全波形反演的一個新目標函數(shù):數(shù)據(jù)域中的微分相似優(yōu)化

高福春,Paul WILLIAMSON,R. Gerhard PRATT

(1.道達爾勘探和生產(chǎn)研究技術中心,休士頓77002,美國;2.西安大略大學,倫敦N6A 3K7,加拿大)

全波形反演的一個新目標函數(shù):數(shù)據(jù)域中的微分相似優(yōu)化

高福春1,Paul WILLIAMSON1,R. Gerhard PRATT2

(1.道達爾勘探和生產(chǎn)研究技術中心,休士頓77002,美國;2.西安大略大學,倫敦N6A 3K7,加拿大)

全波形反演(full waveform inversion,FWI)目前已有廣泛的工業(yè)實踐,但因其本質(zhì)上的非線性,不如走時層析成像等傳統(tǒng)速度建模技術穩(wěn)健,非線性程度也因目標函數(shù)不同而不同。研究分析了FWI中幾種不同目標函數(shù)的性質(zhì),基于定義在數(shù)據(jù)域中的微分相似概念,提出了一種新的目標函數(shù)。初步試驗表明,這種目標函數(shù)對于比較大范圍的數(shù)據(jù)殘差都有凸狀性質(zhì),基于梯度優(yōu)化法時使用該目標函數(shù)的FWI比傳統(tǒng)FWI更穩(wěn)健,而且波形反演的良好分辨率基本得以保留。

全波形反演;目標函數(shù);非線性;微分相似優(yōu)化

三維全波形反演(FWI)在20年前還被認為不可實現(xiàn),但是目前已經(jīng)在工業(yè)上得到應用,這主要得益于計算機硬件和數(shù)值算法的持續(xù)發(fā)展[1-2]。標準的FWI可以得到比反射層析成像等常規(guī)速度建模方法分辨率更高的速度模型,但因其高度非線性,不如常規(guī)速度建模方法穩(wěn)健。為了確保FWI實現(xiàn)物理意義上的收斂,基于二范數(shù)目標函數(shù)的常規(guī)FWI需要初始模型足夠接近真實模型,或者數(shù)據(jù)中含有可靠的低頻信息[3]。通常,人們通過折射波走時層析成像(如PRATT[4],GAO等[5])或深度域成像建模流程(包括反射波層析成像,如SIRGUE等[6])來得到接近真實模型的初始模型,然而,當介質(zhì)比較復雜時,這些方法不一定能產(chǎn)生一個足夠接近真實模型的初始模型,從而確保FWI在使用典型地震數(shù)據(jù)帶寬時總是獲得成功。另一方面,不同類型的目標函數(shù)對相同模型的擾動以及數(shù)據(jù)中噪聲的敏感度不同。為替代傳統(tǒng)的二范數(shù)目標函數(shù),人們提出多種新的目標函數(shù),主要有以下幾種:①人工合成數(shù)據(jù)和實際觀測數(shù)據(jù)之間基于其它范數(shù)的差異性(如TARANTOLA[7],PRATT等[8])的衡量。HUGUES等[9]使用一范數(shù)目標函數(shù)反演了一個墨西哥灣的二維數(shù)據(jù)。GUITTON等[10]推導出基于休伯范數(shù)的反演算法,并證明該算法對數(shù)據(jù)中的異常值不太敏感。KIM等[11]設計出一種FWI算法,旨在頻率域里減少人工合成數(shù)據(jù)和觀測數(shù)據(jù)之間基于二范數(shù)的相位差,而不是振幅和相位的整體差。ARAVKIN等[12]推導出基于學生-T范數(shù)(休伯范數(shù)的一個變種)的FWI算法,該算法對集中在某一區(qū)域的異常值不很敏感。②人工合成數(shù)據(jù)和實際觀測數(shù)據(jù)相似性的度量。ROUTH等[13]設計了基于相關性目標函數(shù)的FWI算法,通過最大限度地提高計算數(shù)據(jù)和觀測數(shù)據(jù)之間的相似性來更新速度模型,而不是最小化它們之間的差別。這兩類替代目標函數(shù)的提出大都是為了達到某一個特定的目標。③人工合成數(shù)據(jù)和觀測數(shù)據(jù)之間經(jīng)過相同變換或處理以后的差異衡量。HA等[14]推導出拉普拉斯域中的FWI算法,指出該方法通過最小化拉普拉斯域中的差異,使其比傳統(tǒng)FWI對數(shù)據(jù)中低頻成分的依賴性相對減弱。LUO等[15]推導并實現(xiàn)了地震波波形包絡(而不是波形本身)反演方法,可重建速度模型的大尺度成分。BAEK等[16]推導出地震波事件記錄引導下的最小二乘反演方法,以減輕傳統(tǒng)FWI在全域最小附近收斂半徑較窄的問題。這類方法的共性是先進行數(shù)據(jù)變換或處理使得信號在其它域中具有相對低頻成分,然后再比較人工合成數(shù)據(jù)和觀測數(shù)據(jù)之間的差異使其發(fā)生周跳的可能性降低,減弱對初始模型的敏感性。上述方法都有其各自的缺點。定義在像域中的微分相似優(yōu)化(differential semblance optimization,簡稱DSO)概念最初由SYMES等[17]提出,作為一種在像域中為更新速度模型而建立的目標函數(shù),它具有凸狀特性優(yōu)勢,受到格外重視。PRATT等[18]分析了一系列基于相似性的目標函數(shù),初步證明DSO概念可以擴展到數(shù)據(jù)域中,并指出這些函數(shù)可作為替代目標函數(shù),使FWI更加穩(wěn)健。雖然PRATT等指出了這些目標函數(shù)對于模型的擾動比常規(guī)目標函數(shù)更線性,然而他們沒有實現(xiàn)用這些函數(shù)作為目標函數(shù)的FWI算法,所以這些函數(shù)在理論上擁有的優(yōu)點還有待進一步證實。本文提出了與PRATT等相似但更加線性化的目標函數(shù),分析并選擇了其中一種作為FWI的目標函數(shù)并予以實現(xiàn)。該FWI算法更加穩(wěn)健,分辨率較傳統(tǒng)的全波形反演有所降低,但較常規(guī)射線類反演方法有明顯提高。將這種目標函數(shù)與傳統(tǒng)的目標函數(shù)聯(lián)合使用能得到分辨率較高的模型。

1 方法原理

1.1 不同的目標函數(shù)

傳統(tǒng)FWI定義在頻率和空間域中基于二范數(shù)的目標函數(shù)(如PRATT等[8])為:

(1)

式中:u和d分別是人工正演模擬和觀察到的波場;“*”表示復共軛。PRATT等[18]提出了如下相干性的衡量:

其中,s被定義為震源因子:

(4)

式中:“T”表示轉置。公式(4)中的震源因子可用m個地震道集來計算,1≤m≤N,N為共震源道集中所有的地震道數(shù)量。公式(3)中的導數(shù)可以是對震源空間位置或檢波器空間位置的導數(shù)。

在本次研究工作中,我們提出了另外一種新的目標函數(shù),稱之為窗口式相似性目標函數(shù)(簡稱WS目標函數(shù))。該函數(shù)在原理上和公式(3)中描述的目標函數(shù)沒有本質(zhì)的區(qū)別,但是更便于實現(xiàn),在實際應用中更具有靈活性。

為了比較由方程(1)—(3)所定義的目標函數(shù)和WS目標函數(shù)對于模型擾動的不同敏感程度,利用PRATT等[18]使用的二維人工合成模型,采用井間數(shù)據(jù)采集方式對4種目標函數(shù)進行了測試。兩口井水平距離150m,井深300m,分別布置51個震源和51個檢波器。模型的速度只是深度的函數(shù),即一維結構模型。其速度分布是一個常背景速度加上一恒定的梯度(v=1000+5z),所以速度隨深度增加,從1000m/s上升到底部的2000m/s,見圖1a。模型中使用的梯度擾動了20次,每次改變5%,從0一直變化到正確梯度的95%,便得到了20個不同的擾動模型,參見圖1b。

圖2顯示了4種不同目標函數(shù)隨速度擾動的演化過程，可以看到：除了常規(guī)的二范數(shù)目標函數(shù),其它3種目標函數(shù)的變化中都不存在區(qū)域最小。常規(guī)目標函數(shù)的變化中存在區(qū)域最小,而且不止一個,如圖2a 所示。采用梯度法時,只有當初始模型的梯度值是正確值的85%時,常規(guī)FWI才會收斂到全域最小。從圖2b至圖2d中3條曲線可以看到,如果它們被當作目標函數(shù),即使速度模型離正確值很遠,也可以收斂到全域最小。但這3個函數(shù)的演化過程不同:圖2b中由公式(2)定義的函數(shù)變化一直緩慢,直到離正確模型很近時才加快收斂,其它兩個函數(shù)的收斂則比較均勻。由于WS函數(shù)相對比較容易實施,我們這里選擇WS函數(shù)作為目標函數(shù)。

圖1 速度模型(a)和各種擾動后的速度隨深度的線性變化(b)

圖2 四種目標函數(shù)收斂性隨速度擾動的變化a 公式(1); b 公式(2); c 公式(3); d WS目標函數(shù)

1.2 使用WS目標函數(shù)的FWI

與主流FWI的實施過程一樣,我們也是用迭代法使得目標函數(shù)變小來更新模型,即:

(5)

2 應用實例

為了驗證WS目標函數(shù)的有效性,我們用Marmousi模型數(shù)據(jù)進行了一系列測試。該數(shù)據(jù)為典型的淺海采集數(shù)據(jù),具有復雜的沉積構造和很強的橫向不均勻性(圖3a)。在2.0km深、9.2km寬的區(qū)域內(nèi),速度變化范圍為1500～5500m/s,共有93個共炮點道集,檢波器在地表均勻分布,位置固定不變,間距為50m,共185個,炮點間距100m。

首先使用一個較好的初始模型進行了一系列測試。該初始模型為真實模型的平滑結果,見圖3b。依次反演了3個頻率分量,分別為3.0,3.5和4.0Hz。

圖3 Marmousi模型數(shù)據(jù)測試結果a 真實模型; b 初始模型; c 常規(guī)FWI反演結果; d 基于WS目標函數(shù)(窗口寬度3)的FWI反演結果; e 基于WS目標函數(shù)(窗口寬度21)的FWI反演結果

采用WS目標函數(shù)時,我們選擇了兩個不同的窗口寬度,分別為3和21。圖3c為常規(guī)FWI反演的最終速度模型,圖3d和圖3e均是用WS目標函數(shù)FWI反演的最終結果,前者窗口寬度為3,后者窗口寬度為21?？梢钥吹?用WS目標函數(shù)時,不管窗口長度如何,在反演頻率允許的分辨率范圍內(nèi),都可以成功地重建真實模型,但它們的分辨率,尤其是采用的窗口寬度比較大時得到的反演模型,要稍低于傳統(tǒng)FWI得到的速度模型。從沉積層和模型邊緣尤其能看出分辨率的不同。沉積層分辨率的差異可歸因于目標函數(shù)的不同,邊緣分辨率的差異可歸因于窗口寬度的不同。由上述分析可以看出,WS目標函數(shù)對模型的擾動敏感性降低了,所以FWI反演的結果分辨率會低一些。窗口寬度越大,邊緣效應就越強,所以模型的邊緣部分受到的影響較大。

以上測試只是確定了WS函數(shù)可以作為FWI的目標函數(shù),但是還沒有證實它的收斂半徑較大。為了進一步測試該目標函數(shù)的性質(zhì),我們對不同目標函數(shù)做了一系列測試(圖4,圖5)。采用的兩個初始模型(模型1和模型2)都是一維結構模型,即速度只是深度的函數(shù),離真實模型比較遠。這兩個模型雖然都是一維模型,但是差異比較大。模型1速度從地表的1500m/s以常梯度變化到底部的2625m/s,而模型2 則從地表1500m/s變化到底部的2000m/s。顯而易見,模型2偏離真實模型更遠,幾乎是常速度模型。用模型1(圖4b)反演了2～6Hz的數(shù)據(jù),常規(guī)FWI得到的最終模型(圖4c)顯然陷入了區(qū)域最小,而用WS目標函數(shù)得到的最終模型(圖4d)則比較合理地重建了真實模型。將圖4d對應的速度模型作為初始模型,我們可以通過常規(guī)FWI得到分辨率較高的速度模型(圖4e)。

圖4 初始模型用模型1的測試結果a 真實模型; b 初始模型; c 常規(guī)FWI反演的最終模型; d 用WS目標函數(shù)FWI得到的最終模型; e 常規(guī)FWI用WS反演模型作為初始模型得到的最終模型(2～6Hz)

圖5 初始模型用模型2的測試結果a 初始模型; b 常規(guī)FWI反演的最終模型; c 用WS目標函數(shù)FWI得到的最終模型; d 常規(guī)FWI用WS反演模型作為初始模型得到的最終模型(2～6Hz)

圖5為用模型2測試的結果。因為模型2離真實模型更遠,可以看到常規(guī)FWI結果誤差很大(圖5b),用WS目標函數(shù)的FWI則可以在一定區(qū)域內(nèi)重建真實模型(圖5c)。這里是指左右下角之外的區(qū)域,原因在于這兩個區(qū)域受到地震波照明的限制。用WS目標函數(shù)的FWI結果作為初始模型,我們又能通過常規(guī)FWI得到分辨率比較高的模型(圖5d)。

3 討論和結論

本文展示并實現(xiàn)了一個FWI常規(guī)目標函數(shù)的替代函數(shù),即WS目標函數(shù)。它的提出受到了PRATT等于2002年提出的幾個目標函數(shù)的啟發(fā)。WS目標函數(shù)和PRATT等提出的目標函數(shù)有相似的性質(zhì),但相對比較容易實現(xiàn),而且可以通過窗口的寬度調(diào)節(jié)目標函數(shù)的性質(zhì)。Marmousi模型測試證實WS目標函數(shù)具有比較大的收斂半徑,較少的局部最小值,能使FWI變得穩(wěn)健。下一步研究工作是對該算法作進一步測試,并將其用于實際地震數(shù)據(jù)。

WS目標函數(shù)會使FWI的分辨率稍微變低,這可以理解而且不是問題,因為該目標函數(shù)可以和常規(guī)目標函數(shù)結合起來使用,使得FWI既穩(wěn)健又具有很好的分辨率。

致謝:感謝道達爾勘探研究和生產(chǎn)技術中心允許發(fā)表這項研究成果,感謝Laurent Sirgue為本研究中的測試提供了所需要的初始模型。

[1] SIRGUE L,ETGEN J,ALBERTIN U.3D full waveform inversion:wide-versus narrow-azimuth acquisitions[J].Expanded Abstracts of 77thAnnual Internat SEG Mtg,2007:1760-1763

[2] PLESSIX R E,GUIDO B,DE MAAG J W,et al.Application of acoustic full waveform inversion to a low-frequency large-offset land data set[J].Expanded Abstracts of 80thAnnual Internat SEG Mtg,2010:930-934

[3] PRATT G R,GAO F,ZEL T C,et al.A comparison of ray-based and waveform tomography:implications for migration[R].Florence,Italy:64thEAGE Conference and Exhibition,2002

[4] PRATT G R.Velocity models from frequency-domain waveform tomography:past,present and future[R].Paris,France:66thEAGE Conference and Exhibition,2004

[5] GAO F C,LEVANDER A R,PRATT R G,et al.Waveform tomography at a ground water contamination site:VSP-surface dataset[J].Geophysics,2006,71(1):H1-H11

[6] SIRGUE L,DENEL B,GAO F.Integrating 3D full waveform inversion into depth imaging projects[J].Expanded Abstracts of 81stAnnual Internat SEG Mtg,2011:2354-2358

[7] TARANTOLA A.Inversion of seismic reflection data in the acoustic approximation[J].Geophysics,1984,49(8):1259-1266

[8] PRATT G R,SHIN C,HICKS G J.Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion[J].Geophysical Journal International,1998,133(2):341-362

[9] HUGUES A D,TARANTOLA A.Multiparameter l1 norm waveform fitting:interpretation of Gulf of Mexico reflection seismograms[J].Geophysics,1999,64(4):1023-1035

[10] GUITTON A,SYMES W W.Robust inversion of seismic data using the Huber norm[J].Geophysics,2003,68(4):1310-1319

[11] KIM W,SHIN C.Phase inversion of seismic data with unknown source wavelet:synthetic examples[J].Expanded Abstracts of 75thAnnual Internat SEG Mtg,2005:1685-1689

[12] ARAVKIN A,VAN LEEUWEN T,HERMANN F.Robust full-waveform inversion using the student’s t-distribution[J].Expanded Abstracts of 81stAnnual Internat SEG Mtg,2011:2669-2673

[13] ROUTH P,KREBS J,LAZARATOS S,et al.Encoded simultaneous source full-wavefield inversion for spectrally shaped marine streamer data[J].Expanded Abstracts of 81stAnnual Internat SEG Mtg,2011:2433-2437

[14] HA W,SHIN C.Laplace-domain full-waveform inversion of seismic data lacking low frequency information[J].Geophysics,2012,77(5):R199-R206

[15] LUO J,WU R.Envelope inversion—some application issues[J].Expanded Abstracts of 83rdAnnual Internat SEG Mtg,2013:1042-1047

[16] BAEK H,CALANDRA H,DEMANET L.Velocity estimation via registration-guided least-squares inversion[J].Geophysics,2014,79(2):R79-R89

[17] SYMES W W,CARAZZONE J J.Velocity inversion by differential semblance optimization[J].Geophysics,1991,56(5):654-663

[18] PRATT G R,SYMES W W.Semblance and differential semblance optimization for waveform tomography:a frequency domain implementation[J].Journal of Conference Abstracts,2002,7(2):183-184

(編輯：戴春秋)

A new objective function for full waveform inversion:differential semblance optimization in data domain

GAO Fuchun1,Paul WILLIAMSON1,R. Gerhard PRATT2

(1.TotalE&PResearchandTechnology,Houston77002,USA;2.UniversityofWesternOntario,LondonN6A3K7,Canada)

Full Waveform Inversion (FWI),while now widely practiced industrially,is less robust than many conventional velocity model building techniques,such as travel time tomography,due to its high non-linearity.Different objective functions in FWI have different degrees of non-linearity.In this study,we investigate the behavior of FWI with different objective functions and propose a new objective function based on differential semblance defined in the data domain.Preliminary tests suggest that this objective function is convex for a large range of data residuals.Gradient-based optimization schemes are therefore more robust than for the standard least-squares formulation;however,the good resolving power of waveform inversion is mostly retained.

full waveform inversion,objective function,non-linearity,differential semblance optimization

2016-10-08;改回日期：2016-11-18。

高福春(1968—),男,道達爾勘探和生產(chǎn)研究技術中心地球物理研究員,研究方向為反演及其在地震波成像中的應用、地震波偏移等。

P631

A

1000-1441(2017)01-0026-05

10.3969/j.issn.1000-1441.2017.01.003