CN113534140A - 基于波场互相关的探地雷达三维成像方法 - Google Patents

Abstract

本发明提出了一种基于波场互相关的探地雷达三维成像方法，能够实现对地下场景进行快速高分辨三维成像处理，计算复杂度低，满足工程实际需求。本发明探地雷达系统采用共中心模式采集回波数据，通过沿俯仰方位步进移动收发天线对得到成像场景C扫描回波数据省去了后续在时间维作傅里叶变换，降低了算法复杂度。考虑到分层介质格林函数的平移不变性，得到收发天线位置分层介质格林函数共轭相乘结果，然后进行二维傅里叶变换，降低了算法复杂度。

Description

基于波场互相关的探地雷达三维成像方法

技术领域

本发明涉及探地雷达技术领域，具体涉及一种基于波场互相关的探地雷达三维成像方法。

背景技术

探地雷达(Ground Penetrating Radar)技术是一种较新的地球物理探测方法，该技术根据高频电磁脉冲波在地下不同物性介质之间得反射、透射和绕射等来探测地下介质分布规律，在国内外已得到广泛的应用。相比于电阻率法、低频电磁感应法等常规无损探测方法，探地雷达技术具有探测效率高，可连续测量，分辨率高，操作方便灵活，探测费用低等优越性，已经成为目前浅表层目标探测的热门技术之一。目前，探地雷达技术的应用领域已经远远超出了“探地”的范畴，作为一种重要的无损探测方法，在隐藏危险物探查、工程质量检测、环境监测、地质勘查、考古以及军事反恐等领域中发挥着重要的作用。

由于GPR探测的根本目标是为了对地下目标或场景进行检测和识别，而最初人们直接使用一维或二维回波来完成检测或识别工作，这种工作方式既不方便直观，又对操作人员的专业技能和先验知识提出了较高层次的要求。为了降低目标检测和识别的难度，GPR成像技术应运而生并不断发展。通过改变探地雷达回波中能量的分布情况，GPR成像技术可以更为直观地表现地下场景及其中目标的几何或物理特征，便于目标解译工作的开展。由于GPR成像技术可使目标解译变得更加容易和准确，已成为当今非破坏性对地探测领域中应用最广泛、最具发展前景的技术之一。

根据成像过程中处理回波数据时所在域的不同，可将GPR成像算法分为时域算法和频域算法两类，其中最典型和最具代表性的算法有时域的后向投影(Back Projection,BP)成像算法和频域的距离偏移(Range Migration,RM)成像算法。BP算法起源于计算机辅助层析技术，是SAR成像中原理最直观简单的一种成像算法。由于原理类似，BP算法可直接应用到GPR成像之中，而且可以克服GPR分层介质带来的不利影响。比如[Halman J I,Shubert K A,Ruck G T.SAR processing of ground-penetrating radar data forburied UXO detection:results from a surface-based system[J].IEEE Transactionson Antennas and Propagation,1998,46(7):1023～1027.]使用BP成像算法对地下未爆炸武器进行了三维成像，结果证明了该算法在对地下目标定位和识别时的有效性。Tanyer使用时域有限差分法对非均匀介质中的BP成像问题，通过在能量求和中引入加权步骤，改善了BP算法的成像精度。[雷文太,刘立业,粟毅.平面分层媒质中目标的TAM-BP成像算法[J].信号处理,2007,23(5):680～685.]在考虑折射和衰减效应的基础上，增加了对近场条件下球面波的传播衰减的考虑，使得BP算法可以适应有耗媒质场景。从已有算法可以发现：BP算法在采样孔径不均匀、运动误差较大时鲁棒性较强，可以获得比频域算法更优的效果，但是其自身具有运算速度慢、旁瓣和杂波能量高等缺点。RM算法起源于地震信号处理，[Stolt RH.Migration by transform[J].Geophysics,1978,43(1):23～48.]首先提出了RM算法。由于原理类似，基于波方程的RM算法同样也适合于GPR成像场景为均匀介质时的情况。比如[Mast J E.Microwave pulse-echo radar imaging for the nondestructiveevaluation of civil structures[D].Urbana,Illinois:University of Illinois atUrbana-Champaign,1993.]将SAR成像中的RM算法引入GPR成像中，由标量波动方程计算出电磁波在各项均匀介质中的传播函数，提出了适合于GPR成像的RM算法；[Gu K,Wang G,LiJ.Migration based SAR imaging for ground penetrating radar systems[J].IEEProceedings of Radar Sonar Navigation,2004,151(5):317～325.]于2004年讨论了使用RM算法对下视和前视GPR进行成像的方法，通过空间补零的方法克服了假频问题。但是，以上的研究都是基于波速不变模型进行推导的，即总是假设成像场景中只存在一种介质，而在实际GPR应用时，受人员安全或设备工作环境的考虑，天线总是距离地面一定的高度，从而使得上述基于单一场景模型的RM算法失效。

逆时偏移(Reverse-time Migration,RTM)最早在地震勘探领域提出，由于其成像精度高等优势逐渐被应用于探地雷达成像[Leuschen C J,Plumb R G.A matched-filter-based reverse-time migration algorithm for ground-penetrating radar data[J].IEEE Trans.Geo.Remote Sens,2001,39(5):929-936.]。探地雷达逆时偏移的基本原理是将地表记录到的接收点电磁波场在时间轴上进行逆向传播，当电磁波场逆推至零时刻，则所有反射波和绕射波的能量都回到最初被反射和绕射的空间位置，然后应用成像条件可获得最终的偏移剖面。逆时偏移基于双程波动方程，且不存在对波动方程的近似，允许电磁波在任意方向进行传播，能够对任意类型的电磁波场进行成像，因此逆时偏移方法相比其他偏移方法具有无法比拟的优越性。但是传统逆时偏移算法需要多次的正向和反向波场外推，然后应用成像条件将每个炮点的成像结果进行累加。因此传统时域逆时偏移需要很大的计算与存储资源，导致其并不能在实际场景中对大数据量进行实时成像处理，故并未得到广泛应用。

BP和RM成像算法往往只适用于远场环境成像，在近场条件下其性能将大幅恶化。传统时域RTM算法计算复杂度过高。在地探地雷达工程实践中，使用场景非常复杂，采集数据往往不均匀。

发明内容

有鉴于此，本发明提出了一种基于波场互相关的探地雷达三维成像方法，能够实现对地下场景进行快速高分辨三维成像处理，计算复杂度低，满足工程实际需求。

为实现上述目的，本发明的技术方案为：

本发明的一种基于波场互相关的探地雷达三维成像方法，包括如下步骤：

步骤1，将探地雷达系统C扫描回波数据在俯仰方位时间维作三维傅里叶变换，得到三维傅里叶变换后的探地雷达扫描回波；

根据成像场景先验信息，计算探地雷达数据采集平面中心点处收发天线分层介质格林函数并共轭相乘，得到共轭相乘结果，然后对共轭相乘结果在俯仰、方位维作傅里叶变换，得到二维傅里叶变换后的共轭相乘结果；

步骤2，将探地雷达发射信号频谱、所述三维傅里叶变换后的探地雷达扫描回波以及所述二维傅里叶变换后的共轭相乘结果相乘，得到三者相乘结果；

将所述三者相乘结果作俯仰以及方位逆傅里叶变换得到各单频点的三维复图像；

将各单频点的三维复图像相参累加得到最终成像结果。

其中，所述步骤1中，所述探地雷达系统采用共中心模式采集回波数据，通过沿俯仰方位步进移动收发天线对得到成像场景C扫描回波数据R([x,y,z＝0],t)，其中[x,y,z＝0]表示收发天线对中点坐标，t表示时间；当探地雷达系统为频率步进体制，则接收信号为R([x,y,z＝0],ω)，ω表示角频率。

其中，所述步骤1中，对采集得到的C扫描回波数据R([x,y,z＝0],t)在俯仰方位时间维作三维傅里叶变换的具体方式为：

其中

为采集C扫描回波数据俯仰方位时间维作三维傅里叶变换的结果。

其中，所述步骤1中，所述分层介质格林函数为G(r,r₀,ω)，其中r＝(x,y,z)为三维成像场景的坐标，r₀＝(x₀,y₀,z＝0)为探地雷达数据采集平面中心点坐标；

利用分层介质格林函数的平移不变性，r_s＝(x_s,y_s,z＝0)位置处的格林函数分布为G(r,r_s,ω)＝G(r-r_s,r₀,ω)；将收发天线位置分层介质格林函数共轭相乘的表达式为：M(r,r₀,ω)＝G(r,r₀,ω)G^*(r+r_d,r₀,ω)

其中*表示复共轭，r_d表示收发天线对间距；

将M(r,r₀,ω)作俯仰方位二维傅里叶变换的表达式为：

其中

为收发天线位置分层介质格林函数共轭相乘后作二维傅里叶变换结果。

其中，所述步骤2中，探地雷达数据采集中心位置处分层介质格林函数共轭相乘结果与作三维傅里叶变换后的探地雷达扫描回波、发射波形频谱三者相乘的表达式为：

其中S(ω)为探地雷达发射信号s(t)的频谱，B([k_x,k_y,z],ω)为三维复图像作俯仰方位傅里叶变化结果。

其中，对B([k_x,k_y,z],ω)作俯仰方位逆傅里叶变换，得到地下场景各频点三维复图像，所述表达式为：

其中

为地下场景各频点的三维复图像。

其中，各频点复图像相参累加得到最终成像结果的表达式为：

其中，I(x,y,z)为最终地下场景三维成像结果。

其中，所述步骤1中，采用非均匀傅里叶变换将回波数据变换到均匀波数域，然后再进行三维傅里叶变换。

其中，所述步骤1中，所述成像场景先验信息包括被探测地下介质的相对介电常数以及天线口面离地面高度。

其中，所述步骤1中，在俯仰、方位维作傅里叶变换的计算根据成像场景先验信息提前计算预存在磁盘中。

有益效果：

本发明考虑到分层介质格林函数的平移不变性，得到收发天线位置分层介质格林函数共轭相乘结果，然后进行二维傅里叶变换，降低了算法复杂度。本发明成像分辨率高，算法稳健性强；适用于远近场场景以及分层介质场景成像；算法复杂度低，收敛性稳定，能对探地雷达实测数据实时处理。

本发明探地雷达系统采用共中心模式采集回波数据，通过沿俯仰方位步进移动收发天线对得到成像场景C扫描回波数据省去了后续在时间维作傅里叶变换，进一步降低了算法复杂度。

本发明采用非均匀傅里叶变换将回波数据变换到均匀波数域，只需精确记录收发天线位置的坐标即可，进一步降低了算法复杂度。

本发明中，二维傅里叶变换后的共轭相乘结果的计算根据成像场景先验信息提前计算预存在磁盘中，大幅提升所提算法数据处理效率，满足工程实际需求。

附图说明

图1是本发明的探地雷达系统数据采集方式示意图

图2是本发明的算法流程图；

图3是本发明探地雷达仿真回波数据的三维场景示意图；

图4是本发明图3三维模型剖面图；

图5是本发明对仿真场景回波数据三维成像结果。

具体实施方式

下面结合附图并举实施例，对本发明进行详细描述。

本发明提出了一种基于波场互相关的探地雷达三维成像方法，技术原理为：

本发明使用共偏移距模式采集GPR数据，即收发天线之间距离恒定。基于波场零延时互相关成像条件，探地雷达成像结果可以表示为：

其中r＝(x,y,z)为成像网格点坐标，对于二维垂直截面和水平切片成像则网格坐标分别取r＝(x,z)，r＝(x,y)；r_d＝(x_d,y_d,z_d＝0)是收发天线之间的坐标差值，并且天线仅在xoy平面以r_d步进移动；u_r(·)与u_s(·)分别是收发天线位置的后向与前向递推波场。对每个天线收发位置的数据进行时域逆时偏移并进行累加便得到最终成像结果I(r)。

根据卷积定理：

其中U_s(·)与U_r(·)分别为前向与后向时域波场值u_r(·)、u_s(·)的频域形式。根据电磁场理论，波场分布是分层介质格林函数与相应收发位置处激励源频谱的乘积，即：

U_s(r,r_s,ω)＝G(r,r_s,ω)S(ω)

其中G(·)频域分层介质格林函数，上标*表示复共轭，S(ω)为探地雷达发射信号s(t)的频谱，

为接收信号时间反转后的傅里叶换。

考虑到分层介质格林函数的平移不变性，即：

G(r,r_s,ω)＝G(r-r_s,r₀,ω)

其中r₀＝(0,0,0)为数据采集平面的中心坐标，G(r,r₀,ω)为对应坐标处的频域分层介质格林函数。则：

其中，M(r,r₀,ω)＝G(r,r₀,ω)G^*(r+r_d,r₀,ω)

在某些远场探测场景中，收发天线间距很小，式(12)可简化为M(r,r₀,ω)＝G(r,r₀,ω)G^*(r,r₀,ω)

根据卷积定理：

其中，

分别为M(r,r₀,ω)、

对x、y方向作二维傅里叶变换的结果。I(r)为最终三维成像结果。

本发明算法流程如图2所示，本发明包括如下步骤：

步骤1，将探地雷达系统C扫描回波数据在俯仰方位时间维作三维傅里叶变换，得到三维傅里叶变换后的探地雷达扫描回波；

所述探地雷达系统采用共中心模式采集回波数据，数据采集方式示意图如图1所示。通过沿俯仰方位步进移动收发天线对得到成像场景C扫描回波数据R([x,y,z＝0],t)，其中[x,y,z＝0]表示收发天线对中点坐标，t表示时间；当探地雷达系统为频率步进体制，则接收信号为R([x,y,z＝0],ω)，ω表示角频率，省去了后续在时间维作傅里叶变换；

对采集得到的C扫描回波数据R([x,y,z＝0],t)在俯仰方位时间维作三维傅里叶变换的表达式为：

其中

为采集C扫描回波数据俯仰方位时间维作三维傅里叶变换的结果，在实际数据采集过程中，数据采集点往往非均匀分布，可以采用非均匀傅里叶变换将回波数据变换到均匀波数域，只需精确记录收发天线位置的坐标即可。

根据成像场景先验信息，计算探地雷达数据采集平面中心点处收发天线分层介质格林函数并共轭相乘，得到共轭相乘结果，然后对共轭相乘结果在俯仰、方位维作傅里叶变换，得到二维傅里叶变换后的共轭相乘结果。

本实施例中，根据被探测地下介质的相对介电常数，以及天线口面离地面高度等成像场景的先验信息，计算探地雷达数据采集平面中心点处分层介质格林函数G(r,r₀,ω)，其中r＝(x,y,z)为三维成像场景的坐标，r₀＝(x₀,y₀,z＝0)为探地雷达数据采集平面中心点坐标；

利用分层介质格林函数的平移不变性，r_s＝(x_s,y_s,z＝0)位置处的格林函数分布为G(r,r_s,ω)＝G(r-r_s,r₀,ω)；将收发天线位置分层介质格林函数共轭相乘的表达式为：

M(r,r₀,ω)＝G(r,r₀,ω)G^*(r+r_d,r₀,ω)

其中*表示复共轭，r_d表示收发天线对间距；

将M(r,r₀,ω)作俯仰方位二维傅里叶变换的表达式为：

其中

为收发天线位置分层介质格林函数共轭相乘后作二维傅里叶变换结果，该步骤

的计算可以根据成像场景先验信息提前计算预存在磁盘中，大幅提升所提算法数据处理效率，满足工程实际需求。

步骤2，将探地雷达发射信号频谱、所述三维傅里叶变换后的探地雷达扫描回波以及所述二维傅里叶变换后的共轭相乘结果相乘，得到三者相乘结果；具体地，探地雷达数据采集中心位置处分层介质格林函数共轭相乘结果与作三维傅里叶变换后的探地雷达扫描回波、发射波形频谱三者相乘的表达式为：

其中S(ω)为探地雷达发射信号s(t)的频谱，B([k_x,k_y,z],ω)为三维复图像作俯仰方位傅里叶变化结果；

将所述三者相乘结果作俯仰以及方位逆傅里叶变换得到各单频点的三维复图像；对B([k_x,k_y,z],ω)作俯仰方位逆傅里叶变换，得到地下场景各频点三维复图像，所述表达式为：

其中

为地下场景各频点的三维复图像。

将各单频点的三维复图像相参累加得到最终成像结果，各频点复图像相参累加得到最终成像结果的表达式为：

其中，I(x,y,z)为最终地下场景三维成像结果。

本发明的效果可以通过以下仿真结果进一步说明：

仿真条件：采用探地雷达正演软件gprMax来模拟真实场景回波数据，激励信号为中心频率为1.5GHz，-3dB带宽约为1.3GHz的Ricker脉，仿真场景为2m*2m*0.7m的三维半空间模型，地下介质的相对介电常数为6，相对磁导率为1，电导率为0.001S/m。

仿真实验一：图3给出了仿真三维场景示意图，目标为间距不同的3对半径为3cm的金属圆柱，其具体参数如图4所示。

图5为所提算法对仿真场景回波数据进行处理得到的三维成像结果，结果证明所提方法具有突出的地下目标重构能力，很好还原了地下目标的几何形态与物理尺寸，能够为地下目标探测提供技术支撑。此外由于所提算法基于波场传播特性进行推导，能够很好适用于近场环境成像。通过与BP三维成像算法进行比较，所提算法具有更低的算法复杂度，适用于工程实际中探地雷达数据的实时成像处理。

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

Claims (10)

1.一种基于波场互相关的探地雷达三维成像方法，其特征在于，包括如下步骤：

步骤1，将探地雷达系统C扫描回波数据在俯仰方位时间维作三维傅里叶变换，得到三维傅里叶变换后的探地雷达扫描回波；

根据成像场景先验信息，计算探地雷达数据采集平面中心点处收发天线分层介质格林函数并共轭相乘，得到共轭相乘结果，然后对共轭相乘结果在俯仰、方位维作傅里叶变换，得到二维傅里叶变换后的共轭相乘结果；

步骤2，将探地雷达发射信号频谱、所述三维傅里叶变换后的探地雷达扫描回波以及所述二维傅里叶变换后的共轭相乘结果相乘，得到三者相乘结果；

将所述三者相乘结果作俯仰以及方位逆傅里叶变换得到各单频点的三维复图像；

将各单频点的三维复图像相参累加得到最终成像结果。

2.如权利要求1所述的方法，其特征在于，所述步骤1中，所述探地雷达系统采用共中心模式采集回波数据，通过沿俯仰方位步进移动收发天线对得到成像场景C扫描回波数据R([x,y,z＝0],t)，其中[x,y,z＝0]表示收发天线对中点坐标，t表示时间；当探地雷达系统为频率步进体制，则接收信号为R([x,y,z＝0],ω)，ω表示角频率。

3.如权利要求2所述的方法，其特征在于，所述步骤1中，对采集得到的C扫描回波数据R([x,y,z＝0],t)在俯仰方位时间维作三维傅里叶变换的具体方式为：

其中

为采集C扫描回波数据俯仰方位时间维作三维傅里叶变换的结果。

4.如权利要求2或3所述的方法，其特征在于，所述步骤1中，所述分层介质格林函数为G(r,r₀,ω)，其中r＝(x,y,z)为三维成像场景的坐标，r₀＝(x₀,y₀,z＝0)为探地雷达数据采集平面中心点坐标；

利用分层介质格林函数的平移不变性，r_s＝(x_s,y_s,z＝0)位置处的格林函数分布为G(r,r_s,ω)＝G(r-r_s,r₀,ω)；将收发天线位置分层介质格林函数共轭相乘的表达式为：

M(r,r₀,ω)＝G(r,r₀,ω)G^*(r+r_d,r₀,ω)

其中*表示复共轭，r_d表示收发天线对间距；

将M(r,r₀,ω)作俯仰方位二维傅里叶变换的表达式为：

其中

为收发天线位置分层介质格林函数共轭相乘后作二维傅里叶变换结果。

5.如权利要求4所述的方法，其特征在于，所述步骤2中，探地雷达数据采集中心位置处分层介质格林函数共轭相乘结果与作三维傅里叶变换后的探地雷达扫描回波、发射波形频谱三者相乘的表达式为：

其中S(ω)为探地雷达发射信号s(t)的频谱，B([k_x,k_y,z],ω)为三维复图像作俯仰方位傅里叶变化结果。

6.如权利要求5所述的方法，其特征在于，对B([k_x,k_y,z],ω)作俯仰方位逆傅里叶变换，得到地下场景各频点三维复图像，所述表达式为：

其中

为地下场景各频点的三维复图像。

7.如权利要求6所述的方法，其特征在于，各频点复图像相参累加得到最终成像结果的表达式为：

其中，I(x,y,z)为最终地下场景三维成像结果。

8.如权利要求1-7任意一项所述的方法，其特征在于，所述步骤1中，采用非均匀傅里叶变换将回波数据变换到均匀波数域，然后再进行三维傅里叶变换。

9.如权利要求1-7任意一项所述的方法，其特征在于，所述步骤1中，所述成像场景先验信息包括被探测地下介质的相对介电常数以及天线口面离地面高度。

10.如权利要求1-7任意一项所述的方法，其特征在于，所述步骤1中，所述二维傅里叶变换后的共轭相乘结果的计算根据成像场景先验信息提前计算预存在磁盘中。

Images