王官超,杜啟振

(中國(guó)石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院,山東青島266580)

?

基于包絡(luò)目標(biāo)函數(shù)的彈性波波形反演

王官超,杜啟振

(中國(guó)石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院,山東青島266580)

摘要:彈性波場(chǎng)的復(fù)雜性致使反演問題非線性性增強(qiáng),尋優(yōu)過程極易陷入局部極值。多尺度反演能在一定程度上提高反演的穩(wěn)定性及分辨率,但多尺度反演策略無(wú)法解決實(shí)際數(shù)據(jù)的低頻缺失問題。針對(duì)該問題,提出了一種新的基于包絡(luò)目標(biāo)函數(shù)的彈性波多尺度反演策略。首先將包絡(luò)反演方法應(yīng)用到彈性波,恢復(fù)地下彈性介質(zhì)參數(shù)中的低波數(shù)信息；然后將包絡(luò)反演結(jié)果作為常規(guī)多尺度反演的輸入,進(jìn)一步提高反演的穩(wěn)定性。Marmousi2模型試算結(jié)果表明該策略能夠更加準(zhǔn)確地反演縱、橫波速度以及密度信息。

關(guān)鍵詞：包絡(luò)反演；低頻缺失；多尺度；彈性波全波形反演

全波形反演技術(shù)(Full Waveform Inversion,FWI)是一種從全波場(chǎng)地震記錄中定量提取地下介質(zhì)參數(shù)的數(shù)據(jù)匹配技術(shù),該技術(shù)利用波場(chǎng)的走時(shí)、振幅以及相位等多種信息,通過不同優(yōu)化方法獲得地下彈性介質(zhì)參數(shù),其分辨率明顯高于傳統(tǒng)基于射線理論的反演方法,反演結(jié)果不僅能為疊前深度偏移提供更為可靠的偏移速度,而且能直接提供精度較高的地下介質(zhì)彈性參數(shù)信息。但是,該技術(shù)存在計(jì)算量大、非線性性強(qiáng)、初始速度場(chǎng)構(gòu)建困難以及波場(chǎng)之間的不確定性等缺點(diǎn),導(dǎo)致波形反演容易陷入循環(huán)跳動(dòng)和局部極小值等問題,當(dāng)涉及到彈性多參數(shù)反演時(shí)這些問題表現(xiàn)得尤為嚴(yán)重,其中初始模型的準(zhǔn)確與否直接影響反演結(jié)果的正確性。

通常情況下,背景宏觀速度模型可以利用初至走時(shí)層析[1]、立體層析[2-3]、反射/折射層析、疊加速度以及偏移速度分析等方法獲得。當(dāng)初始模型遠(yuǎn)離真實(shí)模型時(shí),全波形反演目標(biāo)函數(shù)存在多極值,往往導(dǎo)致最優(yōu)化方法收斂到局部極小[4-5]。多尺度反演策略[6-8]能在一定程度上解決該問題。但是,由于采集技術(shù)的制約，實(shí)際地震勘探中很難記錄到5Hz以下的頻率成分,準(zhǔn)確的含低波數(shù)信息的初始速度模型的構(gòu)建仍然是FWI中亟需解決的關(guān)鍵問題。針對(duì)初始模型構(gòu)建困難這一瓶頸,許多學(xué)者進(jìn)行了大量研究,SHIN等[9-10]提出了利用Laplace域阻尼波場(chǎng)零頻分量反演長(zhǎng)波長(zhǎng)宏觀速度模型作為頻率域全波形反演的初始速度模型。事實(shí)上,將Laplace域和Fourier域的波形反演集合成一個(gè)過程,既可以提高計(jì)算效率,同時(shí)也可以提高算法穩(wěn)定性以及降低對(duì)初始速度模型的依賴性。XU[11]在反演過程中引入反偏移技術(shù)，進(jìn)行了基于反射波的波形反演研究,并在實(shí)際資料應(yīng)用中取得較好效果。胡光輝等[12]采用基于早至波的特征波波形反演建模方法解決收斂過程中由于初始模型不準(zhǔn)確導(dǎo)致的局部極小問題,使基于特征波的全波形反演逐漸受到人們的關(guān)注。地震信號(hào)經(jīng)過包絡(luò)變換后的能量主要集中在極低頻區(qū)域,這種含有極低頻信息的包絡(luò)信號(hào)可以有效反映地下介質(zhì)的大尺度信息。因此,WU等[13-14]在聲波假設(shè)下提出基于包絡(luò)目標(biāo)函數(shù)的全波形反演方法建立低頻速度模型。通過對(duì)低頻信息進(jìn)行補(bǔ)償、頻率外推等方式,減少低頻信息缺失帶來(lái)的不利影響,同樣是一種可行策略。

基于前人的研究成果,本文提出了一種基于彈性波包絡(luò)反演方法的多尺度反演策略。首先將聲波包絡(luò)反演方法推廣應(yīng)用到彈性波,恢復(fù)地下彈性介質(zhì)參數(shù)中的大尺度信息,降低彈性多參數(shù)反演的非線性性；其次將彈性波包絡(luò)反演獲得的含有低頻信息的反演結(jié)果作為常規(guī)多尺度反演的初始模型,進(jìn)一步解決反演過程中存在的不穩(wěn)定性問題。Marmousi2模型試算結(jié)果表明,結(jié)合彈性波包絡(luò)反演和多尺度反演的反演策略能夠更加準(zhǔn)確地反演縱、橫波速度以及密度信息,實(shí)現(xiàn)彈性波多參數(shù)反演。

1基于包絡(luò)目標(biāo)函數(shù)的低頻反演方法

常規(guī)全波形反演通過匹配觀測(cè)數(shù)據(jù)與模擬數(shù)據(jù),使其數(shù)據(jù)殘差δu(t)=u(t)-d(t)達(dá)到最小,定量提取地下介質(zhì)彈性參數(shù)m,目標(biāo)函數(shù)定義為:

(1)

在彈性介質(zhì)中,多分量地震數(shù)據(jù)可通過速度應(yīng)力方程進(jìn)行數(shù)值模擬合成,即:

式中:ρ(x)是密度；u(x,t)和T(x,t)分別為速度分量和應(yīng)力分量；F(x,t)為震源函數(shù)；C(x)為彈性常數(shù)。利用高階交錯(cuò)網(wǎng)格有限差分方法進(jìn)行波動(dòng)方程求解[15],采用PML吸收邊界條件。

包絡(luò)反演與常規(guī)時(shí)間域反演的區(qū)別在于其匹配的是地震數(shù)據(jù)的包絡(luò)而非地震數(shù)據(jù)。根據(jù)聲波包絡(luò)反演的目標(biāo)函數(shù)[12]可以直接給出彈性波包絡(luò)反演目標(biāo)函數(shù):

(3)

式中:d(t)=[dx(t),dz(t)]和u(t)=[ux(t),uz(t)]分別為觀測(cè)和合成的多分量地震記錄；eobs(t)和esyn(t)分別為多分量地震記錄d(t)和u(t)的包絡(luò)；dh(t)和uh(t)為相應(yīng)波形的希爾伯特變換；E為包絡(luò)殘差；參數(shù)p控制包絡(luò)信號(hào)的能量,可以是任意正數(shù),p值越大包絡(luò)信號(hào)中深層反射的能量損失越嚴(yán)重,導(dǎo)致包絡(luò)反演無(wú)法刻畫深層結(jié)構(gòu),在本文的數(shù)值測(cè)試中p=2。

通過推導(dǎo)彈性波包絡(luò)目標(biāo)函數(shù)((3)式)相對(duì)于模型參數(shù)m的一階偏導(dǎo)可得:

(4)

(4)式可表示為雅克比矩陣J以及有效殘差向量η的乘積:

(5)

(6a)

(6b)

(6c)

式中:v=(vx,vz),w=(wx,wz,wxx,wzz,wzx)分別代表正傳波場(chǎng)和反傳殘差波場(chǎng)。因?yàn)椴煌瑓?shù)之間存在著相互關(guān)系,目標(biāo)函數(shù)對(duì)于其它相關(guān)參數(shù)m=[vP,vS,ρ]的梯度可以通過鏈?zhǔn)椒▌t進(jìn)行求解。另外,迭代過程采用拋物線插值方法估計(jì)更新步長(zhǎng)。

包絡(luò)反演假設(shè)地震波形為調(diào)幅信號(hào),包絡(luò)算子是一種非線性的解調(diào)算子,可反映隱藏在調(diào)幅信號(hào)中的低頻調(diào)制信號(hào)的波形變化。圖1a和圖1b分別展示了兩個(gè)余弦衰減信號(hào)疊加的調(diào)幅信號(hào)和單道地震記錄的波形及其包絡(luò),圖1c和圖1d為對(duì)應(yīng)的頻譜特征。從頻譜上可看出,包絡(luò)變換后信號(hào)的能量主要集中在低頻區(qū)域,橫坐標(biāo)越接近0，振幅譜能量就越強(qiáng)(圖1c和圖1d),包絡(luò)反演正是利用這些極低頻信息來(lái)恢復(fù)地下介質(zhì)參數(shù)的長(zhǎng)波長(zhǎng)分量。

圖1　調(diào)幅信號(hào)及其包絡(luò)(a),單道地震記錄及其包絡(luò)(b),對(duì)應(yīng)圖1a的頻譜特征(c)以及對(duì)應(yīng)圖1b的頻譜特征(d)

2基于包絡(luò)目標(biāo)函數(shù)的多尺度反演策略

BUNKS等[6]采用Hamming窗低通濾波器對(duì)地震子波和數(shù)據(jù)進(jìn)行時(shí)間尺度分解,首次實(shí)現(xiàn)了時(shí)間域多尺度全波形反演策略。BOONYASIRIWAT等[19]指出利用Hamming窗進(jìn)行濾波會(huì)有高頻泄漏。為解決Hamming窗濾波的高頻泄漏問題,可以選用其它無(wú)泄漏或少泄漏的濾波窗函數(shù)。本文中采用維納濾波器進(jìn)行低通濾波,實(shí)現(xiàn)多尺度反演策略,維納濾波器表達(dá)式為:

(7)

式中:fWiener是維納濾波器；WOriginal是原始信號(hào)；WTarget為期望得到的信號(hào)；ω是角頻率；ε是為了避免不穩(wěn)定現(xiàn)象發(fā)生而設(shè)定的一個(gè)很小的值。利用維納濾波器將波場(chǎng)數(shù)據(jù)分解成不同頻帶,由低頻至高頻帶進(jìn)行逐頻帶反演，從而實(shí)現(xiàn)多尺度反演策略,降低反演過程陷入局部極小值的風(fēng)險(xiǎn)。圖2a為利用20Hz Ricker子波生成的單炮Z分量記錄,圖2b 為利用5Hz Ricker子波對(duì)圖2a的記錄進(jìn)行維納濾波后的地震記錄。由于實(shí)際地震數(shù)據(jù)往往缺失低頻信息,基于低通濾波的常規(guī)多尺度反演策略無(wú)法反演模型參數(shù)中的大尺度信息,導(dǎo)致反演陷入局部極值或循環(huán)跳動(dòng)。

圖2　利用20Hz Ricker子波生成的單炮記錄(a)以及利用5Hz Ricker子波對(duì)圖2a進(jìn)行維納濾波后的單炮記錄(b)

由于包絡(luò)算子具有解調(diào)地震數(shù)據(jù)中隱含低頻信息的作用，因此，利用包絡(luò)算子解調(diào)出的低頻信息可以有效恢復(fù)地下介質(zhì)的長(zhǎng)波長(zhǎng)成分,進(jìn)而解決常規(guī)多尺度反演無(wú)法克服的低頻缺失或不足問題。同時(shí),包絡(luò)變換后的信號(hào)在富含低頻信息的同時(shí)也存在相應(yīng)的高頻信息。圖3是主頻為7Hz(實(shí)線和虛線)、20Hz(點(diǎn)線和點(diǎn)劃線)Ricker子波及其包絡(luò)的頻譜響應(yīng),可以看出包絡(luò)信號(hào)攜帶的高頻信息與原始信號(hào)主頻帶的高低是相對(duì)的。如果原信號(hào)主頻較高，變換后的信號(hào)頻率也相對(duì)較高，反之亦然。

圖3　20Hz Ricker子波及其包絡(luò)的頻譜響應(yīng)以及7Hz Ricker子波及其包絡(luò)的頻譜響應(yīng)

我們通過數(shù)值算例來(lái)說(shuō)明包絡(luò)中高頻信息對(duì)包絡(luò)反演結(jié)果的影響。圖4是以線性模型作為初始模型分別利用圖3中兩種震源子波進(jìn)行彈性包絡(luò)反演獲得的縱波速度,可以看出,當(dāng)高頻信息較多時(shí)，包絡(luò)反演結(jié)果不穩(wěn)定(圖4b黑色矩形框)。為提高包絡(luò)反演的穩(wěn)定性,在進(jìn)行包絡(luò)反演之前利用類似于時(shí)間域多尺度反演對(duì)原始地震數(shù)據(jù)進(jìn)行低通濾波處理,這樣就把彈性包絡(luò)反演與常規(guī)多尺度反演有機(jī)結(jié)合在一起,自適應(yīng)的實(shí)現(xiàn)了一種從極低頻開始的多尺度反演策略。

因此,基于包絡(luò)目標(biāo)函數(shù)的多尺度反演策略的具體實(shí)現(xiàn)過程為:首先,對(duì)地震數(shù)據(jù)進(jìn)行時(shí)間尺度分解；其次,采用低通濾波后的地震數(shù)據(jù)實(shí)施包絡(luò)反演,恢復(fù)地下彈性介質(zhì)參數(shù)中的大尺度信息；最后,將彈性波包絡(luò)反演結(jié)果作為初始模型進(jìn)行常規(guī)多尺度反演以獲取精確的最終反演結(jié)果。

圖4　彈性波包絡(luò)反演得到的縱波速度a 7Hz震源子波; b 20Hz震源子波

3Marmousi2模型測(cè)試

在反演測(cè)試過程中采用兩種不同類型的震源子波合成多分量地震數(shù)據(jù)(主頻為20Hz的全頻帶Ricker子波和低頻缺失(<7Hz)Ricker子波)來(lái)驗(yàn)證彈性波包絡(luò)反演恢復(fù)地下模型參數(shù)大尺度信息的能力及其在地震數(shù)據(jù)缺失低頻情況下的適用性。為保證包絡(luò)反演的穩(wěn)定性,首先利用7Hz Ricker子波對(duì)兩種類型震源子波生成的地震記錄進(jìn)行維納低通濾波,然后進(jìn)行包絡(luò)反演。圖6和圖7 為以線性模型作為初始模型進(jìn)行包絡(luò)反演迭代30次后得到的平滑背景場(chǎng)。可以看到,針對(duì)兩種不同震源子波包絡(luò)反演獲得的背景模型構(gòu)造非常相近,這說(shuō)明了包絡(luò)反演對(duì)低頻信息的依賴性很弱,即使在低頻缺失的情況下同樣可以反演得到模型中的低頻信息,為常規(guī)多尺度波形反演提供含有低頻信息的初始速度模型,降低全波形反演的非線性性。

圖5　真實(shí)Marmousi2縱波速度模型(a)及其線性模型(b)

圖6　線性模型作為初始模型的彈性波包絡(luò)反演30次迭代后的結(jié)果(全頻帶震源子波)a 縱波速度; b 橫波速度; c 密度

圖7　線性模型作為初始模型的彈性波包絡(luò)反演30次迭代后的結(jié)果(低頻缺失震源子波)a 縱波速度; b 橫波速度; c 密度

圖8和圖9分別為以圖6和圖7中包絡(luò)反演結(jié)果作為初始模型,選取5,10,20Hz的Ricker子波對(duì)地震記錄進(jìn)行維納濾波進(jìn)而實(shí)施多尺度反演得到的最終結(jié)果。相比于直接利用線性模型作為初始模型進(jìn)行多尺度反演的結(jié)果(圖10和圖11),利用包絡(luò)反演加多尺度反演得到的最終結(jié)果在構(gòu)造(黑色橢圓處)上得到了明顯提高,構(gòu)造假象明顯減少。對(duì)比圖10和圖11可以看出,利用低頻缺失震源子波的反演結(jié)果相比于全頻帶震源子波反演結(jié)果的構(gòu)造假象更為嚴(yán)重,這說(shuō)明傳統(tǒng)全波形反演對(duì)震源的頻譜特性更為敏感。

圖8　以圖6所示包絡(luò)反演結(jié)果作為初始模型得到的多尺度反演結(jié)果(全頻帶震源子波)a 縱波速度; b 橫波速度; c 密度

圖9　以圖7所示包絡(luò)反演結(jié)果作為初始模型得到的多尺度反演結(jié)果(低頻缺失震源子波)a 縱波速度; b 橫波速度; c 密度

圖10　以線性模型作為初始模型得到的多尺度反演結(jié)果(全頻帶震源子波)a 縱波速度; b 橫波速度; c 密度

圖11　以線性模型作為初始模型得到的多尺度反演結(jié)果(低頻缺失震源子波)a 縱波速度; b 橫波速度; c 密度

4結(jié)論

通過理論研究與模型試算得到以下認(rèn)識(shí):①包絡(luò)算子可以解調(diào)出地震數(shù)據(jù)中隱含的低頻信息。包絡(luò)反演利用解調(diào)出的低頻信息能有效恢復(fù)模型中的大尺度信息,在地震數(shù)據(jù)缺失低頻的情況下同樣適用。②包絡(luò)信號(hào)在富含低頻信息的同時(shí)也存在相應(yīng)的高頻信息。高頻信息較多時(shí)會(huì)導(dǎo)致包絡(luò)反演不穩(wěn)定,通過將包絡(luò)反演與基于低通濾波的多尺度反演相結(jié)合,自適應(yīng)的實(shí)現(xiàn)了一種從極低頻開始的多尺度反演策略,最大限度地降低了反演對(duì)初始模型的依賴性。

實(shí)際應(yīng)用是全波形反演的挑戰(zhàn),野外數(shù)據(jù)低頻缺失或不足是影響波形反演獲得全波數(shù)帶參數(shù)估計(jì)結(jié)果的主要因素。包絡(luò)反演技術(shù)可恢復(fù)隱藏在地震資料中的低頻信息,建立低頻初始模型,具有很好的發(fā)展前景。

參考文獻(xiàn)

[1]NOLET G.Seismic tomography:with applications in global seismology and exploration geophysics [M].Bonston:Kluwer Academic Publishers,1987,1-24

[2]WOODWARD M J.Wave-equation tomography[J].Geophysics,1992,57(1):15-26

[3]PRIEUX V,OPERTO S,BROSSIER R,et al.Application of acoustic full waveform inversion to the synthetic Valhall model[J].Expanded Abstracts of 77thAnnual Internat SEG,Mtg,2007:2268-2272

[4]VIRIEUX J,OPERTO S.An overview of full-waveform inversion in exploration geophysics[J].Geophysics,2009,74(6):WCC1-WCC26

[5]楊勤勇,胡光輝,王立歆.全波形反演研究現(xiàn)狀及發(fā)展趨勢(shì)[J].石油物探,2014,53(1):77-83

YANG Q Y,HU G H,WANG L X.Research status and development trend of full waveform inversion [J].Geophysical Prospecting for Petroleum,2014,53 (1):77-83

[6]BUNKS C,SALECK F M,ZALESKI S,et al.Multiscale seismic waveform inversion[J].Geophysics,1995,60(5):1457-1473

[7]SIRGUE L,PRATT R G.Efficient waveform inversion and imaging:a strategy for selecting temporal frequencies[J].Geophysics,2004,69(1):231-248

[8]李媛媛,李振春,張凱.頻率域多尺度彈性波全波形反演[J].石油物探,2015,54(3):317-323

LI Y Y,LI Z C,ZHANG K.Multi-scale elastic full waveform inversion in frequency domain [J].Geophysical Prospecting for Petroleum,2015,54 (3):317-323

[9]SHIN C,HA W.A comparison between the behavior of objective functions for waveform inversion in the frequency and Laplace domains[J].Geophysics,2008,73(5):VE119-VE133

[10]SHIN C,CHA Y H.Waveform inversion in the Laplace domain [J].Geophysical Journal International,2008,173(3):922-931

[11]XU S,WANG D,CHEN F,et al.Full waveform inversion for reflected seismic data[J].Expanded Abstracts of 74thEAGE Annual Conference,2012:114-117

[12]胡光輝,王立歆,王杰,等.基于早至波的特征波波形反演建模方法[J].石油物探,2015,54 (1):71-76

HU G H,WANG L X,WANG J,et al.Characteristics waveform inversion based on early arrival waves [J].Geophysical Prospecting for Petroleum,2015,54 (1):71-76

[13]WU R S,LUO J R,WU B.Seismic envelope inversion and modulation signal model[J].Geophysics,2014,79(3):WA13-WA24

[14]LUO J R,WU R S.Seismic envelope inversion:reductionoflocal minima and noise resistance [J].Geophysical Prospecting,2015,63(3):597-614

[15]董良國(guó),馬在田,曹景忠.一階彈性波交錯(cuò)網(wǎng)格高階有限差分解法[J].地球物理學(xué)報(bào),2000,43(3):411-419

DONG L G,MA Z T,CAO J Z.Higher order difference method for first order elastic wave equation with staggered grid [J].Chinese Journal of Geophysics,2000,43(3):411-419

[16]TARANTOLA A.Inversion of seismic reflection data in the acoustic approximation[J].Geophysics,1984,49(8):1259-1266

[17]PLESSIX R E.A review of the adjoint-state method for computing the gradient of a functional with geophysical applications [J].Geophysical Journal International,2006,167(2):495-503

[18]WANG J,ZHOU H,TIAN Y,et al.A new scheme for elastic full waveform inversion based on velocity-stress wave equations in time domain[J].Expanded Abstracts of 82ndAnnual Internat SEG,Mtg,2012:1-5

[19]BOONYASIRIWAT C,VALASEK P,ROUTH P,et al.An efficient multiscale method for time-domain waveform tomography[J].Geophysics,2009,74(6):WCC59-WCC68

(編輯：朱文杰)

Elastic full waveform inversion based on envelope objective function

WANG Guanchao,DU Qizhen

(SchoolofGeosciences,ChinaUniversityofPetroleum,Qingdao266580,China)

Abstract:The complexity of elastic wavefield will increase the nonlinearity of inversion and make the inverse problem fall into local minima.Multi-scale inversion algorithms can improve the stability and resolution of inversion to a certain extent,but it cannot solve the problem of low frequency missing in real data.To alleviate this problem,we propose a new multi-scale scheme based on envelope inversion to obtain the elastic parameters accurately.The main procedures are as follows:Firstly,we extend the envelope inversion method to elastic medium to estimate the low wavenumber components of elastic parameters; then,the envelope inversion results are used as the initial model of conventional multi-scale inversion to further improve the stability of the inversion.The Marmousi2 model test results show that the new multi-scale scheme can estimate P- and S- wave velocities more accurately,as well as the density information.

Keywords:envelope inversion,low frequency missing,multi-scale,elastic full waveform inversion

文章編號(hào)：1000-1441(2016)01-0133-09

DOI:10.3969/j.issn.1000-1441.2016.01.017

中圖分類號(hào)：P631

文獻(xiàn)標(biāo)識(shí)碼：A

基金項(xiàng)目：國(guó)家自然科學(xué)基金(41174100和41074087)和中國(guó)石化地球物理重點(diǎn)實(shí)驗(yàn)室開放基金聯(lián)合資助。

作者簡(jiǎn)介：王官超(1990—),男,碩士在讀,主要從事全波形反演研究工作。通訊作者：杜啟振(1969—),男,教授,博士生導(dǎo)師,主要從事地震彈性波傳播理論與成像等研究工作。

收稿日期：2015-07-07;改回日期：2015-08-24。

王官超,杜啟振.基于包絡(luò)目標(biāo)函數(shù)的彈性波波形反演[J].石油物探,2016,55(1):-141

WANG Guanchao,DU Qizhen.Elastic full waveform inversion based on envelope objective function[J].Geophysical Prospecting for Petroleum,2016,55(1):-141

This research is financially supported by the National Science Foundation of China (Grant Nos.41174100,41074087) and the SINOPEC Key Laboratory of Geophysics Project.