本发明涉及钢筋网质量检测领域，尤其是一种基于sar的gpr多层钢筋网参数检测方法。
背景技术：
探地雷达(gpr)最早起源于1904年德国人hulsmeyer用电磁波信号探测地下的金属物体。1910年，letmbach和lowy在其专利中提出了利用电磁波对地下掩埋物体的探测。到了1926年hulsenbech第一次使用脉冲波技术探测埋藏介质的结构，并发现了介电常数不同的介质交界面会产生电磁波反射，成为了探地雷达研究的基本理论依据。1929年stern将探地雷达应用到实际中，用于测量冰川的厚度。此后，探地雷达主要应用于深层探测。随着现代数字信号处理技术的发展和应用，70年代后，探地雷达从冰层、盐矿等弱耗介质应用扩展到土层、煤层、岩层等有耗介质，出现了用于探测地下坑道的宽带雷达，此后这种雷达快速发展，在中浅层目标探测中得到了广泛的应用。随着探地雷达技术的高速发展和广泛应用，高分辨率数据处理和目标成像技术成为了探地雷达技术发展的关键。近十多年来合成孔径成像技术在探地雷达中逐步得到了应用，同时由于雷达波和地震波在传播过程中存在相似性，探地雷达数据处理和成像可采用地震学中的波动偏移成像技术。
在钢筋网结构越来越普及的今天，钢筋网自身的质量问题引起的病害层出不穷。由于钢筋本身具有易腐蚀的缺点，导致建筑中经常出现钢筋损坏、断裂的病害；又由于在建筑的不同部位对钢筋强度及规格的要求也不同，错误的使用钢筋将导致严重的质量问题以及安全隐患。传统的钢筋检测方法为电流法(给钢筋通电并通过检测电流的大小判断钢筋的腐蚀程度等属性)，但这种方法检测速度极为缓慢而且无法大面积使用。探地雷达在国内外都被广泛应用于地质勘探、工程检测等诸多领域，并且已经被证明是一种行之有效的无损探测方法。相对于传统方法，雷达探测的成像效果以及探测精度都有明显进步，但是目前对钢筋网的探测也仅仅局限于单层、结构较为简单的钢筋网中，对于目前普遍采用的多层钢筋网结构的探测，由于其结构的复杂性，图像往往由于干扰较多而只能手动进行感兴趣区域的抓取和病害判断。由于钢筋网在实际应用中往往铺设面积巨大，例如，高铁所采用的无砟轨道大量使用钢筋混凝土结构，而高铁轨道动辄上千公里，探测数据可达几百gb，在这雷达数据量往往十分庞大的客观条件下，人工数据操作很难在实际应用中进行。随着钢筋网结构的大范围使用，对于多层钢筋网的无损检测和自动处理的技术需求越来越迫切。在此方面，目前许多学者做了大量的研究，例如，利用曲线拟合技术进行特定图像特征的识别，这种方法运算量十分巨大；运用hog(梯度直方图)进行图像特征的提取并用svm(支持向量机)等分类器进行特征分类，或对图像进行小波变换后用分类器进行分类。这些算法集中在特征提取并使用分类器分类的模式，这种方法的最大缺点在于必须事先知道各种物体或病害的图像模式，而这种固定的模式在一些复杂情况下并不存在，并且需要大量实测并确认正确的数据进行模型训练，这些条件都极难实现。
技术实现要素：
发明目的：为克服现有探测技术需要人工参数计算而速度极为缓慢的问题，实现参数自动提取，本发明提出一种基于sar的gpr多层钢筋网参数检测方法。
技术方案：本发明提出的技术方案为：
一种基于sar的gpr多层钢筋网参数检测方法，包括步骤：
(1)获取多层钢筋网的gpr扫描图像，对gpr扫描图像进行预处理，并采用空间域中值滤波法抑制扫描图像数据中的椒盐噪声；
(2)对经过预处理后的gpr扫描图像进行频率波数偏移成像，得到偏移图像；
(3)从偏移图像中自顶层向下逐层提取钢筋网的参数，提取下层钢筋网参数时，将偏移图像中的上层钢筋网截断，提取包括步骤：
(3-1)对偏移图像做像素点灰度值统计，得到每一个灰度值对应的像素点数量占总像素点数量的比例值；
(3-2)将得到的比例值按照其对应的灰度值升序排列形成比例值序列；从比例值序列的最小灰度值一端开始逐个累加比例值，直至比例值的累加和超过设定的阈值t，选取最后累加的一个比例值，将其对应的灰度值作为全局阈值对偏移图像进行二值化处理，得到二值图像；
(3-3)从二值图像中提取钢筋位置坐标，根据钢筋位置坐标计算钢筋水平方向间隔和钢筋深度。
进一步的，所述对于gpr扫描图像进行预处理的方法为：
对gpr扫描图像数据二维矩阵中的每行数据求方差，选出方差小于预设阈值σ的行，并将选出的行中的每个元素减去所在行的平均值。
进一步的，所述对经过预处理后的gpr扫描图像进行频率波数偏移成像的步骤包括：
1)定义预处理后的gpr扫描图像数据为u(x,z＝0,t)，其中，x表示雷达波在水平方向的位置坐标，z表示深度，t为时刻；
2)对u(x,z＝0,t)在x和t方向做二维傅里叶变换：
式中，kx表示雷达波的波数矢量在水平方向的分量；
3)通过差值计算，将u(kx,z＝0,ω)从(kx,ω)空间映射到(kx,kz)空间，得到偏移图像，偏移图像表示为：
式中，v表示电磁波在介质中的波速，kz表示雷达波的波数矢量在垂直方向的分量；b(kx,kz)为中间函数，b(kx,kz)的表达式为：
进一步的，所述步骤(3-3)中计算钢筋水平方向间隔的方法为：
获取二值图像中水平方向上相邻钢筋之间的迹线数，用迹线数乘以相邻迹线之间的间隔距离，得到钢筋水平方向间隔。
进一步的，所述步骤(3-3)中计算钢筋深度的方法为：
获取钢筋位置坐标(x*,y*)，式中，x*为钢筋水平方向位置坐标，y*为钢筋垂直方向位置坐标；钢筋深度的计算公式为：
h＝y*×δt×v
式中，h表示钢筋深度，δt表示采样间隔。
有益效果：与现有技术相比，本发明具有以下优势：
本发明利用偏移成像技术有效提高了gpr图像分辨率以及可识别度。并且计算过程中只使用了gpr探测系统可以轻易得到的先验知识，最终达到了对钢筋网的自动参数提取，大大提高了探测效率。
附图说明
图1为原始gpr图像；
图2为原始gpr图像的偏移图像。
具体实施方式
下面结合附图对本发明作更进一步的说明。
本发明共分为三部分：一、预处理；二、2-d偏移成像；三、图像参数的自动提取。下面对这三个部分分别予以说明。
1)数据预处理
在探地雷达数据中，由雷达波在不同界面引起的反射波或者环境中的噪声在一维a扫描中总是会出现在同一时间位置，这些杂波会在二维b扫描图像中产生一条水平线，这条水平线有可能会掩盖目标位置。另外，发射天线的电磁波经介质表面反射或直接被雷达接收天线接收的直达波能量最大，会叠加在目标有效信号上，特别是浅地层目标回波与其到达时间很接近，就会造成能量较小的目标回波被淹没，从而形成强烈的直达波干扰，不利于目标的检测。
现有的抑制直达波的方法有均值法滤波、小波变换法和自适应抵消法等，其中最简单有效的是均值法。均值法抑制杂波就是将b扫描数据二维矩阵每一行的每个元素减去所在行的平均值。但这种方法在去掉背景噪声的同时也明显削弱了钢筋反射波强度，因此，本发明对数据的抑制直达波的方法作了改进，即先对b扫描数据二维矩阵中的每行数据求方差，然后只对方差极小(接近于0)的行进行上述去均值的操作，保证钢筋反射波的强度不受影响。此外，空间域中值滤波也被用来消除椒盐噪声。
2)2-d偏移成像
首先介绍孔径合成技术基本原理：
探地雷达合成孔径成像基于爆炸反射模型(explodingreflectormodel，简称erm)，主要就是对b扫描图像进行合成孔径处理，使其更接近于目标的真实尺寸和位置，有效地提高探地雷达图像的分辨率，提高图像的可解释性，从而有利于对目标进行精确定位。
在gpr系统中使用的基于波动方程偏移的合成孔径算法主要包括频率波数偏移法(fkmigration)，相移偏移法(phaseshift)，基尔霍夫积分偏移法和有限差分偏移法。频率波数偏移算法要求速度为常数，这对于浅地层环境来说满足条件，且采用快速傅里叶变换(fft)，算法速度和效率远高于其他算法，所以在实际工程应用中多采用频率波数偏移法。
fk偏移假设在介质中的波速为常数。雷达波场传播时满足如下标量赫姆霍兹波动方程：
其中，x表示雷达波沿雷达扫描方向(即水平方向)的坐标，z表示深度方向(即垂直方向)，t为时刻，v表示电磁波在介质中的波速。
对上式进行傅里叶变换可得到其频率-波数(ω,k)域描述：
式中，w是频率，k＝w/v是波数矢量的幅度，kx和kz分别表示波数矢量在水平和垂直方向的分量。
由于发射源到反射点和从反射点到接收点的路径重叠，电磁波总传播时间2倍于单程的从发射源到反射点的时间，所以可以在偏移中把波速看作波在实际介质中的一半。
对雷达接收数据在x和t方向做傅里叶变换：
对上式作深度方向扩展可得到：
u(kx,z,ω)关于kx和ω的逆傅里叶变换如下：
由(3)式可得到：
推导可得：
上式u(x,z,t＝0)即为偏移后的图像，令：
则公式(7)可化为：
由公式(9)可得到目标源的波场，那么偏移成像就是计算t＝0时的一个二维逆傅里叶变换，即实现了gpr数据的偏移成像。图1为原始钢筋网gpr图像，图中钢筋反射波呈现相互交叉影响的双曲线形；图2为偏移成像后的gpr图像，可以看到双曲线经过偏移后被汇聚到接近于圆形的较小区域中，极大增加了图像可识别度，有利于后续的钢筋网参数自动提取。
3)钢筋参数自动提取
经过偏移之后的gpr图像中，钢筋网反射的双曲线被汇聚到一个互不干扰的较小区域，使用图像分割技术可以较好的提取出钢筋。考虑到运算量的制约，这里使用自适应阈值化进行图像分割。但能根据不同图像自动选取自适应阈值的方法是个难题。
通过观察和大量数据统计，我们发现了两个现象：1)在钢筋网探测图像中，钢筋反射波的强度总是远大于其他物体反射波；2)在偏移后的钢筋网gpr图像中，钢筋区域的总面积和gpr图像总面积的比例总是在一个固定的范围内，大约为0.5％～1％之间，且不管gpr参数(时窗，采样率等)变化与否。
利用这些结论，我们提出了基于偏移成像算法的自适应阈值化方法，具体步骤为：
从偏移图像中自顶层向下逐层提取钢筋网的参数，提取下层钢筋网参数时，将偏移图像中的上层钢筋网截断，提取包括步骤：
(1)对偏移图像做像素点灰度值统计，得到每一个灰度值对应的像素点数量占总像素点数量的比例值；
(2)将得到的比例值按照其对应的灰度值升序排列形成比例值序列；从比例值序列的最小灰度值(也可能是最大灰度值，这和雷达信号相位有关)一端开始逐个累加比例值，直至比例值的累加和超过设定的阈值t(t的取值在0.5％～1％之间)，选取最后累加的一个比例值，将其对应的灰度值作为全局阈值对偏移图像进行二值化处理，得到二值图像；
(3)从二值图像中提取钢筋位置坐标，根据钢筋位置坐标计算钢筋水平方向间隔和钢筋深度。
如果是双层钢筋，由于上层钢筋反射强度巨大，导致二值化后只得到上层钢筋位置，下层很难被分割出来。这种情况我们可以先对上层钢筋进行间隔和深度的参数计算，得到结果后将上层钢筋截断(比如数据置零)，再用相同的自适应阈值化方法提取下层钢筋的参数。
钢筋水平方向间隔计算较为简单，首先得到二值图像中钢筋区域灰度值最值(最大值或最小值，从最大值还是最小值开始和雷达信号的相位有关)的位置(x,y)，这里(x,y)就是钢筋的坐标，进而得到水平方向上相邻钢筋之间的迹线数n，钢筋水平间隔等于n*dx(dx为每条迹线之间的间隔，是测量时即可得到的radar参数)。
计算钢筋深度也只需要钢筋坐标(x,y)，钢筋深度的计算公式为：
h＝y*×δt×v
式中，h表示钢筋深度，δt表示采样间隔，v为波速，v＝c/sqrt(a),a为介电常数。
但由于gpr系统发射信号时有一固定延迟(不同的雷达发射延迟也不同)，这就导致计算的电磁波传播时间比实际情况要大。gpr雷达直达波为接收端第一次接收到的信号，我们可以根据发射端和接收端的距离以及直达波到达的时间来计算发射延迟，在计算深度时减去发射延迟后再乘以波速才可得到深度。
本发明利用偏移成像技术有效提高了gpr图像分辨率以及可识别度。并且计算过程中只使用了gpr探测系统即可轻易得到的先验知识，最终达到了对钢筋网的自动参数提取，大大提高了探测效率。
以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。