CN101865758B - 基于多重信号分类算法的冲击载荷定位方法 - Google Patents

Abstract

本发明公布了一种基于多重信号分类算法的结构冲击载荷定位方法，包括如下步骤：步骤一：冲击响应信号的采集；步骤二：利用多重信号分类(MUSIC)算法进行冲击源的波达方向(DOA)估计；步骤三：根据Lamb波传播特性进行冲击源距离的估计。本发明方法简便可靠，可对冲击进行实时、在线监测，冲击定位精度高、实时性强。

Description

基于多重信号分类算法的冲击载荷定位方法

技术领域

本发明基于阵列信号处理技术中的MUSIC算法的结构冲击载荷位置识别方法，属工程结构健康监测方法。

背景技术

复合材料以其比强度高、比刚度大、抗疲劳性能好及材料性能可设计等一系列优点，在航空、航天、汽车等工程领域得到了日益广泛的应用，尤其是在军用、民用飞机上已开始越来越多地使用先进复合材料结构。然而层合板复合材料在服役过程中不可避免的要承受各种能量物体的冲击，比如，冰雹的撞击、飞行器翼面与空中飞鸟的碰撞、飞行器受到枪击及工具经常在维护过程中掉落在飞行器表面等。这些冲击极易造成复合材料结构的内部分层、基体开裂和纤维断裂等损伤。这些内部损伤将使层合结构的力学性能严重退化，强度可削弱35％～40％，导致承载能力大大降低，对结构的整体破坏和失效形成潜在的威胁。而且这些损伤多发生在材料内部不易从表面发现，留下严重的隐患，使得具有损伤的复合材料结构具有突发性和灾难性失效的潜在能力，并且它在损伤、失效等方面的表现却是机理复杂，现象多样，判别困难。因此很有必要对复合材料结构进行全寿命的监测，以确保结构的稳定性和安全性。目前已有许多传统的已被广泛应用的无损检测技术，例如敲击、超声、X射线、电涡流射线、电位测量以及热应力场等方法。但是这些检测方法一般设备复杂、耗时耗力，需要对损伤的位置有初步的了解，使用不方便，局限性大，不易做到服役环境下的实时在线监测，不适合未来大型航空、航天飞行器结构的健康监测与诊断。对于在服役过程中的飞行器，重要的是如何来实时监测冲击试件以保护飞行器的安全运行。而利用基于压电智能结构的健康监测技术可以实现低速冲击的实时、在线监测，评估结构的损伤程度，从而保证结构的稳定性和安全性以及降低结构的维护成本。

目前冲击载荷识别的研究主要集中在实验室针对小试件进行，多是在模拟冲击载荷的情况下进行的实验研究，而且实验对象也多以结构简单的梁或板结构为主，很少在复合材料结构上进行真实的冲击实验。尽管已有各种冲击载荷识别方法，但是目前还没有形成一种能够高精度、实时在线识别冲击载荷的方法。对于复合材料结构的冲击定位问题，原理上可以利用声发射技术进行定位，但由于复合材料的各向异性和应力波在传播过程中的反射、折射、散射及噪声的影响，很难准确确定波速及波达时间，因此单纯采用常规的声发射技术很难对冲击源进行精确定位。而基于优化技术的识别方法在反演过程中需要不断地进行结果迭代，求解时间过长。因此，对于范围稍大的监测区域来说，搜寻最优位置所需的时间更长，单纯运用最优化方法进行实时在线监测已不大可能。并且基于优化技术的识别方法通常直接将时域内许多离散点的冲击力时间历程本身当作未知参数进行求解。由于要求解的未知参数过多，计算任务重，也很难满足在线监测的实时性要求。基于人工智能的荷载识别方法可以解决荷载识别问题，但用于荷载识别的模型需要以实际加载方式或解析解求出响应样本的方法来训练，只能用于已知结构而且对结构的训练工作量较大，训练后的结构在实际荷载识别时亦存在不稳定现象，因此，应用于冲击监测这一实际问题仍有一定的难度。

阵列信号处理技术作为信号处理领域一个重要的分支，近几十年来已在雷达、声纳、通信、地震勘探、射电天文及医学诊断等多种国民经济和军事领域广泛应用。其与传统的单个或少量定向传感器相比，阵列信号处理具有灵活的波束控制、高的信号增益、极强的抗干扰能力及高的空间超分辨能力等优点。MUSIC是Schmidt和Bienvenu及Kopp于1979年提出，由Schmdit在1986年重新发表。MUSIC算法的提出促进了特征结构分类算法的兴起和发展，该算法已成为空间谱估计理论体系中的标志性算法。由于它能对多个空间信号进行识别，故这种方法称为多重信号分类法。由于MUSIC方法发展成熟，分辨率高，对阵列结构要求低等特点，以及超高速超大运算能力信号处理芯片的出现使其在阵列信号处理中的应用非常广泛。因此，将阵列信号处理技术中的MUSIC算法引入到结构健康监测领域中，研究冲击源的定位问题。

发明内容

技术问题：本发明的技术问题是提供一种利用阵列信号处理技术中的MUSIC算法和Lamb波频散特性的结构冲击载荷位置识别方法。本方法利用MUSIC方法实现冲击源波达方向的估计，根据Lamb波的频散特性实现冲击源距离的估计，可以对冲击实时、在线监测，准确、快速地识别出冲击载荷位置。

技术问题：本发明为实现上述目的，采用如下技术方案：

本发明基于多重信号分类算法的结构冲击载荷定位方法，其特征在于包括如下步骤：

步骤一：冲击响应信号的采集

(1)在所述结构上布置由多个相同型号的压电元件组成的压电阵列；

(2)将2块4通道的数采卡0通道并联在一起，保证两块数采卡同步触发采集；

(3)设定采集参数，利用电荷放大器和数据采集卡采集压电阵列输出的冲击响应信号；

步骤二：基于MUSIC算法进行冲击源的波达方向估计

首先用快速傅里叶变换对所述冲击响应信号进行频谱分析，然后根据选定的方向估计中心频率利用小波变换从所述冲击响应信号中提取出中心频率下的窄带信号，最后取窄带信号的幅值最大的信号运用多重信号分类(MUSIC)算法的进行DOA估计；

步骤三：基于Lamb波传播特性进行冲击源的距离估计

首先用小波变换从任一压电元件接收到的冲击响应信号中提取出距离估计中心频率下的Lamb波信号，获得同一频率的对称S₀模式的Lamb波和反对称A₀模式的Lamb波，用最大峰值法求出对称S₀和反对称A₀这两种模式的Lamb波到所述达压电元件的时间差，最后由测得的这两种模式的Lamb波的传播速度求出冲击源到达压电元件的距离，将各压电元件估计出的距离值转化为相对于压电阵列参考点的距离，取所有压电元件测得的距离的平均值作为距离估计值。

优选地，步骤一中，所述压电阵列中压电元件间的距离小于冲击源信号波长。

有益效果：本发明与现有技术相比，本方法简便可靠，可对冲击进行实时在线监测，实时性强、定位精度高，能实际应用于工程结构中冲击载荷的监测及定位，确保工程结构的运行安全，防止事故的发生。

附图说明

图1是一均匀线阵在窄带远场条件下是我信号模型的示意图。它是由M个阵元等间距排列成一条直线的均匀线阵，以左边第一个阵元为参考点，以垂直于阵列的方向为法线方向。N个远场的窄带信号以θ角入射到此均匀线阵上。

图2是铝板上的一组Lamb波频散曲线。纵坐标为群速度，横坐标为板厚与频率的乘积(频厚积)。

图3是压电阵列测得的复合材料结构上的冲击响应信号。

图4是截取出的0.4-0.6ms时间段的小波信号波形。

图5是用MUSIC方法估计出的方位谱。

图6是用小波变换提取出的Lamb波，小波变换的中心频率为25kHz。

图7是本发明方法流程图。

具体实施方式

下面结合附图对发明的技术方案进行详细说明：

如图7所示，本发明的方法包括如下步骤：

(1)冲击响应信号的采集：

包括如下步骤：

(1.1)在结构上布置一个由7个相同型号的压电元件组成的均匀线阵；

(1.2)将2块4通道的数采卡0通道并联在一起，保证两块数采卡同步触发采集。

(1.3)设定采集参数，如采样频率、触发预置、预采点数、采集点数等。利用电荷放大器和数据采集卡采集冲击响应信号。

步骤二：基于MUSIC算法进行冲击源的波达方向估计。

MUSIC方法在窄带远场条件下可实现信号源的DOA估计。但冲击响应信号频率成分十分复杂，直接测得的信号不符合窄带信号模型要求，因此本发明实现DOA估计的基本思路就是：首先利用小波变换对由压电传感器阵列单元直接测得的冲击响应信号进行预处理，提取出利用MUSIC算法进行DOA估计所需要的窄带信号。具体做法就是：首先用快速傅里叶变换对由压电传感器阵列直接测得的冲击响应信号进行频谱分析，然后选定某一中心频率，利用小波变换从压电传感器阵列直接测得的冲击响应信号中提取出此中心频率下的窄带信号，最后取此提取出的响应波信号的某一时间段内的幅值较大的信号运用多重信号分类(MUSIC)算法的进行DOA估计。下面结合MUSIC方法的基本原理对本发明实现DOA估计的方法进行说明。

MUSIC算法的基本思想就是则是将任意阵列输出数据的协方差矩阵进行特征分解，从而得到与信号分量相对应的信号子空间和与信号分量相正交的噪声子空间，然后利用这两个子空间的正交性来估计信号的参数(入射方向、极化信息及信号强度等)。

在窄带远场条件下，N个信号源以入射角θ入射到由M个压电阵元组成的均匀线阵上，如图1所示，DOA估计的数学模型为：

x(t)＝Bs(t)+Noise(t)             (1)

式中，x(t)为阵列的M×1维数据矢量，Noise(t)为阵列的M×1维噪声数据矢量，s(t)为空间信号的N×1维矢量，M和N分别为阵元数目和信号源数目，B为空间阵列的M×N维的流型矩阵(导向矢量阵)，且

B＝[b(θ₁)b(θ₂)…b(θ_N)]          (2)其中，导向矢量

b(θ_i)＝[1 exp(-jω₀τ_2i)…exp(-jω₀τ_Mi)]^T，i＝1，2，…，N    (3)其中，ω₀＝2πf，f为入射信号频率。

下面根据窄带远场信号的数学模型对MUSIC算法实现DOA估计的过程进行分析。由窄带远场信号的数学模型式(1)可知阵列输出数据的协方差矩阵为：

R＝E[xx^H]＝BE[ss^H]B^H+BE[sNoise^H]+E[Noises^H]B+E[NoiseNoise^H]   (4)

由于假设信号与噪声是不相关的，且噪声为平稳的加性高斯白噪声，因此上式中的二、三项为零，上式可化简为：

R＝E[xx^H]＝BE[ss^H]B^H+E[NoiseNoise^H]＝BR_SB^H+R_N  (5)

式中，R_S，R_N分别为信号协方差矩阵和噪声协方差矩阵，设噪声功率为σ²，则有下式：

R＝BR_SB^H+R_N＝BR_SB^H+σ²I                        (6)

由于信号与噪声相互独立，数据协方差矩阵可分解为与信号、噪声相关的两部分，其中R_S是信号的协方差矩阵，BR_SB^H是信号部分。对R进行特征值分解有：

$R=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\λ}}_i{V}_i{V}_i^H+\underset{j=N+1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\λ}}_j{V}_j{V}_j^H=U_S{\mathrm{\Σ}}_S{U}_S^H+U_N{\mathrm{\Σ}}_N{U}_N^H---\left(7\right)$

式中，特征值满足关系λ₁≥λ₂≥…λ_N≥λ_N+1＝…＝λ_N＝σ²。

信号协方差矩阵R特征值分解得到的特征值即为信号功率，在信号功率大于噪声功率的条件下大特征值就对应信号子空间，小特征值对应噪声子空间。因此U_S是由大特征值对应的特征矢量张成的信号子空间，而U_N是由小特征值对应的特征矢量张成的子空间即噪声子空间。

M个特征值对应的特征向量分别为V₁，V₂，…，V_N，V_N+1，…，V_M，其中前N个对应大特征值，后M-N个对应小特征值。对于M-N个小特征值对应的特征向量，有

(R-λ_iI)V_i＝0     i＝N+1，…，M    (8)

即：

$(R-{\stackrel{\‾}{\mathrm{\σ}}}_N^{2}I)V_i=({\mathrm{BR}}_S{B}^H+{\mathrm{\σ}}_N^{2}I-{\mathrm{\σ}}_N^{2}I)V_i---\left(9\right)$

$={\mathrm{BR}}_S{B}^H{V}_i=0,i=N+1,\·\·\·,M$

因为B满秩，R_S非奇异，则有

B^HV_i＝0     i＝N+1，…，M          (10)

上式表明，矩阵B的各个列向量与M-N个最小特征值对应的特征向量正交，也即是矩阵B的各个列向量与噪声子空间正交，故有(11)式。由于R_S非奇异，则阵列方向矩阵B，与阵列输出向量的协方差矩阵的N个大特征值对应的特征向量所张成的信号子空间U_S相同，也就是信号子空间U_S和噪声子空间U_N正交。

$U_N^{H}b\left(\mathrm{\θ}\right)=0---\left(11\right)$

MUSIC算法正是基于上述性质提出的，但在实际应用中，理想的协方差矩阵是未知的，并且无法得到，所谓的协方差矩阵R是利用有限个数据向量X(1)，X(2)，…，X(L)估计得到的，数据协方差矩阵的最大似然函数估计为

$\hat{R}=\frac{1}{L}{\mathrm{XX}}^H---\left(12\right)$

式中，X＝[X(1)，X(2)，…，X(L)]，X为数据向量矩阵，L为采样数据的段数。

对

进行特征分解可以计算得到噪声子空间在实际情况下，存在的噪声与信号源并不是完全不相关，信号方向与噪声方向并不完全正交，即b(θ)与

并不能完全正交，也就是式(11)并不严格成立。因此，实际上求DOA是以最小优化搜索实现的，即

${\mathrm{\θ}}_{\mathrm{MUSIC}}=\arg \underset{\mathrm{\θ}}{\min }b{\left(\mathrm{\θ}\right)}^H{\hat{U}}_N{\hat{U}}_N^{H}b\left(\mathrm{\θ}\right)---\left(13\right)$

式中，arg min表示函数值最小时所对应的角度。于是可以定义MUSIC算法的空间谱函数：

$P_{\mathrm{MUSIC}}=\frac{1}{b^H\left(\mathrm{\θ}\right){\hat{U}}_N{\hat{U}}_N^{H}b\left(\mathrm{\θ}\right)}---\left(14\right)$

根据正交性，上式会在信号DOA值附近出现峰值，与“谱峰”对应的θ，即为DOA的估计。

由上述分析可知，利用MUSIC算法实现冲击源的DOA估计的步骤可归纳如下：

(1)根据公式(12)得到阵列接收数据的协方差矩阵估计

(2)对

进行特征值分解；

(3)确定信号子空间

和噪声子空间

(4)根据信号参数范围由式(14)进行谱峰搜索；

(5)找出极大值点对应的角度就是信号源的DOA。

下面对利用MUSIC算法实现冲击源DOA估计的所需的运算量进行分析，主要包括三部分：首先需要估计阵列协方差矩阵，其运算量为O(M²L)；其次对阵列协方差矩阵作特征值分解，其运算量为O(M³)；最后就是“谱峰”搜索需要的运算量。对于本发明的单冲击源的情况，阵元数目为7，由运算量分析可知，由于阵元数目较少，MUSIC算法估计DOA所需的运算量也就很小，因此该方法具有较高的实时性。

步骤三：基于Lamb波传播特性进行冲击源的距离估计。

在无限均匀各向同性弹性介质中，只存在两种波——横波和纵波，二者分别以各自的特征速度传播而无波形耦合，而在板中则不然，在板的某一点上激励声波，由于声波传播到板的上、下界面时会发生波形转换，经过在板内一段时间的传播之后，因叠加而产生波包，即所谓的板中Lamb波模态。Lamb波在板中传播时，存在不同的模式，各种模式的叠加效果即为Lamb波。本质上，Lamb波是在具有两个平行表面的结构中由横波和纵波相互耦合而成的一种弹性波。Lamb技术是目前研究较多的一种板类结构损伤监测技术。板波理论认为，对于工程上大量使用的板状结构，机械波在板中主要激励起扩展波(最低阶对称波)和弯曲波(最低阶反对称波)两种模式。薄板结构在承受法向冲击产生的板波声发射信号从本质上来说是一种Lamb波。

Lamb波在板中传播时，板中质点产生振动，其振动方式十分复杂，随各种参数(频率、板厚等)的变化而改变。根据薄板两表面质点的振动相位关系，可以把Lamb波的模式分为对称波和反对称波两种，对称波又分为S₀，S₁，…，S_n等多个模式，反对称波又分为A₀，A₁，…，A_m等多个模式。

Lamb波传播特性极为复杂，主要体现在其频散特性和多模现象。所谓的频散特性就是Lamb波的速度随着频率的变化而发生改变。而多模现象就是同一个频率值对应着不止一种Lamb波模式。根据板波理论，对于薄板结构，冲击产生的Lamb波信号，主要是S₀和A₀两种模式的Lamb波，而且在较低的频厚积范围内同一频率的S₀模式Lamb波比A₀模式Lamb波传播速度快，一组铝板中传播的Lamb波群速度频散曲线如图2所示。因此，可根据同一频率的这两种模式的Lamb波的传播速度不同求出这两种模式的Lamb波到达同一阵列单元的时间差，再根据测得的这两种模式的Lamb波的传播速度即可实现距离估计。

薄板结构因受到冲击而产生的Lamb波信号极为复杂，含有较多的模式和丰富的频率成分。这就使直接利用冲击响应信号实现距离估计变得非常困难。因此在估计冲击源的距离时，也同样用小波变换提取出第l(l＝0，…，M-1)个压电阵元的某一中心频率下的Lamb波响应信号，可以获得同一频率的对称S₀模式的Lamb波和反对称A₀模式的Lamb波，然后根据同一频率的对称S₀模式和反对称A₀模式这两种模式的Lamb波的传播速度不同，可用最大峰值法求出S₀和A₀这两种模式的Lamb波到达第此阵元的时间差Δt，最后由测得的这两种模式的Lamb波的传播速度c_S和c_A就求出冲击源到达此阵元的距离d_l，即：

$d_l=\frac{\mathrm{\Δ}{\mathrm{tc}}_S{c}_A}{c_S-c_A}---\left(15\right)$

为了提高距离估计的精度，将各阵元估计出的距离值d_l转化为相对于阵列参考点的距离r_l，取所有阵元测得的距离的平均值作为距离估计值R。

$R=\frac{1}{M}\underset{l=0}{\overset{M-1}{\mathrm{\Σ}}}{r}_l---\left(16\right)$

实施例

试件为T-30碳纤维复合材料层合板，尺寸为610×610×2mm，四边通过螺钉紧固在夹具上。碳纤维铺层的材料为T300/QY8911，铺层顺序为[45/0/-45/90/0/45/0/-45/0]_S。采用由7个相同的压电片组成的均匀线阵，相邻两个压电片中心间的距离为1.1cm。阵元的标号均是从左到右依次为0～6，均以0号压电片的中心为坐标原点，以压电阵列的轴线为实轴，以垂直于阵列的方向为法线方向，如图3所示。传感元件为PZT-5型压电片，直径为8mm，厚度为0.48mm。冲击工具为弹簧机械冲击杆，共有两档冲击能量，冲击头顶部为半球形，直径为10mm，如图3所示。由于碳纤维复合材料层合板为各向异性，各方向的波速不同，所以在进行DOA估计时波速取为0°、45°、90°及-45°方向上的测量值的均值；在进行距离估计时，S₀和A₀两种模式的波速分别取为测得的各阵元接收到的同一中心频率的S₀模式Lamb波的波速均值和A₀模式Lamb波的波速均值。

本实施例包括如下步骤：

步骤一：冲击响应信号。

当冲击锤头作用在试件上的(29.0°，47.6cm)位置时，每个压电元件上的冲击响应信号经过各自的电荷放大器，输入到数据采集卡上，当信号幅值超过设定的触发预置时，数据采集卡开始对冲击响应信号进行采集，并显示、存储在计算机中，0～6号压电片采集到的传感信号如图3所示，从上到下依次为0～6号压电片的信号波形。

步骤二：首先利用快速傅里叶变换对步骤一中获得的冲击响应信号进行频谱分析，然后用小波变换对由压电传感阵列直接测得的冲击响应信号做预处理，小波变换的中心频率为35kHz。然后从提取出的小波信号中截取信号幅值较大的0.4～0.6ms时间段内的数据进行DOA估计，截取获得的小波信号如图4所示。最后对截取得到的信号按照MUSIC算法的步骤进行DOA估计：

(1)根据公式(12)中得到阵列接收数据的协方差矩阵估计

(2)对

进行特征值分解；

(3)确定信号子空间和噪声子空间

(4)根据信号参数范围由式(14)进行谱峰搜索；

(5)找出极大值点对应的角度就是信号源的DOA。

按照上述步骤得到冲击源的方位谱如图5所示，最大谱峰对应的角度为28.7°，角度估计误差为0.3°，由误差可以看出DOA估计的精度非常高。

步骤三：冲击源距离的估计。按照技术路线步骤三中所述的冲击源距离估计方法进行距离估计。以对0号压电片接收到的信号分析为例，首先用小波变换提取出中心频率为25KHz的响应信号如图6所示。从图6可以看出提取出的中心频率为25KHz的响应有3个幅值较大的波包。在较低的频厚积范围内，由于S₀模式的Lamb波传播速度较快，所以先到达传感器，所以第一个波包为对称S₀模式波包，第二个波包为反对称A₀模式波包，第三个较大的波包为由反射造成的两种模式的混叠。然后用最大峰值法可求出S₀模式和A₀模式两个波包到达0号压电片的时间差，再根据测得的两种模式的速度，即可用(15)式估计出冲击源到0号压电片的距离d₀为46.4cm。按照相同的方法可估计出冲击源到达其它压电片的距离。然后运用公式(16)计算出的距离R为48.8cm，即为冲击源到达坐标参考点的距离，从而实现了冲击定位。

Claims (2)

1.一种基于多重信号分类算法的结构冲击载荷定位方法，其特征在于包括如下步骤：

步骤一：冲击响应信号的采集

（1）在所述结构上布置由7个相同型号的压电元件组成的均匀线阵；

（2）将2块4通道的数据采集卡0通道并联在一起，保证两块数据采集卡同步触发采集；

（3）设定采集参数，利用电荷放大器和数据采集卡采集均匀线阵输出的冲击响应信号；

步骤二：基于MUSIC算法进行冲击源的波达方向估计

首先用快速傅里叶变换对所述冲击响应信号进行频谱分析，然后根据选定的方向估计中心频率利用小波变换从所述冲击响应信号中提取出中心频率下的窄带信号，最后取窄带信号的幅值最大的信号运用多重信号分类（MUSIC）算法的进行DOA估计；

步骤三：基于Lamb波传播特性进行冲击源的距离估计

首先用小波变换从任一压电元件接收到的冲击响应信号中提取出距离估计中心频率下的Lamb波信号，获得同一频率的对称S ₀模式的Lamb波和反对称A ₀模式的Lamb波，用最大峰值法求出对称S ₀和反对称A ₀这两种模式的Lamb波到达所述压电元件的时间差，最后由测得的这两种模式的Lamb波的传播速度求出冲击源到达压电元件的距离，将各压电元件估计出的距离值转化为相对于均匀线阵参考点的距离，取所有压电元件测得的距离的平均值作为距离估计值。

2.根据权利要求1所述的基于多重信号分类算法的结构冲击载荷定位方法，其特征在于步骤一中，所述均匀线阵中压电元件间的距离小于冲击源信号波长。

Images