基于頻率-結構雙重加權閾值的地震數據插值

陳思遠，曹思遠*，孫耀光，江雨濛，高君

1 中國石油大學(北京)，北京 102249 2 中國石化石油勘探開發研究院，北京 100083

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

基于稀疏約束的地震數據插值方向，深入討論稀疏變換域結構特性的研究論文較少，通常假設其整體的稀疏性，未考慮信號及變換域的特征，導致變換域中的弱能量被廣義閾值壓制.本文從頻率-波數(FK)域的結構入手，探究沿頻率方向和沿波數方向有效能量及缺失數據分布特征，以減少變換域中的弱有效能量損失.在現有的研究基礎上開展如下工作：(1)以FK變換為例，使用圖示描述基于Lp范數的廣義閾值(下稱廣義閾值)所造成的被插值數據的高低頻缺失現象，探究在FK變換域中缺失數據沿頻率方向的形態，構造頻率加權閾值(FWT)；(2)針對FK變換域的波數方向的聚集性特征，構造結構加權閾值(SWT)以提高稀疏約束能力；(3)模型測試部分，合成簡單模型分別測試“廣義閾值”、“頻率加權閾值(FWT)”、“結構加權閾值(SWT)”、“頻率-結構雙重加權閾值(FSWT)”在插值中的作用，以佐證理論部分的正確性；最后，以實際缺失地震數據測試本文所提出的頻率-結構雙重加權閾值在缺失數據高低頻恢復上的有效性.

需要指出，除FK變換域可作為稀疏變換域之外，Curvelet變換、小波變換、Shearlet變換等均帶有頻率的分解，以滿足不同頻率下數據的稀疏表達；因為地震數據各頻率成分比重不同，所以本研究以FK變換作為稀疏變換旨在證明“頻率-結構雙重加權閾值”相比于“廣義閾值”在數據高低頻恢復中的有效性，其他的帶有頻率分解的類似稀疏變換均可依照該思想構造相應的加權閾值.

1 理論

1.1 基于稀疏反演的地震數據插值理論

時間采樣點為n的m道地震數據插值可以使用正演方程描述為：

Y=LX，

(1)

其中，Y∈mn×1表示觀測數據，即含缺失地震道的數據，L定義為式(2)所示的正演算子，其中，I和0分別表示單位矩陣和零矩陣，可從完整數據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

(5)

式(5)可以通過整形正則化的方法在多項式時間內有效求解(Ge et al., 2011)，給出迭代步長α，對應于式(5)的迭代方程為：

X(k+1)=B-1TλB[X(k)-αLT(LX(k)-Y)]，

(6)

式中，Tλ為式(7)定義的廣義軟閾值處理算子：

(7)

1.2 頻率-結構雙重加權閾值的構造

1.2.1 含缺失道的地震數據在FK域的結構特性

當B為FK變換時，式(5)被稱為FK域的插值反演方程.地震數據因為其頻帶寬度有限，同相軸曲率較小，在FK域中呈聚集性分布；高波數、高頻率位置因地震數據構造的特殊性，也會出現部分弱能量.而由于地震數據不規則缺失現象，缺失地震道的位置產生截斷，如圖1所示，造成沿波數方向的“頻譜泄漏”，隨缺失道的增加，“頻譜泄漏”的能量也隨之增加，缺失地震道的能量將類似于“白噪聲”均勻分布在沿波數方向上(圖1中F1處的虛線)，可使用廣義閾值壓制(圖1下方虛線)；另外，在沿頻率方向上，缺失數據并未造成頻譜泄漏(圖1中K1處的虛線)，故而廣義閾值應單獨作用于沿波數方向，不應作用于沿頻率方向上.

圖1 缺失地震數據FK譜形態示意圖

圖2模擬廣義閾值對白噪聲(頻譜為白譜)和類子波頻譜形態的噪聲(頻譜形態類似為子波頻譜的形態)的壓制效果及所帶來的問題；當噪聲為白噪聲時(圖2a1)，合適的廣義閾值可以有效的壓制噪聲的干擾(圖2d1和圖2d2)，且可以保留高低頻能量；而當噪聲為子波頻譜形態時(圖2a3)，廣義閾值也壓制了噪聲的干擾，但是損失了高頻和低頻的能量(圖2d3和圖2d4).

圖2 不同頻譜形態的噪聲與常數閾值的測試示意圖

1.2.2 基于FK域的頻率-結構雙重加權閾值的構造

在地震數據處理中，低頻和高頻的能量是提高分辨率和儲層反演所需要的，廣義閾值在處理類子波頻譜形態的噪聲中損傷了能量較弱的低頻及高頻信號；對于類子波頻譜形態的噪聲，考慮使用類子波頻譜形態的閾值進行處理，即能量弱的部分減小閾值，聯系上文缺失地震數據FK變換，沿波數方向的“頻譜泄漏”相當于在每個波數上均增加了類子波頻譜形態的噪聲，在使用廣義閾值進行處理時，低頻和高頻所對應頻譜能量較弱的位置，廣義閾值會壓制弱能量，將造成缺失數據高低頻的損失，故而在基于FK域的地震數據插值中需要構造呈類子波頻譜形態的沿頻率方向的閾值加權算子.

缺失數據造成的“頻譜泄漏”沿波數方向呈均勻隨機分布，可使用廣義閾值進行壓制，當缺失數目較多時，“頻譜泄漏”的能量逐漸逼近甚至超過有效信號的能量，這需要沿波數方向構造一個隨FK譜結構變化的加權算子，進行沿波數方向的結構加權；基于頻率-結構的雙重閾值加權算子的構造方式如下.

記觀測數據Y經過沿時間方向的Fourier變換得到的數據為Yf,x，式(1)可寫為：

Yf,x=LXf,x，

(8)

式中，Xf,x表示完成數據X經過沿時間方向的Fourier變換得到的數據，對于數據的每個頻率切片fi,(i=1,2,…,m),式(8)可寫為：

Yfi,x=L0Xfi,x,i=1,2,…,m，

(9)

其中，L0=diag{1,1,…,0,…,1}n×n仍為采樣算子，若Xfi,x在變換域B中為稀疏態，且Yfi,x在變換域B中為非稀疏，則可通過式(10)所示的壓縮感知重建方法得到每個頻率切片對應的完整的數據Xfi,x,遍歷所有頻率切片即完成缺失數據的恢復：

(10)

根據式(3)—(6)，約束方程(10)也可以使用整形正則化方法求解，即：

(11)

(12)

另外，考慮變換域的稀疏性，引入反加權算子hfi,x增強完整數據的每個頻率切片的稀疏約束能力，則式(10)改寫為：

(13)

將頻率-結構雙重加權閾值代入式(7)中，將廣義軟閾值修改為式(14)，結合整形正則化的迭代式(6)，可進行基于頻率-結構雙重加權閾值的缺失地震數據補全：

(14)

2 模型測試

本部分首先使用128×128的合成模型，分別測試廣義閾值，頻率加權閾值(FWT)，結構加權閾值(SWT)和頻率-結構雙重加權閾值(FSWT)各自的作用；如圖3a所示，該模型從左到右Ricker子波主頻由15 Hz遞增到95 Hz；加入高斯白噪后信噪比為24.68 dB.

應用上述四種閾值進行插值，單次插值誤差剖面見圖3b—e，因為SWT作為結構加權閾值，是全頻帶的保幅閾值，誤差變化不明顯；同時，FWT重點傾向于弱頻率能量的保護，在此模型中對應于能量較弱的高頻，故而在全頻帶的誤差剖面上，四種閾值所產生的誤差形態接近.經計算，四種閾值插值結果與原剖面的余弦相似度分別為：廣義閾值：99.27%，FWT：99.44%，SWT：99.41%，FSWT：99.56%.

為了直觀的體現頻率加權閾值對數據高頻補全的重要性，應用時頻分析工具對四種閾值的插值結果及無缺失數據進行頻率分解，分解得到的高頻剖面見圖4.廣義閾值和SWT的插值結果的高頻剖面表明基于廣義閾值的插值難以補全高頻成分(圖4b)，而本文提出的FWT和FSWT的插值結果分頻剖面同相軸連續性較好(圖4c和圖4e)，證明FWT和FSWT在缺失數據高頻恢復上的有效性.

因為模型的主頻是道號的函數，本部分重復試驗100次隨機缺失數據插值，繪制四種閾值插值結果逐道與原剖面(圖3a)的相似度(圖5)，當頻率處于中低頻時(<35 Hz)，相似性FSWT>SWT>FWT=廣義閾值；在證明SWT的保幅效果優異的同時，也說明FWT在能量較強的中低頻部分保幅效果較弱；在高頻部分(>70 Hz)，相似性FSWT>FWT>SWT>廣義閾值，說明FWT在數據高頻部分保幅性優異；而SWT在整個頻帶的相似性排名均優于廣義閾值，證明SWT在全頻帶具有一定保幅作用.

圖3 原始數據及插值誤差

圖4 分頻剖面(高頻)

圖5 四種閾值插值結果與無缺失數據的余弦相似度

繼續合成256道，含256個采樣點的疊前地震數據；如圖6所示為原始數據、缺失50%地震道的地震數據和缺失地震數據的FK譜；添加隨機噪聲后，信噪比為32 dB.本部分將測試廣義閾值和頻率-結構雙重加權閾值(FSWT)對圖6所示復雜模型的適用性.測試環境為Inter i7-9750H處理器，16 GB內存，軟件環境為Matlab2020a.

圖6 合成模型

對于插值反演方程，正則化參數對插值結果影響相對較大.分別選擇不同的正則化參數進行插值測試，記錄插值后與原始數據的余弦相似度、插值后與原始數據低頻剖面以及高頻剖面的余弦相似度和計算時間的曲線如圖7所示：當正則化參數逐漸變小時，任何頻帶范圍內，基于廣義閾值的插值方法和基于FSWT的插值方法所計算的余弦相似度趨于相同(圖7a—c)，但是計算量呈指數型增加(圖7d)；而當正則化參數逐漸變大時，基于廣義閾值的插值方法的相似度快速下降，而基于FSWT的插值方法余弦相似度變化不明顯(圖7a—c)，同時，計算時間趨于穩定(圖7d)，綜合考慮計算成本，如圖8所示取正則化參數等于8時的插值結果進行分析.

圖7 正則化參數對插值性能的影響測試

基于廣義閾值的插值方法的誤差剖面上(圖8b)可見大量同相軸，且誤差剖面的FK譜上有較多聚集性能量(圖8c)，插值結果較差；相比于前者，基于FSWT的插值方法的誤差剖面無明顯可見同相軸，FK域中的聚集能量也較少(圖8d).

將圖6a所示的無缺失數據進行分頻，取低頻(主頻約18 Hz)及高頻(主頻約66 Hz)的剖面如圖9所示；同時將圖8a和圖8d所示的插值結果也進行分頻，繪制高低頻插值結果及誤差剖面如圖10所示；無論是低頻或者高頻部分，基于廣義閾值的插值方法的分頻誤差剖面上(圖10b和圖10d)均可見同相軸，且缺失數據未成功恢復；基于FSWT的插值方法低頻及高頻誤差均較小，無缺失地震道現象出現(圖10f和圖10h).

圖8 基于廣義閾值和FSWT的插值結果

圖9 無缺失數據分頻剖面

3 實際數據測試

實際數據測試部分，選擇M地區實際數據進行插值處理，該數據含有180道，采樣率為4 ms，本文置零約35%的地震道以測試廣義閾值和FSWT的插值效果.如圖所示，圖11a和圖11b分別為原始數據和缺失地震道的數據；圖12為使用廣義閾值和FSWT的插值結果與誤差，從插值剖面上無法看出兩種閾值的區別(圖12a和圖12c)，在相應的誤差剖面上(圖12b和圖12d)FSWT的插值誤差小于使用廣義閾值的插值誤差；同時，本文計算使用兩種閾值插值后剖面缺失位置所占無缺失數據缺失位置能量的比重分別為：廣義閾值：75.7%；頻率-結構雙重加權閾值：89.9%，這證明本文所提出的FSWT在全頻帶上的保幅性優越于傳統廣義閾值.

圖11 實際單炮數據和缺失35%地震道的數據

圖12 基于廣義閾值和FSWT的插值剖面和誤差

圖11a、圖12a和圖12c的缺失地震數據位置的頻譜如圖13所示，基于FSWT插值的結果較好的擬合原始數據的頻譜，而基于廣義閾值的插值結果在低頻和高頻與原始數據頻譜的相似度均較弱；為了得到更準確的測試結果，應用時頻分析工具對原始數據(圖11a)以及插值后的數據圖12a和圖12c進行頻率分解，低頻(主頻約為9 Hz)及高頻剖面(主頻約為90 Hz)如圖14和圖15所示，相比于圖14所示的原始數據分頻結果，本文提出的FSWT無論是低頻(圖15b)或者高頻(圖15d)的插值效果均優異于廣義閾值的插值效果(圖15a和圖15c).

圖13 實際單炮數據和插值后數據的頻譜

圖14 無缺失數據的高低頻剖面

圖15 基于廣義閾值和FSWT插值后的高低頻剖面

4 結論

針對基于傳統廣義閾值插值所造成高低頻難以補全的問題，本文考慮FK域中缺失數據頻率方向的非稀疏，構造類子波形狀的頻率加權算子，同時注意到波數域有效能量聚集性，聯合頻率加權算子，提出頻率-結構雙重加權算子進行缺失地震數據的補全，一定程度上解決了傳統廣義閾值對于缺失位置高低頻難以恢復的問題，保證后續處理環節的進行.算法在適用性方面仍有較大研究空間：

(1)基于頻率-結構雙重加權閾值主要作用在于恢復缺失位置處、能量較弱高低頻信息，由于中頻段(主頻附近)的FSWT與廣義閾值的大小、形態均接近，所以在該頻段內算法改善效果不佳.

(2)頻率域的非離散現象不僅在FK域中出現，包括Curvelet域、Shearlet域等涉及頻率的變換域均有此現象，也可通過構造頻率加權算子提高這些變換域的插值精度.

(3)對于三維數據，可采用頻率-波數-波數(FKK)或者3D-Curvelet作為稀疏變換域，按本文方法可構造三維加權閾值進行高低頻信息的有效恢復.

(4)本文僅在FK域進行加權算子的可行性研究，由于地質構造通常較復雜，可能含有繞射波等大傾角構造，由于變換域限制，大傾角同相軸補全效果不佳，可使用其他稀疏變換域(Curvelet變換)聯合頻率-結構雙重加權閾值對這部分弱能量進行有效恢復.