SSI-COV协方差驱动随机子空间识别：原理、Matlab实现与调参实战

做结构模态测试的人，手里如果已经有几组加速度响应数据，又不想被频域方法的各种窗函数和平均次数搞得心烦，那么SSI-COV（协方差驱动随机子空间识别）是一个非常值得掌握的工具。它直接用环境激励下的响应数据来识别模态频率、阻尼比和振型，不需要人工激励，也不需要知道激励的具体大小。这篇文章就从工程应用的角度，讲讲这套方法的原理、Matlab实现细节、以及我实际调参数踩过的坑。

1. 背景：为什么选SSI-COV而不是频域方法

先聊一个很多人问过的问题：都有了频响函数和峰值拾取法，为什么还要用SSI-COV？原因在于，实际结构测试里，激励往往不是可控的，例如桥梁上的车辆荷载、建筑物上的风荷载、机械运行时的环境振动，这些都是随机激励，输入无法准确测量。频域方法需要把激励和响应做互谱或频响函数，但环境激励下很难拿到高质量的激励信号，强行用峰值法只能看到共振峰，阻尼比算出来误差也大，振型还可能被密集模态污染。

SSI-COV走的是另一条路，直接用响应的协方差序列来构造系统状态空间模型，再把模型转化成模态参数。它的名字里“协方差驱动”指的就是这一步。它的好处非常实际：

• 不需要已知输入，环境激励即可，利用随机激励下输出响应的统计特性
• 能同时识别频率、阻尼比、振型，并且不像FDD那样靠峰值猜测模态阶次
• 抗噪声干扰能力较强，尤其对低频、弱模态比频域法稳定
• 有成熟的状态空间理论基础，模态参数的不确定性可以做后验评估

当然它也有代价。最直观的就是计算量比频域方法高，而且需要人为确定系统阶次。阶次选不好，就很容易出现虚假模态。但这些问题通过稳定性图可以部分缓解，我在后文会专门讲。

我在参与的一个人行天桥振动测试项目里，实测数据只有12个通道的加速度时程，采样率256Hz，总时长20分钟。用PolyMAX和SSI-COV都跑了一遍，SSI-COV识别出的前四阶阻尼比虽然只有0.6%到2.3%，但重复测试间的方差明显比频域方法小，尤其第二阶和第三阶两个频率相距只有3Hz左右的模态，SSI-COV能把振型分开，峰值法则几乎糊成一团。那个项目之后，我就把SSI-COV放进了常规工具箱。

2. 核心原理：状态空间模型和协方差序列的来龙去脉

SSI-COV的理论基础是线性时不变系统的状态空间描述。一个多自由度结构在离散时间下的状态方程和观测方程可以写成：

x(k+1) = A * x(k) + w(k) y(k) = C * x(k) + v(k)

这里的x是状态向量，包含位移和速度（或者对应的离散状态变量），y是实测输出，A是系统矩阵，C是观测矩阵，w和v分别是过程噪声和测量噪声。我们做模态识别的目标，就是从实测输出y(k)中估计系统矩阵A的特征参数，再从中解出振动系统的固有频率、阻尼比和振型。

这里有个关键点随机子空间的假设：w和v都是零均值白噪声，并且与系统状态不相关。这是一个很强的假设，但工程上大多数环境振动激励逼近这个条件。如果数据里有强谐波干扰或者明显的非平稳漂移，请看第5节预处理的方法。

协方差驱动的思路是这样的：先定义输出协方差序列R(i) = E[y(k+i) * y(k)^T]，然后用这些协方差块来构造一个分块Toeplitz矩阵。为什么可以这么做？因为线性系统本身有马尔可夫性，未来输出和过去输入的统计相关性中包含了系统动态的全部信息。具体推导过程不展开，最终结论是Toeplitz矩阵可以分解为可观矩阵和可控矩阵的乘积，也就是：

T = O * G

其中O是观测矩阵和系统矩阵构成的可观矩阵，G是可控矩阵。对T做奇异值分解（SVD），就能得到系统矩阵A的估计。得到A之后，对它做特征值分解，复特征值对应系统的极点，利用极点和采样时间可以换算频率和阻尼比，而观测矩阵C乘以特征向量就得到振型。

我知道很多第一次看到这部分的人会被这些矩阵绕晕，但用一句话总结就是：协方差序列把激励的影响“平均”掉了，剩下的就是系统本身的状态演化规律。SVD把噪声子空间和信号子空间分开，取信号子空间就可以提取模态参数。这就是SSI-COV的全部骨架，其他都是对这个骨架的完善和数值优化。

3. 数据类型和预处理：直接决定识别质量的三个步骤

很多人一开始做SSI-COV识别出来一塌糊涂，第一反应是算法不行，其实八成是数据预处理没做好。数据预处理对SSI-COV的影响，比对频域方法的影响还要大。因为协方差序列对趋势项、直流偏移、高频噪声都非常敏感，任何一个不当处理都会在后续的SVD中变成虚假模态。我一般按下面三步走。

第一步是去趋势和去均值。加速度传感器如果有轻微的温漂或零漂，数据里就会有一个时变趋势，这个趋势在计算协方差时相当于一个低频强分量，会在低阶模态附近造成假峰。建议先用多项式拟合并去除趋势项，同时减去均值。但多项式阶数不用太高，通常1到2阶就够了，太高反而会吃掉真实的低频模态。

第二步是低通滤波。原始数据里高频噪声越少，需要的系统阶次就越低，稳定图上的假模态也越少。不过滤波会改变相位，对频域方法影响很大，但SSI-COV是基于统计特性的，对相位不是特别敏感，所以可以在滤波后直接使用。关键是滤波器要选择零相位版本，例如Matlab的filtfilt函数，可以避免相位偏移。截止频率一般取你关心的最高分析频率的1.5倍左右，例如关心到50Hz，就滤到75Hz，再高只会增加噪声。

第三步是重采样和降采样。这一步常被忽略。SSI-COV的计算量和采样点数直接相关，长数据+高采样率会让汉克尔矩阵和Toeplitz矩阵的维度爆炸，没必要。重采样可以在低通滤波后进行，采样率降到分析最高频率的2.5到3倍即可。例如原采样率1024Hz，关心最高频率30Hz，完全可以直接降采样到100Hz，计算量减少90%，识别结果基本不变。我经常用resample函数，它会自动做抗混叠滤波，比我手动滤波后再抽取靠谱。

预处理之后，最好把数据做一下可视化检查，看所有通道的时程量级是否在一个数量级。如果某个通道的均方根值比其他通道小几十倍，先查传感器，不是每个通道都强制参与识别。在SSI-COV里，通道量级差异过大会让SVD的奇异值偏向大通道，小通道的振型信息容易被淹没。必要时对每个通道做归一化，但振型的绝对值就会丢失，只能从相对值角度使用，所以不是万不得已我不做归一化，而是通过传感器标定和放大倍数设置让所有通道基本一致。

4. Matlab代码实现：从协方差矩阵到频率阻尼振型

下面给出我在工程中使用的SSI-COV主程序框架，代码通俗易懂，核心步骤都加了注释。

4.1 主流程代码框架

function [fn, zeta, phi] = ssi_cov(y, fs, ncols, nrows, order) % y : 输出响应矩阵，每一列为一个测点通道 % fs : 采样频率 % ncols : Toeplitz矩阵的列分块数 % nrows : Toeplitz矩阵的行分块数 % order : 系统阶次（偶数） [N, nch] = size(y); % 计算协方差序列 R(i)，i = 0, 1, ..., nrows+ncols-2 maxlag = nrows + ncols - 2; y_mean = mean(y, 1); y = y - y_mean; R = zeros(nch, nch, maxlag+1); for k = 0:maxlag % 对应样本协方差，滞后k a = y(1:N-k, :); b = y(1+k:N, :); R(:,:,k+1) = (a' * b) / (N-k); end % 构造分块Toeplitz矩阵 T = zeros(nch*nrows, nch*ncols); for i = 1:nrows for j = 1:ncols lag = ncols - j + i - 1; % 需要仔细推导的索引 T((i-1)*nch+1:i*nch, (j-1)*nch+1:j*nch) = R(:,:,lag+1); end end % SVD分解 [U, S, V] = svd(T, 'econ'); % 根据系统阶次order截断 order2 = order / 2; U1 = U(:, 1:order); S1 = S(1:order, 1:order); V1 = V(:, 1:order); % 计算状态矩阵A的估计 O = U1 * sqrt(S1); % 可观测矩阵估计 O_bar = O(1:end-nch, :); % 删掉最后nch行 O_up = O(nch+1:end, :); % 删掉最前nch行 A_est = O_bar \ O_up; % 观测矩阵C估计：取可观测矩阵的前nch行 C_est = O(1:nch, :); % 特征值分解 [Psi, Lambda] = eig(A_est); lambda = diag(Lambda); % 离散极点转连续极点 mu = log(lambda) * fs; % 频率和阻尼比 fn = abs(mu) / (2*pi); zeta = -real(mu) ./ abs(mu) * 100; % 百分比% % 振型 phi = C_est * Psi; % 每个特征向量的列对应一个模态 end

这段代码可以跑通，但我必须提醒几点。

4.2 关于索引和矩阵维度的几个大坑

第一，Toeplitz矩阵的索引是最容易写错的。我上面的代码用的是lag = ncols - j + i - 1，这个公式取决于R(1)对应滞后0。一旦把滞后次序搞反，可能识别出的频率是正确的，但阻尼比符号是反的，振型也会乱。我自己的建议是不要凭记忆写这段，用一个滞后矩阵从左到右、从上到下打出来检查一遍。

第二，系统阶次order必须为偶数，因为每一阶物理模态对应一对共轭复极点。如果是奇数，特征值分解后会出现一个实数极点，对应的“模态”不是物理意义的共振，后续很容易被当成虚假模态误删。

第三，A_est的求解用O_bar \ O_up，Matlab默认的最小二乘。这一步如果O的条件数很差，结果就会不稳。此时可以尝试在SVD截断前进一步做奇异值截断。S值如果出现跳崖式下降，说明截断阶次应该选在跳崖点前的平坦区域。下面我专门讲怎么定阶。

5. 系统阶次怎么定：奇异值曲线和稳定图实战

SSI-COV最让人纠结的就是order取值。order定太小，模态会漏掉，阻尼比估计偏差也大；order定太大，稳定图上全是计算产生的数值极点，筛选起来头痛。我的做法是两层筛选结合。

第一层看奇异值曲线。SVD得到的奇异值S对角元从大到小排列，信号子空间的奇异值远大于噪声子空间的奇异值，所以在对数坐标下，奇异值曲线从陡降变成平缓的转折点就是有效阶次的参考位置。转折点之前的数量乘以2，大致对应系统阶次上限。实际使用时，我在那个点附近再扩大一倍留出余量，比如转折点在25，就可以尝试order在30到50之间变化。

第二层用稳定图。Stabilization Chart我在实际项目里必做，步骤很简单：

• 设定一个order序列，比如从10到80，步长为2
• 对每个order跑一次SSI-COV，得到一组频率、阻尼比、振型
• 按相似容差判断哪些模态在多次识别中稳定

稳定性的标准一般这样设置：

% 频率容差 1% f_tol = 0.01; % 阻尼比容差 10%（阻尼本身识别难度大，可放宽） d_tol = 0.10; % MAC值容差 2% mac_tol = 0.02;

如果某个模态候选在相邻两次order识别中，频率变化小于1%，阻尼比变化小于10%（相对值），MAC大于98%，就可以认定为稳定点。把所有稳定点画在“频率-order”平面上，形成竖线，竖线聚集处就是真实模态。稳定图看起来直观，但筛得很科学。

值得注意的是，阻尼比容差不能太紧。阻尼比本身是识别量里噪声最大的一个，尤其环境激励下，2%的阻尼比和2.2%的阻尼比，在物理上可能没区别，但相对变化超过10%容易把稳定点漏掉。我通常把阻尼的稳定性判据放宽到20%，频率还是1%。

下面给一个画稳定图的常见代码段：

orders = 10:2:60; all_fn = cell(length(orders), 1); all_zeta = cell(length(orders), 1); all_phi = cell(length(orders), 1); for oi = 1:length(orders) [fn_o, zeta_o, phi_o] = ssi_cov(y, fs, 20, 20, orders(oi)); % 只保留正频率且阻尼比在0到10%之间的假想模态 valid = fn_o > 0 & zeta_o > 0 & zeta_o < 10; all_fn{oi} = fn_o(valid); all_zeta{oi} = zeta_o(valid); all_phi{oi} = phi_o(:, valid); end % 循环配对并标记稳定点 % 具体实现可根据数据量优化，这里只给出思虑框架

稳定图的绘制代码不难，但耗时在配对逻辑上。如果每个order识别出20个候选模态，60个order就是1200次配对，用for循环也能跑，但数据量大时要向量化。我在代码里会先用频率排序做预筛，再用MAC矩阵一次算完，这样能省很多时间。

6. 实测踩坑记录：为什么识别的阻尼老是漂

阻尼比对SSI-COV来说是个难点。我整理了自己项目里常见的问题，列成一张表，大家可以直接对照排查。

还有一个我一开始没注意的问题是：数据长度对阻尼识别影响极大。SSI-COV对协方差的估计质量依赖于数据长度。理论上，协方差序列R(i)需要足够多的样本才能收敛。如果数据长度只有几十秒，低频模态的协方差还没被充分平均，阻尼比的估计会偏向于随机游走。我现在的经验是最低确保每个分析频带内至少500个循环周期。比如关注最低频率1Hz，至少采集500秒；如果最低频率0.1Hz，至少采集5000秒。这在实际桥梁测试中往往要熬夜挂机。

另一个容易出错的是重采样导致的虚假模态。我遇到过识别出一个12.5Hz的模态，怎么调参数都不消失，后来发现是原始数据里面有50Hz市电工频干扰，resample到100Hz后混叠到了12.5Hz附近。解决方案是先做50Hz陷波滤波，或者先用零相位低通滤到远低于奈奎斯特频率后再降采样。这个例子充分说明，预处理做的多细致，后续就能多省心。

另外，如果要识别多参考点振型，并希望振型归一化到某一点，建议在输入SSI-COV之前就记录好测点坐标和通道对应关系。识别之后，振型相位是复数域的信息，土木结构小阻尼情况下相位接近0°或180°，但如果识别出的复数振型相位在空间上不是连续分布，往往不是结构问题，而是传感器方向标错了。我在一次试验中吃过亏，后来养成了识别完成后画振型动画的习惯，比单纯看数字可靠得多。

7. 自动筛选真实模态：MAC和模态参与因子的使用

稳定图帮我们缩小了范围，但真实工程数据里，依然会有一些“看起来稳定”但物理上说不通的极点，比如数值模态和真实模态频率太接近，或者谐波分量近似周期信号造成的伪极点。这时候需要另外两个工具：MAC（Modal Assurance Criterion）和模态参与因子。

MAC是判断两个振型向量相关程度的工具，公式是：

MAC_ij = abs(phi_i' * phi_j)^2 / ((phi_i'*phi_i) * (phi_j'*phi_j));

当两个振型来自不同结构模态时，MAC值理论上接近0（数值上一般小于0.2），来自同一模态时接近1。在自动筛选时，我会这样做：

• 计算所有候选模态之间的MAC矩阵
• 把MAC超过0.9且频率也接近的模态合并为一类
• 每一类里选稳定图上出现次数最多的那个频率作为最终模态

这个步骤可以有效去除重复识别。但MAC也有局限，如果传感器数量少或者测点分布不理想，不同振型之间的MAC可能天然很高。比如只有两个测点且对称布置，反对称模态和对称模态在测点上就可能分不开。此时只能靠频率差和阻尼比经验判断，没法完全自动化。

模态参与因子主要用来识别谐波分量。环境激励下的旋转机械或桥梁吊杆，经常会产生近似正弦的窄带激励。这种激励在SSI-COV里会使某个极点的“参与因子”特别大，而且这个极点对应频率的奇异值也很大。参与因子可以从可控矩阵的估计值里算出来，但由于代码篇幅原因不展开。实际项目中，我会在目标频率附近留出一段“警戒区”，把任何频率落在电网工频、转速基频附近的候选模态都标记为待验证，配合现场工况记录判断是否剔除。

8. 实用工具箱推荐和代码调优思路

如果不想从零开始写SSI-COV，Matlab里有几个现成的工具箱可供参考，我在项目中也混用过：

• MACEC：比利时鲁汶大学开发的结构模态分析工具箱，里面有SSI-COV和稳定图实现，学术圈用得很多，可以参考源码学习细节。
• MATLAB System Identification Toolbox：n4sid和ssest函数可以实现部分子空间辨识，但面向一般线性系统，输出的是状态空间模型，需要自己做模态解算。
• 一些GitHub上的开源实现：搜索cov-ssi或stochastic-subspace，质量参差不齐，用之前一定要检查Toeplitz索引和奇异值截断逻辑。

我的建议是第一次做这个研究的人，先自己写一个简化版，比如只识别单自由度系统的模态参数，把每个矩阵打印出来逐步理解。这样比直接拿工具箱跑真实数据，学习效率高得多。理解原理之后再考虑封装成自己的函数库，因为实际项目的通道数、采样频率、阶次范围差异很大，开源代码不一定能直接复用。

调优时有几个思路值得尝试：

1. 分频带处理。如果关心0.1Hz到100Hz很宽的频带，一次做完往往效果不好。可以先低通滤到30Hz识别低频段，再把高通滤掉低频部分识别高频段，但要注意滤波造成的边界效应，数据首尾要预留足够长度的丢弃区。
2. 多段数据平均。如果现场只记录了两次工况，可以把两次数据拼接后进行整体识别吗？不建议直接拼接，因为两次数据的环境激励水平和结构状态可能不同。更好的做法是分别识别，再将频率和振型结果做加权平均，权重和信噪比相关。
3. 通道加权。如果测点较多，某些通道信噪比明显偏低，可以在构建协方差矩阵之前对数据做加权，但加权会影响振型的绝对值，因此我只在严重不平衡时使用，而且识别完会把振型重新缩放到实际测点单位。

工具箱不是重点，重点是你能不能用里面的函数解释清楚每一个输入参数的含义。曾经有朋友问我为什么他用了某个现成代码后阻尼比每天跑结果都不一样，我让他看了代码里对协方差滞后数的选择。那一版代码里ncols和nrows写死了，根本不随数据长度变化，低频部分协方差估计方差极大，阻尼比自然天天变。这种问题只有自己掌握了原理才看得懂。

9. 一次真实数据的完整识别流程复盘

我用一个人行天桥案的例子，把整套流程串起来讲，方便大家照葫芦画瓢。

那次测试一共布了8个加速度传感器，分布在主梁第二跨的八等分点上，采样率200Hz，连续记录20分钟。我按下面的流程走了一遍。

数据预处理：

fs = 200; % 原始数据矩阵 y_raw，每列一个通道 y_raw = detrend(y_raw, 1); % 线性去趋势 y_raw = y_raw - mean(y_raw, 1); % 设计零相位低通，截止40Hz [b, a] = butter(4, 40/(fs/2), 'low'); y_fit = filtfilt(b, a, y_raw); % 降采样到100Hz fs_new = 100; y_new = resample(y_fit, fs_new, fs);

参数选择上，nrows和ncols我取20。注意，这个参数对应Toeplitz矩阵里用到的协方差滞后数，不是越高越好。理论上滞后阶数越大，包含的模态信息越多，但协方差估计误差也随滞后增大而增大。常看到的建议是ncols取最大滞后对应10到30个采样周期，覆盖最关心模态的3个周期以上。以最低模态频率1Hz为例，100Hz采样下对应100个采样点为1周期，ncols=20对应滞后到2s，足够覆盖。但如果是0.2Hz的超低频结构，ncols=20只覆盖1s，远远不够，必须加大。这个关系很好理解。

之后对order从20到80，步长2跑了稳定图，发现三列很明确的稳定线，分别对应约1.9Hz、4.6Hz、8.2Hz。第四阶8.2Hz在较大order时才稳定，阻尼比约1.1%。最终以order=42的频率结果作为最终输出，因为该阶次下所有目标模态都稳定，且没有多余假极点。

振型结果输出前，我把8个测点的相对振幅按测点坐标做了插值，画成侧视图，确认第二阶振型在中跨出现反弯点，和有限元预测一致。这一步不可省，因为只有通过振型空间形态的物理合理性检查，才能确认识别结果不是数值伪模态。

阻尼比的结果在三次重复时间段里分别是：

• 一阶：0.8%、0.9%、0.8%
• 二阶：1.5%、1.6%、1.4%
• 三阶：1.0%、1.1%、1.0%

试想如果不做多段重复，一次测试的阻尼比随机误差可能超过20%。这类统计稳定性信息，在结构健康监测里的阈值设置和模型修正中非常重要。

10. 后续还能怎么扩展

SSI-COV不是一个只能离线处理的算法，现在已经有不少在线识别方案，用滑动窗口实时更新协方差矩阵和模态参数；也有将它与贝叶斯方法结合，输出模态参数的概率分布；还有人在做自动化OoMA（Operational Modal Analysis），把稳定图自动生成、自动筛选、自动报告封装成一套流程。对于Matlab使用者来说，从本文框架出发扩展并不难：

• 把Toeplitz矩阵构造改成稀疏存储，可以支持更大的数据量和更多通道
• 把SVD替换为随机化SVD，计算速度能提升一个量级，适合无人机巡检采集的大量数据
• 与有限元模型结合，通过MAC和频率残差做模型修正，这是结构健康监测领域很热的方向

我个人建议，如果研究方向是结构健康监测，SSI-COV这一套方法论不仅要会调用，更要亲手推导一遍状态空间模型和协方差Toeplitz分解。这个过程能把很多零散知识串起来，也不至于在数据异常时完全摸不着头脑。

最后分享一个小体会：识别结果出来后，先别急着相信计算机给的那几个数。我在实际项目中已经养成习惯，把识别出的频率和理论估算先对一眼，把振型形状画出来和结构几何形状对比一下。凡是和物理直觉冲突很大的结果，九成是输入数据或参数设置出了问题，而不是结构真的出了什么奇异的模态。用这套思路搭配SSI-COV，做起模态分析来会稳很多。