企业官网建设流程全解析

做现场模态测试的工程师大概都有过这种经历：结构明明就在那，环境激励也一直在，但你就是拿不出一份让人信服的阻尼比。频域方法识别频率和振型还说得过去，一碰阻尼比就飘，今天识别出1.2%，明天变成2.8%，谁都不敢签字。后来我把重心转到时域方法，认真跑了SSI-COV（协方差驱动的随机子空间识别），才算把这套参数识别组合拳真正落地。这篇文章我就用Matlab从零实现一次，把原理、代码、参数选择、稳定图判读和现场踩坑一次性讲清楚。目标读者是正在做结构动力学实验、环境激励模态测试，或者写论文需要模态参数识别代码的工程师和研究生，只要你有基本的Matlab和振动理论底子，应该都能跟上。

1. 为什么环境激励场景下，我最终选了SSI-COV

1.1 频域方法在阻尼比识别上的硬伤

在模态测试里，如果还能用锤击或者激振器做人工激励，那问题的确简单很多——你可以直接测频响函数，然后做有理分式拟合，频率、阻尼、振型都能拿到比较稳定的结果。但现场很多工况根本没法人工激励：大跨桥梁靠风、车辆、地脉动，高层建筑靠环境振动，海上平台靠波浪。激励不可控、不可测，只有响应数据，频响函数算不出来，经典的做法是拿功率谱密度曲线凑合当FRF用。

峰值拾取法在这个框架下用得很普遍，功率谱的峰值位置对应固有频率，峰值附近半功率带宽可以估算阻尼比。问题恰恰出在这个“半功率带宽”上，环境激励下的功率谱噪声底很厚，峰值经常被毛刺顶得歪七扭八，带宽稍微一偏，阻尼比能差出一倍。遇到两个模态频率比较靠近的情况，谱峰互相叠加，半功率法直接失效。FDD（频域分解）本质上是把多自由度谱矩阵通过SVD解耦成多个单自由度来曲线拟合，频率和振型精度提升明显，但阻尼比依然要靠谱峰形状拟合，抗噪能力有限，尤其在低阻尼结构上，结果方差非常大。

说白了，频域方法适合“看模态”，不适合“抠阻尼”。

1.2 时域方法的选择：SSI家族为什么更靠谱

时域方法绕过频响函数，直接从时间序列里挖系统矩阵。最经典的ITD、STD和后来的ERA都是做自由响应或者脉冲响应，但环境激励下你拿不到理想自由响应，只能用响应信号自身做统计处理，于是才有了随机子空间识别（Stochastic Subspace Identification，SSI）的用武之地。SSI的核心思路很直接：环境激励虽然不可测，但可以当成白噪声过程建模，这样系统方程就是一个带随机输入的离散状态空间模型，模态参数全部藏在系统矩阵A的特征分解里。

SSI通常分两个分支：一个是直接处理输入输出数据矩阵的SSI-DATA，需要做QR分解和卡尔曼滤波迭代，计算量大且实现复杂；另一个是只处理输出协方差矩阵的SSI-COV，即协方差驱动的随机子空间识别。我在实践中主推SSI-COV，理由很实在：它把目标数据先压缩成协方差块Toeplitz矩阵，再做SVD，算法结构清晰、计算稳定，而且非常容易用Matlab手写实现，不需要额外工具箱。对工程批处理来说，SSI-COV是“不出错”的首选方案。

1.3 我对SSI-COV适用范围的理解

需要说明，SSI-COV的数学前提是激励为白噪声、系统为线性时不变。实际环境激励不可能是纯白噪声——风谱低频分量很重，地脉动里还有人活动、机器运转的窄带成分。但大量实测经验表明，只要把明显趋势和窄带干扰先处理掉，SSI-COV的结果仍然可用，尤其是低阶模态的识别已经是非常成熟的工程手段。

下面我先把SSI-COV的数学骨架讲透，再给可直接跑的Matlab代码。

2. SSI-COV到底在算什么东西：协方差矩阵背后的状态空间重组

2.1 从状态空间模型到输出协方差

离散时间状态空间模型是SSI-COV的起点：

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

x是系统内部状态向量，y是被测响应向量，A是离散状态矩阵，C是输出矩阵，w和v分别代表过程噪声和测量噪声，都假设为零均值白噪声。这里的关键转折是：如果只关注输出y，激励的细节其实可以不管，输出协方差序列本身就携带了系统动力学的全部信息。

定义输出协方差矩阵Ri = E[y(k+i) * y(k)^T]。把状态方程代入，利用w和v的白噪声性质，可以推得：

Ri = C * A^(i-1) * G

其中G = E[x(k+1) * y(k)^T]是下一时刻状态与当前输出的互协方差矩阵。这个式子非常漂亮——不同延时的输出协方差序列，本质上是系统矩阵A的幂次序列被C和G左右夹住。频率、阻尼比信息就藏在A的特征值里，振型藏在与C相关的特征向量里。

2.2 块Toeplitz矩阵的构建：把时间序列重组成矩阵

单看一个Ri不够，SSI-COV把从i=1到2i-1的全部延时协方差块组装成如下分块Toeplitz矩阵：

T(1:1) T(1:2) ... T(1:i) T(2:1) T(2:2) ... T(2:i) ... T(i:1) T(i:2) ... T(i:i)

这个矩阵的行被分成了i个块，每一块对应被测通道数。如果把这个组装过程推到底，T可以严格分解为可观性矩阵O和另一个矩阵Γ的乘积：T = O * Γ。这个因子分解是SSI-COV最核心的一步，它把协方差矩阵“拆开”之后，系统矩阵A就被夹在O里面了。

2.3 SVD分解与系统矩阵提取

对T做奇异值分解：

T = U * S * V^T

取前n个最大奇异值（n就是系统阶次），得到截断后的U1、S1、V1，然后构造：

O = U1 * sqrt(S1)

O的前Nch行就是输出矩阵C，O的位移结构给出了提取A的钥匙。因为可观性矩阵满足位移性质，去掉前Nch行剩下的子矩阵，等于去掉后Nch行的子矩阵左乘A：

A = pinv(O(1:end-Nch, :)) * O(Nch+1:end, :)

这一步在代码里就是两行矩阵运算。拿到A之后，直接做特征值分解，特征值λ与连续时间极点μ的关系是μ = ln(λ) / dt。极点的模对应频率（需要考虑共轭对出现两次），实部对应阻尼比：

f = abs(μ) / (2π) ζ = -real(μ) / abs(μ)

振型则通过C乘特征向量矩阵得到。整个流程下来，没有用到任何频响函数，没有人为选峰，全部由线性代数运算完成。

2.4 为什么用协方差矩阵而不是直接对数据做SVD

这是很多初学的人会问的问题：既然数据矩阵本身也能做SVD，为什么不直接上数据，非要先算协方差？原因有两层。第一，协方差矩阵是对二阶统计量的平均，相当于对数据做了自相关处理，随机噪声在平均过程中被显著压缩，直接对含噪数据做SVD会产生大量噪声主导的奇异值，系统的真实阶次很难辨认。第二，协方差矩阵的尺寸只取决于(2i*Nch)，与采样点数无关，后续SVD的计算负担可控，而数据矩阵会随采样长度线性膨胀，越算越慢。从数学上看，对数据矩阵做QR分解再取R矩阵，经过投影等价后也能得到类似T的Toeplitz结构，所以SSI-COV是更经济的路径。

3. Matlab实现步骤：从仿真数据到模态参数

3.1 仿真系统设计：三自由度质量-弹簧-阻尼模型

算法写完之后必须有一个“真值已知”的算例来验证，所以我先构造一个三自由度集中质量系统。质量矩阵取单位阵，刚度矩阵取三对角形式：

M = diag([1, 1, 1]); k0 = 1000; K = k0 * [2, -1, 0; -1, 2, -1; 0, -1, 1];

阻尼采用比例阻尼C = αM + βK，取α=0.3，β=0.0002。采样率设为50Hz，仿真时长60秒。这个参数设置下，理论模态频率约为2.24Hz、6.27Hz、9.07Hz，阻尼比约为1.2%、0.8%、0.8%。具体的理论值可以在Matlab里用eig(K, M)精确计算，这里给个大致值就能验证代码逻辑。

激励用三个自由度分别施加独立白噪声，响应通过离散状态空间方程生成。我不依赖控制系统工具箱，直接用矩阵指数做离散化：

dt = 1/fs; Ac = [zeros(3), eye(3); -M\K, -M\C]; Bc = [zeros(3); inv(M)]; Ad = expm(Ac * dt); Bd = Ac \ ((expm(Ac * dt) - eye(6)) * Bc);

这段代码在高版本Matlab里可以不用工具箱跑通。生成响应后，前500个采样点作为瞬态段直接丢弃，再给响应叠加约2%的测量噪声，模拟真实采集信号。

3.2 SSI-COV核心代码：一个函数搞定识别

下面是SSI-COV的核心函数，输入是响应矩阵y（每行一个采样点，每列一个通道）、块行数i、系统阶次n、采样步长dt，输出是频率、阻尼比和振型。

function [fn, zeta, phi] = ssi_cov(y, i, n, dt) % SSI-COV: covariance-driven stochastic subspace identification % y : Nsample x Nchannel % i : number of block rows % n : system order (number of eigenvalues to keep) % dt: sampling interval [N, Nch] = size(y); y = y - mean(y, 1); % Build block Hankel matrix with 2*i block rows nb = N - 2*i + 1; H = zeros(2*i*Nch, nb); for k = 1:nb block = y(k:k+2*i-1, :); H(:, k) = block(:); end % Past and future partitions Hp = H(1:i*Nch, :); Hf = H(i*Nch+1:end, :); % Block Toeplitz covariance matrix: future times past divided by nb T = (Hf * Hp') / nb; % SVD and truncation [U, S, ~] = svd(T, 'econ'); U1 = U(:, 1:n); S1 = S(1:n, 1:n); % Observability matrix O = U1 * sqrt(S1); % Shift structure to extract A O_top = O(1:end-Nch, :); O_bot = O(Nch+1:end, :); A = pinv(O_top) * O_bot; % Eigen decomposition [Psi, L] = eig(A); mu = log(diag(L)) / dt; % Sort by frequency ascending [~, idx] = sort(abs(mu)); mu = mu(idx); Psi = Psi(:, idx); fn = abs(mu) / (2*pi); zeta = -real(mu) ./ abs(mu); % Mode shapes: C is the first Nch rows of O C = O(1:Nch, :); phi = C * Psi; % Normalize each mode by its largest absolute value for j = 1:n phi(:, j) = phi(:, j) / max(abs(phi(:, j))); end end

有几个细节需要解释。Hankel矩阵的组装用了循环，这段代码在数据量大的时候不是最快，但胜在直观，工程上足够。T矩阵的构建用Hf * Hp'而不是循环累加，是因为Matlab矩阵乘法本身就会自动并行，远快于任何循环。A矩阵用pinv而不是inv，因为位移矩阵通常不是方阵，伪逆才是最稳健的最小二乘解。

调用方式也很简单：

[fn, zeta, phi] = ssi_cov(y, 60, 6, dt);

块行数i取60，系统阶次n取6（三阶模态对应6个复共轭极点）。跑完后，把识别的频率和阻尼比与理论值对比，结果应该相当接近。

3.3 为什么输出是重复的共轭对

如果直接跑上面的代码，会发现频率数组变成6个值：三个频率各出现两次。这是因为A的复数特征值都是共轭成对出现的，一对共轭极点对应同一阶物理模态。处理方式有两种：一种是保留全部极点在稳定图里用（因为稳定图本身会把共轭对聚合）；另一种是识别完成后只取正频率端，即取imag(mu)大于0（或者按实部符号）的那一半。更稳妥的做法是按模态置信准则做聚合，但在仿真验证里，直接取前一半即可。

3.4 振型归一化：单位化不改变模态形状

振型向量本就不是唯一的，乘以任意非零常数还是同一阶振型。代码里做了“按最大绝对值归一化”，这是输出展示的好习惯。如果有后续的模态质量、模态刚度计算需求，通常会改成按“单位模态质量”归一化，即让φ^T * M * φ = 1。两者只是数值尺度不同，不改变振型形状的相对关系。

4. 算例验证与参数敏感性：i和n到底怎么定

4.1 三自由度系统识别结果对比

仿真参数和理论值如下表（精确值用Matlab计算得到，四舍五入到三位小数）：

频率识别的相对误差基本在0.1%以内，在工程上已经非常好了。阻尼比误差在3%到5%之间，明显优于频域法在噪声工况下的表现。需要坦率地说，阻尼比识别的统计方差本来就远大于频率识别，即便算法正确，单次仿真也有一定随机波动。想要更稳定的阻尼比统计，可以做多次蒙特卡洛仿真取均值，或者做更长的时域数据。

4.2 Hankel块行数i的影响：不是越大越好

块行数i决定了协方差矩阵里能容纳的最大延时，也决定了观测量矩阵的行数。i太小，低频模态的衰减信息还没充分进入协方差矩阵，无法把低频极点“撑住”；i太大，协方差矩阵尺寸膨胀，SVD计算变慢，更重要的是长延时协方差估计因样本量不足而方差增大，反而引入虚假模态。

我做了个简单测试：固定阶次n=6，把i从10扫到120，观察第一阶频率和阻尼比的变化。i在40到80之间时，结果稳定；i小于20时，第一阶频率开始向高频偏移，阻尼比明显飘；i大于100后，出现几条“准稳定”的虚假模态，奇异值衰减变缓。工程上有个实用经验：i乘以采样周期后，覆盖最低阶模态周期的2到3倍以上。本例最低阶频率2.24Hz，周期约0.446秒，采样周期0.02秒，i取60相当于覆盖1.2秒，大约是周期的2.7倍，正好落在合理区间。

4.3 系统阶次n如何初估：SVD奇异值曲线的拐点

n是极点数，等于真实模态数乘以2。但现场不知道真实模态数，所以先看奇异值分布。把T的奇异值画出来，前一段衰减很快，后面进入平缓噪声平台，拐点位置就是系统阶次的良好初估。在仿真里，n=6之后的奇异值会断崖式下降，非常清晰。

实际工程中奇异值拐点不一定锐利，因为噪声大、输出通道多时噪声平台是缓慢下降的。我的做法是把奇异值曲线作为初选，然后跑稳定图定最终阶次，这个方法在下一节展开。

4.4 过阶数的代价：虚假模态都从哪来

如果把n设得过大，比如n=30，等于强行让SVD截断保留30个“主分量”。噪声自由度不够30个也无所谓，算法会硬挤出30个极点，其中真实模态的极点位置变化很小，而噪声极点会在复平面里散布开来，频率随机分布、阻尼比偏低或异常偏高、振型空间不平滑。这类噪声极点特征非常明显：随着n再往上加，它们的位置会继续变化，不会收敛。所以判断真假模态不能只看一次计算结果，必须依赖稳定图。

5. 稳定图：从一堆数学极点里挑出真实模态

5.1 稳定图的思路：让系统阶次“自己招供”

既然n的真值未知，办法就是让n从低到高遍历一遍，把每个n下的所有极点都画在坐标图上，横轴是频率，纵轴是系统阶次。真实模态对应的极点会从某个阶次开始稳定出现：频率基本不动，阻尼比收敛，振型形状不变。噪声极点是“一锤子买卖”，换个阶次就换位置。在图上，真实模态像一条竖直的“稳定线”，噪声极点则像散落的星点，几乎一眼能分辨。

5.2 判稳三准则：频率、阻尼、振型MAC

把每个阶次的极点与上一阶次的极点匹配，用三个准则判稳：

三准则分别记作S（频率稳定）、D（频率+阻尼稳定）、V（频率+阻尼+振型全稳定）。稳定图上通常把完全稳定的点标成最深色，只满足频率的标浅色。理想情况下，真实模态在较高阶次后达到“V”稳定性并沿垂直线堆叠。

5.3 稳定图的最小实现：一次循环加一个匹配逻辑

稳定图实现不复杂，核心就是一个循环：

maxOrder = 20; orders = 2:2:maxOrder; % 偶数阶次对应完整共轭对 f_store = cell(length(orders), 1); z_store = cell(length(orders), 1); phi_store = cell(length(orders), 1); for idx = 1:length(orders) n = orders(idx); [fn, zeta, phi] = ssi_cov(y, i, n, dt); % 只保留正频率侧的极点（取虚部大于0） pos = imag(mu) > 0; % mu是从ssi_cov内部拿到的，实际使用时需要让函数返回mu f_store{idx} = fn(pos); z_store{idx} = zeta(pos); phi_store{idx} = phi(:, pos); end

实际使用时我一般会让ssi_cov额外返回mu和A，方便后续做匹配。匹配策略是：对第k阶的每个极点，找第k-1阶中频率最近的一个极点，然后检查频率、阻尼、MAC三个指标是否都落在阈值内。如果第k-1阶没有可匹配的极点，这个极点在新阶次首次出现，不判稳。这个“由近到远、逐级继承”的过程，本质上是在追踪每条极点在阶次空间的轨迹。

5.4 稳定图判读经验：宁可少选，不要多选

我处理过不少现场数据，稳定图选模态最忌讳“贪多”。真实结构在工作频段内的模态是有限的，稳定图上最漂亮的那几根竖线你直接点选就完事。碰到那些部分稳定、断断续续的点，通常是噪声模态、谐波分量或非线性响应，不要勉强选。选中阶次后，如果振型曲线明显突变、阻尼比超过5%或小于0.1%，请一定回头重新审视这个模态。

6. 实测数据预处理与常见翻车现场

6.1 去均值、去趋势：这是协方差理论的地基

SSI-COV的数学推导建立在零均值平稳随机过程假设上。电涡流位移传感器传出来的信号经常带直流偏置，加速度计有温度漂移，压电传感器还可能有长周期趋势项。这些直流和趋势分量反映在协方差矩阵里就是延时很大的缓变成分，会在极低频制造假模态。预处理必须做：先y = detrend(y, 'constant')去均值，再看趋势项是否明显，明显的话用detrend(y, 'linear')或更平滑的基线去消除。我见过最离谱的一次，把一段带0.1Hz趋势的响应直接喂给算法，识别出一个“完美”的0.05Hz模态，阻尼比0.3%，后来确认纯粹是基线漂移。

6.2 带通滤波必须用零相位，否则阻尼比就毁了

滤波器的设计直接踩过坑。普通IIR滤波（比如butter配filter）有相位延迟，等价于在信号上叠加了一个与频率相关的延时，频率识别影响不大，因为相位偏移不改变峰值位置，但阻尼比是在时域复指数衰减里提取的，相位畸变会直接改变衰减包络，识别出来的阻尼比可能偏大或偏小，没有规律。正确做法是用filtfilt做零相位滤波，在Matlab里就是同一个设计先正向再反向滤波，相位畸变基本消除。滤波带宽也不宜太窄，留出目标模态频率的上下各20%以上的余量，否则滤波器自身的振铃会制造边界处的假稳定点。

6.3 数据时长和测点数量的工程底线

协方差矩阵本质上是对信号做系综平均，样本量直接决定统计精度。我的经验底线是：识别最低阶模态至少要有200个该模态周期的数据长度。比如最低阶1Hz，采样时长至少要200秒；如果要稳定识别阻尼比，300秒以上更稳。测点数量方面，SSI-COV识别的“振型”是测点自由度上的离散振型，不是连续体振型，测点越密空间分辨率越高。如果只有三两个测点，算法依然能算出频率和阻尼，但振型只有少数几个点能描述，评估振型形状时应谨慎。

6.4 环境激励不是白噪声，SSI-COV还扛得住吗

这是现场应用里最常被质疑的问题。严格来说，SSI-COV的白噪声假设现实世界中不成立，但工程上它依然非常抗造，原因在于识别主要依赖于输出协方差序列，宽带的非白噪声激励相当于在输入谱上乘了一个平滑的谱形状，只要激励谱在目标模态附近没有剧烈尖峰或深谷，极点的位置误差通常可以接受。真正危险的是窄带激励，比如旋转机械的转频谐波会在固定频率产生极强的稳定极点，阻尼比趋近于0，判断真实模态时要以工况和先验知识交叉验证。如果风引起的涡振或桥面车流激励有明显窄带分量，建议先对数据做窄带陷波，再进入SSI-COV。

6.5 随机减量（RD）到底要不要做

随机减量技术把平稳随机响应平均成自由衰减响应，早些年常作为SSI的前处理，原因是早期计算能力弱，先降噪再识别能提高信噪比。在我看来，SSI-COV的协方差组合过程本身就已经包含了类随机减量的平均效果，如果直接SSI-COV结果已经足够好，没必要再做RD。只有在信噪比极低、协方差矩阵里噪声成分压过结构成分时，可以先做窄带滤波或RD预处理，然后再识别。多一道处理就多一个参数、多一分人为偏差，能用原数据得到稳定结果就不要画蛇添足。

6.6 同一个坑：时间单位随手写错

最后说一个最不起眼但最容易绊住人的问题。现场采集系统采样率常常是512Hz、1000Hz这类带小数的值，有人写代码时dt习惯性写1，结果识别出来的频率整体差了512倍——几百Hz的模态变成了零点几Hz，所有人都懵了。我的习惯是，在SSI-COV结果出来之后，第一时间把功率谱密度峰值频率和识别频率叠在同一张图上做交叉验证。两个来源如果对得上，说明单位、阶次、截断全都没问题；对不上，就先查dt，再查滤波参数。这个习惯救了我很多次。

写在最后的个人体会

SSI-COV这套方法我在多个项目的实测数据上跑过，从简支梁实验台到大型桥梁环境振动，整体体会是：频率识别可以做到非常自信，振型识别的空间形态取决于测点布置，阻尼比识别则永远是“参考值”——它随数据长度、信噪比、滤波方式、稳定图判据都存在统计波动。但相比频域半功率法，SSI-COV给出的阻尼比可重复性已经强太多。如果你刚接触这个方法，建议先不要直接上现场数据，而是从三自由度仿真算例入手，把i、n、稳定图三个环节都摸透，再去挑战实测信号。代码跑通只是第一步，能解释每一行矩阵运算在做什么、为什么这么做，才算真正掌握它。