基于L1-2正則化的地震波阻抗“塊”反演

耿偉恒 陳小宏* 李景葉 湯 韋 吳 凡 張俊杰

(①中國石油大學(北京)地球物理學院，北京 102249；②油氣資源與探測國家重點實驗室，北京 102249)

0 引言

地震反演技術是提取地層屬性的有效手段之一，其基本過程是根據地表接收的地震數據反推地層的物理參數，從而預測儲層及識別流體[1]。根據所使用的地震數據體，地震反演分為疊后反演和疊前反演兩大類[2-4]。

疊后波阻抗(acoustic impedance, AI)反演是地震反演中的一種重要技術，首次被Lindseth[5]提出，隨后在地球物理反演領域得到快速發展[6-8]。根據實現方式的不同，AI反演主要分為稀疏脈沖反演、基于模型約束的反演、遞歸反演以及有色反演四類，不同方法各有其適用條件及應用范圍[2,9]。盡管AI反演歷經多年且得到很大發展，但仍存在諸多問題，如反問題求解的不適定性、常規反演算法對地層邊界刻畫不清晰以及遞歸反演中的誤差累積等[10-11]。

對于反問題求解的不適定問題，人們進行了研究并提出了相應的解決辦法[12-15]。早在1973年，Claerbout等[16]分析了L2范數和L1范數在反問題求解中的作用。基于L2范數約束，Tikhonov等[17]提出了一套完整的正則化方法解決不適定問題，其基本思想是在滿足數據殘差最小的條件下找到一個使模型的L2范數最小的解。盡管Tikhonov正則化很好地解決了不適定問題，但基于Tikhonov正則化獲得的解較平滑，難以滿足地學界對反演分辨率的要求。Buland等[18]將貝葉斯理論引入地球物理反問題，為基于概率的反演方法提供了數學基礎。Tarantola[19]給出了反問題求解的一些算法細節，包括蒙特卡洛方法以及最小二乘離散問題等。此后，各種范數被提出以在解決不適定問題的同時獲得較高分辨率的反演結果，包括L1范數[20]、Lp范數[21]以及L1-2范數等。L1-2范數作為一種稀疏度量，近年來廣泛用于地球物理反演[11, 22]。

此外，常規的遞歸反演算法首先利用疊后地震數據反演反射系數，然后利用道積分技術計算地層波阻抗信息。這種波阻抗反演算法存在兩個問題：一是對初始波阻抗依賴過高，如果第一層的波阻抗給定不準確，那么所有的波阻抗反演結果都會產生偏差；二是道積分中存在誤差累積，如果某些反演的反射系數存在誤差，那么這種誤差會累積到該目標層之后的所有波阻抗計算結果上。

為了解決波阻抗反演中存在的這些問題，本文提出了基于L1-2正則化的地震波阻抗“塊”反演方法。在前人的基礎上，將L1-2正則化引入基于模型的波阻抗反演，通過借鑒全變分(Total Variation, TV)正則化的思想，利用疊后地震數據直接獲得波阻抗反演結果。首先，推導線性化的波阻抗正演近似公式并分析精度。然后，基于貝葉斯理論，引入L1-2正則化構建波阻抗反演的目標函數，利用迭代重加權最小二乘算法(IRLS)求解目標函數，獲得波阻抗反演結果。由于波阻抗反演為單道反演算法，反演多道數據時道與道之間會產生空間不連續現象，因此文中對反演結果執行f-x域空間預測濾波，改善由噪聲和單道反演算法帶來的空間不連續性[23-24]。最后，利用合成數據和實際資料反演驗證所提方法的有效性和可行性。

1 線性化波阻抗正演公式

基于褶積模型，地震數據可以表示為反射系數和地震子波褶積的形式

(1)

式中：d(t)為地震數據，t為時間；w(τ)為地震子波，τ為時移量；r(t-τ)為地層反射系數；e(t)為噪聲。將式(1)寫為矩陣形式，有

D=WR+e

(2)

式中的變量分別為式(1)變量的矩陣形式。根據定義，當地下為多層水平介質時，任意第i個界面的反射系數為

(3)

式中Zi和Zi+1分別為第i層和第i+1層的波阻抗。將式(3)代入式(2)即可獲得由地層波阻抗得到地震數據的公式。不難看出，式(3)中反射系數與地層波阻抗之間為非線性關系，當與地震子波褶積后，得到的地震記錄與地層波阻抗之間的關系顯然也是非線性的，這種非線性關系會增加反演方程的計算復雜度，也會增強反演的不穩定性。因此需要將式(3)進行線性化近似。根據等價無窮小代換ex～1+x，由式(3)可得[25]

(4)

則有

(5)

將式(5)代入式(2)可得波阻抗正演公式

D=WSZ+e

(6)

式中：Z=ln(Z1,Z2,…,Znt)T為地層波阻抗的自然對數, nt為波阻抗采樣點數；S為差分算子，可表示為

(7)

由式(3)得到式(5)的過程由于使用了等價無窮小代換，因此有必要分析、對比線性化的波阻抗正演近似公式(式(5))和精確公式(式(3))的正演結果。圖1為由式(3)和式(5)獲得的反射系數。由圖可見：當反射系數r的絕對值在一定范圍(|r|≤0.3)時，式(5)與式(3)的結果幾乎完全吻合；當|r|較大時二者的結果才有所差異。因此，當地下反射界面為弱反射界面，即|r|≤0.3時，利用式(5)計算地層的反射系數在理論上不會產生較大誤差[20]。

圖1 由式(3)與式(5)獲得的反射系數

為了進一步了解式(5)的精度，給出了一個基于實際測井資料(圖2a)的正演結果(圖2)。可見，當反射系數的絕對值|r|≤0.3時(圖2b)，式(5)(圖2d)與式(3)(圖2c)的正演結果幾乎沒有差異(圖2e)，進一步證明式(5)作為線性化的波阻抗正演近似公式具有較高精度。

圖2 分別利用式(3)和式(5)對實際測井資料正演得到的地震記錄

2 波阻抗反演

根據貝葉斯理論，已知觀測數據d時，模型m的后驗概率分布P(m|d)等于模型參數的先驗概率分布P(m)與似然函數P(d|m)的乘積再除以邊緣概率密度P(d)，即

(8)

式中P(d)為一個常數，其作用是使后驗概率分布的總和等于1。對于波阻抗反演，假設噪聲服從高斯分布，則似然函數可以表示為

P(D|Z)=

(9)

式中：Ce為噪聲的協方差。阻抗模型Z的先驗分布可表示為

P(Z)=P0Zexp[-R(Z)]

(10)

P(Z|D)∝

(11)

其中λ為超參數。式(11)忽略了邊緣概率密度P(D)。使式(11)中指數項最大的解就是最大后驗概率解，因此可以建立反演目標函數

(12)

引入高斯先驗分布，則R(Z)可表示為

(13)

式中Z0為初始模型。在反演合成數據時，通過平滑真實波阻抗得到初始模型；在反演實際數據時，一般通過插值、外推測井數據得到初始模型。基于高斯先驗分布的反演結果一般較光滑，為此前人利用其他分布提高反演分辨率，如柯西分布、拉普拉斯分布(L1范數)以及L1-2范數。圖3為不同分布的一維概率密度函數對比。由圖可見：①高斯分布的峰值最寬，且沒有“長尾”分布(“長尾”保證了反演結果中允許存在某些“大”值，即保證“稀疏”，如圖中紅色框位置)，所以基于高斯分布的反演結果較光滑且抗噪性差。②柯西分布和拉普拉斯分布在提升反演結果稀疏度方面明顯優于高斯分布，二者的概率密度函數的峰值更尖銳且有明顯的“長尾”。③L1-2范數的概率密度函數的峰值更尖銳，且“長尾”更明顯，因此在理論上可以獲得更稀疏的反演結果。

圖3 不同分布的一維概率密度函數對比

為了更直觀地對比幾種不同的分布，繪制了L0范數、L1范數(對應拉普拉斯分布)、L2范數(對應高斯分布)以及L1-2范數的二維圖像(圖4)。可見，相比于L2范數(圖4c)和L1范數(圖4b)，L1-2范數(圖4d)與L0范數(圖4a)更相似，因此理論上基于L1-2范數約束可以獲得比L2范數和L1范數約束更稀疏、分辨率更高的反演結果。

圖4 L0范數(a)、L1范數(b)、L2范數(c)以及L1-2范數(d)的二維圖像

以二維向量(x，y)為例，幾種范數的定義為

(14)

式中x、y為二維向量的分量。L0范數定義為計算模型中非零項的個數，也被認為是最稀疏的度量。因此可通過對比其他范數圖像與L0范數圖像的相似程度判斷范數的稀疏性。

當引入L1-2正則化時，在式(13)的基礎上，R(Z)可以表示為

(15)

式中：α≥0為權重參數，該值越大，L1-2范數越稀疏；η和β則分別控制高斯項和L1-2正則化項在目標函數中的比例。由于L1-2范數是一種稀疏范數，而波阻抗并不滿足稀疏性假設，因此將式(15)代入式(12)進行反演必定會導致錯誤的反演結果。雖然波阻抗不滿足稀疏假設，但波阻抗的空間導數滿足稀疏條件，因此類似于TV正則化，本文通過對波阻抗的空間導數施加稀疏約束，從而獲得高分辨率的反演結果。這樣，基于L1-2正則化的波阻抗反演的目標函數可定義為

(16)

在反演過程中，通常需要多次試驗獲得參數η和β，以保證最終獲得較滿意的反演結果。利用IRLS算法求解式(16)，將式(16)對Z求導可得

(17)

令?J/?Z=0，可得最終的迭代更新式為

Zn+1=[(WS)T(WS)+ηI+w]-1×

[(WS)TD+ηZn]

(18)

式中：Zn為第n次迭代的解；w為重加權矩陣，可表示為

(19)

利用IRLS算法求解式(16)的算法流程為：

(1)設定超參數η和β、權重參數α以及最大迭代次數K；

(2)設定初始模型Z0；

(3)給定子波矩陣W及差分矩陣S；

(4)令n=1，根據式(19)計算w；

(5)根據式(18)進行求解，并令n=n+1；

(6)重復步驟(4)和步驟(5)，直至達到最大迭代次數。

值得注意的是，K的選擇對于反演結果的收斂以及算法的效率極為重要。一般來說，迭代10次左右即可獲得滿意的反演結果，但K的選擇還與模型的復雜度以及數據的信噪比有關。在反演實際地震數據時，可通過對比井旁道的反演結果與實際測井資料設定較合適的K。

3 模型數據測試

利用單道塊狀模型和二維SEG/EAGE推覆體模型，分別對比、分析基于L2范數、L1范數以及L1-2范數約束的波阻抗反演結果，以驗證所提方法的效果。在反演推覆體模型時，對反演結果執行f-x域空間預測濾波以提高反演結果的橫向連續性。文中信噪比定義為

(20)

式中：s*為無噪地震數據；s為含噪地震數據。

利用式(3)得到單道塊狀模型反射系數，與主頻為30Hz的雷克子波褶積得到合成地震記錄(圖5)，然后利用合成地震記錄反演得到波阻抗(圖6)。從結果來看：①由于沒有施加稀疏約束，且數據含噪，對于塊狀模型，基于L2范數約束的反演效果明顯較差(圖6a、圖6d)。②基于L1范數(圖6b、圖6e)和L1-2范數(圖6c、圖6f)約束的反演結果均呈明顯的“塊狀”特征，反演效果優于L2范數，其中L1-2范數約束的反演結果與真實模型吻合度更高(圖6b、圖6c、圖6e、圖6f的紅色箭頭處)。③隨著噪聲能量增強，反演質量均下降。如基于L2范數約束的反演結果仍然無法精確刻畫地層的“塊狀”特征(圖6d)；基于L1范數約束的反演結果與真實值也存在偏差(圖6e)；基于L1-2范數約束的反演結果雖然在某些層位與真實值稍有差別，但總體上吻合較好，且很好地刻畫了波阻抗的“塊狀”特征(圖6f)，證明了所提方法的有效性。

圖5 基于單道塊狀模型的合成地震數據

圖7為基于圖6數據計算的反射系數。由圖可見：①無論信噪比為10dB還是4dB，基于L2范數約束的反演效果均較差，反演結果與真實的反射系數相差很大。這是由于L2范數為平滑約束，其概率分布函數在均值μ附近的峰值較寬(圖3)，因此反演結果不能很好地壓制一些“小值”，即抗噪性較差，而且概率密度函數“長尾”不明顯，從而使某些離散點被壓制，因此不能反演幅值較大的反射系數。②L1范數和L1-2范數約束的反演效果明顯變好，其中L1-2范數由于更好地壓制了“小值”，同時也保護了“大值”，因此基于L1-2范數約束的反演效果優于基于L1范數約束，尤其是當信噪比為4dB時(圖7b紅色箭頭處)。

圖6 基于圖5數據反演的波阻抗

圖7 基于圖6數據計算的反射系數

利用式(3)得到SEG/EAGE推覆體模型(圖8b)反射系數，與主頻為40Hz的雷克子波褶積并加入信噪比為4dB的隨機噪聲得到合成地震記錄(圖8a)。圖8為基于SEG/EAGE推覆體模型的反演結果。從反演結果可見：①L1-2范數約束的反演結果的分辨率與精度最高(圖8f)，與真實模型(圖8b)吻合很好；基于L1范數約束的反演效果略差(圖8e)，但仍然好于基于L2范數約束(圖8d)。②由于單道算法的限制，反演結果中存在明顯的橫向不連續現象，尤其是當加入稀疏約束時，這種現象更明顯(圖8e、圖8f)。因此，有必要改善反演結果的橫向連續性。③三種范數約束的波阻抗反演結果(圖8d～圖8f)經過f-x域空間預測濾波后，橫向連續性均有提升(圖8g～圖8i)，尤其是圖8h和圖8i。該結果進一步證明了本文所提方法的有效性和可行性。

圖8 基于SEG/EAGE推覆體模型的反演結果

4 實際資料測試

基于M區實際地震資料測試本文方法。圖9為M區疊后地震資料反演結果。可見：①基于L2范數約束的反演結果分辨率較低(圖9d)。②基于L1范數(圖9e)和L1-2范數約束(圖9f)的反演效果得到明顯改善，反演分辨率更高，對地層的刻畫更清晰，且后者優于前者(紅色箭頭處)。圖10為井旁道反演結果與經過Backus平均的實際測井資料對比結果。可見，基于L1-2范數約束的井旁道反演結果與測井曲線更吻合，尤其表現在對“塊狀”特征的反演(紅色箭頭處)。此外，通過求取目的層段(3.0～3.2s)反演結果與真實測井數據的相關系數，可以定量對比不同范數約束的反演結果差異。向量X與Y的相關系數定義為

圖9 M區疊后地震資料反演結果

(21)

式中：c(X，Y)為X與Y的協方差；c(X，X)和c(Y，Y)分別為X和Y的方差。

基于L2范數、L1范數及L1-2范數的反演結果與真實測井數據的相關系數分別為0.6529、0.5197、0.7264，表明：基于L1-2范數的反演結果與真實測井數據更吻合，這是由于L1-2范數約束為稀疏約束(圖10的紅色箭頭處)，可很好地反演地層的“塊狀”信息；基于L2范數和L1范數約束的反演效果則相對較差，其中基于L1范數約束的反演結果的相關系數小于基于L2范數約束，可能原因是前者的反演結果盡管在稀疏性方面要優于后者，但前者的反演結果在大約3.14s處存在一個明顯的層位漂移，從而導致反演結果與真實測井數據的相關系數較低。

圖10 井旁道反演結果與經過Backus平均的實際測井資料對比結果

5 結論

盡管波阻抗反演技術已經相當成熟，但仍然存在一定問題，如反問題的不適定性、反演的分辨率低以及對地層邊界刻畫不清晰等。針對這些問題，本文提出了基于L1-2正則化的地震波阻抗“塊”反演方法。首先推導線性化的波阻抗正演近似公式并進行精度分析，證明近似式的精度較高。然后介紹了基于L1-2范數約束的波阻抗反演方法，與傳統的基于L2范數和L1范數約束的反演方法相比，基于L1-2范數約束的反演結果明顯提高了分辨率，反演結果更具“塊狀”特征。相關系數的定量對比也證明了基于L1-2范數的反演結果優于基于L1范數和L2范數。此外，文中引入f-x域空間預測濾波改善單道反演算法帶來的橫向不連續現象。顯而易見的是，通過對反演結果執行f-x域空間預測濾波，明顯提高了反演結果的橫向連續性。合成和實際地震數據反演證明了所提方法的可行性和有效性。