GCRS与ITRS相互转换，时间系统及转换

大作家佚名

已于 2023-02-14 10:13:03 修改

阅读量6.3k

点赞数 10

分类专栏：算法/方法文章标签： ICRS ITRS

于 2022-10-08 16:52:50 首次发布

本文链接：https://blog.csdn.net/wokaowokaowokao12345/article/details/127211085

版权

算法/方法专栏收录该内容

37 篇文章 6 订阅

订阅专栏

1. 简介

关于坐标系、时间系统转换一直停留在课本上范范而谈的概述上，实际工程的应用较书本上更加详细，涉及到的内容也更多。这里要实现的目标是坐标系统的转换，将涉及到的内容做一些记录。其中，难免有描述不妥或理解不当之处，欢迎指出讨论。

GCRS 是用于地球上观测天体的坐标系，它主要用于地球天文学观测，例如控制望远镜的姿态和观测设备的对准等。它对于地球上的天文学应用非常重要，但不适用于对太阳系的整体研究。

国际地球参考系 ITRS 是我们常用的坐标系，一般性测量涉及的坐标都是这个坐标系。为了更方便地研究卫星问题, 经常需要在这两个坐标系下进行转换。IAU 的标准基本天文程序库 SOFA 给出基本天文运算的库函数, 利用这些库函数, 可以方便地实现这两个坐标系的转换。

本文主要介绍了 ITRS 与 GCRS的发展以及基于 SOFA 的 ITRS 与 GCRS之间不同的转换方法的实现。基础天文标准库 (SOFA) 是国际地球自转服务(IERS)协议提供的关于地球姿态、时间尺度和历法的一系列程序集。

提示：目前广为使用的协议天球坐标系是由 IAU 规定的国际天球坐标系ICRS。依据坐标原点的不同，ICRS可分为太阳系质心天球参考系BCRS和地心天球参考系GCRS两类。前者的坐标原点位于太阳系质心，用于计算行星的运动轨道，编制星表；后者的坐标原点位于地心，用于计算卫星轨道，编制卫星星历。博文中，有些内容想要表达的是GCRS与ITRS之间的转换，写成了ICRS与ITRS之间的转换（一些文献中模糊了ICRS与GCRS的概念，其实应该写清楚的），不够严谨，但是表达的意思是GCRS与ITRS之间的转换！

2. GCRS

2.1 GCRS

GCRS (Geocentric Celestial Reference System) 是一个以地球为原点的天文学坐标系，它与地球相关。在 GCRS 中，天体的位置是相对于地球的，而不是太阳系的中心。

GCRS 具有和 ICRS 相似的坐标系结构，其中一个坐标轴原点对准地球的中心，另一个坐标轴在地球上的正南方向，第三个坐标轴在地球上的正东方向。不同的是，GCRS 是一个动态的坐标系，它随着地球的运动而变化，并且受到地球摄动等因素的影响。

2.2 岁差、章动、极移

在理解坐标系之前需要对地球岁差（进动）、章动和极移有一定的了解，就是要对你所参考的物体有所了解。

岁差： 岁差是由于太阳、月球和行星对地球的吸引力造成的。
（1）日月岁差：太阳和月球对地球赤道隆起部分的吸引力产生不可抵消的力偶作用，从而使得地球自转轴在惯性空间中绕黄道轴进动。北天极绕着黄道轴做进动，进动周期约25700年，每年向西移动50″.37，同时也有章动。
（2）行星岁差：由于太阳系内其它行星引力作用，地球的周年运动并不严格遵守开普勒定律，黄道面的位置是在不断变化的，黄道面的变化也会导致春分点的变化，这使得春分点每年沿着天赤道东进约0″.13。

章动： 在不考虑章动的情况下天极绕黄道轴做圆周运动，此时天极称为平天极。若考虑章动影响则为真天极。
极移： 地球自转轴相对于地球本体是运动的，由此造成的地球极点在地球表面上的位置随时间发生变化，这种现象称为地极移动，简称极移。

极移与岁差章动是两种不同的地球物理现象。岁差和章动是地球自转轴在惯性空间的运动，其在地球内部的相对位置并没有改变，因此它只引起天体坐标的变化，不会引起地球表面经纬度的变化。与之相反，极移是地球自转轴在地球本体内的变化，其在惯性空间中的方向并没有改变，因而会引起地球表面上各地经纬度的变化。

综上所述：
（1）日月岁差会引起地球天极变化，春分点变化；
（2）行星岁差会引起春分点变化；
（3）章动会引起天极变化，春分点变化；
（4）极移会引地表各地经纬度变化；

内容出处：关于J2000.0坐标系与WGS84坐标系的理解。

2.3 什么是赤道坐标系

要理解 ICRS 首先要知道什么是赤道坐标系，因为它是基于赤道坐标系的概念建立的。参考：国际天球参考系统 —— ICRS 简介。

赤道坐标系的极点就是地球自转轴与天球的交点，其赤道就是地球的赤道。下图展示了如何从地球的自转构建赤道坐标系。
在这里插入图片描述

2.4 ICRS简介

国际天球参考系 (International Celestial Reference System，ICRS) 由国际天球参考框架 (ICRF) 实现。1997 年 IAU 第 23 届大会上, 通过并决定自 1998 年 1 月 1 日起, 在天文研究、空间探测、大地测量以及地球动力学等领域中采用 ICRS。ICRS 中的坐标系统与赤道坐标系基本一致，重点在于如何选取参照物 ：ICRS J2000.0 的平均极点在 12 h 的方向是 17.3±0.2 mas，在 18 h 的方向是 5.1±0.2 mas；平均分点在垂直于极轴的方向上从 ICRS 赤经原点移动了 78±10 mas。
参考：国际天球参考系统 —— ICRS 简介、维基百科、百度百科、www.iers.org

下图描述了天球坐标系。
在这里插入图片描述

2.5 ICRS 的实现

我们都知道，确定一个物体的运动（运动事实上也就是位置）需要一个参照物。依据一个特定的参照物，选取一个坐标系，我们会得到一个物理上的参考系。

在通常的赤道坐标系中，我们规定太阳的春分点所在的位置为赤经零点，但是事实上，春分点并不是恒久不变的。由于地球的进动以及章动，春分点无时无刻不在运动着。进动是指地球自转轴绕着空间中的另一个轴旋转，章动则是在进动的过程中更加微小的运动。P 和 N 分别是进动与章动。
在这里插入图片描述
因此，ICRS 选择使用几乎 “恒久不动” 的河外天体作为参照物，这就是 ICRF（International Celestial Reference Frame）—— 国际天球参考框架。

ICRF 是一系列河外源的精确位置，这些源通常是遥远的类星体，它们足够遥远，以保证这一参考框架的 “惯性” 性以及稳定性，这其中就有著名的也是第一个被发现类星体 3C 273。通过 VLBI—— 甚长基线射电干涉测量技术，我们可以将这些射电源的位置精确到数个毫角秒以内，以保证 ICRS 的精确性。

ICRS的原点为地球质心 ，地球质心指向天球北极的方向为Z轴方向，地球质心指向春分点的方向为X轴方向，利用右手法则确定Y轴方向，该坐标系描述卫星在轨飞行时的运动状态。使用该坐标系时一般采用国际组织于1984年启用的协议天球坐标系J2000。下图为空间固定惯性参考系示意图。

在这里插入图片描述

2.6 J2000

J2000 (Julian Year 2000) 是一种天文坐标系，是ICRS (International Celestial Reference System) 的一个参考点。它的原点在太阳系的重心，三个坐标轴分别指向瞬时天球的赤道平面的三个方向。

J2000 是以 2000 年 1 月 1 日午夜 12 时为参考点，以此作为太阳系对象的坐标系。它是天文学研究中常用的参考系，主要用于计算天体位置、运动以及对天体进行分类和分组。

J2000 坐标系非常精确，并且在天文学、天体物理学和天体动力学研究中非常重要。但是，随着时间的推移，天体的位置和运动会有所变化，因此 J2000 坐标系在长期内也需要不断校正和更新。J2000可理解为国际天球惯性参考坐标系ICRS的一种。

此坐标系常被作为地球卫星的惯性坐标系，卫星运动积分等都在此坐标系计算。

3. ITRS

3.1 ITRS简介

1988 年国际大地测量学与地球物理学联合会 (IUGG) 和国际天文学联合会 (IAU) 联合创建了国际地球自转服务局 (IERS), IERS 是国际地球参考框架 (ITRF) 的发布机构, 定期求取 ITRF 的年度解, 这些解由 IERS 年报和技术备忘录公布, 至今 ITRF 已经公布了 12 个版本, 分别为 ITRF88、 ITRF89、ITRF90、ITRF91、ITRF92、ITRF93、 ITRF94、ITRF96、ITRF97、ITRF2000、ITRF2005、 ITRF2008 。其中, ITRF2008 是目前最新的版本。

ITRS （International Terrestrial Reference System）就是我们常说的地固坐标系，其原点在地球质心（包含大气海洋等质量），坐标系xy平面为地球赤道面，地球质心指向地球北极的方向为Z轴方向，x轴指向格林威治子午线与赤道面交点处。此坐标系固定在地球上，地面站测控，以及地球引力场系数等都在此坐标系下定义。该坐标系描述卫星的位置和飞行速度，使用该坐标系时，常采用(World Geodetic System 84,WGS84)坐标系来代替，由于两者十分接近，在使用过程中忽略它们的差别。使用WGS84坐标系可以有效消除卫星运行中地球曲率和自转等因素影响。下图为地固坐标系示意图，其与ICRS明显不同之处是X轴指向，在宇宙空间中ITRS坐标系随地球而自转（相对地球不变），ICRS不会随着地球自转而变化，因为其X轴指向春分点（天球中固定方向）。
在这里插入图片描述

它是大地测量学和地球动力学研究的一种基本坐标系。如果把地球潮汐和地壳运动忽略不计，地球重力场和地面点的位置在这个坐标系中是固定不变的。也就是说这个坐标系仅随地球自转而转动，固定在地球上不变，因而也被称为地固坐标系。地球坐标系的建立已有一百多年的历史，1980年以前主要采用的是光学观测。随着空间大地测量的开展，观测人造的或自然的天体打破了集团或国家独有的观测传统，迫切要求确立与使用公用的地球坐标系。但宇宙间不可能存在绝对固定不动的东西，所以建立这种坐标系只能通过一种协议结果来体现，因而这种坐标系也被称协议地球参考系（CTRS），它是国际上约定统一采用的地球参考系。世界各国经常使用的国际地球参考框架就是这种坐标系的实现。

国际地球自转局（IERS）主要任务是准确及时提供自转参数。它的目的之一，就是建立与保持这个地球参考框架。该框架是国际大地测量学和地球物理学联合会（IUGG）、国际大地测量学协会（IAG）、国际天文学会（IAU）专门决定建立的，有关工作由IERS下属地球参考框架部门负责执行。具体由设在法国巴黎的国家地理院(IGN)大地测量室(LAREG)主持。使用的空间大地测量技术为：激光测月（LLR）；激光测卫（SLR）；甚长基线干涉（VLBI）；全球定位系统（GPS）；多普勒卫星跟踪和无线电定位系统（DORIS）等。

国际地面参考框架（ITRF，International Terrestrial Reference Frame）是ITRS的实现，每隔几年就会产生新的ITRF解决方案，使用最新的数学和测量技术，尽可能精确地实现ITRS。由于实验误差，任何给定的ITRF都将与已有ITRF略有不同。ITRS和ITRF解决方案由国际地球自转和参考系统处（IERS）维护。

4. 时间系统

时间系统满足的两个条件：（1）能作连续均匀的周期性运动, 且运动周期十分稳定;（2）运动周期具有很好的复现性, 即在不同的时期和地点这种周期性的运动都可以通过观测和实验来予以复现。

详细内容可以参考SOFA软件包中的介绍，SOFA提供了Fortran77和C语言的时间转换代码，Matlab论坛提供了对应的Matlab版本。

计时是按照约定的方法来测量时间，在人类历史上，长期以来地球自转是可用的最佳计时器。我们熟悉的公历日期由年、月和日组成，平均年长度为365.25天。天文学家分别使用365.25天和36525天作为 Julian year 和
Julian century。对于更大的间隔，其他天文现象也发挥了作用，特别是地球和月球的轨道周期。

任何科学中的计算都可能涉及精确的时间，但是目前存在多种时间系统。这有几个原因：天文学家必须继续处理过时时间尺度背后的现象，特别是地球自转和行星运动；随着新的时间尺度的引入，与过去的连续性得以继续；在天文应用中，“时钟”的物理背景很重要，无论它是在地球上、移动中还是静止中，还是在航天器上。

常用时间系统如下：

项目	Value
TAI (International Atomic Time)	官方计时标准。
UTC（Coordinated Universal Time）	民用时间的基础。
UT1（Universal Time）	基于地球自转。
TT（Terrestrial Time）	用于太阳系星历查找。
TCG（Geocentric Coordinate Time）	用于以地球为中心的空间计算。
TCB（Barycentric Coordinate Time）	用于地球轨道以外的计算。
TDB（Barycentric Dynamical Time）	TCB的一种标度形式，平均与TT保持一致。

在这里描述的七个时间尺度中，一个是原子时间（TAI），一个为太阳时间（UT1），一种是原子/太阳混合时间（UTC），四个是动力学时间（TT、TCG、TCB、TDB）。每一个都有不同的作用，其中一些之间有几十秒的偏移：当计划天文计算时，选择正确的计算是至关重要的。一个特别常见的错误是，假设只有一种精确的时间，即UTC，它与从望远镜指向（实际上需要UT1）到查找行星位置（需要TDB，可以用TT近似）的一切都兼容。事实上，UTC本身几乎从来都不是天文计算所用的时间刻度，也许除了用于记录之外。下图是一个“路线图”，显示了七个时间尺度之间的关系，以及如何将一个尺度上的时间转换为另一个尺度的同一时间。时间尺度由矩形框表示，其关系显示在圆形框中。
在这里插入图片描述

4.1 时间系统转换关系

∆UT1 = UT1−UTC
∆AT = TAI−UTC
∆T = TT − UT1 = 32.184s + ∆AT − ∆UT1
∆TT = TT−UTC

TAI 1997 Jan 1 00:00:00
TCG 1997 Jan 1 00:00:00 + 32.184 s
TCB 1997 Jan 1 00:00:00 + 32.184 s
TDB 1997 Jan 1 00:00:00 + 32.184 s − 65.5 µs

TT = TAI + 32.184 s.

TCG = TT + LG × (JDTT − TT0)

其中，TT0 = 2443144.5003725, JDTT 是 TT 表示为儒略日, LG = 6.969290134 × $10^{−10}$ 。

TDB = TCB − LB × (JDTCB − T0) × 86400 + TDB0

其中，JDTCB 是TCB 的儒略日, T0 = 2443144.5003725, LB = 1.550519768 × $10^{−8}$ and TDB0 = −6.55 × $10^{−5}$ s.

4.2 历书时(ET－ephemeris time )

由于地球自转的不均匀，以地球自转为基准的时间系统不能满足越来越高的精度要求。于是在1958年国际天文联合会（iau）决定，自1960年起各国开始使用以地球公转周期为基准的历书时编算天文年历。

历书时从太阳平黄经（mean ecliptic longitude）279°41′48．04″的瞬间起算，这个瞬间定义为历书时1900年1月1日12时正。历书时的基本时间单位历书秒（ephemeris second）等于1900年1月1日12时瞬间的回归年长度的31556925．9747分之一。

历书时的观测精度随着观测设备和技术的提高而逐渐增强，从早期的et。、et1发展到后来的et2，其精度可达到 $10^{-11}$ ~ $10^{-9}$ 量级。这个精度高于恒星时、世界时等以地球自转为基准的时间系统，但由于不能实时给出精确的历书时，且精度有限，所以它在使用上有一定的局限性。1976年国际天文联合会决定，自1984年起采用地球力学时和太阳系质心力学时取代历书时。

4.3 恒星时、真太阳时、平太阳时、区时

（1）恒星时：恒星时是以春分点为参考点。由于地球自转使春分点连续两次经过地方上子午圈的时间间隔为一恒星日。

（2）真太阳时：真太阳时是以太阳中心作为参考点的, 太阳中心连续两次通过某地的子午圈的时间间隔称为一个真太阳日。椭圆轨道、角速度不同导致真太阳时的长度不相同。

（3）平太阳时：用一个假太阳来代替真太阳构建的时间系统。它每年和真太阳同时从春分点出发，在天球赤道上从西向东匀速运行，这个速度相当于真太阳在黄道上运行的平均速度，最后和真太阳同时回到春分点。

（4）区时：全球分为 24个标准时区。在同一时区, 统一采用该时区中央子午线上的平太阳时, 称为区时。

4.4 世界时（UT－universal time ）

世界时就是以平子夜为零时起算的格林尼治平太阳时。世界时与平太阳时的单位尺度相同，只是起算点不同，前者比后者晚12个小时。世界时作为一种以地球的周期性自转运动为基准的时间系统，它的精确度主要受地球自转运动速度不均匀和极移的影响。

世界时问题：
地球自转的速度是不均匀的, 它不仅有长期减缓的总赲势, 而且也有季节性的变化以及短周期的变化, 情况较为复杂;
地极在地球上的位置不是固定不变的, 而是在不断移动, 即存在极移现象。这就意味着世界时已不再严格满足作为一个时间系统的基本条件, 因为它实际上不是一个完全均匀的时间系统。

4.5 原子时、国际原子时（TAI-international atomic time）

由于物质内部原子跃迁时辐射和吸收的电磁波频率具有较高的稳定性和复现性，由此建立的原子时，是目前最理想的时间系统。国际上原子时秒长的定义为：位于海平面上的铯原子基态两个超精细能级，在零磁场中跃迁辐射振荡9162631770周所持续的时间，为一原子秒（atomic second）。原子时的原点则由下式定义；

为了统一世界各国的原子时，国际上将100多台原子钟，通过相互比对，经由数据处理推算出统一的原子时系统，称为国际原子时（international atomic time）。TAI的稳定度大约为 $10^{-13}$ 左右。在gps定位中，原子时作为高精度的时间基准，用来测定卫星信号的传播时间。

4.6 协调世界时（UTC－universal coordinated time）

协调世界时是以原子时秒长为基础，在时刻上尽量接近于世界时的一种时间计量系统，是介于原子时和世界时之间一种折衷的时间系统。

国际原子时的准确度为每日数纳秒，而世界时的准确度为每日数毫秒。许多应用部门要求时间系统接近世界时UT，对于这种情况，一种称为协调世界时的折中时标于1972年面世。为确保协调世界时与世界时相差不会超过0.9秒，在有需要的情况下会在协调世界时内加上正或负闰秒。因此协调世界时与国际原子时之间会出现若干整数秒的差别，两者之差逐年积累，便采用跳秒（闰秒）的方法使协调时与世界时的时刻相接近，其差不超过1s。它既保持时间尺度的均匀性，又能近似地反映地球自转的变化。按国际无线电咨询委员会（CCIR）通过的关于UTC的修正案，从1972年1月1日起UTC与UT1（在UT中加入极移改正得到）之间的差值最大可以达到±0.9s。位于巴黎的国际地球自转事务中央局负责决定何时加入闰秒。一般会在每年的6月30日、12月31日的最后一秒进行调整。具体日期由国际地球自转与参考系统服务组织（IERS）安排并通告。

目前，世界各国播发的时间，均以UTC为基准。

4.7 GPS时间系统（GPST - gps time）

为精密导航和测量需要，全球定位系统建立了专用的时间系统，由GPS主控站的原子钟控制。 GPS时属于原子时系统，秒长与国际原子时相同，但与国际原子时的原点不同，即GPST与TAI在任一瞬间均有一常量偏差。

TAI-GPST = 19s

GPST与UTC的时刻，规定在1980年1月6日0时一致，随着时间的积累，两者的差异将表现为秒的整数倍。GPST与UTC之间关系

GPST=UTC+1S×n−19sGPST=UTC+1S×n−19s

到1987年，调整参数n为23，两系统之差为4秒，到1992年调整参数为26，两系统之差已达7秒。

4.8 儒略日（JD）/修正儒略日（MJD）（Julian Date /Modified Julian Date ）

出于许多目的，日历日期很不方便。为此，可以使用儒略日编号系统。JD zero位于7000年前，远早于历史时代，定义为格林威治正午；例如，JD 2449444开始于1994年4月1日正午（12点）。JD 2449443.5是1994年4月1日开始的午夜（0点）。由于儒略日的值太大，加上半天偏移的笨拙，删除前导“24”和“尾0.5”是公认的做法，由此产生所谓的修正儒略日：

MJD = JD − 2400000.5

Matlab 中测试：

% 2018年8月7日0时0分0秒
mjd = mjuliandate(2018,8,7,0,0,0)
jd = juliandate(2018,8,7,0,0,0)

结果

mjd = 58337                 jd = 2458337.5

MJD通常用于计算机应用程序，而不是JD本身，以减少舍入误差的风险。

出于某些天文目的，使用分数年很方便，例如，接近2009年底的给定日期和时间可以写为“2009.93”。

5. ITRS转GCRS理论实现

博文中流程皆为ICRS转换到ITRS的流程。在以下流程的基础上矩阵求逆即可实现ITRS到ICRS的转换。转换采用的式SOFA软件包，基础天文标准库 (SOFA), 是 IAU 赞助的项目, 旨在为天文计算提供权威有效的算法程序和常数数值。 1997 年, SOFA 评审委员会正式创立, 并设置了发布代码的 SOFA 中心。本文中坐标参考系之间转换采用的基础代码来源于该中心。在博文末提供了SOFA链接。

5.1 坐标系统转换

ICRF与ITRF转换方法目前主要有基于春分点的经典转换方法,和天球中间原点(CIO)的新方法。

基于春分点的ICRF与ITRF转换的经典方法
在基于春分点的转换方法中ICRF转化到ITRF要经历了ICRF转化为瞬时平天球坐标系、瞬时平天球坐标系转化为瞬时天球坐标系、瞬时天球坐标系转化为瞬时地球坐标系、瞬时地球坐标系转化为ITRF四个过程。ITRF转换到ICRF经历与之相反的四个转换过程,在实际转换中主要采用直接ICRF转化到ITRF的逆矩阵。基于春分点ICRF转化到ITRF转换的经典方法主要流程如下：

在这里插入图片描述
基于CIO的ICRF与ITRF转换的方法
在基于CIO的系统中，由于不依赖于春分点、采用了新的岁差章动模型以及改进极移改正模型，与传统的基于春分点的系统相比，ICRF与ITRF的相互转换的方法与过程有所不同。新转换方法主要包括ICRF转为天球中间参考系，天球中间参考系转换为瞬时地球坐标系，瞬时地球坐标系转换为ITRF三个过程。基于无旋转原点ICRF转化到ITRF转换的经典方法主要流程如下图所示。
在这里插入图片描述

5.2 坐标系统转换实现

目前最新版SOFA，包括了IAU在2006年公布的岁差模型,主要由两部分组成:天文库和矢量矩阵运算库。其中前者有131个子程序,可以进行天文历法计算、时间计算、星历表计算、岁差章动计算、恒星空间运动计算和主要星表系统的转换等,后者有52个子程序,主要功能是矢量和矩阵的各类操作。

5.3 变换所需参数

[GCRS] = Q(t)R(t)W(t)[ITRS]

其中， [ITRS]和[GCRS]分别对应同一位置向量在ITRS和GCRS坐标系中的坐标。W(t)，R(t)和Q(t)分别对应极移矩阵，地球自转矩阵、岁差和章动矩阵，参考IERS2010第五章。这三个转换矩阵大小都是3X3。

计算转换矩阵需要的参数有5个：

t       - modified julian day referenced to UTC  修正儒略日MJD
dt      - UTC-TAI offset (leapsecs)    [seconds] 通过SOFA软件包可以计算
dut     - UT1-UTC offset               [seconds] 通过SOFA软件包可以计算
xp,yp   - coordinates of the pole      [arcsec ] 根据IERS发布的极移参数线性插值（预测）得到MJD对应的xp和yp

其中时间t用于计算dt、dut、xp和yp四个参数。

5.4 GCRS转ITRS的Matlab程序

计算GCRS转ITRS的转换矩阵。计算W(t)，R(t)和Q(t)矩阵的过程在代码中都有所体现。

%% IERS4m 
% These classes facilitates the CIO based celestial to terrestrial transformation. 
% The main function is GCRS2ITRS which provides the 3x3 celestial to terrestrial transformation matrix.  
% In order to obtain the necessary EOP information use helper object USNO.m.  
%
% There is a detailed write up in the docs/latex folder which explains this transformation in detail as well as explains the matrix formulation implemented in MATLAB.  
% All relevant references are listed in the References section.  Select papers are avaliable in the docs/refs folder. 
%
% Default 2010 IERS convention.
%
%% Example Usage
% This example is taken from David Vallado's paper listed in the References section.  

% date time UTC: 2004/04/06 07:51:28.386 mjd = mjuliandate(2004,4,6,7,51,28.386)
fMJD_UTC = 53101.3274118751;

% init EOP object
eopobj = USNO(); 

% pull latest EOP data from USNO servers
% eopobj = eopobj.initWithFinalsHttp(); 

% interpolate EOP information for date and time
[xp,yp,du,dt] = eopobj.getEOP(fMJD_UTC);

% Vallado et al. 2006, AIAA                        [meters]
X_itrs = [-1033.4793830, 7901.2952754, 6380.3565958]';

% dx,dy = 0                                        [meters]
X_gcrs = [5102.5089592, 6123.0114033, 6378.1369247]';

% compute the 3x3 transformation matrix
GC2IT = IERS.GCRS2ITRS(fMJD_UTC,dt,du,xp,yp);

% perform the coordinate/position conversion: C -> T
X = GC2IT*X_gcrs;

% compute the error in meters
err = sum(sqrt((X-X_itrs).^2))

% X = inv(GC2IT)*X_itrs;
% X = GC2IT\X_itrs;
% err = sum(sqrt((X-X_gcrs).^2))

ITRS转GCRS
ITRS转GCRS与GCRS转ITRS的转换矩阵取逆即可。例如前面得到GCRS转ITRS的转换矩阵为GC2IT 。从GCRS转ITRS可以表示为：X_itrs = GC2IT *X_gcrs，相反从ITRS转GCRS就可以表示为：X_gcrs = inv(GC2IT )*X_itrs 。

完整的计算打开可在Matlab论坛（IERS4m）下载。这里涉及到的技术方法在IERS4m可以配合下面参考中的文档进行理解，因此本博文并没有过于详细的记录。这里欢迎大家留言咨询互相学习。