一种基于下开口扫描的双曲线线性化识别方法

1.本发明涉及雷达图像处理技术领域，具体涉及一种基于下开口扫描的双曲线线性化识别方法。


背景技术：

2.探地雷达(ground penetrating radar)技术是一种用于探测地下目标的地球物理探测方法。传统的裂缝监测方式有目视检查、锤击和取芯等，这些方法要么具有一定的破坏性，要么效率低且成本高。相比于传统检测，使用无损方式进行检测可以降低维护成本和风险。无损检测除了探地雷达方法，还包括电阻率法、低频电磁感应法等。与其他无损检测方法相比，探地雷达检测法探测速度快、探测过程连续、分辨率高、操作方便灵活、探测费用低等优点，被认为是在所有无损检测中连续性、安全性和抗干扰性方面最有效和可靠的。
3.城市道路的管道具有与周围介质不同的介电常数，由于在埋设管道时，管道往往水平埋设在所使用的探地雷达天线的有效深度范围内，而探地雷达的探测方向应垂直于管道方向(在实际应用中，埋地管道的方向可以通过现有的管道图初步确定，然后通过不同方向的多次检测来确定管道的具体方向)。因此，在探测管道时，gpr b扫描图像上会有一条双曲线，并且可以通过拟合生成的双曲线并检查其参数来获得介质中相关管道的位置和尺寸。
4.在实际应用中，由于系统噪声的影响，以及介质的异质性和相互波的相互作用，探地雷达图像可能是有噪声的。因此，使用计算机自动处理探地雷达图像并从中提取双曲线便成为了一项课题。在这一领域出现了许多策略来解决这一问题，包括但不限于基于霍夫变换的方法、最小二乘法、基于机器学习的方法以及这些方法的一些组合。
5.现阶段基于机器学习的方法是主流方法，但是该方法的识别性能十分依赖于训练集的质量，而寻找大量合适的高质量的双曲线图像训练集会是一项挑战。在不使用机器学习方法的情况下，现有方法如最小二乘法等，在信噪比较低，杂波较多的情况下，识别性能仍然有待提升。


技术实现要素：

6.有鉴于此，本发明提供了一种基于下开口扫描的双曲线线性化识别方法，能够解决现有方法在双曲线特征不明显、或双曲线重叠的情况下带来的误差较大，从而影响到识别准确性的问题，提高了识别方法的准确性和鲁棒性。
7.本发明的基于下开口扫描的双曲线线性化识别方法，包括：
8.步骤1，对探地雷达b扫描图像去除直达波后，进行二值分割，得到二值图像；
9.步骤2，对二值图像进行下开口扫描处理，得到各下开口区域点簇；
10.步骤3，裁剪各下开口区域点簇所在图像，并按列扫描，提取该下开口区域点簇的点段信息，所述点段信息包括每一列点段的数量、点段的起始点与终止点的纵坐标，以及点段的长度；取每一点段的中点进行hpl处理，得到折线图h；所述hpl处理具体为：
11.s341，构建新的轨道坐标系o
′
x
′y′
，所述新的轨道坐标系为直角坐标系，原点为设置的参考点；
12.s342，获取s33按列提取的点段的中点并利用将点段中点转换到新的轨道坐标系o
′
x
′y′
下，得到折线曲线图；其中，(x,y)为点段中点在原图像中的坐标；(x
′
,y
′
)为点段中点在新的轨道坐标系下的坐标；其中，yd为s33中提取的点段的中点的最大值；r为探地雷达天线的扫描步进距离，v为电磁波在背景介质中的传播速度，t为时间；
13.s343，基于s33获取的点段长度，将s342的折线曲线图扩充为折线图h；
14.步骤4，构建系列相关矩阵si；相关矩阵s的参数包括：滑块宽度、左步进次数、右步进次数以及单次步进距离；其中，滑块宽度取步骤3折线图顶部宽度的1.5～3倍；左步进次数和右步进次数分别取步骤3折线图h中，顶部中点到左、右图像边界的距离列数；单次步进距离在1～最大单次步进距离之间遍历，即i＝1,2,
…
，最大单次步进距离；所述最大单次步进距离为步骤3折线图中顶部中点到图像底部的距离行数；
15.分别计算系列相关矩阵与步骤3折线图的相关值，得到最大相关值；
16.步骤5，若步骤4获取的最大相关值大于或等于判定门限，则判定该折线图对应的点簇为双曲线点簇。
17.较优的，采用自适应阈值实现二值分割，具体包括：
18.s11，对去除直达波后的图像进行中值滤波，得到图像c1；
19.s12，求取图像c1的全局图像的平均值；将图像中所有小于平均值的像素点的值替换为平均值，得到图像c2；
20.s13，利用基于纵向梯度信息的自适应阈值算法，对图像c2进行阈值分割，得到二值图像c3。
21.较优的，s13中，梯度区域集合定义为：
[0022][0023]
其中，(i,j)为像素点坐标；g为设定的最小梯度阈值；
[0024]
定义自适应阈值t为：
[0025]
t＝mean{g(i,j)}
[0026]
利用自适应阈值t对图像c2进行分割，得到二值化图像c3：
[0027][0028]
较优的，采用均值对消的方式去除直达波，即：对每一个像素点，用该像素点的值减去该行的平均值。
[0029]
较优的，所述步骤1中，对二值图像执行先闭后开操作。
[0030]
较优的，下开口区域点簇的获取方法包括：
[0031]
s21，针对二值图像，按行扫描具有向下开口特征的点段，得到下开口区域并建立点簇；
[0032]
s22，对下开口区域进行一次顶部向上查找和左、右尾部向下查找；向上查找过程中，若上一行的点段与本行点段有纵向重叠，则切割掉上一行中超出本行点段横向区域的像素点，将上一行点段中未超出本行点段横向区域的像素点赋予本行点段的簇类号；继续向上查找，直至上方没有纵向重叠点段；向下查找过程中，若下一行的点段与本行点段有纵向重叠，则切割掉下一行中超出本行点段横向区域内侧的像素点，将下一行点段中未超出本行点段横向区域内侧的像素点赋予本行点段的簇类号；继续向下查找，直至下方没有纵向重叠点段，获得最终的本下开口区域点簇。
[0033]
较优的，所述步骤3中裁剪各下开口区域点簇所在图像，具体为：
[0034]
s31，将二值图像分割成多个子图像，每个子图像中仅包含一个下开口区域点簇；对各子图像中的下开口区域点簇进行图像裁剪，得到裁减后的下开口区域图像；
[0035]
s32，对裁减后的下开口区域图像进行纵向双线性插值；
[0036]
s33，对插值后的图像进行列扫描，按列提取点段信息。
[0037]
较优的，所述s31中，图像裁剪具体为：针对各子图像，标定点簇的4个中点：顶部、左尾部、右尾部以及下开口的中点；基于所述4个中点，对点簇区域进行裁剪，得到裁减后的下开口区域图像。
[0038]
较优的，所述步骤4计算相关值时，先计算相关矩阵s与折线图h的相关度，然后减去冗余矩阵的所有元素的冗余值总和，做为该单次步进距离下的相关值；所述冗余矩阵为：折线图h与相关矩阵s做差之后值为1的区域，根据该区域各点与滑块距离赋予该点冗余值，其与区域为零的矩阵。
[0039]
较优的，所述步骤5中，判定门限为：c3×
滑块宽度
×
(左步进次数+右步进次数+1)，其中，c3为常数。
[0040]
有益效果：
[0041]
(1)本发明将双曲线转化为等间距模型，将目标点簇是否为双曲线的问题转化为点簇是否呈现出等间距特性的问题，使得双曲线的待估计参数a、b减少为一个dt/r的同时，还保留了原本双曲线形状的全部信息，与普通拟合方法相比具有更好的准确性。普通拟合方法(如最小二乘法)有可能会因为双曲线形状的特殊性而产生错误的结果，但该方法可以有效应对这样的情况。
[0042]
(2)本发明基于改进的纵向梯度信息实现自适应阈值二值分割，适用于背景噪声区域较大、背景噪声无法完全滤波的情况，能够实现更好的分割效果，且可自适应阈值分割。
[0043]
(3)下开口扫描过程中，为了达到更好的识别效果，本发明对下开口扫描处理做了适当的修改，不再识别上开口以及“x”点段，并且设定了已标记为下开口的双曲线各个位置的点段范围，对超出范围的像素区域进行合理删除，更贴近实际双曲线的特性，增强算法的鲁棒性。
[0044]
(4)本发明通过引入冗余矩阵，将折线图中折线像素超出滑块范围导致的相关性计算误差考虑在内，降低了识别误差。
附图说明
[0045]
图1为均值对消处理结果；(a)为处理前；(b)为处理后。
[0046]
图2为自适应阈值分割效果；(a)基于传统纵向梯度信息的阈值分割方法示意；(b)为传统方法的实际效果；(c)为本发明自适应阈值分割效果。
[0047]
图3为开闭操作处理结果；(a)为闭操作的处理效果；(b)为开操作的处理效果。
[0048]
图4为osca下开口点段示意；(a)为下开口；(b)为上开口。
[0049]
图5为点段交叠的情况。
[0050]
图6为原判别条件带来的影响；(a)为点簇顶端具有毛刺的情况；(b)为传统判别对于(a)的处理结果。
[0051]
图7为本发明下扫描处理的切割效果示意；(a)为点簇顶端具有毛刺的处理示意；(b)为点簇尾部具有毛刺的处理示意。
[0052]
图8为像素聚类示意。
[0053]
图9为osca处理后的效果。
[0054]
图10为探地雷达场景示意。
[0055]
图11为hpl变换示意。
[0056]
图12为变量替换给点段宽度带来的影响。
[0057]
图13为点段形变示意。
[0058]
图14为需要插值的原因。
[0059]
图15为标定四点。
[0060]
图16为转化后的图像。
[0061]
图17为不同步进下相关值。
[0062]
图18为实验结果；(a)～(g)分别为采用本发明方法识别7张lte雷达图像的识别结果；(h)～(n)分别为采用pfj方法识别7张lte雷达图像的识别结果。
[0063]
图19为本发明流程图。
具体实施方式
[0064]
下面结合附图并举实施例，对本发明进行详细描述。
[0065]
本发明提供了一种基于下开口扫描的双曲线线性化识别方法，首先通过预处理方法(均值对消、自适应阈值、开闭操作等)将b扫描图像转化为一张二值图像；然后通过下开口扫描算法(osca)将具有下开口特征的点簇识别出来，同时对每一个像素点进行聚类，在下开口扫描的过程中应尽量剔除同一个聚类在同一列中有两个点段的情况，以便后续处理；之后通过双曲线点段线性化(hpl)方法，将点簇转化为等间距的线性模型，通过计算机表示出理想等间距图像，通过对比评估转化后的图像的等间距特征，判断该点簇与理想双曲线点簇的相似程度。
[0066]
具体步骤如下：
[0067]
步骤一，对探地雷达b扫描图像进行预处理
[0068]
预处理的过程分为三步：第一步，对图像进行均值对消，去除直达波；第二步，对图像进行自适应阈值处理；第三步，对图像进行开闭操作，去除图像的斑点，连接断点。
[0069]
s1，去除直达波：
[0070]
使用均值对消的方法去除直达波，其具体方式为：对每一个像素点(i,j)，用该像素点的值减去该行的平均值，表达式如下：
[0071][0072]
其中，bm×n表示原始b扫描图像，cm×n表示均值对消后的图像。
[0073]
均值对消处理的对比效果如图1所示，可以看到，直达波部分被很好地去除。
[0074]
s2，自适应阈值二值分割算法：
[0075]
可基于纵向梯度信息的自适应阈值算法，即提取图像中纵向梯度较大的区域，对该区域求取平均值作为图像的阈值t，进行二值分割，如图2(a)所示。但是在背景噪声区域比较大的情况下，由于背景噪声的值往往在正向与负向峰值之间，如果仅对图像中梯度较大的区域取平均，那么背景噪声便无法完全滤除，如图2(b)所示，如果要实现更好的分割效果，需要将图像的负峰去除，因此本实施例做出以下改进：
[0076]
本实施例自适应阈值算法部分由三个环节组成：
[0077]
①
使用中值滤波去除图像噪点，减少平坦区域的梯度值，该步骤是为了避免图像噪点对总体阈值产生影响，经过该步骤得到图像c1；
[0078]
②
求图像c1的全局图像的平均值m，对每一个像素点，若该点的值小于平均值，则令该点的值等于m，该步骤是为了避免负峰区域拉低总体阈值。用mean{}表示取矩阵平均值，用数学方式可表达为：
[0079]
m＝mean{c1(i,j)}
[0080][0081]
③
基于纵向梯度信息的自适应阈值算法进行二值分割；具体的，梯度区域集合定义为：
[0082][0083]
其中g为最小梯度阈值。
[0084]
定义自适应阈值t：
[0085]
t＝mean{g(i,j)}
[0086]
经过阈值t的分割得到一幅二值化图像c3：
[0087][0088]
经过自适应阈值方法处理后的效果如图2(c)所示，经由阈值分割过后，原本的灰度图像转化为一幅二值图像。
[0089]
s3，开-闭操作：
[0090]
开-闭操作分为两步——开操作与闭操作，它们都是形态学操作，表示扩张与侵蚀。假定一结构元dk，开操作定义为：
[0091][0092]
同样的，闭操作定义为：
[0093]
[0094]
表示结构元对图像边缘的侵蚀，侵蚀作用的大小取决于结构元半径的选择。同理，
⊕
表示结构元对图像边缘的扩张。本方案采取先闭后开操作。闭操作的结构元半径设为3，开操作的结构元半径设为2。效果如图3所示。图3(a)图显示了闭操作的处理效果，图3(b)图显示了开操作的处理效果。
[0095]
步骤二，对预处理后的图像进行改进的下开口扫描处理(osca)
[0096]
传统的下开口扫描处理(osca)是通过按行扫描具有向下开口特征的点段来判别下开口点簇的方法。其中，下开口的概念定义为：在第i行的某一个点段(点段1)的下方有两个与该点段在纵向上重叠的点段，分别为点段2与点段3，并且点段1的起点在点段2的起点与终点之间，点段1的终点在点段3的起点与终点之间，三个点段互相连通，如图4(a)所示。在图4(a)中，点段2与点段3存在间隙且每一个点段长度都大于2，这样点段1、点段2与点段3便形成了一处下开口区域。同样的，如果点段2与点段3位于点段1的上方临近一行，三个点段纵向重叠且点段2与点段3之间存在间隙，则这三个点段形成上开口区域，如图4(b)所示。计算机识别出下开口区域后，便为该区域赋予一簇类号，建立一个新的点簇。其中，开口以上的部分称为双曲线的顶部，以下的部分为双曲线的尾部。程序搜索到下开口区域之后会进行一次向上查找与两次向下查找，在向上查找的过程中，若上一行的点段不超出下一行临近点段起止范围，则为该点段的所有像素赋予下方点段的簇类号，并继续向上扫描，直到上方没有重叠点段或重叠的点段超出下方点段所对应的横向的区域为止。两次向下扫描分别沿着点段2与点段3搜索点簇左尾部和右尾部的边界，即，下一行点段存在与上一行临近点段的重叠段，则将该下一行点段的所有像素赋予上方点段的簇类号，并继续向下扫描，直到下方没有重叠点段为止。此外，在将点段加入到已有簇类中时，有时会遇到两个点簇交叠的情况，如图5所示。在遇到两点段交叠的情况时，传统方法采用了将重叠部分删除的策略。
[0097]
本发明为了与后面的步骤衔接，实现更好的准确性，对传统的下开口扫描处理进行了如下改进：
[0098]
其一，对点段超出合理范围的区域做切割
[0099]
在实际图像中，双曲线点簇的顶部区域并不完全遵循上方的点段落在下方临近点段起止区间内的条件，如图6(a)所示，实际图像中双曲线点簇顶部可能会出现微小的毛刺，该条件的不满足会导致上方的所有区域被从原点簇中切除掉，从而严重影响到双曲线的形态。如图6(b)所示，因此本发明做出如下修改：
[0100]
假设点段x被识别为下开口簇的一部分，点段区间为[xb,xe]，且上方存在一点段在纵向上与x存在重叠区域的点段y，点段区间为[yb,ye]，则取区间[yb′
,ye′
]内的像素放入点段x所在的点簇中，其中：
[0101]
y'b＝max(xb,yb)
[0102]
y'e＝min(xe,ye)
[0103]
其中max()和min()分别表示取最大值及最小值。在搜索左右尾部的过程中进行同样的处理，在向下扫描的过程中，对左尾部下方临近点段的终点与右尾部下方临近点段的起点分别做切割。假设左尾部与右尾部的点段为l与r，且分别与下方点段ll与rr有交叠，l与r的横向区间分别为[lb,le]，[rb,re]；ll与rr为[llb,lle]，[rrb,rre]；则将区间[llb,ll
′e]，[rr
′b,rre]内的像素放入l与r所在的点簇中，其中
[0104]
rrb'＝max(rb,rrb)
[0105]
ll'e＝min(le,lle)
[0106]
切割效果如图7所示。切割所带来的影响远不足以影响程序的判断，但是可以大大改善双曲线点簇毛刺过多的情况下的聚类效果。
[0107]
其二，在基于一的聚类条件下，允许交叉区域的像素被添加到多个点簇中，如图8(a)所示。
[0108]
图中的交叉区域视为多个点簇的公共区域，在聚类的过程中，公共区域分别成为所有组成该区域的点簇的组成部分，该区域的像素点包含所有组成该区域点簇的类别号，其在不同的点簇中都会出现。聚类结果如图8(b)、(c)所示。
[0109]
经过下开口扫描处理后的图像如图9所示，可以看到，非向下开口的点簇已被过滤。
[0110]
步骤三，将双曲线点簇转化为等间距点阵
[0111]
在探地雷达中，假设目标所在的横坐标为0，如图10所示，则天线到目标的距离与天线横坐标之间的关系可以表示为：
[0112]
s2＝d2+(r
·
x)2[0113]
r表示天线每一步的步进距离，x表示天线与目标之间间隔的轨道数，d表示目标深度，s表示天线到目标的直线距离，则该式可以化为：
[0114][0115]
在时间轴上进行采样，已知在x处，目标在图像中出现的时间为y，由：
[0116][0117]
代入可以得到：
[0118][0119]
参数yd表示天线位于目标正上方时，回波在图像中出现的时间，即双曲线顶点的纵坐标，如图11中所示。v表示电磁波在背景介质中的传播速度，参数yd已知，由式中可以看到，只需要估计采样时间与速度之积和天线步进的比值v*dt/r，便可唯一确定在图像中某一双曲线区域所对应的双曲线方程。依此，将图像中经过osca处理后的双曲线点簇与理想的双曲线直接进行相关操作，便可以得到双曲线点簇与理想双曲线的相似程度，从而对双曲线进行判断。这也是十分贴近人眼视觉判断的方法。
[0120]
但在实际图像处理中，由于只有特定的雷达采样的文件中才能得到雷达在采样中的步进距离，我们需要一种适用范围广的方法去对双曲线进行判断。本发明将直接从图像处理的角度，从双曲线的特性出发对双曲线点簇进行判断。
[0121]
如图11所示，图中有一参考点与一条可供点段滑动的轨道，如果将双曲线每一列的点段放置于轨道上，同时将点段的纵坐标转化为该点段到参考点的直线距离，通过这种方式得到一个点段集合。由于双曲线是雷达天线在等间距采样下得到的图像，因此将双曲
线的点段通过这种方式重新排列，所有点段在轨道上的距离应呈现出等间距的特征。
[0122]
这也可以通过数学的方式进行解释：
[0123]
如果令则有：
[0124][0125]
化简可得：
[0126][0127]
由此可以看出，y’与x存在线性关系，y’相当于点段在“滑轨”上距离参考点的坐标，如图11中所示，y
’‑
x图像的斜率的绝对值便是参数2r/vdt，因此原有问题便转化为y
’‑
x是否满足一种线性关系。通过该算法避免了在双曲线中进行参数估计的问题。因此，将y-x图像转化为y
’‑
x图像，通过设计评估标准，来判断y’与x是否满足线性关系，便可以得出点簇是否满足双曲线点簇的特质，从而可以避免在双曲线方程中对参数的估计所产生的误差。
[0128]
在实际图像中，点段在纵向上包含一定的宽度，在确定点段在“滑轨”上坐标范围的过程中，需要考虑此变换给点段宽度带来的影响。由图12所示，根据三角形的特性，变换后点段长度总是会大于原点段长度(y
2-y1)，由此可以得到，在转化过程中每一个点段长度都会发生不同程度的扭曲，即原本相同长度的点段，可能会因为位置不同而在转化后变得长度不同。
[0129]
因此在实际构造等间距点段的过程中，必须先提取点段中的某一点进行坐标转换，之后再根据该点在点段中的位置以及点段长度重新构造点段，如图13(b)所示，而不是将点段中的所有点都按照如上方式转化到新的坐标轴上，这样会导致图像的失真。在本发明中，采取提取双曲线中点并进行处理的方式，根据重新排布的中点安排点段位置，将双曲线转化为等间距型折线。
[0130]
假设某一点段中点的坐标为ym，点段长度为l，则点段在“滑轨”上的坐标范围为其中cell()表示向下取整。通过以上方法将双曲线每一列的点段等间距化的过程，本发明中称为双曲线点段线性化(hpl)方法。
[0131]
由于在数字图像中，点段的位移只能取整数，而点段间隔2r/dt有可能是小数，在一些情况下(如图14所示)会导致图像与等间距滑块的相关值较低，因此需要将原图像进行纵向插值处理。
[0132]
图14(a)中是滑块步进为1情况下滑块与实际点段的位置，图14(b)中是滑块步进为2情况下滑块与实际点段的位置，白色区域为实际点段的位置，黑色区域为滑块的位置。可以看到实际点段的间隔应在1与2之间，而滑块间隔只能取整数，在该情况下进行相关操作，会使得相关值较低，并可能导致误判。
[0133]
因此步骤3具体流程如下：
[0134]
依照点簇聚类号对图像进行分割，使得一张图中只有一个类别，此时进行一次扫
描，标定点簇的顶部，左尾部，右尾部，下开口中点四个点，用这四个点确定裁减后图像在原本图像中的位置，如图15所示；
[0135]
假设四个点的坐标分别为：vu＝[v
uc
,v
ur
]；vd＝[v
dc
,v
dr
]；v
l
＝[v
lc
,v
lr
]；vr＝[v
rc
,v
rr
]。
[0136]
根据这四个点切割图像，得到裁减后的双曲线区域；裁减后图像的尺寸为:(v
rc-v
lc
+1)
×
(max(v
rr
,v
lr
)-v
ur
+1)；
[0137]
对裁减后图像进行纵向双线性插值；
[0138]
进行一次列扫描，获取按列分布的列点段信息，包括每一列点段的数量以及点段的长度，点段的起始与终止点的纵坐标等；
[0139]
取每一点段中点进行hpl处理；
[0140]
通过以上步骤，便可以得到一幅由双曲线转化而成的折线图，得到的图像如图16所示，假设该图像为h
mh
×
nh
。
[0141]
步骤四，评估转化后图像的等间距性，判断原图是否为双曲线
[0142]
将滑块在每一列的位置进行拼接，可以得到一个倒“v”型图像，我们称为相关矩阵，相关矩阵具有四个参数，滑块宽度、左步进次数、右步进次数以及单次步进距离，滑块宽度取决于顶点所在点段的宽度；在实际图像中，很多肉眼可识别的双曲线并不是完全对称的，因此在设计滑块滑动次数时，左步进次数与右步进次数应分开设计，左步进次数取决于顶点与左尾部的横向距离；右步进次数取决于顶点与右尾部的横向距离。假设顶点位于该列第k个点段序列上，滑块宽度为z，左右步进次数分别为w
l
，wr：
[0143]
z＝c1·
length(p(v
dc
,k)|v
dc
∈p(v
dc
,k))
[0144]wl
＝v
dc-v
lc
[0145]
wr＝v
rc-v
dc
[0146]
其中p(i,j)表示图像第i列的从上到下第j个点段序列，c1为常数
[0147]
单次步进距离在不超出裁减插值后的范围内遍历，取最大相关值对应的滑块步进，即为所选取步进，表示为：
[0148][0149]
max(r
·wl
,r
·
wr)≤mh[0150]
argmax(cor(r))表示cor(r)取最大值时r的值，cor(r)表示在单次步进距离为r个像素时，相关矩阵与h矩阵的相关值，关于相关值的具体计算会在下一步中说明。
[0151]
例如以下是在不同步进下，矩阵h与相关矩阵的相关程度，对应图17：
[0152]
不同步进下的相关值：
[0153]
step1234567correlate-13563018690311581182642015320051
[0154]
在滑块单次步进距离为6像素时，相关矩阵与矩阵h的相关程度最高，因此本例滑块步进r取6。
[0155]
当h中的像素超出了滑块范围时，与相关矩阵的相关性应该减少，因此，我们引入冗余矩阵。冗余矩阵表示原矩阵h在相关矩阵区域外的像素点分布，即将矩阵h与相关矩阵相减，相减后矩阵的值为一的像素区域。对冗余矩阵中所有值为1的像素点，需要评估该点
的冗余程度，本发明中称为冗余值，冗余值的计算应遵循以下原则：
[0156]
1.冗余像素与本轨道滑块距离较近时，需要视为由噪声带来的误差，冗余值应尽量小；
[0157]
2.冗余像素超过了一定范围后，视觉上已经超出了正常双曲线的区域，足够影响对矩阵h等间距性的判断，因此需要有较大的冗余值；
[0158]
3.计算冗余值的函数理应连续；
[0159]
结合以上条件，定义冗余值计算方法为
[0160]
r(x,y)＝l2/λ
[0161]
λ＝c2z
[0162]
l＝min(|y-r
×
|x-(v
dc-v
lc
+1)|-v
dr
|,|y-r
×
|x-(v
dc-v
lc
+1)|-v
dr
+z|)
[0163]
λ为缓冲区，为滑块宽度的c2倍，c2为常数；l表示该像素点距离所在行滑块(上边界或下边界)的最短距离，z为滑块宽度。l对整体相关值的影响成几何增长。
[0164]
将相关矩阵s与矩阵h做相关，再减去冗余矩阵的所有元素的冗余值总和，即为该步进下的相关值，表示为：
[0165]
corr(step)＝tr(s
·ht
)-sum(r)
[0166]
cor＝max(corr(step))
[0167]
设置判定门限为c3×
滑块宽度
×
(左步进次数+右步进次数+1)，c3为常数，表示为：
[0168]
t＝c3z(w
l
+wr+1)
[0169]
如果cor＞t，则判定该点簇为双曲线点簇。
[0170]
实施例
[0171]
为了验证本发明对地下目标双曲线识别的准确性，本实验采用七张lte雷达数据作为输入，雷达数据的步进r与dt皆视为未知。本实验以识别准确率与召回率作为性能指标，并与pfj方法做对比，采用图像皆为图像质量不佳情况下的雷达探测结果，实验结果如图十八所示
[0172]
本方案识别结果：
[0173] 图像1图像2图像3图像4图像5图像6图像7双曲线个数1519252210识别出的双曲线个数1314262411双曲线被正确识别出的个数131425229召回率13/1514/1911119/10准确率1115/611/29/11
[0174]
pfj方法：
[0175]
pfj(parabolic fitting-based judgment method)，是将点簇通过拟合成抛物线并加以判别的方法，判别标准如下：
[0176]
a)拟合的双曲线需要满足开口向下
[0177]
b)拟合的双曲线需要满足对称轴位于下开口点段的横向坐标范围内
[0178]
以下是pfj方法识别结果：
[0179]
[0180][0181]
(注：若将多个双曲线的连接体识别为一个双曲线，则视为错误识别)
[0182]
通过对比可以看到，在识别低信噪比的双曲线的情况下，本发明中的方法具有明显的优越性。
[0183]
综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

技术特征：
1.一种基于下开口扫描的双曲线线性化识别方法，其特征在于，包括：步骤1，对探地雷达b扫描图像去除直达波后，进行二值分割，得到二值图像；步骤2，对二值图像进行下开口扫描处理，得到各下开口区域点簇；步骤3，裁剪各下开口区域点簇所在图像，并按列扫描，提取该下开口区域点簇的点段信息，所述点段信息包括每一列点段的数量、点段的起始点与终止点的纵坐标，以及点段的长度；取每一点段的中点进行hpl处理，得到折线图h；所述hpl处理具体为：s341，构建新的轨道坐标系o
′
x
′
y
′
，所述新的轨道坐标系为直角坐标系，原点为设置的参考点；s342，获取s33按列提取的点段的中点并利用将点段中点转换到新的轨道坐标系o
′
x
′
y
′
下，得到折线曲线图；其中，(x,y)为点段中点在原图像中的坐标；(x
′
,y
′
)为点段中点在新的轨道坐标系下的坐标；其中，y
d
为s33中提取的点段的中点的最大值；r为探地雷达天线的扫描步进距离，v为电磁波在背景介质中的传播速度，t为时间；s343，基于s33获取的点段长度，将s342的折线曲线图扩充为折线图h；步骤4，构建系列相关矩阵s
i
；相关矩阵s的参数包括：滑块宽度、左步进次数、右步进次数以及单次步进距离；其中，滑块宽度取步骤3折线图顶部宽度的1.5～3倍；左步进次数和右步进次数分别取步骤3折线图h中，顶部中点到左、右图像边界的距离列数；单次步进距离在1～最大单次步进距离之间遍历，即i＝1,2,
…
，最大单次步进距离；所述最大单次步进距离为步骤3折线图中顶部中点到图像底部的距离行数；分别计算系列相关矩阵与步骤3折线图的相关值，得到最大相关值；步骤5，若步骤4获取的最大相关值大于或等于判定门限，则判定该折线图对应的点簇为双曲线点簇。2.如权利要求1所述的方法，其特征在于，采用自适应阈值实现二值分割，具体包括：s11，对去除直达波后的图像进行中值滤波，得到图像c1；s12，求取图像c1的全局图像的平均值；将图像中所有小于平均值的像素点的值替换为平均值，得到图像c2；s13，利用基于纵向梯度信息的自适应阈值算法，对图像c2进行阈值分割，得到二值图像c3。3.如权利要求2所述的方法，其特征在于，s13中，梯度区域集合定义为：其中，(i,j)为像素点坐标；g为设定的最小梯度阈值；定义自适应阈值t为：t＝mean{g(i,j)}利用自适应阈值t对图像c2进行分割，得到二值化图像c3：
4.如权利要求1所述的方法，其特征在于，采用均值对消的方式去除直达波，即：对每一个像素点，用该像素点的值减去该行的平均值。5.如权利要求1～4任一所述的方法，其特征在于，所述步骤1中，对二值图像执行先闭后开操作。6.如权利要求1所述的方法，其特征在于，下开口区域点簇的获取方法包括：s21，针对二值图像，按行扫描具有向下开口特征的点段，得到下开口区域并建立点簇；s22，对下开口区域进行一次顶部向上查找和左、右尾部向下查找；向上查找过程中，若上一行的点段与本行点段有纵向重叠，则切割掉上一行中超出本行点段横向区域的像素点，将上一行点段中未超出本行点段横向区域的像素点赋予本行点段的簇类号；继续向上查找，直至上方没有纵向重叠点段；向下查找过程中，若下一行的点段与本行点段有纵向重叠，则切割掉下一行中超出本行点段横向区域内侧的像素点，将下一行点段中未超出本行点段横向区域内侧的像素点赋予本行点段的簇类号；继续向下查找，直至下方没有纵向重叠点段，获得最终的本下开口区域点簇。7.如权利要求1所述的方法，其特征在于，所述步骤3中裁剪各下开口区域点簇所在图像，具体为：s31，将二值图像分割成多个子图像，每个子图像中仅包含一个下开口区域点簇；对各子图像中的下开口区域点簇进行图像裁剪，得到裁减后的下开口区域图像；s32，对裁减后的下开口区域图像进行纵向双线性插值；s33，对插值后的图像进行列扫描，按列提取点段信息。8.如权利要求1所述的方法，其特征在于，所述s31中，图像裁剪具体为：针对各子图像，标定点簇的4个中点：顶部、左尾部、右尾部以及下开口的中点；基于所述4个中点，对点簇区域进行裁剪，得到裁减后的下开口区域图像。9.如权利要求1所述的方法，其特征在于，所述步骤4计算相关值时，先计算相关矩阵s与折线图h的相关度，然后减去冗余矩阵的所有元素的冗余值总和，做为该单次步进距离下的相关值；所述冗余矩阵为：折线图h与相关矩阵s做差之后值为1的区域，根据该区域各点与滑块距离赋予该点冗余值，其与区域为零的矩阵。10.如权利要求1或9所述的方法，其特征在于，所述步骤5中，判定门限为：c3×
滑块宽度
×
(左步进次数+右步进次数+1)，其中，c3为常数。

技术总结
本发明公开了一种基于下开口扫描的双曲线线性化识别方法。本发明将双曲线转化为等间距模型，将目标点簇是否为双曲线的问题转化为点簇是否呈现出等间距特性的问题，使得双曲线的待估计参数减少为一个的同时，还保留了原本双曲线形状的全部信息，与传统的最小二乘法相比具有更好的准确性。比具有更好的准确性。比具有更好的准确性。


技术研发人员：兰天 陈宏畅 杨小鹏 赵毅 李洪洁 渠晓东 钟世超 曾小路
受保护的技术使用者：北京理工大学
技术研发日：2023.06.05
技术公布日：2023/8/31