CN110348169A - 一种基于压缩感知理论的尾波干涉成像方法 - Google Patents
一种基于压缩感知理论的尾波干涉成像方法 Download PDFInfo
- Publication number
- CN110348169A CN110348169A CN201910674148.6A CN201910674148A CN110348169A CN 110348169 A CN110348169 A CN 110348169A CN 201910674148 A CN201910674148 A CN 201910674148A CN 110348169 A CN110348169 A CN 110348169A
- Authority
- CN
- China
- Prior art keywords
- matrix
- wave interference
- coda wave
- compressive sensing
- disturbance
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)
- Image Processing (AREA)
Abstract
本发明涉及一种基于压缩感知理论的尾波干涉成像方法,由于现有的线性最小二乘差分反演方法存在求解精度低、参数选取困难、计算效率低、难以定位多个扰动点的不足,同时考虑到扰动的空间区域相比于整体介质区域是稀疏的,恰好满足压缩感知理论对于信号稀疏性的要求。本发明所提供的基于压缩感知理论的尾波干涉成像方法,不需要进行复杂的参数选取操作,简单易行,计算效率高,求解结果稳定精确,并且对于同时存在多个扰动区域的情况下,仍能较为准确地获取速度扰动的空间位置和范围,相比于现有的反演方法,在实际的工程应用情况中更为适用。
Description
技术领域
本发明属于尾波干涉成像方法,涉及一种基于压缩感知理论的尾波干涉成像方法,涉及如何对尾波干涉成像的反问题进行精确快速求解的问题。
背景技术
尾波干涉成像是利用尾波时延和扩散近似敏感核来反演散射介质中微小速度扰动空间分布的技术。该技术本质上是对一个矩阵形式的反问题进行求解,该问题是一个欠定问题,有无穷多解,同时具有病态性,难以精确求解。目前,对于该反问题的求解,国内外研究中均采用地球物理反演中常用的线性最小二乘差分反演方法,主要步骤如下:
假设通过尾波干涉技术和扩散近似敏感核理论得到反演模型的线性方程组,并写成矩阵形式如下
T=GV (1)
式中,T∈RM,1表示通过尾波干涉技术获取的M个时延数据;V∈RN,1表示空间各网格点对应的相对波速扰动值,N为网格总数;G∈RM,N表示T到V的映射,其矩阵元素Gij=ΔsKij,其中Kij为敏感核矩阵K的元素,K∈RM,N的每行代表一个敏感核,每列代表一个网格,每个网格的面积为Δs。
利用线性最小二乘差分法对V值进行反演,即
V=V0+CmGT(GCmGT+Cd)-1(T-GV0) (2)
式中,V0为先验初始值,一般为零矩阵;GT表示G的转置;Cd为测得数据的对角协方差矩阵;Cm为介质空间元素平滑矩阵,可用下式计算:
式中,σm为先验标准差;λ0为网格间距;λ为相关长度。σm和λ一般通过L曲线法选取。s1和s2分别为激励点和接收点位置。
在用线性最小二乘差分法反演计算时,还需要循环迭代Cm以寻找最优的V值,迭代公式如下:
该方法存在求解精度不高、参数选取困难、计算效率低的缺陷,并且难以定位多个扰动点,在实际工程中的应用十分困难。
现有的线性最小二乘差分反演方法存在求解精度低、参数选取困难、计算效率低、难以定位多个扰动点的不足。
发明内容
要解决的技术问题
为了避免现有技术的不足之处,本发明提出一种基于压缩感知理论的尾波干涉成像方法,无论是存在单个还是多个扰动区域的情况下,该方法均能较准确地获取速度扰动的空间位置和范围,同时操作简易,具有较高的计算效率,能够更好地应用于工程实践中。
技术方案
一种基于压缩感知理论的尾波干涉成像方法,其特征在于步骤如下:
步骤1:通过尾波干涉技术获取时延数据,并将其存储为列向量;
步骤2:根据观测介质的物理性质构造敏感核矩阵,通过敏感核矩阵将尾波干涉法获取的时延数据与相应的局部变化联系起来,建立所示矩阵形式的反问题:
T=GV
式中,T∈RM,1表示通过尾波干涉技术获取的M个时延数据;V∈RN,1表示空间各网格点对应的相对波速扰动值,N为网格总数;G∈RM,N表示T到V的映射,其矩阵元素Gij=ΔsKij,其中Kij为敏感核矩阵K的元素,K∈RM,N的每行代表一个敏感核,每列代表一个网格,每个网格的面积为Δs;
步骤3:将G、T和V分别视为观测矩阵、测量值矩阵和真实信号,引入稀疏变换矩阵P对V进行离散稀疏变换:
T=GCSVCS, (5)
式中,GCS=GP-1,VCS=PV;
步骤4:将T=GCSVCS转化为最优化问题,压缩感知重构算法求解出稀疏解VCS:
步骤5:通过稀疏逆变换求得原始解V为:
V=P-1VCS, (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;pk(i,j)表示第k次(对应的时间)迭代时在(i,j)处的位移。
为了模拟相对无限大的传播空间,本仿真的速度场边界设置为吸收边界,采用Clayton_Engquist_majda吸收边界算法,反射率约为2.5%,满足本仿真的要求。
尾波干涉理论的实质是通过测量两个不同时刻尾波列的相位差来获取介质在该时段的波速变化。互相关法和延展法是目前基于此理论来提取介质波速变化的主要方法,考虑到仿真信号信噪比较高,本仿真实验采用计算效率更高的互相关法,步骤如下:假设已获取介质变化前后的两列尾波信号uunp和uper,通过两列波的互相关计算获得波速变化,即
式中,2T表示要分析的尾波部分的时间窗长度;t为时间窗中心位置;ts表示互相关中所用的走时差;R的值表示两列波的相似程度。
当ts的某取值tsmax使得R(ts)取最大值时,该走时差tsmax即为所分析尾波部分扰动前后的时间延迟δt。那么根据尾波干涉理论的阐述,可得两列信号的相对波速变化
式中δυ表示波速扰动,是整体介质波速变化的加权平均;υ为扰动前的介质波速。
本仿真实验选择尾波观察时段为1.5~4.7s,每段窗口长度为0.5s,重叠0.2s,每个接收点可以获取扰动前后各10段尾波数据,36对收发对即可计算得到360个时延数据,将其存储为360×1的列向量形式。
(2)引入敏感核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),即可得到敏感核的空间分布。本仿真实验敏感核的构造中散射系数D=8×104m2/s,空间剖分为20×20的离散网格,网格边长为500m。
将步骤(1)中获取的时延数据通过敏感核项K(s1,s2,x0,t)与相应的局部变化联系起来,
式中为介质空间上x0处的相对波速变化,即为需要反演计算的量;δt(t)为x0处相对速度扰动引起的走时变化,在物理意义上等效于尾波干涉求得的介质扰动前后同一接收点波形的时延。
为了使探测范围尽量覆盖整个介质,尾波干涉成像需要大量的敏感核,因此会有许多形如(14)式的方程,为清晰表达,可将线性方程组写成矩阵形式,即(1)式T=GV。
(3)在压缩感知理论框架下求解(1)式。引入稀疏变换矩阵P对V进行离散稀疏变换。从而有(5)式T=GCSVCS,其中,GCS=GP-1,VCS=PV,本仿真采用的稀疏变换矩阵P为离散傅里叶变换矩阵。
(4)将式(5)转化为求解最优化问题式(6)本仿真利用正交匹配追踪算法(OMP)求解出式(5)的稀疏解VCS。
(5)通过稀疏逆变换求得原始解V,即式(7)V=P-1VCS,整个过程无须使用P而仅须用到其逆矩阵P-1,对于傅里叶变换矩阵,P-1=PT。
得到的解V为一个400×1的列向量,用标准的列向量转化为矩阵的方法将其转化为20×20的矩阵形式,并通过MATLAB软件中的surf函数产生二维图像,此图像即为介质在扰动前后相对速度变化的二维空间分布情况。
结果对比:基于以上仿真环境及步骤,图4至8给出针对五个速度扰动算例的线性最小二乘差分成像方法和压缩感知成像方法的成像结果对比,其中,蓝色方框和黄色方框分别表示在该区域内施加+0.5%和-0.5%的速度变化。
两种方法的成像时间如表1所示。
表1线性最小二乘法与压缩感知法计算成像时间对比
结果分析:采用本发明提供的压缩感知方法进行尾波干涉成像,与线性最小二乘差分法相比,在单扰动和多扰动区域情况下都能更好地确定速度扰动区域的空间位置和范围,且避免了复杂的参数确定操作,执行简易,同时在保证反演精度的前提下能够明显减少计算时间,因此采用本发明提供的方法进行尾波干涉成像是可行的。
Claims (2)
1.一种基于压缩感知理论的尾波干涉成像方法,其特征在于步骤如下:
步骤1:通过尾波干涉技术获取时延数据,并将其存储为列向量;
步骤2:根据观测介质的物理性质构造敏感核矩阵,通过敏感核矩阵将尾波干涉法获取的时延数据与相应的局部变化联系起来,建立所示矩阵形式的反问题:
T=GV
式中,T∈RM,1表示通过尾波干涉技术获取的M个时延数据;V∈RN,1表示空间各网格点对应的相对波速扰动值,N为网格总数;G∈RM,N表示T到V的映射,其矩阵元素Gij=ΔsKij,其中Kij为敏感核矩阵K的元素,K∈RM,N的每行代表一个敏感核,每列代表一个网格,每个网格的面积为Δs;
步骤3:将G、T和V分别视为观测矩阵、测量值矩阵和真实信号,引入稀疏变换矩阵P对V进行离散稀疏变换:
T=GCSVCS
式中,GCS=GP-1,VCS=PV;
步骤4:将T=GCSVCS转化为最优化问题,压缩感知重构算法求解出稀疏解VCS:
步骤5:通过稀疏逆变换求得原始解V为:
V=P-1VCS
将二维散射介质沿其长和宽剖分为n×n=N个网格,每个独立网格空间内的相对波速变化视为一致;将求得的N维列向量展开为n×n的矩阵,n×n=N,并将该矩阵导入MATLAB软件中的surf函数产生二维图像,此图像即为介质在扰动前后相对速度变化的二维空间分布情况。
2.根据权利要求1所述基于压缩感知理论的尾波干涉成像方法,其特征在于:所述稀疏变换矩阵为:离散傅里叶变换矩阵、离散余弦变换矩阵或离散小波变换矩阵。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910674148.6A CN110348169B (zh) | 2019-07-25 | 2019-07-25 | 一种基于压缩感知理论的尾波干涉成像方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910674148.6A CN110348169B (zh) | 2019-07-25 | 2019-07-25 | 一种基于压缩感知理论的尾波干涉成像方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN110348169A true CN110348169A (zh) | 2019-10-18 |
CN110348169B CN110348169B (zh) | 2023-03-28 |
Family
ID=68178939
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201910674148.6A Active CN110348169B (zh) | 2019-07-25 | 2019-07-25 | 一种基于压缩感知理论的尾波干涉成像方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN110348169B (zh) |
Cited By (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN111983668A (zh) * | 2020-08-18 | 2020-11-24 | 中国科学技术大学 | 一种获取地震参数估计的方法及系统 |
Citations (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20090196229A1 (en) * | 2008-01-29 | 2009-08-06 | Zukang Shen | ACKNAK and CQI Channel Mapping Schemes in Wireless Networks |
CN103698764A (zh) * | 2013-12-27 | 2014-04-02 | 中国科学院电子学研究所 | 一种稀疏采样条件下的干涉合成孔径雷达成像方法 |
US20150036021A1 (en) * | 2011-11-10 | 2015-02-05 | Centre National De La Recherche Scientifique - Cnrs | Multiple Scattering Medium For Compressive Imaging |
CN105447818A (zh) * | 2015-11-16 | 2016-03-30 | 华东交通大学 | 基于变密度频域稀疏测量的图像重构方法 |
CN106405548A (zh) * | 2016-08-23 | 2017-02-15 | 西安电子科技大学 | 基于多任务贝叶斯压缩感知的逆合成孔径雷达成像方法 |
CN108845316A (zh) * | 2018-06-04 | 2018-11-20 | 中国卫星海上测控部 | 一种基于压缩感知理论的雷达稀疏探测方法 |
-
2019
- 2019-07-25 CN CN201910674148.6A patent/CN110348169B/zh active Active
Patent Citations (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20090196229A1 (en) * | 2008-01-29 | 2009-08-06 | Zukang Shen | ACKNAK and CQI Channel Mapping Schemes in Wireless Networks |
US20150036021A1 (en) * | 2011-11-10 | 2015-02-05 | Centre National De La Recherche Scientifique - Cnrs | Multiple Scattering Medium For Compressive Imaging |
CN103698764A (zh) * | 2013-12-27 | 2014-04-02 | 中国科学院电子学研究所 | 一种稀疏采样条件下的干涉合成孔径雷达成像方法 |
CN105447818A (zh) * | 2015-11-16 | 2016-03-30 | 华东交通大学 | 基于变密度频域稀疏测量的图像重构方法 |
CN106405548A (zh) * | 2016-08-23 | 2017-02-15 | 西安电子科技大学 | 基于多任务贝叶斯压缩感知的逆合成孔径雷达成像方法 |
CN108845316A (zh) * | 2018-06-04 | 2018-11-20 | 中国卫星海上测控部 | 一种基于压缩感知理论的雷达稀疏探测方法 |
Non-Patent Citations (2)
Title |
---|
ZHE BAI 等: "Low-Dimensional Approach for Reconstruction of Airfoil Data via Compressive Sensing", 《AIAA JOURNAL》 * |
朱路 等: "基于变密度稀疏采样的微波辐射干涉测量反演成像方法", 《计算机应用研究》 * |
Cited By (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN111983668A (zh) * | 2020-08-18 | 2020-11-24 | 中国科学技术大学 | 一种获取地震参数估计的方法及系统 |
Also Published As
Publication number | Publication date |
---|---|
CN110348169B (zh) | 2023-03-28 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN107479092B (zh) | 一种基于方向导数的频率域高阶声波方程正演模拟方法 | |
CN107247251A (zh) | 基于压缩感知的三维声源定位方法 | |
KR102021276B1 (ko) | 진폭 보존을 갖는 fwi 모델 도메인 각도 스택들 | |
CN108693403A (zh) | 一种广域虚拟密集化频谱态势生成方法 | |
CN109883532A (zh) | 一种声源识别与声场预报方法 | |
Huang et al. | Three-dimensional GPR ray tracing based on wavefront expansion with irregular cells | |
CN109141614B (zh) | 一种基于网络节点间水声通信信号的声速剖面反演方法 | |
Donelan et al. | A comparison of methods for estimating directional spectra of surface waves | |
Wang et al. | Review of underwater acoustic propagation models. | |
Maurer et al. | Optimized experimental design in the context of seismic full waveform inversion and seismic waveform imaging | |
CN111781639B (zh) | 针对obs多分量数据的炮检互易弹性波全波形反演方法 | |
CN106411438A (zh) | 浅水时变多途水声信道建模 | |
CN113866718B (zh) | 一种基于互质阵的匹配场被动定位方法 | |
CN112444773A (zh) | 基于空域融合的压缩感知二维doa估计方法 | |
CN109407152A (zh) | 基于零均值归一化互相关目标函数的时间域全波形反演方法 | |
CN105572734B (zh) | 一种以逆时偏移算法为引擎的波动方程初至走时层析方法 | |
Kuo et al. | Comparison of three different methods in investigating shallow shear-wave velocity structures in Ilan, Taiwan | |
CN110348169A (zh) | 一种基于压缩感知理论的尾波干涉成像方法 | |
Sun et al. | Estimation of multipath delay-Doppler parameters from moving LFM signals in shallow water | |
CN110146923A (zh) | 一种高效的高精度深度域地震子波提取方法 | |
Chen et al. | Investigating sound speed profile assimilation: An experiment in the Philippine Sea | |
Wapenaar et al. | On the retrieval of the directional scattering matrix from directional noise | |
Andrew et al. | Low-frequency pulse propagation over 510 km in the Philippine Sea: A comparison of observed and theoretical pulse spreading | |
US10317543B2 (en) | Estimation of a far field signature in a second direction from a far field signature in a first direction | |
CN109282841A (zh) | 一种无线无源声表面波传感器的超分辨率测量方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
CB03 | Change of inventor or designer information | ||
CB03 | Change of inventor or designer information |
Inventor after: Hou Hong Inventor after: Zhang Tao Inventor after: Zhu Jianghui Inventor before: Hou Hong Inventor before: Zhang Tao |
|
GR01 | Patent grant | ||
GR01 | Patent grant |