裂隙介質(zhì)中地震波傳播特征近似

何現(xiàn)啟，彭凌星，朱自強，魯光銀

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

賀同江等[1]用迭積微分算子法實現(xiàn)了黏彈性介質(zhì)的地震波場正演模擬,對波傳播特征進行了分析。宮猛等[2]從各向同性介質(zhì)中波場數(shù)值模擬的褶積微分算子法出發(fā),推導(dǎo)出了各向異性雙相介質(zhì)中波場傳播數(shù)值計算的褶積新算法。周銘等[3]應(yīng)用背景噪聲層析成像方法對新疆地區(qū)地殼S波速度結(jié)構(gòu)以及徑向各向異性進行了研究,討論了各向異性產(chǎn)生的機理。Vigen等[4]討論奇異波相鄰方向上近似qS波群速度和極化的問題。Behura等[5]以各向同性黏彈性介質(zhì)為背景介質(zhì),利用攝動法推導(dǎo)出黏彈性 VTI 介質(zhì)中P 波、PS 波的反射系數(shù)近似公式。鄧繼新等[6]利用比奧特(Biot)理論研究了介觀尺度孔隙流體流動介質(zhì)中地震波的傳播規(guī)律。司薌等[7]推導(dǎo)了傾斜裂隙等效TTI介質(zhì)PP波反射系數(shù)公式,建立了反射系數(shù)與裂隙參數(shù)之間的關(guān)系。Tariq[8]利用擾動法求解了橢圓各向異性介質(zhì)中的聲波方程,并用Marmousi模型進行了數(shù)值驗證。張煜等[9]研究了地震勘探頻段內(nèi)針對含流體孔隙介質(zhì)邊界條件的面波的傳播特性。李雨生等[10]運用交錯網(wǎng)格高階有限差分法求解三維單斜各向異性介質(zhì)彈性波方程來模擬彈性波的傳播過程。Alexey等[11]用多參數(shù)擾動法求解了正交各向異性介質(zhì)中的波動方程,大大縮短了求解時間。何現(xiàn)啟等[12]用擾動法推導(dǎo)了極端弱各向異性介質(zhì)中 qP、qS波群速度的一階表達式。何現(xiàn)啟等[13]研究黏彈性各向異性介質(zhì)中地震波的傳播特征,推導(dǎo)出黏彈性垂直定向裂隙的各向異性介質(zhì)中均勻、非均勻波的精確相速度、慢度和群速度公式。盧方超等[14]通過超聲實驗和核磁共振(NMR)實驗,研究了煤孔、裂隙結(jié)構(gòu)應(yīng)力變形各向異性特征。何現(xiàn)啟等[15]通過數(shù)值模型計算,研究了黏彈性HTI(horizontal transverse isotropy)介質(zhì)中地震波的傳播特征。梁鍇等[16]推導(dǎo)了橢球各向異性介質(zhì)是橫向各向同性介質(zhì)中彈性波速度和偏振表達式。趙小龍等[17]用擾動法推導(dǎo)了正交各向異性介質(zhì)巖石力學近似方程。Song等[18]用擾動法多弱各向異性介質(zhì)中P、S波的反射、投射系數(shù)進行了一、二階近似計算。

盡管對各向異性介質(zhì)中地震波傳播有不少報道,但對裂隙介質(zhì)中地震波傳播進行系統(tǒng)研究的甚少,且由于裂隙介質(zhì)中地震波的精確傳播特征表達較復(fù)雜,非線性的隱性表達式不利于參數(shù)反演,因此有必要對其進行線性化近似。

本文為了簡化Christoffel方程,建立了波向量坐標系(x′，o′，z′),其中z′軸方向與波的傳播方向一致,y′為水平方向,其單位向量為e(1)、e(2)、e(3)表示。當波的傳播方向給定時,各向異性介質(zhì)中Christoffel方程的獨立彈性系數(shù)為六個,而在各向同性介質(zhì)中Christoffel矩陣為對角陣,因此通過對各向同性介質(zhì)的擾動便可對弱各向異性介質(zhì)進行攝動分析。

1 裂隙介質(zhì)中相速度精確表達式

由裂隙介質(zhì)的彈性矩陣CEDA,可得到裂隙介質(zhì)Christoffel方程，求解方程可得裂隙介質(zhì)中的qP波、qSV波和SH波的相速度為[15,19-20]

(1)

式(1)中:D=[(c11-c55)F-(c33-c55)E2]2+4(c13+c55)2FE2，E=(-cosφsinθsinφ0+sinφcosφ0)，F(xiàn)=(cosφsinθcosφ0+sinφsinφ0)2+cos2θcos2φ。其中，φ為HTI介質(zhì)對稱軸方位角，φ0為裂隙介質(zhì)對稱軸方位角，θ為介質(zhì)對稱軸傾角，c為裂隙介質(zhì)彈性系數(shù)。

同理，可得到偏振向量及群速度的精確表達式。

2 裂隙介質(zhì)中相速度一階近似

建立野外觀測坐標系xoz，其基本單位向量為(x,y,z)，波的傳播方向n=(nx,ny,nz)=(sinθcosφ, sinθsinφ,cosθ)，其中θ為極化角，φ為方位角。通過擾動法可推導(dǎo)出弱各向異性條件下，極端各向異性介質(zhì)中地震波傳播特征的近似表達式，進一步將HTI介質(zhì)的彈性系數(shù)矩陣代入相關(guān)公式即可推導(dǎo)出HTI介質(zhì)中地震波傳播特征近似公式。通過將HTI介質(zhì)的彈性系數(shù)矩陣進行旋轉(zhuǎn)，可得到裂隙介質(zhì)(extensive dilatancy anisotropy,EDA)的彈性系數(shù)矩陣：

(2)

將式(2)代入極端各向異性介質(zhì)計算公式，即可推導(dǎo)出裂隙弱各向異性介質(zhì)中qP、qS波的一階相下速度[13,15]:

(3)

3 裂隙介質(zhì)中一階偏振向量

將HTI介質(zhì)中的模型坐標轉(zhuǎn)換到測量(勘探)坐標,并用裂隙弱各向異性參數(shù)表示,可得qP、qS波的一階偏振單位向量[13,15-16,20]。

(4)

式(4)中：偏振角

4 裂隙介質(zhì)中一階群速度

將裂隙介質(zhì)的彈性系數(shù)矩陣代入極端各向異性介質(zhì)的計算公式可得[13,16,19]:

(5)

式(5)中：w′1±=F1(θ,φ)±[G1(θ,φ)cos2α+H1(θ,φ)sin2α],w′2±=F2(θ,φ)±[G2(θ,φ)cos2α+H2(θ,φ)sin2α]；并且有:

P(θ,φ)=sin2θ[sin2θ(εxcos4φ+εysin4φ+

δzcos2φsin2φ)+cos2θ(δysin2φ+

δxcos2φ)]+ sin2φsin4θ(ε16cos2φ+

ε26sin2φ)+χzsin2φcos2θsin2θ

(6)

εysin2φ(sin2φcos2θ+cos2φ) -

cos2θ(δysin2φ+δxcos2φ)-

ε26sin2φsin2θ(2cos2θcos2φ+cos2φ)-

(7)

εysin2φ(sin2φcos2θ-cos2φ)-

cos2θ(δysin2φ+δxcos2φ)+

ε26sin2φsin2θ(2cos2θcos2φ-cos2φ)-

2cos2θcos2φ)

(8)

[ε26sin2φ(4cos2φ-1)+ε16cos2φ(1-

(9)

εysin4φ)+cos2θ(δysin2φ+δxcos2φ)]+

2sin2θsin2φsin2θ(ε16cos2φ+ε26sin2φ)

(10)

cos2θ(δy-δx)+δzcos2φsin2θ]+

ε16sin3θcos2φ(2cos2φ-1)+

ε26sin3θsin2φ(2cos2φ+1)]+

(11)

εysin2φ(2sin2φsin2θ-1)+cos2θ×

(12)

(13)

sin2θ)[ε16(cos2φcos2φsinθ-

(14)

cos2φ)(sin2φ+cos2φcos2θ)+

(15)

ε16sinθ(cos2φ-sin22φ)-

ε26sinθcos4φ]-

(16)

cos2φsin2θ)+2εysin2φ(2cos2φ-

sin2φsin2θ)]+δy(2sin2φsin2θ-1)+

ε16sin2θsin2φ(cos2φsin2θ-

cos2φ)ε26sin2θsin2φ(cos2φ-

sin2θsin2φ)]+

(17)

5 算法驗證

下面主要通過模型計算來驗證式(3)的準確性,模型主要參數(shù)見表1,相關(guān)計算結(jié)果見圖1～圖9。

表1 模型參數(shù)

由表1中模型1參數(shù)可計算HTI彈性系數(shù),通過坐標旋轉(zhuǎn)得到裂隙介質(zhì)的彈性系數(shù)矩陣如下:

(18)

由表1中模型2參數(shù)可計算HTI彈性系數(shù),通過坐標旋轉(zhuǎn)可得到裂隙介質(zhì)的彈性系數(shù)矩陣如下:

(19)

圖1 模型1裂隙介質(zhì)精確相速度平視圖Fig.1 Flat view of exact phase velocity of fracture media for model 1

圖2 模型1裂隙介質(zhì)的精確相速度俯視圖Fig.2 Top view of exact phase velocity of fracture media for model 1

圖3 模型1 裂隙介質(zhì)的近似相速度平視圖Fig.3 Flat view of approximation velocity of model 1

圖4 模型1 裂隙介質(zhì)的近似相速度俯視圖Fig.4 Top view of approximation velocity of model 1

圖5 模型1近似相速度相對誤差Fig.5 Relative error of approximation phase velocity for model1

圖6 模型2精確相速度Fig.6 Exact phase velocity of model 2

圖7 模型2近似相速度Fig.7 Approximation velocity of model 2

圖8 模型2 近似相速度的相對誤差Fig.8 Relative error of approximation phase velocity for model 2

圖9 xoy平面內(nèi)P、SH、SV的精確與近似相速度比較Fig.9 Comparation of exact phase velocity and approximation velocity for P、SH、SV in xoy

以上精確相速度通過解christoffel方程得到,近似解由式(3)得到,分析圖1～圖9可知,當各向異性系數(shù)為0.1時,S波(包括SH、SV波)近似相速度相對誤差小于10%,而P波的近似相速度相對誤差小于2%；當各向異性系數(shù)為0.3時,P波相速度相對誤差在10%以內(nèi),而S波(包括SH、SV波)近似相速度相對誤差最大到20%。當各向異性系數(shù)為0.5時,而P波近似相速度相對誤差小于20%,而S波(包括SH、SV波)近似相速度相對誤差最大達50%,由此可見,當各向異性系數(shù)達0.5時,對S波(包括SH、SV波)而言,式(3)計算精度已滿足不了要求,而對P波仍然適用。

6 結(jié)論

為了縮短數(shù)值計算時間,便于進行速度分析及參數(shù)反演,需將地震波的隱性、非線性精確表達式進行近似研究,轉(zhuǎn)換為線性表達式。在擾動法求解弱各向異性條件下qP、qS的相速度、群速度的近似表達式的基礎(chǔ)上,推導(dǎo)了裂隙介質(zhì)中地震波傳播特征的近似表達式。

(1)依據(jù)介質(zhì)模型坐標與觀測坐標,各種介質(zhì)彈性系數(shù)矩陣關(guān)系式推導(dǎo)出裂隙介質(zhì)中地震波相速度、偏振向量及群速度的近似表達。

(2)通過數(shù)值計算,研究了近似計算誤誤差,結(jié)果表明,當各向異性系數(shù)達在0.3以內(nèi)時,對一般裂隙各向異性介質(zhì)都適用,且具有較高的精確,在各向異性系數(shù)為0.5時,P波仍然具有較高的精度,但S波的誤差增大到50%。

(3)因為擾動解是在弱各向異性條件下推導(dǎo)的,相關(guān)近似公式在弱各向異性條件下,精度較高,當各向異性系數(shù)大于0.3時,對S波應(yīng)謹慎使用。

近似計算公式雖然是在弱各向異性條件下推導(dǎo)得出的,但在各向異性系數(shù)達0.3時,對一般各向異性介質(zhì)都適用,且具有較高的精確,為地震波速度分析及參數(shù)反演奠定了基礎(chǔ),對強各向異性介質(zhì)中地震波傳播特征近似研究是下一需要開展的工作。