多波多分量勘探技术是地层岩性识别、陡倾界面成像的有效手段，能够提供更详细的地下介质信息，在油气和地热等能源探测领域具有重要意义。相比于传统的声波逆时偏移，弹性波逆时偏移考虑了地震波在地下传播时的多波耦合特征，可以较为充分地利用弹性波多分量地震数据中丰富的纵横波场信息，实现对复杂地质结构区块的高精度成像^[1]。传统弹性波逆时偏移由于成像时未进行纵横波场分离，导致耦合波场间有串扰的产生，严重影响成像质量^[2]。为了去除串扰伪影和提高成像质量，国内外学者做了许多研究^[3-5]，提出采用波场分离的方法可以进行有效压制。常用的波场分离方法有基于 Helmholtz 分解的波场分离以及基于解耦方程的波场分离。Helmholtz分解自提出以来，就被广泛地应用于弹性波场分解中^[6-8]。Sun等 ^[9]在进行弹性波场外推时，使用散度及旋度算子分离出地表附近的纵波场和横波场，再分别进行PP和转换PS的两次声波逆时偏移。Dellinger等 ^[10]基于Helmholtz理论提出了波数域的类散度算子和类旋度算子，并利用Christoffel方程实现了各向异性均匀弹性介质中的波场分解。Yan等 ^[11]将分离算子从波数域进一步转换到空间域，并制定了一个滤波器分离具有垂直对称轴的横向各向同性介质中的波场。但是以上基于Helmholtz 分解得到的纵波和横波波场通常被理解为标量势与矢量势，得到的振幅以速度的倒数衰减，相位推迟π/2^[1]，不能准确地完成波场分离，需要进行额外的保幅处理和相位校正。为了探寻一种波场分离方法以更贴近于弹性波的传播规律，基于解耦方程进行波场分离的研究也取得了许多进展。马德堂等 ^[12]根据 Helmholtz 分解理论推导出纵波和横波分离的二阶位移方程，通过伪谱法实现了在波场外推过程中完成纵波和横波的分离，这为后续研究提供了一种新思路。李振春等^[13]和Zhang等^[14]进一步推导出一阶速度-应力弹性波场的分离方程，并利用交错网格高阶有限差分算法验证了二维方程的可行性；Xiao 等^[15]引入辅助 P 波和 S 波应力场进行波场分解，实现了VSP局部弹性逆时偏移，使得基于解耦方程的波场分离更加具有矢量传播的物理意义； Wang等 ^[16]证明 Zhang 和 Xiao的公式是形式不一样的等价方程，并给出了新纵波和横波分离的弹性波方程。基于解耦方程得到的是矢量纵波和矢量横波，若仍然采取互相关成像条件进行成像，则得到的反射率模型与传统意义上的不一样^[17]，会给解释工作人员带来困难。为此适用于矢量波场成像的内积成像条件被进一步发展^[18-19]，这种成像条件能较为准确地提取矢量波场各位置的成像值。然而逆时偏移是一种构造成像方法，对地热体的准确模拟和成像还可以考虑温度和物质形态等因素变化^[20-22]，进而发展出更适用的方法。笔者首先介绍波场分离的基本公式和所应用的成像条件；其次将上述提到的公式和成像条件对层状模型进行波场模拟，通过波场快照和单道分析说明其可行性；最后将矢量波场分离公式应用到以侵入花岗型和近代火山型为代表的深部干热岩体模型进行逆时偏移，以检测本方法的优越性、适用性和准确性。

1 方法原理

本文中弹性波逆时偏移可分为3步：①构建震源波场并保存各时刻波场值；②构建检波波场并保存各时刻波场值；③使用震源归一化内积成像条件进行成像。

1.1 二维波动方程

波动方程反映的是受到激励时地下介质每时刻的运动情况，它能表示地震波的传播过程。二维情况下，通常假定地震波没有y分量，即$\frac{\partial }{\partial y}=0$，所以可给出各向同性介质纵横波耦合传播的一阶速度-应力方程 ^[23]为

式中，v_x和v_z分别为地震波的横向速度和纵向速度；τ_xx、τ_zz和τ_xz分别为x方向正应力、z方向正应力和xz方向的剪切应力；λ和μ为拉梅系数；ρ为介质密度。

1.2 弹性波场分离

传统的波场分离方法是基于 Helmholtz 分解对总波场取散度得到纵波场，对波场取旋度得到横波场。此类波场分离方法对原始矢量波场施加了散度和旋度算子，相当于对其求空间导数，导致分离得到的纵波和横波分量产生90°相移，振幅分别产生 $\frac{1}{v_{\mathrm{P}}}$、$\frac{1}{v_{\mathrm{S}}}$ 的改变。Wang等^[16]通过引入一个辅助纵波应力，进一步发展了一种波场矢量分解公式，能较好地解决这种现象，有效地保存了弹性波场的振幅和相位信息。

式中，τ_P为引入的纵波应力；v_xP和v_zP分别为纵波在x和z方向上的速度分量；v_xS和v_zS分别为横波在x和z方向上的速度分量。

为了实现震源波场延拓与检波波场延拓，采用时间二阶、空间八阶的规则网格有限差分格式对式（1）和（2）进行离散。在进行波场模拟时设定完全匹配层（PML）吸收边界，减弱由于人工截断边界而导致的边界反射现象，更有效地压制边界反射和提高目标模拟区域的成像精度。

1.3 内积成像条件

考虑到基于解耦波动方程得出的纵波和横波是矢量形式，若仍采用传统的互相关成像条件来成像，不能反映出反射系数的物理意义，内积成像条件是解决此问题的有效手段^[19]。复杂模型逆时偏移剖面的深部成像通常存在照明不足等问题，使用震源归一化技术可以有效地补偿地震波的能量损失，有利于增强模型深部的能量。为此给出一种震源归一化内积成像条件为

式中，I（x）为各空间位置成像值；S为震源波场；R为检波波场；下标P和S分别表示纵波和横波；符号“·”为内积算子。

例如，I_PS（x）表示震源纵波与检波横波在成像点x处的内积成像值，内积计算格式由S_P（x，t）·R_S（x，t）=|S_P（x，t）||R_S（x，t）|cos〈S_P（x，t），R_S（x，t）〉得出，其中S_P（x，t）表示震源纵波场在t时刻、x位置处的波场值，R_S（x，t）表示检波横波场在t时刻、x位置处的波场值。

通常情况下地震波的主要能量入射角小于40°，然而内积成像条件发生极性反转的条件是入射角大于45°，因此采用内积成像条件不会出现极性反转现象^[17]。针对使用双程波动方程进行波场模拟而导致成像结果存在低频噪音的问题，采用Laplacian滤波对所有逆时偏移剖面进行处理，有效地压制了此现象。

2 数值模拟

为了验证本文提到的波场分离方法的有效性以及对深部高温地热岩体构造成像的适用性，通过对层状模型、侵入花岗型和近代火山型干热岩模型进行测试。层状模型用于验证基于解耦方程进行波场分离的准确性，侵入花岗型和近代火山型用于检验基于解耦方程对深部地热岩体等复杂地质结构成像的优越性和适用性。本文中均使用时间二阶、空间八阶的规则网格有限差分格式完成波场延拓过程，采用集中力震源完成激发（震源增益为1 000）和完全匹配层边界条件进行吸收。横波速度根据泊松体的纵横波比求得，即$V_{\mathrm{S}}=\frac{V_{\mathrm{P}}}{1.732}$，介质密度根据Gardner公式求取，即ρ=0.31V^0.25_P。

2.1 层状模型

为了更为直观地认识基于解耦方程的波场分离方法，采用一个三层模型（图1）测试其有效性和准确性，该模型假定纵波速度分别为3000、4000、3000 m/s，各层厚度分别为600、1800和600 m。模型的网格点数为401×301，网格间距为10 m，时间采样率为1 ms。震源位于模型正中央（201，151），震源函数为25 Hz的雷克子波。

从0.45 s时刻记录到的波场快照中可以看到，基于解耦方程的波场分离效果较好，可以将纯纵波、纯横波和转换波从耦合波场中实现分离，各波型关于震源分别上下、左右对称，符合地震波传播的基本规律。

为了更清楚展示波场分离的振幅和相位特征，在所有波场快照的1.5 km处（图2（a）中白线示意位置）抽取单道。图3（a）~（c）分别为耦合波场、纵波场和横波场水平分量单道信息，图3（d）~（f）分别为耦合波场、纵波场和横波场垂向分量单道信息。通过对比可以看出，分离得到的纵波和横波振幅没有衰减为耦合波场的$\frac{1}{V_{\mathrm{P}}}，\frac{1}{V_{\mathrm{S}}}$，相位也没有延迟，均基本上保留了与耦合波场一样的振幅和相位特性。

由此可见，基于解耦方程的弹性波场可以较为准确地分离出纵波和横波，在传播过程中可以保持其振幅和相位基本不变，无需再额外进行振幅补偿和相位校正。相比于传统基于Helmholtz分解对波场取散度和旋度的方法，该方法更简易、涉及的步骤更少。

2.2 侵入花岗型和近代火山型干热岩

侵入花岗型和近代火山型作为干热岩的两种常见类型，具有较大的研究价值，发展一种高精度成像方法，对地热能的探测和开发具有重要意义。已经对基于解耦方程的弹性波场进行了分析，通过数值模拟论证了其优点，然后使用该方程分别对面积为10 km×7 km的侵入花岗模型和近代火山模型目标区域进行逆时偏移成像研究，所采用的偏移速度模型均由真实速度模型进行5点平滑后所得。

2.2.1 侵入花岗型

侵入花岗型干热岩是由以岩浆岩深层侵入为主的强烈构造活动导致。图4为侵入花岗模型的纵波速度分布。该模型含2条正断层，地层间多为不整合接触关系，较多位置出现尖灭现象，构造比较复杂。模型中最小速度为2000 m/s，最大速度为5500 m/s。花岗岩基岩位于5 km深度以下，有3条主支呈无规则形状向上侵入5套地层可至3 km。

对模型采样1001×701个，横纵空间采样间隔均为10 m，时间采样率为0.3 ms，时间采样21000个，雷克子波主频为20 Hz，以200 m的炮间距在地表共布设51个炮点。采用全接收的观测系统，检波点共1001个，检波点间隔为10 m。

图5（a）、（b）为耦合波场成像结果，图5（c）~（f）为PP、PS、SP、SS成像结果。整体来看，图5（c）~（f）的同相轴能量分布更均衡，偏移假象更少，成像效果较好，可以较为清晰地实现侵入花岗岩体的基本构造成像。

为了方便分析解耦方程对成像结果精度的影响，截取图5中黑色标识框圈定的区域进行分析，将其左侧小黑框圈定区域记为局部1（图6），中间黑框圈定的区域记为局部2（图7）。

图6（a）、（b）和图7（a）、（b）为耦合波场在水平方向和垂直方向的局部成像结果，在白色箭头所指向的位置处，可以看到图6（a）、（b）; 图7（a）中出现目标体同相轴明显错段、能量微弱等现象，图7（b）中出现明显的同相轴合并，偏移假象等现象。都没有很好地对目标体成像，这可能会给后续的解释工作带来困难。而在图6（c）~（f）和图7（c）~（f）的同位置处，尤其是图6（f）和图7（f），可明显地看到同相轴不连续位置更少、同相轴能量更均衡、偏移假象更少、同相轴分辨率更高，可更准确地对目标体构造成像结果进行分析。

2.2.2 近代火山型

近代火山型干热岩主要是由深层高温岩浆向上涌出形成。图8为近代火山型的纵波速度分布情况。该模型左侧出现地堑构造，产生一条正断层和一条逆断层，含多套不规则形状地层，构造复杂。最小速度为2000 m/s，最大速度为5500 m/s。高温火山基岩位于5 km至深处，有2条主支左右入侵其他地层，岩浆房紧挨于火山基岩右旁，其中岩浆向上涌出至1 km处。

对该模型采样也是1001×701个，横纵空间采样间隔均为10 m，时间采样率为0.3 ms，时间采样20000个，子波主频为20 Hz，以200 m的炮间距在地表共布设51个炮点，同样采用的是检波点间隔为10 m的全接收观测系统。

由于岩株的速度相比其附近地层的要大很多，所以在1 km深度附近的岩株顶部发生了强烈反射，使得图9在该位置表现出同相轴的能量较强。相应地，透射波的能量相对减少，该区域至深处地震波能量较弱，目标体同相轴较为不明显。为了更好地分析基于解耦方程进行逆时偏移的成像效果，选择的区域最好是目标体同相轴能量较强的区域，故选择图9中黑色框圈定区域进行成像结果的局部分析，以进一步说明成像精度。

图10（a）、（b）为P-S波耦合的成像结果的局部，在箭头所指向的区域，目标体的同相轴能量较弱，有些地方几乎没有出现同相轴，成像结果不好，这是耦合波成像的常见缺陷。图10（c）~（f）为PP、PS、SP和SS波成像结果局部，相对于图10（a）、（b），同位置成像效果明显变好，同相轴有了较好的呈现，分辨率提高很多。

3 结论

（1）基于解耦方程的波场分解避免了传统Helmholtz分解对波场取散度和旋度导致的振幅减小、相位推迟的问题，无需后续额外校正，可在实现高精度逆时偏移成像的同时，将地震资料处理步骤减少。

（2）侵入花岗型和近代火山型两个典型干热岩模型的PP、PS、SP和SS成像结果说明，较于耦合波场成像，基于解耦方程分离纵横波得到的弹性波成像结果普遍更清晰，同相轴错段位置减少和信息更丰富，深部能量更均衡，分辨率更高。

参考文献