卫星大地测量高程异常编程xWord格式文档下载.docx
《卫星大地测量高程异常编程xWord格式文档下载.docx》由会员分享，可在线阅读，更多相关《卫星大地测量高程异常编程xWord格式文档下载.docx（11页珍藏版）》请在冰点文库上搜索。
%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）;
%地心坐标lanmat=L;
r=sqrt（X*X+Y*Y+Z*Z）;
%----------------------------坐标转换-----------------------------
%----------------------------%计算给定点处的平均正常重力------------
gama=ge*（1+k1*sin（A）*sin（A））/（sqrt（1-e2*sin（A）*sin（A）））*（1-（2*（1+aerfa+0.00344978650684-2*aerfa*sin（A）*sin（A））*H）/R+3*H*H/R*R）;
%EGM96为360阶，若是EGM2008，则修改成2190.p（1,1）=1;
p（2,2）=sqrt（3*（1-sin（fai）*sin（fai）））;
%此处i-1对应公式中
的L
p（i,i）=sqrt（（2*（i-1）+1）/（2*（i-1）））*cos（fai）*p（i-1,i-1）;
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）;
%----------------------------高程异常计算---------------------------
sum=0;
sum2=0;
i
sum=（（R/r）^i）*p（i,j）*（c（i,j）*cos（（j-1）*lanmat）+s（i,j）*sin（（j-1）*lanmat））;
sum2=sum2+sum;
GCYC=GM/r/gama*sum2;
World（mmm）.x=B;
World（mmm）.y=L;
World（mmm）.z=GCYC;
mmm=mmm+1;
%%%%%%%%%%%超大循环%%%%%%%%%%%%%%%%%
%--------------------------写入到文件过程----------------------------
-----
[fname2,fpath2]=uiputfile（'
请选择坐标计算后保存的数据文件'
fn=fopen（strcat（fpath2,fname2）,'
w'
%写入表头fprintf（fn,'
%s\r\n'
B L 高程异常'
fclose（fn）;
fn=fopen（strcat（fpath2,fname2）,'
a'
%a是指在文件末尾添加数据fori=1:
（mmm-1）
fprintf（fn,'
%8d'
World（i）.x）;
%写入坐标X数据（共10位小数5
位）
World（i）.y）;
%写入坐标Y数据（共16位小数5
%16.7f\r\n'
World（i）.z）;
%写入坐标H数据（共16位
小数5位）
%第一行写完，注意在最后\r\n换行
endfclose（fn）end
升级会员 