一種自適應(yīng)選取參考點的Kriging 插值方法研究

鄧顯石

（湖南環(huán)境生物職業(yè)技術(shù)學(xué)院園林學(xué)院，湖南 衡陽421005）

空間插值是利用已知的樣本點來估算其它未知點數(shù)值的過程，常用的空間插值方法主要有IDW、Kriging、趨勢面法、自然鄰域法等。其中Kriging 以變異函數(shù)理論和結(jié)構(gòu)分析為基礎(chǔ)，理論上能夠獲得更高的插值精度，在地質(zhì)、水文、土壤等領(lǐng)域獲得廣泛應(yīng)用[1]。

Kriging 插值算法精度的提高主要受到變異函數(shù)、樣點的分布、參考點的數(shù)量等因素的影響[2]。作為直接影響Kriging 插值誤差的變異函數(shù)，關(guān)于其研究較多。如李莎與徐愛萍等在2012 年提出采用反復(fù)實驗人工擬合變異函數(shù)[3,4]，2016 年徐武平提出一種變異函數(shù)的自動擬合方法[5]。作為影響程度次于變異函數(shù)的參考點的分布及數(shù)量選取，目前卻基本是通過實驗對比方式來進(jìn)行調(diào)優(yōu)。如張仁鐸在2005 年認(rèn)為選擇遠(yuǎn)距離的數(shù)據(jù)點可能會更容易違反二階平穩(wěn)的假設(shè)，一般選擇鄰近的10 個點左右即可[6]；李君在2010 年通過實驗發(fā)現(xiàn)參考點數(shù)在40-60 個點最為合適[7]；olea 認(rèn)為鄰近點的選擇一般應(yīng)介于3-25 點之間[8]。關(guān)于參考樣本點分布均勻能夠提高插值精度已得到多數(shù)學(xué)者的研究證實，如孫益權(quán)通過實驗發(fā)現(xiàn)樣本點分布均勻越能夠提高kriging 的插值精度[9]，但目前顧及樣本點分布不均勻性的自適應(yīng)選取參考樣點的Kriging 插值方法研究較少。

針對此種現(xiàn)狀，本文提出一種自適應(yīng)選取參考樣點的Kriging 方法--OOK。該方法將待插值點的“一階鄰近點”（見后文定義1）作為最終插值的參考點，能夠有效降低

圖7 選擇特征隨吸合次數(shù)增加的變化曲線樣本點采樣分布不均勻?qū)ψ罱K插值結(jié)果造成的影響，同時不需要用戶指定參考點的數(shù)目。

1 普通Kriging 模型

在Kriging 插值中，待插值點的屬性由周邊已知樣本點變量值的線性組合求得[10]，其估計值為：

式中，Y(xi,xj) 為鄰近點xi，xj之間的半方差Y(xi,x0)為鄰近點xi與待插值點x0之間的半方差，可從擬合的半變異函數(shù)模型（球形、高斯、線性模型）中計算得到，φ 為拉格朗日算子。

2 OOK 插值方法

2.1 定義1：將待插值點與研究區(qū)域內(nèi)已知的樣本點建立不規(guī)則三角網(wǎng)，將與待插值點形成鄰接拓?fù)潢P(guān)系的樣本點稱為“一階鄰近點”，如圖（1）所示。由圖可以發(fā)現(xiàn)，采用“一階鄰近點”作為待插值點的參考點，能夠在一定程度上降低樣本點在一側(cè)發(fā)生聚集對參考點選擇造成的影響。

圖1 一階鄰近點

2.2 技術(shù)路線

OOK 插值算法主要包括數(shù)據(jù)清洗，參考點的自適應(yīng)選取，半變異函計算，待求點插值4 個過程，如圖2 所示。

2.2.1 數(shù)據(jù)清洗：現(xiàn)實世界復(fù)雜多變，無論采用何種方式獲取的空間數(shù)據(jù)，均存在一些不可避免的問題[12]，如重復(fù)的空間數(shù)據(jù)點。由于需要將空間點數(shù)據(jù)生成不規(guī)則三角網(wǎng)，因此首先應(yīng)清洗空間數(shù)據(jù)中存在的重復(fù)點。

2.2.2 鄰近點的自適應(yīng)選?。孩偈紫壤脛⒂篮陀?008 年提出的“一種快速生成平面Delaunay 三角網(wǎng)的橫向擴展法”[13]，將n 個已知樣本點建立Delaunay 三角網(wǎng)。對于每個待插值點，不將其與全部樣點重新建立不規(guī)則三角網(wǎng)，而是采取逐步加點法構(gòu)建新的Delaunay 三角網(wǎng)，如此將提高Delaunay 三角網(wǎng)的構(gòu)網(wǎng)效率；②將待插值點q 的“一階鄰近點”作為初始參考點集合；③計算“一階鄰近點”與待插值點之間的距離edgei，當(dāng)edgei＞R(變程[2])，則將此點從初始參考點集合中去除。

圖2 OOK 插值流程

2.2.3 半變異函數(shù)值的計算：本文中半方差函數(shù)模型的擬合采用GS+軟件，其提供的半變異函數(shù)有球狀、高斯等模型，在擬合各類模型時都需設(shè)置步長與“最大滯后距”兩個關(guān)鍵參數(shù)。根據(jù)相關(guān)研究指出，“最大滯后距”一般取研究區(qū)域樣點間最大長度的1/2，步長值的設(shè)置應(yīng)使得各組步長內(nèi)存在不低于30 個樣本點對[7]。

2.2.4 計算待求點的屬性值：根據(jù)步驟2.2.3 求得參考點與待插值點以及參考點之間的半變異值，代入上式（3），計算得到待插值點對應(yīng)的權(quán)重系數(shù)，從而代入公式（1）求得待插值點的屬性值，并對其進(jìn)行精度評估。

3 實驗對比分析

本文采用“Spatial Interpolation Comparison 97 (SIC97)”數(shù)據(jù)集，收集了瑞士1986 年5 月8 日467 個降水站點觀測值，隨機從467 個站點中抽取100 個站點作為訓(xùn)練集，另367 個站點的數(shù)據(jù)作為檢驗數(shù)據(jù)集，467 個站點的分布如圖3 所示，數(shù)據(jù)來源于https://www.researchgate.net/profile/Gregoire_Dubois/publication/281292076_Spatial_Interpolation_Comparison_97_ (SIC97) _dataset）。

圖3 氣象站點分布圖

為檢驗OOK 方法的可行性，將其與OK、AIDW[14]進(jìn)行比較。對比方法為交叉驗證法，比較OK 在不同的選點條件下與OOK、AIDW 的MAE、RMSE、MAX、MIN 值，實驗包含如下兩種情形：

①將上述100 個站點的數(shù)據(jù)作為已知樣本數(shù)據(jù)（單位km2存在2‰個采樣點），估計剩余非邊界點外站點的降水值，下文簡稱為實驗1。

②從467 個站點除邊界點外每次取出一個點作為待插值點（單位km2有12‰個點，實驗2 樣本數(shù)據(jù)的密度較實驗1 大），從剩余樣本點中抽取OK、AIDW 與OOK 插值方法需要的參考點來預(yù)測待插值點的降水值，簡稱為實驗2。

實驗1：研究區(qū)內(nèi)最大降雨量為58.5mm，最低值為1mm，均值為18mm。首先利用GS+擬合出該區(qū)域的半變異函數(shù)模型，在給定“步長”與“最大滯后距”后，對樣本數(shù)據(jù)的Z 值進(jìn)行對數(shù)轉(zhuǎn)換，選擇RSS 值最低的球狀模型（C0=0.071，C0+C=0.554，A0=71800，r2=0.847，RSS=0.0477），可得塊金系數(shù)為12%，低于25%[15]。將OOK、OK（4-16 個參考樣點）、AIDW 應(yīng)用于該數(shù)據(jù)集的插值實驗，得到如圖4。

圖4（左- 右，上- 下為a,b,c,d）

從圖（4-a）可以發(fā)現(xiàn)，OK 插值算法的MAE 誤差隨著參考點的增加而逐步降低，在參考點數(shù)達(dá)到12 以上時變化趨于平緩，在14 個參考點時處于全局最低點。OK 在選擇4-11 個參考點時的MAE 誤差大于OOK，在12 以上略低于OOK 算法，約為0.06mm。

從圖（4-b）可以發(fā)現(xiàn)，OK 插值算法的RMSE 誤差同樣隨著參考點的增多而逐步減少，在參考點數(shù)為12 后呈現(xiàn)平緩的趨勢，在點數(shù)為13 時全局最低。在OK 插值方法選擇5 個參考點后，AIDW 的RMSE 誤差大于OK，同時大于OOK。OOK 插值算法的誤差比OK 選擇4-10 個參考點要低。

從圖（4-c）可以發(fā)現(xiàn)，OK 插值算法的MAX 誤差基本隨參考點數(shù)的增多而下降，在參考點大于5 以后，OK 與AIDW 插值算法的MAX 誤差要優(yōu)于OOK 算法。經(jīng)過分析發(fā)現(xiàn)，OOK、AIDW與OK 插值算法MAX 值最大的待插值點的真實降水量為研究區(qū)域最大-51.8mm，其“一階鄰近點”中存在著近距離的局部最低點-19.2mm，可見其空間異質(zhì)性顯著。由于OK 與OOK 使用的半變異函數(shù)模型不能反映此局部變化，因此造成此點插值誤差較大。但OK 算法隨著搜索范圍的擴大，通過引入更多與待插值點具有弱相關(guān)性的參考點，對該點的插值誤差有輕微的抑制作用，OK 與OOK 的MAX 誤差的極差為1.2mm，發(fā)生在OK 算法選取14 個參考樣點。

從圖（4-d）可知，OK 插值方法的最小值呈現(xiàn)不同振幅的震蕩，其中MIN 值總體上最優(yōu)的插值方法為OOK，其次AIDW，OK最后。

從上述實驗可以看出：1）OOK 算法在OK 算法選取參考點數(shù)為4-11 點時占據(jù)絕對的優(yōu)勢，OK 算法在參考點大于12 時相較于其余兩種插值算法精度更高；2）OOK 插值算法的精度總體上優(yōu)于AIDW。

實驗2：首先將研究區(qū)內(nèi)樣本數(shù)據(jù)降雨量進(jìn)行均方根轉(zhuǎn)換，擬合的半變異函數(shù)模型為球狀模型（C0=0.001，C0+C=1.99，A0=83700，r2=0.959，RSS=0.21），可知塊金系數(shù)為0.05%，相較于實驗1 的12%，空間自相關(guān)性更強。將OOK、OK（4-16 個參考樣點）、AIDW 應(yīng)用于插值實驗，得到如圖5。

從圖（5-a）中可以看出，OK 插值算法參考點為4-6 點時，MAE 誤差逐步降低，其中在參考點為6 時，在3 種插值算法中最低。但隨著參考點的數(shù)量逐漸增多，MAE 值呈現(xiàn)類對數(shù)方式上升，最后在參考點為13 時呈現(xiàn)平穩(wěn)的趨勢，在16 個點時，MAE 在3 種插值算法中誤差最大。OOK 插值方法在OK 參考點數(shù)量為6 時，其MAE 誤差略大于OK。AIDW 在OK 的參考點為9 后其MAE 值低于OK，同時從圖可以發(fā)現(xiàn)AIDW 比OOK 的MAE 誤差大。

圖5（左-- 右，上-- 下為a,b,c,d）

從圖中（5-b）可以發(fā)現(xiàn)，OK 插值算法的RMSE 誤差在插值點數(shù)為6-9 個時呈現(xiàn)"J" 形增長趨勢，隨后趨于穩(wěn)定，在16 個點時誤差在3 種算法中最高，在參考點為6 時，RMSE 值最低。AIDW 在參考點為9 后RMSE 值低于OK。OOK 的誤差比AIDW低，僅在參考點為6 時大于OK。

從圖（5-c）中可以發(fā)現(xiàn)，OK 插值算法的MAX 值5、6 點時全局最低，此后總體趨勢呈現(xiàn)逐步增大并最終趨于穩(wěn)定，AIDW 的MAX 誤差在3 種插值算法中最大。OOK 在OK 選擇參考點為5、6 時，其MAX 值略大，約為0.3mm。

從圖（5-d）可以看出，OK 插值算法的MIN 值隨著點數(shù)的增多呈現(xiàn)上下震蕩，OOK 的MIN 比AIDW 略低，在參考點為8--15 時較OK 算法略低，但總體而言，3 種方法的MIN 相差都較小，極差約為0.03mm。

由可知：a.OOK 插值方法僅在OK 選擇6 個參考點時，其插值精度略低于OK，其余在3 種插值方法中均占據(jù)絕對優(yōu)勢；b.OOK 插值算法的插值精度比AIDW 要高，這是由于該數(shù)據(jù)集空間自相關(guān)性強烈且該區(qū)域的變異函數(shù)模型較好的反映了空間數(shù)據(jù)的分布特征[16]。

綜合以上2 個實驗可以發(fā)現(xiàn)：a. 沒有哪種算法在任何情形下都占據(jù)絕對的優(yōu)勢；b.在空間自相關(guān)性強的數(shù)據(jù)中，顧及樣本分布不均勻性的OOK 插值方法的插值精度（MAE、RMSE、MAX、MIN）較高，尤其在樣本數(shù)據(jù)密度較大時。

4 結(jié)論

本文針對空間樣本點數(shù)據(jù)的分布存在不均勻性的情形，提出了一種自適應(yīng)選取參考樣點的插值方法OOK，通過與經(jīng)典的OK、AIDW 對比發(fā)現(xiàn)：OOK 插值算法無需用戶指定參考樣點的數(shù)量，同時能夠顧及樣點分布的方向不均勻性，能夠提高插值精度。

但自適應(yīng)選取參考點的Kriging 插值方法也存在一定的局限性：①從本文實驗1 的MAX 值誤差分析發(fā)現(xiàn)，方法的抗異質(zhì)性影響的能力弱；②目前Kriging 插值算法的半變異函數(shù)普遍是通過不斷試驗比較得來，加大普通用戶使用Kriging算法的難度。因此未來關(guān)于Kriging 插值方法的研究將聚焦于將ML（機器學(xué)習(xí)）用于空間異質(zhì)性的判別與半變異函數(shù)模型的自動求取，進(jìn)一步提高kriging 插值算法的精準(zhǔn)度與自動化程度。