一、ECEF與ENU坐標系定義

1. ECEF坐標系（地心地固坐標系）

• 原點：地球質心
• X軸：指向本初子午線與赤道交點
• Y軸：在赤道平面內與X軸垂直
• Z軸：指向北極
• 數學表示：$P_{ecef}=(x,y,z)$

2. ENU坐標系（東北天坐標系）

• 原點：本地參考點
• E軸：指向東方
• N軸：指向北方
• U軸：垂直地面向上
• 數學表示：$P_{enu}=(e,n,u)$

二、坐標轉換數學原理

1. 地理坐標轉ECEF：
   $$?\begin{array}{l}x=(N+h)\cos ?\cos \lambda \\ y=(N+h)\cos ?\sin \lambda \\ z=(N(1?e^2)+h)\sin ?\end{array}$$
   其中：

• $N=\frac{a}{\sqrt{1?e^2{\sin }^2?}}$（卯酉圈曲率半徑）
• $e^2=2f?f^2$（橢球偏心率平方）
• $a=6378137m$（WGS84長半軸）
• $f=1/298.257223563$（WGS84扁率）

2. ECEF轉ENU：
   $$\left[\begin{array}{c}e\\ n\\ u\end{array}\right]=\left[\begin{array}{ccc}?\sin \lambda & \cos \lambda & 0\\ ?\sin ?\cos \lambda & ?\sin ?\sin \lambda & \cos ?\\ \cos ?\cos \lambda & \cos ?\sin \lambda & \sin ?\end{array}\right]\left[\begin{array}{c}x?x_0\\ y?y_0\\ z?z_0\end{array}\right]$$

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

五、驗證結果分析

1. 誤差直方圖顯示各方向誤差分布
2. 散點圖展示平面誤差的相關性
3. 典型誤差應小于1e-4米（數值計算誤差）

該實現完整展示了從理論到實踐的全流程，可通過調整測試點數量和分布進行更嚴格的驗證。實際工程應用中需考慮坐標系的旋轉順序、大地水準面模型等更多細節。

（驗證代碼還有BUG，不能直接運行，但是我肝不了了，該天再調了，晚安，嗎喀巴卡）