基于Lp拟范数稀疏约束和交替方向乘子算法的波阻抗反演

2022-10-10张雨强文晓涛

石油物探 2022年5期

张雨强,文晓涛,吴昊,刘军,刘炀

(1.成都理工大学地球物理学院,四川成都 610059;2.成都理工大学油气藏地质及开发工程重点实验室,四川成都 610059;3.中国地质大学构造与油气资源教育部重点实验室,湖北武汉 430074)

地震数据中含有丰富的关于储层地层学、沉积学和岩石性质的信息,地震反演是地质解释中挖掘上述信息的有效途径[1-2]。根据反演所利用的基础数据,可将反演分为叠后反演和叠前反演。相对而言,叠前反演信息量更丰富,反演结果的分辨率更高。但叠后反演也有其自身的优势,如可以相对快速地将地震波形信息转换为具有地质意义的波阻抗信息等,指导储层的识别与刻画[3]。因此,从反演技术的发展趋势来看,叠后地震反演技术在油气勘探中仍然发挥着重要的作用。

在研究反演的过程中,如何提高反演的精度一直都是国内外学者重点关注的问题,比较常用的方法是增加先验信息。基于L2范数约束的Tikhonov正则化反演方法,一定程度上提高了反演的精度。例如,广义线性反演[4]和稀疏尖峰反演[5]。但是,使用L2范数进行正则化约束会导致反演的结果过于光滑,反射界面不明显。鉴于反射系数可看作由大量零值组成的稀疏信号,可利用该稀疏性构造稀疏正则项,进一步对反演结果进行约束[6]。

稀疏正则化的构建,通常是采用L1范数进行正则化约束。在地球物理领域,较早使用L1范数正则化的是在反褶积[7]和数据重构[8]等方面。但由于当时没有合适的算法对L1范数约束的优化问题进行准确求解,使得L1范数正则化的反演算法未能得到进一步的推广。近十多年来,压缩感知理论快速发展,使得L1范数正则化再次引起了学者们的关注。

ZHANG等[9-10]使用基追踪算法求解L1范数约束问题,通过对叠后地震数据反演得到反射系数,为进一步增强横向连续性,引入“Z型”空间导数进行多道反演。LIU等[11]使用L1范数正则化反演方法进行波阻抗反演,并将该方法与阻尼最小二乘阻抗反演方法对比,对比结果显示L1范数正则化反演方法具有更高的反演精度。张繁昌等[12]对基追踪反演算法进行改进,将目标函数中保真项的L2范数约束替换为L1范数约束,并使用对偶对数障碍规划算法进行反演,反演精度相较于传统基追踪算法得到进一步提高。石战战等[13]为了增加在反演中噪声的鲁棒性,在目标函数中引入L1范数拟合项,提出一种L1-L1范数稀疏约束地震反演方法。ZHOU等[14]提出了一种多道基追踪叠后反演算法,该算法在L1范数约束的基础上引入马尔可夫过程描述相邻地震道间的关系,提高了对薄层的识别能力。WU等[15]将L1范数稀疏约束的拓展形式交叠组稀疏用于波阻抗反演,通过交叠组稀疏的结构特性消除了局部异常值,进一步提高反演精度。WANG等[16]提出了一种数据驱动的地震波阻抗反演方法,该方法利用L1范数对地震空间连续信息进行约束并参与反演,进一步提高反演精度。

上述算法大多是基于L1范数对地震数据的稀疏信息进行挖掘,但是L1范数只是L0范数的凸松弛函数,得到的稀疏先验信息有限,于是直接使用L0范数进行反演[17-19]。L0范数虽然是最稀疏的范数,但由于目前主流的求解算法是贪婪迭代策略方法[20],这种方法容易在迭代时陷入局部极值。最近,有学者提出Lp拟范数约束,也就是范数的p次幂(0

选择比L1范数更为稀疏的Lp拟范数(0

1 基本原理

1.1 正演模型

地震信号可以由反射系数与地震子波褶积表示:

S=w*R+N

(1)

式中:S表示地震记录;R表示反射系数;w表示子波;N为随机噪声;“*”代表卷积运算符号。将卷积运算转化为矩阵相乘,得到公式:

S=WR+N

(2)

W的定义为:

(3)

当反射系数较小时,反射系数与波阻抗之间存在以下关系:

(4)

式中:Z为波阻抗;i为采样点序号;j为道数序号。假设L=lnZ,反射系数与波阻抗之间的关系表示为:

R=DL

(5)

式中:D为差分矩阵。D定义为:

(6)

结合以上公式,则公式可以表示为:

S=WDL+N

(7)

1.2 目标函数

在上述正演模型基础上,传统的目标函数可以表示为

(8)

为了得到更高精度反演结果,采用Lp拟范数对反射系数进行稀疏约束。Lp拟范数的具体表达式为:

(9)

式中:0

(10)

(11)

图1 反射系数优化问题a L1范数; b Lp范数

结合(8)式,并在目标函数中加入初始模型约束,得到新的目标函数为:

(12)

式中:L0为初始模型;μ为初始模型约束参数;λ为正则项的约束参数。

1.3 反演方法

在上述目标函数中,不仅存在Frobenius范数,而且还存在Lp拟范数。因此,无法对该问题直接求解。在此引入交替方向乘子算法,将目标函数中的不同范数约束转化为可以直接进行求解的多个子函数,交替进行求解。具体过程为:首先,在(12)式中引入拉格朗日乘子项,用R来代替DL,即

(13)

引入对偶项将上述目标函数转化为无约束优化问题:

(14)

式中:C表示R的对偶项;η为拉格朗日乘子的正则化参数项。

基于交替方向乘子算法(ADMM)流程,将目标函数分解为分别与L,R,C有关的子函数。其中,与L有关的目标函数为:

(15)

(15)式为常规的优化问题,可直接采用凸优化算法的求解方法直接得到L的更新方式,即

Li+1=(DTWTWD+μE+ηDTD)-1·
(DTWTS+μL0+ηDTRi-ηDTCi)

(16)

式中:E表示单位矩阵。在得到更新之后的L时,则可以对R的子函数进行求解,即

(17)

由(17)式中可知,该目标函数为一个基于Lp拟范数的优化问题,在此引入软阈值收缩算法进行求解[25],得到结果为:

(18)

其中,sign函数表示为:

(19)

而关于对偶项C的子函数可以表示为:

(20)

该问题也属于凸优化问题,可利用其求梯度进行求解:

Ci+1=Ci+DLi+1-Ri+1

(21)

综合上述内容,得到反演技术流程框架:

算法:基于Lp稀疏约束和交替方向乘子算法的波阻抗反演输入:地震记录S,初始模型L0,地震子波w,参数项μ,λ,η,p,差分矩阵D,收敛误差ε,道数trace输出:地震波阻抗Z初始化:i=0,t=0,L=L0,R=0,C=01. While t

2 模型测试与实际应用

2.1 模型测试

采用Marmousi2模型的部分数据构成理论模型(图2),对反演方法进行测试。图2a为理论阻抗模型的二维剖面,该模型共有500道地震数据,每道共有450个采样点,在采样点200至250,道数250至450之间存在透镜体形状的含气储层。利用该波阻抗模型,基于公式(4)得到反射系数(图2b),然后将得到的反射系数与30Hz的Ricker子波进行褶积得到地震记录(图2c)。图2d为使用高斯滤波器对原始波阻抗进行低通滤波得到的初始模型。

图2 理论模型a 理论阻抗模型; b 反射系数; c 合成地震记录; d 初始模型

为了验证新方法的可行性及优势,使用基于L1范数约束的基追踪反演算法与Lp拟范数约束的反演算法对上述模型进行反演,并比较两者与理论模型的差异。为了进一步验证新方法的抗噪性,对合成地震记录分别加入20%和50%的高斯随机噪声后,再次反演。

首先,选择L1范数约束与Lp拟范数约束反演的第350道结果进行对比。不同噪声下单道反演结果对比如图3所示,黑色线条为理论阻抗模型;绿色线条为初始模型;红色线条为基于L1范数约束的单道反演结果;蓝色线条为新方法的单道反演结果。由图3可见,在不同噪声情况下的单道反演结果对比中,采用新方法得到的反演结果更接近理论模型,该现象在红色和蓝色椭圆中较为明显。图4,图5和图6 分别展示了不同噪声情况下的反演结果,可以看出,Lp拟范数约束反演结果的横向连续性优于L1范数约束反演结果,在纵向展布上的值也更接近理论模型。为了定量比较,计算了两种反演结果与理论模型的信噪比(signal noise ratio,SNR)与均方误差(root mean square error,RMSE),具体公式为:

图3 不同噪声下单道反演结果对比a 无噪反演结果单道对比; b 含20%噪声单道反演结果对比; c 含50%噪声单道反演结果对比

图4 无噪声反演结果a L1范数约束反演结果; b Lp拟范数约束反演结果

图5 含20%噪声反演结果a L1范数约束反演结果; b Lp拟范数约束反演结果

图6 含50%噪声反演结果a L1范数约束反演结果; b Lp拟范数约束反演结果

(22)

(23)

表1和表2分别为两种约束条件下反演结果的信噪比与均方误差分析结果。由表1可见,新方法得到的反演结果的信噪比大于L1范数得到的反演结果,而其均方误差小于L1范数的反演结果。由此证明了使用新方法得到的反演结果优于L1范数的反演结果。

表1 基于理论模型的信噪比分析结果

表2 基于理论模型的均方误差分析结果

2.2 实际应用

研究区位于四川盆地南缘,目的层为志留系龙马溪组,岩性为页岩。两口钻井(A井和B井)均钻遇优质页岩气储层,但储层厚度较小(小于10m),常规的阻抗反演方法由于受到精度和分辨率的限制,很难对页岩气储层进行有效识别。为此,选择该区域作为试验区,验证新方法的有效性。图7为过A井和B井的连井剖面。本次反演中,在地震层位信息的指导下,选择A井波阻抗数据进行外推插值并做低通滤波处理得到低频模型,选择未进行反演的B井对反演结果进行验证。分别使用基于L1范数约束的反演方法与新方法进行反演,得到的结果如图8和图9所示。

图7 工区地震剖面

由图8和图9可以看到,两种方法反演的纵波阻抗与A、B两井测井曲线吻合较好,表明两种方法在野外资料应用中可行且有效。与基于L1范数约束反演方法的纵波阻抗剖面相比,新方法能够反演出更

图8 基于L1范数约束反演结果

多地层信息(图8和图9中的黄色圆圈)。同时,新方法获得的阻抗信息具有更好的横向连续性(图8和图9 中的红色圆圈)。为了进一步说明新方法的优越性,利用两种方法的反演结果重新合成地震剖面(图10 和图11),再分别与实际地震剖面进行残差分析(图12和图13)。由图12和图13可见,新方法的合成地震剖面与实际地震剖面残差更小,即新方法反演结果与实际数据更吻合。