一種基于CADET的攔截概率快速計(jì)算方法

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

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

1 引言

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

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

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

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

2 攔截概率快速計(jì)算方法

2.1 CADET原理

一般連續(xù)時(shí)域非線性系統(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維隨機(jī)輸入向量，G(t)為確定性的函數(shù)矩陣。

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

(2)

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

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

(3)

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

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

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

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

(4)

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

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

(5)

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

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

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

(6)

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

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

若選取泰勒級(jí)數(shù)展開(kāi)階數(shù)為2階，則

(7)

(8)

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

則式(8)可化簡(jiǎn)為

(9)

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

i(x)=fi(m)

(10)

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

矩陣N為

(11)

2.3 降維攔截概率計(jì)算公式

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

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

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

(12)

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

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

設(shè)動(dòng)能攔截器與目標(biāo)的脫靶量在R以內(nèi)認(rèn)為動(dòng)能攔截器成功攔截目標(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ù)，進(jìn)一步提高計(jì)算效率。

設(shè)

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

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

(15)

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

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

(16)

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

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

(17)

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

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

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

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

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

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

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

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

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

3.3 動(dòng)能攔截器的攔截概率

將動(dòng)能攔截器制導(dǎo)控制系統(tǒng)統(tǒng)計(jì)線性化后，運(yùn)用協(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é)果及分析

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

4.1 仿真條件設(shè)置

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

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

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

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

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

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

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

圖2 均值對(duì)比曲線

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

表3 蒙特卡洛方法和快速計(jì)算方法用時(shí)對(duì)比

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

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

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

5 結(jié)論

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

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