周林, 廖建平, 李景葉, 陳小宏, 劉興業
1 湖南科技大學頁巖氣資源利用湖南省重點實驗室,湘潭 411201 2 中國石油大學(北京)油氣資源與探測國家重點實驗室,北京 102249 3 中國石油大學(北京)海洋石油勘探國家工程實驗室,北京 102249 4 西安科技大學地質與環境學院,西安 710054
巖性和流體識別能夠為儲層含油氣性預測提供關鍵信息,其中流體識別的精度會直接影響儲層油氣分布預測精度.因此,如何合理高效的對儲層流體進行識別一直以來都是油藏勘探開發關注的重點,學者們也在這方面進行了大量研究.Smith和Gidlow(1987)通過聯合縱橫波速度首次定義了流體因子的概念并成功將其應用于含氣砂巖儲層的流體識別.Goodway等(1997)指出拉梅巖石物理常數λρ和μρ在流體檢測和巖性識別方面明顯優于縱橫波速度.基于孔彈性理論(Biot, 1941; Gassman, 1951),Russell等(2003, 2011)定義了具有明確物理意義的流體因子f,即經典的Russell流體因子.Feng等(2007)通過對常用流體因子比較分析指出ρf在流體判別方面表現更好.Yin和Zhang(2014)將有效孔隙流體體積模量Kf定義為新的流體指示因子,有效提升了流體識別的靈敏度.然而,其研究成果顯示,有效孔隙流體體積模量的引入會導致待反演的目標參數數量增加,這會嚴重影響非線性反演算法的穩定性和不確定性.李春寧(2014)構建了一個新的流體指示因子Fρ.相比于經典的Russell流體因子f和孔隙流體體積模量Kf,該流體因子不包含干巖石縱橫波速度比.并且,冷雪梅等(2019)通過研究指出,該流體因子的流體識別效果與Russell流體因子非常接近.從上述研究可以看出,流體因子的種類繁多且各有優缺點,本文將根據后續反演算法的需求,并借助楊培杰等(2016)提出的敏感流體因子定量分析法對現有流體指示因子進行優選.由于儲層的復雜性和多樣性,如果脫離巖性指示信息的約束,僅依靠流體因子反演結果對儲層含油氣性進行預測可能會導致解釋結果存在假象.如果能夠同時得到對儲層流體和巖性指示性明確的參數,共同對儲層含油氣性進行預測將會大大降低預測結果的不確定性.大量研究表明,泊松比具有良好的巖性和含油氣性指示效果(樊長江和王賢, 2006; Ostrander, 1984; 宗兆云等, 2012b; Zong et al., 2013a;Zhou et al., 2017; Zong and Yin, 2017).因此,為了有效提升儲層含油氣性預測精度,本文將以流體因子和泊松比同時作為目標待反演參數.
疊前地震數據中包含著豐富的地下流體和巖性信息,因此疊前AVO/AVA反演技術被廣泛應用于儲層流體檢測和巖性識別.Russell等(2011)推導了包含Russell流體因子、剪切模量和密度的線性近似公式,并采用標準最小二乘反演方法實現了線性反演.宗兆云等(2012a)、Zong等(2012)推導了包含縱波模量和橫波模量的線性近似方程以及彈性阻抗方程,進而基于巖石物理關系實現了Russell流體因子的間接估算.李雷豪等(2019)從減少參數維數入手將上述包含縱橫波模量的三項線性近似方程簡化為兩項式,提升了反演結果的穩定性,并最終通過間接計算實現了高敏流體識別因子的估計.為了避免間接計算造成的累積誤差,Zong等(2013b)利用Russell線性近似方程基于貝葉斯理論實現了Russell流體因子的直接反演.宗兆云等(2012b)、Zong等(2013a)推導了包含楊氏模量和泊松比的近似公式,并實現了泊松比的直接反演.Du和Yan(2013)推導了包含ρf和ρμ的PP和PS波反射系數兩項近似式,并基于這兩個方程實現了縱橫波聯合反演.杜炳毅等(2016)基于Russell近似推導了PS波反射系數近似公式并實現了聯合反演.Zong和Xin(2017)推導了包含楊氏阻抗和泊松阻抗的線性近似方程,并實現了基于該方程的疊前反演,有效提升了非常規儲層流體和巖性的識別精度.Yin和Zhang(2014)、Zong等(2015)、Du等(2019)推導了包含孔隙流體體積模量的線性近似公式并實現了孔隙流體模量的線性反演.冷雪梅等(2019)推導了包含流體因子Fρ的彈性阻抗方程,并實現了基于該方程的聯合反演.吳奎等(2021)從Russell近似表達式出發推導了包含Russell流體因子的彈性阻抗方程,并實現了基于該方程的流體因子反演預測.然而,上述反演方法均是以精確Zoeppritz方程的近似公式為正演算子.眾所周知,近似公式推導過程中的諸多假設條件會極大地限制上述反演方法的應用效果,尤其是在具有強阻抗差異特征的儲層.為了克服近似公式引入的一系列問題,本文將聚焦如何利用精確Zoeppritz方程實現流體因子和泊松比的直接反演.
基于貝葉斯理論構建反演目標函數,并在先驗模型中有針對性地引入合理的先驗分布特征能夠進一步有效提升反演結果的精度(Alemie and Sacchi, 2011).Theune等(2010)、Zhou等(2020)通過研究指出,在先驗模型中引入服從微分拉普拉斯分布的塊約束項能夠通過壓制旁瓣來更好地表征儲層邊界.因此,本文將采用同樣的策略來進一步提升流體因子和泊松比反演結果對儲層邊界的刻畫效果.
本文首先對現有常用流體指示因子利用感流體因子定量分析法進行了優選,然后將傳統形式的精確Zoeppritz方程改寫為包含優選流體因子、泊松比和密度的新形式,并基于新形式的精確Zoeppritz方程構建了貝葉斯框架下的非線性反演目標函數,有效地避免了常規基于近似公式的反演方法對反演結果精度的限制.此外,為了進一步提升反演結果對儲層邊界的刻畫效果,在先驗模型中引入服從微分拉普拉斯分布的塊約束項.合成數據和油田數據應用表明,新方法能夠高精度地估計流體因子和泊松比,且最終能夠基于二者反演結果實現儲層油氣分布特征的準確預測.

(1)


圖1 流體因子敏感系數柱狀圖(第Ⅰ~Ⅲ類AVO含流體砂巖模型)

圖2 流體因子敏感系數柱狀圖(第Ⅳ類AVO含流體砂巖模型)
(2)
傳統形式的精確Zoeppritz方程為(Zoeppritz, 1919; Aki and Richards, 1980):
(3)
式中,VP1、VP2、VS1、VS2、ρ1和ρ2分別表示反射界面兩側的縱橫波速度和密度,θ1和θ2分別表示P波的入射角和透射角,φ1和φ2分別表示PS轉換波的反射角和透射角,RPP和RPS分別表示P波和PS轉換波的反射系數,TPP和TPS分別表示P波和PS轉換波的透射系數.根據Snell定律可得:
(4)
將公式(4)代入公式(3)可得:
(5)
在各向同性介質中,縱波模量、楊氏模量與縱橫波速度以及密度存在如下關系:

(6)
(7)
其中,M1和M2為上下層縱波模量,E1和E2為上下層楊氏模量.將公式(6)和(7)代入公式(2)可得:
Fρ1=ρ1M1-ρ1E1,Fρ2=ρ2M2-ρ2E2.
(8)
在各向同性介質中,縱橫波模量與泊松比、楊氏模量存在如下巖石物理關系:
(9)
(10)
其中,σ1和σ2為上下層泊松比.將公式(9)代入公式(8)可得:
(11)

(12)
將公式(9)和(11)代入公式(10)可得:
(13)

(14)
將公式(12)和(14)代入公式(5),即可得到包含流體因子Fρ、泊松比σ和密度ρ的新形式精確Zoeppritz方程:

(15)
此外,為了便于其他學者對全文研究內容的重現以及后續反演算法的對比,本文在附錄A中給出了包含流體因子Fρ、泊松比σ和密度ρ的傳統線性近似公式的詳細推導過程.
當精確Zoeppritz方程的參數化形式發生改變時,其為反演帶來的不穩定性也會發生變化.因此,在反演之前我們借助殘差函數圖分析法對新方程的反演穩定性進行了分析(Larsen, 1999;Zhou et al., 2020),結果顯示新方程的殘差函數圖都是光滑的,沒有不規則的,且大部分都圍繞單極值呈現閉合的等值線,根據Zhou等(2020)的研究可知,本文推導的新形式精確Zoeppritz方程不會為反演算法帶來明顯的不穩定性,能夠成功地應用于貝葉斯確定性反演中.
以公式(15)所示新形式精確Zoeppritz方程為正演算子,構建貝葉斯理論框架下的非線性反演目標函數.假設似然函數P(d|m)服從高斯分布,則有:
(16)

大量研究表明,如果存在質量可靠且匹配良好的多波數據,多波聯合反演能夠進一步降低反演的不確定性,提升反演結果的精度(Stewart, 1990; Lu et al., 2015; 周林等, 2016; Zhou et al., 2017).
對于多波數據,公式(16)可推廣為:

(17)
基于貝葉斯理論構建目標函數的優勢是能夠通過引入更加有針對性的先驗分布模型進一步提升待反演參數的精度和分辨率.Theune等(2010)、Zhou等(2020)等學者指出,在先驗模型中引入服從微分拉普拉斯分布的塊約束項能夠通過壓制旁瓣來更好地表征儲層邊界.因此,為了進一步提升流體因子反演結果對儲層流體邊界以及泊松比反演結果對巖性邊界的刻畫能力,先驗分布函數P(m)被定義為(Zhou et al., 2020):
(18)
其中,Cm為包含三個待反演參數之間統計相關性的協方差矩陣,μ為均值向量,D為一階微分矩陣算子,kl(l=1,2,3)為尺度參數,每一個待反演參數對應一個尺度參數.
將公式(17)和(18)代入貝葉斯理論中,即可將求解最大后驗概率解的問題轉化為求解如下所示目標函數最小值的問題:
(19)

對于公式(19)所示非線性目標函數,通??梢圆捎酶咚?牛頓法或借助泰勒級數展開對其進行快速穩定求解(周林等, 2016; Zhou et al., 2017),為了避免求解正演算子關于待反演參數的二階偏導數,本文將借助泰勒級數展開的方法對上述目標函數進行求解.根據Zhou等(2017)的研究,對非線性正演方程G(m)進行泰勒級數展開并取其一階近似,然后代入公式(19)可得:
(20)
此時,公式(20)所示目標函數就變成了一個關于目標參數擾動量Δm的線性方程,將公式(20)關于該擾動量求導并令導數等于零可得:
Δm=H-1γ,
(21)
其中:
基于上述求解表達式,我們可以得到目標參數更新項Δm0,從而實現對初始模型m0的更新:
m1=m0+Δm0.
(22)
將m1代入到公式(21),即利用更新后的模型參數m1替代初始模型m0,則又可以得到一個新的目標參數更新項Δm1,這樣我們就可以實現模型參數的連續迭代更新,并最終達到最大迭代次數.因此,公式(19)所示目標函數的解可表示為:
Δmj=(H(mj))-1γ(mj),
(23)
其中:
最終,可得到如下所示待反演參數迭代更新公式:
mj+1=mj+λjΔmj,j=0,1,2,…,
(24)
其中,λj為第j次迭代的步長.

為了驗證基于新形式精確Zoeppritz方程流體因子、泊松比直接反演方法的可行性和穩定性,本文利用圖3所示單井模型進行算法測試.圖中黑線表示真實待反演參數,虛線表示初始模型參數.利用精確Zoeppritz方程基于圖3黑線所示真實參數值計算得到對應的反射系數向量,并與主頻為30 Hz的雷克子波進行褶積,最終得到如圖4a所示不含噪聲的合成PP波和PS波角道集數據.該角道集數據包含4°到40°之間的10個角度(以4°為間隔).為了測試算法的穩定性和抗噪性,在合成數據中加入均方根信噪比為0.5的隨機噪聲,如圖4b所示.圖5和圖6分別為本文提出的基于新形式精確Zoeppritz方程流體因子、泊松比直接反演方法和基于傳統近似公式反演方法的反演結果.其中,圖5a和圖6a為PP波單獨反演的結果,圖5b和圖6b為PP-PS聯合反演的結果.從圖5和圖6的對比可以看出,新方法的反演精度明顯高于基于傳統近似公式的常規方法.此外,從圖5a、b對比可以看出,基于新形式精確Zoeppritz方程的聯合反演能夠進一步提升反演結果的精度.而圖6a、b的對比顯示,由近似公式計算得到的PS波反射系數存在較大誤差,使得基于傳統近似公式的聯合反演不僅未能有效改善反演結果的精度甚至還在一定程度上降低了反演結果的精度.這說明新方法不僅反演精度優于常規方法,而且在PS波信息利用方面也表現的更好,充分驗證了新方法的可行性和有效性.仔細觀察圖5和圖6,我們還可以發現另外一個現象,即不管是傳統方法,還是新方法,密度反演結果精度均低于流體因子和泊松比反演結果精度.對此,我們借助基于單一變量原則的參數變化敏感性分析手段進行了簡要分析,發現包含流體因子、泊松比以及密度項的新方程對密度參數變化的敏感性明顯低于流體因子和泊松比,這也很好地解釋了這種現象的存在.

圖3 單井模型

圖4 合成疊前角道集

圖5 新方法的反演結果

圖6 基于傳統近似公式的常規方法的反演結果
為了測試新方法的穩定性和抗噪性,利用圖4b所示信噪比為0.5的合成數據對研究算法進行測試.圖7a、b分別為含噪數據PP波單獨反演和PP-PS波聯合反演的結果,可以看出,即使在信噪比為0.5的情況下,新方法仍然能夠穩定合理地估計出流體因子和泊松比,說明新方法具有良好的穩定性和抗噪性,同時也很好地驗證了前述關于新方程反演穩定性分析的結論.

圖7 新方法的反演結果(S/N=0.5)
實際資料來自中國東部某勘探工區,圖8為測試區域的疊加剖面.該區域實際采集的角道集數據角度范圍為3°到45°(間隔1°),并已經過一系列常規處理滿足疊前AVA反演的要求.圖8中黑線標示位置為井A所在位置,圖9為井A對應的真實測井曲線.圖10為圖9所示測井數據中流體因子Fρ與泊松比σ曲線交會圖,從該圖可以看出,流體因子Fρ和泊松比σ交會能夠清楚地區分儲層和非儲層區域.

圖8 測試區域的疊加剖面

圖9 井A的真實測井曲線

圖10 圖9所示流體因子和泊松比曲線交會圖
由于工區實際數據的限制,僅有PP波數據參與反演.圖11是利用新方法反演得到的流體因子(圖11a)、泊松比(圖11b)和密度(圖11c)剖面.圖12是利用基于傳統近似公式的常規方法反演得到的流體因子(圖12a)、泊松比(圖12b)和密度(圖12c)剖面.通過對比可以看出,新方法的反演結果在垂向分辨率和橫向連續性方面均明顯優于常規方法.為了進一步展示新方法的優勢,圖13給出了井A位置對應的井旁道反演結果,圖中紅線代表真實測井曲線,藍線代表新方法反演結果,黑色實線代表常規方法反演結果,黑色虛線代表初始模型.通過對比可以看出,新方法流體因子、泊松比反演結果的精度明顯高于常規方法,充分體現了本文方法的優勢.

圖12 基于傳統近似公式的常規方法反演結果剖面

圖13 井旁道反演結果對比
從圖11a所示新方法流體因子反演結果可以看出,區域1和2均顯示明顯的流體因子高值,同時,從圖11b所示泊松比反演結果可以看出,區域1和2均顯示明顯的低泊松比值,根據圖10所示的交會結果很容易判定區域1和區域2位有明顯的含油氣性顯示.此外,從圖11a還可以看出,區域3所示流體因子反演結果雖然呈現高值,但明顯比區域1和2所示流體因子值要低,此時如果僅僅依靠流體因子反演結果很難判定區域3是否為有效儲層.但從圖11b可以看出,區域3所示泊松比值呈現異常低值.大量研究表明,含氣砂巖的泊松比值呈現異常低值,那么我們很容易判定該區域可能是含氣儲層,這也很好解釋了為什么該區域的流體因子呈現高值但卻遠低于區域1和2.基于上述分析,我們可以推測區域1和2 為含油儲層,區域3為含氣儲層,該分析結果與工區實際鉆遇結果一致,充分驗證了本文研究方法的可行性和有效性,對豐富儲層含油氣性高精度預測理論與方法具有十分重要的意義.

圖11 新方法反演結果剖面
本文首先通過流體因子敏感性分析優選出合適的流體因子,然后將傳統形式的精確Zoeppritz方程改寫為包含流體因子、泊松比和密度的新形式,并最終實現了基于新公式的貝葉斯非線性反演,有效地克服了傳統近似公式固有缺陷對流體因子和泊松比反演結果精度的限制.同時,為了進一步提升反演結果對儲層邊界的刻畫能力,有針對性地引入了服從微分拉普拉斯分布的塊約束項.合成數據測試結果表明,新方法能夠穩定合理地反演得到流體因子和泊松比且反演精度與PS轉換波信息的利用均比常規方法表現得更好.二維實際數據測試結果表明,新方法反演結果的精度、垂向分辨率以及橫向連續性均明顯優常規方法,并且,同時利用流體因子和泊松比反演結果對研究區域的含油氣性進行預測,能夠有效降低預測結果的不確定性,充分驗證了本文提出的基于精確Zoeppritz方程的儲層含油氣性預測方法的可行性和有效性.
本文選擇的流體因子雖然很好地避免了干巖石縱橫波速度比對Russell流體因子和流體體積模量流體識別能力的影響,但該流體物理意義并不明確,且在準確給出干巖石縱橫波速度比情況下,其流體識別能力弱于Russell流體因子和流體體積模量.理論上,如果能夠減弱或消除干巖石縱橫波速度比對流體因子識別能力的影響,選擇Russell流體因子或流體體積模量作為流體指示因子應該是最佳選擇.因此,在本文研究內容的基礎上,將隨巖性或深度(走時)變化的干巖石縱橫波速度比同時作為未知待反演參數,開展基于精確Zoeppritz方程的Russell流體因子或流體體積模量直接反演研究,達到消除干巖石縱橫波速度比的影響并同時提升流體識別精度的目的,將是我們下一步計劃研究的內容.但將干巖石縱橫波速度比視為待反演參數或直接反演流體體積模量都將導致待反演參數數量增加,會嚴重影響非線性反演的穩定性同時降低反演結果的不確定性,這也將是我們后續研究中首要解決的難題.
致謝感謝所有合作作者在本文構思、研究及寫作過程中的辛勤付出.感謝兩位審稿專家對本文提出的寶貴建議.
附錄A 包含流體因子Fρ、泊松比σ和密度ρ的線性近似公式推導
優選流體因子Fρ的表達為:
(A1)
泊松比的表達式為:
(A2)
根據多變量微積分鏈式法則可得:
(A3)
(A4)
將公式(A3)兩邊同時除以Fρ,將公式(A4)兩邊同時除以σ,則有:

(A5)
(A6)
其中γ=(VP/VS)sat表示飽和巖石的縱橫速度比.
將公式(A5)和(A6)聯立求解可得:
(A7)
(A8)
將公式(A7)和(A8)代入Aki-Richards近似方程即可得到如下所示包含流體因子Fρ、泊松比σ和密度ρ的傳統線性近似公式:

(A9)
(A10)
附錄B 新形式精確Zoeppritz方程關于反射界面上下層流體因子Fρ、泊松比σ和密度ρ的一階偏導數求解表達式
矩陣形式的新形式精確Zoeppritz方程表達式為:
AR=b,
(B1)

(B2)
其中m=(Fρ1,σ1,ρ1,Fρ2,σ2,ρ2).對式(B2)進行整理可得:
(B3)


(B4)
(B5)
令M=η1/ξ1,N=ξ1/η1,則有:
(B6)

(B7)

(B8)
(B9)

(B10)