删除或更新信息，请邮件至freekaoyan#163.com(#换成@)

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

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

English Abstract


全文HTML

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

2.1.尾波干涉基本原理

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

尾波干涉理论的实质是通过测量两个不同时刻尾波列的相位差来获取介质在该时段的波速变化[1—3]. 互相关法和延展法是目前基于此理论来提取介质波速变化的主要方法[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 $为扰动前的介质波速.
2
2.2.敏感核
-->

2.2.敏感核

为了获取局部扰动的空间分布情况, 需要将上述尾波干涉求得的总体平均时延${\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时的敏感核俯视图
Figure2. 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.

2
2.3.反演成像模型
-->

2.3.反演成像模型

基于扩散近似敏感核, 文献[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}.$

3.基于稀疏重构的尾波干涉成像方法
压缩感知理论指出, 假定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}}}$.
4.数值仿真与结果
为了验证基于稀疏重构算法的成像方法, 现通过数值仿真进行分析. 本数值仿真使用的激励源为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 二维速度场模型
Figure3. 2-D velocity field model.

图 4 激励源及接收点布设
Figure4. 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) 本文方法的反演图像
Figure5. Case 1: (a) Inversion image of linear least squares method; (b) inversion image of the method in this paper.

图 6 算例2　(a) 线性最小二乘法的反演图像; (b) 本文方法的反演图像
Figure6. 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) 本文方法的反演图像
Figure7. 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) 本文方法的反演图像
Figure8. Case 3: (a) Inversion image of linear least squares method; (b) inversion image of the method in this paper.

图 10 算例5　(a) 线性最小二乘法的反演图像; (b) 本文方法的反演图像
Figure10. 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) 本文方法的反演图像
Figure9. 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线性最小二乘法与本文方法计算成像时间对比
Table1.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%.
综合以上数值仿真及实验结果可知, 压缩感知理论下的尾波干涉成像方法的性能优于线性最小二乘差分法, 不仅在多扰动识别能力、定位精度上有所提高, 成像时间也大为缩短.
5.结　论
在尾波干涉成像中, 现有的线性最小二乘差分成像方法定位精度低, 且在多扰动区域的识别中几乎失效. 针对该问题, 充分利用扰动区域的空间稀疏性, 提出了一种基于稀疏重构的尾波干涉成像新方法. 该方法首先根据尾波干涉法获取的时延数据以及扩散近似敏感核矩阵建立反演成像的欠定方程, 其次在压缩感知框架下, 通过信号的稀疏表示进一步构造反演方程, 并利用${l_1}$范数最小化重构稀疏信号, 实现了对速度扰动空间分布的定位. 与线性最小二乘差分法相比, 稀疏重构方法避免了复杂的参数确定操作, 克服了多扰动点识别困难、定位不准确等问题, 在保证反演精度的前提下能够明显减少计算时间. 仿真结果表明, 稀疏重构方法是一个性能优越、稳定可靠的优秀算法, 有助于加快尾波干涉成像技术在地震火山预测、材料无损检测等领域的实际应用.
需要特别指出的是, 线性最小二乘法和稀疏重构方法对于相对扰动大小的定量分析效果均不佳, 敏感核的质量是主要影响因素. 因此, 如何构造能够准确描述尾波时空传播特性的敏感核函数是尾波干涉成像的一个重要方面, 也是后续有待研究的问题. 同时, 本文将敏感核矩阵作为观测矩阵, 其特性对稀疏重构算法性能的影响有待后续着重研究. 另外, 本文中的模拟实验在介质尺寸与边界条件等方面与数值仿真环境差距较大, 后续将设计实施更为合理的实验方案并进一步研究多扰动尾波干涉成像问题.
相关话题/空间 信号 图像 计算 实验
  • 领限时大额优惠券,享本站正版考研考试资料!
    大额优惠券
    优惠券领取后72小时内有效，10万种最新考研考试考证类电子打印资料任你选。涵盖全国500余所院校考研专业课、200多种职业资格考试、1100多种经典教材，产品类型包含电子书、题库、全套资料以及视频，无论您是考研复习、考证刷题，还是考前冲刺等，不同类型的产品可满足您学习上的不同需求。 ...
    考试优惠券 本站小编 Free壹佰分学习网 2022-09-19
  • 交变力磁力显微镜: 在三维空间同时观测静态和动态磁畴
    摘要:传统磁力显微镜(MFM)的磁畴扫描是利用激光束反射探测探针和样品之间的静磁力.因此,对于MFM,直接探测样品在交流磁场作用下的动态磁力仍然是一个挑战.交变力磁力显微镜(A-MFM)使用Co-GdOx超顺磁探针可以实现在交流磁场(频率ωm)作用下探测动态磁力.采集ωm和2ωm信号能准确地表示出样 ...
    中科院物理研究所 本站小编 Free考研考试 2021-12-29
  • 双旋光双反射结构的温度-辐射自稳定性原理和实验研究
    摘要:推导了双旋光双反射结构的反射光偏振态方程,仿真了在环境影响下的反射光偏振态变化情况,发现双旋光双反射结构的偏振无关反射自稳定性.从温度和辐射两个方面实验验证了双旋光双反射结构的偏振无关反射自稳定性.实验结果表明,当温度在–45℃—85℃之间变化时,双旋光双反射结构反射光的平均偏振保持度都能够达 ...
    中科院物理研究所 本站小编 Free考研考试 2021-12-29
  • 基于麦克风的气体超声分子束飞行速度的实验研究
    摘要:超声分子束的膨胀和输运过程是一个较为复杂的分子动力学问题,相关的参数较难准确计算.本文基于麦克风测量方法研究了多种气体(H2,D2,N2,Ar,He,CH4)超声分子束在自由膨胀过程中的平均速度及其沿出射方向在远域空间(喷射距离/喷嘴直径>310)的演变情况,获得了较大范围内分子束平均速度分布 ...
    中科院物理研究所 本站小编 Free考研考试 2021-12-29
  • 钙钛矿太阳能电池研究进展: 空间电势与光电转换机制
    摘要:钙钛矿太阳能电池具有高光电转换效率和低成本制备的特点,是极具希望实现大规模应用的下一代光伏技术.然而,对该类器件的光电转换过程的认知仍然不够清晰,相关研究难以直接观测器件内部的空间电势及其对光生电荷载流子的影响.开尔文探针力显微镜技术能够直接探测出器件空间电势的分布,进而直接反映器件工作的状态 ...
    中科院物理研究所 本站小编 Free考研考试 2021-12-29
  • D-T中子诱发贫化铀球壳内裂变率分布实验
    摘要:中子诱发裂变反应率是表征和检验中子在材料中的输运、裂变放能等过程的重要物理量.贫化铀球壳裂变反应率径向分布数据可为铀核数据宏观检验及研究裂变放能与贫化铀球壳厚度的关系提供数据支持.本文设计了内径为13.1cm,外径分别为18.10,19.40,23.35,25.40,28.45cm的五种不同厚 ...
    中科院物理研究所 本站小编 Free考研考试 2021-12-29
  • 基于有效介质理论的物理性能计算模型的软件实现
    摘要:在改进的有效介质理论的基础上采用C++/Qt混合编程,设计并开发出一套复合材料物理性能模拟计算软件—CompositeStudio.该软件通过格林函数对本构方程进行求解,计算体积分数、颗粒长径比、取向分布、宏观位向对复合材料有效性能的影响.目前软件开发了弹性模量和介电常数两个模块,提供了友好的 ...
    中科院物理研究所 本站小编 Free考研考试 2021-12-29
  • 基于软件定义量子通信的自由空间量子通信信道参数自适应调整策略
    摘要:自由空间量子通信会受到雾霾、沙尘、降雨等自然环境的干扰.为提升环境干扰下量子通信的性能,本文提出了基于软件定义量子通信(softwaredefinedquantumcommunication,SDQC)的自由空间量子通信信道参数自适应调整策略.该策略通过对环境状态实时监测,根据预置在应用层的程 ...
    中科院物理研究所 本站小编 Free考研考试 2021-12-29
  • 基于半导体光纤环形腔激光器的全光广播式超宽带信号源
    摘要:提出一种新型的基于半导体光纤环形腔激光器(semiconductorfiberringlaser,SFRL)全光超宽带(ultra-wideband,UWB)信号源的方案,该方案可以同时产生3路高斯脉冲一阶导数脉冲(monocycle)UWB信号.建立了这种全光UWB信号源的宽带理论模型,通过 ...
    中科院物理研究所 本站小编 Free考研考试 2021-12-29
  • 英国散裂中子源工程材料原位加载衍射实验高温样品环境优化设计
    摘要:英国散裂中子源(ISIS)在工程材料中子衍射领域有着十余年丰富的研究经验,最为典型的衍射谱仪之一的Engin-X在材料、加工等方向有着广泛应用,包括残余应力分布测量、金属相变分析、微观力学研究等.Engin-X通过设置红外加热型高温炉配套材料试验机的样品环境以实现中子衍射原位高温力学实验,目前 ...
    中科院物理研究所 本站小编 Free考研考试 2021-12-29
  • 飞秒脉冲抽运掺镱微结构光纤产生超连续谱的实验研究
    摘要:本文利用钛蓝宝石飞秒激光器抽运自制的掺镱微结构光纤,对微结构光纤中的非线性效应及超连续谱产生机理进行了实验研究.研究发现,当抽运光偏离Yb3+吸收最高峰85nm时,仍具有较高的发光效率.在飞秒脉冲抽运下,位于反常色散区的发射光首先被位于正常色散区的抽运光激发、放大并俘获,然后演化为超短脉冲,随 ...
    中科院物理研究所 本站小编 Free考研考试 2021-12-29

Free考研考试FreeKaoYan.Com
欢迎来到Free考研考试，"为实现人生的Free而奋斗"

© 2020 FreeKaoYan! . All rights reserved.