CN107544095B - 一种地面三维激光点云与探地雷达图像融合的方法 - Google Patents

Abstract

本发明公开了一种地面三维激光点云与探地雷达图像融合的方法，在探地雷达图像与GPS数据同步采集基础上，利用改进的电磁波路径传播模型实现探地雷达图像的时深转换，以地面三维激光数据所在坐标系为参考，结合地理参考绝对定位方程将两不同源异构数据转换到同一坐标系下，实现数据融合。本发明根据地面三维激光扫描仪与探地雷达的数据特点，提出基于三维激光点云与探地雷达图像融合方法，实验结果和误差分析证明，流程比较简单，融合精度较高，具有广泛的适应性；克服两不同源异构数据的分辨率和数据量之间的差异，充分发挥三维激光扫描仪与探地雷达两个传感器各自优势，实现两种不同类型数据之间的有效融合。

Description

一种地面三维激光点云与探地雷达图像融合的方法

技术领域

本发明涉及图像处理领域，具体涉及一种地面三维激光点云与探地雷达图像融合的方法。

背景技术

探地雷达作为利用高频电磁波来确定浅表地层构造的无损地球物理探测技术，由于其操作简便、检测范围广、探测深度深和分辨率高等特点，已被广泛应用于考古探测、地质工程勘察、公路病害检测、道路空洞和地下水污染探测等诸多领域，取得了较好的应用效果。地面三维激光扫描仪以非接触的测量方式快速获取被测物体外部几何形态或三维场景，近年来应用也很广泛。

地面三维激光扫描仪可以获取被测物体外部精细的几何形态或三维场景，但无法穿透物体实现对其内部结构的有效探测，而探地雷达可以穿透物体或地面获取内部异常形态。因此，将地面三维激光扫描仪和探地雷达获取的数据进行融合成为目前研究的热点之一。

探地雷达获取的数据是由连续单道波组成在水平和竖直方向上呈一定比例的二维时间剖面，而地面三维激光扫描仪获取的则是具有拓扑关系的离散点云。地面三维激光点云与探地雷达图像为非同源遥感数据，两种数据源之间既没有统一的坐标系统，又不存在同名特征，目前没有一种方案能克服两种非同源异构特征而实现两数据之间的配准融合。

发明内容

针对地面激光点云与探地雷达图像的数据没有较好的融合实现方法，融合精度不确定的技术问题，本发明提出一种地面三维激光点云与探地雷达图像融合的方法，解决了两异构数据源之间的配准整合问题，实现了两种数据之间的融合显示，很好地发挥了两种技术的各自优势，克服了单一数据源缺陷，促进了两种数据之间的相互解译，提高了探地雷达图像解译的准确性。

为了解决上述技术问题，本发明的技术方案是：一种地面三维激光点云与探地雷达图像融合的方法，在探地雷达图像与GPS数据同步采集基础上，利用改进的电磁波路径传播模型实现探地雷达图像的时深转换，以地面三维激光数据所在坐标系为参考，结合地理参考绝对定位方程将两不同源异构数据转换到同一坐标系下，实现数据融合，其步骤如下：

步骤一：地面三维激光扫描仪采集的点云数据和高分辨率影像直接融合生成彩色激光点云，进行滤波和抽稀处理后，结合三个或三个以上地面已知控制点通过七参数转换法将彩色激光点云由仪器坐标系转换到大地坐标系下；

步骤二：对探地雷达采集的二维时间剖面图像进行相应的数据后处理，后处理包括解震荡滤波、自动增益、背景滤波、带通滤波和图像平滑处理；

步骤三：根据探地雷达和GPS同步采集后的时间同步文件，将处理后的探地雷达图像与GPS坐标进行精确匹配，得到探地雷达图像竖直方向上首个采样点的坐标；

步骤四：在改进的电磁波传播路径模型基础上，根据电磁波在介质中的传播速度实现探地雷达图像纵轴方向上的时深转换，由首个采样点的坐标得到探地雷达图像上所有采样点坐标；

步骤五：根据探地雷达数据采集方式建立了探地雷达探测空间坐标系，以三维激光点云所在坐标系为参考，建立三维激光点云与探地雷达坐标系之间的参数转换模型，将探地雷达图像上所有采样点统一到同一坐标系中，实现探地雷达图像与三维激光点云的融合。

所述七参数转换法将彩色激光点云转换为大地坐标系下的方法为：

其中，(x_j，y_j，z_j)为激光点云j在激光仪器坐标系下的坐标，(X_D，Y_D，Z_D)为激光点云j在大地坐标系下的坐标；ω，γ分别为x,y,z轴对应的旋转参数，Δx，Δy，Δz分别为x,y,z轴的平移参数，m为尺度变换参数；当地面已知控制点多于3个时，通过最小二乘法求7个转换参数的最或然值。

所述探地雷达与GPS同步采集的方法为：

(1).探地雷达采集的二维时间剖面图像e(x_i,t_j)，1≤i≤M,1≤j≤N，M为探地雷达图像道数，i为探地雷达图像的第i道数，N为每道数据上的采样点数，j为第j个采样点，则探地雷达在水平距离x_i＝i·Δx，Δx为采样的道间距，探地雷达在纵轴上的时间往返信号为t_j＝j·Δt，Δt为采样时间间隔；

(2).测距轮的精度Δd＝C/N_d，其中C为测距轮的周长，N_d为测距轮旋转一周的脉冲个数；

(3).由于探地雷达和GPS之间的数据采集是通过测距轮同步触发，探地雷达主机与GPS主机接收的脉冲数应一致，即N_GPS＝M·(Δx/Δd)，N_GPS为GPS打标文件中记录的脉冲事件个数；

(4).GPS接收机同时获取每一个外部脉冲和绝对时刻的空间位置坐标，建立采集的探地雷达数据与GPS数据的对应关系为：

(x_i,y_i,z_i)_GPS＝i·(Δx/Δd)·(x_i,y_i,z_i)_GPR，1≤i≤M,

其中，(x_i,y_i,z_i)_GPS和(x_i,y_i,z_i)_GPR分别为同一时刻探地雷达主机第i道数据和GPS接收机获取的位置信息；从而获取探地雷达图像上每道数据的位置信息；

(5).由于GPS天线的中心与探地雷达天线的中心重合，这坐标转换的过程中x和y轴方向的平移矢量为零，只需计算z轴方向的平移矢量，那么探地雷达图像第i道数据的位置信息为：

(x_s,y_s,z_s)_GPR＝(x_i,y_i,(z_i-h_GPS))_GPS

其中，h_GPS为GPS天线到探地雷达天线中心位置的高度；(x_s,y_s,z_s)_GPR探地雷达天线中心位置的坐标，(x_i,y_i,z_i)_GPS分别为探地雷达天线上的流动站GPS的坐标。

所述探地雷达图像与GPS坐标进行精确匹配的方法是：采用探地雷达与GPS硬件集成的方式来实现探地雷达图像空间位置的获取，将GPS天线固定在探地雷达天线的中心位置进行数据的采集；当探地雷达的主机开始工作时，根据设置的采集间隔，高精度测距轮在行进的过程中会不断触发探地雷达主机采集数据，与此同时也会触发GPS天线主机的I/O口以打标文件的方式记录下此时刻的GPS时间，并存储在记录卡内；数据采集后，根据基站和流动站GPS的差分处理结果，结合道间距和测距轮的精度通过插值实现探地雷达图像上道数据与GPS时间的精确匹配，使每道探地雷达数据都具有精确的位置信息。

所述改进的电磁波传播路径模型包括探地雷达发射天线与接收天线间的电磁波在空气中传播的直达波和地面下的反射波，电磁波的地面下反射面的深度h为：

t₁＝t-t₀

式中，v为电磁波在介质中的传播速度；d为发射天线和接收天线之间的距离，一般地，探地雷达屏蔽天线收发天线之间的距离是固定不变的，不同频率天线的天线距也不相同；t为探地雷达剖面记录下的电磁波传播的双程时间；t₀为地表直达波传播的双程时间；t₁为电磁波从发射天线到遇到异常体的双程传播时间。

将探地雷达图像的数据从探测空间坐标系到与地面三维激光相同的大地参考坐标系的转换的方法是：首先，根据探地雷达图像数据的特点，建立探测空间坐标系o-xyz，以测线的起点o作为坐标原点，x轴表示探地雷达天线沿测线采集数据的方向，y轴表示探地雷达不同测线的方向，z轴表示探地雷达天线探测的方向；

其次，以地面三维激光点云所在大地空间坐标系O-XYZ为参考，将探地雷达探测坐标系o-xyz中任一采样点p(x,y,z)到坐标系O-XYZ下采样点P(X,Y,Z)的转换关系表示为：

即

式中，(X，Y，Z)为探地雷达采样点P在大地空间坐标系O-XYZ下的坐标，(x，y，z)为探地雷达在探测空间坐标系o-xyz下的采样点p的位置坐标，(x₀,y₀,z₀)为探测空间坐标系原点o在大地空间参考坐标系O-XYZ中的坐标，矩阵为两坐标系的空间位置旋转矩阵，九个系数分别为两轴系间的方向余弦。

对探地雷达图像的数据进行属性归一化处理的方法是：激光点云的数据为(X,Y,Z,R,G,B)的形式，其中，X,Y,Z代表激光点云在大地空间参考坐标系下的坐标，R，G，B值代表其纹理颜色，取值范围为0-255之间；转化后探地雷达图像的数据格式为(X,Y,Z,Q)，其中，X,Y,Z表示探地雷达采样点与激光点云对应大地坐标系下的坐标，Q代表采样点的电磁波瞬时振幅值；对探地雷达数据的属性值进行归一化处理，使其与激光点云的颜色属性值处在统一量级上：

其中，Q_max代表探地雷达图像上所有采样点振幅值中的最大值，Q_min表示探地雷达图像上所有采样点振幅值中的最小值，Q表示探地雷达图像上任一采样点的振幅值，Q₀为强度归一化后探地雷达图像上任一采样点的振幅值。

利用空间线性的插值算法实现探地雷达数据的三维显示的方法为：P(x_p,y_p,z_p)为两已知剖面之间对应未知一点的空间坐标为：

x_p＝x_a+K(x_a-x_b)

y_p＝y_a+K(y_a-y_b)

z_p＝z_a+K(z_a-z_b)

其中，P_a(x_a,y_a,z_a)和P_b(x_b,y_b,z_b)分别为两相邻探地雷达二维时间剖面；在求出未知点的空间未知后，该点的强度值Q_p为：

其中，Q_a和Q_b为两相邻探地雷达二维时间剖面上的相对应已知点的强度值；K值则是通过数据采集时采样道间距Δx和二维剖面之间的距离D确定的：(D/Δx)≤K≤D，K为整数。

本发明根据激光点云与探地雷达数据特点，提出基于三维激光点云与探地雷达图像融合方法，实验结果和误差分析证明，流程比较简单，融合精度较高，具有广泛的适应性。本发明地面激光点云与探地雷达图像的数据融合发挥了不同数据源的优势，克服单一数据源的缺陷，促进了两种数据之间的相互解译，具有广阔的应用前景。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明的流程图。

图2为探地雷达与GPS数据时间同步的流程图。

图3为本发明探地雷达电磁波传播路径模型。

图4为本发明探地雷达图像中采样点的坐标。

图5为本发明试验区的图片。

图6为本发明激光点云与探地雷达图像的数据融合显示效果图，其中(a)三维激光点云与探地雷达二维剖面显示图，(b)三维激光点云与探地雷达数据三维显示图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有付出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

如图1所示，一种地面三维激光点云与探地雷达图像融合的方法，在探地雷达与GPS数据同步采集数据的基础上改进电磁波路径传播模型，结合地理参考绝对定位方程和坐标转换将两不同源异构数据统一到相同坐标系下，以实现其数据融合显示，其流程如图1所示。

步骤一：地面三维激光扫描仪采集的点云数据和高分辨率影像直接融合生成彩色激光点云，进行滤波和抽稀处理后，结合三个或三个以上地面已知控制点通过七参数转换法将彩色激光点云由仪器坐标系转换到大地坐标系下。

地面三维激光扫描过程中，不仅可以获取空间复杂目标体表面几何图像数据，也可以通过内置或外置高高分辨率相机获取空间的纹理信息。数据采集之前，通过对相机进行检校和标定以获取其外方位元素，实现地面三维激光扫描仪点云与影像的配准，使点云获取彩色纹理信息，生成三维真彩色点云。

地面三维激光扫描仪获取的数据是仪器坐标系下的点云，可通过七参数转换法将彩色激光点云转换为大地坐标系下，其过程如下：

其中，(x_j，y_j，z_j)为激光点云j在激光仪器坐标系下的坐标，(X_D，Y_D，Z_D)为激光点云j在大地坐标系下的坐标；ω，γ分别为x,y,z轴对应的旋转参数，Δx，Δy，Δz为平移参数，m为尺度变换参数。要求得这7个转换参数，至少需要3个控制点，当控制点多于3个时，可通过最小二乘法求7个转换参数的最或然值。

由于设备系统内部因素和外界环境的影响，地面激光扫描仪在数据采集过程中不可避免地会产生一些噪声点，其中主要是由外界环境因素引起的偶然误差，例如由于车辆、行人、建筑物和树木等因素对扫描目标的遮挡而形成散乱的点云或空洞，通过滤波处理可以剔除点云数据中的噪声点。

地面激光扫描仪扫描速度较快，数据密度大，经过多站拼接后，会出现大量的重叠区域，并且点云是以散乱的形式存在的，占用较大的数据空间，这样不利于后期数据处理、传输及应用。需要在保持一定精度范围内，对点云进行抽稀，以减少冗余数据。

步骤二：对探地雷达采集的探地雷达图像进行相应的数据后处理，后处理包括解震荡滤波、去除地面波、自动增益、背景滤波和带通滤波处理；

探地雷达图像的后处理是对数据进行滤波处理，过程依次包括解震荡滤波、自动增益(AGC)、背景滤波、带通滤波和图像平滑处理。解震荡滤波主要是去除电磁波信号中的直流成分或直流偏移，自动增益的作用是增强电磁波后视信号振幅，避免有效信号的湮没；背景滤波的主要作用是去除水平信号影响，特别是天线的振铃信号；带通滤波作用保持有效信号的截止频率，去除环境或系统噪声；图像平滑处理主要从整体上压制信号散射，提高信噪比，增强可视效果。

步骤三：根据探地雷达和GPS同步采集后的时间同步文件，将后处理的探地雷达图像与GPS坐标进行精确匹配，得到探地雷达图像竖直方向上首个采样点的坐标。

探地雷达与GPS数据之间的时间同步是影响激光点云与探地雷达图像融合精度的主要因素之一。探地雷达与GPS之间的数据一般是通过串口的方式实时通信，这种方式虽然实现过程比较简单，但数据更新速度较慢，数据采集过程中容易发生位置信息的丢失，无法满足探地雷达天线大面积、长距离采集数据的要求。本发明采用探地雷达与GPS硬件集成的方式来实现探地雷达图像空间位置的获取，使探地雷达图像与GPS数据精确匹配。为实现动态实时记录探地雷达天线的位置信息，将GPS天线固定在探地雷达天线的中心位置以组成整体的采集系统。当探地雷达开始工作时，根据设置的采集间隔，高精度测距轮在行进的过程中会不断触发探地雷达主机采集数据，与此同时也会触发GPS天线主机的I/O口以打标文件的方式记录下此时刻的GPS时间，并存储在记录卡内。数据采集后，根据基站和流动站GPS的差分处理结果，结合道间距和测距轮的精度通过插值实现探地雷达图像上道数据与GPS时间的精确匹配，使每道探地雷达数据都具有精确的位置信息，探地雷达与GPS之间的数据同步采集流程，如图2所示。

探地雷达与GPS时间同步的方法如下：

1.探地雷达采集的二维时间剖面图像e(x_i,t_j)，1≤i≤M,1≤j≤N，M为探地雷达图像道数，i为第i道探地雷达图像，N为每道数据上的采样点数，j为第j个采样点，则探地雷达在水平距离x_i＝i·Δx，Δx为采样的道间距，探地雷达在纵轴上的时间往返信号为t_j＝j·Δt，Δt为采样时间间隔。

2.测距轮的精度Δd＝C/N_d，其中C为测距轮的周长，N_d为测距轮旋转一周的脉冲个数。

3.由于探地雷达和GPS之间的数据采集是通过测距轮同步触发，那么探地雷达主机与GPS主机接收的脉冲数应一致，即N_GPS＝M·(Δx/Δd)，N_GPS为GPS打标文件中记录的脉冲事件个数。

4.GPS接收机获取每一个外部脉冲的同时也获取这一绝对时刻的空间位置坐标，在步骤3的基础上建立采集的探地雷达数据与GPS数据的对应关系：

(x_i,y_i,z_i)_GPS＝i·(Δx/Δd)·(x_i,y_i,z_i)_GPR，1≤i≤M,

其中，(x_i,y_i,z_i)_GPS和(x_i,y_i,z_i)_GPR为同一时刻探地雷达主机第i道数据和GPS接收机获取的位置信息。这样就可以获取探地雷达每道数据的位置信息。

5.(x_i,y_i,z_i)_GPR获取的是固定在探地雷达天线正上方的GPS天线位置处的位置信息。由于GPS天线中心与探地雷达天线的中心重合，这样坐标转换的过程中x和y轴方向的平移矢量为零，只需计算z轴方向的平移矢量即可，那么探地雷达图像第i道数据的位置信息为

(x_s,y_s,z_s)_GPR＝(x_i,y_i,(z_i-h_GPS))_GPS

其中，h_GPS为GPS天线到探地雷达天线中心位置的高度，一般为固定距离。

探地雷达与GPS数据之间的同步采集，不仅使采集到的探地雷达图像具有精确的位置信息，而且利用GPS记录下沿测线高精度的高程数据可用于探地雷达图像的地形校正，很大程度上消除了地形因素对探地雷达图像解译的影响。

步骤四：在改进的电磁波传播路径模型基础上，结合电磁波在介质中的传播速度实现探地雷达图像纵轴方向上的时深转换，得到探地雷达图像上所有采样点坐标。

探地雷达图像是以采集起点作为坐标原点的二维时间剖面，横轴表示探地雷达天线在地面上行进的距离，纵轴表示电磁波在介质中传播遇到异常体并反射回来的双程时间。为了实现激光点云与探地雷达图像的精确融合，必须对探地雷达图像竖直方向上进行时深变换，这就涉及到电磁波传播路径模型的选取。实际工程探测中，若电磁波在介质中的传播速度已知，那么探地雷达图像纵轴方向上的深度为：

其中，h为目标体的真实深度，v为电磁波在介质中的传播速度，t为电磁波到达目标体并反射回来的双程时间。由式1可知，此电磁波传播路径模型比较简单，仅反映电磁波在介质中的主要传播路径，忽略了发射天线和接收天线之间的距离和直达波因素对地下目标体精确定位的影响，适用于精度要求不高或收发天线距离较小的情况。

电磁波传播路径模型进行探地雷达图像纵轴方向上的时深转换，必须考虑到探地雷达发射天线和接收天线之间的距离和直达波等因素的影响。综合考虑探地雷达天线采集数据的原理和主要影响因素，即考虑探地雷达的发射天线和接收天线之间的距离和直达因素的影响，改进的电磁波传播路径模型如图3所示。由于探地雷达的收发天线具有极化特点，电磁波在传播过程中的能量散射可以忽略不计，因此电磁波的传播路径主要分为空气中传播的直达波和地面下的反射波，即图3中电磁波的路径1和路径2所示。

那么，通过时深变换目标体的真实深度及地面下反射面的深度h可通过公式计算，即

t₁＝t_j-t₀

式中，v为电磁波在介质中的传播速度；d为发射天线和接收天线之间的距离，一般地，探地雷达屏蔽天线收发天线之间的距离是固定不变的，不同频率天线的天线距也不相同；t_j为探地雷达剖面记录下的电磁波传播的双程时间；t₀为地表直达波传播的双程时间；t₁为电磁波从发射天线到遇到异常体的双程传播时间。当探地雷达图像上的道数据与GPS精确匹配后，结合改进的电磁波传播路径模型可以实现探地雷达图像纵轴方向上的时深转换，进而将每个采样点的深度精确计算出来，其过程如下：

探地雷达图像第i道数据上首个采样点的坐标为：

(x_is,y_is,z_is)_GPR＝(x_i,y_i,(z_i-h_GPS))_GPR。

那么探地雷达图像上第i道数据上第j个采样点的坐标为：

(x_ij,y_ij,z_ij)_GPR＝(x_is,y_is,z_is)_GPR-h_i，

其中，t₁＝t_j-t₀，t_j＝j·Δt，j表示电磁波单道数据上的采样点，Δt为采样时间间隔。

经过时间同步和时深变换后，探地雷达图像上每个采样点都具有精确的位置信息。时深转换为探地雷达图像纵轴上双程传播时间转换为深度的过程。探地雷达原始图像在纵轴上记录的是电磁波到目标体的双程传播时间，如果电磁波的速度已知，就可以计算出异常体的深度。

步骤五：根据探地雷达数据采集方式建立了探地雷达探测空间坐标系，以三维激光点云所在坐标系为参考，建立三维激光点云与探地雷达坐标系之间的转换参数，将探地雷达图像统一到与激光相同的大地坐标系下坐标系中，实现探地雷达图像与地面三维激光点云的融合。

激光点云和探地雷达数据融合实质是求解两个点云数据集之间的对应关系，即数学上的映射关系。理论上，三维激光点云与探地雷达数据的融合就是使地面激光扫描仪和探地雷达两个传感器所有的公有点对(p_s,q_t)满足同一变换矩阵T，即

||Tp_s-q_t＝0||。

实际上，求解上面的方程是比较困难的，因为需要解决点对的查找和变换矩阵T的求解问题。由于三维激光点云与探地雷达之间无重叠部分，无法通过公共点实现两数据融合。三维激光点云与探地雷达图像数据的融合可以转化为两数据坐标系统一的问题，其中最关键的是解决探地雷达探测坐标系下采集的数据到大地空间坐标下的坐标转换。

对于探地雷达数据，要实现从探测空间坐标系到参考大地坐标系的转换，首先是建立探测空间坐标系：根据探地雷达图像数据的特点，建立探测空间坐标系o-xyz，如图4所示，以测线的起点o作为坐标原点，x轴表示探地雷达天线沿测线采集数据的方向，y轴表示探地雷达不同测线的方向，z轴表示探地雷达天线探测的方向。若坐标系O-XYZ代表大地空间坐标系，那么探地雷达探测坐标系o-xyz中任一点p(x,y,z)到坐标系O-XYZ下采样点P的转换关系可通过地理参考绝对定位方程表示：

即

式中，(X，Y，Z)为探地雷达采样点P在大地空间坐标系O-XYZ下的坐标，(x，y，z)为探地雷达在探测空间坐标系o-xyz下的采样点p的位置坐标，(x₀,y₀,z₀)为探测空间坐标系原点o在大地空间参考坐标系O-XYZ中的坐标。矩阵为两坐标系的空间位置旋转矩阵，九个系数分别为两轴系间的方向余弦。

经过坐标转换将激光点云与探地雷达图像的数据统一到相同的大地坐标系下后，为实现两者数据格式的统一，还需对探地雷达图像的数据进行属性归一化处理。激光点云的数据一般为(X,Y,Z,R,G,B)的形式，其中，X、Y、Z代表激光点云在相应大地坐标系下的坐标，R，G，B值代表其纹理颜色，取值范围为0-255之间。而转化后探地雷达图像的数据格式为(X,Y,Z,Q)，其中，X、Y、Z代表激光点云在相应大地坐标系下的坐标，Q代表采样点的电磁波瞬时振幅值，其值一般由介质的属性确定且数值较大。因此，在两者数据融合显示之前，利用式(7)对探地雷达数据的属性值进行归一化处理，使其与点云的颜色属性值处在统一量级上。

其中，Q_max代表探地雷达图像上所有采样点振幅值中的最大值，Q_min表示探地雷达图像上所有采样点振幅值中的最小值，Q表示探地雷达图像上任一采样点的振幅值，Q₀为强度归一化后探地雷达图像上任一采样点的振幅值。

地面三维激光点云与探地雷达图像之间由于无重叠部分，无法通过公共点实现两数据融合。本发明根据两传感器的数据采集方式及特点，提出通过对探地雷达图像与GPS时间同步，将探地雷达图像统一到与激光点云相同的大地坐标系下，两数据在显示软件中同时显示，就可以实现两数据的无缝可视化融合显示。

沿探地雷达测线提取均匀分布的离散特征点，直接计算离散特征点与所对应的激光点云之间差值的均值来验证融合精度：根据地面点和非地面点之间的地形坡度有一定的差别，利用地形坡度的变化关系来实现对非地面点的滤波；假设对于点A(x_A,y_A,z_A)，在给定的半径R的区域内，若该点与任一点之间坡度的最大值小于设定的坡度阀值，则点A为地面点，反之为非地面点；点A(x_A,y_A,z_A)相对于点B(x_B,y_B,z_B)的坡度S_A的计算公式为：通过最大局部坡度滤波的方法提取出激光点云上所有的地面点；沿探地雷达测线上均匀选取离散特征点为P(x_k,y_k,z_k)，那么地面上与之对应的激光特征点通过迭代计算：其中，(x,y,z)为地面上的激光特征点，激光点云与探地雷达数据离散特征点在水平和高程方向的误差可以通过上式计算出来，评估两种数据融合的效果。

根据探地雷达数据的采集方式，三维激光点云与探地雷达图像融合后的数据共有两种显示方式，一种是激光点云与探地雷达二维时间剖面显示的方式，另一种是对于重点检测区域采用等间距采集多道探地雷达二维剖面图像，选择线性的插值算法来实现探地雷达数据的三维显示，将地下三维空间信息更加直观的展示出来，从而实现点云和探地雷达图像的地上地下三维一体化显示。

空间线性的插值算法具体过程：其中P(x_p,y_p,z_p)为两已知剖面之间对应未知一点的空间坐标为

x_p＝x_a+K(x_a-x_b)

y_p＝y_a+K(y_a-y_b)

z_p＝z_a+K(z_a-z_b)

其中，P_a(x_a,y_a,z_a)和P_b(x_b,y_b,z_b)分布为两相邻探地雷达二维时间剖面。在求出未知点的空间未知后，该点的强度值Q_p可通过以下公式计算：

其中，Q_a和Q_b为两相邻探地雷达二维时间剖面上的相对应已知点的强度值。上述方程中K值则是通过根据数据采集时采样道间距Δx和二维剖面之间的距离D确定的：

(D/Δx)≤K≤D，K为整数。

为了验证三维激光点云与探地雷达数据的融合效果，选择良乡校区科技馆旁一处地下有已知管道分布的区域进行实验，实验区进行灰度处理后的图像如图5所示，采用频率为500MHz的探地雷达天线沿测线(图中箭头所示)共采集18道二维时间剖面图，采集的过程中探地雷达与GPS同步采集，剖面之间的距离为20cm。首先按照探地雷达数据后处理的一般流程进行解震荡滤波、自动增益(AGC)、背景滤波、带通滤波和图像平滑处理；对所有的探地雷达剖面进行处理，将探地雷达各道剖面利用数据编辑统一采集长度后，利用线性插值实现探地雷达图像的三维显示。三维激光点云与探地雷达数据融合后的效果图如图6所示，其中图6(a)表示激光点云与探地雷达二维剖面显示的效果图，图6(b)为激光点云与探地雷达图像的地上地下一体化显示效果图。三维激光点云与探地雷达图像的融合，可以将探地雷达图像的位置在三维空间展示出来，并结合地面上的空间信息来对地面下的异常情况进行分析。本发明也适用于车载激光扫描系统或机载激光获取的点云与探地雷达图像的融合。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (6)

1.一种地面三维激光点云与探地雷达图像融合的方法，其特征在于，在探地雷达图像与GPS数据同步采集基础上，利用改进的电磁波传播路径模型实现探地雷达图像的时深转换，以地面三维激光数据所在坐标系为参考，结合地理参考绝对定位方程将两不同源异构数据转换到同一坐标系下，实现数据融合，其步骤如下：

步骤一：地面三维激光扫描仪采集的点云数据和高分辨率影像直接融合生成彩色激光点云，进行滤波和抽稀处理后，结合三个或三个以上地面已知控制点通过七参数转换法将彩色激光点云由仪器坐标系转换到大地坐标系下；

步骤二：对探地雷达采集的二维时间剖面图像进行相应的数据后处理，后处理包括解震荡滤波、自动增益、背景滤波、带通滤波和图像平滑处理；

步骤三：根据探地雷达和GPS同步采集后的时间同步文件，将处理后的探地雷达图像与GPS坐标进行精确匹配，得到探地雷达图像竖直方向上首个采样点的坐标；

步骤四：在改进的电磁波传播路径模型基础上，根据电磁波在介质中的传播速度实现探地雷达图像纵轴方向上的时深转换，由首个采样点的坐标得到探地雷达图像上所有采样点坐标；

步骤五：根据探地雷达数据采集方式建立了探地雷达探测空间坐标系，以三维激光点云所在坐标系为参考，建立三维激光点云与探地雷达坐标系之间的参数转换模型，将探地雷达图像上所有采样点统一到同一坐标系中，实现探地雷达图像与三维激光点云的融合；

所述探地雷达与GPS同步采集的方法为：

(1).探地雷达采集的二维时间剖面图像e(x_i,t_j)，1≤i≤M,1≤j≤N，M为探地雷达图像道数，i为探地雷达图像的第i道数据，N为每道数据上的采样点数，j为第j个采样点，则探地雷达在水平距离x_i＝i·Δx，Δx为采样的道间距，探地雷达在纵轴上的时间往返信号为t_j＝j·Δt，Δt为采样时间间隔；

(2).测距轮的精度Δd＝C/N_d，其中，C为测距轮的周长，N_d为测距轮旋转一周的脉冲个数；

(3).由于探地雷达和GPS之间的数据采集是通过测距轮同步触发，探地雷达主机与GPS主机接收的脉冲数应一致，即N_GPS＝M·(Δx/Δd)，N_GPS为GPS打标文件中记录的脉冲事件个数；

(4).GPS接收机同时获取每一个外部脉冲和绝对时刻的空间位置坐标，建立采集的探地雷达数据与GPS数据的对应关系为：

(x_i,y_i,z_i)_GPS＝i·(Δx/Δd)·(x_i,y_i,z_i)_GPR，1≤i≤M,

其中，(x_i,y_i,z_i)_GPS和(x_i,y_i,z_i)_GPR分别为同一时刻探地雷达天线上的流动站GPS的坐标和探地雷达图像上第i道数据的坐标；从而获取探地雷达图像上每道数据的位置信息；

(5).由于GPS天线的中心与探地雷达天线的中心重合，坐标转换的过程中x轴和y轴方向的平移矢量为零，只需计算z轴方向的平移矢量，那么探地雷达图像第i道数据的位置信息为：

(x_s,y_s,z_s)_GPR＝(x_i,y_i,(z_i-h_GPS))_GPS；

其中，h_GPS为GPS天线到探地雷达天线中心位置的高度；(x_s,y_s,z_s)_GPR为探地雷达天线中心位置的坐标，(x_i,y_i,z_i)_GPS为探地雷达天线上的流动站GPS的坐标；

所述改进的电磁波传播路径模型包括探地雷达发射天线与接收天线间的电磁波在空气中传播的直达波和地面下的反射波，电磁波的地面下反射面的深度h为：

t₁＝t-t₀

式中，v为电磁波在介质中的传播速度；d为发射天线和接收天线之间的距离，探地雷达屏蔽天线收发天线之间的距离是固定不变的，不同频率天线的天线距也不相同；t为探地雷达剖面记录下的电磁波传播的双程时间；t₀为地表直达波传播的双程时间；t₁为电磁波从发射天线到遇到异常体的双程传播时间。

2.根据权利要求1所述的地面三维激光点云与探地雷达图像融合的方法，其特征在于，所述七参数转换法将彩色激光点云转换为大地坐标系下的方法为：

其中，(x_j，y_j，z_j)为激光点云j在激光仪器坐标系下的坐标，(X_D，Y_D，Z_D)为激光点云j在大地坐标系下的坐标；ω，γ分别为x,y,z轴对应的旋转参数，Δx，Δy，Δz分别为x,y,z轴的平移参数，m为尺度变换参数；当地面已知控制点多于3个时，通过最小二乘法求7个转换参数的最或然值。

3.根据权利要求1所述的地面三维激光点云与探地雷达图像融合的方法，其特征在于，所述探地雷达图像与GPS坐标进行精确匹配的方法是：采用探地雷达与GPS硬件集成的方式来实现探地雷达图像空间位置的获取，将GPS天线固定在探地雷达天线的中心位置进行数据的采集；当探地雷达的主机开始工作时，根据设置的采集间隔，高精度测距轮在行进的过程中会不断触发探地雷达主机采集数据，与此同时也会触发GPS天线主机的I/O口以打标文件的方式记录下此时刻的GPS时间，并存储在记录卡内；数据采集后，根据基站和流动站GPS的差分处理结果，结合道间距和测距轮的精度通过插值实现探地雷达图像上道数据与GPS时间的精确匹配，使每道探地雷达数据都具有精确的位置信息。

4.根据权利要求1所述的地面三维激光点云与探地雷达图像融合的方法，其特征在于，将探地雷达图像的数据从探测空间坐标系到与地面三维激光相同的大地参考坐标系的转换的方法是：首先，根据探地雷达图像数据的特点，建立探测空间坐标系o-xyz，以测线的起点o作为坐标原点，x轴表示探地雷达天线沿测线采集数据的方向，y轴表示探地雷达不同测线的方向，z轴表示探地雷达天线探测的方向；

其次，以地面三维激光点云所在大地空间坐标系O-XYZ为参考，将探地雷达探测坐标系o-xyz中任一采样点p(x,y,z)到坐标系O-XYZ下采样点P(X,Y,Z)的转换关系表示为：

即

式中，(X，Y，Z)为探地雷达采样点P在大地空间坐标系O-XYZ下的坐标，(x，y，z)为探地雷达在探测空间坐标系o-xyz下的采样点p的位置坐标，(x₀,y₀,z₀)为探测空间坐标系原点o在大地空间参考坐标系O-XYZ中的坐标，矩阵为两坐标系的空间位置旋转矩阵，九个系数分别为两轴系间的方向余弦。

5.根据权利要求1所述的地面三维激光点云与探地雷达图像融合的方法，其特征在于，对探地雷达图像的数据进行属性归一化处理的方法是：激光点云的数据为(X,Y,Z,R,G,B)的形式，其中，X,Y,Z代表激光点云在大地空间参考坐标系下的坐标，R，G，B值代表其纹理颜色，取值范围为0-255之间；转化后探地雷达图像的数据格式为(X,Y,Z,Q)，其中，X,Y,Z表示探地雷达采样点与激光点云对应大地坐标系下的坐标，Q代表采样点的电磁波瞬时振幅值；对探地雷达数据的属性值进行归一化处理，使其与激光点云的颜色属性值处在统一量级上：

其中，Q_max代表探地雷达图像上所有采样点振幅值中的最大值，Q_min表示探地雷达图像上所有采样点振幅值中的最小值，Q表示探地雷达图像上任一采样点的振幅值，Q₀为强度归一化后探地雷达图像上任一采样点的振幅值。

6.根据权利要求1所述的地面三维激光点云与探地雷达图像融合的方法，其特征在于，利用空间线性的插值算法实现探地雷达数据的三维显示的方法为：P(x_p,y_p,z_p)为两已知剖面之间对应未知一点的空间坐标为：

x_p＝x_a+K(x_a-x_b)

y_p＝y_a+K(y_a-y_b)

z_p＝z_a+K(z_a-z_b)

其中，P_a(x_a,y_a,z_a)和P_b(x_b,y_b,z_b)分别为两相邻探地雷达二维时间剖面；在求出未知点的空间位置后，该点的强度值Q_p为：

其中，Q_a和Q_b为两相邻探地雷达二维时间剖面上的相对应已知点的强度值；K值则是通过数据采集时采样道间距Δx和二维剖面之间的距离D确定的：(D/Δx)≤K≤D，K为整数。