地球坐标系如何定义_经纬度转换为xy坐标

在GNSS的应用开发中，通常会使用到各种坐标系之间的转换，本文将几个经常用到的转换进行了总结，内容如下：

(1)  ECEF 与 站心坐标系（ENU/NED)之间坐标变换矩阵（或称方向余弦矩阵DCM）

(2)  ECEF坐标转LLH(大地经纬高)

(3)  LLH坐标转ECEF坐标

在介绍这些转换之前，先给出需要用到的关于地球参考椭球体的一些物理量，如下：

参考椭球体一些常量的定义：

a -- 地球参考椭球体的长轴（major axis）,在WGS84系下(a = .0)

b -- 地球参考椭球体的短轴（minor axis）,

f --  地球参考椭球体的扁率(flattening coefficient)，通常用f表示，在WGS84系下(f = 1.0/298.)

e -- 参考椭球体的偏心率（eccentricity）

  -- 大地纬度

 -- 大地经度

h -- 椭球体高度

() -- LLH坐标系下的坐标

() -- ECEF坐标系下的坐标

1、常量的值及相互关系

在WGS84坐标体系下，

a = .0

f  = 1.0/298.

2、扁率和偏心率的定义与关系

• 扁率的定义：

• 偏心率的平方 :

3、ENU与ECEF坐标变换矩阵

• ECEF to ENU

• ENU to ECEF

4、NED与ECEF坐标变换矩阵

5、ECEF 转 LLH

     按照下面的顺序计算（下式中，N是参考椭球体卯酉圈曲率半径）：

6、LLH转ECEF

7、 C Code接口函数

#include <math.h> #include <stdint.h> typedef double fp64; typedef float fp32; #define CONV_PI (3.97932) // 3.L #define EARTH_FLAT (3.e-3) // flattening coefficient of reference ellipsoid, WGS84(1.0/298.) #define EARTH_MAJA (.0) // len of semi major axis of reference ellipse const fp64 cst_earth_esqr = (EARTH_FLAT * (2.0-EARTH_FLAT)); // 1st eccentricity squared const fp64 cst_earth_omes = (1.0 - (EARTH_FLAT * (2.0-EARTH_FLAT))); // 1 minus eccentricity squared const fp64 cst_earth_asqr = (EARTH_MAJA * EARTH_MAJA); // semi major axis squared void convert_ecef2llh(fp64 pe[3], fp64 *lat, fp64 *lon, fp64 *alt) { fp64 a0,a1,a2,a3,a4,b0,b1,b2,b3; fp64 b,c0,opqk,qk,qkc,qks,f,fprm; int32_t k; fp64 yt, xt; fp64 earth_efor = cst_earth_esqr * cst_earth_esqr; // square of the 1st eccentricity squared /*-------------------------------------------*/ /* b = (x*x + y*y) / semi_major_axis_squared */ /* c = (z*z) / semi_major_axis_squared */ /*-------------------------------------------*/ b = (pe[0] * pe[0] + pe[1] * pe[1]) / cst_earth_asqr; c0 = pe[2] * pe[2] / cst_earth_asqr; /*-------------------------------------------*/ /* a0 = c * one_minus_eccentricity_sqr */ /* a1 = 2 * a0 */ /* a2 = a0 + b - first_eccentricity_to_fourth*/ /* a3 = -2.0 * first_eccentricity_to_fourth */ /* a4 = - first_eccentricity_to_fourth */ /*-------------------------------------------*/ a0 = cst_earth_omes * c0; a1 = 2.0 * a0; a2 = a0 + b - earth_efor; a3 = -2.0 * earth_efor; a4 = - earth_efor; /*-------------------------------------------*/ /* b0 = 4 * a0, b1 = 3 * a1 */ /* b2 = 2 * a2, b3 = a3 */ /*-------------------------------------------*/ b0 = 4.0 * a0; b1 = 3.0 * a1; b2 = 2.0 * a2; b3 = a3; /*-------------------------------------------*/ /* Compute First Eccentricity Squared */ /*-------------------------------------------*/ qk = cst_earth_esqr; for(k = 1; k <= 3; k++) { qks = qk * qk; qkc = qk * qks; f = a0 * qks * qks + a1 * qkc + a2 * qks + a3 * qk + a4; fprm= b0 * qkc + b1 * qks + b2 * qk + b3; qk = qk - (f / fprm); } /*-------------------------------------------*/ /* Compute latitude, longitude, altitude */ /*-------------------------------------------*/ opqk = 1.0 + qk; if ( pe[0] == 0.0 && pe[1] == 0.0 ) /* on the axis */ { /* We are sitting EXACTLY on the earth's axis.*/ /* Probably at the center or on one of the poles.*/ *lon = 0.0; /* as good as any other value */ if(pe[2] >= 0.0) { *lat = CONV_PI / 2; /* lat above north pole */ } else { *lat = -CONV_PI / 2; /* lat above south pole */ } } else { yt = opqk * pe[2]; xt = sqrt(pe[0] * pe[0] + pe[1] * pe[1]); *lat = atan2(yt, xt); *lon = atan2(pe[1], pe[0]); } *alt = (1.0 - cst_earth_omes / cst_earth_esqr * qk) * EARTH_MAJA * sqrt(b / (opqk * opqk) + c0); } void convert_llh2ecef(fp64 lat, fp64 lon, fp64 alt, fp64 *pe) { fp64 slat; fp64 clat; fp64 r; slat = sin(lat); clat = cos(lat); r = EARTH_MAJA / sqrt(1.0 - cst_earth_esqr * slat * slat); pe[0] = (r + alt) * clat * cos(lon); pe[1] = (r + alt) * clat * sin(lon); pe[2] = (r * cst_earth_omes + alt) * slat; } void calc_cen(fp64 lat, fp64 lon, fp64 cen[3][3]) { fp64 clat, clon, slat, slon; /* compute transcendentals */ slat = sin(lat); clat = cos(lat); slon = sin(lon); clon = cos(lon); #ifdef GEOD_COOR_NED // ECEF to NED cen[0][0] = -clon * slat; cen[0][1] = -slon * slat; cen[0][2] = clat; cen[1][0] = -slon; cen[1][1] = clon; cen[1][2] = 0.0; cen[2][0] = -clon * clat; cen[2][1] = -slon * clat; cen[2][2] = -slat; #else // ECEF to ENU cen[0][0] = -slon; cen[1][0] = -slat*clon; cen[2][0] = clat*clon; cen[0][1] = clon; cen[1][1] = -slat*slon; cen[2][1] = clat*slon; cen[0][2] = 0; cen[1][2] = clat; cen[2][2] = slat; #endif }

8、a simple demo

假设有4个基站（A0，A1，A2，A3），希望得到基站分别站心坐标系下的相对位置坐标、在LLH坐标下的大地坐标、以及它们之间的互转。

站心坐标系的规定：a、选择ENU作为站心系   b、以A0的位置作为站心系的原点

 下面直接给出它们的之间的相互转换的C代码函数接口：

 /** * the coordination of the reference base station */ typedef struct { fp64 pe[3]; fp64 lat,lon,alt; fp64 cen[3][3]; }rbs_coor_t; static rbs_coor_t _rbs_coor; /** * @name : init_rbs_coordination * @brief: initialize the reference coordination * @param(in): fp64 lat , latitude of RBS, unit: rad * @param(in): fp64 lon , longitude of RBS, unit: rad * @param(in): fp64 alt , height of RBS, unit: m * @return: none */ void init_rbs_coordination(fp64 lat, fp64 lon, fp64 alt) { _rbs_coor.lat = lat; _rbs_coor.lon = lon; _rbs_coor.alt = alt; calc_cen(lat, lon, _rbs_coor.cen); convert_llh2ecef(lat, lon, alt, _rbs_coor.pe); } /** * @name : convert_xyz_to_llh * @brief: * @param(in) : const fp64 *xyz, * @param(out): fp64 *lat , pointer of latitude of the object, unit: rad * @param(out): fp64 *lon , pointer of longitude of the object, unit: rad * @param(out): fp64 *alt , pointer of height of the object, unit: m * @return: none */ void convert_xyz_to_llh(const fp64 *xyz, fp64 *lat, fp64 *lon, fp64 *alt) { fp64 prx[3]; fp64 dpe[3]; if (lat == 0 || lon == 0 || alt == 0) return; // dpe = Cne * dpn dpe[0] = _rbs_coor.cen[0][0] * xyz[0] + _rbs_coor.cen[1][0] * xyz[1] + _rbs_coor.cen[2][0] * xyz[2]; dpe[1] = _rbs_coor.cen[0][1] * xyz[0] + _rbs_coor.cen[1][1] * xyz[1] + _rbs_coor.cen[2][1] * xyz[2]; dpe[2] = _rbs_coor.cen[0][2] * xyz[0] + _rbs_coor.cen[1][2] * xyz[1] + _rbs_coor.cen[2][2] * xyz[2]; // prx[0] = _rbs_coor.pe[0] + dpe[0]; prx[1] = _rbs_coor.pe[1] + dpe[1]; prx[2] = _rbs_coor.pe[2] + dpe[2]; // convert_ecef2llh(prx, lat, lon, alt); } /** * @name : convert_llh_to_xyz * @brief: * @param(in): fp64 lat , latitude of the object, unit: rad * @param(in): fp64 lon , longitude of the object, unit: rad * @param(in): fp64 alt , height of the object, unit: m * @param(out) : fp64 *xyz, * @return: none */ void convert_llh_to_xyz(fp64 lat, fp64 lon, fp64 alt, fp64 *xyz) { fp64 dpe[3]; fp64 pe[3]; convert_llh2ecef(lat,lon,alt,pe); dpe[0] = pe[0] - _rbs_coor.pe[0]; dpe[1] = pe[1] - _rbs_coor.pe[1]; dpe[2] = pe[2] - _rbs_coor.pe[2]; xyz[0] = _rbs_coor.cen[0][0] * dpe[0] + _rbs_coor.cen[0][1] * dpe[1] + _rbs_coor.cen[0][2] * dpe[2]; xyz[1] = _rbs_coor.cen[1][0] * dpe[0] + _rbs_coor.cen[1][1] * dpe[1] + _rbs_coor.cen[1][2] * dpe[2]; xyz[2] = _rbs_coor.cen[2][0] * dpe[0] + _rbs_coor.cen[2][1] * dpe[1] + _rbs_coor.cen[2][2] * dpe[2]; } #if 1 #include <stdio.h> int main(char argc, char argv[]) { fp64 lat0 = 39. * CONV_PI / 180.; fp64 lon0 = 116. * CONV_PI / 180.; fp64 alt0 = 31.2; fp64 A1[3] = {-6,0,0}; fp64 A2[3] = {0,-6,0}; fp64 A3[3] = {-6,-6,0}; fp64 lat,lon,alt; init_rbs_coordination(lat0, lon0, alt0); convert_xyz_to_llh(A1, &lat, &lon, &alt); printf("A1(lat,lon,alt) = {%10.6lf,%10.6lf,%10.3lf}\n",lat*180./CONV_PI,lon*180./CONV_PI,alt); convert_xyz_to_llh(A2, &lat, &lon, &alt); printf("A2(lat,lon,alt) = {%10.6lf,%10.6lf,%10.3lf}\n",lat*180./CONV_PI,lon*180./CONV_PI,alt); convert_xyz_to_llh(A3, &lat, &lon, &alt); printf("A3(lat,lon,alt) = {%10.6lf,%10.6lf,%10.3lf}\n",lat*180./CONV_PI,lon*180./CONV_PI,alt); return 0; } #endif

上面程序运行的结果为：

A1(lat,lon,alt) = { 39.,116.,    31.200}
A2(lat,lon,alt) = { 39.,116.,    31.200}
A3(lat,lon,alt) = { 39.,116.,    31.200}