CN105701292A - 一种机动目标转弯角速度的解析辨识技术 - Google Patents

Abstract

本发明公开了一种机动目标转弯角速度的解析辨识技术，通过将机动目标转弯运动的系统状态转移矩阵中每一个关于转弯角速度的非线性函数分量作为一个新的参数，实现新参数与系统状态线性耦合，利用EM框架得到新参数的解析解，并反演解析辨识出转弯角速度。本发明通过将系统状态转移矩阵中每一个关于转弯角速度的非线性函数分量作为一个新的参数，实现新参数与系统状态线性耦合，以突破转弯角速度与系统状态非线性耦合的局限性，利用EM框架得到新参数的解析解，并反演解析辨识出转弯角速度，新技术的实现没有涉及任何近似方法，提高了机动目标转弯角速度的辨识精度，而解析高精度的参数辨识结果又反过来进一步改进了目标状态的估计精度。

Description

一种机动目标转弯角速度的解析辨识技术

技术领域

本发明属于机动目标跟踪领域，涉及一种机动目标转弯角速度的解析辨识技术。

背景技术

目标运动模式的不确定性是机动目标跟踪中面临的最根本的挑战之一，而转弯机动是机动目标的一种很重要的运动形式。由于目标通常具有非合作特性，即为了躲避被跟踪和锁定，目标会进行转弯机动，而对于探测方来说，目标的非合作特性会使得转弯角速度未知且时变，这就迫切要求跟踪技术在估计目标状态的同时，也能够高精度地辨识或估计出未知或时变的目标转弯角速度这一未知参数。

在目标的转弯机动模型中，转弯角速度参数被非线性耦合在线性动态方程的状态转移矩阵中，如果采用传统的状态扩维技术去估计这个参数，那么原来的目标线性模型就变成了非线性系统，自然地，由于现有的非线性动态系统滤波技术，如扩展Kalman滤波(EKF)、容积Kalman滤波(CKF)、无迹Kalman滤波(UKF)等，都是一类近似方法，因此状态和转弯角速度的估计值是近似的而非精确的，估计精度不佳。尽管传统的极大似然(ML)或期望最大化(EM)等辨识技术也可以用来解决这个问题，但转弯角速度参数的非线性耦合特性会使得性能指标的优化变成一个非线性函数最大化过程，这就不得不采用近似技术(如牛顿法)来执行性能指标的优化，进而导致参数的辨识优化结果是近似解而非解析解，这也必然引起较大的状态估计误差。尽管以H-无穷为代表的鲁棒估计技术不需要估计或辨识参数，但它要求已知参数变动所引起的干扰上界，这会导致鲁棒滤波器具有极强的保守性，即估计器的鲁棒性都是以克服这个干扰上界为依据来设计的，而实际上转弯角速度所引起的干扰只会小于等于这个干扰上界，这自然会引起目标状态估计精度变差。

发明内容

本发明的目的在于克服现有机动目标跟踪技术的转弯角速度辨识非解析和状态估计精度不佳等缺点，设计一种机动目标跟踪中转弯角速度的解析辨识方法，以提高目标跟踪的精度。

本发明所采用的技术方案是，一种机动目标转弯角速度的解析辨识技术，其特征在于，通过将机动目标转弯运动的系统状态转移矩阵中每一个关于转弯角速度的非线性函数分量作为一个新的参数，实现新参数与系统状态线性耦合，利用EM框架得到新参数的解析解，并反演解析辨识出转弯角速度。

进一步的，具体按照以下步骤进行：

步骤1、求解转弯角速度解析表达式：

1.1)对机动目标转弯运动的系统模型的动态方程进行非线性参数线性转化；其中，机动目标转弯运动的系统模型如下：即为线性系统(1)，

$x_k=f_{{\mathrm{\Ω}}_k}{x}_{k-1}+v_k---\left(1-1\right),$

y_k=h(x_k)+w_k(1-2)；

其中，公式(1-1)为系统模型的动态方程，公式(1-2)为系统模型量测方程；ξ和η分别表示X和Y方向的位置；和分别表示X和Y方向的速度；T表示采样时间间隔；v_k～N(0,Q)；w_k～N(0,R)；

1.2)将经步骤1.1非线性参数线性转化后的动态方程再转化为x_k＝Φx_k-1+Φ_k-1θ+v_k；

1.3)以步骤1.2中转化后的动态方程以及机动目标转弯运动的系统模型中的量测方程为模型，求解机动目标完整数据对数似然函数；

1.4)根据步骤1.3得到的完整数据对数似然函数，求解机动目标转弯角速度的解析表达式：

a)对完整数据的对数似然函数进行分解，并得到关于状态的对数条件概率密函数，概率密度函数均是待辨识参数θ的线性函数；

b)求步骤a得到的对数条件概率密度关于的条件期望；

c)将Φ_k-1、状态x_i和x_i-1关于状态的各分量展开；

d)求期望关于参数θ的导数，并令导数为零，得到参数θ的解析表达式：

${\widehat{\mathrm{\θ}}}_{k+1}=A^{-1}B---\left(8\right),$

则机动目标转弯角速度的解析表达式为Ω＝θ₄/θ₁；

步骤2、根据步骤1得到的转弯角速度的解析表达式进行转弯角速度辨识及状态估计：

2.1)根据机动目标转弯运动模型要求设置初始状态x₀、初始协方差p₀以及初始转弯角速度Ω₀；

2.2)计算Kalman滤波和RTS平滑获得第t次迭代的状态平滑和协方差；

2.3)根据步骤d得到的A和B的表达式，以及步骤2.2计算的状态平滑和协方差，来计算A和B；

2.4)根据步骤d得到的参数θ的解析表达式计算其估计值

2.5)判断迭代是否满足迭代终止条件，当满足迭代终止条后，即结束迭代本次内层迭代过程，即得到第k次迭代的参数θ的估计值进一步由步骤d中函数关系Ω_k＝θ₄/θ₁得到然后将k取值加1(直到k等于量测长度)，t值重置为1，再返回步骤2.2进行迭代；如果不满足迭代终止条件，将t取值加1，再返回步骤2.2进行迭代。

进一步的，步骤1.1)中非线性参数线性转化的方法具体为，

由于动态方程中的 $f_{{\mathrm{\Ω}}_k}=\left(\begin{array}{cccc}1& \frac{{\mathrm{sin\Ω}}_{k}T}{{\mathrm{\Ω}}_k}& 0& -\frac{1-{\mathrm{cos\Ω}}_{k}T}{{\mathrm{\Ω}}_k}\\ 0& {\mathrm{cos\Ω}}_{k}T& 0& -{\mathrm{sin\Ω}}_{k}T\\ 0& \frac{1-{\mathrm{cos\Ω}}_{k}T}{{\mathrm{\Ω}}_k}& 1& \frac{{\mathrm{sin\Ω}}_{k}T}{{\mathrm{\Ω}}_k}\\ 0& {\mathrm{sin\Ω}}_{k}T& 0& {\mathrm{cos\Ω}}_{k}T\end{array}\right);$

$h\left(x_k\right)=\left(\begin{array}{c}\sqrt{{\mathrm{\ξ}}_k^2+{\mathrm{\η}}_k^2}\\ {\tan }^{-1}\left(\frac{{\mathrm{\η}}_k}{{\mathrm{\ξ}}_k}\right)\end{array}\right),$

其中，T表示采样时间间隔；Ω_k表示k时刻机动目标转弯角速度；

假设θ₁＝sinΩT/Ω，θ₂＝(1-cosΩT)/Ω，θ₃＝cosΩT，θ₄＝sinΩT，

令θ＝[θ₁θ₂θ₃θ₄]^T，则Ω_k和θ的非线性函数关系为Ω_k＝θ₄/θ₁。

进一步的，步骤1.2的转化方法具体为：根据步骤1.1的线性化转化，将动态方程中含有新的参数θ的分量的状态转移矩阵即经步骤1.1转化后的与状态x_k-1按照矩阵乘积乘开，然后转化成Φ_k-1乘以θ的形式。

进一步的，步骤1.3的具体方法为：

以步骤1.2中转化后的动态方程以及机动目标转弯运动的系统模型中的量测方程(1-2)为模型，将目标的运动状态作为缺失数据，将目标镜像距和偏向角的测量数据作为已知数据，构造k-l到k时刻所有状态和量测联合概率密度函数并取对数，即得到完整数据对数似然函数。

进一步的，步骤a的具体方法为：根据步骤1.3构造的完整数据对数似然函数，利用贝叶斯准则和一阶马尔科夫链的性质对其进行分解，可以得到完整数据的对数似然函数分解形式：k-l到k时刻，同一时刻量测关于状态的对数条件概率密度和后一时刻状态关于前一时刻状态的对数条件概率密度。

进一步的，步骤b的具体方法为：对步骤a得到的完整数据对数似然函数的分解形式做关于的条件期望运算，并忽略与参数θ不相关的项，最后剩余k-l到k时刻，后一时刻状态关于前一时刻状态的对数条件概率密度的和式。

进一步的，步骤c的具体方法为：首先根据假设系统动态方程噪声服从高斯分布，将该高斯分布带入到步骤a的期望表达式中并展开，根据展开结果忽略与参数θ不相关的项，然后根据状态平滑的误差协方差的定义推导得到的两个公式如下：

$S_i=\∫x_i{x}_i^T{p}_{{\widehat{\mathrm{\θ}}}_t}\left(x_i\right|Y_{k-l}^k){\mathrm{dx}}_i=P_{i,i|k-l:k}+{\hat{x}}_{i|k-l:k}{\hat{x}}_{i|k-l:k}^T---\left(4\right),$

$T_{i-1,i}=\∫x_{i-1}{x}_i^T{p}_{{\widehat{\mathrm{\θ}}}_t}\left(x_i\right|Y_{k-l}^k){\mathrm{dx}}_{i-1}{\mathrm{dx}}_i=P_{i-1,i-1|k-l:k}{f}_{{\widehat{\mathrm{\Ω}}}_t}^T+{\hat{x}}_{i-1|k-l:k}{\hat{x}}_{i|k-l:k}^T---\left(5\right),$

同时对Φ_k-1和状态x_i和x_i-1关于状态的各分量展开，

${\mathrm{\Φ}}_{i-1}={\mathrm{\Λ}}_1{\stackrel{\·}{\mathrm{\η}}}_{i-1}+{\mathrm{\Λ}}_2{\stackrel{\·}{\mathrm{\ξ}}}_{i-1}---\left(6\right),$

$x_i={\mathrm{\Ψ}}_1{\mathrm{\ξ}}_i+{\mathrm{\Ψ}}_2{\stackrel{\·}{\mathrm{\ξ}}}_i+{\mathrm{\Ψ}}_3{\mathrm{\η}}_i+{\mathrm{\Ψ}}_4{\stackrel{\·}{\mathrm{\η}}}_i---\left(7\right),$

其中， ${\mathrm{\Λ}}_1=\left(\begin{array}{cccc}0& -1& 0& 0\\ 0& 0& 0& -1\\ 1& 0& 0& 0\\ 0& 0& 1& 0\end{array}\right);{\mathrm{\Λ}}_2=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 1\end{array}\right);$ Ψ_j,j＝1,2,3,4表示4×4单位矩阵的第j列。

进一步的，步骤d的具体方法为：根据步骤c得到的期望的表达式以及公式(4)-(7)，求期望关于参数θ的导数，并令导数为零，就能得到参数θ的解析表达式：

θ＝A^-1B(8)，

$A=\underset{i=k-l}{\overset{k}{\mathrm{\Σ}}}\left\{{\mathrm{\Λ}}_1^T{Q}^{-1}{\mathrm{\Λ}}_1{S}_{i-1}\left(4,4\right)+{\mathrm{\Λ}}_2^T{Q}^{-1}{\mathrm{\Λ}}_1{S}_{i-1}\left(2,4\right)+{\mathrm{\Λ}}_1^T{Q}^{-1}{\mathrm{\Λ}}_2{S}_{i-1}\left(2,4\right)+{\mathrm{\Λ}}_2^T{Q}^{-1}{\mathrm{\Λ}}_2{S}_{i-2}\left(2,2\right)\right\},$

$B=\underset{i=k-l}{\overset{k}{\mathrm{\Σ}}}\left\{\underset{j=1}{\overset{4}{\mathrm{\Σ}}}\left\{{\mathrm{\Λ}}_1^T{Q}^{-1}\left({\mathrm{\Ψ}}_j{T}_{i-1,i}\left(4,j\right)-{\mathrm{\Φ\Ψ}}_j{S}_{i-1}\left(4,j\right)\right)\right\}+\underset{j=1}{\overset{4}{\mathrm{\Σ}}}\left\{{\mathrm{\Λ}}_2^T{Q}^{-1}\left({\mathrm{\Ψ}}_l{T}_{i-1,i}\left(2,l\right)-{\mathrm{\Φ\Ψ}}_l{S}_{i-1}\left(2,l\right)\right)\right\}\right\}$

则机动目标转弯角速度的解析表达式为：

Ω＝θ₄/θ₁，

其中，θ₁和θ₄分别为参数θ的第1个分量和第4个分量。

进一步的，步骤2.2中的具体方式为：当k＝l+1时根据步骤2.1的参数初始值计算线性系统(1)k-l到k时刻的状态平滑和协方差；当k＞l+1时，根据上一迭代的参数估计由步骤d中的函数关系Ω_k＝θ₄/θ₁得到结合计算线性系统(1)k-l到k时刻的状态平滑和协方差。

本发明的有益效果是，本发明通过将系统状态转移矩阵中每一个关于转弯角速度的非线性函数分量作为一个新的参数，实现新参数与系统状态线性耦合，以突破转弯角速度与系统状态非线性耦合的局限性，利用EM框架得到新参数的解析解，并反演解析辨识出转弯角速度，新技术的实现没有涉及任何近似方法，提高了机动目标转弯角速度的辨识精度，而解析高精度的参数辨识结果又反过来进一步改进了目标状态的估计精度。

附图说明

图1是本发明步骤1中转弯角速度解析推导流程示意图；

图2是本发明步骤2中转弯角速度辨识及状态估计的算法实现流程示意图；

图3是本发明的解析法与扩维法转弯角速度参数辨识结果；

图4是本发明的解析法与传统EM法转弯角速度参数辨识结果；

图5是本发明的解析法与扩维法的X方向位置估计均方根误差；

图6是本发明的解析法与扩维法的X方向速度估计均方根误差；

图7是本发明的解析法与扩维法的Y方向位置估计均方根误差；

图8是本发明的解析法与扩维法的Y方向速度估计均方根误差；

图9是本发明的解析法与传统EM法的X方向位置估计均方根误差；

图10是本发明的解析法与传统EM法的X方向速度估计均方根误差；

图11是本发明的解析法与传统EM法的Y方向位置估计均方根误差；

图12是本发明的解析法与传统EM法的Y方向速度估计均方根误差；

具体实施方式

下面结合附图和具体实施方式对本发明进行详细说明。

首先，本发明研究的机动目标转弯运动的系统模型如下：即为线性系统(1)，

$x_k=f_{{\mathrm{\Ω}}_k}{x}_{k-1}+v_k---\left(1-1\right),$

y_k=h(x_k)+w_k(1-2)；

其中，公式(1-1)为系统模型的动态方程，公式(1-2)为系统模型量测方程；

$f_{{\mathrm{\Ω}}_k}=\left(\begin{array}{cccc}1& \frac{{\mathrm{sin\Ω}}_{k}T}{{\mathrm{\Ω}}_k}& 0& -\frac{1-{\mathrm{cos\Ω}}_{k}T}{{\mathrm{\Ω}}_k}\\ 0& {\mathrm{cos\Ω}}_{k}T& 0& -{\mathrm{sin\Ω}}_{k}T\\ 0& \frac{1-{\mathrm{cos\Ω}}_{k}T}{{\mathrm{\Ω}}_k}& 1& \frac{{\mathrm{sin\Ω}}_{k}T}{{\mathrm{\Ω}}_k}\\ 0& {\mathrm{sin\Ω}}_{k}T& 0& {\mathrm{cos\Ω}}_{k}T\end{array}\right);h\left(x_k\right)=\left(\begin{array}{c}\sqrt{{\mathrm{\ξ}}_k^2+{\mathrm{\η}}_k^2}\\ {\tan }^{-1}\left(\frac{{\mathrm{\η}}_k}{{\mathrm{\ξ}}_k}\right)\end{array}\right);$

ξ和η分别表示X和Y方向的位置；和分别表示X和Y方向的速度；T表示采样时间间隔；v_k～N(0,Q)；w_k～N(0,R)；Ω_k表示k时刻机动目标转弯角速度。

不难看出，Ω_k非线性耦合于状态转移的某些元素中。如果直接采用状态扩维技术将目标状态与被辨识Ω_k组成一个新的状态，那么由于中某些元素是Ω_k的非线性函数，自然地扩维之后的新状态动态方程就变成了非线性系统，而众所周知，到目前为止所有的非线性状态估计方法均为近似而非解析的技术，Ω_k的估计结果也是近似的而非解析的。同样，由于Ω_k被非线性耦合于中，传统EM技术在计算期望和极大化期望时，都不得不采用近似方法去实现。总之，现有扩维技术和EM辨识技术会导致Ω_k的估计精度损失，引起目标状态估计进度变差，迫切要求探讨新的高精度机动目标转弯角速度辨识技术。

本发明提出了一种机动目标转弯角速度的解析辨识技术，以解决上述问题，如图1和图2所示，具体实现步骤如下：

步骤1、如图1流程图所示的方法，求解转弯角速度解析表达式：

1.1)机动目标转弯运动的系统模型动态方程的非线性参数线性转化：

在目标的转弯机动模型中，转弯角速度参数被非线性耦合在线性动态方程的状态转移矩阵中，而现有转弯角速度参数的估计或辨识技术均属于近似策略，导致目标状态估计估计精度不佳。本发明克服了传统机动目标转弯角速度估计中必须使用近似方法而造成精度损失这一缺点，提出了将非线性参数线性转化的处理方法：假设系统状态转移矩阵中第i个关于系统参数Ω_k的非线性函数分量作为一个新的参数θ_i，即假设θ₁＝sinΩT/Ω，θ₂＝(1-cosΩT)/Ω，θ₃＝cosΩT，θ₄＝sinΩT，令

θ＝[θ₁θ₂θ₃θ₄]^T(2)，

则Ω_k和θ的非线性函数关系如下：

Ω_k＝θ₄/θ₁，

因此只要能解析辨识出θ，就能通过上式的反演计算来解析辨识出Ω_k，以提高目标状态估计精度。

1.2)将动态方程转化为x_k＝Φx_k-1+Φ_k-1θ+v_k：

由于根据步骤1.1)的线性转化重新整理线性系统(1)中的动态方程，

x_k＝Φx_k-1+Φ_k-1θ+v_k(3)，

其中， $\mathrm{\Φ}=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 0\end{array}\right);{\mathrm{\Φ}}_{k-1}=\left(\begin{array}{cccc}{\stackrel{\·}{\mathrm{\ξ}}}_{k-1}& -{\stackrel{\·}{\mathrm{\η}}}_{k-1}& 0& 0\\ 0& 0& {\stackrel{\·}{\mathrm{\ξ}}}_{k-1}& -{\stackrel{\·}{\mathrm{\η}}}_{k-1}\\ {\stackrel{\·}{\mathrm{\η}}}_{k-1}& {\stackrel{\·}{\mathrm{\ξ}}}_{k-1}& 0& 0\\ 0& 0& {\stackrel{\·}{\mathrm{\η}}}_{k-1}& {\stackrel{\·}{\mathrm{\ξ}}}_{k-1}\end{array}\right).$

显然，原来被非线性耦合在中的转弯角速度参数Ω_k就变成了与目标状态Φ_k-1线性耦合的参数θ，这有两点好处：转移矩阵Φ_k-1中的元素均是目标状态变量的线性函数，保障了目标状态估计的精确计算，状态估计精度无损失；Φ_k-1与新参数θ呈线性乘性关系，有利于EM框架对参数θ的解析优化辨识。

1.3)构造完整数据对数似然函数：

EM框架是一种迭代优化算法，每次迭代包含两步，E步和M步。首先构造完整数据的对数似然函数，然后E步计算对数似然函数的期望，M步极大化期望得到参数的估计值，如果两次迭代得到的期望值大于阈值要求，则返回E步继续迭代，否则迭代终止得到参数的估计值。

对于转弯机动目标跟踪问题，目标的运动状态(如位置，速度，加速度等)都是未知的，因此，采用EM算法辨识机动目标的转弯角速度时，将目标的运动状态作为缺失数据，而把目标镜像距，偏向角等可以使用雷达量测到的数据作为已知数据。由于EM框架在计算期望的过程中需要采用Kalman滤波算法获得状态的估计值，为了增强对待辨识参数的跟踪速度，滤波算法的执行是基于某个数据段而不是整个数据区间，因此构造k-l到k时刻所有状态和量测联合概率密度函数并取对数。

1.4)由对数似然函数求解转弯角速度解析表达式：

a)由于机动目标转弯运动的系统模型的动态方程和量测方程为一阶Markov链，根据贝叶斯准则和一阶马尔科夫链的性质对步骤1.3中构造的完整数据的对数似然函数进行分解，可以得到完整数据的对数似然函数分解形式：k-l到k时刻，同一时刻量测关于状态的对数条件概率密度和后一时刻状态关于前一时刻状态的对数条件概率密度，所有这些概率密度函数均是待辨识参数θ的线性函数，这是执行EM解析辨识θ的基础。

b)求步骤a)得到的对数条件概率密度关于的条件期望。由于在M步中，通过求期望关于参数θ的微分极大化期望，因此在期望的求解过程中，为了简化计算可以忽略与参数θ不相关的项，最后剩余k-l到k时刻，后一时刻状态关于前一时刻状态的对数条件概率密度的和式。

c)首先根据假设系统动态方程噪声服从高斯分布，将该高斯分布带入到步骤a)的期望表达式中并展开，根据展开结果忽略与参数θ不相关的项。在期望计算中，由于Φ_k-1包含状态的一些分量，即期望计算中被积函数包含状态某些分量，而积分变量为状态，这对于期望的计算带来了困难。为了解决这一问题，本发明发现可以根据状态平滑的误差协方差的定义来表示积分，从而简化期望的计算。根据状态平滑的误差协方差的定义推导得到的两个公式如下：

$S_i=\∫x_i{x}_i^T{p}_{{\widehat{\mathrm{\θ}}}_t}\left(x_i\right|Y_{k-l}^k){\mathrm{dx}}_i=P_{i,i|k-l:k}+{\hat{x}}_{i|k-l:k}{\hat{x}}_{i|k-l:k}^T---\left(4\right),$

$T_{i-1,i}=\∫x_{i-1}{x}_i^T{p}_{{\widehat{\mathrm{\θ}}}_t}\left(x_i\right|Y_{k-l}^k){\mathrm{dx}}_{i-1}{\mathrm{dx}}_i=P_{i-1,i-1|k-l:k}{f}_{{\widehat{\mathrm{\Ω}}}_t}^T+{\hat{x}}_{i-1|k-l:k}{\hat{x}}_{i|k-l:k}^T---\left(5\right),$

同时需要对Φ_k-1和状态x_i和x_i-1关于状态的各分量展开

${\mathrm{\Φ}}_{i-1}={\mathrm{\Λ}}_1{\stackrel{\·}{\mathrm{\η}}}_{i-1}+{\mathrm{\Λ}}_2{\stackrel{\·}{\mathrm{\ξ}}}_{i-1}---\left(6\right)$

$x_i={\mathrm{\Ψ}}_1{\mathrm{\ξ}}_i+{\mathrm{\Ψ}}_2{\stackrel{\·}{\mathrm{\ξ}}}_i+{\mathrm{\Ψ}}_3{\mathrm{\η}}_i+{\mathrm{\Ψ}}_4{\stackrel{\·}{\mathrm{\η}}}_i---\left(7\right)$

其中， ${\mathrm{\Λ}}_1=\left(\begin{array}{cccc}0& -1& 0& 0\\ 0& 0& 0& -1\\ 1& 0& 0& 0\\ 0& 0& 1& 0\end{array}\right);{\mathrm{\Λ}}_2=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 1\end{array}\right);$ Ψ_j,j＝1,2,3,4表示4×4单位矩阵的第j列。

d)由于期望值是θ的线性表达式，因此求期望关于参数θ的导数，并令导数为零，就能得到参数θ的解析表达式：

${\widehat{\mathrm{\θ}}}_{k+1}=A^{-1}B---\left(8\right)$

$A=\underset{i=k-l}{\overset{k}{\Σ}}\{{\mathrm{\Λ}}_1^T{Q}^{-1}{\mathrm{\Λ}}_1{S}_{i-1}(4,4)+{\mathrm{\Λ}}_2^T{Q}^{-1}{\mathrm{\Λ}}_1{S}_{i-1}(2,4)+{\mathrm{\Λ}}_1^T{Q}^{-1}{\mathrm{\Λ}}_2{S}_{i-1}(2,4)+{\mathrm{\Λ}}_2^T{Q}^{-1}{\mathrm{\Λ}}_2{S}_{i-1}(2,2)\}$

$B=\underset{i=k-l}{\overset{k}{\mathrm{\Σ}}}\left\{\underset{j=1}{\overset{4}{\mathrm{\Σ}}}\left\{{\mathrm{\Λ}}_1^T{Q}^{-1}\left({\mathrm{\Ψ}}_j{T}_{i-1,i}\left(4,j\right)-{\mathrm{\Φ\Ψ}}_j{S}_{i-1}\left(4,j\right)\right)\right\}+\underset{j=1}{\overset{4}{\mathrm{\Σ}}}\left\{{\mathrm{\Λ}}_2^T{Q}^{-1}\left({\mathrm{\Ψ}}_l{T}_{i-1,i}\left(2,l\right)-{\mathrm{\Φ\Ψ}}_l{S}_{i-1}\left(2,l\right)\right)\right\}\right\}$

其中， ${\mathrm{\Λ}}_1=\left(\begin{array}{cccc}0& -1& 0& 0\\ 0& 0& 0& -1\\ 1& 0& 0& 0\\ 0& 0& 1& 0\end{array}\right);{\mathrm{\Λ}}_2=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 1\end{array}\right);$ Ψ_j,j＝1,2,3,4表示4×4单位矩阵的第j列。

则机动目标转弯角速度的解析表达式为：

${\widehat{\mathrm{\Ω}}}_{k+1}={\widehat{\mathrm{\θ}}}_{k+1}^{\left(4\right)}/{\widehat{\mathrm{\θ}}}_{k+1}^{\left(1\right)}$

其中，和分别为参数的第1个分量和第4个分量。

步骤2、根据步骤1得到的转弯角速度的解析表达式以及图2所示的算法进行转弯角速度辨识及状态估计：

2.1)根据步骤1.4得到的参数θ的解析表达式，下一步骤执行迭代计算参数θ的过程。对应图2，其中t表示内层循环指标，k表示外层循环指标，状态和参数初始化：x₀,p₀,Ω₀，

根据步骤1中模型要求设置初始状态和转弯角速度参数x₀,p₀,Ω₀。

2.2)计算Kalman滤波和RTS平滑获得第k次迭代的状态估计：

当k＝l+1时根据步骤2.1的参数初始值计算线性系统(1)k-l到k时刻的状态平滑和协方差；当k＞l+1时，根据上一迭代的参数估计由步骤1中的函数关系Ω_k＝θ₄/θ₁得到结合计算线性系统(1)k-l到k时刻的状态平滑和协方差。

2.3)根据步骤1.4得到的A和B及步骤2.2计算的状态平滑和协方差计算A和B。根据步骤2.2计算的k-l到k时刻状态平滑和协方差可以得到步骤c)中k-l到k时刻对应的S_i和T_i-1,i；结合已知的常数矩阵Ψ_j,j＝1,2,3,4与Λ_j,j＝1,2，根据步骤d)中A和B的公式即可计算得到A和B的值。

2.4)根据步骤1.4中参数θ的解析表达式计算其估计值根据步骤1.4中参数θ的解析表达式：以及步骤2.3计算的A和B的值即可得到

2.5)判断迭代是否满足迭代终止条件，当满足迭代终止条后，即结束迭代本次内层迭代过程，即得到第k次迭代的参数θ的估计值进一步由步骤d)中函数关系Ω_k＝θ₄/θ₁得到然后将k取值加1(直到k等于量测长度)，t值重置为1，再返回步骤2.2进行迭代；如果不满足迭代终止条件，将t取值加1，再返回步骤2.2进行迭代。

本发明效果可通过下面的仿真实例进一步说明。

1.仿真场景及条件

假设机动目标在某一水平面作转弯机动飞行，初始位置为(7000m7000m)，初始速度为(300ms^-1300ms^-1)，转弯角速度Ω为：

系统状态方程方差为：

Q＝diag[q1Mq1M]

$M=\left[\begin{array}{cc}{T}^3/3& T^2/2\\ T^2/2& T\end{array}\right]$

量测方程方差为：

$R=diag\[\begin{array}{cc}{\mathrm{\δ}}_r^2& {\mathrm{\δ}}_{\mathrm{\θ}}^2\end{array}\]$

数据：

T＝1s

q1＝0.1s^-3

δ_r＝10m

δ_θ＝0.01rad

l＝5

状态初始值为：

x₀＝[7000m300m7000m300m]^T

协方差初始值：

P₀＝diag[1000m²100m²s^-21000m²100m²s^-2]

转弯角速度初始值为：

2.仿真内容与结果

对于上面描述的机动目标，采用本发明提出的解析法、传统的EM法和扩维法比较转弯角速度参数辨识的辨识结果以及状态的均方根误差。

由图3可见，相比与扩维法，本发明的解析法对于转弯角速度的辨识具有更高的精度，而扩维法辨识结果波动较大，误差较大。由图4可见，在采样时刻初始值附近以及转弯角速度突变点本发明的解析法具有较快的跟踪速度，说明本发明的方法具有更强的自适应性。由图5-图8可以看出，本发明的解析法对目标位置和速度的估计均方根误差均比扩维法的估计误差小。由图9-图12可以看出，相比与传统的EM法，本发明的解析法对于X方向的位置和速度和Y方向的位置的估计均方根误差在转弯角速度突变点误差较小，对于Y方向速度的估计，在50-150s本发明的解析法具有较小的均方根误差。

因此，本发明的方法具有更高的参数辨识和状态估计精度。

Claims (10)

1.一种机动目标转弯角速度的解析辨识技术，其特征在于，通过将机动目标转弯运动的系统状态转移矩阵中每一个关于转弯角速度的非线性函数分量作为一个新的参数，实现新参数与系统状态线性耦合，利用EM框架得到新参数的解析解，并反演解析辨识出转弯角速度。

2.一种机动目标转弯角速度的解析辨识技术，其特征在于，具体按照以下步骤进行：

步骤1、求解转弯角速度解析表达式：

1.1)对机动目标转弯运动的系统模型的动态方程进行非线性参数线性转化；其中，机动目标转弯运动的系统模型如下：即为线性系统(1)，

$x_k=f_{{\mathrm{\Ω}}_k}{x}_{k-1}+v_k---\left(1-1\right),$

y_k＝h(x_k)+w_k(1-2)；

其中，公式(1-1)为系统模型的动态方程，公式(1-2)为系统模型量测方程； $x={\left(\begin{array}{cccc}\mathrm{\ξ}& \stackrel{\·}{\mathrm{\ξ}}& \mathrm{\η}& \stackrel{\·}{\mathrm{\η}}\end{array}\right)}^T;$ ξ和η分别表示X和Y方向的位置；和分别表示X和Y方向的速度；T表示采样时间间隔；v_k～N(0,Q)；w_k～N(0,R)；

1.2)将经步骤1.1非线性参数线性转化后的动态方程再转化为x_k＝Φx_k-1+Φ_k-1θ+v_k；

1.3)以步骤1.2中转化后的动态方程以及机动目标转弯运动的系统模型中的量测方程为模型，求解机动目标完整数据对数似然函数；

1.4)根据步骤1.3得到的完整数据对数似然函数，求解机动目标转弯角速度的解析表达式：

a)对完整数据的对数似然函数进行分解，并得到关于状态的对数条件概率密函数，概率密度函数均是待辨识参数θ的线性函数；

b)求步骤a得到的对数条件概率密度关于的条件期望；

c)将Φ_k-1、状态x_i和x_i-1关于状态的各分量展开；

d)求期望关于参数θ的导数，并令导数为零，得到参数θ的解析表达式：

${\widehat{\mathrm{\θ}}}_{k+1}=A^{-1}B---\left(8\right),$

则机动目标转弯角速度的解析表达式为Ω＝θ₄/θ₁；

步骤2、根据步骤1得到的转弯角速度的解析表达式进行转弯角速度辨识及状态估计：

2.1)根据机动目标转弯运动模型要求设置初始状态x₀、初始协方差p₀以及初始转弯角速度Ω₀；

2.2)计算Kalman滤波和RTS平滑获得第t次迭代的状态平滑和协方差；

2.3)根据步骤d得到的A和B的表达式，以及步骤2.2计算的状态平滑和协方差，来计算A和B；

2.4)根据步骤d得到的参数θ的解析表达式计算其估计值

2.5)判断迭代是否满足迭代终止条件，当满足迭代终止条后，即结束迭代本次内层迭代过程，即得到第k次迭代的参数θ的估计值进一步由步骤d中函数关系Ω_k＝θ₄/θ₁得到然后将k取值加1(直到k等于量测长度)，t值重置为1，再返回步骤2.2进行迭代；如果不满足迭代终止条件，将t取值加1，再返回步骤2.2进行迭代。

3.如权利要求1所述的一种机动目标转弯角速度的解析辨识技术，其特征在于，所述步骤1.1)中非线性参数线性转化的方法具体为，

由于动态方程中的 $f_{{\mathrm{\Ω}}_k}=\left(\begin{array}{cccc}1& \frac{{\mathrm{sin\Ω}}_{k}T}{{\mathrm{\Ω}}_k}& 0& -\frac{1-{\mathrm{cos\Ω}}_{k}T}{{\mathrm{\Ω}}_k}\\ 0& {\mathrm{cos\Ω}}_{k}T& 0& -{\mathrm{sin\Ω}}_{k}T\\ 0& \frac{1-{\mathrm{cos\Ω}}_{k}T}{{\mathrm{\Ω}}_k}& 1& \frac{{\mathrm{sin\Ω}}_{k}T}{{\mathrm{\Ω}}_k}\\ 0& {\mathrm{sin\Ω}}_{k}T& 0& {\mathrm{cos\Ω}}_{k}T\end{array}\right);$

$h\left(x_k\right)=\left(\begin{array}{c}\sqrt{{\mathrm{\ξ}}_k^2+{\mathrm{\η}}_k^2}\\ {\tan }^{-1}\left(\frac{{\mathrm{\η}}_k}{{\mathrm{\ξ}}_k}\right)\end{array}\right),$

其中，T表示采样时间间隔；Ω_k表示k时刻机动目标转弯角速度；假设θ₁＝sinΩT/Ω，θ₂＝(1-cosΩT)/Ω，θ₃＝cosΩT，θ₄＝sinΩT，令θ＝[θ₁θ₂θ₃θ₄]^T，则Ω_k和θ的非线性函数关系为Ω_k＝θ₄/θ₁。

4.如权利要求1所述的一种机动目标转弯角速度的解析辨识技术，其特征在于，所述步骤1.2的转化方法具体为：根据步骤1.1的线性化转化，将动态方程中含有新的参数θ的分量的状态转移矩阵即经步骤1.1转化后的与状态x_k-1按照矩阵乘积乘开，然后转化成Φ_k-1乘以θ的形式。

5.如权利要求1所述的一种机动目标转弯角速度的解析辨识技术，其特征在于，所述步骤1.3的具体方法为：

以步骤1.2中转化后的动态方程以及机动目标转弯运动的系统模型中的量测方程(1-2)为模型，将目标的运动状态作为缺失数据，将目标镜像距和偏向角的测量数据作为已知数据，构造k-l到k时刻所有状态和量测联合概率密度函数并取对数，即得到完整数据对数似然函数。

6.如权利要求1所述的一种机动目标转弯角速度的解析辨识技术，其特征在于，所述步骤a的具体方法为：根据步骤1.3构造的完整数据对数似然函数，利用贝叶斯准则和一阶马尔科夫链的性质对其进行分解，可以得到完整数据的对数似然函数分解形式：k-l到k时刻，同一时刻量测关于状态的对数条件概率密度和后一时刻状态关于前一时刻状态的对数条件概率密度。

7.如权利要求6所述的一种机动目标转弯角速度的解析辨识技术，其特征在于，所述步骤b的具体方法为：对步骤a得到的完整数据对数似然函数的分解形式做关于的条件期望运算，并忽略与参数θ不相关的项，最后剩余k-l到k时刻，后一时刻状态关于前一时刻状态的对数条件概率密度的和式。

8.如权利要求7所述的一种机动目标转弯角速度的解析辨识技术，其特征在于，所述步骤c的具体方法为：首先根据假设系统动态方程噪声服从高斯分布，将该高斯分布带入到步骤a的期望表达式中并展开，根据展开结果忽略与参数θ不相关的项，然后根据状态平滑的误差协方差的定义推导得到的两个公式如下：

$S_i=\∫x_i{x}_i^T{p}_{{\widehat{\mathrm{\θ}}}_t}\left(x_i\right|Y_{k-l}^k){\mathrm{dx}}_i=P_{i,i|k-l:k}+{\hat{x}}_{i|k-l:k}{\hat{x}}_{i|k-l:k}^T---\left(4\right),$

$T_{i-1,i}=\∫x_{i-1}{x}_i^T{p}_{{\widehat{\mathrm{\θ}}}_t}\left(x_i\right|Y_{k-l}^k){\mathrm{dx}}_{i-1}{\mathrm{dx}}_i=P_{i-1,i-1|k-l:k}{f}_{{\widehat{\mathrm{\Ω}}}_t}^T+{\hat{x}}_{i-1|k-l:k}{\hat{x}}_{i|k-l:k}^T---\left(5\right),$

同时对Φ_k-1和状态x_i和x_i-1关于状态的各分量展开，

${\mathrm{\Φ}}_{i-1}={\mathrm{\Λ}}_1{\stackrel{\·}{\mathrm{\η}}}_{i-1}+{\mathrm{\Λ}}_2{\stackrel{\·}{\mathrm{\ξ}}}_{i-1}---\left(6\right),$

$x_i={\mathrm{\Ψ}}_1{\mathrm{\ξ}}_i+{\mathrm{\Ψ}}_2{\stackrel{\·}{\mathrm{\ξ}}}_i+{\mathrm{\Ψ}}_3{\mathrm{\η}}_i+{\mathrm{\Ψ}}_4{\stackrel{\·}{\mathrm{\η}}}_i---\left(7\right),$

其中， ${\mathrm{\Λ}}_1=\left(\begin{array}{cccc}0& -1& 0& 0\\ 0& 0& 0& -1\\ 1& 0& 0& 0\\ 0& 0& 1& 0\end{array}\right);$ ${\mathrm{\Λ}}_2=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 1\end{array}\right);$ Ψ_j,j＝1,2,3,4表示4×4单位矩阵的第j列。

9.如权利要求8所述的一种机动目标转弯角速度的解析辨识技术，其特征在于，所述步骤d的具体方法为：根据步骤c得到的期望的表达式以及公式(4)-(7)，求期望关于参数θ的导数，并令导数为零，就能得到参数θ的解析表达式：

θ＝A^-1B(8)，

$A=\underset{i=k-l}{\overset{k}{\mathrm{\Σ}}}\left\{{\mathrm{\Λ}}_1^T{Q}^{-1}{\mathrm{\Λ}}_1{S}_{i-1}\left(4,4\right)+{\mathrm{\Λ}}_2^T{Q}^{-1}{\mathrm{\Λ}}_1{S}_{i-1}\left(2,4\right)+{\mathrm{\Λ}}_1^T{Q}^{-1}{\mathrm{\Λ}}_2{S}_{i-1}\left(2,4\right)+{\mathrm{\Λ}}_2^T{Q}^{-1}{\mathrm{\Λ}}_2{S}_{i-1}\left(2,2\right)\right\},$

$B=\underset{i=k-l}{\overset{k}{\mathrm{\Σ}}}\left\{\underset{j=1}{\overset{4}{\mathrm{\Σ}}}\left\{{\mathrm{\Λ}}_1^T{Q}^{-1}\left({\mathrm{\Ψ}}_j{T}_{i-1,i}\left(4,j\right)-{\mathrm{\Φ\Ψ}}_j{S}_{i-1}\left(4,j\right)\right)\right\}+\underset{j=1}{\overset{4}{\mathrm{\Σ}}}\left\{{\mathrm{\Λ}}_2^T{Q}^{-1}\left({\mathrm{\Ψ}}_l{T}_{i-1,i}\left(2,l\right)-{\mathrm{\Φ\Ψ}}_l{S}_{i-1}\left(2,l\right)\right)\right\}\right\}$ 则机动目标转弯角速度的解析表达式为：

Ω＝θ₄/θ₁，

其中，θ₁和θ₄分别为参数θ的第1个分量和第4个分量。

10.如权利要求1所述的一种机动目标转弯角速度的解析辨识技术，其特征在于，所述步骤2.2中的具体方式为：当k＝l+1时根据步骤2.1的参数初始值计算线性系统(1)k-l到k时刻的状态平滑和协方差；当k＞l+1时，根据上一迭代的参数估计由步骤d中的函数关系Ω_k＝θ₄/θ₁得到结合计算线性系统(1)k-l到k时刻的状态平滑和协方差。