CN117741750A - 一种基于Radon变换的多道叠前反褶积方法及系统 - Google Patents
一种基于Radon变换的多道叠前反褶积方法及系统 Download PDFInfo
- Publication number
- CN117741750A CN117741750A CN202410190044.9A CN202410190044A CN117741750A CN 117741750 A CN117741750 A CN 117741750A CN 202410190044 A CN202410190044 A CN 202410190044A CN 117741750 A CN117741750 A CN 117741750A
- Authority
- CN
- China
- Prior art keywords
- deconvolution
- stack
- radon
- prestack
- seismic data
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 229910052704 radon Inorganic materials 0.000 title claims abstract description 101
- SYUHGPGVQRZVTB-UHFFFAOYSA-N radon atom Chemical compound [Rn] SYUHGPGVQRZVTB-UHFFFAOYSA-N 0.000 title claims abstract description 101
- 230000009466 transformation Effects 0.000 title claims abstract description 52
- 238000000034 method Methods 0.000 title claims abstract description 50
- 230000002159 abnormal effect Effects 0.000 claims abstract description 42
- 238000012545 processing Methods 0.000 claims abstract description 31
- 238000012549 training Methods 0.000 claims abstract description 21
- 230000001629 suppression Effects 0.000 claims description 11
- 239000011159 matrix material Substances 0.000 claims description 7
- 238000004364 calculation method Methods 0.000 claims description 5
- 230000008602 contraction Effects 0.000 claims description 2
- 238000012360 testing method Methods 0.000 description 27
- 238000001228 spectrum Methods 0.000 description 7
- 230000008569 process Effects 0.000 description 6
- 230000001737 promoting effect Effects 0.000 description 6
- 230000006872 improvement Effects 0.000 description 4
- 238000005070 sampling Methods 0.000 description 4
- 230000002776 aggregation Effects 0.000 description 3
- 238000004220 aggregation Methods 0.000 description 3
- 230000000694 effects Effects 0.000 description 3
- 238000009499 grossing Methods 0.000 description 3
- 230000036039 immunity Effects 0.000 description 3
- 238000012986 modification Methods 0.000 description 3
- 230000004048 modification Effects 0.000 description 3
- 230000003595 spectral effect Effects 0.000 description 3
- 238000010521 absorption reaction Methods 0.000 description 2
- 230000008901 benefit Effects 0.000 description 2
- 230000015572 biosynthetic process Effects 0.000 description 2
- 230000007423 decrease Effects 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 238000013112 stability test Methods 0.000 description 2
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 1
- 241001148606 Bennia Species 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000012512 characterization method Methods 0.000 description 1
- 230000006835 compression Effects 0.000 description 1
- 238000007906 compression Methods 0.000 description 1
- 238000013170 computed tomography imaging Methods 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000002708 enhancing effect Effects 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 230000010354 integration Effects 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 230000000087 stabilizing effect Effects 0.000 description 1
- 238000010998 test method Methods 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
- 230000002087 whitening effect Effects 0.000 description 1
Classifications
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T10/00—Road transport of goods or passengers
- Y02T10/10—Internal combustion engine [ICE] based vehicles
- Y02T10/40—Engine management systems
Landscapes
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明涉及地震数据处理技术领域，特别是涉及一种基于Radon变换的多道叠前反褶积方法及系统，方法包括：获取待处理叠前地震数据集；将所述待处理叠前地震数据集输入预设的反褶积模型中，输出叠前反褶积后的地震数据集，其中，所述反褶积模型基于训练集训练获得，所述训练集包括满足拉普拉斯分布的异常振幅噪声和满足高斯分布的随机噪声的CMP道集数据；所述反褶积模型通过双曲Radon变换约束结合叠前反褶积算法构建。本发明能够在有效提高分辨率的同时，提高地震道集的信噪比，适用于低信噪比的地震道集处理。
Description
技术领域
本发明涉及地震数据处理技术领域，特别是涉及一种基于Radon变换的多道叠前反褶积方法及系统。
背景技术
地震资料处理中，如何提高地震数据分辨率是一个亟待解决的难点问题。反褶积技术是提高地震数据分辨率的重要算法之一。根据所应用的数据，反褶积可以划分为叠前反褶积（Levin, 1989）以及叠后反褶积（Chen et al., 2022）两类。由于叠后地震数据具有较高的信噪比，因此叠后反褶积通常较为简单。同时，作用在时间域的反褶积算法，包括基于L2范数的反褶积方法（Golub et al., 1999）、考虑地层吸收衰减的非稳态反褶积（Bickel and Natarajan, 1985）、基于L1范数约束的稀疏脉冲反褶积（Taylor, 1979），以及考虑横向连续性的多道反褶积（Ma et al., 2018）等方法。作用在频率域的反褶积算法则包括谱白化（Bennia and Riad., 1990）、谱反演（Chen et al., 2021）和谱模拟（Jianget al., 2020）等算法，这些方法根据设定的振幅谱形态，使地震数据的振幅谱接近预设振幅，从而实现高分辨率的频率域处理。当考虑非稳态数据的特征时，可以使用具备时频定位能力的时频分析方法来实现高分辨率处理。
与叠后数据不同，叠前数据具有明显的时距规律性。因此，在进行反褶积处理时，需要兼顾该特征，避免逐道反褶积导致的横向不连续性，从而影响后续的反演和储层预测工作。对于低信噪比的数据，如何平衡信噪比和高分辨率是一个关键问题。一种增强数据横向连续性的反褶积方法是考虑空间方向或沿同相轴平滑特征的方法，这需要使用局部倾角来（Liu et al., 2010）约束数据。然而，对于低信噪比的数据，很难准确获取到同相轴的倾角，这限制了这类方法的研究进展。另一种思路是假设高分辨率数据在变换域中具有稀疏性，可以使用FK变换（Bai and Wu, 2018）、小波变换（Daubechies, 1996）、Curvelet变换（Starck, 2002）和非局部相似性（Bonar and Sacchi, 2012）等稀疏促进变换来压制随机噪声，并改善高分辨率处理的不适定性问题。
现有反褶积技术在处理叠前数据时存在几个缺陷。首先是地震数据的低信噪比问题，由于面波和高频噪声的存在，反褶积算法容易放大噪声，干扰处理结果。其次，现有算法大多是逐道处理，导致反褶积后结果存在一定的横向不连续性，影响后续的反演和储层预测工作。此外，针对空间方向连续性的问题，现有方法大多考虑平滑特征，即以L2范数作为正则化条件，但在低信噪比数据中难以准确获取同相轴倾角，限制了研究深度。最后，现有反褶积算法在高分辨率处理上存在不适定性，目前采用的稀疏促进变换技术难以得到较好的多道联合处理结果。这些缺陷限制了现有反褶积技术在处理低信噪比数据、保持横向连续性和实现高分辨率处理方面的能力。因此，亟需一种基于Radon变换的多道叠前反褶积方法及系统。
发明内容
本发明的目的是提供一种基于Radon变换的多道叠前反褶积方法及系统，在有效提高分辨率的同时，提高地震道集的信噪比，适用于低信噪比的地震道集处理。
为实现上述目的，本发明提供了如下方案：
一种基于Radon变换的多道叠前反褶积方法，包括：
获取待处理叠前地震数据集；
将所述待处理叠前地震数据集输入预设的反褶积模型中，输出叠前反褶积后的地震数据集，其中，所述反褶积模型基于训练集训练获得，所述训练集包括满足拉普拉斯分布的异常振幅噪声和满足高斯分布的随机噪声的CMP道集数据；所述反褶积模型通过双曲Radon变换约束结合叠前反褶积算法构建。
可选地，所述反褶积模型包括：异常振幅噪声判断单元、第一反褶积单元和第二反褶积单元，所述异常振幅噪声判断单元用于判断所述待处理叠前地震数据集中是否存在异常振幅噪声；所述第一反褶积单元用于当所述待处理叠前地震数据集中存在异常振幅噪声时，基于第一Radon域叠前反褶积算法对所述待处理叠前地震数据集进行叠前反褶积处理；所述第二反褶积单元用于当所述待处理叠前地震数据集中不存在异常振幅噪声时，基于第二Radon域叠前反褶积算法对所述待处理叠前地震数据集进行叠前反褶积处理。
可选地，所述第一Radon域叠前反褶积算法以L1范数作为拟合项，所述第二Radon域叠前反褶积算法以L2范数作为拟合项。
可选地，所述第一Radon域叠前反褶积算法和第二Radon域叠前反褶积算法为：
其中，R为多道反射系数，W为地震子波构成的卷积矩阵，S为多道地震数据，λ为正则化算子，B为双曲Radon变换，|| ||1为数据绝对值的和，|| ||2为数据平方根的和。
可选地，所述第一Radon域叠前反褶积算法和第二Radon域叠前反褶积算法的求解方法为：
设置Radon域数据高斯阈值，以及迭代次数；
基于所述高斯阈值，结合迭代阈值收缩算法对所述待处理叠前地震数据集在Radon域进行噪声压制，直至达到所述迭代次数，获取压制噪声后的Radon域叠前地震数据；
对压制噪声后的Radon域叠前地震数据进行Radon反变换，获取变换域压制噪声后的叠前地震数据集，对所述变换域压制噪声后的叠前地震数据集进行叠前反褶积处理，获取叠前反褶积后的地震数据集。
可选地，所述Radon域数据阈值的计算方法为：
其中，H为阈值，为振幅，h gauss 为当前曲率下构建的高斯窗函数，/>为曲率变量，/>为一次波聚集的曲率位置，/>为窗的大小，/>为当前计算时刻。
为进一步实现上述目的，本发明还提供了一种基于Radon变换的多道叠前反褶积系统，包括：地震数据获取模块、反褶积处理模块；
所述地震数据获取模块，用于获取待处理叠前地震数据集；
所述反褶积处理模块，用于将所述待处理叠前地震数据集输入预设的反褶积模型中，输出叠前反褶积后的地震数据集，其中，所述反褶积模型基于训练集训练获得，所述训练集包括包含异常振幅噪声和随机噪声的CMP道集数据；所述反褶积模型通过双曲Radon变换约束结合叠前反褶积算法构建。
可选地，所述反褶积模型包括：异常振幅噪声判断单元、第一反褶积单元和第二反褶积单元，所述异常振幅噪声判断单元用于判断所述待处理叠前地震数据集中是否存在异常振幅噪声；所述第一反褶积单元用于当所述待处理叠前地震数据集中存在异常振幅噪声时，基于第一Radon域叠前反褶积算法对所述待处理叠前地震数据集进行叠前反褶积处理；所述第二反褶积单元用于当所述待处理叠前地震数据集中不存在异常振幅噪声时，基于第二Radon域叠前反褶积算法对所述待处理叠前地震数据集进行叠前反褶积处理。
本发明的有益效果为：
本发明针对于叠前地震数据，采用Radon变换作为稀疏促进变换，实现叠前数据的反褶积算法，考虑了地震道集的空间方向的可预测性，得到的反褶积结果具备横向连续性，适用于低信噪比的地震道集处理；根据实际数据的噪声类型不同，可灵活选择针对高斯噪声的L2范数作为拟合项、针对拉普拉斯噪声的L1范数作为拟合项，为处理人员提供了参考与指导。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1为本发明实施例的一种基于Radon变换的多道叠前反褶积方法流程图；
图2为本发明实施例的测试一的测试数据，其中，（a）为合成的无噪CMP道集，（b）为加入Laplace噪声的CMP道集，（c）为Radon变换后的域数据，（d）为根据/>域数据构建的时变高斯阈值；
图3为本发明实施例的测试一不同反褶积后的道集，其中，（a）为吉洪诺夫反褶积后的道集，（b）为F-K变换作为稀疏促进变换的反褶积后的道集，（c）为L2范数作为拟合项的Radon变换作为稀疏促进变换的反褶积后的道集，（d）为L1范数作为拟合项的Radon变换作为稀疏促进变换的反褶积后的道集；
图4为本发明实施例的测试一的不同反褶积后的反射系数及地震记录单道对比图；
图5为本发明实施例的测试一的不同反褶积前后振幅谱对比图；
图6为本发明实施例的测试一的不同反褶积稳定性测试结果图，其中，（a）为不同反褶积的信噪估计比，（b）为不同反褶积的结构相似度；
图7为本发明实施例的测试一的不同反褶积对拉普拉斯噪声的抗噪性测试结果图，其中，（a）为不同反褶积的信噪估计比，（b）为不同反褶积的结构相似度；
图8为本发明实施例的测试一的不同反褶积对高斯噪声的抗噪性测试结果图，其中，（a）为不同反褶积的信噪估计比，（b）为不同反褶积的结构相似度；
图9为本发明实施例的测试二的测试数据，其中，（a）为原始地震数据，（b）为原始地震数据第一道的时频谱，（c）为Radon变换后的域数据，（d）为根据/>域数据构建的时变高斯阈值；
图10为本发明实施例的测试二不同反褶积后的数据，其中，（a）为吉洪诺夫反褶积后的数据，（b）为F-K变换作为稀疏促进变换的反褶积后的数据，（c）为L2范数作为拟合项的Radon变换作为稀疏促进变换的反褶积后的数据，（d）为L1范数作为拟合项的Radon变换作为稀疏促进变换的反褶积后的数据；
图11为本发明实施例的测试二的不同反褶积前后振幅谱对比图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。
本实施例提供了一种基于Radon变换的多道叠前反褶积方法，如图1所示，包括：
获取待处理叠前地震数据集；
将待处理叠前地震数据集输入预设的反褶积模型中，输出叠前反褶积后的地震数据集，其中，反褶积模型基于训练集训练获得，训练集包括满足拉普拉斯分布的异常振幅噪声和满足高斯分布的随机噪声的CMP道集数据；反褶积模型通过双曲Radon变换约束结合叠前反褶积算法构建。
反褶积模型包括：异常振幅噪声判断单元、第一反褶积单元和第二反褶积单元，异常振幅噪声判断单元用于判断待处理叠前地震数据集中是否存在异常振幅噪声；第一反褶积单元用于当待处理叠前地震数据集中存在异常振幅噪声时，基于第一Radon域叠前反褶积算法对待处理叠前地震数据集进行叠前反褶积处理；第二反褶积单元用于当待处理叠前地震数据集中不存在异常振幅噪声时，基于第二Radon域叠前反褶积算法对待处理叠前地震数据集进行叠前反褶积处理。
其中，第一Radon域叠前反褶积算法以L1范数作为拟合项，第二Radon域叠前反褶积算法以L2范数作为拟合项。
具体过程如下：
获取原始CMP道集S及子波矩阵W，判断原始CMP道集S中是否存在异常振幅噪声，若不存在异常振幅噪声，则选择对应高斯噪声的L2范数作为拟合项的Radon域叠前反褶积算法，对原始CMP道集S进行叠前反褶积，获取叠前反褶积后的CMP道集Z；若存在异常振幅噪声，则选择对应高斯噪声的L1范数作为拟合项的Radon域叠前反褶积算法，对原始CMP道集S进行叠前反褶积，获取叠前反褶积后的CMP道集Z；
其中，Radon域叠前反褶积算法为：将双曲Radon变换作为稀疏促进变换引入到叠前反褶积算法中，并在Radon域压制叠前数据噪声，得到压制噪声后的数据，并将其进行Radon反变换得到压制噪声后的CMP数据/>，将/>代入到叠前多道反褶积算法中，最终得到反褶积后的CMP道集/>。
具体计算过程如下：
根据褶积原理，叠前地震数据可近似写为地震子波/>同反射系数/>的褶积：
（1）
其中，为时间变量，/>表示随机噪声。将其改写为多道的矩阵形式：
（2）
其中，是由地震子波构成的卷积矩阵；/>，表示多道反射系数，/>为多道地震数据。无论是高分辨率处理，还是反射系数反演，从/>中得到/>都是不适定的过程，因此，考虑增加稳定项，用以提高/>求解的稳定性；依据Tikhonov正则化，易得：
（3）
根据最小二乘原理，最小化目标方程（3），有：
（4）
其中，表示正则化算子，/>为单位矩阵，式（4）的/>为反褶积后的数据，但式（3）存在假设：
a、要求噪声满足高斯分布；
b、要求高分辨率后数据满足高斯分布。
对于信噪比较低的叠前数据，如果数据中含有振幅异常能量，或其它强振幅噪声，直接最小化式（3）会导致不稳定。当假设噪声满足拉普拉斯分布时，修改式（3）为：
（5）
其中，对/>起平滑作用，类似于低通滤波，而实际上，大部分噪声都在地震频带之内，而多道数据之间的噪声具有随机性。因此，有必要考虑多道数据的空间特征，从而一定程度上压制地震频带内的噪声，增强高分辨率数据的稳定性。定义稀疏促进变换/>，认为高分辨率数据经过变换/>后，有效信号稀疏，算法不稳定产生的噪声和数据本身的噪声非稀疏，即：
（6）
式（6）中第一式为常见的Least absolute shrinkage and selection operator，可以应用梯度下降和软阈值交替迭代求解；而式（6）中第二式需用次梯度下降和软阈值交替迭代进行求解，即：
（7）
其中，表示阶跃函数，/>表示叉乘。叠前地震道集大多数符合双曲线形态，而多次波、随机噪声大多不满足该形态。因此，本实施例中选用双曲Radon变换作为稀疏促进变换，下面介绍Radon变换的基本原理。
Radon变换是沿预设轨迹进行线积分，将时空域数据投影到域，使满足该轨迹的信号聚集为速度相关的能量团，不满足该轨迹的信号或者随机噪声被映射到其它位置，基于该思路，Radon变换常被用于图像处理中的线性特征增强、医学CT成像和多次波压制中。预设轨迹为双曲线的2（d）Radon变换为：
（8）
（9）
其中，表示沿偏移距方向的Radon变换算子，/>和/>为Radon域数据。/>表示频率，/>表示虚数单位，/>为最大道数。当已知Radon域系数/>，左乘/>，对频率方向做傅里叶逆变换，可得到时空域地震数据。当已知频率空间域地震数据，可基于式（10）求解Radon域系数最优解/>，即：
（10）
式中，表示离散傅里叶变换矩阵，由于式（10）中/>非满秩，导致/>有无穷解，增加对/>的约束，使/>解唯一，即：
（11）
其中，表示正则化算子。
（12）
基于式（12）可获得域系数/>。
事实上，当式（6）的为Radon变换，式（6）即可实现Radon变换约束的多道叠前反褶积。显然，/>也可以为FK变换、曲波变换等其它稀疏促进变换，其可用迭代软阈值、交替方向乘子法求解式（6）。无论是哪种算法，均需对变换域数据施加阈值处理，以压制算法不稳定导致的噪声和数据本身的噪声，而/>域具有能量聚集性，有效信号聚集在同速度相关的小曲率位置，多次波聚集在相对较大曲率的位置，随机噪声分布在整个Radon域，直接应用常数阈值（软阈值、硬阈值等）会损失部分有效信号，且多次波也具有能量聚集性，也会残余部分多次波。为此，考虑变换域能量聚集性，优化常数阈值，使阈值同曲率、时间相关，即对于每个曲率有：
（13）
（14）
其中，表示振幅，定义为阈值大小，同信噪比相关，一般随时间增加而增大；/>为曲率变量；/>表示一次波聚集的曲率位置，可大致估算，同速度相关，一般随时间增加而增大；/>为窗的大小，一般需要包含一次波的位置，但随深度增加，信噪比降低，该参数随时间增加而增大。
将式（14）带入式（6）中，为了方便并行，应用稳定的交替方向乘子法（ADMM）优化方程，即可得到叠前反褶积后的CMP道集数据。
为进一步优化技术方案，本实施例还提供了一种基于Radon变换的多道叠前反褶积系统，包括：地震数据获取模块、反褶积处理模块；
地震数据获取模块，用于获取待处理叠前地震数据集；
反褶积处理模块，用于将待处理叠前地震数据集输入预设的反褶积模型中，输出叠前反褶积后的地震数据集，其中，反褶积模型基于训练集训练获得，训练集包括满足拉普拉斯分布的异常振幅噪声和满足高斯分布的随机噪声的CMP道集数据；反褶积模型通过双曲Radon变换约束结合叠前反褶积算法构建。
反褶积模型包括：异常振幅噪声判断单元、第一反褶积单元和第二反褶积单元，异常振幅噪声判断单元用于判断待处理叠前地震数据集中是否存在异常振幅噪声；第一反褶积单元用于当待处理叠前地震数据集中存在异常振幅噪声时，基于第一Radon域叠前反褶积算法对待处理叠前地震数据集进行叠前反褶积处理；第二反褶积单元用于当待处理叠前地震数据集中不存在异常振幅噪声时，基于第二Radon域叠前反褶积算法对待处理叠前地震数据集进行叠前反褶积处理。
为对本实施例所提出的一种基于Radon变换的多道叠前反褶积方法及系统进行验证，下面通过测试进行验证：
测试一：
首先对人工合成数据进行测试。该数据为叠前CMP道集数据，如图2所示，包含306个采样点，采样率为2ms，如图2（a）所示，数据中添加了异常能量噪声，如图2（b）所示。对应图2（a）的Radon变换后的域数据如图2（c）所示，图2（d）是依据图2（c）构建的时变高斯阈值，且考虑深层速度不准确的特征，浅层高斯阈值的窗口宽度较小，深层的窗口宽度较大。
分别基于Tikhonov正则化约束的反褶积（Tikhonov-Decon）；fk变换作为稀疏变换域的反褶积（fk-Decon）；L2范数作为拟合项的Radon变换作为稀疏变换的反褶积（L2-Radon-Decon）和L1范数作为拟合项的反褶积方法（L1-Radon -Decon）对道集进行反褶积测试，数据测试结果如图3所示，其中，图3（a）为吉洪诺夫反褶积后的道集，图3（b）为F-K变换作为稀疏促进变换的反褶积后的道集，图3（c）为L2范数作为拟合项的Radon变换作为稀疏促进变换的反褶积后的道集，图3（d）为L1范数作为拟合项的Radon变换作为稀疏促进变换的反褶积后的道集。显然，由于强振幅噪声的存在，Tikhonov-Decon和fk-Decon算法结果均不稳定，这是由于依据L2范数作为反褶积的拟合项，而噪声不符合高斯分布，使反褶积效果变差；另一方面，Tikhonov正则化的平滑特性和fk域的阈值对振幅较大的异常噪声几乎无效果。然而Radon变换在压制该噪声的方面具有得天独厚的优势，并且应用了L1范数为约束项的L1-Radon-Decon的剖面更规整，同相轴的杂乱现象消失。不同反褶积后的反射系数及地震记录单道对比如图4所示，以及不同反褶积前后振幅谱对比如图5所示，从单道对比图上也可发现，本实施例提出的反褶积方法能有效压制异常噪声（图4虚线框），且这几种方法均有效展宽了地震频带。
依据该道集模型重复测试100次，由于高分辨率处理后数据的真实信噪比难以计算。因此，按（Chen et al., 2022）的方法估计信噪比（SNR），同时计算反褶积后数据同反射系数的结构相似度（SSIM ,Wang et al., 2004）。其中，估计的信噪比表征处理后数据质量及规则化程度，信噪比值具有相对意义；结构相似度表示反褶积能力，处理后数据越接近于反射系数，SSIM的值越接近于1。因为不可能得到全频带的反射系数，因此，本参数也具有相对意义。即，信噪比的值越大，表征处理后数据质量越高；结构相似度值越大，表征反褶积效果越好。不同反褶积稳定性测试结果如图6所示，其中，图6（a）为不同反褶积的信噪估计比，图6（b）为不同反褶积的结构相似度，发现重复测试100次，几种方法的稳定性均较高，且建议的方法信噪比绝对占优，几种方法的SSIM值接近，均高于原始数据，说明反褶积的有效性。
依据上述道集，开展算法抗噪性测试，添加现有噪声的1-50倍能量的噪声，测试结果如图7所示，其中，图7（a）为不同反褶积的信噪估计比，图7（b）为不同反褶积的结构相似度。显然，L1范数作为反褶积的拟合项能有效适用于含Laplace噪声的数据，整体信噪比及反褶积能力基本不受噪声水平的影响，其余几种方法的抗噪能力（SNR）和反褶积水平（SSIM）均随噪声强度增加而降低，在黑线位置有明显变化。
同时，也测试了高斯噪声对算法性能的影响，测试结果与理论相同，如图8所示，其中，图8（a）为不同反褶积的信噪估计比，图8（b）为不同反褶积的结构相似度，L2范数作为拟合项的L2-Radon-Decon整体去噪能力和高分辨率能力均强于其它几种方法。
测试二：
将所提出的基于Radon变换约束的多道叠前反褶积算法应用到实际数据中，分别基于Tikhonov正则化约束的反褶积（Tikhonov-Decon）；fk变换作为稀疏变换域的反褶积（fk-Decon）；L2范数作为拟合项的Radon变换作为稀疏变换的反褶积（L2-Radon-Decon）和L1范数作为拟合项的反褶积方法（L1-Radon -Decon）对道集进行反褶积测试。如图9所示，测试数据有2000个采样点，采样率为2ms，如图9（a）所示，覆盖次数为108次，数据的信噪比随深度逐渐下降。数据第一道的时频谱如图9（b）所示，在全部深度上，主频和频带宽度基本保持稳定，说明数据不存在与地层本身相关的吸收衰减；其域数据如图9（c）所示，中浅层聚焦性较好，深层聚焦性较差，信噪比较低；图9（d）是依据式（13）和式（14）计算得到的时变阈值。
应用测试一中的几种测试方法开展反褶积算法的测试，测试结果如图10、图11所示，其中，图10（a）为吉洪诺夫反褶积后的数据，图10（b）为F-K变换作为稀疏促进变换的反褶积后的数据，图10（c）为L2范数作为拟合项的Radon变换作为稀疏促进变换的反褶积后的数据，图10（d）为L1范数作为拟合项的Radon变换作为稀疏促进变换的反褶积后的数据。参与测试的几种算法均实现了频带宽度的提高，进而压缩地震子波，以提高分辨率（如图10的中浅层所示，图11可观测到频带宽度和主频的提升），且Radon变换作为稀疏促进变换的两种反褶积算法，反褶积后数据的信噪比均得到有效改善。
以上所述的实施例仅是对本发明优选方式进行的描述，并非对本发明的范围进行限定，在不脱离本发明设计精神的前提下，本领域普通技术人员对本发明的技术方案做出的各种变形和改进，均应落入本发明权利要求书确定的保护范围内。
Claims (8)
1.一种基于Radon变换的多道叠前反褶积方法，其特征在于，包括：
获取待处理叠前地震数据集；
将所述待处理叠前地震数据集输入预设的反褶积模型中，输出叠前反褶积后的地震数据集，其中，所述反褶积模型基于训练集训练获得，所述训练集包括满足拉普拉斯分布的异常振幅噪声和满足高斯分布的随机噪声的CMP道集数据；所述反褶积模型通过双曲Radon变换约束结合叠前反褶积算法构建。
2.根据权利要求1所述的基于Radon变换的多道叠前反褶积方法，其特征在于，所述反褶积模型包括：异常振幅噪声判断单元、第一反褶积单元和第二反褶积单元，所述异常振幅噪声判断单元用于判断所述待处理叠前地震数据集中是否存在异常振幅噪声；所述第一反褶积单元用于当所述待处理叠前地震数据集中存在异常振幅噪声时，基于第一Radon域叠前反褶积算法对所述待处理叠前地震数据集进行叠前反褶积处理；所述第二反褶积单元用于当所述待处理叠前地震数据集中不存在异常振幅噪声时，基于第二Radon域叠前反褶积算法对所述待处理叠前地震数据集进行叠前反褶积处理。
3.根据权利要求2所述的基于Radon变换的多道叠前反褶积方法，其特征在于，所述第一Radon域叠前反褶积算法以L1范数作为拟合项，所述第二Radon域叠前反褶积算法以L2范数作为拟合项。
4.根据权利要求2所述的基于Radon变换的多道叠前反褶积方法，其特征在于，所述第一Radon域叠前反褶积算法和第二Radon域叠前反褶积算法为：
其中，R为多道反射系数，W为地震子波构成的卷积矩阵，S为多道地震数据，λ为正则化算子，B为双曲Radon变换，|| ||1为数据绝对值的和，|| ||2为数据平方根的和。
5.根据权利要求4所述的基于Radon变换的多道叠前反褶积方法，其特征在于，所述第一Radon域叠前反褶积算法和第二Radon域叠前反褶积算法的求解方法为：
设置Radon域数据高斯阈值，以及迭代次数；
基于所述高斯阈值，结合迭代阈值收缩算法对所述待处理叠前地震数据集在Radon域进行噪声压制，直至达到所述迭代次数，获取压制噪声后的Radon域叠前地震数据；
对压制噪声后的Radon域叠前地震数据进行Radon反变换，获取变换域压制噪声后的叠前地震数据集，对所述变换域压制噪声后的叠前地震数据集进行叠前反褶积处理，获取叠前反褶积后的地震数据集。
6.根据权利要求5所述的基于Radon变换的多道叠前反褶积方法，其特征在于，所述Radon域数据高斯阈值的计算方法为：
其中，H为阈值，/>为振幅，h gauss 为当前曲率下构建的高斯窗函数，/>为曲率变量，/>为一次波聚集的曲率位置，/>为窗的大小，为当前计算时刻。
7.一种基于Radon变换的多道叠前反褶积系统，其特征在于，用于实施如权利要求1-6任一项所述的基于Radon变换的多道叠前反褶积方法，所述系统包括：地震数据获取模块、反褶积处理模块；
所述地震数据获取模块，用于获取待处理叠前地震数据集；
所述反褶积处理模块，用于将所述待处理叠前地震数据集输入预设的反褶积模型中，输出叠前反褶积后的地震数据集，其中，所述反褶积模型基于训练集训练获得，所述训练集包括异常振幅噪声和随机噪声的CMP道集数据；所述反褶积模型通过双曲Radon变换约束结合叠前反褶积算法构建。
8.根据权利要求7所述的基于Radon变换的多道叠前反褶积系统，其特征在于，所述反褶积模型包括：异常振幅噪声判断单元、第一反褶积单元和第二反褶积单元，所述异常振幅噪声判断单元用于判断所述待处理叠前地震数据集中是否存在异常振幅噪声；所述第一反褶积单元用于当所述待处理叠前地震数据集中存在异常振幅噪声时，基于第一Radon域叠前反褶积算法对所述待处理叠前地震数据集进行叠前反褶积处理；所述第二反褶积单元用于当所述待处理叠前地震数据集中不存在异常振幅噪声时，基于第二Radon域叠前反褶积算法对所述待处理叠前地震数据集进行叠前反褶积处理。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202410190044.9A CN117741750B (zh) | 2024-02-21 | 2024-02-21 | 一种基于Radon变换的多道叠前反褶积方法及系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202410190044.9A CN117741750B (zh) | 2024-02-21 | 2024-02-21 | 一种基于Radon变换的多道叠前反褶积方法及系统 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN117741750A true CN117741750A (zh) | 2024-03-22 |
| CN117741750B CN117741750B (zh) | 2024-04-26 |
Family
ID=90277794
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202410190044.9A Active CN117741750B (zh) | 2024-02-21 | 2024-02-21 | 一种基于Radon变换的多道叠前反褶积方法及系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN117741750B (zh) |
Citations (13)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US5920828A (en) * | 1997-06-02 | 1999-07-06 | Baker Hughes Incorporated | Quality control seismic data processing system |
| CN103954992A (zh) * | 2014-04-03 | 2014-07-30 | 中国石油天然气股份有限公司 | 一种反褶积方法及装置 |
| CN104914466A (zh) * | 2015-06-26 | 2015-09-16 | 中国石油大学(华东) | 一种提高地震资料分辨率的方法 |
| CN107607994A (zh) * | 2017-11-07 | 2018-01-19 | 中国海洋石油总公司 | 一种基于高斯平滑的时频域反褶积方法 |
| CN107678062A (zh) * | 2017-09-15 | 2018-02-09 | 上海海洋大学 | 双曲Radon域综合预测反褶积和反馈循环方法压制多次波模型构建方法 |
| CN113341463A (zh) * | 2021-06-10 | 2021-09-03 | 中国石油大学(北京) | 一种叠前地震资料非平稳盲反褶积方法及相关组件 |
| CN114428348A (zh) * | 2020-10-14 | 2022-05-03 | 中国石油化工股份有限公司 | 一种提高地震资料分辨率的方法及系统 |
| US20220291407A1 (en) * | 2020-12-11 | 2022-09-15 | Institute Of Geology And Geophysics, Chinese Academy Of Sciences | Method and device for elastic parameter inversion |
| EP4080250A1 (en) * | 2021-04-22 | 2022-10-26 | Cgg Services Sas | Multiple attenuation and imaging processes for recorded seismic data |
| CN115685340A (zh) * | 2021-07-21 | 2023-02-03 | 中国石油天然气股份有限公司 | 提高薄层结构地震表征精度的方法及装置 |
| US20230251396A1 (en) * | 2022-02-08 | 2023-08-10 | King Fahd University Of Petroleum And Minerals | Event detection and de-noising method for passive seismic data |
| CN116626765A (zh) * | 2023-06-07 | 2023-08-22 | 电子科技大学 | 一种基于二维k-svd和卷积稀疏编码的多道地震反褶积方法 |
| CN117406272A (zh) * | 2023-10-17 | 2024-01-16 | 中海石油(中国)有限公司 | 一种快速多元信息约束的反褶积宽频处理方法及装置 |
-
2024
- 2024-02-21 CN CN202410190044.9A patent/CN117741750B/zh active Active
Patent Citations (13)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US5920828A (en) * | 1997-06-02 | 1999-07-06 | Baker Hughes Incorporated | Quality control seismic data processing system |
| CN103954992A (zh) * | 2014-04-03 | 2014-07-30 | 中国石油天然气股份有限公司 | 一种反褶积方法及装置 |
| CN104914466A (zh) * | 2015-06-26 | 2015-09-16 | 中国石油大学(华东) | 一种提高地震资料分辨率的方法 |
| CN107678062A (zh) * | 2017-09-15 | 2018-02-09 | 上海海洋大学 | 双曲Radon域综合预测反褶积和反馈循环方法压制多次波模型构建方法 |
| CN107607994A (zh) * | 2017-11-07 | 2018-01-19 | 中国海洋石油总公司 | 一种基于高斯平滑的时频域反褶积方法 |
| CN114428348A (zh) * | 2020-10-14 | 2022-05-03 | 中国石油化工股份有限公司 | 一种提高地震资料分辨率的方法及系统 |
| US20220291407A1 (en) * | 2020-12-11 | 2022-09-15 | Institute Of Geology And Geophysics, Chinese Academy Of Sciences | Method and device for elastic parameter inversion |
| EP4080250A1 (en) * | 2021-04-22 | 2022-10-26 | Cgg Services Sas | Multiple attenuation and imaging processes for recorded seismic data |
| CN113341463A (zh) * | 2021-06-10 | 2021-09-03 | 中国石油大学(北京) | 一种叠前地震资料非平稳盲反褶积方法及相关组件 |
| CN115685340A (zh) * | 2021-07-21 | 2023-02-03 | 中国石油天然气股份有限公司 | 提高薄层结构地震表征精度的方法及装置 |
| US20230251396A1 (en) * | 2022-02-08 | 2023-08-10 | King Fahd University Of Petroleum And Minerals | Event detection and de-noising method for passive seismic data |
| CN116626765A (zh) * | 2023-06-07 | 2023-08-22 | 电子科技大学 | 一种基于二维k-svd和卷积稀疏编码的多道地震反褶积方法 |
| CN117406272A (zh) * | 2023-10-17 | 2024-01-16 | 中海石油(中国)有限公司 | 一种快速多元信息约束的反褶积宽频处理方法及装置 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN117741750B (zh) | 2024-04-26 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108267784A (zh) | 一种地震信号随机噪声压制处理方法 | |
| CN110174702B (zh) | 一种海上地震数据低频弱信号恢复的方法和系统 | |
| CN107247290B (zh) | 一种基于时空分数阶滤波的地震资料噪声压制方法 | |
| CN109738951B (zh) | 一种基于地震同相轴子波谱的时变反褶积方法 | |
| Li et al. | Wavelet-based higher order correlative stacking for seismic data denoising in the curvelet domain | |
| CN113887398A (zh) | 一种基于变分模态分解和奇异谱分析的gpr信号去噪方法 | |
| CN113077386A (zh) | 基于字典学习和稀疏表征的地震资料高分辨率处理方法 | |
| CN107607994B (zh) | 一种基于高斯平滑的时频域反褶积方法 | |
| CN110244353A (zh) | 一种基于稀疏范数优化算法的地震数据规则化方法 | |
| CN111736224B (zh) | 一种压制叠前地震资料线性干扰方法、存储介质及设备 | |
| CN104635264B (zh) | 叠前地震数据的处理方法及设备 | |
| Fu et al. | 3-D structural complexity-guided predictive filtering: A comparison between different non-stationary strategies | |
| CN113156514B (zh) | 基于主频波数域均值滤波的地震数据去噪方法及系统 | |
| CN117741750B (zh) | 一种基于Radon变换的多道叠前反褶积方法及系统 | |
| CN112213775B (zh) | 一种高覆盖次数叠前地震数据的保真提频方法 | |
| CN111257931B (zh) | 一种去除海洋地震勘探过船干扰噪音的方法 | |
| Nose-Filho et al. | Algorithms for sparse multichannel blind deconvolution | |
| CN111273351B (zh) | 用于地震资料去噪的结构导向方向广义全变差正则化方法 | |
| CN109782346B (zh) | 一种基于形态成分分析的采集脚印压制方法 | |
| CN114428343A (zh) | 基于归一化互相关的Marchenko成像方法及系统 | |
| Shuchong et al. | Seismic signals wavelet packet de-noising method based on improved threshold function and adaptive threshold | |
| CN110703332B (zh) | 一种鬼波压制方法 | |
| Renfei et al. | Method for wavelet denoising of multi-angle prestack seismic data | |
| CN110749923A (zh) | 一种基于范数方程提高分辨率的反褶积方法 | |
| CN114137606A (zh) | 一种稳健的谱模拟反褶积方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |