• 4.14 MB
  • 61页

基于dem内插的工程土方量计算方法研究

  • 61页
  • 当前文档由用户上传发布,收益归属用户
  1. 1、本文档共5页,可阅读全部内容。
  2. 2、本文档内容版权归属内容提供方,所产生的收益全部归内容提供方所有。如果您对本文有版权争议,可选择认领,认领后既往收益都归您。
  3. 3、本文档由用户上传,本站不保证质量和数量令人满意,可能有诸多瑕疵,付费之前,请仔细先通过免费阅读内容等途径辨别内容交易风险。如存在严重挂羊头卖狗肉之情形,可联系本站下载客服投诉处理。
  4. 文档侵权举报电话:19940600175。
'~00~~~=----~P2~5~8___Wt&:~jf-UDC:--------------*~~%:____,,,J"±".::f!.~II[~-5-:0212523$1w~&~~=i~xm~aAA:2o1s.6.1si"Bx~maAA:2015.6.9 *A~~AA~~~~~~~£*A~~~m~~illff~~~I~~~~~m~d~.-T~~*fJ1!J1JuJ;l.f;f,tl:toJ&iMzJ:H1-,itJtr:f/G"E!.13-;l"tft!!ASt~1t*~mEJufr.J1iJf~~*,ili::pf~!13"1H~~tH3tm:A:¥Meltitt!~JHJL:ftJa1:5000<1实地测量1:5000150m20cm航测成图大比例尺1:1万2.5100m10cm综合程度低,较真航测成图1:2.5万5250m4cm实反映地形地貌航测成图1:5万10500m2cm航测成图1:10万201km1cm编绘成图中比例尺1:25万502.5km4mm一定程度综合,近编绘成图1:50万1005km2mm似反映地形地貌编绘成图小比例尺1:100万10km1mm较高综合程度,反编绘成图映地形大致特征3.2.3实地现场数据采集运用先进的测绘仪器,如GPS、全站仪、三维激光扫描等。在已经控制点的测站上,广策目标点方向、距离和高差、通过计算机软件内业处理即可解算出目标点的x、y、z坐标,适当变换之后可直接生成DEM。实地现场测量可以获取高精度的高程数据,一般用于小范围区域内大比例尺测图同时对高程精度要求较高的工程项目。例如房屋建筑、场地平整、公路铁路前期勘察设计、矿山探测工程、水利枢纽工程等。然而,由于这种方法工作任务繁重、数据采集周期长、数据更新时间跨度长、人力财力耗费较高,通常在大规模的采集中不使用这种方法。总之,这种通过地面实地现场测量自然地表来获取数据的方法,具有数据精度高、适应范围小、费用高等特点。3.2.4既有DEM数据基于测绘与地理信息产业的需要,各个国家基础地理信息数据有了长足的发展和进步,因此各种数字地图产品陆续问世。目前,各个国家都建立了覆盖本国国土并适合本国国情的全系列比例尺DEM。为了满足国家决策和研究的需要,我国也已经逐渐建立了18 兰州交通大学硕士学位论文覆盖整个国土范围的1:5万、1:25万、1:100万数字高程模型,同时为了提高应对复杂自然灾害的应急处理能力,我国还建立了七大江河重点防洪区的1:1万数字高程模型。极大地保障了人民生命和财产安全。各个省也已经分批次开始建立覆盖本省省域的1:1万数字高程模型。在我国,各系列比例尺、各种分辨率DEM产品的建库和生产均由国家基础地理信息中心负责。当运用既有DEM来生产新DEM时,应充分顾及研究目的以及现有DEM的分辨率、数据存储格式、现有数据的精度等因素。3.2.5DEM数据采集方法比较DEM的数据源质量与精度和选择最优的数据获取方式直接影响DEM的最终质量,数字摄影测量是进行数据库更新效率最高的方式之一;以等高线地形图为数据源来生产DEM的方法可以大量推广并应用于实际生产;GPS、激光扫描、雷达干涉在DEM数据采集方面很有发展前景。我们可以从精度、速度、成本、更新程度等方面来全面评价各种采集方法,在实际采集中,也要全面考虑研究目的与需求、精度要求、经济承受能力、硬件软件条件。权衡各种因素后,选择最优最适宜的采集方法。因此对DEM数据采集方法与各自特性进行了比较(表3.2)。表3.2DEM采集方法与特性比较获取方式DEM精度速度成本更新程度应用范围地面测量非常高(cm)耗时很高很困难小范围区域,特别的工程项目摄影测量比较高(cm到m)比较快比较高周期性大工程项目,国家范围数据收集立体遥感低很快低很容易国家范围乃至全球范围数据收集GPS比较高(cm到m)很快比较高容易小范围,特别的项目地形图数字比较低(图上比较耗时低周期性国家范围,军事采集,化0.2到0.4mm)中小比例尺地形图数据获取激光扫描、非常高(cm)很快非常高容易高分辨率,各种范围干涉雷达-19- 基于DEM内插的工程土方量计算方法研究3.3数字高程模型建模的基本方法DEM是地形表面的一个数学模型,可以用一个或多个数学函数来表示DEM,对地形表面进行表达的各种处理可称为表面重建或表面建模。DEM表面用数学表达式可描述为:YX),f(Z(3.4)针对此表达式,总结了最常用的几种多项式函数(表3.3)。表3.3用于表面重建的通用多项式独立项项次表面性质项数aZ0次项平面101次项线性2aXaY12222次项二次抛物面3aXYYaXa34533223次项三次曲面4aYXYaXaYXa67894432234次项四次曲面5aYXYXaYXaYaXa101112131455次项五次曲面6Xa15通常情况下,地形表面建模主要为基于点的建模、基于三角形的建模和基于格网的建模三种,下面将分别对它们进行详细介绍。3.3.1基于点的表面建模在建模方法中基于点的表面建模所使用的逼近函数为表3.3的第一项,即),(ay0xfZ(3.5)其表示一空间小区域水平面,在这里仅需要一个采样点就能确定一个空间小区域水平面,待定系数a0就是该采样点的高程z值,即0ziaz(3.6)整个水平面的高度也就为a0。将建模区域按照采样点来划分子块是因为每个采样点都确定了一个空间水平面,每个字块包含一个采样点,是围绕采样点的一个邻近区域。各个子块拼在一起即可构成一个不连续的数字高程模型。虽然这种基于点的建模操作起来较为简单,可是其难点却在于无法确定各个子块边界和形状。如果单从理论上来分析,这种方法只针对某些独立的点,所以在处理数据时,20 兰州交通大学硕士学位论文对于数据类型是不作限制的。如果采样点是不规则分布的,很难确定各个子块的形状和边界,一般可以采用采样点的Thiessen多边形作为子块的形状和边界。如果采样点是规则分布的,就很容易确定每个子块的形状和边界,而且这些子块可以自然拼接起来。所以基于点的建模仅适合于规则分布的采样点。而这里所建立的模型不是一个连续的模型,因此基于点的建模方法在实际中并不是一种很好的方法。3.3.2基于三角形的表面建模在建模方法中基于三角形的表面建模所使用的逼近函数为表3.3的前三项,即Zf(x,y)a0a1xa2y(3.7)此函数为一个线性函数,他表达的是一个空间水平面或倾斜面,在数学上,确定空间一个平面只需要空间内不在同一条直线的三个点即可。所以我们只要三个采样点,将这三个点的平面坐标和高程代入上式,得到关于待定系数的三个线性方程,解此方程组,可得到三个待定系数a0,a1,a2从而唯一的确定一个空间平面。基于三角形的建模是由三个采样点确定空间中一个平面,而如果该空间平面的覆盖范围只为此三个采样点围成的空间三角平面,那么该建模方法就是将建模区域划分为一系列三角形,每个三角形只由这个三角形的三个顶点即采样点唯一确定,每一个三角形连续拼接起来就构成了一个连续的数字高程模型,该模型就称为不规则三角网模型也可以称为TIN(TriangulatedIrregularNetwork)模型。在几何学中,我们都知道三角形在通常情况下是平面图形中最基本的单元,不论是正方形、菱形这样的简单图形还是其他别的形状的多边形归根结底都是由多个三角形以某种形式组合而成的。无论是用哪种采样方式生成,还是由等高线法生成建模数据。用三角形作为最小单元的表面建模在任何数据结构中都能有非常好的表现。由于三角形不论在整个形状还是在边的大小方面变动都有很大的随意性,所以这种建模方法可以轻松地与陡坎、斜坡、生成线或其他任何较破碎地形数据进行融合,在计算每个空间三角平面时也不复杂,但难点之处在于如何选择三个采样点来确定一个空间平面,同时要保证各个平面之间无缝隙、不重叠地拼接起来,形成一个连续的三角网模型,这是建立三角网模型的难点和关键,也是建立TIN模型的主要工作。3.3.3基于格网的表面建模基于格网的表面建模所使用的逼近函数为表3.3的前四项,即Zf(x,y)a0a1xa2ya3xy(3.8)此函数并非为一线性函数。当x坐标固定时,函数为y的线性函数,当y坐标保持不变时,函数为x的线性函数,所以又称为双线性函数,此函数确定的空间表面也就称-21- 基于DEM内插的工程土方量计算方法研究为双线性曲面。这里需要四个采样点数据即可确定一个双线性曲面,将这四个点的平面坐标和高程代入上式,得到关于待定系数的四个方程,解此方程组,可得到四个待定系数a,a,a,a从而确定一个空间双线性曲面。0123基于格网的建模是由四个采样点确定一个空间双线性曲面,当此空间曲面的覆盖范围为四个采样点所围成的空间四边形,则该建模方法是将建模区域划分为一系列四边形子块,各个子块之间能够无缝隙、不重叠地连续拼接起来,就构成了一个连续的数字高程模型,该模型称为格网模型。从理论的角度来看,此方法适用于任何形式分布的采样点,并且每个空间双线性曲面的计算也不复杂,其难点同三角形方法一样也在于选择哪四个点来确定一个空间双线性曲面,并且各个曲面之间需要无缝隙、不重叠地连接起来形成一个连续的格网模型。对于不规则分布的采样点而言,由采样点直接连接一个连续的四边形网格十分困难,对于规则分布的采样点而言,一个连续的平面四边形网格已经自然形成,也就是说格网模型已经自动建立起来了,这里就称为规则格网模型或GRID,因此在实际应用中,格网模型就是这里所说的规则格网模型,即GRID模型。同时,还可以选择更加多项的高次的多项式来作为逼近函数去拟合实际地表,但考虑到对场地面积较大区域使用高次多项式函数在某些情况下会有可能导致构建的DEM并不符合实际,同时也考虑到实际需要,更高次多项式的逼近函数并不常用。从实用的角度来看,在数据处理过程中格网数据较其他类型数据更加有优势,运用规则格网采样方法和渐进采样方法获取的数据比较适合基于格网的表面建模,而正方形格网数据则更加适合基于格网的表面建模,所以很多软件的DEM模块只接受格网数据。而为了使用这种规则格网数据,首先必须对原始数据进行数据内插,使其成为格网数据来保证输入数据满足软件对数据格式的要求,这种方法用于处理较大范围平缓地区的整体数据,但遇到陡坎和斜坡以及大量断裂线等地表状况极其复杂的情况,如果要想使用这种方法就必须进行增加特征点数量或采集更多外业测量点的处理。3.4不规则三角网的建立3.4.1TIN的简介一系列不规则分布的数据点生成连续的三角面,这一系列三角面共同构成不规则三角网。在数字地形建模中,正是用这个不规则三角网来逼近地形表面。通常用节点(node)、边(edge)、和面(face)作为研究TIN的基本单元。节点:是相邻两个三角形的公共顶点,它同时也作为原始采样数据,用来构建不规则三角网。边:是两节点间的连线,也是两三角形的公共边,边直接表明不规则三角网是不光22 兰州交通大学硕士学位论文滑的。另外,断裂线、特征线以及区域边界也作为边的一部分存在于TIN中。面:是TIN描述实际地形表面形状的最小单元,是由三个节点按照相互距离最近的原则共同组成的一个最优情况下为等边三角形的面。TIN中的每一个三角形都能描述局部地形是水平还是倾斜状态,任意两个三角形都不能交叉也不能重叠。节点、边、面之间还互相存在着某些空间拓扑关系,例如关联、邻接等。3.4.2TIN的体系结构TIN在数据结构上是一种典型的矢量数据结构。这种结构主要通过节点、三角形边和三角形面来表达地形上各个离散点之间的空间拓扑关系,因此,不论对TIN的应用来说还是对数据库的维护来说,如何设计一个兼具结构紧凑、高效、维护方便等优点的TIN的存储与组织结构就变得极其重要。由于地形各个点之间存在着相关关系,相互接近的地形采样点如果它们之间距离越近则它们彼此间的某些性质就越相近。相对规则的三角形对其进行插值时结果要比狭长的三角形对其进行插值的结果更加让人信服。因此,在TIN中对三角形的几何形状有尽可能接近等边三角形、保证形成三角形的三个顶点之间距离最近、三角形网络必须唯一这三条原则。3.4.3TIN的三角剖分算法分类与特点TIN的三角剖分准则是指TIN中三角形生成所必需满足的条件,主要有空外接圆准则、最大最小角准则、最短距离和准则、张角最大准则、面积比准则、对角线准则。这些准则关系到三角形的面积大小和几何形状,更重要的是它们还直接关系到TIN的质量好坏。从DEM数据来源的角度来讲,采样数据一般呈不规则分布或者规则分布,在某些特殊情况下也可能沿等高线进行分布。根据这种分布情况可对三角剖分方法进行分类(图3.1)。Delaunay三角剖分是同时满足空外接圆准则和最大最小角准则的三角剖分,简称为DT。DT三角剖分作为目前应用最频繁和最普遍的三角剖分方法,其优势在于最大限度使三角形接近等边三角形并且从任何采样点构网都能保证最终形成的三角网是相同的。1977年Lawson提出了LOP法则。其基本思想为在由两个有公共边的三角形组成的四边形中,如果一个点落在了另外三个点构成的三角形中,则交换四边形的对角线。只要用LOP法则对在任何三角剖分准则下得到的TIN进行优化处理,就可以得到唯一的DT三角网。3.4.4Delaunay三角网的生成Delaunay三角网是实际中应用最广泛的TIN模型,它能充分反映点的密集程度和地形的复杂状况,最大限度表征地形的结构特征,在地形拟合方面其优势是显而易见的。-23- 基于DEM内插的工程土方量计算方法研究其生成的方法有很多,一般较常见且较为经典的有分割合并算法、三角网生长算法和逐点插入算法。间接DT分割合并算法DT三角剖分空外接圆算法不规则分布数据直接DT逐点插入算法三角形增长算法辐射扫描算法TIN算法类型退火模拟算法数学形态算法IPs算法规则分布数据循环迭代算法层次三角形算法等高线分布数据特征线算法探测优化算法图3.1TIN三角化算法分类(1)分割合并算法分割合并顾名思义就是将复杂问题简单化,整体问题分散化。将采样点集分割成若干个子集,对这些子集分别进行三角剖分,每个子集包含若干个少量的采样点,然后对每个子集进行Delaunay三角剖分,并用LOP算法将每个三角剖分优化为DT三角网,当每个子集各自剖分并且优化完毕后,将各个子集的Delaunay三角网合并,最终形成一个整体Delaunay三角网。具体过程如下:第一步,设采样点集为VP1(x1,y1),P2(x2,y2),,Pn(xn,yn)(3.9)将采样点集按升序进行排序。第二步,如果采样点集中的数据总数超过之前限定的阈值,则用垂直于横轴的分割线将采样点分割成数据个数基本相同的若干个子集,每个规范化的子集包含3~6个采样点。第三步,对每一子集的凸壳分别进行计算并将凸壳作为数据边界,采用任何一种方法对数据子集逐个进行三角剖分并通过LOP算法对其优化,将每一子集优化成为DT三24 兰州交通大学硕士学位论文角剖分。第四步,完成上述步骤之后接着找出连接左右子集的上部和下部公共切线。从下部公共切线开始从左到右为基线边,在基线边的左边找出左右两凸外壳点中夹角最大的点,连接构成新三角形新生成的连接左右两子三角网的边成为新的基线边,新的基线边仍然从左到右,直到基线边推进到上部公共切线为止。在第一步中的工作主要是为了使子三角网不相互重叠和交叉。在第二步中采样点分割是以递归的方式进行的,一般要求每一子集能采用同样的算法进行处理。在第三步中凸壳是采样点的自然极限边界,在所有包含全部数据点的凸多边形中,凸壳是面积最小的,在采样点集中连接任意两点的线段既不能在该凸多边形之外也不能与该凸多边形相交,只能完全位于该凸多边形中。在第四步中合并的方式是一种递归的方式进行的,具体步骤为同层优先,从下至上。这也是分割合并算法中“合并”一词的由来。(2)三角网生长算法三角网生长算法从生长过程角度可以分为收缩生长算法和扩张生长算法。数据点的分布数量和密度对收缩生长算法的运行有很大的影响,而在算法的实际运行过程中还会出现很多不可预知的状况,当边界收缩完成之后,之前一个完整的区域可能会被分解成若干个相互独立的子区域,这使得之后的剖分工作变得极其困难。所以收缩生长算法有很大不确定性,这里主要着重介绍扩张生长算法,此算法是从一个大致处在中心的三个点形成一个三角形,从这个三角形开始向外不断扩展,直到能够形成一个整体三角网最终包含到所有数据点为止。具体过程如下:第一步,生成初始内核三角形。在包含所有数据点区域的几何中心附近选取一点A,寻找距离此点最近的点B,连接后作为初始基线,利用空外接圆准则或者张角最大准则进行三角剖分,在附近区域中得到第三点C,这样,初始Delaunay三角形就形成了(图3.2)。第二步,扩展初始三角形形成三角网。接着以初始三角形的三条边分别为起始边线,利用前面提到的TIN的三角剖分准则中的空外接圆准则或者张角最大准则进行下一个符合条件点的搜寻,于是我们找到了D、E、F三个复合要求的点。在这一步需要注意的是初始边界A、B的连线将整个数据区域分成两个部分,必须要在初始三角形另一顶点的异侧区域内对下一个点进行搜寻,而不能在同侧进行搜寻。也就是说上图中AB为三角形ABC的初始边界,在搜寻下一个点时就要考虑这个点要在以AB为中轴且在C的异侧。同时为了避免三角形的交错和重叠,还要对上一步所选择的第三点进行进一步的确认,保证其一定是我们要的那个点。在这里,我们需要确认的是新三角形的三条边不能被前面所形成的三角形共用过两次或者更多,确认所依据的条-25- 基于DEM内插的工程土方量计算方法研究件是三角形任一条边最多只能与两个三角形所共用。根据这个条件来进行三角形全等比较,如果新三角形的三角形的确被之前的三角形共用过两次,则新三角形无效,如果没有被之前的三角形共用过两次,则三角形有效。图3.2初始三角形第三步,继续不断重复循环第二步的工作,直到初始多边形的每一条边的左边没有点为止。上述三角网生长算法向外扩展的是多边形,而不是三角形,其最终得到的是一个多边形,此多边形的每条边的左边没有采样点,因此该多边形为凸多边形,并且包含了所有的采样点。我们将包含所有采样点的最小凸多边形称为采样点的凸外壳,它既表示采样点的分布区域,同时也表示数字高程模型的覆盖范围。(3)逐点插入算法前面介绍的分割合并算法和三角网生成算法运用静态方法来完成全部构网工作,在整个三角网形成过程中,一旦三角网已经构建完毕,那么此时若再插入新点,这个新点不会破坏原有的已经连接好的三角网结构和连接关系。而这里要介绍的逐点插入算法则恰恰相反,它的构网过程是一个动态过程,每一个新加入的点都会对某些已经连接好的三角形产生影响,迫使他们改变连接方式。具体过程如下:第一步,生成采样点的凸外壳。首先需要确定采样区域的四个角点,得到初始的凸多边形a,b,c,d,其方向为顺时针方向,令采样点坐标为xi,yi,则有左下角点aminxiyi左上角点bminxiyi右上角点cmaxxiyi右下角点dmaxxiyi26 兰州交通大学硕士学位论文依次沿初始凸外壳每一条边的左边,找夹角最小的点,若找到这样的点,则将新点插入初始凸多边形中,同时修改多边形边界,从而形成新的初始凸多边形,重复此步骤,直到初始凸多边形左边没有点为止。此时的凸多边形就为采样点的凸外壳(图3.3)。图3.3生成采样点的凸外壳第二步,凸外壳的三角剖分。由凸外壳三角剖分得到的三角网是一个初始的三角网,在之后的新采样点加入时,还需要不断地进行三角网的调整,所以在此处不必考虑三角形三个采样点的邻近性,也就是将得到的凸外壳多边形划分成一些大的三角形,即仅仅用凸外壳边界点连接三角网(图3.4)。图3.4凸外壳的三角剖分第三步,依次插入凸外壳的每一个采样点。在已经连接好的三角网中,如果某个三角形的外接圆包含了该点,则标记该三角形,找出所有符合条件的三角形。这些三角形都将受到该点的影响,称为该点的影响三角形。所有由影响三角形构成的区域被称为该点的影响区域。删除该区域内所有相邻三角形的公共边,此时原有的三角网就变成一个多边形,称为该点的影响多边形。将该点与新生成的多边形的所有顶点连接,形成了一个内部三角网,然后对此三角网进行LOP处理,重复此步骤,直到所由采样点插入完毕。-27- 基于DEM内插的工程土方量计算方法研究三角网的逐点插入算法对于TIN模型来说是非常适合的,它可以对TIN模型进行增加采样点、删除采样点、移动采样点等的修改操作。3.5数字高程模型内插概述3.5.1DEM内插的概念数字高程模型内插在数学上属于数值分析中的插值问题。内插是二元函数的逼近,这是从数学角度来对其进行的定义。利用地面中已知采样点的三维坐标,形成一数学曲面且保证它是连续的,这些点在内插中不论平面坐标还是高程都是已知的,因此称为内插的参考点。将任一一个需要内插的点的平面坐标代入曲面方程,可以解出该点的高程数据。这就是曲面内插的问题,即数字高程模型内插问题。数字高程模型的内插实质是[26]建立一个与地面对应且双向连续的人造面,它的误差总在一个限度之内。数字高程模型的内插是模型建立、表达、分析、应用中经常遇到的问题,也是核心问题,数字高程模型内插始终伴随于DEM生产、质量控制、精度评定、分析应用的各个步骤。数字高程模型的每一种内插方法都是基于地形表面的空间自相关特性。仅凭上述定义无法真正揭示DEM内插的内涵,事实上,DEM内插是一个复杂的概念,[27]由于地形表通常复杂而多变,因此,某一种内插方法只能适用于特定的地形。一般我们认为,插值的准确性随着采样点数量的增多而提高,其模拟的表面就更加接近实际地[28]形。3.5.2DEM内插的基本原理在待内插点的某个邻近区域内,选择一部分参考点作为内插计算的已知点,再由这些参考点重建地形曲面,也就是建立曲面函数,这个函数称为插值函数或者插值曲面,最后将待内插点的平面坐标代入插值函数求得每个点相对应的高程。其关键环节是如何选择待内插点的邻近区域以及选择符合要求的参考点;如何建立一张曲面使其通过参考点或逼近参考点。目前生成DEM的方法主要是通过区区域内现有的高程点和等高线等信[29]息进行空间内插,以此来得到全区域的DEM。DEM内插的基本过程为:(1)确定待内插点邻近区域,选择满足阈值要求限定的参考点;(2)选择插值函数;(3)求解插值函数:(4)将待插点平面坐标代入插值函数求得待内插点高程。28 兰州交通大学硕士学位论文3.6DEM内插方法分类DEM从其概念提出到现在,经过多年来无数学者的发展和扩充,已经陆续研究出了很多高程内插方法,其分类并没有统一的标准。按参考点与插值曲面的关系可分为:(1)纯曲面内插所有的参考点都位于插值曲面上,即所建立的插值曲面通过所有参考点,这种内插称为纯曲面内插。(2)曲面拟合所建立的插值曲面不通过参考点,即参考点不在插值曲面上,通常是插值曲面逼近参考点,这种内插称为曲面拟合。常用的DEM内插方法如下(图3.5):整体内插纯曲面内插曲面拟合线性内插双线性内插DEM内插分块内插样条函数内插多面函数法内插移动曲面法逐点内插加权平均法图3.5常用内插方法按参考点的分布或者选择范围可分为:(1)整体内插整体内插是将模型区域作为一个整体来考虑,在整个模型区域范围内选择所有的点作为参考点来建立一个统一的插值函数。在这个区域内的任一待内插点都用这同一个插值函数来计算该点的高程。这种方法类似于整体建模方法。(2)分块内插分块内插就是先将模型区域划分成若干个子块,在某个子块内的待内插点选择该子块内的所有采样点作为参考点,用这些参考点建立一个只适用于该子块的插值函数来计算该点高程。当然,对于子块内的任何待内插点,均可以用该子块的插值函数来完成内插计算。(3)逐点内插逐点内插是以待内插点为中心,选取一定阈值为半径形成一圆形面状区域作为待内插点的邻近区域,此时选取该邻近区域内的所有采样数据点作为参考点建立插值函数来计算该点的高程。逐点内插方法因待插点的位置不同而使选择参考点范围也不同,因此-29- 基于DEM内插的工程土方量计算方法研究对于每个待内插点都要分别建立相应的插值函数,当待插点的位置移动时,包含采样数据点的范围也要跟着移动,所以,逐点内插方法也称为移动曲面内插法。3.6.1整体内插整体内插选择整个建模区域的所有采样点作为参考点来进行内插计算,基于整个建模区域来建立插值函数,用此插值函数来对此区域内任何一个待内插点来进行插值计算。插值函数选择一个二元高次多项式,即mmijzf(x,y)cijxy(3.10)i0j0式中,cij(i0,1,2,,m,j0,1,2,,m)为多项式的待定系数,m为多项式的次数。对于次数m的选择,为了确保能解算出待定系数,需要求采样点的个数必须大于或等于多项式的待定系数的数目。为了确定插值函数的具体形式就必须要对多项式的待定系数进行解算。设有n个采样点P1(x1,y1,z1),P2(x2,y2,z2),,Pn(xn,yn,zn)2多项式的待定系数c共有(m1)个,多项式的次数m需满足ij2n(m1)(3.11)才能求解。将n个参考点平面坐标和高程代入(4.2)式,得到mmijzkcijxkyk,(k1,2,,n)(3.12)i0j0这是关于待定系数的n个线性方程,可用矩阵表示如下AXZ(3.13)其中c00mmz11x1y1x1y1x1y1c10z2mm1x2y2x2y2xyc01A22XZz3(3.14)c11mm1xnynxnynxnynzncmm当参考点数n等于待定系数数目时,线性方程组有唯一解,由此可以得到待定系数为30 兰州交通大学硕士学位论文1XAZ(3.15)此多项式曲面通过所有采样点,属于纯曲面内插。当参考点数n大于待定系数数目时,线性方程组不能满足所有参考点,成为一组矛盾方程。可采用最小二乘原理来求解此方程,为此,给每个参考点高程施加一个改正数TV(v1,v2,,vn),得到误差方程VAXZ(3.16)T按最小二乘原理VPVmin解法方程为TTAPAXAPZ(3.17)其中,P为参考点的权阵,解方程得到T1TX(APA)APZ(3.18)此多项式曲面并不完全通过所有采样点,属于曲面拟合内插。整体内插是在整个建模区域上建立一张唯一的、光滑的曲面,能够充分反映宏观地形变化特征。但是如果要对采样点进行增减或移动,那么就要对插值函数多项式的系数作繁杂的改变,这会导致采样点之间出现某种不可预知的振荡现象,且这种振荡无法控制,这也将会直接引起函数不稳定现象的发生。同时,多项式的次数比较高时,对其进行高阶线性方程组求解时出现的计算误差和采样数据误差都会导致待定系数发生较大变化,从而导致高次多项式得到的解不稳定。这些缺点使得整体内插法在实际工作中应用较少。其主要用途是在宏观上对大范围地形的变化趋势进行模拟,同时在地形采样数据的粗差检验方面也有广泛应用,用整体内插来剔除某些在总体趋势中与宏观地形特征不相符的部分。3.6.2分块线性内插由于实际地形通常并不简单,会有许多复杂情况存在,因此我们无法用一个多项式来近似表达整个宏观地形和地貌,基于这种情况,我们在DEM内插过程中通常不采用整体内插法而采用分块内插法,分块内插就是将模型区域细化为若干个子块,各个子块使用不同的内插函数同时各个子块之间还能保持连续或光滑,如果插值函数为线性函数则称为线性插值。这里主要介绍双线性内插。双线性内插的插值函数为zf(x,y)a0a1xa2ya3xy(3.19)为求解此方程,需要有四个参考点,在实际应用中,通常是利用最靠近待内插点P(xp,yp,zp)的四个参考点P1(x1,y1,z1),P2(x2,y2,z2),P3(x3,y3,z3),P4(x4,y4,z4)来求得-31- 基于DEM内插的工程土方量计算方法研究插值函数的四个待定系数,将四个参考点的平面坐标和高程带入上式,则有za0a1xa2ya3xy,i1,2,3,4(3.20)iiiii解线性方程组得1a01x1y1x1y1z1a11x2y2x2y2z2(3.21)a21x3y3x3y3z3a31x4y4x4y4z4由此可知,双线性内插特别适合GRID模型的内插,对于待内插点P,只要找到P点所在方格,利用方格的四个格网点作为参考点确定一个双线性曲面方程,代入待内插点坐标,即可计算出待内插点的高程。对于GRID模型而言,如果将平面坐标系的原点平移到方格的左下角,则四个参考点的坐标可以表示为P1(0,0,z1),P2(d,0,z2),P3(d,d,z3),P4(0,d,z4),其中d为格网间距。相应的待内插点坐标变换为xxpx1yypy1(3.22)将四个参考点的平面坐标和高程全部带入式(3.20)得z1a0z2a0a1d2z3a0a1da2da3dz4a0a2d解出待定系数a0z1a1(z2z1)/da2(z4z1)/d2a3(z1z2z3z4)/d代入式(3.19)得双线性内插的直接计算公式xyyxxyxyzpz1(1)(1)z2(1)()z3()()z4(1)()(3.23)dddddddd如果将式(3.22)修改为x(xpx1)/dy(ypy1)/d(3.24)则基于GRID模型的双线性内插的直接计算公式简化为zpz1(1x)(1y)z2(1y)xz3xyz4(1x)y(3.25)32 兰州交通大学硕士学位论文3.6.3分块曲面内插分块曲面内插中主要介绍样条曲面内插,这里所要介绍的样条曲面如果用具体的形象解释就是将一张整体弹性较好的橡胶薄板压定在各个采样点所在的位置,而其他地方任由薄板自然弯曲。从数学的角度讲,样条曲面就是一个分段的不超过三阶的多项式。我们可以利用这个样条函数来获取在通过各个采样点的一个拟合曲面,并且保证此拟合曲面在相应点上的曲率是最小的。为了保证各个分块之间能够以一个相对平滑的方式过渡,样条函数规定公共边界上的导数连续,以此作为分块之间的连续性条件,这个条件是根据弹性力学条件而设定的。为各分块配置的样条函数一般采用三次多项式的形式,即33ijzf(x,y)aijxy(3.26)i0j0这里有16个待定系数,必须要列出16个线性方程来求解样条函数的具体形式。一般适用于GRID模型的内插,在每个方格建立一个样条函数,用于完成方格内点的内插。为求得方格(i,j)的样条函数并且使方格的样条曲面通过该格网的四个参考点,则样条函数应满足zi,jf(xi,j,yi,j)zi,j1f(xi,j1,yi,j1)(3.27)zi1,jf(xi1,j,yi1,j)zi1,j1f(xi1,j1,yi1,j1)由此可建立四个条件方程。为保持各方格样条曲面的光滑连接,要求在曲面的连接处或在参考点上具有连续的一阶偏导数和混合偏导数,由此可在每个格网点上列出3个条件方程,例如在第i行第j列的格网点上,可列出如下条件33zf(xi,j,yi,j)i1j()i,jaijixyijxxiji1j033zf(xi,j,yi,j)ij1()i,jaijjxyij(3.28)yyiji0j12233zf(xi,j,yi,j)i1j1()i,jaijijxyijxyxyiji1j1各个偏导数可用差商来估计,即-33- 基于DEM内插的工程土方量计算方法研究zzi,j1zi,j1()i,jx2dzzi1,jzi1,j()(3.29)i,jy2dzz()()2i,j1i,j1(zz)(zz)zyyi1,j1i1,j1i1,j1i1,j1()i,j2xy2d4d四个格网点上共列出了12个条件方程,加上之前的四个方程,总共16个方程,由此可解出16个待定系数从而求得此样条曲面函数的具体形式。与整体内插函数相比,样条函数保留了局部地形的细部特征,这在某些情况下是有很重要的意义的。同时样条函数使得构建的DEM既具有连续性又具有光滑性,这对DEM来讲是至关重要的,也是评价DEM质量优劣很重要的一个方面。样条函数在拟合时,由于样条函数多项式阶数较低,不会对数据误差的响应太过于敏感。综上所述,样条函数具有较好的保凸性、逼真性和平滑性。另一方面,样条函数在插值计算时只需要少量的参考点,因此计算速度比整体内插要快得多。3.6.4逐点内插分块内插在内插过程中其分块方式就已经可以确定,而一旦分块方式确定,那么其各个子块的大小、形状和在整个区域中所处的位置都将保持不变,同时包含的参考点的数量也将确定,由此确定子块的插值函数,待内插点落在哪一个子快上,就使用该子块对应的插值函数进行内插计算。逐点内插与分块内插不同,逐点内插是以待内插点为中心,定义一个待内插点的邻近区域,在该邻近区域内选择参考点,由此来确定邻近区域的插值函数。确定的邻近区域的大小、形状、和位置都随待内插点位置的变化而变化,由此确定的插值函数也将随着待内插点的移动而移动,因此逐点内插法又称为移动曲面法。移动曲面法对于每个待内插点可在其邻近区域内选择若干个参考点来分别建立一个插值函数,移动曲面法的插值函数形式如下22zf(x,y)a0a1xa2ya3xya4xa5y(3.30)这是一个二此多项式,有六个待定系数,所以至少需要六个参考点,当参考点多于六个点时,可按最小二乘法求解待定系数。参考点的选择一般要考虑在待内插点的多大范围内选择和在这一范围内选择多少个参考点这两个因素。参考点选择原则是,待内插点为中心,定义一个邻近区域,选择落入邻近区域内的采样点作为参考点,选取的参考点数n应大于等于6。常用邻近区域定义有圆形邻域和正方形邻域。34 兰州交通大学硕士学位论文(1)圆形邻域选择参考点法圆形领域是以待内插点为中心,以R为半径的圆形区域(图3.6)。在邻域范围内的参考点数随半径R的增大而增加,由于参考点的分布一般不是均匀的,所以圆的半径R是动态变化的,其大小取决于采样点的密度。图3.6圆形邻域设在圆形邻域范围内有n个参考点,A为模型区域的总面积,N为模型的总采样点数,则采样点平均密度为Nn(3.31)2AR可选择参考点数的范围为6n15(3.32)如取n=10,则由式(3.31)可得10AR(3.33)N当选择的参考点数n满足式(3.32)的范围时,可直接进行下一步的内插计算;否则,如果参考点数n小于6,则可按一定的步长扩大圆的半径,反之,如果参考点数n大于15,则要按一定的步长缩小圆的半径。在上述圆形邻域选择参考点的方法中,我们不难发现其仅考虑了参考点与圆心的距离,而未考虑参考点分布的方向,有可能出现选择的参考点密集的出现在某个方向上,而在其他方向上分布的较少或者根本没有分布参考点,这将直接影响内插计算的质量。为此我们进行了改进,就是将平面划分成4个或者8个扇形区域,并从每个扇形区域内选取距离中心最近的一两个参考点(图3.7)。-35- 基于DEM内插的工程土方量计算方法研究图3.7按方位选参考点(2)正方形邻域选择参考点法正方形邻域是以待内插点为中心,以d为变长的正方形区域(图3.8)。在正方形邻域范围内的参考点数随着边长d的增大而增加,其边长d大小取决于采样点的密度。图3.8正方形邻域正方形邻域的边长d可用下式估计10Ad(3.34)N如果采样点的坐标为(x,y),待内插点的坐标为(x,y),则当采样点的坐标满足iiddxxix22(3.35)ddyyiy22则采样点可成为内插计算的参考点。当选取的参考点数满足要求时,可直接进行内插计算,如果参考点数不满足要求,可根据情况增大或者缩小正方形的边长。当选取的参考点数大于6个时,需要按最小二乘法计算内插函数,为此需要确定每个参考点的权。不同的参考点在内插计算中起的作用是不同的。由于自然地形具有自相36 兰州交通大学硕士学位论文关性,相距较近的点,其高程相似性较强。。在内插计算中参考点距离待内插点的距离越近,所起的作用就越大,反之若距离越远,则所起的作用就越小。确定权的基本原则是:参考点的权与参考点距待内插点的距离成反比形式。移动曲面的求解是将参考点的坐标和高程代入插值函数。给每个参考点高程加一个改正数v得到误差方程i22zivia0a1xia2yia3xiyia4xia5yi(3.36)引入矩阵符号22z1v11x1y1x1y1x1y1a022z2v21x2y2x2y2x2y2a1ZVAX(3.37)znvn1xnynxnynx2y2a5nn写成矩阵形式VAXZ(3.38)按最小二乘原理求解得到T1TX(APA)APZ(3.39)其中p1P2PPn为参考点的矩阵。当n等于6时,为纯曲面内插,当n大于6时,为曲面拟合。由于需要求解复杂的误差方程,这会导致计算较为复杂。在实际应用中,更为常用的是加权平均法。加权平均法选取参考点和定权的方法跟移动曲面法都相同,所不同的是选择的内插函数为zf(x,y)a0(3.40)T此时,误差方程的系数矩阵为A(111),由此可求得nnT11T(APA)(Pi)APZpizi(3.41)i1i1根据式(3.39)可求得-37- 基于DEM内插的工程土方量计算方法研究nnzpizi/pi(3.42)i1i1上式为加权平均法内插的计算公式。38 兰州交通大学硕士学位论文4不同格网DEM内插方法在土方量计算中对比分析4.1项目区概况兰州新区位于秦王川盆地,东经103°29′22″~103°49′56″,北纬36°17′15″~36°43′29″。地处甘肃省兰州、白银两市之间,这里也是甘肃、青海、宁夏3省省会城市之间的核心位置。南北相距49公里,东西相距23公里,距兰州市38.5公里,白银市79公里,距青海省省会西宁市195公里,距陕西省省会西安市560公里,经景泰到银川有470公里,经国道312线经张掖市、酒泉市可到达新疆维吾尔自治区,距乌鲁木齐1805公里。处于黄土高原、内蒙古高原和青藏高原等三大高原的交汇中心地带,同时也处于陇西盆地与祁连山脉东延之余脉相融合的交错地带。此区域地貌类型属典型的黄土高原丘陵,地貌发育类型为河谷、梁峁、平川及沟壑。土质多为黄绵土,也是黄土高原上分布范围最广的土壤类型,土层较为深厚,表面土层平均厚度在0.5~3米之间,主要由0.25毫米以下颗粒组成,可耕性较好,自然植被稀疏。平均海拔高度1910米。属温带大陆性气候同时也是半干旱气候,春夏秋冬四季区别明显,终年阳光充足,春季多风少雨,夏无酷暑,秋季温凉,冬季寒冷干燥。年日照量1744~2659小时,日照率60%;年平均气温6.9℃,年平均降水量300~350毫米,年蒸发量1880毫米,全年平均无霜期139天,农作物一年一熟;风向方面,兰州新区主导风向为西北风,冬春多西北风,夏秋多东南风,年平均风速2.3米/秒。兰州新区综合保税区位于中川机场东侧的物流产业园区,保税区总规划建设用地面积3.39平方公里,其中办公、仓储、交易区域面积约2.86平方公里,其他附属设施服务区域面积约0.53平方公里。综合保税区区位条件得天独厚,公路、铁路及民用航空运输条件十分便利,距兰州市区也只有45公里路程,且全程为高速公路,特别是新修建的兰秦快速路为双向八车道。与新建成中(川)马(家坪)铁路新区货运站约12公里,与兰州中川机场仅仅只有2公里距离。项目区地块高程最大值1948.41米,高程最小值1930.419米,平均高程1939.89米。地势总体上呈由北向南逐渐降低的趋势,局部地区有坑和土堆。项目区地形图如图4.1所示。-39- 基于DEM内插的工程土方量计算方法研究13砖沥沥沥纬十六路沥沥砖施工区土土土施工区施工区施工区施工区施工区施工区施工区施工区水泥砖施工区施工区施工区土土施工区施工区施工区施工区砖砖经十路经七路施工区砖沥沥土施工区土堆施工区土沥土施工区土堆施工区施工区施工区土堆施工区施工区土砖砖杨5土土堆施工区施工区杨8土堆图号:T1402029-1兰州市1996中川坐标系统兰州市勘察测绘研究院1985国家高程基准,等高距为1.0米2014年2月1996年版图式图4.1项目区地形图4.2实验方法4.2.1实验数据本实验数据来源于兰州市勘察测绘研究院实测地形图同时加密高程点,坐标系为“兰州市1996中川坐标系统”。高程基准为“1985国家高程基准”。运用网络RTK建立控制点,平面精度为一级,高程精度为图根水准,运用全站仪采集碎步点。4.2.2实验用内插方法(1)反距离加权插值法反距离加权插值(InverseDistanceWeighted)以其方法简单容易操作等特点在各种插值方法的研究中应用较为普遍,属于逐点内插法。此方法基于相近相似的原理:即两个物体离得越近,它们的性质就越相似,反之,离得越远则相似性越小。它以插值点与样本点间的距离为权重进行加权平均,离待插值点越近的样本点赋予的权重越大。[30]是从点创建面来测定每个位置的数值。反距离加权插值算法是一种精确性插值算法[31]。设有一系列离散点分布于平面上,已知其坐标和值为X,Y,Z(i1,2,,n),根据已iii40 兰州交通大学硕士学位论文知附近其他离散点的值,用距离加权的方法来估计Z点值,则有nzi2i1diz(4.1)n12i1di其中222dXXYY(4.2)iiiIDW插值通过对包含在某点周围一定范围内的区域的每个采样点值进行取算术平均值,用这个算术平均值作为内插单元值。整个运算过程要求必须保证均分,要求离散点应该尽量均匀分布而不聚集在某个方向上的位置,并且密集程度必须要达到能够真实和最大程度表达局部重要细节特征的变化。在GIS软件中,针对点数据插值的方法有很多,最常见的是IDW插值方法。用这种方法产生的计算结果通常会出现诸如孤立点数据明显高于周围数据点分布模式的情况,这是因为其计算值受到了数据点集群的影响。我们可以运用动态搜索准则在插值过程中[32]对其进行某种程度的改进。(2)样条插值法样条插值法属于分块曲面内插,它使用可最小化整体表面曲率的数学函数来估计[33]值,以生成恰好经过输入点的平滑表面。当使用规则样条函数法来进行插值时能够创建一个逐渐变化的较为平滑的表面,而这里使用的数据可能不包括在样本数据范围中。根据建模现象的不同特性,张力样条函数法能够控制模型表面的硬度,以此来创建不太平滑的表面,其使用的数据受到样本数据范围更为严格的约束。则有nS(x,y)T(x,y)Rr(4.3)iii1式中,i1,2,,n,n为采样点的数目,是通过对线性方程组进行求解而解算出的系i数,r是点X,Y与第i点之间的距离。i对于规则样条函数来说T(x,y)aaxay(4.4)123其中,a是对线性方程组进行求解后算出的系数。i-41- 基于DEM内插的工程土方量计算方法研究21rr2rrR(r)lnc1K0cln(4.5)24222其中,r是点与样本间距离,是权重参数,K是修正贝塞尔函数,c是一个固定0常数,其值为0.577215。对于张力样条函数来说T(x,y)a(4.6)1其中,a是对线性方程组进行求解后算出的系数。11rR(r)2lncK0(r)(4.7)222其中,是权重参数,r是点与样本间距离,c是大小等于0.577215的常数,K是0修正贝塞尔函数。(3)克里金插值法克里金插值法又称空间自协方差最佳插值法,它是以南非工程师克立格(D.GKrige)[34]的名字命名的一种内插方法。此方法属于逐点内插法,根据样本空间几何位置和样本之间彼此相关程度大小的不同,分别对每个样本赋予不同大小的权重,接着运用滑动加权平均的方法来对未知样点上样本平均值进行估计。目前,在前人学者的文献中将克里[35]金插值方法与其它插值方法进行比较研究的依然比较少。克里金法假定空间相关性可以在采样点之间的距离或方向上被反映。它可以被用来解释表面变化。为了最终得到每个位置的最后输出值,克里金法拟合了数学函数与提前设置好数量的点或已经设置好采样半径内的全部点。克里金法是一个多步过程,主要有对于数据的探索性统计研究分析、变异函数表面建模和表面创建,同时还有方差表面的研究。其函数为nZs0iZsi(4.8)i1其中,s为预测位置,n为测量值数,Z(s)为第i个位置处的测量值,为第i个0ii位置处的测量值的权重且这个权重是不确定的。本实验中采用普通克里金法,半变异模型为球面函数。搜索半径点数为12。(4)自然领域插值法自然领域插值法首先寻找距离查询点最近的输入样本子集,然后根据区域的大小不同按照比例对样本赋予一定的权重之后再完成插值操作。其特点是搜索范围具有局部42 兰州交通大学硕士学位论文性,仅使用查询点周围的样本子集作为操作对象,同时保证了插值不会在样本范围之外进行。在自然领域插值法中所有点的自然领域都与与其邻近的Voronoi图(又叫泰森多边形)有关。初始Voronoi图由离散点构成,加入待插点后就会重新建立泰森多边形,前后创建的不同Voronoi图的重叠比例将被用作权重。4.2.3精度评价方法对DEM常用的精度评定方法有检查点法、剖面法、等高线法等。本实验采用检查点法进行精度分析。(1)中误差精度模型中误差(Rootmeansquareerror,RMSE)作为最常用的精度度量指标,其前提假设是误差具有随机性,即误差必须服从均值为0的正态随机分布。通常情况下,误差均值并不一定为0,所以中误差不能揭示误差成分中的系统误差,另一方面,中误差也无法表示单个误差的大小,它只是从整体上描述数字高程模型在多个点处的近似值与其真值之间的离散程度。中误差计算公式为:n2(ZZ)DEMki1(4.9)n其中,Z是通过DEM内插得到的高程值,Z是通过先进的测绘仪器和技术测量DEMk得到的较高精度的高程值,即真值,n为采样点数目。也有学者在一些文献中提出用其他更为全面的精度模型来度量,如平均误差(meanerror,ME)和标准差(standarddeviation,SD):n(ZZ)DEMki1ME(4.10)nn2ZDEMZkMEi1SD(4.11)n但是平均误差和标准差也存在一定的缺陷和不足,因此在本实验中仍旧采用中误差指标对实验结果数据进行评价,用来判断不同内插方法和不同数据网格间距建立的DEM的精度的高低。(2)检查点选择使用中误差模型计算得到的DEM精度值明显受到检查点数的影响。太少的检查点会-43- 基于DEM内插的工程土方量计算方法研究导致不可靠的精度结果,太多的检查点数将引起不必要的成本投入。因此在相应的置信[36]水平下选择适宜的检查点数是必须的。我国国家测绘局1:1万数字高程模型生产技术规定中采用检查点的方式对精度进行[37]检测,用28个检查点对图幅内和图幅边缘进行检测。本实验数据中实验区面积约为0.9平方公里,大约有16000多个外业实测高程点。选取高程点总数3%的点(约490个)作为检查点。4.3实验过程实验选取了反距离加权插值法(IDW)、样条插值法(Spline)、克里金插值法(Kriging)、自然邻域插值法(NaturalNeighbor)的参数中不同数据网格间距作为研究对象,各个插值法网格间距的实验取值及DEM插值精度如表4.1所示。表4.1DEM插值算法中不同数据网格间距精度网格间距(m)IDWSplineKrigingNaturalNeighbor10.4764355310.5676625610.4059572380.40783042420.4829510610.5509105690.4208853480.42673653340.4825952710.598332850.4484277050.44464431360.5922363830.6867082270.5693149490.56829138680.7137032320.8038707050.6982727660.701694974100.7207179360.8472253010.7181344910.716499906120.7252920890.8571853120.7291338970.716583867140.7743769111.2079990570.7812330260.78019592160.8741717211.9149405820.9222439320.916427453四种内插方法的精度对比如图4.2所示。网格间距越小,精度越高。由于在土方计算过程中我们无法获得真实的土方量,所以只能选择相对精确的数值作为真值来评价不同内插方法所计算得到的土方量的精度。由于实验中我们已经发现Spline法构建的DEM精度较低,同时在计算中我们也能发现其土方量与其他三种方法计算得到的结果差距很大,因此在本实验中我们选择IDW法、Kriging法和NaturalNeighbor法三种方法中1m网格所计算出的土方量的平均值作为真值来计算其他间距的网格所构建DEM在土方计算中的相对误差。[38]相对误差计算公式为:44 兰州交通大学硕士学位论文|V-V|真%(4.12)V真其中,V是每种内插方法构建DEM所计算得到的土方量,V是三种方法中1m网格真内插DEM所计算得到的土方量的平均值。0.9IDWSpline0.85KrigingNaturalNeighbor0.80.750.70.65中误差0.60.550.50.450.40246810121416数据网格间距m图4.2不同内插方法构建DEM精度对比分别用不同内插方法得到的DEM与设计高程所建立的DEM叠加计算得出土方量及其精度(不计算Spline法精度)如表4.2所示。三种内插方法计算土方量精度对比如图4.3所示。在本实验中,虽然用三种方法中1m网格所计算出的土方量的平均值作为真值能够对土方计算精度进行评价。但这种方法产生的结果仍然不能完全令人信服。同时,外业测绘人员在数据采集过程中针对项目区场地不同地形而采取不同采集密度,例如在地势较为平坦区域高程点采集较为稀疏,在不平坦区域如坎上砍下、坡顶坡底分别较密集采集并赋予不同点名。在第二章的研究中我们已经发现由于不规则三角网能很好的顾及实际地形特征,因此不规则三角网法在土方计算中的精度很高。考虑以上三个方面的问题,我们首先对外业采集到的16000多个点先进行构建三角网的处理。接着,根据编绘好的地形图以及密集的坎上砍下、坡顶坡底高程点对自动构建的三角网进行详细的人工修改,使自动构建的三角网最大程度地符合实际地形。将此三角网转栅格数据后进行土方-45- 基于DEM内插的工程土方量计算方法研究计算。将此土方量的值作为真实值来计算其他三种方法中不同采样间距的网格所构建DEM在土方计算中的相对误差。表4.2不同插值方法计算土方量及其精度网格IDWSplineKrigingNaturalNeighbor间距(m)3333土方量(m)精度土方量(m)土方量(m)精度土方量(m)精度11062462.560.010031149002.231082403.760.008551074812.780.0014821062423.240.010071149480.231082311.740.008471074796.410.0014641058613.930.013621149224.381078442.920.004861070900.080.0021761062762.250.009751150861.611083479.820.009551075939.550.0025381066164.910.006581164488.591085072.590.011041077704.670.00417101067912.280.004951165569.011088426.770.014161080077.20.00638121061999.240.010461155619.631081773.860.007961074422.760.00111141050336.890.021331248374.631069911.670.003091062209.730.01026161042369.560.028751208805.281057319.130.014821050588.220.021090.03KrigingNaturalNeighbor0.025IDW0.020.015相对误差0.010.00500246810121416数据网格间距m图4.3不同插值方法计算土方量精度对比46 兰州交通大学硕士学位论文表4.3不同插值方法计算土方量及其精度网格IDWSplineKrigingNaturalNeighbor间距(m)3333土方量(m)精度土方量(m)土方量(m)精度土方量(m)精度11062462.560.017691149002.231082403.760.000751074812.780.0062721062423.240.017731149480.231082311.740.000661074796.410.0062941058613.930.021251149224.381078442.920.002921070900.080.0098961062762.250.017411150861.611083479.820.001741075939.550.0052381066164.910.014271164488.591085072.590.003211077704.670.00360101067912.280.012651165569.011088426.770.006311080077.20.00140121061999.240.018121155619.631081773.860.000161074422.760.00663141050336.890.028901248374.631069911.670.010801062209.730.01792161042369.560.036271208805.281057319.130.022451050588.220.02867三种内插方法计算土方量精度对比如图4.4所示。0.04IDW0.035KrigingNaturalNeighbor0.030.0250.02相对误差0.0150.010.00500246810121416数据网格间距m图4.4不同插值方法计算土方量精度对比-47- 基于DEM内插的工程土方量计算方法研究4.4实验分析与结果理论上讲,网格间距越小,所构建DEM精度越高。图4.1中的实验结果可以清楚地证明这一理论。不论是哪一种内插方法,随着网格间距的增大,中误差也越大,说明精度越低。从图4.2可以看出:(1)Spline法不论在哪一种网格间距中的中误差都比其他三种内插方法大,说明Spline法构建的DEM在整体上精度较低,不可取;(2)IDW法在网格间距较小时的中误差高于Kriging法和NaturalNeighbor法,此时其精度低于上述两种方法。之后其中误差与Kriging法和NaturalNeighbor法基本相同,在网格间距较大时,其中误差又高于Kriging法和NaturalNeighbor法,说明IDW法构建的DEM在不同网格间距的精度不稳定;(3)Kriging法和NaturalNeighbor法的中误差在各个数据网格上都非常接近,表明这两种方法精度都非常接近。从图4.3可以看出:(1)IDW法的相对误差在4m网格间距时略有增大,之后相对误差减小,在10m网格间距时达到最小,之后相对误差急剧增大;(2)NaturalNeighbor法的相对误差逐渐增大,在12m网格间距时相对误差达到最小值,之后相对误差也急剧增大;(3)Kriging法的相对误差随着网格间距大小的变化而在某个数值附近连续波动,既没有突增,也没有突减。其相对误差保持相对的稳定;(4)IDW法和NaturalNeighbor法的这种相对误差突增突减的情况表明其在土方计算过程中有不稳定性,选取不同网格间距时会造成土方计算差别较大,因而不可取;(5)Kriging法的相对误差保持相对稳定,在选取不同网格间距时其土方计算量不会相差太大,这表明Kriging法在计算土方量时受网格间距大小的影响不大。从图4.4可以看出:(1)IDW法的相对误差在各个像元网格大小中均高于Kriging法和NaturalNeighbor法,表明其土方计算的精度低于Kriging法和NaturalNeighbor法:(2)Kriging法的相对误差仅在像元网格大小为10m处略高于NaturalNeighbor法,而在整体趋势上始终低于NaturalNeighbor法。这表明Kriging法计算土方的精度高于IDW法和NaturalNeighbor法。综合以上实验可以认定在场地范围较大、地势总体较为平坦并且高程点足够多的情况下,Kriging法在计算土方量的方法中是最优的。48 兰州交通大学硕士学位论文结论伴随着建设工程项目规模越来越大,前期土地整理工作也变得越来越复杂和重要。不同土方计算方法适用性也不同,研究不同地形地貌下土方量计算精度最高的方法对于提高计算效率和控制工程施工成本具有重要的现实意义。本文以兰州新区综合保税区土方计算项目为研究对象,主要完成了以下几个方面的工作和成就:(1)对土方工程的特点和土方计算方法研究的背景和必要性进行了说明。着重分析了国内外专家学者对于常规土方计算方法和基于DEM内插的计算方法的研究现状。提出了常规计算方法的弊端和基于DEM计算的优势,对研究土方量计算精度更高方法的目的和意义进行了说明。(2)对常规方法中的断面法、方格网法、不规则三角网法和等高线法进行了归纳对比,从计算原理和具体计算步骤等方面进行了详细的阐述和分析,对比了各自适用情况,总结了各个方法在土方计算中存在的问题,对下一步基于DEM内插方法的研究提出了展望。(3)研究了数字高程模型的理论,包括概念、特点、分类、数据采集方法以及建模的基本方法。对DEM内插的概念、原理和分类进行了阐述。分析了DEM内插方法的数学模型和内插过程。为之后的实验奠定了理论和数学基础。(4)在前面理论研究的基础上,以兰州新区综合保税区土方计算项目为案例进行了实验分析。应用四种不同DEM内插方法建立DEM,对建立的DEM精度使用中误差精度模型进行评价,剔除了其中一个DEM精度较低的内插方法。同时,对其他三种方法计算得到的土方量精度使用相对误差精度模型进行评价,结果证明在实验区计算土方量的方法中Kriging法是最适宜的。本文虽然通过实验得到了在实验区这种地面高差不大的区域进行土方计算的最优方法。但由于实验条件和实验数据的的限制,本文的研究还存在以下不足:(1)DEM内插是一个复杂的概念,由于地形表面通常较为复杂且多变,同时DEM精度通常与数据点获取的方法、分布和数量等因素有关。不同地形都有适用于它的不同内插方法,没有哪种内插方法是绝对的好与坏。而本实验数据点的获取是通过外业RTK做控制点,全站仪进行碎步测量而得到的,在如此大面积场地中进行高密度碎步点测量是极其耗费人力的,如果现场没有条件进行如此细致的测量,那么实验结论是否会有不同,这一点需要进一步研究。(2)在数字高程模型的建立中,通常会忽略一些特殊地貌。在某些场地中会出现一些比较大的坑,这就牵扯到陡坎的问题,在外业数据采集中对陡坎处的坎上坎下都需要采集。不同内插方法构建DEM在陡坎处是如何处理的,是否会影响最终实验结果仍然需要进一步研究。-49- 基于DEM内插的工程土方量计算方法研究(3)由于无法获取工程的真实土方量,在本实验土方量相对误差的计算中,选择IDW法、Kriging法和NaturalNeighbor法三种方法中1m网格所计算出的土方量的平均值作为真值去评价土方计算精度,这在某些情况下是不合理的,究竟在什么情况下不合理,不合理性在哪,这一点仍然需要讨论。50 兰州交通大学硕士学位论文致谢三年的研究生生活转眼间即将结束,在兰州交通大学学习和在兰州市勘察测绘研究院实习期间的经历令我终生难忘,在这里谨向给予我关心、帮助、支持的老师、同学、同事、朋友致以由衷的感谢!首先特别感谢我在兰州市勘察测绘研究院的导师魏忠邦老师。魏老师,无论是从学术上、工作上还是生活上都是我学习的典范。作为院长,魏老师平时需要主持、安排全院工作。工作非常繁忙。即使在这样的情况下,魏老师仍然在工作之余抽出时间在我的学业上给予悉心指导。从课题的选题、资料的收集、研究方法的确定到论文定稿,魏老师都倾注了大量的时间和精力。在论文的研究和撰写过程中,兰州市勘察测绘研究院的张永忠高工、张剑锋高工、雒喜平高工、槐升团高工、李玉宝高工、王勇军科长都给予了很大的关心和帮助,使本论文的工作得以顺利完成,在此表示衷心的感谢!在兰州市勘察测绘研究院测绘队实习工作期间,我的同事薛道义、郑涛、王祥龙在外业测量工作中给予了我很大的帮助。在地理信息室实习工作期间,张清彦主任、刘小栋副主任、廖文斌、吕江波、张河、唐志龙、马彩云、杨倩也为本论文提出了很多宝贵的建议,在此表示衷心感谢。最后,特别感谢我的父母和亲朋好友,他们的深切关怀和大力支持是我能够顺利完成学业的动力。马龙二零一五年五月-51- 基于DEM内插的工程土方量计算方法研究参考文献[1]刘鸿剑,肖伟红,何建英,等.基于DEM的工程土方量算法研究——以抚州市人工湖项目堤防工程土方量计算为例[J].江西理工大学学报,2009,30(4):17-21.[2]周越轩,刘学军.基于DTM的土方工程计算与精度分析[J].长沙交通学院学报,2000,16(4):39-43[3]吴小燕,陈玉莹.基于复杂区域土方计算的研究[J].矿山测量,2009,5.[4]邹建喜,陈永涛,孙秀丽,等.一种基于样条插值的土方计算方法[J].青岛建筑工程学院学报,2004,25(1):7-9.[5]胡振琪,高永光,李江新.ERDAS在土地整理土方量计算中的运用[J].中国土地科学,2006,20(1):50-54.[6]俞志新,李艳,黄明祥.地统计克立格插值法在工程土方计算中的应用[J].浙江水利科技,2003(4).37-38.[7]梁晓鹤.基于MTA计算土方量的技术研究[J].城市勘测,2010,3.[8]张敏捷.基于EPS平台的土方测算方法探讨——以增城市某人工湖开挖工程为例[J].测绘与空间地理信息,2013,36(5):205-207.[9]柯晓山,张玮,王荣静.采用不规则三角网插值进行土地整理项目前期平整土方量的计算[J].农业工程学报,2004,20(3):243-247.[10]柳长顺,杜丽娟.Arcview在土地整理项目土方量计算中的运用[J].农业工程学报2003,19(2):224-227.[11]GoktepeA.Burak,LavA.Hilmi.MethodforBalancingCut-FillandMinimizingtheAmountofEarthworkintheGeometricDesignofHighways.JournalofTransportationEngineering,Sep.,2003[12]MawdesleyMichaelJ.,AskewBillH.,PattersonDanE.,etal.AcomputerSystemformodelingtheprocessofearthworkPlanninginlinearconstructionProjects.InternationalJournalofComputerApplicationsinTechnology,2004.[13]HarmananiHaidarM.,ZoueinPierretteP.,HajarAouniM..AParallelGeneticAlgorithmfortheGeometricallyConstrainedsitelayoutproblemwithUnequal-sizeFacilities.InternationalJournalofComputationalIntelligence&Application,Dec.2004.[14]KIDNERDB.Higher-orderInterpolationofRegularGridDigitalElevationModels.InternationalJournalofRemoteSensing,2003,14:2981-2987.[15]REESWG.TheAecuracyofDigitalElevationInterpolatedtoHigherResolutions.IntemationalJournalofRemoteSensing,2000,21(l):7-20.[16]VincentC.Accuracyofinterpolationtechniquesforthederivationofdigitalelevationmodelsinrelationtolandformtypesanddatedensity.Geomorphology,2006,77(1):126-141.[17]张光辉.快速计算土方量的方法[J].测绘通报,1997,(5):23-24.[18]李书全.土木工程施工[M].上海:同济大学出版社,2004.[19]TiejunDeng,HanzhiHe.Deviationanalysisofearthworkcalculatedbysquaregrimethod.In:ProceedingofTheSecondInternationalConferenceonManagementScienceandEngineeringManagement.UK:WorldAcademicUnion,2008,594-600.[20]陈黎阳.土方测量计算方法比较研究[J].现代测绘,2010,33(5):36-38.[21]季朝亮,李宗聚,马学民.关于几种土方量计算方法的研究[J].测绘与空间地理信息,52 兰州交通大学硕士学位论文2010,33(3):219-222.[22]罗德仁,邹自力,汤江龙.工程土方量计算比较分析[J].华东理工学院学报,28(1):59-63.[23]李志林,朱庆.数字高程模型[M].武汉:武汉测绘科技大学出版社,2000.[24]周秋生,刘妍,马俊海.数字高程模型及其应用[M].哈尔滨:哈尔滨工程大学出版社,2012.[25]汤国安,刘学军,闾国年.数字高程模型及地学分析的原理与方法[M].北京:科学出版社,2005.[26]胡海,游涟,胡鹏,彭会琨,数字高程模型内插方法的分析和选择[J].武汉大学学报,2011,36(1):82-85.[27]王玲.多种数字高程模型内插方法实验结果分析研究[J].测绘与空间地理信息,2013,36(5):185-188.[28]朱求安,张万昌,余钧辉.基于GIS的空间插值方法研究[J].江西师范大学学报,2004,28(2):183-188.[29]寇程,柯长青.地形平坦地区DEM生成算法的比较研究[J].测绘与空间地理信息,2013,36(7):33-37.[30]廖伟华,徐彬.基于空间内插的场地平整土方工程量计算方法研究[J].黑龙江工程学院学报2007,21(4):14-17.[31]张锦明,游雄,万刚.DEM插值参数优选的试验研究[J].测绘学报,2014,43(2):178-185.[32]彭楠峰.距离反比插值算法与Kriging插值算法的比较[J].大众科技,2008,5.[33]徐潇,谭衢霖,王浩宇,等.复杂地貌地形图等高线内插算法的精度分析[J].遥感信息,2013,28(6):111-115.[34]曹俊茹,刘强,姚吉利,等.基于Kriging插值DEM的计算土方量方法的研究[J].测绘科学,2011,36(3):98-99.[35]靳国栋,刘衍聪,牛文杰.距离加权反比插值法和克里金插值法的比较[J].长春工业大学学报,2003,24(3):53-57.[36]LiZhilin.EffectsofcheckpointsonthereliabilityofDTMaccuracyestimatesobtainedfromexperimentaltests[J].PhotogrammetricEngineeringandRemoteSensing,1991,57(10):1333-1340.[37]周兴华,姚艺强,赵吉先.DEM内插方法与精度评定[J].测绘科学,2005,30(5):86-88.[38]曹俊茹,刘强,姚吉利,等.基于Kriging插值DEM的计算土方量方法的研究[J].测绘科学,2011,36(3):98-99.-53- 基于DEM内插的工程土方量计算方法研究攻读学位期间的研究成果马龙,魏忠邦,王勇军.基于VC++的城市地下管线数据检查程序开发.中国科技人才,2014,10.54'