殷 兵, 吳健平, 章文波

(1.華東師范大學(xué) 地理信息科學(xué)教育部重點實驗室, 上海 200241; 2.北京師范大學(xué) 地理學(xué)與遙感科學(xué)學(xué)院, 北京 100875)

基于GIS的降雨侵蝕力計算器實現(xiàn)研究

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

(1.華東師范大學(xué) 地理信息科學(xué)教育部重點實驗室, 上海 200241; 2.北京師范大學(xué) 地理學(xué)與遙感科學(xué)學(xué)院, 北京 100875)

降雨侵蝕力是土壤流失模型的一個基本因子，在收集和處理降雨數(shù)據(jù)時耗時耗力，且以文本格式保存的降雨數(shù)據(jù)存在諸多缺失問題。為此，將逐日、半月、累月及累年降雨侵蝕力計算模型與組件式GIS平臺ArcGIS Engine及.NET平臺結(jié)合，開發(fā)一種降雨侵蝕力因子計算器，并利用多種空間插值方法對數(shù)據(jù)進行不同的空間插值。結(jié)果表明：該計算器可快速、有效地的管理數(shù)據(jù)；對計算的降雨侵蝕力數(shù)據(jù)，使用插值方法可獲得柵格數(shù)據(jù)，且能進一步進行單一或批量掩膜、重采樣；同時也可對降雨侵蝕力因子進行統(tǒng)計分析、專題圖制作等。為結(jié)果的可視化及保存提供應(yīng)用基礎(chǔ)。

降雨侵蝕力; 空間分析; 統(tǒng)計分析; 空間顯示

降雨侵蝕力(Rainfall Erosivity)是指降雨引起土壤侵蝕的潛在能力，它是評價這種潛在能力的一個動力指標(biāo)[1]。由于降雨過程資料難于獲取，許多研究者提出了基于易獲取的常規(guī)觀測資料計算降雨侵蝕力的簡易算法。在不同的簡易算法中，使用的常規(guī)降雨資料各不相同，周伏建[2]、吳素業(yè)等[3]使用了逐年、月雨量資料；劉秉正[4]、馬至尊[5]、Wischmeier等[6]則同時使用了逐年年雨量和月雨量資料，Arnoldus[7]、Renard等[8]、Yu等[9]使用的是多年平均年雨量和月雨量，還有應(yīng)用逐年年雨量資料、日雨量資料、小時降雨量資料、多年一遇的時段雨量等[10-14]。

降雨侵蝕力的計算及其柵格數(shù)據(jù)的獲取，是基于土壤侵蝕模型預(yù)報土壤流失量的關(guān)鍵所在。降雨數(shù)據(jù)常以文本文件格式保存，數(shù)據(jù)內(nèi)容中存在很多缺失問題，在收集和處理這樣大量的數(shù)據(jù)時非常困難且耗時。而應(yīng)用GIS不僅可以方便地進行降雨侵蝕力計算中各參數(shù)的輸入，還能將模型計算結(jié)果以圖形方式輸出，更形象、直觀的表達計算結(jié)果，顯示降雨侵蝕力的空間分布規(guī)律。范建容等[15]利用GIS空間分析功能獲得長江上游降雨侵蝕力分布圖、降雨侵蝕力年際變化圖和各區(qū)域降雨侵蝕力R值平均年內(nèi)分配曲線；李靜[16]基于組件式GIS開發(fā)在延河流域水土保持效益評價系統(tǒng)中獲取降雨侵蝕力因子的技術(shù)實現(xiàn)。

本文研究的系統(tǒng)以ArcGIS Engine和面向?qū)ο蟪绦蛟O(shè)計語言C#.NET為基礎(chǔ)平臺，進行降雨侵蝕力計算器的開發(fā)。該計算器可有效地解決數(shù)據(jù)的增、刪、改，選擇不同降雨侵蝕力計算模型，計算降雨侵蝕力，使用不同的空間插值方法得到降雨侵蝕力柵格數(shù)據(jù)，并根據(jù)需求對數(shù)據(jù)進行單一或批量掩膜、重采樣；同時，可進一步對降雨侵蝕力進行最大值、最小值、平均值、標(biāo)準(zhǔn)差統(tǒng)計分析，生成曲線圖或?qū)ｎ}圖，為結(jié)果的可視化及保存提供應(yīng)用基礎(chǔ)。

1 降雨侵蝕力計算器的總體設(shè)計

在降雨侵蝕力計算中，各氣象站點數(shù)據(jù)格式往往不統(tǒng)一，需要對各氣象站點數(shù)據(jù)進行統(tǒng)一化管理，并對基礎(chǔ)的矢量數(shù)據(jù)進行管理。針對不同詳細程度的降雨數(shù)據(jù)，本系統(tǒng)提供不同的降雨侵蝕力因子計算模型。

系統(tǒng)主要功能有數(shù)據(jù)庫管理、降雨侵蝕力計算、空間分析、統(tǒng)計制圖及空間顯示，采用自頂向下、逐步求精、可拓展性及模塊化、結(jié)構(gòu)化的設(shè)計原則。系統(tǒng)界面直觀，功能明確，操作簡單。同時，本計算系統(tǒng)開發(fā)基礎(chǔ)為PC機，內(nèi)存2GB，CPU為Inter(R) Core(TM) i3，硬盤100GB以上；應(yīng)用的操作系統(tǒng)為Windows XP+SP3，GIS組件采用ArcGISEngine9.3，系統(tǒng)開發(fā)平臺采用Microsoft Visual Studio 2005，空間數(shù)據(jù)庫采用的是ArcGIS的File Geodatabase。

2 各功能模塊及實現(xiàn)

2.1 數(shù)據(jù)管理

降雨侵蝕力計算器系統(tǒng)以ArcGIS自帶的File Geodatabase數(shù)據(jù)庫方式對空間數(shù)據(jù)、屬性數(shù)據(jù)進行存儲、管理及維護，主要用來管理和導(dǎo)入基礎(chǔ)地理背景資料，以獲取研究區(qū)域所需的地理信息數(shù)據(jù)，包括數(shù)據(jù)的導(dǎo)入、導(dǎo)出、刪除等。由于研究中生成的降雨侵蝕力柵格數(shù)據(jù)的數(shù)據(jù)量大，F(xiàn)ile Geodatabase存儲空間有限，因此柵格數(shù)據(jù)將由文件的方式進行存儲，方便拷貝，并同其他土壤侵蝕模型因子計算。

計算器可導(dǎo)入兩種格式的降雨數(shù)據(jù)，一種是文本文件，另一種是Excel。導(dǎo)入到數(shù)據(jù)庫中的表，表名為R+站點編號，字段分別是站點編號+年(+月+日)+值。數(shù)據(jù)導(dǎo)出的格式為Excel，并可對站點表進行增加、刪除、更新等操作。

2.2 降雨侵蝕力計算

系統(tǒng)中所用的降雨侵蝕力計算模型，針對于不同詳細程度的降雨數(shù)據(jù)：日雨量計算模型用于數(shù)據(jù)比較詳細，精度相對較高的計算；半月計算模型、累月計算模型、累年計算模型應(yīng)用于數(shù)據(jù)詳細程度弱，精度相對偏低的計算。在計算時直接點擊要選用的模型，程序?qū)⒂嬎愕慕涤昵治g力自動保存到數(shù)據(jù)庫。計算完后，在數(shù)據(jù)輸入中點擊相對應(yīng)的模型將產(chǎn)生的降雨侵蝕力寫入到站點圖層。系統(tǒng)中各計算模型所應(yīng)用的公式主要為：

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

α=21.239β-7.3967

(2) 半月計算模型。

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

式中：(EI30)半月——半月降雨侵蝕力[(MJ·mm)/(hm2·h)]；P半月——半月內(nèi)日雨量大于等于12mm的雨量之和(mm)；P半月1440——半月平均最大日雨量(mm)，Ai——區(qū)域i的擬合系數(shù)。

(3) 累月計算模型。

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

式中：(EI30)累月——多年平均月降雨侵蝕力[(MJ·mm)/(hm2·h)]；P累月——多年平均月降雨量(mm)；P累月1440——多年平均月最大日雨量(mm)；Bi——區(qū)域i的擬合系數(shù)。

(4) 累年計算模型。

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

式中：(EI30)累年——多年平均年降雨侵蝕力[(MJ·mm)/(hm2·h·a)]；P累年——多年平均年降雨量(mm)；P累年1440——多年平均年最大日雨量(mm)；Ci——區(qū)域i的擬合系數(shù)。

程序設(shè)計實現(xiàn)流程見圖1。用戶選擇降雨侵蝕力計算模型，并通過數(shù)據(jù)庫中站點圖層矢量數(shù)據(jù)，獲取所需計算的站點編號列表，循環(huán)遍歷獲取每個站點的日降雨資料，如果在數(shù)據(jù)庫中存在該站點，則通過模型進行計算，若不存在將該站點信息保存到日志中。當(dāng)所有站點結(jié)束后，判斷信息保存中信息條數(shù)，如果大于1條記錄，則提示用戶導(dǎo)入缺失站點數(shù)據(jù)或刪除站點圖層中該站點。

圖1 降雨侵蝕力計算流程

2.3 空間分析

空間分析包括空間插值及對插值結(jié)果進行空間處理，對柵格數(shù)據(jù)進行單一裁剪、重采樣，同時也可進行批量裁剪、重采樣。

2.3.1 空間插值

(1) 克里金插值[18]?？死锝?Kriging)插值法又稱空間自協(xié)方差最佳插值法。該方法基于包含自相關(guān)(即測量點之間的統(tǒng)計關(guān)系)的統(tǒng)計模型，包括普通克里金法和泛克里金法。

(2) 反距離權(quán)重插值[19]。反距離權(quán)重(IDW)插值使用一組采樣點的線性權(quán)重組合來確定像元值。權(quán)重是一種反距離函數(shù)，進行插值處理的表面是具有局部因變量的表面。

(3) 樣條函數(shù)插值[20]。利用最小化表面總曲率的數(shù)學(xué)函數(shù)來估計值，從而生產(chǎn)恰好經(jīng)過輸入點的平滑表面。采樣點被拉伸到它們數(shù)量上的高度；樣條函數(shù)折彎一個橡皮頁，該橡皮頁在最小表面總曲率的同時穿過這些輸入點。在穿過采樣點時，它將一個數(shù)學(xué)函數(shù)與指數(shù)量的最近輸入點進行擬合。主要包括規(guī)則樣條函數(shù)和張力樣條函數(shù)。

用戶根據(jù)需求選擇插值方法，輸入所選插值方法中所需參數(shù)，主要有站點圖層、插值字段、輸出像元大小、插值區(qū)域等各種差值方法共有參數(shù)，還有各種方法自帶的參數(shù)，如克里金插值需選擇克里金變異模型，反距離插值需要設(shè)置冪值等。程序?qū)⒄{(diào)用ArcGISEngine中的插值方法，對點數(shù)據(jù)進行插值生成柵格數(shù)據(jù)并保存。

2.3.2 空間處理 基于全國或大區(qū)域的空間插值結(jié)果，其分辨率不可設(shè)置的太高。當(dāng)小范圍研究如3km2以內(nèi)的小流域時，需對大范圍及低分辨率的數(shù)據(jù)進行裁剪、重采樣。針對各個半月所占比例及多年平均年降雨侵蝕力柵格圖并且裁剪范圍不單一的情況，若進行單一裁剪、重采樣，則工作重復(fù)，浪費時間和精力，因此本系統(tǒng)通過對柵格數(shù)據(jù)進行批量處理的方法來實現(xiàn)這些需求。

裁剪柵格和裁剪范圍(重采樣?xùn)鸥窈椭夭蓸臃秶?都有文件和目錄兩種選擇：當(dāng)都選擇文件時，則對文件進行單一裁剪，輸出柵格就保存為文件形式；當(dāng)其中一個為目錄時，這時對數(shù)據(jù)進行批處理，輸出柵格將保存為一個目錄，生成的數(shù)據(jù)名稱將與裁剪柵格及裁剪范圍相聯(lián)系。

2.4 統(tǒng)計制圖

對研究區(qū)各半月降雨侵蝕力所占比例數(shù)據(jù)進行統(tǒng)計分析制圖，可對數(shù)據(jù)的可視化及應(yīng)用提供方便。系統(tǒng)可對各半月降雨侵蝕力所占比例及年平均降雨侵蝕力柵格數(shù)據(jù)進行統(tǒng)計，得出最大值、最小值、平均值及標(biāo)準(zhǔn)差等。并可將統(tǒng)計結(jié)果導(dǎo)出成Excel格式數(shù)據(jù)，同樣也可對各半月降雨侵蝕力所占比例的統(tǒng)計結(jié)果制成曲線圖。

2.5 空間顯示

系統(tǒng)中空間數(shù)據(jù)可以在數(shù)據(jù)視圖中直接顯示，并可實現(xiàn)基本的地圖功能，如數(shù)據(jù)打開、放大縮小、全屏顯示、平移、按比例放大縮小等。在制版視圖中，可對數(shù)據(jù)進行制版設(shè)計，其包含的功能也主要有頁面放大縮小、平移等。

3 降雨侵蝕力因子計算器的應(yīng)用

本研究采用北京及周邊地區(qū)的氣象站點日降雨數(shù)據(jù)資料，運用克里金插值法對北京地區(qū)降雨侵蝕力因子進行計算，統(tǒng)計分析并制作專題圖。

3.1 數(shù)據(jù)入庫

(1) 降雨數(shù)據(jù)和站點數(shù)據(jù)入庫。點擊數(shù)據(jù)管理中批量導(dǎo)入降雨數(shù)據(jù)，選擇相應(yīng)路徑后點擊確定，數(shù)據(jù)庫中將以R+站點編號格式存放降雨數(shù)據(jù)，包含站點編號、名稱、經(jīng)緯度，所在省的站點數(shù)據(jù)將以Station導(dǎo)入到數(shù)據(jù)庫中。

(2) 矢量數(shù)據(jù)入庫。點擊數(shù)據(jù)管理中導(dǎo)入矢量數(shù)據(jù)，選擇中國底圖的存放路徑，填寫導(dǎo)入的名稱。

3.2 降雨侵蝕力計算

(1) 站點侵蝕力表。點擊降雨侵蝕力計算模塊模型選擇中的日雨量計算模型，計算結(jié)束后，在數(shù)據(jù)庫中將生成以R+站點編號+m1命名的站點侵蝕力表。

(2) 數(shù)據(jù)寫入。點擊降雨侵蝕力計算模塊數(shù)據(jù)輸入中的日雨量計算模型，對數(shù)據(jù)庫中站點侵蝕力表進行多年平均半月降雨侵蝕力統(tǒng)計，寫入站點圖層中對應(yīng)字段中。

3.3 空間分析

選擇空間分析模塊空間插值中的克里金插值，相關(guān)參數(shù)設(shè)置為：點要素類為數(shù)據(jù)庫中的站點圖層，插值字段為降雨侵蝕力字段r，插值范圍為數(shù)據(jù)庫中的北京邊界，選擇結(jié)果存放在路徑，及象元大小為30m×30m，其他默認，依次將24個半月降雨侵蝕力所占比例和多年平均降雨侵蝕力進行插值，即可得多年平均降雨侵蝕力。

3.4 統(tǒng)計制圖

選擇上述計算結(jié)果存儲路徑，點擊“確定”，程序?qū)ρ芯繀^(qū)北京24個半月降雨侵蝕力所占比例及多年平均年降雨侵蝕力進行最大值、最小值、平均值、標(biāo)準(zhǔn)差統(tǒng)計，點擊“導(dǎo)出”可將統(tǒng)計結(jié)果導(dǎo)出。

若要制作曲線圖，點擊“曲線圖”按鈕，程序自動對統(tǒng)計值作圖，可選擇顯示的統(tǒng)計值類型，并可導(dǎo)出圖片格式。橫坐標(biāo)代表24個半月的名稱，縱坐標(biāo)代表24個半月所占統(tǒng)計值。將北京市多年平均年降雨侵蝕力進行專題制作，如附圖17所示，顏色從藍到紅降雨侵蝕力值逐漸增加。

4 結(jié)論與探討

本系統(tǒng)在ArcGISEngine9.3+vs2005+C#環(huán)境下，利用ArcGISEngine9.3提供的API，實現(xiàn)了降雨侵蝕力計算的軟件化，為降雨侵蝕力計算提供了技術(shù)支持。系統(tǒng)實現(xiàn)了對降雨數(shù)據(jù)及基礎(chǔ)地理數(shù)據(jù)的管理。將降雨數(shù)據(jù)進行統(tǒng)一化管理，為數(shù)據(jù)的標(biāo)準(zhǔn)化提供了基礎(chǔ)。系統(tǒng)針對不同時段、不同詳細程度的數(shù)據(jù)，提供多種降雨侵蝕力計算方法。同時，系統(tǒng)提供了多種空間插值方法，可對數(shù)據(jù)進行不同的空間插值，用戶可對比不同插值方法結(jié)果精度，以便選擇精度較高插值方法；空間處理可對數(shù)據(jù)進行批處理裁剪、重采樣。另外，本計算系統(tǒng)可對柵格結(jié)果進行統(tǒng)計分析，制作專題圖，生成曲線圖，為結(jié)果的可視化提供基礎(chǔ)。

研究對降雨侵蝕力的計算機可視化呈現(xiàn)進行了實踐，通過軟件的編寫實現(xiàn)證明了系統(tǒng)的可行性。但在具體的實現(xiàn)中仍然存在著如下幾個問題：ArcGISEngine9.3部分功能存在許多bug，處理大量數(shù)據(jù)時，程序可能出現(xiàn)意想不到的錯誤；系統(tǒng)開發(fā)時，運用FileGeodatabase，并將30m×30m的柵格數(shù)據(jù)以文件夾的存儲方式存儲，沒有將此柵格數(shù)據(jù)放入到數(shù)據(jù)庫中。系統(tǒng)沒有對柵格數(shù)據(jù)專題圖的制作進行不同種類選擇，固定了柵格的渲染方式及比例尺，圖例等形式。

以研究區(qū)北京市降雨侵蝕力為例，進行降雨侵蝕力的計算及可視化輸出，雖達到了預(yù)期結(jié)果，但在數(shù)據(jù)庫的選擇上存在的缺陷，數(shù)據(jù)量變大情況下，可能會出現(xiàn)問題，并且沒有對柵格結(jié)果保存到數(shù)據(jù)庫中，沒有實現(xiàn)數(shù)據(jù)統(tǒng)一化管理。因此希望在后來可以將數(shù)據(jù)庫移植到ArcGISSDE空間數(shù)據(jù)庫中，對初始數(shù)據(jù)、中間數(shù)據(jù)、結(jié)果數(shù)據(jù)進行統(tǒng)一管理；并可根據(jù)不同情況，增加降雨侵蝕力的算法，加強對柵格數(shù)據(jù)專題圖制作的功能，最終實現(xiàn)更為完善數(shù)據(jù)管理，侵蝕力計算，空間分析及其可視化。

[1] 王萬忠,焦菊英,郝小品,等.中國降雨侵蝕力R值的計算與分布(Ⅰ)[J].水土保持學(xué)報,1995,9(4):5-18.

[2] 周伏建,黃炎和.福建省降雨侵蝕力指標(biāo)R值[J].水土保持學(xué)報,1995,9(1):13-18.

[3] 吳素業(yè).安徽大別山區(qū)降雨侵蝕力簡化算法與時空分布規(guī)律[J].中國水土保持,1994(4):12-14.

[4] 劉秉正.渭北地區(qū)R值的計算與分布[J].西北林學(xué)院學(xué)報,1993(2):21-29.

[5] 馬至尊.應(yīng)用衛(wèi)星影像估算通用土壤流失方程各因子方法的探討[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].水利與建筑工程學(xué)報,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] 張憲奎.黑龍江省土壤流失預(yù)報中R指標(biāo)的研究[M]:北京:北京林業(yè)出版社,1992:63-68.

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

[16] 李靜.組件式GIS開發(fā)中獲取降雨侵蝕力因子的技術(shù)實現(xiàn):延河流域水土保持效益評價系統(tǒng)為例[J].安徽農(nóng)業(yè)科學(xué),2009,37(15):7278-7280.

[17] 國務(wù)院第一次全國水土普查領(lǐng)導(dǎo)小組辦公室第一次全國水利普查培訓(xùn)教材之六:水土保持情況普查[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標(biāo)段水力侵蝕因子計算分析與制圖

殷兵(1989—),男,江西南昌人,碩士研究生,研究方向為地理信息系統(tǒng)開發(fā)與遙感應(yīng)用。E-mail:bingoyin@ecnu.cn

吳健平(1962—),男,浙江寧波人,博士,教授,主要從事地理信息系統(tǒng)開發(fā)與遙感應(yīng)用。E-mail: jpwu@geo.ecnu.edu.cn

TV121

：A

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