融合点扩散函数的预条件黏声最小二乘逆时偏移

姚振岸 孙成禹 喻志超 马 振

(①东华理工大学放射性地质与勘探技术国防重点学科实验室，江西南昌 330013；②中国石油大学(华东)地球科学与技术学院，山东青岛 266580；③北京大学地球与空间科学学院石油与天然气研究中心，北京 100871)

0 引言

储层孔隙中含有油气水等流体，会引起地震波传播过程中的强振幅衰减和相位畸变[1]，从而导致传统的声波或者弹性波偏移方法分辨率降低或者失效。关于地震波在黏弹性介质中的传播机理，主要有Maxwell模型、Kelvin-Vogit模型、标准线性固体、广义标准线性体、分数阶常Q模型等[2-6]。早期的地震资料衰减补偿主要是使用反Q滤波在数据域实现[7-10]，如Blanch等[11-12]提出用τ方法模拟常Q，并基于伴随状态法在τ-p域中用线性反演方法恢复黏声介质的短波长分量。考虑到地震波衰减是在传播过程中发生的，融合衰减补偿的偏移方法越来越受到重视。Xin等[13]和Xie等[14]利用射线追踪的方法在叠前深度偏移过程中进行衰减补偿；白敏等[15]和代福材等[16]在高斯束叠前深度偏移过程中补偿振幅衰减和校正相位畸变；Dai等[17]、Yu等[18]、Wang[19]和Valenciano等[20]在频率域用单程波动方程偏移进行衰减补偿；陈志德等[21]引入等效Q值实现频率域波场延拓，基于单程波方程和稳定相点原理实现黏声叠前时间偏移，并提高了补偿算法的稳定性。朱峰等[22]推导了黏声VTI介质高阶傅里叶有限差分波场延拓算子，实现了适应高角度入射波场的黏声叠前深度偏移方法。

不同学者提出了多种黏声波动方程，在逆时偏移过程中分别对振幅损失和相位畸变进行补偿和校正[23-27]。李振春等[28]基于广义标准线性体模型对黏声介质波动方程进行线性化，并借助伴随状态法实现了黏声介质最小二乘逆时偏移(Q-LSRTM)的迭代算法，获得了较高的成像分辨率。Dutta等[29]、Dai等[30]基于一阶松弛机制时间域黏声波动方程，用迭代的Q-LSRTM有效补偿和校正了强衰减层导致的振幅损失和相位畸变。徐凯等[31]在反演框架下，构建了Q-LSRTM多约束正则化目标函数，并采用变分方法求解，消除了震源波场与接收点波场互相关引起的噪声，提高了复杂地质体成像的精度、增强了保幅性。相比于传统声波LSRTM，Q-LSRTM能实现反射体的精确归位且振幅更均衡，收敛速率更快[32]。但Q-LSRTM也存在一定问题： 其一，由于在传统Q-LSRTM中用于反传数据残差的伴随Q传播算子也是衰减的，导致成像结果分辨率降低；其二，基于时间域黏声波动的LSRTM计算成本较高，尤其是需要很多次迭代才能达到预期成像效果。在完全弹性情况下，为了提高传统LSRTM的收敛速度，李闯等[33]发展了一种预条件偏移算法，提高了盐下构造成像的分辨率和保幅性。

在实际野外地震勘探中，由于稀疏的震源和检波器，记录偏移孔径有限、速度变化快及构造复杂，以及不均匀的地震照明等原因，偏移剖面当中往往存在多种假象[34]。为了消除偏移假象、提升成像分辨率，偏移反褶积被广泛采用[35]。偏移成像本质上就是真实反射率的一种模糊，模糊算子即为Hessian算子，而偏移反褶积处理就等效于一种近似的Hessian逆处理。Hu等[34-35]在空间—波数域设计了一个偏移反褶积算子压制叠后偏移假象；Yu等[36]将偏移反褶积从叠后偏移拓展到叠前深度偏移。Guitton[37]直接在空间域构建了一系列匹配滤波器近似Hessian矩阵的逆实现偏移反褶积，后来这些匹配滤波器被Aoki等[38]作为预条件化因子用于传统的LSRTM，获得了良好的成像效果，这些匹配滤波器也被称为去模糊滤波器。针对多震源LSRTM，Dai等[39-40]用去模糊滤波器减少串扰噪声并且加速多震源LSRTM的收敛。但目前偏移反褶积的原理和实施都是基于无耗散介质假设，而地下储层往往具有衰减性质。

为了克服Q-LSRTM的缺陷，本文将匹配滤波偏移反褶积与Q-LSRTM相结合。首先基于点扩散函数(PSF)构建黏声去模糊滤波器，然后将其作为Q-LSRTM迭代过程中的预条件化因子，最终实现高分辨的衰减补偿偏移。相比于Q-LSRTM，预条件Q-LSRTM成像结果分辨率更高，反射体振幅更均衡，收敛速度更快。就计算成本来说，尽管预条件Q-LSRTM需要预先计算黏声去模糊滤波器，但其达到预期成像效果需要的迭代次数远比Q-LSRTM少。

1 理论方法

1.1 点扩散函数与去模糊滤波器

在声学介质中，基于Born近似，观测地震记录d的形成过程可以表示为

d=Lm0

(1)

式中：L表示线性模拟算子，与观测系统、震源子波和地下介质模型参数有关；m0是地下反射率模型。对应L，可以获得其伴随偏移算子LT，则偏移结果mmig可表示为

(2)

mmig其实是真实反射率模型m0的模糊版本，LTL就是一个模糊算子，即偏移格林函数[41]，也被称为点扩散函数(PSF)。由于震源子波、观测系统等多种因素的影响，偏移成像中往往存在较多的假象。为了获得更清晰的成像结果，可以取LTL的逆并作用于成像结果

(3)

但在地震成像时，直接计算(LTL)-1是不现实的，为此采用迭代LSRTM方法予以求解，但其计算量往往比标准偏移成像大一个数量级。为了加速收敛，设计去模糊算子的近似逆算子(LTL)-1，并在LSRTM中作为预条件化算子[35-39]。

实际应用中，去模糊滤波器基于背景模型即偏移速度场和参考反射率模型获得，以均匀介质速度背景场和均匀分布点散射体反射率模型为例阐述去模糊滤波器的获取过程。图1a为均匀分布点散射体模型，图1b为对应的叠前深度偏移成像结果，其中震源和检波器均匀分布于模型上方。单个散射体的PSF响应即为其偏移结果。记参考反射率模型为mref，则在Bron近似下利用式(1)可获得正演模拟数据dref，继而可以获得标准的参考偏移结果

mmig-ref=LTLmref=LTdref

(4)

将点散射体模型分解成多个子区域(如图1a中黑色矩形框)，使点散射体位于矩形窗口的中央，窗口大小为wx×wz，足以覆盖其PSF响应的有效能量(图1b中矩形框)。在这个局部窗口内，假设PSF保持不变。在wx×wz局部窗口区域内，去模糊滤波器将参考偏移成像结果和参考反射率模型联系起来[38，40]，即

mmig-ref(x0,z0)idV

(5)

式中：i为局部窗口序号；积分限Vi为整个局部窗口区域。式(5)写成矩阵形式为

[mref]i=Fi⊗[mmig-ref]i

(6)

式中：“⊗”表示空间卷积；[mref]i和[mmig-ref]i分别为第i个局部窗口内参考反射率模型和参考偏移成像结果。与PSF相同，在每个局部窗口内去模糊滤波器保持不变，定义滤波器的空间尺寸为fx×fz。

图1 点散射模型及其偏移成像转换示意图

由卷积的交换律可得

Fi⊗[mmig-ref]i=[mmig-ref]i⊗Fi

(7)

将上述卷积运算写成矩阵运算形式，参考偏移结果[mmin-ref]i变成一个(N+M-1)×N阶卷积矩阵，去模糊滤波器Fi变为N×1阶向量f，其中N=fx×fz且M=wx×wz，即为

[mmig-ref]ifi=[mref]i

(8)

(9)

继而采用三角(LU)分解方法求解式(9)即可获得去模糊滤波器，在第i个局部窗口内近似等于(LTL)-1。在实际应用中，为了保证计算效果，建议局部窗口长度不超过1.5倍地震波长。

1.2 黏声介质中的去模糊滤波器

将去模糊滤波器应用到偏移成像的过程就是一种偏移反褶积处理。前人对偏移反褶积的研究都集中在无耗散介质当中，但强吸收会对地震波振幅和相位产生明显的衰减和畸变[1]。为了补偿衰减效应，Dutta等[29]发展了时间域Q-LSRTM技术，成像结果归位更精确，振幅更均衡。在传统Q-LSRTM中，用于反传数据残差的伴随Q传播算子也是衰减的，因此成像分辨率有一定的损失。

这种分辨率的损失能通过黏声介质中的PSF即偏移格林函数予以解释。给定速度为v0的均匀介质和位于xs处角频率为ω的简谐点源，声学格林函数可以表示为

(10)

如果介质是耗散的，则将完全弹性介质中的实速度v0替换成复速度[1]

(11)

从而得到黏声介质格林函数

(12)

(13)

(14)

(15)

如1.1所述，在偏移成像的不同子区域内，通过匹配参考成像数据和参考反射率模型，得到局部去模糊滤波器，即

(16)

1.3 黏声最小二乘逆时偏移

基于标准线性固体(SLS)黏弹性模型一阶松弛机制，二维黏声介质波动程可写为[11，42]

(17)

式中：ρ为密度；P为压力；K为体积模量；S为震源子波；γp为记忆变量；v为质点速度向量；τ与品质因子Q有关，可表示为

(18)

其中τσ和τε分别表示应力和应变松弛时间。在Born近似的框架下，对体积模量K做一定扰动δK，则扰动场可以写为

(19)

式中忽略了二阶扰动变量。若用压力和记忆变量格林函数GP(xr,t;x0,0)和Gγp(xr,t;x0,0)，式(19)也可表示为

δP(xr,t;xs)

(20)

式中“*”表示时间方向的褶积运算。

在Q-LSRTM的框架下，式(20)的解等效于矩阵—向量运算

dQ=LQm0

(21)

式中：dQ为Born模拟的黏声地震数据；LQ为线性黏声模拟算子；m0为地下介质反射率。用伴随状态法可推导式(19)的伴随方程[11，12，43]为

(22)

式中： (q，u，s)为(P，v，γp)的伴随状态变量；Δd表示每次迭代中实际观测数据和正演数据之间的残差。

成像扰动δm与体积模量扰动δK近似线性相关，而δK可由背景场(通过式(19)计算)和伴随场(通过式(22)计算)之间的零时间延迟互相关得到，即

(23)

上式可用格林函数表示为

[GP(xr,-t;x,0)*ΔP(xr,t;xs)]+

[GτP(xr,-t;x,0)*ΔP(xr,t;xs)]}

(24)

用矩阵可以表示为

(25)

1.4 预条件Q最小二乘逆时偏移

Q-LSRTM的误差泛函[29]为

(26)

(27)

(28)

式中α是更新步长。

2 Marmousi Ⅱ模型测试

为了说明本文提出的预条件Q-LSRTM的有效性，应用Marmousi Ⅱ模型进行测试。图2a和图2b分别为正演速度模型和Q值模型，其中Q只取两个值，Q=20表示衰减层，Q=10000可视为弹性介质。整个模型的网格点数为801×351，模型纵横向采样间隔均为10m。在合成黏声地震记录时，炮间隔为50m，共计150炮均匀分布在模型表面，道间距为10m，震源采用主频为15Hz的Ricker子波。采集系统位于地表，且排列长度固定为整个模型的长度。计算中采用时间二阶精度、空间八阶精度的时间域交错网格有限差分法及C-PML边界条件。本文共用到两套数据，一套是基于声波方程的正演结果，另一套是基于一阶松弛机制黏声波动方程的正演结果。对真实速度模型和Q值模型进行二维平滑，得到偏移速度模型(图2c)和偏移Q值模型(图2d)。分别执行声波LSRTM、Q-LSRTM、预条件Q-LSRTM三种最小二乘逆时偏移，迭代次数均为20次，其中第1次迭代结果即为声波逆时偏移(RTM)、Q-RTM和预条件Q-RTM的成像结果。图3为声波数据声波RTM和LSRTM成像结果，并以此作为基准对比数据。

图2 Marmousi Ⅱ速度及Q模型

图3 声波数据声波RTM(a)和LSRTM(b)成像结果

图4 黏声数据不同方法成像结果对比

图5、图6分别为声波数据LSRTM(图3b)、黏声数据LSRTM(图4b)、Q-LSRTM(图4d)、预条件Q-LSRTM(图4f)成像结果蓝色、红色矩形框局部放大显示，可以看出：预条件Q-LSRTM提高了黏声数据的成像精度和分辨率(蓝色箭头所指)；采用LSRTM对黏声数据进行成像时，不仅存在成像分辨率低的问题，而且由于介质黏弹性导致地震波振幅和相位的畸变，还会导致深层反射体归位出现错乱(图6b红色箭头所指)。相对而言，尽管Q-LSRTM能有效归位主要反射界面，但分辨率较低，且不能呈现构造细节，预条件Q-LSRTM成像与基准声波LSRTM结果比较接近，分辨率较高。

图7为声波数据LSRTM(图5a)、黏声数据LSRTM(图5b)、Q-LSRTM(图5c)、预条件Q-LSRTM(图5d)在x=5.0km处成像结果的中深层归一化波数谱对比，可以看出，传统LSRTM和Q-LSRTM的主波数都相对偏低，只有预条件Q-LSRTM的波数谱与基准声波LSRTM接近，也说明预条件Q-LSRTM提高了黏声数据的成像分辨率。

图5 图3、图4蓝色矩形框局部放大显示

图6 图3、图4红色矩形框局部放大显示

图8为声波数据LSRTM(图6a)、黏声数据LSRTM(图6b)、Q-LSRTM(图6c)、预条件Q-LSRTM(图6d)在x=2.6km处成像结果的深层归一化波数谱对比，可以看出，黏声数据的声波LSRTM波数谱严重偏离基准数据谱，Q-LSRTM有所改善，预条件Q-LSRTM能恢复出高波数细节信息，最接近于基准数据。

图7 x=5.0km处不同方法成像结果的中深层波数谱对比

图8 x=2.6km处不同方法成像结果的深层波数谱对比

图9为黏声数据声波LSRTM、Q-LSRTM、预条件Q-LSRTM(预条件20次)、预条件Q-LSRTM(预条件6次)四种情况下数据残差随着迭代次数的变化曲线，可以看出，四种偏移成像结果的精度是逐渐升高的。就收敛速率来看，预条件Q-LSRTM第4次迭代和第8次迭代时的数据残差接近于Q-LSRTM第8次和第20次迭代的数据残差，因此预条件Q-LSRTM能减少约50%的计算量，尤其在早期的几次迭代当中体现的更为明显。预条件Q-LSRTM收敛速率提升的原因在于黏声去模糊滤波器是Hessian算子逆的一种有效近似估计，尽管需要预先计算去模糊滤波器，但计算成本还是要比Q-LSRTM低得多。兼顾到计算精度和计算效率，参考反射率模型中点散射体的布置应该相对稀疏，由此会引入一些全局噪声，如果预条件在第6次迭代后终止，数据残差会进一步减小，相应地成像精度也有所提高。

为了测试预条件Q-LSRTM对Q值的敏感性，分别应用Q=20、40、80、100四种模型进行偏移成像，结果如图4f、图10所示，图11为对应的参考反射率模型的黏声PSF，可以看出：随着Q值的增大，衰减层及其下部的反射体的成像振幅逐渐减弱，成像分辨率逐渐降低；相反地，当Q值越大时，衰减越弱，参考反射率模型中的点散射体偏移响应即PSF越强，从而用其构建黏声去模糊滤波器时，考虑到的衰减效应越弱、成像精度越低。

图9 数据残差随着迭代次数的变化曲线

图10 不同Q值偏移模型的预条件Q-LSRTM偏移结果

因此，一个相对精确的偏移Q值模型是预条件Q-LSRTM取得良好成像效果的前提。

图11 不同Q值偏移模型的参考反射率模型的黏声PSF

3 结论与认识

(1)在传统Q-LSRTM成像当中，数据残差的反传过程也是衰减的，会导致成像分辨率降低；而基于PSF的预条件Q-LSRTM方法采用黏声去模糊滤波器补偿由于地下强衰减导致的振幅和分辨率降低，能得到较高的成像精度和更均衡的成像振幅。

(2)由于地下模型是未知的，构建黏声去模糊滤波器时一般选用均匀分布的点散射体模型。基于偏移速度场和Q值模型，采用和实际一致的观测系统和地震子波，通过正演模拟和偏移成像获得黏声PSF，最终采用匹配滤波的思想构建黏声去模糊滤波器。兼顾到计算精度和计算效率，参考反射率模型中点散射体密度和局部窗口长度需根据构造复杂度确定。

(3)与Q-LSRTM相比，预条件Q-LSRTM的收敛效率得到大幅度提升，在前几次迭代中更为明显。为了获得更高的成像精度和计算效率，预条件Q-LSRTM只在前几次迭代中使用为宜。就计算成本而言，尽管预条件Q-LSRTM需要预先计算黏声去模糊滤波器，但计算耗时仍然比Q-LSRTM小得多。

(4)与Q-LSRTM类似，为了获得较好的成像效果，预条件Q-LSRTM也需要一个相对精确的Q值模型估计，同时也必须有一个相对精确的偏移速度模型。

(5)本文的方法可拓展到黏弹各向异性等更复杂介质的偏移成像当中，考虑到计算量，建议采用并行加速计算。