具有模态指标的随机子空间识别【包括一致模态指标和模态参与因子】（Matlab代码实现）

👨‍🎓个人主页

💥💥💞💞欢迎来到本博客❤️❤️💥💥

🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。

⛳️座右铭：行百里者，半于九十。

💥1 概述

BSRMWR 可等效为一个多自由度系统，对于多自由度振动系统，其动力学方程为

编辑

其中［M］、［C］、［K］分别表示系统质量矩阵、阻尼矩阵 和 刚 度 矩 阵，三 者 皆 为 实 对 称 矩 阵。{ ·· x } 、 { · x} 、{ x} 分别表示加速度、速度和位移向量。{ F ( t) } 表示系统激励力。令{ F( t) } = { F} e jωt ，{ x} = { X} e jωt ，得:

编辑

由于矩阵［M］、［C］和［K］中存在非零的非对角元素，因此式( 2) 是一组具有耦合项的方程组。运用模态矩阵作为变换矩阵，可进行坐标变换。其中模态矩阵和模态坐标如式( 3) 和( 4) 所示。

编辑

编辑

本文用于识别受高斯白噪声激励影响的2DOF系统，并增加激励和响应的不确定性（也是高斯白噪声）。

指标.EMAC ： 扩展模态幅度相干

指标.MPC ： 模态相位共线性

指标.CMI ： 一致模式指示指标

partfac ： 参与因子

随机子空间识别（Random Subspace Method, RSM）是一种机器学习方法，主要用于增强 ensemble 学习中的多样性，通过在不同的子空间上构建多个基学习器（例如决策树、支持向量机等），然后综合它们的预测来提升整体的预测性能。在处理高维数据时，RSM 能够有效避免过拟合并提高模型的泛化能力。

当涉及到模态指标，特别是“一致模态指标”（Consistent Modal Indicators, CMI）和“模态参与因子”（Modal Participation Factors, MPF）时，我们通常是在讨论模态分析领域，这与结构动力学、振动分析等相关，而非直接与随机子空间识别这一机器学习方法相关联。但我们可以尝试将这两个概念放在一起，探讨如何在结构健康监测、信号处理等领域运用随机子空间识别与模态参数分析的方法。

一致模态指标 (CMI)

在结构动力学中，一致模态指标通常用于评估结构振动响应数据中特定模态的贡献程度。它提供了一种量化特定频率范围内模态活动的方法，帮助识别和分离系统的主要振动模式。CMI 可以用于判断在某个频段内，某个模态是否被清晰地激发或是否与其他模态重叠。

模态参与因子 (MPF)

模态参与因子是描述结构各自由度在某一特定模态振动下的相对重要性的量，反映了结构各个部分对于该模态振动的贡献程度。简单来说，MPF 描述了某一模态下每个节点或自由度的振动幅值相对于整体模态振型的比值。高MPF值表明对应节点在该模态振动中起着主导作用。

随机子空间识别与模态分析的结合

虽然RSM本身不是模态分析的传统工具，但在处理大量传感器数据（如来自桥梁、建筑物的振动监测数据）进行模态参数识别时，RSM可以发挥作用。具体应用方式可能包括：

1. 特征降维与模态提取：通过在不同随机选取的子空间上执行主成分分析（PCA）、奇异值分解（SVD）等，可以有效降低数据维度，并突出显示与特定模态相关的主成分，从而辅助模态参数的识别。
   
2. 多模型集成：构建多个基于不同子集特征的基学习器来识别模态参数（如固有频率、阻尼比、模态形状），最终通过投票或平均等策略综合各个模型的输出，以提高识别精度和稳健性。
   
3. 噪声抑制与特征增强：在高度嘈杂的数据环境中，随机选择的子空间可以帮助过滤掉不相关的噪声特征，强化与感兴趣模态相关的信号，从而提高模态参数估计的准确性。
   

具有模态指标的随机子空间识别方法研究综述

一、随机子空间识别（SSI）方法的基本原理

随机子空间识别（Stochastic Subspace Identification, SSI）是一种基于环境激励的时域模态参数识别方法，由Van Overschee和De Moor于1995年提出。其核心思想是通过构造状态空间模型，结合矩阵分解技术（QR分解、奇异值分解SVD）和最小二乘法，直接从结构响应数据中提取模态参数（频率、阻尼比、振型）。SSI方法适用于环境振动激励下的模态分析，无需人工输入激励信号，且能够处理高噪声数据。

数学模型：

结构系统的动力学行为可表示为离散时间随机状态空间方程：

xk+1=Axk+wkyk=Cxk+vk

其中，A为系统矩阵，C为输出矩阵，wk和vk分别为过程噪声和观测噪声。通过构造Hankel矩阵并对其进行SVD分解，可提取系统矩阵AA和输出矩阵CC，进而通过特征值分解获得模态参数。

算法流程：

1. 数据预处理：采集结构响应数据并进行标准化处理。
2. 构造Hankel矩阵：将时域数据按时间延迟排列为分块矩阵，用于描述系统的动态特性。
3. 投影与分解：通过QR分解降维，再对投影矩阵进行SVD分解，提取系统可观测性子空间。
4. 系统矩阵求解：利用最小二乘法估计AA和CC，并通过特征值分解得到模态频率和阻尼比。
5. 稳定图分析：通过不同阶次模型的参数稳定性筛选真实模态。

二、关键模态指标的定义与作用

1. 一致模态指标（Consistent Modal Indicator, CMI）

CMI用于评估模态参数在不同算法或数据条件下的稳定性，反映模态的真实性。其数学表达式为：

CMI=EMAC×MPC

其中，EMAC（Extended Modal Amplitude Coherence）衡量模态幅值的一致性，MPC（Modal Phase Collinearity）检验振型相位的共线性。CMI值超过80%时，模态被认为是可靠的；低于此阈值需进一步验证。

改进应用：在SSI算法中引入CMI约束，通过优化目标函数最小化不同估计方法间的模态参数差异，例如：

编辑

该方法显著提高了模态密集系统的识别精度。

2. 模态参与因子（Modal Participation Factor, MPF）

MPF描述各阶模态在特定激励方向上的贡献度，物理意义为模态振型与激励向量的投影关系。其计算公式为：

编辑

其中，M为质量矩阵，r为激励方向向量。MPF值越大，表明该模态对系统响应的贡献越显著。工程中通常要求累计有效质量占比超过80%以确保分析可靠性。

应用场景：

• 筛选主导模态，避免冗余计算；
• 评估结构在特定载荷（如地震、风载）下的动力响应特性。

三、SSI算法的改进与联合指标应用

1. 算法优化方向

• 抗噪性提升：采用双协方差Hankel矩阵构造方法，结合聚类分析剔除虚假模态。
• 自动化流程：基于模糊C均值聚类（FCM）和模态相位偏差（MPD）实现模态自动筛选与定阶。
• 计算效率改进：通过滑动窗口技术和并行计算优化Hankel矩阵分解速度。

2. CMI与MPF的联合应用案例

• 桥梁健康监测：在南京长江三桥模态分析中，CMI用于验证模态稳定性，MPF筛选出对车辆荷载敏感的主导模态，两者结合提高了损伤识别的准确性。
• 电力系统低频振荡：改进的双协方差SSI算法利用CMI剔除噪声模态，MPF定位关键振荡模式，实现了类噪声信号下的精确参数辨识。
• 航空航天结构：NASA通过CMI评估ERA算法识别的模态可靠性，结合MPF优化传感器布局，提升结构动态模型精度。

四、工程应用与挑战

典型案例：

• 高土石坝模态识别：基于地震记录的SSI方法结合聚类分析，自动提取Nuozhadu坝体的固有频率和阻尼比，验证了方法的抗干扰能力。
• 悬臂梁仿真分析：通过Marc软件计算MPF，证明低阶模态的有效质量占比超过70%，指导了结构动力设计。

挑战与展望：

• 计算复杂度：Hankel矩阵维度高导致内存消耗大，需进一步优化分解算法。
• 非线性系统适应性：现有SSI方法基于线性假设，如何扩展至非线性时变系统仍需研究。
• 多源数据融合：结合深度学习技术处理多传感器异构数据，提升模态参数估计的鲁棒性。

五、结论

随机子空间识别方法通过引入一致模态指标和模态参与因子，显著提升了模态参数识别的精度与可靠性。CMI为模态验证提供了量化标准，而MPF优化了模态筛选流程，两者联合应用在桥梁、大坝、电力系统中展现出重要价值。未来研究需聚焦计算效率提升、非线性扩展及多学科交叉应用，以推动SSI方法在复杂工程中的更广泛应用。

📚2 运行结果

编辑

编辑

部分代码：

function [Result]=SSID(output,fs,ncols,nrows,cut)

%Input:

%output: output data of size (No. output channels, No. of data)

%fs: Sampling frequency

%ncols: The number of columns in hankel matrix (more than 2/3 of No. of data)

%nows: The number of rows in hankel matrix (more than 20 * number of modes)

%cut: cutoff value=2*no of modes

%Outputs :

%Result  : A structure consist of the below components

%Parameters.NaFreq : Natural frequencies vector

%Parameters.DampRatio : Damping ratios vector

%Parameters.ModeShape : Mode shape matrix

%Indicators.EMAC : Extended Modal Amplitude Coherence

%Indicators.MPC  : Modal Phase Collinearity

%Indicators.CMI  : Consistent Mode Indicator

%Indicators.partfac : Participation factor

%Matrices A,C: Discrete A and C matrices

%--------------------------------------------------------------------------

[outputs,npts]=size(output);       % Computes the size of matrix output

if outputs > npts             % Check if output matrix size is proper or should be transposed

   output=output';

   [outputs,npts]=size(output);

end

%--------------------------------------------------------------------------

% Find block sizes

brows=fix(nrows/outputs);     % brows = how many output blocks.

nrows=outputs*brows;          % Redefine the row numbers.

bcols=fix(ncols/1);           % bcols = how many time steps.

ncols=1*bcols;                % Redefine the column numbers.

m=1;                          % inputs number.

q=outputs;                    % outputs number.

%--------------------------------------------------------------------------

% Form the Hankel matriices Yp,Yf.

Yp=zeros(nrows/2,ncols); %Past Output

Yf=zeros(nrows/2,ncols); %Future Output

for ii=1:brows/2

   for jj=1:bcols

       Yp([(ii-1)*q+1:ii*q] ,[(jj-1)*m+1:jj*m] )=output(:,((jj-1)+ii-1)*m+1:((jj)+ii-1)*m);

   end

end

for ii=brows/2+1:brows

   i=ii-brows/2;

   for jj=1:bcols

       Yf([(i-1)*q+1:i*q] ,[(jj-1)*m+1:jj*m] )=output(:,((jj-1)+ii-1)*m+1:((jj)+ii-1)*m);

   end

end

%--------------------------------------------------------------------------

% Projection

O=Yf*Yp'*pinv(Yp*Yp')*Yp;

%--------------------------------------------------------------------------

% Decompose the data matrix

[R1,Sigma1,S1]=svd(O,0);

sv=diag(Sigma1);

%--------------------------------------------------------------------------

% Truncate the matrices using the cutoff

D=diag(sqrt(sv(1:cut)));      % build square root of the singular values.

Dinv=inv(D);                  % (sigma)^(-1/2)

Rn=R1(:,1:cut);               % use only the principal eigenvectors

Sn=S1(:,1:cut);               % use only the principal eigenvectors

Obs=Rn*D;                     % Observability matrix

%--------------------------------------------------------------------------

% Calculate the realization of A and find eigenvalues and eigenvectors

A=pinv(Obs(1:nrows/2-q,:))*Obs(q+1:nrows/2,:);  % build A

C=Obs(1:q,:);            % build C

clear Yp Yf;

%--------------------------------------------------------------------------

% Extract the modal frequencis , damping ratios and natural frequencis

[Vectors,Values]=eig(A);       % Eigenvalues and Eigenvectors

Lambda=diag(Values);           % roots in the Z-plane

s=log(Lambda).*fs;             % Laplace roots

zeta=-real(s)./abs(s)*100;     % damping ratios

fd=(imag(s)./2./pi);            

                              % damped natural freqs:

shapes=C*Vectors;              % Mode shapes.

%--------------------------------------------------------------------------

%  Calculate  Modal Participation factors

InvV=inv(Vectors);

partfac=std((InvV*D*Sn')')';           %InvV involved in transforming states to the modal amplitudes by uncoupling A matrix

%--------------------------------------------------------------------------

% Calculate EMAC values

partoutput22=(C*A^(brows/2-1)*Vectors)';               %Predicted Last Block in Obs*Vectors Matrix

partoutput2=(Rn(nrows/2-q+1:nrows/2,:)*D*Vectors)';    %Actual Last Block in Obs*Vectors Matrix (Vectors involved in transforming states to the modal amplitudes by uncoupling A matrixx)

for i=1:1:cut

   

sum1=0;

sum2=0;

   

for k=1:1:q

           

Rik=abs(partoutput22(i,k))/abs(partoutput2(i,k));

if Rik>1

Rik=1/Rik;

end

Wik=1-(4/pi)*abs(angle(partoutput22(i,k)/partoutput2(i,k)));

           

if Wik<0

Wik=0;

end

           

EMACout=Rik*Wik;

           

EMACijk=EMACout;

sum1=sum1+EMACijk*(abs(partoutput2(i,k)))^2*100;

sum2=sum2+(abs(partoutput2(i,k)))^2;

end

       

EMACC(i)=sum1/sum2;

end

%--------------------------------------------------------------------------

% Calculate MPC values

🎉3 参考文献

部分理论来源于网络，如有侵权请联系删除。

[1]朱伟明,杨艳,刘泽远,刘程子.基于模态参与因子的无轴承开关磁阻电机振动分析与抑制[J].微电机,2022,55(10):8-15.DOI:10.15934/j.cnki.micromotors.2022.10.001.

[2]熊笑雷,赵鸣.应变模态指标在平板网架结构损伤识别分析中的应用[J].结构工程师,2012,28(02):58-65.DOI:10.15935/j.cnki.jggcs.2012.02.012.

[3] Van Overschee, Peter, and B. L. De Moor. Subspace identification for linear systems: Theory—Implementation—Applications. Springer Science & Business Media, 2012.

[4] R. Pappa, K. Elliott, and A. Schenk, “A consistent-mode indicator for the eigensystem realization algorithm,” Journal of Guidance Control and Dynamics (1993), 1993.

[5] Al Rumaithi, Ayad, "Characterization 资料获取，更多粉丝福利，MATLAB|Simulink|Python资源获取【请看主页然后私信】