搜索

文章查询

x

留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

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

张涛 侯宏 鲍明

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

张涛, 侯宏, 鲍明
PDF
HTML
导出引用
导出核心图
  • 尾波干涉成像是利用尾波时延和扩散近似敏感核来反演散射介质中微小速度扰动空间分布的技术. 该问题本质上是一个欠定问题, 一般无确定解, 导致其难以精确定位介质中微小波速变化的区域. 为解决上述缺陷, 本文利用速度扰动分布在空间上具有稀疏性的特点, 提出了一种基于压缩感知理论中稀疏重构算法的尾波干涉成像方法. 该方法可在散射介质中较准确地获取速度扰动的空间位置和范围, 同时具有较高的计算效率. 数值仿真和实验结果表明: 相比于现有的线性最小二乘差分成像方法, 无论是单个还是多个扰动区域, 该方法均能更精确地进行定位成像, 同时明显减少了计算时间.
      通信作者: 侯宏, houhong@nwpu.edu.cn
    • 基金项目: 国家重点研发计划项目(批准号: 2016YFF0200902)和国家自然科学基金(批准号: 11474230, 11704314)资助的课题
    [1]

    Snieder R 2002 Phys. Rev. E 66 046615

    [2]

    Snieder R, Grêt A, Douma H, Scales J 2002 Science 295 2253

    [3]

    Snieder R 2006 Pure. Appl. Geophys. 163 455

    [4]

    Grêt A, Snieder R, Aster R C, Kyle P R 2005 Geophys. Res. Lett. 32 233

    [5]

    Grêt A, Snieder R, Scales J 2006 J. Geophys. Res. 111 B03305

    [6]

    Larose E, Hall S 2009 J. Acoust. Soc. Am. 125 1853

    [7]

    Niu F, Silver P G, Daley T M, Cheng X, Majer E L 2008 Nature 454 204

    [8]

    宋丽莉, 葛洪魁, 郭志伟, 王小琼 2012 岩石力学与工程学报 31 713

    Song L L, Ge H K, Guo Z W, Wang X Q 2012 Chin. J. Rock. Mech. Eng. 31 713

    [9]

    谢凡, 任雅琼, 王宝善 2017 地球物理学报 60 1470

    Xie F, Ren Y Q, Wang B S 2017 Chin. J. Geophys. 60 1470

    [10]

    Pacheco C, Snieder R 2005 J. Acoust. Soc. Am. 118 1300

    [11]

    Obermann A, Planès T, Larose E, Campillo M 2013 J. Geophys. Res. 118 6285

    [12]

    Lesage P, Reyes G, Arámbula R 2014 J. Geophys. Res. 119 4360

    [13]

    Rossetto V, Margerin L, Planès T, Larose E 2011 J. Appl. Phys. 109 034903

    [14]

    Zhang Y, Planès T, Larose E, Obermann A, Rospars C, Moreau G 2016 J. Acoust. Soc. Am. 139 1691

    [15]

    Donoho D L 2006 IEEE Trans. Inform. Theory 52 1289

    [16]

    Candes E J, Romberg J 2006 IEEE Trans. Inform. Theory 52 489

    [17]

    宁方立, 何碧静, 韦娟 2013 物理学报 62 174212

    Ning F L, He B J, Wei J 2013 Acta Phys. Sin. 62 174212

    [18]

    李龙珍, 姚旭日, 刘雪峰, 俞文凯, 翟光杰 2014 物理学报 63 224201

    Li L Z, Yao X R, Liu X F, Yu W K, Zhai G J 2014 Acta Phys. Sin. 63 224201

    [19]

    时洁, 杨德森, 时胜国, 胡博, 朱中锐 2016 物理学报 65 024302

    Shi J, Yang D S, Shi S G, Hu B, Zhu Z R 2016 Acta Phys. Sin. 65 024302

    [20]

    Duarte M F, Davenport M A, Takhar D, Laska J N, Sun T, Kelly K F, Baraniuk R G 2008 IEEE Sig. Proc. Mag. 25 83

    [21]

    Lingala S G, Hu Y, Dibella E, Jacob M 2011 IEEE Trans. Med. Imaging 30 1042

    [22]

    Tropp J A, Gilbert A C 2007 IEEE Trans. Inform. Theory 53 4655

    [23]

    Mikesell T D, Malcolm A E, Yang D, Haney M M 2015 Geophys. J. Int. 202 347

    [24]

    Hadziioannou C, Larose E, Coutant O, Roux P, Campillo M 2009 J. Acoust. Soc. Am. 125 3688

    [25]

    Hansen P 1992 SIAM Rev. 34 561

    [26]

    Candes E J 2008 Comptes Rendus Mathematique 346 589

    [27]

    Chen J G http://www.paper.edu.cn/releasepaper/content/200606-478 [2006-6-28]

  • 图 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.

    图 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.

    图 3  二维速度场模型

    Fig. 3.  2-D velocity field model.

    图 4  激励源及接收点布设

    Fig. 4.  Layout of the source and receivers.

    图 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.

    图 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.

    图 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.

    图 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.

    表 1  线性最小二乘法与本文方法计算成像时间对比

    Table 1.  The comparison of imaging time between linear least squares method and the method in this paper.

    反演成像方法成像时间/s
    算例1算例2算例3算例4算例5
    线性最小二乘法1.5848501.7935171.6012761.6709322.278217
    本文方法0.2648940.2985830.2535110.2689690.115788
    下载: 导出CSV
  • [1]

    Snieder R 2002 Phys. Rev. E 66 046615

    [2]

    Snieder R, Grêt A, Douma H, Scales J 2002 Science 295 2253

    [3]

    Snieder R 2006 Pure. Appl. Geophys. 163 455

    [4]

    Grêt A, Snieder R, Aster R C, Kyle P R 2005 Geophys. Res. Lett. 32 233

    [5]

    Grêt A, Snieder R, Scales J 2006 J. Geophys. Res. 111 B03305

    [6]

    Larose E, Hall S 2009 J. Acoust. Soc. Am. 125 1853

    [7]

    Niu F, Silver P G, Daley T M, Cheng X, Majer E L 2008 Nature 454 204

    [8]

    宋丽莉, 葛洪魁, 郭志伟, 王小琼 2012 岩石力学与工程学报 31 713

    Song L L, Ge H K, Guo Z W, Wang X Q 2012 Chin. J. Rock. Mech. Eng. 31 713

    [9]

    谢凡, 任雅琼, 王宝善 2017 地球物理学报 60 1470

    Xie F, Ren Y Q, Wang B S 2017 Chin. J. Geophys. 60 1470

    [10]

    Pacheco C, Snieder R 2005 J. Acoust. Soc. Am. 118 1300

    [11]

    Obermann A, Planès T, Larose E, Campillo M 2013 J. Geophys. Res. 118 6285

    [12]

    Lesage P, Reyes G, Arámbula R 2014 J. Geophys. Res. 119 4360

    [13]

    Rossetto V, Margerin L, Planès T, Larose E 2011 J. Appl. Phys. 109 034903

    [14]

    Zhang Y, Planès T, Larose E, Obermann A, Rospars C, Moreau G 2016 J. Acoust. Soc. Am. 139 1691

    [15]

    Donoho D L 2006 IEEE Trans. Inform. Theory 52 1289

    [16]

    Candes E J, Romberg J 2006 IEEE Trans. Inform. Theory 52 489

    [17]

    宁方立, 何碧静, 韦娟 2013 物理学报 62 174212

    Ning F L, He B J, Wei J 2013 Acta Phys. Sin. 62 174212

    [18]

    李龙珍, 姚旭日, 刘雪峰, 俞文凯, 翟光杰 2014 物理学报 63 224201

    Li L Z, Yao X R, Liu X F, Yu W K, Zhai G J 2014 Acta Phys. Sin. 63 224201

    [19]

    时洁, 杨德森, 时胜国, 胡博, 朱中锐 2016 物理学报 65 024302

    Shi J, Yang D S, Shi S G, Hu B, Zhu Z R 2016 Acta Phys. Sin. 65 024302

    [20]

    Duarte M F, Davenport M A, Takhar D, Laska J N, Sun T, Kelly K F, Baraniuk R G 2008 IEEE Sig. Proc. Mag. 25 83

    [21]

    Lingala S G, Hu Y, Dibella E, Jacob M 2011 IEEE Trans. Med. Imaging 30 1042

    [22]

    Tropp J A, Gilbert A C 2007 IEEE Trans. Inform. Theory 53 4655

    [23]

    Mikesell T D, Malcolm A E, Yang D, Haney M M 2015 Geophys. J. Int. 202 347

    [24]

    Hadziioannou C, Larose E, Coutant O, Roux P, Campillo M 2009 J. Acoust. Soc. Am. 125 3688

    [25]

    Hansen P 1992 SIAM Rev. 34 561

    [26]

    Candes E J 2008 Comptes Rendus Mathematique 346 589

    [27]

    Chen J G http://www.paper.edu.cn/releasepaper/content/200606-478 [2006-6-28]

  • [1] 李龙珍, 姚旭日, 刘雪峰, 俞文凯, 翟光杰. 基于压缩感知超分辨鬼成像. 物理学报, 2014, 63(22): 224201. doi: 10.7498/aps.63.224201
    [2] 庄佳衍, 陈钱, 何伟基, 冒添逸. 基于压缩感知的动态散射成像. 物理学报, 2016, 65(4): 040501. doi: 10.7498/aps.65.040501
    [3] 白旭, 李永强, 赵生妹. 基于压缩感知的差分关联成像方案研究. 物理学报, 2013, 62(4): 044209. doi: 10.7498/aps.62.044209
    [4] 李少东, 陈永彬, 刘润华, 马晓岩. 基于压缩感知的窄带高速自旋目标超分辨成像物理机理分析. 物理学报, 2017, 66(3): 038401. doi: 10.7498/aps.66.038401
    [5] 丁亚辉, 孙玉发, 朱金玉. 一种基于压缩感知的三维导体目标电磁散射问题的快速求解方法. 物理学报, 2018, 67(10): 100201. doi: 10.7498/aps.67.20172543
    [6] 王盼盼, 姚旭日, 刘雪峰, 俞文凯, 邱棚, 翟光杰. 基于行扫描测量的运动目标压缩成像. 物理学报, 2017, 66(1): 014201. doi: 10.7498/aps.66.014201
    [7] 仲亚军, 刘娇, 梁文强, 赵生妹. 针对多散斑图的差分压缩鬼成像方案研究. 物理学报, 2015, 64(1): 014202. doi: 10.7498/aps.64.014202
    [8] 刘扬阳, 吕群波, 曾晓茹, 黄旻, 相里斌. 静态计算光谱成像仪图谱反演的关键数据处理技术. 物理学报, 2013, 62(6): 060203. doi: 10.7498/aps.62.060203
    [9] 李广明, 吕善翔. 混沌信号的压缩感知去噪. 物理学报, 2015, 64(16): 160502. doi: 10.7498/aps.64.160502
    [10] 冯丙辰, 方晟, 张立国, 李红, 童节娟, 李文茜. 基于压缩感知理论的非线性γ谱分析方法. 物理学报, 2013, 62(11): 112901. doi: 10.7498/aps.62.112901
    [11] 宁方立, 何碧静, 韦娟. 基于lp范数的压缩感知图像重建算法研究. 物理学报, 2013, 62(17): 174212. doi: 10.7498/aps.62.174212
    [12] 王哲, 王秉中. 压缩感知理论在矩量法中的应用. 物理学报, 2014, 63(12): 120202. doi: 10.7498/aps.63.120202
    [13] 张新鹏, 胡茑庆, 程哲, 钟华. 基于压缩感知的振动数据修复方法. 物理学报, 2014, 63(20): 200506. doi: 10.7498/aps.63.200506
    [14] 时洁, 杨德森, 时胜国, 胡博, 朱中锐. 基于压缩感知的矢量阵聚焦定位方法. 物理学报, 2016, 65(2): 024302. doi: 10.7498/aps.65.024302
    [15] 康荣宗, 田鹏武, 于宏毅. 一种基于选择性测量的自适应压缩感知方法. 物理学报, 2014, 63(20): 200701. doi: 10.7498/aps.63.200701
    [16] 陈明生, 王时文, 马韬, 吴先良. 基于压缩感知的目标频空电磁散射特性快速分析. 物理学报, 2014, 63(17): 170301. doi: 10.7498/aps.63.170301
    [17] 李慧, 赵琳, 李亮. 基于贝叶斯压缩感知的周跳探测与修复方法. 物理学报, 2016, 65(24): 249101. doi: 10.7498/aps.65.249101
    [18] 冷雪冬, 王大鸣, 巴斌, 王建辉. 基于渐进添边的准循环压缩感知时延估计算法. 物理学报, 2017, 66(9): 090703. doi: 10.7498/aps.66.090703
    [19] 康志伟, 吴春艳, 刘劲, 马辛, 桂明臻. 基于两级压缩感知的脉冲星时延估计方法. 物理学报, 2018, 67(9): 099701. doi: 10.7498/aps.67.20172100
    [20] 胡耀华, 刘艳, 穆鸽, 秦齐, 谭中伟, 王目光, 延凤平. 基于多模光纤散斑的压缩感知在光学图像加密中的应用. 物理学报, 2020, 69(3): 034203. doi: 10.7498/aps.69.20191143
  • 引用本文:
    Citation:
计量
  • 文章访问数:  584
  • PDF下载量:  11
  • 被引次数: 0
出版历程
  • 收稿日期:  2019-05-28
  • 修回日期:  2019-06-25
  • 上网日期:  2019-11-26
  • 刊出日期:  2019-10-01

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

  • 1. 西北工业大学航海学院, 海洋声学信息感知工业和信息化部重点实验室, 西安 710072
  • 2. 中国科学院噪声与振动重点实验室(声学研究所), 北京 100190
  • 通信作者: 侯宏, houhong@nwpu.edu.cn
    基金项目: 国家重点研发计划项目(批准号: 2016YFF0200902)和国家自然科学基金(批准号: 11474230, 11704314)资助的课题

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

English Abstract

    • 尾波干涉(coda wave interferometry, CWI), 是一种利用介质中的多次散射波检测介质微小变化的测量方法. 多次散射波是声信号或振动信号在散射介质中传播时, 由于介质的非均匀性(存在如颗粒、裂缝、包体以及孔洞等非均质体)而发生多次散射、反射产生的, 在波形图中显示为直达波后面拖着的长长的“尾巴”, 因此也称尾波. 尾波由于在介质中来回传播, 对介质更密集地采样, 使得介质中的微小变化不断迭加放大, 从而对介质性质的微弱改变具有极高的敏感性, 检测出直达波无法识别的微小变化. 自2002年CWI法首次由Snieder等[13]提出后, 便在地球物理学和材料科学等领域得到了广泛的研究与应用[49]. 文献[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)近年来逐渐兴起并趋于成熟, 在雷达探测、图像处理、声呐定位等领域都有了广泛的应用, 在物理成像方面也表现出强大的生命力[1721]. 文献[17]采用分块压缩感知思想, 提出一种基于${l_p}$范数的压缩感知图像重建算法, 提高重建图像质量的同时缩短了成像时间. 文献[18]提出了基于稀疏测量的超分辨压缩感知鬼成像重建模型, 并搭建了实验装置, 验证了该模型的有效性与优越性, 推动了鬼成像的实用进展. 文献[19]利用声源的空间稀疏性, 在压缩感知理论下建立了矢量阵聚焦定位新方法, 克服了相干声源分辨困难、贡献评价不准确等问题, 且定位精度高, 背景抑制能力强. CS理论指出, 若某信号是可压缩的或经某种变换后稀疏, 则可通过随机观测矩阵(与变换基不相关)将高维信号投影至低维空间中, 最终通过求解优化问题从少量低维数据高概率地重构原始高维信号. 而在尾波干涉成像问题中, 实际扰动通常是空间上聚焦的, 扰动的空间区域相比于监测的整体介质区域往往是稀疏的, 恰好满足CS理论对于信号稀疏性的要求. 本文基于CS理论, 提出了一种尾波干涉成像的新的解算方法. 该方法首先通过CWI和扩散波敏感核构建了用于成像的反演模型; 其次引入CS理论框架下的稀疏变换将欠定方程重新构造; 最终通过压缩感知重构算法——正交匹配追踪算法(orthogonal matching pursuit, OMP)[22]求解该欠定方程, 并获得速度扰动的空间成像. 此方法与先前的线性最小二乘差分成像方法相比, 在单扰动和多扰动情况下对于扰动区域位置和范围的确定都更加准确, 且执行简易, 没有复杂的参数选取, 同时还具有更高的计算效率. 最后结合具体仿真算例, 给出了与线性最小二乘差分成像方法的比较结果, 证明了该成像方法的可行性、有效性和高效率.

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

      图  1  多散射介质中扰动前后波形的比较 (a) 直达波扰动前后的波形; (b)尾波扰动前后的波形

      Figure 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.

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

      ${R^{(t,T)}}({t_{\rm{s}}}) \equiv \frac{{\displaystyle\int_{t - T}^{t + T} {{u_{{\rm{unp}}}}(t'){u_{{\rm{per}}}}(t' + {t_{\rm{s}}}){\rm{d}}t'} }}{{\sqrt {\displaystyle\int_{t - T}^{t + T}\!\!{u_{{\rm{unp}}}^2(t'){\rm{d}}t' \displaystyle\int_{t - T}^{t + T}\!\!{u_{{\rm{per}}}^2(t'){\rm{d}}t'} } } }}, $

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

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

      $\frac{{{\text{δ}}\upsilon }}{\upsilon } = - \frac{{{\text{δ}}t}}{t}, $

      式中${\text{δ}} \upsilon $表示波速扰动, 是整体介质波速变化的加权平均; $\upsilon $为扰动前的介质波速.

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

      $K({s_1},{s_2},{x_0},t) = \frac{{\displaystyle\int_0^t {p({s_1},{x_0},\tau )p({x_0},{s_2},t - \tau ){\rm{d}}\tau } }}{{p({s_1},{s_2},t)}}, $

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

      $p({s_1},{s_2},t) = \frac{1}{{4{\text{π}}Dt}}\exp \left( - \frac{{{{\left| {{s_1} - {s_2}} \right|}^2}}}{{4Dt}}\right), $

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

      图  2  基于扩散近似的二维敏感核示例 (a) t = 1 s时的敏感核空间分布; (b) t = 5 s时的敏感核空间分布; (c) t = 1 s时的敏感核俯视图; (d) t = 5 s时的敏感核俯视图

      Figure 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.

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

      ${\text{δ}}t(t) = - \int_S {K({s_1}} ,{s_2},{x_0},t)\frac{{{\text{δ}}\upsilon }}{\upsilon }({x_0}){\rm{d}}S({x_0}), $

      式中$\dfrac{{{\text{δ}}\upsilon }}{\upsilon }({x_0})$为介质空间上${x_0}$处的相对波速变化, 即为需要反演计算的量; ${\text{δ}}t(t)$${x_0}$处相对速度扰动$\dfrac{{{\text{δ}}\upsilon }}{\upsilon }({x_0})$引起的走时变化, 在物理意义上等效于尾波干涉求得的介质扰动前后同一接收点波形的时延. 因此, 尾波干涉成像的实质是通过实验数据处理得到的时延, 以及根据介质物理性质构造的敏感核, 反解出空间各处的波速变化, 该问题是一个典型的数学物理反问题.

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

      ${{T}} = {{GV}}, $

      式中, ${{T}} \in {R^{M,1}}$表示通过不同的激励接收对和不同的尾波观察窗口获取的M个时延数据; ${{V}} \in {R^{N,1}}$表示空间各网格点对应的相对波速扰动值, N为网格总数; ${{G}} \in {R^{M,N}}$表示${{T}}$${{V}}$的映射, 其矩阵元素${G_{ij}} = \Delta s{K_{ij}}$, 其中${K_{ij}}$为敏感核矩阵${{K}}$的元素, ${{K}} \in {R^{M,N}}$的每行代表一个敏感核, 每列代表一个网格, 每个网格的面积为$\Delta s$.

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

      ${{V}} = {{{V}}_0} + {{{C}}_m}{{{G}}^{\rm{T}}}{({{G}}{{{C}}_m}{{{G}}^{\rm{T}}} + {{{C}}_d})^{ - 1}}({{T}} - {{G}}{{{V}}_0}), $

      式中, ${{{V}}_0}$为先验初始值, 一般为零矩阵; ${{{G}}^{\rm{T}}}$表示${{G}}$的转置; ${{{C}}_d}$为测得数据的对角协方差矩阵; ${{{C}}_m}$为介质空间元素平滑矩阵, 可用下式计算:

      ${{{C}}_m}\left({s_1},{s_2}\right) = {\left({\sigma _m}\frac{{{\lambda _0}}}{\lambda }\right)^2} \cdot \frac{1}{{ch\left(\dfrac{{\left| {{s_1} - {s_2}} \right|}}{\lambda }\right)}},$

      式中, ${\sigma _m}$为先验标准差; ${\lambda _0}$为网格间距; $\lambda $为相关长度. ${\sigma _m}$$\lambda $一般通过L曲线法[25]选取.

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

      ${{C}}_m^{{\rm{new}}} = {{{C}}_m} - {{{C}}_m}{{{G}}^{\rm{T}}}{({{G}}{{{C}}_m}{{{G}}^{\rm{T}}} + {{{C}}_d})^{ - 1}}{{G}}{{{C}}_m}.$

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

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

      ${{x}} = {{{\psi }}^{\rm{T}}}{{a}}, $

      式中, ${{\psi }} \in {R^N}$为稀疏变换矩阵; ${{a}} \in {R^{N,1}}$${{x}}$的稀疏表示.

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

      ${{y}} = {{\phi }}{{a}} = {{\phi }}{{\psi x}} = {{\varTheta x}} ,$

      式中, ${{\phi }}\in {R^{M,N}}$为观测矩阵; ${{y}} \in {R^{M,1}}$${{a}}$在观测矩阵下的观测值; ${{\varTheta }} = {{\phi }}{{\psi }}$称为传感矩阵.

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

      $\mathop {\min }\limits_{{x}} {\left\| {{{\psi x}}} \right\|_0}{\rm{ s}}{\rm{.t}}{\rm{. }}{{y}} = {{\varTheta x}}.$

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

      $\mathop {\min }\limits_{{x}} {\left\| {{{\psi x}}} \right\|_1}{\rm{ s}}{\rm{.t}}{\rm{. }}{{y}} = {{\varTheta x}}.$

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

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

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

      ${{T}} = {{{G}}^{{\rm{CS}}}}{{{V}}^{{\rm{CS}}}}, $

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

      ${{V}} = {{{P}}^{ - 1}}{{{V}}^{{\rm{CS}}}}, $

      整个过程无须使用${{P}}$而仅须用到其逆矩阵${{{P}}^{ - 1}}$, 对于傅里叶变换矩阵, ${{{P}}^{ - 1}} = {{{P}}^{\rm{T}}}$.

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

      图  3  二维速度场模型

      Figure 3.  2-D velocity field model.

      图  4  激励源及接收点布设

      Figure 4.  Layout of the source and receivers.

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

      $\frac{{{\partial ^2}p}}{{\partial {x^2}}} + \frac{{{\partial ^2}p}}{{\partial {y^2}}} + f(x,y,t) = \frac{1}{{{v^2}(x,y)}}\frac{{{\partial ^2}p}}{{\partial {t^2}}}, $

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

      $\begin{split} &{p^{k + 1}}(i,j) \\= & \, 2{p^k}(i,j) - {p^{k - 1}}(i,j)\\ & + \frac{{{v^2}\Delta {t^2}}}{{\Delta {h^2}}}\left\{ { - \frac{1}{{560}}} \right.\left[ {{p^k}(i - 4,j) + {p^k}(i + 4,j)} \right]\\ & + \frac{8}{{315}}\left[ {{p^k}(i - 3,j) + {p^k}(i + 3,j)} \right]\\ & - \frac{1}{5}\left[ {{p^k}(i - 2,j) + {p^k}(i + 2,j)} \right]\\ & +\left. \frac{8}{5}\left[ {{p^k}(i - 1,j) + {p^k}(i + 1,j)} \right] - \frac{{205}}{{72}}{p^k}(i,j)\right\}\\ & + \frac{{{v^2}\Delta {t^2}}}{{\Delta {h^2}}}\left\{ { - \frac{1}{{560}}} \right.\left[ {{p^k}(i,j - 4) + {p^k}(i,j + 4)} \right]\\ & + \frac{8}{{315}}\left[ {{p^k}(i,j - 3) + {p^k}(i,j + 3)} \right]\\ & - \frac{1}{5}\left[ {{p^k}(i,j - 2) + {p^k}(i,j + 2)} \right]\\ & +\left. \frac{8}{5}\left[ {{p^k}(i,j - 1) + {p^k}(i,j + 1)} \right] - \frac{{205}}{{72}}{p^k}(i,j)\right\}\\ & + f(t) * \delta (i - {i_0}) * \delta (i - {j_0})\\ & i = 1,2, \cdots, nx - 1;j = 1,2, \cdots, \\ & ny - 1;k = 2,3, \cdots, nt \end{split}$

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

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

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

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

      图  5  算例1 (a) 线性最小二乘法的反演图像; (b) 本文方法的反演图像

      Figure 5.  Case 1: (a) Inversion image of linear least squares method; (b) inversion image of the method in this paper.

      图  6  算例2 (a) 线性最小二乘法的反演图像; (b) 本文方法的反演图像

      Figure 6.  Case 2: (a) Inversion image of linear least squares method; (b) inversion image of the method in this paper.

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

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

      图  7  实验数据处理结果 (a) 线性最小二乘法的反演图像; (b) 本文方法的反演图像

      Figure 7.  The results of experimental data processing: (a) Inversion image of linear least squares method; (b)inversion image of the method in this paper.

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

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

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

      图  8  算例3 (a) 线性最小二乘法的反演图像; (b) 本文方法的反演图像

      Figure 8.  Case 3: (a) Inversion image of linear least squares method; (b) inversion image of the method in this paper.

      图  10  算例5 (a) 线性最小二乘法的反演图像; (b) 本文方法的反演图像

      Figure 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曲线法获取的${\sigma _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)的结果更为接近真值, 有利于正确分析引起波速变化的物理机制.

      图  9  算例4 (a) 线性最小二乘法的反演图像; (b) 本文方法的反演图像

      Figure 9.  Case 4: (a) Inversion image of linear least squares method; (b) inversion image of the method in this paper.

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

      反演成像方法成像时间/s
      算例1算例2算例3算例4算例5
      线性最小二乘法1.5848501.7935171.6012761.6709322.278217
      本文方法0.2648940.2985830.2535110.2689690.115788

      表 1  线性最小二乘法与本文方法计算成像时间对比

      Table 1.  The comparison of imaging time between linear least squares method and the method in this paper.

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

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

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

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

参考文献 (27)

目录

    /

    返回文章
    返回