资源描述：
《基于matlabgnss水准多项式曲面拟合模型探究》由会员上传分享，免费在线阅读，更多相关内容在学术论文-天天文库。
1、基于MATLABGNSS水准多项式曲面拟合模型探究 摘要：为充分利用GNSS测得的大地高成果，改善GNSS大地高向正常高转换的精度，文中选用多项式曲面拟合法进行高程拟合研究。首先针对最小二乘估计不具备抗粗差的能力，将稳健估计引入到数据预处理中，有效地解决了数学模型失真这一问题；然后结合具体工程，对不同阶次的多项式曲面拟合数据进行对比和精度分析。全部过程利用MATLAB进行程序设计，实现了数据的自动处理。关键词：GNSS高程拟合；粗差剔除；多项式曲面拟合法；MATLAB中图分类号：[TU198+.
2、2]文献标识码：A文章编号：引言GNSS测量得到是WGS-84椭球的大地高，它没有物理意义，而中国采用的是正常高系统，它是以似大地水准面作为参考面的，因此，精确计算GNSS点的正常高，就必须作一些相应的转换。目前工程应用方面，求定地面点正常高的方法主要有GNSS水准方法。7GNSS水准是从几何解析的角度出发，在GNSS网中联测一些已知水准点，再利用这些点上的正常高和大地高求出它们的高程异常值，再据这些点上的高程异常值与坐标的关系，拟合出测区的似大地水准面，内插出其它GNSS点的高程异常，从而求出各
3、个未知点的正常高[1]。虽然GNSS数据采用了先进的采集方式，但是大量的测量数据中难免出现粗差。进行GNSS高程拟合时，多采用基于最小二乘估计来研究不同的拟合方法[2][3]，这种方法要求参与建模的数据不能含有粗差，一旦将粗差引入模型，将会使数学模型歪曲，造成参数的最小二乘估计严重失实。因此处理包含粗差的测量数据，有必要寻找具有排除粗差干扰能力的参数估计方法[4][5]。文中首先将稳健估计引入到数据的预处理中，从而将含有粗差的数据剔除，有效地解决了在精度允许的范围内数学模型失真这一问题。然后结合具
4、体工程，对工程中应用最广泛的多项式曲面拟合法进行了不同阶次的数据对比和精度分析。由于MATLAB突破了高级编程语言要编写大量循环语句的思想，它只需通过一个或几个简单的命令就可以完成矩阵或数组的运算，所以文中利用MATLAB进行程序设计，实现了数据的自动处理[6]。1稳健估计[7]所谓稳健估计，是在粗差不可避免情况下，选择适当的估计方法，使所估计参数尽可能减免粗差的影响，得出正常模式下最佳或最接近最佳的估值。7设有参数向量是未知的非随机量，为了估计，进行次观测得到向量的观测值，由极大似然估计有（1）
5、其中是随机量的密度函数。Huber于1964年提出用代替函数，使其定义广义化，于是得（2）通常残差为未知数的函数，将上式对未知数求一阶导数，并令其等于零，以求出极值点（3）设有参数向量是未知的非随机量，为了估计，进行次观测得到向量的观测值，由极大似然估计有（4）其中是随机量的密度函数。Huber于1964年提出用代替函数，使其定义广义化，于是得（5）通常残差为未知数的函数，将上式对未知数求一阶导数，并令其等于零，以求出极值点（6）法方程，如此反复，直到前后两次解的差值符合确定的精度要求为止。（4）
6、最后结果为7随着函数的选取不同，构成了权函数的多种不同的形式，通常权函数是一个在平差过程中随改正数变化的量，经过多次迭代，从而使含有粗差的奇异观测的权函数为零(或者接近于零)。而相应的残差值在很大程度上反映了其粗差值。常用的稳健估计选权迭代法有Huber法，即：（7）由Huber权函数可以看出，当所有改正数均在-和之间时，Huber估计就是经典最小二乘法估计。而当部分改正数大于时，其与改正数成反比，愈大，对应的愈小，与此相应该观测值对参数估计的影响也愈小。2多项式曲面拟合法多项式曲面拟合法，即对于
7、公共点上的高程异常与平面坐标之间，假定存在如下数学模型（8）式中，为模型待定参数。各高程控制点的己知高程异常与其拟合值之差为：（9）上式称为离差。根据最小二乘法，应在（10）的原则下，解得（8）式中的待定系数。然后再按（8）式求出测区任一点的高程异常值，从而获得点的正常高[8]。3MATLAB程序设计[9][10]3.1粗差剔除程序设计7利用二次多项式曲面拟合方法用稳健迭代法求解参数，权采用Huber法给定，程序如下：functionp=QMWJSY(XY,H,h,b1,c)X=XY(:,1);Y
8、=XY(:,2);n=length(XY);%检测XY的长度e=H-h;%求高程异常X=sym(X);Y=sym(Y);e=sym(e);%符号化%构造YY=[1，x,y,x,y,xy]YY=ones(n,1);YY=sym(YY);%符号化YY(:,2)=X;YY(:,3)=Y;YY(:,4)=X.;YY(:,5)=Y.;YY(:,6)=X.*Y;A=inv(YY’*YY)*YY’*e;V=YY*A-e;%求残差o=sqrt((V’*V)./(n-6));%求单位权中误差
