ECEF与ENU坐标系定义及C语言实现

一、ECEF与ENU坐标系定义

ECEF坐标系（地心地固坐标系）

原点：地球质心
X轴：指向本初子午线与赤道交点
Y轴：在赤道平面内与X轴垂直
Z轴：指向北极
数学表示： $P_{ecef} = (x, y, z)$

ENU坐标系（东北天坐标系）

原点：本地参考点
E轴：指向东方
N轴：指向北方
U轴：垂直地面向上
数学表示： $P_{enu} = (e, n, u)$

二、坐标转换数学原理

地理坐标转ECEF：
$\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}$
其中：

$\frac{a}{\sqrt{1-e^2\sin^2\phi}}$ （卯酉圈曲率半径）
$e^2 = 2f - f^2$ （椭球偏心率平方）
$a = 6378137 m$ （WGS84长半轴）
$f = 1/298.257223563$ （WGS84扁率）

ECEF转ENU：
$\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}$

三、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()