基于尾波干涉法反演结构微小损伤

赵元明，郭怀攀

(中国飞行试验研究院，陕西 西安 710000)

飞行试验中每次飞行前后机务人员必须目视检查飞机的进气道、 挂架、 舱盖等关键部位是否有损伤，但肉眼检查并不能发现某些极小的损伤和暗纹.尾波干涉法为检测这类损伤提供了可能性.多重散射波是振动信号或声音信号在散射介质中传播时，由于介质的非均匀性产生的散射现象.用脉冲信号对被测介质进行激励时，检测到的响应信号是脉冲直达波并且后面拖着长长的“尾巴”，这个“尾巴”称为尾波.与直达波相比，由于多重散射波在介质中来回多次传播，对介质的微小变化进行了重复的采样，可以把微小变化进行放大，因此对检测介质的微小变化具有更好的灵敏度，这就是尾波干涉原理.

Pacheco & Snieder在2005年研究得出多重散射波对于散射介质的微小变化的敏感性并非在空间中均匀分布，并根据振源和监测点的位置、 散射介质物理特性、 尾波的观察时间等变量推导出敏感核的解析式[1].本文对尾波干涉理论及其敏感核构造进行简要介绍，并基于MATLAB建立数值仿真平台，成功地反演出介质的微小变化的大小、 形状及位置分布，验证了算法及程序的有效性，为尾波干涉法在工程实践中的推广及应用奠定了基础.

1 尾波干涉

Snieder在总结前人工作基础上提出尾波干涉(Coda Wave Interferometry，CWI)原理[2]，尾波是振动信号或声音信号在散射介质中经过多次散射的结果，比直达波对介质有更多的采样，所以对介质的微小变化更为敏感，可以识别直达波所不能识别的微小变化.

(1)

式中：δt(t)为介质速度扰动前后同一监测点接收到波形的走时差，s1和s2分别为振(声)源和监测点的位置，K(s1,s2,x0,t)为敏感核，其表达式为

p(x0,s2,t-τ)dτ,

(2)

式中：p为两点之间的波场强度，是时间和位置的函数.在均匀的二维散射介质中，强度p的解析解有两种，一种是扩散方程的解

(3)

式中：l为平均自由程，r为振源到接收点的距离，c为波速，Θ为Heaviside函数.另一种解析解是散射方程的解

(4)

式中：D为扩散系数.把两种解析解带入敏感核进行对比，证实扩散方程的解更能准确地描述敏感核的分布[3], 因此在本文中，采用扩散方程的解构造敏感核.图 1 所示为当t=1 s和t=3 s时, 平均自由程均为500 m，台站间距均为3 km时敏感核的分布.从图 1 中可以看出敏感核呈马鞍形分布，分别在振(声)源和检测点处达到峰值，随着与这两点距离的变大，敏感核相应的降低.对比t=1 s和t=3 s时刻的敏感核分布可以看出，t越大，敏感核分布越广，说明从振(声)源发出的信号到达接收点之前经过的空间范围越广.

图 1 尾波观察时间为1 s和3 s时的敏感核分布Fig.1 Distribution of sensitive nucleus when the observation time of tail wave are 1 s and 3 s

m=m0+CmGT(GCmGT+Cd)-1(d-Gm0),

(5)

式中：Cd为测得的数据的对角协方差矩阵；GT为G的转置，G=Δvk(r′,r,s,t)；m0为初始模型，为零矩阵；Cm为平滑矩阵[4],

式中：λ为相关长度，λ0为网格间距，δm为先验标准差.

2 数值仿真平台

数值仿真平台主要由以下几部分组成： 激励信号、 速度场、 信号传播方式、 监测点的分布等，以此平台为基础，输入激励信号获得响应信号，计算获得速度扰动的空间分布.数值仿真平台结构图如图 2 所示.

图 2 数值仿真平台结构图Fig.2 Structure diagram of numerical simulation platform

2.1 激励源和扰动区域的设定

激励源为一个脉冲信号，如图 3 所示，信号持续时长为6π/100，激励信号函数为

(6)

如图 4 所示，在速度场的正中间施加一个速度+0.5%、 2 km×2 km的矩形速度扰动区域，根据监测点的布设位置，在振(声)源处发出单位脉冲激励，11个监测点连续记录波形.

图 3 激励源波形Fig.3 Excitation source waveform

图 4 在10 km×10 km的监测区域正中央施加一个矩形2 km×2 km的速度扰动区域Fig.4 Add a rectangle with 2 km×2 km velocity perturbation region to the central of 10 km×10 km monitoring area

2.2 速度场和监测点的布设

实际的传播介质是非均匀性的，因此在速度场中，用速度的微小波动模拟介质的非均匀性.为了取得理想的仿真结果，假设速度场是无限大的，速度场的边界应该设置为吸收边界，吸收边界采用Clayton Engquist majda吸收边界条件进行仿真，反射率大约为2.5%，满足本研究的要求[5].

首先，在10 km×10 km的二维速度场中，边界为吸收边界，运用传播速度不同模拟介质的非均匀性，速度场的速度模型为200×200的矩阵，速度矩阵的均值为5 km/s, 方差为200，最小值为4 958.7 km/s，最大值为5 083.5 km/s，如图 5 所示.

监测点的布设原则是确保均匀地覆盖监测区域，通过尝试不同的布设方案，发现监测点布设位置对反演结果有一定的影响，这不在本文讨论范围，因此在本文中监测点的位置固定不变，以消除监测点位置不同对结果的影响，监测测点的分布如图 6 所示，1为振(声)源，2～12为监测点.

图 5 二维速度场模型Fig.5 Two-dimensional velocity field model

图 6 传感器布设分布Fig.6 Sensor layout distribution

2.3 信号传播方式

运用波动方程进行波的传播计算.

(7)

式中：f(x,z)为振源，v为波的传播速度，p为位移.

把等式两边的偏导数用泰勒展开式展开，整理后波动方程转化为差分格式，从泰勒展开式的数学理论得知，展开式的阶数越高越逼近于真实值，在权衡运算效率和精确度后，选择了时间二阶、 空间八阶的有限差分格式.

pk+1(i,j)=2pk(i,j)-pk-1(i,j)+

s(t)*δ(i-i0)*δ(i-j0),

(8)

式中：uk(i,j)为在第k次(对应的时间)迭代时在(i,j)处的位移；v为振动的传播速度，是位置(i,j)的函数； Δh为网格间距，监测区域内网格为均匀分布，因此Δh为固定常数； Δt为采样时间间隔，为采样频率f的倒数.

3 仿真结果分析

图 7(a) 为在2号监测点记录到的波形，最大波峰处为直达波，波峰后拖着长长的尾巴即为散射造成的尾波.由于尾波干涉法观察的是介质的微小变化，因此同一个监测点处扰动前后的波形基本一致，只考虑相位差，忽略幅度的变化[6].

从图 7(a) 中可以看到扰动前后的波形几乎完全重合，这是因为介质发生的是微小变化，所以对波形的幅度及相位的影响非常小.当把波形放大观察时，如图 7(b) 中1～1.4 s的对比波形，仍然发现扰动前后波形高度重合； 而3.1～3.5 s时的波形，发现有很明显的走时差，这是因为尾波把介质的微小变化进行了放大.

图 7 扰动前后的波形及部分波形放大图Fig.7 Waveform before and after disturbance recorded by monitoring point and partial waveform enlargement

3.1 走时差和敏感核的计算

计算走时差常用的方法有两种： 延展法和互相关法[7].延展法主要针对实验数据信噪比较低的情况，通过对波形进行拉伸(压缩)后与扰动前采集到的波形进行互相关计算获得走时差，有较好的稳定性； 互相关法是比较常用的一种方法，相对于延展法计算量较小，对于本文具有较高噪比的仿真数据，采用互相关法[8].选择尾波时间段为3～7.7 s的数据计算走时差，每段数据长度为0.3 s, 重叠0.1 s, 每个检测点可以得到15段尾波数据，用互相关法计算得到的走时差(以振源和2号监测点为例)如图 8 所示[9，10].

图 8 走时差Fig.8 Time lag

以1号点位置为激励源，其它11个监测点记录响应信号，在扰动前后均可记录到11个振动信号，对监测点扰动前后波形进行互相关计算可得到总计11×15=165个走时差.

敏感核表示从振源发出的信号在某一时刻经过空间某一点到达接收点的概率分布.使用散射方程的解计算敏感核.

k(r′,r,s,t)=

(9)

把10 km×10 km的观测区域划分为20×20的网格，网格间隔为500 m.得到振源和2号监测点，时间为3.25 s时刻的敏感核如图 9 分布，从图中可以看出，敏感核为马鞍形，在震源和监测点处敏感核数值较大.

图 9 t=3.25 s时振源和2号监测点的敏感核分布Fig.9 The distribution of sensitive nucleus of vibration source and No.2 monitoring point in the tail wave at t=3.25 s

3.2 速度扰动分布

图 4 所示在监测区域的正中央施加矩形 +0.5% 的速度扰动区域，图 10 为反演结果，黑色方框标注出施加的扰动区域，颜色条表示反演得到速度扰动得百分比，从图中看出反演结果与施加的扰动区域吻合度较好，不仅成功地定位了速度扰动区域的位置，并且反演出扰动的幅度.

这种介质内的振动(声音)传播速度的微小变化可能由于各种原因导致，例如温度、 结构裂纹、 密度变化等，需根据具体介质进行分析.同时也进行了其它区域的反演计算，图 11 是在监测区域的左上方施加+0.5%速度扰动、 右下方施加-0.5% 速度扰动的反演结果.

图 10 正中央施加+0.5%的速度扰动的反演结果Fig.10 Inversion results of velocity perturbation with +0.5% in the center

图 11 施加+0.5%速度扰动和-0.5%速度扰动的反演结果Fig.11 Inversion results of +0.5% velocity perturbation and -0.5% velocity perturbation

4 结束语

本文对尾波干涉理论进行了简要介绍，并基于MATLAB建立数值仿真平台，证明了该算法能够准确监测到结构微小的损伤及定位损伤的具体位置，更深层次的研究以及工程应用需进行继续探究，尾波干涉法是一种较新的探伤技术，尤其是对介质的微小变化具有较高的敏感性，其在飞机试飞、 桥梁工程、 石油开采等方面有广泛的应用.