王向民，王 軍，謝杰濤，郭 治

(1.南京理工大學(xué) 自動(dòng)化學(xué)院，南京 210094;2.中國(guó)人民解放軍 32200 部隊(duì)，遼寧 錦州 121000)

為應(yīng)對(duì)來自空中的威脅，防空反導(dǎo)系統(tǒng)作為最直接有效的反制手段，受到了各國(guó)的高度重視.以高炮為主的近程末端防空武器系統(tǒng)，其具有較高的費(fèi)效比，常作為防御體系的最后屏障，因而它對(duì)來襲目標(biāo)的毀傷關(guān)系到整個(gè)防空作戰(zhàn)的成敗.為獲取高炮武器系統(tǒng)的毀傷概率：一種方法是通過大量的實(shí)彈射擊，統(tǒng)計(jì)出射彈數(shù)目和實(shí)際毀傷目標(biāo)的數(shù)量，二者比值即為武器的毀傷概率,該方法可信度高，但檢驗(yàn)費(fèi)用昂貴，且不能適用于武器系統(tǒng)的設(shè)計(jì)與論證階段；另一種方法是建立基于高炮射擊誤差的毀傷概率數(shù)學(xué)模型，利用數(shù)值計(jì)算的方法獲得高炮武器系統(tǒng)的毀傷概率，該方法通過對(duì)高炮的射擊誤差進(jìn)行適當(dāng)分解，利用誤差模型轉(zhuǎn)換法構(gòu)建毀傷概率的計(jì)算模型.在我國(guó)現(xiàn)行的國(guó)家軍用標(biāo)準(zhǔn)(簡(jiǎn)稱國(guó)軍標(biāo))中，把誤差轉(zhuǎn)換為重復(fù)誤差和不可重復(fù)誤差，雖然降低了計(jì)算復(fù)雜度，但存在較大的模型近似[1-4].射擊誤差中的弱相關(guān)誤差源對(duì)于毀傷概率計(jì)算具有顯著影響，陶德敬等[5]提出將射擊誤差分解為共有分量和非共有分量，運(yùn)用誤差序列預(yù)測(cè)理論給出了確定強(qiáng)相關(guān)誤差和弱相關(guān)誤差下射擊過程的平均毀傷概率表達(dá)式，但該模型只能滿足誤差序列為2階狀態(tài)方程以下的情況；姚志軍等[6]利用弱相關(guān)射擊誤差狀態(tài)方程給出了遞推形式的毀傷概率計(jì)算方法，但同樣對(duì)于高階自相關(guān)誤差序列的相關(guān)性計(jì)算沒有明確給出解析式.隨著現(xiàn)代高炮武器系統(tǒng)射頻的提高，射擊過程中前一發(fā)彈藥的毀傷概率對(duì)其后各發(fā)都有影響，因而每次發(fā)射的彈藥對(duì)毀傷概率的貢獻(xiàn)都是有差別的，而將誤差序列簡(jiǎn)單地近似為1階或2階系統(tǒng)的模型簡(jiǎn)化方法[7]，有悖于高炮射擊的真實(shí)模型.高炮的誤差序列從本質(zhì)上來說，可以視為武器身管的隨機(jī)振動(dòng)、跟蹤誤差以及火控解算誤差等綜合而成的平穩(wěn)時(shí)間序列，有關(guān)時(shí)間序列相關(guān)性的分析方法，包括線性和非線性回歸方法進(jìn)行預(yù)測(cè)[8]、基于數(shù)據(jù)驅(qū)動(dòng)的時(shí)間序列預(yù)測(cè)等[9-10]，需要有足夠的實(shí)驗(yàn)數(shù)據(jù)才能獲得較高的預(yù)測(cè)精度[11-12]，并且隨著時(shí)間序列相關(guān)性階數(shù)的增加，建模的復(fù)雜性和計(jì)算量也會(huì)大幅增加[13].

為了建立更精確的毀傷概率計(jì)算模型，本文以射擊誤差序列所構(gòu)建的狀態(tài)方程為出發(fā)點(diǎn)，運(yùn)用卡爾曼濾波的思想，以遞推方式給出連續(xù)脫靶條件下射擊諸元誤差的密度函數(shù)，將包括由射擊沖擊載荷影響在內(nèi)的射擊誤差的弱相關(guān)分量分解為完全可預(yù)測(cè)與完全不可預(yù)測(cè)的分量，再分別與完全可預(yù)測(cè)的強(qiáng)相關(guān)分量和完全不可預(yù)測(cè)的不相關(guān)分量合并，構(gòu)成2個(gè)隨機(jī)分量，根據(jù)高炮進(jìn)行點(diǎn)射時(shí)的時(shí)空特性，構(gòu)建射擊誤差的隨機(jī)狀態(tài)模型，進(jìn)而推導(dǎo)出相應(yīng)的毀傷概率遞推計(jì)算模型.

1 射擊誤差序列狀態(tài)方程構(gòu)建

x(k)=xb(k)+xr(k)+xq(k)

(1)

式中：xr(k)與xq(k)之和稱為射擊諸元誤差xs(k).射擊誤差序列的狀態(tài)方程可進(jìn)一步表示為

(2)

作為隨機(jī)狀態(tài)方程，還要已知它的初始狀態(tài)方差才能是完備的.設(shè)它的初始狀態(tài)方差為

(3)

考慮到射擊誤差為平穩(wěn)高斯過程，因而可得[5]

(4)

式中：Rn為n階預(yù)測(cè)系數(shù)；w(k)為白噪聲.此時(shí)隨機(jī)誤差序列可分解成可預(yù)測(cè)分量與不可預(yù)測(cè)分量.

2 射擊誤差序列狀態(tài)方程構(gòu)建

射擊誤差序列的一步預(yù)測(cè)特性是影響毀傷概率的一個(gè)決定性因素.對(duì)一個(gè)正常的射擊過程而言，其射擊誤差是由兩部分構(gòu)成：射擊諸元誤差，它是經(jīng)歷了足夠長(zhǎng)的時(shí)間后，暫態(tài)分量可以忽略了的相關(guān)平穩(wěn)序列；射彈散布誤差，它雖然也是平穩(wěn)序列，但是它是在擊發(fā)后才加入到射擊誤差序列之中，并且在有限次數(shù)的發(fā)射后，即行終止的平穩(wěn)序列.

2.1 射擊誤差序列的統(tǒng)計(jì)特性

對(duì)射擊誤差而言，它的不相關(guān)、弱相關(guān)和強(qiáng)相關(guān)3個(gè)平穩(wěn)隨機(jī)序列中的每一個(gè)都是由更多的獨(dú)立的子序列之和構(gòu)成.對(duì)誤差序列的一步預(yù)測(cè)方程,可表示為

(5)

(6)

并且有

(7)

(8)

兩個(gè)獨(dú)立序列和的均值函數(shù)是其k步預(yù)測(cè)值，而協(xié)方差函數(shù)是其k步預(yù)測(cè)方差.兩個(gè)隨機(jī)序列之和存在一個(gè)過渡過程.當(dāng)0<|r|<1且k→∞時(shí)，此過渡過程才能結(jié)束，兩個(gè)獨(dú)立序列的和的均值函數(shù)與方差函數(shù)才能如同隨機(jī)變量一樣求和與方差.顯然，只有兩種情況是例外，即兩個(gè)獨(dú)立序列均為強(qiáng)相關(guān)序列或均為不相關(guān)序列，這是因?yàn)樵谶@兩種情況下，其任何一步的預(yù)測(cè)均值與方差都是不變的，且與其初值相同.

如果射擊誤差序列中的弱相關(guān)分量的過渡過程可以在射擊準(zhǔn)備時(shí)間內(nèi)衰減到可以忽略的條件，射擊過程中的所有弱相關(guān)分量與強(qiáng)相關(guān)分量才可以依隨機(jī)變量求和規(guī)則予以處理.例如：射擊諸元誤差xs(k)有3個(gè)獨(dú)立序列構(gòu)成，即

xs(k)=xg(k)+xw(k)+xq(k)

(9)

(10)

且有xs(k)的一步預(yù)測(cè)系數(shù)表達(dá)式為

(11)

2.2 連續(xù)脫靶條件下射擊諸元誤差的密度函數(shù)

f(xs(M),xs(M-1),…,xs(n),…,xs(1))=

(12)

Xs(k-1)=

[xs(k-1)xs(k-2) …xs(k-n)]T

(13)

若射擊諸元誤差的最高階次為n，根據(jù)式(12)給出的處理方法，xs(k)在當(dāng)k>n時(shí)，則預(yù)測(cè)樣本空間為

{xs(k)|X(k-1)}∈

(14)

此時(shí)xs(k)為一種雙重正態(tài)分布，它作為射擊誤差的動(dòng)態(tài)均值，與前面的各個(gè)動(dòng)態(tài)均值都有關(guān).

由于k∈[1,n]期間，高炮武器控制系統(tǒng)完成了振動(dòng)的過渡過程，當(dāng)k>n，xs(k)進(jìn)入平穩(wěn)序列，表明高炮隨時(shí)可以實(shí)施射擊.

對(duì)首次擊發(fā)的第一枚彈藥，由于射彈散布誤差xb(1)使x(1)=xb(1)+xs(1)的方差產(chǎn)生突變，此時(shí)，x(1)的條件密度函數(shù)可表述為

f(x(1))=f(x(1)|xs(1)f(xs(1)))

(15)

因而在首發(fā)彈藥脫靶條件下，射擊諸元誤差的密度函數(shù)為

g1(xs(1))=

(16)

式中:s1為當(dāng)k=1時(shí)刻靶標(biāo)的迎彈面，其補(bǔ)集為R1.因此首發(fā)彈藥脫靶的概率為

(17)

根據(jù)條件概率密度的定義，連續(xù)兩次發(fā)射均脫靶條件下，射擊諸元誤差的密度函數(shù)為

g1(xs(1))dxs(1)

(18)

而連續(xù)兩發(fā)均脫靶的概率為

P(x(1)?s1,x(2)?s2)=

(19)

依據(jù)上述推導(dǎo)方式，可以得到前N發(fā)彈藥均脫靶條件下，射擊諸元誤差的密度函數(shù)為

gN-1(xs(N-1))dxs(N-1)=

gN-1[xs(N-1)]dxs(N-1)

(20)

上述密度函數(shù)的積分即為連續(xù)發(fā)射N發(fā)彈藥均脫靶的概率.對(duì)于滿足0-1毀傷律的高炮武器系統(tǒng)，在連續(xù)發(fā)射N發(fā)彈藥的條件下，至少命中一發(fā)的概率可表示為

(21)

如果將毀傷定義為命中，則式(21)為射擊諸元誤差存在相關(guān)性時(shí)，高炮武器進(jìn)行N發(fā)連續(xù)射擊的毀傷概率.

3 單管高炮武器系統(tǒng)的毀傷概率

單管高炮武器系統(tǒng)是指在同一瞬時(shí)僅能發(fā)射一發(fā)彈藥.對(duì)于轉(zhuǎn)管火炮，雖然在外形上具有多個(gè)發(fā)射管，但每個(gè)火炮身管只能在轉(zhuǎn)到特定位置后，依單管規(guī)則順序發(fā)射，因此在計(jì)算毀傷概率時(shí)，仍然可按照單管方式加以計(jì)算.

(22)

再考慮到射彈散布誤差序列是在開始時(shí)與射擊準(zhǔn)備誤差序列求和的，故必須以已知射擊諸元Zs(k)條件下的分布特性來表述.基于上述分析可知第一發(fā)彈藥射擊諸元誤差的密度函數(shù)為

(23)

從式(23)出發(fā)，可遞推出前k發(fā)均不毀傷條件下，射擊諸元誤差的密度函數(shù)為

gk(Zs(k))=

gk-1(Zs(k-1))dZs(k-1)

(24)

則單管高炮武器一次N發(fā)的點(diǎn)射的毀傷概率為

(25)

式中：Σb和Σs為協(xié)方差矩陣；det為行列式；I為單位矩陣.

4 數(shù)值仿真與分析

當(dāng)射擊誤差中的弱相關(guān)分量為1階平穩(wěn)序列模型時(shí)，分別采用國(guó)軍標(biāo)、Monte Carlo(MC)模擬、文獻(xiàn)[6]遞推法以及本文方法進(jìn)行計(jì)算，仿真計(jì)算結(jié)果如圖1所示.當(dāng)射擊誤差中的弱相關(guān)分量采用3階平穩(wěn)序列模型時(shí)，誤差序列的各階自相關(guān)系數(shù)為：ri,i=0.81,ri,i+1=0.72,ri,i+2=0.66 (i=1,2,…,21)，仿真計(jì)算結(jié)果如表1所示.當(dāng)射擊誤差中的弱相關(guān)分量采用4階平穩(wěn)序列模型時(shí)，誤差序列的各階自相關(guān)系數(shù)為：ri,i=0.81,ri,i+1=0.72,ri,i+2=0.66，ri,i+3=0.60 (i=1,2,…,20)，仿真計(jì)算結(jié)果如表2所示.

圖1 相關(guān)系數(shù)與點(diǎn)射毀傷概率關(guān)系對(duì)比圖Fig.1 Comparison of bust firing damage probability with correlation coefficient

表1 具有3階自相關(guān)性的射擊誤差序列毀傷概率Tab.1 Damage probabilities of firing error sequence with third-order autocorrelation

表2 具有4階自相關(guān)性的射擊誤差序列毀傷概率Tab.2 Damage probabilities of firing error sequence with fourth-order autocorrelation

上述仿真結(jié)果表明，當(dāng)射擊誤差序列考慮為1階平穩(wěn)序列模型時(shí)，4種計(jì)算方法的結(jié)果毀傷概率結(jié)果差異性很小.但隨著誤差相關(guān)性階數(shù)的增加，國(guó)軍標(biāo)的計(jì)算結(jié)果與其他3種存在較為明顯的差距，這表明將高階自相關(guān)性誤差序列模型近似為1階模型，會(huì)產(chǎn)生一定的模型計(jì)算誤差.而運(yùn)用本文給出的預(yù)測(cè)系數(shù)計(jì)算方法可以較好地適應(yīng)各類高階模型，計(jì)算結(jié)果更加符合射擊誤差的實(shí)際情況.由于仿真中的各弱相關(guān)誤差序列的隨機(jī)特性已經(jīng)給定，采用MC方法模擬射擊過程中的各弱相關(guān)誤差源，可以忽略模型誤差，通過大量的模擬抽樣所得到的數(shù)值計(jì)算結(jié)果是可信的.因此可以將MC方法的計(jì)算結(jié)果作為比較基準(zhǔn)進(jìn)行準(zhǔn)確性驗(yàn)證，通過表1和表2計(jì)算結(jié)果對(duì)比可看出，國(guó)軍標(biāo)方法有20%以上的偏差，文獻(xiàn)[6]方法的偏差為10%以上，而本文計(jì)算結(jié)果的偏差小于3%.

對(duì)于計(jì)算過于復(fù)雜而難以求得解析解的隨機(jī)過程問題，雖然MC模擬可以通過構(gòu)造符合一定規(guī)則的隨機(jī)數(shù)來進(jìn)行數(shù)值求解，但正如前文所述，考慮到高炮武器系統(tǒng)驗(yàn)證的費(fèi)效比問題，高炮武器系統(tǒng)從設(shè)計(jì)到部隊(duì)列裝，不可能通過大量的試驗(yàn)進(jìn)行毀傷概率的驗(yàn)證，在試驗(yàn)數(shù)據(jù)較少的情況下，使用MC模擬射擊打靶過程，存在模型誤差的問題，有可能產(chǎn)生較大的計(jì)算誤差.為進(jìn)一步說明此問題，利用GPML Toolbox version 4.1工具箱在MATLAB R2013a下，生成一個(gè)均值為0，方差為1的高斯隨機(jī)過程，并產(chǎn)生 10 000 組樣本，每個(gè)樣本均勻采樣20個(gè)點(diǎn)，統(tǒng)計(jì)采樣點(diǎn)落在某一區(qū)間的占比，用于模擬高炮打靶發(fā)射彈藥數(shù)與命中目標(biāo)數(shù)的比例，即為命中概率.根據(jù)MC模擬的原則，分別抽取不同數(shù)量的樣本m1，每組總采樣點(diǎn)個(gè)數(shù)為20m1，記錄采樣點(diǎn)落在[-0.5，0.5]之間的個(gè)數(shù)為n1，記p1=n1/(20m1)，稱為抽樣落點(diǎn)占比.而總的采樣點(diǎn)個(gè)數(shù)為 200 000 個(gè)，統(tǒng)計(jì)得到全體采樣點(diǎn)落在[-0.5，0.5]之間的個(gè)數(shù)為 93 662 個(gè)，記p2=93 662/200 000=0.468，稱為全體落點(diǎn)占比.不同采樣數(shù)量下的落點(diǎn)占比與全體落點(diǎn)占比的統(tǒng)計(jì)結(jié)果如圖2所示.

圖2 MC算法抽樣數(shù)量與落點(diǎn)占比的關(guān)系Fig.2 Number of samples versus proportion of dropping points in MC algorithm

從圖2中可知，當(dāng)抽樣數(shù)量較少時(shí)，MC模擬的結(jié)果具有較大的不確定性.模擬過程中還發(fā)現(xiàn)，即使抽樣數(shù)量相同的不同樣本，MC模擬的數(shù)值統(tǒng)計(jì)結(jié)果也存在明顯差異，這對(duì)于試驗(yàn)數(shù)據(jù)較少，且要求毀傷概率計(jì)算具有嚴(yán)格數(shù)學(xué)推理的武器研制與驗(yàn)證單位是不能接受的.而本文的遞推估計(jì)模型具有嚴(yán)格的數(shù)學(xué)推導(dǎo)，只需確定弱相關(guān)誤差的均值、方差和相關(guān)系數(shù)，就可以準(zhǔn)確計(jì)算出毀傷概率.

5 結(jié)語

提出了一種利用高炮武器射擊誤差的狀態(tài)方程估計(jì)毀傷概率的方法，在將射擊誤差序列中的弱相關(guān)部分，分解為可預(yù)測(cè)誤差和不可預(yù)測(cè)誤差的，給出了可預(yù)測(cè)系數(shù)的計(jì)算方法，并利用遞推方法建立了高炮武器連續(xù)點(diǎn)射的毀傷概率計(jì)算模型，可以較好地適用于射擊誤差序列具有多階自相關(guān)性時(shí)高炮武器系統(tǒng)的毀傷概率計(jì)算.本方法也可以作為各類速射火炮進(jìn)行連續(xù)射擊時(shí)毀傷概率計(jì)算的一種參考.