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

夏照華，郭 浩，李瑞平，馮 陽(yáng)，蘇慧敏，趙海雷

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

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

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

1 降雨侵蝕力計(jì)算模型

1.1 水普法計(jì)算模型

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

α=21.239β-7.396 7

(2)

1.2 動(dòng)能方程法計(jì)算模型

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

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

其中，次降雨總動(dòng)能E的計(jì)算公式[6]為

E=ep

(5)

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

(6)

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

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

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

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

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

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

2.2 計(jì)算步驟

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

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

2.2.2 降雨侵蝕力計(jì)算模型選擇

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

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

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

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

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

2.2.5 降雨侵蝕力計(jì)算及匯總

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

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

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

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

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

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

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

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

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

4 結(jié) 論

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