刘百红 李建华 郑四连

(①中国石化石油物探技术研究院，江苏南京 211103； ②中国石油东方地球物理公司研究院地质研究中心，河北涿州 072751)

1 问题的提出

在地震勘探中，叠后地震道的形成可用褶积模型描述

S(t)=r(t)*W(t)+n(t)

(1)

式中：S(t)表示地震道；r(t)表示反射系数序列；W(t)表示子波；n(t)表示噪声；t为时间； “*”表示褶积。当已获得叠后数据体时，往往希望由式(1)得到与数据相应的地下介质的参数，即反射系数序列，进而得到波阻抗，这就是所谓的叠后波阻抗反演[1，2]。这个过程中通常都假设反射系数序列是稀疏的，并通过井震标定确定子波，然后利用某种最优化算法求解设定好的目标函数，获得反射系数序列。

根据反射系数序列是稀疏的这一假设，可将反演问题表述为如下形式的基追踪问题

(2)

式中：d0表示观测地震道；‖·‖0表示L0范数。

由于求L0范数的极值问题很困难，且它对噪声十分敏感，因此在实际计算中常将式(2)近似表示为

(3)

式中‖·‖1表示L1范数。可证明在一定程度上式(3)的解就是式(2)的解，且式(3)比式(2)容易求解。其中，最常见的式(3)解法就是匹配追踪算法。

最初的匹配追踪算法由Mallat等[3]提出，是用于信号分解，随后出现了精度更高的正交匹配追踪，并广泛应用于地球物理领域[4-10]。同时也出现了基于匹配追踪算法的反演，如周东勇等[11]、Wen等[12]提出基于双极子分解匹配追踪算法的反演方法，刘晓晶等[13]提出基于正交匹配追踪算法的反演方法，Zhang等[14-18]和Deborah等[19]提出基于双极子分解基追踪算法的反演方法，印兴耀等[20]则在双极子分解的基础上用梯度投影算法求解反演问题，并开展实际应用[21，22]。与此同时，曹静杰[23]从贝叶斯稀疏反演理论出发，建立了非凸Lp范数正则化的盲反褶积模型，并利用最小二乘法进行反演； 梁东辉等[24]在对反射系数进行L0范数约束的基础上，利用迭代硬阈值法进行反射系数反演。

实际上，在稀疏假设下，对于地球物理反演问题，其表现形式应该是

(4)

式中：k表示稀疏度；p可等于1，也可等于2。当p=1时，式(4)转化为线性规划问题，即可用线性规划算法求解； 当p=2时，式(4)转化为二次规划问题，即可用通用二次规划算法求解。本文提出一种求解方法，依然基于式(2)的形式，即求L0范数的极值问题，但并不直接计算‖r‖0，而是将‖r‖0用一个平滑函数来近似。

2 基本原理

‖r‖0实际上定义为反射系数序列r中非零值的个数。也就是说，如果定义如下函数

(5)

那么‖r‖0就可表示为

(6)

其思路即是将式(6)中的阶跃函数s(ri)用连续的可微函数来近似。这里选用零均值高斯函数

(7)

则有

(8)

或者

(9)

则有

于是

‖r‖0≈n-fσ(r)

其中参数σ为控制精度。σ越小，近似程度越高，σ越大，则函数越平滑。由于n是一个定值，因此求‖r‖0的极小就转化为求fσ(r)的极大值。于是式(2)的基追踪问题就可表示为

(10)

该式可用许多通用的基于导数的优化算法[25](如最速下降法)求解。不同之处是在这个优化问题中要选择σ。由于|r|<1，根据式(9)，可选择一系列σ，其范围为0≤σ≤2，如[σ1,σ2,…,σK]=[2.0,1.9,…,0.1]。然后进行如下迭代计算：

一、固定σ=σ1，选择最小二乘解r0=(WTW+λI)-1·WTd0作为初始解，用最速下降法进行迭代求解：r←r+μ·fσ(r)， 并将其投影到可行域Wr=d0中， 投影方式为r←r-WT·(WWT)-1·(Wr-d0)，从而获得第一次迭代后的最优解r1；

二、令σ=σi，用上一次迭代得到的最优解ri作为初始解，用最速下降法进行下一次迭代，并将其投影到可行域，获得新的最优解ri；

三、令i=i+1，重复第二步；

四、最终获得最优解r=rK。

该过程中，无论是初始化还是在最速下降法中均需要选择一个参数λ，但是其意义不同。在初始化时，最小二乘解中λ的意义是使矩阵求逆稳定，因此这个时候的λ取值不宜太大。而在最速下降法中，λ可通过线性搜索来确定。为了简化计算过程也可指定一个值，但是这个定值需试验确定。本文通过模型试验，取λ为0.01，并且在实际资料计算时，将实际地震资料进行了归一化处理，没有改变λ的值。

对于σ的取值范围可以根据实际选择。在实际地震勘探中由于|r|<1，不妨分别取r=0.5和r=0.05，根据式(7)可得到近似函数的变化趋势，如图1所示。由图可见，即使是反射系数较大的情况，例如r=0.5，取σ最大为2即可使近似函数接近1；而对于反射系数较小的情况，例如r=0.05,σ最大为0.1即可使近似函数接近1。似乎只要取较大的σ(如2或4)，即可很好地近似所有情况，但是实际上其分辨率却不高，即不能区分反射系数大小。因此，在实际应用时，将地震资料归一化以后，在第一层循环中，从大到小、依次按照一定间隔取一系列σ。只要设定最大值、最小值及间隔步长即可。

图1 近似函数随σ的变化趋势

对每一个固定的σ，接下来就是用最优化算法求解式(10)。由于此时的目标函数是连续可微函数，因此，其优化算法就可有很多选择。本文用最速下降法，其终止条件可用最大迭代次数、目标函数的迭代终止误差或者目标函数变量的迭代终止误差进行控制。

3 模型试验

众所周知，如果式(1)中的噪声项为0，子波准确且为已知，则可容易且准确地反演反射系数。但实际情况是噪声项未知，子波虽然能提取但通常不十分准确，此时难以准确地反演反射系数，这对反演算法提出了严峻的考验。为此，本文首先从测井资料获得了准确的反射系数系列； 再用30Hz零相位Ricker子波与之褶积得到模拟地震道； 然后在模拟地震道中加入随机噪声作为观测数据(图2b，信噪比为233)，子波分别用准确的30Hz和28Hz零相位Ricker子波进行反演，反演结果如图3。

从图3a可见，本文提出的反演方法不仅可获得准确的绝对波阻抗值，还可获得较可靠的波阻抗变化趋势及界面位置。同时，弱随机噪声对本文反演方法的影响并不是很大，且实际资料中的随机噪声一般也较小，尤其是叠后或者局部叠加数据中的随机噪声较小。对比图3a与图3b可以看出，子波准确与否对反演结果的影响很大。虽然反演结果也能刻画波阻抗变化趋势以及界面形态，但界面位置及波阻抗值大小都偏离了真实值。这与目前常用的M商用反演软件的情况类似，即子波或者标定的准确与否决定了最终反演结果的可靠性。

用时间域波阻抗楔形模型(图4a)进行试验。时间采样间隔是1ms，共1500道。图4b是用35Hz零相位Ricker子波进行褶积得到的模拟数据。以此模拟数据作为输入，用本文提出的反演方法进行反演，反演时的子波采用准确的子波，即35Hz零相位Ricker子波。图4c是反演得到的波阻抗。对比图4a与图4c可见，当子波准确时，本文的反演方法能十分准确地获得反射界面的真实位置和波阻抗真实值，即使两个反射界面距离非常近，即所谓的薄层情况下，也能准确刻画反射界面的真实位置和波阻抗真实值。

图2 无噪合成道(a)和含噪合成道(b，信噪比为233)

图3 含噪、无噪合成道分别与准确(a)、不准确(b)子波的反演结果与真实值的对比

图4 楔形模型真实波阻抗剖面(a)及其合成地震数据(b)和反演波阻抗剖面(c)

4 实际地震资料处理

对图5实际剖面L1进行处理，图6和图7是反演后的反射系数剖面及相对波阻抗剖面。由于反演所用子波是30Hz零相位Ricker子波，而不是经过标定提取出来的，因此子波未必准确； 另外因未经标定，故波阻抗值也只能是相对的。

利用M商用软件也对图5地震资料进行了相对波阻抗(图8)计算。对比图8和图7可见，本反演方法获得的反射系数剖面和波阻抗剖面能够保持反射系数的稀疏性和横向连续性，从纵向上看，波阻抗的变化趋势以及反射界面还是能够得到清晰的反映，尤其在纵向上，其分辨率并不比商业软件反演的结果低。但是本文的反演方法本质上是稀疏脉冲波阻抗反演的一种实现，且在实施过程中是按道进行的，没有在横向上进行约束。即本文的反演仅仅利用了地震数据，而没有利用包含测井、地质等信息的模型，即没有利用地震信号带限外的高低频信息，尤其是低频信息，因此本文的反演方法得到的结果本质上是带限的相对波阻抗。同时，由于在具体实现时是根据式(4)按道进行的反演，没有规则化项，因此反演结果(波阻抗)仍存在“挂面条”现象(图7)。

图9是实际地震剖面L2，采样间隔为1ms(实际显示0.30～1.80s)。用本文方法反演后加入低频趋势，获得波阻抗剖面(图10)； 再将井旁道反演结果与井口计算的波阻抗(图11)进行对比。从图11可见，加入低频趋势以后，反演结果在总体趋势上与井口的波阻抗吻合很好，既可反映反射界面，也可获得绝对波阻抗，有利于储层描述。另外，从图10可见，反演结果从纵向上反映出波阻抗变化趋势，但横向上不平滑。为了进一步提高反演质量，在做好标定与子波提取的基础上，还需在反演过程中加入低频模型或者横向约束。

图5 实际过井地震剖面L1

图6 本文方法反演的L1反射系数剖面

图7 本文方法反演的L1相对波阻抗剖面

图8 利用M商用软件反演的L1相对波阻抗剖面

图9 实际地震剖面L2

图10 反演后加入低频趋势的L2波阻抗剖面

图11 反演的井旁波阻抗(蓝色)与井口波阻抗(绿色)的对比

5 结论

针对地下地层反射系数是稀疏的这一假设，本文将地震波阻抗反演视为L0范数约束下的最优化问题，并将L0范数用一个平滑函数近似。然后将L0范数约束下的基追踪问题转化为无约束最优化问题，并利用基于导数的最优化方法求解无约束优化问题，求得反射系数，进而计算得到波阻抗。由于目标函数中包含了L0范数约束项，因此反演的反射系数是稀疏的。同时将L0范数用一个平滑函数近似，就可利用基于导数的优化方法，从而使得反演快速收敛。

模型试验表明，当子波准确时，即使有噪声，该方法也可十分准确地获得波阻抗值和反射界面。由于实际地震资料是带限的，因此实际资料计算结果是相对波阻抗，但波阻抗变化趋势及反射界面可清晰刻画，并且在加入低频信息以后波阻抗真值可得到恢复。为了进一步提升反演效果，还需要在具体实现时加入横向约束。