CN107390213B - 一种基于滑动时窗的探地雷达记录剖面的时延曲线提取方法 - Google Patents

Abstract

本发明公开了一种基于滑动时窗的探地雷达记录剖面的时延曲线提取方法，首先计算探地雷达记录剖面中的各孔径点回波信号的能量值，根据一维能量曲线确定各局部峰值点，确定探测目标个数和各目标的中心位置。然后提取某个目标中心位置处回波的时延值，以该目标中心位置处回波的时延值为起始点，设置滑动时窗的起止位置和长度，对中心点左右两侧的孔径点处的邻道数据截取相应的滑动时窗长度内的部分数据，仅对该段数据进行时延估计。遍历整个孔径点，得到该目标在整个孔径点上的时延曲线。再遍历各目标，得到各目标在整个孔径点上的时延曲线提取结果。本方法提取的时延曲线精度更高，根据时延曲线计算得到的成像结果更精确。

Description

一种基于滑动时窗的探地雷达记录剖面的时延曲线提取方法

技术领域

本发明属于探地雷达探测与应用技术领域，具体涉及探地雷达记录剖面中的目标时延曲线提取方法。

背景技术

探地雷达(Ground Penetrating Radar,GPR)是一种有效的浅层隐藏目标探测技术，利用电磁波在媒质电磁特性不连续处产生的反射和散射实现非金属覆盖区域中目标的成像探测。GPR回波的幅度和时延包含目标位置及电磁散射特性等信息。GPR沿一维测线进行空间扫描，在测线的每个孔径点处向地下区域发射电磁波并接收散射回波。每个孔径点处接收一道回波数据，多个孔径点接收的回波数据按列排列，就形成了GPR记录剖面。为进行高精度的层厚估计和地下目标成像，需要对GPR记录剖面中的目标时延曲线进行精确估计。

近些年来，国内外学者提出多种时延估计算法如子空间法，反卷积方法以及压缩感知等算法来处理GPR回波数据，子空间方法包括子空间旋转不变法，多重信号分类算法，最小范数法等，子空间方法应用于GPR地下回波时延估计时需要涉及强相干情况下目标回波的协方差矩阵计算。这类方法能获得高分辨率或超分辨率时延估计，但是不能直接来处理相干信号，所以需要与一些去相关算法进行结合【参考文献：L.Qu,Q.Sun,T.Yang,L.Zhang,Y Sun.Time-delay estimation for ground penetrating radar using ESPRITwith improved spatial smoothing technique.transactions on Geoscience andRemote Sensing,2014；11(8):1315-1319】。机器学习中如神经网络算法，支持向量机回归机算法被应用于GPR的时延估计上，这些方法可以直接处理相干信号，估计精度高。但机器学习算法需要建模与大量的训练序列的支持，计算量较大，运算时间慢【参考文献：Baltzard C L,Wang Y，et al.Time Delay and Permittivity estimation by ground-penetrating radar with support vector regression,IEEE Geoscience and RemoteSensing Letters,2014,11(4):873-877】。压缩感知与稀疏重构的方法近几年已经被应用于雷达领域中参数估计方面，而时延估计是参数估计的一种形式。压缩感知是一个针对信号采样的技术，它通过一些手段，实现了在采样过程中完成了数据压缩的过程。压缩感知的前提条件是信号具有稀疏性或者可压缩性，因此需要利用信号稀疏分解中已有的重构方法在概率意义上实现信号的精确重构或者一定误差下的近似重构。结果表明压缩感知方法能取得较高精确度，并且可直接处理相干信号与重叠回波的情况【参考文献：Jianzhong Li,Gang Wei,Cedric Le Bastard,Yide Wang,Biyun Ma,Meng Sun.Enhanced GPR Signalfor Layered Media Time-Delay Estimation in Low-SNR Scenario,IEEE Geoscienceand Remote Sensing Letters，2016；13(3)：299-303】。

上述方法主要针对的是单孔径点采样回波下的时延估计，并没有考虑GPR记录剖面的时延曲线特征和时延值提取方法。针对多孔径点采样回波，有必要设计一种GPR记录剖面的时延曲线提取方法。

发明内容

本发明所要解决的技术问题是提供一种基于滑动时窗的GPR记录剖面的时延曲线提取方法，该方法能提高GPR时延曲线的提取精度，提高成像效果。以下对各步骤进行具体说明：

一种基于滑动时窗的探地雷达记录剖面的时延曲线提取方法，在探地雷达记录剖面的时延估计中，利用了相邻道信号的相关性，以一种滑动窗口机制来确定邻道数据的预估回波区间，再将截取的信号进行时延估计；该时延曲线提取方法包括以下步骤：

步骤1：计算探地雷达记录剖面的一维能量曲线；根据一维能量曲线上局部峰值点即局部极大值点的个数和各所述局部峰值点的横向位置，确定被测区域中探测目标的个数和各探测目标的中心位置；

步骤2：将各探测目标中心位置对应的一维散射回波信号作为各探测目标的起始估计信号，对其进行时延估计，作为各个探测目标的初始时延估计值；

对一维回波信号进行时延估计为成熟的现有技术，可采用基于压缩感知和稀疏重建的时延估计方法；

对于每个探测目标，以其初始时延估计值为起始点，设置滑动时窗的起止位置和长度，利用滑动时窗对探测目标中心位置左右两侧孔径点处的邻道回波信号进行截取，并对截取的信号进行时延估计；从探测目标中心位置依次向外扩展，分别设置各孔径点处回波信号的滑动时窗的起止位置，进行该孔径点处回波的时延估计；遍历各孔径点，得到该探测目标在各孔径点的时延矢量形成该探测目标的时延曲线；遍历各探测目标，做如上相同处理，得到探地雷达记录剖面的时延曲线。

进一步地，所述步骤1中，计算探地雷达记录剖面的一维能量曲线的步骤为：

记探地雷达沿测线一维扫描得到的二维记录剖面数据为E(x,t)，E(x,t)＝[e(1,t),e(2,t),…,e(N,t)]，E(x,t)中x表示方向维，x＝1,…,N，N表示沿测线方向进行采样的孔径点数；t表示时间维，t＝1,…,M；M表示某孔径点处采集到的一维散射回波信号经数字采样后的采样数点；总采样时长为T ns，采样时间间隔为

对E(x,t)中的每一道数据进行能量值计算：

将第n个孔径点处采集到的一道数据记为e(n,t)，e(n,t)为一维时间信号；通过以下公式计算其能量值：

其中n＝1,…,N；

由P(1)～P(N)形成一维能量曲线图，该一维能量曲线图表示各孔径点处的散射回波信号的能量值随孔径点的变化情况。

进一步地，所述步骤1中，提取一维能量曲线的一阶导数从上到下穿越0的点，作为局部峰值点；设一共有S个局部峰值点，则设定被测区域中探测目标的个数为S，各探测目标的中心位置即为局部峰值点的横向位置，其对应的孔径点序号记为C₁,C₂,…,C_S。

所述步骤2中，设各探测目标中心位置对应的一维散射回波信号分别为第C₁,C₂,…,C_S道信号，即e(C₁,t),e(C₂,t),...,e(C_s,t)，得到的初始时延估计值为τ₁(C₁),τ₂(C₂),…,τ_S(C_S)。

进一步地，所述步骤2中，滑动时窗的长度为：

W＝N_c·ΔT 公式(2)

其中，N_c的取值方法为：

式中ceil(·)表示向大数方向取整，T_w表示回波能量最强的孔径点数据的子波时长，由探地雷达发射信号决定，为已知量。每个滑动窗口的长度是固定的。

根据探测场景的预估厚度和时间窗口长度等先验信息，设置时延估计值的预估范围为(a,b)；

对第s个探测目标的起始估计信号右侧的某道回波信号进行截取时，滑动时窗的起始值t_begin与终止值t_end设置为：

其中，τ_s(n-1)表示该道回波信号上一道信号，即其左侧邻道信号的时延估计值，n＝C_s+1,…,N，s＝1,2,…,S；

对第s个探测目标的起始估计信号左侧的某道回波信号进行截取时，滑动时窗的起始值t_begin与终止值t_end设置为：

其中，τ_s(n+1)表示该道回波信号上一道信号，即其右侧邻道信号的时延估计值，n＝1,L,C_s-1；

当t_end＞b时，设定其中b为时延估计值的预估范围的上限值，为经验取值。

进一步地，基于GPR孔径扫描的连续性和发射/接收天线的宽波束角，某一道时延估计的结果不会偏离上一道估计的时延值很远，因此在对截取的(滑动时窗内的)回波信号进行时延估计后，需要对估计得到的时延点进行判断筛选，得到最终的该道估计结果。具体步骤如下：

探地雷达向地面发射高频电磁波，会在接收端收到各种经过不同界面反射，且具有不同时延的回波信号，这些所有回波信号的总和构成了接收端的信号模型。所以接收端的信号模型由一系列的不同时延不同幅度的回波叠加，基于这个信号模型的求解，时延估计与幅度值应该是可以得到反射序列的幅值估计结果和时延估计结果。

所述步骤3中，已知序号为s的探测目标第n道信号的时延估计结果为τ_s(n)，s＝1,2,...S，采用基于稀疏重建的时延估计方法对下一道信号进行时延估计后，记反射序列的幅值估计结果为对应的时延估计结果

设下一道信号的序号为j；若下一道信号为第n道信号左侧的邻道信号，则j＝n-1,n≥2；若下一道信号为第n道信号右侧的邻道信号，则j＝n+1,n≤N-1；

根据下面所述公式对中的时延点进行判断筛选：

|τ_k-τ_s(n)|≤α 公式(5)

根据探测场景的预估厚度和时间窗口长度等先验信息，设计时延预估范围为τ_i∈(a,b)，在该范围上等距离设置N₀个时延点(包括点a、b以及它们中间的N₀-2个点)，N₀的取值由时延估计算法确定，则表示等分时延点间隔；

基于压缩感知和稀疏重建的时延估计方法【参考文献：Jianzhong Li,Gang Wei,Cedric Le Bastard,Yide Wang,Biyun Ma,Meng Sun.Enhanced GPR Signal for LayeredMedia Time-Delay Estimation in Low-SNR Scenario,IEEE Geoscience and RemoteSensing Letters，2016；13(3)：299-303”】首先要根据先验信息建立一个时延预估范围(a,b)。然后在该范围上，等距离划分间隔N₀时延点，第1个时延点τ₁为a，第N₀个时延点为中间的N₀-2个时延点计算公式为：τ_i＝a+(i-1)·α(1＜i＜N₀)；用τ_i表示重构后划分的各个时延点；用τ_i来构建稀疏矩阵，实现稀疏重构求解时延。

根据公式(5)得到筛选后的时延点τ_select及其对应幅值r_select，根据筛选后的时延点τ_select及其对应幅值r_select确定第j道信号的时延估计值τ_s(j)：

1)当筛选后剩下的时延点为1个的时候，令τ_s(j)＝τ_select；

2)当筛选后剩下的时延点为多个的时候，选取其中对应的幅值最大的时延点τ_m，令τ_s(j)＝τ_m；

3)当筛选后没有剩下时延点时，说明该道信号的时延估计结果与上一道信号的时延估计结果存在较大偏离，没有符合条件的时延点，将该道信号的时延估计结果筛除掉，设置该道信号的时延估计值为空，再下一道信号的滑动时窗的起始值和终止值采用该道信号滑动时窗的起始值和终止值。

有益效果

本发明提出一种基于滑动时窗的GPR记录剖面的时延曲线提取方法，根据一维能量曲线确定地下区域的探测目标个数和各探测目标在探地雷达(GPR)记录剖面的横向维的中心位置，利用了GPR记录剖面邻道信号之间的相关性特征，以探测目标的中心位置的时延值为起始点，依次向外扩展，分别设置各孔径点处回波的滑动时窗的起止位置，利用滑动时窗对回波进行截取并仅对截取的数据进行时延估计。基于GPR孔径扫描的连续性和发射/接收天线的宽波束角特性，某一道时延估计的结果不会偏离邻道信号时延值很远，以此对各道时延估计值进行筛选，去除不符合先验条件的时延点，保证了时延曲线随着孔径点偏离目标的横向中心位置有序向下滑动。本方法求解无需人工干预，有效抑制了非子波段的干扰，提取的时延曲线精度更高，成像更清晰。

附图说明

图1示出了本发明的算法流程图；

图2示出了GPR探测的正演模型；

图3示出了GPR记录剖面预处理后的数据；

图4示出了传统方法对GPR记录剖面的时延估计曲线结果；

图5示出了图3的GPR记录剖面数据的一维能量曲线图；

图6示出了本发明的GPR记录剖面的时延估计曲线结果；

图7示出了图5中的时延估计曲线经成像处理后的成像结果；

图8示出了图6中的时延估计曲线经成像处理后的成像结果。

具体实施方式

以下将结合附图和具体实施例对本发明做进一步详细说明。本发明公开了一种基于滑动时窗的GPR时延估计方法，如图1所示，首先计算探地雷达记录剖面中的各孔径点回波信号的能量值，根据一维能量曲线确定各局部峰值点，确定探测目标的个数和各探测目标在探地雷达(GPR)记录剖面的横向维的中心位置。然后提取某个探测目标的中心位置的时延值，以该探测目标的中心位置的时延值为起始点，设置滑动时窗的终止位置和长度，对中心点左右的邻道数据截取相应的滑动时窗长度内的部分数据，仅对该段数据进行时延估计。依次向外扩展，分别设置各孔径点处回波的滑动时窗的起止位置和长度，进行该孔径点处回波的时延估计。遍历整个孔径点，得到该目标在整个孔径点上的时延曲线。再遍历各目标，得到各目标在整个孔径点上的时延曲线提取结果。

实施例1：

以探测目标的中心位置对应的孔径点序号为C₁，即起始估计信号为第C₁道信号为例，得到初始时延估计值为τ₁(C₁)；向右选择第C₁+1道作为下一道估计信号，则把τ₁(C₁)作为第C₁+1道信号的滑动窗口的起始值t_begin，即第C₁+1道信号的滑动窗口的起始值t_begin与终值t_end为：

滑动窗口的终值t_end存在特殊情况：当t_end＞b即超过预估时间范围的时候。因此对滑动窗口终值t_end做如下处理：t_begin＝b-W，t_end＝b。每一道信号的滑动时窗是变化的，随着每道数据的时延估计值向邻道滑动。

基于滑动时窗对第C₁+1道信号进行加窗截取，截取t_begin-t_end范围的信号段，然后对截取信号段进行时延估计。

如果第C₁+1道信号的时延估计值τ₁(C₁+1)为空，第C₁+2道信号窗口起始值t_begin与终值t_end沿用第C₁+1道的t_begin与t_end值。对第C_s+2,…,N道信号依次做相同处理。向左选择第C₁-1道作为下一道估计信号，滑动时窗位置和长度的处理方式与右侧信道数据处理是相同的。

遍历各孔径点，得到该探测目标在各孔径点的时延矢量形成该目标的时延曲线；遍历各探测目标，做如上相同处理，得到探地雷达记录剖面的时延曲线。

实施例2：

本实例中，用GPRMax2D软件进行正演模拟，正演模拟模型如图2。模拟区域为1.5m×0.4m的地下区域，背景介质的相对介电常数ε₁＝9.0，电导率σ₁＝0.1S/m。目标为2个细长金属棒和1个细长介质棒，半径为0.01m。2个细长金属棒分别埋在(0.7,0.3)和(1.1,0.08)的位置，1个细长介质棒埋在(0.3,0.15)的位置处。发射天线在距地表上方0.005m处从左到右扫描探测，接收天线的高度与发射天线的高度相同，发收天线的间距为0.005m。发收一体天线沿着距地表高度为0.005m的直线在地表上方对地下区域进行合成孔径扫描，整个扫描孔径从0到1.5m，各孔径间距为0.015m，一共有100个孔径点。在每个孔径点处，发射天线向下发射电磁波，接收天线接收从地下散射的回波信号。随着收发一体天线从左到右移动扫描，一共可以获得100道散射数据，经预处理后的数据如图3所示。图4示出了传统方法B-Scan(横向扫描)时延估计曲线图，传统的B-Scan时延估计是各道数据各自估计时延，没有利用邻道信号的互相关信息。下面采用本发明的时延估计算法进行GPR B-Scan的时延曲线提取。图5示出了图3中的GPR B-Scan的一维能量曲线。从图中可见，存在3个局部峰值点，设定该GPR B-Scan中探测目标的个数为3。采用图1的处理流程，分别对这3个探测目标进行滑动时窗估计。探测目标1的中心位置位于序号为20的孔径点处，横向距离为0.3m，估计时延，作为探测目标1的初始时延估计结果；然后从探测目标的中心位置开始，设置滑动时窗的起止位置，并依次对左右两侧邻道数据，提取时窗内的数据，进行时延估计后并进行时延点的筛选；遍历100个孔径点，得到探测目标1在整个孔径点上的时延曲线。探测目标2，探测目标3的中心位置分别位于序号为49,77的孔径点处，对应横向距离为0.735m，1.155m。遍历这3个目标，得到时延提取结果，如图6所示。对比图4和图6，可以看出，本发明的时延曲线提取结果比传统方法的时延曲线提取结果更为精确。采用GPR加窗加权BP算法分别对图4和图6的时延曲线提取结果进行成像处理，结果如图7和图8所示。从成像结果来看，本发明的时延曲线提取后进行成像得到的成像结果比传统方法的时延曲线提取后进行成像得到的成像结果聚焦度更高，成像结果更精确。

对二维成像结果O(x_q,z_p),0≤p≤M-1,0≤q≤L-1，其中，M和L分别为横向和纵向采样点数，O(x_q,z_p)表示点(x_q,z_p)处的能量值；定义聚焦度因子：

该聚焦因子Ω可用来表示二维图像的聚焦程度。理想情况下，一个点目标的聚焦度因子为1。图像聚焦效果越差，该值越小。用上式分别对传统方法的时延曲线提取后进行成像得到的成像结果和本发明的时延曲线提取后进行成像得到的成像结果进行计算，得到二者的聚焦度因子分别为Ω₁＝0.0147，Ω₂＝0.2178。与传统时延曲线方法相比，本发明得到的时延曲线提取结果经过成像处理后的聚焦度提高了13.8倍。

Claims (5)

1.一种基于滑动时窗的探地雷达记录剖面的时延曲线提取方法，其特征在于,包括以下步骤：

步骤1：计算探地雷达记录剖面的一维能量曲线；根据一维能量曲线上局部峰值点即局部极大值点的个数和各所述局部峰值点的横向位置，确定被测区域中探测目标的个数和各探测目标的中心位置；

步骤2：将各探测目标中心位置对应的一维散射回波信号作为各探测目标的起始估计信号，对其进行时延估计，作为各个探测目标的初始时延估计值；

步骤3：对于每个探测目标，以其初始时延估计值为起始点，设置滑动时窗的起止位置和长度，利用滑动时窗对探测目标中心位置左右两侧孔径点处的邻道回波信号进行截取，并对截取的信号进行时延估计；从探测目标中心位置依次向外扩展，分别设置各孔径点处回波信号的滑动时窗的起止位置，进行该孔径点处回波的时延估计；遍历各孔径点，得到该探测目标在各孔径点的时延矢量，形成该目标的时延曲线；遍历各探测目标，做如上相同处理，得到探地雷达记录剖面的时延曲线。

2.根据权利要求1所述的基于滑动时窗的探地雷达记录剖面的时延曲线提取方法，其特征在于，所述步骤1中，计算探地雷达记录剖面的一维能量曲线的步骤为：

记探地雷达沿测线一维扫描得到的二维记录剖面数据为E(x,t)，E(x,t)＝[e(1,t),e(2,t),…,e(N,t)]，E(x,t)中x表示方向维，x＝1,2,…,N，N表示沿测线方向进行采样的孔径点数；t表示时间维，t＝1,2,…,M；M表示某孔径点处采集到的一维散射回波信号经数字采样后的采样数点；总采样时长为T ns，采样时间间隔为

对E(x,t)中的每一道数据进行能量值计算：

将第n个孔径点处采集到的一道数据记为e(n,t)，e(n,t)为一维时间信号；通过以下公式计算其能量值：

其中n＝1,2,…,N；

由P(1)～P(N)形成一维能量曲线图，该一维能量曲线图表示各孔径点处的散射回波信号的能量值随孔径点的变化情况。

3.根据权利要求2所述的基于滑动时窗的探地雷达记录剖面的时延曲线提取方法，其特征在于，所述步骤1中，提取一维能量曲线的一阶导数从上到下穿越0的点，作为局部峰值点；设一共有S个局部峰值点，则设定被测区域中探测目标的个数为S，各探测目标的中心位置即为局部峰值点的横向位置，其对应的孔径点序号记为C₁,C₂,…,C_S。

4.根据权利要求3所述的基于滑动时窗的探地雷达记录剖面的时延曲线提取方法，其特征在于，所述步骤3中，滑动时窗的长度为：

W＝N_c·ΔT 公式(2)

其中，N_c的取值方法为：

式中ceil(·)表示向大数方向取整，T_w表示回波能量最强的孔径点数据的子波时长，由探地雷达发射信号决定，为已知量；

对第s个探测目标的起始估计信号右侧的某道回波信号进行截取时，滑动时窗的起始值t_begin与终止值t_end设置为：

其中，τ_s(n-1)表示该道回波信号上一道信号，即其左侧邻道信号的时延估计值，n＝C_s+1,…,N，s＝1,2,…,S；

对第s个探测目标的起始估计信号左侧的某道回波信号进行截取时，滑动时窗的起始值t_begin与终止值t_end设置为：

其中，τ_s(n+1)表示该道回波信号上一道信号，即其右侧邻道信号的时延估计值，n＝1,…,C_s-1；

当t_end＞b时，设定其中b为时延估计值的预估范围的上限值，为经验取值。

5.根据权利要求4所述的基于滑动时窗的探地雷达记录剖面的时延曲线提取方法，其特征在于，所述步骤3中，已知序号为s的探测目标第n道信号的时延估计结果为τ_s(n)，s＝1,2,...S；采用基于稀疏重建的时延估计方法对下一道信号进行时延估计后，记反射序列的幅值估计结果为对应的时延估计结果

设下一道信号的序号为j；若下一道信号为第n道信号左侧的邻道信号，则j＝n-1,n≥2；若下一道信号为第n道信号右侧的邻道信号，则j＝n+1,n≤N-1；

根据下面所述公式对中的时延点进行判断筛选：

|τ_k-τ_s(n)|≤α 公式(5)

其中，表示等分时延点间隔；(a,b)为根据先验信息设计时延预估范围，N₀为在该范围(a,b)上等距离设置的时延点的个数，N₀的取值由时延估计算法确定；

根据公式(5)得到筛选后的时延点τ_select及其对应幅值r_select，根据筛选后的时延点τ_select及其对应幅值r_select确定第j道信号的时延估计值τ_s(j)：

1)当筛选后剩下的时延点为1个的时候，令τ_s(j)＝τ_select；

2)当筛选后剩下的时延点为多个的时候，选取其中对应的幅值最大的时延点τ_m，令τ_s(j)＝τ_m；

3)当筛选后没有剩下时延点时，说明该道信号的时延估计结果与上一道信号的时延估计结果存在较大偏离，没有符合条件的时延点，将该道信号的时延估计结果筛除掉，设置该道信号的时延估计值为空，再下一道信号的滑动时窗的起始值和终止值采用该道信号滑动时窗的起始值和终止值。