应用m序列及逆重复m序列辨识非线性动态系统

D0I:10.13374/j.issm1001-053x.1981.02.016 北京钢铁学院学报 1981年第2期 应用m序列及逆重复m序列辨识 非线性动态系统 电工教研室幸白男钟延焖 摘 要 本文提出了辨识非线性动态系统的两个方法。非线性系统近似由线性动态系统 和非线性部件组成。在所提出的二个方法中，序列和逆重复m序列用于输入信号 或相关信号，这些方法可用于估算高阶多项式非线性部件的全部系数，并有足钩多 的脉冲响应栅点值。此外，所提出的方法还适于存在缓慢随机漂移干扰时的在线辨 识。由辨识试验的数字模拟可以证实所提方法的正确性。 引 言 本文提出了应用m序列及逆重复m序列辨识非线性动态系统的二个方法。非线性系统近 似由线性动态系统和单值非线性特性部件组成，即所谓维纳模型(W-模型)或Ham mer- stein模型(H模型)。方法一，利用逆重复m序列作被测系统的输入试验信号。对于W- 模型的辨识分二步，第一步，输入不同电平的逆重复序列信号，计算系统输入与输出的互 相关函数，解线性方程组后得到非线性系数a1,脉冲响应函数g(τ)和一组中间函数，第二 步，利用电平相同但偏置不同的逆重复序列信号作输入信号，计算互相关函数并解一组简 单的非线性方程组，最后得到全部非线性系数。对于H-模型的辨识，与上述第二步大致相 同。在上述方法一中，计算互相关函数时都采用了二项式加权平均，所以辨识结果能消除缓 慢漂移的影响。 方法二，输入信号仍采用逆重复m序列信号，不同的是利用与逆重复m序列同步的m序 列信号与被测系统的输出信号进行相关计算。 上述二方法都可用于估算系统中非线性部件的高阶多项式之系数。 辨识试验的数字模拟结果证实了所提出的方法的正确性。 方法 一 1.wiener模型辨识 在图一所示维纳模型中，AL(t)是电平为A的逆重复m序列。L是线性动态系统，NL ◆本文1980年8月4日收到。 142

北 京 铜 铁 学 院 学 报 年 第 期 应用 序列及逆重复 序列辨识 非 线 性 动 态 系 统 ’ 电工 教研 室 李白男 钟延 炯 摘 要 本文提 出了辨识 非线性动 态系统 的两个 方法 。 非线性系统 近 似 由线性动态系统 和 非线性部件组 成 。 在 所提 出的二 个 方法 中 ， 序 列和逆重 复 序 列用 于 输入 信 号 或相关信 号 ， 这 些 方法 可 用 于估 算高阶 多项 式非线性部件 的全 部 系数 ， 并有足 豹 多 的脉 冲响应栅 点值 。 此 外 ， 所提 出的方法还 适 于存在 缓 慢随机漂 移干扰 时的在线辨 识 。 由辨识 试 验 的数 字模 拟可 以证 实所提 方法 的正 确性 。 引 言 本 文 提出 了应用 序 列 及逆重 复 序列 辨识 非线性 动态 系统 的二个方法 。 非线性 系统近 似 由线性动态 系统 和 单值 非线性特性 部件组 成 ， 即所 谓维 纳模 型 一模 型 或 模 型 模 型 。 方法一 ， 利用 逆重 复 序列 作被测 系统 的输入 试验信号 。 对于 模 型 的辨识分二 步 ， 第一 步 ， 输入不 同 电平 的逆重 复 序 列 信号 ， 计 算系统输入 与输 出的互 相 关函数 ， 解线性方程组后得 到非线性 系数 ， 脉 冲响应 函数 幼 和一组 中间函 数 ， 第二 步 ， 利用 电平相同 但偏 置不同 的逆重 复 序列 信号作输入 信号 ， 计 算互相 关函数 并解一组简 单的非线性方程 组 ， 最 后得 到全部非 线性 系数 。 对 于 一模 型 的辨 识 ， 与上述 第二 步大致 相 同 。 在上述 方法一 中 ， 计 算互相 关 函数 时都 采用 了二项 式加 权 平均 ， 所 以 辨 识结果 能消除缓 慢 漂移 的影 响 。 方法二 ， 输入信号仍 采用逆重 复 序列 信号 ， 不同 的是利 用 与逆重复 序列 同步的 序 列 信号与被测 系统 的输 出信号进行相 关计算 。 上述 二方法都可用 于估 算系统 中非线性部件 的高阶多项 式之 系数 。 辨 识 试验的数 字模 拟结果证 实了所提 出的方 法 的正 确性 。 方 法 一 棋型辨识 在图一所示维 纳 模 型 中 ， 是 电平为 的逆 重 复 序 列 。 是线性 动态 系统 ， 本文 年 月 日 收到 。 DOI ：10．13374／j ．issn1001－053x．1981．02．016

是有2n阶多项式的非线性特性，即 n(t》+d(t) A1(t) y(t) z(t) 图一W-模型 2n y=∑au (1) !=1 n(t)是白噪声，d(I)是多项式型缓慢随机漂移。r是偶数。 d(t)= dt 2) 1-0 辨识上述非线性动态系统方法分成两步 第一步： 输入信号是AL(t),系统响应是y(t),有 y(t)=yi(t)+y2(t) (3) 其中： n y(I)= a=A1-6…f0gs-gs1-L(t-5 1=1 L(t-s21-1)dsi.ds21-1 (4) y(=∑A66 ()(-s L(t-s2i)ds1…ds21 (5) AL(t)与z(t)在0至(r+2)T区间权重平均互相关函数可以写成为： r+1(I+1)T $2(t)=1 (6) r+1 (∑K,)T ∑，AK:ztL1-)at 然而： 中Lz(t)=φLy1(T)+中Ly2(t)+中Ld(T)+中Ln(T) (7) 其中， 中Ln(T)=0 (8) +1 (r+1) (∑K)r 143

是有 阶多项 式 的非线性特性 ， 即 图一 一模 型 ， 艺 是 白噪声 ， 是 多项 式 型缓慢 随 机漂移 。 是 偶数 。 幻 。， ‘ ， 艺 ， 辨 识 上述 非线性动 态 系统 方法分 成 两 步 第一 步 输 入 信 号是 ， 系统 响应 是 ， 有 其 中 艺一 二 一 了 … 一 歹 ” ’ 卜 一 … 一 一 … ， ， ‘ 乏一 了 … 了 … 卜 … 一 一 … 与 在 。 至 区 间权 重 平均互相 关函数可 以 写 成为 几 击 ， ， 二 一 皿－ 兮 一 丫 艺 然而 小 小 ， ， 小 ， 小 中 。 丫 其 中 小 小 了 乏 ’ ‘ ’ 了 “ ， ‘ ’ 一 ‘ ’ “ ‘ 乏 ”

K!若按二项式系数选择，可使 φL4(T)=0 注意到 y(t)=-y:(t+T) (11) y2(t)=y2(t+T) (12) 我们有 中Lr2(T)=0 (13) (1+1)T ,()= ∑Ay:(t)L(t-t)dt (14) 最后有 中Lz(t)=中Ly,(t) (15) 代入(4)和(14)到(15)，(15)式可写成下列形式 中L:(t)=∑a-iA中L-(T) (16) 1=1 其中： 4.-(r)=00g(Sg(5,g(521-L (t-51,t-s2,…t-s21-1)ds1…ds21-1 (17) 取A=A,相应有中L(t)=中Lz1(t),i=1,2,…n,这时(16)可写成矩阵形式 a中L,=Aa中Ln (18) 或 a中L=A-中z (19) 从(19)式，我们可找到一组中间函数F,(τ) F,(t)=a21-:中Ln(21-1)(t)i=1,2,…n (20) 当i=1 F(t)=aisg(T) (21) 其中S是AL(t)的自相关函数面积。又假设 · 6g)dA=-1 (22) 从(21)式得到 a-号r(e)dr (23) g(t)=F1(T)/aIs (24) 第二步 输入信号为A。L(t)+B,B是一常偏置量，系统响应是y(t),我们假定 y(t)=y,(t)+y2(t) (25) 144

，若按二项 式 系数选择 ， 可使 小 ， 注 意到 ， 一 我们有 勺 了 ‘汗 小 ， ， 小 ， 丫 级一 一 ” 、 ‘ ， 一 丫 最后 有 小 ， 小 ， 代入 和 到 ， 式可写 成下列形式 小 ‘ ，二 艺 ，一 ” 小 。 卜 ‘ 一 小 。 ‘ 卜 ， · 卜 厂 … 歹了 … 卜 ，小 · · … 丫 一 ， ， ， 丫 一 ， … 一 一 ， ，… 一 相应 有小 ， 丫 小 。 二 ， ， ， 一 ， 这 时 可 写 成矩 阵形 式 小 ， 小 或 从 式 “ 呀， 一 ’ 小 我们可 找 到一组 中间函数 一 小 ， ’ 一 ’ ， ， … 当 ， ， 丫 其 中 是 的 自相 关函数面积 。 又假设 。 入 入 一 从 式得 到 。 ， 、 ， 一 认 丫 第二 步 输入 信号为 。 ， 是一常偏 置 量 ， 系统 响应 是 ， 我们假定 一

其中： y(t)A.fog(s)L(t-s)ds-2a:BA.fog(s)L(t-s)ds+ A.(s:)g(s.)g(5a)L(t-5:)L(t-s2)L(t-5a)dsids:dsa 十… (26) y:(t)--aB+a:A.-(5(52)L(I-5)L(t-5.)dsidsa+ aB-3m,A.BB(S.)E(Sa)L(t-S)L(t-S.)dS:dS. 十… (27) 从以上两式可以得出： y(t)=-y:(t+T) (28) y2(t)=y2(t+T) (29) A。L(t)与Z(t)在0至(r+2)T区间具有权重平均的互相关函数仍然得到 e)=∫rAy,(L-dt (30) 代入(26)到(30)，在0至T区间积分以上公式的两边，给出： n-I -M-CB ·(31) 其中 C1=2aS-4A,F:e)dt-6AF,(t)dr+… C=-a,5+10AF,tat+21AFedt+ C,=4a,5-20;AF.(r)dt-… C4=-5aS+35 A.:()+… (32) a3 Cs=6aaS-… C。=-7a,S+… Ca-1=(-1)"anS M=A,fF()dt+A.fF:()dt+A.F(d 当B=B。=0,相应有中Lz（t)=中Lz。(t),M可从(31)式得到，有： M=中or)dr (33) 145

其中 。 。 户 一 、 、 一 。 少 八 一 、 ， 一 ， 一 一 。 一 产 · … ’ 一 。 · ， 一 一 一 。 。 么 一 一 十 。 。 。 二 从 以 上 两式 可 以得出 ， 一 二 。 与 在 。 至 十 区 间具有权重平均 的互相 关函数仍然得到 小 · · 书 ‘ ’ ‘ ” ’ 。 卜 · 代入 到 ， 在 。 至 区 间积 分 以 上公 式的 两边 ， 给 出 小一 · 一 艺 其 中 、，、于吸 一 吐竺丝一 。 歹 卜 红 。 丫 … ‘ ， 一 一 ‘ 竺一 。 七 一 。 舀 五 。 丫 。 二 一 一 丫 一 … ‘ 一 。 里 互 。 。 厂 “ ， … 。 声 一 … 一 … 一 ， 一 “ 。 。 ’ · · · · 。 。 · · … 当 。 ， 相 应 有 小 小 · 。 ， 可 从 式得 到 ， 有 。 中 ‘ 气 少

当B=B:,相应有中Lz(t)=中Lz1(x),i=1,2,…n-1,(31)式可写成为矩阵形式，有 L:(d+-M)=BC (34) 解以上矩阵方程，我们得到 C=B-'(中a(r)dr-M) (35) 代入中间函数F:(t)以及(35)式的C1值至(32)式，再解一组非线性方程，可得到a2, as,…an各系数。 对于(32)式中的多个解，以及aa,a。…为零的情况，所有非线性系数仍然可以求得。 2.Hammerstein模型的辨识 图2为H模型 n(t)+d(t) A1(t) y(t) z(t) NL 图2H-模型 输入信号为A。L(t)+B,相应输出为y(t),假设 y(t)=yi(t)+y2(t) (36) 其中： y)=(a1+2a:B+aA.2+…0Ag(a)L(t-)d2 (37) y2(t)=-(a2A。2+a2B2+3aA。2B+a3B3+…) (38) 从以上两式可以推导出 y:(t)=-y:(t+T) y:(t)=y2(t+T) 以上从(6)至（5)式的结论仍旧成立，这时，A。L(t)与Z(t)在0至(r+2)T区 间具有权重平均的互相关函数可写成为： 中Lz(t)=A。2Sg(t)(a1+2a2B+a3A,2+3agB2+…) (39) 以上公式从0至T区间积分，并注意到(22)式，得到： n-1 (d-M-CB (30) 其中 M=∑a2n-1A。2i-2 (41) C1=2az+4a4A02+6a8Ag+… C2=3ag+10a6A。2+21a7A。‘+… Cg=4a4+20a。A。2+… C4=5a6+35a7A02+… (42) C。=6a。 Cn-:=nan 146

当 ，， 相 应有 小 二 小 ， ‘ ， ， … 一 式可写 成为矩 阵形式 ， 有 伏丛丝卫 一 卫 卫 解以上 矩 阵方程 ， 我们得 到 少 卫一 如匕竺 一 式的 值 至 式 ， 再解一组非线性方程 ， 代入 中间函数 ， ， 以 及 。 ， … … 。 各系数 。 可得 到 ， 对于 式 中的 多个解 ， 以及 ， 。 … 为零的 情 况 ， 所有非线性 系数仍然可 以求得 。 徽型 的拼识 图 为 模 型 图 一 模型 输入 信号为 。 ， 相 应输 出为 ， 假设 其 中 ， 。 … 广 。 入 一 入 入 “ 一 。 “ 。 … 从 以 上 两 式可 以推 导 出 一 以上从 至 式 的结论仍 旧 成立 ， 这 时 ， 。 与 在 至 区 间具有权 重 平均的 互相 关函数可写 成为 小 。 “ 。 … 以上公 式从 。 至 区 间积分 ， 并注意到 式 ， 得 到 一 命 ‘ 二 ‘ · ’ ’ 一 “ 乏 ， ’ 其 中 乙 ， 。 “ ’ 一 一 一 。 “ 毛 … 。 。 ‘ … 一 … ‘ 。 么 … 。

当B=B。=0,相应有中Lz(t)=中Lz.(r),求得M为： M--Ai(d (43) 当B=B1,相应有中Lz(t)=中Lz1(t),i=1,2,…n-1(42)式可写成矩阵形式如下 -Aiss(dt-M-CB 解以上矩阵方程，我们得到 C=B-Ais()d-M)) (44) 代入所得的C,值到(42)式，解一组线性方程组(42)，系数a2,a,…an便可求得。 代入(43)式中的M值，a2,a3…an到(41)式，得到 nI a1=M-∑a21-1A,2-2 (45) 1=2 在(39)式中，让B=B。=0,代入中Lz(T)到(39)式，g(t)的栅点值由下式得到 g(t)=中Lz(t)/a1A。2S (46) 方法二 1,Hammestein模型辨识 在图二中假设漂移干扰为零，即d(t)=0,输入信号AL(t),Z(t)与AX(t)的互相 关函数为中xz(t),AX(t)为与L(t)同步的m序列。 .(t)=JAz)×t-t)dt (47) 以及 中xz(t)=中xy(t)+中xn(t) (48) 由于 27'(L(t)a1-x(t-t)dt=0i=1,2… (49) 中xn(T)=0 (50) 假设 og(s)ds=1 (51) 得出 中a(t)=N (52) 式中N是m序列的长度。 当A=A1,相应有中xz1(t),i=1,2,…n,(52)式写成下列矩阵形式： 中xz(t)=A·a (53) 或 a=A-1中xz(t) (54) 全部偶次系数a2i,i=1,2,…n可以从(54)式求得。 147

当 。 ， 相应有小 小 乙 。 ， 求得 为 “ 一 万万丐 甲 · 。 、 ’ 。 当 ，， 相 应有小 小 。 ， ， ， … 一 式 可 写 成矩 阵形 式如下 式 了打‘ · 担 一 “ 丝 解以 上矩 阵方程 ， 我们得到 代入所得 的 值 到 式 ， 代入 式 中的 值 ， 一 一 云烤只恤 一 “ ’ 解一组 线性方程 组 ， 。 … 到 式 ， ， 系数 ， 。 ， … 。 便 可求得 。 得 到 一 乏 ，一 。 ’ 一 “ 在 式 中 ， 让 。 ， 代入 小 。 到 式 ， 的 栅点值 由下 式得 到 小 。 ， 方 法 二 橄纽拼识 在图二 中假设 漂移干扰为零 ， 即 二 ， 输入 信号 ， 关函数为中 ， 为与 同 步的 序 列 。 小 乙 下 二仁 一 ‘ 了 丫 卜 “ 、 一 ‘ “ 一 、 一 ‘ ” 、 一 ‘ 一 与 的 互 相 以 及 由于 假设 得 出 小 小 ， 二 小 制 ’ 〔 ‘ ，，〕 “ ’ 一 ” “ 一 ， ‘ 二 ” ‘ “ ‘ ， ， ’ “ 小 、 歹了 ‘ ， “ ‘ 小 票 一 安 。 ，一 、 ‘ 式中 是 序列 的 长度 。 当 ， 相 应 有 小 ， ， ， · ” 一 ， 式写 成下列 矩 阵形 式 小 下 。 或 全都偶次系数 ，， ， ， … 二 一 ’ 小 可 以从 式求得

由AL(t)与Z(t)的互相关函数 )=5gt（∑ (55) 可以求出奇次系数以及g(τ)。 2,Wiener模型辨识 在图1中假设d(t)=0,输入信号为AL(t),AL(t)与Z(t)的互相关函数写成 n L:(t)=∑a2-iA2p2a1-1() (56) 11 其中 a-(e)=U0gsL(t-s)da]21-"L(t-t)dt (57) 从上式我们可以求得奇次系数以及g(τ)。 下面求偶次系数。输入信号仍是AL(t),互相关信号是Z(t)和与L(t)同步的AX(t), 中xz(t)由以下给出 .)=∫AzX-)dt (58) 因为X(t)与L(t)不相关，有 L(-)dt-x(t-)dt =0 (59) 结果 p(t)=∑a1A21+021(t) (60) 其中 中e(r)=27∫[0gsL(t-sds]Xt-r)dt (61) 当A=A1,相应有中xz1(t),i=1,2,…n,(60)式可写成矩阵形式 中xz(t)=Aa中xu(T) (62) 或 a中x(r)=A中xz(r) (63) 代入g(t)至(61)，以及在(61)式求出的中xu21(t)代入(63)，所有偶次系数a2,a4, “a2n可以求出。比较(52)，(60)可知，由中xz（τ)是否为常量可以判定模型结构是H型 还是W型。 方法检验 为了验证以上提出的方法，辨识模拟实验是在FELTXC-256数字计算机上作的。在模 148

由 与 的互相关函数 ‘ 二 ‘ · ， “ · ， 宣一 ， 可 以求出奇次系数 以及 。 橄皿拼识 在图 中假设 “ ， 输入 信号为 与 的互 相 关函数写 成 中 · ‘ 乏 卜 “ ’ 小 。 ‘ “ ，一 ‘ ， ‘ ， 其 中 小 一 ， 二 斋‘ “ ’ ‘广 ，一 、 ‘ “ ’ 一 ” 卜 ‘ ‘ 夕 从上 式我们 可 以求得奇次系数 以 及 。 下面求偶次系数 。 输入 信号仍是 ， 互相 关信号是 和 与 同步的 ， 小 丫 由以下给出 。 · 责歹 ’ ‘ ， “ 一 ， ‘ 因为 与 不相 关 ， 有 。 。 卜 · 二 责 ’ 〔了 “ 一 ， ‘ ‘ “ ’ 一 ‘ ’ ‘ ，一 ， ， 结果 小 · ‘ 艺 ， “ ’ “ 小 · 。 “ ’ ‘ 其 中 当 ， 产 么 厂 产 、 中 “ ” 吸丫 ’ “ 介 一 。 日 “ ’ “ ‘ 一 “ ’ “ “ 入 “ 一 下 ’ “ ‘ 相应 有小 、 ， ， ， … ， 式可写 成矩 阵形式 小 小 或 代入 至 小 。 一 ‘ 小 ， 以 及在 式求 出的 小 代入 ， 所 有偶次系数 ， ， … 。 可 以求出 。 比较 ， 可知 ， 由小、 是否为常量可 以 判定模型结 构是 型 还是 型 。 方 法 检 验 为了验证 以 上提 出的 方法 ， 辨识模 拟实验是在 一 数 字计 算机 上作的 。 在模

拟试验中，非线性部件模型是 V=aip+azp2+a3p3+ap4+aop5 线性动态系统模型是 c1s2+c2s+C3 g(s)=5+b1s+b25+b3 为了简化，假设n(t)=0 模拟试验的辨识结果分别描述如下： 1.试验信号参数 电平，A:1.06,2.22 A。0.29 偏置B:-2,-1,1,1.5 钟周期△t:0.01S 采样区间SI:△t/4 m序列长度N:63 权重系数K1:1,3,3,1 漂移干扰 d(t)=1-2t+3t2 数字模拟辨识结果列举如表1与图3 表1 方法1中W模型辨识结果 模型 非线性部分 参数 a2 a3 理论值 .5 .4 .3 .2 .0 估算值 .503 .405 .306 .205 .0 模型 线性动态系统 参数 Ci c2 C3 bi b3 理论值 .0579 -680.0-681540.49712,9 6815 估算值 ,0468 -666.1-652739.75 702,5 6605 2. 方法1H模型辨识 试验信号参数 电平A:0.1 偏置B:-2,-1,1,1.5 钟周期△t:0.01S m序列长度N:63 采样区间SI:△t/4 数字模拟辨识结果列举如表2与图4。 149

拟试验中 ， 非线性部件模 型是 一 ‘ 。 线性动态系统模型是 名 为了简化 ， 假设 模 拟试 验 的辨 识 结果 分别描述如下 试脸信号今数 电 平 ， 。 偏 置 一 ， 一 ， ， 钟周 期 △ 采 样区 间 △ 序列 长度 权重 系数 ， ， ， ， 漂 移干扰 一 “ 数字模 拟辨识 结果列举如表 与图 表 方法 中 模 型辨 识 结果 参 数 理论值 估 算值 理论值 估 算值 方 法 徽型辨识 试 验信号参数 一 一 一 一 电 平 。 偏 置 一 ， 一 ， ， 钟周 期△ 序列长度 采样区间 △ 数字模拟辨识 结果列举如 表 与图

30V 4, XX估算值 20 o (t) 2 0.4 0.5 0.1 0.2 0.3 t(s) -5 —理论业 节名估募德 -10 图3由数字仿其的V-模型 -15 辨识结果 表2 方法1H模型辨识结果 模 型 非线性部分 线性动态系统 参 数 az a C3 b2 b3 理论值 .5 .4 .3 .2 .1.0579-680.8-681540.49712.96815 估算值 503.403.302 .201.101.0468-666.1-652839.75702.56605 8. 方法2W模型辨识 试验信号参数：电平A:1)1,2,3 2)1.5,2.5 钟周期△t:0.01S m序列长度N:63 采样区间SI:△t/4 150

表 方法 模 型辨识 结果 理 论 值 估 算 值 一 一 一 一 方 法 徽型辨识 试验信号参数 电平 ， ， ， 钟周 期 △ 序 列 长度 采 样区 间 △

50 不 40 30 20 10 0 10 ) 2 0.4 0.5 0 0.1 02 0.3 t(s) -51 避论雀 X×估算雀 10 -15 图4 由数字仿真的1-模型 辨识结果 数字模拟辨识结果列举于表3。 表3 方法2W-模型辨识结果 模 型 非线性部分 线性动态系统 参 数 al az aa a4 aa C2 C3 bi b2 b3 理论值 .5 .3 .2 .1 .0579680.8681540.49712.96815 估算值 .504.406.317.208 .110-.0468666.9653539.75702.56605 4.方法2H-模型辨识 试验信号参数 电平A:1)1,2,3 2)1.5,2.5 钟周期△t:0.01S m序列长度N:63 采样区间SI:△t/4 数字模型辨识结果列举于表4。 151

数字模拟辨识 结果 列 举于表 。 表 方法 一模 型辨识 结果 型 “卜 线 性 部 分 一 线 性 动 态 系 统 数 · 且 · 么 · · 二 。 理 论 值 · · ‘ · 一 。 一 摘竟 一 面恤 一而 一蔽 一压丽爪叫一成 方法 一徽型辨识 试验信号参数 电平 一 ， ， 钟周 期 △ 序列 长度 采样区 间 △ 数字棋型辨识 结果列 举 于表 。