陳思遠(yuǎn),曹思遠(yuǎn)*,孫耀光,江雨濛,高君
1 中國石油大學(xué)(北京),北京 102249 2 中國石化石油勘探開發(fā)研究院,北京 100083
地震采集與處理過程中,受空采區(qū)和壞道的影響,會出現(xiàn)地震道的缺失,需進(jìn)行數(shù)據(jù)補(bǔ)全,為后續(xù)地震數(shù)據(jù)的反褶積、偏移以及儲層反演等工作提供優(yōu)質(zhì)數(shù)據(jù)體.最早使用的插值方法為FX域(Spitz, 1991)和F-XY域(Wang, 2002)的插值方法;考慮相鄰地震道的相似性,出現(xiàn)以差分算子為約束的插值方法(Liu et al., 2015; Zu et al., 2016),并認(rèn)為缺失地震數(shù)據(jù)左右兩側(cè)的平均為缺失處的數(shù)據(jù),由于疊前數(shù)據(jù)傾角較大,該方法效果不盡理想;由此,基于傾角約束(Hedefa et al., 2010)的插值方法便應(yīng)運(yùn)而生,包括基于傾角約束的中值濾波方法等(Huo et al., 2017).另一方面,在波動方程領(lǐng)域,有部分學(xué)者以速度為約束,以波動方程為媒介進(jìn)行插值等(Swindeman and Fomel, 2015).
近年,隨著壓縮感知恢復(fù)算法的發(fā)展,逐漸提出基于稀疏約束的插值算法(劉洋等,2018;孫苗苗等,2019),其假設(shè)無缺失數(shù)據(jù)經(jīng)過某種變換后,在變換域中呈稀疏態(tài),而缺失地震數(shù)據(jù)在變換域內(nèi)非稀疏,相比于無缺失數(shù)據(jù),缺失地震數(shù)據(jù)變換域內(nèi)的幅值較弱,可使用閾值壓制這些非稀疏能量以完成缺失數(shù)據(jù)的補(bǔ)全.目前,可用變換包括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)等.除此之外,低秩約束作為特殊的稀疏約束也被應(yīng)用到地震數(shù)據(jù)的插值中(Chen et al., 2019; Lari et al., 2019; Pan et al., 2017).
在稀疏約束目標(biāo)上,L0范數(shù)是標(biāo)準(zhǔn)的稀疏約束,該范數(shù)的求解是NP難問題(Natarajan, 1995),可使用L0范數(shù)的最佳凸近似L1范數(shù)代替L0范數(shù)求解(Chen et al., 1998),當(dāng)s表示稀疏度,且滿足δ2s 基于稀疏約束的地震數(shù)據(jù)插值方向,深入討論稀疏變換域結(jié)構(gòu)特性的研究論文較少,通常假設(shè)其整體的稀疏性,未考慮信號及變換域的特征,導(dǎo)致變換域中的弱能量被廣義閾值壓制.本文從頻率-波數(shù)(FK)域的結(jié)構(gòu)入手,探究沿頻率方向和沿波數(shù)方向有效能量及缺失數(shù)據(jù)分布特征,以減少變換域中的弱有效能量損失.在現(xiàn)有的研究基礎(chǔ)上開展如下工作:(1)以FK變換為例,使用圖示描述基于Lp范數(shù)的廣義閾值(下稱廣義閾值)所造成的被插值數(shù)據(jù)的高低頻缺失現(xiàn)象,探究在FK變換域中缺失數(shù)據(jù)沿頻率方向的形態(tài),構(gòu)造頻率加權(quán)閾值(FWT);(2)針對FK變換域的波數(shù)方向的聚集性特征,構(gòu)造結(jié)構(gòu)加權(quán)閾值(SWT)以提高稀疏約束能力;(3)模型測試部分,合成簡單模型分別測試“廣義閾值”、“頻率加權(quán)閾值(FWT)”、“結(jié)構(gòu)加權(quán)閾值(SWT)”、“頻率-結(jié)構(gòu)雙重加權(quán)閾值(FSWT)”在插值中的作用,以佐證理論部分的正確性;最后,以實(shí)際缺失地震數(shù)據(jù)測試本文所提出的頻率-結(jié)構(gòu)雙重加權(quán)閾值在缺失數(shù)據(jù)高低頻恢復(fù)上的有效性. 需要指出,除FK變換域可作為稀疏變換域之外,Curvelet變換、小波變換、Shearlet變換等均帶有頻率的分解,以滿足不同頻率下數(shù)據(jù)的稀疏表達(dá);因?yàn)榈卣饠?shù)據(jù)各頻率成分比重不同,所以本研究以FK變換作為稀疏變換旨在證明“頻率-結(jié)構(gòu)雙重加權(quán)閾值”相比于“廣義閾值”在數(shù)據(jù)高低頻恢復(fù)中的有效性,其他的帶有頻率分解的類似稀疏變換均可依照該思想構(gòu)造相應(yīng)的加權(quán)閾值. 時間采樣點(diǎn)為n的m道地震數(shù)據(jù)插值可以使用正演方程描述為: Y=LX, (1) 其中,Y∈mn×1表示觀測數(shù)據(jù),即含缺失地震道的數(shù)據(jù),L定義為式(2)所示的正演算子,其中,I和0分別表示單位矩陣和零矩陣,可從完整數(shù)據(jù)X中采樣得到觀測數(shù)據(jù).式(2)為: (2) 采樣矩陣L是秩虧矩陣,待求的完整數(shù)據(jù)X具有無窮個解,基于壓縮感知理論,可施加不同的正則化約束保證待求X的解唯一,并盡可能獲得最接近真實(shí)數(shù)據(jù)的解.常用的正則化約束包括Tikhonov約束(Fleming, 1990)(L2范數(shù)約束),稀疏約束(L0范數(shù)約束)等,可通過壓縮感知理論建立如下約束方程: (3) 其中,B表示某種變換,完整數(shù)據(jù)X在變換B的作用下可以符合Lq-Norm(0≤q<2)的約束,整形正則化理論表明,可以通過約束待求解數(shù)據(jù)變換域中的形態(tài)以獲得最接近真實(shí)數(shù)據(jù)的解;傳統(tǒng)方法假設(shè)完整數(shù)據(jù)在FK域中稀疏,而缺失后的數(shù)據(jù)非稀疏,通過整形正則化的方法求解該模型得到最接近真實(shí)數(shù)據(jù)的解.將式(3)寫為如式(4)所示的無約束優(yōu)化問題,并使用稀疏約束作為正則化項: (4) 式中,λ定義為正則化算子,由于L0-norm的求解是NP-hard問題,所使用的求解算法是匹配追蹤等(Mallat and Zhang, 1993),該類算法需預(yù)先給定稀疏度,而完整數(shù)據(jù)FK域的數(shù)據(jù)稀疏度無法確定,故匹配追蹤類算法不適用方程(4)的求解;考慮以Lp-norm近似L0-norm,(0 (5) 式(5)可以通過整形正則化的方法在多項式時間內(nèi)有效求解(Ge et al., 2011),給出迭代步長α,對應(yīng)于式(5)的迭代方程為: X(k+1)=B-1TλB[X(k)-αLT(LX(k)-Y)], (6) 式中,Tλ為式(7)定義的廣義軟閾值處理算子: (7) 1.2.1 含缺失道的地震數(shù)據(jù)在FK域的結(jié)構(gòu)特性 當(dāng)B為FK變換時,式(5)被稱為FK域的插值反演方程.地震數(shù)據(jù)因?yàn)槠漕l帶寬度有限,同相軸曲率較小,在FK域中呈聚集性分布;高波數(shù)、高頻率位置因地震數(shù)據(jù)構(gòu)造的特殊性,也會出現(xiàn)部分弱能量.而由于地震數(shù)據(jù)不規(guī)則缺失現(xiàn)象,缺失地震道的位置產(chǎn)生截斷,如圖1所示,造成沿波數(shù)方向的“頻譜泄漏”,隨缺失道的增加,“頻譜泄漏”的能量也隨之增加,缺失地震道的能量將類似于“白噪聲”均勻分布在沿波數(shù)方向上(圖1中F1處的虛線),可使用廣義閾值壓制(圖1下方虛線);另外,在沿頻率方向上,缺失數(shù)據(jù)并未造成頻譜泄漏(圖1中K1處的虛線),故而廣義閾值應(yīng)單獨(dú)作用于沿波數(shù)方向,不應(yīng)作用于沿頻率方向上. 圖1 缺失地震數(shù)據(jù)FK譜形態(tài)示意圖 圖2模擬廣義閾值對白噪聲(頻譜為白譜)和類子波頻譜形態(tài)的噪聲(頻譜形態(tài)類似為子波頻譜的形態(tài))的壓制效果及所帶來的問題;當(dāng)噪聲為白噪聲時(圖2a1),合適的廣義閾值可以有效的壓制噪聲的干擾(圖2d1和圖2d2),且可以保留高低頻能量;而當(dāng)噪聲為子波頻譜形態(tài)時(圖2a3),廣義閾值也壓制了噪聲的干擾,但是損失了高頻和低頻的能量(圖2d3和圖2d4). 圖2 不同頻譜形態(tài)的噪聲與常數(shù)閾值的測試示意圖 1.2.2 基于FK域的頻率-結(jié)構(gòu)雙重加權(quán)閾值的構(gòu)造 在地震數(shù)據(jù)處理中,低頻和高頻的能量是提高分辨率和儲層反演所需要的,廣義閾值在處理類子波頻譜形態(tài)的噪聲中損傷了能量較弱的低頻及高頻信號;對于類子波頻譜形態(tài)的噪聲,考慮使用類子波頻譜形態(tài)的閾值進(jìn)行處理,即能量弱的部分減小閾值,聯(lián)系上文缺失地震數(shù)據(jù)FK變換,沿波數(shù)方向的“頻譜泄漏”相當(dāng)于在每個波數(shù)上均增加了類子波頻譜形態(tài)的噪聲,在使用廣義閾值進(jìn)行處理時,低頻和高頻所對應(yīng)頻譜能量較弱的位置,廣義閾值會壓制弱能量,將造成缺失數(shù)據(jù)高低頻的損失,故而在基于FK域的地震數(shù)據(jù)插值中需要構(gòu)造呈類子波頻譜形態(tài)的沿頻率方向的閾值加權(quán)算子. 缺失數(shù)據(jù)造成的“頻譜泄漏”沿波數(shù)方向呈均勻隨機(jī)分布,可使用廣義閾值進(jìn)行壓制,當(dāng)缺失數(shù)目較多時,“頻譜泄漏”的能量逐漸逼近甚至超過有效信號的能量,這需要沿波數(shù)方向構(gòu)造一個隨FK譜結(jié)構(gòu)變化的加權(quán)算子,進(jìn)行沿波數(shù)方向的結(jié)構(gòu)加權(quán);基于頻率-結(jié)構(gòu)的雙重閾值加權(quán)算子的構(gòu)造方式如下. 記觀測數(shù)據(jù)Y經(jīng)過沿時間方向的Fourier變換得到的數(shù)據(jù)為Yf,x,式(1)可寫為: Yf,x=LXf,x, (8) 式中,Xf,x表示完成數(shù)據(jù)X經(jīng)過沿時間方向的Fourier變換得到的數(shù)據(jù),對于數(shù)據(jù)的每個頻率切片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中為稀疏態(tài),且Yfi,x在變換域B中為非稀疏,則可通過式(10)所示的壓縮感知重建方法得到每個頻率切片對應(yīng)的完整的數(shù)據(jù)Xfi,x,遍歷所有頻率切片即完成缺失數(shù)據(jù)的恢復(fù): (10) 根據(jù)式(3)—(6),約束方程(10)也可以使用整形正則化方法求解,即: (11) (12) 另外,考慮變換域的稀疏性,引入反加權(quán)算子hfi,x增強(qiáng)完整數(shù)據(jù)的每個頻率切片的稀疏約束能力,則式(10)改寫為: (13) 將頻率-結(jié)構(gòu)雙重加權(quán)閾值代入式(7)中,將廣義軟閾值修改為式(14),結(jié)合整形正則化的迭代式(6),可進(jìn)行基于頻率-結(jié)構(gòu)雙重加權(quán)閾值的缺失地震數(shù)據(jù)補(bǔ)全: (14) 本部分首先使用128×128的合成模型,分別測試廣義閾值,頻率加權(quán)閾值(FWT),結(jié)構(gòu)加權(quán)閾值(SWT)和頻率-結(jié)構(gòu)雙重加權(quán)閾值(FSWT)各自的作用;如圖3a所示,該模型從左到右Ricker子波主頻由15 Hz遞增到95 Hz;加入高斯白噪后信噪比為24.68 dB. 應(yīng)用上述四種閾值進(jìn)行插值,單次插值誤差剖面見圖3b—e,因?yàn)镾WT作為結(jié)構(gòu)加權(quán)閾值,是全頻帶的保幅閾值,誤差變化不明顯;同時,F(xiàn)WT重點(diǎn)傾向于弱頻率能量的保護(hù),在此模型中對應(yīng)于能量較弱的高頻,故而在全頻帶的誤差剖面上,四種閾值所產(chǎn)生的誤差形態(tài)接近.經(jīng)計算,四種閾值插值結(jié)果與原剖面的余弦相似度分別為:廣義閾值:99.27%,F(xiàn)WT:99.44%,SWT:99.41%,F(xiàn)SWT:99.56%. 為了直觀的體現(xiàn)頻率加權(quán)閾值對數(shù)據(jù)高頻補(bǔ)全的重要性,應(yīng)用時頻分析工具對四種閾值的插值結(jié)果及無缺失數(shù)據(jù)進(jìn)行頻率分解,分解得到的高頻剖面見圖4.廣義閾值和SWT的插值結(jié)果的高頻剖面表明基于廣義閾值的插值難以補(bǔ)全高頻成分(圖4b),而本文提出的FWT和FSWT的插值結(jié)果分頻剖面同相軸連續(xù)性較好(圖4c和圖4e),證明FWT和FSWT在缺失數(shù)據(jù)高頻恢復(fù)上的有效性. 因?yàn)槟P偷闹黝l是道號的函數(shù),本部分重復(fù)試驗(yàn)100次隨機(jī)缺失數(shù)據(jù)插值,繪制四種閾值插值結(jié)果逐道與原剖面(圖3a)的相似度(圖5),當(dāng)頻率處于中低頻時(<35 Hz),相似性FSWT>SWT>FWT=廣義閾值;在證明SWT的保幅效果優(yōu)異的同時,也說明FWT在能量較強(qiáng)的中低頻部分保幅效果較弱;在高頻部分(>70 Hz),相似性FSWT>FWT>SWT>廣義閾值,說明FWT在數(shù)據(jù)高頻部分保幅性優(yōu)異;而SWT在整個頻帶的相似性排名均優(yōu)于廣義閾值,證明SWT在全頻帶具有一定保幅作用. 圖3 原始數(shù)據(jù)及插值誤差 圖4 分頻剖面(高頻) 圖5 四種閾值插值結(jié)果與無缺失數(shù)據(jù)的余弦相似度 繼續(xù)合成256道,含256個采樣點(diǎn)的疊前地震數(shù)據(jù);如圖6所示為原始數(shù)據(jù)、缺失50%地震道的地震數(shù)據(jù)和缺失地震數(shù)據(jù)的FK譜;添加隨機(jī)噪聲后,信噪比為32 dB.本部分將測試廣義閾值和頻率-結(jié)構(gòu)雙重加權(quán)閾值(FSWT)對圖6所示復(fù)雜模型的適用性.測試環(huán)境為Inter i7-9750H處理器,16 GB內(nèi)存,軟件環(huán)境為Matlab2020a. 圖6 合成模型 對于插值反演方程,正則化參數(shù)對插值結(jié)果影響相對較大.分別選擇不同的正則化參數(shù)進(jìn)行插值測試,記錄插值后與原始數(shù)據(jù)的余弦相似度、插值后與原始數(shù)據(jù)低頻剖面以及高頻剖面的余弦相似度和計算時間的曲線如圖7所示:當(dāng)正則化參數(shù)逐漸變小時,任何頻帶范圍內(nèi),基于廣義閾值的插值方法和基于FSWT的插值方法所計算的余弦相似度趨于相同(圖7a—c),但是計算量呈指數(shù)型增加(圖7d);而當(dāng)正則化參數(shù)逐漸變大時,基于廣義閾值的插值方法的相似度快速下降,而基于FSWT的插值方法余弦相似度變化不明顯(圖7a—c),同時,計算時間趨于穩(wěn)定(圖7d),綜合考慮計算成本,如圖8所示取正則化參數(shù)等于8時的插值結(jié)果進(jìn)行分析. 圖7 正則化參數(shù)對插值性能的影響測試 基于廣義閾值的插值方法的誤差剖面上(圖8b)可見大量同相軸,且誤差剖面的FK譜上有較多聚集性能量(圖8c),插值結(jié)果較差;相比于前者,基于FSWT的插值方法的誤差剖面無明顯可見同相軸,F(xiàn)K域中的聚集能量也較少(圖8d). 將圖6a所示的無缺失數(shù)據(jù)進(jìn)行分頻,取低頻(主頻約18 Hz)及高頻(主頻約66 Hz)的剖面如圖9所示;同時將圖8a和圖8d所示的插值結(jié)果也進(jìn)行分頻,繪制高低頻插值結(jié)果及誤差剖面如圖10所示;無論是低頻或者高頻部分,基于廣義閾值的插值方法的分頻誤差剖面上(圖10b和圖10d)均可見同相軸,且缺失數(shù)據(jù)未成功恢復(fù);基于FSWT的插值方法低頻及高頻誤差均較小,無缺失地震道現(xiàn)象出現(xiàn)(圖10f和圖10h). 圖8 基于廣義閾值和FSWT的插值結(jié)果 圖9 無缺失數(shù)據(jù)分頻剖面 實(shí)際數(shù)據(jù)測試部分,選擇M地區(qū)實(shí)際數(shù)據(jù)進(jìn)行插值處理,該數(shù)據(jù)含有180道,采樣率為4 ms,本文置零約35%的地震道以測試廣義閾值和FSWT的插值效果.如圖所示,圖11a和圖11b分別為原始數(shù)據(jù)和缺失地震道的數(shù)據(jù);圖12為使用廣義閾值和FSWT的插值結(jié)果與誤差,從插值剖面上無法看出兩種閾值的區(qū)別(圖12a和圖12c),在相應(yīng)的誤差剖面上(圖12b和圖12d)FSWT的插值誤差小于使用廣義閾值的插值誤差;同時,本文計算使用兩種閾值插值后剖面缺失位置所占無缺失數(shù)據(jù)缺失位置能量的比重分別為:廣義閾值:75.7%;頻率-結(jié)構(gòu)雙重加權(quán)閾值:89.9%,這證明本文所提出的FSWT在全頻帶上的保幅性優(yōu)越于傳統(tǒng)廣義閾值. 圖11 實(shí)際單炮數(shù)據(jù)和缺失35%地震道的數(shù)據(jù) 圖12 基于廣義閾值和FSWT的插值剖面和誤差 圖11a、圖12a和圖12c的缺失地震數(shù)據(jù)位置的頻譜如圖13所示,基于FSWT插值的結(jié)果較好的擬合原始數(shù)據(jù)的頻譜,而基于廣義閾值的插值結(jié)果在低頻和高頻與原始數(shù)據(jù)頻譜的相似度均較弱;為了得到更準(zhǔn)確的測試結(jié)果,應(yīng)用時頻分析工具對原始數(shù)據(jù)(圖11a)以及插值后的數(shù)據(jù)圖12a和圖12c進(jìn)行頻率分解,低頻(主頻約為9 Hz)及高頻剖面(主頻約為90 Hz)如圖14和圖15所示,相比于圖14所示的原始數(shù)據(jù)分頻結(jié)果,本文提出的FSWT無論是低頻(圖15b)或者高頻(圖15d)的插值效果均優(yōu)異于廣義閾值的插值效果(圖15a和圖15c). 圖13 實(shí)際單炮數(shù)據(jù)和插值后數(shù)據(jù)的頻譜 圖14 無缺失數(shù)據(jù)的高低頻剖面 圖15 基于廣義閾值和FSWT插值后的高低頻剖面 針對基于傳統(tǒng)廣義閾值插值所造成高低頻難以補(bǔ)全的問題,本文考慮FK域中缺失數(shù)據(jù)頻率方向的非稀疏,構(gòu)造類子波形狀的頻率加權(quán)算子,同時注意到波數(shù)域有效能量聚集性,聯(lián)合頻率加權(quán)算子,提出頻率-結(jié)構(gòu)雙重加權(quán)算子進(jìn)行缺失地震數(shù)據(jù)的補(bǔ)全,一定程度上解決了傳統(tǒng)廣義閾值對于缺失位置高低頻難以恢復(fù)的問題,保證后續(xù)處理環(huán)節(jié)的進(jìn)行.算法在適用性方面仍有較大研究空間: (1)基于頻率-結(jié)構(gòu)雙重加權(quán)閾值主要作用在于恢復(fù)缺失位置處、能量較弱高低頻信息,由于中頻段(主頻附近)的FSWT與廣義閾值的大小、形態(tài)均接近,所以在該頻段內(nèi)算法改善效果不佳. (2)頻率域的非離散現(xiàn)象不僅在FK域中出現(xiàn),包括Curvelet域、Shearlet域等涉及頻率的變換域均有此現(xiàn)象,也可通過構(gòu)造頻率加權(quán)算子提高這些變換域的插值精度. (3)對于三維數(shù)據(jù),可采用頻率-波數(shù)-波數(shù)(FKK)或者3D-Curvelet作為稀疏變換域,按本文方法可構(gòu)造三維加權(quán)閾值進(jìn)行高低頻信息的有效恢復(fù). (4)本文僅在FK域進(jìn)行加權(quán)算子的可行性研究,由于地質(zhì)構(gòu)造通常較復(fù)雜,可能含有繞射波等大傾角構(gòu)造,由于變換域限制,大傾角同相軸補(bǔ)全效果不佳,可使用其他稀疏變換域(Curvelet變換)聯(lián)合頻率-結(jié)構(gòu)雙重加權(quán)閾值對這部分弱能量進(jìn)行有效恢復(fù).1 理論
1.1 基于稀疏反演的地震數(shù)據(jù)插值理論
1.2 頻率-結(jié)構(gòu)雙重加權(quán)閾值的構(gòu)造
2 模型測試
3 實(shí)際數(shù)據(jù)測試
4 結(jié)論