高阶次勒让德计算（跨阶次递推公式）

最新推荐文章于 2025-10-29 20:53:40 发布

原创最新推荐文章于 2025-10-29 20:53:40 发布 · 3.7k 阅读

28 ·

CC 4.0 BY-SA版权

文章标签：

#大地测量 #学习交流

学习笔记专栏收录该内容

2 篇文章

订阅专栏

本文介绍如何在卫星大地测量中运用跨阶次递推勒让德公式计算高程异常。通过从ICGEM获取数据，进行坐标转换和重力计算，最终利用Matlab实现计算过程。

跨阶次递推勒让德公式

今天卫星大地测量课的一个问题
问题：已知某一点的大地坐标BLH，求该点的高程异常。

（1）从ICGEM网站获取模型数据，并提取l,m,c,s数据。
（2）在c位系数中减去正常重力位系数
（3）将给定点的大地坐标转化为地心坐标
（4）计算给定点处的平均正常重力
（5）Legendre递推公式的计算(标准前向列递推公式)（跨阶次递推公式）
（6）计算给定点处高程异常

matlab代码：

clc;clear;
% B = input('输入大地纬度B: ');
% L = input('输入大地经度L：');
% H = input('输入大地高： ');
H = 0;
GM = 0.3986004415E+15;
R = 0.6378136300E+07;
Degree = 360;
e2=0.00669437999013;
ge=9.7803253359;
k1=0.00193185265246;
Alpha=1/298.257223563;
Beta=0.00344978650684;
Pi=3.1415926535897932384626433;
% B = B*Pi/180;
% L = L*Pi/180;

%读取LMCS数据（读取前删除数据的前20行介绍）
path = 'C:\Users\Administrator\Desktop\EGM2008.gfc';
Data = importdata(path);
LMCS = Data.data;
C1 = {-4.841667749835220e-04,7.903037335100000e-07,-1.687249611513883e-09,3.460524683954118e-12,-2.650022257475667e-15};

%减去正常重力位系数
z=1;
for i = 2:54
    if(rem(LMCS(i,1),2)==0&&LMCS(i,2)==0)
        LMCS(i,3) = LMCS(i,3) - double(C1{z});        
        z=z+1;
    end
end
%定义高程异常矩阵
Anomaly = zeros(180,360);
%计算全球高程异常（B为-89到90，L为0到359）
for B1 = -89:90
    for L1 = 0:359
        B = B1*Pi/180;
        L = L1*Pi/180;                
        %坐标转换（BLH到XYZ到Phi,Lamda,r）
N = R/sqrt(1-e2*sin(B)^2);
X = (N+H)*cos(B)*cos(L);
Y = (N+H)*cos(B)*sin(L);
Z = (N*(1-e2)+H)*sin(B);
Phi = atan(Z/sqrt(X^2+Y^2));
Lam = L;
r = sqrt(X^2+Y^2+Z^2);
%j计算正常重力
y1 = ge*(1+k1*sin(B)^2)/sqrt(1-e2*sin(B)^2);
y2 = 1 - 2*H*(1+Alpha+Beta-2*Alpha*sin(B)^2)/R+3*H^2/R^2;
y = y1*y2;
P = zeros(Degree+1,Degree+1);

% %跨阶次正规化勒让德系数
P(1,1) = 1;
P(2,1) = sin(Phi)*3^0.5;
P(2,2) = sqrt(3*(1-sin(Phi)^2));

for j = 1:2
   for i = 3:Degree+1
       l=i-1;m=j-1;
       a = sqrt((4*l^2-1)/(l^2-m^2));
       b = sqrt((2*l+1)/(2*l-3))*sqrt(((l-1)^2-m^2)/(l^2-m^2));
       P(i,j) = a*sin(Phi)*P(i-1,j)-b*P(i-2,j);      
   end
end

for j = 3:Degree
    for i = j:j+1
        l=i-1;m=j-1;
      alpha = sqrt((2*l+1)*(l-m)*(l-m-1)/(2*l-3)/(l+m)/(l+m-1));
      if(m==2)
          beta = sqrt(2*(2*l+1)*(l+m-2)*(l+m-3)/(2*l-3)/(l+m)/(l+m-1));
          gama = sqrt(2*(l-m+1)*(l-m+2)/(l+m)/(l+m-1));
      else
          beta = sqrt((2*l+1)*(l+m-2)*(l+m-3)/(2*l-3)/(l+m)/(l+m-1));
          gama = sqrt((l-m+1)*(l-m+2)/(l+m)/(l+m-1));
      end
      P(i,j) = beta*P(i-2,j-2)-gama*P(i,j-2);
    end
    if j+2<Degree+1
        for i = j+2:Degree+1
        l = i-1; m = j-1;
        alpha = sqrt((2*l+1)*(l-m)*(l-m-1)/(2*l-3)/(l+m)/(l+m-1));
        if m ==2
            beta = sqrt(2*(2*l+1)*(l+m-2)*(l+m-3)/(2*l-3)/(l+m)/(l+m-1));
            gama = sqrt(2*(l-m+1)*(l-m+2)/(l+m)/(l+m-1));
        else
            beta = sqrt((2*l+1)*(l+m-2)*(l+m-3)/(2*l-3)/(l+m)/(l+m-1));
            gama = sqrt((l-m+1)*(l-m+2)/(l+m)/(l+m-1));
        end
        P(i,j) = alpha*P(i-2,j)+beta*P(i-2,j-2)-gama*P(i,j-2);
        end
    end
end
l = Degree;m = Degree;
beta = sqrt((2*l+1)*(l+m-2)*(l+m-3)/(2*l-3)/(l+m)/(l+m-1));
gama = sqrt((l-m+1)*(l-m+2)/(l+m)/(l+m-1));
P(l+1,m+1) = beta*P(i-2,j-2)-gama*P(i,j-2);

%求高程异常
z = 2;S = 0;
for i = 2 : Degree
    for j = 0 : i
        S = (R/r)^i*P(i+1,j+1)*(LMCS(z,3)*cos(Lam*j)+LMCS(z,4)*sin(Lam*j))+S;
        z = z+1;
    end
end
Anomaly(B1+90,L1+1) = GM/(y*r)*S;

    end
end
[m,n] = size(Anomaly);
%输出全球高程异常数据为number.txt(当前目录下)
fid = fopen('number.txt','wt');
for i = 1:m
    for j = 1:n
        if j == n
            fprintf(fid,'%g\n',Anomaly(i,j));
        else
            fprintf(fid,'%g\t',Anomaly(i,j));
        end
    end
end
fclose(fid);
%导入数据画图
Ano = importdata('number.txt');
load geoid;
worldmap('World');
geoshow(Ano,geoidrefvec,'DisplayType','texturemap');