各向异性弹性介质方向行波波场分离正演数值模拟

陈可洋

（中国石油大庆油田有限责任公司勘探开发研究院，黑龙江大庆163712）

0 引言

多波多分量弹性波波动方程的地震波正演数值模拟技术一直是国内外地球物理学界的热点。随着弹性波波动理论和计算机技术的不断发展，自20世纪60年代以来这项技术便得到了飞速发展。弹性波正演数值模拟技术能够基本保持弹性波的几何学、运动学和动力学等特征，因此可以达到精确模拟弹性波波场传播规律的目的［1］。目前，弹性波正演数值模拟技术的相关理论已基本成熟，并已形成了一系列经典算法和研究成果。例如，在数值离散方面［2-5］，形成了有限差分法、有限元法和伪谱法等方法，数值频散问题得到了有效压制；在处理人为截断边界方面［6-8］，形成了旁轴近似吸收边界、阻尼吸收边界及最佳匹配层（PML）吸收边界等各类吸收边界条件，提高了计算结果的信噪比；在计算过程稳定性方面［9-10］，形成了特征值分析法、矩阵分析法等平面波解的分析思路；在描述地球介质方面［11-13］，形成了均匀各向同性介质、各向异性介质、黏滞介质、双相介质、裂隙介质及多相孔隙介质等复杂的近似地球介质的相关理论的假设。弹性波正演数值模拟技术的引入不仅大大降低了物理模拟分析和野外多分量地震勘探的成本，而且有效地指导了多分量地震资料的采集、处理和解释，以及相关方法的验证等，因此，在当前的勘探地震学领域发挥着至关重要的作用。

随着地震波数值模拟技术的进步，如何在弹性波正演模拟过程中实现不同类型波场的自动分离逐渐得到广泛重视。其中，最具代表性的为弹性波波场分离数值模拟方法，其次为基于散度算子计算纯纵波和基于旋度算子计算纯横波的波场分离方法［14］，以及构建等效弹性波方程的波场分离数值模拟方法［11，13］。 同时，研究的等效弹性介质也已从目前的单相介质延伸到了双相孔隙介质中，而方向行波分离的方法早已得到广泛应用。例如，单程波叠前深度偏移方法需要对地震波波动方程进行算子分离处理，从而得到单程波偏移算子［15］；地震波数值模拟中的旁轴近似吸收边界也是采用时域单程波算子解决某空间方向的边界吸收问题［6］；还有上行波和下行波分离逆时成像条件可以对目的层有贡献的波场进行选择性成像［16］，从而有效压制逆时偏移噪声的形成。但是，目前尚未见到关于弹性波，特别是各向异性介质中方向行波波场分离方面的相关报道。

在前人研究的基础上，笔者以3个二维理论模型为例，在多波多分量各向异性介质波场模拟中引入波印廷矢量，开展上行波、下行波、左行波和右行波4个方向行波波场分离的弹性波正演数值模拟实验，以期提高对弹性波波场传播规律的认识，并指导多分量弹性波地震资料的数字处理。

1 基本理论

根据弹性波波动力学的相关理论［1］，二维一阶双曲型各向异性介质弹性波波动方程通常有如下的表达式［17］：

式中：vx和vz分别为质点振动速度的水平分量和垂直分量（二维矢量 U={vx，vz}），m/s；τxx，τzz和 τxz分别为应力的3个分量（x和z空间方向的正应力和切应力），N；Cij为各向异性介质弹性参数，109N/m3；ρ为介质密度，kg/m3。式（1）是根据牛顿第二定律和胡克定律来构建的［1］，其具体的数值离散方法、吸收边界条件及其稳定性条件可参见文献［18-20］。

另外，Thomson 参数［21］符号定义如下：

式中：ε和δ均为Thomson参数，无量纲。

同时，泊松比σ公式定义如下：

在此基础上［22-25］对波印廷矢量进行波场数值特征分析，并采用基于该矢量波场对原始多分量弹性波进行波场分离计算，最终实现了多分量各向异性介质弹性波左行波、右行波、上行波和下行波的波场分离数值模拟，并成功地应用于多个数值模拟算例中。

式中：σ为泊松比，无量纲。

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

2 数值模拟算例

2.1 均匀弹性各向异性介质模型

为了验证文中方法的准确性，设计了均匀弹性各向异性介质模型。模型总大小为1 km×1 km，纵横向空间步长均为5m。采用雷克子波波形的涨缩震源，其最大频率为40Hz，并置于模型中央位置处激发。采用的弹性参数C33为4×109N/m3，ε，δ和σ均为无量纲，其值分别为 0.2，0.5和 0.25，密度为103kg/m3，时间步长为0.5ms。满足计算所需的稳定性条件，波场快照记录时间为0.2 s。

图 1（a1）和图 1（a2）分别为采用公式（1）得到的水平分量和垂直分量各向异性介质弹性波波场快照。经分析可知，波场快照中存在2种波场，传播较快的是纵波，其波前面为椭圆形，传播较慢的是由各向异性引起的横波，其波前面形状不规则。在水平分量波场快照中，其左右两侧的波场极性相反；在垂直分量波场快照中，上下两侧的波场极性也相反。图1（b1）为水平方向的波印廷矢量波场快照，蓝色部分其数值为正，红色部分其数值为负，且正好沿着过震源的垂直线分为两瓣，即划分出左行波波场和右行波波场；图1（b2）为垂直方向的波印廷矢量波场快照，正好沿着过震源的水平线分为两瓣，即划分出上行波波场和下行波波场。于是根据图1（b1）和图 1（b2）对图 1（a1）和图 1（a2）的多分量弹性波波场快照进行分离，最终得到了分离后的左行波波场［图 1（c1）和图 1（c2）］、右行波波场［图 1（d1）和图 1（d2）］、上行波波场［图 1（e1）和图 1（e2）］和下行波波场［图 1（f1）和图 1（f2）］，且分离结果准确可靠。图1（g1）和图1（g2）分别为正应力和剪应力的波场快照。在以震源为原点的4个象限中，正应力波场不存在极性反转现象，而剪应力波场存在极性反转现象，对角线方向极性相一致。综合分析认为，文中方法准确实现了均匀、多波多分量弹性各向异性介质中方向行波的波场分离处理。

图1 0.2 s时刻均匀介质弹性波行波波场分离的波场快照及其波印廷矢量Fig.1 Elastic one-waywave field separating snapshots in isotropicmedia and the corresponding Poynting vector at0.2 s

2.2 倾斜界面速度模型

为了验证文中方法在含不同倾角地质情况下的准确性和有效性，以含多个角度的倾斜界面速度模型（图2）为例进行分析。模型大小为2.5 km×1.0 km，纵横向空间步长均为5m；模型顶层Ⅰ左侧部分和底层Ⅵ厚度均为250m，模型顶层Ⅰ右侧部分最大厚度为750m。模型中间部分包含一组不同倾角的倾斜界面，倾斜界面自右向左的倾角依次为30°，45°，60°和 90°，对应的弹性参数 C33依次是：Ⅰ层为 3.24×109N/m3，Ⅱ层为 5.29×109N/m3，Ⅲ层为6.76×109N/m3，Ⅳ层为 9×109N/m3，Ⅴ层为 7.29×109N/m3，Ⅵ层为 12.25×109N/m3，密度均为 103kg/m3，ε，δ和σ均为无量纲，分别为0.2，0.5和0.25。采用雷克子波波形的涨缩震源，其最大频率为40Hz，在模型横向1 km和纵向0.2 km处激发，合成记录的道长为1 s，时间步长为0.5ms，满足计算所需的稳定性条件。地面采集深度为100m，VSP采集位置距离模型左侧100m处，波场快照记录时间为0.2 s。

图2 各向异性弹性介质速度模型Fig.2 Anisotropic elasticmedium velocitymodel

图 3（a1）和图 3（a2）分别为 0.2 s时刻水平分量和垂直分量各向异性介质弹性波波场快照。图3（b1）和图3（b2）分别为水平方向和垂直方向的波印廷矢量波场快照，蓝色部分波场其数值为正，红色部分波场其数值为负。分析图3可知，在含倾斜界面的复杂速度模型中，左行波、右行波、上行波和下行波波场的传播特征较为复杂，但基于波印廷矢量可以对不同类型的行波波场加以自动识别和区分。图 3（c）～（f）分别为根据图 3（b1）和图 3（b2）对图 3（a1）和图3（a2）的多分量弹性波波场快照进行的分离处理结果，最终得到了完全分离的水平分量和垂直分量的左行波［图 3（c1）和图 3（c2）］、右行波［图3（d1）和图 3（d2）］、上行波［图 3（e1）和图 3（e2）］和下行波［图 3（f1）和图 3（f2）］波场快照。 分析分离后的多分量弹性波波场快照可知，文中方法准确实现了较为复杂介质情况下的上行波、下行波、左行波和右行波的波场分离处理。图3（g1）和图3（g2）分别为正应力和剪应力的波场快照。图3（h1）和图3（h2）分别为地面采集的水平分量和垂直分量弹性波模拟记录。图3（i）～（l）分别为对应的水平分量和垂直分量左行波［图 3（i1）和图 3（i2）］与右行波［图 3（j1）和图 3（j2）］。 图 3（k1）和图 3（k2）分别为正应力和剪应力的弹性波模拟记录。图3（l1）和图3（l2）分别为VSP采集的水平分量和垂直分量弹性波数值模拟记录。图3（m）～（n）分别为对应的水平分量和垂直分量上行波［图 3（m1）和图 3（m2）］与下行波［图 3（n1）和图 3（n2）］。 图 3（o1）和图 3（o2）分别为VSP采集的正应力和剪应力的数值模拟记录。综上所述，文中方法在较为复杂的多分量各向异性介质弹性波数值模拟过程中准确实现了上行波、下行波、左行波和右行波的波场分离处理。

图3 0.2 s时刻多分量各向异性介质弹性波行波分离波场快照、模拟记录及其波印廷矢量Fig.3 Elastic one-way wave separating snapshots and numerical records in multi-components anisotropic mediaand the corresponding Poynting vector at 0.2 s

2.3 M armousi模型

为了进一步验证文中方法在更为复杂的弹性各向异性介质模型中的适用性及其应用效果，以Marmousi模型为例进行了分析。模型总大小为3.4 km×1.4 km，纵横向空间网格均为5m；最小弹性参数 C33为 1.06×109N/m3，最大弹性参数 C33为21.81×109N/m3，密度均为 103kg/m3，ε，δ和 σ 均为无量纲，分别为0.2，0.5和0.25。采用雷克子波波形的涨缩震源，其最大频率为40Hz，在横向1.7 km和纵向0.15 km处激发，合成记录的道长为2s，时间步长为0.2ms，满足计算所需的稳定性条件。地面采集深度为100m，波场快照记录时间为0.6 s。

图4 0.6 s时刻M armousi模型多分量各向异性介质弹性波行波分离波场快照、模拟记录及其波印廷矢量Fig.4 Elastic one-way wave separating snapshots and numerical records in multi-components anisotropic mediaMarmousi model and the corresponding Poynting vector at 0.6 s

图 4（a1）和图 4（a2）分别为 0.6 s时刻水平分量和垂直分量各向异性介质弹性波波场快照。图4（b1）和图4（b2）分别为水平方向和垂直方向的波印廷矢量波场快照，蓝色部分波场其数值为正，红色部分波场其数值为负。分析图4可知，在基于Marmousi速度模型的波场快照中，左行波、右行波、上行波和下行波的波波场传播特征更为复杂，存在各种类型的弹性波场转换关系，但基于波印廷矢量可对不同类型的行波波场加以识别和分离。图4（c）～（f）分别为根据图 4（b1）和图 4（b2）对图 4（a1）和图4（a2）的多分量弹性波波场快照进行的分离处理结果，最终得到了完全分离的水平分量和垂直分量左行波［图 4（c1）和图 4（c2）］、右行波［图 4（d1）和图 4（d2）］、上行波［图 4（e1）和图 4（e2）］与下行波［图 4（f1）和图 4（f2）］波场快照。 分析分离后的多分量各向异性介质弹性波波场快照可知，文中方法也准确实现了上行波、下行波、左行波和右行波的波场分离处理。图4（g1）和图4（g2）分别为正应力和剪应力的波场快照。图4（h1）和图4（h2）分别为地面采集的水平分量和垂直分量弹性波模拟记录。图4（i）～（j）分别为对应的水平分量和垂直分量左行波［图 4（i1）和图 4（i2）］与右行波［图 4（j1）和图 4（j2）］波场的数值模拟记录。 ［图 4（k1）和图 4（k2）］分别为地面采集的正应力和剪应力的数值模拟记录。综合分析可知，本文方法准确实现了复杂弹性各向异性介质上行波、下行波、左行波和右行波的波场分离处理。

3 结论

（1）基于多分量弹性波波印廷矢量，实现了多分量各向异性介质弹性波上行波、下行波、左行波和右行波的方向行波波场分离数值模拟。该方法计算量小，计算过程简单，易实现。通过3个不同的理论模型与数值模拟，验证了该方法的准确性和有效性。

（2）通过弹性波方向行波的波场分离数值模拟，进一步提高了对弹性波波场传播规律的认识，同时还能够为多波多分量弹性波资料的方法验证提供指导和帮助。

（3）由于多波多分量弹性波逆时偏移理论仍处于理论研究和探索阶段，该方法目前尚未在逆时成像方面进行应用，但这种方向行波分离的思路构建的优化逆时成像条件在声波波动方程中已取得了较好的应用效果，因此，该方法在多波多分量地震资料处理中具有重要的应用价值。

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

［2］ Dablain M A.Theapplication ofhigh-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，28（4）：1099-1105.

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

［6］ 张秉铭.各向异性介质中弹性波数值模拟与偏移研究［D］.北京：中国科学院地球物理研究所，1997.

［7］ 杨微，陈可洋.加权吸收边界条件的优化设计［J］.石油物探，2009，48（3）：244-246.

［8］ 张会星，宁书年.弹性波动方程叠前逆时偏移［J］.中国矿业大学学报，2002，31（5）：371-375.

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

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

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

［12］ 宁刚，熊章强，陈持逊.波动方程有限差分正演模拟误差来源分析［J］.物探与化探，2008，32（2）：203-206.

［13］ 陈可洋.高精度地震纯波震源数值模拟［J］.岩性油气藏，2012，24（1）：84-91.

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

［15］ DenliH，Huang Lianjie.Elastic-wave reverse-timemigrationwith a wavefield-separation imaging condition［C］.SEG Las Vegas2008 AnnualMeeting，2008：2346-2350.

［16］ 陈可洋，陈树民，李来林，等.弹性波联合叠前逆时偏移数值试验［J］.石油物探，2014，53（1）：8-16.

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

［18］ 裴正林，牟永光.地震波传播的数值模拟［J］.地球物理学进展，2004，19（4）：933-941.

［19］ 陈可洋.地震波数值模拟中差分近似的各向异性分析［J］.石油物探，2010，49（1）：19-22.

［20］ 陈可洋，吴清岭，范兴才，等.高阶高密度三维多波多分量弹性波波场分离正演数值模拟［J］.油气藏评价与开发，2013，3（2）：6-14.

［21］ Thomsen L.Weak elastic anisotropy［J］. Geophysics，1986，51（10）：1954-1966.

［22］ Yoon K，Marfurt Kurt J，Starr W. Challenges in reverse-time migration［G］. SEG 74th Annual International Meeting Expanded Abstracts，2003：1057-1060.

［23］ 陈可洋.几种地震观测方式的逆时成像分析［J］.岩性油气藏，2013，25（1）：95-101.

［24］ 陈可洋，吴沛熹，杨微.扩散滤波方法在地震资料处理中的应用研究［J］.岩性油气藏，2014，26（1）：117-122.

［25］ 陈可洋，杨微，吴清岭，等.地震反射波与散射波波场分离方法初探［J］.岩性油气藏，2013，25（2）：76-81.