WGS84与CGCS2000坐标的精密转换方法和程序实现
The following article is from 溪流之海洋人生 Author 张祥文 陈正伟
自然资源部规定于2018年6月底前完成各类国土资源空间数据向2000国家大地坐标系的转换,并于2018年7月1日后全面使用2000国家大地坐标系;涉及到空间坐标的报部审查和备案项目,全部采用2000国家大地坐标系。高精度GNSS定位技术已经广泛用于各种测绘工作;其中GPS是主要的定位手段,其定位成果采用WGS84坐标系。经实践比对,目前WGS84坐标系与CGCS2000坐标系在我国华东地区相差达到0.6m左右,且随着时间的推移两者的差异会越来越大;海图测量一般认为“CGCS2000坐标系可等同于WGS84坐标系”,但海洋工程测量的测图比例一般都大于1:2000,不能忽略这一差异,因而需要联测CGCS2000坐标的控制点。
在我国大多数沿海岛屿,联测已知CGCS2000坐标的控制点相对比较困难,因此研究不需要联测CGCS2000坐标控制点的WGS84坐标与CGCS2000坐标的转换方法,实现测量成果转换有重要实际应用价值。
一、WGS84坐标系和CGCS2000坐标系
⒈WGS84坐标系
WGS84坐标系属于世界大地坐标系统,由美国国防部制图局建立。WGS84坐标系采用WGS84椭球,其4个基本椭球参数如下:长半轴a=6378137m,扁率f=1/298.257223563,地心引力常数GM=3.986004418×1014m3/s2,地心自转角速度w=7.292115×10-5rad/s。不同时期的WGS84坐标系所采用的参考框架及其参考历元已经历了4次更新,分别为:
⑴1994年的WGS84(G730)与ITRF1991在1994.0历元处一致;
⑵1997年的WGS84(G873)与ITRF1994在1997.0历元处一致;
⑶2002年的WGS84(G1150)与ITRF2000在2001.0历元处一致;
⑷2012年的WGS84(G1674)与ITRF2008在2005.0历元处一致。
其中,ITRF1991至ITRF2008是基于GPS、VLBI、SLR、LLR和DORIS等空间技术在不同年份建立起来的全球参考框架,也是IGS站坐标和速度场的具体体现。
⒉CGCS2000坐标系
CGCS2000国家大地坐标系于2008年7月1日启用,所采用的4个基本参数如下:长半轴a=6378137m,扁率f=1/298.257222101,地心引力常数GM=3.986004418×1014m3/s2,地心自转角速度w=7.292115×10-5rad/s。CGCS2000坐标框架是利用全球47个IGS核心站的ITRF97框架的坐标和速度矢量,以2000.0历元为参考历元,结合我国GNSS观测数据所建立的参考框架。
⒊两种坐标系的差异
WGS84坐标系与CGCS2000坐标系的4个椭球基本参数中只有扁率有微小差异,由此引起的坐标差异约为0.1mm,在当前测量精度下可忽略。由于WGS84坐标框架经过了4次更新,不同时期WGS84坐标框架与CGCS2000坐标框架之间的差异不能忽略,如2012年后的WGS84(G1674)坐标是ITRF2008框架在2005.0参考历元的坐标,与CGCS2000的差别达到半米以上,必须要进行转换。
二、WGS84坐标与CGCS2000坐标的转换方法
采用各地区建成的CORS系统或千寻位置服务的“千寻知寸—FindCM”高级定位服务,或联测已有CGCS2000控制点进行联合平差,可直接求得CGCS2000坐标;但受到CORS或千寻的服务区域,或者CGCS2000控制点的限制。由于GNSS精密点定位不需要联测任何控制点就能获取精确的WGS84坐标,本文研究如何将观测的WGS84坐标直接转换为CGCS2000坐标的转换方法。
由于观测手段的改进和观测精度的提高,ITRF参考框架也在不断精化,不同时期的ITRF框架之间存在系统性差异。且因地球板块运动,各板块之间和板块内部都存在长期漫长的相对运动,引起同一框架不同历元的坐标也有差异。因此,将WGS84坐标转换至CGCS2000坐标,需要进行参考框架转换和历元改正,即利用不同参考框架之间的转换参数进行参考框架转换,利用板块运动速度场模型进行历元改正。理论上,先转换参考框架再进行历元改正与先改正历元再进行参考框架转换是等价的,实际的数据处理结果也验证了这一观点。不同参考框架之间的转换参数由国际地球自转与参考系统服务(IERS)提供,板块运动速度场模型国内、外学者进行了大量研究,可经过比较测试后采用合适的速度场模型。由于板块运动不仅包含线性运动,也包含非线性运动,随着时间的推移非线性运动的累积误差可能逐渐增大。经算例分析统计,从2000年至今近20年时间,累计误差在华东区域约为1~3cm,在可接受范围内。
三、WGS84与CGCS2000坐标的转换步骤
⒈参考框架转换
IERS已发布了ITRF88-94、ITRF96-97、ITRF2000、ITRF2005和ITRF2008、ITRF2014、ITRF2020等全球参考框架。不同参考框架下的三维空间坐标可采用布尔沙模型进行相互转换,其转换公式如下:
式中:Tx,Ty,Tz和Rx,Ry,Rz分别为x,y,z三个坐标轴的平移参数和旋转参数,D为尺度参数。这些转换参数等于基准历元的参数P(t0)加上历元t0到转换历元的变换量:
由于当前的WGS84坐标是ITRF2008框架在2005.0历元的坐标,CGCS2000坐标是ITRF97框架在2000.0历元是坐标,因此WGS84坐标与CGCS2000坐标的框架转换是ITRF2008与ITRF97框架的转换,其历元差为t-t0=5a。两个框架间的转换参数及其变化速率见表1。
表1 ITRF2008转换到ITRF97框架的转换参数及其速率
⒉板块运动改正
地球不是一个刚体,地球板块会有漂移和形变,板块之间还有挤压、抬升、下降等运动,他们的运动趋势从长期分析是一个非线性非匀速运动,但是从局部和短期内可以把它认为是一种线性匀速运动。板块运动改正即根据板块运动速度计算测站的速度,并依据计算速度将站点坐标从某一历元归算到另一历元。
ITRF框架之间转换,历元不同对转换坐标的影响远远大于框架转换系数的影响,这是因为板块运动导致测站的位置变化,累计到当前已达到分米级。板块运动改正的关键是利用合适的板块运动模型计算出测站所在位置的板块运动速度,若基于欧拉矢量方式(有些学者经平差计算后给出的板块速度是空间直角三维坐标的变化速度)表示板块运动模型,则测站速度计算公式为:
计算得到Vx,Vy,Vz后,就可以进行坐标的历元归算,公式如下:
基于当前历元的观测求解得到的WGS84坐标和CGCS2000的框架历元跨度接近20年,如果没有准确的点位速度场,经上面公式改正的点位误差依然可能达到分米级。常用的速度场模型如下:
⑴NNR-NUVEL1A模型
地质学家根据最近百万年的地质学和地球物理资料,推导出板块运动的平均速度模型,目前国际上推荐使用的是NNR-NUVEL1A板块运动模型,该模型将全球划分为14个板块,我国处于欧亚板块的东部。
NNR-NUVEL1A反应的是大时间尺度上板块的稳定性、刚性运动,其采用的数据在中国也比较少,通过NNR-NUVEL1A模型计算得到的中国大陆速度场残差在E方向和N方向最大值都超过30mm/a,整体RMS也接近10mm/a,说明NNR-NUVEL1A模型只扣除了中国大陆速度场的部分运动趋势,因此不能完全反映中国大陆的整体运动。目前国外通用软件在大多采用该模型,这也是通用软件提供的CGCS2000坐标的最大缺陷。
⑵CPM-CGCS2000模型
CPM-CGCS2000模型是基于中国地壳运动观测网络2001—2010年跨度长达10年的观测数据,采用基准优选和变异点数据分段处理等策略,计算获得ITRF2005框架下高精度速度场,同时针对国际上现有的NNR-NUVEL1A、APKIM2005、PB2002等板块模型在中国区域适应性差,基于中国地质构造特性及实际速度场解算结果,构建了中国20个二级板块运动模型,CPM-CGCS2000板块欧拉矢量及板块拟合误差见表2,与国际上现有的几个成熟的模型相比,在整个中国地区,CPM-CGCS2000相较于现有模型更能精确反映站点的水平运动,并且精度提高了2至5倍。其转化精度优于国际上现有比较成熟的速度场模型,同时该模型也是《大地测量控制点坐标转换技术规范》(CH/T2014-2016)规范规定采用的模型。
表2 CPM-CGCS2000板块欧拉矢量及板块拟合误差
板块 | 欧拉矢量/rad·Ma | 拟合误差/mm·a-1 | ||
Ωx | Ωy | Ωz | σ | |
阿尔泰 | 0.000628 | -0.001876 | 0.004746 | 0.74 |
阿拉善 | 0.000410 | -0.005542 | 0.001580 | 0.82 |
巴颜喀拉 | 0.000242 | -0.006253 | 0.002678 | 2.02 |
柴达木 | 0.001663 | -0.010674 | 0.000900 | 2.86 |
华南 | -0.000936 | -0.002695 | 0.004548 | 1.67 |
川滇 | 0.000616 | -0.016341 | -0.002104 | 3.39 |
滇西南 | -0.001726 | -0.002290 | 0.003626 | 2.38 |
拉萨 | 0.002986 | -0.006392 | 0.003822 | 5.02 |
鲁东海 | -0.001957 | -0.000385 | 0.006133 | 0.69 |
羌塘 | 0.002197 | -0.027764 | -0.008042 | 2.48 |
祁连 | 0.000247 | -0.004865 | 0.002917 | 1.94 |
南海 | 0.000102 | -0.004319 | 0.003061 | 1.50 |
天山 | 0.000856 | -0.002689 | 0.004057 | 1.29 |
中蒙 | -0.000743 | -0.001777 | 0.004637 | 0.88 |
塔里木 | 0.000904 | -0.009939 | -0.002243 | 1.76 |
准噶尔 | 0.000748 | -0.000028 | 0.006592 | 1.31 |
中韩 | -0.001021 | -0.001582 | 0.004737 | 0.98 |
华北 | -0.001083 | -0.001761 | 0.005133 | 0.88 |
鄂尔多斯 | -0.001116 | -0.001303 | 0.005514 | 0.98 |
燕山 | -0.000773 | -0.002084 | 0.004447 | 0.81 |
⑶其它模型
NNR-NUVEL1A模型存在残差过大的缺点,CPM-CGCS2000模型的板块划分,一般用户无法获得点位所在的板块,局部小块体处于同一块体的可能性比较大,只能大致认定所在板块。关于速度场模型的研究,魏子卿等构建了中国大陆地区3°×3°和2°×2°速度场模型,以空间直角坐标变化速度方式给出了格网平均速度,方便设计程序自动判断模型参数。苗岳旺等将中国大陆按照省级行政单位划分块体,求解了31组欧拉矢量,一般用户都能方便的判断测区所在省级行政单位。国外其它模型还有APKIM2005、PB2002等。
图1 CPM-CGCS2000板块划分
四、算例
按照上述的参考框架转换方法和转换步骤,本文编制了WGS84和CGCS2000坐标批量相互转换的计算程序,提供了观测文件(*.RAW)或成果文件(*.XYZ)的批量转换。程序内置了本文介绍的几种速度场模型,方便用户根据具体情况进行选择。经比较,相同控制点利用不同速度场模型进行历元归算,CPM-CGCS2000精度最高,NNR-NUVEL1A最差。
图2 WGS84和CGCS2000坐标批量转换程序
2018年8月在浙江某工程中,采用星站差分SeaStar在已知控制点上进行固定点比对(采集时间1h以上,共4232点),测得的坐标(WGS84坐标)与控制点已知的CGCS2000坐标存在偏差(平均差值东向约0.59m,北向-0.25m),将测量成果经程序转换后再与已知CGCS2000的坐标比较,剩余平均偏差约3.5cm,即转换残差为3.5cm,考虑到板块非线性运动和观测误差的综合影响,该误差完全满足工程的要求。观测点相对于已知点CGCS2000坐标的分布如图3所示。
图3 浙江某工程固定点比对观测点分布图
同时本文下载了IGS站SHAO在ITRF2008框架下2000.0-2019.0每年1月1日的坐标,利用程序经框架转换和历元归算转为CGCS2000坐标,并与该点的CGCS2000已知坐标进行比较,从图4(图中X,Y,Z分别代表三个维度的差异,T为总的位置差)可以看出,残差最大不超过3cm。
图4 WGS84坐标转至CGCS2000坐标成果的残差分布
五、结论
通过参考框架转换与历元归算,实现WGS84坐标与CGCS2000坐标的转换,通过前文分析可得出如下几点结论:⑴转换精度约在几个厘米,能满足生产要求;且框架的影响远小于历元的影响,若精度要求为10厘米级时,不进行框架转换也能满足要求;⑵历元的归算受速度场模型的影响较大,经测试CPM-CGCS2000模型精度最高;国际通用软件一般采用NNR-NUVEL1A模型,在中国大陆地区误差较大,需谨慎使用;⑶本文是基于WGS84坐标当前的框架ITRF2008进行的分析,如WGS84再次进行精化改变了GPS卫星星历的框架,需根据具体情况进行相应的转换。
【作者简介】文/张祥文 陈正伟,均来自上海海事测绘中心。第一作者张祥文,男,1967年出生,本科,高级工程师,主要从事海洋测绘研究;通讯作者陈正伟,男,1981年出生,本科,高级工程师,主要从事海洋测绘研究。文章来自《海洋技术学报》(2020年第6期),参考文献略,用于学习与交流,版权归作者及出版社共同拥有。
投影、坐标系、坐标转换...这个坐标系统培训课件给你讲明白(课件可下载)