CN111177917A - 一种基于srtm的坡长提取方法 - Google Patents

Abstract

本发明提出了一种基于SRTM的坡长提取方法，通过建立地理坐标下栅格模型，从该模型中直接推导出SRTM栅格边长的计算方式。计算坡长时考虑坡度截断、河网截断影响，使结果尽可能符合实际情况，最后使用正反遍历计算SRTM的累计坡长。与传统方法相比，本方法不再经过投影坐标变换，从而提高坡长的提取效率；同时也能保证原始数据的精度不下降，坡长提取结果相对精确，坡长结果与现有的DEM计算结果对比显示，坡长差值范围集中，具有规律性。

Description

一种基于SRTM的坡长提取方法

技术领域

本发明涉及侵蚀学、地理学、计算机科学交叉的数字地形分析技术，具体涉及一种基于SRTM的坡长提取方法，设计并且实现基于SRTM地理坐标下坡长的提取方法。

背景技术

坡长是水土保持、土壤侵蚀等研究的重要影响因子。在USLE中将坡长定义为坡面漫流的起点到径流被截断或者流路中断的水平距离。在以下两种情况下被终止(截断)：(1)坡度变缓以致发生沉积时，称为坡度截断(Slope Cutoff Factor，SCF)；(2)明显的河网处，称为河网截断因子(Chanel Networks Cutoff Factor CNCF)。在坡面尺度，通常坡度和坡长可实际测量得到。在流域或者区域土壤侵蚀评价中，坡长通常通过DEM(DigitalElevation Model)进行提取。

SRTM(shuttle radar topography mission,SRTM)数据是指美国“奋进号”航天飞机上的干涉成像雷达系统在全球性作业中,获得了地球纬度范围在60°N至56°S间，面积超过1.19×108km,覆盖陆地表面80％面积的三维雷达数据。其整套数据是以geotif作为基本格式、以1弧度栅格数据作为基本单位存储的。随着侵蚀研究尺度的不断增大，大洲尺度，甚至全球范围的侵蚀估算成为研究的必须。通常使用SRTM数据进行提取，但是基于该数据的坡长直接提取方法尚不成熟。

目前基于SRTM的坡长提取方法主要通过格式转换和投影变换到平面坐标，再以DEM数据的提取方法进行坡长提取，提取坡长的方法相当繁杂，提取效率非常低；在进行坐标变换过程中会使原始数据的精度降低，提取结果误差相对也会增大。

发明内容

针对现有技术存在的不足，本发明的目的在于，提供一种基于SRTM的坡长提取方法，解决现有技术中不能用SRTM数据直接提取坡长的技术问题。

为了解决上述技术问题，本发明采用如下技术方案予以实现：

一种基于SRTM的坡长提取方法，该方法按照以下步骤进行：

步骤1：存储SRTM数据的头信息，同时获取SRTM数据中每一个栅格的坡度值和流向；

步骤2：计算栅格单元坡长值；

步骤2.1：遍历所有栅格，判断当前栅格是否是无值点：若是无值点：则设置该栅格的单元坡长值为0，进入下一个栅格的判断；若是非无值点，执行步骤2.2；

步骤2.2：根据步骤1获取的栅格流向，若流向为东西方向，该当前栅格的单元坡长为LONG；若流向为南北方向，该当前栅格的单元坡长为LAT；若流向为东南、西南、东北、西北中的任意一种，则该当前栅格的单元坡长为diagCellSize；

其中，LAT＝2π×6371000·cellsize/360；

LONG＝LAT·cosθ，θ为当前栅格中任意一点与地球球心连线与赤道的夹角；

步骤2.3：遍历所有栅格，得到每一个栅格的单元坡长值；

步骤3：计算栅格初始坡长值；

步骤3.1：判断当前栅格是否截断，若为不截断，则当前格栅初始坡长值等于单元坡长值；若为截断，则当前栅格的初始坡长值为单元坡长值的一半；

步骤3.2：遍历所有栅格，得到每一个栅格的初始坡长值；

步骤4：计算栅格累积坡长：

步骤4.1：将当前栅格作为初始累积栅格，初始累积栅格初始值设置为0，任意与初始累积栅格相邻的栅格中流向指向初始累积栅格的定义为相邻指向栅格，若其中任意一个相邻指向栅格截断，则取初始值与该相邻指向栅格的初始坡长值的一半之和作为初始累积栅格对该相邻指向栅格的累积坡长；若该相邻指向栅格不截断，则取初始值与该相邻指向栅格的初始坡长值之和作为初始累积栅格对该相邻指向栅格的累积坡长；

依次计算初始累积栅格对其他相邻指向栅格的累积坡长，比较初始累积栅格与所有相邻指向栅格的累积坡长的总和SUM与初始累积栅格的初始坡长值的大小，若SUM＞初始累积栅格的初始坡长值，则将SUM赋值给初始累积栅格作为初始累积栅格的初始坡长值，否则初始累积栅格的初始坡长值不变；

步骤4.2：按照步骤4.1的方法分别正向遍历所有栅格和反向遍历所有栅格，计算每个栅格的SUM并判断是否赋值，至不再出现将SUM赋值给当前栅格的初始坡长值的情况，遍历结束，最后一次赋值时获得的当前栅格的初始坡长值即为提取的该栅格的累积坡长。

具体的，步骤1中，存储SRTM数据的头信息，同时获取SRTM数据中每一个栅格的坡度值和流向，包括：

步骤1.1：创建结构体DemData，用于存放SRTM的头信息cellsize；

步骤1.2：按行读取SRTM的数据，SRTM的头信息cellsize的格式为：前边是名称，接着是空格，后边是头信息cellsize的值；将读取的每行数据存到对应的一个字符串中，然后对该字符串进行空格分割，最后将头信息cellsize的值转换成该值的类型并保存到数据结构DemData对应的属性中；

步骤1.3：采用最陡坡降法计算每个单元栅格的坡度值和流向。

优选的，在步骤2之前，对步骤1的结果进行无值点和洼地填充。

具体的，步骤3中：判断当前栅格是否截断，包括以下方法中的任意一种：

a：若当前栅格是无值点，则当前栅格为截断，否则不是截断；

b：若当前栅格的坡度值小于2.861°，且坡度变化率大于0.7，则当前栅格为截断，否则不是截断；

c：若当前栅格的坡度值大于等于2.861°，且坡度变化率大于0.5，则当前栅格为截断，否则不是截断；

d：若当前栅格的河网值大于设定的河网截断阈值，则当前栅格为截断，否则不是截断。

更具体的，d中，当前栅格的河网值的计算包括：

第一步：将每个栅格的面积S作为该栅格的初始化河网值，其中，S＝LONG×LAT，遍历所有栅格，得到所有栅格的初始化河网值；

第二步：当前栅格的河网值提取：

Step1：计算当前栅格中流向指向当前栅格的相邻栅格的初始化河网值的和，记作sum；

Step2：比较sum与当前栅格的初始化河网值，如果sum＞当前栅格的初始化河网值，则将sum赋值给当前栅格，否则当前栅格的河网值不变；

Step3：正向遍历所有栅格，计算每个栅格的sum并判断是否赋值；

Step4：反向遍历所有栅格，计算每个栅格的sum并判断是否赋值；

Step5：若正反遍历的过程中不再发生赋值，遍历结束，最后一次赋值作为当前栅格的河网值。

本发明与现有技术相比，具有如下技术效果：

本发明基于SRTM的坡长提取方法，有效地避免了SRTM数据格式转换及坐标变换过程，提高提取效率，为从区域地形因子提取到大范围、大尺度土壤侵蚀模型应用起到桥梁作用，为大洲尺度甚至全球范围的地形因子快速提取提供了技术支持。具体来说：

1、本发明的方法建立地理坐标系下的栅格模型，直接对原始的SRTM数据进行分析，计算出坡长值，避免了投影变换过程中导致的原始数据精度损失问题，同时提高了坡长的提取效率。

2、从算法的思想上看，本算法根据栅格所在经纬度，先计算SRTM栅格的边长，使用正反遍历计算累计坡长，同时考虑坡度截断、河网截断，以保证计算的结果符合实际情形。

3、从算法实现的最终结果来看：将计算得到的坡长结果进行投影转换，与平面坐标下的DEM数据得到的坡长结果误差较小，本算法计算结果可接受。

附图说明

图1是基于SRTM的坡长提取流程图。

图2是SRTM栅格模型图。

图3是SRTM栅格切面示意图。

图4水流矩阵示意图。

图5是当前计算栅格流向编码图。

图6是基于DEM的县南沟坡长结果。

图7是基于SRTM的县南沟坡长结果。

图8是SRTM与DEM坡长差值曲线图。

具体实施方式

以下对本发明中涉及到的技术术语进行解释：

其中，DEM数据基于平面坐标，其栅格基于地形的实际长度，一般为米、千米等。而SRTM数据基于地理坐标系，栅格以弧度作为基本单位存储，栅格不规则。SRTM数据，包括头信息cellsize，其格式为：前边是名称，接着是空格，后边是头信息cellsize的值。

此外，本发明所指栅格数据，是以二维矩阵(行和列或格网)的形式来表示空间地物或现象分布的数据组织方式，每个矩阵单位称为一个栅格单元(cell),每个像元都包含一个信息值。栅格数据模型包括数字高程数据、卫星影像、数字正射影像、扫描地图和图形文件等。

本发明所述的流向，是指SRTM中各栅格的流向，应用最多的是D8法，假设单个栅格中的水流只有8中可能的流向，这8个栅格与之相邻，8个栅格相对于该栅格的方位，形成了该栅格的多种可能的流向。其中，以东西方向为例，指的是：该栅格水流流向左右栅格的方向。

本发明的方法，中心思想在于：首先建立地理坐标系下的栅格模型，图2为模拟的地理坐标模型，O为球心，AO为地球半径，以BC为轴的水平方向切面为赤道面，而垂直切面是一个纬线面。如图3(1)，D表示地球表面上的一点，该点与球心O连线与赤道的夹角，弧DC表示夹角为θ时，对应的球面距离，代表该点的纬度，r为这一点纬线圈半径；图3(2)为D的经线圈所在平面，弧DE表示该纬线圈上的一段距离，ω表示这段距离对应的经度差。

用R表示地球的平均半径，则有公式1和公式2：

LAT_DC＝2πR·θ/360 (1)

LONG_DE＝2πr·ω/360 (2)

由r＝R·cosθ有公式3：

LONG_DE＝2πR·cosθ·ω/360 (3)

通过该模型，将SRTM栅格数据的坐标进行一般话计算，若用LONG表示经线方向的像元栅格宽度(距离)，用LAT表示纬线方向的栅格宽度(距离)。栅格南北宽度(度)用A表示，东西宽度(度)B表示。一般在SRTM数据中栅格在经纬度上的跨度相同，即A＝B。

在取R＝6371000m，取β＝1/3600则有公式4和公式5：

LAT＝30.8874791 (4)

LONG＝30.8874791·cosθ (5)

基于单流向D8思想，通过最陡坡降法计算每个单元栅格的坡度和流向，图4为水流矩阵示意图。水流方向规定了每一个栅格的流动路径，而单元坡长则代表一个栅格在流向方向上的长度。坡长计算使用以下方式确定：第一、通过遍历流向矩阵，确定每一个栅格的单元坡长；第二、通过计算截断，让截断点的单元坡长为原来的一半，且不在流路上累加；第三、由于无法确定沿着哪条流路得到的坡长最大，对截断校正过的单元坡长矩阵使用正反遍历方法得到最后的累计坡长。

实施例1：

本实施例提供一种基于SRTM的坡长提取方法，具体流程分为以下步骤：

步骤1：存储SRTM数据的头信息cellsize的值，同时获取SRTM数据中每一个栅格的坡度值和流向。具体包括：

1.1数据输入

(1)读取数据头部信息：创建结构体,按行读取SRTM数据的所有头部信息，存放到结构体对应的属性中。

(2)读取SRTM高程数据。将转换成float类型的高程数据存到高程数据矩阵中，为了计算的方便同时保证每个数据的周围栅格的一致性，将数据扩充两圈，所以数据放进一个扩大了两圈的矩阵中用于计算。

1.2：无值点、洼地填充。

1.3：遍历二维数组，基于单流向思想，计算栅格的坡度、流向。

步骤2：根据流向计算栅格单元坡长值。

(1)在每次遍历时，先判断当前的栅格是否是无值点。如果是无值点：直接设置该栅格的单元坡长值为0并且跳过，进入下一个点的判断。

(2)如果非无值点，则根据栅格的流向确定单元坡长。

设置地球的半径值R＝6371000m。通过获取到的头信息中cellsize这一项的值，计算栅格在经度方向的宽度LAT。计算公式是：

LAT＝2π×6371000·cellsize/360

在纬线方向上的栅格宽度(距离)为LONG，以及：LONG＝2πR·cosθ·ω/360，β为从数据文本的头信息中获取，值为文本数据cellsize这一项的值。即:

LONG＝2π×6371000·cosθ·cellsize/360＝LAT·cosθ

栅格的对角线方向的像元距离diagCellSize，采用勾股定理得：

Step3:单元坡长的确定。如果该点的流向为东西方向，该点的单元坡长为LONG；流向为南北方向，该点的单元坡长为LAT；流向为东南、西南、东北、西北，该点的单元坡长为diagCellSize。

(4)对于平地和洼地坡长值的计算：本发简化了平地流向处理，对于那些由于平地导致不能流出的栅格单独处理。具体做法是：遍历这些栅格的周围八个点，如果周围八个点栅格的高程值大于等于中心点栅格，则将其作为该点栅格的流出方向，同时根据位置关系设置相应的坡长，具体做法是：如果该点位于中心栅格的东西方向，该中心栅格的单元坡长为LONG；如果该点位于中心栅格的南北方向，该中心栅格的单元坡长为LAT；流向为东南、西南、东北、西北，该中心栅格的单元坡长为diagCellSize。

(5)重复以上步骤，直至遍历处理完所有的栅格，完成所有栅格的单元坡长计算。

步骤3：计算栅格初始坡长值；

首先设置截断类型，需先计算河网值，具体包括：

(1)计算汇水面积：

Step1：申请存储汇水的矩阵空间；

Step2：判断栅格是否为无值，如果无值，进行下一个栅格计算，有值进入下一步；

Step3：计算栅格在经度上的长度LONG和和纬度上的长度LAT(计算方法同步骤三)；

Step4：计算栅格的面积S

S＝LONG×LAT

栅格的面积即为该栅格的初始河网值；

Step5：遍历整个数字地图，初始化所有栅格的河网。

(2)河网提取：

Step1：计算栅格c中为8个栅格中流向指向c的栅格的河网值的和，记作sum；

Step2：比较sum与栅格c的河网值，如果sum大于c栅格的河网值，则将sum赋值给c栅格的河网值，否则c栅格的值河网值不变；

Step3：正向遍历整个数字地图，计算整个地图的河网值；

Step4：反向遍历整个数字地图，计算整个地图的河网值；

Step5：如果正反两边的过程中都没有发生将sum的值赋值给c栅格的河网值的操作，说明河网提取完成则结束循环，否则进入step1继续循环计算。

Step6：将计算完成的河网输出到文件中。

初始化栅格坡长需要考虑以下两种截断：(1)坡度变化产生的截断点。

(2)河网截断。

Step1：申请存储截断值的矩阵空间；

step2：判断栅格c的是否为无值，如果是则c点的截断设置成是截断，否则设置成不是截断；

Step3：判断栅格c的坡度是否小于2.861°，如果是而且变化率大于0.7，则栅格c的截断设置成截断；

Step4：判断栅格c的坡度是否大于等于2.861°，如果是而且变化率大于0.5，则c点的截断设置成是截断；

Step5：判断栅格c的河网值是否大于设定的河网截断阈值(河网截断阈值一般是专家经验值)，如果是则c点的截断设置成是截断，否则设置成不截断；

Step6：遍历整个数字地图，完成所有栅格的截断设置。

(3)将初始坡长定义为

其中i为当前栅格，Li代表当前栅格的初始坡长。

Step1：判断栅格c是否截断，如果截断则栅格c的坡长等于原来坡长的一半，如果不截断则栅格c的坡长不发生变化；

Step2：遍历整个数字地图，设置计算每个栅格加入截断后的初始坡长。

步骤4：计算累积坡长。

坡长提取过程是沿着水流方向累加初始坡长的过程，因为计算过程中有累加效果，所以又称为累积坡长，采用正反遍历的方式进行累加，直至栅格坡长值不再发生变化。

Step1：申明临时变量sum初始值设置为0；

Step2：假设栅格a是栅格c周围相邻的一个栅格，而且栅格a的流向指向栅格c，如果栅格a截断，则sum加上栅格a坡长的一半，如果不截断则sum加上栅格a的坡长值；

Step3：用step2的方法计算栅格c周围八个流向栅格c的和；

Step4：比较sum与栅格c的初始坡长值，如果sum大于c栅格的坡长值，则将sum赋值给c栅格的坡长值，否则c栅格的坡长值不变；

Step5：正向遍历整个数字地图，计算整个地图的坡长值；

Step6：反向遍历整个数字地图，计算整个地图的坡长值；

Step7：如果正反两边的过程中都没有发生将sum的值赋值给c栅格的坡长值的操作，说明累计坡长提取完成，结束循环，否则进入step1继续循环计算；

Step8：将计算完成的累计坡长输出到文件中。

至此，算法结束，本发明的方法最终获取到了所有栅格的累积坡长，每一个栅格对应一个累积坡长。

效果验证：

本算法和DEM算法对比：

实验背景：为验证SRTM数据得到坡长结果的正确性，本实验将其计算结果和DEM算法的计算结果作对比、分析。

实验区域：陕北绥德县南沟；实验方法：将县南沟1弧秒SRTM数据(SRTM V1.0)得到的结果进行坐标变换，将结果投影到与30mDEM结果到同一到坐标系下。对比二者的结果。图6基于DEM的县南沟坡长，图7基于SRTM的县南沟坡长。图8是对于SRTM和DEM的坡长结果作差得到的县南沟坡长差值曲线图。

结果分析：由图4和图5可得：基于两种数据得到的LS的整体分布图和坡度值的范围十分相近；通过DEM计算得到的坡长最小值是15m，DEM数字地图的每个栅格的边长都是30m因此当存在截断时，初始的坡长就是15m；而SRTM得到的最小值为12.4m，这是因为通过地球半径计算出了栅格纬度上的跨度较小，通过坡长中的最小值可以推测出，在该测试区域上，栅格的最小纬度上的间距为24.8m，也就是SRTM数据在该纬度上的栅格面积仅为DEM数据同纬度上面积的5/6(按矩形面积计算)，因为地球是一个球体，所以通过SRTM计算的结果应该是更加具有实际意义。通过坡长差值曲线图可以看出，大部分的坡长差值非常小，但也存在少部分差值较大的结果。虽然使用两套数据计算时，将河网截断阈值(一万平方米)设置的相同，但由于不同数据计算得到的栅格面积不同，所以计算出的汇水面积中对应的值也不尽相同，因此会产生截断点不同的现象。假如说a和b是相邻的两个栅格，并且a的流向指向了b，但是在DEM中a和b都由于河网值而截断，但是在SRTM中a和b都没有截断，或者只有b一个截断，那么两个图计算出的坡长结果在b栅格后必然出现巨大的差值。虽然两个数字地图计算的坡长结果中有一些值出现了较大的差值，但是经过综合分析，本方法基于SRTM计算的坡长更为合理。

Claims (5)

1.一种基于SRTM的坡长提取方法，其特征在于，该方法按照以下步骤进行：

步骤1：存储SRTM数据头信息，同时获取SRTM数据中每一个栅格的坡度值和流向；

步骤2：计算栅格单元坡长值；

步骤2.1：遍历所有栅格，判断当前栅格是否是无值点：若是无值点：则设置该栅格的单元坡长值为0，进入下一个栅格的判断；若是非无值点，执行步骤2.2；

步骤2.2：根据步骤1获取的栅格流向，若当前栅格的流向为东西方向，该当前栅格的单元坡长为LONG；若流向为南北方向，该当前栅格的单元坡长为LAT；若流向为东南、西南、东北、西北中的任意一种，则该当前栅格的单元坡长为diagCellSize；

其中，

LAT＝2π×6371000·cellsize/360，

LONG＝2π×6371000·cellsize/360·cosθ，θ为当前栅格中任意一点与地球球心连线与赤道的夹角；

其中，cellsize为SRTM数据中每个栅格对应的头信息中cellsize这一项的取值；

步骤2.3：遍历所有栅格，得到每一个栅格的单元坡长值；

步骤3：计算栅格初始坡长值；

步骤3.1：判断当前栅格是否截断，若为不截断，则当前格栅初始坡长值等于单元坡长值；若为截断，则当前栅格的初始坡长值为单元坡长值的一半；

步骤3.2：遍历所有栅格，得到每一个栅格的初始坡长值；

步骤4：计算栅格累积坡长：

步骤4.1：将当前栅格作为初始累积栅格，初始累积栅格初始值设置为0，任意与初始累积栅格相邻的栅格中流向指向初始累积栅格的定义为相邻指向栅格，若其中任意一个相邻指向栅格截断，则取初始值与该相邻指向栅格的初始坡长值的一半之和作为初始累积栅格对该相邻指向栅格的累积坡长；若该相邻指向栅格不截断，则取初始值与该相邻指向栅格的初始坡长值之和作为初始累积栅格对该相邻指向栅格的累积坡长；

依次计算初始累积栅格对其他相邻指向栅格的累积坡长，比较初始累积栅格与所有相邻指向栅格的累积坡长的总和SUM与初始累积栅格的初始坡长值的大小，若SUM＞初始累积栅格的初始坡长值，则将SUM赋值给初始累积栅格作为初始累积栅格的初始坡长值，否则初始累积栅格的初始坡长值不变；

步骤4.2：将所有剩余栅格中的每个栅格均作为初始累积栅格，按照步骤4.1的方法，遍历所有剩余栅格，计算每个栅格对应的SUM并判断是否赋值，至不再出现赋值情况时遍历结束，最后一次赋值时获得的初始坡长值即为提取的对应栅格的累积坡长。

2.如权利要求1所述基于SRTM的坡长提取方法，其特征在于，步骤1中，存储SRTM数据的头信息，同时获取SRTM数据中每一个栅格的坡度值和流向，包括：

步骤1.1：创建结构体DemData，用于存放SRTM的头信息cellsize；

步骤1.2：按行读取SRTM的数据，SRTM的头信息的格式为：前边是名称，接着是空格，后边是头信息中cellsize这一项；将读取的每行数据存到对应的一个字符串中，然后对该字符串进行空格分割，最后将头信息转换类型并保存到数据结构DemData对应的属性中；

步骤1.3：采用最陡坡降法计算每个单元栅格的坡度值和流向。

3.如权利要求1所述基于SRTM的坡长提取方法，其特征在于，在步骤2之前，对步骤1的结果进行无值点和洼地填充。

4.如权利要求1所述基于SRTM的坡长提取方法，其特征在于，步骤3中：判断当前栅格是否截断，包括以下方法中的任意一种：

a：若当前栅格是无值点，则当前栅格为截断，否则不是截断；

b：若当前栅格的坡度值小于2.861°，且坡度变化率大于0.7，则当前栅格为截断，否则不是截断；

c：若当前栅格的坡度值大于等于2.861°，且坡度变化率大于0.5，则当前栅格为截断，否则不是截断；

d：若当前栅格的河网值大于设定的河网截断阈值，则当前栅格为截断，否则不是截断。

5.如权利要求4所述基于SRTM的坡长提取方法，其特征在于，d中，当前栅格的河网值的计算包括：

第一步：将每个栅格的面积S作为该栅格的初始化河网值，其中，S＝LONG×LAT，遍历所有栅格，得到所有栅格的初始化河网值；

第二步：当前栅格的河网值提取：

Step1：计算当前栅格中流向指向当前栅格的相邻栅格的初始化河网值的和，记作sum；

Step2：比较sum与当前栅格的初始化河网值，如果sum＞当前栅格的初始化河网值，则将sum赋值给当前栅格，否则当前栅格的河网值不变；

Step3：正向遍历所有栅格，计算每个栅格的sum并判断是否赋值；

Step4：反向遍历所有栅格，计算每个栅格的sum并判断是否赋值；

Step5：若正反遍历的过程中不再发生赋值，遍历结束，最后一次赋值作为当前栅格的河网值。

Images