丁 燕，杜啟振，劉力輝，符力耘，冷雪梅，劉子煊

（1.中國石油大學(xué)（華東）深層油氣重點實驗室，山東青島 266580；2.青島海洋科學(xué)與技術(shù)國家實驗室海洋礦產(chǎn)資源評價與探測技術(shù)功能實驗室，山東青島 266580；3.北京諾克斯達石油科技有限公司，北京 100083；4.中國石化國際勘探開發(fā)有限公司，北京 100029）

0 引言

近幾年，隨著疊前地震反演技術(shù)和多波多分量技術(shù)的發(fā)展，儲層流體識別方法的有效性和準確度不斷提高。利用疊前振幅隨偏移距或角度變化信息來進行流體識別是一種重要的方法［1-4］。印興耀等［5-6］通過含有Gassmann 流體項的彈性阻抗方程直接反演來提高流體識別的準確性；Zong 等［7］建立了縱波模量和橫波模量表示的孔隙彈性理論和反射系數(shù)近似方程，并發(fā)展了基于縱波模量和橫波模量AVO 反演的流體識別方法；張世鑫［8］發(fā)展了基于固液解耦的地震流體識別方法，克服了常規(guī)流體識別方法受孔隙度影響產(chǎn)生的流體識別假象，然而，這些方法僅使用了反射縱波地震數(shù)據(jù)，在巖性和流體的共同影響下，利用縱波得到的反演結(jié)果往往存在精度低、反演不穩(wěn)定、邊界刻畫不清晰等問題。

與縱波速度相比，橫波速度受流體影響較小，因此在反演中同時利用縱波和轉(zhuǎn)換橫波地震數(shù)據(jù)，可以降低反演的多解性，提高儲層預(yù)測精度［9-11］。Larsen 等［12］提出了縱波和轉(zhuǎn)換波AVO 聯(lián)合加權(quán)疊加反演方法；陳天勝等［13］基于一種方向加速度優(yōu)化算法，利用縱波和轉(zhuǎn)換橫波反演得到了穩(wěn)定的速度比值；Du 等［14］以彈性孔隙介質(zhì)理論為基礎(chǔ)，將流體因子引入Russell 反射系數(shù)公式，實現(xiàn)了多波聯(lián)合AVO 反演，但常規(guī)流體因子由于受孔隙度和孔隙流體的耦合作用影響，流體識別存在的多解性問題仍未得到有效解決。

綜合前人的研究結(jié)果，首先討論Gassmann 流體因子和流體體積模量隨孔隙度的變化關(guān)系，然后針對復(fù)雜儲層中Gassmann 流體因子受孔隙流體和孔隙度耦合影響而存在的流體識別假象問題，推導(dǎo)基于流體體積模量、剪切模量和孔隙度的PP 波和PS 波Zoeppritz 線性近似方程，建立一種能夠穩(wěn)定獲取等效流體體積模量、剪切模量和孔隙度的縱橫波同步聯(lián)合反演方法，以期解決單縱波反演精度低、反演不穩(wěn)定、邊界刻畫不清晰的問題。

1 理論方法

1.1 孔隙流體體積模量的敏感度分析

Domenico［15］假設(shè)流體的體積模量是氣體和液體體積模量的體積加權(quán)函數(shù)

式中：Kf為流體體積模量，GPa；Sw為含水飽和度，%；Kw為液體體積模量，GPa；Kg為氣體體積模量，GPa。

Han 等［16］提出了Gassmann 流體項經(jīng)驗公式

式中：f為流體因子項為巖石骨架礦物與孔隙度的綜合增益函數(shù)；Φ為孔隙度，%；Kdry為干巖石體積模量，GPa；Km為固體礦物基質(zhì)體積模量，GPa；Kn為干巖石體積模量與基質(zhì)體積模量的比值。

Nur 等［17］給出了臨界孔隙度模型

式中：μdry為干巖石剪切模量，GPa；μm為固體基質(zhì)剪切模量（GPa）；Φc為臨界孔隙度，%。對增益函數(shù)G(Φ)進行簡化，得到

假設(shè)固體骨架礦物相同，孔隙中包含水和氣，則Gassmann 流體因子和流體體積模量隨含水飽和度和孔隙度的變化趨勢（圖1）便可由式（1）與式（2）計算得到。由圖1 可看出，Gassmann 流體因子受孔隙流體和孔隙度的綜合影響，即與孔隙度和含水飽和度呈非線性關(guān)系［圖1（a）］，而流體體積模量在一定孔隙度下隨含水飽和度的變化趨勢則為完全線性［圖1（b）］。因此，如果將流體體積模量與孔隙度分別單獨作為變量，即可消除常規(guī)流體因子受孔隙度和孔隙流體的耦合影響，解決流體識別多解性的問題。

圖1 流體因子隨孔隙度與含水飽和度的變化趨勢Fig.1 Trend of fluid factor with porosity and water saturation

1.2 Kf?μ ?Φ 縱橫波孔隙模量三參數(shù)AVO 近似方程和反演

1.2.1Kf?μ?Φ縱橫波孔隙模量三參數(shù)AVO 近似方程推導(dǎo)

Russell 等［18］建立了彈性孔隙介質(zhì)理論和AVO技術(shù)之間的關(guān)系，推導(dǎo)了關(guān)于流體因子、剪切模量和密度的三參數(shù)AVO 近似式

式中：θ為PP 波的入射角與透射角的平均值，（°）；φ為PS 波反射角與透射角的平均值，（°）；μ為剪切模量，GPa；ρ為密度，g/cm3；γsat為飽和巖石條件下縱橫波速度比均值；γdry為干巖石情況下的縱橫波速度比。

Gassmann 流體項及剪切模量表達式為

式中：vP，vS分別為飽和巖石縱波、橫波速度，m/s。由式（7）和式（8）的流體因子、剪切模量變化率、縱橫波速度變化率與密度變化率的關(guān)系，可以得到

根據(jù)Gardner 等［19］提出的經(jīng)驗關(guān)系式，得到縱波速度與密度變化率之間的關(guān)系

結(jié)合式（9）—（11），可以推導(dǎo)得到

將式（12）分別代入式（5）和式（6），并可將其簡化得到

根據(jù)Han 等［16］總結(jié)得到的流體體積模量Kf與Gassmann 流體項f之間的關(guān)系，聯(lián)立式（2）和式（4）可以得到

將式（15）代入式（13）和式（14），可以得到

式（16）和式（17）即為推導(dǎo)的PP 波和PS 波三參數(shù)Zoeppritz 線性近似公式。

1.2.2 近似方程精度分析

為驗證推導(dǎo)方程的精度［20］，分別設(shè)計了3 種砂巖模型（表1），分別用精確Zoeppritz 方程［21］，Russell方程［18］與式（16），式（17）計算得到地層界面處PP波、PS 波的反射系數(shù)，其中=2.333。當入射角不超過30°時，基于Kf?μ?Φ的縱波近似式與Russell近似式精度相近，而轉(zhuǎn)換波近似方程的精度則略高于Russell 近似式精度［圖2（a）—（b），圖2（e）—（f）］?；贙f?μ?Φ的縱波近似式與Russell 近似精度相近，而轉(zhuǎn)換波近似方程精度在大于30° 時變差［圖2（c）—（d）］。綜合以上分析，新推導(dǎo)的孔隙模量三參數(shù)AVO 近似方程在入射角0°～30°內(nèi)具有較高的精度，因此在最佳入射角0°～30°時利用本文提出的Kf?μ?Φ方程進行縱橫波同步聯(lián)合反演具有可行性。

表1 砂巖模型參數(shù)Table 1 Sandstone model parameters

圖2 3 種模型的PP 波、PS 波反射系數(shù)對比Fig.2 Comparison of PP wave and PS wave reflection coefficient of three kinds of models

1.3 基于貝葉斯理論的疊前同步聯(lián)合反演

對PP 波和PS 波進行同步聯(lián)合反演，考慮M層n個入射角的情況，將式（16）和式（17）寫成矩陣形式

根據(jù)地震褶積模型理論，增加子波矩陣至式（18），則疊前地震道集可以通過矩陣乘法表示為

式中：WP和WS分別為PP 波和PS 波的子波矩陣。

式（19）可以寫成

式中：d為PP 波和PS 波地震記錄矩陣；G為含有子波的系數(shù)矩陣；r為待反演的由等效體積模量、剪切模量和孔隙度組成的列向量系數(shù)矩陣，即

為提高反演的穩(wěn)定性，本文采用基于貝葉斯理論［22］的PP 波和PS 波同步聯(lián)合反演。貝葉斯反演的目的就是在給定測量數(shù)據(jù)d（帶噪聲n）的情況下估計模型參數(shù)r，把似然函數(shù)賦給噪音分布函數(shù)

將（24）式代入到邊緣化公式，取對數(shù)并進一步求導(dǎo)得到目標函數(shù)為

其中

令式（25）等于0，就得到問題的貝葉斯解

式中：μ為數(shù)據(jù)保真項與模型約束項之間權(quán)重的調(diào)節(jié)因子。

反演式（26）存在一定的弱非線性，為得到穩(wěn)定的求解結(jié)果，采用共軛梯度法進行求解。

2 應(yīng)用實例

2.1 Marmousi2 模型試算

為驗證本方法的正確性，本文利用Marmousi2油藏模型［23］進行測試。圖3 為利用縱、橫波速度，密度以及孔隙度等參數(shù)，計算得到模型等效流體體積模量、剪切模量和孔隙度（Kf?μ?Φ）剖面。在此模型參數(shù)基礎(chǔ)上，通過式（16）和式（17）正演得到反射系數(shù)，并與30 Hz 雷克子波褶積生成縱波和轉(zhuǎn)換波角度道集，同時加入信噪比為2 的高斯噪音。圖4是CDP 300 處PP 波與PS 波的角道集，其入射角為0°～30°。

圖3 等效流體體積模量（a）、剪切模量（b）和孔隙度（c）模型Fig.3 Models of equivalent fluid bulk modulus（a），shear modulus（b）and porosity（c）

圖4 CDP 300 處PP 波角道集（a）與PS 波角道集（b）Fig.4 PP wave angle gathers（a）and PS wave angle gathers（b）at CDP 300

利用中心角度為5°，15°，25°部分的疊加數(shù)據(jù)，分別進行單縱波和縱橫波同步聯(lián)合反演。由圖5中CDP 300 處單縱波、縱橫波同步聯(lián)合反演與實際模型單道對比結(jié)果可以得知，本文方法的反演結(jié)果誤差相對較小，穩(wěn)定性高于單縱波反演結(jié)果；圖6為單縱波和縱橫波同步聯(lián)合反演結(jié)果對比，其中，等效流體體積模量剖面能清楚地反映含油氣砂巖特征，剪切模量剖面反映了砂巖與頁巖的巖性差別，孔隙度剖面則準確反映了儲層的孔隙特征。綜合對比結(jié)果，采用本文縱橫波同步聯(lián)合反演方法得到的結(jié)果較單縱波反演結(jié)果具有更高的精確度，且能準確刻畫地質(zhì)體儲層特征和含氣特征。

圖5 CDP 300 處單縱波、縱橫波同步聯(lián)合反演結(jié)果與實際模型單道對比Fig.5 Single channel comparison of single P-wave and PP wave and SS wave with actual model at CDP 300

圖6 單縱波和縱橫波聯(lián)合反演對比Fig.6 Comparison of single P-wave and PP-wave and SS wave joint inversion

2.2 實際數(shù)據(jù)應(yīng)用

利用中國LJ 地區(qū)A 油田M 測線的縱波和轉(zhuǎn)換橫波數(shù)據(jù)進行應(yīng)用效果測試。該區(qū)位于凹陷南緩斜坡，受燕山運動和喜馬拉雅運動的影響，形成了向北傾沒的大型鼻狀構(gòu)造，受構(gòu)造控制，該區(qū)油氣藏類型主要為斷塊油氣藏，兼有巖性、泥巖裂縫等特殊油氣藏［24］。本區(qū)縱波地震資料主頻為35 Hz，轉(zhuǎn)換波地震資料主頻為15 Hz。由于實際縱波和轉(zhuǎn)換橫波記錄的反射振幅相差較大，為得到更可靠的聯(lián)合反演結(jié)果，需要對縱波振幅和轉(zhuǎn)換橫波振幅做匹配。以縱波振幅為標準，每一角度計算一個校正系數(shù)［式（27）］，為了不改變整個角道集的AVO 特征，將各個角度的校正系數(shù)求平均作為最終的校正系數(shù)，圖7 為在時間域內(nèi)匹配一致的縱波和轉(zhuǎn)換波角道集。

圖7 CDP 1000 點PP 波角道集（a）與PS 波角道集（b）Fig.7 PP wave angle gathers（a）and PS wave angle gathers（b）at CDP 1000

利用反演得到的M 測線目標層段局部放大的等效流體體積模量［圖8（a）］、剪切模量［圖8（b）］和孔隙度剖面［圖8（c）］。圖中黑色線柱表示生產(chǎn)井段內(nèi)的油層，紫色線柱表示生產(chǎn)井段內(nèi)的水層，經(jīng)過對比發(fā)現(xiàn)，油層對應(yīng)等效流體體積模量剖面低值［圖8（a）中紅色所示］，孔隙度剖面高值［圖8（c）中紅色所示］，而水層等效流體體積模量值高于油層［圖8（a）中藍色所示］，相對應(yīng)的孔隙度剖面低值［圖8（c）中藍色所示］。油水層所在砂巖與泥頁巖區(qū)在剪切模量剖面上有較明顯的區(qū)分［圖8（b）］。綜上所述，反演結(jié)果中儲層特征清晰，油、水層得到了較好的區(qū)分，且油層范圍與已知生產(chǎn)井的信息較吻合。

3 結(jié)論

（1）新推導(dǎo)的孔隙模量三參數(shù)AVO 近似方程在入射角0°～30°內(nèi)具有較高的精度，因此在最佳入射角0°～30°時利用本文提出的方程進行縱橫波同步聯(lián)合反演具有可行性。

（2）綜合理論分析與數(shù)據(jù)測試結(jié)果，本文提出的縱橫波同步聯(lián)合孔隙模量三參數(shù)AVO 反演方法改善了單縱波反演結(jié)果不穩(wěn)定和邊界模糊性的問題，并在一定程度上壓制了噪聲，提高了反演精度和流體識別的可靠性。