三維隨機(jī)裂隙介質(zhì)建模與GPR響應(yīng)計(jì)算

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

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

0 引言

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

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

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

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

1 建模方法

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

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，將其三等分，得到點(diǎn)P2、P4，設(shè)P1(xp1,yp1,zp1)、P5(xp5,yp5,zp5)，則由定比分點(diǎn)公式(1)可得P2(xp2,yp2,zp2)、P4(xp4,yp4,zp4)。

(1)

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

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

(2)

(3)

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

(4)

P′=P·MT

(5)

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

(6)

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

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

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

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

如圖2所示，設(shè)裂隙深度上、下限為htop和hbot后，將裂隙沿Z軸分n層，共n+1個(gè)界面，每層的底界面即為下一層的頂界面，每段裂隙的底端即為下一段裂隙的頂端。設(shè)裂隙帶在頂界面的中心點(diǎn)位置為Ptop_mid)(Xtop_mid,Ytop_mid,htop)，底界面中心點(diǎn)位置為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層裂隙帶的頂界面中心點(diǎn)坐標(biāo)，P2i_mid(x2i,y2i,z2i)為第i層裂隙帶底界面中心點(diǎn)坐標(biāo)。

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

1.3 裂隙帶位置控制

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

(9)

(10)

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

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

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

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

(11)

式中：εeff為Bruggeman等效介電常數(shù)；εi為第i類(lèi)介質(zhì)的介電常數(shù)，fi為第i類(lèi)介質(zhì)的體積分?jǐn)?shù)，代入裂隙內(nèi)空氣和裂隙周?chē)鷰r土的介電常數(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ì)中的非均勻體看成一個(gè)空間隨機(jī)過(guò)程，用統(tǒng)計(jì)學(xué)方法描述介質(zhì)非均勻性所形成的非均勻分布就是隨機(jī)介質(zhì)模型[13]。由于土壤顆粒的組成及排列方式復(fù)雜，還有水分、有機(jī)質(zhì)和礦物質(zhì)等原因的影響，所以隨機(jī)介質(zhì)模型可以更準(zhǔn)確地描述實(shí)際巖土背景的隨機(jī)非均勻性分布，即：隨機(jī)介質(zhì)模型m(k)為均勻介質(zhì)mi(k)和隨機(jī)擾動(dòng)所δ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ì)的相對(duì)擾動(dòng)的特征，其中f(k)=δm(k)/mi(k)。

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

(14)

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

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

(15)

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

(16)

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

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

2 模型實(shí)例正演

圖6所示為一個(gè)實(shí)例模型，模型大小為5 m×5 m×5 m，網(wǎng)格數(shù)為200×200×200，空間網(wǎng)格步長(zhǎng)為0.025。模型背景分三層，第一層模擬表層土壤，厚2 m，介電常數(shù)平均值為10；第二層模擬濕潤(rùn)粘土層，厚2 m，介電常數(shù)平均值為15；第三層模擬基巖，厚1 m，介電常數(shù)平均值為4。所有層面背景自相關(guān)長(zhǎng)度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種含水率，計(jì)算裂隙的Bruggeman等效介電常數(shù)，如表 1所示。將5種等效介電常數(shù)賦值給裂隙，并對(duì)其進(jìn)行探地雷達(dá)正演模擬。

探地雷達(dá)(GPR)波脈沖激勵(lì)源的中心頻率為100 MHz，時(shí)間步長(zhǎng)為0.05 ns，時(shí)窗長(zhǎng)度為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ī)背景造成的微弱擾動(dòng)，除圖7(e)(含水率20%)外，分界面均較為清晰；在裂隙帶區(qū)域，由于內(nèi)部介質(zhì)比較復(fù)雜，電磁波傳播時(shí)的經(jīng)歷多次反射和繞射過(guò)程，雙曲線波形較為散亂，同相軸更為模糊，但是仍能分辨出裂隙帶的位置和形狀；由于第二層介質(zhì)介電常數(shù)為15，對(duì)比圖7(a)～圖7(e)可發(fā)現(xiàn)，含水率為10%時(shí)(圖7(c)，左)，介電常數(shù)與背景較為接近，裂隙帶造成的波形變化較??；裂隙帶與背景的介電常數(shù)相差越多，波形擾動(dòng)就越大，尤其當(dāng)含水率為20%時(shí)(圖7(e)，左)，也就是裂隙內(nèi)充滿水的情況下，裂隙帶造成的混亂波形掩蓋了第二層分界面。

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

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

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

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

(17)

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

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

3 結(jié)論

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

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