0 引言

弹性波波动方程解耦在弹性波传播理论和实际应用中具有非常重要的意义。在均匀弹性各向同性介质中，对矢量弹性波动方程分别通过求散度和旋度运算，使纵波和横波分离成两个独立传播过程，称为纵波与横波在矢量场中的解耦现象^[1]，即纵波与横波两种波动独立存在，分别满足各自独立传播的波动方程。在大多数各向异性介质中，纵波、横波是耦合在一起传播，描述弹性波传播的方程非常复杂，通常难以解耦，要采用特殊的方法进行处理才能得到各向异性介质中独立传播的纵波、横波波动方程。Alkhalifah^[2-3]提出了基于声学假设的VTI介质和OA介质qP波波动方程；Grechka等^[4]研究了各向异性介质中的横波波动方程；Zhou等^[5]利用二维TTI介质qP波方程进行了模拟和偏移；梁锴等^[6]从频散关系出发研究了三维TTI介质中波动方程的分解；Fowler等^[7]对多个声学假设的qP波方程进行了对比和分析；Duveneck等^[8]和Zhang等^[9]分别提出了TTI介质稳定的qP波方程；Cheng等^[10-11]利用偏振方向投影方法推导了各向异性介质纯qP波和qSV波方程；郭成锋等^[12]利用改进的TTI介质纯P波方程进行了正演模拟和逆时偏移；慕鑫茹等^[13]基于最佳平方逼近方法进行了TTI介质qP波与qSV波的解耦研究；谷一鹏等^[14]利用扩展各向异性线性近似研究了VTI介质qP波和qSV波解耦的相速度近似表征；孙上饶等^[15]基于近似配方法研究了三维TTI介质qP波和qSV波解耦的相速度和群速度近似表征。

横向各向同性(TI)介质是常见的各向异性介质，是薄互层储层或裂缝储层地震资料处理和解释中常用的模型^[16-19]。当TI介质的Thomsen参数满足ε=δ时，被称为椭球各向异性(EA)介质。EA介质的独立弹性常数只有4个。Bakulin等^[20]研究表明含气裂缝介质可以等效为椭球各向异性介质，因此开展椭球各向异性介质的研究具有重要的应用价值。Byun^[21]研究了基于射线理论的椭球各向异性介质的地震参数；Helbig^[22]总结了椭球各向异性介质弹性波波前传播特征和SH波反射透射特征；Thomsen^[23]给出了椭球各向异性介质的Thomsen参数表征；李磊等^[24]讨论了椭球各向异性介质的约束条件；梁锴等^[25]分析了倾斜椭球各向异性介质(TEA)弹性波传播特征，发现在TEA介质中qP波和qSV波的相速度是精确解耦的。

本文在梁锴等^[25]的研究基础上，分析了EA介质弹性波的频散关系和波动方程解耦。从EA介质弹性波波动方程出发，将弹性波频散关系进行因式分解，推导了EA介质完全精确解耦的qP波、qSV波和SH波波动方程，分别描述三种波的独立传播，在此基础上进行了正演模拟并给出了数值示例。

1 EA介质刚度矩阵

各向异性介质中，应力与应变之间线性关系的比例系数由刚度矩阵C确定，而由C表征的弹性波动方程系数的物理意义很不直观。为方便理论研究和实际应用，Thomsen^[23]提出了一套表征TI介质弹性性质的参数，与C的关系可表示为

$ \left\{\begin{array}{l} V_{\mathrm{P} 0}=\sqrt{\frac{C_{33}}{\rho}} \\ V_{\mathrm{S} 0}=\sqrt{\frac{C_{55}}{\rho}} \\ \varepsilon=\frac{C_{11}-C_{33}}{2 C_{33}} \\ \gamma=\frac{C_{66}-C_{44}}{2 C_{44}} \\ \delta=\frac{\left(C_{13}+C_{44}\right)^{2}-\left(C_{33}-C_{44}\right)^{2}}{2 C_{33}\left(C_{33}-C_{44}\right)} \end{array}\right. $

(1)

式中：V_P0、V_S0分别为qP波和qS波沿对称轴方向的相速度；ρ为密度；ε是度量qP波各向异性强度的参数；δ近似表征TI介质对称轴方向qP波相速度对传播角度的二阶导数；γ是度量SH波各向异性强度或横波分裂强度的参数。

当ε=δ或(C₁₃+C₄₄)²=(C₁₁-C₄₄)(C₃₃-C₄₄)时，TI介质退化为EA介质，其Thomsen参数与C的关系为

$ \left\{\begin{array}{l} C_{11}=C_{22}=\rho(1+2 \varepsilon) V_{\mathrm{P} 0}^{2} \\ C_{33}=\rho V_{\mathrm{P0}}^{2} \\ C_{44}=C_{55}=\rho V_{\mathrm{S} 0}^{2} \\ C_{66}=\rho(1+2 \gamma) V_{\mathrm{S0}}^{2} \\ C_{12}=\rho(1+2 \varepsilon) V_{\mathrm{P0}}^{2}-2 \rho(1+2 \gamma) V_{\mathrm{S0}}^{2} \\ C_{13}=C_{23}= \\ \ \ \ \ \rho \sqrt{\left(V_{\mathrm{P} 0}^{2}-V_{\mathrm{S0}}^{2}\right)\left[(1+2 \varepsilon) V_{\mathrm{P} 0}^{2}-V_{\mathrm{S0}}^{2}\right]}-\rho V_{\mathrm{S0}}^{2} \end{array}\right. $

(2)

在本构坐标系下，EA介质的刚度矩阵为

$ \boldsymbol{C}=\left[\begin{array}{cccccc} C_{11} & C_{12} & C_{13} & 0 & 0 & 0 \\ C_{12} & C_{11} & C_{13} & 0 & 0 & 0 \\ C_{13} & C_{13} & C_{33} & 0 & 0 & 0 \\ 0 & 0 & 0 & C_{44} & 0 & 0 \\ 0 & 0 & 0 & 0 & C_{44} & 0 \\ 0 & 0 & 0 & 0 & 0 & C_{66} \end{array}\right] $

(3)

2 EA介质解耦波动方程

波动方程的频散关系在波动方程解耦、单程波方程推导等方面起着非常重要的作用。在各向异性介质中，通过求解Christoffel方程不仅可得到弹性波传播的相速度与群速度，还可求解出EA介质弹性波的精确频散关系方程。

利用EA介质刚度矩阵C，结合弹性动力学的本构方程、运动平衡微分方程和几何方程，可得EA介质弹性波波动方程

$ \left\{\begin{array}{c} C_{11} \frac{\partial^{2} u_{x}}{\partial x^{2}}+C_{66} \frac{\partial^{2} u_{x}}{\partial y^{2}}+C_{44} \frac{\partial^{2} u_{x}}{\partial z^{2}}+\left(C_{12}+C_{66}\right) \times \\ \frac{\partial^{2} u_{y}}{\partial x \partial y}+\left(C_{13}+C_{44}\right) \frac{\partial^{2} u_{z}}{\partial x \partial z}+\rho F_{x}=\rho \frac{\partial^{2} u_{x}}{\partial t^{2}} \\ C_{66} \frac{\partial^{2} u_{y}}{\partial x^{2}}+C_{11} \frac{\partial^{2} u_{y}}{\partial y^{2}}+C_{44} \frac{\partial^{2} u_{y}}{\partial z^{2}}+\left(C_{12}+C_{66}\right) \times \\ \frac{\partial^{2} u_{x}}{\partial x \partial y}+\left(C_{13}+C_{44}\right) \frac{\partial^{2} u_{z}}{\partial y \partial z}+\rho F_{y}=\rho \frac{\partial^{2} u_{y}}{\partial t^{2}} \\ C_{44} \frac{\partial^{2} u_{z}}{\partial x^{2}}+C_{44} \frac{\partial^{2} u_{z}}{\partial y^{2}}+C_{33} \frac{\partial^{2} u_{z}}{\partial z^{2}}+\left(C_{13}+C_{44}\right) \times \\ \frac{\partial^{2} u_{x}}{\partial x \partial z}+\left(C_{13}+C_{44}\right) \frac{\partial^{2} u_{y}}{\partial y \partial z}+\rho F_{z}=\rho \frac{\partial^{2} u_{z}}{\partial t^{2}} \end{array}\right. $

(4)

式中：u=(u_x，u_y，u_z)^T，为位移矢量；F=(F_x，F_y，F_z)^T，为体力项。设弹性波方程(式(4))的平面波解为

$ \boldsymbol{u}=\boldsymbol{p} \exp \left[\mathrm{i}\left(k_{x} x+k_{y} y+k_{z} z-\omega t\right)\right] $

(5)

式中：k_x、k_y、k_z分别为x、y、z方向的波数；p=(p_x，p_y，p_z)^T，为平面波的偏振向量；ω为圆频率。

将平面波解代入EA介质波动方程(式(4))，并忽略体力项的影响，可得EA介质Kelvin-Christoffel方程

$ \left[\begin{array}{ccc} \varGamma_{11}-\rho \omega^{2} & \varGamma_{12} & \varGamma_{13} \\ \varGamma_{12} & \varGamma_{22}-\rho \omega^{2} & \varGamma_{23} \\ \varGamma_{13} & \varGamma_{23} & \varGamma_{33}-\rho \omega^{2} \end{array}\right]\left[\begin{array}{c} p_{x} \\ p_{y} \\ p_{z} \end{array}\right]=0 $

(6)

式中

$ \left\{\begin{array}{l} \varGamma_{11}=\rho\left[(1+2 \varepsilon) V_{\mathrm{P0}}^{2} k_{x}^{2}+(1+2 \gamma) V_{\mathrm{S0}}^{2} k_{y}^{2}+V_{\mathrm{S0}}^{2} k_{z}^{2}\right] \\ \varGamma_{12}=\rho\left[(1+2 \varepsilon) V_{\mathrm{P0}}^{2}-(1+2 \gamma) V_{\mathrm{S0}}^{2}\right] k_{x} k_{y} \\ \varGamma_{13}=\rho \sqrt{\left(V_{\mathrm{P0}}^{2}-V_{\mathrm{S0}}^{2}\right)\left[(1+2 \varepsilon) V_{\mathrm{P0}}^{2}-V_{\mathrm{S0}}^{2}\right]} k_{z} k_{x} \\ \varGamma_{22}=\rho\left[(1+2 \gamma) V_{\mathrm{S0}}^{2} k_{x}^{2}+(1+2 \varepsilon) V_{\mathrm{P0}}^{2} k_{y}^{2}+V_{\mathrm{S0}}^{2} k_{z}^{2}\right] \\ \varGamma_{23}=\rho \sqrt{\left(V_{\mathrm{P0}}^{2}-V_{\mathrm{S0}}^{2}\right)\left[(1+2 \varepsilon) V_{\mathrm{P0}}^{2}-V_{\mathrm{S0}}^{2}\right]} k_{y} k_{z} \\ \varGamma_{33}=\rho\left(V_{\mathrm{S0}}^{2} k_{x}^{2}+V_{\mathrm{S0}}^{2} k_{y}^{2}+V_{\mathrm{P0}}^{2} k_{z}^{2}\right) \end{array}\right. $

(7)

要使式(6)有非零解，必须使其行列式为零，即

$ \begin{aligned} &\left(\rho \omega^{2}\right)^{3}-\left(\varGamma_{11}+\varGamma_{22}+\varGamma_{33}\right)\left(\rho \omega^{2}\right)^{2}+ \\ &\left(\varGamma_{11} \varGamma_{33}+\varGamma_{22} \varGamma_{33}+\varGamma_{22} \varGamma_{11}-\varGamma_{23}^{2}-\varGamma_{13}^{2}-\right. \\ &\left.\varGamma_{12}^{2}\right) \rho \omega^{2}-\left(\varGamma_{11} \varGamma_{22} \varGamma_{33}+2 \varGamma_{12} \varGamma_{13} \varGamma_{23}-\right. \\ &\left.\varGamma_{11} \varGamma_{23}^{2}-\varGamma_{22} \varGamma_{13}^{2}-\varGamma_{33} \varGamma_{12}^{2}\right)=0 \end{aligned} $

(8)

整理式(8)，可得EA介质弹性波精确频散关系

$ \begin{aligned} &{\left[\omega^{2}-(1+2 \varepsilon) V_{\mathrm{P} 0}^{2}\left(k_{x}^{2}+k_{y}^{2}\right)-V_{\mathrm{P0}}^{2} k_{z}^{2}\right] \times} \\ &\left(\omega^{2}-V_{\mathrm{S0}}^{2} k_{x}^{2}-V_{\mathrm{S0}}^{2} k_{y}^{2}-V_{\mathrm{S0}}^{2} k_{z}^{2}\right) \times \\ &{\left[\omega^{2}-(1+2 \gamma) V_{\mathrm{S0}}^{2}\left(k_{x}^{2}+k_{y}^{2}\right)-V_{\mathrm{S0}}^{2} k_{z}^{2}\right]=0} \end{aligned} $

(9)

式(9)可因式分解为三个独立的频散方程，分别是EA介质qP波、qSV波和SH波的解耦频散方程

$ \left\{\begin{array}{l} \omega_{\mathrm{P}}^{2}=(1+2 \varepsilon) V_{\mathrm{P} 0}^{2}\left(k_{x}^{2}+k_{y}^{2}\right)+V_{\mathrm{P0}}^{2} k_{z}^{2} \\ \omega_{\mathrm{SV}}^{2}=V_{\mathrm{S0}}^{2} k_{x}^{2}+V_{\mathrm{S0}}^{2} k_{y}^{2}+V_{\mathrm{S0}}^{2} k_{z}^{2} \\ \omega_{\mathrm{SH}}^{2}=(1+2 \gamma) V_{\mathrm{S0}}^{2}\left(k_{x}^{2}+k_{y}^{2}\right)+V_{\mathrm{S0}}^{2} k_{z}^{2} \end{array}\right. $

(10)

EA介质解耦的波动方程可由EA介质频散关系方程通过Fourier反变换得到，即

$ \frac{\partial^{2} u_{\mathrm{P}}}{\partial t^{2}}=V_{\mathrm{P0}}^{2}\left[(1+2 \varepsilon)\left(\frac{\partial^{2} u_{\mathrm{P}}}{\partial x^{2}}+\frac{\partial^{2} u_{\mathrm{P}}}{\partial y^{2}}\right)+\frac{\partial^{2} u_{\mathrm{P}}}{\partial z^{2}}\right] $

(11)

$ \frac{\partial^{2} u_{\mathrm{SV}}}{\partial t^{2}}=V_{\mathrm{S0}}^{2}\left(\frac{\partial^{2} u_{\mathrm{SV}}}{\partial x^{2}}+\frac{\partial^{2} u_{\mathrm{SV}}}{\partial y^{2}}+\frac{\partial^{2} u_{\mathrm{SV}}}{\partial z^{2}}\right) $

(12)

$ \frac{\partial^{2} u_{\mathrm{SH}}}{\partial t^{2}}=V_{\mathrm{S0}}^{2}\left[(1+2 \gamma)\left(\frac{\partial^{2} u_{\mathrm{SH}}}{\partial x^{2}}+\frac{\partial^{2} u_{\mathrm{SH}}}{\partial y^{2}}\right)+\frac{\partial^{2} u_{\mathrm{SH}}}{\partial z^{2}}\right] $

(13)

式(11)~式(13)由EA介质精确的频散关系解耦导出，没有利用弱各向异性近似或声学假设近似，因此解耦的方程既适合弱各向异性EA介质，也适合强各向异性EA介质。解耦的qP波方程(式(11))只与V_P0、ε有关，与V_S0、γ无关；解耦的qSV波方程(式(12))与各向同性介质的横波方程相同，只与V_S0有关，与V_P0、ε、γ无关，说明在EA介质中qSV波的速度或波前面是各向同性的。解耦的SH波方程(式(13))只与V_S0、γ有关，与V_P0、ε无关。

当EA介质退化为各向同性介质(ε=γ=0)时，式(11)~式(13)也退化为各向同性介质P波和S波的波动方程。

3 数值示例

为了验证EA介质弹性波完全解耦波动方程的正确性和有效性，设计了两个均匀EA介质模型。两个模型的V_P0、V_S0、ρ相同，分别为3000m/s、1732m/s、2000kg/m³；第一个模型ε=γ=0.2，为弱各向异性介质；第二个模型ε=γ=1.0，为强各向异性介质。采用高阶有限差分算法进行正演模拟。

图 1为根据EA介质弹性波群速度公式^[25]计算的250ms时刻理论波前面，可见，在xOz面和yOz面内，qP波和SH波波前为椭圆状，而qSV波为圆形，而在xOy面三种波的波前均为圆形。图 2为弱各向异性EA介质弹性波波动方程的模拟结果，可见，在集中力源激发的弹性波场中包含qP波、qSV波和SH波，qP波和SH波的波前面为椭球状，qSV波波前面为球状，即EA介质中qSV波的速度是各向同性的。图 3为EA介质完全解耦波动方程的模拟结果，由图可见，解耦的qP波场不含有qSV或SH成分，qSV波场不含有qP或SH波成分，SH波波场不含有qP或qSV成分。对于EA介质解耦方程的波场，qP波和SH波的波前面为椭球状，qSV波波前面为球状，说明解耦方程的运动学特征与弹性波方程保持一致，数值模拟结果与理论曲线完全吻合。

图 4为强各向异性EA模型弹性波250ms时刻理论波前面，图 5为强各向异性EA介质弹性波波动方程的模拟结果，图 6为强各向异性EA介质完全解耦波动方程的模拟结果，数值模拟结果与理论曲线完全吻合。图 1~图 6说明这三个完全解耦的波动方程不仅适用于弱各向异性EA介质，也同样适合强各向异性EA介质。同时，EA介质的各向异性越强，qP波和SH波波前面的椭球扁率就越大。

4 结论与认识

一般来说各向异性介质中纵波和横波是耦合传播的，但EA各向异性介质却是一个例外。本文从EA介质精确频散关系方程出发，利用因式分解方法将其分解为三个解耦的频散方程，然后利用傅里叶反变换得到了EA介质完全精确解耦的qP波、qSV波和SH波的波动方程。解耦方程表明，在均匀EA介质中，qP波、qSV波和SH波可以完全解耦，三种波可以独立传播。理论公式和数值模拟均表明，这三个完全解耦方程不仅适用于弱各向异性EA介质，也同样适合强各向异性EA介质。另外，均匀EA介质中qP波和SH波的波前面为椭球状，qSV波的波前面为球状，且qSV波的波前面与各向异性参数无关。

因为倾斜含气裂缝地层可以等效为倾斜椭球各向异性(TEA)介质，本文研究的EA介质是TEA介质对称轴垂直时的特例，因此关于TEA介质弹性波频散关系和解耦波动方程、地震波在倾斜含气裂缝地层中的传播规律，是今后的研究方向。