e-frame与n-frame之间的转换

n-frame下到e-frame下速度向量的转移矩阵

相比n-frame，在e-frame中需要考虑到经纬度因素，即，不再把地球表面近似成一个平面，而是一个球面。速度向量转换关系如下

其中，可分为两次旋转，分别是绕z轴旋转 l（即精度），再绕y轴旋转 -μ-π/2（即纬度）。则可写成如下形式

 将其计算化简可得到如下关系

MATLAB实现代码如下

function R = R_l_mu(l, mu)
    R = [-cos(l) * sin(mu), -sin(l), -cos(l) * cos(mu);
         -sin(l) * sin(mu), cos(l),  -sin(l) * cos(mu);
         cos(mu),           0,       -sin(mu)        ];
end

 

从e-frame坐标得到经纬度

由于地球本身并非一个正圆，故在根据e-frame坐标计算经纬度时需要考虑曲率以及偏心率。计算示意图如下所示

 其中，r为所处位置的球心半径，r0为所处位置到地球表面投影处的球心半径，rn为地球的极轴半径，re为赤道半径，N为曲率半径，μc为投影处纬度而μ为所处位置纬度，h为所处位置到地球表面的高度。其中，N可由下式计算得出

经度可计算为

而纬度以及高度的计算方式如下

其中，e为地球的偏心率，计算方式如下

根据上述参数，如下所示为实现从e-frame坐标到经纬度以及高度的算法

 使用MATLAB代码实现如下

function [l, mu, h] = ECEF_coor_2_l_mu_h(x, y, z)
    alpha = 1e-6;

    re = 6378137;
    rp = 6356752;
    e = 0.0818;
    p = (x.^2 + y.^2).^0.5;
    mu0 = atan((z / p) * (1 - e.^2).^(-1));
    
    N0 = re.^2 / (re.^2 * cos(mu0).^2 + rp.^2 * sin(mu0).^2).^0.5;
    N = N0;
    h = p / cos(mu0) - N0;
    mu = (z / p) * (1 - e.^2 * (N0 / (N0 + h))).^(-1);
    
    while(1)
        if abs(mu-mu0) < alpha
            break;
        else
            mu0 = mu;
            N = re.^2 / (re.^2 * cos(mu0).^2 + rp * sin(mu0).^2).^0.5;
            h = p / cos(mu0) - N0;
            mu = (z / p) * (1 - e.^2 * (N0 / (N0 + h))).^(-1); 
        end
    end
    
    l = atan(y / x);
end

 从经纬度得到e-frame坐标

从经纬度得到e-frame坐标公式如下

 MATLAB代码实现如下

function [x, y, z] = l_mu_h_2_ECEF_coor(l, mu, h)
    re = 6378137;
    rp = 6356752;
    
    N = re.^2 / (re.^2 * cos(mu).^2 + rp.^2 * sin(mu).^2).^0.5;
    x = (N + h) * cos(mu) * cos(l);
    y = (N + h) * cos(mu) * sin(l);
    z = (rp.^2 / re.^2 * N + h) * sin(mu);
end