ITRS/GCRS/J2000坐标系的相互转换
本文主要阐述了目前国际天文界最新规定的岁差章动模型(IAU 2000A/B)，并根据此模型给出ITRS,GCRS和
J2000平赤道地心系相互转换的详细步骤。
本文主要依据《IERS Conventions 2003》(IERS Technical Note No. 32)而来，由于文章为英文，且详细叙述了ITRS,GCRS参考系以及岁差章动等模型，会使得读者初次阅读（或天文背景知识不够）产生阅读障碍。笔者也是在多次反复阅读的基础上并参考其它书籍和文献才稍微弄清楚，因此本文是笔者关于此文献的一个总结，希望能够对读者有所帮助。
本文不打算对坐标系转换的理论（岁差章动，极移等）进行详细讲解，仅仅给出坐标转换的基本步骤和必要常识，以便读者根据此文能够迅速掌握ITRS,GCRS和J2000平赤道地心系相互转换的基本原理，并能够依据IERS提供的Fortran源程序进行实际编程应用。
前言
名词缩写和解释
BCRS barycentric celestial reference system
CEO celestial intermediate origin
CIP celestial intermediate pole
CIRS celestial intermediate reference system
EOP earth orientation parameters
GCRS geocentric celestial reference system
GMST Greenwich mean sidereal time
GAST Greenwich apparent sidereal time
IAU international astronomical union
ICRS international celestial reference system
IERS international Earth rotation and reference systems service
ITRS international terrestrial reference system
J2000 2000 January 1.5
SOFA standards of fundamental astronomy
TEO terrestrial intermediate origin
TIRS terrestrial intermediate reference system
TT terrestrial time
UT universal time UT1
UTC coordinated universal time
MAS milliarcsecond ( )
ITRS
就是我们常说的地固坐标系，其原点在地球质心（包含大气海洋等质量），坐标系xy平面为地球赤道面，z轴指向北极CIO处，x轴指向格林威治子午线与赤道面交点处。
此坐标系固定在地球上，地面站测控，以及地球引力场系数等都在此坐标系下定义。
J2000.0
常被称为J2000平赤道地心坐标系。其原点也是在地球质心，xy平面为J2000时刻的地球平赤道面，x轴指向J2000时刻的平春分点（J2000时刻平赤道面与平黄道面的一个交点）。
此坐标系常被作为地球卫星的惯性坐标系，卫星运动积分等都在此坐标系计算。
GCRS
J2000地心天球坐标系，其定义与J2000平赤道地心坐标系仅有一个常值偏差矩阵B。目前IAU推荐用此坐标系逐渐取代J2000平赤道地心坐标系。
时间量
UT1=UTC+(UT1-UTC) \\ DAT=TAI-UTC\\ TAI=UTC+DAT\\ TT=TAI+32.184s
给定UTC时刻，可求得其距离J2000的世纪数，即
t=(TT-2000 Jan 1d 12h TT) in days/36525
上式中 2000Jan1d12hTT2000 Jan 1d 12h TT对应的儒略日(Julian day)为2451545.0天。
坐标转换
当一坐标系绕其3个轴旋转 角时，则坐标旋转矩阵可表述为：
R_x(\theta)=\begin{pmatrix}1 & 0 & 0 \\0 & cos\theta & sin\theta \\0 & -sin\theta & cos\theta \\\end{pmatrix} \\R_y(\theta)=\begin{pmatrix}cos\theta & 0 & -sin\theta \\0 & 1 & 0 \\sin\theta & 0 & cos\theta \\\end{pmatrix} \\R_z(\theta)=\begin{pmatrix}cos\theta & sin\theta & 0 \\ -sin\theta & cos\theta & 0 \\0 & 0 & 1 \\\end{pmatrix}
其旋转方向符合右手螺旋法则，即逆时针旋转为正方向。另外坐标旋转矩阵具备如下性质：
R^{-1}(\theta)=R^T(\theta)=R(-\theta)
背景知识
我们知道，地球的自转轴在惯性空间中不是固定的，而是不断摆动的。此摆动造成地轴绕北黄级顺时针运动，夹角约为23.5度。于此同时，地轴还在做微小的抖动，见下图。前者的运动称为岁差(Precession)，后者运动成为章动(Nutation)。岁差章动的原因主要有两个方面。其一是太阳系行星对地球绕日轨道所产生的摄动影响；其二是太阳和月球对地球赤道隆起部分的摄动影响。
关于岁差章动的计算，此前一直采用IAU 1976岁差模型和1980章动模型。随着时间的推移，此模型的精度逐渐跟不上需要。因此，IAU又规定，自2003年1月1日起，采用新的岁差章动模型，即IAU 2000A模型（精度到达0.2mas），或者IAU 2000B模型（精度达到1mas）。
岁差章动模型描述了地球自转轴的运动，另外，由于地球表面的海洋，大气运动以及地核内部液体的运动造成地球自转轴并不是相对地球不动的；相对地球北极CIO点来说有个小范围的运动，此种现象成为极移。
此外，地球的自转也不是均匀的，也很复杂。
综上所述，从地固坐标系(ITRS)到地心惯性系(GCRS或J2000平赤道地心系)的坐标转换矩阵由极移，自转和岁差章动组成。
ITRS到GCRS的转换矩阵
本文采用IAU 2000A/B 岁差章动模型，在某历元UTC时刻，ITRS到GCRS的转换矩阵可写成：
\vec r_{GCRS}=Q(t)\cdot R(t)\cdot W(t)\cdot \vec r_{ITRS}=H\cdot G^T(t)\cdot \vec r_{ITRS}
其中， r⃗ ITRS\vec r_{ITRS}和 r⃗ GCRS\vec r_{GCRS}分别对应同一位置向量在ITRS和GCRS坐标系中的坐标。
上式中，W(t)W(t)，R(t)R(t)和Q(t)Q(t)分别对应极移，自转和岁差章动转换矩阵。
在计算R(t)R(t)和Q(t)Q(t)的时候，会有两种计算方法，我们分别称为CEO-based转换方法和equinox-based转换方法。其中前者为IAU提出的新的计算方法。
下面分别给出上述三个转换矩阵的求解过程。
极移矩阵 W(t)W(t)
W(t)=R_z(-s^{'})\cdot R_y(x_p)\cdot R_x(y_p)
上式中， s′s^{'}为：
s^{'}=-0.047mas\cdot t
极移量 (xp,yp)(x_p,y_p)的求解为：
(x_p,y_p)=(x,y)_{IERS}+(\Delta x,\Delta y)_{tidal}+(\Delta x,\Delta y)_{nutation}
极移量主要是由IERS根据天文观测给出的 (x,y)IERS(x,y)_{IERS},每周都有新的观测数据，此外，由于地球潮汐和章动的影响，会对极移有微小的修正 (Δx,Δy)tidal(\Delta x,\Delta y)_{tidal} 和 (Δx,Δy)nutation(\Delta x,\Delta y)_{nutation} 。
上式中， (x,y)IERS(x,y)_{IERS}由IERS给出的观测数据计算求得，(Δx,Δy)tidal+(Δx,Δy)nutation(\Delta x,\Delta y)_{tidal} +(\Delta x,\Delta y)_{nutation} 可由公式计算得到，IERS提供此fortran源程序。
地球自转矩阵R(t)R(t)
R(t)=R_z(-\theta)
地球自转角 θ\theta的求解根据转换方法的不同有不同的求解方式（CEO-based或者equinox-based），具体求解IERS给出了fortran源程序。
岁差章动矩阵 Q(t)Q(t)
前面提到，计算此矩阵有两种方法：
1. CEO-based方法：
这种方法是IAU最新提出的并极力倡导的，其计算公式为：
Q(t)= \begin{pmatrix}1-aX^2 & -aXY & X \\ -aXY & 1-aY^2 & Y \\-X & -Y & 1-a(X^2+Y^2) \\\end{pmatrix} \cdot R_z(s)
其中：
a=1/2+1/8(X^2+Y^2)
上式中：
(X,Y)=(X,Y)_{IAU2000}+(dX+dY)_{IERS}
(X,Y)IAU2000(X,Y)_{IAU2000}和 ss可根据IAU2000A/B岁差章动模型求解出，IERS同样给出求解的fortran源程序，另外，由于IAU2000A/B岁差章动模型没有包含地轴的高频率运动，所以要加上IERS通过观测数据给出的高频率修正项(dX+dY)IERS(dX+dY)_{IERS} 。
2. Equinox-based方法：
Q(t)=B\cdot P(t)\cdot N(t)
其中，常值偏差矩阵 BB ，岁差矩阵P(t)P(t)和章动矩阵 N(t)N(t)如下：
\begin{equation} \left\{ \begin{aligned} B &=R_z(-\delta\alpha_0)\cdot R_y(-\xi_0)\cdot R_x(\eta_0)\\ P(t) &=R_x(-\epsilon_0)\cdot R_z(-\psi_A)\cdot R_x(\omega_A)\cdot R_z(-\chi_A)\\ N(t) &=R_x(-\bar \epsilon)\cdot R_z(-\Delta\psi)\cdot R_x(\bar\epsilon+\Delta\epsilon) \end{aligned} \right. \end{equation}
章动量 (Δψ,Δϵ)(\Delta\psi,\Delta\epsilon)为：
(\Delta\psi,\Delta\epsilon)=(\Delta\psi,\Delta\epsilon)_{IAU2000}+(\delta\Delta\psi,\delta\Delta\epsilon)_{IERS}
上式中， =(Δψ,Δϵ)IAU2000=(\Delta\psi,\Delta\epsilon)_{IAU2000}由IERS 2000A章动模型给出。前面提到过，IAU 2000A/B模型提供的岁差章动不包含高频率项，而是由IERS的观测数据提供（上式右端最后一项），但是在IERS给出的观测数据中仅仅给出 (dX,dY)IERS(dX,dY)_{IERS},我们可以通过IERS提供的fortran源程序将 (dX,dY)IERS(dX,dY)_{IERS}转换为 (δΔψ,δΔϵ)IERS(\delta\Delta\psi,\delta\Delta\epsilon)_{IERS}。
其余参数皆为岁差参数，可以通过公式求出，此处从略。
值得一提的是常值偏差矩阵 中的参数也是给定的，在CEO-based方法求解中，此偏差是包含在(X,Y)IAU2000(X,Y)_{IAU2000}中的。
利用IERS提供的Fortran源程序进行转换
上节中详细讲述了ITRS到GCRS转换矩阵的求解过程，在实际应用中，如果是自己编写源程序的话是件非常琐碎的事情，因为IERS 2000A/B 章动模型的参数多达1000多项。幸而这些基本子程序IERS都提供了，我们所做的就是如何正确的运用这些源程序，并将它们组合起来。
在ftp://maia.usno.navy.mil/conv2000/chapter5/上可下载相关的源程序，子程序列表如下：
| 子程序名 | 说明 |
|---|---|
| BPN2000 | CEO-based intermediate-to-celestial matrix |
| CBPN2000 | equinox-based true-to-celestial matrix |
| EE2000 | equation of the equinoxes (EE) |
| EECT2000 | EE complementary terms |
| ERA2000 | Earth Rotation Angle |
| GMST2000 | Greenwich Mean Sidereal Time |
| GST2000 | Greenwich (apparent) Sidereal Time |
| NU2000A | nutation, IAU 2000A |
| NU2000B | nutation, IAU 2000B |
| POM2000 | form polar-motion matrix |
| SP2000 | the quantity s’ |
| T2C2000 | form terrestrial to celestial matrix |
| XYS2000A | X, Y, s |
| interp.f | Interpolation of IERS polar motion and UT1 time series |
| uai2000.f | IAU 2000 celestial pole offsets conversion (dpsi,deps,dX,dY) |
上表最后可从ftp://hpiers.obspm.fr/iers/models上下载得到。
IERS观测数据的处理
前面一再提到IAU 2000A/B章动岁差模型不包含高频率项，因此在完整的坐标转换过程中，必须考虑到IERS提供观测数据的高频率修正项。
IERS每周发布Bulletins A，每月发布Bulletins B，它们都是描述EOP的参数，下面是综合的EOP参数文件(C 04)，其主要内容如下：
注意上述文件中，dPsi,dEpsilon是针对IAU 1980章动模型的修正项(δΔψ,δΔϵ)IAU1980(\delta\Delta\psi,\delta\Delta\epsilon)_{IAU1980}，不是IAU 2000A章动模型的修正项(δΔψ,δΔϵ)IAU2000(\delta\Delta\psi,\delta\Delta\epsilon)_{IAU2000}，因此不需要使用。
根据UTC时刻的儒略日，加上数据列表MJD,x,y,UT1−UTCMJD,x,y,UT1-UTC，可调用interp.f文件中的interp子程序插值计算出对应UTC时刻的(xp,yp)(x_p,y_p)和(UT1−UTC)IERS(UT1-UTC)_{IERS} 。在子程序interp中，先插值计算出IERS的观测数据(x,y)IERS(x,y)_{IERS}和(UT1−UTC)IERS(UT1-UTC)_{IERS}，然后内部调用子程序PMUT1_OCEANS和PM_GRAVI计算由潮汐和章动引起的高频率修正项(Δx,Δy)tidal,(Δx,Δy)nutation(\Delta x,\Delta y)_{tidal},(\Delta x,\Delta y)_{nutation}和(UT1−UTC)tidal(UT1-UTC)_{tidal}，然后分别相加，给出最后的(xp,yp)(x_p,y_p)和UT1−UTCUT1-UTC 。
另外，根据IERS的观测数据列表MJD,dX,dY,LOD,DATMJD,dX,dY,LOD,DAT，可插值计算出(dX,dY)IERS,LOD(dX,dY)_{IERS},LOD和DATDAT。若为equinox-based方法转换，则需要调用uai2000.f文件中的子程序dXdY_dpsideps将(dX,dY)IERS(dX,dY)_{IERS}转换为(δΔψ,δΔϵ)IERS(\delta\Delta\psi,\delta\Delta\epsilon)_{IERS} 。
有了UT1−UTCUT1-UTC和DATDAT,则可求得UT1,TTUT1,TT和tt 。这些时间量在以后的子程序中都需要。
具体转换步骤
首先调用子程序SP2000求得s′s^{'}，再由上面插值求得的(xp,yp)(x_p,y_p)，调用子程序POM2000即可求得极移矩阵W(t)W(t)。
求地球自转和岁差章动矩阵有两种方法，下面分别叙述：
1. CEO-based transformation
调用ERA2000求得地球自转角θ\theta；
然后调用子程序XYS2000A求得(X,Y)IAU2000(X,Y)_{IAU2000} 和ss，再加上上面观测数据插值的(dX,dY)IERS(dX,dY)_{IERS}，则可求得(X,Y)(X,Y)。根据(X,Y,s)(X,Y,s)，利用子程序BPN2000即可求得岁差章动转换矩阵Q(t)Q(t) 。
2. Equinox-based transformation
调用子程序NU2000A求得章动量(Δψ,Δϵ)IAU2000(\Delta\psi,\Delta\epsilon)_{IAU2000}，再加上由观测数据插值求得的(dX,dY)IERS(dX,dY)_{IERS} 转换后的(δΔψ,δΔϵ)IERS(\delta\Delta\psi,\delta\Delta\epsilon)_{IERS} ，得到最后的(Δψ,Δϵ)(\Delta\psi,\Delta\epsilon)；
由Δϵ\Delta\epsilon，调用GST2000即可求得地球自转角θ\theta;
然后再由(Δψ,Δϵ)(\Delta\psi,\Delta\epsilon)调用子程序CBPN2000，求得岁差章动矩阵Q(t)Q(t)。此处需要对子程序CBPN2000进行简单的说明，其内部进行常值偏差矩阵B和岁差章动矩阵P,N的计算，最后给出矩阵Q(t)Q(t)。
由上述两种方法之一求得W(t),θW(t),\theta，和Q(t)Q(t)，调用子程序T2C2000即可求得ITRS到GCRS的转换矩阵H⋅GT(t)H\cdot G^T(t)。
若采用方法二时，可以用子程序NU2000B替代NU2000A，其它都不变，此种转换的精度稍低(1mas)，但是其计算速度会快很多，在精度要求不是很高的情况下采用此种方法可使计算速度大大提高。
其它的一些说明
1.) 下面给出GCRS和ITRS（包含速度）转换的完整表述：
\begin{equation} \left\{ \begin{aligned} \vec r_{ITRS} &=[Q(t)\cdot R(t)\cdot W(t)]^T\cdot \vec r_{GCRS}=H\cdot G(t)\cdot \vec r_{GCRS}\\ \vec r_{TIRS} &=[Q(t)\cdot R(t)]^T\cdot \vec r_{GCRS}\\ \vec V_{ITRS} &=W^T(t)\cdot \{R^T(t)\cdot Q^T(t)\vec V_{GCRS}-\vec \omega_e\times \vec r_{TIRS}\} \end{aligned} \right. \end{equation}
和
\begin{equation} \left\{ \begin{aligned}\vec r_{GCRS}&=Q(t)\cdot R(t)\cdot W(t)\cdot\vec r_{ITRS} =H\cdot G^T(t)\cdot\vec r_{ITRS}\\ \vec r_{TIRS} &=W(t)\cdot \vec r_{ITRS}\\ \vec V_{GCRS} &=Q(t)\cdot R(t)\cdot \{W(t)\cdot \vec V_{ITRS}+\vec \omega_e\times \vec r_{TIRS}\} \end{aligned} \right. \end{equation}
其中
\omega_e=7.292115146706979\times10^{-5}\{1-LOD/86400\}
上式中的LOD由IERS的观测数据插值求得。
注意上述公式中 r⃗ TIRS\vec r_{TIRS}和 r⃗ ITRS\vec r_{ITRS}的区别， r⃗ TIRS\vec r_{TIRS}为地固系ITRS坐标 r⃗ ITRS\vec r_{ITRS}经过极移转换矩阵后的坐标。
2.) 若需要GCRS和J2000平赤道地心系的相互转换，只需要常值偏差矩阵B即可。其具体的源程序需要读者自行编制。
\begin{equation} \left\{ \begin{aligned}\vec r_{GCRS}&=B\cdot \vec r_{J2000}\\ \vec V_{GCRS} &=B\cdot \vec V_{J2000}\\ \vec r_{J2000} &=B^T\cdot \vec r_{GCRS}\\ \vec V_{J2000} &=B^T \vec V_{GCRS} \end{aligned} \right. \end{equation}
3.) 在IERS提供的fortran源程序中，有部分子程序需要调用IAU SOFA软件包中的子程序，有关IAU SOFA软件包的说明和使用请参见我的文档《IAU SOFA软件包介绍》。
4.) IAU SOFA软件包中也包含有关岁差章动和极移等基本子程序，其主要内容和本文介绍的fortran源程序大同小异。但其软件包中还包含以前的岁差章动模型以及最新的IAU2006岁差模型，读者需要的话可参考其文档说明“sofa_pn.pdf”，并有具体例子，强烈建议读者阅读。
参考文献
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，以及一些其它文件的说明。
ITRS/GCRS/J2000坐标系的相互转换相关推荐
- GNSS-ITRS_GCRS_J2000坐标系的相互转换
背景 转载自:https://blog.csdn.net/u011575168/article/details/52081409 本文主要阐述了目前国际天文界最新规定的岁差章动模型(IAU 2000A ...
- 机械手基坐标系和工具坐标系的相互转换
机械手基坐标系和工具坐标系的相互转换,主要是通过欧拉角来完成的. double D2R(double D) {return D / 180 * M_PI; }double R2D(double R) ...
- ArcGIS坐标系篇之几种常见坐标系的相互转换（二）
几种常见坐标系的相互转换 坐标系的转换是我们在地理数据的使用过程中经常会遇到的问题,像大家经常提到的54转8080转2000WGS84转2000等等,以及ArcGIS经常提到的地理坐标系转投影坐标系. ...
- wgs-84,gcj-02,bd-09的相互转换，高德，世界测量，百度坐标系的相互转换，坐标系转换
高德使用的是gcj-02坐标系,百度使用的是bd09坐标系,注意其间的相互转换 1.首先创建一个GPS对象类 public class Gps {private double wgLat;privat ...
- 浅析相机相关坐标系的相互转换（世界坐标系、相机坐标系、图像坐标系、像素坐标系、内参矩阵、外参矩阵、扭转因子）【相机标定计算机视觉】
最近在网上看到了很多关于坐标系转换的帖子,但是其内容很多都是相互转载(甚至还有一部分是错误的),同时大部分的文章内容都是告诉你四种坐标系间的相互转化的数学公式,看完之后很多时候还是不知所云,本文意在深 ...
- 2D坐标系与3D坐标系的相互转换--python实现
并不是做关于SLAM方向的,但由于某些任务涉及到,故作此笔记- 相机内参矩阵: 不同的的深度摄像头具有不同的特征参数,在计算机视觉里,将这组参数设置为相机的内参矩阵C: [fx0cx0fycy001] ...
- 百度地图（BD-09）坐标系和高德（GCJ-02）坐标系的相互转换（PHP版）
简单介绍一下常见的地理经纬度坐标系: WGS-84坐标系: 是国际标准坐标系,GPS坐标(Google Earth使用.或者GPS模块). GCJ-02坐标系: 火星坐标系,GCJ-02是由中国国家测 ...
- php百度坐标转腾讯坐标,火星坐标系、百度坐标系、腾讯坐标系的相互转换
/** * 火星坐标系 (GCJ-02) 与百度坐标系 (BD-09) 的转换算法 将 GCJ-02 坐标转换成 BD-09 坐标 * * @param gg_lat * @param gg_lon ...
- STK中心天体cb等相关文件说明
STK中的中心天体就是指太阳及太阳系中的行星以及它们的卫星(如月球等).我们在进行相关计算时,有关中心天体的椭球形状.姿态和星历等参数是缺省设置的,在必要的时候我们可以修改相应的参数设置,因此必须清楚 ...
- 电子围栏绘制、校验、不同坐标系（wgs84，jcg02、bd09）的相互转换，看这一篇就够了！
作为一名前端开发者,在工作中经常会用到地图相关的API,比如百度.高德.腾讯等,但是不同地图平台使用的坐标系可能不一样,这就导致定位时会出现几百米到1000左右的误差,比如在小程序平台获取位置的时候我 ...
最新文章
- html 字幕飘动效果,html 滚动字幕 制作滚动字幕效果 参数
- 美轮美奂渐变炫彩海报素材，荣盛回来！
- BERT源码分析PART I
- HandlerThread使用
- win10亮度怎么调_装好原版win10后，没有这样设置过，你的系统怎么可能好用！...
- 简单介绍三级分销系统开发源码
- Verilog学习之异步复位的串联T触发器设计
- 实用好看的个人主页官方源码
- 用户如何制作360度全景图？360度全景图有什么用？
- OpenStack 企业私有云的若干需求（7）：电信行业解决方案 NFV
- theHarvester使用
- 回溯法之活动安排问题
- python中的arange函数,Python中的arange是什么?和range有什么不同？
- python中glob的用途
- 冷暖自知！4年Java小伙收获美团Offer，分享他的社招Java岗4面面经
- PyQT绘制股票K线图
- iOS苹果手机上最好用的3个的小说阅读器
- Python中For循环与While 循环的使用
- git 清除所有历史记录
- 芯片制造22nm制程是什么含义