基于小波變換理論的巖石節理面粗糙度分級表征

元宙昊,葉義成,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所示。

圖24 節理試樣粗糙度系數(部分剖面統計參數平均值) 與試驗反算值對比Fig.24 Comparison between the average value of the roughness coefficient and the inverse calculated value

由圖24可知，與未分級的表征公式計算出的節理面粗糙度系數相比，基于分級表征公式計算出的節理面粗糙度系數與剪切試驗反算出的節理面粗糙度系數之間的差異較小。表明基于節理面粗糙度系數分級表征公式是有效的，且比基于未分級的表征公式更加準確。

5 結 論

(1)提出了能-熵比準則用于小波基優選，采用該準則確立了適用于節理面表面形貌分解的最優小波基為coif5；結合能量保持比和信號標準差建立了具有明確數學定義的臨界分解水平判別準則，確立了第4分解水平為臨界分解水平。

(2)基于小波變換理論，結合小波基優選準則及臨界分解水平判別準則，提出了一種新的節理面表面形貌分解方法。實現了標準JRC曲線中一階、二階起伏體的準確提取。

(3)探明了取樣間隔較低時未分解的標準曲線的統計參數隨粗糙度的增加出現局部凸起現象的原因，當采用未分解的標準曲線的統計參數表征節理面粗糙度時未考慮一階、二階起伏體的統計參數的貢獻差異，從而導致局部凸起現象出現。

(5)采用建立的分級表征公式計算了節理試樣的粗糙度系數，并與剪切試驗反算的節理面粗糙度系數進行了比較，結果表明分級表征公式的計算值與剪切試驗的反算值吻合較好，驗證了粗糙度分級表征公式的可靠性。