地心地固坐标系(ECEF)与站心坐标系(ENU)的转换

您所在的位置:网站首页 ecef坐标系转lla坐标系 地心地固坐标系(ECEF)与站心坐标系(ENU)的转换

地心地固坐标系(ECEF)与站心坐标系(ENU)的转换

2024-05-31 02:36| 来源: 网络整理| 查看: 265

目录

1. 概述

2. 原理

2.1. 平移

2.2. 旋转

 2.3. 总结

3. 实现

4. 参考

1. 概述

      我在《大地经纬度坐标与地心地固坐标的的转换》这篇文章中已经论述了地心坐标系的概念。我们知道,基于地心坐标系的坐标都是很大的值,这样的值是不太方便进行空间计算的,所以很多时候可以选取一个站心点,将这个很大的值变换成一个较小的值。以图形学的观点来看,地心坐标可以看作是世界坐标,站心坐标可以看作局部坐标。

       站心坐标系以一个站心点为坐标原点,当把坐标系定义为X轴指东、Y轴指北,Z轴指天,就是ENU(东北天)站心坐标系。这样,从地心地固坐标系转换成的站心坐标系,就会成为一个符合常人对地理位置认知的局部坐标系。同时,只要站心点位置选的合理(通常可选取地理表达区域的中心点),表达的地理坐标都会是很小的值,非常便于空间计算。

注意站心天向(法向量)与赤道面相交不一定会经过球心 

2. 原理

       令选取的站心点为P,其大地经纬度坐标为 (Bp,Lp,Hp) (Bp,Lp,Hp),对应的地心地固坐标系为 (Xp,Yp,Zp) (Xp,Yp,Zp)。地心地固坐标系简称为ECEF,站心坐标系简称为ENU。

2.1. 平移

       通过第一节的图可以看出,ENU要转换到ECEF,一个很明显的图形操作是平移变换,将站心移动到地心。根据站心点P在地心坐标系下的坐标 (Xp,Yp,Zp) (Xp,Yp,Zp),可以很容易推出ENU转到ECEF的平移矩阵:

 反推之,ECEF转换到ENU的平移矩阵就是T的逆矩阵:

2.2. 旋转

另外一个需要进行的图形变换是旋转变换,其旋转变换矩阵根据P点所在的经度L和纬度B确定。这个旋转变换有点难以理解,需要一定的空间想象能力,但是可以直接给出如下结论:

当从ENU转换到ECEF时,需要先旋转再平移,旋转是先绕X轴旋转 (pi/2−B) (pi/2−B),再绕Z轴旋转 (pi/2+L) (pi/2+L)当从ECEF转换到ENU时,需要先平移再旋转,旋转是先绕Z轴旋转 −(pi/2+L) −(pi/2+L),再绕X轴旋转 −(pi/2−B) −(pi/2−B)

根据我在《WebGL简易教程(五):图形变换(模型、视图、投影变换)》提到的旋转变换,绕X轴旋转矩阵为:

 绕Z轴旋转矩阵为:

从ENU转换到ECEF的旋转矩阵为:

根据三角函数公式:

sin(π/2+α)=cosαsin(π/2−α)=cosαcos(π/2+α)=−sinαcos(π/2−α)=sinα 有:

将(2)、(3)带入(1)中,则有:

 而从ECEF转换到ENU的旋转矩阵为:

 旋转矩阵是正交矩阵,根据正交矩阵的性质:正交矩阵的逆矩阵等于其转置矩阵,那么可直接得:

 2.3. 总结

将上述公式展开,可得从ENU转换到ECEF的图形变换矩阵为:

而从ECEF转换到ENU的图形变换矩阵为:

3. 实现

接下来用代码实现这个坐标转换,选取一个站心点,以这个站心点为原点,获取某个点在这个站心坐标系下的坐标:

#include #include #include using namespace std; const double epsilon = ; const double pi = ; const double d2r = pi / ; const double r2d = / pi; const double a = ; //椭球长半轴 const double f_inverse = ; //扁率倒数 const double b = a - a / f_inverse; //const double b = 6356752.314245; //椭球短半轴 const double e = sqrt(a * a - b * b) / a; void Blh2Xyz(double &x, double &y, double &z) { double L = x * d2r; double B = y * d2r; double H = z; double N = a / sqrt( - e * e * sin(B) * sin(B)); x = (N + H) * cos(B) * cos(L); y = (N + H) * cos(B) * sin(L); z = (N * ( - e * e) + H) * sin(B); } void Xyz2Blh(double &x, double &y, double &z) { double tmpX = x; double temY = y ; double temZ = z; double curB = ; double N = ; double calB = atan2(temZ, sqrt(tmpX * tmpX + temY * temY)); int counter = ; while (abs(curB - calB) * r2d > epsilon && counter < ) { curB = calB; N = a / sqrt( - e * e * sin(curB) * sin(curB)); calB = atan2(temZ + N * e * e * sin(curB), sqrt(tmpX * tmpX + temY * temY)); counter++; } x = atan2(temY, tmpX) * r2d; y = curB * r2d; z = temZ / sin(curB) - N * ( - e * e); } void TestBLH2XYZ() { //double x = 113.6; //double y = 38.8; //double z = 100; // //printf("原大地经纬度坐标:%.10lf\t%.10lf\t%.10lf\n", x, y, z); //Blh2Xyz(x, y, z); //printf("地心地固直角坐标:%.10lf\t%.10lf\t%.10lf\n", x, y, z); //Xyz2Blh(x, y, z); //printf("转回大地经纬度坐标:%.10lf\t%.10lf\t%.10lf\n", x, y, z); double x = ; double y = ; double z = ; //116.9395751953 36.7399177551 printf("地心地固直角坐标:%.10lf\t%.10lf\t%.10lf\n", x, y, z); Xyz2Blh(x, y, z); printf("转回大地经纬度坐标:%.10lf\t%.10lf\t%.10lf\n", x, y, z); } void CalEcef2Enu(Eigen::Vector3d& topocentricOrigin, Eigen::Matrix4d& resultMat) { double rzAngle = -(topocentricOrigin.x() * d2r + pi / ); Eigen::AngleAxisd rzAngleAxis(rzAngle, Eigen::Vector3d(, , )); Eigen::Matrix3d rZ = rzAngleAxis.matrix(); double rxAngle = -(pi / - topocentricOrigin.y() * d2r); Eigen::AngleAxisd rxAngleAxis(rxAngle, Eigen::Vector3d(, , )); Eigen::Matrix3d rX = rxAngleAxis.matrix(); Eigen::Matrix4d rotation; rotation.setIdentity(); rotation.block(, ) = (rX * rZ); //cout


【本文地址】


今日新闻


推荐新闻


CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3