CN106934124A - 一种基于量测变化检测的自适应变划窗方法 - Google Patents

Abstract

本发明公开一种基于量测变化检测的自适应变划窗方法，通过系统噪声或者系统参数变化在量测上的反映以及参数估计值变化量，自适应选取量测窗长：如果前一时刻的参数估计值变化量大于阈值，则选取较大的最大窗长，否则选取较小的最大窗长；如果当前时刻的量测变化率大于前一时刻量测变化率，则增大窗长，否则减小窗长，直到小于最小窗长或者大于最大窗长；然后结合EM框架实现一种基于量测变化检测的自适应变划窗方法。以解决现有划窗法无法自适应选取窗长导致参数辨识精度低或者计算量大的问题。

Description

一种基于量测变化检测的自适应变划窗方法

技术领域

本发明属于参数辨识领域，涉及一种基于量测变化检测的自适应变划窗方法。

背景技术

在基于期望最大化(EM)算法的参数辨识框架中，算法的每一步实现均至少需要计算一个后验概率密度函数,但是如果系统量测的维数比系统状态的维数小，则下一个时刻的量测不足以构造一个有效的后验概率密度函数，相应的估计精度可能就不高。例如水平方向上的机动目标转弯运动，通常量测维数为2而状态的维数为4，因此水平方向的机动目标转弯运动也存在由量测维数小于状态维数引起的误差问题。为了解决这个问题，通常采用方法是划窗法：假设参数在划窗内是恒定不变的，选取一段时间内的量测辨识某一时刻的参数。

然而现有的划窗法选取的量测窗长是恒定的，即恒定划窗法。这种量测窗长选取方法比较容易设计，但是却无法兼顾辨识精度和计算量。如果采用这种恒定窗长法设计参数辨识算法，当选择较大的窗长时，在系统噪声较大或者系统参数发生变化时，辨识算法能得到一个精度较高、跟踪速度较快的辨识结果；但是在噪声较小或系统参数不变的时间段内(通常情况下，在整个量测时间段内，大部分时刻都是这样的)，辨识精度较高，但是时间窗内选择的量测太多而冗余，徒增计算量而对精度的提高贡献有限。当选择较小的窗长时，在系统噪声较大或者系统参数发生突变时，跟踪速度较快，但是辨识精度低；在噪声和系统参数没有发生较大变化的时刻，辨识速度快，但是精度低。

总之，现有的恒定划窗法无法实现自适应选取窗长，无法根据实际需求自适应调整窗长提高参数辨识精度或者降低算法计算量，因此如果能够设计一种自适应变划窗技术，那么就可以平衡参数辨识精度和计算量之间的矛盾。

发明内容

本发明的目的是提供一种基于量测变化检测的自适应变划窗方法，以解决现有划窗法无法自适应选取窗长导致参数辨识精度低或者计算量大的问题。

本发明所采用的技术方案是，一种基于量测变化检测的自适应变划窗方法，通过系统噪声或者系统参数变化在量测上的反映以及参数估计值变化量，自适应选取量测窗长：如果前一时刻的参数估计值变化量大于阈值，则选取较大的最大窗长，否则选取较小的最大窗长；如果当前时刻的量测变化率大于前一时刻量测变化率，则增大窗长，否则减小窗长，直到小于最小窗长或者大于最大窗长；然后结合EM框架实现一种基于量测变化检测的自适应变划窗方法。

进一步的，本发明基于的一类离散时间乘性参数线性系统如下：

x_k+1＝F(x_k)θ+S(x_k)+w_k (1-1)，

y_k＝H_kx_k+v_k (1-2)，

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

x_k+1∈Rⁿ和y_k∈R^m是系统状态和量测向量；F(x_k)∈R^n×p和S(x_k)∈Rⁿ是已知的线性函数，H_k∈R^m×n是量测矩阵；θ∈R^p是与F(x_k)乘性耦合的系统参数；w_k∈Rⁿ和v_k∈R^m分别是均值为零方差为Q和R的高斯白噪声，初始状态x₀服从均值为方差为P₀的高斯分布，k∈{1,2,…}是离散时间；

具体按照以下步骤实施：

步骤1、根据量测变化检测的自适应变划窗技术选取k时刻量测窗长：

步骤1.1、计算k-1时刻量测窗长内参数估计值变化量：

其中，为i时刻参数θ的估计值；表示k-l_k-1-1时刻到k-1时刻参数θ的估计值最大值和最小值的差值的绝对值，即参数θ的估计值在k-l_k-1-1时刻到k-1时刻的变化量；

步骤1.2、根据步骤1.1得到的参数估计值变化量，计算k时刻量测窗长l_k：

步骤a、设置最小窗长为l_min并令l_k＝l_min；

步骤b、计算k时刻和k-1时刻量测变化率和

式中，y_i为i时刻的量测，l_k和l_k-1分别表示k时刻和k-1时刻的窗长；

步骤c、根据步骤a和步骤b计算的最小窗长和k时刻和k-1时刻的量测变化率计算k时刻的窗长l_k；

步骤2、根据步骤1得到的k时刻量测窗长l_k，采用EM算法辨识系统参数θ；

步骤3、根据步骤1得到的k时刻量测窗长l_k、步骤2得到的参数辨识结果进行系统参数和状态的联合估计与辨识。

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

步骤2.1、求解离散时间乘性参数线性系统的完整数据对数似然函数；

步骤2.2、计算步骤2.1得到的完整数据对数似然函数的期望；

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

步骤b、求步骤a得到的对数条件概率密度函数关于的条件期望，其中，为在条件下的概率密度函数，为k时刻第t次迭代参数θ的估计值；

步骤2.3、将x_k，F(x_k)和S(x_k)关于状态的各分量展开；

其中，x_k为k时刻系统状态向量；F(x_k)和S(x_k)为k时刻系统已知的线性函数；

步骤2.4、根据步骤2.3得到的x_k，F(x_k)和S(x_k)的分解形式和步骤b得到的期望的表达式，计算期望，然后取期望关于参数θ的导数，并令导数为零最大化期望；

步骤2.5、根据步骤2.4中使得期望取得最大值的参数θ即可得到参数的表达式。

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

步骤3.1、计算系统状态RTS平滑和协方差获得第k次迭代的状态估计；

步骤3.2、根据步骤3.1计算的状态平滑和协方差计算参数表达式中的变量；

步骤3.3、根据步骤2.5得到的参数θ的表达式计算其估计值

步骤3.4、判断迭代是否满足迭代终止条件，所述迭代终止条件为两次迭代期望误差小于给定阈值ε：

当满足迭代终止条件后，即结束本次内层迭代过程，得到第k次迭代的参数θ的估计值然后将k取值加1，直到k等于量测长度，t值重置为1，返回步骤1进行迭代；

如果不满足迭代终止条件，将t取值加1，返回步骤3.1进行迭代。

本发明的有益效果是，本发明根据系统噪声或者系统参数变化在量测上的反映，自适应选取量测窗长，然后结合EM框架，设计一种基于量测变化检测的自适应变划窗EM算法。新技术通过量测窗长的自适应选取，解决传统恒定划窗法设置较大窗长计算量大、设置较小窗长辨识精度低的问题，平衡基于EM算法参数辨识中辨识精度和计算量的矛盾。

附图说明

图1是本发明一种基于量测变化检测的自适应变划窗方法的流程图；

图2是本发明一种基于量测变化检测的自适应变划窗方法的系统参数辨识流程图；

图3是本发明一种基于量测变化检测的自适应变划窗方法中联合估计与辨识的流程图；

图4是本发明一种基于量测变化检测的自适应变划窗方法与恒定划窗EM算法恒定转弯角速度辨识结果对比图；

图5是本发明一种基于量测变化检测的自适应变划窗方法中不同时刻窗长选取结果图；

图6本发明一种基于量测变化检测的自适应变划窗方法与恒定划窗EM算法位置均方根误差对比图；

图7本发明一种基于量测变化检测的自适应变划窗方法与恒定划窗EM算法速度均方根误差对比图。

具体实施方式

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

本发明研究的一类离散时间乘性参数线性系统如下

x_k+1＝F(x_k)θ+S(x_k)+w_k (1-1)，

y_k＝H_kx_k+v_k (1-2)，

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

x_k+1∈Rⁿ和y_k∈R^m是系统状态和量测向量；F(x_k)∈R^n×p和S(x_k)∈Rⁿ是已知的线性函数，H_k∈R^m×n是量测矩阵；θ∈R^p是与F(x_k)乘性耦合的系统参数；w_k∈Rⁿ和v_k∈R^m分别是均值为零方差为Q和R的高斯白噪声，初始状态x₀服从均值为方差为P₀的高斯分布，k∈{1,2,…}是离散时间。

不难看出，虽然系统参数耦合在系统动态方程中，但是当系统参数发生变化时，通过系统动态方程和量测方程最终会反映在量测上。同理，过程噪声和量测噪声变化时也可以通过量测反映出来。现有的恒定划窗法需要提前设置一个窗长并且在参数辨识过程中始终保持不变，这种方法的缺点是：当选择较大的窗长时，在参数保持恒定不变或者噪声变化较小的时间段存在窗内量测冗余的情况；当选择较小的窗长时，在噪声变化较大或者参数突变时刻辨识精度较低。因此设计一种基于量测变化检测的自适应变划窗方法，弥补现有恒定划窗法的缺点，平衡参数辨识精度和计算量之间的矛盾非常有意义。

本发明提出了一种基于量测变化检测的自适应变划窗方法，并将这种自适应变划窗技术应用在EM算法参数辨识中以解决上述问题，如图1，图2，图3所示，具体实现步骤如下：

步骤1、如图1流程图所示的方法设计基于量测变化检测的自适应变划窗技术选取k时刻量测窗长：

步骤1.1、计算k-1时刻量测窗长内参数估计值变化量：

其中，为i时刻参数θ的估计值，表示k-l_k-1-1时刻到k-1时刻参数θ的估计值最大值和最小值的差值的绝对值，即参数θ的估计值在k-l_k-1-1时刻到k-1时刻的变化量。

在参数辨识过程中，如果上一时刻参数估计值变化量较大(超过阈值δ)，则说明参数估计值波动较大，需要设置较大的窗长提高参数辨识精度。因此，当参数估计值变化量大于阈值δ，则k时刻最大窗长不超过l_max1，否则k时刻最大窗长不超过l_max2(l_max1>l_max2)。其中，阈值δ、最大窗长l_max1和l_max2均为提前设置的参数。

步骤1.2、计算k时刻量测窗长l_k：

步骤a、设置最小窗长为l_min并令l_k＝l_min。

步骤b、计算k时刻和k-1时刻量测变化率和

其中，y_i为i时刻的量测，l_k和l_k-1分别表示k时刻和k-1时刻的窗长；

步骤c、比较和如果即图1中的δ＝0，说明k时刻量测变化较小或者没有发生变化，则k时刻窗长为l_k；如果即图1中的δ＝0，说明k时刻量测变化较大，则l_k＝l_k+1，返回步骤b)，直到δ＝0或者l_k＝l_max。

步骤2、根据步骤1得到的k时刻量测窗长l_k，采取如图2流程图所示的方法求解参数θ表达式：

步骤2.1、构造完整数据对数似然函数

EM算法是一种迭代优化算法，每次迭代包含E步和M步，其中E步计算似然函数的期望，M步最大化期望得到参数的估计值。

通常系统状态为系统内部变化，无法直接获得，而系统量测则可以通过传感器获得，因此将系统状态作为缺失变量，系统量测作为已知数据。在EM框架下计算似然函数的期望需要计算系统状态的估计值，为了增强对待辨识参数的跟踪速度，滤波算法的执行是在某个数据段而不是整个数据区间进行的，因此，构造k-l_k到k时刻所有的状态和量测联合概率密度函数并取对数。

步骤2.2、计算期望

步骤a、根据系统动态方程和量测方程的一阶Markov特性和贝叶斯准则分解步骤2.1得到的完整数据对数似然函数，可以得到其分解之后的和式形式：k-l_k到k时刻，同一时刻量测关于状态的对数条件概率密度函数和后一时刻状态关于前一时刻状态的对数条件概率密度函数。

步骤b、求步骤a得到的对数条件概率密度函数的关于的条件期望，其中，为在条件下的概率密度函数，为k时刻第t次迭代参数θ的估计值；

由于系统参数存在于系统动态方程中，而且在M步的最大化期望是通过求期望关于参数θ的导数极大化期望实现的，因此期望计算过程中可以忽略与参数θ不相关的项：k-l_k到k时刻，同一时刻量测关于状态的对数条件概率密度函数。最后剩余k-l_k到k时刻，后一时刻状态关于前一时刻状态的对数条件概率密度函数的和式。

步骤2.3、根据状态x_k的分量分解x_k，F(x_k)和S(x_k)：

由于F(x_k)和S(x_k)是x_k的函数，即F(x_k)和S(x_k)的某些分量可能是x_k的某些分量的函数，而步骤b期望计算中积分变量为x_k或者x_k-1而不是其某些分量，为了简化期望的计算，根据状态x_k的分量分解x_k，F(x_k)和S(x_k)，将其分解为矩阵乘以状态x_k分量的和式形式。

步骤2.4、根据步骤2.3得到的x_k，F(x_k)和S(x_k)的分解形似和步骤b得到的期望的表达式，计算期望，然后取期望关于参数θ的导数，并令导数为零最大化期望。

步骤2.5、根据步骤2.4中使得期望取得最大值的参数θ即可得到参数的表达式。

步骤3、根据步骤1得到的k时刻量测窗长l_k、步骤2得到的参数辨识结果以及图3所示的算法进行系统参数和状态的联合估计与辨识：

步骤3.1、计算RTS平滑获得第k次迭代的状态估计：

设置状态和参数初始值：状态初始均值和协方差分别为和P₀，参数初始值为当k-l_k≤1时，根据状态初始值和上一次参数估计值计算离散时间乘性参数线性系统k-l到k时刻的状态平滑和协方差；当k-l_k>1时，根据上一次参数估计值结合k-l_k-1时刻的状态估计和协方差计算离散时间乘性参数线性系统k-l到k时刻的状态平滑和协方差。

步骤3.2、根据步骤3.1计算的状态平滑和协方差计算参数表达式中的变量：根据步骤3.1计算的k-l到k时刻的状态平滑和协方差和步骤2.3分解得到的常数矩阵可计算参数表达式中的变量。

步骤3.3、根据步骤2.5得到的参数θ的表达式计算其估计值根据步骤2.5得到的参数θ的表达式以及步骤3.2计算的参数表达式中的变量即可得到

步骤3.4、判断迭代是否满足迭代条件终止条件，迭代终止条件为两次迭代期望误差小于给定阈值ε。

当迭代满足迭代终止条件后，即结束本次内层迭代过程，得到第k次迭代的参数θ的估计值然后将k取值加1，直到k等于量测长度，t值重置为1，返回步骤1进行迭代；如果迭代不满足迭代终止条件，将t取值加1，返回步骤3.1进行迭代。

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

转弯机动目标跟踪模型是离散时间乘性参数线性系统，的一个特例，下面以在水平面上作转弯运动的飞机为例，验证本发明的有效性。

机动目标转弯运动的系统模型为：

y_k＝h(x_k)+v_k，

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

首先通过将状态转移矩阵中关于转弯角速度参数的非线性函数分量看作一个整体作为一个新的参数，然后转化系统模型为离散时间乘性参数线性系统，其次根据步骤2基于量测变化检测的自适应变划窗EM算法即可得到新的参数辨识结果，最后根据新的参数和转弯角速度的函数关系即可得到转弯角速度的辨识结果。

1.仿真场景：

假设机动目标在某一水平面上作转弯机动飞行，初始位置为(7000m 7000m)，初始速度为(300ms^-1300ms^-1)。假设机动目标转弯角速度真实值为Ω＝-5°s^-1。设置转弯角速度初始值为Ω₀＝-1°s^-1。

设置本发明算法的最小窗长为l_min＝2，最大窗长l_max1＝4，l_max2＝5；阈值δ＝0.4。对于上面描述的机动目标，采用本发明的基于量测变化检测的自适应变划窗EM算法(以下简称变划窗法)和基于恒定划窗法的EM的算法(以下简称恒定划窗法)比较转弯角速度参数辨识的辨识结果以及状态的均方根误差(RMSE)。基于恒定划窗法的EM的算法的窗长设置为l＝5。蒙特卡洛仿真次数为100次。

2.仿真结果：

由图4和图5可以看出，在初始时刻附近，变划窗法能够检测出转弯角速度估计值变化超过阈值δ，自适应选取最大窗长l_max2，因此辨识效果与恒定划窗法参数跟踪速度相同。除初始时刻之外，转弯角速度估计值变化较小，最大窗长为l_max1，因此变划窗法窗长均小于恒定划窗法，相应的变划窗法辨识误差大于恒定划窗法。图6和图7状态RMSE也可以说明这一点。

变划窗法位置和速度平均RMSE分别为20.3和9.1。恒定划窗法位置和速度平均RMSE分别为21.7和10.1。变划窗法平均窗长为变划窗法位置和速度平均RMSE分别比恒定划窗法减小6.5％和9.9％。变划窗法以较小的窗长实现甚至超过了恒定窗长法的转弯角速度辨识精度，降低了算法计算量，提高了算法执行效率。

由上面的仿真可以看出，本发明的算法能够以较小的窗长实现甚至超过恒定划窗法的辨识精度；本发明的算法在损失一定的辨识精度条件下能够有效降低算法的计算量。因此，本发明设计自适应变划窗技术能够平衡EM参数辨识精度和计算量之间的矛盾。

本发明方法的步骤1通过对系统量测变化进行检测，实现量测窗长的自适应变化；通过将这种自适应变划窗技术与现有的EM框架进行结合，实现了一种基于量测变化检测的自适应变划窗EM算法。现有的EM算法均是采用恒定划窗法，这种方法无法兼顾算法计算量和系统性能，而根据本发明设计的自适应变划窗技术实现的一种基于量测变化检测的自适应变划窗EM算法能够平衡参数辨识精度和计算量之间的矛盾。

Claims (4)

1.一种基于量测变化检测的自适应变划窗方法，其特征在于，通过系统噪声或者系统参数变化在量测上的反映以及参数估计值变化量，自适应选取量测窗长：如果前一时刻的参数估计值变化量大于阈值，则选取较大的最大窗长，否则选取较小的最大窗长；如果当前时刻的量测变化率大于前一时刻量测变化率，则增大窗长，否则减小窗长，直到小于最小窗长或者大于最大窗长；然后结合EM框架实现一种基于量测变化检测的自适应变划窗方法。

2.如权利要求1所述的一种基于量测变化检测的自适应变划窗方法，其特征在于，

本发明基于的一类离散时间乘性参数线性系统如下：

x_k+1＝F(x_k)θ+S(x_k)+w_k (1-1)，

y_k＝H_kx_k+v_k (1-2)，

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

x_k+1∈Rⁿ和y_k∈R^m是系统状态和量测向量；F(x_k)∈R^n×p和S(x_k)∈Rⁿ是已知的线性函数，H_k∈R^m×n是量测矩阵；θ∈R^p是与F(x_k)乘性耦合的系统参数；w_k∈Rⁿ和v_k∈R^m分别是均值为零方差为Q和R的高斯白噪声，初始状态x₀服从均值为方差为P₀的高斯分布，k∈{1,2,…}是离散时间；

具体按照以下步骤实施：

步骤1、根据量测变化检测的自适应变划窗技术选取k时刻量测窗长：

步骤1.1、计算k-1时刻量测窗长内参数估计值变化量：

${\mathrm{\Δ}}_{{\widehat{\mathrm{\θ}}}_{k-l_{k-1}-1:k-1}}=\left|\underset{k-l_{k-1}-1\≤i\≤k-1}{\max }{\widehat{\mathrm{\θ}}}_i-\underset{k-l_{k-1}-1\≤i\≤k-1}{\min }{\widehat{\mathrm{\θ}}}_i\right|---\left(2\right),$

其中，为i时刻参数θ的估计值；表示k-l_k-1-1时刻到k-1时刻参数θ的估计值最大值和最小值的差值的绝对值，即参数θ的估计值在k-l_k-1-1时刻到k-1时刻的变化量；

步骤1.2、根据步骤1.1得到的参数估计值变化量，计算k时刻量测窗长l_k：

步骤a、设置最小窗长为l_min并令l_k＝l_min；

步骤b、计算k时刻和k-1时刻量测变化率和

${\mathrm{\Δ}}_{y_{k-l_k:k}}=\frac{\left|\underset{k-l_k-1\≤i\≤k}{\max }{y}_i-\underset{k-l_k-1\≤i\≤k}{\min }{y}_i\right|}{l_k+1}---\left(3\right),$

${\mathrm{\Δ}}_{y_{k-l_{k-1}:k-1}}=\frac{\left|\underset{k-l_{k-1}-1\≤i\≤k-1}{\max }{y}_i-\underset{k-l_{k-1}-1\≤i\≤k-1}{\min }{y}_i\right|}{l_{k-1}+1}---\left(4\right),$

式中，y_i为i时刻的量测，l_k和l_k-1分别表示k时刻和k-1时刻的窗长；

步骤c、根据步骤a和步骤b计算的最小窗长和k时刻和k-1时刻的量测变化率计算k时刻的窗长l_k；

步骤2、根据步骤1得到的k时刻量测窗长l_k，采用EM算法辨识系统参数θ；

步骤3、根据步骤1得到的k时刻量测窗长l_k、步骤2得到的参数辨识结果进行系统参数和状态的联合估计与辨识。

3.如权利要求2所述的一种基于量测变化检测的自适应变划窗方法，其特征在于，所述步骤2的具体方法为：

步骤2.1、求解离散时间乘性参数线性系统的完整数据对数似然函数；

步骤2.2、计算步骤2.1得到的完整数据对数似然函数的期望；

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

步骤b、求步骤a得到的对数条件概率密度函数关于的条件期望，其中，为在条件下的概率密度函数，为k时刻第t次迭代参数θ的估计值；

步骤2.3、将x_k，F(x_k)和S(x_k)关于状态的各分量展开；

其中，x_k为k时刻系统状态向量；F(x_k)和S(x_k)为k时刻系统已知的线性函数；

步骤2.4、根据步骤2.3得到的x_k，F(x_k)和S(x_k)的分解形式和步骤b得到的期望的表达式，计算期望，然后取期望关于参数θ的导数，并令导数为零最大化期望；

步骤2.5、根据步骤2.4中使得期望取得最大值的参数θ即可得到参数的表达式。

4.如权利要求2所述的一种基于量测变化检测的自适应变划窗方法，其特征在于，所述步骤3的具体方法为：

步骤3.1、计算系统状态RTS平滑和协方差获得第k次迭代的状态估计；

步骤3.2、根据步骤3.1计算的状态平滑和协方差计算参数表达式中的变量；

步骤3.3、根据步骤2.5得到的参数θ的表达式计算其估计值

步骤3.4、判断迭代是否满足迭代终止条件，所述迭代终止条件为两次迭代期望误差小于给定阈值ε：

当满足迭代终止条件后，即结束本次内层迭代过程，得到第k次迭代的参数θ的估计值然后将k取值加1，直到k等于量测长度，t值重置为1，返回步骤1进行迭代；

如果不满足迭代终止条件，将t取值加1，返回步骤3.1进行迭代。