裂隙介质中地震波传播特征近似

何现启，彭凌星，朱自强，鲁光银

(1.湖南省交通规划勘察设计院有限公司,长沙 410200；2.中南大学地球科学与信息物理学院,长沙 410083)

贺同江等[1]用迭积微分算子法实现了黏弹性介质的地震波场正演模拟,对波传播特征进行了分析。宫猛等[2]从各向同性介质中波场数值模拟的褶积微分算子法出发,推导出了各向异性双相介质中波场传播数值计算的褶积新算法。周铭等[3]应用背景噪声层析成像方法对新疆地区地壳S波速度结构以及径向各向异性进行了研究,讨论了各向异性产生的机理。Vigen等[4]讨论奇异波相邻方向上近似qS波群速度和极化的问题。Behura等[5]以各向同性黏弹性介质为背景介质,利用摄动法推导出黏弹性 VTI 介质中P 波、PS 波的反射系数近似公式。邓继新等[6]利用比奥特(Biot)理论研究了介观尺度孔隙流体流动介质中地震波的传播规律。司芗等[7]推导了倾斜裂隙等效TTI介质PP波反射系数公式,建立了反射系数与裂隙参数之间的关系。Tariq[8]利用扰动法求解了椭圆各向异性介质中的声波方程,并用Marmousi模型进行了数值验证。张煜等[9]研究了地震勘探频段内针对含流体孔隙介质边界条件的面波的传播特性。李雨生等[10]运用交错网格高阶有限差分法求解三维单斜各向异性介质弹性波方程来模拟弹性波的传播过程。Alexey等[11]用多参数扰动法求解了正交各向异性介质中的波动方程,大大缩短了求解时间。何现启等[12]用扰动法推导了极端弱各向异性介质中 qP、qS波群速度的一阶表达式。何现启等[13]研究黏弹性各向异性介质中地震波的传播特征,推导出黏弹性垂直定向裂隙的各向异性介质中均匀、非均匀波的精确相速度、慢度和群速度公式。卢方超等[14]通过超声实验和核磁共振(NMR)实验,研究了煤孔、裂隙结构应力变形各向异性特征。何现启等[15]通过数值模型计算,研究了黏弹性HTI(horizontal transverse isotropy)介质中地震波的传播特征。梁锴等[16]推导了椭球各向异性介质是横向各向同性介质中弹性波速度和偏振表达式。赵小龙等[17]用扰动法推导了正交各向异性介质岩石力学近似方程。Song等[18]用扰动法多弱各向异性介质中P、S波的反射、投射系数进行了一、二阶近似计算。

尽管对各向异性介质中地震波传播有不少报道,但对裂隙介质中地震波传播进行系统研究的甚少,且由于裂隙介质中地震波的精确传播特征表达较复杂,非线性的隐性表达式不利于参数反演,因此有必要对其进行线性化近似。

本文为了简化Christoffel方程,建立了波向量坐标系(x′，o′，z′),其中z′轴方向与波的传播方向一致,y′为水平方向,其单位向量为e(1)、e(2)、e(3)表示。当波的传播方向给定时,各向异性介质中Christoffel方程的独立弹性系数为六个,而在各向同性介质中Christoffel矩阵为对角阵,因此通过对各向同性介质的扰动便可对弱各向异性介质进行摄动分析。

1 裂隙介质中相速度精确表达式

由裂隙介质的弹性矩阵CEDA,可得到裂隙介质Christoffel方程，求解方程可得裂隙介质中的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=(cosφsinθcosφ0+sinφsinφ0)2+cos2θcos2φ。其中，φ为HTI介质对称轴方位角，φ0为裂隙介质对称轴方位角，θ为介质对称轴倾角，c为裂隙介质弹性系数。

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

2 裂隙介质中相速度一阶近似

建立野外观测坐标系xoz，其基本单位向量为(x,y,z)，波的传播方向n=(nx,ny,nz)=(sinθcosφ, sinθsinφ,cosθ)，其中θ为极化角，φ为方位角。通过扰动法可推导出弱各向异性条件下，极端各向异性介质中地震波传播特征的近似表达式，进一步将HTI介质的弹性系数矩阵代入相关公式即可推导出HTI介质中地震波传播特征近似公式。通过将HTI介质的弹性系数矩阵进行旋转，可得到裂隙介质(extensive dilatancy anisotropy,EDA)的弹性系数矩阵：

(2)

将式(2)代入极端各向异性介质计算公式，即可推导出裂隙弱各向异性介质中qP、qS波的一阶相下速度[13,15]:

(3)

3 裂隙介质中一阶偏振向量

将HTI介质中的模型坐标转换到测量(勘探)坐标,并用裂隙弱各向异性参数表示,可得qP、qS波的一阶偏振单位向量[13,15-16,20]。

(4)

式(4)中：偏振角

4 裂隙介质中一阶群速度

将裂隙介质的弹性系数矩阵代入极端各向异性介质的计算公式可得[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)的准确性,模型主要参数见表1,相关计算结果见图1～图9。

表1 模型参数

由表1中模型1参数可计算HTI弹性系数,通过坐标旋转得到裂隙介质的弹性系数矩阵如下:

(18)

由表1中模型2参数可计算HTI弹性系数,通过坐标旋转可得到裂隙介质的弹性系数矩阵如下:

(19)

图1 模型1裂隙介质精确相速度平视图Fig.1 Flat view of exact phase velocity of fracture media for model 1

图2 模型1裂隙介质的精确相速度俯视图Fig.2 Top view of exact phase velocity of fracture media for model 1

图3 模型1 裂隙介质的近似相速度平视图Fig.3 Flat view of approximation velocity of model 1

图4 模型1 裂隙介质的近似相速度俯视图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平面内P、SH、SV的精确与近似相速度比较Fig.9 Comparation of exact phase velocity and approximation velocity for P、SH、SV in xoy

以上精确相速度通过解christoffel方程得到,近似解由式(3)得到,分析图1～图9可知,当各向异性系数为0.1时,S波(包括SH、SV波)近似相速度相对误差小于10%,而P波的近似相速度相对误差小于2%；当各向异性系数为0.3时,P波相速度相对误差在10%以内,而S波(包括SH、SV波)近似相速度相对误差最大到20%。当各向异性系数为0.5时,而P波近似相速度相对误差小于20%,而S波(包括SH、SV波)近似相速度相对误差最大达50%,由此可见,当各向异性系数达0.5时,对S波(包括SH、SV波)而言,式(3)计算精度已满足不了要求,而对P波仍然适用。

6 结论

为了缩短数值计算时间,便于进行速度分析及参数反演,需将地震波的隐性、非线性精确表达式进行近似研究,转换为线性表达式。在扰动法求解弱各向异性条件下qP、qS的相速度、群速度的近似表达式的基础上,推导了裂隙介质中地震波传播特征的近似表达式。

(1)依据介质模型坐标与观测坐标,各种介质弹性系数矩阵关系式推导出裂隙介质中地震波相速度、偏振向量及群速度的近似表达。

(2)通过数值计算,研究了近似计算误误差,结果表明,当各向异性系数达在0.3以内时,对一般裂隙各向异性介质都适用,且具有较高的精确,在各向异性系数为0.5时,P波仍然具有较高的精度,但S波的误差增大到50%。

(3)因为扰动解是在弱各向异性条件下推导的,相关近似公式在弱各向异性条件下,精度较高,当各向异性系数大于0.3时,对S波应谨慎使用。

近似计算公式虽然是在弱各向异性条件下推导得出的,但在各向异性系数达0.3时,对一般各向异性介质都适用,且具有较高的精确,为地震波速度分析及参数反演奠定了基础,对强各向异性介质中地震波传播特征近似研究是下一需要开展的工作。