CN117055116B - 一种强稀疏约束Radon变换多次波压制方法及系统 - Google Patents
一种强稀疏约束Radon变换多次波压制方法及系统 Download PDFInfo
- Publication number
- CN117055116B CN117055116B CN202311315020.3A CN202311315020A CN117055116B CN 117055116 B CN117055116 B CN 117055116B CN 202311315020 A CN202311315020 A CN 202311315020A CN 117055116 B CN117055116 B CN 117055116B
- Authority
- CN
- China
- Prior art keywords
- radon
- representing
- data
- transform
- norm
- 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.)
- Active
Links
- 229910052704 radon Inorganic materials 0.000 title claims abstract description 249
- SYUHGPGVQRZVTB-UHFFFAOYSA-N radon atom Chemical compound [Rn] SYUHGPGVQRZVTB-UHFFFAOYSA-N 0.000 title claims abstract description 249
- 230000009466 transformation Effects 0.000 title claims abstract description 103
- 238000000034 method Methods 0.000 title claims abstract description 33
- 230000001629 suppression Effects 0.000 title claims description 36
- 238000003825 pressing Methods 0.000 claims abstract description 18
- 238000005520 cutting process Methods 0.000 claims abstract description 15
- 239000011159 matrix material Substances 0.000 claims description 33
- 238000004422 calculation algorithm Methods 0.000 claims description 31
- 238000004321 preservation Methods 0.000 claims description 14
- 238000012937 correction Methods 0.000 claims description 10
- 238000010276 construction Methods 0.000 claims description 8
- 230000008569 process Effects 0.000 claims description 7
- 230000002087 whitening effect Effects 0.000 claims description 6
- 238000005070 sampling Methods 0.000 claims description 2
- 238000006243 chemical reaction Methods 0.000 abstract description 4
- 230000006870 function Effects 0.000 description 14
- 238000010586 diagram Methods 0.000 description 11
- 238000012545 processing Methods 0.000 description 9
- 238000004364 calculation method Methods 0.000 description 7
- 230000000694 effects Effects 0.000 description 6
- 230000010354 integration Effects 0.000 description 6
- 230000002776 aggregation Effects 0.000 description 5
- 238000004220 aggregation Methods 0.000 description 5
- 238000003491 array Methods 0.000 description 4
- 238000001914 filtration Methods 0.000 description 4
- 238000011160 research Methods 0.000 description 4
- 239000013598 vector Substances 0.000 description 4
- 238000005056 compaction Methods 0.000 description 3
- 238000000605 extraction Methods 0.000 description 3
- 238000002679 ablation Methods 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 2
- 230000000007 visual effect Effects 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
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 239000002131 composite material Substances 0.000 description 1
- 238000007906 compression Methods 0.000 description 1
- 230000006835 compression Effects 0.000 description 1
- 238000012669 compression test Methods 0.000 description 1
- 238000013170 computed tomography imaging Methods 0.000 description 1
- 230000008878 coupling Effects 0.000 description 1
- 238000010168 coupling process Methods 0.000 description 1
- 238000005859 coupling reaction Methods 0.000 description 1
- 238000013144 data compression Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 230000002427 irreversible effect Effects 0.000 description 1
- 238000013507 mapping Methods 0.000 description 1
- 238000005065 mining Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 238000005498 polishing Methods 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 238000013139 quantization Methods 0.000 description 1
- 238000000926 separation method Methods 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 238000001228 spectrum Methods 0.000 description 1
- 238000003860 storage Methods 0.000 description 1
- 239000000126 substance Substances 0.000 description 1
- 238000003786 synthesis reaction Methods 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 238000011426 transformation method Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/32—Transforming one recording into another or one representation into another
- G01V1/325—Transforming one representation into another
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
- G01V1/362—Effecting static or dynamic corrections; Stacking
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
- G06F17/141—Discrete Fourier transforms
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/40—Transforming data representation
- G01V2210/46—Radon transform
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/40—Transforming data representation
- G01V2210/48—Other transforms
Landscapes
- Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Data Mining & Analysis (AREA)
- Remote Sensing (AREA)
- Life Sciences & Earth Sciences (AREA)
- Theoretical Computer Science (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Optimization (AREA)
- Mathematical Analysis (AREA)
- Computational Mathematics (AREA)
- Environmental & Geological Engineering (AREA)
- Algebra (AREA)
- Geophysics (AREA)
- General Life Sciences & Earth Sciences (AREA)
- Geology (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Acoustics & Sound (AREA)
- Discrete Mathematics (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本申请公开了一种强稀疏约束Radon变换多次波压制方法及系统，方法包括以下步骤：对原始CMP数据做快速傅里叶变换，得到频率‑空间域数据；对所述频率‑空间域数据做Lp‑1范数高分辨率Radon变换，得到Radon域数据最优解；利用多次波能量团切除函数，切除所述Radon域数据最优解中的多次波能量团，得到压制多次波后的Radon域数据；对所述压制多次波后的Radon域数据做保幅Radon反变换，得到压制多次波后CMP数据。本申请相比于传统的基于L1范数的高分辨率Radon变换，其精度更高、稀疏度更强，多次波能量在Radon域的聚焦性更好，提高了Radon域压制多次波的精度。
Description
技术领域
本申请属于地震数据处理技术领域，具体涉及一种强稀疏约束Radon变换多次波压制方法及系统。
背景技术
地震资料处理中，由于地震数据的复杂性，仅时空域的处理方法早已无法满足资料质量要求，因此，众多变换域方法被应用到地震数据处理中，包括Fourier变换、KL变换、小波变换、Radon变换等。其中，Radon变换依据视速度差异分离有效信号和噪声，利用线积分的能量，恢复路径上的信息，该变换常用于医学CT成像，图像处理与增强等领域，在地震数据处理中，考虑CMP道集的多次波和一次波视速度的差异，可完成多次波压制；依据全局变换的特征，也可完成缺失数据重建等工作。
Radon变换已经于1986年引入图像处理中，被用于线性特征的检测与增强（LesleyM）；Douglas J提出双曲Radon变换压制数据的多次波；Sacchi、Maria S和M.M. NurulKabir同时在1995年利用Radon变换完成缺失数据的重建，Herrmann（2000）利用Radon域的稀疏特性，分离混采数据的离散副震源；在这些学者的研究基础上，考虑Radon逆变换计算效率和不同线积分路径问题，Sacchi（1999）提出最小二乘抛物Radon变换，能有效提高计算效率，类似的研究工作还包括Brady（1998）、Daniel（2002）、Amir（2003）、Hu（2013）等。
除计算效率等问题，考虑数据的AVO特征，基于能量重分配的Radon变换已被提出和广泛应用。而变换域的能量聚集性关系到数据的处理精度，因此提出高分辨率稀疏Radon变换，这类变换利用稀疏范数约束Radon变换域，提高射线参数方向的能量聚集性，依据线积分理论，射线参数方向分辨率越高，有效波和多次波的分离效果越高，基于该思想，依据压缩感知反演理论，众多稀疏范数被用于Radon变换的成像工作中，理论上，L0范数约束是压缩感知算法的获得精确解的必要条件，但L0范数的求解是NP-hard问题，因此，该范数的高精度近似是研究的热点问题，目前包括L1范数、L1-L2范数等。
传统的三维抛物Radon变换，大多是采用L2范数或者L1范数作为正则化条件。基于L2范数的三维抛物Radon变换在Radon域中存在严重的能量拖尾现象，多次波能量与有效波能量在对应的曲率上聚焦性较差，部分能量交织在一起难以分离，这严重影响该方法在多次波压制中的应用；基于L1范数的高精度三维抛物Radon变换提高了信号在变换域中的精度，在一定程度上缓解了一次波和多次波由于截距时间较小产生的能量团拖尾现象，但其仍然存在稀疏约束能力不足的情况，当多次波曲率较低时、多次波与有效波混叠时，基于L1范数的高精度三维抛物Radon变换方法在多次波压制中应用效果有限。因此，针对目前广泛使用的Radon变换压制多次波的算法，急需解决普遍存在的变换域稀疏约束能力不足的问题。
发明内容
本申请旨在解决现有技术的不足，提出一种强稀疏约束Radon变换多次波压制方法及系统，针对三维地震数据动校正（NMO）后有效波能量轴与多次波能量轴的差异，采用Lp-1范数强稀疏约束的高分辨率保幅抛物Radon变换，使得动校正后CMP道集中的多次波能量在Radon域聚焦，将聚焦的多次波能量团切除后，对剩余数据做Radon反变换，最终得到压制多次波后的CMP道集数据。
为实现上述目的，本申请提供了如下方案：
一种强稀疏约束Radon变换多次波压制方法，包括以下步骤：
对原始CMP数据做快速傅里叶变换，得到频率-空间域数据；
对所述频率-空间域数据做Lp-1范数高分辨率Radon变换，得到Radon域数据最优解；
利用多次波能量团切除函数切除所述Radon域数据最优解中的多次波能量团，得到压制多次波后Radon域数据；
对所述压制多次波后Radon域数据做保幅Radon反变换，得到压制多次波后CMP数据。
优选的，所述Lp-1范数高分辨率Radon变换为：
其中，表示Radon域数据最优解，/>和/>分别表示不同测线方向的保幅Radon变换算子，/>表示观测地震数据，/>表示傅里叶变换，/>表示离散逆傅里叶变换矩阵，/>表示Radon域数据，/>表示正则化算子，/>表示Lp-1范数中L1范数的权重，/>表示Lp范数，表示L1范数；
所述Radon域数据为：
其中，x和y分别表示沿不同测线方向的炮检距，表示频率，/>表示虚数单位，表示频率空间域地震数据，/>、/>分别表示不同测线方向的Radon变换算子，和/>分别表示不同空间方向的曲率，/>和/>表示不同空间方向的采样点数。
优选的，所述不同测线方向的保幅Radon变换算子为：
其中，、/>分别表示/>、/>的矩阵形式，/>表示x方向保幅算子，/>表示y方向保幅算子，/>、分别表示x、y方向二维时空域地震数据，N表示地震道数，/>、/>分别表示/>、/>方向第/>道地震数据的第j阶正交多项式系数，H表示矩阵的复共轭，/>、/>表示保幅处理前的Radon变换算子，/>、/>分别表示x、y方向的保幅空间矫正算子，为：
其中，表示原始CMP数据的能量阈值，/>与/>构造方式相同。
优选的，通过设置迭代重加权最小二乘算法来求解所述Radon域数据最优解：
其中，k表示迭代次数，、/>分别表示第k次迭代更新的x、y方向的加权矩阵，b表示预白化因子，p表示Lp范数的幂次项，/>表示将括号内的元素排列成一个对角矩阵。
优选的，所述保幅Radon反变换为：
其中，表示反Radon变换得到的频率域地震数据，/>、，分别表示不同测线方向的Radon变换算子，/>表示频率，/>表示虚数单位。
本申请还提供了一种强稀疏约束Radon变换多次波压制系统，应用于上述任一项所述的方法，包括：傅里叶变换模块、Radon变换模块、多次波压制模块和Radon反变换模块；
所述傅里叶变换模块用于对原始CMP数据做快速傅里叶变换，得到频率-空间域数据；
所述Radon变换模块用于对所述频率-空间域数据做Lp-1范数高分辨率Radon变换，得到Radon域数据最优解；
所述多次波压制模块用于利用多次波能量团切除函数切除所述Radon域数据最优解中的多次波能量团，得到压制多次波后Radon域数据；
所述Radon反变换模块用于对所述压制多次波后Radon域数据做保幅Radon反变换，得到压制多次波后CMP数据。
优选的，所述Lp-1范数高分辨率Radon变换为：
其中，表示Radon域数据最优解，/>和/>分别表示不同测线方向的保幅Radon变换算子，/>表示观测地震数据，/>表示傅里叶变换，/>表示离散逆傅里叶变换矩阵，/>表示Radon域数据，/>表示正则化算子，/>表示Lp-1范数中L1范数的权重，/>表示Lp范数，/>表示L1范数；
所述Radon域数据为：
其中，x和y分别表示沿不同测线方向的炮检距，表示频率，/>表示虚数单位，表示频率空间域地震数据，/>、/>分别表示不同测线方向的Radon变换算子，和/>分别表示不同空间方向的曲率，/>和/>表示不同空间方向的采样点数。
优选的，所述不同测线方向的保幅Radon变换算子为：
其中，、/>分别表示/>、/>的矩阵形式，/>表示x方向保幅算子，/>表示y方向保幅算子，/>、分别表示x、y方向二维时空域地震数据，N表示地震道数，/>、/>分别表示/>、/>方向第/>道地震数据的第j阶正交多项式系数，H表示矩阵的复共轭，/>、/>表示保幅处理前的Radon变换算子，/>、/>分别表示x、y方向的保幅空间矫正算子，为：
其中，表示原始CMP数据的能量阈值，/>与/>构造方式相同。
优选的，通过设置迭代重加权最小二乘算法来求解所述Radon域数据最优解：
其中，k表示迭代次数，、/>分别表示第k次迭代更新的x、y方向的加权矩阵，b表示预白化因子，p表示Lp范数的幂次项，/>表示将括号内的元素排列成一个对角矩阵。
优选的，所述Radon反变换模块的工作流程包括：
其中，表示反Radon变换得到的频率域地震数据，/>、，分别表示不同测线方向的Radon变换算子，/>表示频率，/>表示虚数单位。
与现有技术相比，本申请的有益效果为：
本申请相比于传统的基于L1范数的高分辨率Radon变换，其精度更高、稀疏度更强，多次波能量在Radon域的聚焦性更好，提高了Radon域压制多次波的精度。由于正、反Radon变换不可逆，且传统Radon变换不具备保幅性，本申请用频率域迭代重加权最小二乘算法（IRLS）求解Radon变换算子，并将其与基于多项式拟合的地震信号保幅算法相结合，实现基于Lp-1范数的三维高分辨率保幅抛物Radon变换算法。本申请在算法上的集成性较好，可以实现模块化多次波压制，便于地震资料处理人员的使用，具有较强的推广价值。
附图说明
为了更清楚地说明本申请的技术方案，下面对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅是本申请的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1为本申请实施例的方法流程示意图；
图2为本申请实施例的方法流程框图；
图3为本申请实施例的不同范数约束下公式（10）解的相平面图，其中，（a）为L2范数，（b）为L1范数，（c）为Lp范数，（d）为Lp-1范数；
图4为本申请实施例的第一象限内不同正则化算子解区间示意图，其中，（a）为L2范数，（b）为L1范数，（c）为Lp范数，（d）为Lp-1范数；
图5为本申请实施例的正则化算子解区间内面元积分曲线图；
图6为本申请实施例的系统结构示意图；
图7为本申请实施例的三维CMP道集模拟数据示意图；
图8为本申请实施例的模拟数据纵测线方向抽取5个小排列示意图；
图9为本申请实施例的不同算法压制多次波后的小排列示意图，其中，（a）为基于L2范数三维Radon变换算法多次波压制结果，（b）为基于L1范数三维Radon变换算法多次波压制结果，（c）为基于Lp-1范数三维Radon变换算法多次波压制结果；
图10为本申请实施例的不同算法Radon域数据聚焦性示意图，其中，（a）为L2范数Radon域数据，（b）为L1范数Radon域数据，（c）为Lp-1范数Radon域数据，（d）为L2范数Radon域切除的多次波能量团，（e）为L1范数Radon域切除的多次波能量团，（f）为Lp-1范数Radon域切除的多次波能量团；
图11为本申请实施例的图10数据不同算法压制的多次波示意图，其中，（a）为L2范数压制的多次波，（b）为L1范数压制的多次波，（c）为Lp-1范数压制的多次波；
图12为本申请实施例的野外采集三维NMO后的三维数据示意图，其中，（a）为原始三维数据切片，（b）为横测线方向等间距抽取5条小排列，（c）为基于L2范数多次波压制结果，（d）为基于L1范数多次波压制结果，（e）为基于Lp-1范数多次波压制结果。
具体实施方式
下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本申请一部分实施例，而不是全部的实施例。基于本申请中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本申请保护的范围。
为使本申请的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本申请作进一步详细的说明。
实施例一
在本实施例中，如图1、图2所示，一种强稀疏约束Radon变换多次波压制方法，包括以下步骤：
S1.对原始CMP数据做快速傅里叶变换，得到频率-空间域数据。
S2.对频率-空间域数据做Lp-1范数高分辨率Radon变换，得到Radon域数据最优解。
本实施例中，在三维抛物地震数据处理中，三维抛物Radon正变换表示将地震数据沿着抛物面积分，得到域的变换系数。三维抛物Radon正变换公式为：
（1）
其中，是时空域三维地震数据；t表示时间；/>是Radon域零炮检距截距时间；/>为Radon域变换系数；x、y分别表示沿不同测线方向的炮检距；/>和/>为不同方向的曲率。
将式（1）离散化，等式两边延时间方向进行离散傅里叶变换（DFT），并将其改写为矩阵形式，考虑存储量及计算量的问题，将x和y不同测线方向Radon变换算子在频率域解耦，即得到Radon域数据：
（2）
其中，x和y分别表示沿不同测线方向的炮检距，表示频率空间域地震数据，/>和/>表示不同空间方向的采样点数，/>和/>分别表示不同空间方向的曲率，、/>分别表示不同测线方向的Radon变换算子，/>表示频率，/>表示虚数单位。
由上式知，当已知Radon域系数，将其左乘/>和右乘/>，对频率方向做傅里叶逆变换，可得到时空域地震数据。当已知频率空间域地震数据，可基于式（2）求解Radon域系数/>，其最小二乘解为：
（3）
其中，表示离散傅里叶变换矩阵，/>表示离散逆傅里叶变换矩阵，由于式（3）中/>和/>非满秩，导致/>有无穷解，以/>的平滑性作为约束，即以/>范数约束式（4）：
（4）
其中，表示正则化算子，/>越大Radon域系数/>越平滑，反之/>越小，求解稳定性越差。而当考虑/>中能量的聚集性与稀疏性时，一般情况下，基于L1范数约束式（3），可得到高分辨率三维Radon正变换：
（5）
然而，理论上高分辨率Radon变换应基于L0范数约束变换域系数，其公式为：
（6）
式中，表示L0范数。由于L0范数求解困难（NP-hard问题），考虑采用L0范数的凸近似（L1范数）代替L0范数约束求解/>。
在非凸优化问题求解中，虽然大多数情况下L1范数收敛性较好，但其仍然存在稀疏度不足的问题。针对该问题，稀疏性更强的Lp范数是一个较好的解决方案，其表达式为：
（7）
其中，p表示Lp范数的幂次项，取值范围应为[0，1]。当p为0时，式（7）表示为L0范数；当p为1时，式（7）表示为L1范数。
为了分析L2范数、L1范数、Lp范数与Lp-1范数，在Radon变换应用中的稀疏约束能力，本实施例在单位区间内构建一组空间向量：
（8）
此时，应用上述不同范数约束求解向量，可以表示为：
（9）
将式（9）在单位区间内的解投影到坐标轴上，取p=0.5，得到不同范数约束下向量解区间的相平面图，如图3所示。图3可以定性描述不同范数约束条件下目标函数解的稀疏性。其中，图3（a）-（d）表示目标函数解的稀疏性依次增强，即：正则化算子的稀疏性由弱到强依次为：L2范数、L1范数、Lp范数与Lp-1范数。
为了定量分析不同范数约束条件下目标函数解的稀疏性及其参数取值，我们将图3中不同范数约束向量的解空间，映射到2D坐标轴第一象限中的单位区间内，可以得到如图4所示第一象限内不同正则化算子解区间示意图，其解区间对应的函数，如式（10）所示。
（10）
理论上，图4中曲线与坐标轴之间面积越小，对应的正则化算子稀疏性越强。为了定量分析不同正则化算子的稀疏性，在图4中我们分别对不同正则化算子取任意值，计算/>区间内式（10）的面元积分，得到不同正则化算子的解区间内面元积分曲线，如图5所示。为了便于计算，本实施例中给定不同/>取值的离散采样间隔为0.01。
其中，图5横坐标表示取值，纵坐标表示计算/>区间内式（10）的面元积分。分析图5可以从量化的角度确认，正则化算子的稀疏性由弱到强对应图5中四条曲线由上到下，依次为：L2范数、L1范数、Lp范数与Lp-1范数。
考虑Lp范数应用在三维高精度抛物Radon变换中，能够增强地震数据在域中能量的聚集性和稀疏性，有效提高三维抛物Radon变换的分辨率。同时，结合对基于L1范数约束的高分辨率Radon变换的研究，进一步提出了稀疏性更强且精度更高的Lp-1范数，及Lp-1范数高分辨率Radon变换为，并将其应用到三维抛物Radon变换中，如式（11）所示。
（11）
其中，表示Radon域数据最优解，/>和/>分别表示不同测线方向的保幅Radon变换算子，/>表示观测地震数据，/>表示傅里叶变换，/>表示Radon域数据，/>表示正则化算子，/>表示Lp-1范数中L1范数的权重，/>表示Lp范数，/>表示L1范数；
传统Radon变换在计算中未考虑时空域振幅随偏移距的变化特征，导致Radon变换存在不保幅的问题。为了解决这个问题，可以采用正交多项式拟合（OPT），将其应用到三维保幅Radon变换中。
OPT保幅的实现原理是基于这样一个假设：空间上一组地震数据的振幅随偏移距的变化，可以用一个偶次多项式拟合得到。因此，地震数据在偏移坐标上可以展开为一组正交多项式的线性组合，即：
（12）
其中，是二维时空域地震数据；/>是t时刻第j阶多项式拟合系数，N表示地震道数；/>为x方向上第j阶正交多项式。
由于正交多项式满足正交性条件，则可以得到拟合系数的解为：
（13）
通过计算不同时间点上的拟合系数，可以得到OPT谱，即/>随时间和阶次的变化。最后，采用拟合系数/>和正交多项式/>来重构地震数据。
针对NMO校正后的CMP数据，结合OPT保幅思想，创建三维Radon变换保幅算子。考虑到Radon变换后CMP道集数据在空间上会得到补偿，这将在CMP道集中引入假的能量。为了实现CMP道集数据在空间方向的高分辨率保幅，我们在常规OPT保幅算法上进行改进。由于不同测线方向保幅算子构造方式相同，这里以x方向为例，引入空间矫正算子，得到CMP数据空间方向t时刻的高分辨率保幅算子，即：
（14）
其中，表示t时刻原始CMP数据的x空间特征，可以由下式构建：
（15）
其中，表示原始CMP数据的能量阈值，/>与/>构造方式相同，其余参数含义与上述相同。
通常情况下，采用3阶正交多项式系数，即可实现Radon变换中时空域地震数据空间方向的保幅。其中，表示叠加振幅信息、/>表示梯度信息、/>表示曲率信息。对于三维保幅Radon变换，需要分别求取x、y方向的正交多项式系数，将式（14）保幅算子映射到频率域，与Radon变换算子相结合，得到如式（16）不同测线方向的保幅Radon变换算子：
（16）
将式（16）代入式（12），可以实现强稀疏约束Lp-1范数高分辨率保幅三维抛物Radon变换。
在本实施例中，为了求解式（11），将频率域迭代重加权最小二乘算法（IRLS）应用到Lp-1范数约束的三维抛物Radon变换求解中，建立了Lp-1范数约束下的迭代求解方程最优解，如式（17）所示。
（17）
其中，k表示迭代次数，、/>分别表示第/>次迭代更新的x、y方向的加权矩阵，b表示预白化因子，p表示Lp范数的幂次项，H表示矩阵的复共轭，/>表示将括号内的元素排列成一个对角矩阵。
IRLS算法本质上是解决加权最小二乘问题，随着迭代更新加权矩阵，当算法收敛后，可获得三维高精度Radon域系数矩阵/>。对于Lp-1范数约束的三维抛物Radon变换的求解，/>中每个频率索引i下的Radon域系数，其第k+1次迭代的解/>，可以采用式（17）分别迭代求解。
S3.利用多次波能量团切除函数切除Radon域数据最优解中的多次波能量团，得到压制多次波后Radon域数据。
在本实施例中，由于一次波在Radon域表征为零曲率，多次波在Radon域表征为非零曲率，根据此区别，可以人为规避零曲率附近的Radon域数据，对余下数据按照能量大小构造Radon域切除函数。即，
（18）
其中，表示压制多次波能量团后的Radon域数据，A表示切除函数所保留的曲率范围，B表示切除Radon域多次波能量团的能量阈值。最后，令/>，即可实现Radon域切除多次波能量团。
S4.对压制多次波后Radon域数据做保幅Radon反变换，得到压制多次波后CMP数据。
在本实施例中，反变换是通过对域系数反演，最终得到时空域地震数据。反变换公式为
（19）
其中，是时空域三维地震数据。将式（19）离散化，等式两边沿时间方向进行离散傅里叶变换（DFT），考虑存储量及计算量的问题，将式（19）改写为矩阵形式，并将x、y不同方向的Radon变换算子在频率域解耦，得到：
（20）
其中，表示反Radon变换得到的频率域地震数据，/>、，分别表示不同测线方向的Radon变换算子。/>表示频率，/>表示虚数单位。
实施例二
在本实施例中，如图6所示，一种强稀疏约束Radon变换多次波压制系统，包括：傅里叶变换模块、Radon变换模块、多次波压制模块和Radon反变换模块；
傅里叶变换模块用于对原始CMP数据做快速傅里叶变换，得到频率-空间域数据。
Radon变换模块用于对频率-空间域数据做Lp-1范数高分辨率Radon变换，得到Radon域数据最优解。Radon变换模块的工作流程包括：
（21）
其中，表示Radon域数据最优解，/>和/>分别表示不同测线方向的保幅Radon变换算子，/>表示观测地震数据，/>表示傅里叶变换，/>表示离散逆傅里叶变换矩阵，/>表示Radon域数据，/>表示正则化算子，/>表示Lp-1范数中L1范数的权重，/>表示Lp范数，/>表示L1范数。
所述Radon域数据为：
（22）
其中，x和y分别表示沿不同测线方向的炮检距，表示频率，/>表示虚数单位，表示频率空间域地震数据，/>、/>分别表示不同测线方向的Radon变换算子，和/>分别表示不同空间方向的曲率，/>和/>表示不同空间方向的采样点数。
不同测线方向的保幅Radon变换算子为：
（23）
其中，、/>分别表示/>、/>的矩阵形式，表示x方向保幅算子，/>表示y方向保幅算子，/>、/>分别表示x、y方向二维时空域地震数据，N表示地震道数，/>、/>分别表示/>、/>方向第/>道地震数据的第j阶正交多项式系数，H表示矩阵的复共轭，/>、/>表示保幅处理前的Radon变换算子，/>、/>分别表示x、y方向的保幅空间矫正算子。由于/>、/>构造方式相同，这里以/>为例，有：
（24）
其中，表示原始CMP数据的能量阈值，/>与/>构造方式相同，其余参数含义与上文相同。
在本实施例中，通过设置迭代重加权最小二乘算法来求解Radon域数据最优解：
（25）
其中，k表示迭代次数，、/>分别表示第k次迭代更新的x、y方向的加权矩阵，b表示预白化因子，p表示/>范数的幂次项，H表示矩阵的复共轭，/>表示将括号内的元素排列成一个对角矩阵，其余参数含义与上文相同。
多次波压制模块用于利用多次波能量团切除函数切除Radon域数据最优解中的多次波能量团，得到压制多次波后Radon域数据。
在本实施例中，由于一次波在Radon域表征为零曲率，多次波在Radon域表征为非零曲率，根据此区别，可以人为规避零曲率附近的Radon域数据，对余下数据按照能量大小构造Radon域切除函数。即，
（18）
其中，表示压制多次波能量团后的Radon域数据，A表示切除函数所保留的曲率范围，B表示切除Radon域多次波能量团的能量阈值。最后，令/>，即可实现Radon域切除多次波能量团。
Radon反变换模块用于对压制多次波后Radon域数据做保幅Radon反变换，得到压制多次波后CMP数据。Radon反变换模块的工作流程包括：
（27）
其中，表示反Radon变换得到的频率域地震数据，/>、，分别表示不同测线方向的Radon变换算子。/>表示频率，/>表示虚数单位。
实施例三
在本实施例中，通过对合成含多次波的三维动校正后模型数据进行测试，验证本申请的有效性，并与基于L2范数、L1范数的三维保幅抛物Radon算法的结果做对比分析。
如图7所示为合成动校正后的CMP三维数据体，该数据横测线方向（X-Distance）共有60条排列，偏移距范围是-1450至1500m，道间距为50m；纵测线方向（Y-Distance）有51条排列，偏移距范围是0至2500m，道间距为50m。为了便于说明，在纵测线方向（800至1800m）偏移距范围间隔250m抽取5个小道集，如图8所示。
图8中①、⑤为动校正后的一次波，②、③、④、⑥、⑦为多次波，两者之间存在着剩余时差。对图7所示模型分别采用基于L2范数的保幅三维抛物Radon变换、基于L1范数的高分辨率保幅三维抛物Radon变换、基于强稀疏约束Lp-1范数的三维保幅抛物Radon变换进行多次波压制测试，压制结果抽取与图8相同的小排列如图9所示，抽取不同算法纵测线方向偏移距为1250m的一个小排列的Radon域数据，如图10a、b、c所示，不同方法在Radon域切除的多次波能量团如图10d、e、f所示。不同方法对应于图8数据压制的多次波，如图11a、b、c所示。
图9中多次波压制后的小排列及图10d、e、f压制的多次波表明，应用我们提出的强稀疏约束Lp-1范数三维保幅抛物Radon变换算法，在最小伤害有效波的同时，多次波压制效果最好。分析图10a、b、c中不同范数约束下的三维Radon域模型可知，基于L2范数的三维保幅抛物Radon变换，其Radon域数据存在能量拖尾现象，有效波与多次波能量交织在一起，不同能量团很难完整区分，严重影响多次波的压制；基于L1范数的高精度三维保幅抛物Radon变换，其在Radon域中能量团实现了收敛，但是当多次波的曲率较低时（如图8中c能量轴所示）、多次波与有效波叠加或距离较近时（如图8中①、②能量轴及⑤、⑥能量轴所示），其Radon域多次波能量团依然会与有效波能量团产生耦合，不容易区分。因此，基于L1范数的高精度保幅三维抛物Radon变换压制多次波，应用效果有限；基于强稀疏约束Lp-1范数三维保幅抛物Radon变换算法，其Radon域数据能量更聚焦，模型分辨率更高，有效波能量团（应对①、②能量轴）基本上都聚焦在零曲率位置。同时，对应于多次波的①、②能量轴及⑤、⑥能量轴实现了完全分离，更易于多次波的压制。图10d、e、f中对应数据压制掉的多次波，如图11a、b、c所示，从图11中可以发现，由于强稀疏约束Lp-1范数三维保幅抛物Radon变换算法增强了变换域的稀疏约束，提高了变换域的能量聚集性，使得该方法在压制多次波时伤害有效波最少，且更易实现Radon域多次波的压制。
实施例四
将所提出的强稀疏约束Lp-1范数三维保幅抛物Radon变换算法应用到实际数据中，并与基于L2范数、L1范数的三维保幅抛物Radon变换算法的结果做对比分析。该实际数据横测线方向（X-Distance）共有60条排列，偏移距范围是-1450至1500m，道间距为50m；纵测线方向（Y-Distance）有51条排列，偏移距范围是-1250至1250m，道间距为50m。为了便于说明，在纵测线方向（750至1750m）偏移距范围间隔250m抽取5个小道集，如图12所示，图中20、25、30、35、40分别表示偏移距由小到大的排列编号。
图12a为三维动校正后的CMP道集，该实际数据横测线方向（X-Distance）共有60条排列，纵测线方向（Y-Distance）有51条排列。为了便于说明，在纵测线方向间隔250m抽取5个道集，如图12b所示，图中20、25、30、35、40分别表示偏移距由小到大的排列编号。从图12b中可以发现，动校正后一次波信号同相轴基本水平，2s到4s内可明显观测有多次波噪声，如图12b箭头所示，有效波与多次波重合，导致有效信号难以识别。将上述三种方法应用到该数据中发现：基于L2范数的三维Radon变换算法压制多次波效果较差，如图12c方框及箭头处所示，有明显的多次波残留。这是由于多次波与有效波能量团在Radon域无法较好的分离引起的；基于L1范数的三维Radon变换算法精度得到提高，但依然无法较好的分离多次波，如图12d中方框内依然存在明显的多次波干扰，图12d箭头处，多次波与有效波叠加在一起，使能量轴产生断裂现象；从图12中（c、d、e）可以发现，基于Lp-1范数的三维Radon变换算法压制了方框内剩余的多次波，对比箭头处可以发现，应用该方法后，同相轴更为连续，压制多次波效果更好。同时，从图12中可以发现，Radon变换算法具有一定的压制随机噪音的能力。
以上所述的实施例仅是对本申请优选方式进行的描述，并非对本申请的范围进行限定，在不脱离本申请设计精神的前提下，本领域普通技术人员对本申请的技术方案做出的各种变形和改进，均应落入本申请权利要求书确定的保护范围内。
Claims (7)
1.一种强稀疏约束Radon变换多次波压制方法，其特征在于，包括以下步骤：
对原始CMP数据做快速傅里叶变换，得到频率-空间域数据；
对所述频率-空间域数据做Lp-1范数高分辨率Radon变换，得到Radon域数据最优解；
利用多次波能量团切除函数切除所述Radon域数据最优解中的多次波能量团，得到压制多次波后Radon域数据；
对所述压制多次波后Radon域数据做保幅Radon反变换，得到压制多次波后CMP数据；
所述Lp-1范数高分辨率Radon变换为：
其中，/>表示Radon域数据最优解，/>和/>分别表示不同测线方向的保幅Radon变换算子，/>表示观测地震数据，/>表示傅里叶变换，/>表示离散逆傅里叶变换矩阵，/>表示Radon域数据，/>表示正则化算子，/>表示Lp-1范数中L1范数的权重，/>表示Lp范数，/>表示L1范数；
所述Radon域数据为：
其中，x和y分别表示沿不同测线方向的炮检距，/>表示频率，/>表示虚数单位，/>表示频率空间域地震数据，/>、/>分别表示不同测线方向的Radon变换算子，/>和/>分别表示不同空间方向的曲率，/>和/>表示不同空间方向的采样点数；
所述不同测线方向的保幅Radon变换算子为：
其中，/>、/>分别表示/>、/>的矩阵形式，表示x方向保幅算子，表示y方向保幅算子，/>、/>分别表示x、y方向二维时空域地震数据，N表示地震道数，/>、/>分别表示/>、/>方向第/>道地震数据的第j阶正交多项式系数，H表示矩阵的复共轭，/>、/>表示保幅处理前的Radon变换算子，/>、/>分别表示x、y方向的保幅空间矫正算子，/>为：
其中，/>表示原始CMP数据的能量阈值，与/>构造方式相同；
通过设置迭代重加权最小二乘算法来求解所述Radon域数据最优解：
其中，k表示迭代次数，、/>分别表示第k次迭代更新的x、y方向的加权矩阵，b表示预白化因子，p表示Lp范数的幂次项，/>表示将括号内的元素排列成一个对角矩阵。
2.根据权利要求1所述一种强稀疏约束Radon变换多次波压制方法，其特征在于，所述保幅Radon反变换为：
其中，/>表示反Radon变换得到的频率域地震数据，/>、/>，分别表示不同测线方向的Radon变换算子，/>表示频率，/>表示虚数单位。
3.一种强稀疏约束Radon变换多次波压制系统，应用于权利要求1-2任一项所述的方法，其特征在于，包括：傅里叶变换模块、Radon变换模块、多次波压制模块和Radon反变换模块；
所述傅里叶变换模块用于对原始CMP数据做快速傅里叶变换，得到频率-空间域数据；
所述Radon变换模块用于对所述频率-空间域数据做Lp-1范数高分辨率Radon变换，得到Radon域数据最优解；
所述多次波压制模块用于利用多次波能量团切除函数切除所述Radon域数据最优解中的多次波能量团，得到压制多次波后Radon域数据；
所述Radon反变换模块用于对所述压制多次波后Radon域数据做保幅Radon反变换，得到压制多次波后CMP数据。
4.根据权利要求3所述一种强稀疏约束Radon变换多次波压制系统，其特征在于，所述Lp-1范数高分辨率Radon变换为：
其中，表示Radon域数据最优解，/>和/>分别表示不同测线方向的保幅Radon变换算子，/>表示观测地震数据，/>表示傅里叶变换，/>表示离散逆傅里叶变换矩阵，/>表示Radon域数据，/>表示正则化算子，/>表示Lp-1范数中L1范数的权重，/>表示Lp范数，/>表示L1范数；
所述Radon域数据为：
其中，x和y分别表示沿不同测线方向的炮检距，表示频率，/>表示虚数单位，表示频率空间域地震数据，/>、/>分别表示不同测线方向的Radon变换算子，和/>分别表示不同空间方向的曲率，/>和/>表示不同空间方向的采样点数。
5.根据权利要求4所述一种强稀疏约束Radon变换多次波压制系统，其特征在于，所述不同测线方向的保幅Radon变换算子为：
其中，/>、/>分别表示/>、/>的矩阵形式，表示x方向保幅算子，表示y方向保幅算子，/>、/>分别表示x、y方向二维时空域地震数据，N表示地震道数，/>、/>分别表示/>、/>方向第/>道地震数据的第j阶正交多项式系数，H表示矩阵的复共轭，/>、/>表示保幅处理前的Radon变换算子，/>、/>分别表示x、y方向的保幅空间矫正算子，/>为：
其中，/>表示原始CMP数据的能量阈值，与/>构造方式相同。
6.根据权利要求5所述一种强稀疏约束Radon变换多次波压制系统，其特征在于，通过设置迭代重加权最小二乘算法来求解所述Radon域数据最优解：
其中，k表示迭代次数，/>、/>分别表示第k次迭代更新的x、y方向的加权矩阵，b表示预白化因子，p表示Lp范数的幂次项，/>表示将括号内的元素排列成一个对角矩阵。
7.根据权利要求6所述一种强稀疏约束Radon变换多次波压制系统，其特征在于，所述Radon反变换模块的工作流程包括：
其中，/>表示反Radon变换得到的频率域地震数据，/>、/>，分别表示不同测线方向的Radon变换算子，/>表示频率，/>表示虚数单位。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311315020.3A CN117055116B (zh) | 2023-10-12 | 2023-10-12 | 一种强稀疏约束Radon变换多次波压制方法及系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311315020.3A CN117055116B (zh) | 2023-10-12 | 2023-10-12 | 一种强稀疏约束Radon变换多次波压制方法及系统 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN117055116A CN117055116A (zh) | 2023-11-14 |
| CN117055116B true CN117055116B (zh) | 2023-12-19 |
Family
ID=88655828
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202311315020.3A Active CN117055116B (zh) | 2023-10-12 | 2023-10-12 | 一种强稀疏约束Radon变换多次波压制方法及系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN117055116B (zh) |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102401908A (zh) * | 2010-09-07 | 2012-04-04 | 中国石油天然气集团公司 | 一种利用不同模加权稀疏的抛物拉东变换压制多次波的方法 |
| CN106443775A (zh) * | 2016-05-25 | 2017-02-22 | 中国石油集团川庆钻探工程有限公司地球物理勘探公司 | 高分辨率转换波裂缝预测方法 |
| CN107102359A (zh) * | 2017-05-18 | 2017-08-29 | 中国石油集团东方地球物理勘探有限责任公司 | 地震数据保幅重建方法和系统 |
| CN109633741A (zh) * | 2019-01-04 | 2019-04-16 | 吉林大学 | 基于双凸优化稀疏约束的混合震源数据一次波分离方法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6763304B2 (en) * | 2002-05-21 | 2004-07-13 | Pgs Americas, Inc. | Method for processing seismic data to attenuate multiples |
-
2023
- 2023-10-12 CN CN202311315020.3A patent/CN117055116B/zh active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102401908A (zh) * | 2010-09-07 | 2012-04-04 | 中国石油天然气集团公司 | 一种利用不同模加权稀疏的抛物拉东变换压制多次波的方法 |
| CN106443775A (zh) * | 2016-05-25 | 2017-02-22 | 中国石油集团川庆钻探工程有限公司地球物理勘探公司 | 高分辨率转换波裂缝预测方法 |
| CN107102359A (zh) * | 2017-05-18 | 2017-08-29 | 中国石油集团东方地球物理勘探有限责任公司 | 地震数据保幅重建方法和系统 |
| CN109633741A (zh) * | 2019-01-04 | 2019-04-16 | 吉林大学 | 基于双凸优化稀疏约束的混合震源数据一次波分离方法 |
Non-Patent Citations (1)
| Title |
|---|
| 混合域高分辨率双曲Radon变换及其在多次波压制中的应用;巩向博;韩立国;王升超;;石油物探(第05期);全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN117055116A (zh) | 2023-11-14 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Hennenfent et al. | Nonequispaced curvelet transform for seismic data reconstruction: A sparsity-promoting approach | |
| CN105425301B (zh) | 一种频率域三维不规则地震数据重建方法 | |
| Herrmann et al. | Sparsity-and continuity-promoting seismic image recovery with curvelet frames | |
| CN111158049B (zh) | 一种基于散射积分法的地震逆时偏移成像方法 | |
| CN109765616B (zh) | 一种保幅波场延拓校正方法及系统 | |
| CN110456417B (zh) | 一种地震数据多次波压制方法 | |
| CN101614826A (zh) | 三维地震数据处理中实现面元均化的方法和装置 | |
| CN106680796A (zh) | 基于频率干涉的平面全息阵列目标三维表面重构方法 | |
| CN106199698A (zh) | 基于多次波信息的频率域地震数据重构方法 | |
| CN105510975B (zh) | 提高地震数据信噪比的方法及装置 | |
| CN112285709B (zh) | 基于深度学习的大气臭氧遥感激光雷达数据融合方法 | |
| CN113885077A (zh) | 一种基于深度学习的多震源地震数据分离方法 | |
| CN112305612B (zh) | 高分辨率复谱分解时频空间域振幅随偏移距变化校正方法 | |
| CN117055116B (zh) | 一种强稀疏约束Radon变换多次波压制方法及系统 | |
| Liu et al. | Seismic data reconstruction via complex shearlet transform and block coordinate relaxation | |
| Xu et al. | Noniterative least-squares reverse time migration based on an efficient asymptotic Hessian/point-spread function estimation | |
| Yang et al. | A seismic interpolation and denoising method with curvelet transform matching filter | |
| CN110749925A (zh) | 一种宽频逆时偏移成像处理方法 | |
| CN116719086B (zh) | 基于点扩散函数的稀疏海底四分量数据高分辨率成像方法 | |
| CN114185095B (zh) | 一种三维平面波域地震数据多次波压制的方法 | |
| Piedade et al. | Compressive sensing strategy on sparse array to accelerate ultrasonic TFM imaging | |
| Gou et al. | Complex seismic wavefield interpolation based on the Bregman iteration method in the sparse transform domain | |
| CN109917453B (zh) | 基于Shearlet变换的多尺度一次波分离方法 | |
| CN109061564B (zh) | 基于高阶累积量的简化近场定位方法 | |
| CN102998702A (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 |