李 翔

(加拿大英屬哥倫比亞大學(xué),溫哥華V6T1Z4)

基于壓縮感知技術(shù)的全波形反演

李 翔

(加拿大英屬哥倫比亞大學(xué),溫哥華V6T1Z4)

全波形反演技術(shù)雖然已經(jīng)得到了成功應(yīng)用，但其求解一個最小二乘非凸優(yōu)化問題的計算量仍是一個很大的難題。通過隨機降采樣技術(shù)可以減少反演過程中炮數(shù)和頻率數(shù)，從而可以極大程度地降低全波形反演的計算量；然而這種方法受到奈奎斯特采樣定律例證的“維數(shù)災(zāi)難”的限制以及背離‘摩爾定律’現(xiàn)象。為此,研究了基于改進壓縮感知的隨機化降維技術(shù),應(yīng)用壓縮感知理論減少隨機采樣;聯(lián)合隨機采樣和稀疏促進技術(shù),成功減少了地震數(shù)據(jù)的維數(shù),同時保持了有效信息。通過該項技術(shù)的應(yīng)用,牛頓類方法的計算量相當(dāng)于全波場采樣梯度類算法的計算量;將稀疏約束應(yīng)用在反演過程中的模型更新上,不改變波形反演的目標函數(shù),并且能夠壓制由欠采樣產(chǎn)生的虛像噪聲。北海模型數(shù)據(jù)測試結(jié)果證明了該方法的可行性和有效性。

壓縮感知;波形反演;曲波變換;稀疏促進

全波形反演(FWI)[1-5]利用地表或井中觀測的地震信號進行反演以獲得地下介質(zhì)的物理參數(shù),如速度或密度分布等。數(shù)學(xué)上,全波形反演可以看作是一種基于偏微分方程約束的優(yōu)化問題,通過不斷迭代更新模型參數(shù)來減小觀測數(shù)據(jù)和模型數(shù)據(jù)之間的誤差(又稱目標函數(shù),通常用最小二乘誤差表示),當(dāng)目標函數(shù)小于預(yù)定誤差時,就認為所獲得的模型參數(shù)能夠表征地下介質(zhì)的物理參數(shù)分布[3,6]。最近20年,大多是研究一階基于梯度法的優(yōu)化算法在全波形反演問題中的應(yīng)用,如梯度法[3,7]和非線性共軛梯度法[7-11]。然而正如PRATT等[3]和SHIN等[4]的文獻所述,一階方法難以得到可靠的步長,導(dǎo)致一階方法應(yīng)用于FWI問題時收斂較慢,從而影響反演結(jié)果的分辨率。

根據(jù)優(yōu)化理論[3,7,9],當(dāng)初始模型比較接近真實模型時,二階優(yōu)化算法相對于一階算法具有更好的收斂速度,而且二階算法求解的是包含二階信息的海森(Hessian)矩陣的逆,并將其作用在模型更新項,可以補償與震源相關(guān)的模糊效應(yīng)、孔徑不足效應(yīng)和其它與振幅相關(guān)的影響。對于大型反演問題,因為真實海森矩陣不是正定的,顯式構(gòu)建真實的海森矩陣非常困難,求解真實海森矩陣的逆也相當(dāng)困難,故通常利用高斯-牛頓海森矩陣來近似真實的海森矩陣,而高斯-牛頓海森矩陣的逆通常是對稱正定(半正定)的,可以通過迭代法求解[10,12]。因為每次計算海森矩陣作用在某個向量上需要求解多個波動方程,所以求解波形反演中高斯-牛頓海森矩陣逆的計算量非常大。

為了解決求解高斯-牛頓海森矩陣逆計算量巨大的問題,本文采用將模型的更新轉(zhuǎn)入曲波域[13-14]并進行稀疏約束的辦法,對數(shù)據(jù)進行降維采樣,隨機選取炮點及頻率數(shù)[15-16],以減小計算量和提高反演結(jié)果的分辨率。最后用北海模型數(shù)據(jù)測試了方法的可行性和有效性。

1 方法原理

1.1 高斯-牛頓法全波形反演

全波形反演可表示為:

(1)

式中:P為輸入的實際觀測地震數(shù)據(jù),非線性函數(shù)F[m,Q]=DH-1[m]Q為頻率域波動方程正演模擬的地震數(shù)據(jù)。如(1)式所示,求解頻率域波動方程需要對每一頻率求解大型稀疏亥姆霍茲(Helmholtz)矩陣H的逆,其中m為模型參數(shù),Q為震源矩陣。

利用標準高斯-牛頓法實現(xiàn)公式(1)需要求解高斯牛頓的子問題(算法1,具體實現(xiàn)程序見圖1),每一步迭代需要計算傳統(tǒng)波形反演的梯度及海森矩陣的逆,然后將海森矩陣的逆作用到梯度上,最后通過線性搜索選擇合適的步長。

圖1 高斯-牛頓子問題(算法1)的實現(xiàn)程序

然而利用高斯-牛頓法求解上述波形反演問題時計算量非常大,因為每一步迭代都需要求解一個線性的高斯-牛頓子問題(即求解海森矩陣的逆)。而求解子問題也需要多次進行子迭代(此迭代不同于波形反演模型迭代),每次子問題的子迭代需要多次求解波動方程[17-18]。

為了解決計算量巨大的問題,我們結(jié)合降維技術(shù)和壓縮感知技術(shù),隨機降低反演的炮集數(shù)據(jù)和隨機減小頻率數(shù)來減小數(shù)據(jù)量。目前有2種降低炮數(shù)據(jù)的方法,第一種是同時激發(fā)震源技術(shù),該技術(shù)將傳統(tǒng)的單炮震源隨機加權(quán)相加,形成一個超級炮集(一炮有多個震源同時激發(fā));第二種是隨機抽取一部分子集。利用上述降維方法,計算一個高斯-牛頓速度更新雖然需要多次子迭代,但是其計算量可以控制在近似等于一個傳統(tǒng)波形反演梯度迭代的計算量。頻率域超級炮的表達式為:

(2)

式中:FΩ為傅里葉變換。將上述降維算子應(yīng)用于全波形反演,波形反演問題(公式(1))就變?yōu)?

(3)

式中,下劃線代表壓縮降維后的數(shù)據(jù)和算子,即[15,19]:

(4)

數(shù)據(jù)壓縮降維后,波形反演優(yōu)化問題的高斯-牛頓迭代僅需要利用一部分頻率和炮集數(shù)據(jù)。本文方法極大地減少了反演中求解波動方程的數(shù)量。波形反演中需要采用目標函數(shù)的梯度,其梯度表達式為:

(5)

式中:u,v為震源正傳波場和接收點反傳波場。實際上(5)式即為逆時偏移算子表達式,逆時偏移算子的轉(zhuǎn)置即為波動方程的線性近似算子,又稱為Born近似算子。

1.2 稀疏反演迭代

波形反演梯度法利用梯度來更新模型,但梯度并不等于速度更新量,所以需要精確的步長,在實施中尋找精確的步長非常困難[3,5]。二階高斯-牛頓法用海森矩陣的逆來標度梯度,使其擁有相對正確的振幅信息。采用壓縮降維技術(shù)后,波形反演的計算量可以得到極大的減小,然而數(shù)據(jù)量的減小會帶來壓縮產(chǎn)生的人工虛像。所以本文利用壓縮感知技術(shù)求解高斯-牛頓方法的子問題,以此代替?zhèn)鹘y(tǒng)共軛梯度法求解。在本文提出的方法中,我們直接求解Born正演算子(優(yōu)化理論中常稱其為雅可比矩陣)的逆。

對于大多數(shù)反演問題,人們往往避免采用迭代法求解雅可比矩陣的逆,因為迭代法在每次迭代時需要不停地計算其正傳及其轉(zhuǎn)置,因而計算量巨大。此外,在很多優(yōu)化反演問題中,雅可比矩陣往往是病態(tài)的(海森矩陣是奇異的)。然而與很多反演問題不同,地震勘探中的海森矩陣是相對均衡的。因為在初始模型相對理想的情況下,線性Born近似正演能夠很好地匹配非線性波動方程正演,在此情況下海森矩陣的逆可以用簡單的對角矩陣近似代替[18]。

本文中,我們利用壓縮感知技術(shù)降低求解海森矩陣逆的計算量。同時利用L1范數(shù)正則化線性反演來計算“類似高斯牛頓”的迭代更新項,如:

(6)

式中:A∶=RMJS*為Born近似正演矩陣(亦稱雅可比算子,A′A為該問題的海森矩陣);S*為反曲波變換。該稀疏反演技術(shù)在曲波域?qū)ふ乙粋€最稀疏的能夠匹配觀測數(shù)據(jù)的解。曲波變換表征地下層狀介質(zhì)非常有效,所以利用公式(6)求得的最稀疏解可以更容易反演得到接近地下構(gòu)造的模型。降維技術(shù)結(jié)合稀疏反演既能解決二階高斯-牛頓法計算量大的問題,又能壓制由于數(shù)據(jù)壓縮導(dǎo)致的虛像影響。

波形反演存在局部極小點的問題,即所謂的“周期跳躍”現(xiàn)象。人們提出了很多解決這一難題的方法,比如在反演時引入高頻信息和走時信息等[3,5]。

本文中,我們利用一個新的策略讓模型迭代在曲波域中的L1范數(shù)逐漸地變大。這樣讓最接近真實地層構(gòu)造信息先進入數(shù)據(jù)模型結(jié)果,此方法稱為算法2,具體實現(xiàn)程序見圖2。

圖2 算法2的實現(xiàn)程序

綜上所述,通過降維處理,我們可以極大地減少在全波形反演中所需要求解波動方程(偏微分方程)的次數(shù),從而極大地減少計算量,使求解二階高斯-牛頓模型更新的計算量近似等于求一個梯度更新[17,20]。本文方法利用了全波形反演模型更新中結(jié)構(gòu)信息在曲波域是稀疏的,而降維采樣產(chǎn)生的人為虛像在曲波域不稀疏的性質(zhì),在理論上對高斯-牛頓更新進行稀疏約束。在反演迭代過程中,不相干的信息完全被曲波域稀疏約束壓制,該方法也可以用于最小二乘偏移[17.20]。本文給出的稀疏約束全波形反演技術(shù)不僅能夠極大地減小計算量,也能提高反演結(jié)果的分辨率[2,17]。該技術(shù)具有以下優(yōu)點:

1) 高斯-牛頓法的子問題是線性的,而利用稀疏約束的方法求解線性問題正好符合壓縮感知技術(shù)的基本理論。根據(jù)壓縮感知理論[13,21],即壓縮的信號可以顯著低于尼奎斯特采樣率下通過求解稀疏約束的方法恢復(fù)原來的信號。所以壓縮感知技術(shù)結(jié)合隨機選炮以及同時激發(fā)震源技術(shù)可以有效地減小計算量及去除降維采樣帶來的人工虛像。

2) 將曲波域的稀疏約束應(yīng)用于全波形反演問題的高斯-牛頓更新項上,并不會改變?nèi)ㄐ畏囱輧?yōu)化問題的目標函數(shù),所有求解的過程也比較簡單。而且高斯-牛頓法的子問題是線性的,該問題是凸問題并且有全局極小點[22]。

3) 曲波變換表征地下層狀結(jié)構(gòu)非常有效[13-14],所以波形反演以及地震成像中的層狀同相軸形結(jié)構(gòu)在曲波域可以用少量的稀疏表達,該技術(shù)也可用于最小二乘偏移[20-21]。

2 北海模型數(shù)據(jù)測試

為了評估波形反演在解決細節(jié)結(jié)構(gòu)中的能力,以及驗證本文提出的基于壓縮感知的稀疏約束高斯-牛頓法的有效性,我們采用BG北海速度模型數(shù)據(jù)進行了測試。

如圖3所示,該速度模型包含了由實際地震、測井、地質(zhì)等數(shù)據(jù)得到的地質(zhì)構(gòu)造信息。觀測數(shù)據(jù)由時間域有限差分正演模擬產(chǎn)生,采用15Hz雷克子波。數(shù)據(jù)共350炮,炮間距20m,每炮700道,檢波距10m。

正演結(jié)果通過求解亥姆霍茲矩陣得到。圖4為初始速度模型,可看出,速度在橫向上沒有變化。

反演利用頻率域方法從5Hz開始,為了避免局部極小點,共反演了8個頻帶,跨度從5Hz到15Hz,每個頻帶選取20炮以及隨機選取3個頻率(每個頻帶總共10個頻率),進行稀疏約束下高斯-牛頓模型迭代計算(共10次)。對于每個高斯-牛頓法的子問題,利用L1范數(shù)約束下譜梯度投影(SPGL1)[22-24]線性算法進行子迭代(10次)。圖5為傳統(tǒng)高斯-牛頓算法反演結(jié)果,圖6為稀疏約束下高斯-牛頓算法反演結(jié)果。對比圖5和圖6可以明顯看出,如果沒有稀疏約束,在采樣率嚴重不足的情況下,壓縮數(shù)據(jù)會引入很多虛像噪聲,嚴重影響了反演結(jié)果;利用稀疏約束之后,由降維采樣帶來的虛像明顯減少,反演結(jié)果更加真實可靠。圖7和圖8是數(shù)據(jù)匹配的結(jié)果,背景為真實的單炮記錄,其中圖7中的波形圖(綠色波形曲線)為初始模型的正演結(jié)果,圖8為全波形反演模型的正演結(jié)果(綠色波形曲線),可見,反演結(jié)果實現(xiàn)了更好的數(shù)據(jù)匹配。

圖3 真實的BG速度模型

圖4 初始速度模型

圖5 傳統(tǒng)高斯-牛頓法波形反演結(jié)果(每次迭代隨機利用1.7%的數(shù)據(jù))

圖6 稀疏約束高斯-牛頓法反演結(jié)果(每次迭代隨機利用1.7%的數(shù)據(jù))

圖7 初始模型正演紀錄(背景為真實數(shù)據(jù))

圖8 全波形反演結(jié)果正演紀錄(背景為真實數(shù)據(jù))

3 結(jié)論和展望

方法研究和應(yīng)用實例表明:

1) 利用壓縮降維技術(shù),可以極大地減小在反演過程中所用的數(shù)據(jù)量,從而顯著地減小計算量。

2) 將稀疏約束應(yīng)用于反演中的模型更新,不會改變波形反演的目標函數(shù),能夠壓制由欠采樣產(chǎn)生的虛像噪聲。

3) 曲波變換是個多尺度、多角度的可逆變換,在表征地質(zhì)模型上非常有效。所以將反演模型的迭代變換到曲波域,更有利于稀疏約束算法。

壓縮感知技術(shù)打破了傳統(tǒng)尼奎斯特采樣定理,本文將其中的相關(guān)原理應(yīng)用于全波形反演算法,取得了較好的效果。該方法不僅可以用于地震波形反演,也可用于偏移成像和地震數(shù)據(jù)采集。將壓縮感知技術(shù)運用到數(shù)據(jù)采集,可以極大地降低成本,有了壓縮采集的數(shù)據(jù),可直接利用本文方法進行數(shù)據(jù)處理,將會更大地節(jié)省地震勘探成本。

[1] ASKAN A,AKCELIK V,BIELAK J,et al.Full waveform inversion for seismic velocity and anelastic losses in heterogeneous structures[J].Bulletin of the Seismological Society of America,2007,97(6):1990-2008

[2] HERRMANN F J,LI X,ARAVKIN A Y,et al.A modified,sparsity promoting,Gauss-Newton algorithm for seismic waveform inversion[J].Proceedings of SPIE 8138,2011:1-20

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

[4] SHIN C,JANG S,MIN D J.Improved amplitude preservation for prestack depth migration by inverse scattering theory[J].Geophysical Prospecting,2001,49(5):592-606

[5] VIRIEUX J,OPERTO S.An overview of full-waveform inversion in exploration geophysics[J].Geophysics,2009,74(6):WCC1-WCC26

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

[7] MéTIVIER L,BROSSIER R,VIRIEUX J,et al.Full waveform inversion and the truncated Newton method[J].SIAM Journal on Scientific Computing,2013,35(2):B401-B437

[8] GILBERT J C,NOCEDAL J.Global convergence properties of conjugate gradient methods for optimization[J].SIAM Journal on Optimization,1992,2(1):21-42

[9] GRATTON S,LAWLESS A S,NICHOLS N K.Approximate gauss-newton methods for nonlinear least squares problems[J].SIAM Journal on Optimization,2007,18(1):106-132

[10] HESTENES M R,STIEFEL E.Methods of conjugate gradients for solving linear systems[J].Journal of Research of the National Bureau of Standards,1952,49(6):409-436

[11] WARNER M, RATCLIFFE A,NANGOO T,et al.Anisotropic 3d full-waveform inversion[J].Geophysics,2013,78(2):R59-R80

[12] BAUSCHKE H H,COMBETTES P L,LUKE D R.Phase retrieval,error reduction algorithm,and fienup variants:a view from convex optimization[J].Journal of the Optical Society of America A,2002,19(7):1334-1345

[14] KUMAR V.Incoherent noise suppression and deconvolution using curvelet-domain sparsity[D].Vancouver:University of British Columbia,2009

[15] BUNKS C,SALECK F M,ZALESKI S,et al.Multiscale seismic waveform inversion[J].Geophysics,1995,60(5):1457-1473

[16] BURKE J V.A robust trust region method for constrained nonlinear programming problems[J].SIAM Journal on Optimization,1992,2(2):325-347

[17] LI X,ARAVKIN A Y,VAN LEEUWEN T,et al.Fast randomized full-waveform inversion with compressive sensing[J].Geophysics,2012,77(3):A13-A17

[18] HORST R,PARDALOS P,VAN THOAI N.Introduction to global optimization[M].New York:Springer,2000:1-353

[19] ABUBAKAR A,VAN DEN BERG P M.The contrast source inversion method for location and shape reconstructions[J].Inverse Problems,2002,18(2):495

[20] TU N,ARAVKIN A Y,VAN LEEUWEN T.Fast least-squares migration with multiples and source estimation[J]. Expanded Abstracts of 75thEAGE Conference and Exhibition,2013:1-5

[21] HERRMANN F J,LI X.Efficient least-squares imaging with sparsity promotion and compressive sensing[J].Geophysical Prospecting,2012,60(4):696-712

[22] HENNENFENT G,VAN DEN BERG E,FRIEDLANDER M P,et al.New insights into one-norm solvers from the Pareto curve[J].Geophysics,2008,73(4):A23-A26

[23] LI X,HERRMANN F J.Full-waveform inversion from compressively recovered model updates[J].SEG Technical Program Expanded Abstracts,2010:1029-1033

[24] VAN DEN BERG E,FRIEDLANDER M P.Probing the pareto frontier for basis pursuit solutions[J].SIAM Journal on Scientific Computing,2008,31(2):890-912

(編輯:朱文杰)

Full-waveform inversion from compressively recovered updates

LI Xiang

(UniversityofBritishColumbia,VancouverV6T1Z4,Canada)

Although full waveform inversion technique has been successfully applied,the amount of calculation of least-squares non-convex optimization problem is still a challenge.Random sampling technology reduces the number of shot and frequency,and save the full waveform inversion calculation greatly,but it brings curse of dimensionality and departure from Moore’s Law.In this paper with the successful improvement of full-waveform inversion,the current trend of incessantly pushing for higher quality models in increasingly complicated regions of the Earth reveals fundamental shortcomings in our ability to handle increasing problem size numerically.Two main culprits can be identified.First,there is the so-called curse of dimensionality exemplified by Nyquist’s sampling criterion,which puts disproportionate strain on current acquisition and processing systems as the size and desired resolution increases.Secondly,there is the recent departure from Moore’s law that forces us to lower our expectations to compute ourselves out of this.In this paper,we address this situation by randomized dimensionality reduction,which we adapt from the field of compressive sensing.In this approach,we combine deliberate randomized subsampling with structure-exploiting transform-domain sparsity promotion.Our approach is successful because it reduces the size of seismic data volumes without loss of information.With this reduction,we compute Newton-like updates at the cost of roughly one gradient update for the fully-sampled wavefield.Sparsity constrain is employed in the model update in inversion without changing the target function of waveform inversion and suppressing the virtual image noise raised by sub-sampling.The North Sea model testing result proves the feasibility and validity of the method.

compressive sensing,full-waveform inversion,curvelet transform,sparsity promoting

2016-10-11;改回日期:2016-11-03。

李翔(1988—),男,博士,主要從事全波形和最小二乘偏移成像研究。

P631

A

1000-1441(2017)01-0020-06

10.3969/j.issn.1000-1441.2017.01.002