CN110348169B - 一种基于压缩感知理论的尾波干涉成像方法 - Google Patents

Abstract

本发明涉及一种基于压缩感知理论的尾波干涉成像方法，由于现有的线性最小二乘差分反演方法存在求解精度低、参数选取困难、计算效率低、难以定位多个扰动点的不足，同时考虑到扰动的空间区域相比于整体介质区域是稀疏的，恰好满足压缩感知理论对于信号稀疏性的要求。本发明所提供的基于压缩感知理论的尾波干涉成像方法，不需要进行复杂的参数选取操作，简单易行，计算效率高，求解结果稳定精确，并且对于同时存在多个扰动区域的情况下，仍能较为准确地获取速度扰动的空间位置和范围，相比于现有的反演方法，在实际的工程应用情况中更为适用。

Description

一种基于压缩感知理论的尾波干涉成像方法

技术领域

本发明属于尾波干涉成像方法，涉及一种基于压缩感知理论的尾波干涉成像方法，涉及如何对尾波干涉成像的反问题进行精确快速求解的问题。

背景技术

尾波干涉成像是利用尾波时延和扩散近似敏感核来反演散射介质中微小速度扰动空间分布的技术。该技术本质上是对一个矩阵形式的反问题进行求解，该问题是一个欠定问题，有无穷多解，同时具有病态性，难以精确求解。目前，对于该反问题的求解，国内外研究中均采用地球物理反演中常用的线性最小二乘差分反演方法，主要步骤如下：

假设通过尾波干涉技术和扩散近似敏感核理论得到反演模型的线性方程组，并写成矩阵形式如下

T＝GV (1)

式中，T∈R^M,1表示通过尾波干涉技术获取的M个时延数据；V∈R^N,1表示空间各网格点对应的相对波速扰动值，N为网格总数；G∈R^M,N表示T到V的映射，其矩阵元素G_ij＝ΔsK_ij，其中K_ij为敏感核矩阵K的元素，K∈R^M,N的每行代表一个敏感核，每列代表一个网格，每个网格的面积为Δs。

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

V＝V₀+C_mG^T(GC_mG^T+C_d)^-1(T-GV₀) (2)

式中，V₀为先验初始值，一般为零矩阵；G^T表示G的转置；C_d为测得数据的对角协方差矩阵；C_m为介质空间元素平滑矩阵，可用下式计算：

式中，σ_m为先验标准差；λ₀为网格间距；λ为相关长度。σ_m和λ一般通过L曲线法选取。s₁和s₂分别为激励点和接收点位置。

在用线性最小二乘差分法反演计算时，还需要循环迭代C_m以寻找最优的V值，迭代公式如下：

该方法存在求解精度不高、参数选取困难、计算效率低的缺陷，并且难以定位多个扰动点，在实际工程中的应用十分困难。

现有的线性最小二乘差分反演方法存在求解精度低、参数选取困难、计算效率低、难以定位多个扰动点的不足。

发明内容

要解决的技术问题

为了避免现有技术的不足之处，本发明提出一种基于压缩感知理论的尾波干涉成像方法，无论是存在单个还是多个扰动区域的情况下，该方法均能较准确地获取速度扰动的空间位置和范围，同时操作简易，具有较高的计算效率，能够更好地应用于工程实践中。

技术方案

一种基于压缩感知理论的尾波干涉成像方法，其特征在于步骤如下：

步骤1：通过尾波干涉技术获取时延数据，并将其存储为列向量；

步骤2：根据观测介质的物理性质构造敏感核矩阵，通过敏感核矩阵将尾波干涉法获取的时延数据与相应的局部变化联系起来，建立所示矩阵形式的反问题：

T＝GV

式中，T∈R^M,1表示通过尾波干涉技术获取的M个时延数据；V∈R^N,1表示空间各网格点对应的相对波速扰动值，N为网格总数；G∈R^M,N表示T到V的映射，其矩阵元素G_ij＝ΔsK_ij，其中K_ij为敏感核矩阵K的元素，K∈R^M,N的每行代表一个敏感核，每列代表一个网格，每个网格的面积为Δs；

步骤3：将G、T和V分别视为观测矩阵、测量值矩阵和真实信号，引入稀疏变换矩阵P对V进行离散稀疏变换：

T＝G^CSV^CS， (5)

式中，G^CS＝GP^-1，V^CS＝PV；

步骤4：将T＝G^CSV^CS转化为最优化问题，压缩感知重构算法求解出稀疏解V^CS：

步骤5：通过稀疏逆变换求得原始解V为：

V＝P^-1V^CS， (7)

所观测的二维散射介质空间各处的相对波速变化是通过划分网格的方式来离散表达的，将二维散射介质沿其长和宽剖分为n×n＝N个网格，每个独立网格空间内的相对波速变化视为一致；因此，网格划分越细致，即网格划分数量越多，越能精确地表达实际连续的介质空间各处相对波速变化情况。

将求得的N维列向量展开为n×n的矩阵，n×n＝N，并将该矩阵导入MATLAB软件中的surf函数产生二维图像，此图像即为介质在扰动前后相对速度变化的二维空间分布情况。

所述稀疏变换矩阵为：离散傅里叶变换矩阵、离散余弦变换矩阵或离散小波变换矩阵。

有益效果

本发明提出的一种基于压缩感知理论的尾波干涉成像方法，由于现有的线性最小二乘差分反演方法存在求解精度低、参数选取困难、计算效率低、难以定位多个扰动点的不足，同时考虑到扰动的空间区域相比于整体介质区域是稀疏的，恰好满足压缩感知理论对于信号稀疏性的要求。本发明所提供的基于压缩感知理论的尾波干涉成像方法，不需要进行复杂的参数选取操作，简单易行，计算效率高，求解结果稳定精确，并且对于同时存在多个扰动区域的情况下，仍能较为准确地获取速度扰动的空间位置和范围，相比于现有的反演方法，在实际的工程应用情况中更为适用。

附图说明

图1：本发明的流程图

图2：二维散射介质速度场模型

图3：激励源及接收点布设

图4：算例1(a)线性最小二乘法的反演图像；(b)压缩感知方法的反演图像

图5：算例2(a)线性最小二乘法的反演图像；(b)压缩感知方法的反演图像

图6：算例3(a)线性最小二乘法的反演图像；(b)压缩感知方法的反演图像

图7：算例4(a)线性最小二乘法的反演图像；(b)压缩感知方法的反演图像

图8：算例5(a)线性最小二乘法的反演图像；(b)压缩感知方法的反演图像

具体实施方式

现结合实施例、附图对本发明作进一步描述：

本发明主要采用仿真实验的方式验证该方法的可行性，所有步骤都经过实验验证，为实现基于压缩感知的尾波干涉成像，结合实施例，具体实施步骤如下：

(1)本数值仿真使用的激励源为15Hz主频的雷克子波，在10km×10km的非均匀散射介质中心处激发，持续时间为6π/100s，时间采样间隔Δt为0.001s。二维散射介质速度场模型如图2所示，通过传播速度的不同模拟介质的非均匀散射特性，将介质的速度场划分为200×200网格，产生均值为5000m/s，标准差为500m/s的随机数矩阵作为速度矩阵，其中最大速度为7320.95m/s，最小速度为3045.29m/s。接收点的布设原则是确保均匀覆盖介质区域，本仿真中采用的36个接收点和1个激励源的位置分布如图3所示。

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

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

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

为了模拟相对无限大的传播空间，本仿真的速度场边界设置为吸收边界，采用Clayton_Engquist_majda吸收边界算法，反射率约为2.5％，满足本仿真的要求。

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

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

当t_s的某取值t_smax使得R(t_s)取最大值时，该走时差t_smax即为所分析尾波部分扰动前后的时间延迟δt。那么根据尾波干涉理论的阐述，可得两列信号的相对波速变化

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

本仿真实验选择尾波观察时段为1.5～4.7s，每段窗口长度为0.5s，重叠0.2s，每个接收点可以获取扰动前后各10段尾波数据，36对收发对即可计算得到360个时延数据，将其存储为360×1的列向量形式。

(2)引入敏感核K(s₁,s₂,x₀,t)来描述波列在位置s₁发射，途经位置x₀，并在总传播时间t后在位置s₂被接收的概率大小，此概率描述了尾波在位置x₀处传播时间的空间密度，

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

式中，|s₁-s₂|表示s₁和s₂之间的距离；D为介质的散射系数，与其物理性质有关。通过激励源与接收点的位置和尾波的传播时间，对介质空间逐点计算K(s₁,s₂,x₀,t)，即可得到敏感核的空间分布。本仿真实验敏感核的构造中散射系数D＝8×10⁴m²/s，空间剖分为20×20的离散网格，网格边长为500m。

将步骤(1)中获取的时延数据通过敏感核项K(s₁,s₂,x₀,t)与相应的局部变化联系起来，

式中

为介质空间上x₀处的相对波速变化，即为需要反演计算的量；δt(t)为x₀处相对速度扰动/>

引起的走时变化，在物理意义上等效于尾波干涉求得的介质扰动前后同一接收点波形的时延。

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

(3)在压缩感知理论框架下求解(1)式。引入稀疏变换矩阵P对V进行离散稀疏变换。从而有(5)式T＝G^CSV^CS，其中，G^CS＝GP^-1，V^CS＝PV，本仿真采用的稀疏变换矩阵P为离散傅里叶变换矩阵。

(4)将式(5)转化为求解最优化问题式(6)

本仿真利用正交匹配追踪算法(OMP)求解出式(5)的稀疏解V^CS。

(5)通过稀疏逆变换求得原始解V，即式(7)V＝P^-1V^CS，整个过程无须使用P而仅须用到其逆矩阵P^-1，对于傅里叶变换矩阵，P^-1＝P^T。

得到的解V为一个400×1的列向量，用标准的列向量转化为矩阵的方法将其转化为20×20的矩阵形式，并通过MATLAB软件中的surf函数产生二维图像，此图像即为介质在扰动前后相对速度变化的二维空间分布情况。

结果对比：基于以上仿真环境及步骤，图4至8给出针对五个速度扰动算例的线性最小二乘差分成像方法和压缩感知成像方法的成像结果对比，其中，蓝色方框和黄色方框分别表示在该区域内施加+0.5％和-0.5％的速度变化。

两种方法的成像时间如表1所示。

表1线性最小二乘法与压缩感知法计算成像时间对比

结果分析：采用本发明提供的压缩感知方法进行尾波干涉成像，与线性最小二乘差分法相比，在单扰动和多扰动区域情况下都能更好地确定速度扰动区域的空间位置和范围，且避免了复杂的参数确定操作，执行简易，同时在保证反演精度的前提下能够明显减少计算时间，因此采用本发明提供的方法进行尾波干涉成像是可行的。

Claims (2)

1.一种基于压缩感知理论的尾波干涉成像方法，其特征在于步骤如下：

步骤1：通过尾波干涉技术获取时延数据，并将其存储为列向量；

步骤2：根据观测介质的物理性质构造敏感核矩阵，通过敏感核矩阵将尾波干涉法获取的时延数据与相应的局部变化联系起来，建立所示矩阵形式的反问题：

T＝GV

式中，T∈R^M,1表示通过尾波干涉技术获取的M个时延数据；V∈R^N,1表示空间各网格点对应的相对波速扰动值，N为网格总数；G∈R^M,N表示T到V的映射，其矩阵元素G_ij＝ΔsK_ij，其中K_ij为敏感核矩阵K的元素，K∈R^M,N的每行代表一个敏感核，每列代表一个网格，每个网格的面积为Δs；

其中根据观测介质的物理性质构造敏感核矩阵，通过敏感核矩阵将尾波干涉法获取的时延数据与相应的局部变化联系起来的具体过程为：

引入敏感核K(s₁,s₂,x₀,t)来描述波列在位置s₁发射，途经位置x₀，并在总传播时间t后在位置s₂被接收的概率大小，此概率描述了尾波在位置x₀处传播时间的空间密度

式中的p(s₁,s₂,t)表示波列从s₁经过时间t到达s₂的概率；

将步骤1中获取的时延数据通过K(s₁,s₂,x₀,t)与相应的局部变化联系起来：

式中

为介质空间上x₀处的相对波速变化，δt(t)为x₀处相对速度扰动/>

引起的走时变化；

步骤3：将G、T和V分别视为观测矩阵、测量值矩阵和真实信号，引入稀疏变换矩阵P对V进行离散稀疏变换：

T＝G^UV^U

式中，G^U＝GP^-1，V^U＝PV；

步骤4：将T＝G^UV^U转化为最优化问题，压缩感知重构算法求解出稀疏解V^U：

步骤5：通过稀疏逆变换求得原始解V为：

V＝P^-1V^U

将二维散射介质沿其长和宽剖分为n×n＝N个网格，每个独立网格空间内的相对波速变化视为一致；将求得的N维列向量展开为n×n的矩阵，n×n＝N，并将该矩阵导入MATLAB软件中的surf函数产生二维图像，此图像即为介质在扰动前后相对速度变化的二维空间分布情况。

2.根据权利要求1所述基于压缩感知理论的尾波干涉成像方法，其特征在于：所述稀疏变换矩阵为：离散傅里叶变换矩阵、离散余弦变换矩阵或离散小波变换矩阵。

Images