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

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

(①中國(guó)石油大學(xué)(北京)地球物理學(xué)院，北京 102249；②油氣資源與探測(cè)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 102249)

0 引言

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

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

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

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

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

1 線性化波阻抗正演公式

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

(1)

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

D=WR+e

(2)

式中的變量分別為式(1)變量的矩陣形式。根據(jù)定義，當(dāng)?shù)叵聻槎鄬铀浇橘|(zhì)時(shí)，任意第i個(gè)界面的反射系數(shù)為

(3)

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

(4)

則有

(5)

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

D=WSZ+e

(6)

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

(7)

由式(3)得到式(5)的過程由于使用了等價(jià)無窮小代換，因此有必要分析、對(duì)比線性化的波阻抗正演近似公式(式(5))和精確公式(式(3))的正演結(jié)果。圖1為由式(3)和式(5)獲得的反射系數(shù)。由圖可見：當(dāng)反射系數(shù)r的絕對(duì)值在一定范圍(|r|≤0.3)時(shí)，式(5)與式(3)的結(jié)果幾乎完全吻合；當(dāng)|r|較大時(shí)二者的結(jié)果才有所差異。因此，當(dāng)?shù)叵路瓷浣缑鏋槿醴瓷浣缑?，即|r|≤0.3時(shí)，利用式(5)計(jì)算地層的反射系數(shù)在理論上不會(huì)產(chǎn)生較大誤差[20]。

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

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

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

2 波阻抗反演

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

(8)

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

P(D|Z)=

(9)

式中：Ce為噪聲的協(xié)方差。阻抗模型Z的先驗(yàn)分布可表示為

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

(10)

P(Z|D)∝

(11)

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

(12)

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

(13)

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

圖3 不同分布的一維概率密度函數(shù)對(duì)比

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

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

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

(14)

式中x、y為二維向量的分量。L0范數(shù)定義為計(jì)算模型中非零項(xiàng)的個(gè)數(shù)，也被認(rèn)為是最稀疏的度量。因此可通過對(duì)比其他范數(shù)圖像與L0范數(shù)圖像的相似程度判斷范數(shù)的稀疏性。

當(dāng)引入L1-2正則化時(shí)，在式(13)的基礎(chǔ)上，R(Z)可以表示為

(15)

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

(16)

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

(17)

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

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

[(WS)TD+ηZn]

(18)

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

(19)

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

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

(2)設(shè)定初始模型Z0；

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

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

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

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

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

3 模型數(shù)據(jù)測(cè)試

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

(20)

式中：s*為無噪地震數(shù)據(jù)；s為含噪地震數(shù)據(jù)。

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

圖5 基于單道塊狀模型的合成地震數(shù)據(jù)

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

圖6 基于圖5數(shù)據(jù)反演的波阻抗

圖7 基于圖6數(shù)據(jù)計(jì)算的反射系數(shù)

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

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

4 實(shí)際資料測(cè)試

基于M區(qū)實(shí)際地震資料測(cè)試本文方法。圖9為M區(qū)疊后地震資料反演結(jié)果?？梢姡孩倩贚2范數(shù)約束的反演結(jié)果分辨率較低(圖9d)。②基于L1范數(shù)(圖9e)和L1-2范數(shù)約束(圖9f)的反演效果得到明顯改善，反演分辨率更高，對(duì)地層的刻畫更清晰，且后者優(yōu)于前者(紅色箭頭處)。圖10為井旁道反演結(jié)果與經(jīng)過Backus平均的實(shí)際測(cè)井資料對(duì)比結(jié)果?？梢姡贚1-2范數(shù)約束的井旁道反演結(jié)果與測(cè)井曲線更吻合，尤其表現(xiàn)在對(duì)“塊狀”特征的反演(紅色箭頭處)。此外，通過求取目的層段(3.0～3.2s)反演結(jié)果與真實(shí)測(cè)井?dāng)?shù)據(jù)的相關(guān)系數(shù)，可以定量對(duì)比不同范數(shù)約束的反演結(jié)果差異。向量X與Y的相關(guān)系數(shù)定義為

圖9 M區(qū)疊后地震資料反演結(jié)果

(21)

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

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

圖10 井旁道反演結(jié)果與經(jīng)過Backus平均的實(shí)際測(cè)井資料對(duì)比結(jié)果

5 結(jié)論

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