CN1836372A - ***估计方法及程序和记录媒体、***估计装置 - Google Patents

Abstract

本发明确立能最佳地在理论上决定忘却系数的估计方法，同时研究其数值上稳定的估计算法和高速算法。首先，处理部从存储部或输入部读出或输入上限值γ_f (S101)。处理部根据式(15)决定忘却系数ρ(S103)。然后，处理部根据忘却系数ρ执行式(10)～式(13)的超高速H_∞滤波(S105)。处理部101计算式(17)(或下述式(18))的存在条件(S107)，如果其存在条件在所有的时刻得到的满足(S109)，则使γ_f减小Δγ，反复进行相同的处理(S111)。另一方面，在某一γ_f存在条件不能够得到满足时(S109)，对该γ_f加上Δγ作为γ_f的最佳值γ_f ^op输出到输出部并/或存储于存储部(S113)。

Description

***估计方法及程序和记录媒体、***估计装置

技术领域

本发明涉及***估计方法及程序和记录媒体、***估计装置，特别是使用根据H_∞评价基准开发的超高速H_∞滤波器的高速H_∞滤波器算法，同时实现状态估计的稳健化和忘却系数的最佳化的***估计方法及程序和记录媒体、***估计装置。

背景技术

通常，所谓***估计是根据输入输出数据估计***的输入输出关系的数学模型(传递函数、或脉冲响应等)的参数。作为代表性的应用例，有国际通信的回波消除器、数据通信的自动均衡器、音响***的回声消除器和声场重放以及汽车等的主动噪声控制等。详细情况参照非专利文献1.1993年电子情报通信学会的「数字信号处理手册」等。

基本原理

图8表示***估计用的结构图的例子(未知***也可以用IIR(Infinite ImpulseResponse；无限脉冲响应)滤波器表达)

这种***具备未知***1、自适应滤波器2。又，自适应滤波器2具有FIR数字式滤波器3、自适应算法4。

下面对同定未知***1的输出误差方式的一个例子进行说明。其中，u_k是未知***1的输入，d_k是作为所希望的信号的***的输出，d^_k是滤波器的输出。(还有，「^」是估计值的意思，是附在字符正上方的符号，但是由于输入不方便，只是记载于字符的右上方。下同。)

作为未知***的参数，通常使用脉冲响应，因此自适应滤波器利用自适应算法对FIR数字式滤波器3的系数进行调节，使图的评价误差e_k＝d_k-d^_k为最小。

又，向来，为了估计***的参数(状态)，广泛使用以误差协方差矩阵的更新式(Riccati微分方程式)为依据的卡尔门(Kalman)滤波器。详细情况示于非专利文献2、即S.Haykin：Adaptive filter theory，Prentice-Hall(1996)等。

下面对卡尔门滤波器的基本原理进行说明。

如下式所示，以状态空间模型表示的线性***

    x_k+1＝ρ^-1/2x_k，y_k＝H_kx_k+v_k    (1)的状态x_k的最小方差估计值x^_k|k用状态的误差协方差矩阵∑^_k|k-1如下所示得到。

${\hat{x}}_{k|k}={\hat{x}}_{k|k-1}+K_k(y_k-H_k{\hat{x}}_{k|k-1})$

${\hat{x}}_{k+1|k}={\mathrm{\ρ}}^{-\frac{1}{2}}{\hat{x}}_{k|k}---\left(2\right)$

$K_k={\widehat{\mathrm{\Σ}}}_{k|k-1}{H_k}^T{(\mathrm{\ρ}+H_k{\widehat{\mathrm{\Σ}}}_{k|k-1}{H_k}^T)}^{-1}---\left(3\right)$

${\widehat{\mathrm{\Σ}}}_{k|k}={\widehat{\mathrm{\Σ}}}_{k|k-1}-K_k{H}_k{\widehat{\mathrm{\Σ}}}_{k|k-1}$

${\widehat{\mathrm{\Σ}}}_{k+1|k}={\widehat{\mathrm{\Σ}}}_{k|k}/\mathrm{\ρ}---\left(4\right)$

但是，

${\hat{x}}_{0|-1}=0,{\widehat{\mathrm{\Σ}}}_{0|-1}={\mathrm{\ϵ}}_{0}I,{\mathrm{\ϵ}}_0>0---\left(5\right)$

其中，

x_k为状态矢量或只是状态；未知，这是估计的对象。

y_k为观测信号，是滤波器的输入，已知。

H_k为观测矩阵；已知。

v_k为观测噪声；未知。

ρ为忘却系数；通常由尝试错误决定。

K_k为滤波器增益；从矩阵∑^_k|k-1得出。

∑^_k|k对应于x^_k|k的误差的协方差矩阵；利用Riccati微分方程式得出。

∑^_k+1|k对应于x^_k+1|k的误差的协方差矩阵；利用Riccati微分方程式得出。

∑^_1|0对应于初始状态的协方差矩阵；本来是未知的，但是为了方便，使用ε₀I。

又，本发明人已经提出了利用高速H_∞滤波器的***同定算法(参照专利文献1)。这是为了***的同定新决定H_∞评价基准，开发以此为依据的超高速H_∞滤波器的高速算法，同时以该高速H_∞滤波器算法为依据的高速时变***同定方法。高速H_∞滤波器算法能够跟踪每单位时间步骤(step)在计算量O(N)急剧变化的时变***。在上限值的极限与高速卡尔门(Kalman)滤波器算法完全一致。利用这样的***同定能够实现时不变以及时变***的高速实时时间同定估计。

还有，作为***估计的领域中通常知道的方法，应该参照例如非专利文献2、3。

在回波消除器的应用例

在国际电话等长途电话线路中，由于信号的放大等理由，采用4线式线路。另一方面，用户线由于距离比较短，采用2线式线路。

图9是通信***与回波的说明图。在2线式线路与4线式线路的连接部，如图所示引入混合变压器，进行了阻抗匹配。该阻抗匹配如果是完全的，则来自说话人B的信号(声音)只到达说话人A。但是通常实现完全匹配很难，接受信号的一部分泄漏到4线式线路，发生在放大之后再度返回接收者(说话人A)的现象。这就是回波现象(echo)。回波现象随着传输距离的加长(延迟时间的加长)而影响变大，导致通话质量明显变差(在脉冲传输中，即使是近距离对通话质量变坏的影响也很大)。

图10表示回波消除器的原理图。

其中，如图所示导入回波消除器(echo canceller)，用能够直接观察的接受信号与回波逐次估计回波通道的脉冲响应，通过从实际回波减去利用其得到的拟似回波消除回波，谋求将其去除。

回波通道的脉冲响应的估计的进行要使残留回波e_k的均方根误差为最小。这时妨碍回波通道估计的要素是线路噪声和来自说话人A的信号(声音)。通常两个说话人同时开始说话时中断脉冲响应的估计。又，混合变压器的脉冲响应程度为50ms左右，因此取样周期采用125微秒时回波通道的脉冲响应的次数实际上为400左右。

非专利文献1

1993电子情报通信学会“数字信号处理手册”

非专利文献2

S.Haykin：Adaptive filter theory，Prentice-Hall(1996)

非专利文献3

B.Hassibi，A.H.Sayed，and T.Kailath：“Indenfinite-Quadratic Estimation andControl”SIAM(1996)

专利文献1

日本特开2002-135171号公报

发明内容

但是，式(1)～(5)那样的加入已有的忘却系数ρ的卡尔门(Kalman)滤波器中，忘却系数ρ的值必须根据尝试错误决定，非常麻烦。而且忘却系数ρ的值是否果真是最佳值也没有判定手段。

又，在卡尔门滤波器使用的误差协方差矩阵，本来与非零的任意矢量的2次形式通常是正(以下称为“正定”)，但是在利用计算机以单精度进行计算的情况下，其2次形式为负(以下称为“负定”)，可知数值上是不稳定的。而且，计算量为O(N²)(或者O(N³))，因此状态矢量x_k的幂N大的情况下，每1步骤的运算次数急剧增大，对实时处理不适合。

本发明鉴于上述存在问题，其目的在于确立能够在理论上最合适地决定忘却系数的估计方法，同时研究出其数值上稳定的估计算法和高速度的算法。又，本发明的目的在于提供能够适用于通信***和音响***的回波消除器、声场的重放或噪声控制等的***估计方法。

本发明为了解决上述存在问题，应用新考虑的H_∞最佳化方法，导出能够最佳地决定忘却系数的状态估计算法。而且取代应该是正定的误差的协方差矩阵，通过更新其因数矩阵，研究出数值上稳定的估计算法和高速算法。

本发明的第1种解决手段提供一种***估计方法、使计算机执行各步骤用的***估计程序、以及记录该程序的计算机可读记录媒体，

是对下式表示的状态空间模型，

x_k+1＝F_kx_k+G_kw_k

y_k＝H_kx_k+v_k

z_k＝H_kx_k

其中，x_k为状态矢量或只是状态

w_k是***噪声

v_k是观测噪声

y_k是观测信号

z_k是输出信号

F_k是***的动力学

G_k是驱动矩阵

作为评价基准，从用忘却系数ρ加权的干扰决定，而且将向滤波器误差的最大能量增益抑制得比对应于预先被提供的上限值γ_f的项小的估计算法中，同时进行状态估计的可靠性和忘却系数ρ的最佳化用的***估计方法、程序以及记录该程序的计算机可读记录媒体，包含下述步骤，即

处理部从存储部或输入部输入包含上限值γ_f、作为滤波器的输入的观测信号y_k、观测矩阵H_k的值的步骤、

处理部根据上述上限值γ_f决定关于状态空间模型的忘却系数ρ的步骤、

处理部从存储部读取包含初始值或某一时刻的观测矩阵H_k的值，用所述忘却系数ρ执行下式表示的超高速H_∞滤波的步骤、

           x^_k|k＝F_k-1x^_k-1|k-1+K_s，k(y_k-H_kF_k-1x^_k-1|k-1)

其中，

x^_k|k：用观测信号y_o～y_k的时刻k的状态x_k的估计值

F_k：***的动力学

K_s，k：滤波器增益

处理部存储关于超高速H_∞滤波器的求得的值的存储步骤、

处理部根据观测矩阵H_i与滤波器增益K_s，i，计算以所述上限值γ_f及所述忘却系数ρ为依据的存在条件的步骤、以及

处理部使γ_f小下去，反复执行所述超高速H_∞滤波的步骤，以此在各时刻在满足所述存在条件的范围将上限值设定得小，将该值存储于存储部的步骤。

又，本发明的第2种解决手段提供一种***估计装置，是对下式表示的状态空间模型，

x_k+1＝F_kx_k+G_kw_k

y_k＝H_kx_k+v_k

z_k＝H_kx_k

其中，x_k为状态矢量或只是状态

w_k是***噪声

v_k是观测噪声

y_k是观测信号

z_k是输出信号

F_k是***的动力学

G_k是驱动矩阵

作为评价基准，从用忘却系数ρ加权的干扰决定，而且将向滤波器误差的最大能量增益抑制得比对应于预先被提供的上限值γ_f的项小的估计算法中，同时进行状态估计的可靠性和忘却系数ρ的最佳化用的***估计装置，

具备执行估计算法的处理部、以及利用所述处理部进行读出和/或写入，存储与状态空间模型相关的各观测值、设定值、估计值的存储部，

所述处理部从存储部或输入部输入包含上限值γ_f、作为滤波器的输入的观测信号y_k、观测矩阵H_k的值、

所述处理部根据上述上限值γ_f决定关于状态空间模型的忘却系数ρ、

所述处理部从存储部读取包含初始值或某一时刻的观测矩阵H_k的值，用所述忘却系数ρ执行下式表示的的超高速H_∞滤波、

          x^_k|k＝F_k-1x^_k-1|k-1+K_s、k(y_k-H_kF_k-1x^_k-1|k-1)

其中，

x^_k|k：：用观测信号y_o～y_k的时刻k的状态x_k的估计值

F_k：***的动力学

K_s，k：滤波器增益

所述处理部存储关于超高速H_∞滤波器的求得的值、

所述处理部根据观测矩阵H_i与滤波器增益K_s，i，计算以所述上限值γ_f及以所述忘却系数ρ为依据的存在条件、而且

所述处理部使γ_f小下去，反复执行所述超高速H_∞滤波的步骤，以此在各时刻在满足所述存在条件的范围将上限值设定得小，将该值存储于存储部。

本发明的估计方法能够最佳地决定忘却系数，而且算法即使是单精度也能够稳定工作，因此能够以低成本实现高性能。通常，一般的民间通信设备等往往从成本和速度方面以单精度进行计算。因此本发明作为实用性的状态估计算法在各种工业领域都能够带来其效果。

附图说明

图1是本发明实施形态的硬件的结构图。

图2是H_∞滤波器的可靠性和忘却系数最佳化的流程图。

图3是图2中的H_∞滤波器(S105)的算法的流程图。

图4是定理2的平方根阵列算法的说明图。

图5是定理3的数值上稳定的高速算法的流程图。

图6表示脉冲响应{h_i}I＝0²³的值。

图7是定理3的数值上稳定的高速算法得到的脉冲响应的估计结果。

图8是用于***估计的结构图。

图9是关于通信***与回波的说明图。

图10表示回波消除器的原理图。

最佳实施方式

下面对本发明的实施形态进行说明。

1.记号的说明

首先对本发明实施形态中使用的主要记号及其已知或未知进行说明。

x_k：状态矢量或只是状态，未知，这是估计的对象。

x₀：初始状态，未知。

W_k：***噪声，未知。

v_k：观测噪声，未知。

y_k：观测信号；作为滤波器的输入，已知。

z_k：输出信号；未知。

F_k：***的动力学；已知。

G_k：驱动矩阵；执行时已知。

H_k：观测矩阵；已知。

x^_k|k：用观测信号y_o～y_k的时刻k的状态x_k的估计值；由滤波器方程式提供。

x^_k+1|k：用观测信号y_o～y_k的时刻k+1的状态x_k+1的估计值；由滤波器方程式提供。

x^_0|0：状态的初始估计值；本来未知，但是为了方便使用0。

∑^_k|k：对应于x^_k|k的误差的协方差矩阵；由Riccati微分方程式提供。

∑^_k+1|k：对应于x^_k+1|k的误差的协方差矩阵；由Riccati微分方程式提供。

∑^_1|0：对应于初始状态的协方差矩阵；本来未知，但是为了方便使用ε₀I。

K_s，k：滤波器增益；从矩阵∑^_k|k-1得到。

ρ：忘却系数；在定理1～3的情况下，如果γ决定，则根据ρ＝1-χ(γ_f)可自动决定ρ。

e_f，i：滤波器误差

R_e，k：辅助变量

还有，记号上所附加的“^”“ˇ”是估计值的意思。又，“～”、“-”、“∪”等是为了方便附加的记号。这些记号，为了输入方便，记载于文字的右上方，但是如数学式所示，与记载于文字的正上方是一样的。而且x、w、H、G、K、R、∑等是矩阵，如数学式所示用粗体字记载，但是为了输入的方便，用普通的文字记载。

2.***估计的硬件和程序

本***估计方法或***估计装置·***，可以利用使计算机执行其各步骤用的***估计程序、记录***估计程序的计算机可读记录媒体、包含***估计程序，能够下载于计算机的内部存储器的程序产品、包含该程序的服务器等的计算机等提供。

图1是有关本实施形态的硬件的结构图。

该硬件具有作为中央处理器(CPU)的处理部101、输入部102、输出部103、显示部104以及存储部105。又，处理部101、输入部102、输出部103、显示部104以及存储部105用星形连接或总线连接等适当的连接手段连接。存储部105根据需要存储***估计的「1.记号的说明」所示的已知的数据。又，未知·已知的数据和计算出的关于超高速H_∞滤波器的数据·其他数据，利用处理部101根据需要写入和/或读出。

3.能够最合适决定忘却系数的超高速H_∞滤波器

定理1

考虑下式那样的状态空间

x_k+1＝F_kx_k+G_kw_k， w_k，x_k∈R^N                                   (6)

  y_k＝H_kx_k+v_k，   y_k，v_k∈R                                    (7)

  z_k＝H_kx_k，z_k∈R，H_k∈R^1×N，k＝0，1，...，L                  (8)

对这样的状态空间模型，提出下式所示的H_∞评价基准。

满足该H_∞评价基准的状态估计值x^_k|k(或输出估计值zˇ_k|k)，利用下面的电平γ_f的超高速H_∞滤波器提供。

${\hat{x}}_{k|k}=F_{k-1}{\hat{x}}_{k-1|k-1}+K_{s,k}(y_k-H_k{F}_{k-1}{\hat{x}}_{k-1|k-1})---\left(11\right)$

$K_{s,k}={\widehat{\mathrm{\Σ}}}_{k|k-1}{H}_k^T\·{(H_k{\widehat{\mathrm{\Σ}}}_{k|k-1}{H}_k^T+\mathrm{\ρ})}^{-1}---\left(12\right)$

$\left.\begin{array}{c}{\widehat{\mathrm{\Σ}}}_{k|k}={\widehat{\mathrm{\Σ}}}_{k|k-1}-{\widehat{\mathrm{\Σ}}}_{k|k-1}{C}_k^T{R}_{e,k}^{-1}{C}_k{\widehat{\mathrm{\Σ}}}_{k|k-1}\\ {\widehat{\mathrm{\Σ}}}_{k+1|k}=\left(F_k{\widehat{\mathrm{\Σ}}}_{k|k}{F}_k^T\right)/\mathrm{\ρ}\end{array}\right\}---\left(13\right)$

其中，

$R_{e,k}=R_k+C_k{\widehat{\mathrm{\Σ}}}_{k|k-1}{C}_k^T,R_k=\left[\begin{array}{cc}\mathrm{\ρ}& 0\\ 0& -{\mathrm{\ρ\γ}}_f^2\end{array}\right],C_k=\left[\begin{array}{c}{H}_k\\ H_k\end{array}\right]---\left(14\right)$

0＜ρ＝1-χ(γ_f)≤1，γ_f＞1                                   (15)

还有，式(11)表示滤波器方程式，式(12)表示滤波器增益，式(13)表示Riccati微分方程式。

又，驱动矩阵G_k如下所述生成。

$G_k{G}_k^T=\frac{\mathrm{\χ}\left({\mathrm{\γ}}_f\right)}{\mathrm{\ρ}}{F}_k{\widehat{\mathrm{\Σ}}}_{k|k}{F}_k^T---\left(16\right)$

又，为了提高上述高速H_∞滤波器的跟踪能力，上限值γ_f尽可能设定得小，以满足下述存在条件。

${\widehat{\mathrm{\Σ}}}_{i|i}^{-1}={\widehat{\mathrm{\Σ}}}_{i|i-1}^{-1}+\frac{1-{\mathrm{\γ}}_f^{-2}}{\mathrm{\ρ}}{H}_i^T{H}_i>0,i=0,...,k\left(17\right)$

其中，χ(γ_f)是满足χ(1)＝1，χ(_∞)＝0的γ_f的单调衰减函数。

定理1的特征在于，状态估计的可靠性和忘却系数ρ的最佳化同时进行。

图2表示H_∞滤波器的可靠性与忘却系数ρ的最佳化的流程图。在这里，数据块「EXC＞0」：H_∞滤波器的存在条件

Δγ：正实数

首先，处理部101从存储部105或输入部102读出或输入上限值γ_f(S101)。在这一例子中，给予γ_f＞＞1。处理部101利用式(15)决定忘却系数ρ(S103)。其后，处理部101根据忘却系数ρ执行式(10)～式(13)的超高速H_∞滤波(S105)。处理部101对式(17)(或系数式(18))的右边(将其作为EXC)进行计算(S107)，其存在条件如果在所有的时刻得到满足(S109)，则使γ_f只小Δγ并反复进行相同的处理(S111)。另一方面，在某一γ_f存在条件不能够得到满足时(S109)，将该γ_f加Δγ的和作为γ_f的最佳值γ_f ^op输出到输出部103并且/或者存储于存储部105(S113)。还有，在本例子中是加Δγ，但是也可以加除此以外的预先设定的值。该最佳化的过程被称为γ-迭代。还有，处理部101也可以根据需要把H_∞滤波计算步骤S105和存在条件的计算步骤S107等各步骤求得的适当的中间值和最终值适当存储于存储部105，也可以从存储部105读出。

在超高速H_∞滤波器满足存在条件时，式(9)的不等式通常被满足。因此，在使(9)的分母的干扰能量受到限定的情况下，式(9)的分子的平方估计误差的总和是有限的，某一时刻以后的估计误差为0。这意味着如果能够使γ_f更小，则估计值x^_k|k能够快速跟踪状态x_k的变化。

在这里，希望能够注意定理1的超高速H_∞滤波器的算法与通常的H_∞滤波器的算法的不同。而且，在γ_f→∞时，ρ＝1，G_k＝0，定理1的H_∞滤波器的算法与卡尔门(Kalman)滤波器的算法一致。

图3表示图2中的(超高速)H_∞滤波器(S105)的算法的流程图。

超高速H_∞滤波器算法可以摘要说明如下。

〔步骤S201〕处理部101从存储部105读出递归式的初始条件，或者，从输入部102输入初始条件，如图所示那样决定。L表示预先决定的最大数据数。

〔步骤S203〕处理部101将时刻k与最大数据数L加以比较。处理部101在时刻k大于最大数据数时结束处理，如果在最大数据数以下则进入下面的步骤。(如果不要，可以取下条件文。又，也可以根据需要再度开始。)

〔步骤S205〕处理部101利用式(12)计算滤波器增益K_s，k。

〔步骤S207〕处理部101更新式(11)的超高速H_∞滤波器的滤波器方程式。

〔步骤S209〕处理部101用式(13)的Riccati微分方程式计算与误差的协方差矩阵对应的项∑^_k|k、∑^_k+1|k。

〔步骤S211〕使时刻k前进(k＝k+1)，返回步骤S203，只要有数据就继续。

还有，处理部101根据需要适当在存储器105存储在步骤H_∞滤波器计算步骤S205～S209等各步骤求得的合适的中间值和最终值、存在条件等，也可以从存储部105将其读出。

定标器存在条件

而式(17)的存在条件的判定需要O(N²)的计算量。但是如果采用下述条件，则根据计算量O(N)能够验证定理1的H_∞滤波器的存在、即式(9)。

系1.定标器存在条件

如果采用下述存在条件，则能够用计算量O(N)判定超高速H_∞滤波器的存在。

其中，

其中，K_s，i是用式(12)求出的滤波器增益。

证明

以下对系1的证明进行说明。

如果求解2×2的矩阵R_e，k的特征方程式，即

$|\mathrm{\λI}-R_{e,k}|=\left|\begin{array}{cc}\mathrm{\λ}-(\mathrm{\ρ}+H_k{\widehat{\mathrm{\Σ}}}_{k|k-1}{H}_k^T)& -H_k{\widehat{\mathrm{\Σ}}}_{k|k-1}{H}_k^T\\ -H_k{\widehat{\mathrm{\Σ}}}_{k|k-1}{H}_k^T& \mathrm{\λ}-(-{\mathrm{\ρ\γ}}_f^2+H_k{\widehat{\mathrm{\Σ}}}_{k|k-1}{H}_k^T)\end{array}\right|$

则能够得到如下的R_e，k的固有值λ_i。

其中，

如果

则矩阵R_e，k的两个固有值中的一个为正，另一个为负，矩阵R_k与R_e，k具有相同的贯量。因此，如果采用

$H_k{\widehat{\mathrm{\Σ}}}_{k|k-1}{H}_k^T=\frac{H_k{\tilde{K}}_k}{1-\frac{1-{\mathrm{\γ}}_f^{-2}}{\mathrm{\ρ}}{H}_k{\tilde{K}}_k},H_k{\tilde{K}}_k=\frac{H_k^{\mathrm{\ρ}}{K}_{s,k}}{1-{\mathrm{\γ}}_f^{-2}{H}_k{K}_{s,k}}$

则能够得到式(18)的存在条件。其中，H_kK_s，k的计算量为O(N)。

4.数值上稳定的状态估计算法

上述超高速H_∞滤波器为了更新∑^_k|k+1∈R^N×N，每单位时间步骤的计算量为O(N²)、即必须进行正比于N²的算术运算。其中，N为状态矢量x_k的幂。因此随着x_k的幂的增加，本滤波器的执行需要的计算时间急剧增加。而且误差的协方差矩阵∑^_k|k+1根据其性质，通常必须是正定的，但是在数值上也有成为负定的情况。特别是在以单精度计算的情况下，这种倾向显著。已经知道，这时滤波器不稳定。因此，为了算法的实用化和低成本化，最好是研究出在单精度(例如32位)也能够工作的状态估计算法。

在这里，着眼于

R_k＝R^1/2 _kJ₁R^1/2 _k、

R_e，k＝R^1/2 _e，kJ₁R^1/2 _e，k、

∑^_k|k-1＝∑^^1/2 _k|k-1∑^^1/2 _k|k-1

将数值上稳定化的定理1的H_∞滤波器(平方根阵列算法)示于定理2。但是，在这里为了简单，以F_k＝I，但是即使F_k≠I的情况下也同样能够求得。下面表示出实现数值上稳定的状态估计算法用的超高速H_∞滤波器。

定理2

${\hat{x}}_{k|k}={\hat{x}}_{k-1|k-1}+K_{s,k}(y_k-H_k{\hat{x}}_{k-1|k-1})---\left(20\right)$

$K_{e,k}=K_k(:,1)/R_{e,k}\left(\mathrm{1,1}\right),K_k={\mathrm{\ρ}}^{\frac{1}{2}}\left({\mathrm{\ρ}}^{-\frac{1}{2}}{K}_k{R}_{e,k}^{-\frac{1}{2}}{J}_1^{-1}\right)J_1{R}_{e,k}^{\frac{1}{2}}---\left(21\right)$

其中，

$R_k=R_k^{\frac{1}{2}}{J}_1{R}_k^{\frac{1}{2}},R_k^{\frac{1}{2}}=\left[\begin{array}{cc}{\mathrm{\ρ}}^{\frac{1}{2}}& 0\\ 0& {\mathrm{\ρ}}^{\frac{1}{2}}\mathrm{\γf}\end{array}\right],J_1=\left[\begin{array}{cc}1& 0\\ 0& -1\end{array}\right],{\widehat{\mathrm{\Σ}}}_{k|k-1}={\widehat{\mathrm{\Σ}}}_{k|k-1}^{\frac{1}{2}}{\widehat{\mathrm{\Σ}}}_{k|k-1}^{\frac{1}{2}}$

Θ(k)为J-酉矩阵，即满足Θ(k)JΘ(k)^T＝J，J＝(J₁_I)、I是单位矩阵。又，K_k(：，1)表示矩阵K_k的第1列的列矢量。

还有，在式(21)、(22)中，J₁ ^-1和J₁可删除。

图4是定理2的平方根阵列算法的说明图。该计算的算法可以在图2所示的定理1的流程图中的H_∞滤波器的计算(S105)中使用。

本估计算法用以J-酉变换为依据的更新式求其因数矩阵∑^^1/2 _k|k-1∈R^N×N(∑^_k|k-1的平方根矩阵)，代替用Riccati型的更新式求∑^_k|k-1。从这时产生的1-1数据块(block)矩阵和2-1数据块矩阵如图所示求滤波器增益K_s，k。因此，∑^_k|k-1＝∑^^1/2 _k|k-1∑^^1/2 _k|k-1＞0，∑^_k|k-1的正定性得到保证，在数值上能够稳定。还有，定理2的H_∞滤波器的每单位步骤的计算量仍然保持O(N²)不变。

还有，在图4中，J₁ ^-1可删除。

首先，处理部101从存储部105读出或从内部存储器等得到式(22)左边的矩阵的各要素中包含的项，执行J-酉变换(S301)。处理部101根据式(21)从求得的式(22)的右边的矩阵的要素计算***增益K_k、K_s，k(S303、S305)。处理部101根据式(20)计算状态估计值x^_k|k(S307)。

5.状态估计用的数值上稳定的高速算法

如上所述，定理2的H_∞滤波器的每单位步骤的计算量仍然是O(N²)。因此作为计算量的对策，在 H _k ＝ H _k+1 Ψ、 H _k ＝〔u(k)、...、u(o)、o、...、o〕时，考虑利用 x _k＝〔x^T _K、o^T〕^T的一步骤预测的误差的协方差矩阵∑_k+1|k满足下式(24)，即

${\underset{\‾}{\mathrm{\Σ}}}_{k+1|k}-\mathrm{\Ψ}{\underset{\‾}{\mathrm{\Σ}}}_{k|k-1}{\mathrm{\Ψ}}^T=-{\underset{\‾}{L}}_k{R}_{r,k}^{-1}{\underset{\‾}{L}}_k^T,{\underset{\‾}{L}}_k=\left[\begin{array}{c}{\tilde{L}}_k\\ 0\end{array}\right]---\left(24\right)$

更新幂小的 L _k (即L~k)，而不更新 ∑ _k+1|k。在这里，如果注意到 $R_{r,k}=R_{r,k}^{\frac{1}{2}}{\mathrm{SR}}_{r,k}^{\frac{1}{2}},$ 则能够得到如下所述的定理3。

定理3

${\hat{x}}_{k|k}={\hat{x}}_{k-1|k-1}+{\stackrel{\‾}{K}}_{s,k}(y_k-H_k{\hat{x}}_{k-1|k-1})---\left(61\right)$

${\stackrel{\‾}{K}}_{s,k}={\stackrel{\‾}{K}}_k(:,1)/R_{e,k}\left(\mathrm{1,1}\right),{\stackrel{\‾}{K}}_k={\mathrm{\ρ}}^{\frac{1}{2}}\left({\stackrel{\‾}{K}}_k{R}_{e,k}^{-\frac{1}{2}}\right)R_{e,k}^{\frac{1}{2}}---\left(62\right)$

在这里，Θ(k)为任意J-酉矩阵，Cˇ_k＝Cˇ_k+1ψ成立。其中

$R_k=R_k^{\frac{1}{2}}{J}_1{R}_k^{\frac{1}{2}},R_k^{\frac{1}{2}}=\left[\begin{array}{cc}{\mathrm{\ρ}}^{\frac{1}{2}}& 0\\ 0& {\mathrm{\ρ}}^{\frac{1}{2}}\mathrm{\γf}\end{array}\right],J_1=\left[\begin{array}{cc}1& 0\\ 0& -1\end{array}\right],{\widehat{\mathrm{\Σ}}}_{k|k-1}={\widehat{\mathrm{\Σ}}}_{k|k-1}^{\frac{1}{2}}{\widehat{\mathrm{\Σ}}}_{k|k-1}^{\frac{1}{2}}$

还有，定理3的证明将在下面叙述。

上式也可以将K^-，_k(＝P^-1/2K_k)代之以K_k进行整理。

而且如果使用下面的J-酉矩阵

$\mathrm{\Θ}\left(k\right)=(J_1{R}_{e,k}^{\frac{1}{2}}\⊕-R_{r,k}^{\frac{1}{2}})\mathrm{\Σ}\left(k\right)(R_{e,k+1}^{-\frac{1}{2}}{J}_1^{-1}\⊕-R_{r,k+1}^{-\frac{1}{2}})$

，则能够得到定理4的高速化的状态估计算法。其中ψ表示移位矩阵。

定理4

${\hat{x}}_{k|k}={\hat{x}}_{k-1|k-1}+K_{s,k}(y_k-H_k{\hat{x}}_{k-1|k-1})---\left(25\right)$

$K_{s,k}={\mathrm{\ρ}}^{\frac{1}{2}}{\stackrel{\‾}{K}}_k(:,1)/R_{e,k}\left(\mathrm{1,1}\right)---\left(26\right)$

其中，

diag{·}表示对角矩阵，R_e，k+1(1、1)表示矩阵R_e，k+1的1-1分量。又，上式也可以将K^- _k代之以K_k进行整理。

本高速算法利用下面的因式分解，即利用

${\underset{\‾}{\mathrm{\Σ}}}_{k+1|k}-\mathrm{\Ψ}{\underset{\‾}{\mathrm{\Σ}}}_{k|k-1}{\mathrm{\Ψ}}^T=-{\underset{\‾}{L}}_k{R}_{\mathrm{\τ},k}^{-1}{\underset{\‾}{L}}_k^T---\left(32\right)$

中的L~_k∈R(N^+1)×2的更新求滤波器增益K_s，k，因此每单位步骤的计算量为O(N+1)即可。在这里请求、注意下式。

图5表示定理3的数值上稳定的高速算法的流程图。该高速算法加入图2的H_∞滤波器的计算步骤(S105)中，利用γ迭代实现最佳化。因此在存在条件得到满足时，γ_f慢慢减小，而在没有得到满足的时刻如图所示γ_f增加。

H_∞滤波器算法可以摘要叙述如下。

〔步骤S401〕处理部101如图所示那样决定递归式的初始条件。又，L表示最大数据数。

〔步骤S403〕处理部101将时刻k与最大数据数L加以比较。处理部101在时刻k大于最大数据数时结束处理，如果在最大数据数之下则进入下面的步骤。(如果不要，可以取下条件文。又，也可以再度开始。)

〔步骤S405〕处理部101利用式(27)、(31)递归地计算与滤波器增益对应的项K_k+1。

〔步骤S406〕处理部101利用式(29)递归地计算R_e，k+1。

〔步骤S407〕处理部101还利用式(26)、(31)计算K_s，k。

〔步骤S409〕处理部101如果在这里判定存在条件EXC＞0，使存在条件满足，就进入步骤S411。

〔步骤S413〕另一方面，如果在步骤S409不满足存在条件，处理部101就增加γ_f，返回步骤S401。

(步骤S411〕处理部101更新式(25)的H_∞滤波器的滤波器方程式。

〔步骤S415〕处理部101利用式(30)递归地计算R_r，k+1。又，处理部101利用式(28)、(31)递归地计算L~_k+1。

〔步骤S419〕处理部101使时刻k前进(k＝k+1)。返回步骤S403，只要有数据就继续。

还有，处理部101根据需要适当在存储部105存储H_∞滤波器计算步骤S405～S415以及存在条件的计算步骤S409等各步骤求得的适当的中间值以及最终值，也可以将其从存储部105读出。

6.回波消除器

下面建立回波消除问题的数学模型。

首先，如果考虑接收信号{u_k}为向回波通道的输入信号的情况，借助于回波通道的(时变)脉冲响应{h_i〔k〕}，回波{d_k}的观测值用下式表示。

$y_k=d_k+v_k=\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}{h}_i\left[k\right]u_{k-i}+v_k,k=\mathrm{0,1,2},\·\·\·---\left(33\right)$

在这里，u_k、y_k分别表示时刻t_k(＝KT；T是抽样周期)的接收信号与回波，v_k表示时刻t_k的平均值0的线路噪声，h_i〔k〕、i＝0、...、N-1是时变脉冲响应，其抽头数N是已知的。这时，脉冲响应的估计值{h^_i(k〕}如果能够时时刻刻得到，则能够如下所述使用其生成拟似回波。

${\hat{d}}_k=\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}{\hat{h}}_k\left[k\right]u_{k-i},k=\mathrm{0,1,2}\·\·\·---\left(34\right)$

如果将其从回波中减去(y_k-d^_k_0)，则可以消除回波。但是在k-i＜0时，使u_k-1＝0。

根据上面所述，问题归结为从可直接观测的接收信号{u_k}与回波{y_k}逐次估计回波通道的脉冲响应{h_i〔k〕}的问题。

通常，为了把H_∞滤波器使用于回波消除器，首先必须用状态方程式与观测方程式构成的状态空间模型表现式(32)。因此，由于问题是估计脉冲响应{h_i〔k〕}，所以以{h_i〔k〕}为状态变量x_k，如果允许w_k左右的变动，则能够对回波通道使下述状态空间模型成立。

x_k+1＝x_k+G_kw_k，             x_k，w_k∈R^N                   (35)

  y_k＝H_kx_k+v_k，             y_k，v_k∈R                    (36)

  z_k＝H_kx_k，                z_k∈R，H_k∈R^1×N             (37)

其中，

      x_k＝[h₀[k]，…，h_N-1[k]]^T，w_k＝[w_k(1)，…，w_k(N)]^T

      H_k＝[u_k，…，u_k-N+1]

对这样的状态空间模型的超高速和高速H_∞滤波器算法如上所述。又，估计脉冲响应时，通常检测发送信号的发生并且在其间中止估计。

7.对脉冲响应的评价

动作的确认

回波通道的脉冲响应在时间上不变(h_i〔k〕＝h_i)，而且其抽头数N为48的情况下，用仿真对本高速算法进行确认。

$y_k=\underset{i=0}{\overset{47}{\mathrm{\Σ}}}{h}_i{u}_{k-i}+v_k---\left(38\right)$

另外，图6表示脉冲响应{h_i}的值。

在这里，脉冲响应{h_i}_i＝0 ²³采用图示的值，其他{h_i}_i＝24 ⁴⁷采用0。又，v_k是平均值0，方差σ_v ²＝1.0×10^-6是恒定的高斯白噪声，为了方便，抽样周期T取1.0。

又，接收信号{u_k}如下所示用2次的AR模型近似。

u_k＝α₁u_k-1+α₂u_k-2+w_k’

其中α₁＝0.7，α₂＝0.1，w_k’是平均值0，方差σ_w’ ²＝0.04的恒定的高斯白噪声。

脉冲响应的估计结果

图7表示采用定理3的数值上稳定的高速算法得到的脉冲响应的估计结果。在这里，图7(b)的纵轴表示

$\sqrt{}\left\{{{\mathrm{\Σ}}_{i=0}}^{47}{(h_i-{x^{^}}_k(i+1))}^2\right\}$

由此可知，利用本高速算法能够进行良好的估计。但是，采用ρ＝1-χ(γ_f)、χ(γ_f)＝γ_f ^-2、x^_0|0＝0、∑^_1|0＝20I，计算结果以加倍的精度进行。在确认存在条件的同时，设定γ_f＝5.5。

8.定理2的证明

以下关系式成立时，

$\left[\begin{array}{cc}{R}_k^{\frac{1}{2}}& C_k{\widehat{\mathrm{\Σ}}}_{k|k-1}^{\frac{1}{2}}\\ 0& {\mathrm{\ρ}}^{-\frac{1}{2}}{\widehat{\mathrm{\Σ}}}_{k|k-1}^{\frac{1}{2}}\end{array}\right]J\left[\begin{array}{cc}{R}_k^{\frac{1}{2}}& 0\\ {\widehat{\mathrm{\Σ}}}_{k|k-1}^{\frac{1}{2}}{C}_k^T& {\mathrm{\ρ}}^{-\frac{1}{2}}{\widehat{\mathrm{\Σ}}}_{k|k-1}^{\frac{1}{2}}\end{array}\right]$

$=\left[\begin{array}{cc}{R}_{e,k}^{\frac{1}{2}}& 0\\ {\mathrm{\ρ}}^{-\frac{1}{2}}{K}_k{R}_{e,k}^{-\frac{1}{2}}{J}_1^{-1}& {\widehat{\mathrm{\Σ}}}_{k+1|k}^{\frac{1}{2}}\end{array}\right]J\left[\begin{array}{cc}{R}_{e,k}^{\frac{1}{2}}& {\mathrm{\ρ}}^{-\frac{1}{2}}{J}_1^{-1}{R}_{e,k}^{-\frac{1}{2}}{K}_k^T\\ 0& {\widehat{\mathrm{\Σ}}}_{k+1|k}^{\frac{1}{2}}\end{array}\right]---\left(40\right)$

将两边的2×2数据块矩阵的各项加以比较，得到下式。

$R_{e,k}=R_k+C_k{\widehat{\mathrm{\Σ}}}_{k|k-1}{C}_k^T---\left(41\right)$

$K_k={\widehat{\mathrm{\Σ}}}_{k|k-1}{C}_k^T---\left(42\right)$

${\widehat{\mathrm{\Σ}}}_{k+1|k}+{\mathrm{\ρ}}^{-1}{K}_k{R}_{e,k}^{-1}{K}_k^T={\mathrm{\ρ}}^{-1}{\widehat{\mathrm{\Σ}}}_{k|k-1}---\left(43\right)$

这与定理1的F_k＝I时的式(13)的Riccati微分方程式一致。其中，

$J=(J_1\⊕I),J_1\left[\begin{array}{cc}1& 0\\ 0& -1\end{array}\right],C_k=\left[\begin{array}{c}{H}_k\\ H_k\end{array}\right]---\left(44\right)$

另一方面，AJA^T＝BJB^T成立时，B可以利用J-酉矩阵表达为B＝AΘ(k)。因此，根据式(40)，定理1的Riccati微分方程式与下式等价。

$\left[\begin{array}{cc}{R}_{e,k}^{\frac{1}{2}}& 0\\ {\mathrm{\ρ}}^{-\frac{1}{2}}{K}_k{R}_{e,k}^{-\frac{1}{2}}{J}_1^{-1}& {\widehat{\mathrm{\Σ}}}_{k+1|k}^{\frac{1}{2}}\end{array}\right]=\left[\begin{array}{cc}{R}_k^{\frac{1}{2}}& C_k{\widehat{\mathrm{\Σ}}}_{k|k-1}^{\frac{1}{2}}\\ 0& {\mathrm{\ρ}}^{-\frac{1}{2}}{\widehat{\mathrm{\Σ}}}_{k|k-1}^{\frac{1}{2}}\end{array}\right]\mathrm{\Θ}\left(k\right)\·\·\·\left(45\right)$

还有，在式(40)、(45)中，J₁ ^-1可删除。

8-2.定理3的证明

如下所述，假定数据块三角化的J-酉矩阵Θ(k)存在。

这时，如果将上式的两边J＝(J₁_-S)-范数加以比较，可以如下所述决定左边的X、Y、Z。S采用对角分量采取1或-1的对角矩阵。

(1，1)-数据块矩阵

$X{J}_1{X}^T$

$=R_{e,k}+(R_{e,k+1}-R_{k+1})-(R_{e,k}-R_k)=R_{e,k+1}$

因此，从R_e，k+1＝R_e ^1/2 _，k+1 J₁R_e ^1/2 _，k+1，R_k+1＝R_k，得到X＝R_e ^1/2 _，k+1。在这里请注意以下关系成立，即 $J_1^{-1}=J_1(J_1^2=I),S^{-1}=S,R_{e,k+1}^T=R_{e,k+1},R_{r,k}^T=R_{r,k},$

(2，1)-数据块矩阵

$Y{J}_1{X}^T$

$=\left[\begin{array}{c}0\\ {\stackrel{\‾}{K}}_k\end{array}\right]+\left[\begin{array}{c}{\stackrel{\‾}{K}}_{k+1}\\ 0\end{array}\right]-\left[\begin{array}{c}0\\ {\stackrel{\‾}{K}}_k\end{array}\right]$

$=\left[\begin{array}{c}{\stackrel{\‾}{K}}_{k+1}\\ 0\end{array}\right]$

据此得到 $Y=\left[\begin{array}{c}{\stackrel{\‾}{K}}_{k+1}\\ 0\end{array}\right]R_{e,k+1}^{-\frac{1}{2}}{J}_1.$ 其中

(2，2)-数据块矩阵

$-\mathrm{ZS}{Z}^T+Y{J}_1{Y}^T$

$=\left[\begin{array}{c}0\\ {\stackrel{\‾}{K}}_k\end{array}\right]R_{e,k}^{-1}{\left[\begin{array}{c}0\\ {\stackrel{\‾}{K}}_k\end{array}\right]}^T-{\mathrm{\ρ}}^{-1}{\stackrel{\‾}{L}}_k{R}_{r,k}^{-1}{\tilde{L}}_k^T$

因此，

得到 $Z={\tilde{L}}_{k+1}{R}_{r,k+1}^{-\frac{1}{2}}$

8-3.定理4的证明

观测矩阵H_k具有移位特性，而且

                        J＝(J₁_-S)时，利用与定理2相同的方法能够得到下述关系。

$\mathrm{\Θ}\left(k\right)=(J_1{R}_{e,k}^{\frac{1}{2}}\⊕-R_{r,k}^{\frac{1}{2}})\mathrm{\Σ}\left(k\right)(R_{e,k+1}^{-\frac{1}{2}}{J}_1^{-1}\⊕-R_{r,k+1}^{-\frac{1}{2}})$

其中 ，决定R_r，k+1，使下式成立，即

∑(k)^T(R_e，k_-R_r，k)∑(k)＝(R_e，k+1_-R_r，k′+1)

接着，如果在式(46)的第3行新追加R_r，k+1的更新式，则最终得到下式。

从该两边的3×2数据块矩阵的各项的对应得到下面的增益矩阵K^- _k的更新式。

工业应用性

通常的民间的通信设备等，往往从成本和速度面上以单精度进行计算。因此，本发明作为实用性的状态估计算法应该能够给各产业领域带来其效果。又，本发明能够使用于通信***和音响***的回波消除器、音场重放或噪声控制等。

Claims (17)

1.一种***估计方法，是对下式表示的状态空间模型，

x_k+1＝F_kx_k+G_kw_k，

y_k＝H_kx_k+v_k，

z_k＝H_kx_k，

其中，x_k为状态矢量或只是状态，

w_k是***噪声，

v_k是观测噪声，

y_k是观测信号，

z_k是输出信号，

F_k是***的动力学，

G_k是驱动矩阵，

作为评价基准，从用忘却系数ρ加权的干扰决定，而且将向滤波器误差的最大能量增益抑制得比对应于预先被提供的上限值γ_f的项小的估计算法中，同时进行状态估计的可靠性和忘却系数ρ的最佳化用的***估计方法，其特征在于，包含下述步骤，即

处理部从存储部或输入部输入包含上限值γ_f、作为滤波器的输入的观测信号y_k、观测矩阵H_k的值的步骤、

处理部根据上述上限值γ_f决定关于状态空间模型的忘却系数ρ的步骤、

处理部从存储部读取包含初始值或某一时刻的观测矩阵H_k的值，用所述忘却系数ρ执行下式表示的超高速H_∞滤波的步骤、

x^_k|k＝F_k-1x^_k-1|k-1+K_s，k(y_k-H_kF_k-1x^_k-1|k-1)

其中，

x^_k|k：：用观测信号y_o～y_k的时刻k的状态x_k的估计值，

K_s，k：滤波器增益，

处理部存储关于超高速H_∞滤波器的求得的值的存储步骤、

处理部根据求得的观测矩阵H_i或观测矩阵H_i与滤波器增益K_s，i，计算以所述上限值γ_f及所述忘却系数ρ为依据的存在条件的步骤、以及

处理部使γ_f小下去，反复执行所述超高速H_∞滤波的步骤，以此在各时刻在满足所述存在条件的范围将上限值设定得小，将该值存储于存储部的步骤。

2.如权利要求1所述的***估计方法，其特征在于，处理部根据下述公式计算所述存在条件，

${\widehat{\mathrm{\Σ}}}_{i|i}^{-1}={\widehat{\mathrm{\Σ}}}_{i|i-1}^{-1}+\frac{1-{\mathrm{\γ}}_f^{-2}}{\mathrm{\ρ}}{H}_i^T{H}_i>0,i=0,...,k---\left(17\right)$

3.如权利要求1所述的***估计方法，其特征在于，处理部根据下述公式计算所述存在条件，

其中，

4.如权利要求1所述的***估计方法，其特征在于，所述忘却系数ρ与所述上限值γ_f存在下式所述的关系，即

0＜ρ＝1-χ(γ_f)≤1，

其中，χ(γ_f)为满足χ(1)＝1，χ(∞)＝0的γ_f的单调衰减函数。

5.如权利要求1所述的***估计方法，其特征在于，执行所述超高速H∞滤波的步骤，处理部利用下式求得所述滤波器增益K_s，k

${\hat{x}}_{k|k}=F_{k-1}{\hat{x}}_{k-1|k-1}+K_{s,k}(y_k-H_k{F}_{k-1}{\hat{x}}_{k-1|k-1})---\left(11\right)$

$K_{s,k}={\widehat{\mathrm{\Σ}}}_{k|k-1}{H}_k^T\·{(H_k{\widehat{\mathrm{\Σ}}}_{k|k-1}{H}_k^T+\mathrm{\ρ})}^{-1}---\left(12\right)$

$\left.\begin{array}{c}{\widehat{\mathrm{\Σ}}}_{k|k}={\widehat{\mathrm{\Σ}}}_{k|k-1}-{\widehat{\mathrm{\Σ}}}_{k|k-1}{C}_k^T{R}_{e,k}^{-1}{C}_k{\widehat{\mathrm{\Σ}}}_{k|k-1}\\ {\widehat{\mathrm{\Σ}}}_{k+1|k}=\left(F_k{\widehat{\mathrm{\Σ}}}_{k|k}{F}_k^T\right)/\mathrm{\ρ}\end{array}\right\}---\left(13\right)$

其中，

$R_{e,k}=R_k+C_k{\widehat{\mathrm{\Σ}}}_{k|k-1}{C}_k^T,R_k=\left[\begin{array}{cc}\mathrm{\ρ}& 0\\ 0& {-\mathrm{\ρ\γ}}_f^2\end{array}\right],C_k=\left[\begin{array}{c}{H}_k\\ H_k\end{array}\right]---\left(14\right)$

0＜ρ＝1-χ(γ_f)≤1，γ_f＞1                   (15)

还有，式(16)的右边也可以更加一般化，

其中，

x_k为状态矢量或只是状态，

y_k是观测信号，

z_k是输出信号，

F_k是***的动力学，

H_k为观测矩阵，

x^_k|k为用观测信号y_o～y_k的时刻k的状态x_k的估计值，

∑^_k|k：对应于x^_k|k的误差的协方差矩阵，

K_s，k：滤波器增益，

e_f，i：滤波器误差，

R_e，k：辅助变量。

6.如权利要求5所述的***估计方法，其特征在于，执行所述超高速H_∞滤波的步骤包含：

处理部根据初始条件用所述式(12)计算滤波器增益K_s，k的步骤、

处理部更新所述式(11)的H_∞滤波器的滤波器方程式的步骤、

处理部用所述式(13)计算∑^_k|k、∑^_k+1|k的步骤、以及

处理部使时刻k前进，反复执行所述各步骤的步骤。

7.如权利要求1所述的***估计方法，其特征在于，

执行所述超高速H_∞滤波的步骤，

其所述处理部利用增益矩阵K_k，使用下式求所述滤波器增益K_s，k，

${\hat{x}}_{k|k}={\hat{x}}_{k-1|k-1}+K_{s,k}(y_k-H_k{\hat{x}}_{k-1|k-1})---\left(20\right)$

$K_{s,k}=K_k(:,1)/R_{e,k}\left(\mathrm{1,1}\right),K_k={\mathrm{\ρ}}^{\frac{1}{2}}\left({\mathrm{\ρ}}^{-\frac{1}{2}}{K}_k{R}_{s,k}^{-\frac{1}{2}}{J}_1^{-1}\right)J_1{R}_{e,k}^{\frac{1}{2}}---\left(21\right)$

其中，

$R_k=R_k^{\frac{1}{2}}{J}_1{R}_k^{\frac{1}{2}},R_k^{\frac{1}{2}}=\left[\begin{array}{cc}{\mathrm{\ρ}}^{\frac{1}{2}}& 0\\ 0& {\mathrm{\ρ}}^{\frac{1}{2}}{\mathrm{\γ}}_f\end{array}\right],J_1=\left[\begin{array}{cc}1& 0\\ 0& -1\end{array}\right],{\widehat{\mathrm{\Σ}}}_{k|k-1}={\widehat{\mathrm{\Σ}}}_{k|k-1}^{\frac{1}{2}}{\widehat{\mathrm{\Σ}}}_{k|k-1}^{\frac{1}{2}}$

Θ(k)为J-酉矩阵，即满足Θ(k)JΘ(k)^T＝J，J＝(J₁_I)、I是单位矩阵。又，K_k(：，1)表示矩阵K_k的第1列的列矢量，

还有，在这里，在式(21)、(22)中，J₁ ^-1和J₁可删除，

其中，

x^_k|k为用观测信号y_o～y_k的时刻k的状态x_k的估计值，

y_k是观测信号，

F_k是***的动力学，

K_s，k是滤波器增益，

H_k为观测矩阵，

∑^_k|k：对应于x^_k|k的误差的协方差矩阵，

Θ(k)为J-酉矩阵，

R_e，k是辅助变量。

8.如权利要求7所述的***估计方法，其特征在于，执行所述超高速H_∞滤波的步骤包含

处理部用所述式(22)计算K_k、∑^_k+1|k ^1/2的步骤、

处理部根据初始条件用所述式(21)计算滤波器增益K_s，k的步骤、

处理部更新所述式(20)的H_∞滤波器的滤波器方程式的步骤、以及

处理部使时刻k前进，反复执行所述各步骤的步骤。

9.如权利要求1所述的***估计方法，其特征在于，

执行所述超高速H_∞滤波的步骤，

其所述处理部利用增益矩阵K_k，使用下式求所述滤波器增益K_s，k，

${\hat{x}}_{k|k}={\hat{x}}_{k-1|k-1}+{\stackrel{\‾}{K}}_{s,k}(y_k-H_k{\hat{x}}_{k-1|k-1})---\left(61\right)$

${\stackrel{\‾}{K}}_{s,k}={\stackrel{\‾}{K}}_k(:,1)/R_{e,k}\left(\mathrm{1,1}\right),{\stackrel{\‾}{K}}_k={\mathrm{\ρ}}^{\frac{1}{2}}\left({\stackrel{\‾}{K}}_k{R}_{e,k}^{-\frac{1}{2}}\right)R_{e,k}^{\frac{1}{2}}---\left(62\right)$

其中，Θ(k)为任意的J-酉矩阵，C^_k＝C^_k+1Ψ成立。

其中，

$R_k=R_k^{\frac{1}{2}}{J}_1{R}_k^{\frac{1}{2}},R_k^{\frac{1}{2}}=\left[\begin{array}{cc}{\mathrm{\ρ}}^{\frac{1}{2}}& 0\\ 0& {\mathrm{\ρ}}^{\frac{1}{2}}{\mathrm{\γ}}_f\end{array}\right],J_1=\left[\begin{array}{cc}1& 0\\ 0& -1\end{array}\right],{\widehat{\mathrm{\Σ}}}_{k|k-1}={\widehat{\mathrm{\Σ}}}_{k|k-1}^{\frac{1}{2}}{\widehat{\mathrm{\Σ}}}_{k|k-1}^{\frac{1}{2}}$

其中

x^_k|k为用观测信号y_o～y_k的时刻k的状态x_k的估计值，

y_k是观测信号，

K_s，k是滤波器增益，

H_k为观测矩阵，

Θ(k)为J-酉矩阵，

R_e，k是辅助变量。

10.如权利要求9所述的***估计方法，其特征在于，

执行所述超高速H_∞滤波的步骤包含

处理部用所述式(63)计算K_k ^-的步骤、

处理部根据初始条件用所述式(62)计算滤波器增益K^- _s，k的步骤、

处理部更新所述式(61)的H_∞滤波器的滤波器方程式的步骤、以及

处理部使时刻k前进，反复执行所述各步骤的步骤。

11.如权利要求1所述的***估计方法，其特征在于，

执行所述超高速H_∞滤波的步骤，

其所述处理部利用增益矩阵K_k ^-，使用下式求所述滤波器增益K_s，k，

${\hat{x}}_{k|k}={\hat{x}}_{k-1|k-1}+K_{s,k}(y_k-H_k{\hat{x}}_{k-1|k-1})---\left(25\right)$

$K_{s,k}={\mathrm{\ρ}}^{\frac{1}{2}}{\stackrel{\‾}{K}}_k(:,1)/R_{e,k}\left(\mathrm{1,1}\right)---\left(26\right)$

其中，

还有，上式也可以用K_k代替K_k-进行整理，

K，，其中，

y_k是观测信号，

F_k是***的动力学，

H_k为观测矩阵，

x^_k|k为用观测信号y_o～y_k的时刻k的状态x_k的估计值，

K_s，k为滤波器增益，从增益矩阵K_k ^-求得，

R_e，k、L~_k：辅助变量。

2.如权利要求11所述的***估计方法，其特征在于，

执行所述超高速H_∞滤波的步骤包含

处理部根据预先决定的初始条件用所述式(27)递归地计算K^- _k+1的步骤、

处理部用所述式(26)计算***增益K_s，k的步骤、

处理部计算存在条件的步骤、

如果使所述存在条件满足，处理部就更新所述式(25)的H_∞滤波器的滤波器方程式，使时刻k前进，反复执行所述各步骤的步骤、以及

如果没有使所述存在条件满足，就增加上限值γ^f的步骤。

13.如权利要求1所述的***估计方法，其特征在于，而且利用下式从时刻k的状态估计值x^_k|k求输出信号的估计值zˇ_k|k，

zˇ_k|k＝H_kx^_k|k。

14.如权利要求1所述的***估计方法，其特征在于，

使用所述H_∞滤波器方程式，求状态估计值x^_k|k，

对拟似回波如下所述进行估计，

用求得的拟似回波抵消实际回波，以此实现回波消除器，

${\hat{d}}_k=\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}{\hat{h}}_i\left[k\right]u_{k-i,}k=\mathrm{0,1,2},\·\·\·---\left(34\right)$

15.一种***估计程序，是对下式表示的状态空间模型，

x_k+1＝F_kx_k+G_kw_k，

y_k＝H_kx_k+v_k，

z_k＝H_kx_k，

其中，x_k为状态矢量或只是状态，

w_k是***噪声，

v_k是观测噪声，

y_k是观测信号，

z_k是输出信号，

F_k是***的动力学，

G_k是驱动矩阵，

作为评价基准，从用忘却系数ρ加权的干扰决定，而且将向滤波器误差的最大能量增益抑制得比对应于预先被提供的上限值γ_f的项小的估计算法中，同时使计算机执行状态估计的可靠性和忘却系数ρ的最佳化用的***估计程序，其特征在于，使计算机执行下述步骤，即

处理部从存储部或输入部输入包含上限值γ_f、作为滤波器的输入的观测信号y_k、观测矩阵H_k的值的步骤、

处理部根据上述上限值γ_f决定关于状态空间模型的忘却系数ρ的步骤、

处理部从存储部读取包含初始值或某一时刻的观测矩阵H_k的值，用所述忘却系数ρ执行下式表示的的超高速H_∞滤波的步骤、

x^_k|k＝F_k-1x^_k-1|k-1+K_s，k(y_k-H_kF_k-1x^_k-1|k-1)

其中，

x^_k|k：用观测信号y_o～y_k的时刻k的状态x_k的估计值，

F_k：***的动力学，

K_s，k：滤波器增益，

处理部存储关于超高速H_∞滤波器的求得的值的存储步骤、

处理部根据求得的观测矩阵H_i或观测矩阵H_i与滤波器增益K_s，i计算以所述上限值γ_f及所述忘却系数ρ为依据的存在条件的步骤、以及

处理部使γ_f小下去，反复执行所述超高速H_∞滤波的步骤，以此在各时刻在满足所述存在条件的范围将上限值设定得小，将该值存储于存储部的步骤。

16.一种记录使计算机执行用的***估计程序的计算机可读取的记录媒体，是对下式表示的状态空间模型，

x_k+1＝F_kx_k+G_kw_k，

y_k＝H_kx_k+v_k，

z_k＝H_kx_k，

其中，x_k为状态矢量或只是状态，

w_k是***噪声，

v_k是观测噪声，

y_k是观测信号，

z_k是输出信号，

F_k是***的动力学，

G_k是驱动矩阵，

作为评价基准，从用忘却系数ρ加权的干扰决定，而且将滤波器误差的最大能量增益抑制得比对应于预先被提供的上限值γ_f的项小的估计算法中，同时使计算机执行状态估计的可靠性和忘却系数ρ的最佳化同时使计算机执行的***估计程序的记录用的计算机可读取记录媒体，其特征在于，记录使计算机执行下述步骤用的***估计程序的计算机可读取的记录媒体，所述步骤即

处理部从存储部或输入部输入包含上限值γ₁、作为滤波器的输入的观测信号y_k、观测矩阵H_k的值的步骤、

处理部根据上述上限值γ_f决定关于状态空间模型的忘却系数ρ的步骤、

处理部从存储部读取包含初始值或某一时刻的观测矩阵Hk的值，用所述忘却系数ρ执行下式表示的的超高速H_∞滤波的步骤、

${\hat{x}}_{k|k}={\hat{x}}_{k-1|k-1}+K_{s,k}(y_k-H_k{\hat{x}}_{k-1|k-1})---\left(25\right)$

$K_{s,k}={\mathrm{\ρ}}^{\frac{1}{2}}{\stackrel{\‾}{K}}_k(:,1)/R_{e,k}\left(\mathrm{1,1}\right)---\left(26\right)$

其中，

x^_k|k：用观测信号y_o～y_k的时刻k的状态x_k的估计值，

F_k：***的动力学，

K_s，k：滤波器增益，

处理部存储关于超高速H_∞滤波器的求得的值的存储步骤、

，，处理部根据求得的观测矩阵H_i或观测矩阵H_i与滤波器增益K_s，j，计算以所述上限值γ_f及所述忘却系数ρ为依据的存在条件的步骤、以及

处理部使γ_f小下去，反复执行所述超高速H_∞滤波的步骤，以此在各时刻在满足所述存在条件的范围将上限值设定得小，将该值存储于存储部的步骤。

17.一种***估计装置，是对下式表示的状态空间模型，

x_k+1＝F_kx_k+G_kw_k，

y_k＝H_kx_k+v_k，

z_k＝H_kx_k，

其中，x_k为状态矢量或只是状态，

w_k是***噪声，

v_k是观测噪声，

y_k是观测信号，

z_k是输出信号，

F_k是***的动力学，

G_k是驱动矩阵，

作为评价基准，从用忘却系数ρ加权的干扰决定，而且将滤波器误差的最大能量增益抑制得比对应于预先被提供的上限值γ_f的项小的估计算法中，同时进行状态估计的可靠性和忘却系数ρ的最佳化用的***估计装置，其特征在于，具备

执行估计算法的处理部、以及利用所述处理部进行读出和/或写入，存储与状态空间模型相关的各观测值、设定值、估计值的存储部，

所述处理部从存储部或输入部输入包含上限值γ₁、作为滤波器的输入的观测信号y_k、观测矩阵H_k的值、

所述处理部根据上述上限值γ_f决定关于状态空间模型的忘却系数ρ、

所述处理部从存储部读取包含初始值或某一时刻的观测矩阵H_k的值，用所述忘却系数ρ执行下式表示的的超高速H_∞滤波、

x^_k|k＝F_k-1x^_k-1|k-1+K_s，k(y_k-H_kF_k-1x^_k-1|k-1)

其中，

x^_k|k：：用观测信号y_o～y_k的时刻k的状态x_k的估计值，

F_k-1：***的动力学，

K_s，k：滤波器增益，

所述处理部存储关于超高速H_∞滤波器的求得的值、

所述处理部根据求得的观测矩阵H_i或观测矩阵H_i与滤波器增益K_s，i，计算以所述上限值γ_f及以所述忘却系数ρ为依据的存在条件、而且所述处理部使γ_f小下去，反复执行所述超高速H_∞滤波的步骤，以此在各时刻在满足所述存在条件的范围将上限值设定得小，将该值存储于存储部。

Images