第1章:ECEF与ENU坐标转换——从地心到站心的坐标系变换
1.1 为什么需要ECEF和ENU?
做导航定位的人,天天跟坐标系打交道。我个人习惯把坐标系分成两类:一类是“全局坐标系”,比如ECEF(地心地固坐标系);另一类是“局部坐标系”,比如ENU(东北天坐标系)。
ECEF的原点在地球质心,X轴指向本初子午线与赤道的交点,Z轴指向北极。说白了,它是一个跟地球一起转的坐标系。你想想看,卫星的位置、地面站的坐标,通常都用ECEF表示。
ENU就不一样了。它的原点在你脚下,东向为X,北向为Y,天向为Z。我在项目中遇到过这样的场景:无人机起飞后,飞控系统需要知道它相对于起飞点的位置。这时候用ENU就特别直观——往东飞了多少米,往北飞了多少米,高度变化了多少。
所以,ECEF和ENU之间的转换,是导航定位里绕不开的基础操作。
1.2 转换原理:旋转加平移
ECEF到ENU的转换,本质上是一个“旋转+平移”的过程。嗯,这里要注意:ENU的原点(也就是站心)在ECEF中有一个固定的坐标,我们记为(x₀, y₀, z₀),对应的经纬度是(λ₀, φ₀)。
转换分两步走:
- 平移:把ECEF坐标减去站心坐标,得到相对于站心的向量。
- 旋转:把这个向量从ECEF坐标系旋转到ENU坐标系。
旋转矩阵长这样:
R = | -sin(λ₀) cos(λ₀) 0 |
| -sin(φ₀)cos(λ₀) -sin(φ₀)sin(λ₀) cos(φ₀) |
| cos(φ₀)cos(λ₀) cos(φ₀)sin(λ₀) sin(φ₀) |
这个矩阵怎么来的?我简单解释一下:第一行对应东向,第二行对应北向,第三行对应天向。每一行都是ECEF三个轴在对应方向上的投影分量。
核心公式:
ENU坐标 = R × (ECEF坐标 - 站心ECEF坐标)
1.3 代码实现:Python实战
理论讲完了,咱们直接上代码。我习惯用numpy做矩阵运算,简洁又高效。
import numpy as np
def ecef_to_enu(x, y, z, lat0, lon0, h0):
"""
ECEF转ENU
:param x, y, z: 目标点的ECEF坐标(米)
:param lat0, lon0, h0: 站心的纬度(度)、经度(度)、海拔(米)
:return: 东向、北向、天向坐标(米)
"""
# 先把经纬度转成弧度
lat0_rad = np.radians(lat0)
lon0_rad = np.radians(lon0)
# 计算站心的ECEF坐标
a = 6378137.0 # WGS84长半轴
e2 = 0.00669437999014 # 第一偏心率平方
N = a / np.sqrt(1 - e2 * np.sin(lat0_rad)**2)
x0 = (N + h0) * np.cos(lat0_rad) * np.cos(lon0_rad)
y0 = (N + h0) * np.cos(lat0_rad) * np.sin(lon0_rad)
z0 = (N * (1 - e2) + h0) * np.sin(lat0_rad)
# 平移
dx = x - x0
dy = y - y0
dz = z - z0
# 旋转矩阵
R = np.array([
[-np.sin(lon0_rad), np.cos(lon0_rad), 0],
[-np.sin(lat0_rad)*np.cos(lon0_rad), -np.sin(lat0_rad)*np.sin(lon0_rad), np.cos(lat0_rad)],
[ np.cos(lat0_rad)*np.cos(lon0_rad), np.cos(lat0_rad)*np.sin(lon0_rad), np.sin(lat0_rad)]
])
# 应用旋转
enu = R @ np.array([dx, dy, dz])
return enu[0], enu[1], enu[2]
小技巧:我写代码时喜欢把经纬度转弧度放在函数开头,避免后面反复计算。另外,WGS84椭球参数建议写成常量,别硬编码在公式里。
1.4 逆变换:ENU回ECEF
反过来也一样简单。旋转矩阵有个好性质——它是正交矩阵,逆矩阵等于转置矩阵。所以ENU到ECEF的转换就是:
def enu_to_ecef(e, n, u, lat0, lon0, h0):
lat0_rad = np.radians(lat0)
lon0_rad = np.radians(lon0)
# 站心ECEF坐标(同上)
a = 6378137.0
e2 = 0.00669437999014
N = a / np.sqrt(1 - e2 * np.sin(lat0_rad)**2)
x0 = (N + h0) * np.cos(lat0_rad) * np.cos(lon0_rad)
y0 = (N + h0) * np.cos(lat0_rad) * np.sin(lon0_rad)
z0 = (N * (1 - e2) + h0) * np.sin(lat0_rad)
# 旋转矩阵的转置
R_T = np.array([
[-np.sin(lon0_rad), -np.sin(lat0_rad)*np.cos(lon0_rad), np.cos(lat0_rad)*np.cos(lon0_rad)],
[ np.cos(lon0_rad), -np.sin(lat0_rad)*np.sin(lon0_rad), np.cos(lat0_rad)*np.sin(lon0_rad)],
[ 0, np.cos(lat0_rad), np.sin(lat0_rad)]
])
# 先旋转,再平移
xyz = R_T @ np.array([e, n, u]) + np.array([x0, y0, z0])
return xyz[0], xyz[1], xyz[2]
1.5 避坑指南
我曾经踩过的坑:
- 经纬度顺序搞反:很多API里纬度在前、经度在后,但有些库是反的。我建议统一用(lat, lon)顺序,并在注释里写清楚。
- 弧度vs角度:numpy的三角函数默认用弧度。我曾经直接把角度传进去,结果算出来的坐标差了十万八千里。
- 海拔单位:ECEF坐标单位是米,海拔也是米。如果遇到英尺的数据,记得先转换。
1.6 知识体系总览
下面这张图把本章的核心逻辑串起来了。你看一遍就能明白整个转换流程:
1.7 实际应用场景
这个转换在工程里用得特别多。我举几个例子:
- RTK定位:基准站输出ECEF坐标,流动站需要ENU坐标才能知道相对位置。
- 无人机导航:起飞点作为ENU原点,实时计算无人机相对于起飞点的偏移。
- 车载导航:将卫星定位结果转换到车辆所在的局部坐标系,方便路径规划。
你想想看,如果没有这个转换,卫星定位的结果(ECEF)跟实际使用场景(ENU)之间就隔着一道墙。搞懂了转换原理,这道墙就打通了。
好了,这一章就到这里。代码可以直接拿去用,遇到问题欢迎交流。
可可精品资料库