李 佳 霖，樊 子 德，鄧 敏

(中南大學(xué)地理信息系，湖南 長沙 410083)

?

基于空間異質(zhì)分區(qū)的殘差I(lǐng)DW插值方法

李 佳 霖，樊 子 德，鄧 敏

(中南大學(xué)地理信息系，湖南 長沙 410083)

空間插值可以利用已有觀測數(shù)據(jù)修補(bǔ)缺失的觀測數(shù)據(jù)，也可以利用離散數(shù)據(jù)構(gòu)建連續(xù)的表面數(shù)據(jù),但現(xiàn)有的空間插值方法沒有充分考慮空間數(shù)據(jù)的異質(zhì)性。該文提出一種基于空間異質(zhì)分區(qū)的殘差反距離加權(quán)插值方法(RRIDW)。首先根據(jù)采樣點(diǎn)屬性值對研究區(qū)域進(jìn)行空間異質(zhì)分區(qū);為了進(jìn)一步去除不同子區(qū)域內(nèi)的空間趨勢，對每個子區(qū)域計算趨勢面，進(jìn)而計算得到采樣點(diǎn)屬性值的異質(zhì)分區(qū)殘差，利用屬性值殘差進(jìn)行反距離加權(quán)插值;最后結(jié)合趨勢計算得到待求點(diǎn)處的空間插值結(jié)果。實(shí)驗(yàn)采用兩組實(shí)際PM2.5濃度數(shù)據(jù)和降雨量數(shù)據(jù)，運(yùn)用交叉驗(yàn)證方法對RRIDW方法與其他常用空間插值方法進(jìn)行對比分析，驗(yàn)證了該方法的優(yōu)越性和可行性。

空間異質(zhì)性；空間分區(qū)；空間插值；趨勢面分析

0 引言

利用離散點(diǎn)數(shù)據(jù)進(jìn)行空間插值得到連續(xù)表面數(shù)據(jù)，已被廣泛應(yīng)用于各類科學(xué)研究中[1]。為了修補(bǔ)缺失觀測數(shù)據(jù)并且獲得高精度空間表面模型，空間插值利用已有空間數(shù)據(jù)觀測站點(diǎn)的觀測值對待求點(diǎn)處的屬性值進(jìn)行插值估計，最終得到供人們使用的高精度表面數(shù)據(jù)[2]?，F(xiàn)有不同的空間插值方法對于不同領(lǐng)域的空間數(shù)據(jù)具有不同的適應(yīng)性[3]。在實(shí)際應(yīng)用中，通常采用交叉驗(yàn)證法[4]評價不同空間插值方法在不同領(lǐng)域中的插值結(jié)果。

近幾年，部分學(xué)者從理論對經(jīng)典空間插值方法(如反距離加權(quán)方法IDW和普通克里金方法OK)進(jìn)行了優(yōu)化，還有一些學(xué)者針對各種應(yīng)用需求發(fā)展了一些空間插值方法[5-8]。例如，Nalder等[9]將多元線性回歸方法和IDW方法結(jié)合提出一種梯度下降反比法(GIDW)，通過對月均氣溫和降水量的實(shí)驗(yàn)對比，發(fā)現(xiàn)GIDW方法的插值結(jié)果精度高于其他空間插值方法。Lu等[10]根據(jù)待求點(diǎn)鄰域內(nèi)樣本點(diǎn)的空間分布模式確定距離衰減參數(shù)，提出了AIDW方法，并通過實(shí)驗(yàn)證明了可變的距離衰減參數(shù)得到的插值結(jié)果更準(zhǔn)確。Reyes等[11]將LUR模型和BME結(jié)合，對空氣污染物濃度數(shù)據(jù)進(jìn)行建模，實(shí)驗(yàn)結(jié)果優(yōu)于普通克里金方法，但實(shí)際中難以同時獲得多種氣象環(huán)境等數(shù)據(jù)和多源污染物濃度觀測數(shù)據(jù)。Luis等[12]通過四種插值方法對墨西哥城的多種空氣污染物濃度進(jìn)行空間插值，發(fā)現(xiàn)IDW方法和普通克里金方法的插值精度明顯高于另兩種方法。

空間數(shù)據(jù)中通常存在空間異質(zhì)性[13]，但現(xiàn)有的插值方法并沒有充分考慮空間異質(zhì)性對于插值結(jié)果產(chǎn)生的如下主要影響。一方面由于空間異質(zhì)性的影響，研究區(qū)域不同子區(qū)域之間的屬性值可能存在較大的差異，這樣空間上鄰近兩個點(diǎn)的屬性值也有可能相差很多。比如，由于秦嶺對季風(fēng)的阻擋，其北面的陜西省氣候干燥而其南面的四川省潮濕多雨。如果不進(jìn)行空間分區(qū)去除不同子區(qū)域之間空間異質(zhì)性的影響，那么待求點(diǎn)屬性值的插值估計就會受到屬性差異較大的另一子區(qū)域內(nèi)采樣點(diǎn)的影響。另一方面，即使在考慮了空間異質(zhì)性而得到的子區(qū)域內(nèi)部，由于分區(qū)并不能充分去除空間異質(zhì)性的影響，空間數(shù)據(jù)中仍然存在一定的趨勢。比如，秦嶺北面的陜西省年均降雨量也存在南多北少的趨勢。如果不將采樣點(diǎn)對應(yīng)的趨勢值剔除，得到屬性值殘差再進(jìn)行空間插值，那么采樣點(diǎn)屬性值中的趨勢項(xiàng)將會對插值結(jié)果造成不利的影響。

因此，本文提出一種基于空間異質(zhì)分區(qū)的殘差反距離加權(quán)插值方法(RRIDW),以避免將異質(zhì)偏差引入到插值結(jié)果中而對空間插值結(jié)果的精度造成不利影響。為驗(yàn)證本文提出算法的優(yōu)越性和可行性，文中分別以空氣污染物數(shù)據(jù)和氣象要素數(shù)據(jù)進(jìn)行實(shí)驗(yàn)，并使用本文方法與其他常用方法進(jìn)行插值分析和對比驗(yàn)證。

1 基于空間異質(zhì)分區(qū)的殘差I(lǐng)DW插值方法

對于不同類型的空間數(shù)據(jù)，為了消除空間異質(zhì)性對插值的影響，RRIDW插值方法首先根據(jù)采樣點(diǎn)的位置和屬性值進(jìn)行空間分區(qū)；考慮到分區(qū)后每個子區(qū)域內(nèi)部的屬性值仍然存在一定的趨勢，由各個子區(qū)域內(nèi)的采樣數(shù)據(jù)計算其趨勢面方程，并用采樣點(diǎn)屬性數(shù)據(jù)減去對應(yīng)的趨勢值得到采樣點(diǎn)處的殘差值；根據(jù)采樣點(diǎn)處的屬性殘差值，按照IDW方法計算待求點(diǎn)處的屬性殘差值；最后，將待求點(diǎn)處的趨勢值和殘差值相加得到插值結(jié)果。算法流程如圖1。

圖1 算法流程

Fig.1 Flow chart of RRIDW

1.1 空間異質(zhì)分區(qū)

受空間異質(zhì)性的影響，空間上鄰近的兩個點(diǎn)屬性值可能相差很多。本文根據(jù)采樣點(diǎn)數(shù)據(jù)對研究區(qū)域進(jìn)行空間異質(zhì)分區(qū),將大量的空間對象劃分到一定數(shù)量的空間連續(xù)區(qū)域中，使得這些區(qū)域內(nèi)部空間同質(zhì)[14]。使用兼顧空間位置和屬性雙重特性的REDCAP方法進(jìn)行空間分區(qū)。首先，按照一定的鄰近約束規(guī)則將區(qū)域內(nèi)的所有空間對象聚集成為一個空間鄰接樹。根據(jù)研究成果[14,15]，本文選擇聚集效果較好的參數(shù)，聚集方式選擇最小方差增量距離[16]的方法，約束參數(shù)選擇全階策略。對所有連接兩個簇的邊，根據(jù)最小方差增量距離進(jìn)行升序排列，依次將連接兩個簇的邊添加到空間鄰接樹中，直到將所有的空間對象包含在一個簇內(nèi)。然后，根據(jù)空間同質(zhì)性的增量對空間鄰接樹進(jìn)行劃分得到空間分區(qū)?？臻g異質(zhì)性表示方法如下：

(1)

對空間鄰接樹進(jìn)行劃分時，依次將一個空間鄰接子樹劃分為兩個，直到滿足數(shù)量要求。對空間子樹進(jìn)行劃分時，分別計算去除每條子樹的邊之后同質(zhì)性的增量，將其中同質(zhì)性增量最大時所對應(yīng)的邊去除。劃分子樹產(chǎn)生的同質(zhì)性增量計算如下：

(2)

式中：Ra、Rb分別為將空間分區(qū)R分割后的兩個空間分區(qū)。對空間鄰接樹進(jìn)行多次劃分即可得到指定數(shù)量的空間分區(qū)，用于下一步趨勢面方程的計算。

1.2 趨勢面方程的計算

進(jìn)行空間異質(zhì)分區(qū)之后得到的子區(qū)域內(nèi)仍然存在某種特定的趨勢，而傳統(tǒng)的反距離加權(quán)插值方法中并沒有考慮這種趨勢。因而，本文方法針對每一個空間分區(qū)，根據(jù)樣本點(diǎn)屬性值與空間坐標(biāo)的關(guān)系，采用多項(xiàng)式回歸方法，按照最小二乘原理對數(shù)據(jù)點(diǎn)進(jìn)行擬合，得到每個分區(qū)對應(yīng)的趨勢面方程。

首先對采樣數(shù)據(jù)進(jìn)行探索性分析，判斷待插值的數(shù)據(jù)中存在的趨勢類型，然后選擇合適的趨勢面方程進(jìn)行擬合。本文根據(jù)不同分區(qū)內(nèi)采樣點(diǎn)的探索性分析結(jié)果，選擇與其特征較為接近的二元二次趨勢面進(jìn)行擬合，如式(3)所示：

(3)

依次將每一個空間分區(qū)中采樣點(diǎn)的屬性值，減去對應(yīng)位置上的趨勢值，得到該采樣點(diǎn)的屬性值殘差，用于反距離加權(quán)空間插值計算，公式如下：

(4)

1.3 待求點(diǎn)屬性值的計算

采樣點(diǎn)屬性值殘差通過空間異質(zhì)分區(qū)和趨勢剔除的方式，充分顧及了空間異質(zhì)性的影響?，F(xiàn)根據(jù)采樣點(diǎn)屬性值的殘差，按照反距離加權(quán)法計算待求點(diǎn)處屬性值的殘差。反距離加權(quán)法(IDW)是以待求點(diǎn)與采樣點(diǎn)之間的距離為權(quán)重進(jìn)行插值的方法，其權(quán)重貢獻(xiàn)與距離呈反比，公式如下：

(5)

式中：Zi表示待求點(diǎn)鄰域內(nèi)第i個采樣點(diǎn)的屬性值殘差，di為待求點(diǎn)與鄰域內(nèi)第i個采樣點(diǎn)的距離，α為距離衰減參數(shù)。

最后，將計算得到的待求點(diǎn)處屬性值的殘差與其對應(yīng)位置的趨勢面屬性值相加，即為待求點(diǎn)處的空間插值結(jié)果，公式如下：

(6)

式中：Zest即為待求點(diǎn)處插值得到的屬性值。

2 實(shí)驗(yàn)算例分析

為了驗(yàn)證本文方法的優(yōu)越性和可行性，實(shí)驗(yàn)一采用2013年12月至2014年6月北京市35個空氣污染物監(jiān)測站PM2.5濃度均值數(shù)據(jù)；實(shí)驗(yàn)二采用1984-2009年26年中我國554個降雨量觀測站點(diǎn)平均降雨量數(shù)據(jù)。實(shí)驗(yàn)對比經(jīng)典反距離加權(quán)方法(IDW)、普通克立金方法(OK)、基于分區(qū)的反距離加權(quán)法(RIDW)、基于趨勢面的反距離加權(quán)法(TIDW)以及本文提出的基于空間異質(zhì)分區(qū)的殘差反距離加權(quán)插值方法(RRIDW)。由于距離衰減參數(shù)對于插值結(jié)果的影響較大，所以為了保證實(shí)驗(yàn)的準(zhǔn)確性，在用對比方法進(jìn)行插值時，距離衰減參數(shù)分別采用最常用的1、2、3進(jìn)行實(shí)驗(yàn)。實(shí)驗(yàn)結(jié)果中方法名后的數(shù)字(例如IDW_1)表示采用此種方法時的距離衰減參數(shù)α的值。

最后采用留一法進(jìn)行交叉驗(yàn)證，其中平均絕對誤差(MAE)和均方根誤差(RMSE)兩個常用指標(biāo)可用來表達(dá)空間插值結(jié)果的絕對誤差，本文另外采用相對誤差(RE)指標(biāo)表達(dá)空間插值結(jié)果的相對誤差，計算公式為：

(7)

2.1PM2.5濃度數(shù)據(jù)空間插值驗(yàn)證

本實(shí)驗(yàn)中北京市35個空氣污染物監(jiān)測站的分布及PM2.5濃度數(shù)據(jù)空間異質(zhì)分區(qū)情況如圖2所示。由圖2可以發(fā)現(xiàn)，計算得到的北京市空氣污染物濃度空間分區(qū)呈現(xiàn)出南部污染物濃度較高、北部污染物濃度較低的特點(diǎn)，符合實(shí)際情況。

用五種空間插值方法進(jìn)行實(shí)驗(yàn)并得出結(jié)果，進(jìn)行交叉驗(yàn)證得到各項(xiàng)誤差指標(biāo)，如表1所示。

為直觀表現(xiàn)各種空間插值方法實(shí)驗(yàn)結(jié)果的誤差，分別將相對誤差繪制成柱狀圖。由圖3可以發(fā)現(xiàn)，IDW方法插值誤差最大，TIDW方法與OK方法

圖2 北京市空氣污染物濃度異質(zhì)分區(qū)結(jié)果

Fig.2Spatialheterogeneitysub-regionoftheconcentrationofPM2.5inBeijing

表1 五種空間插值方法的實(shí)驗(yàn)結(jié)果

Table 1 Experimental results of five different spatial interpolation methods

MAE(ug/m3)RMS(ug/m3)RE(%)IDW＿17．729．8910．40IDW＿27．829．7110．21IDW＿38．039．6610．16RIDW＿16．067．618．00RIDW＿26．157．758．15RIDW＿36．367．898．30TIDW＿16．047．628．01TIDW＿26．598．078．49TIDW＿37．008．368．79RRIDW＿13．745．125．39RRIDW＿24．095．665．95RRIDW＿34．325．906．21OK6．237．968．37

精度相近，顧及空間異質(zhì)性的RIDW方法的插值效果略優(yōu)于OK方法，均明顯優(yōu)于IDW方法，而本文提出的RRIDW方法的插值結(jié)果相比RIDW方法精度大幅提高。最后給出采用RRIDW方法對北京市監(jiān)測站點(diǎn)平均PM2.5濃度數(shù)據(jù)的插值結(jié)果(圖4)。

圖3 五種空間插值方法的實(shí)驗(yàn)結(jié)果

Fig.3 Results of experiments by five different spatial interpolation methods

2.2 降雨量數(shù)據(jù)的空間插值驗(yàn)證

本實(shí)驗(yàn)中，中國554個氣象監(jiān)測站點(diǎn)的分布情況和降雨量空間異質(zhì)分區(qū)結(jié)果如圖5所示。計算得到的中國平均降雨量空間分區(qū)呈現(xiàn)出從西北到東南降雨量逐漸增大的特點(diǎn)，符合中國氣候分區(qū)的實(shí)際情況。用五種空間插值方法進(jìn)行實(shí)驗(yàn)并得出結(jié)果，進(jìn)行交叉驗(yàn)證得到各項(xiàng)誤差指標(biāo)，如表2所示。

圖4 PM2.5濃度插值結(jié)果

Fig.4 Spatial interpolation results of PM2.5

圖5 降雨量空間異質(zhì)分區(qū)結(jié)果

Fig.5 Spatial heterogeneity sub-region of the precipitation

表2 五種空間插值方法的實(shí)驗(yàn)結(jié)果

Table 2 Experimental results of five different spatial interpolation methods

MAE(ug/m3)RMS(ug/m3)RE(%)IDW＿199．54167．8219．33IDW＿298．42167．9219．34IDW＿399．43170．1419．60RIDW＿193．05157．7918．18RIDW＿291．65158．4118．25RIDW＿392．14160．9618．54TIDW＿197．55166．9219．23TIDW＿297．18166．8919．23TIDW＿398．39168．8819．46RRIDW＿188．52154．6617．82RRIDW＿287．98155．4517．91RRIDW＿388．87157．8918．19OK97．12163．8718．88

為直觀表現(xiàn)各種空間插值方法實(shí)驗(yàn)結(jié)果的誤差，將五種插值方法的誤差繪制柱狀圖，如圖6所示。分析圖6可得到與實(shí)驗(yàn)一相似的空間插值結(jié)果，IDW方法插值誤差最大，TIDW方法插值結(jié)果精度略高于IDW方法，普通克里金方法的空間插值精度比TIDW方法略高，顧及空間異質(zhì)性的RIDW方法的插值效果優(yōu)于普通克里金方法，而本文提出的RRIDW方法的插值精度相比RIDW方法大幅提高。最后給出采用本文RRIDW方法對中國554個氣象站點(diǎn)平均降雨量數(shù)據(jù)的插值結(jié)果(圖7)。

圖6 五種空間插值方法的實(shí)驗(yàn)結(jié)果

Fig.6 Results of experiments by five different spatial interpolation methods

圖7 降水?dāng)?shù)據(jù)插值結(jié)果

Fig.7 Spatial interpolation results of rainfall

3 結(jié)語

針對現(xiàn)有空間插值方法沒有充分考慮空間異質(zhì)性的影響，本文提出一種基于空間異質(zhì)分區(qū)的殘差反距離加權(quán)插值方法(RRIDW)。該方法綜合考慮兩類空間異質(zhì)性對插值的影響，對研究區(qū)域進(jìn)行空間分區(qū)并剔除趨勢值之后，利用采樣點(diǎn)殘差計算待求點(diǎn)處的屬性值，使得插值方法的精度明顯提高。最后通過兩組實(shí)際數(shù)據(jù)驗(yàn)證了本文方法的優(yōu)越性和可行性。未來考慮顧及時間維的數(shù)據(jù)將此插值方法擴(kuò)展至?xí)r空插值研究。

[1] LI J,ANDREW D H.A review of comparative studies of spatial interpolation methods in environmental sciences:Performance and impact factors[J].Ecological Informatics,2011(6):228-241.

[2] SMITH M J D,GOODCHILD M F,LONGLEY P.Geospatial Analysis:A Comprehensive Guide to Principles,Techniques and Software Tools[M].Troubador Publishing Ltd,2007.

[3] LI J,ANDREW D H.Spatial interpolation methods applied in the environmental sciences:A review[J].Environmental Modelling & Software,2014,53(3):173-189.

[4] EFRON B,GONG G.A leisurely look at the bootstrap,the jackknife,and cross-validation[J].The American Statistician,1983,37(1):36-48.

[5] SHIODE N,SHIODE S.Street-level spatial interpolation using network-based IDW and ordinary Kriging[J].Transactions in GIS,2011,15(4):457-477.

[6] 易琳,袁林旺,羅文,等.顧及V-鄰域結(jié)構(gòu)的局部保形插值算法[J].武漢大學(xué)學(xué)報(信息科學(xué)版),2012,37(11):1285-1288.

[7] 李莎,舒紅,徐正全.利用時空Kriging進(jìn)行氣溫插值研究[J].武漢大學(xué)學(xué)報(信息科學(xué)版),2012,37(2):237-241.

[8] 彭思嶺.氣象要素時空插值方法研究[D].長沙:中南大學(xué),2010.

[9] NALDER I A,WEIN R W.Spatial interpolation of climatic normals:Test of a new method in the Canadian boreal forest[J].Agricultural and Forest Meteorology,1998,92(4):211-225.

[10] LU G Y,WONG D W.An adaptive inverse-distance weighting spatial interpolation technique[J].Computers & Geosciences,2008,34(9):1044-1055.

[11] REYES J M,SERRE M L.An LUR/BME framework to estimate PM 2.5 explained by on road mobile and stationary sources[J].Environmental Science & Technology,2014,48(3):1736-1744.

[12] RIVERA-GONZLEZ L O,ZHANG Z,SNCHEZ B N,et al.An assessment of air pollutant exposure methods in Mexico City,Mexico[J].Journal of the Air & Waste Management Association,2015,65:581-591.

[13] WANG J F,CHRISTAKOS G,HU M G.Modeling spatial means of surfaces with stratified nonhomogeneity[J].Geoscience and Remote Sensing,IEEE Transactions on,2009,47(12):4167-4174.

[14] GUO D.Regionalization with dynamically constrained agglomerative clustering and partitioning(REDCAP)[J].International Journal of Geographical Information Science,2008,22(7):801-823.

[15] GUO D,WANG H.Automatic region building for spatial analysis[J].Transactions in GIS,2011,15(s1):29-45.

[16] HAGGETT P,CLIFF A D,FREY A.Locational analysis in human geography[M].London:Edward Arnold,1977.

Residual Inverse Distance Weighting Spatial Interpolation Method Based on Spatial Heterogeneity Sub-region

LI Jia-lin,FAN Zi-de,DENG Min

(DepartmentofGeo-informatics,CentralSouthUniversity,Changsha410083,China)

Spatial interpolation methods,on one hand,can utilize the observation data to complete the missing data,and on the other hand,can also convert the discrete data into a continuous surface.However,the existing spatial interpolation methods lack an adequate consideration of the effect of spatial heterogeneity in spatial data.Therefore,the residual inverse-distance weighting spatial interpolation method based on spatial heterogeneity sub-region was presented in the paper.First of all,according to the sample point attribute value,the research area is divided into spatial heterogeneity sub-regions.And then calculate the trend surface to obtain the residuals of the sample points′ attributes,which are used for inverse distance weighted interpolation in the next step.Finally,the spatial interpolation result is estimated by combining the trend surface values and the interpolated residuals.At the end of this paper,two sets of actual PM2.5 concentrations and rainfall data are adopted to verify the effectiveness of the method proposed by comparing RRIDW with the commonly used spatial interpolation methods.

spatial heterogeneity;space partition;spatial interpolation;trend surface analysis

2015-03-31；

2015-06-22

國家863計劃項(xiàng)目(2013AA122301)；高等學(xué)校博士點(diǎn)專項(xiàng)科研基金項(xiàng)目(20110162110056)；湖南省博士生優(yōu)秀學(xué)位論文資助項(xiàng)目(CX2014B050)

李佳霖(1992-)，男，碩士研究生，主要從事時空插值方法的研究工作。E-mail：garlic_lee@foxmail.com

10.3969/j.issn.1672-0504.2015.05.006

P208

A

1672-0504(2015)05-0025-05