張曉語,杜啟振,符力耘

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

全波形反演(full waveform inversion,FWI)方法是一種基于波動方程正演模擬的反演方法,該方法可在無需任何假設(shè)條件的情況下進(jìn)行波動方程全波場數(shù)值模擬,因此理論上可以提供最高分辨率的速度成像。隨著計算機(jī)技術(shù)的發(fā)展,全波形反演近年來逐漸成為勘探地球物理領(lǐng)域的研究熱點。全波形反演通過建立實際觀測記錄與正演模擬記錄的最小二乘目標(biāo)函數(shù),利用全波場的信息獲取全波數(shù)的速度模型,反演結(jié)果理論上具有最高的分辨率[1-3]。但全波形反演的精度嚴(yán)重依賴于低頻信息、長偏移距地震數(shù)據(jù)以及準(zhǔn)確的初始速度模型。

將速度模型分為淺層速度模型和中、深層速度模型兩部分。淺層模型存在的波主要包括直達(dá)波、潛水波等透射波,該類型的波由于觀測孔徑大,故可以很好地恢復(fù)淺層的低波數(shù)速度[4]。多尺度反演所采用的從低頻逐步反演到高頻的策略可以解決淺層背景速度的恢復(fù)問題,同時多尺度反演還依賴于地震數(shù)據(jù)中的低頻信息,特別是1.5～2.0Hz的地震有效信號在恢復(fù)背景速度中起了非常關(guān)鍵的作用,但在實際情況中,傳統(tǒng)的地震數(shù)據(jù)最低頻約為5Hz,低于2Hz的低頻信息幾乎不存在[5]。包絡(luò)反演將地震信號解調(diào)重構(gòu)得到地震數(shù)據(jù)包絡(luò)用于更新大尺度的地震速度[6-7]。將地震數(shù)據(jù)進(jìn)行包絡(luò)變換后,數(shù)據(jù)的頻帶移至超低頻區(qū)域,該低頻區(qū)域獨立于地震頻帶之外可有效反演大尺度的速度信息[8]。但是地震數(shù)據(jù)求取包絡(luò)后,以淺層直達(dá)波能量為主,包絡(luò)反演對淺層大尺度的速度信息恢復(fù)能力要優(yōu)于對中、深層速度信息恢復(fù)能力。

中、深層速度更新主要依賴于反射波信息,而反射波的觀測孔徑相對較小,幾乎只有高波數(shù)的速度得到了更新[9-10]。為了得到中、深層的背景速度,一方面從梯度優(yōu)化方面入手,對深層能量進(jìn)行補(bǔ)償以提高中、深層反演精度[11-12],另一方面,XU等[13]提出反射波全波形反演的方法旨在利用反射波更新中、深層的低波數(shù)速度參數(shù),并用格林函數(shù)表示反射波全波形反演中的雅可比矩陣及梯度。雅可比矩陣計算量大,故多采用伴隨狀態(tài)法求取梯度[14]。雅可比矩陣被稱為全波形反演的敏感核函數(shù)[15-18],表征速度擾動與波場擾動之間的關(guān)系。將常規(guī)全波形反演的敏感核函數(shù)分解為4個子核,不同子核表征不同位置不同尺度的速度更新量,反射波只有沿著“兔耳”路徑方可對中、深層的低波數(shù)速度進(jìn)行更新[19]。散射理論也可用于表征全波形反演對不同波數(shù)速度的更新,一般認(rèn)為全波形反演包含偏移分量和層析分量[20],對反射波數(shù)據(jù)而言,偏移分量能量遠(yuǎn)大于層析分量能量,而偏移分量對應(yīng)著高波數(shù)速度的更新,因此常規(guī)全波形反演的反射波對高波數(shù)速度的恢復(fù)能力更強(qiáng)。相同方向的正傳波場與殘差波場主要產(chǎn)生層析分量,用于更新低波數(shù)速度,相反方向的正傳波場與殘差波場主要產(chǎn)生偏移分量,用于更新高波數(shù)速度[21]。按照不同的方向分解波場,并更新偏移分量,再計算層析分量也可以達(dá)到反演不同波數(shù)速度參數(shù)的目的[22]。

以上方法均存在兩個問題:①包絡(luò)反演對地震數(shù)據(jù)變換,并未改變雅可比矩陣,且能量集中在淺層,對淺層背景速度的更新能力遠(yuǎn)大于對中、深層背景速度的更新能力;②反射波全波形反演方法可改變雅可比矩陣,且對中、深層背景速度更新能力更強(qiáng),但反演過程需要用到反射波,偏移/反偏移計算可以有效地得到反射波,但是每次迭代需要進(jìn)行6次正演計算,計算量巨大。因此本文提出一種基于包絡(luò)反演的高低波數(shù)同步反演方法,首先將線性初始速度作為一級速度,利用希爾伯特變換將地震數(shù)據(jù)變換為包絡(luò)數(shù)據(jù),重構(gòu)低頻數(shù)據(jù)用于淺層背景速度的恢復(fù);再將包絡(luò)反演的速度作為二級速度進(jìn)行高低波數(shù)同步反演,采用坡印廷矢量波場分解方法分解得到偏移分量及層析分量,根據(jù)偏移分量更新模型的擾動速度,同時根據(jù)層析分量更新模型中、深層的背景速度;最后利用Marmousi2模型數(shù)據(jù)的反演結(jié)果驗證了本文方法的有效性。

1 方法原理

1.1 包絡(luò)反演方法

全波形反演方法受周波跳躍問題的影響,易陷入局部極值。為解決低頻信息缺失情況下的周波跳躍問題,對地震數(shù)據(jù)進(jìn)行非線性變換以重構(gòu)低頻信息,從而為全波形反演提供合理的初始模型。地震數(shù)據(jù)包絡(luò)包含豐富的低頻信息,利用Hilbert變換可得到觀測數(shù)據(jù)和模擬數(shù)據(jù)的包絡(luò)信息,包絡(luò)殘差的L2范數(shù)即為包絡(luò)反演的目標(biāo)函數(shù),可表示為:

式中:u表示波場;?u/?m為雅可比矩陣;y表示模擬地震記錄;H{·}表示希爾伯特變換;yH表示希爾伯特變換后的模擬地震記錄。利用地震數(shù)據(jù)包絡(luò)中所包含的豐富的低頻信息可恢復(fù)介質(zhì)的背景速度,從而為常規(guī)全波形反演提供合理的初始速度模型。如圖1所示,包絡(luò)數(shù)據(jù)的能量集中于淺層的直達(dá)波,因此其只對淺層的低波數(shù)速度具有良好的恢復(fù)能力。

1.2 高低波數(shù)同步反演方法

在全波形反演中,速度模型包括背景速度和擾動速度:

式中:m0為背景速度,即低波數(shù)分量;Δm為擾動速度,即高波數(shù)分量。相應(yīng)的,地震波場可分解為背景速度對應(yīng)的地震波場及擾動速度對應(yīng)的地震波場,格林函數(shù)也可相應(yīng)地分解成兩部分:

圖1 地震數(shù)據(jù)與包絡(luò)數(shù)據(jù)的地震記錄對比(a)和頻譜對比(b)

式中:G為地震波場的格林函數(shù);G0為平滑背景波場的格林函數(shù),用于描述直達(dá)波和潛波;ΔG為擾動波場的格林函數(shù),用于描述反射波等。在一階Born近似假設(shè)下,全波形反演敏感核函數(shù)可表示為[15]:

式中:k(x;xs,xg)又稱為Fréchet導(dǎo)數(shù);u為地震波場;ω為圓頻率;m為速度場;G(x,t;xs)為震源波場的格林函數(shù);G(xg,t;x)為反轉(zhuǎn)波場的格林函數(shù)。因為格林函數(shù)可分解為兩部分,故常規(guī)全波形反演的敏感核函數(shù)可分解為4個子核[19]:

式中:k1(x;xs,xg)=ω2G0(x,t;xs)G0(xg,t;x);k2(x;xs,xg)=ω2G0(x,t;xs)ΔG(xg,t;x);k3(x;xs,xg)=ω2ΔG(x,t;xs)G0(xg,t;x);k4(x;xs,xg)=ω2ΔG(x,t;xs)ΔG(xg,t;x)。不同的子核函數(shù)表征不同位置不同尺度的速度更新量,核函數(shù)k1(x;xs,xg)表征直達(dá)波場之間的關(guān)系,主要依賴于背景速度模型,對于中、深層反射波而言,主要貢獻(xiàn)于高波數(shù)速度場[23];核函數(shù)k2(x;xs,xg)和k3(x;xs,xg)均表征背景波場及散射波場之間的關(guān)系,對沿著“兔耳”路徑即反射波路徑傳播的地震波場而言,即為層析算子的脈沖響應(yīng),對速度模型的低波數(shù)分量有貢獻(xiàn);核函數(shù)k4(x;xs,xg)表征散射波場之間的關(guān)系,主要依賴于擾動速度模型。分離直達(dá)波場和散射波場,即可得到各個子核函數(shù)。

利用坡印廷矢量[21,23]沿地震波傳播方向進(jìn)行地震波場分解,可得到正向直達(dá)波場和背向散射波場。利用坡印廷矢量可將波形反演的核函數(shù)分為偏移分量和層析分量:

式中:kt和km分別表示核函數(shù)的層析分量和偏移分量;GU和GD分別表示向上傳播地震波場和向下傳播地震波場的格林函數(shù)。圖2為層析分量與偏移分量,分別對應(yīng)恢復(fù)速度場中的背景速度與擾動速度。

由于全波形反演的初始速度是平滑的,只能產(chǎn)生較弱的擾動反射波場,因此導(dǎo)致了全波形反演中k1(x;xs,xg)的貢獻(xiàn)遠(yuǎn)大于k2(x;xs,xg)+k3(x;xs,xg)的貢獻(xiàn),也導(dǎo)致了中、深層低波數(shù)速度分量更新能力有限。

圖2 層析分量(a)與偏移分量(b)

圖3為高低波數(shù)同步反演流程,在全波形反演過程中,利用坡印廷矢量計算得到的層析分量和偏移分量,可以實現(xiàn)梯度分解,計算得到相應(yīng)的速度更新量,兩者疊加即為迭代過程中的速度更新量。高低波數(shù)同步反演方法所采用的偏移分量和層析分量同步迭代的方式能夠避免偏移分量和層析分量之間的耦合效應(yīng),在實現(xiàn)擾動速度更新的同時,加強(qiáng)了中、深層低波數(shù)背景速度的更新,一定程度上降低了對初始速度的依賴,提高了反演方法的穩(wěn)定性,反演結(jié)果較為理想。此外,反射波波形反演方法引入偏移/反偏移構(gòu)建層析分量[19],每次迭代需要額外進(jìn)行6次正演運算,其運算消耗及內(nèi)存消耗高;而基于坡印廷矢量的高低波數(shù)同步反演方法,無需偏移/反偏移計算即可得到層析分量,計算效率高。

圖3 高低波數(shù)同步反演流程

1.3 反演策略

本文基于包絡(luò)反演及高低波數(shù)同步反演方法提出了分步多尺度全波形反演策略。圖4為分步多尺度全波形反演策略流程:以線性初始速度作為一級初始速度,通過希爾伯特變換重構(gòu)低頻數(shù)據(jù),并進(jìn)行包絡(luò)反演以恢復(fù)淺層背景速度;將包絡(luò)反演結(jié)果作為二級初始速度,進(jìn)行高低波數(shù)同步反演。主要反演步驟如下:

1) 線性初始速度作為一級初始速度;

2) 將一級初始速度作為一級反演的初始輸入速度模型,采用包絡(luò)反演方法進(jìn)行速度迭代更新,主要實現(xiàn)淺層背景速度更新,得到一級速度反演結(jié)果;

3) 將一級速度反演結(jié)果作為二級反演的初始輸入速度模型,采用高低波數(shù)交替更新迭代的反演策略,實現(xiàn)擾動速度及中、深層背景速度更新,得到最終反演速度。

圖4 分步多尺度全波形反演策略流程

2 模型測試

利用圖5所示的Marmousi2速度模型進(jìn)行測試,模型參數(shù)如下:水平方向網(wǎng)格點數(shù)為301,垂直方向網(wǎng)格點數(shù)為151,水平方向和垂直方向網(wǎng)格間距均為8m,時間采樣間隔為0.5ms,時間采樣點數(shù)為1500,子波采用雷克子波,全波形反演子波主頻為25Hz,震源放置于地表,第1炮的坐標(biāo)為(56m,0),炮間距為56m,檢波點也放置于地表,第1個檢波點坐標(biāo)為(0,0),檢波點間隔為8m,地表全孔徑接收。

將從淺層到深層速度值線性增加的速度模型作為包絡(luò)反演的初始速度,如圖6所示。首先進(jìn)行包絡(luò)反演,地震數(shù)據(jù)的子波主頻為25Hz,對地震數(shù)據(jù)求取包絡(luò)后,頻率向超低頻方向移動,主要集中在0～5Hz附近;能量集中在淺層直達(dá)波處,對淺層的背景速度更新能力強(qiáng)。這一步的包絡(luò)反演為之后的全波形反演提供了背景速度更為準(zhǔn)確的初始速度模型。

圖5 Marmousi2速度模型

圖6 包絡(luò)反演的線性初始速度模型

包絡(luò)反演采用的子波主頻為10Hz,迭代30次后得到的速度模型如圖7所示,可以看出包絡(luò)算子可以解調(diào)出地震數(shù)據(jù)中的超低頻信息并可以有效恢復(fù)速度模型的背景速度信息。分別抽取真實速度模型、線性初始速度模型以及包絡(luò)反演速度模型中第120道,第150道和第200道的速度并對比,結(jié)果如圖8所示,可以看出,恢復(fù)的淺層背景速度接近真實速度,而中、深層速度的更新量較小,導(dǎo)致包絡(luò)反演對背景速度的更新能力集中于淺層。

圖7 包絡(luò)反演(迭代30次)得到的速度模型

將圖7包絡(luò)反演(迭代30次)得到的速度模型作為初始速度模型,進(jìn)行常規(guī)全波形反演,子波主頻25Hz,誤差迭代收斂后得到的速度模型如圖9所示。分別抽取真實速度模型和常規(guī)全波形反演速度模型中的第120道、第150道和第200道的速度進(jìn)行對比,結(jié)果如圖10所示,可以看出,以包絡(luò)反演速度作為初始速度進(jìn)行常規(guī)全波形反演可以較好地更新擾動速度,黃色橢圓區(qū)域所在的淺層(貼近地表處震源會影響干擾速度的更新)背景速度準(zhǔn)確,速度恢復(fù)結(jié)果幾乎與真實速度一致;藍(lán)色橢圓區(qū)域所在的中、深層,誤差收斂后的速度更新結(jié)果并不理想,這主要是由于常規(guī)全波形反演的速度更新能力集中于擾動速度,對中、深層背景速度的更新能力弱。

圖8 包絡(luò)反演(迭代30次)得到的單道速度與真實速度、初始速度對比a 120道; b 150道; c 200道

圖9 常規(guī)全波形反演得到的速度模型

同樣地,將圖7包絡(luò)反演(迭代30次)得到的速度作為初始速度,進(jìn)行高低波數(shù)同步反演,子波主頻25Hz,誤差迭代收斂后得到的速度模型如圖11所示。在反演迭代過程中,利用偏移分量對擾動速度進(jìn)行更新,同步引入層析分量對中、深層的背景速度進(jìn)行迭代更新,迭代40次后根據(jù)優(yōu)化的全波形反演結(jié)果確定反射界面,并利用坡印廷矢量分離出上行波進(jìn)而求取層析分量與偏移分量,結(jié)果如圖12所示。圖13 為采用不同方法反演得到的迭代誤差曲線,可以看出,高低波數(shù)同步反演方法具有良好的收斂性,隨著迭代次數(shù)的增加,誤差穩(wěn)定下降并收斂。分別抽取真實速度模型、高低波數(shù)同步反演速度模型以及常規(guī)全波形反演速度模型中第120道、第150道和第200道的速度進(jìn)行對比,結(jié)果如圖14所示。總體而言,常規(guī)全波形反演及高低波數(shù)同步反演得到的反演速度均與真實速度吻合;相對而言,應(yīng)用高低波數(shù)同步反演方法得到的反演速度與真實速度吻合得更好,且在中、深層區(qū)域(圖11的黑色橢圓區(qū)域和圖14的藍(lán)色橢圓區(qū)域)速度精度更高。

圖10 常規(guī)全波形反演得到的單道速度與真實速度對比a 120道; b 150道; c 200道

圖11 高低波數(shù)同步反演得到的速度模型

圖12 利用坡印廷矢量分解得到的層析分量(a)與偏移分量(b)

圖13 采用不同方法反演得到的迭代誤差曲線

圖14 采用不同方法反演得到的單道速度與真實速度對比a 120道; b 150道; c 200道

3 結(jié)論與認(rèn)識

我們利用坡印廷矢量進(jìn)行梯度分解并構(gòu)建了高低波數(shù)同步反演方法,同時結(jié)合包絡(luò)反演提出了分步多尺度反演策略。該策略引入包絡(luò)反演實現(xiàn)淺層背景速度場更新,然后將更新的背景速度場作為初始速度模型,利用波印廷矢量分解反演梯度的偏移分量和層析分量,交替更新小尺度擾動速度和中、深層處大尺度背景速度,以達(dá)到不同尺度速度同步更新的目的。Marmousi2模型數(shù)據(jù)測試結(jié)果表明,基于包絡(luò)反演的高低波數(shù)同步反演方法可以準(zhǔn)確穩(wěn)定地恢復(fù)中、深層速度,同時避免偏移/反偏移計算帶來的計算消耗和內(nèi)存消耗,計算效率較高。采用基于包絡(luò)反演的高低波數(shù)同步反演方法反演得到的速度模型,為之后地震層位標(biāo)定及油氣預(yù)測奠定了基礎(chǔ)?？紤]到實際資料中初始速度模型不精確,本文方法在三維實際資料中的應(yīng)用還需深入研究。