就地球卫星轨道力学问题,主要涉及地心天球坐标系(GCRS)和国际地球参考系(ITRS),目前仍符合 WGS84 地球引力场系统,考虑沿用的习惯,不必拘泥于学术上的文字定义,在没有特殊需要时,本书仍引用原来的地固坐标系名称。另外,关于地心天球坐标系,一般情况下可忽略参考架偏差,直接引用J2000.0平赤道坐标系作为地心天球坐标系,这意味着将式(1.6)中的参考架偏差矩阵B作为单位阵处理。
分别用 和 表示卫星在地心天球坐标系(J2000.0地心平赤道坐标系)O-xyz和地固坐标系O-XYZ中的同一位置矢量。卫星的位置矢量在这两个坐标系之间的转换关系为
其中坐标转换矩阵(HG)包含四个旋转矩阵,有
这里的(PR),(NR),(ER)和(EP)是岁差、章动、地球自转和极移矩阵,分别由下列各式表达,即
式(1.23)中的x p ,y p 是极移分量;式(1.24)中的格林尼治恒星时S G 由式(1.27)计算得
这里的μ和Δμ是赤经岁差和章动,J2000.0系统中的格林尼治平恒星时 由式(1.28)计算得
式中引数t是 UT1 时间,但计算其他天文量(岁差章动等)时,该引数t则为 TT时间。
式(1.26)中的岁差常数ζ A ,z A 和θ A 的计算公式为:
θ A 是赤纬岁差,相应的赤经岁差μ(或记作m A )为
式(1.25)中的ε是平黄赤交角。IAU(1980)章动序列给出的黄经章动ΔΨ和交角章动Δε的计算公式,包括振幅大于0″.0001的106项。考虑到一般问题涉及的轨道精度要求,只要取振幅大于0″.005的前 20 项(按大小排列)即可,由于是周期项(最快的是月球运动周期项),没有累积效应,故小于0″.005的项引起的误差只相当于地面定位误差,为米级,对于时间而言的差别小于0 S .001。取前 20 项的公式为
相应的赤经和赤纬章动Δμ和Δθ为
其中平黄赤交角的计算公式为
式(1.32)中,章动序列涉及的与太阳和月球位置有关的5个基本幅角α i (i=1,…,5)的计算公式为
其中1 r =360°,章动序列前20项的有关系数如表1.2所示。如果按前面所说的米级精度考虑,式(1.32)右端的A 1j 和B 1j 除表1.2中列出的A 11 和B 11 外也可略去,但在具体工作中的取项多少,不仅需要符合精度要求,还应考虑相应软件的功能和适应性。上述计算章动量的式(1.32)~式(1.35)中出现的t,即式(1.29)中所定义的世纪数T,但对应的是TT时间,即
表1.2 IAU(1980)章动序列的前20项
续表
上述内容中涉及的各旋转矩阵R x (θ),R y (θ),R z (θ)的计算公式为
注意,旋转矩阵R(θ)是正交矩阵,有
该表达形式的上标“T”和“-1”分别表示矩阵的转置和逆。
在IAU2000规范下,地心天球参考系到地球参考系之间的转换过程由式(1.41)表述为
其中[GCRS]和[ITRS]各对应前面IAU1980规范下的地心天球坐标系和地固坐标系。为了表达的连贯性,仍采用同一位置矢量在两个坐标系中的符号 和 来表达变换关系,即
M(t)是岁差、章动矩阵,R(t)是地球自转矩阵,W(t)是极移矩阵。
关于岁差、章动矩阵M(t),基于春分点的转换关系,岁差、章动矩阵可以写为
其中N(t),P(t)和B分别为章动、岁差和参考架偏差矩阵,其中参考架偏差矩阵B已在前面有过说明,它是一个旋转量很小的常数矩阵,在直接引用J2000.0平赤道坐标系作为地心天球坐标系时,就作为单位阵略去,不再考虑,如有特殊需要,考虑它也极其简单,见式(1.6)。于是有
经岁差、章动矩阵M(t)转换后的坐标系即瞬时赤道坐标系,各矩阵的计算方法由下面分别给出。
在第24届国际天文联合会大会(2000年8月,Manchester)上通过决议,从2003年1月1 日起正式采纳IAU2000岁差—章动模型取代IAU1976岁差模型和IAU1980 章动理论。关于岁差章动的计算就是基于 IAU2000 模型,根据不同精度要求,又分为IAU2000A模型和IAU2000B模型两种,前者精度为0.2mas(毫角秒),后者的精度稍低1mas,下面对此作一简要介绍。
关于岁差量的计算
由标准历元J2000.0到计算历元,平赤道坐标系之间转换的三个赤道岁差参数ζ A ,z A ,θ A ,由式(1.45)计算得:
其中t是自标准历元J2000.0起算的计算历元(TT时刻)的儒略世纪数,即
关于章动量的计算
IAU2000章动模型的黄经章动ΔΨ和交角章动Δε的计算公式为
其中ΔΨ p ,Δε p 是行星章动的长周期项,有
式(1.47)中t的意义同前,见式(1.46)。
章动序列中的幅角α i 同样表示成5个基本幅角的线性组合形式,即
式中n i k 是整数,5个与太阳、月亮位置有关的基本幅角F k 由式(1.50)~式(1.54)表达,即
这 5 个与太阳、月球位置有关的基本幅角F k (k=1,…,5),分别为月球的平近点角、太阳的平近点角、月球的平升交点角距、日月平角距和月球轨道升交点平黄经。这5个基本幅角F k (k=1,…,5)就是前面IAU1980章动序列中出现的5个基本幅角α i (i=1,…,5),见式(1.35)。
为了比较,对应IAU1980章动序列的前20项系数,同样列出IAU2000B章动模型前20项系数 的主要部分如表1.3所示。
表1.3 IAU2000B章动序列前20项的主要部分
(1)经典的三次旋转。
其中各旋转角的计算公式见式(1.45)。
(2)岁差矩阵P(t)=(PR)的四次旋转。
其中后三个旋转角的计算公式为
ε 0 的计算下面给出,见式(1.60)。式(1.57)中t的意义同前,见式(1.46)。
其中Δ μ=Δ Ψ cos ε A ,Δ θ=Δ Ψ sin ε A 是赤经和赤纬章动,Δ Ψ是黄经章动,Δ ε是交角章动,ε A 为瞬时平赤道面与黄道面的交角,称为平黄赤交角,计算公式为
ε 0 即历元(J2000.0)平黄赤交角。
这与 IAU1980 规范有差别,涉及时刻t的格林尼治恒星时S G 的计算。在IAU2000规范中引进了中间赤道的概念 [1-3] ,如图1.1所示,图中所画圆周即表示地心参考系中的中间赤道,E为地球质心,γ为春分点。
图1.1 中间赤道示意图
在天球参考系中观察时,中间赤道与 CIO 固结,称为天球中间赤道,地球中间零点TIO 沿赤道逆时针方向运动,周期为1 恒星日。在地球参考系中观察时,中间赤道与TIO 固结,称为地球中间赤道,天球中间零点CIO 以同样周期沿赤道顺时针方向运动。这两种观察所反映的都是地球绕轴的自转运动,CIO 和 TIO 之间的夹角是地球自转角度的度量,称为地球自转角 ERA(Earth Rotating Angle)。格林尼治真恒星时GST与地球自转角ERA有严格区别,在此前提下,矩阵R(t)的计算公式变为
其中GMST,EE分别称为格林尼治平恒星时和二分差(Equation of the Equinoxes),它们的计算公式为
式(1.62)中的θ(UT 1)是地球自转角 ERA,它是世界时UT 1的线型函数,计算公式为
其中d为自标准历元J2000.0起算的世界时儒略日数(对应UT1时刻),即
式(1.63)右端求和项的幅角和振幅如表1.4所示,表中的角α i 与前面式(1.47)中的α i 同类型,其涉及的基本幅角F k (F,D,Ω)见式(1.52)~式(1.54)。
表1.4 α k ,C k 值
与地球自转角不同,格林尼治真恒星时GST 是从真春分点起量的,它与ERA的关系可以写为
其中EO称为零点差,与岁差、章动在中间赤道上的分量有关。根据IAU 2006/2000A岁差、章动模型,零点差的表达式(包含在1975—2025年之间,所有大于0.5μas的项)为
其中t是从J2000.0 起算的儒略世纪数,参数α k 和振幅C k 如表1.4所示。
这里顺便提一句,最早的地球参考框架是国际纬度局(International Latitude Service)根据 1900—1905 五年的观测提出的国际习用原点 CIO(Conventional International Origin)定义了地固坐标系第三轴的平均指向,即原来的地球平均地极指向。该原点的缩写名称CIO已被现在的天球中间零点CIO所占用,请读者注意,不要被这两种称谓所混淆。
经地球自转矩阵R(t)转换后的坐标系即通常所说的准地固坐标系,而再经地球极移矩阵W(t)的转换,才转为地固坐标系。该转换矩阵W(t)与IAU1980规范的转换矩阵(EP)稍有差别,下一小段介绍。
其中x p ,y p 是中间极CIP在地固坐标系ITRS中的两个极移分量,s´称为地球中间零点TIO的定位角,它提供TIO在CIP赤道上的位置,是x p ,y p 的函数,即
s´(t)无法事先获得,但其量级很小,可以用下列线性公式作为其近似,即
其中t的定义同前,即计算时刻距标准历元J2000.0的儒略世纪数。
1.3.3 IAU1980规范与 IAU2000规范之间的对应关系
上述IAU2000新规范与IAU1980规范之间,除岁差章动等有关参数有较小差别外,地球坐标系与天球坐标系的定义与相关符号的采用也与在老规范前提下的习惯用法有所差别,下面厘清两者相关计算公式之间的对应关系。
在IAU2000规范中,天球参考系到地球参考系地固系的转换关系为
其中M(t)是岁差、章动矩阵,R(t)是地球自转矩阵,W(t)是极移矩阵。而在IAU1980规范下的习惯用法,即记卫星在地球参考系和天球参考系的位置矢量分别为 和 ,则上述转换关系式(1.71→)即可类似地写为
从而有如下对应关系
根据上述比较不难看出,关于岁差(按三次旋转考虑)、章动和地球自转(两个规范的格林尼治恒星时S G 稍有差异)三个相关矩阵,除相应的参数稍有微小差别外,转换过程和计算公式完全相同;唯有极移矩阵的表达形式有差别,按IAU2000规范,极移矩阵W(t)应由前面的式(1.68)表达,即
除增加一个微小的转动矩阵R z (-s´)外,剩下的两个转动矩阵R y (x p )R x (y p )却与式(1.77)的转动次序相反。但因为极移量x p ,y p 很小,这样的转动次序改变,引起的差别应该比极移量更小,最终的转换(指地球参考系与天球参考系的转换)结果不会有明显差别,下面将通过具体算例来证实这一点。
比较的算例是采用地球上格林尼治恒星时及相应的坐标转换计算,这一转换过程所涉及的内容基本上涵盖了上述地固坐标系与地心天球坐标系之间转换的全过程。
分别计算世界时2003年1月1日0时,6月1日0时,12月1日0时(UT1,分别对应的简约儒略日为MJD=52640.0,52791.0,52974.0)的格林尼治真恒星时S G 和平恒星时 ,以及赤经章动 Δμ。按前面 IAU1980规范对应的式(1.21)~式(1.28)计算,并考查章动序列的取项精度,结果如表1.5所示。
表1.5 世界时0时对应的格林尼治真恒星时S G 和平恒星时 S G
由比对结果可以看出,章动序列取前20 项,引起的误差不超过0 S .001,相当于米级精度。
分别采用老规范(IAU1980)和新规范(IAU2000)计算进行比较,计算世界时2009年1月1日0时和4月15日0时(UT1,分别对应的简约儒略日为MJD=54832.0,54936.0)的格林尼治真恒星时S G ,平恒星时 和相应的二分差(定义与原赤经章动 Δμ有微小差别,但老规范仍按赤经章动 Δμ计算),以及地心天球坐标系与地固坐标系中位置矢量的转换(1月1日0时对应的是地心天球坐标系到地固坐标系 ,而4月15日0时对应的是地固坐标系到地心天球坐标系 ),结果分别如表1.6和表1.7所示。
从上述两组算例可以看出,关于地球坐标系的 IAU2000 新规范,其有关参数与 IAU1980 老规范相比,并无明显差别,特别是通常所关心的坐标转换,在米级精度前提下,新老规范都可以用,而且对于章动序列只要取到前20项即可。在此处理下,两个规范计算的格林尼治恒星时与标准值(年历所载)之差都不超过0 S .001,相应的位置量转换(包括 IAU2000 新规范的两种岁差计算结果)之差确实都在米级以内。
表1.6 世界时0时对应的格林尼治真恒星时S G 和平恒星时
表1.7 地固坐标系与地心天球坐标系中位置矢量的转换
上述比较结果,可为相关部门在处理卫星轨道资料(特别是涉及历史资料)时提供一个参考信息,即在某些具有多年延续的轨道资料积累中,相应的计算软件系统是否需要跟着IAU2000新规范做较大的改变和如何改变。
1.3.5 地球赤道面的摆动与坐标系的选择问题
对于地球卫星的轨道运动,无论是低轨还是高轨,影响其轨道变化的主要摄动源是地球非球形引力。因此,在考虑地球卫星运动时,无论是理论分析还是具体的定轨或预报,显然是采用地心赤道坐标系,但由于地球自转过程中在外力矩作用下,自转轴在空间摆动(前面提到的岁差章动现象)导致赤道面的变化,从而使地球引力位在确定的地心赤道坐标系中随时间发生改变,这一变化给卫星运动增添附加摄动影响,由此引起与此相关联的地心赤道坐标系的选择问题。早在二十世纪六七十年代,Kozai,Lambeck等人 [7-9] 为了简化上述引力位变化对地球卫星轨道的影响,在建立摄动分析解时,提出了一种混合形式的轨道坐标系O-x´y´z´。其坐标原点仍是地球质心(简称地心),x´-y´坐标面是t时刻的瞬时赤道面,而x´轴方向却是历元(当时取J1950.0)平春分点方向。在这种坐标系中(不考虑极移),地球引力位无变化,其实际影响转化为运动坐标系的“惯性力”影响,但其影响较小,通常可忽略,故这种坐标系在某些领域一直沿用至今。如在空间目标监测中,美国所公开发布的TLE根数即这种坐标系,不过x´轴方向已改为t时刻的平春分点方向。然而在地面测控系统,以及通常的卫星定轨和预报中又宜采用历元平赤道坐标系,即当今的J2000.0地心天球坐标系,特别是高精度的定轨和预报工作,这就导致两种坐标系的互相转换问题,从而引起不必要的麻烦。在标准历元50年改变一次的前提下,作者的一系列工作表明 [10-12] ,上述赤道面在空间变化引起的附加摄动问题可以简单地得到解决,无须采用混合形式的轨道坐标系,在卫星轨道的各类工作中,就采用统一的J2000.0地心天球坐标系。关于这一问题,将在后面第5章中给出全面分析和具体的解决方法,这就使卫星运动在空间坐标系的选择上避免不必要的转换麻烦。