基于FILTERSIM算法的油藏沉積微相模擬研究

孟欣然，梁堰波，孟憲海，楊 欽

(1.北京航空航天大學軟件開發環境國家重點實驗室，北京100191;2.北京航空航天大學計算機學院，北京100191)

0 引言

地質統計學建模方法從傳統上分為兩個流派:基于象元的兩點地質統計學方法和基于目標的地質統計學方法［10］。第一種方法是以變差函數為工具，利用兩點變差函數來反映儲層空間結構，能夠很好地忠實于實際數據的約束，但很難表征復雜的幾何形態 (如彎曲河道)和空間結構，因此更適用于海相而非陸地相的地質結構。而第二種方法則更符合地質工作者的經驗，但是較難與實際數據相結合。多點地質統計學 (Multiple-point Statistics，以下簡稱MPS)［1］模擬方法結合了基于象元的模擬方法和基于對象的模擬方法兩者的優點。首先，多點地質統計學是一種改進的基于象元的方法，繼承了基于象元的模擬方法的靈活性，每一步只對一個像素點進行建模，因此很容易滿足控制點的要求;同時，由于引入了訓練圖像，使得建模結果的曲線連續性有了很大提高，尤其擅長模擬復雜的河道地質結構，更符合實際的需求。

MPS第一個比較重要的算法是SNESIM算法［2］，此處不再贅述。SNESIM算法只能用于油藏沉積微相的建模，而不能用于數值連續的油藏屬性建模，同時，如果沉積相類型的數目過多，SNESIM算法的效率不甚理想。當訓練圖像較復雜或者分類較多時，SNESIM算法對內存的需求也是巨大的。因此，2006年，一種新的 MPS算法——FILTERSIM［3］誕生了，該算法能夠在合理的內存需求下完成沉積微相建模和屬性建模。FILTERSIM算法引入了一系列不同方向的線性濾波器 (linear filter)把訓練圖像進行分類，類似的訓練圖像模式被歸為一類，并以這些模式的均值作為這一類的標志，稱為這個類的“原型” (prototype)，在模擬過程中，原型與當前約束數據最相近的那個類被選為當前網格點的實現。SNESIM算法把所有備選的訓練圖像都存儲在一個搜索樹里，因此對內存要求較高，而FILTERSIM算法只存儲每個訓練圖像的中心點，因此，極大的降低了內存需求。

同時，油藏最重要的屬性參數，如孔隙率、滲透率、含油飽和度等都是連續數據，SNESIM算法無法實現模擬連續數據的需求，而FILTERSIM算法可以分為FILTERSIMCATE和FILTERSIM_CONT，后者可以對連續數據進行精確模擬，彌補了SNESIM算法的不足。

1 FILTERSIM算法概要

FILTERSIM算法主要分為3步:濾波分 (filter score)計算、模式分類和模式模擬。此外，對于已經進行過地震數據處理的油藏，FILTERSIM算法可以綜合地震處理數據進行上述3個過程。

1.1 計算濾波分

FILTERSIM算法一個非常重要的概念是搜索模板(search template)τn，用來搜索訓練圖像。搜索模板的大小用3個整數表示(nx，ny，nz)來表示，并且三個參數必須都是奇數。搜索模板的每一個節點都記錄一個其與中心節點的偏移量，ua=u+ha(a=1，2，…，n)。如圖1所示，ha代表了 (a)中除去中心點u之外的其他80個點較中心點的偏移向量、(b)中除去中心點u之外26個點較中心點u的偏移向量。將搜索模板應用于訓練圖像就獲得了一個圖像模式集。圖2顯示了將搜索模板應用于訓練圖像一個點獲得的一個圖像模式的過程。

FILTERSIM算法里的濾波器 (filter)定義為一系列數據的集合TJ，集合的元數為 J，TJ={u0;hj，j=1，…，J}，集合中的每個節點uj都用與中心節點u0的偏移向量hi=(x，y，z)i來表示，同時賦予該節點一個特定的加權量(weight)fj。偏移向量的元素x、y、z都是整數。為了簡化存儲，FILTERSIM提出了一個濾波分 (filter score)ST(u)的概念，公式是

其中，pat(u+hi)是節點的屬性值，J=nx×ny×nz。一套濾波器是無法捕捉訓練圖像的全部特征的，因此，一個濾波器集合F共同用于發掘訓練圖像的特征。上面的公式可以改寫為

因此，訓練圖像的存儲規模由nx×ny×nz縮減為F。比如說，一個11×11×3像素的3D訓練圖像可以用FILTERSIM算法生成的9個F值來表示。如圖3所示。

圖3 利用濾波器簡化過程

對于連續的訓練圖像 (訓練圖像像素的值代表屬性值)來說，F個濾波分直接用來代表訓練圖像元素本身。對于分類型訓練圖像 (像素值只有幾種類型，代表不同的沉積微相)，假設沉積相種類為K，訓練圖像首先要轉化為一系列二進制指示向量 Ik(u)，k=0，…，K －1，u∈Ti

圖4給出了二維圖形常用的六個濾波器，公式如下:

(a)南北方向均值，中心行權值較重

(b)東西方向均值，中心列權值較重

(c)南北方向梯度，邊緣行權值較重

(d)東西方向梯度，邊緣列權值較重

(e)南北方向曲率

(f)東西方向曲率

圖4 適用于二維訓練圖像的濾波器組合

1.2 圖像模式分類

將F個濾波器應用于K個沉積相的訓練圖像可以獲得F×K個濾波分圖，這個維度遠遠小于圖像模式的像素數，大大減小了維度。相似的圖像模式會得到差不多的濾波分，因此，可以將濾波分進行分類，相似的濾波分為一類，以均值作為這一類的原型 (prototype)，記為prot，公式如下

分類策略最常用的有兩類，一類是交叉分類 (cross partition)，一類是K-均值聚類 (K-mean cluster)，此處不詳細闡述。前者的優點是效率較快，但是分類結果較粗糙。后者是最簡單的非監督學習模型，能夠獲得較好的分類效果。

1.3 圖像模式模擬

圖像模擬過程就是不斷尋找與當前點最相近的模式的過程，具體算法如下:

算法1

(1)利用一定的濾波器獲得濾波分

(2)在濾波分空間 (F×K維)進行模式分類

(3)在模擬網格G中定位控制數據

(4)定義一個訪問G的隨機路徑

(5)for路徑中的每一個節點u，do

(6)提取以u為中心點的約束數據集dev

(7)找到與dev最相近的父原型protp

(8)if protp存在子原型集，then

(9)找到與dev最接近的子原型protc

(10)從protc中隨機選擇一個模式pat

(11)else

(12)從protp中隨機選擇一個模式pat

(13)end if

(14)將pat黏貼到已經模擬的點集中，并凍結dev中的所有點

(15)end for

2 利用FILTERSIM算法進行油藏沉積微相建模

由于SNESIM算法在解決油藏沉積微相方面的能力很強 (前提是沉積相種類較少)，為了評估FILTERSIM算法在沉積微相建模方面的性能，我們對Stanford V數據分別使用SNESIM算法和FILTERSIM進行了沉積微相模擬。如圖5所示，Stanford V數據分為4個部分，其中grid_2d是一個用于進行屬性建模的100*130*1笛卡爾網格，并且提供了三個地震屬性數據，P(crevasse|seis)2、P(mud|seis)2和P(sand|seis)2代表了3種沉積相的分布概率。TI_2D_large提供了算法需要的訓練圖像，well提供約束的井數據。圖6是Standford V數據中的訓練圖像和地震數據。

圖6中，(a)是二維的三相訓練圖像，其中淺色代表砂質相 (sand channel)，深色代表泥相 (mud face)，小點代表斷層 (crevasse);(b)是目標的泥相概率分布圖像。

在本實例中，地震處理提供的3種沉積相分布被用于距離計算，前述的算法可以被改寫成如下過程 (其中tau模型是一個合成概率場的工具):

算法2

(1)對于隨機路徑上的每一個節點u，同時從以實現的模擬數據和地震數據中獲得約束數據集dev和sdev(后者來源于地震數據)

(2)if dev是空集then

(3)用sdev代替dev，并利用新的dev獲取最近的模式

(4)else

(5)用dev獲取最近的模式prot

(6)用tau模型將模式prot和地震數據sdev合成一個新控制數據集dev*，作為本地概率分布模型

(7)尋找與dev*最接近的模式，繼續模擬過程

(8)end if

搜索模板的大小為11×11×1(搜索模板大小可以自行設置，取決于對精度的需求，一般為訓練圖像的1/10-1/5)。圖7是分別利用SNESIM算法和FILTERSIM算法獲得100組結果選取的三個結果。

圖7 利用FILTERSIM算法和SNESIM算法實現二維三相沉積微相分類結果

圖7 中，(a)、(c)、(e)是由FILTERSIM算法實現的3個沉積微相分類結果;(b)、(d)、(f)是由SNESIM算法實現的3個沉積微相分類結果。

圖8是利用FILTERSIM算法和SNESIM算法實現的沉積微相模擬結果生成的mud概率分布圖像?？梢钥闯鯢ILTERSIM算法獲得的概率模型邊緣比較模糊，而SNESIM算法獲得的較清晰，但是兩者都較好地反映了地震數據揭示的地址概率模型。

圖8 利用FILTERSIM算法和SNESIM算法實現的沉積微相模擬結果生成的泥相概率分布圖像

圖8 中，(a)是由100個FILTERSIM算法模擬結果生成的泥相概率分布;(b)是由100個SNESIM算法模擬結果生成的泥相概率分布。

3 利用FILTERSIM算法模擬連續屬性數據

為了驗證FILTERSIM算法在連續屬性數據模擬方面的性能，我們在Stanford VI［4］數據集上進行了實驗。本實驗的訓練圖像是一個80×80×40像素的三維圖像，像素值的范圍是0-1，代表油藏的孔隙度，如圖9所示?？紫抖染禐?.6522，方差為0.078。抽取其中的0.5%數據作為約束數據，如圖10所示。利用連續FILTERSIM算法獲得的模擬結果之一如圖11所示。表1是隨機選取的三個結果中孔隙度均值和方差數據?？梢钥闯?，無論是圖形模式、均值還是方差，FILTERSIM的模擬結果都與訓練圖像十分接近。

圖11 連續FILTERSIM算法獲得的模擬結果之一

表1 FILTERSIM算法隨機選取的3個結果中均值和方差數據

4 結束語

FILTERSIM算法能夠在較少的CPU和內存消耗下獲得與SNESIM算法近似的結果，并且能處理沉積相建模和屬性建模，這是SNESIM算法無法實現的。然而，在沉積相建模質量方面來看，FILTERSIM算法略遜于SNESIM算法，尤其是在邊緣輪廓方面比較模糊。這在井間距離較大或者井數據較少的地質條件下尤其有用。此外，FILTERSIM還能用于連續屬性數據模擬，并且模擬結果能較好地反映訓練圖像的模式 (均值和方差均與訓練圖像十分接近)。

［1］ZHANG T F.Incorporating geological conceptual models and interpretations into reservoir modeling using multiple-point geo-statistics［J］.Earth Science Frontiers，2008，15(1):26-35.

［2］Boucher A.Considering complex training images with search tree partitioning［J］.Computers ＆ Geosciences，2009，35(6):1151-1158.

［3］ZHANG T F.Filter-based training pattern classification for spatial pattern simulation［D］.Ph D dissertation，USA:Stanford University，2006:11-26.

［4］Remy N，Boucher A，WU JB.Applied geostatistics with SGeMS:A Users'guide［M］.Cambridge University Press， 2009:108-134.

［5］WU J B.4D seismic and multiple-point pattern data integration using geostatistics［D］.Ph D dissertation， USA:Stanford University，2007:163-173.

［6］Chugunova T L，HU L Y.Multiple-point simulations constrained by continuous auxiliary data［J］.Mathematical Geosciences 2008，40(2):133-146.

［7］Arpat G B，Caers J.Conditional simulation with patterns［J］.Mathematical Geology，2007，39(2):177-203.

［8］Sanchez-Brea L M，Bernabeu E.Uncertainty estimation by convolution using spatial statistics［J］.IEEE Transactions on Image Processing，2006，15(10):3131-3137.

［9］LU D T，ZHANGT J，YANGQ，et al.A reconstruction method of porous media integrating soft data with hard data［J］.Chinese Science Bulletin，2009，54(11):1876-1885.

［10］LUO Yang.Using MPSin the river reservoir simulation［J］.Geological Science and Technology Information，2008，27(3):68-72(in Chinese).［駱楊.多點地質統計學在河流相儲層建模中的應用［J］.地質科技情報，2008，27(3):68-72.］