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

李 翔

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

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

李 翔

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

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

壓縮感知;波形反演;曲波變換;稀疏促進(jìn)

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

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

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

1 方法原理

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

全波形反演可表示為:

(1)

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

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

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

然而利用高斯-牛頓法求解上述波形反演問(wèn)題時(shí)計(jì)算量非常大,因?yàn)槊恳徊降夹枰蠼庖粋€(gè)線(xiàn)性的高斯-牛頓子問(wèn)題(即求解海森矩陣的逆)。而求解子問(wèn)題也需要多次進(jìn)行子迭代(此迭代不同于波形反演模型迭代),每次子問(wèn)題的子迭代需要多次求解波動(dòng)方程[17-18]。

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

(2)

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

(3)

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

(4)

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

(5)

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

1.2 稀疏反演迭代

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

對(duì)于大多數(shù)反演問(wèn)題,人們往往避免采用迭代法求解雅可比矩陣的逆,因?yàn)榈ㄔ诿看蔚鷷r(shí)需要不停地計(jì)算其正傳及其轉(zhuǎn)置,因而計(jì)算量巨大。此外,在很多優(yōu)化反演問(wèn)題中,雅可比矩陣往往是病態(tài)的(海森矩陣是奇異的)。然而與很多反演問(wèn)題不同,地震勘探中的海森矩陣是相對(duì)均衡的。因?yàn)樵诔跏寄Ｐ拖鄬?duì)理想的情況下,線(xiàn)性Born近似正演能夠很好地匹配非線(xiàn)性波動(dòng)方程正演,在此情況下海森矩陣的逆可以用簡(jiǎn)單的對(duì)角矩陣近似代替[18]。

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

(6)

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

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

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

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

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

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

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

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

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

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

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

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

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

圖3 真實(shí)的BG速度模型

圖4 初始速度模型

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

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

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

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

3 結(jié)論和展望

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

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

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

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

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