基于地震數(shù)據(jù)包絡(luò)的多尺度全波形反演方法

包乾宗,陳俊霓,吳 浩

(1.長(zhǎng)安大學(xué)地質(zhì)工程與測(cè)繪學(xué)院,陜西西安710054;2.四川仁沐高速公路有限責(zé)任公司,四川成都610041)

確定地層速度模型是地震勘探的核心目標(biāo)之一。全波形反演(full waveform inversion,FWI)是在充分利用疊前地震資料中的動(dòng)力學(xué)和運(yùn)動(dòng)學(xué)信息的基礎(chǔ)上反演地層速度模型的主要方法之一[1],是目前地球物理勘探領(lǐng)域的研究熱點(diǎn)。在該領(lǐng)域,TARANTOLA[2]創(chuàng)造性地提出了利用反傳算法來計(jì)算目標(biāo)函數(shù)的梯度,避免了Jacobian矩陣的求取,從而實(shí)現(xiàn)了時(shí)間域全波形反演。GAUTHIER等[3]和MORA[4]實(shí)現(xiàn)了實(shí)際二維地震資料的全波形反演,從實(shí)踐上證明了全波形反演具有精細(xì)刻畫地下構(gòu)造及巖性特征的能力,是一種高精度的速度建模手段;同時(shí),他們指出了全波形反演理論的固有缺陷,如計(jì)算量巨大、解易陷入局部最小值、低頻信息缺失等,限制了全波形反演技術(shù)在實(shí)際資料中的應(yīng)用。

近年來,隨著計(jì)算機(jī)科學(xué)與技術(shù)和優(yōu)化算法的發(fā)展,全波形反演理論和方法的研究有了新的發(fā)展,尤其是GPU并行編碼方式的應(yīng)用使得全波形反演的計(jì)算效率大大提升。馮海新等[5]提出了基于GPU/CPU和震源隨機(jī)編碼技術(shù)的混合域全波形反演,該方法提高了反演精度,減少了炮間串?dāng)_噪聲,在GPU的加速下,計(jì)算效率顯著提升。魏哲楓等[6]使用隨機(jī)邊界技術(shù)反傳、重建震源波場(chǎng),充分發(fā)揮了GPU的計(jì)算能力,提高了全波形反演在三維情況下的計(jì)算效率和穩(wěn)健性。即便如此,基于局部尋優(yōu)算法求解強(qiáng)非線性問題所伴隨的解易陷入局部極小值的問題仍然是全波形反演理論進(jìn)一步發(fā)展的瓶頸。而這種強(qiáng)非線性問題卻對(duì)初始模型有著較高的要求,若初始模型與真實(shí)模型偏差過大,計(jì)算波場(chǎng)與觀測(cè)波場(chǎng)相位差大于半個(gè)周期,會(huì)產(chǎn)生“周波跳躍”現(xiàn)象,導(dǎo)致反演陷入局部極小值而無法得到準(zhǔn)確的結(jié)果。如何為全波形反演提供一個(gè)良好的初始速度模型以緩解反演過程中的“周波跳躍”問題是個(gè)難題。

前人發(fā)展了很多速度建模方法,如偏移速度分析、疊加速度分析以及速度層析成像等方法均可以建立較為準(zhǔn)確的初始速度模型,并能在一定程度上避免全波形反演過程中的局部極小和“周波跳躍”問題[7]。此外,VAN LEEUWEN等[8]提出的波場(chǎng)重構(gòu)反演(wavefield reconstruction inversion,WRI)是近年來較為新穎的反演理論。他們?cè)赪RI方法的目標(biāo)函數(shù)中引入了波動(dòng)方程,使目標(biāo)解的空間擴(kuò)大到數(shù)據(jù)空間和模型空間,然后通過重構(gòu)真實(shí)波場(chǎng)來計(jì)算模型梯度,故而不需要計(jì)算反傳的殘差波場(chǎng),提高了計(jì)算效率,同時(shí)解更為穩(wěn)定,不易陷入局部極小值。段超然等[9]利用波場(chǎng)重構(gòu)反演在優(yōu)化過程中擁有較大自由度的優(yōu)勢(shì),模擬出低頻部分,并以波場(chǎng)重構(gòu)反演的結(jié)果作為初始模型進(jìn)行全波形反演,得到了較好的結(jié)果。李輝等[10]提出了基于高斯束的速度層析方法,將波動(dòng)層析與射線層析的優(yōu)點(diǎn)綜合利用,在Rytov近似和Born近似的基礎(chǔ)上闡述了高斯束層析理論,并應(yīng)用于角度域成像道集偏移速度分析,得到了理想的層析結(jié)果。

在地震數(shù)據(jù)中,不同頻率數(shù)據(jù)對(duì)應(yīng)的周期不同,相比于低頻數(shù)據(jù),高頻數(shù)據(jù)頻率更高、周期更短,在進(jìn)行計(jì)算波場(chǎng)與觀測(cè)波場(chǎng)匹配時(shí)更容易出現(xiàn)“周波跳躍”現(xiàn)象。BUNKS[11]將多尺度反演方法應(yīng)用于時(shí)間域全波形反演,引入濾波技術(shù)對(duì)低頻數(shù)據(jù)進(jìn)行反演以得到大尺度背景結(jié)構(gòu),再用所得結(jié)果作為新的初始模型,對(duì)高頻數(shù)據(jù)進(jìn)行反演得到更為精細(xì)的模型結(jié)構(gòu)。PRATT[12]提出了一種頻率域全波形反演方法,該方法直接在頻率域?qū)崿F(xiàn),因而,很容易實(shí)現(xiàn)從低頻數(shù)據(jù)到高頻數(shù)據(jù)的多尺度反演。多尺度反演方法可以有效降低FWI對(duì)初始模型的依賴,可以緩解計(jì)算波場(chǎng)與觀測(cè)波場(chǎng)間的“周波跳躍”現(xiàn)象,最終獲得較好的反演結(jié)果。但是,當(dāng)數(shù)據(jù)的低頻成分不足時(shí),利用低通濾波技術(shù)難以有效反演出速度模型的長(zhǎng)波長(zhǎng)結(jié)構(gòu),從而導(dǎo)致后續(xù)反演結(jié)果不佳。研究人員發(fā)現(xiàn)利用數(shù)學(xué)變換提取低頻信息,再與全波形反演結(jié)合的方法能取得更好的結(jié)果。SHIN等[13]提出的拉普拉斯域全波形反演方法利用阻尼波場(chǎng)中的零頻率分量成功恢復(fù)出速度模型的長(zhǎng)波長(zhǎng)結(jié)構(gòu),但該方法的反演深度嚴(yán)重依賴于偏移距和拉普拉斯阻尼系數(shù)的選擇。CHI[14]認(rèn)為,雖然地震數(shù)據(jù)自身的低頻信息不足,但在地震數(shù)據(jù)的包絡(luò)中卻含有豐富的低頻成分。在此基礎(chǔ)上,構(gòu)造地震數(shù)據(jù)的包絡(luò)目標(biāo)函數(shù)并進(jìn)行反演,用以恢復(fù)地下速度模型的長(zhǎng)波長(zhǎng)結(jié)構(gòu),再以此作為常規(guī)全波形反演的初始模型,能有效降低反演的非線性程度,并顯著地提高反演結(jié)果的精度。胡勇等[15]提出了利用解調(diào)包絡(luò)來重構(gòu)地震記錄中缺失的低頻數(shù)據(jù),該方法是利用三次樣條插值來提取地震包絡(luò),解決了希爾伯特包絡(luò)在非對(duì)稱信號(hào)情況下出現(xiàn)的振蕩現(xiàn)象。但是,當(dāng)數(shù)據(jù)主頻過高時(shí),單純地采用包絡(luò)反演與常規(guī)全波形反演相結(jié)合的方式,仍會(huì)受“周波跳躍”的干擾,得不到高精度的反演結(jié)果。

受前人工作的啟發(fā),本文先對(duì)高頻數(shù)據(jù)作低通濾波處理,采用低頻數(shù)據(jù)的對(duì)數(shù)包絡(luò)進(jìn)行反演,重構(gòu)出速度模型的長(zhǎng)波長(zhǎng)結(jié)構(gòu),為后續(xù)全波形反演提供一個(gè)較為準(zhǔn)確的初始模型;再采用巴特沃斯低通濾波器實(shí)現(xiàn)從低頻到高頻的分頻多尺度全波形反演。此反演方法降低了反演過程陷入局部極小值的可能,抑制了“周波跳躍”現(xiàn)象對(duì)反演結(jié)果的不良影響。數(shù)值試驗(yàn)結(jié)果證實(shí),該方法可以更好地反演出地下復(fù)雜結(jié)構(gòu)的速度模型。

1 包絡(luò)反演理論

包絡(luò)反演目標(biāo)函數(shù)的定義和梯度的求解與常規(guī)全波形反演所采取的策略極為相似,兩者都是基于最小二乘準(zhǔn)則建立目標(biāo)函數(shù),并用反傳算法計(jì)算梯度。在本節(jié)中首先簡(jiǎn)單回顧了常規(guī)全波形反演理論的目標(biāo)函數(shù)及其梯度求解方法,然后介紹如何利用希爾伯特變換提取地震數(shù)據(jù)的包絡(luò),最后推導(dǎo)出包絡(luò)反演中目標(biāo)函數(shù)的建立及其梯度的求解。

1.1 常規(guī)全波形反演的目標(biāo)函數(shù)建立及其梯度求解

在常規(guī)全波形反演中,目標(biāo)函數(shù)(J)通常根據(jù)計(jì)算波場(chǎng)與觀測(cè)波場(chǎng)之間殘差能量最小來定義:

(1)

式中:ns,nr分別表示激發(fā)點(diǎn)數(shù)和接收點(diǎn)數(shù);ui,j與di,j分別表示第i炮、第j道的計(jì)算波場(chǎng)和觀測(cè)波場(chǎng)[16]。

對(duì)目標(biāo)函數(shù)關(guān)于模型參數(shù)m求導(dǎo)得到其梯度:

(2)

式中:ri,j=ui,j-di,j,即殘差波場(chǎng)[16]。在實(shí)際計(jì)算中,一般使用反傳算法來求解目標(biāo)函數(shù)的梯度[2]。其簡(jiǎn)略步驟為:正向傳播源波場(chǎng),反向傳播殘差波場(chǎng)ri,j,最后計(jì)算二者的零延時(shí)互相關(guān)得到梯度。

1.2 地震數(shù)據(jù)的包絡(luò)

(3)

H{f(t)}表示實(shí)信號(hào)f(t)的希爾伯特變換,希爾伯特變換定義為:

(4)

式中:E(t)為瞬時(shí)振幅,即為信號(hào)f(t)的包絡(luò);Φ(t)為瞬時(shí)相位[17]。圖1顯示了單道地震數(shù)據(jù)及其包絡(luò),可以看出,數(shù)據(jù)的包絡(luò)反映了波形在時(shí)間域上的宏觀變化。圖2顯示了單道地震數(shù)據(jù)及其包絡(luò)的頻譜,可以看出,數(shù)據(jù)包絡(luò)內(nèi)含有豐富的低頻信息。

圖1 地震數(shù)據(jù)及其包絡(luò)

圖2 地震數(shù)據(jù)及其包絡(luò)的頻譜

不同類型的包絡(luò)反映了不同的波形宏觀變化,所提取的低頻成分也不同。設(shè)E(t)p表示數(shù)據(jù)包絡(luò)的p次方,理論上講,p可以取任意正整數(shù),但在實(shí)際反演中,p值一般取1或2,分別表示數(shù)據(jù)包絡(luò)的一次方和二次方。包絡(luò)的指數(shù)形式(圖3)表明:較大的p值(p=2)會(huì)給予早期到達(dá)的數(shù)據(jù)(如直達(dá)波)更大的權(quán)重從而壓制后期到達(dá)的數(shù)據(jù),導(dǎo)致深層的反演效果不佳;較小的p值(p=1)則會(huì)壓制低頻而突出后期的反射波,難以重構(gòu)出速度模型的長(zhǎng)波長(zhǎng)結(jié)構(gòu)。設(shè)logaE(t)表示數(shù)據(jù)包絡(luò)的對(duì)數(shù),取自然數(shù)e為底數(shù)(a=e),得自然對(duì)數(shù)包絡(luò),簡(jiǎn)稱為對(duì)數(shù)包絡(luò)[18]:

(6)

從圖4可以看出,地震數(shù)據(jù)的對(duì)數(shù)包絡(luò)同樣能夠反映出波形的宏觀變化。圖5與圖6分別給出了不同類型的包絡(luò)及其頻譜。對(duì)比對(duì)數(shù)包絡(luò)和指數(shù)包絡(luò)可以看出,對(duì)數(shù)包絡(luò)在反映波形宏觀變化的同時(shí),沒有像包絡(luò)二次方那樣給予早期數(shù)據(jù)過大的權(quán)重,壓制后期數(shù)據(jù);與包絡(luò)一次方相比,對(duì)數(shù)包絡(luò)更為溫和而不尖銳,有效提取了數(shù)據(jù)的低頻成分。因此,本文研究采用對(duì)數(shù)包絡(luò)。

圖3 單道地震數(shù)據(jù)及其包絡(luò)的指數(shù)形式(包絡(luò)二次方已做標(biāo)準(zhǔn)化處理)

圖4 單道地震數(shù)據(jù)及其對(duì)數(shù)包絡(luò)

圖5 單道地震數(shù)據(jù)及其不同類型包絡(luò)(包絡(luò)二次方已做標(biāo)準(zhǔn)化處理)

圖6 單道地震數(shù)據(jù)與不同類型包絡(luò)的頻譜(包絡(luò)二次方頻譜已做標(biāo)準(zhǔn)化處理)

此外,取正演模擬的單炮記錄(采樣間隔0.5ms,采樣長(zhǎng)度1.5s),從更大尺度上來分析地震波場(chǎng)及其包絡(luò)所包含的信息。由Marmousi模型(部分)正演得到的單炮地震記錄及其包絡(luò)(圖7)可以看出,單炮記錄完整地記錄了地下結(jié)構(gòu)的復(fù)雜變化,并且分辨率較高;但包絡(luò)只含有極少的細(xì)節(jié)信息,其主要信息反映的是地下大尺度的背景構(gòu)造。

圖7 單炮記錄(a)及其包絡(luò)(b)

1.3 包絡(luò)反演目標(biāo)函數(shù)的建立及其梯度的求解

與常規(guī)全波形反演類似,包絡(luò)反演的目標(biāo)函數(shù)Jenv是根據(jù)計(jì)算波場(chǎng)包絡(luò)與觀測(cè)波場(chǎng)包絡(luò)之間殘差能量最小來定義:

(7)

對(duì)方程(7)關(guān)于模型參數(shù)m求導(dǎo),可以得到目標(biāo)函數(shù)的梯度:

(8)

式中:renvi,j=lnE(ui,j)-lnE(di,j)表示計(jì)算波場(chǎng)與觀測(cè)波場(chǎng)兩者數(shù)據(jù)包絡(luò)間的殘差。與常規(guī)全波形反演類似,包絡(luò)反演目標(biāo)函數(shù)的梯度也是通過反傳算法進(jìn)行求解。

1.4 基于地震數(shù)據(jù)包絡(luò)的多尺度全波形反演方法

雖然地震數(shù)據(jù)自身的低頻信息不足,但包絡(luò)中卻含有豐富的低頻成分,地震數(shù)據(jù)的包絡(luò)所提取出的低頻信息與地震數(shù)據(jù)的主頻可以相互對(duì)應(yīng)。圖8為不同頻率下地震數(shù)據(jù)及其包絡(luò)的頻譜,可以看出,當(dāng)?shù)卣饠?shù)據(jù)的主頻較高時(shí),其包絡(luò)的頻率也較高。因此,較高頻率的數(shù)據(jù)如果不作任何處理,仍會(huì)受“周波跳躍”現(xiàn)象的干擾而影響反演的穩(wěn)定性。

圖8 不同頻率下地震數(shù)據(jù)及其包絡(luò)的頻譜

本文借鑒BUNKS[11]的分頻多尺度反演思想,對(duì)高主頻的地震數(shù)據(jù)首先做低通濾波處理,再進(jìn)行包絡(luò)反演(低頻包絡(luò)),然后,以包絡(luò)反演的結(jié)果作為全波形反演的初始模型,并在后續(xù)反演過程中,引入巴特沃斯低通濾波器,進(jìn)行分頻多尺度的全波形反演。

2 數(shù)值計(jì)算

選取Marmousi模型的一部分(圖9a)作為真實(shí)模型(圖9b)進(jìn)行反演數(shù)值計(jì)算,采用不含任何構(gòu)造信息的一維線性模型作為初始模型(圖9c),模型的主要參數(shù)為:網(wǎng)格點(diǎn)數(shù)320×120,空間網(wǎng)格大小5m。在模型周圍還設(shè)置了10個(gè)網(wǎng)格點(diǎn)寬度的卷積完全匹配層(convolutional perfectly matched layer,CPML)作為吸收層[19]。采用地表激發(fā)地表接收的觀測(cè)系統(tǒng),在地表均勻分布了19個(gè)炮點(diǎn)和140個(gè)接收點(diǎn),時(shí)間采樣間隔為0.5ms,采樣時(shí)間為1.5s,震源采用主頻為30Hz的Ricker子波。

圖9 原始Marmousi模型(a)、本文所用的真實(shí)模型(b)和線性初始模型(c)

2.1 包絡(luò)反演構(gòu)造初始模型

在進(jìn)行包絡(luò)反演構(gòu)造初始模型的數(shù)值試驗(yàn)里,采用兩種方式進(jìn)行對(duì)比,一種是對(duì)數(shù)據(jù)不作低通濾波處理,直接進(jìn)行包絡(luò)反演;另一種是先對(duì)數(shù)據(jù)作低通濾波處理,濾去其高頻成分再進(jìn)行包絡(luò)反演(低頻包絡(luò))。圖10為普通包絡(luò)、低頻包絡(luò)構(gòu)造的初始速度模型,可以看出:由普通包絡(luò)反演得到的速度模型能夠在一定程度上重構(gòu)速度模型的大尺度構(gòu)造,但是效果較差,其根本原因在于在數(shù)據(jù)主頻較高的情況下,普通包絡(luò)提取的低頻成分不足,反演過程受“周波跳躍”現(xiàn)象干擾;低頻包絡(luò)可以提取數(shù)據(jù)的極低頻成分,反演結(jié)果較好地反映了速度模型的長(zhǎng)波長(zhǎng)結(jié)構(gòu),所展示的背景構(gòu)造更為平滑。

2.2 反演結(jié)果對(duì)比

在6個(gè)尺度(0～5,0～10,0～15,0～20,0～25,0～30Hz)下進(jìn)行分頻多尺度反演,將由低頻數(shù)據(jù)反演得到的結(jié)果作為高頻數(shù)據(jù)反演的初始模型,逐級(jí)反演,提高反演過程穩(wěn)定性,最終得到準(zhǔn)確的反演結(jié)果。按照設(shè)計(jì)的數(shù)值實(shí)驗(yàn)參數(shù),分別進(jìn)行了常規(guī)全波形反演、普通包絡(luò)+多尺度反演和低頻包絡(luò)+多尺度反演,結(jié)果如圖11所示。

常規(guī)全波形反演結(jié)果(圖11a)只能刻畫出真實(shí)模型的大致形態(tài)和輪廓,由于缺乏準(zhǔn)確的初始模型且數(shù)據(jù)主頻較高,因此反演結(jié)果受“周波跳躍”現(xiàn)象的影響較大,反演結(jié)果精度不高,未能刻畫出精細(xì)的構(gòu)造。普通包絡(luò)+多尺度反演結(jié)果(圖11b)顯示,當(dāng)初始模型含有一定的大尺度構(gòu)造(圖10a),并在后續(xù)反演過程中引入分頻多尺度策略時(shí),高頻數(shù)據(jù)的反演過程變得穩(wěn)定,反演精度得到大幅度提升;在模型中部區(qū)域(400～1200m),復(fù)雜結(jié)構(gòu)基本得以重構(gòu),地層界面清晰。低頻包絡(luò)+多尺度反演結(jié)果(圖11c)顯示,一個(gè)良好的初始模型(圖10b)對(duì)反演精度的提升效果極其顯著,此反演方法所重構(gòu)出的速度模型,整體構(gòu)造形態(tài)與真實(shí)模型基本吻合,模型的橫/縱向分辨率高,各地質(zhì)層分界明顯,淺/深部區(qū)域的細(xì)節(jié)構(gòu)造得以準(zhǔn)確刻畫。

圖10 普通包絡(luò)(a)、低頻包絡(luò)(b)構(gòu)造的初始速度模型

圖11 常規(guī)全波形反演(a)、普通包絡(luò)+多尺度反演(b)和低頻包絡(luò)+多尺度反演(c)的結(jié)果

2.3 反演結(jié)果量化分析

在真實(shí)模型(圖9b)中抽取地面距離為x=200m,x=900m和x=1400m處的速度對(duì)圖11所示的3個(gè)反演結(jié)果與真實(shí)模型進(jìn)行速度曲線擬合分析,結(jié)果如圖12所示。

從圖12可以看出:低頻包絡(luò)反演結(jié)果的速度曲線擬合程度最高,其次為普通包絡(luò)反演,常規(guī)全波形反演結(jié)果的速度擬合程度最低;在x=200m處,由于靠近模型邊界,覆蓋次數(shù)較低,低頻成分難以獲得,具有高非線性性,常規(guī)全波形反演與普通包絡(luò)反演得到的速度曲線只能得到速度隨深度的大體變化且與真實(shí)模型相差甚遠(yuǎn),而低頻包絡(luò)反演所得到速度曲線可以準(zhǔn)確反映出速度的宏觀趨勢(shì)和微小變化,并且高度擬合真實(shí)模型的速度變化曲線。

表1給出了x=200m,x=900m和x=1400m處的反演結(jié)果與真實(shí)模型的速度擾動(dòng)比。從表1可以看出,相比另兩種反演方法,低頻包絡(luò)+多尺度反演結(jié)果的速度擾動(dòng)比大幅度降低,凸顯了該方法的優(yōu)越性。

圖12 地面距離為x=200m,x=900m,x=1400m處不同反演結(jié)果與真實(shí)模型的速度擬合曲線

方法x=200mx=900mx=1400m常規(guī)FWI0.1754747630.0960891160.064598744普通包絡(luò)+多尺度0.1685456530.0524345150.044816805低頻包絡(luò)+多尺度0.0345962900.0204759800.022509714

綜合圖11、圖12和表1可以明顯看出3種反演方法的優(yōu)劣差異:相比常規(guī)全波形反演,利用包絡(luò)提取地震數(shù)據(jù)的低頻成分作初步反演,恢復(fù)出大尺度結(jié)構(gòu),為后續(xù)反演提供初始模型,最終反演結(jié)果的精度得以提升且對(duì)細(xì)節(jié)的刻畫能力更強(qiáng);對(duì)于高主頻地震數(shù)據(jù),低頻包絡(luò)的反演結(jié)果比普通包絡(luò)的反演結(jié)果更為平滑,更能反映出速度模型的長(zhǎng)波長(zhǎng)分量,作為后續(xù)反演的初始模型有著更好的結(jié)果;包絡(luò)反演與多尺度反演方法的有機(jī)結(jié)合,極大地避免了“周波跳躍”現(xiàn)象對(duì)反演過程的影響,最終提高了反演精度;與前兩種方法相比,低頻包絡(luò)對(duì)覆蓋次數(shù)較低區(qū)域的反演精度提升十分明顯,減少了對(duì)大偏移距數(shù)據(jù)的依賴。

3 結(jié)論

1) 利用數(shù)據(jù)的對(duì)數(shù)包絡(luò)提取數(shù)據(jù)的低頻信息可以重構(gòu)出模型的大尺度背景,而低頻包絡(luò)反演的結(jié)果則能反映出速度模型的大尺度輪廓和背景構(gòu)造,可作為全波形反演的初始模型。

2) 頻率越高的數(shù)據(jù)在全波形反演中越容易受“周波跳躍”的影響。對(duì)于高頻地震數(shù)據(jù),分頻多尺度反演以低頻數(shù)據(jù)反演結(jié)果作為高頻數(shù)據(jù)的初始模型,逐級(jí)反演,可以有效地解決數(shù)據(jù)主頻過高、反演不收斂、最終結(jié)果精度過低的問題。

3) 本文采用包絡(luò)反演與多尺度反演有機(jī)結(jié)合的反演策略,主要針對(duì)低頻成分占比不高且數(shù)據(jù)主頻較高的地震數(shù)據(jù)難以用常規(guī)全波形反演成功反演的情況,該反演方法可操作性高,能夠明顯提升反演的精度。