GAM IT/GLOBK软件中的坐标参考系及其相互转换

2017-11-03罗文彬

四川水泥 2017年9期

关键词：天球参考系历元

夏飞罗文彬

(1湖北省交通规划设计院，湖北武汉 430051；2中铁第一勘察设计院集团有限公司，陕西西安 710043)

GAM IT/GLOBK软件中的坐标参考系及其相互转换

夏飞1罗文彬2

(1湖北省交通规划设计院，湖北武汉 430051；2中铁第一勘察设计院集团有限公司，陕西西安 710043)

系统介绍了GAM IT/GLOBK软件中使用的协议天球与协议地球参考系、ITRF参考框架、WGS84坐标系、IGS精密星历参考框架、球面坐标系、站心地平坐标系、大地坐标系、空间直角坐标系以及NEU坐标系及其相互转换关系。

GAM IT/GLOBK；坐标系统转换；ITRF；参考框架

1 国际协议天球参考系和国际协议地球参考系

1.1 协议天球和地球参考系的定义

国际协议天球参考系（ICRS）往往用来描述天体和卫星运动，天体和卫星的星历通常都在此系统中表示，技术人员可以利用GAMIT软件来计算卫星的运动方程，以及在实际运动过程中利用ARC的轨道坐标来进行计算，因此我们通常使用的 J2000协议来进行惯性坐标计算。这种计算方式的定义：在实际过程中可以利用地球质心来进行坐标原点的计算，以及选用2000年1月1日质心力学时（TDB）为标准历元，利用瞬时来计算X轴和Y轴来计算，这样X轴、Y轴Z轴就形成一个左右手结构体系，ICRS和瞬时天球坐标系之间可以相互转换。

国际协议地球参考系（ITRS）的定义就是在进行地球力场的计算时，以及对地球表面点的定位计算的时候常常使用，这种计算方式的定义：通常需要以地球质心作为坐标原点，并控制z轴的方向为协议地极方向（CTP），控制好X轴指向格林尼治子午圈与协议地球赤道的交点，这样就就可以构成Y轴X、Z轴构成左右手体系。ITRS与瞬时地球参考系可以通过极移改正进行相互转换。

1.2 协议天球和地球参考系间的相互转换

ICRS和ITRS之间的转换需借助于瞬时真天球参考系与瞬时真地球参考系的坐标轴指向相同来实现。瞬时真天球参考系与协议天球参考系的差别在于地球自转岁差和章动引起的坐标轴指向不同，转换公式如下[4,5]：

式中，P为岁差旋转矩阵，N为章动旋转矩阵。AZ，Aθ和Aξ为岁差参数，、Δ、ψΔ分别为黄赤交角、交角章动和黄经章动。目前一般采用IAU1976年岁差常数和1980章动模型。

协议地球参考系与瞬时真地球参考系的差异是由于极移引起的，将瞬时真地球参考系转换为协议地球参考系时，IERS的惯例是先绕轴反时针旋转 Py ，

再绕轴旋转 Px ，转换公式为[4,5]：

式中，M 为极移旋转矩阵。

因此，协议天球参考系与协议地球参考系的转换关系为：

式中，E为地球自转矩阵，GAST为真春分点的格林尼治时角。

2 ITRF参考框架和WGS84坐标系

2.1 ITRF框架相互变换

ITRF框架属于一种动态性的地球参开框架，其定义就是利用框架的时间演变基准的明确定义、框架定向、尺度以及远点等来实现的。在不同的时期，框架之间的四个基准分量定义也就存在着一定的差异，从而造成框架之间存在着较小的系统性差异，而这些差异通常都可以利用七个参数来进行表示，而且不同的框架之间能够通过坐标系之间的相似变换来进行转换，其转换公式为：

式中 1T、 2T 、 3T 为平移量， 1R 、 2R 、 3R 为旋转量，D为就是尺度改正因子，这就是ITRFxx框架到ITRFyy框架转换的七个转换参数。其中，任一参数P都会在指定时刻t的值都会与基准历元的参数 0( )P t 相等加上基准历元 0t 并转换到历元的变化量就会有：

P( t) = P( t0) + P˙( t - t0) （公式9）

这样利用（8）式和（9）式就可以完成不同参考框架到指定历元t的坐标转换。

假设现在需要将 ITRFxx框架 t1历元下的基准点坐标转换到ITRFyy框架 t2 历元下的，一般有两种方法：（1）先转框架，后转历元：首先将ITRFxx框架 t1历元的基准点坐标转换到ITRFyy框架 t1历

在同一框架下不同历元间转换时，如果已知基准点的速度为 V I T RFxx，则基准点的坐标计算可按下式进行：

以上讨论的是已知速度场的框架点间的转换关系，对于不知道速度场的非框架点来说，可以采用目前国际上推荐使用的NNR-NUVEL-1A板块运动模型进行近似计算。每个板块的角速度分量可从地球物理模型计算得到[5]，因此测站的速度可表示为：

式中 P i为测站所处的板块，Ω为以年为单位的角速度矢量（单位为rad/a），

2.2 ITRF框架与IGS变换

GAMIT/GLOBK高精度数据处理一般采用国际GNSS服务中心提供的IGS精密星历。IGS精密星历选用的是ITRF参考框架，IGS采用的ITRF框架情况如表3-1所示。ITRF97框架之前IGS使用与ITRF相同的参考框架和参考历元，ITRF97框架后，IGS开始使用自己的ITRF实现，以保持一致性，IGS实现的框架与ITRF的差异在1cm精度范围内。

表3 -1 IGS产品对应的ITRF框架[5]

高精度GNSS测量基准的统一，包括基准点坐标基准统一和基线解算时卫星星历基准统一两方面的内容。假设某一GNSS测量中使用的IGS精密星历所在的参考框架为ITRFxx，参考历元为 t1；基准点所在的参考框架为ITRFyy，参考历元为 t0 。采用“先转框架，后转历元”的方法将基准点所在参考框架和参考历元与 IGS精密星历归化到同一参考框架和参考历元下，即此时再对归化后兼容性较好的基准点坐标施加强约束，这样建立的地心独立坐标系属于IGS精密星历所采用的ITRF参考框架，实现了测量基准的统一。

2.3 ITRF框架与WGS84变换

WGS84坐标系属于协议地球参考系，是美国GPS广播星历和美国国防制图局NGS精密星历的参考基准，最初是基于Transit卫星多普勒数据建立的用于GPS广播星历的地球参考系，后来主要是基于GPS观测数据实现。其定义是：原点为地球质心，Z轴指向BIH1984.0定义的协议地极CTP方向，X轴指向BIH1984.0零子午面与CTP对应的赤道的交点，Y轴与Z轴和X轴构成右手系。WGS84参考框架是由一组全球分布的GPS跟踪站的坐标来具体实现的。WGS84系统经过三次精化后，目前与ITRF框架的站坐标差异在1cm以内，在厘米级精度内可认为二者是同一参考框架[5]。一般来说，WGS84（1150GPS周）实用上被认为等同于ITRF2000；WGS84（873 GPS周）实用上被认为等同于ITRF94；WGS84（730 GPS周）实用上被认为等同于ITRF92。

3 球面坐标系、站心地平坐标系、大地坐标系和空间直角坐标系

3.1 球面坐标系与空间直角坐标系变换

GAMIT基线解算前需配置站坐标L文件，即输入各测站点的先验坐标，包括球面坐标（GAMIT的传统格式）和指定历元下的空间直角坐标和速度场（GLOBK的apr文件的格式）两种形式。利用gapr_to_l程序可将GLOBK的apr文件转换为指定历元的站坐标L文件，ITRF框架坐标的apr文件可从MIT的ftp目录获得。ITRF参考框架下的球坐标与空间直角坐标的转换公式如下：

式中，为地心纬度，为地心经度，r为地心向径。

使用球面坐标系能够简化球面坐标与空间的直角坐标之间的转换，并且不需要大地坐标拔秧的迭代运算，降低了运算的难度和复杂程度，地心经度与大地精度一样，并且其经度与纬度也比较接近，其地心向径的变化几乎被认识是大地高的变化。

3.2 站心地平坐标系与空间直角坐标系变换

GAMIT基线解算结果输出文件为基线约束解q文件（详细基线解）和o文件（简略基线解）以及基线松弛解 h文件，主要包括基线解算过程参数和基线结果及其精度信息。在q文件和o文件中，基线解算结果各个分量及其方差协方差阵是以空间直角坐标系和站心地平坐标系两种形式给出的。站心地平坐标系P-NEU定义为：以测站点P为原点，以P点的法线为U轴，指向天顶为正，以子午线方向为N轴，指北为正，E轴垂直于P点的大地子午面，向东为正，构成一个左手参考系。ITRF参考框架下的站心地平坐标系与空间直角坐标系的转换公式如下：

式中，L为大地经度，B为大地纬度。

需要说明的是，q文件和o文件中分别给出了XYZ和NEU形式基线分量的协方差阵，而NEU基线分量的中误差是由XYZ

基线分量的转换关系以及协方差阵并以误差传播定律为依据来进行计算的。通常都是将站心地平坐标中的基线NEU分量中的误差当做基线高程以及水平方向的误差。但这一观念只能够在基线相对较短时可这样认为，若是基线较长就应对其基线NEU分量精度以及测站点中的精度之间的差别进行充分的考虑。

4 GLOBK结果文件中的NEU坐标

GLOBK的输出文件一般为*.prt和*.org，在给出ITRF参考框架下的空间直角坐标的同时还给出了新的NEU坐标，它与站心地平坐标系定义的NEU有关系但又有不同，这种坐标系类似于平面坐标，属于圆锥投影[6]。为了区分我们加用（G）表示GLOBK软件中给出的NEU坐标系，记为 ( )N G 、 ( )E G 、 ( )U G 。NEU( G )坐标系统定义为：

N( G )：WGS84 椭球长半轴 a 与测站纬度之积，显然它是一段弧长，北纬为正。

式中B是以弧度为单位的测站大地纬度。

E( G )：以二万分之一弧度为最小度量单位，测站所在处最靠近的那条平行圈到起始子午线的平行圈弧长。

式中，L是以弧度为单位的测站大地经度。r0为余纬为θ0时的纬圈半径，其计算公式为：

式中，θ0定义为最接近二万分之一弧度的余纬，其计算公式为：

式中， I nt(⋅)表示取整运算。

因此，如果知道了 GLOBK 软件计算得到的 N EU( G )坐标, 我们尚不能严格反算出该测站的 XYZ 坐标，原因是无法得到精确的大地纬度B（因为有了取整运算）。另外从 N ( G )、 E ( G )、 U ( G )的计算过程可以看出，空间一组测站的 N ( G )、 E ( G )、 U ( G )并不构成一个统一的空间直角坐标系，因此也无法用两个测站的 N ( G )、 E ( G )、U ( G )坐标差通过平方和开方的方法得到这两个测站间的空间直线距离，这是它不同于一般真正的NEU坐标系的地方。N EU( G)基线分量的中误差可由XYZ直角坐标系与大地坐标系的转换关系以及大地坐标系与 N EU( G )坐标系的转换关系根据误差传播定律计算得到。