基于粗糙海面的最小二乘殘差法變深度纜接收點(diǎn)鬼波壓制技術(shù)*

李志鵬 何兵壽 楊佳佳 李凱瑞

(1. 中國海洋大學(xué) 山東青島 266100； 2. 青島海洋科學(xué)與技術(shù)國家實(shí)驗(yàn)室海洋礦產(chǎn)資源評價(jià)與探測技術(shù)功能實(shí)驗(yàn)室 山東青島 266071；3. 海底科學(xué)與探測技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室 山東青島 266100； 4. 中國地質(zhì)調(diào)查局青島海洋地質(zhì)研究所 山東青島 266100)

海上拖纜采集的地震數(shù)據(jù)中，不僅包含來自海底反射界面的上行波場，也包含上行波場傳播至海面再由海面反射向下傳播的下行鬼波波場。鬼波的存在會(huì)導(dǎo)致明顯的陷波效應(yīng)，使地震數(shù)據(jù)的頻帶變窄，分辨率降低[1]；同時(shí)，連續(xù)發(fā)育的鬼波也會(huì)在地震剖面上產(chǎn)生虛假同相軸，給地震地質(zhì)解釋造成困擾[2]。

目前，業(yè)界進(jìn)行鬼波壓制的思路主要有2種：一是在野外采集中采用合適的震源組合、電纜形態(tài)或電纜組合減小鬼波的影響[3-9]；二是在資料處理過程中采用合適的算法從實(shí)測數(shù)據(jù)中消除鬼波[10-15]。第一種思路的常用方法包括：①通過立體槍陣的延時(shí)激發(fā)[3-4]壓制震源鬼波；②采用雙檢采集[5]、上下纜采集[6]和變深度纜采集[7-9]等方式壓制接收點(diǎn)鬼波。上述方法均在實(shí)際生產(chǎn)中取得了效果[16]，但實(shí)際采集中往往存在電纜垂向定位難度大、采集成本高[17]等缺陷，而且受采集設(shè)備的嚴(yán)格限制[18]，僅憑野外采集手段無法實(shí)現(xiàn)鬼波的完全壓制。第二種思路的常用方法包括：①對于雙檢采集數(shù)據(jù)，首先進(jìn)行水、陸檢數(shù)據(jù)的一致性處理，然后對兩種檢波器接收到的信息進(jìn)行特征項(xiàng)匹配，進(jìn)而進(jìn)行信息合并[10]，以達(dá)到消除鬼波的目的。②對上下纜數(shù)據(jù)，采用相位校正算法[11]或移位和減去算法[12]進(jìn)行相位修正，然后加權(quán)合并去除鬼波。③對于變深度纜數(shù)據(jù)，常采用聯(lián)合反褶積算法[13-14]或最小二乘反演算法[15]進(jìn)行鬼波壓制。

上述鬼波壓制方法一般以海面水平為假設(shè)前提，此時(shí)海面的反射系數(shù)為-1；而事實(shí)上，受風(fēng)速、風(fēng)向、涌浪、風(fēng)浪和人為活動(dòng)(如航船等)[19]等因素的影響，海面往往是起伏的，此時(shí)海面的反射系數(shù)不一定是-1[20]。這種理論假設(shè)與實(shí)際情況的不相符會(huì)導(dǎo)致上述技術(shù)難以取得理想效果，因此必須研究基于起伏海面的鬼波壓制技術(shù)。本文提出了一種基于粗糙海面的最小二乘殘差法變深度纜接收點(diǎn)鬼波壓制方法，并開展了不同模型的合成數(shù)據(jù)測試，從而實(shí)現(xiàn)了在粗糙海面條件下對接收點(diǎn)鬼波波場的有效壓制。

1 粗糙海面建模與反射系數(shù)計(jì)算

真實(shí)海面是粗糙且不可預(yù)測的，對于海洋數(shù)據(jù)的鬼波壓制，這種海面條件會(huì)導(dǎo)致后期數(shù)據(jù)處理中反射系數(shù)與鬼波延遲時(shí)差等參數(shù)的估算誤差，從而會(huì)造成鬼波記錄時(shí)差、能量及相位的計(jì)算誤差，進(jìn)而降低鬼波壓制效果。理論上，基于起伏海面的鬼波壓制技術(shù)可以消除上述誤差，而實(shí)現(xiàn)起伏海面鬼波壓制的前提是海面起伏量在激發(fā)、接收時(shí)刻的求取與海面成像，而事實(shí)上，受成本、工期和設(shè)備等因素的限制，野外采集時(shí)往往不對海面起伏量進(jìn)行實(shí)時(shí)的密集測量，因此海面起伏的建模需要在室內(nèi)進(jìn)行。

一般來說，某個(gè)時(shí)刻的海面是在一定的高差范圍內(nèi)隨機(jī)起伏的，且起伏量滿足高斯分布[21]，因此可利用統(tǒng)計(jì)學(xué)原理對粗糙海面進(jìn)行建模。

1.1 描述粗糙海面的統(tǒng)計(jì)參量

依據(jù)蒙特卡洛建模方法[22-24]，描述粗糙海面的統(tǒng)計(jì)參量主要包括均方根高度σ、相關(guān)度函數(shù)C、相關(guān)長度l、功率譜密度W、均方根斜率δ和曲率半徑Rc，各參數(shù)的定義如下。

均方根高度σ是描述粗糙海面起伏程度的重要參數(shù)，其值越高，海面起伏程度越大。在拖纜地震勘探中可通過浮標(biāo)數(shù)據(jù)或其它觀測資料近似求取。在野外采集過程的時(shí)間跨度內(nèi)，該值可視為一常數(shù)，若浮標(biāo)等資料缺失，也可以通過最優(yōu)值搜索[25]預(yù)估其值。其定義為

(1)

相關(guān)度函數(shù)C定義為

C(r)=〈h(x),h(x+r)〉

(2)

式(2)中：r是海面上任意兩點(diǎn)之間的距離；〈〉是函數(shù)相關(guān)運(yùn)算符。相關(guān)度函數(shù)是描述海面上任意兩點(diǎn)高程相關(guān)程度的物理量，當(dāng)r=0時(shí)，C(r)最大且有C(r)=σ2。

利用海面均方根高度將相關(guān)度函數(shù)歸一化,可得到相關(guān)系數(shù)ρ

(3)

ρ(r)=1/e時(shí)所對應(yīng)的r值即為相關(guān)長度，記為l。相關(guān)長度是判斷粗糙海面上任意兩點(diǎn)是否相互獨(dú)立的標(biāo)準(zhǔn)之一。

粗糙海面起伏高度的功率譜密度W(k)是相關(guān)度的傅里葉變換結(jié)果[26],表示信號(hào)功率在不同頻率的分布情況，是單位頻率的平均功率量綱,即

(4)

功率譜密度與均方根高度存在以下關(guān)系[26]：

(5)

均方根斜率δ為粗糙海面上各點(diǎn)斜率的均方根值[20],其定義為

(6)

式(6)中：dx為x方向的采樣大小;dh為相距dx兩點(diǎn)對應(yīng)的高度差值。

曲率半徑R反映粗糙面上各點(diǎn)的彎曲程度[27],其定義為

(7)

1.2 粗糙海面成像與反射系數(shù)計(jì)算

粗糙海面可以看作是由大量不同振幅、頻率、相位的諧波疊加而成的，各諧波的振幅可以視作獨(dú)立的高斯隨機(jī)變量[28]，由此可得海面任一點(diǎn)x處起伏量的表達(dá)式為

(8)

由于構(gòu)成粗糙海面的各諧波的振幅是獨(dú)立的高斯隨機(jī)變量，且諧波函數(shù)的相位也滿足正態(tài)分布[29]，因此由式(2)和式(4)可求得高斯密度的相關(guān)函數(shù)和功率譜密度函數(shù)。定義離散波數(shù)kn=2πn/L，可得到長度為L的一維粗糙表面的諧波譜函數(shù)為

(9)

式(9)中：N(0,1)代表均值為0、方差為1的標(biāo)準(zhǔn)正態(tài)分布的隨機(jī)數(shù)；Δk為譜域相鄰的諧波樣本的空間波數(shù)差。再對F(kn)做傅里葉變換,可得到粗糙海面離散化后的函數(shù)關(guān)系為

(10)

式(10)中：xi=iΔx(i=-N/2+1,-N/2+2,…,N/2)，表示粗糙海面上第i個(gè)采樣點(diǎn)。

圖1是按上述理論的粗糙海面建模結(jié)果，其中圖1a為二維建模成像結(jié)果，圖1b是均方根值分別為0.1與0.2時(shí)抽取的y=20 m時(shí)的海面起伏情況。

圖1 粗糙海面建模結(jié)果

在短時(shí)間尺度內(nèi)，采用低階小斜率近似方法[29-30]建立時(shí)變海面復(fù)反射系數(shù)的計(jì)算模型，并由此求出對應(yīng)的粗糙海面反射系數(shù)[31]為

(11)

式(1)中：R0為理想光滑海面的反射系數(shù)，即-1；σ為粗糙海面的均方根振幅；s為一個(gè)和頻率有關(guān)的常數(shù)，其值由入射到海面的地震波頻率成分決定；α為地震波傳播到海面時(shí)的入射角；v為海水聲波速度。

2 基于粗糙海面的變深纜接收點(diǎn)鬼波壓制方法

基于粗糙海面的變深纜接收點(diǎn)鬼波的壓制主要分三步進(jìn)行：①求取粗糙海面的反射系數(shù)；②在τ-p域利用求取出的粗糙海面反射系數(shù)進(jìn)行上下行波波場分離，得到不包含鬼波能量的上行波波場；③將得到的海面上行波場延拓至預(yù)先定義的新的水平基準(zhǔn)面，最終得到壓制鬼波后的記錄。其中，粗糙海面反射系數(shù)的計(jì)算方法已由式(11)給出，在此不再贅述。

2.1 基于最小二乘殘差法的τ -p域波場分離

不考慮船舶反射等背景噪聲，海洋拖纜地震數(shù)據(jù)所記錄的波場p(xi,t)可以視為來自海底地層的上行反射波場u(xi,t)和該上行波傳播到海面后再由海面反射向下傳播的下行波場d(xi,t)(鬼波)之和[32]，即

p(xi,t)=u(xi,t)+d(xi,t)

(12)

如圖2所示，D點(diǎn)在某時(shí)刻接收到的總波場p可視為深部反射波傳播到D點(diǎn)的上行波場u1與由深部傳播至A點(diǎn)再由A點(diǎn)反射至D點(diǎn)的下行波場d1之和。B點(diǎn)為D點(diǎn)在海面的投影，對應(yīng)某炮點(diǎn)接收排列的第i道(i=1,2,3,…,n)。假設(shè)該時(shí)刻在B點(diǎn)同樣存在一個(gè)上行波u2，其傳播方向與u1相同，E點(diǎn)為B點(diǎn)在u1傳播路徑的投影，那么CE段即為u1比u2多傳播的距離。

圖2 接收點(diǎn)波場幾何關(guān)系示意圖

設(shè)u1到達(dá)C點(diǎn)與u2到達(dá)B點(diǎn)的時(shí)差為Δti,m，則有

(13)

式(13)中：αm為第m個(gè)出射角;pm為第m個(gè)射線參數(shù);zi為第i道相對于粗糙海面的沉放深度。

再假設(shè)一維海面在A—C這一范圍內(nèi)無起伏，則u1從D點(diǎn)傳播到C點(diǎn)的時(shí)間tui,m與d1從A點(diǎn)傳播到D點(diǎn)的時(shí)間tdi,m相同，將此時(shí)間記為tsi,m，有

(14)

設(shè)B點(diǎn)的橫坐標(biāo)為xi，則A點(diǎn)與C點(diǎn)的橫坐標(biāo)分別為xi+Δxi,m與xi-Δxi,m。存在如下關(guān)系：t時(shí)刻傳播到D點(diǎn)的上行波u1與t+tsi,m時(shí)刻傳播到C點(diǎn)的上行波場一致，t時(shí)刻傳播到D點(diǎn)的下行波場d1與t-tsi,m時(shí)刻從A點(diǎn)開始傳播的下行波場一致。進(jìn)一步可以得到如下關(guān)系：D點(diǎn)在t時(shí)刻接收到的總波場p可以視作C點(diǎn)在t+tsi,m時(shí)刻接收到的上行波場u1與A點(diǎn)在t-tsi,m時(shí)刻接收到的上行波場u3乘以海面反射系數(shù)之后求和,即

p(xi,t)=u1(xi+Δxi,m,t+tsi,m)+

Rmu3(xi-Δxi,m,t-tsi,m)

(15)

式(15)中：Rm為B點(diǎn)的海面反射系數(shù)。

假設(shè)海水彈性參數(shù)在橫向上均勻不變，則A、C兩點(diǎn)(即橫坐標(biāo)xi+Δxi,m與xi-Δxi,m處)所接收的上行波與B點(diǎn)接收到的上行波存在如下關(guān)系：

uB(xi,t-Δti,m)=uC(xi+Δxi,m,t)

(16)

uB(xi,t+Δti,m)=uC(xi-Δxi,m,t)

(17)

為了使公式簡潔，下文中u(xi,t)均指B點(diǎn)處的波場uB。

聯(lián)立式(15)～(17)可得

p(xi,t)=u(xi,t+tsi,m-

Δti,m)+Rmu(xi,t-tsi,m+Δti,m)

(18)

對式(18)進(jìn)行τ-p變換可得p(pm,τ)，再將p(pm,τ)變換到頻率域得p(pm,ω)，最后在頻率域進(jìn)行線性τ-p反變換得p(xi,ω)[33],即

-jω(pmxi+

exp[-jω(pmxi-Δti,m+tsi,m)]

(19)

式(19)中：U(pm,ω)為上行波場u(xi,t)進(jìn)行τ-p變換后再變換到頻率域的結(jié)果。

令

D1i,j=exp[-jω(pmxi+Δti,m-tsi,m)]

(20)

D2i,j=Rmexp[-jω(pmxi-Δti,m+tsi,m)]

(21)

Di,j=D1i,j+D2i,j

(22)

利用上式將式(19)寫成矩陣形式,可得

P=DU

(23)

其中

通過求解線性方程組(23)即可得到B點(diǎn)上行波波場U。

本文采用最小二乘殘差算法求解上行波波場U。將接收點(diǎn)波場P作為數(shù)據(jù)空間矢量，矩陣D作為模型空間，U作為解空間矢量。選擇正常數(shù)α1、β1使得數(shù)據(jù)空間正交基矢量c1和解空間正交基矢量ω1的范數(shù)為1，由于矩陣的雙對角化迭代過程須滿足[34]：P=β1c1，DTc1=α1ω1，則可采用如下迭代算法[33]建立雙對角矩陣元素與各空間正交基矢量的計(jì)算關(guān)系,即

βk+1ck+1=Dωk-αkck

(24)

αk+1ωk+1=DTck+1-βk+1ωk+1

(25)

k次迭代之后，有

(26)

其中

Wk=(ω1,ω2,…,ωk)

Ck=(c1,c2,…,ck)

式(26)中:ek+1代表單位矩陣中的第k+1列的向量。

令Lk+1=(Bk,αk+1ek+1)，結(jié)合式(26)可得

(27)

在每一步的迭代中計(jì)算殘差rk=P-DUk，假定方程的近似解為

Uk=Wkyk

(28)

其中yk與殘差存在如下關(guān)系，并代入雙對角化關(guān)系P=β1c1，DTc1=α1ω1以及式(27)可得

(29)

2.2 海面起伏影響的消除

利用式(24)～(29)求得的海面上行波場本質(zhì)是粗糙海面上各點(diǎn)的上行波場，需要將其校正至接收點(diǎn)所在的深度。本文采用波場延拓方法進(jìn)行校正,采用式(30)所示的算子[36]進(jìn)行波場延拓。

Pz+Δz=Pze±iθ

(30)

3 模型數(shù)據(jù)測試

3.1 粗糙海面層狀模型

利用圖3所示的層狀模型驗(yàn)證本文方法的有效性，其中海面的起伏量如圖4所示。利用有限差分技術(shù)對圖3進(jìn)行正演模擬，其中模擬子波為主頻35 Hz的雷克子波，震源置于海面，記錄道數(shù)為240道，道間距為12.5 m，記錄長度為0.8 s。變深度纜形態(tài)如圖5，沉放深度在10～20 m逐步加深，變深纜處采用亞網(wǎng)格技術(shù)[37]進(jìn)行空間和時(shí)間的加密采樣，粗糙海面的處理采用雙變網(wǎng)格技術(shù)[38]進(jìn)行。

圖3 層狀模型

圖4 層狀模型頂端一維粗糙海面

圖5 變深度纜形態(tài)

圖6a與圖6b分別為粗糙海面和水平海面條件下正演記錄。對比可看出，圖6a中在接收到海底的一次反射波之后，由于海面起伏不平，一次反射波上行到海面并被粗糙海面反射后，記錄中出現(xiàn)了大量的同相軸斷裂、不連續(xù)等現(xiàn)象，信噪比也較低。利用本文方法進(jìn)行鬼波壓制和海面起伏量校正后的結(jié)果見圖6c，圖6d為圖6b的鬼波壓制結(jié)果。圖6c與圖6d中存在海面下行波導(dǎo)致的多次波干擾，本文利用鬼波壓制過程中的下行波分量進(jìn)行延拓模擬出多次波能量，并從記錄中減去，從而只保留有效反射波的能量記錄。對比表明，本文方法可以取得良好的鬼波壓制效果，消除了粗糙海面反射的下行波帶來的各種干擾。圖6c中存在的少量竄擾噪聲是由于海面的成像精度不足導(dǎo)致的，更為有效的海面成像方法正在研究之中。

圖6 粗糙及水平海面鬼波壓制前后地震記錄

圖7為圖6a與圖6c的振幅譜。從圖7a中可以觀察到明顯的陷波效應(yīng)的存在(如紅色箭頭所示)，而圖7b所示的壓制鬼波后的振幅譜，由于消除了鬼波引起的陷波效應(yīng)，各道頻譜的連續(xù)性增強(qiáng)，并且補(bǔ)償了低頻成分，表明本文方法可以有效壓制粗糙海面的鬼波反射。

圖7 去鬼波前后地震記錄頻譜

3.2 MarmousiII模型

采用二維MarmousiII模型數(shù)據(jù)進(jìn)行試驗(yàn)，模型數(shù)據(jù)由有限差分計(jì)算得到，共50炮，每炮240道接收，接收纜形態(tài)與圖5一致，正演所用的其他參數(shù)與方法和3.1節(jié)相同，模型上部海面的起伏量如圖8所示。

圖8 MarmousiII模型頂端粗糙海面

圖9a與圖9c分別是含鬼波炮集的逆時(shí)偏移剖面與共偏移距剖面(1 000 m偏移距)，圖9b與圖9d分別為采用本文方法壓制鬼波后的逆時(shí)偏移剖面與共偏移距剖面(1 000 m偏移距)，圖9e與圖9f分別為圖9a與圖9b的局部放大，對比發(fā)現(xiàn)，進(jìn)行鬼波壓制前的深度偏移剖面，尤其是淺部與低速體附近分辨率較低，存在由鬼波形成的虛反射成像界面，信噪比遠(yuǎn)低于鬼波壓制后的深度剖面。由于模型中部分反射層位相鄰很近，這些反射層位的鬼波同相軸與相鄰層位的反射波形成的同相軸位置非常接近甚至疊加在一起，因此在圖9a中可見一些相鄰連續(xù)出現(xiàn)的反射層以及圖9b中原本較弱的反射界面在圖9a中能量變強(qiáng)。由于同時(shí)存在多次波的干擾，為更好地確認(rèn)鬼波是否被正確壓制，從所繪制的去鬼波前后的共偏移距剖面(1 000 m偏移距)可以觀察到，圖9c中箭頭①、④處為海底、低速體處的鬼波成像結(jié)果，但在圖9d中得到了良好的壓制；圖9c中箭頭②、⑤處分辨率很低，無法分辨出成像的層位，而鬼波壓制后在圖9d中相應(yīng)位置能看到較為清晰的成像結(jié)果；圖9d中箭頭③處成像與圖9c中相比很弱，因此判斷此處為下行波產(chǎn)生的多次波能量；圖9c中箭頭⑥處約1 600 m和1 800 m兩個(gè)成像界面之間可以觀察到1 600 m處界面鬼波的成像，但在圖9d中得到了壓制，分辨率也得到了提高。由于海底反射層位分布較為復(fù)雜，偏移結(jié)果中多次波和鬼波的成像難以直接進(jìn)行區(qū)分，但依據(jù)海底反射和其鬼波成像的相對位置關(guān)系，配合鬼波壓制前后圖件的對比，可對多次波和鬼波進(jìn)行一個(gè)大致的區(qū)分。圖9c、d中箭頭⑦、⑧所指部分與圖9a、b中紅框處放大區(qū)域(圖9e、f)中箭頭標(biāo)注部分是同一部分，在共偏移距剖面中可以更清晰地看到，箭頭⑧所指部分明顯遠(yuǎn)離上層能量較強(qiáng)的成像界面，因此判斷為多次波；而對于箭頭⑦處，相對于上層界面位置距離與海底和其鬼波的成像位置關(guān)系較為符合，同時(shí)其成像能量的橫向強(qiáng)弱變化也與相鄰上層的變化一致，而多次波則不應(yīng)該存在這樣的關(guān)系，因此判斷箭頭⑦處是鬼波成像。

圖9 去鬼波前后地震剖面

4 結(jié)論

1) 本文方法使用蒙特卡洛原理對粗糙海面建模，并求取對應(yīng)的海面反射系數(shù)，更符合海洋拖纜采集時(shí)的海面條件，具有實(shí)際應(yīng)用價(jià)值。

2) 本文方法使用數(shù)學(xué)原理上更穩(wěn)定、更精確、收斂速度更快的最小二乘殘差算法，在τ-p域?qū)ψ兩疃壤|數(shù)據(jù)進(jìn)行波場分離，從而有效壓制鬼波，消除地震數(shù)據(jù)頻譜中的陷波效應(yīng)。

3) 不同模型的合成數(shù)據(jù)測試結(jié)果表明，本文方法能夠在粗糙海面條件下對接收點(diǎn)鬼波波場實(shí)現(xiàn)準(zhǔn)確、有效壓制。