袁成, 李景葉*, 陳小宏

1 中國石油大學(北京)油氣資源與探測國家重點實驗室， 北京 102249 2 中國石油大學(北京)海洋石油勘探國家工程實驗室， 北京 102249

?

基于概率統(tǒng)計的地震巖相識別不確定性定量評價方法

袁成1,2, 李景葉1,2*, 陳小宏1,2

1 中國石油大學(北京)油氣資源與探測國家重點實驗室， 北京 102249 2 中國石油大學(北京)海洋石油勘探國家工程實驗室， 北京 102249

地震巖相識別能夠提供具有不同儲層特征的巖相分布信息，對巖相識別的不確定性開展定量評價分析可降低后期油藏建模與儲層評價的風險.考慮了地震巖相識別中測井巖相定義、巖石物理建模、井震尺度匹配及地震反演等環(huán)節(jié)的不確定性對巖相識別的影響，基于概率統(tǒng)計方法，引入熵函數(shù)實現(xiàn)了地震巖相識別不確定性定量評價，并結合巖相概率、重建率等多角度綜合定量分析不確定性的構成及傳遞特征，系統(tǒng)地實現(xiàn)了地震巖相識別不確定性評價流程的整體連通.提出了結合屬性交繪特征約束反演參數(shù)空間，提高地震巖相識別運算效率.模擬數(shù)據(jù)分析表明利用熵函數(shù)可精確實現(xiàn)巖相識別不確定性地定量表征，利用屬性交繪特征約束參數(shù)空間既大幅度減少運算量，也可降低地震巖相識別的不確定性.

概率統(tǒng)計； 熵函數(shù)； 巖相識別； 不確定性； 定量評價

1 引言

油田勘探初期測井數(shù)據(jù)稀少，儲層巖相識別主要依靠分布范圍較廣的地震資料.而地震巖相識別面臨較多的不確定因素，對其不確定性信息開展定量評價分析可降低后期油藏建模與儲層評價的風險.

對于地震巖相識別各環(huán)節(jié)不確定性的評價分析，國內(nèi)外學者做了大量地研究.Mukerji等(2001)首先在疊前地震反演中引入統(tǒng)計巖石物理預測儲層物性參數(shù)并表征了巖石物理關系的不確定性.Buland和Omre(2003)結合貝葉斯理論(Tarantola, 2005)與線性AVO反演技術反演彈性參數(shù)信息并通過協(xié)方差矩陣定量表征反演的不確定性.Larsen等(2006)論述了地震反演不確定性對于巖相識別的影響并結合馬爾科夫鏈提高反演精度.Bachrach(2006)聯(lián)合巖石物理、隨機建模方法和貝葉斯理論預測儲層飽和度與孔隙度，采用標準差(STD)分析預測結果的不確定性.Buland等(2008)引入了一種快速貝葉斯反演方法，并論證了準確表征反演不確定性的重要性.Grana和Rossa(2010)闡述了基于地震數(shù)據(jù)預測儲層物性屬性及不確定性分析方法.Grana等(2012a)通過建立儲層物性參數(shù)與巖相的統(tǒng)計關系獲取巖相概率，并定量表征巖相定義的不確定性.針對地震反演的不確定性，國內(nèi)許多學者也在隨機反演(張廣智等，2011；張繁昌等，2014)和多屬性聯(lián)合反演(曹丹平等，2009；楊鍇等，2012)方法研究上做了大量工作.田玉昆等(2013)以馬爾科夫隨機場為先驗模型，按照測井數(shù)據(jù)解釋結果由彈性參數(shù)直接劃分巖相，方法考慮了巖相與彈性參數(shù)間的不確定性，但沒有考慮地震反演的不確定性對巖相識別的影響.騰龍和程玖兵(2014)聯(lián)合統(tǒng)計巖石物理與AVA反演技術預測儲層物性屬性并分析了由于模型和數(shù)據(jù)誤差引入的不確定性.總之，目前關于地震巖相識別的不確定性評價主要開展于地震巖相識別的獨立環(huán)節(jié)，沒有實現(xiàn)整體流程的連通與全面系統(tǒng)地評價.為了實現(xiàn)地震巖相識別不確定性的定量評價，有效降低油藏勘探開發(fā)風險，急需系統(tǒng)地研究其不確定性信息的構成及傳遞特征.

研究在前人工作基礎上，引入熵函數(shù)實現(xiàn)了地震巖相識別不確定性構成及傳遞特征地定量表征，結合巖相概率、重建率等多角度綜合定量分析了其各個環(huán)節(jié)的不確定性對巖相識別的影響，以及地震噪聲對地震巖相識別的影響；并提出采用屬性交繪特征約束反演參數(shù)空間，提高運算效率.最后，基于模型資料詳細分析了地震巖相識別不確定性信息的構成特征，驗證了研究方法的有效性.

2 理論與方法

地震巖相識別包括測井巖相定義、巖石物理建模、井震尺度匹配以及地震數(shù)據(jù)反演等多個環(huán)節(jié)，各環(huán)節(jié)不確定性信息的傳遞和累積嚴重影響巖相識別的精度.因此，實現(xiàn)地震巖相識別不確定性構成和傳播特征的定量評價分析是提高地震巖相識別精度的必要前提.

2.1 測井巖相定義不確定性評價方法

測井巖相定義是地震巖相識別的重要環(huán)節(jié)之一.與傳統(tǒng)沉積相模式不同，在地震巖相識別中巖相是儲層參數(shù)的統(tǒng)計分類，將應用于后續(xù)的儲層建模(Grana et al., 2013).巖相內(nèi)測井物性參數(shù)的統(tǒng)計特征一般滿足高斯分布，公式為

P(R|F)=N(R;μR,ΣR)，

(1)

其中，F(xiàn)表示巖相，R是物性參數(shù)矢量(孔隙度、泥質含量等)，P(R|F)為給定巖相的物性參數(shù)概率分布.N為高斯函數(shù)，其均值和協(xié)方差矩陣分別為μR和ΣR.

根據(jù)貝葉斯理論，給定物性參數(shù)R條件下巖相概率P(F|R)可表示為

P(F|R)∝P(F)P(R|F)，

(2)

其中P(F)是巖相比例先驗信息，ΣP(F)=1，P(F|R)為測井巖相定義的條件概率.給定物性參數(shù)R，可通過P(F|R)獲取巖相F的概率信息.似然函數(shù)P(R|F)與巖相比例先驗信息P(F)的具體參數(shù)通過測井曲線聚類分析獲得(Pirroneetal., 2011).條件概率P(F|R)表征巖相F與物性參數(shù)R的統(tǒng)計關系，也蘊含著測井巖相定義的不確定性.為此，研究引入熵函數(shù)對其不確定性進行定量表征與評價.

熵描述體系的混亂程度，在眾多領域都有重要應用.統(tǒng)計學中熵函數(shù)反映系統(tǒng)的有序程度：狀態(tài)越有序，熵值越小.利用熵函數(shù)可定量表征地震巖相識別的不確定性信息：不確定性越大，熵值越大，反之亦然.傳統(tǒng)的基于反演概率信息或隨機反演獲取多個實現(xiàn)的不確定性評價方法是對不確定性的定性描述，無法對其開展定量評價分析；其次，傳統(tǒng)方法僅能夠對地震巖相識別最終的不確定性開展評價分析，而無法剖析其構成特征以及在地震巖相識別中的傳播規(guī)律.方法引入熵函數(shù)實現(xiàn)了不確定性的定量評價，通過分析對比各環(huán)節(jié)熵值信息，能夠對地震巖相識別中不確定性的構成特征和傳播規(guī)律開展定量表征.與傳統(tǒng)方法相比，引入熵函數(shù)為后期油藏建模與儲層評價風險評估提供了精確的不確定性定量信息.

概率矢量P(F|R)的信息熵(Shannon, 1948)可定義為

(3)

由式(3)可知，P(F|R)與概率均值1/Nf的距離越大，測井巖相定義的不確定性越小，熵函數(shù)值Epy(F|R)越小，反之亦然.因此Epy(F|R)可定量表征測井巖相定義的不確定性，是地震巖相識別不確定性的重要組成部分.

2.2 巖石物理不確定性評價方法

測井巖相定義確定了巖相和儲層物性參數(shù)的統(tǒng)計關系.由于地震數(shù)據(jù)是地層彈性屬性的地球物理響應，地震巖相識別需要建立彈性參數(shù)與物性參數(shù)的巖石物理關系從而進一步獲取巖相信息.在此引入統(tǒng)計巖石物理表征儲層物性與彈性參數(shù)的響應關系.

統(tǒng)計巖石物理是綜合多學科的地球物理方法，結合傳統(tǒng)巖石物理、信息論和統(tǒng)計學以表征傳統(tǒng)巖石物理模型的不確定性(Avseth et al., 2005).統(tǒng)計巖石物理可在給定彈性參數(shù)條件下預測物性參數(shù)概率分布，其中巖石的非均質性及測井數(shù)據(jù)測量誤差等可通過概率信息的不確定性量化表征(Grana et al., 2013).統(tǒng)計巖石物理的概念形式為

Ew=fRPM(R)+ε，

(4)

其中，Ew為測井彈性參數(shù)矢量，fRPM代表傳統(tǒng)巖石物理模型，ε表示誤差.

(5)

(6)

其中P(R|Ew)表示給定測井彈性參數(shù)Ew條件下物性參數(shù)R的概率分布，即統(tǒng)計巖石物理的概率形式.給定測井彈性參數(shù)Ew，通過P(R|Ew)可獲取物性參數(shù)R的概率分布.通過Chapman-Kolmogorov方程(Papoulis,1984)結合測井巖相定義的條件概率P(F|R)計算條件概率P(F|Ew)：

(7)

其中，P(F|Ew)為給定測井彈性參數(shù)Ew條件下巖相F的概率.形如式(3)，計算其熵函數(shù)值Epy(F|Ew).由式(7)可知，獲取巖相概率P(F|Ew)包含測井巖相定義及巖石物理建模兩個環(huán)節(jié)，故其熵函數(shù)值Epy(F|Ew)由測井巖相定義的不確定性Epy(F|R)和巖石物理響應關系P(R|Ew)的不確定性兩個部分構成.由于R為連續(xù)變量，P(R|Ew)不確定性信息無法采用式(3)方法計算，此處可由Epy(F|Ew)與Epy(F|R)之差表征.

2.3 井震匹配不確定性評價方法

地震巖相識別需要聯(lián)合井震信息.由于兩種數(shù)據(jù)測量尺度不同，信息融合過程中會不可避免地產(chǎn)生誤差，從而導致地震巖相識別不確定性的增加.針對其不確定性的評價，需要統(tǒng)計兩類數(shù)據(jù)的概率關系.

地震尺度彈性參數(shù)記為Es.形如式(5)和式(6)，統(tǒng)計其與測井尺度彈性參數(shù)Ew的條件概率P(Ew|Es)，然后結合P(F|Ew)計算P(F|Es)，公式為

(8)

其中P(F|Es)為地震尺度彈性參數(shù)Es條件下的巖相概率，其熵值Epy(F|Es)計算方法參見式(3).由式(8)可知，Epy(F|Es)在Epy(F|Ew)的基礎上進一步考慮了井震尺度匹配P(Ew|Es)的不確定性對巖相識別的影響.

2.4 地震反演不確定性評價方法

地震彈性參數(shù)反演是基于地震數(shù)據(jù)預測儲層彈性屬性的過程，其反演結果往往具有較大的不確定性.在地震彈性參數(shù)反演尋求最優(yōu)解的同時，需要對反演結果的不確定性進行定量評價(Buland and Omre, 2003).

地震彈性參數(shù)反演中，彈性屬性Es的先驗分布P(Es)一般滿足高斯模型：

P(Es)=N(Es;μEs,ΣEs)，

(9)

其中μEs和ΣEs分別為高斯模型的均值和協(xié)方差矩陣.

由線性反演理論可知，地震觀測數(shù)據(jù)S可近似表示為

S=GEs+e，

(10)

其中G是線性算子，e是地震噪聲.

根據(jù)統(tǒng)計學原理，變量概率分布形態(tài)在線性變換前后保持不變.故P(S)亦滿足高斯分布：

P(S)=N(S;μS,ΣS).

(11)

其均值μS與協(xié)方差矩陣ΣS分別滿足：

μS=GμEs，

(12)

ΣS=GΣEsG′+Σe，

(13)

其中Σe為地震噪聲的協(xié)方差矩陣.

設地震數(shù)據(jù)S與彈性參數(shù)Es的聯(lián)合分布滿足高斯模型.給定地震數(shù)據(jù)條件下彈性參數(shù)概率分布P(Es|S)可表示為

P(Es|S)=N(Es;μEs|S,ΣEs|S)，

(14)

其均值μEs|S和協(xié)方差矩陣ΣEs|S分別滿足：

μEs|S=μEs+(GΣEs)′(ΣS)-1(S-μS)，

(15)

ΣEs|S=ΣEs-(GΣEs)′(ΣS)-1(GΣEs)，

(16)

其中μEs|S是反演的最大后驗概率解(MAP)，ΣEs|S定量表征μEs|S的不確定性信息.

聯(lián)合地震反演結果P(Es|S)與條件概率P(F|Es)可得：

(17)

其中，P(F|S)為地震巖相識別的巖相概率，對應熵值Epy(F|S)計算參見式(3).Epy(F|S)在Epy(F|Es)基礎之上考慮了地震反演P(Es|S)的不確定性對巖相識別的影響，其定量表征地震巖相識別總體的不確定性.地震巖相識別各環(huán)節(jié)不確定性信息可結合Epy(F|R)、Epy(F|Ew)、Epy(F|Es)與Epy(F|S)等熵函數(shù)值加以量化分析，從而實現(xiàn)地震巖相識別不確定性構成與傳播特征的定量評價.

3 基于模型數(shù)據(jù)不確定性定量評價

為了驗證方法有效性，更好地認識地震巖相識別不確定性的組成特征，研究以一維砂泥巖儲層為例開展地震巖相識別不確定性定量評價；模型設有三種巖相，包括泥巖相、泥質砂巖相和砂巖相.在實際生產(chǎn)中，方法對儲層巖相間具有地球物理響應差異的儲層類型均具有適用性.模型巖相序列由馬爾科夫矩陣生成，巖相內(nèi)物性參數(shù)(孔隙度與泥質含量)聯(lián)合分布滿足高斯模型且均為水飽和，孔隙度和泥質含量相關系數(shù)為-0.8.各巖相物性參數(shù)統(tǒng)計如表1所示.

表1 巖相物性屬性統(tǒng)計參數(shù)Table 1 Statistical parameters of reservoir properties for lithofacies

假設各巖相物性參數(shù)曲線垂向變化規(guī)律相似可采取相同變差函數(shù)模型擬合.通過變差函數(shù)構建空間協(xié)方差矩陣，與各巖相的物性參數(shù)協(xié)方差矩陣進行克羅內(nèi)克乘積模擬生成孔隙度和泥質含量測井曲線.該方法既保證了物性參數(shù)的相關性，也保持了其垂向連續(xù)性(Grana et al., 2012b).結合硬砂巖巖石物理模型并加入適量噪聲正演獲取彈性參數(shù)曲線，并將其尺度放大到地震數(shù)據(jù)尺度.采用Zoeppritz方程計算反射系數(shù)，選取不同主頻的雷克子波(30 Hz、25 Hz和20 Hz)與對應入射角度(12°、24°和 36°)合成角度道集數(shù)據(jù)并加入適量噪聲(圖1).各巖相物性屬性條件概率P(R|F)與巖相比例P(F)的具體參數(shù)參見表1.通過式(2)計算測井巖相定義的條件概率P(F|R).

設任意巖相的物性參數(shù)分布滿足高斯模型，因此所有巖相總體滿足高斯混合模型(圖2b).高斯混合模型中每一個概率凸峰對應一個巖相.圖2a中各巖相數(shù)據(jù)點基本都落在相應巖相概率等值線內(nèi)，沒有明顯的混疊現(xiàn)象，但這并不代表測井巖相定義的概率關系P(F|R)沒有不確定性.因為對于任意的物性參數(shù)矢量R，巖相概率P(F|R)值總是位于(0,1)內(nèi).結合式(3)可知，其熵值Epy(F|R)是大于零的正數(shù)，代表P(F|R)的不確定性.因此在地震巖相識別中一定要考慮測井巖相定義自身所具有的不確定性.

由熵值Epy(F|R)統(tǒng)計計算過程可知，優(yōu)選巖相敏感屬性是降低測井巖相定義環(huán)節(jié)不確定性的有效方法.實際生產(chǎn)中，通過巖石物理分析確定對巖相較為敏感的測井屬性開展測井巖相定義.選擇巖相敏感屬性越多，測井巖相定義的不確定性相對越低.但應注意的是，屬性過多會嚴重影響算法運算效率，且抑制不確定性的效果趨于平緩.

3.1 巖石物理關系不確定性影響評價

圖1 儲層模型數(shù)據(jù)(a) 巖相序列(淺灰色為泥巖，深灰色為泥質砂巖，黑色為砂巖)； (b) 孔隙度； (c) 泥質含量； (d) 縱波速度； (e) 橫波速度； (f) 密度(實線為測井尺度，虛線為地震尺度)； (g) 合成地震記錄.Fig.1 Reservoir model data(a) Facies sequence (shale in light grey, silty sand in dark grey, sand in black); (b) Porosity; (c) Clay content; (d) P-wave velocity; (e) S-wave velocity; (f) Density (well-log scale in full line, seismic scale in dashed line); (g) Synthetic seismic data.

圖2 孔隙度與泥質含量交繪和概率等值線(虛線為參數(shù)空間邊界)(a)與聯(lián)合概率分布(b)Fig.2 (a) Cross-plot of porosity with clay content and contours of probability. Dashed line is the boundary of parameter space. (b) Joint probability distribution

針對巖石物理建模環(huán)節(jié)中儲層物性與彈性參數(shù)巖石物理響應關系的不確定性，引入統(tǒng)計巖石物理方法對其定量表征.通過統(tǒng)計巖石物理可建立彈性參數(shù)Ew與物性參數(shù)R的概率關系P(R|Ew).首先采用克羅內(nèi)克乘積方法模擬生成多條物性參數(shù)曲線R；再結合傳統(tǒng)巖石物理模型fRPM正演獲取多條彈性曲線Ew并加入適量噪聲ε；然后統(tǒng)計其聯(lián)合分布P(R,Ew)(式(5))進而計算條件概率P(R|Ew)(式(6)).P(R|Ew)為給定測井彈性參數(shù)Ew條件下物性參數(shù)R的概率分布，即統(tǒng)計巖石物理的概率形式.形如式(7)，結合條件概率P(F|R)獲取測井彈性參數(shù)Ew條件下的巖相概率P(F|Ew)，并計算其熵函數(shù)值Epy(F|Ew) (圖3).

圖3中，孔隙度與泥質含量的概率分布信息由測井彈性參數(shù)Ew通過條件概率P(R|Ew)反演獲得，滿足高斯分布.概率均值趨勢與模型參數(shù)基本吻合，驗證了巖石物理模型和統(tǒng)計關系的適用性；高斯分布形態(tài)表征了由數(shù)據(jù)噪聲引起的巖石物理響應關系的不確定性.巖相概率曲線0.33處的虛線表示三個巖相概率的均值，概率曲線與概率均值線的距離越大，巖相識別結果的不確定性越小.結合式(7)可知，熵值Epy(F|Ew)包含測井巖相定義的不確定性Epy(F|R)與巖石物理響應關系P(R|Ew)的不確定性.P(R|Ew)的不確定性信息可通過Epy(F|Ew)與Epy(F|R)之差表征.文中模擬數(shù)據(jù)的熵值Epy(F|R)較小，未在圖中顯示，但在實際評價不確定性的過程中一定要予以考慮.

實際生產(chǎn)中降低統(tǒng)計巖石物理環(huán)節(jié)的不確定性，首先應開展巖石物理診斷確定適用于工區(qū)資料的巖石物理模型，其決定統(tǒng)計概率均值.其次，精確統(tǒng)計分析已選模型對工區(qū)資料的巖石物理響應偏差，響應偏差決定統(tǒng)計方差特征，偏差統(tǒng)計過大會增加反演不確定性，過小則造成不確定性偏低的假象.

3.2 井震匹配不確定性影響評價

地震巖相識別需要聯(lián)合測井和地震信息.由于二者測量尺度的差異，井震信息聯(lián)合之前首先要開展井震尺度匹配.而該環(huán)節(jié)受井震資料品質影響較大，且尺度放大時窗寬度在深度域(λ/4)是隨速度變化的，在時間域時窗寬度(T/4)固定不變(Grana et al., 2013).為了準確評價井震尺度差異所引入的不確定性，而不受尺度放大算法和數(shù)據(jù)誤差等因素的影響，研究在時間域取固定時窗對測井彈性曲線Ew尺度放大到地震尺度彈性參數(shù)Es(圖1).針對巖相序列的尺度放大，引入柱狀圖粗化方法(histogramupscaling)，即首先統(tǒng)計時窗內(nèi)各巖相的樣本數(shù)，然后取統(tǒng)計柱狀高度最大的巖相賦給當前時窗對應位置.

為了準確評價井震尺度匹配環(huán)節(jié)引入的不確定性，將統(tǒng)計巖石物理環(huán)節(jié)中正演獲取的多條彈性參數(shù)Ew曲線尺度放大到地震尺度彈性參數(shù)Es.形如式(5)和式(6)，統(tǒng)計二者的聯(lián)合分布P(Ew,Es)并最終獲得其條件概率P(Ew|Es).結合P(F|Ew)計算地震尺度彈性參數(shù)Es條件下的巖相概率P(F|Es)(式(8))及其熵函數(shù)值Epy(F|Es)(圖4).

巖相概率P(F|Es)在P(F|Ew)基礎上考慮了井震尺度匹配P(Ew|Es)不確定性的影響，圖4中巖相概率較圖3更接近概率均值線.熵函數(shù)值Epy(F|Es)比Epy(F|Ew)顯著增加，反演最大后驗概率解與地震尺度巖相序列更為接近.

對于巖相識別中不確定性的傳遞與累積，亦可結合其數(shù)據(jù)樣點交繪特征加以分析.在物性參數(shù)矢量空間R中，各巖相數(shù)據(jù)樣點基本都落在相應的概率等值線橢圓內(nèi)，沒有明顯混疊現(xiàn)象(圖2a).在測井彈性參數(shù)矢量空間Ew中，各巖相樣點間混疊現(xiàn)象有所增加(圖5a—c).導致這種現(xiàn)象有兩個主要原因：一是巖石物理響應關系的復雜性，不同物性參數(shù)組合R的巖石物理響應Ew可能極其近似；二是數(shù)據(jù)噪聲等因素的影響.這些因素增加了基于測井彈性參數(shù)Ew巖相識別的不確定性.地震尺度彈性參數(shù)矢量空間Es中由于受到井震尺度匹配的不確定性影響，混疊現(xiàn)象更加嚴重(圖5d—f).

圖3 基于測井彈性參數(shù)的巖相識別(a) 巖相序列; (b) 孔隙度一維邊緣概率分布； (c) 泥質含量一維邊緣概率分布; (d) 巖相概率; (e) 最大后驗概率解; (f) 熵函數(shù)值.Fig.3 Facies classification conditioned on elastic data at well-log scale(a) Facies sequence; Unidimensional marginal probability distribution of porosity (b) and clay content (c); (d) Facies probability; (e) Maximum a posteriori; (f) Entropy.

圖4 基于地震尺度彈性參數(shù)Es的巖相識別(a) 巖相序列； (b) 地震尺度巖相序列； (c)巖相概率P(F|Es)； (d) 最大后驗概率解； (e) 熵函數(shù)值.Fig.4 Facies classification conditioned on elastic data at seismic scale Es(a) Facies sequence; (b) Upscaled facies sequence; (c) Facies probability P(F|Es); (d) Maximum a posteriori; (e) Entropy.

圖5 彈性參數(shù)交繪分析虛線為參數(shù)空間邊界，(a—c) 為Ew數(shù)據(jù)交繪，(d—f) 為Es數(shù)據(jù)交繪.Fig.5 Cross-plot of elastic propertiesDashed line is the boundary of parameter space (a—c) Ew data; (d—f) Es data.

井震尺度匹配環(huán)節(jié)不確定性的本質是由于地震資料低分辨率所導致的信息流失.因此對工區(qū)地震資料開展高分辨率處理是降低該環(huán)節(jié)不確定性的有力措施.其次，應選擇適用于工區(qū)資料的尺度粗化算法以及精準的粗化時窗寬度以最大程度地優(yōu)化P(Ew|Es)的統(tǒng)計樣本，從而提高統(tǒng)計精度降低不確定性.

3.3 基于屬性交繪特征約束參數(shù)空間

在構建條件概率P(F|R)、P(R|Ew)和P(Ew|Es)的過程中，由于參數(shù)矢量R、Ew和Es均為多維變量，例如Ew=[vp,vs,ρ]，上述概率矩陣規(guī)模相對較大.以二維概率矩陣P(Ew|Es)為例，其矩陣規(guī)模取決于三維矢量Ew和Es的網(wǎng)格剖分密度.P(Ew|Es)任一維度規(guī)模為三個參量網(wǎng)格剖分長度之積，總體規(guī)模達到六個參量網(wǎng)格剖分長度的乘積.而統(tǒng)計運算量與其矩陣規(guī)模成正相關，因此合理限制概率矩陣的規(guī)模是降低運算成本的重要途徑.

結合圖2和圖5可知，參數(shù)矢量R、Ew和Es參數(shù)樣點都具有一定的交繪特征.根據(jù)其交繪特征限制參數(shù)空間，可降低概率矩陣規(guī)模以減少運算量.圖2a中物性參數(shù)矢量R樣點呈帶狀分布，并沒有完全分布于整個參數(shù)空間.因此根據(jù)交繪特征將R的參數(shù)空間限制在圖中兩條虛線內(nèi)部.根據(jù)彈性參數(shù)Ew和Es中屬性的兩兩交繪特征，分別在圖5(a—c)和圖5(d—f)的二維空間中曲線內(nèi)約束其參數(shù)空間.最終三維矢量Ew和Es的參數(shù)空間分別受到上述曲線在三維空間中所代表的6個曲面約束，極大地壓縮了其參數(shù)空間范圍.對各個參數(shù)矢量(R、Ew和Es)網(wǎng)格剖分，并以網(wǎng)格節(jié)點數(shù)表征參數(shù)空間大小，約束前后參數(shù)空間范圍參見表2.

采用屬性交繪特征約束后，R和Ew的參數(shù)空間分別為全空間的35.54%和14.79%,故P(F|R)與P(R|Ew)的矩陣規(guī)模分別被壓縮為初始規(guī)模的35.54%和5.26%.同理,Es的參數(shù)空間為全空間的11.83%，因此P(Ew|Es)矩陣規(guī)模僅占初始規(guī)模的1.75%.最終極大地降低了概率矩陣P(F|Es)的矩陣規(guī)模和統(tǒng)計運算量，為后續(xù)的地震巖相識別提供了高效的運算途徑；同時剔除了不符合工區(qū)儲層地球物理特征的參數(shù)組合，限制了統(tǒng)計概率的表征范圍，從而提高了統(tǒng)計概率信息強度，降低了巖相識別的不確定性.

表2 各屬性參數(shù)空間約束對比分析Table 2 Comparative analysis of parameter space constraints for each property

3.4 地震反演不確定性影響評價

地震反演往往具有較大的不確定性，影響地震巖相識別精度.為表征地震反演的不確定性對巖相識別的影響，研究采用貝葉斯線性反演理論(Buland and Omre, 2003)由角度道集數(shù)據(jù)S(圖1)反演獲得彈性參數(shù)Es的后驗概率分布P(Es|S)(圖6).通過全貝葉斯分類方法聯(lián)合地震反演結果P(Es|S)與條件概率P(F|Es)獲取地震巖相識別的巖相概率P(F|S)及其熵值Epy(F|S)(圖7).全貝葉斯分類方法考慮了變量之間的相關性，能夠精確地傳遞不確定性信息.

由于考慮了地震反演的不確定性，圖7中巖相概率曲線較圖4和圖3更接近概率均值線；同時概率曲線由于子波濾波作用更為平滑.熵函數(shù)值Epy(F|S)對比Epy(F|Es)明顯增大，體現(xiàn)了地震反演不確定性對巖相識別的具體影響.圖中三類熵函數(shù)值逐漸增加的現(xiàn)象也體現(xiàn)了地震巖相識別不確定性的構成與傳遞特征.

地震噪聲是引起地震反演環(huán)節(jié)不確定性的主要因素，在實際生產(chǎn)中對地震數(shù)據(jù)開展去噪處理提高地震資料信噪比能夠有效降低該環(huán)節(jié)不確定性.其次，提高反演初始模型與地震子波的精度對抑制地震反演不確定性也具有重要意義.

3.5 地震噪聲不確定性影響評價

在實際生產(chǎn)中地震噪聲對巖相識別結果精度也有很大影響.為了定量評價地震噪聲對于巖相識別的影響，研究采用不同信噪比(SNR=1,2,3,5,10)的合成地震數(shù)據(jù)開展地震巖相識別，并定量對比其反演結果精度及不確定性信息，如圖8和圖9所示.

由圖8可知，地震巖相識別的不確定性與地震數(shù)據(jù)噪聲水平成正相關.在2.2、2.4 s附近位置，可明顯觀察到從左到右隨著模擬資料信噪比增加熵值逐步降低.熵函數(shù)清晰地表征了地震噪聲對于巖相識別不確定性的影響.

重建率表征反演巖相與真實巖相的匹配程度，取值范圍為[0,1].重建率越大反演結果與實際資料的匹配程度越高.熵值表征反演的不確定性，重建率描述反演的精度.圖9中，隨著巖相識別各環(huán)節(jié)引入不確定性因素增多，對應熵函數(shù)均值呈增大趨勢而重建率逐步下降，反演不確定性地增加導致反演精度下降.同時，熵函數(shù)均值變化幅度遠大于重建率的變化幅度.這表明對于不確定性信息的表征，熵函數(shù)比重建率更加敏感.圖9中展現(xiàn)了地震巖相識別中各環(huán)節(jié)不確定性的組成及其傳遞特征.通過引入熵函數(shù)，方法實現(xiàn)了地震巖相識別不確定性評價的定量化，從而進一步促進了不確定性分析的系統(tǒng)化以及不確定性傳播的可視化，為后期油藏建模與儲層評價風險評估提供重要指導信息.

圖7 地震巖相識別結果(a) 巖相序列； (b) 地震尺度巖相序列； (c) 巖相概率P(F|S)； (d) 最大后驗概率解； (e) 熵函數(shù)值.Fig.7 Results of seismic facies classification(a) Facies sequence; (b) Upscaled facies sequence; (c) Facies probability P(F|S); (d) Maximum a posteriori; (e) Entropy.

圖8 不同信噪比的地震數(shù)據(jù)巖相識別(a) 最大后驗概率解； (b) 熵函數(shù)值.Fig.8 Seismic facies classification by synthetic seismic data with different SNR(a) Maximum a posteriori; (b) Entropy.

圖9 不同條件數(shù)據(jù)約束下的重建率(虛線)和熵函數(shù)均值(實線)Fig.9 Reconstruction rate (dashed line) and the mean value of entropy (full line) conditioned on different properties

4 結論與建議

研究引入熵函數(shù)全面系統(tǒng)地實現(xiàn)了地震巖相識別不確定性定量評價，提出了結合屬性交繪特征約束反演參數(shù)空間，提高計算效率.研究結果表明：

(1) 與傳統(tǒng)的不確定性定性評價方法相比，引入熵函數(shù)實現(xiàn)了地震巖相識別不確定性構成特征和傳播規(guī)律的定量評價分析，促進了不確定性分析的系統(tǒng)化以及不確定性傳播的可視化，為后期油藏建模與儲層評價風險評估提供了重要指導信息.

(2)隨著不確定性的增加，巖相概率趨近概率均值，熵值增大，最大后驗概率解的重建率降低，熵值與重建率呈負相關；對不確定性信息的表征，熵函數(shù)比重建率更加敏感.

(3)地震資料信噪比越低，地震巖相識別的不確定性越大，熵值增大，重建率降低.

(4)結合屬性交繪特征約束反演參數(shù)空間，極大地減少了統(tǒng)計運算量，提高運算效率；并剔除了不符合工區(qū)地球物理特征的參數(shù)組合，降低了不確定性.

(5)研究方法不僅適用于模擬數(shù)據(jù)中的砂泥巖儲層類型，對于儲層巖相間具有地球物理響應差異的儲層類型均具有適用性.

Avseth P, Mukerji T, Mavko G. 2005. Quantitative Seismic Interpretation. Cambridge: Cambridge University Press.

Bachrach R. 2006. Joint estimation of porosity and saturation using stochastic rock-physics modeling.Geophysics, 71(5): O53-O63.Buland A, Kolbj?rnsen O, Hauge R, et al. 2008. Bayesian lithology and fluid prediction from seismic prestack data.Geophysics, 73(3): C13-C21.

Buland A, Omre H. 2003. Bayesian linearized AVO inversion.Geophysics, 68(1): 185-198.

Cao D P, Yin X Y, Zhang F C, et al. 2009. A study on the method of joint inversion of multiscale seismic data.ChineseJournalofGeophysics(in Chinese), 52(4): 1059-1067, doi: 10.3969/j.issn.0001-5733.2009.04.023.

Grana D, Mukerji T, Dvorkin J, et al. 2012a. Stochastic inversion of facies from seismic data based on sequential simulations and probability perturbation method.Geophysics, 77(4): M53-M72.

Grana D, Paparozzi E, Mancini S, et al. 2013. Seismic driven probabilistic classification of reservoir facies for static reservoir modelling: a case history in the Barents Sea.GeophysicalProspecting, 61(3): 613-629.Grana D, Pirrone M, Mukerji T. 2012b. Quantitative log interpretation and uncertainty propagation of petrophysical properties and facies classification from rock-physics modeling and formation evaluation analysis.Geophysics, 77(3): WA45-WA63.

Grana D, Rossa E D. 2010. Probabilistic petrophysical properties estimation integrating statistical rock physics with seismic inversion.Geophysics, 75(3): O21-O37.

Larsen A L, Ulvmoen M, Omre H, et al. 2006. Bayesian lithology/fluid prediction and simulation on the basis of a Markov-chain prior model.Geophysics, 71(5): R69-R78.

Mukerji T, J?rstad A, Avseth P, et al. 2001. Mapping lithofacies and pore-fluid probabilities in a North Sea reservoir: Seismic inversions and statistical rock physics.Geophysics, 66(4): 988-1001.

Papoulis A. 1984. Probability, Random Variables and Stochastic Processes. New York: McGraw-Hill.

Pirrone M, Grana D, Maffioletti F, et al. 2011. Monte carlo simulation and uncertainty evaluation of log-facies classification based on petroelastic properties. ∥73rd EAGE Conference & Exhibition.

Shannon C E. 1948. A mathematical theory of communication.BellSystemTechnicalJournal, 27: 379-423.

Tarantola A. 2005. Inverse problem theory. ∥ SIAM.

Teng L, Cheng J B. 2014. Bayesian petrophysical-properties estimation with AVA inversion and statistical rock physics.OilGeophysicalProspecting(in Chinese), 49(4): 729-738.

Tian Y K, Zhou H, Yuan S Y. 2013. Lithologic discrimination method based on Markov random-filed.ChineseJournalofGeophysics(in Chinese), 56(4): 1360-1368, doi: 10.6038/cjg20130430.

Yang K, Ai D F, Gen J H. 2012. A new geostatistics inversion and reservoir modeling technique constrained by well-log, crosshole and surface seismic data.ChineseJournalofGeophysics(in Chinese), 55(8): 2695-2704, doi: 10.6038/j.issn.00015733.2012.08.022.

Zhang F C, Xiao Z B, Yin X Y. 2014. Bayesian stochastic inversion constrained by seismic data.OilGeophysicalProspecting(in Chinese), 49(1): 176-182.

Zhang G Z, Wang D Y, Yin X Y, et al. 2011. Study on prestack seismic inversion using Markov Chain Monte Carlo.Chinese

JournalofGeophysics(in Chinese), 54(11): 2926-2932, doi: 10.3969/j.issn.0001-5733.2011.11.022.

附中文參考文獻

曹丹平, 印興耀, 張繁昌等. 2009. 多尺度地震資料聯(lián)合反演方法研究. 地球物理學報, 52(4): 1059-1067, doi: 10.3969/j.issn.0001-5733.2009.04.023.

騰龍, 程玖兵. 2014. 貝葉斯框架下聯(lián)合AVA反演與統(tǒng)計巖石物理的儲層參數(shù)估計. 石油地球物理勘探, 49(4): 729-738.

田玉昆, 周輝, 袁三一. 2013. 基于馬爾科夫隨機場的巖性識別方法. 地球物理學報, 56(4): 1360-1368, doi: 10.6038/cjg20130430.

楊鍇, 艾迪飛, 耿建華. 2012. 測井、井間地震與地面地震數(shù)據(jù)聯(lián)合約束下的地質統(tǒng)計學隨機建模方法研究. 地球物理學報, 55(8): 2695-2704, doi: 10.6038/j.issn.00015733.2012.08.022.

張繁昌, 肖張波, 印興耀. 2014. 地震數(shù)據(jù)約束下的貝葉斯隨機反演. 石油地球物理勘探, 49(1): 176-182.

張廣智, 王丹陽, 印興耀等. 2011. 基于MCMC的疊前地震反演方法研究. 地球物理學報, 54(11): 2926-2932, doi: 10.3969/j.issn.0001-5733.2011.11.022.

(本文編輯 張正峰)

Quantitative evaluation of uncertainties in seismic facies identification based on probabilistic statistics

YUAN Cheng1,2, LI Jing-Ye1,2*, CHEN Xiao-Hong1,2

1StateKeyLaboratoryofPetroleumResourcesandProspecting,ChinaUniversityofPetroleum,Beijing102249,China2NationalEngineeringLaboratoryforOffshoreOilExploration,ChinaUniversityofPetroleum,Beijing102249,China

In the early stage of oilfield exploration, reservoir characterization is always considered to be a risky task, since few well-log data are available at the foremost stage of any production plan in the petroleum industry as well as the ambiguity of seismic data. Uncertainty in seismic reservoir characterization is generally quite large, especially in seismic facies classification which goes through multiple links. Therefore, quantitative uncertainty evaluation is valuable for seismic facies classification. It provides an important guiding sense for risk management as well as decision-making in any petroleum reservoir.For evaluating the uncertainty propagation in seismic facies classification quantitatively, the uncertainty of well-log facies definition, rock physics modeling, scale change and seismic inversion has been taken into consideration in this case. We firstly compute the facies probabilities conditioned on different properties in each step of seismic facies classification. Then, the associated uncertainty and maximum a posteriori (MAP) of facies probabilities are assessed by means of entropy and reconstruction rate respectively, since the variable that represents facies is categorical. The influence of seismic noise on facies classification is also analyzed by synthetic seismic data with different signal noise ratio. In addition, the parameter spaces of well-log and upscaled elastic properties are restricted by the data distribution characters in cross-plot.Synthetic example shows that uncertainty in seismic facies classification could be evaluated quantitatively with this methodology. By introducing entropy, the constitution and propagation of uncertainty can be evaluated quantitatively with the help of facies probability. The total flow chart of uncertainty evaluation in seismic facies classification is connected systematically. Furthermore, the influence of uncertainty propagation and seismic noise on facies classification is illustrated visually by entropy and reconstruction rate. Restriction of parameter spaces by the data distribution characters in cross-plot can not only reduce the computational cost but also the uncertainty in seismic facies classification dramatically since the parameter vectors which fall out of the restricted scopes are precluded.Quantitative uncertainty evaluation brings many details of the uncertainty in seismic facies classification. It makes it possible to assess the propagation and accumulation of uncertainty quantitatively as well as its influence on the result accuracy, since the increment of entropy is considered to be the incoming uncertainty of the current step. It provides us with a unique angle of view to understand the uncertainty in seismic facies classification as well as the great value for risk management and optimal decision-making in the petroleum industry.

Probability statistics； Entropy； Facies identification； Uncertainty; Quantitative evaluation

10.6038/cjg20151032.

Yuan C, Li J Y, Chen X H. 2015. Quantitative evaluation of uncertainties in seismic facies identification based on probabilistic statistics.ChineseJ.Geophys. (in Chinese),58(10):3825-3836,doi:10.6038/cjg20151032.

國家自然科學基金項目(U1262207)，國家科技重大專項課題(2011ZX05019-006)和中國石油大學(北京)?；?2462012KYJJ0508)聯(lián)合資助.

袁成，男，1988年生，在讀博士生，主要從事巖相識別與地質統(tǒng)計學研究.E-mail:yuancheng0124@139.com

*通訊作者 李景葉，男，1978年生，教授，博導，2005年獲得中國石油大學(北京)博士學位，主要從事儲層地球物理學研究.E-mail:ljy3605@sina.com

10.6038/cjg20151032

P631

2014-10-14，2015-07-12收修定稿

袁成, 李景葉, 陳小宏. 2015. 基于概率統(tǒng)計的地震巖相識別不確定性定量評價方法.地球物理學報，58(10):3825-3836,