9.5. DEM 的分析和应用

地形曲面拟合 #

DEM 最基础的应用是求 DEM 范围内任意点的高程，在此基础上进行地形属性分析。 由于已知有限个格网点的高程，可以利用这些格网点高程拟合一个地形曲面，推求区域内任意点的高程。 曲面拟合方法可以看作是一个已知规则格网点数据进行空间插值的特例，距离倒数加权平均方法，克里金插值方法， 样条函数等插值方法均可采用。

立体透视图 #

从数字高程模型绘制透视立体图是 DEM 的一个极其重要的应用。 透视立体图能更好地反映地形的立体形态，非常直观。 与采用等高线表示地形形态相比有其自身独特的优点，更接近人们的直观视觉。 特别是随着计算机图形处理工作的增强以及屏幕显示系统的发展，使立体图形的制作具有更大的灵活性， 人们可以根据不同的需要，对于同一个地形形态作各种不同的立体显示。 例如局部放大，改变高程值Z的放大倍率以夸大立体形态； 改变视点的位置以便从不同的角度进行观察，甚至可以使立体图形转动，使人们更好地研究地形的空间形态。

从一个空间三维的立体的数字高程模型到一个平面的二维透视图，其本质就是一个透视变换。 将“视点”看作为“摄影中心”，可以直接应用共线方程从物点 （X，Y， Z） 计算“像点”坐标 （X,Y） 。 透视图中的另一个问题是“消隐”的问题，即处理前景挡后景的问题。

调整视点、视角等各个参数值，就可从不同方位、不同距离绘制形态各不相同的透视图制作动画。 计算机速度充分高时，就可实时地产生动画 DTM 透视图。

通视分析 #

通视分析有着广泛的应用背景。 典型的例子是观察哨所的设定，显然观察哨的位置应该设在能监视某一感兴趣的区域，视线不能被地形挡住。 这就是通视分析中典型的点对区域的通视问题。 与此类似的问题还有森林中火灾监测点的设定，无线发射塔的设定等。 有时还可能对不可见区域进行分析，如低空侦察飞机在飞行时，要尽可能躲避敌方雷达的捕捉， 飞行显然要选择雷达盲区飞行。 通视问题可以分为五类[Lee,J.（1991）]：

1. 已知一个或一组观察点，找出某一地形的可见区域。

2. 欲观察到某一区域的全部地形表面，计算最少观察点数量。

3. 在观察点数量一定的前提下，计算能获得的最大观察区域。

4. 以最小代价建造观察塔，要求全部区域可见。

5. 在给定建造代价的前提下，求最大可见区。

根据问题输出维数的不同，通视可分为点的通视，线的通视和面的通视。 点的通视是指计算视点与待判定点之间的可见性问题； 线的通视是指已知视点，计算视点的视野问题； 区域的通视是指已知视点，计算视点能可视的地形表面区域集合的问题。 基于格网 DEM 模型与基于TIN模型的 DEM 计算通视的方法差异很大。

1）点对点通视

基于格网 DEM 的通视问题，为了简化问题，可以将格网点作为计算单位。 这样点对点的通视问题简化为离散空间直线与某一地形剖面线的相交问题。 （图9-13）

已知视点V的坐标为（ x_0,y_0,z_0 ）， 以及P点的坐标（ x_1,y_1,z_1 ）。 DEM 为二维数组 Z[M][N] ，则 V 为（ m_0,n_0, Z[m_0,n_0] ），P为（ m_1,n_1,Z[m_1,n_1] ）。计算过程如下：

（1．1）使用Bresenham直线算法，生成V到P的投影直线点集 {x , y}，K=||{x ,y}|| ,并得到直线点集 {x, y} 对应的高程数据 {Z[k], ( k=1,...K-1)} ， 这样形成 V 到 P 的 DEM 剖面曲线。

（1．2）以V到P的投影直线为X轴，V的投影点为原点，求出视线在X-Z坐标系的直线方程：

（0<k<K）

K``为 ``V 到 P 投影直线上离散点数量。

（1．3）比较数组 H[k] 与数组 Z[k] 中对应元素的值，如果

存在 Z[k]>H[k] ，则V与P不可见，否则可见。

2）点对线通视

点对线的通视，实际上就是求点的视野。 应该注意的是，对于视野线之外的任何一个地形表面上的点都是不可见的，但在视野线内的点有可能可见， 也可能不可见。 基于格网 DEM 点对线的通视算法如下：

（2．1）设P点为一沿着 DEM 数据边缘顺时针移动的点，与计算点对点的通视相仿，求出视点到P点投影直线上点集 {x,y} ，并求出相应的地形剖面 {x, y, Z(x, y)} 。

（2．2）计算视点至每个与 Z 轴的夹角：

（2．3）求得。对应的点就为视点视野线的一个点。

（2．4）移动P点，重复以上过程，直至P点回到初始位置，算法结束。

3）点对区域通视

点对区域的通视算法是点对点算法的扩展。 与点到线通视问题相同，P点沿数据边缘顺时针移动。 逐点检查视点至P点的直线上的点是否通视。 一个改进的算法思想是，视点到P点的视线遮挡点，最有可能是地形剖面线上高程最大的点。 因此，可以将剖面线上的点按高程值进行排序，按降序依次检查排序后每个点是否通视， 只要有一个点不满足通视条件，其余点不再检查。 点对区域的通视实质仍是点对点的通视，只是增加了排序过程。

流域特征地貌提取与地形自动分割 #

地形因素是影响流域地貌、水文、生物等过程的重要因子， 地形属性的空间分布特征一直是人们用于描述这些空间过程变化的重要指标。 高精度 DEM 数据和高分辨率、高光谱、多周期的遥感影像， 为人们定量描述流域空间变化过程提供了日益丰富的数据源， 而且人们对流域地貌、水文和生物等过程空间变化机理理解的不断加深， 可以说人类已经进入了一个“ 空间模拟 ”的时代。 基于 DEM 数据自动提取流域地貌特征和进行流域地形自动分割是进行流域 空间模拟 的基础技术。

基于格网 DEM 自动提取流域特征地貌和进行地形自动分割技术主要包括两个方面： 1）流域地貌形态结构定义，定义能反映流域结构的特征地貌，建立格网 DEM 对应的微地貌特征。 2）特征地貌自动提取和地形自动分割算法。格网 DEM 数据是一些离散的高程点数据，每个数据本身不能反映实际地表的复杂性。 为了从格网 DEM 数据中得到流域地貌形态结构，必须采用一个清晰的流域地貌结构模型， 然后针对该结构模型设计自动提取算法。

1）流域结构定义

可以使用一个具有根的树状图来描述流域结构[Shreve]，目前绝大多数算法都沿用这一描述方法。 在此结构中主要包括三个部分，即结点集、界线集和汇流区集。 如图9-14所示。

其具体内容包括几个概念：

1. 沟谷线段：一条具有两侧汇流区的线段；

2. 分水线段：一条具有两侧分水区的线段；

3. 沟谷结点：两条或两条以上沟谷线的交点；

4. 分水线结点：两条或两条以上分水线的交点；

5. 沟谷源点：沟谷的上游起点；

6. 分水线源点：分水线与流域边界的交点；

7. 内部汇流区：汇流区边界不包含流域部分边界的汇流区；

8. 外部汇流区：汇流区边界包括部分流域边界的汇流区。

沟谷结点和沟谷源点共同组成沟谷结点集，所有的沟谷段组成沟谷段集，形成沟谷网络； 所有的分水线组成分水线段集，形成分水线网络。 沟谷段集和分水线段集共同把流域分割成一个汇流区集。

沟谷段是最小的沟谷单位，沟谷段可以分为内部沟谷段和外部沟谷段。 内部沟谷段连接两个沟谷结点，外部沟谷段连接一个沟谷结点和沟谷源点。 同样，分水线段是最小的分水线单位，也分为内部分水线段和外部分水线段。 内部分水线段连接两个分水线结点，外部分水线段连接一个分水线结点和一个分水线源点。

汇流网络中每一沟谷段都有一个汇流区域，这些区域由分水线集控制。 外部沟谷段有一个外部汇流区，内部沟谷段有两个内部汇流区，分布在内部沟谷段两侧。 整个流域被分割成一个个子流域，每个子流域如同树状图上的一片“叶子”。

2）流域特征地貌自动提取和地形自动分割

特征地貌定义与提取： 根据网格点高程与周围高程值的关系，将格网点分为坡地、洼地、分水线、谷地、阶地和鞍部等几类。 先计算中心点与八邻点的高程差，然后对高程差进行排序，再根据高程差序列的特性给中心点格网赋一个特征编码。 然后通过一系列特征码的组合特征，用模式识别的方法，将格网点划分到已知的特征地貌类别。

山脊线和山谷线提取：山脊线和山谷线的自动探测实际上是凹点和凸点的自动搜索。 较为简单的算子是 2*2 的局部算子。 将算子在 DEM 数据中滑动，比较每个格网点与行和列上相邻格网点的高程， 标出其中高程最小（探测山谷线）或高程最大（探测山脊线）的格网点。 对整个 DEM 数据计算一遍后，剩下的未标记格网点就是山脊线或山谷线上的格网点。

流域地形自动分割：流域地形自动分割的目标是将整个流域分割成一个个子汇流区。 大多数算法是利用 3*3 窗口计算流向和基于“溢流跟踪”算法确定汇流网络。 算法过程如下：

（2．1）格网点流向定义

采用 3×3 窗口按8方向搜索计算最大坡向为各网格点的流向。 分别为8方向赋不同的代码，如右图所示。 每个格网有一个从1到9的数值，代表它流向相邻象元的方向，如该象元为凹点，则其值为5（图）。

几种例外情况的处理：

1. 如果一个网格点的最大坡向格网点与之具有相同的高程值，且之前没有其它格网点流向这个相邻格网，则强制流向它。如果还有另外的格网点流向这个相邻格网，则当前格网点为凹点。

2. 当两个或多个相邻格网点的最大坡向相等时，先比较各自相邻格网点坡向，如果仍没解决，继续比较相对格网点的坡向，决定赋一个流向。

3. 对于具有相同高程值的区域则扩大搜索窗口半径，用7×7窗口，如果需要还可以使用更大窗口。

4. 在 DEM 数据的外围加一圈高程值为0的格网点，强制其最大坡向流向研究区之外。

当所有的格网点处理完毕后，生成一个编码1—9的流向图。

（2．2）凹点处理算法

由于凹点的存在，有一些流路不会流向流域出口，而是终止于凹点，所以在进行流域自动分割之前， 还要对凹点进行处理。 流域中凹点既可能是真实的凹点，也可能是由于插值误差造成的，所以不能使用简单的滤波或平滑函数， 将凹点全部去除，目的是将凹点造成的断路连接到主沟谷网络。

搜索所有凹点的相邻最低点（有时可能有多个高程相等的最低点），作为凹点的溢出点， 以溢出点为起点继续搜索比它的高程低或相等的邻点（已经搜索的点忽略），判断是否有比原凹点更低的格网点， 如果没有则以该凹点的溢出点为起点，重复上述搜索过程； 如果搜索到比原凹点低的格网点，将凹点和最低邻点的方向倒转。 如图9-15所示： 高程为48的点为一个凹点，搜索到高程最低的邻点为49，以它为起点继续搜索，找到高程点49， 仍比原凹点高程高，则继续搜索，又找到另一个高程点49，再找到高程47的点，比原凹点高程低，结束搜索， 按搜索方向修改流向，如图9-16中实线箭头方向所示。

（2．3）提取汇流网络

根据修改后的流向图，给定一个点，所有流向它的格网点的总和就是该点的汇流区。 计算方法是给定一个点，搜索8邻点，记录所有流向它的格网点的位置， 然后再以找到的格网点为基点继续搜索记录流向它的格网点，直到没有新的汇流点为止， 所有记录的格网点构成该点的汇流区。

通常沟谷的汇流区面积大于其它格网点的汇流区面积，可以通过设定一个阈值，将汇流区面积大于此阈值的格网点， 标识为沟谷点。 很明显，不同的阈值得到的沟谷网络的复杂性是不同的，这种方法虽然为确定沟谷网络的复杂性提供了灵活性， 但也使得沟谷网络的确定具有太大的随意性。

得到沟谷网络后，可以对沟谷网络进行编码。 首先对沟谷结点编码。 从流域出口开始搜索遍历整个汇流网络，对每个沟谷段的上下游结点进行编码标识，标识值是沟谷段的编码值， 并记录下这些结点的位置。 其次，把沟谷段中的每个格网点标识为沟谷段的编码值。 第三，根据沟谷段上游结点的类型判定沟谷段是内部沟谷段还是外部沟谷段。

（2．4）提取分水网络

递归搜索沟谷段中的每个格网点的汇流区，将汇流区的格网点赋为该沟谷段的标识值，形成各沟谷段的子汇流区。 然后进行边界跟踪，提取子汇流区的边界线为分水线，得到分水线网络。 最后，对沟谷网络和分水线网络及子汇流区进行拓扑编码，以完成流域地形的自动分割。

DEM 计算地形属性 #

由 DEM 派生的地形属性数据可以分为单要素属性和复合属性二种。 前者可由高程数据直接计算得到，如坡度因子，坡向。 后者是由几个单要素属性按一定关系组合成的复合指标，用于描述某种过程的空间变化， 这种组合关系通常是经验关系，也可以使用简化的自然过程机理模型。

单要素地形属性通常可以很容易地使用计算机程序计算得到，包括：

1）坡度、坡向

坡度定义为水平面与局部地表之间的正切值。 它包含两个成分：斜度——高度变化的最大值比率（常称为坡度）； 坡向——变化比率最大值的方向。 地貌分析还可能用到二阶差分凹率和凸率。 比较通用的度量方法是： 斜度用百分比度量，坡向按从正北方向起算的角度测量，凸度按单位距离内斜度的度数测量。

坡度和坡向的计算通常使用 3*3 窗口，窗口在 DEM 高程矩阵中连续移动后， 完成整幅图的计算。 坡度的计算如下：

坡向计算如下：

为了提高计算速度和精度，GIS通常使用二阶差分计算坡度和坡向， 最简单的有限二阶差分法是按下式计算点 i,j 在 x 方向上的斜度：

式中是格网间距（沿对角线时应乘以）。 这种方法计算八各方向的斜度，运算速度也快得多。 但地面高程得局部误差将引起严重得坡度计算误差，可以用数字分析方法来得到更好得结果， 用数字分析方法计算东西方向得坡度公式如下：

同理可以写出其它方向的坡度计算公式。

2）面积、体积

（2．1）剖面积

根据工程设计的线路，可计算其与 DEM 各格网边交点 P_i (X_i，Y_i，Z_i) ， 则线路剖面积为

其中n为交点数； D_i，i+1 为 P_i 与 P_i+1 之距离。 同理可计算任意横断面及其面积。

（2．2）体积

DEM 体积由四棱柱（无特征的格网）与三棱柱体积进行累加得到，四棱柱体上表面用抛物双曲面拟合，三棱柱体上表面用斜平面拟合，下表面均为水平面或参考平面，计算公式分别为

其中 S_3 与 S_4 分别是三棱柱与四棱柱的底面积。

根据两个 DEM 可计算工程中的挖方、填方及土壤流失量。

3）表面积

对于含有特征的格网，将其分解成三角形，对于无特征的格网，可由4个角点的高程取平均即中心点高程， 然后将格网分成4个三角形。 由每一三角形的三个角点坐标 (x_i，y_i， z_i) 计算出通过该三个顶点的斜面内三角形的面积，最后累加就得到了实地的表面积。

三角网内插 #

在建立TIN后，可以由TIN解求该区域内任意一点的高程。 TIN的内插与矩形格网的内插有不同的特点，其用于内插的点的检索比网格的检索要复杂。 一般情况下仅用线性内插，即三角形三点确定的斜平面作为地表面，因而仅能保证地面连续而不能保证光滑。 进行三角网内插，一般要经过以下几个步骤：

1）格网点的检索

给定一点的平面坐标 P(x,y) ，要基于TIN内插该点的高程Z， 首先要确定点P落在TIN的哪个三角形中。 一般的做法是通过计算距离，得到据P点最近的点，设为 Q_1 。 然后就要确定P所在的三角形。 依次取出 Q_1 为顶点的三角形，判断P是否位于该三角形内。 可利用P是否与该三角形每一顶点均在该顶点所对边的同侧（点的坐标分别代人该边直线方程所得的值符号相同）加以判断。 若P不在以 Q_1 为顶点的任意一个三角形中，则取离P次最近的格网点，重复上述处理， 直至取出P所在的三角形，即检索到用于内插P点高程的三个格网点。

2）高程内插

若 P（x，y） 所在的三角形为 ΔQ_1 Q_2 Q_3 ，三顶点坐标为 （x_1， y_1，z_1） ， （x_2，y_2，z_2） 与 （x_3，y_3，z_3 ） ， 则由 Q_1，Q_2 与Q_3 确定的平面方程为

或

令

则P点高程为

等高线追踪 #

基于TIN绘制等高线直接利用原始观测数据，避免了 DTM 内插的精度损失，因而等高线精度较高； 对高程注记点附近的较短封闭等高线也能绘制； 绘制的等高线分布在采样区域内而并不要求采样区域有规则四边形边界。 而同一高程的等高线只穿过一个三角形最多一次，因而程序设计也较简单。 但是，由于TIN的存贮结构不同，等高线的具体跟踪算法跟踪也有所不同。

基于三角形搜索的等高线绘制算法如下：

对于记录了三角形表的TIN，按记录的三角形顺序搜索。 其基本过程如下：

1. 对给定的等高线高程h，与所有网点高程z_i （i=1,2,…，n），进行比较，若z_i =h，则将z_i 加上（或减）一个微小正数ε>0（如ε=10_-4`），以使程序设计简单而又不影响等高线的精度。

2. 设立三角形标志数组，其初始值为零，每一元素与一个三角形对应，凡处理过的三角形将标志置为1，以后不再处理，直至等高线高程改变。

3. 按顺序判断每一个三角形的三边中的两条边是否有等高线穿过。若三角形一边的两端点为P_1 （x_1 ，y_1 ，z_1 ），P_2 （x_2 ，y_2 ，z_2 ）则(z_1-h)(z:_2-h)<0表明该边有等高线点；(z_1-h)(z_2-h)>0表明该边无等高线点。

直至搜索到等高线与网边的第一个交点，称该点为搜索起点， 也是当前三角形的等高线进入边、线性内插该点的平面坐标（x，y）：

1. 搜索该等高线在该三角形的离去边，也就是相邻三角形的进人边，并内插其平面坐标。搜索与内插方法与上面的搜索起点相同，不同的只是仅对该三角形的另两边作处理。

2. 进入相邻三角形，重复第（4）步，直至离去边没有相邻三角形（此时等高线为开曲线）或相邻三角形即搜索起点所在的三角形（此时等高线为闭曲线）时为止。

3. 对于开曲线，将已搜索到的等高线点顺序倒过来，并回到搜索起点向另一方向搜索，直至到达边界（即离去边没有相邻三角形）。

4. 当一条等高线全部跟踪完后，将其光滑输出，方法与前面所述矩形格网等高线的绘制相同。然后继续三角形的搜索，直至全部三角形处理完，再改变等高线高程，重复以上过程，直到完成全部等高线的绘制为止，图中描述了利用三角网生成数值为50的等高线的过程。