一种基于Hesssian矩阵约束的多道双极稀疏谱反演方法与流程

一种基于hesssian矩阵约束的多道双极稀疏谱反演方法
技术领域
1.本发明属于油气及地球物理勘探领域，具体涉及一种基于hessian矩阵约束的多道双极稀疏谱反演方法。


背景技术：

2.由于受调谐效应的影响，应用常规地震资料及传统处理方法无法识别厚度小于调谐厚度的薄层和小尺度地质构造，难以提供高精度的解释成果。谱反演方法不依赖测井资料等先验信息，只需利用地震数据本身的频谱信息就能反演地层稀疏反射系数。由于其能够分辨调谐厚度以下的地震薄层，该方法得到了广泛的研究和讨论。
3.传统的谱反演方法采用单条地震道全局求取频谱，并对目标函数添加1范数约束项进行求解。对单条地震道全局求取频谱，可能造成部分频率域信号重叠，无法凸显局部的有效信息，同时忽略了地震波传播过程中的非平稳性，容易造成误差。
4.对单条地震道目标函数添加1范数约束项，能够对反射系数序列进行稀疏约束，但由于多条地震道分别计算，算法收敛迭代次数不同，各道反射系数没有空间连贯性，在反射系数剖面上出现“错位”、“分叉”等现象，无法反映真实的地层反射系数分布情况。
5.中国专利申请cn111856559a公开了一种基于稀疏贝叶斯学习理论的多道地震谱反演方法，所述方法包括：获取经过预处理的叠后地震资料；从所述叠后地震资料中提取多道地震谱反演参数；从所述叠后地震资料中提取地震子波；基于傅里叶变换提取所述地震子波的频谱信息；根据多道地震谱反演参数构建多道地震谱反演的正演算子矩阵；基于傅里叶变换从所述叠后地震资料中提取多道平均地震道频谱信息；基于所述频谱信息构建多道地震谱反演目标函数，并通过基于最大期望算法的稀疏贝叶斯学习理论求解目标函数以获取多道地震谱反演结果。
6.中国专利申请cn111796325a公开了一种分频迭代约束的随机反演方法，包括：优选反演需要的已钻井；精细井震标定；拓展现有地震资料频带；对拓频后宽频地震资料进行分频处理；高频子波提取；获取三维视砂地比分布体；地震资料驱动为主的反演方法；提升地震资料驱动为主导的随机反演结果的精度，将获得的三维视砂地比分布体作为先验信息加入到高频地震资料驱动的随机反演过程中，从而获得高精度的分频迭代约束随机反演结果；该方法有效提高储层预测精度、更准确地刻画储层横向厚度变化、恢复储层真实的叠置连通关系，为勘探目标评价阶段的钻前井位部署和油田综合调整阶段的井位优化提供重要的参考依据。
7.然而，基于各道反射系数没有空间连贯性，在反射系数剖面上出现“错位”、“分叉”等现象，无法反映真实的地层反射系数分布情况的问题，目前仍未有有效的解决方法。


技术实现要素：

8.本发明主要目的在于提供一种基于hessian矩阵约束的多道双极稀疏谱反演方法，本发明方法在提高地震数据分辨率的同时，还能够保证多道数据中地层的横向连续性，
能够提供更为准确、丰富的地质和储层信息，克服了上述现有技术的不足。
9.为实现上述目的，本发明采用以下技术方案：
10.本发明提供一种基于hessian矩阵约束的多道双极稀疏谱反演方法，包括以下步骤：
11.估计地震子波波长确定时窗长度，将多道地震数据分割；提取时窗内地震子波；在时窗内构建基于hessian矩阵约束的多道双极稀疏谱反演目标函数，计算该时窗内多道地震数据的反射系数；将各时窗中心反射系数结果沿时间顺序组合，得到多道地震数据的谱反演结果。
12.进一步地，估计地震子波波长确定时窗长度的具体方法：
13.估计多道地震数据s(x,t)的子波主频fm：
[0014][0015]
其中s(x,f)表示多道信号s(x,t)的振幅谱，m为地震道数；
[0016]
取主频为fw零相位雷克子波长度作为时窗的窗长，使时窗中心采样点对应零相位雷克子波的主峰位置。
[0017]
进一步地，将多道地震数据分割的方法为：以一个时间采样间隔为移窗步长，利用时窗分割拓展后的各道地震数据；各时窗中心点与原始数据采样点一一对应。
[0018]
进一步地，在时窗内使用统计子波估计方法，提取具有时变特性的地震子波：
[0019]
w(t)＝[1-2(πfmt)2]exp[-(πfmt)2]
[0020]
其中，fm为所在时窗内估计的雷克子波主频。
[0021]
进一步地，构建基于hessian矩阵约束的多道双极稀疏谱反演目标函数：
[0022]
多道地震褶积模型在时间域的褶积运算可以变换为频率域的乘积运算：
[0023][0024]
单个时窗内可以构建谱反演的目标函数：
[0025]
obj＝w(f)
×
r(x,f)-s(x,f)
[0026][0027]
其中，
[0028]
a为子波频谱通过奇偶分解构造的矩阵：
[0029][0030]a11
＝{αore[w(fi)]sin(m
i,k
)sin(ni)-αoim[w(fi)]sin(m
i,k
)cos(ni)}
[0031]a12
＝{αere[w(fi)]cos(m
i,k
)cos(ni) αeim[w(fi)]cos(m
i,k
)sin(ni)}
[0032]a21
＝{αore[w(fi)]sin(m
i,k
)cos(ni) αoim[w(fi)]sin(m
i,k
)sin(ni)}
[0033]a22
＝{-αere[w(fi)]cos(m
i,k
)sin(ni) αeim[w(fi)]cos(m
i,k
)cos(ni)}
[0034]mij
＝πfitj[0035]
ni＝πfi(n-1)δt
[0036]
b为多道地震记录道频谱通过奇偶分解构造的矩阵：
[0037][0038]ro
、re分别为多道反射系数奇、偶分量组成的矩阵，通过奇偶合并能够得到多道反射系数矩阵r(x,t)：
[0039][0040]
αo、αe分别为奇、偶加权系数，n为时窗内采样点数，m为道数，h(r)为反射系数r的hessian矩阵，||g||
2，1
为l2,1范数；||g||f为frobenius或euclidean范数。
[0041]
进一步地，利用最优解算法在时窗内求解谱反演目标函数的最小值。基于hessian矩阵约束的多道双极稀疏谱反演目标函数的结构为典型的线性回归方程，使用常规的最小二乘方法进行迭代拟合可以得到准确的谱反演结果。同时，应用遗传算法、模拟退火法、粒子群算法等全局优化算法也可以得到准确的谱反演结果。
[0042]
进一步地，分割的所有时窗均要计算完毕。
[0043]
进一步地，将各时窗中心点的反射系数值沿时间顺序，最终得到理想的稀疏谱反演结果：
[0044]
r(x,t)＝[r1(x,t)
·
δ(x,t) l r
p
(x,t)
·
δ(x,t)]m×
p
[0045]
其中r
p
(x,t)是第p个时窗内谱反演计算多道反射系数结果，δ(x,t)为狄拉克函数矩阵，m为道数，p为时窗个数。记录所得多道反射系数序列在时窗中心采样点的各值，即对时窗内谱反演多道反射系数结果加狄拉克函数窗，能够消除由时窗分割导致的频谱泄漏影响。
[0046]
与现有技术相比，本发明具有以下优势：
[0047]
本发明方法对多条地震道添加时窗，截取局部信号求取频谱同时计算，时窗内信号近似平稳，更能凸显其局部频谱特征，提取出的子波在具有空变特性的同时具有时变特性。多道同时计算，能够保证各道在计算时参数具有一致性。
[0048]
本发明方法还对谱反演目标函数添加了反射系数的hessian矩阵进行结构约束，hessian矩阵元素由反射系数关于空间二阶偏导数组成，可用于检测三维空间界面奇异性，hessian矩阵的frobenius或euclidean范数可以确定与层边界和断层相对应的超曲面，hessian矩阵约束就是迫使反演出的反射系数沿着3d空间平滑曲面分布，也即反射系数在各道之间具有较好的横向连续性，更符合真实地层反射系数的分布特征。
[0049]
本发明方法能将由于相干形成的复波同相轴分离开来，识别原始地震剖面上无法分辨的、厚度小于调谐厚度的薄层和小尺度地质体，反映地层真实的接触关系，刻画地层真实的尖灭位置，提供高分辨率的、更加丰富的地质信息。在提高地震数据分辨率的同时，该方法能够保证多道数据中地层的横向连续性，能够提供更为准确、丰富的地质和储层信息；
因而，本发明方法可在能够在储层预测和油气检测中发挥重要作用。
附图说明
[0050]
构成本发明的一部分的说明书附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。
[0051]
图1为本发明一具体实施例所述基于hessian矩阵约束的多道双极稀疏谱反演方法流程图；
[0052]
图2为本发明一具体实施例所述砂泥互层楔形模型地震正演剖面；
[0053]
图3为本发明一具体实施例所述砂泥互层楔形模型地震正演剖面使用基于hessian矩阵约束的多道双极稀疏谱反演方法得到的反射系数剖面；
[0054]
图4为本发明一具体实施例所述1砂泥互层楔形模型地震正演剖面使用常规谱反演方法得到的反射系数剖面；
[0055]
图5为本发明一具体实施例所述marmousi模型地震正演剖面；
[0056]
图6为本发明一具体实施例所述marmousi模型地震正演剖面使用基于hessian矩阵约束的多道双极稀疏谱反演方法得到的反射系数剖面；
[0057]
图7为本发明一具体实施例所述marmousi模型地震正演剖面使用常规谱反演方法得到的反射系数剖面。
具体实施方式
[0058]
应该指出，以下详细说明都是示例性的，旨在对本发明提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本发明所属技术领域的普通技术人员通常理解的相同含义。
[0059]
需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本发明的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，当在本说明书中使用术语“包含”和/或“包括”时，其指明存在特征、步骤、操作和/或它们的组合。
[0060]
为了使得本领域技术人员能够更加清楚地了解本发明的技术方案，以下将结合具体的实施例详细说明本发明的技术方案。
[0061]
实施例1
[0062]
如图1所示，所述基于hessian矩阵约束的多道双极稀疏谱反演方法包括以下步骤：
[0063]
s1、输入多道地震数据。
[0064]
s2、估计地震子波以确定时窗窗长，利用时窗分割数据：
[0065]
首先，估计多道地震数据s(x,t)的子波主频fm：
[0066][0067]
其中s(x,f)表示多道信号s(x,t)的振幅谱，m为地震道数。
[0068]
取主频为fw零相位雷克子波长度作为时窗的窗长，使时窗中心采样点对应零相位
雷克子波的主峰位置。
[0069]
之后，将移窗步长设为多道地震数据的一个时间采样间隔，利用时窗分割各道地震数据。此时，各时窗中心点与原始数据各采样点一一对应。
[0070]
s3、在时窗内使用统计子波估计方法，提取具有时变特性的地震子波：
[0071]
w(t)＝[1-2(πfmt)2]exp[-(πfmt)2]，
[0072]
其中，fm为所在时窗内估计的雷克子波主频。
[0073]
s4、构建并求解基于hessian矩阵约束的多道双极稀疏谱反演目标函数：
[0074]
多道地震褶积模型在时间域的褶积运算可以变换为频率域的乘积运算，
[0075][0076]
在单个时窗内可以构建谱反演的目标函数：
[0077]
obj＝w(f)
×
r(x,f)-s(x,f)
[0078][0079]
其中，
[0080]
a为子波频谱通过奇偶分解构造的矩阵：
[0081][0082]a11
＝{αore[w(fi)]sin(m
i,k
)sin(ni)-αoim[w(fi)]sin(m
i,k
)cos(ni)}
[0083]a12
＝{αere[w(fi)]cos(m
i,k
)cos(ni) αeim[w(fi)]cos(m
i,k
)sin(ni)}
[0084]a21
＝{αore[w(fi)]sin(m
i,k
)cos(ni) αoim[w(fi)]sin(m
i,k
)sin(ni)}
[0085]a22
＝{-αere[w(fi)]cos(m
i,k
)sin(ni) αeim[w(fi)]cos(m
i,k
)cos(ni)}
[0086]mij
＝πfitj[0087]
ni＝πfi(n-1)δt
[0088]
b为多道地震记录道频谱通过奇偶分解构造的矩阵：
[0089][0090]ro
、re分别为多道反射系数奇、偶分量组成的矩阵，通过奇偶合并能够得到多道反射系数矩阵r(x,t)：
[0091][0092]
αo、αe分别为奇、偶加权系数，n为时窗内采样点数，m为道数，h(r)为反射系数r的hessian矩阵，||g||
2，1
为l2,1范数；||g||f为frobenius或euclidean范数。
[0093]
基于hessian矩阵约束的多道双极稀疏谱反演目标函数的结构为典型的线性回归方程，使用常规的最小二乘方法进行迭代拟合求解谱反演目标函数的最小值可以得到准确的谱反演结果。同时，应用遗传算法、模拟退火法、粒子群算法等全局优化算法也可以得到准确的谱反演结果。
[0094]
s5、重复步骤s3、s4，直到步骤s2中分割的所有时窗计算完毕。
[0095]
s6、将各时窗中心点的反射系数值沿时间顺序，最终得到理想的稀疏谱反演结果：r(x,t)＝[r1(x,t)
·
δ(x,t) l r
p
(x,t)
·
δ(x,t)]m×
p
[0096]
其中r
p
(x,t)是第p个时窗内谱反演计算多道反射系数结果，δ(x,t)为狄拉克函数矩阵，m为道数，p为时窗个数。记录所得多道反射系数序列在时窗中心采样点的各值，即对时窗内谱反演多道反射系数结果加狄拉克函数窗，能够消除由时窗分割导致的频谱泄漏影响。
[0097]
实施例2
[0098]
对砂泥互层楔形模型地震正演剖面使用实施例1所述基于hessian矩阵约束的多道双极稀疏谱反演方法进行反演，所述砂泥互层楔形模型地震正演剖面如图2所示，计算反射系数剖面如图3所示。
[0099]
实施例3
[0100]
对marmousi模型地震正演剖面使用实施例1所述基于hessian矩阵约束的多道双极稀疏谱反演方法进行反演，所述marmousi模型地震正演剖面如图5所示，计算反射系数剖面如图6所示。
[0101]
将图3、图6所示反演结果与常规谱反演结果(图4、图7)相比，可见：常规谱反演对薄层的分辨能力相对较差，当薄层厚度小于1/8波长时无法分辨，同时反射系数位置发生了错动，缺乏横向连续性；基于hessian矩阵约束的多道双极稀疏谱反演能够分辨厚度小于1/8波长的薄层，反射系数在各道之间具有较好的横向连续性；本发明的谱反演结果具有更高的薄层分辨能力，更符合真实地层反射系数的分布特征，具有较高的应用价值。
[0102]
上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。