多面拟合函数法转换GPS高程

桂林理工大学博文管理学院设计论文XXIII摘要传统的几何水准测量方法,是测绘领域中测定正常高的主要方法。这种方法虽然精度高,但实施起来费时费力,作业效率较低。GPS定位技术自问世以来,就以其精度高,速度快,操作简单等优势对传统的水准测量造成了极大冲击，但是现实测绘和工程建设上需要的具有物理意义的正常高（或正高）高差，因此，如何用GPS高求定正常高，是当前测绘工作者较为热心的课题之一。关键词：GPS水准；大地高：正常高Abstract:Traditionalgeometriclevelingmethod,amajormeanstomeasurethenormalheightinsurveying,produceshighprecisionbutcantaketoomuchtimeandeffortswithlowefficiency.GPStechnologyhashithardthetraditionallevelingwaysviaitshighprecision,fastspeedandsimpleoperationsinceitcametotheworld.However,reallevelingandprojectconstructionrequireanormalaltitudedifferencewithphysicalsignificance,so,howtoidentifythenormalheightwiththeGPSheightisoneofthepreoccupationsforsurveyingworkers.Keywords:GPSleveling;geoidheight;normalheight目录摘要 I1引言 V2高程系统 V2.1大地高系统 V2.2正高系统 VI2.3正常高系统 VI2.4高程系统之间的转换关系 VII3GPS高程测量的基本原理 VII3.1物理大地测量方法 VII3.2几何方法 VII4GPS大地高转化为正常高的方法讨论 VIII4.1平面拟合 VIII4.2多项式曲线拟合 VIII4.3二次多项式曲面拟合法 VIII3.4多面函数曲面拟合 VIII4.5加权平均值法 IX3.6非参数回归法和高程异常变化梯度法 X4.7移动曲面法 X4.8固定边界3次样条插值法 XI4.9非格网GPS散点数据考虑地形改正法 XI4.10BP神经网络算法 XII5算法实例 XIII5.1参与计算的数据 XIII5.2平面拟合法 XIV5.3二次曲面法拟合 XV5.4加权平均值法 XVI5.5多面函数曲面拟合 XVI6结论 XVI参考文献 XVII附录1.平面拟合程序： XIX附录2加权平均值法拟合程序： XXI附录3二次多项式曲面拟合程序 XXII附录4多面函数曲面拟合程序 XXIII致谢 XXIV1引言我国目前采用的高程系统为正常高系统:即以似大地水准面为参考面的高程系统,确定高程通常采用的是几何水准。GPS的广泛应用为确定高程提供了新的途径。但是,用GPS测定的三维坐标,大地经度、大地纬度和大地高都是以WGS84参考椭球面为基准的,GPS测定的大地高是一个几何量,与地球的形状和重力场的分布没有关系。而我们实际在绝大多数情况下使用的为正常高,它是一个物理量。因此,实际使用时,需要将GPS大地高转换为正常高。我们知道,GPS测得是三维坐标,实际上主要利用的是其平面位置信息,高程信息没有充分利用。其原因主要在于现有的各地区的高程异常值精度较低,而且GPS大地高本身的精度也低于平面位置的精度。但若把GPS大地高与现有的水准资料相结合,可对局部地区的似大地水准面进行精化,从而得到较精确的高程异常,获得我们需要的高程。在实际应用的作业的时候，因为区域之间似大地水准面精化并没有普遍被使用，所以我们在考虑范围不大的区域的时候高程异常具有一定几何相关性，几何的高程拟合方法也没有被我们在系统中广泛的被应用。而GPS高程拟合是通过利用这个原理，采用数学的方法，求解GPS观测点的正常高度，在GPS高程拟合模型当中GPS点布成一定区域面的时候我们可以采取数学曲面拟合的方式或方法来求得待定点的正常高。根据观测区中所得知的点平面坐标x与y的高程异常值，采用数值拟合方法，在拟合处测量区域似大地水准面在内插出待求点的高程异常值，从而得求出待求点正常高。。因此,GPS测量正常高的研究不仅具有一定的理论意义,而且具有非常重要的现实意义,有着广阔的应用前景。2高程系统2.1大地高系统大地高是地面上某点沿通过该点的椭球面法线到椭球面的距离,以椭球面为基准的高程系统为大地高系统。通常用H表示。因为椭球面在地球内部的位置，决定于椭球体的定位参数△X0=(△X0△Y0和定向参数ω=（ωXωY所以，当椭球的定位与定向不同时，相应的大地高程系统也是不同的。所以δH表示不同大地高程系统的高程差，则有关系式：δH=+（-12Ne2大地高是一个几何量，它不具有物理上的意义。利用GPS定位技术，可以直接测定观测站在WGS-84中的大地高程。这一高程系统，在工程技术上虽然没有广泛应用，但是它在与水准测量资料、重力测量资料相结合，来研究大地水准面的形状方面，或在结合高程异常资料，来确定点的正常高方面，都具有重要意义。2.2正高系统正高高程系是以大地水准面为高程基准面，地面上任一点的正高高程（简称正高），即该点沿垂线方向至大地水准面的距离。通常以Hg如图一中一点B的正高高程可由下式求得。（2-2-1）设沿垂线BC的重力加速度用gB表示，在垂线BC的不同的点上gB也有不同的值。由位能差相等可得：gBdH=gdh（2-2将dH带入2-1式中，得：（2-2-3）如果取垂线BC上重力加速度的平均值为gmB（2-2-4）式中gmB为常数，gdh为过B点的水准面与大地水准面之间的位能差，也不随路线而异，因此，正高高程不随测量路线的不同而有差异。如果沿着水准路线每隔若干距离测定重力加速度，则上式中的g值是可以得到的。但是由于沿垂线BC的重力加速度gB不但随深入地下深度不同而变化，而且还与地球内部物质密度的分布有关，所以重力加速度的平均值图1-1正高求解因为正高是以大地水准面为基准面，具有重要的物理意义，所以它在水利建设、管道和隧道建设等精密工程技术方面，有着广泛应用。2.3正常高系统为了克服求得正高的值所遇到的困难，莫洛金斯基（М.С.Молоенский）于1945年提出了正常高的概念，将正高系统中不能精确测定的gmB用正常重力γm2.4高程系统之间的转换关系高程位置用大地高H或正高Hg或正常高Hγ表示。如图二所示，大地高H是地面点沿法线到椭球面的距离PP0；正高H正是地面点沿实际重力（垂）线到大地水准面的距离PP1；正常高H正常是地面点沿正常重力（垂）线到似大地水准面的距离PP2H=Hg+N（H=Hγ+ξ（2式中：N为大地水准面差距，ξ为高程异常。图1-2高程系统之间的关系3GPS高程测量的基本原理GPS测量能够精确地给出地面点在WGS-84坐标系中的三维坐标X,Y,Z或B,L,H,经系统变换可以得到地面点在局部坐标系中的大地高。由于各GPS点上的高程异常值无法直接获得,目前还无法直接将大地高精确地转换成实用的海拔高。因此高程异常的确定成为GPS高程转换的关键。目前有物理大地测量方法和几何方法可以确定高程异常。3.1物理大地测量方法利用全球或局部重力模型加上重力测量数据和地形数据按物理大地测量方法计算正常高，即ξ=ξGM+ξ△g+式中ξGM为地球重力模型，ξ△g为重力异常，ξT表示地形起伏对大地水准面的影响。对于一般工程单位,无法获得必要的重力数据,故这种方法难以普及和推广,但很有研究价值。3.2几何方法假设在测区内有若干个既进行了GPS测量又联测了水准高程的GPS点(这样的点称为水准重合点),即可以利用大地高和高程异常之间的关系,推算出各水准重合点上的高程异常,利用这些离散点上的异常值,可以拟合出测区所在局部区域的似大地水准面，进而可以内插出未知点上的高程异常,从而求出其正常高。4GPS大地高转化为正常高的方法讨论4.1平面拟合在很小的范围内,可以认为大地水准面趋于平面,表达式为：ξi=a0+a1xi式中,x、y为高斯平面坐标,a0、a1、a2为待定系数。由已知点的正常高和GPS大地高确定平面的3个参数a0、a1、a2,。求得第i个拟合点的高程异常。采用3点拟合大地水准面时,确定平面的3个已知点与拟合点的距离尽可能近;4.2多项式曲线拟合设测点的ξi和Xi(Yi或拟合坐标)存在如下函数关系：ξx=a0+a1x1+a2采用此种方法拟合似大地水准面,拟合范围越大,高程异常的变化越复杂,削高补低的误差也越大。同时,随着多项式阶次的增高,拟合出的曲面的震荡增大,出现Runger（高次插值龙格现象）现象。4.3二次多项式曲面拟合法二次曲面函数拟合的基本思想是:在区域GPS网内,存在这样一些起算点(既进行了GPS测量又进行了几何水准联测的点),这些点的平面坐标或者地理坐标已知,也知道这些点的大地高和几何水准正常高。这样我可以假设将似大地水准面看成一个二次曲面,将高程异常表示为平面坐标(x,y)或(B,L)的函数。这里只介绍平面坐标的情况,地理坐标与此类似。已知的高程异常确定测区的似大地水准面形状,求出其余各点的高程异常,然后根据高程异常和大地高的关系,求出其他点的正常高。其数学模型为：ξ=fx,y+ε式中,f(x,y)是拟合的似大地水准面;ε是模型误差。将f(x,y)按照次数升序排列展开的具体形式是:fx,y=a0+a式中,a0、a1、a2ξ=a0+a1方程中有6个待定参数,为了求得六个参数,至少需要有六个已知点。当测区的点多于六个的时候,利用最小二乘原理求得最佳拟合参数,即满足ε2=min。利用计算出的拟合参数和未知点的坐标值,计算未知点的高程异常ξ3.4多面函数曲面拟合多面函数法是从几何观点出发,解决根据数据点形成一个平差数学曲面问题。其理论根据是认为“任何数学表面和任何不规则的圆滑表面,总可以利用一系列的有规则的数学表面的总和,以任意精度逼近”。多面函数方程在笛卡尔坐标系中的一般形式为:Z=fx，y=j=1najQ(x,y,xj,式中aj为待定系数，Q(x,y,xj,yj)是x核函数可以任意选择,为了简单,一般都采用具有对称性的距离型,即Qx,y,xj,yj=[x-xjQx,y,xj,yj=[x-(4-4-2)、(4-4-3)两式分别称为正双曲面函数和倒双曲面函数。式中的δ为圆滑因子,也称光滑系数。通常δ可取一小正数或零。将式（4-4-1）中的Z用高程异常ξ代替，选倒双曲面函数为核函数则ξ=fx，y=j=1na应用多面函数拟合内插高程异常,其中有3个问题需要进一步解决:（1）核函数的选择。多面函数是属于纯数学逼近方法,其差值精度与选用的核函数有密切的关系,高程变化大的地区内插的精度往往较低,所以我们要选择一个一成不变的核函数对各种地形进行最小二乘推估,首先应该顾及高程变化大的地区。只要高程变化大的地区保证了精度,对于平坦地区的推估,其核函数影响是不大的。对于待求点的估值实质上就是数据点的加权平均值。（2）圆滑因子δ2的选择。优化选取δ2对提高拟合效果有作用,但比较困难。对于倒双曲线函数,δ2（3）函数结点的选取。如果在一定范围内有较多的GPS水准点,即已知高程异常的点较多,可选取其中部分点为结点,其余点作为拟合高程的检核点。此时拟合结果就与所选结点不同而异,这是需要进一步实验和研究的问题。对于局部GPS网,GPS水准点不会很多,可以全部用来做结点,进行拟合。4.5加权平均值法大地水准面起伏是复杂的，因此高程异常值的变化也是复杂的。但大量的实验表明，相邻点之间高程异常值的变化基本上是一个连续的渐变的过程。也就是说，两个点相距越近，其高程异常值间的联系就越为紧密，亦即:当用一点周围各点的高程异常值估算该点高程异常时，周围点距该点越近，其高程异常值对该点高程异常的影响越大，距该点越远，影响越小。基于以上思考，笔者提出了权的问题，即认为在用某点周围若干点的高程异常推算该点高程异常时，周围各点高程异常值在推算时所占的权重与各点到该点的距离成反比。假设有一点A，其周围有n个点的高程异常值已知为ξi，这些点到A点的水平距离为Li，则各点对A点高程异常影响的权重PPi=1Li并且A点的高程异常ξAξA=PξP以上未加权平均值的基本思想。通过资料分析与野外对比实验表明，式（4-5-1）若改为下式将更为适宜：Pi=1Lin式(4-5-3)中n（1，2，3，…）。一般情况下，随着n值的增大，以的计算效果与实际值的吻合程度将不断提高，通常取Pi=采用加权平均值法推算高程异常进而实现GPS大地高与正常高的转换时，必须使已知点沿控制网外围比较均匀地分布，使推算点位于已知点所围成的多边形以内，否则不能保证推算点计算结果的可靠性。另外还应使多边形范围内有一定数量的已知点，并尽可能地均匀分布。一些实验与分析资料显示，控制网中已知点的数量越多，密度越大，均布情况越好，利用加权平均值法求得的高程异常值的准确度就越高(即GPS大地高转换为正常高的实际精度越高)。一般要求一个网中已知点的数量要维护在10个左右，至少为6个。这是加权平均值法的重要应用条件之一。3.6非参数回归法和高程异常变化梯度法非参数回归法是一种广义的回归方法。其原理为：w（4-6-1）设wnμnx=i=1nwni为回归函数产μnx的权函数沽计，它具有思路直观、模型宽松、计算简单的优点。它的理论基础是概率密度估计的核估计和最邻近估计。该方法的关键在于权函数的适当选取。近邻权是一种具有优良大样本性质的权。一些实践和理论表明,即使是从全部已知点中选用距未知点最近的一个已知点所对应的高程异常值作为该未知点的预测值,其风险也只是最小风险的2倍。高程异常变化梯度法是首先估计出测区范围内高程异常变化的总体趋势,然后选取距待求点最近的已知点上的高程异常值作为待定点高程异常的平滑顶,再考虑待定点的高程异常的波动值,获得未知点上的高程异常。4.7移动曲面法移动曲面拟合法是指用户规定一个有限区域,该区域的位置将随着拟合点的位置变化而移动。选择这种模型的主要原因是可以更好地模拟似大地水准面,而不至于因远离已知点的负面效应损害其精度。在GPS高程拟合中,拟合对象一般是指高程异常。设离散点的高程异常与其平面坐标x、y存在的关系为ξ=fx,y+ε（式中,f(x,y)为高程异常ξ的内插值,ε为误差。设:fx,y=a0+a1x+a写写成矩阵形式有：ξ=BX+ε（4-7-3在进行选点的时候,一般有两种方法。一是基于点数选定法,即在待定点周围选择满足拟合所需的点数即可;二是动态圆半径选点法,所选的点都位于以待定点为圆心,以R为半径的圆内,其中半径R一般通过构造函数来确定。不同的采样点由于相对于待定点的距离不同,地形的变化使得它与待定点的相关程度是不同的。所以,在移动拟合时,加入起控制作用的权函数。常用的权函数有:采用简单的与距离相关的权函数:P=1r2（应用于动态圆半径选点法的距离定权函数:P=（R2-r2）r2p是参考点的权,R是圆的半径,r是待插点与参考点的平面距离。移动曲面拟合法是以每个拟合点为中心，选取周围的点参与拟合并顾忌这些点的分布及地形起伏的影响,并顾及这些点的分布及地形起伏的影响,移动拟合法采用的拟合区域相对较小,可使已知点更好地发挥控制作。其特点是计算灵活,拟合区域越小,精度越高。4.8固定边界3次样条插值法样条函数是一种连续和平滑的组合函数,该函数能在全部结点上计算和微分,这些点通过定义一个闭合区间来确定,此区间还可进一步划分为若干个子区间,对于每个子区间都可通过一组函数表示出该子区间的函数。这些函数在子区间端点处相连,形成一个连续而平滑的组合函数。2个函数相连的点称为节点,为了实现2个函数在公节点处相连,它们必须同时满足某些公共条件。当规定边界条件后,则称为固定边界样条函数。这些函数通常采用低阶多项式的形式。4.9非格网GPS散点数据考虑地形改正法有公式3-1知道高程异常主要由较为光滑部分的中长波项和局部地形起伏引起的短波项组成:ξ=ξ0+ξTC（4ξGM+ξ△g式中，ξ0为高程异常的中长波项;ξTC为高程异常的短波项。由于右了ξξξGMξξξ参考椭球面100Km参考椭球面100Km图4-9-11ookm内大地水准面波长特征基于上述理论，如果能求解分ξTC并在ξ中将扣ξTC除，则可近似认为ξ=ξ0+ξTC为一光滑的几何曲面，这样就可利用测区中已知GPS水准点的正常高和大地高确定一个ξ0的曲面函数(测区高程异常模型)，再利用该模型确定其他待求点的ξ0值。然后利用ξ=ξ04.10BP神经网络算法神经网络是模拟人脑生物神经网络机制而形成的一个并行和分布式的信息处理网络结构。它是由许多并行分布且能自主学习的神经元通过相互连接而成，可以在比较理想的精度内逼近非线性映射规律，能客观地表达输入与输出间的复杂的非线性关系。BP网络学习算法网络结构如图4-10-1所示。网络除输入层、输出层外，还有一层或多层隐含层。对于输入信号，要先向前传播到隐节点，经过激活函数后，再把隐节点的输出信息传播到输出节点，最后输出结果。节点的激活函数通常都选取标准的Sigmoid型函数。误差反传X误差反传X+——输入向量+——输入向量tt期望输出向量J期望输出向量J（导师）信号）（导师）信号）输出层第二隐层第一隐层输入层输出层第二隐层第一隐层输入层图4-10BP网络结构图BP神经网络算法原理：当输出节点的期望输出为ti时，BPyi=f(iwijxi输出节点的计算输出为:Oi=f(iTijyi节点的误差公式为：E==12i（tiBP算法的学习过程由正向和反向传播两部分组成:第一阶段(正向传播过程):给出输入信息通过输入层经隐含层逐层处理并计算每个单元的实际输出值。第二阶段(反向传播过程):若在输出层未能得到期望的输出值，则逐层递归地计算实际输出与期望输出之差值，以便根据此差值调节权值。神经网络模型采用非线性拟合方法完成区域似大地水准面的拟合,其数据结果更贴近多变的似大地水准面,因而可用较少样本值达到较高的拟合精度。而且参加拟合的已知水准点的平面位置和高程精度对预期精度存在明显影响。样本的高程应尽量采用较高等级的水准。同时BP神经网络存在着许多问题,如学习算法的收敛速度很慢,存在局部极小值,网络学习具有不稳定性等。5算法实例5.1参与计算的数据点号x（m）y（m）ξ（m）HGPSH正常13564231.7860499937.723021.074430.44209.367623565858.0800499248.000021.043231.225010.181833564029.5920499613.378021.065430.13409.068643565549.0660498813.558021.028531.110010.081553563826.3180499348.917021.056030.05909.003063566091.4010499632.434021.055430.70909.653673564312.2030500321.498021.085631.206010.120483566375.3460499179.740021.034426.42405.389693564827.1610500392.773021.082629.48708.4044103566324.2510498659.474021.071225.79004.7188113564001.7620500035.270021.078931.338010.2591203567660.2470499189.334021.017031.128010.1110213566814.6690499080.199021.027330.98909.9617223567961.3960498691.018020.997830.16009.1622233566854.8490498567.506021.007630.74109.7334263567524.9090500219.017021.051631.183010.1314283568016.2600499235.369021.013830.33209.3182表5-1-1参与计算的数据其平面分布状况为图5-1-1已知点的平面分布状况5.2平面拟合法当公共点的个数大于3时，则可列出相应的误差方程式为：vi=a0+a1写成矩阵有：V=BX-ξ（5-2-2其中V=v1v2⋮vB=1x1y11X=a0a1a2ξ=ξ1ξ2根据最小二乘原理得X=BTB-1B根据附录中的pmnh.m程序用7、4、6、20号点拟合，用9、23、26号点检验得正常高为：点号已知正常高（m）拟合正常高（m）差值（m）98.40448.40160.0028239.73349.73250.00092610.131410.12940.0320表5-2-1平面拟合所得正常高其计算过程见附录1。5.3二次曲面法拟合为了求得六个参数,至少需要有六个已知点。当测区的点多于六个的时候,利用最小二乘原理求得最佳拟合参数,即满足ε2=min。利用计算出的拟合参数和未知点的坐标值,计算未知点的高程异常vi=a0+a1写成矩阵形式有：V=BX-ξ（5-3-2其中V=(v1v2⋯vB=1x1yX=a0aξ=ξ1有最小二乘原理且根据附录中的ecqm.m程序得点号已知正常高（m）拟合正常高（m）差值（m）1110.259143.289033.02992010.111043.222733.1116219.961743.070533.10885-3-2二次多项式曲面拟合5.4加权平均值法由公式：ξA=PPi=1L根据附录中的jqpjz.m程序用前11个点进行拟合，用20，21,22，好点检验所得正常高为：点号已知正常高（m）拟合正常高（m）差值（m）2010.111010.04910.0619219.96179.91010.0516229.16229.08110.0811表5-3-1加权平均值法所得正常高5.5多面函数曲面拟合根据附录中dmhs.m程序拟合出10,11,20号点的正常高为：点号已知正常高（m）拟合正常高（m）差值（m）1110.259110.2593-0.00022010.111010.09220.0584219.96179.95560.00616结论采用拟合法求定正常高时应注意的问题：(1)采用拟合法确定地面点的正常高是建立在这样2个前提基础上的:①平差后的GPS观测值具有很高的精度,可以看作精确值;②已知高程控制点上的正常高亦可看作精确值。(2)GPS重合水准点的可靠性检验,由于生产实际情况复杂,GPS重合水准点往往等级不同,施测年代、平差精度、沉降情况不同,有的还受到某种程度的破坏,加之联测的粗差,因此必须对测区联测水准点的可靠性进行检验。(3)若使用本地区参考椭球面为基准,还应考虑到本地椭球面与WGS-84椭球面之间的差异。（4）GPS网中水准点的选择与分布。在GPS网中的点除能利用旧点已有的水准高程外,应根据需要适当进行高程联测,在GPS测量规范中有这方面的具体要求。水准重合点的分布对于拟合效果有着重要的影响。原则上要求水准重合点的分布尽可能的均匀,而且在网的边界上布设水准重合点。这样可以大大降低内插出的非重合点上的高程异常的不可靠性。另外,如可能事先根据其它的高程异常资料预测到大地水准面的形状和特征点,适于在特征点上的分布设重合点可以获得更好的拟合效果。在水准测量资料较为充分的情况下,利用GPS测量成果可以确定地面点的正常高,在顾及到以上提到问题情况下,转换后高程的精度可达厘米级,能够满足一般工程建设和大比例尺测图的需要。如有较详细的重力资料,效果更佳。充分利用GPS测量中的高程信息具有重大的经济效益,也有助于GPS定位技术的普及和推广。GPS测定正常高的研究,有利于精化大地水准面的研究,为地学领域的研究提供重要的资料。参考文献【1】刘立龙林文介。GPS测定正常高的方法研究[J]。桂林工学院学报，2001,21（3）：234～238。【2】金时华。多面函数拟合法转换GPS高程[J]。测绘与空间地理信息，2005,28（6）：44～47。【3】王旭刘文生。GPS高程拟合方法的研究[J]。测绘科学，2010,35（增）：28～29【4】王有亮张明。BP神经网络在GPS高程转换中的应用[J]。测绘与空间地理信息，2010,33（6）:124～129。【5】高伟晏磊林沂。GPS大地高转换为正常高的新方法[J]。应用基础与工程科学学报，2008，16（4）：518～523。【6】赵建虎刘经南张