基于贝叶斯公式的无参考信号Lamb波时间反转损伤识别方法

1、基于贝叶斯公式的无参考信号Lamb波时间反转损伤识别方法 缪傲（武汉大学土木建筑工程学院， 武汉 430072）摘 要 现有的主动Lamb波损伤监测大多采用基于参考信号的差信号方法获取损伤散射信号，然后用确定性的方法进行损伤定位，因此在适用性方面受到很大影响。本文引入Lamb波的时间反转过程，提出了一种无参考主动Lamb波损伤识别方法，且考虑损伤识别过程中的不确定性因素，通过贝叶斯公式，推导出损伤位置、波速等未知参数的联合后验概率分布后，利用蒙特卡洛马尔科夫链（MCMC）方法对未知参数进行采样估计，所得到的马尔可夫链的极限分布即为未知参数的后验分布。最后，通过在铝板上的数值仿真结果表明，该方法

2、能够较为准确的识别出损伤的位置、大致范围等特征。 关 键 词 损伤监测；时间反转；无参考；贝叶斯；MCMC1 前言V1V2V12V11V21V21V22Lamb波作为一种超声导波，在结构中衰减慢，能够传播较远的距离，并且对结构中的小损伤敏感，所以在结构健康监测领域中的应用越来越广泛。由于损伤会引起结构中传播的Lamb波的散射，故损伤散射信号直接与损伤位置有关。现有的主动Lamb波损伤监测方法大多基于参考信号，即以健康状态响应信号为基准，通过差信号的方法获取损伤散射信号，但真实结构和外界环境的变化对其影响很大，因此在适用性方面受到了很大的影响，且传统的基于损伤散射信号传播时间的损伤定位方法中，例

3、如，脉冲回波法、椭圆定位法、四点圆弧定位法等，把损伤散射信号传播时间、一定频厚积下的波速都作为确定值，但实际的量测误差和损伤定位中不确定性是不可避免的。例如，在损伤识别前，损伤的尺寸、程度是我们所不知道的，损伤对波速变化产生的影响也是不确定的，这将反过来影响通过理论计算所得到得损伤散射信号传播时间，再者，用做激励和传感的压电片的尺寸也没有考虑，这都给损伤的识别带来误差。除此之外，由于Lamb波的频散特性，损伤散射信号波包在传播过程中会发生畸变，从而也会对损伤散射信号传播时间产生影响。在这种情况下，概率的方法在损伤识别过程中能够考虑各种不确定性因素，因此，比用确定性的方法更合适。针对上述的诸多问

4、题，本文首先根据板结构中主动Lamb波时间反转理论，分析了存在单损伤板结构中Lamb波时间反转聚焦过程，提出了Lamb波无基准损伤散射信号传播时间的提取。然后从概率的角度，考虑损伤识别过程中模型本身和量测所存在的不确定性和误差，把利用时间反转所得到的传感器列阵中损伤散射信号传播时间D作为损伤位置、波速等未知参数的样本信息，再结合其先验信息，采用贝叶斯公式，得出未知参数的联合后验概率分布。最后通过MCMC方法，得到单个未知参数的后验分布。2基于贝叶斯方法的时间反转损伤定位2.1单损伤板结构中Lamb波的时反聚焦时间反转法是法国科学家Fink最先由光学应用引入到Lamb波领域中，并开展了大量的理论

5、和实验研究工作。Lamb波的时反是指将感器所接收到的各个模式的信号在时间域上反转后，在传感器上加载，即所接收到的信号先到后发，后到先发，从而实现了各个模式信号在原始激励处的自适应聚焦。已有的对板结构中传播的Lamb波传感信号时反特性研究表明：不管是在单模式还是双模式（S0，A0模式）下，Lamb波的频散都得到了补偿，同一模式经过前向和反向传播后都在主波峰处下聚焦，而在向前、反向的传播时，不同的群速度将会产生相应的的时间延迟，在主波峰旁辨形成旁瓣。为了降低信号处理的复杂度，选择合适中心频率的窄带信号，采用双面激励方式，以激发出单模式的Lamb波信号。如图a所示，对于单模式的Lamb波信号，由A传

6、播到B存在两条路径，PZTB所接收到的传感信号为直达波信号和损伤散射信号。如图b所示，将B处的传感信号经时域反转后再次加载在PZTB上，在PZTA处得到的时反重构信号可以表示为： 式（7）中，为B点接收到的直达波信号再次以直接波的方式到达A处，、的命名规则一样。根据式（5）可得： 将式（8）带入式（7）可得： 式（9）中基于压电片的激励模型，传递函数、可以简化为由幅值和速度的方程表为： 式中：为激励信号直线传播至传感器的路径，为激励信号经损伤散射后再传播至传感器的路径，为幅值频散方程，为散射系数，为模式的波数。考虑到此确定的模型中，传播距离和窄带激励信号的中心频率为确定值，故式（11）中的幅值

7、项也为确定值，为方便描述，。将式（11）带入式（9），并进行傅立叶反变换，可得时域内的重构信号为： 为了获得这个方程的近似解，通过泰勒级数在中心频率处展开得： 带入角频率、波数、相速度和群速度的关系、，上式变为 其中， 带入公式（14）、（15），公式（12）可表示为由式（17）可以看出：经过时间反转处理后，重构信号主要由三部分组成，第一部分为聚焦信号的主瓣，与初始激励信号的时反信号形状相同，实现了初始激励信号的重构，Lamb波传播过程中的频散得到了补偿；第二、三部分与主瓣部分有的时间延迟，对称的出现在主波峰旁辨形成旁瓣。同时，重构信号主瓣与旁瓣波峰对应的时间差正好为损伤散射信号传播时间于直达

8、波传播时间的时间差。在主动监测的传感器列阵中，直达波传感路径已知，故根据时反响应信号主、旁瓣之间的时间差就可以确定损伤散射信号的传播时间，因而可以实现无参考信号的损伤定位和识别。贝叶斯分析方法贝叶斯分析方法是基于假设的先验概率、以及把量测所得到的数据作为样本信息而得出的。由于考虑了先验信息和样本信息，故未知参数的不确定性减少。随着工程结构的不确定性等问题中概率方法的发展，贝叶斯方法在结构健康监测中起着越来越重要的作用。在损伤的评价、损伤预测、传感器优化布置、监测系统优化设计等工程上都得到了广泛的研究和应用。与一般确定性的方法给出未知参数的精确解相比，贝叶斯方法的基本理念是把未知参数向量的分布看

9、成随机变量联合分布，然后对每个参数做出在区间范围内的估计。具体为：将时间反转所得的Lamb波损伤散射信号传播时间D作为样本信息，综合未知参数的先验分布，再根据贝叶斯公式，得出后验信息，然后根据后验信息去推断未知参数。对于薄板结构，考虑具有条激励传感路径的压电片列阵，对于第（i=1,2）条传感路径的Lamb损伤散射信号传播时间可以示为lamb波从驱动器传播到损伤再到传感器所经历的时间。如图是一个简单的模型用来计算第条传感路径的Lamb损伤散射信号传播时间。在理论上，第i 条传感路径的Lamb损伤散射信号传播时间可以表达为： 式（18）中和分别是第条分别激励传感路径中损伤、驱动器、传感器位置的中心

10、坐标，是在一定的激励频率时lamb波损伤散射信号的传播波速。在式（18）中，压电片的尺寸、损伤的范围都没有考虑，材料属性的不确定性也会使实际波速偏离理论值，故在公式1中没有采用特定的频厚积下所对应的Lamb波波速，而是把波速也作为一个未知参数用来确定损伤的中心位置。因此在该损伤定位方法中，未知参数向量可以表示为。在下文中，为了方便，用来表示中的未知参数。在未知参数被确定后，损伤的定位则被转化为从概率的角度来解决确定未知参数分布的问题。假定建模误差和量测的不确定性分别用变量和表示，那么第条激励传感路径的损伤散射信号的传播时间从概率的角度可以表示为：是在未知参数向量下，利用公式（18）所计算得到的

11、第条传感路径所对应的损伤散射信号的传播时间。为了简单和方便，通常将和假定为平均值为0，标准差分别为和的相互独立的正态分布。由式（19）得： 变量服从正态分布，参数和的似然函数可以表示为： 式（21）中，方差，似然函数为在公式1和给定的未知参数向量下，对于所测得的损伤散射信号传播时间 D的概率分布。 为了方便，将公式3中的求和部分记为 由贝叶斯公式，根据所测得的TOF值，结合未知参数的先验概率密度函数和使然函数，未知参数的后验概率密度函数可以表示为 是 的联合后验分布函数， 是的联合先验分布函数。是一个常数，在后验分布中仅起到一个正则化因子的作用。对于每个参数的边缘概率密度函数 表示参数向量中对

12、除外的其它参数和做一个多维积分根据所测得的损伤散射信号传播时间，公式7给出了每个变量的边缘分布函数之后，根据后验分布函数，可以得到每个未知参数的后验分布，那么损伤的位置和波速度也就能够确定。但是在实际应用中，由于涉及到多个参数的高维积分，不便于用分析的方法计算，使得贝叶斯方法的应用大为受阻。运用MCMC方法，我们可以从复杂的多维积分分布中采样，从而能够对贝叶斯估计方法进行很好的模拟采样由于在贝叶斯方法中为一个常数，在计算后验分布中仅起到一个正则化因子的作用，如果将去掉，可将贝叶斯公式改写为可将贝叶斯公式改写为 其中符号表示左右两边相差一个常数因子式的右边虽然不是正常的概率密度函数，但它确是后验

13、分布的主要部分，故只要式的右边已知，就可以从后验分布中采样，从而避免了计算的问题。MCMC方法蒙特卡洛模拟中，在后验分布中取样，当这些样本独立时，利用大数定律样本均值会收敛到期望值。但在贝叶斯分析领域中，后验分布是复杂的、高维分布，很难得到独立的样本。目前MCMC方法已经成为一种处理复杂统计问题的特别流行的工具，尤其在经常需要复杂的高维积分运算的贝叶斯分析领域，通过合理的定义和实施，MCMC方法总能得到一条具有平稳分布的马尔科夫链，该马尔可夫链的平稳分布就是我们所需要的后验分布。下面给出M-H算法1、构造一个辅助的概率密度函数，通常被称为“提议函数”；2、在第n次迭代时，根据从提议函数中产生一

14、个新状态；3、计算接受概率 概率置我们称之为接受提议，以概率置，我们称之为拒绝提议。重复上述过程，当n足够时，知道得到的马尔可夫链收敛到一个平稳分布。在利用马尔可夫链进行抽样时，在收敛之前的一段时间，各个状态的边际分布还不能认为是稳定分布的，所以在进行估计之前，应该把前面的这n-1次迭代值去掉，这个过程称之为burn-in。提议函数可以是任意的分布形式，但是选择合适的提议函数得到的马尔可夫链能够更好地收敛，通常，选取的提议函数满足对称型，即，这将使接受概率变得更简单些。本文采用基于高斯分布和均匀分布的MCMC算法从后验分布中对未知参数向量进行采样。与传统的MH算法相比，本文所提供算法结合吉布斯

15、采样的理念：如果在其它参数值已知时，那么某个参数的条件分布可以作为后验分布来使用。根据其它变量的当前值，依次对分布的每个变量生成一个样本序列。在此算法中，引入优化参数L，调节调节接受概率，使得马尔可夫链经过适当的迭代次数后收敛。结合未知参数的先验分布、量测所得到得损伤散射信号传播时间，运用MCMC算法，可以得到联合概率分布中各个位置参数的后验分布概率。参数可以直接从伽马分布中进行采样，对于参数的估计，具体算法如下：设定总的迭代次数NT，在收敛之前的一段时间，其分布不能认为是稳定分布，所以在估计时，应该把前面NB的迭代值去掉。1、初始迭当n=0时，； 从中采样 2、n=1时， 从U(0,1)中产

16、生一个随机数u，当u<r时， ，调节参数，；否则 ，调节参数， 式中： k=1时，； k=2时 ， k=3时 ， ， 直接从分布中进行采样。重复步奏直到，当n超过Nb时，调节参数不再变化。舍弃前Nb次的迭代值，剩下的采样 值可以作为一条收敛的马尔可夫链，可以代表位 置参数的后验分布。 数值仿真结果本文设计了相应的验证系统，通过数值仿真方法对上述理论结果进行验证。仿真试件及压电原件布置情况如图所示：采用铝板结构, 杨氏模量72Gpa,密度2700kg/m3，泊松比0.3,其尺寸为580mm* 580mm*2mm，在其相同位置的上下表面布置压电片，构成双面激励方式。以铝板中心作为坐标原点，选

17、用4对压电片作为激励/传感器，其编号和坐标分别为：1#（-90，-90）、2#（90，-90）、3#（90,90）、4#（-90,90）。3数值仿真结果1 2 1 1 1 200mm 1 180mm 1 200mm 1 O仿真结果1 580mm 1 580mm 1 4数值仿真结果1 In the plate view In the thickness view (a)平面视图 (b)侧视图 铝板与PZT布置方案的几何模型 The geometric model of Aluminum plate and PZT layout scheme 应用ANSYS12.0进行建模，结构中损伤是以（20,

18、10）为圆心，直径为8mm的圆形通孔，采用删除有限元单元的方法模拟。划分网格时，单元形状选六面体，基体单元类型采用三维8节点实体单元solid45, 单元尺寸为1mm* 1mm*1mm，损伤周围局部加密，压电原件选用三维耦合场体单元solid5。为了保证压电层与金属面板重合部分传递力和位移的连续性,建模中采用黏结处理，模型的边界全都设置为自由。在仿真实验中，激励信号采用正弦调制窄带信号，表达式为：式中为激励信号的中心频率，取为100kHz,为阶梯函数，为调制信号的波峰数，取为5，激励信号的时频域特性如图所示。激励信号 (a)时域波形（b）频域图Load signal(a)Time domain

19、 input signal （b） Frequency spectrum 实验过程中，采用双面同相激励，垂直于板平面方向加载，此时产生的Lamb波以A0模式为主。逐个选取压电原件为激励源，其余的压电元件作为传感器，根据声波互易性原理，有6个通道的信号是独立的，表示为、 、(下标数字表示激励器和传感器的编号)，然后对采集到的各监测通道结构响应信号施加合适的窗函数进行截取，消除边界反射信号，仅保留直达波信号和损伤散射信号，并对截取后的信号进行时间反转并用作为二次激励，得到结构的时反聚焦信号。最后提取时各监测通道时反重构信号中的损伤损伤散射传播时间用于损伤识别。取通道相应信号进行分析Response signal of sensor波相应信号波形截取并作时间反转后的信号The signal after the interception and time reversal时反聚焦响应信号Reconstructed Time reversedsingal 根据式所表述的时反聚焦信号中主、旁瓣信号对应波峰的时间关系，提取各监测通道下的损伤散射信号传播时间，如表所示。对于监测通道 、，由于损伤距离直达波路径很近，信号散射角度大、传播距离远，加上Lamb波的频散效应导致的波包延拓和幅值衰减，损伤散射