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

衛(wèi)倩倩，陳霄健，閆佳逸，張 寧

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

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

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

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

1 研究區(qū)域和數(shù)據(jù)來源

1.1 研究區(qū)域

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

圖1 山西省地貌圖

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

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

2 研究方法

2.1 空間插值方法

2.1.1 薄盤光滑樣條法

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

局部薄盤光滑樣條的理論統(tǒng)計模型為：

式（1）中：zi為位于空間i點(diǎn)的因變量；f為要估算的關(guān)于xi的未知光滑函數(shù)；xi為d維獨(dú)立變量；b為yi的p維系數(shù)；yi為p維獨(dú)立協(xié)變量；ei為隨機(jī)誤差[13]。

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

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

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

2.1.2 普通克里金法

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

式（3）（4）中：z（x0）為x0點(diǎn)的預(yù)測值；z（xi）為xi點(diǎn)的測量值；λi為測量值對預(yù)測值的權(quán)重系數(shù)。

普通克里金法的優(yōu)點(diǎn)是可估計測試參數(shù)的空間變異分布以及估計參數(shù)的方差分布，但計算步驟煩瑣，計算量大，且變異函數(shù)需要根據(jù)經(jīng)驗(yàn)人為選定[14]。

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

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

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

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

2.2 研究步驟

選用研究區(qū)內(nèi)245個國家級站點(diǎn)2019年的日平均氣溫數(shù)據(jù)作插值運(yùn)算，為評估插值結(jié)果的準(zhǔn)確性，同時對省級氣象觀測站數(shù)據(jù)的準(zhǔn)確性進(jìn)行校驗(yàn)，選取765個省級氣象觀測站作為校驗(yàn)站點(diǎn)。在進(jìn)行插值前，對數(shù)據(jù)的完整性再次進(jìn)行檢驗(yàn)，剔除其中日平均氣溫缺測的站點(diǎn)。

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

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

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

根據(jù)經(jīng)緯度信息，依次將TPS、IDW、OK這3種方法得到的插值表面數(shù)據(jù)提取到相應(yīng)的省級氣象觀測站點(diǎn)，選用相對誤差、絕對相對誤差作為評估指標(biāo)，對3種插值方法的結(jié)果進(jìn)行評估。

日平均氣溫屬于計算值，對于誤差較大的站點(diǎn)，再結(jié)合氣象資料業(yè)務(wù)系統(tǒng)（MDOS）進(jìn)行分析，查看是否有因小時值缺測或質(zhì)控碼標(biāo)識為錯誤導(dǎo)致的日平均氣溫計算疑誤，對省級氣象觀測站的數(shù)據(jù)的準(zhǔn)確性進(jìn)行二次驗(yàn)證，同時也是對插值結(jié)果準(zhǔn)確性的驗(yàn)證。

3 結(jié)果與分析

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

3.1 插值表面對比分析

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

圖2 日平均氣溫插值表面

3.2 TPS法預(yù)測誤差表面分析

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

圖3 TPS插值的預(yù)測誤差表面

3.3 評估指標(biāo)分析

圖4和圖5分別為用TPS、IDW、OK這3種插值方法得到的01-15和07-15的絕對誤差分布圖。經(jīng)統(tǒng)計分析，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方法插值得到的結(jié)果明顯優(yōu)于用IDW和OK方法的插值結(jié)果，從圖6和圖7中，可以得到相同的結(jié)論。

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

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

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

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

3.4 MDOS入庫數(shù)據(jù)比對分析

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

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

4 結(jié)論與討論

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

選用絕對誤差、相對誤差2個評價指標(biāo)，對比分析3種插值方法的插值結(jié)果精度，發(fā)現(xiàn)用TPS方法插值得到的結(jié)果明顯優(yōu)于用IDW和OK方法得到的插值結(jié)果。

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

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