随机子空间识别

第2节：模态参数识别的时域随机子空间方法主讲人：鲍跃全博士、副教授2013.6.6哈尔滨工业大学第6章模态参数识别课程名称《结构损伤识别与健康监测》主要内容3.4.2动力学状态空间模型连续时间状态空间模型离散时间状态空间模型随机状态空间模型3.4.3协方差驱动随机子空间识别方法Hankel矩阵的定义3.4.1随机子空间识别主要思想3.4.1随机子空间识别主要思想随机子空间识别（StochasticSubspaceIdentification-SSI）随机外部激励不可测，假设为白噪声，又称为环境激励。子空间识别模态参数：频率、振型、模态阻尼比时域状态空间方程为基本模型从输出数据的Hankel矩阵的投影的行“子空间”与列“子空间”中获取模态参数协方差驱动随机子空间识别方法数据驱动随机子空间识别方法系统矩阵模态参数Hankel矩阵计算输出数据的协方差矩阵奇异值分解将来行空间投影到过去行空间得到投影矩阵卡尔曼滤波+QR分解3.4.2动力学状态空间模型（1）连续时间状态空间模型n个自由度的系统，其动力特性可用二阶振动微分方程式表示tttttMqCqKqfBu通过定义0,,00tttqCMKxPQqMM0tttBPxQxu将二阶微分方程变成一阶微分方程:两边左乘P-1111110MPMMCMcctttxAxBu可以得到连续时间状态空间方程：式中：1110cIAPQMKMC1100cBBPMB0tttBPxQxu在实际振动测试中，一般不会测量全部自有度的输出，而只检测其中的一部分。假定只在其中lc个自由度处放置传感器，则结构的输出可以表示为：avsttttyCqCqCq11csavaCCCMKCCMC1caDCMBtqtxqtttttMqCqKqBucctttyCxDu观测方程cctttxAxBucctttyCxDu连续时间状态空间模型：110cIAMKMC10cBMB11csavaCCCMKCCMC1caDCMB状态空间矩阵输入矩阵输出矩阵直馈矩阵（2）离散时间状态空间模型cctttxAxBu的一般解：实际测试中，信号的量测都是在离散的时间点，因此需要将连续时间模型转换成离散时间模型。假定信号满足采样定理（Shannon采样定理），设初始时间t0000expexptccctttttttdxAxABu000expexptccctttttttdxAxABu11expexp1ktcccktkttktktdxAxABu设采样时间间隔为，t采样时间序列为0，0tkt1tkt,t2,t,ktTTTkkkktxxqq表示由采样开始时刻位移和速度向量组成的状态向量1kx表示k+1时刻的状态1kt令u在一个采样间隔内是常数假设1kkkxAxBu即expctAA10exptccccdBABAIAB10expexptkckccktdxAxABu类似kkkyCxDucCCcDDkuky表示采样时刻的输入与输出kkkyCxDu1kkkxAxBu离散时间状态空间模型（3）随机状态空间模型kkkkyCxDuv1kkkkxAxBuw实际工程中，数据测量噪声不可避免式中wk是处理过程与建模误差引起的噪声vk传感器误差引起的噪声离散随机状态空间模型kw1kxkxkvkyAC实测过程中，环境激励是不可测量的随机激励（输入），而且强度基本和噪声影响相似，无法将两者区分清楚。将输入项与噪声合并之后，就得到了随机子空间方法的基本模型kkkyCxv1kkkxAxw随机状态空间模型的重要性质假定噪声项为均值为零的平稳过程E0,E0kkwv0500100015002000-4-2024时间序列点幅值举例：实际的白噪声几个前提051015202530354045500246810频率幅度E0,E0TTkkkkxwxv同时wk,vk是零均值白噪声序列，与结构真实状态无关，可知：E,E0Tkkkxxx式中状态协方差矩阵Σ与时间k无关。系统为线性时不变系统，状态序列为平稳随机过程。ETikikRyy定义输出协方差矩阵Ri式中i为任意时刻的时延。0EETTkkkkkkTRyyCxvCxvCΣCR计算状态协方差矩阵Σ，可以得到11=ETkkΣxxEETTTkkkkAxxAwwETkkkkAxwAxwTAΣAQ1iTiRCAG计算协方差矩阵Ri（i=1,2,…）协方差驱动随机子空间识别方法建立的核心表达式。1ETkkGxy定义状态-输出协方差矩阵GEkkkkAxwCxvTAΣCS输出信号y1(t)和对应协方差序列3.4.2协方差驱动随机子空间识别方法0121123112021121123212212221jjiiiijiiiiijiiiijiiiijjyyyyyyyyyyyyYyyyyyyyyyyyy(1)Hankel矩阵定义逆对角线上元素相同yi是l×1维列向量，表示i时刻所有测点的响应0121ipfiipastfutureYYYYj列i块行i块行01211231120121021123121212212221jjiiiijiiiiijpiiiiijfiiiiiijpastfuturejyyyyyyyyyyyyYyyyyYYyyyyYYyyyy另一种分法(2)输出协方差矩阵+=ikikERyyiR是ccll维式中i为任意时刻的时延。-1+=01=limjTikikjkjRyy-11+1212-12-2==iiiiTifpiiiRRRRRRTYYRRR假定输出数据具有各态历经性，则实际测试中，只能得到有限的j个数据点，因此输出协方差可以估计为-1+=01ˆ=jTikikkjRyy-11+1212-12-2ˆˆˆˆˆˆ==ˆˆˆiiTiiifpiiiRRRRRRTYYRRR因此(3)块Toeplitz矩阵分解-1-2-112-22-3-1==iiiiTifpiiiCAGCAGCGCAGCAGCAGTYYCAGCAGCAG2cliniRO可观矩阵2cnliiRΓ可控反转矩阵-1-1==iiiiCCAAGAGGOΓCA对Toeplitz矩阵进行奇异值分解，1111211122==0TTTiTSVTUSVUUUSVSVU是正交矩阵，满足V==TTUUUUI==TTVVVVI奇异值按降序排列12=0SSS1=idiagS120n(4)模态参数识别可将奇异值分解结果分成两个部分，分别用矩阵iOiΓ表示1211=iSOUSΤ1211=iSΓSVΤSΤ是一个非奇异矩阵，可以看成是对原来模型的一种相似变换，即不管怎么取，得到的模型都是等价的1211=iOUS1211=iΓSV=SΤI取-11-1==iiiiiCCATAGAGGOΓCA=1:,:iclCO=:,-1+1:iccliliGΓ根据1iT的定义，同样可以得到2+1iT2-12+1==iiiiiCACATAGAGGOAΓCA1iT2+1iT和具有相同的结构，只是其中包含的协方差kR时延从2到i+1++2+1iiiAOTΓ+表示广义逆1212112+11=iiSUTVS离散状态矩阵A1=AΨΛΨ=nnidiagRΛ是一个对角阵，由离散复特征值组成。Ψ特征向量组成的矩阵离散状态矩阵与连续状态矩阵的关系，=expctAA-1=expccctAΨΛΨ连续状态矩阵cA1=cccAΨΛΨ-1=expccctΨΛΨA与cA具有相同的特征向量，两者特征值的关系为ln=exp=iiiitt系统的复特征值与固有频率，阻尼比之间的关系*2,=-1-iiiiiij22=RIiii=RiiiclnRΦCΨ频率阻尼比振型下弦杆腹杆上弦杆下弦杆腹杆上弦杆例子1002003004005006007008009001000-202x10-3时间(s)加速度(m/s2)0123400.20.40.60.81Frequency(Hz)SpectralDensitychannel5channel6channel7channel10channel11channel12channel17channel19channel21channel23channel26channel27channel29频率(Hz)功率谱0123400.20.40.60.81Frequency(Hz)SpectralDensitychannel5channel6channel7channel10channel11channel12channel17channel19channel21channel23channel26channel27channel290123400.20.40.60.81Frequency(Hz)SpectralDensitychannel5channel6channel7channel10channel11channel12channel17channel19channel21channel23channel26channel27channel29频率(Hz)功率谱510152025303540455000.10.20.30.40.5南北方向长度振型系数第1阶振型（东西向）频率：1.0202阻尼比：0.012001020304050500102030405050-0.4-0.200.20.4南北方向东西方向振型系数东西方向振型系数南北方向01020304050500102030405050-0.4-0.200.20.4南北方向东西方向振型系数01020304050500102030405050-0.4-0.200.20.4南北方向东西方向振型系数01020304050500102030405050-0.4-0.200.20.4南北方向东西方向振型系数01020304050500102030405050-0.4-0.200.20.4南北方向东西方向振型系数东西方向振型系数南北方向第3阶振型（竖向）频率：1.1272阻尼比：0.0184