各向異性介質(zhì)中地震波參數(shù)反演及其在隧道不良地質(zhì)預(yù)報(bào)中的應(yīng)用

何現(xiàn)啟，彭凌星，魯光銀，朱自強(qiáng)，高峰

(1.湖南省交通規(guī)劃勘察設(shè)計(jì)院有限公司，湖南長(zhǎng)沙，410200；2.中南大學(xué)地球科學(xué)與信息物理學(xué)院，湖南長(zhǎng)沙，410083)

EDA(extensive dliatancy anisotropy)介質(zhì)即具有一組定向排列的垂直裂隙的介質(zhì)，也稱廣泛各向異性介質(zhì)、裂隙擴(kuò)容各向異性介質(zhì)和方位各向異性介質(zhì)。通過對(duì)EDA 介質(zhì)的研究，可探明地層裂隙的分布和密度，進(jìn)一步研究裂隙中充填物的性質(zhì)，為災(zāi)害預(yù)報(bào)、地下水探測(cè)等提供重要資料。通過對(duì)地震波速度及其衰減進(jìn)行分析還可以推斷地下冷熱水的運(yùn)移和儲(chǔ)存情況，從而評(píng)價(jià)地下水資源，預(yù)報(bào)地下工程突水、涌泥等地質(zhì)災(zāi)害。地震波振幅隨巖石物理性質(zhì)的變化比地震波速度的變化更大，因此，地震波衰減可能比速度攜帶了更多的表征巖石物理性質(zhì)的信息。在各向異性速度反演及走時(shí)研究方面，盧明輝等[1]采用小生境遺傳算法，對(duì)3條成一定角度的測(cè)線走時(shí)信息進(jìn)行速度和各向異性參數(shù)反演。杜啟振等[2]對(duì)常規(guī)雙曲線時(shí)距方程進(jìn)行了改進(jìn)，給出了VTI(vertical transverse isotropy)介質(zhì)時(shí)距方程。在各向異性地震波衰減估計(jì)方面，BEHURA等[3]研究了VTI介質(zhì)和正交各向異性介質(zhì)中相衰減和群衰減的變化規(guī)律，并利用譜比法進(jìn)行衰減估計(jì)；徐善輝等[4]通過PP 轉(zhuǎn)換縱波和PS 轉(zhuǎn)換橫波的反射系數(shù)反演了各向異性參數(shù)和對(duì)稱軸傾角、方位角等，提出了利用PS 波振幅定性分析TTI(titled transverse isotropy)介質(zhì)對(duì)稱軸傾角的方法。在射線追蹤算法方面，李強(qiáng)等[5]系統(tǒng)綜述了國(guó)內(nèi)外不均勻介質(zhì)中各種主要和實(shí)用的射線追蹤方法。段鵬飛等[6]為了方便、快捷地實(shí)現(xiàn)TI(transverse isotropy)介質(zhì)射線走時(shí)與局部角度的計(jì)算，討論和對(duì)比了2種改進(jìn)的射線追蹤方法。曲英銘等[7]在Vidale 差分算式基礎(chǔ)上推導(dǎo)出不規(guī)則網(wǎng)格差分算式，是適用于不同網(wǎng)格剖分方式的有限差分走時(shí)計(jì)算方法。王東鶴等[8]分析了各種算法的研究現(xiàn)狀，并對(duì)射線追蹤法的發(fā)展趨勢(shì)進(jìn)行了展望。韓松[9]通過1.2D 和3D 的快速掃描法計(jì)算了TTI介質(zhì)中qP 波的走時(shí)，并通過新型快速掃描法(fast scanning method,FSM)，得到走時(shí)。歐陽(yáng)進(jìn)[10]采用射線方法對(duì)復(fù)雜的地質(zhì)模型進(jìn)行了波場(chǎng)正演。HE 等[11]指出含裂縫介質(zhì)不僅具有速度各向異性，而且具有顯著的衰減各向異性，衰減各向異性分析方法可應(yīng)用于裂隙產(chǎn)狀、裂隙密度和發(fā)育范圍等的估計(jì)。居興國(guó)等[12]在程函方程法射線追蹤的基礎(chǔ)上，提出了基于相速度的TTI介質(zhì)射線追蹤新策略。劉瑞合等[13]采用PSV 波分離的射線追蹤方法，為層析反演提供了準(zhǔn)確的旅行時(shí)和射線路徑等信息。龔屹等[14]對(duì)采用射線追蹤方法得到的微地震事件走時(shí)正演值與實(shí)際走時(shí)進(jìn)行了對(duì)比分析。張敏等[15]通過對(duì)運(yùn)動(dòng)學(xué)和動(dòng)力學(xué)追蹤方程進(jìn)行了修改和簡(jiǎn)化，有效地提高了各向異性介質(zhì)射線追蹤算法的計(jì)算效率。洪啟宇等[16]指出利用2個(gè)相互正交的變井源距垂直地震剖面(walkaway vertical seismic profiling)可以完全確定鉆井中TTI 介質(zhì)qP 波的3 個(gè)WA(weak anisotropic)參數(shù)和對(duì)稱軸的2 個(gè)方向參數(shù)。葛子建等[17]基于貝葉斯反演框架建立了P波線性AVAZ(amplitwde versus angle and azimuth)反演方法。國(guó)內(nèi)外研究者在各向異性介質(zhì)參數(shù)反演方面也進(jìn)行了大量研究，如RUSMANUGROHO等[18]發(fā)現(xiàn)Voigt和Tromp參數(shù)能快速、穩(wěn)定地反演陡峭傾斜的背斜結(jié)構(gòu)。XU 等[19]建立了相似性分析框架用以反演最佳各向異性參數(shù)。LU 等[20]為了增強(qiáng)反演結(jié)果，采用層次反演策略來(lái)解決Gauss-Newton 方法中的局部最小解問題，并分別在理想和噪聲的情況下對(duì)合成數(shù)據(jù)進(jìn)行了測(cè)試，發(fā)現(xiàn)該方法可以成功識(shí)別復(fù)雜巖性和流體信息。BOITZ 等[21]研究了橫向各向同性(VTI)介質(zhì)的2種不同效應(yīng)，即源的各向異性和傳播路徑上的各向異性，并用效能張量乘以具有等效剪切模量的各向同性彈性張量，來(lái)探明與巖石各向異性無(wú)關(guān)的構(gòu)造變形情況。HADDEN 等[22]開發(fā)了一種各向異性波形斷層掃描(anisotropic wave tomographe,AWT)方法，利用聲學(xué)近似來(lái)模擬VTI介質(zhì)中的地震波，并將這種方法應(yīng)用于交叉井?dāng)?shù)據(jù)處理，在應(yīng)用于加拿大西部沉積環(huán)境中的井間數(shù)據(jù)處理中具有很高的分辨率。人們對(duì)地震各向異性進(jìn)行了大量研究，但將相關(guān)方法用于裂隙范圍探測(cè)、地下水預(yù)測(cè)及誤差分析尤其是隧道富水區(qū)等地質(zhì)災(zāi)害探測(cè)的研究很少。為此，本文作者在研究彈性HTI(horizontal transverse isotropy)的基礎(chǔ)上，研究黏彈性EDA 介質(zhì)地震波速度分析算法并對(duì)各向異性參數(shù)、裂隙方位、裂隙密度等參數(shù)進(jìn)行反演，對(duì)反演誤差及其影響因素進(jìn)行分析。

1 EDA 介質(zhì)中傾斜界面的動(dòng)校正速度及時(shí)距方程

1.1 方位垂直裂隙(EDA)介質(zhì)中傾斜界面的動(dòng)校正速度

由HTI介質(zhì)精確動(dòng)校正速度表達(dá)式，經(jīng)坐標(biāo)旋轉(zhuǎn)可得到垂直非均勻EDA 介質(zhì)中傾斜界面的精確動(dòng)校正速度[23-25]：

式中：α為共中心點(diǎn)排列的方位角；φ為EDA介質(zhì)的對(duì)稱軸的方位角；?為對(duì)稱軸傾角；vnmo為動(dòng)校正速度。

在弱各向異性介質(zhì)中，同理可得沿界面傾向的P波的動(dòng)校正速度為[25]

由式(2)可知?jiǎng)有Ｕ俣葀nmo(φ,?)主要由垂直各向異性參數(shù)ε(V)與δ(V)的差值和水平界面的動(dòng)校正速度vnmo(0)決定。由HTI介質(zhì)中的動(dòng)校正公式，經(jīng)坐標(biāo)旋轉(zhuǎn)同樣可得沿界面走向的動(dòng)校正速度vnmo((π/2)+φ,?)為

1.2 對(duì)稱面內(nèi)P波動(dòng)校正速度與射線參數(shù)的關(guān)系

由HTI 介質(zhì)中射線動(dòng)校正速度推導(dǎo)[25]，EDA 介質(zhì)中沿界面傾向的動(dòng)校正速度可用射線參數(shù)表示為

式中：q'為一階導(dǎo)數(shù)；q''為二階導(dǎo)數(shù)；p為慢度向量的水平分量；q為慢度向量的垂直分量，

式中：θ=α-φ。

沿界面走向的動(dòng)校正速度可用射線參數(shù)表示為

與HTI介質(zhì)一樣，所有的參數(shù)都為零偏移距射線參數(shù)。動(dòng)校正速度與射線參數(shù)的關(guān)系都隱藏在慢度向量分量和其微商中，為確定動(dòng)校正系數(shù)，便于反演計(jì)算，下面主要討論動(dòng)校正速度與射線參數(shù)的顯性關(guān)系。

2 EDA 介質(zhì)中三維動(dòng)校正速度與射線參數(shù)的關(guān)系

2.1 均勻EDA介質(zhì)中用射線參數(shù)表示的動(dòng)校正速度精確表達(dá)式

由動(dòng)校正速度和慢度向量及2層介質(zhì)中的時(shí)距方程可得[24-25]

式(7)對(duì)所有的非轉(zhuǎn)換波都成立，其中p1和p2分別為零偏移距射線在x和y方向的2 個(gè)水平慢度向量分量，垂直慢度向量分量q=p3，也可由Christoffel方程計(jì)算。

2.2 弱各向異性EDA介質(zhì)中對(duì)稱面內(nèi)射線參數(shù)與動(dòng)校正速度的關(guān)系

EDA 介質(zhì)可由HTI 介質(zhì)繞z旋轉(zhuǎn)角度φ得到，對(duì)垂直方向的參數(shù)沒有影響，因此，計(jì)算EDA 介質(zhì)的垂直慢度仍可采用HTI介質(zhì)中的計(jì)算公式。

將垂直慢度q用η(V)表示，并將q及其微商代入動(dòng)校正方程(7)并用各向異性參數(shù)進(jìn)一步線性化，可得到弱各向異性EDA介質(zhì)中P波的近似動(dòng)校正速度：

其中：vP0為縱波速度；

假設(shè)傾斜方位與介質(zhì)對(duì)稱軸重合,可得界面(α=φ)的傾斜面上的動(dòng)校正速度為[22]

若用代替vP,nmo，EDA 介質(zhì)中的η(V)和分別與HTI介質(zhì)中的η(V)和相等，則式(10)與HTI介質(zhì)中沿傾向的P波動(dòng)校正速度的弱各向異性近似表達(dá)式一樣，vnmo(φ,p1)可由HTI 介質(zhì)中公式計(jì)算。若界面的傾向與各向同性面重合，走向與對(duì)稱面平行，則[22]

3 層狀EDA 介質(zhì)中傾斜界面的似Dix公式

在傾斜EDA 介質(zhì)中，對(duì)動(dòng)校正速度除考慮其測(cè)線方位角外，還要考慮傾角的影響。綜合分析HTI介質(zhì)對(duì)稱面內(nèi)的動(dòng)校正速度公式以及其與坐標(biāo)軸的相互關(guān)系，采用剝層法由相似分析得到[22]：

也可用下式計(jì)算：

其中：f(i)為第i層界面反射波振幅；t0(i)為第i層界面反射波走時(shí)。由式(12)，(13)和(14)即可得到其相關(guān)的等效速度和地層速度。

4 EDA介質(zhì)中的P波參數(shù)反演算法

對(duì)于1個(gè)給定的垂直時(shí)間t0，通過對(duì)水平速度vhor和動(dòng)校正速度vnmo進(jìn)行二維相似分析可以得到1 個(gè)具有最大相似值的模型，相似系數(shù)可由下式計(jì)算：

其中：M為道數(shù)。反射波振幅F(x,t)以及其平方F2(x,t)都是沿著時(shí)距曲線t(t'0,vnmo,x)計(jì)算的(其中雙曲線的頂點(diǎn)位于垂直走時(shí)t'0處，t'0在以時(shí)間t0為中心的平移窗口內(nèi))。2 個(gè)采樣時(shí)間間隔之間的反射波振幅F(x,t)可通過線性插值進(jìn)行計(jì)算。

通過對(duì)多個(gè)排列進(jìn)行速度分析即可確定動(dòng)校正橢圓，由動(dòng)校正速度與各向異性參數(shù)的關(guān)系式即可計(jì)算出各向異性參數(shù)[25]。

所采用的目標(biāo)函數(shù)如下：

其中：(x,n,m)為模型計(jì)算速度；(y,n)為數(shù)據(jù)反演所得速度；x和y為坐標(biāo)軸；m為模型參數(shù)。采用最小二乘法計(jì)算最佳層參數(shù)。

反射界面可用多項(xiàng)式表示，主要參數(shù)為節(jié)點(diǎn)數(shù)、節(jié)點(diǎn)處界面的法向量及x和y方向多項(xiàng)式的維數(shù)。其中，界面法向分量與射線參數(shù)的垂直分量一致，由基于擬牛頓法的最優(yōu)化射線追蹤法確定界面最佳幾何參數(shù)。

5 模型算例分析

5.1 地層模型

模型參數(shù)見表1，模型1示意圖見圖1。其中，第1層界面傾斜10°，第2層和第3層傾斜20°，上面3層為EDA介質(zhì)，下面2層為VTI介質(zhì)。

表1 模型1 彈性參數(shù)Table1 Elastic parameter of Model 1

圖1 模型1示意圖Fig.1 Schematic diagram of model 1

用黏彈性實(shí)射線追蹤法實(shí)現(xiàn)黏彈性EDA 中地震波的數(shù)值模擬，模擬地震波形見圖2。

5.2 速度分析

采用式(16)對(duì)速度進(jìn)行分析，結(jié)果見圖3。從圖3可以看出：0°測(cè)線的疊加速度顯示反演結(jié)果與模型值基本一致，相對(duì)誤差在5%以內(nèi)，動(dòng)校正結(jié)果顯示的界面層數(shù)及界面產(chǎn)狀與實(shí)際情況相符；同時(shí)，0°，10°，20°和30°測(cè)線的疊加速度和動(dòng)校正速度結(jié)果均與實(shí)際情況一致，說(shuō)明速度及動(dòng)校正公式計(jì)算結(jié)果準(zhǔn)確度高。

圖2 模型1 不同方位角的模擬地震波形Fig.2 Synthetic wave shapes of spreads with different azimuth angles for Model 1

圖3 0°測(cè)線的疊加速度和動(dòng)校正速度Fig.3 Stack velocity and NMO velocity for spread with azimuth angle φ=0°

5.3 反演結(jié)果與誤差分析

5.3.1 無(wú)噪聲反演及誤差分析

未添加噪聲的反演界面深度及傾角見圖4。由圖4可知：界面的反射傾角及方位與實(shí)際模型的基本一致，傾角誤差在3°以內(nèi)，方位角反演誤差在5°以內(nèi)。

其他相關(guān)參數(shù)反演結(jié)果見圖5～7，其中，實(shí)線表示參數(shù)的精確值，黑點(diǎn)表示反演值。從圖5～7 可見：迭代300次時(shí)，計(jì)算精度較高收斂較快，誤差分析結(jié)果顯示除方位角誤差稍大外，其余參數(shù)相對(duì)誤差都非常小，均在5%以內(nèi)。

圖4 未加噪聲的模型幾何參數(shù)反演結(jié)果Fig.4 Inversion results of geometric parameter of reflectors without noise

從圖5可見：彈性介質(zhì)各向異性參數(shù)ε及縱波速度反演計(jì)算精度較高，收斂較快，各層參數(shù)的反演值均集中在實(shí)線附近，誤差分析結(jié)果顯示相對(duì)誤差在5%以內(nèi)。

從圖6可見：對(duì)第1 層介質(zhì)方位角反演計(jì)算精度較高，收斂較快，對(duì)第2層和第3層介質(zhì)的反演值相對(duì)較離散，個(gè)別反演值相對(duì)誤差較大。

從圖7可見：對(duì)第1 層和第2 層介質(zhì)密度反演計(jì)算精度較高、收斂較快，對(duì)第3層介質(zhì)的反演值較離散，但反演值均集中在實(shí)線附近，相對(duì)誤差都在5%以內(nèi)。

5.3.2 帶噪聲反演成果及誤差分析

在反演中加入高斯噪聲，加噪后反演的介質(zhì)模型見圖8，反演結(jié)果見圖9～12。從圖9～12 可以看出：加入噪聲后對(duì)反演的速度和精度影響較大，其中各向異性參數(shù)ε和裂隙方位的收斂速度和精度較低，第3層的裂隙方位角反演相對(duì)誤差最大達(dá)10%。其余反演結(jié)果相對(duì)誤差均勻在5%以內(nèi)。

從圖8可見：加入高斯噪聲后，反演的地層界面傾角均與實(shí)際模型的基本一致，噪聲對(duì)反演的結(jié)果影響較小，可忽略。

圖5 各向異性ε及縱波速度反演結(jié)果Fig.5 Inversion results of ε and P wave velocity

圖6 裂隙方位角反演結(jié)果Fig.6 Inversion results of azimuth angle of fracture

圖7 密度反演結(jié)果Fig.7 Inversion results of density

從圖9可見：加入高斯噪聲后，ε反演結(jié)果較離散，速度反演值在實(shí)線附近，相對(duì)誤差在8%以內(nèi)。

從圖10可見：第1層介質(zhì)的方位角反演值在實(shí)線附近，相對(duì)誤差在5%以內(nèi)；第2 層介質(zhì)反演相對(duì)誤差在8%以內(nèi)；第3 層介質(zhì)反演結(jié)果較離散，反演誤差較大。

從圖11可見：第1層和第2層介質(zhì)的密度反演值均在實(shí)線附近，相對(duì)誤差在8%以內(nèi)；第3 層密度反演較離散，反演相對(duì)誤差較大。從圖12可見：第1，2和3層介質(zhì)的各向異性參數(shù)δ反演值均在實(shí)線附近，反演相對(duì)誤差在8%以內(nèi)。

圖9 ε及縱波速度反演結(jié)果Fig.9 Inversion results of ε and P wave velocity

圖8 加噪后反演的介質(zhì)模型幾何參數(shù)Fig.8 Inversion results for geometric parameter of reflectors adding gauss noise

圖10 裂隙方位角反演結(jié)果Fig.10 Inversion results of azimuth angle of fracture

圖11 密度反演結(jié)果Fig.11 Inversion results of density

圖12 各向異性參數(shù)δ反演結(jié)果Fig.12 Inversion results of anisotropy parameter δ

6 工程實(shí)例

為了驗(yàn)證該方法的可行性，將其應(yīng)用于裂隙發(fā)育區(qū)及富水區(qū)預(yù)報(bào)。通過現(xiàn)場(chǎng)測(cè)試，分析地層裂隙走向、節(jié)理密度與速度分布的關(guān)系，研究地震波速度及衰減與地下水賦存狀態(tài)的關(guān)系。

6.1 現(xiàn)場(chǎng)測(cè)試概況

測(cè)區(qū)樁號(hào)為K74+650—YK74+900，隧道走向約340°。依據(jù)洞內(nèi)開挖情況，地下巖體垂直向裂隙極其發(fā)育，地下水極其豐富，為了探明掌子面前方地質(zhì)情況，在地面上共布置4條不同方位的測(cè)線，采集廣角反射法地震勘探數(shù)據(jù)，測(cè)線布置如圖13所示。沿隧道的軸線建立采集坐標(biāo)系，沿軸線設(shè)為x，正方向指向大里程。為確保測(cè)線為長(zhǎng)排列，依據(jù)隧道埋深確定測(cè)線長(zhǎng)度大于270 m，實(shí)際測(cè)線長(zhǎng)度為360 m，測(cè)線2采集地震波數(shù)據(jù)如圖14所示。

圖13 測(cè)線布置圖Fig.13 Schematic diagram of spreads distribution

6.2 數(shù)據(jù)分析

經(jīng)速度分析得動(dòng)校正速度，如圖15所示。從圖15可見動(dòng)校正速度與方位角的關(guān)系可近似為1 個(gè)橢圓，長(zhǎng)軸方向約為335°，由此可見主要裂隙走向約為335°。

圖14 測(cè)線2地震波數(shù)據(jù)Fig.14 Seismic data for spread 2

圖15 不同方位測(cè)線的動(dòng)校正速度Fig.15 NMO velocity for spreads of different azimuths

圖16和圖17所示分別為測(cè)線2 的速度譜和衰減譜。分析圖16和圖17可知：在深度40～65 m處存在1 個(gè)高速度帶，此時(shí)衰減較弱；在65～120 m 時(shí)，雖然速度相對(duì)地表仍然較高，但其衰減系數(shù)較大，地震波衰減較快，由此推測(cè)K74+650—K74+900 范圍內(nèi)，深65～120 m處地下水較豐富。隨后的隧道開挖結(jié)果證明地下水極其豐富，主要裂隙走向約為330°，反演結(jié)果與之相差5°左右。

圖16 測(cè)線2速度譜Fig.16 Velocity spectral of spread 2

圖17 測(cè)線2衰減譜Fig.17 Attenuation spectral for spread 2

7 結(jié)論

1)以HTI 介質(zhì)為基礎(chǔ)，通過坐標(biāo)旋轉(zhuǎn)得到了EDA 裂隙誘導(dǎo)介質(zhì)中地震波的三維動(dòng)校正速度表達(dá)式和走時(shí)方程。以HTI介質(zhì)中的動(dòng)校正速度、慢度向量及走時(shí)方程為基礎(chǔ)，推導(dǎo)出動(dòng)校正速度的射線參數(shù)表達(dá)式。利用最小二乘法對(duì)模擬記錄進(jìn)行非雙曲時(shí)差分析，得到動(dòng)校正速度，然后由動(dòng)校正速度與各向異性參數(shù)及對(duì)稱軸方位角的關(guān)系式進(jìn)行參數(shù)反演。

2)通過對(duì)廣角數(shù)據(jù)進(jìn)行速度分析得到動(dòng)校正速度，采用射線追蹤及多項(xiàng)式擬合法可確定多層傾斜地層界面；由地震波走時(shí)及射線追蹤法，利用多項(xiàng)式擬合法可反演多層傾斜裂隙地層的界面幾何參數(shù)。

3)反演速度與理論計(jì)算速度相對(duì)誤差小于5%；據(jù)各向異性參數(shù)反演可確定EDA 介質(zhì)的裂隙發(fā)育方位。裂隙介質(zhì)參數(shù)反演的相對(duì)誤差均在5%以內(nèi)；加入高斯噪聲后，裂隙方位角和各向異性參數(shù)ε相對(duì)誤差均在10%以內(nèi)，其他相對(duì)誤差均在5%以內(nèi)。

4)速度譜反映的地質(zhì)情況與衰減分析結(jié)果一致，裂隙走向反演結(jié)果相對(duì)誤差為5°。所提出的反演方法可應(yīng)用于隧道裂隙發(fā)育及富水區(qū)預(yù)報(bào)，具有較高的預(yù)報(bào)精度。