一种基于DEM的河谷形态参数自动提取方法与流程

本发明涉及数字地貌分析领域，尤其涉及一种基于DEM的河谷形态参数自动提取方法。

背景技术：
河谷纵剖面和河谷横剖面形态参数计算是河谷地貌研究的重要内容：河谷纵剖面形态参数可以反映河流物理侵蚀过程和各河段流水下切的动力条件，揭示河流地貌对岩性、构造、气候等条件的相应机制；河谷横剖面是由谷底河道及两岸斜坡组成的几何体，其形态沿河流方向不断变化，其形态参数有助于判断河谷类型和发育程度，辅助开展河谷地貌的深入分析。目前河谷地貌分析中河谷纵剖面和横剖面形态参数的自动提取并未引起足够重视。大多现有水文分析软件无法方便地完成河谷纵剖面和横剖面形态参数的自动提取，常需要通过地图量算或繁琐的软件操作完成，具有操作复杂、精度低、难以推广应用等缺点。尤其是河谷横剖面形态参数及其应用更鲜被提及。常用的地貌特征参数，如坡度、坡向、地形粗糙度、起伏度、切割深度等基于DEM通过GIS邻域分析获得，无法从微观上对河谷剖面的几何形态进行定量计算和表征。这在很大程度上影响了河谷地貌演化及相关研究的范围和深度。

技术实现要素：
本发明的目的在于提供一种基于DEM的河谷形态参数自动提取方法，从而解决现有技术中存在的前述问题。为了实现上述目的，本发明所述基于DEM的河谷形态参数自动提取方法，该方法包括：S1，获取需分析河谷流域的DEM数据；S2，通过水文分析获取所述河谷流域的线状河流矢量数据，对属于所述河谷流域的每一条河流进行编号，并从每条河流的河源开始按一定采样间隔a取纵剖面采样点直至河口，并从第一个纵剖面采样点开始，依次给每个点添加索引编号i，i＝1,2,3,…,n；计算各纵剖面采样点的河长值和高程值，所述河长值leni＝a×i，i表示采样点的编号，所述i≥1；所述a为大于1.414×DEM空间分辨率且小于2.828×DEM空间分辨率的正整数；S3，如果计算河谷纵剖面的形态参数，则进入S4；如果计算河谷横剖面的形态参数，则进入S5；S4，按照公式(1)计算河谷纵剖面任意一个纵剖面采样点的标准化河长-坡降指数SLi′；SLi′＝SLi/K(1)；SLi表示任意一个纵剖面采样点的河长-坡降指数，K表示所述纵剖面采样点所在河流的均衡坡降指数；S5，选取一个河流沿岸的任意一个点作为定位点并获得与所述定位点距离最近的纵剖面采样点T及纵剖面采样点T上、下游各三个纵剖面采样点；获得并利用7个采样点的平面坐标，基于最小二乘法计算纵剖面采样点T的河流流向，将所述纵剖面采样点T的河流流向作为所述定位点处的河流流向；获取垂直于所述定位点处河流流向的所述定位点处的河谷横剖面在水平面上的投影线方向；基于DEM通过地形分析并获取与所述河谷相关的两岸山脊的线状矢量数据，利用所述定位点、所述定位点处的河谷横剖面在水平面上的投影线方向、所述定位点处河流流向和所述山脊线状矢量数据进行几何构建，得到河谷横剖面的水平面投影线；按照一定采样间隔b从河谷横剖面的水平面投影线左端点开始提取河谷横剖面采样点直至右端点；所述b为大于1.414×DEM空间分辨率且小于2.828×DEM空间分辨率的正整数；将最接近河流线的河谷横剖面采样点的编号id设定为0，从0点开始分别向其左、向其右给每一个河谷横剖面采样点编号id，设在所述0点左岸的河谷横剖面采样点的编号小于0，记为id＝l，l＝-1,-2,-3……-n，在所述0点右侧的河谷横剖面采样点的编号大于0，记为id＝r，r＝1,2,3……n；计算各个河谷横剖面采样点的高程值和河谷横剖面采样点到0点的水平距离，计算河谷横剖面坡高h、目标深度的河谷宽深比R、计算河谷横剖面的剖面面积A、河谷凹度C、不对称系数M和河宽指数W。优选地，公式(1)中河谷纵剖面上任意一个采样点的河长-坡降指数SLi按照下述实现：设定坡降计算间隔D，其中，D＝2×j×a，j＝1,2,3…,n，从每条河流的第一个纵剖面采样点开始计算长度为D的河段的坡降S，并把得到的S赋值给河段中心点采样点Fi，即设定坡降计算间隔D时，河谷纵剖面上任意一个采样点Fi的坡降Si按照公式(2)计算得到：Si＝(Hi-j-Hi+j)/(2×j×a)(2)；Hi-j、Hi+j分别表示编号为i-j、i+j的纵剖面采样点的高程。更优选地，坡降计算间隔D为最优计算间隔，所述最优坡降计算间隔D通过对比分析确定，具体为：取n个坡降计算间隔d，d＝2×j×a，j＝1,2,3…,n，分别计算d取不同值时，所有河谷纵剖面采样点的坡降Si，并将其作为集合，统计所述集合中Si＜0的纵剖面采样点的数量之和在所有纵剖面采样点中所占的比重p，统计所述集合中Si的取值区间长度q＝Smax-Smin，选择p小于预设值，q大于预设值时的唯一坡降计算间隔，将此唯一坡降计算间隔作为最优坡降计算间隔D，并将采用坡降计算间隔D时计算得到纵剖面采样点坡降Si作为最优坡降；所述预设值的取值根据河谷地形复杂性和DEM数据质量确定；其中，河谷纵剖面上任意一个采样点的河长-坡降指数SLi按照公式(3)计算获得：SLi＝Si×leni(3)；Si表示编号为i的纵剖面采样点的坡降值，leni表示编号为i的纵剖面采样点的河长值。优选地，河流的均衡坡降指数K按照下述计算：设定河谷中任意一条河流的第一个纵剖面采样点编号为1，最后一个纵剖面采样点编号为N，按照公式(4)计算该河流的均衡坡降指数K，K＝(HN-H1)/(lg(lenN)-lg(len1))(4)；HN表示编号为N的的纵剖面采样点的高程值，H1表示编号为1的纵剖面采样点的高程值；lenN表示编号为N的的纵剖面采样点的河长值，len1表示编号为1的纵剖面采样点的河长值。优选地，按照公式(5)计算河谷横剖面坡高h：h＝min(h1，h2)(5)；h1＝Hmax_ld-H0表示左岸斜坡高差，h2＝Hmax_rd-H0表示右岸斜坡高差，H0表示编号为0的河谷横剖面采样点的高程值，Hmax_ld表示左岸河谷横剖面采样点的高程最大值，Hmax_rd表示右岸河谷横剖面采样点的高程最大值。更优选地，计算目标深度h′的河谷宽深比R，按照下述实现：根据河谷横剖面坡高h和目标深度h′，其中，h′<h，遍历河谷横剖面采样点，得到左岸河谷横剖面采样点中的高程最接近h′+H0的河谷横剖面采样点Xl和右岸河谷横剖面采样点中的高程最接近h′+H0的河谷横剖面采样点Xr；获取点Xl到编号为0的河谷横剖面采样点的水平距离len_toRiv_l、点Xr到编号为0的河谷横剖面采样点的水平len_toRiv_r，按照公式(6)计算目标深度为h′处的河谷宽深比R；R＝(len_toRiv_l+len_toRiv_r)/h′(6)；其中，H0表示编号为0的河谷横剖面采样点的高程值。优选地，按照公式(7)计算河谷横剖面的剖面面积A，按照公式(8)计算河谷凹度C，h表示河谷横剖面的坡高，Hid表示编号为id的横剖面采样点的高程值，b表示河谷横剖面水平面投影线的采样间隔，H0表示编号为0的河谷横剖面采样点的高程值，id_min表示左岸河谷横剖面采样点中的高程最接近h+H0的河谷横剖面采样点的编号，id_max表示右岸河谷横剖面采样点中的高程最接近h+H0的河谷横剖面采样点的编号。优选地，遍历得到左岸、右岸河谷横剖面采样点中的高程最接近h+H0的河谷横剖面采样点Yl、Yr，分别获取横剖面采样点Yl、Yr到编号为0的河谷横剖面采样点的水平距离，按照公式(9)计算河谷的不对称系数M，M＝Ll/Lr(9)；Ll表示横剖面采样点Yl到编号为0的河谷横剖面采样点的水平距离；Lr表示横剖面采样点Yr到编号为0的河谷横剖面采样点的水平距离。优选地，河宽指数W按照下述步骤实现：设定角度常数θ，遍历左岸河谷横剖面采样点，得到左岸河谷横剖面采样点Zl，使其与编号为0的河谷横剖面采样点所连直线和水平线所成锐角最接近θ，遍历得到右岸河谷横剖面采样点Zr，使其与编号为0的河谷横剖面采样点所连直线和水平线所成锐角最接近θ；分别获取点Zl、Zr到0点的水平距离，按照公式(10)计算河宽指数W，W＝Wl+Wr(10)；Wl表示横剖面采样点Zl到编号为0的河谷横剖面采样点的水平距离；Wr表示横剖面采样点Zr到编号为0的河谷横剖面采样点的水平距离。本发明的有益效果是：(1)面向对象：以沿河流分布的任意一点为对象，自动提取其标准化河长-坡降指数和所在河谷横剖面的形态特征参数，几何意义表征明确。沿河每一个采样点都可以参与计算，能够从微观上表征河谷地貌及分段差异，有利于分析河谷的整体地形。(2)表征全面：河谷纵剖面的标准化河长-坡降指数和河谷横剖面坡高、宽深比、剖面面积、河谷凹度、不对称系数、河宽指数等形态特征参数，可以共同反映河流的地势变化及河谷整体形态和切割程度等，信息表达全面、完整。(3)干扰因素小，精度更高：本方法所提河谷形态参数基于河谷剖面几何形态进行计算提取，输入数据包括DEM、定位点位置、采样间隔、计算间隔，输入条件少，无假设条件，误差主要来源于DEM质量和空间采样，人为操作干扰程度低，参数计算精度得到提高。(4)运算效率提高：河谷参数计算以沿河分布的定位点为对象，待计算点数目大幅降低，有利于提高计算效率和节约存储空间。同时可以一次性获取所有特征参数的属性，操作简单，具有可维护性。附图说明图1是河流纵剖面参数自动提取流程图；图2是河流横剖面参数自动提取流程图；图3是河谷横剖面形态参数几何示意图；图4是河谷横剖面平面线段示意图，表示定位点，表示河流采样点，表示河谷横剖面端点，表示河流，表示山脊，表示河流流向，表示河谷横剖面的水平面投影线；图5是微地貌点SL’提取示意图，表示河流，表示地貌点，表示普通纵剖面采样点，表示计算纵剖面采样点；图6是研究区河流河长-均衡坡降L-K关系图；图7是麻曲纵剖面高程与坡降分布图；图8是研究区所有河流采样点S值统计图；图9是堵江点、岩性点和断层点的SL’统计盒形图；图10是滑坡堵江自动提取阈值设定示意图；图11是滑坡体堵江预测建模样本数据处理获取图。具体实施方式为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施方式仅仅用以解释本发明，并不用于限定本发明。本发明所述基于DEM的河谷形态参数自动提取方法，参照图1和图2该方法按照以下步骤实现：S1，获取需分析河谷流域的DEM数据；S2，通过基于ArcGIS软件的水文分析模块获取所述河谷流域的线状河流矢量数据，对属于所述河谷流域的每一条河流进行编号，并从每条河流的河源开始按一定采样间隔a取纵剖面采样点直至河口，并从第一个纵剖面采样点开始，依次给每个点添加索引编号i，i＝1,2,3,…,n；计算各纵剖面采样点的河长值和高程值，所述河长值leni＝a×i，i表示采样点的编号，所述i≥1；所述a为大于1.414×DEM空间分辨率，小于2.828×DEM空间分辨率的正整数；S3，如果计算河谷纵剖面的形态参数，则进入S4；如果计算河谷横剖面的形态参数，则进入S5；S4，按照公式(1)计算河谷纵剖面任意一个纵剖面采样点的标准化河长-坡降指数SLi′；SLi′＝SLi/K(1)；SLi表示任意一个纵剖面采样点的河长-坡降指数，K表示所述纵剖面采样点所在河流的均衡坡降指数；S5，选取一个河流沿岸的定位点并获得与所述定位点距离最近的纵剖面采样点T及纵剖面采样点T上、下游各三个纵剖面采样点，利用上述7个采样点的平面坐标，基于最小二乘法计算纵剖面采样点T处的河流流向作为所述定位点处的河流流向，所述定位点处的河谷横剖面的水平面投影线的方向垂直于该点的河流流向，基于DEM通过地形分析获取与所述河谷相关的两岸山脊的线状矢量数据，利用所述定位点、所述河谷横剖面水平面投影线方向和所述线状山脊矢量数据进行几何构建，得到河谷横剖面的水平面投影线，按照一定采样间隔b从河谷横剖面水平面投影线左端点开始提取横剖面采样点直至右端点；所述b为大于1.414×DEM空间分辨率，小于2.828×DEM空间分辨率的正整数；将最接近河流线的横剖面采样点的编号id设定为0，从0点开始分别向其左、向其右给每一个横剖面采样点编号id，设在所述0点左岸的横剖面采样点的编号小于0，记为id＝l，l＝-1,-2,-3……-n，在所述0点右侧的横剖面采样点的编号大于0，记为id＝r，r＝1，2,3……n；计算各个横剖面采样点的高程值和横剖面采样点到0点的距离，计算河谷横剖面坡高、目标深度的河谷宽深比、河谷横剖面的剖面面积、河谷凹度、不对称系数和河宽指数。更详细的解释说明：(一)在步骤S2中，使每条河流都是独立唯一的线要素，并保证线要素的方向与河流流向一致；(二)步骤S4中各个参数的计算方法2.1、公式(1)中河谷纵剖面上任意一个采样点的河长-坡降指数SLi按照下述实现：设定坡降计算间隔D，其中，D＝2×j×a，j＝1,2,3…,n，从每条河流的第一个纵剖面采样点开始计算长度为D的河段坡降S，并把得到的S赋值给河段中心点采样点Fi，即设定坡降计算间隔D时，河谷纵剖面上任意一个采样点Fi的坡降Si按照公式(2)计算得到：Si＝(Hi-j-Hi+j)/(2×j×a)(2)；Hi-j、Hi+j分别表示编号为i-j、i+j的纵剖面采样点的高程。其中每条河流首j个采样点和尾j个纵剖面采样点不参与计算，记Si＝-9999；坡降计算间隔D是纵剖面采样间隔的偶数倍，表示河流的长度。坡降计算间隔D为最优计算间隔，所述最优坡降计算间隔D通过对比分析确定：依次取n个坡降计算间隔d，d＝2×j×a，j＝1,2,3…,n，分别计算d取不同值时，所有河谷纵剖面采样点的坡降Si，并将其作为集合，统计所述集合中Si＜0的纵剖面采样点数的数量之和在所有纵剖面采样点中所占的比重p，统计所述集合中Si的取值区间q＝Smax-Smin，选择p小于预设值，q大于预设值时的唯一计算间隔，将此唯一计算间隔作为最优坡降计算间隔D，并将采用坡降计算间隔D时计算得到纵剖面采样点坡降Si作为最优坡降；所述预设值的取值根据河谷地形复杂性和DEM数据质量确定；表1表示河流纵剖面采样点Fi的坡降计算公式：表1河流纵剖面采样点Fi的坡降计算公式河谷纵剖面上任意一个采样点的河长-坡降指数SLi按照公式(3)计算获得：SLi＝Si×leni(3)；Si表示编号为i的纵剖面采样点的坡降值，leni表示编号为i的纵剖面采样点的河长值，即所述点距离所在河流源头的距离。SL强调河谷的微观地形特征，对河谷坡降变化敏感，同时减轻了采样点坡降随河长增加逐渐减小这一现象对地形监测的干扰，是河谷纵剖面变化识别和形成机制分析的有效工具。SL受多种因素影响，包括基岩性质、构造活动、气候变化、支流汇入和斜坡运动等。沿河流方向若SL在很短的距离内由小到大产生很大增量，说明此处很可能有裂点存在。2.2、河流的均衡坡降指数K按照下述计算：设定每条河流的第一个纵剖面采样点编号为1，最后一个纵剖面采样点编号为N，按照公式(4)计算采样点所在河流的均衡坡降指数K，K＝(HN-H1)/(lg(lenN)-lg(len1))(4)；HN表示编号为N的的纵剖面采样点的高程值，H1表示编号为1的的纵剖面采样点的高程值；lenN表示编号为N的的纵剖面采样点的河长值，len1表示编号为1的的纵剖面采样点的河长值。由于不同河流的形成背景和发育程度各不相同，各河的均衡坡降K值存在较大差异。一般发育程度高、河流长度大的大型河流的K比小型河流大。为了能够在不同河流间开展地貌特征对比分析，需要采用标准化的河长-坡降指标SL’表征河流纵剖面的形态变化。(三)步骤S5中各河谷横剖面形态参数的计算方法，参数的几何意义如图3所示；3.1定位点流向的计算方法：选取一个河流沿岸的定位点并获得与所述定位点距离最近的纵剖面采样点T及纵剖面采样点T上、下游各三个纵剖面采样点，利用上述7个采样点的平面坐标，基于最小二乘法计算纵剖面采样点T处的河流流向作为所述定位点处的河流流向，所述定位点处的河谷横剖面的水平面投影线的方向垂直于该河流流向，具体为:最小二乘法是一种传统的数学估计方法，在求解不同变量依赖关系的经验公式时，要求误差的平方和达到最小。该方法被普遍应用于二维点的直线拟合中，原理如下：假设一组成对二维数据点为(x1,y1；x2,y2；...xn,yn)，它们在x-y直角坐标系中集中分布在一条直线附近，该直线可以表达为y＝kx+b，其中k、b是任意常数。为了对k、b进行参数估计，要求所有二维点的估计值yi’与实际值yi离差的平方和最小，即最小。令当离差的平方和取最小值时，式3.1.1对k、b求偏导的值应为0，可得：得出：基于所述7个纵剖面采样点的坐标利用公式3.1.4计算出的k是纵剖面采样点T处河流流向的斜率。3.2、河谷两岸山脊的线状矢量数据的获取方法：以所述河谷流域的DEM数据和线状河流矢量数据为数据源，借助ArcGIS软件的空间分析模块进行地形分析，识别出山脊的大致分布，在此基础上进行编辑，获取线状山脊矢量数据，为每条山脊线添加两个属性，分别是关联河流编号和相对位置，表示该山脊位于哪条河流的左岸还是右岸。3.3、将所述定位点、所述河谷横剖面的水平面投影线方向和所述线状山脊矢量数据进行几何构建，得到河谷横剖面的水平面投影线，具体按照下述实现：过定位点的河谷横剖面的水平面投影线是过所述定位点连接两岸山脊并垂直于所述定位点处河流流向的线段。参照图4中，若定位点处河流流向的斜率为k，则过所述定位点的河谷横剖面的水平面投影线的斜率为K＝-1/k。过所述定位点以K＝-1/k做直线，交所述左岸山脊线于一点，即河谷横剖面水平面投影线的左端点；交所述右岸山脊线于一点，即河谷横剖面水平面投影线的右端点。连接上述两端点得到过定位点的河谷横剖面的水平面投影线。3.4、按照公式(5)计算河流横剖面坡高h：h＝min(h1，h2)(5)；h1＝Hmax_ld-H0表示左岸斜坡高差，h2＝Hmax_rd-H0表示右岸斜坡高差，H0表示编号为0的河谷横剖面采样点的高程值，Hmax_ld表示左岸河谷横剖面采样点的高程最大值，Hmax_rd表示右岸河谷横剖面采样点的高程最大值。3.5、计算目标深度的河谷宽深比R，按照下述实现：根据河谷横剖面坡高h和目标深度h′，其中h′<h，遍历河谷横剖面采样点，得到左岸河谷横剖面采样点中的高程最接近h′+H0的河谷横剖面采样点Xl和右岸河谷横剖面采样点中的高程最接近h′+H0的河谷横剖面采样点Xr；，其中H0表示编号为0的河谷横剖面采样点的高程值。获取左岸采样点Xl到编号为0的河谷横剖面采样点的水平距离len_toRiv_l、右岸采样点Xr到编号为0的河谷横剖面采样点的水平len_toRiv_r，按照公式(6)计算目标深度为h′处的河流宽深比R；R1＝(len_toRiv_l+len_toRiv_r)/h′(6)；3.6、按照公式(7)计算河谷横剖面的剖面面积A，按照公式(8)计算河谷凹度C，h表示河流坡高，Hid表示编号为id的横剖面采样点的高程值，b表示河谷横剖面采样间隔，H0表示编号为0的横剖面采样点的高程值，id_min表示左岸河谷横剖面采样点中的高程最接近h+H0的河谷横剖面采样点的编号，id_max表示右岸河谷横剖面采样点中的高程最接近h+H0的河谷横剖面采样点的编号。河谷凹度C可表征河谷形态类型，一般认为，C＞1时，所述河谷为“U”型谷；C≤1时，所述河谷为“V”型谷。3.7、计算河谷不对称系数M，按照下述方式实现：遍历得到左岸、右岸河谷横剖面采样点中的高程最接近h+H0的河谷横剖面采样点Yl、Yr；计算横剖面采样点Yl、Yr到编号为0的河谷横剖面采样点的水平距离，按照公式(9)计算山脊的不对称系数M，M＝Ll/Lr(9)；Ll表示横剖面采样点Yl到编号为0的河谷横剖面采样点的水平距离；Lr表示横剖面采样点Yr到编号为0的河谷横剖面采样点的水平距离。不对称系数M主要反映河谷横剖面两岸地形的对称性，当M＝1时，河谷横剖面形态对称；M<1时，河谷横剖面为左陡右缓的右倾坡；当M＞1时，河谷横剖面为右陡左缓的左倾坡。3.8、河宽指数W按照下述方式实现：设定角度常数θ，遍历左岸河谷横剖面采样点，得到左岸河谷横剖面采样点Zl，使其与编号为0的河谷横剖面采样点所连直线和水平线所成锐角最接近θ；同样遍历得到右岸河谷横剖面采样点Zr，使其与编号为0的河谷横剖面采样点所连直线和水平线所成锐角最接近θ；分别获取横剖面采样点Zl、Zr到编号为0的河谷横剖面采样点的水平距离，按照公式(10)计算河宽指数W，W＝Wl+Wr(10)；Wl表示横剖面采样点Zl到编号为0的河谷横剖面采样点的水平距离；Wr表示横剖面采样点Zr到编号为0的河谷横剖面采样点的水平距离，θ的取值根据河谷地形复杂性和DEM数据质量确定。若W＝0说明河道宽度很窄，取W为河谷横剖面的采样间隔b。本发明基于DEM的河谷形态参数自动提取，借助现有GIS软件平台和GIS二次程序编写的方法实现，是GIS技术在河谷地貌学中的综合应用。表2中列举了其中应用的主要GIS功能。ArcGISDesktop等在现有GIS软件提供的相关功能模块可以基于DEM完成水文分析和地形分析，生成线状河流、河流流向和山脊矢量数据。本申请基于软件操作，按照地貌参数提取的原理和方法，对软件输出数据进行编辑修改，使之能够输入后续开发的GIS程序模块中参与运算。表2河谷形态参数自动提取的GIS技术实现GIS程序模块基于ArcEngine二次开发平台和VisualStudio2010开发环境，采用C#高级语言环境，实现现有GIS软件尚不能完成或者操作复杂的功能。ArcGISEngine是ESRI公司为开发嵌入式GIS和独立运行的GIS桌面端应用程序而推出的二次开发组件库，是当前GIS桌面应用程序开发的主要工具之一。该开发工具包是基于ArcObjects组件的软件开发产品，用于构建自定义GIS和制图应用软件，提供了基于COM、.NET、Java和C++四种开发环境的应用程序，能够实现跨平台部署。本申请利用ArcEngine的相关接口，完成了包括采样点生成、河谷横剖面的水平面投影线生成、高程提取以及各类属性数据的读取、编辑等功能，然后基于采样点属性信息建立数组，按照各类参数计算方法编写程序最终实现了河谷纵、横剖面形态参数的自动提取。实施例1利用本方法提取了中国西藏自治区亚东县康布麻曲及其支流的标准化河长-坡降指数SL’的分布情况，并基于该指数开展了滑坡堵江时间的自动识别应用分析。1、数据源与地貌点SL’提取采用于1:5万地形图数字化DEM(空间分辨率约25m)为基础数据，提取并计算了流域内10条主要河流的均衡坡降K，按50m采样间隔计算了康布麻曲及其支流上所有河谷纵剖面采样点的SL’。由于研究区面积相对较小，水系形态相对规则，本申请暂不考虑流域内气候变化、河流袭夺等因素对区内河流地貌演化的影响，只分析滑坡堵江、基岩变化及构造运动三种因素产生的河流地貌效应。基于ArcGIS通过遥感解译和地质图(1:20万)叠加分析获取了研究区沿河分布的滑坡堵江点16个、岩性界线点(下文称“岩性点”)48个和断层点10个，并将它们作为主要地貌点统计SL’的取值情况。地貌点SL’的计算方式如下：由于滑坡堵江事件、基岩变化及断层活动对河谷地貌的改造具有时空演化性，其原始作用位置常随时间产生偏移，同时考虑DEM的采样不确定性，各地貌点上的SL’值未必能真实反映不同因素对河谷纵剖面形态特征的改造程度。针对这一情况，如图5所示，本申请以地貌点为中心在其上、下游各取3个点，然后将6点中SL’取值最大3点的SL’平均值作为地貌点的SL’。2、数据分析(1)河流均衡坡降K河流均衡坡降K反映河流的整体形态、水动力条件和发育程度。研究区河流最短6.8km，最长52.1km，平均长度20.9km；K值最小194，最大1640，平均1020。经曲线拟合发现区内河流的K与L呈对数关系，如图6所示，反映出麻曲流域河流越长，坡降越大的趋势。总体来说，康布麻曲位于流域上游，河流长度短，坡降大，下蚀作用强烈，在新构造运动剧烈、地壳抬升背景下，尚处于河流发育的早期阶段。(2)计算间隔对S的影响图7是取计算间隔为100m时计算得到的康布麻曲支流麻曲纵横面上各点高程H和坡降S的分布图。可以看出S对高程变化响应敏感：上游H值变化不大，S值在0上下小幅波动；L>20000m后，H骤然降低，S随之产生剧烈波动。图7中河流各点S取值虽基本大于0，仍有占总长8％的河段S<0，即横线以下的河段，与实际情况不符，这一误差主要由DEM精度不足产生。为降低这一误差对结果分析的影响，本文按公差100，取100-1000m等差数列的10个项数为计算间隔，分别计算研究区所有河流采样点的S值，进而对比统计不同计算间隔x对应的S<0的采样点所占百分比y1和S取值区间长度y2。如图8所示，两者都随计算间隔增大逐渐减小。曲线拟合发现两组数据均符合对数函数关系，拟合精度接近1，拟合公式分别见公式2-1、2-2。y1＝21.492-2.979ln(x-13.975)(2-1)y2＝1.685-0.185ln(x-13.975)(2-2)为突出反映不同地貌点的SL'特征，计算间隔取值应该使y1偏小，而y2偏大。图8中，当计算间隔x＝300m时，y1＝4.63％(<5％)，y2＝0.65，综合考虑认为该值作为SL'提取的计算间隔比较理想。(3)不同地貌点SL'对比以300m为计算间隔，本申请分别提取并统计了研究区各滑坡堵江点、岩性点和断层点的SL’，结果见表3；并绘制了各地貌点的SL’盒形图，见图9。图9中矩形盒外侧的短边分别表示上、下十分位数，矩形盒的两边表示上、下四分位数，矩形盒内的横线为中位数，方框为平均值，外部的星点为最大、最小值。可知，滑坡堵江点SL'的各统计量显著大于岩性点和断层点的对应值。可以认为研究区中的三类地貌点，根据对河流地貌的改造程度从强到弱排序为滑坡堵江>岩性改变>构造活动。滑坡堵江点的矩形盒在纵轴方向与岩性点和断层点距离较远，具有类间差异大的特征；而后两者的矩形盒宽度小，类内差异小。以上两点有助于滑坡堵江地貌点的提取和识别。表3三类地貌点SL’统计表3、滑坡堵江识别本申请以岩性点、断层点SL'的上十分位数和滑坡堵江点SL'的下十分位数为判定参数规定了滑坡堵江点的提取规则。如图10所示，遍历研究区所有河流采样点，若采样点SL'≥1.518，则认为该点为滑坡堵江点；若0.906≤SL'＜1.518时，且该点50m范围内无岩性点分布，则该点为滑坡堵江点；若0.844≤SL'＜0.906，且该点50m范围内无岩性点或断层点分布，则该点为滑坡堵江点。不满足上述条件之一的点皆不是滑坡堵江点。若干连续采样点同时被提取为滑坡堵江点时，视为一处滑坡堵江事件。利用上述提取规则，本申请提取了研究区滑坡堵江事件33处。经遥感解译和野外验证，确定研究区实际发生滑坡堵江事件28处。提取的滑坡堵江事件中，正确识别24处，错误识别6处，遗漏4处，识别正确率为24/33＝72.73％，能有效识别研究区85.71％(24/28)的滑坡堵江事件，效果比较理想。分析发现被错误识别的滑坡堵江事件基本位于河流转向或支流汇入的河段，而堵江事件漏提的原因主要有两种：其一是堵江坝体规模小，在水力作用下很快溃坝，对河道地形改造程度低，导致SL’偏小；其二是由于距离相近的多处山体滑坡堵江后，被改造河段彼此相连，被误提为1处。通过亚东县康布麻曲及其支流上滑坡堵江事件的提取实验，可以证明在高山峡谷地区利用河长-坡降指数自动识别具有一定规模的滑坡堵江事件是可行的。实施例2利用本方法以中国喜马拉雅山地区已发生的堵江滑坡和非堵江滑坡事件的地理位置为定位点，自动提取各点所在河谷横剖面的形态参数，基于上述参数建立模型，分析滑坡发生后能否堵江的地形敏感性，实现滑坡体堵江能力评价。1、数据源与河谷横剖面形态参数自动提取河谷横剖面形态参数计算的基础数据包括覆盖中国喜马拉雅山地区的高程数据和堵江滑坡、非堵江滑坡事件定位点。本申请选取30m分辨率的ASTERGDEMv2作为高程数据源。堵江滑坡和非堵江滑坡定位点的地理位置，主要通过遥感影像识别获取。考虑到滑坡事件及气候、人为活动等其他因素对河谷横剖面原始形态的改造，定位点多位选在滑坡附近地形变化较小的河道上，而附近地形变化过大的滑坡点则不予采用。本申请最终获取了中国喜马拉雅山地区90个堵江滑坡事件定位点和91个非堵江滑坡事件定位点。基于上述数据，本申请将河流纵剖面采样间隔和横剖面采样间隔皆设置为40m，自动提取了两类定位点的河谷横剖面形态参数，包括坡高、坡高宽深比、3/4坡高宽深比、1/2坡高宽深比、1/4坡高宽深比、横剖面面积、凹度、不对称系数和河宽指数。河宽指数计算中，为了确定θ，分别提取并统计了tanθ取0.87、0.1、0.15、0.2时两类定位点的河宽指数，表4为不同θ设定条件下堵江滑坡定位点的河宽指数统计表。可见，随着θ的增大，堵江滑坡定位点河宽指数的最小值不变，皆为40m，最大值、最小值和标准差逐渐增大，这反映出定位点河宽指数的取值分布随着θ增大趋向离散，呈现差异性，有助于分析不同类型定位点河谷横剖面形态的特征差异。非堵江滑坡样本的河宽指数也具有相同特点。但当θ过大时，河宽指数与实际河宽相差较大，其受斜坡坡形影响的不确定性也增大，因此θ取值不能过大，去折中选择。考虑arctan0.2≈11.3°，对研究区整体属于高山峡谷的地形而言坡度相对偏低，本申请在河宽指数自动提取中设置tanθ为0.2。表4不同θ设定条件下堵江滑坡定位点河宽指数统计表2、两类定位点河谷横剖面形态参数分布特征对比分析本申请统计了中国喜马拉雅山地区堵江滑坡和非堵江滑坡两类定位点的河谷横剖面形态参数，据此对该区堵江滑坡事件与非堵江滑坡事件发生的河谷横剖面形态差异性进行分析，统计结果见表5。表5河谷横剖面形态参数差异性统计表坡高、各宽深比、凹度、河宽指数的样本差异性统计量Z均大于z0.025＝1.96，剖面面积和不对称系数的差异性统计量Z小于该值，据此可以认为在显著性水平0.05条件下，中国喜马拉雅山地区堵江滑坡与非堵江滑坡河谷横剖面的坡高、宽深比、凹度和河宽指数具有显著差异，其中河宽指数差异最大，3/4坡高宽深比次之。堵江滑坡定位点的宽深比普遍偏小，集中在2.5-4.5之间；非堵江滑坡定位点的宽深比分布较离散，大部分大于4。堵江滑坡定位点处凹度平均值为0.958，属于“V”型谷，所有点凹度皆小于1.5，且标准差较小，分布集中；非堵江滑坡定位点的凹度平均值为1.184，最大值为2.655，属于“U”型谷，大于1的样本点占总数的65.9％，标准差大，分布离散，这反映了“V”型谷比“U”型谷更有利生成滑坡堵江事件。河宽指数的最大值、平均值在两类定位点间差异性更明显。此外，研究区两类定位点的横剖面面积和不对称系数不具有明显差异，这在一定程度上说明这两个参数对滑坡体能否堵江关系不大，坡高、宽深比、凹度、坡高和河宽是有可能影响滑坡体堵将能力的主要地形因素。此外，利用ArcGIS表面分析工具提取了两类定位点发生前原始斜坡的平均坡度，统计得出两者平均坡度的差异性统计量Z为2.070，略大于z0.025，其差异显著性小于本申请获取的坡高、各宽深比、凹度和河宽指数，所以河谷横剖面形态参数更适于参与分析河谷地形与滑坡体堵江能力之间的相关性。3、基于逻辑回归的滑坡体堵江能力评价以宽深比、凹度、河宽指数、坡高为主要地形因子，本申请基于逻辑回归方法建立了滑坡体堵江能力评价模型。(1)指标变量选取认为地形、物源和水力是决定山间滑坡发生后滑坡体运动形成堵江的主要条件。本申请基于中国喜马拉雅山地区堵江滑坡定位点与非堵江滑坡定位点河谷横剖面形态参数对比分析，将宽深比、凹度、河宽指数、坡高作为地形条件指标参与模型构建；物源条件指标选取滑坡发生地层的工程地质岩组反映滑坡体的物质结构，滑坡面积代表滑坡规模；水力条件方面，由于滑坡发生当时的河流流速、流量等详细水文数据都无法获取，本申请暂用上游流域面积来粗略表征。综上本申请共采用7个因子参与构建滑坡体堵江能力评价模型，表6为模型变量列表。表6滑坡体堵江预测建模指标变量列表(2)样本数据获取样本数据质量直接影响模型的精度和可靠性。本申请利用中国喜马拉雅山地区90个堵江滑坡点和91个非堵江滑坡点作为样本，其中前者为1-值样本、后者为0-值样本。如图11所示，宽深比、凹度、河宽指数、坡高等地形指标利用本申请编写的河谷形态参数自动提取模块获取，其他参数利用ArcGIS软件数据处理与空间分析实现，其中工程地质岩组通过样本点与岩组图层叠加提取，滑坡面积通过滑坡几何计算提取，上游流域及其面积通过水文分析获取。(3)模型构建利用SPSS软件，构建了中国喜马拉雅山地区滑坡体堵江能力评价模型。在样本数据中，选出80个1-值样本和80个0-值样本参与模型构建，留出11个1-值样本和10个0-值样本用于模型验证。利用SPSS的Logistic回归分析工具，输入样本数据，选出因变量和变量，设置迭代方式Forward:Conditional，构建研究区滑坡体堵江能力评价模型。迭代过程中，逐渐有变量由于显著性指数sig过大而被排除参与模型构建，仅剩宽深比、凹度、河宽指数、工程地质岩组和滑坡面积成为最终参与滑坡体堵江能力预测的变量。中国喜马拉雅山地区滑坡体堵江能力评价模型可表达为：Z＝-1.466-1.468x11-3.586x2-0.10x3+3.157x55+2.117x6(2.3)各变量的回归系数统计结果见表7。表7模型参数列表通过喜马拉雅山地区滑坡体堵江能力评价模型可以看出，该区影响滑坡体堵江能力的主要因素包括地形因素以及滑坡体的规模和物质组成，这基本符合工程地质学对滑坡堵江成因的普遍认识。水力条件，至少上游流域指标对滑坡体堵江能力的影响并不突出，考虑到现有数据获取和采样水平，认为可以忽略水力条件这一参数。根据表7中各变量回归系数Wald值的大小，发现坡高河谷宽深比对滑坡体堵江预测的重要性最大，其他变量根据重要性从高到低分别是滑坡面积、岩性、河宽指数和凹度。其中宽深比、河宽指数、凹度三个主要的地形因子与滑坡堵江成正相关关系，印证了深切河谷地形条件对滑坡堵江的必要性。作为堵江主体，滑坡体规模越大，岩组成分越坚硬，越容易形成滑坡堵江。在喜马拉雅山地区，最容易形成滑坡堵江的岩组类型是较坚硬岩组。2.4模型检验模型验证包括两部分，一部分是利用SPSS在模型构建过程中进行的验证，包括用于单项因素的计分(S.E)验证、Wale验证以及所建模型的样本预测精度。表7中，S.E是Std.Error的缩写，表示回归系数估计量的标准差。Wald是回归系数检验的统计量值；Sig是Wale检验的显著性概率。回归系数Wald值越大，Sig值越小，说明参数显著性高。表中所有回归系数的的Wale>6，sig<0.05，说明预测模型在显示水平0.05条件下，是可以接受的。在此基础上以P＝0.5作为区分滑坡体堵江与否的分类标准，计算并统计了建模样本和验证样本的预测精度，分别见表8和表9。表8建模样本预测精度表9验证样本预测精度从表9中可以看出，研究区滑坡体堵江能力评价模型对1-型样本和0-型样本的预测精度相差不大，1-型样本预测精度略高，总体精度为86.25％。该模型在验证样本的总体预测精度为80.95％，其中对1-型样本的预测精度大于0-型样本。可见该模型预测结果相对保守，能满足灾害防治中滑坡体堵江概率计算的要求。通过采用本发明公开的上述技术方案，得到了如下有益的效果：(1)面向对象：以沿河流分布的任意一点为对象，自动提取其标准化河长-坡降指数和所在河谷横剖面的形态特征参数，几何意义表征明确。沿河每一个采样点都可以参与计算，能够从微观上表征河谷地貌及分段差异，有利于分析河谷的整体地形。(2)表征全面：河谷纵剖面的标准化河长坡降指数和河谷横剖面的坡高、宽深比、剖面面积、河谷凹度、不对称系数、河宽指数等形态特征参数，共同反映河流的地势变化及河谷整体形态和切割程度等，信息表达全面、完整。(3)干扰因素小，精度更高：本方法所提河谷形态参数基于河谷剖面几何形态进行计算提取，输入数据包括DEM、定位点位置、采样间隔、计算间隔，输入条件少，无假设条件，误差主要来源于DEM质量和空间采样，人为操作干扰程度低，参数计算精度得到提高。(4)运算效率提高：河谷参数计算以沿河分布的定位点为对象，待计算点数目大幅降低，有利于提高计算效率和节约存储空间。同时可以一次性获取所有特征参数的属性，操作简单，具有可维护性。以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视本发明的保护范围。