基于栅格GIS土壤侵蚀地形因子的提取算法

226基于栅格基于栅格基于栅格基于栅格GISGISGISGIS土壤侵蚀地形因子的提取算法土壤侵蚀地形因子的提取算法土壤侵蚀地形因子的提取算法土壤侵蚀地形因子的提取算法张照录张照录张照录张照录1崔继红崔继红崔继红崔继红21.山东理工大学资源与环境工程学院淄博2550492.华中农业大学经贸土管学院武汉430070摘摘摘摘要要要要:论述了基于栅格GIS土壤侵蚀地形因子的计算方法和流程在ARCGIS栅格模块中利用AML开发了通用土壤流失方程地形因子的计算程序并给出了程序核心部分的编码实际应用表明计算结果的精度依赖于数据源的精度关键词关键词关键词关键词栅格GIS土壤侵蚀通用土壤流失方程地形因子AMLAlgorithmofTopographicFactorsforSoilErosionBasedonGridGISZHANGZhaolu1,CUIJihong2(1.FacultyofResourcesandEnvironmentalEngineering,ShandongUniversityofTechnology,Zibo255049;2.FacultyofEconomicsandTradeInternational–landManagement,HuazhongAgriculturalUniversity,Wuhan430070)AbstractThisarticlediscussesthealgorithmandtheflowchartoftopographicfactorsforsoilerosionbasedongridGIS.UndertheenvironmentofARCGISgrid,theprogramofcomputingtopographicfactorsforuniversalsoillossequationisdevelopedwitharcmacrolanguage(AML).Atthesametime,thecorecodesoftheprogramaresupplied.Therealapplicationindicatesthattheprecisionoftheresultsdependsontheprecisionofdatasources.KeywordsGridGIS;Soilerosion;Universalsoillossequation(USLE);Topographicfactors;AML计计计计算算算算机机机机工工工工程程程程ComputerEngineering第第第第32卷卷卷卷第第第第5期期期期Vol.3252006年年年年3月月月月March2006工程应用技术与实现工程应用技术与实现工程应用技术与实现工程应用技术与实现文章编号文章编号文章编号文章编号10003428(2006)05022603文献标识码文献标识码文献标识码文献标识码A中图分类号中图分类号中图分类号中图分类号TP306.6在土壤侵蚀研究领域中众多的土壤侵蚀模型如USLERUSLEANSWERSAGNPSWEPP等都需要地形因子作为重要的基础数据地形因子主要是指坡度因子S和坡长因子L由坡度和坡长l来计算传统的坡度和坡长要靠野外实地测量来获取这在研究区域较大时往往是不现实的地理信息系统软件的发展为空间问题的研究提供了强大的支撑平台在ESRI公司ARCGIS地理信息系统软件的栅格模块GRID中利用其开发语言AML进行二次开发理论上可以实现基于栅格数据运算的任一尺度的坡度和坡长的自动提取从而实现坡度因子和坡长因子的自动生成本文以USLE模型中坡度因子和坡长因子的求取过程为例详细介绍整个算法的过程及其实现1算法描述算法描述算法描述算法描述1.1问题的由来问题的由来问题的由来问题的由来土壤侵蚀是发生在地球表面的一种自然现象不同地点的土壤侵蚀等级是不同的土壤侵蚀的强弱程度受多种因子的共同影响其中地形因子就是其中之一坡度越大坡长越长土壤侵蚀的强度就越大这种影响能力通过坡度因子和坡长因子来量化一个地区的坡度和坡长是连续变化的为了问题的简化通常把整个研究区域划分为一系列大小相等的栅格认为每个栅格单元的坡度和坡长是相同的这样一来区域坡度因子和坡长因子的计算转化为求取每个栅格单元的坡度因子和坡长因子1.2基本概念基本概念基本概念基本概念基于栅格的坡度计算方法有很多其结果也有较大差异[1]本文采用最大坡降坡度法[3]计算最大坡降方向的坡度值栅格单元的最大坡降方向定义为该栅格单元标号为9和其邻域栅格单元标号为1-8高度差值最大的方向用12−n表示n为邻域栅格单元的标号图1如果有多个方向的差值相等且最大则进一步扩大邻域范围直到确定唯一的最大坡降方向最大坡降坡度的计算公式为()−=Lzzasi9maxtans为坡度iz为栅格单元的高度值L为两个相邻栅格单元中心点之间的距离当位于1357位置时L=d当位于对角线位置2468时dL×=2d为栅格单元的分辨率1,2,3,...,8i=678591432图图图图1坡度计算窗口及栅格单元编码坡度计算窗口及栅格单元编码坡度计算窗口及栅格单元编码坡度计算窗口及栅格单元编码USLE中坡长的定义为从坡面漫流的起点到下列两类地点之间的水平距离1由于坡度降低出现沉积的地区2坡面漫流汇聚成汇水区的地区[4]因此计算坡长首先找出坡面漫流的起点即地形高点如山脊和漫流汇聚的地区如河道湖泊等集水区坡度降低出现沉积的地区是由截止角度来表征的截止角度定义为沿最大坡降方向相邻像元的坡度变化率其取值范围为0~1如果坡度降低的速率大于截止角度沉积就会出现该单元的坡长值等于0随着坡长的增加侵蚀是不断增加的因此每个栅格单元的坡长值用累积坡长值来表征累积坡长就是从地形高点开始沿着最大坡降方向不断累加形成的1.3算法的流程算法的流程算法的流程算法的流程算法的基础数据为数字高程数据DigitalElevation作者简介作者简介作者简介作者简介张照录1976男博士生讲师主研方向地理信息系统开发与应用崔继红硕士生收稿日期收稿日期收稿日期收稿日期[e-mail]227Model,DEM和水系数据获取DEM的方法有多种最经济最常用的当属数字化地形图由等高线数据转换而来DEM数据要进行预处理即填平处理其主要目的是为了防止坡度计算过程中出现负值负值代表非侵蚀填平后其值为0并不改变该地点的实际地形意义完成预处理之后计算最大坡降方向和最大坡降坡度然后计算每个栅格单元的非累积坡长找出地形高点水系数据的作用是确定汇水区即坡长的终点水系数据经矢栅转换后与地形高点栅格数据进行关系运算得到地形高点和汇水区的叠加数据即确定了坡长起点和终点的栅格数据截止角度的确定一般由对研究区比较熟悉的专家确定默认情况下为0.5具有了以上多种数据之后就可以计算累积坡长并最终生成USLE模型中要求的坡度因子和坡长因子图2DEM生成确定截止角度水系数据截止角度填平DEM水系矢量数据最大坡降方向最大坡降角度水系栅格数据非累积坡长地形高点高点汇水区累积最大坡降坡长坡度坡长因子LS图图图图2基于栅格基于栅格基于栅格基于栅格GIS土壤侵蚀地形因子的算法流程土壤侵蚀地形因子的算法流程土壤侵蚀地形因子的算法流程土壤侵蚀地形因子的算法流程2算法的实现及实例算法的实现及实例算法的实现及实例算法的实现及实例DEM数据的填平处理由GRID模块FILL命令完成最大坡降方向由函数FLOWDIRECTION确定需要说明的是由FLOWDIRECTION计算的水流方向和最大坡降方向是相同的最大坡降坡度的计算按前述公式进行程序较为简单非累积坡长的计算也比较简单如果最大坡降方向为对角线方向则该栅格单元的非累积坡长为栅格单元长度的1.4倍否则为栅格单元的长度2.1地形高点及汇水区栅格数据的生成地形高点及汇水区栅格数据的生成地形高点及汇水区栅格数据的生成地形高点及汇水区栅格数据的生成地形高点定义为没有其它栅格单元的坡面漫流流入或者有其它单元的坡面漫流流入但其高程值相等的栅格单元地形高点的确定是逐个栅格单元进行的计算窗口采用33大小图1用FOCALFLOW函数确定每个栅格单元可能的坡面漫流流入方向如果邻域栅格单元的高度值大于中心栅格单元的高度值则该方向为可能的坡面漫流流入方向之一其标志值为12−nn为邻域栅格单元编号否则为0邻域栅格单元全部标志值之和赋给中心栅格单元作为中心栅格单元的可能的坡面漫流流入标志如果任意一个方位有坡面漫流流入的可能且该方位邻域栅格单元的最大坡降方向指向该单元则该栅格单元不是地形高点其值设为空值如果坡度为0无侵蚀发生其累积坡长值为0如果是地形高点认为坡面漫流的起点在栅格单元中心只有一半距离的栅格单元区存在侵蚀故累积坡长为该栅格单元非累积坡长的一半sl_inflow为坡面漫流流入方向栅格数据sl_flow_buf为带一个栅格宽度缓冲区的最大坡降方向栅格数据sl_ncsl为非累积坡长栅格数据sl_hipts_nodata为地形高点栅格数据sl_slope为最大坡降坡度栅格数据if(%sl_slope%eq0)%sl_hipts_nodata%=0.0elseif((%sl_inflow%&&64)and(sl_flow_buf(0,-1)eq4))%sl_hipts_nodata%=setnull(1eq1)elseif((%sl_inflow%&&128)and(sl_flow_buf(1,-1)eq8))%sl_hipts_nodata%=setnull(1eq1)elseif((%sl_inflow%&&1)and(sl_flow_buf(1,0)eq16))%sl_hipts_nodata%=setnull(1eq1)elseif((%sl_inflow%&&2)and(sl_flow_buf(1,1)eq32))%sl_hipts_nodata%=setnull(1eq1)elseif((%sl_inflow%&&4)and(sl_flow_buf(0,1)eq64))%sl_hipts_nodata%=setnull(1eq1)elseif((%sl_inflow%&&8)and(sl_flow_buf(-1,1)eq128))%sl_hipts_nodata%=setnull(1eq1)elseif((%sl_inflow%&&16)and(sl_flow_buf(-1,0)eq1))%sl_hipts_nodata%=setnull(1eq1)elseif((%sl_inflow%&&32)and(sl_flow_buf(-1,-1)eq2))%sl_hipts_nodata%=setnull(1eq1)else%sl_hipts_nodata%=0.5*%sl_ncsl%endif水系数据与地形高点数据进行关系运算的规则为如果水系栅格数据为空值该栅格单元值为地形高点栅格数据的相应栅格单元值如果非空则为0其表达式为%sl_cum_l%=con(isnull(waterarea),%sl_hipts_nodata%,0)得到的栅格数据sl_cum_l确定了坡长起点和作为坡长终点之一的汇水区2.2累积坡长的计算累积坡长的计算累积坡长的计算累积坡长的计算累积最大坡降坡长所需要的基础数据已全部具备只需要从地形高点开始按照最大坡降方向进行累加计算即可这里有两条假定1在不同方向坡面漫流汇聚地区以最长的累积坡长为主导2当相邻栅格单元沿最大坡降方向的坡度减少速率大于截止角度时处于下方的栅格单元没有侵蚀坡长应设为0累积坡长从此处重新起算以sl_cum_l为基础从8个不同方向分别沿最大坡降方向累加坡长每次累加结束后选择最大值作为栅格单元的累积坡长值直到栅格数据中没有空值即每个栅格单元都有累积坡长为止以下为该程序的核心代码/*测试该单元是否可能有上游单元作为水流来源如果没有那么该单元将被赋0值如果有可能进行下一步测试if(not(%sl_inflow%&&%possible_cell_direction%))%from_cell_grid%=0