0 引言

在各向异性介质中，地震波的传播受多个参数的控制。弹性波波动方程可以模拟地震波的传播过程，但是面临参数多、算法复杂等问题。因为实际地震资料以纵波数据为主^[1-6]，所以许多学者研究了如何准确描述纵波的传播规律。

Alkhalifah^[7-8]从频散关系出发，人为将沿着对称轴方向的qSV波速度设为0，推导了四阶时间—空间域qP波方程以及分裂后的二阶方程组，提出声学近似假设。声学近似假设可以较准确地描述qP波的运动学规律，为各向异性介质纵波勘探提供了理论基础。在此基础上，一些学者开展了关于各向异性介质qP波方程正演方面的研究。Du等^[9]、Zhou等^[10]以及Hestholm^[11]通过简化频散关系，利用不同的降阶处理导出了不同形式的二阶耦合VTI介质qP波方程；吴国忱等^[12-13]在空间—频率域推导了VTI介质qP波波动方程；Duveneck等^[14]从胡克定律及质点方程出发，基于qSV波速度为0的假设，推导了VTI介质二阶耦合qP波波动方程；Xu等^[15-16]基于Alkhalifah的声学近似方程^[7]，推导了各向异性介质的纯qP波波动方程。

虽然Alkhalifah提出的声学近似方法能够较准确地描述qP波，但直接设定沿对称轴方向qSV波波速为0的方法，并不能保证非对称轴方向的qSV波波速也为0，其余方向的残余会在波场中形成菱形的qSV波伪影^[17]。在不改变声学近似原理的基础上，有学者根据菱形伪影产生的原理，在数值求解过程中通过各向异性参数匹配法^[8]、波场滤波法^[18]和辅助波场法^[19]等压制qSV波伪影，但都未能完全消除。而且，声学近似改变了VTI介质原有的弹性介质，只有在Thomsen参数^[20]$ ε \ge \delta $时才能保证声学近似方程计算的稳定性^[17]。

为了解决上述问题，Xu等^[21]将qSV波在各个方向上的波速均设为0，提出一种新声学近似方程。新声学近似方程在理论上可以准确描述纯qP波的传播规律，但形式上却比Alkhalifah的方程更复杂，无法直接用有限差分法求解。Liang等^[22]基于新声学近似的纯qP波频散关系提出一种正演方法，将纯qP波频散关系式通过空间近似^[15]转化为可以直接用有限差分法求解的时间—空间域纯qP波波动方程。但空间近似是利用传播方向矢量代替波数矢量，与Poynting矢量有着类似的缺点，对复杂模型存在模拟精度低的问题^[23]。

本文基于新声学近似的纯qP波频散关系，推导了新的纯qP波波动方程。通过引入中间变量，将频散关系中分数形式的空间—波数域混合算子转化成一个求解中间变量的方程，从而得到时间—空间域纯qP波波动方程组。然后通过有限差分法实现二维VTI介质纯qP波模拟，并与纯qP波等相位面进行对比。数值模拟结果表明新的纯qP波波动方程具有更高的精度。

1 新声学近似

二维VTI介质中qP波和qSV波的精确频散关系分别为^[24]

$ {\omega }_{\mathrm{P}}^{2}=\frac{1}{2}\left\{{V}_{\mathrm{P}0}^{2}\left[(1+2ε ){k}_{x}^{2}+{k}_{z}^{2}\right]+\\{V}_{\mathrm{S}0}^{2}\left({k}_{x}^{2}+{k}_{z}^{2}\right)+\sqrt{D}\right\} $

(1)

$ {\omega }_{\mathrm{S}\mathrm{V}}^{2}=\frac{1}{2}\left\{{V}_{\mathrm{P}0}^{2}\left[(1+2ε ){k}_{x}^{2}+{k}_{z}^{2}\right]+\\{V}_{\mathrm{S}0}^{2}\left({k}_{x}^{2}+{k}_{z}^{2}\right)-\sqrt{D}\right\} $

(2)

式中：$ {\omega }_{\mathrm{P}} $和$ {\omega }_{\mathrm{S}\mathrm{V}} $分别为qP波和qSV波的角频率；$ {k}_{x}=k\mathrm{s}\mathrm{i}\mathrm{n}\theta $和$ {k}_{z}=k\mathrm{c}\mathrm{o}\mathrm{s}\theta $分别是x和z方向的波数，k为波数矢量$ \boldsymbol{k} $的模，θ是$ \boldsymbol{k} $与z轴之间的夹角；$ {V}_{\mathrm{S}0} $和$ {V}_{\mathrm{P}0} $分别为qP波和qSV波沿对称轴方向的相速度；$ ε $和$ \delta $是Thomsen各向异性参数^[20]；而

$ \begin{array}{l}D={\left\{{V}_{\mathrm{P}0}^{2}\left[(1+2ε ){k}_{x}^{2}-{k}_{z}^{2}\right]-{V}_{\mathrm{S}0}^{2}\left({k}_{x}^{2}-{k}_{z}^{2}\right)\right\}}^{2}+\\ 4\left({V}_{\mathrm{P}0}^{2}-{V}_{\mathrm{S}0}^{2}\right)\left[(1+2\delta ){V}_{\mathrm{P}0}^{2}-{V}_{\mathrm{S}0}^{2}\right]{k}_{x}^{2}{k}_{z}^{2}\end{array} $

(3)

Alkhalifah提出的声学近似方法令对称轴方向qSV波速度为0，只能消除对称轴方向的qSV波。为了完全消除各个方向的qSV波，Xu等^[21]将qSV波在各个方向的波速都设为0，即式(2)中$ {\omega }_{\mathrm{S}\mathrm{V}}^{2} $为0，得到$ {V}_{\mathrm{S}0} $关于波数和各向异性参数的函数

$ {V}_{\mathrm{S}0}^{2}=\frac{-2(ε -\delta ){k}_{x}^{2}{k}_{z}^{2}}{(1+2ε ){k}_{x}^{4}+{k}_{z}^{4}+2(1+\delta ){k}_{x}^{2}{k}_{z}^{2}}{V}_{\mathrm{P}0}^{2} $

(4)

然后将式(4)代入式(1)，可得到基于新声学近似的纯qP波频散关系

$ {\omega }_{\mathrm{P}}^{2}={V}_{\mathrm{P}0}^{2}\left[(1+2ε ){k}_{x}^{2}+{k}_{z}^{2}\right]-2(ε -\delta ){V}_{\mathrm{P}0}^{2}G $

(5)

式中

$ G=\frac{{k}_{x}^{2}{k}_{z}^{2}\left({k}_{x}^{2}+{k}_{z}^{2}\right)}{(1+2ε ){k}_{x}^{4}+{k}_{z}^{4}+2(1+\delta ){k}_{x}^{2}{k}_{z}^{2}} $

(6)

为分数形式的空间—波数域混合算子，无法通过有限差分法直接求解。

在推导基于新声学近似的纯qP波频散关系方程时，将各个方向上的qSV波波速均设为0，因此推导出的纯qP波频散关系方程中理论上不包含qSV波。

2 基于空间近似的纯qP波波动方程

为了用有限差分法求解基于新声学近似的纯qP波波动方程，Liang等^[22]构造了波数域标量算子S_k

$ {S}_{k}=\frac{-2\left(ε -\delta \right){k}_{x}^{2}{k}_{z}^{2}}{\left(1+2ε \right){k}_{x}^{4}+{k}_{z}^{4}+2\left(1+\delta \right){k}_{x}^{2}{k}_{z}^{2}} $

(7)

并利用传播方向矢量n代替波数矢量k，即$ \boldsymbol{n}=\boldsymbol{k}/\left|\boldsymbol{k}\right| $，得到空间域标量算子

$ {S}_{n}=\frac{-2\left(ε -\delta \right){n}_{x}^{2}{n}_{z}^{2}}{\left(1+2ε \right){n}_{x}^{4}+{n}_{z}^{4}+2\left(1+\delta \right){n}_{x}^{2}{n}_{z}^{2}} $

(8)

n_x和n_z的计算采用如下近似

$ \left\{\begin{array}{l}{n}_{x}=\frac{\frac{\partial U}{\partial x}}{\sqrt{{\left(\frac{\partial U}{\partial x}\right)}^{2}+{\left(\frac{\partial U}{\partial z}\right)}^{2}}}\\ {n}_{z}=\frac{\frac{\partial U}{\partial z}}{\sqrt{{\left(\frac{\partial U}{\partial x}\right)}^{2}+{\left(\frac{\partial U}{\partial z}\right)}^{2}}}\end{array}\right. $

(9)

式中U为空间域波场。

将算子S_n利用式(9)的近似方法转化到空间域，并与式(5)结合，得到基于空间近似的时间—空间域qP波波动方程

$ \left\{\begin{array}{l}\frac{{\partial }^{2}U}{\partial {t}^{2}}={V}_{\mathrm{P}0}^{2}\left[\left(1+2ε \right)+{S}_{n}\right]\frac{{\partial }^{2}U}{\partial {x}^{2}}+{V}_{\mathrm{P}0}^{2}\left(1+{S}_{n}\right)\frac{{\partial }^{2}U}{\partial {z}^{2}}\\ {S}_{n}=\frac{-2\left(ε -\delta \right){\left(\frac{\partial U}{\partial x}\right)}^{2}{\left(\frac{\partial U}{\partial z}\right)}^{2}}{\left(1+2ε \right){\left(\frac{\partial U}{\partial x}\right)}^{4}+{\left(\frac{\partial U}{\partial z}\right)}^{4}+2\left(1+\delta \right){\left(\frac{\partial U}{\partial x}\right)}^{2}{\left(\frac{\partial U}{\partial z}\right)}^{2}}\end{array}\right. $

(10)

利用空间近似虽然可以将纯qP波频散关系式转化为可以直接用有限差分法求解的时间—空间域纯qP波波动方程, 但由于其波数矢量是用传播方向矢量代替，而传播方向由空间导数计算，因而与利用Poynting矢量计算波传播方向的缺点相似。对于复杂模型，如果有多个波同时到达空间某点，受限于每个空间点上只能表示一个传播方向的原则，导致正演结果精度降低^[23]。

3 新的纯qP波波动方程

为了提高基于新声学近似的纯qP波正演的正演精度，本文提出一种新的纯qP波波动方程。首先将式(5)转换到频率―波数域，为

$ \begin{array}{l}{\omega }^{2}\widetilde{U}={V}_{\mathrm{P}0}^{2}\left[\left(1+2ε \right){k}_{x}^{2}+{k}_{z}^{2}\right]\widetilde{U}-\\ \frac{2\left(ε -\delta \right){V}_{\mathrm{P}0}^{2}\left({k}_{x}^{2}+{k}_{z}^{2}\right){k}_{x}^{2}{k}_{z}^{2}}{\left(1+2ε \right){k}_{x}^{4}+{k}_{z}^{4}+2\left(1+\delta \right){k}_{x}^{2}{k}_{z}^{2}}\widetilde{U}\end{array} $

(11)

式中$ \widetilde{U} $为波数域波场。引入中间变量$ \widetilde{P} $

$ \widetilde{P}=\frac{{k}_{x}^{2}{k}_{z}^{2}}{\left(1+2ε \right){k}_{x}^{4}+{k}_{z}^{4}+2\left(1+\delta \right){k}_{x}^{2}{k}_{z}^{2}}\widetilde{U} $

(12)

则式(11)可表示为

$ \left\{\begin{array}{l}{\omega }^{2}\widetilde{U}={V}_{\mathrm{P}0}^{2}\left[\left(1+2ε \right){k}_{x}^{2}+{k}_{z}^{2}\right]\widetilde{U}-\\ \begin{array}{cc}& \end{array}2{V}_{\mathrm{P}0}^{2}\left(ε -\delta \right)\left({k}_{x}^{2}+{k}_{z}^{2}\right)\widetilde{P}\\ \widetilde{P}\left[\left(1+2ε \right){k}_{x}^{4}+{k}_{z}^{4}+2\left(1+\delta \right){k}_{x}^{2}{k}_{z}^{2}\right]={k}_{x}^{2}{k}_{z}^{2}\widetilde{U}\end{array}\right. $

(13)

转换到时间—空间域，有

$ \left\{\begin{array}{l}\frac{{\partial }^{2}U}{\partial {t}^{2}}={V}_{\mathrm{P}0}^{2}\left[\left(1+2ε \right)\frac{{\partial }^{2}U}{\partial {x}^{2}}+\frac{{\partial }^{2}U}{\partial {z}^{2}}\right]-2{V}_{\mathrm{P}0}^{2}\left(ε -\delta \right)\left(\frac{{\partial }^{2}P}{\partial {x}^{2}}+\frac{{\partial }^{2}P}{\partial {z}^{2}}\right)\\ \frac{{\partial }^{4}U}{\partial {x}^{2}\partial {z}^{2}}=\left(1+2ε \right)\frac{{\partial }^{4}P}{\partial {x}^{4}}+\frac{{\partial }^{4}P}{\partial {z}^{4}}+2\left(1+\delta \right)\frac{{\partial }^{4}P}{\partial {x}^{2}\partial {z}^{2}}\end{array}\right. $

(14)

上式即本文建立的时间—空间域纯qP波波动方程组。与基于空间近似的qP波波动方程(式(10))相比，本文推导过程中未采用任何近似，因此该方程能够精确模拟基于新声学近似的qP波传播规律，理论上具有更高的精度。

采用有限差分法求解式(14)。时间二阶导数采用二阶中心差分格式，空间二阶和四阶导数的有限差分格式为

$ \left\{\begin{array}{l}\frac{{\partial }^{2}f}{\partial {x}^{2}}\approx \frac{1}{{\left(\mathrm{\Delta }x\right)}^{2}}\left({f}_{i+1,j}-2{f}_{i,j}+{f}_{i-1,j}\right)\\ \frac{{\partial }^{2}f}{\partial {z}^{2}}\approx \frac{1}{{\left(\mathrm{\Delta }z\right)}^{2}}\left({f}_{i,j+1}-2{f}_{i,j}+{f}_{i,j-1}\right)\\ \frac{{\partial }^{4}f}{\partial {x}^{4}}\approx \frac{1}{{\left(\mathrm{\Delta }x\right)}^{4}}\left({f}_{i+2,j}-4{f}_{i+1,j}+6{f}_{i,j}-4{f}_{i-1,j}+{f}_{i-2,j}\right)\\ \frac{{\partial }^{4}f}{\partial {z}^{4}}\approx \frac{1}{{\left(\mathrm{\Delta }z\right)}^{4}}\left({f}_{i,j+2}-4{f}_{i,j+1}+6{f}_{i,j}-4{f}_{i,j-1}+{f}_{i,j-2}\right)\\ \frac{{\partial }^{4}f}{\partial {x}^{2}\partial {z}^{2}}\approx \frac{1}{{\left(\mathrm{\Delta }x\right)}^{2}{\left(\mathrm{\Delta }z\right)}^{2}}\left[4{f}_{i,j}-2\left({f}_{i+1,j}+{f}_{i-1,j}+{f}_{i,j+1}+{f}_{i,j-1}\right)\right.+\left.{f}_{i+1,j+1}+{f}_{i+1,j-1}+{f}_{i-1,j-1}+{f}_{i-1,j+1}\right]\end{array}\right. $

(15)

式中：f代表网格点上P或U；Δx和Δz分别是x和z方向的空间采样间隔；i、j分别为两个方向的网格点序号。

式(14)中第二式是求解中间变量P的方程，将其离散后形成一个线性方程组

$ \boldsymbol{A}\boldsymbol{P}=\boldsymbol{U} $

(16)

式中：$ \boldsymbol{A} $为由各向异性参数组成的系数矩阵，大小为(N_x×N_z)×(N_x×N_z)，其中N_x和N_z分别为x和z方向网格点的个数，图 1为其矩阵示意图；$ \boldsymbol{P} $和$ \boldsymbol{U} $分别为每个网格点上P和U组成的向量，大小为N_x×N_z。

通过求解线性方程可以得到中间变量P，将P代入式(14)中第一个公式，通过有限差分法迭代求解U，可得到qP波正演模拟结果。另外，在正演过程中采用完全匹配层(PML)边界条件^[25-26](附录A)。

4 正演数值模拟 4.1 均匀模型试算

为了证明新声学近似能有效消除qSV波伪影并解决$ ε < \delta $情况下计算不稳定的问题，本文将新的纯qP波波动方程分别在$ ε > \delta $(模型A)、$ ε =\delta $(模型B)、$ ε < \delta $(模型C)三种情况下与Alkhalifah的声学近似qP波波动方程及基于空间近似的纯qP波波动方程进行对比。表 1为三个均匀模型的各向异性参数，将qP波的V_P0设为3000 m/s。震源采用主频为20 Hz的Ricker子波，位于模型中心。网格数为301×301，水平和垂直网格间距均设为6 m，时间采样率为0.1 ms，在(200，200)处安置检波器。

图 2为基于Alkhalifah提出的声学近似纯qP波波动方程模拟的三个模型波场快照，可以看出当$ ε > \delta $时波场中会出现明显的qSV波伪影，当$ ε < \delta $时由于计算不稳定不能获得波场图像，只有当$ ε =\delta $时才能有效避免qSV波伪影。

与图 2相比，基于空间近似的纯qP波波动方程模拟的波场快照(图 3)和本文提出的新的纯qP波波动方程模拟的波场快照(图 4)基本消除了qSV波伪影，并且在$ ε < \delta $情况下计算也稳定，能够更精确地描述qP波的传播特征。

为了比较基于空间近似的qP波波动方程与新的纯qP波波动方程的精度，根据Xu等^[21]推出的纯qP波群速度绘出纯qP波的等相位面(图 3、图 4中红线)。由图 3b、图 4b可以看出：当$ ε =\delta $时，两种方法模拟的波场快照与纯qP波的等相位面几乎重合, 并且波场快照中无明显频散。而当$ ε \ne \delta $时，基于空间近似qP波波动方程模拟的波场快照(图 3a、图 3c)与纯qP波等相位面出现了明显的偏差，且波场中出现明显的频散(图 3d、图 3e箭头所示)；而新的纯qP波波动方程模拟的波场快照(图 4a、图 4c)与纯qP波的等相位面几乎一致, 且波场中没有明显的频散(图 4d、图 4e箭头所示)。这是由于纯qP波频散关系中含有分数阶的空间—波数域混合算子G，当$ ε \ne \delta $时，如式(6)~式(8)所示，基于空间近似的纯qP波波动方程对G进行了从波数到空间导数的近似，从而出现了频散现象，同时导致了与纯qP波的等相位面误差。而本文提出的新的纯qP波波动方程则是利用中间变量对G求解，没有采用任何近似处理，故而在波场中并未出现明显的频散现象，且更加接近纯qP波的等相位面。当$ ε =\delta $时，纯qP波频散关系式(5)中G为0，因此两种方法都没有误差。

为了进一步比较三种方法的差异，绘制了均匀介质中检波点处记录的波形，如图 5所示。可以看出：当$ ε \ne \delta $时，基于空间近似方法模拟的波形由于受到频散的影响，波形以及振幅均发生了明显的改变，而本文新的纯qP波波动方程方法模拟的波形以及振幅没有明显改变。综上，本文的纯qP波波动方程在正演模拟中具有更高精度。

4.2 复杂模型试算

应用Sigsbee2a模型和Marmousi-Ⅱ模型对比新的纯qP波波动方程与基于空间近似的纯qP波波动方程对复杂模型的适应能力。将水平和垂直的网格间距均设为10 m，时间采样率为0.1 ms，震源采用主频为15 Hz的Ricker子波，震源位于图 6a的星号处, 并在三角标处设置三个检波点。

图 6a~图 6c为Sigsbee2a模型的速度、各向异性参数。在各向异性参数相同的情况下，基于空间近似波动方程模拟的qP波(图 6d)出现了严重的频散，而本文新的波动方程模拟的qP波(图 6e)几乎没有频散。图 7为两种波动方程在三个检波点处记录的波形对比，可以看出空间近似方法模拟的波形振幅与新的方法相比，其振幅变化更大，原因在于：首先，空间近似只能表示一个传播方向，当多个波重叠时，会降低模拟精度^[23]；其次，空间近似会在一定程度上加重频散现象。而本文新的波动方程在推导过程中没有使用近似，因此对于复杂模型具有更高的模拟精度。

Marmousi-Ⅱ模型的参数如图 8a~图 8c所示。由于受到较为明显的频散干扰，从基于空间近似的纯qP波波动方程模拟的波场快照中(图 8d)无法清晰地观察qP波的传播规律，并且检波点处接收到的波场信息(直达波、反射波)均被频散干扰，使检波点处的记录不够精确(图 9)；而本文提出的新波动方程模拟出的波场快照(图 8e)中未发现明显的频散，可以清楚地观察到qP波在复杂介质中的传播规律，检波点处也能够接收到较准确的波场信息(直达波、反射波)，由此可进一步说明本文提出的纯qP波方程对稍复杂的模型具有更好的适应性，在正演模拟过程中不容易产生频散，正演模拟结果也更加精确。

结合两种模型的各向异性参数的数值分布(图 6b~图 6c、图 8b~图 8c)以及波场快照(图 6d~图 6e、图 8d~图 8e)，可以看出波场快照中出现明显频散的地方，其各向异性参数之间的差值$ \left|ε -\delta \right| $较大，可说明当$ \left|ε -\delta \right| $较大时，本文提出的方法与空间近似的方法相比，具有更高的精度，对复杂的模型具有更好的适应性。

5 结论

本文基于新声学近似的纯qP波频散关系，通过引入中间变量推导了新的纯qP波波动方程，避免了空间近似带来的误差，从理论上提高了正演结果的精度。

均匀VTI介质模型模拟结果完全消除了qSV波伪影，在$ ε < \delta $的情况下也具有稳定性，验证了新的声学近似方法可以准确描述纯qP波的传播规律。与基于空间近似的纯qP波波动方程相比，新的纯qP波波动方程模拟的波场没有明显的频散干扰，与纯qP波的等相位面更接近，模拟精度更高。Sigsbee2a模型和Marmousi-Ⅱ模型模拟结果表明，本文方法与空间近似的方法相比，前者对复杂模型具有更强的适应性，不容易受到频散的干扰，能够准确地描述qP波在复杂介质中的运动规律，具有更高的精度。

附录A PML边界条件

以x方向为例(左、右边界)，通过参考复数伸展坐标系^[27-28]，将新的x方向坐标表示为

$ \tilde{x}\left(x\right)=x-\frac{\mathrm{i}}{\omega }{\int }_{0}^{x}{d}_{x}\left(s\right)\mathrm{d}s $

(A-1)

式中：$ \tilde{x}\left(x\right) $是新的复数空间坐标；$ {d}_{x}\left(s\right) $为边界衰减系数，是一个关于x的实函数，可以表示为

$ {d}_{x}\left(x\right)=\left\{\begin{array}{l}2\mathrm{\pi }\alpha F{\left(\frac{x-{x}_{0}}{{L}_{\mathrm{P}\mathrm{M}\mathrm{L}}}\right)}^{2}\\ 0\end{array}\right. $

(A-2)

其中：α为指数常数，通常取1.79；F是子波主频；L_PML表示PML吸收层的厚度；x₀为边界常数。

对式(A-1)求导可得

$ \frac{\partial \tilde{x}}{\partial x}=1+\frac{{d}_{x}\left(x\right)}{\mathrm{i}\omega } $

(A-3)

记$ \xi =\frac{\partial \tilde{x}}{\partial x} $，$ \frac{1}{\xi }=\frac{\partial x}{\partial \tilde{x}}=\frac{\mathrm{i}\omega }{\mathrm{i}\omega +{d}_{x}\left(x\right)} $，则

$ \frac{\partial }{\partial \tilde{x}}=\frac{\partial }{\partial x}\left(\frac{1}{\xi }\right)=\frac{-\mathrm{i}\omega \cdot {d}_{x}^{\text{'}}\left(x\right)}{{\left[\mathrm{i}\omega +{d}_{x}\left(x\right)\right]}^{2}} $

(A-4)

式中$ {d}_{x}^{\text{'}}\left(x\right) $为$ {d}_{x}\left(x\right) $的一阶导数。

将式(14)中的P和U带入式(A-4)并对其求偏导，可得

$ \left\{\begin{array}{l}\frac{{\partial }^{2}U}{\partial {\tilde{x}}^{2}}=\frac{{\partial }^{2}U}{\partial {x}^{2}}{\left(\frac{1}{\xi }\right)}^{2}+\frac{\partial U}{\partial x}\frac{1}{\xi }\left[\frac{\partial }{\partial x}\left(\frac{1}{\xi }\right)\right]\\ \frac{{\partial }^{2}U}{\partial {\tilde{z}}^{2}}=\frac{{\partial }^{2}U}{\partial {z}^{2}}{\left(\frac{1}{\xi }\right)}^{2}+\frac{\partial U}{\partial z}\frac{1}{\xi }\left[\frac{\partial }{\partial z}\left(\frac{1}{\xi }\right)\right]\\ \frac{{\partial }^{2}P}{\partial {\tilde{x}}^{2}}=\frac{{\partial }^{2}P}{\partial {x}^{2}}{\left(\frac{1}{\xi }\right)}^{2}+\frac{\partial P}{\partial x}\frac{1}{\xi }\left[\frac{\partial }{\partial x}\left(\frac{1}{\xi }\right)\right]\\ \frac{{\partial }^{2}P}{\partial {\tilde{z}}^{2}}=\frac{{\partial }^{2}P}{\partial {z}^{2}}{\left(\frac{1}{\xi }\right)}^{2}+\frac{\partial P}{\partial z}\frac{1}{\xi }\left[\frac{\partial }{\partial z}\left(\frac{1}{\xi }\right)\right]\end{array}\right. $

(A-5)

将波场分裂为$ U={U}_{1}+{U}_{2} $、$ P={P}_{1}+{P}_{2} $，同时引入中间变量Q₁和Q₂，得到频率域PML吸收边界条件

$ \left\{\begin{array}{l}{\left(\mathrm{i}\omega +{d}_{x}\right)}^{2}{U}_{1}=\left(1+2ε \right)\frac{{\partial }^{2}U}{\partial {x}^{2}}+{Q}_{1}\\ {Q}_{1}=-\left(1+2ε \right)\frac{{d}_{x}^{\text{'}}}{\mathrm{i}\omega +{d}_{x}}\frac{\partial U}{\partial x}\\ {\left(\mathrm{i}\omega \right)}^{2}{U}_{2}=\frac{{\partial }^{2}U}{\partial {z}^{2}}\\ {\left(\mathrm{i}\omega +{d}_{x}\right)}^{2}{P}_{1}=-2\left(ε -\delta \right)\frac{{\partial }^{2}P}{\partial {x}^{2}}+{Q}_{2}\\ {Q}_{2}=2\left(ε -\delta \right)\frac{{d}_{x}^{\text{'}}}{\mathrm{i}\omega +{d}_{x}}\frac{\partial P}{\partial x}\\ {\left(\mathrm{i}\omega \right)}^{2}{P}_{2}=-2\left(ε -\delta \right)\frac{{\partial }^{2}P}{\partial {z}^{2}}\end{array}\right. $

(A-6)

则时间域PML吸收边界条件为

$ \left\{\begin{array}{l}{\left({\partial }_{t}+{d}_{x}\right)}^{2}{U}_{1}=\left(1+2ε \right)\frac{{\partial }^{2}U}{\partial {x}^{2}}+{Q}_{1}\\ \left({\partial }_{t}+{d}_{x}\right){Q}_{1}=-\left(1+2ε \right){d}_{x}^{\text{'}}\frac{\partial U}{\partial x}\\ {\left({\partial }_{t}\right)}^{2}{U}_{2}=\frac{{\partial }^{2}U}{\partial {z}^{2}}\\ {\left({\partial }_{t}+{d}_{x}\right)}^{2}{P}_{1}=-2\left(ε -\delta \right)\frac{{\partial }^{2}P}{\partial {x}^{2}}+{Q}_{2}\\ \left({\partial }_{t}+{d}_{x}\right){Q}_{2}=2\left(ε -\delta \right){d}_{x}^{\text{'}}\frac{\partial P}{\partial x}\\ {\left({\partial }_{t}\right)}^{2}{P}_{2}=-2\left(ε -\delta \right)\frac{{\partial }^{2}P}{\partial {z}^{2}}\end{array}\right. $

(A-7)

同理可获得z方向的PML边界条件。