CN105403883A - 一种探地雷达地下目标位置检测方法 - Google Patents

Abstract

本发明涉及一种探地雷达地下目标位置检测方法，包括以下步骤：将探地雷达B-Scan数据进行多尺度的振幅、相位和方向特征分量的提取；提取第一个尺度的三个特征分量，利用HSB通道表示方法结合成彩色图像；找出目标感兴趣区域；对目标感兴趣区域内的彩色图像进行灰度化，并对目标感兴趣区域的振幅图像进行边缘提取；进行目标双曲线的定位；根据目标双曲线实现目标的位置确定。该方法不需要进行数据的训练，容易进行在线检测；能够较完整保留目标信息，从而提升目标定位的精度，对于浅层探测而言,不易产生虚警,不易漏掉非金属管线目标；能够提升杂波抑制效果，有效区分较强的杂波和目标回波；算法较为简单，保证快速进行目标定位。

Description

一种探地雷达地下目标位置检测方法

技术领域

本发明涉及一种探地雷达地下目标位置检测方法，属于探地雷达探测领域。

背景技术

探地雷达是近几十年迅速发展起来的一种有效的浅层地下目标探测技术，它是一种非破坏性探测手段，具有探测速度快、分辨率高、操作方便灵活、探测成本低等诸多优点，已被广泛应用于地下目标，如空洞、管道、地雷等的探测及定位。

探地雷达探测的二维回波数据称为B-Scan数据，它是后续雷达信号处理、目标识别及解译的数据基础，探地雷达目标定位技术也要基于B-Scan数据。对实现目标准确定位影响最大的是探地雷达B-Scan数据中的“杂波”。探地雷达杂波可看作是除了目标回波以外的各种回波，通常包括天线直达波、地表回波、地下非均匀介质产生的回波、以及伪目标所产生的回波等等。探地雷达杂波使得对地下目标的准确探测变得困难，尤其对于浅层埋地目标，目标回波与地表回波相比是较弱的成分，并且目标回波与地表回波间的时延很小，目标回波易被地表强回波这类杂波所淹没。因此有效的抑制杂波方法是实现探地雷达目标准确定位的首要任务。

目前国内外常见的定位方法主要是基于B-Scan图像的双曲线提取，根据提取到的双曲线进行速度估计再计算目标深度。目前常用的地下目标定位方法主要有以下五种：1、基于神经网络对双曲线的提取，需要较多的数据进行训练，不易实现在线检测；2、釆用模糊聚类的模式识别方法，对于金属管线和非金属管线都可能存在的浅层探测而言，容易产生虚警，且容易漏掉非金属管线目标；3、基于图像分割和霍夫变换的方法，应用在浅层探测管线的时候，不能有效区分较强的杂波和目标回波；4、基于图像分割和模板匹配的方法应用在浅层探测管线时候，由于管径的大小可能多变，从而对应的模版也较多，导致算法运算时间较长；5、基于形态学的曲线检测，是根据图像的灰度值进行检测判断，能够判断目标的区域但是得到是多根曲线，进行下一步计算还需对曲线进行处理。

另外，单演小波(Monogenic)信号是一维解析信号二维延拓的结果，该方法是将实系数的小波基和它复值Riesz变换结合起来构成多尺度单演信号分析，从而使得每个小波系数都有振幅，相位和方向信息。

发明内容

本发明的目的是提供一种探地雷达地下目标位置检测方法，用以解决传统的地下目标定位方法存在着诸多弊端的问题。

为实现上述目的，本发明的方案包括一种探地雷达地下目标位置检测方法，所述位置检测方法包括以下步骤：

(1)、将探地雷达B-Scan数据进行多层分解提取，提取多尺度单演信号调制特征分量，包括振幅、相位和方向特征分量；

(2)、提取第一个尺度，即从多层分解中选择第一层的振幅、相位和方向特征分量，利用HSB通道表示方法结合成彩色图像；

(3)、从振幅特征分量中找出目标感兴趣区域；

(4)、对目标感兴趣区域内的彩色图像进行灰度化，并对所述目标感兴趣区域的振幅图像进行边缘提取；

(5)、进行目标双曲线的定位；

(6)、根据所述目标双曲线实现目标的位置确定。

所述步骤(1)具体为：

1)、构造小波域单演信号分析框架：假设f是原始信号，解析小波ψ_i(x)，得到的小波系数ω_i[k]：

${\mathrm{\&omega;}}_i\&lsqb;k\&rsqb;=<f,{\mathrm{\&psi;}}_{i,k}>=({\mathrm{\&psi;}}_i*f)\left(x\right){|}_{x=D^{-(i+1)}k},$

其中，r_1，i[k]是Riesz变换的实部，r_2，i[k]是Riesz变换的虚部；

2)、计算局部结构矩阵：

其中，n，m＝1，2；r₁[k]和r₂[k]分别是Riesz变换的实部和虚部；

基于张量的方法得到方向向量为：

$\mathrm{\&theta;}=\frac{1}{2}\arctan \left(\frac{2J_{12}}{J_{22}-J_{11}}\right),$ u＝(cosθ，sinθ)，

根据选定的方向，相应的方向Hilbert变换是：

q_i[k]＝r_1，i[k]cosθ+r_2，i[k]sinθ或者

3)、计算振幅特征分量和相位特征分量：

$A_i\&lsqb;k\&rsqb;=\sqrt{{\mathrm{\&omega;}}_i{\&lsqb;k\&rsqb;}^2+q_i{\&lsqb;k\&rsqb;}^2},{\mathrm{\&xi;}}_i\&lsqb;k\&rsqb;=\arctan \left(\frac{q_i\&lsqb;k\&rsqb;}{{\mathrm{\&omega;}}_i\&lsqb;k\&rsqb;}\right).$

所述步骤(2)中，所述利用HSB通道表示方法结合成彩色图像具体为：所述彩色图像中的色度分量用方向特征分量表示，饱和度分量用相位特征分量表示，亮度分量用相位特征分量表示。

所述步骤(3)具体为：

1)、针对每一道A-Scan数据，利用能量和方差两个统计量，选取能量和方差都为峰值的区域，并结合第一设定阈值，确定出目标所在的深度范围，用时间窗口表示为：t_a～t_b；

2)、让每道A-scan数据中除t_a～t_b时间窗口之外的数据全为0，然后根据A-scan数据的能量变化情况曲线，选择曲线中最大值处的位置，并结合第二设定阈值，确定目标的目标感兴趣区域。

所述步骤(4)具体为：

1)、对目标感兴趣区域内的彩色图像进行灰度化；

2)、采用高斯滤波进行去噪和剔除虚假目标；

3)、利用Canny算子将所述目标感兴趣区域的正负图像进行边缘提取转换为二值图像。

所述步骤(5)包括以下两个步骤：

1)、利用最小二乘法估计目标双曲线的顶点坐标；

2)、根据估计的目标双曲线的顶点坐标，采用霍夫变换来求取波速，进而得到目标双曲线的定位。

所述步骤(5)中的步骤1)具体为：设定所述二值图像中值为1的点为待拟合的点，优化条件是：

其中，(x₀,t₀)是估计系数的向量，x_i和t_i是二值图像中n个值为1的点数据行向量的横纵坐标值，i＝1,…,n，f(x₀,t₀,x)是的变形， $t=\frac{2}{v}\left(\sqrt{{\left(\frac{{\mathrm{vt}}_0}{2}\right)}^2+{(x-x_0)}^2}\right);$

然后采用对多元函数求极值的方法，对x₀求偏导，让偏导等于0；对t₀求偏导，让偏导等于0；进而得到关于参数(x₀,t₀)的方程，通过解所述关于参数(x₀,t₀)的方程得到拟合的曲线；

最后，取该拟合好的曲线的纵坐标的极大值点为(x₀,t₀)中的坐标x₀，同时得到坐标t₀；

所述步骤(5)中的步骤2)具体为：根据估计出的顶点坐标(x₀,t₀)，采用一维的霍夫变换求取波速v：根据探地区域设定传播速度v在一定范围内变换，在此范围内变换参数v的值，此时参数空间将有聚焦点，聚焦点的集中的值即为波速v。

所述目标的定位包括：

目标深度定位：H＝t₀·v，其中，t₀为双曲线顶点的反射回波时延，v为电磁波在介质中的传播速度，H为目标的深度位置；

目标水平位置定位：S＝x₀·Δx，x₀为双曲线顶点的横坐标，Δx为天线的移动步长，S为目标的水平位置。

首先，本发明提供的探地雷达地下目标位置检测方法，不需要进行数据的训练，所以容易进行在线检测；能够较完整保留目标信息，从而提升目标定位的精度，对于金属管线和非金属管线都可能存在的浅层探测而言,不易产生虚警,所以不易漏掉非金属管线目标；能够提升杂波抑制效果，有效区分较强的杂波和目标回波；基于的算法较为简单，算法运算时间短，保证了快速进行目标定位。

附图说明

图1是本发明提供的探地雷达地下目标位置检测方法的流程图；

图2是探地雷达实测的B-Scan回波图；

图3是彩色图像二值化和ROI算法限定区域后得到的图；

图4是采用雷达天线与目标B-Scan回波的几何位置关系图。

具体实施方式

下面结合附图对本发明做进一步详细的说明。

单演小波(Monogenic)信号是一维解析信号二维延拓的结果，该方法是将实系数的小波基和它复值Riesz变换结合起来构成多尺度单演信号分析，从而使得每个小波系数都有振幅，相位和方向信息，因此，采用单演小波进行目标回波振幅特征的提取，能够在较完整保留目标信息的同时提升杂波抑制效果，从而提升目标定位的精度。

本发明是一种基于多尺度单演(Monogenic)信号分析的方法来进行探地雷达地下目标的位置检测和确定，如图1所示，具体如下：

1、首先将图2所示的探地雷达回波B-Scan图像进行多层分解提取，提取多尺度单演信号调制特征分量，包括振幅、相位和方向特征分量。具体包含以下步骤：

1.1、针对探地雷达的B-Scan探测回波数据f(χ)，其中χ＝(x，y)，求其Riesz变换：

其中即表示f(χ)的Riesz变换，为Riesz变换的实部，为Riesz变换的虚部。

1.2、构造小波域单演信号分析框架。假设f是原始信号，解析小波ψ_i(x)，得到的小波系数ω_i[k]：

${\mathrm{\&omega;}}_i\&lsqb;k\&rsqb;=<f,{\mathrm{\&psi;}}_{i,k}>=({\mathrm{\&psi;}}_i*f)\left(x\right){|}_{x=D^{-(i+1)}k}---\left(2\right)$

其中，r_1，i[k]是ω_i[k]的Riesz变换的实部，r_2，i[k]是ω_i[k]的Riesz变换的虚部。

1.3、计算局部结构矩阵：

${\&lsqb;J\left(k\right)\&rsqb;}_{mn}=J_{mn}=\underset{l\&Element;Z^2}{\&Sigma;}v\&lsqb;l-k\&rsqb;r_m\&lsqb;l\&rsqb;r_n\&lsqb;l\&rsqb;,---\left(4\right)$

其中，n，m＝1，2；r₁[k]和r₂[k]分别是Riesz变换的实部和虚部；

基于张量的方法得到方向向量为：

$\mathrm{\&theta;}=\frac{1}{2}\arctan \left(\frac{2J_{12}}{J_{22}-J_{11}}\right),---\left(5\right)$

u＝(cosθ，sinθ)(6)

根据选定的方向，相应的方向Hilbert变换可以用下式表示：

q_i[k]＝r_1，i[k]cosθ+r_2，i[k]sinθ(7)

或者：

${\mathrm{\&omega;}}_i\&lsqb;k\&rsqb;+q_i\&lsqb;k\&rsqb;j=A_i\&lsqb;k\&rsqb;\&CenterDot;e^{{\mathrm{j\&xi;}}_i\&lsqb;k\&rsqb;}---\left(8\right)$

1.4、计算得到的振幅特征分量和相位特征分量：

$A_i\&lsqb;k\&rsqb;=\sqrt{{\mathrm{\&omega;}}_i{\&lsqb;k\&rsqb;}^2+q_i{\&lsqb;k\&rsqb;}^2},---\left(9\right)$

${\mathrm{\&xi;}}_i\&lsqb;k\&rsqb;=\arctan \left(\frac{q_i\&lsqb;k\&rsqb;}{{\mathrm{\&omega;}}_i\&lsqb;k\&rsqb;}\right)---\left(10\right)$

从而实现了将探地雷达B-Scan数据进行多层分解提取，在小波子频带中由小波振幅提取出多尺度的振幅、相位和方向特征分量。

2、从上述多尺度的信息中，即多层分解中，提取第一个尺度的振幅、相位和方向特征分量，即从多层分解中选择第一层的振幅、相位和方向特征分量，利用HSB通道表示方法结合成彩色图像，该彩色图像中，色度分量用方向特征分量表示，饱和度分量用相位特征分量表示，亮度分量用相位特征分量表示。

3、采用基于能量统计的目标感兴趣区域提取方法从上述振幅分量图像中提取目标区域。具体包括以下两个步骤：

3.1获取时间窗口：

探地雷达二维B-Scan回波图像是由多道一维A-Scan回波数据组成的。在步骤1中提取出的B-Scan图像振幅分量的基础上，针对每一道A-Scan数据，根据其在时间t方向的特点可知，在非目标区域A-Scan在t方向的方差较小，能量较小；在目标区域，A-scan的方差较大，能量较大。利用能量和方差两个统计量，选取能量和方差都为峰值的区域，并结合合适的阈值，确定出目标所在的深度范围，用时间表示t_a～t_b，即时间窗口。(其中，本段中的阀值可由人工预设或迭代计算方式获得，这是常规技术，这里不做赘述)

3.2加窗能量法提取ROI：

选择时间窗口后，对A-scan数据进行归一化处理，让每道A-scan中除t_a～t_b之外的数据全为0。加窗后，在测线方向上每道A-scan数据都减少了许多干扰信号，剩下的信号成分中如果无目标，则每个采样点的幅值接近于零；如有目标，则有较大幅值。由此，可以根据加窗后每个A-scan数据的能量变化情况，得到对应的能量变化情况曲线，选择能量曲线中最大值处的位置，并结合阈值的选取确定目标的ROI位置。(其中，本段中的阀值也可由人工预设或迭代计算方式获得，这是常规技术，这里不做赘述)

上述为本发明提供的一种ROI提取方法，现有技术中的常用的ROI提取方法有时频分析、神经网络、模板匹配法。时频分析方法对信噪比不高的区域，特别是有相关频率的干扰的区域容易产生误判；神经网络法需要有大的已知数据做训练，其计算速度也不能很好地适用于现场探测；模板匹配法需要对测量数据有很多的先验知识，因此限制了它的应用。

4、对目标感兴趣区域内的彩色图像进行灰度化，并利用Canny算子针对上述ROI区域的振幅图像进行边缘提取。具体包含以下三个步骤：

4.1、对目标感兴趣区域(ROI区域)内的彩色图像进行灰度化；

4.2、采用高斯滤波对ROI区域的振幅图像进行去噪和剔除虚假目标；

4.3、利用Canny算子将已通过高斯滤波处理后的ROI区域的振幅图像进行边缘提取转换为二值图像，如图3所示。

5、进行目标双曲线的定位，整体包括以下两个步骤：

5.1、利用最小二乘法估计目标双曲线的顶点坐标；

5.2、根据估计的目标双曲线的顶点坐标，采用霍夫(Hough)变换来求取波速，进而得到目标双曲线的定位。

由探地雷达原理知：

$\frac{t^2}{t_0^2}-\frac{4{(x-x_0)}^2}{v^2{t}_0^2}=1---\left(11\right)$

其中，x表示天线位置，x₀表示目标的水平位置，v表示电磁波在介质中的传播速度，t₀表示目标反射回波时延。因此，求出顶点坐标(x₀,t₀)以及波速v即可准确定位目标。

首先利用最小二乘方法来估计双曲线的顶点坐标(x₀,t₀)，这样做可以减少Hough变换估计的参数。

具体来说将上面步骤4.3中得到的二值图像中值为1的点为待拟合的点，这里的优化条件是：

其中(x₀,t₀)是估计系数的向量，x_i和t_i是二值化后的图像中n个值为1的点数据行向量的横纵坐标值，i＝1,…,n，这里f(x₀,t₀,x)是(11)式的变形， $t=\frac{2}{v}\left(\sqrt{{\left(\frac{{\mathrm{vt}}_0}{2}\right)}^2+{(x-x_0)}^2}\right).$

最小二乘拟合时采用对多元函数求极值的方法，对x₀求偏导，让偏导等于0；对t₀求偏导，让偏导等于0；得到关于参数(x₀,t₀)的方程，通过解此方程得到拟合的曲线。

拟合好曲线后，通过取该曲线纵坐标的极大值点得到顶点的坐标x₀，对应的坐标即为顶点的坐标t₀。

根据估计出的顶点坐标(x₀,t₀)，进一步采用一维的霍夫(Hough)变换求取波速v，并进行目标双曲线的定位：根据探地区域设定传播速度v在一定范围内变换，在此范围内变换参数(波速)v的值，此时参数空间将有聚焦点，聚焦点的值越大波速值越准确，所以取聚焦点的集中的值为波速v。

6、实现目标的定位，如图4所示，具体包含以下两个步骤：

6.1、目标深度定位：根据探地雷达原理知，双曲线顶点的纵坐标表示最短回波时延，即在这一测点探地雷达距离目标最近，因此双曲线的纵坐标就代表目标对应的深度。将提取到的双曲线顶点的反射回波时延乘以电磁波在介质中的传播速度v，就得到目标的深度位置H，即：

H＝t₀·v(10)

6.2、目标水平位置定位：同理，双曲线顶点的横坐标即代表目标对应的水平位置，因此，将提取到的双曲线顶点的横坐标乘以天线的移动步长Δx，就得到目标的水平位置S，即：

S＝x₀·Δx(11)

在顶点坐标(x₀,t₀)处，即在目标处，振幅的能量最大。

说明书中未详细说明的部分属于本领域技术人员公知的现有技术。

以上给出了具体的实施方式，但本发明不局限于所描述的实施方式。本发明的基本思路在于上述基本方案，对本领域普通技术人员而言，根据本发明的教导，设计出各种变形的模型、公式、参数并不需要花费创造性劳动。在不脱离本发明的原理和精神的情况下对实施方式进行的变化、修改、替换和变型仍落入本发明的保护范围内。

Claims (8)

1.一种探地雷达地下目标位置检测方法，其特征在于，所述位置检测方法包括以下步骤：

(1)、将探地雷达B-Scan数据进行多层分解提取，提取多尺度单演信号调制特征分量，包括振幅、相位和方向特征分量；

(2)、提取第一个尺度，即从多层分解中选择第一层的振幅、相位和方向特征分量，利用HSB通道表示方法结合成彩色图像；

(3)、从振幅特征分量中找出目标感兴趣区域；

(4)、对目标感兴趣区域内的彩色图像进行灰度化，并对所述目标感兴趣区域的振幅图像进行边缘提取；

(5)、进行目标双曲线的定位；

(6)、根据所述目标双曲线实现目标的位置确定。

2.根据权利要求1所述的探地雷达地下目标位置检测方法，其特征在于，所述步骤(1)具体为：

1)、构造小波域单演信号分析框架：假设f是原始信号，解析小波ψ_i(x)，得到的小波系数ω_i[k]：

${\mathrm{\&omega;}}_i\&lsqb;k\&rsqb;=<f,{\mathrm{\&psi;}}_{i,k}>=({\mathrm{\&psi;}}_i*f)\left(x\right){|}_{x=D^{-(i+1)}k},$

其中，r_1，i[k]是Riesz变换的实部，r_2，i[k]是Riesz变换的虚部；

2)、计算局部结构矩阵：

其中，n，m＝1，2；r₁[k]和r₂[k]分别是Riesz变换的实部和虚部；

基于张量的方法得到方向向量为：

$\mathrm{\&theta;}=\frac{1}{2}arctan\left(\frac{2J_{12}}{J_{22}-J_{11}}\right),$ u＝(cosθ，sinθ)，

根据选定的方向，相应的方向Hilbert变换是：

q_i[k]＝r_1，i[k]cosθ+r_2，i[k]sinθ或者 ${\mathrm{\&omega;}}_i\&lsqb;k\&rsqb;+q_i\&lsqb;k\&rsqb;j=A_i\&lsqb;k\&rsqb;\&CenterDot;e^{{\mathrm{j\&xi;}}_i\&lsqb;k\&rsqb;};$

3)、计算振幅特征分量和相位特征分量：

$A_i\&lsqb;k\&rsqb;=\sqrt{{\mathrm{\&omega;}}_i{\&lsqb;k\&rsqb;}^2+q_i{\&lsqb;k\&rsqb;}^2},{\mathrm{\&xi;}}_i\&lsqb;k\&rsqb;=arctan\left(\frac{q_i\&lsqb;k\&rsqb;}{{\mathrm{\&omega;}}_i\&lsqb;k\&rsqb;}\right).$

3.根据权利要求1所述的探地雷达地下目标位置检测方法，其特征在于，所述步骤(2)中，所述利用HSB通道表示方法结合成彩色图像具体为：所述彩色图像中的色度分量用方向特征分量表示，饱和度分量用相位特征分量表示，亮度分量用相位特征分量表示。

4.根据权利要求1所述的探地雷达地下目标位置检测方法，其特征在于，所述步骤(3)具体为：

1)、针对每一道A-Scan数据，利用能量和方差两个统计量，选取能量和方差都为峰值的区域，并结合第一设定阈值，确定出目标所在的深度范围，用时间窗口表示为：t_a～t_b；

2)、让每道A-scan数据中除t_a～t_b时间窗口之外的数据全为0，然后根据A-scan数据的能量变化情况曲线，选择曲线中最大值处的位置，并结合第二设定阈值，确定目标的目标感兴趣区域。

5.根据权利要求1所述的探地雷达地下目标位置检测方法，其特征在于，所述步骤(4)具体为：

1)、对目标感兴趣区域内的彩色图像进行灰度化；

2)、采用高斯滤波进行去噪和剔除虚假目标；

3)、利用Canny算子将所述目标感兴趣区域的正负图像进行边缘提取转换为二值图像。

6.根据权利要求5所述的探地雷达地下目标位置检测方法，其特征在于，所述步骤(5)包括以下两个步骤：

1)、利用最小二乘法估计目标双曲线的顶点坐标；

2)、根据估计的目标双曲线的顶点坐标，采用霍夫变换来求取波速，进而得到目标双曲线的定位。

7.根据权利要求6所述的探地雷达地下目标位置检测方法，其特征在于，所述步骤(5)中的步骤1)具体为：设定所述二值图像中值为1的点为待拟合的点，优化条件是：

其中，(x₀,t₀)是估计系数的向量，x_i和t_i是二值图像中n个值为1的点数据行向量的横纵坐标值，i＝1,…,n，f(x₀,t₀,x)是的变形， $t=\frac{2}{v}\left(\sqrt{{\left(\frac{{\mathrm{vt}}_0}{2}\right)}^2+{(x-x_0)}^2}\right);$

然后采用对多元函数求极值的方法，对x₀求偏导，让偏导等于0；对t₀求偏导，让偏导等于0；进而得到关于参数(x₀,t₀)的方程，通过解所述关于参数(x₀,t₀)的方程得到拟合的曲线；

最后，取该拟合好的曲线的纵坐标的极大值点为(x₀,t₀)中的坐标x₀，同时得到坐标t₀；

所述步骤(5)中的步骤2)具体为：根据估计出的顶点坐标(x₀,t₀)，采用一维的霍夫变换求取波速v：根据探地区域设定传播速度v在一定范围内变换，在此范围内变换参数v的值，此时参数空间将有聚焦点，聚焦点的集中的值即为波速v。

8.根据权利要求7所述的探地雷达地下目标位置检测方法，其特征在于，所述目标的定位包括：

目标深度定位：H＝t₀·v，其中，t₀为双曲线顶点的反射回波时延，v为电磁波在介质中的传播速度，H为目标的深度位置；

目标水平位置定位：S＝x₀·Δx，x₀为双曲线顶点的横坐标，Δx为天线的移动步长，S为目标的水平位置。