[发明专利]一种多重测量丢失情形下时变目标跟踪系统的量化滤波方法

权利要求书

1.一种多重测量丢失情形下时变目标跟踪系统的量化滤波方法，其特征在于，该方法的具体步骤为：

步骤一、建立多重测量丢失情形下时变目标跟踪系统的动态模型：

在跟踪目标的过程中，采用具有多个传感器的雷达系统收集测量信息，建立具有多重测量丢失和信号量化情形下、受乘性噪声及加性噪声影响的时变目标跟踪系统模型：

xk+1=f(xk)+Σi=1n1αi,kAi,kxk+ωkyk=ΞkCkxk+Σi=1m1βi,kCi,kxk+vk---(1)]]>

式（1）中：

f(x_k)=Φ_kx_k+G(h(x_k)+H)，

h(xk)=-gρ(x2,k)2βx·1,k2+x·2,k2x·1,kx·2,k,]]>

ρ(x2,k)=θ1e-θ2x2,k]]>

Φk=1T000100001T0001,]]>

G=T220T00T220T,]]>

H=0-g]]>

Ai,k=0.12sin(k)0000-0.0200000.1sin(2k)00000.15]]>

Ck=10000010,]]>

Ci,k=0.15000000.30]]>

x_k+1为k+1时刻的目标跟踪系统的状态变量，y_k为k时刻的目标跟踪系统模型的测量输出；xk=x1,kx·1,kx2,kx·2,kT]]>为k时刻的目标跟踪系统的状态变量，上标“^T”表示矩阵的转置，x_1,k为k时刻的目标横坐标，为x_1,k的导数，x_2,k为k时刻的目标纵坐标，为x_2,k的导数；f(x_k)为非线性函数；T为采样周期；g为重力加速度；β为弹道系数；ρ(·)为空气密度，当x_2,k<9144米时，θ₁=1.227,θ₂=1.093×10^-4，当x_2,k≥9144米时，θ₁=1.754,θ₂=1.49×10^-4；α_i,k，β_i,k，ω_k和ν_k都是零均值的高斯白噪声，方差分别为为1，1，Q_k=c·diag{q,q}，c=0.1m²/s³，diag{·}表示对角矩阵，

q=T33T22T22T,Rk=100·1001;]]>n₁和m₁均是正整数；

矩阵为刻画多重测量丢失的矩阵，

其中，为服从伯努利分布的随机变量，i=1,2,…,m；m表示传感器的个数；

步骤二、对步骤一的目标跟踪系统模型的测量输出进行对数型量化，获得目标跟踪系统模型的测量输出的量化数据

通过公式：

qi(ykj)=ui(j),11+δjui(j)<ykj≤11-δjui(j)0,ykj=0-qj(-ykj),ykj<0,j=1,2,...,m---(2)]]>

对目标跟踪系统模型的测量输出进行对数型量化；

式中，为第j个传感器在k时刻的测量值；

δj=1-χ(j)1+χ(j);]]>

χ^(j)是取值于(0,1)中的常数；

为量化密度，它的值取自水平集Μ_j；

Mj={±ui(j),ui(j)=(χ(j))iu0(j),i=0,±1,±2,...}∪{0},u0(j)>0---(3)]]>

步骤三、计算步骤一中的多重测量数据丢失现象矩阵的数学期望，获得测量数据丢失的概率；

步骤四、根据测量数据丢失概率及跟踪系统模型的测量输出的量化数据，对多重测量丢失情形下时变目标跟踪系统状态进行一步预测和k+1时刻的状态估计；

通过公式：

x^k+1|k=f(x^k|k)---(4)]]>

x^k+1|k+1=x^k+1|k+Gk+1(y~k+1-Ξ‾k+1Ck+1x^k+1|k)---(5)]]>

获得，k时刻的一步预测值和k+1时刻系统状态的估计

式中，为x_k在k时刻的系统的状态估计；为非线性函数在处的函数值；

G_k+1为滤波器增益矩阵；

为k+1时刻滤波器接收到的信息，

其中，q(yk+1)=q1(yk+11)q2(yk+12)...qm(yk+1m)T,]]>

y_k+1为k+1时刻的测量输出，为k+1时刻的第1个传感器的测量输出，为k+1时刻的第2个传感器的测量输出，为k+1时刻的第m个传感器的测量输出；

为k+1时刻随机变量的数学期望，为k+1时刻随机变量的数学期望，为k+1时刻随机变量的数学期望；

步骤五、利用步骤四获得的系统状态的一步预测和k+1时刻的状态估计值对步骤一中的非线性函数f(x_k)进行线性化处理及余项估计，获得系统状态的一步预测误差和k+1时刻的滤波误差；

将公式（1）减去公式（4）获得一步预测误差

x~k+1|k=f(xk)-f(x^k|k)+Σi=1n1αi,kAi,kxk+Dkωk---(6)]]>

式（6）中，为k+1时刻的一步预测误差；

对非线性函数f(x_k)在处做泰勒级数展开：

f(xk)=f(x^k|k)+Akx~k|k+o(|x~k|k|)---(7)]]>

式（7）中，为k时刻的滤波误差，

Ak=∂f(xk)∂xk|xk=x^k|k,∂f(xk)∂xk]]>为函数f(x_k)关于x_k的偏导数，

o(|x~k|k|)=BkΘkLkx~k|k]]>为泰勒级数展开的高阶项，

B_k和L_k是有界矩阵，矩阵Θ_k刻画线性化误差，它满足I为单位矩阵；

步骤六、利用一步预测误差的协方差递推方程及滤波误差的协方差递推方程，获得一步预测误差P_k+1|k和滤波误差P_k+1|k+1；

一步预测误差的协方差通过公式：

Pk+1|k=E{x~k+1|kx~k+1|kT},]]>

获得，式中，为的数学期望，

通过一步预测误差的协方差递推方程：

Pk+1|k=(Ak+BkΘkLk)Pk|k(Ak+BkΘkLk)T+Σi=1n1Ai,kE{xkxkT}Ai,kT+Qk---(8)]]>

获得一步预测误差P_k+1|k，

式（8）中，为的数学期望，

滤波误差的协方差通过公式：

Pk|k=E{x~k|kx~k|kT}]]>

获得，式中，为的数学期望，

通过滤波误差的协方差递推方程：

Pk+1|k+1=(I-Gk+1Ξ‾k+1Ck+1)Pk+1|k(I-Gk+1Ξ‾k+1Ck+1)T+E{Hk+1+Hk+1T}+Gk+1Δk+1Ξ‾k+1Ck+1E{xk+1xk+1T}Ck+1TΞ‾k+1Δk+1TGk+1T+Gk+1(I+Δk+1)(Xk+1+Πk+1+Rk+1)(I+Δk+1)TGk+1T---(9)]]>

获得滤波误差P_k+1|k+1，

式（9）中：

Πk+1=Σi=1m1Ci,k+1E{xk+1xk+1T}Ci,k+1T]]>

Hk+1=-(I-Gk+1Ξ‾k+1Ck+1)x~k+1xk+1TCk+1TΞ‾k+1Δk+1TGk+1T]]>

“ο”为Hadamard（阿达玛）乘积算子，

I为单位矩阵，

Δk+1=diagΔk+1(1)Δk+1(2)...Δk+1(m)]]>为刻画量化不确定性的对角矩阵，

为矩阵△_k+1的第1个对角分量，

为矩阵△_k+1的第2个对角分量，

为矩阵△_k+1的第m个对角分量，

|Δk+1(1)|≤δ1,|Δk+1(2)|≤δ2,|Δk+1(m)|≤δm,]]>

为的范数；

步骤七、计算滤波误差的协方差P_k+1|k+1的上界并构造滤波器增益矩阵G_k+1：

滤波误差的协方差的上界通过公式：

Σk+1|k=Ak(Σk|k-1-γ1,kLkTLk)-1AkT+γ1,kTBkBkT+Qk+Σi=1n1Ai,k[(1+ϵ1)Σk|k+(1+ϵ1-1)x^k|kx^k|kT]Ai,kT---(10)]]>

Σk+1|k+1=(1+ϵ2)(I-Gk+1Ξ‾k+1Ck+1)Σk+1|k(I-Gk+1Ξ‾k+1Ck+1)T+Gk+1{(1+ϵ2-1)tr(ΛΞ‾k+1Ck+1Φk+1|kCk+1TΞ‾k+1Λ)I+tr(Ψk+1|k)[(I-γ2,k+1ΛΛ)-1+γ2,k+1-1I]}Gk+1T---(11)]]>

获得，Σ_k+1|k+1为滤波误差P_k+1|k+1的上界矩阵，

通过公式

Gk+1=Σk+1|kCk+1TΞ‾k+1{(1+ϵ2)Ξ‾k+1Ck+1Σk+1|kCk+1TΞ‾k+1+(1+ϵ2-1)tr(ΛΞ‾k+1Ck+1Φk+1|kCk+1TΞ‾k+1Λ)I+tr(Ψk+1|k)[(I-γ2,k+1ΛΛ)-1+γ2,k+1-1I]}-1---(12)]]>

获得滤波器增益矩阵G_k+1；

式（10）、（11）和（12）中，

Φk+1|k=(1+ϵ3)Σk+1|k+(1+ϵ3-1)x^k+1|kx^k+1kT]]>

Λ=diag{δ₁,δ₂,…,δ_m}

γ_1,k，γ_2,k+1，ε₁，ε₂和ε₃都是正常数，为矩阵Σ_k|k的逆，为ε₁的逆，tr(Ψ_k+1|k)为矩阵Ψ_k+1|k的迹；

步骤八、将步骤七获得的滤波器增益矩阵G_k+1带入到步骤四中的式5，获得多重测量丢失情形下时变目标跟踪系统的状态估计值，实现多重测量丢失情形下时变目标跟踪系统的量化滤波。