頻率—慢度域延拓關鍵參數優化的鬼波壓制方法

青杰,顧漢明*,王建花

(1.中國地質大學(武漢)地球物理與空間信息學院，湖北武漢430074； 2.中海油研究總院有限責任公司，北京100028；3.海洋油氣勘探國家工程研究中心，北京100028)

0 引言

海上勘探中，拖纜通常沉放于水面以下，受水面強反射的影響，海洋地震資料中地層反射波后會伴隨反相位的鬼波。鬼波的陷波效應使得地震記錄有效頻帶變窄，故壓制鬼波是拓寬地震資料頻帶的重要途徑。

為了得到分辨率和信噪比更高的寬頻地震資料，上/下纜采集[1-3]、雙檢接收技術[4-5]、上/下源激發技術[6-7]等相繼提出。此外，斜纜采集也被推出并應用，其通過變化檢波器深度的采集方式使陷波特性多樣化，壓制鬼波后可較好地拓展地震資料的頻帶。Soubaras等[8-9]提出原始數據結合鏡像數據聯合反褶積壓制鬼波；Sablon等[10]將多道聯合反褶積技術應用于同步多級震源與二、三維斜纜數據；Wang等[11]提出了傾斜拖纜“bootstrap”處理方法，用于處理τ-p域中的斜纜數據；Song等[12]在高分辨率Radon變換基礎上進行傾斜拖纜數據處理，拓寬了地震頻帶；張振波等[13]在珠江口盆地使用斜纜采集、處理方法獲得了高品質的海上地震數據；王建花等[14]通過波場延拓生成鏡像數據共軛迭代反演獲得去鬼波地震數據；張威等[15]利用最小平方殘差法反演求解海表面上行波，并延拓得到拖纜處的上行波波場；馬繼濤等[16]將基于波場外推和閾值截斷的方法應用于鬼波壓制；宋建國等[17]提出利用格林函數多次迭代壓制鬼波的方法；張寶金等[18]研究了倍頻視角下的寬頻地震子波優化方法，用于鬼波壓制后寬頻地震數據子波整形；王兆旗等[19]對不同時刻的靜態粗糙海面組合得到了動態海面，用以測試鬼波壓制效果。在復雜采集環境中，假設條件下的鬼波參數與實際參數相比誤差較大，為此，Grion等[20]提出了以相移法延拓波場來壓制鬼波的算法，對澳大利亞近海起伏海面上的數據有著良好的應用效果；Matta等[21]提出了自適應去鬼波方法；許自強等[22]研究了抑制鬼波的優化聯合反褶積算法；王沖等[23-24]分別在頻率—波數域與頻率—慢度域對斜纜進行迭代反演壓制鬼波，擴展了地震頻帶；張威等[25]在平纜數據上對反演得到的海面波場正、反向延拓到拖纜處獲得拖纜深度，消除了起伏海面影響；封強等[26]以負熵度量去鬼波結果的非高斯性來搜索得到更精確的海面反射系數和鬼波延遲時間，繼而通過反褶積壓制鬼波；王永強等[27]基于編碼/解碼思想，在Bayes反演框架下進行了鬼波預測與壓制；陳勝紅等[28]以自舉法τ-p域稀疏反演鬼波壓制等多種方法為核心，形成了寬頻準三維處理集成技術。

鬼波濾波方法需考慮鬼波與一次波的延遲時間和振幅差異系數，但崎嶇海底、起伏海面及目的層深度等因素會使鬼波參數受到一定的擾動，經常規方法處理后記錄中會存在噪聲。本文基于頻率—慢度域鬼波壓制方法構建了波場延拓算子，用原始記錄分別加上、減去延拓一次的記錄，以兩個新記錄乘積的絕對值之和作為度量參數，交替迭代尋找最優平均鬼波延遲時間和振幅差異系數，進而去除鬼波。另外，為降低鬼波參數時變的影響，利用滑動時間窗口將地震記錄離散化，使得鬼波參數在窗口內近似相等，再分別對窗口內數據尋找最優鬼波參數以消除鬼波。斜纜合成數據和實際資料的試算結果表明，該方法能較好地消除鬼波，達到拓寬地震記錄頻帶的目的。

1 頻率—慢度域鬼波壓制原理

在拖纜沉放至某一深度時，由于海水面起伏、海底崎嶇以及復雜構造等因素，拖纜鬼波與一次波的振幅差異系數不再為-1，這里用R表示該差異系數。對于同一水平慢度pj的平面波，由圖1中鬼波波場ug與一次波波場up的幾何關系可得

圖1 檢波器接收記錄示意圖

(1)

式中：τ為波前面的截距時間；Δτ為鬼波與一次波的延遲時間；xi、zi分別為檢波器炮檢距離和深度，其中i=1,2,…,N，且N是記錄的道數；θj為水平慢度pj對應平面波的出射角度,j=1,2,…,M，其中M為慢度個數；VW是海水波場傳播速度。

檢波器的波場u是一次波與鬼波相加的結果，即

u(pj,τ)=up(pj,τ)+uG(pj,τ)

=up(pj,τ)+Rup(pj,τ-Δτ)

(2)

該接收點總波場d(xj,t)為M個水平慢度之和，即

(3)

式中τ=t-pjxi，其中t為時間。對式(3)做關于時間t的傅里葉變換，得

(4)

式中ω為頻率。對于N道數據而言，每一個頻率ω有

D(x,ω)=G(x,p,ω)U(p,ω)

(5)

式中：D(x,ω)為有N個元素的列向量；U(p,ω)是有M個元素的列向量，即

D(x,ω)=[D(x1,ω),D(x2,ω),…,D(xN,ω)]T

(6)

U(p,ω)=[U(p1,ω),U(p2,ω),…,U(pM,ω)]T

(7)

G(x,p,ω)為N×M階矩陣，即

G(xi,pj,ω)=(1+Re-iω Δ τ)e-iωpjxi

(8)

U=(GTG+μI)-1GTD

(9)

式中：μ為阻尼系數；I為單位矩陣。

2 頻率—慢度域延拓原理

將式(8)中G(xi,pj,ω)分解可得

G(xi,pj,ω)=Gp(xi,pj,ω)-Gg(xi,pj,ω)

(10)

式中：Gp(xi,pj,ω)=e-iωpjxi；Gg(xi,pj,ω)=-Re-iω Δ τe-iωpjxi。由式(5)可得

D(x,ω)=[Gp(x,p,ω)-Gg(x,p,ω)]U(p,ω)

(11)

(12)

(13)

(14)

以此延拓類推，則有

(15)

3 鬼波參數擾動分析

由式(8)可見，鬼波算子與鬼波參數(振幅差異系數及延遲時間)有關。實際資料中，鬼波參數會因采集環境不同而有所變化，海底的復雜地形會使鬼波與一次波的路徑相差較大，延遲時間因此不同；此外，由于檢波器中接收的鬼波和上行波不是從同一個點反射，兩者的振幅差異系數也會偏離假定的反射系數。海面是隨時間起伏的，拖纜深度隨之產生變化，鬼波延遲時間也隨空間和時間變化；反射層深度會使得延遲時間時變，也會影響鬼波和上行波的振幅差異系數。當然影響鬼波延遲時間的因素還有海水中波場傳播速度等。

在此以延遲時間的擾動簡述崎嶇海底、起伏海面和反射地層深度等三個環境因素對鬼波參數的影響。采用兩層模型簡單模擬各個因素引起的延遲時間擾動，基礎模型如圖2a所示。海底與海面相距500 m，選取的檢波器位于炮檢距為1000 m、深度為20 m的位置，海水中波場傳播速度固定為1500 m/s，采用50 Hz的Ricker子波，各個模型參數因考察影響因素不同而改變。圖2b是海底情況，海底反射情況復雜，此處僅以鬼波的海底反射點位置抬升簡單演示。圖2c對應于海面起伏情況，檢波器的深度因此變化。圖2d為目的層深度對延遲時間的影響，為操作方便對該情況進行簡化，以海底深度變化查看目的層深度引起的延遲時間擾動。

圖3a～圖3c為圖2b～圖2d模型對應的數值模擬記錄，為同一個檢波器在影響因素改變時接收到的記錄，用于觀察延遲時間的變化，此處已將一次波校平，每個記錄的第一道對應影響因素未改變的情況；圖3d～圖3f為圖3a～圖3c減去各自第一道記錄的結果，這就使得各個因素對延遲時間的影響清晰可見。隨著鬼波在海底的反射點位置抬升，延遲時間隨之減小(圖3a)，當上抬高度到達1.0 m后差值記錄(圖3d)能量就不能忽視；由圖3b可見，當海面上升，鬼波與一次波的延遲時間隨之增加，當海面上升0.5 m后，鬼波就與改變前有較大的相位差(圖3h)；隨著目的層深度的增加，對應的延遲時間也隨之變化(圖3c、圖3f)，且前期變化快，后期漸緩，深度相差較大時，延遲時間差也對應增大。深度越大，反射時間越長，可見不同反射時間的延遲時間會有些許差距，但當目的層越深，鬼波與一次波的路徑就越接近，反射地層深度造成的延遲時間擾動就越小。

圖2 延遲時間影響因素

圖3 海底模型及數值模擬記錄

上述只是一個簡單的模擬，實際可能的情況要復雜得多。如海面是隨時間起伏的，則延遲時間也對應增大或減小，是時變的，在整個拖纜上海面的影響也是變化的，故而該因素對應的延遲時間擾動也是空變的；海底情形多變，并不僅是上抬，還有如傾斜、下陷、凸起等，影響復雜；反射地層深度造成的擾動變化快、慢難以統一度量，各因素交錯作用。另外，炮檢距也會對各個因素造成的延遲時間擾動產生影響。

以上情況說明，各因素作用會讓鬼波參數產生隨時間和空間變化的不容忽視的擾動，這使得鬼波參數在實際采集中難以估算。因此，在處理海洋地震記錄時，需要相對地震記錄求最優平均延遲時間和振幅差異系數，獲得一個整體最優的處理結果。

4 參數尋優原理

4.1 尋優原理

設Gg+(xj,pj,ω)=-Rie-iω Δ τie-iωpjxi,GG-(xj,pj,ω)=-e-iω Δ τie-iωpjxi/Ri，則

(16)

(17)

(18)

(19)

然后，對Q+進行反傅里葉變換得時間—空間域記錄q+，對Q-進行反傅里葉變換得時間—空間域記錄q-。將尋優判定函數設置為q+與q-乘積絕對值之和，即

(20)

式中S(xi)為第i道的判定函數；q+相當于原記錄加上延拓后記錄；q-相當于原記錄減去延拓后記錄。當鬼波與一次波的延遲時間和振幅差異系數最優時，延拓后記錄的一次波移動到原記錄中鬼波的位置，q+中的原記錄鬼波與延拓后記錄的一次波相加抵消，得到如式(17)右邊所表達的記錄結果，即一次波和延拓后的鬼波；當延遲時間或振幅差異系數有誤差時，原記錄的鬼波與延拓后記錄的一次波將不能完全抵消，會有殘余噪聲。q-與q+相對應，在q+中振幅因原記錄鬼波與延拓記錄一次波相加削弱的時間段，q-則因原記錄鬼波與延拓后記錄的一次波相減增強，鬼波參數最優時便如式(18)右邊所示。

當q+與q-相乘并取絕對值相加時，原記錄鬼波所在時間段的振幅變化是影響相乘結果的主要原因。若延遲時間和振幅差異系數最優，q+中原記錄鬼波去除徹底，q-中增強后的原記錄鬼波與q+中極小的原記錄鬼波殘留相乘并取絕對值相加后結果就較??；相對的若延遲時間和振幅差異系數有誤差，則q+中原記錄鬼波所在時間段鬼波殘留較多，與q-中增強的鬼波進行相乘并取絕對值相加時結果就會偏大。由此，通過兩個新構建記錄的波形特征可對不同延遲時間和不同振幅差異系數的記錄處理結果進行判定，找到最優處理結果對應的延遲時間和振幅差異系數。每道最小S對應的延遲時間和振幅差異系數即為所求結果。近炮檢距記錄可適當增加延拓次數以改善效果。

4.2 滑動時窗

從本文的正演地震記錄(圖8a)提取出一道記錄，即圖4a展示的第100道記錄。圖4a標注了該道6個主要由一次波和鬼波構成的反射波組，分別提取了這6個反射波組的鬼波延遲時間(圖4b)，可以看到鬼波的延遲時間是隨時間波動的，前期變化較快，后期平緩，可知鬼波參數是時變的。

圖4 第100道記錄點位圖(a)及對應的延遲時間示意圖(b)

圖5為滑動時窗示意圖，滑動時窗中每一個時窗長度對應一個最優延遲時間和振幅差異系數。分時窗時應該盡量將鬼波和上行波分在同一個時窗，前、后滑動窗口應該有重疊部分，以保證地震數據的全覆蓋。分時窗時窗口不宜太大，太大則前、后鬼波參數相差較大，導致無法求得一個合適的延遲時間和振幅差異系數；窗口也不宜太小，太小則計算量增加，時間方向滑動窗口約為200～500 ms即可，在該時段采集環境對鬼波參數的影響相對一致，空間方向長度根據處理情況調整。開始計算時，可先選擇一炮記錄進行時窗長度試算，在盡量使鬼波與一次波在同一個窗口的前提下，以一定的間隔增加時窗長度，將每一個時窗長度在炮記錄上進行處理測試，分別查看應用效果，選擇其中效果較好的時窗長度中的高值，將之擴散應用到整個接收記錄。在記錄長度的中、后段可以適當增加時窗長度，減小計算量。具體的窗口可根據實際處理情況實時修改。

圖5 滑動時窗示意圖

4.3 鬼波壓制步驟

(1)劃分滑動時間窗口，保證覆蓋整個地震記錄；

(3)在所給范圍內枚舉延遲時間，對于每一個延遲時間Δτ，振幅差異系數不變，按式(16)～式(20)計算對應的S值；

(4)對每個延遲時間的S比較大小，最小值所對應延遲時間即為最優延遲時間；

(5)在振幅差異系數范圍內枚舉鬼波振幅差異系數，保持步驟(4)所得到的延遲時間不變，對于每一個振幅差異系數，由式(16)～式(20)計算對應的S值；

(6)比較每個振幅差異系數的S值，最小值所對應的振幅差異系數即為所求最優鬼波振幅差異系數；

(7)用所得的振幅差異系數，重復步驟(3)和步驟(4)，得到新的延遲時間，再由新延遲時間得到新的鬼波振幅差異系數，據此不斷迭代，直至到達最大迭代次數；

(8)以所得到的延遲時間和振幅差異系數，由式(8)和式(9)計算去鬼波記錄U(p,ω)，再進行頻率域Radon反變換和傅里葉反變換，得到時間—空間域的去鬼波記錄U(x,t)；

(9)對每個窗口內的數據重復步驟(2)～步驟(8)，最后將所有窗口數據整合成一個完整的記錄。

利用延拓尋找最優平均鬼波延遲時間和振幅差異系數進行鬼波壓制方法的流程圖如圖6所示。

圖6 鬼波壓制流程

5 合成數據試算

為驗證本文方法去除鬼波的效果，利用圖7a所示復雜鹽丘模型的正演模擬記錄進行試算。斜纜采用圖7b拋物線式斜纜，炮點位于800 m處，檢波器沉放深度為6～50 m，記錄共250道，道間距為10 m，時間采樣間隔為0.002 s，記錄長度為2.9 s，采用主頻為35 Hz的Ricker子波。

圖7 復雜鹽丘模型(a)與斜纜深度隨炮檢距變化曲線(b)

在原始炮集(圖8a)中，因為海面的強反射，一次波后緊跟著一個反相的鬼波，鬼波與一次波的延遲時間隨著炮檢距的增大呈先增后減的狀態。圖8b是拖纜鬼波壓制后的炮集。由圖可見，經本文方法處理后，一次波后的反相同相軸基本消失，整個記錄的相位變得相對單一。對參數尋優時需將數據轉換到頻率—慢度域進行延拓，再換算到時間—空間域進行處理結果判定。將原始炮集轉換到時間—慢度域時(圖8c)，可發現與圖8a相似的特征，鬼波與一次波對應的同相軸先逐漸分開再靠攏，兩者相位表現相反，各同相軸與時間—空間域的同相軸相對應。圖8d是去鬼波炮集的時間—慢度域記錄，其中與鬼波對應的同相軸得到有效壓制，鬼波去除相對干凈。

圖8 不同域去鬼波前、后的合成記錄

從圖9a的近炮檢距記錄可以看到鬼波與一次波隨著炮檢距的增大逐漸分開的過程，經鬼波壓制后(圖9b)，鬼波同相軸的能量得到極大程度的衰減?？梢姳疚姆椒▽聿ㄓ兄己玫娜コЧ?。鬼波會與一次波相互干涉，使得頻譜出現陷波點，陷波點處頻譜能量會有較大衰弱(圖10a)。在原始炮集的頻率—空間域振幅譜(圖10a)中可以看到陷波頻率隨檢波器深度變化，顯現出了斜纜陷波特征的多樣性。在去除鬼波后，陷波點處的能量得到有效地恢復，頻譜變得連續，拓寬了記錄的有效帶寬(圖10b)，驗證了本文方法可壓制鬼波、彌補由鬼波引起的陷波點能量損失。

圖9 時間—空間域去鬼波前(a)、后(b)合成記錄局部圖

圖10 去鬼波前(a)、后(b)合成記錄振幅譜

6 實際數據試算

為進一步驗證本文方法的有效性，對實際海上斜纜采集數據進行試算，該數據震源深度為5 m，道間距為12.5 m，時間采樣間隔為0.002 s，總計為300道，斜纜大致呈拋物線狀，拖纜沉放深度為5～50 m。圖11a是經過直達波切除、去噪等預處理的斜纜炮集記錄。在實際資料中，因各種采集環境因素造成的擾動，鬼波與一次波的實際延遲時間與理論上有一定誤差，鬼波與一次波的振幅差異系數也不再與假定的相同。

海上斜纜數據的反射波主要出現在2 s后(圖11a)，鬼波緊隨在上行波后，兩者間的距離隨著炮檢距的增大呈先逐漸分開后漸漸靠近的變化，相位表現相反。經過鬼波參數的延拓尋優，擾動的鬼波參數得到修正，經過鬼波處理后，緊隨上行波的反相同相軸得到極大壓制(圖11b)，記錄更加清晰，有效信號更加突出，說明本文方法可較好地消除斜纜的鬼波。實際資料變換到時間—慢度域，同樣出現鬼波與一次波成對出現的特征(圖11c)，兩者相位相反，鬼波與一次波同相軸隨著慢度增加先逐漸分離再慢慢并攏，應用本文方法進行鬼波處理后，時間—慢度域中鬼波對應的同相軸基本消失(圖11d)，表明鬼波得到有效去除。從圖12的局部對比圖也能看到一次波后的同相軸能量得到較大衰減，炮記錄相位更加單一，鬼波得到較好地壓制，進一步說明了基于延拓尋找最優參數的鬼波壓制方法可較好地去除由海面虛反射引起的鬼波，提高資料的分辨率。

圖12 圖11去鬼波前(a)、后(b)實際資料局部放大顯示

從圖13a可以看到實際資料中鬼波引起的陷波特征以及斜纜采集數據中陷波頻率多樣性的特點。經過鬼波處理后(圖13b)，陷波處的能量得到有效補償，振幅譜連續性增加，使得地震記錄的頻帶得到拓寬，較好地去除了鬼波帶來的陷頻影響，說明本文方法在實際數據處理中的有效性。在鬼波去除前后第100道至第150道記錄振幅譜對比圖中(圖14)，可見到多個陷波頻率以及陷波頻率處明顯的能量低值，經本文方法去除鬼波后，振幅譜變得更加連續，低頻和高頻的能量都有所增益，斜纜陷波頻率多樣性的特點得到較好利用，提高了資料的分辨率，表明本文方法對資料的寬頻處理有著良好的作用。

圖13 去鬼波前(a)、后(b)實際資料振幅譜

圖14 鬼波壓制前、后第100～150道振幅譜對比

7 結論

通過對斜纜合成數據以及海上實際斜纜資料的處理分析，可以得到如下結論：

(1)崎嶇海底、起伏海面、地層反射構造及拖纜沉放深度誤差等因素直接影響到鬼波相對于一次波的振幅及延遲時間的擾動，繼而影響到傳統方法鬼波壓制效果；

(2)結合滑動時窗方法，基于頻率—慢度域延拓算子，通過交替迭代求解，得到最優鬼波相對于一次波的振幅差異系數和延遲時間；

(3)基于頻率—慢度域延拓的鬼波參數最優化拖纜鬼波壓制方法應用于合成和實際海上拖纜地震炮集記錄的鬼波壓制，鬼波參數擾動的影響得到消除，鬼波壓制效果得以提升，地震資料的頻帶明顯拓寬。