本发明属于地震勘探，具体涉及一种基于k1峭度约束的地震频域反褶积方法、系统、设备及存储介质。
背景技术：
1、盆地的地表大多是被黄土塬、沙漠、戈壁滩等恶劣条件覆盖，且地震波在传播过程中存在吸收衰减等效应，导致采集到的反射地震信号波形展宽、主频下降，纵向分辨率较差，从而影响后续储层结构特征及内部流体特性的精细描述。因此，地震资料高分辨处理是盆地地震资料处理中一个基础和关键的技术。
2、地震资料高分辨率处理方法的方法有很多，主要包括谱白化方法(陈传仁等)、反q滤波(wang)、反褶积(margrave et al.,2011)等。其中，反褶积方法通过压缩地震子波获取地层反射系数，是地震资料处理中提高纵向分辨率、细致刻画薄储层最常用的方法之一。假设地震反射系数序列的稀疏性是目前地震反褶积常用的方法之一。根据褶积模型和地震反射系数序列的稀疏性假设，地球物理学家们提出了一系列基于反演框架的地震反褶积方法，taylor et al.(1979)利用交替迭代求解地震子波与反射系数对应的反问题，来同时估计地震子波与反射系数。在taylor等工作基础上，wang et al.,(2016)提出一种基于稀疏toeplitz矩阵分解的稀疏盲反褶积方法。上述算法大多是假设地震子波是未知的来实现的。如果地震子波已知，根据褶积模型的原理，可以将上述基于迭代优化算法的反褶积算法进行简化，如脉冲反褶积(潘树林等,2019)等。近年来，随着深度学习的快速发展，该算法也被引入基于反演框架的地震反褶积算法中(chai et al.,2021；chen et al.2023；yang etal.2024)。
3、同时，假设地震反射系数序列是满足独立同分布的随机变量是另外一种常用的反褶积方法。最初的反褶积方法通常假设地层反射系数序列是服从独立同分布的高斯分布，然后利用二阶统计量来实现反褶积算法(treitel and robinson,1964)。但是，实际地层的反射系数序列并不严格满足高斯分布，限制了上述方法的应用范围(蔡连芳和田学民,2012；luo.,2024)。随着地球物理学家们对反射系数序列统计规律的深入研究，发现实际地震反射系数序列更适合用服从独立同分布的非高斯分布来模拟(wang et al.,2018)。wiggins(1978)提出了最小熵反褶积算法，该算法利用最大化地震数据的四阶统计量(即峭度)来估计逆滤波器，从而实现反射系数序列的非高斯性度量。donoho(1981)进一步发展了最小熵反褶积算法，利用概率和信息论等理论阐述了盲反褶积的一般框架。taleb et al.(2001)提出了基于互信息熵的反褶积算法，该算法利用输出信号的互信息熵来调整逆滤波器，从而获得较优的反褶积结果。larue et al.(2006)提出了频率域基于最小化互信息率的反褶积算法。在此工作基础上，baan和pham(2008)对其进行推广，在子波的频带内对数据进行反褶积处理，同时引入维纳滤波器来有效地抑制噪声。此外，地球物理学家们也探讨了rényi熵的性质和度量方法，并提出了最小rényi熵的盲反褶积算法(nadarajah et al.,2003；)。上述方法都是通过评估输出结果的统计特征(如峭度、互信息率、熵等)来判断反褶积的反射系数序列是否满足非高斯分布的。这类方法虽然也可以得到较好的结果，但在一些复杂地形中所得到的地震分辨率仍不足以满足勘探需要。主要存在如下缺点：
4、(1)对反射系数序列的稀疏性假设过于严格，而实际地震反射系数序列过于复杂，难以用稀疏性假设来获得准确的反射序列。
5、(2)针对实际地震数据，统计分布特征比高斯分布更加尖锐，反射稀疏序列的非高斯性无法被很好地表示。在复杂地形中，地震子波类型多变不容易匹配。
技术实现思路
1、为了克服上述现有技术的缺点，本发明的目的在于提供一种基于k1峭度约束的地震频域反褶积方法、系统、设备及存储介质，能够有效解决现有技术中难以用稀疏性假设来获得准确的反射序列，以及复杂地形中地震子波类型多变不容易匹配的技术难题。
2、为了达到上述目的，本发明采用以下技术方案予以实现：
3、本发明公开的一种基于k1峭度约束的地震频域反褶积方法，首先，引入广义beta小波，使其能够匹配各种类型的地震子波。然后，利用k1峭度来更加准确地度量地震反射系数序列的非高斯性，从而提高地震数据的纵向分辨率。
4、本发明公开的基于k1峭度约束的地震频域反褶积方法，具体包括以下步骤：
5、1)对地震数据预处理，引入广义beta母小波函数，并估算对应的地震子波的主频；
6、2)假设其参数α已知，计算相应的地震子波频率表达式逆滤波器表达式
7、3)通过k1峭度估计出最优的子波，从而使得反射系数更满足非高斯性；
8、4)通过最优的地震子波wopt(t)，得到最优的wiener滤波器aopt(ω)并计算出最后的反褶积结果。
9、优选地，步骤1)具体操作包括：
10、引入广义beta母小波函数可以更好的匹配地震子波。广义beta小波的频率域定义为：
11、
12、式中，ω为角频率；α和β是广义beta小波的波形参数，调整这两个参数可以控制广义beta母小波的形态，使其满足不同类型的地震子波；μα,β是归一化常数；u(ω)是单位阶跃信号，保证广义beta小波是解析的。通过调整两个参数α和β，广义beta小波变换表现出很大的灵活性。
13、同时，在处理实际地震数据时，地震数据的主频可以使用如下公式计算：
14、
15、优选地，步骤2)中，具体操作为：
16、首先设置最大迭代次数为nmax和α＝{αk},k＝1,2,....,nmax的范围；
17、其次根据式(1)广义beta小波的定义，可以推导出广义beta母小波的主频与波形参数α、β的关系，表示为：
18、
19、式中，ωp是广义beta母小波的主频。假设参数α已知，在主频可以通过式(2)计算的情况下可以求得震子波的参数β。由此地震子波频率表达式可以通过式(1)求得。
20、最后逆滤波器表达式可以通过下式求得：
21、
22、式中，wa(ω)表示子波wa(t)的频率表达，为wa(ω)的复共轭。ξ表示wiener滤波器的阻尼因子。该滤波器为wiener滤波器，wiener滤波器被认为是在最小二乘假设下的最优滤波器。
23、优选地，步骤3)中，具体操作为：
24、首先，当已知地震记录x(t)和wiener滤波器时，可以得到估计的反射系数频率表达式为
25、
26、其次，为了让反射系数更加满足非高斯性，用k1峭度来度量估计出的反射系数即
27、
28、式中，j(α)是反射系数的k1峭度估计；α是待优化的广义beta子波波形参数，通过优化公式(6)可以获得最优的子波，从而使得反射系数更满足非高斯性。优化公式可以表示为进一步，为了计算方便，令b＝e，公式(6)可以改写为：
29、
30、最后求解上述优化公式，可以得到最优的α，从而得到最优的子波wopt(t)。
31、进一步优选地，步骤3)中，基于k1峭度来度量估计出的反射系数，其中：
32、峭度是描述数据集位置及变异度的统计量，其定义是归一化四阶中心矩。传统峭度k的定义为：
33、
34、式中，z表示随机变量；e(z)表示随机变量z的均值；var(z)表示随机变量z的方差。当随机变量均值为0，且方差为单位方差时，峭度k可以重新表示为：k＝e(z)4。利用峭度的峰值可以表示信号的脉冲分量，如果峭度峰值越高，信号强度分布也越锐利。与正态分布对比，峭度可以用来度量一个分布是轻尾或者是重尾。
35、为了利用概率密度的头部信息指示随机变量的非高斯性，k1峭度定义为：
36、k1(b)＝e(cb(-|z|)). (8-2)
37、式中，z表示随机变量，尺度参数c为归一化常数，对于给定的参数b，尺度参数c的作用是使正态分布的k1峭度等于3；b的界限是固定的，2≤b≤20。k1峭度被定义为具有尺度参数c和形状参数b的标准变量z的连续函数的期望值，且函数是凸函数。显然，k1峭度对概率分布的峰值更加敏感，而传统的峭度k对概率分布的尾部更加敏感。也就是说，当概率密度函数的分布比高斯分布更加尖锐的时候，k1峭度比传统的峭度k更具有优势。根据上述分析，k1峭度是利用概率密度的峰值信息来指示随机变量的非高斯性，而实际地震反射系数序列为尖峰、长尾分布的。
38、因此，利用k1峭度的特性可以更优地度量反射系数的非高斯性。为了便于实际应用，本方法中将k1峭度近似表示成：
39、k1(b)∝e(b(-|z|))(8)式中，b≥2；符号∝表示正交于，也就是说k1峭度与e(b(-|z|))正交。
40、进一步地，步骤4)中，具体操作为：
41、根据最优的子波wopt(t)得到最优的逆滤波器：
42、
43、式中，wopt(ω)是最优子波wopt(t)频率表达形式，是最优子波wopt(ω)的复共轭；ξ为阻尼参数；最后，计算得到反褶积的结果。
44、更进一步地，阻尼参数ξ根据噪声水平选择，当噪声标准差为σ时，阻尼参数ξ近似设置为2.5σ；在实际地震数据处理中，由于噪声水平σ未知，利用最小绝对误差来估计实际数据的噪声误差。
45、本发明还公开了实现上述的基于k1峭度约束的地震频域反褶积方法的装置，包括：
46、地震数据获取单元，用于对地震数据进行处理，包括地震子波主频的获取和数据的预处理；
47、模型构建单元，用于对预处理后的地震数据进行模型构建，获得优化模型；
48、模型求解单元，用于对优化模型进行算法求解，得到输出最优aopt(ω)；
49、反分块处理单元，用于对最优aopt(ω)进行处理，获得重构的地震数据。
50、本发明还公开了一种计算机设备，包括存储器和处理器，所述存储器存储有计算机程序，所述处理器执行计算机程序时实现上述的基于k1峭度约束的地震频域反褶积方法的步骤。
51、本发明还公开了一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序被处理器执行时实现上述的基于k1峭度约束的地震频域反褶积方法的步骤。
52、与现有技术相比，本发明具有以下有益效果：
53、本发明公开的基于k1峭度约束的地震频域反褶积方法，通过引入广义beta母小波函数，能够更好的匹配不同类型的地震子波，同时，k1峭度能够更优地度量尖峰、长尾分布的地震反射系数序列的非高斯性。因此，基于k1峭度的频域反褶积方法，其可以有效度量地震反射系数序列的非高斯性，从而能够提高地震数据的分辨率，尤其是多个相邻的细小反射系数。因此利用本发明的方法可以有效地识别出多个含气薄储层，从而有利于后续的井位部署优化。
1.一种基于k1峭度约束的地震频域反褶积方法，其特征在于，包括以下步骤：
2.根据权利要求1所述的基于k1峭度约束的地震频域反褶积方法，其特征在于，广义beta小波的频率域定义为：
3.根据权利要求1所述的基于k1峭度约束的地震频域反褶积方法，其特征在于，步骤2)中，具体操作为：
4.根据权利要求3所述的基于k1峭度约束的地震频域反褶积方法，其特征在于，步骤3)中，具体操作为：
5.根据权利要求4所述的基于k1峭度约束的地震频域反褶积方法，其特征在于，步骤3)中，基于k1峭度来度量估计出的反射系数，其中：
6.根据权利要求4所述的基于k1峭度约束的地震频域反褶积方法，其特征在于，步骤4)中，具体操作为：
7.根据权利要求6所述的基于k1峭度约束的地震频域反褶积方法，其特征在于，阻尼参数ξ根据噪声水平选择，当噪声标准差为σ时，阻尼参数ξ近似设置为2.5σ；在实际地震数据处理中，由于噪声水平σ未知，利用最小绝对误差来估计实际数据的噪声误差。
8.实现权利要求1～7中任意一项所述的基于k1峭度约束的地震频域反褶积方法的系统，其特征在于，包括：
9.一种计算机设备，包括存储器和处理器，所述存储器存储有计算机程序，其特征在于，所述处理器执行计算机程序时实现权利要求1至7中任意一项所述的基于k1峭度约束的地震频域反褶积方法的步骤。
10.一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现权利要求1至7中任意一项所述的基于k1峭度约束的地震频域反褶积方法的步骤。
