王華忠， 馮波， 劉少勇， 胡江濤， 王雄文， 李輝

同濟(jì)大學(xué)海洋與地球科學(xué)學(xué)院 波現(xiàn)象與反演成像研究組， 上海 200092

?

疊前地震數(shù)據(jù)特征波場分解、偏移成像與層析反演

王華忠， 馮波*， 劉少勇， 胡江濤， 王雄文， 李輝

同濟(jì)大學(xué)海洋與地球科學(xué)學(xué)院 波現(xiàn)象與反演成像研究組， 上海 200092

本文提出了一套疊前地震數(shù)據(jù)稀疏表達(dá)(特征波場合成)、深度偏移成像和層析成像的處理流程.不同于傳統(tǒng)的變換域中的數(shù)據(jù)稀疏表達(dá)理論，本文利用局部平面波的傳播方向(慢度矢量)，在中心炮檢點(diǎn)處同時進(jìn)行波束合成，從而將地震數(shù)據(jù)投影到局部平面波域(高維空間)中.由于波束合成后的地震數(shù)據(jù)描述了局部平面波的方向特征，因此稱之為特征波場.然而波束合成算法需要估計局部平面波的慢度矢量.當(dāng)?shù)卣饠?shù)據(jù)受噪聲干擾時，難以在常規(guī)τ-p譜中自動估計局部平面波的射線參數(shù)(慢度矢量).本文提出了基于反演理論的特征波場合成方法，可以同時反演局部平面波及其傳播方向，從而提高特征波合成的自動化程度并保持方法的穩(wěn)健性.通過特征波場合成，可以將地震數(shù)據(jù)分解為單獨(dú)的震相(波形).這樣的數(shù)據(jù)可以直接用來成像及反演.在局部平面波域中，由于局部平面波的入射與出射射線參數(shù)已知，傳統(tǒng)的Kirchhoff疊前深度偏移(PSDM)和高斯束/控制束PSDM可以實(shí)現(xiàn)從“沿等時面的畫弧”到“向反射點(diǎn)(段)的直接投影”的轉(zhuǎn)變，疊前偏移的效率以及成像質(zhì)量可以同時提高.此外，特征波場與地下反射點(diǎn)(段)的一對一映射關(guān)系使得疊前深度偏移與層析成像融為一體，可以極大地提高速度反演的效率.數(shù)值試驗(yàn)證明了特征波場合成、疊前深度成像以及層析反演的有效性.

地震數(shù)據(jù)稀疏表達(dá); 特征波場合成； 立體數(shù)據(jù)空間； 局部平面波； 特征波場疊前深度偏移成像； 特征波場層析反演

1 引言

地震數(shù)據(jù)中含有豐富的有效信號，然而時空域的信號表達(dá)方式無法直接描述地震數(shù)據(jù)的特征.為此，可以用某種數(shù)學(xué)變換抓取地震信號的主要特征，實(shí)現(xiàn)地震數(shù)據(jù)的稀疏表達(dá).這有利于地震信號處理、波場分析和反演成像.事實(shí)上，地震數(shù)據(jù)可以在變換域中稀疏表達(dá)，常用的變換域有Fourier域、小波域、Curvelet域、Radon域、Gabor域等.廣義上講，可以通過選取一組正交基、或框架函數(shù)、甚至基函數(shù)字典，使得變換域中的地震數(shù)據(jù)具有稀疏性質(zhì).利用數(shù)據(jù)在變換域中的稀疏性，可以實(shí)現(xiàn)地震數(shù)據(jù)的壓縮與重構(gòu)插值等.Xu等(2005)在Fourier域反演非規(guī)則地震數(shù)據(jù)頻譜，實(shí)現(xiàn)數(shù)據(jù)重構(gòu).Herrmann和Hennenfent(2008)、白蘭淑等(2014)，利用地震數(shù)據(jù)在Curvelet域中的稀疏性，實(shí)現(xiàn)地震數(shù)據(jù)重建.

然而，Curvelet變換無法直觀描述地震波的傳播特征，因而沒有明確的地球物理含義.本文中，我們把疊前地震數(shù)據(jù)抽象為不同出射角度(即出射慢度矢量)以及不同到達(dá)時的局部平面波的線性疊加，而接收到的地震記錄可以看作局部平面波波前在空間離散的檢波器處的采樣.在立體層析理論(Billette和Lambaré，1998)中，地震記錄被簡化為炮點(diǎn)、檢波點(diǎn)坐標(biāo)，炮點(diǎn)、檢波點(diǎn)處的射線參數(shù)以及射線的雙程走時.上述參數(shù)構(gòu)成了立體層析中的數(shù)據(jù)空間，但忽略了地震子波的波形信息.

基于線性Radon變換(LRT)的局部平面波分解的方法是主要的地震信號的線性表達(dá)方法，以其明確的物理含義(可以描述局部平面波的傳播方向)而被廣泛應(yīng)用.非線性的Radon變換，如拋物Radon變換(Sacchi和Ulrych，1995)，雙曲Radon變換(Liu和Sacchi，2002)，多項(xiàng)式Radon變換(牛濱華等，2001)等，則主要用于多次波壓制、去噪、數(shù)據(jù)重構(gòu)等方面.不同于Fourier變換和正交小波變換，LRT的基函數(shù)是非完備正交的，因而導(dǎo)致LRT結(jié)果難以實(shí)現(xiàn)稀疏表達(dá).王雄文和王華忠(2014)將常規(guī)的LRT轉(zhuǎn)化為壓縮感知問題，在壓縮感知的理論框架下實(shí)現(xiàn)了高分辨率的LRT.在Radon域中，可以測量局部平面波在炮點(diǎn)處的入射射線參數(shù)Ps及檢波點(diǎn)處的出射射線參數(shù)Pr.若引入了局部平面波的方向信息，疊前地震記錄d=d(xs,xr,t)表現(xiàn)為立體地震數(shù)據(jù)，可描述為d=d(xs,ps,xr,pr,t).

本文中，首先考慮地震數(shù)據(jù)在局部平面波域中的稀疏表達(dá)，并引入地震數(shù)據(jù)的特征分解方法，隨后給出基于特征波數(shù)據(jù)的偏移成像與反演成像方法.

2 特征波場合成

在天然地震學(xué)中，通常用檢波器排列的某種組合來分離相干信號和噪聲.最基本的方法就是波束合成方法(BeamForming，Rost和Thomas，2002；Krüger等，1993)，它利用局部平面波波前到達(dá)不同檢波器的時差，對檢波器接收到的地震記錄進(jìn)行時移疊加，使局部平面波相干疊加而其他信號相互抵消，從而提取某種特定的波型(震相).

波束合成的思想也可應(yīng)用于勘探地震學(xué)中，如高斯束偏移方法(Hill，1990；Hale，1992)，在炮集中對不同檢波器中心進(jìn)行加窗局部傾斜疊加(即波束合成)，利用炮集的τ-p譜進(jìn)行偏移.在射線束偏移(Sun等，2001)和控制束偏移(Gao，2006)中，通過選擇τ-p譜的能量，實(shí)現(xiàn)有選擇性的波場傳播與成像.

波束合成算法需要估計局部平面波的方向信息.本文中，我們考慮Radon譜約束下的局部平面波反演方法(胡江濤等，2013)：

(1)

上述Radon譜約束下的局部平面波反演方法可以提高Radon譜的精度并壓制泄漏的噪聲.更進(jìn)一步地，為了得到更稀疏的局部平面波解(王雄文和王華忠，2014)，可以把式(1)改寫為壓縮感知問題，在L0范數(shù)意義下求解如下反問題：

(2)

利用匹配追蹤方法求解(2)式，可以得到局部空間窗內(nèi)的稀疏平面波解.

通過傳統(tǒng)的線性Radon變換或者高分辨率的平面波反演方法，可以得到局部平面波及其慢度矢量(或射線參數(shù)p).此時，我們可以考慮對炮集的波束合成.

在共炮點(diǎn)集中，對中心檢波點(diǎn)xr0的波束合成可以寫為

dBF(xr0,pr,t,xs)=

(3)

我們把在中心檢波點(diǎn)處進(jìn)行一次波束合成之后的地震數(shù)據(jù)記為dBF(xr0,pr,t,xs).顯然，經(jīng)過波束合成后，實(shí)現(xiàn)了出射方向?yàn)閜r的局部平面波的同相疊加，而其他信號相互干涉.

若將地震數(shù)據(jù)分選為共檢波點(diǎn)道集，此時可以再次在中心炮點(diǎn)xs0處進(jìn)行波束合成：

dDBF(xr0,pr,xs0,ps,t)=

(4)

同公式(3)，經(jīng)過中心炮點(diǎn)處的再次波束合成后，實(shí)現(xiàn)了入射方向?yàn)閜s的局部平面波的同相疊加.相應(yīng)的，我們把兩次波束合成之后的地震數(shù)據(jù)記為dDBF(xr0,pr,xs0,ps,t).其中，(xs0,xr0)為波束合成的中點(diǎn)坐標(biāo)，而(pr,ps)為局部平面波的慢度矢量.

將(3)式中的波束合成代入(4)式中，可得同時在中心炮點(diǎn)和檢波點(diǎn)處的波束合成算法：

(5)

在三維時，上述公式中的坐標(biāo)與射線參數(shù)均為矢量，表示為xs=(sx,sy)，xr=(gx,gy)，ps=(psx,psy)，pr=(prx,pry).經(jīng)過特征波場合成，五維地震數(shù)據(jù)d=d(xs,xr,t)被投影到九維的立體數(shù)據(jù)空間d=d(xs,ps,xr,pr,t)中.由于(5)式波束合成的地震數(shù)據(jù)可以描述局部平面波的傳播方向特征，因此我們稱(5)式合成的地震數(shù)據(jù)為特征波場.

從公式(5)可以看出，特征波場合成與局部平面波的入射與出射射線參數(shù)有關(guān)，其本質(zhì)是利用局部平面波的入射和出射慢度計算道間時差(Δt=p·Δx)，并將波束合成孔徑內(nèi)的地震道按照道間時差進(jìn)行時移疊加.在時移疊加的過程中，只有給定慢度矢量(局部平面波入射與出射射線參數(shù))的局部平面波波前實(shí)現(xiàn)了同相疊加，而其他的信號相互抵消，從而實(shí)現(xiàn)不同震相(波型)的分離，即特征波場合成.

事實(shí)上，上述方法也可以推廣到任意線性排列的觀測系統(tǒng)中(如VSP或井間觀測系統(tǒng)等)，只要滿足局部的炮點(diǎn)/檢波點(diǎn)排列具有線性特征.相應(yīng)的，局部平面波的慢度矢量(射線參數(shù))定義為地震波走時對地震排列(炮點(diǎn)或檢波點(diǎn)坐標(biāo))的方向?qū)?shù).

此外，在頻率域特征波場合成公式(5)對應(yīng)的相移矩陣為Ai,j=exp(iωpj·(xi-x0))，其中，射線參數(shù)向量p=(ps,pr)，炮檢坐標(biāo)向量x=(xs,xr).此時，基于反演理論的特征波場合成方法可以在L2或L0范數(shù)意義下求解(形式與公式(1)和(2)相同)，無需預(yù)先估計射線參數(shù)再進(jìn)行特征波場合成.

3 特征波場偏移成像

地震數(shù)據(jù)偏移成像的過程可以描述為地表波場的反傳播加上成像條件.以運(yùn)動學(xué)射線追蹤為傳播算子的Kirchhoff疊前深度偏移(PSDM)技術(shù)在橫向變速劇烈情況下存在多路徑和振幅焦散等問題，基于射線束的成像方法(劉少勇等，2012)在一定程度上解決了這些問題.高斯束偏移技術(shù)(Hill，1990；Hale，1992)已經(jīng)成為工業(yè)界一個重要的疊前深度偏移成像工具.Sun和Schuster(2001)提出菲涅爾帶控制的波路徑偏移方法，將積分法偏移在成像空間中沿等時面的投影轉(zhuǎn)化為沿胖射線路徑的投影，在計算效率上有一定優(yōu)勢.Liu和Palacharla(2011)對高斯束偏移進(jìn)行簡化，只保留其運(yùn)動學(xué)特征，推導(dǎo)出簡單的解析表達(dá)式來控制射線束的寬度和振幅，理論上說該實(shí)現(xiàn)方式效率要優(yōu)于高斯束偏移，成像精度基本有保障.快速射線束(Fast-beam)偏移(Gao等，2006)便是該思想的一種實(shí)現(xiàn).特征波偏移成像利用射線束波傳播算子來反傳局部特征波場，既保留了射線類成像的靈活性，又保證了成像精度.

特征波場是地震波場在高維空間中的稀疏表達(dá)，且局部平面波的入射與出射慢度矢量已知.利用這個性質(zhì)，射線束類的偏移方法可以實(shí)現(xiàn)從“畫弧”向“搬家”的轉(zhuǎn)變，進(jìn)一步提高成像效率.與上節(jié)的特征波場合成方法結(jié)合，匹配合適的局部射線束傳播算子，可以實(shí)現(xiàn)特征波疊前深度偏移成像(characteristicwaveimaging，簡記為CWI).更進(jìn)一步的，“搬家”或一對一映射的偏移方法與基于角度道集的層析速度反演可以自然的融合成一體，共享同樣的數(shù)據(jù)空間和成像空間的數(shù)據(jù)體.成像與層析速度反演一體化也是特征波成像的一大優(yōu)勢.

傳統(tǒng)基于炮道集的射線束(包括高斯束)成像條件可以描述為

I(x)= ∫Ω∫ps∫prWBF(x)DBF(τ,pr;xs(ps),xr)

×dprdpsdΩ|τ=τs+τr.

(6)

其中：x表示成像點(diǎn)坐標(biāo)，WBF表示由動力學(xué)射線追蹤得到的振幅加權(quán)函數(shù)，W表示成像孔徑，即對地下成像點(diǎn)有貢獻(xiàn)的地表炮檢點(diǎn)范圍.類似的，特征波成像條件可以表示為

I(x)=∫ΩWDBF(x)DDBF(τ,ps,pr;xr,xs)dΩ|τ=τs+τr.

(7)

其中，WDBF表示一個更為簡單的振幅加權(quán)函數(shù)，它是成像空間坐標(biāo)的函數(shù)，在不同的成像方法中可有不同的表達(dá).不論在傳統(tǒng)射線束成像還是特征波成像過程中，由于可以方便的根據(jù)走時梯度計算射線參數(shù)，進(jìn)而通過射線參數(shù)求出張角(三維時是張角和方位角).張角θ可以由下式表示：

(8)

傳統(tǒng)成像條件和特征波成像條件的幾何解釋可由圖1表示.圖中紅色反射層為地層真實(shí)位置，傳統(tǒng)射線束成像方法是將局部平面波投影到整個橢圓(三維情況為橢球)上.由于特征波場在地表的入射與出射慢度已知，因而特征波成像只對特征波場進(jìn)行“能量搬家”.紅色小橢圓的位置代表一對特征波射線束在成像域的菲涅爾帶范圍.從圖1可以看出，基于特征波場表達(dá)，實(shí)現(xiàn)了成像方式從“畫弧”到“搬家”的轉(zhuǎn)變.傳統(tǒng)積分類成像的畫弧成像過程必然會帶來畫弧噪聲，而搬家的方式可最大限度減少此類噪聲同時能兼顧菲涅爾帶的成像分辨率，其成像效率也會得到很大的提高.

圖1 特征波場偏移成像示意圖Fig.1 A sketch of characteristic wave imaging

4 特征波場走時反演

傳統(tǒng)的波動方程走時反演方法(Luo和Schuster，1991)利用互相關(guān)方法估計觀測數(shù)據(jù)與模擬數(shù)據(jù)的時差，需要對兩者進(jìn)行加窗處理，從而提取初至波或者折射波.該方法需要拾取特定的波型并在時間域加窗，繁重的拾取工作限制了該方法的應(yīng)用.

利用特征波場合成，單道地震記錄被分離為一系列單獨(dú)的震相.此時，可直接利用獨(dú)立的震相進(jìn)行反演，無需對地震記錄加窗處理.更重要的是，特征數(shù)據(jù)使得多尺度、逐層反演成為可能，如先利用折射波反演淺層速度，再利用由淺到深的反射波反演中深層速度.

由于特征波場只含有單獨(dú)的波型(震相)，可用互相關(guān)函數(shù)(Luo和Schuster，1991)測量觀測數(shù)據(jù)與模擬數(shù)據(jù)的時差，無需對地震數(shù)據(jù)加窗.本文中，我們將特征波場合成算法應(yīng)用初至波場分解，并應(yīng)用于透射波走時反演中.至于特征波場分解在反射波走時反演中的應(yīng)用，可以參考馮波(2015)，本文不再贅述.

若觀測數(shù)據(jù)d與模擬數(shù)據(jù)u經(jīng)過特征波場合成后，得到的初至震相分別可以表示為dDBF與uDBF.此時，觀測數(shù)據(jù)dDBF與模擬數(shù)據(jù)uDBF的互相關(guān)函數(shù)可以表示為

f(xr,xs;τ)=

∫dDBF(xr,t+τ;xs)uDBF(xr,t;xs)dt.

(9)

(10)

基于走時殘差

L2范數(shù)的目標(biāo)函數(shù)可以表示為

(11)

經(jīng)過推導(dǎo)，目標(biāo)泛函梯度(馮波，2015)為

(12)

其中，g0(x,-t|xr,0)代表逆時傳播的格林函數(shù)，u0為入射場，m為慢度平方，y(t)為伴隨源：

(13)

對于初至走時反演，dDBF與uDBF代表對初至波的特征波場合成，因此無需對數(shù)據(jù)進(jìn)行加窗處理.

梯度類的迭代算法為：

(14)

其中，αk為步長，P為梯度預(yù)條件算子.

5 數(shù)值試驗(yàn)

根據(jù)上文的分析，我們可以利用特征波場合成實(shí)現(xiàn)不同波型(震相)的分離，而分離后的單獨(dú)震相可以用于后續(xù)的成像及反演.

在數(shù)值實(shí)驗(yàn)中，我們首先分析不同信噪比的數(shù)據(jù)對射線參數(shù)估計以及特征波場合成的影響，然后，用一個簡單模型(凹陷模型)展示特征波場成像的效果及其優(yōu)勢.最后，用特征波場合成方法提取初至波場，并用于波動方程初至走時反演中.

5.1 特征波場合成

我們用Marmousi模型正演觀測數(shù)據(jù).速度模型為：X方向網(wǎng)格數(shù)目nx=737，Z方向網(wǎng)格數(shù)目nz=250；網(wǎng)格間隔均為12.5 m.觀測系統(tǒng)設(shè)置為：第一炮的坐標(biāo)為625 m，炮間隔25 m，總共301炮.檢波器為固定排列，第一道的坐標(biāo)為625 m，檢波器間隔12.5 m，總共601道.用10 Hz主頻的Ricker子波作為震源.

首先，選取中心地震道的炮點(diǎn)和檢波點(diǎn)坐標(biāo)為(1850, 3000).圖2a中，從左到右分別是xs0=1850 m處的共炮點(diǎn)道集，局部空間窗內(nèi)(xr0=3000 m)的地震道及其傾斜疊加得到的τ-p譜.從τ-p譜中可以拾取局部平面波的出射射線參數(shù)pr，它代表了局部平面波在地表的出射方向.我們拾取了能量最大的4個包絡(luò)對應(yīng)的出射射線參數(shù).

為了測量局部平面波的入射射線參數(shù)，需要抽取共檢波點(diǎn)道集.圖2b中，從左到右分別是xr0=3000 m處的共檢波點(diǎn)道集，局部空間窗內(nèi)(xs0=1850 m)的地震道及其傾斜疊加得到的τ-p譜.從τ-p譜中可以拾取局部平面波的出射射線參數(shù)ps，它代表了局部平面波在地表的入射方向.與圖2a中對應(yīng)，我們拾取能量最大的四組入射射線參數(shù).

利用圖2中拾取的入射和出射射線參數(shù)，根據(jù)公式(5)，可進(jìn)行特征波場合成.圖3a中，第1—4道為合成的特征波場.其中，第1道代表對直達(dá)波的特征波場合成，它將直達(dá)波從地震記錄中提取.第2—4道代表對反射波的特征波場合成，它將反射波從地震記錄中提取.第5道為重構(gòu)的中心道(第1—4道疊加)，與原始中心地震道(第6道)相比，只保留了我們提取的震相，提高了特征波的信噪比.同時，考慮到特征波場(第1—4道)的時間局部特征，可以(在變換域)進(jìn)一步壓縮存儲.

接著，為了測試噪聲對特征波場合成的影響，我們對原始數(shù)據(jù)加入了高斯白噪聲(S/N=100).圖4中展示了相應(yīng)的道集及其τ-p譜(所有參數(shù)含義同圖2).可以看出，高斯隨機(jī)噪聲(S/N=100)降低了τ-p譜的分辨率，因而可能會對τ-p譜的自動拾取有一定的影響，但對人工拾取的影響較小.圖3b中，展示了對圖4中的數(shù)據(jù)進(jìn)行特征波場合成的結(jié)果(所有參數(shù)同圖2a).顯然，特征波場合成之后，仍然能夠得到高信噪比的單獨(dú)震相.這是由于特征波場合成的過程中，同相的信號相互增強(qiáng)，而其他信號相互干涉.若進(jìn)一步增強(qiáng)噪聲強(qiáng)度(S/N=20)，此時，隨機(jī)噪聲完全壓制了有效信號，即使在τ-p譜也無法拾取射線參數(shù).若我們?nèi)匀挥谜鎸?shí)的射線參數(shù)(圖2中拾取的結(jié)果)，特征波場合成結(jié)果如圖4c所示.其中，直達(dá)波的波形仍然具有較為保真的形態(tài)，而后續(xù)的幾個反射波則受到了一定程度的噪聲干擾，略有震蕩.

上述試驗(yàn)說明，本文中的特征波場合成算法十分穩(wěn)健，隨機(jī)噪聲對其影響較小.當(dāng)噪聲能量逐漸增強(qiáng)時，τ-p譜的分辨率會逐漸降低，直至無法拾取準(zhǔn)確的射線參數(shù)，進(jìn)而導(dǎo)致特征波場合成算法失效.

5.2 特征波成像

本文用洼陷模型測試特征波場成像方法的有效性，其速度模型如圖5a所示.模擬數(shù)據(jù)共201炮，炮間距20 m，每一炮301道接收，道間距10 m.本節(jié)分別對該模型數(shù)據(jù)進(jìn)行傳統(tǒng)的Kirchhoff PSDM和特征波成像來進(jìn)行效率效果的對比.圖5b是Kirchhoff PSDM的成像剖面，由于其成像過程是通過等時面疊加收斂繞射波，成像剖面上留下了畫弧噪聲.圖5c是特征波場成像得到的成像剖面，由于其通過對特征域數(shù)據(jù)進(jìn)行“搬家”，其畫弧范圍控制在局部范圍內(nèi)，其畫弧噪聲明顯弱于傳統(tǒng)Kirchhoff PSDM.由于特征波的成像過程中可以得到射線方向，張角道集可以方便的輸出.傳統(tǒng)Kirchhoff PSDM一般只能輸出地表偏移距道集，在對后續(xù)的速度分析或?qū)游龀上竦倪m應(yīng)性方面上來說，張角道集質(zhì)量要優(yōu)于地表偏移距道集.圖5(d，e)分別為洼陷模型3500 m處的Kirchhoff PSDM地表偏移距道集和特征波成像的張角道集，對比兩圖可見，張角道集的成像質(zhì)量要高于地表偏移距道集.

圖2 (A)和(B)分別為共炮點(diǎn)/檢波點(diǎn)道集，局部窗內(nèi)的地震道以及局部平面波的入射/出射線參數(shù)譜Fig.2 (A) and (B) is the common shot/receiver gather, local seismic traces and the ray parameter spectrum, respectively

圖3 Marmousi模型正演記錄的特征波場合成對比(a) 無噪聲; (b) 加入高斯白噪聲，S/N=100; (c) 加入高斯白噪聲，S/N=20Fig.3 Comparison of characteristic wavefield decomposition

圖4 (A)和(B)分別為原始數(shù)據(jù)加入高斯隨機(jī)噪聲后的共炮點(diǎn)/檢波點(diǎn)道集，局部窗內(nèi)的地震道以及局部平面波的入射/出射線參數(shù)譜Fig.4 (A) and (B) is the common shot/receiver gather with Gaussian random noise,local seismic traces and the ray parameter spectrum, respectively

特征波成像由于只對特征數(shù)據(jù)進(jìn)行反投影(特征數(shù)據(jù)是對原地震數(shù)據(jù)的一個壓縮處理)，其成像效率要遠(yuǎn)遠(yuǎn)高于傳統(tǒng)的Kirchhoff PSDM或者傳統(tǒng)的射線束成像.本文對比了Kirchhoff PSDM和特征波成像在洼陷模型上的成像效率，具體數(shù)據(jù)見表1，數(shù)據(jù)得到了有效的壓縮，同時成像效率也有近一個數(shù)量級的提高.由于特征波場包含了數(shù)據(jù)方向信息，特征波場成像方法可以方便的推廣到目標(biāo)體成像，圖5f展示了僅僅對最后一個反射層的面向目標(biāo)成像結(jié)果.對特定目標(biāo)成像將會進(jìn)一步提高效率，縮短處理周期.

5.3 特征波場反演

在本節(jié)中，我們將特征波場合成算法分離初至波型，并用于波動方程初至走時層析中.速度模型(BP2004速度模型的右端部分，橫向約26 km，深度約12 km)如圖6a所示.觀測系統(tǒng)為：炮點(diǎn)深度12.5 m，炮點(diǎn)間隔為200 m，共107炮.檢波器深度12.5 m，間隔400 m，每炮54道.觀測系統(tǒng)總共5778道.初始速度模型采用一維模型，其中，水體速度以及水底深度已知.總共進(jìn)行18次反演迭代，反演得到的速度模型見圖6b.反演得到了整體背景速度的變化趨勢，其中模型右側(cè)的高速隆起反演得較好，但沒有反演出淺層小尺度的低速異常.圖7中為初至波走時對比曲線.其中，紅色、綠色和藍(lán)色曲線分別是真實(shí)模型、初始模型和第18次反演的速度模型中的初至走時.可以看出，真實(shí)的初至走時與反演得到的速度模型中的初至走時基本吻合，實(shí)現(xiàn)了初至走時反演的目標(biāo).

表1 K-PSDM 和 CWI 效率對比

6 討論與結(jié)論

本文提出了一種地震數(shù)據(jù)的稀疏表達(dá)方法.不同于傳統(tǒng)的變換域中的數(shù)據(jù)稀疏表達(dá)理論，該方法利用局部平面波的傳播方向(慢度矢量)，在中心炮檢點(diǎn)處同時進(jìn)行波束合成，從而將地震數(shù)據(jù)投影到局部平面波域(高維空間，二維觀測對應(yīng)5D數(shù)據(jù)；三維觀測對應(yīng)9D數(shù)據(jù))中.由于波束合成后的地震數(shù)據(jù)描述了局部平面波的方向特征，因此我們稱之為特征波場.相應(yīng)的，我們稱這種波場合成(分離)算法為特征波場合成方法.

特征波場合成的思想是利用局部平面波到達(dá)不同檢波器的時差，對地震記錄進(jìn)行時移疊加，使得給定的入射/出射慢度的局部平面波相干加強(qiáng)而其他信號相互抵消，從而分離地震記錄中的不同的波型(震相).因此，該方法對含有噪聲的地震數(shù)據(jù)有較好的適應(yīng)能力.文中的特征波場合成正問題需要估計射線參數(shù)，當(dāng)?shù)卣饠?shù)據(jù)受噪聲干擾嚴(yán)重時，難以在常規(guī)τ-p譜中自動估計局部平面波的射線參數(shù)(慢度矢量).此時，可以考慮基于反演理論的特征波場合成方法.在壓縮感知理論框架下同時反演局部平面波及其傳播方向，從而提高特征波合成的自動化程度并保持方法的穩(wěn)健性.

然而，特征波場合成正問題要求局部炮/檢排列在同一水平面上.因此，對于地表高程劇烈變化的起伏地表，該方法不再成立.此外，若炮點(diǎn)(檢波器)采樣過于稀疏，也會影響射線參數(shù)估計的精度.此時，可以僅考慮單次波束合成.

通過特征波場合成，可以將地震數(shù)據(jù)分解為單獨(dú)的震相.此時，傳統(tǒng)的Kirchhoff及高斯束偏移方法可以實(shí)現(xiàn)從“沿等時面的畫弧”到向“向反射點(diǎn)歸位搬家投影”的轉(zhuǎn)變.相應(yīng)的，常速介質(zhì)中的偏移脈沖響應(yīng)由“橢圓/球”退化為點(diǎn)響應(yīng).同時，由于在偏移過程中不需要對炮/檢射線參數(shù)掃描，因而疊前偏移效率將會大幅提高.若炮點(diǎn)的慢度矢量(入射射線參數(shù))難以估計，特征波場的偏移成像可以自動退化為高斯束類的成像方法.

更進(jìn)一步地，特征波偏移成像這種一對一映射的偏移方法非常適合與基于角度道集的速度層析反演方法融合成一體.它們具有同樣的數(shù)據(jù)空間和成像空間.這種特征波場驅(qū)動的高效偏移成像和速度層析反演融為一體的技術(shù)體系會改變目前的地震成像流程，與大規(guī)模計算機(jī)集群和現(xiàn)代圖像處理方法結(jié)合，在很大程度上可以實(shí)現(xiàn)半自動或自動的地震波成像.

圖5 速度模型(a)，Kirchhoff PSDM的成像剖面(b)，特征波場成像剖面(c), Kirchhoff PSDM的成像點(diǎn)道集(d)，特征波場成像點(diǎn)道集(e)及面向目標(biāo)成像結(jié)果(f)Fig.5 Velocity model (a), image produced by Kirchhoff PSDM (b),characteristic wavefield PSDM (c), Kirchhoff PSDM common imaging point gathers in offset domain, characteristic wavefield PSDM angle-domain CIGs, target-oriented migration result (f)

圖6 速度模型(a)及反演得到的速度模型(b)Fig.6 The true velocity model (a) and the inverted result (b)

圖7 初至波走時對比，背景為觀測的炮集((a)和(b)分別為兩個不同位置的炮集).紅色、綠色和藍(lán)色曲線分別是真實(shí)模型、初始模型和第18次反演的速度模型中的初至走時Fig.7 The comparison of first-arrival travel-times. The red, green and blue curves are the true, initial and inverted first-arrival travel-time. The background of (a) and (b) is the observed shot gather from different location

Bai L S, Liu Y K, Lu H Y, et al. 2014. Curvelet-domain joint iterative seismic data reconstruction based on compressed sensing.ChineseJ.Geophys. (in Chinese), 57(9): 2937-2945, doi: 10.6038/cjg20140919.

Billette F, Lambaré G. 1998. Velocity macro-model estimation from seismic reflection data by stereotomography.GeophysicalJournalInternational, 135(2): 671-690.

Herrmann F J, Hennenfent G. 2008. Non-parametric seismic data recovery with curvelet frames.Geophys.J.Int., 173(1): 233-248.

Feng B. 2015. Data domain wave equation tomographyfor velocity inversion (in Chinese)[Ph. D. thesis]. Shanghai: Tongji University.

Gao F, Zhang P, Wang B, Dirks V. 2006. Fast beam migration—A step toward interactive imaging. ∥76th Annual International Meeting, SEG, Expanded Abstracts, 2470-2473.

Geng Y, Wu R S, Gao J H. 2011. Efficient Gaussian packets representation and seismic imaging. ∥81st Annual International Meeting, SEG, Expanded Abstracts, 3398-3403. Hale D. 1992. Migration by the Kirchhoff, Slant Stack and Gaussian beam Methods. Center for Wave Phenomena, Colorado School of Mines, Research Report CWP-126.Hill N R. 1990. Gaussian beam migration.Geophysics, 55(11): 1416-1428.

Hu J T, Wang H Z, Wang X W. 2013. Beam forming by Radon spectrum constraint and its application.GeophysicalProspectingforPetroleum(in Chinese), 52(4): 335-338.Iturbe I, Roux P, Virieux J, et al. 2009. Travel-time sensitivity kernels versus diffraction patterns obtained through double beam-forming in shallow water.TheJournaloftheAcousticalSocietyofAmerica, 126(2): 713-720.Krüger F, Weber M, Scherbaum F, et al. 1993. Double beam analysis of anomalies in the core-mantle boundary region.Geophys.Res.Lett., 20(14): 1475-1478.

Li H, Wang H Z, Feng B, et al. 2014. Efficient pre-stack depth migration method using characteristic Gaussian packets.ChineseJ.Geophys.(in Chinese), 57(7): 2258-2268, doi: 10.6038/cjg20140720.

Liu J, Palacharla G. 2011. Multiarrival Kirchhoff beam migration.Geophysics, 76(5) WB109-WB118.

Liu S Y, Cai J X, Wang H Z, et al. 2012. Technology and application of beam ray prestack imaging method in foothill areas.GeophysicalProspectingforPetroleum(in Chinese), 51(6): 598-605.

Liu Y, Sacchi M D. 2002. De-multiple via a fast least-squares hyperbolic radon transform. ∥72nd Ann. Internat. Mtg., Soc. Expl. Geophys.. SEG Expanded Abstracts, 2182-2185.

Luo Y, Schuster G T. 1991. Wave-equation traveltime inversion.Geophysics, 56(5): 645-653.

Niu B H, Sun C Y, Zhang Z J, et al. 2001. Polynomial radon transform.ChineseJ.Geophys. (in Chinese), 44(2): 263-271.Rost S, Thomas C. 2002. Array seismology: Methods and applications.Rev.Geophys., 40(3): 2-1-2-27.Sacchi M D, Ulrych T J. 1995. High-resolution velocity gathers and offset space reconstruction.Geophysics, 60(4): 1169-1177.

Skarsoulis E K, Cornuelle B D. 2004. Travel-time sensitivity kernels in ocean acoustic tomography.J.Acoust.Soc.Am., 116(1): 227-238.Sondergaard P L. 2007. Finite discrete gabor analysis. http:∥www2.mat.dtu.dk/people/P.Soendergaard/articles/soender-phd.pdf.

Sun H, Schuster G T. 2001. 2-D wavepath migration.Geophysics, 66(5): 1528-1537.

Wang X W, Wang H Z. 2014. A research of high-resolution plane-wave decomposition based on compressed sensing.ChineseJ.Geophys. (in Chinese), 57(9): 2946-2960, doi: 10.6038/cjg20140920.

Xu S, Zhang Y, Pham D, et al. 2005. Antileakage Fourier transform for seismic data regularization.Geophysics, 70(4): V87-V95.

附中文參考文獻(xiàn)

白蘭淑, 劉伊克, 盧回憶等. 2014. 基于壓縮感知的Curvelet域聯(lián)合迭代地震數(shù)據(jù)重建. 地球物理學(xué)報, 57(9): 2937-2945, doi: 10.6038/cjg20140919.

馮波. 2015. 數(shù)據(jù)域波動方程層析速度反演方法研究 [博士論文]. 上海: 同濟(jì)大學(xué).

胡江濤, 王華忠, 王雄文. 2013. 拉東譜約束的射線束形成方法及應(yīng)用. 石油物探, 52(4): 335-338.

李輝, 王華忠, 馮波等. 2014. 特征高斯波包疊前深度偏移方法. 地球物理學(xué)報, 57(7): 2258-2268, doi: 10.6038/cjg20140720.

劉少勇, 蔡杰雄, 王華忠等. 2012. 山前帶地震數(shù)據(jù)射線(束)疊前成像方法研究與應(yīng)用. 石油物探, 51(6): 598-605.

牛濱華, 孫春巖, 張中杰等. 2001. 多項(xiàng)式Radon變換. 地球物理學(xué)報, 44(2): 263-271.王雄文, 王華忠. 2014. 基于壓縮感知的高分辨率平面波分解方法研究. 地球物理學(xué)報, 57(9): 2946-2960, doi: 10.6038/cjg20140920.

(本文編輯 汪海英)

Characteristic wavefield decomposition, imaging and inversion with prestack seismic data

WANG Hua-Zhong, FENG Bo*, LIU Shao-Yong, HU Jiang-Tao, WANG Xiong-Wen, LI Hui

WavePhenomenaandInversionImagingGroup(WPI),TongjiUniversity,Shanghai200092,China

The sparse expression of seismic data is important to signal processing, migration and parameter inversion. The conventional linear Radon transform is generalized to invert for both the local plane-wave and shot/receiver ray parameters simultaneously, which can compress the prestack seismic data in the local plane-wave domain. Consequently, it will be convenient to the subsequent ray-based imaging and inversion.Based on the combination of shot/receiver ray parameters, a new plane-wave prediction model is presented which can predict the local plane-wave within the spatial window of local shot/receiver center. The basic assumption is the local linear property of the wave-front. The seismic traces within the spatial window of local shot/receiver center are time-shifted and stacked according to the shot/receiver ray parameters and spatial distance, and the seismic phase with the given shot/receiver ray parameter can be separated naturally. Since the presented time-shift and stack procedure relies on the propagation direction of local plane-wave, we call it characteristic wave decomposition (CWD). Besides, the inversion-based CWD under the framework of compressed sensing (CS) theory can invert for both the local plane-wave and shot/receiver ray parameters simultaneously. Due to the pre-estimated ray-parameters for shot/receiver center, the loop for shot/receiver ray-parameters is avoided. Therefore the efficiency of conventional ray-based migration can be increased remarkably, and the ellipsoid migration response in constant media can be simplified to a point response. The new imaging scheme is called the characteristic wave imaging (CWI).First, the proposed CWD method is tested with the Marmousi synthetic data. Numerical examples show that the seismic phases can be separated completely. Meanwhile, even if the original seismic traces are contaminated by severe random noise, the numerical results demonstrate the robustness of CWD method.Then, the CWI method is tested with the Waxian synthetic data. Due to the high S/N ratio, the seismic data are compressed considerably. Compared with Kirchhoff PSDM method, the imaging efficiency is increased over one order of magnitude. Due to the pre-estimated ray-parameters, the angle domain CIGs can be produced directly.Finally, the CWD method is applied to the separation of first-arrival waveform, serving as the input data for wave-equation based transmission traveltime tomography.The proposed CWD method can compress the prestack seismic data in the local plane-wave domain. Under the framework of compressed sensing theory, both the local plane-wave and shot/receiver ray parameters can be inverted simultaneously. However, the local linear property of wave-front is the basic assumption of CWD. Therefore, it is not suitable for complex topography. Besides, in case of sparse shot/receiver sampling, it is difficult to invert for the shot/receiver ray-parameter simultaneously. The double beam forming procedure in CWD will be degenerated to the classical beam forming method.The presented CWI scheme can map the local plane-wave into the model space with the pre-estimated shot/receiver ray-parameters, so that the imaging efficiency is increased and the migration noise is attenuated. Meanwhile, the angle domain CIGs can be produced easily. Combined with the angle-domain tomography theory, it may be beneficial to the automation of velocity model building procedure.

Sparse expression of seismic data; Characteristic wavefield decomposition; Stereo data space; Local plane wave; Characteristic wavefield imaging; Characteristic wavefield inversion

10.6038/cjg20150617.

國家重點(diǎn)基礎(chǔ)研究計劃項(xiàng)目(2011CB201002)，國家自然科學(xué)基金(41374117)，國家重大專項(xiàng)(2011ZX05005-005-008HZ, 2011ZX05006-002, 2011ZX05023)以及中國石化地球物理重點(diǎn)實(shí)驗(yàn)室開放基金(33550006-14-FW2099-0026)資助.

王華忠，男，1964生，教授，博士生導(dǎo)師，主要從事地震波傳播、偏移成像與反演、近地表與中深層速度分析(反演)和速度模型建立、地震數(shù)據(jù)的規(guī)則化處理等方面的研究工作.E-mail:herbhuak@vip.163.com

*通訊作者 馮波，男，1984年生，博士，主要從事地震波偏移成像與反演成像等理論與應(yīng)用研究.E-mail:ancd111@163.com

10.6038/cjg

P631

2014-02-19，2015-05-21收修定稿

王華忠， 馮波， 劉少勇等.2015.疊前地震數(shù)據(jù)特征波場分解、偏移成像與層析反演.地球物理學(xué)報，58(6):2024-2034,

Wang H Z, Feng B, Liu S Y, et al. 2015. Characteristic wavefield decomposition, imaging and inversion with prestack seismic data.ChineseJ.Geophys. (in Chinese),58(6):2024-2034,doi:10.6038/cjg20150617.