GPS水准高程拟合报告
GPS水准高程拟合报告
实验目的：
1掌握GPS水准高程拟合的基本原理，了解高精度GPS水准的研究意义；
2能够利用Matlab编程实现几何内插法拟合GPS水准高程；
实验内容：
利用Matlab编程实现几何内插法拟合GPS水准高程，并作内插结果分析
实验原理：
1大地水准面，参考椭球面，正高，大地高之间的几何关系
A 正高的定义是：由地面点沿通过该点的铅垂线到大地水准面的距离。一般用符号Hg表示。
B 大地高的定义是：由地面点沿通过该点的椭球面法线到椭球面的距离。也称为椭球高，一般用符号H表示。大地高是一个纯几何量，不具有物理意义。同一个点，在不同的基准下，具有不同的大地高。利用GPS，可以测定地面点的WGS-84中的大地高。
C 大地水准面差距：大地水准面到椭球面的距离，称为大地水准面差距，记为hg （或N）。
如上图可以看出大地高和正高之间的关系：H=Hg+ hg
2几何内插法原理
几何内插法是通过一些既进行了GPS观测又具有水准资料的点上的大地水准面差距，采用平面或曲面拟合，配置三次样条等内插方法，得到其他点上的大地水准面差距从而反算这些点上的正高。
3二多项式拟合
N=a0+a1*dB+a2*dL+a3*dB2+a4*dL2+a5*dL*dB 公式一
式中dB=B-B0;dL=L-L0;B0=1/n∑B;L0=1/n∑L,n为GPS观测点的数量。
利用其中一些具有水准观测资料的公共点上的的大地高和正高可以计算出这些点的大地水准面差距。利用这些公共点的观测资料求得公式一的参数，再利用求得的公式进行其他点的大地水准面差距内插，和正高的拟合；
实验步骤：
1输入已知点的GPS观测值和相应的正常高构成矩阵B,L,H,h,分别是纬度矩阵，经度矩阵，大地高矩阵，正高矩阵；
2计算dB=B-B0;dL=L-L0;B0=1/n∑B;L0=1/n∑L，构成矩阵矩阵dB,dL和大地水准面差距矩阵N=H-h;
3将以上计算得到的矩阵代入公式一经过间接平差求得相应的参数a i,这样就能构成一个确定的多项式二；
4输入未知待求点的GPS观测值构成矩阵BB,LL,HH，计算相应的dBB,dLL；
5将dBB,dLL矩阵代入多项式二，解算出对应点的大地水准面差距NN矩阵；
6反算各点的正高h=H-NN;
7对计算得到的正高，大地水准面差距做对比分析；
实验分析：
1本实验中可以选择两种差值公式算法
（1）N=a0+a1*dB+a2*dL+a3*dB2+a4*dL2+a5*dL*dB
（2）N=a0+a1*B+a2*L+a3*B2+a4*L2+a5*L*B
采用公式（1）的插值结果如下：
Δh(dB)散点图
注：Δh(dB)是插值点的水准资料与插值结果的差值
采用公式（2）的插值结果如下：
Δh(dB)散点图
由以上散点图可以看出以下特点：
A 采用d
B dL作为变量值求得的插值函数所插值的结果与真是值差距普遍集中在-0.4到-0.6之间，显然这种拟合结果的精度并不是很高，另外我们不得不怀疑这种拟合结果存在某种系统误差因为其散点图具有一定的偏向性，同时作为公共点已知数据的十个点的差值结果亦如想象的那样理想，这证明我们的差值函数参数的解算没用问题。
B采用 B L作为变量值求得的插值函数所插值的结果与真是值差距普遍集中在-0.4到0.4之间，显然这种拟合结果的精度也不是很高，但较第一种差值有了一定的提高，这种插值结果的散点图分布较为随机不太可能存在系统误差，同时作为公共点已知数据的十个点的差值结果亦如想象的那样理想，这证明我们的差值函数参数的解算没用问题。
总结以上拟合结果尽管其精度不尽相同但是两种方法的精度显然不能与传统的水准测量相媲美，完全不能满足实际工程水准的需要。
2 下面采用拉格朗日n次插值对原数据在进行处理
拉格朗日三次差值残差分布
拉格朗日四次差值残差分布
拉格朗日五次差值残差分布
拉格朗日六次差值残差分布
注：1由于五次和六次插值的效果较差，五次插值仅选取了拟合效果相对较好的前39项作图，六次插值仅选取前25项；
2拟合多项式参数求解并未进行平差解算；
分析拉格朗日插值的几幅残差图发现在四次插值时其残差分布与Δh(dB)散点图比较接近，其他的几种高次插值效果都比较差，而且随着插值次数的提高，拟合效果也来越差其残差成指数增加，由此可见多项式插值对GPS水准拟合并不是一种好的拟合方式；对此我们可以尝试用空间插值，例如克里金插值法。
实验总结：
与常规的水准测量而言GPS水准具有效率高费用低等特点，可以在大范围内进行高程测量，对于大型工程而言传统的大地水准将耗费大量的人力物力财力，且观测周期非常长，如果GPS可行那么将为大范围水准测量开辟崭新的空间时代，同时高精度的GPS水准也将改变传统大地测量在高程基准面确定上的难题，为测绘作业提供高精度的大地水准面。
但是就目前而言GPS水准高程的精度问题是其进行工程应用的难题，其原因主要是大地水准面差距的精度低，目前我们还难以获得较高精度的全球大地水准面，解决GPS水准的主要方法主要是几何插值法，就本次实验可以看出多项式插值的精度并不能满足实际需要，但是综合目前各方面的实际情况，GPS水准可以满足三等水准的要求。对于一些特殊的工程而言要获取高精度的GPS水准采用的处理方式是，在局部区域内建立若干GPS长期观测基站利用这些基站的观测值和IGS精密星历解算该区域的大地水准面，以获得精度较高的大地水准面差距从而能获得精度较高的高程，这方面的例子主要是港珠澳大桥项目，不过由于成本太高耗时太长一般的项目难以企及，不具有工程实用性。
附Matlab程序代码
%原始数据
B=[4160726,4154896,4159270,4148735,4157177,4153413,4157949,4151556,41 55573,4157027]'
L=[469370,473772,473622,474714,474356,476134,
476233,477565,476918,477814]'
H=[43.085,47.677,37.002,56.859,43.904,78.045, 41.697,51.519,33.57, 32.715]' h=[15.395,19.758,9.137,28.900,15.993,50.071,13.721,23.492,5.574,4.685]'
%定义dB=B-B0
dB=zeros(10,1);
dL=zeros(10,1);
B0=0;
L0=0;
for i=1:10
B0=B0+B(i,1);
L0=L0+L(i,1);
end
B0=B0/10.0;
L0=L0/10.0;
for i=1:10
dB(i,1)=B(i,1)-B0;
dL(i,1)=L(i,1)-L0;
end
%V=XA-N;求相应数据
X=zeros(10,6);
%初始化X，赋值
for i=1:10
X(i,1)=1.0;
X(i,2)=dB(i,1);
X(i,3)=dL(i,1);
X(i,4)=dB(i,1)*dB(i,1);
X(i,5)=dL(i,1)*dL(i,1);
X(i,6)=dL(i,1)*dB(i,1);
end
%求N
N=zeros(10,1);
for i=1:10
N(i,1)=H(i,1)-h(i,1);
end
%间接平差求A；X'PV=0 V=XA-N P=E；
A=zeros(6,1);
s=X'*X;
t=inv(s);
u=X'*N;
A=t*u;
%读入要计算的原始数据
BB=[4149569,4152803,4155546,4159369,4157612,4147678,4155194,4159612,4 151539,4157232,4150583,4158508,4160739,4156273,4154803,4157419,4163806,4 163116,4146497,4153030,4143589,4148280,4150278,4156103,4146715,4149128,4 154168,4155429,4144712,4148524,
4150861,4154013,4146430,4145664,4148161,4143343,4145375,4144379,4142039, 4139017]'
LL=[479940,479819,479795,479113,479861,482166,482154,482603,482291,483 485,485202,485355,485636,486046,487445,488532,488447,490930,486936,487566 ,488117,487969,488166,490412,490503,490823,491023,491935,492338,494523,49 4172,493828,495124,498357,497419,500731,500943,504035,504085,504544]' HH=[35.852,35.099,33.292,32.411,32.857,40.737,34.026,33.772,35.689,33.418, 65.076,39.009,32.066,105.537,128.872,32.045,37.03,34.798,49.149,84.576,
44.983,58.323,86.432,36.14,59.754,88.67,100.893,76.64,49.964,39.018,52.246,71.3 9,33.804,32.367,32.36,64.272,58.95,31.911,37.269,37.935]'
BB0=0;
LL0=0;
for i=1:40
BB0=BB0+BB(i,1);
LL0=LL0+LL(i,1);
end
BB0=BB0/40.0;
LL0=LL0/40.0;
dBB=zeros(40,1);
dLL=zeros(40,1);
XX=zeros(40,6);
%计算各点大地水准面差距NN
for i=1:40
dBB(i,1)=BB(i,1)-BB0;
dLL(i,1)=LL(i,1)-LL0;
XX(i,1)=1.0;
XX(i,2)=dBB(i,1);
XX(i,3)=dLL(i,1);
XX(i,4)=dBB(i,1)*dBB(i,1);
XX(i,5)=dLL(i,1)*dLL(i,1);
XX(i,6)=dLL(i,1)*dBB(i,1);
end
%MM是10个已知点的大地水准面差距的插值矩阵，o是10个已知值的正高矩阵
NN=zeros(40,1);
NN=XX*A;
MM=X*A;
hh=HH-NN;
o=H-MM;