基于頻率-波數(shù)濾波的聯(lián)合多尺度全波形反演方法及其實現(xiàn)策略

黃建平, 崔 超, 劉夢麗

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

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

1 原 理

1.1 時間域全波形反演基本原理

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

(1)

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

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

(2)

Sadj=ucal-uobs.

(3)

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

mk+1=mk-αkHkmχ.

(4)

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

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

利用原始信號與目標(biāo)信號的頻譜信息,Wiener濾波器能夠?qū)⒃葱盘栟D(zhuǎn)化為十分接近目標(biāo)信號的形式。Boonyasiriwate等[5]提出了基于Wiener濾波器的高效時間域波形反演方法,其中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]實現(xiàn)了在時間域的多尺度全波形反演方法,其計算效率明顯高于傳統(tǒng)方法。然而通過對公式(5)分析可得,與一般的基于頻率濾波的多尺度反演方法相同,Wiener濾波器不能對地震數(shù)據(jù)中缺失的低頻信息進(jìn)行恢復(fù),因此其適用性在很大程度上依賴于地震數(shù)據(jù)的最低有效頻率。在地震數(shù)據(jù)質(zhì)量較差、低頻信息缺失嚴(yán)重時,Wiener并不能緩解全波形反演的非線性問題。

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

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

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

(7)

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

(8)

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

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

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

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

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

P(u)?-

(9)

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

(10)

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

(11)

W(θ)形式可以表示為

(12)

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

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

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

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

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

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

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

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

2 模型試算

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

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

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

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

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

2.1 時間域多尺度反演方法對低頻信息的依賴性測試

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖9 速度場局部放大圖Fig.9 Partial enlarged details of velocity field

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

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

圖11 橫向1 km處抽得的縱向速度剖面對比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ù)據(jù)中獲得足夠低的波數(shù)成分,以保證后續(xù)反演過程中觀測數(shù)據(jù)與模擬數(shù)據(jù)誤差不大于半個周期。

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

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

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

參考文獻(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] 敖瑞德,董良國,遲本鑫.不依賴子波、基于包絡(luò)的FWI初始模型建立方法研究[J].地球物理學(xué)報,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)交錯網(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].中國石油大學(xué)學(xué)報(自然科學(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] 李振春,雍鵬,黃建平,等.基于矢量波場分離彈性波逆時偏移成像[J].中國石油大學(xué)學(xué)報(自然科學(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.