石油地球物理勘探  2024, Vol. 59 Issue (2): 230-237  DOI: 10.13810/j.cnki.issn.1000-7210.2024.02.005
0
文章快速检索     高级检索

引用本文 

张天悦, 林凯, 文晓涛, 赵炼, 张雨强, 雷扬. 应用Lp拟范数稀疏约束的纵横波速比直接反演. 石油地球物理勘探, 2024, 59(2): 230-237. DOI: 10.13810/j.cnki.issn.1000-7210.2024.02.005.
ZHANG Tianyue, LIN Kai, WEN Xiaotao, ZHAO Lian, ZHANG Yuqiang, LEI Yang. Direct inversion of P-wave to S-wave velocity ratio by Lp quasi-norm sparse constraints. Oil Geophysical Prospecting, 2024, 59(2): 230-237. DOI: 10.13810/j.cnki.issn.1000-7210.2024.02.005.

本项研究受国家自然科学青年基金项目“含流体孔隙介质倾斜入射地震平面波一阶近似解析理论研究”(42104131)和四川省自然科学基金青年基金项目“四川盆地深部碳酸盐岩储层含油气预测方法研究”(2022NSFSC1140)联合资助

作者简介

张天悦  硕士研究生,2000年生;2021年获河南理工大学地球信息科学与技术专业学士学位。现于成都理工大学攻读地质工程专业硕士学位,主要从事油气储层解释方面的学习和研究

林凯, 四川省成都市成华区二仙桥东三路1号成都理工大学地球物理学院,610059。Email:linkai02102@163.com

文章历史

本文于2023年6月9日收到,最终修改稿于同年12月28日收到
应用Lp拟范数稀疏约束的纵横波速比直接反演
张天悦1,2 , 林凯1,2 , 文晓涛1,2 , 赵炼1,2 , 张雨强1,2 , 雷扬3     
1. 成都理工大学油气藏地质及开发工程重点实验室, 四川成都 610059;
2. 成都理工大学地球物理学院, 四川成都 610059;
3. 中国石化石油工程地球物理公司南方分公司, 四川成都 610041
摘要:纵横波速比(vP/vS)是识别气藏、描述储层特征和判别岩性的重要解释工具。目前主要是通过反射系数近似方程反演得到纵、横波速度,再进一步计算纵横波速比,但是这种间接计算方法会产生累积误差。为了直接从叠前地震数据反演纵横波速比,文中提出了一种新的广义弹性阻抗方程,再进一步推导出一个与纵横波速比、纵波速度、密度相关的纵波反射系数近似方程。为了得到精度较高的反演结果,基于推导出的反射系数近似方程,提出一种基于Lp拟范数稀疏约束的叠前地震反演方法,并通过交替方向乘子算法求解。将提出的直接反演方法应用于理论模型和实际数据,并与间接反演方法相对比,结果表明该直接反演方法的反演结果精度较高,对含气储层的边界刻画更清晰。
关键词反演    纵横波速比    广义弹性阻抗    Lp拟范数    交替方向乘子算法    
Direct inversion of P-wave to S-wave velocity ratio by Lp quasi-norm sparse constraints
ZHANG Tianyue1,2 , LIN Kai1,2 , WEN Xiaotao1,2 , ZHAO Lian1,2 , ZHANG Yuqiang1,2 , LEI Yang3     
1. State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation, Chengdu University of Technology, Chengdu, Sichuan 610059, China;
2. College of Geophysics, Chengdu University of Technology, Chengdu, Sichuan 610059, China;
3. South Branch of SINOPEC Petroleum Engineering Geophysical Co., Ltd, Chengdu, Sichuan 610041, China
Abstract: The P-wave to S-wave velocity ratio(vP/vS) is a vital tool for gas reservoir identification, reservoir characterization, and lithology recognition. At present, the P-wave velocity and S-wave velocity are mainly obtained through the inversion of the reflection coefficient approximate equation, and then the vP/vS is calculated. However, this indirect calculation method creates a cumulative error. To obtain the vP/vS directly from pre-stack seismic data, this paper proposes a new generalized elastic impedance equation and further derives an approximate equation of the P-wave reflection coefficient, which is linked to the vP/vS, P-wave velocity, and density. To achieve high-precision inversion results, this paper proposes a prestack seismic inversion method based on the sparse constraint of the Lp quasi-norm utilizing the derived approximation equation of reflection coefficient, which is solved by the alternating direction multiplier algorithm. The proposed direct inversion method is applied to theoretical models and practical data and compared with the indirect inversion method. The results demonstrate that the direct inversion method exhibits higher inversion accuracy and clearer boundary characterization of gas-bearing reservoirs.
Keywords: inversion    P-wave to S-wave velocity ratio(vP/vS)    generalized elastic impedance    Lp quasi-norm    alternating direction multiplier algorithm    
0 引言

地震反演是地震勘探中重要的解释方法,它可以利用地震资料解释地下弹性参数,进一步进行储层预测,其中弹性参数包括纵波速度、横波速度和密度等。地震反演可以分为叠后反演和叠前反演,叠后反演具有简便、高效的特点,在储层预测中得到了广泛的应用,但是只能得到少数类型的参数,不能得到反映物理性质和流体特性的参数。与叠后反演相比,叠前地震反演利用地震振幅和炮检距的信息,为研究岩性和储层变化提供了充足有效的数据量[1-2],因此叠前地震反演是油气储层定量解释的一种重要手段[3-5]

纵横波速度比(vP/vS)是表征储层和区分岩性的重要参数。各种实验测量表明,vP/vS在岩性识别和异常孔隙压力预测方面具有广泛的应用[6]。特别是,由于该比值与泊松比直接相关[7-8],对流体性质高度敏感[9-10],已广泛用为有效油气鉴别器。

鉴于vP/vS的重要性,已开发了多种估计方法。较常用的是基于振幅随炮检距的变化(AVO)理论,使用AVO近似式从地震数据中反演纵、横波速度或模量或阻抗,以间接方式计算vP/vS[11]。这种间接的弹性参数转换会引起累积误差。因此,本文在前人研究的基础上提出一种从地震数据中直接获得vP/vS的反演方法,避免了间接方法的不稳定性。

Connolly[12]首次提出了弹性阻抗(EI)的概念。EI是在改进的Zoeppritz方程基础上推导而来,包含了纵波速度、横波速度和密度等信息,较常规波阻抗能反映更多的信息。然而,常规的EI方程存在一个不良特征,数值随入射角有较大尺度的变化。Whitcombe[13]对EI方程进行了标准化,使EI在一定角度范围内可以直观对比。Whitcombe等[14]对常规的EI定义进行了扩展,使其能够满足对流体和岩性的预测。为了平衡反演精度和稳定性,EI反演中嵌入的vP/vS通常是一个恒定值。马劲风[15]考虑了vP/vS随深度的变化,导出了递推形式的广义弹性阻抗(GEI)表达式。根据马劲风的研究,GEI的精度高于EI,但是GEI表达式较复杂,限制了反演的适用性。本文对GEI方程进行了改进,导出了关于vP/vS的反射系数近似方程,将vP/vS作为一个模型参数,可直接由叠前地震数据反演。

由于地球物理反演问题的多解性,导致反演结果精度较低。较常用的方法是通过正则化方法对反演问题添加先验信息。印兴耀等[16]在基于L1范数约束的基追踪反演算法基础上,加入低频模型约束,对波阻抗进行反演,发现该方法能够消除子波旁瓣对反演结果的影响。刘晓晶等[17]基于反射系数奇偶分解实现了关于流体因子在深部储层的基追踪EI反演方法。She等[18]在一阶纵向差分的基础上,提出了基于L1范数约束的高阶差分叠前地震反演方法,可以有效应对地质结构更复杂的情况。文晓涛等[19]在L1范数稀疏约束的基础上,加入低频模型残差约束项,提出一种低频稀疏双约束的波阻抗反演方法,有效提高反演结果的分辨率。虽然L1范数正则化方法提高了反演的稀疏性,但是对稀疏先验信息挖掘程度不够。Woodworth等[20]指出Lp拟范数比L1范数更接近L0范数。Chen等[21]基于Lp拟范数的稀疏性约束,提出了一种新的时频分析方法。张雨强等[22]将Lp拟范数稀疏约束项引入叠后反演,得到了精度更高的反演结果。综合来看,与L1范数相比,Lp拟范数是一种更为稀疏的L0范数近似。

基于以上分析,本文提出了一种基于Lp拟范数稀疏约束的叠前vP/vS直接反演方法。

1 基本原理 1.1 改进的GEI公式

EI由于含有地震AVO信息,在地震解释和流体识别中得到广泛的应用。但常规弹性阻抗中的vP/vS通常设为常数,未考虑到随深度的变化。马劲风[15]将GEI表示为

$ \mathrm{G}\mathrm{E}{\mathrm{I}}_{i}=\frac{{\rho }_{i}{{v}_{\mathrm{P}}}_{i}}{\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{i}}{\left(1-\frac{\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}{\theta }_{i}}{{r}_{i}^{2}}\right)}^{4} $ (1)

式中:ri=vPi/vSi,为第i层纵横波速比;ρ为岩石密度;θ为P波入射角。式(1)是关于vP/vS、纵波速度和密度的表达式。P波反射系数可以表示为

$ R\left(\theta \right)=\frac{1}{2}\frac{\mathrm{\Delta }\mathrm{G}\mathrm{E}\mathrm{I}}{\mathrm{G}\mathrm{E}\mathrm{I}}\approx \frac{1}{2}\mathrm{l}\mathrm{n}\frac{\mathrm{G}\mathrm{E}{\mathrm{I}}_{i+1}}{\mathrm{G}\mathrm{E}{\mathrm{I}}_{i}} $ (2)

对式(1)进行对数运算,简化为

$ \mathrm{l}\mathrm{n}\left(\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{i}\mathrm{G}\mathrm{E}{\mathrm{I}}_{i}\right)=\mathrm{l}\mathrm{n}{\rho }_{i}+\mathrm{l}\mathrm{n}{{v}_{\mathrm{P}}}_{i}+4\mathrm{l}\mathrm{n}\left(1-\frac{\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}{\theta }_{i}}{{r}_{i}^{2}}\right) $ (3)

$ \theta $小于30°时,$ 1/{r}_{i}^{2} $最大值为0.5,$ \mathrm{s}\mathrm{i}{\mathrm{n}}^{2}{\theta }_{i} $最大值为0.25,所以$ \mathrm{s}\mathrm{i}{\mathrm{n}}^{2}{\theta }_{i}/{r}_{i}^{2} $的最大值为0.125,是一个较小的值,将$ \mathrm{l}\mathrm{n}(1-\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}{\theta }_{i}/{r}_{i}^{2}) $用泰勒公式展开,有

$ \mathrm{l}\mathrm{n}\left(1-\frac{\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}{\theta }_{i}}{{r}_{i}^{2}}\right)\approx -\frac{\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}{\theta }_{i}}{{r}_{i}^{2}} $ (4)

将式(4)代入式(3),可得到新的GEI表达式

$ \mathrm{G}\mathrm{E}{\mathrm{I}}_{i}=\frac{{\rho }_{i}{{v}_{\mathrm{P}}}_{i}}{\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{i}}{\mathrm{e}}^{\left(-\frac{4\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}{\theta }_{i}}{{r}_{i}^{2}}\right)} $ (5)

应用Goodway等[23]根据实际资料给出的含气砂岩模型(表 1)进一步说明改进的GEI公式准确性。改进前、后计算的反射系数与精确Zoeppritz方程及Aki-Richards近似方程[24]计算的反射系数对比如图 1所示。该模型包含两个反射系数:上两层之间的负反射系数和下两层之间的正反射系数。由图 1可以看出,改进后计算的反射系数比改进前更接近Zoeppritz方程,因为改进前忽略密度反射率而产生了较大的误差;在入射角小于25°时,改进后的反射系数计算结果与Zoeppritz方程及Aki-Richards方程基本一致。

表 1 三层含气砂岩模型参数

图 1 四种方法计算的反射系数对比 (a)负反射系数;(b)正反射系数
1.2 反射系数公式

在上、下介质弹性参数差异小的情况下,使用改进的GEI建立线性化的纵波反射系数方程。对式(5)进行微分,可得

$ \begin{array}{l}\mathrm{\Delta }\mathrm{G}\mathrm{E}\mathrm{I}=\left[\frac{{v}_{\mathrm{P}}\mathrm{\Delta }\rho +\mathrm{s}\mathrm{e}{\mathrm{c}}^{2}\theta \mathrm{t}\mathrm{a}{\mathrm{n}}^{2}\theta \rho \mathrm{\Delta }{v}_{\mathrm{P}}}{\mathrm{c}\mathrm{o}\mathrm{s}\theta }-8\mathrm{s}\mathrm{i}\mathrm{n}\theta \times \right.\\ \begin{array}{cc}& \left.\mathrm{t}\mathrm{a}\mathrm{n}\theta {\left(\frac{1}{r}\right)}^{2}\rho \mathrm{\Delta }{v}_{\mathrm{P}}-\frac{8\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\theta }{\mathrm{c}\mathrm{o}\mathrm{s}\theta }\rho {v}_{\mathrm{S}}\mathrm{\Delta }\left(\frac{1}{r}\right)\right]{\mathrm{e}}^{\left(\frac{-4\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\theta }{{r}^{2}}\right)}\end{array}\end{array} $ (6)

式中$ \mathrm{\Delta }\left(1/r\right)=-\mathrm{\Delta }r/{r}^{2} $。结合式(2)、式(6)可得纵波反射系数的线性公式

$ \begin{array}{l}R\left(\theta \right)=4{\left(\frac{{v}_{\mathrm{S}}}{{v}_{\mathrm{P}}}\right)}^{2}\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\theta \frac{\mathrm{\Delta }r}{r}+\\ \begin{array}{cc}& \left[\frac{1}{2}\mathrm{s}\mathrm{e}{\mathrm{c}}^{2}\theta -4{\left(\frac{{v}_{\mathrm{S}}}{{v}_{\mathrm{P}}}\right)}^{2}\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\theta \right]\frac{\mathrm{\Delta }{v}_{\mathrm{P}}}{{v}_{\mathrm{P}}}+\frac{1}{2}\frac{\mathrm{\Delta }\rho }{\rho }\end{array}\end{array} $ (7)

式中$ \mathrm{\Delta }r/r $$ \mathrm{\Delta }{v}_{\mathrm{P}}/{v}_{\mathrm{P}} $$ \mathrm{\Delta }\rho /\rho $分别为纵横波速比、纵波速度和密度的反射系数,具体分别为

$ \frac{\mathrm{\Delta }r}{r}=\frac{2\left({r}_{i+1}-{r}_{i}\right)}{{r}_{i+1}+{r}_{i}} $ (8)
$ \frac{\mathrm{\Delta }{v}_{\mathrm{P}}}{{v}_{\mathrm{P}}}=\frac{2\left({{v}_{\mathrm{P}}}_{i+1}-{{v}_{\mathrm{P}}}_{i}\right)}{{{v}_{\mathrm{P}}}_{i+1}+{{v}_{\mathrm{P}}}_{i}} $ (9)
$ \frac{\mathrm{\Delta }\rho }{\rho }=\frac{2\left({\rho }_{i+1}-{\rho }_{i}\right)}{{\rho }_{i+1}+{\rho }_{i}} $ (10)
1.3 反演算法

式(7)可以写为矩阵形式

$ \boldsymbol{R}\left(\theta \right)=a\left(\theta \right)\frac{\mathrm{\Delta }\boldsymbol{r}}{\boldsymbol{r}}+b\left(\theta \right)\frac{\mathrm{\Delta }{\boldsymbol{v}}_{\mathrm{P}}}{{\boldsymbol{v}}_{\mathrm{P}}}+c\left(\theta \right)\frac{\mathrm{\Delta }\boldsymbol{\rho }}{\boldsymbol{\rho }} $ (11)

式中$ a\left(\theta \right) $$ b\left(\theta \right) $$ c\left(\theta \right) $为纵横波速比、纵波速度和密度反射系数的参数项,表示为

$ \left\{\begin{array}{l}a\left(\theta \right)=4{\left(\frac{{v}_{\mathrm{S}}}{{v}_{\mathrm{P}}}\right)}^{2}\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\theta \\ b\left(\theta \right)=\frac{1}{2}\mathrm{s}\mathrm{e}{\mathrm{c}}^{2}\theta -4{\left(\frac{{v}_{\mathrm{S}}}{{v}_{\mathrm{P}}}\right)}^{2}\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\theta \\ c\left(\theta \right)=\frac{1}{2}\end{array}\right. $ (12)

$ \mathrm{\Delta }\boldsymbol{r}/\boldsymbol{r} $的元素$ \mathrm{\Delta }{r}_{i, j}/{r}_{i, j} $表示为

$ \frac{\mathrm{\Delta }{r}_{i, j}}{{r}_{i, j}}=\mathrm{\Delta }\mathrm{l}\mathrm{n}{r}_{i, j}=\mathrm{l}\mathrm{n}{r}_{i+1, j}-\mathrm{l}\mathrm{n}{r}_{i, j} $ (13)

式中j为道号。于是

$ \frac{\mathrm{\Delta }\boldsymbol{r}}{\boldsymbol{r}}=\boldsymbol{D}{\boldsymbol{L}}_{r} $ (14)

式中:$ {\boldsymbol{L}}_{r}=\mathrm{l}\mathrm{n}\boldsymbol{r} $,表示对矩阵r中每一个元素都取对数;D为差分矩阵,即

$ \boldsymbol{D}=\left(\begin{array}{ccccc}-1& 1& & & \\ & -1& 1& & \\ & & \ddots & \ddots & \\ & & & -1& 1\end{array}\right) $ (15)

类似地,有

$ \frac{\mathrm{\Delta }{\boldsymbol{v}}_{\mathrm{P}}}{{\boldsymbol{v}}_{\mathrm{P}}}=\boldsymbol{D}{\boldsymbol{L}}_{{v}_{\mathrm{P}}} $ (16)
$ \frac{\mathrm{\Delta }\boldsymbol{\rho }}{\boldsymbol{\rho }}=\boldsymbol{D}{\boldsymbol{L}}_{\rho } $ (17)

则式(11)可以表示为

$ \boldsymbol{R}\left(\theta \right)=a\left(\theta \right)\boldsymbol{D}{\boldsymbol{L}}_{r}+b\left(\theta \right)\boldsymbol{D}{\boldsymbol{L}}_{{v}_{\mathrm{P}}}+c\left(\theta \right)\boldsymbol{D}{\boldsymbol{L}}_{\rho } $ (18)

式中:$ {\boldsymbol{L}}_{{v}_{\mathrm{P}}}=\mathrm{l}\mathrm{n}{\boldsymbol{v}}_{\mathrm{P}} $$ {\boldsymbol{L}}_{\rho }=\mathrm{l}\mathrm{n}\boldsymbol{\rho } $。根据褶积模型可以得到地震记录和反射系数的关系表达式

$ \boldsymbol{S}\left(\theta \right)=\boldsymbol{W}\boldsymbol{R}\left(\theta \right) $ (19)

式中$ \boldsymbol{W} $表示子波矩阵。如果有N个入射角,S可写为

$ \boldsymbol{S}={\left[\begin{array}{cccc}\boldsymbol{S}\left({\theta }_{1}\right)& \boldsymbol{S}\left({\theta }_{2}\right)& \cdots & \boldsymbol{S}\left({\theta }_{N}\right)\end{array}\right]}^{\mathrm{T}} $ (20)

式(19)可进一步表示为

$ \boldsymbol{S}=\boldsymbol{A}\boldsymbol{L} $ (21)

式中:$ \boldsymbol{L}={\left[\begin{array}{ccc}{{\boldsymbol{L}}_{r}}^{\mathrm{T}}& {{\boldsymbol{L}}_{{v}_{\mathrm{P}}}}^{\mathrm{T}}& {{\boldsymbol{L}}_{\rho }}^{\mathrm{T}}\end{array}\right]}^{\mathrm{T}} $$ \boldsymbol{A} $可以表示为

$ \boldsymbol{A}=\left[\begin{array}{ccc}a\left({\theta }_{1}\right)\boldsymbol{W}\boldsymbol{D}& b\left({\theta }_{1}\right)\boldsymbol{W}\boldsymbol{D}& c\left({\theta }_{1}\right)\boldsymbol{W}\boldsymbol{D}\\ a\left({\theta }_{2}\right)\boldsymbol{W}\boldsymbol{D}& b\left({\theta }_{2}\right)\boldsymbol{W}\boldsymbol{D}& c\left({\theta }_{2}\right)\boldsymbol{W}\boldsymbol{D}\\ & ⋮& \\ a\left({\theta }_{N}\right)\boldsymbol{W}\boldsymbol{D}& b\left({\theta }_{N}\right)\boldsymbol{W}\boldsymbol{D}& c\left({\theta }_{N}\right)\boldsymbol{W}\boldsymbol{D}\end{array}\right] $ (22)

反演的目标函数可表示为

$ J\left(\boldsymbol{L}\right)=\underset{\boldsymbol{L}}{\mathrm{m}\mathrm{i}\mathrm{n}}{‖\boldsymbol{S}-\boldsymbol{A}\boldsymbol{L}‖}_{2}^{2} $ (23)

为了提高反演的稳定性,加入了纵横波速比、纵波速度、密度的初始模型约束项,同时采用Lp拟范数构造反射稀疏约束项

$ \begin{array}{l}J\left(\boldsymbol{L}\right)\mathrm{ }=\underset{\boldsymbol{L}}{\mathrm{m}\mathrm{i}\mathrm{n}}\left({‖\boldsymbol{S}-\boldsymbol{A}\boldsymbol{L}‖}_{2}^{2}+\alpha {‖{\boldsymbol{L}}_{r}-\mathrm{l}\mathrm{n}{\boldsymbol{r}}_{0}‖}_{2}^{2}+\right.\\ \left.\mu {‖{\boldsymbol{A}}_{1}\boldsymbol{L}‖}_{p}^{p}+\beta {‖{\boldsymbol{L}}_{{v}_{\mathrm{P}}}-\mathrm{l}\mathrm{n}{\boldsymbol{v}}_{\mathrm{P}0}‖}_{2}^{2}+\gamma {‖{\boldsymbol{L}}_{\rho }-\mathrm{l}\mathrm{n}{\boldsymbol{\rho }}_{0}‖}_{2}^{2}\right)\end{array} $ (24)

式中:αβγ分别为纵横波速比、纵波速度和密度的初始模型约束项参数;μ为稀疏约束项参数;r0vP0ρ0分别为纵横波速比、纵波速度和密度的初始模型;$ {‖·‖}_{p}^{p} $表示Lp拟范数,其中p的取值范围为0~1,p值越小,反演结果会越稀疏,可根据实际情况取最优值;矩阵$ {\boldsymbol{A}}_{1} $满足反射系数$ \boldsymbol{R}={\boldsymbol{A}}_{1}\boldsymbol{L} $,可表示为

$ {\boldsymbol{A}}_{1}=\left[\begin{array}{ccc}a\left({\theta }_{1}\right)\boldsymbol{D}& b\left({\theta }_{1}\right)\boldsymbol{D}& c\left({\theta }_{1}\right)\boldsymbol{D}\\ a\left({\theta }_{2}\right)\boldsymbol{D}& b\left({\theta }_{2}\right)\boldsymbol{D}& c\left({\theta }_{2}\right)\boldsymbol{D}\\ & ⋮& \\ a\left({\theta }_{N}\right)\boldsymbol{D}& b\left({\theta }_{N}\right)\mathbf{W}\mathbf{D}& c\left({\theta }_{N}\right)\boldsymbol{D}\end{array}\right] $ (25)

因为目标函数包括了Lp$ {\mathrm{L}}_{2} $范数,因此采用交替方向乘子法求解[25]。引入拉格朗日项$ {\boldsymbol{L}}_{\mathrm{y}} $代替Lp拟范数的运算,将式(24)转换为可以直接求解的优化问题

$ \begin{array}{l}J\left(\boldsymbol{L}\right)=\underset{\boldsymbol{L}, {\boldsymbol{L}}_{y}}{\mathrm{m}\mathrm{i}\mathrm{n}}\left({‖\boldsymbol{S}-\boldsymbol{A}\boldsymbol{L}‖}_{2}^{2}+\alpha {‖{\boldsymbol{L}}_{r}-\mathrm{l}\mathrm{n}{\boldsymbol{r}}_{0}‖}_{2}^{2}+\right.\\ \mu {‖{\boldsymbol{L}}_{\mathrm{y}}‖}_{p}^{p}+\left.\beta {‖{\boldsymbol{L}}_{{v}_{\mathrm{P}}}-\mathrm{l}\mathrm{n}{\boldsymbol{v}}_{\mathrm{P}0}‖}_{2}^{2}+\gamma {‖{\boldsymbol{L}}_{\rho }-\mathrm{l}\mathrm{n}{\boldsymbol{\rho }}_{0}‖}_{2}^{2}\right)\\ \mathrm{s}\mathrm{t}{\boldsymbol{L}}_{\mathrm{y}}={\boldsymbol{A}}_{1}\boldsymbol{L}\end{array} $ (26)

再引入对偶项$ \boldsymbol{C} $,将式(26)转换成无约束的最优化问题

$ \begin{array}{l}J\left(\boldsymbol{L}\right)=\underset{\boldsymbol{L}, {\boldsymbol{L}}_{y}, \boldsymbol{C}}{\mathrm{m}\mathrm{i}\mathrm{n}}\left({‖\boldsymbol{S}-\boldsymbol{A}\boldsymbol{L}‖}_{2}^{2}+\mu {‖{\boldsymbol{L}}_{\mathrm{y}}‖}_{p}^{p}+\right.\\ \alpha {‖{\boldsymbol{L}}_{r}-\mathrm{l}\mathrm{n}{\boldsymbol{r}}_{0}‖}_{2}^{2}+\beta {‖{\boldsymbol{L}}_{{v}_{\mathrm{P}}}-\mathrm{l}\mathrm{n}{\boldsymbol{v}}_{\mathrm{P}0}‖}_{2}^{2}+\\ \left.\gamma {‖{\boldsymbol{L}}_{\rho }-\mathrm{l}\mathrm{n}{\boldsymbol{\rho }}_{0}‖}_{2}^{2}+\lambda {‖{\boldsymbol{L}}_{\mathrm{y}}-{\boldsymbol{A}}_{1}\boldsymbol{L}-\boldsymbol{C}‖}_{2}^{2}\right)\end{array} $ (27)

式中λ为对偶项的约束参数。将含有多参数的目标函数分解为多个单参数问题,其中关于$ \boldsymbol{L} $的子目标函数为

$ \begin{array}{l}{J}_{1}\left(\boldsymbol{L}\right)=\underset{\boldsymbol{L}}{\mathrm{m}\mathrm{i}\mathrm{n}}\left({‖\boldsymbol{S}-\boldsymbol{A}\boldsymbol{L}‖}_{2}^{2}+\alpha {‖{\boldsymbol{L}}_{r}-\mathrm{l}\mathrm{n}{\boldsymbol{r}}_{0}‖}_{2}^{2}+\right.\\ \beta {‖{\boldsymbol{L}}_{{v}_{\mathrm{P}}}-\mathrm{l}\mathrm{n}{\boldsymbol{v}}_{\mathrm{P}0}‖}_{2}^{2}+\gamma {‖{\boldsymbol{L}}_{\rho }-\mathrm{l}\mathrm{n}{\boldsymbol{\rho }}_{0}‖}_{2}^{2}+\\ \left.\lambda {‖{\boldsymbol{L}}_{\mathrm{y}}-{\boldsymbol{A}}_{1}\boldsymbol{L}-\boldsymbol{C}‖}_{2}^{2}\right)\end{array} $ (28)

拉格朗日项$ {\boldsymbol{L}}_{\mathrm{y}} $的子目标函数为

$ {J}_{2}\left({\boldsymbol{L}}_{\mathrm{y}}\right)=\underset{{\boldsymbol{L}}_{\mathrm{y}}}{\mathrm{m}\mathrm{i}\mathrm{n}}\left(\mu {‖{\boldsymbol{L}}_{\mathrm{y}}‖}_{p}^{p}+\lambda {‖{\boldsymbol{L}}_{\mathrm{y}}-{\boldsymbol{A}}_{1}\boldsymbol{L}-\boldsymbol{C}‖}_{2}^{2}\right) $ (29)

关于对偶项$ \boldsymbol{C} $的子目标函数为

$ {J}_{3}\left(\boldsymbol{C}\right)=\underset{\boldsymbol{C}}{\mathrm{m}\mathrm{i}\mathrm{n}}\lambda {‖{\boldsymbol{L}}_{\mathrm{y}}-{\boldsymbol{A}}_{1}\boldsymbol{L}-\boldsymbol{C}‖}_{2}^{2} $ (30)

引入交替方向乘子算法及软阈值收缩算法对上述不同目标函数交替求解

$ \begin{array}{l}{\boldsymbol{L}}^{(k+1)}={\left({\boldsymbol{A}}^{\mathrm{T}}\boldsymbol{A}+\boldsymbol{\eta }\boldsymbol{E}+\lambda {\boldsymbol{A}}_{1}^{\mathbf{T}}{\boldsymbol{A}}_{1}\right)}^{-1}\left[{\boldsymbol{A}}^{\mathrm{T}}\boldsymbol{S}+\right.\\ \begin{array}{cc}& \left.\boldsymbol{\eta }{\boldsymbol{L}}_{0}+\lambda {\boldsymbol{A}}_{1}^{\mathbf{T}}{\boldsymbol{L}}_{\mathrm{y}}^{\left(k\right)}-\lambda {\boldsymbol{A}}_{1}^{\mathbf{T}}{\boldsymbol{C}}^{\left(k\right)}\right]\end{array}\end{array} $ (31)
$ \begin{array}{l}{{\boldsymbol{L}}_{\mathrm{y}}}^{(k+1)}=\mathrm{m}\mathrm{a}\mathrm{x}\left[\left|{\boldsymbol{A}}_{1}{\boldsymbol{L}}^{(k+1)}+{\boldsymbol{C}}^{\left(k\right)}\right|-{(\lambda /\mu )}^{p-2}\times \right.\\ \left.{\left|{\boldsymbol{A}}_{1}{\boldsymbol{L}}^{(k+1)}+{\boldsymbol{C}}^{\left(k\right)}\right|}^{p-1}, 0\right]\times \mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left[{\boldsymbol{A}}_{1}{\boldsymbol{L}}^{(k+1)}+{\boldsymbol{C}}^{\left(k\right)}\right]\end{array} $ (32)
$ {\boldsymbol{C}}^{(k+1)}={\boldsymbol{C}}^{\left(k\right)}+{\boldsymbol{A}}_{1}{\boldsymbol{L}}^{(k+1)}-{{\boldsymbol{L}}_{\mathrm{y}}}^{(k+1)} $ (33)

式中:$ {\boldsymbol{L}}_{0}={\left[\begin{array}{ccc}\mathrm{l}\mathrm{n}{\boldsymbol{r}}_{0}^{\mathrm{T}}& \mathrm{l}\mathrm{n}{\boldsymbol{v}}_{\mathrm{P}0}^{\mathrm{T}}& \mathrm{l}\mathrm{n}{\boldsymbol{\rho }}_{0}^{\mathrm{T}}\end{array}\right]}^{\mathrm{T}} $k为迭代次数;E为单位矩阵;$ \left|•\right| $表示对矩阵每个元素取绝对值;max(·,0)表示比较矩阵各个元素与0,取最大值,即对矩阵中值为负的元素置零;$ {\left|•\right|}^{p-1} $表示对矩阵每个元素取绝对值后计算p$ - $1次方;sign(·)表示对矩阵每个元素取正、负号;而

$ \boldsymbol{\eta }=\left(\begin{array}{ccc}\alpha & 0& 0\\ 0& \beta & 0\\ 0& 0& \gamma \end{array}\right) $ (34)
2 模型测试与实际应用 2.1 模型测试

为验证本文提出反演方法的效果,根据A井测井数据设计了一维地震模型。该模型是利用A井纵波速度、横波速度和密度测井数据通过Zoeppritz方程计算得到纵波反射系数,并与主频为30 Hz的Ricker子波褶积得到不同角度的合成地震道集(图 2a)。为了验证本文方法的抗噪性,对合成地震数据加上其最大振幅30%的高斯噪声(图 2b)。

图 2 合成地震道集数据 (a)加噪前;(b)加噪后

将Lp拟范数稀疏约束的纵横波速比直接反演方法应用于该模型。图 3为无噪条件下的反演结果,与理论模型很接近,验证了本文方法的可行性。图 4为含噪条件下的反演结果,可以发现:纵波速度和密度由于噪声的影响,反演精度降低,但仍能反映基本趋势;纵横波速比反演结果受噪声干扰较小,仍具有较高的精度,说明本文方法具有较强的抗噪性。

图 3 无噪合成数据本文方法的反演结果(红色)与理论模型(黑色)及初始模型(绿色)的对比

图 4 含噪合成数据本文方法反演结果(红色)与理论模型(黑色)及初始模型(绿色)的对比

为了进一步验证直接反演方法的优势,将利用Aki-Richards近似方程间接得到的纵横波速比与直接得到的反演结果进行对比(图 5)。由图 5可见,无论是无噪条件下还是含噪条件下,直接反演结果均更接近理论模型,尤其是黑色虚线圈处更加明显。

图 5 直接(红色)和间接(蓝色)反演结果与理论模型(黑色)及初始模型(绿色)的对比 (a)无噪;(b)含噪

为了定量分析提出方法与间接反演方法的反演效果,使用信噪比(Signal-to-noise Ratio,SNR)定量评估不同反演方法的反演效果。具体定义为

$ \mathrm{S}\mathrm{N}\mathrm{R}=10\times \mathrm{l}\mathrm{g}\frac{{‖\boldsymbol{r}-\stackrel{-}{\boldsymbol{r}}‖}_{2}^{2}}{{‖{\boldsymbol{r}}^{\mathrm{i}\mathrm{n}\mathrm{v}}-\stackrel{-}{\boldsymbol{r}}‖}_{2}^{2}} $ (35)

式中:r为理论纵横波速度比模型;$ \stackrel{-}{\boldsymbol{r}} $表示理论模型的平均值;$ {\boldsymbol{r}}^{\mathrm{i}\mathrm{n}\mathrm{v}} $表示反演结果。

表 2为两种方法反演结果的SNR,可见,无论在无噪或含噪的条件下,本文方法都具有更高的信噪比,从而验证了方法的有效性。

表 2 两种方法反演结果的SNR
2.2 实际数据应用

将本文的直接反演方法应用于M工区实际地震数据。该区目的层主要为砂岩含气储层。图 6是入射角为10°、17°和24°的部分角度叠加剖面,共500道,每道含有241个采样点,采样间隔为1 ms。该地震剖面过A、B两口井,均钻遇含气储层,因此具有较高的研究价值。图 7为本文方法反演的纵横波速比、纵波速度和密度剖面,从图中可以看出,反演剖面与A、B两井测井曲线趋势吻合较好,验证了本文方法实用性。

图 6 不同角度的叠加地震剖面 (a)10°;(b)17°;(c)24°

图 7 本文方法反演的vP/vS(a)、vP(b)和ρ(c)剖面

对比本文方法得到的vP/vS反演结果与基于Aki-Richards近似方程得到的间接反演结果(图 8)可进一步说明本文方法的优越性。气层的纵横波速比反映为异常低值,与测井解释结论对比可知,两种方法得到的反演剖面与测井解释结论基本吻合,但是由于间接反演方法存在累积误差,对部分含气储层不能有效识别(黑色椭圆处)。并且由于直接反演方法得到的反演结果不存在累积误差,其对储层边界刻画更为清晰(黑色矩形框处)。

图 8 本文方法(a)与间接方法(b)反演的vP/vS剖面对比 井柱上绿色表示含气储层,蓝色代表非气藏层。

为了进一步验证本文方法应用于实际数据的可行性,对该工区进行三维反演。数据由500条Xline和200条Inline组成。图 9为本文方法的三维直接反演结果,可见,横向连续性较高、气藏分布较清晰。图 10是提取三维反演结果的沿层切片,可以看出,低值区(红色)显示为气藏异常,与A井和B井吻合较好,验证了本文方法的有效性。

图 9 本文方法的三维反演结果

图 10 本文方法三维反演结果的沿层切片
3 结论

(1) 本文在广义弹性阻抗基础上提出了新的广义弹性阻抗方程,该方程与纵横波速比、纵波速度和密度有关,计算的纵波反射系数更接近Zoeppritz方程。为了直接反演纵横波速比,由新的广义弹性阻抗方程导出了一个与纵横波速比、纵波速度、密度相关的纵波反射系数近似方程。

(2) 为了提高反演精度,在稀疏反演理论的基础上引入Lp拟范数对反演参数进行稀疏约束,应用交替方向乘子算法求解。模型数据测试及实际数据应用结果表明本文方法反演精度较高。

(3) 本文的直接反演方法避免了间接反演方法产生的累积误差,具有一定的抗噪能力,提高了反演的稳定性,同时提高了对含气储层的识别能力。

参考文献
[1]
李建华, 刘百红, 张延庆, 等. 叠前AVO反演在储层含油气性预测中的应用[J]. 石油地球物理勘探, 2016, 51(6): 1180-1186.
LI Jianhua, LIU Baihong, ZHANG Yanqing, et al. Oil-bearing reservoir prediction with prestack AVO inversion[J]. Oil Geophysical Prospecting, 2016, 51(6): 1180-1186.
[2]
郭强, 雒聪, 刘红达, 等. 自适应优化参数模拟退火的叠前地震联合反演方法[J]. 石油地球物理勘探, 2023, 58(3): 670-679.
GUO Qiang, LUO Cong, LIU Hongda, et al. Prestack seismic hybrid inversion based on simulated annealing algorithm with adaptive optimization parameters[J]. Oil Geophysical Prospecting, 2023, 58(3): 670-679.
[3]
徐斌, 陈学华, 张杰, 等. 依赖频率的纵、横波衰减参数叠前反演方法[J]. 石油地球物理勘探, 2022, 57(6): 1427-1435.
XU Bin, CHEN Xuehua, ZHANG Jie, et al. Pre-stack inversion method for frequency-dependent P-wave and S-wave attenuation parameters[J]. Oil Geophysical Prospecting, 2022, 57(6): 1427-1435.
[4]
李坤, 印兴耀. 混合概率模型驱动的叠前地震反演方法[J]. 石油地球物理勘探, 2020, 55(4): 839-853.
LI Kun, YIN Xingyao. Prestack seismic inversion driven by mixture probabilistic models[J]. Oil Geophysical Prospecting, 2020, 55(4): 839-853.
[5]
陈建江, 印兴耀, 张广智. 层状介质AVO叠前反演[J]. 石油地球物理勘探, 2006, 41(6): 656-662.
CHEN Jianjiang, YIN Xingyao, ZHANG Guangzhi. Prestack AVO inversion of layered medium[J]. Oil Geophysical Prospecting, 2006, 41(6): 656-662.
[6]
CHEN F B, ZONG Z Y, YANG Y M. Amplitude-variation-with-offset inversion using P- to S-wave velocity ratio and P-wave velocity[J]. Geophysics, 2022, 87(4): N63-N74. DOI:10.1190/geo2021-0623.1
[7]
WANG X Q, SCHUBNEL A, FORTIN J, et al. High vP/vS ratio: Saturated cracks or anisotropy effects?[J]. Geophysical Research Letters, 2012, 39(11): L11307.
[8]
YIN X Y, ZHANG S X. Bayesian inversion for effective pore-fluid bulk modulus based on fluid-matrix decoupled amplitude variation with offset approximation[J]. Geophysics, 2014, 79(5): R221-R232. DOI:10.1190/geo2013-0372.1
[9]
张金强, 刘振峰, 刘喜武, 等. 致密砂岩储层自适应岩石物理建模方法[J]. 石油地球物理勘探, 2021, 56(4): 792-800, 808.
ZHANG Jinqiang, LIU Zhenfeng, LIU Xiwu, et al. Self-adaptive rock physics modeling method for tight sandstone reservoirs[J]. Oil Geophysical Prospecting, 2021, 56(4): 792-800, 808.
[10]
程冰洁, 徐天吉, 梁群, 等. 多波速度比参数含气性识别研究与应用[J]. 石油地球物理勘探, 2014, 49(2): 307-315.
CHENG Bingjie, XU Tianji, LIANG Qun, et al. Gas identification using multi-wave velocity ratio parameters[J]. Oil Geophysical Prospecting, 2014, 49(2): 307-315.
[11]
张凌远, 张宏兵, 尚作萍, 等. 基于Zoeppritz方程的叠前和叠后混合多参数非线性地震反演[J]. 石油地球物理勘探, 2021, 56(1): 164-171.
ZHANG Lingyuan, ZHANG Hongbing, SHANG Zuoping, et al. Nonlinear multi-parameter hybrid inversion of pre-stack and post-stack seismic data based on Zoeppritz equation[J]. Oil Geophysical Prospecting, 2021, 56(1): 164-171.
[12]
CONNOLLY P. Elastic impedance[J]. The Leading Edge, 1999, 18(4): 438-452. DOI:10.1190/1.1438307
[13]
WHITCOMBE D N. Elastic impedance normalization[J]. Geophysics, 2002, 67(1): 60-62. DOI:10.1190/1.1451331
[14]
WHITCOMBE D N, CONNOLLY P A, REAGAN R L, et al. Extended elastic impedance for fluid and lithology prediction[J]. Geophysics, 2002, 67(1): 63-67. DOI:10.1190/1.1451337
[15]
马劲风. 地震勘探中广义弹性阻抗的正反演[J]. 地球物理学报, 2003, 46(1): 118-124.
MA Jinfeng. Forward modeling and inversion method of generalized elastic impedance in seismic exploration[J]. Chinese Journal of Geophysics, 2003, 46(1): 118-124.
[16]
印兴耀, 刘晓晶, 吴国忱, 等. 模型约束基追踪反演方法[J]. 石油物探, 2016, 55(1): 115-122.
YIN Xingyao, LIU Xiaojing, WU Guochen, et al. Basis pursuit inversion method under model constraint[J]. Geophysical Prospecting for Petroleum, 2016, 55(1): 115-122.
[17]
刘晓晶, 印兴耀, 吴国忱, 等. 基于基追踪弹性阻抗反演的深部储层流体识别方法[J]. 地球物理学报, 2016, 59(1): 277-286.
LIU Xiaojing, YIN Xingyao, WU Guochen, et al. Identification of deep reservoir fluids based on basis pursuit inversion for elastic impedance[J]. Chinese Journal of Geophysics, 2016, 59(1): 277-286.
[18]
SHE B, WANG Y J, ZHANG J S, et al. AVO inversion with high-order total variation regularization[J]. Journal of Applied Geophysics, 2019, 161: 167-181. DOI:10.1016/j.jappgeo.2018.12.014
[19]
文晓涛, 杨吉鑫, 李雷豪, 等. 低频稀疏双约束宽频带地震阻抗反演[J]. 天然气工业, 2019, 39(5): 45-52.
WEN Xiaotao, YANG Jixin, LI Leihao, et al. Low-frequency sparse double-constrained broadband seismic impedance inversion[J]. Natural Gas Industry, 2019, 39(5): 45-52.
[20]
WOODWORTH J, CHARTRAND R. Compressed sensing recovery via nonconvex shrinkage penalties[J]. Inverse Problems, 2016, 32(7): 075004. DOI:10.1088/0266-5611/32/7/075004
[21]
CHEN Y P, PENG Z M, GHOLAMI A, et al. Seismic signal sparse time-frequency representation by Lp-quasi norm constraint[J]. Digital Signal Processing, 2019, 87(4): 43-59.
[22]
张雨强, 文晓涛, 吴昊等. 基于Lp拟范数稀疏约束和交替方向乘子算法的波阻抗反演[J]. 石油物探, 2022, 61(5): 856-864.
ZHANG Yuqiang, WEN Xiaotao, WU Hao, et al. Seismic acoustic impedance inversion using Lp quasi-norm sparse constraint and alternating direction multiplier algorithm[J]. Geophysical Prospecting for Petroleum, 2022, 61(5): 856-864.
[23]
GOODWAY B, CHEN T W, DOWNTON J. Improved AVA fluid detection and lithology discrimination using Lamé petrophysical parameters "λρ", "μρ", & "λ/μ fluid stack" from P and S inversions[C]. SEG Technical Program Expanded Abstracts, 1997, 16: 183-186.
[24]
AKI K, RICHARDS P G. Quantitative Seismology: Theory and Methods[M]. San Francisco: W H Freeman and Company, 1980.
[25]
BOYD S, PARIKH N, CHU E, et al. Distributed optimization and statistical learning via the alternating direction method of multipliers[J]. Foundations and Trends in Machine Learning, 2011, 3(1): 1-122.