一种高精度的纵横波阻抗反演方法

摘要

本发明公开了一种高精度的纵横波阻抗反演方法，包括：首先，在Fatti近似公式的基础上进一步优化了纵、横波阻抗反射系数近似公式；然后，根据反射系数奇偶分解理论和拉格朗日乘子法，构建出纵、横波阻抗低频模型同步约束的L1‑2最小化反演目标函数；最后，采用凸函数差异算法DCA和交替方法乘子法ADMM联合对其求解，得到纵、横波阻抗。

权利要求

1.一种高精度的纵横波阻抗反演方法，其特征在于，包括：通过对纵、横波阻抗和密度的地震反射系数近似关系进行优化推导处理，得到优化后的纵、横波阻抗反射系数近似关系，其中，所述优化后的纵、横波阻抗反射系数近似关系为：利用地震褶积模型和所述优化后的纵、横波阻抗反射系数近似关系，得到反射系数分解后的地震褶积模型；

利用所述反射系数分解后的地震褶积模型及拉格朗日乘子法，构建纵、横波阻抗低频模型同步约束的反演目标函数；

利用L1-2最小化算法对所述纵、横波阻抗低频模型同步约束的反演目标函数进行反演，得到奇偶分量系数，其包括：通过对所述反演目标函数进行化简处理，得到L1-2最小化反演目标函数；

利用增广拉格朗日乘子法对所述L1-2最小化反演目标函数进行约束，得到最终目标函数利用交替方法乘子法ADMM对所述最终目标函数进行求解，得到奇偶分量系数；

利用所述奇偶分量系数，得到纵、横波阻抗；

其中，θ为入射角，b为实际工区测井资料拟合的参数，分别为t时刻对应的纵波阻

抗和横波阻抗，α为权重参数。

2.根据权利要求1所述的方法，其特征在于，所述利用地震褶积模型和所述优化后的纵、横波阻抗反射系数近似关系，得到反射系数分解后的地震褶积模型包括：通过对所述优化后的纵、横波阻抗反射系数近似关系中的反射系数进行奇偶分解，得到反射系数分解后的纵横、波阻抗反射系数近似关系；

根据地震褶积模型和所述反射系数分解后的纵、横波阻抗反射系数近似关系，得到反射系数分解后的地震褶积模型。

3.根据权利要求2所述的方法，其特征在于，所述反射系数分解后的纵、横波阻抗反射系数近似关系包括：r＝Dm；

其中，所述r是指纵横波阻抗反射系数；所述D是指反射系数分解算子；所述m是指奇偶分量系数。

4.根据权利要求3所述的方法，其特征在于，所述根据地震褶积模型和所述反射系数分解后的纵、横波阻抗反射系数近似关系，得到反射系数分解后的地震褶积模型包括：通过将所述反射系数分解后的纵、横波阻抗反射系数近似关系r＝Dm代入所述地震褶积模型s＝WCr中，得到反射系数分解后的地震褶积模型：s＝WCDm；

其中，所述s是指多角度叠前地震记录；所述W是指w(θ)组成的多角度子波矩阵；所述C是指由Cp(θ)、Cs(θ)组成的系数矩阵。

5.根据权利要求4所述的方法，其特征在于，所述纵、横波阻抗同步约束的反演目标函数包括：其中，所述β1是指纵波阻抗低频约束的权重系数；所述β2是指横波阻抗低频约束的权重系数；所述 所述 所述 以及所述B均是指积分算子；所述Ip_low是指纵波阻抗低频模型；所述Is_low是指横波阻抗低频模型；所述LLP为低通滤波矩阵。

6.根据权利要求5所述的方法，其特征在于，所述利用所述反射系数分解后的地震褶积模型及拉格朗日乘子法，构建纵、横波阻抗低频模型同步约束的反演目标函数包括：根据测井资料，分别建立纵波低频模型和横波低频模型；

根据所述拉格朗日乘子法、所述纵波低频模型以及所述横波低频模型对基追踪标准公式进行同步约束，得到所述反演目标函数。

7.根据权利要求1所述的方法，其特征在于，所述利用所述奇偶分量系数，得到纵、横波阻抗包括：根据所述奇偶分量系数，得到纵、横波反射系数；

通过对所述纵、横波反射系数进行积分处理，得到纵、横波阻抗。

说明书

一种高精度的纵横波阻抗反演方法

技术领域

[0001] 本发明涉及地震叠前反演技术领域，特别涉及一种高精度的纵横波阻抗反演方法。

背景技术

[0002] 地震叠前反演是储层流体识别和评价的有效途径。叠前AVO(Amplitude Variation with Offset，带偏移量的振幅变化)反演方法，是获取流体识别因子的主要技术之一，该反演方法能够综合利用测井、岩石物理等资料，因此，在实际应用中备受关注。

[0003] 叠前AVO反演方法的核心包括反射系数近似公式和反演算法两部分。许多学者通过对Zoeppritz方程进行简化，得到不同形式的反射系数近似公式。然而，所有的反演方法均是基于L1范数的凸优化算法，虽然得到了较为广泛的应用，但考虑到地震反演问题本身固有的“病态解”和“多解性”等问题，并指出凸优化算法会导致反演结果陷入次优稀疏解，为解决这一问题，非凸优化算法引起了众多学者的关注。其中，利用L1-2最小化算法对地震衰减进行补偿，并取得较好的实际应用效果。但是，在常用的AVO反演方法存在稀疏性差、垂向分辨率低以及反演精度低的问题依然没有得到良好的改善。

发明内容

[0004] 根据本发明实施例提供的方案解决的技术问题是常用的AVO反演方法存在稀疏性差、垂向分辨率低以及反演精度不高。

[0005] 根据本发明实施例提供的一种高精度的纵横波阻抗反演方法，包括：

[0006] 通过对纵、横波阻抗和密度的地震反射系数近似关系进行优化推导处理，得到优化后的纵、横波阻抗反射系数近似关系；

[0007] 通过对所述优化后的纵、横波阻抗反射系数近似关系中的反射系数进行奇偶分解，得到反射系数分解后的纵横、波阻抗反射系数近似关系；

[0008] 根据地震褶积模型和所述反射系数分解后的纵、横波阻抗反射系数近似关系，得到反射系数分解后的地震褶积模型；

[0009] 利用所述反射系数分解后的地震褶积模型及拉格朗日乘子法，构建纵、横波阻抗低频模型同步约束的反演目标函数；

[0010] 利用L1-2最小化算法对所述纵、横波阻抗低频模型同步约束的反演目标函数进行反演，得到奇偶分量系数；

[0011] 利用所述奇偶分量系数，得到纵、横波阻抗。

[0012] 优选地，所述优化后的纵、横波阻抗反射系数近似关系包括：

[0013]

[0014] 其中，θ为入射角，b为实际工区测井资料拟合的参数，
分别为t时刻对应的纵波阻抗和横波阻抗。

[0015] 优选地，所述反射系数分解后的纵、横波阻抗反射系数近似关系包括：

[0016] r＝Dm；

[0017] 其中，所述r是指纵横波阻抗反射系数；所述D是指反射系数分解算子；所述m是指奇偶分量系数。

[0018] 优选地，所述根据地震褶积模型和所述反射系数分解后的纵、横波阻抗反射系数近似关系，得到反射系数分解后的地震褶积模型包括：

[0019] 通过将所述反射系数分解后的纵、横波阻抗反射系数近似关系r＝Dm代入所述地震褶积模型s＝WCr中，得到反射系数分解后的地震褶积模型：s＝WCDm；

[0020] 其中，所述s是指多角度叠前地震记录；所述W是指w(θ)组成的多角度子波矩阵；所述C是指由Cp(θ)、Cs(θ)组成的系数矩阵。

[0021] 优选地，所述纵、横波阻抗同步约束的反演目标函数包括：

[0022]

[0023] 其中，所述β1是指纵波阻抗低频约束的权重系数；所述是指横波阻抗低频约束的权重系数；所述 所述 所述 以及所述B均是指积分算子；所述Ip_low是指纵波阻抗低频模型；所述Is_low是指横波阻抗低频模型；所述LLP为低通滤波矩阵。

[0024] 优选地，所述利用所述反射系数分解后的地震褶积模型及拉格朗日乘子法，构建纵、横波阻抗低频模型同步约束的反演目标函数包括：

[0025] 根据测井资料，分别建立纵波低频模型和横波低频模型；

[0026] 根据所述拉格朗日乘子法、所述纵波低频模型以及所述横波低频模型对基追踪标准公式进行同步约束，得到所述反演目标函数。

[0027] 优选地，所述利用L1-2最小化算法对所述纵横波阻抗低频模型同步约束的反演目标函数进行反演，得到奇偶分量系数包括：

[0028] 通过对所述反演目标函数进行化简处理，得到L1-2最小化反演目标函数；

[0029] 利用增广拉格朗日乘子法对所述L1-2最小化反演目标函数进行约束，得到最终目标函数；

[0030] 利用交替方法乘子法ADMM对所述最终目标函数进行求解，得到奇偶分量系数。

[0031] 优选地，所述利用所述奇偶分量系数，得到纵、横波阻抗包括：

[0032] 根据所述奇偶分量系数，得到纵、横波反射系数；

[0033] 通过对所述纵、横波反射系数进行积分处理，得到纵、横波阻抗。

[0034] 根据本发明实施例提供的方案，在稀疏性、反演精度、抗噪性等方面明显优于常用的基追踪算法(BasisPursuit，BP)，同时实际资料的应用证实了新方法的可行性，其反演结果为储层流体识别提供了更可靠的依据。

附图说明

[0035] 图1是本发明实施例提供的一种高精度的纵横波阻抗反演方法流程图；

[0036] 图2是本发明实施例提供的不同方程下反射系数对比图；

[0037] 图3是本发明实施例提供的反射系数奇偶分解示意图；

[0038] 图4是本发明实施例提供的LP范数单位球示意图。

具体实施方式

[0039] 以下结合附图对本发明的优选实施例进行详细说明，应当理解，以下所说明的优选实施例仅用于说明和解释本发明，并不用于限定本发明。

[0040] 图1是本发明实施例提供的一种高精度的纵横波阻抗反演方法流程图，如图1所示，包括：

[0041] 步骤S101：通过对纵、横波阻抗和密度的地震反射系数近似关系进行优化推导处理，得到优化后的纵、横波阻抗反射系数近似关系；

[0042] 步骤S102：通过对所述优化后的纵、横波阻抗反射系数近似关系中的反射系数进行奇偶分解，得到反射系数分解后的纵横、波阻抗反射系数近似关系；

[0043] 步骤S103：根据地震褶积模型和所述反射系数分解后的纵、横波阻抗反射系数近似关系，得到反射系数分解后的地震褶积模型；

[0044] 步骤S104：利用所述反射系数分解后的地震褶积模型及拉格朗日乘子法，构建纵、横波阻抗低频模型同步约束的反演目标函数；

[0045] 步骤S105：利用L1-2最小化算法对所述纵、横波阻抗低频模型同步约束的反演目标函数进行反演，得到奇偶分量系数；

[0046] 步骤S106：利用所述奇偶分量系数，得到纵、横波阻抗。

[0047] 其中，所述优化后的纵、横波阻抗反射系数近似关系包括：

[0048]

[0049] 其中，θ为入射角，b为实际工区测井资料拟合的参数，
分别为t时刻对应的纵波阻抗和横波阻抗。

[0050] 其中，所述反射系数分解后的纵、横波阻抗反射系数近似关系包括：r＝Dm；其中，所述r是指纵横波阻抗反射系数；所述D是指反射系数分解算子；所述m是指奇偶分量系数。

[0051] 其中，所述根据地震褶积模型和所述反射系数分解后的纵、横波阻抗反射系数近似关系，得到反射系数分解后的地震褶积模型包括：通过将所述反射系数分解后的纵、横波阻抗反射系数近似关系r＝Dm代入所述地震褶积模型s＝WCr中，得到反射系数分解后的地震褶积模型：s＝WCDm；其中，所述s是指多角度叠前地震记录；所述W是指w(θ)组成的多角度子波矩阵；所述C是指由Cp(θ)、Cs(θ)组成的系数矩阵。

[0052] 其中，所述纵、横波阻抗同步约束的反演目标函数包括：其中，所述β1是指纵波阻抗
低频约束的权重系数；所述是指横波阻抗低频约束的权重系数；所述 所述所述 以及所述B是指积分算子；所述Ip_low是指纵波阻抗低频模型；所
述Is_low是指横波阻抗低频模型；所述LLP为低通滤波矩阵。

[0053] 其中，所述利用所述反射系数分解后的地震褶积模型及拉格朗日乘子法，构建纵、横波阻抗低频模型同步约束的反演目标函数包括：根据测井资料，分别建立纵波低频模型和横波低频模型；根据所述拉格朗日乘子法、所述纵波低频模型以及所述横波低频模型对基追踪标准公式进行同步约束，得到所述反演目标函数。

[0054] 其中，所述利用L1-2最小化算法对所述纵横波阻抗低频模型同步约束的反演目标函数进行反演，得到奇偶分量系数包括：通过对所述反演目标函数进行化简处理，得到L1-2最小化反演目标函数；利用增广拉格朗日乘子法对所述L1-2最小化反演目标函数进行约束，得到最终目标函数；利用交替方法乘子法ADMM对所述最终目标函数进行求解，得到奇偶分量系数。

[0055] 其中，所述利用所述奇偶分量系数，得到纵、横波阻抗包括：根据所述奇偶分量系数，得到纵、横波反射系数；通过对所述纵、横波反射系数进行积分处理，得到纵、横波阻抗。

[0056] 纵、横波阻抗反射系数近似方程推导

[0057] 纵、横波阻抗和密度的地震反射系数近似关系，称为Fatti三项式，其形式如公式(1)：

[0058]

[0059] 式中: 同时，Fatti等人发现在 入射角度不大的情况下，可将公式(1)简化为：

[0060]

[0061] 公式(2)一般称作Fatti两项式，该公式适用于缺乏大角度地震数据的情况。为了能够将叠前多角度数据均用于反演，同时考虑到密度项是一个难以准确反演的参数，根据Gardner经验公式：

[0062]

[0063] 式中，a,b为根据实际工区测井资料拟合的参数。由公式(3)可得密度反射系数与纵波速度反射系数的关系如下：

[0064]

[0065] 考虑均匀各向同性完全弹性无孔介质中，纵波阻抗与纵波速度、密度的关系式Ip＝ρVp，并对其求导可得：

[0066]

[0067] 结合公式(4)和(5)，可推导出密度和纵波阻抗的关系：

[0068]

[0069] 将公式(6)带入公式(1)可推导出纵、横波阻抗反射系数近似方程：

[0070]

[0071] 表1 GoodWay含气砂岩与页岩模型(Goodway，1997)

[0072]

[0073] 利用给出的含气砂岩模型，如表1所示，对公式(7)、精确Zoeppritz方程和Fatti近似式计算的反射系数进行了对比分析，结果如图2所示，由此可见，本发明方法得到的反射系数近似方程精度高于忽略密度项的Fatti两项式。

[0074] 正演模型

[0075] 上述内容推导了叠前多角度地震记录与纵、横波阻抗反射系数之间的关系。考虑多个界面，N个角度的情况下，反射系数近似方程写成矩阵形式，如公式(8)：

[0076]

[0077] 式中，rp、rs分别为纵波、横波阻抗反射系数；Cp(θ)、Cs(θ)分别由 B(θ)组成的对角阵系数矩阵。根据地震褶积模型，考虑多个界面、N个角度情况下，推导出地震记录与反射系数之间的关系，如公式(9)：

[0078]

[0079] 式中，w(θ)为角度子波矩阵。将公式(8)、(9)结合，并简写成矩阵形式，如公式(10)：

[0080] s＝WCr  (10)

[0081] 式中， W为由w(θ)组成的多角度子波矩阵，C是由Cp(θ)、Cs(θ)组成的系数矩阵，具体形式见公式(8)、(9)。

[0082] L1-2最小化反演方法

[0083] 采用反射系数奇偶分解理论，有效提高了反演结果的垂向分辨率。因此，文中采用该方法对纵、横波阻抗反射系数进行奇偶分解，并通过拉格朗日乘子法构建出纵、横波阻抗同步约束的反演目标函数，然后采用L1-2最小化算法进行了纵、横波阻抗的精确反演。

[0084] (1)反射系数奇偶分解

[0085] 根据奇偶分解理论：一个反射系数对可以唯一分解为一个反射系数偶分量和一个反射系数奇分量，如图3所示。将反射系数r进行奇偶分解，如公式(11)：

[0086] r＝are+bro  (11)

[0087] 式中，re为反射系数偶分量；ro为反射系数奇分量；a、b为奇偶分量分解系数。

[0088] 考虑同步反演纵、横波阻抗反射系数，结合公式(11)，将纵、横波反射系数进行奇偶分解，如公式(12)：

[0089]

[0090] 将公式(12)代入(10)中，简写为公式(13)：

[0091] s＝WCDm  (13)

[0092] 式中，D表示反射系数分解算子，m表示奇偶分量对应的系数，即待反演量。D、m具体形式见公式(12)。

[0093] (2)纵、横波低频模型约束的反演目标函数

[0094] 采用基追踪对公式(13)进行反演，其求解形式定义如下：

[0095]

[0096] 式中，λ为正则化调节参数，用来控制待反演参数的稀疏程度。若采用公式(14)进行直接反演，会引起缺乏实际地质背景而导致反演结果横向连续性差。为此根据拉格朗日乘子法，利用测井资料建立纵、横波低频模型对公式(14)进行同步约束，其形式如公式(15)：

[0097]

[0098] 式中，β1、β2分别为纵波阻抗、横波阻抗低频约束的权重系数，B、Ip_low、Is_low分别是积分算子和根据测井资料建立的纵波阻抗、横波阻抗低频模型；LLP为低通滤波矩阵。

[0099] (3)L1-2最小化反演方法

[0100] 近年来，提出的L1-2最小化算法被应用于稀疏反演，且其性能明显优于传统的L1范数算法，图4给出了L1和L1-2范数单位球示意图，由图4可知，L1-2最小化算法更接近最优稀疏解。

[0101] 令 则将公式(15)进一步化简为L1-2最小化算法求解形式，如公式(16)所示：

[0102]

[0103] 式中：α为权重参数。采用凸函数差异算法(Difference of Convex Algorithm，DCA)将公式(16)分解为F(m)＝G(m)-H(m)，其中：

[0104]

[0105] 根据DCA迭代公式(Yin et al.,2015)，将L1-2最小化反演目标函数化简为：

[0106]

[0107] 式中，yk为H(m)在mk处的梯度，其定义为：

[0108]

[0109] 引入增广拉格朗日乘数法对公式(18)进行约束，可得：

[0110]

[0111] 式中，z为辅助中间变量，w为拉格朗日乘子，ρ为惩罚参数。公式(20)由交替方法乘子法(Alternating Direction Method of Multipliers，ADMM)求解，其迭代递推公式如：

[0112]

[0113] 最后，利用反演结果进行道积分，可得出对应的纵、横波阻抗，如公式(22)：

[0114]

[0115] 综上所述，首先，在Fatti近似公式的基础上进一步优化了纵、横波阻抗反射系数近似公式；然后，根据反射系数奇偶分解理论和拉格朗日乘子法，构建出纵、横波阻抗低频模型同步约束的L1-2最小化反演目标函数；最后，采用DCA(DifferenceofConvexAlgorithm，凸函数差异算法)和ADMM(Alternating Direction Method of Multipliers，交替方法乘子法)联合对其求解。

[0116] 根据本发明实施例提供的方案，在稀疏性、反演精度、抗噪性等方面明显优于常用的基追踪算法(BasisPursuit，BP)，同时实际资料的应用证实了新方法的可行性，其反演结果为储层流体识别提供了更可靠的依据。

[0117] 尽管上文对本发明进行了详细说明，但是本发明不限于此，本技术领域技术人员可以根据本发明的原理进行各种修改。因此，凡按照本发明原理所作的修改，都应当理解为落入本发明的保护范围。