Orthogonal Rotation-invariant V Moments and Application to Image Reconstruction
-
摘要: 定义在单位圆盘上的正交旋转不变矩函数(如Zernike矩) 具有非常广泛的应用. 本文基于一类正交分段多项式函数系--V系统, 构造了一种新型的矩函数, 称之为正交旋转不变V矩(简称为V矩). 除了正交性、旋转不变性之外, 由于V系统具有次数低、表达式简单的优点, V矩能够避免传统矩函数中高阶多项式的计算, 从而能够保证数值稳定性, 降低计算复杂度. 实验结果表明, V矩比传统的正交旋转不变矩具有更好的图像重建与图像检索结果.Abstract: Orthogonal rotation-invariant moments (ORIMs) such as Zernike moments introduced and defined on a continuous unit disk have been proven to be powerful tools in many applications. A new kind of ORIMs named V moments (VMs) is introduced in this paper, which is based on a class of orthogonal piecewise polynomial functions--V system. In addition to orthogonality and rotation-invariant, VMs have two special advantages of numerical stability and low computational complexity, which originate from the low degrees and simple expressions of the V system. Experimental results are presented to show that VMs can have a better image reconstruction and image retrieval results than the traditional moments.
-
子空间辨识算法由于能对多输入多输出系统采用统一框架建立状态空间模型,在系统辨识和控制工程领域受到广泛关注 [1]. 有一些子空间算法被提出用于开环辨识工业过程,得到一致估计结果 [2]. 但是,由于过程操作安全性和稳定性的需要,许多工业过程限制在闭环条件下进行辨识实验,由于反馈控制作用的影响,使得过程噪声和输入存在相关性,使得传统开环子空间方法产生辨识偏差 [3]. 闭环系统辨识因此在近年来受到很多关注和探讨 [4],一些闭环子空间辨识算法被相继提出 [5]. 这些算法可被归纳为三类,第一类方法 [6-7]采用辅助变量策略消除噪声影响,保证估计结果的一致性;第二类方法 [8]采用最小二乘法对高阶VARX模型(Vector autoregressive with exogenous inputs model)进行计算得到马尔科夫估计参数,由于VARX模型只包括当前时刻的不可测噪声,该噪声和VARX模型的过去时刻输入无关,从而可保证所得参数的一致性;第三类算法 [9]用噪声预估值代替真实值进行计算保证得到一致估计结果.
基于奇偶空间的闭环子空间辨识算法SIMPCAwc [6]采用过去时刻输入输出数据作为辅助变量来消除噪声,以得到无噪声的输入输出数据,然后从无噪声影响的输入输出数据奇偶空间中提取得到扩展可观测矩阵和下三角形Toeplitz矩阵,从而求得系统矩阵,该方法取得较好辨识精度.然而,文献[7]指出当闭环系统设定点输入激励为不相关白噪声序列时,虽然引入辅助变量与噪声不相关,可以有效地消除噪声,但由于该辅助变量和系统设定点输入激励也不相关,导致从无噪声输入输出数据奇偶空间中提取参数可能同时含有过程模型参数信息和控制器参数信息,因而无法对它们进行区分,从而致使过程模型估计出现偏差.针对SIMPCAwc [6]辨识方法存在的问题,本文通过将输入输出数据正交投影到新息数据的正交补空间来消除噪声,以得到新的无噪声数据矩阵,进而从其对应的奇偶空间中提取得到扩展可观测矩阵和下三角形Toeplitz矩阵. 由于新息数据的正交补空间数据和噪声无关,且同时与系统设定输入激励相关,确保本文方法从新的无噪声奇偶空间中提取的参数只包含过程模型参数,有效地保证估计结果的一致无偏性.由于新息数据为不可测量数据,本文通过模型推导,得到和待辨识状态空间模型等价的VARX模型. 在此基础上,采用最小二乘法对VARX模型进行计算以得到新息的一致估计值. 采用新息一致估计值代替真实值,以完成模型参数估计.为了论证说明本文方法的有效性,严格分析和给出了本文算法保证一致估计的条件.
1. 问题描述
本文研究如下线性离散状态空间过程模型:
$S:\left\{ \begin{align} & x(t+1)=Ax(t)+Bu(t)+w(t) \\ & y(t)=Cx(t)+Du(t)+v(t) \\ \end{align} \right.$
(1) 其中,$x(t)\in {{R}^{{{n}_{x}}}}$,$u(t)\in {{R}^{{{n}_{u}}}}$,$y(t)\in {{R}^{{{n}_{y}}}}$分别为系统状态和过程输入和输出. $v(t)\in {{R}^{{{n}_{y}}}}$和$w(t)\in {{R}^{{{n}_{x}}}}$分别为过程测量噪声. A,B,C,D 分别为相应维数的系统矩阵.本文研究系统在闭环工作条件下,利用系统输入和输出观测数据,辨识对象状态空间(亦称子空间)模型.
由于闭环反馈控制作用的影响,使得过程测量噪声和输入存在相关性,若直接通过模型(1)来辨识系统矩阵,很难消除噪声对辨识结果的不利影响.因此,采用卡尔曼滤波原理 [10],将系统模型(1)等价表示为新息形式
${{S}_{I}}:\left\{ \begin{array}{*{35}{l}} \begin{align} & x(t+1)=Ax(t)+Bu(t)+Ke(t) \\ & y(t)=Cx(t)+Du(t)+e(t) \\ \end{align} \\ \end{array} \right.$
(2) 其中,K为卡尔曼滤波增益,新息$e(t)$为零均值白噪声,当$i＜j$时,新息$e(j)$和输入输出$\{u(i),y(i)\}$不相关.
进一步定义$\bar{A}=A-KC$和$\bar{B}=B-KD$,模型(2)可被等价描述为如下预测形式:
${{S}_{P}}:\left\{ \begin{array}{*{35}{l}} \begin{align} & x(t+1)=\bar{A}x(t)+\bar{B}u(t)+Ky(t) \\ & y(t)=Cx(t)+Du(t)+e(t) \\ \end{align} \\ \end{array} \right.$
(3) 其中,假设$\bar{A}$的特征值严格位于单位圆内.
定义过去和将来水平数分别为p和f,过去和将来输入向量分别为$u_p(t)=[u(t-p)^{\textrm T}$ $\cdots$ $u(t-2)^{\textrm T}$ $u(t-1)^{\textrm T}]^{\textrm T}$和$u_f(t)=[u(t)^{\textrm T}$ $\cdots$ $u(t+f-2)^{\textrm T}$ $u(t+f-1)^{\textrm T}]^{\textrm T}$,定义过去和将来输入Hankel 矩阵Up $=$ $[u_p(t)^{\textrm T}$ $\cdots$ $u_p(N)^{\textrm T}]^{\textrm T}$和$U_f=[u_f(t)^{\textrm T}$ $\cdots$ $u_f(N)^{\textrm T}]^{\textrm T}$,输出和新息数据做类似定义.
对式(3)进行迭代可得:
$x(t)={{\bar{A}}^{p}}x(t-p)+{{\bar{L}}_{1}}{{u}_{p}}(t)+{{\bar{L}}_{2}}{{y}_{p}}(t)$
(4) 其中,扩展可观性矩阵分别表示为 $\bar{L}_1=[\bar{A}^{p-1}\bar{B}$ $\cdots$ $\bar{A}\bar{B}$ $\bar{B}]$,$\bar{L}_2=[\bar{A}^{p-1}K$ $\cdots$ $\bar{A}K$ $K]$.初始状态为$x(t$ $-$ $p)$.当p充分大时,$x(t-p)$可被忽略,将式(4)带入式(3})得到等价VARX模型
${{S}_{V}}:y(t)=C{{\bar{L}}_{1}}{{u}_{p}}(t)+C{{\bar{L}}_{2}}{{y}_{p}}(t)+e(t)$
(5) 本文将采用模型(2)和(5)对系统矩阵进行辨识.定义$X_p=[x(t-p)$ $\cdots$ $x(t-p+N-1)]$ 和Xf $=$ $[x(t)$ $\cdots$ $x(t+N-1)]$.通过对式(2)进行迭代可得:
$Y(t)={{\Gamma }_{f}}{{X}_{f}}+{{H}_{f}}{{U}_{f}}+{{G}_{f}}{{E}_{f}}$
(6) 其中,扩展可观测矩阵为$\Gamma_f=[C^{\textrm T}$ $\cdots$ $(CA^{f-1})^{\textrm T}]^{\textrm T}$.下三角形Toeplitz矩阵分别为
$\begin{align} & {{H}_{f}}=\left[ \begin{matrix} D & \cdots & \cdots & 0 \\ CB & \cdots & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ C{{A}^{f-2}} & \cdots & CB & D \\ \end{matrix} \right] \\ & {{G}_{f}}=\left[ \begin{matrix} 0 & \cdots & \cdots & 0 \\ C & \cdots & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ C{{A}^{f-2}} & \cdots & C & 0 \\ \end{matrix} \right] \\ \end{align}$
2. 闭环子空间辨识算法和一致性分析
2.1 闭环子空间辨识算法
在式(6)的基础上,通过同时计算扩展可观测矩阵$\Gamma_f$和下三角形Toeplitz矩阵Hf实现对系统矩阵的辨识.首先将输入数据移至式(6)的左侧,得到:
$[I-{{H}_{f}}]{{W}_{f}}={{\Gamma }_{f}}{{X}_{f}}+{{G}_{f}}{{E}_{f}}$
(7) 其中,$W_f=[Y_f^{\textrm T}U_f^{\textrm T}]^{\textrm T}$,对Wp做同样定义.
为求解式(7)得到$\Gamma_f$和Hf的估计值,需要消除未知状态和新息的影响.通过在式(7)的左右侧同时引入$\Gamma_f$的正交补向量$\Gamma_f^{\bot}$,由于$(\Gamma_f^{\bot})^{\textrm T}\Gamma_f=0$
$\left[ {{(\Gamma _{f}^{\bot })}^{\text{T}}}-{{(\Gamma _{f}^{\bot })}^{\text{T}}}{{H}_{f}} \right]{{W}_{f}}={{(\Gamma _{f}^{\bot })}^{\text{T}}}{{G}_{f}}{{E}_{f}}$
(8) 通过将式(8)正交投影到Ef的正交补空间来消除新息噪声,
$\left[ {{(\Gamma _{f}^{\bot })}^{\text{T}}}-{{(\Gamma _{f}^{\bot })}^{\text{T}}}{{H}_{f}} \right]{{W}_{f}}\Pi _{{{E}_{f}}}^{\bot }=0$
(9) 其中,$\Pi_{E_f}^{\bot}=I_N-E_f^{\textrm T}(E_fE_f^{\textrm T})^{-1}E_f$是Ef的正交补.由于$\lim_{N \to \infty }E_f\Pi_{{E}_f}^{\bot}=0$和$\lim_{N \to \infty }R_f\Pi_{{E}_f}^{\bot}=R_f$ $\neq$ $0$.根据文献[7]结论可知,无噪声数据块$W_f\Pi_{E_f}^{\bot}$的奇偶空间只包含过程模型参数信息,不会包括控制器模型参数信息.因此,通过$W_f\Pi_{E_f}^{\bot}$的奇偶空间可得到过程模型参数的一致估计值.
由于新息Ef未知,不能进行模型参数估计.这里利用等价辅助模型(5)计算Ef的估计值$\hat{E}_f$,再采用估计值代替真实值进行后续计算.定义新的过去输入输出Hankel 矩阵$U_p(t,N+f)=[u_p(t)$ $\cdots$ $u_p(N+f)]$,$Y_p(t,N+f)=[y_p(t)$ $\cdots$ $y_p(N+f)]$,新息数据做同样定义.定义$W_p(t,N+f)=$ $[U_p^{\textrm T}(t,$ $N+f)$ $Y_p^{\textrm T}(t,N+f)]^{\textrm T}$ 和$Y(t,N+f)=$ $[y(t)$ $\cdots$ $y(N+f)]$.同时定义 $\theta= [C\bar{L}_1$ $C\bar{L}_2]$.采用最小二乘法对式(5)进行计算,可得:
$\hat{\theta }=Y(t,N+f)W_{p}^{\dagger }(t,N+f)$
(10) 其中,
$\begin{align} & W_{p}^{\dagger }(t,N+f)= \\ & W_{p}^{\text{T}}(t,N+f){{[{{W}_{p}}(t,N+f)W_{p}^{\text{T}}(t,N+f)]}^{-1}} \\ \end{align}$
由于估计值$\hat{\theta}$和真实值的误差为
$\Delta \theta =\hat{\theta }-\theta =E(t,N+f)W_{p}^{\dagger }(t,N+f)$
(11) 其中,
$E(t,N+f)=[e(t)\cdots (N+f)]$
由于$\lim_{N \to \infty }E(t,N+f)W_p^{\textrm T}(t,N+f)=0$,若$\lim_{N \to \infty }W_p(t,N+f)W_p^{\textrm T}(t,N+f)>0$ (可作为默认成立条件),则$\hat{\theta}$是一致估计值.
将一致估计值$\hat{\theta}$带入估计值$\hat{Y}(t,N+f)$,可以得到一致估计值$\hat{E}(t,N+f)$ (证明见 第2.2节).
$\begin{align} & \hat{E}(t,N+f)=Y(t,N+f)-\hat{Y}(t,N+f)= \\ & Y(t,N+f)-\hat{\theta }{{W}_{p}}(t,N+f)= \\ & Y(t,N+f)[{{I}_{N+f}}-W_{p}^{\dagger }(t,N+f){{W}_{p}}(t,N+f)]= \\ & Y(t,N+f)\Pi _{{{W}_{p}}(t,N+f)}^{\bot } \\ \end{align}$
(12) 其中,
$\begin{align} & \Pi _{{{W}_{p}}(t,N+f)}^{\bot }={{I}_{N+f}}-W_{p}^{\text{T}}(t,N+f)\times \\ & {{\left[ {{W}_{p}}(t,N+f)W_{p}^{\text{T}}(t,N+f) \right]}^{-1}}{{W}_{p}}(t,N+f) \\ \end{align}$
定义
$\Pi _{{{W}_{p}}(t,N+f)}^{\bot }=[{{\beta }_{1}}\cdots {{\beta }_{N+f}}]$
(13) 则未来噪声Hankel矩阵的一致估计值可通过以下方式重构得到:
${{{\hat{E}}}_{f}}=Y(t,N+f)\left[ \begin{matrix} {{\beta }_{1}} & \cdots & {{\beta }_{N}} \\ {{\beta }_{2}} & \cdots & {{\beta }_{N+1}} \\ \vdots & \ddots & \vdots \\ {{\beta }_{f}} & \cdots & {{\beta }_{f+N}} \\ \end{matrix} \right]$
(14) 进一步将Ef用其一致估计值代替,可得一致估计值$W_f\Pi_{\hat{E}_f}^{\bot}$.通过SVD分解得到:
$\begin{align} & {{W}_{f}}\Pi _{{{{\hat{E}}}_{f}}}^{\bot }= \\ & \left[ {{{\hat{U}}}_{1}}\hat{U}_{1}^{\bot } \right]\left[ \begin{matrix} {{{\hat{\Sigma }}}_{1}} & 0 \\ {{{\hat{\Sigma }}}_{2}} & 0 \\ \end{matrix} \right]\left[ \begin{matrix} \begin{matrix} \hat{V}_{1}^{\text{T}} \\ {{(\hat{V}_{1}^{\bot })}^{\text{T}}} \\ \end{matrix} \\ \end{matrix} \right]={{{\hat{U}}}_{1}}{{{\hat{\Sigma }}}_{1}}\hat{V}_{1}^{\text{T}} \\ \end{align}$
(15) 其中,$\hat{U}_1$是$W_f\Pi_{\hat{E}_f}^{\bot}$的前$n_x+fn_u$个特征向量,则$[(\Gamma_f^{\bot})^{\textrm T}-(\Gamma_f^{\bot})^{\textrm T}H_f]$的一致估计值如下(证明见第2.2节):
$\left[ {{(\hat{\Gamma }_{f}^{\bot })}^{\text{T}}}-{{(\hat{\Gamma }_{f}^{\bot })}^{\text{T}}}{{{\hat{H}}}_{f}} \right]={{(\hat{U}_{1}^{\bot })}^{\text{T}}}$
(16) 定义$\hat{U}_1^{\bot}=[P_1^{\textrm T}P_2^{\textrm T}]^{\textrm T}$,其中,$P_1=\hat{U}_1^{\bot}(1:fn_y$,$:)$,$P_2=\hat{U}_1^{\bot}(1+fn_y:end,:)$ (Matlab表示).可得:
$\hat{\Gamma }_{f}^{\bot }={{P}_{1}}$
(17) $-{{(\hat{\Gamma }_{f}^{\bot })}^{\text{T}}}{{{\hat{H}}}_{f}}=P_{2}^{\text{T}}$
(18) 由于$(\Gamma_f^{\bot})^{\bot}=\Gamma_f$,根据式(17),可得:
${{{\hat{\Gamma }}}_{f}}={{(\hat{\Gamma }_{f}^{\bot })}^{\bot }}={{I}_{f{{n}_{y}}}}-{{P}_{1}}{{(P_{1}^{\text{T}}{{P}_{1}})}^{-1}}P_{1}^{\text{T}}$
(19) 为了得到Hf的估计值,做如下定义:
$-P_{1}^{\text{T}}=[{{\phi }_{1}}\cdots {{\phi }_{f}}]$
(20) $P_{2}^{\text{T}}=[{{\varphi }_{1}}\cdots {{\varphi }_{f}}]$
(21) 其中,${{\phi }_{i}}\in {{R}^{(i{{n}_{y}}-{{n}_{x}})\times {{n}_{y}}}},{{\varphi }_{i}}\in {{R}^{(i{{n}_{y}}-{{n}_{x}})\times {{n}_{u}}}}$.
将式(20)和式(21)带入式(18),可得:
$\left[ \begin{matrix} {{\phi }_{1}} & \cdots & {{\phi }_{f-1}} & {{\phi }_{f}} \\ {{\phi }_{2}} & \cdots & {{\phi }_{f}} & 0 \\ \vdots & \vdots & \ddots & \vdots \\ {{\phi }_{f}} & 0 & \cdots & 0 \\ \end{matrix} \right]{{H}_{f1}}=\left[ \begin{matrix} {{\varphi }_{1}} \\ \phi {{i}_{2}} \\ \vdots \\ {{\phi }_{f}} \\ \end{matrix} \right]$
(22) 其中,$H_{f1}=[D^{\textrm T}(CB)^{\textrm T}\cdots(CA^{f-2}B)^{\textrm T}]$.
采用最小二乘法可得$H_{f1}$的一致估计如下:
${{{\hat{H}}}_{f1}}={{\left[ \begin{matrix} {{\phi }_{1}} & \cdots & {{\phi }_{f-1}} & {{\phi }_{f}} \\ {{\phi }_{2}} & \cdots & {{\phi }_{f}} & 0 \\ \vdots & \vdots & \ddots & \vdots \\ {{\phi }_{f}} & 0 & \cdots & 0 \\ \end{matrix} \right]}^{\dagger }}\left[ \begin{matrix} {{\varphi }_{1}} \\ \phi {{i}_{2}} \\ \vdots \\ {{\phi }_{f}} \\ \end{matrix} \right]$
(23) 系统矩阵估计值$\hat{A}$和$\hat{C}$可直接从$\hat{\Gamma}$中提取,即
$\hat{C}=\hat{\Gamma }(1:{{n}_{y}},1:{{n}_{x}})$
(24) $\hat{A}={{{\hat{\Gamma }}}^{\dagger }}(1+{{n}_{y}}:f{{n}_{y}},1:{{n}_{x}})\hat{\Gamma }(1:(f-1){{n}_{y}},1:{{n}_{x}})$
(25) 由于
${{H}_{f1}}=\left[ \begin{matrix} {{I}_{{{n}_{y}}}} & 0 \\ 0 & \Gamma (1:(f-1){{n}_{y}},1:{{n}_{x}}) \\ \end{matrix} \right]\left[ \begin{matrix} \begin{matrix} D \\ B \\ \end{matrix} \\ \end{matrix} \right]$
(26) 系统矩阵估计值$\hat{B}$和$\hat{D}$可采用最小二乘法从$\hat{H}_{f1}$中计算得到:
$\left[ \begin{matrix} \begin{matrix} {\hat{D}} \\ {\hat{B}} \\ \end{matrix} \\ \end{matrix} \right]={{\left[ \begin{matrix} {{I}_{{{n}_{y}}}} & 0 \\ 0 & \hat{\Gamma }(1:(f-1){{n}_{y}},1:{{n}_{x}}) \\ \end{matrix} \right]}^{\dagger }}{{{\hat{H}}}_{f1}}$
(27) 本文基于新息估计和正交投影的闭环子空间辨识方法(Closed-loop subspace identification method using innovation estimation and orthogonal projection,CSIMIEOP)可总结如下:
步骤 1. 由式(10)求解$\hat{\theta}$,再由式(12)计算$\hat{E}(t,N+f)$,采用式(14)构造$\hat{E}_f$.
步骤 2. 通过SVD分解对式(15)进行计算.
步骤 3. 通过式(19})和式(23)计算估计值$\hat{\Gamma}$和$\hat{H}_{f1}$.
步骤 4. 通过式(24),式(25)和式(27)求解系统矩阵$\hat{A}$,$\hat{B}$,$\hat{C}$ 和 $\hat{D}$.
2.2 闭环一致条件分析
由第2.1节闭环子空间辨识算法可知,本文采用新息估计和正交投影消除噪声,得到一致估计结果. 辨识结果是否一致取决于$\hat{E}(t,N+f)$ 和$[(\hat{\Gamma}_f^{\bot})^{\textrm T}$ $-$ $(\hat{\Gamma}_f^{\bot})^{\textrm T}\hat{H}_f]$是否一致,本文对它们的一致估计条件进行分析和说明,给出如下定理.
定理 1. 若$$\lim_{N \to \infty }W_p(t,N+f)W_p^{\textrm T}(t,N+f)>0$$ 则$\hat{E}(t,N+f)$ 为一致估计值.
证明. 由式(12)可知,估计值$\hat{E}(t,N+f)$和真实值的误差为
$\begin{align} & \Delta E(t,N+f)=\hat{E}(t,N+f)-E(t,N+f)= \\ & Y(t,N+f)\Pi _{{{W}_{p}}(t,N+f)}^{\bot }-E(t,N+f)= \\ & [Y(t,N+f)-\hat{Y}(t,N+f)]-E(t,N+f)= \\ & [\theta -\hat{\theta }]{{W}_{p}}(t,N+f)= \\ & \Delta \theta {{W}_{p}}(t,N+f)= \\ & E(t,N+f){{W}_{p}}{{(t,N+f)}^{\dagger }}{{W}_{p}}(t,N+f) \\ \end{align}$
(28) 由于
$\underset{N\to \infty }{\mathop{\lim }}\,E(t,N+f)W_{p}^{\text{T}}(t,N+f)=0$
若
$\underset{N\to \infty }{\mathop{\lim }}\,{{W}_{p}}(t,N+f)W_{p}^{\text{T}}(t,N+f)＞0$
则
$\underset{N\to \infty }{\mathop{\lim }}\,\Delta E(t,N+f)=0$
从而可知$\hat{E}(t,N+f)$为一致估计值.
定理 2. 若系统可控可观,且$E_fE_f^{\textrm T}>0$ 和
$\bar{E}\left\{ \left[ \begin{matrix} {{u}_{p}}(t) \\ {{u}_{f}}(t) \\ {{e}_{p}}(t) \\ {{e}_{f}}(t) \\ \end{matrix} \right]{{\left[ \begin{matrix} {{u}_{p}}(t) \\ {{u}_{f}}(t) \\ {{e}_{p}}(t) \\ {{e}_{f}}(t) \\ \end{matrix} \right]}^{\text{T}}} \right\}＞0$
则$[(\hat{\Gamma}_f^{\bot})^{\textrm T} -(\hat{\Gamma}_f^{\bot})^{\textrm T}\hat{H}_f]$为一致估计值.
证明. 令$\Omega =\left[ \begin{matrix} {{\Gamma }_{f}} & {{H}_{f}} \\ 0 & I \\ \end{matrix} \right]$,由式(9)可知:
$\Omega \left[ \begin{matrix} \begin{matrix} {{U}_{f}} \\ {{X}_{f}} \\ \end{matrix} \\ \end{matrix} \right]\left[ \begin{matrix} {{I}_{N}}-E_{f}^{\text{T}}{{({{E}_{f}}E_{f}^{\text{T}})}^{-1}}{{E}_{f}} \\ \end{matrix} \right]$
(29) 由式(29)可得:
$\begin{array}{l} \mathop {\lim }\limits_{N \to 0} {W_f}\Pi _{{{\hat E}_f}}^ \bot W_f^{\rm{T}} = \\ \Omega \left[ {\begin{array}{*{20}{c}} \begin{array}{l} {U_f}\\ {X_f} \end{array} \end{array}} \right]\left[ {{I_N} - E_f^{\rm{T}}{{({E_f}E_f^{\rm{T}})}^{ - 1}}{E_f}} \right]{\left[ {\begin{array}{*{20}{c}} \begin{array}{l} {U_f}\\ {X_f} \end{array} \end{array}} \right]^{\rm{T}}}{\Omega ^{\rm{T}}} = \\ \Omega \left\{ {{R_1} - {R_2}{{({E_f}E_f^{\rm{T}})}^{ - 1}}R_2^{\rm{T}}} \right\}{\Omega ^{\rm{T}}} \end{array}$
(30) 其中,
$\begin{array}{l} {R_1} = \left[ {\begin{array}{*{20}{c}} {{U_f}{X_f}} \end{array}} \right]{\left[ {\begin{array}{*{20}{c}} {{U_f}{X_f}} \end{array}} \right]^{\rm{T}}}\\ = \bar E\left\{ {\left[ {\begin{array}{*{20}{c}} {{u_f}(t)}\\ {x(t)} \end{array}} \right]{{\left[ {\begin{array}{*{20}{c}} {{u_f}(t)}\\ {x(t)} \end{array}} \right]}^{\rm{T}}}} \right\}\\ {R_2} = \left[ {\begin{array}{*{20}{c}} \begin{array}{l} {U_f}\\ {X_f} \end{array} \end{array}} \right]E_f^{\rm{T}} = \bar E\left\{ {\left[ {\begin{array}{*{20}{c}} {{u_f}(t)}\\ {x(t)} \end{array}} \right]e_f^{\rm{T}}(t)} \right\} \end{array}$
若系统为可控可观系统,则$\Omega$为满秩矩阵
${\rm{rank}}(\Omega ) = {n_x} + f{n_u}$
(31) 由于新息为平稳零均值白噪声系列,可知$E_fE_f^{\rm T}$ $>$ $0$,据文献[11]定理2可知,式(31)中第2项的秩可通过如下计算得到:
$\text{rank}\left\{ \bar{E}\left\{ \left[ \begin{matrix} {{u}_{f}}(t) \\ x(t) \\ {{e}_{f}}(t) \\ \end{matrix} \right]{{\left[ \begin{matrix} {{u}_{f}}(t) \\ x(t) \\ {{e}_{f}}(t) \\ \end{matrix} \right]}^{\text{T}}} \right\} \right\}-f{{n}_{y}}$
(32) 进一步将式(4)带入式(32),可得:
$\begin{align} & \bar{E}\left\{ \left[ \begin{matrix} {{u}_{f}}(t) \\ x(t) \\ {{e}_{f}}(t) \\ \end{matrix} \right]{{\left[ \begin{matrix} {{u}_{f}}(t) \\ x(t) \\ {{e}_{f}}(t) \\ \end{matrix} \right]}^{\text{T}}} \right\}= \\ & \Upsilon \bar{E}\left\{ \left[ \begin{matrix} {{u}_{p}}(t) \\ {{u}_{f}}(t) \\ {{e}_{p}}(t) \\ {{e}_{f}}(t) \\ \end{matrix} \right]{{\left[ \begin{matrix} {{u}_{p}}(t) \\ {{u}_{f}}(t) \\ {{e}_{p}}(t) \\ {{e}_{f}}(t) \\ \end{matrix} \right]}^{\text{T}}} \right\}{{\Upsilon }^{\text{T}}} \\ \end{align}$
(33) 其中,$\Upsilon =\left[ \begin{matrix} 0 & {{I}_{p{{n}_{u}}}} & 0 & 0 \\ {{L}_{1}} & 0 & {{L}_{1}} & 0 \\ 0 & 0 & 0 & {{I}_{p{{n}_{y}}}} \\ \end{matrix} \right]$.
由于系统为可控可观测系统,采用文献[11]定理2可知$\Upsilon$为满秩矩阵.同时如果 式(33)第2项为正定满秩矩阵,根据文献[12]给定的秩条件可得:
$\begin{align} & \text{rank}\left\{ \overline{\text{E}}\left\{ \left[ \begin{matrix} {{u}_{f}}(t) \\ x(t) \\ {{e}_{f}}(t) \\ \end{matrix} \right]{{\left[ \begin{matrix} {{u}_{f}}(t) \\ x(t) \\ {{e}_{f}}(t) \\ \end{matrix} \right]}^{\text{T}}} \right\} \right\}= \\ & f({{n}_{u}}+{{n}_{y}})+{{n}_{x}} \\ \end{align}$
(34) 因此,矩阵(32)为满秩矩阵.
$\text{rank}\left( \left[ {{R}_{1}}-{{R}_{2}}{{({{E}_{f}}E_{f}^{\text{T}})}^{-1}}R_{2}^{\text{T}} \right] \right)=f{{n}_{u}}+{{n}_{x}}$
(35) 由式(30)、式(32)和式(35)可知\vskip0.1mm\noindent
$\text{rank}\left( \underset{N\to 0}{\mathop{\lim }}\,{{W}_{f}}\Pi _{{{{\hat{E}}}_{f}}}^{\bot }W_{f}^{\text{T}} \right)=f{{n}_{u}}+{{n}_{x}}$
(36) 则
$\begin{align} & \text{rank}\left( \underset{N\to 0}{\mathop{\lim }}\,{{W}_{f}}\Pi _{{{{\hat{E}}}_{f}}}^{\bot } \right)= \\ & \text{rank}\left( \underset{N\to 0}{\mathop{\lim }}\,{{W}_{f}}\Pi _{{{{\hat{E}}}_{f}}}^{\bot }W_{f}^{\text{T}} \right)=f{{n}_{u}}+{{n}_{x}} \\ \end{align}$
(37) 由以上秩条件可知,$\lim_{N \to \infty }W_f\Pi_{\hat{E}_f}^{\bot}$的非零特征向量个数为$f(n_u+n_y)-n_x-fn_u=fn_y-n_x$.因此$\lim_{N \to \infty }\hat{U}_1^{\bot}$ 是$\lim_{N \to \infty }W_f\Pi_{\hat{E}_f}^{\bot}$的零特征向量.则
$\underset{N\to 0}{\mathop{\lim }}\,\hat{U}_{1}^{\bot }{{W}_{f}}\Pi _{{{{\hat{E}}}_{f}}}^{\bot }=0$
(38) 由于${\rm rank}(\Gamma_f^{\bot})=fn_y-n_x$ [13],因此
$\text{rank}\left( [{{(\Gamma _{f}^{\bot })}^{\text{T}}}-{{(\Gamma _{f}^{\bot })}^{\text{T}}}{{H}_{f}}] \right)=f{{n}_{y}}-{{n}_{x}}$
(39) 同时,由于
$\underset{N\to 0}{\mathop{\lim }}\,\left[ {{(\Gamma _{f}^{\bot })}^{\text{T}}}-{{(\Gamma _{f}^{\bot })}^{\text{T}}}{{H}_{f}} \right]{{W}_{f}}\Pi _{{{{\hat{E}}}_{f}}}^{\bot }=0$
(40) 所以$[(\Gamma_f^{\bot})^{\textrm T}-(\Gamma_f^{\bot})^{\textrm T}H_f]$和$\lim_{N \to \infty }\hat{U}_1^{\bot}$满足
$\begin{align} & rowspace\left( [{{(\Gamma _{f}^{\bot })}^{\text{T}}}-{{(\Gamma _{f}^{\bot })}^{\text{T}}}{{H}_{f}}] \right)= \\ & rowspace\left( \underset{N\to 0}{\mathop{\lim }}\,\hat{U}_{1}^{\bot } \right) \\ \end{align}$
(41) 因此,$[(\hat{\Gamma}_f^{\bot})^{\textrm T} -(\hat{\Gamma}_f^{\bot})^{\textrm T}\hat{H}_f]=(\hat{U}_1^{\bot})^{\textrm T}$,由于$\hat{E}_f$是一致估计值,可知$(\hat{U}_1^{\bot})^{\textrm T}$为一致估计值,则$[(\hat{\Gamma}_f^{\bot})^{\textrm T}$ $-$ $(\hat{\Gamma}_f^{\bot})^{\textrm T}\hat{H}_f]$为一致估计值.
由于$[(\hat{\Gamma}_f^{\bot})^{\textrm T} -(\hat{\Gamma}_f^{\bot})^{\textrm T}\hat{H}_f]$为一致估计值,根据平移变换法求取系统矩阵的一致不变性,可知系统矩阵估计值也为一致估计值.
2.3 闭环一致条件合理性分析
第2.2节定理2的两个假设条件涉及未知新息信息,以上假设是否合理,将直接决定本文方法能否取得一致估计结果.
由于模型(1)可表示为模型(2),则新息的协方差可表示为$\bar{E}[e(t)e^{\textrm T}(t)]=CP^{\textrm T}C+R_3$ (具体表述可参考文献[10]第5章),其中P为半正定矩阵,对于系统噪声有$R_3=\bar{E}[v(t)v^{\textrm T}(t)]>0$,则$\bar{E}[e(t)e^{\textrm T}(t)]$ $>$ $0$,可知$E_fE_f^{\textrm T}>0$.
对于假设2,由于
$\begin{align} & \underset{N\to 0}{\mathop{\lim }}\,{{W}_{p}}(t,N+f)W_{p}^{\text{T}}(t,N+f)= \\ & \bar{E}\left\{ \left[ \begin{matrix} {{u}_{p}}(t) \\ {{y}_{p}}(t) \\ \end{matrix} \right]{{\left[ \begin{matrix} {{u}_{p}}(t) \\ {{y}_{p}}(t) \\ \end{matrix} \right]}^{\text{T}}} \right\} \\ \end{align}$
(42) 进一步,由式(2)可知
$\begin{align} & {{y}_{p}}(t)={{\Gamma }_{p}}[{{L}_{1}}u_{p}^{\text{T}}(t-p)+{{L}_{2}}e_{p}^{\text{T}}(t-p)]+ \\ & {{H}_{p}}{{u}_{p}}(t)+{{{\bar{G}}}_{p}}{{e}_{p}}(t) \\ \end{align}$
(43) 其中,扩展可观性矩阵分别表示为 ${L}_1=[{A}^{p-1}{B}$ $\cdots$ ${A}{B}$ ${B}]$,${L}_2=[{A}^{p-1}K$ $\cdots$ ${A}K$ $K]$,下三角形Toeplitz矩阵为
${{H}_{f}}=\left[ \begin{matrix} I & \cdots & \cdots & 0 \\ C & \cdots & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ C{{A}^{f-2}} & \cdots & C & I \\ \end{matrix} \right]$
将式(43)带入式(42),可得:
$\begin{align} & \bar{E}\left[ \begin{matrix} {{u}_{p}}(t) \\ {{y}_{p}}(t) \\ \end{matrix} \right]{{\left[ \begin{matrix} {{u}_{p}}(t) \\ {{y}_{p}}(t) \\ \end{matrix} \right]}^{\text{T}}}= \\ & \Xi \bar{E}\left[ \begin{matrix} {{u}_{p}}(t-p) \\ {{u}_{p}}(t) \\ {{e}_{p}}(t-p) \\ {{e}_{f}}(t) \\ \end{matrix} \right]{{\left[ \begin{matrix} {{u}_{p}}(t-p) \\ {{u}_{p}}(t) \\ {{e}_{p}}(t-p) \\ {{e}_{f}}(t) \\ \end{matrix} \right]}^{\text{T}}}{{\Xi }^{\text{T}}} \\ \end{align}$
(44) 其中,$\Xi=\left[\begin{array}{cccc} 0&I_{pn_u}&0&0\\Gamma_{P}L_1&H_p&\Gamma_{P}L_2&\bar{G}_p \end{array}\right]$.
若系统可控可观,由文献[11]定理2可知,${\rm rank}\left(\Xi\right) =p(n_u+n_y)$,$\Xi$是满秩矩阵,则式(42)正定的充分必要条件为式(44)中第2项正定,由于变量p和f可取任何值,则假设2成立的充分必要条件为$\lim_{N \to \infty }W_{p+f}W_{p+f}^{\textrm T}>0$.可通过验证$\lim_{N \to \infty }W_{p+f}W_{p+f}^{\textrm T}>0$来确定假设2是否成立.即若$\lim_{N \to \infty }W_{p+f}W_{p+f}^{\textrm T}>0$,则假设2成立.
3. 仿真研究
考虑文献[6]中研究的闭环系统
$y(t)-0.9y(t-1)=u(t-1)+e(t)+0.9e(t-1)$
(45) 其中,反馈控制结构为$u(t)=- 0.6y(t)+r(t)$.过程噪声设置为方差为0.2的白噪声序列.
对系统设定输入激励$r(k)$为单位方差白噪声序列和相关序列两种情况进行研究. 其中相关序列设为$r(t)=(1+0.8q^{-1}+0.6q^{-2})r_0(t)$,$r_0(t)$为单位方差白噪声序列. 过去和未来水平数均设置为10.在不同数据长度$N\in[200,8 000]$的情况下进行1 000次Monte Carlo仿真. 本文方法的辨识结果与SIMPCAwc算法 [6]进行比较,当$r(t)$为单位方差白噪声系列时,系统极点(真实值为0.9)的估计平均值如图 1所示,系统极点的估计标准方差如图 2 所示.当系统设定输入$r(t)$为相关系列时,系统极点的平均值如图 3所示,系统极点的估计标准方差如图 4所示.
图 1 系统设定输入激励为白噪声的系统极点估计平均值Fig. 1 Mean value of estimated poles for white noise setpoint excitation
图 2 系统设定输入激励为白噪声的系统极点估计标准方差Fig. 2 Standard deviation of estimated poles for white noise setpoint excitation
图 3 系统设定输入激励为相关序列的系统极点估计平均值Fig. 3 Mean value of estimated poles for correlated quasi-stationary setpoint excitation
图 4 系统设定输入激励为相关序列的系统极点估计标准方差Fig. 4 Standard deviation of estimated poles for correlated quasi-stationary setpoint excitation从以上结果可以看出,当系统设定输入为单位方差白噪声系列时,SIMPCAwc算法得出有偏估计结果;只有当系统设定输入激励为相关序列时,SIMPCAwc算法才能保证无偏估计结果.本文CSIMIEOP算法对系统设定输入激励为无关序列和相关序列时均可得到一致无偏结果,并且估计精度优于SIMPCAwc算法.
4. 结论
本文提出一种基于新息估计和正交投影的闭环子空间辨识算法,对系统设定点输入激励为白噪声无关序列和相关序列的情况均可得到一致无偏估计结果,并且相对于近期有关文献给出的方法如SIMPCAwc算法 [6],能进一步提高辨识精度.同时,严格分析和证明了本文算法保证一致估计的条件.最后通过仿真实例验证了本文方法的有效性和优越性.
-
[1] Flusser J, Suk T, Zitova B. Moments and Moment Invariants in Pattern Recognition. UK: John Wiley Sons, 2009. [2] [2] Teague M R. Image analysis via the general theory of moments. Journal of the Optical Society of America, 1980, 70(8): 920-930 [3] [3] Teh C H, Chin R T. On image analysis by the methods of moments. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1988, 10(4): 496-513 [4] [4] Sheng Y L, Shen L X. Orthogonal Fourier-Mellin moments for invariant pattern recognition. Journal of the Optical Society of America A, 1994, 11(6): 1748-1757 [5] [5] Khotanzad A, Hong Y H. Invariant image recognition by Zernike moments. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1990, 12(5): 489-497 [6] [6] Kim H K, Kim J D, Sim D G, Oh D I. A modified Zernike moment shape descriptor invariant to translation, rotation and scale for similarity-based image retrieval. In: Proceedings of the 2000 IEEE International Conference on Multimedia and Expo. New York, USA: IEEE, 2000. 307-310 [7] [7] Novotni M, Klein R. 3D Zernike descriptors for content based shape retrieval. In: Proceedings of the 8th ACM Symposium on Solid Modeling and Applications. New York, USA: ACM, 2003. 216-225 [8] [8] Gope C, Kehtarnavaz N, Hillman G. Zernike moment invariants based photo-identification using fisher discriminant model. In: Proceedings of the 26th Annual International Conference of the IEEE Engineering in Medicine and Biology Society. San Francisco, USA: IEEE, 2004. 1455-1458 [9] [9] Li S, Lee M C, Pun C M. Complex Zernike moments features for shape-based image retrieval. IEEE Transactions on Systems, Man, and Cybernetics, Part A: Systems and Humans, 2009, 39(1): 227-237 [10] Chen Z, Sun S K. A Zernike moment phase-based descriptor for local image representation and matching. IEEE Transactions on Image Processing, 2010, 19(1): 205-219 [11] Liao S X, Pawlak M. On image analysis by moments. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1996, 18(3): 254-266 [12] Liao S X, Pawlak M. On the accuracy of Zernike moments for image analysis. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1998, 20(12): 1358-1364 [13] Xin Y Q, Pawlak M, Liao S. Accurate computation of Zernike moments in polar coordinate. IEEE Transactions on Image Processing, 2007, 16(2): 581-587 [14] Kotoulas L, Andreadis I. Accurate calculation of image moments. IEEE Transactions on Image Processing, 2007, 16(8): 2028-2037 [15] Lin H B, Si J, Abousleman P G. Orthogonal rotation-invariant moments for digital image processing. IEEE Transactions on Image Processing, 2008, 17(3): 272-282 [16] Ma H, Qi D X, Song R X, Wang T J. The complete orthogonal v-system and its applications. Communications on Pure and Applied Analysis, 2007, 6(3): 853-871 [17] Huang C, Yang L H, Qi D X. A new class of multi-wavelet bases: V-system. Acta Mathematica Sinica, English Series, 2012, 28(1): 105-120 [18] Liang Yan-Yan, Song Rui-Xia, Wang Xiao-Chun, Qi Dong-Xu. Complete orthogonal V-system and its application in geometrical information reconstruction. Journal of Computer-Aided Design Computer Graphics, 2007, 19(7): 871-875, 883 (梁延研, 宋瑞霞, 王小春, 齐东旭. 完备正交V-系统及其在几何信息重构中的应用. 计算机辅助设计与图形学学报, 2007, 19(7): 871-875, 883) [19] Li Jian, Song Rui-Xia, Ye Meng-Jie, Liang Yan-Yan, Qi Dong-Xu. Orthogonal reconstruction of 3D model based on V-System over triangular domain. Chinese Journal of Computers, 2009, 32(2): 193-202(李坚, 宋瑞霞, 叶梦杰, 梁延研, 齐东旭. 基于三角域上V-系统的三维几何模型的正交重构. 计算机学报, 2009, 32(2): 193-202) [20] Song Rui-Xia, Chen Xi, Sun Hong-Lei, Yao Dong-Xing, Xue Guan-Chen. A novel algorithm of classification and retrieval for shape group. Journal of Computer-Aided Design Computer Graphics, 2011, 23(12): 1981-1986 (宋瑞霞, 陈曦, 孙红磊, 姚东星, 薛冠辰. 形状群组的分类和检索算法. 计算机辅助设计与图形学学报, 2011, 23(12): 1981-1986) [21] Qi Dong-Xu, Song Rui-Xia, Li Jian. Discontinuous Orthogonal Functions. Beijing: Scientific Press, 2011. (齐东旭, 宋瑞霞, 李坚. 非连续正交函数. 北京: 科学出版社, 2011.) 期刊类型引用(12)
1. 王亚军，白翱，张博，郭超. 基于TOGAF的工艺知识全生命周期管理架构研究. 知识管理论坛. 2024(06): 519-532 . 
百度学术2. 刘拥民，罗皓懿，谢铁强. 基于XGBoost-ARIMA方法的PM_(2.5)质量浓度预测模型的研究及应用. 安全与环境学报. 2023(01): 211-221 . 百度学术3. 张旭，张亮，金博，张红哲. 基于不确定性的多元时间序列分类算法研究. 自动化学报. 2023(04): 790-804 . 本站查看4. 缪燕子，王志铭，李守军，代伟. 基于背景值和结构相容性改进的多维灰色预测模型. 自动化学报. 2022(04): 1079-1090 . 本站查看5. 徐任超，阎威武，王国良，杨健程，张曦. 基于周期性建模的时间序列预测方法及电价预测研究. 自动化学报. 2020(06): 1136-1144 . 本站查看6. 李晓理，张博，杨旭. 基于图像混合核的列生成PM_(2.5)预测. 工程科学学报. 2020(07): 922-929 . 百度学术7. 杨博帆，张琳，张搏，宋亚飞，丁尔启. 动态多模型指数平滑法融合的在线预测方法. 系统工程与电子技术. 2020(09): 2013-2021 . 百度学术8. 乔俊飞，丁海旭，李文静. 基于WTFMC算法的递归模糊神经网络结构设计. 自动化学报. 2020(11): 2367-2378 . 本站查看9. 郝广涛，林清华，李晓梅. 超短期负荷预测中指数平滑法平滑系数的确定方法. 莆田学院学报. 2020(05): 80-86 . 百度学术10. 杨亚莉，李智伟，钟卫军. 基于二向注意力循环神经网络的PM_(2.5)浓度预测. 空军工程大学学报(自然科学版). 2020(06): 101-106 . 百度学术11. 慕凯，张祥，余士龙，李俊，李宇，戴笑俊. 滑动加权马尔科夫模型在降水量预测中的应用. 气象水文海洋仪器. 2019(04): 34-36 . 百度学术12. 周强，肖强宏，王浩然，高乐乐. 基于BEADS-ESMC组合算法的三相光伏并网逆变柜触点红外温度预测方法. 变频器世界. 2018(11): 72-78 . 百度学术其他类型引用(16)
-
下载:
计量
- 文章访问数: 1664
- HTML全文浏览量: 68
- PDF下载量: 847
- 被引次数: 28
下载:
