基于疊前地震道集的流體因子同步反演方法

劉厚裕， 白俊雨

(1.中國石油化工股份有限公司 華東油氣分公司，南京 200019;2.中國石油化工股份有限公司 石油物探技術(shù)研究院，南京 211103)

0 引言

準(zhǔn)確識別儲層流體是高效地震勘探中非常重要的一環(huán)。疊前地震數(shù)據(jù)較疊后地震數(shù)據(jù)含有更為豐富的信息，從疊前地震數(shù)據(jù)中提取儲層巖性及流體參數(shù)是理論研究和實際應(yīng)用的熱點(diǎn)。其中AVO分析技術(shù)利用地震振幅隨偏移距變化關(guān)系，對地下介質(zhì)性質(zhì)進(jìn)行分析和預(yù)測，Ostrander[1]通過實際勘探和成功探井資料所證實AVO技術(shù)的應(yīng)用效果，之后該項技術(shù)在工業(yè)界取得極為廣泛的應(yīng)用；Russell等[2-4]基于巖石物理理論推導(dǎo)建立的Gassmann流體因子和Russell流體因子，在流體識別中比其他常規(guī)流體因子更加敏感，尤其在致密儲層中具有良好的識別能力；Gang Chen等[5]推導(dǎo)了縱波速度、流體因子、密度的縱波反射系數(shù)公式，通過兩步法實現(xiàn)了流體因子和密度的反演；Shengjun Li等[6]推導(dǎo)了泊松比、楊氏模量、Gassmann流體因子與縱波反射系數(shù)的關(guān)系，實現(xiàn)了泊松比、楊氏模量、Gassmann流體因子的直接反演方法。

杜炳毅等[7]以Russell近似理論為基礎(chǔ)，推導(dǎo)了轉(zhuǎn)換波AVO近似公式，通過PP波和PS波聯(lián)合反演獲取了流體因子、剪切模量和密度；高剛等[8]根據(jù)淺層砂巖儲層彈性特征設(shè)計了砂巖敏感識別因子SF4；蔡生娟等[9]針對非均質(zhì)儲層的流體識別問題，消除孔隙度對流體因子的影響，構(gòu)建孔隙度非敏感流體因子PINF；李超等[10]提出了一種在缺少遠(yuǎn)偏移距數(shù)據(jù)情況下精確估計流體的疊前地震反演方法，該方法同步反演Gassmann流體項和剪切模量；劉軍等[11]以珠江口盆地東部中深層碎屑巖儲層為例，基于巖石物理測試數(shù)據(jù)對比分析了流體因子對不同流體的敏感性；楊勤林等[12]對松遼盆地某區(qū)塊的致密砂巖儲層開展巖石物理分析和含氣性預(yù)測，結(jié)合測井和巖石物理測試數(shù)據(jù)分析表明含氣砂巖表現(xiàn)為低縱橫速度比和低流體因子f等特征，交會分析顯示流體因子f對含水和含氣砂體的區(qū)分能力相對最強(qiáng)；劉道理等[13]將流體因子f的頻散程度作為深層儲層含氣性預(yù)測的識別因子,實現(xiàn)了基于疊前地震資料的頻變流體參數(shù)反演。

對彈性參數(shù)進(jìn)行轉(zhuǎn)換而獲取流體因子的方法，存在累積誤差且不能充分發(fā)揮疊前地震資料振幅變化信息。對此筆者對流體因子的地震反演方法進(jìn)行研究，構(gòu)建了低頻約束的流體因子疊前地震道集同步反演目標(biāo)函數(shù)，并給出了梯度信息的解析解，實現(xiàn)了直接從疊前地震道集反演Gassmann流體因子、剪切模量和密度等三項參數(shù)反演。

1 基本原理

1.1 Gassmann流體因子反射系數(shù)近似方程

Russell[2-4]等基于巖石物理基礎(chǔ)理論，推導(dǎo)并提出了孔隙且飽含流體巖石的Russell流體因子(ρs和ρf)，對Aki-Richards方程進(jìn)行變換得到了f-μ-ρ反射系數(shù)方程[3]為式(1)。

(1)

(2)

將式(1)改寫成如下簡潔形式：

Rpp(θ)=A(θ)Rf+B(θ)Rμ+C(θ)Rρ

(3)

其中：

(4)

1.2 疊前道集正演方程

根據(jù)上述Russell的f-μ-ρ反射系數(shù)方程和褶積理論，地震記錄可以表示成地震子波與反射系數(shù)的褶積，即：

S=W*R

(5)

其中：W為地震子波矩陣;R為反射系數(shù)矩陣，具體由上述公式(3)中的Rpp(θ)描述；S為地震記錄(角度道集)。

1.3 目標(biāo)函數(shù)建立和求解

在流體因子疊前反演問題中，疊前角度道集數(shù)據(jù)y是已知的，即實測地震數(shù)據(jù)，未知待反演求解參數(shù)x，其中x=(Rf,Rμ,Rρ)。地震子波作為已知參數(shù)給定，速度比參數(shù)γdry、γsat由具體工區(qū)巖石物理分析結(jié)果確定，此處作為已知參數(shù)給定。筆者把實測地震數(shù)據(jù)與待反演流體因子參數(shù)都作為隨機(jī)變量，所要建立的目標(biāo)函數(shù)由已知實測數(shù)據(jù)條件下關(guān)于流體因子參數(shù)的條件概率的形式來表達(dá)，待反演參數(shù)在條件概率取得最大值的時候體現(xiàn)為最優(yōu)解。根據(jù)貝葉斯理論，定義目標(biāo)函數(shù)F(x)的具體形式為式(6)。

(6)

(7)

在式(7)兩邊取對數(shù)，得到進(jìn)一步簡化的目標(biāo)函數(shù)：

F(x)=(f(x)-y)T(f(x)-y)+

(8)

其中：式(8)右邊第一項是數(shù)據(jù)匹配項；第二項是低頻約束項；λ是調(diào)節(jié)權(quán)值，用于調(diào)整兩項在目標(biāo)函數(shù)中的比重。目標(biāo)函數(shù)F(x)對未知參數(shù)x=(Rf,Rμ,Rρ)求偏導(dǎo)，得到梯度F(x)=(?F(x)/?Rf,?F(x)/?Rμ,?F(x)/?Rρ)：

(9)

用上述定義的目標(biāo)函數(shù)與其梯度信息，通過共軛梯度法或阻尼最小二乘法求解上述優(yōu)化問題.然后對反演結(jié)果進(jìn)行積分，得到最終的反演結(jié)果，積分公式為式(10)。

(10)

其中：t0為初始時間；f(t0)為反演參數(shù)的初始值。

2 算法應(yīng)用

2.1 模型試驗

選用如表1所示的致密砂巖模型參數(shù)，模型一描述了砂巖在完全飽水狀態(tài)和完全飽氣狀態(tài)的參數(shù)，模型中將完全飽水狀態(tài)設(shè)置為上層，將完全飽氣狀態(tài)設(shè)置為下層，上層和下層的骨架及孔隙參數(shù)一致，僅是飽含流體不同,Kdry/μ為0.893。模型二中上層為泥巖，下層為含水砂巖，令上層泥巖和下層水砂的Kdry/μ均為0.998。利用f-μ-ρ反射系數(shù)方程計算彈性邊界處的反射系數(shù)，同時與Zoeppritz方程的精確解、Bortfeld-Aki-Richards、Shuey近似式進(jìn)行比較。

圖1中Russell近似計算結(jié)果與Zoeppritz方程精確解計算結(jié)果雖然存在差異，但是反射系數(shù)隨入射角變化趨勢以及與精確解的誤差滿足應(yīng)用于地震資料進(jìn)行反演的要求。

表1 砂巖巖石物理性質(zhì)

圖1 反射系數(shù)Fig.1 Reflection coefficient(a)模型一的反射系數(shù)計算結(jié)果；(b)模型二的反射系數(shù)計算結(jié)果

圖2 地震子波Fig.2 Seismic wavelet

2.2 井旁道反演

圖3 測井曲線和角度道集Fig.3 Logging curves and angle gathers(a)測井曲;(b)角度道集

圖4 井旁角度道集和殘差道集Fig.4 Angle gather and residual gather(a)輸入的角度道集；(b)反演結(jié)果正演的角度道集；(c)殘差道集

為檢驗本文方法的效果，選用典型含流體模型進(jìn)行正演得到地震道集數(shù)據(jù)，然后對該道集數(shù)據(jù)進(jìn)行反演。圖2是主頻為30 Hz的地震子波，采樣率為1 ms。圖3是測井曲線和合成角度道集，在圖3(a)中，紅色虛線圈中部分(626 ms ～632 ms)的縱波速度和密度表現(xiàn)為降低趨勢、而橫波速度并未與縱波速度和密度同步發(fā)生降低趨勢，這于流體的剪切模量為“0”的物理基礎(chǔ)是相符的，在各向同性介質(zhì)速度計算公式中，飽含流體巖石的體積模量降低，剪切模量保持不變，密度變低，在速度上的響應(yīng)為縱波速度降低，橫波速度不變。雖然受密度降低影響，但橫波速度的變化與縱波速度變化相比，可近似為“0”。在各向同性介質(zhì)假設(shè)前提下，紅色虛線圈中部分的縱波速度、橫波速度和密度與上下圍巖相比，能夠有效地描述巖石飽含流體后的狀態(tài)。圖3(b)是利用地震子波、流體因子f、剪切模量、密度曲線正演得到的地震角度道集，該地震角度道集采樣率為1 ms，該角度道集共包含45個地震道，入射角為從0°到44°，間隔為1°。

圖5 單道波形匹配Fig.5 Single channel waveform matching

給定圖3中角度道集數(shù)據(jù)，利用本文方法從該道集數(shù)據(jù)中反演流體因子f、剪切模量和密度。圖4中(a)為輸入的角度道集，(b)為反演結(jié)果正演的角度道集，(c)為殘差道集。圖5中曲線為從圖4(a)和(b)中抽取的第15道，波形匹配良好，共計546個樣點(diǎn)，平均誤差為0.000 3，標(biāo)準(zhǔn)差為0.002 046，誤差頂?shù)變蓚?cè)較大，由正演數(shù)據(jù)截取所致。

圖6中，藍(lán)色線為未經(jīng)平滑處理的實際測井曲線，黑色線是低頻模型，紅色線是反演結(jié)果，其中626 ms～632 ms處為圖3輸入測井曲線紅色虛線圈中部分所描述的優(yōu)質(zhì)含氣層；圖6(a)中反演得到的流體因子在626 ms-632 ms處表現(xiàn)為低值特征，與測井曲線計算的流體因子特征一致，兩者基本吻合；圖6(b)反演的剪切模量與輸入剪切模量曲線在數(shù)值范圍和形態(tài)上基本一致；圖6(c)中反演的密度在626 ms-632 ms處與測井?dāng)?shù)據(jù)具有較好地匹配。對比反演預(yù)測結(jié)果與輸入測井曲線可以看出，兩者吻合度較高，驗證了本方法的有效性。

2.3 實際資料反演

選用一條二維線進(jìn)行反演，該數(shù)據(jù)是驗證AVO屬性提取、疊前反演算法效果的通用測試數(shù)據(jù)，由HRS軟件的測試數(shù)據(jù)提供。圖7為超道集和疊加剖面，其中疊加剖面有131道，每道250個采樣點(diǎn)，采樣率為2 ms，氣層段，疊后地震剖面為一明顯的亮點(diǎn)，圖7中紅色曲線為測井?dāng)?shù)據(jù)計算的流體因子f。圖8為原始測井曲線，依次為縱波速度、橫波速度、密度、GR、速度比、流體因子f、電阻率、含水飽和度，最右側(cè)Tos在630 ms處標(biāo)記的Gas范圍內(nèi)有8 m左右氣層，氣層段測井曲線為電阻率相對較高，約為8 Ω·m，密度相對較低，約為2.2 g/cm3,聲波速度相對較低，約為2 300 m/s。圖9為GR曲線與流體因子f交會，從圖9中可以看出，流體因子f能夠區(qū)分出氣層與非氣層，為基于彈性參數(shù)預(yù)測氣層提供了巖石物理基礎(chǔ)。圖10是由常規(guī)疊前反演縱波阻抗、橫波阻抗和密度換算的流體因子f剖面，圖11是本文方法直接反演的流體因子f，top1層為優(yōu)質(zhì)含氣層，bot為含氣層，采用本方法得到的流體因子f，剖面含氣特征跟非含氣特征差異更加明顯，氣藏的邊界更加清晰，表明了該方法的可靠性和優(yōu)越性。

圖6 反演結(jié)果Fig.6 Inversion results(a)為流體因子反演結(jié)果；(b)為剪切模量反演結(jié)果；(c)為密度反演結(jié)果

圖7 超道集和疊加剖面Fig.7 Super gather and post stack section(a)超道集；(b)疊加剖面

圖8 井1Fig.8 Logging curve

圖9 GR與F流體因子交會Fig.9 Crossplot of GR and F

3 結(jié)論和討論

圖10 由公式計算的流體因子FFig.10 Gassmann fluid factor

圖11 直接反演的流體因子FFig.11 Gassmann fluid factor

筆者提出了一種流體因子直接反演方法，利用疊前地震道集資料反演儲層流體因子屬性，該方法以流體因子疊前反射原理為基礎(chǔ)，將低頻約束引入反演目標(biāo)泛函，構(gòu)建了流體因子疊前同步反演目標(biāo)函數(shù)，并給出了梯度信息的解析解，正演精度對比分析顯示了Russell近似方程滿足疊前反演的精度要求。典型流體模型實驗結(jié)果顯示了反演曲線與實測曲線匹配良好，表明流體因子同步反演方法能夠有效地從疊前地震資料中反演預(yù)測流體敏感屬性。本文方法綜合了巖石物理基礎(chǔ)理論、AVO基礎(chǔ)理論、數(shù)學(xué)最優(yōu)化分析等多種理論實現(xiàn)，其中干巖石速度比、飽和巖石速度比等參數(shù)需要由靶區(qū)巖石物理測試數(shù)據(jù)統(tǒng)計分析獲取。

致謝

感謝評審專家對本文提供的建設(shè)性建議和編輯部的大力支持！