宋 豹，于劍橋，楊 迪，郭斐然

(北京理工大學(xué)宇航學(xué)院，北京 100081)

1 引言

攔截概率是武器系統(tǒng)攔截效能的重要評價標(biāo)準(zhǔn)，也是制定多彈協(xié)同攔截策略的重要依據(jù)。在多智能體攻防對抗場景下，只有準(zhǔn)確地得到單攔截器對單目標(biāo)的攔截概率，才能在此基礎(chǔ)上為攔截器合理地選擇攔截對象，充分發(fā)揮攔截器效能，保證攔截方整體的攔截效果達到最優(yōu)[1]。在攻防對抗過程中，攔截策略的制定一般是動態(tài)的，這對攔截概率的計算速度提出了較高的要求[2]。

對于導(dǎo)彈這類非線性時變系統(tǒng)，一般情況下采用蒙特卡羅法(Monte-Carlo)進行精度分析和概率計算[3，4]。蒙特卡洛法是一種基于隨機數(shù)的方法，為了提高仿真精度需要進行大量的模擬計算，時間成本非常高，不適合在戰(zhàn)場上即時運算。

協(xié)方差分析描述函數(shù)技術(shù)(Covariance Analysis Describing Function Technique，CADET)是美國ASC公司提出的一種快速分析非線性制導(dǎo)系統(tǒng)制導(dǎo)精度的方法。自提出以來，該方法得到了廣泛應(yīng)用。Wang X S等[5]利用CADET對初始擾動下的側(cè)向短周期運動進行了仿真分析。蔣瑞民等[6]提出了一種分析存在內(nèi)部參數(shù)攝動的導(dǎo)彈姿態(tài)控制系統(tǒng)的方法。張遠等[7]利用CADET進行落點預(yù)報，提供落點的誤差散布。Zhang T T等[8]基于CADET建立了多種影響因素下的視軸角速度誤差傳播模型。Guo J G等[9]針對一般的高超聲速飛行器姿態(tài)跟蹤系統(tǒng)進行了控制精度分析。以上研究結(jié)果表明，CADET運算速度及準(zhǔn)確度較高，但所分析的物理系統(tǒng)狀態(tài)量較少，結(jié)構(gòu)較為簡單，在非線性強、耦合程度高的復(fù)雜系統(tǒng)上的應(yīng)用效果仍需進一步驗證。

本文首次提出將CADET應(yīng)用于動能攔截器攔截概率的計算這一復(fù)雜非線性問題之中，在保證計算精確度的前提下提高計算效率。首先介紹CADET的基本原理，基于泰勒級數(shù)展開提出高維隨機變量函數(shù)的統(tǒng)計線性化近似求解方案，據(jù)此快速獲得動能攔截器的三維落點均值矩陣和協(xié)方差矩陣；其次跟據(jù)這兩個矩陣建立攔截概率計算公式，利用分部積分法降低了積分維數(shù)；最后將該方法應(yīng)用于動能攔截器六自由度制導(dǎo)控制系統(tǒng)。本文提出的方法在保證計算準(zhǔn)確度的前提下可以大幅縮短攔截概率的計算時間，為多彈協(xié)同攔截策略的制定提供有效依據(jù)。

2 攔截概率快速計算方法

2.1 CADET原理

一般連續(xù)時域非線性系統(tǒng)可以表示為

x(t)=f(x，t)+G(t)w(t)

(1)

其中，f(x，t)為n維狀態(tài)向量x(t)的非線性矢量函數(shù)，x(t)是系統(tǒng)的n維狀態(tài)向量，w(t)為p維隨機輸入向量，G(t)為確定性的函數(shù)矩陣。

隨機狀態(tài)向量x(t)包括確定性分量m(t)和隨機分量r(t)

(2)

其中p(t)為狀態(tài)向量的協(xié)方差矩陣。

同樣，隨機輸入向量w(t)也由均值向量b(t)和u(t)組成

(3)

其中，u(t)是具有譜密度矩陣Q(t)的高斯白噪聲。δ(t-τ)為狄拉克函數(shù)。

對于非線性系統(tǒng)應(yīng)用CADET，在利用描述函數(shù)理論將非線性系統(tǒng)統(tǒng)計線性化后，再結(jié)合協(xié)方差分析的方法進行統(tǒng)計性能分析。

對于非線性矢量函數(shù)f(x，t)，其擬線性化形式可以表示為+Nr。利用描述函數(shù)理論求取合適的N值使得這個近似表達式與f(x，t)的均方差E[eTSe]達到最小，則產(chǎn)生最好近似。其中e=f--Nr，S為任意半正定矩陣。

若假設(shè)x(t)呈聯(lián)合正態(tài)分布，描述函數(shù)和準(zhǔn)線性增益矩陣N可以用下式計算

(4)

其中，g(x)為x(t)的概率密度函數(shù)。

在將非線性系統(tǒng)統(tǒng)計線性化后，便可以得到系統(tǒng)狀態(tài)均值和協(xié)方差矩陣的傳播方程

(5)

式中，只要m(0)和p(0)已知便可通過一次積分得到任意時刻的狀態(tài)均值和協(xié)方差矩陣。

2.2 高維隨機變量函數(shù)的統(tǒng)計線性化近似求解方案

若函數(shù)fi(x)(i=1，2，…，n)關(guān)于狀態(tài)向量x(t)高階可微，將fi(x)在m處進行多元函數(shù)的泰勒展開，可以得到

(6)

其中，s代表泰勒展開階數(shù)。

針對不同的fi(x)適當(dāng)選取泰勒級數(shù)展開階數(shù)并忽略皮亞諾余項可以簡化描述函數(shù)i(x)的計算。一般選取泰勒級數(shù)展開階數(shù)為1～2階即可滿足精度要求。

若選取泰勒級數(shù)展開階數(shù)為2階，則

(7)

(8)

根據(jù)概率論相關(guān)知識，式(8)中

則式(8)可化簡為

(9)

類似的，若選取泰勒級數(shù)展開階數(shù)為1階，則

i(x)=fi(m)

(10)

考慮到六自由度動能攔截器制導(dǎo)控制系統(tǒng)維數(shù)較高，本文選取1階泰勒級數(shù)展開進行近似求解。

矩陣N為

(11)

2.3 降維攔截概率計算公式

利用式(5)可以快速得到彈目相對位置的均值矩陣和協(xié)方差矩陣，但落點散布的均值和協(xié)方差矩陣并不能作為制定多彈協(xié)同策略的直接依據(jù)，為此，需要得到攔截器對對應(yīng)目標(biāo)的攔截概率。

本小節(jié)根據(jù)相對位置的均值矩陣和協(xié)方差矩陣給出攔截概率公式，并將其降維以提高計算效率。

假設(shè)n維隨機變量x服從正態(tài)分布N(m，p)，其概率密度函數(shù)[10]為

(12)

其中，p為隨機向量的協(xié)方差矩陣，m為隨機向量的均值向量。

協(xié)方差矩陣中的非對角元素表征著狀態(tài)變量之間的相關(guān)程度。若p為對角矩陣，則x為獨立的隨機向量，若p為非對角矩陣，則x為非獨立的隨機向量。一般地，對于有控導(dǎo)彈系統(tǒng)，相對位置向量x=[Δx，Δy，Δz]T為非獨立隨機向量。

設(shè)動能攔截器與目標(biāo)的脫靶量在R以內(nèi)認(rèn)為動能攔截器成功攔截目標(biāo)，則攔截概率可以表示為

(13)

其中，Ω為一球型區(qū)域，方程為

(Δx-mΔx)2+(Δy-mΔy)2+(Δz-mΔz)2≤R2

(14)

將概率密度函數(shù)g(x)沿Δz方向積分，則可以將三維概率密度函數(shù)轉(zhuǎn)化為二維概率密度函數(shù)，進一步提高計算效率。

設(shè)

則概率密度函數(shù)g(x)可以降維表示為式(15)。具體推導(dǎo)過程見附錄A。

(Φ(h(mΔz+R))-Φ(h(mΔz-R)))×

(15)

式(15)中Φ(h(Δz))表示累積分布函數(shù)，可通過讀取標(biāo)準(zhǔn)正態(tài)分布表獲得，且

攔截概率相應(yīng)地可以表示為

(16)

其中，Ω′為一圓型區(qū)域，方程為

(Δx-mΔx)2+(Δy-mΔy)2≤R2

(17)

3 快速計算方法在動能攔截器制導(dǎo)控制模型中的應(yīng)用

將本文提出的攔截概率快速計算方法應(yīng)用于動能攔截器制導(dǎo)控制模型中，以驗證該方法在復(fù)雜系統(tǒng)上的可行性。

3.1 動能攔截器制導(dǎo)控制模型

動能攔截器制導(dǎo)控制模型內(nèi)部關(guān)系可由圖1表示。

圖1 動能攔截器制導(dǎo)控制模型

在動能攔截器制導(dǎo)控制模型中，主要隨機狀態(tài)變量如表1所示，共計26項，組成高維隨機狀態(tài)向量x(t)。

表1 動能攔截器制導(dǎo)控制模型中主要隨機變量

3.2 動能攔截器制導(dǎo)控制系統(tǒng)的統(tǒng)計線性化

為保證統(tǒng)計分析的準(zhǔn)確性，需對動能攔截器制導(dǎo)控制系統(tǒng)中的非線性方程均進行統(tǒng)計線性化。為便于計算，本文對于低維函數(shù)方程進行完全統(tǒng)計線性化，采用2.2節(jié)中提出的近似求解方案對高維函數(shù)方程統(tǒng)計線性化。

3.3 動能攔截器的攔截概率

將動能攔截器制導(dǎo)控制系統(tǒng)統(tǒng)計線性化后，運用協(xié)方差分析技術(shù)可以得到系統(tǒng)響應(yīng)的狀態(tài)矢量x(t)的均值和協(xié)方差的方程式(18)。

(18)

其中，m(t)為狀態(tài)向量x(t)的均值

m(t)=[mωy4，mωz4，mamy2，mamz2，mm，mVm，mθm，mψvm，

mXm，mYm，mZm，maTx4，maTy4maTz4，mVTx，mVTy，

mVTz，mxT，myT，mzT，mR，mq1，mq2，mΔq1，mΔq2]T

4 仿真結(jié)果及分析

彈目相對位置的均值矩陣和協(xié)方差矩陣是攔截彈效能評估和命中概率計算的主要依據(jù)。本文分別通過快速計算方法和蒙特卡洛方法對典型攔截場景進行計算機仿真，對比兩種方法得到的彈目相對位置的均值和方差，以分析快速計算方法的準(zhǔn)確性；同時給出兩種方法的仿真時間，以說明快速計算方法的有效性。

4.1 仿真條件設(shè)置

初始值設(shè)置：動能攔截器初始位置(2360km，0km，0km)，初始速度大小6000m/s，初始速度方向指向目標(biāo)，目標(biāo)初始位置(2235.5km，79km，0.5km)，目標(biāo)初始速度(3795m/s，-1881m/s，20m/s)。

擾動量設(shè)置：高低角偏差Δq1均值設(shè)置為0，標(biāo)準(zhǔn)差設(shè)置為5e-5rad；方位角偏差Δq2均值設(shè)置為0，標(biāo)準(zhǔn)差設(shè)置為5e-5rad。

步長設(shè)置：1ms，仿真末端步長減小至0.01ms。

硬件環(huán)境：利用同一臺工作站仿真。

4.2 仿真結(jié)果對比與分析

圖2和圖3分別為目標(biāo)不采取機動策略時彈目相對位置的均值和方差對比曲線。從圖2可以看出，采用快速計算方法得到的相對位置均值和蒙特卡洛方法得到的結(jié)果是一致的，在制導(dǎo)控制系統(tǒng)的作用下逐漸減小。從圖3可以看出采用快速計算方法得到的相對位置方差和蒙特卡洛方法得到的結(jié)果也是一致的，在測角誤差的影響下相對位置的散布逐漸變大。

將落點均值和協(xié)方差代入式(16)中并進行積分運算便可以得到脫靶量特定范圍以內(nèi)的攔截概率。從表2中可以看出，根據(jù)快速計算方法得到的命中概率和進行500次蒙特卡洛仿真得到的結(jié)果非常接近，且從表2中可以看出，采用本文方法對命中概率的計算的快速性較好，大約是500次蒙特卡洛仿真所用時間的1/90，且10秒的運行時間基本可以滿足彈上實時運算的要求。

圖2 均值對比曲線

表2 蒙特卡洛方法和快速計算方法在不同脫靶量限制下的攔截概率

表3 蒙特卡洛方法和快速計算方法用時對比

圖4 t(aT≠0)=2s時攔截概率對比圖

圖5 t(aT≠0)=4s時攔截概率對比圖

圖4和圖5分別為在目標(biāo)機動力持續(xù)作用時間為2s和4s的條件下，根據(jù)快速計算方法得到的攔截概率對比圖。從圖4和圖5可以看出，對于不同的目標(biāo)機動力大小、機動力持續(xù)時間和突防起始位置，攔截概率有很大的區(qū)別。這表明只有在探測到目標(biāo)的機動方式時，才能獲得準(zhǔn)確地攔截概率。而在實際戰(zhàn)場環(huán)境下，對于目標(biāo)的探測時長是有限的，這也進一步說明了提出攔截概率快速計算方法的必要性。

5 結(jié)論

對于動能攔截器攔截概率的統(tǒng)計計算，一般通過蒙特卡洛方法實現(xiàn)，時間成本很高。應(yīng)用本文所提方法對動能攔截器制導(dǎo)控制系統(tǒng)進行處理可以一次性地得到狀態(tài)變量的均值矩陣和協(xié)方差矩陣，將其代入降維攔截概率計算公式即可解算攔截概率，在保證計算準(zhǔn)確度的前提下可以大幅提高計算效率，能夠為戰(zhàn)場實時決策提供有效信息。文中給出的仿真實例說明了所提方法的可行性和有效性。

本文提出的攔截概率計算方法的創(chuàng)新點主要有：為滿足多彈協(xié)同策略的制定對于獲取攔截概率的實時性要求，首次結(jié)合協(xié)方差分析描述函數(shù)法進行攔截概率的計算；基于泰勒級數(shù)展開對傳統(tǒng)CADET進行改進，對復(fù)雜系統(tǒng)的統(tǒng)計線性化采用近似求解方案，能夠快速獲得狀態(tài)向量的均值矩陣和協(xié)方差矩陣；采用分部積分法對三維攔截概率計算公式進行降維，能夠進一步減少計算時間。