蔡志成,顧漢明,曹靜杰

(1.河北地質(zhì)大學(xué)勘查技術(shù)與工程學(xué)院,河北石家莊050031;2.中國(guó)地質(zhì)大學(xué)地球物理與空間信息學(xué)院地球內(nèi)部多尺度成像湖北省重點(diǎn)實(shí)驗(yàn)室,湖北武漢430074)

逆時(shí)偏移(RTM)技術(shù)相對(duì)其它偏移方法具有成像精度高的優(yōu)勢(shì)[1-3],能夠處理速度橫向變化劇烈、陡傾角構(gòu)造等復(fù)雜地質(zhì)情況[4]。近年來(lái),隨著計(jì)算機(jī)技術(shù)的快速發(fā)展,RTM方法逐漸走向?qū)嵱没痆5-7]?；趶椥越橘|(zhì)假設(shè)的逆時(shí)偏移方法在理論上考慮了實(shí)際介質(zhì)中波傳播多波多分量問(wèn)題,能得到較為準(zhǔn)確的成像結(jié)果,因而備受關(guān)注。

彈性波逆時(shí)偏移問(wèn)題的研究主要包括震源、檢波點(diǎn)波場(chǎng)重建和成像條件應(yīng)用。不少學(xué)者在彈性波RTM過(guò)程中先進(jìn)行縱橫波波場(chǎng)分離然后再成像[8-10]。采用的主要方法是基于Helmholtz分解的波場(chǎng)分離[11],可分為直接波場(chǎng)分離法[12-14]和間接波場(chǎng)分離法。直接波場(chǎng)分離法是在空間域?qū)Σ▓?chǎng)進(jìn)行散度和旋度運(yùn)算分離得到縱、橫波波場(chǎng),但存在振幅相位校正問(wèn)題。AKI等[15]選擇空間域矢量分解方法解決了成像時(shí)極性反轉(zhuǎn)問(wèn)題,但其波場(chǎng)的物理意義發(fā)生了改變;ZHU[16]嘗試在波數(shù)域進(jìn)行波場(chǎng)分離,但該方法需要進(jìn)行數(shù)據(jù)域轉(zhuǎn)換以及波數(shù)域歸一化處理等才可保證其物理意義。間接波場(chǎng)分離技術(shù)則是改造波動(dòng)方程,在波場(chǎng)延拓過(guò)程中分解得到縱橫波分量。馬德堂等[17]利用各向同性介質(zhì)縱橫波分離方程實(shí)現(xiàn)波場(chǎng)完全分解。XIAO等[18]通過(guò)縱波方程得到縱波波場(chǎng)再與彈性波波場(chǎng)作差得到分離的橫波波場(chǎng),此類方法能夠得到與原波場(chǎng)一致的縱橫波波場(chǎng)。WANG等[19]將基于縱橫波矢量解耦的波動(dòng)方程應(yīng)用于彈性波逆時(shí)偏移成像,得到獨(dú)立的縱橫波偏移結(jié)果。吳瀟等[20]對(duì)比了5種波場(chǎng)分離方法在彈性波逆時(shí)偏移中的應(yīng)用效果。與直接波場(chǎng)分離法相比,間接波場(chǎng)分離法成像質(zhì)量更好,但這類方法需要改造波動(dòng)方程,增加了計(jì)算量和內(nèi)存存儲(chǔ)。

應(yīng)用準(zhǔn)確成像條件是彈性波逆時(shí)偏移成像的關(guān)鍵環(huán)節(jié)。激發(fā)時(shí)間成像條件是利用提取震源與檢波點(diǎn)波場(chǎng)對(duì)應(yīng)成像時(shí)間波場(chǎng)值,得到最終偏移結(jié)果[21-22];還有通過(guò)提取震源、檢波點(diǎn)波場(chǎng)作比值進(jìn)行成像[23-24]。該類方法計(jì)算效率高,且有一定的保幅特性,但轉(zhuǎn)換波成像物理意義有待研究。較為常用的互相關(guān)成像條件能夠利用全波場(chǎng)信息,有學(xué)者通過(guò)選擇性成像[25]和濾波的方法壓制背向散射噪聲干擾[26],但該壓制技術(shù)同時(shí)也會(huì)對(duì)有效信號(hào)造成損傷。ROCHA等[26-27]利用彈性波波動(dòng)方程及能量守恒原則,提出新的成像條件,成像結(jié)果表示反射系數(shù)能量,無(wú)需對(duì)整個(gè)波場(chǎng)進(jìn)行縱橫波分離,也避免了極性反轉(zhuǎn)等傳統(tǒng)彈性波偏移成像的問(wèn)題。張曉語(yǔ)等[28-29]在ROCHA等提出方法的基礎(chǔ)上發(fā)展了基于能量密度自解耦成像條件,此成像條件能夠有效壓制低頻散射成像干擾,主要應(yīng)用了二階方程,對(duì)于具有較高效率的交錯(cuò)網(wǎng)格一階方程尚無(wú)具體實(shí)現(xiàn)方式。本文發(fā)展了一階速度應(yīng)力彈性波能量范數(shù)成像條件的逆時(shí)偏移方法,通過(guò)構(gòu)造中間應(yīng)變量給出了一階方程情況下各分量轉(zhuǎn)換規(guī)則,提出了具體波場(chǎng)更新方式。最后通過(guò)數(shù)值算例測(cè)試了該方法對(duì)彈性波逆時(shí)偏移中低頻噪聲的壓制效果。

1 方法原理

1.1 彈性波逆時(shí)偏移理論

逆時(shí)偏移理論框架可按其實(shí)現(xiàn)步驟來(lái)描述:①震源波場(chǎng)正向外推;②檢波點(diǎn)波場(chǎng)逆時(shí)延拓;③成像條件應(yīng)用。震源波場(chǎng)正向外推可簡(jiǎn)化為如下邊值問(wèn)題[1]:

(1)

式中:Us(x,z,t)為震源波場(chǎng)值;f(t)為加載震源函數(shù)。波場(chǎng)值更新采用數(shù)值離散化的波動(dòng)方程求取。檢波點(diǎn)波場(chǎng)逆時(shí)延拓則與正向外推相反[1],可表達(dá)為:

(2)

式中:g(t)為地震記錄值;Ur(x,z,t)為檢波點(diǎn)波場(chǎng)值;Tmax為最大記錄時(shí)間。公式(2)中將檢波點(diǎn)地震記錄作為初始條件,利用波場(chǎng)延拓算子求取進(jìn)行逆時(shí)延拓。逆時(shí)偏移是記錄每一計(jì)算遞推時(shí)刻的震源波場(chǎng)和檢波點(diǎn)波場(chǎng),然后選取相應(yīng)成像條件得到最終偏移剖面。

本文方法基于完全彈性假設(shè),波場(chǎng)延拓過(guò)程采用各向同性介質(zhì)一階速度應(yīng)力彈性波方程[3]:

(3)

式中:σxx、σzz、σxz為應(yīng)力分量;vx、vz為速度分量;vP、vS分別為介質(zhì)中縱、橫波速度;ρ為密度。本文采用偽譜法對(duì)此方程進(jìn)行數(shù)值離散。

1.2 能量范數(shù)成像條件

逆時(shí)偏移中常用的互相關(guān)成像條件[3]表達(dá)為:

(4)

式中:I為成像結(jié)果;Us為震源波場(chǎng);Ur為檢波點(diǎn)波場(chǎng);e表示炮索引號(hào);t為時(shí)間。在應(yīng)用該成像條件時(shí),存在檢波點(diǎn)波場(chǎng)與震源波場(chǎng)相同路徑傳播的互相關(guān)噪聲。

能量范數(shù)成像條件首先從能量角度以一維方程聲波定義一個(gè)空間域(Ω)內(nèi)能量表達(dá)[26]:

(5)

式中:U為波場(chǎng)值;v為介質(zhì)中波速。定義波場(chǎng)能量范數(shù)為:

(6)

類似兩個(gè)波場(chǎng)U、V的內(nèi)積能量范數(shù)可表達(dá)為:

(7)

借助(4)式,U、V為正反傳波場(chǎng)時(shí),能量范數(shù)成像條件[26]為:

(8)

在彈性介質(zhì)假設(shè)條件下,公式(8)中兩個(gè)波場(chǎng)對(duì)應(yīng)為彈性波場(chǎng)。

1.3 能量范數(shù)成像條件噪聲壓制分析

下面進(jìn)一步分析上述成像條件噪聲壓制原理。將檢波點(diǎn)波場(chǎng)與震源波場(chǎng)看作是多方向的平面波合成,分別將震源和檢波點(diǎn)波場(chǎng)變換至頻率域,互相關(guān)成像條件可以表示為:

(9)

(10)

式中:ks、kr分別為震源波場(chǎng)和檢波點(diǎn)波場(chǎng)的波數(shù)向量;θ為波場(chǎng)波數(shù)向量夾角。利用頻散關(guān)系:

(11)

得到能量范數(shù)成像條件[26]為:

(12)

在ω2項(xiàng)前加入一個(gè)cos2θc項(xiàng),將其近似表示為:

(13)

當(dāng)選擇θc與θ相等時(shí),則該夾角的能量范數(shù)為0,將其反變換回空間域可以得到:

(14)

因此,可以選擇θc來(lái)壓制震源波場(chǎng)與檢波點(diǎn)波場(chǎng)任意角度互相關(guān)造成的干擾噪聲?；ハ嚓P(guān)成像條件的后散射干擾可以近似看作兩個(gè)波場(chǎng)成180°互相關(guān)干擾,這里設(shè)置θc為90°入射,即得到最終壓制反向散射干擾成像條件為:

(15)

將其應(yīng)用于二維彈性波逆時(shí)偏移成像中得到垂直分量和水平分量成像條件:

(16)

(17)

式中:(IE)x、(IE)z分別為能量范數(shù)成像條件偏移結(jié)果的水平分量和垂直分量;Usx、Usz分別為震源波場(chǎng)x分量和z分量;Urx、Urz分別為檢波點(diǎn)波場(chǎng)x分量和z分量。

1.4 能量范數(shù)成像條件彈性波逆時(shí)偏移

公式(16)和公式(17)利用一階速度應(yīng)力彈性波方程求取時(shí)需引入中間變量,首先假設(shè):

(18)

對(duì)(18)式取時(shí)間一階導(dǎo)數(shù),變換得到:

(19)

由(19)式假設(shè)的中間分量可以通過(guò)速度分量空間偏導(dǎo)求取,對(duì)(19)式利用時(shí)間二階差分和空間偽譜法進(jìn)行數(shù)值離散,有:

(20)

1.5 偽譜法數(shù)值離散求解方式

圖1 能量反射成像條件一階彈性波逆時(shí)偏移各分量更新示意

2 模型試算分析

為驗(yàn)證算法正確性,設(shè)計(jì)兩層介質(zhì)模型進(jìn)行試算分析,模型參數(shù)如圖2所示,密度由一般經(jīng)驗(yàn)公式計(jì)算得到。對(duì)單炮正演得到地震記錄進(jìn)行偏移,震源采用30Hz雷克子波,位置為(1000m,0);201道接收,道間距10m,檢波點(diǎn)深度為0,0.5ms采樣,1.0s記錄長(zhǎng)度。

圖2 兩層模型速度參數(shù)a 縱波速度; b 橫波速度

圖3為利用偽譜法彈性波互相關(guān)成像條件和能量范數(shù)成像條件對(duì)單炮記錄進(jìn)行偏移的結(jié)果。從圖3a 和圖3b可以看出,由于直達(dá)波近似水平方向傳播與檢波點(diǎn)波場(chǎng)繞射產(chǎn)生的P波和S波產(chǎn)生的干擾被較好壓制(箭頭處),在界面上方的大部分繞射產(chǎn)生的干擾很大程度上得到了消除;從圖3c和圖3d可以看出,橢圓部分低頻噪聲部分被壓制。多炮偏移疊加結(jié)果如圖4所示。從圖4a和圖4b可以看出,界面上方的低頻散射干擾明顯降低;圖4c和圖4d驗(yàn)證了能量范數(shù)成像條件的有效去噪效果。

圖3 偽譜法彈性波互相關(guān)成像條件和能量范數(shù)成像條件單炮記錄偏移結(jié)果a vx分量互相關(guān)成像條件偏移結(jié)果; b vz分量互相關(guān)成像條件偏移結(jié)果; c 能量范數(shù)成像條件x分量偏移結(jié)果; d 能量范數(shù)成像條件z分量偏移結(jié)果

圖4 偽譜法彈性波互相關(guān)成像條件和能量范數(shù)成像條件多炮偏移疊加結(jié)果a vx分量互相關(guān)成像條件偏移結(jié)果; b vz分量互相關(guān)成像條件偏移結(jié)果; c 能量范數(shù)成像條件x分量偏移結(jié)果; d 能量范數(shù)成像條件z分量偏移結(jié)果

設(shè)計(jì)凹陷模型(圖5)進(jìn)行多分量逆時(shí)偏移試算,分別測(cè)試互相關(guān)成像條件的偏移、縱橫波解耦的彈性波偏移、能量范數(shù)成像條件的彈性波逆時(shí)偏移和多炮偏移疊加結(jié)果,如圖6所示。與直接利用vx、vz分量進(jìn)行互相關(guān)成像(圖6a,圖6b)相比,基于縱橫波分離的成像PP波與SS波偏移結(jié)果更為清晰,第二界面上方的低頻干擾明顯降低(圖6c,圖6f)。PS波與SP波的偏移效果則較差(圖6d,圖6e),說(shuō)明不同類型波場(chǎng)的直接互相關(guān)成像方法帶來(lái)了更多噪聲,甚至有效信息被覆蓋,而關(guān)于轉(zhuǎn)換波的成像方法則還需要利用其它技術(shù)來(lái)完成?；谀芰糠稊?shù)成像條件的偏移結(jié)果x分量與z分量存在的低頻噪聲明顯降低(圖6g,圖6h),正反傳波場(chǎng)的非界面處互相關(guān)噪聲得到了較好的壓制。

圖5 凹陷模型速度參數(shù)a 縱波速度模型; b 橫波速度模型

圖6 凹陷模型彈性波多炮偏移疊加結(jié)果a vx分量互相關(guān)成像條件偏移; b vz分量互相關(guān)成像條件偏移疊加; c 縱橫波波場(chǎng)分離PP波偏移多炮疊加; d 縱橫波波場(chǎng)分離PS波偏移多炮疊加; e 縱橫波波場(chǎng)分離SP波偏移多炮疊加; f 縱橫波波場(chǎng)分離SS波偏移多炮疊加; g 能量范數(shù)成像條件x分量偏移多炮疊加; h 能量范數(shù)成像條件z分量偏移多炮疊加

3 復(fù)雜模型試算

截取Marmousi模型的一部分進(jìn)行本文方法的噪聲壓制測(cè)試,模型大小為2000m×700m,網(wǎng)格大小5m×5m,速度參數(shù)如圖7所示,密度由一般經(jīng)驗(yàn)公式計(jì)算得到。正演模擬共41炮,炮點(diǎn)位置從0到2000m,炮間距50m。單炮401道接收,道間距5m,分布于整個(gè)模型地表附近,0.5ms采樣,2s記錄長(zhǎng)度。炮點(diǎn)、檢波點(diǎn)深度均為0,采用30Hz雷克子波點(diǎn)震源(同時(shí)激發(fā)縱波和橫波)。從最終偏移疊加剖面(圖8)可以看出,與互相關(guān)成像條件偏移結(jié)果(圖8a,圖8b)相比,基于能量范數(shù)成像條件的結(jié)果(圖8g,圖8h)中低頻噪聲得到了較好的壓制;與縱橫波分離后的成像結(jié)果(圖8c,圖8d,圖8e,圖8f)相比,圖8g和圖8h的信噪比較高。

圖7 復(fù)雜模型速度參數(shù)a 縱波速度模型; b 橫波速度模型

圖8 復(fù)雜模型彈性波多炮偏移疊加結(jié)果a vx分量互相關(guān)成像條件偏移多炮疊加; b vz分量互相關(guān)成像條件偏移多炮疊加; c 縱橫波波場(chǎng)分離PP波偏移多炮疊加; d 縱橫波波場(chǎng)分離PS波偏移多炮疊加; e 縱橫波波場(chǎng)分離SP波偏移多炮疊加; f 縱橫波波場(chǎng)分離PS波偏移多炮疊加; g 能量范數(shù)成像條件x分量偏移多炮疊加; h 能量范數(shù)成像條件z分量偏移多炮疊加

4 結(jié)論

本文從能量范數(shù)的角度詳細(xì)分析了能量范數(shù)成像條件及其噪聲壓制原理,將這一成像條件發(fā)展到了一階速度應(yīng)力彈性波逆時(shí)偏移方法中。通過(guò)引入中間應(yīng)變分量進(jìn)行所需分量更新,實(shí)現(xiàn)了一階速度應(yīng)力彈性波方程逆時(shí)偏移方法能量范數(shù)成像條件的應(yīng)用。通過(guò)復(fù)雜模型試算,對(duì)比vx分量、vz分量互相關(guān)成像條件偏移結(jié)果,本文方法能夠壓制低頻噪聲;對(duì)比縱橫波波場(chǎng)分離偏移,本文方法無(wú)需考慮轉(zhuǎn)換波極性反轉(zhuǎn)問(wèn)題,避免進(jìn)行解耦波動(dòng)方程延拓。