彈性阻抗反演的后驗正則化方法

楊 曉 祝厚勤 王彥飛*

(①中國科學院油氣資源研究重點實驗室，中國科學院地質與地球物理研究所，北京100029； ②中國科學院大學，北京100049； ③中國科學院地球科學研究院，北京100029； ④中國石油勘探開發研究院，北京100083)

1 引言

頁巖氣作為一種清潔、高效的非常規能源，其開發和利用對發展低碳經濟、優化能源結構具有重要意義[1-3]。巖石脆性作為遴選高品質頁巖的重要評價指標，是頁巖氣儲層描述的主要內容之一[4]，而脆性指數則是巖石脆性的主要評價指標。通常，脆性指數的選取是基于模型正演敏感性分析，利用縱、橫波阻抗、拉梅系數、楊氏模量、泊松比等彈性參數，獲得不同方法的脆性指數結果，最終優選出敏感的脆性指數[5-8]。

目前，疊后聲波阻抗反演、AVO反演和彈性阻抗反演是提取地層彈性參數的主要方法。疊后聲波阻抗反演方法結果單一，未充分利用疊前地震信息[9]。AVO反演方法對數據信噪比要求較高，數據信噪比低可能導致地震子波隨炮檢距變化，必然會降低AVO反演結果的準確性[10]。與疊后聲波阻抗反演和AVO反演相比，Connolly[11]提出的彈性阻抗反演方法，考慮了AVO效應，利用部分疊加提高數據信噪比，能對近、遠炮檢距資料進行有效標定并獲得更為豐富、穩定、可靠的彈性參數反演結果[12-14]，有利于頁巖敏感性脆性指數的選取。

彈性阻抗反演是通過求解與離散線性算子方程相關的最小二乘問題實現。彈性阻抗反演是Ha-damard意義下的不適定問題[15]，即不能同時滿足解的存在性、唯一性和穩定性的要求。因此直接求解離散線性算子方程相關的最小二乘問題會帶來不穩定的計算結果[16,17]，而正則化方法是求解不適定問題的有效方法[15]。正則化方法是用一簇與原問題相“鄰近”的適定問題的解去逼近原問題的真實解，即通過構造正則算子得到原問題穩定的近似解。其中，正則化參數控制與原問題的“鄰近程度”。正則化方法有Tikhonov方法、Landweber迭代法、Newton型方法、共軛梯度法等[18]，其中Tikhonov正則化方法應用最為廣泛。

最優正則化參數的選取是正則化方法的核心問題之一。正則化參數取值過大會導致借助構造正則算子所獲得的解與原問題偏差過大，而取值太小會過多地“繼承”原問題的不適定性而導致解的不穩定。正則化參數的選取方法通常分為先驗和后驗兩大類[18]。選取正則參數的先驗方法雖具有理論分析意義，但是在實際中常常難以實現。后驗方法有如下幾種形式：偏差原理[19]、廣義偏差原理[20]、誤差極小化準則[21]、Arcangeli準則[22]等，這些方法需要預先估計原始數據的誤差水平。針對原始資料誤差水平未知的情況，Tikhonov等[23]提出擬最優準則、Hanke等[24]提出L-曲線準則、Golub等[25]提出廣義交叉校驗準則。其中，利用L-曲線準則確定最優化正則參數的方法在圖像恢復、生命科學、遙感技術等領域已經得到廣泛應用[15]。

為克服彈性阻抗反演問題的不適定性，本文利用Tikhonov正則化方法和L-曲線準則，將解的L2范數最小作為約束條件加入原問題并建立最優化模型。在數值實現中，利用雙濾波因子技巧，即Tikhonov濾波因子與截斷奇異值閾值濾波因子相結合的方式，獲得較好的彈性阻抗反演結果。

2 雙濾波因子正則化

2.1 雙濾波因子正則化

實測地震記錄可視為地層反射系數與震源子波的卷積[26]，即

(1)

式中：ti表示時間延遲；Ri為地層反射系數；w為震源子波；d(t)表示地震信號；N(t)表示干擾噪聲。

式(3)可以寫成矩陣方程的形式

d=WR+N

(2)

式中：d=(d1,d2,…,dm)T;R=(R1,R2,…,Rn)T;N=(N1,N2,…,Nm)T以及

(3)

W是將模型空間投影到觀測空間的算子。式(2)的線性反問題的求解是彈性阻抗反演的關鍵。

Tikhonov正則化方法根據Lagrange乘子法或優化理論中構造泛函的思想，通過引入光滑泛函構造正則算子，進而求解不適定性反問題[18]。本文利用Tikhonov正則化思想，將式(2)的線性反問題的求解簡化為數據擬合問題，即‖WR-d‖2→min的形式，并建立基于L2范數約束的最優化模型

Jα(R)=‖WRα-d‖22+α‖K(Rα-R0)‖22→min

(4)

式中：K是規范化矩陣，通常是一個(半)正定矩陣，本文令K為單位矩陣；R0表示初始解(提供解空間的先驗信息)，本文令R0=0；α>0是正則化參數，平衡殘差項‖WRα-d‖22與附加約束項‖K(Rα-R0)‖22，以避免過正則化或欠正則化。

W=UΣV

(5)

式中：U、V為酉矩陣，U=(u1,u2,…,um)，V=(v1,v2,…,vn)；Σ為半正定矩陣，形式如下

(6)

其中，Σr=diag(σ1,σ2,…,σr)，σ1≥σ2≥…≥σr，σi是W的奇異值。矩陣算子W的奇異系統可記為{σi,ui,vj}。基于Moore-Penrose廣義逆，根據Picard定理[18,22]可將最優化問題(式(4))的極小解寫成

(7)

(8)

式中fi(α)是濾波因子。為了獲得原問題的穩定近似解，需要選擇合適的濾波因子。本文針對彈性阻抗反演問題考慮兩種濾波因子。

(1)Tikhonov濾波因子

在K是單位矩陣的情況下，Tikhonov濾波因子為

(9)

由于矩陣W的病態性，不適定問題(式(2))的解可以利用W的數值秩rδ表示。設Bδ為W的擾動矩陣，則數值秩rδ滿足

rδ=min{rank(Bδ+A)∶Bδ∈Rm×n，‖Bδ‖2≤δ}

(10)

(2)閾值濾波因子

根據以上討論，式(8)的濾波函數fi(α)也可以取截斷奇異值分解的形式

(11)

在數值實現中，把上述Tikhonov濾波因子與閾值濾波因子結合起來，形成雙濾波因子正則化。由于引起彈性阻抗反演問題求解不穩定的因素主要是正演算子小奇異值以及地震數據中的噪聲干擾，因此，在確定最優正則參數α時，借助閾值濾波因子把小奇異值濾掉，以提高算法的穩定性。然后，再利用Tikhonov濾波因子，求解最優化模型(式(4))的極小解。

2.2 后驗選取正則參數

在正則化方法中，最優的正則參數α起到平衡殘差項和有關解約束條件的作用。Hansen等[27,28]正是基于這種考慮，提出L-曲線準則，即在對數坐標系統下借助殘差項‖WRα-d‖2和約束條件‖K(Rα-R0)‖2隨正則參數α變化的曲線確定最優正則參數的方法。該方法因曲線在對數尺度中會呈現“L”形狀而得名，其關鍵參數是L-曲線隅角。Hanke等[24]定義L-曲線隅角為其在對數尺度下的最大曲率。

令ρ(α)=lg‖WRα-d‖2,θ(α)=lg‖Rα‖2，定義一個有關正則參數α的非線性函數

(12)

式中：ρ′、ρ″和θ′、θ″分別為函數ρ(α)和θ(α)的一階和二階導數；φ(α)稱作參數α的曲率函數。

下面給出利用Tikhonov方法和截斷奇異值方法相結合，以及L-曲線準則求取彈性阻抗反演問題的雙濾子正則化算法。

(1)對算子矩陣W作奇異值分解

[U,Σ,V]=svd(W)

(13)

式中Σ=diag(σ1,σ2,…，σn)是矩陣奇異值，呈降序排列。

(2)基于σr≤α≤σ1(r≤n)，r是W的數值秩，借助閾值濾波因子，計算反問題逼近解的范數

(14)

以及殘差的范數

(15)

(3)以圖形形式(lg‖Rα‖2, lg‖WRα-d‖2)繪制L-曲線、確定曲線隅角，隅角處對應的正則參數α就是最優正則參數。

(4)借助Tikhonov濾波因子，獲得反問題的正則解

(16)

3 數值試驗

3.1 理論模型正演

首先，利用Castagna等[29]提出的第三類經典AVO模型進行試算。如圖1所示，構建一個三層楔狀構造模型，時間范圍為0～250ms，共140個地震道，地層巖性從上至下分別為頁巖、含氣砂巖、含水砂巖，各地層的彈性參數如表1所示。圖2是地層彈性阻抗隨入射角的變化曲線。由圖可見，頁巖層和含水砂巖層的彈性阻抗隨入射角而增大，含氣砂巖層的彈性阻抗值隨入射角而減小。也就是說，當入射角增大時，各反射界面上、下地層的彈性阻抗差異變大。

圖1 彈性阻抗模型

圖2 地層的彈性阻抗隨入射角的變化曲線

巖性縱波速度m/s橫波速度m/s密度g/cm3頁巖18304002.02含氣砂巖14405801.53含水砂巖21306701.90

圖3 正演地震剖面(無噪聲)

圖4 加噪正演地震剖面(S/N=5)

為便于分析對比，本文主要研究入射角為30°的情況。圖2是彈性阻抗模型的正演剖面。正演過程中采用零相位雷克子波

(17)

式中：fm=30Hz是子波主頻;k是時間采樣間隔，k=1ms。子波樣點個數是101。

由圖3、圖4的模型正演剖面可見，頁巖和含氣砂巖之間的反射同相軸為負極性、弱反射；當數據中存在噪聲時，難以在地震剖面上準確定位頁巖和含氣砂巖地層的分界面。

3.2 理論模型反演

利用雙濾波因子正則化算法對圖4的地震數據進行反演。圖5是該反演過程的L-曲線，橫、縱坐標分別是對數尺度下反問題逼近解的范數‖Rα‖2和殘差項的范數‖WRα-d‖2。根據Hanke等[24]的定義，利用式(8)確定L-曲線的隅角(紅色標記)以及它所對應的正則參數α=0.8846。

圖5 理論模型的L-曲線

圖6是加噪正演數據的彈性阻抗反演剖面。為了驗證雙濾波因子正則化算法的精度，圖7給出反演結果與真實模型的誤差剖面。從圖7可以看出，頁巖地層反演結果的平均誤差約為7%，而含氣砂巖、含水砂巖的反演結果誤差相對較高，分別約為10%和13%。

圖6 彈性阻抗反演剖面

圖7 誤差剖面

3.3 實際數據反演

利用某油田的實際三維地震資料對雙濾波因子正則化算法進行測試。本次反演目的層為泥盆系頁巖氣層，構造較為平緩，頂面埋深為2500～4500m，沉積于開闊海域環境，最大海進盆地快速充填期，上覆和下伏地層均為泥灰巖。圖8所示為L80主測線入射角為25°的地震剖面，可以看出，目的層為一低頻弱振幅波谷反射，反射特征清楚，橫向連續性好。

通過巖石物理分析，頁巖與灰巖阻抗差異明顯，并且其差異隨入射角度的增大而增大。如圖9所示，在25°入射時，頁巖的彈性阻抗范圍是9×103～1.1×104g/cm3·m/s，泥灰巖彈性阻抗范圍是1.1×104～1.7×104g/cm3·m/s，說明彈性阻抗能有效地區分頁巖與泥巖。

圖10是實際資料的L-曲線，描述實際資料反演過程中的殘差項和約束條件隨正則參數的變化情況。圖中的橫、縱坐標分別是對數尺度下反問題逼近解的范數‖Rα‖2 和殘差項的范數‖WRα-d‖2。根據Hanke等[24]的定義，利用式(8)確定L-曲線隅角(曲線紅色標記處)及其對應的正則參數α=0.5469。

圖8 L80主測線25°入射角的地震剖面

圖9 頁巖與灰巖彈性阻抗特征分析

圖10 實際資料的L-曲線

圖11是彈性阻抗反演剖面。由圖可見，借助于彈性阻抗剖面可以容易地區分頁巖與灰巖。目的層頁巖是剖面1700～1750ms范圍的深藍色部分，表現為明顯的低阻抗特征，具備較好的追蹤解釋條件，厚度為15～20m。剖面上1650～1700ms的淺藍色區域解釋為泥灰巖，其彈性阻抗范圍為1.15×104～1.3×104g/cm3·m/s。圖11的彈性阻抗反演結果與圖9所示的該地區巖石物理分析結果相吻合。

圖11 本文方法反演的彈性阻抗剖面

基于彈性阻抗反演結果可獲得縱、橫波阻抗、楊氏模量、泊松比等彈性參數，進而求取脆性指數[5-8]。圖12是根據下式計算的數據敏感脆性指數

Brit=Eρ

(18)

式中：E是楊氏模量；ρ是密度。如圖12所示，在頁巖層段內部，脆性指數較高的區域，即黃色到紅色區域，為頁巖甜點區，內部細節較為豐富，具有很好的指示作用。根據反演結果與測線附近井資料的對比分析，目的層段彈性阻抗均值與測井數據的誤差約為11%，脆性指數的分布也與鉆井基本吻合。

圖12 基于彈性反演結果計算的敏感脆性指數剖面

4 結論

為克服彈性阻抗反演問題的不適定性，本文提出基于Tikhonov濾波因子和閾值濾波因子相結合的雙濾波因子正則化算法，輔以后驗最優正則化參數選取方法(L-曲線準則)，獲得了較好的彈性阻抗反演結果，取得如下幾點認識。

(1)建立基于L2范數約束的彈性阻抗反演最優化模型。

(2)提出先借助閾值濾波因子和L-曲線準則確定最優正則化參數，然后利用Tikhonov濾波因子求解最優化模型的極小解的方法，能有效地提高計算精度以及解的穩定性。

(3)通過低信噪比理論模型和實際地震資料對算法進行測試。低信噪比理論模型的反演結果與模型正演情況基本吻合；基于實際地震資料的彈性阻抗反演結果與該地區巖石物理分析結果也基本符合，進一步獲得的脆性指數分布的預測結果也與鉆井吻合。這表明算法具有良好的穩定性和抗噪性，具有推廣應用前景。