基于GIS的降雨侵蚀力计算器实现研究

殷 兵, 吴健平, 章文波

(1.华东师范大学 地理信息科学教育部重点实验室, 上海 200241; 2.北京师范大学 地理学与遥感科学学院, 北京 100875)

基于GIS的降雨侵蚀力计算器实现研究

殷 兵1, 吴健平1, 章文波2

(1.华东师范大学 地理信息科学教育部重点实验室, 上海 200241; 2.北京师范大学 地理学与遥感科学学院, 北京 100875)

降雨侵蚀力是土壤流失模型的一个基本因子，在收集和处理降雨数据时耗时耗力，且以文本格式保存的降雨数据存在诸多缺失问题。为此，将逐日、半月、累月及累年降雨侵蚀力计算模型与组件式GIS平台ArcGIS Engine及.NET平台结合，开发一种降雨侵蚀力因子计算器，并利用多种空间插值方法对数据进行不同的空间插值。结果表明：该计算器可快速、有效地的管理数据；对计算的降雨侵蚀力数据，使用插值方法可获得栅格数据，且能进一步进行单一或批量掩膜、重采样；同时也可对降雨侵蚀力因子进行统计分析、专题图制作等。为结果的可视化及保存提供应用基础。

降雨侵蚀力; 空间分析; 统计分析; 空间显示

降雨侵蚀力(Rainfall Erosivity)是指降雨引起土壤侵蚀的潜在能力，它是评价这种潜在能力的一个动力指标[1]。由于降雨过程资料难于获取，许多研究者提出了基于易获取的常规观测资料计算降雨侵蚀力的简易算法。在不同的简易算法中，使用的常规降雨资料各不相同，周伏建[2]、吴素业等[3]使用了逐年、月雨量资料；刘秉正[4]、马至尊[5]、Wischmeier等[6]则同时使用了逐年年雨量和月雨量资料，Arnoldus[7]、Renard等[8]、Yu等[9]使用的是多年平均年雨量和月雨量，还有应用逐年年雨量资料、日雨量资料、小时降雨量资料、多年一遇的时段雨量等[10-14]。

降雨侵蚀力的计算及其栅格数据的获取，是基于土壤侵蚀模型预报土壤流失量的关键所在。降雨数据常以文本文件格式保存，数据内容中存在很多缺失问题，在收集和处理这样大量的数据时非常困难且耗时。而应用GIS不仅可以方便地进行降雨侵蚀力计算中各参数的输入，还能将模型计算结果以图形方式输出，更形象、直观的表达计算结果，显示降雨侵蚀力的空间分布规律。范建容等[15]利用GIS空间分析功能获得长江上游降雨侵蚀力分布图、降雨侵蚀力年际变化图和各区域降雨侵蚀力R值平均年内分配曲线；李静[16]基于组件式GIS开发在延河流域水土保持效益评价系统中获取降雨侵蚀力因子的技术实现。

本文研究的系统以ArcGIS Engine和面向对象程序设计语言C#.NET为基础平台，进行降雨侵蚀力计算器的开发。该计算器可有效地解决数据的增、删、改，选择不同降雨侵蚀力计算模型，计算降雨侵蚀力，使用不同的空间插值方法得到降雨侵蚀力栅格数据，并根据需求对数据进行单一或批量掩膜、重采样；同时，可进一步对降雨侵蚀力进行最大值、最小值、平均值、标准差统计分析，生成曲线图或专题图，为结果的可视化及保存提供应用基础。

1 降雨侵蚀力计算器的总体设计

在降雨侵蚀力计算中，各气象站点数据格式往往不统一，需要对各气象站点数据进行统一化管理，并对基础的矢量数据进行管理。针对不同详细程度的降雨数据，本系统提供不同的降雨侵蚀力因子计算模型。

系统主要功能有数据库管理、降雨侵蚀力计算、空间分析、统计制图及空间显示，采用自顶向下、逐步求精、可拓展性及模块化、结构化的设计原则。系统界面直观，功能明确，操作简单。同时，本计算系统开发基础为PC机，内存2GB，CPU为Inter(R) Core(TM) i3，硬盘100GB以上；应用的操作系统为Windows XP+SP3，GIS组件采用ArcGISEngine9.3，系统开发平台采用Microsoft Visual Studio 2005，空间数据库采用的是ArcGIS的File Geodatabase。

2 各功能模块及实现

2.1 数据管理

降雨侵蚀力计算器系统以ArcGIS自带的File Geodatabase数据库方式对空间数据、属性数据进行存储、管理及维护，主要用来管理和导入基础地理背景资料，以获取研究区域所需的地理信息数据，包括数据的导入、导出、删除等。由于研究中生成的降雨侵蚀力栅格数据的数据量大，File Geodatabase存储空间有限，因此栅格数据将由文件的方式进行存储，方便拷贝，并同其他土壤侵蚀模型因子计算。

计算器可导入两种格式的降雨数据，一种是文本文件，另一种是Excel。导入到数据库中的表，表名为R+站点编号，字段分别是站点编号+年(+月+日)+值。数据导出的格式为Excel，并可对站点表进行增加、删除、更新等操作。

2.2 降雨侵蚀力计算

系统中所用的降雨侵蚀力计算模型，针对于不同详细程度的降雨数据：日雨量计算模型用于数据比较详细，精度相对较高的计算；半月计算模型、累月计算模型、累年计算模型应用于数据详细程度弱，精度相对偏低的计算。在计算时直接点击要选用的模型，程序将计算的降雨侵蚀力自动保存到数据库。计算完后，在数据输入中点击相对应的模型将产生的降雨侵蚀力写入到站点图层。系统中各计算模型所应用的公式主要为：

(1) 逐日计算模型[17]。

α=21.239β-7.3967

(2) 半月计算模型。

(EI30半月)=Ai·P半月·P半月1440

式中：(EI30)半月——半月降雨侵蚀力[(MJ·mm)/(hm2·h)]；P半月——半月内日雨量大于等于12mm的雨量之和(mm)；P半月1440——半月平均最大日雨量(mm)，Ai——区域i的拟合系数。

(3) 累月计算模型。

(EI30累月)=Bi·P累月·P累月1440

式中：(EI30)累月——多年平均月降雨侵蚀力[(MJ·mm)/(hm2·h)]；P累月——多年平均月降雨量(mm)；P累月1440——多年平均月最大日雨量(mm)；Bi——区域i的拟合系数。

(4) 累年计算模型。

(EI30)累年=Ci·P累年·P累年1440

式中：(EI30)累年——多年平均年降雨侵蚀力[(MJ·mm)/(hm2·h·a)]；P累年——多年平均年降雨量(mm)；P累年1440——多年平均年最大日雨量(mm)；Ci——区域i的拟合系数。

程序设计实现流程见图1。用户选择降雨侵蚀力计算模型，并通过数据库中站点图层矢量数据，获取所需计算的站点编号列表，循环遍历获取每个站点的日降雨资料，如果在数据库中存在该站点，则通过模型进行计算，若不存在将该站点信息保存到日志中。当所有站点结束后，判断信息保存中信息条数，如果大于1条记录，则提示用户导入缺失站点数据或删除站点图层中该站点。

图1 降雨侵蚀力计算流程

2.3 空间分析

空间分析包括空间插值及对插值结果进行空间处理，对栅格数据进行单一裁剪、重采样，同时也可进行批量裁剪、重采样。

2.3.1 空间插值

(1) 克里金插值[18]。克里金(Kriging)插值法又称空间自协方差最佳插值法。该方法基于包含自相关(即测量点之间的统计关系)的统计模型，包括普通克里金法和泛克里金法。

(2) 反距离权重插值[19]。反距离权重(IDW)插值使用一组采样点的线性权重组合来确定像元值。权重是一种反距离函数，进行插值处理的表面是具有局部因变量的表面。

(3) 样条函数插值[20]。利用最小化表面总曲率的数学函数来估计值，从而生产恰好经过输入点的平滑表面。采样点被拉伸到它们数量上的高度；样条函数折弯一个橡皮页，该橡皮页在最小表面总曲率的同时穿过这些输入点。在穿过采样点时，它将一个数学函数与指数量的最近输入点进行拟合。主要包括规则样条函数和张力样条函数。

用户根据需求选择插值方法，输入所选插值方法中所需参数，主要有站点图层、插值字段、输出像元大小、插值区域等各种差值方法共有参数，还有各种方法自带的参数，如克里金插值需选择克里金变异模型，反距离插值需要设置幂值等。程序将调用ArcGISEngine中的插值方法，对点数据进行插值生成栅格数据并保存。

2.3.2 空间处理 基于全国或大区域的空间插值结果，其分辨率不可设置的太高。当小范围研究如3km2以内的小流域时，需对大范围及低分辨率的数据进行裁剪、重采样。针对各个半月所占比例及多年平均年降雨侵蚀力栅格图并且裁剪范围不单一的情况，若进行单一裁剪、重采样，则工作重复，浪费时间和精力，因此本系统通过对栅格数据进行批量处理的方法来实现这些需求。

裁剪栅格和裁剪范围(重采样栅格和重采样范围)都有文件和目录两种选择：当都选择文件时，则对文件进行单一裁剪，输出栅格就保存为文件形式；当其中一个为目录时，这时对数据进行批处理，输出栅格将保存为一个目录，生成的数据名称将与裁剪栅格及裁剪范围相联系。

2.4 统计制图

对研究区各半月降雨侵蚀力所占比例数据进行统计分析制图，可对数据的可视化及应用提供方便。系统可对各半月降雨侵蚀力所占比例及年平均降雨侵蚀力栅格数据进行统计，得出最大值、最小值、平均值及标准差等。并可将统计结果导出成Excel格式数据，同样也可对各半月降雨侵蚀力所占比例的统计结果制成曲线图。

2.5 空间显示

系统中空间数据可以在数据视图中直接显示，并可实现基本的地图功能，如数据打开、放大缩小、全屏显示、平移、按比例放大缩小等。在制版视图中，可对数据进行制版设计，其包含的功能也主要有页面放大缩小、平移等。

3 降雨侵蚀力因子计算器的应用

本研究采用北京及周边地区的气象站点日降雨数据资料，运用克里金插值法对北京地区降雨侵蚀力因子进行计算，统计分析并制作专题图。

3.1 数据入库

(1) 降雨数据和站点数据入库。点击数据管理中批量导入降雨数据，选择相应路径后点击确定，数据库中将以R+站点编号格式存放降雨数据，包含站点编号、名称、经纬度，所在省的站点数据将以Station导入到数据库中。

(2) 矢量数据入库。点击数据管理中导入矢量数据，选择中国底图的存放路径，填写导入的名称。

3.2 降雨侵蚀力计算

(1) 站点侵蚀力表。点击降雨侵蚀力计算模块模型选择中的日雨量计算模型，计算结束后，在数据库中将生成以R+站点编号+m1命名的站点侵蚀力表。

(2) 数据写入。点击降雨侵蚀力计算模块数据输入中的日雨量计算模型，对数据库中站点侵蚀力表进行多年平均半月降雨侵蚀力统计，写入站点图层中对应字段中。

3.3 空间分析

选择空间分析模块空间插值中的克里金插值，相关参数设置为：点要素类为数据库中的站点图层，插值字段为降雨侵蚀力字段r，插值范围为数据库中的北京边界，选择结果存放在路径，及象元大小为30m×30m，其他默认，依次将24个半月降雨侵蚀力所占比例和多年平均降雨侵蚀力进行插值，即可得多年平均降雨侵蚀力。

3.4 统计制图

选择上述计算结果存储路径，点击“确定”，程序将对研究区北京24个半月降雨侵蚀力所占比例及多年平均年降雨侵蚀力进行最大值、最小值、平均值、标准差统计，点击“导出”可将统计结果导出。

若要制作曲线图，点击“曲线图”按钮，程序自动对统计值作图，可选择显示的统计值类型，并可导出图片格式。横坐标代表24个半月的名称，纵坐标代表24个半月所占统计值。将北京市多年平均年降雨侵蚀力进行专题制作，如附图17所示，颜色从蓝到红降雨侵蚀力值逐渐增加。

4 结论与探讨

本系统在ArcGISEngine9.3+vs2005+C#环境下，利用ArcGISEngine9.3提供的API，实现了降雨侵蚀力计算的软件化，为降雨侵蚀力计算提供了技术支持。系统实现了对降雨数据及基础地理数据的管理。将降雨数据进行统一化管理，为数据的标准化提供了基础。系统针对不同时段、不同详细程度的数据，提供多种降雨侵蚀力计算方法。同时，系统提供了多种空间插值方法，可对数据进行不同的空间插值，用户可对比不同插值方法结果精度，以便选择精度较高插值方法；空间处理可对数据进行批处理裁剪、重采样。另外，本计算系统可对栅格结果进行统计分析，制作专题图，生成曲线图，为结果的可视化提供基础。

研究对降雨侵蚀力的计算机可视化呈现进行了实践，通过软件的编写实现证明了系统的可行性。但在具体的实现中仍然存在着如下几个问题：ArcGISEngine9.3部分功能存在许多bug，处理大量数据时，程序可能出现意想不到的错误；系统开发时，运用FileGeodatabase，并将30m×30m的栅格数据以文件夹的存储方式存储，没有将此栅格数据放入到数据库中。系统没有对栅格数据专题图的制作进行不同种类选择，固定了栅格的渲染方式及比例尺，图例等形式。

以研究区北京市降雨侵蚀力为例，进行降雨侵蚀力的计算及可视化输出，虽达到了预期结果，但在数据库的选择上存在的缺陷，数据量变大情况下，可能会出现问题，并且没有对栅格结果保存到数据库中，没有实现数据统一化管理。因此希望在后来可以将数据库移植到ArcGISSDE空间数据库中，对初始数据、中间数据、结果数据进行统一管理；并可根据不同情况，增加降雨侵蚀力的算法，加强对栅格数据专题图制作的功能，最终实现更为完善数据管理，侵蚀力计算，空间分析及其可视化。

[1] 王万忠,焦菊英,郝小品,等.中国降雨侵蚀力R值的计算与分布(Ⅰ)[J].水土保持学报,1995,9(4):5-18.

[2] 周伏建,黄炎和.福建省降雨侵蚀力指标R值[J].水土保持学报,1995,9(1):13-18.

[3] 吴素业.安徽大别山区降雨侵蚀力简化算法与时空分布规律[J].中国水土保持,1994(4):12-14.

[4] 刘秉正.渭北地区R值的计算与分布[J].西北林学院学报,1993(2):21-29.

[5] 马至尊.应用卫星影像估算通用土壤流失方程各因子方法的探讨[J].中国水土保持,1989(3):24-27.

[6]WischmeierWH,SmithDD.Predictingrainfallerosionlosses-Aguidetoconservationplanning[R].Predictingrainfallerosionlosses:Aguidetoconservationplanning,1978.

[7]ArnoldusHMJ.MethodologyusedtodeterminethemaximumpotentialaverageannualsoillossduetosheetandrillerosioninMorocco[R].FAOSoilsBulletins(FAO),1977.

[8]RenardKG,FreimundJR.UsingmonthlyprecipitationdatatoestimatetheRfactorintherevisedUSLE[J].Journalofhydrology,1994,157(1):287-306.

[9]YuB,RosewellCJ.AnassessmentofadailyrainfallerosivitymodelforNewSouthWales[J].SoilResearch,1996,34(1):139-152.

[10] 李玉泉.基于日雨量的降雨侵蚀力模型研究[J].水利与建筑工程学报,2007,5(2):12-15.

[11] 谢云,章文波.用日雨量和雨强计算降雨侵蚀力[J].水土保持通报,2001,21(6):53-56.

[12] 殷水清.用小时降雨资料估算降雨侵蚀力的方法[J].地理研究,2007,26(3):541-547.

[13] Mikhailova E A, Bryant R B, Schwager S J, et al. Predicting rainfall erosivity in Honduras[J]. Soil Science Society of America Journal,1997,61(1):273-279.

[14] 张宪奎.黑龙江省土壤流失预报中R指标的研究[M]:北京:北京林业出版社,1992:63-68.

[15] 范建容,严冬,郭祥.GIS支持下的长江上游降雨侵蚀力时空分布特征分析[J].水土保持研究,2010,17(1):92-96.

[16] 李静.组件式GIS开发中获取降雨侵蚀力因子的技术实现:延河流域水土保持效益评价系统为例[J].安徽农业科学,2009,37(15):7278-7280.

[17] 国务院第一次全国水土普查领导小组办公室第一次全国水利普查培训教材之六:水土保持情况普查[M].北京:中国水利水电出版社,2010:209-210.

[18] Oliver M A, Webster R. Kriging: a method of interpolation for geographical information systems[J]. International Journal of Geographical Information System,1990,4(3):313-332.

[19] Watson D F, Philip G M. A refinement of inverse distance weighted interpolation[J]. Geo-processing,1985,2(4):315-327.

[20] Mitáš L, Mitášová H. General variational approach to the interpolation problem[J]. Computers & Mathematics with Applications,1988,16(12):983-992.

DevelopmentofRainfallErosivityCalculatorBasedonGIS

YIN Bing1, WU Jian-ping1, ZHANG Wen-bo2

(1.KeyLaboratoryofGeographicInformationScience,MinistryofEducation,EastChinaNormalUniversity,Shanghai200241,China; 2.SchoolofGeography,BeijingNormalUniversity,Beijing100875,China)

Rainfall erosivity is one of the essential factors in soil erosion models. Collection and processing of rainfall data are usually time-consuming. These data are generally stored in text format, which causes a common issue on data missing. This study developed a rainfall erosivity calculator by integrating the daily, half-monthly, monthly and yearly rainfall erosivity calculation models with ArcGIS Engine and .NET platform. Moreover, various spatial interpolation methods were implemented in this calculator in order to generate different interpolation results based on the collected rainfall data. The results demonstrate that this developed calculator is able to manage rainfall data quickly and effectively. Different interpolation methods can be employed to generate raster data for the spatial distribution of rainfall erosivity. These raster data can be masked or resampled individually or in batch by request. Statistics and thematic mapping of the rainfall erosivity can also be conducted by using this calculator, which provides the application basis for the visualization and storage of results.

rainfall erosivity; spatial analyst; statistical analysis; spatial visualization

：2014-04-22

第一次全国水利普查水土保持情况普查全国土壤侵蚀影响因子计算分析与制图第I标段水力侵蚀因子计算分析与制图

殷兵(1989—),男,江西南昌人,硕士研究生,研究方向为地理信息系统开发与遥感应用。E-mail:bingoyin@ecnu.cn

吴健平(1962—),男,浙江宁波人,博士,教授,主要从事地理信息系统开发与遥感应用。E-mail: jpwu@geo.ecnu.edu.cn

TV121

：A

：1005-3409(2014)04-0123-04