基于栅格数据的降雨侵蚀力计算工具及其应用

夏照华，郭 浩，李瑞平，冯 阳，苏慧敏，赵海雷

(北京地拓科技发展有限公司，北京 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年全国逐月、年降雨侵蚀力，计算结果表明该工具具有良好的运算速度和计算精度。基于栅格数据的降雨侵蚀力计算工具开发，能为区域土壤侵蚀相关因子的快速更新与土壤侵蚀模型的推广奠定基础。