基于Matlab的ARMA模型时间序列分析法仿真

基于Matlab的ARMA模型时间序列分析法仿真

ARMA模型时间序列分析法简称为时序分析法，是一种利用参数模型对有序随机振动响应数据进行处理，从而进行模态参数识别的方法。参数模型包括AR自回归模型、MA滑动平均模型和ARMA自回归滑动平均模型。1969年Akaike H首次利用自回归滑动平均ARMA模型进行了白噪声激励下的模态参数识别。 N个自由度的线性系统激励与响应之间的关系可用高阶微分方程来描述，在离散时间域内，该微分方程变成由一系列不同时刻的时间序列表示的差分方程，即ARMA时序模型方程：

k02Nakxtkbkftk (1)

k02N 式(1)表示响应数据序列xt与历史值xtk的关系，其中等式的左边称为自回归差分多项式，即AR模型，右边称为滑动平均差分多项式，即MA模型。2N为自回归模型和滑动均值模型的阶次，ak、bk分别表示待识别的自回归系数和滑动均值系数，ft表示白噪声激励。当k＝0时，设a0b01。 由于ARMA过程{xt}具有唯一的平稳解为

xthifti (2)

i0式中：hi为脉冲响应函数。 xt的相关函数为

RE[xtxt]i0k0hihkE[ftiftk] (3)

ft是白噪声，故

2kiE[ftiftk] (4)

0other式中：2为白噪声方差。

将此结果代人式(3)，即可得

R2i0hihi (5)

因为线性系统的脉冲响应函数hi，是脉冲信号，激励该系统时的输出响应，故由ARMA过程定义的表达式为

k02Nakhtkbktkbt (6)

k02N2N 利用式(5)和式(6)，可以得出：

k02NakRlkhiakilki0k02i0hibil (7)

对于一个ARMA过程，当是大于其阶次2N时，参数bk＝0。故当l>2N时，式(7)恒等于零，于是有

-1

- 好好学习，天天向上

RiakRlk0,l2N (8)

k02N或写成

k02NakRlkRi,l2N (9)

设相关函数的长度为L，并令M＝2N。对应不同的l值，由代人以上公式可得一组方程：

a1RMa2RM1aMR1RM1aR1M1a2RMaMR2RM2 (10) aRaRaR2L2MLMRL1L1将式(10)方程组写成矩阵形式，则有

RMRM1R1a1RM1RaRRRM2M1M22 (11) RRRL2LMaML1RL或缩写为

[R](LM)M{a}M1{R}(LM)1 (12)

式(12)为推广的Yule-walker方程。一般情况下，由于L比2N大得多，采用伪逆法可求得方程组的最小二乘解，即

{a}([R]T[R])1([R]T{R}) (13)

由此求得自回归系数ak(k1,2,,2N)。

滑动平均模型系数bk(k1,2,,2N)可通过以下非线性方程组来求解：

2b02b12bMc0b0b1bM1bMc1 (14)

b0bMcM其中

cki02Nj02NaiajCkij,k0,1,2,,2N (15)

式中：Ck为响应序列xt的自协方差函数。

滑动平均模型MA系数bk的估算方法很多，主要的有基于Newton-Raphson算法的迭代最优化方法和基于最小二乘原理的次最优化方法。 当求得自回归系数ak和滑动均值系数bk后，可以通过ARMA模型传递函数的表达式计算系统的模态参数，ARMA模型的传递函数为

H(z)k0k02N2Nbkzk (16)

akzk用高次代数方程求解方法计算分母多项式方程的根：

-2

- 好好学习，天天向上

k02Nakzk1a1za2z2a2Nz2N0 (17)

或表示成以下形式的方程：

z2Na1z2N1a2N1za2N0 (18)

求解得到的根为传递函数的极点，它们与系统的模态频率k，和阻尼比k的关系为

2zkexp(skt)exp[(kkjk1k)t] (19) 2zkexp(skt)exp[(kkjk1k)t]并且由式(19)可求得模态频率k，和阻尼比k，即

Rklnzksktkr|Rk|/t (20) 1k21(Im(R)/Re(R))kk为计算模态振型，需要先求出留数。设q点处激励p点响应的传递函数Hpq(s)的第是阶留数为Akpq，可用下式计算留数：

AkpqlimHpq(z)(zzk)zzkk0k02N2Nbkzk(zzk)akzkzzk (21)

振型向量可以通过对一系列响应测点求出的留数处理得到。对于一个有n个响应测点的结构，首先需要从n个对应同一阶模态的留数中找出绝对值最大的测点，假设该点是测点m，对应第k阶模态的归一化复振型向量可由下式求出：

{k}[Ak1qAk2qAknq]T/Akmq (22)

-3