基于ATV正則化與初始模型約束的波阻抗反演*

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

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

1 波阻抗反演方法研究現狀

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

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

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

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

2 波阻抗反演方法

2.1 正演模型

2.2 反演方法

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

(1)

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

(2)

(3)

(4)

(5)

(6)

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

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

While

k=k+1

End

End

3 仿真實驗

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

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

圖2 反射系數剖面Fig. 2 Reflectivity Section

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

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

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

圖5 反演結果Fig. 5 Inversion Result

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

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

4 結語

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