三維空間屬性體克里金插值方法的研究

房鵬 陳麗鈞

摘要：該文以克里金插值法為基礎(chǔ)，實現(xiàn)了一種滑動鄰域克里金方法，并針對該方法無法應(yīng)用于三維空間屬性體建模的缺陷，提出了一種改進(jìn)算法。從屬性體建模的角度，該文詳細(xì)介紹了滑動鄰域克里金方法基本原理，分析這種方法的優(yōu)勢和缺點(diǎn)，并提出一種改進(jìn)的滑動鄰域克里金方法，實現(xiàn)了其在三維空間中的屬性插值。最后通過體繪制以切片形式展示插值結(jié)果，說明本文提出的改進(jìn)算法可以應(yīng)用到三維空間屬性體建模，并具有較高的效率。

關(guān)鍵詞：克里金插值;滑動鄰域克里金;三維空間屬性體;體繪制

中圖分類號：TP3 文獻(xiàn)標(biāo)識碼：A

文章編號：1009-3044（2020）01-0231-02

克里格法是國際上公認(rèn)的空間插值方法，也是地質(zhì)統(tǒng)計學(xué)的主要方法。它具有線性，無偏和最小方差估計的特點(diǎn)?？死锝鹚惴ㄔ蜑槠胀死锝?，由于該算法實現(xiàn)比較簡單，本文在實現(xiàn)普通克里金的基礎(chǔ)上又實現(xiàn)了一種高效的滑動鄰域克里金方法，該方法由梅鋼等人提出，其原理為：以變程為邊長將研究區(qū)域劃分為正方形網(wǎng)格，構(gòu)造該區(qū)域內(nèi)的克里金方程組并求解得到插值結(jié)果。本文在實現(xiàn)該算法后又研究了大規(guī)模四維屬性數(shù)據(jù)體插值，因采用普通克里金和滑動鄰域克里金方法均無法對其進(jìn)行插值，故對梅鋼等人提出的滑動鄰域克里金算法進(jìn)行了改進(jìn)，實現(xiàn)了一種局部求變程并通過漫水法進(jìn)行鄰域內(nèi)點(diǎn)集填充的滑動鄰域克里金方法。

傳統(tǒng)的克里格法選擇所有的已知點(diǎn)進(jìn)行插值，計算時間較長，尤其在大規(guī)模采集數(shù)據(jù)的情況下，甚至無能為力;并且在實際中，真實的協(xié)方差是估計出來的，而數(shù)據(jù)測量時又不可避免產(chǎn)生誤差，這就有可能導(dǎo)致在大的鄰域上插值產(chǎn)生的均方差比在相對較小的鄰域上插值時的均方差要大。故在計算前對每個待插點(diǎn)分別選擇其鄰域內(nèi)的一部分已知點(diǎn)作為原始估值計算數(shù)據(jù)。因此，如何選取鄰域點(diǎn)成了克里格法亟需解決的問題。孫立雙等提出了基于空間分布權(quán)系數(shù)鄰域選點(diǎn)算法，該算法通過建立插值點(diǎn)鄰域點(diǎn)空間分布權(quán)系數(shù)的函數(shù)，確定參估值計算的鄰域點(diǎn)。劉永社根據(jù)大量插值試驗結(jié)果經(jīng)驗認(rèn)為搜索數(shù)據(jù)采集點(diǎn)時搜索半徑通常取1.3-2倍的變程值。牛文杰提出：利用Delaunay三角網(wǎng)格剖分，建立各個三角形和不同數(shù)據(jù)點(diǎn)的包含關(guān)系，利用三角網(wǎng)格的拓?fù)潢P(guān)系，檢索以三角形某一頂點(diǎn)為共同頂點(diǎn)的所有三角形圍成的凸多邊形內(nèi)包含的數(shù)據(jù)點(diǎn)，直到檢索到的數(shù)據(jù)點(diǎn)滿足要求。杜宇健結(jié)合溫度場計算的實際應(yīng)用詳細(xì)分析了克里格法鄰近點(diǎn)選擇中需要考慮的原則，并采用Delaunay三角劃分搜索和固定距離搜索相結(jié)合的鄰近點(diǎn)搜索策略，提出一種利用變程的Delauuay-固定距離滑動鄰域算法。梅鋼等提出了一種基于變程的滑動鄰域克里金方法，該方法是一種全局求變程，并根據(jù)已知變程進(jìn)行局部克里金方程組計算的插值方法，彌補(bǔ)了Delaunay三角劃分方法等需要對每個待插點(diǎn)必須求一次變差函數(shù)系數(shù)矩陣的逆矩陣的缺點(diǎn)。但是該方法有一個缺陷，即用該方法進(jìn)行大規(guī)模數(shù)據(jù)插值時，進(jìn)行全局變程計算量巨大因而存在無法進(jìn)行插值計算的可能。為解決這一問題，本文對該方法進(jìn)行了改進(jìn)，實現(xiàn)了一種局部求變程，并通過漫水法進(jìn)行鄰域內(nèi)點(diǎn)集填充的克里金插值方法。下面分別對滑動鄰域克里金方法和改進(jìn)的克里金屬性體插值方法進(jìn)行介紹。

1滑動鄰域克里金插值

變差函數(shù)能反映或刻畫區(qū)域化變量的許多重要性質(zhì)，例如可以通過“變程”反映變量相關(guān)性的范圍?；瑒余徲蚩死锝鸩逯到梃b滑動鄰域克里格法的思想，總體思路為：對所有的已知點(diǎn)求變程;確定矩形網(wǎng)格劃分區(qū)域范圍;將矩形區(qū)域劃分為覆蓋上述區(qū)域的正方形網(wǎng)格，網(wǎng)格邊長即為全局變程;計算落在每個正方形單元內(nèi)的已知點(diǎn)和待插點(diǎn);以正方形單元為中心，以其相鄰的正方形單元和自身單元為該單元內(nèi)所有待插點(diǎn)的滑動鄰域。值得注意的是：若鄰域內(nèi)已知點(diǎn)太少，根據(jù)正方形單元在網(wǎng)格中的拓?fù)潢P(guān)系向四周延展，直到滿足條件。本文以三維曲面高程插值為例對滑動鄰域克里金插值方法的鄰域劃分方式和實現(xiàn)步驟進(jìn)行介紹：

1.1鄰域劃分

1.2算法實現(xiàn)步驟

（1）將已知采樣點(diǎn)數(shù)據(jù)加入原始點(diǎn)集合。

（2）對已知所有點(diǎn)進(jìn)行全局范圍求變程a。

（3）根據(jù)求解得到的變程a，用鄰域劃分方法進(jìn)行正方形網(wǎng)格單元劃分。劃分后的網(wǎng)格起點(diǎn)坐標(biāo)為mGridXO和mGridYO。

（4）將已知原始點(diǎn)分配到其各自所屬的網(wǎng)格單元。各點(diǎn)所在的鄰域位置用公式（xO-mGridXO）/a和（yo-mGridYo）la確定。

（5）對所有待插點(diǎn)進(jìn)行遍歷確定其所在鄰域，以全局變程a構(gòu)建局部范圍內(nèi)的克里金方程組，求解方程組得到待插點(diǎn)結(jié)果。

2改進(jìn)的滑動鄰域克里金屬性體插值

本文實現(xiàn)了基于變程的滑動鄰域克里金插值方法，并在三維空間的曲面高程數(shù)據(jù)插值中得到應(yīng)用，但實驗證明該方法無法應(yīng)用到四維屬性數(shù)據(jù)體插值，一是因為屬性數(shù)據(jù)體規(guī)模大，在全局范圍內(nèi)求解變程值會導(dǎo)致計算困難甚至無法得到插值結(jié)果;二是因為當(dāng)采樣點(diǎn)分布不均勻時，通過該方法得到的插值結(jié)果誤差較大。

針對以上兩個問題，本文對以原滑動鄰域克里金插值方法做了改進(jìn)，改進(jìn)方案是在原算法劃分網(wǎng)格的基礎(chǔ)上進(jìn)行更細(xì)致的二層網(wǎng)格劃分，并且限制了鄰域內(nèi)最大點(diǎn)數(shù)和最小點(diǎn)數(shù)。經(jīng)實驗證明該方法既縮短了運(yùn)算時間，又保證了插值效果的合理性。下面同樣以三維曲面高程插值為例對該改進(jìn)算法的步驟進(jìn)行介紹：

（1）獲取已知采樣點(diǎn)集合，用步長step進(jìn)行第一層網(wǎng)格單元的劃分，鄰域劃分方法與原滑動鄰域算法所述一致。之后以step/3為步長對一層大網(wǎng)格單元進(jìn)行二次劃分。

（2）將已知采樣點(diǎn)坐標(biāo)分配在二層網(wǎng)格中，確定每個采樣點(diǎn)在二層網(wǎng)格中所屬的具體單元。

（3）使用漫水法，首先搜索小網(wǎng)格內(nèi)的采樣點(diǎn)，然后將搜索得到的點(diǎn)填充到其所屬的一層大網(wǎng)格，直到大網(wǎng)格內(nèi)達(dá)到規(guī)定的最少點(diǎn)數(shù)要求。

（4）根據(jù)克里金方法，對填充得到的一層網(wǎng)格單元內(nèi)的數(shù)據(jù)點(diǎn)進(jìn)行求解，得到該層網(wǎng)格內(nèi)的局部變程a和理論變差。

（5）遍歷所有待插點(diǎn)，分別確定每一待插點(diǎn)所在的大網(wǎng)格單元，根據(jù)步驟（4）求得的該網(wǎng)格的變程a和變差函數(shù)在當(dāng)前網(wǎng)格單元構(gòu)建克里金方程組，求解方程組得到插值結(jié)果。算法結(jié)束。

3實驗結(jié)果與分析

因四維屬性數(shù)據(jù)體規(guī)模龐大，普通克里金和梅鋼等人提出的滑動鄰域克里金方法無法完成插值運(yùn)算，如圖1所示為改進(jìn)的滑動鄰域克里金插值算法對屬性數(shù)據(jù)體的插值實現(xiàn)結(jié)果。

此次參與運(yùn)算的屬性樣本點(diǎn)共有21365個，普通克里金和梅鋼等人提出的滑動克里金方法經(jīng)長時間計算無法完成屬性數(shù)據(jù)體插值。本文改進(jìn)算法耗時358s，最大誤差為20.37%，最小誤差為0.85%，平均誤差為9.27%。因而可以看出，在數(shù)據(jù)量規(guī)模大的情況下，普通克里金和滑動克里金方法均無法完成四維屬性數(shù)據(jù)體的插值運(yùn)算，而本文改進(jìn)的滑動克里金算法可以在上述情況下正常工作，完成四維屬性數(shù)據(jù)體插值，且插值速度在可接受范圍內(nèi)。