胡光輝,賀偉光

(中石化石油物探技術研究院有限公司,江蘇南京211103)

全波形反演(FWI)將疊前地震數(shù)據(jù)信息映射為地下結構或介質(zhì)參數(shù)分布,是目前地震學中最強大的建模成像工具之一。在勘探地球物理學中,早在1975年就建立了從地面地震數(shù)據(jù)推斷內(nèi)部結構的概念[1-2],但FWI的理論建立者仍被認為是LAILLY[3]和TARANTOLA[4]。對于全波形反演一詞的理解存在兩個觀點:第1個觀點堅持認為只有利用整個地震炮集數(shù)據(jù)的波形信息才能稱為全波形反演(full指的是整個地震炮集數(shù)據(jù));第2個觀點更寬泛一些,把所選擇的地震事件信息完全利用即可成為全波形反演(full暗含的意思是fully)。無論是哪種觀點,都認為波的傳播需要經(jīng)過求解波動方程來模擬,利用其它退化的模擬方法(比如射線方程法、單程波方程法)都不能稱為全波形反演。在全波形反演的早期發(fā)展中,由于計算能力有限以及缺乏低頻和長偏移距觀測數(shù)據(jù),因而其發(fā)展十分緩慢。自2000年以來,由于全波形反演理論上的發(fā)展和計算能力的進步,該方法逐步發(fā)展成為現(xiàn)代工業(yè)生產(chǎn)中的常規(guī)地震資料處理方法之一,尤其是海上地震資料。然而,無論是在全波形反演的早期發(fā)展階段還是如今的工業(yè)應用階段,周跳(cycle skipping)問題依然沒有完全解決[5-6]。盡管有一些精巧的反演策略[7](頻率逐漸增加和偏移距逐漸增加),但是傳統(tǒng)的最小二乘函數(shù)(L2函數(shù))在許多情況下依然會陷入局部極值。

從數(shù)學上講,全波形反演是在彈性動力學方程約束下的局部優(yōu)化問題。數(shù)學上的處理有兩種:第1種思路是構建拉格朗日函數(shù);第2種思路是建立增廣拉格朗日函數(shù)。增廣拉格朗日函數(shù)是將彈性動力學方程作為二次項來構建優(yōu)化函數(shù)[8]。從數(shù)學界來看,增廣拉格朗日函數(shù)是凸性良好的優(yōu)化函數(shù)。波動方程以二次型的方式引入增廣拉格朗日函數(shù)并不要求波動方程被精確地滿足。這種處理方式相當于增大了波形反演的解的搜索空間,同時增強了方程的凸性。但是,增廣拉格朗日函數(shù)的計算要求很高,原因是涉及模擬算子和伴隨算子的乘積。在頻率域,計算代價可以接受,但在時間域則相當困難。因此,增廣拉格朗日函數(shù)在工業(yè)生產(chǎn)中很少使用。全波形反演的主流框架一直是第1個思路,PLESSIX[9]對該框架進行了綜述。在拉格朗日框架中,作為約束的彈性動力學方程被線性地引入到優(yōu)化函數(shù)中,之后可以方便地得到伴隨場震源和梯度的計算。拉格朗日框架要求在每次迭代時都能準確地求解彈性動力學方程。等價地說,在全波形反演的整個迭代過程中,約束方程(彈性動力學方程)始終得到滿足。這限制了未知參數(shù)的搜索空間,因此導致拉格朗日框架的凸性不如增廣拉格朗日框架。但是,拉格朗日框架因為計算代價較低而成為學術界和工業(yè)界的主流方式。在拉格朗日函數(shù)框架方案中,避免周跳的主要方法是設計凸目標函數(shù),如L1范數(shù)[10]、混合或Huber函數(shù)[11-13]、零延遲互相關[14-17](相當于歸一化的L2函數(shù))、包絡函數(shù)[18-24]、懲罰互相關函數(shù)[25]、基于反褶積的函數(shù)[26]和最優(yōu)傳輸函數(shù)[27-33]。其中,最優(yōu)化傳輸函數(shù)的優(yōu)勢是其比較的是全局地震數(shù)據(jù),而不是逐點比較地震數(shù)據(jù)。這種特性保持了地震數(shù)據(jù)的事件間和檢波點端記錄數(shù)據(jù)的相干性,從而使得較傳統(tǒng)L2函數(shù)能更穩(wěn)定地反演。本研究采用拉格朗日函數(shù)的框架,綜述了幾種最優(yōu)傳輸函數(shù)在全波形反演中的應用。

最優(yōu)化傳輸起源于法國的一個工程問題,即如何將沙子從源頭位置運輸出來,以最少的能耗來填補施工區(qū)域的洞。這個問題被MONGE[34]用數(shù)學方法表述為一個L1最小化問題。在一些研究中,它可能被命名為EMD(Earth Mover Distance)。方程的答案可以通過求解非線性方程,這種非線性方程也稱為Monge-Ampere方程。將1范數(shù)推廣到任何其它范數(shù)稱為p-Wasserstein函數(shù),其中,p表示范數(shù)。證明Monge方程解的存在性,以及求解都是很困難的。KANTOROVICH[35]取得了突破,將Monge公式擴展到更高維空間會得到線性問題;具體來說,這是線性分配問題。FERRADANS等[36]和LELLMANN等[37]在圖像處理中引入最優(yōu)化傳輸函數(shù),該函數(shù)是一種自動且正確的匹配模式。最優(yōu)化傳輸函數(shù)目前已廣泛應用于機器學習或人工智能(AI)等領域。最優(yōu)傳輸問題根據(jù)分布是連續(xù)的還是離散的[38]分為3類,即連續(xù)(兩個分布都是連續(xù)的)、半離散(一個是連續(xù)的,另一個是離散的)和離散的(兩者都是離散的)。連續(xù)問題可以用Monge-Ampere方程或流體動力學方程來解決。半離散問題利用計算幾何和牛頓優(yōu)化技術可以快速解決。離散問題用正則化最優(yōu)傳輸求解器或標準線性賦值求解器求解。

直接計算最優(yōu)傳輸問題的運算量巨大。例如,當將n個點移動到n個目標位置時,每個點都可以移動到任何位置。因此,未知變量的數(shù)目為n2。KANTOROVICH[35]將原始的最優(yōu)化傳輸問題轉(zhuǎn)換到對偶空間中,未知量的個數(shù)變?yōu)閚,極大地減少了離散情形下未知變量的個數(shù)。在對偶空間中,這可以用標準的線性方程求解器來方便地求解。同時,也有一些專門針對線性分配的算法加速收斂,比如匈牙利算法[39]或拍賣算法[40]。在圖像比較領域,在傳輸代價函數(shù)中添加關于傳輸算子的正則化被廣泛采用[41]。為了提高計算效率,人們發(fā)展了更為精巧和復雜的計算技術,如多層次采樣策略、多尺度正則化策略。在多層次采樣策略中,將整個數(shù)據(jù)集采樣到多個層次。將原始數(shù)據(jù)集定義為零級,沿每個維度以2的比率進行降采樣得到第1級(對于2D問題,這個數(shù)字現(xiàn)在是原始數(shù)據(jù)集的0.25倍)。第1級再次重新采樣以獲得下一級,最后一級是最小的數(shù)據(jù)集合。數(shù)據(jù)分級采樣之后,求解器從最小的數(shù)據(jù)集開始(也是最終的數(shù)據(jù)集)。最后一級的傳輸平面計算出來之后,進行插值,成為前一級的初始會傳輸平面。這個過程一直持續(xù)到零級。在每個層面上,都采用多尺度的正則化策略來計算傳輸平面。在各級計算中,正則化方法相應地從相對較大的值開始,逐漸減小到接近于零。在傳輸平面代表的映射中,每個點會被分割到幾個位置上。在零級數(shù)據(jù)集上,傳輸平面上會有一些非常集中的能量分布,每個能量團的中心對應于一個一一映射的解。正則化的引入使得無法得到單點對單點的映射,但是可以定義為質(zhì)心坐標強制傳輸平面成為一對一的映射。

最優(yōu)化傳輸并不能直接應用到地震數(shù)據(jù)處理中。地震信號是上、下震蕩,有正、有負的。而最優(yōu)化傳輸處理的是正信號,且兩個比較對象的質(zhì)量要相等。因此需要對數(shù)據(jù)進行預處理,使得地震信號滿足非負且質(zhì)量相等的條件。依據(jù)數(shù)據(jù)預處理的方式,目前有3種最優(yōu)化傳輸公式在全波形反演中得到了很好的應用,即2-Wasserstein函數(shù)[27-29],KR-OT函數(shù)[30-32]和圖空間的OT函數(shù)GSOT[33]。2-Wasserstein函數(shù)是基于范數(shù)2準則的最優(yōu)傳輸函數(shù)的公式。它是第一個被引入到地球物理學的最優(yōu)傳輸函數(shù)[27](W22函數(shù),第一個2意味著2范數(shù),第二個2意味著平方)。它是一維的算法,每次對比一對地震道信號。每一個合成地震道和觀測地震道都應該經(jīng)過預處理為正的并歸一化(保證質(zhì)量相等)。最簡單的預處理是減去最小值并除以其積分。該公式的獨特優(yōu)點是存在解析解,計算速度快且穩(wěn)定。擴展到其它規(guī)范(p-Wasserstein函數(shù))是很直接的。缺點是預處理對理論凸性略有破壞。在極端情形下,可以假設地震信號減去負無窮大,在此情形下,最優(yōu)化傳輸函數(shù)退化為常規(guī)的最小二乘誤差。

第2個引入地球物理學的最優(yōu)化傳輸是由METIVIER等[30]提出的KR-OT。它是基于1范數(shù)(1-Wasserstein函數(shù))的最優(yōu)化傳輸函數(shù),并且是在對偶空間中表示。KR-OT在一維、二維和三維空間中都具有相似的數(shù)學結構。二維全波形反演的案例表明,2D KR-OT明顯優(yōu)于1D KR-OT。正如METIVIER所解釋的那樣,1D KR-OT反演的退化是忽略了接收器間的相干性。在KR-OT的公式中,最優(yōu)化傳輸函數(shù)只依賴于兩個數(shù)據(jù)殘差(地震數(shù)據(jù)可以直接發(fā)送到求解器,而不需要移動或歸一化),這是優(yōu)點也是缺點。弊端在于所有的地震事件從一開始就混合在一起,當初始合成數(shù)據(jù)與觀測數(shù)據(jù)相差較遠時,映射可能無法正確匹配地震事件。KR-OT的局限之一是采樣在時間軸和空間軸上應該是均勻的。至于將KR-OT擴展到其它領域,如超聲成像[42],原則上是不合適的,因為接收器并不是均勻分布的。

上述兩個最優(yōu)化傳輸公式都屬于連續(xù)傳輸問題?？碧降厍蛭锢韺W中的第3種最優(yōu)傳輸函數(shù)是圖空間最優(yōu)傳輸函數(shù)GS-OT,這是離散傳輸問題。將一維離散信號擴展為時間-振幅二維圖中的點集[33,43-44]。這種處理方法自然滿足了數(shù)據(jù)非負的要求。同時,如果用相同數(shù)量的點離散地震數(shù)據(jù),質(zhì)量守恒也是自然滿足的。從原理上將沿著任何一個維度都可以任意采樣,去除了均勻采樣的要求。GS-OT僅取決于點的數(shù)量,即在一個地震道中的離散采樣點。所以直接擴展到高維空間必然面對計算量巨大的挑戰(zhàn)。

最近,XU等[45]開發(fā)了基于圖空間公式的高維最優(yōu)傳輸函數(shù)。雖然這項研究是在全波形反演中的應用,但作者最初的研究動機是要匹配兩個彩色點云。每個顏色點都與一個向量(R,G,B)相關聯(lián),但最優(yōu)化傳輸只針對標量函數(shù)。從數(shù)學上講,XU等[45]開發(fā)的是比較兩個向量場的最優(yōu)化傳輸函數(shù)。該研究進行了二維和三維實驗。正如作者解釋的那樣,當點的數(shù)量少于1000時,計算非?？?但對于超過10000個點的彩色點云來說,計算時間過長。因為該研究將整個地震圖分成互相不重疊的小區(qū)域(二維地震圖是長方形,三維是長方體),大幅度降低了計算成本。依據(jù)在于,地震圖的對比中,合成數(shù)據(jù)的地震事件是不應當與觀測事件相差太遠。所以劃成小區(qū)域后,數(shù)據(jù)匹配依然存在合理性。但是在圖像對比中,整張圖可能會旋轉(zhuǎn)或者倒轉(zhuǎn),在這種情形下不能劃分小區(qū)域。經(jīng)過數(shù)據(jù)區(qū)域劃分處理后,將高維最優(yōu)傳輸函數(shù)納入全波形反演中變得可行和現(xiàn)實。該研究還將熵變換應用于空間坐標和時間坐標。雖然在一些模型中是否采用熵變換并不影響反演,但有趣的是,在Chevron 2014盲測中,應用熵變換的反演結果優(yōu)于沒有熵變換的反演。

最優(yōu)化傳輸函數(shù)是在原始空間或?qū)ε伎臻g中求解。在數(shù)學領域,一種新的基于混合算法的最優(yōu)傳輸方法被提出。該方法是基于2-Wasserstein函數(shù)的二維算法,不再是一維單道處理。這是由JACOBS等[46]開發(fā)的一種新穎的方法來解決最優(yōu)化傳輸函數(shù)。首先,將原始形式轉(zhuǎn)換到對偶空間,在數(shù)學上就得到了兩個等價的方程:一個在原始空間,另一個在對偶空間。新穎的部分在于第2步。JACOBS等[46]只迭代一步原始解;然后跳到對偶空間迭代一步;之后再回到原始空間繼續(xù)循環(huán),直到得到一個收斂的結果。JACOBS等[46]將其命名為往復法(back and forth method,BFM)。

最初的最小二乘目標函數(shù)是逐點比較的算法,后來興起的從互相關到最優(yōu)傳輸?shù)哪繕撕瘮?shù)更多是基于全局比較算法。因此,可以更好地利用相鄰點之間甚至地震事件之間的相干性。這很好地提高了目標函數(shù)的凸性和地下結構的一致性。本研究旨在提供對全波形反演中最優(yōu)化傳輸函數(shù)的介紹,解釋了目前引入全波形反演的幾種最優(yōu)化傳輸公式,包括W2(1D版本,2-Wasserstein函數(shù))、BFM、KR-OT、GS-1D和GS-2D(高維最優(yōu)傳輸函數(shù))。

1 方法理論

本文重點在于討論最優(yōu)化傳輸函數(shù),但為了說明與全波形反演的結合應用,先簡單介紹全波形反演的實現(xiàn)。之后重點介紹5個最優(yōu)化傳輸函數(shù)的基本原理。W2是1D 2-Wasserstein函數(shù),BFM是2D或者3D的2-Wasserstein函數(shù),KR-OT是2D或者3D的1-Wasserstein函數(shù),GS-1D是一維信號的圖空間最優(yōu)化傳輸,GS-2D是多維信號的圖空間最優(yōu)化傳輸。

1.1 全波形反演

FWI已經(jīng)發(fā)展了幾十年,現(xiàn)在是地震資料處理業(yè)界的常規(guī)技術。FWI包含一個目標函數(shù)和一個約束方程(又稱為控制方程)。目標函數(shù)是刻畫合成地震數(shù)據(jù)與觀測地震數(shù)據(jù)直接偏差的函數(shù),控制方程是描述聲波或彈性地震波傳播的方程。最早也是最常見的目標函數(shù)是最小二乘函數(shù)(L2函數(shù)),即:

(1)

式中:p(t)是合成數(shù)據(jù);d(t)是觀測數(shù)據(jù),求和針對所有的震源和接收器。目標函數(shù)中應用了加權函數(shù)W(t),因為并不是所有的地震事件都被利用。例如,在應用聲波方程求解器處理陸地地震數(shù)據(jù)時,需要預先去除觀測數(shù)據(jù)中的面波和橫波。目標函數(shù)的優(yōu)化通常是基于局部優(yōu)化求解器,模型更新的搜索方向是基于梯度的。梯度的計算通常是采用伴隨方[9]。正傳波場與伴隨波場之間的累積求和為梯度。計算梯度只需要花費2～3個波傳播的計算時間。伴隨波場是通過反向傳播伴隨源,伴隨源由目標函數(shù)關于合成數(shù)據(jù)的一階導數(shù)給出。目標函數(shù)的一階擾動為:

(2)

由此可確定伴隨源s(t)=W(t)[p(t)-d(t)]。確定了伴隨源之后,即可通過正傳波場與反傳波場的相關求取模型更新的梯度,繼而通過不同的迭代求解算法逐步更新模型。

1.2 最優(yōu)化傳輸函數(shù)的模式匹配特性

最優(yōu)傳輸函數(shù)基于模式的匹配。這個性質(zhì)可以用插值的例子來說明。圖1在左下角顯示圓盤,在右上角顯示正方形盤。傳統(tǒng)的插值(線性插值)將它們與兩個加權因子分別相乘,之后相加。圖2中的第1行顯示了線性插值的演化,其中圓盤上的加權因子從1逐漸減小到0,同時正方形盤上的加權因子從0增加到1?？梢钥吹?插值的效果是此消彼長的變化。但有時,希望可以像水流一樣把物體形態(tài)光滑地變化成為另一個。圖2中的第2行顯示的正是這種變化,這一種基于最優(yōu)化傳輸?shù)牟逯?首先求出圓盤與正方形盤直接的傳輸平面,然后對平面插值即可?？梢钥吹?圓盤的形狀和位置均平滑地演化成正方形盤和其所在的位置。

最優(yōu)傳輸函數(shù)的這種插值特性對FWI非常有幫助。傳統(tǒng)的FWI經(jīng)常被周跳問題困擾。最小二乘函數(shù)只是點對點地比較兩個地震道。一旦兩個地震事件相距較遠,在伴隨場震源上兩個事件就不重疊,也就是說,它們之間就沒有相互作用了。這就解釋了為什么最小二乘函數(shù)只在兩個地震事件誤差在半周期內(nèi)時起作用,同時在誤差大于半周期的時候會陷入局部極值。而波動方程層析建模的凸性較FWI更好,因為兩個事件總是可以通過互相關聯(lián)系起來。因此,我們期望具有全局對比功能的最優(yōu)化傳輸函數(shù),它比最小二乘函數(shù)具有更好的凸性。

1.2.1W2函數(shù)

(3)

應遵守質(zhì)量守恒原則,即:

(4)

圖3 計算W2的示意

基于方程(3)和方程(4),可以計算出伴隨源。伴隨源s(t)由W2相對于合成的p(t)的偏導數(shù)給出[47]。引入中間變量,有:

(5)

最后的伴隨源為:

(6)

數(shù)值實現(xiàn)過程可總結為:

1) 對每個地震道進行預處理,使其為正。建議將每個處理過的地震道除以其自身的質(zhì)量,使總質(zhì)量為1。

移動和歸一化屬于線性變換,把數(shù)據(jù)變?yōu)檎盘栠€包括非線性變換。非線性變換包括平方、指數(shù)或雙曲正切(tanh)。CHEN等[47]在地震定位的研究中采用平方法對地震數(shù)據(jù)進行預處理。作者設計了一個兩層模型,并根據(jù)模型層位速度的不同計算了相對應的合成地震圖,最終做了目標函數(shù)的凸性分析,如圖4 所示。令人驚訝的是,將地震數(shù)據(jù)進行平方這個預處理優(yōu)于移動和歸一化策略。采用平方預處理得到的目標函數(shù)分布圖更加光滑,凸性更好。

圖4 目標函數(shù)誤差分布[47]

與此同時,體波的振幅和弱事件的振幅會得到增強。振幅平衡是彈性全波形反演中重要的步驟之一,對于提升體波對于深部的反演貢獻巨大。HE等[31-32]在各向異性介質(zhì)中詳細論述了振幅問題。

1.2.2BFM函數(shù)

一維問題的最優(yōu)化傳輸存在解析解。與一維問題不同,高維問題不存在最優(yōu)化傳輸?shù)慕馕鼋?。多年?數(shù)學家們一直在努力尋找有效的求解方法。來回跳躍法BFM屬于有效求解的方法之一。首先,對原問題構造拉格朗日函數(shù),將其轉(zhuǎn)化為對偶空間的等價問題。至此,有兩種等價的最優(yōu)化傳輸問題,一種在原始空間,另一種在對偶空間。BFM的獨特優(yōu)點是,最優(yōu)解是在原始空間和對偶空間之間迭代產(chǎn)生的。在原始空間迭代一次,然后跳到對偶空間,在對偶空間同樣只迭代一次,然后再跳回原始空間;一直循環(huán)直到問題收斂。這種方法的數(shù)學基礎比較復雜[46]。BFM仍然要求輸入是非負的,用于W2的預處理也適用于BFM。

1.2.3KR-OT函數(shù)

KR-OT是基于1-Wasserstein的函數(shù);但并不是求解原始問題,而是轉(zhuǎn)化到對偶空間中求解。KR-OT的代價函數(shù)是基于1范數(shù)c(x,y)=|x-y|,在對偶空間中的表達式為[30]:

(7)

式中:φ(x)在有界的1-Lipschitz空間中,而λ是一個用戶定義的參數(shù)。第1個不等式保證了φ不能太大的要求,第2個不等式排除了突然變化的情況。HE等[32]仔細研究了不等式約束的影響,其結論是存在較大的λ取值空間,都可以給出合理的解決方案φ。KR-OT的顯著特性是,該公式只依賴于兩個分布之間的差異,所以不要求輸入的地震信號是正的。對于任何負分布或非正分布,可以將兩個分布加上同一個常數(shù)變?yōu)檎?但是從KR-OT的表達式看,這個常數(shù)是會被消去而不起任何作用。

在全波形反演中,地震數(shù)據(jù)是以離散形式存儲的。在三維情況下,對應的KR-OT公式為:

(8)

式中:常數(shù)因子Δx1,Δx2,Δx3在公式(8)中的第1個方程中被省略,uijk是合成數(shù)據(jù)與觀測數(shù)據(jù)的差異。公式(8)可以用SDMM(simultaneous direction method of multipliers)高效求解。METIVIER等[30]詳細說明了如何解決KR-OT。其中,求解矩陣逆是一個沉重的計算負擔。后來,METIVIER等[48]證明了其中的矩陣運算與泊松方程是等價的,從而可以在頻率域用FFT方法或在時域用多重網(wǎng)格方法有效地求解泊松方程。雖然如此,三維計算仍然非常昂貴,主要原因是每次求解KR-OT都需要多次求解三維泊松方程。在二維空間中,從公式(8)中去掉一維,可以大大減少計算時間。因此,偽三維處理也是折衷的方法。將三維地震數(shù)據(jù)體沿地震數(shù)據(jù)線方向分割成切片(地震數(shù)據(jù)切片形成二維矩陣),并使用2D KR-OT求解器對每個切片進行處理。

KR-OT的輸入信號是殘差,故所有的地震事件從一開始就是混合的。這會影響KR-OT的凸性。與其它OT公式相比,KR-OT的獨特優(yōu)勢是KR-OT將解約束為1-Lipschitz函數(shù)。這種優(yōu)勢對于包含面波的陸地地震數(shù)據(jù)具有很大的潛力。圖5為KR-OT能平衡振幅,展示的是一個陸地數(shù)據(jù)全波形反演的伴隨場。L2伴隨源以面波為主,因此,重建由體波采樣的深層區(qū)域是困難的。相比之下,KR-OT伴隨源的振幅分布平衡,反演確實顯示了從淺到深的重構模型。

圖5 KR-OT能平衡振幅[32]

1.2.4GS-1D函數(shù)

圖空間最優(yōu)化傳輸(Graph Space Optimal Transport)處理的是離散化的問題。將一維地震道擴展為時間-振幅的二維空間,最優(yōu)化傳輸?shù)拇鷥r函數(shù)為:

cij=(pi-dj)2+η(ti-tj)2

(9)

式中:pi是在ti時刻采樣的合成地震道;dj是在tj時刻采樣的觀測地震道;cij是點(ti,pi)至點(tj,dj)的移動成本;η是用戶定義的參數(shù),用來控制時間軸的權重。最優(yōu)化傳輸要求質(zhì)量守恒,在圖空間中,質(zhì)量守恒約束指的是合成地震道與觀測地震道包含相同數(shù)量的點。配對一對地震道的總成本是對所有個體成本的求和,即:

(10)

(11)

伴隨源si看起來像傳統(tǒng)的L2的情況,但其中的觀測數(shù)據(jù)是重排列之后的觀測數(shù)據(jù)。在L2的情況下,兩個不重疊地震事件是沒有相互作用的,這也是周跳的根源。GS 1D依然能夠通過求解最優(yōu)的排列來聯(lián)系這些事件。

圖6顯示了L2、KR-OT和GS-1D的伴隨源。圖6a中,觀測數(shù)據(jù)包含多個地震事件,而合成數(shù)據(jù)只有兩個。圖6b中,常規(guī)的最小二乘函數(shù)無法平衡振幅,而KR-OT和GS-OT都在很好的平衡振幅。輸入信號由Ricker子波組成,模擬的是事件錯位和事件缺失的情況。L2伴隨源具有較大的振幅差異,而OT伴隨源確實平衡了振幅。然而,這并不意味著在這兩個OT函數(shù)中都忽略了振幅信息。在全波形反演中,隨著迭代的繼續(xù),觀測數(shù)據(jù)的振幅會得到越來越好的利用。最后,所有的事件都是完全擬合的。

圖6 L2、KR-OT和GS-OT的伴隨源[31]

1.2.5GS-2D(高維OT)

(12)

圖7和圖8均引自XU 等[45]。圖7顯示了兩個二維質(zhì)量分布之間的映射;兩個質(zhì)量分布都是由高斯函數(shù)生成,但是中心位置不同。二維的GS-OT算法很好地求解出了匹配關系(黑色實線)。若是采用一維算法,則只能沿著橫向或者縱向匹配,而無法正確求解出真實的匹配關系。在全波形反演中,地震事件經(jīng)常是在時間軸方向和空間軸方向都有移動,這也是為何要發(fā)展高維OT算法的原因。圖8顯示的是三維質(zhì)量分布之間的三維映射。兩個球體分別對應的是兩個三維的質(zhì)量分布,可以看出,三維的GS-OT算法正確地匹配到了兩個三維質(zhì)量分布。XU 等[45]將高維最優(yōu)化傳輸GS 2D應用于聲學和彈性全波形反演試驗,包括彈性SEG/EAGE逆沖模型和Chev-ron 2014的數(shù)據(jù)盲測試驗。在逆沖模型中,L2反演重構出了橫波速度參數(shù)vS的淺層結構,GS 1D重建了vS的淺部和深部。相應地,體波和面波都比L2得到了更好的利用。GS-2D得到了最佳的反演,即vP和vS都顯示了高分辨率和高保真度的結構。在XU等[45]展示的Chevron 2014數(shù)據(jù)盲測實驗中,使用熵變換的代價函數(shù)x=xlogx得到了更干凈的速度。

圖7 二維質(zhì)量分布之間的最優(yōu)化映射

圖8 三維質(zhì)量分布之間的最優(yōu)化映射

2 反演實驗

我們已經(jīng)實現(xiàn)了上述所有的最優(yōu)化傳輸函數(shù)W2,BFM,KR-OT,GS-1D和GS-2D的FWI。本節(jié)將展示兩個模型測試結果。第1個是透射例子;第2個是選取的SEAM Ⅱ Foothill模型的切片。

2.1 透射模型

真實模型是在均勻背景vP=3.0km/s的模型中添加了vP=3.6km/s的圓形異常體。圓形異常體的直徑為3.6km。假設密度為常數(shù)ρ=1g/cm3。以均勻背景速度作為初始模型。圖9顯示了真實的模型和初始的模型。在頂部等距離布置了32個爆炸源,在底部布置了300個壓力接收器。震源子波為10Hz的Ricker子波。以10Hz為參考計算,圓形異常體的大小是主波長的10倍。

圖9 全波形反演測試

圖10顯示了觀測數(shù)據(jù)以及觀測數(shù)據(jù)與合成數(shù)據(jù)的交錯比較。對比圖清楚地說明了周期跳躍問題。L2反演進行了6次迭代陷入了局部極值,反演結果在圖11a中。KR-OT反演也被局限在某個局部最小值中(圖11d)。在這個特定的例子中,初始模型與實際模型相差太遠。正如作者所解釋的,KR-OT將L2誤差作為輸入,它的性能是受到影響的。其它OT反演(W2,BFM,GS-1D,GS-2D)都正確地恢復了模型。

圖10 觀測數(shù)據(jù)、合成數(shù)據(jù)與觀測數(shù)據(jù)對比

2.2 SEAM Ⅱ Foothill模型

所用模型是從原始的3D SEAM Ⅱ FOOTHILL模型中抽取的一個剖面,同時去除了近地表低速區(qū)。真實模型如圖12所示,縱、橫向采樣均是20m,橫向有720個格點,長約14km;縱向有200個格點,深約4km。可以明顯看到地層發(fā)生扭曲和位移,但是模型沿著縱向的整體變化率并不大。采用3Hz的Ricker子波作為震源子波。在頂部布置了48個震源和270個接收器。模型正演采用了吸收邊界條件。圖13為觀測數(shù)據(jù)的炮集。圖中的初至大致沿著一條直線,因為模型縱向變化不大,所以潛水波(diving wave)穿透深度較淺,從而導致走時跟距離接近線性關系。本文創(chuàng)建了4個初始模型,如圖14 所示。前3個由真實模型采用高斯平滑得到,平滑半徑分別為10,20,30個網(wǎng)格長度。最后一個是一維線性遞增模型。相應的合成地震數(shù)據(jù)如圖15所示。初始模型越簡單,合成數(shù)據(jù)也越簡單。對前面兩個初始模型,尚能識別反射波;后面兩個不易辨認反射波。對每個初始模型都進行了6次反演測試,包括1次常規(guī)反演和5次OT反演?？偣灿肧EAM Ⅱ Foothill模型進行了24次反演實驗。

圖12 實際速度模型

圖13 觀測數(shù)據(jù)

圖14 初始模型

圖15 根據(jù)圖14計算出來的合成地震數(shù)據(jù)

第1個初始模型最接近真實模型,且其初始合成數(shù)據(jù)已經(jīng)包含了大量反射事件。圖16顯示了使用第1個初始模型的反演結果。所有6個反演都給出了正確的結構,且反演結果基本相當。常規(guī)全波形反演中的L2目標函數(shù)的吸引域(attraction basin)較小,從另一個角度來說,L2的分辨率較高。圖16的反演結果顯示了OT函數(shù)也同樣具有較高的分辨率。與此相比,若是采用波動方程層析(wave equation tomography),通常需要用L2接著往下反演;而OT則具有一步到位的潛力。

第2個初始模型比第1個更為光滑,大部分清晰的反射結構已經(jīng)被平滑掉。其合成數(shù)據(jù)僅能識別一個大的反射事件。初至波則與第1個模型的初至波較為接近。原因為:①初至波穿透較淺;②初至波沿著淺地表傳播,其到時本身就是模型結構的加權平均,所以對模型光滑不會導致初至太大的變化。圖17 展示的是6個目標函數(shù)對應的反演結果。常規(guī)反演沒有完全重構速度模型,其最深處能識別的結構在中間2km深處。與真實模型相比,其淺部的速度缺乏背景場的更新,更接近于偏移模式。在早期的全波形反演發(fā)展中,地球物理學家已經(jīng)意識到反射波對于中間波長的速度模型是缺乏反應的(intermediate-wavelength gap),所以用L2函數(shù)會導致全波形反演向著偏移的模式演化,并陷入局部極值。圖17所示結果實際上揭示了OT函數(shù)對于走時信息的提取能力。OT函數(shù)是基于模式匹配的,模式上的匹配是隱含著運動學信息(kinematic),如圖17和18所示。一維OT提取的是走時信息,高維OT除了提取走時信息外,還提取了事件沿著空間軸方向的移動信息。這些運動學信息的提取對于模型低波數(shù)成分的恢復起了決定性作用。從圖17的5個OT結果上看,3個二維OT函數(shù)的反演結果略微優(yōu)于兩個一維OT函數(shù)的結果。在模型的左下方,兩個一維OT的結果未能清晰地重構出低速帶。需要說明的是,一維GS-OT的用戶輸入?yún)?shù)不大容易選取最優(yōu)輸入?yún)?shù),在特定模型上累計的經(jīng)驗推廣到另一個完全不同的模型上往往并不會得到同樣的結果。從原理上講,GS-OT是具有良好的凸性的,比W2更優(yōu)越,圖17e中展示的反演結果或許還有進一步提高的可能。

在作者積累的全波形反演經(jīng)驗中,模型陷入局部極值后可能會有一些難以解釋的現(xiàn)象。圖18展示的是用第3個初始模型得到的6個反演結果。相比于圖17中的常規(guī)反演,圖18中的常規(guī)反演得到的模型更新更深。雖然兩者都是因全波形反演的偏移模式而陷入局部極值,但無法解釋為何圖18a會有更深的模型更新。與圖17類似,圖18中同樣是3個二維OT函數(shù)的反演結果略微優(yōu)于兩個一維OT函數(shù)的結果。并且GS-1D的反演結果確實比W2更好。

圖18 以圖14c作為初始模型得到的反演結果

圖19展示的是初始模型為速度隨深度線性增加模型對應的6個反演結果,對于所有這些函數(shù),線性增加的模型可能離真實的模型太遠。W2的反演可能是因為在陷入局部極值的情形下步長不合適而出現(xiàn)了特別大和特別小的速度結構。其它5個反演結果應當是全波形反演的偏移模式,背景模型難以更新。值得注意的是,KR-OT對應的反演結果模型迭代得最深,這可能跟KR-OT獨特的平衡能力相關。KR-OT的輸入信號即是數(shù)據(jù)殘差,在經(jīng)過KR-OT平衡后,來自于深部的反射信號得到了加強,所以導致較深的模型迭代。

圖19 以圖14d作為初始模型得到的反演結果

3 討論

圖11中的常規(guī)全波形反演陷入局部極值,即使沒有真實模型參考也能識別是局部極值。但在實際的全波形反演中,應用常規(guī)全波形反演后遇到的情形并不像圖11那么清晰。關于實際數(shù)據(jù)中出現(xiàn)的局部極值現(xiàn)象,一個更為準確的推斷是,在部分區(qū)域模型恢復得很好,部分區(qū)域與真正的結構差異較大。所以導致了全波形反演與某些測井資料對應得很好,而與其它測井信息則差異較大。也正是因為這個原因,才凸顯出凸性更好的目標函數(shù)的重要性,這樣即可在絕大部分區(qū)域都能成功避免局部極值。從目標函數(shù)的研究歷史看,從最初的點對點的比較,逐步發(fā)展為道與道的比較?；诘辣容^的目標函數(shù)包括:包絡函數(shù)、互相關、反卷積和互相關走時差等一系列函數(shù)。最優(yōu)化傳輸是基于道的整個數(shù)據(jù)體的模式比較,支持整個數(shù)據(jù)體之間的比較。從研究歷史看,基于數(shù)據(jù)整體的比較是一個發(fā)展方向,其優(yōu)勢在于可以在全局比較中正確匹配地震事件的偏差,是全局優(yōu)化策略。

最優(yōu)化傳輸函數(shù)的凸性取決于:①范數(shù)。通常會選擇1或2范數(shù)(p-Wasserstein函數(shù))。通?？梢哉J為2范數(shù)比1范數(shù)更凸。此外,在大多數(shù)情況下,2范數(shù)更容易操作。相比之下,1范數(shù)在零點附近不是平滑的,比如,|x|在零附近的導數(shù)就是間斷的。不恰當?shù)貙?范數(shù)合并到全波形反演中可能會導致伴隨源在某些點或者區(qū)域的不連續(xù)性。KR-OT函數(shù)雖然是基于1范數(shù),但并沒有出現(xiàn)這種不連續(xù)的現(xiàn)象,因為KR-OT是在對偶空間中表示的,1范數(shù)被轉(zhuǎn)移到了約束方程中。②預處理。對于數(shù)據(jù)預處理需要謹慎對待。比如W2常用到的移動加規(guī)則化,其移動量不宜太大?；诜蔷€性變換的預處理有很多,在作者的經(jīng)驗中,tanh表現(xiàn)不錯,而且用戶輸入?yún)?shù)容易調(diào)節(jié)。③時間空間坐標是否引入最優(yōu)化傳輸目標函數(shù),圖空間最優(yōu)化傳輸是將時間坐標和空間坐標顯式地引入,其結果也變得更凸,同時,避免了負信號的處理。

具有全局優(yōu)化潛力的最優(yōu)化傳輸函數(shù)的計算量比常規(guī)函數(shù)要大。在實用化方面,一維的最優(yōu)化傳輸函數(shù)的計算量很小,可以應用在二維和三維全波形反演中。若將三維數(shù)據(jù)分割成二維的矩陣,那么,本文提到的3個二維最優(yōu)化傳輸函數(shù)同樣可以應用,計算代價在可接受范圍內(nèi)。若是直接處理三維數(shù)據(jù)體,計算量較大,可能需要一些加速策略。從目標函數(shù)的性能來看,沒有哪一個函數(shù)可以處理所有的情形,但是有些函數(shù)會在大部分地質(zhì)環(huán)境中表現(xiàn)良好。使用SEAM Ⅱ Foothill模型的全波形反演實驗證明,總體趨勢是2D最優(yōu)化傳輸函數(shù)優(yōu)于1D版本。BFM功能優(yōu)于W2功能;GS-2D技術優(yōu)于GS-1D技術。從理論上說,二維函數(shù)能提取空間軸方向的信息,比一維更具潛力。這與作者在其它各種模型上測試得到的認識是一致的。綜合作者在各種模型上的測試經(jīng)驗看,W2在數(shù)據(jù)模式較為簡單的情形下表現(xiàn)良好,可能是因為預處理影響到了函數(shù)性能。KR-OT因為輸入數(shù)據(jù)是數(shù)據(jù)誤差,導致地震事件從一開始就混雜在一起,所以合成數(shù)據(jù)跟真實數(shù)據(jù)不能相差太遠。但是其具有獨特且優(yōu)越的振幅平衡能力。在陸地數(shù)據(jù)的彈性全波形反演中KR-OT是首選。從原理上說,KR-OT比L2更凸,但會比GS-1D略遜一籌。GS-OT(包括GS-1D和GS-2D)的優(yōu)勢在于將坐標引入計算代價中,從根本上避免了數(shù)據(jù)的預處理。但是,凸性良好的函數(shù)具有共同的小瑕疵,即在反演模型接近真實模型時,反映在數(shù)據(jù)上的偏差可能不會太大,所以反演收斂速度會很慢。基于此,在GS-OT反演后,可以酌情在后續(xù)加入常規(guī)反演。BFM的性能原則上比KR-OT更凸,但遜于GS-2D。BFM是基于2范數(shù)的最優(yōu)化傳輸,而KR-OT是基于1范數(shù)的,2范數(shù)比1范數(shù)更凸。GS-2D也是基于2范數(shù),同時,避免了數(shù)據(jù)預處理,應當比BFM更凸。

4 結論

本文綜合對比了目前學術界和工業(yè)界新提出的最優(yōu)化傳輸目標函數(shù)。常規(guī)的最小二乘目標函數(shù)是點對點的比較方式,對于勘探地震數(shù)據(jù)地震道中的時間相關性和地震道之間的空間相關性沒有顯式地提取。當初始模型足夠準確,合成地震事件與真實地震事件偏差在半個波長以內(nèi)的時候,常規(guī)的最小二乘目標函數(shù)是可以反演出足夠接近真實的地下速度模型。然而在實際數(shù)據(jù)處理中,相對于采集中設定的低頻范圍,初始模型常常是不夠準確的。這也是學術界和工業(yè)界尋求其它目標函數(shù)的動力。從目標函數(shù)的研究歷史看,的確是朝著高維方向發(fā)展。在目前的研究中,最優(yōu)化傳輸通常采用1范數(shù)和2范數(shù);對于給定的范數(shù),又可以分為單道處理和高維處理。綜合我們的測試結果看,高維的最優(yōu)化傳輸函數(shù)優(yōu)于一維的版本。在實際應用中,可以綜合考慮計算成本和效果。對于工區(qū)不太復雜的情形,可以采用W2,因為其存在解析解,計算速度最快。在復雜地區(qū),可以采用BFM,KR-OT,GS-1D或者GS-2D;在這4種最優(yōu)化傳輸函數(shù)中,以GS-1D計算成本最小。其中,只有KR-OT是基于1范數(shù),其凸性可能遜于BFM;但是KR-OT獨有的振幅平衡能力使得KR-OT在處理存在強振幅差異地震數(shù)據(jù)的時候具有特殊的優(yōu)勢。