各向異性介質中地震波參數反演及其在隧道不良地質預報中的應用

何現啟，彭凌星，魯光銀，朱自強，高峰

(1.湖南省交通規劃勘察設計院有限公司，湖南長沙，410200；2.中南大學地球科學與信息物理學院，湖南長沙，410083)

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

1 EDA 介質中傾斜界面的動校正速度及時距方程

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

由HTI介質精確動校正速度表達式，經坐標旋轉可得到垂直非均勻EDA 介質中傾斜界面的精確動校正速度[23-25]：

式中：α為共中心點排列的方位角；φ為EDA介質的對稱軸的方位角；?為對稱軸傾角；vnmo為動校正速度。

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

由式(2)可知動校正速度vnmo(φ,?)主要由垂直各向異性參數ε(V)與δ(V)的差值和水平界面的動校正速度vnmo(0)決定。由HTI介質中的動校正公式，經坐標旋轉同樣可得沿界面走向的動校正速度vnmo((π/2)+φ,?)為

1.2 對稱面內P波動校正速度與射線參數的關系

由HTI 介質中射線動校正速度推導[25]，EDA 介質中沿界面傾向的動校正速度可用射線參數表示為

式中：q'為一階導數；q''為二階導數；p為慢度向量的水平分量；q為慢度向量的垂直分量，

式中：θ=α-φ。

沿界面走向的動校正速度可用射線參數表示為

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

2 EDA 介質中三維動校正速度與射線參數的關系

2.1 均勻EDA介質中用射線參數表示的動校正速度精確表達式

由動校正速度和慢度向量及2層介質中的時距方程可得[24-25]

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

2.2 弱各向異性EDA介質中對稱面內射線參數與動校正速度的關系

EDA 介質可由HTI 介質繞z旋轉角度φ得到，對垂直方向的參數沒有影響，因此，計算EDA 介質的垂直慢度仍可采用HTI介質中的計算公式。

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

其中：vP0為縱波速度；

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

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

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

在傾斜EDA 介質中，對動校正速度除考慮其測線方位角外，還要考慮傾角的影響。綜合分析HTI介質對稱面內的動校正速度公式以及其與坐標軸的相互關系，采用剝層法由相似分析得到[22]：

也可用下式計算：

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

4 EDA介質中的P波參數反演算法

對于1個給定的垂直時間t0，通過對水平速度vhor和動校正速度vnmo進行二維相似分析可以得到1 個具有最大相似值的模型，相似系數可由下式計算：

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

通過對多個排列進行速度分析即可確定動校正橢圓，由動校正速度與各向異性參數的關系式即可計算出各向異性參數[25]。

所采用的目標函數如下：

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

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

5 模型算例分析

5.1 地層模型

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

表1 模型1 彈性參數Table1 Elastic parameter of Model 1

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

用黏彈性實射線追蹤法實現黏彈性EDA 中地震波的數值模擬，模擬地震波形見圖2。

5.2 速度分析

采用式(16)對速度進行分析，結果見圖3。從圖3可以看出：0°測線的疊加速度顯示反演結果與模型值基本一致，相對誤差在5%以內，動校正結果顯示的界面層數及界面產狀與實際情況相符；同時，0°，10°，20°和30°測線的疊加速度和動校正速度結果均與實際情況一致，說明速度及動校正公式計算結果準確度高。

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

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

5.3 反演結果與誤差分析

5.3.1 無噪聲反演及誤差分析

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

其他相關參數反演結果見圖5～7，其中，實線表示參數的精確值，黑點表示反演值。從圖5～7 可見：迭代300次時，計算精度較高收斂較快，誤差分析結果顯示除方位角誤差稍大外，其余參數相對誤差都非常小，均在5%以內。

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

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

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

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

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

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

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

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

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

圖7 密度反演結果Fig.7 Inversion results of density

從圖9可見：加入高斯噪聲后，ε反演結果較離散，速度反演值在實線附近，相對誤差在8%以內。

從圖10可見：第1層介質的方位角反演值在實線附近，相對誤差在5%以內；第2 層介質反演相對誤差在8%以內；第3 層介質反演結果較離散，反演誤差較大。

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

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

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

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

圖11 密度反演結果Fig.11 Inversion results of density

圖12 各向異性參數δ反演結果Fig.12 Inversion results of anisotropy parameter δ

6 工程實例

為了驗證該方法的可行性，將其應用于裂隙發育區及富水區預報。通過現場測試，分析地層裂隙走向、節理密度與速度分布的關系，研究地震波速度及衰減與地下水賦存狀態的關系。

6.1 現場測試概況

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

圖13 測線布置圖Fig.13 Schematic diagram of spreads distribution

6.2 數據分析

經速度分析得動校正速度，如圖15所示。從圖15可見動校正速度與方位角的關系可近似為1 個橢圓，長軸方向約為335°，由此可見主要裂隙走向約為335°。

圖14 測線2地震波數據Fig.14 Seismic data for spread 2

圖15 不同方位測線的動校正速度Fig.15 NMO velocity for spreads of different azimuths

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

圖16 測線2速度譜Fig.16 Velocity spectral of spread 2

圖17 測線2衰減譜Fig.17 Attenuation spectral for spread 2

7 結論

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

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

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

4)速度譜反映的地質情況與衰減分析結果一致，裂隙走向反演結果相對誤差為5°。所提出的反演方法可應用于隧道裂隙發育及富水區預報，具有較高的預報精度。