全局约束反演多道吸收补偿方法

2021-05-15刘金涛王小卫刘梦丽

石油地球物理勘探 2021年2期

刘金涛王孝王小卫苏勤刘梦丽袁焕

(中国石油勘探开发研究院西北分院，甘肃兰州 730030)

0 引言

众所周知，地震波在地下介质中传播时，由于大地滤波作用，会导致地震波能量衰减、速度频散，降低分辨薄层的能力。目前已提出了许多地层吸收衰减模型(Q模型)描述这一现象，其中常Q模型是地震勘探领域使用最多的地层吸收衰减模型[1]。地震波的衰减效应破坏了地震子波在时间方向的稳态特征，导致传统的反褶积方法很难从非稳态地震记录中准确地恢复反射系数序列。因此，基于吸收补偿的提高分辨率方法具有明确的物理意义和较大应用潜力[2-8]。

反Q滤波是吸收补偿的主要方法，有多种实现方式，其中，波场延拓和稀疏反演是两类应用最为广泛的方法。由于波场延拓方法在外推过程中存在放大函数，极易放大高频噪声，结果不稳定。因此，稀疏反演类反Q滤波方法应运而生[9-10]。Zhang等[11]将衰减补偿看作基于反演的最小二乘问题，并用贝叶斯理论进行正则化约束；Braga等[12]在小波域实现了L2范数约束的反Q滤波；Oliveira等[13]在频率域实现了L1范数约束的衰减补偿，Wang等[14]对其进行改进，在时间域实现了L1范数约束的衰减补偿。以上方法都只考虑了地震道纵向上的信息，而没有考虑相邻道之间的相关性，即都是单道补偿方法，因而结果往往会出现横向不连续问题。为了解决这一问题，将相邻道之间的相关性作为先验信息加入反演系统进行多道反演补偿的方法逐渐发展起来。

Zhang等[15]提出了一种多道反演方法，该算法将横向约束融入反演算法中，反演的结果较常规单道算法具有更好的横向连续性。Kazemi等[16]假设反射系数序列稀疏分布，提出了一种非线性多道盲源反褶积方法，主要优点是抗噪性较强且不需要子波的先验信息。Yuan等[17]先将数据转化到变换域，然后在变换域数据上加一个约束项，实现了多道波阻抗反演，在一定程度上缓解了单道反演算法的横向不连续性问题；在此基础上，Yuan等[18]提出了一种在频率—波数域的多道衰减补偿方法，增强了地震同相轴横向连续性，提高了吸收补偿的稳定性和抗噪性。

f-x域预测滤波是描述线性事件空间预测特性的有用工具，主要被用于数据插值和去噪领域[19-22]。为提高地层吸收补偿的稳定性和抗噪性，本文提出了一种基于反演的全局约束吸收补偿方法。该方法首先在f-x域计算信号预测算子[23-25]，用于表征地震信号在空间上的结构关系和能量变化；然后，将该算子作为横向正则化约束引入到纵向吸收补偿系统，建立具有双向(纵向和横向)即全局约束的反演目标函数；最后应用模型数据和实际数据验证了方法的有效性。

1 地层吸收模型

1.1 Q的定义

当平面波通过均匀的黏弹性介质时，其振幅衰减和速度频散可以用一个无量纲的参数Q描述

(1)

式中：ω为圆频率；α(ω)和v(ω)分别为衰减系数和相速度。由于品质因子Q必须为正数，所以

(2)

进一步假设

(3)

则式(1)可以近似为

(4)

1.2 修正的Kolsky-Futterman模型

Kolsky[26]假设在测量频带范围内衰减系数α(ω)与频率呈严格的线性关系，即

(5)

式中ω0为参考圆频率。同时，Kolsky定义其相速度为

(6)

将式(5)和式(6)代入式(4)，可得

(7)

对于较大的品质因子，即Q(ω0)≫1，有

(8)

式中γ=1/[πQ(ω0)]。则相速度和品质因子可近似为

(9)

(10)

其中自由选择的参考频率ω0小于最低的测量频率。

Wang等[27]指出，式(6)相速度仅仅是对于ω≫ω0的一个渐近表达式。由于地震勘探资料的频率范围相对较低(一般小于500Hz)，其修正的相速度公式为

(11)

式中：h是一个频率不相关的常数；ωh为重新定义的调谐参数。

对于修正的相速度模型，其相应的修正Q模型为

(12)

本文所用的地层吸收模型即为上述模型。

2 基于反演的纵向约束吸收补偿方法

在完全弹性介质中，可以用褶积模型简单描述地震数据采集过程，即地震记录可由地震子波和反射系数序列褶积外加噪声得到

d(t)=r(t)*w(t)+n(t)

(13)

式中：d(t)为地震记录；r(t)为反射系数；w(t)为震源子波；n(t)为噪声。

在黏弹性介质中，由于地层的吸收衰减作用，导致震源子波会随着传播时间而变化。因此，需要将稳态的褶积模型推广到非稳态

(14)

式中w(τ，t)为时变子波，且有w(0，t)=w(t)。

(15)

将式(14)离散可改写为矩阵向量形式，有

(16)

式中N为地震道的样点数。

式(16)的加入纵向约束的L2范数最优化问题可用目标函数表示为

(17)

式中：W为时变子波矩阵；r为反射系数向量；d为地震记录向量；μt为纵向约束因子，主要作用为调节数据残差项(第一项)与正则化项(第二项)的权重，使得反演结果更合理。通过极小化目标函数，式(17)的解可表示为

r=(WTW+μtI)-1WTd

(18)

3 基于反演的全局约束吸收补偿方法

假设地震剖面S(t，x)只有一个振幅为常数、斜率为p的线性同相轴，转换到频率域，可表示为

S(f，x)=A(f)exp(j2πfxp)

(19)

式中A(f)是同相轴子波的频谱。假设道间隔Δx为常数，有x=(k-1)Δx，k=1，2，…，K，K为地震剖面的道数。对于任意频率f，有

Sk(f)=Aexp(jθk)

(20)

式中θ=2πfpΔx。则第k道可用第k-1道表示为

Sk(f)=a1(f)Sk-1(f)

(21)

式中a1(f)=exp(jθ)。式(21)递归方程是一阶差分方程，也称为一阶自回归(AR)模型。通过该模型，可以沿着空间方向递归地预测有效信号。

类似地，假设在时空域有M条线性同相轴，则可以用一个M阶的差分方程表示为

(22)

式中ai(f)表示信号预测算子，主要用于预测无噪声的叠加复谐波。式(22)也可写成预测误差形式

(23)

式中：g0=1；gi=-ai。

在实际数据中，噪声是不可避免的。假设地震数据由有效信号和噪声叠加而成，即

Dk=Sk+Nk

(24)

将上式代入式(23)，可得

(25)

上式为自回归滑动平均(ARMA)模型，其自回归分量和滑动平均分量相等。

(26)

为了更加准确地预测有效信号，AR模型中的参数M一般需大于实际同相轴的条数，而且M的确定需要应用统计准则。

假设M=2，式(26)可以写成线性方程组

(27)

其矩阵表达式为

D-Ya=N

(28)

为了最小化噪声能量，建立目标函数

(29)

对目标函数求偏导，并令导数为零，有

YHYa=YHD

(30)

式中上标“H”表示共轭转置。需要注意的是，矩阵YHY为Toeplitz矩阵，因此可用Levinson递归法快速求解。为了获得稳定的滤波器系数，通常需要在Toeplitz矩阵的对角线上增加一个小的扰动μ1，则式(30)改写为

(YHY+μ1Ι)a=YHD

(31)

那么待求的信号预测算子为

a=(YHY+μ1I)-1YHD

(32)

每一个频率的预测算子是相互独立的，将每个频率的一维预测算子依次排列，可得一维预测算子

(33)

式中I为频点数。将每个频率的一维预测算子作用于含随机噪声的数据，可预测出有效信号

(34)

(35)

可见，该信号预测误差算子除了在输出位置的系数为-1以外，其他系数与信号预测算子的系数完全相同。地震信号的空间可预测性主要来源于地震同相轴的空间展布特征(即构造特征)，信号预测误差算子直接反映了地震反射结构特征。需要注意的是，预测误差算子作用于地震数据的结果不是有效信号，而是预测误差。

现有的纵向约束吸收补偿并未考虑到地震道横向之间的关系。为了弥补这一缺陷，将地震信号的反射结构特征(即信号预测误差算子)也作为一种约束引入到吸收补偿算法。基于反演的全局约束吸收补偿方法具体计算过程如下。

(1)通过解最小二乘的问题

(36)

(2)由预测算子可得到预测误差算子

P=(a-M，…，a-2，a-1，-1，a1，a2，…，aM)

(37)

(3)通过反演

(38)

计算去噪地震数据。式中μ2为权衡参数。通过多次迭代求解式(38)，可得频率域去噪数据。

(39)

式中μx为横向约束因子。新的目标函数的解可表示为

(40)

本文提出的基于反演的全局约束多道吸收补偿方法的横向约束算子由原始的地震剖面求得，包含了原始数据的反射结构信息，能够较明显地提高吸收补偿结果的稳定性和抗噪性，同时还保持了同相轴的横向连续性。

4 应用

4.1 Marmousi模型试算

为验证本文方法相对于纵向约束吸收补偿方法的优越性，利用Marmousi模型进行测试。

图1a是原始反射系数与主频为30Hz的Ri-cker子波褶积得到的Marmousi模型模拟数据，可以看出，Marmousi模型构造复杂，既有平滑的同相轴，也有大型的背斜构造和一些小型断裂。图1b是Marmousi模型的衰减数据，其中Q设为50。图1c和图1d是信噪比(SNR)分别为5和20的加噪衰减数据，可以看出，由于地层吸收的影响，地震资料的分辨率降低，原有的构造反射特征变模糊。另外，随机噪声污染了地震反射信号。

分别采用传统的纵向约束方法和本文的全局约束方法对图1b所示的无噪衰减数据进行补偿，结果如图2所示。其中纵向约束吸收补偿过程中μt=0.1，本文方法中μt=0.1，μx=1.0。比较图2a与图2b可以看出，在无噪情况下，纵向约束吸收补偿方法和本文方法均能很好地恢复地震波能量，提高地震资料的分辨率。对比图中红色箭头处可知，本文方法的结果横向更加连续、平滑，构造结构更清晰，弱反射得到了较好的恢复，证明了本文方法的效果更好。图3为无噪衰减数据两种方法补偿前、后单道频谱对比，可以看出，两种方法都有效补偿了中高频的振幅衰减，并且拓宽了原始数据的地震频带，但本文方法补偿结果的频带略宽，包含更丰富的高频有效信息。

图1 Marmousi模型数据(a)原始模拟数据； (b)衰减模拟数据； (c)SNR=5的加噪衰减数据； (d)SNR=20的加噪衰减数据

图2 无噪衰减数据两种方法补偿结果(a)纵向约束补偿； (b)本文方法

对于含噪衰减数据，纵向约束吸收补偿结果虽然恢复了地震波能量，提高了地震资料的分辨率，但是吸收补偿结果的信噪比较低，地震同相轴的横向连续性较差(图4a)。本文方法的补偿结果不仅能够提高地震资料的分辨率，使得构造反射结构更清晰，同时也能保证补偿结果的信噪比，横向上更加连续(图4b)，如图中红色箭头所示，本文方法恢复出了更多的细节，表明本文方法比纵向约束吸收补偿方法具有更强的稳定性和抗噪性。

图5为SNR=5含噪衰减数据利用两种方法补偿结果的单道对比，可以看出，本文方法补偿后与原始数据更相似，无论振幅还是相位，纵向补偿效果均较好(红色箭头所示)。图6为图5结果的频谱对比，可以看出，在含噪情况下，两种方法均能补偿中高频的振幅衰减，并且拓宽了补偿前数据的地震频带。此外，纵向约束补偿结果的高频部分略微强于本文方法补偿结果，有可能是前者对高频噪声的过度放大所致。

图3 无噪数据两种方法补偿前、后单道频谱对比(a)第100道； (b)第103道； (c)第123道

图4 SNR=5(左)和SNR=20(右)含噪数据两种方法补偿结果对比(a)纵向约束方法； (b)本文方法

以均方根振幅为衡量地震剖面横向连续性，图7为SNR=5含噪衰减数据两种方法补偿结果的横向连续性对比，可见本文方法补偿结果横向上与补偿前数据更加接近，横向上抖动较小(红色箭头所示)，体现了本文方法在保持剖面横向连续性和稳定性方面的优势。

图5 SNR=5含噪数据两种方法补偿结果单道对比(a)第100道； (b)第103道； (c)第123道

图6 SNR=5含噪数据两种方法补偿结果单道频谱对比(a)第100道； (b)第103道； (c)第123道

图7 SNR=5含噪数据两种方法补偿结果的横向连续性对比

4.2 实际数据应用

图8a为二维实际地震叠加剖面，采用纵向约束吸收补偿方法和本文方法的处理结果分别如图8b、图8c所示。可以看出，纵向约束补偿方法虽然可以提高地震资料的分辨率，但补偿后地震剖面中的噪声被放大，信噪比较低；经本文方法处理后，地震资料的能量衰减得到有效补偿，地震分辨率得到一定提升，而且补偿后地震记录信噪比较高，噪声得到有效压制，地震剖面的同相轴更连续。由此可见，本文提出的基于反演的全局约束吸收补偿方法比传统的纵向约束吸收补偿方法具有更强的稳定性和抗噪性，其结果更利于后续的解释和反演。

图9为图8数据的局部放大对比，可以看出，相较于原始剖面(图9a)，本文方法补偿剖面(图9c)同相轴增多、更加清晰，提高了分辨率和信噪比，构造更加清楚可辨；相较于纵向约束方法补偿结果(图9b)，本文方法补偿结果信噪比较高，横向更连续。

图10为实际地震剖面两种方法补偿后的单道频谱对比，可以看出，两种方法均能有效补偿地震数据中高频振幅衰减，一定程度上拓宽了数据的地震频带，其中纵向约束方法处理结果高频部分相对较强，是高频噪声放大效应的体现。