第1章：WGS-84与ECEF坐标转换——从经纬高到地心地固的正算与反算

各位同学，今天我们来聊聊导航定位里最基础、也最绕不开的一个环节——WGS-84坐标系和ECEF坐标系之间的转换。

说白了，WGS-84就是我们平时说的经纬高（LLA），而ECEF是地心地固坐标系。你手机里的GPS芯片拿到的是WGS-84坐标，但做惯性导航、卫星轨道计算时，大家更习惯用ECEF。这两者怎么互转？嗯，这就是本章要讲的核心。

1.1 两个坐标系长什么样？

先简单回顾一下：

WGS-84（LLA）：用纬度φ、经度λ、海拔高度h表示。纬度范围-90°到90°，经度-180°到180°。地球被近似为一个椭球，赤道半径约6378137米，极半径约6356752米。
ECEF（地心地固）：原点在地球质心，X轴指向本初子午线与赤道交点，Y轴指向东经90°方向，Z轴指向北极。单位是米。

我个人习惯把ECEF想象成一个“固定在地球上的直角坐标系”，地球自转时它也跟着转。而WGS-84更像是一个“贴在地球表面的网格”。

1.2 正算：从LLA到ECEF

正算公式其实不复杂，但有个关键参数——椭球的第一偏心率平方e²。WGS-84椭球的参数是：

参数	符号	值
长半轴	a	6378137.0 m
短半轴	b	6356752.3142 m
第一偏心率平方	e²	0.00669437999014

公式如下：

N = a / sqrt(1 - e² * sin²(φ))
X = (N + h) * cos(φ) * cos(λ)
Y = (N + h) * cos(φ) * sin(λ)
Z = (N * (1 - e²) + h) * sin(φ)

这里N是卯酉圈曲率半径。你想想看，如果h=0，这个点就在椭球面上。我在项目中遇到过有人直接把地球当球体算，结果在赤道附近误差能到几十米——嗯，做高精度定位时千万别这么干。

1.3 反算：从ECEF到LLA

反算就有点意思了。因为公式里φ和h互相依赖，没法直接解。常用的方法是迭代法，或者用Bowring算法。

我推荐用迭代法，稳定且容易理解：

# 初始值
p = sqrt(X² + Y²)
φ = atan2(Z, p * (1 - e²))  # 初始猜测
h = 0

# 迭代
for i in range(10):
    N = a / sqrt(1 - e² * sin²(φ))
    h = p / cos(φ) - N
    φ = atan2(Z, p * (1 - e² * N / (N + h)))

一般迭代3-5次就收敛了。我曾经在写嵌入式代码时，为了省计算量只迭代了2次，结果在极地附近误差偏大——后来老老实实改成了5次。

1.4 Python实现

下面给出完整的正反算代码。我习惯把WGS-84参数定义成常量，避免硬编码。

import math

# WGS-84 椭球参数
A = 6378137.0          # 长半轴
E2 = 0.00669437999014  # 第一偏心率平方

def lla_to_ecef(lat_deg, lon_deg, h):
    """正算：经纬高 -> ECEF"""
    lat = math.radians(lat_deg)
    lon = math.radians(lon_deg)
    
    sin_lat = math.sin(lat)
    cos_lat = math.cos(lat)
    cos_lon = math.cos(lon)
    sin_lon = math.sin(lon)
    
    N = A / math.sqrt(1 - E2 * sin_lat * sin_lat)
    
    X = (N + h) * cos_lat * cos_lon
    Y = (N + h) * cos_lat * sin_lon
    Z = (N * (1 - E2) + h) * sin_lat
    
    return X, Y, Z

def ecef_to_lla(X, Y, Z, max_iter=10):
    """反算：ECEF -> 经纬高（迭代法）"""
    p = math.sqrt(X*X + Y*Y)
    if p < 1e-12:
        # 在极点附近
        lat = math.copysign(math.pi/2, Z)
        lon = 0.0
        h = abs(Z) - A * math.sqrt(1 - E2)
        return math.degrees(lat), math.degrees(lon), h
    
    lat = math.atan2(Z, p * (1 - E2))  # 初始猜测
    h = 0.0
    
    for _ in range(max_iter):
        sin_lat = math.sin(lat)
        N = A / math.sqrt(1 - E2 * sin_lat * sin_lat)
        h = p / math.cos(lat) - N
        lat = math.atan2(Z, p * (1 - E2 * N / (N + h)))
    
    lon = math.atan2(Y, X)
    return math.degrees(lat), math.degrees(lon), h

# 示例：北京天安门坐标（约）
lat, lon, h = 39.9042, 116.4074, 44.0
X, Y, Z = lla_to_ecef(lat, lon, h)
print(f"ECEF: {X:.2f}, {Y:.2f}, {Z:.2f}")

lat2, lon2, h2 = ecef_to_lla(X, Y, Z)
print(f"反算回LLA: {lat2:.6f}, {lon2:.6f}, {h2:.3f}")

小技巧：反算时如果p值很小（比如在南北极附近），直接用atan2(Z, p)会出问题。我一般会加个判断，p小于1米时单独处理。

1.5 知识体系结构图

下面这张图展示了本章的核心逻辑：

核心要点：正算有闭式解，直接代入公式即可。反算需要迭代，但收敛很快。记住WGS-84椭球参数不要用错，a和e²是固定的。

1.6 避坑指南

我总结几个实际项目中容易踩的坑：

角度单位：Python的math库三角函数默认弧度，别忘了转。我曾经见过同事直接把度数传进去，结果算出来的坐标偏了十万八千里。
高度h的单位：LLA里的h是海拔高度（椭球高），单位是米。有些数据源给的是大地高，要确认清楚。
反算收敛性：在赤道附近迭代很快，但在极地附近可能需要更多迭代次数。我一般设最大迭代10次，实际3-5次就够。
精度问题：单精度float32在计算大数时误差明显。做高精度定位时，建议用float64（Python默认就是双精度，没问题）。

注意：ECEF坐标的数值通常很大（地球半径约6378公里），X、Y、Z都是百万米量级。做数值计算时注意浮点精度，尤其是减法操作可能产生灾难性抵消。

好了，这一章的内容就到这里。正算和反算是导航定位的基石，后面讲惯性导航、卫星定位时都会用到。建议你把代码跑一遍，看看天安门、自由女神像、悉尼歌剧院这些地标的ECEF坐标是多少——嗯，亲手算一遍比看十遍都管用。

可可精品资料库

📚 课程目录

第1章：WGS-84与ECEF坐标转换——从经纬高到地心地固的正算与反算

1.1 两个坐标系长什么样？

1.2 正算：从LLA到ECEF

1.3 反算：从ECEF到LLA

1.4 Python实现

1.5 知识体系结构图

1.6 避坑指南