疊前重加權L1范數稀疏約束的地震反演方法

趙云，文曉濤*，尹川，韓文明，李陳龍

（1.成都理工大學油氣藏地質及開發工程重點實驗室，四川成都 610059；2.成都理工大學地球物理學院，四川成都 610059；3.中海石油國際能源服務（北京）有限公司，北京 100010）

0 引言

與疊后地震反演方法相比，疊前地震反演方法可獲得更多的彈性參數，能更有效地進行儲層預測和流體識別[1-3]。為了提高疊前地震反演結果的準確性，可通過約束方法[4-6]將先驗信息引入反演。Berkhout[7]利用反射系數的稀疏信息和初始模型的趨勢信息，應用L2范數約束添加先驗約束信息，提高了反演結果的準確性。然而，Zhang 等[8]發現采用L2范數約束先驗信息會出現低稀疏性，導致反演結果分辨率變低，相鄰巖層邊界位置呈現非聚焦性平滑過渡。由于地層反射系數向量通常具有稀疏性，向量中存在稀疏非零值，因此在稀疏表示過程中如何提高先驗約束信息的稀疏性非常重要。

Taylor等[9]提出了利用L1范數稀疏約束地震反演方法，隨后Thueune 等[10]、Kong 等[11]、Zhang 等[12]均驗證了L1范數的稀疏性優于L2范數。Guitton 等[13]優化了L1范數，提出了結合柯西約束的混合范數。李昕潔等[14]將該方法用于全波形反演，取得了較好的效果。Zhang 等[8]提出了基于L1范數優化方法的基追蹤反演(BPI)算法，傳統疊前AVO 反演邊界尖銳問題得到一定程度的解決，同時證實了L1范數約束的稀疏性更好。張雨強等[15]、Sun 等[16]在地球物理反演中引入了自適應Lp范數稀疏表示方法，可反演出“塊狀”地層，反演結果比L1范數縱向分辨率更高。Pérez 等[17]提出了一種加權混合L2,1范數的疊前AVO 反演方法，優化了反演函數約束項的稀疏性，可獲得更多的塊狀解。王治強等[18]采用了全變分(TV)約束波阻抗的方法，既可以有效壓制隨機噪聲，還能較好地刻畫地層邊界。Li等[19]基于分裂Bregman 迭代算法，將L1范數與全變分正則化結合，實現了更高橫向分辨率的疊前地震多道同時反演。Huang 等[20]、Wang 等[21]、耿偉恒等[22]將非凸L1-2正則化方法拓展到疊前AVO 反演，證明了L1-2范數比L1范數更稀疏。然而，上述稀疏約束忽略了巖層邊界中的地震反射振幅信息，導致反演結果出現速度偽層，為后續參數預測積累了誤差。

為此，本文提出了一種疊前重加權L1范數稀疏約束的地震反演方法(簡稱“本文方法”)，在反演目標函數稀疏項的構建中引入了基于反射系數表示的重加權矩陣，即“重加權L1范數”。相比L1范數稀疏約束(簡稱“傳統方法”)，重加權L1范數的優越性主要體現在：其稀疏約束表示可以通過重加權矩陣補充地層邊界信息。加權地層反射邊界信息的引入能夠提高L1范數約束的稀疏性，本文將重加權L1范數稀疏約束和初始低頻模型約束同時用于構造反演目標函數，采用交替方向乘子法(ADMM)和軟閾值收縮算法(ISTA)進行迭代求解[23-25]。理論模型測試和實際工區應用結果均表明，本文方法明顯提高了疊前多參數同時反演的分辨率，對層速度邊界刻畫更清晰，極大減弱了速度偽層和密度偽層現象。

1 方法原理

1.1 正演模擬

地震信號一般可以表示為

式中：S為地震記錄；ω為地震子波；R為地層反射系數；N為噪聲；“*”表示褶積運算。

引入子波卷積矩陣W，將式(1)中的褶積運算轉化為簡單線性運算，其中W定義為

式中m為地震子波的長度。式(1)可簡化為

基于Aki和Richard反射系數近似公式[26]

其中

針對式(5)，定義變量

式中：RVP、RVS、Rρ分別為縱波速度、橫波速度、密度反射系數；a(θ)、b(θ)、c(θ)皆為角度系數。

近似地將RVP表示為

將MVP記作VP的自然對數，即用矩陣形式表示

根據式(7)、式(8)，將縱波的反射系數表示成矩陣形式

式中D0為單參數一階差分矩陣，可定義為

同理可得

式中：MVS為橫波速度的自然對數矩陣；Mρ為密度的自然對數矩陣。

最終推導出正演模型的矩陣表達式為

式中：B為角度系數矩陣；D為差分矩陣；M和D分別定義為

1.2 反演目標函數構建

用于反演縱波速度、橫波速度和密度的自然對數與地震記錄之間的關系為

式中||·||2為L2范數。

基于傳統L1范數稀疏約束包含低頻背景約束項，傳統方法的反演目標函數為

式中：M0為低頻模型矩陣；λ為低頻模型約束項系數；α為L1范數稀疏約束項系數；||·||1為L1范數。

重加權矩陣數學表示形式[27]為

式中：Q為自適應重加權矩陣；qi為重加權矩陣對角線第i個元素；ri為反射系數向量R的第i個元素；N為單道采樣點。矩陣Q的具體形式為

為了將重加權矩陣元素大小與反射系數建立聯系，利用邊界反射振幅信息將qi定義為

式中ξ為一個很小的常量，本文稱為加權因子。qi與反射系數的絕對值呈負相關。當上、下兩層參數變化大時，即為巖層邊界時，反射系數較大，此時權重減小，使稀疏約束項的影響減小，故可以在巖層邊界處產生較大反射系數；而當上、下兩層參數變化小或不變時，即代表同一層內，此時反射系數小或無，將導致權重增大，間接增大了稀疏約束的比重。通過這種動態調節約束項權重，本文方法可以得到比傳統方法(不重加權)更稀疏的反演結果。

參照張雨強等[15]對比范數稀疏性的方法，將本文重加權L1范數和傳統L1范數的稀疏問題分別表示為

式中：H為預測值；HR是實際值；η為極小正數。

由圖1 可見，交點位置即為問題的最優解，圖1b 交點(藍點)位置比圖1a 更趨近坐標縱軸，這說明重加權L1范數比傳統L1范數更具稀疏性，可獲得更稀疏的反演結果。因此，本文構建的反演目標函數為

圖1 L1范數（a）與重加權L1范數（b）的稀疏性對比

1.3 反演算法

對于傳統方法的目標函數(式(17))和本文方法的目標函數(式(23))，本文采用交替方向乘子方法(ADMM)將含多個未知參數的目標函數分解成多個單一參數的子問題進行求解。

以本文方法為例，首先采用拉格朗日乘子P代替L1范數中的QBDM矩陣進行運算，將上述非線性問題轉化為有約束的線性問題，得到目標函數

進一步引入二次懲罰項和對偶項[28]，將式(24)有約束的目標函數轉化為無約束的目標函數，得到

于是，將式(25)拆解成分別關于M、P和C的三個子優化問題，即

①與M相關的目標函數為

②與P相關的目標函數為

③與C相關的目標函數為

針對式(26)采用函數梯度方法求解，可得

式中：“·”為標量乘法運算符號；k為迭代次數；I為單位矩陣。

針對式(27)，引入軟閾值收縮算法求解同時含有L1范數和L2范數的目標函數，可得

式中sign(·)為符號函數，其定義為

針對式(28)，可根據式(26)相同的求解方法，得到

基于式(20)，本文引入的重加權矩陣對角線元素的迭代表達式為

通過對目標函數的詳細求解，本文方法的地震反演流程框架如表1所示。

表1 疊前重加權L1范數的反演流程

2 模型測試

選擇經典Marmousi Ⅱ模型測試本文方法和傳統方法。該模型為橫向500道、縱向500個采樣點，其縱波速度、橫波速度和密度的剖面分別如圖2所示。該模型中包含多種不同厚度、不同形狀的砂巖地層，其中多組薄層將作為驗證本文算法的重要地層。該模型還發育邊界清晰的斷層和背斜構造，適合驗證本文方法在刻畫地層邊界、減弱偽層現象方面的優越性。

圖2 不同參數的理論模型

2.1 無噪測試

首先基于式(4)，利用理論模型數據計算入射角為10°、20°和30°的地層反射系數；然后分別與頻率為30 Hz的雷克子波褶積合成對應的近、中、遠角度的地震記錄(圖3)。

圖3 無噪情況下不同角度數據的合成地震記錄

對反演目標函數添加可以補充背景信息的低頻模型，對理論模型的縱波速度、橫波速度和密度數據進行高斯低通濾波，即可獲得疊前縱波速度、橫波速度和密度的低頻模型(圖4)。

圖4 不同參數低頻初始模型

由式(23)可知，正則化參數α用于控制稀疏約束項的權重，α越大，稀疏約束的力度越大；參數λ用于控制低頻模型約束的權重，λ越大，低頻模型補充的先驗背景信息的影響越大。在反演迭代過程中，低頻模型約束項的主要作用是保證迭代的穩定性，而不對反演效果起主導作用，因此本文引入L曲線[29]確定本文方法和傳統方法的最優正則化參數α。經過大量啟發式參數測試，在無噪情況下，傳統方法的最優參數為λ1=5.1×10－4、λ2=5.0×10－4、λ3=5.0×10－3、α=1.5×10－8、μ=1.5×10－8；本文方法的最優參數為λ1=5.1×10－6、λ2=5.0×10－6、λ3=9.0×10－5、α=8.5×10－10、μ=1.5×10－9、ξ=1.0×10－5。

分別采用傳統方法與本文方法對縱波速度、橫波速度和密度進行同時反演。對比圖5 與圖6 可知，本文方法反演的地層邊界更清晰，縱向分辨率明顯更高(黑色箭頭處)，尤其針對薄氣砂巖儲層刻畫能力更強(黑色虛線圓圈)。

圖5 無噪情況下傳統方法不同參數的反演結果

圖6 無噪情況下本文方法不同參數的反演結果

為了更準確地分析兩種方法的反演效果，圖7 展示了無噪情況下第100 道(單道)不同參數反演結果。由圖可見，在傳統方法的反演結果中，同一地層出現較強的“正弦形”擾動(紅色箭頭處)，導致同層中產生許多偽速度層和偽密度層(小波浪狀)，反演結果不夠準確；而在本文方法反演結果中，邊界更清晰，很大程度地減弱了偽層現象，保證了同層反演的穩定性。傳統方法反演整體呈脈沖式(黑色箭頭處)，局部平滑了相鄰地層邊界，這是降低縱向分辨率的主要原因，而本文方法的反演結果與理論模型吻合度更高。傳統方法識別薄層的分辨率低、精度低(藍色箭頭處)，而本文方法通過加權反射系數的稀疏性，對薄層的識別能力明顯提高，同時也在一定程度上降低了密度反演結果的不穩定性。

2.2 抗噪性測試

為測試本文方法的抗噪性，對遠、中、近角度的地震記錄添加20%的高斯白噪聲(圖8)，再分別采用傳統方法和本文方法反演。在調參數過程中，隨著加噪程度的增大，建議適當增大加權因子ξ。同樣，得出加噪情況下傳統方法的最優參數為λ1=1.0×10－1、λ2=5.5×10－2、λ3=7.0×10－1、α=5.7×10－1、μ=2.0×10－4；本文方法的最優參數為λ1=1.9×10－2、λ2=1.5×10－2、λ3=2.4×10－1、α=1.5×10－1、μ=2.9×10－4、ξ=1.0×10－5。

對比圖9與圖10可知，本文方法反演的地層邊界分辨率更高(黑色箭頭處)；薄氣砂巖儲層邊界更清晰，偽層現象大幅減少(黑色虛線圈處)。由此證明，在加噪情況下本文方法仍然能夠獲得比傳統方法分辨率更高的縱波速度、橫波速度和密度剖面。

圖9 添加20%的高斯白噪聲的傳統方法不同參數反演結果

圖10 添加20%的高斯白噪聲的本文方法不同參數反演結果

圖11展示了加噪情況下第100道(單道)反演效果對比。由圖可見，相比傳統方法，本文方法反演的縱波速度、橫波速度和密度邊界更清晰，且對薄層反演精度更高，識別更準確(黑色箭頭處)。

圖11 添加20%的高斯白噪聲的第100 道不同參數反演結果

為了進一步測試本文方法的抗噪性，對地震記錄加50%的高斯白噪聲，得到的地震記錄如圖12所示。分別采用本文方法與傳統方法反演，通過啟發式確定傳統方法的最優參數為λ1=1.0×10－1、λ2=5.5×10－2、λ3=7×10－1、α=5.7×10－1、μ=7.0×10－1；本文方法的最優參數為λ1=1.0×10－2、λ2=1.4×10－2、λ3=1×10－1、α=8.7×10－1、μ=5.9×10－4、ξ=1.0×10－2。

傳統方法和本文方法得到的反演結果分別如圖13 和圖14所示。由圖可見，本文方法對細薄氣砂巖儲層的形態刻畫更清晰，識別能力更強(黑色虛線圈內)；同時，本文方法對地層邊界刻畫更清晰，而傳統方法模糊了地層邊界。

圖13 添加50%的高斯白噪聲的傳統方法不同參數反演結果

圖14 添加50%的高斯白噪聲的本文方法不同參數反演結果

通過對比第100 道(單道)不同參數反演結果(圖15)，本文方法對邊界刻畫能力更強(黑色箭頭處)，反演結果更準確(藍色箭頭處)。

圖15 添加50%的高斯白噪聲的第100 道不同參數反演結果

綜上所述，通過添加不同的噪聲，本文方法仍然能夠獲得比較合理、可靠的反演結果。

2.3 魯棒性測試

為測試本文算法的魯棒性，對地震記錄添加大小為地震振幅一半的異常值，異常值數量為地震信號長度的3%。此時，傳統方法的最優參數為λ1=8.0×10－2、λ2=5.0×10－2、λ3=4.0×10－1、α=1.5×10－6、μ=1.5×10－4；本文方法的最優參數為λ1=1.0×10－2、λ2=1.0×10－2、λ3=9.0×10－1、α=9.5×10－4、μ=1.0×10－4、ξ=1.0×10－5。

分別采用本文方法和傳統方法對第100道進行單道反演，結果如圖16所示。由圖可見，兩種方法均受到異常值的較大影響，傳統方法產生大量不準確的脈沖解，偽層現象和非聚焦平滑過渡現象加重(黑色箭頭)；而本文方法表現出更強的穩定性，邊界刻畫更準確，這說明本文方法魯棒性更強，反演結果更加穩定、可靠。

2.4 反演效果評價

引入信噪比(SNR)和歸一化均方根誤差(NRMSE)[30]用于評價反演效果。與均方根誤差(RMSE)相比，NRMSE 更能反映出反演結果與理論模型之間的平均差異(或偏差)和差異的可變性，以從橫向和縱向的角度更直觀、準確地對比兩種方法的反演效果。它們表達式分別為

式中：Y為理論參數值；為理論參數的均值；X為反演得到的參數；Ymax為理論參數的最大值；Ymin為理論參數的最小值；Tr為剖面總道數。

由表2可知，本文方法反演的縱波速度、橫波速度和密度的信噪比均高于傳統方法；從縱向分析，本文方法的歸一化均方根誤差均低于傳統方法；從橫向分析，本文方法與傳統方法獲得的縱波速度、橫波速度和密度中，縱波速度、橫波速度的歸一化均方根誤差均小于密度，說明縱波速度、橫波速度與模型差異的可變性更小，反演能力更強，反演結果質量更高。

表2 無噪情況下不同方法反演結果評價

表3和表4分別為加噪20%和50%的評價結果。與表2 相似，均證明了本文方法比傳統方法反演精度更高，抗噪性更強。

表3 加噪20%情況下不同方法反演結果評價

表4 加噪50%情況下不同方法反演結果的定量評價

3 實際資料應用

實際資料來源于中國西部地區，目的層儲層主要為頁巖。選擇A 井(直井)測井數據參與實際應用。

首先，對1.770 s 處目的層按優勢響應角度篩選出三組中心入射角的地震數據(圖17)，并用于縱波速度、橫波速度和密度的疊前三參數反演，它們依次為3°(1°～5°)、15°(13°～17°)、27°(25°～29°)。然后，對A井的縱波速度、橫波速度和密度井數據進行內插、外推和低通濾波，構建出各參數的低頻模型。再分別采用傳統方法和本文方法進行反演，此時傳統方法的最優參數為λ1=4.0×109、λ2=1.5×1010、λ3=5.0×1011、α=8.7×1013、μ=2.5×1010；本文方法的最優參數為λ1=4.0×109、λ2=5.0×1010、λ3=8.0×1011、α=8.7×1012、μ=2.5×105、ξ=5.0×10－2。

圖17 過A 井的共角度疊加地震剖面

圖18為過A 井的單道反演結果，可見本文方法比傳統方法與實際數據吻合得更好。

圖18 過A 井的單道不同方法、不同參數反演結果

同樣，由圖19 可見，本文方法反演結果與實際數據吻合較好(黑色箭頭處)。對于縱波速度反演結果來說，傳統方法受偽層干擾，導致地層分辨率低(藍色虛線圈內)，出現混層現象；本文方法對薄層識別能力明顯增強(紅色箭頭處)，地層的橫向連續性更好(黑色虛線圈內)，對地層邊界刻畫更清晰(藍色虛線圈內)，縱向分辨率明顯提高。類似地，橫波速度和密度反演結果也能證明以上結論，這表明本文方法在實際應用中的可行性和優越性。

4 結論

本文提出的疊前重加權L1范數稀疏約束的地震反演方法，充分利用了疊前地震記錄的邊界振幅信息。通過理論模型測試和實際工區應用，與傳統基于L1范數稀疏約束的疊前反演方法對比，得出如下結論：

(1)本文方法能產生更多稀疏解，提高了反演結果的稀疏性，有利于刻畫地層邊界，彌補了L1范數在地層邊界非聚焦平滑過渡的缺陷；

(2)本文方法加權利用了地震反射系數的稀疏性，反演結果分辨率更高，減弱了L1范數約束反演中嚴重的偽層現象，使疊前地震反演結果更加準確；

(3)本文方法具有較強的抗噪性，可為后續參數預測奠定更準確的基礎。

由于本文方法采用的是單道反演形式，后續可以研究如何實現疊前重加權L1范數稀疏約束的多道同時反演，以進一步提升預測地層的橫向連續性。