CN115239894A - 一种适用于大尺度地理坐标系栅格数据的ls因子提取方法 - Google Patents

Abstract

本发明涉及一种适用于大尺度地理坐标系栅格数据的LS因子提取方法，适用于地理坐标系栅格数据的LS因子提取方法设计和实现，能提升LS因子提取效率，并且设计的针对大尺度LS因子的计算流程，完善了大尺度LS因子提取的方式；该方法可快速有效地对基于高分辨率全球尺度SRTM1进行LS因子的提取与计算；包括如下步骤：步骤1、数据分块：步骤2、添加缓冲区并将栅格数据格式转化为ASCII数据：步骤3、带缓冲区数据的LS因子提取：步骤4、ASCII结果数据格式转化并提取不带缓冲区各LS因子数据：步骤5、数据融合提取所需尺度的地理坐标系栅格LS因子；将步骤4得到的LS因子数据通过Arcmap软件中的镶嵌工具融合即可得到所需尺度的地理坐标系下栅格LS因子。

Description

一种适用于大尺度地理坐标系栅格数据的LS因子提取方法

技术领域

本发明涉及计算机科学、侵蚀学与地理学交叉的数字地形分析技术，具体地说是一种适用于大尺度地理坐标系栅格数据的LS因子提取方法。

背景技术

地形是影响土壤侵蚀的重要因子，在美国的通用土壤流失方程和修正通用土壤流失方程式(Universal soil loss equation/revised universal soil loss equation，USLE/REUSLE)以及中国土壤流失方程(Chinese Soil Loss Equation， CSLE)中，地形(LS)因子是用于土壤流失定量计算的重要指标。目前用于区域尺度LS因子的提取算法主要基于投影坐标系下的DEM(Digital Elevation Model)栅格数据进行，而用于大尺度范围基于地理坐标系下栅格数据(如SRTM、ASTER GDEM、AW3D30 DEM等)进行LS因子提取的研究鲜有报道。数据以1弧秒分辨率(约30米)的栅格数据作为基本单位，采用常用的“WGS_1984”地理坐标系，文件格式为hgt，分辨率高约30m。整套SRTM1数据文件的每张图幅数据覆盖经纬度各1°，全球数据量极大，总面积超过1.19×10⁸km²。使用以SRTM1为主的地理坐标系栅格数据进行大尺度范围LS因子的提取十分必要。

目前的研究多将SRTM1等地理坐标系栅格数据进行格式转化与投影转换后再利用常用GIS工具(如ArcGIS、ENVI和SAGA等软件)逐步完成 LS因子的计算。如杨勤科等人使用SRTM1和30米分辨率的ASTER GDEM 高程数据在东北样区和黄土样区的相关研究就使用了投影转换的方式进行。基于地理坐标系的SRTM1数据，以经纬度作为基本单位，由于数据投影为球面距离，因此SRTM1的栅格形状近似为梯形。如何直接通过SRTM1数据而不通过投影转换方式进行地形因子相关参数(如坡度、坡长、坡长(L) 因子、坡度(S)因子及地形(LS)因子等)的计算，从而完成大尺度乃至全球LS因子的提取成为研究的难点。

目前LS因子提取算法多基于DEM，主要通过遍历栅格，比较中心栅格与周围栅格的高程值，根据流向算法思想来确定相应的坡度、流向、单元坡长以及汇水面积，从而确定L因子和S因子来计算LS因子。流向算法分单流向和多流向，其中以单流向D8算法应用最为广泛。由于DEM栅格是边长相等的正方形且边长已知，所以计算地形参数相对容易。而大尺度范围LS 因子的提取由于数据量太大，计算十分耗时，有研究针对欧洲尺度水土流失问题使用到SRTM进行相关地形参数的计算，但研究范围也只限于欧洲尺度；也有就全球尺度土壤侵蚀问题使用高分辨率空间分布数据进行了基于 RUSLE模型的LS因子的计算的研究，但其使用的3弧秒分辨率的SRTM数据难以保证提取精度的准确性。当前也缺少通过大尺度SRTM1数据集提取 LS因子的具体计算流程与方式，对大尺度乃至全球LS因子的估算并不快捷高效。

如何通过单流向D8算法完成地理坐标系下LS因子的提取并通过大尺度范围地理坐标系栅格数据实现海量LS因子的提取是一个有待解决的技术问题。

适用于地理坐标系栅格数据的LS因子提取方法设计和实现，有效的避免了对地理坐标系下栅格数据的坐标转换的过程，能提升LS因子提取效率，并且设计的针对大尺度LS因子的计算流程，完善了大尺度LS因子提取的方式。该方法可快速有效地对基于高分辨率全球尺度SRTM1进行LS因子的提取与计算，进而应用于全球土壤侵蚀制图与分析中。

发明内容

本发明的目的是克服现有技术中存在的不足，提供一种适用于大尺度地理坐标系栅格数据的LS因子提取方法，适用于地理坐标系栅格数据的LS因子提取方法设计和实现，有效的避免了对地理坐标系下栅格数据的坐标转换的过程，能提升LS因子提取效率，并且设计的针对大尺度LS因子的计算流程，完善了大尺度LS因子提取的方式；该方法可快速有效地对基于高分辨率全球尺度SRTM1进行LS因子的提取与计算，进而应用于全球土壤侵蚀制图与分析中。

为实现以上技术目的，本发明的技术方案是：一种适用于大尺度地理坐标系栅格数据的LS因子提取方法，包括如下步骤：

步骤1、数据分块：

步骤1.1、将所需尺度地理坐标系栅格各块数据通过Arcmap镶嵌工具融合为一整个大块栅格；

步骤1.2、按实际需要将大块栅格数据通过Arcmap裁剪工具规则划分为小块；

步骤2、添加缓冲区并将栅格数据格式转化：

步骤2.1、通过Arcmap软件中的镶嵌工具，将划分后的每一单块地理坐标系栅格数据，考虑当前数据块周围四个方向，当存在栅格数据时，添加1°缓冲区范围数据进行数据融合；

步骤2.2、通过Arcmap软件中的栅格转文本工具将融合的栅格数据转为 ASCII数据格式；

步骤3、带缓冲区数据的LS因子提取：

步骤3.1、创建日志文件；

步骤3.2、读取ASCII数据头文件和高程及LS因子提取中的参数信息；

步骤3.3、对地理坐标系栅格数据进行无值点和洼地的填充并更新高程值；

步骤3.4、遍历高程二维数组计算坡度、流向和单元坡长；申请存储坡度、流向和单元坡长的矩阵空间，每次遍历高程矩阵时，先判断当前栅格是否为无值点；如果是无值点直接将坡度记为“0”并跳过该点进入下一个点的判断；如果不为无值点则根据D8流向算法思想，按公式(1)方式计算最大坡度值作为该栅格的坡度记录在坡度矩阵中：其中E_c代表中心栅格的高程值，E_i代表当前栅格的高程值，cellsize为当前栅格与中心栅格的距离，其中栅格的南北方向距离记为h_x，栅格的东西方向距离记为h_y，栅格的对角线方向距离记为diagcellsize，h_x和h_y由式(2)(3)得出，diagcellsize根据h_x和 h_y由勾股定理计算得出，其中θ为地理坐标系栅格像元南北宽度；将平地与洼地情况下的坡度设定为0.1，重复上述步骤直至所有有值点坡度计算完成；将坡度最大值所在方向定为该栅格的流向，依据对应方向流向编码将值记录在流向矩阵中，重复上述步骤直至所有有值点流向计算完成；将celsize值作为单元坡长值记录在单元坡长矩阵中；

slope＝max(deg·arctan((E_c-E_i)/cellsize)) (1)

hx＝30.8874791 (2)

hy＝30.8874791·cosθ (3)

步骤3.5、遍历对应二维数组，计算初始汇水面积和初始坡长以及汇水面积:

Step1：计算初始汇水面积，申请存储初始汇水面积的矩阵空间，每次遍历高程矩阵和流向矩阵时，先判断当前栅格是否为无值点；如果无值进行下一个栅格计算，有值则根据流向矩阵和累计栅格个数记录初始化汇水面积；其中栅格的面积等于长与宽的乘积，即h_x·h_y，此面积即初始汇水面积值；重复上述步骤直至所有有值点初始汇水面积赋值完成；Step2：计算初始坡长，申请存储初始坡长的矩阵空间，每次遍历高程矩阵和流向矩阵时，先判断当前栅格是否为无值点；如果无值进行下一个栅格计算，有值则根据流向矩阵和单元坡长矩阵记录初始化坡长；其中单元坡长值为初始坡长值；重复上述步骤直至所有有值点初始坡长赋值完成；Step3：计算汇水面积，申请存储汇水面积的矩阵空间，遍历流向矩阵和初始汇水面积矩阵，将流向当前栅格的汇水面积值的和记为amount；比较amount与当前栅格的汇水面积值，选择较大值作为当前栅格的汇水面积值；正向遍历整个初始化汇水面积数组，计算整个栅格数据的汇水面积值；反向遍历整个初始化汇水面积数组，计算整个栅格数据的汇水面积值；如果正反两边的过程中都没有发生将amount的值赋值给当前栅格的汇水面积值的操作，说明汇水面积提取完成结束循环，否则重头重复继续循环计算；重复上述步骤直至所有点汇水面积赋值完成；

步骤3.6、遍历对应二维数组，根据坡度截断和沟道截断计算累计坡长以及出入口坡长：

Step1：设置坡道截断和沟道截断；申请存储截断值的矩阵空间；每次遍历高程矩阵时，判断当前栅格是否为无值点，并考虑以下截断情况：

如果无值则将该点设置为截断，否则设置成不是截断；遍历坡度矩阵，以坡度的5％作为分界点，小于5％，截断因子设置为0.7；大于等于5％时，截断因子值设置为0.5；当栅格坡度与截断因子的乘积大于流出方向的栅格坡度时，将栅格设置为截断；

遍历汇水面积矩阵，判断栅格的汇水面积值是否大于设定的河网阈值，如果是则将栅格设置为截断，否则不设为截断；

重复上述步骤设置每个栅格的截断；

Step2：计算累计坡长，申请存储累计坡长的矩阵空间，遍历截断矩阵和初始坡长矩阵，先判断当前栅格是否为截断，如果截断则当前栅格的坡长等于初始坡长的一半，如果不截断则当前栅格的初始坡长不变；重复上述步骤设置计算每个栅格加入截断后的初始坡长到初始坡长数组中；再申明临时变量total初始值设置为0；遍历流向矩阵，假设栅格a是当前栅格c周围相邻的一个栅格，而且栅格a的流向指向栅格c，如果栅格a截断则total加上栅格a坡长的一半，如果不截断则total加上栅格a的坡长值；用此方法计算流向当前栅格的坡长值的和记为total；比较total与当前栅格的坡长值，选择较大值作为当前栅格的坡长值；正向遍历整个初始化坡长数组，计算整个栅格数据的坡长值；反向遍历整个初始化坡长数组，计算整个栅格数据的坡长值；如果正反两边的过程中都没有发生将total的值赋值给当前栅格的坡长值的操作，说明累计坡长提取完成结束循环，否则重头重复继续循环计算；重复上述步骤直至所有点累计坡长赋值完成；

Step3：计算出入口坡长，申请存储出口和入口坡长的矩阵空间，每次遍历高程矩阵时，先判断当前栅格是否为无值点；如果无值进行下一个栅格计算，有值则根据流向矩阵、初始坡长矩阵和累计坡长矩阵计算当前栅格的入口坡长和出口坡长；重复上述步骤直至所有有值点出、入口坡长赋值完成；

步骤3.7、遍历对应二维数组，计算S因子与L因子：

Step1：计算S因子，申请存储坡度因子(S)的矩阵空间；每次遍历高程矩阵时，先判断当前栅格是否为无值点；如果无值进行下一个栅格计算，有值则根据坡度矩阵，依据公式(4)计算S因子值；重复上述步骤遍完成每一个栅格的S因子计算；式中θ为坡度；

Step2：计算L因子，申请存储坡长因子(L)的矩阵空间；每次遍历高程矩阵时，先判断当前栅格是否为无值点；如果无值进行下一个栅格计算，有值则根据出入口坡长矩阵，依据公式(5)和(6)计算分段坡L值；其中当入口坡长小于出口坡长值时，使用公式(6)，否则使用公式(6)；重复上述步骤遍完成每一个栅格的L因子计算；式中λ为坡长，m为坡长指数，λout、λin分别为栅格出口及入口的坡长(m)；

L＝(λ/22.13)^m (5)

步骤3.8、遍历对应二维数组，计算LS因子：申请存储坡度坡长因子(LS) 的矩阵空间；每次遍历高程矩阵时，先判断当前栅格是否为无值点；如果无值进行下一个栅格计算，有值则根据L因子和S因子矩阵，依据公式(7) 计算LS因子值；重复上述步骤遍完成每一个栅格的LS因子计算；

LS＝L·S (7)

步骤4、ASCII结果数据格式转化并提取不带缓冲区各LS因子数据：

步骤4.1，通过Arcmap软件中的ASCII转栅格工具将各块带缓冲区的LS因子数据转化回栅格数据；

步骤4.2，通过Arcmap软件中的裁剪工具，使用各未带缓冲区的分块栅格数据裁剪上述所得的LS因子栅格数据；

步骤5、数据融合提取所需尺度的地理坐标系栅格LS因子；将步骤4得到的LS数据通过Arcmap软件中的镶嵌工具融合即可得到所需尺度的地理坐标系下栅格LS因子。

作为优选，步骤1中，针对大尺度范围数据使用分块方式处理数据。

作为优选，在步骤3之前，添加1°缓冲区范围数据并转换栅格格式为 ASCII文本；该方式使用Arcmap软件通过镶嵌和栅格转ASCII完成缓冲区添加和栅格数据转换。

作为优选，步骤3.2中，读取地理坐标系栅格数据的ASCII文件和LS 因子提取中各参数信息的过程为：

Step1：创建一个名为DemData的结构体存放ASCII头信息和LS因子提取中设置的参数信息，申请二维数组保存高程值；

Step2：打开栅格数据文本文件，如果打开失败写入日志并停止执行；

Step3：先按行读取ASCII文件中的内容，文件头部分的格式中以“名称-空格-值”的形式记录；之后将读取的每行数据存到每个字符串中，然后用空格对该字符串进行分割，将得到的值转换成该值的类型并保存到创建的结构体对应的属性中，该过程重复到头文件读取完毕；再将LS因子提取中设置的流向编码、截断因子和河网阈值等参数信息按相同形式记录在结构体对应属性中。

作为优选，步骤3.4中，D8流向算法的栅格编码方式以中心栅格为例，根据不同栅格对应的高程值判断中心栅格周围八个方向的流向，方向自东、东南、南、西南、西、西北、北到东北分别记为1、2、4、8、16、32、64 和128。

作为优选，步骤3.6的Step2中，初始坡长值的更新根据栅格是否截断由单元坡长重新进行赋值：对于基础栅格，初始坡长为原单元坡长值；对于截断栅格，初始坡长为原单元坡长值的一半。

作为优选，步骤5之前，将计算得到的ASCII结果数据转换回栅格数据并去除缓冲区，该方式使用Arcmap软件通过ASCII转栅格和裁剪完成栅格数据转换和无重叠区域LS因子的构建。

从以上描述可以看出，本发明的有益效果在于：

1.针对地理坐标系栅格数据特点，建立以SRTM1为主的地理坐标系栅格模型，基于单流向D8算法思想，直接根据栅格特点逐步计算并提取LS因子，算法提取的LS因子符合实际研究，提取效率较高；

2.针对大尺度范围栅格数据量大的特点，设计大尺度LS因子提取流程，基于“分块计算再合并”的思想，计算每块数据时设置缓冲区以保证计算的准确性，通过ARCGIS等地理信息系统完成较大范围LS因子的提取，流程具有可操作性，为大尺度乃至全球LS因子的提取提供了方案；

3.根据方法实现的最终结果可以看出：将该方法计算得到的LS因子结果与投影坐标系下DEM转化数据得到的LS因子结果统一在同一坐标系统下，99％的LS差值集中在±1之间；通过该算法和流程完成的全球LS因子的提取结果显示，全球LS因子点图过渡平滑且无缝，符合LS因子计算的值域范围，并且欧洲区域计算结果与已有欧洲区域LS因子提取研究结果分布相似且地形特征更加明显，本方法计算结果合理，可实施性强。

附图说明

图1是适用于大尺度地理坐标系栅格(SRTM1为例)的LS因子提取方法流程图；

图2是基于地理坐标系栅格(SRTM1为例)的LS因子提取算法流程图；

图3是SRTM1栅格数据经纬线面图

图4是SRTM1栅格流向与编码方式图

图5是基于DEM数据的县南沟LS因子结果图

图6是基于SRTM1数据的县南沟LS因子结果图

图7是SRTM1与DEM的LS因子差值频率统计图。

具体实施方式

下面结合具体实施例对本发明作进一步说明。

土壤侵蚀已被确定为主要的土壤威胁之一。为应对土壤侵蚀状况复杂的挑战，进行全球土壤侵蚀评估并制作全球土壤侵蚀图来解决土壤侵蚀带来的风险十分必要。地形(LS)因子是常用的土壤侵蚀模型(USLE/RUSLE/CSLE) 中的重要因素。本发明设计并实现了适用于大尺度地理坐标系栅格数据的LS 因子提取方法，为计算全球范围LS因子提供了技术支持和解决方案。

基于GIS的LS因子提取方法，需要先对地理坐标系下栅格进行坐标变换，然后再使用坡度、流向、栅格计算器等常用水文分析工具，才能逐步实现LS 因子的提取。并且在提取大尺度LS因子时，数据量的增多会导致操作起来十分复杂，需要较高的专业性和严谨程度。针对上述现有研究技术中存在的缺陷或不足，设计实现适用于地理坐标系栅格数据的LS因子提取方法及大尺度LS因子提取流程，将极大方便用户的使用，为大尺度乃至全球范围LS因子的提取提供技术支撑与解决方法。

参见图1和图2，以全球范围地理坐标系下栅格数据(SRTM1为主，30 米分辨率ASTERGDEM为空洞地区补充数据的融合数据)为例，适用于大尺度地理坐标系栅格数据的LS因子提取方法的具体流程分为以下步骤：

步骤1：将全球尺度计算的地理坐标系栅格数据进行分块；

使用Arcmap软件中的裁剪工具，将全球拼接的大块SRTM1和ASTER GDEM融合数据按实际需要进行规则划分，划分后的每块数据结构一致，均为grid栅格，全球栅格数据被划分为228个大块。

步骤2：为每个大块数据添加1°的缓冲区并将数据块转化为ASCII文本；

使用Arcmap软件中的融合工具，将划分后的每一单块栅格数据添加其中，并在栅格数据周围四个方向(存在栅格数据时)添加1°缓冲区范围数据用于减少数据边界效应，通过此方式将数据融合为各个较大块的栅格数据；然后使用Arcmap软件中的栅格转文本工具将融合的栅格数据块转为ASCII 数据，各块数据均为含有头部信息和高程值。

步骤3：基于LS因子提取算法完成带缓冲区的各块栅格数据的LS因子提取；

(1)创建日志文件

考虑到算法执行时会出现问题，尤其是在运算栅格数据量大时，在算法运行时间较长的情况下，通过在算法执行时在输出位置建立一个日志文件并创建一个日志文件的类用以输出必要的日志信息。

(2)读取数据

i.输入带缓冲区的栅格融合数据并读取文件头部和LS因子提取中的参数信息：先申请创建结构体用于存放地理坐标系栅格数据的头文件信息；再打开栅格数据文本文件，如果打开失败写入日志并停止执行；再按行读取文件中的内容，文件头部分的格式中以“名称-空格-值”的形式记录；之后将读取的每行数据存到每个字符串中，然后用空格对该字符串进行分割，将得到的值转换成该值的类型并保存到创建的结构体对应的属性中，该过程重复到头文件读取完毕；最后将LS因子提取中设置的流向编码、截断因子和河网阈值等参数信息按相同形式记录在结构体对应属性中。

ii.读取数据高程值：申请创建高程二维数组。由于高程数据中的每一行的每个数据都是以空格间隔的。算法执行过程中按行读取高程值，以字符串形式读取每一行的数据并利用空格进行分割，然后将分割出的每一个数据的字符串转换成float类型并存储到高程数据矩阵中。为了计算方便高程矩阵被扩充了两圈，所以读取数据会将数据置于一个扩大了两圈的矩阵中。

(3)无值点、洼地填充。

i.对高程数据中的无值点和洼地进行填充：无值点是由于测量方式本身的误差导致高程数据中的错误值。洼地在地形数据表现为周围栅格数据均高于中心栅格数据的现象。算法执行时需要将这些栅格填充为周围八个栅格的最小高程值。

(4)遍历高程二维数组计算坡度、流向和单元坡长。

i.申请存储坡度、流向和单元坡长的矩阵空间。每次遍历高程矩阵时，先判断当前栅格是否为无值点。如果是无值点直接将坡度记为“0”并跳过该点进入下一个点的判断；如果不为无值点则按以下方式计算栅格单元边长 cellsize：

由公式(2)(3)得到的结果设置地理坐标系栅格数据经纬线方向上的栅格宽度(距离)h_x与h_y，对于h_y有公式(8)：

cosθ＝cos((yllcorner+((nrows-1)-(i-2))·cellsize)/deg (8)

其中，yllcorner是从栅格头文件信息中获取的表示该数据所处空间的坐标信息；nrows是在栅格数据原有的行数基础上扩充了两圈后的行数；i表示当前遍历的栅格相对于起始栅格在经度方向的偏移量。对于栅格的对角线方向的像元距离diagcellsize，根据勾股定理由公式(9)计算得出：

ii.根据D8流向算法思想计算坡度、流向及单元坡长：

①：依次将当前栅格与周围的8个栅格的高程值进行对比，计算周围栅格与中心栅与中心栅格所成角度，计算方式如公式(10)所示：

angle＝deg·arctan((E_c-E_i)/cellsize) (10)

其中E_c代表中心栅格的高程值，E_i代表当前栅格的高程值。cellsize为当前栅格与中心栅格的距离，其中栅格东西方向为h_y，栅格南北方向为h_x，栅格对角线方向上为diagcellsize。

②：选取周围角度最大值作为中心栅格的坡度值，计算方式如公式(11)所示：

slope＝max(deg·arctan((E_c-E_i)/cellsize)) (11)

③：选取周围角度最大值所在方向为该中心栅格的流向，依据流向编码对其进行赋值。

④：根据栅格流向，将栅格单元边长值celsize作为单元坡长值进行赋值。

iii.平地与洼地情况下坡度的设定：由于平地、洼地无流向，本算法对此作简化处理，为保证栅格与河网链接，将其坡度设置为最小值0.1。重复以上步骤直至所有栅格遍历完成，并将计算完成的坡度结果输出到坡度数组中。(5)遍历对应二维数组，计算初始汇水面积和初始坡长以及汇水面积。

i.初始汇水面积提取：申请存储初始汇水面积的矩阵空间。每次遍历高程矩阵和流向矩阵时，先判断当前栅格是否为无值点。如果无值进行下一个栅格计算，有值则根据流向矩阵和累计栅格个数记录初始化汇水面积。其中栅格的面积等于长与宽的乘积，即h_x·h_y，由公式(2)(3)(8)可得到此面积即初始汇水面积值。重复上述步骤直至所有有值点初始汇水面积赋值完成，并将计算完成的初始汇水面积值输出到初始汇水面积数组中；

ii.初始坡长提取：申请存储初始坡长的矩阵空间。每次遍历高程矩阵和流向矩阵时，先判断当前栅格是否为无值点。如果无值进行下一个栅格计算，有值则根据流向矩阵和单元坡长矩阵记录初始化坡长。将单元坡长值赋值给初始坡长。重复上述步骤直至所有有值点初始坡长赋值完成，并将计算完成的初始坡长输出到初始坡长数组中。

iii.汇水面积提取：申请存储汇水面积的矩阵空间，遍历流向矩阵和初始汇水面积矩阵，计算流向当前栅格的汇水面积值的和记为amount；比较 amount与当前栅格的汇水面积值，选择较大值作为当前栅格的汇水面积值；正向遍历整个初始化汇水面积数组，计算整个栅格数据的汇水面积值；反向遍历整个初始化汇水面积数组，计算整个栅格数据的汇水面积值；如果正反两边的过程中都没有发生将amount的值赋值给当前栅格的汇水面积值的操作，说明汇水面积提取完成结束循环，否则重头重复继续循环计算。重复上述步骤直至所有点汇水面积赋值完成，最后将计算完成的汇水面积输出到汇水面积数组中。

(6)遍历对应二维数组，根据坡度截断和沟道截断计算累计坡长以及出入口坡长。

i.申请存储截断值的矩阵空间。每次遍历高程矩阵时，判断当前栅格是否为无值点。如果无值则将该点设置为截断，否则设置成不是截断；遍历坡度矩阵，判断栅格的坡度是否小于5％(约为2.861°)，由于此时截断因子为 0.7，当栅格坡度与截断因子的乘积大于流出方向的栅格坡度时，将栅格设置为截断；如果坡度大于等于5％(约为2.861°)，此时截断因子为0.5，当栅格坡度与截断因子的乘积大于流出方向的栅格坡度时，将栅格也设置为截断；再遍历汇水面积矩阵，判断栅格的汇水面积值是否大于设定的河网阈值，如果是则将栅格设置为截断，否则不设为截断；重复上述步骤设置每个栅格的截断并将计算完成的截断设置输出到截断数组中。

ii.累计坡长提取：申请存储累计坡长的矩阵空间，遍历截断矩阵和初始坡长矩阵，先判断当前栅格是否为截断，如果截断则当前栅格的坡长等于初始坡长的一半，如果不截断则当前栅格的初始坡长不变；重复上述步骤设置计算每个栅格加入截断后的初始坡长到初始坡长数组中。

再申明临时变量total初始值设置为0；遍历流向矩阵，假设栅格a是当前栅格c周围相邻的一个栅格，而且栅格a的流向指向栅格c，如果栅格a 截断则total加上栅格a坡长的一半，如果不截断则total加上栅格a的坡长值。用此方式计算流向当前栅格的坡长值的和记为total；比较total与当前栅格的坡长值，选择较大值作为当前栅格的坡长值；正向遍历整个初始化坡长数组，计算整个栅格数据的坡长值；反向遍历整个初始化坡长数组，计算整个栅格数据的坡长值；如果正反两边的过程中都没有发生将total的值赋值给当前栅格的坡长值的操作，说明累计坡长提取完成结束循环，否则重头重复继续循环计算。重复上述步骤直至所有点累计坡长赋值完成，最后将计算完成的累计坡长输出到累计坡长数组中；

iii.出、入口坡长提取：申请存储出口和入口坡长的矩阵空间。每次遍历高程矩阵时，先判断当前栅格是否为无值点。如果无值进行下一个栅格计算，有值则根据流向矩阵、初始坡长矩阵和累计坡长矩阵计算当前栅格的入口坡长和出口坡长；重复上述步骤直至所有有值点出、入口坡长赋值完成，并将计算完成的出入口坡长输出到出、入口坡长数组中。

(7)遍历对应二维数组，计算S因子与L因子。

i.S因子提取：申请存储坡度因子(S)的矩阵空间。每次遍历高程矩阵时，先判断当前栅格是否为无值点。如果无值进行下一个栅格计算，有值则根据坡度矩阵，依据公式(4)计算S因子值；重复上述步骤遍完成每一个栅格的S因子计算；

ii.L因子提取：申请存储坡长因子(L)的矩阵空间。每次遍历高程矩阵时，先判断当前栅格是否为无值点。如果无值进行下一个栅格计算，有值则根据出入口坡长矩阵，依据公式(5)和(6)计算分段坡L值。其中当入口坡长小于出口坡长值时，使用公式(6)，否则使用公式(5)；重复上述步骤遍完成每一个栅格的L因子计算。

(8)遍历对应二维数组，计算LS因子。

申请存储坡度坡长因子(LS)的矩阵空间。每次遍历高程矩阵时，先判断当前栅格是否为无值点。如果无值进行下一个栅格计算，有值则遍历L因子和S因子矩阵，根据公式(7)计算LS因子值；重复上述步骤遍完成每一个栅格的LS因子提取。

步骤4：将LS因子提取算法得到的ASCII形式的LS因子数据转化为栅格数据，并使用不带缓冲区的各块地形数据对已得到的LS因子数据块进行裁剪，提取不带缓冲区的各LS因子块数据；

用本发明的LS因子提取算法完成各局部块的计算结果为ASCII数据，通过使用Arcmap软件中的ASCII转栅格工具将各块LS因子数据转化为栅格数据。再使用Arcmap软件中的裁剪工具，通过各分块栅格数据裁剪上述所得的LS因子栅格数据块，从而得到不带缓冲区(各块之间无重叠)的局部块LS因子栅格数据。

步骤5：将步骤四得到的LS因子栅格数据融合，最终得到1弧秒分辨率的全球LS因子。

使用Arcmap软件中的融合工具，将步骤4得到的各局部块LS因子进行融合操作，即可得到全球1弧秒分辨率分辨率的无缝LS因子栅格数据。

实验部分：

将本地理坐标系栅格数据LS因子提取算法和传统投影坐标系栅格数据 LS因子提取算法进行对比，并使用该计算流程与算法完成全球LS因子提取，对比已有欧洲区域研究结果：

实验背景：

地理坐标系下栅格数据得到的计算结果的正确性与可实施性验证，最直接有效的测试方式就是将其计算结果和传统投影坐标系下的LS因子提取算法对比分析，并通过实施该方法得到的大尺度乃至全球LS因子结果，参照与已有研究分析其合理性。

实验区域：

陕北绥德县南沟流域的SRTM1与30米分辨率DEM数据、全球地理坐标系栅格数据(SRTM1为主，30米分辨率ASTERGDEM为空洞地区补充数据的融合数据)。

实验方法：

1.将1弧秒分辨率县南沟SRTM1数据通过上述LS因子提取算法计算得到的结果进行投影转换，与将30m分辨率县南沟DEM数据通过传统投影坐标系下LS因子提取算法计算得到的结果置于同一坐标系下，对比二者结果的差异。图5是基于DEM的投影坐标系下LS因子提取算法的县南沟LS值，图6是基于SRTM1的地理坐标系下LS因子提取算法的县南沟LS因子值。图7是对SRTM1和DEM的LS因子结果作差扩大100倍得到的县南沟LS 因子差值频率统计图；

2.将全球地理坐标系栅格数据通过上述流程和算法计算来完成1弧秒分辨率全球LS因子提取，查看其分布和值域范围是否合理，并对比已有欧洲区域LS因子数据看结果分布。

结果分析：

DEM数据为30米的数字高程模型数字地图，SRTM1为1弧秒分辨率的遥感高程地图(即30米分辨率SRTM数据)。由图4和图5可得：基于两种数据得到的LS因子的整体空间分布图和LS因子值的范围十分相近，由两种不同的数据地图提取出的LS因子的趋势是一致的。由于LS因子是由S因子和L因子计算得到的，其中S因子是由坡度计算得到，L因子是由坡长计算得到，而LS因子主要受坡度影响较大，因此LS因子的计算值会受到一定影响。

STRM1数据值是栅格在经度上的跨越距离，通过地球的半径和几何公式算的以1弧秒为单位的SRTM数据在经度上的跨越距离约等于30米。由于 SRTM1数据栅格在经度上的跨越长度固定，其在纬度上的跨越长度随纬度的增大而减小，因此不在赤道位置上的SRTM1数据的栅格大小略小于30m分辨率的DEM栅格。LS因子结果表明，由于DEM数据与STRM1数据在相同位置的高程差相等，当栅格的流向不南北方向时，SRTM1的栅格间距会小于DEM栅格，所以由SRTM1计算得出的LS因子值会大于DEM计算得出的LS因子值；当栅格流向为南北方向时，SRTM1的栅格间距会略大于DEM 栅格间距，所以STRM1提取的LS因子值有时会大于DEM提取的LS因子值。从LS因子差值频率统计图可知：99％的差值集中在±1之间，结果与传统DEM方法结果保持高度一致。从全球LS因子提取结果图可知：LS因子点图过渡平滑且无缝，符合LS因子计算的值域范围，值分布合理。与已有欧洲区域提取结果相比，本方法LS因子结果地形分布特征更加明显，在LS因子较大值和较小值区域分层现象更为突出，易于发现局部LS因子高低点。虽然部分区域LS因子值分布有偏差，但这是基于不同的LS因子计算公式(实际坡度及其对应LS因子计算公式的参数设置不同)所致，误差在合理范围内。

县南沟区域SRTM1计算的LS因子结果与DEM计算的结果略有不同，但两种数据产生的结果高度相似验证了本LS因子提取算法的准确性。而基于全球地理坐标系栅格数据得到的全球LS因子数据分布比较“平滑”，值域范围正确，空间分布合理，与已有欧洲区域LS因子结果分布相似且地形特征更加明显。形成的1弧秒分辨率全球无缝隙地图能有效解决全球范围LS 因子点图分辨率过低的缺点，也能较好地反映了地形对坡面侵蚀的影响。该方法能够快速、高效、准确地提取大尺度乃至全球LS因子。

以上为本发明较佳的实施方式，本发明所属领域的技术人员还能够对上述实施方式进行变更和修改，因此，本发明并不局限于上述的具体实施方式，凡是本领域技术人员在本发明的基础上所作的任何显而易见的改进、替换或变型均属于本发明的保护范围。

Claims (7)

1.一种适用于大尺度地理坐标系栅格数据的LS因子提取方法，其特征在于，包括如下步骤：

步骤1、数据分块：

步骤1.1、将所需尺度地理坐标系栅格各块数据通过Arcmap镶嵌工具融合为一整个大块栅格；

步骤1.2、按实际需要将大块栅格数据通过Arcmap裁剪工具规则划分为小块；

步骤2、添加缓冲区并将栅格数据格式转化：

步骤2.1、通过Arcmap软件中的镶嵌工具，将划分后的每一单块地理坐标系栅格数据，考虑当前数据块周围四个方向，当存在栅格数据时，添加1°缓冲区范围数据进行数据融合；

步骤2.2、通过Arcmap软件中的栅格转文本工具将融合的栅格数据转为ASCII数据格式；

步骤3、带缓冲区数据的LS因子提取：

步骤3.1、创建日志文件；

步骤3.2、读取ASCII数据头文件和高程及LS因子提取中的参数信息；

步骤3.3、对地理坐标系栅格数据进行无值点和洼地的填充并更新高程值；

步骤3.4、遍历高程二维数组计算坡度、流向和单元坡长；申请存储坡度、流向和单元坡长的矩阵空间，每次遍历高程矩阵时，先判断当前栅格是否为无值点；如果是无值点直接将坡度记为“0”并跳过该点进入下一个点的判断；如果不为无值点则根据D8流向算法思想，按公式(1)方式计算最大坡度值作为该栅格的坡度记录在坡度矩阵中：其中E_c代表中心栅格的高程值，E_i代表当前栅格的高程值，cellsize为当前栅格与中心栅格的距离，其中栅格的南北方向距离记为h_x，栅格的东西方向距离记为h_y，栅格的对角线方向距离记为diagcellsize，h_x和h_y由式(2)(3)得出，diagcellsize根据h_x和h_y由勾股定理计算得出，其中θ为地理坐标系栅格像元南北宽度；将平地与洼地情况下的坡度设定为0.1，重复上述步骤直至所有有值点坡度计算完成；将坡度最大值所在方向定为该栅格的流向，依据对应方向流向编码将值记录在流向矩阵中，重复上述步骤直至所有有值点流向计算完成；将cellsize值作为单元坡长值记录在单元坡长矩阵中；

slope＝max(deg·arctan((E_c-E_i)/cellsize)) (1)

hx＝30.8874791 (2)

hy＝30.8874791·cosθ (3)

步骤3.5、遍历对应二维数组，计算初始汇水面积和初始坡长以及汇水面积:

Step1：计算初始汇水面积，申请存储初始汇水面积的矩阵空间，每次遍历高程矩阵和流向矩阵时，先判断当前栅格是否为无值点；如果无值进行下一个栅格计算，有值则根据流向矩阵和累计栅格个数记录初始化汇水面积；其中栅格的面积等于长与宽的乘积，即h_x·h_y，此面积即初始汇水面积值；重复上述步骤直至所有有值点初始汇水面积赋值完成；Step2：计算初始坡长，申请存储初始坡长的矩阵空间，每次遍历高程矩阵和流向矩阵时，先判断当前栅格是否为无值点；如果无值进行下一个栅格计算，有值则根据流向矩阵和单元坡长矩阵记录初始化坡长；其中单元坡长值为初始坡长值；重复上述步骤直至所有有值点初始坡长赋值完成；Step3：计算汇水面积，申请存储汇水面积的矩阵空间，遍历流向矩阵和初始汇水面积矩阵，将流向当前栅格的汇水面积值的和记为amount；比较amount与当前栅格的汇水面积值，选择较大值作为当前栅格的汇水面积值；正向遍历整个初始化汇水面积数组，计算整个栅格数据的汇水面积值；反向遍历整个初始化汇水面积数组，计算整个栅格数据的汇水面积值；如果正反两边的过程中都没有发生将amount的值赋值给当前栅格的汇水面积值的操作，说明汇水面积提取完成结束循环，否则重头重复继续循环计算；重复上述步骤直至所有点汇水面积赋值完成；

步骤3.6、遍历对应二维数组，根据坡度截断和沟道截断计算累计坡长以及出入口坡长：

Step1：设置坡道截断和沟道截断；申请存储截断值的矩阵空间；每次遍历高程矩阵时，判断当前栅格是否为无值点，并考虑以下截断情况：

如果无值则将该点设置为截断，否则设置成不是截断；遍历坡度矩阵，以坡度的5％作为分界点，小于5％，截断因子设置为0.7；大于等于5％时，截断因子值设置为0.5；当栅格坡度与截断因子的乘积大于流出方向的栅格坡度时，将栅格设置为截断；

遍历汇水面积矩阵，判断栅格的汇水面积值是否大于设定的河网阈值，如果是则将栅格设置为截断，否则不设为截断；

重复上述步骤设置每个栅格的截断；

Step2：计算累计坡长，申请存储累计坡长的矩阵空间，遍历截断矩阵和初始坡长矩阵，先判断当前栅格是否为截断，如果截断则当前栅格的坡长等于初始坡长的一半，如果不截断则当前栅格的初始坡长不变；重复上述步骤设置计算每个栅格加入截断后的初始坡长到初始坡长数组中；再申明临时变量total初始值设置为0；遍历流向矩阵，假设栅格a是当前栅格c周围相邻的一个栅格，而且栅格a的流向指向栅格c，如果栅格a截断则total加上栅格a坡长的一半，如果不截断则total加上栅格a的坡长值；用此方法计算流向当前栅格的坡长值的和记为total；比较total与当前栅格的坡长值，选择较大值作为当前栅格的坡长值；正向遍历整个初始化坡长数组，计算整个栅格数据的坡长值；反向遍历整个初始化坡长数组，计算整个栅格数据的坡长值；如果正反两边的过程中都没有发生将total的值赋值给当前栅格的坡长值的操作，说明累计坡长提取完成结束循环，否则重头重复继续循环计算；重复上述步骤直至所有点累计坡长赋值完成；

Step3：计算出入口坡长，申请存储出口和入口坡长的矩阵空间，每次遍历高程矩阵时，先判断当前栅格是否为无值点；如果无值进行下一个栅格计算，有值则根据流向矩阵、初始坡长矩阵和累计坡长矩阵计算当前栅格的入口坡长和出口坡长；重复上述步骤直至所有有值点出、入口坡长赋值完成；步骤3.7、遍历对应二维数组，计算S因子与L因子：

Step1：计算S因子，申请存储坡度因子(S)的矩阵空间；每次遍历高程矩阵时，先判断当前栅格是否为无值点；如果无值进行下一个栅格计算，有值则根据坡度矩阵，依据公式(4)计算S因子值；重复上述步骤遍完成每一个栅格的S因子计算；式中θ为坡度；

Step2：计算L因子，申请存储坡长因子(L)的矩阵空间；每次遍历高程矩阵时，先判断当前栅格是否为无值点；如果无值进行下一个栅格计算，有值则根据出入口坡长矩阵，依据公式(5)和(6)计算分段坡L值；其中当入口坡长小于出口坡长值时，使用公式(6)，否则使用公式(6)；重复上述步骤遍完成每一个栅格的L因子计算；式中λ为坡长，m为坡长指数，λout、λin分别为栅格出口及入口的坡长(m)；

L＝(λ/22.13)^m (5)

步骤3.8、遍历对应二维数组，计算LS因子：申请存储坡度坡长因子(LS)的矩阵空间；每次遍历高程矩阵时，先判断当前栅格是否为无值点；如果无值进行下一个栅格计算，有值则根据L因子和S因子矩阵，依据公式(7)计算LS因子值；重复上述步骤遍完成每一个栅格的LS因子计算；

LS＝L·S (7)

步骤4、ASCII结果数据格式转化并提取不带缓冲区各LS因子数据：

步骤4.1，通过Arcmap软件中的ASCII转栅格工具将各块带缓冲区的LS因子数据转化回栅格数据；

步骤4.2，通过Arcmap软件中的裁剪工具，使用各未带缓冲区的分块栅格数据裁剪上述所得的LS因子栅格数据；

步骤5、数据融合提取所需尺度的地理坐标系栅格LS因子；将步骤4得到的LS数据通过Arcmap软件中的镶嵌工具融合即可得到所需尺度的地理坐标系下栅格LS因子。

2.如权利要求1所述的适用于大尺度地理坐标系栅格数据的LS因子提取方法，其特征在于，步骤1中，针对大尺度范围数据使用分块方式处理数据。

3.如权利要求1所述适用于大尺度地理坐标系栅格数据的LS因子提取方法，其特征在于，在步骤3之前，添加1°缓冲区范围数据并转换栅格格式为ASCII文本；该方式使用Arcmap软件通过镶嵌和栅格转ASCII完成缓冲区添加和栅格数据转换。

4.如权利要求1所述的适用于大尺度地理坐标系栅格数据的LS因子提取方法，其特征在于，步骤3.2中，读取地理坐标系栅格数据的ASCII文件和LS因子提取中各参数信息的过程为：

Step1：创建一个名为DemData的结构体存放ASCII头信息和LS因子提取中设置的参数信息，申请二维数组保存高程值；

Step2：打开栅格数据文本文件，如果打开失败写入日志并停止执行；

Step3：先按行读取ASCII文件中的内容，文件头部分的格式中以“名称-空格-值”的形式记录；之后将读取的每行数据存到每个字符串中，然后用空格对该字符串进行分割，将得到的值转换成该值的类型并保存到创建的结构体对应的属性中，该过程重复到头文件读取完毕；再将LS因子提取中设置的流向编码、截断因子和河网阈值等参数信息按相同形式记录在结构体对应属性中。

5.如权利要求1所述的适用于大尺度地理坐标系栅格数据的LS因子提取方法，其特征在于，步骤3.4中，D8流向算法的栅格编码方式以中心栅格为例，根据不同栅格对应的高程值判断中心栅格周围八个方向的流向，方向自东、东南、南、西南、西、西北、北到东北分别记为1、2、4、8、16、32、64和128。

6.如权利要求1所述的适用于大尺度地理坐标系栅格数据的LS因子提取方法，其特征在于，步骤3.6的Step2中，初始坡长值的更新根据栅格是否截断由单元坡长重新进行赋值：对于基础栅格，初始坡长为原单元坡长值；对于截断栅格，初始坡长为原单元坡长值的一半。

7.如权利要求1所述的适用于大尺度地理坐标系栅格数据的LS因子提取方法，，其特征在于，步骤5之前，将计算得到的ASCII结果数据转换回栅格数据并去除缓冲区，该方式使用Arcmap软件通过ASCII转栅格和裁剪完成栅格数据转换和无重叠区域LS因子的构建。

Images