数字地面模型的应用.doc

第五章 数字地面模型的应用数字地面模型的应用是很广泛的。在测绘中可用于绘制等高线、坡度、坡向图、立体透视图，制作正射影像图、立体景观图、立体匹配片、立体地形模型及地图的修测。在各种工程中可用于体积、面积的计算，各种剖面图的绘制及线路的设计。军事上可用于导航(包括导弹与飞机的导航)、通讯、作战任务的计划等。在遥感中可作为分类的辅助数据。在环境与规划中可用于土地利用现状的分析、各种规划及洪水险情预报等。本章重点介绍数字地面模型在测绘中的应用。151 基于矩形格网的DEM多项式内插 DEM最基础的应用(也是各种应用的基础)是求DEM范围内任意一点P(X，Y)的高程。由于此时已知该点所在的DTM格网各个角点的高程，因此可利用这些格网点高程拟合定的曲面，然后计算该点的高程。所拟合的曲面一般应满足连续乃至光滑的条件。一、双线性多项式(双曲面)内插根据最邻近的4个数据点，可确定一个双线性多项式或用矩阵形式表示为 (152)利用4个已知数据点求出4个系数，然后根据待定点的坐标(X，Y)与求出的系数内插出该点的高程。双线性多项式的特点是：当坐标X (或Y)为常数时，高程Z与坐标Y (或X)成线性关系，故称其为“双线性”。 当数据点规则排列组成矩形或正方形时，可导出直接计算的内插公式。例如当4个已知点组成矩形时(如图151)，可由双线性多项式的定义、直积运算或阵列代数推导出内插公式。1矩阵直积法由P点所在格网4个角点高程Zij(i，j0,1)，列出关于系数的方程或写为 因此由于其中为矩阵的直积运算符号。根据直积的运算规则(参考附录四)：将四阶矩阵的求逆转变为两个二阶矩阵的求逆。由X1- X0= Y1- Y0L可得出内插公式：(158)2阵列代数法 阵列代数是一个新的有力的数学工具，用以扩展线性代数，解决多级数据问题。当原始数据点呈规则格网排列，且方程的系数矩阵能够进行直积分解成x和y方向的两个矩阵时，则矩阵形式的方程可变换为阵列代数形式的阵列方程。由于矩阵F中的一个行阵可按下式分解： 则因此可列出阵列方程：或因此解算的待定参数为将A代人式(152)整理即得双线性内插公式(158)。这两种方法与直接求逆的工作量之比近似地为，当待定参数数目以及维数增多时，计算工作量的节省更加显著。 双线性多项式内插只能保证相邻区域接边处的连续，不能保证光滑。但因其计算量较小，因而是最常用的方法。二、双三次多项式(三次曲面内插)三次曲面方程为因此三次曲面方程可以利用阵列代数方式表达为 若数据点呈方格网分布(图152)，糖坐标原点平移至待定点P所在方格网之左下角则P点坐标(x，y)满足0xL，0yL，其中L为格网边长。为了简单令L=1，则0x1，0y1由于待定系数共有16个，因而除了P所在格网四顶点高程外，还需要已知其点处的一阶偏导数与二阶混合导数，其值可按下式计算：因为且xi-1y j-1-1，xiyi0，xi+1y i+1=1，xi+2y i+2=2，所以不难求出x与y之逆： 从而可求得16个待定的参数值为三次多项式内插虽然属于局部函数内插，即在每一个方格网内拟合一个三次曲面，但由于考虑了一阶偏导数与二阶混合导数，因而它能保证相邻曲面之间的连续与光滑。 设相邻两格网曲面函数为，公共格网边为AB(图152)，相应的参数短阵为，因此则在接边线AB上任意一点(在左右格网平面坐标分别为(1，y)与(0，y)在左右曲面的函数值分别为因为在两个端点A，B处的高程及y方向的斜率相等，即，可得因此，即分界线AB上两相邻曲面的连续性得证。同理，利用A，B上x方向的斜率与扭曲(二阶混合导数)对应相等可证明边界线朋上两曲面沿x方向的斜率均相等、所以是光滑的。若采用不完整的三次多项式内插，即令，解算时只考虑四角数据点的高程与一阶导数。由于未考虑扭曲条件，根据以上证明可知，它只能保证相邻曲面在分界线上的连续，但不能保证光滑。152 等高线的绘制根据规则格网DTM自动绘制等高线，主要包括以下两个步骤：利用DTM的矩形格网点的高程内插出格网边上的等高线点，并将这些等高线点按顺序排列(即等高线的跟踪)；利用这些顺序排列的等高线点的平面坐标X，Y进行插补，即进一步加密等高线点并绘制成光滑的曲线(即等高线的光滑)。一、等高线跟踪 在数字地面模型格网边上内插并排列等高线点的方法很多但总的来说可以分为两种方式：一是对每条等高线边内插边排序；另一种方式是对同一高程的等高线先内插出所有等高线点，再逐一排列每条等高线的点。 1按每条等高线的走向顺序插点 这是一种按逐条等高线的走向边搜索边插点的方法，因此内插等高线点及其排列是同时完成的，其主要过程如下： (1)确定等高线高程。 为了在整个绘图范围中绘制出全部等高线，首先要根据DTM中的最低点高程与最高点高程计算最低等高线高程与最高等高线高程：当则。其中Z为等高距，INT为取整运算，即截去小数部分，则各等高线高程为(2)计算状态矩阵。为了记录等高线通过DTM格网的情况，可设置两个状态矩阵H（K）与V（K）序列 分别表示等高线穿过DTM网格水平边与竖直边的状态：m+1为DTM的行数，n+1为DTM的列数。其中hij与vij的值也可用逻辑值“真”(true)与“假”(false)。 由于格网(i，j)水平边有高程为zk的等高线通过的条件为：等高线高程介于DTM网格水平边两端点高程之间，即这两个条件等价于同理格网(i，j)竖直边有高程为zk的等高线通过的条件为:则状态矩阵H（K）与V（K的元素为为了避免上述判别式为零的情况可将所有等于等高线高程的格网点上的高程加(或减)上个微小的数o： 若Zi,jzk，则在实际进行软件设计时，H（K），V（K）(K1，2，)与H（0）V（0）占据相同的存贮单元。(3)搜索等高线的起点。 与边界相交的等高线为开曲线，而不与边界相交的等高线为闭曲线。通常首先跟踪开曲线，即沿DTM的四边搜索，所有的元素均对应着一条开曲线的一个起点(或终点)。在搜索到一个开曲线的起点后，要将其相应的状态矩阵元素置零。处理完开曲线后，再处理闭曲线。此时可按先列(行)后行(列)的顺序搜索DTM内部格网的水平边(或竖直边)，所遇到的第一个等高线通过的边即闭曲线的起点边。闭曲线起点对应的矩阵元素仍保留原值1。（4)内插等高线点。等高线点的坐标一般采用线性内插。网格(i，j)水平边上等高线点坐标(xPyP)为(1540)其中xix0+ix；yjy0+jy；(x0, y0)为DTM起点坐标；x，y为DTM方向与方向的格网间隔。格网(i，j)竖直边上等高线点的坐标(xq，yq)为(1541) (5)搜索下一个等高线点。 在找到等高线起点后，即可顺序跟踪搜索等高线点。为此可将每一DTM格网边进行编号为1，2，34(如图l53所示)，则等高线的进入边号IN有4种可能。设进入边号为1，即IN1，按固定的方向(顺时针或逆时针)搜索等高线穿过此格网的离去边号OUT；如按逆时针方向搜索，则首先判断2号边，其次3号边，最后4号边，即对于上述三种情况，下一格网的编号与进入边号分别为同理可分析处理进入边号IN23，4的情况。在搜索到下一个等高B线点后即按式(1540)或式(1541)计算该点坐标。将每一搜索到的等高线点对应的状态矩阵元素置零是必要的，它表明该等高线点已被处理过了。当状态矩阵H（K），V（K）变为零矩阵时，高程zk的等高线就全部被搜索出来了。以上是格网中仅有高程为zk的一条等高线通过时的情况。若格网中有高程为zk的两条等高线通过，即该格网中4条格网边都有等高线穿过，这种特殊情况可按上述统一的逆时针方向(或顺时针方向)搜索下一等高线点，也可以借助于格网中心点高程Zc：进行判断离去边。仍考虑该格网进入边号IN1的情况，此时离去边号只可能是OUT2或OUT4，而不可能为3号边，即OUT3。另一种处理办法是根据在格网中的地面可用双曲抛物面近似，用Zzk平面切割该双曲抛物面得到一双曲线，该双曲线即该格网内高程为zk的两条等高线。双曲面方程为设以DTM格网(i，j)左下角为原点，边长为单位长度1，由格网四角点坐标计算系数:则双曲线方程为(1545)再将坐标系平移到双曲线的中心(x0，y0)(1546)其中()为新坐标系坐标。将式(1546)代人(1545)得(15-47)令的系数为零得(15-48)将式(1548)代人式(15-47)整理得因此当0时，等高线通过I，象限(图154)；当0时，等高线通过，象限(图155)。地形特征线是表示地貌形态、特征的重要结构线。若在等高线绘制过程中不考虑地形特征线，就不能正确地表示地貌形态；降低精度，就不能完整地表达山脊山谷的走向及地貌的细部。因此，必须在DTM数据采集、建立及应用的整个过程中考虑地形特征线。 在搜索等高线时，为了考虑地形特征线必须注意以下几点： 若在某一条格网边上有地形特征线(如山脊线)穿过。如图l56所示，必须采用特征线与格网线之交点(如固a，b,c)与相应的格网点(如图中LM，N)内插等高线点，而不能直接用格网点内插等高线。例如由aN一线性内插点1bM线性内插点2bL一线性内插点4 此时就可能在同一条格网边上出现两个等高线点，例如图中格网边ML上就出现了等高线点2和4。这样就个能仅仅以一个逻辑值(“真”或“假”)来简单地判断该格网线上是否存在等高线，因为特征线已将该格网线分成两个线段。为此，在一个计算机字中不能简单地赋予一个逻辑值，而应将一个计算机字分成相应的“字段”使用，最简单的“字段”是“位”(即bit)。例如，当特征线将格网线分成两段时，则在计算机中取2bit，分别以高位的“1”或“0”表示格网线的上段(或水平格网线之左段)有、无等高线通过，低位表示下段(或右段)有、无等高线通过。在跟踪按索等高线时，当等高线穿过山脊线(或山谷线)，还必须在山脊线(山谷线)上补插等高线点，例如图中由特征点b，c内插等高线点3。由图156的例子可以看出，当考虑了特征线时，内插出等高线点1，2，3，4，5，从而保证了山脊线的走向，正确地表示了地貌。否则，不考虑特征线时，只能内插得到两个等高线点1，5，因此难以保证地形特征与精度；当等高线遇到断裂线或边界时则等高线必须“断”在断裂线或边界线上。 特征线穿过一个DTM格网边共有6种情况(如图157所示)，它将一个格网分成2个多边形(三角形、四边形或五边形)，下一个等高线点的搜索应在由进入边(或进入的半边)与特征线及原格网边(或半边)组成的多边形内进行。当等高线穿过特征线时，则应继续在该格网的另外一部分(也是一多边形) 搜索离去边，此时，该格网含有该等高线的三个点两个是与格网边的交点，一个是与特征线的交点。对于断裂线，则不存在离去边，等高线就终止在断裂线上对闭曲线应从其起点向另一方向搜索。(6)搜索等高线终点。 对于开曲线，当一个点是DTM边界上的点时，该点即为此等高线的终点。对于闭曲线当一个点也是该等高线第一点时，起点即为其终点。由于在搜索闭曲线起点时，保留其对应的状态矩阵元素为l，这就保证了能够搜索到闭曲线的终点。 2整体内插出整个数字地形模型内全部等高线穿越格网边的交点，然后按每条等高线将这些点分别排列并存贮 这种方法是按数字地形模型的格网边的顺序(例如先按行、后按列)内插计算出全部等高线穿越格网边交点的坐标x，y，然后按等高线的顺序将属于每一条等高线的点找出来，并按等高线的走向将它们顺序排列，并存贮在磁带或磁盘上。在格网边上内插等高线点的坐标可以按上述的线性内插的方法进行，在格网的密度较稀或精度要求较高时也可用高次多项式内插。将离散的等高线点顺序地排列起来，可按如下两个条件进行：方向条件。要求从已经排列好的两个相邻等高线点出发，至下一个等高线点的方向变化为最小；距离条件。要求从一个已经排列好的等高线点到下一个点之间的距离为最小。一般来说，在排列时应综合考虑此两个条件。因为若仅仅考虑一个条件有时候会导致等高线的错误排列。如图158(a)表示仅仅按方向条件排列，图(b)表示仅仅按距离条件排列而导致错误的例子。但在实际计算程序中，为了尽可能节省计算时间，主要是按距离条件排列。因为所谓“距离条件”通常可以化为“坐标增量条件”，也就是要求相邻两个等高线点(Xl，Y2)，(X2，Y2)之间的坐标增量无论是X= X2Xl或YY2Y1，均不得超过DTM的格网边长L，即 XL及YL同时要求相邻两个等高线点之连线不得穿越格网线。假定X2Xl，Y2Y1时，此条件可用下列不等式表示。是一个任意小值，如10-5(注意此时整个DTM是以格网起点为平面坐标的原点) 。如图158中所示的两个错误都是由于“相邻两个等高线点”之连线穿越格网线所产生的。二、等高线光滑（曲线内插)由上述步骤获得的是系列离散的等高线点，即等高线与DTM格网边及特征线的交点， 显然，若将这些离散点依次相连只能获得不光滑的由一系列折线组成的“等高线”。为了获得一条光滑的等高线，在这些离散的等高线点之间还必须插补(加密)。插补的方法很多，一般来说对于插补的方法有以下的要求：(1)曲线应通过已知的等高线点(常称为节点)；(2)曲线在节点处光滑，即其一阶导数(或二阶导数)是连续的(3)相邻两个节点间的曲线没有多余的摆动；(4)同一等高线自身不能相交。目的，常用的些插补方法都能严格满足上述的条件保证，特别是当节点分布不均匀或较稀疏时，问题更为突出。“张力样条”函数的插补方法，主要是针对解决曲线的多余摆动而提出来的。但是，由于数字地形模型的格网较密，离散等高线点分布比较均匀，而且比较密集，因此，一般来说利用分段三次多项式插补方法也能满足后两个条件。 经过上述的等高线跟踪与光滑处理，即可将等高线图经数控绘图仪绘出或显示在计算机屏幕上。153 立体透视图 从数字高程模型绘制透视立体图是DEM的一个极其重要的应用。透视立体图能更好地反映地形的立体形态，非常直观。与采用等高线表示地形形态相比有其自身独特的优点，更接近人们的直观视觉。特别是随着计算机图形处理工作的增强以及屏幕显示系统的发展，使立体图形的制作具有更大的灵活性，人们可以根据不同的需要，对于同一个地形形态作各种不同的立体显示。例如局部放大，改变Z的放大倍率以夸大立体形态；改变视点的位置以便从不同的角度进行观察，甚至可以使立体图形转动，使人们更好地研究地形的空间形态。 从一个空间三维的立体的数字高程模型到一个平面的二维透视图，其本质就是一个透视变换。我们可以将“视点”看做为“摄影中心”，因此我们可以直接应用共线方程从物点(X，Y，Z)计算“像点”坐标(x，y)，这对于摄影测量工作者来说是一个十分简单的问题。透视图中的另一个问题是“消除”的问题，即处理前景挡后景的问题。 从三维立体数字地面模型至二维平面透视图的变换方法很多，利用摄影原理的方法是较简单的一种基本分为以下几步进行： (1)选择适当的高程Z的放大倍数m与参考面高程Z。这对夸大地形之立体形态是十分必要的，令Zijm(ZijZ0)。(2)选择适当的视点位置XS，YS，ZS；视线方位t(视线方向)， (视线的俯视角度)。如图159所示，S为视点，SO(y)轴是中心视线(相当于摄影机主光轴)，为了在视点S与视线方向SO上获得透视图，先要将物方坐标系旋转至“像方”空间坐标系Sx1y1z1：或(1551) 式中 在通过平移旋转将物方坐标X，Y，Z换算到像方空间坐标x1，y1，z1以后，怎样通过“缩放”，投影到透视平面(相当于像面)上，即怎样设置透视平面到视点S的距离像面主距f比较合理的方法是通过被观察的物方数字高程模型的范围Xmax，Xmin，Y max，Ymin以及像面的大小(设像面宽度为w高度为H)，自动确定像面主距f，其算法如下： 计算DEM四个角点的视线投射角，之几何意义见图1510，i=1，2，3，4；xli，yli，z1i是由DEM四个角点坐标例如(Xmin，Ymin，Z1)通过公式(1551)所求得的四个角点的像方空间坐标。 从中选取max，min，max，min，即从而再由像面曲大小求主距：(3)根据选定的或计算所获得的参数XS，YS，ZS，以及主距f计算物方至像方之透视变换，得DEM各节点之“像点”坐标x，y：(4)隐藏线的处理。 在绘制立体图形时，如果前面的透视剖面线上各点的z坐标大于(或部分大于)后面某一条透视剖面线上各点的z坐标，则后面那条透视剖面线就会被隐藏或部分被隐藏，这样的隐藏线就应在透视图上消去，这就是绘制立体透视图的“消隐”处理，如图1511所示。 欲根本上解决这一问题是比较困难的，主要是计算量太大，一般经常使用的种近似方法被称为“峰值法”或“高度缓冲器算法”，名称虽各不相同，但其基本思想是相同的。 基本思想是将“像面”的宽度划分成m个单位宽度x0,例如对于一个分辨率为1024个像素的图形显示终端，则可以将整个幅面分成1024个像素，即单位宽度为像素，又如在图解绘图时，可令单位宽度x00.1mm(或0.2 mm)，则将绘图范围划分为m列，定义一个包含m个元素的缓冲区zbufm，使zbuf的每一元素对应一列。 在绘图的开始将缓冲区zbuf全部赋值zmin(或零)，即以后在绘制每一线段时，首先计算该线段上所有“点”的坐标。设线段的两个端点为Pi(Xi,Yi)与Pi+1(Xi+1,Yi+1)，则该线段上端点对应的绘图区列号即缓冲区zbuf的对应单元号为Pi与Pi+1之间各“点”对应的缓冲区单元号为它们的z坐标由线性内插计算为当绘每一“点”时，就将该“点”的z坐标z(k)与缓冲区中的相应单元存放的z坐标进行比较，当时，该“点”被前面已绘过的点所遮挡，是隐藏点，则不予绘出。否则，当时，该“点”是可视点，这时应将该“点”绘出，并将新的该绘图列的最大高度值赋予相应缓冲区单元在整个绘图过程中，缓冲区各单元始终保存相应绘图列的最大高度值。 (5)从离视点最近的DTM剖面开始，逐剖面地绘出，对第一条剖面的每一格网点，只需与它前面的一个格网点相连接；对以后的各剖面的每一格网点，不仅要与其同一剖面的前一格网点相连接还应与前一剖面的相邻格网点相连接(当然，被隐藏的部分是不绘出的)。 (6)调整各个参数值，就可从不同方位、不同距离绘制形态各不相同的透视图制作动画。省计算机速度充分高时就可实时地产生动画DTM透视图。154 DEM的其他应用一、坡度、坡向的计算 1斜平面与地平面之夹角 如图1512所示，一斜平面与水平面之夹角为， 且其在X方向与y方向之夹角为x与Y，显然 又因为所以2坡向QO与x铀之夹角T为坡向角3由4个格网点拟合一平面之坡度设平面方程为 以(0，0)点为原点，可列误差方程法方程为解为因为平面的法矢量为n=A B -1T，所以坡度角之余弦为Z方向单位矢量0 0 -1T 与n之数积所以坡向角之正切为所得公式与斜平面情况的公式相同。二、面积、体积的计算1剖面积根据工程设计的线路，可计算其与DEM各格网边交点Pi(Xi,Yi, Zi)，则线路剖面积为其中n为交点数；Di,i+1为Pi与Pi+1之距离同理可计算任意横断面及其面积。 2体积 DEM体积由四棱柱(无特征的格网)与三棱柱体积进行累加得到，四棱柱体上表面用双曲抛物面拟合，三棱柱体上表面用斜平面拟合，下表面均为水平面或参考平面，计算公式分别为其中S3与S4分别是三棱柱与四棱柱的底面积。 根据新老DEM可计算工程中的挖方、填方及土壤流失量。 3表面积对于含有特征的格网，将其分解成三角形，对