应变功率谱密度传递比与工作应变模态参数识别

王世东，任伟新

（1.合肥工业大学土木与水利工程学院，安徽合肥 230009；2.深圳大学土木与交通工程学院，广东深圳 518060）

引 言

模态参数是结构固有的动力特性，通过识别模态参数可以对结构的工作状态进行评估。传统的结构模态参数识别基于结构的输入（激励）和输出（响应），对于处于运营（工作）状态下的桥梁、建筑、核反应堆、大坝、海上平台等大型土木工程结构，要准确测得结构所受到的激励非常困难，而且做不到实时监测。工作模态参数分析（Operational Modal Anal⁃ysis，OMA）是仅基于响应数据，即只需测试结构响应信息，完成模态参数识别的过程，此时一般将结构工作期间所受的激励（环境激励）假定为白噪声。基于环境激励的结构工作模态参数识别方法快速发展，并形成了一系列较为经典的方法，如峰值拾取法（Peak⁃picking，PP）、频域分解法（Frequency Domain Decomposition，FDD）、多参考最小二乘复频域法、NExT 法（Natural Excitation Technique）、特征系统实现算法（Eigensystem Realization Algorithm，ERA）、随机子空间识别（Stochastic Subspace Identi⁃fication，SSI）等。大量实践表明：结构正常工作期间所受激励难以避免地包含有非白噪声如谐波成分等，OMA 方法对激励所做的白噪声假定与实际激励有差别，会影响模态参数的识别结果，甚至导致错误的识别结果［1⁃2］。

Yan 等［3］提出了功率谱密度传递比（Power Spectrum Density Transmissibility，PSDT）的概念。PSDT 定义为两响应测点i，j与一参考测点k的互功率谱密度函数Sik和Sjk的比值。从理论上证明了在系统的极点处，PSDT 收敛于两测点i和j的振型系数之比φi/φj，这一特性与激励类型和参考点位置的选择无关。因此，在一种激励工况下，对相同的测点选择不同的参考点来计算多组PSDT，在系统极点处，各组PSDT 均相等且收敛于对应测点的振型系数。随后，张昱等［4］和Li 等［5］细致地论述了以上特性。基于PSDT 的概念，Yan 等［6］进一步提出了改进的PSDT 方法，即PSDT⁃driven PP 和EPSDT 方法。Araujo 等［7⁃8］基于PSDT 相继提出了PSDTM⁃SVD方法和改进的PSDTM⁃SVD 方法。结果表明，基于PSDT 的工作模态参数识别方法具有理论价值和实际应用优势［9］。

经典的模态参数识别指位移模态参数识别，所得到的振型为位移振型［10］。和位移相比，应变对结构的局部损伤（如裂缝、孔洞等）更加敏感，能更好地反映出结构局部特性的变化，在损伤识别领域应用较多［11⁃12］。早期的模态试验通过测试位移模态然后借助中心差分法间接得到应变模态［13⁃14］，为了避免中心差分所产生的数值误差，学者们研究直接利用应变测试数据来建立应变响应模型。伊立言［15］将应变计用于模态试验，提出了应变模态的概念。Yam等［16］推导了应变频响函数（Strain Frequency Re⁃sponse Functions，SFRFs），给出了应变模态参数（频率、阻尼、应变振型）识别的实验测试方法。随后，许多经典的位移模态分析方法也被引入到应变模态参数识别中，如随机子空间方法、特征系统实现算法、频域空间域分解法等［17⁃20］。

为了实现仅基于应变响应测试的结构应变工作模态参数识别，本文定义了响应应变功率谱密度传递比，从理论上证明了其在系统的极点处为应变振型系数之比。利用这一性质，选取一系列不同的参考点构造响应应变功率谱密度传递比矩阵，在系统的极点处对该矩阵进行奇异值分解，分解所得左奇异矩阵的第一列向量即为应变振型，从而实现结构工作应变模态参数的识别。该方法不需要对激励做白噪声假设，也不需要改变激励类型，直接通过测试一种工况下的应变响应数据就能识别出结构的应变工作模态参数。

1 应变模态分析理论

由模态分析理论，结构位移响应u由位移振型向量φr和模态坐标qr叠加而成。类似地，结构的应变响应ε也可以用应变振型φεr和模态坐标q′r叠加而成，如下式所示：

位移和应变是结构同一能量状态的两种表达，对于第r阶模态，应变模态坐标q′r和位移模态坐标qr相等。对于线性时不变系统，模态坐标可以表示为：

式中kr，mr和cr分别代表第r阶模态刚度、模态质量和模态阻尼；N表示模态的数目；F表示激励力的幅值；ω表示圆频率；j 表示虚数符号，j2=−1；t表示时间。由方程（2）和（3）得到如下式所示的应变模态表达式：

由方程（3）和（4）得到j点激励，i点响应的应变频响函数（Strain Frequency Response Function，SFRF）表达式如下：

式（5）表示第r阶应变模态振型向量在i测点的系数，φjr由表示第r阶位移模态振型向量在i测点的系数，应变频响函数(ω)组成应变频响矩阵Hε如下：

将方程（6）写成如下式所示的极点留数形式：

式中s表示一个复数，表示系统的极点；表示留数矩阵；为第r阶模态参与系数向量，其含义是第r阶模态在振动响应中所占的比重；和分别为和λr的复共轭，ζr表示第r阶模态阻尼比。当不考虑阻尼的影响时有ζr=0，λr=jωr为系统的极点。此时结构的第r阶频率如下式所示：

2 应变功率谱密度传递比

对于平稳随机过程，系统的输入与输出有如下关系：

式中S(s)m×m表示响应的功率谱密度矩阵；m为位移响应的测试点数；G(s)n×n表示激励的功率谱密度矩阵；n表示激励的作用点数；H(s)为m×n阶频响函数矩阵；H*(s)T表示H(s)的共轭转置。

由方程（9），点i与点j的互功率谱可以表示为以下方程：

式中表示H(s) 第i行；（s）j表示矩阵H*(s)T的第j列。方程（10）可以写成如下式所示的累加求和形式：

式中Hik(s)表示i点响应和k点（k=1，2，…，N）激励之间的传递函数；Gkn(s)表示在点k和点n（n=1，2，…，N）上激励的互功率谱密度；(s)表示j点的位移响应和n点激励的传递函数的复共轭。

选定任一响应测试点p为参考点，i点和j点的功率谱密度传递比PSDT 定义为：i点响应yi（t）与参考点p点响应yp（t）之间的互功率谱Sip（s），j点响应yj（t）与参考点p点响应yp（t）之间的互功率谱Sjp（s），两者之比表示为：

借鉴位移响应功率谱密度传递比的概念，定义响应应变功率谱密度传递比（Strain Power Spectrum Density Transmissibility，SPSDT），即i点应变响应εi（t）与参考点p点应变响应εp（t）之间的互功率谱，j点应变响应εj（t）与参考点p点响应εp（t）之间应变响应的互功率谱密度两者之比表示为：

由方程（6）和（9），对于平稳随机过程，可以建立用应变表示的系统的输入与输出关系：

式中Sε(s)m×m表示结构应变响应的功率谱矩阵，m为响应测点；G(s)n×n表示激励的功率谱密度矩阵，n为激励作用点数；Hε(s)为m×n阶应变传递函数矩阵；Hε*(s)T为Hε(s)的共轭转置。由方程（9），（10），（11）和（14）可得：

则由方程（15）和（13），响应应变功率谱密度传递比可以表示为：

式中(s)表示k点激励，i点响应的应变传递函数；Gkn(s)表示在k点和n点的激励的互功率谱密度。在方程（16）中，取分母与分子有相同下标的项得到εtkn(s)：

当s→λr，εtkn(s)的极限：

式中和分别代表i测点和j测点在第r阶应变模态振型中对应的振型系数。在复平面中，当fn（s）和gn（s）（n=1，2，…，N）满足下式时：

方程（19）中η表示同一常数，可得：

结合方程（16），（18），（19）和（20）可以得到：

同理，当参考点选择为q点时，可以得到下式：

至此证明了在系统的极点处，SPSDT 等于测点i和测点j的应变振型系数和的比值，这一特性与参考点的选择无关，只与测点i和j的位置有关。

3 基于SPSDT 的应变模态参数识别

对测点i和j，选择不同参考点p和q，对SPSDT做差，由方程（21）和（22）可以得到：

当s→λr时：

对方程（24）取倒数记为Δ-1[(s)］，显然系统的极点是Δ-1[(s)]的极点。需要注意的是，由方程（25）所得到的极点中可能包含与系统极点无关的虚假极点。为了减少虚假极点的出现，融合所有参考点的信息构造如下方程，确定系统的极点：

选定k（k=1，2，…，N）测点为基准点，然后结合其他N个测点和N个参考点的信息，构成SPSDT矩阵εT(s)如下式所示：

由方程（21），（22）和（27），当s=λr时：

方程（28）中应变功率谱密度传递比矩阵εT(λr)的秩为1。在λ r处对εT(λr)矩阵进行奇异值分解：

方程（29）中所得到的左奇异矩阵U的第一列向量为第r阶应变振型。依次改变k（k=1，2，…，N）的取值，S表示对角线位置的元素由奇异值组成的对角矩阵，VH表示右奇异矩阵。由方程（27），（28）和（29）计算N个向量（k=1，2，…，N）。对以上所得N个向量按照最大值归一化以后取平均，即可得到应变振型如下式所示：

综上，基于响应应变功率谱密度传递比的结构应变工作模态参数识别过程如图1所示。

图1 基于SPSDT 的应变模态参数识别流程图Fig.1 Flow chart of strain modal parameter identification based on SPSDT

4 数值算例

如图2所示为三跨连续梁模型，连续梁总长为60 m，跨度为20 m+20 m+20 m。截面为矩形，宽0.6 m，高0.2 m。梁的材料特性如下：密度为7900 kg/m3，弹性模量为210 GPa。采用有限元方法，全桥划分为60 个单元，每个单元的长度为1 m，单元类型为BEAM188，通过自振分析可得到结构固有频率。

图2 三跨连续梁模型Fig.2 Model of 3-span continuous beam

用Matlab 生成61 组服从标准正态分布的随机数，每组数据包含5000 个数，模拟随机激励同步加载到各个节点，对梁做动力时程分析，时间步长为0.01 s，计算节点的应变时程数据。采用建立的基于响应应变功率谱密度传递比方法，仅基于应变响应识别梁的应变模态参数。按照最大值归一化，得到应变模态振型如图3所示，频率识别结果比较如表1所示。

表1 频率识别结果Tab.1 Recognition results of frequencies

值得注意的是，应变模态和位移模态是结构同一种状态的两种表达。通常位移模态振型φ和应变模态振型φε是不同的，如图3所示，图中位移模态振型由有限元计算直接得到，而应变振型为用本文方法识别得到的。

图3 位移振型和应变振型Fig.3 Displacement mode shapes and strain mode shapes

对于受弯为主的梁而言，应变是位移的二阶导数。有限元计算无法直接计算出应变模态振型，一般由位移模态振型差分2 次得到应变模态振型。本文所建立的应变功率谱密度传递比方法，基于结构的应变响应，可以直接识别出结构的应变模态振型。

5 实验验证

实验简支梁模型的梁长为3.05 m，净跨为3 m。矩形截面高为22.5 mm，宽为125 mm，材料为铝合金，密度为2.7×103kg/m3，弹性模量为70 GPa。梁等间距分为20 段，共设置19 个应变测点，测点布置如图4所示。实验采用的应变片为电阻式应变计，型号为BX120⁃3AA，由浙江黄岩测试仪器厂生产，电阻值为120 Ω，灵敏系数为2.08。实验数据采集使用NI 动态数据采集仪，该仪器包括的组件有NI PXIe⁃1082、PXIe⁃8840 控制器、PXIe⁃4330 应变测试模块、TB4330 接线盒。简支梁模型、动态数据采集仪、支座模型和应变片分别如图4～7 所示。

图4 简支梁模型Fig.4 Model of simply supported beam

图5 NI 数据采集仪Fig.5 NI data acquisition instrument

图6 支座模型Fig.6 Support model

图7 应变片Fig.7 Strain gauge

本文给出的SPSDT 模态参数识别法不需要对激励类型做任何假设，仅利用响应数据就能识别模态参数。因此实验室做试验时用小锤敲击作为加载，加载位置和敲击力大小随机。采样时长为150 s，采样频率为1000 Hz。测点10 的时长30 s 的应变时程曲线如图8所示。

图8 测点10 时长30 s 的应变时程Fig.8 30 seconds strain history of testing point 10

针对模型实验梁实测应变响应数据，采用本文SPSDT 方法进行了应变模态参数识别，并采用传统的频域分解法（FDD）和随机子空间方法（SSI）识别，同时给出有限元结果和理论解，五种方法所得固有频率结果比较如表2所示，结果吻合良好。

表2 频率结果比较Tab.2 Comparison of frequency results

由SPSDT 方法直接识别出的前四阶应变模态振型如图9所示。各振型按照最大值归一化处理，表明本文建立的基于响应应变功率谱密度传递比结构工作应变模态参数识别方法，仅基于应变响应，不仅能方便地识别出结构的模态频率，还可以识别出结构的工作应变模态振型。图9对比了解析应变振型和识别应变振型；表3给出了前四阶解析应变振型和识别的应变振型的MAC 值，其中对角项的MAC 均大于0.99，非对角项小于0.01。结果表明两种方法所得到的结果是一致的。

表3 解析应变振型与识别应变振型的MACTab.3 MAC of analysis strain mode shape and identifi‑cation strain mode shape

图9 简支梁应变振型Fig.9 Strain mode shapes of simply supported beam

为探讨数据长度（采样时长）对识别结果的影响，选取一段时长4 min，采样频率1000 Hz 的应变响应数据，将其划分成时长为1，2，3 和4 min 的四段数据，进行应变模态参数识别。4 种不同时长的应变数据识别出实验梁的固有频率比较如表4所示，所得前四阶频率结果相同。

表4 采样时长与频率Tab.4 Sampling time and frequency

对于应变模态振型的比较，采用模态置信准则（Modal Assurance Criterion，MAC）：

式（31）中，当MAC 为0 时，表示两向量完全无关；当MAC 为1 时，表示两向量完全相关。以采样时长60 s 的数据识别出的前四阶应变模态振型作为基准，和其他三组不同时长数据做比较。分别取2，3，4 min 时长数据所识别出的前四阶应变模态振型计算MAC，各组数据所识别结果不同阶次的应变模态振型之间的MAC 均小于0.01，相同阶次的应变振型MAC 接近于1。在表5中列出不同时长数据所识别的前四阶相同阶次应变模态振型之间的MAC值，其中第二列数值为采样时长1 min 的数据所识别的前四阶应变模态振型的MAC。

表5 采样时长与MACTab.5 Sampling time and MAC

6 结 论

（1）选取一个应变响应测点作为参考点，定义了一个响应应变功率谱密度传递比（Strain Power Spectrum Density Transmissibility，SPSDT），从理论上证明了SPSDT 在系统的极点处为应变振型系数之比；

（2）选取一系列不同的响应参考点构造响应应变功率谱密度传递比矩阵，在系统的极点处对该矩阵进行奇异值分解，分解所得左奇异矩阵的第一列向量即为应变振型，从而实现结构工作（仅基于响应）应变模态参数的识别；

（3）数值模拟算例和实验室模型试验验证了所提出方法的有效性，并分析了采样时长对识别结果的影响，结果表明该方法具有较好的鲁棒性；

（4）与传统的工作模态分析方法相比，SPSDT方法不需要对激励做白噪声假定，不需要多种激励类型，仅在一种激励下，即可识别出结构的工作应变模态参数。