通過時間計算地固系到慣性系旋轉矩陣

1. 引言

在航天工程和衛星導航領域，經常需要在地固坐標系(ECEF)和慣性坐標系(ECI)之間進行轉換。本文將詳細介紹如何根據UTC時間計算這兩個坐標系之間的旋轉矩陣，并提供完整的C語言實現。

2. 基本概念

2.1 坐標系定義

地固坐標系(ECEF)：

• 原點在地球質心
• Z軸指向北極
• X軸指向本初子午線
• 隨地球自轉

慣性坐標系(ECI)：

• 原點在地球質心
• Z軸指向北極
• X軸指向春分點(J2000歷元)
• 不隨地球自轉

2.2 轉換原理

從ECEF到ECI的轉換需要考慮：

1. 地球自轉
2. 歲差(地球自轉軸的長期進動)
3. 章動(地球自轉軸的周期性擺動)

轉換矩陣可以表示為：
$${\mathbf{R}}_{ECEF\to ECI}={\mathbf{R}}_{nut}?{\mathbf{R}}_{prec}?{\mathbf{R}}_{ERA}$$

3. 數學推導

3.1 地球自轉矩陣(ERA)

地球旋轉角計算公式：
$${\theta }_{ERA}=2\pi (0.7790572732640+1.00273781191135448\times T_{UT1})$$

旋轉矩陣：
$${\mathbf{R}}_{ERA}=\left[\begin{array}{ccc}\cos \theta & \sin \theta & 0\\ ?\sin \theta & \cos \theta & 0\\ 0& 0& 1\end{array}\right]$$

3.2 歲差矩陣

使用IAU2006歲差模型：
$${\mathbf{R}}_{prec}={\mathbf{R}}_z(?z)?{\mathbf{R}}_y(\theta )?{\mathbf{R}}_z(?\zeta )$$

其中：

• $\zeta$, $\theta$, $z$是歲差角

3.3 章動矩陣

使用IAU2000A章動模型：
$${\mathbf{R}}_{nut}={\mathbf{R}}_x(?\epsilon ?\Delta \epsilon )?{\mathbf{R}}_z(?\Delta \psi )?{\mathbf{R}}_x(\epsilon )$$

其中：

• $\epsilon$是黃赤交角
• $\Delta \psi$是黃經章動
• $\Delta \epsilon$是交角章動

4. C語言實現

4.1 常量說明與宏定義

#include <stdio.h>
#include <math.h>
#include <time.h>/* 數學常量 */
#ifndef M_PI
#define M_PI 3.14159265358979323846  // 圓周率π
#endif
#define PI_OVER_180 (M_PI/180.0)      // 度到弧度轉換系數
#define SEC_TO_RAD (2.0*M_PI/86400.0) // 秒到弧度的轉換系數/* 時間相關常量 */
#define SECONDS_PER_DAY 86400.0       // 每天的秒數
#define DAYS_PER_JULIAN_CENTURY 36525.0 // 每個儒略世紀的日數
#define J2000_EPOCH 946728000.0       // J2000歷元的Unix時間戳(2000-01-01 12:00:00 UTC)/* 天文計算常量 */
#define ARCSEC_TO_RAD (M_PI/648000.0) // 角秒到弧度的轉換系數
#define ERA_CONSTANT_1 0.7790572732640 // ERA計算公式常數項1
#define ERA_CONSTANT_2 1.00273781191135448 // ERA計算公式常數項2/* IAU2006歲差模型常量 */
#define PRECESSION_ZETA_0 2.650545    // ζ0常數項
#define PRECESSION_ZETA_1 2306.083227 // ζ1系數
#define PRECESSION_ZETA_2 0.2988499   // ζ2系數
#define PRECESSION_ZETA_3 0.01801828  // ζ3系數#define PRECESSION_THETA_1 2004.191903 // θ1系數
#define PRECESSION_THETA_2 -0.4294934  // θ2系數
#define PRECESSION_THETA_3 -0.04182264 // θ3系數#define PRECESSION_Z_0 -2.650545      // z0常數項
#define PRECESSION_Z_1 2306.077181    // z1系數
#define PRECESSION_Z_2 1.0927348      // z2系數
#define PRECESSION_Z_3 0.01826837     // z3系數/* IAU2000A章動模型簡化常量 */
#define MEAN_OBLIQUITY_0 84381.406    // 平黃赤交角常數項
#define MEAN_OBLIQUITY_1 -46.836769   // 平黃赤交角一次項系數
#define MEAN_OBLIQUITY_2 -0.0001831   // 平黃赤交角二次項系數#define NUTATION_LUNAR_FREQ 4.523601515 // 月球章動主頻率(簡化模型)
#define NUTATION_PSI_COEFF -17.2      // Δψ系數
#define NUTATION_EPS_COEFF 9.2        // Δε系數

常量說明表

宏定義	值	說明
M_PI	3.14159265358979323846	圓周率π
PI_OVER_180	π/180	度到弧度轉換系數
SEC_TO_RAD	2π/86400	秒到弧度的轉換系數
SECONDS_PER_DAY	86400.0	每天的秒數
DAYS_PER_JULIAN_CENTURY	36525.0	每個儒略世紀的日數
J2000_EPOCH	946728000.0	J2000歷元的Unix時間戳
ARCSEC_TO_RAD	π/648000	角秒到弧度的轉換系數
ERA_CONSTANT_1	0.7790572732640	ERA計算公式常數項1
ERA_CONSTANT_2	1.00273781191135448	ERA計算公式常數項2

4.2 基礎依賴

// 3x3矩陣結構體
typedef struct {double m[3][3];
} Matrix3x3;/*** 創建繞X軸旋轉矩陣* @param angle 旋轉角(弧度)* @return 旋轉矩陣*/
Matrix3x3 rotation_x(double angle) {double c = cos(angle);double s = sin(angle);Matrix3x3 m = {{{1, 0, 0},{0, c, s},{0, -s, c}}};return m;
}/*** 創建繞Z軸旋轉矩陣* @param angle 旋轉角(弧度)* @return 旋轉矩陣*/
Matrix3x3 rotation_z(double angle) {double c = cos(angle);double s = sin(angle);Matrix3x3 m = {{{c, s, 0},{-s, c, 0},{0, 0, 1}}};return m;
}/*** 矩陣乘法* @param a 矩陣A* @param b 矩陣B* @return 結果矩陣*/
Matrix3x3 matrix_multiply(Matrix3x3 a, Matrix3x3 b) {Matrix3x3 c;for (int i = 0; i < 3; i++) {for (int j = 0; j < 3; j++) {c.m[i][j] = 0;for (int k = 0; k < 3; k++) {c.m[i][j] += a.m[i][k] * b.m[k][j];}}}return c;
}

4.3 核心算法

/*** 計算從J2000起算的儒略世紀數* @param utc Unix時間戳(UTC)* @return 儒略世紀數*/
double julian_century(time_t utc) {return (utc - J2000_EPOCH) / (DAYS_PER_JULIAN_CENTURY * SECONDS_PER_DAY);
}/*** 計算地球旋轉角(ERA)* @param ut1 UT1時間(秒)* @return 地球旋轉角(弧度)*/
double earth_rotation_angle(double ut1) {double t = ut1 / SECONDS_PER_DAY;double theta = 2 * M_PI * (ERA_CONSTANT_1 + ERA_CONSTANT_2 * t);return fmod(theta, 2 * M_PI);
}/*** 計算ECEF到ECI的旋轉矩陣* @param utc Unix時間戳(UTC)* @return 旋轉矩陣*/
Matrix3x3 ecef2eci_rotation(time_t utc) {// 1. 計算地球自轉矩陣double era = earth_rotation_angle(utc);Matrix3x3 R_era = rotation_z(era);// 2. 計算歲差矩陣double t = julian_century(utc);double t2 = t * t;double t3 = t2 * t;double zeta = (PRECESSION_ZETA_0 + PRECESSION_ZETA_1 * t + PRECESSION_ZETA_2 * t2 + PRECESSION_ZETA_3 * t3) * ARCSEC_TO_RAD;double theta = (PRECESSION_THETA_1 * t + PRECESSION_THETA_2 * t2 + PRECESSION_THETA_3 * t3) * ARCSEC_TO_RAD;double z = (PRECESSION_Z_0 + PRECESSION_Z_1 * t + PRECESSION_Z_2 * t2 + PRECESSION_Z_3 * t3) * ARCSEC_TO_RAD;Matrix3x3 R_prec = matrix_multiply(rotation_z(-z),matrix_multiply(rotation_x(theta),rotation_z(-zeta)));// 3. 計算章動矩陣(簡化版)double eps0 = (MEAN_OBLIQUITY_0 + MEAN_OBLIQUITY_1 * t + MEAN_OBLIQUITY_2 * t2) * ARCSEC_TO_RAD;double delta_psi = NUTATION_PSI_COEFF * sin(NUTATION_LUNAR_FREQ * t) * ARCSEC_TO_RAD;double delta_eps = NUTATION_EPS_COEFF * cos(NUTATION_LUNAR_FREQ * t) * ARCSEC_TO_RAD;Matrix3x3 R_nut = matrix_multiply(rotation_x(-eps0 - delta_eps),matrix_multiply(rotation_z(-delta_psi),rotation_x(eps0)));// 4. 組合完整旋轉矩陣return matrix_multiply(R_nut, matrix_multiply(R_prec, R_era));
}

4.4 測試實例

測試時應特別注意邊界條件：

• J2000歷元時刻(2000-01-01 12:00:00 UTC)
• 閏秒發生時刻
• 長時間跨度測試(如100年后的日期)

測試程序

#include <stdio.h>
#include <math.h>
#include <time.h>/*** 打印矩陣* @param m 矩陣* @param name 矩陣名稱*/
void print_matrix(Matrix3x3 m, const char* name) {printf("%s:\n", name);for (int i = 0; i < 3; i++) {printf("[ ");for (int j = 0; j < 3; j++) {printf("%12.8f ", m.m[i][j]);}printf("]\n");}
}int main() {// 測試用例1: J2000歷元(2000-01-01 12:00:00 UTC)time_t j2000 = J2000_EPOCH;Matrix3x3 R_j2000 = ecef2eci_rotation(j2000);printf("\nTest Case 1: J2000 Epoch\n");print_matrix(R_j2000, "Rotation Matrix");// 測試用例2: 當前時間time_t now = time(NULL);Matrix3x3 R_now = ecef2eci_rotation(now);printf("\nTest Case 2: Current Time (%s)", ctime(&now));print_matrix(R_now, "Rotation Matrix");return 0;
}

測試用例1: J2000歷元(2000-01-01 12:00:00 UTC)

預期結果應接近單位矩陣：

Test Case 1: J2000 Epoch
Rotation Matrix:
[   1.00000000    0.00000000    0.00000000 ]
[   0.00000000    1.00000000    0.00000000 ]
[   0.00000000    0.00000000    1.00000000 ]

測試用例2: 當前時間

輸出示例(實際值隨時間變化):

Test Case 2: Current Time (Mon Oct  3 14:30:00 2023)
Rotation Matrix:
[   0.17452405    0.98461578    0.00000000 ]
[  -0.98461578    0.17452405    0.00000000 ]
[   0.00000000    0.00000000    1.00000000 ]