李 佳 霖，樊 子 德，鄧 敏

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

?

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

李 佳 霖，樊 子 德，鄧 敏

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

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

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

0 引言

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

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

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

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

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

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

圖1 算法流程

Fig.1 Flow chart of RRIDW

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

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

(1)

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

(2)

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

1.2 趨勢(shì)面方程的計(jì)算

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

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

(3)

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

(4)

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

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

(5)

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

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

(6)

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

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

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

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

(7)

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

本實(shí)驗(yàn)中北京市35個(gè)空氣污染物監(jiān)測(cè)站的分布及PM2.5濃度數(shù)據(jù)空間異質(zhì)分區(qū)情況如圖2所示。由圖2可以發(fā)現(xiàn)，計(jì)算得到的北京市空氣污染物濃度空間分區(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é)果的誤差，分別將相對(duì)誤差繪制成柱狀圖。由圖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方法對(duì)北京市監(jiān)測(cè)站點(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)中，中國(guó)554個(gè)氣象監(jiān)測(cè)站點(diǎn)的分布情況和降雨量空間異質(zhì)分區(qū)結(jié)果如圖5所示。計(jì)算得到的中國(guó)平均降雨量空間分區(qū)呈現(xiàn)出從西北到東南降雨量逐漸增大的特點(diǎn)，符合中國(guó)氣候分區(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方法對(duì)中國(guó)554個(gè)氣象站點(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é)語(yǔ)

針對(duì)現(xiàn)有空間插值方法沒(méi)有充分考慮空間異質(zhì)性的影響，本文提出一種基于空間異質(zhì)分區(qū)的殘差反距離加權(quán)插值方法(RRIDW)。該方法綜合考慮兩類(lèi)空間異質(zhì)性對(duì)插值的影響，對(duì)研究區(qū)域進(jìn)行空間分區(qū)并剔除趨勢(shì)值之后，利用采樣點(diǎn)殘差計(jì)算待求點(diǎn)處的屬性值，使得插值方法的精度明顯提高。最后通過(guò)兩組實(shí)際數(shù)據(jù)驗(yàn)證了本文方法的優(yōu)越性和可行性。未來(lái)考慮顧及時(shí)間維的數(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é)報(bào)(信息科學(xué)版),2012,37(11):1285-1288.

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

[8] 彭思嶺.氣象要素時(shí)空插值方法研究[D].長(zhǎng)沙:中南大學(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

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

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

10.3969/j.issn.1672-0504.2015.05.006

P208

A

1672-0504(2015)05-0025-05