一、ECEF与ENU坐标系定义
- ECEF坐标系(地心地固坐标系)
- 原点:地球质心
- X轴:指向本初子午线与赤道交点
- Y轴:在赤道平面内与X轴垂直
- Z轴:指向北极
- 数学表示: P e c e f = ( x , y , z ) P_{ecef} = (x, y, z) Pecef=(x,y,z)
- ENU坐标系(东北天坐标系)
- 原点:本地参考点
- E轴:指向东方
- N轴:指向北方
- U轴:垂直地面向上
- 数学表示: P e n u = ( e , n , u ) P_{enu} = (e, n, u) Penu=(e,n,u)
二、坐标转换数学原理
- 地理坐标转ECEF:
{ x = ( N + h ) cos ϕ cos λ y = ( N + h ) cos ϕ sin λ z = ( N ( 1 − e 2 ) + h ) sin ϕ \begin{cases} x = (N + h)\cos\phi\cos\lambda \\ y = (N + h)\cos\phi\sin\lambda \\ z = (N(1-e^2)+h)\sin\phi \end{cases} ⎩ ⎨ ⎧x=(N+h)cosϕcosλy=(N+h)cosϕsinλz=(N(1−e2)+h)sinϕ
其中:
- N = a 1 − e 2 sin 2 ϕ N = \frac{a}{\sqrt{1-e^2\sin^2\phi}} N=1−e2sin2ϕa(卯酉圈曲率半径)
- e 2 = 2 f − f 2 e^2 = 2f - f^2 e2=2f−f2(椭球偏心率平方)
- a = 6378137 m a=6378137m a=6378137m(WGS84长半轴)
- f = 1 / 298.257223563 f=1/298.257223563 f=1/298.257223563(WGS84扁率)
- ECEF转ENU:
[ e n u ] = [ − sin λ cos λ 0 − sin ϕ cos λ − sin ϕ sin λ cos ϕ cos ϕ cos λ cos ϕ sin λ sin ϕ ] [ x − x 0 y − y 0 z − z 0 ] \begin{bmatrix} e \\ n \\ u \end{bmatrix} = \begin{bmatrix} -\sin\lambda & \cos\lambda & 0 \\ -\sin\phi\cos\lambda & -\sin\phi\sin\lambda & \cos\phi \\ \cos\phi\cos\lambda & \cos\phi\sin\lambda & \sin\phi \end{bmatrix} \begin{bmatrix} x-x_0 \\ y-y_0 \\ z-z_0 \end{bmatrix} enu = −sinλ−sinϕcosλcosϕcosλcosλ−sinϕsinλcosϕsinλ0cosϕsinϕ x−x0y−y0z−z0
三、C语言实现(保存为ecef2enu.c)
#include <stdio.h>
#include <math.h>#define WGS84_A 6378137.0
#define WGS84_F 1/298.257223563typedef struct { double x, y, z; } ECEF;
typedef struct { double lat, lon, alt; } Geodetic;
typedef struct { double e, n, u; } ENU;ECEF geodetic_to_ecef(Geodetic geo) {double a = WGS84_A;double f = WGS84_F;double e2 = 2*f - f*f;double sinphi = sin(geo.lat);double cosphi = cos(geo.lat);double N = a / sqrt(1 - e2*sinphi*sinphi);ECEF ecef;ecef.x = (N + geo.alt) * cosphi * cos(geo.lon);ecef.y = (N + geo.alt) * cosphi * sin(geo.lon);ecef.z = (N*(1-e2) + geo.alt) * sinphi;return ecef;
}ENU ecef_to_enu(ECEF target, Geodetic ref_geo) {ECEF ref_ecef = geodetic_to_ecef(ref_geo);double dx = target.x - ref_ecef.x;double dy = target.y - ref_ecef.y;double dz = target.z - ref_ecef.z;double sinphi = sin(ref_geo.lat);double cosphi = cos(ref_geo.lat);double sinlam = sin(ref_geo.lon);double coslam = cos(ref_geo.lon);ENU enu;enu.e = -sinlam*dx + coslam*dy;enu.n = -sinphi*coslam*dx - sinphi*sinlam*dy + cosphi*dz;enu.u = cosphi*coslam*dx + cosphi*sinlam*dy + sinphi*dz;return enu;
}int main() {// 北京参考点(39.9042°N, 116.4074°E,海拔43m)Geodetic ref = {.lat = 39.9042 * M_PI/180,.lon = 116.4074 * M_PI/180,.alt = 43};// 生成测试数据(实际应用时可从文件读取)for(int i=0; i<5; i++) {Geodetic target_geo = {.lat = ref.lat + 0.01*i,.lon = ref.lon + 0.01*i,.alt = ref.alt + 10*i};ECEF target_ecef = geodetic_to_ecef(target_geo);ENU enu = ecef_to_enu(target_ecef, ref);printf("ENU: E=%.3fm, N=%.3fm, U=%.3fm\n", enu.e, enu.n, enu.u);}return 0;
}
编译执行:
gcc ecef2enu.c -lm -o ecef2enu && ./ecef2enu
四、Python验证代码
import numpy as np
import matplotlib.pyplot as plt
from pyproj import Transformer
from subprocess import check_output# 生成100个测试点
np.random.seed(42)
ref_lat, ref_lon = 39.9042, 116.4074
d_lats = np.random.uniform(-0.1, 0.1, 100)
d_lons = np.random.uniform(-0.1, 0.1, 100)
alts = np.random.uniform(0, 1000, 100)# 运行C程序获取结果
c_output = check_output(["./ecef2enu"]).decode().strip().split('\n')
c_enu = [list(map(float, line.split()[2::2])) for line in c_output]# 使用pyproj计算
ecef_trans = Transformer.from_crs(4326, 4978)
enu_trans = Transformer.from_crs(4326, 4467, authority="EPSG",lon_0=ref_lon, lat_0=ref_lat, h_0=0)errors = []
for i in range(100):# Python计算结果lat = ref_lat + d_lats[i]lon = ref_lon + d_lons[i]alt = alts[i]# pyproj计算ENUecef = ecef_trans.transform(lat, lon, alt)enu_py = enu_trans.transform(lat, lon, alt)# 计算误差enu_c = c_enu[i] if i <5 else [0,0,0] # 仅示例前5个errors.append(np.array(enu_py) - np.array(enu_c))# 可视化
errors = np.array(errors)
plt.figure(figsize=(12,4))plt.subplot(131)
plt.hist(errors[:,0], bins=20)
plt.title('East Error Distribution')
plt.xlabel('Meters')plt.subplot(132)
plt.hist(errors[:,1], bins=20)
plt.title('North Error Distribution')plt.subplot(133)
plt.scatter(errors[:,0], errors[:,1], alpha=0.6)
plt.xlabel('East Error')
plt.ylabel('North Error')
plt.tight_layout()
plt.savefig('enu_errors.png')
plt.show()
五、验证结果分析
- 误差直方图显示各方向误差分布
- 散点图展示平面误差的相关性
- 典型误差应小于1e-4米(数值计算误差)
该实现完整展示了从理论到实践的全流程,可通过调整测试点数量和分布进行更严格的验证。实际工程应用中需考虑坐标系的旋转顺序、大地水准面模型等更多细节。
(验证代码还有BUG,不能直接运行,但是我肝不了了,该天再调了,晚安,吗喀巴卡)