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

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

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

0 引言

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

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

1 材料與方法

1.1 研究區(qū)與數(shù)據(jù)來源

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

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

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

1.2 數(shù)據(jù)處理

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

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

1.3 模型檢驗方法

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

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

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

2 結(jié)果與分析

2.1 最優(yōu)模型選取

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

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

表3 平均降水數(shù)據(jù)ANUSPLIN 插值輸出統(tǒng)計校驗Table 3 Comparisons of the ANUSPLIN output statistics for rainfall interpolation

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

2.2 插值結(jié)果

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

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

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

2.3 插值精度分析

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

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

圖4 冬季降水插值結(jié)果和觀測值的散點分布（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 夏季降水插值結(jié)果和觀測值的散點分布（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 結(jié)論

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