一種基于MPS和ISOMAP的空間數據重建方法

杜 奕 張 挺 黃 濤

1(上海第二工業大學工學部 上海 201209)2(上海電力學院計算機科學與技術學院 上海 200090)3(中國科學技術大學近代力學系 合肥 230027)(duyi@sspu.edu.cn)

?

一種基于MPS和ISOMAP的空間數據重建方法

杜 奕1張 挺2黃 濤3

1(上海第二工業大學工學部 上海 201209)2(上海電力學院計算機科學與技術學院 上海 200090)3(中國科學技術大學近代力學系 合肥 230027)(duyi@sspu.edu.cn)

在空間數據重建過程中,條件數據對重建結果影響較大，在僅有少量條件數據的情況下,重建結果常常出現較多的不確定性,此時適合采用不確定性插值方法重建空間數據.作為目前不確定性插值的主流方法之一,多點信息統計法(multiple-point statistics, MPS)可以從訓練圖像提取模式的本質特征,然后將這些特征復制到待模擬區域.由于傳統采用線性降維的MPS方法無法有效處理非線性數據,因此將等距特征映射(isometric mapping, ISOMAP)應用到MPS方法,以實現對非線性數據的降維.提出基于MPS和ISOMAP的空間數據重建方法,通過模式庫構建、模式降維、模式分類、模式提取等步驟能夠較為準確地重構出未知的空間數據,為MPS處理非線性空間數據提供了新思路.實驗結果表明:該方法所重建的空間數據具有與訓練圖像相似的結構特征.

訓練圖像;多點信息統計法;模式;等距特征映射;降維

Fig. 1 Well data(hard conditional data) and seismic data(soft conditional data).圖1 井位數據(硬條件數據)和地震數據(軟條件數據)

空間數據是指包含空間特征的數據集合.現實生活中的空間數據,既有來自科學試驗的數據,也有來自生產實踐的數據,其格式與特征是包羅萬象[1-3].與一般數據相比,空間數據具有一定的特點,主要表現在2個方面：1)空間數據除包含屬性信息外,還含有空間信息；2)空間數據具有空間自相關性,即鄰近數據間的相關性比距離較遠的數據間的相關性要大得多[4-6].目前獲取某些領域真實有效的空間數據存在一定的難度,出現這種情況的主要原因是科學實驗難度大或勘探開發費用高.

數據插值成為重建空間數據的1個有效手段[7-8].插值方法可分為“確定”性插值方法和“不確定”性插值方法.“確定”性插值方法的插值形式、插值函數參數以及插值結果基本都是確定的[9-10].不確定性插值方法的不確定性一方面表現在選用的插值方式具有隨機性,另一方面表現在插值參數的選取和確定依賴概率統計原則[11-13].

不確定性插值方法主要包括克里金(Kriging)方法和隨機模擬方法.多點信息統計法(multiple-point statistics, MPS)是目前隨機模擬的主流.通過再現高階統計量,MPS能夠從訓練圖像中捕捉復雜的特征樣式并把它們復制到待模擬區域中,從而獲得最終模擬結果[14-15].

在MPS重建數據過程中,訓練圖像的模式特征決定了模擬結果,但這些模式往往維數較高,數據處理較為困難.因此,模式降維問題成為MPS中的研究熱點.一部分主流MPS方法,如過濾器隨機模擬(filter-based stochastic simulation, FILTERSIM)[15]和距離模式隨機模擬(distance-based pattern simulation, DISPAT)[16],將數據降維引入到其隨機模擬過程中.但是,這些方法都采用線性降維處理模式數據,而真實情況下的數據并非都以線性形式存在著.如果空間數據處于非線性狀態卻使用線性方法降維,那么將會嚴重降低隨機模擬的質量.另外,如果把線性降維當作非線性降維的一種特殊情況[17-18],那么基于空間數據非線性降維的隨機模擬具有更廣泛的一般性,可同時用于線性和非線性的空間數據重建.

條件數據對模擬結果的影響很大,而條件數據又可以分為硬條件數據和軟條件數據.硬條件數據是已知的精確采樣點數據,但可能數量很少.與硬條件數據相對應的是軟條件數據.軟條件數據總體上能夠比較準確地反映所研究變量的變化趨勢,但精度較差.在許多領域里,由于受到客觀條件或技術水平限制,所能得到的硬條件數據非常有限,但是可以獲得相對比較豐富的軟條件數據.例如在石油勘探過程中,所能獲得的硬條件數據(井位數據)往往非常少,而有關所研究對象的軟條件數據(如地質解釋和地震資料等)卻相對較為豐富.如圖1(a)表示勘探得到的井位數據,其中黑色和白色點分別表示砂巖(sand)和頁巖(shale)井位的分布,是一種硬條件數據.圖1(b)為對應區塊的地震資料,是頁巖分布概率,可作為軟條件數據.可見硬條件數據的數量相對較少,如果能充分利用較為豐富的軟條件數據,那么必然會提高所建模型的精度.

本文基于訓練圖像的非線性降維,提出一種利用多點信息統計法重建空間數據的方法,可以分別在結合使用軟硬條件數據,只有硬條件數據和無條件數據情況下實現空間數據重建.由于本方法通過捕獲訓練圖像的內在結構特征進行數據重建,因此適用于條件數據較少或沒有條件數據情況下的數據重建.雖然重建結果具有隨機性和不確定性,但是它們能在已知條件數據和先驗模型基礎上提供未來數據發展趨勢的預測,具有一定實際意義.

Fig. 2 Data template.圖2 數據模板

1 MPS中的一些基本概念

設數據模板為τD,它是由D個向量組成的幾何形態,τn={hα;α=1,2,…,D}.設模板中心位置為u,模板其他位置uα=u+hα(α=1,2,…,D).圖2(a)是一個9×9節點組成的二維模板,uα由中心點u和80個向量hα所組成,各向量用網格節點表示.而在三維空間中數據模板的定義也是適用的,圖2(b)是由3×3×3體素組成的三維模板,模板中心點u用黑色表示,其周圍的各三維節點表示各個向量的位置,用淺灰色表示.注意本文的“節點”表示二維平面的像素或者三維空間的體素.假設變量S可以取K個狀態值,即有狀態值的集合{sk,k=1,2,…,K}.令d(uα)表示在uα的狀態值.為了獲得在某個向量u的狀態值,選取離其最近的D個數據作為模擬時的條件數據.D個條件數據和它們的幾何結構組成了“數據事件”,記為dD[19]:

dD={d(uα)=skα,α=1,2,…,D,kα∈[1,K]},

(1)

式(1)代表D個向量在uα位置的S(u1),S(u2),…,S(uD)分別為某個狀態值sk.

圖3(a)(b)所示為圖2(a)的二維數據模板所捕獲的2個數據事件,數據模板中包含3種灰度的節點:深灰、淺灰和白色節點.深灰和淺灰節點表示2種已知數據，白色節點為未知數據.同理,圖3(c)(d)所示為圖2(b)的三維數據模板所捕獲的2個數據事件.其中,待模擬點u仍然位于3×3×3數據模板的中心,用黑色表示.深灰色節點表示已知數據點，淺灰色節點表示未知數據.

Fig. 3 2D data event and 3D Data event.圖3 二維數據事件和三維數據事件

2 方法核心思想

Fig. 4 Capture a pattern in a training image.圖4 捕獲訓練圖像中的一個模式

本文方法主要由4個部分組成:建立訓練圖像模式庫、模式的降維、模式的分類和模式的提取.各部分詳細內容分別介紹于2.1至2.4節.

2.1 建立訓練圖像的模式庫

利用數據模板掃描訓練圖像來捕獲數據事件,這些數據事件被稱為“模式”,而這些模式包含了訓練圖像的結構特征.例如圖4(a)是一個3×3的數據模板；圖4(b)是利用數據模板掃描一個二維訓練圖像；圖4(c)是捕獲的一個模式.如果按照每次移動一個節點,從左到右、從上到下進行掃描,那么就可以獲得訓練圖像全部的模式庫,如圖5所示.對于三維訓練圖像同理可以建立類似的模式庫.

Fig. 5 Build a pattern database from a training image.圖5 獲取訓練圖像的模式庫

設TI(u)表示以u為中心的模式,ti(u+hi)表示u+hi位置的屬性值,那么有：

TI(u)=(ti(u+h1),ti(u+h2),…,ti(u+hD)).

(2)

可見模式TI(u)和數據模板尺寸相同,如果將各個模式從訓練數據中抽取出來形成模式庫,那么與具體位置u無關,此時可設模式庫中的模式個數為Npat,那么第k個模式可以表示為

patk=(patk(h1),patk(h2),…,patk(hD)),

k=1,2,…,Npat,

(3)

其中,patk(hi)與ti(u+hi)(i=1,2,…,D)一一對應.整個模式庫patDb可以表示為

patDb=(pat1,pat2,…,patNpat),

(4)

其中，每個patk(k=1,2,…,Npat)可以視為模式空間中的一個點.

2.2 模式的降維

降維方法可分為線性降維和非線性降維.線性降維方法主要用來對線性數據降維,代表性方法有主成分分析(principal component analysis, PCA)和多維尺度變換(multi-dimensional scaling, MDS)等.然而現實世界的許多數據處于非線性結構,采用線性降維方法難以處理非線性數據.本文采用等距特征映射(isometric mapping, ISOMAP)實現模式的非線性降維[20-21].

每個模式patk(k=1,2,…,Npat)含有D個分量,可以視其維數為D.本方法通過ISOMAP將patDb中的各個模式patk降維到d(d?D)維空間.降維后得到的低維數據集設為Y=(y1,y2,…,yNpat)(其中yk∈d,k=1,2,…,Npat),即樣本個數不變,每個樣本的維數從D變為d.利用ISOMAP降維的步驟如下[21]:

1) 構造模式近鄰圖G.計算模式庫patDb中所有模式對(pati,patj)間的歐氏距離,記為dE(pati,patj)：

(5)

比較歐氏距離,確定每個模式的k個最相似模式；構造模式近鄰圖G,圖G中的節點和各模式一一對應.邊代表近鄰關系,即如果pati和patj是相似模式,G中對應節點i和j用邊連接,否則斷開.邊權值用模式間歐氏距離dE(pati,patj)表示.

2) 計算模式間的最短路徑dM(pati,patj).當圖G中存在邊(i,j)時,設最短路徑dM(pati,patj)=dE(pati,patj),否則dM(pati,patj)=∞,對于p=1,2,…,Npat：

dM(pati,patj)=min{dM(pati,patj),

dM(pati,patp)+dM(patp,patj)}.

(6)

最短路徑距離陣可表示為DM=(dM(pati,patj))2,i,j=1,2,…,Npat.

3) 計算低維映射.計算對稱陣:

(7)

其中，I是Npat階單位陣,l是元素全為1的Npat維列向量.低維空間的維度為d,則求取矩陣B的d個最大的特征值λ1≥λ2≥…≥λd,對應的特征向量是α1,α2,…,αd,則可得降維后的各個模式：

Y=(y1,y2,…,yNpat)T=

(8)

至此獲得patDb={patk}(patk∈D,k=1,2,…,Npat)的低維映射Y={yk}(yk∈d,k=1,2,…,Npat),數據從高維空間D映射到低維空間d.低維空間的維數d可以通過最大似然估計法獲得,具體步驟可見文獻[22].

2.3 模式的分類

在完成訓練圖像中模式的降維后,需要對獲取的低維模式進行分類.對這些模式采用一種基于密度的聚類方法(density-based spatial clustering of applications with noise, DBSCAN)進行分類[23].其主要思想是：只要鄰近區域的密度(對象的數目)超過某個閾值,就繼續執行聚類操作.

在利用DBSCAN劃分模式空間后,獲得若干“子空間”,每個子空間命名為Cell.圖6所示是一個二維空間被劃分為若干個Cell,其中包含一個空的Cell,每個黑點表示一個降維后的模式.對于每個Cell,可以定義一個與數據模板相同形狀的“平均模板”,稱為Prototype. Prototype是屬于該Cell的所有模式在原高維空間中各節點位置的均值.例如,圖7(a)是一幅訓練圖像(150×150像素),包含頁巖和砂巖2種狀態值,圖7(b)表示屬于一個Cell的所有15個模式(13×13像素),圖7(c)表示該Cell對應的Prototype(13×13像素).本質上Prototype可以視為每個Cell中所有模式的均值,是該Cell所有模式的“代表”.

Fig. 6 Illustration of separating pattern space into some Cells.圖6 劃分模式空間成各個Cell

Fig. 7 Illustration of a Prototype.圖7 Prototype示意圖

Fig. 8 A patch and its inner part and outer part.圖8 patch及其inner和outer部分

Prototype的值prot(l)(hα)表示屬于一個Cell中的所有圖案在各個hα位置的均值,定義為

(9)

protl=(prot(l)(h1),prot(l)(h2),…,prot(l)(hD)),

(10)

其中，l=1,2,…,L,protl可以視為第l個模式類的平均模式.

2.4 模式的提取

利用掃描訓練圖像時使用的數據模板掃描待重建區域,檢索以未知點u為中心的數據事件內的所有已知節點,這些點視為硬條件數據.設已知硬條件數據點數目為n′≤D.提取的訓練圖像的模式(稱為patch)本質上即訓練圖像內的圖像塊.將patch分為2個部分：inner部分和outer部分.它們可以作為其他后繼節點重建時的硬條件數據.不同之處在于,inner部分會被“固定”為重建結果,不會再對其進行模擬；而patch的outer部分不僅作為其他節點模擬的條件數據,而且該部分的所有節點會被當作未知點重新模擬.如圖8(a)是一個patch,圖8(b)(c)分別是其inner部分(虛線內部分)和outer部分(虛線外部分).如果一個patch的inner部分被設定得比較大,那么重建速度會較快,但是容易產生局部的不連續性；反之,較小的inner部分會導致重建速度比較慢,但是重建效果會較好.

設t表示數據類型,則patch中一共可能有3種數據：

1)t=1,原始硬條件數據,它們被分配到各節點位置上；

2)t=2,已經模擬過的節點,它們被固定為硬條件數據,即patch的inner部分(不包含inner部分內的原始硬條件數據)；

3)t=3,通過“粘貼”Cell中的patch而已知的節點,但這些點會被重新模擬,即patch的outer部分(不包含outer部分內的原始硬條件數據).

另外定義“距離函數”如下：

(11)

其中，t=1,2,3;l=1,2,…,L；dis(d(uα),protl)表示求取數據事件d(uα)和protl中對應的已知節點間的距離.每種節點會根據其類型而給定一個權值ω(t),表示其在求取距離函數中的重要性,注意有ω(3)≤ω(2)≤ω(1).可見原始硬條件數據點在距離函數中的作用最大.

假設以u為中心的硬條件數據事件hd(uα)與各個protl之間的距離可以表示為

(12)

(13)

如果在位置u處也對應存在一個以u為中心的軟條件數據事件sd(uα),它與各個protl之間的距離可以表示為

(14)

(15)

由式(15)可得整合軟硬條件數據的數學模型,即求取軟硬條件數據情況下數據事件與protl之間的距離向量：

Distotal=perhd×Dishd+μ×Dissd,

(16)

其中,perhd表示在數據模板中已知硬條件數據的節點數占模板總節點數的百分比,因為當硬條件數據較多時,其能提供的輔助判斷信息會增加,權值相應變大；而μ(0≤μ≤1)是表示軟數據準確性的權值,它的取值取決于人們對軟數據準確性的主觀意見.當μ=0,表示完全忽略軟數據；而當μ=1,表示完全信任軟數據.另外,Dishd和Dissd要進行歸一化處理,即令其取值范圍為[0,1].

假設有:

(17)

則Distotal可以表示為

(18)

3 重建過程的完整流程

根據初始條件數據的不同,節點u的重建可以分為4種情況:1)無條件數據模擬(稱為情況C1);2)只有硬條件數據情況下的模擬(稱為情況C2);3)軟條件數據和硬條件數據同時存在時的模擬(稱為情況C3);4)只有軟條件數據情況下的模擬(稱為情況C4).第2節對本方法各核心部分作了重點介紹,本節將第2節部分整合起來,形成本方法的完整流程:

1) 利用數據模板掃描訓練圖像，構建模式庫patDb.

2) 利用ISOMAP對模式庫中的模式進行降維.

3) 采用DBSCAN對降維后的模式進行分類.

4) 定義一條隨機訪問路徑，對重建區域內的未知節點進行訪問.

5) 檢查訪問路徑上的待模擬節點u是已知的硬條件數據還是已模擬節點.如果是已知的硬數據或已模擬節點，則對隨機路徑上的下個節點進行判定；否則轉向步驟6).

6) 檢索以u為中心的數據模板內的硬條件數據(設數目為n′)和對應的軟條件數據.根據條件數據的不同,提取patch的方法分為4種情況:

① 如果條件數據屬于情況C1(即n′=0,無軟條件數據)，那么隨機地從訓練圖像模式庫patDb中隨機提取一個patch，然后將該patch直接復制到以當前模擬節點u為中心的重建區域.

② 如果條件數據屬于情況C2(即n′>0,無軟條件數據)，即以u為中心的硬數據事件非空，那么就在訓練圖像所有的protl中利用式(16)尋找與當前數據事件最接近的protl.一旦該protl被確定,則從與之對應的Cell中隨機提取一個patch,然后將它復制到以當前模擬節點u為中心的重建區域.注意在計算式(16)時,只能獲得硬距離Dishd,軟距離Dissd不存在.

③ 如果條件數據屬于情況C3(即n′>0,且有軟條件數據)，那么同時有軟硬條件數據,此時直接利用式(16)獲得與當前數據事件最接近的protl,然后從與該protl對應的Cell中隨機提取一個patch復制到重建區域.

④ 如果條件數據屬于情況C4(即n′=0,僅有軟條件數據)，也是利用式(16)獲得與當前數據事件最接近的protl,然后從與該protl對應的Cell中隨機提取一個patch復制到重建區域.在重建的初始階段時,式(16)中僅存在軟距離Dissd,但是隨著重建過程的進行,已模擬節點不斷增加,這些節點成為后繼節點模擬的硬條件數據.那么此時變為軟硬條件數據結合情況下的重建.

7) 將patch復制到重建區域后,再將patch的部分節點固定為inner部分,該patch的其他節點作為outer部分.

8) 重復步驟5)～7)，繼續對其他節點進行模擬，直到隨機路徑上的所有節點被模擬完畢.

本方法的流程圖如圖9所示.值得注意的是,由于情況C3下硬條件數據通常較為稀疏,即在大多數情況下重建區域內只有軟條件數據,因此本文將情況C3和C4視為同一種情況,并且在實驗部分僅對情況C1,C2,C3進行測試.

Fig. 9 Flow chart of the proposed method.圖9 本文方法流程圖

4 實驗結果及分析

4.1 重建質量的評價標準

對于評價重建結果優劣的標準,隨機重建與一般重建方法有較大的不同.一般重建方法中的重建數據與原始數據的相似度是由兩者在每個相同位置具有相同狀態值元素的比例決定的.在一般重建方法中,該比例越高,重建相似度就越高,但是這個判定標準對于隨機重建并不適用.隨機重建并不要求隨機模擬結果與訓練圖像完全一致,它的真正作用體現在捕捉訓練圖像的結構特征,并將這些特征復制到重建數據中.隨機重建可以給出多個模擬結果,這些結果是對訓練圖像結構特征的反映.重建效果的優劣在于這些結構特征是否被復制到重建數據中.

因此,本文評價重建結果時并不直接比較重建結果和訓練圖像是否相同,而是通過多個重建結果的均值和方差分布以及變差函數來衡量.變差函數能夠反映變量在某個方向上空間結構變化的相關性和變異性.如果2幅圖像中的某個屬性值在同一個方向上具有相似的變差函數曲線,那么可以說明這2幅圖像中的該屬性值在此方向上具有相似的結構特征；反之,該屬性值在此方向上的結構差異較大.變差函數[24]的定義為

(19)

其中,Z(u)和Z(u+h)分別表示位置u和u+h處的變量狀態值,h表示距離向量,E表示求數學期望值.

4.2 重建實驗

為了驗證本方法有效性,分別對三維離散型數據和三維連續型數據進行重建實驗.所有實驗都是在Intel core i3-4160T (3.1 GHz CPU)和4 GB DDR3內存環境下運行的.式(11)中的權值ω(1)=0.5,ω(2)=0.3,ω(3)=0.2.如第3節所述,由于將C3和C4視為同一種情況,故以下實驗只分別對C1,C2,C3這3種情況下的重建進行測試.

4.2.1 三維離散型數據重建

首先實驗三維離散型數據的重建效果.圖10(a)是一幅多孔介質(碳酸巖)體數據訓練圖像(80×80×80體素),包含孔隙和骨架2種狀態值,孔隙所占百分比(孔隙度)為0.201;圖10(b)是訓練圖像的截面圖(3個截面分別為X=40,Y=40和Z=40);圖10(c)是訓練圖像的孔隙結構圖.圖11是硬條件數據,分別是孔隙點和骨架點,其中孔隙點比例為0.182,設重建區域為80×80×80體素.由于多孔介質內部只存在孔隙和骨架2種狀態值,因此重建區域內孔隙和骨架在某個體素位置的概率之和為1.圖12(a)(b)分別是孔隙和骨架分布的概率圖,作為重建時的軟條件數據.根據最大似然估計法,設置低維空間的維數d=9[22].情況C1,C2,C3下重建結果分別如圖13～15所示.可見,3種情況下基本都可以重建孔隙和骨架的結構特征.對情況C1～C3下各實現10次重建,設孔隙狀態值為1,骨架狀態值為0,則可得3種情況下重建結果均值和方差,如圖16和圖17所示.3種情況下均值(即孔隙度)見表1所示.求取的方差只有6種取值:0,0.09,0.16,0.21,0.24,0.25,各體素位置的方差值分布如圖18所示.可以看出C3重建結果的方差值最小,說明C3重建結果的波動性最小.C1～C3重建圖像均值和訓練圖像的變差函數曲線如圖19所示,可見訓練圖像與C3重建結果的變差函數最為接近.情況C1～C3下10次重建的時間和占用的最大內存如表2所示,可見C3時使用的時間和內存最大.以上實驗表明,結合使用軟硬數據能夠提高重構質量，但是會占用較多的CPU時間和內存.

Fig. 10 Training image of carbonatite.圖10 碳酸巖訓練圖像

Fig. 11 Hard conditional data of carbonatite.圖11 碳酸巖硬條件數據

Fig. 12 Probability distribution of pore spaces and grains.圖12 孔隙和骨架分布的概率圖

Fig. 13 Reconstructed results of C1.圖13 C1重建結果

Fig. 14 Reconstructed results of C2.圖14 C2重建結果

Fig. 15 Reconstructed results of C3.圖15 C3重建結果

Fig. 16 Average of 10 times of reconstructed carbonatite results under three conditions.圖16 3種情況下的碳酸巖重建結果均值

Fig. 17 Variance of 10 times of reconstructed carbonatite results under three conditions.圖17 3種情況下的碳酸巖重建結果方差

Table 1 Average of the Training Image and Reconstructed Carbonatite Images with C1, C2 and C3

表1 訓練圖像、C1,C2,C3情況下重建的碳酸巖圖像均值

Fig. 18 Variance distribution at each voxel location under three conditions of C1, C2 and C3.圖18 C1,C2,C3情況下各體素位置的方差值分布

Table 2 Maximum Time and Memory for Reconstructing Ten Carbonatite Images with C1, C2 and C3

表2 C1,C2,C3情況下重建10幅碳酸巖圖像的時間和占用的最大內存

ItemMaximumTime∕sMaximumMemory∕MBC16891239C26932243C37353286

Fig. 19 Variogram curves of the training image and average reconstructed images of carbonatite.圖19 碳酸巖訓練圖像和重建圖像均值的變差函數

4.2.2 三維連續型數據重建

進一步實驗三維連續型數據的重建效果.圖20(a)是一幅滲透率體數據的訓練圖像(50×50×20體素,均值為1.634 Darcy),設重建區域也為50×50×20體素.重建所使用的硬數據與軟數據分別如圖20(b)(c)所示.硬條件數據占全部重建區域體素點的1%.情況C1～C3下重建結果分別如圖21～23所示.在C1～C3這3種情況下各實現10次重建,則可得3種情況下重建結果均值和方差,如圖24和圖25所示.各情況下均值見表3所示.各體素位置的方差值分布如圖26所示,可見C3重建結果的方差值最小,說明C3情況下重建結果的波動性最小.C1～C3重建圖像均值和訓練圖像的變差函數曲線如圖27所示,說明C3的重建結果最接近訓練圖像.C1～C3情況下10次重建的時間和占用的最大內存如表4所示,可見3種情況的差異不大,但是C3情況下的重建時間和內存占用量最多.以上實驗再次證明,結合使用軟硬數據有助于提高重構質量,但是也會占用較多的CPU時間和內存.

Fig. 20 Training image, hard and soft conditional data for reconstructing permeability image.圖20 滲透率圖像重建所使用的訓練圖像和軟硬條件數據

Fig. 21 Reconstructed results of C1.圖21 C1重建結果

Fig. 22 Reconstructed results of C2.圖22 C2重建結果

Fig. 23 Reconstructed results of C3.圖23 C3重建結果

Fig. 24 Averages of 10 times of reconstructed permeability results under the conditions of C1, C2 and C3.圖24 C1,C2,C3情況下的滲透率重建結果均值

Fig. 25 Variance of 10 times of reconstructed permeability images under the conditions of C1, C2 and C3.圖25 C1,C2,C3情況下的滲透率重建結果方差

Table 3 Average of the Training Image and the Reconstructed Permeability Images with C1, C2 and C3

表3 訓練圖像、C1,C2,C3情況下重建的滲透率圖像均值

ItemAveragePermeability∕DarcyTrainingImage1．634C11．751C21．723C31．613

Fig. 26 Variance distribution at each voxel location under the conditions of C1, C2 and C3.圖26 C1,C2,C3情況下各體素位置的方差值分布

Table 4 Maximum Memory and Time for Reconstructing Images with C1, C2 and C3

表4 C1,C2,C3情況下重建圖像的時間和最大內存

4.3 與FILTERSIM和DISPAT的比較

Fig. 28 Reconstructed results of FILTERSIM.圖28 FILTERSIM重建結果

如引言所述,FILTERSIM和DISPAT是典型的基于線性降維的MPS方法.現利用FILTERSIM和DISPAT與本方法進行對比.實驗樣本即4.2.1節采用的碳酸巖.對FILTERSIM和DISPAT進行結合軟硬條件數據情況下(即C3情況)的模擬,即采用圖10的訓練圖像、圖11的硬條件數據和圖12的軟條件數據.分別利用FILTERSIM和DISPAT進行10次重建,每個重建結果為80×80×80體素,一些重建結果的孔隙結構如圖28和圖29所示.可見FILTERSIM和DISPAT能較好地再現孔隙結構,只是FILTERSIM生成的孔隙的連通性較好,易于形成較大范圍的連通孔隙空間;而DISPAT生成的孔隙空間連通性稍差,因此大范圍連通的孔隙較少.訓練圖像和3種方法(本文方法、FILTERSIM和DISPAT)重建結果均值的變差函數參見圖30所示.可見本文方法重建結果的變差函數和訓練圖像變差函數最為接近.FILTERSIM和DISPAT生成上述10個結果的時間和內存占用如表5所示.與表2的C3(即同時使用軟硬條件數據)情況相比,可見本方法具有計算速度和內存使用的優勢.

Fig. 29 Reconstructed results of DISPAT.圖29 DISPAT重建結果

Fig. 30 Variograms of the training image and average reconstructed images using three methods.圖30 訓練圖像和3種方法重建圖像的平均變差函數

Table 5 Maximum Time and Memory for Reconstructing Ten Carbonatite Images using FILTERSIM and DISPAT

表5 FILTERSIM和DISPAT重建圖像的時間和最大內存

5 總 結

重建區域內的條件數據對重建結果的準確性影響很大.當已知條件數據較少甚至沒有條件數據時,非確定性插值方法是重建空間數據的有效手段.多點信息統計法是目前非確定性插值方法的重要分支.它從訓練圖像提取模式特征,并將其復制到待模擬區域,以完成對缺失空間數據的重建.但是訓練圖像往往包含大量非線性數據,傳統多點信息統計法處理線性數據效果較好,但是并沒有針對非線性數據提出具體的處理策略.一些基于降維分類思想的多點信息統計法只是采用線性降維方法來處理包括非線性數據在內的所有數據,這可能會造成空間數據重建質量的下降.

本文采用ISOMAP對訓練圖像中的模式數據進行降維處理.由于ISOMAP是一種非線性降維方法,因此適用于降低非線性數據的維度.如果把線性數據視為非線性數據的特殊情況,那么ISOMAP同樣可以用于線性數據降維,這大大拓展了多點信息統計法的應用范圍.基于MPS和ISOMAP提出在無條件數據,只有硬條件數據和結合軟硬條件數據情況下的空間數據重建方法,實驗驗證了方法的有效性.又與傳統采用線性降維的多點信息統計法(即DISPAT和FILTERSIM)進行比較,顯示出采用非線性降維方法處理空間數據的優勢.值得指出的是：雖然本方法采用的隨機模擬會使得重建結果具有一定的不確定性并產生多種可能結果,但是這些結果都是在當前條件數據和先驗模型下對未知區域數據的合理預測,對風險評測研究和分析決策均具有指導意義.

[1]Wang Wencheng. Research and realization of volume rendering in 3D data fields[D]. Beijing: Institute of Software, Chinese Academy of Sciences, 1998 (in Chinese)(王文成. 三維數據場體繪制技術的研究和實現[D]. 北京: 中國科學院軟件研究所, 1998)

[2]Vedadi F, Shirani S. A MAP-based image interpolation method via viterbi decoding of Markov chains of interpolation functions[J]. IEEE Trans on Image Processing, 2014, 23(1): 424-438

[3]Delibasis K K, Kechriniotis A I, Assimakis N D. New closed formula for the univariate hermite interpolating polynomial of total degree and its application in medical image slice interpolation[J]. IEEE Trans on Signal Processing, 2012, 60(12): 6294-6304

[4]Alpay D, Lewkowicz I. Interpolation by polynomials with symmetries[J]. Linear Algebra and Its Applications, 2014, 456: 64-81

[5]Wang Yunsen, Yang Dongsheng, Liu Yinzhong. A real-time look-ahead interpolation algorithm based on Akima curve fitting[J]. International Journal of Machine Tools & Manufacture, 2014, 85: 122-130

[6]Matej S, Kazantsev I G. Fourier-based reconstruction for fully 3-D PET: Optimization of interpolation parameters[J]. IEEE Trans on Medical Imaging, 2006, 25(7): 845-854

[7]Angelo L D, Stefano P D. C1 continuities detection in triangular meshes[J]. Computer-Aided Design, 2010, 42: 828-839

[8]Chen K, Lorenz D A. Image sequence interpolation based on optical flow, segmentation, and optimal control[J]. IEEE Trans on Image Processing, 2012, 21(3): 1020-1030

[9]Pan Rijing. Quadratic B-spline interpolation curves with tangent constraints on data points[J]. Chinese Journal of Computers, 2007, 30(12): 2132-2141 (in Chinese)(潘日晶. 滿足數據點切向約束的二次B樣條插值曲線[J].計算機學報, 2007, 30(12): 2132-2141)

[10]Paiement A, Mirmehdi M, Xie X, et al. Integrated segmentation and interpolation of sparse data[J]. IEEE Trans on Image Processing, 2014, 23(1): 110-125

[11]Jing Minggang, Wu Jitao. Fast image interpolation using directional inverse distance weighting for real-time applications[J]. Optics Communications, 2013, 286: 111-116

[12]Du Yi, Zhang Ting, Lu Detang, et al. An interpolation method using an improved Markov model[J]. Journal of Computer Research and Development, 2012, 49(3): 565-571 (in Chinese)(杜奕, 張挺, 盧德唐, 等. 一種基于改進的Markov模型的插值方法[J].計算機研究與發展, 2012, 49(3): 565-571)

[13]Wang Guozhao, Fang Lincong. On control polygons of quartic pythagorean hodograph curves[J]. Computer Aided Geometric Design, 2009, 26: 1006-1015

[14]Guardiano F, Srivastava M R. Multivariate geostatistics: Beyond bivariate moments[EB/OL]. (1993-03-21)[2014-03-02]. https://pangea.stanford.edu/ researchgroups/ scrf /sites /default /files

[15]Zhang Tuanfeng. Filter-based training pattern classification for spatial pattern simulation[D]. Palo Alto: Stanford University, 2006

[16]Honarkhah M. Stochastic simulation of patterns using distance-based pattern modeling[D]. Palo Alto: Stanford University, 2011

[17]Lei Yingke. The study of manifold learning algorithms and their applications[D]. Hefei: University of Science and Technology of China, 2011 (in Chinese)(雷迎科. 流形學習算法及其應用研究[D]. 合肥: 中國科學技術大學, 2011)

[18]Wang Yong. Manifold learning based classification and clustering approaches with their applications[D]. Changsha: National University of Defense Technology, 2011 (in Chinese)(王勇. 基于流形學習的分類與聚類方法及其應用研究[D]. 長沙: 國防科學技術大學, 2011)

[19]Zhang Ting, Lu Detang, Li Daolun, et al. Research on statistical information reconstruction of images based on multiple-point geostatistics integrating soft data with hard data[J]. Journal of Computer Research and Development, 2010, 47(1): 43-52 (in Chinese)(張挺, 盧德唐, 李道倫, 等. 基于軟硬數據的多點信息統計法在圖像統計信息重構中的應用研究[J]. 計算機研究與發展, 2010, 47(1): 43-52)

[20]Kim K, Lee D. Inductive manifold learning using structured support vector machine[J]. Pattern Recognition, 2014, 47: 470-479

[21]Tenenbaum B J, Silva V, Langford C J. A global geometric framework for nonlinear dimensionality reduction[J]. Science, 2000, 290(22): 2319-2323

[22]Levina E, Bickel P J. Maximum likelihood estimation of intrinsic dimension[C] //Proc of Neural Information Processing Systems (NIPS 2005). Cambridge, MA: MIT Press, 2005: 777-784

[23]Ester M, Kriegel H P, Sander J, et al. A density-based algorithm for discovering clusters in large spatial databases[C] //Proc of Knowledge Discovery and Data Mining(KDD’96). Menlo Park: AAAI Press, 1996: 226-231

[24]Goovaerts P. Geostatistics for Natural Resources Evaluation[M]. New York: Oxford University Press, 1997:45-60

Du Yi, born in 1977. PhD and associate professor in Shanghai Polytechnic University. Member of China Computer Federation. Her main research interest is data mining.

Zhang Ting, born in 1979. PhD and associate professor in Shanghai University of Electric Power. His main research interest is image reconstruction.

Huang Tao, born in 1979. PhD and associate professor in the University of Science and Technology of China. His main research interest is high performance computing.

A Reconstruction Method of Spatial Data Using MPS and ISOMAP

Du Yi1, Zhang Ting2, and Huang Tao3

1(College of Engineering, Shanghai Polytechnic University, Shanghai 201209)2(CollegeofComputerScienceandTechnology,ShanghaiUniversityofElectricPower,Shanghai200090)3(DepartmentofModernMechanics,UniversityofScienceandTechnologyofChina,Hefei230027)

Conditional data influence the reconstructed results greatly in the reconstruction of spatial data. Reconstructed results often show a number of uncertainties when only sparse conditional data are available, so it is suitable to use indefinite interpolation to reconstruct spatial data. As one of the main indefinite interpolation methods, multiple-point statistics (MPS) can extract the intrinsic features of patterns from training images and copy them to the simulated regions. Because the traditional MPS methods using linear dimensionality reduction are not suitable to deal with nonlinear data, isometric mapping (ISOMAP) is combined with MPS to address the above issues. A method using MPS and ISOMAP for the reconstruction of spatial data is proposed for the accurate reconstruction of unknown spatial data by constructing pattern dataset, dimensionality reduction of patterns, classification of patterns and extraction of patterns, which has provided a new idea for dealing with nonlinear spatial data by MPS. The experimental results show that the structural characteristics of reconstructed results using this method are similar to those of training images.

training image; multiple-point statistics (MPS); pattern; isometric mapping (ISOMAP); dimensionality reduction

2015-05-21；

2015-12-22

中國科學院戰略性先導科技專項課題(XDB10030402);中國石油與中國科學院重大戰略合作項目(2015A-4812);國家自然科學基金項目(41672114);上海市自然科學基金項目(16ZR1413200);上海第二工業大學校級重點學科建設計算機科學與技術項目(XXKZD1604) This work was supported by the State Priority Research Program of the Chinese Academy of Sciences (XDB10030402), CNPC-CAS Strategic Cooperation Research Program (2015A-4812), the National Natural Science Foundation of China (41672114), the Natural Science Foundation of Shanghai (16ZR1413200), and the Key Discipline of Computer Science and Technology of Shanghai Polytechnic University (XXKZD1604).

張挺(tingzh@shiep.edu.cn)

TP391.41