matlab-convert-aerospace-coordinates

航空航天基础知识

航空航天工具箱的核心函数，用于单位转换、时间转换、坐标变换和旋转表示。

适用场景

• 在航空航天单位制（SI、英制、海里制）之间进行转换
• 计算儒略日、修正儒略日、十进制年或 TDB 时间
• 在 ECEF、ECI、LLA 或平地坐标系之间变换位置
• 在大地纬度和地心纬度之间转换
• 构建用于坐标系变换的 DCM（ECI↔ECEF、ECEF↔NED、机体↔风轴、机体↔稳定轴）
• 在欧拉角、DCM、四元数和 Rodrigues 向量之间转换
• 使用 quaternion 对象进行旋转运算、插值或组合
• 确定数据所采用的坐标系

不适用场景

• 大气模型（atmosisa、atmoscoesa）——请使用 matlab-compute-aerospace-environment
• 空速修正（correctairspeed）——此处未涵盖
• 轨道传播或卫星机动——请使用 Aerospace Toolbox、Aerospace Blockset 或 Satellite Communications Toolbox
• 导航专用变换（lla2enu、lla2ned）——请使用 Navigation Toolbox
• Simulink 模块——请使用 /model-based-design-core:building-simulink-models 和 Aerospace Blockset

工作流

1. 确定坐标系——确定数据所在的坐标系以及目标坐标系。有关坐标系定义和选择指南，请参阅 references/coordinate-systems.md。
2. 先转换单位——调用变换函数前，确保输入单位符合函数要求。
3. 应用变换——使用合适的函数，并确保参数顺序正确。
4. 验证——将结果转换回原始坐标系；误差应小于 1e-10。

主要函数

单位转换

函数	转换量	单位
convlength	长度	'ft'、'm'、'km'、'in'、'mi'、'naut mi'
convvel	速度	'ft/s'、'm/s'、'km/s'、'in/s'、'km/h'、'mph'、'kts'、'ft/min'
convang	角度	'deg'、'rad'、'rev'
convacc	加速度	'ft/s^2'、'm/s^2'、'km/s^2'、'in/s^2'、'km/h-s'、'mph/s'、'G''s'
convangacc	角加速度	'deg/s^2'、'rad/s^2'、'rpm/s'
convangvel	角速度	'deg/s'、'rad/s'、'rpm'
convforce	力	'N'、'lbf'
convmass	质量	'kg'、'lbm'、'slug'
convpres	压力	'Pa'、'psi'、'psf'、'atm'
convtemp	温度	'K'、'R'、'F'、'C'
convdensity	密度	'kg/m^3'、'slug/ft^3'、'lbm/ft^3'、'lbm/in^3'

所有转换函数使用相同的签名：output = convXXX(value, fromUnit, toUnit)

时间转换

函数	用途	起始版本
juliandate	日历日期 → 儒略日	R2006b
mjuliandate	日历日期 → 修正儒略日（JD − 2400000.5）	R2006b
decyear	日历日期 → 十进制年	R2006b
leapyear	判断年份是否为闰年	R2006b
tdbjuliandate	地球时 → TDB 儒略日	R2015a

坐标变换

函数	起点	终点	起始版本
lla2ecef	LLA（大地坐标）	ECEF	R2006b
ecef2lla	ECEF	LLA（大地坐标）	R2006b
lla2eci	LLA	ECI	R2014a
eci2lla	ECI	LLA	R2014a
ecef2eci	ECEF（位置/速度/加速度）	ECI	R2019a
eci2ecef	ECI（位置/速度/加速度）	ECEF	R2019a
eci2aer	ECI	AER（方位角、仰角、距离）	R2015a
lla2flat	LLA	平地坐标系	R2011a
flat2lla	平地坐标系	LLA	R2011a
geod2geoc	大地纬度	地心纬度	R2006b
geoc2geod	地心纬度	大地纬度	R2006b
ijk2keplerian	ECI 位置/速度	开普勒轨道根数	R2019a
keplerian2ijk	开普勒轨道根数	ECI 位置/速度	R2019a
siderealTime	UTC 儒略日	GMST（度）/GAST（秒）	R2021a
dcmeci2ecef	—	ECI 到 ECEF 的 DCM	R2013b
dcmecef2ned	—	ECEF 到 NED 的 DCM	R2006b
dcm2latlon	ECEF 到 NED 的 DCM	纬度/经度	R2006b
dcmbody2wind	Alpha、Beta	机体到风轴的 DCM	R2006b
dcm2alphabeta	机体到风轴的 DCM	Alpha、Beta	R2006b
dcmbody2stability	Alpha	机体到稳定轴的 DCM	R2022a

旋转表示

函数	起点	终点	起始版本
angle2dcm	欧拉角	DCM	R2006b
dcm2angle	DCM	欧拉角	R2006b
angle2quat	欧拉角	四元数（1×4）	R2006b
quat2angle	四元数（1×4）	欧拉角	R2007b
dcm2quat	DCM	四元数（1×4）	R2006b
quat2dcm	四元数（1×4）	DCM	R2006b
angle2rod	欧拉角	Rodrigues 向量	R2017a
rod2angle	Rodrigues 向量	欧拉角	R2017a
dcm2rod	DCM	Rodrigues 向量	R2017a
rod2dcm	Rodrigues 向量	DCM	R2017a
quat2rod	四元数（1×4）	Rodrigues 向量	R2017a
rod2quat	Rodrigues 向量	四元数（1×4）	R2017a

四元数对象

方法	用途
quaternion(E,'eulerd',RS,PF)	根据欧拉角（度）创建
quaternion(E,'euler',RS,PF)	根据欧拉角（弧度）创建
quaternion(RM,'rotmat',PF)	根据旋转矩阵创建
quaternion(RV,'rotvec')	根据旋转向量（弧度）创建
compact(q)	提取 [w x y z] 数组
eulerd(q,RS,PF)	转换为欧拉角（度）
euler(q,RS,PF)	转换为欧拉角（弧度）
rotmat(q,PF)	转换为旋转矩阵
rotvec(q) / rotvecd(q)	转换为旋转向量（弧度/度）
rotatepoint(q,pts)	旋转点（主动旋转）
rotateframe(q,pts)	旋转坐标系（被动旋转）
normalize(q)	归一化为单位四元数
slerp(q1,q2,t)	球面线性插值
meanrot(q)	数组的平均旋转
dist(q1,q2)	角距离（弧度）
angvel(q,dt,PF)	根据四元数数组计算角速度
randrot(n)	均匀随机旋转

四元数数学运算（基于数组）

函数	用途	起始版本
quatmultiply(q,r)	四元数乘积（组合旋转）	R2006b
quatconj(q)	共轭（向量部分取反）	R2006b
quatinv(q)	逆（共轭 / 范数²）	R2006b
quatnormalize(q)	归一化为单位四元数	R2006b
quatnorm(q)	平方范数（q·q）	R2006b
quatmod(q)	模（范数的平方根）	R2006b
quatrotate(q,v)	用四元数旋转向量	R2006b
quatdivide(q,r)	四元数除以四元数	R2006b
quatinterp(p,q,f,method)	插值（'slerp'、'lerp'、'nlerp'）	R2016a
quatexp(q)	四元数的指数	R2016a
quatlog(q)	四元数的自然对数	R2016a
quatpower(q,pow)	四元数的幂	R2016a

模式

单位换算

% Always: convXXX(value, 'from', 'to')
alt_m = convlength(35000, 'ft', 'm');        % 10668.0 m
speed_ms = convvel(250, 'kts', 'm/s');       % 128.61 m/s
angle_rad = convang(45, 'deg', 'rad');       % 0.7854 rad
accel_g = convacc(9.81, 'm/s^2', 'G''s');   % 1.0 G
omega_rpm = convangvel(360, 'deg/s', 'rpm'); % 60 rpm

时间换算

% Julian Date from components (year, month, day, hour, min, sec)
jd = juliandate(2024, 6, 15, 12, 0, 0);     % 2460477.0

% Modified Julian Date
mjd = mjuliandate(2024, 6, 15, 12, 0, 0);   % 60476.5

% From datetime objects
dt = datetime(2024, 6, 15, 12, 0, 0);
jd = juliandate(dt);

% Decimal year
dy = decyear(2024, 6, 15);                   % 2024.4536

% TDB Julian Date from Terrestrial Time [yr mo day hr min sec]
tt = [2024 6 15 12 0 0];
jdTDB = tdbjuliandate(tt);

LLA ↔ ECEF

% LLA is [latitude_deg, longitude_deg, altitude_m]
lla = [40, -74, 0];
ecef = lla2ecef(lla);           % [1348613.0, -4703172.4, 4077985.6] m

% Round-trip verification
lla_check = ecef2lla(ecef);     % [40.0, -74.0, 0.0]

LLA ↔ ECI（与时间相关）

% ECI transforms require UTC time
lla = [40, -74, 1000];
utc = [2024 6 15 12 0 0];
posECI = lla2eci(lla, utc);

% Back to LLA
lla_check = eci2lla(posECI, utc);

开普勒根数 ↔ ECI (IJK)

% Convert ECI position/velocity to Keplerian orbital elements
rijk = [-2981784; 5207055; 3161595];   % ECI position (m), 3x1
vijk = [-3384; -4887; 4843];           % ECI velocity (m/s), 3x1
[a, ecc, incl, RAAN, argp, nu] = ijk2keplerian(rijk, vijk);
% a=m, angles in degrees

% Convert Keplerian elements back to ECI position/velocity
[rijk_back, vijk_back] = keplerian2ijk(a, ecc, incl, RAAN, argp, nu);

% For a different central body (e.g., Mars)
[a_mars, ecc_mars, incl_mars, RAAN_mars, argp_mars, nu_mars] = ...
    ijk2keplerian(rijk, vijk, CentralBody="Mars");

恒星时

% Compute Greenwich Mean Sidereal Time (GMST) and Apparent (GAST)
utcJD = juliandate(2019, 1, 4, 12, 0, 0);  % UTC as Julian date
dUT1 = 0;    % UT1-UTC offset (seconds)
dAT  = 37;   % TAI-UTC offset (leap seconds)
[thGMST, thGAST] = siderealTime(utcJD, dUT1, dAT);
% thGMST in degrees, thGAST in seconds

平地近似

% Good for short-range simulations (< ~100 km from reference)
llo = [40, -74];    % reference lat/lon (deg)
psio = 0;           % angular direction of flat Earth x-axis (rad, 0=North)
href = 0;           % reference height (m)

lla_point = [40.01, -73.99, 100];
flatPos = lla2flat(lla_point, llo, psio, href);  % [x, y, z] in meters
lla_back = flat2lla(flatPos, llo, psio, href);

欧拉角 ↔ DCM

% CRITICAL: angle2dcm expects RADIANS, not degrees
yaw = deg2rad(30); pitch = deg2rad(10); roll = deg2rad(5);
dcm = angle2dcm(yaw, pitch, roll, 'ZYX');

% Extract angles back (returns radians)
[y, p, r] = dcm2angle(dcm, 'ZYX');
fprintf('Yaw=%.1f, Pitch=%.1f, Roll=%.1f deg\n', rad2deg(y), rad2deg(p), rad2deg(r));

DCM ↔ 四元数（Aerospace Toolbox 格式）

% Aerospace Toolbox quaternion format: [q0 q1 q2 q3] = [scalar, vector]
q = dcm2quat(dcm);       % 1x4, scalar-first
dcm_back = quat2dcm(q);  % 3x3

% Verify orthogonality
err = max(abs(dcm' * dcm - eye(3)), [], 'all');
assert(err < 1e-14, 'DCM is not orthogonal');

四元数对象（现代方法）

% Create from Euler angles — specify sequence AND frame type
q = quaternion([30 10 5], 'eulerd', 'ZYX', 'frame');

% Rotate a point (active rotation)
pt = [1 0 0];
pt_rotated = rotatepoint(q, pt);   % [0.8529, 0.4924, -0.1736]

% Compose rotations by multiplication
q_total = q2 * q1;  % applies q1 first, then q2

% Always normalize after arithmetic accumulation
q = normalize(q);

四元数插值（SLERP）

% Smooth interpolation between two orientations
q1 = quaternion([0 0 0], 'eulerd', 'ZYX', 'frame');
q2 = quaternion([90 0 0], 'eulerd', 'ZYX', 'frame');
t = linspace(0, 1, 5)';
q_interp = slerp(q1, q2, t);

% Verify: yaw progresses linearly for pure yaw rotation
e = eulerd(q_interp, 'ZYX', 'frame');
% e(:,1) = [0, 22.5, 45, 67.5, 90]

四元数运算（基于数组的函数）

% Compose two rotations: q2 applied after q1
q1 = angle2quat(deg2rad(30), 0, 0, 'ZYX');
q2 = angle2quat(0, deg2rad(10), 0, 'ZYX');
q_total = quatmultiply(q2, q1);  % q2 * q1 (apply q1 first)

% Inverse rotation
q_inv = quatinv(q1);
q_identity = quatmultiply(q1, q_inv);  % [1 0 0 0]

% Rotate a vector
v = [1 0 0];
v_rot = quatrotate(q1, v);  % rotate v by q1

% Normalize after accumulation
q_accumulated = quatnormalize(q_total);

% Interpolation (slerp, lerp, nlerp)
q_mid = quatinterp(q1, q2, 0.5, 'slerp');

% Relative rotation: what rotation takes q1 to q2?
q_rel = quatdivide(q2, q1);  % q_rel such that q2 = q_rel * q1

% Exponential/logarithm (useful for angular velocity integration)
q_log = quatlog(q1);         % maps to tangent space
q_back = quatexp(q_log);     % back to quaternion

% Fractional rotation (half the rotation of q1)
q_half = quatpower(q1, 0.5);

Rodrigues 向量

% Compact 3-element representation (singular at 180 deg)
rod = dcm2rod(dcm);           % 1x3 vector
dcm_back = rod2dcm(rod);

% Convert between all representations
rod = quat2rod(q_array);      % q is M×4 (scalar-first)
q_back = rod2quat(rod);

机体坐标系、风轴系和稳定坐标系

% Body-to-wind DCM from angle of attack and sideslip
alpha = deg2rad(5);
beta = deg2rad(2);
dcm_bw = dcmbody2wind(alpha, beta);

% Extract alpha/beta from a DCM
[alpha_out, beta_out] = dcm2alphabeta(dcm_bw);

% Body-to-stability (alpha only, no sideslip)
dcm_bs = dcmbody2stability(alpha);

ECEF 到 NED 坐标系 DCM

% Get the DCM to rotate vectors from ECEF to local NED
lat = deg2rad(40);
lon = deg2rad(-74);
dcm_ecef2ned = dcmecef2ned(lat, lon);

% Transform an ECEF velocity to NED
v_ecef = [10; 20; 30];
v_ned = dcm_ecef2ned * v_ecef;

ECI 到 AER（方位角、仰角、斜距）

% Compute look angles from a ground station to a satellite in ECI
posECI = [-2981784, 5207055, 3161595];  % satellite ECI position (m)
utc = [2019 1 4 12 0 0];               % observation time
lla0 = [28.5, -80.5, 0];               % ground station [lat, lon, alt] (deg, deg, m)

aer = eci2aer(posECI, utc, lla0);
% aer = [azimuth_deg, elevation_deg, slant_range_m]

约定

• angle2dcm 和 angle2quat 要求输入弧度——请务必先使用 deg2rad() 或 convang 转换
• 默认旋转顺序为 'ZYX'（偏航-俯仰-横滚）——请务必明确指定，以免产生歧义
• Aerospace Toolbox 四元数格式为 [scalar, i, j, k]——1×4 数组函数采用标量优先格式
• quaternion 对象也采用标量优先格式（compact 返回 [w x y z]）
• rotatepoint 与 rotateframe——rotatepoint 旋转点（主动旋转）；rotateframe 旋转坐标系（被动旋转）。两者互为逆运算。
• LLA 的顺序为 [lat, lon, alt]，单位分别为度和米
• geod2geoc 和 geoc2geod 的纬度参数要求弧度，高度/半径单位为米
• ECI 函数需要 UTC 时间——结果会随地球自转而变化
• 四元数算术运算后应进行归一化，以防止漂移
• Rodrigues 向量在 180° 处存在奇异性——任意旋转请使用四元数

常见错误

错误	错误原因	正确做法
向 angle2dcm 传入角度值	函数要求弧度	使用 deg2rad() 或 convang(val,'deg','rad')
省略旋转顺序	默认的 ZYX 顺序可能不符合你的约定	始终明确传入 'ZYX'、'XYZ' 等
认为四元数格式为 [x y z w]	Aerospace Toolbox 使用 [w x y z]（标量优先）	检查格式；使用 compact(q) 进行验证
需要 rotateframe 时却使用 rotatepoint	两者互为逆运算——用错会导致旋转方向相反	主动旋转使用 rotatepoint；被动旋转使用 rotateframe
四元数运算后跳过归一化	四元数漂移会导致模长不为 1，从而使旋转失真	累积旋转后调用 normalize(q)
向 geod2geoc 传入角度值	纬度参数要求弧度	转换为弧度：geod2geoc(deg2rad(lat), alt)
ECI 变换时忘记提供 UTC 时间	ECI 位置取决于该时刻地球的自转状态	始终提供 [yr mo day hr min sec]
在接近 180° 旋转时使用 Rodrigues 向量	Rodrigues 向量在 π 处存在奇异性	改用四元数表示
忽略接近 ±90° 俯仰角时的万向节锁	欧拉角提取会丢失一个自由度	计算时直接使用四元数或 DCM

万向节锁

当第二次旋转达到 ±90° 时，欧拉角表示会丢失一个自由度（对于 ZYX：俯仰角 = ±90°）。其表现包括：

• 俯仰角接近 ±90° 时，dcm2angle 会返回意外的偏航角/横滚角
• 在奇异点附近插值欧拉角会产生不稳定的轨迹

解决方案： 使用四元数进行计算和插值。仅在显示或便于人理解时才转换为欧拉角。

% Detect gimbal lock risk
[~, pitch, ~] = dcm2angle(dcm, 'ZYX');
if abs(pitch) > deg2rad(85)
    warning('Near gimbal lock — use quaternion representation');
end

参考资料

• 有关详细的坐标系定义（ECEF、ECI、NED、机体、风轴、稳定）、轴约定以及坐标系选择指南，请参阅 references/coordinate-systems.md。
• 有关详细的旋转数学内容，请参阅 references/rotation-conventions.md：DCM 属性、四元数代数、万向节锁理论，以及各种表示形式之间的转换路径。
• 有关完整的单位字符串表和链式转换，请参阅 references/unit-conversion-guide.md。

---

版权所有 2026 The MathWorks, Inc.

---