基础测绘控制测量工作中坐标框架与历元的转换

0、 引言

随着2000国家大地坐标系的颁布和推行, 原国家测绘地理信息局要求在现行测绘工作中要以2000国家大地坐标系作为测绘基准, 特别是在国家和省级测绘项目中强制要求使用该坐标系。

基础测绘作为省级测绘项目, 必须使用2000国家大地坐标系。在基础测绘项目的控制测量实施中, 特别是涉及长基线解算时, 通常会遇到收集到的控制点的坐标框架不一致以及起算数据与观测数据历元不一致的问题。本文拟通过编程的方法探索如何快速有效地解决基础测绘控制测量工作中坐标框架与历元不一致的问题。

1、 坐标框架与历元问题

1.1、 坐标框架问题

国际地球参考框架 (ITRF) 是以甚长基线干涉测量 (VLBI) 、卫星激光测距 (SLR) 、全球定位系统 (GPS) 和多里斯系统 (DORIS) 等空间大地测量技术构成的全球大地测量框架。目前, ITRF是国际公认的应用前景最广泛、精度最高的地心坐标系框架。ITRF通过国际地球自转服务 (IERS) 组织分布于全球的跟踪站的坐标和速度场来维持并提供用户使用。IERS每年将全球各站的观测数据进行综合处理和分析, 得到ITRF框架, 以IERS年报和IERS技术备忘录的形式发布。到目前为止, 已发布的ITRF系列有ITRF88、ITRF、ITRF90、ITRF91、ITRF92、ITRF93、ITRF94、ITRF96、ITRF97、ITRF2000、ITRF2005、ITRF2008、ITRF2014共13个框架。

我国2000国家大地坐标系就是定义在ITRF97地心坐标系统中的区域性地心坐标框架。当使用2000国家大地坐标系作为平面坐标基准进行测量时, 需要注意起算点坐标和目标成果数据的坐标框架是否是ITRF97框架。当出现不同于ITRF97框架的坐标时, 要将其他框架的坐标成果统一归算到ITRF97框架下进行使用。

1.2、 历元问题

2000国家大地坐标系作为一种着眼于全球、考虑动态因素的坐标系, 其定义与1980西安坐标系相比有很大的不同。例如, 当我们要用某一观测历元的数据以IGS的站点作为起算点计算待求大地控制点的2000国家大地坐标系的坐标时, 就要将IGS站点的坐标, 从2000国家大地坐标系所定义的历元推算到观测时的历元, 进行平差或者坐标推算, 而后再进行到2000国家大地坐标系所定义的历元的坐标推算。

2、 ITRF框架与历元转换公式

虽然ITRF参考框架是随着时间而变化的, 但无论处于哪个框架哪个历元, 都可以根据一定的方法向任何一个框架及相应的历元进行转换。

2.1、 同一参考框架下不同历元的转换

每个ITRF参考框架都是由空间大地测量观测站的坐标和运动速度定义的, 本文中同一ITRF坐标框架下的不同历元的转换按如下公式计算:

式 (1) 中, Vx、Vy、Vz为观测站在该参考框架下的运动速度;t0为参考历元;t为待求历元。如果该点为IGS跟踪站点, 则其一定历元和一定框架下的坐标成果及速度场数据可以从ITRF官方网站等有关网站上下载。如果该点为待求点, 则其坐标和速度场数据无法从相关网站上获得, 需要通过从相关网站上获得的已知样本点以及其与样本点间的观测数据通过一定的方法插值求得。本文拟采用反距离加权法 (IDW) 计算, 样本点从ITRF官方网站上获得。如公式 (2) 所示。

式 (2) 中, Z' (x0) 为待估值, x0为观测的待估值点;Z (xi) 为区域内位于xi的观测值;Di是待求点距离样各本点的距离, n为参与插值的样本点个数;k为距离的幂。

2.2、 同一历元不同参考框架的转换

同一历元不同ITRF框架之间的转换采用七参数空间相似变换的布尔莎模型实现, 如公式 (3) 所示。

式 (3) 中, 有3个平移参数:Tx、Ty、Tz;3个旋转参数:Rx、Ry、Rz;1个尺度转换因子:D;同时还有七参数的变化率:X1为框架转换前的坐标, X2为框架转换后的坐标。

转换时, 首先应获得转换公式的七参数, 七参数的变化率乘以转换历元与标准历元的历元差得到转换历元下七参数的变化大小, 与七参数标准值相加得到转换公式所需的七参数。其中:

为了求变化率, 对式 (1) 求导并简化, 得到:

这样, 利用式 (3) 和式 (4) , 就可以实现从一个框架到另一个框架指定历元的转换了。

2.3、 不同参考框架不同历元的转换

采用框架和历元分开转换的方法。可以先进行历元转换, 再进行框架转换;也可以先进行框架转换, 再进行历元转换。具体实现同2.2节和2.3节所述的方法。

3、 程序实现及验证

依照上述原理, 用C#语言编制了ITRF框架转换与历元推算专用程序。应用从ITRF官方网站上下载的IGS站点任意历元的多种坐标框架的成果和速度场数据进行验证。

3.1、 框架转换的验证

选取BJFS、URUM、KIT3、POL2、SHAO、LHAZ和WUHN等7个IGS站点作为验证对象, 下载这7个站点2000.0历元的不同框架的坐标成果, 与用程序对框架进行转换的成果进行对比。结果最大差值为-32.2 mm, 中误差为±9.4 mm。

3.2、 有无历元推算时的验证

1) 以BJFS、URUM、KIT3、POL2、SHAO、LHAZ和WUHN等7个IGS站点作为验证对象, 下载这7个站点2000年1月1日和2016年6月14日两个时间点的ITRF2005框架坐标成果和速度场数据, 用程序将2000年1月1日成果用历元推算至2016年6月14日, 而后用推算成果与下载成果进行对比。结果推算坐标成果与下载坐标成果各坐标进行比对, 各坐标分量差值均不超过1 mm。

2) 以BJFS、URUM、KIT3、POL2、SHAO、LHAZ和WUHN等7个IGS站点作为验证对象。采用ITRF官方网站上分别以2001、2003、2006、2009、2012、2013、2014、2016等8个年份中各选取一天下载的观测数据, 以URUM为待求点, 其他6个IGS站点为已知点, 用GAMIT10.6软件进行基线解算, 用GPSNET软件进行平差计算。

当以6个IGS站点的ITRF 97框架、2000.0历元的成果作为起算数据时, 用2001年观测数据解算的基线, 平差后基线向量的残差最小, 其最大值只有4.5 cm;其他年份的观测数据随着年份的增加, 平差后基线残差呈逐年增大的趋势, 用2016年观测数据解算的基线, 平差后基线向量的残差最大, 其最大值已经达到了41.5 cm。

当以6个IGS站点的ITRF 97框架、与观测数据相同历元的成果作为起算数据时, 平差后各坐标分量的基线残差值没有出现随时间变化的趋势, 且各基线向量的残差值都在6 cm之内, 大部分在4 cm之内。

4、 结束语

经过验证, 运用本文所述的原理和程序进行框架转换时, 坐标分量的最大误差不超过3.5 cm, 中误差不超过1 cm, 可以满足C级及以下等级的控制测量起算数据的精度要求。

如不进行历元推算, 直接采用2000.0历元、ITRF 97框架的IGS站点成果作为起算数据进行平差计算, 其误差将随着时间的推移逐渐增大, 到2016年, 其最大值已经达到41.5 cm;如以与观测数据同历元的IGS站点的成果作为起算数据, 其误差可以提高一个数量级。