尹舒祚， 喻思羽， 李少華， 杜 川， 金 浩

(長江大學(xué)地球科學(xué)學(xué)院， 武漢 430100)

1962年，法國數(shù)學(xué)家Matheron[1]提出了區(qū)域化變量的概念，創(chuàng)立了地質(zhì)統(tǒng)計學(xué)。傳統(tǒng)的地質(zhì)統(tǒng)計建模方法[2]，如克里金法使用變差函數(shù)對空間場進行估值，但是克里金法作為確定性建模方法，對具有結(jié)構(gòu)性和隨機性的復(fù)雜地質(zhì)特征表達存在一定的局限性。斯坦福大學(xué)Journel將隨機性特征引入地質(zhì)統(tǒng)計學(xué)儲層建模，提出了隨機模擬方法。1992年多點地質(zhì)統(tǒng)計學(xué)(multi-point geostatistical simulation, MPS)被提出，并且以2002年提出的Snesim算法為關(guān)鍵點，迅速成為隨機建模的研究熱點[3]。使用“訓(xùn)練圖像”的多點地質(zhì)統(tǒng)計建模方法可很好地表達復(fù)雜地質(zhì)變量的空間結(jié)構(gòu)性，綜合了基于象元和基于目標方法的優(yōu)點[4]，克服了傳統(tǒng)基于變差函數(shù)地質(zhì)統(tǒng)計學(xué)在刻畫儲層形態(tài)結(jié)構(gòu)特征方面的不足，也可以很好忠于硬數(shù)據(jù)，展現(xiàn)了良好的應(yīng)用前景。目前多點地質(zhì)統(tǒng)計學(xué)的建模方法有迭代法和非迭代法兩類[5]，迭代法存在難以滿足收斂條件、計算速度過于緩慢等問題，導(dǎo)致難以在實際建模過程中得到應(yīng)用，故逐漸被非迭代法取代。非迭代法包含了Strebelle[6]提出的基于概率的Snesim算法，Arpat[7]提出的基于數(shù)據(jù)樣式的Simpat算法，以及之后發(fā)展的Filtersim[8]、DS-MPS[9]等多種多點地質(zhì)統(tǒng)計算法。

多點地質(zhì)統(tǒng)計建模方法已經(jīng)在眾多領(lǐng)域得到應(yīng)用，研究學(xué)者[10-14]應(yīng)用不同的多點算法對三角洲前緣微相、河流相儲層分布、儲層沉積相及數(shù)字巖心等領(lǐng)域進行地質(zhì)模型的建立和探討，證實了多點地質(zhì)統(tǒng)計建模方法的實用性和可行性。使用MPS的一個重要問題是建模質(zhì)量不僅取決于所使用的訓(xùn)練圖像，還取決于一系列算法參數(shù)[15]。建模參數(shù)作用于空間多點數(shù)據(jù)事件的提取與重建完整流程，因此選用建模參數(shù)是否合適不僅影響建模效果，而且存在建模效率及內(nèi)存開銷[16-17]問題。為此，中外相關(guān)學(xué)者做了大量的研究，包括以彎曲河道作為訓(xùn)練圖像，利用Snesim算法對重要的建模參數(shù)(條件數(shù)據(jù)點最大數(shù)量、搜索鄰域、多重網(wǎng)格等)進行了敏感性分析[18-20]，并對如何合理設(shè)置相關(guān)參數(shù)提出了建議。Straubhaar等[21]在Impala中用列表結(jié)構(gòu)及索引樹取代搜索樹，實現(xiàn)了算法計算并行化，提高了模擬效率；喻思羽等[22]通過深入剖析Simpat算法原理及特點，針對其在平衡計算效率和內(nèi)存空間問題上存在的不足，提出了基于樣式降維聚類的多點地質(zhì)統(tǒng)計建模算法；Yang[23]和Honarkhah等[24]對最佳模板大小，利用熵提出了不同的模板優(yōu)選方法。Meerschman等[25]對并查集(disjoint-set，DS)算法的綜合敏感性分析，提供了DS參數(shù)設(shè)置的指南；Baninajar等[26]提出了MPS自動參數(shù)優(yōu)化器，這是一種基于隨機優(yōu)化的通用方法，可以快速逼近任何MPS方法和不同類型設(shè)置的最優(yōu)參數(shù)。Bai等[27]提出了一種估算Snesim最小多重網(wǎng)格數(shù)量的定量方法，使用空間關(guān)聯(lián)度來分析尋找訓(xùn)練圖像具有統(tǒng)計意義的最大規(guī)模的結(jié)構(gòu)信息，以此來估計多重網(wǎng)格的合適取值。對建模者在訓(xùn)練圖像優(yōu)選上存在的主觀標準，近幾年有學(xué)者[28-30]提出了客觀的訓(xùn)練圖像評估和優(yōu)選方法。

每個MPS算法都有一組特定的參數(shù)，這些參數(shù)設(shè)置對多點地質(zhì)建模中模擬效果的好壞至關(guān)重要，需在建模之前進行選擇。雖然MPS算法在最佳參數(shù)選擇上已有一些篩選方法，但包含了人為主觀因素，缺乏一定的準確客觀評價原則。到目前為止，尋找最佳參數(shù)的常見做法是進行靈敏度分析，但是該方法過于煩瑣。為了解決這一問題，現(xiàn)提出基于平均熵值的多點地質(zhì)統(tǒng)計建模參數(shù)優(yōu)選方法。在地質(zhì)統(tǒng)計學(xué)中，熵是一種隨機性的統(tǒng)計度量，可用來表征訓(xùn)練圖像的紋理，通過計算不同尺寸模式分解下每個模型的平均熵序列，并將訓(xùn)練圖像和模擬圖像的平均熵曲線差異用Hsim相似度進行量化表示，得到基于平均熵值的空間相關(guān)性評價指標與建模參數(shù)的關(guān)系曲線，通過對量化數(shù)據(jù)進行分析，找到Hsim相似度曲線趨于平穩(wěn)的拐點值，將該值作為最佳模擬參數(shù)。

1 基于平均熵值的多點地質(zhì)統(tǒng)計建模參數(shù)優(yōu)選方法

1.1 熵值計算原理

熵的概念由德國物理學(xué)家Clausius于1865年提出，最初是用來描述“能量退化”的物質(zhì)狀態(tài)參數(shù)之一。隨著統(tǒng)計物理、信息論等一系列學(xué)科的發(fā)展和完善，熵在不同領(lǐng)域得到了重要應(yīng)用。在地質(zhì)統(tǒng)計學(xué)中，熵作為一種隨機性的統(tǒng)計度量，可用來表征訓(xùn)練圖像的紋理。根據(jù)香農(nóng)第三定理——信源編碼定理，香農(nóng)的熵是衡量消息中所包含的信息，將該概念應(yīng)用于訓(xùn)練圖像上，可以確定所需的最小信息，以便可靠地表示由二進制數(shù)字編碼的模式。所以熵在模板大小優(yōu)選方面，通過靜態(tài)增加訓(xùn)練圖像中的模板大小，可以觀察到熵值會先急劇增加(編碼或壓縮圖像的平均信息數(shù)在增加)，之后隨著模板大小增加到表示訓(xùn)練圖像平穩(wěn)特征的最佳大小，熵值會以很慢的速度趨于平穩(wěn)(編碼具有一些平穩(wěn)特征的圖像所需信息近似訓(xùn)練圖像本身的平穩(wěn)性特征)。

熵在模板大小選擇中，假設(shè)某一特定對象(如泥巖相)在訓(xùn)練圖像上大量重復(fù)，最佳模板大小將是捕獲泥巖對象的模板大小。當(dāng)模板大小逐漸增加時，會有如下情況：①模板大小低于最佳值，提取模式中包含的信息量將增加；②模板大小達到最佳值，作為泥巖對象的訓(xùn)練圖像特征被很好地編碼在圖像中；③模板大小超過最佳模板，提取模式中的信息量不再增加。

根據(jù)上述概念，隨著模板大小增加，以熵值為特征的信息量逐漸增加，當(dāng)信息量增加到與訓(xùn)練圖像相似時，熵值趨于穩(wěn)定。尺寸大小為nx×ny×nz的圖像熵值計算公式為

(1)

式(1)中：K為隨機變量的可能結(jié)果數(shù)；pi為概率質(zhì)量函數(shù)(即直方圖)。熵值越高，意味著更多的隨機性。

1.2 實現(xiàn)步驟

對于任意隨機模擬方法而言，建模參數(shù)值發(fā)生遞增或者遞減變化時，生成隨機模型的表征變化屬于漸變，而不是突變。提出的建模參數(shù)優(yōu)選方法建立在此前提條件下，以多點地質(zhì)建模參數(shù)——數(shù)據(jù)樣板尺寸為例，隨著數(shù)據(jù)樣板尺寸逐漸增大，生成隨機模型和訓(xùn)練圖像的空間相關(guān)性和幾何形態(tài)越來越相似。

采用有序遞增排列的建模參數(shù)集模擬一組隨機模型集，計算在不同尺寸模式分解方式下訓(xùn)練圖像和所有隨機模型的平均熵，并對相同建模參數(shù)值實現(xiàn)的100個隨機模型模式均熵求均值，利用基于相似性度量的高維數(shù)據(jù)聚類Hsim算法描述隨機模型與訓(xùn)練圖像(training image, TI)的平均熵曲線差異，選取擬合曲線上開始趨于平穩(wěn)的拐點所對應(yīng)參數(shù)值作為優(yōu)選參數(shù)。以建模參數(shù)“數(shù)據(jù)模板尺寸”為例，圖1為基于平均熵值的建模參數(shù)優(yōu)選方法具體流程，具體步驟如下。

步驟一輸入訓(xùn)練圖像TI，給定MPS關(guān)鍵參數(shù)——數(shù)據(jù)模板尺寸的集合C={1,2,…,K},集合C為連續(xù)單調(diào)遞增數(shù)值，K表示數(shù)據(jù)模板大小。

步驟二采用MPS，使用TI和C生成隨機模型集M，其中第k個參數(shù)Ck對應(yīng)隨機模型集定義為Mk，每個Mk中包含m個隨機模型，模型集M中共計有K×m個隨機模型。

步驟三選取參數(shù)Ck的模型集Mk和TI，計算不同尺寸模式分解下的平均熵序列，計算公式為

(2)

步驟四計算M模型集K組模型的模式平均熵與訓(xùn)練圖像的平均熵曲線差異，依次比對模型和訓(xùn)練圖像平均熵值曲線和Hsim相似性，空間相似性計算公式為

(3)

高維數(shù)據(jù)相似性度量函數(shù)Hsim()計算公式為

(4)

式(4)中：S(X,Y)為對象X和Y在空間上的相似程度，其值越大，表明數(shù)據(jù)之間的相似性程度越高；d為空間維度。

步驟五k每次增加1，若k≤K-1，則進入步驟三、步驟四，否則進入步驟六。

圖1 優(yōu)選方法流程圖Fig.1 Flow chart of the preferred method

步驟六建立參數(shù)k與相似度的擬合曲線，分析并選取曲線變化趨于平穩(wěn)時的拐點參數(shù)值k作為最優(yōu)參數(shù)。隨著建模參數(shù)值的遞增，基于建模參數(shù)值的隨機模型與TI相似度逐漸增加，并達到難以人為區(qū)分的程度，即當(dāng)建模參數(shù)大于最佳值后，隨機模型的空間結(jié)構(gòu)特征變化趨于穩(wěn)定。

2 Simpat建模參數(shù)優(yōu)選實例

多點地質(zhì)統(tǒng)計儲層隨機建模方法依賴于訓(xùn)練圖像，訓(xùn)練圖像作為地質(zhì)先驗?zāi)Ｐ?，其質(zhì)量的好壞對模擬結(jié)果至關(guān)重要。使用基于樣式的多點地質(zhì)統(tǒng)計算法Simpat以channel模型作為訓(xùn)練圖像進行非條件模擬，對Simpat的重要建模參數(shù)——數(shù)據(jù)模板尺寸開展研究。

圖2 channel訓(xùn)練圖像Fig.2 The channel training image

圖2中，channel為河流相數(shù)字模型，包含砂巖和泥巖兩種相類型，訓(xùn)練圖像大小為101×101。使用Simpat對TI進行非條件模擬，針對模擬模板尺寸參數(shù)的優(yōu)選問題，設(shè)置數(shù)據(jù)樣板尺寸大小從1遞增到15，在每個數(shù)據(jù)樣板參數(shù)下，模擬生成100個隨機模型。如圖3所示為不同模板大小生成的隨機模型，隨著模板尺寸逐漸增加，生成隨機模型的河道連續(xù)性快速改善，形態(tài)表征與TI越來越相似，當(dāng)模板大小增加到3或者4后，河道隨機模型的紋理表征改善越來越低，人為很難區(qū)分模擬結(jié)果的紋理表征，若通過人為主觀視覺來確定改善效果的臨界值，給定一個最優(yōu)參數(shù)，會增加建模過程中的不確定性。

根據(jù)圖1的參數(shù)優(yōu)選方法流程，采用式(2)計算每個隨機模型和訓(xùn)練圖像在不同尺寸分解模式下的平均熵序列，圖4為模板大小為1～15的隨機模型與TI的平均熵序列擬合曲線，隨著模板大小增加到3后，隨機模型平均熵序列和訓(xùn)練圖像擬合曲線逐漸接近，根據(jù)地質(zhì)統(tǒng)計學(xué)中熵值的概念，表示隨機生成模型的空間形態(tài)結(jié)構(gòu)達到穩(wěn)定，即模擬結(jié)果和TI越來越相似。

再求相同模板參數(shù)模擬的100個隨機模型平均熵序列均值，由式(4)計算訓(xùn)練圖像channel和隨機模型平均熵的Hsim相似度，來量化表示兩者之間的差異度。如圖5所示，隨著模擬樣板尺寸從1開始增加，隨機模型和TI的相似度逐漸增加并趨于穩(wěn)定，即100個隨機模擬結(jié)果的平均熵和TI的平均熵曲線差異越來越小并趨于穩(wěn)定。當(dāng)模擬樣板尺寸達到3時，空間相關(guān)性評價指標逐漸平穩(wěn)，模擬結(jié)果質(zhì)量不再有明顯提升，從而確定出本例最優(yōu)模板參數(shù)大小。

圖3 基于Simpat不同模板尺寸的channel隨機模型Fig.3 The channel random model based on different template sizes of Simpat

圖4 隨機模型和TI的平均熵序列擬合曲線Fig.4 Random model and TI mean entropy sequence fitting curves

圖5 基于平均熵值的Hsim相似度變化曲線(channel)Fig.5 Hsim similarity change curves based on average entropy(channel)

使用多點地質(zhì)統(tǒng)計算法Simpat對另一訓(xùn)練圖像circle進行非條件模擬，同時對Simpat的重要建模參數(shù)——數(shù)據(jù)模板尺寸開展研究。圖6為circle數(shù)字模型，包含砂巖和泥巖兩種相類型，訓(xùn)練圖像大小為101×101。

圖6 circle訓(xùn)練圖像Fig.6 The circle training image

使用Simpat對TI進行模擬，模式實現(xiàn)流程和參數(shù)設(shè)置同河道訓(xùn)練圖像，圖7為模板大小1-15生成的circle隨機模型。根據(jù)圖1參數(shù)優(yōu)選方法流程，隨機模型、熵值計算、相似度計算及量化擬合熵值差異曲線同河道訓(xùn)練圖像，得到Hsim相似度曲線如圖8所示，同channel，隨著模板大小增加，隨機模型和訓(xùn)練模型空間相似度越來越相似，由空間相似度評價指標可得模板大小3為最優(yōu)參數(shù)。

圖8 基于平均熵值的Hsim相似度變化曲線(circle)Fig.8 Hsim similarity change curve based on average entropy (circle)

圖7 基于Simpat不同模板尺寸的circle隨機模型Fig.7 The circle random model based on different template sizes of Simpat

3 結(jié)論

(1)選取最佳參數(shù)適用于離散的訓(xùn)練圖像，模式平均熵值計算采用單點計算代替兩點熵值計算，節(jié)約了計算時間,提高了參數(shù)優(yōu)選效率。

(2)結(jié)合多點地質(zhì)統(tǒng)計隨機模型空間特征和建模參數(shù)的相關(guān)性認識，建模參數(shù)值在遞增或遞減時，生成隨機模型特征屬于漸變而非突變，提出一種基于平均熵值的多點地質(zhì)統(tǒng)計建模參數(shù)優(yōu)選方法。以建模參數(shù)——模板大小為例，隨著模板尺寸增加，隨機模型與TI的空間結(jié)構(gòu)、形態(tài)特征越來越相似，利用模式均熵來量化其空間特征，結(jié)合Hsim相似度度量方法對隨機模型和TI模式均熵差異進行定量化分析，給出客觀的優(yōu)選參數(shù)值，對多點建模參數(shù)選擇和優(yōu)化建模效果具有實際意義。