基于分形損失理論的粗糙節理巖體中應力波波速研究

李業學，謝向東，秦 麗

(1.襄樊學院 建筑工程學院，湖北 襄樊 410053;2.四川大學 水利水電學院，成都 610065)

爆破仍然是當今世界在地下開采過程中采用的最為廣泛的方法之一。爆炸所釋放的能量部分用于破巖，部分則以應力波形式自爆源向外傳播，遭遇裂縫、節理等不連續結構面后，將產生復雜的應力狀態，對工程體與地質體的安全可能造成隱患，因而，為了達到防患于未然的目的，深入探討節理對應力波傳播的影響規律就顯得很有必要。

許多研究者［1－5］探討了巖石特性和聲波波速的相互關系，發現這兩者密切相關，存在大量影響巖石中應力波波速的因素，這些因素大致可以包括:巖石種類、密度、顆粒尺度和形狀、空隙率、各項異性、空隙水、裂隙斷層不連續結構面、圍壓和溫度。其中，巖石中包含的節理是影響應力波波速的主要因素之一。

學者們通過分析巖石中微裂隙、結構等特點，探討了描述節理的各參量對應力波波速的影響規律，比較有代表性的有:Oda［6－7］從巖體節理、裂隙的角度提出了用裂隙張量來表示巖體中裂隙的幾何形態(裂隙的密度、尺寸和方向)的方法，并建立了裂隙張量與波速之間的關系。文獻［8］提出了考慮裂紋寬度影響的波速計算方法。文獻［9］將節理劃分為微觀節理與宏觀節理，分別探討了它們對應力波波速的影響規律，類似的研究還有［10－12］。研究者們已經就節理對波速的影響這一問題做了大量有意義的工作，并取得了豐碩的研究成果，但縱觀以上研究，他們的研究工作均基于節理面光滑這一假定，事實上，經過億萬年成巖過程形成的巖石，所包含的節理可以說必然粗糙不平，為了研究成果更貼近工程實際，力學家們繼續探討了節理面幾何構型對應力波傳播的影響規律。比較有代表性的有，文獻［13］通過在巖石表面縱橫刻畫刻痕模擬節理面的粗糙性，探討了應力波波速與節理面粗糙性的關系;但是，該文獻描述粗糙節理構型參數(節理粗糙度系數FRC)是作者規定的，存在著大量經驗方法，這些在理論是上不完善的至少是不可接受的，隨著分形幾何的誕生，其幾何構型通常采用分形維數來描述。

分形維數計算理論自分形幾何由Mandelbrot［14］創立以來，其研究在短短數十年間就取得了累累碩果，從三角形棱柱法［15］，投影覆蓋法［16］，立方體覆蓋法［17］直至改進立方體覆蓋法［18］，計算方法在理論上可以說已經發展到了近乎完美的地步，但是維數計算過程中存在瓶頸—激光表面儀的使用。該儀器的掃描視野有限，相對巨大尺度的工程體經常顯得無能為力，因而，學者們隨后探討了粗糙表面維數的間接計算法——圖像維數，如:二值化圖像維數［19］，灰度圖像維數［20］，在一定程度上緩解了這一困難，但問題的關鍵在于現實生活中的圖像是真彩圖，二值化圖像與灰度圖像僅僅是真彩圖像的理想狀態與簡化結果。因而，為了使研究成果更切合工程實際，探討真彩圖像的維數計算方法勢在必行。

本文擬定探討真彩圖像維數的計算方法，提出相應計算理論，并基于損傷理論推導應力波在節理巖石中傳播波速的解析解，建立圖像維數與應力波波速間的關系。

1 應力波穿越分形節理時的波速

在三維一般情況下，當忽略體力影響時，運動方程可表示為:

其中:Vi(i=X、Y、Z)為質點速度，σmn(m，n=X、Y、Z)為應力;ρ是體積密度。對于均勻，各向同性，線彈性介質，可利用Lame形式的Hooke定律代入后消除一些參數，然而，眾所周知，通過億萬年成巖過程形成的巖石內部包含大量微孔隙，裂紋，夾雜甚至宏觀裂縫等一系列不連續結構面，它不可能是一種理想的均質各向同性體，顯然并不滿足Hooke定律的應用條件，一條可行的措施是將包含節理的巖石材料從幾何構型與物理性質上進行等效，使之滿足Hooke定律應用條件，必須注意的是:在理論推導過程中，所有參數必須代之以對應的等效物理參數?；谏鲜龅刃Ъ俣傻?

其中:μ、λ 為 lame常數;εij(i，j=X，Y，Z)為應變;ν為泊松比。uX、uY、uZ為對應位移分量。

對式(2)中三個式子分別對X，Y，Z微分再相加則可得:

這是對體積膨脹Δ的線性雙曲型偏微分方程，表示體積膨脹Δ以波速Vp傳播，即:

上式也可寫為:

其中:Vp為應力波波速;Ed為有效彈模。

巖體內包含大量不連續結構面，或者說巖石材料存在一定程度的損傷，且材料損傷演化具有分形性質［16，17］，因而本文結合損傷力學原理與分形理論，給出了一種兼顧損傷細觀特征描述和宏觀損傷力學分析需要的分形損傷變量ω(d，δ)，可以表示為:

其中ω0是歐氏空間表觀損傷變量;d為節理面分形維數;de為典型域的Euclidean維數;δ是量測尺度。由于實際粗糙巖石表面在Euclidean空間中非二維的面，而是三維的體，故其歐式維數應該是3，即de=3，所以有:

引入損傷變量后，有效動態彈性模量可表示為:

其中Ed0為初始彈模。

將式(7)代入式(8)得:

將式(8)代入式(5)得:

依據上式，我們考慮兩種極限情況，其一，如果巖石材料處于理想的完整狀態，或者說沒有出現損傷狀態，此時有:

另外，依據分形理論，當材料內沒有孔隙、裂紋等不連續結構面存在時，此時的巖石材料在分形空間中的維數為3，即:

將式 (11)、式(12)代入式(10)可得:

分析式(13)可知，此式正是通過經典彈性波理論導出的應力波波速計算公式。

考慮另外一個極限情況，巖石完全開裂，即發生了完全損傷，此時巖石損傷程度為:

將式(14)代入式(10)可得:

式(15)表明，應力波不能穿越完全開裂的巖石，這與實際情況是完全相符的。

通過極限狀態分析可知，本文基于分形損傷理論導出的應力波波速的解析解是經典理論公式的推廣，在實際工程中將更具有普適性，而在經典彈性理論框架內推導出的應力波速計算公式是分形空間中波速公式在理想狀態下的特例。

2 參數研究

為了揭示波速隨著表面維數變化而改變的規律性，從式(10)不難看出，必須分別探討粗糙表面分形維數與表觀損傷變量的計算方法。以下各自給出上述兩參量(維數、表觀損傷變量)的計算過程，并由此導出考慮節理面分形效應的應力波波速。

2.1 節理面圖像分形維數計算

工程尺度上的節理大小通常大大超出激光表面儀的掃描視野，使得計算節理面維數存在很大困難，這往往可能成為探討粗糙節理對應力波波速影響規律的瓶頸。為了克服上述困難，本文在描述節理粗糙度時采用節理面圖像維數替換“普通”節理面維數，以下探討真彩圖像維數計算的相關理論。

2.1.1 顏色表面的模型構建

如圖1所示，RGB顏色模型是一個加色模型，多種基色的強度加在一起生成一種顏色。顏色單元體邊界中的每一個顏色點可表示一個三元組(R，G，B)，三分量R、G、B取值范圍均在區間［0，1］內。依據這一模型，任何一種顏色 Cλ在 RGB坐標中可表示為如下矢量:

圖1 顏色單元體Fig.1 Color units

一副圖象由多個象素點構成，每一象素點對應一種RGB模式的顏色，由于任何一種RGB顏色可表示為形如式(16)的矢量，所以，在每個象素點處均可構建一個三維矢量，如:矢量 Aa、Bb、Cc、Dd(見圖 2)。連接各矢量的終端組成2～3維的粗糙曲面(如圖2)，即為所構建的顏色表面。

2.1.2 真彩圖象分維計算方法

由上述的顏色表面構建程序可以看到，每個象素點與該象素點R、G、B值所構建的矢量都是以各自象素點所在位置為坐標原點，各矢量所在的坐標空間不統一，因而有必要將所有矢量統一至同一坐標系下。具體算法如下:

(1)對圖象所在平面進行網格剖分，網格間距選為象素間距δ;

(2)建立如圖1所示坐標系統，坐標原點設在A點，其中，詳圖1中A、B、C、D四點在該坐標系下的坐標分別為(iδ，jδ，0)、((i+1)δ，jδ，0)、(iδ，(j+1)δ，0)、((i+1)δ，(j+1)δ，0)，所以，顏色表面上對應的 a、b、c、d 四點坐標分別是:(Ra+iδ，Ga+jδ，Ba)、(Rb+(i+1)δ，Gb+jδ，Bb)、(Rc+iδ，Gc+(j+1)δ，Bc)、(Rd+(i+1)δ，Gd+(j+1)δ，Bd)。由此可得:顏色表面上任一點的坐標在統一坐標系下可表示為:(R+i·δ，G+j·δ，B)(i=0，1，2…m，j=0，1，2…n)(m 與 n 分別為網格點的行列數)。依據上述方法，我們獲取了粗糙顏色表面的三維坐標。

圖2 顏色表面示意圖Fig.2 Schematic figure on color surface

2.1.3 圖像維數計算

由顏色表面構建過程不難知道，顏色表面上的點分布不均勻，在分析該顏色表面的分形特性前，有必要利用泛克里金法對三維坐標數據網格化。基于周宏偉［21］提出的分形維數計算理論－立方體覆蓋法，計算其分形維數，該分維即為真彩圖象分形維數。具體過程如下:如圖3所示，在平面XOY上存在一正方形網格，網格尺寸為δ，正方形的四個角點處分別對應四個高度 h(i，j)，h(i，j+1)、h(i+1，j)和 h(i+1，j+1)(1≤i，j≤n－1，n為每個邊的量測點數)。用邊長為δ的立方體對粗糙顏色表面進行覆蓋，計算覆蓋區域δ×δ內的立方體個數，即在第i，j網格內，覆蓋粗糙面的立方體個數 Gi，j為:

式中INT為取整函數。

則覆蓋整個粗糙顏色表面所需的立方體總數為:

改變觀測尺度再次覆蓋粗糙顏色表面，計算覆蓋整個表面所需的立方體總數，若粗糙顏色表面具有分形性質，按分形理論，立方體總數G(δ)與尺度δ之間存在如下關系:

式中D為粗糙顏色表面自相似分形維數。

圖3 立方體覆蓋法Fig.3 Cubic covering method

基于上述原理，本文計算了所有節理面圖像維數(見表2)，限于篇幅，本文僅給出了一組典型的巖石節理表面圖像及其對應的粗糙顏色表面形貌圖、維數計算的log－log圖(見圖4)。其中δ指量測尺度;G是覆蓋盒子總數。

2.2 表觀損傷變量

依據損傷力學的相關理論，本文的表觀損傷變量采用缺陷體積與名義體積之比這一定義，即:

其中:V為缺陷體積，V0為名義體積，V為損傷后的體積。

對三維坐標數據中的高度坐標進行縱橫方向搜索，找出高度值中的最大者與最小者，即:

圖4 7#巖石節理面圖像、對應顏色表面圖、維數計算log-log圖Fig.4 Image of 7#rock joint surface，picture of corresponding color surface and bi－ logarithm plot used to computer dimension

其中p，q分別為縱橫方向上的掃描點數。

名義體積可表示為:

由掃描坐標可求得巖石損傷后的體積:

將式(22)、式(23)代入式(20)可得:

2.3 應力波波速的計算

在除預制裂紋外巖樣內其它部分致密且均質這一假定近似滿足的前提下，一組巖樣損傷程度不斷增加，可等效認為是某一個巖樣損傷的動態持續增大，因此存在一個沒有損傷或損傷度較小的初始狀態，取此狀態的彈性模量作為彈模的計算初始值。為了計算便利，并考慮計算精度要求，本文采用沒有預制裂紋的致密花崗巖的彈模，作為巖樣損傷演化過程中的初始值Ed0。其它計算參數(如:體積密度)通過常規試驗測定，具體見表1，將表中參數代入公式(10)，即可求得考慮節理面圖像分形效應后的應力波波速(見表2)。

表1 波速計算參數表Tab.1 Table on parameter for computering wave speed

表2 圖像分形維數與應力波波速Tab.2 Image dimension and velocity of stress wave

3 節理巖石的超聲波試驗

為驗證節理巖石中波速理論推導正確與否，本節通過節理巖石的超聲波試驗，測定超聲波穿越分形節理的波速，建立波速與圖像維數間的關系，并與理論結果進行對比分析。

本試驗采用四川大學MTS實驗室超聲波儀，該設備由超聲波信號激發探頭，信號接收探頭以及顯示波形的示波器組成，當超聲波信號由激發探頭產生后，經巖石介質傳至預制節理(見圖5)，由于節理與巖石的波阻抗差異，此時超聲波將發生復雜的透反射現象，部分反射回發射端，部分超聲波穿越節理傳至巖樣的另一端，被接收探頭感知，由示波器采集并保持下來，限于篇幅，本文僅給出了巖樣的采集結果(見圖6)。分析保存的波形，測量巖樣的長度，可分別求取超聲波傳播的時間與距離，由此計算超聲波在節理巖樣中的傳播速度，計算成果見表2。

圖5 超聲波儀示意圖Fig.5 Schematic figure on ultrasonic apparatus

圖6 超聲波波形Fig.6 Ultrasonic form

4 分析與討論

依據式(10)計算出的波速(見表2)，我們繪制了波速隨維數變化關系曲線圖(見圖7)，它揭示了應力波穿越不同分形節理時的波速變化規律，并清楚顯示:

(1)隨著巖石節理面圖像分形維數增加，或者說節理面的粗糙度增大，應力波穿越該節理時的波速相應減小。反之，如果節理面圖像維數減小，則應力波波速隨之增大。兩者之間呈現顯著的非線性變化關系。這一規律能佐證兩點:① 節理面與節理面真彩圖像具有一致的分形特性，即:節理面越粗糙，分形維數值越大，則節理面圖像從“顏色”角度體現的粗糙度越高，對應的圖像分形維數也越大，這與文獻［20］的“節理面與其灰度圖具有一致的分形性”結論是一致的，所不同的是，本文研究的是真彩圖像分形特性，而文獻［20］僅探討了灰度圖的分形特性。② 對于兩參數間所表現的非線性關系，節理面的不規則性應該“功不可沒”，至少可以認為上述規律的出現部分是來自它的“貢獻”。

(2)盡管應力波波速隨著維數增大持續減小，但波速隨維數減小的速率在不同的維數區間里而有所差異。當圖像分維在區間［2.01203，2.01690］時，在節理面圖像維數增大過程中，應力波波速減小速率相對較小;而當維數大于2.01690且小于2.04116時，應力波波速減小速率加快，約為上一區間的波速變化速率的6.7倍。隨著圖像維數的進一步增大，當處于區間［2.04116，2.04626］時，波速減小速率反而減緩。

(3)通過與文獻［13］的對比分析知，文獻［13］與本文得出了類似的研究結論，即:隨著節理面粗糙度的增加，應力波波速相應減小。但文獻［13］描述節理面粗糙度是文中作者自定義的經驗參數——斷面粗糙系數(fracture roughness coefficient)，且該文中的結論為探索性的試驗成果，而本文通過理論推導得出了應力波在節理巖體中波速的解析解，將描述粗糙性的參數由圖像維數取代，克服了文獻［13］中的經驗主義問題，并將試驗研究結論上升到理論高度，是節理巖體中應力波波速研究的一點改進。

圖7 波速與圖像維數間的關系曲線圖Fig.7 Relation curve between wave velocity and image dimension

5 結論

綜合運用分形幾何學、圖形圖像學、損傷力學的相關理論，導出節理面圖像維數與應力波波速的解析解，并通過數值方法分析了圖像維數改變時穿越節理面應力波波速的變化規律，主要結論包括:

(1)基于真彩圖像的像素點顏色可以表示為一個三維空間中矢量的基本假定，構建了節理面圖像的“粗糙顏色表面”，依據維數計算原理——立方體覆蓋法，計算出該“顏色表面”的分形維數，此維數即為節理面圖像的分形維數。并通過實例對該方法進行了驗證，結果表明該圖像維數計算方法合理可行。

(2)采用分形損傷的基本原理，從構型與參數等效角度推導了應力波在分形節理巖石中傳播波速的解析解，從理論上建立了應力波波速與節理面圖像維數間的定量關系，即:

(3)通過數值方法，探討了應力波波速與節理面圖像維數之間的相互關系。結果表明:隨著節理面圖像維數增大，或者說節理面的粗糙度增大，則在節理巖石中傳播波速隨之減小，兩者間呈現非線性變化關系，且在不同的維數區間，波速隨維數的變化速率存在差異，上述規律至少能說明兩點:① 節理面與其真彩圖像具有一致的分形性;② 節理面的不規則性是造成兩者間呈現非線性變化關系的重要因素之一。

［1］Gurevich B，Lopatnikov S L.Velocity and attenuation of elastic waves in finely layered porous rocks［J］.Geophysical Journal International，1995，121(3):933 －947.

［2］Wright T W.Elastic wave propagation through a material with voids ［J］. J.Mech.Phys.Solids， 1998，46(10):2033－2047.

［3］鄧慶田，楊智春.導波在多損傷板結構中的散射［J］.振動與沖擊，2010，29(4):40 －44.

［4］楊宏峰，施行覺.軸向壓力下砂巖波速的實驗研究［J］.地球物理學進展，2004，19(2):481 －485.

［5］韓 嵩，蔡美峰.節理巖體物理模擬與超聲波試驗研究［J］.巖石力學與工程學報，2007，26(5):1026－1033.

［6］Oda M.Permeability tensor for discontinuous rock masses［J］.Geotechnique，1985，35(4):483 －495.

［7］Oda M.A crack tensor and its relation to wave velocity anisotropy in jointed rock masses［J］.International Journal of Rock Mechanics and Mining sciences and Geomechanics Abstracts，1986，23(6):387 －397.

［8］胡 剛，郝傳波，景海河.爆炸作用下巖石介質應力波傳播規律研究［J］.煤炭學報，2001，226(3):270－273.

［9］Ju Y，Sudak L，Xie H P.Study on stress wave propagation in fractured rocks with fractal joint surfaces［J］.International Journal of Solids and Structures，2007，44(13):4256 －4271.

［10］逯靜洲，林 皋，王 哲.混凝土經歷三向受壓荷載歷史后強度劣化及超聲波探傷方法的研究［J］.工程力學，2002，19(5):52 －57.

［11］劉天云，劉光廷.拱壩地震動隨機響應分析［J］.工程力學，2002，17(6):20 －25.

［12］朱金穎，陳龍珠，嚴細水.混凝土受力狀態下超聲波傳播特性研究［J］.工程力學，1998，15(3):111－117.

［13］Kahraman S.The effects of fracture roughness on P-wave velocity［J］.Engineering Geology，2002，63(3):347 －350.

［14］Mandelbrot B B.How long is the coastline of Britain?Statistical self-similarity and fractal dimension ［J］.Science，1967，155(4):636 －638.

［15］Clarke K C.Computation of the fractal dimension of topographic surfaces using the triangular prism surface area method［J］.Computers ＆ Geoscience，1986，12(5):713－722.

［16］Xie H P，Wang J A.Direct fractal measurement of fracture surfaces［J］.International Journal of Solids and Structures，1999，36(20):3073 －3084.

［17］Zhou H W，Xie H P.Direct estimation of the fractal dimensions of a fracture surface of rock ［J］.Surface Review and Letters，2003，10(5):751 －762.

［18］張亞衡，周宏偉，謝和平.粗糙表面分形維數估算的改進立方體覆蓋法［J］.巖石力學與工程學報，2005，24(17):3192－3196.

［19］彭瑞東，謝和平，鞠 楊.二維數字圖像分形維數的計算方法［J］.中國礦業大學學報，2004，33(1):19－24.

［20］Pentland A P.Fractal—based description of nature scenes［J］.IEEE Transactions on PAMI，1984，6(6):315 －326.