劉煥莉，范增祿，韓明稚，田國強

(河北省氣象信息中心，河北 石家莊 050021)

引言

空間插值是基于已知采樣點的數(shù)據(jù)預(yù)測未知空間(點)數(shù)據(jù)值的過程[1]。目前，常用的空間插值方法有反距離加權(quán)法、趨勢面法、克里金插值法、薄盤樣條插值法等, 以及形式多樣的多元線性回歸方法等[2]，上述插值方法在本質(zhì)上均基于一定的統(tǒng)計學(xué)模型，但在空間化模擬過程中各有不足，其差別主要表現(xiàn)在適用性、插值原理、空間分辨率等方面[3]。因此，通過比較而選擇一個合適的、適用于數(shù)據(jù)空間分布特點的插值方法是空間插值的關(guān)鍵[4]。

20世紀(jì)90年代以來，以ANUSPLIN專用空間插值軟件為代表的區(qū)域精細(xì)化插值方法在國內(nèi)外得到廣泛應(yīng)用[5-14]。如錢永蘭等[15]利用ANUSPLIN軟件對中國境內(nèi)667個基本、基準(zhǔn)地面氣象觀測站1961—2006年逐日氣溫和降水?dāng)?shù)據(jù)進行插值，并將插值結(jié)果分別與反向距離權(quán)重法和克里金法插值結(jié)果進行對比，結(jié)果證明ANUSPLIN插值法對降水的插值誤差較另外兩種方法略大，而氣溫的插值效果則優(yōu)于反向距離權(quán)重法和克里金法。劉志紅等[16]選擇黃土高原中北部多沙粗砂區(qū)為研究區(qū)，利用ANUSPLIN方法對研究區(qū)內(nèi)52個氣象站近21 a月平均多種氣象要素進行插值，結(jié)果證明在所有氣象要素中，溫度的空間插值誤差最小，風(fēng)速、水汽壓的誤差中等，日照時數(shù)和降雨量的誤差較大?？臻g內(nèi)插的精度與臺站分布密度息息相關(guān)，趙煜飛等[17]分析了不同臺站密集度情況下插值的精度大小，得出臺站密集程度越高，插值效果越好的結(jié)論。廖榮偉等[18]利用ANUSPLIN、SHERPAD和OI插值方法對中國地面氣溫數(shù)據(jù)進行了網(wǎng)格化并對插值精度進行了比較，得出了考慮高程的ANUSPLIN插值的綜合評價較好。

以上研究多是基于長時間序列的國家級地面氣象觀測站氣象要素進行插值，數(shù)據(jù)集空間分辨率多為0.5°×0.5°和0.1°×0.1°，高空間分辨率的格點數(shù)據(jù)具有規(guī)范的格式、易于分析、能更精細(xì)地表征空間分布特征等優(yōu)點[19]。鑒于溫度、降水等要素具有明顯的局地小氣候特征[20]，需要高空間密度的觀測站點來捕獲數(shù)據(jù)，本文選取京津冀區(qū)域以及臨近省共3 974個國家級及區(qū)域氣象觀測站，利用ANUSPLIN插值軟件對2018年汛期6—8月逐日氣溫資料進行插值，得到了京津冀區(qū)域逐日氣溫0.01°×0.01°格點數(shù)據(jù)集，并對數(shù)據(jù)集的質(zhì)量進行檢驗評估。

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

1.1 研究區(qū)域

研究區(qū)域選擇京津冀地區(qū)(圖1)，京津冀地區(qū)位于華北平原北部，行政區(qū)域包括北京、天津、雄安新區(qū)及河北，地域面積21.6萬km2，位于113.5°～120°E，36°～42.3°N。西倚太行山與山西省為鄰；北部與內(nèi)蒙古自治區(qū)接壤；南部與河南、山東兩省毗鄰；東與遼寧省相接。京津冀區(qū)域地貌類型多樣，平原、丘陵、盆地、山地、高原等地貌類型盡在其中，整體地勢呈現(xiàn)西北高東南低的地形特點。區(qū)域內(nèi)海拔2 000 m以上的山峰有10座，第一高峰位于張家口小五臺山東臺，海拔2 882 m[21-22]。

圖1 研究區(qū)示意圖(a.國家級和區(qū)域氣象觀測站分布；b.高程，填色，單位：m)

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

本文研究所使用的數(shù)據(jù)為113.4°～120°E，36°～42.5°N范圍內(nèi)2018年6—8月京津冀區(qū)域2 402個站點的逐日氣溫數(shù)據(jù)，為減弱邊界效應(yīng)造成的插值誤差，引入臨近省區(qū)(內(nèi)蒙古、山西、山東、河南、遼寧)1 572個站點數(shù)據(jù)。數(shù)據(jù)來源于河北省氣象信息共享平臺，數(shù)據(jù)內(nèi)容包括臺站經(jīng)緯度、海拔高度和日平均氣溫。數(shù)據(jù)的準(zhǔn)確性對插值結(jié)果影響顯著，因此原始數(shù)據(jù)均經(jīng)過格式檢查、缺測檢查、界限值檢查、內(nèi)部一致性檢查、時變檢查、持續(xù)性檢查等質(zhì)量控制。

由于氣溫的空間分布與地形特征緊密相連，采用ANUSPLIN軟件對研究區(qū)2018年6—8月逐日氣溫進行空間插值，為提高空間插值精度，減弱地形條件對氣溫空間插值精度的影響，真實反映研究區(qū)域氣溫的時空分布規(guī)律，引入高程數(shù)據(jù)DEM[23]，分辨率為1 km，通過數(shù)字高程模型提取地形特征因子(如海拔高度、坡向、坡度、地形遮蔽度等)，利用ArcGIS將DEM數(shù)據(jù)轉(zhuǎn)化成ASCII格式，作為協(xié)變量參與空間插值運算。

2 研究方法

2.1 空間插值方法

本文采用ANUSPLIN專用空間插值軟件對京津冀區(qū)域氣溫站點資料進行空間插值，為評估ANUSPLIN軟件的插值效果，以日平均氣溫為研究對象，綜合比較基于多種插值方法的插值精度，確定最適合京津冀地區(qū)氣溫的插值方法。

2.1.1 ANUSPLIN插值

ANUSPLIN軟件是基于普通薄盤和局部薄盤樣條函數(shù)對多變量數(shù)據(jù)進行內(nèi)插的工具，能同時進行多個表面的空間插值，適用于時間序列的氣象數(shù)據(jù)。ANUSPLIN軟件提供了完整的統(tǒng)計分析、數(shù)據(jù)診斷等功能，同時也支持多種數(shù)據(jù)輸入和表面查詢功能。

局部薄盤光滑樣條(partial thin plate smoothing splines)是對薄盤光滑樣條原型的擴展[24], 它除普通的樣條自變量外允許引入線性協(xié)變量子模型, 如溫度和海拔之間、降水和海岸線的關(guān)系等[25]。

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

zi=f(xi)+bTyi+ei(i=1,…,N)

(1)

式中：zi是位于空間i點的因變量；xi是獨立變量，f是要估算的關(guān)于xi的未知光滑函數(shù)；yi為獨立協(xié)變量；b為yi的系數(shù)；ei為隨機誤差。

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

(2)

其中Jm(f)是函數(shù)f(xi)的粗糙度測算函數(shù)，定義為f的m階偏導(dǎo)，ρ為光滑參數(shù)，在數(shù)據(jù)保真度與曲面的粗糙度之間起平衡作用，在ANUSPLIN中通常用廣義交叉驗證GCV的最小化以及最大似然法GML的最小化來確定。

ANUSPLIN程序共分為8個程序模塊(SPLINA，SPLINB，SELNOT，ADDNOT，DELNOT，GCVGML，LAPPNT，LAPGRD)。通過輸入固定格式的數(shù)據(jù)文件，包含站號、經(jīng)度、緯度、高程信息(DEM數(shù)據(jù))，再執(zhí)行ANUSPLIN中的內(nèi)置模塊函數(shù)，最終輸出插值后的數(shù)據(jù)結(jié)果，同時包含光滑參數(shù)RHO、擬合數(shù)據(jù)誤差列表文件(均值、標(biāo)準(zhǔn)差、方差等)、擬合表面系數(shù)的協(xié)方差矩陣、誤差表面等。薄盤光滑樣條總共有18種插值模型選擇 (獨立變量、協(xié)變量和樣條次數(shù)多種組合) 見表1，根據(jù)不同氣象要素選擇合適插值模型。

ANUSPLIN最佳模型判斷標(biāo)準(zhǔn)：GCV或GML最??；信噪比最小；信號自由度小于站點的一半，但不能過小以至于尋找不到最優(yōu)光滑參數(shù)；模型成功率判斷中無“*”。

表1 薄盤光滑樣條函數(shù)模型

2.1.2 反距離權(quán)重方法

反距離權(quán)重插值基于相似相近原理，即兩個物體離得越近，它們的值越相似，反之，則相似性越小。反距離權(quán)重(inverse distance weighted，IDW)方法以插值點和樣點的距離為權(quán)重進行加權(quán)平均[26]，算法如下：

(3)

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

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

2.1.3 普通克里金法

普通克里金(ordinary kriging,OK)法是從地統(tǒng)計學(xué)的克里格法中演化出的一種插值方法，它以區(qū)域性變量理論為基礎(chǔ)，半變異函數(shù)為分析工具，其實質(zhì)是利用半變異函數(shù)對未知樣點進行線性無偏、最優(yōu)估計。它不僅考慮預(yù)測點位置與已知數(shù)據(jù)位置的相互關(guān)系，而且還考慮變量的空間相關(guān)性[8]。算法如下：

(4)

(5)

式中：z(x0)為x0點的預(yù)測值，z(xi)為xi點的測量值，λi為測量值對預(yù)測值的權(quán)重系數(shù)。

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

2.1.4 樣條函數(shù)法

樣條函數(shù)(spline function， SPLINE)法主要通過估計方差，利用一些特征節(jié)點，用多項式擬合的方法來產(chǎn)生平滑的插值曲線[27]。算法如下：

(6)

樣條函數(shù)插值法具有易操作，插值速度快的優(yōu)點，但存在插值后誤差不能直接估計，采樣點稀少時效果不好的缺點。

2.2 評估指標(biāo)

為評估ANUSPLIN的插值效果，選取空間分布均勻的15%站點作為檢驗站點(檢驗站點不參與插值)，剩余站點作為插值站點，檢驗站點分布見圖2。分別使用IDW、OK、SPLINE和ANUSPLIN進行插值，使用最近距離法將插值的高分辨率資料匹配到京津冀區(qū)域各站點上，得到了逐日氣溫空間插值數(shù)據(jù)集及實測數(shù)據(jù)集，采用相關(guān)系數(shù)(Corr)、平均絕對誤差(mean absolute error，MAE)、平均相對誤差(mean relative error，MRE)等作為評估插值效果的指標(biāo)，討論不同方法對氣溫插值結(jié)果的影響。

圖2 插值站點與檢驗站點分布

3 結(jié)果與分析

3.1 插值結(jié)果與分析

依據(jù)ANUSPLIN最佳插值模型判斷標(biāo)準(zhǔn)，篩選出模型7、模型8、模型9為最優(yōu)待用模型，再分別用三個待用模型進行插值，分析其誤差分布，確定氣溫的最優(yōu)插值模型為模型8。

圖3為采用4種空間插值方法構(gòu)建的6、7、8月每月1日的京津冀區(qū)域氣溫空間分布圖，通過比較可以發(fā)現(xiàn)，ANUSPLIN體現(xiàn)了引入DEM預(yù)測氣溫的優(yōu)勢，插值結(jié)果光滑，清晰反映出氣溫過渡變化的特點，而IDW、OK、SPLINE在曲線平滑程度以及局部地區(qū)的空間分布具有一定差異，其中IDW、SPLINE插值結(jié)果容易產(chǎn)生“牛眼”現(xiàn)象。因此，在選取插值方法時，應(yīng)充分考慮空間特異性，引入適當(dāng)?shù)挠绊懸蜃樱浞址治鰯?shù)據(jù)，才能對特定環(huán)境下的氣象要素進行精確、合理的插值[28]。

圖3 日平均氣溫空間分布(a. ANUSPLIN, b. IDW, c. OK, d. SPLINE；單位：℃)

綜上所述，4種插值方法獲得的格點化氣溫數(shù)據(jù)均能直觀體現(xiàn)京津冀區(qū)域氣溫由北向南遞增的空間分布特征，氣溫分界線沿太行山脈、燕山山脈，分界線以北區(qū)域溫度較以南平原地區(qū)氣溫低，日氣溫最大值均出現(xiàn)在低海拔平原地區(qū)，而最小值均出現(xiàn)在崇禮圍場高海拔地區(qū)，日氣溫整體呈現(xiàn)東南高西北低的趨勢。究其原因是在影響氣溫分布的諸多因素中，以海拔高度和地形影響最顯著[29-30]。

3.2 評估結(jié)果與分析

相關(guān)系數(shù)(Corr)可以用來表征預(yù)測值與實測值之間的相似程度，用誤差占比、MAE、MRE作為評估指標(biāo)來檢驗插值精度。相關(guān)系數(shù)檢驗結(jié)果如下：由圖4可以看出，4種插值方法得出的預(yù)測值與實測值相關(guān)性均較好，相關(guān)系數(shù)均大于0.86，ANUSPLIN相關(guān)系數(shù)高于其他3種插值法，相關(guān)系數(shù)平均達0.97，IDW、OK次之，相關(guān)系數(shù)平均達0.96，而SPLINE略低，相關(guān)系數(shù)平均達0.95。

圖4 基于不同插值方法的相關(guān)系數(shù)分布(紅色線為ANUSPLIN，藍色線為OK，黃色線為IDW，綠色線為SPLINE)

基于IDW、OK、SPLINE和ANUSPLIN軟件插值的結(jié)果有一定差距，通過計算預(yù)測值與實測值的誤差占比、MAE、MRE作為檢驗插值精度的標(biāo)準(zhǔn)，可以看出，不同插值方法的MAE、MRE、誤差占比并不相同(圖5、圖6、表2)。由表2可見，4種插值方法中，基于ANUSPLIN軟件的插值結(jié)果最優(yōu)，其次為IDW、OK，基于SPLINE的氣溫插值結(jié)果最差。

圖5 平均絕對誤差空間分布(a. ANUSPLIN, b. IDW, c. OK, d. SPLINE；單位：℃)

圖6 平均相對誤差空間分布(a. ANUSPLIN, b. IDW, c. OK, d. SPLINE；單位：%)

表2 誤差參數(shù)分析

在氣溫插值誤差空間分布方面，插值誤差較大的站點分布在冀北高原燕山丘陵及太行山脈高海拔地區(qū)，而平原地區(qū)則誤差較小，其原因可能是冀北高原燕山丘陵及太行山脈一帶站點稀少且海拔較高，賈洋和崔鵬[11]的研究也證明ANUSPLIN插值精度主要受海拔影響，同時站點密度也是造成插值誤差的原因之一[31]。

綜上所述，4種插值方法在京津冀區(qū)域插值精度為ANUSPLIN>IDW>OK>SPLINE,基于ANUSPLIN軟件的插值方法更加適合京津冀區(qū)域的逐日氣溫空間插值。

4 結(jié)論與討論

本文使用京津冀區(qū)域以及臨近省區(qū)2018年國家級及區(qū)域氣象觀測站日氣溫數(shù)據(jù)，采用ANUSPLIN研制京津冀區(qū)域逐日氣溫格點數(shù)據(jù)集，并與IDW、OK、SPLINE插值方法進行對比分析，得出以下結(jié)論：

1)基于4種插值方法對京津冀區(qū)域逐日氣溫數(shù)據(jù)進行插值的結(jié)果表明，ANUSPLIN插值結(jié)果與數(shù)字高程模型極為相似，插值效果優(yōu)于IDW、OK、SPLINE，精度更高、光滑度更好。

2)通過計算預(yù)測值與實測值的相關(guān)系數(shù)(Corr)、MAE、MRE等作為檢驗插值精度的標(biāo)準(zhǔn)，4種插值方法得出的預(yù)測值與實測值相關(guān)性均較好，在站點密度足夠高的前提下，各種插值方法用于同一組數(shù)據(jù)的插值結(jié)果差異不大，而插值誤差則存在差異，綜合評估，ANUSPLIN的插值結(jié)果最優(yōu)，其次為IDW，再次為OK，基于SPLINE的氣溫插值結(jié)果最差，ANUSPLIN最適用于京津冀區(qū)域逐日氣溫插值。

3)京津冀區(qū)域逐日氣溫插值誤差較大的站點分布在冀北高原、燕山丘陵及太行山脈一帶，平原地區(qū)誤差較小，海拔及站點密度是影響ANUSPLIN插值精度的因素之一。

綜合分析，引入DEM的ANUSPLIN插值效果最優(yōu)，但由于受模型本身的局限，本文只選擇了經(jīng)度、緯度和高程(m)三項作為自變量，沒有考慮更多影響局地氣溫變化的因素，建議在后續(xù)的研究中，針對不同長度的時間序列、海拔、地形、DEM誤差等因素開展更深入的研究。同時在驗證插值精度時，僅采用了相關(guān)系數(shù)(Corr)、MAE、MRE幾個誤差指標(biāo)進行簡單驗證，考慮到插值過程中會平滑掉一些極端值，在后續(xù)的研究中，應(yīng)盡可能選取更多的誤差指標(biāo)及氣溫特征值進行驗證，全面衡量插值精度。