李 曙，李 冉，金鵬程，曹繼華，陳澤平，周元玲，楊 喜

(吉首大學(xué)信息科學(xué)與工程學(xué)院，湖南 吉首 416000)

1 波阻抗反演方法研究現(xiàn)狀

地震波阻抗是一種重要的巖石屬性，它與巖性、孔隙度等巖石物理參數(shù)密切相關(guān)，因此如何獲取一直是勘探地震學(xué)領(lǐng)域的研究熱點.疊后地震反演是得到波阻抗的有效方法，傳統(tǒng)的基于Tikhonov正則化的波阻抗反演方法，通過求解L2范數(shù)最優(yōu)化問題來得到波阻抗.這類方法具有簡單和理論成熟的優(yōu)點，但是存在反演結(jié)果易受噪聲影響、垂直分界面過于光滑等問題，從而導(dǎo)致反演得到的波阻抗剖面的垂向分辨率不高、多解性較大等.為了克服這些問題，地球物理學(xué)家們嘗試引入新的理論.近年來，以壓縮感知和稀疏最優(yōu)化為代表的新理論引起了他們的注意.基于壓縮感知的算法，在只有少數(shù)觀測數(shù)據(jù)的前提下，能夠以極高的精度對稀疏(或變換域稀疏)信號進行重構(gòu).[1]基于此，大量地球物理研究人員提出了基于壓縮感知理論的地震信號噪聲衰減、地震反褶積、疊前彈性參數(shù)反演和波阻抗反演等方法[2-6].在波阻抗反演方面，2010年，Wang Y F[7]提出了一種基于L1范數(shù)正則化和梯度下降的方法，該方法驗證了L1范數(shù)正則化在波阻抗反演方面的可行性.L1范數(shù)正則化可以提高反演結(jié)果的抗噪聲性能，但是沒有充分利用地下巖層的先驗信息，如一般認(rèn)為巖層是呈塊狀分布的，即同一巖層內(nèi)部的巖石屬性近似.為了將這一信息引入波阻抗反演以得到更符合地質(zhì)先驗的結(jié)果，2014年，Zhang F C等[8]提出了基于L1范數(shù)誤差函數(shù)和全變分(Total Variation,TV)正則化的波阻抗反演方法，將數(shù)據(jù)保真項用L1范數(shù)表征的同時將波阻抗的全變分作為約束項，利用迭代重加權(quán)最小二乘方法求出反射系數(shù)后，借助波阻抗和反射系數(shù)之間的遞推關(guān)系求出波阻抗.實驗結(jié)果表明：在已知0時刻的波阻抗值情況下，該方法對高斯噪聲和脈沖噪聲具有較好的魯棒性；全變分正則化的加入使得該方法能更清晰地刻畫地層的邊界.與Zhang F C等用L1范數(shù)表征數(shù)據(jù)保真項的做法不同，2015年，Liu C等[9]提出了用L2范數(shù)表征數(shù)據(jù)保真項并對反射系數(shù)項施加L1范數(shù)約束的波阻抗反演方法.與傳統(tǒng)的阻尼最小二乘反演方法相比，基于L1范數(shù)正則化的波阻抗反演方法具有明顯的抗噪聲優(yōu)勢.

上述波阻抗反演方法只考慮了單道情況，沒有考慮地震數(shù)據(jù)各道之間的聯(lián)系.H Hamid等[10]發(fā)現(xiàn)了這一問題，提出了具有橫向約束的多道波阻抗同時反演方法.研究結(jié)果顯示，多道同時反演能充分利用地震數(shù)據(jù)的橫向連續(xù)性信息，提高反演結(jié)果的橫向連續(xù)性和抗噪聲性能，從而揭示更多的地質(zhì)結(jié)構(gòu)特征.A Gholami[11]提出了一種基于全變分正則化的非線性多道波阻抗同時反演方法，在進行多道同時反演時，利用地震數(shù)據(jù)在時空信息方面的特點和全變分正則化在獲得塊狀波阻抗方面的優(yōu)勢，并采用DCT變換提升計算速度.Yuan S Y等[12]提出了具有變換域稀疏提升的多道波阻抗同時反演方法，利用信號和噪聲在變換域中系數(shù)值不同的特點來抑制噪聲.實驗結(jié)果表明，利用該方法可以穩(wěn)定反演過程，提高反演結(jié)果的橫向連續(xù)性，對高波數(shù)噪聲有較強的抑制能力.

前面介紹的反演方法都只是聯(lián)合使用1～2種正則化方法，沒有充分利用先驗信息.在這種背景下，Li S等[4]提出了具有多種正則化約束的地震波阻抗反演方法，綜合了反射系數(shù)的L1范數(shù)稀疏約束、波阻抗的各向同性全變分約束、初始模型約束和多道同時反演等手段.實驗結(jié)果顯示，與基于L1范數(shù)正則化或全變分正則化的多道同時反演方法相比，該方法具有更優(yōu)的抗噪聲性能和更小的反演誤差，且能更清晰地刻畫地層的分界面.研究結(jié)果同時表明：與基于L1范數(shù)正則化多道同時反演方法相比，基于全變分正則化多道同時反演的反演效果有明顯提升；與基于全變分正則化多道同時反演方法相比，基于L1范數(shù)正則化和全變分正則化的多道同時反演方法的反演效果只是略有提升，卻面臨更復(fù)雜的正則化參數(shù)選擇問題.

上述反演方法中用到的全變分方法都是各向同性全變分(Isotropic Total Variation,ITV)方法，與ITV正則化方法相比，基于L1范數(shù)最小化的各向異性全變分正則化(Anisotropic Total Variation,ATV)方法對地震信號的稀疏性描述能力更強且實現(xiàn)簡單.由于地球的濾波效應(yīng)將導(dǎo)致地表接收到的地震信號的低頻成分缺失，而在反演算法中加入初始模型約束可有效緩解這一問題，因此筆者將ATV正則化和初始模型約束相結(jié)合并采用多道同時反演的方式，實現(xiàn)基于ATV正則化和初始模型約束的多道波阻抗同時反演.

2 波阻抗反演方法

2.1 正演模型

2.2 反演方法

筆者將ATV正則化與初始模型約束相結(jié)合，提出了基于ATV正則化與初始模型約束的多道波阻抗同時反演方法(簡稱新方法).新方法的目標(biāo)函數(shù)為數(shù)據(jù)保真項、ATV正則化項和初始模型約束項，即

(1)

等式右邊分別為數(shù)據(jù)保真項、ATV正則化項和初始模型約束項.為了求解該目標(biāo)函數(shù)，筆者采用split Bregman迭代方法[13]，(1)式可進一步表示為

(2)

(3)

(4)

(5)

(6)

至此，目標(biāo)函數(shù)的求解已完成.為了更清晰地說明新方法的執(zhí)行過程，下面給出算法的偽代碼：

Input：μ，λ，α，tol，L，S，L′

While

k=k+1

End

End

3 仿真實驗

為了測試新方法的效果，筆者利用Marmousi2模型[14]中的部分?jǐn)?shù)據(jù)進行仿真實驗，真實的波阻抗模型如圖1所示.從圖1可以看出，深度方向有450個采樣點，距離方向有500道.利用圖1所示的真實波阻抗剖面先求得反射系數(shù)剖面(圖2)，然后用反射系數(shù)剖面和一個主頻為40 Hz的雷克子波褶積產(chǎn)生合成地震記錄，再向該合成記錄中添加20%的高斯隨機噪聲，產(chǎn)生反演所需的合成地震剖面(圖3)，最后對真實波阻抗模型進行高斯低通濾波，并將濾波后的結(jié)果作為初始波阻抗模型(圖4).

圖1 真實的波阻抗模型Fig. 1 True Acoustic Impedance Model

圖2 反射系數(shù)剖面Fig. 2 Reflectivity Section

圖3 合成地震剖面Fig. 3 Synthetic Seismic Section

圖4 初始波阻抗模型Fig. 4 Initial Acoustic Impedance Model

設(shè)置參數(shù)μ=4×10-5,λ=4.5×10-3,α=0.6×10-5，反演結(jié)果如圖5所示.為便于觀察單道的反演情況，筆者抽取第200道的反演結(jié)果，并與相應(yīng)位置的真實波阻抗曲線、初始模型波阻抗曲線作對比，結(jié)果如圖6所示.

圖5 反演結(jié)果Fig. 5 Inversion Result

圖6 第200道的反演結(jié)果Fig. 6 Inversion Result of Trace 200

將圖5的反演結(jié)果剖面與圖1的真實波阻抗剖面進行對比，就可發(fā)現(xiàn)即使在較強噪聲干擾的情況下，新方法仍然能得到較好的反演效果：不僅能清晰地勾勒出大多數(shù)層位，而且反演結(jié)果與真實值之間的相對誤差也較小.圖6所示的單道反演結(jié)果，從微觀層面表明新方法的反演結(jié)果與真實波阻抗值之間的誤差較小.

4 結(jié)語

提出了一種基于ATV正則化與初始模型約束的多道波阻抗同時反演方法，利用該方法既能得到塊狀的反演結(jié)果，又能緩解地震信號缺失低頻信息的問題.此外，為了求解目標(biāo)函數(shù)，筆者提出了一種基于split Bregman的迭代算法，并選用Marmousi 2模型中的部分?jǐn)?shù)據(jù)進行測試.測試結(jié)果表明，即使在噪聲干擾較大的情況下仍然能獲得較好的反演結(jié)果，反演結(jié)果與真實值之間的相對誤差較小.