杜 奕，周雪萍，李遠剛

基于序貫指示模擬的信息預測方法

杜 奕1，周雪萍2，李遠剛3

（1. 上海第二工業(yè)大學計算機與信息學院，上海 201209；2. 新疆油田公司采氣一廠，克拉瑪依 834000；3. 新疆油田公司采油二廠，克拉瑪依 834000）

信息預測在醫(yī)藥、地質、軍事等眾多工程應用領域中發(fā)揮著重要作用。在模擬過程中僅利用少量的已知數(shù)據(jù)預測未知信息非常困難。為了獲得更好的模擬效果，在估計或預測未知區(qū)域的過程中采用內插和外插的方法。作為一種非參數(shù)的估計算法，指示克里格方法雖然能夠在給定周圍條件數(shù)據(jù)的情況下估計任何位置的條件累積分布函數(shù)，但并不能保證結果分布的有效性。序貫指示模擬(Sequential indicator simulation, SISIM)綜合了指示克里格方法和序貫方法在模擬非參數(shù)連續(xù)和離散分布函數(shù)上各自的優(yōu)點。該方法根據(jù)一系列自定義閾值確定非參數(shù)累積分布函數(shù)，模擬過程中利用底部和頂部外插函數(shù)給出分布函數(shù)的形狀。實驗結果表明模擬結果具有與真實數(shù)據(jù)相似的分布特征，證明了該方法的實用性。

插值；指示克立格；序貫指示模擬；外插法；信息預測

0 引言

未知區(qū)域的信息預測在很多科學領域發(fā)揮著非常重要的作用。例如，在醫(yī)藥、軍事、地質、氣象和礦業(yè)等領域，信息預測技術被廣泛應用。在預測未知信息的過程中，插值方法可以利用一些離散點來估計未采樣點的未知屬性，從而根據(jù)一些數(shù)學、物理等規(guī)則建立正確完整的數(shù)學模型。盡管目前已經(jīng)提出了大量的插值方法，但精確的信息模擬和預測仍然難以實現(xiàn)。插值方法主要分為兩類：“確定性”方法和“不確定性”方法。這里的“確定性”是指插值函數(shù)的形式、參數(shù)和結果總是確定的?！按_定性”方法包括距離反比加權法、基于三角網(wǎng)格的方法以及基函數(shù)法等?！安淮_定性”插值方法的“不確定”性，一方面表現(xiàn)在插值形式的隨機性上，另一方面表現(xiàn)在插值參數(shù)的選取和確定需要依賴于概率統(tǒng)計原則[1-2]。不確定性插值方法主要包括地質統(tǒng)計學領域的各種克里格方法和隨機模擬方法?；谧儾詈瘮?shù)的克里格方法和隨機模擬方法都屬于兩點地質統(tǒng)計法[3-4]，它們僅描述空間兩點間的關系，難以重構如曲線形狀這樣的復雜圖形。

不確定性方法一般分為兩類[5]。一類是生成唯一插值結果的插值方法。這種插值方法通常作為低通過濾器用來平滑模擬變量空間中可變的局部細節(jié)。該方法提供一種不確定性的測量手段，如克里格方差。第二類不確定性方法是隨機模擬方法，可以生成多種可能的空間分布結果。這些隨機算法屬于典型的全通過濾器，能夠再現(xiàn)空間可變性的所有范圍。多種隨機結果之間的波動可以對結果中隱含的不確定性特征進行可視的定量測量。

指示克里格是一種非參數(shù)估計算法，該方法在給定條件數(shù)據(jù)的前提下可以估計任意位置上條件累積分布函數(shù)(conditional cumulative distribution function, ccdf)。該方法適用于連續(xù)型和離散型變量的信息預測。指示克里格方法將輸入的數(shù)據(jù)轉換成二進制的指示數(shù)據(jù)，然后對這些轉換過的指示數(shù)據(jù)應用簡單或普通克里格方法。對于連續(xù)型變量，通過克里格方法對在閾值以下的概率確定條件累積概率函數(shù)。對于離散型變量，可以估計屬于某一特定分類的概率分布[4]。

序貫指示模擬(SISIM)將指示克里格方法和序貫方法相結合，實現(xiàn)了非參數(shù)的連續(xù)型與離散型的分布模擬。根據(jù)連續(xù)型序貫指示模擬(SISIM)的理論，本文提出了一種利用外插方法的連續(xù)型序貫指示模擬的信息預測方法，實現(xiàn)結果驗證了該方法的實用性與有效性。

1 概念與方法

1.1序貫模擬

序貫模擬的基本思想是：某一位置u鄰域內的所有已知數(shù)據(jù)(原始數(shù)據(jù)和已模擬的數(shù)據(jù))都可作為條件數(shù)據(jù)，在這一前提下進行模擬?？紤]N個隨機變量Zi的聯(lián)合分布。Zi可以代表：

● 某一區(qū)域內離散在N個網(wǎng)格節(jié)點上的同一屬性；

● 同一點處的N個不同屬性；

● N’個節(jié)點上的K個屬性的聯(lián)合分布，N =KN’。

這N個隨機變量的n個數(shù)據(jù)，其相應的N元的條件累積分布函數(shù)ccdf則可表示成：

這里Zi是第i個屬性值，F(xiàn)N(·)是條件概率分布函數(shù)(ccdf)。為了得到一個來自(1)式的N元樣本，可以由N個相繼的步驟來完成。每一步都是ccdf中的一個抽樣，這樣先前已模擬的數(shù)據(jù)可作為下一個抽樣的條件數(shù)據(jù)，條件數(shù)據(jù)不斷增加，已知信息集由(n)更新為(n+1)，序次考慮所有N個隨機變量，重復上述過程。

1.2克里格

克里格在本質上是一種廣義線性退化算法，它可以在最小方差情況下提供最優(yōu)估計值。該估計值是由硬數(shù)據(jù)的線性組合實現(xiàn)的[5-7]，計算公式表示如下：

式(2)中的z*(u)是信息預測估計值。z(u)表示待預測區(qū)域A的硬數(shù)據(jù)(主要信息)，uA∈是位置向量。z(uα)(α = 1,2, …, n)是位于uα位置的第α個信息采樣數(shù)據(jù)。n是區(qū)域A中的采樣數(shù)據(jù)點數(shù)目。權值λα由克里格方程獲得，uα與u之間的關系式如下：

式(3)中的h是描述uα與u之間距離關系的向量。實際上，式(3)要受限于uα與u之間的關系，例如兩者的變差函數(shù)和協(xié)方差。

1.3指示克里格

指示克里格方法可以應用于離散與連續(xù)型變量。但是，對于連續(xù)型變量，指示克里格方法不能保證概率的單調增長，也無法保證連續(xù)性變量的取值在0和1之間[4]。

連續(xù)型指示變量的定義如下：

式(4)中z(u)為待模擬的變量，zk(k = 1, …, K)為閾值。

對于任意閾值zk，指示克里格可以在已知n個數(shù)據(jù)的條件下估計條件累積概率函數(shù)。公式(5)給出了條件累積概率函數(shù)的定義：

由公式(5)可知，根據(jù)不同的閾值zk, k = 1, …, K，可以計算生成不同的條件累積概率函數(shù)I(u, zk)。

對于指示克里格，實驗變差函數(shù)γI?(h,z)定義為

式(6)中N(h)為所有變差函數(shù)的數(shù)目。本文僅討論連續(xù)型SISIM，不涉及離散型指示克里格方法。

1.4 SISIM

在連續(xù)型SISIM過程中，可以沿著模擬路徑的每個位置，對于每個閾值利用基于指示條件數(shù)據(jù)的指示克里格方法估計ccdf值[4,8]。

SISIM模擬過程中不需要任何高斯假設。在連續(xù)型SISIM過程中，通過估計不超過已限定閾值集的概率計算生成ccdf。給定的閾值越多，對應的條件累積概率函數(shù)越精確。SISIM模擬過程如下[4]：

步驟1 確定閾值：z1, z2,…, zK；

步驟2 定義一條隨機模擬路徑；

步驟3 在模擬路徑的每個位置u上，尋找其周圍的條件數(shù)據(jù)z(u1), z(u2) ,…, z(un)，將每個數(shù)據(jù)z(uα) (α = 1,2,…, n)轉換為指示值表示的向量；

步驟4 對于K個閾值，根據(jù)克里格法估計每個閾值對應的指示隨機向量I(u, zk)，I(u, zk) = Prob{z(u)＜ zk| z(u1), z(u2),…, z(un) }，估算對應于變量z(u)的條件累積概率函數(shù)Fz(u)；

步驟5 從上述條件累積概率函數(shù)中抽取一個值，將其作為位置u上的數(shù)據(jù)值；

步驟6 重復步驟3至步驟5直到模擬完所有位置。

條件累積概率函數(shù)Fz(u)的定義如下：

式(7)中Fzk(u)根據(jù)指示克里格法計算獲得。φlti(z) 和 φuti(z)分別為上，下尾外推(the lower and upper tail extrapolation)。

1.5 上下尾外推

非參數(shù)條件累積概率函數(shù)F(u)根據(jù)一系列閾值z1, z2…, zK進行定義，并以相同的增量1/(K + 1)實現(xiàn)遞增，如F(z1) = 1/(K + 1)， 而F(zK) = K/(K+1)。如果z1和zK不是最小和最大值，則需要對分布函數(shù)F(u)的尾部進行建模。如果z1和 zK分別為最小和最大值則不需要進行任何外推估計。下尾外推函數(shù)提供了最小值zmin和第一個值z1的分布形狀[4,9]。

下尾外推函數(shù)如下定義[4,10]：

(1) Z被限定：F(zmin) = 0。F的下尾外推函數(shù)為

式(8)中的1ω≥用于調節(jié)上述函數(shù)的增加或減少。如果1ω=, 那么zmin和 z1之間的所有值均相等。

(2) Z不被限定：F的下尾外推函數(shù)為

式(9)中的exp為指數(shù)函數(shù)。

上尾外推函數(shù)如下定義：

(1) Z被限定：F(zmin) =1。上尾外推函數(shù)為

式(10)中的01ω≤≤用于調節(jié)上述函數(shù)的增加或減少。如果1ω=, 那么zmax和 zK之間的所有值均相等。

(2) Z不被限定：F的上尾外推函數(shù)為

式(10)中的1ω≥.

2 實驗結果和分析

實驗所用的二維樣本數(shù)據(jù)為實際測得的某區(qū)域的真實海撥值。該數(shù)據(jù)集為包含10 000個數(shù)據(jù)點的點集網(wǎng)格，用于測試評價SISIM的模擬效果。用于條件數(shù)據(jù)的樣本數(shù)據(jù)如圖1所示。背景設為黑色用于突出顯示樣本數(shù)據(jù)。工具條中的不同顏色表示區(qū)域中已知數(shù)據(jù)的不同取值?？梢钥闯鰞H有少量結點為已知數(shù)據(jù)。

圖1 樣本數(shù)據(jù)Fig.1 Sample data

圖1 中的樣本數(shù)據(jù)是從圖2所示的真實海撥數(shù)據(jù)中抽樣獲得。因此可以通過對比模擬結果和圖2的相似性來評估本文算法的有效性。相似度越高，則說明方法越有效。原始圖像數(shù)據(jù)的直方圖如圖3所示。

圖2 原始圖像數(shù)據(jù)Fig.2 Original data

圖3 原始圖像數(shù)據(jù)的直方圖Fig.3 Histogram of original data

下面給出了本文方法的實驗結果。兩個利用SISIM和樣本數(shù)據(jù)生成的隨機模擬結果分別如圖4(a)和(b)所示。結果顯示模擬結果和原始圖像數(shù)據(jù)的結構特征非常相似。

圖4 使用SISIM和樣本數(shù)據(jù)獲得的模擬結果. (a)結果1 (b)結果2Fig.4 Simulated results using SISIM and sample data. (a) result 1; (b) result 2.

上述兩個模擬結果的直方圖分別如圖5(a)和5(b)所示。可以看出模擬結果和原始圖像數(shù)據(jù)直方圖中的數(shù)值分布非常相似，證明本文方法有效。

圖5 序貫指示模擬結果的直方圖 (a)結果1;(b)結果2Fig.5 The histograms of simulated results using SISIM and sample data. (a) result 1; (b) result 2

表1給出了模擬結果與初始圖像數(shù)據(jù)的平均值及方差值，不難看出模擬結果與初始數(shù)據(jù)的平均值和方差值也都非常接近。

表1 初始數(shù)據(jù)和模擬結果的平均值和方差值Tab.1 The average and variance of original data and simulated results

3 結束語

基于序貫指示模擬的信息預測方法可以有效實現(xiàn)未知信息的預測。本文討論了連續(xù)型的序貫指示模擬方法。在不需要任何高斯假設的前提下，SISIM將指示條件數(shù)據(jù)和序貫模擬相結合實現(xiàn)每個待模擬結點上未知數(shù)據(jù)的有效預測。模擬過程中使用了指示克里格和外推法用于計算上下尾部函數(shù)。由實驗結果可以看出，模擬結果與初始數(shù)據(jù)非常相似，從而證明了本文方法的實用性。

[1] ZHANG T , LU D T , LI D L . Porous media reconstruction using a cross-section image and multiple-point geostatistics[C]// Proceedings of ICACC, 2009, Singapore, 2009:24-29.

[2] ZHANG T, LU D T, LI D L . A statistical information reconstruction method of images based on multiple-point geostatistics integrating soft data with hard data[C]//Proceedings of ISCSCT 2008. Shanghai, China, 2008:573-578.

[3] LU D T, ZHANG T, YANG J Q , LI D L , KONG X Y. A reconstruction method of porous media integrating soft data with hard data[J]. Chinese Science Bulletin, 2009, 54(11): 1876-1885.

[4] RENY N, BOUCHER A, WU J B. Applied Geostatistics with SGeMS: A Users’ Guide[M]. Cambridge: Cambridge University Press, 2009:108-134.

[5] ZHANG T, DU Y. The study of cokriging using a Markov model[C]// Proceedings of ICETC 2010. Shanghai, China, 2010:6-11.

[6] ZHANG T, DU Y. A novel interpolation method using a Markov model and COSGSIM[C]// Proceedings of ICCAE 2010. Singapore, 2010: 691-695.

[7] ZHANG T, DU Y. A novel method for information prediction[C]//Proceedings of ICSAP 2010. Bangalore, India, 2010:3-7.

[8] MA X L, JOURNEL A G. An expanded GSLIB cokriging program allowing for two Markov models[J]. Computers & Geosciences, 1999, 25: 627-639.

[9]GOOVAERT SP. Geostatistics for Natural Resources Evaluation[M]. New York: Oxford University Press. 1997: 54-80.

[10] ALMEIDA A, JOURNEL A G. Joint simulation of multiple variables with a Markov-type coregionalization model[J]. Mathematical Geology, 1994, 26(5): 565-588.

An Information Prediction Method Using SISIM

DU Yi1, ZHOU Xue-ping2, LI Yuan-gang3
(1. School of Computer and Information, Shanghai Second Polytechnic University, Shanghai 201209, P. R. China; 2. NO.1 gas production plant of Xinjiang oilfield company, Karamay 834000, P. R. China; 3. NO.2 oil production plant of Xinjiang oilfield company, Karamay 834000, P. R. China）

Information prediction plays an important role in many fields such as medical, geological and military fields. However, it is quite difficult to predict the unknown information only by some sparse known data in the process of simulation. Therefore, some interpolation or extrapolation methods are used to estimate or predict the unknown region for better simulated results. Indicator kriging is a non-parametric estimation algorithm used to estimate the conditional cumulative distribution function at any location given neighboring conditioning data, which cannot guarantee that the resulting distributions are valid. Sequential indicator simulation (SISIM) combines the indicator kriging with the sequential method to simulate non-parametric continuous or categorical distributions. In SISIM, a non-parametric cumulative distribution function is determined from a set of threshold values. The lower tail extrapolation function and the upper tail extrapolation function provide the shapes of the distribution functions during the simulation of SISIM. Experimental results show that the simulated results have the similar characteristics of distributed values with the original one, demonstrating that the proposed method is practical.

interpolation; indicator kriging; SISIM; extrapolation; information prediction

TP391

A

1001-4543(2010)04-0285-06

2010-09-15;

2010-10-28

杜奕（1977－），女，江蘇吳江人，副教授，博士，主要研究方向為數(shù)據(jù)分析與處理，電子郵件：duyi@it.sspu.cn