楊 曉 祝厚勤 王彥飛*

(①中國(guó)科學(xué)院油氣資源研究重點(diǎn)實(shí)驗(yàn)室，中國(guó)科學(xué)院地質(zhì)與地球物理研究所，北京100029； ②中國(guó)科學(xué)院大學(xué)，北京100049； ③中國(guó)科學(xué)院地球科學(xué)研究院，北京100029； ④中國(guó)石油勘探開發(fā)研究院，北京100083)

1 引言

頁(yè)巖氣作為一種清潔、高效的非常規(guī)能源，其開發(fā)和利用對(duì)發(fā)展低碳經(jīng)濟(jì)、優(yōu)化能源結(jié)構(gòu)具有重要意義[1-3]。巖石脆性作為遴選高品質(zhì)頁(yè)巖的重要評(píng)價(jià)指標(biāo)，是頁(yè)巖氣儲(chǔ)層描述的主要內(nèi)容之一[4]，而脆性指數(shù)則是巖石脆性的主要評(píng)價(jià)指標(biāo)。通常，脆性指數(shù)的選取是基于模型正演敏感性分析，利用縱、橫波阻抗、拉梅系數(shù)、楊氏模量、泊松比等彈性參數(shù)，獲得不同方法的脆性指數(shù)結(jié)果，最終優(yōu)選出敏感的脆性指數(shù)[5-8]。

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

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

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

為克服彈性阻抗反演問(wèn)題的不適定性，本文利用Tikhonov正則化方法和L-曲線準(zhǔn)則，將解的L2范數(shù)最小作為約束條件加入原問(wèn)題并建立最優(yōu)化模型。在數(shù)值實(shí)現(xiàn)中，利用雙濾波因子技巧，即Tikhonov濾波因子與截?cái)嗥娈愔甸撝禐V波因子相結(jié)合的方式，獲得較好的彈性阻抗反演結(jié)果。

2 雙濾波因子正則化

2.1 雙濾波因子正則化

實(shí)測(cè)地震記錄可視為地層反射系數(shù)與震源子波的卷積[26]，即

(1)

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

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

d=WR+N

(2)

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

(3)

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

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

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

(4)

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

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的奇異系統(tǒng)可記為{σi,ui,vj}?；贛oore-Penrose廣義逆，根據(jù)Picard定理[18,22]可將最優(yōu)化問(wèn)題(式(4))的極小解寫成

(7)

(8)

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

(1)Tikhonov濾波因子

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

(9)

由于矩陣W的病態(tài)性，不適定問(wèn)題(式(2))的解可以利用W的數(shù)值秩rδ表示。設(shè)Bδ為W的擾動(dòng)矩陣，則數(shù)值秩rδ滿足

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

(10)

(2)閾值濾波因子

根據(jù)以上討論，式(8)的濾波函數(shù)fi(α)也可以取截?cái)嗥娈愔捣纸獾男问?/p>

(11)

在數(shù)值實(shí)現(xiàn)中，把上述Tikhonov濾波因子與閾值濾波因子結(jié)合起來(lái)，形成雙濾波因子正則化。由于引起彈性阻抗反演問(wèn)題求解不穩(wěn)定的因素主要是正演算子小奇異值以及地震數(shù)據(jù)中的噪聲干擾，因此，在確定最優(yōu)正則參數(shù)α?xí)r，借助閾值濾波因子把小奇異值濾掉，以提高算法的穩(wěn)定性。然后，再利用Tikhonov濾波因子，求解最優(yōu)化模型(式(4))的極小解。

2.2 后驗(yàn)選取正則參數(shù)

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

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

(12)

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

下面給出利用Tikhonov方法和截?cái)嗥娈愔捣椒ㄏ嘟Y(jié)合，以及L-曲線準(zhǔn)則求取彈性阻抗反演問(wèn)題的雙濾子正則化算法。

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

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

(13)

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

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

(14)

以及殘差的范數(shù)

(15)

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

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

(16)

3 數(shù)值試驗(yàn)

3.1 理論模型正演

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

圖1 彈性阻抗模型

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

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

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

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

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

(17)

式中：fm=30Hz是子波主頻;k是時(shí)間采樣間隔，k=1ms。子波樣點(diǎn)個(gè)數(shù)是101。

由圖3、圖4的模型正演剖面可見，頁(yè)巖和含氣砂巖之間的反射同相軸為負(fù)極性、弱反射；當(dāng)數(shù)據(jù)中存在噪聲時(shí)，難以在地震剖面上準(zhǔn)確定位頁(yè)巖和含氣砂巖地層的分界面。

3.2 理論模型反演

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

圖5 理論模型的L-曲線

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

圖6 彈性阻抗反演剖面

圖7 誤差剖面

3.3 實(shí)際數(shù)據(jù)反演

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

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

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

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

圖9 頁(yè)巖與灰?guī)r彈性阻抗特征分析

圖10 實(shí)際資料的L-曲線

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

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

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

Brit=Eρ

(18)

式中：E是楊氏模量；ρ是密度。如圖12所示，在頁(yè)巖層段內(nèi)部，脆性指數(shù)較高的區(qū)域，即黃色到紅色區(qū)域，為頁(yè)巖甜點(diǎn)區(qū)，內(nèi)部細(xì)節(jié)較為豐富，具有很好的指示作用。根據(jù)反演結(jié)果與測(cè)線附近井資料的對(duì)比分析，目的層段彈性阻抗均值與測(cè)井?dāng)?shù)據(jù)的誤差約為11%，脆性指數(shù)的分布也與鉆井基本吻合。

圖12 基于彈性反演結(jié)果計(jì)算的敏感脆性指數(shù)剖面

4 結(jié)論

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

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

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

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