中国科技核心期刊      中国指挥与控制学会会刊     军事装备类重点期刊
Theory & Research

Adaptive Interacting Multi-model Algorithm Based on Target Moving Angular Velocity

  • LIU Yong-pu 1 ,
  • WANG Pei-cheng 2
Expand
  • 1. Jiangsu Automation Research Institute, Lianyungang 222061
  • 2. Chang’an University, Xi’an 710061, China

Received date: 2020-05-06

  Request revised date: 2020-06-12

  Online published: 2022-05-07

Copyright

Copyright reserved © 2022

Abstract

In view of the mobile target tracking problem, an adaptive multi-model filtering algorithm based on the target line coordinate system is proposed while the tracking sensor can provide the target moving angular velocity. Based on the classical interactive multi-model algorithm, the algorithm introduces the target speed measurement information, and uses the azimuth method to transform the parameters between the target line coordinate systems at different moments, and makes full use of the target speed measurement information and the noise measurement error information to track the maneuvering target. The simulation results show that it can effectively shorten the filter convergence time and improve the filtering accuracy of the mobile target comparing with the single-model Kalman filtering algorithm.

Cite this article

LIU Yong-pu , WANG Pei-cheng . Adaptive Interacting Multi-model Algorithm Based on Target Moving Angular Velocity[J]. Command Control and Simulation, 2020 , 42(5) : 25 -31 . DOI: 10.3969/j.issn.1673-3819.2020.05.005

随着科学技术的发展,反舰导弹升级换代的速度明显加快。各国为获取前所未有的打击优势,研制的新一代反舰导弹都朝着超声速、高机动的方向发展,如俄罗斯的3M-5E俱乐部、印度/俄罗斯的布拉莫斯、中国台湾的雄风III等超声速导弹等[1,2]。同时,新型雷达、光电等多种传感器在机动目标跟踪目标领域得以应用,新的更多的量测数据如目标径向速度、角速率等信息得以获取。美国的“密集阵”系统是利用角速率信息实现坐标系转换、角速率跟踪方面的典范,它极大地提高了跟踪性能,减少了反应时间,在实战和演练中都取得了不俗的成绩。如何快速准确地估计机动目标[3]状态信息,提高跟踪精度,缩短反应时间,已成为防空反导作战中一个重要的研究课题。
传统线性最小方差意义下的α-β滤波、Kalman滤波等线性滤波算法由于目标运动模型单一,无法解决算法模型与机动目标真实运动模型的匹配问题,不能较好地对目标运动状态进行估计。对此,国内外学者提出了交互式多模型滤波算法[4,5,6,7],通过多个模型的交互作用实现对机动运动状态的估计,可较好地克服单一滤波模型本身的误差,被认为是解决机动目标跟踪问题的主流方法。
本文在经典交互多模型的基础上,提出了一种基于目标运动角速度的自适应交互多模型滤波算法。该算法在瞄准线坐标系下对目标运动进行建模,充分利用目标速度量测和噪声量测误差信息进行机动目标跟踪,克服了直角坐标系中噪声的相关性,增强了算法运动模型对目标实际机动运动的匹配性,通过仿真分析证明了算法的正确性。

1 目标运动角速度的引入

理论计算与实践证明,在传统目标量测信息的基础上,目标速度信息的引入可有效提高对目标的滤波精度。Fitzgerald[8]采用三状态模型研究了上述问题,结果发现速度量测的引入对位置估计精度提高可达一个数量级以上。周宏仁[9]从理论和仿真方面通过增加观测矩阵的秩对改善跟踪精度问题进行了相关研究。张怀根[10]将径向速度引入目标跟踪,提高了跟踪精度,缩短了收敛时间。王勇[11]将目标运动角速度引入火控滤波,定量分析了目标运动角速度对滤波性能的影响。下面以传递函数思想,在连续时间域内,对速度量测的引入做一定理论分析。
假定目标做一维匀速运动,设目标状态方程为

X ˙(t)=ΦX(t)+w(t)

式中X(t)= x ( t ) x ˙ ( t )为目标状态向量,系统状态矩阵Φ= 0 1 0 0,w(t)是均值为零的高斯白噪声。
目标观测方程为

Y(t)=HX(t)+υ(t)

式中H= 1 0,υ(t)是均值为零,方差为R的高斯观测噪声。
根据卡尔曼滤波,状态估计方程为
$\begin{array}{c}\hat{\dot{X}}(t)=\boldsymbol{\Phi} \hat{\boldsymbol{X}}(t)+\boldsymbol{K}(t)[\boldsymbol{Y}(t)-\boldsymbol{H} \hat{\boldsymbol{X}}(t)]=[\boldsymbol{\Phi}- \\\boldsymbol{K}(t) \boldsymbol{H}] \hat{\boldsymbol{X}}(t)+\boldsymbol{K}(t) \boldsymbol{Y}(t)\end{array}$
式中K(t)为卡尔曼增益矩阵。
假定卡尔曼滤波器达到稳态,K(t)趋于常数K,对式(3)两边进行拉普拉斯变换,并设初始状态为零,则有

X ^(s)=(sI-Φ+KH)-1KY(s)

式中,s为拉普拉斯算子,X(s)=L{X(t)},Y(s)=L{Y(t)}。
对式(2)两边进行拉普拉斯变换,并设初始状态为零,则有

Y(s)=HX(s)+υ(s)

式中,υ(s)=L{υ(t)}。
由式(4)、(5)可得

X ^(s)= X ^ d(s)+ X ^ r(s)

X ^ d(s)=(sI-Φ+KH)-1KHX(s)

X ^ r(s)=(sI-Φ+KH)-1KHυ(s)

式中, X ^ d(s)相当于滤波器稳态部分, X ^ r(s)相当于滤波器输出暂态部分。由卡尔曼滤波器特性可知,当滤波器达到稳定状态时, X ^ r(s)输出为零。
K= k 1 k 2 T,将ΦH代入式(7)并整理得
X ^ d (s)= 1 s ( s + k 1 ) + k 2 k 1 s + k 2 0 k 2 s 0
由上式,易推知

x ^ d x ^ · d= 1 s ( s + k 1 ) + k 2 k 1 s + k 2 0 k 2 s 0 x ˙ s x ^ · s=

k 1 s + k 2 s ( s + k 1 ) + k 2 x ( s ) k 2 s ( s + k 1 ) + k 2 x . ( s )
当采用目标速度量测数据时,观测矩阵的秩由1变为2,观测矩阵变为

H= 1 0 0 0

此时,卡尔曼增益矩阵为
K= k 11 k 12 k 21 k 22
将式(11)、(12)代入式(7)得
$\left[\begin{array}{c}\hat{x}_{d} \\\hat{\dot{x}}_{d}\end{array}\right]=\frac{1}{\left(s+k_{11}\right)\left(s+k_{12}\right)-k_{21}\left(k_{12}-1\right)}\left[\begin{array}{cc}k_{11}\left(s+k_{22}\right)-k_{21}\left(k_{12}-1\right) & k_{12}\left(s+k_{22}\right)-k_{22}\left(k_{12}-1\right) \\-k_{11} k_{21}+k_{21}\left(s+k_{11}\right) & -k_{12} k_{21}+k_{22}\left(s+k_{11}\right)\end{array}\right]\left[\begin{array}{l}x(s) \\\dot{x}(s)\end{array}\right]$
进一步整理得
$\left[\begin{array}{c}\hat{x}_{d} \\\hat{\dot{x}}_{d}\end{array}\right]=\left[\begin{array}{c}\frac{k_{12} s\left(s+k_{22}\right)-k_{22} s\left(k_{12}-1\right)+k_{11}\left(s+k_{22}\right)-k_{21}\left(k_{12}-1\right)}{\left(s+k_{11}\right)\left(s+k_{12}\right)-k_{21}\left(k_{12}-1\right)} x(s) \\\frac{k_{22}\left(s+k_{11}\right)-k_{21}\left(k_{12}-1\right)}{\left(s+k_{11}\right)\left(s+k_{12}\right)-k_{21}\left(k_{12}-1\right)} \dot{x}(s)\end{array}\right]$
比较式(10)、(13)可发现,两式传递函数的分母多项式中,s的幂都为2,但方程(13)的分子多项式增加了一次幂,这意味着目标速度量测数据的引入相当于在关于位置、速度的传递函数中增加了一个零点,滤波器的带宽将会增大,动态误差将减小,从而减小滤波的收敛时间,提高量测信息的滤波精度。
类似匀速直线运动,对匀加速等其他运动模型而言,引入速度量测同样能够提高滤波器的跟踪精度,在此不再详细论述。

2 目标线速度计算

瞄准线坐标系OXaYaZa,原点O是跟踪旋转轴与俯仰轴的交点,OXa轴为跟踪器俯仰轴,OYa为跟踪器瞄准线,OZa垂直于OXaYa的方向,轴OXa、轴OYaOZa组成右手直角坐标系。如图1q、Δε为方向角和高低角跟踪偏差。设ω为瞄准线坐标系旋转角速度。
图1 瞄准线坐标系与跟踪误差角
假定目标距离为D,目标的直角坐标分别为xayaza,则在瞄准线坐标系下,目标位置ra

ra= x a y a z a= D cosΔ ε sinΔ q D cosΔ ε cosΔ q D sinΔ ε

根据Coriolios定理可知,目标在瞄准线坐标系下运动线速度为
v a =( D ˙cosΔεsinΔq-DΔ ε ˙sinΔεsinΔq+DΔ q ˙cosΔεcosΔq i a+( D ˙cosΔεcosΔq-DΔ ε ˙sinΔεcosΔq-DΔ q ˙cosΔεsinΔq j a+( D ˙sinΔε+DΔ ε ˙cosΔε j a+ ω × r a
其中,ω=(ωxa,ωya,ωza)为瞄准线坐标系旋转角速度,

ω × r a=ra[(sinΔε·ωya-cosΔεcosΔq·ωza i a+(-sinΔε·ωxa+cosΔεsinΔq·ωza j a+(cosΔεcosΔq·ωxa-cosΔεsinΔq·ωya k a]

在传感器角伺服系统具有理想的稳定跟踪状态时,即跟踪偏差Δqε=0时,将式(14)、(16)代入式(15)可得目标线速度为
v xa v ya v za = - D ω za D ˙ D ω xa
其中, D ˙为目标的径向速度,目标距离的微分。
在不考虑舰艇摇摆的情况下,从瞄准线到甲板坐标系的运动角速率可表示成如下形式
ω = ε ˙ + q ˙

Mq= cos q sin q 0 sin q cos q 0 0 0 1,Mε= 1 0 0 0 cos ε sin ε 0 - sin ε cos ε

是坐标系转换的方向余弦矩阵,其中,q为甲板坐标系中跟踪瞄准线的舷角,右舷为正;ε为甲板坐标系中跟踪器瞄准线的俯仰角,向上为正;由公式(18)可得瞄准线坐标系相对于稳定平台角速率为

ωa= ω xa ω ya ω za=MεMq 0 0 - q ˙+Mε ε ˙ 0 0

简化可得

ωa= ω xa ω ya ω za= ε ˙ - q ˙ sin ε - q ˙ cos ε

式(20)给出了瞄准线坐标系相对于稳定平台的角速率公式,可通过安装在瞄准线坐标系的速率陀螺装置测出其相对于惯性空间中的角速率值。
在舰艇摇摆情形下,瞄准线坐标系要保持对目标的稳定跟踪,需平台稳定装置抵消舰艇摇摆影响,此时安装在瞄准线坐标系中的速率陀螺装置测出的瞄准线坐标系运动角速率包含摇摆引起的滚转角速率和目标运动引起的角速率两部分。为研究方便起见,本文考虑舰艇无摇摆情形。

3 引入目标运动速度的自适应交互多模型算法

基于角速率的自适应交互多模型算法是经典多模型算法的进一步延伸。经典交互多模型算法在稳定的直角坐标系中进行滤波,不需要对量测和目标状态等进行转换,而本文提出的算法是在自适应坐标系中进行滤波。自适应坐标系是一种假想坐标系,其沿雷达天线每一坐标轴(瞄准线坐标系)建立目标加速度模型,假设天线框是瞬时固定的(即在量测间隔时间内不相对惯性空间旋转),该坐标系并不时刻与瞄准线坐标系重合,只有在新的量测输入(跟踪雷达和速率陀螺)情况下,自适应坐标系才经过变换与瞄准线坐标系完全重合,其他时刻自适应坐标系在空间是固定不变的。用交互多模型算法处理目标的位置、速度的量测,就像在最后获得目标量测的瞬间,瞄准线坐标系被固定一样。目标状态估计器就是相对于自适应坐标系来计算目标的位置、速度和加速度的。在交互多模型的每个滤波过程结束后,应将目标现在的状态估计与误差协方差等转换到下一个相应的新指向。
对舰艇近程防御而言,目标近程运动往往以直线运动为主。模型集选择以CV模型、CA模型和CS模型等直线运动模型为主比较合理。在此以CV模型、CA模型组合对基于角速率的自适应交互多模型算法进行介绍,模型1为CV模型,模型2为CA模型。
假设滤波器方程在三个空间坐标轴上完全独立,状态空间和量测空间通过测量方程Yk=HkXk+Vk联系起来。滤波在自适应坐标系的三个轴上进行,以x轴方向为例,在步进坐标系中对目标运动模型进行建模,x为位置, x ˙为速度, x ¨为加速度, x 为加加速度,目标的状态向量包括距离量测和速度量测两项,建立目标运动模型。
模型1目标运动方程为

x ˙ x ¨= 0 1 0 0 x x ˙+ 0 1w(t)

模型2目标运动方程为

x ˙ x ¨ x = 0 1 0 0 0 1 0 0 0 x x ˙ x ¨+ 0 0 1w(t)

量测方程为

x x ˙= 1 0 0 0 1 0 x x ˙ x ¨+ 0 0 1υ(t)

其中,w(t)是均值为零,方差为σ2的高斯白噪声。υ(t)是测量误差,是均值为零的高斯白噪声。
同理,建立yz方向上的目标运动模型。
在经典多模型算法基础上,本文提出算法的一个递推循环主要由输入交互、滤波、概率更新、输出交互和坐标转换五个部分。
1)输入交互
此过程主要完成对输入的交互部分。将模型1与模型2的滤波输出通过式(24)的交互原则,转换为滤波器的输入部分。则各滤波器的输入为

X ^ 0 j(k-1/k-1)= i = 1 2 X ^ i(k-1/k-1)uij(k-1/k-1),j=1,2

P0j(k-1/k-1)=∑uij(k-1/k-1){Pi(k-1/k-1)+

[ X ˙ i(k-1/k-1)- X ^ 0 j(k-1/k-1)][ X ^ i(k-1/k-1)- X ^ 0 j(k-1/k-1)]T} (25)

其中,uij(k-1/k-1)=pijui(k-1)/ c ̅ j, c ̅ j= i = 1 2pijui(k-1)。
2)滤波
对应模型j进行卡尔曼滤波,状态预测值为

X ^ j(k/k-1)=Φj X ^ 0 j(k-1/k-1)

状态预测误差协方差为

Pj(k/k-1)=ΦjP0j(k-1/k-1) Φ j T+Qj (27)

卡尔曼增益为

Kj=Pj(k/k-1) H j T S j - 1(k)

k时刻的滤波值为

X ^ j(k/k)= X ^ j(k/k-1)+Kj(k)[Y(k)-Hj X ˙ j(k/k-1)]

滤波协方差为

Pj(k/k)=[I-Kj(k)Hj]Pj(k/k-1)

3)概率更新
模型概率更新:

uj(k)= 1 cΛi(k) i = 1 2Pijui(k-1)=Λi(k) c ̅ j/c

以上式中,

Sj(k)=HjPj(k/k-1) H j T+R

式中,c为归一化常数,且c= j = 1 2Λj(k)/ c ̅ j
Λj(k)为观测Y(k)的似然函数:

Λj(k)= 1 2 π | S j ( k ) | exp - 1 2 γ T j ( k ) S j - 1 ( k ) γ j ( k )

4)输出交互
此部分主要对多模型滤波器的输出量进行交互。自适应交互滤波模型中,交互后的滤波输出为相应时刻自适应坐标系中的滤波输出值,以待输入交互滤波,进入下一轮循环。此过程中应注意,交互前的滤波值应有两种流向,用于相应时刻多模型滤波器的交互输出和转换为下一时刻自适应坐标系中的多模型滤波器的交互输入。

X ^(k/k)= j = 1 2 X ^ j(k/k)uj(k) (34)

P(k/k)= i = 1 2uj(k){Pj(k/k)+[ X ^ j(k/k)- X ^(k/k)][ X ^ j(k/k)- X ^(k/k)]T}

5)坐标转换
任一时刻滤波器的输出转换包括目标状态的转换和误差协方差的转换。滤波后的目标状态变换为

x 1 k + 1 ( k / k ) y 1 k + 1 ( k / k ) z 1 k + 1 ( k / k )= C k k + 1 x 1 k ( k / k ) y 1 k ( k / k ) z 1 k ( k / k )

x ˙ 1 k + 1 ( k / k ) y ˙ 1 k + 1 ( k / k ) z ˙ 1 k + 1 ( k / k )= C k k + 1 x ˙ 1 k ( k / k ) y ˙ 1 k ( k / k ) z ˙ 1 k ( k / k )

x ¨ 1 k + 1 ( k / k ) y ¨ 1 k + 1 ( k / k ) z ¨ 1 k + 1 ( k / k )= C k k + 1 x ¨ 1 k ( k / k ) y ¨ 1 k ( k / k ) z ¨ 1 k ( k / k )

上述转换为将模型1滤波器的状态输出转换到下一时刻自适应坐标系中。类似地,可将模型2滤波器的状态输出转换到下一时刻的自适应坐标系中,此处不再详述。
滤波后的误差协方差转换到下一时刻自适应坐标系中仅作为下一时刻滤波器输入交互使用,在该时刻的滤波器误差协方差输出交互过程中应使用该时刻的误差协方差。下面以模型1滤波器为例对误差协方差的转换进行介绍,模型2的误差协方差转换与模拟1类似。
对于x通道误差协方差转换,有

P 11 P 12 P 13 P 21 P 22 P 23 P 31 P 32 P 33 x 1 k + 1= P 11 P 12 P 13 P 21 P 22 P 23 P 31 P 32 P 33 x 1 k C 11 2+ P 11 P 12 P 13 P 21 P 22 P 23 P 31 P 32 P 33 y 1 k C 12 2+ P 11 P 12 P 13 P 21 P 22 P 23 P 31 P 32 P 33 z 1 k C 13 2

对于y通道误差协方差转换,有

P 11 P 12 P 13 P 21 P 22 P 23 P 31 P 32 P 33 y 1 k + 1= P 11 P 12 P 13 P 21 P 22 P 23 P 31 P 32 P 33 x 1 k C 21 2+ P 11 P 12 P 13 P 21 P 22 P 23 P 31 P 32 P 33 y 1 k C 22 2+ P 11 P 12 P 13 P 21 P 22 P 23 P 31 P 32 P 33 y 1 k C 23 2

对于z通道误差协方差转换,有

P 11 P 12 P 13 P 21 P 22 P 23 P 31 P 32 P 33 z 1 k + 1= P 11 P 12 P 13 P 21 P 22 P 23 P 31 P 32 P 33 x 1 k C 31 2+ P 11 P 12 P 13 P 21 P 22 P 23 P 31 P 32 P 33 y 1 k C 32 2+ P 11 P 12 P 13 P 21 P 22 P 23 P 31 P 32 P 33 z 1 k C 33 2

其中 C k k + 1= C 11 C 12 C 13 C 21 C 22 C 23 C 31 C 32 C 33k时刻瞄准线坐标系转换到k+1时刻瞄准线坐标系的转换矩阵,其转换过程为先将目标状态由k时刻的瞄准线坐标系转换到甲板坐标系,然后再变换为k+1时刻的瞄准线坐标系,故 C k k + 1= C d a(k)T C d a(k+1)。
其中,k时刻甲板坐标系到瞄准线坐标系变换由跟踪传感器瞄准线的方位和俯仰角变换两个步骤组成,具体的变换步骤为
1) 甲板坐标系OXdYdZd经由方位角变换(变换矩阵Mq:绕OZd轴逆向旋转方位角q)为坐标系OXcYcZd;
2) 坐标系OXcYcZd经由高低角变换(变换矩阵Mε:绕OXc轴正向旋转高低角ε)为坐标系OXaYaZa
综合可得k时刻从甲板坐标系到瞄准线坐标系变换矩阵为
C d a (k)=Cε·Cq= cos q - sin q 0 cos ε sin q cos ε cos q sin ε - sin ε sin q - sin ε cos q cos ε

4 仿真计算

设定目标航路参数:目标距观测点的初始距离约为10 000 m,海拔高度为50 m,航路捷径为300 m,采样周期T=0.02 s。目标沿北偏东60°以300 m/s的速度匀速直线飞行15 s,随后以20 m/s2的加速度直线飞行,量测噪声σD=5 m,σq=σε=0.5 mrad, σ ω xa= σ ω ya= σ ω za=0.1 mrad/s。在已知理论值的情况下,统计本文算法与标准的单一模型卡尔曼滤波算法的滤波收敛时间及滤波误差,横轴表示时间,单位为s,纵轴表示位置,单位为m。图2为两算法在舰艇地理坐标系X轴、Y轴的位置滤波效果图。由图2可以看出,单一模型卡尔曼滤波算法收敛时间约为1.98 s,本文算法收敛时间约为0.22 s,后者算法收敛速度明显优于前者。通过表1可知,后者算法滤波精度优于前者。图3为目标机动过程中的两种算法滤波效果图。由图3可以看出,目标第15 s开始机动后,前者算法已开始滤波发散,后者保持对机动目标的平稳滤波。由图2图3表1对比可知,在机动目标跟踪过程中,本文算法的滤波效果明显优于单一模型卡尔曼滤波算法。
图2 两种模型算法滤波对比图
图3 目标机动情形下两种模型算法滤波对比图
表1 单一模型卡尔曼滤波算法与自适应交互多模型算法对比
X轴位置误差均
方差/m
Y轴位置误差均
方差/m
单一模型卡尔曼滤波算法 2.436 7 1.623 6
自适应交互多模型算法 0.626 8 0.364 5

5 结束语

针对新型反舰导弹等高速机动目标跟踪难题,本文在分析噪声特性和挖掘目标量测信息的基础上,推导了目标运动角速度转化为线速度方法,并在瞄准线坐标系中对机动目标进行建模,提出了一种基于目标角速度的自适应交互多模型滤波算法。仿真结果表明,与单一模型卡尔曼滤波算法相比,本文算法可有效提高机动目标滤波性能,缩短滤波收敛时间。不过,本文算法也有缺点,每次循环过程中均需对目标状态和误差协方差进行变换,增加了计算量。
[1]
徐国亮, 王勇. 舰炮反导火控原理[M]. 北京: 北京理工大学出版社, 2018.

[2]
姜青山, 张福光, 汤波. 超音速反舰导弹装备发展现状及使用特点分析[J]. 海军航空工程学院学报, 2002, 17(6):681-682.

[3]
周宏仁, 敬忠良, 王培德. 机动目标跟踪[M]. 2版. 北京: 国防工业出版社, 1994.

[4]
Mazor E, Averbuch A, Bar-Shalom Y. Interacting Multiple Model Methods in Target Tracking: a Survey[J]. IEEE Transactions on Aerospace and Electronic Systems, 1998, 34(1):103-123.

DOI

[5]
Blom H A, Bar-Shalom Y. The Interactive Multiple Model Algorithm for System with Markov Switching Coefficients[J]. IEEE Transactions on Automatic Control, 1988, 33(8):780-783.

DOI

[6]
李辉, 沈莹, 张安, 等. 交互式多模型目标跟踪的研究现状及发展趋势[J]. 火力与指挥控制, 2006, 31(11):1-4.

[7]
陈利斌, 佟明安. 机动目标跟踪的交互式多模型自适应滤波算法[J]. 火力与指挥控制, 2000, 25(4):36-38,42.

[8]
Fitzgerald R J. Simple Tracking Filters Position and Velocity Measurements[J]. IEEE Transactions on Aerospace and Electronic Systems, 1982, 18(5):531-537.

[9]
Zhou Hongren. Tracking of Maneuvering Targets[D]. Minneapolis:University of Minnesota, 1984.

[10]
张怀根, 张林让, 吴君顺. 利用径向速度观测值提高目标跟踪性能[J]. 西安电子科技大学学报(自然科学版), 2005, 32(5):667-670.

[11]
王勇, 梁燊, 徐国亮. 基于目标运动角速度的火控滤波技术[J]. 现代防御技术, 2015, 43(5):124-128.

Outlines

/