基于線性貝葉斯和改進的典型相關分析法的地震巖石物理反演

曹亞梅, 周輝*, 于波, 王玲謙,2

1 油氣資源與工程全國重點實驗室, CNPC物探重點實驗室, 中國石油大學(北京), 北京 102249 2 中國石油大學(北京)理學院, 北京 102249 3 東北石油大學地球科學學院, 大慶 163318

0 引言

隨著復雜巖性油氣藏或非常規油氣藏開發的深入,傳統的彈性參數反演方法已經不能完全滿足儲層預測和油藏描述的需求,這推動了巖石物理反演的廣泛應用,以獲取地震數據中包含的豐富的彈性和物性參數信息(Gallop, 2006; González et al., 2008; Spikes et al., 2007; 撒利明等, 2015).聯合地震與巖石物理反演是定量描述油氣藏的主要技術之一,在儲層研究中占據重要地位.反演方法通常有兩類,分別為一步法和兩步法.一步法將巖石物理模型融入到地震模型中形成地震巖石物理模型,實現從地震響應直接推斷儲層物性參數(Mazzotti and Zamboni, 2003; Bosch et al., 2007; Chiappa and Mazzotti, 2009; Aleardi et al., 2017;劉軍等,2022).兩步法根據AVO近似模型從地震響應估計彈性參數(Russell, 1988;Avseth et al., 2010; 印興耀等,2020;趙容容等,2022),然后基于巖石物理模型將彈性參數轉化為儲層物性參數(Doyen, 2007; 劉興業等, 2018; Mavko et al., 2020).對一步法來說,反演結果直接受到地震數據的約束,而實際地震響應可能不支持井位處儲層性質的測量結果.此外,在復雜地質條件下,彈性性質和巖石物理性質之間的關系是不同地質情況相聯合作用的產物,根據地震數據直接獲得儲層物性參數的地下空間展布難以推廣到整個研究區域(Azevedo and Soares, 2017).而兩步法利用反演獲得的彈性性質與感興趣的巖石物理性質直接校準標定,保證導出的巖石物理關系受地震觀測數據的約束,從而降低了反演結果對儲層參數初始模型的依賴性.

由于觀測地震數據的采集誤差及噪聲等問題,地震反演存在很強的不適定性(Varela et al., 2006; 王寧等, 2022).貝葉斯理論可以融合地震、測井等多尺度信息,在井震同時約束下反演目標參數,量化反演結果的不確定性,在地球物理領域中應用廣泛(Buland and More, 2003; Hansen et al., 2006).Grana等(2017)、Grana(2018)利用貝葉斯理論,結合先驗混合高斯分布反演復雜巖相儲層的彈性性質.Yu等(2021)在線性貝葉斯反演的基礎上提出去相關方法,減少了彈性參數間的串擾,提高反演結果的精確度和反演過程的穩定性.袁成等(2015)、Luo等(2018)也針對線性貝葉斯反演進行了大量的研究和廣泛的實際資料應用.

巖石物理模型是巖石物理反演的核心,其建立方法主要有兩種:一是依賴于理論巖石物理方程(De Figueiredo et al., 2017, 2018; Liu et al., 2018; Lang and Grana, 2018; 周家雄等, 2019;劉仕友等,2022),二是利用統計巖石物理關系(Mukerji et al., 2001; Bachrach, 2006; Grana and Rossa, 2010; 印興耀等, 2014; 張冰等,2018).當地下介質復雜時,理論巖石物理方程難以準確描述研究區物性參數與彈性參數之間的關系,且大多理論巖石物理方程是非線性的,這大幅增加了巖石物理反演的計算難度.Grana(2016)提出一種針對理論巖石物理模型的線性反演方法,采用泰勒一階近似對巖石物理模型線性處理,將線性化的模型納入貝葉斯框架,實現彈性與儲層物性性質的高效轉化.但線性化巖石物理模型僅適用于非線性程度較低的巖石物理模型,對高度非線性的巖石物理模型強行線性化可能會失效(張佳佳等, 2020).統計巖石物理關系主要以統計分析的手段建立彈性參數與物性參數之間的關系,在誤差可接受范圍內獲得二者間的線性近似關系,大大提高了計算效率,同時統計巖石物理建模可應用于不同工區,具有更廣泛的適用性(Bachrach, 2006; Dubrule, 2003; Gui et al., 2015; 印興耀等, 2014).

目前應用的統計巖石物理關系大多基于多元回歸分析,可能會產生過擬合現象,從而降低預測儲層物性參數的準確性.為克服彈性參數與儲層物性參數之間的過擬合問題,本文引入典型相關分析(Canonical Correlation Analysis, CCA)來構建統計巖石物理模型.典型相關分析通過構建多組典型相關變量,利用典型相關變量反映原數據間的耦合特征,在地震儲層表征(Alvarez et al., 2015)、儲層建模預測(Satija et al., 2017)及地質建模(印興耀等, 2019)中均有應用.Grana等(2022)將典型相關分析應用于巖石物理隨機反演,并獲得了高精度的巖石物理反演結果.然而,直接利用原數據提取典型相關變量來構建巖石物理關系容易出現計算不穩定的現象,為此,本文采用Barnett和Preisendorfer(1987)改進的典型相關分析方法(BP-CCA),有機結合主成分分析(Principal Component Analysis,PCA)和CCA,提取原變量間的主成分,對主成分應用CCA,從而提取原數據間的隱性相關特征,建立穩定的統計巖石物理關系,實現高效、準確的巖石物理反演.

為實現適用于不同復雜地質條件下的高精度儲層參數預測,本文提出了一種結合線性貝葉斯理論與BP-CCA的地震巖石物理反演方法.在線性貝葉斯框架下,從疊前地震數據中反演縱、橫波速度及密度等彈性參數,通過BP-CCA算法計算測井資料中彈性參數與儲層參數之間的隱式線性關系,將構建的統計巖石物理關系應用于巖石物理反演,獲得地下物性參數的反演結果.本文分別對基于BP-CCA的巖石物理反演和基于多元線性回歸的巖石物理反演進行了合成數據實驗和實際數據應用測試,驗證了本方法可以獲得更準確的結果.

1 基本原理

1.1 地震正演模型

在地震褶積模型框架下,觀測地震數據d可以表示為地震子波W與反射系數rpp的褶積:

d=Wrpp+e,

(1)

其中,e為誤差項.反射系數的計算十分復雜,為簡化計算,常用Zoeppritz方程的近似式來求取反射系數,其中Aki-Richards公式的應用最為廣泛(Aki and Richards, 2002),具體表達式為

(2)

式中,vp、vs和ρ分別代表巖石的縱波速度、橫波速度和密度,θ為入射角度.因此,式(1)所示的地震正演表達式可以寫為

d=Gm+e=WADm+e,

(3)

其中,G為地震巖石物理正演算子,m=((lnvp)T,(lnvs)T,(lnρ)T)T為對數彈性參數向量,A為反射系數矩陣,D為一階差分算子.

1.2 疊前線性貝葉斯反演

在貝葉斯框架下,假設彈性參數m服從高斯分布:

m～N(μm,Σm),

(4)

其中,μm和Σm分別代表m的先驗均值向量和先驗協方差矩陣.假設地震數據的噪聲也服從高斯分布:

e～N(0,Σe).

(5)

根據線性高斯理論和式(3),觀測地震數據服從的高斯分布函數為

d～N(Gμm,GTΣmG+Σe),

(6)

其中,Gμm和GTΣmG+Σe分別對應地震記錄的均值向量和協方差矩陣.根據貝葉斯理論,在地震數據約束下,彈性參數的后驗概率分布也服從高斯分布,后驗均值向量μm|d和后驗協方差矩陣Σm|d可表示為

(7)

1.3 基于BP-CCA的巖石物理反演

為獲得孔隙度、泥質含量等儲層物性參數的空間展布,本文從已知的彈性參數和儲層物性參數(如測井數據)中獲得兩者之間的多元統計關系,從而建立起能夠最大程度逼近準確彈性參數與儲層參數關系的統計巖石物理模型.

1.3.1 統計巖石物理關系構建

傳統的統計巖石物理模型大多通過多元回歸構建,然而大部分巖石物理關系通常是非線性的,這對巖石物理關系的整合表征提出了挑戰.本文采用BP-CCA算法來定義彈性參數與物性參數之間的隱性關系.BP-CCA是一種結合了PCA和CCA的統計預測方法,原始變量通過PCA可轉化為用少數主分量表征主要特征的近似變量,降低信息冗余,而CCA可以將彈性與物性參數這兩組變量整合為多組具有最大相關性的典型變量,此時彈性與物性參數之間的統計關系就轉化為多組典型變量之間的相關關系.

m=μm+σmm′,

(8)

z=μz+σzz′,

(9)

其中,μm和μz分別為m和z的均值向量,σm和σz為標準差,m′和z′為標準化后的參數.假設m′和z′之間存在線性關系:

z′=Fm′+ε,

(10)

其中,F是一個表征m′與z′之間統計關系的多元回歸矩陣.

為獲得矩陣F的最優值,引入殘差函數g(F)并滿足以下條件:

g(F)=z′-Fm′.

(11)

根據最小二乘法,當殘差g(F)最小時,回歸矩陣F表達式為

F=(z′m′T)(m′m′T)-1.

(12)

通常情況下,實際數據中存在噪聲干擾,F的求取具有較強的不適定性,直接求解式(12)將引起求解不穩定,降低結果的可靠性.因此,針對式(12)進行正則化處理,結合CCA的思想,提高求解F過程的穩定性.具體做法如下:

首先,對標準化后的彈性參數m′與儲層參數z′應用PCA:

(13)

(14)

其中,Um和Uz為表征相同參數空間相關性的特征矩陣,Vm與Vz為表征不同參數之間相關性的特征矩陣,以上均為正交矩陣.將m′映射到前k個主分量上,將z′映射到前l個主分量上,即:

(15)

(16)

(17)

(18)

其中,R和P分別為左、右奇異矩陣,奇異值Scca的對角元素為典型相關系數,用于表征典型相關變量的相關性.此時,矩陣F可進一步寫為

(19)

(20)

(21)

此時,式(19)表示的是通過測井數據構建的統計巖石物理關系.

1.3.2 儲層物性參數預測

根據構建的巖石物理關系進行儲層參數預測的過程可分為以下5步:

(1)將線性貝葉斯反演獲得的彈性參數mo作為巖石物理反演的觀測數據,標準化獲得m′o:

(22)

(23)

(24)

(25)

(5)根據儲層參數的均值和協方差矩陣將預測的標準儲層參數恢復,獲得實際預測儲層物性參數zp:

zp=μzp+σzpz′p,

(26)

其中,μzp和σzp分別為從先驗信息中獲得的物性參數均值和標準差.

2 合成數據測試

為了驗證本文提出方法的有效性,采用Stanford模型的一部分進行測試.圖1a、c、e展示了縱、橫波速度及密度等彈性參數的三維真實模型,圖2a、c為孔隙度和泥質含量等物性參數的三維真實模型.該三維模型在時間方向上包含80個采樣點,采樣間隔為1 ms,主測線和聯絡測線方向上分別包含150和200道數據,測線間隔為20 m.以50 ms處的時間等時切片作為目的層,圖1b、d、f分別為準確模型中縱、橫波速度及密度的目的層等時切片,圖2b、d分別為孔隙度和泥質含量的目的層切片.基于Aki-Richards公式,利用主頻為30 Hz的Ricker子波與圖1a—c所示的彈性參數真實模型合成角道集正演記錄,入射角度分別為9°、15°和21°.添加隨機噪聲使信噪比為6 dB,形成如圖3所示的含噪角道集,用于后續的一維、二維及三維反演驗證測試.

圖1 三維真實彈性參數模型及目的層切片

圖2 三維真實物性參數模型及目的層切片

圖3 三維含噪角道集,入射角分別為(a)9°、(b)15°、(c)21°

為直觀分析BP-CCA方法構建統計巖石物理關系的過程,沿主測線80道與聯絡測線40道處抽取一道地震數據進行一維反演.用于獲取彈性參數與物性參數之間多元統計關系的偽井數據是通過隨機抽取三維真實模型中任意三道數據得到的,共計240個樣本點,彈性和物性參數對應的最大PCA主分量個數分別為k=3和l=2,在應用BP-CCA時,k和l越大,統計巖石物理關系越準確,但計算量也隨之增大,因此需要確定最優k和l值.根據偽井數據,以相關系數和均方根誤差來定量評價不同主分量組合(l,k)預測物性參數與真實物性參數的關系,如圖4所示.可以看出,k=2和l=1分別為彈性和物性參數的最優主分量,獲得圖5所示的典型相關變量對,其中灰色散點代表典型相關變量值,黑線為散點趨勢線,典型相關系數為0.9987,表征典型相關變量對的相關性,該統計巖石物理關系達到了用于巖石物理反演的標準.

圖4 偽井數據獲得的預測與真實物性參數間的相關系數與均方根誤差直方圖

圖5 BP-CCA獲得的第一典型相關變量對散點圖

圖6a、b分別為彈性參數和物性參數的一維反演結果,其中,黑色實線為反演結果,灰色點線為真實模型,灰色虛線為初始模型.對比發現,一維彈性和物性參數的反演結果均達到了較高的準確度,基本與真實模型吻合.利用均方根誤差對反演結果進行量化評價,圖6a的縱波速度、橫波速度及密度的反演結果與真實模型之間的均方根誤差分別為0.0179 km·s-1、0.0281 km·s-1及0.0478 g·cm-3,圖6b的孔隙度、泥質含量的反演結果與準確模型之間的均方根誤差分別為0.0127和0.0206,量化評價結果表明,噪聲干擾并未掩蓋彈性和物性參數的有效信息,本文提出方法可以有效表征儲層特性.

圖6 一維反演結果

接下來,沿主測線80道處切取二維地震數據進行二維反演測試.切得的二維地震數據如圖7所示,對應的彈性參數和物性參數真實模型分別如圖8a—c和圖8d—e所示.圖9a—c為基于線性貝葉斯理論的彈性參數反演結果,圖9d—e為基于BP-CCA的物性參數反演結果.二維測試是以單道反演的方式實現的,通過計算單道的均方根誤差來定量評價二維反演結果的精度,如圖10所示,其中圖10a—e分別為縱波速度、橫波速度、密度、孔隙度和泥質含量的反演結果與準確模型之間的均方根誤差.可以看出,隨機噪聲對反演結果造成一定程度的干擾,使得反演結果與真實模型之間存在些許差異,但從均方根誤差來看,噪聲干擾造成的反演誤差仍在可接受范圍內,且均方根誤差呈現出標準正態分布的狀態,與前文假設印證,從而證實了本方法的有效性和抗噪性.

圖7 二維含噪角道集,入射角分別為(a)9°、(b)15°、(c)21°

圖8 二維真實模型

圖9 二維反演結果

圖10 二維反演結果的均方根誤差

進一步以單道反演的方式測試本方法在三維合成數據上的表現,圖11a、c、e和圖12a、c所示為三維彈性參數和物性參數反演結果,對應地,圖11b、d、f和圖12b、d分別展示彈性和物性參數反演結果的目的層切片.對比圖1和圖2中真實模型的目的層切片可知,反演結果的切片形態與之基本吻合,且整體連續性較好,地質規律明顯.圖13為三維反演結果的單道均方根誤差分布直方圖,圖14為反演結果目的層切片與真實目的層切片間的相對誤差,誤差在可接受范圍內,達到目標精度.

圖11 三維彈性參數反演結果及目的層切片

圖12 三維物性參數反演結果及目的層切片

圖13 三維反演結果的均方根誤差

圖14 反演結果目的層切片相對誤差

3 實際數據測試

合成數據測試證明了本方法的有效性,現采用一個實際野外數據集測試本方法的應用效果.圖15為二維部分疊加角道集,其中角道集的入射角度分別為7.5°、17.5°和27.5°,圖中所示地震記錄位于1.4～1.8 s之間,時間采樣間隔為1 ms.二維剖面橫向上有340道數據,道間距為25 m,橫向跨度8.5 km.用1口測井數據參與本次測試,該井位于第144道,即3.6 km處.參與獲取多元統計關系的實原始樣本點為400個,為了避免過擬合,采用核密度估計(Silverman, 1986)對原始樣本點進行擴充.核密度估計是一種從一組數據點估計其非參數概率密度函數的統計方法,可以通過高斯核隨機抽取新樣本從而達到擴充原始樣本集的目的,使得新樣本點繼承原始樣本集的特征,擴充樣本點后孔隙度與泥質含量的核密度散點圖如圖16所示.將擴充后的樣本點應用BP-CCA算法,主分量個數分別為k=2和l=1.圖17為根據樣本集中計算的典型相關變量對散點圖,灰色散點代表典型相關變量值,黑色實線為散點趨勢線,典型相關系數為0.9927.圖18為基于BP-CCA算法預測的標準化物性參數與標準化測井數據的散點交會圖,其中黑色實線為散點趨勢線,孔隙度的預測值與真實值之間的相關系數為0.8453,泥質含量的預測值與真實值之間的相關系數為0.8012,說明獲得的統計巖石物理關系可以較好地表征彈性參數與物性參數之間的關系.

圖15 實際疊前數據集,入射角分別為

圖16 孔隙度與泥質含量的核密度散點圖

圖17 基于BP-CCA的第一典型相關變量對散點圖

圖18 基于BP-CCA預測物性參數與真實測井數據的散點圖

將構建的統計巖石物理關系應用于井旁地震道數據的反演,并與測井數據進行對比展示,以進一步檢驗本方法的準確性.彈性和物性參數的反演結果分別如圖19a、b所示,其中黑色實線為反演結果,灰色點線為測井數據,灰色虛線為測井數據濾波平滑得到的初始模型.可以看出,反演結果與測井數據之間的吻合程度較高,縱波速度、橫波速度、密度、孔隙度和泥質含量的反演結果與測井數據之間的均方根誤差分別為0.2103 km·s-1、0.1464 km·s-1、0.0897g·cm-3、0.0467和0.0590,表征結果精度較高.為進一步驗證BP-CCA方法相較傳統統計巖石物理反演的優越性,對井旁地震數據進行基于統計回歸的巖石物理反演(Grana and Rossa, 2010),結果如圖19c所示.反演結果與測井數據之間有較好的擬合,然而二者之間的均方根誤差分別為0.0652和0.0622,反演精度小于本文提出的方法.

圖19 井旁道反演結果

接下來實現工區二維反演,圖20a—c和圖20d—e分別展示了彈性及物性參數的反演結果,為便于對比,將已知井的彈性參數及物性參數顯示于二維反演剖面中.可以觀察到,井旁道反演結果與測井數據匹配良好,剖面整體反演結果連續性較高.

圖20 實際數據集反演結果

4 結論

本文基于改進的典型相關分析(BP-CCA)方法構建了穩健的統計巖石物理關系,并結合線性貝葉斯地震反演發展了一種兩步地震巖石物理反演方法,可實現從地震數據精確、高效地預測目標儲層的物性參數.本方法的新穎之處在于,從降維處理的原始彈性及物性參數中求取典型相關變量,利用典型相關變量獲得簡約且穩定的最大線性相關關系,進而建立彈性參數與儲層物性參數間的隱式統計巖石物理模型.本方法的優勢是,可實現從測井數據中直接獲得統計巖石物理模型,有效避免巖石物理模型的校準.此外,BP-CCA作為一種統計分析方法,可拓展到構建任意彈性與物性性質之間的統計巖石物理關系,具有廣泛的適用性.對比傳統的基于多元回歸的巖石物理反演方法,本方法反演精度更高,可實現對儲層參數的精準描述.