基于八叉樹的修正克里金空間插值算法

王金鑫， 秦子龍， 曹澤寧， 陳藝航， 石 焱

(1.鄭州大學(xué) 地球科學(xué)與技術(shù)學(xué)院，河南 鄭州 450001； 2.鄭州大學(xué) 水利科學(xué)與工程學(xué)院，河南 鄭州 450001)

0 引言

空間插值是一種根據(jù)有限的采樣點(diǎn)數(shù)據(jù)來推求同一區(qū)域其他未知點(diǎn)數(shù)據(jù)的算法。它是基于“地理學(xué)第一定律”的假設(shè)：空間位置上距離越近的點(diǎn)，具有相似特征的可能性越大，反之，就越小[1]。實(shí)際應(yīng)用中，人們所獲取的樣本數(shù)據(jù)往往是有限的。利用插值算法能夠很好地解決由點(diǎn)到面的賦值問題，是地學(xué)科學(xué)研究常用的方法。尤其在地質(zhì)學(xué)研究中，由于地質(zhì)空間錯(cuò)綜復(fù)雜，勘查得到的都是離散的、空間分布不均勻的樣本數(shù)據(jù)，難以表達(dá)真實(shí)的地質(zhì)情況[2]。因此，地質(zhì)現(xiàn)象的空間分布就需要通過插值來實(shí)現(xiàn)[3]。

目前，空間插值主要分為確定性插值與空間統(tǒng)計(jì)插值?？死锝鸩逯底鳛榭臻g統(tǒng)計(jì)學(xué)插值的主要方法[4]，能夠很好地表達(dá)復(fù)雜地質(zhì)空間結(jié)構(gòu)特征[5]，也是目前應(yīng)用最廣泛的一種空間插值方法。孫宗良[6]提出了一種基于變程的滑動(dòng)鄰域克里金插值方法，結(jié)果表明，該方法在三維地質(zhì)建模中效率得到了明顯提升。黃蕾蕾[7]對不同插值方法在地質(zhì)建模中的適應(yīng)性進(jìn)行了分析。王長鵬等[8]基于離散網(wǎng)格曲率對克里金插值所使用的樣本點(diǎn)進(jìn)行選取，極大地提高了地質(zhì)曲面的建模效率。馮波等[9]對反距離加權(quán)插值與自然鄰域法在地質(zhì)建模上的適用性進(jìn)行了比較，結(jié)果表明，反距離加權(quán)插值的精度更高。隨著研究的深入，空間插值方法被應(yīng)用于多個(gè)領(lǐng)域，更多的影響因素也被考慮進(jìn)來。李慧晴等[10]將坡向與地形作為權(quán)重修正指數(shù)對反距離加權(quán)插值進(jìn)行了改進(jìn)，有效減小了降水插值的誤差。曹端廣等[11]在對氣溫空間插值時(shí)，考慮了風(fēng)向與風(fēng)速的影響，提高了氣溫空間插值的準(zhǔn)確性。莫躍爽等[12]使用不同的插值方法對喀斯特地區(qū)降水量進(jìn)行了插值計(jì)算，結(jié)果表明，使用克里金插值得到的效果最好。

克里金插值是一種定義在空間有限域上的算法，在實(shí)際應(yīng)用中常會(huì)涉及計(jì)算區(qū)域選取的問題，即待插值點(diǎn)的鄰域搜索問題。該問題是影響克里金插值效率與精度的重要因素之一[13]?？死锝鸩逯掂徲蛩阉鞯乃惴ㄖ饕泄潭ň嚯x法[14]、固定數(shù)目法[15-16]、Delaunay固定鄰域選擇算法等[17]，但這些算法對于樣本點(diǎn)的空間均勻分布有較強(qiáng)的依賴性，對于分布不均勻的樣本點(diǎn)，插值結(jié)果會(huì)受到一定的影響[4,13]。同時(shí)，在構(gòu)建三維地質(zhì)模型時(shí)，基于傳統(tǒng)的插值方法往往會(huì)存在插值點(diǎn)數(shù)據(jù)冗余、空間分布不均勻等問題。

針對以上問題，本文在普通克里金插值方法的基礎(chǔ)上，利用八叉樹網(wǎng)格剖分改進(jìn)插值點(diǎn)的鄰域搜索算法，優(yōu)化參考點(diǎn)的選取策略，使得插值后的數(shù)據(jù)空間分布更均勻，算法精度與效率進(jìn)一步提高。

1 基于八叉樹的修正克里金插值算法

1.1 普通克里金插值原理

普通克里金插值是地質(zhì)統(tǒng)計(jì)學(xué)中廣泛使用的插值方法，也是地質(zhì)統(tǒng)計(jì)學(xué)的重要組成部分。該方法用來估算未采樣位置的屬性值，是一種最優(yōu)無偏估計(jì)方法：

(1)

1.2 插值點(diǎn)精度飽和條件

空間插值的精度主要與距離相關(guān)[18]，當(dāng)已知點(diǎn)與待插值點(diǎn)距離過遠(yuǎn)時(shí)，則對待插值點(diǎn)的影響很??；同時(shí)，當(dāng)待插值點(diǎn)附近的已知點(diǎn)達(dá)到一定數(shù)量時(shí)，再增加已知點(diǎn)個(gè)數(shù)對插值結(jié)果精度提高的作用很小[19]；且已有的研究表明，插值計(jì)算的鄰域點(diǎn)至少要超過3個(gè)[20]。因此，本文采用交叉驗(yàn)證的方法，即每次從采樣數(shù)據(jù)中選取一些點(diǎn)作為未知點(diǎn)，利用普通克里金插值的方法預(yù)測該點(diǎn)的值。實(shí)驗(yàn)時(shí)，以高程插值為例，已知點(diǎn)的個(gè)數(shù)以5為間隔，由5到90依次遞增。每次從已知點(diǎn)中隨機(jī)選取5個(gè)點(diǎn)作為未知點(diǎn)進(jìn)行驗(yàn)證，最后求出這5個(gè)點(diǎn)的均方根誤差。計(jì)算結(jié)果如圖1(a)所示。

圖1 交叉驗(yàn)證插值結(jié)果圖Figure 1 Interpolation results of cross-validation

由圖1(a)可知，當(dāng)已知點(diǎn)的數(shù)量達(dá)到一定值時(shí)，插值點(diǎn)的變化很小或幾乎不再變化。當(dāng)已知點(diǎn)數(shù)量大于30時(shí)，插值結(jié)果基本上不再發(fā)生變化，即插值點(diǎn)的精度達(dá)到飽和。

為驗(yàn)證上述結(jié)論對于其他屬性值的適用性，本文對某塊土壤中某元素的質(zhì)量分?jǐn)?shù)也進(jìn)行了插值計(jì)算。實(shí)驗(yàn)時(shí)仍以5個(gè)樣本點(diǎn)作為起始樣本點(diǎn)，以5為間隔，依次遞增，在驗(yàn)證結(jié)果時(shí)同樣采用交叉驗(yàn)證的方式。如圖1(b)所示，實(shí)驗(yàn)結(jié)果與上述結(jié)論基本相同，說明這種插值精度規(guī)律具有一般性。

1.3 八叉樹鄰域點(diǎn)搜索策略

根據(jù)克里金插值的基本原理，隨著距離的增大，數(shù)據(jù)的相關(guān)性會(huì)降低。當(dāng)鄰域過大時(shí)會(huì)導(dǎo)致相關(guān)性低的點(diǎn)參與計(jì)算，從而影響插值的效率；而當(dāng)鄰域過小時(shí)則會(huì)導(dǎo)致插值點(diǎn)不足，使得插值結(jié)果精度不能得到保證?；诖耍疚奶岢鲆环N基于八叉樹的空間插值算法，對普通克里金插值中的鄰域搜索策略進(jìn)行改進(jìn)，主要步驟如下。

步驟1構(gòu)建樣本點(diǎn)數(shù)據(jù)的最小外包圍盒。通過確定樣本點(diǎn)數(shù)據(jù)中的xmin、ymin、zmin與xmax、ymax、zmax來獲取最小外包圍盒，即對于任意的樣本點(diǎn)數(shù)據(jù)均滿足下式：

(2)

步驟2對構(gòu)建好的最小外包圍盒進(jìn)行八叉樹剖分，將剖分后的包圍盒進(jìn)行編碼標(biāo)記，從上到下順時(shí)針依次記為0～7。為了對剖分后的子立方體進(jìn)行區(qū)分，第2次剖分得到的子立方體采用“父立方體編碼+子立方體編碼”的方式進(jìn)行編碼，依此類推。以第1次剖分后編碼為“1”的立方體為例，第2次剖分后的8個(gè)小立方體依次記為10～17(圖2)。

圖2 八叉樹剖分示意圖Figure 2 Diagram of octree division

步驟3在對外包圍盒進(jìn)行剖分時(shí)，須有約束條件對剖分次數(shù)限制。根據(jù)上文得到的插值點(diǎn)精度飽和條件，當(dāng)已知點(diǎn)數(shù)量超過15時(shí)，插值點(diǎn)的精度可達(dá)到8 m以下，能夠滿足基本的需求；當(dāng)已知點(diǎn)數(shù)量超過30時(shí)，插值點(diǎn)的精度提高很小。因此，剖分后的立方體內(nèi)樣本點(diǎn)數(shù)量若大于30則繼續(xù)剖分；若大于15小于30則不再剖分；若小于15，則在插值時(shí)采用鄰近(面鄰近中已知點(diǎn)多的)立方體內(nèi)的樣本點(diǎn)數(shù)據(jù)補(bǔ)充已知點(diǎn)數(shù)量至15；若剖分后立方體內(nèi)樣本點(diǎn)數(shù)量為0，則刪除該立方體所占的空間。

步驟4在進(jìn)行鄰域搜索時(shí)，待插值點(diǎn)利用各自歸屬的小包圍盒中的樣本點(diǎn)進(jìn)行插值(需補(bǔ)充的除外)，這樣能夠最大程度保證插值點(diǎn)的精度，同時(shí)，能夠減少計(jì)算量，提高插值效率。

1.4 自適應(yīng)空間插值

在進(jìn)行實(shí)際的插值計(jì)算時(shí)，原始樣本數(shù)據(jù)通常是稀疏且分布不均的。插值的實(shí)質(zhì)是由點(diǎn)的特征推測面的特征，對大多數(shù)的地學(xué)應(yīng)用而言，插值后空間數(shù)據(jù)點(diǎn)的相對均勻分布是較好的情況。一般方法是基于一定的數(shù)學(xué)規(guī)則進(jìn)行加密，但由于樣本點(diǎn)的非均勻性會(huì)導(dǎo)致最終空間點(diǎn)的分布不均。本文從建模數(shù)據(jù)出發(fā)，通過定義“點(diǎn)密度”來約束插值點(diǎn)的密度與分布，從而達(dá)到自適應(yīng)加密的效果。

設(shè)總的樣本點(diǎn)數(shù)量為N，剖分后的每個(gè)小立方體內(nèi)樣本點(diǎn)與插值點(diǎn)的數(shù)量之和為Ni(Ni的初始值為各個(gè)小立方體內(nèi)樣本點(diǎn)的數(shù)量，即Ni0)，剖分后總的立方體個(gè)數(shù)為n(除去不包含任何數(shù)據(jù)的空立方體)，研究區(qū)域點(diǎn)的總數(shù)量為T(包括樣本點(diǎn)與插值點(diǎn)，T的初始值為樣本點(diǎn)的數(shù)量，即T0)，則初始時(shí)的點(diǎn)密度di0=Ni0/T0。各個(gè)小立方體加密點(diǎn)數(shù)量ti與點(diǎn)密度di的關(guān)系滿足：

(3)

插值計(jì)算的具體步驟如下。

步驟1構(gòu)建樣本點(diǎn)的最小外包圍盒，分別在x、y和z方向上根據(jù)實(shí)際需求進(jìn)行等間距剖分，將其剖分為i×j×k大小的格網(wǎng)，將位于研究區(qū)域范圍內(nèi)的格點(diǎn)作為待估點(diǎn)進(jìn)行插值計(jì)算[30-31]。

步驟2根據(jù)式di=Ni/T計(jì)算第1次插值后各個(gè)小立方體內(nèi)的點(diǎn)密度，判斷di是否符合要求(當(dāng)di=1/n時(shí)，認(rèn)為達(dá)到要求)。

步驟3若未達(dá)到要求，則繼續(xù)進(jìn)行加密計(jì)算。對于未達(dá)到要求的小立方體，對其內(nèi)部點(diǎn)繼續(xù)加密。在繼續(xù)加密時(shí)，根據(jù)式(3)與步驟1中的方法來確定加密點(diǎn)的數(shù)量與位置。

步驟5重復(fù)步驟3、4，直到滿足要求為止(各個(gè)小立方體內(nèi)的點(diǎn)密度di=1/n時(shí))。

圖3 基于八叉樹的空間插值流程圖Figure 3 Flow chart of spatial interpolation based on octree

2 實(shí)驗(yàn)驗(yàn)證

實(shí)驗(yàn)數(shù)據(jù)采用河南省鄭州市航空港經(jīng)濟(jì)區(qū)某地層的實(shí)測鉆孔數(shù)據(jù)，*.obj格式，Beijing54坐標(biāo)系。鉆孔數(shù)據(jù)在三維空間中的離散分布見圖4，共計(jì)69個(gè)樣本點(diǎn)，經(jīng)度方向的極差為14 889.82 m，緯度方向的極差為19 805.86 m，高程的極差為72.37 m。

圖4 原始樣本點(diǎn)空間分布Figure 4 Spatial distribution of original sample points

對上述數(shù)據(jù)，分別應(yīng)用基于固定距離、固定數(shù)目鄰域搜索策略的普通克里金插值[21-22]，反距離加權(quán)插值，本文方法進(jìn)行加密計(jì)算。本次實(shí)驗(yàn)以Visual Studio 2017為開發(fā)平臺(tái)，以O(shè)SG(open scene graph)作為圖形引擎，調(diào)用Eigen庫參與計(jì)算，以標(biāo)準(zhǔn)C++作為開發(fā)語言，構(gòu)建實(shí)驗(yàn)系統(tǒng)并進(jìn)行計(jì)算。實(shí)驗(yàn)計(jì)算機(jī)配置：DELL Tower 5810，8核E5-1620V33.5 GHz 處理器，16 G內(nèi)存，256 G SSD，Windows10 專業(yè)版。

加密計(jì)算時(shí)，待估點(diǎn)的選擇采用常規(guī)的方法，用70 m×70 m×25 m的網(wǎng)格覆蓋整個(gè)研究區(qū)域，將落在研究區(qū)域內(nèi)的格點(diǎn)作為待插值點(diǎn)。針對本次實(shí)驗(yàn)數(shù)據(jù)，固定距離法中的距離設(shè)置為樣本點(diǎn)間距離最大值的1/4(800 m)和1/8(400 m)，這樣設(shè)置可以保證待插值點(diǎn)鄰域內(nèi)有較為合適數(shù)量的樣本點(diǎn)進(jìn)行計(jì)算，保證計(jì)算結(jié)果的精度；在設(shè)置固定數(shù)目法中的鄰域點(diǎn)數(shù)時(shí)，依據(jù)本文1.2節(jié)計(jì)算結(jié)果將其設(shè)為15、30，以保證插值結(jié)果的精度。然后采用交叉驗(yàn)證的方法(依次將每個(gè)樣本點(diǎn)作為未知點(diǎn))進(jìn)行檢驗(yàn)計(jì)算，以誤差絕對值的最大值、最小值與均方根誤差作為評定指標(biāo)，并對不同方法耗時(shí)進(jìn)行比較，結(jié)果如表1所示。

表1 不同插值方法的比較Table 1 Comparison of different interpolation methods

本次實(shí)驗(yàn)參與計(jì)算的插值點(diǎn)數(shù)共106 126個(gè)。由表1可知，本文方法在絕對誤差的最大值、最小值、均方根誤差以及耗時(shí)方面均優(yōu)于使用固定數(shù)目法(15個(gè)樣本點(diǎn))、固定距離法(400 m)、固定距離法(800 m)鄰域搜索策略的普通克里金插值和反距離加權(quán)插值。與使用固定數(shù)目法(30個(gè)樣本點(diǎn))鄰域搜索策略的普通克里金插值相比，本文方法在結(jié)果精度上較低，但在效率上要遠(yuǎn)優(yōu)于該方法(耗時(shí)僅為固定數(shù)目法(30個(gè)樣本點(diǎn))的15%)。

將上述采用固定數(shù)目法(30個(gè)樣本點(diǎn))和固定距離法(800 m)的普通克里金插值和本文方法的結(jié)果可視化，如圖5所示。為更直觀地比較不同方法的結(jié)果，分別計(jì)算不同方法得到的結(jié)果中相鄰點(diǎn)的最大、最小距離，結(jié)果如表2所示。由圖5與表2可知，使用固定距離法(400 m)和固定數(shù)目法(30個(gè)樣本點(diǎn))鄰域搜索策略的普通克里金插值得到的結(jié)果，其插值點(diǎn)在空間分布上存在過疏或過密的情況，使用本文的方法得到的結(jié)果，能夠自適應(yīng)插值加密樣本點(diǎn)，使得插值點(diǎn)在空間上分布更為均勻。

圖5 不同插值方法結(jié)果可視化比較Figure 5 Results visualization based on different interpolation methods

表2 不同方法得到的插值結(jié)果的空間分布比較Table 2 Interpolation results obtained by different methods compared in spatial distribution

為驗(yàn)證本文方法在地質(zhì)真三維建模中的應(yīng)用，采用球體離散網(wǎng)格中的球體測地線八叉樹網(wǎng)格[23](sphere geodesic octree grid，SGOG)構(gòu)建上述地層的真三維模型。SGOG網(wǎng)格利用大圓弧和半徑中分的剖分規(guī)則，以 0°～180°首子午圈、東西經(jīng) 90°子午圈、赤道3條大圓弧為界線，將地球體剖分成8個(gè)相同的球面三棱錐(八分體)，對每個(gè)棱錐進(jìn)行遞歸的橫向和徑向剖分，直到滿足精度為止。SGOG采用“圈層碼(十六進(jìn)制)+八分體碼(八進(jìn)制)+球面位置碼(四進(jìn)制)+徑向位置碼(二進(jìn)制)”的網(wǎng)格編碼模型，其網(wǎng)格體元具有結(jié)構(gòu)簡單、變形適中、拓?fù)潢P(guān)系一致、幾何特征明晰等特點(diǎn)。球體網(wǎng)格機(jī)制為基于體元的地質(zhì)真三維模型的構(gòu)建提供了一種新方法，其原理是將離散點(diǎn)的坐標(biāo)映射到相應(yīng)層次的網(wǎng)格編碼，然后進(jìn)行網(wǎng)格繪制即可。實(shí)驗(yàn)采用SGOG橫向剖分16層(三角形網(wǎng)格邊長約150 m)的網(wǎng)格，在盡可能減少數(shù)據(jù)量的前提下，充分體現(xiàn)地層的完整性及結(jié)構(gòu)特征，地層結(jié)構(gòu)渲染圖如圖6所示，整個(gè)地層模型共有網(wǎng)格32 520個(gè)。

圖6 地層結(jié)構(gòu)渲染圖Figure 6 Stratum structure rendering

分別采用基于固定數(shù)目法(30個(gè)樣本點(diǎn))、固定距離法(800 m)的普通克里金插值與本文方法進(jìn)行插值計(jì)算。實(shí)驗(yàn)發(fā)現(xiàn)：采用基于固定數(shù)目法(30個(gè)樣本點(diǎn))的普通克里金插值法構(gòu)建上述層次的地層網(wǎng)格需要120 104個(gè)插值點(diǎn)，基于固定距離法(800 m)的普通克里金插值需要122 004個(gè)插值點(diǎn)，本文方法只需88 132個(gè)插值點(diǎn)。本文方法減少將近1/3的冗余點(diǎn)。為更直觀地展示不同方法的建模效果，以88 132個(gè)插值點(diǎn)為限，分別使用基于固定數(shù)目法(30個(gè)樣本點(diǎn))、固定距離法(800 m)的普通克里金插值得到的數(shù)據(jù)進(jìn)行建模，結(jié)果如圖7所示?？梢钥闯?，使用常規(guī)搜索策略建立的模型，在相同剖分層次下，會(huì)出現(xiàn)一些漏洞(圖7)，而本文方法則不然(圖6)。

圖7 固定數(shù)目法與固定距離法建模結(jié)果Figure 7 Modeling results of fixed number method and fixed distance method

3 結(jié)論

本文提出了一種基于八叉樹的修正克里金空間插值算法。依據(jù)普通克里金插值原理，設(shè)計(jì)相應(yīng)的精度驗(yàn)證實(shí)驗(yàn)，得出其精度飽和的閾值；通過八叉樹空間剖分，優(yōu)化其鄰域搜索策略；顧及加密點(diǎn)的空間分布均衡，提出基于“點(diǎn)密度”的自適應(yīng)加密算法，并與傳統(tǒng)的插值算法進(jìn)行對比實(shí)驗(yàn)。結(jié)果表明：本文方法在插值精度與效率上均優(yōu)于大多數(shù)傳統(tǒng)插值算法，同時(shí)保證了加密點(diǎn)在空間分布的相對均勻性，有效減少了數(shù)據(jù)冗余。本文方法能夠應(yīng)用于多種基于離散點(diǎn)的空間插值問題，可將該方法用于對其他屬性插值的場景中(如利用空間離散點(diǎn)數(shù)據(jù)插值水體中微量元素含量、水質(zhì)重金屬污染擴(kuò)散程度等)。同時(shí)，對于不同的數(shù)據(jù)，可根據(jù)實(shí)際情況將本文所優(yōu)化的鄰域搜索策略應(yīng)用于其他插值算法當(dāng)中。本研究提出的方法僅針對插值算法中的鄰域搜索進(jìn)行了改進(jìn)，對于插值過程中其他步驟的改進(jìn)、在插值計(jì)算時(shí)考慮更多的影響因素以及減小運(yùn)算復(fù)雜度是下一步研究的方向。本文方法為插值模型中的鄰域搜索和地球系統(tǒng)空間的建模與表達(dá)提供了一種有效的技術(shù)手段。