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

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

0 引言

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

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

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

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

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

1.1 普通克里金插值原理

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

(1)

1.2 插值點精度飽和條件

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

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

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

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

1.3 八叉樹鄰域點搜索策略

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

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

(2)

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

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

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

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

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

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

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

(3)

插值計算的具體步驟如下。

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

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

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

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

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

2 實驗驗證

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

圖4 原始樣本點空間分布Figure 4 Spatial distribution of original sample points

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

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

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

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

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

圖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

為驗證本文方法在地質(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個相同的球面三棱錐(八分體)，對每個棱錐進(jìn)行遞歸的橫向和徑向剖分，直到滿足精度為止。SGOG采用“圈層碼(十六進(jìn)制)+八分體碼(八進(jìn)制)+球面位置碼(四進(jìn)制)+徑向位置碼(二進(jìn)制)”的網(wǎng)格編碼模型，其網(wǎng)格體元具有結(jié)構(gòu)簡單、變形適中、拓?fù)潢P(guān)系一致、幾何特征明晰等特點。球體網(wǎng)格機(jī)制為基于體元的地質(zhì)真三維模型的構(gòu)建提供了一種新方法，其原理是將離散點的坐標(biāo)映射到相應(yīng)層次的網(wǎng)格編碼，然后進(jìn)行網(wǎng)格繪制即可。實驗采用SGOG橫向剖分16層(三角形網(wǎng)格邊長約150 m)的網(wǎng)格，在盡可能減少數(shù)據(jù)量的前提下，充分體現(xiàn)地層的完整性及結(jié)構(gòu)特征，地層結(jié)構(gòu)渲染圖如圖6所示，整個地層模型共有網(wǎng)格32 520個。

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

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

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

3 結(jié)論

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