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

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

CGCS2000大地問題解算工程運(yùn)用研究

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

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

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

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

0　引言

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

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

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

1　大地測量主題正算

1.1定義及理論公式

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

圖1　大地問題解算示意圖

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

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

式中：

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

b——參考橢球短半軸；

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

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

σ——球上的角距，從點(diǎn)P1到點(diǎn)P2；

σ1——球上的角距，從赤道到點(diǎn)P1；

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

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

1.2工程實(shí)踐經(jīng)驗(yàn)

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

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

1.2.2算法修正

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

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

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

1.2.3P2與P1重合的處理

當(dāng)大地線長s為0 m時，P2點(diǎn)與P1點(diǎn)重合，直接給出P2點(diǎn)的坐標(biāo)（B2=B1，L2=L1），大地反方位角A21為A12的反方向角，兩者相差π。

1.2.4穿越極點(diǎn)問題處理

當(dāng)大地方位角A12為π或者0時，有可能會跨越南北極點(diǎn)，一旦s大于P1點(diǎn)與南極點(diǎn)或者北極點(diǎn)的大地線長，大地主題正算無解，因此，需要進(jìn)行特殊處理。

第1步是需要對是否跨越北極點(diǎn)或者南極點(diǎn)進(jìn)行計算和判斷。如果大地方位角A12為0，調(diào)用大地測量主題反算來計算P1點(diǎn)到北極點(diǎn)的大地線長SJ；如果大地方位角A12為π，調(diào)用大地測量主題反算來計算P1點(diǎn)到南極點(diǎn)的大地線長SJ。在計算得到SJ后，如果SJ＜s，則一定會跨越北極點(diǎn)或者南極點(diǎn)，需要執(zhí)行第2步操作；

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

1.2.5大地線長超長問題處理

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

第1步是判斷大地線長是否超長，如果Ss＜s，則需要執(zhí)行第2步操作；

第2步是以P1點(diǎn)為起點(diǎn)，大地方位角為A12，大地線長為Ss，計算出分段終點(diǎn)Ps坐標(biāo)及大地反方位角As1；

第3步是Ps點(diǎn)為起點(diǎn)，大地方位角為As1的反向角，大地線長為s-Ss，計算出P2點(diǎn)坐標(biāo)和大地反方位角A2s，P2點(diǎn)到P1點(diǎn)的大地方位角A21就等于A2s。

1.3驗(yàn)證題目

表2　大地測量主題正算驗(yàn)證題目

2　大地測量主題反算

2.1定義及理論公式

已知P1點(diǎn)坐標(biāo)（B1、L1），P2點(diǎn)的坐標(biāo)（B2、L2），求大地方位角A12、大地反方位角A21、大地線的長度s，叫大地測量主題反算或第2主題［1］，如圖1所示。

GJB6304-2008《2000中國大地測量系統(tǒng)》給出的標(biāo)準(zhǔn)算法［3］如下。

開始迭代計算如下諸式：

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

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

2.2工程實(shí)踐經(jīng)驗(yàn)

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

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

2.2.2歸化經(jīng)差

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

2.2.3點(diǎn)和點(diǎn)重合情況處理

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

2.2.4判斷A12、A21的改進(jìn)算法

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

①計算輔助變量

②判斷起點(diǎn)是否為南極點(diǎn)或北極點(diǎn)，若是，令a1=0；

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

如果a1＞0，那么

如果a1＜0，那么

2.2.5對趾問題處理

2.2.5.1問題提出和解決思路

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

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

2.2.5.2改進(jìn)算法

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

間接方法計算步驟如下：

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

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

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

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

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

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

⑦輸出計算結(jié)果

2.3驗(yàn)證題目

表3　大地測量主題反算驗(yàn)證題目

3　結(jié)論

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

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

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

［3］中國人民解放軍總裝備部.GJB6304-2008 2000中國大地測量系統(tǒng)［S］.北京：中國人民解放軍總裝備部，2008.

［4］史國友，趙慶濤，王玉梅，等.貝塞爾大地主題反解的改進(jìn)算法［J］.交通運(yùn)輸工程學(xué)報，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-），男，江西余江人，碩士。研究方向：指控裝備。