徐洪群，何文東，劉斌，孫洪旗

（1.炮兵防空兵裝備技術(shù)研究所，北京100012；2.西南計算機有限責任公司，重慶400060）

CGCS2000大地問題解算工程運用研究

徐洪群1，何文東2，劉斌2，孫洪旗2

（1.炮兵防空兵裝備技術(shù)研究所，北京100012；2.西南計算機有限責任公司，重慶400060）

軍事信息系統(tǒng)研制中經(jīng)常涉及CGCS2000坐標系下的大地問題解算，但直接采用韋森特公式進行軟件設(shè)計存在一些問題。通過對韋森特公式在工程運用中的問題進行了分析，對正解、反解算法進行了補充完善并給出了相關(guān)的計算實例供編程檢查。改進后的算法簡單實用，適合軍事信息系統(tǒng)軟件的實現(xiàn)。

CGCS2000，大地問題解算，韋森特公式，工程化

0　引言

軍事信息系統(tǒng)研制中經(jīng)常涉及CGCS2000坐標系下的大地問題解算，目前通常采用6°帶高斯平面坐標，并采用高斯平面坐標計算距離和方向等參數(shù)。由于高斯平面坐標存在分帶的問題，不適合長距離大地計算，因此，有必要采用大地坐標方式，并采用大地問題解算算法來適應(yīng)新形勢下的作戰(zhàn)和裝備發(fā)展需求。

大地問題解算分為大地主題正算和大地主題反算，是研制軍事信息系統(tǒng)的基本計算。目前由于還沒有直接解算橢球面三角形的封閉公式，許多數(shù)學家和測量工作者研究得出了50余種解算方法，其中大部分適用于短距離（達120 km），一部分適用于中距離（達400 km），只有幾種可適用于長達20 000km的距離［1］。這些算法中，比較典型的有高斯-赫耳默特方法、貝塞爾方法、陳俊勇方法、嵌套系數(shù)法、韋森特公式等等。關(guān)于這些算法的詳細介紹可以參考《橢球大地測量學》、《2000國家大地坐標系實用寶典》、GJB6304-2008《2000中國大地測量系統(tǒng)》等書籍和相關(guān)論文。但值得注意的是，部分資料給出的公式存在印刷錯誤，進行編程運用時需要進行推敲。

我國正在推廣CGCS2000坐標系，并頒布了GJB6304-2008《2000中國大地測量系統(tǒng)》，該標準推薦使用韋森特公式進行大地問題解算。韋森特公式適用于線長從數(shù)厘米到近20000km，大地線長的精度為毫米級，大地方位角的精度為0.000 1 s［3］。工程設(shè)計中直接采用韋森特公式，存在不能實現(xiàn)任意距離的計算以及某些其他特殊情況的計算等問題。因此，本文將在CGCS2000坐標系下對韋森特公式加以詳細介紹，并結(jié)合工程實際對韋森特公式進行補充和改進，以期彌補原公式的不足，使之更能符合工程化運用需要。

1　大地測量主題正算

1.1定義及理論公式

已知P1點坐標（B1、L1）以及其大地方位角A12、大地線的長度s，求大地線的終點P2的坐標（B2、L2）及其大地反方位角A21，叫大地測量主題的正算或第一主題［1］，如圖1所示。

圖1　大地問題解算示意圖

GJB6304-2008《2000中國大地測量系統(tǒng)》給出的標準算法［3］如下：

直至σ無顯著變化之后，計算：

式中：

e——子午橢圓第一偏心率；

b——參考橢球短半軸；

α——在赤道處大地線的方位角；

λ——在輔助球上的經(jīng)差；

σ——球上的角距，從點P1到點P2；

σ1——球上的角距，從赤道到點P1；

σm——球上的角距，從赤道到點P1和點P2之間的線長中點。

表1　CGCS2000地球橢球參數(shù)表［2］

1.2工程實踐經(jīng)驗

1.2.1參數(shù)范圍及單位

計算時所有的角度均以弧度為單位，緯度范圍［-π/2，π/2］，經(jīng)度范圍［0，2π］，大地方位角范圍［0，2π］，大地線長以米為單位。

1.2.2算法修正

①可能是印刷錯誤，原算法的式（7）不正確。經(jīng)編程證實，正確的計算式應(yīng)為：

②經(jīng)過大量數(shù)據(jù)試驗，發(fā)現(xiàn)求取式（8）的λ（代表P2與P1點在輔助球上的經(jīng)差）時應(yīng)結(jié)合大地方位角A12進行計算和修正：

③在進行大地反方位角A21計算時，需要進行象限判斷和修正：一是需要依據(jù)sinA12和tanA21的取值來共同決定A21的象限；二是三角函數(shù)起始0°與地球大地反方位角起始0°重合，但是兩者遞增方向相反，三角函數(shù)計算結(jié)果沿逆時針方向遞增，地球大地反方位角沿順時針方向遞增，因此，需要按照不同象限對三角函數(shù)計算結(jié)果進行修正。

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，調(diào)用大地測量主題反算來計算P1點到北極點的大地線長SJ；如果大地方位角A12為π，調(diào)用大地測量主題反算來計算P1點到南極點的大地線長SJ。在計算得到SJ后，如果SJ＜s，則一定會跨越北極點或者南極點，需要執(zhí)行第2步操作；

第2步是在明確跨越南北極點的情況下，調(diào)整P1點位置到南極點或者北極點，將大地方位角A12反向，大地線長調(diào)整為s-SJ，再次執(zhí)行大地測量主題正算，其計算結(jié)果為最終結(jié)果。

1.2.5大地線長超長問題處理

經(jīng)過大量數(shù)據(jù)測試，發(fā)現(xiàn)當大地線長s在［19 900 000，20 003 931.46］區(qū)間內(nèi)時，大地主題正算的計算結(jié)果錯誤。因此，在實際運用中，需要采用分段計算方法來解決大地線長超長問題，可以將分段長度Ss設(shè)定在19000000m，大地線長超過該長度則采用分段計算方法。計算步驟如下：

第1步是判斷大地線長是否超長，如果Ss＜s，則需要執(zhí)行第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中國大地測量系統(tǒng)》給出的標準算法［3］如下。

開始迭代計算如下諸式：

直至λ無顯著變化之后（可以讓差值小于0.00000001），計算：

上述各式中的參數(shù)含義和取值參見大地測量主題正算部分。

2.2工程實踐經(jīng)驗

2.2.1參數(shù)范圍及單位

計算時所有的角度均以弧度為單位，緯度范圍［-π/2，π/2］，經(jīng)度范圍［0，2π］，大地方位角范圍［0，2π］，大地線長以米為單位。

2.2.2歸化經(jīng)差

采用式（15）計算經(jīng)差時，其計算結(jié)果位于［-2π，2π］范圍，參考橢球上任意兩點的經(jīng)差的絕對值不會超過π，因此，需對式（15）的計算結(jié)果進行歸化，將經(jīng)差歸化到［-π，π］之間。

2.2.3點和點重合情況處理

當P1點和P2點重合（經(jīng)度和緯度均相同）時，直接給出大地線長s為0m，大地方位角A12為0rad，大地反方位角A21為0rad。

2.2.4判斷A12、A21的改進算法

計算A12、A21時，原算法存在缺陷。一是未針對特殊情況（如兩點均在赤道上或起點為極點等）處理，二是對象限的判定未給出可行的辦法，因此，實際編程時應(yīng)改進。結(jié)合《貝塞爾大地主題反解的改進算法》一文中提出的算法［4］，改進后的算法如下。

①計算輔助變量

②判斷起點是否為南極點或北極點，若是，令a1=0；

③如果a1=0（兩點在同一經(jīng)線上，或穿越極點，或起點在南極點或北極點），那么

如果a1＞0，那么

如果a1＜0，那么

2.2.5對趾問題處理

2.2.5.1問題提出和解決思路

對趾點是指穿過參考橢球球心直線與參考橢球面相交的兩點。對趾點的經(jīng)度差絕對值為π，緯度相反或者均等于0。南極點和北極點為特殊對趾點，其經(jīng)度差絕對值可以為0～2π區(qū)間的任何數(shù)值。例如P1點（北緯40°，東經(jīng)120°）的對趾點為P2點（南緯40°，西經(jīng)60°）。

當P1點和P2點為對趾點或者接近對趾時，韋森特反解公式可能無解［3］。因此，這種情況下不能直接采用該公式進行計算。由于P1、P2最短大地線之間必定與二者的中間經(jīng)線有交點Pc（或者兩線重疊），因此，可以通過分別求取P1與Pc、P2與Pc的大地線長和大地方位角的思路來實現(xiàn)間接計算。間接計算方法的核心就是尋找Pc的坐標。

2.2.5.2改進算法

為了解決P1點和P2點為對趾點或者接近對趾時，韋森特反解公式可能無解的問題，需要將韋森特公式反解公式編寫成子模塊供調(diào)用，并對迭代計算循環(huán)進行了次數(shù)限制（最多400次）以避免死循環(huán)。當獲得P1點和P2點的輸入后，直接調(diào)用該子模塊，如果不能獲得計算結(jié)果，采用間接方法計算。

間接方法計算步驟如下：

①計算P1點和P2點的中間經(jīng)線的經(jīng)度；

②依據(jù)P1點緯度B1和P2點緯度B2，計算出中央經(jīng)線的緯度范圍Bm：

③將中間經(jīng)線沿南北方向4等分，計算出5個等分點（Pm）的坐標，分別記為；

④調(diào)用子模塊分別求取P1與各等分點的大地線長SP1m和大地方位角AP1m、P2與各等分點的大地線長SP2m和大地方位角AP2m，進而求出P1點經(jīng)各等分點至P2點的大地線長，分別記為；

⑤在大地線長Sm組中，按以下算法搜索出最小值Smin及其相鄰值SL和 SR。并記錄中間經(jīng)線的緯度段：

⑥分別計算Smin與SL、Smin和SR的差值。如果某個差值大于允許值（如0.000 1 m），則將中間經(jīng)線緯度范圍收縮到區(qū)間，并跳轉(zhuǎn)到步驟③繼續(xù)迭代搜索；如果兩個差值均不大于允許值，則可認為Smin已是P1、P2間的最短大地線長s，Smin所對應(yīng)的等分點是Pc點，結(jié)束搜索；

⑦輸出計算結(jié)果

2.3驗證題目

表3　大地測量主題反算驗證題目

3　結(jié)論

本文對GJB6304-2008《2000中國大地測量系統(tǒng)》推薦的韋森特公式工程化過程中的相關(guān)問題進行了補充討論，對算法所未考慮的特殊情況提出了解決措施。本文取得的成果適用于在全球范圍內(nèi)任意兩點的大地測量主題正算和主題反算，成果較為簡單實用，易于軟件實現(xiàn)，目前已經(jīng)成功運用到某型指控裝備研制中。

［1］陳健，晁定波.橢球大地測量學［M］.北京：測繪出版社，1989.

［2］程鵬飛.2000國家大地坐標系實用寶典［M］.北京：測繪出版社，2008.

［3］中國人民解放軍總裝備部.GJB6304-2008 2000中國大地測量系統(tǒng)［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-），男，江西余江人，碩士。研究方向：指控裝備。