基于分治算法的谷值因子斜纜數(shù)據(jù)鬼波壓制方法*

崔明曉 王守東 王豆豆 齊子威

(1.中國石油大學(北京)油氣資源與探測國家重點實驗室 北京 102249； 2.中國石油大學(北京)海洋石油勘探國家工程實驗室 北京 102249 )

海洋地震數(shù)據(jù)拖纜采集時通常是將采集拖纜放置于海水面以下，由于海水面的反射，檢波器除了接收到地下反射波外，還會接收到反射波到達海水面后向下產(chǎn)生的虛反射，該虛反射也被稱為鬼波[1]。由于鬼波陷波效應的影響，使得地震數(shù)據(jù)頻帶變窄，從而降低了地震剖面的分辨率[2]，增加了地震地質(zhì)解釋的困難，因此必須在地震資料采集和處理時設法消除鬼波的影響[3]。以往常規(guī)海洋地震數(shù)據(jù)采集一般采用等深度纜采集技術[4]。為了壓制鬼波，國際上先后發(fā)展了OBC[5]、上下纜[6-7]和雙檢等采集方式[8]，提出了上下纜采集數(shù)據(jù)處理方法[9]和雙檢采集數(shù)據(jù)處理方法[10]，如PGS公司提出了雙傳感器采集系統(tǒng)中上行波場分離的方法壓制虛反射[11],這些方法對于消除鬼波影響，保持地震數(shù)據(jù)的低頻和高頻信息都有較好的效果。

近幾年，CGG公司推出了變深度纜采集技術[12]。針對變深度纜采集數(shù)據(jù)，Sablon等[13]提出了鏡像偏移壓制虛反射的方法，該方法能夠有效壓制虛反射和拓寬地震頻帶[14]；Soubaras等[15]提出了聯(lián)合反褶積壓制虛反射的方法，該方法在墨西哥灣取得了良好的應用效果；Wang等[16]提出了一種在F-XY域利用實際炮集和鏡像炮集聯(lián)合反褶積來壓制鬼波的算法；許自強 等[17]提出了海上變深度纜數(shù)據(jù)最優(yōu)化壓制虛反射的方法；葉林 等[18]提出了峰值因子最大化壓制鬼波方法。本文推導了利用鏡像變深度纜地震數(shù)據(jù)聯(lián)合反褶積去鬼波的算法，針對該方法中因鬼波延遲時間不準確而導致不能得到最優(yōu)壓制效果的問題，提出了通過谷值因子來判斷最精確鬼波延遲時間的技術思路，并進行了模型試算及實際數(shù)據(jù)處理，均取得了較好的去鬼波效果。

1 方法原理

1.1 數(shù)據(jù)最優(yōu)化聯(lián)合反褶積去鬼波方法

對于沉放深度為zi處的變深度纜所接收到的實際波場p(x,zi,t)，可以分解為上行波場u(x,zi,t)和下行波場d(x,zi,t)之和的形式，見式(1)，其中上行波場為地震波自激發(fā)后向下傳播經(jīng)地下界面反射直接到達檢波器的波場，下行波場則是地下反射波先到達海水面后經(jīng)海水面反射再到達檢波器的波場。

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

(1)

式(1)中：x為偏移距變量;t為時間變量。

把式(1)變換到頻率域,可得

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

(2)

波場傳播示意圖見圖1，其中深度為z1=zr處的實際變深度纜某接收點接收到的總波場為p(x,z1,t)，上行波場為u(x,z1,t)，下行波場為d(x,z1,t)；首先在深度為2zr的位置定義一個假想的海水界面，然后放置一個關于此界面與實際接收點對稱的鏡像接收點，其深度為z2，且z2=3zr，則假設鏡像接收點接收到的總波場為p(x,z2,t)，上行波場為u(x,z2,t)，下行波場為d(x,z2,t)。由式(1)、(2)可得

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

(3)

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

(4)

同時使得鏡像接收點接收到的下行波場等于實際拖纜接收到的上行波場，即

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

(5)

圖1 波場傳播示意圖Fig.1 Schematic of wave propagation

通常情況下，我們可以把海水面的反射系數(shù)取值為-1，此時就可以把下行波場寫成上行波場和波場延拓算子的形式,即

(6)

將式(5)代入式(6)，可得

U(x,z2,ω)=-U(x,z1,ω)ejk(z2-z1 )

(7)

將式(5)～(7)式代入式(3)、(4)，可得

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

(8)

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

(9)

從式(8)、(9)可以看出，實際拖纜和鏡像接收點處接收到的總波場都可以用實際拖纜處接收到的上行波場和一個算子計算得到，聯(lián)立式(8)、(9)，可以得出鏡像記錄是在實際原始記錄的基礎上作用了一個算子S(x,ω)，即

P(x,z2,ω)=P(x,z1,ω)S(x,ω)

(10)

其中，鏡像算子S(x,ω)為

S(x,ω)=-ejk(z2-z1)

(11)

因為原始記錄和鏡像記錄又可以用實際拖纜z1處的上行波場U(x,z1,ω)和鬼波算子G(x,ω)表示，即

(12)

其中

(13)

式(13)中：Δt為上行波與下行波之間的延遲時間。從而頻率空間域的鏡像記錄可以表示為

P(x,z2,ω)=P(x,z1,ω)(-ejωΔt)

(14)

這樣即可推導出用原始記錄和帶有鬼波延遲時間Δt的算子來求取鏡像記錄的公式。

利用式(14)可以由原始地震記錄計算出鏡像地震記錄，壓制鬼波可以通過原始記錄與鏡像記錄聯(lián)合反褶積法來實現(xiàn)，具體方法可以通過求解如下最小平方問題來實現(xiàn)，即

U(x,ω)=

(15)

通過求解式(15)可以得到上行波場頻率域數(shù)據(jù)U(x,ω)，再通過反傅里葉變換即可得到時間域上行波場u(x,t)，也就是壓制鬼波后的地震數(shù)據(jù)。

1.2 谷值因子提出

以上推導中，鏡像記錄求解時需要已知鬼波延遲時間Δt，但實際中Δt的影響因素有很多，包括檢波器深度、海水面、地震波傳播角度等，因此，往往不能事先得出精確的Δt，使得計算得到的鏡像記錄不夠準確，進而導致通過數(shù)據(jù)聯(lián)合反褶積方法得不到最優(yōu)的鬼波壓制效果。為了求得最精確的鬼波延遲時間Δt，本文采用了谷值因子的思路。

當鬼波延遲時間不準確時，原始記錄和鏡像記錄做單道反褶積得到的波形會產(chǎn)生震蕩，與精確結果相比波形復雜。為了刻畫波形的形態(tài)，用得到的單道地震數(shù)據(jù)最大值和最小值之差和該道記錄的均方根作商，假設代入在第i個鬼波延遲時間求得的鬼波數(shù)據(jù)為di(x,z,t)，由此定義了1個參數(shù)V的計算方法，即

(16)

其中，rms[di(x,z,t)]為該下行波地震道數(shù)據(jù)的均方根，則代入的每一個Δt均可得到一個V，把所有得到的V值中最小的一個稱為谷值因子。谷值因子所對應的Δt即為該地震道最精確的鬼波延遲時間，用該延遲時間得到的去鬼波數(shù)據(jù)是最精確的。因此,搜索最精確的鬼波延遲時間Δt的過程即為搜索谷值因子的過程。谷值因子可通過建立如下目標函數(shù)來確定,即

s=minV

(17)

1.3 分治算法搜索谷值因子

根據(jù)以上方法，通過搜索谷值因子來確定最精確鬼波延遲時間。為了提高搜索速度，本文采用分治算法搜索，其思想是：當在搜索一組數(shù)據(jù)中的最小值時，如果參與搜索的數(shù)據(jù)非常多，根本無法直接找出其中的最小值，使得計算過程非常復雜，可以把參與搜索的數(shù)據(jù)分成幾組，即把一個大的問題分解成幾個子問題[19]，如果這些子問題可以直接得出結果，也就是分解產(chǎn)生的數(shù)據(jù)能夠直接判斷最小值，則問題就很容易解決了；如果被劃分成的每個數(shù)據(jù)組依然不容易判斷最小值，可以把子問題繼續(xù)劃分為更小的子問題，以此類推繼續(xù)劃分數(shù)據(jù)，直至可以判斷最小值，并得出谷值因子為止。

通過分治算法可以有效搜索出谷值因子，得出最精確的鬼波延遲時間，從而得出最優(yōu)的去鬼波結果?；诜种嗡惴ㄋ阉鞴戎狄蜃觼硐聿ǖ倪^程可以分成以下4個步驟：

1) 確定鬼波延遲時間的搜索范圍和步長，把該范圍內(nèi)的每一個鬼波延遲時間代入式(14)，計算對應的鏡像數(shù)據(jù)。

2) 由原始數(shù)據(jù)和鏡像數(shù)據(jù)，應用數(shù)據(jù)最優(yōu)化聯(lián)合反褶積方法求得上行波場和下行波場，利用時間域下行波場分別計算對應的參數(shù)V。

3) 在得到的所有參數(shù)V中通過分治算法快速尋找最小值Vmin，即為谷值因子，其對應的Δt即為最精確的鬼波延遲時間。

4) 根據(jù)最精確的鬼波延遲時間計算最準確的鏡像記錄，進而得到最精確的去鬼波地震數(shù)據(jù)。

2 模型試算

2.1 驗證谷值因子的有效性

為了驗證谷值因子對應著最佳鬼波延遲時間的正確性，在單層介質(zhì)的情況下給定一個主頻為50 Hz的單道信號(圖2a)，再以50 ms為標準延遲時間給定一個鬼波信號(圖2b)，則2個信號相加即可認為是帶有鬼波的原始信號(圖2c)。設置鬼波延遲時間的搜索范圍Δt∈(25 ms，75 ms)，搜索步長Δτ=1 ms，則可知在此范圍內(nèi)共可搜索50次，分別得到Δt從25 ms到75 ms的鬼波記錄和每一個鬼波波記錄對應的參數(shù)V。圖2d、e、f、g、h分別為Δt=30、45、50、55、75 ms時得到的鬼波信號，可以看出圖2f即Δt=50 ms時的鬼波信號與真實的鬼波信號最接近。

把搜索50次中每次得到的V值計算出來，并繪制變化曲線圖(圖3)，發(fā)現(xiàn)當Δt=50 ms，即Δt為精確的鬼波延遲時間時所對應的V值最小，即為谷值因子，這時提取的鬼波正好最接近真實鬼波。

圖2單層介質(zhì)某單道信號不同鬼波延遲時間得到的鬼波波形

Fig.2Ghostwaveformwithdifferentghostwavedelayofsinglechannelsignalinsinglelayermedia

圖3 單層介質(zhì)參數(shù)V值隨鬼波延遲時間Δt的變化曲線Fig.3 Curve of the V change with the ghost wave delay Δt in single layer media

以上結果驗證了谷值因子在單層介質(zhì)信號上的有效性。再對多層介質(zhì)的情況進行驗證，給出1個三層介質(zhì)下的一次波某單道信號(圖4a)，然后再以鬼波延遲時間為50 ms得到鬼波信號(圖4b)，2個信號相加即可得到帶有鬼波的原始信號(圖4c)。與單道信號一樣，在搜索范圍Δt∈(25 ms，75 ms)以步長為1 ms依次代入其中的每一個鬼波延遲時間，分別得到下行波場，并由其計算每一個對應的參數(shù)V。圖4d、e、f、g、h分別為Δt=30、45、50、55、75 ms時得到的鬼波信號，可以看出圖4f即Δt=50 ms時的鬼波信號與真實的鬼波信號最接近。

圖4 多層介質(zhì)某單道信號不同鬼波延遲時間得到的鬼波波形Fig.4 Ghost waveform with different ghost wave delay of single channel signal in multilayer media

把每次搜索對應的參數(shù)V值隨延遲時間Δt變化曲線畫出來(圖5)，可以看出多層介質(zhì)情況下最精確的鬼波延遲時間對應的參數(shù)V也是最小的，即多層介質(zhì)中谷值因子同樣可以作為最優(yōu)鬼波延遲時間Δt的判斷標準。

圖5 多層介質(zhì)參數(shù)V值隨鬼波延遲時間Δt變化曲線Fig.5 Curve of the V change with the ghost wave delay Δt in multilayer media

由此可見，谷值因子所對應的鬼波延遲時間Δt即是最精確的鬼波延遲時間，因此，為了得到最優(yōu)化的去鬼波結果，計算搜索谷值因子是一種比較有效的方法。

2.2 簡單模型試算

建立的簡單速度模型如圖6a所示。橫向采樣間隔1.5 m，縱向采樣間隔1.5 m，水層速度1 500 m/s，采用70 Hz雷克子波。聲波方程有限差分正演模擬，單邊接收，每炮645道。震源放于左側(cè)水面上，檢波器深度4.5～50.0 m，圖6b為纜深變化，圖6c為第100道V值大小隨鬼波延遲時間的變化曲線。圖7a為帶有鬼波的原始地震記錄，圖7b為對應的鏡像地震記錄，圖7c為運用本文方法去除鬼波后的地震記錄。圖8a、b分別為對應的去鬼波前、后的頻譜圖。由簡單模型的地震記錄和頻譜圖可以看出，運用本文方法對帶有鬼波的原始記錄進行去鬼波處理，在剖面上可以有效的壓制鬼波，提高地震數(shù)據(jù)的分辨率；在頻率域可以消除陷波效應，提高頻帶寬度。

圖6 簡單速度模型(a)、纜深變化(b)和第100道V值隨鬼波延遲時間Δt變化曲線(c)Fig.6 Simple speed model(a),cable depth change(b)and the curve of V change with the ghost wave delay Δt in the 100th trace(c)

圖7 簡單模型原始地震記錄局部(a)、鏡像記錄局部(b)和去鬼波后地震記錄局部(c)Fig.7 Part of original seismograms(a),mirror seismograms(b) and seismograms after deghost(c)of the simple model

圖8 簡單模型處理前(a)與處理后(b)頻譜圖對比Fig.8 Spectrum before(a)and after(b)processing of the simple model

2.3 復雜模型試算

建立的復雜速度模型如圖9a所示。橫向采樣間隔1.5 m，縱向采樣間隔1.5 m，水層速度1 500 m/s，采用70 Hz雷克子波。聲波方程有限差分正演模擬，單邊接收，每炮1 345道。震源放于左側(cè)水面上，檢波器深度4.5～50.0 m，圖9b為纜深變化，圖9c為第150道V值隨鬼波延遲時間的變化曲線。圖10a為帶有鬼波的原始地震記錄，圖10b為最優(yōu)鬼波延遲時間計算得到的鏡像地震記錄，圖10c為運用本文方法去除鬼波后的地震記錄。圖11a、b分別為對應的去鬼波前、后的頻譜圖。由復雜模型的地震記錄和頻譜圖可以看出，本文方法對去除復雜模型地震記錄的鬼波和消除陷波效應也具有良好的效果。

圖9 復雜速度模型(a)、纜深變化(b)和第150道V值隨鬼波延遲時間Δt變化曲線(c)Fig.9 Complex speed model(a),cable depth change(b)and the curve of V change with the ghost wave delay Δt in the 150th trace(c)

圖10 復雜模型原始地震記錄鏡局部(a)、鏡像記錄局部(b)和去鬼波后地震記錄局部(c)Fig.10 Part of original seismograms(a),mirror seismograms(b)and seismograms after deghost(c)of the complex model

圖11 復雜模型處理前(a)與處理后(b)頻譜圖對比Fig.11 Spectrum before(a)and after(b)processing of the complex model

3 實際數(shù)據(jù)處理

用本文方法對某海上地震數(shù)據(jù)進行處理。圖12a為原始地震記錄，圖12b為鏡像地震記錄，圖12c為用本文方法進行處理后得到的地震剖面，圖13a、b、c分別為實際資料相應的局部顯示。由實際地震數(shù)據(jù)的處理結果可以看出，本文壓制鬼波方法取得了良好的效果。

圖12 實際資料原始地震記錄(a)、鏡像地震記錄(b)和壓制鬼波后地震記錄(c)Fig.12 Original seismograms(a),mirror seismograms(b) and seismograms after deghost(c)of the actual data

圖13 實際資料原始局部顯示(a)、鏡像局部顯示(b)和去鬼波后局部顯示(c)Fig.13 Part of original seismograms(a),mirror seismograms(b) and seismograms after deghost(c)of the actual data

4 結論

1) 在應用數(shù)據(jù)最優(yōu)化聯(lián)合反褶積方法壓制斜纜數(shù)據(jù)鬼波時，通過搜索谷值因子可以得到最精確的鬼波延遲時間，從而解決鬼波延遲時間不準確的問題，進而得到更優(yōu)化的斜纜數(shù)據(jù)鬼波壓制效果。

2) 在搜索谷值因子時，分治算法可以作為一種有效的搜索方法，通過該方法可以減少搜索次數(shù)，并準確地找到谷值因子，進而確定鬼波延遲時間，得到去鬼波的地震記錄。

3) 應用本文方法對模擬數(shù)據(jù)和實際數(shù)據(jù)進行處理均得到了較好的去鬼波效果，在剖面上提高了地震資料的分辨率，拓寬了地震頻帶。

[1] 陸敬安，伍忠良，曾憲軍.海洋地震勘探中地震波、鬼波綜合效應分析與應用[J].海洋技術，2006,25(4):76-78，98.

LU Jing’an,WU Zhongliang,ZENG Xianjun.The synthesized effect of seismic wave and ghost reflection and its application in marine seismic survey[J].Ocean Technology,2006,25(4):76-78，98.

[2] 李套山，高洪強，劉振夏，等.虛反射信息的采集及其形成機制、頻率響應的理論討論[J].石油物探，1997,36(4)：38-44.

LI Taoshan,GAO Hongqiang,LIU Zhenxia,et al.The acquisition of ghost information and the study of its generation and frequency response[J].GPP,1997,36(4):38-44.

[3] 劉建磊，王修田.海上地震虛反射相位剔除法反褶積[J].海洋地質(zhì)與第四紀地質(zhì)，2002,22(4):117-121.

LIU Jianlei,WANG Xiutian.Deconvolution by phase spectrum deghosting in marine seismic survey[J].Marine Geology & Quaternary Geology,2002,22(4):117-121.

[4] 王振峰，李緒深，孫志鵬，等．瓊東南盆地深水區(qū)油氣成藏條件和勘探潛力[J].中國海上油氣，2011,23(1):7-13.

WANG Zhenfeng,LI Xushen,SUN Zhipeng,et al.Hydrocarbon accumulation conditions and exploration potential in the deep-water region，Qiongdongnan basin[J].China Offshore Oil and Gas,2011,23(1):7-13.

[5] 王守君.海底電纜地震技術優(yōu)勢及在中國近海的應用效果[J].中國海上油氣，2012,24(2):9-12.

WANG Shoujun.Technical advantages of OBC seismic survey and its application effects offshore China[J].China Offshore Oil and Gas,2012,24(2):9-12.

[6] HILL D,COMBED L,BACON J.Over/under acquisition and data processing:the next quantum leap in seismic technology[J].First Break,2006,24(6):81-95.

[7] MOLDOVEANU N,COMBEE J.Over/under towed-streamer acquisition:a method to extend seismic bandwidth to both higher and lower frequencies[J].The Leading Edge,2007,26(1):41-58.

[8] TENGHAMN R,VAAGE S,BORRESEN C.A dual-sensor,towed marine streamer:its viable implementation and initial results[J].SEG Expanded Abstracts,2007(1):989.

[9] 劉春成，劉志斌，顧漢明，等.利用上下纜合并算子確定海上上下纜采集的最優(yōu)沉放深度組合[J].石油物探，2013,52(6)：623-629.

LIU Chuncheng,LIU Zhibin,GU Hanming,et al.The determination of optimal sinking depths of over/under streamers in offshore survey by merge operator[J].GPP,2013,52(6):623-629.

[10] SOLLNER W,BROX E,WIDMAIER M,et al.Surface related multiple suppression in dual senor towed steamer data[J].SEG Technical Program Expanded Abstracts,2007,26:2540-2544.

[11] CARLSON C,LONG A,SOLLNER W,et al.Increased resolution and penetration from a towed dual-sensor streamer[J].First Break,2007,25(12):71-77.

[12] SOUBARAS R,DOWLE R.Variable-depth streamer,a broadband marine solution[J].First Break,2010,28:89-96.

[13] SABLON R,RUSSIER D,ZURITA O,et al.Multiple attenuation for variable-depth streamer data:from deep to shallow water[C].81st Annual Meeting,SEG,Expanded Abstracts,2011:3505-3509.

[14] SOUBARAS R.Pre-stack deghosting for variable-depth streamer data[C].82nd Annual Meeting,SEG,Expanded Abstracts,2012:780-789.

[15] SOUBARAS R.Deghosting by joint deconvolution of a migration and a mirror migration[C].80th Annual Meeting,SEG,Expanded Abstracts,2010:3406-3410.

[16] WANG P,PENG C.Premigration deghosting for marine streamer data using a bootstrap approach[C].82nd Annual Meeting,SEG,Expanded Abstracts,2012:1-5.

[17] 許自強，方中于，顧漢明，等.海上變深度纜數(shù)據(jù)最優(yōu)化壓制鬼波方法及其應用[J].石油物探，2015,54(4):404-413.

XU Ziqiang,FANG Zhongyu,GU Hanming,et al,The application of optimal deghosting algorithm on marine variable-depth streamer data[J].GPP,2015,54(4):404-413.

[18] 葉林,韓立國.基于峰值因子最大化的鬼波壓制方法[J].世界地質(zhì),2016,35(4):1102-1107.

YE Lin,HAN Liguo.Deghosting method based on maximization of peak factor[J].Global Geology,2016,35(4):1102-1107.

[19] 李紹華，王建新，馬振宇，等.基于加權分治技術的setpacking精確算法[J].小型微型計算機系統(tǒng)，2010,31(6)：1180-1184.

LI Shaohua,WANG Jianxin,MA Zhenyu,et al.Exact algorithm for set packing based on measure and conquer[J].Chinese Computer Systems,2010,31(6):1180-1184.