0 引言

反射系数由地下介质的速度和密度共同决定，地震数据可表达为反射系数和地震子波的卷积^[1]，因此反射系数被视为地震数据与地下地质结构之间的桥梁^[2]。近些年，随着稀疏约束优化算法的发展，使反射系数估计成为一种常规处理方法，从而促进反射系数反演在薄层识别^[3-4]、强反射屏蔽消除^[5]、非平稳地震信号解释^[6]、小构造检测^[7]及块化波阻抗反演^[8]等领域得到广泛应用。

三维反射系数数据体可以通过各类稀疏优化算法进行逐道反演获得，如：Taylor等^[9]借助线性规划方法提出了L₁范数反褶积方法；Velis^[10]利用快速模拟退火算法提出了随机稀疏脉冲反褶积方法以估计反射系数序列；Dondurur^[2]提出的黄金反褶积算法可以获得混合相位子波时的反射系数估计；Zhang等^[11]基于基追踪提出了反射系数反演方法；Yuan等^[12]将褶积问题转化到频率域，利用稀疏贝叶斯优化方法实现了反射系数谱反演。上述逐道反射系数反演方法只采用了各种纵向稀疏约束策略，没有利用相邻地震道之间的高度相关性，因此在低信噪比情况下获得的反射系数反演结果常常呈抖动假象，横向连续性差，与实际情况不符(实际沉积岩地层阻抗界面在横向上呈平滑起伏特征)，一些对油气储运至关重要的小构造(如小隆起、薄层等)很容易被掩盖。此外，逐道反演策略在处理三维地震数据时需要逐道分别经历一次迭代优化过程，期间反复进行时间域的矩阵乘法和求逆运算，计算效率相对较低。

为降低噪声引起的抖动效应，近年来，通过在反演目标函数中引入与邻道信息相关的约束项，发展出一些有效的多道反射系数反演方法，如：Heimer等^[13]利用Markov-Bernoulli随机场模拟反射系数剖面的局部连续性；Kazemi等^[14]假设各道对应的子波相同，建立了多道反射系数稀疏反演目标函数；Yuan等^[15]利用时间—空间域中线性同相轴在f-k域中的稀疏性，提出了f-k域多道反射系数反演方法；Wang等^[16]在反演目标函数中引入一种预测算子，提高了反射系数横向连续性；Ma等^[17]利用协方差矩阵描述地震道之间的空间相关性，在稀疏贝叶斯反演框架下压制了随机噪声引起的抖动效应；Du等^[18]利用地震数据中获得的结构特征设计了横向平滑正则化算子并用于稀疏反褶积；纪永祯等^[19]通过引入相邻道数据作为约束，将Velis^[10]提出的单道随机反射系数反演方法发展到多道；印兴耀等^[20]基于反射系数与地震数据在结构上具有相似性的特点，利用地震数据互相关描述地层反射系数的空间结构特征，提高了反演结果的横向连续性。但是，加入各种横向约束后的反演目标函数在优化过程中需要在时间—空间域处理更大尺度的矩阵乘法和求逆运算^[20]，求解效率会进一步降低。而且这些多道反演方法只针对二维地震剖面，大多难以直接拓展到三维地震数据体，只能通过逐线反演方式获得三维反射系数，线与线之间缺乏联系，易在联络测线方向产生抖动假象。

为了消除逐道或逐线反演获得的三维反射系数中的抖动假象，应将地震道在三维空间中与相邻道横向上的高度相关性引入反演目标函数，进而开发能够高效求解的三维反演算法。全变分(Total Variation，TV)正则化正是借助待求参数中相邻数据间的高度相似性来提高反演结果横向连续性的方法，由于其实质是一种稀疏约束，因此在提高横向连续性的同时能够保持断点、边界点等稀疏的“尖锐”信息^[21]。但常规TV正则化采用相邻地震道间的一阶导数构造约束项(即前、后两道之差)，在最小化TV函数时会压制相邻道之间微小的振幅差异，造成反演结果中的“阶梯”效应^[22-23]。本文采用合成地震数据横向二阶导数构造TV正则化项，有效信号横向二阶导数体现振幅横向变化的快慢，而地震振幅在横向上呈现局部线性变化特征，因此横向二阶导数主要突出不规则噪声，对其TV函数进行最小化时的保幅性优于常规TV正则化。此外，在分裂Bregman优化框架下给出了三维反演目标函数在频率—波数域中的快速求解方法，借助频率—波数域点乘避开了三维卷积和求逆运算，因此本文方法在压制抖动假象的同时能够兼顾运算效率。

1 方法原理 1.1 三维反射系数反演目标函数

利用符号“*”表示地震子波列向量w与三维反射系数R的三维卷积，则三维地震数据可表示为

$ \boldsymbol{S}=\boldsymbol{w} * \boldsymbol{R}+\boldsymbol{n} $

(1)

式中n为噪声。在最小二乘意义下建立如下三维反射系数反演目标函数

$ \mathit{\Phi}_1(\boldsymbol{R})=\frac{1}{2}\|\boldsymbol{S}-\boldsymbol{w} * \boldsymbol{R}\|_{\mathrm{F}}^2+\mu\|\boldsymbol{R}\|_1 $

(2)

式中：‖·‖_F表示Frobenius范数，即三维数据体中各元素平方和的算术平方根；‖·‖₁表示L₁范数，即三维数据体中各元素绝对值之和；μ为稀疏约束正则化参数。将三维数据体替换为单道地震信号，式(2)就退化为常规逐道反射系数反演的目标函数。

1.2 横向二阶导数TV正则化目标函数

尽管式(2)将各道数据包含在同一个反演目标函数中，可以直接求解获得三维反射系数，但该目标函数缺乏横向约束信息，无法压制反射系数的抖动，因此本文在式(2)中加入合成数据的横向二阶导数构造TV正则化项，具体原因可用含噪地震数据(图 1a)进行说明。图 1b为该地震剖面3.0 s时的振幅曲线，红色框中CDP号为150~165，可以发现振幅在横向上呈现局部线性化特征，因此横向一阶导数剖面中(图 1c)除了噪声信息外还存在大量黑色箭头所指的体现前、后地震道振幅差异的规则信息。横向二阶导数剖面(图 1d)中主要突出不规则噪声，基本不含规则的有效信息。构造横向二阶导数TV正则化项的方式为：首先令U=w*R表示合成的三维地震数据，则U在xOy平面内的二阶导数包括：x方向二阶偏导数U_xx=D_xx*U、y方向二阶偏导数U_yy=D_yy*U、混合偏导U_xy=U_yx=D_xy*U，其中D_xx、D_yy和D_xy分别表示对应的二阶差分算子(图 2)。

构造TV正则化项前，首先定义如下组合算子

$ \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{ \boldsymbol{\varDelta} }} = {{\boldsymbol{D}}_{xx}} + {{\boldsymbol{D}}_{yy}}}\\ {\mathit{\boldsymbol{ \boldsymbol{\varLambda} }} = {{\boldsymbol{D}}_{xx}} - {{\boldsymbol{D}}_{yy}}}\\ {\mathit{\boldsymbol{ \boldsymbol{\varGamma} }} = 2{{\boldsymbol{D}}_{xy}}} \end{array}} \right. $

(3)

为了全方位体现相邻道在三维空间中的相关性，构造的TV正则化项为

$ T({\boldsymbol{R}}) = \frac{{\sqrt 2 }}{2}\sum {\sqrt {{{\left( {{\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}_1}*{\boldsymbol{U}}} \right)}^2} + {{\left( {{\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}_2}*{\boldsymbol{U}}} \right)}^2} + {{\left( {{\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}_3}*{\boldsymbol{U}}} \right)}^2}} } $

(4)

式中：Ψ₁=Δ；Ψ₂=Λ；Ψ₃=Γ。将式(3)代入式(4)后，可得

$ T(\boldsymbol{R})=\sum \sqrt{\left(\boldsymbol{D}_{x x} * \boldsymbol{U}\right)^2+\left(\boldsymbol{D}_{y y} * \boldsymbol{U}\right)^2+2\left(\boldsymbol{D}_{x y} * \boldsymbol{U}\right)^2} $

(5)

当只在单一方向上进行横向约束时，式(5)的TV正则化函数退化为横向二阶导数的L₁范数，即T(R)=‖D_xx*U‖₁。

将式(5)的横向约束项加入式(2)，可得到最终的三维反射系数反演目标函数

$ \mathit{\Phi}_2(\boldsymbol{R})=\frac{1}{2}\|\boldsymbol{S}-\boldsymbol{w} * \boldsymbol{R}\|_{\mathrm{F}}^2+\mu\|\boldsymbol{R}\|_1+\lambda T(\boldsymbol{R}) $

(6)

式中λ为横向约束正则化参数。μ值越大反演结果越稀疏，λ值越大横向平滑程度越高，一般利用小样本数据经过多次试验获得最佳参数组合，因此正式反演时将二者视为已知量。

1.3 分裂Bregman优化框架下的快速求解方法

式(6)中的时间—空间域三维卷积运算使其很难被直接优化求解，因此基于分裂Bregman优化框架在频率—波数域推导了快速求解算法^[24-25]。

首先定义中间变量d≈R和v_m≈Ψ_m*U=Ψ_m*w*R(m=1、2、3)，则在分裂Bregman优化框架下求解式(6)等价于交替求解下列反演问题^[24]

$ \begin{gathered} \boldsymbol{R}^{(k+1)}=\min _{\boldsymbol{R}}\left[\frac{1}{2}\left\|\boldsymbol{S}-\boldsymbol{w}^* \boldsymbol{R}\right\|_{\mathrm{F}}^2+\frac{\alpha}{2} \| \boldsymbol{d}^{(k)}-\boldsymbol{R}-\right. \\ \left.\boldsymbol{b}^{(k)}\left\|_{\mathrm{F}}^2+\frac{\beta}{2} \sum\limits_{m=1}^3\right\| \boldsymbol{v}_m^{(k)}-\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}_m * \boldsymbol{U}-\boldsymbol{q}_m^{(k)} \|_{\mathrm{F}}^2\right] \end{gathered} $

(7)

$ \boldsymbol{d}^{(k+1)}=\min\limits _{\boldsymbol{d}}\left[\mu\|\boldsymbol{d}\|_1+\frac{\alpha}{2}\left\|\boldsymbol{d}-\boldsymbol{R}^{(k+1)}-\boldsymbol{b}^{(k)}\right\|_{\mathrm{F}}^2\right] $

(8)

$ \begin{aligned} \boldsymbol{v}_m^{(k+1)}= & \min _{\boldsymbol{v}_m}\left[\lambda\left\|\boldsymbol{v}_1, \boldsymbol{v}_2, \boldsymbol{v}_3\right\|_2+\frac{\beta}{2} \sum\limits_{m=1}^3 \| \boldsymbol{v}_m-\right. \\ & \left.\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}_m * \boldsymbol{U}-\boldsymbol{q}_m^{(k)} \|_{\mathrm{F}}^2\right] \end{aligned} $

(9)

$ \boldsymbol{b}^{(k+1)}=\boldsymbol{b}^{(k)}-\left[\boldsymbol{d}^{(k+1)}-\boldsymbol{R}^{(k+1)}\right] $

(10)

$ \boldsymbol{q}_m^{(k+1)}=\boldsymbol{q}_m^{(k)}-\left[\boldsymbol{v}_m^{(k+1)}-\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}_m * \boldsymbol{U}^{(k+1)}\right] $

(11)

式中：α和β皆为使反演过程稳定的正则化参数；k为迭代次数；$\left\|\boldsymbol{v}_1, \boldsymbol{v}_2, \boldsymbol{v}_3\right\|_2=\sqrt{\boldsymbol{v}_1^2+\boldsymbol{v}_2^2+\boldsymbol{v}_3^2}$表示三个数据体v₁、v₂和v₃相同位置处的三个数的L₂范数；三维中间变量b和q_m初值设置为0。尽管式(7)中反演问题的目标函数是可导的，但在时间域处理该三维问题需要涉及三维卷积和大型矩阵求逆运算。因此基于三维时间—空间域卷积对应频率—波数域点乘运算，本文将其目标函数转化到频率—波数域

$ \begin{aligned} \mathit{\Phi}_3(\hat{\boldsymbol{R}})= & \frac{1}{2}\|\hat{\boldsymbol{S}}-\hat{\boldsymbol{w}} \circ \hat{\boldsymbol{R}}\|_{\mathrm{F}}^2+\frac{\alpha}{2}\left\|\hat{\boldsymbol{d}}^{(k)}-\hat{\boldsymbol{R}}-\hat{\boldsymbol{b}}^{(k)}\right\|_{\mathrm{F}}^2 \\ & \frac{\beta}{2} \sum\limits_{m=1}^3\left\|\hat{\boldsymbol{v}}_m^{(k)}-\hat{\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}}_m \circ \hat{\boldsymbol{w}} \circ \hat{\boldsymbol{R}}-\hat{\boldsymbol{q}}_m^{(k)}\right\|_{\mathrm{F}}^2 \end{aligned} $

(12)

式中：“ $\hat{\boldsymbol{\bullet}}$”表示在三维傅里叶变换域；“$\circ$”表示两个三维数据体的点乘运算。对式(12)求导并令导数为0，可得

$ \left\{\begin{aligned} & \boldsymbol{R}^{(k+1)}= \mathrm{IFFT}_{3 \mathrm{D}}\left[\boldsymbol{P}^{-1} \circ \boldsymbol{Q}^{(k)}\right] \\ & \boldsymbol{P}=\overline{\hat{\boldsymbol{w}}} \circ \hat{\boldsymbol{w}}+\alpha \boldsymbol{I}+\beta \overline{\hat{\boldsymbol{w}}} \circ \hat{\boldsymbol{w}} \circ \sum\limits_m\left(\overline{\hat{\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}}}_m \circ \hat{\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}}_m\right) \\ & \boldsymbol{Q}^{(k)}= \overline{\hat{\boldsymbol{w}}} \circ \hat{\boldsymbol{S}}+\alpha\left[\hat{\boldsymbol{d}}^{(k)}-\hat{\boldsymbol{b}}^{(k)}\right]+ \\ & \beta \overline{\hat{\boldsymbol{w}}} \circ \sum\limits_m\left\{\overline{\hat{\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}}}_m \circ\left[\hat{\boldsymbol{v}}_m^{(k)}-\hat{\boldsymbol{q}}_m^{(k)}\right]\right\} \end{aligned}\right. $

(13)

式中：“IFFT_3D”表示三维傅里叶逆变换；“$\overline{\boldsymbol{\bullet}}$”表示复共轭；I为元素全为1的三维数据体。式(8)和式(9)中的反演问题可分别采用软阈值函数直接求解^[21-22]

$ \begin{aligned} \boldsymbol{d}^{(k+1)}= & \max \left\{\operatorname{abs}\left[\boldsymbol{R}^{(k+1)}+\boldsymbol{b}^{(k)}\right]-\frac{\mu}{\alpha} \boldsymbol{I}, \bf{0}\right\} \circ \\ & \operatorname{sign}\left[\boldsymbol{R}^{(k+1)}+\boldsymbol{b}^{(k)}\right] \end{aligned} $

(14)

$ \boldsymbol{v}_m^{k+1}=\max \left(1-\frac{\lambda / \beta}{\boldsymbol{s}^{(k)}}, \bf{0}\right) \circ\left[\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}_m * \boldsymbol{U}^{(k+1)}+\boldsymbol{q}_m^{(k)}\right] $

(15)

式中：三维参数$\boldsymbol{s}^{(k)}=\sqrt{\sum\limits_{m=1}^3\left[\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}_m * \boldsymbol{U}^{(k+1)}+\boldsymbol{q}_m^{(k)}\right]^2}$(开平方运算表示对三维数据体中的各个元素开平方根)；abs(∙)和sign(∙)分别表示对三维数据体中的各个元素取绝对值和符号；max(∙, 0)表示将三维数据体中小于0的元素用0替换。

实验表明，交替计算式(13)~式(15)、式(10)、式(11)即可在几十次迭代后直接获得优质的三维反射系数反演结果。利用分裂Bregman方法对式(6)进行优化的详细步骤如下。

(1) 输入三维地震数据S，地震子波向量w，正则化参数μ、λ、α和β，复合横向差分算子Ψ_m(m=1、2、3)，最大迭代次数K；定义初始值d⁽⁰⁾=b⁽⁰⁾=v₁⁽⁰⁾=v₂⁽⁰⁾=v₃⁽⁰⁾=q₁⁽⁰⁾=q₂⁽⁰⁾=q₃⁽⁰⁾=0和k=0；

(2) 当k≤K时，依次根据式(13)~式(15)、式(10)、式(11)计算三维数据体R^(k)、d^(k)、v_m^(k)、b^(k)和q_m^(k)的值，并令k+1⇒k；

(3) 重复步骤(2)直到k＞K，并输出三维反射系数反演结果R。

2 方法验证

图 3a为一个断层模型反射系数剖面，该模型上盘顶、底反射系数随CDP号增大而减小，下盘规律相反，因此可用该模型模拟地震反射振幅的横向缓慢变化以及断点。利用主频为30 Hz的Ricker子波与该反射系数褶积并加入20%高斯白噪合成地震剖面(图 3b)。由于噪声的存在，使逐道反演剖面(图 3c)中存在大量随机的小反射系数，且红色箭头处反射系数出现微小错动，反演结果横向连续性差。利用合成数据横向二阶导数的L₂范数(即式(6)中横向约束项为$\left\|\boldsymbol{D}_{x x} * \boldsymbol{U}\right\|_2^2$)进行反演的结果如图 3d所示，尽管随机抖动现象得到较好的压制，但红色箭头处本应是尖锐断点的位置反射系数呈平滑过渡到零的特性，与事实不符。图 3e和3f分别为横向一阶导数和二阶导数TV正则化反演所得反射系数剖面，可以看出，TV正则化稀疏约束不仅可以压制随机抖动，而且断点保持准确。为了比较横向一阶导数和二阶导数TV正则化的保幅性，图 4展示了两种TV正则化方法反演所得上盘和下盘顶、底界面反射系数值与真值的对比。横向一阶导数TV正则化方法所得结果具有明显的“阶梯”效应，有些相邻道之间的差异被抹平，横向二阶导数TV正则化方法所得结果反射系数平滑变化，与真实值吻合度高，表现出更强的保幅性。

分别对三维推覆体模型含噪数据(图 5a)进行逐道反演、逐线反演和三维反演，分析三种反演方法的抖动假象及算法效率。其中逐道反演只在纵向上采用L₁范数稀疏约束，逐线反演在测线方向上增加了合成数据横向二阶导数的L₁范数约束。剖面上，由于随机噪声的干扰使逐道反演结果(图 5b)中水平界面位置出现微小错位，造成明显的抖动假象。逐线反演结果(图 5c)中抖动现象得到一定压制，但压制效果无法与本文三维反演结果(图 5d)媲美。平面上，红色箭头所指位置存在一个小隆起构造，逐道反演结果由于抖动假象严重，如果不仔细观察很容易在平面上忽略该小隆起；逐线反演结果中该处小隆起构造的成像质量得到提高，但隆起边界的连续性相对较差；本文三维反演结果在成像清晰度和边界连续性方面，都明显优于其他两种反演方式。此外，黄色箭头处存在一个波阻抗差异较小的连续反射界面，在逐道反演结果上该界面与抖动干扰十分类似，识别、解释困难，在逐线反演结果上可识别、解释，但横向连续性无法达到本文三维反演的效果。

迭代次数设置为100时，分别用三种反演方法对该三维推覆体模型数据进行反演，处理时间见表 1(处理器型号：Inter(R) Core(TM) i7-9700 CPU @3.00 GHz；RAM：16.0 GB)。逐道反演时不需要涉及与横向约束相关的矩阵运算，且涉及到的矩阵运算维度较小，因此运行效率尚可接受。逐线反演时每次迭代都需要消耗大量机时求解大型矩阵的逆，因此反演效率较低。本文三维反演算法用频率域点乘实现大型矩阵的乘法运算和求逆运算，在12 min内即可完成该三维地震数据的反演。

3 实际应用

以两个实际三维地震数据反射系数反演为例，对比逐道、逐线及本文三维反演算法所得结果，进一步说明本文方法在抖动假象压制和兼顾算法效率方面的优势。图 6a为加拿大阿尔伯塔Erskine地区的小规模三维地震数据体，该数据存在明显的随机干扰，分别采用逐道(图 6b)、逐线(图 6c)和本文三维反演方式(图 6d)获得反射系数数据体。从反演结果可以看出，本文三维反演方法获得的反射系数横向连续性(图 6d)明显好于其他两种反演方法(图 6b、图 6c)，而且粉色箭头处的同相轴错断现象得到较好的保持。此外，尽管本文方法对该实际三维数据中的所有地震道进行同步反演，但由于在频率域中计算避免了大型矩阵乘积和求逆运算，计算效率反而最高(表 1)。

为了进一步说明利用本文方法压制抖动后对小构造的成像优势，在图 7中给出了Erskine三维数据0.41 ms的时间切片对比。粉色箭头位置存在两个隆起构造，由于左侧隆起规模小，在逐道反演时间切片上被随机抖动淹没(图 7b)，逐线反演时间切片上隆起边界刻画模糊(图 7c)，但本文三维反演时间切片上隆起边界成像清晰(图 7d)，有助于人工或机器自动识别和解释。

图 8a展示了东海一斜坡带上X工区的三维叠后地震数据(地层倾向与主测线方向一致)，分别采用时间域逐道(图 8b)、逐线(图 8c)及本文三维反演方式(图 8d)进行反射系数反演。整体上图 8d抖动程度最小，而且粉色箭头处的断点清晰。

图 9为第100线的反演结果，地震剖面上粉色箭头处出现振幅异常，在对应反射系数剖面上显示存在水退引起的地层下超，因此振幅异常是由地层厚度变化引起的。对比反射系数反演结果可以发现，本文三维反演算法对该地质现象成像最清晰(图 9d)。地震剖面上蓝色箭头处为薄层引起的复合波，本文三维反演方法对该薄层成像的横向可追踪性最好，其他两种方法的结果由于存在抖动效应，出现上、下层面相交的假象(图 9b、图 9c)，效果欠佳，不利于薄层厚度的解释。算法效率方面，本文方法可以在5 min左右完成该三维数据的反演，其他两种方法每次迭代都需要在时间域求解矩阵乘法及逆运算，而且不同的道(或线)需要分别经历一次完整的迭代过程，因此二者的计算效率较低(表 1)。

4 讨论

尽管本文提出的三维反演方法在频率—波数域中避免了耗时的三维卷积和求逆运算，但同时引入了三维快速傅里叶正、反变换，对于大规模三维地震数据进行三维快速傅里叶变换时效率很低，当计算机内存不大时，甚至会出现因内存不足而无法计算的情况。为了避免这种情况的发生，可以将大型三维地震数据分为若干小规模三维地震数据块后再分别采用本文方法进行三维反演。从而既可以获得三维反演对抖动假象的全方位压制效果，又不至于引起计算困境。

5 结束语

为了压制三维反射系数反演结果中的抖动假象，本文将合成数据横向二阶导数TV正则化项加入反演目标函数，进而提出了反射系数三维反演方法。通过合成数据和实际数据分析发现，地震数据有效波振幅在局部范围内呈现横向缓慢线性变化的特征，用二阶导数构造TV正则化项的保幅性优于常规一阶导数TV正则化。为解决时间—空间域三维反演的算法效率问题，在分裂Bregman优化框架下将三维反演目标函数转换到频率—波数域，用点乘代替了耗时的三维卷积和求逆运算，实现了快速三维反演。合成地震数据试算和实际地震数据处理结果表明，与时间—空间域逐道反演、逐线反演方法相比，本文三维反演算法能够全方位压制抖动假象，并保持同相轴错断、中断等“尖锐”信息，有利于三维数据体的小构造和薄层的边界刻画和解释。