Home Journals Progress in Geophysics
Progress in Geophysics

Abbreviation (ISO4): Prog Geophy      Editor in chief:

About  /  Aim & scope  /  Editorial board  /  Indexed  /  Contact  / 

Seismic data regularization and interpolation approach based on compressive sensing principle

  • LieQian DONG ,
  • Heng ZHOU ,
  • YunYun SANG ,
  • QingQin ZENG ,
  • HongGuang FAN ,
  • YongQing TIAN
Expand
  • Bureau of Geophysical Prospecting INC., China National Petroleum Corporation(BGP), Zhuozhou 072750, China

Received date: 2024-03-07

  Online published: 2025-03-13

Copyright

Copyright ©2025 Progress in Geophysics. All rights reserved.

Abstract

Seismic samples are typically designed on a perfect Cartesian grid. However, field constructions can disrupt the sampling geometry, resulting in the samples missing or off-the-grid. Our research goals are to simultaneously regularize off-the-grid samples and interpolate missing data for 3D seismic data under the framework of compressive sensing, which combines a 3D curvelet transform, a fast iterative threshold algorithm, and a merging sampling operator. The new sampling operator combines a binary mask with a barycentric Lagrangian operator for simultaneous interpolation and regularization. The fast iterative threshold algorithm is helpful to improve the interpolation accuracy and efficiency while solving the ill-posed problem. Finally, we demonstrate the effectiveness of the proposed approach by simulated and field datasets.

Cite this article

LieQian DONG , Heng ZHOU , YunYun SANG , QingQin ZENG , HongGuang FAN , YongQing TIAN . Seismic data regularization and interpolation approach based on compressive sensing principle[J]. Progress in Geophysics, 2025 , 40(1) : 276 -284 . DOI: 10.6038/pg2025HH0428

0 引言

地震勘探中,采集的地震数据都被期望在规则的网格点分布(张慕刚等,2017).但是受采集成本和采集因素(复杂的地质条件、障碍物等)的限制,实际采集的地震数据在空间方向上往往呈现不规则和不完整分布(Naghizadeh and Sacchi, 2010).对于大的障碍物会导致采集样点的缺失,对于小的障碍物,采集样点会被偏移到附近的地方进行数据的采集,引起采集数据的不规则分布,这在数据处理中需要对缺失和偏点的数据进行数据重构和规则化处理.此外,基于压缩感知的地震勘探技术(李成博和张宇,2018),在观测系统设计中采用非规则样点设计采集数据,可以降低作业成本和改善地下成像质量,也引起工业界的关注,该技术的核心也是非规则数据的重构恢复.因此,地震数据重构在地震数据处理环节具有非常重要的影响,是后期地震资料数据处理工作取得良好效果的前提.发展有效的地震数据重构方法,对实现地震数据高信噪比、高保真的重构具有重要的意义.
地震数据重构问题最早是由Larner等(1981)所提出的,他们对不完整地震道的插值及野外地震观测数据的设计进行了深入的研究.经过多年的发展,许多国内外学者就地震数据插值的问题进行了更深层次的探讨,提出了许多行之有效的重构算法.1991年,Spitz(1991)提出了一种新的插值算法,即将地震数据从t-x域上变换到f-x域上,利用变换域上数据的特点实现地震数据插值的目的.自Spitz之后,许多学者开始研究地震数据在不同变换域上的特点,以期获得新的插值算法.从t-x域算法发展到f-t域,然后到f-k域算法(宜明理等,2001).Liu和Sacchi(2004)提出最小加权范数重构算法(MWNI), 并通过选择自适应谱加权范数作为优化问题的正则化项来限制反问题的解,该方法虽然不抗假频,但是重构效果要好于常规Fourier插值算法.Xu等(2005)提出抗泄露傅里叶变换的完全不规则缺失地震数据的正则化方法,并经过不断的迭代实现了地震数据的频谱重正交,减少了频谱泄露和混叠的现象.随着压缩感知技术的兴起,Curvelet变换(Candès and Donoho, 2004)在地震数据重建和去噪领域得到了广泛地应用,且取得了良好的重建效果.刘国昌等(2011)将Curvelet变换和稀疏约束相结合,并引入L1范数进行约束,最后嵌入凸集投影(POCS)算法(Abma and Kabir, 2006)来求解反演模型.Gong等(2018)对双曲Radon变换进行改进,并将其与Curvelet变换结合,提出稀疏性框架理论,实现了缺失地震数据的重建.Hennenfent和Herrmann(2008)设计了一种抖动随机欠采样方案,以曲波变换进行数据重建,分析了该随机采集方法的性能.针对存在采样点偏离预设网格点的情况,可以采用非均一傅里叶变换(Yu et al., 2015)和线性插值的方法(Li et al., 2012)去规则化采样点偏离预设网格点的数据,然后将规则化和重构统一纳入反演框架,实现偏点数据和缺失数据的同时规则化和重构恢复.随着人工智能机器学习和深度学习的飞速发展,基于机器学习的地震数据重建法也随之兴起(俞国庆等, 2019; Wang et al., 2019; Zhang et al., 2020; 郑浩和张兵, 2020).机器学习的方法通过对地震数据自身的学习获取特征分布规律从而实现数据重建,取得了不错的重构效果.
本文主要是针对地震采集中数据缺失和由于障碍物导致的采样点偏离预设网格点的情况,基于压缩感知理论提出了一种新的数学模型,实现了地震数据的同时规则化和数据重构处理.该方法是基于一种新的融合式采样算子,三维曲波变换以及快速迭代阈值算法.融合采样算子包含了规则网格点的二值采样(mask)算子与纠正偏点的重心拉格朗日(BL)算子,mask算子实现了对缺失数据的重构,BL算子实现了偏点数据到期望输出网格点的规则化处理,两者结合实现了对地震数据的同时规则化和重构处理;快速迭代阈值算法可有效解决反问题并提高算法的计算效率.应用该技术对模型和实际地震数据进行测试,验证该方法可以有效改善地震数据品质.

1 方法原理

为了清楚的表示规则网格数据和实际采集观测的数据,令 $X, D_{\text {on }} \in \mathbb{R}^{n_t} \times \mathbb{R}^{n_s} \times \mathbb{R}^{n_r}$分别表示规则网格数据和实际观测的在规则网格点的数据. $D_{\text {of }} \in \mathbb{R}$ ${ }^{n_t} \times \mathbb{R}^m$表示实际采集的偏离规则网格点的数据. $n_{\mathrm{t}}$, $n_{\mathrm{s}}, n_{\mathrm{r}} \in \mathbb{R}$分别为时间采样点数、炮点数和接收点数. $m \in \mathbb{Z}$表示偏离规则网格点的采样点数.令 $\boldsymbol{x}=\mathrm{vec}$$(X) \in \mathbb{R}^{n_{\mathrm{t}} \times n_{\mathrm{s}} \times n_{\mathrm{r}}}, \boldsymbol{d}_{\text {on }}=\operatorname{vec}\left(D_{\text {on }}\right) \in \mathbb{R}^{n_{\mathrm{t}} \times n_{\mathrm{s}} \times n_{\mathrm{r}}}, \boldsymbol{d}_{\text {of }}=\operatorname{vec}$$\left(\boldsymbol{D}_{\text {of }}\right) \in \mathbb{R}^{n_{\mathrm{t}} \times m}$, 其中vec表示矩阵化因子,xdondof分别为XDonDof的矩阵化形式.
地震数据采集过程在数学上可以通过以下形式来表述:
$\boldsymbol{d}_{\text {on }}=\operatorname{diag}(\operatorname{vec}(M)) \boldsymbol{x} ,$
$\boldsymbol{d}_{\mathrm{of}}=\tilde{\boldsymbol{B}} \boldsymbol{x} ,$
其中,$M \in \mathbb{R}^{N_1} \times \mathbb{R}^{N_s} \times \mathbb{R}^{N_r}$ 表示采样矩阵,是一个对角线元素是0或1的对角阵,0表示数据缺失道,1表示实际采集的在规则网格点的数据. $\operatorname{diag}(\cdot)$ 表示对角阵因子.$\widetilde{\boldsymbol{B}}$ 表示BL算子.
下面介绍下$\tilde{\boldsymbol{B}}$的构建过程:
如果实际观测网格与期望输出规则网格是不相关的,可以在获得观测数据后确定期望输出规则网格.我们将实际观测网格上的真实位置表示为p1, …, pm,期望输出规则网格上的位置表示为l1, …, ln.对于j=1, …, m, $k \ll n$,给定:
$s_j \triangleq \underset{s \in\{1, \cdots, n-k\}}{\operatorname{argmin}} \prod\limits_{i=0}^k\left|p_j-l_{s+i}\right|,$
其中sj表示距离pj位置最近的规则网格点l的下标位置,即lsj是距离pj最近的规则网格点.$\operatorname{argmin} \|$ 表示取最小值,$\Pi \|$ 表示数相乘运算.则一维拉格朗日算子可以表示为:
$L_j=\frac{p_j-l_{s_j}}{l_{s_j}-l_{s_j+1}},$
式(4)是拉格朗日算子的一阶形式,也可以延拖到k阶.
在实际地震数据采集过程中,有可能存在两个方向采样点都偏离预设期望输出规则网格点的情况,因此在一维拉格朗日算子的基础上,构建了二维拉格朗日算子(Berrut and Trefethen, 2004),它由一维拉格朗日算子的重心和张量组成.假设实际采集数据和采集数据的位置分别为dij和(ai, bj),期望输出数据和期望输出的规则网格点位置分别为χ(v, z)和(vi, zj),i=1,…,n, j=1,…,m,则有:
$L_{i j}(a, b)=L_i(a) L_j(b),$
$\chi(v, z)=\sum\limits_{i=0}^n \sum\limits_{j=0}^m l_{i j}(a, b) d_{i j},$
其中Lij为二维拉格朗日算子,上述二维拉格朗日插值运算的计算复杂度为O(n2),其中n为一个方向上采样点个数,且该算子存在计算的数值不稳定的缺陷.
鉴于此,本文优化了二维拉格朗日算子,构建了二维重心拉格朗日算子.它可以降低计算复杂度,同时也增加数值计算的稳定性.重心拉格朗日权因子ωi, j的定义如下:
$\omega_{i j}=\frac{1}{\prod\limits_{k=0, k \neq i}^n\left(a_i-a_k\right) \prod\limits_{k=0, k \neq j}^n\left(b_j-b_k\right)}$
则公式(6)可以重新写成:
$\chi(v, z)=\frac{\sum\limits_{i=0}^n \sum\limits_{j=0}^m \frac{\omega_{i j}}{\left(v-x_i\right)\left(z-y_j\right)} d_{i j}}{\sum\limits_{i=0}^n \sum\limits_{j=0}^m \frac{\omega_{i j}}{\left(v-x_i\right)\left(z-y_j\right)}},$
其中:
$\tilde{B}_{i, j}=\frac{\omega_{i j}}{\left(v-x_i\right)\left(z-y_j\right)},$
$\tilde{B}_{i, j}$算子,该算子的计算复杂度为O(n),需要注意的是插值坐标(v, z)不能与给定的位置重叠,或者在实际处理中可以在分母中加上控制数$\epsilon$,避免零值点的出现.
BL算子建立了偏点数据到规则网格的映射关系.令 $\boldsymbol{d}=\left[\boldsymbol{d}_{\mathrm{on}}^{\mathrm{T}}, \boldsymbol{d}_{\mathrm{of}}^T\right]^{\mathrm{T}}, \boldsymbol{A}=\left[\operatorname{diag}(\operatorname{vec}(\boldsymbol{M})), \tilde{B}^{\mathrm{T}}\right]^{\mathrm{T}}$,则实际地震采集数据可以表示为:
$\boldsymbol{d}=\boldsymbol{A x},$
A表示融合采样算子,该算子包含了缺失数据的位置信息和偏点数据到规则网格点的信息,将两者同时纳入非规则数据重构的最优化问题中,实现了非规则数据的同时规则化和重构.
由于求解非规则数据重构的问题(9)通常是欠定的,所以数学上通常考虑将式(9)转化为带有先验约束的最优化问题然后去求解.最常用的方法之一是迭代收缩阈值算法,该最优化求解一般表达式如下所示:
$\boldsymbol{Y}_{k+1}=\left(\boldsymbol{I}-\boldsymbol{A}^*\right) C^{\mathrm{T}} T_\lambda C \boldsymbol{Y}_k+\boldsymbol{A}^* \boldsymbol{d},$
其中,Yk为第k次迭代的重构结果,Tλ为阈值函数,通常为硬阈值函数或软阈值函数,λ为选取的阈值,A *为A的逆,C表示稀疏变换,本文选取的为3D曲波变换,CT为C的逆变换.
为了进一步提高迭代收缩阈值算法计算速度,本文采用了一种优化的加速迭代阈值算法(Beck and Teboulle, 2009; 董烈乾等,2023),提高了求解最优化问题的计算效率.
$\boldsymbol{Y}_{k+1}=T_\lambda\left(\boldsymbol{Z}_k\right)+\frac{\theta_k-1}{\frac{1+\sqrt{1+4 \theta_k^2}}{2}}\left(C^{\mathrm{T}} T_\lambda C\left(\boldsymbol{Z}_k-\boldsymbol{Z}_{k-1}\right)\right),$
其中,$\boldsymbol{Z}_k=\left(\boldsymbol{I}-\boldsymbol{A}^*\right) \boldsymbol{Y}_k+\boldsymbol{A}^* \boldsymbol{d},\theta_k$是控制因子,取值为从1逐渐减小到0.

2 数据测试

为了验证方法的有效性,首先选取大小为三维模拟炮集数据进行方法验证,该正演模型数据时间采样点为1001,时间采样间隔为0.004 s,接收线为151条,每条接收线有151道接收,接收线距为100 m,道距为25 m.图 1a为规则采样的三维模拟数据,图 1b为对图 1a进行了50%非规则缺失,然后又对25%的采样点加入了半个道距的随机抖动因子,即进行了偏点处理,处理后的数据有25%的采样点在原始设计的网格点上,25%的采样点偏离原始采样网点,50%的采样点非规则缺失.图 2展示了两种采样算子重构结果对比,融合采样算子重构结果(图 2c,SNR=10.29 dB)明显好于mask采样算子重构结果(图 2a,SNR=6.41 dB),图 2bd分别为图 2ac相对于原始数据的残差.因为在mask采样算子重构方法中,采样点近似地被分配到就近的规则网格上,换句话说,mask采样方法给定的数据位置是不准确的,因此导致重建结果不理想.而融合采样算子考虑了偏点数据的位置信息,可以更好的处理非规则数据存在偏点的情况,重构效果也能够得到改善.图 3分别为模拟时间时间切片对比图.通过对比可以看出,融合采样算子重构结果重构效果更好,与原始数据差值最小.图 4为原始数据、非规则欠采样后数据、mask方法重构结果以及融合采样算子重构结果的FK谱图对比,通过对比看出融合采样算子重构结果FK谱更接近于原始数据的FK谱,主要的相干能量得到了加强,信噪比也更高,而mask采样方法重构结果相干能量存在模糊的现象,信噪比也较融合采样算子重构结果要低.
图1 三维模拟数据

(a)规则采样数据;(b)非规则欠采样数据.

Fig 1 3D simulated data

(a) Raw data; (b) Irregularly decimated data.

图2 三维模拟数据重构结果对比

(a)mask采样算子重构及(b)重构误差;(c)融合采样算子重构及(d)重构误差.

Fig 2 Reconstruction result comparisons of 3D simulated data

(a) Reconstruction result by using mask sampling operator and (b) reconstruction error; (c) Reconstruction result by using merging sampling operator and (d) reconstruction error.

图3 三维模拟数据时间切片对比图

(a)原始数据;(b)非规则欠采样;(c)mask采样算子重构及(d)重构误差;(e)融合采样算子重构及(f)重构误差.

Fig 3 Time slice comparisons of 3D simulated data

(a) Raw data; (b) Irregularly decimated data; (c) Reconstruction result by using mask sampling operator and (d) reconstruction error; (e) Reconstruction result by using merging sampling operator and (f) reconstruction error.

图4 FK谱图对比

(a)原始数据;(b)非规则欠采样后数据;(c)mask重构结果;(d)融合采样算子重构结果.

Fig 4 FK spectrum comparisons

(a) Raw data; (b) Irregularly decimated data; (c) Reconstruction result by using mask sampling operator; (d) Reconstruction result by using merging sampling operator.

接着选取实际数据进行方法测试.该采集的实际数据接收线距和炮线距分别为200 m,初始设计的规则炮点距和接收点距分别为50 m.实际采集中,炮点和接收点进行了25%的非规则缺失,由于障碍物的影响,同时存在采样点偏离预设网格点的情况.为了实现缺失采样点的重构和偏点数据的规则化,本文在十字道集进行数据的重构和规则化处理.图 5a为待处理的十字道集数据的点位展布图,图 5b图 5a红框区域的放大图.其中蓝色的实心点表示期望输出规则采样点,红色的圆圈表示缺失的采样点,蓝色的钻石形状表示实际采集的在网格点上的采样点,红色的三角形表示实际采集的偏离期望输出网格的偏点.目的是将实际采集数据中的缺失数据进行重构和偏点数据进行规则化恢复到期望输出的网格点上(即蓝色的实心点位置上).图 6a展示的是实际采集的数据,可以看出由于采样点的缺失,存在有效信息的缺失,图 6bc分别为mask算子和融合采样算子重构的结果,对比可以看出两种方法都可以将缺失的数据进行恢复,但由于实际数据存在偏点的影响,mask重构方法没有包含偏点的位置信息,并不能很好地将偏点的数据进行规则化处理,存在同相轴的扭曲,而融合采样算子重构包含了偏点到规则网格点的信息,能够很好地处理存在偏点数据的重构,重构后同相轴也更加连续(红色箭头所指),整体信噪比也得到提升.图 7为实际数据重构前后时间切片的对比,对比可以清楚的看出,缺失的数据得到了恢复,与mask重构方法对比,本文基于融合采样算子的重构方法在解决存在偏点的地震数据数据重构上更优优势,重构效果更好.
图5 (a) 实际数据采样点分布图;(b)图(a)红框区域采样点分布放大图

Fig 5 (a) Sampling point distribution map of field data; (b) Enlarged view of rectangle area in Fig.(a)

图6 实际数据重构结果对比

(a)原始采集的数据;(b)mask采样算子重构结果;(c)融合采样算子重构结果.

Fig 6 Reconstruction result comparisons of field data

(a) Raw data; (b) Reconstruction result by using mask sampling operator; (c) Reconstruction result by using merging sampling operator.

图7 实际数据时间切片效果对比

(a)原始采集的数据;(b)mask采样算子重构结果;(c)融合采样算子重构结果.

Fig 7 Time slice comparisons of field data

(a) Raw data; (b) Reconstruction result by using mask sampling operator; (c) Reconstruction result by using merging sampling operator.

3 结论

本文在压缩感知理论基础上,提出了一种新的地震数据重构模型,主要亮点在于采用了一种新的融合采样算子.该融合采样算子结合了规则网格点的二值采样算子与纠正偏点的重心拉格朗日算子,可以将缺失数据和偏点数据的位置信息包含在新的采样算子中,然后利用基于压缩感知的重构技术实现对地震数据的同时重构和规则化处理.应用该方法对模拟数据和实际数据进行重构测试,并与基于规则网格点的mask采样算子重构方法对比,本文方法可以在恢复缺失的数据的同时,能够改善由于偏点对地震数据重构的影响,重构后数据同相轴连续性更好,数据整体信噪比也更高.

感谢两位匿名审稿专家给予的宝贵修改建议,感谢哈尔滨工业大学于四伟教授对本文方法实现给与的建设性建议和帮助.

Abma R , Kabir N . 3D interpolation of irregular data with a POCS algorithm. Geophysics, 2006, 71 (6): E91- E97.

DOI

Beck A , Teboulle M . A fast iterative shrinkage-thresholding algorithm for linear inverse problems. SIAM Journal on Imaging Sciences, 2009, 2 (1): 183- 202.

DOI

Berrut J P , Trefethen L N . Barycentric lagrange interpolation. SIAM Review, 2004, 46 (3): 501- 517.

DOI

Candès E J , Donoho D L . New tight frames of curvelets and optimal representations of objects with piecewise C2 singularities. Communion Pure and Applied Mathematics, 2004, 57 (2): 219- 266.

DOI

Dong L Q , Zhang M G , Wang C H , et al. Reconstruction of non-uniformly sampled seismic data based on fast POCS algorithm. Oil Geophysical Prospecting, 2023, 58 (2): 334- 339.

Gong X B , Wang S C , Du L Z . Seismic data reconstruction using a sparsity-promoting apex shifted hyperbolic Radon-curvelet transform. Studia Geophysica et Geodaetica, 2018, 62 (3): 450- 465.

DOI

Hennenfent G , Herrmann F J . Simply denoise: wavefield reconstruction via jittered undersampling. Geophysics, 2008, 73 (3): V19- V28.

DOI

Larner K , Gibson B , Rothman D . Trace interpolation and the design of seismic surveys. Geophysics, 1981, 46 (4): 407- 415.

Li C B, Mosher C C, Kaplan S T. 2012. Interpolated compressive sensing for seismic data reconstruction. //SEG Technical Program Expanded Abstracts 2012. SEG, 1 -6.

Li C B , Zhang Y . CSI: an efficient high-resolution seismic acquisition technology based on compressive sensing. Geophysical Prospecting for Petroleum, 2018, 57 (4): 537- 542.

Liu B , Sacchi M D . Minimum weighted norm interpolation of seismic records. Geophysics, 2004, 69 (6): 1560- 1568.

DOI

Liu G C , Chen X H , Guo Z F , et al. Missing seismic data rebuilding by interpolation based on Curvelet transform. Oil Geophysical Prospecting, 2011, 46 (2): 237- 246.

Naghizadeh M , Sacchi M D . Beyond alias hierarchical scale curvelet interpolation of regularly and irregularly sampled seismic data. Geophysics, 2010, 75 (6): WB189- WB202.

DOI

Spitz S . Seismic trace interpolation in the F-X domain. Geophysics, 1991, 56 (6): 785- 794.

DOI

Wang B , Zhang N , Lu W , et al. Deep-learning-based seismic data interpolation: a preliminary result. Geophysics, 2019, 84 (1): V11- V20.

DOI

Xu S , Zhang Y , Pham D , et al. Antileakage fourier transform for seismic data regularization. Geophysics, 2005, 70 (4): V87- V95.

DOI

Yi M L , Yan Y S , Wei X , et al. Using anti-alias trace interpolation in F-K domain. Geophysical Prospecting for Petroleum, 2001, 40 (2): 36- 41.

Yu G Q , Jia R S , Sun Y Y , et al. Reconstruction of seismic data based on over-complete dictionary learning. Progress in Geophysics, 2019, 34 (1): 229- 235.

Yu S W , Ma J W , Zhang X Q , et al. Interpolation and denoising of high-dimensional seismic data by learning a tight frame. Geophysics, 2015, 80 (5): V119- V132.

DOI

Zhang H , Yang X Y , Ma J W . Can learning from natural image denoising be used for seismic data interpolation?. Geophysics, 2020, WA115- WA136.

Zhang M G , Luo F , Wang C H , et al. Industrialization progress of wide-azimuth, broadband, and high-density (WBH) seismic acquisition technique. Natural Gas Industry, 2017, 37 (11): 1- 8.

Zheng H , Zhang B . Intelligent seismic data interpolation via convolutional neural network. Progress in Geophysics, 2020, 35 (2): 721- 727.

烈乾 , 慕刚 , 长辉 , 等. 应用快速POCS算法的非均匀地震数据重建. 石油地球物理勘探, 2023, 58 (2): 334- 339.

成博 , . CSI: 基于压缩感知的高精度高效率地震资料采集技术. 石油物探, 2018, 57 (4): 537- 542.

国昌 , 小宏 , 志峰 , 等. 基于Curvelet变换的缺失地震数据插值方法. 石油地球物理勘探, 2011, 46 (2): 237- 246.

明理 , 又生 , , 等. F-K域抗假频道内插. 石油物探, 2001, 40 (2): 36- 41.

国庆 , 瑞生 , 圆圆 , 等. 基于超完备字典学习的缺失地震数据重构方法. 地球物理学进展, 2019, 34 (1): 229- 235.

DOI

慕刚 , , 长辉 , 等. "两宽一高"地震采集技术工业化应用的进展. 天然气工业, 2017, 37 (11): 1- 8.

, . 基于卷积神经网络的智能化地震数据插值技术. 地球物理学进展, 2020, 35 (2): 721- 727.

DOI

Outlines

/