封 強(qiáng) 韓立國(guó) 楊 帆

(吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院，吉林長(zhǎng)春130026)

0 引言

海洋地震勘探中，受海水面虛反射影響，采集的原始地震數(shù)據(jù)包括地下反射波和鬼波。由于鬼波的陷頻特性，限制了地震勘探的頻帶寬度，降低了地震資料的分辨率[1]。因此，要獲得寬頻帶地震資料，必須消除鬼波的影響。

在此背景下，業(yè)界陸續(xù)推出并應(yīng)用了上下纜[2-4]、上下源[3]、雙檢[4]等海洋地震勘探技術(shù)。近年來(lái)，變深度纜采集技術(shù)及相關(guān)的去鬼波算法受到學(xué)者關(guān)注。即摒棄常規(guī)拖纜采集的將拖纜固定在某一深度的方式，而采用CGG公司首創(chuàng)的檢波器沉放深度隨炮檢距增大而增大的變深度纜采集方式[5]。變深度纜中不同檢波器沉放深度具有陷頻多樣性，疊加不同道集，能同時(shí)在低頻端和高頻端兩個(gè)方向拓寬頻帶[6]。針對(duì)鬼波壓制，Soubaras[7]提出聯(lián)合反褶積方法； 張振波等[8]將斜纜寬頻地震勘探技術(shù)應(yīng)用于珠江口盆地； 許自強(qiáng)等[9]采用最優(yōu)化聯(lián)合反褶積算法壓制鬼波，在實(shí)際數(shù)據(jù)處理中取得良好效果； Wang等[10-11]基于Bootstrap方法做加窗處理和鏡像偏移，估計(jì)鬼波延遲時(shí)間，在f-x-y域和τ-p域壓制鬼波； 王沖等[12]采用最小二乘反演迭代算法，獲得準(zhǔn)確的鬼波延遲時(shí)間； 葉林等[13]提出利用峰值因子搜尋最優(yōu)延遲時(shí)間壓制鬼波。為了提高計(jì)算的穩(wěn)定性，張威等[14]基于LSMR算法反演求解海面觀測(cè)的上行波場(chǎng)，再延拓到拖纜以壓制斜纜鬼波； 吳娟等[15]在高斯束偏移中加入穩(wěn)定求解的鬼波補(bǔ)償算子，更利于消除鬼波影響； 王建花等[16]基于波場(chǎng)延拓和反演，增加約束條件壓制變深度纜鬼波； 馬繼濤等[17]基于波場(chǎng)外推和閾值截?cái)鄩褐乒聿ā?/p>

大多數(shù)鬼波壓制方法假設(shè)海面反射系數(shù)為-1，但實(shí)際上海面反射系數(shù)并不等于-1，同時(shí)，不同時(shí)間海面各處反射系數(shù)并不完全相同。Masoomzadeh等[18]在估計(jì)鬼波延遲時(shí)間時(shí)，假設(shè)海面反射系數(shù)是頻率和慢度的函數(shù)。王沖等[19]基于理論下行波與實(shí)際下行波之間的平方誤差最小理論，在頻率慢度域計(jì)算了鬼波與一次反射波振幅系數(shù)和延遲時(shí)間。

在重點(diǎn)調(diào)研了變深度纜鏡像記錄及聯(lián)合反褶積算法相關(guān)文獻(xiàn)的基礎(chǔ)上，本文提出一種基于非高斯性最大化的變深度纜鬼波壓制方法，以負(fù)熵為非高斯性度量，采用滑動(dòng)時(shí)空數(shù)據(jù)窗口克服海面反射系數(shù)隨時(shí)間變化和延遲時(shí)間隨地層深度變化問(wèn)題，獲取最優(yōu)海面反射系數(shù)和延遲時(shí)間，用聯(lián)合反褶積算法較徹底地壓制了鬼波。

1 變深度纜聯(lián)合反褶積算法

變深度纜某一深度zi所接收到的實(shí)際波場(chǎng)p(x,zi,t)是由深部地層反射的上行波場(chǎng)u(x,zi,t)與該上行波場(chǎng)因海面反射產(chǎn)生的下行波場(chǎng)d(x,zi,t)之和構(gòu)成，即

p(x,zi,t)=u(x,zi,t)+d(x,zi,t)

(1)

對(duì)式(1)沿時(shí)間t做一維傅里葉變換，得到

P(x,zi,ω)=U(x,zi,ω)+D(x,zi,ω)

(2)

在圖1中：沉放深度z1=zr處拖纜接收到的總波場(chǎng)為p(x,z1,t)，上、下行波場(chǎng)分別為u(x,z1,t)和d(x,z1,t)； 拖纜位置往下zr處設(shè)定為“假想海面”，以此為對(duì)稱面，往下zr處即是相對(duì)于變深度纜接收點(diǎn)的鏡像接收點(diǎn)，其沉放深度z2=3zr，則該假想拖纜接收到的總波場(chǎng)為p(x,z2,t)，對(duì)應(yīng)的上、下行波場(chǎng)分別為u(x,z2,t)和d(x,z2,t)。各波場(chǎng)對(duì)時(shí)間t的一維傅里葉變換分別記為P(x,z1,ω)、P(x,z2,ω)、U(x,z1,ω)、U(x,z2,ω)、D(x,z1,ω)和D(x,z2,ω)。

圖1 波場(chǎng)傳播示意圖

要求鏡像拖纜接收的下行波等于實(shí)際拖纜接收的上行波，即

D(x,z2,ω)=U(x,z1,ω)

(3)

假設(shè)海面反射系數(shù)為-r，下行波場(chǎng)是上行波場(chǎng)在海面的反射而產(chǎn)生，因此下行波場(chǎng)可以由單程波場(chǎng)延拓算子得到

(4)

將式(3)代入式(4)得

(5)

再將式(3)～式(5)代入式(2)，得到

P(x,z1,ω)=U(x,z1,ω)[1-re-jk(z2-z1)]

(6)

(7)

聯(lián)立式(6)與式(7)，可得

P(x,z2,ω)=-rejk(z2-z1)P(x,z1,ω)

(8)

通過(guò)式(8)可得到鏡像記錄，鏡像記錄等于原始記錄施加一個(gè)算子，鏡像算子可以表示為

J(x,ω)=-rejk(z2-z1)=-rejω τ

(9)

式中τ為上行波與下行波之間的延遲時(shí)間。

式(6)、式(7)可改寫成

P(x,z1,ω)=U(x,z1,ω)G1(x,ω)

(10)

P(x,z2,ω)=U(x,z2,ω)G2(x,ω)

(11)

式中

G1(x,ω)=1-re-jk(z2-z1)=1-re-jω τ

(12)

(13)

鏡像記錄也可表示為

P(x,z2,ω)=-re-jω τP(x,z1,ω)

(14)

上述問(wèn)題的最小平方解是

U(x,z1,ω)

(15)

上標(biāo)“*”表示共軛運(yùn)算。在實(shí)際地震記錄與其鏡像記錄已知的情況下，式(15)可求得變深度纜的上行波場(chǎng)U(x,z1,ω)，即鬼波壓制后的頻率域波場(chǎng)。

2 非高斯性最大化的參數(shù)搜索

實(shí)際地震信號(hào)具有非高斯分布的特性[20]。由中心極限定理可知：多個(gè)相互獨(dú)立的隨機(jī)信號(hào)的和比單個(gè)信號(hào)更趨于高斯分布。實(shí)際地震記錄包括一次波和鬼波，即一次波的非高斯性大于實(shí)際地震記錄的非高斯性。因此，鬼波壓制效果最好的地震記錄的非高斯性最大。

在所有方差相等的隨機(jī)變量中，高斯變量具有最大熵，非高斯性越強(qiáng)熵越小。負(fù)熵J可作為一種非高斯性度量。為便于計(jì)算，選用負(fù)熵近似式

J(y)=[E(G(y))-E(G(v))]2

(16)

式中：v是具有零均值、單位方差的高斯信號(hào)；y是具有零均值、單位方差的隨機(jī)信號(hào)；E(·)表示求期望值；G(·)是可選非線性函數(shù)，這里選擇G(x)=-e0.5x2。負(fù)熵近似式對(duì)噪聲不敏感，具有良好的魯棒性，計(jì)算快速。

在上述聯(lián)合反褶積算法中，海面反射系數(shù)和鬼波延遲時(shí)間是影響上行波準(zhǔn)確性的兩個(gè)最重要參數(shù)，但實(shí)際上并不能得到準(zhǔn)確的海面反射系數(shù)和鬼波延遲時(shí)間。建立基于非高斯性最大化的目標(biāo)函數(shù)

(17)

其中

式中：u0(n)是對(duì)u(n)零均值化和方差歸一化的結(jié)果；v(n)是均值為0、方差為1的高斯信號(hào)。當(dāng)u0(n)或u(n)為高斯信號(hào)時(shí)，目標(biāo)函數(shù)等于0； 而當(dāng)u0(n)或u(n)越趨向于非高斯信號(hào)，目標(biāo)函數(shù)越大于0； 反之亦然。因此，當(dāng)目標(biāo)函數(shù)最大時(shí)，代表此時(shí)的參數(shù)為最優(yōu)的反射系數(shù)和延遲時(shí)間。

3 滑動(dòng)數(shù)據(jù)窗口處理

同一檢波點(diǎn)接收到的鬼波與反射波之間的延遲時(shí)間隨地層深度變化(圖2)，同一道的鬼波延遲時(shí)間不完全相同。海面反射系數(shù)并不是恒定值，而是隨時(shí)間發(fā)生變化。單道地震記錄只選定固定的延遲時(shí)間和反射系數(shù)進(jìn)行去鬼波處理，很難取得理想的壓制效果。因此，必須沿時(shí)間和空間進(jìn)行窗口劃分，離散化地震記錄，使同一窗口內(nèi)的地震記錄的延遲時(shí)間近似相等，同時(shí)確定海面反射系數(shù)。

圖2 鬼波與一次反射波間延遲時(shí)隨炮檢距和旅行時(shí)變化示意圖

變深度纜各接收點(diǎn)的鬼波延遲時(shí)間與拖纜沉放深度有關(guān)，隨著沉放深度增大，鬼波延遲時(shí)間也逐漸增大。因此，窗口在空間方向上不能太大，否則很難找到一個(gè)較準(zhǔn)確的延遲時(shí)間。實(shí)際海面是粗糙起伏的，反射系數(shù)也不是一成不變的； 同樣，在時(shí)間域也需選取適當(dāng)窗口。在本文中，選擇滑動(dòng)窗口為3～5道、300～400ms(圖3)。在此基礎(chǔ)上，需擴(kuò)大搜索鬼波延遲時(shí)間范圍的上、下限，選擇適當(dāng)步長(zhǎng)，以便能覆蓋窗口內(nèi)的各道的準(zhǔn)確的延遲時(shí)間?；瑒?dòng)數(shù)據(jù)窗口在地震記錄上有規(guī)律地移動(dòng)，后續(xù)窗口與前一窗口有重疊部分，保證地震記錄被各窗口完全覆蓋。對(duì)地震記錄做滑動(dòng)窗口處理，使窗口內(nèi)去鬼波記錄非高斯性最大化，從而獲得窗口內(nèi)的最優(yōu)延遲時(shí)間和反射系數(shù)，克服鬼波參數(shù)變化問(wèn)題?；瑒?dòng)窗口完全覆蓋整條記錄后，將各窗口鬼波延遲時(shí)間和海面反射系數(shù)代入式(12)和式(13)，繼而通過(guò)式(15)得到去鬼波后的頻率—空間域地震記錄。

圖3 地震記錄滑動(dòng)窗口示意圖

4 數(shù)據(jù)測(cè)試

4.1 驗(yàn)證非高斯性最大化方法的有效性

為了驗(yàn)證基于非高斯性最大化方法的有效性，首先模擬一次波和鬼波，選擇鬼波參數(shù)：反射系數(shù)r為-0.960，鬼波延遲時(shí)間τ為45.00ms，設(shè)置反射系數(shù)范圍[-1.0∶0.001∶-0.6]、 延遲時(shí)間(ms)范圍[40∶0.01∶50]，利用聯(lián)合反褶積計(jì)算上行波場(chǎng)； 再計(jì)算得到二維負(fù)熵圖譜(圖4)，利用上述方法搜索負(fù)熵最大值，搜尋最優(yōu)參數(shù)。估算參數(shù)re=-0.961和τe=45.00ms很接近初始設(shè)定參數(shù)。

圖5顯示部分參數(shù)對(duì)的去鬼波結(jié)果，第一道為原始鬼波信號(hào)，第二道為獲得的最優(yōu)參數(shù)去鬼波結(jié)果，其余道是參數(shù)空間中隨機(jī)選取的部分參數(shù)去鬼波結(jié)果。通過(guò)對(duì)比，可見不同的鬼波參數(shù)對(duì)鬼波壓制效果產(chǎn)生了很大影響，利用負(fù)熵估計(jì)出來(lái)的參數(shù)能很好地壓制鬼波。因此，基于非高斯性最大化能夠選擇最優(yōu)的鬼波參數(shù)。

4.2 簡(jiǎn)單模型測(cè)試

為了檢測(cè)基于非高斯性最大化的變深度纜聯(lián)合反褶積算法壓制鬼波的效果，對(duì)水平層狀模型的合成變深度纜單炮記錄做鬼波壓制試處理。模擬時(shí)采用主頻為50Hz的雷克子波，記錄道數(shù)為101，道間距為6m，采樣間隔為0.5ms，記錄時(shí)間為0.9s，炮檢距為0，單邊放炮，拖纜沉放深度為6～50m。

圖6a為正演模擬的變深度纜單炮記錄，對(duì)其加入隨機(jī)噪聲以測(cè)試本文方法的抗噪性。隨著炮檢距增大，拖纜沉放深度逐漸增大，鬼波與一次波之間的延遲時(shí)間相應(yīng)增大，在地震記錄上表現(xiàn)為同相軸由開始時(shí)的相互重疊到慢慢分離，且兩者距離越來(lái)越大。在壓制鬼波過(guò)程中，首先設(shè)定反射系數(shù)搜索下限為-1.0，搜索上限為-0.6，步長(zhǎng)為0.01；延遲時(shí)間τ(=2d/c)的搜索上限為1.5τ，下限為0.5τ，步長(zhǎng)為0.1ms。滑動(dòng)窗口在時(shí)空域地震記錄上有規(guī)律地移動(dòng)，獲得各窗口的二維負(fù)熵圖譜，搜尋其中的最大值，分別得到各窗口內(nèi)最優(yōu)海面反射系數(shù)和鬼波延遲時(shí)間； 基于式(15)，對(duì)模擬的原始地震記錄和鏡像記錄(圖6b)應(yīng)用聯(lián)合反褶積法求得去鬼波波場(chǎng)(圖6c)。

鏡像記錄與原始地震記錄存在明顯的同相軸相位差異： 鏡像記錄首先接收到的是鬼波，因此其同相軸相位最初是負(fù)的；而原始地震記錄首先接收到的是一次反射波，其同相軸相位最初是正的。

對(duì)比圖6a與圖6c，可見運(yùn)用本文方法對(duì)原始地震記錄進(jìn)行處理后，在含噪情況下鬼波已被較徹底地去除，只剩下一次反射波； 同時(shí)表明噪聲的存在對(duì)本文方法去鬼波效果的影響極小。

總之，模型測(cè)試結(jié)果表明：本文方法在含噪情況下能有效壓制鬼波，且抗噪性強(qiáng)、精度高。

4.3 實(shí)際數(shù)據(jù)處理

用實(shí)際采集的變深度纜數(shù)據(jù)進(jìn)行非高斯性最優(yōu)化去鬼波算法測(cè)試。施工參數(shù)為：拖纜沉放深度5～50m，道間距12.5m，采樣間隔2ms，記錄長(zhǎng)度7.2s。先對(duì)該變深度纜單炮記錄做預(yù)處理(圖7a，選取2.0～4.0s，1～110道)并得到對(duì)應(yīng)鏡像記錄(圖7b)；再應(yīng)用本文方法做濾波處理，得到去鬼波后地震記錄(圖7c)。

由于鬼波的存在，在反射波同相軸之后會(huì)出現(xiàn)與其相位相反的虛假同相軸(圖7a)，降低了地震剖面的分辨率。因鬼波延遲時(shí)間和反射系數(shù)都能較準(zhǔn)確地獲取，故在利用本文方法獲得的鏡像記錄(圖7b)中對(duì)應(yīng)位置處都會(huì)呈現(xiàn)明顯的相位相反。

圖6 模擬的原始(a)、鏡像(b)及去鬼波后(c)的地震記錄

從實(shí)際地震記錄及其鬼波壓制后記錄的局部放大(圖8)可見：經(jīng)去鬼波處理后地震記錄(圖8b)中鬼波被徹底壓制，只剩下一次波同相軸，且與原始記錄中的反射波基本重合。充分說(shuō)明本文提出的基于非高斯性最大化的最優(yōu)化聯(lián)合反褶積算法對(duì)實(shí)際海上變深度纜數(shù)據(jù)中的鬼波具有較強(qiáng)壓制作用。

圖7 實(shí)際(a)、鏡像(b)及壓制鬼波后(c)的地震記錄

圖8 實(shí)際地震記錄鬼波壓制前(a)、后(b)的局部放大

圖9為實(shí)際地震數(shù)據(jù)鬼波壓制前、后的頻譜。對(duì)比發(fā)現(xiàn)用本文方法壓制鬼波后(圖9b)頻譜中的陷波點(diǎn)已基本消除，陷波頻率的能量得到了補(bǔ)償，低頻和高頻能量都有顯著提高。

圖9 實(shí)際地震數(shù)據(jù)鬼波壓制前(a)、后(b)的頻譜

5 結(jié)束語(yǔ)

文中分析、對(duì)比了模擬的和實(shí)際的變深度纜單炮記錄鬼波壓制方法及效果，并得到如下結(jié)論：

(1)以負(fù)熵度量去鬼波結(jié)果的非高斯性，可得到精確的海面反射系數(shù)和鬼波延遲時(shí)間，從而獲得理想的鬼波壓制效果。

(2)采用滑動(dòng)時(shí)空數(shù)據(jù)窗口可克服海面反射系數(shù)和鬼波延遲時(shí)間的變化問(wèn)題，獲得精確鏡像記錄和去鬼波記錄，提高了變深度纜資料處理精度。

(3)在同相軸干涉嚴(yán)重的復(fù)雜構(gòu)造區(qū)及大炮檢距情況下，時(shí)空域壓制鬼波存在局限性，鬼波壓制效果欠佳。這是后續(xù)研究的課題。