一种探地雷达记录剖面中目标顶点的检测方法与流程

本发明属于探地雷达探测与应用技术领域，具体涉及运用探地雷达图像检测目标顶点的方法。

背景技术：

探地雷达(groundpenetratingradar，gpr)是一种有效的浅层隐藏目标探测技术，利用电磁波在媒质电磁特性不连续处产生的反射和散射实现非金属覆盖区域中目标的成像探测。gpr回波的幅度和时延包含目标位置及电磁散射特性等信息。gpr沿一维测线进行空间扫描，在测线的每个孔径点处向地下区域发射电磁波并接收散射回波。每个孔径点处接收一道回波数据，多个孔径点接收的回波数据按列排列，就形成了gpr记录剖面。通过对gpr记录剖面进行处理以实现对地下未知区域进行探测和参数反演是一种无损检测方式。gpr可以在岩石、土壤、冰、淡水、公共道路以及各种结构物等介质使用，还可以探测到地表下的物质、材质变化、空隙和裂隙等，广泛应用于道路工程、建筑工程、考古等无损探测领域。

gpr在扫描的过程中，当电磁波遇到与背景媒质电磁特性不同的目标时会产生电磁散射，部分散射波被gpr接收天线接收。随着天线的移动，地下目标的散射回波在gpr记录剖面中就形成了一条时延曲线。时延曲线的顶点与目标的真实空间位置有一一对应关系。从图形图像学角度分析，只需确定gpr记录剖面中时延曲线的顶点位置，便可以获得地下目标的空间分布。因此，gpr记录剖面中顶点位置的检测算法就成了gpr对地探测中目标定位的关键问题。国外对gpr记录剖面中的顶点检测问题进行了大量的研究。文献1“liu,m.wang,andq.cai,“thetargetdetectionforgprimagesbasedoncurvefitting,”inproc.3rdint.congr.imagesignalprocess.(cisp),vol.6.yantai,china,oct.2010,pp.2876–2879”提出了通过对gpr图像中的曲线进行曲线拟合和边缘检测，从而对图像中的目标进行提取的方法。文献2“z.w.wang,m.zhou,g.g.slabaugh,j.zhai,andt.fang,“automaticdetectionofbridgedeckconditionfromgroundpenetratingradarimages,”ieeetrans.autom.sci.eng.,vol.8,no.3,pp.633–640,jul.2011”提出了通过对gpr图像中的曲线使用偏微分方程式的数学方法和模板匹配的方式，确定图像中的双曲线顶点位置，以此来对钢筋混凝土桥面腐蚀情况进行检测。文献3“parneetkaur；kristinj.dana；franciscoa.romero；nenadgucunski,“automatedgprrebaranalysisforroboticbridgedeckevaluation,”ieeetransactionsoncybernetics.vol.46,2016,pp.2265–2276”提出了基于图像的梯度特性和曲线拟合方式的双曲签名，通过集成机器学习分类的方法，来检测桥面钢筋腐蚀情况。文献4“laurencemertens；raffaelepersico；loredanamatera；sébastienlambot,“automateddetectionofreflectionhyperbolasincomplexgprimageswithnoaprioriknowledgeonthemedium,”ieeetransactionsongeoscienceandremotesensing.vol.54,2016,pp.580-596”提出了基于gpr图像进行边缘检测结果中的边缘点生成双曲线拟合函数的目标检测方法。上述方法，在信噪比较高的情况下，有较高的检测精度，但在信噪比较低的情况下，gpr剖面记录图像中目标双曲线边缘轮廓发生变化，其平滑性、对称性等曲线特征会遭到破坏，因此会造成目标点的漏检或误检的情况，直接影响目标的检测精度。而且，边缘检测方法需要手动调整阈值，实时性较差。因此，有必要针对现有技术中存在的问题设计一种检测精度高，且能实现自动化的顶点检测方法。

技术实现要素：

本发明所要解决的技术问题是，针对现有技术的不足，提供一种基于自动阈值边缘检测和模板匹配的gpr记录剖面中目标顶点的检测方法，可以在有噪声的情况下，有更好的检测性能，且能实现自动化的顶点检测。

本发明的技术方案如下：

一种探地雷达记录剖面中目标顶点的检测方法，包括以下步骤：

步骤1：对原始探地雷达记录剖面进行预处理，并计算预处理后的探地雷达记录剖面沿测线方向的一维能量曲线，提取该能量曲线的局部极大值个数，据此设定该探地雷达记录剖面的目标顶点个数的预设区间(取值范围)；

步骤2：对预处理后的探地雷达记录剖面进行自动阈值匹配的边缘检测，得到一幅二值图像，根据二值图像估计可能存在的目标顶点；

步骤3：根据估计出的可能存在的目标顶点坐标建立匹配模板，使用匹配模板与边缘检测得到的二值图像进行模板匹配，舍弃匹配相似度低的顶点，得出匹配后的目标顶点；

步骤4：对匹配后的目标顶点进行聚类分析，去除虚假顶点，得出检测结果。

进一步地，所述步骤1中包括以下步骤：

步骤1.1、记原始的探地雷达记录剖面为矩阵in×m，其第x行第t列元素记为in×m[x,t]，x表示采集数据的方向维，其取值范围为[1,n]，n表示沿测线方向进行采样的孔径点数；t表示时间维，其取值范围为[1,m]，m表示某孔径点处采集到的一维散射回波经数字采样后的采样数点；总采样时长为tns，时间采样间隔为in×m[x,t]表示原始的探地雷达记录剖面上第x行第t列像元的像素值对in×m进行预处理，即去直达波和幅度归一化处理，结果记为矩阵i″n×m；

步骤1.2、根据公式计算每道数据的能量，由此得出探地雷达记录剖面的一维能量曲线p(x)；其中i″n×m[x,t]表示矩阵i″n×m第x行第t列元素；

步骤1.3、寻找一维能量曲线p(x)上的局部极大值点，记为(x,p(x))；

步骤1.4、统计所有局部极大值点的个数，记为s；据此设定该探地雷达记录剖面的目标顶点个数的预设区间为[smin,smax]，其中smin为round(u·s)，smax为round(v·s)，u和v为经验参数，round(·)表示四舍五入取整。

进一步地，u＝0.8，v＝1.2。

进一步地，所述步骤1.1中对in×m进行去直达波和幅度归一化处理的步骤如下：

定义一个元素全为1的矩阵one1×m，从矩阵in×m中取出第一列元素，记为in×1，根据公式(1)对矩阵in×m进行去直达波处理，得到矩阵i′n×m；再根据公式(2)对矩阵i′n×m进行幅值归一化处理，得到i″n×m；

i′n×m＝in×m-in×1·one1×m(1)

i″n×m＝i′n×m/max{i′n×m[x,t]}(2)。

进一步地，所述步骤2包括以下步骤：

步骤2.1：设定初始阈值α＝α0；

步骤2.2：采用canny边缘检测算法，根据设定的阈值对预处理之后的探地雷达记录剖面进行边缘检测，得出对应阈值下的二值图像bα(x,t)；

步骤2.3：根据二值图像bα(x,t)估计出可能存在的目标顶点：

首先，对该二值图像bα(x,t)先进行行遍历，再进行列遍历，找到二值图像bα(x,t)中所有像素值为1的点；定义二值图像bα(x,t)中所有像素值为1的点为二值图像的边缘点，记为e(xm,tm),m＝1,...,m′，其中m为边缘点的标号，m′为边缘点的总个数；将连续的边缘点连接成边缘曲线；

然后，对所有边缘点逐点进行处理，估计可能存在的顶点：

先针对任一个边缘点e(xm,tm)，计算出该点右边连续且具有相同梯度的边缘点数量，记为n′(即边缘点右边连续的n′个像素点的像素值都为1，即它们都为边缘点，且具有相同的梯度；计算边缘点的梯度即求图像梯度；求图像梯度可以把图像看成二维离散函数，图像梯度就是这个二维离散函数求导；图像梯度：g(x,t)＝dx(i,j)+dt(i,j)；dx(x,t)＝iα(i+1,j)-iα(i,j)；dy(x,t)＝iα(i,j+1)-iα(i,j)，其中，iα(i,j)是二值图像bα(x,t)上位置为(i,j)处的像素点的像素值)；

再基于边缘点e(xm,tm)设定子图像[xm:(xm+n′),(tm-n′):(tm+n′+n′)]，即取二值图像bα(x,t)中第xm行到xm+n′行，第tm-n′列到第tm+n′+n′列之间的区域作为子图像，其中n′用来设定子图像大小，在本方法中取值为5；

若子图像的左边[xm:(xm+n′),(tm-n′):tm]边缘点的梯度为上升趋势，且子图像的右边[xm:(xm+n′),(tm+n′):(tm+n′+n)]边缘点的梯度为下降趋势，则说明在这个子图像中可能存在目标顶点，计算该子图像中可能存在的目标顶点坐标为(xm,tm+t′)，其中t′为并记录该目标顶点坐标，否则说明这个子图像中不存在目标顶点；

针对每一个边缘点进行上述处理后，得到二值图像bα(x,t)中所有可能存在的目标顶点，记为eα(xk,tk),k＝1,2,…,k，(xk,tk)为第k个可能存在的目标顶点的坐标，k为可能存在的目标顶点的个数；

步骤2.4：判断k是否在目标顶点个数的预设区间[smin,smax]内，若是，则结束步骤2，取当前估计得到的eα(xk,tk),k＝1,2,…,k进行后续步骤3的处理；否则递增阈值，设定阈值α＝α+δα，并返回步骤2.2重复上述处理过程，直到α＝α1为止。若α＝α1，而k还是不在预设区间[smin,smax]内，说明探地雷达记录剖面图像的信噪比特别小，图像中的虚假顶点较多(噪声干扰太严重)，一般不会出现这种情况；如果出现了这种极端的情况，则无法检测出图像中的目标顶点。

进一步地，所述α0＝0.01，δα＝0.02，α1＝0.9。

进一步地，所述步骤3包括以下步骤：

步骤3.1：设置标准双曲线模板方程为其中，a＝t0，其中t0表示探地雷达到地下目标的最小时延，x0表示最小时延对应的横向探测距离，εr表示目标的相对介电常数；

步骤3.2：初始化k＝1；

步骤3.3：基于目标顶点eα(xk,tk)建立子二值图像；

基于步骤2中找出的可能存在的目标顶点eα(xk,tk)，建立子二值图像b′α(x,t),x∈[xkmin,xkmax],t∈[tkmin,tkmax]，其中xkmin、xkmax、tkmin和tkmax分别为xk、xk+x′、tk-b和tk+b，x′＝0.2·t0，t0是天线的周期，f为探地雷达天线的频率；

步骤3.4：基于目标顶点eα(xk,tk)建立匹配模板；

在有噪声的情况下，实际数据图像的边缘曲线与标准双曲线模板对比会有一定的偏差，因此，定义一个参数δb表示偏差的大小，来对标准双曲线模板进行修正。随着图像噪声的增大，对应的边缘检测的阈值也随之增大，因此，边缘检测之后的边缘曲线与标准双曲线模板对比的偏差也会随之增大，检测阈值与参数δb具有线性关系；设定δb＝10×α，其中α表示边缘检测中的阈值。修正后的双曲线为将[tkmin,tkmax]内的整数t依次代入修正后的双曲线得到相应的x，满足x∈[xkmin,xkmax]的所有解为修正后的双曲线对应的边缘点；将该修正后的双曲线对应的边缘点绘制在与上述子二值图像b′α(x,t),x∈[xkmin,xkmax],t∈[tkmin,tkmax]大小相同的二值图像中，得到匹配模板，记为mk(x,t),x∈[xkmin,xkmax],t∈[tkmin,tkmax]，匹配模板中的两条双曲线即分别由修正后的双曲线和对应的边缘点构成的双曲线，称为理论双曲线【由于修正后的双曲线对中存在±δb，因此会得到两条理论双曲线，它们具有相同的顶点，由附图6可以看出，其方框中，中间的较长的曲线为需要检测的边缘曲线，而在其上下的两条曲线则为对应的±δb的标准双曲线】；

步骤3.5：使用匹配模板mk(x,t),x∈[xkmin,xkmax],t∈[tkmin,tkmax]与子二值图像b′α(x,t),x∈[xkmin,xkmax],t∈[tkmin,tkmax]进行模板匹配，计算匹配相似度ρk，并将ρk与预设的匹配相似度阈值ρ0相比；如果ρk>ρ0，则认为e(xk,tk)为匹配的顶点，否则认为e(xk,tk)不是匹配的顶点；

步骤3.6：设置k＝k+1，若k≤k，则重复步骤3.3～步骤3.5对下一个顶点建立子二值图像、匹配模板、进行匹配相似度计算及与匹配相似度阈值的对比判断；若k>k则结束，将得到的所有匹配的顶点记为v(xβ,tβ),β＝1,...,b，(xβ,tβ)为第β个匹配的目标顶点的坐标，b为匹配的目标顶点的总个数。

进一步地，所述步骤3.5中，通过下述公式计算匹配相似度ρk：

ρk＝((ρ1+ρ3+ρ4+ρ5)·2+ρ2+ρ6)·0.1(3)

其中ρ1、ρ2、ρ3、ρ4、ρ5和ρ6分别表示模板匹配约束参数dreli、sym、drti、drxi、disi和sur的匹配相似度，计算公式如下：

各个模板匹配参数的计算方法如下：

将缓冲区，即子二值图像b′α(x,t),x∈[xkmin,xkmax],t∈[tkmin,tkmax]以直线x＝x0为界分成b′α(x,t),x∈(xkmin,x0),t∈[tkmin,tkmax]和b′α(x,t),x∈(x0,xkmax),t∈[tkmin,tkmax]左右两个部分，即左右缓冲区；

1)根据以下公式计算左右缓冲区的相对边缘点密度：

其中，drel1和drel2分别表示左缓冲区和右缓冲区的相对边缘点密度；ndata1和ndata2分别表示左缓冲区和右缓冲区中的边缘点【即缓冲区中像素值为1的点】的数量；ntheo是匹配模板mk(x,t),x∈[xkmin,xkmax],t∈[tkmin,tkmax]上的边缘点【即修正后的双曲线对应的边缘点】数量；在理想情况下：dreli应该是趋近100％，但其可接受范围为50％到270％之间。

2)根据左缓冲区和右缓冲区的相对边缘点密度计算边缘曲线的对称性规范化系数：

在理想情况下，sym趋近于1，其可接受的范围为sym≥0.6，这时可以从对称性的角度说，可以接受其为双曲线。

3)计算左缓冲区和右缓冲区相对理论双曲线在t方向上与x方向上的数据范围；

左缓冲区和右缓冲区相对理论双曲线在t方向上的数据范围计算方法如下：

其中，drt1和drt2分别为左、右缓冲区相对理论双曲线在t方向上的数据范围；trt1和trt2分别是指将匹配模板匹配至缓冲区上时，匹配模板上的理论双曲线在左、右缓冲区中t方向上的数据范围(最小值和最大值的差值的绝对值)；adrt1和adrt2是分别是指左、右缓冲区在t方向上的数据范围。

同理，左缓冲区和右缓冲区相对理论双曲线在x方向上的数据范围计算方法如下：

其中，drx1和drx2分别为左缓冲区和右缓冲区相对理论双曲线在x方向上的数据范围；trx1和trx2分别是指将匹配模板匹配至缓冲区上时，匹配模板上的理论双曲线在左缓冲区和右缓冲区中x方向上的数据范围(最小值和最大值的差值的绝对值)；adrx1和adrx2是分别是指左缓冲区和右缓冲区在x方向上的数据范围。在理想条件下，drti和drxi的值趋近于1；实际情况下，drti和drxi的可接受范围为：drti≥0.7，drxi≥0.7。

4)计算左缓冲区和右缓冲区相对于理论双曲线在t方向上的数据分布：

其中，dis1和dis2分别表示左缓冲区和右缓冲区相对于理论双曲线在t方向上的数据分布，nmatch1和nmatch2分别表示左缓冲区和右缓冲区中的边缘点在落在匹配模板中两条理论双曲线之间的个数【记和分别表示模板中两条理论双曲线上的边缘点，且xm1≤xm2，对于左右缓冲区中的边缘点e(xm,tm)，如果满足tm1＝tm＝tm2且xm1≤xm≤xm2，则说明边缘点e(xm,tm)落在匹配模板中两条理论双曲线之间】；

如果只考虑相对数据范围和相对边缘分布密度还是不足以判定双曲线，如，在相互干扰的情况下，产生的图形，这时还得考虑数据分布情况。判定的条件：disi≤0.7，就可以认为存在双曲线。

5)缓冲区中边缘曲线外的边缘点的相对数量sur，计算方法如下。

即使数据分布情况满足disi≤0.7，但有可能在缓冲区内部的数量比较高，会影响对双曲线的判断。其中：nsur表示在缓冲区(包括左右两个缓冲区)中边缘曲线外的离散的边缘点的数量【对于左缓冲区内离散的边缘点e(xd,td)，若其相对于该缓冲区中边缘曲线上的某个边缘点e(xc,tc)，满足xd<xc且td<tc；对于右缓冲区内离散的边缘点e(xd,td)，若其相对于该缓冲区中边缘曲线上的某个边缘点e(xc,tc)，满足xd>xc且td<tc，则说明这些边缘点在边缘曲线外部】；nbf表示在缓冲区(包括左右两个缓冲区)中边缘曲线内的离散的边缘点的数量【对于左缓冲区内离散的边缘点e(xd,td)，若其相对于该缓冲区中边缘曲线上的某个边缘点e(xc,tc)，满足xd>xc且td>tc；对于右缓冲区内离散的边缘点e(xd,td)，若其相对于该缓冲区中边缘曲线上的某个边缘点e(xc,tc)，满足xd<xc且td>tc，则说明这些边缘点在边缘曲线内部】；在理想情况下，sur的值比较小，即nbf变小，则对应的sur变小。因此，在sur≤0.7的情况下，是可以接受的。

设定模板匹配约束参数ρ1～ρ6均等于1时的模板匹配相似度作为匹配相似度阈值，即ρ0＝1。

进一步地，所述步骤4包括以下步骤：

步骤4.1：对步骤3中得出的匹配的目标顶点v(xβ,tβ),β＝1,...,b进行聚类分析得到簇cn，n＝1,...,n1，其中n为簇序号，n1为簇数，各个簇包含的目标顶点记为个数qn，n＝1,...,n1；

步骤4.2：设定初始值n＝1，γ＝1；

步骤4.3：若qn＝1成立，则将簇cn中包含的目标顶点作为该簇中过滤后的目标顶点vγ；若qn＝1不成立，则比较簇cn中各个目标顶点的能量值大小，得到该簇中能量值最大的目标顶点，若该簇中存在多个能量值最大且相等的目标顶点，则取这些目标顶点中对应的时间点最小的目标顶点作为该簇中过滤后的目标顶点；若该簇中只存在一个能量值最大的目标顶点，则将该目标顶点作为该簇中过滤后的目标顶点vγ；

步骤4.4：判断n＝n1是否成立，若不成立，则设定返回步骤4.3对下一个簇中的顶点进行过滤；若n＝n1成立，则结束；遍历n1个簇，得到原始探地雷达记录剖面中的顶点估计序列vγ，γ＝1,...,n1，作为顶点检测结果。

有益效果：

本发明提出一种基于自动阈值边缘检测和模板匹配的gpr记录剖面中目标顶点的检测方法。对gpr记录剖面进行边缘检测，得出图像中的曲线边缘轮廓。利用边缘点的梯度变化关系，预估出顶点的坐标位置；再基于顶点的位置生成相应的子图像窗口，结合双曲线拟合生成标准的双曲匹配模板；针对边缘点设定上述六个匹配条件，得到模板匹配的相似度，去除相似度较低的顶点；最后，对所得顶点进行聚类分析，去除虚假顶点，进而得到要检测的目标顶点。本发明中自动匹配最佳阈值，不必手动调整，提高了算法的灵活性。在有噪声的情况下，实际数据图像的边缘点与标准双曲线模板对比会有一定的偏差，因此，定义一个参数δb表示偏差的大小；通过对δb参数的控制，使得在有噪声的情况，也有较好的检测性能，可应用于实际含噪gpr记录剖面中的目标检测。

附图说明

图1示出了本发明的处理流程图；

图2示出了gpr原始b-scan图像；

图3示出了原始数据经预处理之后的图像；

图4示出了预处理之后图像的一维能量曲线；

图5示出了边缘检测后预估的顶点；

图6示出了根据各顶点的坐标生成的匹配模板曲线；

图7示出了进行模板匹配舍弃相似度较低的顶点后的检测结果；

图8示出了运用聚类分析去除虚假顶点得到的顶点检测结果；

图9示出了运用传统方法对图2所示的gpr图像进行顶点检测的结果。

具体实施方式

以下将结合附图和具体实施例对本发明做进一步详细说明。

本发明主要自动检测gpr图像目标顶点，如图1所示，采用gprmax2.0仿真软件生成原始图像，进行加噪处理模拟产生含噪的实测gpr图像。运用本发明的算法对含噪gpr图像进行处理，首先对含噪的原始图像进行预处理，然后对预处理后的探地雷达记录剖面采用自动匹配阈值的方式进行边缘检测，设定初始阈值，进行边缘检测和顶点估计，判断顶点个数是否位于预设区间。如是，则取当前的边缘检测结果进行下一步。如否，则递增阈值，直到预估计的目标顶点个数在预置区间为止。接着提取边缘检测得到的二值图像中的目标顶点，基于这些目标顶点的坐标位置建立匹配模板，运用该模板与边缘检测得出的二值图像进行模板匹配并计算匹配相似度，设定匹配阈值，舍弃匹配相似度低于该阈值的顶点，得出匹配后的目标顶点。再基于聚类分析方法对目标顶点进行过滤，去除虚假顶点，最后得到探地雷达记录剖面中的目标顶点的检测结果。

实施例1：

本实例中，通过gprmax2.0仿真软件生成探地雷达原始图像，gpr扫描的横向区间长度为2.5m，总的横向维采样点数为920点，相邻两孔径点的采样间隔为0.25mm。纵向维为时间维，总的采样点数为1357点，总的采样时长为12ns，该数据的尺寸为1357×920的二维矩阵。仿真生成的gpr数据为不含噪声的理想数据，将其归一化。为模拟真实的gpr探测场景，加入噪声，使得信噪比snr等于10db，加噪后的gpr记录剖面如图2所示。

对gpr原始图像进行预处理，利用公式(1)、(2)将图像中的直达波去除并对幅值进行归一化，得出预处理的图像如图3所示；

绘制出原始图像的一维能量曲线，如图4所示；根据原始图像的一维能量曲线特征，在极值点处可能存在顶点，求出极值点的个数为23，设定可能存在顶点的数量取值范围[18,28]；

接下来对预处理的图像进行边缘检测找出最佳的边缘检测阈值，采用canny边缘检测函数，设定初始阈值α0为0.01，α以0.02的步长递增，到0.9结束。根据每个阈值获取边缘检测的二值图像bα(x,t)，再根据二值图像估计出可能存在的目标顶点。对该二值图像bα(x,t)先进行行遍历，再进行列遍历，找到二值图像bα(x,t)中所有像素值为1的点；定义二值图像bα(x,t)中所有像素值为1的点为二值图像的边缘点，记为e(xm,tm),m＝1,...,m′，其中m为边缘点的标号，m′为边缘点的总个数；

然后，对所有边缘点逐点进行处理，估计可能存在的顶点：针对任一个边缘点e(xm,tm)，计算出该点右边连续且具有相同梯度的边缘点数量记为n，再基于边缘点e(xm,tm)，设定子图像[xm:(xm+n′),(tm-n′):(tm+n′+n′)]，其中n′表示像素点的数目，用来设定子图像大小，在本方法中取值为n′＝5。若子图像的左边[xm:(xm+n′),(tm-n′):tm]的边缘点的梯度为上升趋势，且子图像的右边[xm:(xm+n′),(tm+n′):(tm+n′+n)]的边缘点的梯度为下降趋势；则可以说明在这个子图像中可能存在顶点，此时计算顶点坐标为(xm,tm+t′)，其中t′为并记录该顶点坐标；根据上述子图像中边缘点的梯度关系找出图像中所有的顶点估计结果记为eα(xk,tk),k＝1,2,…,k，k为得出顶点的数量。在边缘检测二值图像中标出所有的估计顶点如图5所示，根据上述处理过程，在本实例中找到23个估计顶点，即k＝23。

设置标准的双曲线模板方程为其中，a＝x0,根据上述找出的可能存在的顶点eα(xk,tk),k＝1,2,…,k，然后遍历每一个顶点，当k＝1时，再基于顶点(226,150)建立[xk:(xk+x^·),(tk-b):(tk+b)]大小的子图像记为bα(x,t),x∈[145,155],t∈[226,232]，在有噪声的情况下，定义δb＝10·α，来对标准双曲线模板进行修正，修正后的双曲线方程为将该修正后的双曲线对应的边缘点绘制在设定的匹配模板中，就建立好了匹配的模板mk(x,t),x∈[145,155],t∈[226,232]，如图6所示。

下面进行模板匹配的过程，上述的步骤7，根据设定了6个匹配约束参数，先判断第一个参数。根据公式(4)，分别求出drel1＝2.75和drel2＝5.37，根据公式(5)，drel1和drel2，得到对称性规范化系数sym＝0.51；根据公式(6)得出drt1＝0.71，同理得出drt2＝0.71、dry1＝0.92、dry2＝2.53；接下来计算disi：根据公式(7)得出dis1＝0.50、dis2＝0.37；根据公式(8)得出sur＝-4.54，综合上述结果，通过公式(3)得出模板匹配相似度ρk＝0.6，小于相似度阈值。因此，认为(226,150)不是目标顶点。重复上述的判断方式，对所有的顶点进行上述的判断过程得到所有匹配之后顶点记为v(xβ,tβ),β＝1,...,b，并在原始图像上标出这些顶点的位置，如图7所示。

对匹配之后的目标顶点进行过滤：

对上述步骤得出的目标顶点进行聚类分析得到簇cn，n＝1,...,3，及各簇包含的顶点个数分别为q1＝3,q2＝2,q3＝3，设定初始值n＝1，γ＝1；对于簇cn，如果qn＝1成立，则vγ＝vn，其中vγ为过滤后的目标顶点，vn为簇cn中包含的唯一一个顶点；如果qn＝1不成立，则针对该簇提取其中能量值最大的顶点，顶点序号记为[m1,...,mp]，此时若该簇中只存在一个能量值最大的顶点，即p＝1成立的话，则直接认为为该簇中过滤后的目标顶点vγ；否则，若p≠1，说明该簇中存在多个能量值最大且相等的顶点，则提取其中对应时间点最小的顶点作为该簇中过滤后的目标顶点vγ，即如果n＝n1不成立，则设定对下一个顶点进行上述处理过程，遍历n1个簇，得到gpr记录剖面中的顶点估计序列vγ，也就是顶点检测结果，如图8所示。

运用传统方法对图2所示的gpr记录剖面进行顶点提取，结果如图9所示，采用手动调节阈值方式，共检测出7个目标顶点。标记1处实际只有一个目标顶点，而传统方法检测出了2个目标点，包含了1个虚假顶点。在图中标记2，3，4处没有目标顶点，而这3处都检测出了目标顶点，出现了误检的情况。本发明可以自动匹配阈值，提高了算法的灵活性，并且在低信噪比情况下有更高的检测精度，没有出现误检、漏检的情况。