CN106931993A - 一种量测窗长自适应滑动的时变参数估计方法 - Google Patents

Abstract

本发明公开了一种量测窗长自适应滑动的时变参数估计方法，用自适应滑动的窗长选择更加有用的量测进行参数估计，通过传感器测量含有未知参数的系统以获得量测，再将量测传给第一级事件驱动量测选择器，其量测选择条件与现有事件驱动策略完全相同；经第一级事件驱动量测选择器选择后，将量测通过通信网络传递给第二级事件驱动量测选择器，以自适应选取靠近当前时刻的一个滑动窗长内的所有量测，第二级事件驱动的触发条件是后一时刻量测来决定前一时刻量测的取舍，依次类推，直至事件驱动条件不满足，舍掉不满足事件驱动条件时刻以前的所有量测；最后将窗长内的量测传递给估计器对参数进行估计。以解决FIR线性系统时变或突变参数的快速高精度估计与辨识问题。

Description

一种量测窗长自适应滑动的时变参数估计方法

技术领域

本发明属于参数估计或辨识技术领域，涉及一种量测窗长自适应滑动的时变参数估计方法。

背景技术

在估计与辨识技术中，量测信息越丰富，越有助于提升状态或参数的估计精度，但从算法的执行效率角度来说，在保障估计高精度或至少精度不会损失太多的前提下，如何有效快速地处理量测信息一直是工程技术关注的热点。基于事件驱动的量测选取策略为解决工程应用(特别是网络化多传感器系统)中算法快速执行问题提供了一种很好的解决方案，该方法通过在传感器测量与网络通讯之间设置传递触发或驱动条件，如图2所示，目的是选择对估计或辨识相对重要的量测，如果量测满足传递触发条件，则传感器传递量测的事件被驱动，否则传感器即使采集到量测，信息也不会被传递。事件驱动的量测选取策略实现了估计或辨识在精度与计算量之间的有效折衷，一方面被选择的相对重要的量测保障了估计或辨识精度，而事件触发又降低了需要被处理的量测信息量，降低了通信负担及算法计算负担。

本发明面向单传感器在不同采样时刻的量测时间驱动策略设计问题，主要解决的是FIR线性系统参数时变或突变情况下的快速高精度跟踪估计问题。如图1所示，现有事件驱动策略仅针对当前时刻量测进行传递事件判断，从传感器视角来看，不同时刻量测的重要性因时间驱动而不同，但从估计器或辨识角度来说，被传感器传递至估计器的量测具有同等重要性，这对于辨识常值参数是比较有效的，但在辨识时变参数时，则精度不佳。因为对于时变参数来说，越靠近当前时刻的量测越能及时准确地反映参数的变化特征，如果借鉴事件驱动的思路来设计一种量测窗长的自适应选取或滑动机制，即根据参数的当前变化特征，自适应选取靠近当前时刻的一段时间窗内的所有量测，而舍掉该时间窗以前的所有量测，这不仅可以实现对时变参数的快速估计与辨识，而且可以降低计算量。

发明内容

本发明的目的是提供一种量测窗长自适应滑动的时变参数估计方法，以解决FIR线性系统时变或突变参数的快速高精度估计与辨识问题。

本发明所采用的技术方案是，一种量测窗长自适应滑动的时变参数估计方法，用自适应滑动的窗长选择更加有用的量测进行参数估计，通过传感器测量含有未知参数的系统以获得量测，再将量测传给第一级事件驱动量测选择器，其量测选择条件与现有事件驱动策略完全相同；经第一级事件驱动量测选择器选择后，将量测通过通信网络传递给第二级事件驱动量测选择器，以自适应选取靠近当前时刻的一个滑动窗长内的所有量测，第二级事件驱动的触发条件是后一时刻量测来决定前一时刻量测的取舍，依次类推，直至事件驱动条件不满足，舍掉不满足事件驱动条件时刻以前的所有量测；最后将窗长内的量测传递给估计器对参数进行估计。

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

第一步、建立FIR线性系统如下：

其中k是当前采样时刻，y∈R^m是量测，观测矩阵H_k∈R^n×m满足θ∈R^m是待估计的未知时变参数，包括待估计参数Iθ₁和待估计参数IIθ₂，v_k∈R^m是零均值高斯白噪声，噪声的方差是且是正定矩阵；

第二步、选择由第一步中FIR线性系统产生的量测，分别确定一级与二级事件驱动量测选择器的规则和含有待优化参数的使用概率；

第三步、由事件驱动量测选择器规则计算非迭代形式的对数极大似然函数，进而得到初步参数估计表达式；

第四步、优化第三步中得到的初步参数估计解析表达式，得到精确参数估计解析表达式。

进一步的，第二步的具体方法为，设计二级事件驱动量测选择器的迭代形式规则如下：

z_k＝{y_k,γ_kz_k-1} (3)，

其中{z₁}＝{y₁}，k为当前时刻，y_k是k时刻传感器测量得到的量测值，γ_k为k时刻由使用概率为成功概率的二项分布产生的值，z_k是二级事件驱动量测选择器k时刻的输出；

使用概率如下：

其中{τ_k∈R^m|1≤k≤N}，正定矩阵两个待优化参数，通过后面的设计实现理想的估计效果，y_k是k时刻传感器测得的量测；

通过当前时刻传递过来的量测，计算出使用概率，将此概率作为二项分布成功的概率，根据此二项分布的输出为0还是1决定是否使用前一时刻量测，若为0，则此前时刻量测均舍弃，若为1则使用前一时刻量测，依次类推直到二项分布输出为0。

进一步的，第三步的具体方法为，结合事件驱动量测选择器迭代规则和FIR线性系统，写出含待估计参数的对数极大似然函数如下：

其中：

I_n为n阶单位矩阵

，n为量测的维数；

对对数极大似然函数求一阶二阶导数，由于对数极大似然函数为凸函数，令其一阶导函数为零，则可以得到初步参数估计的解析表达式如下，其中含有待定参数：

其中

为k时刻对待估计参数的估计值，u为哑元，Δ_k和τ_k为优化参数，∑为噪声方差。

进一步的，第四步的具体方法为，在得到了初步的参数估计解析表达式之后，计算参数估计期望的表达式如下：

将代入初步参数估计解析表达式(8)，得到优化后的参数估计解析表达式如下：

其中

为k时刻对待估计参数的估计值，θ_k-1为k-1时刻对待估计参数的估计值，u为哑元，Δ_k为优化参数，Σ为噪声方差。

本发明的有益效果是，通过设置第二级事件驱动来实现量测窗长的自适应滑动，相当于量测被第二次进行了取舍，算法计算量进一步降低，而在窗长内的量测可以认为是最能及时准确反映参数时变或突变特性的，如果用这些量测去估计参数，突变参数辨识精度会大大提高。

附图说明

图1是本发明一种量测窗长自适应滑动的时变参数估计方法中第一级事件驱动策略框图；

图2是本发明一种量测窗长自适应滑动的时变参数估计方法中第二级事件驱动策略框图；

图3是本发明一种量测窗长自适应滑动的时变参数估计方法中在两级事件驱动策略下，FIR参数辨识技术框图；

图4是本发明一种量测窗长自适应滑动的时变参数估计方法中待估计参数I突变时辨识性能比较；

图5是本发明一种量测窗长自适应滑动的时变参数估计方法中待估计参数II突变时辨识性能比较；

图6是本发明一种量测窗长自适应滑动的时变参数估计方法中第二级事件驱动量测选取直观图；

图7是本发明一种量测窗长自适应滑动的时变参数估计方法中参数估计框图；

图8是本发明一种量测窗长自适应滑动的时变参数估计方法中二级事件驱动量测选择器框图；

图9是本发明一种量测窗长自适应滑动的时变参数估计方法中优化参数估计表达式框图。

具体实施方式

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

本发明采用两级事件驱动策略，如图3所示。第一级事件驱动策略设置在传感器与通信带宽之间，其量测传递条件与现有事件驱动策略完全相同，第二级事件驱动策略设置在通信网络与估计器之间，经第一级传递后的量测在被用于估计参数时，第二级事件驱动的原则是越靠近当前时刻，量测的重要性越强，目的是自适应选取靠近当前时刻的一个窗长内的所有量测，第二级事件驱动的触发条件是后一时刻量测来决定前一时刻量测的取舍，依次类推，直至事件驱动条件不满足，舍掉不满足事件驱动条件时刻以前的所有量测。由于第二级事件驱动的存在，每个时刻量测窗长是根据事件驱动策略自适应选取的，从整个时间空域来说，量测窗长在事件驱动策略下自适应滑动。

本发明提供了一种量测窗长自适应滑动的时变参数估计方法，如图7参数估计框图所示，具体按照以下步骤实施：

第一步：本发明考虑的FIR线性系统如下：

其中k是当前采样时刻，y∈R^m是量测，观测矩阵H_k∈R^n×m满足θ∈R^m是待估计的未知时变参数，包括待估计参数Iθ₁和待估计参数IIθ₂，v_k∈R^m是零均值高斯白噪声，噪声的方差是且是正定矩阵；

第二步：选择由第一步中系统产生的量测，分别确定一级与二级事件驱动量测选择器的规则和含有待优化参数的使用概率。

一级事件驱动量测选择器与二级事件驱动量测选择器协同工作，如图3所示，传感器k时刻测出量测y_k，将量测传给一级事件驱动量测选择器，其k时刻输出u_k通过通信网络传给二级事件驱动量测选择器，最后估计器k时刻接收到的是二级事件驱动量测选择器k时刻的输出z_k。估计器通过z_k对参数进行估计与辨识。

我们的目的是精确并且快速的估计未知的时变参数，为了达到目的现在有多种基于事件驱动的量测选择策略，以选择相对重要的量测。这样能够保证或者至少不会损失太多估计精度。但是由于需要被处理的量测数量大幅度减少，计算效率会大大提高。

一级事件驱动量测选择器的规则如下：

u_N＝{u₁,...,u_k,...,u_N} (2)，

并且

其中y_k为k时刻一级事件驱动量测选择器的输入，u_N为时刻1到时刻N的一级事件驱动量测选择器的输出序列，u_k为k时刻的输出，λ_k是仅能决定k时刻是否传输量测y_k的决定变量。这里λ_k的取值与时刻无关，即每一个时刻量测权重都是一样的。

如图1所示，在现有事件驱动规则下，即本发明中的一级事件驱动量测选择器，被传到通信网络中的量测的权重都是一样的，意味着在估计器估计未知的时变参数时，被传送的所有量测重要性都是一样的。然而事实上量测采样时刻越靠近当前时刻，就能够更高的反映出参数特性，比如定常，时变或者突变特性。因此即使现有规则可以估计未知参数，但是尤其在估计时变参数时，其跟踪速度和估计精度很难达到实际估计精度和效率的要求。这是因为所有量测在同样的权重下，采样时间远离当前时刻的量测与靠近当前时刻的量测被选择的概率一样，但是前者也许会降低甚至污染参数估计的精度和速度。

为了克服上述困难并且实现在保持高估计精度的前提下，估计器的跟踪性能大幅度提升，在通信网络之后加上了二级事件驱动量测选择器，其规则如下：

z_k＝{y_k,γ_kz_k-1} (3)，

其中{z₁}＝{y₁}。如图2所示，y_k是k时刻传感器测量得到的量测值。γ_k为k时刻由使用概率为成功概率的二项分布产生的值，其等于0或1。用于决定k时刻的量测y_k是否被用于参数估计与辨识。z_k是基于事件驱动的量测选取器k时刻的输出，被估计器用于参数估计与辨识。

如图8所示，k时刻二级事件驱动量测选择器将前一时刻，即k-1时刻的输出z_k-1通过反馈方式和k时刻产生的决定变量γ_k结合起来，即通过γ_k的值来决定反馈回来的z_k-1是否使用，在将决定变量γ_k与z_k-1结合起来后，加上k时刻通信网络传递过来的最新量测y_k构成新的序列z_k,即为二级事件驱动量测选择器k时刻的输出，输出给估计器用于参数估计与辨识。

用于通信网络之后的二级事件驱动量测选择器加上通信网络之前的一级事件驱动量测选择器，二者的协同工作会更加智能的选择更加靠近当前时刻，更加重要的量测用于参数估计，相反的丢弃那些不重要的，可能会对参数估计值带来污染的量测信息，这样能够在保证精度的前提下，大大提高参数估计速度，在对时变参数进行估计时尤其有效。

使用概率方如下：

其中{τ_k∈R^m|1≤k≤N}，正定矩阵两个待优化参数，通过后面的设计实现理想的估计效果，y_k是k时刻传感器测得的量测。我们借鉴了独立多维高斯分布概率密度函数的模型，将分布的期望和方差设为两个待优化参数，这样可以通过调整这两个参数进一步优化估计表达式的精度。

第三步：由事件驱动量测选择器规则计算非迭代形式的对数极大似然函数，进而得到初步参数估计表达式：

结合事件驱动量测选择器迭代规则和FIR线性系统，写出含待估计参数的对数极大似然函数如下：

其中：

I_n为n阶单位矩阵

n为量测的维数。

求得对数极大似然函数之后，要确定参数估计的表达式，即确定带估计参数为何值时，对数极大似然函数取得最大值。这样问题就转化为对数极大似然函数凸优化问题。通过对对数极大似然函数求一阶二阶导数，可以发现对数极大似然函数为凸函数，令其一阶导函数为零，则可以得到参数估计的解析表达式，但其中含有待定参数。

对式(5)求一阶二阶导数如下：

其中

l_k(θ)为对数极大似然函数。u为哑元。Δ_k和τ_k为优化参数。

∑为噪声方差。

式(7)说明对数极大似然函数的二阶导数恒小于零，则对数极大似然函数为凸函数，我们令一阶导函数式(6)为0，可以得到估计参数的表达式如下：

式(7)和式(8)中

为k时刻对待估计参数的估计值。u为哑元。Δ_k和τ_k为优化参数。

∑为噪声方差。

第四步：优化第三步中得到的初步参数估计解析表达式，得到更为精确的参数估计解析表达式：

如图9所示，在得到了初步的参数估计解析表达式之后，为了优化估计效果并且提高精度，让初步参数估计为无偏估计。

式(9)中算出了参数估计期望的表达式，我们可以看出只有当时，参数估计才是无偏估计。将代入初步参数估计表达式(8)得到优化后的参数估计解析表达式如下：

其中

为k时刻对待估计参数的估计值。θ_k-1为时刻对带估计参数的估计值。u为哑元。Δ_k为优化参数。∑为噪声方差。

在本发明中事件驱动量测选择器实际上等效于每一个不同时刻都会生成长度不同的窗，在这个窗内的量测就会被选择用于状态估计与参数辨识。而窗的长度，即窗长就决定了使用量测数量的多少，在有的实际情况中，为了估计的速度，会降低量测使用数量，甚至设置最高上限。这时就对窗长的大小有一个约束，这个约束就与另一个待优化变量的约束形成一一映射的关系，即通过约束窗长的大小从而可以约束另一个待优化参数，但无法具体确定其大小。

这样就得到了估计效果更佳理想的参数估计解析表达式。

参数估计效果如图4、图5所示，可以明显看出本发明的参数估计精度和速度明显优于现有参数估计系统。

图6是每个时刻窗长的大小，可以明显看出本发明估计参数时都只用了数量很小，但很能反映参数特性的有用量测。

在现有的参数估计系统中，量测调度器都在通信网络之前，其通用性和可扩展性都受到了不小的限制。而且事件驱动策略设计都没有将不同时刻量测的重要性加以区分，每个量测取到的概率均相等。这就导致了采样时刻远离当前时刻的量测也被选择从而大幅度降低甚至污染了参数估计的性能，如跟踪性能和估计精度等。

我们设计的事件驱动策略在选择量测时，把量测的采样时间和量测自身的重要性联系了起来，可以更好的去除那些对参数估计有污染，不能反映当前参数特性的量测。所以在对时变的参数进行估计时，可以智能的选择靠近当前时刻的量测而抛弃那些无用的远离当前时刻的量测。这样一来理论上估计器的跟踪性能和估计的准确性会得到很大的提升。

Claims (5)

1.一种量测窗长自适应滑动的时变参数估计方法，其特征在于，用自适应滑动的窗长选择更加有用的量测进行参数估计，通过传感器测量含有未知参数的系统以获得量测，再将量测传给第一级事件驱动量测选择器，其量测选择条件与现有事件驱动策略完全相同；经第一级事件驱动量测选择器选择后，将量测通过通信网络传递给第二级事件驱动量测选择器，以自适应选取靠近当前时刻的一个滑动窗长内的所有量测，第二级事件驱动的触发条件是后一时刻量测来决定前一时刻量测的取舍，依次类推，直至事件驱动条件不满足，舍掉不满足事件驱动条件时刻以前的所有量测；最后将窗长内的量测传递给估计器对参数进行估计。

2.如权利要求1所述的一种量测窗长自适应滑动的时变参数估计方法，其特征在于，具体按照以下步骤实施：

第一步、建立FIR线性系统如下：

$y_k=H_k^T\mathrm{\θ}+v_k---\left(1\right),$

其中k是当前采样时刻，y∈R^m是量测，观测矩阵H_k∈R^n×m满足θ∈R^m是待估计的未知时变参数，包括待估计参数Iθ₁和待估计参数IIθ₂，v_k∈R^m是零均值高斯白噪声，噪声的方差是且是正定矩阵；

第二步、选择由第一步中FIR线性系统产生的量测，分别确定一级与二级事件驱动量测选择器的规则和含有待优化参数的使用概率；

第三步、由事件驱动量测选择器规则计算非迭代形式的对数极大似然函数，进而得到初步参数估计表达式；

第四步、优化第三步中得到的初步参数估计解析表达式，得到精确参数估计解析表达式。

3.如权利要求2所述的一种量测窗长自适应滑动的时变参数估计方法，其特征在于，所述第二步的具体方法为，设计二级事件驱动量测选择器的迭代形式规则如下：

z_k＝{y_k,γ_kz_k-1} (3)，

其中{z₁}＝{y₁}，k为当前时刻，y_k是k时刻传感器测量得到的量测值，γ_k为k时刻由使用概率为成功概率的二项分布产生的值，z_k是二级事件驱动量测选择器k时刻的输出；

使用概率如下：

$p({\mathrm{\γ}}_k=1|y_k=y)=1-\exp (-\frac{1}{2}||y_k-{\mathrm{\τ}}_k|{|}_{{\mathrm{\Δ}}_k^{-1}}^2)---\left(4\right),$

其中{τ_k∈R^m|1≤k≤N}，正定矩阵两个待优化参数，通过后面的设计实现理想的估计效果，y_k是k时刻传感器测得的量测；

通过当前时刻传递过来的量测，计算出使用概率，将此概率作为二项分布成功的概率，根据此二项分布的输出为0还是1决定是否使用前一时刻量测，若为0，则此前时刻量测均舍弃，若为1则使用前一时刻量测，依次类推直到二项分布输出为0。

4.如权利要求2所述的一种量测窗长自适应滑动的时变参数估计方法，其特征在于，所述第三步的具体方法为，结合事件驱动量测选择器迭代规则和FIR线性系统，写出含待估计参数的对数极大似然函数如下：

$l_k\left(\mathrm{\θ}\right)=\ln p_{y_k}\left(y\right)+(1-{\mathrm{\γ}}_k)\ln p\{{\mathrm{\γ}}_k=0\}+\underset{u=2}{\overset{k}{\Σ}}\underset{l=u}{\overset{k}{\Π}}{\mathrm{\γ}}_l\ln p_{y_{u-1}}\left(y\right)+\underset{u=3}{\overset{k}{\Σ}}\underset{l=u}{\overset{k}{\Π}}{\mathrm{\γ}}_l(1-{\mathrm{\γ}}_{u-1})p\{{\mathrm{\γ}}_{u-1}=0\}---\left(5\right),$

其中：

$p_{y_k}\left(y\right)=\frac{\exp (-\frac{1}{2}||y_k-H_k^T\mathrm{\θ}|{|}_{{\mathrm{\Σ}}^{-1}}^2)}{\sqrt{{\left(2\mathrm{\π}\right)}^n\det \left(\mathrm{\Σ}\right)}},$

I_n为n阶单位矩阵，n为量测的维数；

对对数极大似然函数求一阶二阶导数，由于对数极大似然函数为凸函数，令其一阶导函数为零，则可以得到初步参数估计的解析表达式如下，其中含有待定参数：

$\begin{array}{c}{\widehat{\mathrm{\θ}}}_k={\[H_k{\Σ}^{-1}{H}_k^T+(1-{\mathrm{\γ}}_k)H_k{({\mathrm{\Δ}}_k+\Σ)}^{-1}{H}_k^T+\underset{u=2}{\overset{k}{\Σ}}{m}_u+\underset{u=3}{\overset{k}{\Σ}}{q}_u\]}^{-1}\\ \×\[H_k{\Σ}^{-1}{y}_k+(1-{\mathrm{\γ}}_k)H_k{({\mathrm{\Δ}}_k+\Σ)}^{-1}{\mathrm{\τ}}_k+\underset{u=2}{\overset{k}{\Σ}}{s}_u+\underset{u=3}{\overset{k}{\Σ}}{c}_u\]\end{array}---\left(8\right),$

其中

为k时刻对待估计参数的估计值，u为哑元，Δ_k和τ_k为优化参数，∑为噪声方差。

5.如权利要求2所述的一种量测窗长自适应滑动的时变参数估计方法，其特征在于，所述第四步的具体方法为，在得到了初步的参数估计解析表达式之后，计算参数估计期望的表达式如下：

$\begin{array}{c}E\left({\widehat{\mathrm{\θ}}}_k\right)\\ =P\{{\mathrm{\γ}}_k=1\}E\left({\widehat{\mathrm{\θ}}}_k\right|{\mathrm{\γ}}_k=1)+P\{{\mathrm{\γ}}_k=0\}E\left({\widehat{\mathrm{\θ}}}_k\right|{\mathrm{\γ}}_k=0)\\ =P\{{\mathrm{\γ}}_k=1\}{\[H_k{\Σ}^{-1}{H}_k^T+\underset{u=2}{\overset{k}{\Σ}}{m}_u+\underset{u=3}{\overset{k}{\Σ}}{q}_u\]}^{-1}\×\left[\begin{array}{c}{H}_k{\Σ}^{-1}E\left(y_k\right)+\underset{u=2}{\overset{k}{\Σ}}\underset{l=u}{\overset{k}{\Π}}{\mathrm{\γ}}_l{H}_{u-1}{\Σ}^{-1}E\left(y_{u-1}\right)\\ +\underset{u=3}{\overset{k}{\Σ}}\underset{l=u}{\overset{k}{\Π}}{\mathrm{\γ}}_l(1-{\mathrm{\γ}}_{u-1})H_{u-1}{({\mathrm{\Δ}}_{u-1}+\Σ)}^{-1}E\left({\mathrm{\τ}}_{u-1}\right)\end{array}\right]\\ +P\{{\mathrm{\γ}}_k=0\}{(H_k{\Σ}^{-1}{H}_k^T+H_k{\left({\mathrm{\Δ}}_k+\Σ\right)}^{-1}{H}_k^T)}^{-1}\×\left(H_k{\Σ}^{-1}\right(E\left(y_{k-1}\right))+H_k{\left({\mathrm{\Δ}}_k+\Σ\right)}^{-1}E({\mathrm{\τ}}_{k-1}\left)\right)\end{array}---\left(9\right)$

将代入初步参数估计解析表达式(8)，得到优化后的参数估计解析表达式如下：

$\begin{array}{c}{\widehat{\mathrm{\θ}}}_k={\[H_k{\Σ}^{-1}{H}_k^T+(1-{\mathrm{\γ}}_k)H_k{({\mathrm{\Δ}}_k+\Σ)}^{-1}{H}_k^T+\underset{u=2}{\overset{k}{\Σ}}{m}_u+\underset{u=3}{\overset{k}{\Σ}}{q}_u\]}^{-1}\\ \×\[H_k{\Σ}^{-1}{y}_k+(1-{\mathrm{\γ}}_k)H_k{({\mathrm{\Δ}}_k+\Σ)}^{-1}{H}_k^T{\mathrm{\θ}}_{k-1}+\underset{u=2}{\overset{k}{\Σ}}{s}_u+\underset{u=3}{\overset{k}{\Σ}}{c}_u\]\end{array}---\left(10\right),$

其中

为k时刻对待估计参数的估计值，θ_k-1为k-1时刻对待估计参数的估计值，u为哑元，Δ_k为优化参数，∑为噪声方差。