系统辨识中不同噪声情况的研究

D0I:10.13374/i.issn1001-053x.1990.03.010 北京科技大学学报 第12卷第3期 Vo1.12No.3 1990年5月 Journal of University of Science and Technology Beijing May 1990 系统辨识中不同噪声情况的研究 沃利拉M·刘宏才· 摘要：对实际过程中的量测噪声和系统噪出进行了分析，非做了仿直研究，仿真结果 了解到它们的差别，这对利用系统辨识的方法来建模和参数估计是有益的，还提出一种交替 矩阵方法，它不仅可用来进行参数估计还可用来求解病态矩阵问题。 关键词：系统辨识，参数估计，脉冲响应函数 Different Noises in System Identification Olehla Miroslov Liu Hongcai ABSTRACT:The measurement noise and system noise in the practical process are studied.The analysis and the simulation results show that the impulse responses of them are different.It is profit for modelling and system parameter estimation.The Rotation Matrix method is also proposed.It not only estimates system parameters but also finds the solution of ill-matrix. KEY WORDS:system identification,parameter estimation,impulse response 在利用系统辨识的方法对实际过程建模与参数估计时，常以一个等效的噪声来代表量测 噪声、系统噪声和模型误差等。本文对量测噪声和系统噪声进行了分析和仿真研究，仿真结 果显示出两者之间的差别，认清这一点对于利用系统辨识的方法来建模是有益处的。此外提 出一种参数估计的算法，也可用来对病态矩阵求解。 1问题陈述 考虑过程模型为 A(z-I)y(k)=B(z-1)u(k-d)+2()+f(k) (1) 1989-12-23收稿 ·为高级进修生，捷克斯洛伐克里培瑞兹技术大学博士、副教授 …自动化系(Dcpt,of Autom,and Control Eng,. 257

、 、 第径 卷第 期 年 弓 月 北 京 科 技 大 学 学 报 系统辨识 中不 同噪声情况的研究 沃利拉 ’ 刘宏 才 ’ ‘ 、 、 捅 要 对 实际 过 柳扣的 量测 噪 声和 系统 噪 声进 行 了分 析 ， 井 做 了 仿 真 研究 ， 仿 真结 果 了解 到它 们 的差 别 ， 这对 利用 系统 辨识 的方法 来建模和 参数 估计是 有益 的 ， 还 提 出一 种交 替 矩 阵 方法 ， 它 不 仅可 用 来进行 参 数 估计 还可用 来求 解病态 矩 阵问题 。 关键词 系 统 辨识 ， 参数 佑计 ， 脉 冲 响应 函数 几 、 ，、 、 ， ， · 现 · ， 了 一 ， ， 、 在利 用 系统辨识的 方 法对 实际过 程建模与参数估计时 ， 常以 一个等效 的噪 声来代表量测 噪 声 、 系统噪声和 模型 误差等 。 本文对 量测噪声和 系统噪声进行 了分析和仿 真研究 ， 仿 真结 果显示 出两 者之 间的差 别 ， 认 清这 一点对于利 用 系统辨识 的方 法 来建 模是有益 处的 。 此 外提 出一种参数估计 的算法 ， 也 可用来对 病态 矩阵求解 。 问 题 陈 述 考虑过 程模 型为 一 ‘ 夕 一 ’ 一 一 一 收 稿 为高 级进 修 生 ， 捷 克 斯洛伐克里 培瑞兹技 术大学 博 士 、 副教授 自动化系 DOI ：10．13374／j ．issn1001－053x．1990．03．010

其中 y(k)和u()为可测的输出变量和控制变量；d为过程的时滞时间：(k)为不可测量的扰 动变量，{2(k)}为独立同分布的随机序列，且满足E{z(k)}=0,E(2(k)2T(k)}=σ2， f(k)为漂移噪声， 2-1为单位延迟算子， A(21)、B(21)分别为输出变量和控制变量的加权多项式，即 A(2-1)=1+a121+…+a.2- B(-1)=b。+b12-1+…+bm2-m 1.1考患噪声在输出端的情况（设f(k)=0) 由图1可写出过程模型为 y()+∑1y(质-i)=∑bu(质-i)+2()+∑12(k-i) (2) i=1 i=0 i1 2) +2k) u(k) + B(2 xky(k) )B(3 A(z-) +Oyk) AE B(( z+ b 图1有物出噪声的过程模型 用2有输人噪声的过程模型 Fig.1 Process model with output Fig.2 Process model with input noisc noise 其中 x()为过程的理想输出变量， 2(k)为过程的输出噪声。 式(2)也可写成 A(z-1)y(k)=B(=-1)u(k)+A(-1)z(k) (3) 其中 符号说明同前。 这种模型参数的估计值，可采用最小二乘法或其他方法获得。 1.2考忠噪声在输入端的情况（图2a) 由于是线性系统，故图2a可画成图2b。其模型为 A2-)y(®)=B2-)u()+B(2-12-）z(®) A(2-1) (4) 这时最优滤波器为(1/B(z-1))。 由式(3)和式(4)看出，噪声在输出端和输入端时的最优滤被器是不同的。当噪声未知 258

其 中 ， 幻 和 。 幻 为可测 的输 出变量和控 制 变 量， 为过 程的时 滞时 间 句 为 不可侧 量的扰 动 变 量 ， 幻 为独 立 同分 布的随机序 列 ， 且满足 笼 ， 以 幻 。 ， 为 漂移噪 声 一 ‘ 为单位延 迟算子 一 ’ 、 二 ’ 分别 为输 出变量和 控 制变量的加权 多项式 ， 即 一 ’ 一之 一 ‘ … 二 二 一 “ 一 。 之 一 ’ … 二 之 一 口 考虑嗓 声在翰出端 的情况 设 二 由图 可写 出过 程模型 为 ， “ 名 ， ‘ 一 ‘ 名 。 一 ‘ “ 名 掩一 ‘ 吻羽 卫孵 图 有粉 出嗓声的 过程模 型 图 有 物入嗓 声的 过程摸型 其 中 为过程的理 想输出变量， 以 为过 程的输出噪 声 。 式 也可写 成 一 ’ 二 一 一 二 掩 其 中 符号说明 同 前 。 这种模型 参 数的估 计值 ， 可采用 最小 二乘法或其他方法获得 。 。 考虑 嗓声在 入端 的情 况 田 由于 是线 性 系统 ， 故图 可 画成图 。 其模型为 一 二 二 一 吞 一 之 一 之 一 ‘ 庵 口 这 时 最优滤 波器 为 “ 一 ’ 。 由式 和 式 看出 ， 噪声在输出端和输入 端时的最优滤波器是不同的 。 当噪声未知

时，噪声滤波器不是采用文献〔1)中的1/A(2-1),而是用1/C(=1)或D(z1)或C(z-1)1 D(z-1)。如果采用1/C(2-1)的话，其模型为熟知的形式 A(z-1)y()=B(z-1)u(k)+C(2-i)z(k) (5) 其中 C(2-1)为2~的多项式（滑动平均） C(21)=1+C12-1+…+Cn2-p 同样，D(21)也可写成：~1的多项式（自回归）。当然最好是用C(e-)/D(e-1)形式的滤波 器，如图3所示。 于是过程模型为 A(z-1)y(k)=B(-1)u(k)+e() (6) 其中 e()=4(2)C-)z®)△H(2-1)z(k) D(z-1) 4 ·R) g() -y() 2 C(2) z1 2() n(k) L、 (3) B(2-) A(2)+ -y火R) )+2(k) 图3有输出噪声的过程模型 图4有输人量测噪声的过程模型 Fig,3 Process model with output noise Fig.4 Process model with input measu- rement noisc n(k)为随机噪声，在大多数情形，n()都可表示具有有理谐密度的平稳随机序列。它总 可表示为一个白噪声序列{()}为输人的线性系统输出。 D(s-1)=d12-1+…+dg2-g 以上讨论的是参数模型，当然也可用非参数模型来表示，如果图1用脉冲响应函数g() 表示的话，那么模型写为 m-1 y(i)=∑9(G》u(i-i)+z(i) (7) j=0 其中 i=1,2,…,W;N>m 1.3考虑量测噪声在输入端的情况（图4） 这时过程模型为 259

时 ， 噪声滤 波器 不是采 用文献 〔 〕 中的 二 ’ ， 而是 用 厂 二 ’ 或 刀 二 ’ 或 二 一 ， 犷 ’ 。 如果采 用 犷 ’ 的 话 ， 其 模型为熟 知的形式 之 一 ’ 夕 二 之 一 ’ 之 一 ’ 及 其 中 一 ’ 为 一 ’ 的 多项式 滑 动平均 犷 ’ 之 一 ’ … 艺 一 ， 同样 ， 二 ’ 也 可写 成 一 ’ 的 多项式 自回 归 。 当然 最 好是 用 一 ’ 二 ‘ 形式的滤波 器 ， 如图 所示 。 于 是过程模型 为 之 一 二 一 掩 勺一 其 中 。 月 之 之 一 二 叠 之 一 ’ 二 才 夕 花 夕 瓦 城九卜之 花 图 有输出噪 声的过 程模型 通 图 呈 有输入量 测噪 声的 过 程模型 爪 王 爪 城句 为随机噪 声 ， 在大多数情形 ， 。 句 都可表示具 有有理 谱密度的 平稳随机序列 。 它总 可表示为 一个 白噪 声序列 以 幻 为输人 的线 性系统输 出 。 一 ’ 声 一 … 十 沪 一 以上 讨论的是 参数模 型 ， 当然 也 可用非参数模型 来表示 ， 如果图 用脉 冲响应 函数 表 示 的话 ， 那 么模型写 为 一 ， ‘ 云 。 ‘ 一 ‘ 二 其 中 ， ， … ， 》 州 考虑 侧 噪 声在输 入端 的情 况 图 这时过程模型 为

每-1 y(i)=∑g(j)〔u(i-)+z(i-i)〕 (8) j=0 如果噪声：()为相关噪声，那么就不能得到一致性估计。 1,4考虑过程中存在零漂f()的情况 由于f()是趋势项，它可表示为： f)=∑w(j)i (9) i=0 其中 j=1,2,…,N 通过实验就可按式(9)来拟合趋势项f(k)。 1.5交替矩阵(Rotation Matrix)法 上述参数模型有时可简写成 y=60+e (10) 其中 y=〔y(1)y(2)…y(N)]T多 pT(1) 中= pT(2) T(N)) -y(0) -y(-1)…-y(1-n),4(0)u(-1)…u(1-m) -y(1) -y(0)… -y(2-n),4(1)u(0)…u(2-m) :;: -y(N-1)-y('-2)…-y(N-n),u(N-1)u(N-2)…u(N-M)5 0T=〔a:…0.b1…baJ eT=〔e(1)e(2)e(N)〕。 式(10)的未知参数可用多种参数估计方法求得。如用最小二乘法的话，参数估计值为 8=(中T中)-'中Ty=(φ*)-1Ty (11) 要使上式有解，必须中T阵为非奇异。如果中T中阵为奇异，则称中？中为病态矩阵。这时 可用Householder?变换的方法求解未知参数，但这方法运算较麻烦。这里介绍一种简单的方 法一交替矩阵T法，其主要思想是将φ·阵经过多次选得的交替矩阵左乘后，变为一个三角阵。 而T?T=【(单位阵)，说明参数估计的误差平方和是不变的。 T阵的一般形式为 260

伍 一 ， ‘ 习 。 〔 ‘ 一 ‘ 一 〕 如 果噪 声 旬 为 相关噪 声 ， 那 么就 不能 得到 一 致 性估计 。 考虑过 程 中存在零 漂 的情 况 由干厂 是趋势项 ， 它 可表示 为 护 吞 二 习 却 二 其 中 二 ， ， … ， 通过 实验就可 按式 来拟 合 趋势项 。 交替 矩阵 法 上述参数模 型有时可简写 成 价 其中 二 〔夕 一 〕 甲 伊 卯 一 一 一 一 … 一 … 一 一 ” ， 赵 一 又 一 ， “ 一 · ” “ … ‘ 一 优 ‘ 一 用 一 少 刃 一 一 夕 入 尹 一 · 一 少 一 ， “ 一 一 一 一 、‘ ‘… 苗︸一一 口 二 〔 口 一 。 … 二 〕 〔 … 〕 。 式 的未知参数可 用多种 参数估计方法求 得 。 如 用 最小 二乘 法的 话 ， 参数估计值 为 口 二 价 沪 一 ’ 功 夕 二 功 一 ‘ 价 要使 上 式 有解 ， 必须 沪 功阵为非 奇异 。 如果 功 功阵为 奇异 ， 则称 价 必为病 态 矩 阵 。 这 时 可 用 变换 的方法求解未知参数 ， 但 这 方 法运算较麻 烦 。 这 里介绍 一种 简 单 的方 法一交替矩 阵 法 ， 其主 要 思想是 将功 阵经 过多次 选 得的 交替矩 阵左乘后 ， 变为一 个三角阵 。 而 单 位阵 ， 说 明参数 估 计的误 差平 方和 是不变的 。 阵的 一般形 式为

c 其中c2+s2=1 用矩炬阵T左乘中·阵后，则中·阵中的相应的元素值会改变。如果选定c和s为某值，可以使 T中“矩阵的某一元素值为零，然后再选定c1和51又可使T中“矩阵的另一元素值为零，如此重 复就可使“阵的下三角全为零元素。这样6·阵变为上三角阵，就可求？的逆了，举个例子 说明具体做法。 例： 设 （p11P12P13） 0*= P21P22P23 P31P32pa3) 现在希望p21,p31及P32为零，做法分三步。 第一步 选择T21,并左乘·，得 C12-S12 0 T21= S12C12 0 001)； 「c12p11-912p21=pc12p12-512p12=p2 T21= s12p11+c12p21=p1S12p12+c12p23=p2) 931 P32 C12p13-S12p13=pg s12p13+C12p23=p} =功*1) P33 设c12+s2=1 S12p11+C12p21=p}=0 则得 -P21 P11 s12=√91+p31； c12=/p1+p1 第二步 选择Tg1为 261

… 尸 、 、 … 二 一 。 点 … 了 ，土 了 … ‘ 炜、 其 中 ’ 用矩阵 左 乘必 ‘ 阵后 ， 则价 ‘ 阵中的 相 应的元 素值会改 变 。 如果选定。 和 为某值 ， 可 以使 少 矩阵的某 一元 素 值 为零 ， 然后再选 定 。 ，和 ， 又可使 势 矩阵的 另 一元 素值为 零 ， 如此重 复就可使功 阵的 下 三 角全为零元素 。 这 样娇 阵变为上 三 角阵 ， 就可求价 势的逆 了 ， 举 个例子 说明具体做 法 。 例 设 功 二 切 “ 甲 ‘ ’ 沪 ’ … 尹 ， ‘ 甲， ， ， 贬 甲 甲 甲 少… 现在希望 尹 ， 沪 及 卯 为零 ， 做 法分三步 。 第一步 选择 ， 并左乘 功 ’ ， 得 、 阮 一 厂 一 … 。 ， … 甲 一 甲 甲 ， 甲 尹 、 ， 二 尹 委、 ， 尹 一 ， 叨 沪 飞 ， 中 甲 尹 飞， 甲 ， 甲 ， 。 一 。 甲 ， 尹 犷孟 ’ 、 ， 毋 十 卯 二 尹 飞 ， 甲 。 ， ‘ 一 设 子 二 尹 ， ， 切 甲 毛 则得 “ ， ， 二 护势 切 十 甲 呈 “ ，， 二 护乒乳下币 第 二步 选择 。 为