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

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

0 引言

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

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

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

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

1.1 數(shù)據(jù)來源

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

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

1.2 插值方法

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

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

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

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

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

1.2.3 普通克里金

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

1.2.4 協(xié)同克里金

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

1.2.5 薄板樣條函數(shù)

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

1.3 檢驗(yàn)方法

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

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

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

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

2 結(jié)果與分析

2.1 空間分布對比

基于圖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é)果相近。對其他月份的插值結(jié)果均繪圖顯示，結(jié)果表明，Cokriging和TPS插值的結(jié)果更加平滑，而精確度有待進(jìn)一步評估，以下闡述五種插值方法的誤差檢驗(yàn)

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

2.2 逐月對比

之后交叉驗(yàn)證得到的MAE與RMSE作為判斷標(biāo)準(zhǔn)，判斷每個月份五種方法插值的精度高低。如圖3顯示，IDW法在10月時RMSE最大，為0.90 ℃，且MAE最大，為0.49 ℃。3月時RMSE與MAE最小，分別為0.67和0.36 ℃。LPI法的RMSE與MAE分別于10和9月達(dá)到峰值，為0.87和0.47 ℃。3月時LPI插值法的RMSE和MAE最小，為0.65和0.35 ℃。OK法的RMSE與MAE在10月達(dá)到最大，分別為0.90和0.49 ℃。3月時RMSE與MAE較小，分別為0.66和0.37 ℃。加入了高程作為協(xié)變量的Cokriging法的RMSE和MAE在10月達(dá)到最大，分別為0.83和0.56 ℃，3月時RMSE與MAE較小，為0.68和0.37 ℃。與IDW、LPI、OK、Cokriging法相比，TPS法的RMSE與MAE明顯較小，在9月達(dá)到峰值，分別為0.37和0.29 ℃，3月時RMSE與MAE較小，分別為0.25和0.18 ℃。在1—12月中每個月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é)。從誤差范圍來看，冬季的插值結(jié)果在負(fù)誤差方向表現(xiàn)比其他時間尺度明顯，表明利用Anusplin對研究區(qū)冬季氣溫數(shù)據(jù)插值時，會更大概率出現(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估算氣溫時出現(xiàn)低估的現(xiàn)象。這可能是夏季氣溫插值誤差小于冬季的原因。此外，不同站點(diǎn)下墊面性質(zhì)的不同，對氣溫插值的結(jié)果也有一定影響，會造成一定的誤差。

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

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

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

2.3 影響因素分析

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

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

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

3 結(jié)論

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

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

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

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