### WGS84与东北天坐标系的关系及转换方法
#### 1. WGS84与东北天坐标系简介
WGS84（World Geodetic System 1984）是一种全球性的大地测量坐标系统，用于定义地球表面的位置。它由经度、纬度和高度组成[^1]。
东北天坐标系（ENU，East-North-Up Coordinate System），是以某一点为中心建立的地方直角坐标系。其三个轴分别指向东（X轴）、北（Y轴）以及垂直向上方向（Z轴）。该坐标系通常被应用于局部区域内的导航和定位计算中[^3]。
---
#### 2. 转换原理
从WGS84到东北天坐标系的转换涉及两个主要步骤：
- **第一步**：将WGS84地理坐标（经纬度和海拔）转换为地心地固坐标系（ECEF，Earth-Centered Earth-Fixed Coordinate System）中的三维笛卡尔坐标。
- **第二步**：通过旋转矩阵将ECEF坐标投影至以指定点为中心的东北天坐标系下。
具体过程如下：
---
#### 3. 数学模型与算法实现
##### (1) WGS84 → ECEF 的转换公式
给定一个WGS84位置 \( (\phi, \lambda, h) \)，其中：
- \( \phi \): 纬度（单位为弧度）
- \( \lambda \): 经度（单位为弧度）
- \( h \): 高程（相对于椭球面的高度）
可以利用以下公式将其转化为ECEF坐标 \( (x, y, z) \):
\[
N(\phi) = \frac{a}{\sqrt{1-e^2\sin^2{\phi}}}
\]
\[
x = (N(\phi)+h)\cos{\phi}\cos{\lambda}
\]
\[
y = (N(\phi)+h)\cos{\phi}\sin{\lambda}
\]
\[
z = ((1-e^2)N(\phi)+h)\sin{\phi}
\]
其中：
- \( a \approx 6378137m \): 椭球半径；
- \( e^2 \approx 0.00669438 \): 扁率平方；
---
##### (2) ECEF → ENU 的转换公式
设中心点的地心坐标为 \( C(x_c, y_c, z_c) \)，对应的地理坐标为 \( (\phi_c, \lambda_c, h_c) \)。对于任意目标点P，在完成上述WGS84→ECEF变换后得到的目标点ECEF坐标为\( P(x_p, y_p, z_p) \)。
构建旋转矩阵R来表示从ECEF到ENU的方向变化关系：
\[
R =
\begin{bmatrix}
-\sin{\lambda} & -\sin{\phi}\cos{\lambda} & \cos{\phi}\cos{\lambda}\\
\cos{\lambda} & -\sin{\phi}\sin{\lambda} & \cos{\phi}\sin{\lambda}\\
0 & \cos{\phi} & \sin{\phi}
\end{bmatrix}
\]
最终ENUCoordinates可通过向量乘法获得：
\[
V_{enu}= R(V_{ecef}- V_{ref})
\]
这里，
- \( V_{enu}=[v_e,v_n,v_u]^T\) 表示东北天上坐标分量；
- \( V_{ecef}=[x,y,z]^T\) 是待求解点在ECEF下的表达形式；
- \( V_{ref}=[x_c,y_c,z_c]^T\) 则代表参考原点所在处对应于整个体系里的固定数值。
---
#### 4. 实现代码示例
以下是基于Python的一个简单实现案例:
```python
import numpy as np
def wgs84_to_ecef(lat, lon, alt):
"""
Convert from geodetic coordinates to ECEF.
Parameters:
lat(float): Latitude in radians
lon(float): Longitude in radians
alt(float): Altitude above the ellipsoid
Returns:
tuple: X,Y,Z coordinates of point in meters
"""
a = 6378137.0 # Semi-major axis length (equatorial radius), m
f = 1/298.257223563 # Flattening factor
b = (1-f)*a # Semi-minor axis length (polar radius)
N = a / np.sqrt(1-(f*(2-f))*np.sin(lat)**2)
x = (N + alt) * np.cos(lat) * np.cos(lon)
y = (N + alt) * np.cos(lat) * np.sin(lon)
z = (((1-f)**2*N) + alt) * np.sin(lat)
return x, y, z
def ecef_to_enu(ecef_point, ref_lat, ref_lon, ref_alt):
"""
Transform an ECEF coordinate into local East-North-Up frame.
Args:
ecef_point(tuple): Point's Cartesian coords [X,Y,Z], unit:meters
ref_lat(float): Reference latitude(rad)
ref_lon(float): Reference longitude(rad)
ref_alt(float): Reference altitude(meter)
Return:
enu_coords(list): List containing east,north and up components respectively.
"""
xecef, yecef, zecef = ecef_point
x_ref, y_ref, z_ref = wgs84_to_ecef(ref_lat, ref_lon, ref_alt)
delta_x = xecef-x_ref;delta_y=yecef-y_ref;delta_z=zecef-z_ref;
sin_lambda=np.sin(ref_lon); cos_lambda=np.cos(ref_lon);
sin_phi=np.sin(ref_lat); cos_phi=np.cos(ref_lat);
rotation_matrix=[
[-sin_lambda , cos_lambda , 0],
[-cos_lambda*sin_phi,-sin_lambda*sin_phi,cos_phi ],
[ cos_lambda*cos_phi,sin_lambda*cos_phi,sin_phi ]
]
xyz_local=np.dot(rotation_matrix,[delta_x,delta_y,delta_z])
return list(map(lambda v:"%.6f"%v,xyz_local))
if __name__ == "__main__":
phi_deg, lambda_deg, height_m = 39.914266, 116.403988, 0
phi_rad = np.radians(phi_deg)
lambda_rad = np.radians(lambda_deg)
# Step1: From WGS84 -> ECEF
ecef_result=wgs84_to_ecef(phi_rad, lambda_rad, height_m )
print(f"ECEF Coordinates:{list(map('%.6f'.format,ecef_result))}")
# Step2: Define reference location for transformation origin
center_phi_deg,center_lambda_deg,height_center=39.914266,116.403988,0;
center_phi_rad=np.radians(center_phi_deg);
center_lambda_rad=np.radians(center_lambda_deg);
# Perform final conversion step
enu_output=ecef_to_enu(ecef_result,center_phi_rad,center_lambda_rad,height_center)
print(f"ENU Coordinates:[East={enu_output[0]}m,North={enu_output[1]}m,Up={enu_output[2]}m]")
```
---
超级会员免费看
订阅专栏 解锁全文
扫一扫
2410
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
