葛瑞星,周青松,張劍云,諶詩娃

(國防科學(xué)技術(shù)大學(xué) 電子對抗學(xué)院,安徽 合肥 230037)

基于迭代凸優(yōu)化的導(dǎo)彈探測時間估算方法

葛瑞星,周青松,張劍云,諶詩娃

(國防科學(xué)技術(shù)大學(xué) 電子對抗學(xué)院,安徽 合肥 230037)

為準(zhǔn)確估算遠(yuǎn)程預(yù)警雷達(dá)對彈道導(dǎo)彈探測的時間窗口,需要聯(lián)立導(dǎo)彈軌跡與雷達(dá)探測范圍的數(shù)學(xué)模型,求解非線性二次方程組。提出了一種非線性二次方程組的求解方法,首先將核問題轉(zhuǎn)化為二次等式的優(yōu)化問題,引入誤差范數(shù)約束非凸優(yōu)化問題,將二次等式進(jìn)行一階和二階泰勒展開,將原問題在局部近似為一個最優(yōu)化問題并利用迭代凸優(yōu)化求解,得到雷達(dá)探測范圍與導(dǎo)彈彈道的交匯點坐標(biāo)的解析解,實現(xiàn)對雷達(dá)探測時間窗口的準(zhǔn)確估算。通過場景模擬與求解分析證明了該方法的正確性,并表明該方法具有精度高、通用性強(qiáng)的特點。

彈道導(dǎo)彈;探測;雷達(dá)探測窗口;非線性二次方程組;迭代凸優(yōu)化

在彈道導(dǎo)彈攻防中,為提高彈道導(dǎo)彈突防成功率、減少導(dǎo)彈防御系統(tǒng)預(yù)警時間,常在彈道導(dǎo)彈上布置電子戰(zhàn)裝備以干擾預(yù)警雷達(dá)探測[1]。導(dǎo)彈在飛行過程中為規(guī)避探測常處于電磁靜默狀態(tài),延遲或滯后釋放電磁干擾易暴露目標(biāo),故估算預(yù)警雷達(dá)對彈道導(dǎo)彈探測的時間窗口,選擇合適的時間窗口對遠(yuǎn)程預(yù)警雷達(dá)進(jìn)行干擾具有重要的研究意義[2]。

近年來,圍繞導(dǎo)彈彈道模型與雷達(dá)探測模型的研究不在少數(shù),文獻(xiàn)[3～5]等研究了簡化的彈道導(dǎo)彈的通用性模型,文獻(xiàn)[6]研究了常見導(dǎo)彈發(fā)射參數(shù)對導(dǎo)彈彈道的影響,文獻(xiàn)[7]根據(jù)彈道理論建立了雷達(dá)探測彈道導(dǎo)彈的仿真模型,文獻(xiàn)[8]建立了探測雷達(dá)對彈道目標(biāo)跟蹤的數(shù)學(xué)模型并進(jìn)行了仿真,文獻(xiàn)[9]結(jié)合SPY雷達(dá)具體特性,分析了雷達(dá)對導(dǎo)彈目標(biāo)的探測效能。但這些研究很少從宏觀角度將導(dǎo)彈彈道與雷達(dá)探測的數(shù)學(xué)模型進(jìn)行結(jié)合研究探測時間估算的通用方法。在求解非線性方程時,一般采用迭代最小二乘法、梯度下降法、擬牛頓法等方法[10],這類方法容易導(dǎo)致解的發(fā)散,缺乏魯棒性。本文將雷達(dá)對導(dǎo)彈的探測時間估算問題轉(zhuǎn)化為非線性方程組求解的數(shù)學(xué)問題,通過迭代凸優(yōu)化的算法進(jìn)行求解,算法具有較高的精度與較強(qiáng)的通用性。

1 探測模型建立與場景分析

估算遠(yuǎn)程預(yù)警雷達(dá)對彈道導(dǎo)彈的探測時間,首先要分析雷達(dá)對導(dǎo)彈的探測場景,在綜合考慮導(dǎo)彈軌跡模型與雷達(dá)探測范圍模型的基礎(chǔ)上進(jìn)行分析。

1.1 坐標(biāo)系建立

本文中導(dǎo)彈發(fā)射點與落點的選定、彈道軌道根數(shù)的求解和模型建立均在地心慣性坐標(biāo)系(ECI)下進(jìn)行,該坐標(biāo)系以地心O為原點,zE軸沿地球自轉(zhuǎn)軸指向正北方向,xE軸為赤道面與本初子午面交線方向由地心向外,yE軸與xE、zE軸構(gòu)成右手坐標(biāo)系。

為簡化導(dǎo)彈軌道模型與雷達(dá)探測模型聯(lián)立后的解析表達(dá)式,定義坐標(biāo)系Oxyz(I系),該坐標(biāo)系以地心O為原點,以發(fā)射點A、目標(biāo)點B連線的中垂線為y軸,方向由地心向外;z軸垂直于OAB所在平面,方向與OA、OB外積方向相同;x軸與y、z軸構(gòu)成右手坐標(biāo)系。I系示意圖如圖1所示。

1.2 探測模型的簡化

由于本文主要解決非線性方程組的求解問題,不失一般性,可對導(dǎo)彈運(yùn)動模型與雷達(dá)探測模型作如下簡化:根據(jù)彈道導(dǎo)彈的飛行時間較短,可忽略其軌道攝動與其自由段運(yùn)動時間長的特點,將導(dǎo)彈運(yùn)動模型近似為繞地球質(zhì)點的二體橢圓運(yùn)動[11];雷達(dá)對彈道導(dǎo)彈進(jìn)行探測時,在忽略導(dǎo)彈運(yùn)動過程中目標(biāo)雷達(dá)截面積(RCS)S起伏效應(yīng)的情況下,可認(rèn)為雷達(dá)對其探測的最大距離保持不變,探測范圍可近似看作一個以雷達(dá)為中心的球;遠(yuǎn)程預(yù)警雷達(dá)大多為地基或海基雷達(dá),部署海拔較低,故可將雷達(dá)的探測視線范圍近似看作過雷達(dá)部署點與地球相切的平面。

1.3 探測場景分析

在上述簡化的基礎(chǔ)上,將導(dǎo)彈軌跡看作I系Oxy平面內(nèi)的橢圓曲線;雷達(dá)探測范圍看作R3中的球體;雷達(dá)視線看作與地球相切的平面,在彈道平面Oxy內(nèi)對雷達(dá)的探測場景進(jìn)行討論。

1)無法被探測的場景。

在平面Oxy中,雷達(dá)探測范圍球在該平面內(nèi)的截面為一個圓。易得當(dāng)該圓與導(dǎo)彈軌跡相離或相切時,導(dǎo)彈無法被雷達(dá)探測,如圖2、圖3所示。

2)能夠被探測的場景。

當(dāng)雷達(dá)探測范圍與彈道軌跡存在交點時,位于交點之間的彈道軌跡能被探測,考慮雷達(dá)視線對探測范圍的遮擋效應(yīng)與發(fā)動機(jī)脫離對探測范圍的影響,對雷達(dá)探測場景進(jìn)行討論。

②發(fā)動機(jī)脫離對雷達(dá)探測的影響。 導(dǎo)彈在主動段結(jié)束時發(fā)動機(jī)脫離,彈頭繼續(xù)按橢圓軌跡運(yùn)動,此時導(dǎo)彈的雷達(dá)截面積(RCS)發(fā)生較大變化,雷達(dá)探測的最大距離改變,對應(yīng)探測球半徑發(fā)生改變,此時雷達(dá)探測模型在Oxy平面中截面是兩段圓心相同的圓弧,記發(fā)動機(jī)脫離點為P點,脫離前探測圓弧為l1,脫離后為l2,將兩段圓弧與導(dǎo)彈軌跡相交的情況分為以下幾類進(jìn)行討論。

2 探測時間估算的數(shù)學(xué)模型

估算遠(yuǎn)程預(yù)警雷達(dá)對彈道導(dǎo)彈的探測時間問題等效于計算導(dǎo)彈進(jìn)入與離開雷達(dá)探測范圍的時刻,該問題的求解需要得到上述關(guān)鍵點的具體位置與導(dǎo)彈運(yùn)動的時間-位置參數(shù)方程。

為方便模型化簡計算,在自定義坐標(biāo)系(I系)下討論導(dǎo)彈運(yùn)動模型與雷達(dá)探測模型。

簡化后的彈道軌跡在自定義坐標(biāo)系I系中可看作三維橢球在彈道平面上的截線軌跡:

(1)

視線平面在I系中表現(xiàn)為一個與地球相切的平面,被Oxy平面所截得的直線表示為

(2)

雷達(dá)探測范圍模型是R3中一個球體,表示為

(3)

式中:X∈R3為三維空間內(nèi)任意點;Rmax為雷達(dá)最大作用距離,參考相控陣?yán)走_(dá)搜索方程[12]:

(4)

式中:Pav為雷達(dá)平均功率;d為雷達(dá)孔徑;ts為雷達(dá)搜索時間;S為目標(biāo)雷達(dá)截面積;L為衰減因子;Ni=kTeBnFn為接收機(jī)輸出的等效噪聲功率;Ωs為搜索空域立體角;D0為檢測因子,表示雷達(dá)檢測所需的最小信噪比。

考慮到導(dǎo)彈在主動段結(jié)束時發(fā)動機(jī)脫離,此時導(dǎo)彈的RCS發(fā)生較大變化,雷達(dá)探測最大距離改變。對于工作頻率10 GHz的雷達(dá),分離前的導(dǎo)彈彈體的RCS量級約為1 m2;分離后的導(dǎo)彈彈頭RCS量級約為0.1 m2。記導(dǎo)彈發(fā)動機(jī)脫離點坐標(biāo)XP=(xPyP0)T,對應(yīng)的雷達(dá)探測距離分別記做R1max和R2max。故在考慮發(fā)動機(jī)脫離時,雷達(dá)探測范圍模型表示為

(5)

將計算導(dǎo)彈在雷達(dá)探測范圍內(nèi)的邊界點問題轉(zhuǎn)化為一系列方程組的求解問題。

計算視線平面與導(dǎo)彈軌道切點:

(6)

計算發(fā)動機(jī)脫離前后探測球與導(dǎo)彈軌跡交點:

(7)

(8)

解得各交點后,結(jié)合1.3節(jié)分析的場景求得實際的探測范圍。

得到導(dǎo)彈在雷達(dá)探測范圍內(nèi)的邊界點后,代入導(dǎo)彈運(yùn)動的時間-位置參數(shù)方程求解雷達(dá)探測時間,本文采用常見的橢圓軌道根數(shù)模型進(jìn)行計算。

在ECI坐標(biāo)中橢圓軌道模型可根據(jù)軌道根數(shù)(a,e,i,ψ,ω,tP)確定,其中,a為軌道長半軸;e為偏心率;i為軌道傾角,表示軌道平面與赤道平面的夾角;ψ為升交點角距,表示升交點與ECI系中x軸的夾角;ω為近地點角距,表示升交點與近地點之間的夾角;tz為導(dǎo)彈起飛至近地點用時。

確定軌道根數(shù)后,可以計算給定時刻t的偏近地點角E(t),進(jìn)而求得該點的地心距離r(t)。

(9)

式中:μ為引力常數(shù)與地球質(zhì)量的乘積。在已知探測范圍邊界點位置后可由上式反推導(dǎo)彈經(jīng)過該點的時刻。

3 迭代凸優(yōu)化求解算法

從上一節(jié)中的分析不難看出,估算雷達(dá)預(yù)警時間的關(guān)鍵在于式(6)～式(8)的非線性方程組求解問題,本文采取的求解方式是將非線性方程的定值求解問題轉(zhuǎn)化為變量優(yōu)化問題進(jìn)行計算。

3.1 彈道軌跡與雷達(dá)探測范圍交點的計算

(10)

該優(yōu)化問題中,約束條件中表示的三維曲面上的點集是非凸集合,由于三維曲面在任意點附近可用過該點的超平面進(jìn)行近似逼近,將其局部松弛為平面后可對優(yōu)化模型采取凸優(yōu)化的算法求解。

對三維曲面函數(shù)進(jìn)行泰勒展開,取其一階近似:令X=X0+ΔX,f(X)在X=X0處的Taylor展開為f(X)=f(X0)+(f(X0))TΔX+o(ΔX),其中f(X0)為f(X)在X0處的微分;o(ΔX)為在ΔX→0時比ΔX高階的無窮小。故一階近似后擬合的超平面方程為f*(X)=f(X0)+(f(X0))TΔX。

不難看出,f*(X)對f(X)的逼近程度與ΔX的取值大小呈負(fù)相關(guān),在選取的X0點與所求點相差較大時,一次近似后的運(yùn)算精度有限,考慮采用迭代算法進(jìn)行多次迭代優(yōu)化求解。

在選取初始點X01后,將ΔX1作為變量對問題進(jìn)行優(yōu)化求解,將求得的X1=X01+ΔX1作為下一次運(yùn)算的初始點X02將ΔX2作為變量進(jìn)行迭代運(yùn)算,第k步運(yùn)算表示為

(11)

由第1節(jié)中對雷達(dá)探測情景的分析可知,在導(dǎo)彈能夠被雷達(dá)探測到的情況下,雷達(dá)探測球與導(dǎo)彈軌跡存在2個交點,即優(yōu)化問題存在2個全局最優(yōu)解,在設(shè)置迭代初始值X01時,分別在下方交點的下行位置和上方交點的上行位置取初始點以控制優(yōu)化方向,得到對應(yīng)最優(yōu)解。

3.2 視線平面與彈道軌跡交點計算

在確定探測球與橢圓軌道交點后,需進(jìn)一步求解雷達(dá)視線平面與橢圓軌道的交點得到實際探測范圍,同理,取式(6)中線性函數(shù)的二范數(shù)為目標(biāo)函數(shù),非線性部分為約束條件,將問題轉(zhuǎn)化為優(yōu)化模型:

(12)

令X=X0+ΔX,采用相同方法對約束條件松弛后進(jìn)行迭代凸優(yōu)化求解,第k步運(yùn)算表示為

(13)

在迭代初始值選取上,由于只關(guān)心靠近視平面線與導(dǎo)彈軌跡靠近雷達(dá)一側(cè)的交點,故選擇靠近該點的雷達(dá)坐標(biāo)Xr為優(yōu)化初始值。

以下給出上述采用迭代凸優(yōu)化方法求解非線性方程組算法的正式描述。

1)任務(wù)。

求解非線性方程組。

2)問題轉(zhuǎn)化。

①分離方程組中的線性函數(shù)與非線性函數(shù),取線性函數(shù)的二范數(shù)作為目標(biāo)函數(shù),非線性函數(shù)作為約束條件,將方程組求解問題轉(zhuǎn)化為優(yōu)化問題:

②令X=X0+ΔX,將非線性函數(shù)在X0點處進(jìn)行泰勒展開,取其一階近似:

3)迭代求解。

①初始值選取:設(shè)定迭代次數(shù)k=0,根據(jù)求解目標(biāo)位置,選擇該點附近的已知點作為迭代運(yùn)算的初始點X01;

②迭代次數(shù)k增加1,對凸優(yōu)化問題進(jìn)行求解,得到ΔXk;

③更新解:計算Xk=X0k+ΔXk

④更新初始值:將新解作為新的優(yōu)化初始值X0k+1=Xk;

4 場景模擬與仿真分析

為驗證上述算法,給出3種場景下的導(dǎo)彈軌跡模型與雷達(dá)探測模型,對導(dǎo)彈探測時間進(jìn)行估算,具體參數(shù)如表1所示。表中λ0、φ0、h0分別為導(dǎo)彈發(fā)射點北緯緯度、東經(jīng)經(jīng)度和高度;λf、φf、hf為落點對應(yīng)量;tP為發(fā)動機(jī)脫離時導(dǎo)彈飛行時間;S1、R1max分別為彈頭脫離前的雷達(dá)截面面積和探測距離;λr、φr、hr為雷達(dá)的位置;S2和R2max為彈頭脫離后的雷達(dá)截面積和探測距離。

表1 不同場景下探測的關(guān)鍵參數(shù)

根據(jù)表1參數(shù)建立導(dǎo)彈軌跡模型與雷達(dá)探測模型,設(shè)定迭代終止誤差ε=1×10-10,采用迭代凸優(yōu)化求解算法得到發(fā)動機(jī)脫離前雷達(dá)探測球與導(dǎo)彈軌跡交點X1、X2;發(fā)動機(jī)脫離后雷達(dá)探測球與導(dǎo)彈軌跡交點X3、X4;視線平面與導(dǎo)彈軌跡交點Xc;再通過軌道根數(shù)模型得到發(fā)動機(jī)脫離點位置P,各點在Oxy平面橢圓軌道中的位置關(guān)系分別在圖10～圖12中表示。

以交匯點與導(dǎo)彈彈道的距離Δ=norm(Xk-B)TA-1(Xk-B)-1作為衡量誤差的標(biāo)準(zhǔn),算法求得各交匯點X1、X2、X3、X4、P的計算誤差由表2所示。

表2 不同場景下算法計算誤差

定義導(dǎo)彈發(fā)射時間為零時刻,結(jié)合導(dǎo)彈根數(shù)模型,計算不同場景下雷達(dá)對導(dǎo)彈的探測時間窗口,結(jié)果如表3所示。表中tb為探測起始時間,te為探測終止時間,Δt為探測總時間。

表3 不同場景下雷達(dá)對導(dǎo)彈的探測時間窗口

場景1中的探測時間窗口為導(dǎo)彈發(fā)射后105～120 s,總時間為15 s;場景2中由2個探測時間窗口,分別為導(dǎo)彈發(fā)射后102～120 s與導(dǎo)彈發(fā)射后133～245 s,第一段窗口總時間為18 s,第二段窗口總時間為112 s;場景1中的探測時間窗口為導(dǎo)彈發(fā)射后39～206 s,總時間為167 s。

從表1的計算誤差與表2對雷達(dá)探測時間窗口的計算結(jié)果可以看出,本文提出的方法具有較高的計算精度,同時能夠在不同的雷達(dá)探測情景下實現(xiàn)對探測時間窗口的估計。

5 結(jié)束語

分析仿真計算結(jié)果,本文提出的迭代凸優(yōu)化算法能夠有效求解非線性二次方程組,適用于不同的探測場景且具有較高的精度。有效地解決了遠(yuǎn)程預(yù)警雷達(dá)對彈道導(dǎo)彈探測時間的估算問題。同時該方法對于三維空間目標(biāo)運(yùn)動軌跡交點的求解問題具有較強(qiáng)的通用性,在航天器探測跟蹤、GPS定位等領(lǐng)域同樣具有應(yīng)用前景。

[1] 金林.彈道導(dǎo)彈防御系統(tǒng)綜述[J].現(xiàn)代雷達(dá),2012,34(12):1-7.

JIN Lin.Overview of ballistic missile defense system[J].Modern Radar,2012,34(12):1-7.(in Chinese)

[2] 袁可亮,施自勝,羅金亮,等.彈道導(dǎo)彈突防中的電磁威脅分析及對策[J].艦船電子對抗,2013,36(2):33-34.

YUAN Ke-liang,SHI Zi-sheng,LUO Jin-liang,et al.Analysis and countermeasures of electromagnetism threats in the ballistic missile defense penetration[J].Shipboard Electronic Countermeasure,2013,36(2):33-34.(in Chinese)

[3] LI X R,JILKOV V P.Survey of maneuvering target tracking.Part Ⅱ:motion models of ballistic and space targets[J].IEEE Transactions on Aerospace & Electronic Systems,2010,46(1):96-119.

[4] 曹敏,劉永祥,黎湘.彈道導(dǎo)彈軌跡參數(shù)生成技術(shù)研究[J].電光與控制,2004,11(2):35-37.

CAO Min,LIU Yong-xiang,LI Xiang.Study on generation of trajectory parameters for ballistic missiles[J].Electronics Optics & Control,2004,11(2):35-37.(in Chinese)

[5] 霍偉,王曉鈞.基于彈道仿真的反導(dǎo)系統(tǒng)分析[J].中國電子科學(xué)研究院學(xué)報,2009,4(5):507-510.

HUO Wei,WANG Xiao-jun.Analysis of antimissile system based on trajectory simulation[J].Journal of China Academy of Electronics and Information Technology,2009,4(5):507-510.(in Chinese)

[6] XU H,CHEN W,NING X.Simulation research on missile tracking under the guidance of online real radar[C]//Asian Simulation Conference.Singapore:Springer,2016:519-531.

[7] 朱力,王盛利,于立.彈道導(dǎo)彈的雷達(dá)探測系統(tǒng)仿真[J].計算機(jī)仿真,2000,17(1):46-48.

ZHU Li,WANG Cheng-li,YU Li.The simulation of radar detection system of ballistic missile[J].Computer Simulation,2000,17(1):46-48.(in Chinese)

[8] PEDERSEN R N,MOUNTCASTLE P D,VISS S J.Method and system for analyzing ballistic trajectories:IEEE US9342907[R].2016.

[9] 張海成,楊江平,王晗中.艦載SPY雷達(dá)對彈道導(dǎo)彈探測效能分析[J].現(xiàn)代雷達(dá),2012,34(1):1-4.

ZHANG Hai-cheng,YANG Jiang-ping,WANG Han-zhong.Analysis of detection performance of shipborne SPY radar for ballistic missile[J].Modern Radar,2012,34(1):1-4.(in Chinese)

[10] 谷曉輝,龔春林,谷良賢.導(dǎo)彈運(yùn)動方程求解策略研究[J].計算機(jī)仿真,2013,30(7):92-95.

GU Xiao-hui,GONG Chun-lin,GU Liang-xian.Strategy study on missile motion equations[J].Computer Simulation,2013,30(7):92-95.(in Chinese)

[11] 肖濱,郭鵬程,衡軍,等.戰(zhàn)術(shù)彈道導(dǎo)彈的彈道仿真[J].系統(tǒng)仿真技術(shù),2008,4(4):213-217.

XIAO Bin,GUO Peng-cheng,HENG Jun,et al.Ballistic simulation of tactical ballistic missile[J].System Simulation Technology,2008,4(4):213-217.(in Chinese)

[12] 張光義,趙玉潔.相控陣?yán)走_(dá)技術(shù)[M].北京:電子工業(yè)出版社,2007.

ZHANG Guang-yi,ZHAO Yu-jie.Phased-array radar techn-ology[M].Beijing:Publishing House of Electronics Industry,2007.(in Chinese)

EstimationMethodofMissileDetectionTime-windowBasedonIterativeConvexOptimization

GE Rui-xing,ZHOU Qing-song,ZHANG Jian-yun,CHEN Shi-wa

(School of Electronic Countermeasure,National University of Defense Technology,Hefei 230037,China)

To accurately estimate the time-window of long-range early-warning radar dectecting ballistic missile,the simultaneous mathematical model of missile trajectory and radar detection-range was built,and a method of solving the nonlinear quadratic equations was proposed.The core problem was transformed into an optimization problem,and the error norm was introduced to constrain the non-convex optimal problem.The first and second order Taylor expansion of the quadratic equation was carried out,and the original problem was approximated to an optimization problem locally,and the iterative convex optimization was used to solve the problem.The analytical solution of radar detection-range and trajectory joint-point were obtained to realize the accurate estimation of radar detection time-window.Scene simulation and solution analysis prove that the method is correct,and this method has the characteristics of high precision and versatility.

ballistic missile;detection;radar detection time-window;nonlinear quadratic equations;iterative convex optimization

TJ765.4

A

1004-499X(2017)04-0022-07

2017-05-25

葛瑞星(1993- ),男,碩士研究生,研究方向為導(dǎo)彈實際過程中的電子對抗技術(shù)。E-mail:geruixing1993@sina.com。