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

房鵬 陳麗鈞

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

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

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

文章編號(hào)：1009-3044（2020）01-0231-02

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

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

1滑動(dòng)鄰域克里金插值

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

1.1鄰域劃分

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

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

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

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

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

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

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

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

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

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

（2）將已知采樣點(diǎn)坐標(biāo)分配在二層網(wǎng)格中，確定每個(gè)采樣點(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ù)克里金方法，對(duì)填充得到的一層網(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實(shí)驗(yàn)結(jié)果與分析

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

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