地震波动方程方向行波波场分离正演数值模拟与逆时成像

陈可洋，陈树民，李来林，吴清岭，范兴才，刘振宽

（1.中国石油大庆油田有限责任公司 勘探开发研究院，黑龙江大庆163712；2.中国石油大庆油田有限责任公司勘探事业部，黑龙江 大庆163453）

0 引言

地震波正演数值模拟技术一直是国际地球物理学界的热点内容之一。随着地震波动理论和计算机技术的不断发展，自20世纪60年代以来该项技术便得到了飞速发展，其中采用波动方程的地震波数值模拟技术能更好地保持地震波的几何学、运动学和动力学等特征，因此可达到精确模拟地震波传播规律的目的［1］。目前，地震波正演数值模拟技术的相关理论已基本成熟，并已形成了一系列经典算法和研究成果，如在数值离散方面［2-5］，以交错网格有限差分法最具代表性，有效压制了数值频散并提高了计算效率；在吸收边界方面［6-8］，以最佳匹配层（PML）吸收边界条件最具代表性，在减小计算量的同时减少了虚假波场的引入，有效提高了计算结果的信噪比。地震波正演数值模拟技术的引入不仅大大降低了物理模拟技术和实际野外地震勘探的成本，同时还可有效指导实际地震资料的采集、处理、解释及方法验证等，因此，该技术在当前勘探地震学领域中发挥着至关重要的作用。

随着地震波数值模拟技术的进步，如何在地震波正演模拟过程中实现不同类型波场的自动识别与分离问题逐渐得到广泛重视，有利于对复杂波场传播规律的认识和应用。其中，最具代表性的是弹性波波场分离数值模拟方法，目前该方法包括基于散度算子计算纯纵波和基于旋度算子计算纯横波的波场分离数值模拟方法［9］，也有构建等效弹性波方程的纵横波波场分离数值模拟方法［10］，它们均可实现纯纵波和纯横波波场的自动分离。同时，研究介质也从单相介质延伸到双相孔隙介质。截至目前，未见方向行波波场分离数值模拟方面的相关报道，但针对地震方向行波波场分离的相关方法早已得到应用，如利用单程波偏移算子进行地震成像［11］和利用单程波吸收边界条件进行地震波正演中的边界吸收处理［12］等。在前人研究的基础上，笔者采用3种模型实例，基于波印廷矢量，开展地震波场中上行波、下行波、左行波和右行波共4种类型方向行波的波场分离数值模拟实验，以期提高对地震波波场传播规律的认识，并指导实际地震资料数字处理。

1 基本理论

地震波波动方程通常可表述为一阶双曲型形式［6］

式中：P为质点的应力波场，无量纲；u和w分别为2个空间方向的质点振动速度，m/s；v为介质速度，m/s；f为震源函数，无量纲；xs和zs分别为震源在模型x和z空间方向的坐标，m；t为波场传播时间，s；ρ为密度，kg/m3。

式（1）的具体数值离散方法、吸收边界条件及其稳定性条件可参见文献［1］。

Yoon等［11］给出了用于地震波场的波印廷矢量计算公式，即

该波印廷矢量代表地震波在介质中传播时其波前面的法线传播方向。笔者在此基础上对波印廷矢量进行了波场数值特征分析，并基于该矢量波场实现对原始波场的分离，最终实现了左行波、右行波、上行波和下行波的波场分离正演数值模拟，且计算过程简单易行。

2 方向行波分离数值算例

2.1均匀介质模型

模型大小为1 km×1 km，横、纵向空间步长均为5m。采用最大频率为60Hz的雷克子波作为震源，并置于模型中央位置处激发。地震波速度为2 000 m/s，介质密度为 1.1×103k/m3，时间步长为0.5ms，满足计算所需的稳定性条件［7］，波场快照记录时间为0.2 s。

分析图1（a）可知，在各向同性介质中，地震波波场快照的波前面为圆形。图1（b）为水平方向波印廷矢量波场快照，蓝色波场数值为正，红色波场数值为负，且正好沿过震源的垂直线分为2瓣，即划分为左行波波场和右行波波场；同理，图1（c）为垂直方向波印廷矢量波场快照，正好沿过震源的水平线分为2瓣，即划分为上行波波场和下行波波场。根据图 1（b）和图 1（c）对图 1（a）完整的地震波波场快照进行筛选，最终得到了完全分离出来的左行波波场［图 1（d）］、右行波波场［图 1（e）］、上行波波场［图 1（f）］和下行波波场［图 1（g）］。根据计算结果分析可知，在各向同性介质中本文所提的方向行波波场分离数值模拟方法准确可靠。

图1 0.2 s时刻均匀介质地震行波分离的波场快照及其波印廷矢量Fig.1 Seismic one-way wave separating snapshots in isotropic media and corresponding Poynting vector at 0.2 s

2.2 倾斜界面速度模型

图2 复杂介质模型Fig.2 Complex velocity model

图2为含有多个角度的倾斜界面速度模型，模型大小为2.5 km×1.0 km，横、纵向空间步长均为5m。模型顶层①和底层⑦厚度分别为150m和250m，对应的地震波速度分别为1 800m/s和3 500m/s。模型中间部分包含一组不同倾角的倾斜界面，而且倾斜界面自右向左的倾角依次为 30°，45°，60°和90°，模型中介质②～⑥对应的地震波速度分别为2 000m/s，2 300m/s，2 600m/s，3 000m/s 和 2 700m/s，介质密度均为1.1×103k/m3。采用最大频率为60Hz的雷克子波作为震源，在模型（1.0 km，0.2 km）位置处激发，合成记录长度为1 s，时间步长为0.5ms，满足计算所需的稳定性条件［7］。地面采集深度为100m，VSP采集位置距模型左侧100m，波场快照记录时间为0.2 s。

图3（a）和图3（b）分别为水平方向和垂直方向的波印廷矢量波场快照，其中蓝色部分波场数值为正，红色部分波场数值为负。分析图3（a）和图3（b）可知，在这种较为复杂的速度模型中，不同方向地震行波波场的传播特征非常复杂，无法有效识别哪些波场属左行波波场，哪些波场属右行波波场，哪些波场属上行波波场以及哪些波场属下行波波场，而基于波印廷矢量可以对不同类型的地震方向行波波场实现自动识别与分离。图3（c）为0.2 s时刻的地震波全波场快照，图 3（d）～（g）分别为对图 3（c）的全波场快照进行方向行波波场分离的结果，最终得到了完全分离的左行波、右行波、上行波和下行波波场快照。图3（h）和图3（j）分别为100m深度位置采集和VSP观测方式采集的模拟记录。从这2种模拟记录中无法有效识别与分离不同类型的方向行波波场。 图 3（i）和图 3（k）分别为图 3（h）和图3（j）对应的波印廷矢量数值模拟记录，图 3（l）～（o）分别为该方法分离出的左行波、右行波、上行波和下行波波场数值模拟记录。综合分析可知，该方法适合较为复杂的介质模型，而且可以获得准确的左行波、右行波、上行波和下行波波场分离结果。

图3 复杂介质模型中地震方向行波波场分离的波场快照、模拟记录及其波印廷矢量Fig.3 Seismic one-way wave separating snapshots, numerical records and corresponding Poyntingvector in complex velocity model

2.3 M armousi模型

为了进一步验证方向行波波场分离数值模拟方法在更为复杂的介质模型中的应用效果，以Marmousi模型进行模拟。模型大小为3.4 km×1.4 km，横、纵向空间步长均为5m。地震波最小速度为1 028m/s，最大速度为4 670m/s，采用最大频率为60Hz的雷克子波作为震源，在（1.70 km，0.15 km）处激发，合成记录长度为2 s，时间步长为0.2ms，满足计算所需的稳定性条件。地面采集深度为100m，VSP采集位置距模型左侧1 km，波场快照记录时间为0.6 s。

图 4（a）和图 4（b）分别为 Marmousi模型水平方向和垂直方向的波印廷矢量波场快照。分析图4（a）和图4（b）可知，由于介质的复杂性，在该介质的波场快照中左行波、右行波、上行波和下行波波场的传播特征变得更为复杂，而基于波印廷矢量的地震波波场仍可有效地分离出不同类型的方向行波波场。图 4（c）为 0.6 s时刻的地震波全波场快照，图 4（d）～（g）分别为对图4（c）的完整地震波波场快照进行方向行波波场的分离处理，最终得到了完全分离的左行波、右行波、上行波和下行波波场快照。图4（h）和图4（j）分别为Marmousi模型100m深度位置采集和VSP观测方式采集的地震波正演模拟记录。从这2种模拟记录中也无法有效识别和分离不同类型的方向行波波场。图4（i）和图4（k）分别为图4（h）和图 4（j）对应的波印廷矢量模拟记录，图 4（l）～（o）分别为该方法分离出的左行波、右行波、上行波和下行波的波场模拟记录。Marmousi模型数值分析表明，该方法的应用效果也准确有效。

图4 M armousi模型中地震方向行波分离的波场快照、模拟记录及其波印廷矢量Fig.4 Seismic one-way wave separating snapshots, numerical records and the corresponding Poynting vector in Marmousi model

2.4 行波分离方法在逆时成像中的应用

以图2速度模型为例，震源采用最大频率为60 Hz的雷克子波，从模型最左侧5m处开始向右侧激发，震源深度均为10m，炮间距为50m，共激发50炮，检波器埋置深度为10m，道间距为5m，共501个检波器。合成记录长度为2 s，时间步长为0.5ms，满足计算所需的稳定性条件［7］。采用交错网格有限差分法，其差分精度分别为时间2阶和空间10阶，在边界处采用20个网格点的PML吸收边界条件来压制边界反射波，可保证有效模拟区域内具有较高的信噪比和人工截断处具有较好的边界吸收效果，然后采用相关法逆时成像条件对采集到的所有炮记录进行逆时成像处理。

图5（a）为炮点和检波点处上行波和下行波传播方向示意图，S代表炮点波场，R代表检波点波场，原因在于100m深度位置检波器仅能采集这2种方向行波波场。分析波场传播机理可知，检波点的下行波波场（逆向传播）其实与炮点上行波波场的传播路径一致，检波点的上行波波场与炮点下行波波场的传播路径一致，将这些波波场两两作互相关法逆时成像得到的是逆时偏移噪声，而检波点的下行波波场与炮点下行波波场的传播路径以及检波点的上行波波场与炮点上行波波场的传播路径不同，因此将这些波波场两两作互相关法逆时成像得到的是真实的地层有效信号［16］。图 5（b）为常规逆时偏移成像剖面，其成像结果为

式中：*代表相关成像。

图5 方向行波示意图及其波场分离方法在逆时成像中的应用Fig.5 Directional one-way wave maps and its application in reverse-time migration

分析图 5（b）和图 5（c）可知，应用方向行波波场分离后的逆时成像剖面，其偏移噪声得到了有效压制，信噪比更高。由此可见，基于方向行波的地震波场分离方法在实际地震资料复杂构造成像处理中具有较好的应用效果。

3 结论

（1）根据波印廷矢量在地震波场中的特点，将其应用于地震方向行波的波场分离正演数值模拟中，实现对常规地震波场中的左行波、右行波、上行波和下行波波场进行自动识别与分离，而且该方法计算量小，计算过程简单易行。

（2）通过3种不同模型的数值模拟，验证了该方法的准确有效性，并且得到了完全分离的左行波、右行波、上行波和下行波4种方向行波的波场快照和数值模拟记录。通过地震方向行波的波场分离处理，可进一步提高对地震波场的传播过程及规律的认识，同时还能为实际地震资料和VSP资料的方法验证提供指导与帮助。

（3）通过地震方向行波波场分离数值模拟在逆时成像中的应用，其成像效果得到显著改善，偏移噪声得到有效压制，因此，该方法在实际地震资料成像等处理中具有重要的研究和应用价值。

［1］ 陈可洋.高阶弹性波波动方程正演模拟及逆时偏移成像研究［D］.大庆：大庆石油学院，2009.

［2］ Dablain M A.The application of high-order differencing to the scalar waveequation［J］.Geophysics，1986，51（1）：54-66.

［3］ 张中杰，腾吉文，杨顶辉.声波与弹性波数值模拟中的褶积微分算子法［J］.地震学报，1996，18（1）：63-69.

［4］ 程冰洁，李小凡.2.5维地震波场褶积微分算子法数值模拟［J］.地球物理学进展，2008，23（4）：1099-1105.

［5］ 朱生旺，魏修成.波动方程非规则网格任意阶精度差分法正演［J］.石油地球物理勘探，2005，40（2）：149-153.

［6］ 陈可洋.边界吸收中镶边法的评价［J］.中国科学院研究生院学报，2010，27（2）：170-175.

［7］ 廖振鹏，周正华，张艳红.波动数值模拟中透射边界的稳定实现［J］.地球物理学报，2002，45（4）：533-545.

［8］ 马在田.地震成像技术有限差分法偏移［M］.北京：石油工业出版社，1989.

［9］ 刘洪林，陈可洋，杨微，等.高阶交错网格有限差分法纵横波波场分离数值模拟［J］.地球物理学进展，2010，25（3）：877-884.

［10］ 陈可洋.一阶速度-应力Biot双相各向同性介质弹性波波场分离数值模拟［J］.计算物理，2011，28（3）：404-412.

［11］ YYoon K，Marfurt K J，Starr W. Challenges in reverse-time migration［C］. 74th SEG Annual International Meeting，Denver，Colorado，2004：1057-1060.

［12］ Sava P C，Formel S. Angle-domain common-image gathers by wavefieldcontinuation methods［J］. Geophysics，2003，68（14）：1065- 1074.

［13］ Fomel S. Theory of 3-D angle gathers in wave-equation seismic imaging［J］. Journal of Petroleum Exploration and Production Technology， 2011，1：11-16.

［14］ Sava P，Formel S. Time-shift imaging condition in seismic migration ［J］. Geophysics，2006，71（6）：209-217.

［15］ Xu Sheng，Zhang Yu，Tang Bing. 3D angle gathers from Reversetime migration［J］. Geophysics，2011，76（2）：77-92.

［16］ 陈可洋，吴清岭，范兴才，等.地震波叠前逆时偏移脉冲响应研究与应用［J］.石油物探，2013，52（2）：163-170.