頻率—慢度域延拓關(guān)鍵參數(shù)優(yōu)化的鬼波壓制方法

青杰,顧漢明*,王建花

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

0 引言

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

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

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

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

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

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

(1)

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

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

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

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

(2)

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

(3)

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

(4)

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

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)

式中：μ為阻尼系數(shù)；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 鬼波參數(shù)擾動分析

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

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

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

圖2 延遲時間影響因素

圖3 海底模型及數(shù)值模擬記錄

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

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

4 參數(shù)尋優(yōu)原理

4.1 尋優(yōu)原理

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

(16)

(17)

(18)

(19)

然后，對Q+進(jìn)行反傅里葉變換得時間—空間域記錄q+，對Q-進(jìn)行反傅里葉變換得時間—空間域記錄q-。將尋優(yōu)判定函數(shù)設(shè)置為q+與q-乘積絕對值之和，即

(20)

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

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

4.2 滑動時窗

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

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

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

圖5 滑動時窗示意圖

4.3 鬼波壓制步驟

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

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

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

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

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

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

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

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

利用延拓尋找最優(yōu)平均鬼波延遲時間和振幅差異系數(shù)進(jìn)行鬼波壓制方法的流程圖如圖6所示。

圖6 鬼波壓制流程

5 合成數(shù)據(jù)試算

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

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

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

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

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

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

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

6 實際數(shù)據(jù)試算

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

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

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

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

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

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

7 結(jié)論

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

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

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

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