冲击源定位是现代结构健康监测中的主要部分之一，通过定位可以确定工程结构中异常冲击源的位置．通常情况下，结构中的冲击事件是由外部物体与结构的撞击、结构构件的失效、裂缝发展等引起的．例如，高速运行的飞行器可能会与碎片等不明物件发生撞击，须要尽快锁定撞击位置，以确保飞行器安全稳定运行[1]．核电站工程中应用的松动件定位系统(loose parts monitoring system，LPMS)[2]可用来监测核反应堆一回路的螺钉、销钉等零部件脱落，以免导致反应堆安全性能下降．以上典型的冲击源定位方法是通过在结构的某些测点位置安装传感器，以捕捉和分析冲击源事件发生后弹性波经过结构传播产生的信号．此外，声发射领域(acoustic emission，AE)[3-4]会通过分析物体变形及破裂等过程产生的应力波信号的特征，以实现声源发射区域位置的定位．梅兆池等[5]针对管道泄漏这一问题，结合经验模态分解方法和小波分解重构的方法，对信号进行降噪及信号频域特征提取处理，对比分析不同孔径、不同管道压力下该方法的效果，检验了基于声发射及融合信号频域提取方法的管道泄漏检测定位技术的有效性．童国炜等[6]以位置坐标、声源发射时间作为未知变量，将其他含参量作为误差项的基础上，提出一种融合Geiger定位方法、同伦算法及牛顿梯度方法的源定位算法，并结合试验验证其算法对声发射源位置的高精度定位．目前，针对结构冲击源定位的方法主要是基于传感器获得各种特征信息，再利用信号处理手段实现定位，例如采用信号的幅值、信号的到达时间(或到达时间差)等常见的信号特征作为输入．通过设定目标函数，该目标函数涉及多个传感器处弹性波的到达时间差(time difference of arrival，TDOA)和方向相关波速，得到目标函数最小值下对应的源位置．Yin等 [7]提出了“Z”形星团作为“I”形星团的替代方案，利用连续小波变换估计弹性波的到达时间(time of arrival，TOA)，在未知材料性质的平板中实现了冲击源定位．同时，针对该问题还可利用希尔伯特包络线的第1个峰值作为信号到达点，来计算敲击平板时不同传感器间的时延，通过将小波降噪与希尔伯特包络算法相结合，实现平板上的碰撞定位．Park等[8]提出了基于时间反转的冲击源定位方法；Simone等[9]提出了由四个传感器组成传感器簇群的线性化定位技术；Sabzevari等[10]提出了一种基于波衰减分析的定位技术；Gollob等[11]提出了一种基于改进迪克斯特拉算法的定位方法；Park等[12]提出了一种基于波前形状的冲击源定位方法，用来求解材料性质未知的各向异性板，且对比频谱比方法，根据声源信号到达前后的波形数据来判断波的性质差别，从而定位声源位置．Seno等[13]开发一种基于人工神经网络(ANN)的传感定位方法，可模拟实际环境和操作条件下复合材料板中的冲击事件，使得在不同冲击条件下具有较高的准确度．以上研究中所提及的方法主要通过确定传感器接收到应力波的时间来进行冲击源定位，但是如果无法准确求解波的到达时间或者到达时间差，那么这些基于到达时间(到达时间差)的冲击源定位方法将导致定位结果失效．本研究根据传感器信号能量与传感器相对于冲击源位置的相关性，建立传感器信号能量与传感器相对于冲击源位置间的多参数模型，从而解算出目标函数最小值对应的冲击源位置坐标．用该方法进行结构冲击源定位时，可避免对波的到达时间(到达时间差)的高精度计算，且无须知道结构材料力学性质的先验信息．1 信号能量与传感器、冲击源间的相对位置关系图1给出了结构模型中冲击源和布置在结构模型上距离冲击源d处的传感器示意图．根据图1显示的坐标轴，冲击源位置(Dx，Dy)与传感器位置(x，y)之间的相对距离为d，相对夹角为θ，有：d=(x-Dx)2+(y-Dy)2；(1)θ=arctan(y-Dy)/(x-Dx)．(2)10.13245/j.hust.241215.F001图1冲击源位置与传感器相对位置图当外界的激励作用于结构模型时，会产生兰姆波[14]．当兰姆波在结构中传播时，其振幅与距离的平方根成反比[15]，可给出振幅A与距离d之间的关系式A(d)∝1/d    (d≠0)．(3)当兰姆波在结构中传播时，波在时间t和距离d处的振幅关系式[16]可表示为    A(d, t)=A0/d∙exp[i(kxx+ikyy-wt)]    (d≠0),(4)式中：A0为波源处的振动峰值，i2=-1；kx=kcos θ；ky=ksin θ；k=Re(k)+iIm(k)，Re(k)为波矢量k实数部，Im(k)为波矢量k虚数部．因此，式(4)可进行转换为    A(d,t)=A0/dexp{i[Re(k)(xcos θ+ysin θ)-ωt]-Im(k)(xcos θ+ysin θ)}    (d≠0).(5)当实际计算距离d处的峰值振幅a时，取A(d，t)的实部作为计算模量，    a(d)=|A(d, t)|=A0/d∙exp[-Im(k)(xcos θ+ysin θ)].(6)由于式(6)中x=Dx+dcos θ，y=Dy+dsin θ，因此式(6)中d处的a可进行转换为    a=A0/dexp[-Im(k)∙(Dxcos θ+Dysin θ+d)].(7)1.1　当角度恒定时信号能量与距离的关系当角度θ恒定时，兰姆波的最大振幅可表示为    a=a(d)=A0/dexp[-Im(k)(Dxcos θ0+Dysin θ0+d)]=A0/dek1exp[-Im(k)d].(8)由于θ为恒定值，可令常数θ=θ0，且Dx和Dy也为常数，因此可令k1=-Im(k)(Dxcos θ0+Dysin θ0)．根据式(8)，当d值较大，且大于5倍波长[17]时，兰姆波的最大振幅可表示为：a(d)=A0exp(k1)exp[-Im(k)'d];(9)exp[-Im(k)'d]=1/dexp[-Im(k)d].传感器可接收到从声源处传播的信号，从而计算出第i个传感器获得的能量，所以N个传感器信号的总能量为：E(d)=∑i=1NEi(d)=∑i=1Nai(d)2;(10)E(d)≈p1'exp(-q2'd),(11)式中p1'和q2'为常数．对式(11)两边取对数，有lnE(d)=lnp1'-q2'd=p1+qd，(12)式中：p1=lnp1';q=-q2'．1.2　当距离恒定时信号能量与角度的关系当d恒定时，对于sin θ和cos θ在θ=0°处进行泰勒展开，则对式(7)可进行转换，有：cos θ=∑m=0∞(-1)mθ2m(2m)!,sin θ=∑n=0∞(-1)nθ2n+1(2n+1)!;(13)    a=A0/dexp[-Im(k)(Dxcos θ+Dysin θ+d)]=A0dexp-Im(k)[Dx∑m=0∞(-1)m∙θ2m(2m)!+Dy∑n=0∞(-1)nθ2n+1(2n+1)!+d].(14)由于距离d为恒定值，可令常数d =d0，因此式(14)可转换为    a=a(θ)=A0/d0exp[-Im(k)d0]∙exp-Im(k)Dx∑m=0∞(-1)mθ2m(2m)!+Dy∑n=0∞(-1)nθ2n+1(2n+1)!.(15)传感器可接收到从冲击源处传播的信号，从而计算出第i个传感器获得的能量，所以N个传感器信号的总能量E(θ)可表示为E(θ)=∑i=1NEi(θ)=∑i=1N[ai(θ)]2．(16)结合式(15)和式(16)，可以发现能量E(θ)与常数θ部分的关系与式(9)和式(10)类似．此外，E(θ)中的指数部分涉及到关于θ的多项式，说明E(θ)是一个高度非线性的方程．为进一步优化该冲击源定位方程涉及的非线性问题，在E(θ)的表达式中只保留指数多项式的线性项，用常量参数分别替换常系数和指数项，从而可使须要优化的维度最小．最终，得到E(θ)的表达式E(θ)≈p2'exp(r1'+rθ)．(17)对式(17)两边取对数，有lnE(θ)=lnp2'+r1'+rθ=p2+rθ，(18)式中p2=lnp2'+r1'．1.3　信号能量与角度和距离之间的关系在实际的定位事件中，能量函数E会同时包含θ和d这两个变量．结合式(12)和式(18)，可得lnE=p+qd+rθ，(19)式中p=p1+p2．根据式(1)和式(2)，将d和θ分别用Dx和Dy代入式(19)，将其转化为    lnE=p+q(x-Dx)2+(y-Dy)2+rarctan(y-Dy)/(x-Dx).(20)2 基于灰狼优化算法的定位方法式(20)是含有5个未知参数的非线性方程，未知参数分别为Dx，Dy，p，q和r．因此，利用在结构模型上不同位置安装的N个传感器，可建立如下N个非线性方程组，其中每个方程对应唯一的一个传感器，    lnEi-p-q(xi-Dx)2+(yi-Dy)2-rarctan(yi-Dy)/(xi-Dx)=0    (i=1,2,…,N).(21)对于传感器i(i=1，2，…，N)，(xi，yi)为第i个传感器的坐标值，Ei为第i个传感器的能量值．通过建立方程组(21)，不断优化参数p，q，Dx，Dy和r，求解出目标函数σ的最小值，有    σ=∑i=1nlnEi-p-q(xi-Dx)2+(yi-Dy)2-rarctan(yi-Dy)/(xi-Dx)2.(22)根据目标函数σ的最小值结果，得到冲击源位置坐标(Dx，Dy)．2.1　建立灰狼优化算法模型灰狼优化(grey wolf optimization，GWO)算法[17] 通过模仿自然界中狼群的捕食行为模式，并且在此基础上依据灰狼群体的共同协作机制实现优化目标，可用于求解多目标优化问题．该算法具备结构简单、调节参数少、易于快速求解目标结果等优点．本研究利用灰狼优化算法，针对求解目标函数σ最小值对应的参数Dx和Dy．灰狼优化算法的基本思路如下．灰狼群体被严格划分为四层金字塔型的社会等级层次，分别为α狼、β狼、δ狼和ω狼，其中：α狼作为狼群的最顶层，是狼群的领导者，对β狼、δ狼和ω狼具有控制权；β狼作为狼群的次顶层，负责协助α狼，对δ狼和ω狼具有控制权；δ狼作为狼群的次底层，负责协助α狼和β狼，对ω狼具有控制权；ω狼作为狼群的最底层，所占比例最大，服从于α狼、β狼和δ狼．灰狼优化算法中的α狼、β狼和δ狼，代表最优的3个解，即最优解、次优解和较优解．α狼、β狼和δ狼指导ω狼向着目标搜索，具体搜索过程包含围困猎物、追剿猎物、攻击猎物和寻找猎物．2.2　灰狼优化算法的应用针对本研究的冲击源定位问题，灰狼优化算法的具体计算流程步骤如下．步骤1 初始化灰狼优化算法参数，包含最大迭代次数、种群数量、变量维度(未知参数为p，q，Dx，Dy和r)等；步骤2 根据实际应用场景给定的上下界条件，通过随机赋予灰狼群体个体的位置初始化X的值；步骤3 计算每一种狼的适应度σ值，并将狼群中适应度值最小(最优)的灰狼位置信息保存为Xα，将狼群中适应度值较小(次优)的灰狼位置信息保存为Xβ，将狼群中适应度值最差(较差)的灰狼位置信息保存为Xδ；步骤4 更新狼群中个体X的位置；步骤5 不断更新参数；步骤6 计算灰狼种群中灰狼个体的适应度值，同时更新三种头狼(α狼、β狼和δ狼)的最佳位置；步骤7 判断是否达到优化算法的最大迭代次数tmax，若符合最大迭代次数条件，则结束优化过程返回Xα值作为最优解，若不满足最大迭代次数条件，则转到步骤4～7．当达到最大迭代次数时，目标函数σ值将逼近最小值，此时可求得参数Dx和Dy，作为预测的冲击源位置坐标．目标函数σ值(适应度值)的变化曲线如图2所示，可以发现：随着迭代次数的增加，适应度值先是迅速降低，然后缓慢减小，随后收敛曲线保持较平稳的状态．具体表现为：当迭代次数为1～5次时，目标值呈现数量级式快速下降；当迭代次数为6～10次时，随着迭代次数的增加，目标值缓慢地降低，基本保持在同一数量级；迭代次数大于10次后，目标值的变化较为平稳，波动幅度较低．10.13245/j.hust.241215.F002图2适应度值变化曲线3 冲击源定位试验3.1　试验一(平板模型敲击试验)试验一通过在平板模型上安装加速度传感器，利用加速度传感器来接收力锤敲击平板的信号．本试验平台由测试对象、加速度传感器、数据采集设备及计算机等组成．测试对象钢板的材质为Q235，钢板的结构型式为正方体，其中钢板的具体尺寸为1.300 m×0.800 m×0.008 m．通过在钢板的四个角分别垫上橡胶隔振块以达到缓冲减振的效果．平面钢板上共布置8个加速度传感器(型号为B&K 4513B)．数据采集模块型号为B&K 3053-B-120，数据采样率设为65.536 kHz．图3为加速度传感器测点位置及敲击位置，其中红色方块为传感器测点，“X”形标记点为敲击位置．10.13245/j.hust.241215.F003图3平板测点及敲击点位置示意图试验一(平板模型敲击试验)中加速度传感器的位置坐标及敲击点位置坐标如下：S1~S8为传感器测点位置；1号～10号为敲击点位置．3.2　试验二(大尺度舱段模型敲击试验)试验二通过在大尺度舱段模型上布置传感器测点，利用加速度传感器来接收力锤敲击舱段模型的信号，如图4所示．试验二也由测试对象、加速度传感器、数据采集设备及计算机等组成．测试对象舱段的材质为Q235，舱段的结构型式为柱壳结构，其中舱段的具体尺寸为Φ3.00 m×8.80 m×0.01 m(即直径3.00 m，长度8.80 m，厚度0.01 m)．与试验一相比，本试验使用的舱段模型与平板的材质相同，结构型式不同．10.13245/j.hust.241215.F004图4舱段模型传感器位置及敲击位置示意图图4为舱段试验模型布置的平面展开图，其中大尺度舱段模型右下侧靠近舱门的加速度传感器命名为“S01”，其余传感器的编号自下而上、从右往左分别递增，舱段模型的外壳壳体上共布置53个加速度传感器(S01～S53)．图4中X轴方向的长度为8.65 m，Y轴方向的长度为9.42 m，肋骨间距为0.80 m．图4中红色方块(S01～S53)为传感器测点位置，蓝色“X”形标记点(1～11号)为敲击点位置．4 试验结果分析与讨论4.1　定位结果根据式(22)可知冲击源位置的理论计算坐标为(Dx，Dy)，但实际应用时目标值σ不一定完全等于0．因此，本研究利用灰狼优化算法，求解当目标值σ最小时，解算出式(22)中对应的参数Dx和Dy．本研究用定位误差(e)、定位误差度(f)、定位准确度(g)和相关系数(r)衡量定位结果．a．定位误差通过计算实际冲击源位置与定位位置之间的距离来评估定位误差，e=(Dx,act-Dx)2+(Dy,act-Dy)2．(23)b．定位误差度定位误差度即定位误差e与定位区域最大尺寸Lmax的比值，f=(e/Lmax)×100%．(24)c．定位准确度定位准确度即满足定位误差度要求(误差度＜20%)的敲击点数Np与总敲击点数N的比值，g=(Np/N)×100%．(25)d．相关系数[18]由式(20)可知ln E与θ和d之间存在较强线性关系，因此代入实际试验数据可建立因变量为ln E，自变量为θ和d的线性模型y'=ln E'=b0θ+b1d+b2,(26)式中y'(ln E')为该线性关系式的拟合值，利用线性回归得到系数b0，b1和b2的值．因此，可得到相关系数表达式r=∑(yi'-y¯)2/∑(yi-y¯)2．(27)对于n组数据(θ1，d1，ln E1)，(θ2，d2，ln E2)，…，(θn，dn，ln En)，y'为第i组数据的线性拟合值，y¯为n组数据的均值，yi为第i组数据的实测值．利用r建立ln E关于θ和d的线性方程来描述该线性方程与原始数据点之间的关联程度．本研究在平板试验模型上布置8个加速度传感器测点，共敲击10处位置，每个位置重复敲击2次．截取长度为500个采样点的数据信号．由于振动数据信号为加速度信号，因此对加速度信号进行二次积分，可得到位移历时曲线，再通过累计求和位移的平方计算出能量E的值，从而可求解出ln E的值．平板试验模型每个敲击位置的定位误差、相关系数、定位误差度及定位准确度结果见表1．10.13245/j.hust.241215.T001表1平板试验模型定位结果(Lmax=1.526 m)第1次敲击平板第2次敲击平板敲击点位置定位误差/m定位误差度/%敲击点位置定位误差/m定位误差度/%1号0.24816.231号0.22814.922号0.0523.402号0.0614.023号0.20313.293号0.22314.624号0.1127.374号0.15810.345号0.25016.375号0.27417.966号0.17311.356号0.21714.247号0.34022.277号0.29219.128号0.0946.138号0.1006.579号0.30419.939号0.31720.7710号0.0654.2810号0.1429.33由表1计算可知：第一次敲击平板(1～10号敲击点位置)，r=0.75，e=0.184 m，f=12.05%，g=90%；第二次敲击平板(1～10号敲击点位置)，r=0.65，e=0.201 m，f =13.19%，g=90%．本研究在大尺度舱段试验模型上布置53个加速度传感器测点，共敲击11处位置，每个位置重复敲击2次，截取长度为500个采样点的数据信号．由于振动数据信号为加速度信号，因此对加速度信号进行二次积分，可得到位移历时曲线，再通过累计求和位移的平方计算出能量E的值，从而可求解出ln E的值．大尺度舱段试验模型每个敲击位置的定位误差、相关系数、定位误差度及定位准确度结果见表2．10.13245/j.hust.241215.T002表2大尺度舱段试验模型定位结果(Lmax=12.79 m)第1次敲击舱段第2次敲击舱段敲击点位置定位误差/m定位误差度/%敲击点位置定位误差/m定位误差度/%1号4.21332.941号3.47427.162号4.27133.392号3.30225.823号1.43611.233号1.2259.584号1.88214.714号1.75113.695号1.85514.505号1.76213.786号0.6915.406号0.9127.137号1.31810.307号1.38810.858号0.6825.338号0.6384.999号1.0017.839号1.1118.6910号1.1869.2710号1.34010.4811号4.76137.2211号3.57127.92由表2计算可知：第一次敲击舱段(1～11号敲击点位置)，r=0.41，e=2.118 m，f =16.56%，g=72.73%；第二次敲击舱段(1～11号敲击点位置)，r=0.43，e=1.861 m，f=14.55%，g=72.73%．4.2　定位误差分析通过选取不同长度的数据点，可得到表3中关于平板和舱段的定位误差结果．10.13245/j.hust.241215.T003表3定位结果分析试验模型数据长度re/mf/%第1次敲击平板5000.750.18412.061 0000.530.21914.352 0000.420.28318.554 0000.430.23315.278 0000.380.29319.20第2次敲击平板5000.650.20113.171 0000.510.22714.882 0000.460.27818.224 0000.450.29619.408 0000.400.31720.77第1次敲击舱段5000.412.11816.561 0000.382.33518.262 0000.362.49619.524 0000.352.53219.808 0000.332.59720.30第2次敲击舱段5000.431.86114.551 0000.412.08916.332 0000.362.38818.674 0000.372.39618.738 0000.342.46319.264.2.1　定位误差结果与相关系数的关系结合表3中的数据，可得到平板、舱段模型的定位结果数据，绘制如图5和图6所示定位误差、定位误差度与相关系数的拟合结果图．从图5和图6可以看出：平板模型和舱段模型的定位误差、定位误差度与相关系数整体表现为线性负相关，随着相关系数的增大，定位误差及定位误差度逐渐减小；定位误差(定位误差度)与相关系数的线性拟合度较高，说明定位误差(定位误差度)与相关系数的关联程度较高．10.13245/j.hust.241215.F005图5平板试验模型定位误差和定位误差度10.13245/j.hust.241215.F006图6舱段试验模型定位误差和定位误差度4.2.2　相关系数与数据长度的关系结合表3中的数据，得到平板、舱段定位结果数据，绘制如图7所示相关系数与截取数据长度的关系曲线(数据长度即为用于定位计算的数据采样点数)．10.13245/j.hust.241215.F007图7相关系数与数据长度的关系曲线根据图7可知：随着数据长度的增加，平板及舱段试验模型的相关系数总体上呈现降低趋势，且当数据长度在2 000个采样点以内时，相关系数的下降程度显著高于数据长度为2 000～8 000个采样点．这是由于波在边界会发生反射，随着数据信号长度的增加，用于计算的数据点并不完全是入射的兰姆波信号，从而造成冲击源定位的误差增大，导致相关系数降低．结合图7，整体上来看，平板试验模型的相关系数高于舱段试验模型的相关系数．这可能是由于舱段试验模型实际是圆柱形曲面结构，而相关系数的计算是基于平面结构模型展开的，从而造成舱段试验模型的相关系数偏低．4.2.3　定位误差与数据长度的关系结合表3中的数据，绘制定位误差、定位误差度与截取数据长度的关系曲线如图8所示．10.13245/j.hust.241215.F008图8定位误差、定位误差度与数据长度关系曲线根据图8可知：平板及舱段试验模型的定位误差会随着数据长度的增加呈现先明显上升后趋于平缓的变化趋势，平板试验模型的定位误差小于舱段试验模型的定位误差；平板及舱段试验模型定位误差度与数据长度的变化曲线较为一致，随着数据长度的增加，定位误差度均表现为先上升后平缓变化的趋势；平板及舱段试验模型的定位误差度区间较为一致，均处于10%～22%之间．5 结论本研究基于波在结构中的传播特性，从理论上推导出信号能量与传感器相对于冲击源位置存在线性相关性，结合平板模型和大尺度舱段模型敲击试验，运用灰狼优化算法预测出冲击源位置，分析结构冲击源定位误差、定位误差度与相关系数、数据信号长度之间的联系，得到结论如下．a．通过计算传感器信号能量与传感器相对于冲击源位置的相关性，建立传感器信号能量与传感器相对于冲击源位置之间的多参数模型，从而解算出目标函数最小值对应的冲击源位置坐标．b．运用灰狼优化搜索算法，可较精准求解目标函数最小值及对应的冲击源位置坐标．灰狼优化算法具有较强的收敛性能、参数较少、容易实现等特点．c．平板和舱段试验模型的定位误差、定位误差度与相关系数呈现出较强的线性负相关关系，随着相关系数的增大，定位误差及定位误差度总体上逐渐减小．d．由于波的反射，平板及舱段试验模型的相关系数随着数据信号长度的增加表现为下降趋势．当数据长度在2 000个采样点以内时，相关系数的下降程度显著高于数据长度为2 000～8 000个采样点的相关系数，平板试验模型的相关系数高于舱段试验模型的相关系数．e．冲击源定位的预测精度随数据信号长度的增加而降低，平板和舱段试验模型的定位误差及定位误差度会随数据信号长度的增加，先上升后趋于平缓变化；平板试验模型的定位误差小于舱段试验模型的定位误差；平板和舱段试验模型定位误差度的变化区间保持一致，均位于10%～22%之间．