卫星大地测量高程异常编程xWord文档下载推荐.docx

1、只找到了 EGM96 数据，下面以 EGM96 数据为例，说明编程的思路和步骤。二、利用 Matlab 对 EGM96 数据文件进行读取，为了方便编程，对 EGM 数据进行了预处理。首先利用结构体存储每一列，然后再分别利用数组，对文件中的 l,m,c,s进行存储。其具体的代码如下：%c31相当于c20位系数endl=;m=;c=;s=;for i=1:yinzil（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;c（3,1

2、）=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 位系数在数组中表示

3、为了 361*361 的下三角矩阵。三、将大地坐标（B，L，H）转化为地心坐标（𝜑，𝜆，𝑟）。根据上述转换模型，设计 Matlab 代码如下：for L=-180:10:180H=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）

4、;fai=atan（Z/（sqrt（X*X+Y*Y）;%地心坐标lanmat=L;% 超大循环体 %for B=-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

5、（1+ ksin2 B）2 （1+a+ 0.00344978650684 - 2asin2 B）HH 2 1- e2 sin2 Bg =e11-+ 3RR2 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;

6、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对应公式中的Lp（i,i）=sqrt（2*（i-1）+1）/（2*（i-1）*cos（fai）*p（i-1,i-1）;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）

7、;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; endP 矩阵的运行结果如下：六、计算N值由上述模型计算 N 值。并写入到文本中，利用 Matlab 的绘图函数进行绘制。最终结果为：3. 总结和体会通过此次作业，我有以下收获：通过编程计算某点的高程异常值，使得对高程异常有了更深的理解。学会了使用 Matlab 进行高程异常的计算。扩充了自己利用 Matlab 对数据进行读取的知识，同时也更加熟悉了 Matlab 中的矩阵运算和也存在这一些不足之处

8、，比如限于对编程语言的熟悉程度，没有开发一套界面出来，另外只用了 EGM96 的数据，相关的代码还有一些不足，计算的结果精度还不够等问题。这也是我以后改进的动力。另外，在日后的学习中，也希望自己能够加强编程语言的学习，复习回顾基础知识，多看一些专业相关的论文，弥补自身的不足。4. 附录：高程异常计算的 Matlab 代码function readEGM%读取 EGM96 文件clear all;clc;format long g;fname,fpath=uigetfile（*.txt,选择 EGM96 文件）; fp=fopen（strcat（fpath,fname）,ri=1;zb=;% 坐

9、标矩阵为空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,%fzb（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，所以+1 s（l（i）+1,m（i）+1）=zb

10、（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 结构体循环变量for B=0:5:90% 超大循环体 % for L=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）