CN106199580A - 一种基于模糊推理系统的Singer模型改进算法 - Google Patents

Abstract

本发明公开了一种基于模糊推理系统的Singer模型改进算法，属于雷达机动目标跟踪领域。本发明的方法可改善强机动条件下Singer模型跟踪性能不高的问题，有效提高雷达对机动目标的跟踪性能。本发明的方法包括以下步骤：（一）对雷达接收信号进行采样，利用Wigner‑Hough变换方法得到目标径向加速度和径向速度；（二）在数据处理阶段将径向加速度和径向速度引入模糊推理系统；（三）根据模糊推理系统的输出对Singer模型的状态向量预测协方差进行实时调整，从而实现机动目标的精确跟踪。本发明能够精确实时地反映目标的机动情况，提高了目标跟踪的精度，在速度和加速度估计精度方面也有明显改善，工程实现容易，具有较强的工程应用价值和推广前景。

Description

一种基于模糊推理系统的Singer模型改进算法

一、技术领域

本发明隶属于雷达机动目标跟踪处理领域，适用于高/中脉冲重复频率雷达(如机载脉冲多普勒雷达等)对机动目标的精确跟踪。

二、背景技术

随着科学技术的发展，现代飞机和导弹的机动能力大大增强，对雷达的检测和跟踪性能提出了新的挑战，对高机动目标跟踪也成为雷达目标跟踪领域的难点和重点。现有的机动目标跟踪技术是在雷达测量信息(位置、多普勒速度)的基础上，利用各种机动模型实现对机动目标的跟踪，其中，Singer模型是一种常用的有效机动目标跟踪算法，主要由以下3个步骤实现：

(1)将雷达接收机输出的目标回波信号进行A/D变换，送雷达数据处理计算机执行以下步骤；

(2)建立机动目标模型，假设机动加速度的“当前”概率密度用修正的瑞利分布描述，均值为“当前”加速度预测值，随机机动加速度在时间轴上符合一阶时间相关过程；

(3)采用扩展卡尔曼滤波算法对机动目标进行自适应跟踪，估计出目标的位置、速度和加速度状态。

这种方法具有以下缺陷：

(1)当目标加速度发生急剧变化时，存在着跟踪精度低、动态时延明显较大的问题；

(2)速度和加速度估计的误差较大，不能精确实时地反映目标的机动情况。

三、发明内容

本发明的目的是提出一种基于模糊推理系统的Singer模型改进算法，解决现有Singer模型算法对目标突发机动的自适应跟踪能力不强，以及速度和加速度估计误差较大的问题。

本发明提出的基于模糊推理Singer模型机动目标跟踪算法的技术方案包括以下步骤：

步骤1：将雷达接收机接收到的线性调频信号s(t)通过采样器以采样间隔T_s进行采样，变为离散信号s(nT_s)，其中n表示采样点序号；将s(nT_s)送入雷达信号处理计算机；

在雷达信号处理计算机中执行以下步骤：

步骤2：初始化

λ设为雷达波长；

T_s设为信号处理阶段采样间隔；

T′设为信号观测时间；

SNR设为信噪比；

T_s′设为数据处理阶段采样间隔；

σ_ρ设为距离量测误差；

σ_θ设为角度量测误差；

σ_v设为径向速度量测误差；

σ_a设为径向加速度量测误差；

X(0)设为目标起始状态；

u₁、u₂分别设为模糊推理系统中目标径向速度和径向加速度的变化率；

α设为Singer模型中的机动时间常数；

步骤3：将采样后的离散信号送入雷达信号处理计算机求取其Wigner-ville分布，并对该分布进行Hough变换，搜索变换图峰值，根据峰值坐标确定信号的径向速度和径向加速度

将目标的位置、径向速度和径向加速度信息送到雷达的数据处理计算机，在雷达数据处理计算机中执行以下步骤：

步骤4：建立目标运动方程

考虑二坐标雷达数据处理问题，直角坐标系下目标运动模型表示为：

X(k)＝Φ(k,k-1)X(k-1)+w(k-1)

其中，状态向量为：

$X\left(k\right)={\left[\begin{array}{cccccc}{x}_k& {\stackrel{\·}{x}}_k& {\stackrel{\·\·}{x}}_k& y_k& {\stackrel{\·}{y}}_k& {\stackrel{\·\·}{y}}_k\end{array}\right]}^T$

式中，x_k、和y_k、分别表示k时刻目标x方向和y方向的位置、速度、加速度。

Φ为状态转移矩阵，w(k-1)为零均值方差为Q的高斯白噪声。

状态转移矩阵Φ为：

$\left\{\begin{array}{c}\mathrm{\Φ}(k,k-1)=diag\left({\mathrm{\Φ}}_{11}\right(k,k-1),{\mathrm{\Φ}}_{22}(k,k-1\left)\right)\\ {\mathrm{\Φ}}_{11}(k,k-1)={\mathrm{\Φ}}_{22}(k,k-1)=\left[\begin{array}{ccc}1& T_s& (e^{-\mathrm{\α}T}+\mathrm{\α}T-1)/{\mathrm{\α}}^2\\ 0& 1& (1-e^{-\mathrm{\α}T})/\mathrm{\α}\\ 0& 0& e^{-\mathrm{\α}T}\end{array}\right]\end{array}\right.$

式中，α为机动时间常数；T_s为采样周期。

离散时间过程噪声w具有协方差:

$Q=2{\mathrm{\α\σ}}_a^2\left[\begin{array}{ccc}{q}_{11}& q_{12}& q_{13}\\ q_{12}& q_{22}& q_{23}\\ q_{13}& q_{23}& q_{33}\end{array}\right]$

式中，Q精确表达式为：

$q_{11}=\frac{1}{2{\mathrm{\α}}^5}\[1-e^{-2\mathrm{\α}T}+2\mathrm{\α}T+\frac{2{\mathrm{\α}}^3{T}^3}{3}-2{\mathrm{\α}}^2{T}^2-4{\mathrm{\αTe}}^{-\mathrm{\α}T}\]$

$q_{12}=\frac{1}{2{\mathrm{\α}}^4}\[e^{-2\mathrm{\α}T}+1-2e^{-\mathrm{\α}T}+2{\mathrm{\αTe}}^{-\mathrm{\α}T}-2\mathrm{\α}T+{\mathrm{\α}}^2{T}^2\]$

$q_{13}=\frac{1}{2{\mathrm{\α}}^3}\[1-e^{-2\mathrm{\α}T}-2{\mathrm{\αTe}}^{-\mathrm{\α}T}\]$

$q_{22}=\frac{1}{2{\mathrm{\α}}^3}\[4e^{-\mathrm{\α}T}-3-e^{-2\mathrm{\α}T}+2\mathrm{\α}T\]$

$q_{23}=\frac{1}{2{\mathrm{\α}}^2}\[e^{-2\mathrm{\α}T}+1-2e^{-\mathrm{\α}T}\]$

$q_{33}=\frac{1}{2\mathrm{\α}}\[1-e^{-2\mathrm{\α}T}\]$

步骤5：建立目标量测方程

由于量测向量和状态向量之间存在非线性强度高的缺点，因此本专利优先处理位置量测值，对其进行无偏转换卡尔曼滤波，然后使径向速度和径向加速度依据状态向量的估计值进行EKF处理，这样既可以降低了线性化误差，又可以保证滤波的稳定性。

普通雷达在极坐标系下可获得位置量测(距离和方位角)，利用Wigner-Hough变换方法可以得到径向速度量测和径向加速度量测，因此测量向量表示为

$\begin{array}{c}{Z}_k^m={\left[\begin{array}{cccc}{r}_k^m& {\mathrm{\θ}}_k^m& v_{rk}^m& a_{rk}^m\end{array}\right]}^T\\ ={\left[\begin{array}{cccc}r& \mathrm{\θ}& v& a\end{array}\right]}^T+{\left[\begin{array}{cccc}\tilde{r}& \widetilde{\mathrm{\θ}}& \tilde{v}& \tilde{a}\end{array}\right]}^T\end{array}$

其中

$\left\{\begin{array}{c}{r}_k=\sqrt{x_k^2+y_k^2}\\ {\mathrm{\θ}}_k=arctan(y_k/x_k)\end{array}\right.$

上式中，为k时刻从雷达信号处理阶段获得的测量向量，其中包括距离量测方位角量测径向速度量测和径向加速度量测r、θ、v和a分别为距离、方位角、径向速度和径向加速度的真值。和为相应的量测误差，假设它们都是均值为零的高斯白噪声序列且相互统计独立，方差分别为和

(1)位置转换测量方程

从极坐标到直角坐标的位置转换测量方程可表示为

$\begin{array}{c}{Z}_1\left(k\right)=H_L\left(k\right)X\left(k\right)\\ =\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\end{array}\right]\left[\begin{array}{c}x\\ y\\ \stackrel{\·}{x}\\ \stackrel{\·}{y}\\ \stackrel{\·\·}{x}\\ \stackrel{\·\·}{y}\end{array}\right]\end{array}$

Z₁(k)为k时刻目标的量测，转换测量协方差为：

$R_1=\left[\begin{array}{cc}{\mathrm{\σ}}_x^2& {\mathrm{\σ}}_{xy}\\ {\mathrm{\σ}}_{xy}& {\mathrm{\σ}}_y^2\end{array}\right]$

其中，

${\mathrm{\σ}}_x^2=\frac{1}{2}(r^2+{\mathrm{\σ}}_{\mathrm{\ρ}}^2)(1+{\mathrm{\λ}}_{\mathrm{\θ}}^{\′}cos2\mathrm{\θ})+({\mathrm{\λ}}_{\mathrm{\θ}}^{-2}-2)r^2{\cos }^2\mathrm{\θ}$

${\mathrm{\σ}}_y^2=\frac{1}{2}(r^2+{\mathrm{\σ}}_{\mathrm{\ρ}}^2)(1-{\mathrm{\λ}}_{\mathrm{\θ}}^{\′}cos2\mathrm{\θ})+({\mathrm{\λ}}_{\mathrm{\θ}}^{-2}-2)r^2{\sin }^2\mathrm{\θ}$

${\mathrm{\σ}}_{xy}={\mathrm{\σ}}_{yx}=\frac{1}{2}(r^2+{\mathrm{\σ}}_{\mathrm{\ρ}}^2){\mathrm{\λ}}_{\mathrm{\θ}}^{\′}sin2\mathrm{\θ}+({\mathrm{\λ}}_{\mathrm{\θ}}^{-2}-2)r^{2}sin\mathrm{\θ}cos\mathrm{\θ}$

式中，

${\mathrm{\λ}}_{\mathrm{\θ}}^{\′}=e^{-2{\mathrm{\σ}}_{\mathrm{\θ}}^2}={\mathrm{\λ}}_{\mathrm{\θ}}^4$

(2)径向加速度和径向速度量测方程

从极坐标到直角坐标的径向加速度和径向速度量测方程进行线性化处理得到：

$\begin{array}{c}{Z}_2\left(k\right)=h(k,X(k))+W\left(k\right)\\ =\left[\begin{array}{cccccc}{S}_1\left(k\right)& S_2\left(k\right)& 0& S_3\left(k\right)& S_4\left(k\right)& 0\\ S_5\left(k\right)& 0& S_6\left(k\right)& S_7\left(k\right)& 0& S_8\left(k\right)\end{array}\right]\left[\begin{array}{c}x\\ y\\ \stackrel{\·}{x}\\ \stackrel{\·}{y}\\ \stackrel{\·\·}{x}\\ \stackrel{\·\·}{y}\end{array}\right]+W\left(k\right)\end{array}$

式中，W(k)为量测噪声，假定为加性零均值白噪声，

$S_1\left(k\right){|}_{X=X\left(k\right|k)}=\frac{\stackrel{\·}{x}ysin\mathrm{\θ}}{x^2+y^2}-\frac{y\stackrel{\·}{y}cos\mathrm{\θ}}{x^2+y^2}$

S₂(k)|_X＝X(k|k)＝cosθ

$S_3\left(k\right){|}_{X=X\left(k\right|k)}=-\frac{\stackrel{\·}{x}xsin\mathrm{\θ}}{x^2+y^2}+\frac{x\stackrel{\·}{y}cos\mathrm{\θ}}{x^2+y^2}$

S₄(k)|_X＝X(k|k)＝sinθ

$\begin{array}{c}{S}_5\left(k\right){|}_{X=X\left(k\right|k)}=\frac{\stackrel{\·\·}{x}ysin\mathrm{\θ}}{x^2+y^2}-\frac{y\stackrel{\·\·}{y}cos\mathrm{\θ}}{x^2+y^2}-\frac{x\stackrel{\·}{x}{(\stackrel{\·}{x}sin\mathrm{\θ}-\stackrel{\·}{y}cos\mathrm{\θ})}^2}{(x^2+y^2)\sqrt{x^2+y^2}}\\ +\frac{2(\stackrel{\·}{x}\sin \mathrm{\θ}-\stackrel{\·}{y}\cos \mathrm{\θ})(-y\stackrel{\·}{x}\cos \mathrm{\θ}-y\stackrel{\·}{y}\sin \mathrm{\θ})}{(x^2+y^2)\sqrt{x^2+y^2}}\end{array}$

S₆(k)|_X＝X(k|k)＝cosθ

$\begin{array}{c}{S}_7\left(k\right){|}_{X=X\left(k\right|k)}=-\frac{\stackrel{\·\·}{x}xsin\mathrm{\θ}}{x^2+y^2}+\frac{x\stackrel{\·\·}{y}cos\mathrm{\θ}}{x^2+y^2}-\frac{y\stackrel{\·}{y}{(\stackrel{\·}{x}sin\mathrm{\θ}-\stackrel{\·}{y}cos\mathrm{\θ})}^2}{(x^2+y^2)\sqrt{x^2+y^2}}\\ +\frac{2(\stackrel{\·}{x}\sin \mathrm{\θ}-\stackrel{\·}{y}\cos \mathrm{\θ})(x\stackrel{\·}{x}\cos \mathrm{\θ}+x\stackrel{\·}{y}\sin \mathrm{\θ})}{(x^2+y^2)\sqrt{x^2+y^2}}\end{array}$

S₈(k)|_X＝X(k|k)＝sinθ

步骤6：模糊推理系统

F(x)＝[u₁,u₂]，u₁,u₂分别表示目标的径向速度和径向加速度信息，由于量测向量是二维的，u₁,u₂分别定义为：

$u_1=\left|{\mathrm{\Δv}}_{rk}^m\right|=|v_{rk}^m-v_{rk-1}^m|,u_2=\left|{\mathrm{\Δa}}_{rk}^m\right|=|a_{rk}^m-a_{rk-1}^m|$

其中，和分别代表k时刻信号处理阶段得到的径向速度和径向加速度值，u₁,u₂分别代表k时刻目标径向速度和径向加速度的变化率，它们的值越大说明目标机动性越强，反之，机动性越小。模糊推理系统的输入可以划分为多个模糊子集，模糊子集的数量越多，输出的精度越高，但相应的计算量也会增加。综合考虑，本专利将u₁,u₂划分为七个子集，分别为[A₁,A₂,A₃,A₄,A₅,A₆,A₇]、[B₁,B₂,B₃,B₄,B₅,B₆,B₇]，其隶属度函数可以通过模糊统计方法来确定，本专利选用常用论域[-6,-5,-4,-3,-2,-1,0,1，2,3,4,5,6]对两个输入变量u₁,u₂进行模糊化处理，定义上有7个语言变量值分别为NB、NM、NS、0、PS、PM和PB(负大、负中、负小、零、正小、正中和正大)，公式如下：

$U=\frac{12}{b-a}\[u-\frac{a+b}{2}\]$

其中，a，b代表输入变量的取值范围，将区间为[a，b]的输入量u₁，u₂转化为[-6，6]区间的变化量U₁，U₂。u₁，u₂的隶属度函数μ(U₁)，μ(U₂)采用三角函数，该函数形状简单、计算量小，当输入量变化时具有较好的灵活性。假设u₁，u₂的取值范围分别为[0，60]和[0，90]，其隶属度函数分别如图2和3所示，模糊子集的计算公式为：

$A={\∫}_{U_1}\mathrm{\μ}\left(U_1\right){U_1}^{-1}{\mathrm{dU}}_1$

模糊推理系统输出f是一个调整因子，可以根据目标的机动情况自动调节滤波器的预测协方差。理论上f的范围可以很大，根据经验及试验统计分析，将f的取值范围定义在区间[0，8]，将其划分为五个模糊子集[C₁，C₂，C₃，C₄，C₅]，定义上有其个语言变量值分别为为NM、NS、0、PS、和PM，隶属度函数也采用三角形函数，如图4所示。

当所有输入特征量处于低值时，认为目标处于弱机动目标状态，f取较小的值；当所有特征量处于极大值时，则认为目标处于强机动目标状态，f取较大的值，以保证跟踪精度，从而本专利建立模糊推理规则如图5，其中×为不存在项，这些状态在实际目标运动过程中不出现，因此规则数目相对减少。计算输入量和输出量的模糊关系R，公式如下：

R＝(A₁×B₁×C₁)∪(A₁×B₂×C₁)∪…∪(A₇×B₇×C₅)

对某一时刻的输入u₁，u₂，与模糊关系R进行合成运算，得到输出C：

$\begin{array}{c}C=A^{\′}\×B^{\′}\×R\\ =\left[\begin{array}{c}\mathrm{\μ}\left(U_{11}\right)\\ \mathrm{\μ}\left(U_{12}\right)\\ .\\ .\\ .\\ \mathrm{\μ}\left(U_{17}\right)\end{array}\right]\×\left[\begin{array}{cccc}\mathrm{\μ}\left(U_{21}\right)& \mathrm{\μ}\left(U_{22}\right)& \mathrm{...}& \mathrm{\μ}\left(U_{27}\right)\end{array}\right]\×R\end{array}$

其中，U_1i，U_2i(i＝1，2，…，7)分别代表某一时刻输入量u₁，u₂在其模糊子集上的取值。最后将模糊量精确化，把经过模糊推理得到的以模糊子集表示的调节量C，根据隶属函数计算出最终输出精确值f，这里采用面积重心法进行反模糊化，重心法的计算式为：

$f=\frac{\mathrm{\Σ}\mathrm{\μ}\left(f_i\right)f_i}{\mathrm{\Σ}\mathrm{\μ}\left(f_i\right)}$

步骤6：自适应滤波算法

在滤波算法中，先将极坐标系下位置量测无偏转换到直角坐标系下，采用标准的卡尔曼滤波算法进行跟踪滤波，再利用状态滤波值对带径向加速度和径向速度的量测方程进行扩展卡尔曼滤波，并将模糊推理系统得到的f对状态向量预测协方差进行实时调整，其具体流程可以表示为：

步骤1：时间更新滤波估计

1)状态向量一步预测为

X(k)＝Φ(k|k-1)X(k-1)+W(k-1)

2)状态向量预测协方差为

P(k|k-1)＝f×Φ(k|k-1)P(k-1|k-1)+Q(k-1)

步骤2：位置量测更新滤波估计

1)量测向量一步预测为

Z₁(k|k-1)＝H₁(k)X(k|k-1)

式中，

2)新息为

In(k)＝Z₁(k)-Z₁(k|k-1)

3)新息协方差为

S(k)＝H₁(k)P(k|k-1)H₁′(k)+R₁(k)

4)卡尔曼增益为

K(k)＝P(k|k-1)H₁(k)′S^-1(k)

5)状态向量一步更新为

X(k|k)＝X(k|k-1)+K(k)In(k)

6)状态向量协方差更新为

P(k|k)＝P(k|k-1)-K(k)S(k)K′(k)

步骤3：径向速度和径向加速度量测更新滤波估计

1)量测向量一步预测为

Z₂(k|k-1)≈H₂(k)X(k-1)

2)量测雅克比矩阵为

$H_2\left(k\right)=\left[\begin{array}{cccccc}{S}_1\left(k\right)& S_2\left(k\right)& 0& S_3\left(k\right)& S_4\left(k\right)& 0\\ S_5\left(k\right)& 0& S_6\left(k\right)& S_7\left(k\right)& 0& S_8\left(k\right)\end{array}\right]$

3)新息为

In(k)＝Z₂(k)-Z₂(k|k-1)

4)新息协方差为

S₂(k)＝H₂(k)P₂(k|k-1)H₂′(k)

5)卡尔曼增益为

K₂(k)＝P(k|k-1)H₂(k)′(k)S^-1(k)

6)状态向量一步更新为

X₂(k|k)＝X₂(k|k-1)+K₂(k)In₂(k)

7)状态向量协方差更新为

P₂(k|k)＝P(k|k-1)-K₂(K)S₂(k)K′₂(k)

步骤4：最终滤波估计

1)状态向量更新为

X(k|k)＝X₂(k|k)

2)状态向量协方差更新为

P(k|k)＝P₂(k|k)

当输出量f小于1时，说明目标处于弱机动状态，这时令f为1，机动模型退化为Singer模型。

和现有技术相比，本发明的有益效果说明：

(1)本发明方法相比现有技术相比，在目标机动运动阶段跟踪过程更加稳定，收敛速度和跟踪精度都有所提高；

(2)本发明方法相比现有技术在滤波后X方向速度和加速度估计的精度也有所改善。

附图说明

附图1是本发明提出的基于模糊推理系统的Singer模型机动目标跟踪算法整体流程图；

附图2是本发明模糊推理系统中输入量u₁的隶属度函数图；

附图3是本发明模糊推理系统中输入量u₂的隶属度函数图；

附图4是本发明模糊推理系统中输出量f的隶属度函数图；

附图5是本发明模糊推理规则；

附图6是目标运动轨迹；

附图7是基于传统Singer模型和基于本发明算法的目标径向距离误差曲线对比图；

附图8是基于传统Singer模型和基于本发明算法的目标X轴位置误差曲线；

附图9是基于传统Singer模型和基于本发明算法的目标Y轴位置误差曲线；

附图10是基于传统Singer模型和基于本发明算法的目标X轴速度误差曲线；

附图11是基于传统Singer模型和基于本发明算法的目标Y轴速度误差曲线；

附图12是基于传统Singer模型和基于本发明算法的目标X轴加速度误差曲线；

附图13是基于传统Singer模型和基于本发明算法的目标Y轴加速度误差曲线。

五、具体实施方式

下面结合附图对本发明的方法进行详细描述。

实施例条件：假设雷达波长λ＝8mm，雷达在信号处理阶段采样间隔为T_s＝5×10^{- 5}s，信号观测时间为T＝20ms，信噪比为SNR＝0dB；雷达在数据处理阶段采样间隔为T_s′＝1s，距离量测误差σ_r＝100m，角度量测误差径向速度量测误差σ_v＝1m/s，径向加速度量测误差σ_a＝1m/s²，目标起始状态：X(0)＝[120000m,-426m/s,0m/s²,2000m,0m/s,0m/s²]^T，目标运动过程历时90s，模糊推理系统中设u₁,u₂的取值范围分别为u₁∈[0,60]，u₂∈[0,90]，目标发生机动时刻及加速度大小如表1所示，目标直角坐标系下运动轨迹如图6所示。Singer模型参数设定：α＝1/20，a_max＝100m/s²，P_max＝0.95，P_min＝0.05。

表1目标机动运动情况

目标发生机动时刻(s)	t＝31	t＝38	t＝49	t＝61	t＝65	t＝66	t＝81
								X方向加速度值(m/s2)	5	-8	10	0	-10	-5	5
Y方向加速度值(m/s2)	-10	18	-20	30	-8	0	-10

将雷达接收机接收到的模拟信号送入雷达处理计算机中执行以下步骤(参照说明书附图1)：

步骤1：将雷达接收机接收到的线性调频信号s(t)通过采样器以采样间隔T_s进行采样，变为离散信号s(nT_s)，其中n表示采样点序号；将s(nT_s)送入雷达信号处理计算机；

步骤2：初始化

雷达波长λ设为8mm；

信号处理阶段采样间隔T_s设为5×10^-5s；

信号观测时间T设为20ms；

信噪比SNR设为0dB；

数据处理阶段采样间隔T_s′设为1s；

距离量测误差σ_r设为100m；

角度量测误差设为0.5°；

径向速度量测误差σ_v设为1m/s；

径向加速度量测误差σ_a设为1m/s²；

目标起始状态X(0)设为[120000m,-426m/s,0m/s²,2000m,0m/s,0m/s²]^T；

目标运动过程设为90s；

模糊推理系统u₁,u₂取值范围分别设为u₁∈[0,60]，u₂∈[0,90]；

Singer模型参数α设为1/20；

步骤3：将采样后的离散信号送入雷达信号处理计算机求取其Wigner-ville分布，并对该分布进行Hough变换，搜索变换图峰值，根据峰值坐标确定信号的径向速度和径向加速度

将目标的位置、径向速度和径向加速度信息送到雷达的数据处理计算机，在雷达数据处理计算机中执行以下步骤：

步骤4：建立目标运动方程

考虑二坐标雷达数据处理问题，直角坐标系下目标运动模型表示为：

X(k)＝Φ(k,k-1)X(k-1)+w(k-1)

其中，状态向量为：

$X\left(k\right)={\left[\begin{array}{cccccc}{x}_k& {\stackrel{\·}{x}}_k& {\stackrel{\·\·}{x}}_k& y_k& {\stackrel{\·}{y}}_k& {\stackrel{\·\·}{y}}_k\end{array}\right]}^T$

式中，x_k、和y_k、分别表示k时刻目标x方向和y方向的位置、速度、加速度。

Φ为状态转移矩阵，w(k-1)为零均值方差为Q的高斯白噪声。

状态转移矩阵Φ为：

$\left\{\begin{array}{c}\mathrm{\Φ}(k,k-1)=diag\left({\mathrm{\Φ}}_{11}\right(k,k-1),{\mathrm{\Φ}}_{22}(k,k-1\left)\right)\\ {\mathrm{\Φ}}_{11}(k,k-1)={\mathrm{\Φ}}_{22}(k,k-1)=\left[\begin{array}{ccc}1& T_s& (e^{-\mathrm{\α}T}+\mathrm{\α}T-1)/{\mathrm{\α}}^2\\ 0& 1& (1-e^{-\mathrm{\α}T})/\mathrm{\α}\\ 0& 0& e^{-\mathrm{\α}T}\end{array}\right]\end{array}\right.$

式中，α为机动时间常数；T_s为采样周期。

离散时间过程噪声w具有协方差:

$Q=2{\mathrm{\α\σ}}_a^2\left[\begin{array}{ccc}{q}_{11}& q_{12}& q_{13}\\ q_{12}& q_{22}& q_{23}\\ q_{13}& q_{23}& q_{33}\end{array}\right]$

式中，Q精确表达式为：

$q_{11}=\frac{1}{2{\mathrm{\α}}^5}\[1-e^{-2\mathrm{\α}T}+2\mathrm{\α}T+\frac{2{\mathrm{\α}}^3{T}^3}{3}-2{\mathrm{\α}}^2{T}^2-4{\mathrm{\αTe}}^{-\mathrm{\α}T}\]$

$q_{12}=\frac{1}{2{\mathrm{\α}}^4}\[e^{-2\mathrm{\α}T}+1-2e^{-\mathrm{\α}T}+2{\mathrm{\αTe}}^{-\mathrm{\α}T}-2\mathrm{\α}T+{\mathrm{\α}}^2{T}^2\]$

$q_{13}=\frac{1}{2{\mathrm{\α}}^3}\[1-e^{-2\mathrm{\α}T}-2{\mathrm{\αTe}}^{-\mathrm{\α}T}\]$

$q_{22}=\frac{1}{2{\mathrm{\α}}^3}\[4e^{-\mathrm{\α}T}-3-e^{-2\mathrm{\α}T}+2\mathrm{\α}T\]$

$q_{23}=\frac{1}{2{\mathrm{\α}}^2}\[e^{-2\mathrm{\α}T}+1-2e^{-\mathrm{\α}T}\]$

$q_{33}=\frac{1}{2\mathrm{\α}}\[1-e^{-2\mathrm{\α}T}\]$

步骤5：建立目标量测方程

由于量测向量和状态向量之间存在非线性强度高的缺点，因此本专利优先处理位置量测值，对其进行无偏转换卡尔曼滤波，然后使径向速度和径向加速度依据状态向量的估计值进行EKF处理，这样既可以降低了线性化误差，又可以保证滤波的稳定性。

普通雷达在极坐标系下可获得位置量测(距离和方位角)，利用Wigner-Hough变换方法可以得到径向速度量测和径向加速度量测，因此测量向量表示为

$\begin{array}{c}{Z}_k^m={\left[\begin{array}{cccc}{r}_k^m& {\mathrm{\θ}}_k^m& v_{rk}^m& a_{rk}^m\end{array}\right]}^T\\ ={\left[\begin{array}{cccc}r& \mathrm{\θ}& v& a\end{array}\right]}^T+{\left[\begin{array}{cccc}\tilde{r}& \widetilde{\mathrm{\θ}}& \tilde{v}& \tilde{a}\end{array}\right]}^T\end{array}$

其中

$\left\{\begin{array}{c}{r}_k=\sqrt{x_k^2+y_k^2}\\ {\mathrm{\θ}}_k=arctan(y_k/x_k)\end{array}\right.$

上式中，为k时刻从雷达信号处理阶段获得的测量向量，其中包括距离量测方位角量测径向速度量测和径向加速度量测r、θ、v和a分别为距离、方位角、径向速度和径向加速度的真值。和为相应的量测误差，假设它们都是均值为零的高斯白噪声序列且相互统计独立，方差分别为和

(1)位置转换测量方程

从极坐标到直角坐标的位置转换测量方程可表示为

$\begin{array}{c}{Z}_1\left(k\right)=H_L\left(k\right)X\left(k\right)\\ =\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\end{array}\right]\left[\begin{array}{c}x\\ y\\ \stackrel{\·}{x}\\ \stackrel{\·}{y}\\ \stackrel{\·\·}{x}\\ \stackrel{\·\·}{y}\end{array}\right]\end{array}$

Z(k)为k时刻目标的量测，转换测量协方差为：

$R_1=\left[\begin{array}{cc}{\mathrm{\σ}}_x^2& {\mathrm{\σ}}_{xy}\\ {\mathrm{\σ}}_{xy}& {\mathrm{\σ}}_y^2\end{array}\right]$

其中，

${\mathrm{\σ}}_x^2=\frac{1}{2}(r^2+{\mathrm{\σ}}_{\mathrm{\ρ}}^2)(1+{\mathrm{\λ}}_{\mathrm{\θ}}^{\′}cos2\mathrm{\θ})+({\mathrm{\λ}}_{\mathrm{\θ}}^{-2}-2)r^2{\cos }^2\mathrm{\θ}$

${\mathrm{\σ}}_y^2=\frac{1}{2}(r^2+{\mathrm{\σ}}_{\mathrm{\ρ}}^2)(1-{\mathrm{\λ}}_{\mathrm{\θ}}^{\′}cos2\mathrm{\θ})+({\mathrm{\λ}}_{\mathrm{\θ}}^{-2}-2)r^2{\sin }^2\mathrm{\θ}$

${\mathrm{\σ}}_{xy}={\mathrm{\σ}}_{yx}=\frac{1}{2}(r^2+{\mathrm{\σ}}_{\mathrm{\ρ}}^2){\mathrm{\λ}}_{\mathrm{\θ}}^{\′}sin2\mathrm{\θ}+({\mathrm{\λ}}_{\mathrm{\θ}}^{-2}-2)r^{2}sin\mathrm{\θ}cos\mathrm{\θ}$

式中，

${\mathrm{\λ}}_{\mathrm{\θ}}^{\′}=e^{-2{\mathrm{\σ}}_{\mathrm{\θ}}^2}={\mathrm{\λ}}_{\mathrm{\θ}}^4$

(2)径向加速度和径向速度量测方程

从极坐标到直角坐标的径向加速度和径向速度量测方程进行线性化处理得到：

$\begin{array}{c}{Z}_2\left(k\right)=h(k,X(k))+W\left(k\right)\\ =\left[\begin{array}{cccccc}{S}_1\left(k\right)& S_2\left(k\right)& 0& S_3\left(k\right)& S_4\left(k\right)& 0\\ S_5\left(k\right)& 0& S_6\left(k\right)& S_7\left(k\right)& 0& S_8\left(k\right)\end{array}\right]\left[\begin{array}{c}x\\ y\\ \stackrel{\·}{x}\\ \stackrel{\·}{y}\\ \stackrel{\·\·}{x}\\ \stackrel{\·\·}{y}\end{array}\right]+W\left(k\right)\end{array}$

式中，W(k)为量测噪声，假定为加性零均值白噪声，

$S_1\left(k\right){|}_{X=X\left(k\right|k)}=\frac{\stackrel{\·}{x}ysin\mathrm{\θ}}{x^2+y^2}-\frac{y\stackrel{\·}{y}cos\mathrm{\θ}}{x^2+y^2}$

S₂(k)|_X＝X(k|k)＝cosθ

$S_3\left(k\right){|}_{X=X\left(k\right|k)}=-\frac{\stackrel{\·}{x}xsin\mathrm{\θ}}{x^2+y^2}+\frac{x\stackrel{\·}{y}cos\mathrm{\θ}}{x^2+y^2}$

S₄(k)|_X＝X(k|k)＝sinθ

$\begin{array}{c}{S}_5\left(k\right){|}_{X=X\left(k\right|k)}=\frac{\stackrel{\·\·}{x}ysin\mathrm{\θ}}{x^2+y^2}-\frac{y\stackrel{\·\·}{y}cos\mathrm{\θ}}{x^2+y^2}-\frac{x\stackrel{\·}{x}{(\stackrel{\·}{x}sin\mathrm{\θ}-\stackrel{\·}{y}cos\mathrm{\θ})}^2}{(x^2+y^2)\sqrt{x^2+y^2}}\\ +\frac{2(\stackrel{\·}{x}\sin \mathrm{\θ}-\stackrel{\·}{y}\cos \mathrm{\θ})(-y\stackrel{\·}{x}\cos \mathrm{\θ}-y\stackrel{\·}{y}\sin \mathrm{\θ})}{(x^2+y^2)\sqrt{x^2+y^2}}\end{array}$

S₆(k)|_X＝X(k|k)＝cosθ

$\begin{array}{c}{S}_7\left(k\right){|}_{X=X\left(k\right|k)}=-\frac{\stackrel{\·\·}{x}xsin\mathrm{\θ}}{x^2+y^2}+\frac{x\stackrel{\·\·}{y}cos\mathrm{\θ}}{x^2+y^2}-\frac{y\stackrel{\·}{y}{(\stackrel{\·}{x}sin\mathrm{\θ}-\stackrel{\·}{y}cos\mathrm{\θ})}^2}{(x^2+y^2)\sqrt{x^2+y^2}}\\ +\frac{2(\stackrel{\·}{x}\sin \mathrm{\θ}-\stackrel{\·}{y}\cos \mathrm{\θ})(x\stackrel{\·}{x}\cos \mathrm{\θ}+x\stackrel{\·}{y}\sin \mathrm{\θ})}{(x^2+y^2)\sqrt{x^2+y^2}}\end{array}$

S₈(k)|_X＝X(k|k)＝sinθ

步骤6：模糊推理系统

F(x)＝[u₁,u₂]，u₁,u₂分别表示目标的径向速度和径向加速度信息，由于量测向量是二维的，u₁,u₂分别定义为：

$u_1=\left|{\mathrm{\Δv}}_{rk}^m\right|=|v_{rk}^m-v_{rk-1}^m|,u_2=\left|{\mathrm{\Δa}}_{rk}^m\right|=|a_{rk}^m-a_{rk-1}^m|$

其中，和分别代表k时刻信号处理阶段得到的径向速度和径向加速度值，u₁，u₂分别代表k时刻目标径向速度和径向加速度的变化率，它们的值越大说明目标机动性越强，反之，机动性越小。模糊推理系统的输入可以划分为多个模糊子集，模糊子集的数量越多，输出的精度越高，但相应的计算量也会增加。综合考虑，本专利将u₁，u₂划分为七个子集，分别为[A₁，A₂，A₃，A₄，A₅，A₆，A₇]、[B₁，B₂，B₃，B₄，B₅，B₆，B₇]，其隶属度函数可以通过模糊统计方法来确定，本专利选用常用论域[-6，-5，-4，-3，-2，-1，0，1，2，3，4，5，6]对两个输入变量u₁，u₂进行模糊化处理，定义上有7个语言变量值分别为NB、NM、NS、0、PS、PM和PB(负大、负中、负小、零、正小、正中和正大)，公式如下：

$U=\frac{12}{b-a}\[u-\frac{a+b}{2}\]$

其中，a，b代表输入变量的取值范围，将区间为[a，b]的输入量u₁，u₂转化为[-6，6]区间的变化量U₁，U₂。u₁，u₂的隶属度函数μ(U₁),μ(U₂)采用三角函数，该函数形状简单、计算量小，当输入量变化时具有较好的灵活性。假设u₁,u₂的取值范围分别为[0,60]和[0,90]，其隶属度函数分别如图2和3所示，模糊子集的计算公式为：

$A={\∫}_{U_1}\mathrm{\μ}\left(U_1\right){U_1}^{-1}{\mathrm{dU}}_1$

模糊推理系统输出f是一个调整因子，可以根据目标的机动情况自动调节滤波器的预测协方差。理论上f的范围可以很大，根据经验及试验统计分析，将f的取值范围定义在区间[0,8]，将其划分为五个模糊子集[C₁,C₂,C₃,C₄,C₅]，定义上有其个语言变量值分别为为NM、NS、0、PS、和PM，隶属度函数也采用三角形函数，如图4所示。

当所有输入特征量处于低值时，认为目标处于弱机动目标状态，f取较小的值；当所有特征量处于极大值时，则认为目标处于强机动目标状态，f取较大的值，以保证跟踪精度，从而本专利建立模糊推理规则如图5，其中×为不存在项，这些状态在实际目标运动过程中不出现，因此规则数目相对减少。计算输入量和输出量的模糊关系R，公式如下：

R＝(A₁×B₁×C₁)∪(A₁×B₂×C₁)∪…∪(A₇×B₇×C₅)

对某一时刻的输入u₁,u₂，与模糊关系R进行合成运算，得到输出C：

$\begin{array}{c}C=A^{\′}\×B^{\′}\×R\\ =\left[\begin{array}{c}\mathrm{\μ}\left(U_{11}\right)\\ \mathrm{\μ}\left(U_{12}\right)\\ .\\ .\\ .\\ \mathrm{\μ}\left(U_{17}\right)\end{array}\right]\×\left[\begin{array}{cccc}\mathrm{\μ}\left(U_{21}\right)& \mathrm{\μ}\left(U_{22}\right)& \mathrm{...}& \mathrm{\μ}\left(U_{27}\right)\end{array}\right]\×R\end{array}$

其中，U_1i,U_2i(i＝1,2,…,7)分别代表某一时刻输入量u₁,u₂在其模糊子集上的取值。最后将模糊量精确化，把经过模糊推理得到的以模糊子集表示的调节量C，根据隶属函数计算出最终输出精确值f，这里采用面积重心法进行反模糊化，重心法的计算式为：

$f=\frac{\mathrm{\Σ}\mathrm{\μ}\left(f_i\right)f_i}{\mathrm{\Σ}\mathrm{\μ}\left(f_i\right)}$

步骤6：自适应滤波算法

在滤波算法中，先将极坐标系下位置量测无偏转换到直角坐标系下，采用标准的卡尔曼滤波算法进行跟踪滤波，再利用状态滤波值对带径向加速度和径向速度的量测方程进行扩展卡尔曼滤波，并将模糊推理系统得到的f对状态向量预测协方差进行实时调整，其具体流程可以表示为：

步骤1：时间更新滤波估计

1)状态向量一步预测为

X(k)＝Φ(k|k-1)X(k-1)+W(k-1)

2)状态向量预测协方差为

P(k|k-1)＝f×Φ(k|k-1)P(k-1|k-1)+Q(k-1)

步骤2：位置量测更新滤波估计

1)量测向量一步预测为

Z₁(k|k-1)＝H₁(k)X(k|k-1)

式中，

2)新息为

In(k)＝Z₁(k)-Z₁(k|k-1)

3)新息协方差为

S(k)＝H₁(k)P(k|k-1)H₁′(k)+R₁(k)

4)卡尔曼增益为

K(k)＝P(k|k-1)H₁(k)′S-¹(k)

5)状态向量一步更新为

X(k|k)＝X(k|k-1)+K(k)In(k)

6)状态向量协方差更新为

P(k|k)＝P(k|k-1)-K(k)S(k)K′(k)

步骤3：径向速度和径向加速度量测更新滤波估计

1)量测向量一步预测为

Z₂(k|k-1)≈H₂(k)X(k-1)

2)量测雅克比矩阵为

$H_2\left(k\right)=\left[\begin{array}{cccccc}{S}_1\left(k\right)& S_2\left(k\right)& 0& S_3\left(k\right)& S_4\left(k\right)& 0\\ S_5\left(k\right)& 0& S_6\left(k\right)& S_7\left(k\right)& 0& S_8\left(k\right)\end{array}\right]$

3)新息为

In(k)＝Z₂(k)-Z₂(k|k-1)

4)新息协方差为

S₂(k)＝H₂(k)P₂(k|k-1)H₂′(k)

5)卡尔曼增益为

K₂(k)＝P(k|k-1)H₂(k)′(k)S^-1(k)

6)状态向量一步更新为

X₂(k|k)＝X₂(k|k-1)+K₂(k)In₂(k)

7)状态向量协方差更新为

P₂(k|k)＝P(k|k-1)-K₂(K)S₂(k)K′₂(k)

步骤4：最终滤波估计

1)状态向量更新为

X(k|k)＝X₂(k|k)

2)状态向量协方差更新为

P(k|k)＝P₂(k|k)

当输出量f小于1时，说明目标处于弱机动状态，这时令f为1，机动模型退化为Singer模型。

为了比较对现有技术的基于Singer模型的机动目标跟踪算法在相同条件下进行了仿真，附图7是基于传统Singer模型和基于本发明算法的目标径向距离误差曲线对比图，附图8是基于传统Singer模型和基于本发明算法的目标X轴位置误差曲线，附图9是基于传统Singer模型和基于本发明算法的目标Y轴位置误差曲线，附图10是基于传统Singer模型和基于本发明算法的目标X轴速度误差曲线，附图11是基于传统Singer模型和基于本发明算法的目标Y轴速度误差曲线，附图12是基于传统Singer模型和基于本发明算法的目标X轴加速度误差曲线，附图13是基于传统Singer模型和基于本发明算法的目标Y轴加速度误差曲线。可见，本发明和现有技术相比，具有过程稳定，收敛速度快，跟踪精度高的优势，同时滤波后X方向的速度和加速度估计精度也有了很大改善。

Claims (3)

1.一种基于模糊推理系统的Singer模型改进算法，其特征在于包括以下技术措施：

(1)信号处理阶段通过Winger-Hough变换得到每个时刻的目标径向速度和径向加速度信息，并送到雷达数据处理计算机；

(2)在数据处理阶段将得到的径向速度和径向加速度的时间变化率作为模糊推理系统输入，经过模糊推理后输出调整因子f，用来自动调节Singer模型滤波器的预测协方差；

(3)将极坐标系下位置量测无偏转换到直角坐标系下，采用标准的卡尔曼滤波算法进行跟踪滤波，再利用状态滤波值对带径向加速度和径向速度的量测方程进行扩展卡尔曼滤波，并将模糊推理系统得到的调整因子f对状态向量预测协方差进行实时调整，以实现机动目标的精确跟踪。

2.根据权利要求1所述的基于模糊推理系统的Singer模型改进算法，其特征在于：在数据处理阶段将得到的径向速度和径向加速度的时间变化率作为模糊推理系统输入，经过模糊推理后输出调整因子f，用来自动调节滤波器的预测协方差：

设F(x)＝[u₁,u₂]，u₁,u₂分别表示目标的径向速度和径向加速度信息，由于量测向量是二维的，u₁,u₂分别定义为：

$u_1=\left|{\mathrm{\Δv}}_{rk}^m\right|=|v_{rk}^m-v_{rk-1}^m|,u_2=\left|{\mathrm{\Δa}}_{rk}^m\right|=|a_{rk}^m-a_{rk-1}^m|$

其中，和分别代表k时刻信号处理阶段得到的径向速度和径向加速度值，u₁,u₂分别代表k时刻目标径向速度和径向加速度的变化率；将u₁,u₂划分为七个子集，分别为[A₁,A₂,A₃,A₄,A₅,A₆,A₇]、[B₁,B₂,B₃,B₄,B₅,B₆,B₇]，其隶属度函数可以通过模糊统计方法来确定，本专利选用常用论域[-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6]对两个输入变量u₁,u₂进行模糊化处理，定义上有7个语言变量值分别为NB、NM、NS、0、PS、PM和PB(负大、负中、负小、零、正小、正中和正大)，公式如下：

$U=\frac{12}{b-a}\[u-\frac{a+b}{2}\]$

其中，a,b代表输入变量的取值范围，将区间为[a,b]的输入量u₁,u₂转化为[-6,6]区间的变化量U₁,U₂；u₁,u₂的隶属度函数μ(U₁),μ(U₂)采用三角函数，该函数形状简单、计算量小，当输入量变化时具有较好的灵活性；假设u₁,u₂的取值范围分别为[0,60]和[0,90]，模糊子集的计算公式为：

$A={\∫}_{U_1}\mathrm{\μ}\left(U_1\right){U_1}^{-1}{\mathrm{dU}}_1$

模糊推理系统输出f是一个调整因子，可以根据目标的机动情况自动调节滤波器的预测协方差；根据经验及试验统计分析，将f的取值范围定义在区间[0,8]，将其划分为五个模糊子集[C₁,C₂,C₃,C₄,C₅]，定义上有其个语言变量值分别为为NM、NS、0、PS、和PM，隶属度函数也采用三角形函数；

当所有输入特征量处于低值时，认为目标处于弱机动目标状态，f取较小的值；当所有特征量处于极大值时，则认为目标处于强机动目标状态，f取较大的值，以保证跟踪精度，计算输入量和输出量的模糊关系R，公式如下：

R＝(A₁×B₁×C₁)∪(A₁×B₂×C₁)∪…∪(A₇×B₇×C₅)

对某一时刻的输入u₁,u₂，与模糊关系R进行合成运算，得到输出C：

$\begin{array}{c}C=A^{\′}\×B^{\′}\×R\\ =\left[\begin{array}{c}\mathrm{\μ}\left(U_{11}\right)\\ \mathrm{\μ}\left(U_{12}\right)\\ .\\ .\\ .\\ \mathrm{\μ}\left(U_{17}\right)\end{array}\right]\×\[\begin{array}{cccc}\mathrm{\μ}\left(U_{21}\right)& \mathrm{\μ}\left(U_{22}\right)& \mathrm{...}& \mathrm{\μ}\left(U_{27}\right)\end{array}\]\×R\end{array}$

其中，U_1i,U_2i(i＝1,2,…,7)分别代表某一时刻输入量u₁,u₂在其模糊子集上的取值；最后将模糊量精确化，把经过模糊推理得到的以模糊子集表示的调节量C，根据隶属函数计算出最终输出精确值f，这里采用面积重心法进行反模糊化，重心法的计算式为：

$f=\frac{\Σ\mathrm{\μ}\left(f_i\right)f_i}{\Σ\mathrm{\μ}\left(f_i\right)}.$

3.根据权利要求1所述的基于模糊推理系统的Singer模型改进算法，其特征在于：将极坐标系下位置量测无偏转换到直角坐标系下，采用标准的卡尔曼滤波算法进行跟踪滤波，再利用状态滤波值对带径向加速度和径向速度的量测方程进行扩展卡尔曼滤波，并将模糊推理系统得到的调整因子f对状态向量预测协方差进行实时调整，以实现机动目标的精确跟踪：

步骤1：时间更新滤波估计

1)状态向量一步预测为

X(k)＝Φ(k|k-1)X(k-1)+W(k-1)

2)状态向量预测协方差为

P(k|k-1)＝f×Φ(k|k-1)P(k-1|k-1)+Q(k-1)

步骤2：位置量测更新滤波估计

1)量测向量一步预测为

Z₁(k|k-1)＝H₁(k)X(k|k-1)

式中，

2)新息为

In(k)＝Z₁(k)-Z₁(k|k-1)

3)新息协方差为

S(k)＝H₁(k)P(k|k-1)H′₁(k)+R₁(k)

4)卡尔曼增益为

K(k)＝P(k|k-1)H₁(k)′S^-1(k)

5)状态向量一步更新为

X(k|k)＝X(k|k-1)+K(k)In(k)

6)状态向量协方差更新为

P(k|k)＝P(k|k-1)-K(k)S(k)K′(k)

步骤3：径向速度和径向加速度量测更新滤波估计

1)量测向量一步预测为

Z₂(k|k-1)≈H₂(k)X(k-1)

2)量测雅克比矩阵为

$H_2\left(k\right)=\left[\begin{array}{cccccc}{S}_1\left(k\right)& S_2\left(k\right)& 0& S_3\left(k\right)& S_4\left(k\right)& 0\\ S_5\left(k\right)& 0& S_6\left(k\right)& S_7\left(k\right)& 0& S_8\left(k\right)\end{array}\right]$

3)新息为

In(k)＝Z₂(k)-Z₂(k|k-1)

4)新息协方差为

S₂(k)＝H₂(k)P₂(k|k-1)H′₂(k)

5)卡尔曼增益为

K₂(k)＝P(k|k-1)H₂(k)′(k)S^-1(k)

6)状态向量一步更新为

X₂(k|k)＝X₂(k|k-1)+K₂(k)In₂(k)

7)状态向量协方差更新为

P₂(k|k)＝P(k|k-1)-K₂(K)S₂(k)K′₂(k)

步骤4：最终滤波估计

1)状态向量更新为

X(k|k)＝X₂(k|k)

2)状态向量协方差更新为

P(k|k)＝P₂(k|k)

当输出量f小于1时，说明目标处于弱机动状态，这时令f为1，机动模型退化为Singer模型。