反演近地表物性參數(shù)的地震層析成像方法綜述

劉玉柱,吳世林,劉偉剛,黃鑫泉,伍 正

(1.同濟大學(xué)海洋地質(zhì)國家重點實驗室,上海200092;2.同濟大學(xué)海洋與地球科學(xué)學(xué)院,上海200092)

1 層析算法

基于射線理論的地震走時層析成像在工業(yè)界獲得了成功應(yīng)用,現(xiàn)在仍是商業(yè)軟件的核心建模模塊。其基本原理建立在對正問題的一階Taylor級數(shù)展開的基礎(chǔ)上:

t=f(s)≈f(s0)+f′(s0)Δs

(1)

式中:t為射線的走時;f(s)為表達走時與慢度s之間非線性關(guān)系的函數(shù);s0為背景慢度。在射線理論框架下,正問題可以用矩陣公式表達為:

f(s)=L(s)·s

(2)

式中:L(s)表示慢度模型s下的射線路徑長度矩陣。根據(jù)Fermat原理,射線路徑是相對穩(wěn)定的,因此:

f′(s0)=L(s0)

(3)

將(3)式代入(1)式即得傳統(tǒng)的射線走時層析方程組:

L(s0)Δs=Δt

(4)

式中:Δs為慢度擾動向量;Δt為真實走時與理論合成走時之間的走時殘差向量。傳統(tǒng)層析成像通過求解上述方程組實現(xiàn)局部慢度向量的更新。如此,利用局部線性化+迭代方式實現(xiàn)非線性走時層析反演。

受計算能力的限制,20世紀(jì)70—80年代以來,地震層析成像一直采用經(jīng)驗性的求解大規(guī)模線性方程組(4)的方式進行模型的更新[1]。典型的算法包括反投影法(BPT)、代數(shù)重建法(ART)和聯(lián)合代數(shù)重建法(SIRT)等。這類方法基于加權(quán)平均的方式對模型進行更新,計算效率高但反演精度低。另一類反演類算法包括最小二乘QR分解(LSQR)和截斷奇異值分解(SVD)等,基于奇異矩陣的廣義逆算法,反演效果更好,但計算量比較大,算法不穩(wěn)定[2]。但無論經(jīng)驗性解法還是廣義逆解法都需要存儲層析矩陣,對計算機內(nèi)存要求很高。

借鑒全波形反演的實現(xiàn)過程[3],在目標(biāo)函數(shù)的局部最優(yōu)化基礎(chǔ)上[4]建立了一種走時層析成像實現(xiàn)方式,該方式通過尋找某一慢度向量s*,使得如下目標(biāo)函數(shù)極小:

(5)

式中:t0為觀測走時;tc為滿足(2)式的理論合成走時。結(jié)合(3)式得到上述目標(biāo)函數(shù)的梯度[5]:

(6)

不難推而廣之,只要采用(5)式的目標(biāo)函數(shù)形式,任何反演應(yīng)用都可以得到如(6)式的目標(biāo)函數(shù)梯度表達式,區(qū)別在于L取其一般形式為依賴于正問題的敏感核函數(shù)K,Δt取其一般形式為數(shù)據(jù)殘差Δd。因此,按照下列公式沿著最速下降方向(也可以采用其它更優(yōu)的更新方向)即可實現(xiàn)模型的更新。

sk+1=sk-αg

(7)

式中:k為反演中的迭代次數(shù);α為步長,可以利用線性搜索法計算得到[6]。不難看出,目標(biāo)函數(shù)梯度的計算是關(guān)鍵,理論上可以根據(jù)(6)式算得,但需要計算并存儲核函數(shù)矩陣,在實際應(yīng)用中不現(xiàn)實。

隨著伴隨狀態(tài)法在全波形反演(FWI)中的成功應(yīng)用[7],越來越多的學(xué)者將伴隨狀態(tài)法應(yīng)用到層析成像中[4,8-10]。伴隨狀態(tài)法走時層析將程函方程作為傳統(tǒng)目標(biāo)函數(shù)的約束項,利用拉格朗日乘子法進行求解,進而間接得到目標(biāo)函數(shù)梯度的表達式:

(8)

式中:v為地下介質(zhì)速度分布;λ為伴隨旅行時場,在地表處滿足:

(9)

在介質(zhì)內(nèi)部滿足:

(10)

式中:n為地表法向量。層析成像無需射線追蹤即可直接求得目標(biāo)函數(shù)的梯度,配合快速掃描正演算法,計算效率高,內(nèi)存占用小。但無法實施預(yù)條件,為此給出了一種等效的預(yù)條件算子[10]。

LIU等[11]在全波形反演研究中提出了一種改進的散射積分算法,之后又將該算法應(yīng)用于地震走時層析成像[12-13]。該算法與伴隨狀態(tài)法皆用于局部最優(yōu)化算法中梯度的求取,但改進的散射積分法將梯度的計算表達為核函數(shù)—殘差乘的累積求和形式:

pSD=-g=LTΔt=

(11)

式中:pSD為最速下降方向,即負(fù)梯度方向;lij為敏感核函數(shù)矩陣L中的元素,表示第i條射線在第j個模型網(wǎng)格里的長度;Δti為Δt向量中的元素,表示第i個數(shù)據(jù)的走時殘差。可見,該方法盡管需要計算核函數(shù),但無需存儲核函數(shù),內(nèi)存占用極小、計算效率高、物理意義明確。更重要的是,由于計算了核函數(shù),預(yù)條件與Hessian矩陣的利用非常方便。如利用公式(12),獲得梯度的同時就可以得到準(zhǔn)確的預(yù)條件算子,基本不需額外的計算量。這是伴隨狀態(tài)法層析無法實現(xiàn)的。

(12)

式中:H0為Hessian矩陣的對角線元素,即準(zhǔn)確的預(yù)條件算子。圖1展示了預(yù)條件算子對層析成像反演的影響,尤其當(dāng)初始模型梯度嚴(yán)重小于真實模型梯度時,無預(yù)條件的層析在淺層形成了高速屏蔽層,導(dǎo)致深部的速度結(jié)構(gòu)很難被反演出來。

圖1 預(yù)條件算子對層析成像的影響a 真實速度模型; b 初始速度模型; c 帶預(yù)條件層析反演結(jié)果; d 無預(yù)條件層析反演結(jié)果

2 基于不同理論的地震層析成像方法

根據(jù)(4)式,地震層析成像包含3個要素:核函數(shù)、物性參數(shù)與數(shù)據(jù)。這3個要素相互關(guān)聯(lián),某些情況下又可以相互獨立。核函數(shù)依賴于表達物性參數(shù)與數(shù)據(jù)之間關(guān)系的方程(正問題),因此根據(jù)所基于的正演理論不同,層析成像可分為射線層析、菲涅爾體層析與波動層析。傳統(tǒng)的地震層析成像基于射線理論[14-16],射線本身即為層析的核函數(shù)。由于其簡單、高效,在層析靜校正[17]、背景模型建立[18-19]中應(yīng)用廣泛。但受其高頻假設(shè)的限制,其反演結(jié)果分辨率不高[20]。

LUO等[21]提出了基于波動方程的走時層析成像方法。該方法利用波形互相關(guān)構(gòu)造走時殘差,然后利用伴隨狀態(tài)法求解最小二乘目標(biāo)函數(shù),迭代更新模型。隨后,MARQUERING等[22]提出了有限頻地震層析成像方法。該方法基于波動理論構(gòu)建核函數(shù),同樣利用波形互相關(guān)構(gòu)造走時殘差,在天然地震學(xué)研究中得到了成功應(yīng)用。但這些方法需要波動方程的正演模擬計算,計算量龐大。

為此,SPETZLER等[23]和LIU等[24]對走時菲涅爾體的性質(zhì)進行了研究,劉玉柱等[25]進而提出了帶限菲涅爾體地震層析成像方法。該方法基于波動理論顯式構(gòu)建核函數(shù)((13)式),但只考慮對單一震相貢獻最大的菲涅爾體范圍內(nèi)的波路徑(圖2)。

(13)

式中:ω為圓頻率;Im表示取復(fù)數(shù)的虛部;G(g,r)為無擾動速度場v0(r)中r點在檢波點g處的格林函數(shù);u(r,s)為炮點s點在空間任意一點r處的格林函數(shù)。

后續(xù)LIU等[13]又將改進的散射積分法引入梯度的計算,將射線走時殘差替代為波形互相關(guān)殘差,計算效率顯著提高,內(nèi)存占用量大幅降低,反演精度有所提高。菲涅爾體層析與射線層析的理論模型反演效果如圖3所示,對比圖3b和圖3c可見,反演中考慮波動理論核函數(shù)的確可以提高反演精度。圖4和圖5展示了西部山區(qū)某實際資料的散射積分法菲涅爾體層析與射線層析反演結(jié)果及其對應(yīng)的層析靜校正疊加剖面[13],可以看到菲涅爾體層析的優(yōu)勢明顯。

圖2 射線(黑色實線)、菲涅爾體(白色虛線)與走時波路徑(整個模型空間)

圖3 復(fù)雜起伏地表模型(a)以及射線層析(b)與菲涅爾體走時層析(c)反演結(jié)果

圖4 陸上實際資料的散射積分法菲涅爾體層析反演結(jié)果(a)與射線層析反演結(jié)果(b)

圖5 陸上實際資料的散射積分法菲涅爾體層析靜校正疊加剖面(a)與射線層析靜校正疊加剖面(b)

3 利用不同初至信息的地震層析成像方法

傳統(tǒng)的走時層析成像利用的是射線到達時,即初至走時信息。該走時與射線路徑核函數(shù)相匹配。隨著有限頻理論和波動理論在層析成像中的應(yīng)用,波形互相關(guān)走時替代了射線走時,被用于構(gòu)建走時殘差。但無論是射線走時還是互相關(guān)走時,都是單一走時,反演精度和分辨率提高有限。

為了提高反演分辨率,SHENG等[26]提出了初至波形反演,由于相對于單一走時信息,初至波形包含的信息量更多,因此反演精度和分辨率顯著提高。但該方法需要進行波動方程的正演模擬,計算量與傳統(tǒng)的全波形反演相當(dāng),因此難以應(yīng)用到實際當(dāng)中。劉玉柱等[27]提出了利用高斯束合成初至波形,以提高初至波形反演計算效率的方法,取得了一定的效果。

為了避開時間域波動方程的正演模擬,同時利用初至波形中更多的信息,CHOI等[28]提出了瞬時走時層析成像方法,劉玉柱等[29]提出了初至波相位走時層析方法。這些方法基于波動理論,利用初至中頻率依賴的相位信息實現(xiàn)初至波多頻率走時反演,反演精度介于射線層析與波形層析之間,但計算效率較波形反演明顯降低。圖6給出了利用射線走時、瞬時走時與初至波形的層析反演結(jié)果,可以看出隨著利用的信息量的增加,反演精度逐步提升。相對于射線層析,波形層析的反演精度有很大的提升。但值得一提的是,隨著信息量的增加,層析反演的非線性程度也在增強,即對初始模型的依賴性也逐次增強。因此,實際應(yīng)用中可以通過逐步增加信息量來逐級提高反演的精度。

圖6 利用不同信息的初至層析反演結(jié)果a Overthrust真實模型; b 常梯度初始模型; c 射線層析反演結(jié)果; d 瞬時走時層析反演結(jié)果; e 初至波形層析反演結(jié)果

圖7a展示了某海上OBS實際資料的初至波形層析反演結(jié)果。圖7b是與該反演速度無關(guān)的商業(yè)軟件疊前時間偏移(PSTM)剖面,T23軸和T25軸與圖7a 中深度為2.0km和4.7km的界面匹配得非常好[30]。

圖7 海上OBS實際資料的初至波形層析反演結(jié)果(a)以及某商業(yè)軟件海面拖攬數(shù)據(jù)PSTM剖面(b)

4 從單參數(shù)到多參數(shù)

傳統(tǒng)的層析成像方法以反演P波速度為主,但這已經(jīng)不能滿足成像、儲層識別等對地下介質(zhì)多參數(shù)反演的需求。因此,多參數(shù)層析成像研究逐漸多起來,包括各向異性層析[31-33]、Q層析[34-36]、縱橫波速度與密度層析[37-38]等。但多參數(shù)層析的主要難點在于多參數(shù)耦合,即不同的參數(shù)對觀測數(shù)據(jù)都有影響。雖然理論上可以導(dǎo)出并分析不同參數(shù)的敏感核函數(shù),但同時反演多參數(shù)仍然很困難,尤其存在強弱參數(shù)之分的情況下,弱參數(shù)往往很難被反演出來。研究表明,除利用必要的先驗信息作為約束外,還需要考慮以下3個方面。

1) 參數(shù)化方法[39-42]。采用合適的參數(shù)化方式可以有效降低多參數(shù)之間的強弱差異,有望實現(xiàn)多參數(shù)的同時反演。如VTI介質(zhì)走時層析成像中,速度屬于強參數(shù),各向異性參數(shù)屬于弱參數(shù),但如果采用三速度參數(shù)化方式則可以使得對應(yīng)的核函數(shù)在相同的數(shù)量級上,這對反演很有利[5]。

2) 反演策略[43-45]。盡管多個參數(shù)對地震數(shù)據(jù)都有影響,但不同的參數(shù)對不同的數(shù)據(jù)子集可能有不同的影響權(quán)重,采用針對性的反演策略往往可以收到事半功倍的效果。如劉玉柱等[5]在各向異性走時層析反演中采用了兩步法反演策略,即第一步同時反演速度和各向異性參數(shù),然后將反演的速度和最開始的各向異性參數(shù)作為初始模型再進行第二輪反演,得到可靠的各向異性參數(shù)(如圖8所示)。

3) Hessian矩陣的利用[11,46-49]。Hessian矩陣是目標(biāo)函數(shù)對模型參數(shù)的二階偏導(dǎo)數(shù)構(gòu)成的系數(shù)矩陣。在單參數(shù)反演中,Hessian矩陣的對角線元素反映了波場的幾何擴散效應(yīng),非主對角線元素反映了空間不同散射體對地震波傳播的干涉效應(yīng);在多參數(shù)反演中,如果將Hessian矩陣視為塊狀矩陣,那么Hessian矩陣的主對角塊元素反映的是單參數(shù)對波場的以上影響,非主對角塊元素反映的則是包含幾何擴散和干涉效應(yīng)在內(nèi)的多參數(shù)耦合作用。如在變密度聲波方程全波形反演中,Hessian矩陣可以表達為[50]:

(14)

式中:Kv,Kρ分別表示速度和密度的核函數(shù)矩陣;T表示復(fù)矩陣的共軛轉(zhuǎn)置。因此,在多參數(shù)反演中考慮Hessian矩陣可以實現(xiàn)多參數(shù)的同時反演,提高多參數(shù)的反演精度。但遺憾的是,目前在近地表物性參數(shù)層析反演(尤其是多參數(shù)反演)中較少使用完整的Hessian矩陣。

圖8 VTI介質(zhì)兩步法多參數(shù)走時層析反演a 二維復(fù)雜起伏地表模型地表以下200m深度范圍內(nèi)的真實速度結(jié)構(gòu); b 第一輪反演結(jié)果(地表以下200m深度范圍內(nèi)的速度結(jié)構(gòu)); c 第二輪反演結(jié)果(地表以下200m深度范圍內(nèi)的速度結(jié)構(gòu)); d 二維復(fù)雜起伏地表模型地表以下200m深度范圍內(nèi)的ε模型; e 第一輪反演結(jié)果(地表以下200m深度范圍內(nèi)的ε值); f 第二輪反演結(jié)果(地表以下200m深度范圍內(nèi)的ε值)

5 結(jié)論與展望

關(guān)于地震層析成像的研究持續(xù)了三十多年,從射線理論到波動理論,從走時反演到波形反演,從單參數(shù)到多參數(shù),從方程求解到局部最優(yōu)化,對其理論、方法與應(yīng)用研究從簡單到復(fù)雜,反演精度與分辨率逐步提高。從以上綜述不難看出,影響反演精度最主要的因素是利用的信息,波形反演要遠(yuǎn)比走時反演獲得的反演分辨率高;其次是理論基礎(chǔ),由于地震波傳播的有限頻特性,相較于射線層析,波動理論反演可以獲得更高的反演精度;最后是反演算法、反演策略、參數(shù)化等對反演質(zhì)量也具有一定的影響。

多參數(shù)反演是目前層析成像的研究熱點,除Hessian矩陣的有效利用外,反演策略的制定也很關(guān)鍵,值得深入研究。另外,本文提出的改進的散射積分算法方便預(yù)條件、Hessian矩陣?yán)贸浞?、?nèi)存占用小、反演穩(wěn)定等諸多優(yōu)點使其尤其適合于求解地震層析成像反演問題。相信它會有更加廣泛的實際應(yīng)用價值。