CN109492194B - 一种基于数学向量几何的dem二阶地形因子计算方法 - Google Patents

Abstract

本发明公开一种基于数学向量几何的DEM二阶地形因子计算方法，该方法首先将坡度或者坡向数据使用坡向极坐标系进行坐标表达，然后将坐标转化为普通极坐标系进行表达，接着再将普通极坐标转化为平面直角坐标系下进行表达，最终化为平面直角坐标系下的向量表达；再依据向量的运算法则实现向量之间的运算，计算得到东西和南北方向上的坡度或坡向变化率，最终求解得到准确的二阶地形因子计算结果。本发明充分考虑到坡度和坡向数据具有的方向属性，在数据表达和计算方式上都依据向量的基本规则，求解出更为科学准确的二阶地形因子结果。

Description

一种基于数学向量几何的DEM二阶地形因子计算方法

技术领域

本发明涉及数字地形分析算法设计领域，具体涉及一种基于数学向量几何的DEM二阶地形因子计算方法。

背景技术

坡度与坡向是两个极为重要的一阶地形因子，是相互联系的两个参数，坡度反映斜坡的倾斜程度，坡向反映斜坡所面对的方向。坡度大小直接影响着地表物质流动与能量转换的规模与强度，是制约生产力空间布局的重要因子。坡向是决定地表面局部地面接收阳光和重新分配太阳辐射量的重要地形因子之一，直接造成局部地区气候特征的差异，同时，也直接影响到诸如土壤水分、地面无霜期以及作物生长适宜性程度等多项重要的农业生产指标。

坡度变率是地面坡度在微分空间的变化率，是依据坡度的求算原理，在所提取的坡度值的基础上对地面每一点再求算一次坡度，即坡度之坡度(slope of slope，SOS)。坡度是地面高程的变化率的求解，因此，坡度变率表征了地表面高程相对于水平面变化的二阶导数。坡向变率是指在地表的坡向提取基础之上，对坡向变化率值的二次提取，亦即坡向之坡度(slope of aspect，SOA)。坡向变率在所提取的地表坡向矩阵的基础上沿袭坡度的求算原理，提取地表局部微小范围内坡向的最大变化情况。坡度变化率、坡向变化率是对地形基本因子坡度和坡向变化情况进行度量的指标，这两个因子在地貌形态结构研究中具有重要的意义。

以坡向变率的求解为例。最初提取SOA是直接采用对坡向数值矩阵求取坡度的方法，这种方法很快被证实在北坡有极大的误差。汤国安等提出基于正反DEM结合的方式消除误差，并得到了广泛的应用。这种方法优点是易于理解，直接在GIS软件中建立模型即可达到快速求解的目的。但是这种方法操作繁琐，求解步骤冗余，且存在着误差特征概括不全面以及高SOA区域计算误差大的缺陷。谢轶群等提出六分法的方式来改进正反地形法存在的误差。这种方法完成了对高SOA区域的数值修正，简化了正反DME法求解SOA的操作步骤。但是六分法算法本身较为复杂，不易理解，且对低SOA区间的数值并没有修正。

另一方面，传统方法求解坡度变率会夸大坡度变化信息。传统方法求解坡度变率同样是参考坡度的计算方法，对坡度数据再求一次坡度。坡度数据不同于坡向数据，其角度大小规定是0到90度，因此不会存在角度差大于180度的情况，也就不会出现如直接求解坡向变率这样显而易见的错误，传统求解方法也可以得到看似满意的结果，然而这样的结果是不科学的。但正是由于这样看似满意的结果，造成了对坡度也具有方向性的忽视。

发明内容

发明目的：针对现有的技术存在的上述问题，提供一种基于数学向量几何的DEM二阶地形因子计算方法。

技术方案：本发明所述的包括以下步骤：S1：获得研究区域的栅格DEM数据，从中提取每一地面点的坡度或坡向值；S2：基于每一坡度或坡向值计算在坡度或坡向极坐标系下每一地面点的坐标向量；S3：将每一地面点的坐标向量从坡度或坡向极坐标系转换为普通极坐标系；S4：将每一地面点的坐标向量从普通极坐标系转换为平面直角坐标系；S5：基于平面直角坐标系中的每一地面点的坐标向量计算东西和南北方向的坡度或坡向变化率；S6：基于东西和南北方向的坡度或坡向变化率计算二阶地形因子。

步骤S2中，坡度或坡向极坐标系是以正北方向为起始方向且以顺时针为旋转方向；且步骤S2通过式(1)实现：

A＝(r,α) (1)

其中，A为每一地面点的坐标在坡度或坡向极坐标系中的表示，r是栅格DEM数据中栅格单元的大小，α是坡向值。

步骤S3通过式(2)和(3)实现：

A'＝(r,f(α)) (2)

其中，A’为每一地面点的坐标在普通极坐标系中的表示。

步骤S4通过式(4)和(5)实现：

A″＝(x,y) (4)

其中，A″为每一地面点的坐标在平面直角坐标系中的表示，β是A′的旋转角度，且β＝f(α)。

步骤S5中，通过式(6)计算东西和南北方向的坡度或坡向变化率：

其中，h代表中心栅格单元的坡度或坡向向量，i和j分别代表中心栅格单元的行号和列号，Δx代表东西方向坡度或坡向变化率，Δy代表的是南北方向坡度或坡向变化率。

步骤S6中，通过式(6)计算二阶地形因子：

其中，SOS为坡度变率，SOA为坡向变率。

有益效果：本发明与现有技术相比，其优点为：(1)能够有效地将坡度或坡向的方向属性用于计算，避免了传统标量方法无视方向属性造成坡向变率计算结果夸大的误差；(2)弱化的坡度或坡向变率计算结果以及显著的频率变化分布特征在相当程度上是更加利于对坡度或坡向突变区域的提取；(3)坡度或坡向变率计算结果分布特征更加符合地表空间中坡面转折的渐变、突变、以及不变的地貌形态分布的基本认知，即不变及渐变的地形区域占据主要地位。因此，本发明所提出的基于向量几何的二阶地形因子计算方法为精准数字地形分析提供参考。

附图说明

图1(a)和1(b)分别为向量做差运算和向量视角下的标量做差运算；

图2(a)和2(b)分别为坡度或坡向极坐标系和普通极坐标系；

图3为平面直角坐标系下向量表示特征；

图4为本发明的基于向量几何法的二阶地形因子计算方法算法流程图；

图5为本发明实施例中采用的实验样区与DEM数据；

图6(a)至6(c)分别为对延安样区采用正反DEM法、六分法和本发明的向量几何法计算坡向变率结果；

图7为基于延安样区的不同方法坡向变率计算结果频率分布；

图8(a)和8(b)为延安样区标量方法和向量方法计算坡向变率结果；

图9为基于延安样区的不同方法坡度变率计算结果频率分布。

具体实施方式

以下通过一个具体的实施例并结合附图对本发明进行详细说明。

如图4，本发明的基于数学向量几何的DEM二阶地形因子计算方法包括如下步骤：

步骤一：获得研究区域较高分辨率栅格DEM数据，可以在ArcGIS软件中进行坡度或者坡向提取，作为计算二阶地形因子的基础数据。

延安市位于黄河中游，属黄土高原丘陵沟壑区。延安地貌以黄土高原、丘陵为主。地势西北高东南低，平均海拔1200米左右。北部的白于山海拔1600－1800米，最高点在吴旗县五谷城乡的白于山顶，海拔1809.8米；最低点在宜川县集义乡猴儿川，海拔388.8米，相对高差1421米。北部以黄土梁峁、沟壑为主，占全区总面积72％；南部以黄土塬沟壑为主，占总面积19％；全区石质山地占总面积9％。西部子午岭，南北走向，构成洛河与泾河的分水岭，是高出黄土高原的基岩山地之一，海拔1500－1600米，主峰1687米；东部黄龙山，大致呈南北方向延伸，海拔1500米，主峰(大岭)海拔1788.7米；中部劳山，呈西北——东南走向，平均海拔1400米，主峰(大墩梁)海拔1464米。黄龙山和劳山统称为梁山山脉，形成延安地区地形的骨架。由于样区地形复杂，因此选择陕西省延安市实验样区5m分辨率DEM数据作为原始数据为例进行说明，样区位置如图5所示。对该实验样区首先使用ArcGIS进行坡度、坡向提取，获得计算二阶地形因子所需要的实验数据准备。

步骤二：基于每一坡度或坡向值计算在坡度或坡向极坐标系下每一地面点的坐标向量。

如图2(a)所示，在坡度或坡向极坐标系下，每一地面点的坐标向量可以认为处在一个以正北方向为起始方向，顺时针为旋转方向的特殊极坐标系下。坡度或坡向值只有首先转为坡向极坐标系统，使用坐标表达，才能进一步转为向量，其表示为：

A＝(r,α) (1)

其中r是长度，这里是栅格DEM数据中栅格单元的大小；α是旋转角度，这里是坡度或坡向值。

步骤三：将每一地面点的坐标向量从坡度或坡向极坐标系转换为普通极坐标系。

如图2(b)所示，普通极坐标系是较为常用的一种坐标系，并且可以与平面直角坐标系进行快速转换。

从坡度或坡向极坐标系到普通极坐标系的转换可以通过式(2)和(3)实现：

A'＝(r,f(α)) (2)

其中，A’为每一地面点的坐标在普通极坐标系中的表示。

步骤四：将每一地面点的坐标向量从普通极坐标系转换为平面直角坐标系。关于直角坐标，如图3所示，分别取与X轴、Y轴方向相同的两个单位向量

作为基底。任作一个向量/>

由平面向量基本定理可知，有且只有一对实数x、y，可以使得等式：/>

成立，于是把(x,y)叫做向量/>

的直角坐标表示，记作：/>

平面直角坐标系更适合于向量的表示，因此需要将极坐标系下的坡向坐标转到平面直角坐标系下。这样一来，就可以运用平面向量基本定理，使用向量的平面直角坐标表达法来表达坡度或坡向数据。采用式(4)和(5)可以实现从普通极坐标系到平面直角坐标系的转换：

A″＝(x,y) (4)

其中，A″为每一地面点的坐标在平面直角坐标系中的表示，β是A′的旋转角度，且β＝f(α)。

步骤五：基于平面直角坐标系中的每一地面点的坐标向量计算东西和南北方向的坡度或坡向变化率。

坡度变化率、坡向变化率是对地形基本因子坡度和坡向变化情况进行度量的指标，坡度变化率是地形单元中坡度变化的描述或者说坡度的坡度，而坡向变化率则是坡向变化程度的量化表达即坡向的坡度。这里重要的是变化率的计算。本发明使用三阶反距离平方权差分算法计算坡度或坡向在东西和南北方向上的变化率，如式(6)所示：

其中，h代表中心栅格单元的坡度或坡向向量，i和j分别代表中心栅格单元的行号和列号，Δx代表东西方向坡度或坡向变化率，Δy代表的是南北方向坡度或坡向变化率。

需要注意的是，这里应当使用向量的运算法则进行公式运算。

二阶地形因子的计算中很重要的部分是邻域栅格单元内格网数值的代数做差运算。由于传统方法进行差分运算使用的是未经向量化转换的坡度或坡向数值数据，运算法则使用的是标量之间的代数运算法则，忽视了方向对运算的影响。该方法得到的结果与地表真实的坡度或坡向变率存在偏差。在平面直角坐标系中，传统方法是将坡向矩阵中的每个数值看作是二维空间中起点在Y轴上，方向指向X正半轴，模长为坡度或坡向值的向量。由于方向均指向X正半轴，此时，向量之间的运算可以使用标量的代数运算法则代替。但是，该运算方法使得向量之间的运算丧失二维空间意义，在方向相同时，该运算实质仍是一维标量运算(图1(b))。

在本发明中，坡度或坡向隐式表达的坡度或坡向数值矩阵数据将进行向量化表达，坡度或坡向矩阵中每个坡度或坡向值都有其方向。此时，坡度或坡向做差运算则不是一维标量运算，而是符合坡向定义的二维向量运算(如图1(a))。在平面直角坐标系下，平面向量间的运算法则如下：

已知

则/>

已知

则/>

已知

和实数λ,则/>

已知

则/>

的模长是/>

通过以上四个向量之间的运算法则，可完整、准确地将三阶反距离平方权差分算法采用向量运算法则进行计算，并应用于二阶地形因子的计算。

步骤六：基于东西和南北方向的坡度或坡向变化率计算二阶地形因子。将步骤五中计算得到的结果代入式(6)，即可计算得到二阶地形因子SOS或SOA：

其中，SOS为坡度变率，SOA为坡向变率。

针对坡向变率的特殊性，同时采用正反DEM法、六分法以及向量几何法计算坡向变率，并对所计算结果进行定量分析。坡向变率计算结果如图6(a)至(c)所示。可以看出，正反DEM法和六分法的结果差异相对较小，如图6(a)和(b)。而本文所提向量法得到结果明显不同于传统的标量计算方式结果。具体表现为，坡向变率值较高的区域呈现出较强的空间结构特征，而坡向变化平缓区域的坡向变率值较小(图6(c))。可以看出，在不考虑坡向特征在数学上存在向量方向特征时，传统的基于数学标量规则的坡向变率计算方法在相当程度上夸大了地表面的坡向变化特征，使得地表空间中的不变和渐变的地形特征所表现出来的坡向变化难以得到合理的结果。

图7为延安样区采用六分法、正反DEM法、以及向量法计算得到的坡向变率频率分布结果。综合图6和图7可知，对于黄土地貌类型区域，山脊、山谷、山顶等地形特征部位在样区中表现为地形骨架信息，处于量少而又重要的地位。这些地形特征部位坡向变化较大，属于地表空间地形变化的突变区域。因此，相应的坡向变率计算结果较大。而黄土地貌类型区域地形中大部分的区域仍然是坡向转折较小，坡面形态改变相对平缓的区域。因此，地形特征仍属于渐变区域，相应的坡向变率计算结果较小。向量法得到的坡向变率结果在相当程度上更好的符合黄土地貌地貌坡向变化的认知，坡向变率低值区域占大部分，坡向变率高值区域占小部分。而基于标量运算的六分法和正反DEM法恰恰相反，高坡向变率值栅格频率远高于低坡向变率值栅格频率，在相当程度上夸大了地表坡面转变的信息。

对于坡度变率，恰恰由于其没有坡向变率那种特殊的误差，才使得人们忽略了方向属性对二阶地形因子求解的重要性。然而对于复杂地形，如延安实验样区，坡度变率求解过程中忽视方向属性仍会造成极大的误差，这将会对以坡度变率为基础的后续研究带来难以估计的严重后果。图8给出了延安样区采用传统标量法、向量几何法的计算结果展示，可以看出对于坡度变率的求解，标量运算同样存在夸大地表坡度转换信息的错误。图9为采用传统标量法、向量几何法的坡度变率频率分布结果，可以发现传统方法确实出现了极大的偏差。

Claims (6)

1.一种基于数学向量几何的DEM二阶地形因子计算方法，其特征在于，包括以下步骤：

S1：获得研究区域的栅格DEM数据，从中提取每一地面点的坡度或坡向值；

S2：基于每一坡度或坡向值计算在坡度或坡向极坐标系下每一地面点的坐标向量；

S3：将每一地面点的坐标向量从坡度或坡向极坐标系转换为普通极坐标系；

S4：将每一地面点的坐标向量从普通极坐标系转换为平面直角坐标系；

S5：基于平面直角坐标系中的每一地面点的坐标向量计算东西和南北方向的坡度或坡向变化率；

S6：基于东西和南北方向的坡度或坡向变化率计算二阶地形因子。

2.根据权利要求1所述的DEM二阶地形因子计算方法，其特征在于，步骤S2中，坡度或坡向极坐标系是以正北方向为起始方向且以顺时针为旋转方向；且步骤S2通过式(1)实现：

A＝(r，α) (1)

其中，A为每一地面点的坐标在坡度或坡向极坐标系中的表示，r是栅格DEM数据中栅格单元的大小，α是坡向值。

3.根据权利要求2所述的DEM二阶地形因子计算方法，其特征在于，步骤S3通过式(2)和(3)实现：

A′＝(r，f(α)) (2)

其中，A′为每一地面点的坐标在普通极坐标系中的表示。

4.根据权利要求3所述的DEM二阶地形因子计算方法，其特征在于，步骤S4通过式(4)和(5)实现：

A″＝(x，y) (4)

其中，A″为每一地面点的坐标在平面直角坐标系中的表示，β是A′的旋转角度，且β＝f(α)。

5.根据权利要求4所述的DEM二阶地形因子计算方法，其特征在于，步骤S5中，通过式(6)计算东西和南北方向的坡度或坡向变化率：

其中，h代表中心栅格单元的坡度或坡向向量，i和j分别代表中心栅格单元的行号和列号，Δx代表东西方向坡度或坡向变化率，Δy代表的是南北方向坡度或坡向变化率。

6.根据权利要求5所述的DEM二阶地形因子计算方法，其特征在于，步骤S6中，通过式(6)计算二阶地形因子：

其中，SOS为坡度变率，SOA为坡向变率。

Images