自2008年三峡工程175 m试验性蓄水以来，在高程145~175 m之间形成了高差达30 m的水库消落带，导致岸坡岩体节理增生并向失稳方向发展（黄波林等，2019; 刘新荣等，2020; 殷跃平等，2020）。三峡库区岸坡灰岩溶蚀岩体劣化发展诱发地质灾害引起学者广泛关注，典型案例有龚家坊滑坡、箭穿洞危岩体（Wang Luqi et al.，2020）、板壁岩危岩体、棺木岭危岩体、龙门寨危岩体（黄波林等，2021）等。由于灰岩岩体的劣化时效性较长，如何对原位灰岩岩体劣化进行有效的定量刻画和表征是极为重要的工作。前人通过室内试验，研究了灰岩标准岩样在水岩耦合（Wang Luqi et al.，2020）、不同溶液环境（Dochez et al.，2014）、不同应力环境（冯夏庭等，2005）、不同温度（Mutlutürk et al.，2004）下的劣化特征，取得了大量成果。但是室内试验不能揭示原位节理岩体的劣化程度（Hoek et al.，2019）。针对灰岩岸坡岩体劣化，殷跃平等（2020）通过多期次原位跨孔声波测试结合广义H-B准则，深入研究了“完整岩块+岩体结构面”的二元劣化规律，并指出岩体结构面劣化对岸坡的长期稳定性不利。黄波林等（2019）、殷跃平等（2021）、Yan Guoqiang et al.（2022）进一步采用原位消落带典型灰岩岩体测窗影像进行节理统计对比分析，给出了消落带岩体的劣化程度分级标准。

基于岩体节理高清影像进行岩体劣化程度评估是一种直接高效的方法（Huang Bolin et al.，2020）。岩体劣化最直观的表现是节理的新增、延长、变宽（黄波林等，2019; 张永双等，2019; Yan Guoqiang et al.，2022）。传统的基于节理数J_v的统计往往具有较大的局限性，受控于测线位置以及测窗分布，对高度随机分布的节理信息统计极容易产生“偏差”和“纰漏”（卢波等，2005），而且岩体节理长度、条数等统计受观测人员的主观因素影响较高（徐黎明等，2011）。分形理论采用分形维数（fractal dimension）则可以很好地对自然界中不规则线、面、体等天然随机分布物件测度（Mandelbrot，1967），在岩石损伤过程、破碎分形等多个方面取得成功应用（涂新斌等，2005; 蔡建超等，2015; Ju Yang et al.，2018; Kou Miaomiao et al.，2019; 张艳博等，2021）。本文结合前期研究成果，采用分形理论对多个典型消落带岩体测窗节理进行了统计分析，在量化评估岩体劣化程度的基础上，对三峡库区未来表层岩体的劣化发展模式及破碎相变时间进行了预测，对于库区岸坡劣化程度调查及岩体结构劣化的演化趋势预测具有一定指导作用。

1 消落带岩体节理劣化分形计算方法

分形理论中最重要的两个参数为分维数与缺项值。分维数可以反映分形集（岩体节理）所占空间（测窗）大小以及节理的不规则性（平直度）。缺项值可以表征岩体节理集群化程度、分布的均匀性、各向异性。联合分维数与缺项值可以较好地描述消落带岩体节理的生长过程与节理分布的集群化程度。分维数的计算有多种方法（如码尺法、沙盒法等），目前最为常用的方法为计盒法（Gagnepain et al.，1986）。缺项值的计算则多采用滑盒法进行统计分析（Roy et al.，2010）。下面主要针对这两种方法在消落带岩体节理劣化中的应用进行介绍和推导。

1.1 分维数的计盒法运算

计盒法是指采用不同尺度r的盒子去覆盖分形体（如节理曲线、面、孔隙体等），所需盒子数目N（r）将随盒子尺度r变化，且N（r）-r在统计意义上满足：

式中，D_b即为计盒法得到的分形维数。对式（1）两边取对数可得岩体节理lnN-lnr拟合回归线，该回归线的斜率即为分维数D_b。

数学上的严格自相似分形无论选取多大尺度进行计算，其分形维数应是严格一致的。但是在实际中，由于不同统计尺度的精准性以及拍照像素的限制，观测窗只是满足统计意义上的自相似性和无标度特性。在工程应用中，选取的测窗满足工程尺度大小即可，本文中统计的测窗长宽尺度是在5~19 m之间，可以满足代表性岩体的统计分析。由公式（1）、（2）不难求得欧氏几何中线、面、体的分维数D_b分别为1、2、3（图1）。对于消落带岩体节理来讲，初始仅含一条或若干条平直光滑节理时（图1a），测窗岩体节理的分维数D_b=1。随着节理的生长、发育，整个测窗区域极限状态为全面密布节理（图1b），理想状态下此刻分维数D_b=2。根据相关文献研究（Mandelbrot，1982; Gibson et al.，2006），对于同一均匀节理分布区，二维与三维空间对应节理的分维数在统计上近似相差1，即有：

式中，D_b3、D_b2代表三维、二维空间对应的节理分维数。本文采用由美国NIH实验室开发的Image J分形计算程序进行计盒法运算（Smith Jr et al.，1996）。

1.2 缺项值的滑盒法运算

以5×5的节理二元图像缺项值计算为例进行展示说明（图2）。图中仅包含黑色（1）和白色（0）两种数据单元，分别表示研究对象的分形集（如节理）以及其他区域（如固体基质以及非量测对象），黑色像元的不同排列代表了现实岩体节理的不同长度、产状的排列组合。采用一个r×r大小的滑盒从图像左上角开始计数。图2a中滑盒大小r=2，统计每个盒子包含分形集（节理/孔隙）的个数（记为M），盒子每次向右移动一个单元（一般为图像最小像元），当盒子平移到每一行尾部时，需向下移动一个单元，直到滑盒移动轨迹覆盖整个图像/测窗区域；增加滑盒大小r（图2b），重复上述移动，可以获得不同r值对应缺项值。每个滑盒包含分形集个数（记为M）出现的频次记为n（M，r），n（M，r）除以全部滑盒个数得概率分布函数Q（M，r），其中全部滑盒数记为（r_总-r_b+1）²， r_总为图像总长，r_b为滑盒尺寸。

（a）—r=2滑盒法计算示意图，其中红框为滑盒大小，小黑框为节理分形集；（b）—r=4滑盒法计算示意图，其中红框为滑盒大小，小黑框为节理分形集

(a) —schematic diagram of gliding box method (r=2) , in which the red box is the gliding box and the small black box is the joints fractal; (b) —schematic diagram of gliding box method (r=4) , in which the red box is the gliding box and the small black box is the joints fractal

为方便分析概率函数Q（M，r），构建统计动差函数$Z_Q^{（q）}（r）$如下（蔡建超等，2015）：

缺项值Λ（r）可定义为q=2时的统计代差函数除以q=1时的统计代差函数的平方，如下式（蔡建超等，2015）：

在具有平移不变性（理想均匀分布的分形集）的测窗对象分布中，由于Q（M，r）为狄拉克函数，则有$Z_Q^{（2）}（r）={\left[Z_Q^{（1）}（r）\right]}^2$，此时缺项值恒为1，而非平移不变性的分形集缺项要大于1，且分形集分布越不均匀对应缺项值相对越大。如图2所示，当r=2时，$Z_Q^{（2）}（2）=\sum _M{M}^{2}Q（M，2）=2.125$和$Z_Q^{（1）}（2）=\sum _{M}MQ（M，2）=1.375$，代入式（5）得Λ（2）=1.124；当r=4时，$Z_Q^{（2）}（4）=\sum _M{M}^{2}Q（M，4）=30.5$和$Z_Q^{（1）}（2）=\sum _{M}MQ（M，2）=5.5$，可得Λ（4）=1.0083。本文采用Vadakkan（2009）公布的开源程序进行滑盒法计算，且对图2所示算例进行计算、验证，与公式（4）、（5）手算结果一致。

2 消落带岩体劣化分形计算结果

选取典型测窗岩体进行劣化分形计算。这里仅选取岩体劣化强烈、较强烈、一般的3个测窗案例进行计算展示，其他典型测窗采用类似的计算分析步骤。首先处理照片获得二值化节理图像（图3~5），然后采用前文所述程序进行计算分析，获得消落带测窗上下岩体节理的分维数D_b以及缺项值Λ（r）。消落带上下测窗像素精度大小完全相同，为1605×1605左右。根据比拟法，消落带之上A测窗统计数据可视为2008年的B测窗统计数据，消落带处B测窗为2019年拍摄的岩体（殷跃平等，2021）。

（a）—轿顶峰1^#消落带上下测窗岩体节理二值化图像（据殷跃平等，2021改）；（b）—对应测窗岩体节理分维数拟合曲线；（c）—对应测窗岩体节理缺项值曲线

(a) —binary image of rock mass joints in upper and lower measuring windows of Jiaodingfeng 1^# drawdown zone (modified from Yin Yueping et al., 2021) ; (b) —fitting curve of fractal dimension of rock mass joints corresponding to measuring window; (c) —lacunarity curve of rock mass joints corresponding to measuring window

（a）—箭穿洞2^#消落带上下测窗岩体节理二值化图像（据殷跃平等，2021改）；（b）—对应测窗岩体节理分维数拟合曲线；（c）—对应测窗岩体节理缺项值曲线

(a) —binary image of rock mass joints in upper and lower measuring windows of Jianchuandong 2^# drawdown zone (modified from Yin Yueping et al., 2021) ; (b) —fitting curve of fractal dimension of rock mass joints corresponding to measuring window; (c) —lacunarity curve of rock mass joints corresponding to measuring window

（a）—巴雾峡口消落带上下测窗岩体节理二值化图像（据殷跃平等，2021改）；（b）—对应测窗岩体节理分维数拟合曲线；（c）—对应测窗岩体节理缺项值曲线

(a) —binary image of rock mass joints in upper and lower measuring windows of Bawu Gorge drawdown zone (modified from Yin Yueping et al., 2021) ; (b) —fitting curve of fractal dimension of rock mass joints corresponding to measuring window; (c) —lacunarity curve of rock mass joints corresponding to measuring window

图3中轿顶峰消落带岩体可视为岩体劣化强烈的代表；岩体经过11年水位变动后，岩体结构发生明显降级由完整块状向碎裂状发展。由分维拟合曲线可以明显看出，随着岩体劣化，岩体节理逐渐延展、增生，节理分维数发生了明显的升维现象。具体升维数值为：ΔD_b=D_b（2019）-D_b（2008）=0.2146，且对应的缺项值Λ（r）曲线发生明显下降。这表明随着岩体节理的生长、加密，节理分布相对更加均匀，消落带处B测窗节理分布的集群性明显要低于消落带之上的A测窗节理。图4、图5中分维数D_b以及缺项值Λ（r）呈现类似的规律，只是具体的变化量值不同。如图4中岩体劣化较强烈的箭穿洞2号测窗ΔD_b=0.0984，缺项值Λ（r）曲线呈现微弱下降。图5中岩体劣化一般的巴雾峡口测窗ΔD_b=0.0054，缺项值Λ（r）曲线呈现极其微弱下降。

联合图3、图4、图5分维数D_b曲线以及缺项值Λ（r）曲线可以明显看到，分维数拟合曲线相关系数R²都高达0.99，表明岩体测窗节理劣化分形特征明显。其次由于岩体节理劣化的存在，分维数皆呈现“升维”趋势，而缺项值皆呈下降趋势，且劣化程度越大，这种变化趋势愈发明显。表明库岸岩体表层节理的劣化生长、发育，整体是向“均匀密布”测窗的趋势发展的。缺项值Λ（r）曲线显示，随着计盒尺寸r的增加，缺项值Λ（r）对应下降。这其实也体现了随着统计尺寸的增加，岩体节理均匀性增大的普遍客观规律。比如标准岩芯宏观尺度上可视为均质的，但是在微观尺度上则呈现劣化损伤、微观孔隙的高度不均质性，与统计尺度有很大的关系。

Λ（r）除了与分形集本身的离散性、集群性有关外，与计盒尺寸r也密切相关。可见采用可以代表测窗岩体的尺度r来刻画、描述缺项值 Λ（r）是有必要的。对于本文测窗缺项值Λ（r）曲线来看（图3c，图4c，图5c），选取不大于200的像元数进行描述是较为合理的（对应实际尺寸为2~2.5 m）。采用缺项值 Λ（r）曲线陡缓转折点来反推岩体REV（表征单元）尺寸的下限值是一种较为可靠的量化手段。对于本文岩体节理测窗REV尺寸下限值约为2~2.5 m左右，小于测窗整体尺寸5~19 m，可见采用测窗尺寸5~19 m 是可以有效代表岸坡岩体节理的分布情况。

对文献（殷跃平等，2021）中所有典型测窗进行分维数D_b以及升维数ΔD_b值计算、统计如表1所示。在原有定性分级的基础上，发现采用分维数D_b以及升维值ΔD_b可以很好的刻画、划分岩体劣化程度。当ΔD_b≥0.15，岩体劣化强烈；当0.08≤ΔD_b＜0.15，岩体劣化较强烈；当ΔD_b＜0.08，岩体劣化一般（表1）。利用这一分类标准进行岩体劣化程度划分的结果，与采用劣化变量Det、ΔJ_ov、ΔJ_v、ΔGSI值进行划分基本一致。

从表1来看，分维数升维值ΔD_b与ΔJ_ov与ΔGSI等具有很好的相关性，都呈现随岩体劣化程度加剧，数值递增趋势。将表1中的10个案例测窗数据进行相关统计拟合分析，可得分维数升维值ΔD_b与ΔJ_ov、ΔGSI拟合函数如下：

结合公式（6）、（7）与图6、图7可知，分维数升维值ΔD_b与ΔJ_ov以及ΔGSI之间可采用指数函数拟合，相关系数R²分别达0.967、0.923。可见随着岩体劣化分维数的升高，岩体测窗的劣化速度可能会呈现一个加速递增趋势，劣化变量ΔJ_ov以及ΔGSI增加明显。岩体劣化分维数的升维值ΔD_b可以很好地刻画岸坡岩体的劣化程度。

3 讨论

由图1、表1可知，岩体节理的分维数介于1~2之间。理想状态下起始岩体仅含有若干条平直光滑的节理，对应分维数D_b=1；当岩体节理逐渐新生、延展，整个测窗区域极限状态为全面密布曲折节理，此时对应极限分维数D_b=2。这一过程是由完整结构岩块破裂演变为碎裂结构岩体，可将其定义为“劣化破碎相变”，对应时间即为“劣化破碎相变时间”。通过D_b（2019）、D_b（2008）两期分维数可以大体推测岩体节理分维数升维速率公式如下：

式中，λ为节理分维数升维速率；则当岩体节理逐渐发展达到极限分维数D_b=2时，可得所需的劣化破碎相变时间为：

采用公式（8）、（9）对10个测窗计算达到极限分维数所需的劣化破碎相变时间统计汇总见表2。由表2可知，对于不同劣化程度岩体来讲，呈现出岩体劣化愈强烈，岩体劣化至极限破碎相变的时间最短、最快；这和现场调查中的定性认知是一致的。对于岩体劣化为强烈程度的测窗，极限破碎相变时间T＜50年，建议开展专业监测并关注岩体劣化对岸坡稳定性的影响程度。

结合现场调查并借鉴岩石强度损伤演化过程，岩体节理分维劣化演变模式可概化为两种（图8）。A型节理分维数先快速劣化升维之后缓慢升维并不断逼近极限分维数D_b=2。部分较致密岩体如灰岩整体呈现B型类“S”状的劣化升维形式，早期劣化衰变较慢，历经一定的损伤劣化周期后出现一段加速劣化过程，在后期劣化速率再次放缓（Wang Luqi et al.，2020）。随着岩体劣化的发展，节理分维数虽然在逐渐抬升，但是现实中岩体节理的生长、延展是不可能真正达到极限分维数D_b=2。通过公式（8）、（9）结合两期次节理实测分维数进行极限破碎相变时间的计算，实际上是图8中展示的割线估计方法。割线与极限分维数交点计算的极限破碎相变时间是偏小的，计算结果偏保守。由图可知，割线估计法计算的劣化破碎相变时间，已逼近两种演变模式的节理分维曲线拐点，是基本满足工程实际应用的。

在三峡库区中受岩体劣化作用消落带岩体节理逐渐扩展，岩体完整性进一步下降，岸坡稳定性趋于下降。特别是消落带岩体对应的极限劣化破碎相变时刻，岸坡稳定性剧烈下降，可能使得消落带岸坡处于一种几近失稳的临界状态。本文提出的基于岩体节理图像学的分维计算法来预测岩体的极限劣化破碎相变时刻可作为传统定性调查、评估岸坡稳定状态的一种有益补充。极限分维计算有望通过简单、大量的图像调查，如区域高精度无人机、三维激光扫描、航空摄影调查等，对岩体的劣化状态以及可能的破坏时间进行大范围的快速评估。

4 结论

采用分形理论中的分形维数以及缺项值对三峡库区巫峡段灰岩岩体劣化程度、标准划分、可能的劣化破碎相变时间以及发展趋势进行了初步研究，得到如下主要结论：

（1）消落带岩体节理劣化分布具有分形属性。消落带之下的B测窗分维数D_b大于其上部A测窗，呈现升维特性；对应缺项值Λ（r）则小于上部，且随着计盒尺寸增加，节理分布的“集群性”降低，缺项值Λ（r）曲线陡缓转折点可视为REV岩体尺寸下限；整体来看，劣化使得岩体节理生长向更加均匀、密布的趋势发展。

（2）消落带岩体节理分维数的升维值ΔD_b可以定量刻画岩体的劣化程度。当ΔD_b≥0.15，岩体劣化程度为强烈等级；当0.08≤ΔD_b＜0.15，岩体劣化较强烈；当ΔD_b＜0.08，岩体劣化程度为一般。且分维数升维值ΔD_b与 ΔJ_ov、ΔGSI满足指数函数关系。

（3）讨论提出了灰岩岩体可能的劣化破碎相变时间计算公式。采用割线法计算了岩体劣化强烈的消落带极限破碎相变时间为T＜50年。探讨了2种典型节理分维数劣化演变模式，其中A模式岩体劣化分维数持续快速增长，更值得关注。

（4）分维数D_b受观测尺度影响较小，应用范围广，对于岩体劣化调查具有较高的统计稳定性。建议在岩体劣化调查中进行多期次的分形维数计算、调查，为三峡库区岸坡劣化演化预测提供更好的支撑。

参考文献