CN109655064B - 一种基于相对转动的地球固连系-惯性系快速高精度转换方法 - Google Patents

Abstract

本发明提出了一种基于相对转动的地球固连系‑惯性系快速高精度转换方法。该方法仅在事先生成星历插值表时，调用SOFA的相关子程序。在生成星历表后，只需要查找表中包含对应待求时刻的区间，读取相关的参数并进行线性插值，最后代入相应的矩阵，完成三次乘法运算即可。该方法只消耗少量的计算、不需要过多的天文理论背景、成本低廉易于实现，在航天的制导导航与控制等方面都具有重要的意义，可用于固连系坐标到惯性系坐标的快速变换。与传统的基于SOFA程序方法相比，这种新型的坐标转换方法大大降低了算法上的复杂度，在保证一定精度的同时能够以很小的计算代价完成相同的功能。

Description

一种基于相对转动的地球固连系-惯性系快速高精度转换 方法

技术领域

本发明提供一种基于相对转动的地球固连系-惯性系快速高精度转换方法，它涉及一种通过计算两个时刻坐标变换矩阵的相对转角，对其间固连系-惯性系坐标变换矩阵进行快速计算的方法，属于导航技术领域。

背景技术

考虑到地球固连参考系并不是惯性参考系，牛顿运动定律不适用于这种参考系，然而在航天器的制导导航与控制的计算中，通常需要在以地球质心为原点的惯性参考系(如天球参考系)中。目前在航天工业中，大部分的卫星应用(如卫星导航、任务规划)需要同时得到航天器在固连参考系和惯性参考系中的坐标，因此就避免不了在地球参考系和天球参考系之间频繁地进行坐标转换。

国际天文学联合会(International Astronomical Union，简称IAU)一直致力于固连系与惯性系之间坐标变换的理论研究(如岁差章动模型等)并定期发布相关参数。IAU在本世纪初提出的基于天球中间点(Celestial Intermediate Origin，简称CIO)或称为无旋转点的理论比原有的基于春分点的转换理论简化了不少，但是对于工程应用还是显得较为复杂。

国际大地测量学与地球物理学联合会(International Union of Geodesy andGeophysics，简称IUGG)提供了极移的计算理论，并定期发布极移的观测数据。

由这两个机构的岁差、章动和极移相关理论，能够建立完整的坐标变换体系。

IAU和IUGG共同组建了国际地球自转服务局(International Earth RotationService，简称IERS)。IERS提供了一套用于计算从固连系到惯性系坐标变换的参考程序SOFA(Standards of Fundamental Astronomy)。但是SOFA程序十分繁杂，共有200多个子程序，在计算调用过程中付出的代价很大，在航天工业领域使用显得极为不便。例如，实现一次较为简单的时间系统转换，仍然需要调用十几个相关的子程序；实现一次从固连系到惯性系的坐标变换，更是需要数十个相互依赖的子程序。用户需要有一定的背景知识才能理解固连系到惯性系坐标变换的理论来使用程序，这就要付出额外的时间代价。而且在坐标系的变换过程中，大多数的中间变量和临时变量对于用户而言是无需了解的。

本专利提出的基于相对转动的地球固连系-惯性系快速高精度转换方法是一种通过数值计算的快速坐标变换方法。这种新型的坐标变换方法仅需要拟合日期和后一天的相关数据即可插值计算得到当天任意时刻的从固连系到惯性系的坐标变换矩阵。

在坐标转换的计算中，需要付出的代价仅仅是简单的线性插值和矩阵乘法。因此，这种方法的效率很高，能够大大缩减程序执行的时间。同时，这种转换方法也不需要用户拥有相关的背景理论知识，在工程应用中能够提高使用的便利性。

综上，基于相对转动的地球固连系-惯性系快速高精度转换方法具有巨大的研究与应用价值，它不依赖于SOFA的繁杂程序，只需通过事先计算好固定时刻的相关参数，即可计算得到在所需时间范围内，任意时刻对应的从地球固连系到地球惯性系的坐标变换矩阵。

发明内容

(一)发明目的：本发明创新性地利用相对转动的思想，结合插值计算方法实现了快速计算从地球固连系到地球惯性系的坐标变换矩阵。与传统的基于SOFA程序计算坐标转换矩阵的方法相比，基于这种新型坐标转换方法的程序复杂度低，计算量较小，也不涉及较深的理论背景，极大地提高了航天工业领域中相关应用的效率，实现成本低廉，具有广泛的应用前景。

(二)技术方案

本发明是一种基于相对转动的地球固连系-惯性系快速高精度转换方法，其步骤如下：

步骤一：准备工作

地球固连参考系定义：

常用的地球固连参考系有国际地球参考系、协议地球参考系、1984世界大地坐标系和2000国家大地坐标系。国际地球参考系(International Terrestrial ReferenceSystem，简称ITRS)是目前最常用的地球参考系，它由IERS发布。IERS一直在对框架进行不断地改进和修正，至今IERS已经公布了13个版本，目前最新的版本是ITRF2014。该参考系的z轴被定义为与地球自转轴平行或重合，x轴则是指向本初子午线与地球赤道的交点。本发明中实现参考系变换时使用的地球参考系是ITRS。

地心惯性参考系定义：

常用的地心惯性参考系有J2000.0参考系、质心天球参考系和地心天球参考系(Geocentric Celestial Reference System，简称GCRS)。GCRS的坐标轴指向由甚长干涉(VLBI)测量的一组河外射电源在J2000.0天球赤道坐标系确定，是一个准空间固定参考系，可以便于表示天体在空间中的方位。J2000.0参考系是用J2000.0时刻的天球赤道与赤经章动来确定的天球参考系，即平春分点平赤道地心参考系，该参考系的坐标原点位于地球质心，其xoy平面是J2000.0时刻的平赤道面，x轴指向J2000.0时刻的平春分点。由于该坐标系只考虑了岁差修正而忽略了章动部分的修正，因此在航天器精密定轨中一般不使用该参考系。本发明中实现参考系变换时使用的地心惯性参考系是GCRS。

基础天文标准库(SOFA)：是由国际天文联合会(IAU)赞助的项目，旨在为天文计算过程中提供具有权威性的有效算法以及部分常数的数值。SOFA程序的第一版代码于2001年10月底发布，在此之后基本是每18个月都会进行版本维护。SOFA的算法最初是用Fotran77编写的，由于C语言的广泛应用及其高效性，在2009年2月IAU又推出了ANSI C版本的SOFA。目前，最新版SOFA(截止2018年10月)主要包含了天体测量学、日历、时间尺度、黄道坐标、地球自转和恒星时间、银河坐标、地心/大地转换、岁差，章动和极移以及星历转换等内容。

两部分儒略日：

本发明中使用的儒略日由两个变量(JD1与JD2)组成。JD1表示的是儒略日的整数部分，JD2表示的是儒略日的小数部分。以两个长双精度的变量来存储日期有直观、精度高和方便计算的优点。

地球旋转角(Earth Rotation Angle，简称ERA)：

地球旋转角ERA₀是t₀时刻天球中间极赤道上CIO与TIO之间的角距离，它是UT1的线性函数。当t₁>t₀时，t₁时刻的地球旋转角ERA₁可以表示为：

上式中，

为t₀至t₁时刻的地球旋转平均角速度。

地球旋转角矩阵[ERA]是一个绕y轴旋转ERA角度的3×3矩阵。

固连系-惯性系转换矩阵：

[GCRS]＝[T][ITRS]＝[W][ERA][CIO][ITRS] (2)

上式中，[T]表示整个坐标变换过程的总体矩阵，[CIO]表示CIO理论中的岁差章动矩阵，[ERA]表示地球旋转角矩阵，[W]表示极移矩阵，它们都可以通过调用SOFA的子程序进行计算。

相对转角理论：

对于给定的初始时刻t₀和下一时刻t₁，它们对应的变换总体矩阵分别为T₀、T₁，那么从t₀到t₁时刻的变换矩阵Q可以通过下式得到：

T₁＝QT₀

Q＝T₁T′₀＝([W(t₁)][ERA(t₁)][CIO(t₁)])([W(t₀)][ERA(t₀)][CIO(t₀)])′ (3)

＝[W(t₁)][ERA(t₁)][CIO(t₁)][CIO(t₀)]′[ERA(t₀)]′[W(t₀)]′

上式及以后的公式中，”矩阵’”(具有上撇符号的矩阵)均代表原矩阵的逆矩阵。

考虑到在短时间内(1天以内)，岁差章动矩阵[CIO]的变化很小：

[CIO(t₀)]′[CIO(t₁)]≈I

上式中，I为单位矩阵。

因此，

t时刻的总体变换矩阵T可以通过以下公式得到：

上式中，T₀可以通过固连系-惯性系转换矩阵计算得到。

通过这种方法就可以得到t₀至t₁时刻之间，任意时刻t所对应的坐标变换矩阵。

步骤二：初始参数计算

通过调用SOFA的子程序iauUtcut1计算给定IERS数据发布时刻t₀(UTC)、后一天的对应时刻t₁(即t₀+1)(UTC)对应的UT1时标：

iauUtcut1(t.UTCJD1,t.UTCJD2,&t.UT1JD1,&t.UT1JD2)) (6)

上式中，t.UTCJD1和t.UTCJD2表示在t时刻以UTC时标表示的两部分儒略日，同样地t.UT1JD1、t.UT1JD2为t时刻以UT1时标表示的两部分儒略日。需要分别计算发布时刻t₀和后一时刻t₁对应的UT1时标。

得到两个时刻对应的UT1时标后，再通过调用SOFA的子程序iauEra00计算初始时刻t₀、后一时刻t₁对应的地球旋转角ERA₀、ERA₁：

iauEra00(t.UT1JD1,t.UT1JD2) (7)

上式中，t.UT1JD1、t.UT1JD2为t时刻以UT1时标表示的两部分儒略日。需要分别计算发布时刻t₀和后一时刻t₁对应的地球旋转角ERA，并记

因此，可以计算得到t₀和t₁时刻之间的地球旋转平均角速度

调用SOFA子程序计算极移矩阵W(t₀)和岁差章动矩阵CIO(t₀)，在计算得到这两个3×3矩阵的相应数值后，通过以下公式将它们分解为欧拉角：

θ₀＝asin(-CIO[3])

ψ₀＝asin((CIO[2]/cosθ₀))

P_y＝asin(-W[3]) (9)

P_z＝asin((W[2]/cosP_y))

P_x＝asin((W[6]/cosP_y))

步骤三：输出星历插值表

执行一次步骤二，可以得到t₀、

θ₀、ψ₀、

P_x0、P_y0和P_z0的数值。使用C语言的fprint函数或C++的标准输出流函数将它们按照顺序存储至逗号分隔符文件(扩展名为.csv)中，这样就完成了一条星历的计算。

若将时间间隔(t₁-t₀)设置为1天，使用

…进行迭代循环计算N次，就可以得到一个条目数为N的星历插值表。通过这个插值表，可以计算t₀时刻至t₀+N天的时间内，任意时刻对应的固连系-惯性系的坐标变换矩阵。

在生成星历插值表后，计算上述时间段内的坐标变换矩阵时，无需调用任何SOFA的子程序，而且星历插值表只需要生成一次保存即可。

步骤四：读取星历插值表

使用C语言的fopen函数或C++的标准输入流函数，按行循环读取星历插值表中的数值，并将这些数值保存至多个一维数组或一个二维数组中，即可完成星历插值表的读取。

步骤五：初始化t时刻的相关矩阵

首先，使用二分法(或者其他高效的下标搜索方法)查找t所在的时间区间(t∈(t₀,t₁),t₁-t₀＝1)对应的t₀,t₁在数组中的下标。

获取t₀,t₁对应的下标后，再读取t₀对应的

θ₀、ψ₀、

P_x0、P_y0和P_z0的数值以及t₁对应的P_x1、P_y1和P_z1的数值。通过对P_x0、P_y0、P_z0以及P_x1、P_y1和P_z1线性插值(例如：

)，得到Px、P_y和Pz，从而完成参数的初始化。

通过式(2)和下面的公式，可以完成矩阵T₀、[W(t₀)]、

以及[W(t₁)]的初始化：

[W(t₀)]＝Rx(P_x0)Ry(P_y0)Rz(P_z0)

[W(t)]＝Rx(P_x)Ry(P_y)Rz(P_z)

上式中，Rx，Ry和Rz分别表示绕x、y、z轴旋转相应欧拉角生成的矩阵。

例如：

通过此步骤，即可完成相关矩阵的初始化工作。

步骤六：求解t时刻对应的坐标变换矩阵

步骤五已经初始化了[W(t₀)]矩阵，但是在式(5)中，需要它的逆矩阵。无论是岁差章动矩阵CIO还是极移矩阵W，它们的逆矩阵都与它们的转置矩阵相同。这个结论可以通过旋转矩阵的性质简单证明，此处不再赘述。

对[W(t₀)]进行转置，得到[W(t₀)]′。然后将它和其他相关矩阵的数值代入式(5)，通过三次矩阵乘法可以获得t时刻的总体变换矩阵T。

通过以上步骤，提出了一种基于相对转动的地球固连系-惯性系快速高精度转换方法；该方法仅在事先生成星历插值表时，调用SOFA的相关子程序。在生成星历表后，只需要查找表中包含对应待求时刻的区间，读取相关的参数并进行线性插值，最后代入相应的矩阵，完成三次乘法运算即可。

该方法只消耗少量的计算、不需要过多的天文理论背景、成本低廉易于实现，在航天的制导导航与控制等方面都具有重要的意义，可用于固连系坐标到惯性系坐标的快速变换。与传统的基于SOFA程序方法相比，这种新型的坐标转换方法大大降低了算法上的复杂度，在保证一定精度的同时能够以很小的计算代价完成相同的功能。

(三)优点

本发明提供的一种基于相对转动的地球固连系-惯性系快速高精度转换方法，相较于传统的坐标变换的优点在于：

1)本发明创新性地应用了相对转动理论，避免了大量SOFA子程序被反复调用造成的执行时间过长。仅仅需要计算三次线性插值和三次矩阵乘法即可完成大部分的运算，从而很好的提高了程序计算的效率。

2)本发明提出的利用星历插值表计算对应时刻的固连系－惯性系坐标变换矩阵的方法，理论简明易懂，使用方便；不同于在使用SOFA程序时，需要了解大量的参数以及调用的子程序对应的天文背景和复杂的定义。由于程序复杂度低、计算代价较小、执行效率高，本发明具有广泛的应用前景。

附图说明

图1是本发明所述方法流程图

具体实施方式

下面将结合图1和技术方案对本发明的具体实施过程做进一步的详细说明。

本发明提供的一种基于相对转动的地球固连系-惯性系快速高精度转换方法，见图1所示，其步骤如下：

步骤一：初始参数计算

通过调用SOFA的子程序iauUtcut1计算给定IERS数据发布时刻t₀(UTC)、后一天的对应时刻t₁(即t₀+1)(UTC)对应的UT1时标：

iauUtcut1(t.UTCJD1,t.UTCJD2,&t.UT1JD1,&t.UT1JD2)) (11)

上式中，t.UTCJD1和t.UTCJD2表示在t时刻以UTC时标表示的两部分儒略日，同样地t.UT1JD1、t.UT1JD2为t时刻以UT1时标表示的两部分儒略日。需要分别计算发布时刻t₀和后一时刻t₁对应的UT1时标。

得到两个时刻对应的UT1时标后，再通过调用SOFA的子程序iauEra00计算发布时刻t₀、后一时刻t₁对应的地球旋转角ERA₀、ERA₁：

iauEra00(t.UT1JD1,t.UT1JD2) (12)

上式中，t.UT1JD1、t.UT1JD2为t时刻以UT1时标表示的两部分儒略日。需要分别计算初始时刻t₀和后一时刻t₁对应的地球旋转角ERA，并记

因此，可以计算得到t₀和t₁时刻之间的地球旋转平均角速度

调用SOFA子程序计算极移矩阵W(t₀)和岁差章动矩阵CIO(t₀)，在计算得到这两个3×3矩阵的相应数值后，通过以下公式将它们分解为欧拉角：

θ₀＝asin(-CIO[3])

ψ₀＝asin((CIO[2]/cosθ₀))

P_y＝asin(-W[3]) (14)

P_z＝asin((W[2]/cosP_y))

P_x＝asin((W[6]/cosP_y))

该步骤如图1中第一个方框所示。

步骤二：输出星历插值表

执行一次步骤二，可以得到t₀、

θ₀、ψ₀、

P_x0、P_y0和Pz0的数值。使用C语言的fprint函数或C++的标准输出流函数将它们按照顺序存储至逗号分隔符文件(扩展名为.csv)中，这样就完成了一条星历的计算。

若将时间间隔(t₁-t₀)设置为1天，使用

…进行迭代循环计算N次，就可以得到一个条目数为N的星历插值表。通过这个插值表，可以计算t₀时刻至t₀+N天的时间内，任意时刻对应的固连系-惯性系的坐标变换矩阵。

在生成星历插值表后，计算上述时间段内的坐标变换矩阵时，无需调用任何SOFA的子程序，而且星历插值表只需要生成一次保存即可。

该步骤如图1中第二个方框所示。

步骤三：读取星历插值表

使用C语言的fopen函数或C++的标准输入流函数，按行循环读取星历插值表中的数值，并将这些数值保存至多个一维数组或一个二维数组中，即可完成星历插值表的读取。

该步骤如图1中第三个方框所示。

步骤四：获取t时刻所在区间

使用二分法(或者其他高效的下标搜索方法)查找t所在的时间区间(t∈(t₀,t₁),t₁-t₀＝1)对应的t₀,t₁在数组中的下标。

该步骤如图1中第四个方框所示。

步骤五：初始化t时刻的相关矩阵

获取t₀,t₁对应的下标后，再读取t₀对应的

θ₀、ψ₀、

P_x0、P_y0和P_z0的数值以及t₁对应的P_x1、P_y1和P_z1的数值。通过对P_x0、P_y0、P_z0以及P_x1、P_y1和P_z1线性插值(例如：

)，得到P_x、P_y和P_z，从而完成参数的初始化。

通过式(2)和下面的公式，可以完成矩阵T₀、[W(t₀)]、

以及[W(t₁)]的初始化：

[W(t₀)]＝Rx(P_x0)Ry(P_y0)Rz(P_z0)

[W(t)]＝Rx(P_x)Ry(P_y)Rz(P_z)

上式中，Rx，Ry和Rz分别表示绕x、y、z轴旋转相应欧拉角生成的矩阵。

例如：

通过此步骤，即可完成相关矩阵的初始化工作。

该步骤如图1中第五个方框所示。

步骤六：求解t时刻对应的坐标变换矩阵

步骤五已经初始化了[W(t₀)]矩阵，但是在式(5)中，需要它的逆矩阵。无论是岁差章动矩阵CIO还是极移矩阵W，它们的逆矩阵都与它们的转置矩阵相同。这个结论可以通过旋转矩阵的性质简单证明，此处不再赘述。

对[W(t₀)]进行转置，得到[W(t₀)]′。然后将它和其他相关矩阵的数值代入式(5)，通过三次矩阵乘法可以获得t时刻的总体变换矩阵T。

该步骤如图1中第六个方框所示。

通过以上步骤，提出了一种基于相对转动的地球固连系-惯性系快速高精度转换方法；该方法仅在事先生成星历插值表时，调用SOFA的相关子程序。在生成星历表后，只需要查找表中包含对应待求时刻的区间，读取相关的参数并进行线性插值，最后代入相应的矩阵，完成三次乘法运算即可。

该方法只消耗少量的计算、不需要过多的天文理论背景、成本低廉易于实现，在航天的制导导航与控制等方面都具有重要的意义，可用于固连系坐标到惯性系坐标的快速变换。与传统的基于SOFA程序方法相比，这种新型的坐标转换方法大大降低了算法上的复杂度，在保证一定精度的同时能够以很小的计算代价完成相同的功能。

Claims (1)

1.一种基于相对转动的地球固连系-惯性系快速高精度转换方法，其特征在于：其步骤如下：

步骤一：准备工作

地球固连参考系定义：

常用的地球固连参考系有国际地球参考系、协议地球参考系、1984世界大地坐标系和2000国家大地坐标系；国际地球参考系(International Terrestrial Reference System，简称ITRS)是目前最常用的地球参考系，它由IERS发布；IERS一直在对框架进行不断地改进和修正，至今IERS已经公布了13个版本，目前最新的版本是ITRF2014；该参考系的z轴被定义为与地球自转轴平行或重合，x轴则是指向本初子午线与地球赤道的交点；本发明中实现参考系变换时使用的地球参考系是ITRS；

地心惯性参考系定义：

常用的地心惯性参考系有J2000.0参考系、质心天球参考系和地心天球参考系(Geocentric Celestial Reference System，简称GCRS)；GCRS的坐标轴指向由甚长干涉(VLBI)测量的一组河外射电源在J2000.0天球赤道坐标系确定，是一个准空间固定参考系，可以便于表示天体在空间中的方位；J2000.0参考系是用J2000.0时刻的天球赤道与赤经章动来确定的天球参考系，即平春分点平赤道地心参考系，该参考系的坐标原点位于地球质心，其xoy平面是J2000.0时刻的平赤道面，x轴指向J2000.0时刻的平春分点；本发明中实现参考系变换时使用的地心惯性参考系是GCRS；

基础天文标准库(SOFA)：是由国际天文联合会(IAU)赞助的项目，旨在为天文计算过程中提供具有权威性的有效算法以及部分常数的数值；SOFA程序的第一版代码于2001年10月底发布，在此之后基本是每18个月都会进行版本维护；SOFA的算法最初是用Fotran77编写的，由于C语言的广泛应用及其高效性，在2009年2月IAU又推出了ANSI C版本的SOFA；目前，最新版SOFA(截止2018年10月)主要包含了天体测量学、日历、时间尺度、黄道坐标、地球自转和恒星时间、银河坐标、地心/大地转换、岁差，章动和极移以及星历转换等内容；

两部分儒略日：

本发明中使用的儒略日由两个变量(JD1与JD2)组成；JD1表示的是儒略日的整数部分，JD2表示的是儒略日的小数部分；

地球旋转角(Earth Rotation Angle，简称ERA)：

地球旋转角ERA₀是t₀时刻天球中间极赤道上CIO与TIO之间的角距离，它是UT1的线性函数；当t₁>t₀时，t₁时刻的地球旋转角ERA₁可以表示为：

上式中，

为t₀至t₁时刻的地球旋转平均角速度；

地球旋转角矩阵[ERA]是一个绕y轴旋转ERA角度的3×3矩阵；

固连系-惯性系转换矩阵：

[GCRS]＝[T][ITRS]＝[W][ERA][CIO][ITRS] (2)

上式中，[T]表示整个坐标变换过程的总体矩阵，[CIO]表示CIO理论中的岁差章动矩阵，[ERA]表示地球旋转角矩阵，[W]表示极移矩阵，它们都可以通过调用SOFA的子程序进行计算；

相对转角理论：

对于给定的初始时刻t₀和下一时刻t₁，它们对应的变换总体矩阵分别为T₀、T₁，那么从t₀到t₁时刻的变换矩阵Q可以通过下式得到：

上式及以后的公式中，“矩阵’”(具有上撇符号的矩阵)均代表原矩阵的逆矩阵；

考虑到在短时间内(1天以内)，岁差章动矩阵[CIO]的变化很小：

[CIO(t₀)]′[CIO(t₁)]≈I

上式中，I为单位矩阵；

因此，

t时刻的总体变换矩阵T可以通过以下公式得到：

上式中，T₀可以通过固连系-惯性系转换矩阵计算得到；

通过这种方法就可以得到t₀至t₁时刻之间，任意时刻t所对应的坐标变换矩阵；

步骤二：初始参数计算

通过调用SOFA的子程序iauUtcut1计算给定IERS数据发布时刻t₀(UTC)、后一天的对应时刻t₁(即t₀+1)(UTC)对应的UT1时标：

iauUtcut1(t.UTCJD1,t.UTCJD2,&t.UT1JD1,&t.UT1JD2)) (6)

上式中，t.UTCJD1和t.UTCJD2表示在t时刻以UTC时标表示的两部分儒略日，同样地t.UT1JD1、t.UT1JD2为t时刻以UT1时标表示的两部分儒略日；需要分别计算发布时刻t₀和后一时刻t₁对应的UT1时标；

得到两个时刻对应的UT1时标后，再通过调用SOFA的子程序iauEra00计算发布时刻t₀、后一时刻t₁对应的地球旋转角ERA₀、ERA₁：

iauEra00(t.UT1JD1,t.UT1JD2) (7)

上式中，t.UT1JD1、t.UT1JD2为t时刻以UT1时标表示的两部分儒略日；需要分别计算发布时刻t₀和后一时刻t₁对应的地球旋转角ERA，并记

因此，可以计算得到t₀和t₁时刻之间的地球旋转平均角速度

调用SOFA子程序计算极移矩阵W(t₀)和岁差章动矩阵CIO(t₀)，在计算得到这两个3×3矩阵的相应数值后，通过以下公式将它们分解为欧拉角：

步骤三：输出星历插值表

执行一次步骤二，可以得到t₀、

θ₀、ψ₀、

P_x0、P_y0和P_z0的数值；使用C语言的fprint函数或C++的标准输出流函数将它们按照顺序存储至逗号分隔符文件(扩展名为.csv)中，这样就完成了一条星历的计算；

若将时间间隔(t₁-t₀)设置为1天，使用

进行迭代循环计算N次，就可以得到一个条目数为N的星历插值表；通过这个插值表，可以计算t₀时刻至t₀+N天的时间内，任意时刻对应的固连系-惯性系的坐标变换矩阵；

在生成星历插值表后，计算上述时间段内的坐标变换矩阵时，无需调用任何SOFA的子程序，而且星历插值表只需要生成一次保存即可；

步骤四：读取星历插值表

使用C语言的fopen函数或C++的标准输入流函数，按行循环读取星历插值表中的数值，并将这些数值保存至多个一维数组或一个二维数组中，即可完成星历插值表的读取；

步骤五：初始化t时刻的相关矩阵

首先，使用二分法(或者其他高效的下标搜索方法)查找t所在的时间区间(t∈(t₀,t₁),t₁-t₀＝1)对应的t₀,t₁在数组中的下标；

获取t₀,t₁对应的下标后，再读取t₀对应的

θ₀、ψ₀、

P_x0、P_y0和P_z0的数值以及t₁对应的P_x1、P_y1和P_z1的数值；通过对P_x0、P_y0、P_z0以及P_x1、P_y1和P_z1线性插值(例如：

得到P_x、P_y和P_z，从而完成参数的初始化；

通过式(2)和下面的公式，可以完成矩阵T₀、[W(t₀)]、

以及[W(t₁)]的初始化：

上式中，Rx，Ry和Rz分别表示绕x、y、z轴旋转相应欧拉角生成的矩阵；

例如：

通过此步骤，即可完成相关矩阵的初始化工作；

步骤六：求解t时刻对应的坐标变换矩阵

步骤五已经初始化了[W(t₀)]矩阵，但是在式(5)中，需要它的逆矩阵；无论是岁差章动矩阵CIO还是极移矩阵W，它们的逆矩阵都与它们的转置矩阵相同；这个结论可以通过旋转矩阵的性质简单证明，此处不再赘述；

对[W(t₀)]进行转置，得到[W(t₀)]′；然后将它和其他相关矩阵的数值代入式(5)，通过三次矩阵乘法可以获得t时刻的总体变换矩阵T。

Images