基于迭代誤差補(bǔ)償?shù)陌l(fā)射點(diǎn)坐標(biāo)精確求解方法

溫永祿，馬 駿，黃 蓓，李 強(qiáng)，許 諾

（中國(guó)運(yùn)載火箭技術(shù)研究院，北京 100071）

0 引言

隨著現(xiàn)代化戰(zhàn)爭(zhēng)的動(dòng)態(tài)化和復(fù)雜化發(fā)展，軍事信息系統(tǒng)將是提升現(xiàn)代化戰(zhàn)爭(zhēng)下指揮員指揮能力、發(fā)揮武器裝備作戰(zhàn)效能的重要抓手［1-2］，而軍事信息系統(tǒng)經(jīng)常涉及到發(fā)射點(diǎn)精確計(jì)算問(wèn)題，該問(wèn)題屬于CGCS2000 坐標(biāo)系下的大地主題解算。對(duì)于小射程的導(dǎo)彈武器，目前通常采用6°帶高斯平面坐標(biāo)，并采用高斯平面坐標(biāo)計(jì)算距離和方向等參數(shù)，但是由于高斯平面坐標(biāo)存在分帶處理，對(duì)于大射程的導(dǎo)彈武器，在計(jì)算距離和方向等參數(shù)時(shí)，與真實(shí)情況存在較大偏差［3-4］。因此，為適應(yīng)新形勢(shì)下的作戰(zhàn)和裝備發(fā)展需求，通常采用大地主題解算算法，解決不同射程下的發(fā)射點(diǎn)坐標(biāo)計(jì)算問(wèn)題。

目前針對(duì)橢球模型下的大地主題解算問(wèn)題，還沒(méi)有可用于直接解算橢球面三角形的解析求解方法。國(guó)內(nèi)外眾多學(xué)者一直關(guān)注該問(wèn)題并提出了多種解算方法，分別適用于短距離、中距離和任何距離。對(duì)于短距離（小于120 km），比較典型是史賴(lài)伯（O. Schreiber）公式和高斯平均引數(shù)公式［5-6］，該公式是在勒讓德（Legendre）的基礎(chǔ)上，利用大地線(xiàn)終點(diǎn)在起始子午圈上的垂足，用以解算大地主題。此后，高斯（Gauss）提出把勒讓德技術(shù)改化成大地線(xiàn)起點(diǎn)和終點(diǎn)的平均緯度和平均方位角為引數(shù)的冪級(jí)數(shù)，這樣勒讓德級(jí)數(shù)式中所有的偶次冪都消失了，可使項(xiàng)數(shù)減少，從而大大加速了級(jí)數(shù)的收斂，形成高斯平均引數(shù)公式，最終達(dá)到大地主題問(wèn)題解算的目的，但是其解算精度與距離有關(guān)，距離越長(zhǎng)，收斂越慢，因此，只適用于較短的距離。對(duì)于中距離（小于400 km）大地主題解算，常用的是巴烏曼投影公式，巴烏曼采用橢球面對(duì)球面的等距離投影方法求解大地主題［7］。針對(duì)中長(zhǎng)距離解算問(wèn)題，貝塞爾（Bessel）提出將橢球面上的參數(shù)按一定條件投影到輔助球面上，然后由邊長(zhǎng)和經(jīng)差投影公式求出輔助球面參數(shù)表示的大地要素［8］，實(shí)現(xiàn)從球面到橢球面上的相應(yīng)轉(zhuǎn)換。貝塞爾大地投影的方法，解算的距離相比高斯平均引數(shù)法要遠(yuǎn)得多，依據(jù)白塞爾的這種解法，派生出許許多多的公式［9-10］，有的是逐漸趨近的解法，有的是直接解法，還有的是其簡(jiǎn)化公式，例如導(dǎo)航使用的大地線(xiàn)長(zhǎng)近似計(jì)算公式，白塞爾公式的直接解法——陳俊勇公式，保持緯度不變的大地投影——張志新公式［11］。韋森特（T.Vincenty）以貝塞爾公式為基礎(chǔ)，推導(dǎo)出嵌套系數(shù)公式［12-14］。1985年，張學(xué)廉以該公式為基礎(chǔ)，采用嵌層約化的方法對(duì)冪級(jí)數(shù)及系數(shù)進(jìn)行改化，推導(dǎo)出3 個(gè)嵌套系數(shù)和兩個(gè)迭代改正數(shù)，進(jìn)而求解距離和經(jīng)差改正項(xiàng)，該方法統(tǒng)稱(chēng)為韋森特公式，也稱(chēng)為嵌套系數(shù)法，可適用于任何距離的大地主題解算，而且計(jì)算量小。韋森特公式能夠適用于線(xiàn)長(zhǎng)從數(shù)厘米到近20 000 km的大地問(wèn)題解算，大地線(xiàn)長(zhǎng)的精度為毫米級(jí)，大地方位角的精度為0.000 1 s，并且我國(guó)于2008 年頒布了GJB6304-2008《2000 中國(guó)大地測(cè)量系統(tǒng)》，推廣了韋森特公式的使用。

針對(duì)軍事信息系統(tǒng)火力籌劃運(yùn)用需求，為實(shí)現(xiàn)對(duì)目標(biāo)的精確打擊，通常需要預(yù)先根據(jù)導(dǎo)彈武器的射程、發(fā)射點(diǎn)方位角（射向）、目標(biāo)點(diǎn)坐標(biāo)，精確計(jì)算發(fā)射點(diǎn)坐標(biāo)，供載機(jī)航線(xiàn)規(guī)劃和作戰(zhàn)規(guī)劃使用，該問(wèn)題是一種非典型的大地主題解算問(wèn)題，目前尚未形成統(tǒng)一的計(jì)算方法解決該問(wèn)題。因此，針對(duì)該問(wèn)題，本文在大地主題正解的基礎(chǔ)上，提出一種基于迭代誤差補(bǔ)償?shù)姆堑湫痛蟮刂黝}解算方法，根據(jù)已知的發(fā)射點(diǎn)方位角、目標(biāo)點(diǎn)坐標(biāo)和射程，初步計(jì)算目標(biāo)點(diǎn)方位角，采用射向迭代法，再調(diào)用大地主題正解模型求得新的發(fā)射點(diǎn)方位角，求取兩次發(fā)射點(diǎn)方位角之間的偏差，利用偏差對(duì)目標(biāo)點(diǎn)方位角進(jìn)行補(bǔ)償修正，進(jìn)一步減小發(fā)射點(diǎn)方位角誤差，直至計(jì)算得到的發(fā)射點(diǎn)方位角與已知的發(fā)射點(diǎn)方位角誤差小于預(yù)設(shè)的極小值，從而解決發(fā)射點(diǎn)坐標(biāo)精確計(jì)算問(wèn)題，該方法計(jì)算適用射程范圍廣、精度高、計(jì)算量小，能夠滿(mǎn)足軍事信息系統(tǒng)工程化運(yùn)用需求。

1 典型大地主題解算方法

假設(shè)發(fā)射點(diǎn)大地經(jīng)緯度為B1，L1，發(fā)射點(diǎn)至目標(biāo)點(diǎn)的大地方位角（射向）為A12和射程S，目標(biāo)點(diǎn)大地經(jīng)緯度為B2，L2，目標(biāo)點(diǎn)至發(fā)射點(diǎn)的大地方位角為A21，各角度及距離示意圖如圖1 所示。

圖1 大地問(wèn)題解算示意圖

CGCS2000 地球橢球參數(shù)表如表1 所示：

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

1.1 大地主題正解

大地主題正解問(wèn)題是根據(jù)發(fā)射點(diǎn)的大地經(jīng)緯度B1、L1和發(fā)射點(diǎn)至目標(biāo)點(diǎn)的大地方位角A12和射程S，推求目標(biāo)點(diǎn)大地經(jīng)緯度B2、L2和目標(biāo)點(diǎn)至發(fā)射點(diǎn)的大地方位角A21。GJB6304-2008《2000 中國(guó)大地測(cè)量系統(tǒng)》［12］給出的標(biāo)準(zhǔn)算法計(jì)算步驟如下：

發(fā)射點(diǎn)（B1、L1）的地心緯度u1正切三角函數(shù)：

從赤道到發(fā)射點(diǎn)連線(xiàn)在球上的角距σ1的正切三角函數(shù)：

從發(fā)射點(diǎn)到目標(biāo)點(diǎn)連線(xiàn)在球上的角距σ 近似值：

依次迭代計(jì)算如下3 個(gè)公式：

直至σ 無(wú)顯著變化之后（變化值小于MIN=10-10），計(jì)算：

目標(biāo)點(diǎn)（B2、L2）的大地緯度B2正切三角函數(shù)：

發(fā)射點(diǎn)（B1、L1）和目標(biāo)點(diǎn)（B2、L2）在輔助球上的經(jīng)差正切三角函數(shù)：

發(fā)射點(diǎn)（B1、L1）和目標(biāo)點(diǎn)（B2、L2）的經(jīng)差ΔL

目標(biāo)點(diǎn)（B2、L2）的大地經(jīng)度L2：

目標(biāo)點(diǎn)（B2、L2）的大地方位角A21的正切三角函數(shù)：

該模型同樣適用于已知目標(biāo)點(diǎn)坐標(biāo)、射程、目標(biāo)點(diǎn)方位角，計(jì)算發(fā)射點(diǎn)坐標(biāo)和發(fā)射點(diǎn)方位角的情況。

1.2 大地主題反解

大地主題反解問(wèn)題是根據(jù)已知發(fā)射點(diǎn)、目標(biāo)點(diǎn)的大地經(jīng)緯度B1、L1、B2、L2，求兩點(diǎn)之間射程S 和正反大地方位角A12、A21。GJB6304-2008《2000 中國(guó)大地測(cè)量系統(tǒng)》［5］給出的標(biāo)準(zhǔn)算法如下。

發(fā)射點(diǎn)（B1、L1）的地心緯度正切三角函數(shù)：

目標(biāo)點(diǎn)（B2、L2）的地心緯度正切三角函數(shù)：

依次迭代計(jì)算如下公式：

發(fā)射點(diǎn)（B1、L1）和目標(biāo)點(diǎn)（B2、L2）之間的射程：

發(fā)射點(diǎn)（B1、L1）的大地方位角A12的正切三角函數(shù)：

目標(biāo)點(diǎn)（B2、L2）的大地方位角A21的正切三角函數(shù)：

當(dāng)兩點(diǎn)接近對(duì)趾時(shí)，反解公式可能無(wú)解。（對(duì)趾?jiǎn)栴}是指兩點(diǎn)的直線(xiàn)連線(xiàn)穿過(guò)參考橢球球心的情況，主要指射程在20 000 000 m 左右，具體解決方法見(jiàn)文獻(xiàn)［8］）。

2 基于迭代誤差補(bǔ)償?shù)姆堑湫痛蟮刂黝}解算方法

根據(jù)射程、射向（發(fā)射點(diǎn)方位角）和目標(biāo)點(diǎn)坐標(biāo)，計(jì)算發(fā)射點(diǎn)坐標(biāo)、目標(biāo)點(diǎn)方位角，該問(wèn)題是非典型的大地正解問(wèn)題［13］，已知目標(biāo)點(diǎn)大地經(jīng)緯度B2、L2和發(fā)射點(diǎn)至目標(biāo)點(diǎn)的大地方位角A12和射程S，推求發(fā)射點(diǎn)大地經(jīng)緯度B1、L1和目標(biāo)點(diǎn)至發(fā)射點(diǎn)的大地方位角A21。具體求解流程如圖2 所示。

1）首先，根據(jù)A12，估算得出初步的A'21；

2）根據(jù)初步計(jì)算的A'21、目標(biāo)點(diǎn)大地緯經(jīng)度B2、L2、射程S，調(diào)用大地主題正算模型，計(jì)算發(fā)射點(diǎn)大地方位角A'12；

圖2 非典型的大地正解計(jì)算流程

4）采用迭代法進(jìn)行求解，循環(huán)計(jì)算以下步驟，直至ΔA 小于給定的極小誤差值；

②根據(jù)更新后的A'21、輸入的目標(biāo)點(diǎn)大地緯經(jīng)度B2、L2和射程S，調(diào)用大地主題正算模型，計(jì)算發(fā)射點(diǎn)大地方位角A'12；

④判斷ΔA' 是否小于給定的極小誤差值（MIN=10-10）；如果滿(mǎn)足，則退出迭代；如果迭代次數(shù)超過(guò)最大迭代次數(shù)（TIMES=400），則逐步放大迭代誤差要求（乘以10 倍，MIN=10-9，10-8，…，10-3）直至迭代誤差不大于10-3；

5）根據(jù)最終計(jì)算得到的A'21、輸入的目標(biāo)點(diǎn)大地緯經(jīng)度B2、L2和射程S，調(diào)用大地主題正算［14］，得出發(fā)射點(diǎn)緯經(jīng)度B1、L1。

3 算例分析

為驗(yàn)證本文提出方法的正確性和適用性，在Win7 i7-6700 3.4 GHz Matlab 2012b 的環(huán)境下進(jìn)行仿真測(cè)試，極小值取值10-10，計(jì)算時(shí)間約0.1 s。為驗(yàn)證算法的正確性，選用國(guó)軍標(biāo)GJB6304-2008《2000中國(guó)大地測(cè)量系統(tǒng)》中的算例進(jìn)行計(jì)算對(duì)比；為驗(yàn)證算法的適用性，在20 000 000 m 的射程范圍內(nèi)，隨機(jī)在全球范圍內(nèi)生成發(fā)射點(diǎn)和目標(biāo)點(diǎn)坐標(biāo)，形成計(jì)算樣例，對(duì)計(jì)算結(jié)果的適用性進(jìn)行分析。

3.1 與GJB 中的算例對(duì)比

根據(jù)GJB6304-2008《CGCS2000 中國(guó)大地測(cè)量系統(tǒng)》中提供的解算算例，以目標(biāo)點(diǎn)經(jīng)度L2、目標(biāo)點(diǎn)緯度B2、射程S、發(fā)射點(diǎn)方位角A12為輸入計(jì)算發(fā)射點(diǎn)經(jīng)度L1、發(fā)射點(diǎn)緯度B1和目標(biāo)點(diǎn)方位角A21，以標(biāo)準(zhǔn)中的算例結(jié)果為參考對(duì)計(jì)算結(jié)果進(jìn)行對(duì)比分析。

通過(guò)與標(biāo)準(zhǔn)算例中的結(jié)果進(jìn)行對(duì)比，計(jì)算偏差如圖3 所示，計(jì)算得到的發(fā)射點(diǎn)坐標(biāo)緯度最大偏差為0.000 333''，經(jīng)度最大偏差為0.000 121''，發(fā)射方位角偏差最大為0.000 196''，滿(mǎn)足發(fā)射點(diǎn)坐標(biāo)的計(jì)算精度要求。

表2 仿真計(jì)算結(jié)果

圖3 與標(biāo)準(zhǔn)算例的計(jì)算偏差

3.2 隨機(jī)測(cè)試算例

在目標(biāo)點(diǎn)坐標(biāo)緯度范圍［-90，90］°，經(jīng)度范圍［0，360）°，射程［0，20 000 000］m，射向［0，360）°范圍內(nèi)，隨機(jī)生成測(cè)試樣例，測(cè)試次數(shù)設(shè)置為5 000 次，對(duì)結(jié)果進(jìn)行統(tǒng)計(jì)與分析。首先，根據(jù)目標(biāo)點(diǎn)坐標(biāo)、射程、射向（發(fā)射點(diǎn)方位角），采用非典型的大地主題正解模型計(jì)算發(fā)射點(diǎn)坐標(biāo)、目標(biāo)點(diǎn)方位角［15］；然后再根據(jù)得到的發(fā)射點(diǎn)坐標(biāo)，結(jié)合目標(biāo)點(diǎn)坐標(biāo)，采用國(guó)軍標(biāo)中的大地主題反解公式計(jì)算射程、發(fā)射點(diǎn)和目標(biāo)點(diǎn)方位角，并與已知射程、發(fā)射點(diǎn)方位角和非典型大地主題解算計(jì)算得到的發(fā)射點(diǎn)坐標(biāo)進(jìn)行對(duì)比，分析計(jì)算誤差如下頁(yè)圖4～圖6 所示。

根據(jù)仿真計(jì)算結(jié)果可看出，在設(shè)定的數(shù)據(jù)取值范圍內(nèi)，隨機(jī)生成的5 000 個(gè)測(cè)試樣例中，射程計(jì)算偏差絕對(duì)值最大為6.08×10-3m，偏差量級(jí)為10-3m；發(fā)射點(diǎn)方位角偏差絕對(duì)值最大為7.23×10-7°m，偏差量級(jí)為10-7°（10-4''），目標(biāo)點(diǎn)方位角偏差絕對(duì)值最大為7.78×10-7°，偏差量級(jí)為10-7°（10-4''），計(jì)算精度能夠滿(mǎn)足發(fā)射點(diǎn)精確計(jì)算要求，而且算法具有很好的魯棒性。

圖4 目標(biāo)點(diǎn)坐標(biāo)分布

圖5 射程偏差

圖6 目標(biāo)點(diǎn)方位角偏差

圖7 發(fā)射點(diǎn)方位角偏差

4 結(jié)論

本文通過(guò)分析發(fā)射點(diǎn)精確計(jì)算特點(diǎn)，在大地主題正算（韋森特公式）、大地主題反算求解方法的基礎(chǔ)上，提出一種基于迭代誤差補(bǔ)償?shù)姆堑湫痛蟮刂黝}正解方法，采用國(guó)軍標(biāo)GJB6304-2008 中的算例和隨機(jī)測(cè)試樣例，對(duì)算法的適用性和精確性進(jìn)行驗(yàn)證，通過(guò)仿真結(jié)果進(jìn)行分析可知，本文提出的方法適用于20 000 km 內(nèi)任意兩點(diǎn)的發(fā)射點(diǎn)精確計(jì)算，方法簡(jiǎn)單實(shí)用，能夠滿(mǎn)足軍事指揮控制信息系統(tǒng)的發(fā)射點(diǎn)坐標(biāo)精確計(jì)算要求。