地球物理学报  2021, Vol. 64 Issue (10): 3730-3741   PDF    
基于频率-结构双重加权阈值的地震数据插值
陈思远1, 曹思远1, 孙耀光1, 江雨濛1, 高君2     
1. 中国石油大学(北京), 北京 102249;
2. 中国石化石油勘探开发研究院, 北京 100083
摘要:地震数据的插值是地震数据处理的关键环节,插值效果的优劣关系到后续反演、去噪、反褶积、偏移等工作的进行.基于稀疏反演的常规频率-波数(FK)域插值方法往往使用广义阈值,而由于缺失地震数据仅沿波数方向呈稀疏态,在沿频率方向呈类子波频谱形态分布的非稀疏态,直接应用基于Lp范数的广义阈值会难以补全缺失数据的高低频信息,不利于后续依赖低频信息的反演及依赖高频的高分辨率处理工作的开展;基于该问题,本文设计沿频率方向类子波频谱形态变化的频率加权阈值,同时考虑完整地震数据在波数域的聚集分布,引入沿波数方向的结构加权阈值,提出基于频率-结构的双重加权阈值,在补全缺失地震道的同时,保证其高低频信息的有效恢复.模型和实际数据表明,相比于广义阈值,基于频率-结构的双重加权阈值可以减少插值误差,最大限度保证缺失地震道位置插值结果高低频信息的完整性.
关键词: 插值      结构加权      频率加权      FK变换      稀疏反演     
Seismic data interpolation based on the frequency-structure dual weighted threshold
CHEN SiYuan1, CAO SiYuan1, SUN YaoGuang1, JIANG YuMeng1, GAO Jun2     
1. China University of Petroleum(Beijing), Beijing 102249, China;
2. Sinopec Petroleum Exploration andProduction Research Institute, Beijing 100083, China
Abstract: Seismic data interpolation is a crucial part of seismic data processing. The quality of the interpolation is related to the subsequent inversion, denoising, deconvolution, migration, and other work. Conventional FK interpolation methods based on sparse inversion often use generalized thresholds. However, the missing seismic data in these methods only shows a sparse state along the wavenumber direction and a non-sparse state with a wavelet-like spectrum distribution along the frequency direction. The direct application of the generalized threshold based on the Lp norm will make it difficult to complete the high- and low-frequency information of the missing data, which is not conducive to the subsequent inversion that relies on low-frequency information and the development of high-resolution processing that relies on high frequency. To solve this problem, we design a frequency-weighted threshold that changes along the frequency direction and introduces a structure-weighted threshold along the wavenumber direction considering the aggregation distribution of complete seismic data in the wavenumber domain. Combining the above two points, we propose a frequency-structure dual weighted threshold to ensure the effective recovery of high- and low-frequency information while completing missing traces. By comparing with the generalized threshold, the model and actual data show that the frequency-structure dual weighted threshold can reduce the interpolation error and ensure the integrity of the high- and low-frequency information to the largest extent.
Keywords: Interpolation    Structure weighted    Frequency weighted    FK transform    Sparse inversion    
0 引言

地震采集与处理过程中,受空采区和坏道的影响,会出现地震道的缺失,需进行数据补全,为后续地震数据的反褶积、偏移以及储层反演等工作提供优质数据体.最早使用的插值方法为FX域(Spitz, 1991)和F-XY域(Wang, 2002)的插值方法;考虑相邻地震道的相似性,出现以差分算子为约束的插值方法(Liu et al., 2015; Zu et al., 2016),并认为缺失地震数据左右两侧的平均为缺失处的数据,由于叠前数据倾角较大,该方法效果不尽理想;由此,基于倾角约束(Hedefa et al., 2010)的插值方法便应运而生,包括基于倾角约束的中值滤波方法等(Huo et al., 2017).另一方面,在波动方程领域,有部分学者以速度为约束,以波动方程为媒介进行插值等(Swindeman and Fomel, 2015).

近年,随着压缩感知恢复算法的发展,逐渐提出基于稀疏约束的插值算法(刘洋等,2018孙苗苗等,2019),其假设无缺失数据经过某种变换后,在变换域中呈稀疏态,而缺失地震数据在变换域内非稀疏,相比于无缺失数据,缺失地震数据变换域内的幅值较弱,可使用阈值压制这些非稀疏能量以完成缺失数据的补全.目前,可用变换包括Fourier变换(Gülünay, 2003王锦妍等,2020)、FK变换(Wang and Lu, 2017)、小波变换(Guo et al., 2015)、Curvelet变换(Herrmann and Hennenfent, 2008; Yang et al., 2012)、Shearlet变换(Liu et al., 2018)、Radon变换(Yu et al., 2007)及Seislet变换(Gan et al., 2015邓盾等,2019)等.除此之外,低秩约束作为特殊的稀疏约束也被应用到地震数据的插值中(Chen et al., 2019; Lari et al., 2019; Pan et al., 2017).

在稀疏约束目标上,L0范数是标准的稀疏约束,该范数的求解是NP难问题(Natarajan, 1995),可使用L0范数的最佳凸近似L1范数代替L0范数求解(Chen et al., 1998),当s表示稀疏度,且满足δ2s<sqrt(2)-1(Candès et al., 2006)时两者解相同,基于整形正则化(Beck and Teboulle, 2009; Fomel, 2007; Goldstein and Stanley, 2009)等方法可有效求解;近年,相比于L1范数,非凸拟范数Lp(0<p<1) 的出现更好近似了L0范数,该范数稀疏约束能力较高的原因是其本质上修改了整形正则化中的阈值,增大阈值之间的差距,进而增强变换域的稀疏约束能力,基于该拟范数的稀疏恢复问题,可在多项式时间内有效求解(Ge et al., 2011).上述稀疏约束假设数据在变换域中满足严格k稀疏,所使用的阈值作用于整个变换域,但实际数据不满足该条件;同时地震数据波数、频率、时间等在变换域中具有非随机分布特征,即稀疏变换域的结构特性表现为非均匀分布,使得变换域的数据表现为近似k稀疏,使上述算法在稀疏重建过程中存在重建误差,所损失部分为变换域中有效波的弱能量(Foucart and Rauhut, 2017).

基于稀疏约束的地震数据插值方向,深入讨论稀疏变换域结构特性的研究论文较少,通常假设其整体的稀疏性,未考虑信号及变换域的特征,导致变换域中的弱能量被广义阈值压制.本文从频率-波数(FK)域的结构入手,探究沿频率方向和沿波数方向有效能量及缺失数据分布特征,以减少变换域中的弱有效能量损失.在现有的研究基础上开展如下工作:(1)以FK变换为例,使用图示描述基于Lp范数的广义阈值(下称广义阈值)所造成的被插值数据的高低频缺失现象,探究在FK变换域中缺失数据沿频率方向的形态,构造频率加权阈值(FWT);(2)针对FK变换域的波数方向的聚集性特征,构造结构加权阈值(SWT)以提高稀疏约束能力;(3)模型测试部分,合成简单模型分别测试“广义阈值”、“频率加权阈值(FWT)”、“结构加权阈值(SWT)”、“频率-结构双重加权阈值(FSWT)”在插值中的作用,以佐证理论部分的正确性;最后,以实际缺失地震数据测试本文所提出的频率-结构双重加权阈值在缺失数据高低频恢复上的有效性.

需要指出,除FK变换域可作为稀疏变换域之外,Curvelet变换、小波变换、Shearlet变换等均带有频率的分解,以满足不同频率下数据的稀疏表达;因为地震数据各频率成分比重不同,所以本研究以FK变换作为稀疏变换旨在证明“频率-结构双重加权阈值”相比于“广义阈值”在数据高低频恢复中的有效性,其他的带有频率分解的类似稀疏变换均可依照该思想构造相应的加权阈值.

1 理论 1.1 基于稀疏反演的地震数据插值理论

时间采样点为nm道地震数据插值可以使用正演方程描述为:

(1)

其中,Ymn×1表示观测数据,即含缺失地震道的数据,L定义为式(2)所示的正演算子,其中,I0分别表示单位矩阵和零矩阵,可从完整数据X中采样得到观测数据.式(2)为:

(2)

采样矩阵L是秩亏矩阵,待求的完整数据X具有无穷个解,基于压缩感知理论,可施加不同的正则化约束保证待求X的解唯一,并尽可能获得最接近真实数据的解.常用的正则化约束包括Tikhonov约束(Fleming, 1990)(L2范数约束),稀疏约束(L0范数约束)等,可通过压缩感知理论建立如下约束方程:

(3)

其中,B表示某种变换,完整数据X在变换B的作用下可以符合Lq-Norm(0≤q<2)的约束,整形正则化理论表明,可以通过约束待求解数据变换域中的形态以获得最接近真实数据的解;传统方法假设完整数据在FK域中稀疏,而缺失后的数据非稀疏,通过整形正则化的方法求解该模型得到最接近真实数据的解.将式(3)写为如式(4)所示的无约束优化问题,并使用稀疏约束作为正则化项:

(4)

式中,λ定义为正则化算子,由于L0-norm的求解是NP-hard问题,所使用的求解算法是匹配追踪等(Mallat and Zhang, 1993),该类算法需预先给定稀疏度,而完整数据FK域的数据稀疏度无法确定,故匹配追踪类算法不适用方程(4)的求解;考虑以Lp-norm近似L0-norm,(0<p≤1);p值越小,所约束的BX越稀疏,改写式(4)为:

(5)

式(5)可以通过整形正则化的方法在多项式时间内有效求解(Ge et al., 2011),给出迭代步长α,对应于式(5)的迭代方程为:

(6)

式中,Tλ为式(7)定义的广义软阈值处理算子:

(7)

该算子返回广义阈值λ2-p|X|p-1处理后的结果,其中0<p≤1,|X|越大,阈值λ2-p|X|p-1越小;反之,|X|越小,阈值λ2-p|X|p-1越大;当p=1时,广义阈值退化为常数阈值.

1.2 频率-结构双重加权阈值的构造 1.2.1 含缺失道的地震数据在FK域的结构特性

B为FK变换时,式(5)被称为FK域的插值反演方程.地震数据因为其频带宽度有限,同相轴曲率较小,在FK域中呈聚集性分布;高波数、高频率位置因地震数据构造的特殊性,也会出现部分弱能量.而由于地震数据不规则缺失现象,缺失地震道的位置产生截断,如图 1所示,造成沿波数方向的“频谱泄漏”,随缺失道的增加,“频谱泄漏”的能量也随之增加,缺失地震道的能量将类似于“白噪声”均匀分布在沿波数方向上(图 1中F1处的虚线),可使用广义阈值压制(图 1下方虚线);另外,在沿频率方向上,缺失数据并未造成频谱泄漏(图 1中K1处的虚线),故而广义阈值应单独作用于沿波数方向,不应作用于沿频率方向上.

图 1 缺失地震数据FK谱形态示意图 Fig. 1 Sketch showing FK spectrum of missing data

图 2模拟广义阈值对白噪声(频谱为白谱)和类子波频谱形态的噪声(频谱形态类似为子波频谱的形态)的压制效果及所带来的问题;当噪声为白噪声时(图 2a1),合适的广义阈值可以有效的压制噪声的干扰(图 2d1图 2d2),且可以保留高低频能量;而当噪声为子波频谱形态时(图 2a3),广义阈值也压制了噪声的干扰,但是损失了高频和低频的能量(图 2d3图 2d4).

图 2 不同频谱形态的噪声与常数阈值的测试示意图 Fig. 2 Schematic diagram showing noise of different spectrum forms and the test of constant threshold
1.2.2 基于FK域的频率-结构双重加权阈值的构造

在地震数据处理中,低频和高频的能量是提高分辨率和储层反演所需要的,广义阈值在处理类子波频谱形态的噪声中损伤了能量较弱的低频及高频信号;对于类子波频谱形态的噪声,考虑使用类子波频谱形态的阈值进行处理,即能量弱的部分减小阈值,联系上文缺失地震数据FK变换,沿波数方向的“频谱泄漏”相当于在每个波数上均增加了类子波频谱形态的噪声,在使用广义阈值进行处理时,低频和高频所对应频谱能量较弱的位置,广义阈值会压制弱能量,将造成缺失数据高低频的损失,故而在基于FK域的地震数据插值中需要构造呈类子波频谱形态的沿频率方向的阈值加权算子.

缺失数据造成的“频谱泄漏”沿波数方向呈均匀随机分布,可使用广义阈值进行压制,当缺失数目较多时,“频谱泄漏”的能量逐渐逼近甚至超过有效信号的能量,这需要沿波数方向构造一个随FK谱结构变化的加权算子,进行沿波数方向的结构加权;基于频率-结构的双重阈值加权算子的构造方式如下.

记观测数据Y经过沿时间方向的Fourier变换得到的数据为Yf, x,式(1)可写为:

(8)

式中,Xf, x表示完成数据X经过沿时间方向的Fourier变换得到的数据,对于数据的每个频率切片fi, (i=1, 2, …, m), 式(8)可写为:

(9)

其中,L0=diag{1, 1, …, 0, …, 1}n×n仍为采样算子,若Xfi, x在变换域B中为稀疏态,且Yfi, x在变换域B中为非稀疏,则可通过式(10)所示的压缩感知重建方法得到每个频率切片对应的完整的数据Xfi, x遍历所有频率切片即完成缺失数据的恢复:

(10)

根据式(3)—(6),约束方程(10)也可以使用整形正则化方法求解,即:

(11)

(12)

若保证迭代式(11)的有效求解,正则化系数λ∝‖BYfi, xB为沿空间方向的Fourier变换时,作为广义阈值λ2-p|Xfi, x|p-1,其正比于每个频率切片沿波数方向“频谱泄漏”振幅的均值.基于此,对每个频率切片进行插值恢复时,需要赋以不同的正则化系数,且该正则化系数与每个频率切片沿波数方向“频谱泄漏”振幅的均值成正比,即呈类子波频谱形态的沿频率方向的频率加权阈值.

另外,考虑变换域的稀疏性,引入反加权算子hfi, x增强完整数据的每个频率切片的稀疏约束能力,则式(10)改写为:

(13)

反加权算子hfi, x作用于正则化系数λ上,其对变换域的数据进行加权处理,减小|BXfi, x|能量较强部分的阈值,增强|BXfi, x|能量较弱部分的阈值,从而增强对变换域的稀疏约束,Candès等(2008)已证明式(13)所对应求解方法的收敛性,可有效提高稀疏约束.由于反加权算子hfi, x数据驱动,具备自适应性,考虑B为沿空间方向的Fourier变换,对应着数据X的结构特征,故而称(λhfi, k)2-p|Xfi, K|p-1为沿波数方向的结构加权阈值.

综上所述,假设缺失的地震记录X(f, k)经过FK变换为M(f, k),其中,f表示频率,k为波数,首先沿波数k方向取平均值W(f),对W(f)沿f方向平滑,得到类子波谱(f),沿波数方向复制该谱,得到(f, k),联合广义阈值λ2-p|X|p-1,引入⊗表示Hadamard product,得到频率加权阈值(λ)2-p⊗|X|p-1(FWT). 然后设定阈值θ,并令其中μ<1,构造FK谱的反加权算子H(f, k),得到结构加权阈值(λH)2-p⊗|X|p-1(SWT),基于FK变换的频率和结构的双重加权阈值(FSWT)可写为(λF)2-p⊗|X|p-1=[(λ)2-p⊗|X|p-1]⊗[(λH)2-p⊗|X|p-1],当F ≡1时,双重加权阈值退化为广义阈值.

将频率-结构双重加权阈值代入式(7)中,将广义软阈值修改为式(14),结合整形正则化的迭代式(6),可进行基于频率-结构双重加权阈值的缺失地震数据补全:

(14)

2 模型测试

本部分首先使用128×128的合成模型,分别测试广义阈值,频率加权阈值(FWT),结构加权阈值(SWT)和频率-结构双重加权阈值(FSWT)各自的作用;如图 3a所示,该模型从左到右Ricker子波主频由15 Hz递增到95 Hz;加入高斯白噪后信噪比为24.68 dB.

图 3 原始数据及插值误差 (a) 原始数据;(b) 基于广义阈值的插值误差;(c) 基于FWT的插值误差;(d) 基于SWT的插值误差;(e) 基于FSWT的插值误差. Fig. 3 Original data and interpolation errors (a) Original data; (b) Interpolation error using generalized threshold; (c) Interpolation error using FWT; (d) Interpolation error using SWT; (e) Interpolation error using FSWT.

应用上述四种阈值进行插值,单次插值误差剖面见图 3be,因为SWT作为结构加权阈值,是全频带的保幅阈值,误差变化不明显;同时,FWT重点倾向于弱频率能量的保护,在此模型中对应于能量较弱的高频,故而在全频带的误差剖面上,四种阈值所产生的误差形态接近.经计算,四种阈值插值结果与原剖面的余弦相似度分别为:广义阈值:99.27%,FWT:99.44%,SWT:99.41%,FSWT:99.56%.

为了直观的体现频率加权阈值对数据高频补全的重要性,应用时频分析工具对四种阈值的插值结果及无缺失数据进行频率分解,分解得到的高频剖面见图 4.广义阈值和SWT的插值结果的高频剖面表明基于广义阈值的插值难以补全高频成分(图 4b),而本文提出的FWT和FSWT的插值结果分频剖面同相轴连续性较好(图 4c图 4e),证明FWT和FSWT在缺失数据高频恢复上的有效性.

图 4 分频剖面(高频) (a) 原始数据分频剖面;(b) 基于广义阈值的插值分频剖面;(c) 基于FWT的插值分频剖面;(d) 基于SWT的插值分频剖面;(e) 基于FSWT的插值分频剖面. Fig. 4 Frequency-division profiles (high frequency) (a) Frequency-division profile of original data; (b) Interpolated frequency-division profile using generalized threshold; (c) Interpolated frequency-division profile using FWT; (d) Interpolated frequency-division profile using SWT; (e) Interpolated frequency-division profile using FSWT.

因为模型的主频是道号的函数,本部分重复试验100次随机缺失数据插值,绘制四种阈值插值结果逐道与原剖面(图 3a)的相似度(图 5),当频率处于中低频时(<35 Hz),相似性FSWT>SWT>FWT=广义阈值;在证明SWT的保幅效果优异的同时,也说明FWT在能量较强的中低频部分保幅效果较弱;在高频部分(>70 Hz),相似性FSWT>FWT>SWT>广义阈值,说明FWT在数据高频部分保幅性优异;而SWT在整个频带的相似性排名均优于广义阈值,证明SWT在全频带具有一定保幅作用.

图 5 四种阈值插值结果与无缺失数据的余弦相似度 Fig. 5 Cosine similarity between the original data and interpolated data using four thresholds

继续合成256道,含256个采样点的叠前地震数据;如图 6所示为原始数据、缺失50%地震道的地震数据和缺失地震数据的FK谱;添加随机噪声后,信噪比为32 dB.本部分将测试广义阈值和频率-结构双重加权阈值(FSWT)对图 6所示复杂模型的适用性.测试环境为Inter i7-9750H处理器,16 GB内存,软件环境为Matlab2020a.

图 6 合成模型 (a) 原始数据;(b) 缺失数据;(c) 图 6a的FK谱. Fig. 6 Synthesis model (a) Original data; (b) Missing data; (c) FK spectrum of Fig. 6a.

对于插值反演方程,正则化参数对插值结果影响相对较大.分别选择不同的正则化参数进行插值测试,记录插值后与原始数据的余弦相似度、插值后与原始数据低频剖面以及高频剖面的余弦相似度和计算时间的曲线如图 7所示:当正则化参数逐渐变小时,任何频带范围内,基于广义阈值的插值方法和基于FSWT的插值方法所计算的余弦相似度趋于相同(图 7ac),但是计算量呈指数型增加(图 7d);而当正则化参数逐渐变大时,基于广义阈值的插值方法的相似度快速下降,而基于FSWT的插值方法余弦相似度变化不明显(图 7ac),同时,计算时间趋于稳定(图 7d),综合考虑计算成本,如图 8所示取正则化参数等于8时的插值结果进行分析.

图 7 正则化参数对插值性能的影响测试 (a) 全频带余弦相似度;(b) 低频余弦相似度;(c) 高频余弦相似度;(d) CPU计算时间. Fig. 7 Interpolation performance of regularization parameters (a) Cosine similarity in full frequency band; (b) Cosine similarity in the low frequency range; (c) Cosine similarity in high-frequency band; (d) CPU calculation time.
图 8 基于广义阈值和FSWT的插值结果 (a) 基于广义阈值的插值剖面;(b) 图 6a和图 8a的差;(c) 图 8b的FK谱;(d) 基于FSWT的插值剖面;(e) 图 6a和图 8e的差;(f) 图 8e的FK谱. Fig. 8 Interpolated results using generalized threshold and FSWT (a) Interpolated profile using generalized threshold; (b) Difference between Fig. 6a and Fig. 8a; (c) FK spectrum of Fig. 8b; (d) Interpolated profile using FSWT; (e) Difference between Fig. 6a and Fig. 8e; (f) FK spectrum of Fig. 8e.

基于广义阈值的插值方法的误差剖面上(图 8b)可见大量同相轴,且误差剖面的FK谱上有较多聚集性能量(图 8c),插值结果较差;相比于前者,基于FSWT的插值方法的误差剖面无明显可见同相轴,FK域中的聚集能量也较少(图 8d).

图 6a所示的无缺失数据进行分频,取低频(主频约18 Hz)及高频(主频约66 Hz)的剖面如图 9所示;同时将图 8a图 8d所示的插值结果也进行分频,绘制高低频插值结果及误差剖面如图 10所示;无论是低频或者高频部分,基于广义阈值的插值方法的分频误差剖面上(图 10b图 10d)均可见同相轴,且缺失数据未成功恢复;基于FSWT的插值方法低频及高频误差均较小,无缺失地震道现象出现(图 10f图 10h).

图 9 无缺失数据分频剖面 (a) 低频剖面;(b) 高频剖面. Fig. 9 The frequency-division profile of complete data (a) Profile in the low-frequency band; (b) Profile in the high-frequency band.
图 10 基于广义阈值和FSWT的插值结果分频剖面 (a) 基于广义阈值的插值结果低频剖面;(b) 图 9a和图 10a的差;(c) 基于广义阈值的插值结果高频剖面;(d) 图 9b和图 10c的差;(e) 基于FSWT的插值结果低频剖面;(f) 图 9a和图 10e的差;(g) 基于广义阈值的插值结果高频剖面;(h) 图 9b和图 10g的差. Fig. 10 Interpolated frequency-division profiles using generalized threshold and FSWT (a) Interpolated profile in the low frequency range using generalized threshold; (b) Difference between Fig. 9a and Fig. 10a; (c) Interpolated profile in the high frequency range using generalized threshold; (d) Difference between Fig. 9b and Fig. 10c; (e) Interpolated profile in the low frequency range using FSWT; (f) Difference between Fig. 9a and Fig. 10e; (g) Interpolated profile in the high frequency range using FSWT; (h) Difference between Fig. 9b and Fig. 10g.
3 实际数据测试

实际数据测试部分,选择M地区实际数据进行插值处理,该数据含有180道,采样率为4 ms,本文置零约35%的地震道以测试广义阈值和FSWT的插值效果.如图所示,图 11a图 11b分别为原始数据和缺失地震道的数据;图 12为使用广义阈值和FSWT的插值结果与误差,从插值剖面上无法看出两种阈值的区别(图 12a图 12c),在相应的误差剖面上(图 12b图 12d)FSWT的插值误差小于使用广义阈值的插值误差;同时,本文计算使用两种阈值插值后剖面缺失位置所占无缺失数据缺失位置能量的比重分别为:广义阈值:75.7%;频率-结构双重加权阈值:89.9%,这证明本文所提出的FSWT在全频带上的保幅性优越于传统广义阈值.

图 11 实际单炮数据和缺失35%地震道的数据 (a) 原始无缺失数据;(b) 缺失35%地震道的数据. Fig. 11 Real single shot and missing seismic data (a) Original complete data; (b) Missing seismic data with 35% missing traces.
图 12 基于广义阈值和FSWT的插值剖面和误差 (a) 基于广义阈值的插值剖面;(b) 图 12a和图 11a的差;(c) 基于FSWT的插值剖面;(d) 图 12c和图 11a的差. Fig. 12 Interpolated profiles and errors using generalized threshold and FSWT (a) Interpolated profile using generalized threshold; (b) Difference between Fig. 12a and Fig. 11a; (c) Interpolated profile using FSWT; (d) Difference between Fig. 12c and Fig. 11a.

图 11a图 12a图 12c的缺失地震数据位置的频谱如图 13所示,基于FSWT插值的结果较好的拟合原始数据的频谱,而基于广义阈值的插值结果在低频和高频与原始数据频谱的相似度均较弱;为了得到更准确的测试结果,应用时频分析工具对原始数据(图 11a)以及插值后的数据图 12a图 12c进行频率分解,低频(主频约为9 Hz)及高频剖面(主频约为90 Hz)如图 14图 15所示,相比于图 14所示的原始数据分频结果,本文提出的FSWT无论是低频(图 15b)或者高频(图 15d)的插值效果均优异于广义阈值的插值效果(图 15a图 15c).

图 13 实际单炮数据和插值后数据的频谱 Fig. 13 The spectrum of real single shot and interpolated data
图 14 无缺失数据的高低频剖面 (a) 无缺失数据的低频剖面;(b) 无缺失数据的高频剖面. Fig. 14 Profiles of complete data in high- and low-frequency range (a) Profile of complete data in low-frequency range; (b) Profile of complete data in high-frequency range.
图 15 基于广义阈值和FSWT插值后的高低频剖面 (a) 基于广义阈值插值后的低频剖面;(b) 基于FSWT插值后的低频剖面;(c) 基于广义阈值插值后的高频剖面;(d) 基于FSWT插值后的高频剖面. Fig. 15 Interpolated profiles in the high- and low-frequency range using generalized threshold and FSWT (a) Interpolated profile in the low-frequency range using generalized threshold; (b) Interpolated profiles in the low-frequency range using FSWT; (c) Interpolated profile in the high-frequency range using generalized threshold; (d) Interpolated profile in the high-frequency range using FSWT.
4 结论

针对基于传统广义阈值插值所造成高低频难以补全的问题,本文考虑FK域中缺失数据频率方向的非稀疏,构造类子波形状的频率加权算子,同时注意到波数域有效能量聚集性,联合频率加权算子,提出频率-结构双重加权算子进行缺失地震数据的补全,一定程度上解决了传统广义阈值对于缺失位置高低频难以恢复的问题,保证后续处理环节的进行.算法在适用性方面仍有较大研究空间:

(1) 基于频率-结构双重加权阈值主要作用在于恢复缺失位置处、能量较弱高低频信息,由于中频段(主频附近)的FSWT与广义阈值的大小、形态均接近,所以在该频段内算法改善效果不佳.

(2) 频率域的非离散现象不仅在FK域中出现,包括Curvelet域、Shearlet域等涉及频率的变换域均有此现象,也可通过构造频率加权算子提高这些变换域的插值精度.

(3) 对于三维数据,可采用频率-波数-波数(FKK)或者3D-Curvelet作为稀疏变换域,按本文方法可构造三维加权阈值进行高低频信息的有效恢复.

(4) 本文仅在FK域进行加权算子的可行性研究,由于地质构造通常较复杂,可能含有绕射波等大倾角构造,由于变换域限制,大倾角同相轴补全效果不佳,可使用其他稀疏变换域(Curvelet变换)联合频率-结构双重加权阈值对这部分弱能量进行有效恢复.

References
Beck A, Teboulle M. 2009. A fast iterative shrinkage-thresholding algorithm for linear inverse problems. Siam Journal on Imaging Sciences, 2(1): 183-202. DOI:10.1137/080716542
Candès E J, Romberg J K, Tao T T. 2006. Stable signal recovery from incomplete and inaccurate measurements. Communications on Pure and Applied Mathematics, 59(8): 1207-1223. DOI:10.1002/cpa.20124
Candès E J, Wakin M B, Boyd S P. 2008. Enhancing Sparsity by Reweighted 1 Minimization. Journal of Fourier Analysis and Applications, 14(5-6): 877-905. DOI:10.1007/s00041-008-9045-x
Chen SS, Donoho D L, Saunders M A. 1998. Atomic decomposition by basis pursuit. SIAM Journal on Scientific Computing, 20(1): 33-61. DOI:10.1137/S1064827596304010
Chen Y K, Chen X H, Wang Y F, et al. 2019. The interpolation of sparse geophysical data. Surveys in Geophysics, 40(1): 73-105. DOI:10.1007/s10712-018-9501-3
Deng D, Pei J X, Deng Y, et al. 2019. Key techniques research for improving mid-deep layer structure imaging in deep water seismic data of South China Sea. Progress in Geophysics (in Chinese), 34(2): 702-708. DOI:10.6038/pg2019BB0551
Fleming H E. 1990. Equivalence of regularization and truncated iteration in the solution of Ⅲ-posed image reconstruction problems. Linear Algebra and Its Applications, 130: 133-150. DOI:10.1016/0024-3795(90)90210-4
Fomel S. 2007. Shaping regularization in geophysical-estimation problems. Geophysics, 72(2): R29-R36. DOI:10.1190/1.2433716
Foucart S, Rauhut H. 2017. A mathematical introduction to compressive sensing. Bull. Am. Math, 54: 151-165.
Gan S W, Wang S D, Chen Y K, et al. 2015. Dealiased seismic data interpolation using seislet transform with low-frequency constraint. IEEE Geoscience and Remote Sensing Letters, 12(10): 2150-2154. DOI:10.1109/LGRS.2015.2453119
Ge D D, Jiang X Y, Ye Y Y. 2011. A note on the complexity of Lp minimization. Mathematical Programming, 129(2): 285-299. DOI:10.1007/s10107-011-0470-2
Goldstein T, Osher S. 2009. The split Bregman method for L1-regularized problems. Siam Journal on Imaging Sciences, 2(2): 323-343. DOI:10.1137/080725891
Guo N M, Chen M, Cui Y F. 2015. Prestack seismic data reconstruction method based undecimated discrete wavelet transform. //85th Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts, 1-5.
Gülünay N. 2003. Seismic trace interpolation in the Fourier transform domain. Geophysics, 68(1): 355-369. DOI:10.1190/1.1543221
Hedefa M, Huang M Z, Zhu W H. 2010. Time domain trace interpolation for irregularly sampled seismic data by using a localized dip search approach. //80th Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts, 3569-3573.
Herrmann F J, Hennenfent G. 2008. Non-parametric seismic data recovery with curvelet frames. Geophysical Journal International, 173(1): 233-248. DOI:10.1111/j.1365-246X.2007.03698.x
Huo S D, Zhu W H, Shi T K. 2017. Iterative dip-steering median filter. Journal of Applied Geophysics, 144: 151-156. DOI:10.1016/j.jappgeo.2017.05.012
Lari H H, Naghizadeh M, Sacchi M D, et al. 2019. Adaptive singular spectrum analysis for seismic denoising and interpolation. Geophysics, 84(2): V133-V142. DOI:10.1190/geo2018-0350.1
Liu J C, Chou Y X, Zhu J J. 2018. Interpolating seismic data via the POCS method based on shearlet transform. Journal of Geophysics and Engineering, 15(3): 852-876. DOI:10.1088/1742-2140/aaa5d1
Liu Y, Fomel S, Liu C. 2015. Signal and noise separation in prestack seismic data using velocity-dependent seislet transform. Geophysics, 80(6): WD117-WD128. DOI:10.1190/geo2014-0234.1
Liu Y, Zhang P, Liu C, et al. 2018. Seismic data interpolation based on Bregman shaping iteration. Chinese Journal of Geophysics (in Chinese), 61(4): 1400-1412. DOI:10.6038/cjg2018L0422
Mallat S G, Zhang Z F. 1993. Matching pursuits with time-frequency dictionaries. IEEE Transactions on Signal Processing, 41(12): 3397-3415. DOI:10.1109/78.258082
Natarajan B K. 1995. Sparse approximate solutions to linear systems. Siam Journal on Computing, 24(2): 227-234. DOI:10.1137/S0097539792240406
Pan X, Cao S Y, Zu S H, et al. 2017. SVD-constrained MWNI with shaping theory. IEEE Geoscience and Remote Sensing Letters, 14(6): 856-860. DOI:10.1109/LGRS.2017.2683537
Spitz S. 1991. Seismic trace interpolation in the F-X domain. Geophysics, 56(6): 785-794. DOI:10.1190/1.1443096
Sun M M, Li Z C, Li Z N, et al. 2019. Reconstruction of seismic data with weighted MCA based on compressed sensing. Chinese Journal of Geophysics (in Chinese), 62(3): 1007-1021. DOI:10.6038/cjg2019L0336
Swindeman R, Fomel S. 2015. Seismic data interpolation using plane-wave shaping regularization. //85th Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts, 3853-3858.
Wang B F, Lu W K. 2017. Accurate and efficient seismic data interpolation using FK-curvelet transform. //87th Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts, 280-283.
Wang J Y, Liu Z H, Dai Z G. 2020. Logarithmic function-based sparse constraints for irregular seismic data Reconstruction. Progress in Geophysics (in Chinese), 35(6): 2228-2238. DOI:10.6038/pg2020DD0444
Wang Y H. 2002. Seismic trace interpolation in the f-x-y domain. Geophysics, 67(4): 1232-1239. DOI:10.1190/1.1500385
Yang P L, Gao J H, Chen W C. 2012. Curvelet-based POCS interpolation of nonuniformly sampled seismic records. Journal of Applied Geophysics, 79: 90-99. DOI:10.1016/j.jappgeo.2011.12.004
Yu Z, Ferguson J, McMechan G, et al. 2007. Wavelet-Radon domain dealiasing and interpolation of seismic data. Geophysics, 72(2): V41-V49. DOI:10.1190/1.2422797
Zu S H, Zhou H, Chen Y K, et al. 2016. Interpolating big gaps using inversion with slope constraint. IEEE Geoscience and Remote Sensing Letters, 13(9): 1369-1373. DOI:10.1109/LGRS.2016.2587301
邓盾, 裴健翔, 邓勇, 等. 2019. 南海深水盆地改善中深层地震成像关键技术研究. 地球物理学进展, 34(2): 702-708. DOI:10.6038/pg2019BB0551
刘洋, 张鹏, 刘财, 等. 2018. 地震数据Bregman整形迭代插值方法. 地球物理学报, 61(4): 1400-1412. DOI:10.6038/cjg2018L0422
孙苗苗, 李振春, 李志娜, 等. 2019. 基于压缩感知的加权MCA地震数据重构方法. 地球物理学报, 62(3): 1007-1021. DOI:10.6038/cjg2019L0336
王锦妍, 刘智慧, 代志刚. 2020. 基于对数函数稀疏约束的随机缺失地震数据的重建. 地球物理学进展, 35(6): 2228-2238. DOI:10.6038/pg2020DD0444