国家测绘地理信息局2008年6月18日发布公告:我国自2008年7月1日起,启用2000国家大地坐标系(CGCS2000)。2000国家大地坐标系与现行国家大地坐标系转换、衔接的过渡期为8~10年。自然资源部、国家测绘地理信息局2017年3月10日发布《关于加快使用2000国家大地坐标系的通知(国土资发[2017]30号)》,明确规定2018年6月底前完成全系统各类国土资源空间数据向2000国家大地坐标系转换,2018年7月1日起全面使用2000国家大地坐标系,届时,国家测绘地理信息局将停止提供非2000国家大地坐标系下的测绘成果。

从公告和文件可以看出,CGCS2000系统必将代替现阶段地质系统大量使用的WGS84、1954北京(BJ54)、1980西安坐标系(Xi’an80)。WGS84、BJ54和Xi’an80坐标系与CGCS2000坐标系在地球椭球参数、坐标框架及参考历元等方面均存在差异^[1,2,3]。现有地质、物化探成果需要进行坐标转换后,才不影响继续使用^[4,5,6]。各坐标系向CGCS2000坐标系转换,相关的研究较多^{[7,8,9,10,11]},但大多是基本理论的公式推导,由于坐标转换模型较多,且各数据管理软件采用的坐标转换模型均不一致,这为实际使用带来了诸多不便。

ArcGIS系列软件是美国环境研究所(environment system research institute,ESRI)开发的GIS软件,主要应用软件是Desktop,包含ArcMap、ArcCatalog、ArcToolbox等用户界面组件,能够进行数据管理和分析,功能强大。ArcGIS在测绘、地质、物化探成果管理中均有着广泛的应用^[12,13,14]。分析和研究ArcGIS软件的坐标转换方法,将更加便捷有效地将已有的各类成果数据的坐标系转换为CGCS2000坐标系。

建国以来,我国先后使用的坐标系统有BJ54坐标系、Xi’an80坐标系以及随着GPS系统的兴起,广泛使用的WGS84坐标系和我国最新建立的具有现代意义的CGCS2000坐标系。除了CGCS2000坐标系刚刚开始推广,其余坐标系统伴随着我国经济建设的发展,在各行各业得到了广泛的应用。

但上述坐标系统在坐标系类型、定向、椭球参数、参考历元上均有显著不同(见表1),这就导致了各参考椭球之间需要进行严格的转换计算,方能保证点的坐标在不同的坐标系统下正确地显示。

在ArcGIS中,定义坐标系统时,将平面和高程的坐标系统分别定义,即XY坐标系和Z坐标系,这表示在坐标系变换前后,Z值保持不变。根据坐标类型,XY坐标又分为地理坐标系和投影坐标系两类,地理坐标系是以经纬度表示点位置的坐标形式,投影坐标系则是以采用投影计算后的直角坐标值表示点位置的坐标形式,我国大中比例尺地图均采用高斯—克吕格投影,大于1∶5万比例尺地形图采用3度分带投影,1∶5万及以下比例尺地形图采用6度分带投影。

ArcGIS中可以实现同一椭球内地理坐标系与高斯投影坐标系之间以及不同椭球之间任意坐标系的相互转换,只需准确选择待转数据转换前后的坐标系统即可。同一椭球坐标系之间的坐标类型转换有严密的数学公式,计算过程中不需要输入转换参数,而不同椭球坐标系之间的椭球变换需要计算坐标转换参数。坐标转换参数应根据实际采用的转换模型利用数据范围内有两套坐标系的公共点坐标计算得出,不同工区、不同模型计算的坐标转换参数不能混用。

同椭球之间的地理坐标系与投影坐标系之间的变换,采用高斯正反算公式直接计算。同一组坐标数据,ArcGIS中进行同椭球之间的坐标变换得出的结果与常见的其他软件得出的结果完全一致,这里不再赘述。

针对不同椭球之间的坐标转换,ArcGIS软件提供的坐标转换模型较多,包含了我国常用的坐标转换模型,如“COORDINATE_FRAME”和“MOLODENSKY”,即布尔莎七参数模型和莫洛金斯基三参数模型。七参数模型的公式推导过程和三参数的相类似,下面以莫洛金斯基模型(MOLODENSKY)为例,推导该模型坐标转换参数的计算过程。

考虑到ArcGIS中不能定义空间直角坐标系,那么根据广义大地坐标微分公式^[15]:

(1)X=(N+H)cosBcosL,Y=(N+H)cosBsinL,Z=[N(1-e2)+H]sinB。

将式(1)全微分:

(2)dXdYdZ=JdBdLdH+CΔAΔF,

其中:J和C均为系数矩阵,略去旋转参数项和尺度变化项,并整理得:

(3)dBdLdH=-sinBcosLM+H-sinBsinLM+HcosBM+H-sinL(M+H)cosBcosL(N+H)cosB0cosBcosLcosBsinLsinBΔXΔYΔZ+N(M+H)Ae2sinBcosBM(2-e2sin2B)(M+H)(1-F)sinBcosBsinBcosB00-NA(1-e2sin2B)M1-F(1-e2sin2B)sin2BΔAΔF,

(4)dBdLdH=B2L2H2-B1L1H1。

以从BJ54坐标系向CGCS2000坐标系转换为例,式(3)右侧矩阵中的大地纬度B、大地经度L、大地高H、子午圈曲率半径M、卯酉圈曲率半径N、第一偏心率平方e²、长半轴A、扁率F均为BJ54坐标及其对应的椭球参数;ΔA和ΔF为BJ54对应椭球(Krasovsky_1940)和CGCS2000坐标对应椭球的长轴差和扁率差。

式(4)中的B₁、L₁、H₁为公共点的BJ54地理坐标,B₂、L₂、H₂为CGCS2000地理坐标,B和L的单位为度,H的单位为m。一般情况下,BJ54坐标以高斯坐标的形式给出,那么所对应的B和L可以根据高斯反算公式计算得出。考虑到ArcGIS坐标转换前后高程保持不变的特点,BJ54坐标的大地高H可以利用该点的正常高h代替。

根据需要转换的数据范围内3个公共点的两套坐标值,可以列出9个式(3)的方程,采用矩阵运算即可求得其中的3个平移参数ΔX、ΔY、ΔZ;若有多个公共点的坐标值,则利用间接平差原理^[16],将BJ54坐标系下的坐标视为观测值,设WGS84坐标系下的坐标为无误差,由于各点的坐标可视为同精度独立观测值,因此观测值的权P=1,根据式(3)列出误差方程的矩阵形式:

(5)V=BX-L,

利用最小二乘准则求解误差方程的待定参数X,得:

(6)X=(BTB)-1(BTL)。

最终,根据式(6)即可求解得出有多个公共点的坐标值时的坐标转换三参数。

将X代入式(5),求的坐标改正数V,利用改正数进行坐标转换的精度评定:

(7)σ0=VTPVn-t,

其中:n为总观测个数,t为必要观测数。

ArcGIS软件的坐标转换菜单均集合在Arctoolbox的“投影和变换”工具中。按待转换文件的类型分为栅格文件和矢量文件,栅格文件的坐标变换主要涉及“平移”和“投影栅格”两类,“平移”对Z值无影响,而“投影栅格”需要对Z值进行重采样;矢量文件的坐标变换主要指“投影”。不同椭球之间转换时需预先创建自定义地理(坐标)变换方法,即选择地理坐标变换模型。ArcGIS软件中的各坐标系向CGCS2000坐标系转换的流程如图1。

WGS84和CGCS2000同属现代意义的大地测量基准,由表1可以看出,这两个坐标系统之间只有椭球扁率、参考框架和参考历元有些许差别,根据相关研究^[2,11],椭球扁率和参考框架引起的坐标差异基本忽略不计,参考历元不同引起的坐标差异,涉及到板块运动的速度场,根据2000坐标中国大陆速度场的最大值计算,参考历元引起的坐标位移约5 cm/a^[17],即WGS84坐标与CGCS2000坐标最大差0.6 m。所以在进行WGS84至CGCS2000成果转换时,应根据精度要求具体分析。

已有资料中,若是手持式GPS测量所得的点类成果,该类成果坐标系统属于WGS84,若精度要求较低时,可以认为WGS84坐标与CGCS2000坐标一致;若对精度要求较高时,首先需要确定两套坐标之间的北、东坐标位移量,对栅格数据可以采用“平移”工具向CGCS2000转换,矢量数据则投影为高斯直角坐标系后直接加减计算就得出相应的CGCS2000坐标。

BJ54和Xi’an80坐标系均属于参心坐标系,在椭球参数上与CGCS2000有很大的差异。在北纬18°~54°、东经72°~135°的范围内,由BJ54转换到CGCS2000时,北坐标平移量为-22~-74 m,东坐标平移量为-61~+104 m^[7];由Xi’an80转换到CGCS2000时,北坐标平移量为+9~+30 m,东坐标平移量为+104~+117 m^[4]。从上述可以看出,这两类坐标向CGCS2000转换时不能以简单的加减坐标位移量而得出相应的坐标,必须利用坐标转换模型转换。

在ArcGIS中进行BJ54或Xi’an80坐标向CGCS2000转换时,应利用数据范围内具有两套坐标的公共点坐标值,根据式(6)计算得出该区内的坐标转换三参数,若区内有转换参数,则首先应查清转换参数是何种模型计算得出,与ArcGIS软件采用的模型相一致时方可使用。在ArcGIS中创建自定义地理坐标变换方法,将转换参数输入,利用“投影和变换”工具,就可将BJ54或Xi’an80成果准确的转换为CGCS2000成果。

笔者采用2016年青海东昆仑地区实施的某矿产地质调查项目开展1∶5万重力调查工作时获得的测量坐标成果,在ArcGIS软件中进行坐标转换计算,对上述的转换方法进行进一步说明及验证转换效果。

1∶5万重力调查方法对测量精度的要求很高^[18],经规范、系统的联测工作区内的控制点,同时利用国家B级控制点的CGCS2000坐标,约束平差得出其余控制点的CGCS2000坐标,经统计控制网中误差为平面±0.013 m,高程±0.017 m,工区内控制点的分布情况如图2所示。由于重力勘探中,为方便与已有的资料融合拼接,布格重力异常求取过程中参与各项改正的测点坐标数据都一般采用为BJ54坐标、正常高程系统,所以在实际生产中,利用该区内CGCS2000向BJ54坐标转换参数,得出控制点的BJ54坐标,进而在GPS-RTK测量时,将这些控制点的BJ54坐标作为基准站坐标值,即测得所有重力测点的BJ54坐标,精度统计采用同精度检查的方式进行,经统计测点平面位置中误差±0.118m,高程中误差±0.047m。为了进行转换后坐标值对比,同时利用控制点的CGCS2000坐标作为基准站坐标值,测得部分测点的CGCS2000坐标。

利用7个控制点的BJ54坐标和CGCS2000坐标,按照式(3)计算得出BJ54至CGCS2000转换参数,并根据式(7)统计转换参数的精度为±0.035 m。在ArcGIS软件中,输入待转换的BJ54坐标测点,同时,创建一个新地理坐标变换并将上述转换参数输入。即可利用Arctoolbox组件的“投影和变换”工具将上述测点的坐标系转换为CGCS2000坐标。由于ArcGIS不进行高程转换,则将转换得出的测点CGCS2000平面坐标与实测的坐标值进行比较,比较结果见表2,由表2可以看出:

1) ArcGIS的转换误差完全满足我们的工作需求;

2) 较差呈现随机误差的特征,是转换参数精度和测量精度的综合反映,而转换参数的精度取决于计算参数的控制点个数及在工区的分布情况。

笔者主要探讨了ArcGIS软件的坐标系统管理、投影方式设置及坐标转换方法等内容,特别是针对不同椭球坐标系之间的转换,通过理论公式的推导,得出ArcGIS软件中参数的计算方式及参数的精度评定方法,在ArcGIS软件中实现了不同坐标系向CGCS2000坐标系地精确转换。在此基础上,提出了已有地质、物化探成果各类坐标向CGCS2000坐标转换的解决方法,并将该方法应用于实际生产中,取得了理想的效果。同时,希望这些方法能够为需要利用ArcGIS软件进行坐标转换的其他领域人员提供一些参考。