屠 寧

(同濟(jì)大學(xué)海洋與地球科學(xué)學(xué)院,上海200092)

壓縮感知理論的出現(xiàn)[6-7]為解決規(guī)模龐大的地震反演問題提供了新思路。在地震反演偏移成像問題中,對于一個壓縮的物理觀測模型,壓縮感知理論提出了:①精確求解原始觀測模型所對應(yīng)解的充分條件,從而指導(dǎo)我們最優(yōu)化設(shè)計(jì)壓縮觀測系統(tǒng);②如何利用稀疏優(yōu)化的方法從壓縮觀測的數(shù)據(jù)中恢復(fù)真實(shí)的地震像。本文將介紹作者與HERRMANN等在這個領(lǐng)域共同研究取得的一些重要成果[8-11]:①最小二乘逆時偏移的快速計(jì)算;②消除逆時偏移中未知的子波參數(shù)對反演結(jié)果的影響;③包含表面多次波的地震數(shù)據(jù)的成像。同時展示了合成和實(shí)際數(shù)據(jù)的實(shí)驗(yàn)結(jié)果,以驗(yàn)證本文方法的可行性和有效性。

1 方法和合成數(shù)據(jù)實(shí)驗(yàn)

1.1 快速最小二乘偏移

HERRMANN等[8]提出了基于壓縮感知的最小二乘偏移方法,TU等[9]在此基礎(chǔ)上,分析了波動方程線性化過程中引入的線性誤差對該方法的影響。與傳統(tǒng)的炮編碼一樣,基于壓縮感知的最小二乘偏移也是通過對數(shù)據(jù)降維來實(shí)現(xiàn)波動方程模擬量的減少。與傳統(tǒng)的最小二乘偏移不同,求解壓縮感知問題的Basis Pursuit De-Noise(BPDN)模型[6]引入了稀疏約束,其在頻率域求解的形式如下:

(1)

式中:Ω表示對數(shù)據(jù)降維時選擇的頻率子集;Σ表示選擇的編碼炮集;i和j分別表示頻率和炮集的編號;d表示隨機(jī)降采樣的觀測數(shù)據(jù)(在本文中,根據(jù)上下文可判斷出是僅包含一次波的觀測數(shù)據(jù),還是同時包含一次波和多次波的觀測數(shù)據(jù));m0表示光滑的背景模型;wi表示地震子波w在某個頻率上的頻譜值;s表示不包含子波效應(yīng)的隨機(jī)降采樣后的震源波場(即主要包含與觀測數(shù)據(jù)相對應(yīng)的震源的位置信息,子波效應(yīng)包含在w中,降采樣方式和觀測數(shù)據(jù)d一致);C表示曲波變換(C*表示曲波反變換);x是模型擾動δm在曲波域的系數(shù)(δm=C*x);σ是可調(diào)節(jié)參數(shù),用以控制匹配數(shù)據(jù)的精度,其選擇往往反映數(shù)據(jù)的信噪比。(1)式的意義是,在所有可能的曲波系數(shù)中,尋找一個最稀疏的系數(shù)x,經(jīng)過曲波反變換后得到模型擾動的一個估計(jì),再經(jīng)過反偏移,可以在由σ指定的誤差范圍內(nèi)匹配觀測數(shù)據(jù)d。引入曲波變換的意義在于,由于地質(zhì)模型在空間域并非稀疏,曲波變換可將其稀疏化[12-13]。直觀上解釋,上述問題的含義是:在所有可能的曲波系數(shù)中,尋找最稀疏(一范數(shù)最小)的解,這個解經(jīng)過反曲波變換(變換到空間域中),得到地下介質(zhì)的擾動模型,該模型經(jīng)過反偏移算子的作用,能在σ指定的誤差范圍內(nèi)匹配觀測到的地震數(shù)據(jù)。HERRMANN等[8]利用譜梯度投影算法(SPGL1)求解上述問題[14]。

對比傳統(tǒng)的最小二乘偏移算法,壓縮感知引入了稀疏約束(公式(1)第2行)。然而壓縮感知對最小二乘偏移的影響還不止于此:壓縮感知理論對數(shù)據(jù)如何進(jìn)行降采樣(同時對應(yīng)對炮集如何進(jìn)行編碼)也有指導(dǎo)作用。壓縮感知理論指出,在采樣過程中引入隨機(jī)因素,是最終恢復(fù)原始信號的關(guān)鍵條件。信號恢復(fù)的充分條件由Restricted Isometry Property(RIP)條件給出[6];在地震信號的處理中,HENNENFENT等[15]對隨機(jī)采樣的重要性亦作了直觀的闡述。因此,在對觀測數(shù)據(jù)進(jìn)行降采樣時,應(yīng)隨機(jī)選取頻率的子集;在編碼炮集時,也應(yīng)引入足夠的隨機(jī)因素。詳細(xì)的降采樣方案請參考文獻(xiàn)[8]和文獻(xiàn)[11]。

另一個問題就是波動方程正演的線性化(即用反偏移代替正演)帶來的誤差。在數(shù)據(jù)觀測較多,成像系統(tǒng)是超定系統(tǒng)時,最小二乘的方法可以自動地將正交于反偏移算子的列空間的噪聲排除在解空間之外;然而,隨著降采樣帶來的系統(tǒng)欠定性增強(qiáng),噪聲可能會在解空間投影,從而給真實(shí)解帶來干擾。TU等[9]詳細(xì)分析了線性誤差對成像系統(tǒng)帶來的影響,并提出,根據(jù)特定的規(guī)則在不同的迭代之間引入不同的降采樣算子(即對數(shù)據(jù)和炮集重新進(jìn)行采樣)會顯著提高成像系統(tǒng)對線性噪聲的魯棒性。

利用SEG/EAGE鹽丘模型的一個二維剖面,對上述算法的有效性作了驗(yàn)證[10]。實(shí)驗(yàn)?zāi)Ｐ偷膮?shù)如下:模型3.9km深,15.7km長,網(wǎng)格間距約為24m。共設(shè)置323個炮點(diǎn)和檢波點(diǎn),利用互易原理,將數(shù)據(jù)設(shè)置為每一炮都有共同的檢波器來接收。真實(shí)速度模型及其平滑之后的背景模型如圖1所示。圖2展示了一系列實(shí)驗(yàn)結(jié)果。首先,利用線性數(shù)據(jù)(即利用反偏移算子作用到真實(shí)模型擾動上生成的地震數(shù)據(jù))分析無任何噪聲的狀況下所提算法的有效性。圖2a為對線性數(shù)據(jù)做單次逆時偏移的結(jié)果;圖2b為對線性數(shù)據(jù)做快速最小二乘偏移的結(jié)果。在快速最小二乘偏移中,所用合成炮的數(shù)量約為原炮數(shù)的5%,所用頻率子集的大小約為總頻率集合的15%,我們運(yùn)行60個SPGL1的迭代,在這樣的反演參數(shù)設(shè)置下,總的波動方程的求解數(shù)量相當(dāng)于利用所有數(shù)據(jù)進(jìn)行單次逆時偏移,即圖2a和圖2b所展示的結(jié)果包含幾乎相同的波動方程求解。接著,驗(yàn)證存在相干噪聲的情況下該算法的效果。用時間域的iWave波場正演引擎來模擬觀測波場[16],用我們的頻率域波場線性反偏移引擎來進(jìn)行成像,結(jié)果如圖2c所示。這兩種波場模擬引擎間存在數(shù)值算法上的差別,同時,iWave正演的觀測波場含有高階散射。在這種情況下,快速最小二乘成像方法依然取得了較好的效果。和圖2a相比,圖2b和圖2c均展現(xiàn)出了較高的分辨率。

圖1 快速最小二乘逆時偏移合成數(shù)據(jù)實(shí)驗(yàn)所用的真實(shí)速度模型(a)和光滑速度模型(b)

圖2 線性數(shù)據(jù)的單次逆時偏移結(jié)果(a)和線性數(shù)據(jù)(b)、iWave正演數(shù)據(jù)(c)的快速最小二乘偏移結(jié)果

1.2 震源子波估計(jì)

無論是單次逆時偏移還是最小二乘逆時偏移,都需要預(yù)先知道地震子波。如果地震子波估算不準(zhǔn),那么正傳波場就會有誤,從而影響地震成像以及最小二乘偏移中梯度的求取。求取可靠的地震子波在技術(shù)上或者計(jì)算上面臨較大的挑戰(zhàn):依賴統(tǒng)計(jì)信息[17]或者大量的計(jì)算[18]。最小二乘偏移本身是一種反演算法,因此我們將子波的反演融入其中,在不給出子波先驗(yàn)信息的情況下,得到可靠的地震成像結(jié)果。

借助在非線性最小二乘問題中變量投影(Variable Projection)的方法[19-20],并將其理論推廣到包含稀疏約束的非線性最小二乘問題中,提出了不依賴于子波信息的最小二乘逆時偏移算法[10]。對于子波w來說,一旦x確定,存在下述解析表達(dá):

(2)

使用交替求解的方式,每更新一次x,就更新一次w的估計(jì),即可在未給定子波的情況下求解最小二乘偏移的梯度[10]。

繼續(xù)利用前文所用的模型進(jìn)行合成數(shù)據(jù)的實(shí)驗(yàn),主要結(jié)果如圖3所示。圖3a展示了利用線性化數(shù)據(jù)進(jìn)行最小二乘逆時偏移的結(jié)果,但正演中使用的地震子波和真實(shí)的子波存在一個相位差,其影響是反射層發(fā)生錯位(如箭頭所示)和整體地震像散焦。圖3b展示了利用相同的線性化數(shù)據(jù),在未給定子波(第一個迭代中的子波起始估計(jì)是白譜子波，即時間域的脈沖)的情況下,利用本文包含子波估計(jì)的成像方法得到的地震像。圖3c展示了利用iWave正演的數(shù)據(jù)進(jìn)行實(shí)驗(yàn)的結(jié)果(其它參數(shù)與圖3b所用的一致)。從結(jié)果來看,本文方法可以在未知真實(shí)子波的情況下,得到可靠的地震像。對比圖3b,圖3c和圖2b,圖2c可以看出,采用本文方法得到的地震像與知道真實(shí)子波情況下得到的地震像具有相似的高精度和分辨率。

1.3 表面多次波成像

海洋表面多次波是海洋地震資料處理時面臨的強(qiáng)干擾。逆時偏移本身基于單次散射的假設(shè),并不能對多次波進(jìn)行有效成像。我們從描述自由表面反射的自由表面多次波消除公式出發(fā),研究出對多次波進(jìn)行成像以及一次波和多次波聯(lián)合成像的方法[10-11]。海洋表面多次波的正確成像依賴于:①正確識別表面多次波對應(yīng)的震源波場。一次波的成像中,正傳波場為點(diǎn)震源波場;多次波的成像中,正傳波場是海洋表面的整體下行波場,可以看作一個面震源[11,21-23]。②新的成像條件。對于含有多階海洋表面反射的多次波來說,廣泛應(yīng)用于一次波成像的互相關(guān)成像條件會造成各階次反射之間產(chǎn)生非物理的串?dāng)_噪聲[24-25],形成對正確地震像的干擾。應(yīng)用反褶積成像條件可以一定程度緩解這個問題[24,26],而真正衰減串?dāng)_噪聲則需要利用最小二乘反演的方法[11,27]。

基于我們提出的快速最小二乘逆時偏移的框架,實(shí)現(xiàn)對表面多次波的正確成像。通過對正傳震源波場做修正,用自由表面的下行波場代替點(diǎn)震源波場,可實(shí)現(xiàn)對表面多次波的正確成像[11]以及一次波和多次波的聯(lián)合成像[10,28]。數(shù)據(jù)的預(yù)處理細(xì)節(jié)參考文獻(xiàn)[10]。本文主要介紹一次波和多次波聯(lián)合成像的方法。這時,數(shù)據(jù)匹配的目標(biāo)變?yōu)?

(3)

需要注意的是,此處的觀測數(shù)據(jù)d同時包含了上行的一次波波場和多次波波場。對于預(yù)測數(shù)據(jù),可以看到,在反偏移算子中,震源波場wisj-di,j不只是包含點(diǎn)震源s(經(jīng)過波場傳播后預(yù)測一次波),同時包含廣義的面震源波場即觀測數(shù)據(jù)d(前面的“-”表示水面的反射系數(shù)為-1,上行的觀測數(shù)據(jù)經(jīng)過水面反射,再經(jīng)過和一次波相同的波場傳播后,預(yù)測表面多次波)。除去震源波場的變化,其它成像過程與前文所述一次波的反演成像過程相同。另外,這種情況下,由于震源波場的變化,子波的估計(jì)方法則為:

(4)

選用專為驗(yàn)證與表面多次波有關(guān)的算法而設(shè)計(jì)的Sigsbee 2B模型,對包含海洋表面多次波的數(shù)據(jù)進(jìn)行了測試。真實(shí)和背景速度模型見圖4。數(shù)據(jù)正演使用iWave,使用自由表面邊界條件來正演包含多次波和鬼波的觀測數(shù)據(jù)。假設(shè)真實(shí)地震子波已知,利用本文算法對一次波和多次波進(jìn)行聯(lián)合快速最小二乘偏移成像,結(jié)果如圖5a所示,其整體的波動方程模擬次數(shù)依然控制在單次逆時偏移的水平上。圖5b 展示了通過子波估計(jì)對一次波和多次波進(jìn)行聯(lián)合成像的結(jié)果。可見,圖5a和圖5b結(jié)果具有很高的一致性,與真實(shí)模型對比可以看出,這兩個結(jié)果中多次波帶來的串?dāng)_噪聲均得到了很好的壓制。文獻(xiàn)[10]中詳細(xì)討論了本文方法對子波初始估計(jì)的魯棒性,本文方法在子波具有錯誤的初始振幅和相位的情況下,均能得到高精度的地震像。

圖4 包含多次波數(shù)據(jù)偏移實(shí)驗(yàn)所用真實(shí)的Sigsbee 2B速度模型(a)和光滑后的速度模型(b)

圖5 使用真實(shí)子波(a)和子波估計(jì)方法(b)得到的一次波與多次波聯(lián)合快速最小二乘偏移成像結(jié)果

2 實(shí)際數(shù)據(jù)實(shí)驗(yàn)

利用英國北海海域某工區(qū)實(shí)際數(shù)據(jù),對文中所述算法進(jìn)行驗(yàn)證。圖6a為利用互相關(guān)成像條件對多次波進(jìn)行偏移成像的結(jié)果(成像過程中,震源為面震源,由表面位置的下行波場構(gòu)成,反傳數(shù)據(jù)為多次波),可

以看到不同階次多次波導(dǎo)致的串?dāng)_噪聲。圖6b為利用本文快速最小二乘偏移算法對多次波成像的結(jié)果。圖7和圖8分別為圖6中區(qū)域A和區(qū)域B放大后的結(jié)果?？梢钥闯?圖6a中出現(xiàn)的極強(qiáng)的串?dāng)_噪聲在圖6b中得到了很好的壓制。

圖9為不同數(shù)據(jù)的偏移成像結(jié)果。圖9a是一次波成像結(jié)果;圖9b是多次波單獨(dú)成像結(jié)果;圖9c是一次波與多次波聯(lián)合成像的結(jié)果?？梢钥闯?多次波單獨(dú)成像和一次波成像存在較好的可比性,較強(qiáng)的連續(xù)反射層(紅色框中的反射層)均能得到較高精度的成像,但是整體信噪比較低,有些反射層(黑色箭頭所指)沒有成像;對于某些一次波成像沒能揭示出的結(jié)構(gòu)(白色箭頭所指),多次波成像能較好的刻畫出來。從圖9c可以看出,一次波與多次波的聯(lián)合成像,不僅在數(shù)據(jù)處理環(huán)節(jié)避免了多次波和一次波的分離,而且改進(jìn)了多次波成像本身信噪比偏低的缺點(diǎn)。

圖6 利用互相關(guān)成像條件(a)和快速最小二乘逆時偏移(b)對多次波進(jìn)行成像的結(jié)果

圖7 圖6a(a)和圖6b(b)中區(qū)域A放大后的結(jié)果

圖8 圖6a(a)和圖6b(b)中區(qū)域B放大后的結(jié)果

圖9 不同數(shù)據(jù)的偏移成像結(jié)果a 僅用一次波; b 僅用多次波;c 一次波和多次波聯(lián)合成像

3 結(jié)論

本文將壓縮感知應(yīng)用于地震成像中,通過人工壓縮觀測系統(tǒng)和觀測數(shù)據(jù),大幅降低最小二乘逆時偏移的成本?？焖僮钚《四鏁r偏移是個靈活的框架,在這個框架下,介紹了如何利用變量投影技術(shù)解決地震偏移中未知子波的問題,并利用廣義面震源實(shí)現(xiàn)了海洋多次波的成像以及一次波與多次波的聯(lián)合成像。模型數(shù)據(jù)以及海上實(shí)際數(shù)據(jù)應(yīng)用結(jié)果表明,本文方法在未知地震子波信息的情況下,利用相當(dāng)于單次逆時偏移的計(jì)算量,可以得到高分辨率的最小二乘偏移成像結(jié)果,并對多次波進(jìn)行正確的歸位和成像,為今后更加有效和多樣化地利用多次波信息提供了思路。

[1] LAILLY P.The seismic inverse problem as a sequence of before stack migrations[C]∥BEDNAR J B.Conference on inverse scattering:theory and application,Philadelphia:SIAM,1983:206-220

[2] NEMETH T,WU C,SCHUSTER G T.Least-squares migration of incomplete reflection data[J].Geophysics,1999,64(1):208-221

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

[4] ROMERO L,GHIGLIA D,OBER C,et al.Phase encoding of shot records in prestack migration[J].Geophysics,2000,65(2):426-436

[5] SCHUSTER G T,WANG X,HUANG Y,et al.Theory of multisource crosstalk reduction by phase-encoded statics[J].Geophysical Journal International,2011,184(3):1289-1303

[6] DONOHO D.Compressed sensing[J].IEEE Transactions on Information Theory,2006,52(4):1289-1306

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

[9] TU N,LI X,HERRMANN F J.Controlling linearization errors in l_1 regularized inversion by rerandomization[J].Expanded Abstracts of 83rdAnnual Internat SEG Mtg,2013:4640-4644

[10] TU N,ARAVKIN A Y,LEEUWEN T,et al.Source estimation with surface-related multiples-fast ambiguity-resolved seismic imaging[J].Geophysical Journal International,2016,205(3):1492-1511

[11] TU N,HERRMANN F J.Fast imaging with surface-related multiples by sparse inversion[J].Geophysical Journal International,2015,201(1):304-317

[13] HERRMANN F J,MOGHADDAM P P,STOLK C.Sparsity- and continuity-promoting seismic image recovery with curvelet frames[J].Applied and Computational Harmonic Analysis,2008,24(2):150-173

[14] 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

[15] HENNENFENT G,HERRMANN F J.Simply denoise:wavefield reconstruction via jittered undersampling[J].Geophysics,2008,73(3):V19-V28

[16] TERENTYEV I S,VDOVINA T,SYMES W W,et al.iWave:a framework for wave simulation[EB/OL].[2017-01-10].http:∥www.trip.caam.rice.edu/software/iwave/doc/html/index.html

[17] ULRYCH T J,VELIS D R,SACCHI M D.Wavelet estimation revisited[J].The Leading Edge,1995,14(11):1139-1143

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

[19] GOLUB G,PEREYRA V.The differentiation of pseudo-inverses and nonlinear least squares problems whose variables separate[J].SIAM Journal on Numerical Analysis,1973,10(2):413-432

[20] GOLUB G,PEREYRA V.Separable nonlinear least squares:the variable projection method and its applications[J].Inverse Problems,2003,19(2):R1-R26

[21] BERKHOUT A J.Migration of multiple reflections[J].Expanded Abstracts of 63rdAnnual Internat SEG Mtg,1993:1022-1025

[22] GUITTON A.Shot-profile migration of multiple reflections[J].Expanded Abstracts of 72ndAnnual Internat SEG Mtg,2002:1296-1299

[23] ZUBERI A,Alkhalifah T.Imaging by forward propagating the data:theory and application[J].Geophysical Prospecting,2013,61(S1):248-267

[24] MUIJS R,ROBERTSSON J O A,HOLLIGER K.Prestack depth migration of primary and surface-related multiple reflections:part Ⅰ—imaging[J].Geophysics,2007,72(2):S59-S69

[25] LIU Y,CHANG X,JIN D,et al.Reverse time migration of multiples for subsalt imaging[J].Geophysics,2011,76(5):WB209-WB216

[26] WHITMORE N D,VALENCIANO A A,SOLLNER W.Imaging of primaries and multiples using a dual-sensor towed streamer[J].Expanded Abstracts of 80thAnnual Internat SEG Mtg,2010:3187-3192

[27] WONG M,BIONDI B,RONEN S.Imaging with multiples using least-squares reverse time migration[J].The Leading Edge,2014,33(9):970-976

[28] TU N,HERRMANN F J.Fast least-squares imaging with surface-related multiples:application to a North-Sea data set[J].The Leading Edge,2015,34(7):788-794