本发明涉及探地雷达图像处理方法，用于地下柱状反射体倾角大小的检测，尤其是对探测场景进行数学模型构建，得到柱状反射体双曲线信号对应的函数方程，并进行曲线拟合的过程，属于检测与图像处理领域。
背景技术：
探地雷达是利用电磁波探测地下物质性质与分布规律的地球物理技术，其发射天线向地下发射一定频率的电磁波，并被地下目标物反射至接收天线，通过分析所接收到电磁波的振幅与形状等信息，即可推知地下目标物的物理特性。
地下管线是保障城市平稳运行的重要基础设施，负责实现城市内部的供水、供气等功能，城市的快速发展和承灾能力与地下管线体系的完善性密切相关，故寻找一种准确快速的地下管线检测方法至关重要。目前，地下管线检测的主要方法有充电法、电磁感应法等。探地雷达作为一种有效的探测手段，同样广泛应用于地下管线探测中，相比于其它方法而言，探地雷达能够完成非金属管线的探测，并具有无损探测、对地下管线干扰小的优点。
在现有的研究中，通常利用探地雷达对地下管线的直径与深度进行估计。就常用的收发共置型探地雷达而言，电磁波从发射天线发出，到达地下目标物，被目标物反射后再被接收天线所接收，即在雷达图像上形成反射信号，该反射信号所对应的是电磁波从发射天线到目标物再回到接收天线所消耗的时间，与雷达装置在测线上位置的函数关系。利用这一函数关系，雷达图像上的信号即可用数学方程式表达出来。由于方程式中包含诸如管线直径、管线深度等未知参数，故可采用曲线拟合的方法求出这些未知参数的估计值。
然而，到目前为止，利用探地雷达对地下管线倾角的检测却很少被关注。了解地下管线的倾斜角度，可以有效获取管线系统中存在的异常状况，有助于专业人员对地下管线进行更加精确的检查与维护。以柱状反射体为例，考虑到其在空间中的倾斜角度对其雷达信号函数解析式的影响，柱状反射体倾角会影响其双曲线反射信号的特征，致使很难利用雷达信号对其进行精确定位。故在对地下管线进行检测时，必须对其倾角信息加以了解，才能更加精准地获取地下管线的其它信息。
技术实现要素：
为实现本发明之目的，采用以下技术方案予以实现：
一种检测地下柱状反射体倾角大小的方法，其中，所述方法包括如下步骤：
步骤a，获取地下柱状物的雷达扫描图像；
步骤b，从步骤a所得到的雷达扫描图像中，在柱状物的双曲线反射信号上提取多个采样点坐标；
步骤c，建立数学模型，利用相应算法进行曲线拟合，计算得到柱状物倾角的测量值。
所述的方法，其中：步骤a中，获取地下柱状物的雷达扫描图像的具体方式是：使用收发天线一体的便携式雷达系统，通过布设多条测线，对地下柱状反射体进行探测。
所述的方法，其中：步骤b中，在柱状物的双曲线反射信号上提取每个采样点坐标之前，先对步骤a获取的原始雷达数据进行预处理，再在经过预处理的雷达图像中，提取双曲线反射信号上振幅绝对值最大的点，记录该采样点的坐标值。
所述的方法，其中：步骤c包括：将一倾斜埋于地下的柱状反射体置于空间直角坐标系中，根据雷达天线图像上的反射信号即电磁波的双程走时t与天线在测线上位置l的函数关系建立数学模型。
所述的方法，其中所述数学模型为：
其中，
所述步骤c包括，在建立了上述数学模型后，将步骤b中所选采样点的坐标值与上述方程进行曲线拟合；将
所述的方法，其中，所述方法还包括在所述步骤a中，在被探测柱状反射体的走向已知的情况下，在垂直于反射体的走向上布设多条测线的过程。
所述的方法，其中，所述方法还包括在所述步骤b中，在图像中寻找均匀分布于双曲线信号上待选点，从待选点所在的a扫描中找到振幅绝对值最大的点，并将待选点的位置调整至振幅绝对值最大点处，记录该采样点的坐标。
所述的方法，其中，所述方法所述步骤c中建立数学模型的过程具体为：
建立包含一倾斜埋于地下的柱状反射体的场景数学模型，在此场景中，柱状反射体的走向是已知的，天线测线与柱状反射体所在平面垂直，以天线测线所在的直线为x轴，以柱状反射体所在平面与地面的交线为y轴，z轴垂直于x轴与y轴构成空间右手直角坐标系。
所述的方法，其中，所述方法还包括在所述步骤c中推导柱状反射体双曲线信号所对应的函数关系式的过程，具体过程为：
首先，设柱状反射体的倾角为
其中h表示柱状反射体和z轴的交点与原点之间的距离；天线的测线与x轴重合，点p即反射体距离天线最近的点，就是反射体所在直线与测线所在直线的公垂线与反射体的交点，其坐标为：
其次，设天线从起始点出发所移动的距离为l，在实际情况下，雷达天线的起始位置与模型的起始位置之间存在位置差，此处引入偏移量l0，故天线的坐标为(l-l0,0,0)，天线与点p的距离d即可表示为：
然后，将此距离转化成电磁波的走时t：
其中，v为电磁波在介质中传播的速度，其与介质的介电常数之间的关系表示为：
c为真空中的光速，ε为介质的介电常数；
最后，根据上述关系式，得到雷达图像中双曲线信号对应的方程：
所述的方法，其中，所述方法还包括在所述步骤c中基于雷达图像中的双曲线信号进行曲线拟合的过程，即输入拟合方程以及步骤b中所选的多个采样点的坐标，得到包括反射体倾角
附图说明
以下附图仅旨在于对本发明做示意性说明和解释，并不限定本发明的范围。其中，
图1为根据本发明的一个具体实施例的利用探地雷达进行地下柱状反射体倾角估测的流程示意图；
图2为根据本发明提供的估测地下柱状反射体倾角大小的方法，在步骤a中利用雷达装置对地下反射体进行探测的示意图；
图3为根据本发明提供的估测地下柱状反射体倾角大小的方法，在步骤b中对雷达图像进行预处理的实例，(a)为探地雷达采集的原始图像，(b)为经过预处理后的雷达图像；
图4为根据本发明提供的估测地下柱状反射体倾角大小的方法构建的数学模型，描述了包含一倾斜埋于地下的柱状反射体的场景；
图5为根据本发明设计实验的实地场景。
具体实施方式
为了对本发明的技术特征、目的和效果有更加清楚的理解，现说明本发明的具体实施方式。但本领域的技术人员应该知道，以下实施例并不是对本发明技术方案作的唯一限定，凡是在本发明技术方案精神实质下所做的任何等同变换或改动，均应视为属于本发明的保护范围。
图1为利用探地雷达进行地下柱状反射体倾角估测的流程示意图；下面详细说明根据本发明提供的利用探地雷达估测地下柱状反射体倾角的原理，所述方法包括如下三个步骤：
步骤a，获取地下柱状物的雷达扫描图像；
步骤b，从步骤a所得到的雷达扫描图像中，在柱状物的双曲线反射信号上提取多个采样点的坐标；
步骤c，建立数学模型，利用相应算法进行曲线拟合，计算得到柱状物倾角的测量值。
对于上述三个步骤，步骤a为雷达数据收集和获取的过程，需在待测量地点完成；步骤b和c为雷达数据数据处理过程，可根据实时情况和需求进行安排完成。在步骤a中，本发明要求在柱状物上方布设多条平行测线，以保证实现对一段连续柱状物的倾角估测，便于分析柱状物倾角的变化情况。在步骤b中，本发明对双曲线信号上待选点所在的a扫描通道进行查看，并在每个a扫描上统一对最大振幅值出现的点进行选择，以保证曲线拟合的精度。在步骤c中，根据探地雷达工作原理，建立数学模型，并依据模型所反映的双曲线信号的函数关系式，进行曲线拟合，实现柱状反射体倾角参数的估测，这是本发明的一个重要优势。
下面分别对步骤a，b和c予以详细介绍：
步骤a，获取地下柱状物的雷达扫描图像
本发明使用收发天线一体的便携式雷达系统，对地下柱状反射体进行探测，测量过程如图2所示。在所要探测的柱状反射体走向已知的情况下，应在垂直于反射体的走向上，布设多条测线。在图2中，雷达天线沿其中一条测线(如图2中的l1)移动时，可得到带有地下反射体反射信号的一幅图像，根据图像中的反射信号，本发明所提供的方法将按照步骤b与步骤c，估测相应测线正下方位置处反射体的倾角(如图2中的
步骤b，从步骤a所得到的雷达扫描图像中，在柱状物的双曲线反射信号上提取多个采样点的坐标
本发明在提取采样点信号上点的坐标之前，首先需要对步骤a获取的原始雷达数据进行预处理，主要操作包括解振荡滤波、初至时间订正、背景去除和信号增益。图像预处理前后对比如图3(a)和3(b)所示，上述预处理操作能够达到突出双曲线信号、去除杂波的作用。
然后，为实现步骤c中的曲线拟合，需要从经过预处理的雷达图像中，提取双曲线反射信号上的点。由于反射体的双曲线反射信号是由正相与负相多层信号组成，在信号的不同层上取点，会影响步骤c中曲线拟合的结果。根据探地雷达的工作原理，电磁波走时与天线位置之间的函数关系对应的正是柱状物反射信号中，由各个a扫描内振幅最大点组成的双曲线信号条带。为了避免出现较大的拟合误差，本发明所提供的方法首先在图像上找到若干待选点作为采样点，这些待选点均匀分布于双曲线信号上，然后再对每个待选点所在的a扫描进行查看，并找到该a扫描上振幅绝对值最大的点，最后将待选点的位置调整至振幅绝对值最大点处。采样点的数量为10-20个，优选的是15个。
最后，记录所选点的坐标值。根据步骤a所获取的雷达数据，所述雷达数据上的每个采样点pi的坐标包括该点对应的天线在测线上的位置li和电磁波从天线到目标物再回到天线的双程走时ti。其中i是大于1的整数，表示采样点的序号。
步骤c，建立数学模型，利用相应算法进行曲线拟合，计算得到柱状物倾角的估测值
图4为本发明所提供的方法构建的简化数学模型，描述了包含一倾斜埋于地下的柱状反射体的场景。在此场景中，柱状反射体的走向是已知的，为保证雷达信号质量，天线测线(图4中的l轴：l1,l2或l3)应与柱状反射体所在平面垂直。以天线测线所在的直线为x轴，以柱状反射体所在平面与地面的交线为y轴，z轴垂直于x轴与y轴构成空间右手直角坐标系。
雷达天线在探测过程中，需沿测线匀速缓慢移动，并越过柱状反射体上方空间，根据空间解析几何知识，对于天线在测线上的每一个位置，总能在柱状物上找到一点p，使得点p到天线的距离比该柱状物上其它点到点p的距离都要近。考虑到探地雷达的成像原理，雷达图像上反射体的反射信号所表示的，即为电磁波从天线处发射并传播至点p再被反射回天线处所用的时间与天线在测线上位置的函数关系。
首先，设柱状反射体的倾角为
其中h表示柱状反射体和z轴的交点与原点之间的距离；
其次，设天线从起始点出发沿x轴所移动的距离为l，在实际情况下，雷达天线的起始位置与模型的起始位置之间存在位置差，此处引入偏移量l0，故天线的坐标为(l-l0,0,0)，天线与点p的距离d即可表示为：
然后，将此距离转化成电磁波的走时t：
其中，v为电磁波在介质中传播的速度，其与介质的介电常数之间的关系可表示为：
c为真空中的光速，ε为介质的介电常数。
最后，根据上述关系式，可得雷达图像中双曲线信号对应的方程：
根据数学模型得到上述方程后，本发明提供的方法将利用该方程作为拟合方程，结合步骤b中所选点的坐标值进行曲线拟合。曲线拟合需要借助计算机软件来完成，所采用的算法为levenberg-marquardt算法。拟合方程中，t与l分别为因变量与自变量，由于介质的介电常数ε需预先进行测量或估计，并作为先验知识代入方程，故待拟合的参数仅包括h、l0和反射体倾角
为了更好地说明本发明的技术效果，在实地设计控制实验，利用本发明提供的方法，用探地雷达对以一定倾角埋入沙土中的钢筋进行扫描，并从所得到图像的反射信号中提取钢筋倾角参数。控制实验选择在中国内蒙古自治区新巴尔虎左旗(48°15’16”n,118°25’30”e)的一处沙地。首先，选择一块地势平坦、无植被干扰的区域，并在该区域随机选择若干点，测量其土壤含水量，并取其平均值θ，根据topp公式：
θ＝-5.3×10-2+2.92×10-2ε-5.5×10-4ε2+4.3×10-6ε3
即可求解出实验区沙土的平均介电常数。之后在实验区挖7个长度为1.2m、宽度为0.6m的沙坑，用抹刀将沙坑底部修整为倾角为0°、10°、15°、20°、30°、40°、45°的斜坡，并保证底部斜面的中心点距离地面0.35m。挖好沙坑后，将长度为0.5m的钢筋置于底部斜面上，并保持与沙坑的长边垂直，如图5所示，再用原沙土将沙坑小心填满，将表面的沙土铺平。而后利用天线中心频率为1600mhz的雷达装置对钢筋进行扫描，扫描时测线应经过钢筋在地面上投影的中心点，并与钢筋在地面上的投影垂直，以保证获取雷达信号的质量。
在获得雷达图像，并对图像进行预处理后，即可利用步骤c对图像中的双曲线进行曲线拟合，在双曲线上确定待选点，记录点的坐标值，然后将测得的实验区沙土介电常数代入拟合公式中，并将坐标值与公式输入数学软件中，得到最终的拟合结果如下表所示，此处仅列出反射体倾角的拟合结果：
从结果可以看出，反射体倾角的拟合值与对应的真实值接近，估测的绝对误差不超过5°。本发明所提供的方法能够精准地估测地下柱状反射体的倾角，从反射体的反射信号出发，利用数学模型与曲线拟合算法，可以准确反演反射体的倾角，这将有助于快速、无损地获取地下目标物的走向情况，对地下钢筋、管道结构的探测与检修提供便利，对基础设施建设具有重要意义。
以上所述仅为本发明示意性的具体实施方式，并非用以限定本发明的范围。任何本领域的技术人员，在不脱离本发明的构思和原则的前提下所作的等同变化、修改与结合，均应属于本发明保护的范围。