基于改进模态柔度灵敏度和l

基于改进模态柔度灵敏度和l
p
正则化的结构损伤识别方法
技术领域
1.本发明属于结构健康监测领域，具体涉及一种基于改进模态柔度灵敏度和l
p
正则化的结构损伤识别方法。


背景技术：

2.由于环境侵蚀、过载、疲劳和其他因素，工程结构的局部区域在正常使用年限内可能会遭受损坏。结构损伤或劣化会导致其材料和几何特征的改变，这些会降低结构刚度并导致不利的振动响应。因此，发展有效的损伤检测方法，对结构损伤位置和程度进行可靠评估是非常重要的。
3.在过去的几十年里，基于振动的损伤识别技术得到了广泛的发展。模态参数与结构固有特性有关，不受外界激励的影响，研究人员利用模态参数建立了多种损伤识别策略。以往研究多采用结构的固有频率和模态振型进行损伤识别。然而，固有频率是全局动力特性参数，难以提供充足的局部损伤信息，现实中，结构的损伤却往往发生在极少数的局部位置。模态振型可表征结构的空间信息，从而对局部损伤更为敏感，但是，需要在结构上布置众多传感器来获取满足损伤识别要求的模态振型，且高阶模态振型在实际中难以测量。考虑到上述模态参数的局限性，另一种策略是利用模态柔度，模态柔度是基本模态参数的变体，并且对局部损伤更敏感。此外，利用模态柔度矩阵进行损伤识别的一个重要优点是只需测量少量的模态参数。
4.结构损伤识别问题通常被认为是不适定的，输入数据的微小波动可能导致解的巨大突变。而且，相对于整个结构的所有单元，少量的损伤位置体现出稀疏性。因此，损伤单元的识别同样也是一个稀疏恢复问题。根据压缩感知理论，l
p
正则化(0《p《1)可通过较少的测量数据就可得到更稀疏解的恢复，且噪声鲁棒性好。因此，为了解决以上两个问题，亟待提出一种基于模态柔度和正则化的结构损伤识别方法。


技术实现要素：

5.本发明通过扩展特征值灵敏度函数，进一步改进模态柔度灵敏度函数，引入l
p
正则化处理噪声环境下损伤检测的不适定问题，建立了一种基于改进模态柔度灵敏度和l
p
正则化的结构损伤识别方法。该方法较于传统的l1和l2正则化技术能够更精确的定位和量化结构的单损伤和多损伤；
6.本发明的目的是通过如下技术方案来实现的：该基于改进模态柔度灵敏度和l
p
正则化的结构损伤识别方法，包括如下步骤：
7.step1：通过在特征值导数中引入振型的变化改进了模态柔度灵敏度公式；
8.step2：通过数值模拟和试验获得结构的模态参数分别计算模态柔度矩阵，对结构损伤前后的模态柔度矩阵求残差建立目标函数，引入l
p
正则化进行优化求解；
9.step3：计算刚度损伤因子向量α，进而定位和量化结构损伤；
10.所述step1包括以下步骤：
11.step1-1，通过测试分别获得梁结构损伤前后的模态参数,由频率和振型分别计算模态柔度矩阵f；
[0012][0013]
式中，k为结构的刚度矩阵，λi＝(2πfi)2是结构第i阶特征值，fi是结构第i阶固有频率，是第i阶模态振型；
[0014]
step1-2，第i阶模态柔度基本灵敏度公式可表示为：
[0015][0016]
式中，p为结构参数，ω＝(k-λim)-1
，m为结构的质量矩阵，忽略损伤对结构质量的影响，令质量矩阵对损伤参数的一阶导数为0；
[0017]
step1-3，通过引入模态振型的变化对特征值灵敏度公式进行改进：
[0018][0019]
step1-4，令：
[0020][0021]
step1-5，最终，改进的模态柔度灵敏度公式可表示为：
[0022][0023]
如果只测量了前m阶模态参数，则改进模态柔度灵敏度公式可表示为：
[0024][0025]
具体的，step2测试获得的模态参数m不少于3阶。
[0026]
具体的，step2中l
p
正则化技术，选择p＝1/2，既可解决识别过程的不适定问题又可获得更精确的稀疏解。
[0027]
具体的，step3中包括如下步骤：
[0028]
step3-1，构造最小化目标函数，并在其中引入l
p
正则化项：
[0029][0030]
式中，是刚度损伤因子向量α的解范数，表示目标函数的残差范数，r
·
α＝sα+r0，其中，f表示结构的模态柔度矩阵；
[0031]
step3-2，基于l
p
正则化利用模态柔度进行损伤识别，优化公式为：
[0032][0033]
式中，fa和fe分别表示由数值模拟和试验得到的模态柔度矩阵，β是正则化参数，β＞0；
[0034]
step3-3，选取最优正则化参数β：
[0035]
绘制残差范数和解范数相对于与正则化参数β的曲线，选取两条“l”形曲线拐点所对应的β值作为最优参数；
[0036]
刚度损伤因子向量α的解范数：
[0037][0038]
目标函数的残差范数：
[0039][0040]
step3-4利用加权l1正则化求解公式，公式修改如下：
[0041][0042]
其中，权重ωi需要利用上一步刚度损伤因子αi迭代计算，即ωi＝|(αi|+ε)-1/2
；
[0043]
step3-5，结构损伤检测的l
p
算法实现如下：
[0044]
1).定义容差准则tol＝0.00001和迭代截止次数l
max
＝10；
[0045]
2).初始化权重：
[0046]
3).求解加权l
p
最小化问题：
[0047][0048]
4).更新权重：
[0049][0050]
式中，ε＝0.001，保证中的零值项不具有太大的权重；
[0051]
5).程序迭代终止于收敛准则满足或l＝l
max
；否则，设置l＝l+1，重复2～5步。
[0052]
具体的，step3中刚度损伤因子αi可定义为：
[0053][0054][0055]
式中，ne是结构单元数目，γi是刚度衰减系数，αi＝0表示结构中第i个单元未发生损伤，αi＝-1表示结构中第i个单元完全损坏，αi的值可表示结构的损伤位置和损伤程度。
[0056]
本发明的有益效果在于：本发明通过在特征值导数中引入振型的变化，建立了对
局部损伤更为敏感的模态柔度灵敏度表达公式，考虑到结构损伤分布的稀疏性，引入l
p
正则化算法用于结构稀疏损伤识别。该方法能有效对结构的单损伤和多损伤进行定位和定量分析，为结构的无损检测与评估提供了一种有效的新方法。
附图说明
[0057]
图1是本发明方法的流程框图。
[0058]
图2是本发明实例一中桁架结构示意图。
[0059]
图3是本发明实例一中噪声等级1下ds1残差和解范数相对于正则化参数变化曲线图。
[0060]
图4是本发明实例一中噪声等级1下ds2残差和解范数相对于正则化参数变化曲线图。
[0061]
图5是本发明实例一中噪声等级2下ds1残差和解范数相对于正则化参数变化曲线图。
[0062]
图6是本发明实例一中噪声等级2下ds2残差和解范数相对于正则化参数变化曲线图。
[0063]
图7是本发明实例一中无噪环境下ds1平面桁架损伤识别图。
[0064]
图8是本发明实例一中无噪环境下ds2平面桁架损伤识别图。
[0065]
图9是本发明实例一中噪声等级1下ds1平面桁架损伤识别图。
[0066]
图10是本发明实例一中噪声等级1下ds2平面桁架损伤识别图。
[0067]
图11是本发明实例一中噪声等级2下ds1平面桁架损伤识别图。
[0068]
图12是本发明实例一中噪声等级2下ds2平面桁架损伤识别图。
[0069]
图13是本发明实例一中ds1-e3下改进灵敏度公式与基本公式识别误差对比图。
[0070]
图14是本发明实例一中ds2-e10下改进灵敏度公式与基本公式识别误差对比图。
[0071]
图15是本发明实例一中ds2-e18下改进灵敏度公式与基本公式识别误差对比图。
[0072]
图16是本发明实例二中试验悬臂梁结构。
[0073]
图17是本发明实例二中悬臂梁的有限元模型。
[0074]
图18是本发明实例二中ds1残差和解范数相对于正则化参数变化曲线图。
[0075]
图19是本发明实例二中ds2残差和解范数相对于正则化参数变化曲线图。
[0076]
图20是本发明实例二中ds3残差和解范数相对于正则化参数变化曲线图。
[0077]
图21是本发明实例二中ds1悬臂梁损伤识别图。
[0078]
图22是本发明实例二中ds2悬臂梁损伤识别图。
[0079]
图23是本发明实例二中ds3悬臂梁损伤识别图。
[0080]
图24是本发明实例二中改进灵敏度公式与基本公式识别误差对比图。
[0081]
具体实施方法
[0082]
下面结合附图和实施例子对本发明做进一步的详细的描述。
[0083]
参见图1，是本发明基于改进模态柔度灵敏度和l
p
正则化的结构损伤识别方法的流程框图，其具体步骤如下：
[0084]
step1：通过在特征值导数中引入振型的变化改进了模态柔度灵敏度公式；
[0085]
step2：通过数值模拟和试验获得结构的模态参数分别计算模态柔度矩阵，对结构
损伤前后的模态柔度矩阵求残差建立目标函数，引入l
p
正则化进行优化求解；
[0086]
step3：计算刚度损伤因子向量α，进而定位和量化结构损伤；
[0087]
所述step1包括以下步骤：
[0088]
step1-1，通过测试分别获得梁结构损伤前后的模态参数,由频率和振型分别计算模态柔度矩阵f；
[0089][0090]
式中，k为结构的刚度矩阵，λi＝(2πfi)2是结构第i阶特征值，fi是结构第i阶固有频率，是第i阶模态振型；
[0091]
step1-2，第i阶模态柔度基本灵敏度公式可表示为：
[0092][0093]
式中，p为结构参数，ω＝(k-λim)-1
，m为结构的质量矩阵，忽略损伤对结构质量的影响，令质量矩阵对损伤参数的一阶导数为0；
[0094]
step1-3，一般的特征值灵敏度公式可以写成：
[0095][0096]
通过引入模态振型的变化对灵敏度方程进行改进：
[0097][0098]
基于振动理论可知，结构振型之间是相互独立的，一组完整模态振型可以用来定义第r阶振型的灵敏度系数，振型灵敏度函数定义如下：
[0099][0100]
将式(5)代入式(4)可得改进的特征值灵敏度函数：
[0101][0102]
step1-4，令：
[0103][0104]
step1-5，将式(6)、式(7)和式(2)组合，改进的模态柔度灵敏度公式可表示为：
[0105][0106]
如果只测量了前m阶模态参数，则改进模态柔度灵敏度公式可表示为：
[0107][0108]
具体的，step2测试获得的模态参数不少于3阶。
[0109]
具体的，step2中l
p
正则化技术，选择p＝1/2，既可解决识别过程的不适定问题又可获得更精确的稀疏解。
[0110]
具体的，step3中包括如下步骤：
[0111]
step3-1，构造最小化目标函数，并在其中引入l
p
正则化项：
[0112][0113]
式中，是刚度损伤因子向量α的解范数，表示目标函数的残差范数，r
·
α＝sα+r0，其中，f表示结构的模态柔度矩阵；
[0114]
step3-2，基于l
p
正则化利用模态柔度进行损伤识别，优化公式为：
[0115][0116]
式中，fa和fe分别表示由数值模拟和试验得到的模态柔度参数,β是正则化参数，β＞0；
[0117]
step3-3，选取最优正则化参数β：
[0118]
绘制残差范数和解范数相对于与正则化参数β的曲线，选取两条“l”形曲线拐点所对应的β值作为最优参数；
[0119]
刚度损伤因子向量α的解范数：
[0120][0121]
目标函数的残差范数：
[0122][0123]
step3-4利用加权l1正则化求解公式，公式修改如下：
[0124][0125]
其中，权重ωi需要利用上一步刚度损伤因子αi迭代计算，即ωi＝(|ai|+ε)-1/2
；
[0126]
step3-5，结构损伤检测的l
p
算法实现如下：
[0127]
1).定义容差准则tol＝0.00001和迭代截止次数l
max
＝10；
[0128]
2).初始化权重：
[0129]
3).求解加权l
p
最小化问题：
[0130][0131]
4).更新权重：
[0132][0133]
式中，ε＝0.001，保证中的零值项不具有太大的权重；
[0134]
5).程序迭代终止于收敛准则满足或l＝l
max
；否则，设置l＝l+1，重复2～5步。
[0135]
具体的，step3中刚度损伤因子αi可定义为：
[0136][0137][0138]
式中，ne是结构单元数目，γi是刚度衰减系数，αi＝0表示结构中第i个单元未发生损伤，αi＝-1表示结构中第i个单元完全损坏，αi的值可表示结构的损伤位置和损伤程度。
[0139]
具体的，在目标函数中引入传统l1正则化、l2正则化同l
p
正则化进行对比；
[0140]
在目标函数中加入l2正则化项：
[0141][0142]
其中，是参数向量α的解范数；表示目标函数的残差范数；β最优正则化参数；
[0143]
利用l曲线的拐点来选取最优正则化参数β；
[0144]
在目标函数中加入l1正则化项：
[0145][0146]
式中，β最优正则化参数；
[0147]
下面是应用本发明方法的具体工程实施例。以下的描述涉及附图时，除非另有表示，不同附图的相同数字表示相同或相似的要素。
[0148]
实施例一：参见图二，桁架结构共六跨，各个杆件长为2m。杆件截面尺寸为0.10*0.10m，弹性模量和材料密度分别为2.0
×
10
11
n/m2和2850kg/m3。在matlab中采用桁架单元建立有限元模型，有限元模型由13个节点和23个单元组成。表1给出了包括不同损伤的位置和严重程度在内的详细信息。
[0149]
表1平面桁架的损伤工况：
[0150][0151]
假定在各种损伤状态下采用完整的模态振型进行数值模拟。首先利用前三阶无噪声的模态参数用于损伤识别。然后，在模态参数上施加一些随机噪声电平以模拟真实的测量条件，并将噪声模态参数简单地模拟为：
[0152]
[0153][0154]
其中，κ是噪音水平；μ是一个均值为零、方差为1.0的随机量，服从正态分布；
[0155]
表2不同噪音水平：
[0156][0157]
具体实施步骤如下：
[0158]
步骤1：通过在特征值导数中引入振型的变化改进了模态柔度灵敏度公式。
[0159]
步骤2：通过数值模拟和试验获得结构的模态参数分别按式(1)计算模态柔度矩阵。
[0160]
步骤3：确定每种损伤条件的最优正则化参数，无噪声条件下，选择β＝0.005，其他噪声水平，通过式(12)、式(13)绘制残差范数和解范数与正则化参数β的曲线，来选取合适的正则化参数，对式(15)采用matlab中的序列二次规划(sqp)算法求解。不同正则化参数的残差和解范数曲线如图3～6所示。利用无噪声数据对基于l1、l2和l
p
正则化算法的损伤识别结果进行了对比，结果如图11～14所示。基于l1正则化和l
p
正则化的损伤模拟值与实际损伤严重程度接近，l1正则化利用无噪声模态数据具有更高的精度。但l2正则化方法的误差较大，特别是ds1中，损伤单元3(e3)的识别误差为13.09％。对基于l1、l2和l
p
正则化算法分别给出在两种不同噪声水平下各种损伤工况的识别结果，结果如图7～12所示。对于噪声水平1的情况，虽然l1正则化算法仍能识别损伤单元，但l1正则化算法对噪声更敏感，输出结果不太稀疏，特别是在多损伤情况下。随着测量数据中噪声的增加(噪声等级2)，l1和l2正则化算法已经无法准确识别严重程度较低的损伤单元。在工况2中，对于损伤严重程度为10％的单元3，l1算法产生12.08％的识别误差，而l2更是难以对单元3进行识别，误判单元2为损伤单元。相比之下，所提出的l
p
正则化性能相当好，它能对四种损伤情况下的损伤单元进行精确定位和定量，最大识别误差仅为5.16％。而且，与其它两种算法相比，l
p
正则化算法对噪声不敏感，在噪声干扰下仍能输出稀疏解。
[0161]
通过对比损伤工况的计算值与实际值间的相对误差，分析改进的灵敏度函数的优越性。如图13～15所示，各个工况下，使用基本公式识别损伤的相对误差都大于改进灵敏度函数。特别的是，当损伤程度较大时，这种误差更为明显。且可观察到使用改进的灵敏度识别误差基本都低于4％。对比结果证明了改进后的灵敏度函数的优越性。
[0162]
实施例二：梁的一端通过g形夹具固定在钢板上，模拟悬臂梁的固定边界条件，如图16所示。通过切割悬臂梁设置了不同的损伤严重程度和位置，损伤工况描述在表3中。
[0163]
表3悬臂梁损伤工况：
[0164][0165]
具体实施步骤如下：
[0166]
步骤1：通过在特征值导数中引入振型的变化改进了模态柔度灵敏度公式。
[0167]
步骤2：通过数值模拟和试验获得结构的模态参数分别按式(1)计算模态柔度矩阵。使用10个等距安装在横梁上的加速度计对结构进行振动测试，测试悬臂梁在三种不同损伤工况下的前六阶频率和相关振型。采用100个等长euler-bernoulli梁单元对悬臂梁进行模拟。每个单元的长度为10mm，与实验中的切割长度一致，可反映损伤位置的稀疏性且有效地量化损伤，悬臂梁试验概况如图17所示。
[0168]
步骤3：计算刚度损伤因子向量α，进而定位和量化结构损伤。通过式(12)、式(13)，确定每种损伤条件的最优正则化参数，对式(15)进行求解。图18-20示出了基于l
p
算法求解的残差范数和解范数相对于正则化参数β的曲线，双轴曲线拐点对应的横坐标值为最佳正则化参数。比较和分析了l1、l2和l
p
正则化算法在三种损伤工况下的差异，如图21-23所示。从这些图中可以看出，利用l2正则化识别出的损伤单元分布在许多位置，虽然可以近似检测到损伤位置，但损伤严重程度与实际情况相差较大。与l2正则化相比，l1和l
p
正则化性能更好。它们能够以更高的精度定位和量化悬臂梁的损伤。此外，所提算法具有更稀疏的识别结果，更符合实际情况。
[0169]
通过比较相对误差，证明了改进后的灵敏度函数的优越性。图24揭示了两种损伤工况的比较评估结果。结果表明，改进后的灵敏度公式比基本公式提供了更可靠、更准确的识别结果。从图中可以看出，改进的函数的最大识别误差不超过6％，而基本函数的相对误差明显较大。而且，随着损伤严重程度的增加，相对误差变化更加显著。
[0170]
以上所述仅为本发明的两个实施例,凡依本发明申请专利范围所做的均等变化与修饰,皆属于本发明的涵盖范围。

技术特征：
1.一种基于改进模态柔度灵敏度和l
p
正则化的结构损伤识别方法，其特征在于，包括如下步骤：step1：通过在特征值导数中引入振型的变化改进了模态柔度灵敏度公式；step2：通过数值模拟和试验获得结构的模态参数分别计算模态柔度矩阵，对结构损伤前后的模态柔度矩阵求残差建立目标函数，引入l
p
正则化进行优化求解；step3：计算刚度损伤因子向量α，进而定位和量化结构损伤；所述step1包括以下步骤：step1-1，通过测试分别获得梁结构损伤前后的模态参数,由频率和振型分别计算模态柔度矩阵f；式中，k为结构的刚度矩阵，λ
i
＝(2πf
i
)2是结构第i阶特征值，f
i
是结构第i阶固有频率，是第i阶模态振型；step1-2，第i阶模态柔度基本灵敏度公式可表示为：式中，p为结构参数，ω＝(k-λ
i
m)-1
，m为结构的质量矩阵，忽略损伤对结构质量的影响，令质量矩阵对损伤参数的一阶导数为0；step1-3，通过引入模态振型的变化对特征值灵敏度公式进行改进：step1-4，令：step1-5，最终，改进的模态柔度灵敏度公式可表示为：如果只测量了前m阶模态参数，则改进模态柔度灵敏度公式可表示为：2.根据权利要求1中所述基于改进模态柔度灵敏度和l
p
正则化的结构损伤识别方法，其特征在于：step2测试获得的模态参数m不少于3阶。3.根据权利要求1中所述基于改进模态柔度灵敏度和l
p
正则化的结构损伤识别方法，其特征在于：step2中l
p
正则化技术，选择p＝1/2，既可解决识别过程的不适定问题又可获得更精确的稀疏解。
4.根据权利要求1中所述基于改进模态柔度灵敏度和l
p
正则化的结构损伤识别方法，其特征在于，step3中包括如下步骤：step3-1，构造最小化目标函数，并在其中引入l
p
正则化项：式中，是刚度损伤因子向量α的解范数，表示目标函数的残差范数，r
·
α＝sα+r0，其中，f表示结构的模态柔度矩阵；step3-2，基于l
p
正则化利用模态柔度进行损伤识别，优化公式为：式中，f
a
和f
e
分别表示由数值模拟和试验得到的模态柔度矩阵，β是正则化参数，β＞0；step3-3，选取最优正则化参数β：绘制残差范数和解范数相对于与正则化参数β的曲线，选取两条“l”形曲线拐点所对应的β值作为最优参数；刚度损伤因子向量α的解范数：目标函数的残差范数：step3-4利用加权l1正则化求解公式，公式修改如下：其中，权重ω
i
需要利用上一步刚度损伤因子α
i
迭代计算，即ω
i
＝(|α
i
|+ε)-1/2
；step3-5，结构损伤检测的l
p
算法实现如下：1).定义容差准则tol＝0.00001和迭代截止次数l
max
＝10；2).初始化权重：3).求解加权l
p
最小化问题：4).更新权重：式中，ε＝0.001，保证中的零值项不具有太大的权重；5).程序迭代终止于收敛准则满足或l＝l
max
；否则，设置l＝l+1，重复2～5步。5.根据权利要求1中所述基于改进模态柔度灵敏度和l
p
正则化的结构损伤识别方法，其
特征在于：step3中刚度损伤因子α
i
可定义为：可定义为：式中，ne是结构单元数目，γ
i
是刚度衰减系数，α
i
＝0表示结构中第i个单元未发生损伤，α
i
＝-1表示结构中第i个单元完全损坏，α
i
的值可表示结构的损伤位置和损伤程度。

技术总结
本发明公开了一种基于改进模态柔度灵敏度和l


技术研发人员：殷新锋 晏万里 叶航艇 刘扬 唐盛华
受保护的技术使用者：长沙理工大学
技术研发日：2023.07.17
技术公布日：2023/9/20