系统辨识与滤波

第二讲辨识三要素一、数据本节介绍辨识数据的特点及获得适宜辨识的数据的方法。随机过程x(t):在每一个时间点t=t0上，x(t)都是一个随机变量，其概率密度函数p(x,t)随时间 变化。平稳随机过程：在所有时间点上，概率分布都相同的随机过程，其概率密度函数p(x)不随时 间变化。各态遍历平稳随机过程：从整个时间轴上看，每个随机事件都会发生的平稳随机过程。其谱 密度函数与概率密度函数类似。时间平均等于集合平均。数字特征特征随机过程平稳随机过程各态遍历平稳随机过程均值 (数学 期望值)H (t) = j 8 xp (x, t) dxX-8口 E(x(t)口 x )三 口 xp, (t)三口 = x=limj x (t) dt-T=lim L 义 x (k)N T8k =1均方 值W 2 (t) = j 8 x2 p (x, t) dx x-8 E x 2( t)W 2(t)三W 2xxW 2( t)三W 2 = x 2=limj x2 (t) dtL- T=limU x2 (k)N T8k =1方差。2(t) = E x(t)-日x(t)2=j8 x -日(t )2 p (x, t) dx-8* Ex(t) -|!x(t)2 Var x(t)Q 2 (t)三。2Q 2 =limj x(t) - p 2 dtL- T1 v-=lim 乙x(k) x2N Tsk =1自相 关函 数R (t , t ) Ex(t )x(t )X 1 212=j J x x p (x , x ; t , t ) dx dx 12 2 12 121 2一s 一s2Rx (t1,12)=R. (t2,.)=R (0, t - t )=R .(，)R (t ) = R ( -t )=limJ x(t)x(t + t )dtT T- T1羽=lim 乙 x (k) x (k + l)N Tsk=11N-1=lim乙 x (k) x (k + l)N ts N lk =1协方 差函 数C (t , t ) E x (t ) f (t )* x 121x 1x(12) f x(12) Cov x(t ), x(t )=J Jx(t ) -日(t ) x(t )-1x 12-s-s日(t ) p (x , x ; t , t ) dx dx x 22121212C x (t1,12)=C (12,t=C (0,12 - 11)=R 3)-日2=Cx (T )Cx (T ) = Cx (-t )limJ x (t) - 口 *Tt- Tx(t +t ) - 口 dt1 V _ 一 .lim 乙x (k) - x *=N ts Nk =1x (k + l) - x 1N l-lim乙x (k) - x *=Nts N lk =1x (k + l) - x 一些 关系w 2( t) = Rx (t, t)b 2(t) =w 2(t) - |! 2(t)xxx=R (t, t)-日 2( t)C (t , t ) = R (t , t ) x 12x 12日(t川(t )x 1x 2w 2 = R (0)b 2 = W 2 日 2=R (0)-日 2Cx (T)=R 印)口2w 2 = R (0)1乎、=lim 乙 x2 (k)NTsk=1b 2 = R (0) - x2C (t ) = R (t ) - /互相关函 数R (t ,t ) = Ex(t )y(t ) xy 1212=J J xy p (x, y; t , t ) dxdy 212-s-sR (t12)=R (t = t t )R , 3)=limJ x(t)y(t +t )dtTt- T1 v=宫 Xx (k”(小l1-1=limX x (k) y (k + l)N T8k =1互协 方差 函数C (t ,t ) = Ex(t )- xy 121日(t ) y (t )-日(t )x 12y 2带带=J J x -日(11) y -8-8日(t ) p (x, y; t , t ) dxdy y 2212=R (t , t )一日(t 川(t )x!1 2x 1y 2C (t12)=C (t = 12 - t =R (t ) -口口气3)=limJ x (t) - 口 *Tt- Ty (t + t ) - 口 dt一. 1 V 一 一 .lim X x (k) - x *=n * Nk =1y (k + l) - y 1n- ilimX x (k) - x *=NT8 N lk =1y (k + l) - y 其它R (t ) = R (-t ) Cxy 3) = Cyx (-T)说明：离散计算时假设采样时间间隔为七，则时延T= l * T0。相关函数的性质1 V 2 = R (0) N02 R(T ) = R (-T )3 R (0) I R(T ) |4若中)是周期为T的信号，则其自相关函数也是周期为T的信号。即:x(t) = x(t + T) n R (t ) = R (t + T)5 若 x(t)=y(t)+z(t)，且y(t)与 z(t)互不相关(R (t )三 0 )，则yzR (t ) = R (t ) + R (t )6若x(t)=y(t)+z，其中y = 0，z是一个常数，则R (t ) = R (t ) + z xy7若如)均值为零，且不含有周期性成分，则当r很大时，迫)与x(t+T )必然是互相独立的(不 相关)，因此，R (t ) = 0 , T充分大。8若x(t)均值为零，则C (T ) = R(T )。这是因为在通常情况下，C(T )等于R(T )向下平移日2，因此，当日2=0时，两者相等。9 对于线性系统 y(k)=G(z)u(k)，有R , (t ) = G (z)R (t )10虽然x(t)是个随机过程，但R (t )却不是随机过程，而是一个确定性的时间函数。js II X ( j)II2 d 2兀-s帕萨瓦尔Parseval)定理与功率谱Parseval定理：确定性信号x(t)的总能量为：Js x2 (t)dt =确定性信号x(t)的平均功率：s lim II X (j)l|2 -s 2 T Tlim j T x2 (t) dt = j2 T - t2兀,确定性信号x(的平均谱密度：1S ()=lik XI j(2 ) | |x 2 T t随机性信号x(t)的平均谱密度：1S ()=libE XI I j(2 ) I I xTts 2TT维纳一肯塔金(Wiener-Khintchine关系式:随机过程x(t)的谱密度Sx(o )与自相关函数Rx (t )构成一组傅立叶变换对:S (o ) = js R (t )e - jot dtx-s xR (t ) = 2-js S (o ) e j ot d o定义互谱密度为互相关函数的傅立叶变换:S (jo ) = js R (t )e - jot dtR (t ) = 土卜 S (o ) e jot d o应用维纳一肯塔金关系式，可以证明，对于频率响应为G(jo )的线性系统，在随机输入下 的输出谱密度和互谱密度分别为：s ()=| G (j)|2 s ()七(加)=G (加)S,()输出谱密度关系告诉我们：要充分激励系统，就要使输入信号的频谱“宽”于系统频谱。白噪声白噪声的概念：白噪声过程(一系列不相关的随机变量组成的理想化随机过程)相关函数：R (T ) = Q 25 (T )谱密度：S 3 ) =b 2(一8 W +8 )近似白噪声过程谱密度：s ()=b ，一 o(为给定的远大于过程的截止频率)相关函数：R (T )讨论白噪声时，还要涉及到白噪声的概率分布，服从正态分布的白噪声称为高斯白噪声。n维白噪声：一个n维随机过程W (t)满足:E(W (t) = 0、Cov W (t), W (t +t ) = EW (t)W (t + t ) = Q5 (t )其中Q为正定常数矩阵，则称W (t)为n维白噪声过程。白噪声特点：如果一个零均值、平稳随机过程的谱密度为常数，我们称之为白噪声(由白色光联想而 得)。一个平稳随机序列w化)具有恒定功率谱密度，若这个过程是离散的，称白噪声序列。白噪声有以下特点：1 E (w )= w = 0。2 S (少)=b 2，频谱宽度无限。3 . T = 0 一3 R(T) = Q 25(T)，其中，s(T )为脉冲函数(Dirac函数)，即8 (t ) = (且j 3 8 (t ) d t = 1。4无记忆性，即t时刻的数值与t时刻以前的过去值无关，也不影响t时刻以后的将来 值。从另一意义上说，即不同时刻的随机信号互不相关。白噪声的用途：1作为系统输入时，有R (T ) = g(T )，=0,1,2,.，即为系统的单位脉冲响应。2作为被辨识系统输入时，可以激发系统的所有模态，可对系统充分激励。3作为被辨识系统输入时，可防止数据病态，保证辨识精度。4在辨识过程中，以输出估计误差是否具有白色性来判断辨识方法的优劣，也可用来判 断模型的结构和参数是否合适。5产生有色噪声。白噪声的产生方法：1 (0，1)均匀分布白噪声：化=取小数179* & 1，i=1,2,3，初值可取为：曰 0 = 11 2 35。2正态分布白噪声n(口,q 2): m (k) = + b (i2 & 6)，其中e为服从(0，i=11)均匀分布的白噪声(教材51页)。例题：白噪声生成利用乘同余法生成U0,1均匀分布的随机数x.= Ax , (mod M ) x& = r U 0,1i M其中,M = 215 = 32768 ,(循环周期2S2)A = 179 三 3 (mod 8)x0 = 11,(奇数)利用U0,1均匀分布的随机数生成正态分布的白噪声v(k) = b丈 & 6 N (0,b 2 )飞.=1 )其中，标准差b分别取0,0.1,0.5。编程语句for (k = 1; k 252 , k + + ) ksai = 0;for (i = 1; i 12; i + + )xi = A * xi;xi = MOD( xi, M );ksai = ksai + FLOAT (xi/M );v(k) = Sigma * (ksai 一 6.0);有色噪声有色噪声是指每一时刻的噪声和另一时刻的噪声相关，因而其谱密度也不再是常数。在 工业生产实际中，白噪声在物理上是不存在的，常见的往往是有色噪声。表示定理：设平稳噪声序列e(k)的谱密度S ()是o的实函数，或是cos o的有理函数，那么必 定存在一个渐近稳定的线性环节，使得如果环节的输入是白噪声序列，则环节的输出是谱密 度为Se(o)的平稳噪声序列e(k)。成形滤波器：表示定理中所涉及到的线性环节称为成型滤波器。定理表明：有色噪声序列可以看成 由白噪声序列驱动的线性环节的输出(如下图所示)。白噪声线性环节，(成形滤波器)H(z-1)有色噪声w(k)e(k)可以证明：如果(e(k)的谱密度S (o )是cos o的有理函数，那么一定存在一个成型滤 e波器，它的脉冲传递函数为：D (z T)1 + c1 z t + + c ： z-气C (z -1)1 + d z -1 + + d z - nd且C (z -1), D (z -1)的根都在z平面的单位圆内。M序列(二位式最大长度伪随机序列PRBS-Pseudo Random BinarySignal)M序列序列的概念：1辨识输入信号要求具有白噪声的统计特性。2 M序列是PRBS的一种形式，它的自相关函数接近脉冲函数，M序列具有近似的白噪 声性质，有持续激励，可保证有好的辨识精度，即1 . T Const , t = 0R (t ) = j M (t) M (t + t ) dt I。3 M序列净扰动小，幅度、周期、时拍易控制，实现简单，不影响系统正常工作，可 进行在线辨识。例：4 阶 M 序列 1111010110010001 1 1 1。1 1 11 -11-111 -1-11 -1 -1-11 1 1 10 0 0 0 0M序列的产生:M序列的定义:x = a xa xa x , i = P + 1, P + 2, ，注意求和为模2i 1 i-12 i - 2P i - P求和（异或），P表示级数，且a? = 1，适当选取a 1?a2,-工八（反馈通道的选择）就可以使序列以（2P - 1） bit的最长周期循环，这种最长周期二元序列就称M序列。M序列的生成结构图（P55图2.29）M序列生成实例:用M序列作为辨识的输入信号，M序列的循环周期取七=21 = 15 ,时钟节拍 t = 1秒，幅度a = 1,逻辑0”为a,逻辑1”为-a，特征多项式自选，如尸（s） = s4s31，见P59表 2.11 中（P,n,k）=（4,2,2）或（,n2,%，）=（4,3,0）。M 序列的生成结构图为（P55 图 2.30）：在CP移位脉冲的作用下，P55图2.30中各级寄存器状态的变化如表所示（P55表2.10）。CP123456789101112131415161718 C1111100010011010111 C011110001001101011 2C3101111000100110101 C 4010111100010011010 用M序列作为辨识的输入信号M序列的循环周期取N=26 - 1 = 63 ,时钟节拍PAt = 1秒，幅度a = 1,逻辑0”为a,逻辑1”为-a，特征多项式自选，如尸(s) = s6s51 , 见 P59 表 2.11 中(P,n,k)=(6,2,2)或(,%,%，)=(6,5,0)。生成M序列的结构图各级寄存器状态的变化可自行推导。编程语句for (k = 1; k 0; i) M (i) = M (i - 1);ifM ( 0 ) = 0 thenu ( k ) = a;ifM ( 0 ) = 1thenu ( k ) = 一 a;注：若用MATLAB，可参看Idinput函数(用Help命令查看有关函数说明)特征多项式解决如何选取反馈通道的问题，以保证生成M序列。定义多项式：G (s) = xs，(无限阶)i重要结果：G ( s)=-F ( s )F (s) = 1 a sj (有限阶)j=1称F (s)为M序列的特征多项式，称F (s) = 0为M序列的特征方程。注意1：此时选取M序列初始状态为：x = 1, x = 0, , x = 0。注意2：生成M序列的结构图完全由特征多项式F (s)确定。注意3：不是任何多项式都可以作为生成M序列的特征多项式，它必须满足一下条件。必要条件：特征多项式F(s)是不可约多项式，即：不可分解的多项式。充分必要条件：特征多项式F ( s )是本原多项式，即F ( s )是多项式snp 1的一个因子，其中Np = (2 p - 1)，P是F (s)的阶次。满足以上两个条件的部分特征多项式见表2.11，注意表的使用，P确定时，构成M序 列的方案不唯一，其多项式F (s)也不唯一。M序列的性质周期性：M序列的循环周期Np = (2 p - 1)位，P (级)个移位寄存器所能表示的最多 状态数。出现0的次数为(Np-1)/2次，出现1的次数为(N疽1)/2次，一个周期中1的 个数比0的个数多1。例如：4阶M序列中有8个1，7个0。M序列的游程：M序列中某种状态连续出现的段称为游程。一个P级M序列的游 程总数为2p 1，其中0游程与1游程各占一半。长度为i位(1 i P - 2)的游程占1/2，即有2P-I个，但长度为(P - 1)位的游程只有一个，为0游程，长度为P位游程也 只有一个，为1游程。M序列的可加性：所有M序列都具有移位可加性(模2和)，模2和后的结果仍为M 序列，并与原M序列移位等价。M序列的自相关函数幅度的选取：作变换：M （，）= 0（1 - 2x.）,幅度变为G和-1。计算式：R+其中，是位移脉冲周期（时钟节拍）。M N At 0 p结果：RM/a 2 1 -离散形式为：R (t) = (k)M (k + x) M N p k =0-Ar t ArO据自相关函数的性质知,Ar t (N l)ArR （T）也是N ，的偶函数，如图所示（P64图2.34）0MP当心充分大时，M序列自相关函数近似脉冲函数，故M序列是较理想的辨识的输入信 号。然而，用M序列作为辨识的输入信号时，需要根据辨识对象的先验知识，适当选择M序列的循环周期N 3、移位脉冲周期（时钟节拍）函和M序列的幅值s才能获得比较P理想的辨识结果。M序列的谱密度：R (T ) = 7?(1； ) + A (2)(T )MMM(j27?( 2 饥)=MNpS (co) = S(i)(oo) +S(2)(co)MM005(i)(co) = 2k E c 8 (co - ZrcoMkk = -s2兀 co =0 N pS(2)(co)=M8 (co)2兀0 2 (N + 1) 1 一 sin (/ 2200 8 (cd kco )0k o以。0-rN 2P1CD Af_ 2_S (cd )=M+271(2 2N 28 (co)M序列的频谱：设，为一步M序列信号的持续时间，心勇是一个周期的持续时间。M序列的频谱为：必t2 兀 2sinS ()=-d ()+ (N + 1) 0；2 2 Z 5 (kro )，其中: pk = s2k丰0上式告诉我们:1M序列的频谱不是光滑的曲线，而是线条谱。2 M序列的直流分量（S（ro）k=红竺）与N；成反比，因此，加大Np，可减少M序P列中的直流分量。3 M序列的频带为B = 1 Hz，因此，减少At，可增加带宽。3 A t4谱线密度与NpAt成正比，M序列应用：将M序列作为扰动信号有以下好处：1幅值可取a,-a，容易选择，且当Np充分大时，均值约等于02在一个周期内，自相关函数匕近似为函数，因此，以M序列为输入的线性系统，其互相关函数序列等于脉冲响应序列(*大于过渡过程)。3对于多输入单输出系统，可将同一 M序列的不同移位序列(例：M(k)、M(k+j)、 M(k+2j)等)作为各输入信号的扰动信号，当输出对各输入的过渡过程小于j时，可认为输 入信号之间是互相正交的。逆重复M序列：逆重复M序列的提出：谱分析表明：M序列含有的直流分量(Sm S) k=空虹 与Np2P成反比)，这会造成对辨识对象的“净扰动”，这通常是不希望的。而逆重复M序列可克服 这一缺点。逆重复M序列的产生：将M序列M(k)与0或1的方波序列S(k)按模2和，可得逆重复 M 序列IM(k)，即：(IM (k)=( M (k)(S(k)。如：M(k)S(k)IM (k)11110001001101011110001001101010 1010101010101010101010101010正M序列逆M序列(即正M序列对应位按位求反)0 10 110111001111101001000110000逆重复M序列的特点：1周期=2Np，一个周期中1的个数与0的个数相等，EIM(k)=0，即IM(k)是零均值序列。 具有白噪声的优点，而其幅值可控。2逆M序列的前半周期与后半周期为逆重复，即：IM(k)= -IM(k + Np)。3逆M序列IM(k)与原M序列M(k)不相关，即它们的互相关函数等于零。辨识试验设计设计原则：在安全的前提下，尽可能地激励系统；保持输入输出关系；适当解耦；明确目的、 要求和模型用途。了解辨识对象划清要辨识系统的边界，选好输入/输出，从边界外连入的其它信号尽量保持稳定，并 作为被辨识系统的噪声。整体/局部：确定哪些输入需要叠加扰动信号，哪些输入要保持稳定。输入/输出/噪声：确定是多输入还是单输入(耦合关系)、确定过渡过程是否有明显差异(时间 常数)、了解噪信比的大小(滤波)及噪声类型(白色、有色)。值域范围：确定信号采集时是否需要零迁(加、减直流分量)、放大。安全工况：可叠加扰动信号的类型与幅值选择工况：生产负荷、试验时间、系统区域隔离、地理区域隔离、安全措施。扰动信号设计要点：扰动信号频带应宽于过程的工作频带；持续时间为3-5个扰动信号周期；幅值由安 全工况确定，对可中断生产的系统，且试验不引起原材料浪费的，可实施单独试验。扰动信号类型：M序列、逆M序列、白噪声序列等。扰动信号幅值：由安全工况确定。M序列设计一级数P、步宽At、正交化设系统过渡过程时间为乙最高工作频率为fmax，(通常，fmax=(315)/T)则以1.5*T3* fmax为原则，选择P和At。实际应用过程：1确认系统过渡过程T；2选择最高工作频率fmax=(315)/T)；(设为fmax =k/T)；3 根据 1/At 3* fmax，得出 At =T/(3k)；4 根据 1.5*T1.5*(输出对各输入的过渡过程时间)。白噪声序列：1/A/ 3* fmax数据采集设计采样间隔：满足香农定理，设采样间隔为T0, 1/ T02*fmax。对于有扰动信号的系统，T0 At , 且At应是T0的整数倍。实验时间：数据长度为3-5个扰动信号周期，预实验1-2个扰动信号周期。物理信号处理：信号零点选择、滤波器选择、放大器选择。*现场准备：(扰动信号发生器：测试、加限幅器、选择叠加方式、接线。信号采集装置：测试、时间间隔标定、接线试验设计报告：上述内容及人员组织，报领导批准)数据预处理传感器非线性校正：根据传感器生产厂提供的静态非线性曲线，编制程序对相应的采样数据 进行校正；滤波：对于输入输出加相同的滤波器，输入输出关系不变；重抽样：对于过渡过程较长的子系统，可用重抽样的方法，重新选择采样间隔，以降低模型 阶次。标准化：零均值、单位方差目的：降低舍入误差，防止病态方程，提高计算精度。方法：y(k )+ay(k-1)=bu(k -1)+cu(k-2)+d(1)y(k)/ N+ Eay(k -1)/N= Ebw(k-1)/N+ Ecw(k-2)/N+dY =y(k )/N= y(k-1)/NU =u(k-1)/N= E u(k-2)/N,d=(1+a) Y -(b+c)UY(k)= y(k)- YU(k)=u(k )-U(2)Y(k)+aY(k -1)=bU(k -1)+cU(k-2)a 2y= E Y(k)2/N , a 2u= E U(k)2/NY。(k)=Y(k)/ a y , Ua(k)=U(k)/ a uY (k )+a Y(k-1)=b Ua(k-1)* a J a +c U。(k-2) * a J a yB= b* a u/a y,C= c* a u/a yY (k)+a Y (k-1)=B Ua (k-1)+C Ua (k-2)b= B * a y/ a u, c= C*a y/ a u -反演公式二、辨识准则误差定义参数误差：J=E a 2 /n适合于算法研究;、一、一、，L_八_一八、一、一方程误差：J=E e(k)2/Ne(k)= y(k)+ a y(k-1)-b u(k-1)-c u(k-2), J 是参数的线性函数；只能作一步预测；适合模型为：AY=Bv UU2；、I 、，一 L一-.人一.I输出误差：J=E e(k)2/N, e(k)= y(k)- y(k), y(k)=-ay(k-1)+bu(k-1)+ cu(k-2), J 是 参数的非线性函数；当输入序列已知时，可作多步预测；适合模型为：Y = / U1 + A2- U 2白化误差 加权准则：J= e(k) TWe(k) /N。模型结构对误差的影响模型阶次越高，J值越小，导致模型过估计。采样间隔对模型结构及误差的影响采样间隔增大一倍，模型阶次大致减少一半；要进行不同采样步长下模型的 比较，误差公式应改为：J= e(k) TWe(k) *T2 /N , T0为采样步长。