林景宜， 曾昭發(fā)， 李 靜， 張 領(lǐng)， 槐 楠， 胡志鵬

(1.吉林大學(xué) 國土資源部應(yīng)用地球物理重點實驗室，長春 130026；2.吉林大學(xué) 地球探測科學(xué)與技術(shù)學(xué)院，長春 130026)

0 引言

土壤和巖石顆粒具有明顯的塑性、脹縮性、吸濕性和黏結(jié)性，是影響收縮特征的重要因素，巖土的收縮能力與土壤黏粒含量呈顯著正相關(guān)[1]，在遇水和失水后，會在上述性質(zhì)差異較大的區(qū)域，產(chǎn)生有自相似性的裂隙分布，裂隙具有明顯的分維特征，且裂隙率、含水量、滲透系數(shù)和變形量之間具有一定的關(guān)系[2]；巖土中的裂隙帶是巖土中的軟弱結(jié)構(gòu)和水的優(yōu)先流動途徑[3]；坡體上裂隙的存在能夠?qū)е聝?yōu)先流的發(fā)生，使地表水很快到達(dá)土壤深層，增大了入滲面積，提高了入滲速度[4]，同時，裂隙受雨水滲透力及滲透侵蝕的作用，不斷地延伸、擴(kuò)大，最終形成滑坡[5]?？梢妼r土介質(zhì)中的裂隙帶進(jìn)行探測，是災(zāi)害預(yù)警和災(zāi)害評價重要的基礎(chǔ)。

構(gòu)成巖土的顆粒和孔隙大小不同，形狀各異，并且可能以各種方式連接，結(jié)構(gòu)形狀十分復(fù)雜，通過常規(guī)的歐式幾何很難模擬巖土內(nèi)部的結(jié)構(gòu)。但是巖土結(jié)構(gòu)狀況的一些參量，比如粒徑分布、裂隙發(fā)育程度、聯(lián)通狀況都表現(xiàn)出分形特征[6]，所以近年來，不少研究采用分形理論模擬巖土結(jié)構(gòu)，Perfect等[7]、Perrier等[8]運用QSF分形理論對含聚合物顆粒的碎裂土壤結(jié)構(gòu)進(jìn)行了定量化模擬；Kravchenko等[9]、Medina等[10]計算出了土壤粒分形模型來模擬土壤結(jié)構(gòu)和不同尺度下質(zhì)量構(gòu)成；Shepherd[11]用Koch曲線模擬孔隙通道連通狀況問題；Jiang等[12]建立基于混合型自相關(guān)函數(shù)的多尺度隨機(jī)介質(zhì)模型來模擬介質(zhì)的非均一性。然而，上述研究偏重于對巖土整體結(jié)構(gòu)的建模，未能模擬出巖土中軟弱帶的局部特征。

巖土介質(zhì)建模后，對模型進(jìn)行數(shù)值模擬的方法眾多，其中探地雷達(dá)是對淺層介質(zhì)探測的重要方法。探地雷達(dá)(GPR)數(shù)據(jù)的分析和解釋基于地下介質(zhì)的電性參數(shù)，主要包括介電常數(shù)和電導(dǎo)率，應(yīng)用在介質(zhì)界面和目標(biāo)形狀的探測方面[13]，在裂隙帶探測中也能發(fā)揮重要作用。Godio等[14]用探地雷達(dá)和其他方法綜合探測了含水層的裂隙帶和滲透層的幾何形態(tài)；鄧國文等[15]在隧道超前預(yù)測中用探地雷達(dá)對斷層破碎帶和裂隙帶進(jìn)行了探測和解釋。在GPR數(shù)值模擬方面，Li等[16]采用CFS-RIPML邊界導(dǎo)出高階FDTD方法，應(yīng)用到了GPR數(shù)值模擬中，有效降低了數(shù)值色散引起的誤差，提高了結(jié)果的準(zhǔn)確性。在結(jié)果處理和解釋過程中，時頻分析方法起到了重要作用，Stockwell等[17]在以Morlet小波為基本小波的基礎(chǔ)上提出了可逆的S變換(ST)時頻分析方法，S變換采用高斯窗函數(shù)，繼承了短時傅立葉變換(STFT)和小波變換(WT)的優(yōu)點，可以有效壓制噪聲，獲得更高分辨率的響應(yīng)結(jié)果；鄧攻等[18]使用S變換對深反射地震弱信號進(jìn)行了分解，提取出了被噪聲湮滅的低頻細(xì)節(jié)特征，提高了剖面的分辨率和同相軸的連續(xù)性。

筆者在前人研究的基礎(chǔ)上，使用三維Koch分形方法和Bruggeman等效介質(zhì)理論，來模擬真實裂隙帶分布，并使用隨機(jī)等效介質(zhì)理論建立巖土背景。建立一個三層背景下的裂隙帶模型，根據(jù)不同含水率計算裂隙的等效介電常數(shù)。對模型進(jìn)行GPR數(shù)值模擬，定性研究其波場特征，再對響應(yīng)進(jìn)行S變換得到時頻特征，定量擬合含水率與響應(yīng)的關(guān)系，最后進(jìn)行裂隙帶深度的誤差分析。

1 建模方法

本模型的構(gòu)建過程主要包括三維Koch曲線構(gòu)造、裂隙帶形態(tài)控制、裂隙帶位置控制、裂隙成分等效、巖土背景隨機(jī)化幾個步驟。

1.1 三維Koch曲線構(gòu)造

Koch曲線由瑞典數(shù)學(xué)家科赫(Koch)于1904年提出，是一種具有典型自相似性的折線，應(yīng)用三維Koch曲線進(jìn)行分形建模，可以更好地模擬出裂隙結(jié)構(gòu)的自相似性和聯(lián)通狀況。

圖1 滿足廣義能量極小值原理的三維Koch曲線Fig.1 Three-dimensional Koch curve of the principle of generalized energy minimal value

三維Koch曲線如圖1所示，取一條線段P1P5，將其三等分，得到點P2、P4，設(shè)P1(xp1,yp1,zp1)、P5(xp5,yp5,zp5)，則由定比分點公式(1)可得P2(xp2,yp2,zp2)、P4(xp4,yp4,zp4)。

(1)

在三維空間中，ΔP2P3P4可以圍繞軸P2P4任意旋轉(zhuǎn)，所以點P3(xp3,yp3,zp3)不唯一。在自然界中非線性系統(tǒng)屬于強(qiáng)迫耗散動力系統(tǒng)，由廣義能量極小值原理和熱力學(xué)第一、第二定律可知，這種系統(tǒng)在經(jīng)過足夠長的時間后，總是趨向一種不可逆過程最弱、系統(tǒng)廣義能量最小的有序狀態(tài)[19]。裂隙結(jié)構(gòu)屬于不規(guī)則的非線性復(fù)雜系統(tǒng)，也遵循廣義的能量最小值原理。

設(shè)三維Koch曲線整體為一剛體系統(tǒng)，則其彈性勢能為“0”，系統(tǒng)總能量等于重力勢能，在系統(tǒng)質(zhì)量一定時，重力勢能與系統(tǒng)重心位置有關(guān)，重心位置越低，則重力勢能越小，系統(tǒng)總能量也越小[19]。如圖1所示，當(dāng)點P3(xp3,yp3,zp3)在Z方向上最小時，重心位置最低，此時ΔP2P3P4所在平面垂直于平面XOY，過點做一垂線L，由幾何關(guān)系知，L平行于平面XOY。以L為軸線，將線段P2P4順時針旋轉(zhuǎn)60，旋轉(zhuǎn)后點P4即為點P3所在位置。

(2)

(3)

M=+cosθ·(I-)+sinθ·A*

(4)

P′=P·MT

(5)

式(2)～式(5)中:A=[ax,ay,az]為單位長旋轉(zhuǎn)軸;θ為旋轉(zhuǎn)角;I為三階單位矩陣;P為待旋轉(zhuǎn)點的坐標(biāo);P′為旋轉(zhuǎn)后坐標(biāo)。

(6)

可求得旋轉(zhuǎn)后的坐標(biāo)點P3，此時整個系統(tǒng)的總能量最小。將P1到P5的位置都求出后，每兩個相鄰的點之間再進(jìn)行相同計算過程，重復(fù)幾層以后，可遞歸出一條滿足廣義能量最小值定理的Koch曲線。

1.2 裂隙帶形態(tài)控制

在遇水和失水后，會在收縮性質(zhì)差異較大的區(qū)域，產(chǎn)生有眾多自相似性的裂隙分布，形成裂隙帶。

圖2 3維Koch分形裂隙的構(gòu)造過程Fig.2 The construction process of 3D Koch fractal fissure

如圖2所示，設(shè)裂隙深度上、下限為htop和hbot后，將裂隙沿Z軸分n層，共n+1個界面，每層的底界面即為下一層的頂界面，每段裂隙的底端即為下一段裂隙的頂端。設(shè)裂隙帶在頂界面的中心點位置為Ptop_mid)(Xtop_mid,Ytop_mid,htop)，底界面中心點位置為Pbot_mid(Xbot_mid,Ybot_mid,hbot)，裂隙帶傾角為α，則由幾何關(guān)系可得：

tanα=

(7)

當(dāng)Xtop_mid、Ytop_mid和α給定后，可求得一組Xbot_mid,Ybot_mid，確定了Ptop_mid和Pbot_mid的位置(圖2)。

(8)

式中:P1i_mid(x1i,y1i,z1i)為第i層裂隙帶的頂界面中心點坐標(biāo)，P2i_mid(x2i,y2i,z2i)為第i層裂隙帶底界面中心點坐標(biāo)。

設(shè)裂隙帶的長和寬為dx和dy，為了裂隙帶在內(nèi)部體現(xiàn)出每條裂隙位置的隨機(jī)性，在(xni mid±dx/2,yni mid±dy/2)范圍內(nèi)隨機(jī)取值，即可得到第i層裂隙的頂端點坐標(biāo)(x1i,y1i,z1i)和底端點坐標(biāo)(x2i,y2i,z2i)通過式(7)、式(8)可畫出第i層的分形裂隙，將各層連接即為一條裂隙，多次遞歸上部步驟即可畫出一個具有傾角的裂隙帶。

1.3 裂隙帶位置控制

經(jīng)過處理后得到固定位置處的裂隙帶，還需將裂隙帶旋轉(zhuǎn)平移到指定位置。將裂隙帶各個轉(zhuǎn)折點位置記為一個矩陣[X,Y,Z]，X、Y、Z均為列向量，

(9)

(10)

式中:β為裂隙帶在xoy面上的旋轉(zhuǎn)角度;R為二維旋轉(zhuǎn)矩陣;[Dx,Dy]表示裂隙帶在x和y方向上平移;[Xnew,Ynew]為旋轉(zhuǎn)平移變換后的，在z方向上，裂隙帶位置沒有發(fā)生變化;Znew與Z相同。如圖3所示，先旋轉(zhuǎn)β改變裂隙帶傾斜方向，再平移D改變裂隙帶位置?？傻昧严稁У男D(zhuǎn)和平移后的位置[Xnew,Ynew,Znew]。

圖3 裂隙帶的旋轉(zhuǎn)平移過程Fig.3 The rotation process of the fissure belt

1.4 應(yīng)用等效介質(zhì)方法

經(jīng)過上幾步處理，實際得到的是裂隙帶的轉(zhuǎn)折點骨架，還要給裂隙帶賦予實際的物理意義。將模型中的每一條裂隙看作是周圍眾多細(xì)小裂隙的等效，不僅能突出需要研究的尺度，還能提高計算效率，應(yīng)用Bruggeman等效介質(zhì)理論：

(11)

式中：εeff為Bruggeman等效介電常數(shù)；εi為第i類介質(zhì)的介電常數(shù)，fi為第i類介質(zhì)的體積分?jǐn)?shù)，代入裂隙內(nèi)空氣和裂隙周圍巖土的介電常數(shù)和體積分?jǐn)?shù)，即可得到等效介電常數(shù)εeff，為了更好地體現(xiàn)出隨機(jī)效果，將εeff加上隨機(jī)變化δε(δε滿足高斯分布)：

εfra=εeff+δε

(12)

式中:εfra即為裂隙的等效介電常數(shù)，賦予介電常數(shù)后繪制成的裂隙帶(圖4),含有35條裂隙，分四層，傾角為45°,εfra=5。

圖4 Bruggeman等效處理后的裂隙帶Fig.4 Fissure belt after satisfying Bruggeman equivalent theory

1.5 裂隙帶巖土背景隨機(jī)化處理

將復(fù)雜非均勻介質(zhì)中的非均勻體看成一個空間隨機(jī)過程，用統(tǒng)計學(xué)方法描述介質(zhì)非均勻性所形成的非均勻分布就是隨機(jī)介質(zhì)模型[13]。由于土壤顆粒的組成及排列方式復(fù)雜，還有水分、有機(jī)質(zhì)和礦物質(zhì)等原因的影響，所以隨機(jī)介質(zhì)模型可以更準(zhǔn)確地描述實際巖土背景的隨機(jī)非均勻性分布，即：隨機(jī)介質(zhì)模型m(k)為均勻介質(zhì)mi(k)和隨機(jī)擾動所δm(k)組成(δm(k)滿足高斯分布)，描述隨機(jī)介質(zhì)模型的基本思路描述如下：

m(k)=mi(k)+δm(k)=mi(k)[1+f(k) ]

(13)

式中：三維空間位置矢量表示為k=(x,y,z)；用f(k)代表這種隨機(jī)介質(zhì)的相對擾動的特征，其中f(k)=δm(k)/mi(k)。

對于f(k)的選擇，采用了小尺度非均勻性的橢球式混合型自相關(guān)函數(shù)：

(14)

式中：r是代表模糊度因子(當(dāng)r=0時為高斯型橢圓自相關(guān)函數(shù)，當(dāng)r=1時為指數(shù)型橢圓自相關(guān)函數(shù)，r介于兩者之間的，則為混合型自相關(guān)函數(shù))；a、b、c分別代表x、y、z方向的自相關(guān)長度；再用式(15)作傅里葉變換，得到能量譜密度函數(shù)：

R′ (k)=‖F(xiàn)(k)‖2

(15)

式中:k=(kx,ky,Kkz)再加入隨機(jī)相位φ(k)，變換后得到功率譜函數(shù)：

(16)

對公式(16)中F′ (k)進(jìn)行傅立葉逆變換，從波數(shù)域變換到空間域的R(x,y,z),得到隨機(jī)等效背景模型(圖 5),均勻介質(zhì)mi(k)=5，x、y、z方向的自相關(guān)長度A=5、B=5、C=5,模糊度r=0.5。

圖5 隨機(jī)等效處理后的背景介質(zhì)Fig.5 The medium after satisfying the stochastic equivalent media theory

2 模型實例正演

圖6所示為一個實例模型，模型大小為5 m×5 m×5 m，網(wǎng)格數(shù)為200×200×200，空間網(wǎng)格步長為0.025。模型背景分三層，第一層模擬表層土壤，厚2 m，介電常數(shù)平均值為10；第二層模擬濕潤粘土層，厚2 m，介電常數(shù)平均值為15；第三層模擬基巖，厚1 m，介電常數(shù)平均值為4。所有層面背景自相關(guān)長度A=5,B=5,C=5,模糊度r=0.5，背景電導(dǎo)率為0.01 s/m。第二層內(nèi)的裂隙帶大小為5 m×1.25 m×2 m，傾角為38°，裂隙占裂隙帶體積的20%，背景介質(zhì)占裂隙帶體積的80%。

圖6 在隨機(jī)背景中的分形等效裂隙帶模型Fig.6 Fractal equivalent fracture belt model in random background

表1 不同含水率下裂隙帶的Bruggeman等效介電常數(shù)

從0%到20%取5種含水率，計算裂隙的Bruggeman等效介電常數(shù)，如表 1所示。將5種等效介電常數(shù)賦值給裂隙，并對其進(jìn)行探地雷達(dá)正演模擬。

探地雷達(dá)(GPR)波脈沖激勵源的中心頻率為100 MHz，時間步長為0.05 ns，時窗長度為120 ns。正演模擬得到響應(yīng)結(jié)果后去除直達(dá)波，如圖7所示，

圖7 不同含水率模型的探地雷達(dá)模擬響應(yīng)結(jié)果Fig.7 Numerical simulation of GPR with different moisture content

左側(cè)圖7(a)～圖7(e)依次為0%、5%、10%、15%、20%含水率的探地雷達(dá)正演響應(yīng)結(jié)果。在圖7(a)～圖7(e)中均可觀察到隨機(jī)背景造成的微弱擾動，除圖7(e)(含水率20%)外，分界面均較為清晰；在裂隙帶區(qū)域，由于內(nèi)部介質(zhì)比較復(fù)雜，電磁波傳播時的經(jīng)歷多次反射和繞射過程，雙曲線波形較為散亂，同相軸更為模糊，但是仍能分辨出裂隙帶的位置和形狀；由于第二層介質(zhì)介電常數(shù)為15，對比圖7(a)～圖7(e)可發(fā)現(xiàn)，含水率為10%時(圖7(c)，左)，介電常數(shù)與背景較為接近，裂隙帶造成的波形變化較??；裂隙帶與背景的介電常數(shù)相差越多，波形擾動就越大，尤其當(dāng)含水率為20%時(圖7(e)，左)，也就是裂隙內(nèi)充滿水的情況下，裂隙帶造成的混亂波形掩蓋了第二層分界面。

取2 m位置處單道記錄進(jìn)行S變換，得到5種含水率響應(yīng)的時頻圖，如圖7所示，右側(cè)圖7(a)～圖7(e)依次為0%、5%、10%、15%、20%含水率響應(yīng)結(jié)果的S變換時頻圖。由于S變換壓制了噪聲，減弱了隨機(jī)背景的影響，使得兩層分界面和裂隙帶的異常均清晰可見。對比圖7中圖7(a)～圖7(e)，可以發(fā)現(xiàn)含水率為10%時(圖7(c)，右)，由于裂隙帶介電常數(shù)與背景較為接近，兩層分界面的異常明顯強(qiáng)于裂隙帶異常；含水率變低或變高，即介電常數(shù)與背景差異越大，裂隙帶異常就越明顯；且含水率較高時(圖7(d)和圖7(e)，右)，裂隙帶異常明顯強(qiáng)于分界面異常值。

不同含水率模型響應(yīng)的2 m位置處單道記錄(圖7(f))也佐證了上述規(guī)律。在前50 ns，不同含水率模型響應(yīng)基本一致；在含水率為10%時，裂隙帶位置(68 ns附近)振幅變化相對較小，而含水率越低或越高，振幅變化越大，異常越明顯。

由Bruggeman等效介質(zhì)理論計算可得出介電常數(shù)與第二層背景相等(ε=15)時，含水率為8.63%，此時裂隙帶基本無響應(yīng)。將不同含水率和其時頻圖中裂隙帶處的最大幅值進(jìn)行三次多項式擬合處理，得到的擬合曲線和公式如圖8和公式17所示。通過時頻變換，我們得到了該模型下的裂隙帶含水率與正演響應(yīng)的定量關(guān)系，根據(jù)此擬合公式，可用響應(yīng)估算含水率。

圖8 不同含水率時頻圖在裂隙帶位置最大幅值的曲線擬合Fig.8 Curve fitting of the maximum amplitude of the time - frequency pattern with different water content

(17)

根據(jù)不同含水率時頻圖(圖7(a)～圖7(e)，右)中裂隙帶異常最大值位置的走時，可以計算出該道裂隙帶中心深度，誤差估計如表2所示，可以發(fā)現(xiàn)，響應(yīng)結(jié)果對裂隙帶深度的估計較為準(zhǔn)確。

表2 不同含水率模型響應(yīng)結(jié)果的裂隙帶深度誤差估計Tab. 2 The errors estimation of fissure belt for response results of different moisture model

3 結(jié)論

三維Koch分形曲線具有聯(lián)通性、隨機(jī)性和自相似性，通過Koch曲線模擬的骨架表現(xiàn)了裂隙分布的分維特征；對裂隙中空氣和水應(yīng)用Bruggenman等效介質(zhì)方法，得到多種含水率的巖土介電常數(shù)，滿足了不同情況和計算效率的需求；應(yīng)用隨機(jī)等效介質(zhì)方法使巖土背景呈現(xiàn)多尺度的隨機(jī)性，模擬了真實巖土介質(zhì)的復(fù)雜狀況。

對不同含水率下的裂隙帶模型進(jìn)行GPR數(shù)值模擬后，可以發(fā)現(xiàn)裂隙帶含水率對響應(yīng)波形的影響較大。裂隙等效介電常數(shù)與背景相差越多，裂隙帶的響應(yīng)越明顯，當(dāng)裂隙內(nèi)充滿水后，噪聲甚至可以掩蓋分界面波形。對響應(yīng)進(jìn)行S變換(ST)后分析其時頻特征，既壓制了混雜的噪聲，又可擬合出該模型下的含水率-響應(yīng)曲線。根據(jù)此曲線的定量關(guān)系，可通過響應(yīng)估算含水率。最后進(jìn)行的裂隙帶深度的誤差分析,證實了數(shù)值模擬及時頻分析結(jié)果的可信性。此建模、探地雷達(dá)響應(yīng)數(shù)值模擬及信號時頻分析方法可以為巖土體中軟弱帶及含水裂隙帶的探測提供參考。