一种地震子波和反射系数的反演方法及反演系统转让专利
申请号 : CN202210922138.1
文献号 : CN115327624B
文献日 : 2023-07-14
发明人 : 高静怀 , 徐文豪 , 田亚军
申请人 : 西安交通大学
摘要 :
本发明公开一种地震子波和反射系数的反演方法和反演系统。该方法假设地震子波具有紧支撑且光滑，同时假设反射系数是相对稀疏的，并构造了相应的反演地震子波和反射系数序列的优化问题。通过使用交替迭代，本发明将基于紧光滑性和相对稀疏度的地震子波和反射系数的联合反演问题拆分成子波反演子问题和反射系数反演子问题，并通过近端算法求解这两个子问题。相对于已有的地震子波和反射系数反演方法，本发明的优势包括最优参数易于选择、反演得到的反射系数横向连续性好等优点。
权利要求 :
1.一种地震子波和反射系数的反演方法，其特征在于，包括以下步骤：
S1、采集叠前地震记录数据，并对所述叠前地震记录数据进行预处理以形成叠后地震记录数据，其中，所述叠后地震记录数据包括多道叠后地震记录子数据；
S2、对所述多道叠后地震记录子数据取平均值以生成叠后地震记录平均道数据，对所述叠后地震记录平均道数据的振幅谱作光滑处理以形成零相位的初始子波，同时将初始反射系数序列取为零，在子波和反射系数的迭代开始时，将所述初始子波和所述初始反射系数序列分别设为最新的子波和最新的反射系数序列；
S3、使用所述最新的子波，通过反射系数近端算法求带有相对稀疏度约束的反射系数优化问题的最优解以形成最新的反射系数序列；
S4、使用所述最新的反射系数序列，通过子波近端算法求带有紧支撑约束和全变差光滑正则化项的子波优化问题的最优解以形成最新的子波；
S5、判断所述子波和所述反射系数的迭代是否满足终止条件；若是，则输出最终形成的子波和反射系数序列，否则返回步骤S3。
2.根据权利要求1所述的地震子波和反射系数的反演方法，其特征在于，所述对所述叠后地震记录平均道数据的振幅谱作光滑处理以形成零相位的初始子波具体包括以下步骤：对所述多道叠后地震记录子数据取平均值以生成叠后地震记录平均道数据，对所述叠后地震记录平均道数据的傅里叶振幅谱做光滑处理，以得到所述初始子波的振幅谱，然后再通过令所述初始子波的振幅谱为该初始子波的傅里叶变换进行反傅里叶变换以得到所述零相位的初始子波，其具体表达式如下：其中，Ntr表示叠后地震记录的道数，di表示第i道的叠后地震记录，n表示地震记录的时k间采样点数， 表示取上界函数，smooth (·)表示对振幅谱连续做k次七点三次多项式光滑。
3.根据权利要求1所述的地震子波和反射系数的反演方法，其特征在于，所述使用最新的子波，通过反射系数近端算法求带有相对稀疏度约束的反射系数优化问题的最优解以形成新的反射系数序列具体包括以下步骤：构造带有相对稀疏度约束的反射系数优化问题，其具体表达式如下：
其中，Tw表示子波w对应的托普利茨矩阵，ri表示第i道的反射系数序列， 和分别表示第i道非零反射系数中绝对幅度最小和最大的分量，cr表示相对稀疏度参数。
4.根据权利要求3所述的地震子波和反射系数的反演方法，其特征在于，所述cr的取值范围为0.001≤cr≤0.01。
5.根据权利要求3所述的地震子波和反射系数的反演方法，其特征在于，所述使用最新的反射系数序列，通过子波近端算法求带有紧支撑约束和全变差光滑正则化项的子波优化问题的最优解以形成所述最新的子波具体包括以下步骤：构造带有紧支撑约束和全变差光滑正则化项的子波优化问题，其具体表达式如下：
其中，λw表示通过地震道数规范化的全变差光滑正则化项的系数，Δt表示时间采样间隔，f0是将地震数据各道的峰值频率作平均后得到的地震数据主频估计值，||w||0表示w的零范数， 表示通过一个主周期内的采样点数规范化后的子波非零分量个数参数，子波的长度为2m‑1。
6.根据权利要求5所述的地震子波和反射系数的反演方法，其特征在于，所述λw的取值范围为0.001≤λw≤0.01， 取1‑5之间的正整数。
7.根据权利要求1‑6中任一项所述的地震子波和反射系数的反演方法，其特征在于，对所述叠前地震记录数据进行预处理的方法包括保幅去噪、真振幅恢复、静校正处理中的至少一种。
8.根据权利要求1‑6中任一项所述的地震子波和反射系数的反演方法，其特征在于，所述叠后地震记录数据包括叠后地震记录的每一道数据。
9.一种地震子波和反射系数的反演系统，其特征在于，所述反演系统通过权利要求1‑8中任一项所述的地震子波和反射系数的反演方法实现。
说明书 :
一种地震子波和反射系数的反演方法及反演系统
技术领域
[0001] 本发明属于地球物理勘探技术领域，尤其涉及一种地震子波和反射系数的反演方法及反演系统。
背景技术
[0002] 子波和反射系数反演在地震资料处理中有着重要的应用，包括高分辨率处理和地震波阻抗反演等方面。子波和反射系数反演的主要困难在于地震记录的非平稳性以及子波和反射系数反演的病态性。
[0003] 目前克服地震记录非平稳性的常用方法包括有重叠的分段、品质因子补偿和深度学习补偿等。关于克服子波反演的病态性，目前已有的方法包括相位校正、压缩映射算子、L1范数稀疏正则化、全变差正则化、深度学习反演等。关于克服反射系数反演的病态性，目前已有的方法包括固定非零反射系数的个数后反演反射系数的位置和振幅、L1范数稀疏正则化、基于多道反射系数构造齐次线性方程组、基于包络峰值个数的L0范数约束、深度学习反演等。
[0004] 目前的子波和反射系数反演方法具有以下不足：
[0005] (1)基于深度学习的反演方法一般需要有合适的标签数据才能发挥比较好的效果，但实际中测井数据相对较少，且测井数据和地震数据的匹配存在困难；
[0006] (2)正则化项和约束项往往需要选择合适的参数才能有效发挥作用，基于人工调试的参数选取方式不利于大规模的工业应用，而基于L曲线和广义交叉验证的选取方式计算成本较高；
[0007] (3)对反射系数反演，施加横向约束会导致强间断的区域失真，而在不施加横向约束的情况下，基于L1范数正则化或L0范数约束得到的反射系数横向连续性有待提高。
发明内容
[0008] 本发明旨在至少在一定程度上解决上述技术问题之一或至少提供一种有用的商业选择。为此，本发明的一个目的在于提出一种地震子波和反射系数的反演方法，本发明的地震子波和反射系数的反演方法假设地震子波具有紧支撑且光滑，同时假设反射系数是相对稀疏的；通过使用交替迭代，本发明将基于紧光滑性和相对稀疏度的子波和反射系数联合反演问题拆分成子波反演子问题和反射系数反演子问题，并使用近端算法求解这两个子问题，本发明的优势包括最优参数易于选择、反演得到的反射系数分辨率高且横向连续性好等。本发明的另一个目的在于提供一种地震子波和反射系数的反演系统。
[0009] 根据本发明的地震子波和反射系数的反演方法包括以下步骤：
[0010] S1、采集叠前地震记录数据，并对所述叠前地震记录数据进行预处理以形成叠后地震记录数据，其中，所述叠后地震记录数据包括多道叠后地震记录子数据；
[0011] S2、对所述多道叠后地震记录子数据取平均值以生成叠后地震记录平均道数据，对所述叠后地震记录平均道数据的振幅谱作光滑处理以形成零相位的初始子波，同时将初始反射系数序列取为零，在子波和反射系数的迭代开始时，将所述初始子波和所述初始反射系数序列分别设为最新的子波和最新的反射系数序列；
[0012] S3、使用所述最新的子波，通过反射系数近端算法求带有相对稀疏度约束的反射系数优化问题的最优解以形成最新的反射系数序列；
[0013] S4、使用所述最新的反射系数序列，通过子波近端算法求带有紧支撑约束和全变差光滑正则化项的子波优化问题的最优解以形成最新的子波；
[0014] S5、判断所述子波和所述反射系数迭代是否满足终止条件；若是，则输出最终形成的子波和反射系数序列，否则返回步骤S3。
[0015] 本发明的地震子波和反射系数的反演方法假设地震子波具有紧支撑且光滑，同时假设反射系数是相对稀疏的；通过使用交替迭代，本发明将基于紧光滑性和相对稀疏度的子波和反射系数联合反演问题拆分成子波反演子问题和反射系数反演子问题，并使用近端算法求解这两个子问题，本发明的优势包括最优参数易于选择、反演得到的反射系数分辨率高且横向连续性好等。
[0016] 另外，根据本发明的地震子波和反射系数的反演方法还可以具有以下的技术特征：
[0017] 所述对所述叠后地震记录平均道数据的振幅谱作光滑处理以形成零相位的初始子波具体包括以下步骤：
[0018] 对所述多道叠后地震记录子数据取平均值以生成叠后地震记录平均道数据，对所述叠后地震记录平均子数据的傅里叶振幅谱做光滑处理，以得到所述初始子波的振幅谱，然后再通过令所述初始子波的振幅谱为该初始子波的傅里叶变换进行反傅里叶变换以得到所述零相位的初始子波，其具体表达式如下：
[0019]
[0020] 其中，Ntr表示叠后地震记录的道数，di表示第i道的叠后地震记录，n表示地震记录k的时间采样点数， 表示取上界函数，smooth (·)表示对振幅谱连续做k次七点三次多项式光滑。
[0021] 所述使用最新的子波，通过反射系数近端算法求带有相对稀疏度约束的反射系数优化问题的最优解以形成新的反射系数序列具体包括以下步骤：
[0022] 构造带有相对稀疏度约束的反射系数优化问题，其具体表达式如下：
[0023]
[0024] 其中，Tw表示子波w对应的托普利茨矩阵，ri表示第i道新的反射系数序列，和 分别表示第i道非零反射系数中绝对幅度最小和最大的分量，cr表示相对稀疏度参数。
[0025] 所述cr的取值范围为0.001≤cr≤0.01。
[0026] 所述使用最新的反射系数序列时，通过子波近端算法求带有紧支撑约束和全变差光滑正则化项的子波优化问题的最优解以形成所述最新的子波具体包括以下步骤：
[0027] 构造带有紧支撑约束和全变差光滑正则化项的子波优化问题，其具体表达式如下：
[0028]
[0029] 其中，λw表示通过地震道数规范化的全变差光滑正则化项的系数，Δt表示时间采样间隔，f0是将地震数据各道的峰值频率作平均后得到的地震数据主频估计值，||w||0表示w的零范数， 表示通过一个主周期内的采样点数规范化后的子波非零分量个数参数，子波的长度为2m‑1。
[0030] 所述λw的取值范围为0.001≤λw≤0.01， 取1‑5之间的正整数。
[0031] 对所述叠前地震记录数据进行预处理的方法包括保幅去噪、真振幅恢复、静校正处理中的至少一种。
[0032] 所述叠后地震记录数据包括叠后地震记录的每一道数据。
[0033] 本发明还提供了一种地震子波和反射系数的反演系统，所述反演系统通过上述任意的地震子波和反射系数的反演方法实现。
[0034] 本发明具有以下有益效果：
[0035] 本发明在子波和反射系数反演中，对子波引入紧光滑约束，同时对反射系数引入相对稀疏约束，使得本发明的反演方法参数易于选取且所得的反射系数横向连续性好。与目前常用的TSMF子波和反射系数反演方法对比，本发明所提的方法可以在更高效的参数选取方法下得到反射系数横向连续性更好、更符合地质构造特点的反演结果。进一步，更高效的参数选取可以为子波和反射系数反演的大规模工业化应用提供技术支撑，而反射系数横向连续性更好、更符合地质构造特点的反演结果可以为更精确的油气勘探提供技术支撑。
[0036] 本发明的附加方面和优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本发明的实践了解到。
附图说明
[0037] 本发明的上述和/或附加的方面和优点从结合下面附图对实施例的描述中将变得明显和容易理解，其中：
[0038] 图1是本发明一个实施例的流程示意图；
[0039] 图2是本发明的一个实施例的某海上实际地震记录的二维剖面图；
[0040] 图3是本发明的一个实施例的对比方法得到的反射系数序列图；
[0041] 图4是本发明的一个实施例的采用本发明方法得到的反射系数序列图；
[0042] 图5是本发明的一个实施例的原始地震记录、对比方法得到的反射系数、本发明得到的反射系数的频谱对比图；
[0043] 图6是本发明的一个实施例的对比方法和本发明得到的子波对比图。
[0044] 具体实施方法
[0045] 下面详细描述本发明的实施例，所述实施例的示例在附图中示出，其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施例是示例性的，旨在用于解释本发明，而不能理解为对本发明的限制。
[0046] 地震子波和反射系数反演通过对叠后地震数据进行处理，得到地震子波和反射系数序列，是地震资料高分辨处理和地震波阻抗反演的重要手段。本发明假设地震子波具有紧支撑且光滑，同时假设反射系数是相对稀疏的，从而构造相应的子波和反射系数反演优化问题，并通过交替迭代方法和近端算法求解相应优化问题，最终得到了参数容易选择且反射系数横向连续性好的地震子波和反射系数反演方法。
[0047] 图1是本发明的一个实施例的流程示意图；图2是本发明的一个实施例的某海上实际地震记录二维剖面图；图3是本发明的一个实施例的对比方法得到的反射系数序列图；图4是本发明的一个实施例的采用本发明方法得到的反射系数序列图；图5是本发明的一个实施例的原始地震记录、对比方法得到的反射系数、本发明得到的反射系数的频谱对比图；图
6是本发明的一个实施例的对比方法和本发明得到的子波对比图。参考图1‑图6，本发明提供了一种地震子波和反射系数的反演方法，所述反演方法包括以下步骤：
6是本发明的一个实施例的对比方法和本发明得到的子波对比图。参考图1‑图6，本发明提供了一种地震子波和反射系数的反演方法，所述反演方法包括以下步骤：
[0048] S1、采集叠前地震记录数据，并对所述叠前地震记录数据进行预处理以形成叠后地震记录数据，其中，所述叠后地震记录数据包括多道叠后地震记录子数据。
[0049] 具体的，可通过保幅去噪、真振幅恢复、静校正处理等方法来对采集到的地震数据做预处理，以得到高信噪比的叠后地震记录数据。在具体实施中，既可以采用保幅去噪、真振幅恢复、静校正处理等方法中的一种来对采集到的地震数据做预处理，也可以采用保幅去噪、真振幅恢复、静校正处理等方法中的两种或者三种进行叠加处理，从而得到更加精确的叠后地震记录数据。本具体实施中，叠后地震记录数据记为{di},i＝1,2,…,Ntr，其中Ntr表示叠后地震记录的道数。在一个实施例中，叠后地震记录数据可包括叠后地震记录的每一道数据，即叠后地震记录的每一道数据即形成了叠后地震记录的数据集。在其他实施例中，叠后地震记录数据也可以采取叠后地震记录的每一道数据中的一部分。
[0050] S2、对所述多道叠后地震记录子数据取平均值以生成叠后地震记录平均道数据，对所述叠后地震记录平均道数据的振幅谱作光滑处理以形成零相位的初始子波，同时将初始反射系数序列取为零，在子波和反射系数的迭代开始时，将初始子波和初始反射系数序列设为最新的子波和最新的反射系数序列。
[0051] 具体的，对所述叠后地震记录平均道数据的傅里叶振幅谱做光滑处理，以得到所述初始子波的振幅谱，然后再通过令所述初始子波的振幅谱为该初始子波的傅里叶变换进行反傅里叶变换以得到所述零相位的初始子波，其具体表达式如下：
[0052]
[0053] 其中，n表示地震记录的时间采样点数， 表示取上界函数，smoothk(·)表示对振幅谱连续做k次七点三次多项式光滑，初始化反射系数序列 初始化迭代次数k＝0。
[0054] 在具体实施中，在子波和反射系数的迭代开始时，将初始子波和初始反射系数分别设为最新的子波、最新的反射系数序列，从而便于后序的迭代操作。
[0055] S3、使用所述最新的子波，通过反射系数近端算法求带有相对稀疏度约束的反射系数优化问题的最优解以形成最新的反射系数序列。
[0056] 具体的，假设最新的子波为w(k)，构造如下带有相对稀疏度约束的反射系数优化问题：
[0057]
[0058] 其中，Tw表示子波w对应的Toeplitz(托普利茨)矩阵， 和 分别表示第i道非零反射系数中绝对幅度最小和最大的分量，cr表示相对稀疏度参数。对进行保幅去噪处理了的高信噪比地震数据，cr一般可取[0.001,0.01]，在本实施例中，cr可取0.005。由于在子波固定时，各道的反射系数反演是独立的，因此优化问题(2)可简化为如下形式：
[0059]
[0060] 其中，d和r分别表示任一道的地震记录和相应的反射系数序列。对优化问题(4)中的 在搜索点s处做一阶泰勒展开，其中搜索点s可基于迭代变量通过Nesterov方法确定，并加入近端正则化项 其中，L是迭代步长的倒数，从而可将优化问题(4)转化为如下形式：
[0061]
[0062] 根据反射系数近端算法，优化问题(5)可转化为如下的向量逼近问题：
[0063]
[0064] 其中， 参考已有的关于固定稀疏度问题的解法，本发明给出了如下关于向量逼近问题的解：
[0065]
[0066] 其中， 表示v中绝对幅值大于 的分量个数，[v]K+1表示v中绝对幅值第K+1大的分量。
[0067] S4、使用所述最新的反射系数序列，通过子波近端算法求带有紧支撑约束和全变差光滑正则化项的子波优化问题的最优解以形成最新的子波。
[0068] 具体的，假设最新的反射系数为 子波的长度为2m‑1，可以理解的是，子波的长度大于0,构造如下带有紧支撑约束和全变差(Total Variation,TV)光滑正则化项的子波优化问题：
[0069]
[0070] 其中，λw表示通过地震道数规范化的TV光滑正则化项的系数，Δt表示时间采样间隔，f0是将地震数据各道的峰值频率作平均得到的估计的地震数据主频，||w||0表示w的零范数也即子波中非零分量的个数，Nf0表示通过一个主周期(对应f0)内的采样点数规范化后的子波非零分量个数参数。特别地，对进行保幅去噪处理的高信噪比地震数据，λw一般可取[0.001,0.01]， 一般可取1至5的正整数。在本实施例中，λw取0.005, 取3。将m固定为同时利用卷积的对称性，可将优化问题(3)简化为如下形式：
[0071]
[0072] 其中， 是如下n×(2n‑1)矩阵 的中间2m‑1列：
[0073]
[0074] 类似步骤S3的做法，在搜索点s处对lossw(w)做一阶泰勒展开并加入近端正则化项 可将优化问题(8)转化为如下形式
[0075]
[0076] 根据子波近端算法，优化问题(10)可转化为如下的向量逼近问题：
[0077]
[0078] 其中， 类似于Fused Lasso问题的求解，我们将优化问题(11)转化为如下的嵌套优化问题：
[0079]
[0080] 其中
[0081]
[0082] 优化问题(13)是典型的TV正则化问题，本发明采用经典的迭代重加权最小二乘(Iterative Reweighted Least Square,IRLS)方法求解。参考固定稀疏度约束优化问题的解法，优化问题(12)可采用如下形式的解：
[0083]
[0084] 其中，
[0085] S5、判断子波和反射系数迭代是否满足终止条件；若是，则输出最终反演得到的子波和反射系数序列，否则返回步骤S3。
[0086] 具体的，更新迭代次数k＝k+1，判断迭代是否满足终止条件(通常选为最新的子波与上次迭代的子波的二范数误差小于一定阈值)，若满足，则输出最终反演得到的子波和反射系数序列；否则，返回步骤S3。总体而言，本发明的地震子波和反射系数的反演方法，其通(k)过步骤S3和步骤S4不断持续循环迭代，在一个实施例中，使用最新的子波w 通过反射系数近端算法求带有相对稀疏度约束的反射系数优化问题的最优解形成新的反射系数
此时该反射系数 即为形成的最新的反射系数序列；然后使用该最新的反射系数序列为 并通过子波近端算法求带有紧支撑约束和全变差光滑正则化项的子波优化问题(k+1) (k+1)
的最优解以形成最新的子波w ；紧接着，该形成的最新的子波w 又通过反射系数近端算法的求解形成最新的反射系数 然后又使用该最新形成的反射系数序列 通
(k+2)
过子波近端算法的求解形成最新的子波w ……；如此循环往复的迭代，直到最新的子波与上次迭代的子波的二范数误差小于一定阈值后停止迭代。
此时该反射系数 即为形成的最新的反射系数序列；然后使用该最新的反射系数序列为 并通过子波近端算法求带有紧支撑约束和全变差光滑正则化项的子波优化问题(k+1) (k+1)
的最优解以形成最新的子波w ；紧接着，该形成的最新的子波w 又通过反射系数近端算法的求解形成最新的反射系数 然后又使用该最新形成的反射系数序列 通
(k+2)
过子波近端算法的求解形成最新的子波w ……；如此循环往复的迭代，直到最新的子波与上次迭代的子波的二范数误差小于一定阈值后停止迭代。
[0087] 本发明还提供了一种地震子波和反射系数的反演系统，所述反演系统通过上述的各反演方法实现。
[0088] 实际地震资料测试
[0089] 这一节使用实际地震资料，并通过与目前常用的Toeplitz稀疏矩阵分解(Toeplitz‑Sparse Matrix Factorization,TSMF)地震子波和反射系数反演方法的对比，来验证本发明提出的紧光滑且相对稀疏(Compact Smoothness and Relative Sparsity,CSRS)方法的有效性。这里，本发明所提CSRS方法中的参数均使用上述的固定默认参数。图2给出了某海上实际叠后地震数据的一个二维剖面，该剖面共801道，截取的时间范围为1498～2096毫秒，时间采样间隔为2毫秒。测井数据位于第275道的位置。图3、4分别给出了通过TSMF和CSRS方法得到的反射系数剖面。图3和图4的对比表明，在使用固定默认参数的情况下，CSRS可以得到分辨率高(即地质构造更清晰)且比TSMF横向连续性更好(即更符合地质构造特点)的反射系数剖面。图5给出了原始地震数据、TSMF所得反射系数剖面、CSRS所得反射系数剖面的振幅谱比较。图5表明，在有效拓宽高频的同时，CSRS方法相对TSMF方法可以更好地保留地震数据的低中频部分(即所得反射系数的保真性更好)。图6给出了相应的子波比较，图6表明CSRS方法可以得到类似雷克子波形状的合理地震子波。
[0090] 需要指出的是，横向连续性更好是指反射系数的空间横向连续性更好。由于实际数据没有真实地震子波作为对比，一般实际数据的子波反演往往是要求能得到类似雷克子波形状的合理地震子波。另外，虽然TSMF方法得到的子波比CSRS更加光滑，但由于实际震源子波的复杂性以及地层对子波的衰减效应，不一定越光滑的子波就越合理。实际上，由于CSRS得到的子波对应的反射系数横向连续性更好(更符合地质构造特点)，通过本发明的方法得到的子波可能更接近于实际的地震子波。
[0091] 本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。
[0092] 尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制，本领域的普通技术人员在不脱离本发明的原理和宗旨的情况下在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。