基于CGCS2000椭球的子午弧长正算

甘梦仙，胡 鹏，亢瑞红

（1.池州学院 资源环境与旅游系，安徽池州 247000；2.安徽新华学院 土木与环境工程学院，安徽合肥230088；3.安徽大学 资源与环境工程学院，安徽合肥230601）

大地测量外业的成果是以铅垂线和大地水准面作为基准线和基准面，内业计算时需要的参考基准是法线和参考椭球面。因此，要得到内业的成果，需要将外业测量数据计算和转换，这一工作称为大地测量计算。子午线弧长在测量计算时经常会用到，如研究高斯投影计算、弧度测量计算。

2000国家大地坐标系是我国新一代地心坐标系统，自2008年开始启用。目前的大地测量学和控制测量学等相关教科书上，椭球子午弧长等公式、算例均是基于Krassovsky椭球和IAG-75国际椭球参数。原有的教学内容已不能满足2000国家大地坐标系下椭球面上的测量计算要求。因此，为了与原有坐标承接，确保测绘教学内容实时更新，给出新坐标系下的子午弧长正反算公式十分有必要。李海祥[1]将子午线弧长积分公式按e2展开至e16项，实现高精度的子午线弧长解算；刘仁钊[2]根据二项式定理及三角函数积分递推关系导出弧长积分任意精度的公式；过家春等[3]采用椭圆第二类积分方法实现了弧长的解算，完善了子午线弧长理论。本文在国家高精度子午弧长正算公式的基础上，给出的新坐标系下子午弧长正算实用公式，原理简单、计算简便。计算结果与文献[1]和文献[2]（以下简称迭代算法和递归算法）方法的结果进行比较，得出结论。

1 椭球基本参数简介

2000国家大地坐标系的地球椭球参数与克拉索夫斯基椭球、IAG-75国际椭球、WGS-84椭球体的参数列表如下：

表1 常用地球椭球参数

2 子午弧长正算实用公式

根据国家测绘局《现有测绘成果转换到2000国家大地坐标系技术指南》（下面简称“指南”）的规定，子午弧长计算公式各系数精度有所变动。

设有子午线上一点p，其大地纬度为B，则该点至赤道的子午线弧长X计算公式：

式中

为方便编程计算，令t=tanB，将上式倍角函数化简为幂级数：

其中：

第一偏心率是小于1的数，上式e的10次项以上对结果的影响小于0.0001m，因此可将上式保留8次项即可：

上式（3）为CGCS2000椭球子午弧长计算的实用公式。

3 实例计算与分析

将格拉索夫斯基椭球参数代入式（3），以MATLAB编程计算，结果与已有两种算法（迭代与递归）作对比分析[5-8]，解算结果如表2。

表2 子午弧长计算结果

由上表可以看出，以子午线弧长实用公式计算的子午线弧长可保留0.001项的精度，可达到mm级，满足工程控制测量的精度要求。MATLAB编程计算时，在同等运算精度条件下，e的10次及以上的项对结果的精度影响可以忽略不计。为便于学生手工计算或电算，将适于CGCS2000的子午线弧长实用计算式表示如下：

我国的大地测量计算中还会用到原来采用过的椭球参数，以下是控制测量学教科书中未列出CGCS2000椭球的实用计算式：

WGS-84椭球子午线弧长计算式与我国CGCS2000计算式在主项上基本一致，本文不列出；如有需要计算的，请参照CGCS2000计算式解算。

4 结论

计算表明，该方法简便实用，方便学生和测绘工作者手算和电算，可作为教学和学习的参考。文中特色之处在于将子午弧长正算的精密公式改进为方便实用的计算式，并补充了相关教材缺少的CGCS2000子午弧长计算的示例。至于子午弧长的反算及其在高斯投影垂足纬度的计算问题将是下一步研究的内容。

[1]李海祥，张伟国.高精度任意元素椭球面子午线长度的正反算[J].测绘与空间地理信息，2013（5）:179-181，184.

[2]刘仁钊，伍吉仓.任意精度的子午线弧长递归计算[J].大地测量与地球动力学，2007（5）:59-62.

[3]过家春，赵秀侠，徐丽，等.基于第二类椭圆积分的子午线弧长公式变换及解算[J].大地测量与地球动力学，2011（4）:94-98.

[4]孔祥元，郭际明，刘宗泉.大地测量学基础[M].武汉:武汉大学出版社，2010：114-120．

[5]覃军，骆奇峰.在Excel中实现统一椭球的高斯投影正反算[J].测绘与空间地理信息，2009（2）:205-208.

[6]赵英志，刘永涛，郑玉军.利用VB6.0实现2000国家大地坐标系高斯正反算程序的编写[J].测绘通报，2010（5）:38-41.

[7]刘修善.计算子午线弧长的数值积分法[J].测绘通报，2006（5）:4-6.

[8]过家春.基于第二类椭圆积分的子午线弧长反解新方法[J].大地测量与地球动力学，2012（3）:116-120.