绕射波与一次波联合黏声最小二乘逆时偏移

doi:10.11720/wtyht.2023.2636

绕射波与一次波联合黏声最小二乘逆时偏移

徐雷良¹, 赵国勇¹, 张剑², 钟天淼³, 谷佳莹³, 游剑³, 曲英铭^,³

1.中石化石油工程地球物理有限公司科技研发中心,江苏南京 211100

2.中石油化工股份有限公司石油工程地球物理公司胜利分公司,山东东营 257100

3.中国石油大学(华东) 地球科学与技术学院,山东青岛 266580

Joint Q-compensated least-squares reverse time migration using primary and diffracted waves

XU Lei-Liang¹, ZHAO Guo-Yong¹, ZHANG Jian², ZHONG Tian-Miao³, GU Jia-Ying³, YOU Jian³, QU Ying-Ming^,³

1. Science and Technology Research and Development Center,Sinopec Petroleum Engineering Geophysical Limited Corporation,Nanjing 211100,China

2. Shengli Branch Company,Sinopec Petroleum Engineering Geophysical Limited Corporation,Dongying 257100,China

3. School of Geosciences,China University of Petroleum (East China),Qingdao 266580,China

通讯作者: 曲英铭(1990-),男,博士,中国石油大学(华东)副教授,博士生导师,主要从事地震波传播、成像与反演等方面的研究工作。Email:quyingming@upc.edu.cn

责任编辑: 叶佩

收稿日期: 2021-12-11 修回日期: 2022-05-27

基金资助:

国家自然科学基金(42174138)
国家自然科学基金(42074133)
中国石化科技项目(P22165)

Received: 2021-12-11 Revised: 2022-05-27

作者简介 About authors

徐雷良(1983-),男,高级工程师,主要从事地震资料采集与处理技术研究工作。Email:sl-xull.osgc@sinopec.com

摘要

由于照明不足,小尺度断层和孔洞的成像是一个难题,地下衰减导致地震波振幅损失和相位畸变,在成像过程中忽略这种衰减会造成偏移振幅模糊。基于黏声最小二乘逆时偏移(QLSRTM)能够优化小尺度构造的成像,但这需要大量的迭代和计算成本。为了提高小尺度构造的成像效果,提出了一种充分使用绕射波的面向地质目标的的黏声LSRTM(J-QLSRTM)。在该方法中,构建了新的目标函数和梯度公式,并基于反演理论和伴随理论,推导了一次波和绕射波的Q补偿波场传播算子、Q补偿伴随算子和Q衰减反偏移算子。数值实例证明了J-QLSRTM比传统QLSRTM和声波J-LSRTM更有优势。

关键词： 最小二乘逆时偏移; 黏声; 绕射波

Abstract

Poor illumination poses great challenges to the imaging of small-scale faults and pores.Subsurface attenuation leads to amplitude loss and phase distortion of seismic waves,and ignoring such attenuation during imaging will blur migration amplitudes.The Q-compensated least-squares reverse time migration (QLSRTM) can improve the imaging of these small-scale structures,but it requires a huge amount of iterations and computational cost.To improve the imaging effect of these small-scale structures,this study proposed a geological-target-oriented joint QLSRTM (J-QLSRTM) that fully utilizes diffracted waves.In this method,a new objective function and gradient formula was constructed.Moreover,the Q-compensated wavefield propagation operators,Q-compensated adjoint operators,and Q-attenuated demigration operators were derived for both primary and diffracted waves based on the inversion and adjoint theories.The numerical examples verified that the proposed J-QLSRTM is superior to the conventional QLSRTM and the acoustic J-LSRTM.

Keywords： least-squares reverse time migration; viscoacoustic; diffracted waves

PDF (4481KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

徐雷良, 赵国勇, 张剑, 钟天淼, 谷佳莹, 游剑, 曲英铭. 绕射波与一次波联合黏声最小二乘逆时偏移[J]. 物探与化探, 2023, 47(1): 91-98 doi:10.11720/wtyht.2023.2636

XU Lei-Liang, ZHAO Guo-Yong, ZHANG Jian, ZHONG Tian-Miao, GU Jia-Ying, YOU Jian, QU Ying-Ming. Joint Q-compensated least-squares reverse time migration using primary and diffracted waves[J]. Geophysical and Geochemical Exploration, 2023, 47(1): 91-98 doi:10.11720/wtyht.2023.2636

0 引言

小尺度断层和孔洞给高分辨率成像带来了难题。在复杂勘探区域的成像中,基于双程波方程的逆时偏移(RTM)比基于射线方程和单程波方程的偏移更有优势。尽管如此,由于观测系统的限制,传统RTM很难对这些复杂地下构造进行成像。为了克服传统偏移方法中的这些问题,基于最小二乘反演和伴随状态方法的最小二乘逆时偏移法(LSRTM)被提出并迅速发展^[1]。然而,该方法需要大量的迭代和计算成本才能生成高分辨率的小尺度断层和孔洞。

地下任何小于或能与地震波波长相比拟的不均匀体(如断块、断棱)均可看作是绕射点,当地震波通过绕射点时产生的波称为绕射波,地震波通过大尺度构造产生的波称为反射波,通常利用一次反射波进行成像,而将绕射波进行压制。但绕射波中蕴含着丰富的与断层、尖灭、盐丘以及缝洞型储集体等非均质性构造有关的信息,因此,含有特殊传播路径的绕射波应该能够被充分用于提高一次波成像而不是压制。绕射成像方法有两种:直接绕射成像和间接绕射成像。直接成像方法主要包括Kirchhoff法^[2⇓-4]、叠加法^[5]、角域法^{[6⇓⇓⇓-10]}。该方法根据绕射波和反射波的能量差异,在成像过程中直接分离绕射的能量。间接方法实现绕射成像的步骤如下:首先从叠前道集(例如共炮点道集、平面波道集、共偏移距道集等)中分离出绕射波,然后将绕射波场单独成像。该方法主要包括共偏移距道集法^[11]、共绕射点剖面法^[12]、共炮点记录法^[3,13-14]、平面波记录法^[15]、倾斜叠加法^[16]以及Radon变换法^[17]。

然而,地球具有严重的粘滞性特性。由粘弹性引起的地震波衰减会在成像中导致振幅减弱和反射层错位^[18-19]。如果在成像中忽略地下衰减,将会在偏移成像中导致更严重的振幅损失和相位畸变。为了补偿衰减,在成像中主要考虑两种方法。一种为反Q滤波法,该方法有效且计算效率高^[20⇓-22],但它不能在复杂的地质构造中正确校正衰减;另一种为Q补偿偏移法^[23⇓-25]。近年来,一些基于各种黏声拟微分方程的Q补偿LSRTM(QLSRTM)逐渐发展起来^[26⇓-28]。

为了提高小尺度断层和孔洞的成像效果,提出一种面向目标结构的绕射波和一次波联合方法QLSRTM(J-QLSRTM),并通过典型模型来检验绕射波场传播算符的精度,然后应用衰减Sigsbee2B模型验证了它相对于QLSRTM和无补偿J-LSRTM的优势,并通过野外数据实例,对J-QLSRTM方法的应用灵活性进行了评价。

1 方法原理

1.1 J-QLSRTM的目标函数

包含一次波和绕射波的观测衰减地震数据( $d_{Q}^{o b s}$ )可由下式表示:

(1)

d_{Q}^{o b s}

d_{Q_{0}}^{o b s}

d_{Q_{r}}^{o b s}

式中: $d_{Q_{0}}^{o b s}$ 和 $d_{Q_{r}}^{o b s}$ 分别为一次波和绕射波的观测衰减数据。对于常规LSRTM,目标函数E为:

(2)E(m)=

\frac{1}{2} {({L^{S}}_{Q_{0}} m - d_{Q}^{o b s})}_{2}^{2}

\frac{1}{2} {(({L^{S}}_{Q_{0}} m - d_{Q_{0}}^{o b s}) - d_{Q_{r}}^{o b s})}_{2}^{2}

式中: ${L^{S}}_{Q_{0}}$ 为一次波衰减反偏移算子(线性模型算子); $‖ ‖_{2}$ 为L₂范数;m为反射模型。

因此,由式(2)可导出梯度公式

(3)g=

\underset{(Ⅰ)}{\underset{⏟}{{L^{R}}_{Q_{0}} ({L^{S}}_{Q_{0}} m - d_{Q_{0}}^{o b s})}}

\underset{(Ⅱ)}{\underset{⏟}{{L^{R}}_{Q_{0}} d_{Q_{r}}^{o b s}}}

。

式中: ${L^{R}}_{Q_{0}}$ 为一次波伴随算子。由式(3)可知,梯度由两部分组成:真梯度部分(I)和绕射波观测数据的噪声梯度部分(II),随着LSRTM的迭代,这些噪声会变得更强。为了克服这一问题,QLSRTM采用了绕射波和一次波相结合的方法,通过提高小尺度孔洞的成像效果以获得更好的分辨率。

将所提出的J-QLSRTM新目标函数构造为

(4)E(m)=

\frac{1}{2} {(a_{0} ({L^{S}}_{Q_{0}} m - d_{Q_{0}}^{o b s}) + a_{1} ({L^{S}}_{Q_{r}} m - d_{Q_{r}}^{o b s}))}_{2}^{2}

式中:a₀、a₁为系数; ${L^{S}}_{Q_{r}}$ 为绕射波的衰减偏移算子。从式(4)可推导出梯度公式(推导过程见附录A)

(5)g=

\frac{1}{δ v} (w_{R_{0}} F_{C_{0}} {U^{R}}_{C_{0}} + w_{R_{1}} F_{C_{r}} {U^{R}}_{C_{r}})

式中: ${U^{R}}_{C_{0}}$ 和 ${U^{R}}_{C_{r}}$ 分别为衰减介质中的反传一次波和绕射波(本文上标S和R分别表示震源正向传播和检波器反向传播,下标C和Q分别表示Q补偿和Q衰减); $F_{C_{0}}$ 、 $F_{C_{r}}$ 分别为用来产生波场扰动δ ${U^{S}}_{C_{0}}$ 、δ ${U^{S}}_{C_{r}}$ 的源矩阵; ${U^{S}}_{C_{0}}$ 、 ${U^{S}}_{C_{r}}$ 为相应的正向传播波场;δv为v的扰动; $w_{R_{0}}$ 、 $w_{R_{1}}$ 为加权系数; ${U^{R}}_{C} ({U^{R}}_{C_{0}}, {U^{R}}_{C_{r}})$ 和F_C $(F_{C_{0}}, F_{C_{r}})$ 将在下一部分中给出。

1.2 粘声反偏移与伴随拟微分方程

Qu等人用下式描述粘声拟微分方程:

(6)

\begin{array}{l} [\frac{1}{v^{2}} \frac{\partial^{2}}{\partial t^{2}} - \frac{τ ε}{2} \frac{1}{v} \frac{\partial}{\partial t} \sqrt{- \nabla} - \nabla + σ \frac{τ (ε + 1)}{4} \frac{\partial}{\partial t} \nabla] \cdot \\ {U^{S}}_{0} = F \end{array}

式中: ${U^{S}}_{0}$ 为正向传播波场;F为源矩阵;Ñ= $\frac{\partial^{2}}{\partial x^{2}}$ + $\frac{\partial^{2}}{\partial z^{2}}$ ,为Laplacian算子;τ= $\frac{τ_{ε}}{τ_{σ}}$ -1,τ_ε和τ_σ分别为应变和应力松弛时间,并能够通过τ-Q关系τ_σ= $\frac{\sqrt{Q^{2} + 1} - 1}{ω}$ )和τ_ε= $\frac{\sqrt{Q^{2} + 1} + 1}{ω Q}$ 计算得到,Q为质量因数,ω为角频率;t为时间;x和z是空间坐标;σ为正则化参数,本文设为0.01;v为速度,其与参考速度v₀的关系为v=v₀ ${|ω / ω_{0}|}^{(π Q) - 1}$ ,其中ω₀为参考角频率,本文设为10 Hz;ε是控制Q衰减和Q补偿的参数,可以通过下式设定:

\{\begin{array}{l} ε = - 1, {U^{S}}_{0} = {U^{S}}_{Q_{0}}, Q 衰 减; \\ ε = 1, {U^{S}}_{0} = {U^{S}}_{C_{0}}, Q 补 偿 。 \end{array}

在Q补偿的情况下,式(6)可简化为

(7)L^S

{U^{S}}_{0}

(x,t)=F,

式中:x为坐标;L^S为波场正向传播算子,

(8)L^S=

\{\begin{array}{l} {L^{S}}_{Q} = \frac{1}{v^{2}} \frac{\partial^{2}}{\partial t^{2}} + \frac{τ}{2} \frac{1}{v} \frac{\partial}{\partial t} \sqrt{- \nabla} - \nabla, \\ {L^{S}}_{C} = \frac{1}{v^{2}} \frac{\partial^{2}}{\partial t^{2}} + \frac{τ}{2} \frac{1}{v} \frac{\partial}{\partial t} \sqrt{- ▽} - \nabla + σ \frac{τ}{2} \frac{\partial}{\partial t} \nabla 。 \end{array}

根据Born近似理论,可以得到反偏移方程(推导过程见附录B):

(9)L^Sδ

{U^{S}}_{0}

(x,t)=F₀(I) 。

根据伴随状态理论,可以得到伴随方程:

(10)

{L^{R}}_{C} {U^{R}}_{C_{0}}

(x,t)=

d_{Q_{0}}^{c a l}

d_{Q_{0}}^{o b s}

。

式中: ${L^{R}}_{C}$ 为Q补偿波场反向传播算子; $d_{Q_{0}}^{c a l}$ 为一次波衰减合成数据,可用 $d_{Q_{0}}^{c a l}$ =δ ${U^{S}}_{Q_{0}}$ (x_r,t)计算,其中x_r为检波器的坐标。

式(6)的伴随方程是它们的正向传播算子的自动伴随^[25]:

(11)

\begin{array}{l} [\frac{1}{v^{2}} \frac{\partial^{2}}{\partial t^{2}} - \frac{τ ε}{2} \frac{1}{v} \frac{\partial}{\partial t} \sqrt{- \nabla} - \nabla + \\ σ \cdot \frac{τ (ε + 1)}{4} \frac{\partial}{\partial t} \nabla] {U^{R}}_{0} = d_{Q_{0}}^{c a l} - d_{Q_{0}}^{o b s} \end{array}

。

1.3 J-QLSRTM的正向和反向传播波场

附录B的反偏移方程(式(B-5))中的F'_C为

(12)

F_{C}^{’}

\frac{2 δ v}{v^{3}} (\frac{\partial^{2} {U^{S}}_{C}}{\partial t^{2}} - \frac{τ v}{4} \sqrt{- \nabla^{2}} \frac{\partial {U^{S}}_{C}}{\partial t})

则式(5)可以变换为:

(13)

\begin{array}{l} g = w_{R_{0}} [\frac{2}{v^{3}} (\frac{\partial^{2} {U^{S}}_{C_{0}}}{\partial t^{2}} - \frac{τ v}{4} \sqrt{- \nabla^{2}} \frac{\partial {U^{S}}_{C_{0}}}{\partial t}) \cdot {U^{R}}_{C_{0}}] + \\ w_{R_{1}} [\frac{2}{v^{3}} (\frac{\partial^{2} {U^{S}}_{C r}}{\partial t^{2}} - \frac{τ v}{4} \sqrt{- \nabla^{2}} \frac{\partial {U^{S}}_{C_{r}}}{\partial t}) \cdot {U^{R}}_{C_{r}}] 。 \end{array}

式中:正向传播波场 ${U^{S}}_{C_{0}}$ (x,t)可以由

(14)

{L^{S}}_{C} {U^{S}}_{C_{0}}

(x,t)=F(x_s,t)

生成,其中x_s为震源坐标; ${U^{R}}_{C_{r}}$ (x,t)可以由

(15)

{L^{R}}_{C} {U^{R}}_{C_{r}}

(x,t)=

d_{Q_{r}}^{c a l}

d_{Q_{r}}^{o b s}

生成,其中 $d_{Q_{r}}^{c a l}$ 和 $d_{Q_{r}}^{o b s}$ 可以通过使用倾角滤波器或平面波解构滤波器抑制 $d_{Q_{0}}^{c a l}$ 和 $d_{Q_{0}}^{o b s}$ 得到^[29]。

采用最速下降法的k次迭代的步长(β⁽^k⁾)为

(16)β⁽^ik⁾=

\frac{g^{(k)} g^{(k)}}{g^{(k - 1)} g^{(k - 1)}}

用g⁽^k⁾和β⁽^k⁾计算共轭梯度方向z⁽^k⁾和相应的步长α⁽^k⁾。第k次迭代的成像结果可以用下式更新:

(17)I⁽^k⁾=I⁽^k^-1)-α⁽^k⁾z⁽^k⁾ 。

J-QLSRTM中波场延拓算子、反偏移算子、伴随算子的离散格式见附录C。

本方法采用的是复频移完美匹配层边界条件^[30],算子的稳定性条件为^[31]

(18)

\frac{Δ t}{Δ x}

v_p<

\frac{6}{7 \sqrt{2}}

。

2 数值实例

2.1 含小孔洞的3层衰减模型的绕射波计算

第1种模型为3层衰减模型(图1),模型尺寸为2 400 m×2 400 m,水平方向和垂直方向网格间距分别为8 m×8 m。图1给出了不同层位和孔洞的速度和Q值。在第2层反射面上有6个8 m×8 m的孔洞,速度和Q值分别为4 000 m/s和20。

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 三层衰减模型与小尺度孔洞

Fig. 1 Three layer attenuation model with small-scale caves

数据由沿地面排布的301个检波器生成,相邻检波器间隔8 m,中间放炮,震源函数为主频30 Hz的Ricker子波,采样间隔为0.6 ms,记录时间为1.8 s。图2a和图2b比较了黏声介质和声波介质中的直接波、反射波和绕射波的炮记录,图2b中绕射波由于Q衰减明显变弱;图2c和2d分别显示了从衰减炮记录中分离出来的绕射波炮记录和反射波炮记录(切除了直达波),图2c、2d与图2a对比可以看出,绕射波和反射波被准确地分离开。

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 黏声介质(a)和声波介质(b)的炮记录以及从衰减炮记录中分离出的绕射波的炮记录(c)和分离出的反射波(d)

Fig.2 The shot records from the viscoelastic medium(a), acoustic medium(b),and the shot records of diffracted waves separated from the attenuating shot records(c) and separated reflected waves by cutting off the direct wave(d)

图3给出了600 ms波场补偿试验结果。与绕射波的反向波场(图3a和3b)相比,绕射能量在全波形波场中是不可见的(图3c和3d),并且Q衰减造成的能量损失得到了很好的补偿。

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 600 ms波场补偿试验

Fig.3 Wave field compensation test with 600 ms

2.2 衰减Sigsbee2B模型

采用本文提出的J-QLSRTM对Sigsbee2B数据集进行演示。一些小尺度孔洞很难精确成像。该示例中对标准Sigsbee2B模型进行了修改,通过重新采样,x的范围变为12 000 m,z的范围变为4 000 m,x和z方向的空间网格间隔分别为15 m×10 m。衰减Sigsbee2B的速度模型和Q模型分别如图4a和图4b所示。

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 衰减Sigsbee2B速度模型(a)和Q模型(b)

Fig. 4 The attenuating Sigsbee2B velocity model(a)and Q model(b)

衰减数据集由以下观测系统和正演模拟参数合成:衰减数据由201个检波器生成,相邻检波器间隔15 m,中间激发两端接收,共120炮,炮间距75 m,源函数为主频25 Hz的Ricker子波,时间步长为0.6 ms,记录时间为6 s。图5a所示为炮记录,图5b和5c所示的分别为分离出的绕射波和一次波炮记录。在本例中,J-QLSRTM迭代30次。图6a为Q-RTM成像结果。相比之下,如图6b所示的常规无补偿RTM成像结果经过衰减补偿后,分辨率更好,深部能量更强。

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 炮记录(a)、分离的绕射波炮记录(b)和分离的一次波炮记录(c)

Fig. 5 Shot record(a),separated diffractive wave shot records(b) and separated primary wave shot records(c)

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 QRTM (a)和无补偿RTM (b)进行拉普拉斯滤波后的Sigsbee2B模型图像

Fig. 6 Images of the Sigsbee2B model after Laplacian filtering from QRTM(a) and noncompensated RTM(b)

图7a为30次迭代后的成像结果。为了比较效果,还给出了传统的一次波QLSRTM的图像(图7b)和棱柱波QLSRTM结果(图7c)。与一次波QLSRTM图像相比,图7a所示的一次波与绕射波联合成像结果具有更清晰的小尺度孔洞(椭圆圈示)。此外,J-QLSRTM图像具有更好的分辨率、更宽的照度和更平衡的振幅,与棱柱波QLSRTM成像相比,虽然棱柱波QLSRTM结果小尺度孔洞和断层得到了很好的成像,但无法得到大尺度构造的成像结果。采用本文提出的J-QLSRTM方法,平均一次迭代花费的时间约等于常规LSRTM的3.45倍,约等于QLSRTM的1.52倍。

图7

新窗口打开| 下载原图ZIP| 生成PPT

图7 J-QLSRTM(a)、常规的一次波QLSRTM(b)和棱柱波QLSRTM(c)迭代30次后的成像结果

Fig.7 The imaging results after 30 iterations using the proposed J-QLSRTM(a) and the conventional QLSRTM(b) and prism wave(c)

图8为反射系数和使用声波数据的J-LSRTM声波成像,图9为4.5 km处的波数谱。从图9可以看出,J-QLSRTM补偿了幅值损失并对相位进行了校正,成像结果更接近于参考值。

图8

新窗口打开| 下载原图ZIP| 生成PPT

图8 反射系数(a)和使用声波数据的声波LSRTM图像(b)

Fig.8 The reflectivity(a) and acoustic LSRTM image using acoustic data(b)

图9

新窗口打开| 下载原图ZIP| 生成PPT

图9 4 500 m处的波数谱

Fig. 9 Wavenumber spectra at the distance of 4 500 m

实际工作中,速度和Q场通常难以准确获得,因此,本文对速度场和Q场存在的误差情况进行了测试。

图10a所示为存在明显误差的速度场,为真实速度场加入10%的随机误差再经过高斯平滑50次后产生的速度场;采用该速度场得到的J-QLSRTM成像结果如图11a所示。可以看出,当速度场存在明显误差时,J-QLSRTM成像结果较差,存在明显的成像噪声,且同相轴不连续,成像位置不准确。

图10

新窗口打开| 下载原图ZIP| 生成PPT

图10 存在明显误差的速度场(a)和Q场(b)

Fig.10 Velocity field(a) and Q field(b) with obvious errors

图11

新窗口打开| 下载原图ZIP| 生成PPT

图11 采用图10a速度场(a)和图10bQ场(b)得到的J-QLSRTM成像结果

Fig.11 J-QLSRTM imaging results obtained by using the velocity field in Fig.10a(a) and by using the Q field in Fig.10b(b)

图10b所示为存在明显误差的Q场,该模型为真实Q模型增大20%,并采用高斯平滑100次得到的Q模型;采用该Q模型得到的J-QLSRTM成像结果如图11b所示。从图中可以看出,虽然输入的Q极不准确,但依旧得到了可接受的成像结果。因此,可以说明,该J-QLSRTM对速度模型比Q模型的敏感性更高。

3 结论

本文通过联合应用一次波和绕射波构建新的目标函数和梯度公式,研发了一种Q补偿的一次波和绕射波联合LSRTM方法,以提高小尺度断层和孔洞的成像效果。基于反演理论提出的方法补偿了沿一次波和绕射波传播路径的Q衰减。数值实例表明,与传统的QLSRTM和未补偿的J-QLSRTM相比,该方法通过对一次波和绕射波的联合成像和补偿衰减,获得了具有更高质量、更清晰的复杂结构信息、更强深部能量和更高分辨率的图像。

附录A:J-QLSRTM梯度方程的推导

目标函数δE的扰动可以表示为

(A-1)

$\begin{array}{c}\delta E=\frac{1}{2} \delta\left\|_{a_{0}}\left(\boldsymbol{L}_{Q_{0}}^{\mathrm{S}} \boldsymbol{m}-\boldsymbol{d}_{Q_{0}}^{\mathrm{obs}}\right)+a_{1}\left(\boldsymbol{L}_{Q_{\mathrm{r}}}^{\mathrm{S}} \boldsymbol{m}-\boldsymbol{d}_{Q_{\mathrm{r}}}^{\mathrm{obs}}\right)\right\|_{2}^{2}= \\\left\langle\left(a_{0} \boldsymbol{Y} \boldsymbol{U}_{Q_{0}}^{\mathrm{S}}+a_{1} \boldsymbol{Y} \boldsymbol{U}_{Q_{\mathrm{r}}}^{\mathrm{S}}\right),\left[a_{0}\left(\boldsymbol{Y} \boldsymbol{U}_{Q_{0}}^{\mathrm{S}}-\boldsymbol{d}_{Q_{0}}^{\mathrm{obs}}\right)+a_{1}\left(\boldsymbol{Y} \boldsymbol{U}_{Q_{\mathrm{r}}}^{\mathrm{S}}-\boldsymbol{d}_{Q_{\mathrm{r}}}^{\mathrm{obs}}\right)\right]\right\rangle= \\\left\langle\left[a_{0} \boldsymbol{Y}\left(\boldsymbol{L}_{C_{0}}^{\mathrm{S}} \boldsymbol{F}_{C_{0}}\right)+a_{1} \boldsymbol{Y}\left(\boldsymbol{L}_{C_{\mathrm{r}}}^{\mathrm{S}} \boldsymbol{F}_{C_{\mathrm{r}}}\right)\right],\left[a_{0}\left(\boldsymbol{Y} \boldsymbol{U}_{Q_{0}}^{\mathrm{S}}-\boldsymbol{d}_{Q_{0}}^{\mathrm{obs}}\right)+\right.\right. \\\left.\left.a_{1}\left(\boldsymbol{Y} \boldsymbol{U}_{Q_{\mathbf{r}}}^{\mathrm{s}}-\boldsymbol{d}_{Q_{\mathrm{r}}}^{\mathrm{obs}}\right)\right]\right\rangle=\left\langle a_{0} \boldsymbol{Y}\left(\boldsymbol{L}_{C_{0}}^{\mathrm{S}} \boldsymbol{F}_{C_{0}}\right), a_{0}\left(\boldsymbol{Y} \boldsymbol{U}_{Q_{0}}^{\mathrm{S}}-\boldsymbol{d}_{Q_{0}}^{\mathrm{obs}}\right)\right\rangle+ \\\left\langle a_{1} \boldsymbol{Y}\left(\boldsymbol{L}_{C_{\mathrm{r}}}^{\mathrm{S}} \boldsymbol{F}_{C_{\mathrm{r}}}\right), a_{1}\left(\boldsymbol{Y} \boldsymbol{U}_{Q_{\mathrm{r}}}^{\mathrm{S}}-\boldsymbol{d}_{Q_{\mathrm{r}}}^{\mathrm{obs}}\right)\right\rangle+\text { crosstalks, }\end{array}$

式中:Y为将波场限制在检波器位置处的算子。忽略串扰项,式(A-1)可以简化为

(A-2)$\begin{array}{c}\delta E=a_{0}^{2}\left\langle\boldsymbol{Y}\left(\boldsymbol{L}_{C_{0}}^{\mathrm{s}} \boldsymbol{F}_{C_{0}}\right),\left(\boldsymbol{Y} \boldsymbol{U}_{Q_{0}}^{\mathrm{s}}-\boldsymbol{d}_{Q_{0}}^{\mathrm{obs}}\right)\right\rangle+a_{1}^{2}\left\langle\boldsymbol{Y}\left(\boldsymbol{L}_{C_{\mathrm{r}}}^{\mathrm{s}} \boldsymbol{F}_{C_{\mathrm{r}}}\right),\right. \\\left.\left(\boldsymbol{Y} \boldsymbol{U}_{Q_{\mathrm{r}}}^{\mathrm{s}}-\boldsymbol{d}_{Q_{\mathrm{r}}}^{\mathrm{obs}}\right)\right\rangle=a_{0}^{2}\left\langle\boldsymbol{F}_{C_{0}}, \boldsymbol{L}_{C_{0}}^{\mathrm{R}} \boldsymbol{Y}^{\mathrm{R}}\left(\boldsymbol{Y} \boldsymbol{U}_{Q_{0}}^{\mathrm{s}}-\boldsymbol{d}_{Q_{0}}^{\mathrm{obs}}\right)\right\rangle+ \\a_{1}^{2}\left\langle\boldsymbol{F}_{C_{\mathrm{r}}}, \boldsymbol{L}_{C_{\mathrm{r}}^{\mathrm{R}}}^{\mathrm{R}} \boldsymbol{Y}^{\mathrm{R}}\left(\boldsymbol{Y} \boldsymbol{U}_{Q_{\mathrm{r}}}^{\mathrm{s}}-\boldsymbol{d}_{Q_{\mathrm{r}}}^{\mathrm{obs}}\right)\right\rangle 。\end{array}$

式中:Y^R为将接收器数据扩展到整个模型空间的算子,并且

(A3)

{L^{R}}_{C_{r}}

Y^R(Y

{U^{S}}_{Q_{r}}

d_{Q_{r}}^{o b s}

{U^{R}}_{C_{r}}

。

令 $w_{R_{0}}$ = $a_{0}^{2}$ , $w_{R_{1}}$ = $a_{1}^{2}$ ,则有

(A4)g=

\frac{δ E}{δ v}

\frac{1}{δ v} (w_{R_{0}} F_{C_{0}} {U^{R}}_{C_{0}} + w_{R_{1}} F_{C_{r}} {U^{R}}_{C_{r}})

。

附录B:基于粘声拟微分方程的反偏移算子

基于Born近似,参数扰动δv=v-v₀导致扰动波场δU=U-U₀,U和U₀由下式确定:

(B-1)

[\frac{1}{(v_{0} {+ δ v)}^{2}} \frac{\partial^{2}}{\partial t^{2}} - \frac{τ ε}{2} \frac{1}{(v_{0} + δ v)} \frac{\partial}{\partial t} \sqrt{- \nabla} - \nabla + σ \frac{τ (ε + 1)}{4} \frac{\partial}{\partial t} \nabla]

(U₀+δU)=0,

(B-2)

[\frac{1}{v_{0}^{2}} \frac{\partial^{2}}{\partial t^{2}} - \frac{τ ε}{2} \frac{1}{v_{0}} \frac{\partial}{\partial t} \sqrt{- \nabla} - \nabla + σ \frac{τ (ε + 1)}{4} \frac{\partial}{\partial t} \nabla]

U₀=0 。

基于泰勒展开,近似有

(B-3)

\{\begin{array}{l} \frac{1}{v^{2}} = \frac{1}{v_{0}^{2}} - \frac{2}{v_{0}^{3}} + O (δ v), \\ \frac{1}{v} = \frac{1}{v^{0}} - \frac{δ v}{v_{0}^{2}} + O (δ v) 。 \end{array}

式中:O(δv)对应于v的高次项。

式(B-2)减去式(B-1),并忽略δv的高次项,得到

(B-4)

[\frac{1}{v_{0}^{2}} \frac{\partial^{2}}{\partial t^{2}} - \frac{τ ε}{2} \frac{1}{v_{0}} \frac{\partial}{\partial t} \sqrt{- \nabla} - \nabla + σ \frac{τ (ε + 1)}{4} \frac{\partial}{\partial t} \nabla]

·
δU=

(\frac{2 δ v}{v_{0}^{3}} \frac{\partial^{2}}{\partial t^{2}} - \frac{τ ε}{2} \frac{δ v}{v_{0}^{2}} \frac{\partial}{\partial t} \sqrt{- \nabla})

U₀ 。

因此,可以构建源矩阵 $F_{C_{0}}$ 、 $F_{C_{r}}$ 以生成波场扰动δ ${U^{S}}_{C_{0}}$ 、δ ${U^{S}}_{C_{r}}$ :

(B-5)

F_{C}^{'}

\frac{2 δ v}{v^{^{3}}} (\frac{\partial^{2} {U^{^{S}}}_{C}}{\partial t^{^{2}}} - \frac{τ v}{4} \sqrt{- \nabla^{^{2}}} \frac{\partial {U^{^{S}}}_{C}}{\partial t})

。

将反射系数模型定义为

(B-6)I(x)=

\frac{2 δ v}{v_{0}}

将式(B-6)代入式(B-4),得到粘声反偏移方程:

(B-7)

[\frac{1}{v_{0}^{2}} \frac{\partial^{2}}{\partial t^{2}} - \frac{τ ε}{2} \frac{1}{v_{0}} \frac{\partial}{\partial t} \sqrt{- \nabla} - \nabla + σ \frac{τ (ε + 1)}{4} \frac{\partial}{\partial t} \nabla]

δU=I(x)

(\frac{2}{v_{0}^{2}} \frac{\partial^{2}}{\partial t^{2}} - \frac{τ ε}{4} \frac{1}{v_{0}} \frac{\partial}{\partial t} \sqrt{- \nabla})

U₀ 。

附录C:算子的离散格式

本文采用混合空间差分格式进行数值离散,即衰减项在时间-波束域进行离散,其他项在时间-空间域进行离散。

(C-1)

{\frac{\partial^{2} U}{\partial t^{2}}|}_{t = n}

≈

\frac{1}{Δ t^{2}}

(Uⁿ⁺¹+Uⁿ^-1-2Uⁿ),

式中:Δt为时间采样间隔,n表示时间离散点数。

衰减项在时间-波束域的离散格式为:

(C-2)

{\sqrt{- (\frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial z^{2}})} \frac{\partial U}{\partial t}|}_{t = n}

≈

\frac{1}{Δ t}

Φ^-1

\sqrt{k_{x}^{2} + k_{z}^{2}}

·[Φ(Uⁿ-F(Uⁿ^-1)],

式中:Φ表示傅里叶变换,Φ^-1表示傅里叶逆变换,Δt为时间采样间隔,k_x为水平方向的波束,k_z为垂直方向的波束。

其他项在空间域的离散格式为

(C-3)

{\frac{\partial^{2} U}{\partial x^{2}}|}_{x = i}

{D^{2}}_{x}

U≈

\frac{1}{Δ x^{2}}

[c₀U_i+

\overset{M}{\sum_{m = 1}}

c_m(U_i₊_m+U_i_-_m)],

(C-4)

{\frac{\partial^{2} U}{\partial z^{2}}|}_{z = j}

D_{z}^{2}

U≈

\frac{1}{Δ z^{2}}

[c₀U_j+

\overset{M}{\sum_{m = 1}}

c_m(U_j₊_m+U_j_-_m)],

式中:M为差分精度(2M表示差分阶数),c₀、c_m表示差分稀疏,i和j表示空间水平分量和垂直分量的离散点数,m为求和计数单位。

反偏移算子的递推公式为

(C-5)

U_{0}^{S, n + 1}

U_{0}^{S, n}

U_{0}^{S, n - 1}

\frac{τ ε v}{2}

ΔtΦ^-1

\sqrt{k_{x}^{2} + k_{z}^{2}}

·[Φ(

U_{0}^{S, n}

)-Φ(

U_{0}^{S, n - 1}

)]+Δt²v²

{D^{2}}_{x} U_{0}^{S, n}

+Δt²v²

{D^{2}}_{z} U_{0}^{S, n}

+F(t)·v²,

(C-6)$\begin{array}{c}\delta U^{n+1}=2 \delta U^{n}-\delta U^{n-1}+\frac{\tau \varepsilon v}{2} \Delta t \Phi^{-1} \sqrt{k_{x}^{2}+k_{z}^{2}} \cdot \\{\left[\Phi\left(\delta U^{n}\right)-\Phi\left(\delta U^{n-1}\right)\right]-\frac{\sigma \tau(\varepsilon+1) v^{2}}{4} \Delta t \cdot} \\\left(\delta U^{n}-\delta U^{n-1}\right)+\Delta t^{2} v^{2} D_{x}^{2} \delta U^{n}+\Delta t^{2} v^{2} D_{z}^{2} \delta U^{n}+I(x) \cdot \\\left\{2\left(U_{0}^{S, n+1}-2 U_{0}^{S, n}+U_{0}^{S, n-1}\right)-\frac{\tau \varepsilon v}{4} \Delta t \Phi^{-1} \sqrt{k_{x}^{2}+k_{z}^{2}} \cdot\right. \\\left.\left[\Phi\left(U_{0}^{S, n}\right)-\Phi\left(U_{0}^{S, n-1}\right)\right]\right\}\end{array}$ 。

伴随算子的递推公式为

(C-7)

U_{0}^{R, n - 1}

U_{0}^{R, n}

U_{0}^{R, n + 1}

\frac{τ ε v}{2}

ΔtΦ^-1

\sqrt{k_{x}^{2} + k_{z}^{2}}

·[Φ(

U_{0}^{R, n}

)-Φ(

U_{0}^{R, n - 1}

)]+Δt²v²

{D^{2}}_{x} U_{0}^{R, n}

+Δt²v²

{D^{2}}_{z} U_{0}^{R, n}

(C-8)δ

U^{^{n + 1}}

=2δ

U^{^{n}}

-δUⁿ^-1+

\frac{τ ε v}{2}

ΔtΦ^-1

\sqrt{k_{x}^{2} + k_{z}^{2}}

·[Φ(δ

U^{^{n}}

)-Φ(δ

U^{^{n - 1}}

)]-

\frac{σ τ (ε + 1) v^{2}}{4}

Δt·(

U^{_{0}^{R, n}}

U^{_{0}^{R, n - 1}}

d_{Q_{0}}^{c a l}

d_{Q_{0}}^{o b s}

。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

Schuster

, Wang

, Huang

, et al.

Theory of multisource crosstalk reduction by phase-encoded statics

[J]. Geophysical Journal International, 2011, 184(3):1289-1303.

DOI:10.1111/j.1365-246X.2010.04906.x URL [本文引用: 1]

[2]

Kozlov

, Barasky

, Korolev

Imaging scattering objects masked by specular reflections

[J]. Society of Exploration Geophsicists, 2004, 23(1):1131-1134.