基于探地雷达的地下异常体三维建模方法及装置

1.本发明涉及地下异常体检测领域，具体而言涉及一种基于探地雷达的地下异常体三维建模方法及装置。


背景技术：

2.随着城市地下空间开发利用的规模和范围的不断扩大，对于地下空间的勘探需求日益增加。在多种地球物理勘探方法中，探地雷达技术通过发射天线向地下发射电磁波，电磁波信号穿过地下介质，当遇到不同介电特性的物体时，如空洞、金属、水等，由于界面分界处介质的介电常数差异，会产生的幅值明显的回波信号，回波信号能够反应地下结构或者目标物的特征。探地雷达技术因定位准确、无损检测、探测快速、操作灵活、精度高等特点，被广泛应用于各领域的探测与成像。
3.目前，基于探地雷达的地下异常体三维建模方法通常需要先基于机器学习或者深度学习模型对异常区域进行自动识别，由于地下介质情况复杂，训练样本有限，训练的模型多适用于单一场景，且多基于回波的双曲线特征提取地下目标，对于其他任意形状的地下回波目标自动识别能力不足，普适性较差；另外，地下异常体通常通过三维切片展示，展示度不高，不能直观地进行观测，不利于地下灾害风险的实时监控和预警。


技术实现要素：

4.（一）发明目的本发明的目的是提供一种基于探地雷达的地下异常体三维建模方法及装置，用于解决现有的地下异常体检测展示度低，对非双曲线回波特征的异常体自动解译建模能力较差的问题。
5.（二）技术方案为解决上述问题，本发明的第一方面提供了一种基于探地雷达的地下异常体三维建模方法，包括：沿着多个测线获取待测区域的雷达测线数据；获取每个所述雷达测线数据对应的异常点云；获取每个异常点云对应的异常区的边缘轮廓；将全部所述边缘轮廓拼接，得到异常体的三维轮廓；获取落入所述三维轮廓内的插值点；为所述插值点赋予属性值，所述属性值根据所述插值点与所述三维轮廓的距离设置；基于所述插值点的位置与属性值构建三维模型。
6.在一些实施例中，获取每个异常点云对应的异常区的边缘轮廓包括：通过半径为设定值的球在异常点云上滚动；将所述球滚过的异常点连接得到初始边缘轮廓；
判断每个所述初始边缘轮廓对应的异常区是否为空洞，若是，则将所述初始边缘轮廓作为所述边缘轮廓；其中，所述设定值根据预估目标异常体的尺寸以及雷达中心频率进行预设。
7.在一些实施例中，判断每个所述初始边缘轮廓对应的异常区是否为空洞，包括：选取所述初始边缘轮廓对应的异常区中间道波形，在第一预设值至第二预设值区间内搜索振幅绝对值的最大值，对应的位置即为异常区的顶面子波主瓣位置，根据顶面子波主瓣位置获得顶面子波主瓣振幅值，其中，第一预设值为异常区上轮廓点，第二预设值为异常区上轮廓点+子波持时/采样间隔；在第二预设值至第三预设值区间内搜索振幅绝对值的最大值，对应的位置为异常区的底面子波主瓣位置，根据底面子波主瓣位置获得底面子波主瓣振幅值，其中，第三预设值为异常区下轮廓点位置；判断顶面子波主瓣振幅值与底面子波主瓣振幅值的乘积是否小于0，若是，则该异常区为空洞，若否，则该异常区不是空洞。
8.在一些实施例中，获取每个所述雷达测线数据对应的异常点云前，还包括：获取所述雷达测线数据能量横向叠加值最大的峰值采样点；基于所述峰值采样点生成均值作差区间；通过均值作差法对所述均值作差区间内的采样点进行直达波滤除。
9.在一些实施例中，所述雷达测线数据能量横向叠加值的计算公式如下：其中，x
i,j
为雷达测线数据的第j道回波数据的第i个采样点的雷达波振幅值，ei为雷达测线数据的n道回波数据在第i个采样点的雷达波能量横向叠加值。
10.在一些实施例中，所述直达波滤除的计算公式如下：其中，fc为探地雷达天线中心频率（mhz），t
ricker 为子波的周期（ns）,子波持时为2倍子波周期，dt为采样间隔，n
max
为雷达波能量横向叠加值峰值位置。
11.在一些实施例中，获取每个所述雷达测线数据对应的异常点云前，还包括：dewow滤波, 去除直流漂移；通过带通滤波消除所述雷达测线数据的噪音；获取每个所述雷达测线数据对应的异常点云包括：采用频率波数偏移方法根据指定波速对所述雷达测线数据进行偏移成像，压制绕射波的干扰；通过最大类间方差法对偏移成像后的所述雷达测线数据进行图像分割得到异常区域；根据所述指定波速对所述雷达测线数据进行时深转换，将异常区各点的双程回波走时转换为异常区各点所处深度，得到实际坐标尺寸的异常点云。
12.在一些实施例中，所述指定波速为预设波速区间内，对应的阈值分割区域面积最小的波速。
13.在一些实施例中，通过最大类间方差法对偏移成像后的所述雷达测线数据进行图像分割前，还包括：对偏移成像后的所述雷达测线数据进行深度补偿以增强所述雷达测线数据的信号强度。
14.在一些实施例中，所述属性值为插值点与所述三维轮廓距离的平方。
15.本发明第二方面提供了一种基于探地雷达的地下异常体三维建模装置，包括：第一获取模块，用于沿着多个测线获取待测区域的雷达测线数据；第二获取模块，用于获取每个所述雷达测线数据对应的异常点云；第三获取模块，用于获取每个异常点云对应的异常区的边缘轮廓；拼接模块，用于将全部所述边缘轮廓拼接，得到异常体的三维轮廓；第四获取模块，用于获取落入所述三维轮廓内的插值点；赋值模块，用于为所述插值点赋予属性值，所述属性值根据所述插值点与所述三维轮廓的距离设置；数据处理模块，用于基于所述插值点的位置与属性值构建三维模型。
16.（三）有益效果本发明提供了一种基于探地雷达的地下异常体三维建模方法，通过预先设定好参数，自动提取雷达测线数据的异常区域，将全部异常点云的边缘轮廓拼接能够得到异常体的三维轮廓，根据插值点与三维轮廓的距离对插值点赋予属性值并构建三维模型，能够实现雷达测线数据的自动实时处理、解释，能够对异常体是否为空洞有初步判断，异常体展现更加直观，有助于实现地下灾害的智能化管控。
附图说明
17.图1是本发明实施例提供的一种基于探地雷达的地下异常体三维建模方法的流程图；图2是本发明实施例提供的探地雷达扫描方式示意图；图3是本发明实施例提供的一种获取每个异常点云对应的异常区的边缘轮廓的流程图；图4是本发明实施例提供的点云滤波和异常体轮廓提取示意图；图5是本发明实施例提供的一种基于探地雷达的地下异常体三维建模方法的详细流程图；图6是本发明实施例提供的偏移成像与图像分割示意图；图7是本发明实施例提供的异常区属性判别示意图；图8是本发明实施例提供的异常体的三维轮廓的拼接图以及建模效果图；图9是本发明实施例提供的一种基于探地雷达的地下异常体三维建模装置的结构示意图。
具体实施方式
18.在下文的描述中，给出了大量具体的细节以便提供对本发明更为彻底的理解。然而，对于本领域技术人员而言显而易见的是，本发明可以无需一个或多个这些细节而得以实施。在其他的例子中，为了避免与本发明发生混淆，对于本领域公知的一些技术特征未进行描述。
19.应予以注意的是，这里所使用的术语仅是为了描述具体实施例，而非意图限制根据本发明的示例性实施例。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式。此外，还应当理解的是，当在本说明书中使用术语“包含”和/或“包括”时，其指明存在特征、整体、步骤、操作、元件和/或组件，但不排除存在或附加一个或多个其他特征、整体、步骤、操作、元件、组件和/或它们的组合。
20.现在，将参照附图更详细地描述根据本发明的示例性实施例。然而，这些示例性实施例可以多种不同的形式来实施，并且不应当被解释为只限于这里所阐述的实施例。应当理解的是，提供这些实施例是为了使得本发明的公开彻底且完整，并且将这些示例性实施例的构思充分传达给本领域普通技术人员。
21.图1是本发明实施例提供的一种基于探地雷达的地下异常体三维建模方法的流程图。如图1所示，本发明实施例提供了一种基于探地雷达的地下异常体三维建模方法包括如下步骤：s101：沿着多个测线获取待测区域的雷达测线数据。如图2所示，多个测线可平行间隔设置，探地雷达以测量轮触发的连续测量模式沿着测线水平移动，接收的回波数据为雷达测线数据。
22.s102：获取每个雷达测线数据对应的异常点云。探地雷达通过发射天线向地下发射电磁波，电磁波信号穿过地下介质。当遇到不同介电特性的物体时，如空洞、金属、水等，由于界面分界处介质的介电常数差异，会产生幅值明显的回波信号，回波信号经过处理能够得到异常点云。异常点云用于表征异常体的剖面特征。
23.s103：获取每个异常点云对应的异常区的边缘轮廓。本实施例中，边缘轮廓可以为一个或多个，边缘轮廓用于表征异常体剖面的边界线。
24.s104：将全部边缘轮廓拼接，得到异常体的三维轮廓。本实施例中，三维轮廓为一个或多个，三维轮廓用于表征异常体的轮廓特征。
25.s105：获取落入三维轮廓内的插值点。
26.s106：为插值点赋予属性值，属性值根据插值点与三维轮廓的距离设置。本实施例中，插值点与三维轮廓的距离越近，属性值的绝对值越小。
27.s107：基于插值点的位置与属性值构建三维模型。
28.本发明实施例提供的一种基于探地雷达的地下异常体三维建模方法，通过预先设定好参数，自动提取雷达测线数据的异常区域，能够逐个将单测线的所有异常区域逐一进行属性判别，筛去非空洞区域，将多测线全部异常点云的边缘轮廓拼接能够得到地下异常体的三维轮廓，根据插值点与三维轮廓的距离对插值点赋予属性值并构建三维模型，能够实现雷达测线数据的自动实时处理、解释，能够对异常体是否为空洞有初步判断，异常体展现更加直观，有助于实现地下灾害的智能化管控。
29.图3是本发明实施例提供的一种获取每个异常点云对应的异常区的边缘轮廓的流
程图。如图3所示，获取每个异常点云对应的异常区的边缘轮廓包括如下步骤：s1031：通过半径为设定值的球在异常点云上滚动；s1032：将球滚过的异常点连接得到初始边缘轮廓；s1033：选取初始边缘轮廓对应的异常区中间道波形，在第一预设值至第二预设值区间内搜索振幅绝对值的最大值，对应的位置即为异常区的顶面子波主瓣位置，根据顶面子波主瓣位置获得顶面子波主瓣振幅值，其中，第一预设值为异常区上轮廓点，第二预设值为异常区上轮廓点+子波持时/采样间隔；s1034：在第二预设值至第三预设值区间内搜索振幅绝对值的最大值，对应的位置为异常区的底面子波主瓣位置，根据底面子波主瓣位置获得底面子波主瓣振幅值，其中，第三预设值为异常区下轮廓点位置；s1035：判断顶面子波主瓣振幅值与底面子波主瓣振幅值的乘积是否小于0，若是，则该异常区为空洞，将初始边缘轮廓作为边缘轮廓，若否，则该异常区不是空洞，删除该异常区对应的初始边缘轮廓；其中，设定值根据预估目标异常体的尺寸以及雷达中心频率进行预设。
30.电磁波从介电系数大的介质入射介电系数小的介质,回波相位不变,从介电系数小的介质到介电系数大的介质反向，空洞底层反射波必然与空洞顶端反射波的相位相反，而通过希尔伯特变换获得的瞬时相位变化信息不易识别。本实施例中，以特征明显易识别的雷达子波主瓣的正负作为相位判别的依据，以异常区子波顶底界面主瓣正负不一致初判为空洞，将主瓣相位一致的不是空洞特征的异常区从整体异常点云中剔除。
31.具体地，获取每个异常点云的边缘轮廓包括：通过半径为设定值的球在异常点云上滚动；将球滚过的异常点连接得到边缘轮廓；其中，可根据预估目标异常体的尺寸以及雷达中心频率进行预设。基于剖面数据自动生成异常体外轮廓存在以下难点：实测数据中单个剖面数据常常存在多个异常体区域，直接对单剖面多个异常体求取外轮廓可能使多个异常体连接一个封闭整体，不符合地下异常体实际的分布情况，同时，不关心的小尺寸异常体也会对目标异常体的轮廓生成带来干扰。为克服以上难点，本实施例，利用matlab中alphashape函数的特点，通过将alpha值设为2，自动提取异常体外轮廓，本实施例还包括对经过上述步骤得到的外轮廓shp结构体利用shp.regionthreshold删除小于设定剖面区域阈值的边缘轮廓，排除小尺寸区域干扰，剖面区域阈值可以为0.2平方米，最终利用boundaryfacets函数获取异常点云的边缘轮廓坐标。图4为经过偏移处理后，利用alphashape函数求取异常区域外轮廓的效果。
32.图5是本发明另一实施例提供的一种基于探地雷达的地下异常体三维建模方法的流程图。如图5所示，本发明另一实施例提供的一种基于探地雷达的地下异常体三维建模方法包括数据采集，数据自动化处理和三维建模与可视化表达。
33.本实施例中，数据采集包括待测区域规划，测线布置和雷达参数选择。数据自动化处理包括单测线数据处理和多测线数据处理。
34.本实施例中，单测线数据处理包括：s501：直达波滤除。当地下深层反射信号较为强烈时，对整体雷达剖面图像进行均值作差会导致反射信号所处深度出现横状条纹，对异常体识别带来不必要的干扰。本技术采用能量分布法确定直达波位置，从而根据雷达中心频率限定均值作差区间。由于直达波
到达时能量会产生明显突变，且能量幅值最大，通过计算雷达图像的横向能量叠加分布图，将最大横向能量位置定位为直达波峰值时刻。水平方向的能量分布及直达波峰值时刻计算如下式所示：其中，x
i,j
为第j道雷达测线数据的第i个采样点的雷达波振幅值，ei为多道雷达测线数据在第i个采样点的雷达波能量叠加值。通过读取雷达数据头文件中的中心频率，计算雷达子波持时，设定记录开始点至（峰值时刻+2倍子波周期）作为均值作差区间，其中雷达子波周期计算如下：fc为探地雷达天线中心频率（mhz），t
ricker 为子波的周期（ns）。则改进后的均值法去直达波公式为：其中，fc为探地雷达天线中心频率（mhz），子波持时为2倍子波周期，dt为采样间隔，n
max
为雷达波能量横向叠加值峰值位置。
35.本技术优势在于通过分析雷达图像特征，自动找出直达波位置，限制均值作差的区间，从而最大程度减小均值作差对雷达信号的影响，保证地下深层异常体的准确识别。
36.s502：滤波。雷达信号经常会产生一些低频和高频噪音，干扰目标信号的识别。本技术通过读取雷达数据头文件的天线中心频率fc，选取0.25倍fc至2.5倍fc作为带通滤波器滤波区间，对测线的每道数据进行带通滤波，消除雷达信号中的低频和高频噪音。
37.s503：波速搜索。在步骤s504中，当设定波速过大，会发现数据偏移过度，雷达图像上会呈现倒置双曲线，若波速过小，则偏移不够，恰好的波速会将双曲线收缩到点状或圈状。基于这一特点，本技术将偏移成像和阈值分割加入到波速搜索当中，通过设定电磁波可能波速变化范围，依据偏移成像后阈值分割得出的异常体区域面积最小，确定出最佳波速。探地雷达波速标定通常需要获取已知深度的地下目标物回波数据，进而反算出探测区域的波速，本方法的优势在于省去测量人员实地标定波速的环节，直接根据采集得到的数据遍历计算即可得到所测区域波速，方法快捷简便，在仅有雷达数据的情况下即可快速获得波速。
38.s504：偏移成像；当天线由远及近扫描有限尺寸的物体时，该物体在雷达图像上表现为双曲线形态。绕射会模糊一些有用的信息，导致对地下目标体的尺寸和形状作出错误的解译。本实施例采用基于波动方程的偏移方法—频率波数偏移方法，又称为f-k方法，该方法精度高、稳定性好、运算速度快，适应于波传播速度变化较小的情形，通过偏移将绕射波回归到真实位置。
39.s505：自动增益；本实施例通过采用自动增益方法，对偏移后的雷达图像进行深度补偿，增强深层信号的强度。
40.s506：阈值分割。为了进一步提取异常体目标，本实施例采用最大类间方差法（otsu法）自动对图像进行分割。通过求取经偏移、自动增益后的雷达剖面数据的绝对值，将
雷达数据转换至0~255范围内，得到剖面数据的灰度图。采用otsu法对灰度图进行阈值计算，基于阈值对偏移处理后的图像进行阈值分割，提取得到异常体区域。图6为通过gprmax得到的bscan剖面图以及经过偏移成像、自适应阈值分割得到的结果图。
41.s507：判断异常区域是否最小，若否则返回s503。
42.s508：波速锁定，剖面区域提取。s507的判断结果为是，则根据确定的波速，获取异常区域。
43.s509：对异常区域的雷达测线数据进行时深转换，将双程走时转换为深度，提取经偏移成像、阈值分割后的剖面异常点云，数据第一列为沿着测线方向长度，第二列为所处深度得到异常点云。
44.s510：对该测线异常区的属性逐个进行判断。自动选取异常区中间道波形，在第一预设值至第二预设值区间内搜索振幅绝对值的最大值，对应的位置即为异常区的顶面子波主瓣位置，根据顶面子波主瓣位置获得顶面子波主瓣振幅值，其中，第一预设值为异常区上轮廓点，第二预设值为异常区上轮廓点+子波持时/采样间隔；在第二预设值至第三预设值区间内搜索振幅绝对值的最大值，对应的位置为异常区的底面子波主瓣位置，根据底面子波主瓣位置获得底面子波主瓣振幅值，其中，第三预设值为异常区下轮廓点位置；根据顶面子波主瓣振幅值与底面子波主瓣振幅值乘积是否小于0，判断该异常区是否为空洞，若大于0，排除此异常区为空洞的可能。图7为异常区中间道波形拾取以及中间道波形顶底面子波波形主瓣位置。
45.s511：多测线数据处理包括：根据多测线经纬度，将经纬度进行坐标系转换，对经过上述单测线处理的异常点云进行循环储存，得到包含水平坐标（x,y）和深度坐标(z)的多测线异常点云数据；s512：可视化表达包括设定插值控制点、赋值构造三维数组、isosurface提取点、面信息和三维obj模型生成。
46.相关技术中，三维模型构建通过先将剖面轮廓线连接成三角面片，接着将多个面片拼接成表面模型完成建模，但是多测线间剖面求取外轮廓时，两个剖面异常体之间的连接、包含、分离等关系常常需要人工结合常识和规律进行判断，当剖面轮廓线复杂，难以取得较好效果。为加快建模速度，本实施例采用采用设置控制点的方法构造三维立方体数组，利用等值面法求取三维模型点、面信息。利用numregions(shp)求取每个测线剖面异常区域个数，对异常区域外轮廓逐一按照如下方式进行赋值：首先，按照100x100的数组对剖面进行插值，得到控制点坐标，接着逐个判断控制点是否在异常体区域内部，若在区域内，计算该控制点距该异常区域轮廓点最近距离d
in
，将该距离的平方（d
in2
）作为该点的属性值；若在区域外，计算该控制点距所有异常区域轮廓最近的距离d
out
，将该距离的平方的负值（-d
out2
）作为该点的属性值；若在异常区域边界上，该点的属性值设为0。将三维轮廓上的异常点增设为控制点，并将属性值赋值为0；寻找剖面所有异常区域内的任意一点作为该异常区域内部控制点，计算该增补控制点距该异常区域的三维轮廓最近的距离，将该距离的平方作为该点的属性值。对合并后的控制点的属性值进行插值。利用scatteredinterpolant函数，设定测线最大长度、最大深度、多测线距离作为插值区间，综合考虑模型精度需求与建模效率，合理设置插值步长，对多测线剖面控制点的属性值进行进一步三维插值，得到插值
后三维区域属性值。采用matlab中的isosurface函数，对构造得到三维数组v提取参数值为0的等值面fv，利用isonormals函数计算等值面顶点的法向量，得到写入obj所需的点信息fv.vertices以及面信息fv.faces。最后利用write_wobj函数生成异常体三维模型。图8是发明实施例提供的异常体三维模型建模效果图，左侧为异常体的点云模型，右侧为采用控制点构造三维数组生成的三维模型，如图8所示，本实施例提供的异常体三维建模方法能够准确地建立地下空间三维模型自动载入地理信息系统。
47.本实施例提供的基于探地雷达的地下异常体三维建模方法能够代替人工的推测与分析，结合剖面测线数据，实现自动判断地下结构的连接与否以及异常体的合理分割，受剖面轮廓线影响较小，自动生成的模型可以直接导入地理信息系统进行展示，实现地下三维场景的一体化构建与动态更新，对地下探测中的异常体展现更加直观。
48.图9是本发明实施例提供的一种基于探地雷达的地下异常体三维建模装置的结构示意图。如图9所示，基于同一发明构思，本发明实施例提供了一种基于探地雷达的地下异常体三维建模装置，包括：第一获取模块11，用于沿着多个测线获取待测区域的雷达测线数据；第二获取模块12，用于获取每个雷达测线数据对应的异常点云；第三获取模块13，用于获取每个异常区的边缘轮廓；属性判别模块14，用于识别每个异常区的属性，筛去非空洞异常区域；拼接模块15，用于将多测线全部边轮廓拼接，得到异常体的三维轮廓；第四获取模块16，用于判断获取落入三维轮廓内的插值点；赋值模块17，用于为插值点赋予属性值，属性值根据插值点与三维轮廓的距离设置；数据处理模块18，用于基于插值的位置与属性值构建三维模型。
49.应当理解的是，本发明的上述具体实施方式仅仅用于示例性说明或解释本发明的原理，而不构成对本发明的限制。因此，在不偏离本发明的精神和范围的情况下所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。此外，本发明所附权利要求旨在涵盖落入所附权利要求范围和边界、或者这种范围和边界的等同形式内的全部变化和修改例。