編碼與解碼框架下的局部平面波域淺水多次波壓制方法

徐 鵬,劉亞輝,王華忠

(1.波現象與智能反演成像研究組(WPI),同濟大學海洋與地球科學學院,上海200092;2.中國石油大學(華東)地球科學與技術學院,山東青島266580)

隨著地震采集技術的不斷發展,“兩寬一高”地震勘探日漸常態化,但是地震數據的品質與成像結果的分辨率仍然受到多種因素的制約,多次波是其中重要因素之一。當前勘探地震學的核心問題是地震波反演成像,其目標在于利用低波數的光滑速度場估計高波數散射體。現有成像方法多數假設觀測數據中的反射波僅在地下經歷一次反射,而地震波在真實介質中的傳播通常包含多個波阻抗界面的反射,即包含多次波,常規針對一次波的成像結果會含有假象,嚴重影響后續的數據分析與解釋。

對于多次波的處理可以分為“壓制”和“利用”兩類。“壓制”類是將數據中的多次波視為噪聲去除后,利用常規的一次波成像方法處理;“利用”類則將多次波也視為有效信號,使用更復雜的成像方法對一次波、多次波進行聯合成像。后者能更多地挖掘出數據中的地下介質信息,是更理想的成像方式。例如劉伊克等[1]對南海深水區域含多次波的數據進行偏移成像,成功彌補了目標區域的照明不足;王永強[2]利用最小二乘逆時偏移(LS-RTM)對一次波-多次波進行聯合成像,驗證了一次波-多次波聯合成像方法的可行性。

但是,就目前研究現狀而言,多次波成像方法適用范圍有限,因為它對偏移速度場的要求過于嚴格,同時為了消除不同階次的反射波之間的錯誤相關值,引入最小二乘(LS)迭代流程,使得總體計算量增加。因此,先去除多次波,再利用成熟的常規算法進行成像是目前更為合理的選擇。

勘探地震技術發展至今,已有多種多次波壓制方法在實際生產中得到應用。根據其基本原理,可將這些方法分為兩類:第一類是濾波方法,基于多次波與一次波在道集中的形態特性差異,將多次波剔除;第二類是褶積方法,基于多次波的物理成因,由波動理論建立一次波與多次波的高維褶積模型,預測出多次波的位置、形態,實現對多次波的精確壓制。第一類方法主要利用一次波與多次波(在變換域)的可分離性壓制多次波,具有成本低、計算量小等特點,實際應用廣泛,代表方法有最優疊加、特征值譜濾波、f-k域濾波等。然而,由于多次波與一次波并非完全可分,因此壓制效果差強人意。第二類方法采用“預測-減去”的策略,預測階段主要考慮多次波的運動學信息,減去階段則主要用動力學信息進行匹配。代表方法有經典預測反褶積、純數據驅動的自由表面多次波消除(SRME)[3-6]、基于模型的淺水多次波壓制(簡稱SWD)[7]等。該方法將多次波壓制分解為多次波模型預測與自適應相減兩步,降低了實現難度,但是引入了額外的中間數據,降低了自動化程度[8-9]。

淺海環境為多次波壓制問題引入了額外的難點。淺海環境主要的多次波是水體相關多次波,由于海底-海面路徑短,且海面、海底都是強反射界面,因此該類多次波的特點是周期短、能量強,對有效信號干擾極強[10-11]。然而,由于現有的觀測系統大多是針對中深層目標體設計的,用于淺海底的反射波觀測時采樣間隔偏大,同時更易造成近偏移距數據的缺失,故常規的多次波壓制方法不能直接用于淺水數據的多次波壓制。例如,經典的SRME作為一種數據驅動的高維褶積類方法,對數據采樣間距、完備性有較高要求,因而無法在淺水數據中取得令人滿意的效果[12-15]。SWD方法引入了海底形態模型,對海底反射的刻畫更準確,但其本質仍是Rayleigh積分的數值實現,故對積分孔徑、積分網格都有較高要求,在近偏移距數據缺失時仍有局限性。因此,需要發展針對淺海環境的多次波壓制技術。

本文針對淺海環境中多次波問題的特點,引入信號領域的“編碼-解碼”框架,設計一種針對淺水多次波的“預測-減去”類壓制方法,借助稀疏τ-p變換克服不規則空間采樣問題,再利用射線理論高效生成多次波模型,最終在局部平面波域完成觀測數據與多次波模型的匹配相減。

1 水體多次波問題描述

水體中存在海水-空氣、海水-海底兩個強反射界面,因此海上地震數據中多次波往往更為發育。海上勘探中多次波種類主要有:鬼波、水體相關多次波、其它自由表面相關多次波以及層間多次波。前三者成因與海上觀測系統緊密相關,層間多次波基本與陸上勘探相同。

鬼波可分為源鬼波、檢鬼波和源檢鬼波。震源激發后,地震波上行至水面發生反射繼而下行,稍遲于震源直接激發的下行波傳播,即源鬼波。地震波與地下介質相互作用后上行至海面反射一次后,再傳至檢波點,即檢鬼波。若同時符合源鬼波、檢鬼波的路徑特點,則稱為源檢鬼波。

水體相關多次波是自由表面相關多次波的一種特例,以海底向上反射開始或/和結束,可分為源端、檢端和源檢兩端三種。源端水體相關多次波是指地震波自震源激發后,先在海底-海面間震蕩一次,再自海面向深層介質傳播的波現象。檢端水體相關多次波是指地震波穿過地層介質后,最終被檢波器接收之前,在海面-海底間反射一次所形成的波現象。同時滿足源端水體相關多次波、檢端水體相關多次波路徑特點的波現象被稱為源檢兩端水體相關多次波。另外,僅在水層中震蕩傳播的波現象被稱為純水體鳴震[16-17]。

除鬼波、水體相關多次波之外的自由表面多次波統稱為其它自由表面相關多次波。此類多次波含義較廣,例如在海底之下波阻抗界面與水面之間的多次反射、多次折射等。

在海上勘探中(尤其是淺海探區),首先必須面對、同時也是干擾最強的多次波是水體相關多次波。水體相關多次波的存在,掩蓋了深層反射信號,影響了微斷層、鹽丘下部構造的成像[12,18]。

2 正問題:用編碼觀點看待多次波形成正過程

從各種多次波成因來看,多次波實質上是一次波的延拓結果。因此可從多個角度對多次波進行表述,如將多次波表述為波動方程的積分解、鏡像虛震源激發的波場等。本文從信號分析角度出發,將觀測數據視為輸入信號,將地下介質整體視為線性濾波器(即“地層濾波系統”),多次波視為輸出信號。即將多次波視為一次波的編碼結果,而編碼過程是一次波在地下介質中震蕩傳播的過程,或者說取“地層濾波系統”的脈沖響應作為編碼算子。將多次波形成的正過程表述為編碼過程,則多次波的消除過程可理解為解碼過程。引入信號分析領域的“編碼-解碼”框架使得輸入信號、編碼算子、輸出信號等概念更加直觀,簡化整個多次波形成過程的描述,有利于多次波壓制算法的設計與實現。

圖1a為真實地下介質(含海水和空氣)示意圖,為方便討論,引入參考介質,如圖1b所示,即去除海面的自由邊界并在原海面及以上空間填充海水的假想介質。基于圖1a介質的地震波傳播情況如圖1c所示,可見由于海水自由表面的存在,一次反射波又被向下反射,繼續與地下介質發生作用,進而產生一階多次波;一階多次波繼續產生二階多次波,直至無窮階多次波。由于反射系數絕對值小于1,且地震波傳播會發生衰減,高階多次波可以忽略不計。

在參考介質中,震源直接激發不會產生自由表面多次波。真實觀測數據中的自由表面多次波可看作將真實觀測數據作為邊值條件加載在深度z=0平面上所引起的波場(圖1d)。可用經典的波動方程積分解解釋:以所有上行波在海面的反射R0U(x′;xs,ω)為邊值條件,以海面+半徑無窮大的下半球為封閉面,最終在檢波點引起的波場。即:

(1)

式中:M(xr;xs,ω)表示自由表面相關多次波;U(x′;xs,ω)表示觀測數據(含一次波、多次波);Gref(x′;xr,ω)表示參考介質中的格林函數;ω表示頻率;xs,xr,x′分別表示震源點、檢波點以及積分微元的坐標;R0表示海面反射系數;積分曲面S取為海平面;n表示積分微元處的單位法向量。公式(1)即水體相關多次波在“編碼-解碼”框架下的表達式。

等價地,可以用向量形式描述上述過程:

式中,輸入信號為所有上行波,編碼算子為海水層的脈沖響應,輸出信號為檢波點觀測到的波場。應該指出,引入編碼算子A是為了簡潔表達(1)式的線性關系,相當于離散化表示的積分核。后續計算中不需要將該算子顯式表達為矩陣,僅需要實現A的矩陣向量乘。

地震數據中的一次波及各階多次波表達式可以拆分,得到更直觀的編碼預測關系:以一次波在海面的反射作為邊值條件,在參考介質中傳播得到的波場即為一階檢端多次波;以一階檢端多次波在海面的反射作為邊值條件,在參考介質中傳播得到的波場即為二階檢端多次波;以此類推,以k階檢端多次波在海面的反射作為邊值條件,在參考介質中傳播得到的波場即為k+1階檢端多次波。

圖1 地下介質及地震波傳播a 真實地下介質; b 參考介質; c 地震波在真實地下介質中的傳播; d 在參考介質中借助面源描述地震波傳播

(4)

式中,Mi表示第i階自由表面相關多次波,P表示一次反射波。總之,算子A的作用是海水層對海面處上行波的編碼“升階”。

分別累加各式等號兩端,得:

(5)

式中,Mtotal為觀測到的總體多次波(各階多次波之和),右側為觀測數據。

同樣地,若將各階多次波都寫為關于一次波P的函數:

(6)

可得總體多次波與一次波的關系:

(7)

其中,總體多次波不是一個直接觀測的量。為了建立一次波P與直接觀測結果Uobs之間的關系,在(7)式左右兩端都加上P,因此得到描述一次波與含多次波的實際觀測數據之間的線性編碼關系:

(8)

其中,I表示單位矩陣。

3 反問題:用解碼的觀點看待多次波壓制

(9)

其中,L為迭代次數。

令:

(10)

將公式(10)代入公式(9),得:

(11)

(12)

公式(12)即為水體相關多次波的求解公式,它與SRME的“預測+減去”公式有類似之處。相對SRME使用觀測數據來代替(1)式的Rayleigh積分,本文使用了特征波域的射線束傳播算子更為精確地實現編碼算子A,完成了(12)式中的A(ω)Uobs(ω)計算,即總體多次波的預測。

公式(12)描述了理想情況下的多次波壓制或一次波估計的數學模型。但事實上,由于對水體介質各屬性(如海面和海底的反射系數分布、水體速度模型、海底起伏形態)不能進行精確描述,實際數值實現的A算子都是近似的。

多次波壓制過程具體可描述為:

1) 利用公式(13)計算或預測多次波模型。

(13)

2) 利用公式(14)估計一個整形濾波器f,使得輸出的一次波能量最小。

(14)

3) 采用公式(15)輸出一次波。

(15)

借鑒SRME方法,公式(12)的減法是自適應減法,且在局部平面波域進行。相對于空間域減法,局部平面波域自適應減法能更好地區分多次波與一次波的斜率差別,也更能容忍兩者的交叉,使得輸出的去噪結果質量更好。

4 在平面波域實現多次波的編碼預測與一次波的解碼

公式(1)是將海面每個微元作為參考介質中激發的偶極震源,積分計算所有微元在檢波點處所產生的波場。作為波動方程的積分解,理想情況下公式(1)在數學上是精確成立的,并有明確的物理意義。然而在實際應用中,有限的積分孔徑、不規則的面元分布、近偏移距數據的缺失等問題,使得公式(1)的精確性下降。編碼算子A是公式(1)中積分所描述的線性關系的緊湊表達。在局部平面波域實現編碼算子A可以緩解上述問題。在每個局部區域,作為輸入信號的入射波上行波場都可分解為不同方向的平面波[19-20]。局部區域激發的二次波場的傳播過程可以看作一組不同方向的平面波傳播過程。該方法在地震成像領域應用廣泛,如馮波[21]通過局部平面波分解解決波形反演中的“周期跳躍”問題;劉少勇[22]通過局部平面波分解解決常規偏移中的畫弧問題;孫維薔[23]將局部平面波分解用于鬼波與自由表面多次波的統一壓制。

本文將局部平面波分解用于描述多次波的形成過程。單個平面波分量的傳播過程如圖2a所示,相比傳統預測方式(圖2b),該方法考慮了各個平面波分量的方向,同時簡化了地震波在水體中的傳播過程。局部平面波域水體相關多次波預測方法的實現步驟如下:

1) 將計算區域的海面分解為若干局部窗。

2) 通過稀疏求解局部窗τ-p系數,將局部窗內的入射波上行波場分解為不同方向的平面波。

3) 從局部窗中心向檢波面進行反射射線追蹤(從窗中心出發,在海底反射,終點在檢波面,見圖2a),并記錄各條射線出射方向、走時(T)。

4) 將局部平面波當作有方向、有寬度的射線來傳播,各個τ-p系數按照與其方向對應的射線的走時計算延遲,然后“放置”到射線路徑終點對應小鄰域內檢波點的數據上(需考慮相對于各檢波點的時差)。

5) 對所有局部窗重復步驟2)～4),完成所有計算。

預測過程主要利用運動學信息構建了多次波模型,其波形、振幅等動力學信息與真實數據仍有差異,因此精細的“匹配相減”是不可或缺的措施。得益于本文方法中的“平面波分解”步驟,可以將“匹配相減”過程放在局部平面波域,從而能額外考慮方向信息,使得匹配過程更為精細。

圖2 觀測數據作為二次源對于多次波的不同貢獻方式a 本文方法中的局部平面波源方式; b 積分解形式中的偶極點源方式

5 實際應用

為了驗證本文方法的有效性及實用性,設計、實施了淺海環境的OBC數據實驗和拖纜數據實驗。

5.1 淺海OBC數據

淺海OBC數據共含有8條檢波線,每條檢波線有240個檢波器。炮點、檢波點間距為25m,炮線間距250m,檢波線間距400m,觀測系統如圖3所示。采集環境平均水深約37m,水體相關多次波周期約為50ms(該數據中單個震相的波形時長約40ms),嚴重干擾了有效信號。

圖3 OBC數據觀測系統

采用以下步驟對多次波進行了壓制:

1) 通過“近道自相關”提取工區網格內的水深,并插值為光滑曲面作為海底模型。該海底深度模型是深度關于坐標的二階光滑函數。

2) 以單炮道集為獨立單位進行多次波預測。將單炮數據按檢波點坐標劃分為若干局部空間窗口,并利用稀疏τ-p變換進行平面波分解。

3) 在每個空間窗內,將所有平面波分量“傳播”到檢波陣列中。即將每個平面波分量當作一束有寬度的射線,按照其起飛角(由參數p計算),沿“海底-海面”反射路徑,貢獻到檢波陣列(圖2a),最終得到該單炮道集的多次波預測結果。

4) 進行“局部平面波域匹配相減”,減去多次波,輸出結果即為一次波。

考慮各個單炮以及各個空間窗之間的獨立性,本文在各個單炮道集間采用MPI多進程并行計算,每炮內各空間窗之間采用OpenMP多線程并行的計算策略,從而提高計算效率。

本文方法中的稀疏估計平面波系數方法可以利用有限樣點還原出局部波前面,從而緩解采樣不足問題,因此比常規方法具有明顯優勢。其中局部平面波傳播在三維坐標系下進行,能夠自動適應實際三維OBC觀測系統。OBC數據實驗結果表明,本文方法能夠準確地預測并消除水體相關多次波。圖4和圖5 對比了多次波壓制處理前后OBC炮集和OBC數據體,圖6為原始OBC數據及其平面波系數,圖7為多次波壓制處理后OBC數據及其平面波系數,可以看出多次波已被成功分離。圖8對比了OBC炮集多次波壓制處理前后自相關結果,可見數據自相關中的周期性成分已被成功消除。

圖5 原始OBC數據體(a)與多次波壓制處理后OBC數據體(b)對比

圖6 原始OBC數據(a)及其局部平面波系數(b)

圖7 多次波壓制處理后OBC數據(a)及其局部平面波系數(b)

圖8 OBC炮集多次波壓制處理前(a)、后(b)自相關結果對比

5.2 淺海拖纜數據

淺海拖纜數據共1149炮,炮間距25m,每炮4條拖纜,檢波器間距12.5m,拖纜間距100m。采集環境水深90m,多次波周期約120ms。該數據中單個震相的波形時長約40ms,水體相關多次波的存在使得有效數據同相軸下方多出若干“虛影”。相比OBC數據,拖纜數據的觀測系統更具挑戰性:近偏移距數據缺失、拖纜發生漂移甚至彎曲(圖9)。得益于稀疏τ-p反演,本文方法在面對非規則網格時,能根據采樣點真實坐標將數據分解為不同方向的平面波,解決了觀測系統帶來的問題。

拖纜數據多次波壓制的步驟與OBC數據實驗相同。圖10對比了拖纜數據多次波壓制處理前后的疊加剖面,可以看到多次波的影響已被成功消除。圖11對比了拖纜數據多次波壓制處理前后自相關結果,可以看出零延遲之外的旁瓣已被消除,這意味著數據中的周期性成分已被成功壓制。

圖9 拖纜數據中炮檢點分布

圖10 拖纜數據多次波壓制處理前(a)、后(b)疊加剖面對比

圖11 拖纜數據多次波壓制處理前(a)、后(b)自相關結果對比

6 討論與結論

本文分析了淺海環境中多次波問題的特點,引入信號分析領域的“編碼-解碼”框架,提出了針對性的多次波壓制方法。將多次波的產生與壓制用一個線性系統及其逆進行描述,即將上行波在水層震蕩產生多次波的過程視為編碼過程,將多次波壓制過程視為解碼過程,以精簡的方式從線性系統求解角度給出了壓制算法。該方法將水體相關多次波的預測轉移到局部平面波域,利用稀疏τ-p反演緩解了采樣不規則、近偏移距數據缺失造成的影響。利用淺海OBC數據和拖纜數據驗證了方法的有效性和高效性,展示了方法在淺海環境地震數據處理中的應用前景。

本文方法相對傳統方法具有如下優勢:①稀疏求解τ-p系數,能在非均勻網格下得到各個窗口的平面波分解,克服了非均勻采樣問題,緩解了近道缺失的影響;②所需的海底模型為光滑曲面,不需要離散建模,因此能應對崎嶇海底,且不會引入模型的鋸齒和波場的數值頻散現象;③傳播過程化簡為從窗中心到海面的反射射線追蹤更高效,由于射線追蹤計算次數等于局部窗口數目,遠小于檢波點數目,故總體效率也高于MWD(model-based water-layer demultiple)算法中從檢波點到檢波點的射線追蹤;④各個單炮道集以及各個局部窗的計算過程相互獨立,具有良好的可并行性;⑤將匹配相減放在局部平面波域進行,能夠區分一次波、多次波的不同視速度,提高多次波壓制精度。

但是,本文方法在多次波預測過程中將每個平面波分量當作一束有寬度的射線(即胖射線)傳播,對于射線寬度的選擇尚缺理論依據,需要進行手動調試。射線過寬時,所構建的多次波模型的頻帶將偏低;射線過窄時,構建的多次波模型會出現一些空缺點。可能的解決方法是根據“菲涅爾半徑”確定射線寬度,需要開展進一步的探索研究。

展望未來多次波相關技術的發展,“壓制”與“利用”應當并重。勘探地震學的重要方法如層析、FWI、成像等,目前仍依賴于一次波假設。在實際應用中,多次波壓制的效率與效果,將與最終處理、解釋結果緊密相關。因此,發展有效的多次波壓制技術具有重要意義。同時,多次波成像應予以重視,一種能夠應對多次波的成像方法將大大簡化地震數據的預處理,并能更多地挖掘出數據中的地下介質信息,最終呈現出更準確的成像結果。故對待多次波,應當立足“壓制”,展望“利用”。

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