0 引言

研究表明，地下介质广泛存在地震各向异性，传统的各向同性假设不能准确描述地震波运动学特征^[1-2]。因此，基于各向同性的逆时偏移算法对实际地震资料成像时会产生构造假象、偏移噪声等问题，不能满足当前高精度地震勘探成像要求^[3-5]。为使地震资料准确成像，Crampin^[6]提出一种具有垂直对称轴的横向各向同性(VTI)介质模型近似周期性排列的薄互层和水平排列的平行裂缝等地质体。然而，VTI介质假设不能对陡倾地质体准确成像。如在盐丘侧翼覆盖页岩情况下，对称轴不在垂直方向，若忽略此倾角信息，成像时会产生构造假象和噪声干扰。因此，人们提出了具有倾斜对称轴的横向各向同性(TTI)介质模型弥补VTI介质模型的不足^[7]。

尽管各向异性全弹性波方程可以较准确地描述实际地下介质中的地震波传播特征，但是巨大的计算量、计算内存和波场分离困难等问题限制了其在业界的应用^[8-9]。此外，由于目前地震勘探主要利用纵波反射信息，前人研究了可极大减小计算量与内存需求的拟声波波动方程，以期更好地处理纵波勘探数据^[10-17]。Thomsen^[1]首次提出弱各向异性理论，并用Thomsen各向异性参数表征地下各向异性强度。Tsvankin^[18]推导了由Thomsen参数表示的VTI介质精确相速度公式。为便于纵波成像，Alkhalifah^[19]假设沿着垂直方向的横波速度为零，简化了VTI介质相速度公式，可以降低频散方程的复杂性，且不会明显影响qP波运动学信息的准确性。随后，Alkhalifah^[20]从简化的精确相速度公式出发推导了可以由有限差分法(FDM)求解的四阶qP波波动方程。之后，人们从Alkhalifah的简化频散关系出发，发展了多种不同的耦合拟声波方程^[21-22]。另一种耦合拟声波方程的推导思路是从声学近似的胡克定律出发，得到的拟声波方程可以较准确地描述波场的运动学、动力学特征，方程具有明确的物理意义^[23-25]。由于拟声波近似仅假设沿着垂直对称轴方向的横波速度为零，而其他方向的横波速度并不为零，从而造成模拟的波场中存在伪横波干扰^[26]。为消除伪横波干扰，利用伪横波不能在各向同性介质和椭圆各向异性介质中传播的特点，可在震源周围加一个震源环^[20]或者滤波每一个时刻的波场^[27]。此外，当ε < δ(ε和δ是Thomsen弱各向异性参数)时，耦合拟声波方程模拟的波场存在不稳定问题。

人们通过推导解耦的纯qP波方程解决上述耦合拟声波方程出现的问题。Du等^[11]对VTI介质精确相速度公式进行一阶泰勒展开近似，推导了解耦的TTI介质qP波与qSV波方程，并由伪谱法(PSM)波动方程实现TTI介质逆时偏移。纯qP波方程虽然解决了耦合拟声波方程的一些问题，但是近似相速度公式精度低、PSM效率低、倾角剧烈变化的复杂各向异性介质波场传播不稳定等问题仍制约着其在生产中的应用^[28]。近年来，相继提出了高精度、高效、稳定的纯qP波方程与算法。Chu等^[29]、Zhan等^[30]利用平方根近似替代精确相速度公式或其波数域形式，并分别使用FDM和快速展开法(REM)得到纯qP波方程，大大提高了波场模拟效率。Song等^[31]首先提出使用最佳平方逼近(Optimal Quadratic Approximation，OQA)近似精确相速度公式，随后Mu等^[32]在OQA基础上发展了较高精度的VTI介质四阶近似相速度公式，并证明其精度高于由一阶泰勒近似得到的相速度公式。最近，Li等^[33]发展了一种新的隐式FDM替代传统PSM，实现了TTI介质纯qP波高效、稳定的数值模拟。

本文借鉴前人经验，首先使用OQA近似VTI介质精确相速度公式，分别得到较高精度的近似qP波、qSV波相速度公式。然后从近似相速度公式出发，推导了二维TTI介质解耦qP波、qSV波波动方程。为实现高效、稳定的数值模拟，使用新的有限差分—伪谱混合法(Hybrid Finite-Difference/Pseudo-Spectral Method，H-FD/PSM)求解波动方程。最终，建立了高精度、高效率、稳定的TTI介质正演模拟与逆时偏移成像算法。通过复杂VTI和TTI介质正演模拟与逆时偏移测试，展示了文中方法在复杂各向异性介质地震数据处理中的优势与应用前景。

1 方法原理 1.1 qP波传播方程

二维VTI介质精确相速度公式为^[18]

$ \begin{array}{l} \frac{{v_{\rm{P}}^2\left( \theta \right)}}{{v_{{\rm{P}}0}^2}} = 1 + \varepsilon {\sin ^2}\theta - \frac{f}{2} + \frac{f}{2} \times \\ \;\;\;\;\;\;\;\;\;\sqrt {{{\left( {1 + \frac{{2\varepsilon {{\sin }^2}\theta }}{f}} \right)}^2} - \frac{{8\left( {\varepsilon - \delta } \right){{\sin }^2}\theta {{\cos }^2}\theta }}{f}} \end{array} $

(1)

式中：v_P(θ)为qP波相速度，其中θ为相位角；f=1- $\frac{{v_{{\rm{So}}}^2}}{{v_{{\rm{P}}0}^2}}$，其中v_P0和v_S0分别为沿着对称轴方向的qP波和qSV波速度；ε和δ为Thomsen弱各向异性参数。根据Alkhalifah^[20]提出的拟声学近似(即假设沿对称轴方向的横波速度为零)，可得f≡1，因此式(1)变为

$ \begin{array}{l} \frac{{v_{\rm{P}}^2\left( \theta \right)}}{{v_{{\rm{P}}0}^2}} = \frac{1}{2} + \varepsilon {\sin ^2}\theta + \frac{1}{2} \times \\ \;\;\;\;\;\;\sqrt {{{\left( {1 + 2\varepsilon {{\sin }^2}\theta } \right)}^2} - 8\left( {\varepsilon - \delta } \right){{\sin }^2}\theta {{\cos }^2}\theta } \end{array} $

(2)

Song等^[31]首先使用OQA的思想逼近式(2)的平方根部分，并作以下假设

$ \begin{array}{*{20}{c}} {{f_1}\left( \theta \right) = {{\left[ {{{\left( {1 + 2\varepsilon {{\sin }^2}\theta } \right)}^2} - 8\left( {\varepsilon - \delta } \right){{\sin }^2}\theta {{\cos }^2}\theta } \right]}^{\frac{1}{2}}}}\\ {\theta \in \left( {0,2{\rm{ \mathsf{ π} }}} \right)} \end{array} $

(3)

OQA的思想(详见附录A)就是找到一组系数λ₁, λ₂, …, λ_n，求

$ \mathop {\min }\limits_{\theta \in \left( {0,2{\rm{ \mathsf{ π} }}} \right)} \left\{ {\int_a^b {{{\left[ {{f_1}\left( \theta \right) - \sum\limits_{i = 1}^n {{\lambda _i}{\varphi _i}\left( \theta \right)} } \right]}^2}} {\rm{d}}\theta } \right\} $

(4)

Song等^[31]使用函数最高次幂为2的二阶函数系{sin²θ, 1}近似f₁(θ)，但是仅在椭圆各向异性时精度很高，当Thomsen各向异性参数相差较大时误差较大。这是因为二阶函数系本质上假设ε=δ，从而忽略了式(3)中的非椭圆各向异性信息。为了补偿被忽略的非椭圆各向异性信息，本文使用四阶函数系{sin⁴θ, sin²θ, 1}和六阶函数系{sin⁶θ, sin⁴θ, sin²θ, 1}近似f₁(θ)，以期得到更高的精度。将VTI介质的垂直对称轴旋转变为TTI介质，假设新的对称轴与垂直方向的夹角为ϕ，分别得到TTI介质二阶、四阶、六阶近似qP波相速度公式

$ \begin{array}{l} \frac{{2v_{\rm{P}}^2\left( {\theta ,\phi } \right)}}{{v_{{\rm{P}}0}^2}} = 1 + 2\varepsilon {\sin ^2}\left( {\theta - \phi } \right) + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\lambda _1}{\sin ^2}\left( {\theta - \phi } \right) + {\lambda _2} \end{array} $

(5)

$ \begin{array}{l} \frac{{2v_{\rm{P}}^2\left( {\theta ,\phi } \right)}}{{v_{{\rm{P}}0}^2}} = 1 + 2\varepsilon {\sin ^2}\left( {\theta - \phi } \right) + \\ \;\;\;\;\;\;\;\;\;{\lambda _1}{\sin ^4}\left( {\theta - \phi } \right) + {\lambda _2}{\sin ^2}\left( {\theta - \phi } \right) + {\lambda _3} \end{array} $

(6)

$ \begin{array}{l} \frac{{2v_{\rm{P}}^2\left( {\theta ,\phi } \right)}}{{v_{{\rm{P}}0}^2}} = 1 + 2\varepsilon {\sin ^2}\left( {\theta - \phi } \right) + {\lambda _1}{\sin ^6}\left( {\theta - \phi } \right) + \\ \;\;\;\;\;\;\;\;\;{\lambda _2}{\sin ^4}\left( {\theta - \phi } \right) + {\lambda _3}{\sin ^2}\left( {\theta - \phi } \right) + {\lambda _4} \end{array} $

(7)

众所周知，高阶近似的精度更高，但将消耗更多的计算时间和计算内存，特别是TTI介质。图 1为TTI介质qP波相速度曲线。由图可见：当Thomsen参数相差较大时，二阶近似与精确值存在较大差异(图 1b~图 1d)；四阶、六阶近似与精确值差异都非常小，且六阶近似的精度略高于四阶近似，但是六阶近似的计算量和计算内存更大、推导过程复杂。综合考虑计算精度与计算效率，本文采用四阶近似qP波相速度公式。

为了证明对于具有不同方向对称轴的TTI介质四阶近似qP波相速度公式仍具有较高精度，利用式(6)计算具有不同方向对称轴的TTI介质qP波相速度曲线(图 2)。可见，具有不同方向对称轴的4种TTI介质的qP波相速度四阶近似与精确值只有微小的差异，可以较准确地模拟qP波传播。

已知TTI介质平面波具有如下关系

$ \left\{ {\begin{array}{*{20}{l}} {\sin \theta = \frac{{v\left( {\theta ,\phi } \right){k_x}}}{\omega }}\\ {\cos \theta = \frac{{v\left( {\theta ,\phi } \right){k_z}}}{\omega }} \end{array}} \right. $

(8)

式中：k_x和k_z分别为x、z方向的波数成分；ω为角频率。将式(8)代入式(6)并乘以傅里叶域波场U_P(k_x, k_z, ω)，同时仅在频率域应用傅里叶反变换，可得二维TTI介质qP波波动方程

$ \begin{array}{l} \frac{{{\partial ^2}{U_{\rm{P}}}\left( {{k_x},{k_z},t} \right)}}{{\partial {t^2}}} = - \frac{{v_{{\rm{P}}0}^2}}{2}\left[ {k_x^2 + k_z^2 + \left( {{\lambda _3} + } \right.} \right.\\ \;\;\;\;\left. {{\lambda _1}{{\cos }^4}\phi + 2\varepsilon {{\cos }^2}\phi + {\lambda _2}{{\cos }^2}\phi } \right)\frac{{k_x^4}}{{k_x^2 + k_z^2}} + \\ \;\;\;\;\left( {{\lambda _3} + {\lambda _1}{{\sin }^4}\phi + 2\varepsilon {{\sin }^2}\phi + {\lambda _2}{{\sin }^2}\phi } \right)\frac{{k_z^4}}{{k_x^2 + k_z^2}} + \\ \;\;\;\;\left( {2{\lambda _3} + {\lambda _2} + 2\varepsilon + 1.5{\lambda _1}{{\sin }^2}2\phi } \right)\frac{{k_x^2k_z^2}}{{k_x^2 + k_z^2}} + \\ \;\;\;\;\left( { - 2\varepsilon \sin 2\phi - {\lambda _2}\sin 2\phi - 2{\lambda _1}\sin 2\phi {{\cos }^2}\phi } \right)\frac{{k_x^3{k_z}}}{{k_x^2 + k_z^2}} + \\ \;\;\;\;\left. {\left( { - 2\varepsilon \sin 2\phi - {\lambda _2}\sin 2\phi - 2{\lambda _1}\sin 2\phi {{\sin }^2}\phi } \right)\frac{{{k_x}k_z^3}}{{k_x^2 + k_z^2}}} \right] \times \\ \;\;\;\;{U_{\rm{P}}}\left( {{k_x},{k_z},t} \right) \end{array} $

(9)

1.2 qSV波传播方程

同样地，从VTI介质qSV波精确相速度公式出发^[18]，推导TTI介质qSV波波动方程。qSV波精确相速度公式为

$ \begin{array}{l} \frac{{v_{\rm{S}}^2\left( \theta \right)}}{{v_{{\rm{P}}0}^2}} = 1 + \varepsilon {\sin ^2}\theta - \frac{f}{2} - \frac{f}{2} \times \\ \;\;\;\;\;\;\sqrt {{{\left( {1 + \frac{{2\varepsilon {{\sin }^2}\theta }}{f}} \right)}^2} - \frac{{8\left( {\varepsilon - \delta } \right){{\sin }^2}\theta {{\cos }^2}\theta }}{f}} \end{array} $

(10)

此处，首先定义

$ \begin{array}{l} {f_2}\left( \theta \right) = \left[ {{{\left( {1 + \frac{{2\varepsilon {{\sin }^2}\theta }}{f}} \right)}^2} - } \right.\\ \;\;\;\;\;\;\;\;\;{\left. {\frac{{8\left( {\varepsilon - \delta } \right){{\sin }^2}\theta {{\cos }^2}\theta }}{f}} \right]^{\frac{1}{2}}}\quad \theta \in \left( {0,2{\rm{ \mathsf{ π} }}} \right) \end{array} $

(11)

分别用二阶函数系{sin²θ, 1}、四阶函数系{sin⁴θ, sin²θ, 1}、六阶函数系{sin⁶θ, sin⁴θ, sin²θ, 1}近似f₂(θ)，分别得到TTI介质二阶、四阶、六阶近似qSV波相速度公式

$ \begin{array}{*{20}{c}} {\frac{{v_{\rm{S}}^2\left( {\theta ,\phi } \right)}}{{v_{{\rm{P0}}}^2}} = 1 + \varepsilon {{\sin }^2}\left( {\theta - \phi } \right) - \frac{f}{2} - }\\ {\frac{f}{2}\left[ {{\lambda _1}{{\sin }^2}\left( {\theta - \phi } \right) + {\lambda _2}} \right]} \end{array} $

(12)

$ \begin{array}{l} \frac{{v_{\rm{S}}^2\left( {\theta ,\phi } \right)}}{{v_{{\rm{P0}}}^2}} = 1 + \varepsilon {\sin ^2}\left( {\theta - \phi } \right) - \frac{f}{2} - \\ \;\;\;\;\;\;\;\frac{f}{2}\left[ {{\lambda _1}{{\sin }^2}\left( {\theta - \phi } \right) + {\lambda _2}{{\sin }^2}\left( {\theta - \phi } \right) + {\lambda _3}} \right] \end{array} $

(13)

$ \begin{array}{l} \frac{{v_{\rm{S}}^2\left( {\theta ,\phi } \right)}}{{v_{{\rm{P0}}}^2}} = 1 + \varepsilon {\sin ^2}\left( {\theta - \phi } \right) - \frac{f}{2} - \\ \;\;\;\;\;\;\;\;\frac{f}{2}\left[ {{\lambda _1}{{\sin }^6}\left( {\theta - \phi } \right) + {\lambda _2}{{\sin }^4}\left( {\theta - \phi } \right) + } \right.\\ \;\;\;\;\;\;\;\;\left. {{\lambda _3}{{\sin }^2}\left( {\theta - \phi } \right) + {\lambda _4}} \right] \end{array} $

(14)

图 3为TTI介质qSV波相速度曲线。由图可见：当Thomsen参数相差较大时，二阶近似与精确值存在较大差异(图 3b~图 3d)；四阶、六阶近似与精确值差异都非常小，且六阶近似的精度略高于四阶近似，但是六阶近似的计算量和计算内存更大、推导过程复杂。综合考虑计算精度与计算效率，本文采用四阶近似qSV波相速度公式。

为了证明对于具有不同方向对称轴的TTI介质四阶近似qSV波相速度公式仍具有较高精度，利用式(13)得到具有不同方向对称轴的TTI介质qSV波相速度曲线(图 4)。可见，具有不同方向对称轴的4种TTI介质的qSV波相速度四阶近似精度略低于qP波，但也可以较准确地模拟qSV波传播。

将式(8)代入式(13)并乘以傅里叶域波场U_S(k_x, k_z, ω)，同时仅在频率域应用反傅里叶变换，可得二维TTI介质qSV波波动方程

$ \begin{array}{l} \frac{{{\partial ^2}{U_{\rm{S}}}\left( {{k_x},{k_z},t} \right)}}{{\partial {t^2}}} = v_{{\rm{P}}0}^2\left[ {\left( {1 - \frac{f}{2} + \varepsilon {{\cos }^2}\phi - } \right.} \right.\\ \;\;\;\;\;\;\left. {\frac{f}{2}{\lambda _1}{{\cos }^4}\phi - \frac{f}{2}{\lambda _2}{{\cos }^2}\phi - \frac{f}{2}{\lambda _3}} \right)\frac{{k_x^4}}{{k_x^2 + k_z^2}} + \\ \;\;\;\;\;\;\left( {1 - \frac{f}{2} + \varepsilon {{\sin }^2}\phi - \frac{f}{2}{\lambda _1}{{\sin }^4}\phi - \frac{f}{2}{\lambda _2}{{\sin }^2}\phi - } \right.\\ \;\;\;\;\;\;\left. {\frac{f}{2}{\lambda _3}} \right)\frac{{k_z^4}}{{k_x^2 + k_z^2}} + \left( {2 - f + \varepsilon - } \right.\\ \;\;\;\;\;\;\left. {3f{\lambda _1}{{\sin }^2}\phi {{\cos }^2}\phi - \frac{f}{2}{\lambda _2} - f{\lambda _3}} \right)\frac{{k_x^2k_z^2}}{{k_x^2 + k_z^2}} + \\ \;\;\;\;\;\;\left( { - 2\varepsilon \sin \phi \cos \phi + 2f{\lambda _1}\sin \phi {{\cos }^3}\phi + } \right.\\ \;\;\;\;\;\;\left. {f{\lambda _2}\sin \phi \cos \phi } \right)\frac{{k_x^3{k_z}}}{{k_x^2 + k_z^2}} + ( - 2\varepsilon \sin \phi \cos \phi + \\ \;\;\;\;\;\;\left. {\left. {2f{\lambda _1}{{\sin }^3}\phi \cos \phi + f{\lambda _2}\sin \phi \cos \phi } \right)\frac{{{k_x}k_z^3}}{{k_x^2 + k_z^2}}} \right] \times \\ \;\;\;\;\;\;{U_{\rm{S}}}\left( {{k_x},{k_z},t} \right) \end{array} $

(15)

图 5为图 2和图 4的极坐标表示形式，反映了在弱各向异性假设某一时刻均匀介质(速度场为常数)波场快照的波前面形状^[11]，qP波(图 5上)和qSV波(图 5)的四阶近似相速度和真实相速度几乎没有差异，间接说明本文推导的qP波和qSV波波动方程可以准确地描述TTI介质中波的传播特征。

1.3 H-FD/PSM

式(9)和式(15)可方便地采用PSM求解，但是PSM计算量大且易出现数值不稳定。为提高复杂TTI介质数值模拟稳定性和计算效率，本文借鉴前人经验，提出利用H-FD/PSM求解波动方程^[32-33]。此时，式(9)可重写为

$ \begin{array}{l} \frac{{{\partial ^2}{U_{\rm{P}}}\left( {x,z,t} \right)}}{{\partial {t^2}}} = \frac{{v_{{\rm{P}}0}^2}}{2}\left[ {\left( {1 + {\lambda _3} + {\lambda _1}{{\cos }^4}\phi + } \right.} \right.\\ \;\;\;\left. {2\varepsilon {{\cos }^2}\phi + {\lambda _2}{{\cos }^2}\phi } \right)\frac{{{\partial ^4}{U_q}\left( {x,z,t} \right)}}{{\partial {x^4}}} + \\ \;\;\;\left( {1 + {\lambda _3} + {\lambda _1}{{\sin }^4}\phi + 2\varepsilon {{\sin }^2}\phi + } \right.\\ \;\;\;\left. {{\lambda _2}{{\sin }^2}\phi } \right)\frac{{{\partial ^4}{U_q}\left( {x,z,t} \right)}}{{\partial {z^4}}} + \\ \;\;\;\left( {2 + 2{\lambda _3} + {\lambda _2} + 2\varepsilon + 1.5{\lambda _1}{{\sin }^2}2\phi } \right)\frac{{{\partial ^2}{U_q}\left( {x,z,t} \right)}}{{\partial {x^2}\partial {z^2}}} + \\ \;\;\;\left( { - 2\varepsilon \sin 2\phi - {\lambda _2}\sin 2\phi - 2{\lambda _1}\sin 2\phi {{\cos }^2}\phi } \right) \times \\ \;\;\;\;\frac{{{\partial ^2}{U_q}\left( {x,z,t} \right)}}{{\partial {x^3}\partial z}} + \left( { - 2\varepsilon \sin 2\phi - {\lambda _2}\sin 2\phi - } \right.\\ \;\;\;\;\;\left. {\left. {2{\lambda _1}\sin 2\phi {{\sin }^2}\phi } \right)\frac{{{\partial ^2}{U_q}\left( {x,z,t} \right)}}{{\partial x\partial {z^3}}}} \right] \end{array} $

(16)

其中

$ {\nabla ^2}{U_q}\left( {x,z,t} \right) = {U_{\rm{P}}}\left( {x,z,t} \right) $

(17)

对于泊松方程(式(17))，可用PSM、快速泊松变换、LU分解、多重网格法等方法求解^[34-35]。为便于计算，本文用PSM简单实现高精度求解，那么

$ {U_q}\left( {x,z,t} \right) = {\rm{FF}}{{\rm{T}}^{ - 1}}\left\{ {\frac{{{\rm{FFT}}\left[ {{U_{\rm{P}}}\left( {x,z,t} \right)} \right]}}{{ - k_x^2 - k_z^2}}} \right\} $

对于式(16)，可用时间二阶、空间高阶FDM求解。因此，可用H-FD/PSM求解以下TTI介质qSV波控制方程

$ \begin{array}{l} \frac{{{\partial ^2}{U_{\rm{S}}}\left( {x,z,t} \right)}}{{\partial {t^2}}} = v_{{\rm{P}}0}^2\left[ {\left( {1 - \frac{f}{2} + \varepsilon {{\cos }^2}\phi - \frac{f}{2}{\lambda _1}{{\cos }^4}\phi - } \right.} \right.\\ \;\;\;\;\;\;\;\;\left. {\frac{f}{2}{\lambda _2}{{\cos }^2}\phi - \frac{f}{2}{\lambda _3}} \right)\frac{{{\partial ^4}{U_k}\left( {x,z,t} \right)}}{{\partial {x^4}}} + \\ \;\;\;\;\;\;\;\;\left( {1 - \frac{f}{2} + \varepsilon {{\sin }^2}\phi - \frac{f}{2}{\lambda _1}{{\sin }^4}\phi - } \right.\\ \;\;\;\;\;\;\;\;\left. {\frac{f}{2}{\lambda _2}{{\sin }^2}\phi - \frac{f}{2}{\lambda _3}} \right)\frac{{{\partial ^4}{U_k}\left( {x,z,t} \right)}}{{\partial {x^4}}} + \\ \;\;\;\;\;\;\;\;\left( {2 - f + \varepsilon - 3f{\lambda _1}{{\sin }^2}\phi {{\cos }^2}\phi - } \right.\\ \;\;\;\;\;\;\;\;\left. {\frac{f}{2}{\lambda _2} - f{\lambda _3}} \right)\frac{{{\partial ^4}{U_k}\left( {x,z,t} \right)}}{{\partial {x^2}\partial {z^2}}} + \left( { - 2\varepsilon \sin \phi \cos \phi } \right. + \\ \;\;\;\;\;\;\;\;\left. {2f{\lambda _1}\sin \phi {{\cos }^3}\phi + f{\lambda _2}\sin \phi \cos \phi } \right)\frac{{{\partial ^4}{U_k}\left( {x,z,t} \right)}}{{\partial {x^3}\partial z}} + \\ \;\;\;\;\;\;\;\;\left( { - 2\varepsilon \sin \phi \cos \phi + 2f{\lambda _1}{{\sin }^3}\phi \cos \phi + } \right.\\ \;\;\;\;\;\;\;\;\left. {\left. {f{\lambda _2}\sin \phi \cos \phi } \right)\frac{{{\partial ^4}{U_k}\left( {x,z,t} \right)}}{{\partial x\partial {z^3}}}} \right] \end{array} $

(18)

$ {\nabla ^2}{U_k}\left( {x,z,t} \right) = {U_{\rm{s}}}\left( {x,z,t} \right) $

(19)

对式(18)用FDM求解，式(19)用PSM求解。

图 6为二维TTI Wedge模型，存在模型参数剧烈变化界面，用该模型验证本文提出的H-FD/PSM的稳定性和计算效率。根据上述模拟参数，分别采用H-FD/PSM和PSM进行数值模拟。图 7为二维TTI Wedge模型波场快照。由图可见：①H-FD/PSM波场快照(图 7a~图 7c)较稳定，没有出现奇异值，其中t＝4.0s时刻(图 7c)的残余波场是FDM处理突变界面时产生的数值噪声，这是由于波传播一段时间后大数值的有效界面反射被吸收边界完全吸收，残留了小数值噪声，但这种数值噪声不影响偏移成像；②PSM波场快照(图 7d)不稳定。表 1为两种数值模拟方法的耗时。可见，PSM耗时约为H-FD/PSM耗时的2.5倍，且由于PSM需要使用多次正、反傅里叶变换，也会增大耗时。

1.4 逆时偏移成像

逆时偏移基于双程波动方程，具有成像精度高、不受倾角限制、适用于横向速度变化等特点^[37-39]。叠前逆时偏移过程包括三个步骤：

(1) 地震波场沿时间方向正向延拓至最大时刻，为了节省TTI介质波场传播计算内存，本文采用波场重构的策略，只保留最后两个时刻的全波场和每个时刻的边界区域波场^[40-41]。

(2) 检波点波场在时间方向从最大时刻结合炮记录逆时延拓至零时刻。

(3) 为补偿地震波传播的深层能量损失，提高成像分辨率，本文采用震源归一化互相关成像条件对两次正、反向传播得到的同一时刻波场成像，成像条件为

$ I\left( {x,z} \right) = \frac{{\int_0^{{T_{\max }}} {S\left( {x,z,t} \right)R\left( {x,z,t} \right){\rm{d}}t} }}{{\sum\limits_0^{{T_{\max }}} {{S^2}\left( {x,z,t} \right)} }} $

(20)

式中：S(x, z, t)为震源波场在时间方向的正向延拓；R(x, z, t)为检波点波场在时间方向的反向延拓。

图 8和图 9分别为qP波和qSV波脉冲响应。由图可见：①图 8无伪横波干扰；图 8和图 9在ϕ=90°时仍具有清晰的能量，表明逆时偏移没有构造倾角限制，能对陡倾构造准确成像。②考虑到各向异性影响，随着地层倾角变化，各向异性对称轴也随之改变(图 8、图 9箭头所示)，导致qP波和qSV波的波前面不同于各向同性介质(图 5)。

2 模型测试

为了验证本文提出的TTI介质逆时偏移算法的有效性，对Overthrust TTI和BP2007 TTI模型进行逆时偏移试算。

2.1 Overthrust TTI模型

模型参数如图 10所示，使用H-FD/PSM得到偏移炮记录，然后分别利用不同方法对炮记录进行逆时偏移。图 11为OverthrustTTI模型逆时偏移结果。由图可见：各向同性(图 11a)和VTI(图 11b)逆时偏移导致成像剖面中TTI异常体下伏地层的同相轴连续性变差(图中椭圆框区域)，且异常体内部倾斜界面模糊、界面归位不准(图中黑色箭头所示)，并引入了噪声；TTI逆时偏移可对模型精确成像(图 11c)。

2.2 BP2007 TTI模型

鉴于原始BP2007 TTI模型较大，本文仅截取其中一部分进行成像测试，速度模型与各向异性参数模型如图 12所示。图 13为基于OQA的纯qP波叠前逆时偏移成像结果。由图可见：各向同性逆时偏移(图 13a)、VTI逆时偏移(图 13b)造成盐丘边界成像位置出现较大偏差(图中箭头与实线框所示)；TTI逆时偏移(图 13c)的同相轴振幅横向均衡性变好。这主要是由于各向同性逆时偏移忽略了各向异性影响、VTI逆时偏移忽略了TTI介质对称轴倾角的变化，从而不同程度地造成目标体成像位置偏差，甚至由焦散问题造成振幅不均衡。此外，由于高速盐丘的遮挡和单边接收观测系统导致盐丘左、右侧成像能量不均衡，致使左侧深部部分同相轴未能成像。因此，本文提出的TTI介质逆时偏移算法可对复杂模型准确成像。

3 结束语

本文利用最佳平方逼近近似VTI介质精确相速度公式，得到高精度的qP波和qSV波相速度公式；从相速度公式出发，推导了新的TTI介质解耦qP波和qSV波波动方程。与常规耦合拟声波方程相比，本文推导的方程无伪横波干扰，可避免伪横波干扰产生的成像噪声。同时采用有限差分—伪谱混合法求解波动方程，与传统伪谱法求解方式相比，具有更高计算效率与稳定性。数值模拟实验表明，本文方法可实现复杂TTI介质偏移成像。

附录A 最佳平方逼近原理

假设f(x)是区间x∈[a, b]上的连续函数，φ_j(x)(j=0, 1, …, n)是x∈[a, b]上的一个线性无关函数系，且在x∈[a, b]上都是连续的，确定广义多项式

$ \varphi \left( x \right) = \sum\limits_{j = 0}^n {{\lambda _j}{\varphi _j}\left( x \right)} $

(A-1)

的系数λ₀, λ₁, …, λ_n，使

$ \min \left\{ {\int_a^b {{{\left[ {f\left( x \right) - \varphi \left( x \right)} \right]}^2}{\rm{d}}x} } \right\} $

(A-2)

这样得到的函数φ(x)是f(x)在x∈[a, b]上的最佳平方逼近。

最佳平方逼近的系数解可表示为

$ \left( {\begin{array}{*{20}{c}} {{\lambda _0}}\\ {{\lambda _1}}\\ {{\lambda _2}}\\ \vdots \\ {{\lambda _n}} \end{array}} \right) = {\left( {\begin{array}{*{20}{c}} {s_0^0}&{s_1^0}&{s_2^0}& \cdots &{s_n^0}\\ {s_0^1}&{s_1^1}&{s_2^1}& \cdots &{s_n^1}\\ {s_0^2}&{s_1^2}&{s_2^2}& \cdots &{s_n^2}\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ {s_0^n}&{s_1^n}&{s_2^n}& \cdots &{s_n^n} \end{array}} \right)^{ - 1}}\left( {\begin{array}{*{20}{c}} {{m_0}}\\ {{m_1}}\\ {{m_1}}\\ \cdots \\ {{m_n}} \end{array}} \right) $

(A-3)

其中

$ s_j^i = \int_a^b {\left[ {{\varphi _i}(x) \times {\varphi _j}(x)} \right]{\rm{d}}x} \;\;\;i = 0,1, \cdots ,n $

(A-4)

$ {m_j} = \int_a^b {\left[ {{\varphi _j}(x) \times f(x)} \right]{\rm{d}}x} $

(A-5)