基于ANUSPLIN 的降水空間插值方法研究

李任君，高懋芳，李 強，李百壽

（1. 桂林理工大學測繪地理信息學院，廣西桂林541006；2. 中國農業科學院農業資源與農業區劃研究所/農業農村部農業遙感重點實驗室，北京100081；3. 江蘇師范大學地理測繪與城鄉規劃學院，徐州221000）

0 引言

降水是氣象觀測的重要因素之一，高精度的降水數據是氣象、農業干旱分析的重要基礎［1］。氣象數據表現為空間內的不均勻性、不連續性，通過站點測得的氣象數據，以一定的空間插值算法求解未知點的降水量是獲得降水空間數據的重要途徑［2-3］。插值過程中，由于降水具有間斷性和空間不連續性的特點［4］，影響因素眾多，且存在很強的季節變化，與氣溫空間插值相比，降水數據空間插值不確定性更大，更加復雜困難［5］。在所有插值方法中，基于統計插值技術的克里金法和薄盤光滑樣條函數法最為適用，這些技術建模只將空間分布作為觀測數據的函數而不需要其先驗知識和物理過程，能提高插值的準確度［6-8］。一些研究表明，利用平滑樣條法對降水的插值效果明顯優于其他方法［9-10］。錢永蘭將局部薄盤樣條插值結果分別與反向距離權重法和普通克里金法的插值結果進行對比，顯示ANUSPLIN 的插值誤差最小，適當增加站點數量和提高DEM 精度可進一步提高ANUSPLIN 軟件的插值精度［11］。劉正佳利用ANUSPLIN 軟件，以高程、距海岸線距離及高程—海岸線距離分別為協變量對多年平均降水量進行空間插值，發現以高程為協變量比以距海岸線距離和高程—海岸線距離為協變量所獲得的插值結果更加精確［12］。譚劍波采用克里金和ANUSPLIN 空間插值軟件對青藏高原東南緣年均溫度和年累計降水進行插值，發現降水插值中將高程作為協變量并不能體現ANUSPLIN 插值的優勢，但仍然優于克里金空間插值方法［13］。徐翔利用ANUSPLIN 軟件對中國康滇區多年月平均溫度和降水進行空間插值，并與克里金法和反距離加權法比較得出，ANUSPLIN 更適用于復雜山地環境下氣象要素的空間插值［14］。

雖然目前已有大量的研究在不同時空尺度上分析了我國不同地區和流域降水要素的插值，但是還沒有一種插值方法適用于所有的氣象要素，應根據研究目的和研究區域的自然地理地形特征選取適宜的插值方法及參數［15］。由于對平原地區的降水插值方法研究甚少，基于此，文章對黃淮海平原458 個氣象站點的多年平均降水數據進行處理，將連續365 d 的30 年平均降水數據每16 d 的降水累積量為一個試驗樣本，共得到23 個降水量數據樣本。以分辨率分別為90 m 和1 km 的數字高程模型（DEM）作為協變量，采用三變量局部薄盤樣條函數，通過ANUSPLIN 軟件進行降水插值。使用預留驗證集對該種插值方法在黃淮海地區降水插值的精度進行評價，并將插值結果與克里金插值方法進行對比，以期為相關科學研究提供理論依據和數據支撐。

1 材料與方法

1.1 研究區與數據來源

黃淮海地區是我國重要的糧、棉、油生產基地，對國民經濟發展起著舉足輕重的作用［16］。該地區屬溫帶大陸性季風氣候，冬季寒冷干燥，夏季高溫多雨，年平均降水量為500～600 mm，全年降水主要集中在夏季，秋、冬、春3 季均為水分虧缺期，經常出現先澇再旱，甚至連年干旱，是全國受旱面積最大的區域之一，嚴重制約著該地區農業生產力。

研究從中國氣象數據共享網（http：//cdc.cma.gov.cn/）收集了研究區內457 個氣象站點 1981 年 1 月至 2010 年 12 月每日（08：00- 次日 08：00）的降水量數據，選取 428 個氣象站點作為試驗站點輸入模型進行插值，預留29 個均勻分布在研究區的氣象站點作為驗證站點，研究區地理位置以及樣點分布見圖1。試驗中采用的90 m 和1 km 分辨率的數字高程模型數據（Digital Elevation Model，DEM）來源于CGIAR-CSI SRTM 中國區域數據（http：//srtm.csi.cgiar.org/），其具有精度高、應用廣泛、獲取免費等優點。

圖1 研究區DEM 及氣象站點分布Fig.1 DEM of the study area and distribution of meteorological stations

1.2 數據處理

應用ANUSPLIN 軟件，首先需要將降水數據通過SPSS 軟件處理成程序要求的標準格式，包括站點代碼、經度、緯度、高程信息和與氣象有關的屬性數據集，對數據進行檢查修正處理，剔除重復站點、缺值、錯值，數據格式詳見表1，以9 月30 日—10 月16 日樣本數據其中5 個站點的數據為例。ANUSPLIN 軟件所需要的DEM 數據需要轉換成ASCII 形式書寫，利用ArcGIS 軟件對DEM 數據進行拼接裁剪、投影變換為ALBERS投影，然后利用轉換工具轉為ASCII 格式文件，同時檢查和保證轉換過程中數據的正確性。

表1 1981—2010 年平均降水數據格式（rain.dat）Table 1 Data format of average precipitation from 1981 to 2010

1.3 模型檢驗方法

為了體現ANUSPLIN 對氣象數據插值的優勢，通過對氣象站點進行驗證，計算均方根誤差（RMSE）來評價ANUSPLIN 和常用的克里金插值方法的精度，公式為：

式（1）中，Zx和Zy分別代表降水量實測值和預測值，n為測站數。

ANUSPLIN 軟件提供了用于判別誤差來源和插值質量的統計參數，該研究選取廣義交叉驗證均方根（GCV）、殘差均方根（MSR）和方差估計均方根（VAR）等指標進行最優模型選取，具體公式詳見ANUSPLIN 的用戶指南［17］。

2 結果與分析

2.1 最優模型選取

由數據預處理部分得到ANUSPLIN 軟件特定格式的降水數據和DEM 數據。根據數據和ANUSPLIN 的用戶指南編寫運行ANUSPLIN 的SPLINA 和LAPGRD 這兩個程序模塊的批處理命令文件。APLINA 程序模塊負責生成曲面擬合結果的系數以及誤差統計文件，LAPGRD 程序模塊通過輸入SPLINA 程序模塊生成的表面系數文件得到插值表面。通過對前人工作經驗的總結［18-19］，該試驗選取了適用于降水插值的6 個模型（變量、協變量和樣條次數多種組合），以經緯度、高程作為變量對多年平均降水數據進行試驗，6 種模型參數設置見表2。ANUSPLIN 在日志文件（Log file and List file）中提供了一系列用于判別誤差來源和診斷插值質量的統計參數，如均值、方差、標準差、擬合曲面參數的有效數量估計Signal（信號自由度）、剩余自由度Error、光滑參數RHO、廣義交叉驗證GCV、期望真實均方誤差MSE、均方殘差MSR、方差估計VAR 及其平方根等，這些統計參數和模型診斷參數用來判斷站點是否充分，數據是否存在錯誤，是否存在合適的模型可以成功插值，幫助用戶輕松選擇出最佳空間分布模型［20-21］。利用黃淮海平原428 個站點的降水數據運行6 種模型，對產生的判別誤差來源和插值質量的統計參數進行整理選取最優模型，以9 月30 日至10 月16 日的數據為例，誤差統計參數見表3。

表2 六種待選薄盤光滑樣條函數模型Table 2 Six optional smooth spline function models for thin disks

表3 平均降水數據ANUSPLIN 插值輸出統計校驗Table 3 Comparisons of the ANUSPLIN output statistics for rainfall interpolation

依據最佳模型判斷標準：基于廣義交叉驗證或最大似然法，統計誤差最小，信噪比（SNR，信號自由度與剩余自由度之比）最小，信號自由度小于站點的1/2，SPLINA 生成的誤差文件中無*符號。Signal 表示擬合曲面的復雜程度，大于觀測站點的一半表示數據存在顯著錯誤或者是存在不適合表面模型的數據，擬合過程找不到最優光滑參數，這些情況在誤差統計文件中以*符號標出。可以根據插值產生的誤差統計文件，檢驗并消除原始數據在位置和數值上的錯誤，發現程序執行過程中的問題，并做出及時的調整［22-23］。對運行SPLINA 與LAPGRD 的6 種模型產生的誤差統計數據進行分析發現：高程作為協變量比高程作為獨立變量的插值精度要高；其他幾種模型的統計誤差數據都相差不大，DEM 為90 m 的插值精度普遍高于分辨率為1 km 插值。當DEM 分辨率為90 m 時，以經緯度為變量、高程為協變量，樣條次數設為4 的模型精度最高，所有數據都可以找到適合的表面模型，精度明顯高于樣條次數為2 或者3 的結果，其中模型6 為最優待用模型。

2.2 插值結果

應用ANUSPLIN 軟件，以90 m 分辨率的高程數據作為第三變量，選用模型6 為最優待用模型，對1981—2010 連續30 年23 個時段的平均降水數據樣本進行插值，在ArcGIS中將輸出表面文件轉換成tiff 圖像格式文件，進行顯示操作，結果見圖2。

圖2 1981 年1 月1 日至2010 年1 月4 日研究區全年降水分布圖（單位：mm）Fig.2 Distribution map of precipitation in the study area from January 1,1981 to January 4，2010

從圖2 空間分布上可看出西北部地區全年降水量較低，降水主要分布在南部地區，由南到北降水量呈遞減趨勢分布，5 月份開始，降水區域慢慢向北部擴大，8 月份時全區降水主要分布在北部和中部，南部降水量減少；之后降水分布主要在東部，西部降水減少，容易發生北旱南澇和東澇西旱，和實際的降水空間分布一致；時間分布上可看出黃淮海平原1 月份降水最少，7 月份降水最多，全年降水過分集中在7—8 月份，9 月份開始，全區降水量逐漸減少，發生先澇后旱的概率比較高，經常造成連年干旱。

2.3 插值精度分析

為了更好地體現ANUSPLIN 對氣象數據插值的優勢，利用常用的克里金插值方法對研究區6 月10—25 日期間的降水數據進行空間插值，如圖3 所示，通過ANUSPLIN 獲得的插值表面細節更加突出，能反映局部地形特征，使降水隨高程的梯度變化更加直觀，并且有比較好的平滑度。對冬季（1 月1—16 日）和夏季（6 月26—7 月12 日）的降水插值結果進行驗證，以均方根誤差（RMSE）為驗證幾種插值模型的標準，RMSE 通過插值獲得的驗證站點的降水數據與其觀測值計算所得，圖4 所示為冬季降水插值，采用ANUSPLIN 插值，DEM 分辨率分別為90 m、1 km 作為高程數據和克里金方法獲得的插值結果的RMSE 分別為：0.379 mm，0.376 mm，0.451 mm；圖5 為夏季降水插值，DEM 分辨率分別為90 m、1 km 作為高程數據和克里金方法3 種插值結果的RMSE 分別為：4.186 mm，4.237 mm，4.310 mm。

圖3 基于克里金和ANUSPLIN 的降水插值柵格（a：克里金插值；b：ANUSPLIN 插值）Fig.3 Precipitation interpolation grid based on Kriging and ANUSPLIN

根據降水插值的誤差比較結果得出：（1）無論是夏季降水還是冬季降水，ANUSPLIN的插值精度都高于克里金插值方法，夏季降水插值精度明顯低于冬季降水插值；（2）90 m分辨率和1 km 分辨率的DEM 分別作為ANUSPLIN 的協變量輸入時，在冬季降水量較少的數據插值結果中差別不是很大，在夏季降水插值中，以90 m 分辨率的DEM 作為高程輸入模型獲得的插值精度更高，DEM 柵格特征更明顯。

圖4 冬季降水插值結果和觀測值的散點分布（a）ANUSPLIN 插值-90 mDEM；（b）ANUSPLIN 插值-1 kmDEM；（c）克里金插值Fig.4 The scatterplot of true value and predictive value of the precipitation interpolation in winter

圖5 夏季降水插值結果和觀測值的散點分布（a）ANUSPLIN 插值-90 mDEM（b）ANUSPLIN 插值-1 kmDEM（c）克里金插值Fig.5 The scatterplot of true value and predictive value of the precipitation interpolation in summer

3 結論

由于不同區域地形復雜多樣，氣候差異性大，氣象站數量有限，影響降水數據插值精度的因素眾多，插值使用的模型和方法都需要根據研究區域實際情況進行選擇。該研究利用國內外廣泛應用的專用氣象插值軟件ANUSPLIN，以高程作為插值變量因子對黃淮海平原連續30 年地降水數據進行插值，得到與遙感數據的時空分辨率相匹配的柵格數據，對插值數據進行驗證，得出結論：（1）與克里金方法相比，ANUSPLIN 插值精度更高，可以更細致地描繪出黃淮海平原地形因素對降雨空間分布的影響，分辨率為90 m 的DEM 作為高程輸入數據所獲得的插值結果比分辨率為1 km 的插值結果精度相對高，但是前者試驗時間明顯長于后者；（2）夏季降水量集中的6—8 月份的插值精度比降水量較少的冬季降水插值精度低，雖然提高DEM 的分辨率，可以使夏季降水插值精度有所提升，但是并沒有改善降雨量的增加與精度誤差呈正相關的關系。通過使用ANUSPLIN 對黃淮海平原的降水插值方法的研究，探討插值精度與DEM 分辨率的關系，獲得較高精度的平均降水柵格數據，為黃淮海平原干旱分析并指導當地農作物灌溉生產提供重要的科學依據。