利用格林函數的多次迭代鬼波壓制方法

宋建國 馬 安 黃 晟 劉 炯 陳飛旭

(①山東省深層油氣重點實驗室，山東青島 266580； ②中國石油大學(華東)地球科學與技術學院，山東青島 266580； ③中海石油海南分公司，海南海口 570311； ④中國石化石油勘探開發研究院，北京 100083； ⑤中國石油塔里木油田公司勘探開發研究院，新疆庫爾勒 841000)

0 引言

在海洋地震勘探中，因海水表面這一強反射界面的存在，檢波器采集的地震數據會包含反射波和鬼波[1]。鬼波的存在會嚴重限制地震勘探的頻帶寬度、降低地震資料的信噪比。因此，為提高地震資料品質，必須研發能有效壓制鬼波的適用技術。經過幾十年的研發，已涌現出多種鬼波壓制技術，這些技術主要分為兩大類。

第一類是通過改變采集方式，采集更寬頻帶地震資料以限制鬼波。這些地震采集技術包括上下纜采集、雙檢采集、斜纜采集和海底電纜采集等[2]。上下纜采集技術最早由Parrack[3]提出，隨后S?nneland等[4]將波場分離方法應用于上下纜采集地震資料，并有效地壓制了鬼波。但上下雙纜必須在同一垂直剖面內，對電纜的定位精度要求較高，且雙纜方式的施工成本也較高。隨著科技的發展，雙檢采集技術被提出并應用于海上地震勘探，利用水檢和陸檢接收波場的極性特點，將水檢與陸檢資料能量匹配合并處理以壓制鬼波[5-6]，該技術在使用前需先對陸檢資料進行標定，且合并壓制效果受陸檢資料信噪比影響。斜纜采集技術是將電纜傾斜布置于海水中，不同炮檢距的檢波器因沉放深度不同，其采集的地震資料波場也不同，即采集的波場更豐富，鬼波時差也具有多樣性[7-8]，有利于后續鬼波壓制處理。近年來，海底電纜采集技術越來越受到重視，海底電纜將檢波器布設在海底，位置相對固定，受海浪影響小，具有較大的炮檢入射角，故在海洋地震勘探中具有明顯優越性[9]。

第二類是鬼波壓制處理技術，是指在資料處理階段通過反演手段計算一次波或通過波場外推方式計算出鬼波，再通過匹配相減壓制鬼波。Fokkema等[10]是最早通過頻率—波數域波場外推法實現鬼波壓制的。Zhang[11]提出一種基于格林函數理論的鬼波壓制技術，該方法引入參考介質與擾動介質，通過選擇不同的閉合曲面和格林函數壓制經海面反射的下行波，即壓制了鬼波。賀兆全等[12]論證了雙檢合并理論及處理技術，由于水檢和陸檢的自身結構和工作原理不同，二者對上行波響應極性相同，而對鳴震、鬼波等下行波極性相反，將水檢與陸檢合并可得到較滿意的鬼波壓制效果，但陸檢資料的信噪比往往較低，這也成為制約該方法應用效果的主要因素。Soubaras[13]提出利用偏移和鏡像偏移數據進行聯合反褶積的鬼波壓制方法，即首先生成相對海水表面對稱的鏡像數據，再對數據和鏡像數據分別做疊前時間偏移，最后利用偏移后的兩個數據道集進行聯合反褶積，成功壓制了鬼波。Wang等[14]基于Bootstrap方法，在f-x-y域實現了鬼波壓制。

王芳芳等[15]詳細論述了逆散射級數法鬼波壓制原理，并實現了基于該理論的鬼波壓制。Sun等[16]在MWD方法基礎上提出EMWD(Enhanced model-based water-layer related demultiple)算法，與MWD相比該方法減少了對自適應減法的依賴性，降低了計算成本。Song等[17]提出頻率域高分辨率Radon變換鬼波壓制方法，并取得了較好的實際應用效果。王沖等[18]推導出頻率—慢度域鬼波濾波算子及自適應迭代反演求解上行波的算法。張威等[19]基于平面波傳播理論，導出τ-p域拖纜觀測的總波場與海水表面觀測的上行波波場之間的新關系式； 通過建立線性τ-p方程，利用最小平方殘差法反演得到海水表面觀測的上行波波場； 再通過延拓得到拖纜觀測的上行波波場。馬繼濤等[20]提出一種基于波場外推和閾值截斷的鬼波壓制算法，該算法對含鬼波的地震數據進行兩次波場外推，并對結果進行求和及閾值截斷處理，得到海面位置的地震數據； 然后利用該數據再次進行波場外推，分別得到鬼波和一次波，通過自適應相減得到不含鬼波的一次波波場。封強等[21]推導出包含海面反射系數和鬼波延遲時間的鏡像記錄生成公式和聯合反褶積去鬼波公式，并基于非高斯性最大化，以負熵為非高斯性度量，采用滑動時空數據窗口解決鬼波參數變化問題，獲得最優化的海面反射系數和鬼波延遲時間，進而通過聯合反褶積壓制鬼波。畢麗飛等[22]通過簡化逆散射級數公式預測層間多次波，引入二維卷積盲分離方法實現層間多次波的自適應匹配相減，模型和實際資料測試均證明所提方法能在壓制層間多次波的同時，更好地保護一次波。李鐘曉等[23]引入偽地震數據算法，實現了一次波與預測多次波的分離。李繼偉等[24]提出一種基于自適應相減濾波與Curvelet變換組合的壓制面波方法。

基于格林函數理論的鬼波預測方法不需已知地層速度信息，就能準確預測鬼波； 同時格林函數的構建靈活、便捷，適用于各種觀測系統[25-26]。本文從克希霍夫衍射理論出發，根據地震波傳播路徑構建鬼波傳播的格林函數，將其與原始地震記錄褶積從而預測鬼波，再結合曲波域匹配相減技術實現對鬼波的壓制。但格林函數與原始地震數據褶積時，一次波會預測出鬼波(一階多次波)，低階多次波預測出高階多次波，使預測處理后的波場更復雜； 在曲波域匹配相減采用閾值處理方法時，很難設置合適閾值，即會削減一次波能量，造成一次波損傷和鬼波壓制不徹底等問題。本文提出的多次迭代格林函數鬼波壓制方法是先將原始波場與預測鬼波進行匹配相減得到一次波波場，利用該一次波波場與格林函數進行褶積計算新的鬼波波場，再從原始波場中將鬼波匹配相減得到更干凈的一次波波場。模型和實際資料處理結果表明，該方法能在有效保護一次波的同時較好地壓制鬼波，且顯著拓寬地震資料頻寬。

1 方法原理

1.1 基于格林函數的鬼波預測方法

基于格林函數的鬼波預測方法具體過程為：根據已知海水深度和地震波在海水的傳播速度等參數，構建鬼波傳播的格林函數； 將其與原始地震記錄褶積得到預測波場，將原始波場與預測波場進行匹配相減得到一次波波場； 將該一次波場與格林函數褶積計算新的預測波場，再從原始波場中減去預測波場得到新的一次波。通過上述迭代運算，逐漸實現對鬼波的準確預測及壓制。具體表達式如下

(1)

式中：M1表示預測波場；D表示原始地震記錄；G表示格林函數；P1表示通過一次匹配相減得到的一次波波場；Mi+1表示第i次迭代計算的預測波場； ?表示褶積運算；Ξ表示自適應相減過程。構建的格林函數為

(2)

式中：R代表海水表面反射系數；ω為角頻率；r為鬼波相對于一次波多傳播的距離；t為對應的傳播時間。但這樣構建的格林函數僅考慮了地震波傳播時能量的球面擴散和時移等運動學特征，并未考慮地震波的動力學特征，而動力學特征對于鬼波的精準預測也十分重要。

本文基于Kirchhoff衍射理論，閉合曲面SC外一個點源P0對曲面內一點P的波場值可由曲面邊界上的波場值求出。即一個閉合曲面SC對曲面內某一點P的波場值總貢獻u為

(3)

式中：A表示曲面SC上的波場振幅；rP0-Q為點源P0到曲面SC上一點Q的距離；rQ-P為點Q到P之間的距離；k=ω/vw，為波數，其中vw為海水中波速；θ為rP0-Q與rQ-P之間夾角；K(θ)表示角度因子，其值如下

(4)

在海底電纜勘探中，將電纜以上半空間看成閉合曲面SC的內部，其上半空間趨于無窮遠處的曲面上波場為零，因此可由電纜上的波場值計算出SC曲面內某一點P的波場值，從而實現波場延拓。

鬼波的形成可看成是檢波器記錄的地震波波場延拓后的克希霍夫疊加[27]。對于拖纜觀測系統，檢波器鬼波傳播路徑如圖1a所示，可見其海底一階多次波(圖1b)相對于檢波器鬼波多經過了一次海底反射，且多傳播了一段距離，這個多傳播的距離會隨電纜深度而變化。對于海底電纜觀測系統，相當于將拖纜觀測系統的電纜深度增大到海底。此時海底一階多次波(圖1d)相對于鬼波(圖1c)多傳播的路徑很短，其傳播時間可忽略不計，因此海底電纜的一階多次波與檢波器鬼波會混疊在一起，這兩種波可在處理中同時進行壓制。

圖1 海洋地震勘探地震波傳播示意圖(a)拖纜觀測系統鬼波傳播路徑； (b)拖纜觀測系統一階多次波傳播路徑；(c)海底電纜觀測系統鬼波傳播路徑； (d)海底電纜觀測系統一階多次波傳播路徑

本文以海底電纜觀測系統為例，置于海底的檢波器在接收到海底反射波的同時也會作為一個新的球面波震源，產生二次元波向外擴散，二次源波經海面反射后產生下行波場，再次被檢波器接收。按照該原理對海底電纜采集波場進行延拓，可實現對海底電纜鬼波和多次波的預測，但在延拓過程中只考慮海水表面的強反射，因此預測的只是與水層有關的多次波。經延拓后，海底地層一次反射波變為檢波器鬼波或一階多次波(圖2紅色路徑)； 海底一階多次波經延拓后變為二階多次波(圖2綠色路徑)，依此類推，低階多次波延拓后變為更高階多次波。

圖2 海底電纜波場延拓示意圖

根據地震波的傳播路徑，進一步計算經海水表面反射后最終被檢波器接收到的波場值

(5)

式中：F(t)=R·δ(t)，表示海面反射波的理論脈沖；rex為波場延拓距離。由該式可將鬼波預測格林函數改寫為新形式，即在式(2)中加入角度算子以控制不同入射角地震道的疊加權重，則格林函數變化為

(6)

在三維情況下利用格林函數方法理論預測鬼波時，需對大量數據進行疊加求和，對計算機運算性能有一定要求。為了在不影響預測效果的前提下提高計算效率，采用限制預測孔徑的方法，只使用設定孔徑內的地震數據進行鬼波預測(圖3)，這樣可在保證鬼波預測精度的同時大幅度減少運算量。

圖3 有限孔徑鬼波預測示意圖

1.2 曲波域匹配相減技術

傳統的匹配相減技術是在最小平方理論基礎上發展起來的，該方法在能量匹配過程中會干擾原始地震數據并引入新的噪聲。因此，為了更好地適應各種復雜地震數據，在盡可能消除鬼波的同時保留有效信號，本文引入曲波變換法以實現匹配相減。曲波變換由于其多尺度和多方向等特性，在圖像處理中可很好地對細節信息進行刻畫和提取，在地球物理領域得到廣泛應用[28]。另外，相比于小波變換等其他稀疏變換，曲波變換對二維分段光滑函數具有更好的非線性逼近性，同時其稀疏性也更強[29]，可取得更好的匹配相減效果。

曲波變換是將原始信號在給定的函數或矢量集上進行分解，得到原始信號在曲波域的稀疏表示。曲波變換通過基函數與信號的內積形式實現信號的稀疏表示，連續的曲波變換的系數可表示為

(7)

(8)

在得到曲波系數c(j，l，p)之后，采用閾值處理方法壓制鬼波。常見的曲波閾值處理方法有兩種，即硬閾值法和軟閾值法[30]。

硬閾值法的準則是，設定一個閾值T，將大于T的曲波系數予以保留，小于T的曲波系數賦值為零，即

(9)

軟閾值法的準則是，同樣設定一個閾值，將小于閾值的曲波系數賦值為零，大于閾值的系數減去給定閾值，即

(10)

式中： sgn(·)為階躍函數；c(j,l,p)、C(j，l，p)分別為閾值處理前、后的曲波系數。

硬閾值法和軟閾值法都有各自優點和局限性。硬閾值法在保留有效波同相軸及弱小信號方面有較好效果，在曲率較大處去噪效果較好； 軟閾值法在壓制噪聲方面優勢較明顯，在相同閾值情況下，噪聲壓制比硬閾值法更徹底，但同時也會在一定程度上損害一次波。在曲波變換壓制鬼波方法中，閾值的設定至關重要，因此本文在上述兩種方法基礎上，采用改進的相對強度閾值設定方法。

在預測的鬼波與原始地震數據匹配較好的情況下，將原始地震數據的曲波系數與預測的鬼波的曲波系數做對比，得到相對強度，即有

(11)

式中：cr(j，l，p)為原始地震數據的曲波系數；cm(j，l，p)為匹配后預測的鬼波曲波系數。相對強度的閾值法準則為： 若q(j，l，p)T，則將其視為有效波。然后再采取硬閾值法計算曲波系數，即小于閾值的系數賦值為零，大于閾值的系數被保留。

在預測鬼波過程中，由于子波變化和能量差異等原因，造成地震資料在相位和振幅等信息上存在誤差，所以可根據地震數據能量隨時間的變化用不同的時窗截取地震數據。把每個時窗內的數據進行最小二乘匹配； 再將原始數據與匹配后的鬼波數據進行曲波正變換，在曲波域采用閾值處理法壓制鬼波； 最后將處理后的數據反變換到時間域，便可得到壓制鬼波后的地震數據。

2 模型試算

為了驗證本文論述方法的有效性，使用射線追蹤法對海底電纜進行正演模擬，通過模型數據分析本文方法有效性。設計層狀模型海底深度為50m，海水速度為1500m/s，海底以下地層速度依次增大，震源深度為5m，檢波器放置于海底，共1201道，道間距為5m，震源采用主頻為30Hz的雷克子波，采樣間隔為1ms，接收5s的地震記錄。正演模擬時，對模型上邊界不做任何處理，以模擬海平面的自由界面； 對左、右、下界面使用PML吸收邊界條件，以模擬地震波在半無限介質中的傳播。數值模擬采用的速度模型如圖4所示，圖中橫縱坐標(m)可轉化為網格點數，網格單元尺寸為5m×5m。

從利用速度模型(圖4)正演模擬得到的單炮記錄(圖5a)可見，在每一個地層界面的反射波之后都跟隨一個與其波形相似的波，它們之間存在一定的時差，該“跟隨”波即是鬼波。

圖4 速度模型

圖5b為一次迭代鬼波壓制后單炮記錄，從中可見鬼波雖得到很好壓制，但同時對一次波也造成了一定損傷。為此，調整曲波域匹配相減的閾值，所得結果如圖5c所示，可見調整匹配相減參數后一次波得到了很好保留，但其鬼波尚有一定程度的殘留。于是利用本文多次迭代法對壓制后單炮記錄進行第二次迭代預測，將預測的多次波再從原始數據中匹配相減后得到圖5d所示結果。第二次迭代匹配相減所用閾值與第一次圖5b中設定閾值完全相同。對比這兩次迭代結果，可見在圖5d中不僅一次波得到充分保留，而且在圖5c中殘留的鬼波也得到了較徹底的壓制。

圖5 模擬地震記錄的一次、二次迭代結果(a)正演炮集記錄； (b)一次迭代壓制鬼波后記錄(匹配相減損害一次波)；(c)一次迭代調整匹配閾值后記錄(匹配相減殘留鬼波)； (d)二次迭代鬼波壓制后炮集記錄

通過單道數據分析展示鬼波壓制效果。圖6a和圖6b分別為第800和第1000道的單道曲線鬼波壓制前(紅線)、后(藍線)的結果，對比可見在鬼波壓制后緊跟在一次反射波之后的鬼波已被消除。

圖6 第800 道(a)和第1000道(b)鬼波壓制前、后單道對比

進一步對鬼波壓制前、后單炮記錄做頻譜分析，驗證本文方法對陷頻點的補償作用。陷頻點公式為

(12)

式中：fnotch為陷頻點處頻率；n代表第n個陷頻點；h為電纜沉放深度； 模型數據的電纜深度為50m，可知陷頻點的頻率分別為15、30、45Hz。

圖7a和圖7b分別為鬼波壓制前、后地震數據的f-k譜，可見鬼波壓制后，陷波頻率得到較好補償。對鬼波壓制前、后地震數據的頻譜(圖8)進行分析，可見鬼波壓制前(紅線)位于15、30、45Hz等處的陷頻點，在鬼波壓制后(藍線)都得到有效補償。

圖7 鬼波壓制前(a)、后(b)的f-k譜分析

圖8 鬼波壓制前、后模擬數據頻譜分析

3 實例分析

M工區水深約為18m，使用海底電纜進行采集，共部署5條測線，每炮650道，道間距為50m，采樣點為3500個，采樣間隔為2ms，震源深度為5m。對比鬼波壓制前(圖9a)、后(圖9b)的炮集記錄，可見處理后的單炮記錄上鬼波得到較好壓制，消除了一些假反射軸，提高了資料信噪比。

圖9 原始(a)和用本文方法壓制鬼波后(b)的單炮記錄

分析從該實際炮集記錄得到的頻譜(圖10)，可見鬼波壓制后陷頻點振幅得到很好補償，地震資料的高頻成分增強，頻帶拓寬。

圖10 鬼波壓制前、后實際地震數據頻譜分析

觀察鬼波壓制前(圖11a)、后(圖11b)該工區實際疊加剖面，可見鬼波壓制后消除了(原始疊加剖面中的)一些虛假同相軸(紅色箭頭指示處)，地震資料的信噪比得到提高。

圖11 鬼波壓制前(a)、后(b)實際疊加剖面

圖12為另一實際工區鬼波壓制前、后的對比疊加剖面，可見鬼波壓制前的一些虛假同相軸受到壓制(圖12b)，進一步證明了本文所提方法對鬼波的有效壓制作用。

圖12 鬼波壓制前(a)、后(b)局部疊加剖面對比

4 結論

鬼波壓制是海洋寬頻地震勘探核心問題，采用基于格林函數的多次迭代鬼波預測與匹配相減鬼波壓制方法，對模擬和實際海底電纜數據進行處理，有效提高了數據的信噪比，得出如下結論。

(1)根據鬼波傳播路徑構建鬼波的格林函數，將其與原始地震數據褶積計算出鬼波，通過匹配相減得到一次波后，將該一次波再與格林函數褶積得到新的鬼波波場，通過如此迭代可逐步提高鬼波的預測精度。

(2)曲波域匹配相減是采用相對強度的閾值處理實現鬼波壓制，但在壓制鬼波的同時，務必兼顧最大程度地保留一次波有效信息。

(3)本文方法屬于模型驅動方法，需用準確的檢波器深度和海水速度計算準確的波場延拓距離，從而實現高精度的鬼波預測，且適用于拖纜和海底電纜等多種觀測系統。

(4)實際海上地震勘探中，海水界面反射系數常常不是嚴格的-1.0，因此需考慮粗糙海水表面對波場延拓的影響。海底沉積物的軟硬對地震波反射振幅也有一定影響。下一步擬將這兩個影響因素添加到算法中，從而進一步改進、完善算法。