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

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

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

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

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

1 數(shù)據(jù)與方法

1.1 研究區(qū)概況

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

1.2 數(shù)據(jù)來源與處理

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

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

1.3 研究方法

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

薄板樣條估計z（xi）由一個適當(dāng)光滑的函數(shù)計算，該函數(shù)最小為:

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

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

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

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

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

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

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

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

2 結(jié)果與分析

2.1 插值結(jié)果

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

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

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

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

2.2 插值結(jié)果比較

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

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

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

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

3 結(jié)論與討論

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

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

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

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