CN111241904B - 一种基于盲源分离技术的欠定情况下运行模态识别方法 - Google Patents

Abstract

本发明涉及一种基于盲源分离技术的欠定情况下运行模态识别方法，属于大型陆基平台、航空航天、船舶和建筑检测领域。本发明提高基于盲源分离的运行模态分析方法在欠定情况、近频模态等情况下的模态参数辨识精度、提高计算效率、降低计算成本。此外，降低该方法在模态参数辨识过程中对模态阶数等先验知识的依赖，降低对噪声信号的敏感性，使得该方法计算量小、鲁棒性强、方便使用。本发明即使在缺乏专业知识背景的情况下也能进行操作，能够在结构动力学工程应用中广泛应用于线性结构欠定情况下的模态参数辨识。

Description

一种基于盲源分离技术的欠定情况下运行模态识别方法

技术领域

本发明涉及一种基于盲源分离技术的欠定情况下运行模态识别方法，尤其涉及一种基于系数成分分析和基于密度聚类的运行模态参数辨识方法，属于大型陆基平台、航空航天、船舶和建筑检测领域。

背景技术

随着机械设备逐步趋于大型化、高效化和复杂化，大型陆基平台、航空航天、船舶、建筑等各个领域对机械系统的综合性能提出了更高的要求。实际机械系统在运行过程中，产生的有害振动不仅影响自身以及其他设备的安全性、稳定性以及可靠性，同时会危害操作人员的健康安全。目前，一方面，设计人员常常通过调整结构参数使其固有频率远离共振区，从而降低设备有害振动；另一方面，依赖于外界激励和结构响应信息，通过设备提供的“控制力”对设备进行反向补偿来实现对有害振动的主动控制。然而，上述手段均需要提供对象系统较为准确的固有频率、阻尼比和模态振型等特性参数。显然，相较于理论建模的高误差，难求解，实时性差等缺点来说，最直接可靠的系统参数识别途径是通过实验测试手段，直接获得系统实际的振动特性参数和结构参数。

对于某些复杂的机械设备，各组件装配过程中存在一定的预留间隙，设备运行过程中会产生部件变形，这些原因使得系统实际参数和设备运行参数相差较大，因此需要对系统进行实际工作状态下的参数识别。而复杂机械系统振源较多，难以实时获得激励数值，不能满足传统实验模态参数识别技术，且在实际运行过程中由于控制成本，体积不允许等原因只能通过有限的传感器采集设备振动信号，因此有必要研究在欠定情况下的运行模态参数识别技术。

运行模态分析与有限元计算模态分析法、依赖输入输出的传统模态分析法共同构成了现在的模态分析领域。运行模态分析技术经过几十年的发展，产生了很多的方法，这些算法大多是基于参数模型的，这些方法为传统运行模态分析方法。传统运行模态分析方法分为时域法和频域法两类。时域方法的优点在于无须对信号进行傅里叶变换，减少计算量，避免了频谱泄露，但其对噪声较为敏感。而频域方法通过经典的谱平均方法来降低噪声的影响，但同时将面临频谱泄露等问题。无论是时域法还是频域法，存在对噪声敏感度高、不能区分近频或者重频模态、计算量大等问题。

总之，传统试验模态参数辨识法中需要明确的输入激励和输出响应，在对运行中的设备进行模态参数辨识的应用存在较大困难，而现有的运行模态参数辨识法通常需要操作人员有较强的领域知识，需要明确模态阶数等先验知识，且算法理论复杂，计算资源占用较高，在模态近频、高阻尼比、和欠定条件下的辨识精度不能得到保证。

发明内容

本发明的目的是为了解决机械系统工作多年，产生变形，危害的问题，提供一种基于盲源分离技术的欠定情况下运行模态识别方法。使用该方法辨识得到的参数来指导机械系统的评估，维护，提高安全系数。

针对运行过程中的结构体模态参数辨识中存在的上述问题，本发明公开一种基于盲源分离技术的欠定情况下运行模态识别方法要解决的技术问题是：提高基于盲源分离的运行模态分析方法在欠定情况、近频模态等情况下的模态参数辨识精度、提高计算效率、降低计算成本。此外，降低该方法在模态参数辨识过程中对模态阶数等先验知识的依赖，降低对噪声信号的敏感性，使得该方法计算量小、鲁棒性强、方便使用。本发明即使在缺乏专业知识背景的情况下也能进行操作，能够在结构动力学工程应用中广泛应用于线性结构欠定情况下的模态参数辨识。

本发明公开的一种基于盲源分离技术的欠定情况下运行模态识别方法，首先使用时频分析法连续小波变换将多自由度振动系统的自由响应信号进行稀疏表示，将信号从时域变换到时频域满足盲源信号分离的稀疏限定要求；采用单源点筛选法SSP对稀疏域内的观测信号进行筛选，得到严格满足稀疏性的观测样点，这些样点的实部向量和虚部向量共线，使用能量筛选法筛选出能量占比高的点，剔除部分噪声点；对筛选后的点进行标准化处理后通过基于密度的聚类法DBSCAN提取各阶模态振型方向向量；通过压缩感知中的基追踪技术在欠定情况下进行源信号重构，得到稀疏域内的单频源信号，从所述源信号中提取各阶模态固有频率并通过半功率带宽法进行一次阻尼比估计；使用连续小波变换的逆变换将重构的源信号从时频域信号转换为时域信号，使用对数减幅值率法对阻尼比进行二次估计，最后加权融合两次估计的阻尼比得到最终估计的各阶阻尼比，完成所有模态参数的辨识。

根据所述的一种基于盲源分离技术的欠定情况下运行模态识别方法得到的设备运行过程中的固有频率、阻尼比、模态振型等结构模态参数，能够检测设备是否符合工程结构标准及隔振等要求，能够指导工程结构的设计。另外，还能为结构的健康监测、结构故障诊断、结构振动控制等方面的应用提供有力的支持，具有广泛的应用前景与效益。

本发明目的是通过下述技术方案实现的。

一种基于盲源分离技术的欠定情况下运行模态识别方法，包括如下步骤：

步骤1：对待辨识参数结构的自由振动加速度观测数据进行稀疏表示，将观测数据从时域转换到稀疏域；所述待辨识参数结构包括：机械结构的固有频率、阻尼比和模态振型。

步骤1.1：观测数据处理。将安置在待测结构体上的多个传感器测量得到的振动加速度数据整理、合并得到观测信号矩阵X^m×T，即混合信号矩阵。所述混合信号矩阵的上角标中，m表示观测信号个数，T表示每个信号的采样点个数。

步骤1.2：Morlet小波变换参数选取。由于Morlet复小波具有较优的时频局部化特性，其波形和振动数据形状相似，且在Heisenberg测不准原理问题中表现很好，因此选用Morlet小波对观测信号进行连续小波变换。Morlet小波的中心频率和带宽分别设置为fc-fb。

步骤1.3：观测数据的稀疏表示。考虑到各阶模态响应是一种在不同频率下的一组相互正交的基，可看作是不同频率下的源信号，因此模态响应在时频域内满足稀疏性。混合信号矩阵X^m×T的每一行代表一个传感器测量的时域数据，通过式(1)对X^m×T逐行进行连续小波变换。

上式中，b是时移因子，代表了小波基的平移距离，a是尺度因子，代表了小波基的缩放程度，ψ^*(·)是小波基的复共轭；t表示时间；x(t)表示t时刻结构的振动响应的幅值；ψ_a,b(t)表示t时刻小波基的幅值。得到的复数W_x(a,b)是小波变换系数，可以看作是观测信号在某一时移因子和尺度因子组合(a,b)下的能量。通过式(2)将小波变换系数的变量域从时移因子和尺度因子组合(a,b)转换到时间和频率组合(t,k)有

将所有观测数据变换完成后得到三维矩阵X^F×T×m，表示第m个观测信号在频率k、t时刻下的小波变换系数。至此，完成了响应信号的稀疏化预处理步骤。

式中，f_c为Morlet小波中心频率，T_s为信号采样周期，T为采样总时长，b_max为最大平移因子。

步骤2：对步骤1得到的稀疏域观测信号的三维矩阵X^F×T×m进行二次筛选，选出严格满足稀疏性的样点，同时剔除一部分噪声信号。

步骤2.1：对稀疏域观测信号进行单源点(SSP)筛选。对步骤1得到的三维矩阵X^{F ×T×m}取实部和虚部分别得到矩阵

和

若点是一个SSP，那么该点表示的向量与振型矩阵的列向量共线，则该点的实部和虚部向量需要对齐。满足式(3)的点即为SSP。

式中

向量表示一个在时间和频率组合(t,k)位置的散点的虚部坐标，该点的维度为m。

向量表示一个在时间和频率组合(t,k)位置的散点的实部坐标；|·|表示绝对值。令SSP阈值Δcosθ＝0.1。

步骤2.2：对稀疏信号进行能量筛选。由于时频域中还可能存在很多由噪音等生成的能量点，这些点将影响聚类精度和混合矩阵的估计精度，因此通过式(4)剔除这些低能量点。

||X(t,k,:)||＜λmax||X(t,k,:)||,0＜λ＜1 (4)

||X(t,k,:)||表示向量X(t,k,:)的长度，其中向量X(t,k,:)由时刻，频率k下的m个元素组成。将能量低于最大能量和筛选因子λ乘积的点剔除，至此完成了对稀疏空间的观测信号点的二次筛选，保留下来的点集记为X^SSP-ENG。

步骤3：对步骤2筛选得到的点集进行聚类，并估计待测结构体的估计振型矩阵W^{m ×n}。

步骤3.1：对点集X^SSP-ENG进行归一化预处理。通过式(5)对点集中每个点进行标准化处理，使得各点表示的向量模长为1，此时各点在m维空间中与原点距离相同，各点分布在以原点为圆心，半径为1的半球面上。

其中

是归一化前点集中第i个点向量的第j个元素；

是归一化后点集中第i个点向量的第j个元素。

步骤3.2：对归一化后的点集X^Nor聚类。基于密度的聚类算法(DBSCAN)有无需输入模态阶数，且该算法具有空间复杂度低、可以区分不规则形状的簇、对噪声鲁棒性高等优势，因此用该算法进行聚类。指定领域的半径大小Eps，领域中最小点的个数Min_pts，在聚类过程中余下的噪音点会被删除，有效点被分成许多簇，第h个簇包含的点集表示为X_ζ ^h＝{X₁ ^h,X₂ ^h,...,X_s ^h}。

步骤3.3：为了进一步提高估计精度，依据式(6)将点集中大于均方差的点剔除。

m是点的空间维度，即传感器个数；s是第h个簇中点的个数；

表示第h个簇中点_ζ的第i个元素；

表示第h个簇中所有点空间向量第i个元素的均值。

步骤3.4：得到估计振型矩阵。经过步骤3.3得到了较高精度的n个列向量，通过点集X^SSP-ENG记录的频率信息从小到大排序为向量K，并将振型列向量按照对应频率从小到大的顺序从左到右排列并形成估计矩阵W^m×n，矩阵W^m×n即为要求的振型矩阵。

步骤3.5：估计振型矩阵精度验证。若有该结构的理论模态振型A^m×n，则可以根据式(7)，计算估计模态振型的精度。

上式得到的结果为向量MAC＝{mac₁,mac₂,...,mac_n}，其中各元素代表该阶估计模态向量的精度。W_i表示估计振型矩阵的第i个列向量；A_i表示理论阵型矩阵的第i个列向量；mac的取值范围为0到1，越接近1表示精度越高。

步骤4：模态坐标重构，即源信号重构。通过步骤3得到的估计振型矩阵W^m×n和步骤1得到的稀疏域观测信号的三维矩阵X^F×T×m重构稀疏域源信号的三维矩阵S^F×T×n，之后通过逆小波变换将稀疏域源信号变换为时域源信号，实现模态坐标的重构。

步骤4.1：正定情况下的模态坐标重构。在正定的情况下，即m＝n，估计振型矩阵W^{m ×n}是一个方阵，此时只需要直接使用时域源信号的观测信号矩阵X^m×T通过式(8)来重建模态响应S^n×T。

S^n×T＝(W^m×n)^-1X^m×T (8)

步骤4.2：欠定情况下的模态坐标重构。对于欠定情况，即m<n，估计振型矩阵W^m×n可以从式(9)表示的线性方程组的无穷解中找出n个最具有稀疏性的解作为重构的模态响应S^n×T。

X^m×T＝W^m×nS^n×T (9)

因此在该情况下使用压缩感知领域的基追踪算法进行求解。算法如式(10)。

s(ξ)是重构后稀疏域信号；||s(ξ)||₁是信号的1范数；

表示实数空间；X是观测信号矩阵；σ是阈值；用l¹最小化方法来重构源信号，将线性系统等式约束用松弛约束替代，从而更好地应对有噪声或不完全数据的情况。使用每一个稀疏域的观测信号向量X^F×T×m(t,k,:)可以通过式(10)计算出m个长度为n的源信号向量s_t,k，将所有s_t,k存在三维矩阵中，得到重构的稀疏空间源信号矩阵S^F×T×n，其中t表示时间，k表示频率，S^F×T×n表示S是一个子空间排序为频率空间F，时间空间T和长度空间n的三维张量空间。

步骤4.3：对重构的源信号进行逆稀疏化。对步骤4.2得到的稀疏空间源信号矩阵S^F×T×n的n个子空间S^F×T(a,b)进行连续小波变换逆变换，得到模态响应S^n×T。至此完成了模态响应的重建。

步骤5：固有频率和阻尼比的估计。根据步骤4重构得到的稀疏域和时域源信号估计结构体的固有频率和阻尼比。

步骤5.1：估计固有频率。步骤4.2得到了重构的稀疏空间源信号矩阵的每一个子空间S^F×T×n(:,:,i),i＝1,2,…,n，每一个子空间上只存在一种信号频率。为了提高精度，通过式(11)将子空间行向量所有元素相加。

通过上式可以得到n个F维列向量，提取每个列向量的最大值，根据式计算出该值对应的频率，得到估计的固有频率列向量F＝{f₁,f₂,…,f_n}。

步骤5.2：通过半功率带宽法估计阻尼比。根据步骤5.1得到的固有频率F＝{f₁,f₂,...,f_n}和重构的稀疏空间源信号矩阵实部

其中i＝1,...,n,n是源信号的个数。使用半功率带宽法估计阻尼比。提取源信号矩阵实部

每个i维度的最大值Q，以

为高画一条平行于频率轴f的直线与曲线

相交于两点，这两点为半功率点，对应的频率分别为f¹，f²。则通过式(12)计算出该固有频率对应的阻尼比，列为向量ξ¹。

其中

是通过半功率带宽法得到的阻尼比向量ξ¹中的第i个元素，表示第i阶阻尼比的值，系统共有n阶，n的大小与源信号个数相同。步骤5.3：通过对数减幅法估计阻尼比。步骤4.3得到了模态响应矩阵S^n×T为单频衰减信号，可以用该信号估计每个频率对应的阻尼比。经过推导可得阻尼比和对数减幅率之间的关系如式(13)。

其中ξ表示阻尼比，δ表示对数减幅率。

单频衰减信号的对数减幅率可以通过每个周期衰减幅度来估计，为了提高精度，将q个周期的对数减幅率δ取平均，如式(14)。

上式中A_i为信号S^n×T第i个周期的峰值,A_i+1是i+1周期的峰值。由此可求得每阶固有频率对应的阻尼比，列为向量ξ²。

步骤5.4：加权融合得到阻尼比。为了提高算法对阻尼比估计的鲁棒性，通过

对步骤5.2和步骤5.3得到的阻尼比ξ¹和ξ²进行加权平均，求得最终的阻尼比。至此本文所提出的运行模态识别方法运行结束，得到了模态频率K，阻尼比ξ和振型矩阵W。

步骤6：应用步骤5辨识的结构模态参数指导结构动力学领域的结构分析与设计。

根据步骤5辨识的结构模态参数，能得到工程结构的固有频率、阻尼比和振型模态参数，检测是否符合工程结构标准及隔振等要求，能够指导工程结构的设计。另外，得到的结构模态参数还能为结构的健康监测、结构故障诊断、结构振动控制等方面的应用提供有力的支持，具有广泛的应用前景与效益。

有益效果：

1、本发明公开的一种基于盲源分离技术的欠定情况下运行模态识别方法，采用连续小波变换对时域振动信号进行稀疏转换，在满足盲源分离方法准许条件的同时保证了较高的时频分辨率；

2、本发明公开的一种基于盲源分离技术的欠定情况下运行模态识别方法，采用单源点选择法和能量选择法两次筛选，保证了信号的严格稀疏性的同时去除了部分大部分噪声，时候采用基于密度的聚类法DBSCAN提取模态振型，提高了模态振型的估计精度和算法的鲁棒性；

3、本发明公开的一种基于盲源分离技术的欠定情况下运行模态识别方法，使用基追踪算法并通过l¹最小化以及将线性系统等式约束用松弛约束替代的方法进行源信号的重构，能更好地应对有噪声或不完全数据的情况；

4、本发明公开的一种基于盲源分离技术的欠定情况下运行模态识别方法，通过用半功率带宽法和对数减幅率法进行两次阻尼比估计并将两次估计结果加权融合，提高了阻尼比的估计精度。

使用盲源分离这种无参数方法进行模态辨识相比于传统运行模态分析方法来说具有无需先验信息、计算简单等优势。且相对于比较成熟的独立分量分析方法ICA和基于二阶统计法SOBI，稀疏成分分析法SCA最大的优点在于能够求解欠定问题，除此之外使用对信号的稀疏表示能够较大程度提高计算效率。

附图说明

图1为本发明一种基于盲源分离技术的欠定情况下运行模态识别方法的流程图；

图2为具体实施方式实施例1中的五自由度弹簧-阻尼器-质量仿真系统图；

图3为具体实施方式实施例1中的五自由度仿真系统每个质量体的响应曲线；

图4为具体实施方式实施例1中使用连续小波变换对时域观测信号进行稀疏转换的结果图；

图5为具体实施方式实施例1中二次筛选各阶段过程图。(a)为二次筛选前稀疏域观测信号点在三维空间中的分布图；(b)为稀疏域的观测信号点满足系数条件的分布图；(c)为二次筛选后稀疏域观测信号点在三维空间中的分布图；(d)为二次筛选后稀疏域观测信号点关于原点堆成的分布图；

图6为具体实施方式实施例1中二次筛选前稀疏域观测信号点进行标准化在三维空间中半球面上的分布图；

图7为具体实施方式实施例1中使用基追踪算法对稀疏域源信号进行重构得到的5个源信号时频图；

图8为具体实施方式实施例1中固有频率计算图；

图9为具体实施方式实施例1中使用对数减幅率法计算阻尼比图。

具体实施方式

下面结合附图与实施例对本发明作进一步说明。

建筑领域是本发明主要应用领域之一，为了更好地说明本发明的目的和优点，下面通过对一个桥梁等效模型进行分析，该等效模型为自由振动下的五自由度弹簧阻尼结构，以此对本发明做出详细解释。

实施例1：

本实施例的五自由度弹簧-阻尼器-质量系统，如图2所示。本实施例的五自由度弹簧-阻尼器-质量系统包含质量块、阻尼器、弹簧各五个，这些弹簧阻尼器均是定常的，不随时间变化。

将动力学方程式(15)转化为状态空间方程，使用SIMULINK根据状态空间方程构建状态空间仿真模型如图2。上式中，M，C，K分别为质量、刚度、阻尼矩阵，为系统响应，如式(16)所示。F(t)在本实例中施加0.01s大小为120N的瞬态冲击。

分别对五个质量块施加Gauss白噪声激励，得到系统中五个质量块的响应。采用五个自由度位移为辨识所用的响应信号样本，将求解结果以f＝100Hz重新采样，记录时间为10s(t∈[0,10]),信号长度N＝1000，五个自由度的位移响应曲线如图3所示。

本实施例公开的一种基于盲源分离技术的欠定情况下运行模态识别方法，包括以下步骤：

步骤1：对待辨识参数结构的自由振动加速度观测数据进行稀疏表示，将其从时域转换到稀疏域。对采集到的五个通道的时域数据进行连续小波变换，将其转换到时频域，即稀疏域，得到五个通道的时频图如图4所示。图4中各子图里横坐标表示时间，纵坐标表示频率，蓝色颜色越趋近于红色表示该位置小波变换系数越大。

在本实施方式中，CWT参数设置为：选择Morlet复小波，小波中心频率和带宽参数Fc-Fb设置为3-3，尺度的数量设置为6144。

步骤2：将第一次筛选前得到的点绘制到3维空间中。从图5中(a)可以看到这些点形成了5条过原点的直线。每条直线周边有一些发现的点，这些点一部分是噪声点，另一部分是满足特殊条件，即式(4)的点。

对步骤1得到的稀疏域观测信号X^F×T×m进行二次筛选，选出严格满足稀疏性的样点，同时剔除一部分噪声信号。对稀疏信号进行单源点(SSP)筛选，公式在各个点出出计算结果如图5中(b)，蓝色为满足SSP条件的点，黄色为不满足SSP条件的点。对前步筛选得到的稀疏信号进行能量筛选。如图5中(c)所示为经过能量筛选后的散点，可以明显看出散点形成的5条直线方向更加清晰，直线周边发散的点少了很多，得到的点集为X^SSP-ENG。

在本实施方式中，SSP筛选和能量筛选的阈值分别设置为：Δcosθ＝0.1，ΔENG＝0.01。

步骤3：对步骤2筛选得到的点集进行聚类，并估计待测结构体的振型矩阵W。由图5中(c)可得散点狙击形成的直线关于原点对称，因此同一条直线拥有两种方向向量表示方法。在这里将这些散点依据某一坐标轴以原点为中心对称，得到的散点如图5中(d)。可以看到对称后的直线有唯一的方向向量。之后为了使用DBSCAN聚类法，对点集X^SSP-ENG进行归一化预处理。将这些点变换在以原点为圆心，半径为1的半球面上。变换后的点如图6所示，得到点集X^Nor。

使用基于密度的聚类方法DBSCAN对归一化后的点集X^Nor聚类，将这些分布在半球面上的点划分为一个个簇，每一簇代表一阶模态。剔除每一簇中偏差较大的点，得到了较高精度的n个振型向量。

通过点集X^SSP-ENG记录的频率信息将振型列向量按照频率从小到大的顺序从左到右排列并形成估计振型矩阵W^m×n，如式(17)。

在本实施方式中，DBSCAN聚类算法的参数设定为：领域的半径大小Eps＝0.1，领域中最小点的个数Min_pts＝150。

步骤4：模态坐标重构，即源信号重构。通过步骤3估计得到的振型矩阵W^m×n和步骤1得到的稀疏域观测信号X^F×T×m重构稀疏域源信号S^F×T×n，之后通过逆小波变换将稀疏域源信号变换为时域源信号，实现模态坐标的重构。由于本例子为正定情况，即m＝n，振型矩阵W^m×n是一个方阵，此时只需要直接使用时域源信号X^m×T通过式来重建模态响应，重建的模态响应在时频域如图7所示。

步骤5：固有频率和阻尼比的估计。对步骤4得到的的时域源信号进行连续小波变换转换到稀疏空间，对源信号矩阵的每一个子空间S^F×T×n(:,:,i),i＝1,2,...,n，将子空间行向量所有元素相加，可以得到n个F维列向量，提取每个列向量的最大值，计算出该值对应的频率如图8所示，红色点为识别的固有频率，得到估计的固有频率列向量F＝{1.5360,3.9083,5.9284,7.9746,10.8011}。

通过半功率带宽法估计阻尼比，根据式得到阻尼比向量ξ¹＝[0.02702,0.01488,0.01390,0.01550,0.01680]。之后通过对数减幅法估计阻尼比如图9所示。图中红点为每个周期振动源信号的波峰，求得对数减幅率δ，之后可以得出阻尼比，列为向量ξ²＝[0.02799,0.01504,0.01448,0.01582,0.01622]。

通过

对加权平均，求得最终的阻尼比ξ＝[0.02751,0.01496,0.01419,0.01566,0.01651]。

本实施方法案例在正定情况下识别到的结构模态参数和理论值比较如表所示。

表1在正定情况下辨识到的三种模态参数和理论值对比表

由表可知使用本发明的一种方法在只知道输出不知道输入且正定的情况下识别到的结构模态参数精度较高。

实施例2：

本实施使用与实例1相同的五自由度弹簧-阻尼器-质量仿真系统，如图2所示。系统组成及各参数设定与实例1相同。

由于该实例针对欠定情况，因此该实例选取三组样本，每组样本包含三个自由度的加速度为辨识所用的响应信号样本，分别为X_1-3-5,X_2-4-5,X_1-2-4，角标表示质量块的编号。将求解结果以f＝100Hz重新采样，记录时间为10s(t∈[0,10])，信号长度N＝1000。

本实施例公开的一种基于盲源分离技术的欠定情况下运行模态识别方法，包括以下步骤：

步骤1：对待辨识参数结构的自由振动加速度观测数据进行稀疏表示，将其从时域转换到稀疏域。对采集到的三组三个通道的时域数据进行连续小波变换，将其转换到时频域，即稀疏域。

在本实施方式中，CWT参数设置为：选择Morlet复小波，小波中心频率和带宽参数Fc-Fb设置为3-3，尺度的数量设置为6144。

步骤2：对步骤1得到的稀疏域观测信号X^F×T×m进行二次筛选，选出严格满足稀疏性的样点，同时剔除一部分噪声信号。首先对稀疏信号进行单源点(SSP)筛选。将第一次筛选后得到的点绘制到3维空间中。根据式对前步筛选得到的稀疏信号进行能量筛选。

在本实施方式中，SSP筛选和能量筛选的阈值分别设置为：Δcosθ＝0.1，ΔENG＝0.01。

步骤3：对步骤2筛选得到的点集进行聚类，并估计待测结构体的振型矩阵W。将散点依据某一坐标轴以原点为中心对称，得到的散点如图。使用DBSCAN聚类法，对点集进行归一化预处理，使用基于密度的聚类方法DBSCAN对归一化后的点集聚类。剔除每一簇中偏差较大的点，得到了较高精度的n个振型向量。通过点集X^SSP-ENG记录的频率信息将振型列向量按照频率从小到大的顺序从左到右排列并形成三组估计振型矩阵W^m×n，如表2。

在本实施方式中，DBSCAN聚类算法的参数设定为：领域的半径大小Eps＝0.1，领域中最小点的个数Min_pts＝150。

本实施方法案例在欠定情况下识别到的结构模态参数和理论值比较如表所示。

表2三种欠定情况下识别到的结构模态参数和理论值比较表

由上表可知使用本发明的一种方法在只知道输出不知道输入且欠定的情况下识别到的结构模态参数精度较高。

由表1和表2能够看出，本实施例公开的一种基于盲源分离技术的欠定情况下运行模态识别方法，能够在仅知道输出结果不知道输入激励的情况下，使用有限的传感器更加精确地辨识出线性结构固有频率、阻尼比和振型三种模态参数，辨识精度高，计算成本低，方便使用，在结构动力学领域具有良好的应用前景。

以上所述的具体描述，对发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例，用于解释本发明，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (4)

1.一种基于盲源分离技术的欠定情况下运行模态识别方法，其特征在于：包括如下步骤：

步骤1：对待辨识参数结构的自由振动加速度观测数据进行稀疏表示，将观测数据从时域转换到稀疏域；所述待辨识参数结构包括：机械结构的固有频率、阻尼比和模态振型；

步骤1.1：观测数据处理；将安置在待测结构体上的多个传感器测量得到的振动加速度数据整理、合并得到观测信号矩阵X^m×T，即混合信号矩阵；所述混合信号矩阵的上角标中，m表示观测信号个数，T表示每个信号的采样点个数；

步骤1.2：Morlet小波变换参数选取；由于Morlet复小波具有较优的时频局部化特性，其波形和振动数据形状相似，且在Heisenberg测不准原理问题中表现很好，因此选用Morlet小波对观测信号进行连续小波变换；Morlet小波的中心频率和带宽分别设置为fc-fb；

步骤1.3：观测数据的稀疏表示；考虑到各阶模态响应是一种在不同频率下的一组相互正交的基，可看作是不同频率下的源信号，因此模态响应在时频域内满足稀疏性；混合信号矩阵X^m×T的每一行代表一个传感器测量的时域数据，通过式(1)对X^m×T逐行进行连续小波变换；

上式中，b是时移因子，代表了小波基的平移距离，a是尺度因子，代表了小波基的缩放程度，ψ^*(·)是小波基的复共轭；t表示时间；x(t)表示t时刻结构的振动响应的幅值；ψ_a,b(t)表示t时刻小波基的幅值；得到的复数W_x(a,b)是小波变换系数，可以看作是观测信号在某一时移因子和尺度因子组合(a,b)下的能量；通过式(2)将小波变换系数的变量域从时移因子和尺度因子组合(a,b)转换到时间和频率组合(t,k)有

将所有观测数据变换完成后得到三维矩阵X^F×T×m，表示第m个观测信号在频率k、t时刻下的小波变换系数；至此，完成了响应信号的稀疏化预处理步骤；

式中，f_c为Morlet小波中心频率，T_s为信号采样周期，T为采样总时长，b_max为最大平移因子；

步骤2：对步骤1得到的稀疏域观测信号的三维矩阵X^F×T×m进行二次筛选，选出严格满足稀疏性的样点，同时剔除一部分噪声信号；

步骤2.1：对稀疏域观测信号进行单源点(SSP)筛选；对步骤1得到的三维矩阵X^F×T×m取实部和虚部分别得到矩阵

和

若点是一个SSP，那么该点表示的向量与振型矩阵的列向量共线，则该点的实部和虚部向量需要对齐；满足式(3)的点即为SSP；

式中

向量表示一个在时间和频率组合(t,k)位置的散点的虚部坐标，该点的维度为m；

向量表示一个在时间和频率组合(t,k)位置的散点的实部坐标；|·|表示绝对值；令SSP阈值Δcosθ＝0.1；

步骤2.2：对稀疏信号进行能量筛选；由于时频域中还可能存在很多由噪音等生成的能量点，这些点将影响聚类精度和混合矩阵的估计精度，因此通过式(4)剔除这些低能量点；

||X(t,k,:)||＜λmax||X(t,k,:)||,0＜λ＜1 (4)

||X(t,k,:)||表示向量X(t,k,:)的长度，其中向量X(t,k,:)由时刻，频率k下的m个元素组成；将能量低于最大能量和筛选因子λ乘积的点剔除，至此完成了对稀疏空间的观测信号点的二次筛选，保留下来的点集记为X^SSP-ENG；

步骤3：对步骤2筛选得到的点集进行聚类，并估计待测结构体的估计振型矩阵W^m×n；

步骤3.1：对点集X^SSP-ENG进行归一化预处理；通过式(5)对点集中每个点进行标准化处理，使得各点表示的向量模长为1，此时各点在m维空间中与原点距离相同，各点分布在以原点为圆心，半径为1的半球面上；

其中

是归一化前点集中第i个点向量的第j个元素；

是归一化后点集中第i个点向量的第j个元素；

步骤3.2：对归一化后的点集X^Nor聚类；基于密度的聚类算法(DBSCAN)有无需输入模态阶数，且该算法具有空间复杂度低、可以区分不规则形状的簇、对噪声鲁棒性高等优势，因此用该算法进行聚类；指定领域的半径大小Eps，领域中最小点的个数Min_pts，在聚类过程中余下的噪音点会被删除，有效点被分成许多簇，第h个簇包含的点集表示为X_ζ ^h＝{X₁ ^h,X₂ ^h,...,X_s ^h}；

步骤3.3：为了进一步提高估计精度，依据式(6)将点集中大于均方差的点剔除；

m是点的空间维度，即传感器个数；s是第h个簇中点的个数；

表示第h个簇中点_ζ的第i个元素；

表示第h个簇中所有点空间向量第i个元素的均值；

步骤3.4：得到估计振型矩阵；经过步骤3.3得到了较高精度的n个列向量，通过点集X^{SSP -ENG}记录的频率信息从小到大排序为向量K，并将振型列向量按照对应频率从小到大的顺序从左到右排列并形成估计矩阵W^m×n，矩阵W^m×n即为要求的振型矩阵；

步骤3.5：估计振型矩阵精度验证；若有该结构的理论模态振型A^m×n，则可以根据式(7)，计算估计模态振型的精度；

上式得到的结果为向量MAC＝{mac₁,mac₂,...,mac_n}，其中各元素代表模态的阶数估计模态向量的精度；W_i表示估计振型矩阵的第i个列向量；A_i表示理论阵型矩阵的第i个列向量；mac的取值范围为0到1，越接近1表示精度越高；

步骤4：模态坐标重构，即源信号重构；通过步骤3得到的估计振型矩阵W^m×n和步骤1得到的稀疏域观测信号的三维矩阵X^F×T×m重构稀疏域源信号的三维矩阵S^F×T×n，之后通过逆小波变换将稀疏域源信号变换为时域源信号，实现模态坐标的重构；

步骤4.1：正定情况下的模态坐标重构；在正定的情况下，即m＝n，估计振型矩阵W^m×n是一个方阵，此时只需要直接使用时域源信号的观测信号矩阵X^m×T通过式(8)来重建模态响应S^n×T；

S^n×T＝(W^m×n)^-1X^m×T (8)

步骤4.2：欠定情况下的模态坐标重构；对于欠定情况，即m<n，估计振型矩阵W^m×n可以从式(9)表示的线性方程组的无穷解中找出n个最具有稀疏性的解作为重构的模态响应S^n×T；

X^m×T＝W^m×nS^n×T (9)

因此在该情况下使用压缩感知领域的基追踪算法进行求解；算法如式(10)；

s(ξ)是重构后稀疏域信号；||s(ξ)||₁是信号的1范数；

表示实数空间；X是观测信号矩阵；σ是阈值；用l¹最小化方法来重构源信号，将线性系统等式约束用松弛约束替代，从而更好地应对有噪声或不完全数据的情况；使用每一个稀疏域的观测信号向量X^F×T×m(t,k,:)可以通过式(10)计算出m个长度为n的源信号向量s_t,k，将所有s_t,k存在三维矩阵中，得到重构的稀疏空间源信号矩阵S^F×T×n，其中t表示时间，k表示频率，S^F×T×n表示S是一个子空间排序为频率空间F，时间空间T和长度空间n的三维张量空间；

步骤4.3：对重构的源信号进行逆稀疏化；对步骤4.2得到的稀疏空间源信号矩阵S^F×T×n的n个子空间S^F×T(a,b)进行连续小波变换逆变换，得到模态响应S^n×T；至此完成了模态响应的重建；

步骤5：固有频率和阻尼比的估计；根据步骤4重构得到的稀疏域和时域源信号估计结构体的固有频率和阻尼比；

步骤5.1：估计固有频率；步骤4.2得到了重构的稀疏空间源信号矩阵的每一个子空间S^{F ×T×n}(:,:,i),i＝1,2,...,n，每一个子空间上只存在一种信号频率；为了提高精度，通过式(11)将子空间行向量所有元素相加；

通过上式可以得到n个F维列向量，提取每个列向量的最大值，根据式计算出该值对应的频率，得到估计的固有频率列向量F＝{f₁,f₂,...,f_n}；

步骤5.2：估计阻尼比；

步骤6：应用步骤5辨识的结构模态参数指导结构动力学领域的结构分析与设计；

根据步骤5辨识的结构模态参数，能得到工程结构的固有频率、阻尼比和振型模态参数，检测是否符合工程结构标准及隔振等要求，能够指导工程结构的设计；另外，得到的结构模态参数还能为结构的健康监测、结构故障诊断、结构振动控制等方面的应用提供有力的支持。

2.如权利要求1所述的方法，其特征在于：步骤5.2所述估计阻尼比的方法为：通过半功率带宽法估计阻尼比；根据步骤5.1得到的固有频率F＝{f₁,f₂,...,f_n}和重构的稀疏空间源信号矩阵实部

其中i＝1,...,n,n是源信号的个数；使用半功率带宽法估计阻尼比；提取源信号矩阵实部

每个i维度的最大值Q，以

为高画一条平行于频率轴f的直线与曲线

相交于两点，这两点为半功率点，对应的频率分别为f¹，f²；则通过式(12)计算出该固有频率对应的阻尼比，列为向量ξ¹；

其中

是通过半功率带宽法得到的阻尼比向量ξ¹中的第i个元素，表示第i阶阻尼比的值，系统共有n阶，n的大小与源信号个数相同。

3.如权利要求1所述的方法，其特征在于：步骤5.2所述估计阻尼比的方法为：通过对数减幅法估计阻尼比；步骤4.3得到了模态响应矩阵S^n×T为单频衰减信号，可以用该信号估计每个频率对应的阻尼比；经过推导可得阻尼比和对数减幅率之间的关系如式(13)；

其中ξ表示阻尼比，δ表示对数减幅率；

单频衰减信号的对数减幅率可以通过每个周期衰减幅度来估计，为了提高精度，将q个周期的对数减幅率δ取平均，如式(14)；

上式中A_i为信号S^n×T第i个周期的峰值,A_i+1是i+1周期的峰值；由此可求得每阶固有频率对应的阻尼比，列为向量ξ²。

4.如权利要求1所述的方法，其特征在于：步骤5.2所述估计阻尼比的方法为：

步骤1)，通过半功率带宽法估计阻尼比；根据步骤5.1得到的固有频率F＝{f₁,f₂,...,f_n}和重构的稀疏空间源信号矩阵实部

其中i＝1,...,n,n是源信号的个数；使用半功率带宽法估计阻尼比；提取源信号矩阵实部

每个i维度的最大值Q，以

为高画一条平行于频率轴f的直线与曲线

相交于两点，这两点为半功率点，对应的频率分别为f¹，f²；则通过式(15)计算出该固有频率对应的阻尼比，列为向量ξ¹；

其中

是通过半功率带宽法得到的阻尼比向量ξ¹中的第i个元素，表示第i阶阻尼比的值，系统共有n阶，n的大小与源信号个数相同；

步骤2)，通过对数减幅法估计阻尼比；步骤4.3得到了模态响应矩阵S^n×T为单频衰减信号，可以用该信号估计每个频率对应的阻尼比；经过推导可得阻尼比和对数减幅率之间的关系如式(16)；

其中ξ表示阻尼比，δ表示对数减幅率；

单频衰减信号的对数减幅率可以通过每个周期衰减幅度来估计，为了提高精度，将q个周期的对数减幅率δ取平均，如式(17)；

上式中A_i为信号S^n×T第i个周期的峰值,A_i+1是i+1周期的峰值；由此可求得每阶固有频率对应的阻尼比，列为向量ξ²；

步骤3)，加权融合得到阻尼比；为了提高算法对阻尼比估计的鲁棒性，通过

对步骤1)和步骤2)得到的阻尼比ξ¹和ξ²进行加权平均，求得最终的阻尼比；至此本文所提出的运行模态识别方法运行结束，得到了模态频率K，阻尼比ξ和振型矩阵W。

Images