基于柵格數(shù)據(jù)的降雨侵蝕力計算工具及其應(yīng)用

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

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

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

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

1 降雨侵蝕力計算模型

1.1 水普法計算模型

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

α=21.239β-7.396 7

(2)

1.2 動能方程法計算模型

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

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

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

E=ep

(5)

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

(6)

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

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

2 降雨侵蝕力計算工具的設(shè)計與實現(xiàn)

2.1 總體設(shè)計及計算流程

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

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

2.2 計算步驟

2.2.1 輸入數(shù)據(jù)類型選擇

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

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

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

2.2.3 數(shù)據(jù)集規(guī)范性檢測

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

2.2.4 數(shù)據(jù)整理

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

2.2.5 降雨侵蝕力計算及匯總

輸入數(shù)據(jù)集與各模型計算所需條件匹配后，程序?qū)⒆詣佑嬎銓?yīng)的降雨侵蝕力。對于計算時段為日、月降雨侵蝕力的，程序會自動調(diào)用降雨侵蝕力匯總模塊，自動匯總計算年降雨侵蝕力。

2.3 技術(shù)實現(xiàn)

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

降雨侵蝕力計算工具編譯后為ENVI+IDL的擴(kuò)展模塊，使用時需借助ENVI軟件，不能獨立使用。該工具支持32位和64位操作系統(tǒng)，能一次性處理全國范圍的 TRMM 數(shù)據(jù)集或 CMORPH 數(shù)據(jù)集的數(shù)據(jù)?？紤]到用戶使用的方便性，設(shè)計時盡量做到錄入界面友好、清晰，輸出結(jié)果規(guī)范。

3 降雨侵蝕力計算工具的應(yīng)用

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

本次計算采用水普法計算模型，先對逐小時降雨量進(jìn)行求和獲取日降雨量柵格數(shù)據(jù)，日降雨量數(shù)據(jù)以月為單位存放，逐月計算月降雨侵蝕力，最后匯總為年降雨侵蝕力?；谀杲涤昵治g力柵格，在GIS軟件中采用空間統(tǒng)計分析方法，匯總?cè)珖饕》萁涤昵治g力最小值、最大值、均值等，統(tǒng)計結(jié)果見表1。

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

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

4 結(jié) 論

基于ENVI+IDL遙感軟件開發(fā)了降雨侵蝕力計算工具，支持水普法計算模型與動能方程法計算模型。該工具直接基于柵格數(shù)據(jù)對降雨侵蝕力進(jìn)行計算，所得結(jié)果可直接在ENVI、ArcMap等常見的RS與GIS軟件中使用?；谠摴ぞ撸褂?015年CMORPH逐小時降雨數(shù)據(jù)集，計算出2015年全國逐月、年降雨侵蝕力，計算結(jié)果表明該工具具有良好的運(yùn)算速度和計算精度。基于柵格數(shù)據(jù)的降雨侵蝕力計算工具開發(fā)，能為區(qū)域土壤侵蝕相關(guān)因子的快速更新與土壤侵蝕模型的推廣奠定基礎(chǔ)。