吳新宏，張紫琪，王 磊，高 帆，仇理寬

（1.中國(guó)人民解放軍92942 部隊(duì)，北京 100161；2.上海機(jī)電工程研究所，上海 201109）

0 引言

現(xiàn)代反艦導(dǎo)彈普遍采用多變彈道技術(shù)［1］結(jié)合超低空飛行進(jìn)行突防，這大大提升了防空反導(dǎo)彈的攔截難度，對(duì)導(dǎo)彈制導(dǎo)律的快速響應(yīng)性和魯棒性提出了更高的要求。現(xiàn)已有諸多成熟的控制理論為制導(dǎo)律的設(shè)計(jì)提供了新的方向，例如，包含目標(biāo)機(jī)動(dòng)補(bǔ)償?shù)谋壤龑?dǎo)引律、滑模變結(jié)構(gòu)制導(dǎo)律、自抗擾控制導(dǎo)引律等。這些制導(dǎo)律在實(shí)際使用中需要獲取目標(biāo)的機(jī)動(dòng)信息，但現(xiàn)有的導(dǎo)引頭并不具備直接測(cè)量目標(biāo)加速度的能力。若不對(duì)導(dǎo)彈加以補(bǔ)償會(huì)造成較大脫靶量［4］，所以需對(duì)目標(biāo)加速度進(jìn)行估計(jì)。在多變彈道技術(shù)中，蛇形機(jī)動(dòng)是一種典型的目標(biāo)機(jī)動(dòng)方式。蛇形機(jī)動(dòng)是指導(dǎo)彈進(jìn)行平面內(nèi)的S 形運(yùn)動(dòng)，其運(yùn)動(dòng)軌跡形狀似蛇。文獻(xiàn)［5］探討了蛇形機(jī)動(dòng)提升導(dǎo)彈機(jī)動(dòng)能力的有效性。文獻(xiàn)［6］經(jīng)過(guò)仿真驗(yàn)證了小周期高頻率的蛇形機(jī)動(dòng)目標(biāo)更有利于提升反艦導(dǎo)彈的突防能力。在現(xiàn)役導(dǎo)彈中，俄羅斯的“白蛉”已具備蛇形機(jī)動(dòng)的能力。

在卡爾曼濾波中，Singer 模型［7］和“當(dāng)前”統(tǒng)計(jì)模型［8］常用來(lái)構(gòu)建機(jī)動(dòng)目標(biāo)的運(yùn)動(dòng)模型，這兩種模型的加速度方差受預(yù)設(shè)加速度極值的影響。目標(biāo)加速度估計(jì)是由目標(biāo)機(jī)動(dòng)頻率和加速度方差共同決定的。當(dāng)目標(biāo)的實(shí)際加速度大于預(yù)設(shè)加速度極值時(shí)，加速度估計(jì)精度將會(huì)下降［9］。另外，由于蛇形機(jī)動(dòng)目標(biāo)具有高機(jī)動(dòng)性和復(fù)雜性，采用單運(yùn)動(dòng)模型難以描述目標(biāo)的運(yùn)動(dòng)，因此，可用交互多模型（Inter?acting Multiple Model，IMM）［11］算法提升模型與真實(shí)目標(biāo)的匹配性。IMM 采用多個(gè)目標(biāo)運(yùn)動(dòng)模型，在線調(diào)整各模型匹配度，最終輸出各模型濾波值的加權(quán)和為其結(jié)果。

針對(duì)以上為問(wèn)題，本文改進(jìn)了卡爾曼濾波算法中目標(biāo)加速度方差的計(jì)算公式，建立了基于擴(kuò)展卡爾曼濾波（EKF）算法的自適應(yīng)交互多模型的目標(biāo)加速度估計(jì)方法，實(shí)現(xiàn)對(duì)不同機(jī)動(dòng)過(guò)載的水平蛇形機(jī)動(dòng)目標(biāo)的加速度估計(jì)。

1 數(shù)學(xué)模型的建立

1.1 彈-目相對(duì)運(yùn)動(dòng)方程

導(dǎo)彈和目標(biāo)的相對(duì)運(yùn)動(dòng)方程，可用彈體執(zhí)行坐標(biāo)系Ox1ay1az1a和彈上視線坐標(biāo)系Oxs1ys1zs1之間的轉(zhuǎn)換關(guān)系來(lái)描述，彈目相對(duì)運(yùn)動(dòng)方程為

式中：R為彈目相對(duì)距離為彈目相對(duì)距離變化率；εs1為彈體視線高低角；βs1為彈體視線方位角為彈體視線高低角速率為彈體視線方位角速率；aTr、aTε、aTβ分 別為目標(biāo)加 速度在坐標(biāo)系Oxs1ys1zs13 個(gè)軸上的投影；aMr、aMε、aMβ分別為導(dǎo)彈加速度在坐標(biāo)系Oxs1ys1zs13 個(gè)軸上的投影。

1.2 狀態(tài)方程和量測(cè)方程

采用EKF，取狀態(tài)向量為

取控制向量為

目標(biāo)加加速度在彈上視線坐標(biāo)系3 個(gè)軸上的投影可分別用函數(shù)far、faε、faβ表示，則狀態(tài)方程可寫(xiě)為

式中：w(t)為均值為零的白噪聲向量。

對(duì)連續(xù)非線性狀態(tài)方程進(jìn)行離散、線性化。假設(shè)系統(tǒng)無(wú)輸入，取采樣時(shí)間為ΔT，將x(t+ΔT)在t進(jìn)行泰勒展開(kāi)，省略二階以上的高階項(xiàng)，可得

令t+ΔT?k+1，t?k，則有

對(duì)φ(k)線性化，可得

將Φ(k)記為狀態(tài)轉(zhuǎn)移矩陣，離散線性化狀態(tài)方程為

式中：w(k)為零均值的過(guò)程白噪聲序列。

w(k)統(tǒng)計(jì)特性如下：

式中：Qk為系統(tǒng)噪聲的誤差方差陣；σkj為克羅尼克函數(shù)。σkj特性如下：

半主動(dòng)式導(dǎo)引頭可提供彈目相對(duì)距離變化率、視線角和視線角速率信息，因此，取觀測(cè)向量為

則離散形式的量測(cè)方程為

式中，

v(k)表示零均值的過(guò)程白噪聲序列，其統(tǒng)計(jì)特性如下：

式中：Rk為噪聲誤差方差陣。

2 自適應(yīng)加速度目標(biāo)估計(jì)方法

2.1 交互多模型算法

蛇形機(jī)動(dòng)目標(biāo)具有機(jī)動(dòng)性強(qiáng)、機(jī)動(dòng)發(fā)生時(shí)刻不確定、機(jī)動(dòng)持續(xù)時(shí)間不確定的特點(diǎn)。選用單一的目標(biāo)運(yùn)動(dòng)形式不能準(zhǔn)確描述蛇形目標(biāo)的運(yùn)動(dòng)狀態(tài)，因此，本文選用多個(gè)運(yùn)動(dòng)模型描述目標(biāo)運(yùn)動(dòng)形式，基于IMM 算法實(shí)現(xiàn)目標(biāo)加速度的估計(jì)。IMM 算法運(yùn)用兩個(gè)及以上的模型來(lái)描述目標(biāo)運(yùn)動(dòng)形式，實(shí)時(shí)調(diào)整各目標(biāo)運(yùn)動(dòng)模型的匹配概率，根據(jù)匹配概率對(duì)濾波結(jié)果進(jìn)行加權(quán)融合，最終輸出機(jī)動(dòng)目標(biāo)的狀態(tài)估計(jì)。

本文選用Singer 模型和“當(dāng)前”統(tǒng)計(jì)模型構(gòu)建IMM 算法的模型集。Singer 模型認(rèn)為目標(biāo)的加速度可以用一階時(shí)間相關(guān)過(guò)程描述，且加速度均值為零，可有效描述目標(biāo)的加速度發(fā)生突變的情況?！爱?dāng)前”統(tǒng)計(jì)模型假設(shè)目標(biāo)的加速度服從均值為“當(dāng)前”時(shí)刻加速度的一階時(shí)間相關(guān)過(guò)程，該模型認(rèn)為下一時(shí)刻目標(biāo)加速度取在“當(dāng)前”時(shí)刻加速度的鄰域內(nèi)，可有效描述目標(biāo)加速度連續(xù)變化的情況。該算法的流程如圖1 所示。

圖1 交互多模型算法流程Fig.1 Flow chart of the IMM algorithm

在該算法中，采用Singer 模型的EKF 濾波器（M1）和采用“當(dāng)前”統(tǒng)計(jì)模型的EKF 濾波器（M2）并行工作。交互多模型算法每次遞推循環(huán)可分為4個(gè)步驟：1）將k?1 時(shí)刻得到的兩個(gè)濾波器的模型概率μj(k?1)和狀態(tài)估計(jì)值協(xié)方差陣作為輸入進(jìn)行計(jì)算，分別得到各模型的混合狀態(tài)估計(jì)值和協(xié)方差 陣作為EKF 的初始條 件；2）根據(jù)k時(shí)刻得到的觀測(cè)值進(jìn)行EKF，得到k時(shí)刻各濾波器的狀態(tài)估計(jì)值(k/k)、協(xié)方差陣(k/k)和觀測(cè)殘差Λj；3）根據(jù)Λj計(jì)算似然函數(shù)，更新濾波器的模型概率μj(k)；4）基于模型概率μj(k)，對(duì)個(gè)濾波器的狀態(tài)估計(jì)值求和得到最終的估計(jì)結(jié)果(k/k)。IMM 算法每次循環(huán)遞推的詳細(xì)步驟見(jiàn)文獻(xiàn)［12］。

2.2 加速度方差自適應(yīng)算法

當(dāng)目標(biāo)運(yùn)動(dòng)模型為Singer 模型時(shí)，目標(biāo)加加速度分量可描述為

式中：α為目標(biāo)的機(jī)動(dòng)頻率。

目標(biāo)加速度方差為

式中：amax為目標(biāo)加速度最大值；pmax為目標(biāo)機(jī)動(dòng)加速度最大的概率值；p0為目標(biāo)不機(jī)動(dòng)的概率。

當(dāng)目標(biāo)運(yùn)動(dòng)模型為“當(dāng)前”統(tǒng)計(jì)模型時(shí)，目標(biāo)加加速度分量可描述為

“當(dāng)前”統(tǒng)計(jì)模型的加速度方差計(jì)算公式如下：

式中：a?max為目標(biāo)的加速度最小值，通常為先驗(yàn)常值。

從式（15）和式（16）中可以看出：Singer 模型和“當(dāng)前”統(tǒng)計(jì)模型的加速度方差的計(jì)算結(jié)果受預(yù)先設(shè)定的加速度最大值的影響。在實(shí)際應(yīng)用中，目標(biāo)的最大加速度通常是未知的。當(dāng)目標(biāo)實(shí)際加速度超出了預(yù)先設(shè)定的加速度最大值時(shí)，加速度估計(jì)精度就會(huì)下降。因此，本文基于預(yù)測(cè)殘差對(duì)加速度方差的計(jì)算方法進(jìn)行了改進(jìn)。

在未獲取量測(cè)值時(shí)，k時(shí)刻的目標(biāo)最優(yōu)速度預(yù)測(cè)值可由k?1 時(shí)刻的速度和加速度表示如下：

在得到量測(cè)值以后，經(jīng)過(guò)殘差修正的速度估計(jì)值為

式中：Δa(k)表示k時(shí)刻的加速度變化量。通過(guò)修正，(k/k)相 較于(k/k?1)考慮了加 速度變 化量Δa(k)的影響，因此Δa(k)還可表示為

式（19）說(shuō)明目標(biāo)加速度變化量和速度預(yù)測(cè)值殘差是線性相關(guān)，文獻(xiàn)［13］研究表明目標(biāo)加速度方差同加速度變化量的絕對(duì)值也是線性相關(guān)的，即

式中：C1為比例系數(shù)，取C1=1。

因此，可得k時(shí)刻目標(biāo)加速度方差為

式中：C為比例系數(shù)。

本文的彈-目相對(duì)運(yùn)動(dòng)方程是建立在彈上視線坐標(biāo)系下的，狀態(tài)量不能直接描述目標(biāo)加速度變化量。根據(jù)彈目相對(duì)速度vs與狀態(tài)量之間關(guān)系，目標(biāo)的加速度方差可表示為

3 仿真分析

本文選用交互多模型算法和方差自適應(yīng)的交互多模型算法，對(duì)水平蛇形機(jī)動(dòng)目標(biāo)的加速度進(jìn)行估計(jì)。假設(shè)導(dǎo)彈保持靜止?fàn)顟B(tài)，則濾波模型的控制向量為0，導(dǎo)彈的初始位置為XM0=[ 0 m 0 m 0 m]，初始速度為VM0=[ 0 m/s 0 m/s 0 m/s]。目標(biāo)的初始位置為XT0=[ 40 000 m 10 m 2 000 m]，初始速度為VT0=[ 800 m/s 0 m/s 0 m/s]，運(yùn)動(dòng)周期為4 s，目標(biāo)機(jī)動(dòng)過(guò)載如下：

式中：ntx、nty、ntz分別為目標(biāo)機(jī)動(dòng)過(guò)載在地面坐標(biāo)系3 個(gè)軸的分量；A為機(jī)動(dòng)過(guò)載最大值。

本文將針對(duì)最大過(guò)載為10、15 和20 的蛇形機(jī)動(dòng)目標(biāo)進(jìn)行仿真。對(duì)彈目視線角速率引入0.005 rad/s 的噪聲，設(shè)定濾波模型中目標(biāo)加速度最大值為150 m2/s。

圖2～圖4 分別為目標(biāo)過(guò)載最大值為10、15 和20 的目標(biāo)加速度估計(jì)曲線，圖5～圖7為經(jīng)過(guò)50 次蒙特卡洛仿真得到的加速度均方根誤差（RMSE）曲線。兩種估計(jì)方法的加速度均方根誤差均值見(jiàn)表1。

圖2 過(guò)載最大值為10 時(shí)加速度估計(jì)曲線Fig.2 Estimated curves of the acceleration when the maximum overload is 10

圖3 過(guò)載最大值為15 時(shí)加速度估計(jì)曲線Fig.3 Estimated curves of the acceleration when the maximum overload is 15

圖4 過(guò)載最大值為20 時(shí)加速度估計(jì)曲線Fig.4 Estimated curves of the acceleration when the maximum overload is 20

圖5 過(guò)載最大值為10 時(shí)加速度均方誤差Fig.5 RMSE curves of the acceleration when the maximum overload is 10

圖6 過(guò)載最大值為15 時(shí)加速度均方根誤差Fig.6 RMSE curves of the acceleration when the maximum overload is 15

圖7 過(guò)載最大值為20 時(shí)加速度均方根誤差Fig.7 RMSE curves of the acceleration when the maximum overload is 20

表1 兩種算法的加速度均方根誤差均值Tab.1 Average values of the acceleration RMSEs obtained from two algorithms m/s2

由圖2～圖4 可知：自適應(yīng)IMM 算法的加速度估計(jì)方法的對(duì)三種不同過(guò)載的目標(biāo)加速度估計(jì)結(jié)果更為良好，響應(yīng)速度更為迅速。由圖5～圖7 可以看出，自適應(yīng)IMM 算法的加速度估計(jì)結(jié)果均方根誤差小于IMM 算法，隨著目標(biāo)機(jī)動(dòng)過(guò)載的增加，IMM算法受限于預(yù)先設(shè)定的目標(biāo)最大加速度值，估計(jì)誤差大大增加。由圖4 可知，當(dāng)目標(biāo)過(guò)載最大值達(dá)到20 的時(shí)候，IMM 算法估計(jì)加速度極大值約為16。表1 顯示了自適應(yīng)IMM 算法的加速度估計(jì)結(jié)果均方根誤差均值小于IMM 算法。而本文通過(guò)改善目標(biāo)加速度方差的計(jì)算方法，實(shí)現(xiàn)了對(duì)目標(biāo)加速度方差的動(dòng)態(tài)調(diào)整，提高了目標(biāo)加速度的估計(jì)精度。

4 結(jié)束語(yǔ)

本文基于交互多模型算法建立了水平蛇形機(jī)動(dòng)目標(biāo)加速度的估計(jì)方法，針對(duì)Singer 模型和“當(dāng)前”統(tǒng)計(jì)模型原有的加速度方差計(jì)算方法受限于預(yù)設(shè)加速度極值的缺陷，建立了新的加速度方差的計(jì)算方法，形成了蛇形機(jī)動(dòng)目標(biāo)加速度估計(jì)的自適應(yīng)IMM 算法。通過(guò)仿真驗(yàn)證，該算法可有效提高目標(biāo)的加速度估計(jì)精度，具有一定的工程應(yīng)用價(jià)值，為遭遇點(diǎn)預(yù)測(cè)［14］和制導(dǎo)律的設(shè)計(jì)奠定了基礎(chǔ)。