CN105044785A - 一种基于模糊聚类与Hough变换的探地雷达地下管线探测方法 - Google Patents

Abstract

本发明提供一种地下管线反演方法，属于计算机数字图像处理领域。其特点在于综合采用双曲线渐近线和模糊C均值聚类以及聚类Hough变换反演地下管线位置及半径等参数。本发明基于共偏置模式探地雷达，通过B扫描得到雷达数据，再通过目标追踪获取供反演的一系列坐标。利用这一系列坐标所蕴含的变量特征，提取出雷达波速。结合模糊C均值聚类方法，达到获取准确雷达波波速的目的。在获取雷达波波速以后，再通过三角形相似得到一系列管线边缘坐标值，利用聚类Hough变换圆检测，进一步准确地反演出管线的主要参数。本发明具有预测准确、鲁棒性强的特点。

Description

一种基于模糊聚类与Hough变换的探地雷达地下管线探测方法

技术领域

本发明属于计算机数字图像处理领域，针对探地雷达图像的特点，利用双曲线渐近线、模糊C均值聚类以及聚类Hough变换反演地下管线参数。

背景技术

探地雷达发射机向地下介质发射超高频宽带短脉冲电磁波，遇到地表和不同介质的地下目标和界面时，部分脉冲波将反射回地面，被接收天线接收，然后进行数据处理，当目标物在天线信号范围之内并且信噪比适当时，隐蔽物可以被探地雷达探测出。探地雷达的测深和分辨率与以下几个因素有关：天线频率、发射功率、传播介质的电磁特性、目标物的形状和大小。本专利采用共偏置模式探地雷达，共偏置模式是最常见的探测方式，在进行数据采集时，雷达收发天线以固定间距沿测线移动探测，当收发天线间隔为零时称为无偏置模式。当收发天线同时移动一次获得一个雷达A扫描记录，收发天线沿测线同步移动得到一个B扫描雷达数据记录，横坐标记录接收天线的水平位移，纵坐标记录雷达波双程时间差，根据双程时差与地下介质中的波速即可估算出地下目标掩埋深度，再结合横坐标信息就能确定目标位置。

1957年，HugoSteinhaus首次提出K均值聚类的思想，直到1967年，JamesMacQueen才实现了该算法。K均值聚类，对于数据的划分是硬性的，每个样本数据都被严格地划分为属于某一类。但在实际问题中事物往往不满足“非此即彼”的条件，那就要考虑有模糊存在的问题，即有些事物或特征不是仅仅属于某一特定的类，而是“亦此亦彼”，只是属于不同类的程度有所不同。因此将模糊数学理论引入聚类分析中，采用模糊聚类分析可以取得更好的效果。模糊C均值聚类算法，是一种基于模糊划分的聚类算法，它的思想是使得被划分到同一簇的对象之间相似度最大，而不同簇之间的相似度最小，用隶属度确定每个数据点属于某个聚类的程度。1973年，J.C.Dunn和J.C.Bezdek提出了该算法，作为早期K均值聚类方法的一种改进。

发明内容

本发明的目的是通过探地雷达B扫描得到的双曲线渐近线和模糊C均值聚类算法准确反演地下管线参数。本方法具有良好的适应性、准确性和稳定性。

本方法按照以下步骤进行雷达波速预测：

步骤(A1)：雷达数据对地下管线检测成像具有双曲线特征，该双曲线仅有半支且开口向下，双曲线上共有n个元素，其横坐标为水平位置x_i，i＝1,2,……,n，纵坐标为信号往返时间t_i，i＝1,2,……,n。为减少误差数据影响，去除该双曲线顶点(x₀,t₀)附近的数据，余下m个数据；

步骤(A2)：求这m个数据构成的斜率，利用双曲线上元素连线的斜率逐渐接近渐近线的特点，将该双曲线上隔c个点连线的斜率作为渐近线的估计斜率k_j，j＝1,2,……,m-c；

步骤(A3)：对于一系列渐近线的估计斜率k_j，j＝1,2,……,m-c，我们可以借助公式v＝2λ/k得到一系列波速的估计值，其中为矫正系数；

步骤(A4)：利用模糊C均值聚类算法对m-c个波速估计值进行聚类，通过合理选取初始聚类数量s和校正系数λ，将元素最多的类别中心作为实际波速的结果

反演管线参数具体步骤如下：

步骤(B1)：首先需要粗略计算管线中轴坐标(xm,ym)，其中，xm为使得t_i取得最小值时x_i的值，即x₀，或者几个使得t_i取得最小值时x_i的均值，其中i＝1,2,……,n，(x_f,t_f)为曲线上距离顶点(x₀,t₀)较远的一点；

步骤(B2)：依据三角形相似则对应边成比例，我们可以得到管线边缘点坐标(x_ei,y_ei)，为了克服部分数据提取时产生的误差，这里要先对双曲线坐标(x_i,y_i)进行最小二乘拟合，得出拟合结果，代入方程，其中 $r=\sqrt{{(x_i-\mathrm{xm})}^2+{\left(\mathrm{ym}\right)}^2}-\frac{{\hat{v}t}_i}{2},x_{\mathrm{ei}}=x_i+(\mathrm{xm}-x_i)\frac{\hat{v}{t}_i}{r+{\hat{v}t}_i},$ $y_{\mathrm{ei}}=\mathrm{ym}\×\frac{\hat{v}{t}_i}{r+\hat{v}{t}_i},i=\mathrm{1,2},......,n;$

步骤(B3)：管线边缘点坐标(x_ei,y_ei)可能在管线边缘上，也可能不在管线边缘上，因此我们需要通过聚类Hough变换排除误差数据的影响，进而获取管线半径的估计值

步骤(B4)：管线半径的估计值作为已知，我们再结合之前计算得到的ym，便可以得到管线的中轴新坐标(xm,ym)，其中xm＝xm，综上，我们便得到了管线探测的参数 $\{\mathrm{xm},\mathrm{ym},\hat{r},\hat{v}\}.$

本发明中采用的双曲线最小二乘拟合步骤如下：

步骤(C1)：通过双曲线坐标(x_i,t_i)得到线性方程坐标(X_i,T_i)，其中

步骤(C2)：依据最小二乘原理计算线性方程的参数(A,B)，其中 $B=\frac{1}{m}({\mathrm{\Σ}}_1^n{Y}_i-A{\mathrm{\Σ}}_1^n{X}_i);$

步骤(C3)：计算拟合双曲线的估计参数其中

本发明中采用的聚类Hough变换步骤如下：

步骤(D1)：管线边缘点坐标(x_ei,y_ei)，i＝1,2,……,n，作为Hough变换圆检测的输入参数，同时还要定义半径的变化范围(r_min,r_max)和每次改变的步长Δr，另外还要确定模糊C均值聚类的初始聚类的数目s；

步骤(D2)：Hough变换圆检测将管线边缘点坐标(x_ei,y_ei)，i＝1,2,……,n作为圆心，任意两个边缘点(x_e1,y_e1)，(x_e2,y_e2)构造的圆的交点(x_c1,y_c1)和(x_c2,y_c2)可以通过以下公式计算， $x_{c1}=\frac{-b+\sqrt{b^2-4\mathrm{ac}}}{2a},x_{c1}=\frac{-b-\sqrt{b^2-4\mathrm{ac}}}{2a},$ 同时， $y_{c1}=\frac{y_{e1}+y_{e2}}{2}-\left(\frac{x_{e1}-x_{e2}}{y_{e1}-y_{e2}}\right)(x_{c1}-\frac{x_{e1}+x_{e2}}{2}),$ $y_{c2}=\frac{y_{e1}+y_{e2}}{2}-\left(\frac{x_{e1}-x_{e2}}{y_{e1}-y_{e2}}\right)(x_{c2}-\frac{x_{e1}+x_{e2}}{2}),$ 其中 $c=\frac{{(x_{e1}^2+x_{e2}^2+{(y_{e1}-y_{e2})}^2)}^2-4x_{e1}^2{x}_{e2}^2}{4{(y_{e1}-y_{e2})}^2}-{(r_{\min }+\mathrm{\Δr})}^2,$ $b=(x_{e1}+x_{e2})[1+{\left(\frac{x_{e1}-x_{e2}}{y_{e1}-y_{e2}}\right)}^2],a=1+{\left(\frac{x_{e1}-x_{e2}}{y_{e1}-y_{e2}}\right)}^2,$ 圆的交点可能有1个、2个或者0个；

步骤(D3)：利用模糊C均值聚类算法对圆的交点的集合(x_ck,y_ck)，k＝1,2,……,n(n-1)/2进行聚类，由于通常x_ck误差较小，因此主要对y_ck进行聚类，我们将元素最多的类别中心作为管线中轴y坐标的估计值再利用得到管径的估计值

本发明就有以下优点：

1、雷达波速预测准确，能够满足管线参数计算的需要。

2、有较强的鲁棒性，适合多种探地雷达数据的处理。

3、内存需求较低，避免了内存不足的情况。

附图说明

图1本发明预测雷达波速流程图

图2本发明反演管线参数流程图

图3本发明聚类Hough变换流程图

图4本发明简明流程图

具体实施方式

本发明采用共偏置模式探地雷达，利用B扫描得到的雷达数据，结合双曲线渐近线和模糊C均值聚类预测雷达波速。利用波速初步确定管线轴心坐标，通过三角形相似对应边成比例得到管线边缘坐标，再通过聚类Hough变换圆检测进一步确定管线轴心，从而反演出地下管线参数。

预测雷达波速流程如下：

(1)如图1所示，通过提取B扫描得到的雷达数据，得到一系列坐标(x_i,t_i)。这一系列坐标具有双曲线特征，在坐标中寻找顶点坐标(x₀,t₀)，去除顶点坐标附近数据。

(2)利用余下数据(x_i,t_i)估计斜率k_j，并通过计算进一步得到一系列雷达波速v_j。采用模糊C均值聚类计算雷达波速v_j的估计值

反演管线参数流程如下：

(1)如图2所示，首先利用一系列坐标(x_i,t_i)以及雷达波速的估计值粗略计算管线中轴坐标(xm,ym)，然后再进一步提取管线边缘点坐标(x_ei,y_ei)。

(2)对管线边缘点坐标(x_ei,y_ei)做聚类Hough变换，得到较为准确的管线中轴坐标(xm,ym)。从而进一步确定管线半径的估计值这样就得到了管线反演的参数

聚类Hough变换流程如下：

(1)将管线边缘点坐标(x_ei,y_ei)作为圆心，计算其交点集合(x_ck,y_ck)。

(2)利用模糊C均值聚类算法对交点y_ck聚类，从而得到管线中轴y坐标的估计值再利用进一步求出

Claims (4)

1.一种基于模糊聚类与Hough变换的探地雷达地下管线探测方法，该方法建立在雷达数据对地下管线检测成像的理论基础上，其特征在于，与理论推导的双曲线图形相对照，采用渐近线斜率作为雷达波速的计算依据，利用矫正参数和模糊C均值聚类算法获取准确的雷达波速，在计算雷达波速过程中，依次含有以下步骤：

步骤(A1)：雷达数据对地下管线检测成像具有双曲线特征，该双曲线仅有半支且开口向下，双曲线上共有n个元素，其横坐标为水平位置x_i，i＝1,2,……,n，纵坐标为信号往返时间t_i，i＝1,2,……,n。为减少误差数据影响，去除该双曲线顶点(x₀,t₀)附近的数据，余下m个数据；

步骤(A2)：求这m个数据构成的斜率，利用双曲线上元素连线的斜率逐渐接近渐近线的特点，将该双曲线上隔c个点连线的斜率作为渐近线的估计斜率k_j，j＝1,2,……,m-c；

步骤(A3)：对于一系列渐近线的估计斜率k_j，j＝1,2,……,m-c，我们可以借助公式v＝2λ/k得到一系列波速的估计值，其中λ为矫正系数；

步骤(A4)：利用模糊C均值聚类算法对m-c个波速估计值进行聚类，通过合理选取初始聚类数量s和校正系数λ，将元素最多的类别中心作为实际波速的结果

2.如权利要求1所述方法，其特征在于，依据三角形相似得到管线边缘坐标，再利用聚类Hough变换圆检测进一步确定管线的中轴坐标，从而获取管线半径在内的一系列参数，步骤如下：

步骤(B1)：首先需要粗略计算管线中轴坐标(xm,ym)，其中，xm为使得t_i取得最小值时x_i的值，即x₀，或者几个使得t_i取得最小值时x_i的均值，其中i＝1,2,……,n，(x_f,t_f)为曲线上距离顶点(x₀,t₀)较远的一点；

步骤(B2)：依据三角形相似则对应边成比例，我们可以得到管线边缘点坐标(x_ei,y_ei)，为了克服部分数据提取时产生的误差，这里要先对双曲线坐标(x_i,y_i)进行最小二乘拟合，得出拟合结果，代入方程，其中

步骤(B3)：管线边缘点坐标(x_ei,y_ei)可能在管线边缘上，也可能不在管线边缘上，因此我们需要通过聚类Hough变换排除误差数据的影响，进而获取管线半径的估计值

步骤(B4)：管线半径的估计值作为已知，我们再结合之前计算得到的ym，便可以得到管线的中轴新坐标其中综上，我们便得到了管线探测的参数

3.对于权利要求2所述方法中提到的双曲线计算最小二乘拟合，我们要通过降维的方法，对线性函数实现最小二乘拟合，然后得出双曲线的估计参数具体步骤如下：

步骤(C1)：通过双曲线坐标(x_i,t_i)得到线性方程坐标(X_i,T_i)，其中

步骤(C2)：依据最小二乘原理计算线性方程的参数(A,B)，其中

步骤(C3)：计算拟合双曲线的估计参数其中

4.对于权利要求2所述方法中提到的聚类Hough变换，其特征在于，利用模糊C均值聚类对Hough变换的交点进行聚类，其中元素最多的类别看做空间交点最集中的区域，该类别中心作为实际交点的估计位置，这就使得参数空间的统计问题转换为空间交点的聚类问题，以时间换取空间，大幅度减少了计算的内存消耗，具体步骤如下：

步骤(D1)：管线边缘点坐标(x_ei,y_ei)，i＝1,2,……,n，作为Hough变换圆检测的输入参数，同时还要定义半径的变化范围(r_min,r_max)和每次改变的步长Δr，另外还要确定模糊C均值聚类的初始聚类的数目s；

步骤(D2)：Hough变换圆检测将管线边缘点坐标(x_ei,y_ei)，i＝1,2,……,n作为圆心，任意两个边缘点(x_e1,y_e1)，(x_e2,y_e2)构造的圆的交点(x_c1,y_c1)和(x_c2,y_c2)可以通过以下公式计算，同时， 其中 圆的交点可能有1个、2个或者0个；

步骤(D3)：利用模糊C均值聚类算法对圆的交点的集合(x_ck,y_ck)，k＝1,2,……,n(n-1)/2进行聚类，由于通常x_ck误差较小，因此主要对y_ck进行聚类，我们将元素最多的类别中心作为管线中轴y坐标的估计值再利用得到管径的估计值