深圳市氣溫空間插值方法比較

何輪凱 孫立群 李晴嵐 謝坤 王書欣 王德立 徐倩

（1 中國(guó)科學(xué)院深圳先進(jìn)技術(shù)研究院，深圳 518055；2 深圳市氣象局，深圳518040；3 香港大學(xué)，香港）

0 引言

氣象數(shù)據(jù)作為環(huán)境因子，是地學(xué)和氣候模型等相關(guān)研究的重要參數(shù)[1]。理論上氣象要素的空間分布可以通過(guò)高密度站網(wǎng)采集，然而實(shí)際上氣象觀測(cè)站點(diǎn)的空間分布稀疏異質(zhì)，觀測(cè)序列長(zhǎng)短不一，為網(wǎng)格化的模型應(yīng)用帶來(lái)了許多不確定性[2-3]。因此，利用有限的站點(diǎn)進(jìn)行合理的空間分布插值計(jì)算不僅可以提高模型運(yùn)行效率，更有助于氣象預(yù)報(bào)結(jié)果的精細(xì)化顯示[4]。常見(jiàn)的插值方法包括反距離權(quán)重插值法（Inverse Distance Weighted，簡(jiǎn)稱IDW）、局部多項(xiàng)式插值法（Local Polynomial Interpolation，簡(jiǎn)稱LPI）、克里金插值法和薄板樣條函數(shù)法（Thin Plate Spline，簡(jiǎn)稱TPS）等[5]。由于氣象要素的變化較為復(fù)雜，不同插值方法插值得出的結(jié)果有較大的差別，不存在絕對(duì)最優(yōu)的插值方法，但存在在特定環(huán)境區(qū)域中的最優(yōu)插值方法[6-7]。如何根據(jù)一個(gè)區(qū)域的地形、氣候變化特征、站點(diǎn)分布等選擇最合適的氣溫插值方法對(duì)于氣溫?cái)?shù)據(jù)進(jìn)行格點(diǎn)化表達(dá)是一個(gè)很重要的問(wèn)題。

目前，薄板樣條函數(shù)法與克里金法在考慮到地形因素影響的氣象數(shù)據(jù)空間插值中具有較高的準(zhǔn)確率和效率[8]。薄板樣條函數(shù)法和克里金法只將空間分布作為觀測(cè)數(shù)據(jù)的函數(shù)，不需要其物理過(guò)程與先驗(yàn)知識(shí)，有效地提高了插值的準(zhǔn)確度[9]。Hutchinson等[10]從理論上對(duì)比了薄板樣條函數(shù)法和克里金法的不同，認(rèn)為當(dāng)樣條函數(shù)法的粗糙度穩(wěn)定，且克里金法選擇了適當(dāng)?shù)陌胱儺惡瘮?shù)時(shí)，兩者均能獲得較好的插值效果，但是實(shí)際應(yīng)用中數(shù)據(jù)結(jié)構(gòu)更加簡(jiǎn)單、誤差估算更加精確的薄板樣條函數(shù)法更受推崇[11]。為了方便薄板樣條函數(shù)法的使用，Hutchinson等開(kāi)發(fā)了氣候數(shù)據(jù)薄板樣條函數(shù)法插值的專用軟件Anusplin[12]，該軟件可以引進(jìn)與氣象要素相關(guān)的協(xié)變量線性子模型，如高程等，基于不依賴于模型的廣義交叉驗(yàn)證（GCV）法，相對(duì)于克里金法更加穩(wěn)健，更能反映氣象要素在空間中自然分布的特征規(guī)律，已在國(guó)際上得到廣泛的應(yīng)用[13-14]。

目前，深圳市共有近200個(gè)可以正常工作的自動(dòng)氣象觀測(cè)站。在廣東省乃至全國(guó)都屬于自動(dòng)站點(diǎn)較為密集地區(qū)。然而，這些站點(diǎn)相對(duì)于深圳近2000 km2的面積來(lái)講仍然無(wú)法滿足某些氣象模式大區(qū)域、無(wú)縫隙、網(wǎng)格化的數(shù)據(jù)要求，并且在氣象預(yù)報(bào)結(jié)果的表達(dá)上也無(wú)法做到千米級(jí)的精細(xì)化呈現(xiàn)。基于以上現(xiàn)狀，選用較好的插值方法，建立深圳市氣溫格點(diǎn)化數(shù)據(jù)集十分有必要。本研究將以深圳市128個(gè)氣象站點(diǎn)的月平均氣溫?cái)?shù)據(jù)為研究對(duì)象，對(duì)包括薄板樣條函數(shù)法（TPS）、反距離插值法（IDW）、局部多項(xiàng)式插值法（LPI）、普通克里金和協(xié)同克里金在內(nèi)的5種插值法進(jìn)行插值效果比較，對(duì)比各類插值方法的優(yōu)劣及插值效果的影響因素，為深圳市氣溫插值的方法選擇提供科學(xué)依據(jù)，優(yōu)化插值效果。

1 數(shù)據(jù)與方法

1.1 數(shù)據(jù)來(lái)源

本研究使用的氣溫?cái)?shù)據(jù)集來(lái)源于深圳市氣象局，時(shí)間序列為2017年的1—12月，原始數(shù)據(jù)集共有195個(gè)站點(diǎn)，每個(gè)站點(diǎn)對(duì)應(yīng)的數(shù)據(jù)集包含了2017年每天的平均氣溫，為保證氣溫插值的質(zhì)量，站點(diǎn)中異常數(shù)據(jù)、缺失數(shù)據(jù)超過(guò)4%的站點(diǎn)均予以剔除。剔除有部分缺失數(shù)據(jù)、數(shù)據(jù)異常等的67個(gè)站點(diǎn)，剩余128個(gè)站點(diǎn)，由站點(diǎn)逐日的數(shù)據(jù)推算出每個(gè)月的平均氣溫后整理。每個(gè)站點(diǎn)數(shù)據(jù)包括站點(diǎn)的代號(hào)、經(jīng)度、緯度、海拔和月平均氣溫。

由于邊緣效應(yīng)，本研究使用的插值范圍比深圳市的實(shí)際范圍略大，為113.7°—114.7°E，22.3°—22.9°N，但插值時(shí)不增加深圳市外的站點(diǎn)。考慮到深圳市高程圖精度過(guò)高不利于計(jì)算效率的提高，精度過(guò)低插值效果不夠精確，調(diào)整精度為0.02°×0.02°（約為2 km×2 km）。

1.2 插值方法

本研究使用了五種插值方法，分別是IDW、LPI、普通克里金（OK）、協(xié)同克里金（Cokriging）和基于Anusplin軟件的TPS，以下分別對(duì)這五種方法進(jìn)行簡(jiǎn)單的介紹。

1.2.1 反距離權(quán)重法

IDW基于地理學(xué)第一定律-相似相近[15]，假設(shè)彼此距離較近的要素要比彼此距離較遠(yuǎn)的要素更相似。當(dāng)為任何未測(cè)量的位置預(yù)測(cè)值時(shí)，反距離權(quán)重法會(huì)采用預(yù)測(cè)位置周圍的測(cè)量值。與距離預(yù)測(cè)位置較遠(yuǎn)的測(cè)量值相比，距離預(yù)測(cè)位置最近的測(cè)量值對(duì)預(yù)測(cè)值的影響更大。反距離權(quán)重法假定每個(gè)測(cè)量點(diǎn)都有一種局部影響，而這種影響會(huì)隨著距離的增大而減小。由于這種方法為距離預(yù)測(cè)位置最近的點(diǎn)分配的權(quán)重較大，而權(quán)重卻作為距離的函數(shù)而減小，因此稱之為反距離權(quán)重法。該方法計(jì)算方式較為簡(jiǎn)單，但易受極值的影響。本研究使用的IDW插值法基于ArcGIS中Geostatistical Analyst模塊中確定性方法提供的IDW算法進(jìn)行插值計(jì)算。

1.2.2 局部多項(xiàng)式插值法

全局多項(xiàng)式插值法可以根據(jù)整個(gè)表面擬合多項(xiàng)式[16]，而LPI可以對(duì)位于指定重疊鄰域內(nèi)的多個(gè)多項(xiàng)式進(jìn)行擬合。通過(guò)使用大小和形狀、鄰域數(shù)量和部分配置，可以對(duì)搜索鄰域進(jìn)行定義。LPI僅使用既定鄰域內(nèi)的點(diǎn)對(duì)指定階數(shù)（0 階、1 階、2 階、3 階，等等）的多項(xiàng)式進(jìn)行擬合。鄰域相互重疊，每次預(yù)測(cè)所使用值是位于鄰域中心的擬合的多項(xiàng)式的值。本研究使用的LPI插值法基于ArcGIS中Geostatistical Analyst模塊中確定性方法提供的LPI算法進(jìn)行插值計(jì)算。

1.2.3 普通克里金

普通克里金是最早被提出和系統(tǒng)研究的克里金法[17]，并隨著地統(tǒng)計(jì)學(xué)的發(fā)展衍生出一系列變體和改進(jìn)算法。普通克里金是一個(gè)線性估計(jì)系統(tǒng)，適用于任何滿足各向同性假設(shè)的固有平穩(wěn)隨機(jī)場(chǎng)。各向同性假設(shè)下的固有平穩(wěn)隨機(jī)場(chǎng)中，數(shù)學(xué)期望與其位置無(wú)關(guān)，且協(xié)方差僅是點(diǎn)間距離的函數(shù)。通常隨機(jī)場(chǎng)的協(xié)方差函數(shù)是未知的，需要使用變異函數(shù)作為近似，此時(shí)變異函數(shù)也僅與點(diǎn)間距離有關(guān)。本研究使用的普通克里金插值法基于ArcGIS中Geostatistical Analyst模塊中地統(tǒng)計(jì)方法提供的Kriging算法進(jìn)行插值計(jì)算。

1.2.4 協(xié)同克里金

協(xié)同克里金是在普通克里金的基礎(chǔ)上發(fā)展起來(lái)的一種插值方法[18]，其將與變量高度相關(guān)的輔助信息作為協(xié)變量加入到克里金插值過(guò)程中，既考慮自變量的空間自相關(guān)性，又考慮了協(xié)變量與自變量之間的相關(guān)性。該方法建立于協(xié)同區(qū)域化變量理論的基礎(chǔ)上，通過(guò)建立變異函數(shù)模型和趨勢(shì)移除，利用主變量自相關(guān)性和主變量與其他所有變量的互相關(guān)性進(jìn)行更好的預(yù)測(cè)。本研究使用的協(xié)同克里金插值法基于ArcGIS中Geostatistical Analyst模塊中地統(tǒng)計(jì)方法提供的Cokriging算法進(jìn)行插值計(jì)算。

1.2.5 薄板樣條函數(shù)

薄板樣條函數(shù)插值法最早由Wahba[19]于1979年提出，1984年Hutchinson等對(duì)其進(jìn)行了改進(jìn)，以適用于更大的數(shù)據(jù)集[10]。Bates等將其拓展為局部薄板光滑樣條法[20]。薄板樣條函數(shù)法在插值過(guò)程中利用一個(gè)平面去擬合所有的站點(diǎn)，使得它可以通過(guò)站點(diǎn)組成的樣條，得到逼近所有控制點(diǎn)且曲率最小的光滑曲面。光滑參數(shù)由GCV的最小化或廣義最大似然的最小化確定，使用了最優(yōu)的光滑參數(shù)實(shí)現(xiàn)了逼真度和光滑度的最佳平衡，保證精度可靠且曲面光滑連續(xù)[21]。本研究中使用GCV的最小化確定光滑參數(shù)，評(píng)估誤差后使用優(yōu)化較好的二次樣條函數(shù)進(jìn)行插值研究。本研究使用的薄板樣條函數(shù)插值法基于Anusplin4.3版本的算法。

1.3 檢驗(yàn)方法

本研究采用交叉驗(yàn)證法中的留一驗(yàn)證法對(duì)5種插值方法的精度進(jìn)行評(píng)價(jià)，每次使用氣象站點(diǎn)中的一個(gè)站點(diǎn)作為驗(yàn)證資料，剩余的站點(diǎn)用作訓(xùn)練資料，進(jìn)行插值處理，通過(guò)測(cè)算實(shí)際值與預(yù)測(cè)值的誤差評(píng)價(jià)不同方法的插值精度。采用平均絕對(duì)誤差（MAE）和均方根誤差（RMSE）作為5種插值方法誤差的評(píng)價(jià)指標(biāo)，直接評(píng)價(jià)插值方法得出的結(jié)果與實(shí)際溫度的溫度差值：

式中，Z0(xi)和Z(xi)分別是第i個(gè)站點(diǎn)的預(yù)測(cè)值與實(shí)際值。平均絕對(duì)誤差與均方根誤差越小，說(shuō)明插值方法的準(zhǔn)確度越高。

基于反距離權(quán)重法、局部多項(xiàng)式法、普通克里金法、協(xié)同克里金法等4種方法，利用ArcGIS中Geostatistical Analyst模塊提供的交叉驗(yàn)證方法，使用所有數(shù)據(jù)對(duì)趨勢(shì)和自相關(guān)模型進(jìn)行估計(jì)。它會(huì)每次移除一個(gè)數(shù)據(jù)位置，然后預(yù)測(cè)關(guān)聯(lián)的數(shù)據(jù)值。在本研究中，交叉驗(yàn)證每次會(huì)省略一個(gè)點(diǎn)，然后使用剩余的127個(gè)點(diǎn)計(jì)算此位置的值。將省略點(diǎn)位置的預(yù)測(cè)值與實(shí)際值相比較。然后對(duì)第二個(gè)點(diǎn)重復(fù)此過(guò)程，以此類推。交叉驗(yàn)證會(huì)對(duì)所有點(diǎn)的測(cè)量值和預(yù)測(cè)值進(jìn)行比較。

基于Anusplin自帶的薄板樣條函數(shù)法，對(duì)每個(gè)月的數(shù)據(jù)集進(jìn)行分割，每次生成一個(gè)數(shù)據(jù)量為127個(gè)點(diǎn)的子集，以及剩余1個(gè)點(diǎn)的驗(yàn)證集，將子集輸入Anusplin軟件的插值模塊生成插值表面，使用Anusplin帶有的驗(yàn)證功能輸入驗(yàn)證集，得到驗(yàn)證集中的點(diǎn)對(duì)應(yīng)的預(yù)測(cè)值，重復(fù)此過(guò)程128次，可得到所有數(shù)據(jù)點(diǎn)的驗(yàn)證結(jié)果。

2 結(jié)果與分析

2.1 空間分布對(duì)比

基于圖1 中深圳市氣象站點(diǎn)分布與高程圖，以2017年7月為例，應(yīng)用TPS、IDW、LPI、OK、Cokriging等5種方法進(jìn)行平均氣溫的插值（圖2）。與IDW、LPI、OK等相比，將高程作為協(xié)變量進(jìn)行插值的TPS與Cokriging法得出的插值結(jié)果更加平滑。TPS法插值得出的溫度分布為22.34～29.51 ℃，Cokriging法插值得出的溫度分布為21.96～30.63 ℃，兩者結(jié)果相近。對(duì)其他月份的插值結(jié)果均繪圖顯示，結(jié)果表明，Cokriging和TPS插值的結(jié)果更加平滑，而精確度有待進(jìn)一步評(píng)估，以下闡述五種插值方法的誤差檢驗(yàn)

圖1 深圳市氣象站點(diǎn)分布（a）與高程圖（b） Fig. 1 Distribution of meteorological stations (a) and elevation map in Shenzhen (b)

2.2 逐月對(duì)比

之后交叉驗(yàn)證得到的MAE與RMSE作為判斷標(biāo)準(zhǔn)，判斷每個(gè)月份五種方法插值的精度高低。如圖3顯示，IDW法在10月時(shí)RMSE最大，為0.90 ℃，且MAE最大，為0.49 ℃。3月時(shí)RMSE與MAE最小，分別為0.67和0.36 ℃。LPI法的RMSE與MAE分別于10和9月達(dá)到峰值，為0.87和0.47 ℃。3月時(shí)LPI插值法的RMSE和MAE最小，為0.65和0.35 ℃。OK法的RMSE與MAE在10月達(dá)到最大，分別為0.90和0.49 ℃。3月時(shí)RMSE與MAE較小，分別為0.66和0.37 ℃。加入了高程作為協(xié)變量的Cokriging法的RMSE和MAE在10月達(dá)到最大，分別為0.83和0.56 ℃，3月時(shí)RMSE與MAE較小，為0.68和0.37 ℃。與IDW、LPI、OK、Cokriging法相比，TPS法的RMSE與MAE明顯較小，在9月達(dá)到峰值，分別為0.37和0.29 ℃，3月時(shí)RMSE與MAE較小，分別為0.25和0.18 ℃。在1—12月中每個(gè)月TPS法的MAE與RMSE均遠(yuǎn)小于其余四種方法。

圖2 基于5種插值方法的深圳市7月平均氣溫空間分布 （a）OK法插值結(jié)果；（b）LPI法插值結(jié)果；（c）IDW法插值結(jié)果；（d）Cokriging法插值結(jié)果；（e）TPS法插值結(jié)果 Fig. 2 Spatial distribution of July mean temperature in Shenzhen based on five interpolation methods(a) Interpolation results of OK; (b) Interpolation results of LPI; (c) Interpolation results of IDW; (d) Interpolation results of Cokriging; (e) Interpolation results of TPS

圖3 基于5種插值方法交叉驗(yàn)證誤差的RMSE、MAE逐月分析 Fig. 3 Monthly analysis of RMSE and MAE based on cross-validation errors of five interpolation methods

在插值精度上，RMSE與MAE的變化具有季節(jié)特征，在3月左右的插值精度明顯優(yōu)于10月左右的插值精度。RMSE結(jié)果顯示，TPS插值法在春夏季節(jié)的插值精度優(yōu)于秋冬季節(jié)。從誤差范圍來(lái)看，冬季的插值結(jié)果在負(fù)誤差方向表現(xiàn)比其他時(shí)間尺度明顯，表明利用Anusplin對(duì)研究區(qū)冬季氣溫?cái)?shù)據(jù)插值時(shí)，會(huì)更大概率出現(xiàn)低估的情況[22-23]。某些地區(qū)由于局地氣候的影響，在季節(jié)氣溫方面表現(xiàn)出特殊性，如在冬季這些地區(qū)的站點(diǎn)氣溫遠(yuǎn)高于同海拔區(qū)域[24]。本研究也發(fā)現(xiàn)，9月—次年2月的某些站點(diǎn)氣溫高于同海拔地區(qū)，導(dǎo)致Anusplin估算氣溫時(shí)出現(xiàn)低估的現(xiàn)象。這可能是夏季氣溫插值誤差小于冬季的原因。此外，不同站點(diǎn)下墊面性質(zhì)的不同，對(duì)氣溫插值的結(jié)果也有一定影響，會(huì)造成一定的誤差。

大氣層的主要直接熱源是地面，離地面越遠(yuǎn)，得到的地面輻射越少，氣溫也就越低。根據(jù)氣溫垂直遞減率[21]，每上升100 m，氣溫下降0.6 ℃，將高程作為協(xié)變量加入插值模型中，可以更好地模擬地形對(duì)氣溫的影響。

以深圳市最高峰梧桐山及周邊的站點(diǎn)為例（圖1中紅框內(nèi)的5個(gè)站點(diǎn)，5個(gè)站點(diǎn)的信息如表1所示），大梧桐及小梧桐站點(diǎn)地勢(shì)較高，沙頭角、蓮塘、明珠等站點(diǎn)地勢(shì)較低，高程圖顯示，該地區(qū)地形變化較大，其中小梧桐站點(diǎn)介于最高峰和平原之間，是地形變化最明顯的站點(diǎn)。將該站移除，以剩余站點(diǎn)的數(shù)據(jù)進(jìn)行插值，得到小梧桐站點(diǎn)的插值結(jié)果（圖4）。該站點(diǎn)的插值結(jié)果顯示，在12個(gè)月的氣溫逐月插值中基于Anusplin軟件的TPS方法準(zhǔn)確度明顯優(yōu)于其余方法。IDW、LPI、OK結(jié)果基本一致，偏高3～4 ℃。Cokriging方法結(jié)果誤差稍小，TPS方法的插值結(jié)果與實(shí)際觀測(cè)結(jié)果基本一致，平均每月的誤差為0.27 ℃。

表1 梧桐山及周邊站點(diǎn)信息 Table 1 Information of Wutong Mountain and its surrounding sites

2.3 影響因素分析

高程可能影響插值的最終結(jié)果，將每次插值的誤差絕對(duì)值（即觀測(cè)值減去預(yù)測(cè)值再取絕對(duì)值）與深圳市高程圖進(jìn)行對(duì)比，發(fā)現(xiàn)分布特征與海拔有一定的聯(lián)系。對(duì)每個(gè)月的插值誤差絕對(duì)值與該觀測(cè)點(diǎn)的海拔進(jìn)行Pearson相關(guān)分析（圖5），結(jié)果顯示，IDW、LPI、OK的插值誤差絕對(duì)值與海拔之間的相關(guān)系數(shù)較高，12個(gè)月的相關(guān)系數(shù)均大于0.75，呈現(xiàn)顯著正相關(guān)關(guān)系。Cokriging將高程作為協(xié)變量進(jìn)行插值，一定程度上避免了高程對(duì)誤差的影響，因此Cokriging方法的誤差與高程相關(guān)性較小。TPS方法的插值誤差與高程相關(guān)性較小，可基本認(rèn)為TPS方法可最大程度減少由高程帶來(lái)的溫度插值誤差。

圖4 不同插值方法對(duì)小梧桐站氣溫插值結(jié)果逐月分析 Fig. 4 Monthly analysis of temperature interpolation results at Xiaowutong station by different interpolation methods

圖5 5種插值方法交叉驗(yàn)證結(jié)果誤差絕對(duì)值與高程相關(guān)性分析 Fig. 5 Analysis of absolute error and elevation correlation of five interpolation methods for cross-validation results

3 結(jié)論

基于2017年深圳市128個(gè)氣象站點(diǎn)的氣溫監(jiān)測(cè)數(shù)據(jù)與深圳市數(shù)字高程模型，對(duì)比了反距離權(quán)重法、局部多項(xiàng)式插值法、普通克里金法、協(xié)同克里金法和Anusplin軟件中的薄板樣條函數(shù)法，對(duì)深圳市2017年全市月平均氣溫進(jìn)行了空間插值，通過(guò)氣象站點(diǎn)逐個(gè)交叉驗(yàn)證的方法，對(duì)比了不同空間插值方法在不同月份的插值結(jié)果。研究顯示：

1）與IDW、LPI、OK法相比，Cokriging法與基于Anusplin軟件的TPS法插值得到的擬合曲面均較為平滑。

2）在插值交叉驗(yàn)證結(jié)果的逐月對(duì)比中，與IDW、LPI、OK、Cokriging法相比，TPS法的均方根誤差與平均絕對(duì)誤差較小，可有效提高氣溫?cái)?shù)據(jù)的空間插值的精度。2017年深圳市氣溫空間插值的交叉驗(yàn)證誤差顯現(xiàn)明顯的季節(jié)性特征，TPS法的春夏兩季精度明顯高于秋冬兩季。

3）高程是IDW、LPI、OK法插值重要的誤差來(lái)源，其他方法插值與協(xié)變量為高程的TPS法的主要差別在于不考慮高程所帶來(lái)的系統(tǒng)誤差。TPS法在氣溫插值上的應(yīng)用充分考慮到了海拔的影響，因此避免了高程帶來(lái)的誤差，Anusplin軟件在氣溫空間插值領(lǐng)域值得推薦。