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

（1.中國人民解放軍92942 部隊，北京 100161；2.上海機電工程研究所，上海 201109）

0 引言

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

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

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

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

1.1 彈-目相對運動方程

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

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

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

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

取控制向量為

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

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

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

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

對φ(k)線性化，可得

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

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

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

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

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

則離散形式的量測方程為

式中，

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

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

2 自適應(yīng)加速度目標估計方法

2.1 交互多模型算法

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

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

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

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

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

當目標運動模型為Singer 模型時，目標加加速度分量可描述為

式中：α為目標的機動頻率。

目標加速度方差為

式中：amax為目標加速度最大值；pmax為目標機動加速度最大的概率值；p0為目標不機動的概率。

當目標運動模型為“當前”統(tǒng)計模型時，目標加加速度分量可描述為

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

式中：a?max為目標的加速度最小值，通常為先驗常值。

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

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

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

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

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

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

因此，可得k時刻目標加速度方差為

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

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

3 仿真分析

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

式中：ntx、nty、ntz分別為目標機動過載在地面坐標系3 個軸的分量；A為機動過載最大值。

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

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

圖2 過載最大值為10 時加速度估計曲線Fig.2 Estimated curves of the acceleration when the maximum overload is 10

圖3 過載最大值為15 時加速度估計曲線Fig.3 Estimated curves of the acceleration when the maximum overload is 15

圖4 過載最大值為20 時加速度估計曲線Fig.4 Estimated curves of the acceleration when the maximum overload is 20

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

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

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

4 結(jié)束語

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