黃建平, 崔 超, 劉夢(mèng)麗

(1.中國(guó)石油大學(xué)地球科學(xué)與技術(shù)學(xué)院,山東青島 266580;2.海洋國(guó)家實(shí)驗(yàn)室海洋礦產(chǎn)資源評(píng)價(jià)與探測(cè)技術(shù)功能實(shí)驗(yàn)室,山東青島 266071)

由于固有的非線性特征,全波形反演對(duì)初始模型和地震數(shù)據(jù)中的低頻成分具有較強(qiáng)的依賴性[1-2]?；陬l率由低到高的多尺度反演(multi-scale full waveform inversion, MFWI)是提高全波形反演穩(wěn)定性的有效方法[3-5],然而該方法仍然受限于地震數(shù)據(jù)中最低有效頻率的大小[6]。受采集能力和處理技術(shù)的限制,實(shí)際地震數(shù)據(jù)中的低頻成分往往不可靠。為此國(guó)內(nèi)外學(xué)者提出利用包絡(luò)[7-11]、小波變換[12]等方法恢復(fù)地震數(shù)據(jù)中的低頻信息。Alkhalifah[13]指出,波形反演問(wèn)題的關(guān)鍵在于保證反演波數(shù)的連續(xù)性。全波形反演能夠恢復(fù)的波數(shù)成分不僅與地震數(shù)據(jù)的頻率有關(guān),還與地震波傳播的散射角相關(guān)[6,14-16]?；诖?筆者在傳統(tǒng)的時(shí)間域多尺度反演方法基礎(chǔ)上,發(fā)展一種基于頻率-波數(shù)濾波的聯(lián)合多尺度反演方法(hybrid multi-scale full waveform inversion, HMFWI),利用大散射角信息彌補(bǔ)地震數(shù)據(jù)中缺失的低頻成分,以降低傳統(tǒng)全波形反演方法對(duì)低頻數(shù)據(jù)的依賴性。

1 原 理

1.1 時(shí)間域全波形反演基本原理

基于L2范數(shù)的全波形反演方法的目標(biāo)泛函[1,17]為

(1)

式中,uobs和ucal分別為由震源s激發(fā)在r接收點(diǎn)位置處t時(shí)刻的觀測(cè)數(shù)據(jù)和正演數(shù)據(jù);m為模型參數(shù)。

目標(biāo)泛函關(guān)于模型參數(shù)的梯度項(xiàng)可以通過(guò)伴隨狀態(tài)法[18-19]求得:

(2)

Sadj=ucal-uobs.

(3)

全波形反演通過(guò)對(duì)模型參數(shù)的迭代更新求得最優(yōu)解,其迭代公式可以表示為

mk+1=mk-αkHkmχ.

(4)

式中,k為迭代次數(shù);Hk為對(duì)梯度的修正項(xiàng),可采用共軛梯度法、L-BFGS方法等,本文中采用共軛梯度法;αk為步長(zhǎng),可由拋物插值法獲得[20]。

1.2 基于Wiener濾波器的多尺度反演方法

利用原始信號(hào)與目標(biāo)信號(hào)的頻譜信息,Wiener濾波器能夠?qū)⒃葱盘?hào)轉(zhuǎn)化為十分接近目標(biāo)信號(hào)的形式。Boonyasiriwate等[5]提出了基于Wiener濾波器的高效時(shí)間域波形反演方法,其中Wiener濾波器可以表示為

(5)

式中,Wtarget(f)為目標(biāo)子波的頻譜;Worigional(f)為原始子波的頻譜;f為頻率;ε為保證穩(wěn)定性的極小值。在地震數(shù)據(jù)子波已知的情況下,利用Wiener濾波器可以將原高頻地震數(shù)據(jù)轉(zhuǎn)化為主頻較低的數(shù)據(jù),轉(zhuǎn)化后的地震數(shù)據(jù)可以表示為

u′=F-1(WWiener(f)U(f)) .

(6)

式中,U(f)為原始地震記錄的頻譜;F-1為傅里葉反變換。

利用濾波器(5),并結(jié)合高效的反演頻率選擇策略[21],Boonyasiriwate等[5]實(shí)現(xiàn)了在時(shí)間域的多尺度全波形反演方法,其計(jì)算效率明顯高于傳統(tǒng)方法。然而通過(guò)對(duì)公式(5)分析可得,與一般的基于頻率濾波的多尺度反演方法相同,Wiener濾波器不能對(duì)地震數(shù)據(jù)中缺失的低頻信息進(jìn)行恢復(fù),因此其適用性在很大程度上依賴于地震數(shù)據(jù)的最低有效頻率。在地震數(shù)據(jù)質(zhì)量較差、低頻信息缺失嚴(yán)重時(shí),Wiener并不能緩解全波形反演的非線性問(wèn)題。

1.3 基于Poynting矢量的多尺度反演方法

由地震波的散射理論可知[6,14],全波形反演對(duì)地下介質(zhì)的分辨率主要與散射角θ(本文中參照Z(yǔ)hou 等[22]定義散射角)和地震波波長(zhǎng)λ相關(guān)。因此對(duì)于低波數(shù)成分(背景場(chǎng))的反演不僅僅依賴于地震數(shù)據(jù)的頻率成分,也依賴于地震波的散射角信息。前人在此基礎(chǔ)上提出了多種全波形反演的優(yōu)化算法:由于大偏移距一般對(duì)應(yīng)于大散射角信息,Shipp和Singh[23]提出利用數(shù)據(jù)加窗的方法,首先在早期反演大偏移距數(shù)據(jù)獲得模型低波數(shù)組分,然后進(jìn)行近偏移距數(shù)據(jù)反演的多尺度反演方法,然而該方法難以避免大偏移距數(shù)據(jù)的周波跳躍問(wèn)題[24];Xie[6]以及Alkhalifah[15]等提出通過(guò)散射角信息對(duì)梯度濾波以控制模型更新的波數(shù)成分,但是在地震數(shù)據(jù)低頻信息缺失或者數(shù)據(jù)采集質(zhì)量不高時(shí)該方法的適用性有限[6];Alkhalifah和Wu[16]利用散射角濾波實(shí)現(xiàn)了梯度中背景速度成分與散射成分的分離,結(jié)合多次散射,以獲得更加豐富有效的低波數(shù)成分。

如圖1所示(其中KS、KR分別表示與震源、檢波點(diǎn)相關(guān)的波數(shù)成分[28]),全波形反演模型波數(shù)與散射角和波長(zhǎng)的關(guān)系可以表示為

(7)

式中,λ為固有波長(zhǎng),由速度和頻率決定;θ為散射角,主要取決于觀測(cè)系統(tǒng)的分布;n表示波數(shù)K的方向。波長(zhǎng)的定義為

(8)

式中,v為模型速度。由式(7)和(8)可得:一方面,在觀測(cè)系統(tǒng)一定的條件下,反演所得到模型的波數(shù)成分主要取決于地震數(shù)據(jù)的頻率成分,頻率越低,反演所得到的波數(shù)成分越低;另一方面,在地震數(shù)據(jù)低頻信息不足時(shí),波形反演的低波數(shù)成分可以通過(guò)大散射角信息獲得,如回折波,大偏移距的反射波等。

圖1 全波形反演的波數(shù)成分與觀測(cè)系統(tǒng)分布之間關(guān)系Fig.1 Relationship between wavenumber recovered by FWI and acquisition geometry

基于此,本文中通過(guò)設(shè)計(jì)與散射角相關(guān)的濾波器[6,15-16],對(duì)梯度項(xiàng)按照散射角信息進(jìn)行濾波。在反演早期提取波場(chǎng)傳播中的大散射角信息,首先獲得較為準(zhǔn)確的低波數(shù)模型,然后逐漸減小散射角,以增加模型中的高波數(shù)組分,并與傳統(tǒng)的時(shí)間域多尺度反演相結(jié)合,來(lái)有效提高全波形反演在低頻成分缺失情況下的穩(wěn)定性和適應(yīng)性。

本文中利用Poynting矢量實(shí)現(xiàn)對(duì)梯度根據(jù)散射角信息的濾波。Poynting矢量通過(guò)波場(chǎng)數(shù)據(jù)以及波場(chǎng)的時(shí)間、空間導(dǎo)數(shù)獲得地震波的傳播方向,在地震數(shù)據(jù)處理中被應(yīng)用于去除逆時(shí)偏移中的低頻噪音[25-26]、進(jìn)行地震勘探的照明分析[27]、抽取共成像點(diǎn)道集[28-30]以及對(duì)全波形反演梯度的預(yù)處理[6]等方面。

在二維情況下,Poynting矢量[25]可以表示為

P(u)?-

(9)

利用Poynting矢量可以獲得震源波場(chǎng)與伴隨波場(chǎng)之間的夾角余弦值:

(10)

利用震源波場(chǎng)與伴隨波場(chǎng)的余弦值設(shè)計(jì)濾波器W(θ),將原有的梯度項(xiàng)表示為

(11)

W(θ)形式可以表示為

(12)

式中，λ1、λ2表示所期望的散射角余弦值范圍,其取值范圍為-1<λ2<λ1<1。

利用公式(11)能夠有效控制梯度中的波數(shù)成分:在反演早期采用較小的λ1、λ2值,只保留大散射角信息,以保證背景速度場(chǎng)的有效更新,然后逐漸放寬λ1、λ2的范圍,逐步對(duì)反演結(jié)果進(jìn)行細(xì)化,實(shí)現(xiàn)基于波數(shù)濾波的多尺度反演方法。必須說(shuō)明的是,濾波器W(θ)中參數(shù)λ1、λ2的確定是決定基于波數(shù)濾波的多尺度反演方法成敗的關(guān)鍵因素。對(duì)于波數(shù)濾波中閾值的確定在諸多論文中有所涉及:Alkhalifah和Wu[16]利用波數(shù)信息控制背景場(chǎng)和擾動(dòng)場(chǎng)的更新量;Xie[6]通過(guò)波數(shù)域的濾波實(shí)現(xiàn)了大尺度速度場(chǎng)的更新;Tang[31]等利用散射角信息實(shí)現(xiàn)了對(duì)全波形反演梯度中層析成分與偏移成分的加權(quán)重組以提高波形反演對(duì)中低波數(shù)組分的反演能力。然而前人并沒(méi)有給出有效的濾波閾值的控制方法,僅僅說(shuō)明波數(shù)閾值應(yīng)該逐漸變化或者通過(guò)人為控制確定。由于波形反演中的關(guān)鍵問(wèn)題是控制波數(shù)的連續(xù)性[13],而利用散射角對(duì)梯度中的波數(shù)成分進(jìn)行控制時(shí)不存在類似于基于頻率濾波的多尺度反演的高效算法[5,21]。因此本文認(rèn)為,為了能夠保證反演的準(zhǔn)確性,應(yīng)該在計(jì)算量允許的情況下盡可能地增加反演階段。除此之外,還可以根據(jù)不同階段目標(biāo)泛函的收斂情況進(jìn)行濾波閾值的確定,若前一階段目標(biāo)泛函收斂較慢,則可以適當(dāng)增大濾波范圍的變化量,反之則應(yīng)當(dāng)適當(dāng)減小變化量。

1.4基于頻率-波數(shù)濾波的聯(lián)合多尺度全波形反演流程

綜上所述,波形反演的低波數(shù)成分來(lái)源包括低頻數(shù)據(jù)和大散射角信息。常規(guī)的基于頻率濾波的多尺度反演方法對(duì)地震數(shù)據(jù)中的低頻信息具有較強(qiáng)的依賴性,在低頻數(shù)據(jù)缺失情況下往往反演失敗;而僅僅利用波數(shù)濾波在地震數(shù)據(jù)主頻較高時(shí)存在嚴(yán)重的周波跳躍問(wèn)題。因此本文中提出兩種基于頻率-波數(shù)濾波的聯(lián)合多尺度反演策略:一方面,通過(guò)首先對(duì)地震數(shù)據(jù)做低通濾波,可以有效避免單一波數(shù)濾波方法存在的周波跳躍問(wèn)題,提高反演穩(wěn)定性;另一方面,通過(guò)大散射角信息彌補(bǔ)地震低頻數(shù)據(jù)缺失所不能恢復(fù)的低波數(shù)成分,能有效降低傳統(tǒng)全波形反演方法對(duì)低頻數(shù)據(jù)的依賴性,增強(qiáng)多尺度反演方法的適應(yīng)性。

假設(shè)地震子波已知的情況下,策略一:將波數(shù)濾波與頻率濾波分離,首先進(jìn)行基于波數(shù)濾波的多尺度反演,在其提供較為準(zhǔn)確的初始模型時(shí),開(kāi)始進(jìn)行基于頻率濾波的多尺度反演。其具體流程為:首先通過(guò)式(5)所示的Wiener濾波器對(duì)原始地震數(shù)據(jù)進(jìn)行濾波,以獲得其中相對(duì)低頻的信息。然后利用公式(11)獲得不同散射角下對(duì)應(yīng)的梯度信息,散射角由大到小逐次迭代直到獲得該頻率下對(duì)應(yīng)的所有波數(shù)成分(-1≤cosθ≤1)。最后利用傳統(tǒng)的時(shí)間域多尺度反演方法對(duì)模型進(jìn)行更新,在后續(xù)的頻率成分反演中不對(duì)梯度項(xiàng)進(jìn)行波數(shù)濾波,簡(jiǎn)記為HMFWI1。

策略二:將波數(shù)濾波與頻率濾波交替進(jìn)行,將整個(gè)反演分為內(nèi)、外兩個(gè)循環(huán)。將波數(shù)濾波作為外循環(huán),內(nèi)循環(huán)做頻率由低到高的多尺度反演,然后逐漸增大外循環(huán)反演所用的散射角信息,直到獲得所有的頻率和散射角信息。相同的可以將頻率濾波作為外循環(huán),將波數(shù)濾波作為內(nèi)循環(huán),簡(jiǎn)稱為HMFWI2。

具體的反演流程如圖2所示,其中條件1表示目標(biāo)泛函收斂或者反演次數(shù)達(dá)到上限,條件2表示頻率濾波迭代完成,條件3表示波數(shù)濾波迭代完成。

本文中所采用的方法與Alkhalifah和Wu[16]所給波數(shù)濾波方法在思想上是類似的,即通過(guò)直接利用散射角信息對(duì)梯度中的波數(shù)成分進(jìn)行控制,以達(dá)到控制模型更新的波數(shù)成分的目的。然而,除了所采用的濾波方法不同之外(本文中采用Poynting矢量方法,Alkhalifah和Wu[16]采用擴(kuò)展成像方法),一方面,兩者在低波數(shù)成分的來(lái)源上是不同的:Alkhalifah和Wu[16]利用波數(shù)濾波的目的在于對(duì)梯度場(chǎng)的分解,波數(shù)濾波起到閾值作用,以獲得有效的背景場(chǎng)與散射場(chǎng),其反演所獲得的低波數(shù)成分來(lái)自于模型擾動(dòng)場(chǎng)所得到的一次散射和多次散射成分。相比而言,本文中研究的重點(diǎn)在于給出一種頻率-波數(shù)聯(lián)合濾波的多尺度反演方法,在解決單一波數(shù)濾波周波跳躍問(wèn)題的基礎(chǔ)上彌補(bǔ)由于低頻缺失不能恢復(fù)的低波數(shù)成分,通過(guò)兩種濾波方法的聯(lián)合作用直接獲得低波數(shù)組分。另一方面,由于Alkhalifah和Wu[16]的方法需要構(gòu)造較為有效的擾動(dòng)場(chǎng)獲得低波數(shù)梯度,如果在時(shí)間域進(jìn)行會(huì)極大地增加計(jì)算量,相比而言本文中方法在計(jì)算量上具有一定的優(yōu)勢(shì)。

2 模型試算

為了驗(yàn)證本文方法的適用性,采用如圖3所示的國(guó)際標(biāo)準(zhǔn)SEG/EAGE推覆體模型進(jìn)行測(cè)試。首先驗(yàn)證傳統(tǒng)多尺度全波形反演對(duì)低頻信息的依賴性,然后將本文方法在低頻數(shù)據(jù)缺失情況下進(jìn)行測(cè)試。

具體的反演參數(shù)為:模型縱橫向尺寸為1 km×4 km, 每200 m做一次正演模擬,共20炮,在海水表面從0開(kāi)始每10 m放置一個(gè)檢波器,共400個(gè)接收點(diǎn)。使用時(shí)間2階,空間8階有限差分正演模擬及PML邊界條件做正演模擬[32-34]。地震子波采用主頻為15 Hz的雷克子波,反演過(guò)程中給定震源子波和水層速度。對(duì)水層以下進(jìn)行三角平滑并乘以0.9獲得反演的初始模型(圖4),初始模型與真實(shí)模型存在較大差異。

圖2 兩種聯(lián)合多尺度全波形反演策略流程Fig.2 Workflows of two inversion strategies of hybrid multi-scale inversion

圖3 真實(shí)模型:SEG/EAGE推覆體模型Fig.3 True model:SEG/EAGE overthrust model

圖4 反演所采用的初始模型Fig.4 Initial model for inversion

2.1 時(shí)間域多尺度反演方法對(duì)低頻信息的依賴性測(cè)試

為了驗(yàn)證傳統(tǒng)時(shí)間域多尺度反演方法對(duì)低頻信息的依賴性,首先利用完整的雷克子波進(jìn)行反演測(cè)試。完整的雷克子波的時(shí)間域波形與頻率域頻譜如圖5中藍(lán)色曲線所示,為了保證反演精度,所用的反演頻率分別為5、8、11、14 Hz和全頻段(利用Boonyasiriwate等[5]給出的頻率選擇標(biāo)準(zhǔn),反演頻率應(yīng)該為5、15 Hz)。每個(gè)頻率迭代10次,采用全頻段時(shí)迭代直至誤差泛函收斂。

對(duì)原始數(shù)據(jù)做時(shí)間域多尺度的全波形反演所得的最終反演結(jié)果如圖6所示?？梢?jiàn)由于觀測(cè)系統(tǒng)偏移距較大,回折波作用較深,在回折波的作用區(qū)域反演結(jié)果與真實(shí)模型較為接近。尤其是在0.5 km以上,速度層位刻畫較為精細(xì),兩個(gè)高陡斷層反演較為準(zhǔn)確。然而在中深層,由于初始模型與真實(shí)模型相差較遠(yuǎn),并且回折波作用較弱,速度層間差異較小,反演效果較為模糊。在深度0.8 km以下范圍內(nèi),尤其是在模型兩側(cè),由于不能接收到來(lái)自于該區(qū)域的反射波信息,全波形對(duì)該區(qū)域速度更新及其微弱。總體而言,時(shí)間域多尺度反演方法較為成功地解釋了地下的構(gòu)造情況。

為了驗(yàn)證低頻信息缺失對(duì)傳統(tǒng)時(shí)間域多尺度反演的影響,采用高通濾波后的雷克子波作為震源子波,其頻譜如圖5(b)中紅色實(shí)線所示,濾去6 Hz以下的低頻信息,并在6～12 Hz之間加過(guò)渡帶以保證穩(wěn)定性,由其頻譜可見(jiàn)低頻信息損失嚴(yán)重。低頻信息的缺失會(huì)導(dǎo)致時(shí)間域波形的高頻抖動(dòng)(圖5(a)中紅色實(shí)線),由于全波形局部極小問(wèn)題主要來(lái)源于波形的復(fù)雜性,因此意味著局部極小的增多,非線性增強(qiáng),反演難度增大。

圖5 高通濾波前后雷克子波對(duì)比Fig.5 Comparison of Ricker wavelets before and after using high pass filter

圖6 全頻帶地震數(shù)據(jù)的時(shí)間域多尺度反演結(jié)果Fig.6 Inversion result of multi-scale full waveform inversion in time domain using full frequency banddata

采用與完整數(shù)據(jù)時(shí)完全相同的反演參數(shù),低頻數(shù)據(jù)缺失的時(shí)間域全波形反演結(jié)果如圖7所示。對(duì)比分析可知,由于缺少低頻信息所提供的低波數(shù)成分,傳統(tǒng)的時(shí)間域多尺度反演方法僅在橫向大于2 km的范圍內(nèi)對(duì)模型進(jìn)行了有效更新,原因在于該區(qū)域構(gòu)造較為簡(jiǎn)單并且回折波的穿透范圍較深。在橫向小于2 km范圍內(nèi),由于模型構(gòu)造較為復(fù)雜,反演結(jié)果在0.2 km以下的深層范圍內(nèi)幾乎不能提供任何有效的構(gòu)造信息,證明了在低頻缺失情況下,傳統(tǒng)的時(shí)間域多尺度全波形反演非線性問(wèn)題嚴(yán)重,不能對(duì)模型進(jìn)行有效更新,反演失效。

圖7 低頻數(shù)據(jù)缺失的時(shí)間域多尺度反演結(jié)果Fig.7 Inversion result of multi-scale full waveform inversion in time domain without low frequency data

2.2 基于頻率-波數(shù)濾波的聯(lián)合多尺度反演方法測(cè)試

在驗(yàn)證傳統(tǒng)時(shí)間域多尺度反演方法對(duì)數(shù)據(jù)中低頻信息具有較強(qiáng)依賴性的基礎(chǔ)上,為了說(shuō)明本文中方法的有效性,分別采用本文中提出的兩種策略進(jìn)行低頻數(shù)據(jù)缺失情況下的反演測(cè)試。所采用的反演參數(shù)如下:①對(duì)于HMFWI1,為了保證反演結(jié)果的準(zhǔn)確性,采用較為密集的濾波閾值采樣。在散射角多尺度反演過(guò)程中將λ2固定為-1,λ1分別取-0.8,-0.5,-0.2,0.1,0.4,0.7,1。由于多尺度反演的每個(gè)階段不一定要求目標(biāo)泛函收斂,而僅需為下一階段提供較為準(zhǔn)確的初始速度場(chǎng)[16],且在反演早期更新背景場(chǎng)情況下目標(biāo)泛函收斂較快,因此為了保證較高的計(jì)算效率,令每個(gè)階段迭代5次。其次采用傳統(tǒng)的時(shí)間域多尺度反演方法,為了保證反演精度,本文中采用等間隔的頻率采樣方法,所用的反演頻率分別為5、8、11、14 Hz和全頻段(符合Boonyasiriwate等[5]給出的頻率選擇標(biāo)準(zhǔn)),同樣為了保證計(jì)算效率,每個(gè)頻率迭代10次,采用全頻段時(shí)迭代直至誤差泛函收斂。②對(duì)于HMFWI2,將基于波數(shù)濾波的多尺度反演作為外循環(huán),所采用的λ2、λ1值與HMFWI1相同,將傳統(tǒng)時(shí)間域多尺度反演作為內(nèi)循環(huán),所采用的反演頻率與HMFWI1相同。為了保證計(jì)算效率,每個(gè)頻率段迭代5次,采用全頻段和全散射角時(shí)迭代直至目標(biāo)泛函收斂。

圖8為HMFWI1和HMFWI2的反演結(jié)果。由圖8可見(jiàn):在低頻數(shù)據(jù)缺失的情況下,相比于傳統(tǒng)的時(shí)間域多尺度反演方法(圖7),基于頻率-波數(shù)濾波的聯(lián)合多尺度全波形反演結(jié)果有明顯提升,即使在低頻信息缺失情況下也能得到較為準(zhǔn)確的反演結(jié)果。無(wú)論是橫向大于2 km的水平層位置還是較為復(fù)雜的褶皺區(qū)域,速度層位的刻畫都較為精確。橫向小于2 km處的高陡斷層位置準(zhǔn)確。

圖8 HMFWI1和HMFWI2的反演結(jié)果Fig.8 Inversion results of HMFWI1 and HMFWI2

對(duì)比HMFWI1和HMFWI2兩種方法反演結(jié)果可以發(fā)現(xiàn),HMFWI2的反演結(jié)果(圖8(b))要明顯優(yōu)于HMFWI1(圖8(a)),尤其是在推覆體下部位置(圖9)。究其原因可能是在以基于波數(shù)濾波的多尺度反演作為外循環(huán)時(shí)每個(gè)散射角循環(huán)可以利用所有地震數(shù)據(jù)的頻率組分,因此能夠恢復(fù)的波數(shù)成分也較為豐富,而HMFWI1在進(jìn)行散射角多尺度反演時(shí)只能用到有限的頻率成分(本例中為高通濾波后5 Hz雷克子波),因此能夠恢復(fù)的波數(shù)成分有限。圖10給出了所有散射角信息反演完成時(shí)HMFWI1和HMFWI2所得到的最終反演結(jié)果。相比而言HMFWI2的反演結(jié)果更為精確,模型中的主要構(gòu)造的形態(tài)及反演絕對(duì)數(shù)值都更加接近真實(shí)速度場(chǎng)。

為了更加清晰的對(duì)比幾種反演方法的差異,圖11給出了位于1 km位置處的縱向速度剖面。通過(guò)觀察可知:由于初始速度場(chǎng)與真實(shí)速度場(chǎng)相差較遠(yuǎn),傳統(tǒng)的時(shí)間域全波形反演方法僅能更新淺層(深度小于100 m)部分構(gòu)造,而本文中提出的兩種HMFWI方法能夠進(jìn)一步更新中、深部模型,極大提高了傳統(tǒng)全波形反演的適應(yīng)范圍。同時(shí),經(jīng)過(guò)分析還發(fā)現(xiàn),在模型深層反演效果方面,HMFWI2反演結(jié)果要略優(yōu)于HMFWI1。

最后,為了證明本文中給出的反演策略相對(duì)于Xie[6]提出的基于角度域波數(shù)濾波器的多尺度全波形反演方法的優(yōu)勢(shì),圖12給出了通過(guò)對(duì)梯度進(jìn)行波數(shù)濾波的多尺度全波形反演在低頻數(shù)據(jù)缺失情況下的反演結(jié)果?？梢?jiàn)反演結(jié)果沒(méi)有體現(xiàn)任何有效的速度場(chǎng)信息,說(shuō)明了在低頻信息缺失的情況下,僅僅采用波數(shù)濾波不能避免周波跳躍問(wèn)題。

圖9 速度場(chǎng)局部放大圖Fig.9 Partial enlarged details of velocity field

需要說(shuō)明的是,由于本方法要求在反演早期有較為可靠的大散射角信息,因此主要適用于大偏移距、速度場(chǎng)梯度較大的模型。在近偏移距情況下,可以利用偏移得到的偏移剖面作為先驗(yàn)信息獲得大散射角信息[16,24]。在計(jì)算效率上,本文中方法要明顯低于傳統(tǒng)的時(shí)間域全波形反演,因此需要進(jìn)一步探討高效的反演策略以提高其對(duì)復(fù)雜介質(zhì)[35]適用性。

圖10 所有散射角反演后得到的速度模型Fig.10 Velocity models obtained after inversions using all scattering angles

圖11 橫向1 km處抽得的縱向速度剖面對(duì)比Fig.11 Comparison of vertical velocity profiles at 1 km

圖12 由Xie[6]提出的“基于角度域波數(shù)濾波的多尺度全波形反演”結(jié)果Fig.12 Inverted result using “angle-domain wavenumber filter based full-waveform inversion (AWFWI)” proposed by Xie[6]

3 結(jié) 論

(1)波形反演的關(guān)鍵在于保證反演波數(shù)的連續(xù)性,當(dāng)初始模型較差時(shí),需要從地震數(shù)據(jù)中獲得足夠低的波數(shù)成分,以保證后續(xù)反演過(guò)程中觀測(cè)數(shù)據(jù)與模擬數(shù)據(jù)誤差不大于半個(gè)周期。

(2)傳統(tǒng)的時(shí)間域多尺度反演方法通過(guò)由低頻到高頻遞進(jìn)反演,在反演早期獲得較為準(zhǔn)確的低波數(shù)組分,以提高波形反演問(wèn)題的收斂性。但是該方法在本質(zhì)上依賴于地震數(shù)據(jù)中的低頻信息,在低頻信息缺失或質(zhì)量較差以及初始模型不能保證最低有效頻率下的觀測(cè)數(shù)據(jù)與模擬數(shù)據(jù)誤差小于半個(gè)周期時(shí),該方法反演失敗。

(3)根據(jù)地震波散射理論,全波形反演能夠反演的模型低波數(shù)成分不僅僅來(lái)源于低頻數(shù)據(jù),也可以來(lái)源于大散射角信息。因此當(dāng)?shù)卣饠?shù)據(jù)中低頻數(shù)據(jù)缺失時(shí),通過(guò)設(shè)計(jì)與散射角相關(guān)的濾波器,提取梯度中的大散射角成分,以彌補(bǔ)由于低頻缺失不能恢復(fù)的低波數(shù)成分,能夠保證后續(xù)反演收斂性。

(4)將基于波數(shù)濾波的多尺度全波形反演方法與傳統(tǒng)的時(shí)間域多尺度全波形反演方法相結(jié)合,一方面可以有效彌補(bǔ)低頻信息缺失時(shí)傳統(tǒng)多尺度反演方法不能恢復(fù)的低波數(shù)成分;另一方面可以保證波數(shù)濾波反演方法能夠避免周波跳躍現(xiàn)象。與傳統(tǒng)的時(shí)間域多尺度全波形反演方法相比,本文方法即使在低頻數(shù)據(jù)缺失情況下也能得到較為準(zhǔn)確的反演結(jié)果,有效地降低了波形反演對(duì)低頻信息的依賴性。

參考文獻(xiàn):

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

[2] FICHTNER A. Full seismic waveform modelling and inversion [M]. Heidelberg:Spinger-Verlag,2010.

[3] FICHTNER A. Multiscale full waveform inversion[J]. Geophysical Journal International, 2013,194(1):534-556.

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

[5] BOONYASIRIWATC, VALASEK P, ROUTHP,et al. An efficient multiscale for time-domain wave form tomography[J].Geophysics, 2009,74(6):WCC59-WCC68.

[6] XIE X B. An angle-domain wavenumber filter for multi-scale full-waveform inversion: SEG International Exhibition and Annual Meeting, New Orleans, October 18-23, 2015 [C].Tulsa: Society of Exploration Geophysicists,2015.

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

[8] CHI B, DONG L, LIU Y. Full waveform inversion method using envelope objective function without low frequency data[J].J Appl Geophys, 2014,109:36-46.

[9] LUO J, WU R S. Seismic envelope inversion: reduction of local minima and noise resistance [J]. Geophysical Prospecting, 2015,63:597-614.

[10] 敖瑞德,董良國(guó),遲本鑫.不依賴子波、基于包絡(luò)的FWI初始模型建立方法研究[J].地球物理學(xué)報(bào),2015,50(6):1998-2010.

AO Ruide, DONG Liangguo, CHI Benxin. Source-independent envelope-based FWI to build an initial model[J]. Chinese Journal of Geophysics, 2015,50(6):1998-2010.

[11] HUANG C, DONG L G, CHI B X. Elastic envelope inversion using multicomponent seismic data with filtered-out low frequency [J]. Applied Geophysics, 2015,12(3):362-377.

[12] YANHUA O Y, FREDERIK J S. Multiscale adjoint waveform-difference tomography using wavelets [J]. Geophysics, 2014,79(3):WA79-WA95.

[13] ALKHALIFAH T. Full-model wavenumber inversion: an emphasis on the appropriate wavenumber continuation[J].Geophysics, 2016,81(3):R89-R98.

[14] XIE X B, JIN S W, WU R S. Wave equation based seismic illumination analysis [J].Geophysics,2006,71(5):S169-S177.

[15] ALKHALIFAH T. Scattering-angle based filtering of the waveform inversion gradients[J] Geophysical Journal International, 2015,200:363-373.

[16] ALKHALIFAH T, WU Z D. Multi scattering inversion for low-model wave numbers[J]. Geophysics, 2016,81(6):R417-R428.

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

[18] PLESSIX RE. 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.

[19] FICHTNER A, TRAMPERT J. Hessian kernels of seismic data functionals based upon adjoint techniques[J]. Geophysical Journal International, 2011,185(2):775-798.

[20] K?HN D. Time domain 2D elastic full waveform tomography[D]. Kiel: Kiel University, 2011.

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

[22] ZHOU W, BROSSIER R, OPERTO S, et al. Full waveform inversion of diving & reflected waves for velocity model building with impedance inversion based on scale separation[J].Geophysical Journal International, 2015,202(3):1535-1554.

[23] SHIPP R M, SINGH S C. Two dimensional full wavefield inversion of wide aperture marine seismic streamer data [J]. Geophysical Journal International, 2002,151:325-344.

[24] WANG F, DANIELA D, HERVé C, et al. Wave form inversion based on wave field decomposition[J]. Geophysics, 2016,81(6):R457-R470.

[25] YOON K, MARFURT K J. Reverse-time migration using the Poynting vector[J]. Exploration Geophysics, 2006,37(1):102-107.

[26] EDVALDO S A, REYNAM C P, ADRIANO W G S. Symplectic scheme and the Poynting vector in reverse-time migration [J].Geophysics, 2014,79(5):S163-S172.

[27] YAN Z, XIE X B. Full-wave seismic illumination and resolution analyses: a Poynting-vector-based method[J]. Geophysics,2016,81(6):S447-S458.

[28] THOMAS A D, GRAHAM A W. RTM angle gathers using Poynting vectors:SEG International Exhibition and Annual Meeting, San Antonio, October 18-23, 2011[C].Tulsa: Society of Exploration Geophysicists, 2011.

[29] JIA Y, WARREN R. Improving the stability of angle gather computation using Poynting vectors:SEG International Exhibition and Annual Meeting, Texas, September 22-27, 2013[C].Tulsa: Society of Exploration Geophysicists,2013.

[30] JIA Y, THOMAS A D. Reverse time migration angle gathers using Poynting vectors [J].Geophysics, 2016,81(6):S511-S522.

[31] TANG Y X, LEE S, ANATOLY B, et al. Tomographically enhanced full wave field inversion:SEG International Exhibition and Annual Meeting, Texas, September 22-27,2013[C].Tulsa: Society of Exploration Geophysicists,2013.

[32] ALAN R L. Fourth-order finite-difference P-SV seismograms[J]. Geophysics, 1988,53(11):1425-1436.

[33] 李娜,李振春,黃建平,等.Lebedev網(wǎng)格與標(biāo)準(zhǔn)交錯(cuò)網(wǎng)格耦合機(jī)制下的復(fù)雜各向異性正演模擬[J].石油地球物理勘探,2014,49(1):121-131.

LI Na, LI Zhenchun, HUANG Jianping, et al. Numerical simulation with coupling Lebedev and standard staggered grid schemes for complex anisotropic media[J]. Oil Geophysical Prospecting, 2014,49(1):121-131.

[34] 黃建平,楊宇,李振春,等.基于 M-PML 邊界的Lebedev網(wǎng)格起伏地表正演模擬方法及穩(wěn)定性分析[J].中國(guó)石油大學(xué)學(xué)報(bào)(自然科學(xué)版),2016,40(4):47-56.

HUANG Jianping, YANG Yu, LI Zhenchun, et al. Lebedev grid finite difference modeling for irregular free surface and stability analysis based on M-PML boundary condition[J]. Journal of China University of Petroleum (Edition of Natural Science),2016,40(4):47-56.

[35] 李振春,雍鵬,黃建平,等.基于矢量波場(chǎng)分離彈性波逆時(shí)偏移成像[J].中國(guó)石油大學(xué)學(xué)報(bào)(自然科學(xué)版),2016,40(1):42-48.

LI Zhenchun, YONG Peng, HUANG Jianping, et al. Elastic wave reverse time migration based on vector wave field seperation[J]. Journal of China University of Petroleum(Edition of Natural Science), 2016,40(1):42-48.