一种基于克里金法的数据驱动方法，用于确定复合材料中同时发生的损伤位置

《MECHANICAL SYSTEMS AND SIGNAL PROCESSING》：A kriging-based data-driven approach for localization of simultaneous impacts in composites

　　•基于人工神经网络的单次及两次同时冲击在复合材料中的定位方法。•同时冲击产生的波场叠加会带来到达时间信息的不确定性，从而限制基于机器学习的定位精度。•一种基于奇异值分解的到达时间输入特征集，用于训练用于冲击定位的人工神经网络。•即便训练数据有限，也可通过四传感器配置实现对两次同时冲击的定位。•基于克里金的方法可提升实验场景下的预测精度。

引言
由于高强度重量比和设计灵活性，复合材料在现代飞机中的应用日益广泛。然而，这类结构极易受到低速冲击的影响，比如技术员手中工具滑落、飞行中遭遇鸟击，或是恶劣天气下的冰雹袭击。高速冲击通常会造成明显的损伤，而低速冲击则往往仅留下几乎难以察觉的冲击痕迹[1]。长期而言，即便是一些细微且未被发现的损伤，也会缩短飞机的使用寿命并增加维护成本[2]。

为检测复合材料中的低速冲击，人们广泛采用了多种无损检测技术，如X射线断层扫描、热成像、射线检测、超声波检测以及声发射技术[3]。其中，声发射技术无需任何驱动装置即可实现实时冲击检测[4]。

在基于声发射的定位方法中，三角测量法最为简单。该方法通过测量由声源产生的信号被安装在结构上的传感器网络捕获时的到达时间差，再结合已知的波速分布，在假设波速均匀且与方向无关的前提下，通过求解线性方程组来估算冲击位置[5]。但对于复合材料这类各向异性且材质不均匀的材料，由于群速会随传播方向变化，这一假设就不再适用。为解决这一问题，Ciampa等人[6]提出了一种使用六个换能器阵列的方法，通过双曲定位技术求解非线性方程组，从而提高在复合材料中定位声发射源的精度。

在现实场景中，飞机结构常常会同时遭受多次冲击，这意味着可能有多个声发射源同时活跃。这给上述传统定位方法带来了巨大挑战。已有研究多集中于在铝板、钢板等各向同性材料中定位两个同时存在的冲击源。Sai等人[7]在铝合金板上安装了四个光纤布拉格光栅传感器来捕捉冲击信号，他们利用香农小波变换提取窄带成分，并通过时间反演聚焦技术来定位两个冲击源。但该方法依赖于精确的群速估算和小波提取，因此较为复杂且耗时。He等人[8]提出了一种将奇异值分解与波束成形相结合的混合方法，即SVD-BFM方法，用于定位钢板中的双重冲击源。Tai等人[9]则提出了一种频域陷波加权波束成形算法，通过将陷波理论与频域加权波束成形相结合，来定位钢板中的两个同时冲击源。该算法假设传到传感器上的信号幅度相同，仅有时间差，这一假设在传感器间距较近的各向同性材料中较为成立。但在复合材料这类各向异性材料中，不同传播路径上的波色散和衰减会导致幅度变化和时间偏移，使得该假设不再适用[10]。Zong等人[11]展示了如何使用带有七个传感器的近场二维音乐算法对复合板上的两次同时冲击进行定位。Ren等人[12]则开发了一种基于空间波数过滤技术的复合结构多冲击成像与定位方法，该框架使用了十三个传感器，成功实现了对监测区域内多次冲击事件的精准定位。

对于这类逆向定位问题，人们通常会采用机器学习方法。其原理是从已知冲击位置的信号中学习规律，进而预测未知冲击事件的位置。Liu等人[13]结合小波包能量、多域特征以及蝙蝠算法支持向量回归模型，用于定位碳纤维增强塑料板中的单次冲击。Damm等人[14]则使用卷积神经网络，以传感器信号的短时傅里叶变换作为输入，来确定单次冲击的位置。Wang等人[15]通过时间序列特征、递归量化特征以及递归绘图特征来训练卷积神经网络模型，用于定位复合飞机结构上的单次冲击。Sun等人[16]提出采用贝叶斯更新方法，构建最优的人工神经网络模型用于冲击数据集的训练。Deng等人[17][18]则在贝叶斯反演框架内，结合虚拟时间反演、到达时间以及信号包络特征来实现冲击位置的定位，他们采用马尔可夫链蒙特卡洛采样方法来估计与冲击位置相关的参数后验分布，随后再通过加权平均技术对多频率后的验分布进行决策级融合。Khodaei等人[19]比较了最大信号幅度、信号包络以及到达时间三种特征，认为到达时间是定位单次冲击最关键的信号特征。Balasubramanian等人[20]系统地评估了包括人工神经网络、卷积神经网络以及长短期记忆网络在内的多种神经网络架构，研究发现，当使用到达时间特征对人工神经网络模型进行训练时，其在预测单次冲击位置方面的准确率最高。Kossel等人[21]发现，当有两个声发射源同时存在时，相互重叠的信号会使得到达时间无法用于精准定位，他们在铝制试样上通过独立成分分析实现盲源分离，以便在定位前分离信号，但该方法仍缺乏有效的信号分组策略，需进一步研究。

当多个冲击同时发生时，相互重叠的信号会掩盖后发生的冲击源的到达时间信息，这使得像人工神经网络这样的元模型难以训练[22]。解决这一问题要么可以通过增加传感器数量来实现，要么可以通过分析足够长的信号中的到达与反射模式来完成。这一局限性表明，仅依靠到达时间特征是不够的，要实现多次冲击的精准定位，还需要更丰富的信号特征。

在数据驱动型问题中，收集足够的训练数据至关重要。而生成多次冲击的训练数据既耗时又需要多种不同的实验组合。由于航空航天领域的训练数据极为匮乏，因此元模型应当能够适应各种冲击场景下多次同时冲击的定位需求。对于波传播特性不均匀的复合结构，有限元建模可用于生成训练数据[23][24][25][26]，不过要在所有环境及运行条件下进行完整的兰姆波物理仿真，成本极高且耗时极长。

现有的多源分离技术，如波束成形、盲源分离、多信号分类等方法，都基于波传播各向同性的假设，通常还需要精确的波速模型以及密集的传感器布局，因此并不适用于复合材料中多次冲击的检测。为应对这些挑战，本研究提出了通过奇异值分解进行特征工程，并结合基于克里金的后处理方法。本研究的主要贡献如下：
i) 提出了一种新的特征表示方法，利用奇异值分解处理到达时间信息，从多次冲击产生的重叠信号中提取可靠的到达时间特征，从而克服传统到达时间特征在多次冲击场景下的局限性。与传统将奇异值分解用于去噪或明确源分离不同，该方法将分解得到的分量视为互补的到达时间假设，而非真实的声源。
ii) 构建了一个统一的框架，该框架可使用相同的传感器配置和特征集来检测单次及多次冲击，从而无需增加传感器密度，也无需事先进行信号分离。
iii) 训练过程仅依赖基线的有限元仿真数据，且通过包含不同冲击力及冲击时间差的未见仿真试验与实验进行了验证，这证明了该模型即便仅在数值仿真数据上训练，也能在真实且不确定的运行条件下保持稳定性，因此无需花费大量成本构建实验训练数据集。
iv) 为进一步提升上述模型的稳定性，本研究引入了基于克里金的后处理方法，该方法无需重新训练或修改信号特征，即可提高实验条件下的定位精度。

本文的其余部分结构如下：第2节描述了实验冲击装置，第3节详细介绍了基于奇异值分解和克里金的定位方法的理论基础。第4节阐述了实验装置的数值建模方法，第5节说明了如何利用单次冲击信号生成两次同时冲击信号。第6节讨论了利用奇异值分解进行特征提取的方法，第7节则明确了使用仿真数据训练人工神经网络时的输入特征与输出目标。第8节展示了研究结果与讨论，包括使用实验信号对所提方法的验证。最后，第9节对整个研究进行了总结。

材料与试样
本研究在堆叠顺序为[0/45/90/-45]s的复合层压板上进行了低速冲击实验。该层压板由环氧树脂和E玻璃纤维制成，纤维体积分数约为50%，是通过图1(a)所示的真空辅助树脂灌注工艺制造而成的。随后按照树脂灌注系统的要求进行了后固化处理。所有制造工作均在环境温度为25℃的受控条件下进行。

奇异值分解
假设复合材料表面同时发生了位于第i位和第j位的两次冲击。传感器采集到的信号可以用矩阵形式表示，记为X∈R(M×N)，其中M代表传感器的数量（或测量通道数），N代表时间样本的数量（且M>N）。应用奇异值分解后，X可以表示为公式(1)：
X=UΣVT
其中，U∈RM×M和V∈RN×N是分别包含左奇异向量和右奇异向量的正交矩阵，Σ∈RM×N则是対角矩阵。

数值建模
冲击仿真是在ABAQUS显式动力学软件中进行的，模拟对象为带有表面粘贴的压电传感器的层压板，这些传感器的尺寸和材料属性与2.2.1节中描述的相同，如图7所示。每层复合材料都被单独建模，其材料属性则根据各自的局部坐标系进行设定。GFRP和压电材料的力学性能见表1[32]。为减少计算时间，复合材料层被...

单次冲击下的传感器响应
对于每次冲击事件，都会从四个压电传感器s1、s2、s3和s4中记录长度为N的兰姆波信号。与第i位冲击位置相关的传感器响应si可以用公式(9)表示：
si=s1is2is3is4i
冲击作用点位于结构上预定义的35个网格点上，每个冲击网格点对应的记录信号可以用矩阵形式表示，如公式(10)所示：
Si35×4=s1s2::s35=s11s21s31s41s12s22s32s42::::::::s135s235s335s435
图9展示了...

利用归一化阈值法从冲击信号中提取到达时间
假设如图11所示，有四个传感器（S1、S2、S3和S4）捕获到了冲击引发的信号。信号的到达时间可以被定义为信号超过特定阈值的时间，该阈值取决于传感器的灵敏度、最小可检测应变以及系统噪声水平[35]。在一些研究中，人们会通过阈值处理来抑制干扰精确到达时间测定的噪声。初始阈值设置为最大值的1.4倍。

使用仿真数据训练人工神经网络的输入特征与输出
为确定复合材料表面坐标分别为(xi,yi)和(xj,yj)的两个同时产生的冲击源的位置，本研究采用了人工神经网络。输入向量的结构为[Z]1×36=[ΔT1×4T1×32]，这样每个训练样本就有36个输入特征。用于训练的数据集包含595组双重冲击样本和35组单次冲击样本，每组样本都对应一个输出向量，其中包含了两个冲击源的坐标，表示为x...

算法1. 信号预处理
输入：sa(n)，其中a=1:4，n=1:N
输出：ya(n)
流程如下：
对于f=1到N：
计算差值：d(n)=sa(n)-sa(n-1)，且sa(0)=0
进行零填充：d(1)=0
进行积分：z(n)=∑k=2nd(k)+d(k-1)2，且z(1)=0
进行归一化：y(n)=∣z(n)∣/∣z(n)∣max
更新值：sa(n)=y(n)
结束循环
返回结果：ya(n)

使用实验信号与仿真信号验证所提方法
为验证所提出的方法，仿真信号必须与实验信号具有高度相关性。图25对比了在位置i=1、2、3和4处发生的单次冲击时，传感器S1、S2、S3和S4所采集的实验信号与其对应的仿真信号。可以看出，球体冲击产生的数值信号在波形特征上与实验记录的信号十分相似。

结论
本研究提出了一种新的方法，用于定位复合层压板中的多次低速同时冲击，该方法针对波场叠加效应以及各向异性波传播对传统到达时间特征的削弱问题进行了改进。这里的主要创新点在于利用SVD来增强TDoA特征。此处使用SVD并非用于传统的去噪或信号源分离，而是为了在稀疏感知条件下，通过数据驱动的方式区分重叠的冲击信号。研究采用了通用的ANN框架。CRediT作者贡献声明：V. Mishra：写作——审阅与编辑、写作——初稿撰写、可视化、验证、软件、方法论、研究分析、正式分析、数据整理、概念构建；S. Kumar：可视化；M.R. Sunny：写作——审阅与编辑、可视化、监督、软件、资源、方法论；A. Mukherjee：监督。利益冲突声明：作者们声明自己没有已知的可能会影响本文所述工作的财务利益或个人关系。V. Mishra | S. Kumar | M.R. Sunny | A. Mukherjee