技术领域
[0001] 本
发明属于计算机应用领域，涉及坡度提取，具体涉及一种一种基于SRTM数据的坡度提取方法。
背景技术
[0002] 地形是影响
土壤侵蚀的重要因子，相关研究已经有多年历史。目前的研究已经实现了区域尺度上基于DEM(Digital Elevation Model)的地形因子提取
算法，而高
分辨率遥感影像数据进行大范围的地形因子提取的研究少有报道。通常DEM数据栅格尺寸大小一致，以规整的正方形作为基本单位。SRTM V1.0是美国太空总署(NASA)和国防部国家测绘局(NIMA)以及德国与意大利航天机构共同合作完成联合完成地球表面所进行的遥感测量数据。
覆盖全球80％以上的陆地面积，同时
精度高达约30m。整套数据以1弧度的栅格数据作为基本单位，为geotif格式，每景数据覆盖经纬度各1°，数据量巨大。目前的研究多将SRTM进行格式转化和投影转换后，再利用常用工具(如：ArcGIS等)进行坡度计算。地理坐标下的SRTM数据，拍摄采用经纬度作为单位，投影过程为球面距离，因此SRTM的栅格形状为不规则边长，如何直接使用SRTM数据进行坡度提取，不再进行投影转换，进而为其他地形要素提取提供参考，成为研究的难点。
[0003] 目前基于DEM的坡度提取算法较多，核心思想是对栅格进行遍历，比较中心栅格的高程值和周围八个格子的高程值，通过最陡坡降、三阶反距离平方权、三阶不带权差分等方法确定坡度值。由于栅格都是大小相等的规整正方形，边长相等，计算坡度相对容易。常用的坡度计算方法如：最陡坡降法，
[0004] 全球土壤侵蚀问题是关乎全人类共同命运的严峻生态问题，为了从宏观上把控地形的侵蚀变化，高分辨率全球土壤侵蚀制图亟待解决。本发明设计并实现了基于SRTM的坡度提取算法，其为计算全球SRTM地形因子提供了解决方案和技术支持。
[0005] 基于ArcGIS的坡度提取算法，需要先经过坐标变换，然后再使用坡度分析工具，才能实现坡度的提取，其对操作人员专业性的要求较高。
发明内容
[0006] 针对
现有技术存在的不足，本发明的目的在于，提供一种基于SRTM数据的坡度提取方法，解决现有技术中的坡度提取算法对操作人员专业性的要求较高，且无法在大空间范围进行坡度提取的技术问题。
[0007] 为了解决上述技术问题，本发明采用如下技术方案予以实现：
[0008] 一种基于SRTM数据的坡度提取方法，该方法按照以下步骤进行：
[0009] 步骤一，读取数据：
[0010] 步骤S11，读取文件头部；
[0011] 步骤S12，读取SRTM高程数据；
[0012] 步骤二，无值点和洼地填充：
[0013] 步骤S21，对SRTM高程数据进行无值点填充：
[0014] 步骤S22，对经过无值点填充后的SRTM高程数据进行洼地填充：
[0015] 步骤三，遍历二维数组，计算坡度值：
[0016] 步骤S31，对经过无值点和洼地填充的SRTM高程数据遍历二维数组，在每次遍历时，先判断SRTM高程数据中当前栅格是否是无值点，如果是无值点，则直接设置坡度值为0并且跳过，进入下一个栅格的判断；若果是非无值点，则执行步骤S32；
[0017] 所述的栅格为SRTM高程数据的基本单元；
[0018] 步骤S32，如果是非无值点：则进行以下计算。
[0019] 设置地球的半径值R＝6371000m，通过获取到的头信息cellsize的值；
[0020] 计算当前栅格在经线方向上的像元宽度XCellsize；
[0021] XCellsize＝2π×6371000·cellsize/360；
[0022] 计算当前栅格在
纬线方向上的像元宽度YCellsize，
[0023] YCellsize＝2πR·cosα·cellsize/360；
[0024] 计算当前栅格的对
角线方向的像元距离diagCellSize；
[0025]
[0026] 其中，cellsize为SRTM数据中文件头部中用于标识栅格大小的数据；α为当前栅格的纬度；
[0027] 步骤S33，根据最陡坡降法确定坡度值：
[0028] Step1：依次将当前栅格和周围的8个栅格的高程值进行对比，如果周围某个栅格的高程值不是无值点，并且高程值大于当前栅格，则计算该周围某个栅格与当前栅格所成角度angle＝deg·arctan((DC-DO)/CellStep)；
[0029] 其中：
[0030] DC代表当前栅格的高程值；
[0031] DO代表周围某个栅格的高程值；
[0032] CellStep为周围某个栅格与当前栅格的距离，若周围某个栅格与当前栅格的
位置关系为东西方向，则CellStep为YCellsize，若周围某个栅格与当前栅格的位置关系为南北方向，则CellStep为XCellsize，若周围某个栅格与当前栅格的位置关系为南北方向，则CellStep为diagCellSize；
[0033] Step2：将周围8个栅格计算出的angle最大值确定为当前栅格的坡度值；
[0034] 步骤S34，依次将当前栅格和周围的8个栅格的高程值进行对比，如果周围某个栅格的高程值不是无值点，并且8个栅格的高程值均小于等于当前栅格，则当前栅格对应的地形为平地或洼地，设置当前栅格的坡度最小值为0.1；
[0035] 步骤S35，重复以上步骤S31至步骤S34，直至遍历处理完所有的栅格，完成基于SRTM数据的坡度提取。
[0036] 本发明还具有如下技术特征：
[0037] 步骤S11中读取文件头部的具体过程为：
[0038] Step1：创建一个结构体DemData，用于存放SRTM的文件头信息；
[0039] Step2：打开SRTM文本文件，如果打开失败写入日志并退出程序；
[0040] Step3：按行读取SRTM文本文件中的内容，SRTM文本文件的文件头部分的格式为：前边是名称，接着是空格，后边是值；
[0041] Step4：将读取的一行数据存到一个字符串中，然后用对该字符串进行以空格分割，然后将值转换成该值的类型并保存到数据结构DemData中对应的属性中；
[0042] Step5：重复step4的步骤读完SRTM文本文件的文件头部。
[0043] 步骤S12中读取SRTM高程数据的具体过程为：所述的SRTM高程数据的每一行的每个数据之间是以空格分开的。读取的过程中按行读取，读取每一行的字符串形式的数据，利用空格进行分割，将分割出的每一个数据部分的字符串，转换成float类型的数据存到高程数据矩阵中。
[0044] 本发明与现有技术相比，具有如下技术效果：
[0045] (Ⅰ)本发明的方法直接对原始的SRTM数据进行分析，计算出坡度值，这样提取坡度的效率更高。
[0046] (Ⅱ)从算法的思想上看，本算法根据栅格所在经纬度，先计算获得格网四条边的长度，进而再以最陡坡降的方向，提取坡度值，坡度提取符合理论常识，具有可操作性。
[0047] (Ⅲ)从算法实现的最终结果来看：将计算得到的坡度结果进行投影转换，与平面坐标下的DEM数据得到的坡度结果99％以上误差在0.5°以内，本算法计算结果可接受。
附图说明
[0048] 图1是基于SRTM的坡度提取
流程图。
[0049] 图2是SRTM栅格的俯视图示意图。
[0050] 图3是图2为栅格中心像元与8邻域像元示意图。
[0051] 图4是基于DEM的县南沟坡度结果。
[0052] 图5是基于SRTM的县南沟坡度结果。
[0053] 图6是SRTM与DEM坡度差值
频率统计图。
[0054] 以下结合
实施例对本发明的具体内容作进一步详细解释说明。
具体实施方式
[0055] 如何在SRTM的地理
坐标系下实现成为本
专利的研究创新。本发明设计实现基于SRTM的坡度提取算法，将极大的方便用户使用，为大尺度地形因子提取提供技术
支撑。基于SRTM的坡度提取算法的设计和实现，有效的避免了SRTM数据格式及投影的变换过程，能够提升计算效率，对于在大空间范围应用土壤侵蚀模型起到积极的促进作用。
[0056] SRTM数据中，栅格并非传统的规整正方形，每个栅格的俯视图类似梯形，而梯形的每条边都只能通过地球半径和经纬度进行计算。如图2所示：
[0057] 图2(a)为地球过地轴的旋转面，O为球心，R为赤道半径，α为地球表面某一点与球心的连线与赤道面的夹角，其值表示这一点的纬度，r为这一点纬线圈的半径。X表示夹角为α时，对应的球面距离。O’表示纬线圈所在平面与地轴的交点。图2(b)为纬线圈所在平面，Y表示纬线圈上的一段球面距离，β表示这段距离对应的经度差。
[0058] 用R表示地球的平均半径，则有公式1和公式2：
[0059] X＝2πR·α/360 (1)
[0060] Y＝2πR·β/360 (2)
[0061] 由r＝R·cosα有公式3
[0062] Y＝2πR·cosα·β/360 (3)
[0063] 对于地理坐标下的SRTM数据，若X表示经线方向上的像元宽度(距离)，即X＝XCellsize，Y表示纬线方向上的像元宽度(距离)，即Y＝YCellsize，那么α与β分别表示像元南北宽度(度)A和东西宽度(度)B。在STRM数据中像元栅格所跨越的经度和纬度是相同的，即A＝B。在取R＝6371000m，取β＝1/3600则有公式4和公式5：
[0064] XCellsize＝30.8874791 (4)
[0065] YCellsize＝30.8874791·cosα (5)
[0066] 图2为栅格中心像元与8邻域像元示意图，a,b,c,d,e,f,g,h,i分别表示相应像元的高程值。坡度公式为：
[0067]
[0068] 其中：
[0069]
[0070]
[0071] 以下给出本发明的具体实施例，需要说明的是本发明并不局限于以下具体实施例，凡在本
申请技术方案
基础上做的等同变换均落入本发明的保护范围。
[0072] 实施例1：
[0073] 遵从上述技术方案，如图1至图3所示，本实施例给出一种基于SRTM数据的坡度提取方法，该方法按照以下步骤进行：
[0074] 步骤一，读取数据：
[0075] 步骤S11，读取文件头部；
[0076] 步骤S11中读取文件头部的具体过程为：
[0077] Step1：创建一个结构体DemData，用于存放SRTM的文件头信息；
[0078] Step2：打开SRTM文本文件，如果打开失败写入日志并退出程序；
[0079] Step3：按行读取SRTM文本文件中的内容，SRTM文本文件的文件头部分的格式为：前边是名称，接着是空格，后边是值；
[0080] Step4：将读取的一行数据存到一个字符串中，然后用对该字符串进行以空格分割，然后将值转换成该值的类型并保存到数据结构DemData中对应的属性中；
[0081] Step5：重复step4的步骤读完SRTM文本文件的文件头部。
[0082] 步骤S12，读取SRTM高程数据；
[0083] 步骤S12中读取SRTM高程数据的具体过程为：所述的SRTM高程数据的每一行的每个数据之间是以空格分开的。读取的过程中按行读取，读取每一行的字符串形式的数据，利用空格进行分割，将分割出的每一个数据部分的字符串，转换成float类型的数据存到高程数据矩阵中；为了计算的方便将数据扩充两圈，所以在读完数据后就是把数据放进一个扩大了两圈的矩阵中。
[0084] 步骤二，无值点和洼地填充：
[0085] 步骤S21，对SRTM高程数据进行无值点填充：
[0086] 步骤S22，对经过无值点填充后的SRTM高程数据进行洼地填充：
[0087] 无值点一般是在卫星飞行测量产生的错误数据。根据地形
水文学的实际情况，在地形中不存在周围都有数据而中心没有数据的栅格，因此需要将这些栅格填充为周围八个栅格的最小高程值。
[0088] 所述的洼地填充的方式是填充周围最低栅格的高程。
[0089] 步骤三，遍历二维数组，计算坡度值：
[0090] 步骤S31，对经过无值点和洼地填充的SRTM高程数据遍历二维数组，在每次遍历时，先判断SRTM高程数据中当前栅格是否是无值点，如果是无值点，则直接设置坡度值为0并且跳过，进入下一个栅格的判断；若果是非无值点，则执行步骤S32；
[0091] 所述的栅格为SRTM高程数据的基本单元；
[0092] 步骤S32，如果是非无值点：则进行以下计算。
[0093] 设置地球的半径值R＝6371000m，通过获取到的头信息cellsize的值；
[0094] 计算当前栅格在经线方向上的像元宽度XCellsize；
[0095] XCellsize＝2π×6371000·cellsize/360；
[0096] 计算当前栅格在纬线方向上的像元宽度YCellsize，
[0097] YCellsize＝2πR·cosα·cellsize/360；
[0098] 计算当前栅格的对角线方向的像元距离diagCellSize；
[0099]
[0100] 其中，cellsize为SRTM数据中文件头部中用于标识栅格大小的数据；α为当前栅格的纬度；
[0101] 步骤S33，根据最陡坡降法确定坡度值：
[0102] Step1：依次将当前栅格和周围的8个栅格的高程值进行对比，如果周围某个栅格的高程值不是无值点，并且高程值大于当前栅格，则计算该周围某个栅格与当前栅格所成角度angle＝deg·arctan((DC-DO)/CellStep)；
[0103] 其中：
[0104] DC代表当前栅格的高程值；
[0105] DO代表周围某个栅格的高程值；
[0106] CellStep为周围某个栅格与当前栅格的距离，若周围某个栅格与当前栅格的位置关系为东西方向，则CellStep为YCellsize，若周围某个栅格与当前栅格的位置关系为南北方向，则CellStep为XCellsize，若周围某个栅格与当前栅格的位置关系为南北方向，则CellStep为diagCellSize；
[0107] Step2：将周围8个栅格计算出的angle最大值确定为当前栅格的坡度值；
[0108] 步骤S34，依次将当前栅格和周围的8个栅格的高程值进行对比，如果周围某个栅格的高程值不是无值点，并且8个栅格的高程值均小于等于当前栅格，则当前栅格对应的地形为平地或洼地，设置当前栅格的坡度最小值为0.1；
[0109] 步骤S35，重复以上步骤S31至步骤S34，直至遍历处理完所有的栅格，完成基于SRTM数据的坡度提取。
[0110] 实验验证：
[0111] 本发明的基于SRTM数据的坡度提取方法和DEM算法对比：
[0112] 实验背景：SRTM数据得到的计算结果的正确性验证，最直接有效的测试方式就是将其计算结果和目前已经成熟的DEM算法作对比、分析。
[0113] 实验区域：陕北绥德县南沟；实验方法：将县南沟1弧秒SRTM数据(SRTM V1.0)得到的结果进行坐标变换，将结果投影到与30mDEM结果到同一到坐标系下。对比二者的结果。图4基于DEM的县南沟坡度，图5基于SRTM的县南沟坡度。图6是对于SRTM和DEM的坡度结果作差扩大100倍得到的县南沟坡度差值频率统计图。
[0114] 结果分析：由图4和图5可得：基于两种数据得到的LS的整体分布图和坡度值的范围十分相近；从坡度差值频率统计图可见：99％的差值集中在±0.5之间。SRTM数据栅格在经度上的跨越长度是固定的，而在纬度上的跨越长度会随着纬度的增大而减小，所以不在赤道上的SRTM数据的栅格大小是略小于DEM30米栅格的。在结果中，DEM数据与STRM数据在相同位置的高程差是相等的；但是如果流向不是指向南北，那么SRTM数据中栅格间距是小于DEM的，所以坡度就自然大于DEM数据的坡度；而当流向是指向南北的话，那么SRTM的栅格间距略大于DEM栅格的间距，所以也存在STRM中的坡度大于DEM中的坡度值的情况。所以会出现SRTM计算的坡度结果与DEM结果有略微不同。但是二者整体结果无限接近足以表明本发明的正确性、可用性。
