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

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

（1 中國科學院深圳先進技術研究院，深圳 518055；2 深圳市氣象局，深圳518040；3 香港大學，香港）

0 引言

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

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

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

1 數據與方法

1.1 數據來源

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

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

1.2 插值方法

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

1.2.1 反距離權重法

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

1.2.2 局部多項式插值法

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

1.2.3 普通克里金

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

1.2.4 協同克里金

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

1.2.5 薄板樣條函數

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

1.3 檢驗方法

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

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

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

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

2 結果與分析

2.1 空間分布對比

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

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

2.2 逐月對比

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

圖2 基于5種插值方法的深圳市7月平均氣溫空間分布 （a）OK法插值結果；（b）LPI法插值結果；（c）IDW法插值結果；（d）Cokriging法插值結果；（e）TPS法插值結果 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種插值方法交叉驗證誤差的RMSE、MAE逐月分析 Fig. 3 Monthly analysis of RMSE and MAE based on cross-validation errors of five interpolation methods

在插值精度上，RMSE與MAE的變化具有季節特征，在3月左右的插值精度明顯優于10月左右的插值精度。RMSE結果顯示，TPS插值法在春夏季節的插值精度優于秋冬季節。從誤差范圍來看，冬季的插值結果在負誤差方向表現比其他時間尺度明顯，表明利用Anusplin對研究區冬季氣溫數據插值時，會更大概率出現低估的情況[22-23]。某些地區由于局地氣候的影響，在季節氣溫方面表現出特殊性，如在冬季這些地區的站點氣溫遠高于同海拔區域[24]。本研究也發現，9月—次年2月的某些站點氣溫高于同海拔地區，導致Anusplin估算氣溫時出現低估的現象。這可能是夏季氣溫插值誤差小于冬季的原因。此外，不同站點下墊面性質的不同，對氣溫插值的結果也有一定影響，會造成一定的誤差。

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

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

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

2.3 影響因素分析

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

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

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

3 結論

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

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

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

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