郭 旭，陸 港，周 瀟

(空軍預(yù)警學(xué)院，湖北 武漢 430019)

0 引 言

彈道導(dǎo)彈預(yù)警是戰(zhàn)略反導(dǎo)作戰(zhàn)的前提，其發(fā)展水平關(guān)乎國家命脈，因而受到了廣泛重視。由于起步較早，國外軍事強國的預(yù)警體系建設(shè)已較為完善，具備了很強的預(yù)警探測能力[1-2]。要對導(dǎo)彈的預(yù)警過程及預(yù)警參數(shù)進(jìn)行分析，需要使用多種專業(yè)功能突出的仿真軟件[3]。其中STK功能強大，可用于航天、情報、雷達(dá)、電子對抗、導(dǎo)彈防御等各方面[4-5]。我國目前已經(jīng)開始彈道導(dǎo)彈預(yù)警體系的開發(fā)，但由于起步較晚，很多技術(shù)尚處于起步階段，還需要很多探索。

探索研究預(yù)警參數(shù)與導(dǎo)彈各項預(yù)警參數(shù)之間的關(guān)系，從而找到最合理的優(yōu)化配置方法，選擇合理的探測精度，得到最為優(yōu)化的預(yù)警結(jié)果，在一定程度上可以提高導(dǎo)彈預(yù)警系統(tǒng)的總體設(shè)計質(zhì)量，為今后導(dǎo)彈仿真優(yōu)化設(shè)計以及對我國的彈道預(yù)警體系提供參考。

1 導(dǎo)彈落點預(yù)報及誤差分析

通過衛(wèi)星觀測導(dǎo)彈，并將信息提供給雷達(dá)，雷達(dá)就能對目標(biāo)導(dǎo)彈進(jìn)行觀測，得到目標(biāo)導(dǎo)彈在一段時間內(nèi)的位置、速度信息。根據(jù)這些信息，就可以對目標(biāo)導(dǎo)彈的落點進(jìn)行預(yù)報。

1.1 航天器坐標(biāo)系簡單介紹

為了描述運動點的規(guī)律，確定點在空間中的位置，必須確定恰當(dāng)?shù)淖鴺?biāo)系，包括確定坐標(biāo)原點、參考平面和參考平面上的主方向3個要素。這里首先對常用的坐標(biāo)系進(jìn)行以下簡單介紹：

(1) 地心慣性坐標(biāo)系OeXIYIZI

地心慣性坐標(biāo)系是一個慣性參考基準(zhǔn)，用來描述航天器在慣性空間中的絕對運動。地心慣性坐標(biāo)系OeXIYIZI的3個要素定義如下：

(a) 坐標(biāo)原點為地球質(zhì)心。

(b) 參考平面OeXIYI為地球赤道面。

(c) 參考平面上的主方向：OeXI軸指向春分點方向；OeZI軸與地球南北軸重合，指向北極；OeYI軸與OeXI軸和OeZI軸構(gòu)成右手坐標(biāo)系。

(2) 地心固連坐標(biāo)系OeXEYEZE

地心固連坐標(biāo)系是隨地球一同轉(zhuǎn)動的動參考基準(zhǔn)，用以描述航天器相對地球的視運動，即航天器相對地球的相對運動。地心固連坐標(biāo)系OeXEYEZE的3個要素定義如下：

(a) 坐標(biāo)原點Oe為地球質(zhì)心。

(b) 參考平面OeXEYE為地球赤道面。

(c) 參考平面上的主方向：OeXE軸在赤道面指向格林尼治子午線(零子午線)方向；OeZE軸與地球南北軸重合，指向北極；OeYE軸與OeZE軸構(gòu)成右手坐標(biāo)系。

(3) 導(dǎo)彈發(fā)射極坐標(biāo)系o1x1y1z1

導(dǎo)彈極坐標(biāo)系是以導(dǎo)彈的發(fā)射點為參考點，用來描述導(dǎo)彈相對與發(fā)射點的視運動。導(dǎo)彈極坐標(biāo)系o1x1y1z1的3個要素定義如下：

(a) 坐標(biāo)原點o1為導(dǎo)彈在地球表面的發(fā)射點。

(b) 參考平面為與發(fā)射點所在位置相切的平面。

(c) 參考平面上的主方向：o1y1指向?qū)椀陌l(fā)射方向；o1x1與o1z1構(gòu)成右手坐標(biāo)系。

(4) 瞬時地心坐標(biāo)系oxyz

瞬時地心坐標(biāo)系也稱為真地球坐標(biāo)系，它的3個要素定義如下：

(a) 坐標(biāo)原點位于地球質(zhì)心。

(b) 參考平面OeXEYE為地球赤道面。

(c) 參考平面上的主方向：oz軸指向瞬時地球自轉(zhuǎn)軸方向，ox軸指向瞬時赤道面和包含瞬時地球自轉(zhuǎn)軸與平均天文臺參考點的子午面的交點。

(5) J2000慣性坐標(biāo)系oAxAyAzA

由于恒星赤經(jīng)和赤緯會因歲差(與恒星的自行)改變，所以天文學(xué)家們經(jīng)常指定某一特定的紀(jì)元作參考點。J2000坐標(biāo)系是指由J2000時刻的天赤道與二分點來定義天球參考坐標(biāo)系。它是在天文學(xué)上使用的歷元，前綴“J”代表這是一個儒略紀(jì)元法，它指的是儒略日期TT時2 451 545.0，即相對于世界協(xié)調(diào)時間2000年1月1日11∶58∶55.816。

1.2 常見坐標(biāo)系之間的轉(zhuǎn)換

1.2.1 地心慣性坐標(biāo)系、地心固連坐標(biāo)系之間的轉(zhuǎn)換

航天器在地心慣性坐標(biāo)系OeXIYIZI和地心固連坐標(biāo)系OeXEYEZE的三維位置表示是不相同的，但可以相互轉(zhuǎn)換。2種坐標(biāo)系的原點以及OeZI軸與OeZE軸始終是重合的，地心固連坐標(biāo)系OeXE軸和OeYE軸卻始終繞地球南北軸轉(zhuǎn)動，而地心慣性坐標(biāo)系在慣性空間中固定不動。如果知道某初始時刻t0地心固連坐標(biāo)系的OeXE和地心慣性坐標(biāo)系OeXI軸之間的夾角θ=θ0+Ωe(t-t0)，Ωe為地球自轉(zhuǎn)速度，其值約為7.292 115 ×10-5rad/s。

設(shè)在t0時刻之后的任意時刻t，航天器在地心慣性坐標(biāo)系和地心固連坐標(biāo)系中的三維位置用直角坐標(biāo)表示分別為(xI,yI,zI) 和(xE,yE,zE)，兩者之間可以采用下式進(jìn)行相互轉(zhuǎn)換，即：

(1)

式中：θ為此時刻地心固連坐標(biāo)系的OeXE軸與地心慣性坐標(biāo)系的OeXI軸之間的夾角。

1.2.2 地心固連坐標(biāo)系OeXEYEZE與極坐標(biāo)o1x1y1z1的轉(zhuǎn)換

極坐標(biāo)系轉(zhuǎn)換到地心固連坐標(biāo)系，可以先把極坐標(biāo)系o1x1y1z1轉(zhuǎn)換到地心慣性坐標(biāo)系OeXIYIZI，然后再把地心慣性坐標(biāo)系OeXIYIZI轉(zhuǎn)換成地心固連坐標(biāo)系OeXEYEZE。

1.3 導(dǎo)彈落點預(yù)報算法

1.3.1 將雷達(dá)探測到的導(dǎo)彈數(shù)據(jù)進(jìn)行坐標(biāo)系轉(zhuǎn)換

(1) 導(dǎo)彈極坐標(biāo)系o1x1y1z1轉(zhuǎn)換到導(dǎo)彈發(fā)射坐標(biāo)系oxyz

利用STK生成的導(dǎo)彈彈道選擇2個時刻TQ和TP，記錄2個時刻下對應(yīng)的雷達(dá)測量的導(dǎo)彈參數(shù)[6](方位角αQ、αP,俯仰角βQ、βP,徑向距離rQ、rP)，這實際上相當(dāng)于導(dǎo)彈在發(fā)射極坐標(biāo)系下的坐標(biāo)(αQ,βQ,rQ)、(αP,βP,rP)。

通過坐標(biāo)系轉(zhuǎn)換，可以得到導(dǎo)彈在發(fā)射坐標(biāo)系oxyz下的坐標(biāo)(xQ,yQ,zQ)、(xP,yP,zP)。坐標(biāo)系轉(zhuǎn)換過程如下：xQ=rQcosβQcosαQ,yQ=rQsinβQ,zQ=rQcosβQsinαQ;xP=rPcosβPsinαP,yP=rPsinβP,zP=rPcosβPsinαP。

(2) 導(dǎo)彈發(fā)射坐標(biāo)系oxyz轉(zhuǎn)換到地心坐標(biāo)系OeXEYEZE

得到導(dǎo)彈在發(fā)射坐標(biāo)系oxyz下的2點的坐標(biāo)(xQ,yQ,zQ)及(xP,yP,zP)后，繼續(xù)進(jìn)行坐標(biāo)系轉(zhuǎn)化，得到這2點在地心坐標(biāo)系OeXEYEZE下的坐標(biāo)。這2個坐標(biāo)系的轉(zhuǎn)化方法如下：

這里已知雷達(dá)地理坐標(biāo)(λ0,φ0,h)，其中λ0是雷達(dá)方位角，φ0是雷達(dá)的俯仰角，h是雷達(dá)的高度。坐標(biāo)系的關(guān)系如圖1所示。

圖1 坐標(biāo)系轉(zhuǎn)換示意圖

根據(jù)圖1，可以得到2個坐標(biāo)系之間的轉(zhuǎn)化關(guān)系：

(2)

GE=M2[-(90°+α0)]M1[φ0]M3[λ0-90°]

(3)

這樣就得到了導(dǎo)彈在地心坐標(biāo)系下的坐標(biāo)。

(3) 地心慣性坐標(biāo)系OeXIYIZI轉(zhuǎn)換成地心固連坐標(biāo)系OeXEYEZE

將地心坐標(biāo)系下的坐標(biāo)(XE,YE,ZE)轉(zhuǎn)換到發(fā)射瞬時地心坐標(biāo)系下(XI,YI,ZI)，坐標(biāo)轉(zhuǎn)化公式為：

(4)

IE=M2[ΩG+we·t]·M3[ΩG]

(5)

這樣就得到了導(dǎo)彈在發(fā)射瞬時坐標(biāo)系下的坐標(biāo)。

最終要將發(fā)射瞬時地心坐標(biāo)系下的坐標(biāo)(XI,YI,ZI)轉(zhuǎn)換到地心慣性坐標(biāo)系J2000下(XJ,YJ,ZJ)，轉(zhuǎn)換過程需考慮地球自轉(zhuǎn)與攝動的影響。

1.3.2 彈道導(dǎo)彈落點預(yù)報算法

首先以彈道所在平面為坐標(biāo)平面，以地球球心O(橢圓軌道1個焦點)為坐標(biāo)原點，以長軸所在直線為Yd軸建立彈道平面坐標(biāo)系，P、Q兩點為導(dǎo)彈飛行過程中的兩點，B為導(dǎo)彈落點，如圖2所示。

圖2 導(dǎo)彈飛行模型

根據(jù)STK直接讀出P、Q兩點的vP,vQ及rP,rQ，計算相關(guān)軌道根數(shù)和P、Q點坐標(biāo)，計算步驟如下：

計算導(dǎo)彈在J2000坐標(biāo)系下落點時，利用線性空間內(nèi)積不變原理：

(6)

上式可簡化為：

(7)

將對應(yīng)的坐標(biāo)代入，會發(fā)現(xiàn)上式是一個關(guān)于(x1,y1,z1)的三元二次方程，求解即可得到導(dǎo)彈在J2000坐標(biāo)系下的落點坐標(biāo)(x1,y1,z1)。

2 導(dǎo)彈飛行時間的計算及誤差分析

2.1 導(dǎo)彈飛行時間的計算模型

首先以彈道所在平面為坐標(biāo)平面，以地球球心O(橢圓軌道一個焦點)為坐標(biāo)原點，以長軸所在直線為Yd軸建立彈道平面坐標(biāo)系，Xd軸為與Yd軸垂直的坐標(biāo)軸。A為導(dǎo)彈發(fā)射點，B為導(dǎo)彈彈著點，如圖3所示。

圖3 導(dǎo)彈飛行時間計算模型

根據(jù)STK直接讀出B點的vB及rB，計算B點對應(yīng)的真近點角，計算過程如下：

2.2 導(dǎo)彈飛行時間算法的穩(wěn)定性及有效性仿真分析

以上是在理想狀況下，之后對vP,vQ及rP,rQ分別引入均值為0的高斯白噪聲。進(jìn)行100次獨立重復(fù)試驗，計算出每次試驗對應(yīng)的落點,預(yù)報平均值和均方根誤差，并畫出對應(yīng)曲線。

2.2.1 仿真分析流程

圖4 導(dǎo)彈飛行時間算法流程圖

導(dǎo)彈飛行時間算法的穩(wěn)定性仿真驗證步驟如圖4所示,主要包括以下步驟：

(1) 利用STK進(jìn)行雷達(dá)探測導(dǎo)彈仿真，得到導(dǎo)彈在飛行過程中每一時刻的位置、速度信息。

(2) 通過坐標(biāo)系轉(zhuǎn)換，將導(dǎo)彈在極坐標(biāo)系下的位置信息轉(zhuǎn)化到J2000坐標(biāo)系下。

(3) 調(diào)用向量內(nèi)積不變算法，進(jìn)行導(dǎo)彈落點預(yù)報。

(4) 在雷達(dá)的探測數(shù)據(jù)中引入測量誤差，并在不同測量誤差等級下進(jìn)行多次獨立重復(fù)試驗，求得各噪聲等級下落點預(yù)報均值和均方根誤差。

(5) 繪制曲線圖。

(6) 分析落點預(yù)報算法對于測量誤差的穩(wěn)定性。

2.2.2 對導(dǎo)彈的位置測量數(shù)據(jù)加入誤差，驗證算法

首先在導(dǎo)彈的位置測量信息中加入均值為0、方差為σ=0∶50∶500的高斯白噪聲，進(jìn)行100次獨立重復(fù)試驗，計算出每次試驗對應(yīng)的落點預(yù)報平均值和均方根誤差，并畫出對應(yīng)曲線，結(jié)果如圖5所示。

圖5 加入速度測量誤差后的導(dǎo)彈飛行時間預(yù)測效果圖

通過圖6可以看出，加入均值u=0、方差為σ=0∶50∶500(單位m)的位置測量誤差后，導(dǎo)彈預(yù)測飛行時間的均值比較穩(wěn)定，變化不大，說明飛行時間的預(yù)測值對位置測量的穩(wěn)定性較好。

圖6 加入速度測量誤差后的導(dǎo)彈飛行時間均值

通過圖7可以看出，導(dǎo)彈預(yù)測飛行時間的均方根誤差(RMSE)并未隨著誤差等級σ的增大而增大，穩(wěn)定在5 s附近，整體波動幅度較小，說明位置測量誤差對導(dǎo)彈飛行時間的預(yù)測影響較小。對于導(dǎo)彈飛行時間的預(yù)測，實際應(yīng)用中可以放寬對導(dǎo)彈位置測量精度的要求。

圖7 加入速度測量誤差后的導(dǎo)彈飛行時間RMSE

2.2.3 對導(dǎo)彈的速度測量數(shù)據(jù)加入誤差，驗證算法

在導(dǎo)彈的速度測量信息中引入均值為0、方差為0∶50∶500的高斯白噪聲，同樣進(jìn)行100次獨立重復(fù)試驗，計算出每次試驗對應(yīng)的落點預(yù)報平均值和均方根誤差，并畫出對應(yīng)曲線，結(jié)果如圖8所示。

圖8 軟件進(jìn)行導(dǎo)彈飛行時間預(yù)測的效果圖

通過圖9可以看出，加入均值u=0、方差為σ=0∶50∶500 m/s的速度測量誤差后，導(dǎo)彈預(yù)測飛行的均值穩(wěn)定在1 810 s附近，有小范圍的變化。即使速度測量誤差的均值u=0，但導(dǎo)彈預(yù)測落點均值曲線不是一條直線，而是一條曲線，這說明導(dǎo)彈速度測量誤差對于飛行時間的預(yù)測有一定影響。

圖9 加入速度誤差后的導(dǎo)彈飛行時間均值

通過圖10可以看出，導(dǎo)彈預(yù)測飛行時間的RMSE隨著測量誤差的增大而增大，并且增速較快。

與圖7相比，僅在速度測量誤差σ=50 m/s時，導(dǎo)彈預(yù)測落點的RMSE就已達(dá)到34.096 6 s；而當(dāng)測量誤差σ=500 m/s時，導(dǎo)彈預(yù)測落點的RMSE將達(dá)到347.204 4 s。這說明導(dǎo)彈速度測量誤差對最終的導(dǎo)彈飛行時間預(yù)報精度影響較大，在實際應(yīng)用中，應(yīng)盡可能地控制速度測量誤差，使其位于較小的范圍內(nèi)。

圖10 加入速度誤差后的導(dǎo)彈飛行時間RMSE

2.3 落點預(yù)報算法穩(wěn)定性仿真分析

2.3.1 仿真分析流程

導(dǎo)彈落點預(yù)報算法的仿真驗證步驟如圖11所示[7]。

圖11 導(dǎo)彈落點預(yù)報算法流程

主要步驟如下：

(1) 利用STK進(jìn)行雷達(dá)探測導(dǎo)彈仿真，得到導(dǎo)彈在飛行過程中每一時刻的位置、速度信息。

通過坐標(biāo)系轉(zhuǎn)換，將導(dǎo)彈在極坐標(biāo)系下的位置信息轉(zhuǎn)化到J2000坐標(biāo)系下。

調(diào)用向量內(nèi)積不變算法，進(jìn)行導(dǎo)彈落點預(yù)報。

(2) 在雷達(dá)的探測數(shù)據(jù)中引入測量誤差，并在不同測量誤差等級下進(jìn)行多次獨立重復(fù)試驗，求得各噪聲等級下落點預(yù)報均值和均方根誤差。

(3) 繪制曲線圖。

(4) 分析落點預(yù)報算法對于測量誤差的穩(wěn)定性。

2.3.2 對導(dǎo)彈的位置測量數(shù)據(jù)加入誤差，驗證算法

驗證算法穩(wěn)定性的步驟為：

首先用VC調(diào)用STK建立場景[8-9]，場景的起始時刻為“2007年7月1日12∶00”，結(jié)束時刻為“2007年7月2日12∶00”。

隨后建立導(dǎo)彈與雷達(dá)坐標(biāo)系。這里選擇更貼合實際的仿真場景，由于印度始終是我國在西南方向的一個有力對手，其最新研制的“烈火-5”型彈道導(dǎo)彈的打擊范圍已經(jīng)可以覆蓋北京。

因此假設(shè)導(dǎo)彈的發(fā)射點位于印度的新德里，其經(jīng)緯度、海拔為(78.10°，13.37°，0 km)，彈著點假設(shè)

位于我國北京，其經(jīng)緯度、海拔為(116°,39.9°，0 km)。監(jiān)測雷達(dá)位于云南，其經(jīng)緯度、海拔為(103.34°，34.15°，2 km)。至此場景參數(shù)設(shè)定完畢。利用雷達(dá)監(jiān)測導(dǎo)彈，得到來襲導(dǎo)彈在不同時刻的位置、速度信息。通過這些數(shù)據(jù)，調(diào)用落點預(yù)報算法，就可以進(jìn)行導(dǎo)彈落點預(yù)報。

最后，在雷達(dá)所監(jiān)測到的導(dǎo)彈位置信息中分別引入均值為0、方差為0∶50∶500的高斯白噪聲，進(jìn)行100次獨立重復(fù)試驗，計算出每次試驗對應(yīng)的預(yù)測落點平均值和均方根誤差，并畫出對應(yīng)曲線，計算結(jié)果如圖12～圖14所示。

圖12 加入位置測量誤差后導(dǎo)彈落點預(yù)報的效果圖

圖13 加入位置測量誤差后的導(dǎo)彈預(yù)測落點均值

圖14 加入位置測量誤差后預(yù)測導(dǎo)彈預(yù)測落點的RMSE

通過圖13可以看出，加入均值u=0、方差為σ=0∶50∶500(單位m)的位置測量誤差后，導(dǎo)彈預(yù)測落點的均值比較穩(wěn)定，變化不大，說明導(dǎo)彈預(yù)測落點均值對位置測量的穩(wěn)定性較好。

通過圖14可以看出，在位置測量誤差σ<150 m時，導(dǎo)彈預(yù)測落點的RMSE位于0.043～6.3 km范圍內(nèi)，比較穩(wěn)定。并且由正態(tài)分布的特性可知，當(dāng)位置測量誤差σ>150 m時，導(dǎo)彈預(yù)測落點的RMSE將隨著測量誤差的增大而驟然增大，說明當(dāng)測量誤差大于某一等級后，落點預(yù)測算法的穩(wěn)定性開始降低，將無法滿足落點預(yù)測精度。因此應(yīng)該將位置測量誤差σ控制在150 m以內(nèi)。

2.3.3 對導(dǎo)彈的速度測量數(shù)據(jù)加入誤差，驗證算法

在得到的導(dǎo)彈在不同時刻的速度信息中引入均值為0、方差為0∶50∶500的高斯白噪聲，同樣進(jìn)行100次獨立重復(fù)試驗，計算出每次試驗對應(yīng)的落點預(yù)報平均值和均方根誤差，并畫出對應(yīng)曲線，結(jié)果如圖15～圖17所示。

圖15 加入速度測量誤差后的導(dǎo)彈落點預(yù)報效果示意圖

圖16 加入速度測量誤差后的導(dǎo)彈預(yù)測落點均值

圖17 加入速度測量后的導(dǎo)彈預(yù)測落點RMSE

通過圖16可以看出，加入均值u=0、方差為σ=0∶50∶500(單位m/s)的速度測量誤差后，導(dǎo)彈預(yù)測落點的均值有小范圍的變化，說明導(dǎo)彈的落點預(yù)測與測量速度有一定的關(guān)系，即使速度測量誤差的均值u=0，但導(dǎo)彈預(yù)測落點均值曲線不是一條直線，而是一條曲線。

通過圖17可以看出，導(dǎo)彈預(yù)測落點的RMSE隨著測量誤差的增大而增大，并且增速較快。與圖14相比，僅在速度測量誤差σ=50 m/s時，導(dǎo)彈預(yù)測落點的RMSE就已超過10 km，而當(dāng)測量誤差σ=500 m/s時，導(dǎo)彈預(yù)測落點的RMSE將達(dá)到126.10 km。說明導(dǎo)彈速度測量誤差對最終的導(dǎo)彈落點預(yù)報精度影響較大，在實際應(yīng)用中，我們應(yīng)該盡可能地控制速度測量誤差，使其位于一個較小的范圍內(nèi)。

3 結(jié)束語

介紹了導(dǎo)彈的預(yù)警參數(shù)中的彈著點位置及飛行時間這2個參數(shù)的計算模型。對于彈著點進(jìn)行求解計算，本文采用的是線性空間不變法。這種算法，相較于傳統(tǒng)的解析式法來說，計算量小，運算速度塊，實時性好；在導(dǎo)彈的飛行時間計算方面，本文采用的是橢圓軌道的交點法；在介紹了導(dǎo)彈的預(yù)警參數(shù)的計算模型之后，本文又對這種算法的穩(wěn)定性及有效性進(jìn)行了分析。

通過分析得到了不同類別的測量誤差對導(dǎo)彈預(yù)警參數(shù)預(yù)報精度的影響程度。在計算資源有限的情況下，優(yōu)先保證對預(yù)報精度影響較大的誤差的測量精度，從而能夠?qū)ΜF(xiàn)有的預(yù)警過程進(jìn)行優(yōu)化，提高導(dǎo)彈預(yù)警的準(zhǔn)確性。