基于ANUSPLIN的黃土丘陵區氣象要素插值方法對比與評估

肖 旭， 鄭 誠， 丁成琴， 范陳哲， 白悅江， 林隆超，閆 婷， 高 宇， 史海靜,5

（1.陜西省延安市氣象局，陜西 延安 716000；2.西北農林科技大學草業與草原學院，陜西 楊凌 712100；3.西北農林科技大學水土保持研究所，陜西 楊凌 712100；4.西北農林科技大學資源環境學院，陜西 楊凌 712100；5.中國科學院水利部水土保持研究所，陜西 楊凌 712100）

氣象因子是影響農業生產、人類活動以及生態環境的重要要素，尤其是在干旱半干旱地區。氣象要素在農業和環境等相關領域研究較多，主要應用于模型預測和作物生產等方面［1］。關于氣象的研究，最初人們依賴于氣象站的雨量計監測，但由于氣象站點數目較少，空間不完全覆蓋，且時效性較差，在完整捕捉氣候事件和描述空間變化等方面能力不足。隨著空間處理和分析技術的發展，空間插值的方法被用來描述氣象因子的空間分布，一定程度上實現了氣象空間覆蓋，提高了人們對氣象空間變化的認識。

目前，插值方法較多，且在系統模擬和反演領域都有不同程度的應用，常見的方法有局部多項式法(Local polynomial interpolation)、樣 條 函 數 法(Spline)、克里金插值法(Kriging)和反向距離權重法(Inverse Distance Weight, 簡稱IDW)等［2-4］。ANUSPLIN 是由Hutchinson［5］開發用于不同空間和時間尺度水文氣象插值的方法，其主要特點是可包含線性協變量例如地形變量等，廣泛應用于高山丘陵等地形變化起伏較大的地區［6-7］。此外，ANUSPLIN 具有簡單的操作流程，不需要對其參數進行單獨的預先校準，生成的溫度降水時間序列曲面具有較高的精度和可靠性［8］。近年來，ANUSPLIN 插值方法已被國內外學者廣泛應用于氣候數據的插值［9-10］。黃土丘陵區位于黃土高原中部，是我國生態恢復的重點區域之一，該區溝壑縱橫，地形起伏多變，溫度降水等氣象因子受地形影響較大［11］，傳統的插值方法多通過氣象站空間距離進行估算，而黃土丘陵區多以山地為主，地形變化足以影響空間溫度和降水的差異性分布［12］，單從距離的遠近獲得氣象數據并不能真實的反映該地區的氣象變化，從而無法滿足植被恢復和重建工作的需求。因此，針對地形復雜的丘陵區，探索包含不同分辨率地形變化的氣象插值方法對于進一步完善氣象空間數據和相關研究具有重要意義。

本文采用ANUSPLIN 空間插值方法，對延河流域2010—2020年的溫度降水進行空間插值，同時選取25 m、90 m、1 km 三個不同精度的數字高程模型(Digital Elevation Model，DEM)作為協變量，對比分析各個尺度溫度降水插值結果的精度，評估ANUSPLIN 插值方法在黃土丘陵溝壑區的適用性，并選擇區域溫度降水最優數據集，補充丘陵區氣象插值精度不夠的短板，為準確獲取地形起伏較大的丘陵區氣象空間插值數據提供可參考的理論依據和方法支持。

1 數據與方法

1.1 研究區概況

本研究以典型黃土丘陵區的延河流域為研究區域(36°23′～37°17′N，108°45′～110°28′E)，該區覆蓋面積為7687 km2，海拔454～1765 m，該流域從東南向西北依次經過延長、延安、安塞和志丹，是典型的丘陵區（圖1）。氣候全年較干，夏季多雨，冬季降水較少，處于暖溫帶大陸性季風半濕潤氣候向溫帶半干旱氣候的過渡區，全年降雨量496 mm，年平均氣溫8.8 ℃。從東南向西北，降雨、溫度呈現出明顯的梯度變化特征。土壤以黃綿土為主, 植被有森林(刺槐和遼東櫟)、灌木叢和草地等。

1.2 數據來源與處理

1.2.1 氣象數據 基礎氣象數據來自延河流域及周邊站點。數據獲取時間為2010—2021年，數據類型為日值數據，包含日平均溫度、最高溫度、最低溫度、日降水量等氣候要素，本研究通過多年（11 a）的年均溫度和年均降水進行插值。由于是長時間序列數據，部分數據存在異?；蛉笔У?，因此，數據經過篩選、剔除異常數據，最終保留105個氣象站（圖1）。

1.2.2 DEM 數據 延河流域25 m 分辨率DEM 是通過對日本ALOS 衛星搭載的PALSAR 傳感器25 m DEM 數據進行處理得到的。ALOS 25 m DEM 可從NASA 的ASF DAAC (https://search.asf.alaska.edu/#/)免費下載。90 m 和1 km 分辨率DEM 數據來自CGIAR-CSI SRTM中國區域數據。

1.3 研究方法

1.3.1 ANUSPLIN 模型基本原理 ANUSPLIN 是一種使用薄板平滑樣條對噪聲多變量數據進行透明分析和插值的工具。薄板平滑樣條函數法實際上是利用一個平面去擬合所有的氣象站點，使得該曲面能通過站點構成“樣條”，得到逼近所有控制點彎曲最小的光滑曲面，即利用M 次多項式，對曲線進行分段擬合，用局部的薄板來繪制分段連續的曲線。TPS 是對樣條函數法的曲面擴展，插值過程中利用最優的光滑參數實現逼真度和光滑度最佳平衡，保證了精度可靠，同時插值曲面光滑連續［13］。此外，ANUSPLIN 還允許轉換自變量和因變量，并允許處理具有缺失數據值的數據集。 當對因變量應用變換時，ANUSPLIN 允許擬合曲面的反向變換，計算相應的標準誤差，并糾正這些變換引起的小偏差。研究發現，這在將表面擬合到降水數據和其他自然為正或非負的數據時特別方便［14］。

薄板樣條估計z（xi）由一個適當光滑的函數計算，該函數最小為:

式中：yi為觀測i的觀測數據值；?(xi)為被估計的樣條函數；λ為一個正數，稱為平滑參數；Jm（?）是用m階偏導數定義的函數?的粗糙度的度量［m階偏導數稱為SPLINA 輸入中的樣條的階(它也被稱為粗糙度階）］［10］。

局部薄盤光滑樣條（Partial thin plate smoothing splines）的理論統計模型公式：

式中：Zi為位于空間i點的因變量；T為迭代次數；?(xi)為估算關于xi的未知光滑函數，xi為獨立變量；yi為p維獨立協變量；b為yi的p維系數；ei為隨機誤差［10］。

上式中，函數?和系數b是通過最小二乘估計來確定的：

式中：wi為權重的已知局部相對變異系數；Jm(?)為函數?(xi)的粗糙度測度函數，定義為?的偏導（稱為樣條次數，也叫粗糙次數）；λ為正的光滑參數，在數據保真度與曲面的粗糙度之間起平衡作用。在AUNSPLIN中通常用廣義交叉驗證GCV的最小化以及最大似然法GML的最小化來確定［15-16］。

1.3.2 插值精度評價指標 本研究選取數據的15%作為驗證集，選取過程采用隨機抽樣，將隨機站點數據與插值結果進行比較，以評估氣象模擬的準確性和偏差［17］。三種分辨率的插值精度比較適合使用交叉驗證法。本研究插值的效果采用均方根誤差RMSE來評估，表達式如下：

式中：Zx和Zy分別代表氣候數據實測值和插值結果；n為計算序列長度。RMSE 的值越小表示氣象模擬效果越好。

此外，將實測值與預測值進行相關性分析，用確定性系數R2表示一元多項式回歸方程擬合度的高低。R2取值分布在0～1之間，R2值越大，表示模擬效果越好。

2 結果與分析

2.1 插值結果

基于ANUSPLIN 分別獲取了研究區2010—2021 年年均溫和累計降水的插值曲面見圖2，從圖中可以看出，相較于1 km 的插值，90 m 與25 m 分辨率的插值結果細節表現得更加突出，能反映局部地形特征，也有較好的平滑度。

由圖2 還可以看出，東部延長縣逐日年均溫度較高，西部靖邊、志丹以及中部的安塞區逐日年均溫度較低，三種分辨率的逐日年均溫度整體呈“東高西低”的趨勢。25 m 分辨率的溫度插值范圍為7.725～13.329 ℃（圖2a），90 m 分 辨 率 為7.786～13.354 ℃（圖2b），1 km分辨率為7.872～13.313 ℃（圖2c）。

從降雨插值可以看出，低降雨的空間分布主要集中在延河流域西部的靖邊地區、中部的寶塔區以及東部延長的少部分區域。25 m 分辨率的逐日降雨插值范圍為1.062～1.847 mm，90 m 分辨率的插值范圍為1.046～1.852 mm，兩者最低降雨量相差0.016 mm，上限相差0.05 mm，較為接近，但1 km 分辨率的插值結果與25 m、90 m 的插值結果相差過大，下限為1.117 mm，比90 m 的結果高0.071 mm、比25 m 的插值結果高0.055 mm；上限為1.73 mm，比兩者都低，與25 m 的插值上限相差0.117 mm，整體插值結果范圍縮小。

從插值的溫度降雨空間圖中可以看出，細節特征比較突出，且有明顯的“馬賽克”紋路，與過去的Cokriging 插值相比，沒有出現較大的斑塊與“牛眼”，說明ANUSPLIN 有比較好的平滑度，且很好的擬合了地形特征，這與劉林等［18］的研究結論一致。此外，插值結果顯示延河流域的中部、西北部降雨較低，東部降雨量較高，這在劉春利等［19-20］的研究結果中得以驗證。且土壤水分變化也與之一致［21］。溫度插值特征也符合以往氣象站公布的數據規律。延河流域的西北部為草原地帶［22］，整體呈降水低，溫度低，東部為森林草原地帶與森林帶，降水高，溫度高，與插值結果一致，說明ANUSPLIN 模型能較好的適應黃土高原丘陵溝壑區。

2.2 插值結果比較

本研究以25 m、90 m、1 km 三種分辨率的海拔數據為協變量對黃土丘陵溝壑區的溫度降雨進行插值，相較于前期研究中采用的1 km×1 km 或500 m×500 m 分辨率［23］，本研究在一定程度上提高了精度，減小了插值誤差。從圖3 可以得出，對于不同DEM分辨率下溫度插值結果，根據計算的均方根誤差（RMSE），25 m 分辨率的插值數據精度較好，90 m 分辨率插值數據精度其次，1 km插值數據精度較差。而在確定性系數（R2）方面，25 m 分辨率的R2為0.65，擬合效果較好，1 km分辨率擬合效果其次，90 m分辨率擬合效果稍差。綜上所述，在三種不同精度的DEM分辨率模擬情景下，溫度插值精度排序為：25 m＞90 m＞1 km。賈洋等［24］的研究表明，DEM 分辨率越高，ANUSPLIN 模型插值精度也越高，溫度插值結果也遵循這一規律，高精度的地形數據可以有效反映空間差異，也可提高溫度空間插值的精確度。

圖3 基于不同分辨率DEM溫度降水觀測值與預測值分析結果Fig.3 Analysis results of observed and predicted values of temperature rainfall based on DEM with different resolutions

降雨插值結果精度檢驗按照年均溫來檢驗，空間數據取365 d。對于降雨插值結果的RMSE，90 m<25 m<1 km，即90 m 分辨率的插值精度最高，1 km 的插值精度較低。同時，90 m 分辨率的插值結果與實測結果擬合時，90 m 分辨率也是三種分辨率中最大的，R2=0.833，擬合效果良好，其次是25 m 分辨率，1 km 分辨率的擬合效果最差。這與溫度插值規律并不相同，筆者推測出現這一結果的原因可能是由于黃土丘陵區山地較多［25-27］，降水容易在溝地匯集，導致站點數據存在誤差。尤其是黃土丘陵區普遍坡長在70～200 m之間，90 m的空間精度與地形變化精度恰好一致，確定性檢驗仍需后續工作進一步驗證。但是，該結果也說明在丘陵區設置氣象站點要考慮地形變化。

高分辨率數據集(包括遙感數據和長時間序列的氣候數據)是目前黃土丘陵區相關研究迫切需要的。且近年來干旱區氣候變化較大［28］，實時氣象數據對于氣候-植被研究和人類活動具有重要價值。這一氣象數據和空間分布研究對于干旱區氣候變化研究非常寶貴，同時高精度分辨率的數據集可以為該區域生態恢復和氣候變化研究提供基礎數據支撐［29］。

3 結論與討論

本研究選取黃土高原典型丘陵區延河流域作為研究區，以研究區內及周邊105 個氣象站2010—2021 年逐日氣象數據為基礎資料，以1 km、90 m、25 m 三種分辨率的數字高程數據為協變量，利用ANUSPLIN 插值方法對研究區不同地形分辨率氣象數據進行空間插值，通過氣象站站點交叉驗證方法，對比分析各分辨率數字高程數據的插值精度水平。研究發現數字高程數據的分辨率影響氣象空間插值的結果，具體結果如下：

（1）溫度插值顯示，延河流域東部延長地區溫度較高、西部溫度較低符合以往氣象站數據規律［30］；降水量空間插值顯示中部與西北部較低，與該區域蒸散發時空變化趨勢基本特征一致［31］，延安市周邊本研究降雨量較少，其主要原因可能是人為干擾導致氣象站點數據存在誤差，以往研究也表明，延河流域的延安市是人類活動較為頻繁的區域［32］。在趙美亮等［33］在西北地區水文空間分布研究中表明，土地利用改變會影響水文要素分布?？偟膩碚f，ANUSPLIN 模型對黃土丘陵溝壑區溫度降雨空間插值有較好的適應性。

（2）在三種不同的DEM 分辨率模擬場景下，溫度插值精度排序為：25 m＞90 m＞1 km，符合數字高程數據分辨率越高，插值精度越高［34］；降雨插值精度排序為：90 m>25 m>1 km,不完全符合以上結論，表明降雨數據空間插值存在最佳分辨率的數字高程數據，最佳數字高程數據與該區域地形特征有關。在野外觀測中發現，黃土丘陵區山脈坡長多60～300 m之間，與該結果具有一致性，坡度和坡長等地形數據可能更適合描述降雨數據信息，而不是25 m 的海拔數據。

本研究基于ANUSPLIN 插值方法，結合數字高程數據進行多種空間尺度的空間精度驗證，該結果與觀測結果擬合較好，說明了結果具有可靠性，所得結論能夠為丘陵區插值的應用提供參考。然而，本研究僅考慮了海拔數據作為協變量，然而坡向、坡位等地形因子對插值結果影響也不可忽視，將多種協變量考慮在ANUSPLIN 插值方法中，可以為氣象插值精度優化方法提供更有意義的參考依據。