074-RTKlib中ECEF坐标转换为大地坐标问题

最新推荐文章于 2024-07-17 17:06:57 发布

魔方的块

最新推荐文章于 2024-07-17 17:06:57 发布

阅读量2.8k

点赞数 3

CC 4.0 BY-SA版权

分类专栏： RTKLIB 文章标签： RTKLIB 坐标变换大地坐标系 ECEF

本文链接：https://blog.youkuaiyun.com/Pro2015/article/details/87940439

RTKLIB 专栏收录该内容

27 篇文章

订阅专栏

RTKlib中ECEF坐标转换为大地坐标的方法跟平时使用的不同，特此记录。
常用的大地坐标转换为ECEF坐标的公式有：

$\tag{1} \begin{cases} x = (R_n+h)cosLcos\lambda \\ y = (R_n+h)cosLsin\lambda \\ z = [R_n(1-e^2)+h]sinL \end{cases}$

各符号具体含义不再赘述。

整理（1）式中第三式可得：

$\tag{2} (R_n+h)sinL = z + R_n e^2 sinL$

为直观表现，参考下图：

为与图中字母对应，改写公式（2）为：

$\tag{3} (N+h)sin\phi = z + N e^2 sin\phi$

RTKlib程序：

extern void ecef2pos(const double *r, double *pos)
{
    double e2=FE_WGS84*(2.0-FE_WGS84),r2=dot(r,r,2),z,zk,v=RE_WGS84,sinp; 
    for (z=r[2],zk=0.0;fabs(z-zk)>=1E-4;) {
        zk=z;
        sinp = z / sqrt(r2 + z*z);
        v=RE_WGS84/sqrt(1.0-e2*sinp*sinp); 
        z=r[2]+v*e2*sinp;
    }
    pos[0]=r2>1E-12?atan(z/sqrt(r2)):(r[2]>0.0?PI/2.0:-PI/2.0);
    pos[1]=r2>1E-12?atan2(r[1],r[0]):0.0;
    pos[2]=sqrt(r2+z*z)-v;
}

由于 $sinϕ=AC/(N+h)sin\phi = AC/(N+h)$ ，所以对于sinp = z / sqrt(r2 + z*z);，得到的为小于 $ϕ\phi$ 的角度的正弦值，相应的，程序中的卯酉圈曲率半径v，即公式中的 $N$ 也是小于实际值。
循环中的最后一行z=r[2]+v*e2*sinp;为z向图中 $A C$ 的逼近。注意：z越逼近 $A C$ ，得到的纬度越准确。最终逼近 $A C$ ，从而得到较为准确的 $ϕ\phi$ ，便可进行下面的解算。