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

南望山小雷达

已于 2023-03-30 20:24:54 修改

阅读量5.7k

点赞数 11

文章标签：算法 matlab

于 2021-02-25 21:31:44 首次发布

本文链接：https://blog.youkuaiyun.com/FQ_QIU/article/details/114106853

版权

地球重力场模型是地球重力位的数学表达式，通常用球谐函数的级数形式来表示。地球外部任意一点P地心求坐标（r,θ,λ），其重力位V的球谐表达式为：
在这里插入图片描述
根据上式可知，只需知道地球外部一点的球坐标，以及相应重力场模型球谐系数，即可解算该位置的重力位。

重力矢量和重力梯度分别是重力位对球坐标的一阶、二阶导数，还涉及球坐标与直角坐标的转换问题，计算较为复杂。球坐标与局部指北坐标的转换关系如下：
在这里插入图片描述
EGM2008是近来由美国国家地理空间情报局（NGA）发布的全球超高阶地球重力场模型，模型球谐系数阶次扩展至2190次，相当于模型的空间分辨率约为5′(约9 km)。该模型采用了GRACE卫星跟踪数据、卫星测高数据和地面重力数据等数据。

读取EGM2008模型球谐系数

fid = fopen(filename);
hasErrors = true;

%读取头文件
while 1
  line = fgetl(fid);
  if ~ischar(line), break, end
  if isempty(line), continue, end
  keyword = 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

%获取地心引力参数GM
  if(strcmp(keyword{1}, 'earth_gravity_constant'))
    cells = textscan(line,'%s%f',1);
    GM = cells{2};
  end

  if(strcmp(keyword{1}, 'errors'))
    cells = textscan(line,'%s%s',1);
    if(strcmp(cells{2}, 'no'))
      hasErrors = false;
    end
  end

  if(strcmp(keyword{1}, 'end_of_head'))
    break
  end
end

% 声明球谐系数阵
cnm = zeros(maxDegree+1, maxDegree+1);
snm = zeros(maxDegree+1, maxDegree+1);

% 读取球谐系数
while 1
  line = fgetl(fid);
  if ~ischar(line), break, end
  if isempty(line), continue, end
  if(hasErrors)
    cells = textscan(line,'%s%d%d%f%f%f%f');
  else
    cells = textscan(line,'%s%d%d%f%f');
  end
  cnm(cells{2}+1, cells{3}+1) = cells{4};
  snm(cells{2}+1, cells{3}+1) = cells{5};
end

fclose(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次
         %组合法计算勒让德级数的一阶导、二阶导
      %计算Vrr
      Vq(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==0
           Vq(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==1
           Vq(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));
      else
           Vq(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==0
          Vq(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==1
          Vq(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));
      else
          Vq(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==0
          Vq(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==1
          Vq(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));
      else
          Vq(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
      
      %%  为了后续坐标转换，先算好三个一阶导
      %计算Vr
      Vq(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==0
          Vq(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==1
          Vq(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));
      else
          Vq(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);
     end
    end

% legendreFunctions函数

function Pnm = legendreFunctions(theta1, maxDegree)

theta=theta1*pi/180; %角度转弧度
global legendreFactor1;
global legendreFactor2;

if (size(legendreFactor1,1) < maxDegree+1)
  legendreFactor1 = zeros(maxDegree+1, maxDegree+1);
  legendreFactor2 = zeros(maxDegree+1, maxDegree+1);
  legendreFactor1(2,2) = sqrt(3);
  for n=2:maxDegree
    legendreFactor1(n+1,n+1) = sqrt((2*n+1)/(2*n));
  end
  for m=0:maxDegree-1
    for n=m+1:maxDegree
      f=(2*n+1)/((n+m)*(n-m));
      legendreFactor1(n+1,m+1) = sqrt(f*(2*n-1));
      legendreFactor2(n+1,m+1) = sqrt(f*(n-m-1)*(n+m-1)/(2*n-3));
    end
  end
end

cosTheta = cos(theta);
sinTheta = sin(theta);

Pnm = zeros(maxDegree+1, maxDegree+1);
Pnm(1,1)=1;
for n=1:maxDegree
  Pnm(n+1,n+1) = legendreFactor1(n+1,n+1)*sinTheta*Pnm(n,n);
end
for m=0:maxDegree-1
  Pnm(m+2,m+1) = legendreFactor1(m+2,m+1)*cosTheta*Pnm(m+1,m+1);
end
for m=0:maxDegree-1
  for n=m+2:maxDegree
    Pnm(n+1,m+1) = legendreFactor1(n+1,m+1)*cosTheta*Pnm(n,m+1)-legendreFactor2(n+1,m+1)*Pnm(n-1,m+1);
  end
end

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

Vcal=zeros(size(Vq,1),6);
for n=1:size(Vq,1)
    %Vxx
    Vcal(n,1)=Vq(n,7)/data(n,4)+Vq(n,3)/power(data(n,4),2);
    %Vyy
    Vcal(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));
    %Vzz
    Vcal(n,3)=Vq(n,1);
    %Vxy
    Vcal(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));
    %Vxz
    Vcal(n,5)=Vq(n,8)/power(data(n,4),2)-Vq(n,5)/data(n,4);
    %Vyz
    Vcal(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:71
    for j=1:73
        
         Plm = legendreFunctions(90-(2.5*i-90), 182); %地心余纬
         for l=0:180   %l阶，取到180阶
           for m=0:l    %m次  
           %计算Vrr
           Vglobe(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);
     
           end
        end

    end
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