Matlab||EGM2008模型计算GOCE沿轨重力梯度及全球重力梯度分布

地球重力场模型是地球重力位的数学表达式，通常用球谐函数的级数形式来表示。地球外部任意一点P地心求坐标（r,θ,λ），其重力位V的球谐表达式为：

根据上式可知，只需知道地球外部一点的球坐标，以及相应重力场模型球谐系数，即可解算该位置的重力位。

重力矢量和重力梯度分别是重力位对球坐标的一阶、二阶导数，还涉及球坐标与直角坐标的转换问题，计算较为复杂。球坐标与局部指北坐标的转换关系如下：

EGM2008是近来由美国国家地理空间情报局（NGA）发布的全球超高阶地球重力场模型，模型球谐系数阶次扩展至2190次，相当于模型的空间分辨率约为5′(约9 km)。该模型采用了GRACE卫星跟踪数据、卫星测高数据和地面重力数据等数据。

读取EGM2008模型球谐系数

fid = fopen(filename);
hasErrors = true;%读取头文件
while 1line = fgetl(fid);if ~ischar(line), break, endif isempty(line), continue, endkeyword = textscan(line,'%s',1);%获取最大阶数if(strcmp(keyword{1}, 'max_degree'))cells = textscan(line,'%s%d',1);maxDegree = cells{2};end%获取参考地球半径if(strcmp(keyword{1}, 'radius'))cells = textscan(line,'%s%f',1);R = cells{2};end%获取地心引力参数GMif(strcmp(keyword{1}, 'earth_gravity_constant'))cells = textscan(line,'%s%f',1);GM = cells{2};endif(strcmp(keyword{1}, 'errors'))cells = textscan(line,'%s%s',1);if(strcmp(cells{2}, 'no'))hasErrors = false;endendif(strcmp(keyword{1}, 'end_of_head'))breakend
end% 声明球谐系数阵
cnm = zeros(maxDegree+1, maxDegree+1);
snm = zeros(maxDegree+1, maxDegree+1);% 读取球谐系数
while 1line = fgetl(fid);if ~ischar(line), break, endif isempty(line), continue, endif(hasErrors)cells = textscan(line,'%s%d%d%f%f%f%f');elsecells = textscan(line,'%s%d%d%f%f');endcnm(cells{2}+1, cells{3}+1) = cells{4};snm(cells{2}+1, cells{3}+1) = cells{5};
endfclose(fid);

计算球坐标系下的重力梯度分量

因EGM2008是超高阶重力场模型，为了减轻计算量，取到180阶。data采用上篇GOCE卫星重力梯度观测数据，计算沿轨梯度值。

maxDegree = size(cnm,1)-1;  %0阶开始
Vq=zeros(2678400,9);        %声明一个球坐标系下重力梯度矩阵for  n=1:2678400  %matlab的矩阵是从1开始的，不是0开始。Plm = legendreFunctions(90-data(n,2), 182); %地心余纬for l=0:180   %l阶，取到1000阶for m=0:l          %m次%组合法计算勒让德级数的一阶导、二阶导%计算VrrVq(n,1) =Vq(n,1)+ GM/R*(l+1)*(l+2)/power(data(n,4),2) *power(R/data(n,4),l+1)*(cnm(l+1,m+1)*cosd(m*data(n,3)) + snm(l+1,m+1)*sind(m*data(n,3)))* Plm(l+1,m+1);%计算VλλVq(n,2)=Vq(n,2)+GM/R*(-1)*power(R/data(n,4),l+1)*(cnm(l+1,m+1)*cosd(m*data(n,3)) + snm(l+1,m+1)*sind(m*data(n,3)))* Plm(l+1,m+1)*m*m;%计算Vθθif m==0Vq(n,3)=Vq(n,3)+GM/R*power(R/data(n,4),l+1)*(cnm(l+1,m+1)*cosd(m*data(n,3)) + snm(l+1,m+1)*sind(m*data(n,3)))*(-l*(l+1)/2*Plm(l+1,1)+sqrt(l*(l-1)*(l+1)*(l+2)/8)*Plm(l+1,3));elseif m==1Vq(n,3)=Vq(n,3)+GM/R*power(R/data(n,4),l+1)*(cnm(l+1,m+1)*cosd(m*data(n,3)) + snm(l+1,m+1)*sind(m*data(n,3)))*((2*l*(l+1)+(l-1)*(l+2))/-4*Plm(l+1,2)+sqrt((l-2)*(l-1)*(l+2)*(l+3))/4*Plm(l+1,4));elseVq(n,3)=Vq(n,3)+GM/R*power(R/data(n,4),l+1)*(cnm(l+1,m+1)*cosd(m*data(n,3)) + snm(l+1,m+1)*sind(m*data(n,3)))*(sqrt((l+m)*(l-m+1)*(l+m-1)*(l-m+2))/4*Plm(l+1,m-1)-((l+m)*(l-m+1)+(l-m)*(l+m+1))/4*Plm(l+1,m+1)+sqrt((l-m)*(l+m+1)*(l-m-1)*(l+m+2))/4*Plm(l+1,m+3));           end%计算VrλVq(n,4)=Vq(n,4)+GM/R*(l+1)/-data(n,4)*power(R/data(n,4),l+1)*(snm(l+1,m+1)*cosd(m*data(n,3))-cnm(l+1,m+1)*sind(m*data(n,3)))*m*Plm(l+1,m+1);%计算Vrθif m==0Vq(n,5)=Vq(n,5)+GM/R*(l+1)/-data(n,4)*power(R/data(n,4),l+1)*(cnm(l+1,m+1)*cosd(m*data(n,3)) + snm(l+1,m+1)*sind(m*data(n,3)))*sqrt(1*(l+1)/2)*Plm(l+1,2);elseif m==1Vq(n,5)=Vq(n,5)+GM/R*(l+1)/-data(n,4)*power(R/data(n,4),l+1)*(cnm(l+1,m+1)*cosd(m*data(n,3)) + snm(l+1,m+1)*sind(m*data(n,3)))*(sqrt(2*l*(l+1))/-2*Plm(l+1,1)+sqrt((l-1)*(l+2))/2*Plm(l+1,3));elseVq(n,5)=Vq(n,5)+GM/R*(l+1)/-data(n,4)*power(R/data(n,4),l+1)*(cnm(l+1,m+1)*cosd(m*data(n,3)) + snm(l+1,m+1)*sind(m*data(n,3)))*(sqrt((l+m)*(l-m+1))/-2*Plm(l+1,m)+sqrt((l-m)*(l+m+1))/2*Plm(l+1,m+2));end%计算Vλθif m==0Vq(n,6)=Vq(n,5)+GM/R*power(R/data(n,4),l+1)*(snm(l+1,m+1)*cosd(m*data(n,3))-cnm(l+1,m+1)*sind(m*data(n,3)))*m*sqrt(1*(l+1)/2)*Plm(l+1,2);elseif m==1Vq(n,6)=Vq(n,5)+GM/R*power(R/data(n,4),l+1)*(snm(l+1,m+1)*cosd(m*data(n,3))-cnm(l+1,m+1)*sind(m*data(n,3)))*m*(sqrt(2*l*(l+1))/-2*Plm(l+1,1)+sqrt((l-1)*(l+2))/2*Plm(l+1,3));elseVq(n,6)=Vq(n,5)+GM/R*power(R/data(n,4),l+1)*(snm(l+1,m+1)*cosd(m*data(n,3))-cnm(l+1,m+1)*sind(m*data(n,3)))*m*(sqrt((l+m)*(l-m+1))/-2*Plm(l+1,m)+sqrt((l-m)*(l+m+1))/2*Plm(l+1,m+2));end%%  为了后续坐标转换，先算好三个一阶导%计算VrVq(n,7)=Vq(n,7)+GM/R*(l+1)/-data(n,4)*power(R/data(n,4),l+1)*(cnm(l+1,m+1)*cosd(m*data(n,3)) + snm(l+1,m+1)*sind(m*data(n,3)))* Plm(l+1,m+1);%计算Vθif m==0Vq(n,8)=Vq(n,8)+GM/R*power(R/data(n,4),l+1)*(cnm(l+1,m+1)*cosd(m*data(n,3)) + snm(l+1,m+1)*sind(m*data(n,3)))*sqrt(1*(l+1)/2)*Plm(l+1,2);elseif m==1Vq(n,8)=Vq(n,8)+GM/R*power(R/data(n,4),l+1)*(cnm(l+1,m+1)*cosd(m*data(n,3)) + snm(l+1,m+1)*sind(m*data(n,3)))*(sqrt(2*l*(l+1))/-2*Plm(l+1,1)+sqrt((l-1)*(l+2))/2*Plm(l+1,3));elseVq(n,8)=Vq(n,8)+GM/R*power(R/data(n,4),l+1)*(cnm(l+1,m+1)*cosd(m*data(n,3)) + snm(l+1,m+1)*sind(m*data(n,3)))*(sqrt((l+m)*(l-m+1))/-2*Plm(l+1,m)+sqrt((l-m)*(l+m+1))/2*Plm(l+1,m+2));end%计算VλVq(n,9)=Vq(n,9)+GM/R*power(R/data(n,4),l+1)*(snm(l+1,m+1)*cosd(m*data(n,3))-cnm(l+1,m+1)*sind(m*data(n,3)))*m*Plm(l+1,m+1);endend

坐标转换
将球坐标系下的重力梯度值转换成局部指北坐标系。

Vcal=zeros(size(Vq,1),6);
for n=1:size(Vq,1)%VxxVcal(n,1)=Vq(n,7)/data(n,4)+Vq(n,3)/power(data(n,4),2);%VyyVcal(n,2)=Vq(n,7)/data(n,4)+cotd(90-data(n,2))*Vq(n,8)/power(data(n,4),2)+Vq(n,2)/(power(data(n,4),2)*power(sind(90-data(n,2)),2));%VzzVcal(n,3)=Vq(n,1);%VxyVcal(n,4)=Vq(n,6)/(power(data(n,4),2)*sind(90-data(n,2)))-cosd(90-data(n,2))*Vq(n,9)/(power(data(n,4),2)*power(sind(90-data(n,2)),2));%VxzVcal(n,5)=Vq(n,8)/power(data(n,4),2)-Vq(n,5)/data(n,4);%VyzVcal(n,6)=Vq(n,9)/(power(data(n,4),2)*sind(90-data(n,2)))-Vq(n,4)/data(n,4)/sind(90-data(n,2));
end

将EGM2008模型计算的重力梯度值与GOCE卫星观测值求差，得到以下残差图：
可见，GOCE卫星观测值还需要进一步滤波去噪，改正其他扰动因子。

附：

笔者尝试绘制全球2.55°分辨率的地心向径重力梯度分量Vrr分布图，因缺少全球地表地心距，选取地心距6.410^6m参考球面的梯度分布图。


Vglobe=zeros(71,73);
for i=1:71for j=1:73Plm = legendreFunctions(90-(2.5*i-90), 182); %地心余纬for l=0:180   %l阶，取到180阶for m=0:l    %m次  %计算VrrVglobe(i,j) =Vglobe(i,j)+ GM/R*(l+1)*(l+2)/power(6.4*10^6,2) *power(R/(6.4*10^6),l+1)*(cnm(l+1,m+1)*cosd(m*(5*j-180)) + snm(l+1,m+1)*sind(m*(5*j-180)))* Plm(l+1,m+1);endendend
end%画图---等值线
[X,Y] = meshgrid(-180:5:180,-87.5:2.5:87.5);
contourf(-180:5:180,-87.5:2.5:87.5,Vglobe,100,'LineStyle','none');
hold on;% 海岸线绘制：coast.dat为海岸线坐标文件
load('coast.dat');
plot(coast(:,1),coast(:,2),'k');colorbar;
shading flat;
title('重力梯度Vrr(r=6.4*10^6)')
xlabel('Longitude','FontSize',12)
ylabel('Latitude','Fontsize',12)
axis tight

分布图如下：（正确性有待考察）

Matlab||EGM2008模型计算GOCE沿轨重力梯度及全球重力梯度分布相关推荐

matlab iri模型,iri-model IRI模型计算电离层延迟的matlab 算法 265万源代码下载- www.pudn.com...
文件名称: iri-model下载收藏√ [ 5 4 3 2 1 ] 开发工具: matlab 文件大小: 4182 KB 上传时间: 2013-05-29 下载次数: 85 详细说明: ...
零维水温模型计算——基于matlab建立计算模型
零维水温模型计算--基于matlab建立计算模型零维水温模型基本方程计算参数准备 matlab模型的建立四阶龙格-库塔法的matlab实现水面的净流通量φ~s~的计算函数湖水与岩土的换热φ ...
e-006 matlab,基于MATLAB进行潮流计算
基于MATLAB 进行潮流计算学生:王仕龙 2011148213 指导老师:李咸善摘要:电力系统潮流计算方法有两类,即手算潮流和计算机潮流计算.手算潮流主要借助于形成简化的等值电路来实现,这种方法 ...
【Python】SEBS模型计算蒸散发
教程照片及其他详细信息请关注微信公众号:夫也的笔记公众号内容包含:GEE.ArcGIS.ENVI.MATLAB.Python和R语言教程和实际案例分享 SEBS模型介绍近年来出现了许多遥感反演蒸散 ...
matlab三相短路电流计算程序_基于MATLAB的短路电流计算程序编制.pdf
基于MATLAB的短路电流计算程序编制维普资讯 2008年第4期 <贵州电力技术> (总第 l06期) 基于 MATLAB的短路电流计算程序编制武汉大学电气工程学院周冬旭向俊杰陈 ...
matlab与科学计算王沫然,MATLAB与科学计算（第3版）王沫然著电子工业出版社 9787121180521...
商品描述: 基本信息书名:MATLAB与科学计算(第3版) 定价:49.80元作者:王沫然编著出版社:电子工业出版社出版日期:2012-10-01 ISBN:9787121180521 字数 ...
基于MATLAB梁模型振型动画程序设计
基于MATLAB梁模型振型动画程序设计梁模型振型动画 %该模型为一端固定梁模型振型动画 clear; close all; clc % 系统参数 E = 1e7; A = 1.5; rho = 2. ...
Matlab软件使用讲解(5)，Matlab数学工程计算2023a中文版下载安装教程
Matlab是一款被广泛应用于科学.技术和工程等领域的数学软件,在数据分析.算法研究.信号处理等方面具有很高的实用性.在本文中,我将通过举例的方式向大家介绍Matlab软件的常用功能和使用技巧,并结合 ...
matlab klobuchar模型,区域似大地水准面精化模型算法的优选
区域似大地水准面精化模型算法的优选第20卷第1期 2011年2月 ENGINEERING 测绘工程 Vol.20l.1OFSURVEYINGANDMAPPING Feb.,2011 区域似大地 ...

Matlab||EGM2008模型计算GOCE沿轨重力梯度及全球重力梯度分布

Matlab||EGM2008模型计算GOCE沿轨重力梯度及全球重力梯度分布相关推荐

最新文章

热门文章