基于小波变换理论的岩石节理面粗糙度分级表征

2022-08-18元宙昊叶义成罗斌玉阎要锋

煤炭学报 2022年7期

元宙昊,叶义成,2,罗斌玉,阎要锋

(1.武汉科技大学资源与环境工程学院,湖北武汉 430081;2.武汉科技大学湖北省工业安全工程技术研究中心,湖北武汉 430081;3.武汉理工大学土木与建筑工程学院,湖北武汉 430070)

随着“一带一路”国家倡议的提出，在地质条件复杂、地质灾害频发的西部地区开建了一批重大基础工程。这些基础工程的建设在给岩体工程稳定性计算提出新挑战同时还对节理岩体剪切力学特性的准确评估提出了新要求。粗糙度作为影响节理面剪切强度的重要因素之一，节理面粗糙度对节理岩体剪切力学特性的准确评估具有重要意义，为此如何实现节理面粗糙度的准确表征逐渐成为节理岩体力学领域研究的热点问题。

节理面形貌是多尺度的，可将其分为一阶起伏体(初级波纹起伏结构)和二阶起伏体(次级细微粗糙结构)。因一阶、二阶起伏体表面形貌特征存在差异，2者对节理试样剪切力学行为的影响是不同的。陈世江等、齐豫等在未对节理表面形貌进行分解的前提下，提出采用起伏度与起伏角表征一阶起伏体，采用分形维数、轮廓线伸长率表征二阶起伏体，并建立了节理面粗糙度分级表征公式。LIU等通过2种采样间隔(0.5 mm及5.0 mm)将标准JRC曲线分解为一阶、二阶起伏体，通过对一阶、二阶起伏体的统计参数进行计算，建立了粗糙度分级表征公式，并提出分级表征公式可以更准确地描述节理面的粗糙程度。然而，LIU等通过固定采样间隔来实现节理面形貌的分解方法主观性较强。

实现节理面粗糙度分级表征的关键在于准确分离节理面形貌中的一阶、二阶起伏体。一阶起伏体为起伏度较大的低频成分，二阶起伏体为起伏高度较低的高频成分。小波变换理论作为一种多尺度分析方法可以实现信号中的近似分量(低频分量)和细节分量(高频分量)的分类提取。该方法可以为节理面表面形貌中一阶、二阶起伏体的准确分离提供一种新途径。基于小波变换理论，以db8小波基为母小波，以近似分量的标准差开始突变为判据来确定临界分解水平，以Shapiro-Wilk检验(W检验)为检验方法对二阶起伏体的分布进行检验，可实现节理面形貌中一阶、二阶起伏体的分离。该方法多被用于探究粗糙节理面对流体流动的影响，是否适用于节理面粗糙度分级表征仍有待进一步研究。同时，在采用小波变换方法对节理面进行分解时，小波基函数的选取及最优分解层数的确定均直接关系到一阶、二阶起伏体的准确提取，从而对节理面粗糙度分级表征的准确性造成影响。ZOU等通过试错法来选取最优小波基，该方法比较耗时。当信号中存在噪声时较难比较信号与小波基之间的相似性。通过近似分量的标准差开始发生突变为判据来确定临界分解水平时对于不同形貌的表面得到的临界分解水平不同。同时，ZOU等在对二阶起伏体分布特征判别时采用的W检验的适用范围为样本数量小于50的情况，而在采样间隔较小时会出现样本数量较大情况。因此，在小波基优选及临界分解水平判别方面仍需进一步研究。

综上所述，节理面表面形貌是多尺度的，在对未分解的标准JRC曲线的统计参数进行计算时会出现统计参数随粗糙度的增加存在局部凸起现象，该现象会导致部分标准曲线的粗糙程度被低估。该现象的发生主要由于未充分考虑一阶、二阶起伏体的形貌贡献导致。基于此，以小波变换理论为基础，提出能-熵比准则对小波基进行优选，确立适合于标准JRC曲线的最优小波基，结合基于能量保持百分比和信号标准差定义的临界分解水平准则对临界分解水平进行确定，实现节理面中一阶、二阶起伏体的准确分离。然后计算出一阶、二阶起伏体的统计参数，基于这些统计参数建立节理面粗糙度的分级表征公式。最后开展节理面剪切试验，对粗糙度分级表征公式进行验证。研究为节理面粗糙度准确表征提供了新方向。

1 基于小波变换理论的节理面形貌分解方法

1.1 小波变换理论

任何复杂的曲线均可看成不同波长、振幅的正弦波、余弦波的组合体，Barton的10条标准曲线也是如此。小波变换作为一种信号的时间-尺度分析方法，通过伸缩和平移等运算功能可实现对信号进行多尺度的细化分析，其具体表达式为

(1)

式中，wt为小波变换；为伸缩因子；为平移因子；为小波基函数,的复共轭；()为表征信号的函数；为信号的采样间隔。

任何信号均可以通过离散小波变换(DWT)分解为高频(细节)和低频(近似)分量。整个分解过程可以看作一个滤波及卷积的过程，通过将第-1分解水平的近似系数与低通滤波器系数进行卷积，获得第分解水平的近似分量。同理，通过将第-1分解水平的近似系数与高通滤波器系数进行卷积，亦可以获得第分解水平的细节分量。

1.2 小波基优选准则

在对信号进行小波变换时，小波基函数选取的好坏直接决定了信号分解的成败。在小波基的优选准则中，能量准则与熵准则均可用于小波基优选，而能-熵比准则是一个将能量准则与熵准则相结合的小波基优选准则。该准则既考虑了近似分量的能量还考虑了能量的谱分布。基于巴什瓦尔定理，小波变换得到的离散小波分解系数反映了信号的能量特征。能量准则是以近似分量的能量为度量准则提出的一种基于最大能量匹配准则的小波基优选准则。在对信号进行小波变换时，近似分量包含的能量越高，选择的小波基越优。近似分量中包含的能量可表示为

(2)

式中，()为第个离散小波分解系数(第分解水平)；为第分解水平的离散小波分解系数的总数。

为保证信号特征的有效提取，还应考虑能量的谱分布。在信息论中，香农熵是与随机变量相关的不确定性的度量。小波系数的能量分布可以由香农熵描述，即

(3)

=|()|()

(4)

式中，为小波系数的能量概率分布。

式(3)，(4)中，当=0时，log=0，且0≤()≤log。若小波系数只有一个不为0，则熵为0；若所有系数的能量分配概率相同，则熵值为log。节理面的一阶起伏体具有显著的低频特征，当节理面形貌被分解为多层时，最佳小波基应在低频处产生大的幅度系数，在其他频率处产生可忽略的幅度系数。从信号中提取的能量越高，信号的小波变换越有效。同时，熵越低，能量密度越高。因此，将能量准则与熵准则结合起来，提出一个新的小波基优选准则—能-熵比准则，表达式为

()=()()

(5)

最优小波基从被分析的信号中提取最大的能量，经变换得到的离散小波分解系数的香农熵最小。能-熵比越大，表明选取的小波基越优。

1.3 临界分解水平判别准则

在信号处理领域，临界分解水平的判别问题就是最优分解层数的确定。从信号处理角度来讲，当分解层数过多，会造成有用信息的丢失，即节理面形貌中的一阶起伏体被当作二阶起伏体提取；分解层数过少则降噪的目的很难达到，即节理面形貌中的二阶起伏体不能被完全提取，从而不能实现节理面表面形貌中一阶、二阶起伏体的准确分离。节理面形貌中的一阶起伏体反映了宏观起伏特征，而二阶起伏体可视为服从均值接近0的高斯分布的白噪声。能量保持百分比′可用于表征分解后的近似分量保持原信号特征的能力，信号标准差为分解后细节分量偏离原信号数据平均值的度量。通过分析′及随分解层数的变化规律可实现临界分解水平的判别。当达到临界分解水平前，一阶起伏体仍包含在近似分量中，近似分量仍可反映宏观起伏特征，因此′与不会出现显著变化，当超过临界分解水平后，部分一阶起伏体被作为二阶起伏体提取，近似分量将难以准确表征宏观起伏特征，从而导致′与开始发生显著变化。因此，将′与开始发生显著变化时作为临界分解水平判别准则。其中，能量保持百分比′定义为

′=

(6)

式中，为近似分量的能量；为初始信号的能量。

信号标准差定义为

(7)

式中，′()为近似分量的信号序列；()为原信号序列；为信号序列包含信号的个数。

2 基于小波变换的标准JRC曲线的分解结果

2.1 标准JRC曲线轮廓坐标获取与标准JRC曲线重构

Barton提出的10条标准JRC曲线是以图片形式给出，计算其统计参数的前提是进行数字化处理，即获取各条曲线的坐标值。目前，灰度图像处理方法在标准JRC曲线坐标值获取方面已被广泛应用。与其他方法相比，该方法的人为误差较少，其主要步骤为：① 标准JRC曲线图片的获取；② 曲线分割；③图片灰度处理；④ 二值化处理；⑤ 曲线中心线提取；⑥ 坐标计算。具体流程如图1所示。

注:xk,yk为沿水平方向第k个像素点的横、纵坐标；μ为每一个像素点的尺寸；n0为标准JRC曲线沿水平方向包含的像素点的个数； L为标准JRC曲线沿水平方向的长度；l为灰阶矩阵中第k列中灰阶较小的像素点最大行号与最小行号的平均值。图1 灰度图像处理流程Fig.1 Gray scale image processing flow

根据灰度图像处理结果，对10条标准JRC曲线进行重构，重构结果如图2所示。

图2 重构后的标准JRC曲线Fig.2 Standard roughness profile after reconstruct

为验证重构的标准JRC曲线准确性，对曲线的波幅()、粗糙度统计参数和-1进行了统计，并与文献[14,17,39-40]中给出的结果进行对比，结果如图3所示。

由图3可知，参数，，-1均表现出相同的变化趋势。表明采用灰度图像处理方法来实现标准JRC曲线的数字化是合理有效的。

2.2 确定最优小波基

Matlab小波工具箱中预定义的离散小波系包括：Daubechies(db)小波系、Coiflet(coif)小波系、SymletsA(sym)小波系、双正交小波Biorthogonal(biorNr.Nd)小波系。采用较为常用的db(1～10)，coif(1～5)，sym(2～10)三种小波系来选取最优小波基。基于1.2节提出的小波基优选准则，采用式(2)，(3)计算了各小波基的能量和熵，结果分别如图4，5所示。

在获得各小波基的能量和熵后，采用式(5)计算各小波基的能-熵比，计算结果见表1。

根据能-熵比优选准则，将表1中的各小波基的能-熵比从大到小进行排序，并根据其序次给出相应的分值(排序第1得10分，第2得9分，依此类推，第10得1分)，排序结果见表2。从表2中可以看出，coif5小波基的得分最高，因此将coif 5作为标准JRC曲线分解的最优小波基。

图3 统计参数Fig.3 Statistical parameters

图4 各小波基的能量Fig.4 Energy values of each wavelet basis

图5 各小波基的熵Fig.5 Entroppy values of each wavelet basis

表2 基于能-熵比准则(R(j))的小波基排序Table 2 Sorting wavelet basis based on energy-entropy ratio criterion(R(j))

2.3 标准JRC曲线的分解结果

2.3.1 标准JRC曲线的小波变换结果

基于2.2节确定的最优小波基，应用Matlab的小波工具箱对10条标准JRC曲线进行小波变换处理。图6为1号标准JRC曲线经过8个分解水平后的小波变换结果。其中a为近似分量，d为细节分量(=1,2，…,8)。其他标准JRC曲线的小波变换结果均与1号标准JRC曲线类似，此处只对1号标准曲线进行展示。

由图6可知，当分解水平较低时，近似分量a中仍包含曲线的宏观起伏/波动，细节分量d为频率较大、幅度较小的波动。a和d的形貌特征随分解水平的增加而不断变化。当分解水平超过一定值后，d的频率开始变小，幅度开始变大。

2.3.2 临界分解水平的确定

根据式(7)，(8)对各条标准曲线在各分解水平下的近似分量的′和进行了计算，计算结果如图7所示。

图6 1号标准JRC曲线小波变换结果(8水平分解后)Fig.6 Wavelet transform result of No.1 standard JRC profile(After 8 level decomposition)

图7 各分解水平下标准JRC曲线近似分量的能量保持比(E′)及信号标准差(S)Fig.7 E′ and S of the approximate components of the standard JRC profile at each decomposition level

由图7可知，当分解水平超过4时，各条标准曲线的近似分量的能量保持比开始快速下降，而信号标准差开始快速上升，根据1.3节提出的临界分解水平判别准则，确定第4分解水平为临界分解水平。

为进一步验证1.3节提出的临界分解水平判别准则的有效性，对各条标准曲线的二阶起伏体的分布特征进行了统计，结果如图8所示。各条标准曲线的均值基本处于0附近，相关系数基本大于0.95，表明各条曲线的二阶起伏体均满足高斯分布，这与以往研究的结果是一致的。综合图7，8的结果，确定10条标准JRC曲线的临界分解水平均为第4分解水平。

2.3.3 标准JRC曲线分解结果

在确定最优小波基与临界分解水平后，将10条标准JRC曲线分解为一阶起伏体和二阶起伏体，如图9所示。由图9可知，各条标准曲线的一阶起伏体可以反映原始曲线的宏观起伏特征，而二阶起伏体可以反映原始曲线中的幅值较低频率较高的小尺度的波动。一阶、二阶起伏体的形貌特征满足ISRM给出的定性描述，表明提出的节理面表面形貌分解方法是可行的。

图8 二阶起伏体高度分布统计Fig.8 Frequency count of unevenness

图9 标准JRC曲线分解结果Fig.9 Standard JRC profile decomposition result

3 节理面粗糙度分级表征

3.1 统计参数变化规律

(8)

(9)

式中，为标准JRC曲线离散点的个数；Δ为标准JRC曲线沿水平方向的取样间隔。

采用式(8)，(9)对原始标准JRC曲线及其各分解水平下近似分量和细节分量的统计参数(,-1)进了计算，计算结果如图10所示。

在此基础上，对各标准JRC曲线中的一阶、二阶起伏体的统计参数进行了计算，结果如图11所示。

图10 标准JRC曲线各分解水平下的统计参数Fig.10 Statistical parameters at each decomposition level of the standard JRC profile

图11 标准JRC曲线及其一阶、二阶起伏体的统计参数Fig.11 Statistical parameter for standard roughness profile and the waviness and unevenness

3.2 节理面粗糙度分级表征公式

为实现节理面粗糙度的定量表征，一系列的未分级的表征公式被提出。JANG等在对已有公式分析的基础上，提出已有公式可以表示为

JRC=()+

(10)

式中，为统计参数；，，为拟合参数。

为准确表征节理面粗糙程度，采用一阶、二阶起伏体的统计参数对节理面粗糙度进行分级表征。在对陈世江、LIU、齐豫等提出的表征公式进行总结的基础上，提出新的节理面粗糙程度的分级表征公式为

JRC=()+()+

(11)

式中，，为拟合参数。

利用Matlab按式(10)，(11)形式对节理面粗糙度系数进行拟合，拟合结果如图12，13所示，拟合公式分别为式(12)～(15)。

JRC=8377()0333 5-4418=0928 9

(12)

JRC=9149(-1)0207 2-355=0931 0

(13)

(14)

(15)

图12 基于标准JRC曲线的统计参数的拟合结果Fig.12 Fitting results by the statistical parameter of standard roughness profile

由图12，13可知，用未分解的标准曲线的统计参数，-1表征节理面粗糙度系数时的分别为0.928 9 和0.931 0。采用分解后的一阶、二阶起伏体的统计参数来表征节理面粗糙度系数时的分别为0.956 1和0.955 6。对比式(12)与(14)、式(13)与(15)可知，式(12)与(13)对节理面上所有的起伏体都采用同一权重来表示，式(14)与(15)中分别考虑了一、二阶起伏体的权重。同时因一阶、二阶起伏体对节理面剪切力学行为的贡献存在差异，与采用同一权重相比，分别考虑一阶、二阶起伏体的权重会更为合理。通过分析式(14)，(15)可以发现，一阶起伏体的统计参数的指数部分均小于1，而二阶起伏体的统计参数的指数部分均大于1，同时一、二阶起伏体的统计参数均小于1。从数学角度可知，当底数小于1时，指数越大，数值越小，指数越小，数值越大。为此，根据式(14)，(15)对一阶、二阶起伏体的统计参数代表的粗糙度系数进行了计算，如图14所示。

图13 基于一阶、二阶起伏体统计参数的拟合结果Fig.13 Fitting result by the statistical parameters of waviness and unevenness

图14 一阶、二阶起伏体的粗糙度系数Fig.14 Roughness coefficient for waviness and unevenness

由图14可知，基于一阶起伏体的统计参数计算出的粗糙度远大于二阶起伏体的统计参数的计算值。表明通过一阶、二阶起伏体的统计参数的权重系数控制可以对其形貌贡献进行调整。同时，由式(14)，(15)对标准曲线的粗糙度进行了计算，计算结果如图15所示。由图15可知，JRC随着参考值的增加也呈现出递增趋势，未出现局部凸起现象。表明通过分级表征公式可以避免采用未分解的标准JRC曲线的统计参数表征时出现的粗糙度系数被低估的现象发生。

图15 JRC与JRCc对比Fig.15 Comparison between JRC and JRCc

4 节理面粗糙度分级表征公式有效性验证

为验证节理面粗糙度分级表征公式的有效性，针对节理试样开展了不同法向应力条件下的压剪试验。根据得到的剪切试验结果，通过对Barton剪切强度公式逆向展开求出节理面粗糙度系数。通过对比JRC与分级表征公式的计算值对分级表征公式的有效性进行验证。

4.1 节理试样的制备

制备具有相同表面形貌的节理面试样是在不同条件下展开剪切试验的基础。目前，“三维扫描+3D打印+相似材料浇筑”技术被广泛应用于具有相同形貌的节理面试样制备。该方法主要包括初始节理面试样制取，节理面表面形貌获取，节理面模板制备以及节理面试样浇筑这4个步骤，具体制备流程如图16所示。

4.2 剪切试验

试验采用的节理面试样尺寸为100 mm(长)×100 mm(宽)×100 mm(高)，采用浇筑相似材料的方式来制备节理面试样。试验采用的石膏为高强石膏，与传统建筑石膏相比，该石膏材料易于搅拌及脱模，凝固速度快，力学性质稳定，均质且各向同性好。制备完成的试样如图16(d)所示。剪切试验在微机控制电子式岩石直剪仪(图17)上展开试验,采用的边界条件为恒定法向应力(CNL)。由于锚固工程支护多处于地下巷道等围岩2～3 m以浅的范围内。根据相关经验，选取法向应力分别为1，2，4 MPa。加载速率根据ISRM的建议值设定，首先以0.3 mm/min的加载速率在上剪切盒顶部施加法向荷载至目标值，当加载至目标值后保持恒定。然后以0.3 mm/min的加载速率在下剪切盒上施加剪切力，当剪切位移达到8 mm后停止加载。

上述4种表面形貌及3个法向应力下的剪切试验结果如图18所示。随着法向应力的增加，节理试样的剪切强度逐渐增大。

此外，采用单轴压缩试验对水膏质量比为0.3∶1.0的石膏试样(50 mm×100 mm)的强度进行了测试，采用测斜试验对表面平滑的节理试样的基本摩擦角进行了测试，获得试样的基本物理力学参数见表3。

图16 节理面试样制备流程Fig.16 Preparation process for joint specimen

图17 微机控制电子式岩石直剪仪Fig.17 Computer controlled electronic rock direct shear apparatus

4.3 节理面粗糙度系数反算

基于4种形貌的节理试样在1，2，4 MPa作用下的剪切强度，结合Barton剪切强度公式(式(16))，可反算出节理面的各节理面的粗糙度系数(式(17))。

(16)

(17)

式中，为节理面剪切强度；为节理面法向应力；JCS为节理壁面强度，一般按试样的单轴抗压强度取值；为基本摩擦角。

采用式(17)与石膏试样的基本物理力学参数(表3)对4个节理面的表面粗糙度系数进行反算，结果见表4。

表3 石膏试样基本物理力学参数Table 3 Basic physical and mechanical parameters of gypsum specimen

表4 各节理面试样表面粗糙度系数反算值Table 4 Inverse calculated value of surface roughness coefficient of each joint specimen

图18 剪切应力-剪切位移曲线Fig.18 Shear stress vs shear displacement curves

4.4 节理面分级表征公式验证

为实现节理面分级表征公式的验证，首先通过“小波变换+临界分解水平判别准则”将节理面表面形貌进行分解，分解过程如图19所示。此处只展示1号节理面表面形貌的小波变换过程。

图19 1号节理面小波变换过程Fig.19 Wavelet transform process of No.1 joint specimen

以4 mm为间距，沿节理试样剪切方向取26个剖面计算节理面粗糙度系数，如图20所示。

基于未分级表征公式(12)，(13)对4个节理面的104条剖面的粗糙度系数展开计算，计算结果如图21所示。从图21可以看出，除3号节理面的粗糙度系数分布不太均匀外，其余剖面的粗糙度系数均较为均匀。

基于分级表征公式(14)，(15)对4个节理面的104条剖面的粗糙度系数展开计算，计算结果如图22所示。从图22可以看出，除3号节理面外，其他节理面的粗糙度系数分布均较为均匀。

图20 1号节理面剖面布置Fig.20 Section layout for No.1 joint specimen

图21 节理面粗糙度计算结果(Z2,RP-1)Fig.21 Joint surface roughness calculation result by parameter(Z2,RP-1)

图22 节理面粗糙度计算结果Fig.22 Joint surface roughness calculation result by the

对比图21，22可知，基于节理面粗糙度分级表征公式的计算结果大于未分级表征公式的计算结果。为比较试验反算值与未分级表征公式(12)，(13)、分级表征公式(14)，(15)的计算结果之间的差异，对26个剖面的粗糙度系数求平均值，结果如图23所示。

图23 节理面粗糙度系数的平均值与试验反算值对比Fig.23 Comparison between the average value of the roughness coefficient and the inverse calculated value

由图23可知，1～4号节理面粗糙度系数的平均值均小于试验反算值，其中未分级表征公式的节理面粗糙度计算值与试验反算值之间的差异大于基于分级表征公式的节理面粗糙度计算值。LIU等在比较剪切试验反算出的节理面粗糙度系数与基于统计参数的计算值时也发现该问题，并提出节理面剪切过程中只有部分较为粗糙的剖面来提供抵抗力，因此应该对较为粗糙的剖面的粗糙度系数进行比较，而非通过平均值进行比较。因此提出从节理面中选取部分较为粗糙的剖面来表征节理面的粗糙系数，统计结果如图24所示。