弹性波逆时偏移的角度域共成像点道集提取方法

2022-08-02徐蔚亚朱成宏魏哲枫张晓语杜启振

石油地球物理勘探 2022年4期

徐蔚亚朱成宏魏哲枫张晓语杜启振

(①页岩油气富集机理与有效开发国家重点实验室，北京 100083；②中国石化弹性波理论与探测技术重点实验室，北京 100083；③中国石化石油勘探开发研究院，北京 100083；④齐鲁工业大学(山东省科学院)山东省科学院海洋仪器仪表研究院，山东青岛 266061；⑤中国石油大学(华东)深层油气重点实验室，山东青岛 266580)

0 引言

逆时偏移(RTM)方法[1-4]于20世纪80年代提出，该方法基于波动方程实现波场延拓，以其无倾角限制、能对多种波(反射波、多次波、棱柱波)成像而受到广泛关注[5-6]。弹性波逆时偏移(ERTM)算法[7-8]由McMechan研究组提出，采用更接近真实的弹性介质和弹性波场，与传统的声波RTM成像方法相比，理论假设与实际更接近，并能同时实现纵波和横波的深度域成像。

地震波波场分离是ERTM方法中的重要环节，也是得到纯波模式成像结果的关键一步。除了基于Helmholtz定理的波场分离方法[9-12]以外，目前常用的是解耦延拓方程法[13-16]。解耦延拓方程法可在波场延拓过程中实现均匀介质的纵、横波质点位移的解耦或者质点振动速度的矢量解耦，但在应力解耦时，横波应力存在能量串扰问题[17-19]。为解决横波应力中的纵波串扰问题，Du等[16-17]提出了伪横波应力构建方法，实现了纵波和横波应力的解耦。

共成像点道集(CIG)作为叠前深度偏移的重要输出道集，可以用于改善成像质量，提供反映地下岩性的振幅和相位等信息，同时能用于速度建模。CIG[20]主要分为共炮检距道集(Common-offset CIG，COCIG)和角度域道集(Angle-domain CIG，ADCIG)。Sava等[21]利用单程波偏移提取地下局部COCIG，通过转换获得ADCIG。三维情况下，Xu等[22-23]在频率—波数域通过传播方向求取反射角信息和方位角信息，利用合适的成像条件提取角道集，优点是角道集质量较好，缺点是计算量和存储量均较大。王保利等[24]在实现过程中利用波场延拓获得的质点振动速度和应力，提取了成像点处的Poynting矢量，从而确定入射波和反射波的传播方向，最后通过构建成像条件提取了ADCIG，该方法简单易行；汪天池等[25]将该方法拓展到倾角域进行反射波与绕射波分离以实现绕射波RTM。Zhao等[26]通过提取地下COCIG进而提取了PP波和PS波的ADCIG，并分析了二者的运动学特征。Wang等[27]通过构建Poynting矢量提取了PP波和PS波的ADCIG，但该方法没有考虑分离得到的横波应力存在的纵波能量串扰问题。

针对横波应力串扰问题，本文基于解耦延拓方程及伪横波应力构建方法，开展了基于Poynting矢量提取横波反射角的研究，构建了基于ERTM的PP波和PS波ADCIG求取方法；针对PP波ADCIG不同角度的成像特性，提出了把ADCIG按角度衰减的叠加成像策略。最后应用模型数据验证了本文方法的正确性。

1 方法原理

1.1 解耦延拓方程

ERTM是利用弹性波动方程进行波场延拓，并结合波场分离方法得到解耦的震源端和检波端的纵/横波波场，然后运用弹性波成像条件得到最终的PP波和PS波成像结果。解耦延拓方程[14-19]能实现弹性波场的构建和解耦。其中，一阶弹性波波场方程[18]为

(1)

纵波方程为

(2)

横波波场可表示为

(3)

式中：vx、vz分别为总弹性波场中x和z方向的质点振动速度分量；vP,x、vP,z为纵波波场的两个方向速度分量；vS,x、vS,z为横波波场的两个方向速度分量；τxx、τzz、τxz为总弹性波场应力分量；τP为纵波体应力；λ、μ和ρ分别为弹性介质的拉梅常数和密度。

总弹性波场应力分量减去纵波体应力τP，可以构建横波应力分量[15，17]

(4)

式中τS,xx、τS,zz和τS,xz为横波应力分量。对上式进行时间积分，可得由弹性应变表示的横波应力

τS,ij=2μ(εij-bδij)

(5)

式中：εij为弹性应变分量；b=εxx+εzz，为弹性体应变。由于弹性应变既包含纵波应变εP,ij，又包含横波应变εS,ij，则有

εij=εP,ij+εS,ij

(6)

同时横波应力可以改写为

τS,ij=2μ[(εP,ij+εS,ij)-(bP+bS)δij]

(7)

式中：bP=εP,xx+εP,zz，为纵波体应变；bS=εS,xx+εS,zz，为横波体应变。可见横波应力除了横波应变作用2μ(εS,ij-bSδij)外，还有纵波应变作用2μ(εP,ij-bPδij)，因此，横波应力中含有纵波串扰。

为压制纵波能量在横波应力中的影响，Du等[16-17]基于横波质点振动速度提出了伪横波应力构建方程

(8)

式中：τqS,xx和τqS,zz分别为x方向和z方向的伪横波正应力分量；τqS,xz为伪横波切应力。基于横波质点振动速度构建的伪横波应力与横波应力在形式上保持一致。不同于横波应力，伪横波应力只与横波质点振动速度相关，避免了纵波应变对横波应力的影响。由于横波质点振动速度满足∂vS,z/∂z+∂vS,x/∂x=0，则伪横波应力满足τqS,xx=-τqS,zz。因此，基于解耦延拓方程和伪横波应力构建方法可以实现纵横波应力场的解耦，进而可以实现弹性应力和质点振动速度的解耦。

1.2 基于Poynting矢量提取PP波和PS波反射角

Poynting将能流密度矢量用于描述电磁能传播方向，即通过电场强度与磁场强度点乘得到电磁波的传播方向，后人把能流密度矢量称为Poynting矢量。Yoon等[28]给出了地震波场Poynting矢量表达式

S=-v·τ

(9)

式中：v=(vx，vz)为质点振动速度矢量；τ为应力张量。

利用解耦的质点振动速度及应力，可以确定纵波和横波的Poynting矢量，进而确定纵波和横波能量传播方向。

图1 纵波Poynting矢量与反射角的关系示意图

(10)

由于纵波反射角等于入射角，因此开角的一半即为纵波反射角。根据Snell定律，纵波入射角β和横波反射角α之间满足sinβ/VP=sinα/VS，其中VP、VS为纵、横波传播速度。

图2 纵波入射角β与横波反射角α关系示意图

(11)

式中

(12)

由VP、VS和开角θPS就可以确定PS波反射角α。

由于Poynting矢量的非稳健性，往往采用在一定空间窗上进行平滑的方式求取θPP或θPS，即

(13)

式中Ω表示平滑的空间窗口。

1.3 PP波和PS波角度域道集提取方法

获得纵波入射角β和横波反射角α后，便可以提取ADCIG以对成像过程进行质控，从而进一步提高成像精度。常规成像条件是对角道集所有角度的叠加，故不能直接用于提取ADCIG。为此，需对常规成像条件进行修改。借鉴王保利等[24]的方法，引入高斯采样函数，将常规成像条件拓展为PP波角度域共成像点成像条件

IPP_ADCIG(x,βk)

(14)

利用成像域PP波成像条件可以得到PP波的ADCIG，其中ERTM中PP波ADCIG的计算流程如下：

(1)在炮域构建震源波场，计算每一时刻的纵波质点振动速度场，并将其存储至硬盘；

(2)在炮域构建检波点波场，计算每一时刻的纵波质点振动速度场及应力场；

(3)从硬盘读取纵波质点振动速度场，利用式(9)计算当前时刻的Poynting矢量，并利用式(10)计算纵波的开角和入射角；

(4)利用式(13)进行角度平滑；

(5)以离散角作为角道集输出的角度变量，利用式(14)输出PP波ADCIG。

类似地，利用得到的横波反射角α，本文给出了PS波ADCIG成像条件

IPS_ADCIG(x,βk)

(15)

利用成像域PS波成像条件，可以得到PS波的ADCIG。其中ERTM中PS波ADCIG的计算流程如下：

(1)在炮域构建震源波场，计算每一时刻的纵波质点振动速度场，并将其存储至硬盘；

(2)在炮域构建检波点波场，计算每一时刻的横波质点振动速度场及应力场；

(3)从硬盘读取纵波质点振动速度场，利用式(9)计算当前时刻的Poynting矢量，并求取纵波入射及横波反射之间的开角θPS；

(4)结合Snell定律及波场传播关系，利用式(11)计算PS波反射角；

(5)应用式(13)进行角度平滑；

(6)以离散角作为角道集输出的角度变量，利用式(15)输出PS波ADCIG。

1.4 PP波ADCIG角度衰减+Laplace滤波的组合叠加成像策略

受双程波动方程影响，ERTM方法进行PP波成像过程中引入了低波数噪声。低波数噪声反射角接近90°，主要存在于角道集的大角度区域，因此对角道集进行叠加时，可通过对大角度切除获得消除低波数噪声的RTM结果。但是，道集中大角度道也含有部分有效信息。

借鉴Laplace滤波的物理机制，本文提出了在PP波角道集中小角度按纵波入射角度余弦cosβk的方式进行叠加成像、大角度实施Laplace滤波的组合叠加成像的策略，即

(16)

式中：Nβ2为离散角度个数；Nβ1为小角度与大角度的分界点序号。分界角取值范围一般为30°～40°。

根据式(16)可知，本文的ADCIG叠加成像的策略利用了角道集中所有的角度成像，小角度cosβk值更大，因此最终成像结果中既突出了角度道集中小角度道的贡献，又避免了大角度道有效信息的损失。

需要指出的是，对于PS波成像，由于震源端P波与检波端S波传播路径不同，且振动方向也不一致，因此PS波成像不存在低波数噪声。为此，PS波角道集叠加成像时可以忽略低波数噪声影响，不需要大角度道集噪声压制。

2 模型验证

2.1 水平层状模型

构建水平层状模型验证本文的叠加成像策略的正确性。模型如图3所示，网格数为400×350，尺寸为10m×10m，纵横波速度比为1.73。加载主频为15Hz的Ricker子波作为爆炸源，道间距为10m，炮间距为100m，在地表全覆盖，共接收40炮；时间采样间隔为1ms，记录长度为3.2s。采用时间二阶、空间十阶的有限差分算法进行数值模拟。

图3 水平层状模型

用式(14)和式(15)分别提取出水平层状模型的PP波角道集和PS波角道集，如图4所示。其中，每个角道集角度范围为0°～90°，间隔为1°；每隔50个CMP抽取一个角道集显示。从PP波成像道集可以看出，同相轴水平，埋深1500m、2500m界面在ADCIG中的同相轴连续性较好，成像噪声较少。由红色箭头处可以看出，PP低波数噪声主要位于角道集的浅层大角度的位置。从PS波成像角道集可看出，虽然有多次波干扰，但是埋深1500m界面在ADCIG中同相轴是水平的，而且连续性较好。这证明了PP波角道集和PS波角道集提取方法的正确性。

图4 水平层状模型的PP波(a)与PS波角道集(b)

为便于AVO/AVA分析，利用Zoeppritz方程[29]计算深度1500m界面的PP波和PS波的理论反射系数曲线。提取深度1500m、CMP200处PP波角道集和PS波角道集振幅曲线并归一化，并其与理论反射系数对比(图5)显示，提取的振幅随角度变化曲线与理论反射系数曲线基本一致，说明本文角道集提取方法基本可行。

图5 1500m处界面PP波(a)和PS波(b)角道集归一化振幅曲线与理论反射系数曲线对比

分别抽取水平层状模型PP波角道集中0°～40°、41°～70°及71°～90°的道进行叠加成像，分析角度对成像结果的影响(图6)。0°～40°的叠加成像结果(图6a)中几乎没有任何低波数噪声，41°～70°叠加成像结果(图6b)中低波数噪声较少，而71°～90°的成像结果(图6c)几乎全部是低频噪声，但有部分有效成像信息。

图6 水平层状模型PP波不同角度范围的叠加成像结果

采用本文方法提出的PP波角度衰减+Laplace滤波的组合叠加成像策略(式(16))，取分界角度为40°，得到组合叠加成像结果(图7a)，并抽取CMP200处的0°～40°和40°～90°角道集(图7b、图7c)，可见低波数噪声得以压制，且同相轴连续。

图7 水平层状模型PP波组合叠加成像结果(a)及CMP200处

2.2 Salt模型

利用Salt模型(图8)测试本文算法对复杂模型的成像精度和效果。模型网格数为1500×700，网格尺寸为10m×10m，纵横波速度比为1.73，密度取常数2.0g/cm3。数值模拟时，震源采用主频为30Hz的Ricker子波以爆炸形式加载至P波应力分量上；

图8 Salt纵波速度模型

时间采样间隔为1ms；炮点间隔为100m，共150炮；在地表以10m为间隔均匀分布的检波点进行接收。采用在时间二阶、空间十阶的有限差分算法进行正向和反向波场延拓。

提取的Salt模型的PP波的ADCIG如图9a所示，角度范围为0°～90°，角度间隔为2°，每隔100个CMP抽取一个角道集进行展示。利用本文提出的组合叠加策略得到的PP波成像结果如图9b所示。由图9a可以看出，每个ADCIG的同相轴都是水平的，且存在低波数噪声。采用组合叠加策略后，成像结果的低波数噪声得到有效压制。

图9 Salt模型PP波共成像点道集(a)和成像结果(b)

用本文提出的基于伪横波应力角道集方法提取的Salt模型PS波ADCIG如图10a所示。对0°～90°角道集进行叠加得到的PS成像结果如图10b所示。由图可以看出，每个ADCIG上的同相轴都是水平的，验证了本文PS波角道集提取方法的正确性。在浅层PS波成像结果优于PP波，在深层PP波成像结果优于PS波。

图10 Salt模型PS波共成像点道集(a)和成像结果(b)

图11a和图11b分别为基于伪横波应力(本文方法)提取的PS波部分CMP角道集及0°～45°叠加成像结果；图11c和图11d分别为基于横波应力提取的PS波部分CMP角道集及0°～45°的叠加成像结果。图11c中同相轴不连续，图11d中含有低波数噪声；图11a中同相轴较连续而且图11b中不存在低波数噪声。原因在于图11c和图11d中构造的横波应力(式(4))中含有纵波成分，导致提取的PS波角道集也含有纵波成分，因此，基于横波应力提取的PS波角道集中存在纵波能量干扰和低波数噪声；相比之下，图11a中基于伪横波应力(式(8))提取的PS波角道集较好压制了能量串扰。