基于稀疏重构的尾波干涉成像方法*

张涛 侯宏† 鲍明

1) (西北工业大学航海学院, 海洋声学信息感知工业和信息化部重点实验室, 西安 710072)

2) (中国科学院噪声与振动重点实验室(声学研究所), 北京 100190)

尾波干涉成像是利用尾波时延和扩散近似敏感核来反演散射介质中微小速度扰动空间分布的技术.该问题本质上是一个欠定问题, 一般无确定解, 导致其难以精确定位介质中微小波速变化的区域.为解决上述缺陷, 本文利用速度扰动分布在空间上具有稀疏性的特点, 提出了一种基于压缩感知理论中稀疏重构算法的尾波干涉成像方法.该方法可在散射介质中较准确地获取速度扰动的空间位置和范围, 同时具有较高的计算效率.数值仿真和实验结果表明:相比于现有的线性最小二乘差分成像方法, 无论是单个还是多个扰动区域,该方法均能更精确地进行定位成像, 同时明显减少了计算时间.

1 引 言

尾波干涉 (coda wave interferometry, CWI),是一种利用介质中的多次散射波检测介质微小变化的测量方法.多次散射波是声信号或振动信号在散射介质中传播时, 由于介质的非均匀性(存在如颗粒、裂缝、包体以及孔洞等非均质体)而发生多次散射、反射产生的, 在波形图中显示为直达波后面拖着的长长的“尾巴”, 因此也称尾波.尾波由于在介质中来回传播, 对介质更密集地采样, 使得介质中的微小变化不断迭加放大, 从而对介质性质的微弱改变具有极高的敏感性, 检测出直达波无法识别的微小变化.自2002年CWI法首次由Snieder等[1−3]提出后, 便在地球物理学和材料科学等领域得到了广泛的研究与应用[4−9].文献[7]在SAFOD(San Andreas Fault Observatory at Depth) 进行了一次主动震源井间实验, 观测到距试验点3 km外的一个3级地震发生前数小时, 井内波速增加了0.8%.文献[8]运用CWI测得水泥样品水饱和进程和楼房承载柱的尾波波速随环境湿度的变化,水饱和使得尾波波速相对变化0.7%, 楼房承载柱随相对湿度每变化1%, 尾波波速变化0.213%.文献[9]基于CWI的测量方法, 利用超声尾波在实验室观测了1.5 m岩石断层的黏滑过程, 获得了精度高达 10–6的相对波速变化, 相当于~10 kPa 的应力变化.以上研究都取得了一定成果, 均通过CWI的测量方法获取了介质的相对波速变化量, 但是这些变化量仅仅是观测介质各位置相对波速变化的加权平均, 并不能反映波速变化的空间分布情况.然而, 正确定位波速变化区域对于分析引起波速变化的物理机制十分重要, 因此对波速变化的空间分布成像近年来成为了CWI法的研究前沿与难点之一.早在 2005年, Pacheco 和 Snieder[10]便 对CWI测量方法进行了延伸, 利用多次散射波场的扩散近似, 提出了一种通过扩散近似敏感核将平均走时延迟与介质波速局部变化联系起来的技术, 使定位局部微小变化成为了可能.文献[11]分析了2010年6月至 12月La Reunion 岛上 Piton de la Fournaise活火山的连续环境噪声记录, 利用10月的大爆发前观察到的速度变化确定了爆发点速度显著下降(高达0.6%)的区域.文献[12]使用环境噪声互相关和延展法处理了墨西哥Volcán de Colima地区超过15年的连续记录地震数据, 最明显的速度变化是2003年Tecomán附近的7.4级地震期间下降了2.6%, 并采用了基于辐射传输近似的敏感核反演定位水平面上的扰动.文献[13]描述了一种LOCADIFF技术, 利用CWI法以及结合扩散传播模型的最大似然法反演定位非均匀强散射介质中的微小变化, 并通过二维有限差分进行了数值模拟.文献[14]将LOCADIFF这种超声波成像技术成功应用在混凝土结构无损检测与评价中.对于CWI结合敏感核成像介质波速变化的空间分布以定位局部微小改变的研究, 国内外目前涉及较少, 反演问题是这种成像技术发展与应用的阻力之一.尾波干涉成像技术归根到底是对一个矩阵形式的反问题进行求解, 该问题是一个欠定问题, 有无穷多解, 同时具有病态性, 难以精确求解.前述研究中均采用地球物理反演中常用的线性最小二乘差分反演方法, 该方法存在求解精度不高、参数选取困难的缺陷, 并且难以定位多个扰动点, 在实际工程中的应用十分困难.

另一方面, 由 Donoho[15]以 及 Candes 和Romberg[16]于2006年正式提出的压缩感知理论(compressed sensing, CS)近年来逐渐兴起并趋于成熟, 在雷达探测、图像处理、声呐定位等领域都有了广泛的应用, 在物理成像方面也表现出强大的生命力[17−21].文献 [17]采用分块压缩感知思想, 提出一种基于lp范数的压缩感知图像重建算法, 提高重建图像质量的同时缩短了成像时间.文献[18]提出了基于稀疏测量的超分辨压缩感知鬼成像重建模型, 并搭建了实验装置, 验证了该模型的有效性与优越性, 推动了鬼成像的实用进展.文献[19]利用声源的空间稀疏性, 在压缩感知理论下建立了矢量阵聚焦定位新方法, 克服了相干声源分辨困难、贡献评价不准确等问题, 且定位精度高, 背景抑制能力强.CS理论指出, 若某信号是可压缩的或经某种变换后稀疏, 则可通过随机观测矩阵(与变换基不相关)将高维信号投影至低维空间中, 最终通过求解优化问题从少量低维数据高概率地重构原始高维信号.而在尾波干涉成像问题中, 实际扰动通常是空间上聚焦的, 扰动的空间区域相比于监测的整体介质区域往往是稀疏的, 恰好满足CS理论对于信号稀疏性的要求.本文基于CS理论, 提出了一种尾波干涉成像的新的解算方法.该方法首先通过CWI和扩散波敏感核构建了用于成像的反演模型; 其次引入CS理论框架下的稀疏变换将欠定方程重新构造; 最终通过压缩感知重构算法——正交匹配追踪算法 (orthogonal matching pursuit,OMP)[22]求解该欠定方程, 并获得速度扰动的空间成像.此方法与先前的线性最小二乘差分成像方法相比, 在单扰动和多扰动情况下对于扰动区域位置和范围的确定都更加准确, 且执行简易, 没有复杂的参数选取, 同时还具有更高的计算效率.最后结合具体仿真算例, 给出了与线性最小二乘差分成像方法的比较结果, 证明了该成像方法的可行性、有效性和高效率.

2 尾波干涉成像基本原理

2.1 尾波干涉基本原理

尾波干涉理论指出, 多散射介质中的微小改变将导致波速的微弱变化, 主要表现为尾波信号在扰动前后的相位偏移.图1显示了微小结构变化发生前后的多散射介质中某处记录的波形信号, 插图显示了直达波(图1(a))和尾波(图1(b))波形的详细信息.不难看出, 直达波于扰动前后在幅值和相位上几乎没有变化, 而尾波则表现出明显的走时延迟, 该时延即为尾波对于介质中微弱变化多次采样后的体现, 是直达波所不能包含的物理信息.

尾波干涉理论的实质是通过测量两个不同时刻尾波列的相位差来获取介质在该时段的波速变化[1−3].互相关法和延展法是目前基于此理论来提取介质波速变化的主要方法[23,24], 考虑到仿真信号信噪比较高, 本文采用计算效率更高的互相关法[23].步骤如下:假设已获取介质变化前后的两列尾波信号uunp和uper, 通过两列波的互相关计算获得波速变化, 即

图1 多散射介质中扰动前后波形的比较 (a) 直达波扰动前后的波形; (b)尾波扰动前后的波形Fig.1.Comparison between typical time traces of a wave propagating in a multiple scattering medium before and after a small perturbation:(a) The first arrival waves before and after a small perturbation; (b) the coda waves before and after a small perturbation.

式中, 2T表示要分析的尾波部分的时间窗长度;t为时间窗中心位置;ts表示互相关中所用的走时差;R的值表示两列波的相似程度.

当ts的某取值tsmax使得R(ts) 取最大值时, 该走时差tsmax即为所分析尾波部分扰动前后的时间延迟δt.那么根据Snieder[1]对尾波干涉的系统阐述, 可得两列信号的相对波速变化

式中 δυ表示波速扰动, 是整体介质波速变化的加权平均;υ为扰动前的介质波速.

2.2 敏感核

为了获取局部扰动的空间分布情况, 需要将上述尾波干涉求得的总体平均时延δt与介质波速局部变化联系起来.然而, 由于介质的复杂性, 多散射效应产生的尾波的传播路径长而无序, 因此难以将每个时延与一组特定的传播轨迹配对, 确定变化的具体来源也就十分困难.一些学者[10,13]并没有尝试确定每个尾波序列的传播路径, 而是研究了一种基于时空传播统计描述的替代方法, 引入敏感核,K(s1,s2,x0,t) 来描述波列在位置s1发射, 途经位置x0, 并在总传播时间t后在位置s2被接收的概率大小, 此概率描述了尾波在位置x0处传播时间的空间密度.

式中的p(s1,s2,t) 表示波列从s1经过时间t到达s2的概率.在实际应用中, 此概率可以用波场强度等效, 波场强度是位置和时间的函数, 可以通过扩散方程的解来近似.在无限大二维介质中, 扩散方程的解为

式中, |s1−s2| 表示s1和s2之间的距离;D为介质的散射系数, 与其物理性质有关.通过激励源与接收点的位置和尾波的传播时间, 对介质空间逐点计算K(s1,s2,x0,t) , 即可得到敏感核的空间分布.图2所示为尾波观察时间t分别为1 s和5 s, 散射系数为 5.8×105m2/s , 激励源 (S)与接收点 (R)间距5000 m时的敏感核空间分布图.从图2中可以看出, 敏感核在三维空间上呈马鞍形分布, 在激励点和接收点处达到峰值, 敏感核的值随着与这两点的距离增大而降低, 敏感核的空间分布也随着尾波观察时间的增加而展宽.t= 1 s时, 以直达波和单散射波为主, 对介质的采样密度和范围较小,t= 5 s时, 以多次散射产生的尾波为主, 对介质进行密集采样的同时增加了探测范围.

2.3 反演成像模型

基于扩散近似敏感核, 文献[10,11,13]提出了一个线性模型, 将尾波干涉法获取的平均时延数据通过敏感核项K(s1,s2,x0,t) 与相应的局部变化联系起来.

图2 基于扩散近似的二维敏感核示例 (a) t= 1 s时的敏感核空间分布; (b) t= 5 s 时的敏感核空间分布; (c) t= 1 s 时的敏感核俯视图; (d) t= 5 s时的敏感核俯视图Fig.2.Examples of sensitivity kernel based on the diffusion approximation in 2-D:(a) Spatial representation of the sensitivity kernel when t= 1 s; (b) spatial representation of the sensitivity kernel when t= 5 s; (c) vertical view of the sensitivity kernel when t=1 s; (d) vertical view of the sensitivity kernel when t= 5 s.

为了使探测范围尽量覆盖整个介质, 尾波干涉成像需要大量的敏感核, 因此会有许多形如(5)式的方程, 为清晰表达, 可将线性方程组写成矩阵形式:

式中,T∈RM,1表示通过不同的激励接收对和不同的尾波观察窗口获取的M个时延数据;V∈RN,1表示空间各网格点对应的相对波速扰动值,N为网格总数;G∈RM,N表示T到V的映射,其矩阵元素Gij=∆sKij, 其中Kij为敏感核矩阵K的元素,K∈RM,N的每行代表一个敏感核, 每列代表一个网格, 每个网格的面积为 ∆s.

利用线性最小二乘差分法对V值进行反演, 即

式中,V0为先验初始值, 一般为零矩阵;GT表示G的转置;Cd为测得数据的对角协方差矩阵;Cm为介质空间元素平滑矩阵, 可用下式计算:

式中,σm为先验标准差;λ0为网格间距;λ为相关长度.σm和λ一般通过 L 曲线法[25]选取.

在反演计算时, 还需要循环迭代Cm以寻找最优的V值, 迭代公式如下:

3 基于稀疏重构的尾波干涉成像方法

压缩感知理论指出, 假定N维真实信号x∈RN在某变换域下是L稀疏的(即则可利用信号x的任意M(M≈Lln(N/L)) 个线性测量值高概率精确重构x.压缩感知原理主要包括三个过程:信号的稀疏表示、测量矩阵构建以及算法重构.

首先, 对信号进行稀疏表示如下:

式中,ψ∈RN为稀疏变换矩阵;a∈RN,1为x的稀疏表示.

其次, 构建一个与变换基底ψ不相关的观测矩阵对a进行线性观测

式中,ϕ∈RM,N为观测矩阵;y∈RM,1是a在观测矩阵下的观测值;Θ=ϕψ称为传感矩阵.

最后, 对信号进行重构, 是压缩感知理论中极为关键的一环.(11) 式中y=Θx称为欠定方程,y中的未知量有N个, 而方程仅有M个, 理论上该反问题有无穷多个解, 但如果Θ满足有限等距性质 (restricted isometry property, RIP), 则有唯一最稀疏解.该过程等效为求解一个最优化问题

求解 (12)式是最小l0范数优化问题, 属于NP-hard问题, 直接求解此问题数值极不稳定, 且计算复杂度高, 实际工程中实施困难.Candes[26]提出,在一定条件下, 可以利用l1范数代替l0范数, 即

因此, (12)式的问题转化为(13)式的凸优化问题,可以通过线性规划理论求解, 本文采用应用广泛的OMP法作为重构算法.

可以观察到(6)式和(11)式在数学形式上相似, 且扰动的空间区域相比于整体介质区域是稀疏的, 恰好满足压缩感知理论对于信号稀疏性的要求, 因此, 如果将G,T和V分别视为观测矩阵、测量值矩阵和真实信号, 则可在压缩感知理论框架下, 通过信号稀疏表示和稀疏信号重构的方法求解(6)式.

本文引入稀疏变换矩阵P对V进行离散稀疏变换.从而有

式中,GCS=GP−1,VCS=PV.常用的稀疏变换矩阵有离散傅里叶变换矩阵、离散余弦变换矩阵、离散小波变换矩阵等, 图像信号在二维小波变换域是稀疏的, 声信号在傅里叶变换域是稀疏的, 本文本质上是对一维声信号的处理, 因此采用离散傅里叶变换矩阵.利用OMP等凸优化算法求解出方程(14)的稀疏解VCS, 最终通过稀疏逆变换即可求得原始解V, 即

整个过程无须使用P而仅须用到其逆矩阵P−1,对于傅里叶变换矩阵,P−1=PT.

4 数值仿真与结果

为了验证基于稀疏重构算法的成像方法, 现通过数值仿真进行分析.本数值仿真使用的激励源为 15 Hz 主频的雷克子波, 在 1 0km×10km 的非均匀散射介质中心处激发, 持续时间为 6π/100s ,时间采样间隔 ∆t为0.001 s.通过传播速度的不同模拟介质的非均匀散射特性, 将介质的速度场划分为 2 00×200 网格, 产生均值为 5 000m/s , 标准差为500m/s的随机数矩阵作为速度矩阵, 其中最大速度为 7 320.95m/s , 最小速度为 3 045.29m/s , 如图3所示.接收点的布设原则是确保均匀覆盖介质区域, 本仿真中采用的36个接收点和1个激励源的位置分布如图4所示.

本数值仿真运用二维声波方程模拟波的传播,

式中,p(x,y,t) 为质点的位移函数,f(x,y,t) 为激励源函数,v(x,y) 为介质在 (x,y) 处的速度.采用时间二阶、空间八阶的有限差分法对(16)式进行数值离散, 即有

图3 二维速度场模型Fig.3.2-D velocity field model.

图4 激励源及接收点布设Fig.4.Layout of the source and receivers.

式中,nx,ny分别为划分在x,y方向的网格个数;nt为时间间隔个数; ∆h,∆t分别为空间、时间网格上的划分步长, ∆h=50m , ∆t=0.001s ;pk(i,j)表示第k次(对应的时间)迭代时在 (i,j) 处的位移.

为了模拟相对无限大的传播空间, 本文仿真的速度场边界设置为吸收边界, 采用Clayton_Engquist_majda吸收边界算法, 反射率约为2.5%, 满足本研究的要求[27].

下面基于以上仿真环境, 根据第2和第3节的原理, 给出针对五个速度扰动算例的线性最小二乘差分成像方法和本文成像方法的性能比较.其中,速度扰动的空间分布通过将介质剖分为 2 0×20 的网格来离散表达, 则有400个待反演的未知量, 网格单元边长为500 m; 选择尾波观察时段为1.5—4.7 s, 每段窗口长度为 0.5 s, 重叠 0.2 s, 每个接收点可以获取扰动前后各10段尾波数据, 36对收发对即可计算得到360个时延数据; 敏感核的构造中散射系数D=8×104m2/s ; 线性最小二乘差分方法中的相对长度λ=750m.

对于单点扰动, 通过算例1和算例2进行了两种算法的成像对比, 如图5和图6所示.其中, 算例1在速度场的左下方(具体位置见图5蓝色方框)施加一个大小为+0.5%, 面积为 1 km×2km 的速度扰动区域; 算例2在图6右上方蓝色方框区域添加一个大小为+0.5%, 面积为 2 km×2km 的速度波动.

图5(a)和图6(a)分别为算例1和算例2利用线性最小二乘差分方法迭代10次的结果, 通过L曲线法获取的σm分别为 0.00328和 0.00266,图5(b)和图6(b)分别为算例1和算例2在压缩感知理论框架下利用OMP算法迭代5次和7次的结果.从算例1的结果可以看出, 两种算法均能准确地突显出扰动的位置, 而在扰动范围的确定方面本文方法更为接近真实值; 但从图6中可以清晰地看出, 线性最小二乘方法反演出的算例2扰动区域与真实设置相距甚远, 而本文方法依然能得到较好的反演结果.可见, 本文方法在单扰动情况下较线性最小二乘法有更高的精确性与稳定性.

图7展示了实测数据的处理结果.该实验是一个小型的模拟实验, 利用激振器发出前述雷克子波, 作为对 1.5m×1.2m 的石膏板小样品的激励信号, 通过石膏板上布设的五个加速度传感器采集微小扰动发生前后的波形信号.石膏板的波速与石膏板的含水率有关, 含水率增加, 板中超声波传播速度相应减小, 因此, 本实验采用加湿的方式对样品介质添加扰动.

将石膏板剖分为 1 0×10 的网格, 则待反演量的个数为100, 网格单元长宽分别为 0.15 m和0.12 m, 在石膏板的中心位置加水以增加中心点处的含水率作为局部扰动; 尾波分析时段为0.2—2.8 s,窗口长度均为 0.6 s, 重叠 0.4 s, 共 11 段尾波数据,则可获得 5 ×11=55 个时延数据; 敏感核构造中的散射系数D=140m2/s ; 线性最小二乘法中的相对 长度λ=0.1m , 网格 间距λ0取0.134 m.

图5 算例 1 (a) 线性最小二乘法的反演图像; (b) 本文方法的反演图像Fig.5.Case 1:(a) Inversion image of linear least squares method; (b) inversion image of the method in this paper.

图6 算例 2 (a) 线性最小二乘法的反演图像; (b) 本文方法的反演图像Fig.6.Case 2:(a) Inversion image of linear least squares method; (b) inversion image of the method in this paper.

图7 实验数据处理结果 (a) 线性最小二乘法的反演图像; (b) 本文方法的反演图像Fig.7.The results of experimental data processing:(a) Inversion image of linear least squares method; (b)inversion image of the method in this paper.

在该模拟实验中, 仅仅对板的中心处进行了加湿, 扰动的空间分布已足够稀疏, 因此在使用稀疏重构方法时, 可以直接求解(6)式的最小l1范数,而无需进行稀疏变换, 图7(b)即为稀疏重构法的成像结果, 图7(a)为线性最小二乘法的迭代结果,L曲 线 法 求 得σm为 0.33698.对 比 7(a)和 图7(b)可以看出, 两者对于扰动位置均有不错的定位效果, 但稀疏重构法的反演结果中, 扰动更为聚焦,范围更接近于实际扰动的施加情况.

前面讨论了单点扰动算例, 并验证了所提方法的优越性能.在实际中, 往往会出现多个扰动的情况, 因此算例3至5考虑了多个扰动区域的仿真,图8—图10分别给出了两种算法反演结果的对比,图中的蓝色方框和黄色方框分别表示在该区域内施加+0.5%和–0.5%的速度变化.

图8 算例 3 (a) 线性最小二乘法的反演图像; (b) 本文方法的反演图像Fig.8.Case 3:(a) Inversion image of linear least squares method; (b) inversion image of the method in this paper.

图9 算例 4 (a) 线性最小二乘法的反演图像; (b) 本文方法的反演图像Fig.9.Case 4:(a) Inversion image of linear least squares method; (b) inversion image of the method in this paper.

图10 算例 5 (a) 线性最小二乘法的反演图像; (b) 本文方法的反演图像Fig.10.Case 5:(a) Inversion image of linear least squares method; (b) inversion image of the method in this paper.

图8(a), 图9(a)和图10(a)分别是算例 3, 4,5用线性最小二乘法迭代10次的结果, L曲线法获 取 的σm分 别 为 0.01244, 0.00201和 0.00208,图8(b), 图9(b)和图10(b)分别为算例 3, 4, 5 用本文方法迭代5次、5次和4次的结果.从图9和图10可以直观地看出, 本文方法对于多扰动区域的反演质量明显优于线性最小二乘法, 线性最小二乘法在算例4和算例5的结果中, 仅能定位出部分扰动, 且其位置与范围也有一定程度的偏离真值,而本文方法对于设定的多扰动区域, 均能较为准确地反演出其位置与范围; 而从算例3的反演图像来看, 图8(a)确定的扰动区域与真值相差甚远, 甚至失去了参考价值, 相比之下, 图8(b)的结果更为接近真值, 有利于正确分析引起波速变化的物理机制.

表1 线性最小二乘法与本文方法计算成像时间对比Table 1.The comparison of imaging time between linear least squares method and the method in this paper.

此外, 两种方法在5个算例下反演成像时间对比如表1所列.

从表1可以看出, 在5个算例中, 本文方法的成像时间均小于线性最小二乘法, 平均减少了86%左右.同时, 模拟实验中两种方法的成像时间分别为 1.075540 s和 0.122640 s, 稀疏重构法的成像时间比线性最小二乘法少了88.6%.

综合以上数值仿真及实验结果可知, 压缩感知理论下的尾波干涉成像方法的性能优于线性最小二乘差分法, 不仅在多扰动识别能力、定位精度上有所提高, 成像时间也大为缩短.

5 结 论

在尾波干涉成像中, 现有的线性最小二乘差分成像方法定位精度低, 且在多扰动区域的识别中几乎失效.针对该问题, 充分利用扰动区域的空间稀疏性, 提出了一种基于稀疏重构的尾波干涉成像新方法.该方法首先根据尾波干涉法获取的时延数据以及扩散近似敏感核矩阵建立反演成像的欠定方程, 其次在压缩感知框架下, 通过信号的稀疏表示进一步构造反演方程, 并利用l1范数最小化重构稀疏信号, 实现了对速度扰动空间分布的定位.与线性最小二乘差分法相比, 稀疏重构方法避免了复杂的参数确定操作, 克服了多扰动点识别困难、定位不准确等问题, 在保证反演精度的前提下能够明显减少计算时间.仿真结果表明, 稀疏重构方法是一个性能优越、稳定可靠的优秀算法, 有助于加快尾波干涉成像技术在地震火山预测、材料无损检测等领域的实际应用.

需要特别指出的是, 线性最小二乘法和稀疏重构方法对于相对扰动大小的定量分析效果均不佳,敏感核的质量是主要影响因素.因此, 如何构造能够准确描述尾波时空传播特性的敏感核函数是尾波干涉成像的一个重要方面, 也是后续有待研究的问题.同时, 本文将敏感核矩阵作为观测矩阵, 其特性对稀疏重构算法性能的影响有待后续着重研究.另外, 本文中的模拟实验在介质尺寸与边界条件等方面与数值仿真环境差距较大, 后续将设计实施更为合理的实验方案并进一步研究多扰动尾波干涉成像问题.