CN111368247B

CN111368247B - 基于快速正交字典的稀疏表征正则化叠前avo反演方法

Info

Publication number: CN111368247B
Application number: CN202010168773.6A
Authority: CN
Inventors: 王峣钧; 刘宇; 厍斌; 胡光岷
Original assignee: University of Electronic Science and Technology of China
Current assignee: University of Electronic Science and Technology of China
Priority date: 2020-03-12
Filing date: 2020-03-12
Publication date: 2021-11-30
Anticipated expiration: 2040-03-12
Also published as: CN111368247A

Abstract

本发明公开一种基于快速正交字典的稀疏表征正则化叠前AVO反演方法，应用于地球物理勘探解释技术领域，针对现有基于KSVD字典反演的运算速率慢且受控参数多的问题，本发明根据已知的角道集地震数据、子波序列、入射角信息和初始模型参数，构建目标函数；然后加入正则化约束；最后基于正交字典迭代求解正则化约束目标函数，得到最终的模型参数；本发明的方法在保证与KSVD字典方法具有相当反演效果的同时，能有效提升稀疏反演速率。

Description

基于快速正交字典的稀疏表征正则化叠前AVO反演方法

技术领域

本发明属于地球物理勘探解释技术领域，特别涉及一种利用正交字典学习结合稀疏表征正则化约束的叠前AVO反演技术。

背景技术

AVO(Amplitude variation with offset,振幅随偏移距的变化)反演试图将叠前振幅转化为反映地下介质岩性和含油气信息的有效弹性参数，是当前最为重要的储层预测手段之一。由于先验信息不足，数据受噪声干扰严重、反演参数多等原因，使得AVO反演具有多解性和不稳定性的特点。减少反演结果非唯一性和不稳定性是反演算法的核心问题(Samet al，2015)，正则化技术便是最常用的提高反演过程稳定性的方法之一。正则化方法首先由Tikhonov提出(Tikhonov and Arsenin,1977)。后来，许多科学家在Tikhonov框架的基础上发展了各种正则化方法。通常这类正则化方法假设地层模型服从特定的分布，比如高斯分布或柯西分布，这样我们可以通过解析化的手段对模型参数进行表征，将解析化的先验信息用于约束反演求解过程，即可降低反演多解性。这类方法已经被证明是行之有效的技术手段，但是随着我们对储层分析的深入，这类对储层假设统一数学模型的方法已经很难适应地质特征不同区域的反演问题(Bin She et al,2018,2019)。

近年来有学者在地震反演中引入了数据驱动先验信息学习及其约束反演方法，取得了显著的效果(Bin She et al,2018,2019)。该方法采用KSVD字典学习方法从测井数据中学习地层特征，利用字典的稀疏表征作为先验信息实现对地震反演的约束，该方法使得反演先验信息获取从数学分布假设变为从数据中获取，避免了工区单一数学模型假设问题，提升了先验信息对地层特征刻画的全面性和有效性。这类数据驱动AVO反演方法的核心是字典的学习和基于字典的稀疏表征过程。目前基于KSVD字典反演方法假设其字典原子冗余，即原子之间存在信息重叠，在该算法中使用冗余度较高的字典原子能够提高表征系数的稀疏性，但也造成了字典学习过程迭代的增加和小块冗余重叠，因此反演过程存在运算效率较慢且受控参数多等问题，极大的限制了算法的实际应用。

发明内容

为解决上述技术问题，本发明提出一种基于快速正交字典的稀疏表征正则化叠前AVO反演方法，在能够保证反演效果与原有基于KSVD字典学习叠前AVO反演方法相当的条件下，本发明的方法参数调节更加便捷，运算效率更高，更适合实际地震数据反演。

本发明采用的技术方案为：一种基于快速正交字典的稀疏表征正则化叠前AVO反演方法，包括：

A、根据已知的角道集地震数据、子波序列、入射角信息和反演模型参数，构建目标函数；

B、加入正则化约束；

C、基于正交字典迭代求解正则化约束目标函数，得到最终的模型参数。

步骤A中所述的模型参数包括：纵波速度、横波速度以及密度；所述反演模型参数为：纵波速度的初始值、横波速度的初始值以及密度的初始值。

步骤B所述加入正则化约束的目标函数为：

其中，

用于计算反演过程的均方误差，函数R(x)是关于模型参数x的正则化项，argmin表示使目标函数取最小值时的变量值，μ₁、μ₂和μ₃表示正则化系数，R_i表示取块矩阵；x_P、x_S和x_ρ分别表示纵波速度、横波速度和密度；D_p表示由训练样本学习得到的纵波速度字典，D_s示由训练样本学习得到的横波速度字典，D_ρ表示由训练样本学习得到的密度字典；

为x_P经D_p稀疏表示的稀疏系数，

为x_S经D_s稀疏表示的稀疏系数，

为x_ρ经D_ρ稀疏表示的稀疏系数。

步骤C具体为：

C1、采用坐标下降策略求解步骤B得到的目标函数，得到三个子问题：

子问题(a)为基于零阶Tikhonov正则化约束的地震反演问题；

k表示迭代次数；

子问题(b)是一个稀疏逼近问题；

子问题(c)为一个常规最小化问题；

C2、通过迭代求解三个子问题，得到最终的模型参数。

步骤C2中采用梯度下降算法求解子问题(a)。

步骤C2中所述求解子问题(b)，具体为：采用正交字典对数据进行稀疏表征，从而迭代求解稀疏表征与字典更新，根据迭代收敛的正交字典求解得到稀疏系数。

本发明的有益效果：本发明在保证与KSVD字典方法具有相当反演效果的同时，能有效提升稀疏反演速率；并且相比于KSVD字典方法需要考虑字典大小与稀疏度多个参数对反演效果的影响下，本发明方法无须过多考虑字典大小对反演的影响，本发明主要通过调节硬阈值一个参数来控制反演的效果，噪声测试与实际资料测试也验证了本发明所提出的方法的有效性。综上，本发明提出的方法能有效应对复杂的地震反演问题，并且运算效率更高，参数调节更加便捷，能够有效促进字典学习和稀疏表征理论在处理大规模地震数据反演领域的发展与应用。

附图说明

图1为本发明的方案流程图；

图2(a)-图2(i)为本发明实施例提供的本发明方法与现有的KSVD方法在4dB高斯噪声环境下的叠前三参数反演结果对比图；

其中，图2(a)为真实纵波速度，图2(b)为横波速度，图2(c)为密度剖面，图2(d)为KSVD字典方法的纵波速度反演结果，图2(e)为KSVD字典方法的横波速度反演结果，图2(f)为KSVD字典方法的密度剖面反演结果，图2(g)为正交字典方法的纵波速度反演结果，图2(h)为正交字典方法的横波速度反演结果，图2(i)为正交字典方法的密度剖面反演结果；

图3(a)-图3(f)为本发明实施例提供的本发明方法与现有的KSVD方法在国内某工区叠前实际资料反演结果对比图；

其中，图3(a)为KSVD字典约束反演的纵波，图3(b)为KSVD字典约束反演的横波，图3(c)为KSVD字典约束密度反演剖面，图3(d)为ORTD字典约束反演的纵波，图3(e)为ORTD字典约束反演的横波，图3(f)为ORTD字典约束密度反演剖面；

图4(a)-图4(f)为本发明实施例提供的国内某工区叠前井道处反演结果对比图；

其中，图4(a)为KSVD字典约束反演的纵波，图4(b)为KSVD字典约束反演的横波，图4(c)为KSVD字典约束密度反演剖面，图4(d)为ORTD字典约束反演的纵波，图4(e)为ORTD字典约束反演的横波，图4(f)为ORTD字典约束密度反演剖面。

具体实施方式

为便于本领域技术人员理解本发明的技术内容，下面结合附图对本发明内容进一步阐释。

如图1所示，本发明的基于快速正交字典的稀疏表征正则化叠前AVO反演方法，包括以下步骤：

S1、提取地震子波。本发明采用确定性子波提取方法，确定性子波提取方法是在有测井数据的情况下，利用测井资料计算出的反射系数序列，然后基于褶积模型结合井旁地震道，采用最小平方法估算出地震子波的方法；

S2、井震标定。井震标定是联系地震数据和测井数据的桥梁，其方法是由原始井旁地震道来提取相应时间窗口范围内的地震子波，然后再利用子波褶积得到的合成井旁地震道和原始井旁地震道进行对比；

S3、建立初始模型。建立初始模型首先要根据地震层位解释信息，按照一定地质沉积规律在解释层位间细分出多个小层，搭建出相对精细的地质框架结构。然后在精细地质框架模型下，运用合适的数学插值方法(克里金算法等)，在单井井震联合标定的基础上，将测井数据沿地震解释层位进行外推和内插，构建平滑、闭合的初始模型；

S4、目标函数构建。在基于字典学习方法的叠前AVO反演方程的基础上，采用正交字典学习与稀疏表征技术。步骤S4的详细过程如下：

S41、基于字典学习方法的叠前AVO反演方程

本发明的AVO反演基于Fatti于1994年改进的Aki&Richards近似方程：

R_PP(θ)＝c₁R_P+c₂R_S+c₃R_ρ (1)

式中：R_pp为纵波反射系数，c₁＝1+tan²θ，c₂＝-8γ²tan²θ，c₃＝-0.5tan²θ+2γ²sin²θ；R_p表示纵波阻抗反射系数，R_s表示横波阻抗反射系数，R_ρ表示密度反射系数；θ为纵波入射角。

对于任意入射角的叠前角道集数据，对应的反演模型可以表示成如下形式：

式中，s(θ)表示任意入射角度的角道集地震数据，L_P表示纵波波阻抗的自然对数，L_S表示横波波阻抗的自然对数，L_ρ表示密度的自然对数，W(θ)表示子波矩阵，D表示导数矩阵。

叠前角道集反演模型需要求解纵波速度、横波速度和密度三个弹性参数，Hampson(2005)经过研究发现纵波阻抗的对数、横波阻抗的对数和密度的对数之间存在着如下线性关系：

其中，i表示解释层位序号；

和

分别表示纵横波波阻抗的对数和密度的对数；

和

是线性拟合后的残差；a,b,c,d为常线性拟合系数，可以从测井曲线中学习得到。因此，当纵波速度、横波速度和密度的初始值已知时，可以很容易地转换成模型参数：

根据上述线性关系，将叠前角道集反演的输入设置成如下模型参数的形式，n代表采样点的数量：

将式(5)代入到式(2)的叠前反演模型，可表示为：

当已知角道集地震数据、子波序列、入射角信息和初始模型参数，便可以通过式(6)构建目标函数得到对应的模型参数，进而转换为纵波速度、横波速度和密度。

叠前反演的模型可以表示成s＝Gx，其中s表示叠前角道集地震数据，G表示正演矩阵。反演过程可以被描述为最小二乘最小化问题。由于这个反问题是病态的，因此通常采用正则化方案来求解该问题。加入正则化约束的反演目标函数为：

式(7)中第一项

用于计算反演过程的均方误差，第二项函数R(x)是关于模型参数的正则化项，argmin表示使目标函数取最小值时的变量值。基于字典学习的稀疏表征地震反演也是通过在反演目标函数中添加正则化项进行先验信息约束，其对应的目标函数表示为：

其中，x＝[L_P△L_S△L_ρ]表示模型参数；μ₁、μ₂和μ₃表示正则化系数，代表了地震约束项和稀疏字典约束之间的相对权重大小；R_i表示取块矩阵；x_P、x_S和x_ρ分别表示纵波速度、横波速度和密度；D_P表示由训练样本学习得到的纵波速度字典，D_S表示由训练样本学习得到的横波速度字典，D_ρ表示由训练样本学习得到的密度字典；

为x_P经D_P稀疏表示的稀疏系数，

为x_S经D_S稀疏表示的稀疏系数，

为x_ρ经D_ρ稀疏表示的稀疏系数。

对于上述叠前AVO反演算法的具体求解过程，本发明为了解决这个复杂的方程(8)，考虑添加中间变量y，得到式(9)的表达式，然后再对式(9)进行恒等变换：

其中，y表示对模型参数x的约束，其值等于模型参数的初始值；y_P、y_S和y_ρ分别表示纵波速度、横波速度和密度；β为正则化系数。

求解式(9)通常采用坐标下降策略(Coordinate Descent Strategy)每次只更新一个变量，可以将上面的方程分成以下三个子问题的交替求解：

k表示迭代次数；

其中子问题(a)将其看成是基于零阶Tikhonov正则化约束的地震反演问题，可以利用梯度下降算法求解；问题(b)是一个稀疏逼近问题，问题(c)为一个常规最小化问题。通过迭代三个子问题来解决原来的问题(9)，直到满足一定的收敛准则。

其中问题(b)的求解需要通过学习得到的字典来对数据进行稀疏表征从而求解得到稀疏系数，不同字典学习方法具有不同的稀疏系数求解策略。本发明采用正交字典学习方法，问题(b)的求解过程采用硬阈值操作，有关正交字典学习方法的详细介绍如下。

S42、正交字典学习

相比于基于解析变换的字典而言，从一组训练样本中学习的字典可以实现更好的稀疏表示，例如KSVD算法(Aharon et al，2006)作为一种众所周知的字典学习方法，它采用贪婪算法(如：正交匹配追踪(OMP))和奇异值分解(SVD)进行稀疏编码和字典更新过程，其在训练字典前需要提前设置好字典原子长度、原子个数以及稀疏度三个参数，并且由于KSVD字典学习通常基于冗余假设，因此其生成的字典原子个数通常大于原子长度(非方阵)。基于学习的过完备字典虽然比解析变换字典具有更好的稀疏逼近特性(Rubinsteinet al,2010)，但有研究表明过度冗余的字典并不一定具有最佳的表征能力(Ravishankarand Bresler，2011)，并且过完备字典中的原子由于缺少相关性约束，造成了字典学习过程迭代的增加和小块冗余重叠，因此过完备的字典学习具有更高的计算复杂度。

本发明采用(Bao et al，2013)提出的正交字典学习方法，其字典原子之间存在正交性约束(字典原子个数等于原子长度(方阵))，在训练字典过程只需提前设置好字典原子长度以及阈值系数两个参数。相比于KSVD字典该方法计算速度快，性能好，参数少，具有以下目标函数：

s.t.D^TD＝I_n

其中，y_i表示用于训练的样本数据块，

v_i表示y_i对应的稀疏系数，D表示由学习得到的n个原子组成的正交字典。采用正交字典将大大简化字典更新和稀疏编码的计算。

与传统的KSVD方法在交替极小化框架下学习过完备字典一样，正交字典学习也采用交替迭代两个子问题：稀疏编码和字典更新求解式(13)。具体而言，首先随机初始化正交字典D⁽⁰⁾，对于k＝0,1,…,K-1,进行以下交替迭代：

1、稀疏编码：给定正交字典D^(k)，通过求解式(14)得到稀疏系数V^(k)

2、字典更新：给定稀疏系数V^(k)，通过求解式(15)来更新字典；

s.t.D^TD＝I_n

由于字典D具有正交性特点，在进行稀疏编码和字典更新时，每一次迭代都不需要解决任何最小化问题且具有更低的计算复杂度，如表1展示了KSVD字典与本发明的基于正交字典训练迭代一次过程的理论计算复杂度公式(Bao et al，2013)，表1中m、n分别表示字典原子大小和原子个数，K表示稀疏系数稀疏度。其中稀疏编码通过硬阈值操作完成，T_λ(·)表示硬阈值算子，字典更新通过单个SVD完成。有关算法的完整描述，详见后文的算法1。

表1迭代一次KSVD字典与正交字典训练过程计算复杂度(Bao et al，2013)

其中，符号“：＝”表示定义为。

S5、选取合适的正则化参数。如果正则化参数选取得过小，则反演的不适定性得不到较好地改善，正则化方法无法发挥作用；如果正则化参数选取得过大，则会出现反演结果过拟合而偏离真实情况的现象，需要根据工区盲井测试所得误差最小化结果进行调整测试。

S6、选取合适的硬阈值系数。如果硬阈值系数选取得过小，则容易造成反演结果的过拟合效应；如果硬阈值系数选取的过大，则容易造成反演结果的欠拟合效应，需要根据工区盲井测试所得误差最小化结果进行调整测试。

S7、迭代更新模型。从初始模型出发，利用正演合成记录与叠前地震记录的残差计算出的模型修正量来不断更新模型，使模型正演合成记录与实际地震资料达到最佳吻合，最后输出反演结果。

以下结合图2-4和具体实施例来说明本发明的技术效果：

实施例1

为了评价本发明的基于正交字典稀疏促进正则化约束的叠前AVO反演速率，并检验其反演效果。通过采样和时深转换处理后，选取了时间域的大小为200×400，采样间隔为4ms的Marmousi模型以及主频25Hz雷克子波来测试所提出的算法，如图2(a)所示为真实纵波速度，如图2(b)所示为横波速度，如图2(c)所示为密度剖面。与基于冗余假设的KSVD字典稀疏约束的反演方法一致，我们从Marmousi模型400道数据中等间距的选取相同的20道井数据取小块后作为训练样本。本是实施例中用于反演约束的纵波速度、横波速度和密度的正交字典分别为D_P、D_S和D_ρ，字典大小均设为50×50，分别由对应的训练样本经算法1学习得到。

如图2(a)-图2(i)为Marmousi模型加入4dB高斯噪声环境下的叠前三参数反演结果对比图，其中图2(d)为KSVD字典(字典大小为50×1500)方法的纵波速度反演结果，图2(e)为KSVD字典(字典大小为50×1500)方法的横波速度反演结果，图2(f)为KSVD字典(字典大小为50×1500)方法的密度剖面反演结果，图2(g)为正交字典方法的纵波速度反演结果，图2(h)为正交字典方法的横波速度反演结果，图2(i)为正交字典方法的密度剖面反演结果。从图2中可以看出，相比于KSVD字典，本发明提出的正交字典约束下的反演结果其纵波速度，横波速度反演效果相当，但密度反演效果有明显提升。同时对比测试了两种反演算法的训练字典耗时及反演耗时情况，其结果如表2所示，采用本发明的正交字典在保证反演效果的同时，体现出了明显的效率优势。

表2Marmousi模型中KSVD方法与正交字典(ORTD)方法运行时间比较

字典类型	字典大小	训练字典耗时	叠前Marmousi剖面反演耗时
				KSVD字典	50×1500	43.92秒	629.17秒
正交字典	50×50	0.16秒	224.91秒

实施例2

为测试所提出方法的实际应用性，将本发明的方法应用到国内某工区的叠前实际资料反演中，该实际资料时间域的大小为80×853，其中训练字典过程采用该工区的13口井数据经低通滤波后取小块组成训练样本。如图3(a)-图3(f)所示为基于两种字典的反演结果对比图，其中图3(a)为KSVD字典(字典大小为50×700)约束反演的纵波，图3(b)为KSVD字典(字典大小为50×700)约束反演的横波，图3(c)为KSVD字典(字典大小为50×700)约束密度反演剖面，图3(d)为ORTD字典(字典大小为50×700)约束反演的纵波，图3(e)为ORTD字典(字典大小为50×700)约束反演的横波，图3(f)为ORTD字典(字典大小为50×700)约束密度反演剖面。从图3(a)-图3(f)中可以看出本发明提出的方法(ORTD)其反演结果的横向连续性，分辨率都稍好于KSVD反演结果。为更直观对比其反演效果，选取了井道处数据进行盲井测试(该井未参与训练字典)，图4(a)为KSVD字典(字典大小为50×700)约束反演的纵波，图4(b)为KSVD字典(字典大小为50×700)约束反演的横波，图4(c)为KSVD字典(字典大小为50×700)约束密度反演剖面，图4(d)为ORTD字典(字典大小为50×700)约束反演的纵波，图4(e)为ORTD字典(字典大小为50×700)约束反演的横波，图4(f)为ORTD字典(字典大小为50×700)约束密度反演剖面；图4(a)-图4(f)的定量化对比结果如表3，可以发现本发明的方法反演结果精度与KSVD方法相当，因此本发明方法能有效应用于实际资料反演中。

表3井道处反演结果均方根误差(RMSE)对比

我们用叠前Marmousi模型作为参考模型，生成叠前地震记录，验证本发明的有效性。已知真实的叠前Marmousi纵横波速度、密度剖面，如图2(a)、图2(b)、图2(c)所示。我们采用叠前Marmousi模型与主频为25Hz的Ricker子波褶积并加入4dB的高斯噪声得到合成地震记录进行反演测试。如图2(d)、图2(e)、图2(f)为KSVD字典方法(KSVD)叠前三参数反演结果。图2(g)、图2(h)、图2(i)为正交字典方法(ORTD)叠前三参数反演结果。从图2(a)-图2(i)中可以发现，两种方法反演结果其纵波速度，横波速度反演效果相当，但正交字典方法密度反演效果有明显提升。同时对比测试了两种反演算法的训练字典耗时及反演耗时情况，其结果如表2所示，发现采用正交字典方法(ORTD)在保证反演效果的同时，并体现出明显的效率优势。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的权利要求范围之内。

Claims

1.一种基于快速正交字典的稀疏表征正则化叠前AVO反演方法，其特征在于，包括：

B、加入正则化约束；步骤B所述加入正则化约束的目标函数为：

其中，

是用于计算反演过程的均方误差，函数R(x)是关于模型参数x的正则化项，argmin表示使目标函数取最小值时的变量值，μ₁、μ₂和μ₃表示正则化系数，R_i表示取块矩阵；x_P、x_S和x_ρ分别表示纵波速度、横波速度和密度；D_P表示由训练样本学习得到的纵波速度字典，D_S表示由训练样本学习得到的横波速度字典，D_ρ表示由训练样本学习得到的密度字典；