基于包絡反演的高低波數同步反演方法

張曉語,杜啟振,符力耘

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

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

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

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

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

1 方法原理

1.1 包絡反演方法

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

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

1.2 高低波數同步反演方法

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

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

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

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

式中:k(x;xs,xg)又稱為Fréchet導數;u為地震波場;ω為圓頻率;m為速度場;G(x,t;xs)為震源波場的格林函數;G(xg,t;x)為反轉波場的格林函數。因為格林函數可分解為兩部分,故常規全波形反演的敏感核函數可分解為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)。不同的子核函數表征不同位置不同尺度的速度更新量,核函數k1(x;xs,xg)表征直達波場之間的關系,主要依賴于背景速度模型,對于中、深層反射波而言,主要貢獻于高波數速度場[23];核函數k2(x;xs,xg)和k3(x;xs,xg)均表征背景波場及散射波場之間的關系,對沿著“兔耳”路徑即反射波路徑傳播的地震波場而言,即為層析算子的脈沖響應,對速度模型的低波數分量有貢獻;核函數k4(x;xs,xg)表征散射波場之間的關系,主要依賴于擾動速度模型。分離直達波場和散射波場,即可得到各個子核函數。

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

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

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

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

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

圖3 高低波數同步反演流程

1.3 反演策略

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

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

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

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

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

2 模型測試

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

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

圖5 Marmousi2速度模型

圖6 包絡反演的線性初始速度模型

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

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

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

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

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

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

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

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

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

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

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

3 結論與認識

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