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

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

(1. 中國海洋大學 山東青島 266100； 2. 青島海洋科學與技術國家實驗室海洋礦產資源評價與探測技術功能實驗室 山東青島 266071；3. 海底科學與探測技術教育部重點實驗室 山東青島 266100； 4. 中國地質調查局青島海洋地質研究所 山東青島 266100)

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

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

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

1 粗糙海面建模與反射系數計算

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

一般來說，某個時刻的海面是在一定的高差范圍內隨機起伏的，且起伏量滿足高斯分布[21]，因此可利用統計學原理對粗糙海面進行建模。

1.1 描述粗糙海面的統計參量

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

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

(1)

相關度函數C定義為

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

(2)

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

利用海面均方根高度將相關度函數歸一化,可得到相關系數ρ

(3)

ρ(r)=1/e時所對應的r值即為相關長度，記為l。相關長度是判斷粗糙海面上任意兩點是否相互獨立的標準之一。

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

(4)

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

(5)

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

(6)

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

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

(7)

1.2 粗糙海面成像與反射系數計算

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

(8)

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

(9)

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

(10)

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

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

圖1 粗糙海面建模結果

在短時間尺度內，采用低階小斜率近似方法[29-30]建立時變海面復反射系數的計算模型，并由此求出對應的粗糙海面反射系數[31]為

(11)

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

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

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

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

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

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

(12)

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

圖2 接收點波場幾何關系示意圖

設u1到達C點與u2到達B點的時差為Δti,m，則有

(13)

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

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

(14)

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

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

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

(15)

式(15)中：Rm為B點的海面反射系數。

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

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點處的波場uB。

聯立式(15)～(17)可得

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

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

(18)

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

-jω(pmxi+

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

(19)

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

令

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點上行波波場U。

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

β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)，結合式(26)可得

(27)

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

Uk=Wkyk

(28)

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

(29)

2.2 海面起伏影響的消除

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

Pz+Δz=Pze±iθ

(30)

3 模型數據測試

3.1 粗糙海面層狀模型

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

圖3 層狀模型

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

圖5 變深度纜形態

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

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

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

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

3.2 MarmousiII模型

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

圖8 MarmousiII模型頂端粗糙海面

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

圖9 去鬼波前后地震剖面

4 結論

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

2) 本文方法使用數學原理上更穩定、更精確、收斂速度更快的最小二乘殘差算法，在τ-p域對變深度纜數據進行波場分離，從而有效壓制鬼波，消除地震數據頻譜中的陷波效應。

3) 不同模型的合成數據測試結果表明，本文方法能夠在粗糙海面條件下對接收點鬼波波場實現準確、有效壓制。