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

薄 楊，黃存東

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

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

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

1 Kriging算法

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

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

(1)

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

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

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

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

(2)

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

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

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

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

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

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

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

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

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

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

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

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

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

(3)

(4)

W =

(5)

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

本文系統(tǒng)中，選用最簡單的理論變異函數(shù)模型——線性模型，并取一個整形定值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;

// 通過GetDistance()函數(shù)獲得兩點間距離

}

}

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;//通過指針從點集合中獲取該點z坐標(biāo)

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

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

4.1 工程應(yīng)用

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

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

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

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

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

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

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

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

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

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

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

5 小結(jié)

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

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

參考文獻：

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

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

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

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

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

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