山西省逐日日平均氣溫插值方法適用性評估＊

衛倩倩，陳霄健，閆佳逸，張 寧

（山西省氣象信息中心，山西 太原 030006）

在全球變化研究中，地學模型和氣候模型等相關研究的重要參數之一，是以高分辨率、格點化的氣候數據作為環境因子[1]。而氣溫是氣候變化研究中最基本指標之一，它不但在氣候變化檢測研究中有著基礎地位，同時又與氣候預測、氣候模式以及其他一系列的氣候變化問題研究有十分緊密的關系[2]。由于經緯度、海陸分布以及地貌特征與下墊面特性的差異等不同，體現在氣溫的空間分布上，會表現為具有明顯的區域特征[3]。通過地面氣象站觀測到的氣溫數據在空間上是離散的、孤立的點，只能代表觀測站點有限區域內的氣溫狀況，難以反映氣候空間變化連續過渡的基本特征，站網密度的增加也不能有效地將站點氣溫數據連接成面。因此，為獲得格點化的氣候數據，通常會選擇進行空間插值。

空間插值的實質是通過已知的離散樣點的數據，按照某種數學關系估算出未知點的數據。近年來，研究人員針對空間插值方法展開了大量的研究[4-6]，常用的空間插值方法主要有克里金法、反距離加權法、泰森多邊形法、多項式回歸法、樣條函數法等。張海平等[7]針對ArcGIS軟件中7種常用的空間插值方法，從不同視角進行適用性分類，比較各插值方法的差異，為研究人員如何選擇合適的插值模型提供參考。使用薄板局部光滑樣條插值方法的ANUSPLIN空間插值軟件，相較于傳統的插值方法，只將空間分布作為觀測數據的函數，同時考慮了相關因素的影響，插值結果也更符合氣象要素特點，插值結果的準確度也有所提高。

目前使用ANUSPLIN軟件進行氣候數據插值的應用較多，對較大時間尺度（月、年）或較小時間尺度（日），或不同區域面積上的不同氣象要素的插值研究均有所涉及[1，8-11]。本研究利用ANUSPLIN軟件，對山西省2019年國家級地面氣象站點的逐日日平均氣溫進行插值，并與采用普通克里金法、反距離權重法插值的結果進行對比，選取絕對誤差、相對誤差2個評估指標，分析比較3種插值方法的結果差異，找出更加適合山西省的逐日平均氣溫空間插值的方法，同時對省級氣象觀測站日平均氣溫的缺測數據的插補進行思考，提高其數據完整性。

1 研究區域和數據來源

1.1 研究區域

本文的研究區域為山西省，地貌如圖1所示。山西省地處黃河流域中部，因地屬太行山以西，故名山西。東與河北省為鄰，太行山作為天然屏障；西、南部以黃河為塹，與陜西省、河南省相望；北與內蒙古自治區毗連。地理位置范圍在北緯34°34′～40°44′，東經110°14′～114°33′之間，面積15.67萬km2。山西境內是典型的由廣泛黃土覆蓋的山地高原，地勢東北高西南低，地貌類型復雜多樣，有山地、丘陵、臺地、平原，山多川少，山地、丘陵面積占全省總面積的80.1%，大部分地區海拔在1 500 m以上，最高點為五臺山主峰葉斗峰，海拔3 061.1 m，為華北最高峰（信息來源于山西省人民政府網站省情概貌）。

圖1 山西省地貌圖

1.2 數據來源

本文研究所使用的數據為山西省2019年全年共874個站點的日平均氣溫資料，其中包括109個國家級地面氣象站（基準站、基本站、原一般站）和764個省級氣象觀測站（國家級業務考核站）。同時，為盡可能地減小由于地緣邊界站點減少造成的插值誤差[12]，選取周邊省份的136個國家站數據同時參與插值運算，經緯度范圍為北緯34°～41°，東經110°～115°，共計245個國家級站點數據。以上數據均來源于山西省氣象信息共享平臺，數據內容包括站點的經度、緯度、海拔高度和日平均氣溫。數據的準確性對插值結果影響較大，因此調取的數據均已通過氣象資料業務系統（MDOS）質量控制檢查。由于氣溫的分布與地形、海拔高度等密切相關，在插值過程中引入高程數據作為協變量，所需DEM數字高程數據下載自地理空間數據云網站（www.gscloud.cn），選取SRTMDEMUTM 90 M分辨率數據高程數據產品。

2 研究方法

2.1 空間插值方法

2.1.1 薄盤光滑樣條法

薄盤光滑樣條（Thin Plate Smoothing Spline，TPS）法是對樣條函數法的曲面擴展，常用于不規則分布數據的多變量平滑內插，利用光滑參數來達到數據保真度和擬合曲面光滑度之間的優化平衡[11]。局部薄盤光滑樣條是對薄盤光滑樣條原型的擴展[12]，它除普通的樣條自變量外，允許引入線性協變量子模型，如溫度和海拔之間、降水和海岸線的關系等。

局部薄盤光滑樣條的理論統計模型為：

式（1）中：zi為位于空間i點的因變量；f為要估算的關于xi的未知光滑函數；xi為d維獨立變量；b為yi的p維系數；yi為p維獨立協變量；ei為隨機誤差[13]。

函數f和系數b通過公式（2）的值最小來獲得，即由最小二乘法估計來確定：

式（2）中：ρ為光滑參數；Jm（f）為函數f（xi）的粗糙度測算函數，定義為f的m階偏導，通常由廣義交叉驗證GCV的最小化、最大似然法GML的最小化、或期望真實平方誤差MSE的最小化來確定。

ANUSPLIN軟件是基于普通薄盤和局部薄盤樣條函數對多變量數據進行內插的工具，能同時進行多個表面的空間插值，包括通過提供全面的統計分析、數據診斷等過程，同時也支持靈活的數據輸入和表面查詢功能。該軟件包含了18個插值模型，研究人員可以根據模型診斷結果，為不同氣象要素選擇最佳的插值模型。

2.1.2 普通克里金法

普通克里金（Ordinary Kriging，OK）法是從地統計學的克里格法中演化出的一種插值方法，它以區域性變量理論為基礎，半變異函數為分析工具，根據待估樣點有限鄰域內若干已測定的樣點數據，對未知樣點進行線性無偏、最優估計。它不僅考慮預測點位置與已知數據位置的相互關系，而且還考慮變量的空間相關性，是目前應用最廣的克里金插值方法[14-15]。其算法如下：

式（3）（4）中：z（x0）為x0點的預測值；z（xi）為xi點的測量值；λi為測量值對預測值的權重系數。

普通克里金法的優點是可估計測試參數的空間變異分布以及估計參數的方差分布，但計算步驟煩瑣，計算量大，且變異函數需要根據經驗人為選定[14]。

2.1.3 反距離權重法

反距離權重法（Inverse Distance Weighted，IDW）是一個均分的過程，通過對采樣點進行線性加權來決定輸出柵格值的大小，基于相似相近原理，加權與距離成反比[14-15]。該方法以插值點和樣點的距離為權重進行加權平均，算法如下：

式（5）中：z為估算值；n為參與插值的樣點數目；zi為第d i個樣點的觀測值；di為插值點到第i個樣點的距離；p為用于計算距離權重的冪指數，一般p=2。

IDW是根據點數據生成柵格圖層的常用方法，但它易受樣本點極值的影響，計算中常出現一種孤立點數據明顯高于周圍數據點的情況，即“牛眼”現象[14]。

2.2 研究步驟

選用研究區內245個國家級站點2019年的日平均氣溫數據作插值運算，為評估插值結果的準確性，同時對省級氣象觀測站數據的準確性進行校驗，選取765個省級氣象觀測站作為校驗站點。在進行插值前，對數據的完整性再次進行檢驗，剔除其中日平均氣溫缺測的站點。

利用ANUSPLIN軟件進行TPS插值。首先，將參與插值的國家級站點的日平均氣溫數據、DEM高程數值進行格式轉換，編寫軟件運行所需的cmd文件，運行splina模塊，得到相應的日志文件。對照最佳模型判斷標準，通過分析選擇，從ANUSPLIN軟件的18個待選插值模型中，確定模型9為最佳模型，即選用經度、緯度、高程（m）作為自變量，樣條次數為4，進行插值運算，得到輸出殘差文件、光滑參數文件、表面文件、列表文件、擬合表面系數的誤差協方差文件等。其次，根據得到的表面文件、擬合表面系數的誤差協方差文件，設定相應的輸出格式，編寫運行lapgrd模塊所需的cmd文件，得到最終的插值表面文件及插值預測誤差表面等相關文件。

利用ArcGIS軟件進行IDW和OK這2種插值。將參與插值的國家級站點的日平均氣溫數據，經格式轉換，保存為shapefile格式文件，再依次選擇IDW和OK模塊工具，進行相應的插值運算，得到插值表面文件。

分析對比3種插值方法得到的插值表面文件，同時對使用TPS方法得到的插值預測誤差表面文件進行分析。

根據經緯度信息，依次將TPS、IDW、OK這3種方法得到的插值表面數據提取到相應的省級氣象觀測站點，選用相對誤差、絕對相對誤差作為評估指標，對3種插值方法的結果進行評估。

日平均氣溫屬于計算值，對于誤差較大的站點，再結合氣象資料業務系統（MDOS）進行分析，查看是否有因小時值缺測或質控碼標識為錯誤導致的日平均氣溫計算疑誤，對省級氣象觀測站的數據的準確性進行二次驗證，同時也是對插值結果準確性的驗證。

3 結果與分析

選用典型月份的每月15日，即01-15、04-15、07-15、10-15這4 d的逐日平均氣溫數據進行插值結果分析。

3.1 插值表面對比分析

圖2中的（a）（b）（c）（d）依次為使用TPS、IDW、OK這3種插值方法得到的01-15、04-15、07-15、10-15這4 d的平均氣溫插值表面。從圖中可以明顯看出，使用TPS方法進行插值，得到的插值表面溫度分布與地形走勢更為一致，同時紋理清晰，分辨率更高、細節更突出，對于一些地勢起伏地帶的溫度變化也有很好地體現，可以更好地反映氣溫隨海拔高度的變化。使用IDW和OK這2種插值方法得到的插值表面，溫度變化趨勢不明顯，不能很好地體現山西境內地勢的高低不同對溫度產生的影響，其中，使用IDW方法進行插值時，“牛眼”現象尤為明顯，但對于小范圍的地勢起伏還是可以較好地反映出溫度隨海拔高度的變化。而使用OK插值方法得到的插值表面，對于小范圍的地勢起伏地帶，高低溫區別不明顯，得到的插值表面趨向于平滑的連續曲面。但三者都對海拔較高、溫度偏低的東北部一帶和海拔較低、溫度也相應偏高的中南部臨運盆地一帶的平川地區，這種大范圍的、地勢高低有明顯不同的高低溫區有很好的指示。

圖2 日平均氣溫插值表面

3.2 TPS法預測誤差表面分析

圖3中的（a）（b）（c）（d）依次為01-15、04-15、07-15、10-15用TPS插值方法得到的預測誤差表面圖。從圖中可以看出，由于引入了周邊省份的國家級站點數據進行插值，邊緣地帶的插值結果的準確性也得到了保證。插值誤差較大的區域主要分布在中部偏西一帶，對比站點分布圖及DEM數字高程圖，該范圍內海拔較高，且站點分布間隔相對較遠，同時參與插值的周邊省份站點也偏少，因此，預測誤差相對較高，這點也與前人研究結論相同[1]。

圖3 TPS插值的預測誤差表面

3.3 評估指標分析

圖4和圖5分別為用TPS、IDW、OK這3種插值方法得到的01-15和07-15的絕對誤差分布圖。經統計分析，01-15的插值絕對誤差小于等于2℃的比例依次為91.49%、78.40%、79.97%，平均絕對誤差依次為0.93℃、1.35℃、1.29℃；07-15的插值絕對誤差小于等于2℃的比例依次為92.37%、78.55%、70.66%，平均絕對誤差依次為1.06℃、1.49℃、1.75℃。由此可見，用TPS方法插值得到的結果明顯優于用IDW和OK方法的插值結果，從圖6和圖7中，可以得到相同的結論。

圖4 01-15絕對誤差分布圖

圖5 07-15絕對誤差分布圖

圖6 01-15相對誤差分布圖

圖7 07-15相對誤差分布圖

3.4 MDOS入庫數據比對分析

選用01-15利用TPS插值方法得到的插值數據與MDOS入庫小時數據進行分析判斷。

01-15有10個省級氣象觀測站日平均氣溫為缺測，依次查詢MDOS入庫的小時氣溫數據情況，其中有3個臺站僅有一個定時時次數據未正常入庫，導致日平均氣溫缺測。還有2個臺站缺測時次也小于6個時次，但根據省級氣象觀測站的日平均氣溫計算規則，只要有1個定時時次缺測，則日平均氣溫缺測，因此，這樣會造成一些臺站由于個別定時時次的缺測造成日平均氣溫的缺測，因此，這時可以考慮用插值數據來進行缺測補值，當然也要考慮對插值誤差的訂正。

4 結論與討論

利用局部薄盤光滑樣條函數法進行日平均氣溫空間插值，與反距離權重法和普通克里金法相比，由于考慮了溫度同海拔高度變化的關系，在插值過程中，引入高程作為協變量，插值結果更符合氣溫隨地勢高低的變化規律，精度更高、光滑度更好。但3種插值方法對于大范圍的、地勢高低有明顯不同的高低溫區的插值精度區別不大。

選用絕對誤差、相對誤差2個評價指標，對比分析3種插值方法的插值結果精度，發現用TPS方法插值得到的結果明顯優于用IDW和OK方法得到的插值結果。

省級氣象觀測站的日平均氣溫計算方法為4次定時值的平均，只要有一次定時缺測，則日平均氣溫相應缺測。由于傳輸不穩定或儀器故障等原因，引起的部分時次缺測現象比較常見，可以考慮利用插值數據進行缺測補值，提高氣象服務的數據完整性，對于如何進行誤差訂正，這部分還有待研究。

在下一步的研究當中，將以此研究為基礎，選用TPS方法，利用ANUSPLIN軟件完成山西省2019年逐日平均氣溫格點數據集的建立，并逐步擴展到長時間序列的高分辨率的其他氣象要素格點數據集的建立，并盡可能選取更多的誤差指標及各氣象要素相應的特征值進行驗證[6，16]，提高插值精度，為氣候研究提供參考。