基于柵格數據的降雨侵蝕力計算工具及其應用

夏照華，郭 浩，李瑞平，馮 陽，蘇慧敏，趙海雷

(北京地拓科技發展有限公司，北京 100084)

我國是世界上土壤侵蝕最嚴重的國家之一[1]。根據第一次全國水利普查數據[2]，我國現有水力侵蝕面積為129.32萬km2，水土流失直接影響到社會經濟的可持續發展。土壤侵蝕的產生是多重自然、社會因素相互作用的結果。降雨是自然因素之一，而降雨侵蝕力是評價降雨引起土壤侵蝕的潛在能力的一個動力指標[1]，也是土壤侵蝕模型的基礎因子之一。降雨侵蝕力的計算方法主要有兩種：一種是依據日、月、年等時段的降雨量采用經驗公式推算降雨侵蝕力，以第一次全國水利普查降雨侵蝕力計算模型(以下簡稱“水普法計算模型”)為典型代表，采用這種方法的大多采用氣象站點數據，通過空間插值獲取研究區域的降雨侵蝕力，如謝云等[3-4]的研究成果。另一種是采用通用土壤流失方程中確定的EI30指標來計算(以下簡稱“動能方程法計算模型”)，由于降水過程資料不易獲取，因此采用這種方法一般使用時段最大降雨強度來替代最大30 min降雨強度進行計算。關于這種方法的研究目前多集中于基于TRMM[5]，采用動能方程法計算降雨侵蝕力，如朱強等[6]用最大3 h降雨強度替代最大30 min降雨強度計算了大凌河流域的降雨侵蝕力，靳秋桐等[7]對朱強的計算公式進行了改進，提出了降雨侵蝕力的校正方法，并將其應用于薊運河上游地區。

在降雨侵蝕力計算模型的軟件化方面，部分研究人員嘗試開發降雨侵蝕力的計算工具，在一定程度上提高了降雨侵蝕力的計算速度，避免了人工計算的煩瑣，如殷兵等[8]的降雨侵蝕力計算器。但目前公開發表的文獻中，大多是基于氣象站點矢量數據或者局部區域的柵格數據的降雨侵蝕力計算工具，很少涉及關于全國大范圍柵格數據的降雨侵蝕力計算工具的研究。基于上述學者的研究成果，本研究介紹了基于柵格數據的降雨侵蝕力計算工具的設計思路及其實現方法，并運用該工具基于“中國自動站與CMORPH融合的逐時降水量0.1°網絡數據集(1.0版)”(以下簡稱“CMORPH數據集”)計算了2015年全國逐月、年降雨侵蝕力。

1 降雨侵蝕力計算模型

1.1 水普法計算模型

根據逐日降雨量數據，剔除日降雨量小于12 mm的非侵蝕性降雨后，分別計算日、月、年降雨侵蝕力。月降雨侵蝕力計算公式[9]為

α=21.239β-7.396 7

(2)

1.2 動能方程法計算模型

RUSLE手冊建議采用30min降雨數據計算降雨侵蝕力，計算公式[6-7]為

式中：R為年平均降雨侵蝕力，MJ·mm/(hm2·h·a)；n為總年數；j為參與計算的年份序數；m為第j年內的暴雨次數；k為參與計算的次暴雨序數；E為次降雨總動能，MJ/hm2；I30為30min時間分辨率資料統計的次降雨最大30min降雨強度，mm/h。

其中，次降雨總動能E的計算公式[6]為

E=ep

(5)

e=0.29×[1-0.72exp(-0.082I)]

(6)

上兩式中：e為單位降雨動能，MJ/(hm2·mm)；p為降雨量，mm；I為降雨強度，mm/h。

采用柵格數據計算降雨侵蝕力時，可使用次降雨最大時段(1h、3h)雨強替代最大30min雨強進行計算。

2 降雨侵蝕力計算工具的設計與實現

2.1 總體設計及計算流程

本研究設計降雨侵蝕力計算工具的目的是實現基于柵格數據的降雨侵蝕力因子的快速計算，便于用戶使用。輸入數據為逐日、月時間序列的降雨強度或降雨量柵格數據集。計算流程如圖1所示。

圖1 降雨侵蝕力計算工具總體設計及計算流程

2.2 計算步驟

2.2.1 輸入數據類型選擇

根據輸入柵格數據的類型選擇對應的入口，降雨侵蝕力計算工具能支持TRMM數據集、CMORPH數據集，以及通過氣象站點數據插值生成的日、月、年時間序列數據集，所有數據集均應按指定要求命名與存放。

2.2.2 降雨侵蝕力計算模型選擇

降雨侵蝕力計算工具設計由用戶手動選擇所要采用的降雨侵蝕力模型，支持水普法計算模型與動能方程法計算模型，還預留有接口，可后期擴展支持常用土壤侵蝕模型。

2.2.3 數據集規范性檢測

根據用戶輸入的數據集與選擇的計算模型，程序將自動根據內置規則判斷是否滿足計算要求。水普法計算模型需要輸入的數據為日降雨量。長時間序列的日降雨量的命名必須滿足“數據集標識符+降雨量時段標識符+日期”，如CMORPH的日降雨量命名為 CMORPH_DAY_20150410.tif。柵格數據格式可以為 *.tif 或*.dat等 ENVI 支持的格式，輸入時間序列的文件擴展名與命名方式建議保持一致。動能方程法計算模型需要輸入的數據為時段降雨量。TRMM數據與CMORPH數據的時段降雨量文件的命名直接采用各數據集的命名規范，命名時應包含有時間、數據集標識等信息，統一轉換為 ENVI軟件能支持的格式。程序能自動識別不滿足計算要求的數據，為避免無效數據參與計算，當用戶輸入數據不滿足要求時，程序會自動終止計算過程，并提示用戶對數據進行整理后重新輸入計算。

2.2.4 數據整理

對于不滿足要求的數據，根據提示，用戶可使用數據整理模塊對相關數據集進行整理：格式未統一的，可使用格式轉換模塊統一對指定文件夾下的數據進行轉換；選擇模型與輸入數據時段不一致的，如選擇水普法計算模型，但輸入的是時段降雨量數據，會提示用戶“輸入的數據集不滿足模型要求”，由用戶基于輸入的時段降雨量數據改用日降雨量計算模塊批量計算日降雨量后，再次輸入日降雨量調用水普法計算模型進行降雨侵蝕力計算；數據存放不規范的，如采用動能方程法計算模型批量計算月降雨侵蝕力的時候，數據必須按月存放，若發現某一月份的文件夾中存有其他月份的文件，則會提示用戶使用數據整理模塊對降雨數據進行自動整理。

2.2.5 降雨侵蝕力計算及匯總

輸入數據集與各模型計算所需條件匹配后，程序將自動計算對應的降雨侵蝕力。對于計算時段為日、月降雨侵蝕力的，程序會自動調用降雨侵蝕力匯總模塊，自動匯總計算年降雨侵蝕力。

2.3 技術實現

本研究降雨侵蝕力計算工具的開發是基于 ENVI+IDL環境，柵格文件的讀寫可通過調用 ENVI 的柵格數據讀寫函數實現，一次性將所需數據讀入到數組，將長時間序列的多個柵格圖層逐像元讀入到三維數組(分別為行、列、時間維度)。計算降雨侵蝕力的時候，以行、列作為循環控制變量，逐像元計算降雨侵蝕力，計算結果寫入柵格數據存儲的二維數組，計算結果調用ENVI 軟件的柵格寫入程序，輸出為*.tif或*.dat等常用的GIS軟件與RS軟件支持的格式。

降雨侵蝕力計算工具編譯后為ENVI+IDL的擴展模塊，使用時需借助ENVI軟件，不能獨立使用。該工具支持32位和64位操作系統，能一次性處理全國范圍的 TRMM 數據集或 CMORPH 數據集的數據。考慮到用戶使用的方便性，設計時盡量做到錄入界面友好、清晰，輸出結果規范。

3 降雨侵蝕力計算工具的應用

以CMORPH數據集作為輸入數據，使用降雨侵蝕力計算工具計算了2015年全國逐月、年降雨侵蝕力。

本次計算采用水普法計算模型，先對逐小時降雨量進行求和獲取日降雨量柵格數據，日降雨量數據以月為單位存放，逐月計算月降雨侵蝕力，最后匯總為年降雨侵蝕力。基于年降雨侵蝕力柵格，在GIS軟件中采用空間統計分析方法，匯總全國主要省份降雨侵蝕力最小值、最大值、均值等，統計結果見表1。

表1 2015年全國主要省(區、市)降雨侵蝕力統計值 MJ·mm/(hm2·h·a)

計算過程中，CMORPH數據集使用地面和衛星兩個來源的降水數據。地面觀測降水資料來自全國3萬多個自動觀測站(包括國家級自動站和區域自動站)逐小時降水量。衛星反演降水產品選用由美國環境預測中心的氣候預測中心開發的實時衛星反演CMORPH降水產品，原始CMORPH資料的空間分辨率為8 km，時間分辨率為30 min。將原始時空分辨率的CMORPH數據空間重采樣得到1 h、0.1°×0.1°的衛星反演降水產品。

4 結 論

基于ENVI+IDL遙感軟件開發了降雨侵蝕力計算工具，支持水普法計算模型與動能方程法計算模型。該工具直接基于柵格數據對降雨侵蝕力進行計算，所得結果可直接在ENVI、ArcMap等常見的RS與GIS軟件中使用。基于該工具，使用2015年CMORPH逐小時降雨數據集，計算出2015年全國逐月、年降雨侵蝕力，計算結果表明該工具具有良好的運算速度和計算精度。基于柵格數據的降雨侵蝕力計算工具開發，能為區域土壤侵蝕相關因子的快速更新與土壤侵蝕模型的推廣奠定基礎。