基于CGCS2000椭球的坐标转换程序实现与精度分析

纪志刚，余 锐，宣 伟，郭剑琴

（1.天津市市政工程设计研究院，天津300051 2.武汉大学 测绘学院，湖北 武汉430079）

在许多大型线状工程中，地球曲率引起的长度变形对项目影响显著，需要采用高斯投影的方式对数据进行处理。现有程序多基于克拉索夫斯基椭球或1975国际椭球[1,2]，不利于解算。笔者将CGCS2000椭球参数代入原始方程推出新公式，并将公式编制成程序，利用计算机进行解算，满足了需求。

1 2000国家大地坐标系

2000国家大地坐标系的原点为包括海洋和大气的整个地球的质量中心。坐标系的Z轴由原点指向历元2000.0的地球参考极方向，X轴指向格林尼治参考子午线与地球赤道面（历元2000.0）的交点，Y轴与Z轴、X轴构成右手正交坐标系[3]。采用广义相对论意义下的尺度，2000国家大地坐标系采用的地球椭球为CGCS2000椭球[4],其基本参数如表1所示。

表1 CGCS2000椭球基本参数表[5]

2 程序实现

2.1 程序过程概况

该程序的实现过程如图1所示。

2.2 程序概述

图1 程序运算流程图

该程序基于VS2008平台用C#语言编写，主要模块为6个函数，分别为solve_x0()、conversion_degree()、conversion_Text()、Gauss_ZhengSuan()、Gauss_FanSuan()和Gauss_HuanDai()。

1）Double solve_x0(doubleB)，该函数的功能是利用给定的参数B，去求解赤道到纬度B所对应平行圈的子午线弧长：

式中，M为子午圈曲率半径。将M按级数展开，并进行积分得弧长计算公式为：

将CGCS2000椭球参数代入得各系数数值，如表2所示。

表2 弧长计算公式系数表

将参数代入式（2）得弧长计算公式为：

2）conversion_degree()，该函数的功能是将rad转换为（°）。由于坐标反算的结果为大地坐标，坐标换带时也需要将投影坐标首先反算得到大地坐标。

3）conversion_Text()，该函数的功能是将大地坐标（只考虑经纬度，不考虑高程）转换为rad。程序中大地坐标的输入形式用下例说明：如大地纬度为95°23′18.3″，输入形式为95.231 83。

4）Gauss_ZhengSuan()，该函数的功能是高斯正算，即将已知大地坐标转换为目标带的高斯投影坐标。Gauss_ZhengSuan()所解算出来的大地坐标，其y坐标在原来的基础上增加了500 km。基于CGCS2000椭球元素，适合电算的高斯正算公式为[6]：

式中,

将CGCS2000坐标系参数代入式（7），然后汇编入Gauss_ZhengSuan()函数中。

5）Gauss_FanSuan()，该函数的功能是将高斯投影坐标反算为大地坐标。其所接收的大地坐标的y坐标也是在原来的基础上增加了500 km。Gauss_FanSuan()函数的核心为高斯反算公式，基于CGCS2000椭球元素的高斯反算公式为：

接下来求解Bf。用迭代法，首先令

式中，X为已知投影点X坐标，6 367 449.145 71为式（5）的第1个系数。每次按式（11）迭代计算，直至

6）Gauss_HuanDai()，该函数的功能是高斯换带，即将已知带投影坐标换算为目标带坐标。换带类型不同，所得结果会有所不同。常规的分带有6°带和3°带2种，函数根据分带的不同，又细分为4种类型：6°带换算为 6°带、6°带换算为 3°带、3°带换算为 3°带、3°带换算为6°带[7]。其中6°带和3°带之间的相互转换具有相对复杂性，下面重点介绍二者的相互转换过程。

① 6°带换算为3°带。以图2为例，首先在6°带内利用已知点A的坐标，采用Gauss_FanSuan()函数计算出A点的纬度B以及经度与中央子午线的经差l。若l<1.5°，在3°带内，A点位于以117°为中央子午线的投影带内，此时l不变，再利用Gauss_ZhengSuan()函数计算出A点在该带内的投影坐标；若l>1.5°，在3°带内，A点位于以120°为中央子午线的投影带内，此时l=6°-l，再利用Gauss_ZhengSuan()函数计算出A点在该带内的投影坐标。

② 3°带换算为6°带。首先利用Gauss_FanSuan函数计算出A点的纬度B和经度与中央子午线的经差l。判断A点所在3°带中央子午线是否满足L0=6n-3。若满足，则A点坐标不变；若不满足，l=1.5°+l，再利用Gauss_ZhengSuan函数计算出A点在该带内的投影坐标。程序实例图如图3、图4所示。

图3 程序实例图（1）

图4 程序实例图（2）

3 算例分析

为检验转换结果的准确性，下面均采用Powercoor软件转换结果作参考值，同类型的转换结果进行比较，以获取转换精度。为了同软件中的显示格式相同，各表中经纬度均采用如下表示形式：30°30′15.2″表示为30.301 52。以下表中初始数据均为模拟数据。

表3 正算结果分析表/m

3.1 正反算数据表及结果分析

坐标分量中误差计算公式分别为：

点位中误差为：

对表3各点分别作中误差分析，得出最大点位中误差为0.011 49 m。

表3和表4中，原始大地坐标经过正反算后又得到新的大地坐标，将二者进行比较，结果如表5所示。

表4 反算结果分析表

表5 新旧大地坐标对照表

表3～表5中数据显示，在X方向和B方向转换结果较参考值均偏大，Y方向和L方向转换结果较参考值均偏小，存在一定的系统误差。

3.2 换带计算分析表及结果分析

将坐标（5 635 146.875 0，629 151.817 9）分别以6°→ 6°，6°→ 3°，3°→ 6°，3°→ 3°4 种方式进行换带计算，并将结果与参考值进行比较，见表6。对表6结果同样利用式（11）进行点位中误差分析，经计算得各点转换后最大中误差为1.047 cm。

表6 换带计算表

4 结 语

综上所述，坐标正算结果中误差最大差值可达到1.149 cm,一般情况其精度在mm级；反算得到的大地坐标精度可以达到10-4（″）量级；坐标换带精度一般也可以达到mm级，基本满足了工程建设的需要。

[1]林晓静,张小红.ITRF2005与CGCS2000坐标转换方法与精度分析[J].大地测量与地球动力学,2010,30(2):117-119

[2]唐运海,何诚.坐标转换及换带计算的研究与实验分析[J].测绘与空间信息,2011,34（2）:1-5

[3]李征航,黄劲松.GPS测量与数据处理［M］.武汉:武汉大学出版社,2005

[4]程鹏飞,文汉江.2000国家大地坐标系椭球参数与GRS80和WGS84的比较[J].测绘学报,2009,38(3):89-94

[5]孔祥元,郭际明,刘宗泉.大地测量学基础[M].武汉:武汉大学出版社,2007

[6]刘飞,周琳琳.GPS大地坐标向地方坐标转换的实用方法研究[J].华东师范大学学报,2005(1):73-77

[7]胡圣武.对高斯投影分带的研究[J].地理空间信息, 2012,10(1):54-56