CGCS2000大地问题解算工程运用研究

徐洪群，何文东，刘斌，孙洪旗

（1.炮兵防空兵装备技术研究所，北京100012；2.西南计算机有限责任公司，重庆400060）

CGCS2000大地问题解算工程运用研究

徐洪群1，何文东2，刘斌2，孙洪旗2

（1.炮兵防空兵装备技术研究所，北京100012；2.西南计算机有限责任公司，重庆400060）

军事信息系统研制中经常涉及CGCS2000坐标系下的大地问题解算，但直接采用韦森特公式进行软件设计存在一些问题。通过对韦森特公式在工程运用中的问题进行了分析，对正解、反解算法进行了补充完善并给出了相关的计算实例供编程检查。改进后的算法简单实用，适合军事信息系统软件的实现。

CGCS2000，大地问题解算，韦森特公式，工程化

0　引言

军事信息系统研制中经常涉及CGCS2000坐标系下的大地问题解算，目前通常采用6°带高斯平面坐标，并采用高斯平面坐标计算距离和方向等参数。由于高斯平面坐标存在分带的问题，不适合长距离大地计算，因此，有必要采用大地坐标方式，并采用大地问题解算算法来适应新形势下的作战和装备发展需求。

大地问题解算分为大地主题正算和大地主题反算，是研制军事信息系统的基本计算。目前由于还没有直接解算椭球面三角形的封闭公式，许多数学家和测量工作者研究得出了50余种解算方法，其中大部分适用于短距离（达120 km），一部分适用于中距离（达400 km），只有几种可适用于长达20 000km的距离［1］。这些算法中，比较典型的有高斯-赫耳默特方法、贝塞尔方法、陈俊勇方法、嵌套系数法、韦森特公式等等。关于这些算法的详细介绍可以参考《椭球大地测量学》、《2000国家大地坐标系实用宝典》、GJB6304-2008《2000中国大地测量系统》等书籍和相关论文。但值得注意的是，部分资料给出的公式存在印刷错误，进行编程运用时需要进行推敲。

我国正在推广CGCS2000坐标系，并颁布了GJB6304-2008《2000中国大地测量系统》，该标准推荐使用韦森特公式进行大地问题解算。韦森特公式适用于线长从数厘米到近20000km，大地线长的精度为毫米级，大地方位角的精度为0.000 1 s［3］。工程设计中直接采用韦森特公式，存在不能实现任意距离的计算以及某些其他特殊情况的计算等问题。因此，本文将在CGCS2000坐标系下对韦森特公式加以详细介绍，并结合工程实际对韦森特公式进行补充和改进，以期弥补原公式的不足，使之更能符合工程化运用需要。

1　大地测量主题正算

1.1定义及理论公式

已知P1点坐标（B1、L1）以及其大地方位角A12、大地线的长度s，求大地线的终点P2的坐标（B2、L2）及其大地反方位角A21，叫大地测量主题的正算或第一主题［1］，如图1所示。

图1　大地问题解算示意图

GJB6304-2008《2000中国大地测量系统》给出的标准算法［3］如下：

直至σ无显著变化之后，计算：

式中：

e——子午椭圆第一偏心率；

b——参考椭球短半轴；

α——在赤道处大地线的方位角；

λ——在辅助球上的经差；

σ——球上的角距，从点P1到点P2；

σ1——球上的角距，从赤道到点P1；

σm——球上的角距，从赤道到点P1和点P2之间的线长中点。

表1　CGCS2000地球椭球参数表［2］

1.2工程实践经验

1.2.1参数范围及单位

计算时所有的角度均以弧度为单位，纬度范围［-π/2，π/2］，经度范围［0，2π］，大地方位角范围［0，2π］，大地线长以米为单位。

1.2.2算法修正

①可能是印刷错误，原算法的式（7）不正确。经编程证实，正确的计算式应为：

②经过大量数据试验，发现求取式（8）的λ（代表P2与P1点在辅助球上的经差）时应结合大地方位角A12进行计算和修正：

③在进行大地反方位角A21计算时，需要进行象限判断和修正：一是需要依据sinA12和tanA21的取值来共同决定A21的象限；二是三角函数起始0°与地球大地反方位角起始0°重合，但是两者递增方向相反，三角函数计算结果沿逆时针方向递增，地球大地反方位角沿顺时针方向递增，因此，需要按照不同象限对三角函数计算结果进行修正。

1.2.3P2与P1重合的处理

当大地线长s为0 m时，P2点与P1点重合，直接给出P2点的坐标（B2=B1，L2=L1），大地反方位角A21为A12的反方向角，两者相差π。

1.2.4穿越极点问题处理

当大地方位角A12为π或者0时，有可能会跨越南北极点，一旦s大于P1点与南极点或者北极点的大地线长，大地主题正算无解，因此，需要进行特殊处理。

第1步是需要对是否跨越北极点或者南极点进行计算和判断。如果大地方位角A12为0，调用大地测量主题反算来计算P1点到北极点的大地线长SJ；如果大地方位角A12为π，调用大地测量主题反算来计算P1点到南极点的大地线长SJ。在计算得到SJ后，如果SJ＜s，则一定会跨越北极点或者南极点，需要执行第2步操作；

第2步是在明确跨越南北极点的情况下，调整P1点位置到南极点或者北极点，将大地方位角A12反向，大地线长调整为s-SJ，再次执行大地测量主题正算，其计算结果为最终结果。

1.2.5大地线长超长问题处理

经过大量数据测试，发现当大地线长s在［19 900 000，20 003 931.46］区间内时，大地主题正算的计算结果错误。因此，在实际运用中，需要采用分段计算方法来解决大地线长超长问题，可以将分段长度Ss设定在19000000m，大地线长超过该长度则采用分段计算方法。计算步骤如下：

第1步是判断大地线长是否超长，如果Ss＜s，则需要执行第2步操作；

第2步是以P1点为起点，大地方位角为A12，大地线长为Ss，计算出分段终点Ps坐标及大地反方位角As1；

第3步是Ps点为起点，大地方位角为As1的反向角，大地线长为s-Ss，计算出P2点坐标和大地反方位角A2s，P2点到P1点的大地方位角A21就等于A2s。

1.3验证题目

表2　大地测量主题正算验证题目

2　大地测量主题反算

2.1定义及理论公式

已知P1点坐标（B1、L1），P2点的坐标（B2、L2），求大地方位角A12、大地反方位角A21、大地线的长度s，叫大地测量主题反算或第2主题［1］，如图1所示。

GJB6304-2008《2000中国大地测量系统》给出的标准算法［3］如下。

开始迭代计算如下诸式：

直至λ无显著变化之后（可以让差值小于0.00000001），计算：

上述各式中的参数含义和取值参见大地测量主题正算部分。

2.2工程实践经验

2.2.1参数范围及单位

计算时所有的角度均以弧度为单位，纬度范围［-π/2，π/2］，经度范围［0，2π］，大地方位角范围［0，2π］，大地线长以米为单位。

2.2.2归化经差

采用式（15）计算经差时，其计算结果位于［-2π，2π］范围，参考椭球上任意两点的经差的绝对值不会超过π，因此，需对式（15）的计算结果进行归化，将经差归化到［-π，π］之间。

2.2.3点和点重合情况处理

当P1点和P2点重合（经度和纬度均相同）时，直接给出大地线长s为0m，大地方位角A12为0rad，大地反方位角A21为0rad。

2.2.4判断A12、A21的改进算法

计算A12、A21时，原算法存在缺陷。一是未针对特殊情况（如两点均在赤道上或起点为极点等）处理，二是对象限的判定未给出可行的办法，因此，实际编程时应改进。结合《贝塞尔大地主题反解的改进算法》一文中提出的算法［4］，改进后的算法如下。

①计算辅助变量

②判断起点是否为南极点或北极点，若是，令a1=0；

③如果a1=0（两点在同一经线上，或穿越极点，或起点在南极点或北极点），那么

如果a1＞0，那么

如果a1＜0，那么

2.2.5对趾问题处理

2.2.5.1问题提出和解决思路

对趾点是指穿过参考椭球球心直线与参考椭球面相交的两点。对趾点的经度差绝对值为π，纬度相反或者均等于0。南极点和北极点为特殊对趾点，其经度差绝对值可以为0～2π区间的任何数值。例如P1点（北纬40°，东经120°）的对趾点为P2点（南纬40°，西经60°）。

当P1点和P2点为对趾点或者接近对趾时，韦森特反解公式可能无解［3］。因此，这种情况下不能直接采用该公式进行计算。由于P1、P2最短大地线之间必定与二者的中间经线有交点Pc（或者两线重叠），因此，可以通过分别求取P1与Pc、P2与Pc的大地线长和大地方位角的思路来实现间接计算。间接计算方法的核心就是寻找Pc的坐标。

2.2.5.2改进算法

为了解决P1点和P2点为对趾点或者接近对趾时，韦森特反解公式可能无解的问题，需要将韦森特公式反解公式编写成子模块供调用，并对迭代计算循环进行了次数限制（最多400次）以避免死循环。当获得P1点和P2点的输入后，直接调用该子模块，如果不能获得计算结果，采用间接方法计算。

间接方法计算步骤如下：

①计算P1点和P2点的中间经线的经度；

②依据P1点纬度B1和P2点纬度B2，计算出中央经线的纬度范围Bm：

③将中间经线沿南北方向4等分，计算出5个等分点（Pm）的坐标，分别记为；

④调用子模块分别求取P1与各等分点的大地线长SP1m和大地方位角AP1m、P2与各等分点的大地线长SP2m和大地方位角AP2m，进而求出P1点经各等分点至P2点的大地线长，分别记为；

⑤在大地线长Sm组中，按以下算法搜索出最小值Smin及其相邻值SL和 SR。并记录中间经线的纬度段：

⑥分别计算Smin与SL、Smin和SR的差值。如果某个差值大于允许值（如0.000 1 m），则将中间经线纬度范围收缩到区间，并跳转到步骤③继续迭代搜索；如果两个差值均不大于允许值，则可认为Smin已是P1、P2间的最短大地线长s，Smin所对应的等分点是Pc点，结束搜索；

⑦输出计算结果

2.3验证题目

表3　大地测量主题反算验证题目

3　结论

本文对GJB6304-2008《2000中国大地测量系统》推荐的韦森特公式工程化过程中的相关问题进行了补充讨论，对算法所未考虑的特殊情况提出了解决措施。本文取得的成果适用于在全球范围内任意两点的大地测量主题正算和主题反算，成果较为简单实用，易于软件实现，目前已经成功运用到某型指控装备研制中。

［1］陈健，晁定波.椭球大地测量学［M］.北京：测绘出版社，1989.

［2］程鹏飞.2000国家大地坐标系实用宝典［M］.北京：测绘出版社，2008.

［3］中国人民解放军总装备部.GJB6304-2008 2000中国大地测量系统［S］.北京：中国人民解放军总装备部，2008.

［4］史国友，赵庆涛，王玉梅，等.贝塞尔大地主题反解的改进算法［J］.交通运输工程学报，2009，9（1）：17-18.

Rescarch of Engineering Application of Geodesic Calculation Based on CGCS2000 Coordinate System

XU Hong-qun1，HE Wen-dong2，LIU Bin2，SUNHong-qi2
（1.Equipment And Technologies Research Institute of FA&ADA，Beijing 100012，China；2.Chongqing Southwest Computer Limited Liability Company，Chongqing 400060，China）

The geodesic calculation based on the CGCS2000 coordinate system is always referred in the research and development of the military information system.However，directly desingning softwares with the Vicente Formula will result in some problems.According to the analyses of the problems of the Vicente Formula exerted in engineering，this paper improves the normal，inverse solution algorithms and gives out the correlative calculation instances for the programming inspection. The improved algorithms are simple and applied，fits the research and development of the Military Information System well.

CGCS2000，solvegeodetic problem，vicenteformula，engineering

TP311.1

A

1002-0640（2016）10-0168-06

2015-08-15

2015-09-22

徐洪群（1966-），男，江西余江人，硕士。研究方向：指控装备。