卫星大地测量高程异常编程x.docx
《卫星大地测量高程异常编程x.docx》由会员分享，可在线阅读，更多相关《卫星大地测量高程异常编程x.docx（11页珍藏版）》请在冰豆网上搜索。
卫星大地测量 编程作业
利用EGM模型计算高程异常
利用EGM模型计算高程异常
1.前言
高程异常是似大地水准面至地球椭球面的高度，在GPS水准中，求解高程异常的方法是用水准测量的方法联测GPS网中的若干GPS点的正常高，然后根据GPS点的大地高求出各公共点的高程异常，然后由公共点的平面坐标和高程异常用数值拟合的方法拟合出区域的似大地水准面，最后即可求得各点的高程异常值，进而求出各点的正常高。
这种方法在实际运用中操作不便，而且耗时耗力，精度不高。
如果用EGM96或EGM2008 求得的大地水准面差距之差能以较高的精度将GPS测定的大地高差转换为正常高差,那么,测区内既不需要联测
GPS水准点,也不需要测量正常高差,只需1～2个已知水准点就能够解决测区所有GPS点的高程转换问题。
因此，利用EGM模型计算高程异常在实际应用中是非常重要的步骤。
EGM96模型是美国NASA/GSFC和国防制图局（DMA）联合研制的360阶全球重力场模型,被公认为是同阶次模型中最好的一个。
美国国家地理空间情报局
（NGA）最新给出了EGM2008重力场模型,该新一代地球重力场模型达到2159
阶次（球谐系数的阶扩展至2190,次为2159）,空间分辨率约为5′。
高精度
、高分辨率局部或区域大地水准面不仅为大地测量、地球物理、地球动力学及海洋学等地球科学的研究和应用提供基础地球空间信息,而且也是当今构建数字地球必不可少的信息之一。
下文就利用EGM模型计算高程异常的方法进行了总结。
2.步骤
一、从ICGEM网站获取EGM96和EGM2008模型数据。
发现网站上并没有数据，于是求助同学。
只找到了EGM96数据，下面以EGM96数据为例，说明编程的思路和步骤。
二、利用Matlab对EGM96数据文件进行读取，为了方便编程，对EGM数据进行了预处理。
首先利用结构体存储每一列，然后再分别利用数组，对文件中的l,m,c,s
进行存储。
其具体的代码如下：
%c31相当于c20位系数
end
l=[];m=[];c=[];s=[];
fori=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，所以+1s（l（i）+1,m（i）+1）=zb（i）.xs;
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。
另外，还要对C20至C100减去正常重力位系数。
这一步，使C，S系数对应了l,m，方便后来的调用。
从下图中可以看到，C和S位系数在数组中表示为了361*361的下三角矩阵。
三、将大地坐标（B，L，H）转化为地心坐标（𝜑，𝜆，𝑟）。
根据上述转换模型，设计Matlab代码如下：
forL=-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;%λ
%%%%%%%%%%%%%%%%%超大循环体%%%%%%%%%%%%%%
forB=-90:
10:
90
%----------------------------坐标转换-------------------------
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结构体循环变量
注意的地方是，B在后面应转化为弧度制再带入计算。
四、计算给定点处的平均正常重力
g（1+k
sin2B）æ 2（1+a+0.00344978650684-2asin2B）H H2ö
1-e2sin2B
g= e 1
ç1- +3 ÷
ç R R2÷
è ø
k1=0.00193185265246
根据上述公式编制对应的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）））;fori=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
fori=3:
（jieshu+1）forj=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的绘图函数进行绘制。
最终结果为：
3.总结和体会
通过此次作业，我有以下收获：
通过编程计算某点的高程异常值，使得对高程异常有了更深的理解。
学会了使用Matlab进行高程异常的计算。
扩充了自己利用Matlab对数据进行读取的知识，同时也更加熟悉了Matlab中的矩阵运算和
也存在这一些不足之处，比如限于对编程语言的熟悉程度，没有开发一套界面出来，另外只用了EGM96的数据，相关的代码还有一些不足，计算的结果精度还不够等问题。
这也是我以后改进的动力。
另外，在日后的学习中，也希望自己能够加强编程语言的学习，复习回顾基础知识，多看一些专业相关的论文，弥补自身的不足。
4.附录：
高程异常计算的Matlab代码
functionreadEGM %读取EGM96文件clearall;clc;
formatlongg;
[fname,fpath]=uigetfile（'*.txt','选择EGM96文件'）;fp=fopen（strcat（fpath,fname）,'r'）;
i=1;zb=[]; %坐标矩阵为空
while~feof（fp） %文件结尾时feof（文件句柄）=0，其它为1,
%此时也可以用feof（fp）==0代替zb（i）.xl=fscanf（fp,'%d',1）;zb（i）.xm=fscanf（fp,'%d',1）;
zb（i）.xc=fscanf（fp,'%f',1）;zb（i）.xs=fscanf（fp,'%f',1）;
zb（i）.sc=fscanf（fp,'%f',1）;zb（i）.ss=fscanf（fp,'%f',1）;i=i+1;
end
fclose（fp）; %关闭句柄为fp的文件
%------------------------------------------------------------
yinzi=i-1;%总行数因子l=[];m=[];c=[];s=[];
fori=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，所以+1s（l（i）+1,m（i）+1）=zb（i）.xs;
end
c（3,1）=c（3,1）+0.484166774985E-03;%减去正常重力位系数c31相当于c20位系数
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;
%----------------------------弄好C和S系数----------------------
%----------------------------坐标转换-------------------------
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结构体循环变量
forB=0:
5:
90 %%%%%%%%%%%%%%%%%超大循环体%%%%%%%%%%%%%%forL=0:
5:
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）
升级会员 