ITRS/GCRS/J2000坐标系的相互转换
李云飞 (blitheli@gmail.com)
南京大学天文系/上海航天技术研究院
2008-12-6
本文主要阐述了目前国际天文界最新规定的岁差章动模型(IAU 2000A/B),并根据此模型给出
ITRS,GCRS 和 J2000 平赤道地心系相互转换的详细步骤。
本文主要依据《IERS Conventions 2003》(IERS Technical Note No. 32)而来,由于文章为英文,
且详细叙述了 ITRS,GCRS 参考系以及岁差章动等模型,会使得读者初次阅读(或天文背景知识不够)
产生阅读障碍。笔者也是在多次反复阅读的基础上并参考其它书籍和文献才稍微弄清楚,因此本文
是笔者关于此文献的一个总结,希望能够对读者有所帮助。
本文不打算对坐标系转换的理论(岁差章动,极移等)进行详细讲解,仅仅给出坐标转换的基
本步骤和必要常识,以便读者根据此文能够迅速掌握 ITRS,GCRS 和 J2000 平赤道地心系相互转换的
基本原理,并能够依据 IERS 提供的 Fortran 源程序进行实际编程应用。
0 前言
0.1 名词缩写和解释
barycentric celestial reference system
celestial intermediate origin
celestial intermediate pole
celestial intermediate reference system
earth orientation parameters
geocentric celestial reference system
Greenwich mean sidereal time
Greenwich apparent sidereal time
international astronomical union
international celestial reference system
international Earth rotation and reference systems service
international terrestrial reference system
2000 January 1.5
standards of fundamental astronomy
terrestrial intermediate origin
terrestrial intermediate reference system
terrestrial time
universal time UT1
coordinated universal time
milliarcsecond (
BCRS
CEO
CIP
CIRS
EOP
GCRS
GMST
GAST
IAU
ICRS
IERS
ITRS
J2000
SOFA
TEO
TIRS
TT
UT
UTC
MAS
ITRS
就是我们常说的地固坐标系,其原点在地球质心(包含大气海洋等质量),坐标系 xy 平面为地
0.001
''
)
0.2
球赤道面,z 轴指向北极 CIO 处,x 轴指向格林威治子午线与赤道面交点处。
此坐标系固定在地球上,地面站测控,以及地球引力场系数等都在此坐标系下定义。
0.3 J2000.0 Mean equator, mean equinox dynamical system
常被称为 J2000 平赤道地心坐标系。其原点也是在地球质心,xy 平面为 J2000 时刻的地球平赤
道面,x 轴指向 J2000 时刻的平春分点(J2000 时刻平赤道面与平黄道面的一个交点)。
此坐标系常被作为地球卫星的惯性坐标系,卫星运动积分等都在此坐标系计算。
0.4 GCRS
J2000 地心天球坐标系,其定义与 J2000 平赤道地心坐标系仅有一个常值偏差矩阵 B。目前 IAU
推荐用此坐标系逐渐取代 J2000 平赤道地心坐标系。
0.5 一些时间量
1
(
−
UT UTC UT UTC
1
+
DAT TAI UTC
−
TAI UTC DAT
+
TT TAI
32.184
=
=
=
=
+
s
)
(1)
给定 UTC 时刻,可求得其距离 J2000 的世纪数,即
t = (TT - 2000 January 1d 12h TT) in days/36525.
上式中 2000 January 1d 12h TT 对应的儒略日(Julian day)为 2451545.0 天。
0.6 坐标转换
当一坐标系绕其 3 个轴旋转θ角时,则坐标旋转矩阵可表述为:
R
x
( )
θ
R
y
( )
θ
R
z
( )
θ
0
cos
θ
sin
θ
−
0
θ
1
0
1
0
0
cos
0
sin
θ
cos
θ
sin
θ
0
⎛
⎜
= ⎜
⎜
⎝
⎛
⎜
= ⎜
⎜
⎝
⎛
⎜
= −⎜
⎜
⎝
sin
cos
0
0
⎞
⎟
sin
θ
⎟
⎟
cos
θ
⎠
sin
θ
−
⎞
⎟
0
⎟
⎟
cos
⎠
θ
⎞
⎟
θ
⎟
⎟
⎠
θ
0
0
1
其旋转方向符合右手螺旋法则,即逆时针旋转为正方向。另外坐标旋转矩阵具备如下性质:
R
1( )
−
θ
=
T
R
( )
θ
=
R
(
)
θ
−
1 背景知识
我们知道,地球的自转轴在惯性空间中不是固定的,而是不断摆动的。此摆动造成地轴绕北黄
级顺时针运动,夹角约为 23.5 度。于此同时,地轴还在做微小的抖动,见下图。前者的运动称为岁
差(Precession),后者运动成为章动(Nutation)。岁差章动的原因主要有两个方面。其一是太阳系行星
对地球绕日轨道所产生的摄动影响;其二是太阳和月球对地球赤道隆起部分的摄动影响。
关于岁差章动的计算,此前一直采用 IAU 1976 岁差模型和 1980 章动模型。随着时间的推移,
此模型的精度逐渐跟不上需要。因此,IAU 又规定,自 2003 年 1 月 1 日起,采用新的岁差章动模型,
即 IAU 2000A 模型(精度到达 0.2mas),或者 IAU 2000B 模型(精度达到 1mas)。
岁差章动模型描述了地球自转轴的运动,另外,由于地球表面的海洋,大气运动以及地核内部
液体的运动造成地球自转轴并不是相对地球不动的;相对地球北极 CIO 点来说有个小范围的运动,
此种现象成为极移。
此外,地球的自转也不是均匀的,也很复杂。
综上所述,从地固坐标系(ITRS)到地心惯性系(GCRS 或 J2000 平赤道地心系)的坐标转换矩阵由
极移,自转和岁差章动组成。
23.5
0
图表 1 地轴的岁差章动现象
2
ITRS 到 GCRS 的转换矩阵
本文采用 IAU 2000A/B 岁差章动模型,在某历元 UTC 时刻,ITRS 到 GCRS 的转换矩阵可写成:
r
GCRS
=
Q t R t W t
( )
( )
( )
⋅
⋅
⋅
=
T
HG t
( )
⋅
(2)
r
ITRS
r
ITRS
分别对应同一向量在 ITRS 和 GCRS 坐标系中的坐标。
( )R t 和
分别对应极移,自转和岁差章动转换矩阵。
( )Q t
其中,
r
和
r
ITRS
GCRS
,( )W t
上式中,
在计算 ( )R t 和
( )Q t
的时候,会有两种计算方法,我们分别称为 CEO-based 转换方法和
equinox-based 转换方法。其中前者为 IAU 提出的新的计算方法。
下面分别给出上述三个转换矩阵的求解过程。
2.1 极移矩阵
( )W t
W t
( )
=
R
z
(
−
上式中, 为:
's
s R x
')
(
⋅
y
)
⋅
R yp
x
(
)
p
(3)
极移量(
x
,
p
)
y 的求解为:
p
x
,
(
p
s
'
=
-0.047
mas t
⋅
y
p
)
=
x y
( ,
)
IERS
y
+ Δ Δ
x
(
,
)
tidal
y
+ Δ Δ
x
(
,
)
nutation
极移量主要是由 IERS 根据天文观测给出的
(
球潮汐和章动的影响,会对极移有微小的修正
,每周都有新的观测数据,此外,由于地
和 (
由 IERS 给出的观测数据计算求得,(
y
x
,
Δ Δ
x
y
,
)
Δ Δ
x y
( ,
)IERS
x
y
,
)tidal
Δ Δ
y
+ Δ Δ
x y
上式中,( ,
)IERS
)nutation
(
tidal
可由公式计
。
x
,
nutation
)
算得到,IERS 提供此 fortran 源程序。
2.2 地球自转矩阵 ( )R t
R t
( )
=
R θ
)
−
z
(
(4)
地球自转角θ的求解根据转换方法的不同有不同的求解方式(CEO-based 或者 equinox-based),
具体求解 IERS 给出了 fortran 源程序。
2.3 岁差章动矩阵
( )Q t
前面提到,计算此矩阵有两种方法:
1.) CEO-based 方法:
这种方法是 IAU 最新提出的并极力倡导的,其计算公式为:
上式中:
2
Q t
( )
aXY
aX
1
⎛
−
−
⎜
aXY
aY
1
2
−
= −
⎜
⎜ −
Y
X
−
⎝
a
X
1 2 1 8(
=
其中:
+
2
+
X
Y
a X
(
2
+
Y
2
)
⎞
⎟
⎟
⎟
⎠
⋅
R s
( )
z
1
−
Y
)
2
(
X Y
,
)
=
(
X Y
,
)
IAU
2000
+
(
dX dY
,
)
IERS
(5)
(6)
s
(
X Y
,
)IAU
2000
和 可根据 IAU 2000A/B 岁差章动模型求解出,IERS 同样给出求解的 fortran 源程
序,另外,由于 IAU2000A/B 岁差章动模型没有包含地轴的高频率运动,所以要加上 IERS 通过观测
(
数据给出的高频率修正项
2.) Equinox-based 方法:
dX dY
)IERS
。
,
Q t
( )
= ⋅
B P t N t
( )
( )
⋅
( )P t
其中,常值偏差矩阵 B ,岁差矩阵
( )N t
)
(
(
)
⋅
ξ
−
δα
−
η
0
0
0
R
R
)
(
)
(
(
(
ε ψ
ω
⋅
−
A
x
z
A
0
R
(
)
(
)
(
ε
ψ ε ε
+ Δ
Δ
−
x
和章动矩阵
R
R
)
⋅
x
y
R
R
)
⋅
⋅
x
z
R
R
)
⋅
⋅
z
x
B R
=
z
P t
( )
=
N t
( )
=
(
⎧
⎪
⎨
⎪
⎩
如下:
χ
−
A
)
章动量 (
)ψ ε
Δ
,
Δ 为:
(
)
ψ ε ψ ε
Δ
(
Δ = Δ
Δ
)
,
,
IAU
2000
+ Δ
(
)
δ ψδ ε
Δ
,
IERS
(7)
(8)
(9)
(
,
Δ
)IAU
2000
ψ ε
Δ
上式中,
由 IERS 2000A 章动模型给出。前面提到过,IAU 2000A/B 模型提供
的岁差章动不包含高频率项,而是由 IERS 的观测数据提供(上式右端最后一项),但是在 IERS 给
转
出的观测数据中仅仅给出
,我们可以通过 IERS 提供的 fortran 源程序将
dX dY
,dX
dY
(
(
,
)IERS
)IERS
换为
)IERS
Δε 。
(
,δ ψδΔ
其余参数皆为岁差参数,可以通过公式求出,此处从略。
值得一提的是常值偏差矩阵 B 中的参数也是给定的,在 CEO-based 方法求解中,此偏差是包含
X Y
,
中的。
(
)IAU
在
3 利用 IERS 提供的 Fortran 源程序进行转换
2000
上节中详细讲述了 ITRS 到 GCRS 转换矩阵的求解过程,在实际应用中,如果是自己编写源程
序的话是件非常琐碎的事情,因为 IERS 2000A/B 章动模型的参数多达 1000 多项。幸而这些基本子
程序 IERS 都提供了,我们所做的就是如何正确的运用这些源程序,并将它们组合起来。
在ftp://maia.usno.navy.mil/conv2000/chapter5/上可下载相关的源程序,子程序列表如下:
图表 2 IERS 提供有关 IAU 2000 的 fortran 源程序
子程序名
BPN2000
CBPN2000
EE2000
EECT2000
ERA2000
GMST2000
GST2000
NU2000A
NU2000B
POM2000
SP2000
T2C2000
XYS2000A
interp.f
uai2000.f
说明
CEO-based intermediate-to-celestial matrix
equinox-based true-to-celestial matrix
equation of the equinoxes (EE)
EE complementary terms
Earth Rotation Angle
Greenwich Mean Sidereal Time
Greenwich (apparent) Sidereal Time
nutation, IAU 2000A
nutation, IAU 2000B
form polar-motion matrix
the quantity s’
form terrestrial to celestial matrix
X, Y, s
Interpolation of IERS polar motion and UT1 time series
IAU 2000 celestial pole offsets conversion (dpsi,deps,dX,dY)
3.1
上表最后可从 ftp://hpiers.obspm.fr/iers/models 上下载得到。
IERS 观测数据的处理
前面一再提到 IAU 2000A/B 章动岁差模型不包含高频率项,因此在完整的坐标转换过程中,必
须考虑到 IERS 提供观测数据的高频率修正项。
IERS 每周发布 Bulletins A,每月发布 Bulletins B,它们都是描述 EOP 的参数,下面是综合的
EOP 参数文件(C 04),其主要内容如下:
注意上述文件中,dPsi,dEpsilon 是针对 IAU 1980 章动模型的修正项
图表 3 IERS 观测数据 C 04 文件
(
δ ψδ ε
Δ
Δ
,
)IAU
,不是 IAU
1980
,
(
δ ψδ ε
Δ
Δ
,
)IAU
y
p
和
,p
)IERS
(UT1-UTC)
2000A 章动模型的修正项
,因此不需要使用。
MJD x y
, UT1-UTC
,
和
)
2000
根据 UTC 时刻的儒略日,加上数据列表
x
(
)
T1-UTC
子程序插值计算出对应 UTC 时刻的
x y
(U
IERS 的 观 测 数 据 ( ,
PM_GRAVI 计算由潮汐和章动引起的高频率修正项
UT1-UTC
然后分别相加,给出最后的
MJD dX dY LOD DAT ,可插值计算出
,可调用 interp.f 文件中的 interp
。在子程序 interp 中,先插值计算出
, 然 后 内 部 调 用 子 程 序 PMUT1_OCEANS 和
(
,
。
,
y 和
)p
另外,根据 IERS 的观测数据列表
,LOD
和 DAT。若为 equinox-based 方法转换,则需要调用 uai2000.f 文件中的子程序 dXdY_dpsideps 将
dX dY
(
,
有了
(1)求得UT1,TT和t 。这些时间量在以后的子程序中都需
(
δ ψδ ε
Δ
和DAT,则可根据式
)IERS
转换为
UT1-UTC
(UT1-UTC)
dX dY
)nutation
x
y
Δ Δ
x
y
Δ Δ
)IERS
)IERS
。
px
IERS
IERS
,(
tidal
tidal
,
Δ
(
,
,
)
(
,
和
,
,
,
,
要。
3.2 具体转换步骤
首先调用子程序 SP2000 求得 ,再由上面插值求得的 (
's
x
,
y ,调用子程序 POM2000 即可求
p
)
p
得极移矩阵
( )W t
。
求地球自转和岁差章动矩阵有两种方法,下面分别叙述:
1.) CEO-based transformation
调用 ERA2000 求得地球自转角θ;
X Y
,
然后调用子程序XYS2000A求得 (
,
nox-based transformation
)X Y ,参见式(6)。根据 (,
则可求得(
2.) E
qui
调用子程序NU2000A求得章动
后的 (
由 εΔ
)IERS
δ ψδ ε
Δ
,调用 GST20
转换
Δ
,
(
量
)IAU
,
ψ ε
Δ
Δ
)ψ ε
,得到最后的 (
,
Δ ,参见式(9);
Δ
自转角θ;
00 即可求得地球
2000
)IAU
,
, )X Y s ,利用子程序BPN2000 即可求得岁差章动转换矩阵 ( )Q t 。
s 再加上上面观测数据插值的 (
dX dY
)IERS
和 ,
2000
,
,再加上由观测数据插值求得的
(
dX dY
,
)IERS
然后再由 (
)ψ ε
Δ
Δ
,
调用子程序CBPN2000 ,求得岁差章动矩阵 Q t 。此处需要对子程序
(8)进行常值偏差矩阵B和岁差章动矩阵P,N的计算,最后
( )
CBPN2000 进行简单的说明,其内部利用式
给出矩阵 ( )Q
t
。
由上述两种方法之一求得
,( )W t θ和
( )Q t
,调用子程序 T2C2000 即可求得 ITRS 到 GCRS 的
转换矩阵
THG t
( )
。
若采用方法二时,可以用子程序 NU2000B 替代 NU2000A,其它都不变,此种转换的精度稍低
(1mas),但是其计算速度会快很多,在精度要求不是很高的情况下采用此种方法可使计算速度大大
提高。
4 其它的一些说明
1.) 下面给出 GCRS 和 ITRS(包含速度)转换的完整表述:
⎧
r
ITRS
⎪⎪
r
⎨
TIRS
⎪
V
⎪⎩
ITRS
r
GCRS
⋅
Q t R t W t
( )
[
( )
( )]
T
=
⋅
⋅
Q t R t
r
[
( )]
( )
T
=
⋅
⋅
GCRS
{
R t Q t V
W t
( )
( )
( )
T
T
T
=
⋅
GCRS
⋅
=
HG t
( )
⋅
r
GCRS
rω⊕
×
TIRS
−
}
⎧
r
GCRS
⎪⎪
r
⎨
TIRS
⎪
V
⎪⎩
GCRS
r
ITRS
=
⋅
⋅
( )
Q t R t W t
( )
( )
⋅
=
r
W t
( )
=
⋅
ITRS
Q t R t W t V
( )
( )
⋅
⋅
=
ITRS
( )
⋅
{
T
HG t
( )
⋅
r
ITRS
rω⊕
×
TIRS
+
}
(10)
(11)
和
其中:
ω
⊕
=
7.292115146706979 10
×
5
−
⎧
1
−⎨
⎩
LOD
86400
⎫
⎬
⎭
上式中的 LOD 由 IERS 的观测数据插值求得。
注意上述公式中
的区别,
和
r
TIRSr
TIRSr
经过极移转换矩阵后的坐标。
2.) 若需要GCRS和J2000 平赤道地心系的相互转换,只需要常值偏差矩阵B即可(式(8))。其具体的
为地固系 ITRS 坐标 ITRS
ITRS
r
源程序需要读者自行编制。
r
⎧
GCRS
⎪
V
⎪
GCRS
⎨
r
⎪
J
⎪
V
⎩
J
2000
2000
B r
= ⋅
J
2000
B V
= ⋅
J
B r
T
⋅
=
GCRS
B V
T
⋅
=
GCRS
2000
(12)
3.) 在 IERS 提供的 fortran 源程序中,有部分子程序需要调用 IAU SOFA 软件包中的子程序,有关
IAU SOFA 软件包的说明和使用请参见我的文档《IAU SOFA 软件包介绍》。
4.) IAU SOFA 软件包中也包含有关岁差章动和极移等基本子程序,其主要内容和本文介绍的 fortran
源程序大同小异。但其软件包中还包含以前的岁差章动模型以及最新的 IAU2006 岁差模型,读
者需要的话可参考其文档说明“sofa_pn.pdf”,并有具体例子,强烈建议读者阅读。
5 参考文献
1.
IERS Conventions(2003)
可从IERS网站上下载(http://www.IERS.org)。此文详细叙述了IAU 2003A/B岁差章动模型,以及
ITRS和GCRS坐标系的定义和详细转换过程,也是本文档的主要英文依据。在尽可能的情况下,读
者可以多阅读几遍。
2. 《Fundamentals of Astrodynnamics and Applications》, Third Edition. Microcosm Press, 2007
此书为 David A. Vallado 所著,此书中第 3 章“Coordinate and Systems”详细介绍了时间,坐标
系系统,以及详细的 ITRS 到 GCRS 的转换过程。
3.
ftp://maia.usno.navy.mil/conv2000/chapter5/和ftp://hpiers.obspm.fr/iers/models
此两 ftp 网站上含有 IAU 2000 岁差章动模型的所有 Fortran 源程序。也是本文中子程序的来源处。
4. http://www.iau-sofa.rl.ac.uk/
从此网站上可下载 IAU SOFA 软件包,里面同样包含 IAU 2000 岁差章动模型的所有 Fortran 源
程序,另外还包括最新的 IAU 2006 模型,以及儒略日计算,行星历表等常用基本子程序。
5. http://hpiers.obspm.fr/eop-pc/
此网站包含 EOP 各种类型数据,包含 IERS 的 Bulletin A/B 和 C 04,以及一些其它文件的说明。