CN104020495B

CN104020495B - 一种基于探地雷达的地下管线参数自识别方法

Info

Publication number: CN104020495B
Application number: CN201410283348.6A
Authority: CN
Inventors: 郑晶; 杨峰; 杜翠; 乔旭; 夏云海
Original assignee: China University of Mining and Technology Beijing CUMTB
Current assignee: China University of Mining and Technology CUMT; China University of Mining and Technology Beijing CUMTB
Priority date: 2014-06-24
Filing date: 2014-06-24
Publication date: 2015-05-06
Anticipated expiration: 2034-06-24
Also published as: CN104020495A

Abstract

本发明专利提供一种基于探地雷达的地下管线参数自识别方法，该方法通过检测和处理由地下管线的回波信号，自动识别地下管线目标的位置、尺寸等信息，其实现步骤为：1.预处理，利用SVD对探地雷达接收的原始剖面图像进行变换分析，分离背景和目标信号；2.类双曲线边缘提取，对预处理图像进行边缘检测，提取管线反射形成的类双曲线同相轴的位置及形状；3.目标初始特征参数提取，对步骤2产生的类双曲线同相轴图像，进行管线空间位置、上覆地层介质平均速度和管径等参数进行提取，获得管径反演计算的初始参数；4.速度扫描及管线参数反演迭代；5.Hough变换提取管线参数。本发明是基于探地雷达探测技术，对地下管线目标体进行管径参数的自识别方法。

Description

一种基于探地雷达的地下管线参数自识别方法

技术领域

本发明属于探地雷达地下目标识别的技术领域，涉及一种基于探地雷达的地下管线参数的自识别方法，该方法通过处理由介质不连续处散射回来的高频脉冲电磁波，达到检测地下管线目标，并对其位置、尺寸等信息进行自识别的目的。

背景技术

地质雷达是通过发射高频脉冲电磁波进行地下目标探测的技术方法。电磁波的产生及其传播规律是进行地质雷达数据处理及解释的研究基础。电磁波是通过变化的电场与变化的磁场相互激发在介质中进行传播并形成电磁场。介质中（尤其是土壤中）电磁场、电磁波的传播——波速、衰减、反射与折射的理论是地质雷达探测检测的理论基础。

地质雷达是一种快速、高效、无损探测的物探方法。常用的雷达探测方法主要以反射式探测和透射式探测这两种。反射式雷达探测原理：反射探测是雷达的发射天线和接收天线都放置于被测介质表面，通过向介质中发射高频宽带电磁波脉冲信号，并接收介质中有差异地方所反射回的电磁波信号实现探测的过程。本发明专利就是基于该探测原理开展的。

发明内容

本发明的目的是对浅层地下管线目标实现自识别，提供目标的定位及管径等信息。

本发明提出了一种基于反射回波的浅层地下管线目标探测预处理技术。由于地表的强反射回波、设备天线间的耦合、媒质内电导率的不均匀及其中大尺寸颗粒的反射等杂波信号的存在，使得管径等目标反射信号不明显，因此对回波信号进行预处理，从中滤除背景噪声及各类杂波信号对于后续的目标检测和目标识别非常关键，预处理技术包括：零偏矫正、数字滤波、小波软阈值去噪、测线间去均值、SVD变换分离背景和目标信号等。

本发明还给出了一种类双曲线的检测技术，该技术采用小波模极大值方法（但不局限于此，还可用方向算子、梯度算子、canny算子及Sobel等方法进行）对预处理后雷达剖面进行边缘提取，在得到类双曲线同时，最大的保留了有用信号，实现了对不连续的双曲线进行自适应拾取、检测，并采用线连接方法将较近的类双曲线相连，得到连续的双曲线。

本发明还给出了一种基于波动反射理论的目标定位、管径尺寸及介质速度的计算方法，该方法利用已知的接收天线位置和到达时间等信息，计算出目标的位置、管径尺寸及介质速度等参数，作为初始参考值，计算方法参照下面的公式（2）。

本发明还给出了一种空间散射点真实位置的反演方法，由于读取误差的存在，利用公式(2)提取的参数与实际值存在一定误差，只能作为反演计算的初始值，对速度在一定范围内进行扫描，反演出各扫描参数对应的散射点位置，当参数值与真实值一致时，反演出的散射点位置组合成真实的管线外壁。

本发明还给出了一种基于Hough变换的圆检测技术，结合峰值检测技术，得到管线目标位置及管径等参数的最佳值：当参数值与真实值一致时，反演出的曲线形状最接近于真实的管线外壁，也最类似于圆，而其他情况则类似于椭圆或双曲线，以反演的结果作为Hough的输入时，当曲线越接近于圆时，参数域中的‘亮点’值越大，通过峰值检测技术，在速度扫描范围内，找到最大值，并将该值对应的速度值记录为介质速度的估计值，从而对应得到管线的目标位置及管径等参数。

本发明的特征在于，同时对预处理后的剖面进行边缘提取，在消除杂波的基础上，最大限度的保留了原始剖面中的目标双曲线边缘信号。

本发明的特征在于，给出了一种基于波动反射理论的目标定位、管径及介质速度的计算方法，结合多个已知接收信号的到达时间和接收天线位置进行运算，即可计算出目标的位置、管径的尺寸以及介质的速度等参数，无需已知介质的传输速度，计算方法见公式（2），由于并对参数在一定的范围内进行扫描计算，反演出一定参数范围内的类双曲线上信号对应空间真实散射点的位置范围，由于读取误差的存在，利用公式(2)提取的参数与实际值存在一定误差，只能作为反演计算的初始值。

本发明的特征在于，以基于波动反射理论的计算结果为参考值，给出了目标真实散射点的反演方法，计算方法见公式（3），对速度在一定范围内进行扫描时，反演出各扫描参数对应的散射点位置，当参数值与真实值一致时，反演出的散射点位置可以组合成真实的管线外壁。

本发明的特征在于，以速度扫描范围内的反演结果作为Hough变换圆检测的输入时，只有当参数值与真实值一致时，反演出的曲线形状才最接近于圆，而其他情况则类似于椭圆或双曲线，曲线越接近于圆时，参数域中的‘亮点’积累的值越大，通过峰值检测技术，可以得到管线目标位置及管径等参数的最佳估计值，检测结果误差小，可信度高。

为了达到上述目的，本发明采取的步骤如下：

1）利用对接收的雷达剖面进行预处理；

2）对预处理后的剖面同时进行边缘检测，找出类双曲线信号；

3）结合目标的特征和波动反射理论，根据已知的接收天线位置及信号到达时间，求取目标的位置、管径以及介质的速度等参数的初始值；

4）进行速度扫描，并针对各扫描值，反演出各空间散射点位置信息；

5）基于Hough变换，对步骤4）中的反演结果做圆检测，得到变换后参数域的值，并对其进行峰值检测，记录峰值对应的介质速度，并得到目标位置以及管径等参数的最佳估计值。

在上述技术方案中，步骤1）所述预处理是从中接收信号中滤除背景噪声及各类杂波信号。预处理技术包括：零偏矫正、数字滤波、小波软阈值去噪、测线相关去除直达波等。

步骤2）利用小波模极大值方法对剖面边缘进行检测，也可以借助其他图像边缘检测方法，包括方向算子、梯度算子、canny算子及Sobel算子等。

步骤3）所述对地下管径目标进行雷达探测时，在雷达剖面内，已知的是各个接收天线的坐标和它所接收到信号的时刻,不知道的是目标的位置和尺寸以及介质传播速度。如图2，设目标位置的空间坐标为，管径为R，天线第个接收点坐标为，检测到类双曲线信号的时刻为，介质传播的速度为，在目标散射点和第个接收天线之间的走时方程见公式（1）。

（1）

方程中有四个位置量，利用不同接收信号的已知量，做差分，我们可以计算出目标的位置、尺寸和速度的信息。计算方法见公式（2）。

（2）

在剖面中定位一个未知点的坐标位置和管径尺寸，需要在双曲线上获取四个以上的坐标，为了提高精度，可以滑动取点，再做平均。

步骤4）中所述的真实散射点位置的反演方法见公式（3），假设真实的散射点坐标为，其中介质传播速度为步骤3）中求得的初始值在一定范围内进行扫描的结果，地下管线目标的位置和管径可以依赖公式（2）得到不同速度扫描值下的对应值。

（3）

步骤5）利用Hough变换将目标雷达图像变换到了参数空间，在坐标平面上确定一个圆需要三个参数——圆的半径、圆心的坐标，因此圆的Hough变换是一个以圆的半径和圆心坐标为参数的三维空间。检测图像平面上的圆的问题就转换为检测参数空间上三维锥面的交点位置问题。图像平面的方程转化为参数平面上的示意图如图3。采用基于最小距离原则的递归方法来找到所有的局部极大值，这样就可以检测出所有可能的圆心。将圆心定位之后，采用阈值化的方法就能找到可能的半径。只有当参数值与真实值一致时，反演出的曲线形状才最接近于圆，而其他情况则类似于椭圆或双曲线，椭圆和双曲线在Hough变换后，不能形成有效的交点，那么圆心就会迅速扩散，圆心的累积值将会降低，记录下扫描结果对应的圆心累计值，采用峰值检测技术，就可以找到反演后最接近圆时对应的速度值，进而得到对应的管线目标位置及管径等参数的最佳估计值。

附图说明

图1为本发明方法的流程图。

图2为目标在雷达探测中的示意图。

图3为圆的参数空间表示。

图4为预处理前后的雷达剖面图。

图5为边缘检测后提取的目标双曲线。

具体实施方式

下面结合附图和具体实施方式对本发明作进一步详细描述：在一个长为2.5m，高为2m的空间内，填满水泥，中间放置一个金属管，中心位置为(1m, 0.8m), 天线频率为1GHz，相对介电常数为9。管径从0.35m~0.6m，每次变化量为0.05m。天线位置为离地表0.05m，初始横坐标为0.0875m，每次位移量为0.02m，共采200道。

在软件中将采集的信号，按照如下步骤进行处理：

1）利用FIR滤波器对信号进行带通滤波，并进行道内和道间零偏矫正，小波软阈值去噪和SVD变换，滤除随机噪声和直达波以及地表反射信号的干扰，图4为预处理前后的雷达剖面图；

2）对预处理后的剖面利用小波模极大值方法进行边缘检测，并将边缘检测结果进行相关，找出目标曲线对应在剖面上的位置，保留原始剖面中相应位置的边缘检测结果，其它滤除，图5为边缘检测后的目标双曲线提取结果；

3）对边缘上的点进行循环取点，结合公式（1）和（2），求取目标的位置、管径的尺寸以及介质的速度等参考参数；

4）一定的范围内扫描参数；

5）对每个扫描结果利用公式（3）进行目标真实的空间散射点位置信息反演；

6）基于Hough变换，对所有的反演结果做圆检测，记录变换后参数域的‘圆心’累积值，并对其进行峰值检测，得到对应的介质速度值，并计算出目标的位置以及管径等参数，作为最佳的估计值。

利用前面阐述的地下管线目标参数自识别方法，计算目标的位置和管径特征，如表1所示。其中坐标Y的误差来源：测线间隔为0.02m，起始坐标为0.0875m，离真实坐标最近的值为0.9875m。从管径估计结果来看，误差随着管径的增大有所减小，误差控制在5%以内，算法精度高。

表 1

Claims

1.本发明涉及一种基于探地雷达的地下管线参数自识别方法，其特征在于，可以根据探地雷达接收到管线的散射信号反演出管径目标的位置及管径，包括如下步骤：

步骤一：预处理，首先采用零偏矫正、数字滤波、小波软阈值去噪、从接收信号中滤除背景噪声及各类杂波信号，然后利用SVD变换的优势，分离背景和目标信号；

步骤二：双曲线边缘提取，利用小波模极大值方法对预处理后的剖面进行边缘检测，找出目标曲线位置并最大限度的保留类双曲线信号在原始剖面中的形状特征；

步骤三：目标初始特征参数求取，结合目标的特征和波动反射理论，求取介质的速度、目标的位置以及管径等参数，由于读取误差的存在，本步骤提取参数与实际值存在一定误差，只能作为反演计算的初始参数；

步骤四：以3中计算得到的参数为初始参数，在一定范围内对速度参数值进行扫描，针对每次的扫描值，对地下管线目标的真实散射点位置进行反演，并计算出相应的管径参数，反演的具体数学表达式为：

Y_{i}^{'} = Y_{i} + \sqrt{\frac{{({Vt}_{i} / 2)}^{2}}{1 + {(1 / (Y_{0} - Y_{i}))}^{2} Z_{0}^{2}}}

公式中(Y′_i,Z′_i)真实的散射点坐标，V为以3中求得的参数为初始参数并在一定范围内进行扫描的速度值、(Y₀,Z₀)为3中求取的初始目标位置、(Y_i,0)是天线第i个接收点坐标，t_i是检测到类双曲线信号的时刻；

步骤五：基于Hough变换，对反演结果做圆检测，记录每个反演信息在参数域中对应的峰值，进行峰值检测，作为识别管径参数的依据，得到介质速度、目标位置以及管径等参数的最佳估计值，从而达到管径检测的目的。

2.如权利要求1所述的方法，其特征在于，所述步骤三中可以在未知介质速度的情况下，对目标位置、尺寸和速度信息进行计算，为步骤四的反演提供初始参数，计算的具体数学表达式为：

V = 2 \sqrt{\frac{t_{i} (Y_{k} - Y_{j}) (Y_{l} - Y_{k}) (Y_{l} - Y_{j}) + t_{j} (Y_{l} - Y_{k}) (Y_{i} - Y_{l}) (Y_{k} - Y_{i}) - t_{k} (Y_{j} - Y_{i}) (Y_{i} - Y_{l}) - t_{l} (Y_{j} - Y_{i}) (Y_{k} - Y_{j}) (Y_{k} - Y_{i})}{t_{i} t_{j} (Y_{l} - Y_{k}) (t_{j} - t_{i}) - t_{k} t_{j} (t_{k} - t_{j}) (Y_{i} - Y_{l}) - t_{l} t_{k} (t_{l} - t_{k}) (Y_{j} - Y_{i}) - t_{i} t_{l} (t_{i} - t_{l}) (Y_{k} - Y_{j}) - t_{i} t_{k} (t_{k} - t_{i}) (Y_{i} - Y_{j}) - t_{j} t_{l} (t_{l} - t_{j}) (Y_{k} - Y_{i})}}

R = - \frac{(Y_{j} - Y_{i}) (Y_{k} - Y_{j}) (Y_{i} - Y_{k}) + (t_{i}^{2} (Y_{k} - Y_{j}) + t_{j}^{2} (Y_{i} - Y_{k}) + t_{k}^{2} (Y_{i} - Y_{k}) + t_{k}^{2} (Y_{j} - Y_{i})) {(V / 2)}^{2}}{V (t_{i} (Y_{k} - Y_{j}) + t_{j} (Y_{i} - Y_{k}) + t_{k} (Y_{j} - Y_{i}))}

Y = - \frac{1}{2} \frac{RV (t_{j} - t_{i}) + (Y_{k} - Y_{j}) (Y_{i} - Y_{k}) + {(V / 2)}^{2} (t_{j}^{2} - t_{i}^{2}) - Y_{j}^{2} + Y_{i}^{2}}{(Y_{j} - Y_{i})}

Z = \sqrt{(R + {Vt}_{j} / 2) - {(Y_{j} - Y)}^{2}}

其中，V为速度值，R为管径，(Y_i,0)、(Y_j,0)、(Y_k,0)分别是天线第i、j、k个接收点坐标，t_i、t_j、t_k是分别是天线第i、j、k个接收点接收信号中检测到类双曲线信号的时刻，(Y,Z)为管线目标圆心位置对应的坐标。

3.如权利要求1所述的方法，其特征在于，所述步骤五中，采用Hough变换对反演结果进行圆检测，采用基于最小距离原则的递归方法来找到所有参数域内的局部极大值，对圆心进行定位，当参数值与真实值一致时，反演出的曲线形状最接近于圆，其他情况类似于椭圆或双曲线，椭圆和双曲线在Hough变换后，不能形成有效的交点，圆心的累积值会降低，记录下扫描结果对应的圆心累计值，采用峰值检测技术，找到反演后最接近圆时对应的速度值，进而得到对应的管线目标位置及管径等参数的最佳估计值。