高維地震數據Wiener中心濾波方法

張力起,張 猛,2,王華忠,秦廣勝

(1.波現象與智能反演成像研究組(WPI),同濟大學海洋與地球科學學院,上海200092;2.中國石油化工股份有限公司勝利油田分公司物探研究院,山東東營257000;3.中國石油化工股份有限公司中原油田分公司物探研究院,河南濮陽457001)

近十余年,油氣地震勘探技術的顯著進展主要體現在“兩寬一高”地震數據采集技術的普遍使用和地震波反演成像方法技術的逐漸實用化。噪聲在地震波反演成像中影響巨大,它決定了反演成像方法的收斂性及成像質量?！皟蓪捯桓摺钡卣饠祿杉玫搅?盡可能)無假頻且高維的疊前地震數據,為在高維數據(信號)空間中進行徹底的噪聲壓制提供了良好的數據基礎。

勘探地震中的地震波場,可以表達成信號、相干噪聲和隨機噪聲3部分的疊加,可記為D=S+Nc+Nr,其中,D為觀測數據,S為地震信號,Nc為相干噪聲,Nr為隨機噪聲。噪聲壓制的思想框架通常是將含噪信號視為Hilbert空間中的一個平方可積函數,選取合適的基函數對信號進行稀疏、近似的表達(信號預測的正問題);對上述信號預測正問題,在一定的估計理論基礎上,求解最佳的信號預測器,從而實現信號的最佳估計。在這樣的函數逼近思想下,地震數據的去噪、插值、壓縮和特征提取等正問題是統一的。

在Bayes估計理論框架下,假定信號可以用既定的預測器進行預測,噪聲滿足某種統計分布,信號自身滿足一些先驗信息(例如具有稀疏性),可以通過后驗概率最大化估計出最佳預測信號。若假定信號是線性可預測的且噪聲滿足高斯分布,那么可以在均方誤差或L2范數誤差最小準則下對信號實現最佳預測。常見的方法包括f-x濾波[1]、奇異譜分析(依據觀測數據的自相關函數分解)[2]。信號在一些變換空間如小波基、曲波基、物理小波基和字典基[3-5]中,具有更好的稀疏性,將其與L0/L1范數下稀疏反演理論[6]結合可以更好地實現信號預測,達到信號分析的目的。這是當前現代信號分析的重點發展方向。此外,按照地震數據(隨機信號)的統計特征設計濾波器是另一種重要的濾波方法,例如高斯濾波器、中值濾波器。將此類統計濾波方法擴展到高維數據的情況是必然的,例如非局部均值濾波器[7]、三維塊匹配濾波器[8-9],此時數據或圖像中結構信息在構建統計濾波器時是非常關鍵的。據此,基于信號逼近理論的濾波方法和基于信號結構信息的統計濾波方法逐漸有融合的趨勢。

基于AR模型的預測濾波方法通常是在頻率空間域(f-x)采用前向或后向預測濾波方法,或利用前向(后向)預測算子的共軛對稱性同時進行后向(前向)預測濾波,從而達到壓制噪聲的目的。GüLüNAY[10]將非因果濾波(類似高維Wiener中心濾波)應用于三維疊后數據的去噪;WANG[11]將f-x域地震數據插值推廣至三維的f-x-y域進行數據插值;LIU等[12]利用ARMA模型和非因果空間預測算子進行f-x-y濾波。

“兩寬一高”地震勘探中采集的巨量地震數據(一塊上百平方千米的工區疊前數據達到十幾到幾十TB),對高維空間、高效且高質量的噪聲壓制技術的研發有著迫切的需求。本文未采用更復雜的信號預測器(譬如高維小波變換),而是將AR模型的f-x濾波器(前向或后向預測濾波器)修改為Wiener中心預測濾波器,同時將高維數據排成Hankel矩陣的形式使其與Wiener中心預測濾波器對應,在預測誤差的L2范數定義下求解對應的法方程獲得最佳Wiener中心濾波器系數來構建高維線性預測濾波器,實現最佳濾波。對于含有復雜波場的地震數據,為滿足局部線性假設,采用對數據局部取窗的方式進行去噪。數據窗長度以及濾波算子長度的選取會明顯地影響去噪效果和執行效率。為了提高本文方法的計算效率,利用MPI進程級并行和OpenMP線程級并行來加速計算,將本文方法應用于合成數據和實際數據,結果表明其有良好的去噪效果。

1 Bayes框架下的地震數據分析理論

地震數據預處理(地震信號處理)中包含的大多數問題,譬如去噪、插值、多次波壓制等正問題的構建都可以在Bayes反演理論框架下統一。實際采集的地震數據被認為是隨機信號,滿足一定的概率分布。地震數據去噪的正問題是建立對數據中所蘊含信號的預測理論,通常利用函數逼近理論建立包含參數的信號預測模型。預測結果與實測數據的殘差是隨機的,滿足一定的先驗概率分布。因此,隨機地震觀測數據可表示為:

(1)

(2)

式中:P(m|dobs)=P(dobs|m)P(m)/P(dobs),P(m|dobs)為后驗概率密度函數,其中,P(dobs|m)為實測數據dobs對模型參數m的先驗概率密度函數,P(m)為模型參數m的先驗概率密度函數,P(dobs)為實測數據dobs的先驗概率密度函數。P(m|dobs)可以理解為當觀測數據為實測數據dobs時,對應不同的參數m時的概率。在高維情況下,后驗概率密度函數P(m|dobs)幾乎無法由實際計算得到。因此,常取后驗概率密度函數P(m|dobs)最大值時對應的參數m為最終的估計結果:

(3)

一般地,在假設P(dobs|m)和P(m)這兩個概率密度函數都是高斯函數且正問題是線性的情況下,有:

(4)

式中:C為常數;S(m)為代價函數。

其中,代價函數S(m)為:

(5)

式中:mprior是模型參數m的先驗值;CD是數據自相關矩陣;CM是模型參數自相關矩陣。采用各種數值優化算法求解公式(5),均可得到估計的模型參數m。

信號預測算子G(·)可記為G,它由各種選定的基函數族中的基函數線性疊加產生。基函數族可分為兩類,第1類是預先選定的基函數簇,包括Fourier基函數、余弦基函數、Radon基函數、小波基函數和曲波基函數等;第2類是由數據驅動獲得的基函數簇,包括K-L變換、奇異譜分析和字典學習獲得的基函數等。模型參數m為線性信號預測器的系數?；瘮祵嵸|上是信號的特征成分?；瘮颠x擇恰當,則可以用很少的特征成分組合很好地表達信號,這是信號(與圖像)分析中最根本的基礎。在上述Bayes框架下,對模型參數進行最佳估計后,即可實現對數據中所蘊含信號的最佳表達。

上述抽象的理論框架建立了數據分析的理論基礎,原則上適用于任何信號和圖像分析。本文中提出的高維Wiener中心濾波方法則是基于信號是線性可預測且噪聲滿足高斯分布的假設,在L2范數誤差最小準則下對信號進行建模和最佳預測。

2 f-x域中線性信號模型的建立與信號預測

Wiener濾波可實現對線性信號或相干噪聲的最佳預測,達到壓制隨機噪聲(也包括相干噪聲)和提高信噪比的目的。在頻率空間域,這些線性信號可以表示成具有線性相位移的諧波疊加。時間空間域中的二維信號s(t,x)在f-x域表示為:

(6)

式中:p為水平射線參數;Δx為道間距;m為道數;Δt為相鄰兩道間的時移量;A(f)為子波的頻譜;φ=2πfpΔx,φ代表由p確定的線性同相軸的線性相位移。

對于含有一個線性同相軸的任意一個單頻,相鄰道之間存在如下關系:

(7)

式中:a=ei2πfpΔx;Si-1為頻率域第i-1道的值;Si為頻率域第i道的值。(7)式說明各道信號之間存在可預測性。類比AR模型,當含有p個線性同相軸時,可以得到前向線性預測濾波器{gi},i∈[1,p]:

(8)

(9)

(10)

圖1 算子長度為3時前向預測濾波(圓圈)和后向預測濾波(方框)示意

(11)

(12)

為求解濾波器系數g,需要建立如下誤差泛函J(g):

(13)

(13)式的法方程為:

(14)

(15)

式中:Df為前向預測濾波器構建的矩陣;gf為前向預測濾波器系數;Db為后向預測濾波器構建的矩陣;gb為后向預測濾波器系數。

3 高維Wiener中心濾波

根據實際數據各自構建法方程,然后求解得到的前向預測濾波器和后向預測濾波器(基于AR模型構建),只適用于噪聲滿足高斯分布、振幅不隨空間變化的局部線性信號預測,若不滿足該條件,則在壓制噪聲的同時會破壞信號的振幅和相位信息。Wiener中心預測方法適用于振幅緩變的局部線性信號。利用數據的多級Hankel矩陣排布來構建Wiener中心預測濾波器的法方程,進而實現高維數據的預測濾波,該方法稱為高維Wiener中心濾波方法。

3.1 Hankel矩陣的構建

對于一維數據(向量),其構建的Hankel矩陣是指每一條交叉對角線上的元素都相等的矩陣,即Hankel矩陣的項hi,j=hi-1,j+1。對于高維數據,也同樣能夠構建出Hankel矩陣,即塊狀Hankel矩陣(block Hankel matrix)。多維信號以二維數據為例,D∈Rm×n所代表的矩陣為:

(16)

(17)

(18)

式中:l2+k2-1=m。圖2為將3×3的二維數據排成塊狀Hankel矩陣結構。

圖2 由二維數據構建的塊狀Hankel矩陣結構

3.2 Wiener濾波中心預測濾波

將二維數據變換到f-x域后,單個頻率片數據為一維數據。采用圖3中的一維Wiener中心預測濾波方法進行預測濾波。假定數據長度為6,Wiener中心濾波算子長度為2,一維Wiener中心濾波器的構建過程如圖4所示,可以看出,只需要對一維數據進行Hankel排布就可以得到公式(14)中的D和S。

構建多級Hankel矩陣可以將一維Wiener中心濾波方法推廣到高維。以二維Wiener中心預測濾波

s

s

k,l

(19)

圖5 二維Wiener中心預測濾波示意(算子大小為5×5)

圖6 根據單頻空間二維數據矩陣構建D和向量d

若是按行優先的方式進行數據的Hankel矩陣排布,首先對窗內(紅色框)數據按行優先方式排成一行;然后將窗沿行方向滑動一個單位長度;再將窗內數據按行優先方式排成一行,依照上述方式滑動窗并取數據進行排布;直到窗滑動到數據末端;再將窗按列滑動一個單位長度,重復上述按行方向滑動窗取數據并排布的操作;重復以上過程直到窗滑動到數據末端。

高維Wiener中心濾波的計算步驟與f-x濾波計算步驟類似,具體如下:

1) 將數據從時間空間域變換到頻率空間域;

2) 對每個頻率切片數據進行Hankel矩陣排布,構造法方程,估計濾波器系數;

3) 利用濾波器系數對變換后的復數域數據進行褶積運算;

4) 將濾波后的結果變換回時間空間域。

為了滿足信號的空間局部線性假設條件,對數據進行取窗處理。因三維地震數據的單個頻率片數據為二維復數域(如圖7藍色框所示),故對二維復數域數據取窗數據(圖7中黑色框所示),對缺少的數據進行補零后,將其排成Hankel矩陣,最終濾波后輸出數據的區域如圖7紅色框所示,具體算法流程如圖8所示。

圖7 二維單個頻率片數據的I/O示意

圖8 實際數據處理流程

為了滿足局部線性假設條件,需對地震數據進行取窗處理。不同的數據窗長度和不同的算子長度會造成去噪結果和效率的明顯差異。長數據窗會破壞局部線性假設,無益于提高去噪能力,同時還會增加計算量;而短數據窗則導致該算法去噪能力下降,無法適應復雜波場的情況。因此,需要根據數據波場的復雜程度選取合適的數據窗長度和算子長度。一般情況下,算子長度選取為3或5,空間窗長度選取范圍為15～30,時間窗長度選取范圍為100～200個采樣點。由于采用了分塊處理數據的方式并且預測濾波是在單個頻率片上進行的,故可以采用并行策略加速運算。并行主要利用MPI和OpenMP進行加速。利用MPI多進程分塊讀取高維窗數據體,各個進程對各自的數據體獨立進行高維去噪;在每個進程均利用OpenMP多線程分別處理單個頻率片數據。

4 測試應用及分析

對三維數據體取窗后獲得的數據體同樣是三維的,而頻率空間域中濾波算子的維度則是二維的。公式(20)為信噪比(SNR)的定義。

(20)

式中:{xi},i∈[1,N]為無噪聲數據;{ni},i∈[1,N]為含噪聲數據;N為數據長度。

4.1 三維合成地震數據

圖9中合成的三維數據大小為600×60×60,時間采樣間隔為4ms,空間采樣間隔為1m,子波主頻為30Hz。圖9a為三維合成數據中某條測線的地震記錄,圖9b為三維合成數據加入-6dB高斯白噪聲后該測線的地震記錄(對應圖9a同一條測線)。采用不同長度的一維前向和后向預測濾波器沿兩個方向各做一次一維預測濾波,其濾波結果分別如圖9c、圖9e 所示,數據殘差剖面中均出現了強同相軸(圖9d、圖9f),說明這種濾波方法破壞了信號結構,壓制噪聲的能力差。分別采用不同大小的數據窗和二維Wiener中心濾波算子進行預測濾波,結果如圖10 所示?？梢钥闯?當數據窗(100×10×10)和Wiener濾波算子(3×3)較小時,濾波效果較差;當數據窗(100×30×30)和Wiener濾波算子(5×5)較大時,其數據殘差剖面包含的有效信號較小且信噪比較大,濾波效果較好。

圖11中合成的三維合成單炮炮集數據大小為1500×120×120,時間采樣間隔為4ms,空間采樣間隔為10m,子波主頻為30Hz。圖11a為三維合成炮集數據中的某條測線的地震記錄,圖11b為加入-10dB高斯白噪聲后三維合成炮集數據的地震記錄(對應圖11a同一條測線)。采用二維Wiener中心濾波方法進行預測濾波,圖11c顯示了數據時窗大小為300×30×30,濾波算子大小為3×3時的二維Wiener中心濾波后的地震記錄,圖11d為采用二維Wiener

圖9 三維合成數據及一維前向和后向預測濾波結果a 三維合成數據中某條測線的地震記錄; b 加入-6dB高斯白噪聲后三維合成數據中某條測線的地震記錄(對應圖9a同一條測線); c 數據窗大小為100×10×10,算子長度均為3時前向預測和后向預測濾波后的結果; d 圖9c與圖9a的數據殘差剖面(SNR為15.69dB); e 數據窗大小為100×10×10,算子長度均為5時前向預測和后向預測濾波后的地震記錄; f 圖9e與圖9a的數據殘差剖面(SNR為12.98dB)

中心濾波方法得到的圖11c與圖11a的數據殘差剖面,圖中無明顯的連續的同相軸,并且噪聲幅值小?？梢钥闯鯳iener中心濾波方法在壓制噪聲的同時較好地保留了有效信號。

4.2 三維實際地震資料

圖12a為某地區實際采集的四維數據(疊前CDP道集)中某個CDP道集(三維數據體),其數據大小為1501×120×34。當數據窗大小為200×30×20,濾波算子為3×3時,采用二維Wiener中心濾波方法對數據進行預測濾波,得到的CDP道集和殘差剖面分別如圖12b和圖12c所示?？梢钥闯龆SWiener中心濾波方法能夠壓制隨機噪聲和部分相干噪聲,在振幅存在局部變化時(圖12a和圖12b中的紅圈)能較好地保留信號。

圖11 三維合成炮集數據二維Wiener中心濾波結果a 三維合成炮集數據中某條測線的地震記錄; b 加入-10dB高斯白噪聲后三維合成炮集數據中某條測線的地震記錄(對應圖11a同一條測線); c 數據窗大小為300×30×30,濾波算子大小為3×3時二維Wiener中心濾波后的地震記錄; d 圖11c與圖11a的數據殘差剖面(SNR為12.24dB)

圖12 實際數據二維Wiener中心濾波結果a 野外采集得到的某個CDP道集; b 數據窗大小為200×30×20,濾波算子大小為3×3時二維Wiener中心濾波后的CDP道集; c 圖12b與圖12a的數據殘差剖面

某地區的三維實際數據疊后剖面如圖13a所示,數據大小為1501×567×10。當數據窗大小為300×20×12,濾波算子為3×3時的采用二維Wiener中心濾波方法對數據進行預測濾波,得到的疊后剖面和殘差剖面分別如圖13b和圖13c所示?？梢钥闯?二維Wiener中心濾波方法能夠壓制中、深層雜亂的噪聲,并且很好地保留淺、中層連續的同相軸以及斷層信息。

圖13 實際數據二維Wiener中心濾波結果a 三維野外實際數據疊后剖面; b 數據窗大小為300×20×12,濾波算子大小為3×3時二維Wiener中心濾波后的疊后剖面; c 圖13b與圖13a的數據殘差剖面

5 討論與結論

“兩寬一高”地震勘探針對的是小尺度的精細油藏和復雜油藏,相應的成像方法也逐漸進入以FWI+LS_RTM為標志的反演成像階段。噪聲對成像的影響需要被重新審視。偏移成像中的隨機噪聲可通過高覆蓋的成像疊加被壓制,對成像結果影響較小;而反演成像過程中噪聲則會嚴重影響反演成像方法的收斂性及反演結果的精度(分辨率)。另外,針對當前“兩寬一高”地震數據采集得到的巨量數據體,在實際生產中實用的、高維空間中可實現的且效率和效果達到均衡的去噪方法還是相當缺乏的。

本文所討論的高維Wiener中心濾波方法,是由基于AR模型中的前向預測濾波器和后向預測濾波器發展得到的,其要求信號是線性可預測的并且噪聲滿足高斯分布的。對于非線性的信號,可采用對數據局部取窗的方式滿足局部線性可預測的假設條件;對于噪聲不滿足高斯分布假設的信號,需要修正Bayes理論框架下的信號預測理論和方法以達到更好的去噪效果。地震數據基本都是高維數據體,采用低維空間中的預測算子很難利用其高維空間信息,也無法達到滿意的去噪效果。本文通過對數據體進行Hankel矩陣排布解決了這一問題,結合Wiener中心預測濾波方法,實現了三維或者更高維數據空間中的噪聲壓制。實驗結果表明,高維Wiener中心濾波方法能夠在有效壓制噪聲的同時較好地保留信號,但存在計算效率低的問題。本文采用MPI和OpenMP并行計算可以部分地提高計算效率。但是本文方法也存在一些問題,如對于非規則數據、復雜波場(振幅變化較為強烈)以及在信噪比相對較低的情況下則很難最佳地預測信號并壓制噪聲。

到目前為止,地震數據去噪(插值)真正的問題依然是所建立的地震信號預測算子不能很好地預測實測數據中的信號,以及噪聲的概率分布是未知的。當前,信號建模的合理性、參數估計的精度、算法計算效率等方面遠未達到理想的程度,這使得去噪結果達不到反演成像的要求。“兩寬一高”的數據采集方式提供了更好的數據基礎,但(尤其在復雜近地表探區)地震數據在空間時間上的復雜變化使得當前已有的去噪技術并沒有達到令人滿意的、滿足反演成像要求的程度,未來仍需要繼續深入研究。

致謝:感謝中石油勘探開發研究院及西北分院、中海油研究院和湛江分公司、中國石油化工股份有限公司石油物探技術研究院和勝利油田分公司對波現象與智能反演成像研究組(WPI)研究工作的資助與支持。