简介：这份PDF资料聚焦ARMA模型时间序列分析法在模态参数识别中的应用，面向结构动力学、振动测试与信号处理方向的学习者与工程技术人员，帮助读者从有序随机振动响应数据中提取自然频率、阻尼比与振型等动态特性。内容共5页，系统讲解AR模型与MA模型的组成逻辑、ARMA时序模型方程、脉冲响应函数与相关函数推导，并给出推广Yule-Walker方程、伪逆法最小二乘求解自回归系数、Newton-Raphson迭代估算滑动平均系数，以及由传递函数极点反推模态频率与阻尼比的完整公式链条，还涉及留数与归一化复振型向量的计算思路。资源包为1个PDF文件，约199KB，轻量便于随时查阅。目前已有819人学习，适合希望快速掌握时序分析法原理推导与模态参数识别流程的读者作为公式速查与入门参考。

1. 从振动数据到模态参数：这份 5 页 PDF 到底解决了什么问题

做结构模态测试的人多半遇到过这种场景：锤击法或者激振器试验采回来一堆响应数据，频响函数曲线看着挺漂亮，可一到密集模态或者环境激励（只有响应没有激励）的工况，峰值拾取法就开始翻车——两个挨得很近的峰根本分不开，阻尼比更是估得离谱。这份《ARMA模型时间序列分析法 时序分析法 模态参数识别的方法 原理讲解 公式推导 共5页.pdf》讲的，就是绕开频域峰值拾取、直接在时域里用参数模型把模态参数抠出来的路子。它的核心思路是：把有序随机振动响应数据看作一个 ARMA 过程，用差分方程去拟合这段序列，再从拟合出的传递函数极点里反解出模态频率和阻尼比。适合已经懂一点结构动力学、手头有实测响应数据、想补上时域识别这一环的工程师；如果你连自功率谱和互功率谱都还没算过，建议先把频域基础打牢再回来啃这 5 页。

2. ARMA 时序模型的数学骨架：差分方程、Yule-walker 与传递函数

2.1 从微分方程到差分方程：AR 和 MA 各自管什么

N 个自由度的线性系统，激励与响应之间本来是连续时间域的高阶微分方程。到了离散时间域，微分变成差分，响应序列 $x_t$ 就被写成当前值和历史值、以及白噪声激励历史值的线性组合：

$$x_t = \sum_{k=1}^{2N} a_k x_{t-k} + \sum_{k=0}^{2N} b_k f_{t-k}$$

等号左边那串 $\sum a_k x_{t-k}$ 是自回归部分（AR），它管的是“当前响应和过去响应之间的关系”，本质上是系统自身惯性和弹性的记忆效应；右边 $\sum b_k f_{t-k}$ 是滑动平均部分（MA），它管的是“外部白噪声激励经过系统之后留下的痕迹”。$2N$ 是阶次，$a_k$、$b_k$ 是待识别系数，$f_t$ 是白噪声。当 $k=0$ 时约定 $a_0=b_0=1$。这里有个容易看漏的点：ARMA 里的阶次用的是 $2N$ 而不是 $N$，因为一个 N 自由度系统的特征方程是 N 阶的，写成差分方程后对应 2N 个系数，后面解极点时也是从 2N 次代数方程里出根。

为什么非要拆成 AR 和 MA 两块？只留 AR 的话，模型等价于把响应看成自身历史的线性外推，对宽带白噪声激励下的响应拟合没问题，但对激励本身有色、或者测量噪声混进来的情况就偏；只留 MA 的话，参数估计会变成非线性问题，迭代容易发散。ARMA 把线性部分（AR）和非线性部分（MA）分开处理，先线性解 AR，再非线性解 MA，这是它能落地工程的关键。

2.2 相关函数与 Yule-walker 方程：把参数估计变成线性代数

直接对差分方程做最小二乘是行不通的，因为 $f_t$ 是未知的白噪声。这份 PDF 走的是相关函数路线：先算响应序列的自相关函数 $R_\tau$，利用白噪声自相关只在 $\tau=0$ 处有值（方差 $\sigma^2$）这个性质，把式 (3) 代入后得到：

$$R_\tau = \sigma^2 \sum_{i=0}^{\infty} h_i h_{i+\tau}$$

再结合脉冲响应函数满足的差分关系，当滞后 $l > 2N$ 时 MA 系数 $b_k$ 全部为零，于是得到一组只含 AR 系数 $a_k$ 的方程：

$$R_l = \sum_{k=1}^{2N} a_k R_{l-k}, \quad l > 2N$$

把不同的 $l$ 值代进去，就凑成推广的 Yule-walker 方程。写成矩阵形式是 $[R]{a} = {R'}$，其中 $[R]$ 是自相关矩阵，${R'}$ 是右端向量。因为相关函数长度 $L$ 通常远大于 $2N$，方程个数多于未知数个数，属于超定方程组，用伪逆法求最小二乘解：

$${a} = ([R]^T[R])^{-1}[R]^T{R'}$$

这一步是整个流程里最“线性”的部分，也是最好写代码的部分。AR 系数解出来之后，MA 系数 $b_k$ 要靠非线性方程组 (14) 来求，PDF 里提到两类方法：基于 Newton-Raphson 的迭代最优化，和基于最小二乘原理的次最优化。工程上我一般先用次最优方法拿初值，再用 Newton-Raphson 精修，直接上迭代很容易因为初值太远而发散。

2.3 从传递函数极点到模态频率与阻尼比

AR 和 MA 系数都拿到之后，ARMA 模型的传递函数是：

$$H(z) = \frac{\sum_{k=0}^{2N} b_k z^{-k}}{\sum_{k=0}^{2N} a_k z^{-k}}$$

分母多项式等于零就是特征方程：

$$\sum_{k=0}^{2N} a_k z^{-k} = 0 \Rightarrow z^{2N} + a_1 z^{2N-1} + \cdots + a_{2N} = 0$$

解这个 2N 次代数方程，得到 2N 个根 $z_k$，它们就是传递函数的极点。极点一般是共轭成对出现的，每一对对应一阶模态。极点和模态参数的换算关系是：

$$z_k = \exp(s_k \Delta t), \quad s_k = -\xi_k \omega_k \pm j\omega_k\sqrt{1-\xi_k^2}$$

从极点反解模态频率和阻尼比：

$$\omega_k = \frac{|\ln z_k|}{\Delta t}, \quad \xi_k = \frac{-\text{Re}(\ln z_k)}{|\ln z_k|}$$

这里 $\Delta t$ 是采样间隔。实际写代码时要注意：$\ln z_k$ 是多值的，取主值就行，因为采样定理保证了 $|\text{Im}(\ln z_k)| < \pi$。算完频率和阻尼比，振型还得靠留数。留数 $A_{pqk}$ 用极点处的极限求：

$$A_{pqk} = \lim_{z \to z_k} (z - z_k) H_{pq}(z) = \frac{\sum b_k z_k^{-k}}{\prod_{i \neq k}(z_k - z_i)}$$

对同一阶模态，把 n 个测点的留数都求出来，找绝对值最大的那个测点作为参考，归一化之后就得到复振型向量。这一步是 ARMA 法能出振型的关键，也是很多人只算频率阻尼、不算振型的原因——留数对噪声比极点敏感得多。

3. 把公式落成代码：从响应数据到模态参数的完整实现

3.1 数据预处理与自相关函数估计

拿到实测响应数据，第一步不是直接套公式，而是去均值、去趋势。振动信号里如果混了直流分量或者温度漂移，自相关函数在 $\tau$ 大时会翘起来，Yule-walker 方程直接解歪。常见做法是先做一阶差分或者多项式拟合去趋势，再用 Welch 法估计自相关。

import numpy as np from scipy.signal import detrend, correlate def estimate_autocorr(x, max_lag): """ 估计响应序列的自相关函数 x: 一维响应数组 max_lag: 最大滞后点数，一般取 2N 的 3~5 倍 """ x = detrend(x, type='linear') # 去线性趋势 x = x - np.mean(x) # 去均值 n = len(x) # 用 FFT 加速自相关，比直接循环快一个量级 acf = correlate(x, x, mode='full')[n-1:] / n return acf[:max_lag+1]

这段代码里detrend用线性去趋势，是因为实测数据里最常见的干扰就是传感器温漂带来的慢变基线。correlate用mode='full'之后取后半段，得到的是无偏估计的近似。max_lag取 2N 的 3 到 5 倍，是因为 Yule-walker 方程要用到 $l > 2N$ 的滞后，留够余量才能让最小二乘稳定。如果数据长度不够，自相关尾部噪声大，解出来的 AR 系数会飘。

3.2 用伪逆法解 Yule-walker 方程求 AR 系数

有了自相关序列，构造自相关矩阵和右端向量，直接上伪逆。这里的关键是矩阵的构造方式：第 $i$ 行第 $j$ 列的元素是 $R_{M+i-j}$，其中 $M=2N$。

def solve_ar_coeff(acf, order): """ 解 Yule-walker 方程求 AR 系数 acf: 自相关序列，长度至少 2*order+1 order: 模型阶次 2N """ M = order L = len(acf) - 1 # 构造自相关矩阵 R，尺寸 (L-M) x M R = np.zeros((L - M, M)) for i in range(L - M): for j in range(M): R[i, j] = acf[M + i - j] # 右端向量 r = acf[M+1:L+1] # 伪逆法最小二乘解 a = np.linalg.pinv(R) @ r return a

np.linalg.pinv内部走的是 SVD，比直接求逆稳，因为自相关矩阵在滞后大时接近奇异。order就是 PDF 里的 $2N$，如果你有 3 个自由度，order 取 6。L是自相关序列长度，一般取 order 的 5 到 10 倍。这里有个血泪经验：如果R的条件数超过 $10^{12}$，解出来的 $a$ 会完全不可信，这时候要么降阶，要么加正则化项。

3.3 MA 系数求解与极点提取

MA 系数求解是非线性的，PDF 里给了 Newton-Raphson 和次最优两条路。工程上我一般先用次最优方法：把式 (14) 的非线性方程组在初值附近线性化，迭代几步拿到粗略的 $b_k$，再用 Newton-Raphson 精修。

from scipy.optimize import fsolve def solve_ma_coeff(acf, a_coeff, order): """ 解 MA 系数，用 fsolve 做非线性方程组求解 acf: 自相关序列 a_coeff: 已求出的 AR 系数 order: 2N """ M = order # 计算中间量 C_k = sum_i sum_j a_i a_j R_{k+i-j} def compute_C(k): s = 0.0 for i in range(M+1): for j in range(M+1): ai = 1.0 if i == 0 else a_coeff[i-1] aj = 1.0 if j == 0 else a_coeff[j-1] s += ai * aj * acf[abs(k + i - j)] return s def equations(b): eqs = [] for k in range(M+1): lhs = sum(b[i] * b[i+k] for i in range(M+1-k)) eqs.append(lhs - compute_C(k)) return eqs b0 = np.ones(M+1) * 0.1 # 初值给小的正数 b_sol = fsolve(equations, b0) return b_sol

fsolve默认用 MINPACK 的 hybrd 算法，本质就是拟牛顿法。初值给 0.1 而不是 0，是因为 $b_0$ 在式 (14) 里出现在分母位置，给 0 会直接除零。compute_C里的双重循环在 order 不大时（一般不超过 20）完全够用，order 再大就得改成矩阵运算。解出 $b$ 之后，构造分母多项式系数，用np.roots求根：

def extract_modal_params(a_coeff, dt): """ 从 AR 系数提取模态频率和阻尼比 a_coeff: AR 系数，长度 2N dt: 采样间隔 """ # 分母多项式：z^{2N} + a1 z^{2N-1} + ... + a_{2N} coeffs = np.concatenate([[1.0], a_coeff]) poles = np.roots(coeffs) # 只取模大于 1 的根（因果系统极点应在单位圆外） poles = poles[np.abs(poles) > 1.0] modal_params = [] for z in poles: ln_z = np.log(z) omega = np.abs(ln_z) / dt xi = -np.real(ln_z) / np.abs(ln_z) modal_params.append((omega, xi)) return modal_params

np.roots对 2N 次多项式用的是伴随矩阵特征值法，数值稳定性比直接求根公式好。极点筛选那一步很关键：理论上因果系统的极点应该在单位圆外（$|z|>1$），但实测数据里总有几个根落在圆内，那是噪声或者数值误差产生的虚假模态，直接扔掉。算出来的 $\omega$ 是圆频率，除以 $2\pi$ 才是赫兹。

4. 避坑与排查：ARMA 模态识别里最容易翻车的五个地方

4.1 现象：解出的阻尼比是负数或者大得离谱

原因：极点位置对 AR 系数的误差极其敏感，而 AR 系数又是从自相关矩阵的最小二乘解里来的。如果自相关序列尾部噪声大，或者矩阵条件数太高，解出来的 $a_k$ 会有微小扰动，反映到极点上是实部符号翻转，阻尼比就成负的。另一个常见原因是采样频率选得太高，模态频率对应的归一化频率接近 0 或 0.5，极点挤在一起分不开。

解决：先检查自相关序列在最大滞后处是否已经衰减到接近零，如果还在振荡说明数据里有未去除的周期成分。把采样频率降到模态最高频率的 5 到 10 倍，别盲目追求高采样率。解出极点后，对阻尼比做物理约束：$\xi$ 在 0 到 0.2 之间是结构模态的常见范围，超出这个范围的根直接标记为可疑。

4.2 现象：阶次选高了出现一堆虚假模态，选低了真实模态被吞掉

原因：ARMA 模型的阶次 $2N$ 需要事先知道或者估计。实际结构自由度是连续的，离散成 N 个模态只是近似，阶次选不对，模型要么过拟合噪声，要么欠拟合真实动态。

解决：用 AIC 或者 BIC 准则扫一遍阶次。AIC 的公式是 $\text{AIC} = \ln(\hat{\sigma}^2) + 2p/N$，其中 $p$ 是参数个数，$N$ 是数据长度。具体做法是从低阶往高阶扫，画 AIC 随阶次变化的曲线，取曲线拐点。我一般还会配合稳定图：同一个阶次下，频率和阻尼比随模型阶次变化很小的极点才认为是真实模态。

4.3 现象：MA 系数迭代不收敛，fsolve报“迭代次数超限”

原因：式 (14) 的非线性方程组对初值很敏感，如果初值离真解太远，Newton-Raphson 会发散。另一个原因是 $b_0$ 的约束没有加进去，解出来的 $b$ 不满足 $b_0=1$ 的约定。

解决：先用次最优方法拿初值——把式 (14) 在 $b_k$ 的零附近做一阶泰勒展开，解一个线性最小二乘问题，得到的解作为fsolve的初值。如果还不行，把fsolve的xtol放宽到 $10^{-6}$，maxfev加到 5000。实在不收敛就退回纯 AR 模型，虽然 MA 部分丢了，但频率估计通常还能用。

4.4 现象：算出来的频率和频域峰值拾取对不上，差了百分之几

原因：ARMA 法识别的是离散时间模型的极点，频域峰值拾取找的是功率谱的局部极大值。两者在阻尼小的时候应该一致，但阻尼大或者模态密集时，功率谱峰值会被相邻模态“拉偏”，而 ARMA 极点不受这个影响。所以对不上不一定是 ARMA 错了，可能是频域方法偏了。

解决：用半功率带宽法从频响函数上单独估一个阻尼比，和 ARMA 的结果对比。如果 ARMA 的频率落在两个频域峰之间，而且阻尼比更合理，那大概率是频域方法分辨不开。反过来，如果 ARMA 频率跑到频域峰外面很远，先检查采样间隔 $\Delta t$ 有没有代错。

4.5 现象：振型算出来相位乱跳，归一化之后符号对不上

原因：留数对噪声的敏感度比极点高一个量级，尤其是响应测点信噪比低的时候，留数的相位误差会直接传到振型上。另外，归一化时选的参考测点如果正好在某阶模态的节点附近，留数绝对值很小，归一化会放大误差。

解决：算留数之前先对响应数据做带通滤波，只保留目标模态附近的频带。参考测点不要选节点位置，可以选留数绝对值最大的测点，或者干脆用所有测点的留数做整体最小二乘拟合。振型符号统一用参考测点的相位做基准，相位差超过 90 度就翻转。

5. 进阶技巧：用稳定图和留数拟合把识别结果钉死

把 ARMA 跑通只是第一步，真正让结果可信的是后处理。我一般会做两件事：稳定图和留数整体拟合。

稳定图的做法是：把模型阶次从 2 扫到 20，每个阶次都跑一遍 ARMA，把识别出的频率和阻尼比画在同一张图上。真实模态的特征是——随着阶次增加，频率和阻尼比基本不变，形成一条竖直的“稳定轴”；虚假模态则到处乱跳。判断稳定的阈值我一般设频率变化小于 1%、阻尼比变化小于 5%。这一步能把 4.2 里说的阶次选择问题直接可视化。

留数整体拟合是针对振型的。单个测点的留数误差大，但同一阶模态下所有测点的留数应该满足同一个传递函数结构。把所有测点的留数堆成一个矩阵，做一次全局最小二乘：

def global_residue_fit(residues, poles, mode_idx): """ 对同一阶模态的所有测点留数做整体拟合 residues: 形状 (n_points, n_modes) 的留数矩阵 poles: 极点数组 mode_idx: 目标模态索引 """ zk = poles[mode_idx] # 每个测点的留数除以参考测点留数，得到归一化振型 ref_idx = np.argmax(np.abs(residues[:, mode_idx])) phi = residues[:, mode_idx] / residues[ref_idx, mode_idx] # 用相位一致性做符号校正 phase_ref = np.angle(phi[ref_idx]) for i in range(len(phi)): if np.abs(np.angle(phi[i]) - phase_ref) > np.pi/2: phi[i] = -phi[i] return phi

ref_idx选留数绝对值最大的测点，避免选到节点。相位校正那一步是必须的，因为np.angle返回的是主值，相邻测点的相位差如果超过 180 度会被折叠，导致振型符号跳变。做完这一步，振型向量就可以直接拿去和有限元结果做 MAC 对比了。

从那以后我每次跑 ARMA 模态识别，都强制走一遍“去趋势 → 自相关检查 → 阶次扫描 → 稳定图 → 留数整体拟合”这五步，少一步结果都不敢往外发。希望帮到你。

本文还有配套的精品资源，点击获取