基于奇异熵定阶的ERA方法优化

特征系统实现法（ERA），首先由美国国家航空与宇航局（NASA）所属的Langley研究中心于1984年提出的。它是一种属于多输入多输出（MIMO）的时域整体模态参数辨识方法。它移植了自动控制理论中的最小实现理论，利用脉冲响应数据，采用奇异值分解的方法，求得系统的特征值与特征向量，从而求得模态参数。该方法于1984年提出后，当年即在美国伽利略航天器的模态分析中应用，次年又在航天飞机载巨型太阳能帆板的太空模态试验中应用，均取得良好效果。与同时用的其他几个方法比较，该方法具有最佳精度，是目前比较先进的一种时域参数识别方法。

特征系统实现算法是利用实测的脉冲响应函数，运用奇异值分解方法，确定系统的阶次和状态方程中系统矩阵A，输入矩阵B和输出矩阵C，进而求解系统矩阵A的特征值问题，求得极点与留数，从而确定系统的模态参数。当矩阵A、B、C的阶次为最小时，即为最小实现，此时系统是可控的、可观的。

5.1.2.1　系统状态方程的描述

对于一个N自由度系统，若在P个点激励，在L个点上测量响应，可用下列状态方程描述：

式中：K为采样点序号；X（K）为在KΔ时刻系统的状态向量（2N×1），Δ为采样时间间隔；Y（K）为在KΔ时刻的实测响应向量（L×1）；F（K）为在KΔ时刻系统的输入向量（P×1）；A为系统矩阵（2N×2N）；B为输入矩阵（2N×P），又称控制矩阵；C为输出矩阵（L×2N），又称观测矩阵。

对一线性定常系统，自由响应可用脉冲响应函数来代替。因此自由响应的最小实现问题常用脉冲响应的最小实现问题来代替。

5.1.2.2　系统脉冲响应函数（或自由振动响应）的提取

对于传统的实验模态分析，系统的脉冲响应可由实测传递函数的拉氏变换求得，但大型水工结构在工作状态下的激励荷载很难测得，无法得到传递函数。考虑到诱发水工结构振动的激励荷载主要是水流脉动荷载，一般情况下都属于各态历经的平稳随机过程，且具有一定的激励频带宽度，可以认为是一种有色白噪声荷载。而结构在白噪声荷载激励作用下的脉冲响应函数或自由振动响应可由自然激励技术（NEx T法）或随机减量法提取。

1.随机减量法

随机减量法是从线性振动系统的一个或多个平稳随机响应样本中，消除或减少随机成分而获得一定初始激励下的自由响应信号。该方法的主要思想是利用平稳随机振动信号的平均值为零的性质，将包含有确定性振动信号和随机信号两种成分的实测振动响应信号进行辨别，将确定性信号从随机信号中分离出来，得到自由衰减振动响应信号。图5.1.1所示为单测点随机信号的样本平均，取起始采样幅值xs。将随机响应信号分成K个长度相等、可重叠的样本。每个样本段的起始采样幅值均取为xs=x（tk），其中k=1，2，…，K。各段样本起始点的斜率x·（tk）正负交替出现，如图5.1.1中所示。当k=1，3，5，…，奇数时，x·（tk）＞0；当k=2，4，6，…，偶数时，x·（tk）＜0。tk为第k个样本的起始采样时刻，而往后的时间以τ表示。

将K个样本进行平均，得

图5.1.1　单测点随机信号的样本平均

对于不同的τ，δ（τ）有不同的值。δ（τ）称为随机减量特征信号。实际上它是由初始位移激励而引起的自由响应，可解释为：①初始条件（初始激励）引起的自由响应；②激励而引起的自由响应（受迫暂态响应中的一部分）；③激励而引起的受迫响应。当激励为随机激励时，后两部分响应亦为随机的。当进行多次样本平均后，后两部分响应趋于零。对于第一部分响应来说，初始激励有两种：一种为初始位移xs激励；另一种为初始速度x·（tk）激励。由于在选取各段样本的起始条件时，取x·（tk）正负交替，因此在多次平均后，亦趋向于零。因此在响应中剩下的就只是由初始位移激励而引起的自由响应，亦称为阶跃自由响应。

当样本起始点x（tk）及x·（tk）取值不同时，可得到不同的自由响应曲线。即：

（1）当x（tk）=xs，而x·（tk）正负交替时，k=1，2，…，K，经多次样本平均后，可得由初始位移而引起的阶跃自由衰减响应。

（2）当取x（tk）=0，x·（tk）＞0，k=1，2，…，K时，经多次样本平均后，得到由正初速度引起的正脉冲自由衰减响应。

（3）当取x（tk）=0，x·（tk）＜0，k=1，2，…，K时，经多次样本平均后，得到负脉冲自由衰减响应。

对于多自由度系统多个测点的情况，采用下列方法提取自由振动响应具有更高的精度。

以两个自由度系统，两个测点随机响应信号为例，图5.1.2为两个测点的随机响应信号x1（t）、x2（t），这两个信号是由某些已知的或未知的随机输入引起的，并在同一时刻记录下来。

它们的随机减量信号（自由衰减响应信号）为

图5.1.2　两个测点随机响应信号样本

现对其中一个信号的起始值作如下规定：

当t=tk时

而对另一测点的信号x2（t）不作任何规定。这样，经过样本多次平均后，对第一测点来说，得到的是对初始位移x1（tk）=xs的自由响应。而对第二测点来说，得到的是对第一测点上初始位移xs在第二测点上引起的自由响应。因为对第二测点的初始条件未加任何限制，它们亦是随机的，因此在多次平均后亦消失了。

对于具有N个自由度的结构，若取N个测点，仅对其中第j个测点的初始采样条件加以限制，而对其他各测点则无任何限制。在进行多次样本平均后，得到的各点的随机减量特征信号δi（τ），i=1，2，…，N，就表示各点对j点初始位移xj（tk）=xs的自由衰减响应。所谓初始位移响应即相当于在第j点作用一个静荷载，使结构在j点产生一个静位移xs，然后将静荷载突然去掉，此时结构产生自由振动响应。此种方法要求同一结构上的各个测点必须同时采样。由于此法以多个测点中的某个测点为参考，与单独将各个测点进行随机减量相比，考虑了各个测点自由振动衰减信号之间的相关联系，提取的自由振动衰减信号更加贴切实际，有利于提高模态参数整体识别法的精度。

以某悬臂梁在水荷载作用下的结构3号点和4号点应变测量测点同时采样信号为例，如图5.1.3所示。

图5.1.3　某悬臂梁在水荷载作用下的结构3号点和4号点应变测量采样信号

以3号点为例，对比将3号点取初始幅值后进行随机减量得到的自由衰减信号与仅将4号取初始幅值，而其他点未做任何限制，进行随机减量后得到的自由衰减信号，如图5.1.4所示。

图5.1.4　利用多测点获得的衰减信号

显然，对3号点进行限制初始幅值后，得到的自由衰减信号质量不如未加任何限制而得到的自由衰减信号。可见，对于同一结构上的不同点，同时采样的条件下，仅对其中一个点取初始幅值，而其他点未加任何限制条件下得到的各个点的自由衰减信号质量会大大提高，因其考虑了测点与测点结构之间的相关关系。

进行随机减量特征提取时，一个关键性的问题是信号截取初始幅值xs的选取，信号长度一般是一定的，当截取的初始幅值取得较大时，截取的子信号段数将减少，会使有效的平均次数减少，平均效果变差；相反，若截取的初始幅值过小，虽然平均次数增多了，但由于小幅值产生的位移激振量值小，效果亦较差。实际应用时，既要保证适当大的截取初始幅值，又要保证一定数量的平均次数，显然，它们之间是相互矛盾的，这就给应用带来一定的困难。

2.NEx T法

NEx T法由美国SADIA国家实验室于1995年提出，基本原理是认为白噪声环境激励下结构两点之间响应的互相关函数和脉冲响应函数有相似的表达式，对于自由度为N的线性系统，当系统的k点受到力fk（t）的激励，系统i点的响应xik（t）可表示为

式中：φiξ为第i测点的第ξ阶模态振型；akξ为仅同激励点k和模态阶次ξ有关的常数项。

当系统的k点受到单位脉冲力激励时，就得到系统i点的脉冲响应hik（t），可表示为

对系统k点输入力fk（t）进行激励，系统i点和j点测试得到的响应分别为xik（t）和xjk（t），这两个响应的互相关函数的表达式为

假定激励f（t）是理想白噪声，根据相关函数的定义，则有

式中：δ（p-q）为脉冲函数；ak为仅同激励点k有关的常数项。将式（5.1.30）代入式（5.1.29）并积分，得到

对式（5.1.31）的积分部分进行计算并化简，得

将式（5.1.32）代入式（5.1.31）得

对式（5.1.33）做进一步的化简，整理后可得

式中，，仅为同参考点j和阶次ξ有关的常数项。

比较式（5.1.34）和式（5.1.28），可以看出，线性系统在白噪声激励下两点的响应的互相关函数和脉冲激励下的脉冲响应函数的数学表达式是完全一致的，互相关函数确实可以表征为一系列复指数叠加的形式。在这方面，互相关函数具有和系统的脉冲响应函数同样的性质。同时结构各测点的同阶模态振型乘以同一因子时，并不改变模态振型的特征，因此互相关函数可以用来代替脉冲响应函数。

5.1.2.3　脉冲响应矩阵的建立

系统的脉冲响应可由实测传递函数的拉氏逆变换求得，对各点的脉冲响应函数h（t）进行离散采样后，便可得离散的脉冲响应函数序列h（K）（K=1，2，…）。设采样的时间间隔为Δ，在KΔ时刻，各测量点的脉冲响应可构成下列脉冲响应矩阵：

式中：hij为j点激励、i点测量的脉冲响应函数；L、P分别为测量点与激励点的数目。

脉冲响应的最小实现问题是已知及求A、B、C，并使三重矩阵[A、B、C]的阶次最小。在求得系统矩阵A后，再由其特征值与特征向量确定模态参数。

5.1.2.4　构造Hankel矩阵

脉冲响应的最小实现一般是从生成Hankel分块矩阵开始，Hankel矩阵具有如下形式：

式中：ji（i=1，2，…，α-1）、tj（j=1，2，…，β-1）均为任意正整数。

从理论上讲，矩阵有不变的秩，且此秩正好是系统的阶次，系统的极点数可由它确定。但由于噪声干扰，由实测数据生产的必然有秩的亏损。当式（5.1.36）中的α、β增加到一定程度后，其秩才趋于不变。tj、ji、α、β选择的准则是能获得这个不变的秩，且矩阵的尺寸最小。常见的一种选择是令tj=ji=iΔ（i=1，2，…）。此时矩阵的逆对角元素均相等，即可得

5.1.2.5　脉冲响应矩阵与系统矩阵、输入矩阵和输出矩阵之间的关系

对线性定常系统脉冲响应与矩阵A、B、C之间有如下关系：

式（5.1.38）可推导如下：

系统的传递函数可表示为

依据式（5.1.39），并考虑脉冲响应函数与传递函数之间的Z变换关系：

可得传递函数：

式中：hlp（i）为Markov系数，即脉冲响应函数。对传递函数矩阵而言，则有

式中：为脉冲响应矩阵。

此外，系统的状态方程可写为

对式（5.1.44）进行Z变换：

将式（5.1.45）代入式（5.1.46）可得

可得传递函数

令Z-1A=X，则

将式（5.1.49）展开成泰勒级数，可得

则式（5.1.49）可写成

将式（5.1.51）代入式（5.1.48），得

将式（5.1.52）与式（5.1.43）逐一对照，可得

对式（5.1.38）进行递推，可得

继续递推，并代入式（5.1.37），可得

式中：P为可观性矩阵；Q为可控性矩阵；α，β为可观、可控性指数，且有(www.xing528.com)

由此即可导出特征系统实现算法的主要计算公式。

5.1.2.6　算法实现

由式（5.1.56），当K=1时，有

显然亦有

对）进行奇异值分解，有

式中：U为Lα×2N维左奇异向量；V为Pβ×2N维右奇异值向量矩阵；Σ=diag（λ1，λ2，…，λ2N）为2N×2N维奇异值对角矩阵，λi为奇异值，并有矩阵的秩即为系统的阶次，可由不为零的奇异值的个数来确定。

U、V是正交归一化矩阵，且满足

假设存在一个矩阵H#，满足下式：

则由式（5.1.58）可得

可见，H#是的一种广义逆矩阵。

按照矩阵广义逆的定义，H#可写成如下两种形式：

（1）若矩阵的秩等于其列数，则

（2）若矩阵的秩等于其行数，则

按上述两种逆矩阵表达式之一可推导出H#的计算公式，现以第一种形式进行推导：

由于U及V为正交矩阵，故有U-T=U，V-T=V，因此

下面确定矩阵A、B、C与U、V、Σ之间的关系，令

类似式（5.1.37）可得

对式（5.1.71）左乘右乘EP，可得

将式（5.1.59）代入式（5.1.72），得

在式（5.1.73）矩阵A的两边各插入一个单位矩阵QH#P，可得

考虑到式（5.1.67）及式（5.1.60），式（5.1.74）可改写为

式（5.1.75）与式（5.1.38）具有极其相似的表达式，即

因此系统的状态方程规范型可写为

由式（5.1.78）可见，A矩阵的阶次取决于Σ的阶次，而矩阵Σ∈R2N×2N。尽管矩阵的阶次很高（Lα×Pβ），经过奇异值分解后Σ属于2N阶。由于

由此可知，系统矩阵A为2N阶方阵；相应的状态矢量X（K）的阶次必为2N阶，它是描述2N阶系统的最小阶次，因此是最小实现。

对作奇异值分解的含义，可理解如下：

（1）从逼近理论来看，UΣVT是在所有子空间中的最佳逼近。

（2）从信号处理角度来看，用UΣVT代替相当于对数据进行了一次维纳滤波。被滤掉的是对应于奇异值为零的那些与输入、输出无关的随机噪声。因此状态方程无需再为噪声提供出口，无需再进行扩阶。

（3）以最少的参数、最小的阶次来描述系统的特征和进行运算，从而减小了运算量和内存占用量。

（4）提高了算法的抗噪声干扰能力，避免了在模态转化过程中产生计算误差，即出现虚假模态。

5.1.2.7　模态参数辨识

系统的模态参数可由系统矩阵A的特征值及特征向量来确定。系统矩阵A可由式（5.1.78）确定：

对矩阵A进行特征值分解，求出特征值矩阵，然后求出特征向量矩阵：

式中：Z为特征值矩阵；Ψ为特征向量矩阵。

由此便可确定系统的模态振型。

矩阵A的特征值与系统特征值有如下关系：

式中：Zr和λr均为复数，可写成

然后由下列公式可求得系统的模态频率ωr及模态阻尼ξr：

模态矩阵（或振型矩阵）：

模态贡献因子矩阵（或振幅矩阵）：

5.1.2.8　模态置信因子

当测量噪声不足够低且观测数据的确定性又不高时，则截断奇异值确定的秩大于实际系统能控能观的模态数，表明噪声仍然存在于缩减分块矩阵中。为此需引入模态置信度因子。

1.模态幅值相关因子γ

模态空间的实现系统包含有结构模态分量和噪声模态分量。在模态空间，理想线性系统的特征值和特征向量具有确定性，而噪声模态具有随机性，两模态空间一般不能完全分离。某阶实测模态的幅值与该阶真实模态的幅值在时间历程上的相关程度，反映了两者之间的逼近程度。真实模态可用实测的系统特征值外推构造。

令控制输入矩阵为

式中：“*”表示复共轭转置，P维向量br与系统特征值λr（r=1，2，…，2N）相对应；以做幅值系数，按指数外推得真实模态幅值时间历程：

定义复向量qr为实测模态的幅值时间历程，可表示为

故第r阶模态的幅值相关因子为

若第r阶实测模态受噪声干扰较小，由真实模态的确定性可推知，随时间的推移，模态幅值不应有明显变异，即γr趋近于1，表明该阶模态接近真实模态；反之该阶模态为非真实模态。

2.模态相位共线度因子μ

对于小阻尼结构，一般的模态形态是可以观察的。理想线性的模态振型相位角是0°或180°，即振型实部和虚部的回归直线是水平坐标轴线。由于结构非线性，模态间的耦联和信噪比低等因素，导致振型向量的相位角偏离0°或180°。据此，可用模态相位共线度因子μ描述振型精度。

设cr（r=1，2，…，2N）为相应于第j阶辨识模态的振型列向量。计算下面的值：

式中：sign（·）表示取符号。

由此定义第r阶模态相位共线度因子：

μr越趋近于1，表明第r阶模态振型精度越高；反之，μr越趋近于0，表明振型精度很低，其模态属性难以确定。

上述两个模态置信因子，γ反映模态纯度大小，μ反映模态振型精度大小。因此在筛选结构真实模态时，以γ为主，以μ为辅。

5.1.2.9　基于奇异熵的系统阶次确定

在对原始信号经奇异值分解降噪处理后，可利用NEx T法求得结构某点的系统脉冲响应，再利用ERA算法识别结构模态参数。利用ERA算法识别时首先需要构造Hankel矩阵，求出系统矩阵A，对A进行特征值分解可求得系统特征值。理论上系统的特征值个数为2N，虽然原始信号已进行降噪处理，但仍有部分噪声存在。对式（5.1.60）中的进行奇异值分解时，奇异值对角矩阵元素往往大于2N。

单凭奇异值大小确定系统的阶次往往比较困难。通过比较式（5.1.60）中）和式（4.3.8）中的吸引子轨道矩阵D可以发现，的构造和D的构造形式完全相同，反映的信息是结构脉冲响应信息矩阵，D反映的是结构原始信号信息矩阵。因此，的定阶可用基于矩阵奇异熵的方法来确定。通过式（4.3.6）、式（4.3.7）计算经奇异值分解后的奇异谱和奇异熵增量。对于同一脉冲响应信号而言，无论信号受噪声干扰的严重程度如何，对其有效特征信息进行完整抽取所需的奇异谱阶次是一定的（也即系统的阶次是一定的）。选取奇异熵增量开始降低到渐近值时的阶次作为系统结构定阶是非常合理的，这一系统定阶过程是基于奇异熵消噪与定阶的两个过程同时进行，即对脉冲响应信号再一次进行了消噪，同时也确定了结构系统的阶次。

5.1.2.10　基于奇异熵定阶降噪的ERA法的基本步骤及算法流程图

在获取泄流激励下水工结构动力响应后，首先对信号进行降噪，再利用NEx T法计算泄流结构系统脉冲响应函数参数矩阵，通过构造Hankel矩阵及利用奇异熵定阶技术，利用特征系统实现法（ERA）寻找系统的一个最小实现，得出最小阶次的系统矩阵，对系统矩阵进行特征值分解并剔除虚假模态，即得结构最终的模态参数，其基本步骤如下：

（1）获得结构测点的动力响应数据X，并构造矩阵D，利用奇异熵或者小波降噪技术进行噪声剔除，得到重构信号

（2）利用NEx T法计算结构测点的脉冲响应函数。

（3）利用脉冲响应函数构造Hankel矩阵并对进行奇异值分解。

（4）计算矩阵奇异值分解后的奇异熵，并确定奇异谱阶次，也即结构系统的阶次。

（5）最后根据已确定的阶次确定系统矩阵A，输入矩阵B和输出矩阵C。

（6）求解系统矩阵A的特征值问题，求得极点与留数，从而确定系统的模态参数。

基于奇异熵定阶降噪的ERA法识别流程图如图5.1.5所示。

图5.1.5　基于奇异熵定阶降噪的ERA法识别流程图