利用EGM96和EGM2008模型，分别计算计算高程异常、重力异常和垂线偏差。

1.前言

高程异常是似大地水准面至地球椭球面的高度，在GPS水准中，求解高程异常的方法是用水准测量的方法联测GPS网中的若干GPS点的正常高，然后根据GPS点的大地高求出各公共点的高程异常，然后由公共点的平面坐标和高程异常用数值拟合的方法拟合出区域的似大地水准面，最后即可求得各点的高程异常值，进而求出各点的正常高。这种方法在实际运用中操作不便，而且耗时耗力，精度不高。如果用EGM96或EGM2008求得的大地水准面差距之差能以较高的精度将 GPS测定的大地高差转换为正常高差, 那么测区内既不需要联测GPS水准点 ,也不需要测量正常高差,只需1～2个已知水准点就能够解决测区所有GPS点的高程转换问题。因此，利用EGM模型计算高程异常在实际应用中是非常重要的步骤。

EGM96模型是美国 NASA/GSFC和国防制图局 (DMA)联合研制的 360 阶全球重力场模型 , 被公认为是同阶次模型中最好的一个。美国国家地理空间情报局 (NGA)最新给出了EGM2008重力场模型 , 该新一代地球重力场模型达到2159阶次(球谐系数的阶扩展至2190,次为2159),空间分辨率约为5′。高精度、高分辨率局部或区域大地水准面不仅为大地测量、地球物理、地球动力学及海洋学等地球科学的研究和应用提供基础地球空间信息,而且也是当今构建数字地球必不可少的信息之一 。下文就利用EGM2008模型计算高程异常的方法进行了总结。

2.步骤

一、从ICGEM网站获取EGM96和EGM2008模型数据。发现网站上并没有数据，于是求助同学。只找到了EGM96数据，下面以EGM96数据为例，说明编程的思路和步骤。

二、利用Matlab对EGM96数据文件进行读取，为了方便编程，对EGM数据进行了预处理。

首先利用结构体存储每一列，然后再分别利用数组，对文件中的l,m,c,s进行存储。其具体的代码如下：

l=[];m=[];c=[];s=[];

for i=1:yinzi

l(i)=zb(i).xl;%将l下标变成矩阵

m(i)=zb(i).xm;%将m下标变成矩阵

c(l(i)+1,m(i)+1)=zb(i).xc;%矩阵下标不能为0，所以+1

s(l(i)+1,m(i)+1)=zb(i).xs;

end                               %c31相当于c20位系数

c(3,1)=c(3,1)+0.484166774985E-03; %减去正常重力位系数

c(5,1)=c(5,1)-0.790303733511E-06;

c(7,1)=c(7,1)+0.168724961151E-08;

c(9,1)=c(9,1)-0.34605248394E-11;

c(11,1)=c(11,1)+0.265002225747E-14;

值得注意的地方是：Matlab中数组的下标不能为0，故在上述第5行代码中加了1。另外，还要对C2 0至C10 0减去正常重力位系数。这一步，使C，S系数对应了l,m，方便后来的调用。从下图中可以看到，C和S位系数在数组中表示为了361*361的下三角矩阵

三、将大地坐标（B，L，H）转化为地心坐标（φ，λ，r）。

根据上述转换模型，设计Matlab代码如下：

%----------------------------坐标转换-------------------------

R=6378136.460;GM=3.986004415E+14;PI=3.1415926535897932384626433;

e2=0.00669437999013;ge=9.7803253359;aerfa=1/298.257223563;%扁率

k1=0.00193185265246;B=-90;L=-180;H=0;World=[];mmm=1; 

%mmm是world结构体循环变量

for B=-90:10:90    %%%%%%%%%%%%%%%%% 超大循环体 %%%%%%%%%%%%%%

    for L=-180:10:180

    H=0;

PP=fix(B);  AA=(B-PP)*100; QQQ=fix(AA);degs=AA-QQQ; %化弧度

A=(PP+QQQ/60.00+degs/36.0)*pi/180.0;

N=R/(sqrt(1-e2*sin(A)*sin(A)));

X=(N+H)*cos(A)*cos(L);

Y=(N+H)*cos(A)*sin(L);

Z=(N*(1-e2)+H)*sin(A);

fai=atan(Z/(sqrt(X*X+Y*Y)));%地心坐标

lanmat=L;%λ

注意的地方是，B在后面应转化为弧度制再带入计算。

四、计算给定点处的平均正常重力

根据上述公式编制对应的Matlab代码。

五、Legendre 递推公式的计算

根据前面推算的C和S系数阵以及l,m参数，利用上述公式可以计算P值并综合到矩阵中。P矩阵也是361*361的下三角矩阵。其具体的代码已经运行结果如下：

%-------------------------Legendre 递推公式的计算-------------------

p=[];

jieshu=360;%EGM96为360阶，若是EGM2008，则修改成2190.

p(1,1)=1;

p(2,1)=sqrt(3)*sin(fai);

p(2,2)=sqrt(3*(1-sin(fai)*sin(fai)));

for i=3:(jieshu+1)

    p(i,i-1)=sqrt(2*(i-1)+1)*sin(fai)*p(i-1,i-1); %此处i-1对应公式中的L

    p(i,i)=sqrt((2*(i-1)+1)/(2*(i-1)))*cos(fai)*p(i-1,i-1);

end

for i=3:(jieshu+1)

    for j=1:i-2

        %式子太长，避免括号错位，分开计算

    ddd=sqrt((2*(i-1)*(i-1))/((i-1)*(i-1)-(j-1)*(j-1)))*sin(fai)*p(i-1,j);

dddd=sqrt(((2*(i-1)+1)/(2*(i-1)-3))*((i-1-1)*(i-1-1)-j*j)/((i-1)*(i-1)-j*j))*p(i-2,j);

p(i,j)=ddd-dddd;

    end

end

P矩阵的运行结果如下

六、计算N值

由上述模型计算N值。并写入到文本中，利用Matlab的绘图函数进行绘制。最终结果为：

附录：Matlab代码

https://download.csdn.net/download/qq_38437948/19400765?spm=1001.2014.3001.5503