简介自回归滑动平均ARMA模型时间序列分析法是结构动力学与工程领域中常用的模态参数识别手段特别适合从白噪声激励下的随机振动响应数据中提取自然频率、阻尼比与振型等动态特性。文档共5页系统梳理自回归AR模型、滑动平均MA模型及两者结合的ARMA模型数学表达并围绕推广的Yule-Walker方程、系数估计、传递函数极点求解等关键环节给出公式推导。内容还说明如何由极点换算模态频率与阻尼比并通过留数处理得到归一化复振型向量。文档虽短但步骤紧密适合力学、土木、机械等专业学生与振动测试工程师快速建立自回归滑动平均时序分析框架也可作为课程复习或项目实践的公式速查手册。包体仅含一个PDF文档压缩后约199KB轻量便携已有819人浏览学习是一份兼顾原理讲解与工程应用的简明资料。1. ARMA 模型为什么能让响应信号直接通向模态参数做模态测试的人都清楚这类场景结构太大太重力锤敲下去能量根本覆盖不了低频激振器能激励起来但安装和标定耗时太长环境激励下激励又是不可测的。频响函数法卡在最基本的前提上必须同时拿到输入和输出。反直觉的地方在于ARMA 模型可以完全绕开激励测量把响应序列本身当作分析对象用时间序列里的统计规律把系统的频率、阻尼和振型参数解出来。这就是标题里的“时序分析法”也是 ARMA 模型在工况模态分析里的核心价值。它适合桥梁、风力发电机、大型机械这类激励不可测的平稳随机响应场景。下面按“模型原理 → 公式推导 → 代码落地 → 工程排错”这条主线把整条路径讲透。2. ARMA 模型的表达方式与时序分析法的前提假设2.1 从结构振动方程到 ARMA(p, q) 差分方程多自由度结构的离散振动方程可以写成矩阵形式M·x C·x K·x f(t)其中 M、C、K 分别是质量、阻尼和刚度矩阵f(t) 是外部激励。对响应做等间隔采样后单测点的位移或加速度响应 y_t 在离散域满足一个线性差分方程y_t a1·y_{t-1} a2·y_{t-2} ... ap·y_{t-p} e_t b1·e_{t-1} ... bq·e_{t-q}等号左边由过去的响应值递归决定右边是白噪声序列 e_t 及其延迟项的线性组合。p 和 q 分别是自回归阶数和滑动平均阶数这个模型就是 ARMA(p, q)。结构动力学方程被离散化后自然呈现出这种递归关系而不是显式的传递函数形式这正是时序分析法区别于频响函数法的起点。2.2 AR 部分决定结构极点MA 部分处理激励成形ARMA 模型把系统响应分成了两层含义。AR 部分即等号左边的过去响应加权它完全由系统自身的 M、C、K 决定承载着模态频率和模态阻尼的核心信息。MA 部分则刻画激励的成形作用当外部激励不是理想白噪声而是带有一定频带特征的有色噪声时MA 项相当于一个成形滤波器对白噪声进行平滑后再驱动结构。理解这一点对后续求模态参数非常重要。模态识别的落点几乎全在 AR 系数上MA 系数虽然不直接影响频率和阻尼但参与了响应能量的分配决定各阶模态在实测信号里的激励权重。忽略 MA 部分退化成纯 AR 模型当然也能算但在激励是有色噪声时AR 系数会产生系统性偏差。2.3 平稳性、可逆性与 ARX 模型的边界ARMA 时序分析法成立的前提比频响函数法多两个硬性条件。第一是平稳性响应的均值和自相关函数只能与时间差有关结构出现明显退化、温度漂移或慢变载荷时数据必须预先处理。第二是可逆性MA 部分的特征根必须落在单位圆内否则白噪声序列无法由响应反解出来。下表把几种相近模型放在一起对比便于明确 ARMA 的使用边界。模型表达式要点激励假设识别目标AR(p)y_t Σa_i·y_{t-i} e_t激励为白噪声极点频率、阻尼MA(q)y_t e_t Σb_i·e_{t-i}描述有色激励成形激励谱ARMA(p,q)上述两者组合白噪声通过成形滤波器极点 留数ARXy_t Σa_i·y_{t-i} Σc_i·u_{t-i} e_t激励 u_t 可直接测量频响函数实际工程里使用 ARMA 而不是 ARX 的最常见场景是环境激励脉动风、地面微振、行车荷载都可以视为平稳随机激励测不到具体的 u_t只能把输入建模成白噪声驱动的成形过程。3. 模态参数识别的那条推导主线从 ARMA 系数到频率和阻尼3.1 由差分方程写出离散传递函数将 ARMA 差分方程两边做 z 变换忽略初始条件后得到离散传递函数H(z) B(z) / A(z)其中分子和分母分别是B(z) 1 b1·z^-1 b2·z^-2 ... bq·z^-q A(z) 1 a1·z^-1 a2·z^-2 ... ap·z^-p频响函数对应 z e^(jωTs) 时 H(z) 的取值而共振峰出现在分母 A(z) 趋近于零的位置。所以结构的模态信息被压缩到了分母多项式 A(z) 的根里面这就是为什么后面所有的参数识别都围绕 AR 系数展开。MA 部分只通过 B(z) 零点改变各频率的激励幅值。3.2 核心推导极点映射与频率阻尼换算A(z) 0 是一个 p 次代数方程两边乘以 z^p 后可以用标准的求根方法处理A(z) 1 a1·z^-1 a2·z^-2 ... ap·z^-p 0 等价形式 z^p a1·z^(p-1) a2·z^(p-2) ... a(p-1)·z ap 0 解出 p 个离散极点 z_r其中复极点成共轭对出现。离散极点 z_r 与连续域系统极点 s_r 之间存在指数映射关系。设采样周期为 Ts 1/fs则s_r (1/Ts) · ln(z_r) fs · [ln|z_r| j·(arg(z_r) 2πk)]这里只取 k 0 的主值分支。k 非零时对应频率折叠后的镜像极点一般直接丢弃。得到连续域极点 s_r σ_r j·ω_dr 后第 r 阶模态的固有频率和阻尼比由下面两式计算ω_nr |s_r| sqrt(σ_r² ω_dr²) ζ_r -σ_r / |s_r|注意 σ_r 是负数所以阻尼比为正数。工程上更常用的形式是用有阻尼频率 f_d 表述f_d ω_dr / (2π)且 f_d 与无阻尼固有频率 f_n 满足 f_d ≈ f_n·√(1 - ζ_r²)小阻尼时二者几乎一致。提示离散极点 z_r 的模长如果大于 1取对数后 σ_r 会变成正数对应阻尼比为负这是系统发散的不稳定极点在识别结果中必须剔除。3.3 复模态极点配对与多测点振型提取多自由度系统的每一阶模态对应一对共轭复极点例如 s_r σ jω_d 和 s_r* σ - jω_d 是同一阶模态。识别时只需要取虚部为正的上半平面极点并按频率远近将它们配成一对防止不同阶模态的极点交叉配对。振型的提取则依赖留数。对每个测点的响应序列分别拟合 ARMA 模型AR 部分理论上对所有测点一致因为系统本身的极点是全局属性但各测点在同一频率处的残差幅值不同这个幅值比就是振型分量。实际操作中逐测点拟合的 AR 系数会有微小差异常见做法是取所有测点 AR 系数的平均值或做全局 ARMA 拟合再用每测点的 MA 系数计算留数从而得到完整振型。4. 把公式落到代码两阶段最小二乘实现 ARMA 参数识别4.1 为什么不用直接法估计 ARMA 系数ARMA(p, q) 的参数估计是非线性问题因为等号右边的 e_t 不可观测。直接做最大似然估计需要迭代优化收敛性和初值选择都比较麻烦。工程中最稳妥的替代方案是两阶段最小二乘第一阶段用高阶 AR 模型对响应做预白化把残差当作对 e_t 的估计第二阶段把残差序列代入回归矩阵用普通最小二乘同时解出 AR 和 MA 系数。4.2 预白化与线性回归的完整实现import numpy as np def ar_prewhiten(y, L20): 一阶段用 AR(L) 预白化返回残差序列估计的 e_t n len(y) y_mat np.stack([y[i:n-Li] for i in range(1, L1)], axis1) target y[L:] coef, *_ np.linalg.lstsq(y_mat, target, rcondNone) resid np.zeros_like(y, dtypefloat) resid[L:] target - y_mat coef return coef, resid def arma_fit(y, p, q, pre_order20): 二阶段线性回归估计 ARMA(p,q) 系数返回 a、b、残差 _, e ar_prewhiten(y, pre_order) n len(y) rows, targets [], [] start max(pre_order, p, q) for t in range(start, n): row [-y[t-1-i] for i in range(p)] [e[t-1-i] for i in range(q)] rows.append(row) targets.append(y[t]) coef, *_ np.linalg.lstsq(np.asarray(rows), np.asarray(targets), rcondNone) a np.r_[1.0, coef[:p]] b np.r_[1.0, coef[p:]] return a, b, e回归矩阵中 y 列取负号是为了让估计出来的系数直接匹配 A(z) 1 a1·z^-1 ... 的形式。e 列不加负号MA 系数直接对应 B(z)。lstsq 返回的最小二乘解在数据量几千点时已经足够稳定数据更长时建议换用递推最小二乘减小内存占用。4.3 从 a 系数换算频率和阻尼的代码片段def arma_modes(a, fs, f_max50.0, zeta_max0.2): 由 AR 系数求极点并换算模态频率与阻尼比 z np.roots(a) s fs * np.log(z) modes [] for sk in s: if sk.imag 0: continue wn abs(sk) freq sk.imag / (2 * np.pi) zeta -sk.real / wn if 0 freq f_max and 0 zeta zeta_max: modes.append((freq, zeta, wn)) return modes核心是 np.roots(a)。a 的最高次项为 1roots 返回 z 平面极点np.log 完成从 z 平面到 s 平面的映射。实际测试中你会发现频率越高极点越靠近单位圆log 映射对极点位置误差越敏感因此高频段的阻尼比往往比低频段更分散。f_max 过滤掉超出关注频带的极点zeta_max 防止把奇异解当作真实模态。4.4 阶次 p、q 怎么定AIC 与奇异值联合判断ARMA 阶次没有先验值时最常用的做法是扫描 p 和 q对每组阶次计算信息准则AIC N · ln(σ_e²) 2·(p q) BIC N · ln(σ_e²) (p q)·ln(N)其中 σ_e² 是残差 e_t 的方差N 是样本点数。AIC 随着阶次增加先快速下降再缓慢爬升拐点对应的阶次就是推荐值。由于模态数是未知的工程习惯是让 p 至少覆盖可能模态数的两倍再配合叠加图或稳定图做最终确认。仅用 AIC 容易选偏低的阶次建议将 AIC 最小值附近的多个阶次全部计算模态参数观察哪些频率和阻尼随阶次变化保持稳定。5. ARMA 模态识别落地时的 4 个细节与稳定图检验5.1 采样频率与识别频带要有明确边界采样率不是越高越好。采样率过高会让关注频带只占整个频谱的一小部分ARMA 模型的拟合精力被大量分配到带外噪声上过低则高频模态折返混叠。一般让目标最高模态频率落在 0.1 到 0.25 倍 fs 之间。例如目标最高模态 20 Hz采样率设定在 80 到 100 Hz 更合适而不是直接使用采集仪默认的 1 kHz 档位。5.2 趋势项和低频漂移必须先处理环境激励数据里常常叠加大幅值低频漂移这直接破坏平稳性前提。进入 ARMA 拟合前用一次差分或高通滤波把直流和极低频分量滤掉是标准操作。需要注意的是差分会放大高频噪声数据信噪比不高时优先选择高通滤波器截止频率设在关注频带最低频率的 1/3 左右避免把真实低频模态一同滤掉。5.3 数据长度不足时优先缩短阶次ARMA 的统计特性依赖足够的样本数经验上每阶待估模态至少需要覆盖 20 个以上振动周期。数据只有几秒钟而最低模态只有 0.5 Hz 时p 和 q 取到 30 阶以上容易产生虚假极点。这时优先缩减阶次把 p 设定在 2 倍模态数左右再用残差白噪声检验判断模型是否充分提取了系统信息。5.4 用稳定图区分物理模态与数学极点ARMA 识别出的极点数永远等于 AR 阶数其中混杂着大量用来拟合噪声的数学极点。区分方法是逐阶扫描 p把每一阶下识别出的频率和阻尼比画在同一张图上真实物理模态的频率和阻尼会随着阶次升高保持稳定数学极点则四处漂移。稳定判据一般取频率变化小于 1%、阻尼比变化小于 5% 或绝对差小于 0.005。p_list range(10, 50, 2) stable [] for p in p_list: a, _, _ arma_fit(y, p, p, pre_order40) for freq, zeta, _ in arma_modes(a, fs100.0, f_max20.0): stable.append((p, freq, zeta))输出后以 p 为横轴、频率为纵轴散点绘图垂直方向几乎排成一条直线的点带就是可信的物理模态。阻尼比散乱但频率稳定的极点通常是真实存在的弱激励模态频率和阻尼都稳定的极点则是整条识别流程中最值得留用的结果。本文还有配套的精品资源点击获取