经纬高(LLA)坐标系地心地固(ECEF)坐标系与东北天(ENU)坐标系转换

        前段时间在做水下机器人项目,添加了RTK,读取到的数据为经纬高坐标系中度分形式的经纬度信息,无法直接用于定位,还需要进行坐标系的转换。后来在学习Cartographer时处理GPS数据也有提到这方面的知识,于是决定汇总到一起进行学习,下文将对相关内容进行介绍。

一.LLA坐标系与ECEF坐标系之间的转换

        其中lon为经度信息,lat为纬度信息,alt为高度信息,f为极扁率,N为基准椭球体的曲率半径,将传入的信息先转换为ECEF坐标系下的X Y Z

        当已知ECEF坐标系下点为(x,y,z)时,需转换为LLA坐标系下的(lon,lat,alt):

       最初lon是未知的,可以假设为0,经过计策迭代之后就能收敛。

二.ECEF坐标系转ENU坐标系

        ECEF坐标系原点为地球的质心,x轴延伸通过本初子午线和赤道的交点。 z轴延伸通过的北极(即与地球旋转轴重合)。 y轴完成右手坐标系,穿过赤道和90度经度。

         坐标原点P_{0}(x_{0},y_{0},z_{0}),目标点P(x,y,z),以P_{0}为坐标原点的东北天坐标系坐标(e,n,u),LLA坐标系下的坐标LLA_{0} (lon_{0},lat_{0},alt_{0}),其中S为正交矩阵:

        而从ENU坐标系转换到ECEF坐标系为:

       以下Cartographer代码中,为计算第一帧GPS数据指向ECEF坐标系下原点的坐标变换:

cartographer::transform::Rigid3d ComputeLocalFrameFromLatLong(
    const double latitude, const double longitude) {
  const Eigen::Vector3d translation = LatLongAltToEcef(latitude, longitude, 0.);
  const Eigen::Quaterniond rotation =
      Eigen::AngleAxisd(cartographer::common::DegToRad(latitude - 90.),
                        Eigen::Vector3d::UnitY()) *
      Eigen::AngleAxisd(cartographer::common::DegToRad(-longitude),
                        Eigen::Vector3d::UnitZ());
  return cartographer::transform::Rigid3d(rotation * -translation, rotation);
}

        首先调用的是LatLongAltToEcef函数,从名字不难猜到这个函数是用来将得到的经纬高坐标系下的经纬度坐标转换成ECEF坐标系下的坐标:

Eigen::Vector3d LatLongAltToEcef(const double latitude, const double longitude,
                                 const double altitude) {
  // note: 地固坐标系(Earth-Fixed Coordinate System)也称地球坐标系, 
  // 是固定在地球上与地球一起旋转的坐标系.
  // 如果忽略地球潮汐和板块运动, 地面上点的坐标值在地固坐标系中是固定不变的.

  // https://en.wikipedia.org/wiki/Geographic_coordinate_conversion#From_geodetic_to_ECEF_coordinates
  
  constexpr double a = 6378137.;  // semi-major axis, equator to center.
  constexpr double f = 1. / 298.257223563;
  constexpr double b = a * (1. - f);  // semi-minor axis, pole to center.
  constexpr double a_squared = a * a;
  constexpr double b_squared = b * b;
  constexpr double e_squared = (a_squared - b_squared) / a_squared;
  const double sin_phi = std::sin(cartographer::common::DegToRad(latitude));
  const double cos_phi = std::cos(cartographer::common::DegToRad(latitude));
  const double sin_lambda = std::sin(cartographer::common::DegToRad(longitude));
  const double cos_lambda = std::cos(cartographer::common::DegToRad(longitude));
  const double N = a / std::sqrt(1 - e_squared * sin_phi * sin_phi);
  const double x = (N + altitude) * cos_phi * cos_lambda;
  const double y = (N + altitude) * cos_phi * sin_lambda;
  const double z = (b_squared / a_squared * N + altitude) * sin_phi;

  return Eigen::Vector3d(x, y, z);
}

        回到ComputeLocalFrameFromLatLong函数,rotation为ECEF坐标系下的坐标到 local 坐标系的旋转矩阵,translation 表示ECEF坐标系原到 (latitude, longitude) 位置的平移,该旋转矩阵表示先绕Y轴调整经度,然后再绕Z轴调整纬度。

  • 6
    点赞
  • 55
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
转换 ENU 坐标系ECEF 坐标系: ```c++ #include <cmath> // 定义常量 const double pi = 3.14159265358979323846; const double a = 6378137; // WGS84 椭球体长半轴 const double f = 1.0 / 298.257223563; // WGS84 椭球体扁率 const double b = a * (1 - f); // WGS84 椭球体短半轴 const double e2 = 1 - pow((b / a), 2); // WGS84 椭球体第一偏心率的平方 // 转换函数 void enu2ecef(double& x, double& y, double& z, const double& lat0, const double& lon0, const double& h0, const double& east, const double& north, const double& up) { double sinp = sin(lat0 * pi / 180); double cosp = cos(lat0 * pi / 180); double sinl = sin(lon0 * pi / 180); double cosl = cos(lon0 * pi / 180); double dx = -sinl * east - cosl * north; double dy = cosl * east - sinl * north; double dz = up; x = -sinl * dx + cosl * dy + cosp * dz + a * cosp * e2 * sinp; y = -cosl * dx - sinl * dy + cosp * dz + a * cosp * e2 * sinp; z = + sinp * dz + b * sinp * sinp * sinp + h0; } ``` 转换 ECEF 坐标系经纬坐标系: ```c++ #include <cmath> // 定义常量 const double pi = 3.14159265358979323846; const double a = 6378137; // WGS84 椭球体长半轴 const double f = 1.0 / 298.257223563; // WGS84 椭球体扁率 const double b = a * (1 - f); // WGS84 椭球体短半轴 const double e2 = 1 - pow((b / a), 2); // WGS84 椭球体第一偏心率的平方 // 转换函数 void ecef2lla(const double& x, const double& y, const double& z, double& lat, double& lon, double& h) { double r = sqrt(x * x + y * y); double E = sqrt(a * a - b * b); double F = 54 * b * b * z * z; double G = r * r + (1 - e2) * z * z - e2 * E * E; double c = (e2 * e2 * F * r * r) / (G * G * G); double s = pow(1 + c + sqrt(c * c + 2 * c), 1.0 / 3.0); double P = F / (3 * pow((s + 1 / s + 1), 2) * G * G); double Q = sqrt(1 + 2 * e2 * e2 * P); double r0 = -(P * e2 * r) / (1 + Q) + sqrt(0.5 * a * a * (1 + 1.0 / Q) - P * (1 - e2) * z * z / (Q * (1 + Q)) - 0.5 * P * r * r); double U = sqrt(pow((r - e2 * r0), 2) + z * z); double V = sqrt(pow((r - e2 * r0), 2) + (1 - e2) * z * z); double zo = (b * b * z) / (a * V); h = U * (1 - b * b / (a * V)); lat = atan((z + e2 * zo) / r); lon = atan2(y, x); } ``` 转换经纬坐标系ECEF 坐标系: ```c++ #include <cmath> // 定义常量 const double pi = 3.14159265358979323846; const double a = 6378137; // WGS84 椭球体长半轴 const double f = 1.0 / 298.257223563; // WGS84 椭球体扁率 const double b = a * (1 - f); // WGS84 椭球体短半轴 const double e2 = 1 - pow((b / a), 2); // WGS84 椭球体第一偏心率的平方 // 转换函数 void lla2ecef(const double& lat, const double& lon, const double& h, double& x, double& y, double& z) { double sinp = sin(lat); double cosp = cos(lat); double sinl = sin(lon); double cosl = cos(lon); double N = a / sqrt(1 - e2 * sinp * sinp); x = (N + h) * cosp * cosl; y = (N + h) * cosp * sinl; z = (N * (1 - e2) + h) * sinp; } ``` 转换 ENU 坐标系经纬坐标系: ```c++ #include <cmath> // 定义常量 const double pi = 3.14159265358979323846; const double a = 6378137; // WGS84 椭球体长半轴 const double f = 1.0 / 298.257223563; // WGS84 椭球体扁率 const double b = a * (1 - f); // WGS84 椭球体短半轴 const double e2 = 1 - pow((b / a), 2); // WGS84 椭球体第一偏心率的平方 // 转换函数 void enu2lla(const double& lat0, const double& lon0, const double& h0, const double& east, const double& north, const double& up, double& lat, double& lon, double& h) { double x, y, z; enu2ecef(x, y, z, lat0, lon0, h0, east, north, up); ecef2lla(x, y, z, lat, lon, h); } ``` 转换经纬坐标系ENU 坐标系: ```c++ #include <cmath> // 定义常量 const double pi = 3.14159265358979323846; const double a = 6378137; // WGS84 椭球体长半轴 const double f = 1.0 / 298.257223563; // WGS84 椭球体扁率 const double b = a * (1 - f); // WGS84 椭球体短半轴 const double e2 = 1 - pow((b / a), 2); // WGS84 椭球体第一偏心率的平方 // 转换函数 void lla2enu(const double& lat0, const double& lon0, const double& h0, const double& lat, const double& lon, const double& h, double& east, double& north, double& up) { double x0, y0, z0; double x, y, z; lla2ecef(lat0, lon0, h0, x0, y0, z0); lla2ecef(lat, lon, h, x, y, z); double dx = x - x0; double dy = y - y0; double dz = z - z0; double sinp = sin(lat0 * pi / 180); double cosp = cos(lat0 * pi / 180); double sinl = sin(lon0 * pi / 180); double cosl = cos(lon0 * pi / 180); east = -sinl * dx + cosl * dy; north = -sinp * cosl * dx - sinp * sinl * dy + cosp * dz; up = cosp * cosl * dx + cosp * sinl * dy + sinp * dz; } ```
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值