楊 曉 祝厚勤 王彥飛*
(①中國科學(xué)院油氣資源研究重點實驗室,中國科學(xué)院地質(zhì)與地球物理研究所,北京100029; ②中國科學(xué)院大學(xué),北京100049; ③中國科學(xué)院地球科學(xué)研究院,北京100029; ④中國石油勘探開發(fā)研究院,北京100083)
頁巖氣作為一種清潔、高效的非常規(guī)能源,其開發(fā)和利用對發(fā)展低碳經(jīng)濟、優(yōu)化能源結(jié)構(gòu)具有重要意義[1-3]。巖石脆性作為遴選高品質(zhì)頁巖的重要評價指標(biāo),是頁巖氣儲層描述的主要內(nèi)容之一[4],而脆性指數(shù)則是巖石脆性的主要評價指標(biāo)。通常,脆性指數(shù)的選取是基于模型正演敏感性分析,利用縱、橫波阻抗、拉梅系數(shù)、楊氏模量、泊松比等彈性參數(shù),獲得不同方法的脆性指數(shù)結(jié)果,最終優(yōu)選出敏感的脆性指數(shù)[5-8]。
目前,疊后聲波阻抗反演、AVO反演和彈性阻抗反演是提取地層彈性參數(shù)的主要方法。疊后聲波阻抗反演方法結(jié)果單一,未充分利用疊前地震信息[9]。AVO反演方法對數(shù)據(jù)信噪比要求較高,數(shù)據(jù)信噪比低可能導(dǎo)致地震子波隨炮檢距變化,必然會降低AVO反演結(jié)果的準(zhǔn)確性[10]。與疊后聲波阻抗反演和AVO反演相比,Connolly[11]提出的彈性阻抗反演方法,考慮了AVO效應(yīng),利用部分疊加提高數(shù)據(jù)信噪比,能對近、遠(yuǎn)炮檢距資料進行有效標(biāo)定并獲得更為豐富、穩(wěn)定、可靠的彈性參數(shù)反演結(jié)果[12-14],有利于頁巖敏感性脆性指數(shù)的選取。
彈性阻抗反演是通過求解與離散線性算子方程相關(guān)的最小二乘問題實現(xiàn)。彈性阻抗反演是Ha-damard意義下的不適定問題[15],即不能同時滿足解的存在性、唯一性和穩(wěn)定性的要求。因此直接求解離散線性算子方程相關(guān)的最小二乘問題會帶來不穩(wěn)定的計算結(jié)果[16,17],而正則化方法是求解不適定問題的有效方法[15]。正則化方法是用一簇與原問題相“鄰近”的適定問題的解去逼近原問題的真實解,即通過構(gòu)造正則算子得到原問題穩(wěn)定的近似解。其中,正則化參數(shù)控制與原問題的“鄰近程度”。正則化方法有Tikhonov方法、Landweber迭代法、Newton型方法、共軛梯度法等[18],其中Tikhonov正則化方法應(yīng)用最為廣泛。
最優(yōu)正則化參數(shù)的選取是正則化方法的核心問題之一。正則化參數(shù)取值過大會導(dǎo)致借助構(gòu)造正則算子所獲得的解與原問題偏差過大,而取值太小會過多地“繼承”原問題的不適定性而導(dǎo)致解的不穩(wěn)定。正則化參數(shù)的選取方法通常分為先驗和后驗兩大類[18]。選取正則參數(shù)的先驗方法雖具有理論分析意義,但是在實際中常常難以實現(xiàn)。后驗方法有如下幾種形式:偏差原理[19]、廣義偏差原理[20]、誤差極小化準(zhǔn)則[21]、Arcangeli準(zhǔn)則[22]等,這些方法需要預(yù)先估計原始數(shù)據(jù)的誤差水平。針對原始資料誤差水平未知的情況,Tikhonov等[23]提出擬最優(yōu)準(zhǔn)則、Hanke等[24]提出L-曲線準(zhǔn)則、Golub等[25]提出廣義交叉校驗準(zhǔn)則。其中,利用L-曲線準(zhǔn)則確定最優(yōu)化正則參數(shù)的方法在圖像恢復(fù)、生命科學(xué)、遙感技術(shù)等領(lǐng)域已經(jīng)得到廣泛應(yīng)用[15]。
為克服彈性阻抗反演問題的不適定性,本文利用Tikhonov正則化方法和L-曲線準(zhǔn)則,將解的L2范數(shù)最小作為約束條件加入原問題并建立最優(yōu)化模型。在數(shù)值實現(xiàn)中,利用雙濾波因子技巧,即Tikhonov濾波因子與截斷奇異值閾值濾波因子相結(jié)合的方式,獲得較好的彈性阻抗反演結(jié)果。
實測地震記錄可視為地層反射系數(shù)與震源子波的卷積[26],即
(1)
式中:ti表示時間延遲;Ri為地層反射系數(shù);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)的線性反問題的求解是彈性阻抗反演的關(guān)鍵。
Tikhonov正則化方法根據(jù)Lagrange乘子法或優(yōu)化理論中構(gòu)造泛函的思想,通過引入光滑泛函構(gòu)造正則算子,進而求解不適定性反問題[18]。本文利用Tikhonov正則化思想,將式(2)的線性反問題的求解簡化為數(shù)據(jù)擬合問題,即‖WR-d‖2→min的形式,并建立基于L2范數(shù)約束的最優(yōu)化模型
Jα(R)=‖WRα-d‖22+α‖K(Rα-R0)‖22→min
(4)
式中:K是規(guī)范化矩陣,通常是一個(半)正定矩陣,本文令K為單位矩陣;R0表示初始解(提供解空間的先驗信息),本文令R0=0;α>0是正則化參數(shù),平衡殘差項‖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的奇異系統(tǒng)可記為{σi,ui,vj}?;贛oore-Penrose廣義逆,根據(jù)Picard定理[18,22]可將最優(yōu)化問題(式(4))的極小解寫成
(7)
(8)
式中fi(α)是濾波因子。為了獲得原問題的穩(wěn)定近似解,需要選擇合適的濾波因子。本文針對彈性阻抗反演問題考慮兩種濾波因子。
(1)Tikhonov濾波因子
在K是單位矩陣的情況下,Tikhonov濾波因子為
(9)
由于矩陣W的病態(tài)性,不適定問題(式(2))的解可以利用W的數(shù)值秩rδ表示。設(shè)Bδ為W的擾動矩陣,則數(shù)值秩rδ滿足
rδ=min{rank(Bδ+A)∶Bδ∈Rm×n,‖Bδ‖2≤δ}
(10)
(2)閾值濾波因子
根據(jù)以上討論,式(8)的濾波函數(shù)fi(α)也可以取截斷奇異值分解的形式
(11)
在數(shù)值實現(xiàn)中,把上述Tikhonov濾波因子與閾值濾波因子結(jié)合起來,形成雙濾波因子正則化。由于引起彈性阻抗反演問題求解不穩(wěn)定的因素主要是正演算子小奇異值以及地震數(shù)據(jù)中的噪聲干擾,因此,在確定最優(yōu)正則參數(shù)α?xí)r,借助閾值濾波因子把小奇異值濾掉,以提高算法的穩(wěn)定性。然后,再利用Tikhonov濾波因子,求解最優(yōu)化模型(式(4))的極小解。
在正則化方法中,最優(yōu)的正則參數(shù)α起到平衡殘差項和有關(guān)解約束條件的作用。Hansen等[27,28]正是基于這種考慮,提出L-曲線準(zhǔn)則,即在對數(shù)坐標(biāo)系統(tǒng)下借助殘差項‖WRα-d‖2和約束條件‖K(Rα-R0)‖2隨正則參數(shù)α變化的曲線確定最優(yōu)正則參數(shù)的方法。該方法因曲線在對數(shù)尺度中會呈現(xiàn)“L”形狀而得名,其關(guān)鍵參數(shù)是L-曲線隅角。Hanke等[24]定義L-曲線隅角為其在對數(shù)尺度下的最大曲率。
令ρ(α)=lg‖WRα-d‖2,θ(α)=lg‖Rα‖2,定義一個有關(guān)正則參數(shù)α的非線性函數(shù)
(12)
式中:ρ′、ρ″和θ′、θ″分別為函數(shù)ρ(α)和θ(α)的一階和二階導(dǎo)數(shù);φ(α)稱作參數(shù)α的曲率函數(shù)。
下面給出利用Tikhonov方法和截斷奇異值方法相結(jié)合,以及L-曲線準(zhǔn)則求取彈性阻抗反演問題的雙濾子正則化算法。
(1)對算子矩陣W作奇異值分解
[U,Σ,V]=svd(W)
(13)
式中Σ=diag(σ1,σ2,…,σn)是矩陣奇異值,呈降序排列。
(2)基于σr≤α≤σ1(r≤n),r是W的數(shù)值秩,借助閾值濾波因子,計算反問題逼近解的范數(shù)
(14)
以及殘差的范數(shù)
(15)
(3)以圖形形式(lg‖Rα‖2, lg‖WRα-d‖2)繪制L-曲線、確定曲線隅角,隅角處對應(yīng)的正則參數(shù)α就是最優(yōu)正則參數(shù)。
(4)借助Tikhonov濾波因子,獲得反問題的正則解
(16)
首先,利用Castagna等[29]提出的第三類經(jīng)典AVO模型進行試算。如圖1所示,構(gòu)建一個三層楔狀構(gòu)造模型,時間范圍為0~250ms,共140個地震道,地層巖性從上至下分別為頁巖、含氣砂巖、含水砂巖,各地層的彈性參數(shù)如表1所示。圖2是地層彈性阻抗隨入射角的變化曲線。由圖可見,頁巖層和含水砂巖層的彈性阻抗隨入射角而增大,含氣砂巖層的彈性阻抗值隨入射角而減小。也就是說,當(dāng)入射角增大時,各反射界面上、下地層的彈性阻抗差異變大。
圖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。子波樣點個數(shù)是101。
由圖3、圖4的模型正演剖面可見,頁巖和含氣砂巖之間的反射同相軸為負(fù)極性、弱反射;當(dāng)數(shù)據(jù)中存在噪聲時,難以在地震剖面上準(zhǔn)確定位頁巖和含氣砂巖地層的分界面。
利用雙濾波因子正則化算法對圖4的地震數(shù)據(jù)進行反演。圖5是該反演過程的L-曲線,橫、縱坐標(biāo)分別是對數(shù)尺度下反問題逼近解的范數(shù)‖Rα‖2和殘差項的范數(shù)‖WRα-d‖2。根據(jù)Hanke等[24]的定義,利用式(8)確定L-曲線的隅角(紅色標(biāo)記)以及它所對應(yīng)的正則參數(shù)α=0.8846。
圖5 理論模型的L-曲線
圖6是加噪正演數(shù)據(jù)的彈性阻抗反演剖面。為了驗證雙濾波因子正則化算法的精度,圖7給出反演結(jié)果與真實模型的誤差剖面。從圖7可以看出,頁巖地層反演結(jié)果的平均誤差約為7%,而含氣砂巖、含水砂巖的反演結(jié)果誤差相對較高,分別約為10%和13%。
圖6 彈性阻抗反演剖面
圖7 誤差剖面
利用某油田的實際三維地震資料對雙濾波因子正則化算法進行測試。本次反演目的層為泥盆系頁巖氣層,構(gòu)造較為平緩,頂面埋深為2500~4500m,沉積于開闊海域環(huán)境,最大海進盆地快速充填期,上覆和下伏地層均為泥灰?guī)r。圖8所示為L80主測線入射角為25°的地震剖面,可以看出,目的層為一低頻弱振幅波谷反射,反射特征清楚,橫向連續(xù)性好。
通過巖石物理分析,頁巖與灰?guī)r阻抗差異明顯,并且其差異隨入射角度的增大而增大。如圖9所示,在25°入射時,頁巖的彈性阻抗范圍是9×103~1.1×104g/cm3·m/s,泥灰?guī)r彈性阻抗范圍是1.1×104~1.7×104g/cm3·m/s,說明彈性阻抗能有效地區(qū)分頁巖與泥巖。
圖10是實際資料的L-曲線,描述實際資料反演過程中的殘差項和約束條件隨正則參數(shù)的變化情況。圖中的橫、縱坐標(biāo)分別是對數(shù)尺度下反問題逼近解的范數(shù)‖Rα‖2 和殘差項的范數(shù)‖WRα-d‖2。根據(jù)Hanke等[24]的定義,利用式(8)確定L-曲線隅角(曲線紅色標(biāo)記處)及其對應(yīng)的正則參數(shù)α=0.5469。
圖8 L80主測線25°入射角的地震剖面
圖9 頁巖與灰?guī)r彈性阻抗特征分析
圖10 實際資料的L-曲線
圖11是彈性阻抗反演剖面。由圖可見,借助于彈性阻抗剖面可以容易地區(qū)分頁巖與灰?guī)r。目的層頁巖是剖面1700~1750ms范圍的深藍色部分,表現(xiàn)為明顯的低阻抗特征,具備較好的追蹤解釋條件,厚度為15~20m。剖面上1650~1700ms的淺藍色區(qū)域解釋為泥灰?guī)r,其彈性阻抗范圍為1.15×104~1.3×104g/cm3·m/s。圖11的彈性阻抗反演結(jié)果與圖9所示的該地區(qū)巖石物理分析結(jié)果相吻合。
圖11 本文方法反演的彈性阻抗剖面
基于彈性阻抗反演結(jié)果可獲得縱、橫波阻抗、楊氏模量、泊松比等彈性參數(shù),進而求取脆性指數(shù)[5-8]。圖12是根據(jù)下式計算的數(shù)據(jù)敏感脆性指數(shù)
Brit=Eρ
(18)
式中:E是楊氏模量;ρ是密度。如圖12所示,在頁巖層段內(nèi)部,脆性指數(shù)較高的區(qū)域,即黃色到紅色區(qū)域,為頁巖甜點區(qū),內(nèi)部細(xì)節(jié)較為豐富,具有很好的指示作用。根據(jù)反演結(jié)果與測線附近井資料的對比分析,目的層段彈性阻抗均值與測井?dāng)?shù)據(jù)的誤差約為11%,脆性指數(shù)的分布也與鉆井基本吻合。
圖12 基于彈性反演結(jié)果計算的敏感脆性指數(shù)剖面
為克服彈性阻抗反演問題的不適定性,本文提出基于Tikhonov濾波因子和閾值濾波因子相結(jié)合的雙濾波因子正則化算法,輔以后驗最優(yōu)正則化參數(shù)選取方法(L-曲線準(zhǔn)則),獲得了較好的彈性阻抗反演結(jié)果,取得如下幾點認(rèn)識。
(1)建立基于L2范數(shù)約束的彈性阻抗反演最優(yōu)化模型。
(2)提出先借助閾值濾波因子和L-曲線準(zhǔn)則確定最優(yōu)正則化參數(shù),然后利用Tikhonov濾波因子求解最優(yōu)化模型的極小解的方法,能有效地提高計算精度以及解的穩(wěn)定性。
(3)通過低信噪比理論模型和實際地震資料對算法進行測試。低信噪比理論模型的反演結(jié)果與模型正演情況基本吻合;基于實際地震資料的彈性阻抗反演結(jié)果與該地區(qū)巖石物理分析結(jié)果也基本符合,進一步獲得的脆性指數(shù)分布的預(yù)測結(jié)果也與鉆井吻合。這表明算法具有良好的穩(wěn)定性和抗噪性,具有推廣應(yīng)用前景。