包乾宗,陳俊霓,吳 浩
(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)造性地提出了利用反傳算法來(lái)計(jì)算目標(biāo)函數(shù)的梯度,避免了Jacobian矩陣的求取,從而實(shí)現(xiàn)了時(shí)間域全波形反演。GAUTHIER等[3]和MORA[4]實(shí)現(xiàn)了實(shí)際二維地震資料的全波形反演,從實(shí)踐上證明了全波形反演具有精細(xì)刻畫(huà)地下構(gòu)造及巖性特征的能力,是一種高精度的速度建模手段;同時(shí),他們指出了全波形反演理論的固有缺陷,如計(jì)算量巨大、解易陷入局部最小值、低頻信息缺失等,限制了全波形反演技術(shù)在實(shí)際資料中的應(yīng)用。
近年來(lái),隨著計(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)非線性問(wèn)題所伴隨的解易陷入局部極小值的問(wèn)題仍然是全波形反演理論進(jìn)一步發(fā)展的瓶頸。而這種強(qiáng)非線性問(wèn)題卻對(duì)初始模型有著較高的要求,若初始模型與真實(shí)模型偏差過(guò)大,計(jì)算波場(chǎng)與觀測(cè)波場(chǎng)相位差大于半個(gè)周期,會(huì)產(chǎn)生“周波跳躍”現(xiàn)象,導(dǎo)致反演陷入局部極小值而無(wú)法得到準(zhǔn)確的結(jié)果。如何為全波形反演提供一個(gè)良好的初始速度模型以緩解反演過(guò)程中的“周波跳躍”問(wèn)題是個(gè)難題。
前人發(fā)展了很多速度建模方法,如偏移速度分析、疊加速度分析以及速度層析成像等方法均可以建立較為準(zhǔn)確的初始速度模型,并能在一定程度上避免全波形反演過(guò)程中的局部極小和“周波跳躍”問(wèn)題[7]。此外,VAN LEEUWEN等[8]提出的波場(chǎng)重構(gòu)反演(wavefield reconstruction inversion,WRI)是近年來(lái)較為新穎的反演理論。他們?cè)赪RI方法的目標(biāo)函數(shù)中引入了波動(dòng)方程,使目標(biāo)解的空間擴(kuò)大到數(shù)據(jù)空間和模型空間,然后通過(guò)重構(gòu)真實(shí)波場(chǎng)來(lái)計(jì)算模型梯度,故而不需要計(jì)算反傳的殘差波場(chǎng),提高了計(jì)算效率,同時(shí)解更為穩(wěn)定,不易陷入局部極小值。段超然等[9]利用波場(chǎng)重構(gòu)反演在優(yōu)化過(guò)程中擁有較大自由度的優(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ì)初始模型的依賴(lài),可以緩解計(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)重依賴(lài)于偏移距和拉普拉斯阻尼系數(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ò)來(lái)重構(gòu)地震記錄中缺失的低頻數(shù)據(jù),該方法是利用三次樣條插值來(lái)提取地震包絡(luò),解決了希爾伯特包絡(luò)在非對(duì)稱(chēng)信號(hào)情況下出現(xiàn)的振蕩現(xiàn)象。但是,當(dāng)數(shù)據(jù)主頻過(guò)高時(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)從低頻到高頻的分頻多尺度全波形反演。此反演方法降低了反演過(guò)程陷入局部極小值的可能,抑制了“周波跳躍”現(xiàn)象對(duì)反演結(jié)果的不良影響。數(shù)值試驗(yàn)結(jié)果證實(shí),該方法可以更好地反演出地下復(fù)雜結(jié)構(gòu)的速度模型。
包絡(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ù)的建立及其梯度的求解。
在常規(guī)全波形反演中,目標(biāo)函數(shù)(J)通常根據(jù)計(jì)算波場(chǎng)與觀測(cè)波場(chǎng)之間殘差能量最小來(lái)定義:
(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ì)算中,一般使用反傳算法來(lái)求解目標(biāo)函數(shù)的梯度[2]。其簡(jiǎn)略步驟為:正向傳播源波場(chǎng),反向傳播殘差波場(chǎng)ri,j,最后計(jì)算二者的零延時(shí)互相關(guān)得到梯度。
(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ò)的頻譜
不同類(lèi)型的包絡(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)稱(chēng)為對(duì)數(shù)包絡(luò)[18]:
(6)
從圖4可以看出,地震數(shù)據(jù)的對(duì)數(shù)包絡(luò)同樣能夠反映出波形的宏觀變化。圖5與圖6分別給出了不同類(lèi)型的包絡(luò)及其頻譜。對(duì)比對(duì)數(shù)包絡(luò)和指數(shù)包絡(luò)可以看出,對(duì)數(shù)包絡(luò)在反映波形宏觀變化的同時(shí),沒(méi)有像包絡(luò)二次方那樣給予早期數(shù)據(jù)過(guò)大的權(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ù)及其不同類(lèi)型包絡(luò)(包絡(luò)二次方已做標(biāo)準(zhǔn)化處理)
圖6 單道地震數(shù)據(jù)與不同類(lèi)型包絡(luò)的頻譜(包絡(luò)二次方頻譜已做標(biāo)準(zhǔn)化處理)
此外,取正演模擬的單炮記錄(采樣間隔0.5ms,采樣長(zhǎng)度1.5s),從更大尺度上來(lái)分析地震波場(chǎng)及其包絡(luò)所包含的信息。由Marmousi模型(部分)正演得到的單炮地震記錄及其包絡(luò)(圖7)可以看出,單炮記錄完整地記錄了地下結(jié)構(gòu)的復(fù)雜變化,并且分辨率較高;但包絡(luò)只含有極少的細(xì)節(jié)信息,其主要信息反映的是地下大尺度的背景構(gòu)造。
圖7 單炮記錄(a)及其包絡(luò)(b)
與常規(guī)全波形反演類(lèi)似,包絡(luò)反演的目標(biāo)函數(shù)Jenv是根據(jù)計(jì)算波場(chǎng)包絡(luò)與觀測(cè)波場(chǎng)包絡(luò)之間殘差能量最小來(lái)定義:
(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ī)全波形反演類(lèi)似,包絡(luò)反演目標(biāo)函數(shù)的梯度也是通過(guò)反傳算法進(jìn)行求解。
雖然地震數(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ù)反演過(guò)程中,引入巴特沃斯低通濾波器,進(jìn)行分頻多尺度的全波形反演。
選取Marmousi模型的一部分(圖9a)作為真實(shí)模型(圖9b)進(jìn)行反演數(shù)值計(jì)算,采用不含任何構(gòu)造信息的一維線性模型作為初始模型(圖9c),模型的主要參數(shù)為:網(wǎng)格點(diǎn)數(shù)320×120,空間網(wǎng)格大小5m。在模型周?chē)€設(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)
在進(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ò)提取的低頻成分不足,反演過(guò)程受“周波跳躍”現(xiàn)象干擾;低頻包絡(luò)可以提取數(shù)據(jù)的極低頻成分,反演結(jié)果較好地反映了速度模型的長(zhǎng)波長(zhǎng)結(jié)構(gòu),所展示的背景構(gòu)造更為平滑。
在6個(gè)尺度(0~5,0~10,0~15,0~20,0~25,0~30Hz)下進(jìn)行分頻多尺度反演,將由低頻數(shù)據(jù)反演得到的結(jié)果作為高頻數(shù)據(jù)反演的初始模型,逐級(jí)反演,提高反演過(guò)程穩(wěn)定性,最終得到準(zhǔn)確的反演結(jié)果。按照設(shè)計(jì)的數(shù)值實(shí)驗(yàn)參數(shù),分別進(jìn)行了常規(guī)全波形反演、普通包絡(luò)+多尺度反演和低頻包絡(luò)+多尺度反演,結(jié)果如圖11所示。
常規(guī)全波形反演結(jié)果(圖11a)只能刻畫(huà)出真實(shí)模型的大致形態(tài)和輪廓,由于缺乏準(zhǔn)確的初始模型且數(shù)據(jù)主頻較高,因此反演結(jié)果受“周波跳躍”現(xiàn)象的影響較大,反演結(jié)果精度不高,未能刻畫(huà)出精細(xì)的構(gòu)造。普通包絡(luò)+多尺度反演結(jié)果(圖11b)顯示,當(dāng)初始模型含有一定的大尺度構(gòu)造(圖10a),并在后續(xù)反演過(guò)程中引入分頻多尺度策略時(shí),高頻數(shù)據(jù)的反演過(guò)程變得穩(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)確刻畫(huà)。
圖10 普通包絡(luò)(a)、低頻包絡(luò)(b)構(gòu)造的初始速度模型
圖11 常規(guī)全波形反演(a)、普通包絡(luò)+多尺度反演(b)和低頻包絡(luò)+多尺度反演(c)的結(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é)的刻畫(huà)能力更強(qiáng);對(duì)于高主頻地震數(shù)據(jù),低頻包絡(luò)的反演結(jié)果比普通包絡(luò)的反演結(jié)果更為平滑,更能反映出速度模型的長(zhǎng)波長(zhǎng)分量,作為后續(xù)反演的初始模型有著更好的結(jié)果;包絡(luò)反演與多尺度反演方法的有機(jī)結(jié)合,極大地避免了“周波跳躍”現(xiàn)象對(duì)反演過(guò)程的影響,最終提高了反演精度;與前兩種方法相比,低頻包絡(luò)對(duì)覆蓋次數(shù)較低區(qū)域的反演精度提升十分明顯,減少了對(duì)大偏移距數(shù)據(jù)的依賴(lài)。
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ù)主頻過(guò)高、反演不收斂、最終結(jié)果精度過(guò)低的問(wèn)題。
3) 本文采用包絡(luò)反演與多尺度反演有機(jī)結(jié)合的反演策略,主要針對(duì)低頻成分占比不高且數(shù)據(jù)主頻較高的地震數(shù)據(jù)難以用常規(guī)全波形反演成功反演的情況,該反演方法可操作性高,能夠明顯提升反演的精度。