基于Kriging插值算法的三維地形構(gòu)造研究

薄 楊，黃存東

(安徽國(guó)防科技職業(yè)學(xué)院信息技術(shù)學(xué)院，安徽 六安 237011)

Kriging算法是現(xiàn)今三維可視化建模中一種進(jìn)行數(shù)據(jù)插值的優(yōu)秀算法.該算法是一種基于線性回歸分析方法的算法，并且在線性回歸分析方法的基礎(chǔ)上進(jìn)行了改進(jìn).該算法由線性回歸以及非參數(shù)兩個(gè)部分組成，線性回歸的不同選擇對(duì)所模擬模型的性質(zhì)不會(huì)產(chǎn)生很大影響，而隨機(jī)過(guò)程的實(shí)現(xiàn)通過(guò)非參數(shù)部分得以實(shí)現(xiàn)，隨機(jī)過(guò)程與高斯分布一致，同時(shí)通過(guò)極大似然估計(jì)可以確定協(xié)方差矩陣的系數(shù)[1,2].

Kriging在進(jìn)行某一未知點(diǎn)的計(jì)算時(shí)，是基于周圍點(diǎn)的一系列已知信息來(lái)完成的，在一定的區(qū)域內(nèi)通過(guò)線性組合中的不同信息確定不同的權(quán)值來(lái)估計(jì)未知點(diǎn)的信息，最小化估計(jì)值誤差的方差用來(lái)確定不同的權(quán)值.由此可見，Kriging算法是一種對(duì)線性無(wú)偏估計(jì)的最優(yōu)化表達(dá).

1 Kriging算法

Kriging模型由兩個(gè)部分構(gòu)成：多項(xiàng)式部分和隨機(jī)分布部分.模型的具體形式為：

y(x) =F(β，x) +z(x) =fT(x)β+z(x)

(1)

其中β為回歸系數(shù)，f(x)為變量x的多項(xiàng)式函數(shù).f(x)提供模擬空間中的全局近似，局部近似由z(x)進(jìn)行模擬.

一般情況下，出于簡(jiǎn)化模型的考慮， f(x)通常取常量1，而回歸量并不能對(duì)模擬的精度起到?jīng)Q定性作用.z(x)是一個(gè)不獨(dú)立存在的隨機(jī)過(guò)程，服從正態(tài)分布N(0，σ2)，且協(xié)方差非零.

z(x)的協(xié)方差矩陣為：

Cov[Z(xi),Z(xj)]=σ2R[R(xi,xj)]

(2)

任意兩個(gè)點(diǎn)xi和xj在空間內(nèi)的相關(guān)性由R(xi,xj)來(lái)確定.

2 Kriging估值算法的優(yōu)勢(shì)

Kriging 算法先確定一個(gè)可能影響插值點(diǎn)數(shù)據(jù)計(jì)算結(jié)果的區(qū)域范圍，在這個(gè)范圍內(nèi)選擇已知點(diǎn)信息進(jìn)行數(shù)據(jù)插值計(jì)算，該算法考慮了空間內(nèi)已知點(diǎn)對(duì)未知點(diǎn)的影響.Kriging算法考慮了已知點(diǎn)與未知點(diǎn)在空間分布的特性，并結(jié)合多種點(diǎn)屬性信息后，為達(dá)到最小估計(jì)方差、無(wú)偏差和線性的估計(jì)，對(duì)不同屬性賦與不同的系數(shù)，在此基礎(chǔ)上進(jìn)行加權(quán)平均[3].

與樣條插值法等傳統(tǒng)的估值方向相比，Kriging估值算法有明顯的優(yōu)勢(shì)：

(1)考慮了網(wǎng)格化空間中被描述點(diǎn)的空間相關(guān)度，估值計(jì)算更接近于實(shí)際，也更加準(zhǔn)確；

(2)給出了Kriging方差，數(shù)值估算的可靠度清晰.

從這個(gè)意義上說(shuō)，只有Kriging法才是一種真正的插值方法，Kriging估值技術(shù)是一種國(guó)際上公認(rèn)的空間估值技術(shù)[4].

3 基于Kriging插值的地形構(gòu)造方法設(shè)計(jì)

有限的已知點(diǎn)信息是進(jìn)行三維地形模擬的基礎(chǔ)，本文三維地形構(gòu)造系統(tǒng)數(shù)據(jù)以工程區(qū)域的等高線圖例為依據(jù)，進(jìn)行人工讀點(diǎn).在人工讀點(diǎn)過(guò)程中可以靈活地對(duì)數(shù)據(jù)點(diǎn)進(jìn)行修改和修正，以保障最終的成圖效果.

在三維地形構(gòu)造的過(guò)程中，系統(tǒng)首先依據(jù)已知的數(shù)據(jù)信息對(duì)離散數(shù)據(jù)點(diǎn)進(jìn)行插值，即根據(jù)已知點(diǎn)對(duì)投影面中的Z坐標(biāo)進(jìn)行插值，該系統(tǒng)采取普通Kriging法進(jìn)行插值；其次，使用一定的規(guī)則對(duì)離散點(diǎn)進(jìn)行擬合，來(lái)形成一個(gè)連續(xù)的曲面.最終對(duì)曲面進(jìn)行色彩渲染.系統(tǒng)結(jié)構(gòu)設(shè)計(jì)如圖1所示.

圖1 三維地形構(gòu)造系統(tǒng)結(jié)構(gòu)設(shè)計(jì)

普通克里格法作為Kriging算法中最廣泛的一種估值方法，是本文模型構(gòu)造過(guò)程中數(shù)值估算所采用的算法.假設(shè)X0為待估計(jì)點(diǎn)，首先計(jì)算第n點(diǎn)和第m點(diǎn)的變異量，變異量cov(xn,xm)與兩點(diǎn)距離相關(guān)，所以cov(xn,xm)= cov(xm,xn).

普通克立格法中點(diǎn)X0處的估計(jì)值Z(X0)=Y′W-1B，其中

(3)

(4)

W =

(5)

可以看出，由于cov(xn,xm)= cov(xm,xn)，因此可知W矩陣為對(duì)稱矩陣.

本文系統(tǒng)中，選用最簡(jiǎn)單的理論變異函數(shù)模型——線性模型，并取一個(gè)整形定值N作為線性斜率，則相對(duì)于兩點(diǎn)間距離，變異函數(shù)Z(h)=N*h.

以下分別是創(chuàng)建W、B、Y矩陣程序代碼及注釋：

W: for(int j=0; j

{

for(int i=0; i< Size; i++){

if(i == Size -1 || j == Size e-1)

{

W[i][j] = 1;

if(i == Size -1 && j == Size -1)

W[i][j] = 0;

continue;

}

W[i][j] = ::GetDistance(first, i, j) * N;

// 通過(guò)GetDistance()函數(shù)獲得兩點(diǎn)間距離

}

}

B: for(int i=0; i< Size-1; i++) {

double distance = ::GetDistance(x, y, first, i);

B[i] = distance * N;

}

B[Size -1] = 1;

Y: for(i=0; i< Size -1; i++)

Y[i] = (*(first+i)).z;//通過(guò)指針從點(diǎn)集合中獲取該點(diǎn)z坐標(biāo)

通過(guò)以上矩陣可計(jì)算，Z(X0)=Y’W-1B,在系統(tǒng)中對(duì)逆矩陣的求解采取了矩陣LU分解法以加快計(jì)算速度，該部分的輸出結(jié)果Z(X)為繪制曲面模型提供數(shù)據(jù)支持[5,6].

4 三維地形構(gòu)造

4.1 工程應(yīng)用

本文基于一個(gè)地形相對(duì)復(fù)雜的工程區(qū)域進(jìn)行三維面模型構(gòu)造，在地形構(gòu)模中采取等高線讀點(diǎn)方式，此方式可以精確地描述地形的細(xì)節(jié)特點(diǎn)，相比較網(wǎng)格剖分的讀點(diǎn)方式有著自身的優(yōu)勢(shì).在讀點(diǎn)前首先剔除一些對(duì)模型構(gòu)造無(wú)用的等高線，這樣在初始階段就人為地簡(jiǎn)化了模型構(gòu)造的難度，減少了軟件運(yùn)算的時(shí)間.

該工程中的三維地形構(gòu)造所取等高線如圖2所示，原始地形點(diǎn)為圖中圓點(diǎn)所示，可見在平緩區(qū)域所讀點(diǎn)少，而在地形復(fù)雜、變化多樣區(qū)域，所讀點(diǎn)數(shù)明顯增多.此操作的目的是增加復(fù)雜地形區(qū)域的點(diǎn)相關(guān)性.

圖2 原始點(diǎn)數(shù)據(jù)采集示意圖

圖3 原始數(shù)據(jù)

點(diǎn)坐標(biāo)以文本形式存儲(chǔ)于文件當(dāng)中，形式如圖3，每行三個(gè)數(shù)據(jù)分別為X、Y和高程Z.插值后構(gòu)建模型如圖4所示.

圖4 基于Kriging插值的三維地形構(gòu)造效果

由此可見，基于Kriging 插值的三維面模型構(gòu)造細(xì)致地展示了工作區(qū)域的三維地形地貌.

4.2 傳統(tǒng)三維地形構(gòu)造比較

在傳統(tǒng)地形展示中，人們更多地在CAD的基礎(chǔ)上基于傳統(tǒng)插值方法進(jìn)行二次開發(fā)，實(shí)現(xiàn)對(duì)三維地形的構(gòu)造.下圖5是實(shí)際工作中依托CAD平臺(tái)開發(fā)的三維地形軟件成圖效果，該軟件使用LISP語(yǔ)言基于樣條插值算法和網(wǎng)格剖分進(jìn)行CAD二次開發(fā).

此時(shí)，基于CAD的二次開發(fā)構(gòu)模只實(shí)現(xiàn)了地形輪廓的簡(jiǎn)單構(gòu)造，且成圖速度較慢.而以同樣圖例，使用Kriging 插值算法進(jìn)行模型構(gòu)造如圖6所示，構(gòu)造速度快，構(gòu)造效果較好.

圖6 Kriging 插值算法地表構(gòu)模

5 小結(jié)

本文分析了Kriging算法的計(jì)算過(guò)程，以普通Kriging插值算法為例，基于OPENGL技術(shù)，在C++平臺(tái)上設(shè)計(jì)并實(shí)現(xiàn)了三維地表模型構(gòu)造軟件，針對(duì)某一地區(qū)的復(fù)雜地形進(jìn)行了三維地形的建模和可視化顯示，并取得了較好的效果.

Kriging算法考慮了空間中空間屬性的分布變異特性，該算法首先確定對(duì)待估值點(diǎn)會(huì)有影響的一個(gè)空間區(qū)域，采集這個(gè)區(qū)域內(nèi)的已知點(diǎn)來(lái)進(jìn)行未知點(diǎn)的數(shù)值估計(jì).估值過(guò)程中考慮了空間分布中多種屬性之間的相互關(guān)系，是一種最優(yōu)的插值算法.

參考文獻(xiàn)：

[1]伊堯國(guó),劉慧平,齊建超,等.基于改進(jìn)Kriging插值模型的城市地面沉降變形趨勢(shì)面模擬[J].大地測(cè)量與地球動(dòng)力學(xué),2017,(9):898-902,927.

[2]張娟.有約束線的普通kriging插值算法的研究[J].黑龍江科技信息,2017,(12):14.

[3]李靜思,潘潤(rùn)秋,范馥麟.基于Kriging模型的地面氣溫空間插值研究[J].西南師范大學(xué)學(xué)報(bào)(自然科學(xué)版),2016,(5):21-27.

[4]李莎,舒紅,董林.基于時(shí)空變異函數(shù)的Kriging插值及實(shí)現(xiàn)[J].計(jì)算機(jī)工程與應(yīng)用,2011,(23):25-26,38.

[5]雷能忠,王心源,蔣錦剛,等.不同地形與取樣數(shù)的Kriging插值精度對(duì)比研究——以舒城縣土壤全氮空間分異為例[J].水文地質(zhì)工程地質(zhì),2008,(5):86-91.

[6]曾懷恩,黃聲享.基于Kriging方法的空間數(shù)據(jù)插值研究[J].測(cè)繪工程,2007,(5):5-8,13.