一種自適應選取參考點的Kriging 插值方法研究

鄧顯石

（湖南環境生物職業技術學院園林學院，湖南 衡陽421005）

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

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

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

圖7 選擇特征隨吸合次數增加的變化曲線樣本點采樣分布不均勻對最終插值結果造成的影響，同時不需要用戶指定參考點的數目。

1 普通Kriging 模型

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

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

2 OOK 插值方法

2.1 定義1：將待插值點與研究區域內已知的樣本點建立不規則三角網，將與待插值點形成鄰接拓撲關系的樣本點稱為“一階鄰近點”，如圖（1）所示。由圖可以發現，采用“一階鄰近點”作為待插值點的參考點，能夠在一定程度上降低樣本點在一側發生聚集對參考點選擇造成的影響。

圖1 一階鄰近點

2.2 技術路線

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

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

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

圖2 OOK 插值流程

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

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

3 實驗對比分析

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

圖3 氣象站點分布圖

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結論

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

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