0 引言

地震反演技术是提取地层属性的有效手段之一，其基本过程是根据地表接收的地震数据反推地层的物理参数，从而预测储层及识别流体^[1]。根据所使用的地震数据体，地震反演分为叠后反演和叠前反演两大类^[2-4]。

叠后波阻抗(acoustic impedance, AI)反演是地震反演中的一种重要技术，首次被Lindseth^[5]提出，随后在地球物理反演领域得到快速发展^[6-8]。根据实现方式的不同，AI反演主要分为稀疏脉冲反演、基于模型约束的反演、递归反演以及有色反演四类，不同方法各有其适用条件及应用范围^{[2, 9]}。尽管AI反演历经多年且得到很大发展，但仍存在诸多问题，如反问题求解的不适定性、常规反演算法对地层边界刻画不清晰以及递归反演中的误差累积等^[10-11]。

对于反问题求解的不适定问题，人们进行了研究并提出了相应的解决办法^[12-15]。早在1973年，Claerbout等^[16]分析了L₂范数和L₁范数在反问题求解中的作用。基于L₂范数约束，Tikhonov等^[17]提出了一套完整的正则化方法解决不适定问题，其基本思想是在满足数据残差最小的条件下找到一个使模型的L₂范数最小的解。尽管Tikhonov正则化很好地解决了不适定问题，但基于Tikhonov正则化获得的解较平滑，难以满足地学界对反演分辨率的要求。Buland等^[18]将贝叶斯理论引入地球物理反问题，为基于概率的反演方法提供了数学基础。Tarantola^[19]给出了反问题求解的一些算法细节，包括蒙特卡洛方法以及最小二乘离散问题等。此后，各种范数被提出以在解决不适定问题的同时获得较高分辨率的反演结果，包括L₁范数^[20]、L_p范数^[21]以及L_1-2范数等。L_1-2范数作为一种稀疏度量，近年来广泛用于地球物理反演^{[11, 22]}。

此外，常规的递归反演算法首先利用叠后地震数据反演反射系数，然后利用道积分技术计算地层波阻抗信息。这种波阻抗反演算法存在两个问题：一是对初始波阻抗依赖过高，如果第一层的波阻抗给定不准确，那么所有的波阻抗反演结果都会产生偏差；二是道积分中存在误差累积，如果某些反演的反射系数存在误差，那么这种误差会累积到该目标层之后的所有波阻抗计算结果上。

为了解决波阻抗反演中存在的这些问题，本文提出了基于L_1-2正则化的地震波阻抗“块”反演方法。在前人的基础上，将L_1-2正则化引入基于模型的波阻抗反演，通过借鉴全变分(Total Variation, TV)正则化的思想，利用叠后地震数据直接获得波阻抗反演结果。首先，推导线性化的波阻抗正演近似公式并分析精度。然后，基于贝叶斯理论，引入L_1-2正则化构建波阻抗反演的目标函数，利用迭代重加权最小二乘算法(IRLS)求解目标函数，获得波阻抗反演结果。由于波阻抗反演为单道反演算法，反演多道数据时道与道之间会产生空间不连续现象，因此文中对反演结果执行f-x域空间预测滤波，改善由噪声和单道反演算法带来的空间不连续性^[23-24]。最后，利用合成数据和实际资料反演验证所提方法的有效性和可行性。

1 线性化波阻抗正演公式

基于褶积模型，地震数据可以表示为反射系数和地震子波褶积的形式

$ d(t) = \int_{ - \infty }^{ + \infty } w (\tau )r(t - \tau ){\rm{d}}\tau + e(t) $

(1)

式中：d(t)为地震数据，t为时间；w(τ)为地震子波，τ为时移量；r(t－τ)为地层反射系数；e(t)为噪声。将式(1)写为矩阵形式，有

$ \mathit{\boldsymbol{D}} = \mathit{\boldsymbol{WR}} + \mathit{\boldsymbol{e}} $

(2)

式中的变量分别为式(1)变量的矩阵形式。根据定义，当地下为多层水平介质时，任意第i个界面的反射系数为

$ {r_i} = \frac{{{Z_{i + 1}} - {Z_i}}}{{{Z_{i + 1}} + {Z_i}}} $

(3)

式中Z_i和Z_i+1分别为第i层和第i+1层的波阻抗。将式(3)代入式(2)即可获得由地层波阻抗得到地震数据的公式。不难看出，式(3)中反射系数与地层波阻抗之间为非线性关系，当与地震子波褶积后，得到的地震记录与地层波阻抗之间的关系显然也是非线性的，这种非线性关系会增加反演方程的计算复杂度，也会增强反演的不稳定性。因此需要将式(3)进行线性化近似。根据等价无穷小代换e^x~1+x，由式(3)可得^[25]

$ \frac{{{Z_{i + 1}}}}{{{Z_i}}} = \frac{{1 + {r_i}}}{{1 - {r_i}}} \approx \frac{{{{\rm{e}}^{{r_i}}}}}{{{{\rm{e}}^{ - {r_i}}}}} = {{\rm{e}}^{2{r_i}}} $

(4)

则有

$ {r_i} = \frac{1}{2}\left( {\ln {Z_{i + 1}} - \ln {Z_i}} \right) $

(5)

将式(5)代入式(2)可得波阻抗正演公式

$ \mathit{\boldsymbol{D}} = \mathit{\boldsymbol{WSZ}} + \mathit{\boldsymbol{e}} $

(6)

式中：Z=ln(Z₁, Z₂, …, Z_nt)^T为地层波阻抗的自然对数, nt为波阻抗采样点数；S为差分算子，可表示为

$ \mathit{\boldsymbol{S}} = \frac{1}{2}\left( {\begin{array}{*{20}{c}} { - 1}&1&{}&{}&{}\\ {}&{ - 1}&1&{}&{}\\ {}&{}& \ddots & \ddots &{}\\ {}&{}&{}&{ - 1}&1 \end{array}} \right) $

(7)

由式(3)得到式(5)的过程由于使用了等价无穷小代换，因此有必要分析、对比线性化的波阻抗正演近似公式(式(5))和精确公式(式(3))的正演结果。图 1为由式(3)和式(5)获得的反射系数。由图可见：当反射系数r的绝对值在一定范围(|r|≤0.3)时，式(5)与式(3)的结果几乎完全吻合；当|r|较大时二者的结果才有所差异。因此，当地下反射界面为弱反射界面，即|r|≤0.3时，利用式(5)计算地层的反射系数在理论上不会产生较大误差^[20]。

为了进一步了解式(5)的精度，给出了一个基于实际测井资料(图 2a)的正演结果(图 2)。可见，当反射系数的绝对值|r|≤0.3时(图 2b)，式(5)(图 2d)与式(3)(图 2c)的正演结果几乎没有差异(图 2e)，进一步证明式(5)作为线性化的波阻抗正演近似公式具有较高精度。

2 波阻抗反演

根据贝叶斯理论，已知观测数据d时，模型m的后验概率分布P(m|d)等于模型参数的先验概率分布P(m)与似然函数P(d|m)的乘积再除以边缘概率密度P(d)，即

$ P(\mathit{\boldsymbol{m}}\mid \mathit{\boldsymbol{d}}) = \frac{{P(\mathit{\boldsymbol{d}}\mid \mathit{\boldsymbol{m}})P(\mathit{\boldsymbol{m}})}}{{P(\mathit{\boldsymbol{d}})}} $

(8)

式中P(d)为一个常数，其作用是使后验概率分布的总和等于1。对于波阻抗反演，假设噪声服从高斯分布，则似然函数可以表示为

$ \begin{array}{l} P(\mathit{\boldsymbol{D}}\mid \mathit{\boldsymbol{Z}}) = \\ \;\;\;\;{P_0}\exp \left[ { - \frac{1}{2}{{(\mathit{\boldsymbol{D}} - \mathit{\boldsymbol{WSZ}})}^{\rm{T}}}\mathit{\boldsymbol{C}}_e^{ - 1}(\mathit{\boldsymbol{D}} - \mathit{\boldsymbol{WSZ}})} \right] \end{array} $

(9)

式中：C_e为噪声的协方差。阻抗模型Z的先验分布可表示为

$ P(\mathit{\boldsymbol{Z}}) = {P_{0z}}\exp [ - R(\mathit{\boldsymbol{Z}})] $

(10)

式中：R(Z)为Z的函数，不同的先验分布对应不同的R(Z)；P₀和P_0Z分别为两个常数，其作用是使概率分布的总和为1。联立式(9)和式(10)并假设噪声独立分布(即C_e=σ_e²I，σ_e为噪声的标准差，I为单位矩阵)，即可得到波阻抗的后验概率分布

$ \begin{array}{l} P(\mathit{\boldsymbol{Z}}\mid \mathit{\boldsymbol{D}}) \propto \\ \exp \left[ { - \frac{1}{{2\sigma _e^2}}{{(\mathit{\boldsymbol{D}} - \mathit{\boldsymbol{WSZ}})}^{\rm{T}}}(\mathit{\boldsymbol{D}} - \mathit{\boldsymbol{WSZ}}) - \lambda R(\mathit{\boldsymbol{Z}})} \right] \end{array} $

(11)

其中λ为超参数。式(11)忽略了边缘概率密度P(D)。使式(11)中指数项最大的解就是最大后验概率解，因此可以建立反演目标函数

$ J(\mathit{\boldsymbol{Z}}) = \frac{1}{{2\sigma _e^2}}{(\mathit{\boldsymbol{D}} - \mathit{\boldsymbol{WSZ}})^{\rm{T}}}(\mathit{\boldsymbol{D}} - \mathit{\boldsymbol{WSZ}}) + \lambda R(\mathit{\boldsymbol{Z}}) $

(12)

引入高斯先验分布，则R(Z)可表示为

$ R(\mathit{\boldsymbol{Z}}) = \frac{1}{2}\left\| {\mathit{\boldsymbol{Z}} - {\mathit{\boldsymbol{Z}}_0}} \right\|_2^2 $

(13)

式中Z₀为初始模型。在反演合成数据时，通过平滑真实波阻抗得到初始模型；在反演实际数据时，一般通过插值、外推测井数据得到初始模型。基于高斯先验分布的反演结果一般较光滑，为此前人利用其他分布提高反演分辨率，如柯西分布、拉普拉斯分布(L₁范数)以及L_1-2范数。图 3为不同分布的一维概率密度函数对比。由图可见：①高斯分布的峰值最宽，且没有“长尾”分布(“长尾”保证了反演结果中允许存在某些“大”值，即保证“稀疏”，如图中红色框位置)，所以基于高斯分布的反演结果较光滑且抗噪性差。②柯西分布和拉普拉斯分布在提升反演结果稀疏度方面明显优于高斯分布，二者的概率密度函数的峰值更尖锐且有明显的“长尾”。③L_1-2范数的概率密度函数的峰值更尖锐，且“长尾”更明显，因此在理论上可以获得更稀疏的反演结果。

为了更直观地对比几种不同的分布，绘制了L₀范数、L₁范数(对应拉普拉斯分布)、L₂范数(对应高斯分布)以及L_1-2范数的二维图像(图 4)。可见，相比于L₂范数(图 4c)和L₁范数(图 4b)，L_1-2范数(图 4d)与L₀范数(图 4a)更相似，因此理论上基于L_1-2范数约束可以获得比L₂范数和L₁范数约束更稀疏、分辨率更高的反演结果。

以二维向量(x，y)为例，几种范数的定义为

$ \left\{ {\begin{array}{*{20}{l}} {{L_0} = \left\{ {\begin{array}{*{20}{c}} 0&{(x, y) = 0}\\ 2&{x \ne 0, y \ne 0}\\ 1&{{\rm{其他}}} \end{array}} \right.}\\ {{L_2} = \sqrt {|x{|^2} + |y{|^2}} }\\ {{L_1} = |x| + |y|}\\ {{L_{1 - 2}} = |x| + |y| - \sqrt {|x{|^2} + |y{|^2}} } \end{array}} \right. $

(14)

式中x、y为二维向量的分量。L₀范数定义为计算模型中非零项的个数，也被认为是最稀疏的度量。因此可通过对比其他范数图像与L₀范数图像的相似程度判断范数的稀疏性。

当引入L_1-2正则化时，在式(13)的基础上，R(Z)可以表示为

$ R(\mathit{\boldsymbol{Z}}) = \frac{\eta }{2}{\left\| {\mathit{\boldsymbol{Z}} - \mathit{\boldsymbol{Z}}} \right\|_0}_2^2 + \beta \left( {{{\left\| \mathit{\boldsymbol{Z}} \right\|}_1} - \alpha {{\left\| \mathit{\boldsymbol{Z}} \right\|}_2}} \right) $

(15)

式中：α≥0为权重参数，该值越大，L_1-2范数越稀疏；η和β则分别控制高斯项和L_1-2正则化项在目标函数中的比例。由于L_1-2范数是一种稀疏范数，而波阻抗并不满足稀疏性假设，因此将式(15)代入式(12)进行反演必定会导致错误的反演结果。虽然波阻抗不满足稀疏假设，但波阻抗的空间导数满足稀疏条件，因此类似于TV正则化，本文通过对波阻抗的空间导数施加稀疏约束，从而获得高分辨率的反演结果。这样，基于L_1-2正则化的波阻抗反演的目标函数可定义为

$ \begin{array}{l} J(\mathit{\boldsymbol{Z}}) = \frac{1}{2}\left\| {\mathit{\boldsymbol{D}} - \mathit{\boldsymbol{WSZ}}} \right\|_2^2 + \frac{\eta }{2}\mathit{\boldsymbol{Z}} - {\mathit{\boldsymbol{Z}}_0}_2^2 + \\ \;\;\;\;\;\;\;\;\;\beta \left( {{{\left\| {\mathit{\boldsymbol{SZ}}} \right\|}_1} - \alpha {{\left\| {\mathit{\boldsymbol{SZ}}} \right\|}_2}} \right) \end{array} $

(16)

在反演过程中，通常需要多次试验获得参数η和β，以保证最终获得较满意的反演结果。利用IRLS算法求解式(16)，将式(16)对Z求导可得

$ \begin{array}{l} J(\mathit{\boldsymbol{Z}}) = \frac{1}{2}\left\| {\mathit{\boldsymbol{D}} - \mathit{\boldsymbol{WSZ}}} \right\|_2^2 + \frac{\eta }{2}\mathit{\boldsymbol{Z}} - {\mathit{\boldsymbol{Z}}_0}_2^2 + \\ \;\;\;\;\;\;\;\;\;\beta \left( {{{\left\| {\mathit{\boldsymbol{SZ}}} \right\|}_1} - \alpha {{\left\| {\mathit{\boldsymbol{SZ}}} \right\|}_2}} \right) \end{array} $

(17)

令$\partial J/\partial \mathit{\boldsymbol{Z}} = 0$，可得最终的迭代更新式为

$ \begin{array}{l} {\mathit{\boldsymbol{Z}}^{n + 1}} = {\left[ {{{(\mathit{\boldsymbol{WS}})}^{\rm{T}}}(\mathit{\boldsymbol{WS}}) + \eta \mathit{\boldsymbol{I}} + \mathit{\boldsymbol{w}}} \right]^{ - 1}} \times \\ \;\;\;\;\;\;\;\;\;\;\left[ {{{(\mathit{\boldsymbol{WS}})}^{\rm{T}}}\mathit{\boldsymbol{D}} + \eta {\mathit{\boldsymbol{Z}}^n}} \right] \end{array} $

(18)

式中：Zⁿ为第n次迭代的解；w为重加权矩阵，可表示为

$ w = \beta {\mathit{\boldsymbol{S}}^{\rm{T}}}\left[ {\frac{I}{{{\rm{diag}}\left( {\left| {\mathit{\boldsymbol{S}}{\mathit{\boldsymbol{Z}}^{n - 1}}} \right|} \right)}} - \alpha \frac{I}{{\mathit{\boldsymbol{S}}{\mathit{\boldsymbol{Z}}^{n - 1}}_2}}} \right]\mathit{\boldsymbol{S}} $

(19)

利用IRLS算法求解式(16)的算法流程为：

(1) 设定超参数η和β、权重参数α以及最大迭代次数K；

(2) 设定初始模型Z₀；

(3) 给定子波矩阵W及差分矩阵S；

(4) 令n=1，根据式(19)计算w；

(5) 根据式(18)进行求解，并令n=n+1；

(6) 重复步骤(4)和步骤(5)，直至达到最大迭代次数。

值得注意的是，K的选择对于反演结果的收敛以及算法的效率极为重要。一般来说，迭代10次左右即可获得满意的反演结果，但K的选择还与模型的复杂度以及数据的信噪比有关。在反演实际地震数据时，可通过对比井旁道的反演结果与实际测井资料设定较合适的K。

3 模型数据测试

利用单道块状模型和二维SEG/EAGE推覆体模型，分别对比、分析基于L₂范数、L₁范数以及L_1-2范数约束的波阻抗反演结果，以验证所提方法的效果。在反演推覆体模型时，对反演结果执行f-x域空间预测滤波以提高反演结果的横向连续性。文中信噪比定义为

$ {\rm{SNR}} = 10\lg \frac{{\left\| {{\mathit{\boldsymbol{s}}^ * }} \right\|_2^2}}{{\left\| {\mathit{\boldsymbol{s}} - {\mathit{\boldsymbol{s}}^ * }} \right\|_2^2}} $

(20)

式中：s^*为无噪地震数据；s为含噪地震数据。

利用式(3)得到单道块状模型反射系数，与主频为30Hz的雷克子波褶积得到合成地震记录(图 5)，然后利用合成地震记录反演得到波阻抗(图 6)。从结果来看：①由于没有施加稀疏约束，且数据含噪，对于块状模型，基于L₂范数约束的反演效果明显较差(图 6a、图 6d)。②基于L₁范数(图 6b、图 6e)和L_1-2范数(图 6c、图 6f)约束的反演结果均呈明显的“块状”特征，反演效果优于L₂范数，其中L_1-2范数约束的反演结果与真实模型吻合度更高(图 6b、图 6c、图 6e、图 6f的红色箭头处)。③随着噪声能量增强，反演质量均下降。如基于L₂范数约束的反演结果仍然无法精确刻画地层的“块状”特征(图 6d)；基于L₁范数约束的反演结果与真实值也存在偏差(图 6e)；基于L_1-2范数约束的反演结果虽然在某些层位与真实值稍有差别，但总体上吻合较好，且很好地刻画了波阻抗的“块状”特征(图 6f)，证明了所提方法的有效性。

图 7为基于图 6数据计算的反射系数。由图可见：①无论信噪比为10dB还是4dB，基于L₂范数约束的反演效果均较差，反演结果与真实的反射系数相差很大。这是由于L₂范数为平滑约束，其概率分布函数在均值μ附近的峰值较宽(图 3)，因此反演结果不能很好地压制一些“小值”，即抗噪性较差，而且概率密度函数“长尾”不明显，从而使某些离散点被压制，因此不能反演幅值较大的反射系数。②L₁范数和L_1-2范数约束的反演效果明显变好，其中L_1-2范数由于更好地压制了“小值”，同时也保护了“大值”，因此基于L_1-2范数约束的反演效果优于基于L₁范数约束，尤其是当信噪比为4dB时(图 7b红色箭头处)。

利用式(3)得到SEG/EAGE推覆体模型(图 8b)反射系数，与主频为40Hz的雷克子波褶积并加入信噪比为4dB的随机噪声得到合成地震记录(图 8a)。图 8为基于SEG/EAGE推覆体模型的反演结果。从反演结果可见：①L_1-2范数约束的反演结果的分辨率与精度最高(图 8f)，与真实模型(图 8b)吻合很好；基于L₁范数约束的反演效果略差(图 8e)，但仍然好于基于L₂范数约束(图 8d)。②由于单道算法的限制，反演结果中存在明显的横向不连续现象，尤其是当加入稀疏约束时，这种现象更明显(图 8e、图 8f)。因此，有必要改善反演结果的横向连续性。③三种范数约束的波阻抗反演结果(图 8d~图 8f)经过f-x域空间预测滤波后，横向连续性均有提升(图 8g~图 8i)，尤其是图 8h和图 8i。该结果进一步证明了本文所提方法的有效性和可行性。

4 实际资料测试

基于M区实际地震资料测试本文方法。图 9为M区叠后地震资料反演结果。可见：①基于L₂范数约束的反演结果分辨率较低(图 9d)。②基于L₁范数(图 9e)和L_1-2范数约束(图 9f)的反演效果得到明显改善，反演分辨率更高，对地层的刻画更清晰，且后者优于前者(红色箭头处)。图 10为井旁道反演结果与经过Backus平均的实际测井资料对比结果。可见，基于L_1-2范数约束的井旁道反演结果与测井曲线更吻合，尤其表现在对“块状”特征的反演(红色箭头处)。此外，通过求取目的层段(3.0~3.2s)反演结果与真实测井数据的相关系数，可以定量对比不同范数约束的反演结果差异。向量X与Y的相关系数定义为

$ C = \frac{{c(\mathit{\boldsymbol{X}}, \mathit{\boldsymbol{Y}})}}{{\sqrt {c(\mathit{\boldsymbol{X}}, \mathit{\boldsymbol{X}}) \times c(\mathit{\boldsymbol{Y}}, \mathit{\boldsymbol{Y}})} }} $

(21)

式中：c(X，Y)为X与Y的协方差；c(X，X)和c(Y，Y)分别为X和Y的方差。

基于L₂范数、L₁范数及L_1-2范数的反演结果与真实测井数据的相关系数分别为0.6529、0.5197、0.7264，表明：基于L_1-2范数的反演结果与真实测井数据更吻合，这是由于L_1-2范数约束为稀疏约束(图 10的红色箭头处)，可很好地反演地层的“块状”信息；基于L₂范数和L₁范数约束的反演效果则相对较差，其中基于L₁范数约束的反演结果的相关系数小于基于L₂范数约束，可能原因是前者的反演结果尽管在稀疏性方面要优于后者，但前者的反演结果在大约3.14s处存在一个明显的层位漂移，从而导致反演结果与真实测井数据的相关系数较低。

5 结论

尽管波阻抗反演技术已经相当成熟，但仍然存在一定问题，如反问题的不适定性、反演的分辨率低以及对地层边界刻画不清晰等。针对这些问题，本文提出了基于L_1-2正则化的地震波阻抗“块”反演方法。首先推导线性化的波阻抗正演近似公式并进行精度分析，证明近似式的精度较高。然后介绍了基于L_1-2范数约束的波阻抗反演方法，与传统的基于L₂范数和L₁范数约束的反演方法相比，基于L_1-2范数约束的反演结果明显提高了分辨率，反演结果更具“块状”特征。相关系数的定量对比也证明了基于L_1-2范数的反演结果优于基于L₁范数和L₂范数。此外，文中引入f-x域空间预测滤波改善单道反演算法带来的横向不连续现象。显而易见的是，通过对反演结果执行f-x域空间预测滤波，明显提高了反演结果的横向连续性。合成和实际地震数据反演证明了所提方法的可行性和有效性。