基于連通性函數(shù)的多點地質(zhì)統(tǒng)計建模參數(shù)敏感性分析方法

杜 川, 喻思羽, 李少華, 方 紅

(1.長江大學(xué) 地球科學(xué)學(xué)院,武漢 430100；2.中國石油 天然氣遼河油田分公司 勘探開發(fā)研究院,盤錦 124010)

0 引言

隨著多點地質(zhì)統(tǒng)計學(xué)(multi-point statistics， MPS)受到廣泛的關(guān)注，涌現(xiàn)出許多相對成熟的MPS建模算法。相對于傳統(tǒng)的兩點地質(zhì)統(tǒng)計學(xué)，多點地質(zhì)統(tǒng)計學(xué)利用訓(xùn)練圖像代替變差函數(shù)表達地質(zhì)變量的空間結(jié)構(gòu)性[1]，從中獲取更高階的樣本統(tǒng)計量[2],MPS建模優(yōu)點在于既能夠模擬地質(zhì)體復(fù)雜的幾何結(jié)構(gòu)，也能忠實于條件數(shù)據(jù)。Guardiano等[2]提出的首個MPS算法每模擬一個網(wǎng)格節(jié)點都需要遍歷整個訓(xùn)練圖像，計算成本非常高，使得該算法難以實際應(yīng)用。之后，更好的算法被開發(fā)出來改善這些缺點。多點統(tǒng)計技術(shù)包括基于象元的建模算法SNESIM[3]、基于樣式的建模算法FILTERSIM[4]、SIMPAT[5]、DISPAT[6]、基于直接采樣的算法DS[7]、IMPACA[8]、HOSIM[9]、CCSIM[10]、基于圖像縫合技術(shù)的Image Quilting[11]。這些MPS算法使用不同的方法來讀取訓(xùn)練圖像、存儲信息和復(fù)制模式。Mariethoz[12]對這些多點地質(zhì)建模方法進行了定性比較。多點地質(zhì)統(tǒng)計學(xué)建模已用于許多應(yīng)用之中，例如常見的在儲層地質(zhì)建模研究領(lǐng)域，利用Simpat模擬河流相儲層分布[13]；多點地質(zhì)統(tǒng)計學(xué)在河流相儲層建模中的應(yīng)用[14]；在三角洲相儲層建模領(lǐng)域的應(yīng)用，如多點地質(zhì)統(tǒng)計學(xué)建模方法在復(fù)雜疊置樣式砂體表征中的應(yīng)用[15]。

在兩點地質(zhì)統(tǒng)計學(xué)建模中，通過變差函數(shù)考慮統(tǒng)計數(shù)據(jù)之間的空間相關(guān)性[16]。類似的在多點地質(zhì)統(tǒng)計學(xué)建模過程，需要通過給定的數(shù)據(jù)事件對訓(xùn)練圖像進行掃描來獲取一個統(tǒng)計先驗?zāi)Ｐ停嚓P(guān)的統(tǒng)計數(shù)據(jù)取決于給定數(shù)據(jù)事件的參數(shù)，通過人為的參數(shù)調(diào)整，找到合適的參數(shù)設(shè)置是獲取良好建模效果的前提，不同的參數(shù)不僅影響建模的質(zhì)量，也涉及CPU及內(nèi)存開銷問題[17]。為此前人也進行了大量的探究，常見的參數(shù)敏感性分析，如Liu[18]提供了關(guān)于SNESIM參數(shù)以及每個參數(shù)對模型再現(xiàn)質(zhì)量和CPU成本的完整指南；Meerschman[19]也為DS直接采樣算法做了類似的工作；Honarkhah[6]提出了DISPAT，它對用戶提供的參數(shù)不太敏感。在該方法中，通過計算不同樣板尺寸的圖案的平均熵，然后使用所謂的肘形圖的拐點，自動選擇最佳樣板尺寸。Strebelle[20]提出了一種通過優(yōu)化數(shù)據(jù)樣板的大小來最小化SNESIM的計算成本，同時保持模式再現(xiàn)質(zhì)量的方法。在該方法中，最小可接受的數(shù)據(jù)樣板大小是通過閾值來選擇的，超過該閾值，額外的條件數(shù)據(jù)不會提高估計的條件概率。Kolbj?rnsen[21]開發(fā)了一種自動方法來定量確定馬爾可夫網(wǎng)格建模中使用的多個網(wǎng)格的最佳數(shù)量，該方法使用從訓(xùn)練圖像獲得的方向相關(guān)函數(shù);Bai[22]還提出了一種量化方法來估計SNESIM算法的多個網(wǎng)格的最小數(shù)量。它使用由改進的聯(lián)合計數(shù)統(tǒng)計量化的空間關(guān)聯(lián)度來尋找最粗略的尺度，并因此基于模擬尺度應(yīng)不小于目標尺度的標準來估計多個網(wǎng)格的最小數(shù)量。

傳統(tǒng)建模參數(shù)的評價主要依靠人工識別，即首先給定一組(有序)參數(shù)集，使用該參數(shù)集里的每個參數(shù)模擬一組隨機模型，然后通過人工判別優(yōu)選與訓(xùn)練圖像相似度較為合適的模型，以該模型對應(yīng)的參數(shù)作為最優(yōu)參數(shù)。人工識別的精確性取決于建模工作者的經(jīng)驗，具有較強主觀性，同時人工識別效率低，不適應(yīng)于現(xiàn)代自動化生產(chǎn)的需要。為了解決這一問題，我們引入了一種基于連通性函數(shù)[23-24]的多點地質(zhì)統(tǒng)計學(xué)建模參數(shù)敏感性分析方法，通過計算條件數(shù)據(jù)下生成的模擬圖像的連通性概率曲線，并把訓(xùn)練圖像與模擬圖像連通性概率曲線之間的相似性用余弦相似度來量化，得到基于連通性函數(shù)的空間相關(guān)性評價指標與建模參數(shù)的關(guān)系曲線，通過對量化數(shù)據(jù)進行分析，找到相似度曲線趨于平緩的拐點值，以此拐點值作為最優(yōu)參數(shù)，大于最優(yōu)值參數(shù)將變得不敏感，模型效果不會隨著參數(shù)值增加而得到明顯改善。

1 連通性計算原理

連通性函數(shù)被定義為網(wǎng)格上的兩個點屬于同一個連接區(qū)域的概率[24]，利用連通性函數(shù)計算區(qū)域化變量指示圖的連通性，可以反映區(qū)域化變量在某一方向某一距離上連通性概率的變化程度。理論上隨機模型與訓(xùn)練圖像的連通性差異越小，則表明二者沉積相的空間相關(guān)性及位置分布特征越相近，也說明了用于生成該隨機模型的建模參數(shù)越可靠。Eulogio[25]詳細地描述了連通性概率計算原理，其計算原理如下：

圖1定義工區(qū)網(wǎng)格為矩形G，G的網(wǎng)格節(jié)點(或象元)屬于子集S(砂巖相)或補集Sc(泥巖相)

G=SUSc

(1)

G的象元取值函數(shù)定義為指示函數(shù)I(u)

(2)

如果u和u′=u+h屬于相同的S，用u?u′表示，連通性函數(shù)τ(h)定義為屬于S的兩點連通的概率：

τ(h)=P(u?u+h│u,u+h∈S)

(3)

其近似計算公式為式(4)。

(4)

以砂巖相連通性為例，如圖1所示u和u′有三種情況，①u和u+h有一個或者兩個都屬于泥巖相，例如u3與u3+h或者u4與u4+h，則不參與砂巖相連通計算；②u和u+h都是砂巖，但不屬于相同砂巖體，例如u2與u2+h，則增加N(u,u+h∈S)數(shù)值；③u和u+h都是砂巖，屬于相同砂巖體，例如u1與u1+h則增加#N(u?u+h│u,u+h∈S)和#N(u,u+h∈S)值。

圖1 連通性示意圖Fig.1 Schematic diagram of connectivity

2 參數(shù)敏感性分析方法

2.1 實現(xiàn)步驟

基于連通性函數(shù)的多點地質(zhì)統(tǒng)計學(xué)建模參數(shù)優(yōu)選方法是建立在一個前提條件下：以經(jīng)典多點算法Snesim的一個重要建模參數(shù)--“搜索節(jié)點數(shù)”為例，基于連通性函數(shù)優(yōu)選搜索節(jié)點數(shù)；隨著搜索節(jié)點的增加，隨機模型與訓(xùn)練圖像的空間相關(guān)性、幾何形態(tài)越來越相似。

在此前提條件下采用有序排列的建模參數(shù)集建立一組隨機模型集，然后計算基于連通性函數(shù)和余弦相似度CosSim函數(shù),評價建模參數(shù)集的多點地質(zhì)統(tǒng)計隨機模型與訓(xùn)練圖像的空間相關(guān)性及結(jié)構(gòu)特征相似性，進而建立基于連通性函數(shù)的空間相關(guān)性評價指標與建模參數(shù)的關(guān)系曲線，選取評價指標開始平穩(wěn)、進入平臺區(qū)域時的拐點所對應(yīng)的參數(shù)值作為最優(yōu)參數(shù)。

以建模參數(shù)“搜索節(jié)點數(shù)”為例，基于連通性函數(shù)的建模參數(shù)優(yōu)選方法具體流程及流程圖(圖2)如下：

圖2 優(yōu)選方法流程圖Fig.2 Flow chart of optimization method

1)輸入訓(xùn)練圖像TI，給定MPS的關(guān)鍵建模參數(shù)——搜索節(jié)點的K個連續(xù)單調(diào)遞增數(shù)值(例如：1，2，…，K)組成的集合S={1，2，…，K}。

2)采用MPS，使用TI和S生成隨機模型集M，其中第k個參數(shù)Sk相應(yīng)的隨機模型集定義為Mk，每個Mk包含n個隨機模型，則M中共計有N=K×n個模型。

(5)

(6)

計算Mk所有模型與訓(xùn)練圖像TI關(guān)于連通性函數(shù)參數(shù)組的余弦相似度，然后取其數(shù)學(xué)期望作為建模參數(shù)Sk對應(yīng)的空間相關(guān)性評價指標ξk

ξk=E[cosθk(Ai,Bi)]

(7)

5)k增加1，如果k小于K-1，則進入步驟4)，否則進入步驟6)。

6)建立參數(shù)k與相似度的關(guān)系曲線，分析曲線變化趨勢，把曲線變化趨于平緩時的拐點建模參數(shù)k(i)作為最優(yōu)參數(shù)。隨著建模參數(shù)值的遞增，基于建模參數(shù)值的隨機模型與訓(xùn)練圖像之間相似度逐步增加，并最終達到難以區(qū)分的程度，即距離也來越小，相似度越來越高，同時對應(yīng)的參數(shù)越不敏感。

2.2 Snesim建模參數(shù)敏感性分析實例

訓(xùn)練圖像是地質(zhì)先驗?zāi)Ｐ停錁邮健⒋笮〉榷紩δM結(jié)果產(chǎn)生重大影響，以一個網(wǎng)狀河訓(xùn)練圖像channel為例，圖3訓(xùn)練圖像大小為101×101的網(wǎng)格，黑色部分代表河道砂，白色部分代表背景泥，訓(xùn)練圖像數(shù)字化的結(jié)果是一個101×101的矩陣，“1”代表河道，“0”代表背景泥。其中背景泥占比0.7，河道占比0.3。

圖3 訓(xùn)練圖像channelFig.3 Training image channel

2.2.1 搜索節(jié)點數(shù)

利用Snesim基于像元的多點地質(zhì)建模對訓(xùn)練圖像進行非條件模擬，為了找出最優(yōu)建模參數(shù)值，這個過程中其他參數(shù)保持不變，只對Snesim重要建模參數(shù)“搜索節(jié)點”數(shù)進行更改，設(shè)置搜索節(jié)點數(shù)5到75，步長10。其他參數(shù)分別設(shè)置為：分類變量比例0.7、0.3，伺服器系統(tǒng)參數(shù)0.5;搜索半徑hx=50，hy=50，搜索角度angel_1=0, angel_2=0，多重網(wǎng)格層數(shù)為3。對應(yīng)每個搜素節(jié)點數(shù)下進行模擬生成300個隨機模型。圖4是不同搜索節(jié)點數(shù)所模擬出來的一個模型(300個實現(xiàn)模型中的一個)，隨著搜索節(jié)點數(shù)量的增加，所模擬出來的河道形態(tài)越來越好，與訓(xùn)練圖像越來越相似，當增加到25和35以后，設(shè)置不同搜索節(jié)點數(shù)所模擬出來的模型的空間特征相似用肉眼就無法直接區(qū)分開來，同時隨著搜索節(jié)點增加，該參數(shù)敏感性降低，所生成的隨機模型效果改善不明顯甚至幾乎難以區(qū)分，往往在考慮計算效率以及實現(xiàn)效果上，通過人為肉眼觀察給定一個最優(yōu)參數(shù)值往往具有誤差，會給建模中帶來不確定性。

圖4 不同搜索節(jié)點數(shù)對應(yīng)的實現(xiàn)模型Fig.4 Implementation models corresponding to different node numbers(a)Node=5;(b)Node=15;(c)Node=25;(d)Node=35;(e)Node=45;(f)Node=55;(g)Node=65;(h)Node=75

根據(jù)圖2提出的優(yōu)選方法流程，計算隨機模型與訓(xùn)練圖像砂巖泥巖相在x、y方向上不同距離下的連通概率，基于連通性概率函數(shù)計算公式(4)計算訓(xùn)練圖像channel砂巖相泥巖相在x、y方向上的連通概率，圖5是訓(xùn)練圖像在x和y方向上砂巖泥巖在不同滯后距下的連通概率，隨著滯后距的增加砂巖泥巖相在x、y方向上的連通性總體上不斷趨于減小，其中砂巖的變化有周期化的上下波動，并且其變化規(guī)律能大致反映砂巖的空間分布特征。

圖5 訓(xùn)練圖像channel在x、y方向上的連通性概率曲線Fig.5 Connectivity probability curve of training image channel in x and y directions

不同搜索節(jié)點所實現(xiàn)的隨機模型與訓(xùn)練圖像在空間型態(tài)結(jié)構(gòu)特征上存在差異，隨著搜索節(jié)點數(shù)量的增加，隨機模型的形態(tài)結(jié)構(gòu)等特征與訓(xùn)練圖像越來越相似，差異會越來越小，這種差異也同時體現(xiàn)在了連通性概率曲線上如圖6所示，其中mod1代表搜索節(jié)點數(shù)等于5的隨機模型(綠色)，Mod2是搜索節(jié)點數(shù)等于75的隨機模型(藍色)，實線與虛線分別代表泥巖相與砂巖相在x或y方向上的連通概率，由圖6可知，模型與訓(xùn)練圖像的形態(tài)特征和空間結(jié)構(gòu)越相似，二者之間的連通性概率曲線就越相近相似。(圖6曲線紅色和藍色相近，紅色與綠色差異大)

圖6 訓(xùn)練圖像及其隨機模型的連通性概率曲線對比圖Fig.6 Comparison of the connectivity probability curves of the training image and its random model(a)x方向；(a)y方向

利用余弦相似度cosSimilarity函數(shù)計算訓(xùn)練圖像與隨機模型連通性概率曲線的相似度，得到隨機模型(搜索節(jié)點數(shù)從5～75)與訓(xùn)練圖像關(guān)于連通性概率的余弦相似度，然后取其數(shù)學(xué)期望作為建模參數(shù)“搜索節(jié)點”對應(yīng)的空間相關(guān)性評價指標ξk建立相似度與搜索節(jié)點數(shù)的關(guān)系曲線(圖7)，隨著搜索節(jié)點數(shù)從5開始增加，300個隨機模型與訓(xùn)練圖像的相似度期望不斷增加，表明模型與訓(xùn)練圖像越來越相似(根據(jù)余弦相似度定義，余弦值越接近1，表示越相似)，當搜索節(jié)點數(shù)達到35時，空間相關(guān)性評價指標進入一個平臺區(qū)域，參數(shù)變得不敏感，模型質(zhì)量將不再有明顯提升，是本例中的優(yōu)選參數(shù)。

圖7 建模參數(shù)(搜索節(jié)點數(shù))基于連通性函數(shù)的余弦相似度變化曲線Fig.7 Cosine similarity change curve based on connectivity function for modeling parameters (node numbers)

2.2.2 多重網(wǎng)格層數(shù)

多重網(wǎng)格層數(shù)對模擬的效果影響很大，但是隨之而來的是消耗更大的計算性能，利用上述分析方法，可以找到一個最優(yōu)值，既能達到良好建模效果，同時也節(jié)約計算性能。對“多重網(wǎng)格重數(shù)”這一參數(shù)利用相同方法測試，設(shè)置多重網(wǎng)格重數(shù)從1到8。

保持其他參數(shù)不變設(shè)置搜索節(jié)點數(shù)35、分類變量比例0.7、0.3、伺服系統(tǒng)參數(shù)0.5，搜索半徑hx=50，hy=50、搜索角度angel_1=0, angel_2=0。利用snesim多點地質(zhì)建模對訓(xùn)練圖像進行非條件模擬，在每個多重網(wǎng)格重數(shù)下進行模擬生成300個隨機模型(圖8顯示對應(yīng)網(wǎng)格層數(shù)下所模擬的300個實現(xiàn)模型中的一個)。

圖8 不同網(wǎng)格重數(shù)對應(yīng)的實現(xiàn)模型Fig.8 Implementation models corresponding to different multigrids(a)Multigrids=1;(b)Multigrids=2;(c)Multigrids=3;(d)Multigrids=4;(e)Multigrids=5;(f)Multigrids=6;(g)Multigrids=7;(h)Multigrids=8

利用余弦相似度cosSimilarity函數(shù)計算訓(xùn)練圖像與隨機模型連通性概率曲線的相似度，得到隨機模型(多重網(wǎng)格重數(shù)1-8)與訓(xùn)練圖像關(guān)于連通性概率的余弦相似度，然后取其相似度均值建立其對應(yīng)多重網(wǎng)格重數(shù)變化的曲線(圖 9)。

圖9 建模參數(shù)(網(wǎng)格重數(shù))基于連通性函數(shù)的余弦相似度變化曲線Fig.9 Cosine similarity change curve based on connectivity function for modeling parameters (multigrids)

隨著多重網(wǎng)格數(shù)從1開始增加，300個隨機模型與訓(xùn)練圖像的相似度期望不斷增加，表明模型與訓(xùn)練圖像越來越相似(根據(jù)余弦相似度定義，余弦值越接近1，表示越相似)，當網(wǎng)格重數(shù)到3時，空間相關(guān)性評價指標進入一個平臺區(qū)域，模型質(zhì)量將不再有明顯提升，參數(shù)敏感性降低。這是因為當網(wǎng)格很粗時，相應(yīng)的用來掃描訓(xùn)練圖像的樣板規(guī)模會大于訓(xùn)練圖像的規(guī)模，導(dǎo)致大的數(shù)據(jù)事件很難被找到。此時數(shù)據(jù)事件就會相應(yīng)縮小，而這就等同于使用較小的數(shù)據(jù)樣板進行掃描，一般來說，網(wǎng)格重數(shù)設(shè)置多少取決于最大尺度結(jié)構(gòu)信息的大小，在本例中101×101的網(wǎng)格訓(xùn)練圖像中，利用連通性方法定量化地給出最佳網(wǎng)格重數(shù)為3，超過3以后，模擬效果不會有明顯改善，該實驗結(jié)果也符合前人所提出的認識。

3 結(jié)論

采用基于連通性函數(shù)的多點地質(zhì)統(tǒng)計學(xué)建模參數(shù)敏感性分析方法，分別對兩個關(guān)鍵建模參數(shù)(搜索節(jié)點數(shù)、多重網(wǎng)格層數(shù))進行敏感性測試，得出以下認識：

1)隨著搜索節(jié)點數(shù)和多重網(wǎng)格層數(shù)增加，模型與訓(xùn)練圖像的形態(tài)視覺特征越來越相似，利用連通性函數(shù)來量化這一空間相關(guān)性特征，結(jié)合余弦相似度度量方法對其進行定量的分析，找出價指標開始平穩(wěn)、進入平臺區(qū)域時的拐點所對應(yīng)的參數(shù)值作為最佳參數(shù)。

2)針對本例101×101的河道訓(xùn)練圖像，在snesim模擬過程中搜索節(jié)點數(shù)不能過少，過少捕捉不到河道的結(jié)構(gòu)特征，過多并不能明顯改善建模效果，參數(shù)在到達一定值后變得不敏感，反而會增加計算量，降低效率，多重網(wǎng)格層數(shù)也是如此。利用連通性的方法定量化的分析參數(shù)敏感性，找到搜索節(jié)點數(shù)設(shè)置在35左右，多重網(wǎng)格層數(shù)設(shè)置為3層，能獲取良好模擬效果，也能減少不必要的運算，提高效率。