基于Kriging方法的負偏態(tài)分布數(shù)據(jù)的空間插值研究

黃 靜，王愛倩，翟世龍

（1.新疆農(nóng)業(yè)大學 水利與土木工程學院，新疆 烏魯木齊 830052；2.新疆克孜爾水庫管理局，新疆 庫車 842000；3.新疆地震局，新疆 烏魯木齊830011）

Kriging插值方法又稱空間局部插值法，是以變異函數(shù)理論和結構分析為基礎，在有限區(qū)域內(nèi)對區(qū)域化變量進行無偏最優(yōu)估計的一種方法［1-2］。該理論和方法建立在二階平穩(wěn)或內(nèi)蘊假設的基礎上，所以要求進行分析計算的數(shù)據(jù)服從正態(tài)分布，但在實際應用中，數(shù)據(jù)常具有異常值，高偏度以及非正態(tài)分布性質(zhì)的特點，這對于變異函數(shù)擬合及插值的穩(wěn)健性有著極大的影響［3］。目前在對實際數(shù)據(jù)進行穩(wěn)健處理方面做了很多研究，最為有效的方法是將非正態(tài)的原始數(shù)據(jù)轉(zhuǎn)化為正態(tài)或近似正態(tài)分布的變換［4］。變換的方法有直接變換，如取倒數(shù)、取平方根、取對數(shù)，以及變換能夠較強的冪變換和Johnson變換。但實際中，正態(tài)變換并不是對所有的數(shù)據(jù)都奏效，比如存在趨勢的拖尾負偏態(tài)數(shù)據(jù)，無法進行正態(tài)轉(zhuǎn)換時，就要先對數(shù)據(jù)進行趨勢分析，采用全局多項式擬合配合Kriging插值方法來達到穩(wěn)健性插值的目的。

1 數(shù)據(jù)基本情況

本文選取的樣本數(shù)據(jù)是某水庫庫底高程數(shù)據(jù)，庫區(qū)范圍：南北方向約長13km，東西方向約長8km，樣本數(shù)共計466個，采樣間隔不等，間隔最小約150m，最大約1 000m左右。

2 數(shù)據(jù)基本統(tǒng)計分析

對466個樣本數(shù)據(jù)進行直方圖統(tǒng)計分析，如圖1所示，在直方圖中現(xiàn)實樣本數(shù)據(jù)是存在拖尾的負偏態(tài)數(shù)據(jù)，峰值為2.55，偏度為-0.74，對其進行KS正態(tài)檢驗，P值小于0.01，說明該樣本數(shù)據(jù)不符合正態(tài)分布，見圖2。

3 數(shù)據(jù)趨勢分析和全局趨勢擬合

對該樣本數(shù)據(jù)進行正態(tài)變換，仍無法通過K-S正態(tài)檢驗。對數(shù)據(jù)進行分析，采用趨勢分析來判斷樣本數(shù)據(jù)在特定方向是否存在趨勢，若存在趨勢，先選擇一個擬合效果最好的多項式對散點進行內(nèi)插，然后再對殘差隨機短程變異進行Kriging插值建立殘差模型，見圖3。

在圖3中，X軸代表實際地理位置的南北方向，Y軸代表實際地理位置的東西方向，因此從該圖中可以看出，該樣本數(shù)據(jù)南北方向呈現(xiàn)南高北低趨勢，在東西方向呈現(xiàn)兩頭高中間低的趨勢，因此，采用二次多項式對該數(shù)據(jù)進行全局趨勢擬合，并得到殘差數(shù)據(jù)。

圖1 樣本高程數(shù)據(jù)頻率直方圖

圖2 樣本高程數(shù)據(jù)正態(tài)分布概率圖

圖3 樣本高程數(shù)據(jù)南北-東西方向投影趨勢

4 殘差數(shù)據(jù)的隨機短程變異函數(shù)模型參數(shù)估計

對殘差數(shù)據(jù)計算變異函數(shù)值，通過變異函數(shù)模型進行模擬，目的在于確定一個最佳的擬合模型，最終計算出變異函數(shù)模型的幾個主要參數(shù)：塊金值、基臺值和變程。要想獲得精度較高的Kriging插值結果，變異函數(shù)模型的參數(shù)估計至關重要，首先要確定離散數(shù)據(jù)計算變異函數(shù)值的步長、最大步長和步數(shù)，其次要根據(jù)不同方向上的變異函數(shù)值分布情況，判斷數(shù)據(jù)是否存在各項異性，將此影響加入到模型中。

最大步長的設定一般是研究區(qū)域沿某一方向的最大尺度［5］，最大步長等于步長×步長組數(shù)。本文采用平均最鄰近法確定步長大小，步長為203m，步長組數(shù)為30，最大步長是6 090m，約是最大方向長度的一半。

在圖4中，變異函數(shù)表面圖上每個柵格單元用顏色進行編碼，從圖4中可以看出，在西北-東南方向的變異性要比南西-北東方向增加的略快。此時，選用的變異函數(shù)模型為tetraspherical，變程為2.5 km，基臺值為17.682，不存在塊金效應。下面對數(shù)據(jù)進一步進行方向效應分析，見圖5。

圖4 未考慮方向效應的變異函數(shù)散點圖及表面

圖5 考慮方向效應的變異函數(shù)圖及表面

在圖5中，散點圖顯示的是在某一方向的變異函數(shù)分布情況，在變異函數(shù)表面圖中橢圓的長軸表示變異函數(shù)在該方向有最大變程，其值為5.2km，超過這個變程殘差數(shù)據(jù)空間不相關，長軸方向表示變異函數(shù)在此方向有最小變程，其值為1.9km。當變異函數(shù)值達到穩(wěn)定水平時的值相同時，這就是函數(shù)模型的基臺值為18.507，塊金值為0.023 57。

5 Kriging插值結果及交叉驗證

得到變異函數(shù)擬合模型，并且該模型中考慮了趨勢面的剔除和方向效應，從而計算Kriging插值系數(shù)，對數(shù)據(jù)進行插值。

為了選擇最佳的變異函數(shù)模型來提高精度，分別采用不同的變異函數(shù)模型進行計算，并通過交叉驗證［6］結果進行比較，本文把最佳插值精度結果列出，選用的模型是Tetraspherical，見圖6。

圖6 Kriging插值交叉驗證結果

從圖6可以看出，選用Tetraspherical作為變異函數(shù)擬合模型進行Kriging插值，其交叉驗證預測精度結果是：平均預測誤差為-0.029 53，該值應接近于0；均方根標準誤差為1.001，其值應接近于1；均方根誤差為1.479，平均標準差為1.417，其值越小，說明精度越高，是選取最佳變異函數(shù)模型的依據(jù)。

6 結 論

本文對拖尾的負偏態(tài)數(shù)據(jù)進行深入分析，將數(shù)據(jù)剝離為趨勢項和殘差兩部分，分別采用二次多項式曲面擬合及Kriging插值的方法進行插值，在變異函數(shù)模型擬合時考慮方向效應，對模型參數(shù)計算進行優(yōu)化，其插值結果采用交叉驗證進行精度評定，最終選定Tetraspherical作為最終的變異函數(shù)模型，其插值精度最高。

通過對Kriging插值結果進行交叉驗證發(fā)現(xiàn)，預測精度較低的點主要分布在測區(qū)的邊緣地帶，其原因有二，一是樣本數(shù)據(jù)選用的是某水庫庫底的高程數(shù)據(jù)，其地形變化在庫邊變化較快；二是在采樣時，庫邊數(shù)據(jù)采樣間隔大，數(shù)據(jù)量少。在對部分空間數(shù)據(jù)進行插值時，邊緣數(shù)據(jù)的處理至關重要，采用什么模型和方法來對空間數(shù)據(jù)進行插值來提高邊緣數(shù)據(jù)的插值精度，是以后研究的重點。

［1］曾懷恩，黃聲享.基于Kriging方法的空間數(shù)據(jù)插值研究［J］.測繪工程，2007，16（5）：5-8.

［2］張梁，林韜，汪慶鋒.使用克里金插值法進行CGCS2000到海南地方坐標系的轉(zhuǎn)換［J］.測繪與空間地理信息，2014，37（9）：219-222.

［3］李曉輝，袁峰，白曉宇，等.典型礦區(qū)非正態(tài)分布土壤元素數(shù)據(jù)的正態(tài)變換方法對比研究［J］.地理與地理信息科學，2010，26（6）：102-105.

［4］陳道貴，胡乃聯(lián)，李國清.區(qū)域化變量非正態(tài)分布的穩(wěn)健性［J］.北京科技大學學報，2009，31（4）：412-417.

［5］吳學文，晏路明.普通Kriging法的參數(shù)設置及變異函數(shù)模型選擇方法［J］.地球信息科學，2007，9（3）：104-108.

［6］陳華鑫，姜藝，李碩，等.基于ArcGIS的交通量預測模型［J］.同濟大學學報，2010，38（8）：104-108.