地表面积估算与资源分配的DEM模型.docx

地表面积估算与资源分配的DEM模型摘要本文针对三维地形的地表面积的大小以及监测器的合理布置以求达到最有 效的监测范围问题进行了模型分析，建立了 DEM模型和三维空间最佳面积覆盖模 型，考虑该矩形区域内地形表面的实际情况，综合影响模型准确程度的一系列因 素，对其进行模型精确度的分析问题一：根据已知测量点数据，利用二维差值法生成格网结构的数字高程模型， 运用DEM三维地形地表面积递推算法求解地表面积，其精度用残差总和的大小与 回归分析法进行综合评估问题二：对于矩形区域的无盲区覆盖问题，建立了三维空间的随机最佳覆盖模型， 运用分布式启发算法完成整个矩形区域内无盲区的覆盖问题文章最后，对该文所应用的模型在其他一些重要领域进行了简单的介绍，扩 展了模型从而加强了它的普遍适用性关键词：二维差值法DEM模型残差总和回归函数随机最佳模型1 问题的提出2500200016001000主12米以上部分的500问题一:试 地表面积要在长4000 米，宽 3000 米的某矩形区域内进行一项野外作业，需要布设一 些监控器附件 data5B.txt 中给出了该区域的步长为 5 米的网格节点对应的海拔高 度值及等高线图，其概貌如下图所示：_ ” /|丿、, 7 I yj j/ I tZLi VJ # I ^-rzJ N /乂 , y问题二：如果要在该矩形区域内布设某型号的监控器，已知在无障碍情况下 最远可观测距离为500米，为实现整个矩形区域的无盲区全覆盖监视，需研究在 矩形区域内设置多少个监视器？其最佳位置在何处？若在该区设置 28 个监控 器，这些监控器应该设置在什么位置？可以覆盖矩形区域内的多大面积。

2 问题的分析2.1 问题一分析：已知该矩形区域的等高线图及相应海拔高度值（见附录），要对海拔在12 米以上的部分 地表面积进行计算在计算实际地形的地表面积时，因其地形是随机起伏的，很难准确的运 用函数来描述这种起伏状况，而用传统的数学方法又不易得到准确值，故运用一种基于逼近 思想的递归算法来求解DEM模型，从而算出三维地形的地表面积对其模型精度的分析评 估中，影响其精度因素主要有地形表面的特征、原始数据的分布及密度、DEM表面的建模 方法以及其表面特性等，而原始数据的内插法为直接影响因素，因此运用残差总和的大小进 行比较，再者结合回归函数分析与之比较，可知其精度较高，模型较为准确2.2 问题二分析：要实现对整个矩形区域内无盲区全覆盖监测，可以通过设计一个基于位置的 算法来找到一条连接网络中任意 S 源节点和 Sink 指令中心的路径，并使该路径 上所有点的最小被观察能力最大化，转化为无线传感器网络中的“最佳覆盖问 题”，运用分布式启发算法进行求解3 模型假设1）问题的描述涉及到地表的规模，当规模发生变化后，问题性质不发生变化2）问题的出口条件定为判断各个小三角形面积是否小于阀值3）监测器在任何情况下均可以达到最大范围的监测效果。

4 符号系统x ：表示 x 轴坐标 ;y ：表示 y 轴坐标 ;z ：表示 z 轴坐标 ;s ：三角形 1 的面积 ;1s ：三角形 2 的面积 ;2s ：三角形面积 ;p : 第一点的坐标 ;1:三角形边长和的一半;p ：第二点的坐标 ;2p ：第三点的坐标 ;3p ：第四点的坐标 ;4p : 一点二点之间的距离 ;12p : 一点三点之间的距离 ;13p : 二点三点之间的距离 ;23p : 一点四点之间的距离 ;14p : 二点四点之间的距离 ;24p : 三点四点之间的距离 ;34p : i与j两点之间的距离ijS :残差总和ey :线性函数回归1y :二次函数回归2y 3 :负指数函数回归5 模型的建立与求解5.1 问题一 针对问题一，在计算地表面积时，首先圈定一个多边形，因三维地形随机起 伏，故此多边形内部的三维地形表面也是随机起伏的在已知平面内建立空间直 角坐标系O-xyz，将多边形的各个顶点通过鼠标在三维空间中进行三维拾取，然 后从目标多边形ABCDE（如图一）的某一个固定的顶点出发，依次连接各个顶点， 这样便可把多边形分解为一个一个的三角形ABC （如图二），为了能够更加精确 的计算出三角形固定的三维地形的表面积，再对三角形依照地形特征进行拆分， 关键是要求出三角形内的一个点，这个点能够代表出这个地形的地表特征，将其 与其余的点相连接，得到三个更小的三角形（如图三），以此类推，直到三角形 的面积小于一个阀值或者拆分位于三角形的平 面内， 最后利用海伦公式 S = Jp（p - a）（p - b）（p - c）， p = " + " + °，可得三角形面积。

将逼近思想运用^2于本题，分割成n*n个目标多边形，把每个多边形分成三个三角形，运用海伦公 式求此三角形的面积，利用循环解把n*n个目标多边形面积相加可得面积总和一〜 \ /图二 未分割的三角形示意图图三 分割后的三角形示意图图一目标多边形示意图求解如下：建立空间直角坐标系O-xyz，取四个峰顶为四个顶点，连接四个 点可得到一个四边形，任意连接一对角线，将四边形分割成两个三角形（如下图 四）利用海伦公式及逼近的思想，可近似的得到地表面积的大小若把本题表面分解成 3000*4000 个矩形，把每个矩形分成 2 个三角形，利用海伦公式和MATLAB得出循环解求得总面积为S = 9.7633*106 + 0.06 .因此地表面积为（9.7633*106 + 0.06）m2对运用 DEM 模型求得的三维地表面积大小的精确度进行残差总和和回归函 数模型进行综合评估，公式内容如下：残差总和公式S丄£ e 2，其中残差总和越小，模型的精确度越高e 飞 n - 2 ii=1线性函数回归y二P + 0 x+r101二次函数回归 y = 0 + 0 x+ 0 x 2 +r2 0 1 2负指数函数回归y =0 e-x +r30式中， r 为回归方程的函数，是判断回归函数拟合效果的重要指标，运用 MATLAB 软件分别求出线性、二次、负指数函数回归三个函数各自的残差总和， 由算法结果可知：负指数函数的残差总和最小，因此采用负指数函数回归作为对 DEM 模型精确度的分析。

经 MATLAB 软件求得，该模型的计算误差约为 1.02%，精确度相对较高5.2 问题二 针对问题二，因为所有节点的监测半径是相同的，当找到一个所有节点覆盖 半径的最小值 r 满足条件：只要经过监测器节点为圆心， r 为半径所确定的多个 圆内部就可以找到一条连接网络中 S 源节点与 Sink 指令中心的路径 求解过程如下：根据现有的 3 维 WSN 拓扑图生成化简但不失地理特性的 3 维德洛 内三角图 D （x ）（如下图五）留映射为对应的 2 维德洛内三角图 D （x ）（如3 3 2 2下图六），顶点集就可以尽可能多的保持原有 3 维 WSN 拓扑图的基本地理特征 在保持 3 维 WSN 拓扑图一定级别的地理特征条件下， 3 维 WSN 拓扑图的覆盖即可 转化为 D （x ）的覆盖问题，故只需保证其构成的三角形可被传感器节点覆 33盖监测即可图五 3 维德洛内三角图图六 2 维德洛内三角形示意图为了实现整个矩形区域的无盲区全覆盖监视，在该矩形区域内应设置监视器 的数量如下图所示，图中所有圆的个数即为所要求得的监测器的数量，共有 30 个监测器，最有可能的位置便是在图中所示圆的中心处设置。

图七 区域内监控器个数示意图如果在该区域内设置 28 个监测器，最有可能的设置位置便是在图六 2 维德洛 内三角形的各个顶点处，那样便可以得到三维地形的最佳覆盖面积6 模型的评价和推广6.1 优点：1 ） DEM 模型求解面积具有客观性、实际性，在现实生活中易被采纳，具有可 操作性实用价值2 ）递归算法能够灵活调整阀值，更好的控制模型的计算精度3 ）考虑较全面，运用较广泛，实际性强6.2 改进之处：1 ）选择恰当的高程点来确定拆分三角形较为困难2 ）模型是在理想的情况下建立的，并没有考虑野外的动态因素，因此模型仍 存在着不足，但总的来说，还是具有很好的借鉴意义和指导价值的6.3 模型的推广应用DEM 模型已经在各领域均有广泛的应用，比如线路勘测设计、水利建设工程、 军事工程等等方面，对于自然灾难而形成的洪水、泥石流淹没区的估算，山区日 照分析和 GIS 中的三维可视化，以及本文涉及的平面等高线图、立体等高线图， 只要是与地理位置相联系的空间定位，均可以用到数字高程模型7 模型改进本文所建模型是从地势表面宏观方向考虑的，从而忽略了较为微小的一些真 实地势概貌，若将插值单元的长度与宽度尽可能的缩小以至趋于无穷小，那么所 建模型变更能反映出山势地貌的真实情况，计算出的地表面积以及监测器的合理 布置的精确性也会有大大的提高。

参考文献1］李志林，朱庆.数字高程模型 Im ］ . 武汉：武汉测绘科技大学出版社， 2004 b］孙家广•计算机图形学（第三版） Im ］ . 北京：清华大学出版社， 2003b］任彦，张思东•线传感器网络三维空间最佳覆盖路由协议•电子学报，2006, 2B］朱衡君.MATLAB语言及实践教程•北京：清华大学出版社，北京交通出版社， 2005k］姜启源•数学模型•北京：高等教育出版社，1987k］杨晓云，顾利亚.字高程模型的误差处理与精度评估•铁道勘查，2004，5侯飞，韩奎峰.多分辨率DEM的插值残差计算及精度评价模型.工程勘察2010,附录一:（ 步长为 5m 的三维立体图形的程序）30002500200015001000500°0 500 1000 1500 2000 2500 3000 3500 4000图为该区域等高线示意图clc, clear a=load('data1.txt');amin=min(min(a)), amax=max(max(a)) x0=0:5:4000; y0=0:5:3000;[xx0,yy0]=meshgrid(x0,y0); mesh(xx0,yy0,a) v=[-49,1,12,44,75,107,171,184:80:357,357:100:614,614];figure, c=contour(x0,y0,a,v); clabel(c)附录二：(地行地表面积的算法程序) clc, clear %该程序计算最小面积a=load('data1.txt');a=a';x0=0:5:4000; y0=0:5:3000; pp=csape({x0,y0},a) x=0:4000; y=0:3000;z=fnval(pp,{x,y});[m,n]=size(z);s=0;for i=1:m-1for j=1:n-1p1=[x(i),y(j),z(i,j)]; p2=[x(i+1),y(j),z(i+1,j)]; p3=[x(i+1),y(j+1),z(i+1,j+1)];p4=[x(i),y(j+1),z(i,j+1)];p12=norm(p1-p2); p23=norm(p3-p2); p13=norm(p3-p1); p14=norm(p4-p1); p34=norm(p4-p3);if all([z(i,j),z(i+1,j),z(i+1,j+1)]>=12) z1=(p12+p23+p13)/2;s1=。