卫星大地测量高程异常编程xWord文档下载推荐.docx
《卫星大地测量高程异常编程xWord文档下载推荐.docx》由会员分享,可在线阅读,更多相关《卫星大地测量高程异常编程xWord文档下载推荐.docx(11页珍藏版)》请在冰豆网上搜索。
只找到了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:
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;
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);
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
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,'
zb(i).xc=fscanf(fp,'
%f'
zb(i).xs=fscanf(fp,'
zb(i).sc=fscanf(fp,'
zb(i).ss=fscanf(fp,'
i=i+1;
fclose(fp);
%关闭句柄为fp的文件
%------------------------------------------------------------
yinzi=i-1;
%总行数因子l=[];
%将l下标变成矩阵
%将m下标变成矩阵
%矩阵下标不能为0,所以+1s(l(i)+1,m(i)+1)=zb(i).xs;
%减去正常重力位系数c31相当于c20位系数
c(5,1)=c(5,1)-0.790303733511E-06;
c(9,1)=c(9,1)-0.34605248394E-11;
c(11,1)=c(11,1)+0.265002225747E-14;
%----------------------------弄好C和S系数----------------------
GM=3.986004415E+14;
PI=3.1415926535897932384626433;
e2=0.00669437999013;
ge=9.7803253359;
%mmm是world结构体循环变量
forB=0:
5:
90 %%%%%%%%%%%%%%%%%超大循环体%%%%%%%%%%%%%%forL=0:
%取出度AA=(B-PP)*100;
%取出分QQQ=fix(AA);
%取出秒A=(PP+QQQ/60.00+degs/36.0)*pi/180.0;
X=(N+H)*cos(A)*cos(L);
fai=atan(Z/(sqrt(X*X+Y*Y)