黃 靜，王愛倩，翟世龍

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖4 未考慮方向效應(yīng)的變異函數(shù)散點(diǎn)圖及表面

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

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

5 Kriging插值結(jié)果及交叉驗(yàn)證

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

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

圖6 Kriging插值交叉驗(yàn)證結(jié)果

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

6 結(jié) 論

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

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

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

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

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

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

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

［6］陳華鑫，姜藝，李碩，等.基于ArcGIS的交通量預(yù)測(cè)模型［J］.同濟(jì)大學(xué)學(xué)報(bào)，2010，38（8）：104-108.