最小二乘法的状态转移方程（2）

矢志航天的阿洪

已于 2024-04-07 11:44:16 修改

阅读量700

点赞数 12

文章标签： matlab 开发语言算法

于 2024-04-07 11:42:40 首次发布

本文链接：https://blog.csdn.net/weixin_57997461/article/details/137432795

版权

文章详细阐述了地球指向参数的文本格式，包括世界时UT1与协调世界时UTC的关系，以及格林尼治地方时、真恒星时和章动的计算方法，涉及Nutation模型、月球和太阳影响、进动和章动理论。同时介绍了从真日期坐标系到地球赤道和格林尼治子午线系统的转换矩阵。

摘要由CSDN通过智能技术生成

一、文本格式

通常给出的地球指向参数文本格式为：

1962 01 01 37665 -0.012700 0.213000 0.0326338 0.0017230 0.064261 0.006067 0.000000 0.000000 2

1962 01 01：公历日期

37665 ：儒略日

-0.0127 ：x极移，单位：角秒

0.213 ：y极移，单位：角秒

0.0326338：世界时UT1与协调世界时UTC的差

0.0017230：LOD 地球日长变化单位：s

0.064261 ：x向的章动值单位：角秒

0.006067 ：y向的章动值单位：角秒

0.000000 ：x向极移的变化量单位：角秒

0.000000 ：y向极移的变化量单位：角秒

2 ：国际原子时TAI与协调世界时UTC之间相差整数秒单位：s

角秒即为1°/3600，又已知pi=180°，求得1角秒约等于4.848136811e-6

二、应用

格林尼治地方平时M称为世界时，记为UT。天文经度 $\lambda$ 处的地方平时与世界时的关系为

$m-\textup{UT}=\lambda/15$

时间的单位是小时（h），经度的单位是度（°）。通过测站直接观测天体获得的世界时记为UT0，对应于瞬时极子午圈。修正极移引起子午圈变化后得到的世界时记为UT1，对应于平均极子午圈。

根据瞬时真天球坐标系与瞬时地球坐标系的转换关系，两者的差别仅是X轴的指向差一个格林尼治真恒星时。真恒星时的计算由平恒星时和赤经岁差引起。平恒星时的计算与UT1有关

$\textup{GMST}_{s}=24110.54841+8640184.812886T_0+\\1.002737909350795\textup{UT1}+(0.093104-6.2\times 10^{-6}T)T^2$

其中 $T_0$ 是该天UT1时刻的简略儒略日的整数部分与2000年1月1日12：00的简略儒略日之差除以100年即36525天， $T$ 为该天的简略儒略日与2000年1月1日12：00的简略儒略日之差除以100年即36525天

$T_0=\frac{\textup{[UT}1]-51544.5}{36525}\\ T=\frac{\textup{UT}1-51544.5}{36525}$

则格林尼治地方平时的求解即写为

$\textup{GMST}_{\textup{rad}}=2\pi\frac{\textup{GMST}_{s}}{86400}-[2\pi\frac{\textup{GMST}_{s}}{86400}]$

代码如下：

%--------------------------------------------------------------------------
%
% gmst: Greenwich Mean Sidereal Time
%
% Input:
%  Mjd_UT1    Modified Julian Date UT1
%
% Output:
%  gmstime	   GMST in [rad]
%
% Last modified:   2015/08/12   M. Mahooti
%
%--------------------------------------------------------------------------
function gmstime = gmst(Mjd_UT1)

Secs = 86400;                       % Seconds per day
MJD_J2000 = 51544.5;

Mjd_0 = floor(Mjd_UT1);
UT1   = Secs*(Mjd_UT1-Mjd_0);       % [s]
T_0   = (Mjd_0  -MJD_J2000)/36525;
T     = (Mjd_UT1-MJD_J2000)/36525;

gmst  = 24110.54841 + 8640184.812866*T_0 + 1.002737909350795*UT1...
        + (0.093104-6.2e-6*T)*T*T;  % [s]

gmstime = 2*pi*Frac(gmst/Secs);     % [rad], 0..2pi

地球自转轴的方向还受到被称为章动的小周期扰动的影响。它们是由于在进差处理中平均的月球和太阳扭矩的月和年变化。对章动的主要贡献来自于月球轨道相对于地球赤道的不同方向，这一点由月球的升交点黄经表示，它引起了周期性的变化。在本文给出IAU1980章动模型，其表达式如下式所示

$\Delta\mathit{\Psi}=\sum_{i=1}^{106}(\Delta\mathit{\Psi})_i\cdot \sin(\phi_i )\\ \Delta\mathit{\varepsilon }=\sum_{i=1}^{106}(\Delta\mathit{\varepsilon })_i\cdot \cos(\phi_i )$

共106项，其常数计入如下形式。

每个术语都描述了月球和太阳轨道的平均元素的周期函数

$\phi _i=p_{l,i}l+p_{l',i}l'+p_{F,i}F+p_{D,i}D+p_{\Omega,i}\Omega$

整数系数为 $p_{l,i},p_{l',i},p_{F,i},p_{D,i},p_{\Omega,i}$ 。其它参数是月球的平近地点角 $l$ ，太阳平近地点角 $l'$ ，月球的升交点赤经 $F$ ，太阳和月球的平均经度之间的差异 $D$ ，和月球轨道的升交点赤经。计算公式如下：

式中：T为该天的简略儒略日与2000年1月1日12：00的简略儒略日之差除以100年即36525天。其计算代码如下：

%--------------------------------------------------------------------------
%
% NutAngles: Nutation in longitude and obliquity
%
% Input:
%   Mjd_TT    Modified Julian Date (Terrestrial Time)
%
% Output:
%   dpsi, deps
%
% Last modified:   2018/01/04   M. Mahooti
%
%--------------------------------------------------------------------------
function [dpsi, deps] = NutAngles(Mjd_TT)

global const

T  = (Mjd_TT-const.MJD_J2000)/36525;
T2 = T*T;
T3 = T2*T;
rev = 360*3600;  % arcsec/revolution

N_coeff = 106;
C = [
  % l  l' F  D Om    dpsi    *T     deps     *T
    0, 0, 0, 0, 1,-1719960,-1742,  920250,   89 ;   %   1
    0, 0, 0, 0, 2,   20620,    2,   -8950,    5 ;   %   2
   -2, 0, 2, 0, 1,     460,    0,    -240,    0 ;   %   3
    2, 0,-2, 0, 0,     110,    0,       0,    0 ;   %   4
   -2, 0, 2, 0, 2,     -30,    0,      10,    0 ;   %   5
    1,-1, 0,-1, 0,     -30,    0,       0,    0 ;   %   6
    0,-2, 2,-2, 1,     -20,    0,      10,    0 ;   %   7
    2, 0,-2, 0, 1,      10,    0,       0,    0 ;   %   8
    0, 0, 2,-2, 2, -131870,  -16,   57360,  -31 ;   %   9
    0, 1, 0, 0, 0,   14260,  -34,     540,   -1 ;   %  10
    0, 1, 2,-2, 2,   -5170,   12,    2240,   -6 ;   %  11
    0,-1, 2,-2, 2,    2170,   -5,    -950,    3 ;   %  12
    0, 0, 2,-2, 1,    1290,    1,    -700,    0 ;   %  13
    2, 0, 0,-2, 0,     480,    0,      10,    0 ;   %  14
    0, 0, 2,-2, 0,    -220,    0,       0,    0 ;   %  15
    0, 2, 0, 0, 0,     170,   -1,       0,    0 ;   %  16
    0, 1, 0, 0, 1,    -150,    0,      90,    0 ;   %  17
    0, 2, 2,-2, 2,    -160,    1,      70,    0 ;   %  18
    0,-1, 0, 0, 1,    -120,    0,      60,    0 ;   %  19
   -2, 0, 0, 2, 1,     -60,    0,      30,    0 ;   %  20
    0,-1, 2,-2, 1,     -50,    0,      30,    0 ;   %  21
    2, 0, 0,-2, 1,      40,    0,     -20,    0 ;   %  22
    0, 1, 2,-2, 1,      40,    0,     -20,    0 ;   %  23
    1, 0, 0,-1, 0,     -40,    0,       0,    0 ;   %  24
    2, 1, 0,-2, 0,      10,    0,       0,    0 ;   %  25
    0, 0,-2, 2, 1,      10,    0,       0,    0 ;   %  26
    0, 1,-2, 2, 0,     -10,    0,       0,    0 ;   %  27
    0, 1, 0, 0, 2,      10,    0,       0,    0 ;   %  28
   -1, 0, 0, 1, 1,      10,    0,       0,    0 ;   %  29
    0, 1, 2,-2, 0,     -10,    0,       0,    0 ;   %  30
    0, 0, 2, 0, 2,  -22740,   -2,    9770,   -5 ;   %  31
    1, 0, 0, 0, 0,    7120,    1,     -70,    0 ;   %  32
    0, 0, 2, 0, 1,   -3860,   -4,    2000,    0 ;   %  33
    1, 0, 2, 0, 2,   -3010,    0,    1290,   -1 ;   %  34
    1, 0, 0,-2, 0,   -1580,    0,     -10,    0 ;   %  35
   -1, 0, 2, 0, 2,    1230,    0,    -530,    0 ;   %  36
    0, 0, 0, 2, 0,     630,    0,     -20,    0 ;   %  37
    1, 0, 0, 0, 1,     630,    1,    -330,    0 ;   %  38
   -1, 0, 0, 0, 1,    -580,   -1,     320,    0 ;   %  39
   -1, 0, 2, 2, 2,    -590,    0,     260,    0 ;   %  40
    1, 0, 2, 0, 1,    -510,    0,     270,    0 ;   %  41
    0, 0, 2, 2, 2,    -380,    0,     160,    0 ;   %  42
    2, 0, 0, 0, 0,     290,    0,     -10,    0 ;   %  43
    1, 0, 2,-2, 2,     290,    0,    -120,    0 ;   %  44
    2, 0, 2, 0, 2,    -310,    0,     130,    0 ;   %  45
    0, 0, 2, 0, 0,     260,    0,     -10,    0 ;   %  46
   -1, 0, 2, 0, 1,     210,    0,    -100,    0 ;   %  47
   -1, 0, 0, 2, 1,     160,    0,     -80,    0 ;   %  48
    1, 0, 0,-2, 1,    -130,    0,      70,    0 ;   %  49
   -1, 0, 2, 2, 1,    -100,    0,      50,    0 ;   %  50
    1, 1, 0,-2, 0,     -70,    0,       0,    0 ;   %  51
    0, 1, 2, 0, 2,      70,    0,     -30,    0 ;   %  52
    0,-1, 2, 0, 2,     -70,    0,      30,    0 ;   %  53
    1, 0, 2, 2, 2,     -80,    0,      30,    0 ;   %  54
    1, 0, 0, 2, 0,      60,    0,       0,    0 ;   %  55
    2, 0, 2,-2, 2,      60,    0,     -30,    0 ;   %  56
    0, 0, 0, 2, 1,     -60,    0,      30,    0 ;   %  57
    0, 0, 2, 2, 1,     -70,    0,      30,    0 ;   %  58
    1, 0, 2,-2, 1,      60,    0,     -30,    0 ;   %  59
    0, 0, 0,-2, 1,     -50,    0,      30,    0 ;   %  60
    1,-1, 0, 0, 0,      50,    0,       0,    0 ;   %  61
    2, 0, 2, 0, 1,     -50,    0,      30,    0 ;   %  62
    0, 1, 0,-2, 0,     -40,    0,       0,    0 ;   %  63
    1, 0,-2, 0, 0,      40,    0,       0,    0 ;   %  64
    0, 0, 0, 1, 0,     -40,    0,       0,    0 ;   %  65
    1, 1, 0, 0, 0,     -30,    0,       0,    0 ;   %  66
    1, 0, 2, 0, 0,      30,    0,       0,    0 ;   %  67
    1,-1, 2, 0, 2,     -30,    0,      10,    0 ;   %  68
   -1,-1, 2, 2, 2,     -30,    0,      10,    0 ;   %  69
   -2, 0, 0, 0, 1,     -20,    0,      10,    0 ;   %  70
    3, 0, 2, 0, 2,     -30,    0,      10,    0 ;   %  71
    0,-1, 2, 2, 2,     -30,    0,      10,    0 ;   %  72
    1, 1, 2, 0, 2,      20,    0,     -10,    0 ;   %  73
   -1, 0, 2,-2, 1,     -20,    0,      10,    0 ;   %  74
    2, 0, 0, 0, 1,      20,    0,     -10,    0 ;   %  75
    1, 0, 0, 0, 2,     -20,    0,      10,    0 ;   %  76
    3, 0, 0, 0, 0,      20,    0,       0,    0 ;   %  77
    0, 0, 2, 1, 2,      20,    0,     -10,    0 ;   %  78
   -1, 0, 0, 0, 2,      10,    0,     -10,    0 ;   %  79
    1, 0, 0,-4, 0,     -10,    0,       0,    0 ;   %  80
   -2, 0, 2, 2, 2,      10,    0,     -10,    0 ;   %  81
   -1, 0, 2, 4, 2,     -20,    0,      10,    0 ;   %  82
    2, 0, 0,-4, 0,     -10,    0,       0,    0 ;   %  83
    1, 1, 2,-2, 2,      10,    0,     -10,    0 ;   %  84
    1, 0, 2, 2, 1,     -10,    0,      10,    0 ;   %  85
   -2, 0, 2, 4, 2,     -10,    0,      10,    0 ;   %  86
   -1, 0, 4, 0, 2,      10,    0,       0,    0 ;   %  87
    1,-1, 0,-2, 0,      10,    0,       0,    0 ;   %  88
    2, 0, 2,-2, 1,      10,    0,     -10,    0 ;   %  89
    2, 0, 2, 2, 2,     -10,    0,       0,    0 ;   %  90
    1, 0, 0, 2, 1,     -10,    0,       0,    0 ;   %  91
    0, 0, 4,-2, 2,      10,    0,       0,    0 ;   %  92
    3, 0, 2,-2, 2,      10,    0,       0,    0 ;   %  93
    1, 0, 2,-2, 0,     -10,    0,       0,    0 ;   %  94
    0, 1, 2, 0, 1,      10,    0,       0,    0 ;   %  95
   -1,-1, 0, 2, 1,      10,    0,       0,    0 ;   %  96
    0, 0,-2, 0, 1,     -10,    0,       0,    0 ;   %  97
    0, 0, 2,-1, 2,     -10,    0,       0,    0 ;   %  98
    0, 1, 0, 2, 0,     -10,    0,       0,    0 ;   %  99
    1, 0,-2,-2, 0,     -10,    0,       0,    0 ;   % 100
    0,-1, 2, 0, 1,     -10,    0,       0,    0 ;   % 101
    1, 1, 0,-2, 1,     -10,    0,       0,    0 ;   % 102
    1, 0,-2, 2, 0,     -10,    0,       0,    0 ;   % 103
    2, 0, 0, 2, 0,      10,    0,       0,    0 ;   % 104
    0, 0, 2, 4, 2,     -10,    0,       0,    0 ;   % 105
    0, 1, 0, 1, 0,      10,    0,       0,    0     % 106
 ];

% Mean arguments of luni-solar motion

%   l   mean anomaly of the Moon
%   l'  mean anomaly of the Sun
%   F   mean argument of latitude
%   D   mean longitude elongation of the Moon from the Sun 
%   Om  mean longitude of the ascending node  
l  = mod (  485866.733 + (1325.0*rev +  715922.633)*T ...   
                               + 31.310*T2 + 0.064*T3, rev );
lp = mod ( 1287099.804 + (  99.0*rev + 1292581.224)*T ...            
                               -  0.577*T2 - 0.012*T3, rev );
F  = mod (  335778.877 + (1342.0*rev +  295263.137)*T ...   
                               - 13.257*T2 + 0.011*T3, rev );
D  = mod ( 1072261.307 + (1236.0*rev + 1105601.328)*T ...   
                               -  6.891*T2 + 0.019*T3, rev );
Om = mod (  450160.280 - (   5.0*rev +  482890.539)*T ...   
                               +  7.455*T2 + 0.008*T3, rev );

% Nutation in longitude and obliquity [rad]  
dpsi = 0;
deps = 0;

for i=1:N_coeff
  arg  =  ( C(i,1)*l+C(i,2)*lp+C(i,3)*F+C(i,4)*D+C(i,5)*Om )/const.Arcs;
  dpsi = dpsi + ( C(i,6)+C(i,7)*T ) * sin(arg);
  deps = deps + ( C(i,8)+C(i,9)*T ) * cos(arg);
end
    
dpsi = 1e-5 * dpsi/const.Arcs;
deps = 1e-5 * deps/const.Arcs;

格林尼治真地方时（GAST）与格林尼治平地方时（GMST）的关系为

$\textup{GAST}=\textup{GMST}+d\psi\cos{\varepsilon }$

%--------------------------------------------------------------------------
%
% gast: Greenwich Apparent Sidereal Time
%
% Input:
%   Mjd_UT1   Modified Julian Date UT1
%
% Output:
%   gstime    GAST in [rad]
%
% Last modified:   2018/01/04   M. Mahooti
%
%--------------------------------------------------------------------------
function gstime = gast(Mjd_UT1)

global const

gstime = mod(gmst(Mjd_UT1) + EqnEquinox(Mjd_UT1), const.pi2);

%--------------------------------------------------------------------------
%
% EqnEquinox: Computation of the equation of the equinoxes
%
% Input:
%   Mjd_TT    Modified Julian Date (Terrestrial Time)
% 
% Output:
%    EqE      Equation of the equinoxes
%
% Notes:
%   The equation of the equinoxes dpsi*cos(eps) is the right ascension of
%   the mean equinox referred to the true equator and equinox and is equal
%   to the difference between apparent and mean sidereal time.
%
% Last modified:   2015/08/12   M. Mahooti
% 
%--------------------------------------------------------------------------
function EqE = EqnEquinox (Mjd_TT)

% Nutation in longitude and obliquity 
[dpsi, deps] = NutAngles (Mjd_TT);

% Equation of the equinoxes
EqE = dpsi * cos(MeanObliquity(Mjd_TT));

由于行星的进动，黄道的倾角略有下降，相当于

代码如下：

%--------------------------------------------------------------------------
%
% MeanObliquity: Computes the mean obliquity of the ecliptic
%
% Input:
%  Mjd_TT    Modified Julian Date (Terrestrial Time)
% 
% Output:
%  MOblq     Mean obliquity of the ecliptic
%
% Last modified:   2018/01/04   M. Mahooti
% 
%--------------------------------------------------------------------------
function MOblq = MeanObliquity(Mjd_TT)

global const

T = (Mjd_TT-const.MJD_J2000)/36525;

MOblq = const.Rad*(23.43929111-(46.8150+(0.00059-0.001813*T)*T)*T/3600);

最后得到了真实日期的坐标系（采用的进动章动理论定义）与地球赤道和格林尼治子午线对齐的系统之间的转换矩阵

$\boldsymbol{\Theta} (t)=R_z(\textup{GAST})$

代码如下：

%--------------------------------------------------------------------------
%
% GHAMatrix: Transformation from true equator and equinox to Earth equator
%            and Greenwich meridian system 
%
% Input:
%   Mjd_UT1   Modified Julian Date UT1
% 
% Output:
%   GHAmat    Greenwich Hour Angle matrix
%
% Last modified:   2015/08/12   M. Mahooti
%
%--------------------------------------------------------------------------
function GHAmat = GHAMatrix(Mjd_UT1)

GHAmat = R_z(gast(Mjd_UT1));

%--------------------------------------------------------------------------
%  Input:
%    angle       angle of rotation [rad]
%
%  Output:
%    rotmat      rotation matrix
%--------------------------------------------------------------------------
function [rotmat] = R_z(angle)

C = cos(angle);
S = sin(angle);
rotmat = zeros(3,3);

rotmat(1,1) =      C;  rotmat(1,2) =   S;  rotmat(1,3) = 0.0;
rotmat(2,1) = -1.0*S;  rotmat(2,2) =   C;  rotmat(2,3) = 0.0;
rotmat(3,1) =    0.0;  rotmat(3,2) = 0.0;  rotmat(3,3) = 1.0;

矢志航天的阿洪

关注

12
点赞
踩
21

收藏

觉得还不错? 一键收藏
0
评论
最小二乘法的状态转移方程（2）

通常给出的地球指向参数文本格式为：1962 01 01：公历日期37665 ：儒略日-0.0127 ：x极移，单位：角秒0.213 ：y极移，单位：角秒0.0326338：世界时UT1与协调世界时UTC的差0.0017230：LOD 地球日长变化单位：s0.064261 ：x向的章动值单位：角秒0.006067 ：y向的章动值单位：角秒0.000000 ：x向极移的变化量单位：角秒。
复制链接

扫一扫