探地雷达成像算法研究报告

1、-. z.探地雷达成像算法研究摘要探地雷达(Ground Penetrating Radar,简称GPR)集无损检测、穿透能力强、分辨率高等众多优点而成为检测和识别地下目标的一种有效技术手段。性能优良的探地雷达成像方法有助于准确定位地下目标，同时提高对目标的检测和识别能力，从而推动探地雷达在城市质量监控、地质灾害、考古挖掘、高速公路无损检测、地雷探测等各个方面得到更广泛的应用。本文以中国电波传播研究所的探地雷达LD-2000为实验设备，从中读取探测数据。以MATLAB为软件平台，实现了探地雷达数据的显示、处理、成像几个局部。其中数据显示方式包括数据的波形堆积图，剖面面色阶图以及带数据波形图；数

2、据处理局部包括直达波的去除、背景噪声的去除、振幅增益等；雷达成像算法局部主要采用波前成像算法和投影层析成像算法。Imaging Algorithm of Ground Penetrating Radar ABSTRACTGPR (Ground Penetrating Radar, referred GPR) set of non-destructive testing, penetration ability, many advantages of high resolution detection and identification of underground and bee the t

3、arget of an effective technical means. E*cellent performance GPR imaging approach helps pinpoint undergroundtargets, while increasing the target detection and identification capabilities, thereby promoting the quality of ground penetrating radar surveillance in the city, geological disasters, archae

4、ological e*cavation, highway nondestructive testing, mine detection, etc. aspects to be more widely used.In this paper, China Institute of Radiowave Propagation GPR LD-2000 for the e*perimental apparatus, reads probe data. MATLAB as the software platform to achieve a ground-penetrating radar data di

5、splay, processing, imaging several parts. Wherein the data includes a data waveform display stacked, with a cross-sectional side view and a gradation data waveform; data processing section includes the removal of the direct wave, the background noise removal, the amplitude gain, etc.; radar imaging

6、algorithm some of the major imaging algorithm and the wavefront projection tomography algorithms.-. z.1 绪论11 选题的背景及意义雷达是利用电磁波探测目标的电子设备。发射电磁波对目标进展照射并接收其回波，由此获得目标至电磁波发射点的距离、距离变化率径向速度、方位、高度等信息。而探地雷达则是利用电磁波在地下媒质中的传播与反射特性进展地下不可见目标体或界面的探查与定位分辨。探地雷达实质是向地下发射电磁波，通过电磁波的回波来分析地下的情况2。由于地下的物质与空气相比反射情况会复杂的多，因此，与雷达

7、相比，探地雷达的回波很大程度上会受到各种噪声的干扰。探地雷达在工作的过程中不断的接收回波的数据信息，数据量十分庞大，同时，接收的回波信息中，不可防止的会收到外界或者仪器本身产生的干扰2。所以在一般情况下，未经过处理的信号几乎无法判别，因此探地雷达处理准确的可视化成像图形分析越来越受到使用者的重视。再者，由于探地雷达现在的应用已经开拓了很大的应用领域，尤其是在工程地质领域的应用。探地雷达开拓应用领域的优点在于: 探地雷达是一种非破坏性的探测技术，可以平安地用于城市和正在建立中的工程现场。工作场地条件宽松，适应性强 (对于轻便类的仪器)；抗电磁干扰能力强，可在城市各种噪声环境下工作，环境干扰影响小

8、；具有工程上较满意的探测深度和分辨率。现场直接提供实时剖面记录图，图像清晰直观；便携微机控制数字采集、记录、存储和处理。轻便类仪器现场仅需3人或更少人员工作，工作效率高。本文通过系统研究和阅读大量文献资料，着重的研究了从数据的提取到，到探地雷达回波信号的预处理，再到探地雷达成像算法的研究。12 探地雷达成像算法的开展探地雷达Ground Penetrating Radar是利用电磁波探测地下目标，通过分析电磁信号与地下目标的相互作用，提取目标的性质、形状等信息1。由于地球表层和人们生产生活的密切关系，电磁波在地下媒质中的传播与反射比空间更为复杂，有许多新的课题需要研究解决，所以探地雷达技术愈来

9、愈受到人们的重视，得到了迅速的开展与提高，电磁波探测地下目标成为电磁波应用的一个重要领域。国外兴旺国家自20世纪60年代以后，探地雷达技术得到迅速开展，国自20世纪90年代以来也开场重视和开展探地雷达技术的研究和应用，并开发出了实用的产品。但是国产品在分辨率、使用方便性、对雷达信号成像和图像解译技术等方面与国外产品存在差距，因此，国探地雷达的应用中绝大多数采用国外的产品。当前随着探地雷达技术的飞速开展，先进的高分辨数据处理和成像技术成为探地雷达技术开展的关键，成像方法也趋于多样化。探地雷达合成孔径聚焦成像技术自上个世纪90年代初以来已得到逐步应用，同时基于雷达波和地震波在运动学上的相似性，反射

10、地震学上的波动偏移成像技术也逐步应用于探地雷达数据处理和成像1。上个世纪70年代初Standford大学的JClaerbout教授首先提出了用有限差分法解单程波动方程的近似式，用地面观测的地震数据重建地震波在地下传播过程中的波场，从这些传播过程的波场中提取地震界面剖面像的数据，组成地震偏移剖面，这种成像方法即为有限差分偏移成像。在上个世纪70年代后期，RStolt和JGazdag等人又先后提出了基于波动方程偏移的Stolt偏移成像方法和PhaseShift偏移成像方法，由于此两种方法在计算中可以充分利用快速傅立叶交换，因此计算简单，效率高，很快得到推广。波动方程偏移成像在最近十年间迅速开展并不

11、断完善，许多人对此做出了有益的奉献。其中，Loewenthal等人的爆炸反射面的概念，Hubral和Lamer等人的深度偏移的概念，AJBerkhout提出的偏移过程是一个空间卷积的概念，我国的马在田院士的高阶方程的分裂算法等都为波动方程偏移成像技术的开展作出了奉献，同时促进了探地雷达成像技术的进一步开展。当前仍有许多学者还在探索波动方程偏移成像，以期更加完善该方法，这也必将为探地雷达成像技术的开展注入新的活力。最近几年，应用于随机不均匀介质中的时逆(Time Reverse)方法成为国际上研究的热点。时逆方法最早是由法国巴黎大学的Mathias Fink教授提出，并在超声医学中进展研究，其可

12、以分为物理时逆和虚拟或计算时逆3。在物理时逆中，逆传播场是由发射天线向未知的真实介质辐射的；而在虚拟或计算时逆中，逆传播是在一个虚拟的参考介质中进展数值仿真。物理时逆在医学上已有许多应用，如肾结石的粉碎等，最近也出现了把物理时逆应用到如地雷引爆、通信等方面的探讨。而虚拟时逆主要应用于成像，最近两年，由Rice大学的Liliana Borcea教授、Standford大学GPapanicolaou教授、美国Lawrence Livermore国家实验室的JBerryman等学者组成的研究小组对随机介质中时逆成像方法的理论进展了大量深入的研究，并基于超声对此理论进展了大量的数值仿真，为时逆成像方法

13、在各个领域的应用提供了坚实的根底。最近，Northeastern大学的Anthony JDevaney教授在Standford大学作报告时首先开场了时逆成像方法在探地雷达对地下目标成像中的探讨，为探地雷达成像技术开辟了新的方向4。当前，探地雷达成像技术由于受到各自数学模型的影响，都有自己的缺乏和优点。成孔径聚焦成像方法适用条件广，但计算复杂；Stolt偏移成像方法和PhaseShift偏移成像方法计算简单，但无法适用于复杂的地下构造，复杂的速度模型；时逆成像方法要求地下介质有随机变化，但变化幅度必须限制在一个很小的围5。另外，探地雷达像根本上是一个二维截面图，探地雷达三维成像技术还处于研究之中

14、。探地雷达成像方法作为探地雷达的关键技术之一，其开展方向是高分辨力、高效率、高精度成像。在此要求之下，在国外大量科研人员的不懈努力下，探地雷达技术将不断提高，其应用围也将逐步拓展。在我国，探地雷达除了用于各种建筑和公路质量监控外，还将在国防和国家平安部门有很大的应用前景，现实的或潜在的应用还包括：地雷探测、地下掩体的探测、货物平安检查、打击毒品走私等方面7。13 本文的研究容本文的主要研究容如下：绪论本章介绍了课题的选题背景和原因。介绍了探地雷达开展过程，探地雷达成像算法的国外开展状况以及今后的开展趋势。 探地雷达根底本章是从电磁波的根底理论出发，讲述了电磁波在岩石介质中的传播特性。同时分析了

15、探地雷达的工作原理和探地雷达回波模型，研究了探地雷达的三种数据显示方式A-scan,B-scan,C-scan，三种显示方式Wiggle曲线，灰度图，彩色图。第三章 探地雷达数据预处理探地雷达所接收到的信号十分复杂，脉冲在通过地下介质的过程中，波形和波幅将发生较大的变化，而脉冲余振、系统部干扰、地表不光滑或地下介质不均匀等引起的散射以及剖面旁侧的绕射等干扰，均使得实时记录图像多变和不易分辨。为此我们需要先对我们所测量的数据进展一些处理。做奇异值分解以去除探地雷达最大的干扰成分直达波；做道增益以补偿电磁波在传播的过程的衰减；对探地雷达数据做中值滤波以去除背景噪声的影响。第四章 探地雷达成像算法的

16、研究本章开场介绍的两种探地雷达的扫描方式合成孔径，实孔径。然后主要介绍基于合成孔径扫描方式下两种成像算法波前成像算法，投影层析成像。投影层析成像与CT类似，通过接收空间各方位的散射场，就可以通过层析技术实现目标的反演。只是在CT中接收的是前向透射场，而在雷达扫描中接收的是后向散射场。目标函数可以用逆拉顿变换得到8。波前重建最早作为地震波迁徙技术出现在地震学中。基于地震波与电磁波传播的相似性，Cafforio和Soumekh提出了两个根本一样的SAR模型。在系统模型中Cafforio和Soumekh运用一阶波恩近似，假定目标电磁散射可以模型化为目标各个离散散射中心相应的矢量和8。2 探地雷达根底

17、21 探地雷达理论根底211 电磁场根底理论探地雷达采用高频电磁波进展探测，电磁波的传播满足麦克斯韦方程组，即2.1式中，E为电场强度(),B为磁感应强度(T), H为磁场强度()J为电流密度(), D为电位移( ),p为电荷密度()。式2.1中，各个方程依次分别为微分形式的法拉第电磁感应定律，安培电流环路定律，磁场无源定律以及电场高斯定律。当求解E, B, H和D四个未知矢量时，仅己知J和P是不够的，需要引入介质的本构关系。本构关系就是场量之间的关系，它决定于电磁场所在介质中的性质。对于均匀、线性和各向同性介质而言，它的本构关系可以简化为2.2式中，为电导率(S/m)，为介电常数(F/m)

18、,磁导率(H/m ) 。结合介质的本构关系，麦克斯韦方程组可以改写为只含两个矢量场的形式即2.3麦克斯韦方程组说明，随着时间变化的磁场会产生随时间变化的电场，随时间变化的电场会产生随时间变化的磁场。麦克斯韦方程组提醒的现象简单的说就变化的磁场和变化的电场相互激发，并且变化的磁场和变化的电场以一定的速度向外传播，这就形成了电磁波9。212 电磁波在岩土介质中的传播电磁波根据其波面的形状可以分为平面波、柱面波和球面波，其中平面波是最根本、最具有电磁波普遍规律的电磁波类型。探地雷达所发射的的电磁波可以经傅立叶变换换算一系列的谐波，这些谐波近似为平面波，则探地雷达电磁波传播以平面谐波的传播规律为根底。

19、在实际应用中，介质为有损介质，即电导率。根据欧姆定律，麦克斯韦方程组可写成如下复数形式2.4其中，为介质中的等效复介电常数。由此可以得出如下波动方程: 2.5其中，为复数，令，则可得到衰减常数、相位常数。2.6在探地雷达应用中，我们通常比拟关心电磁波的传播速度和衰减因子。其中电磁波的相速度并不是一个常数，而是频率的函数，即，即2.7在描述导电介质中均匀平面波特性的公式中，称为损耗角正切，表示了导电电流和位移电流幅度之比。当时，只有位移电流则为理想介质;当时，传导电流远远小于位移电流为低损耗介质;当时，传导电流远远大于位移电流为良导体。损耗角正切与三者都密切的关系，即便在同一介质中，对于不同频率

20、的电磁波，其表达出来的介电常数与电导率会随之变化。假设介质为低损耗介质，即，此时2.8此时，平面波的电场强度近似等于磁场强度。大多数岩石介质为非磁性、非导电介质，其传播速度为，其中为光速，为相对介电常数。此时，电磁波的速度主要取决于介质的介电常数。衰减常数与电导率成正比，与介电常数的平方根成反比，电磁波能量的衰减主要是由于感生涡流损失引起的。假设介质为良导体，即，此时2.9此时，随着电导率、磁导率增加，以及电磁波频率升高，电磁波的衰减越快。波速与频率的平方根成正比，与电导率的平方根成反比，波速是频率和电导率的函数7。213 探地雷达工作原理探地雷达利用主频为数十兆赫兹至上千兆赫兹的高频电磁波，

21、以高频带短脉冲的形式由T*向地下发射，经由地下目标体或者地层反射回来并由接收线圈R*所接收。当相邻介质的电磁特性有差异时，电磁波在接触面附近发生电磁波的反射与折射等现象，每次遇到存在电磁特性差异的接触界面电磁波均发生反射与折射。原理如下列图2.1所示图2.1 探地雷达探测示意脉冲波双程走时为，对于雷达波速，可以采用近似求解。电磁波垂直入射的情况下，反射系数R由接触界面两侧介质的相对常数，确定，如下2.10探地雷达移动发射和接收天线的同时，接收到反射电磁波的双程走时相应变化，如下列图2.2所示为同步移动发射和接收天线，对地下目标体反射波接收的原理图及回波曲线图。图2.2 探地雷达探测原理及对应的

22、回波曲线214 探地雷达回波模型根据电磁波传播和反射理论，对于近距离的探地雷达，接收信号常常由四个局部组成:天线串扰、地面反射、目标反射和白噪声。模型如下列图所示：图2.3 探地雷达接收信号模型图因此，接收信号可公式化为:。公式中各个变量分析：C 是发射天线和接收天线之间的串扰信号。一般情况下，当发射天线和接收天线之间位置固定后，信号C可以认为是一个确定的信号。对于这个干扰信号可以通过天线校正或者在发射天线和接收天线之间加屏蔽隔板进展消除。b信号是地面反射信号及地下各种杂质干扰带来的噪声信号。由于天线和地面之间的阻抗不匹配，所以引起电磁波的屡次反射。地面的第一次反射信号最强，以后逐渐衰减。对于

23、探测浅地层的探地雷达来说，这种屡次反射造成的影响是较大的，地下各种杂质带来的干扰统称为背景噪声。都是在杂波抑制中要去除的信号。S信号是目标反射的信号，也是真正需要的信号。这是一个与目标体相关确实定信号，由于与到达的时间和目标体的位置有关，因此它也是一个变化的信号；n是噪声项，它由测量噪声和模型噪声组成，可以假定它是一个加性高斯白噪声4。本文中的杂波抑制算法局部即是根据上述这些杂波的特性，完成一个信号处理的算法对上述的杂波干扰给予有效的去除。22 探地雷达数据显示探地雷达的工作方式是在水平方向的一个位置上，发射天线对地下发射脉冲信号，接收天线接收在该位置上的回波信号，以此作为一列的回波数据;然后

24、收发天线移动到下一位置，发射天线再对地下发射脉冲信号，接收天线再接收在该位置上的回波信号并作为新的一列回波数据;以此类推，等间隔的移动N个位置，便可接收到N列的回波数据，即得到一次扫描的数据。如下列图所示，是在一个位置上接收到一列回波信号的情况，称之为A扫描(A一Scan)。图2.4 探地雷达A-scan示意图图2.5 探地雷达A-scan数据图下列图所示的是在等间隔的N个位置上接收到的N列回波信号的情况，称之为B扫描(B-Scan)。图2.6探地雷达B-scan示意图由上图所示可知，一个B扫描是由N个A扫描组成，当对一个B扫描的回波信号进展抽样量化后，便得到探地雷达的原始数据，它是一个的矩阵

25、。(M是对每个A扫描数据的采样点数，N是包含的A扫描的个数)。假设对一个B扫描的图像直接成像，将在所成的像中看到一条由目标反射所形成的围很宽的双曲线带。本文的算法针对的是B-Scan的数据进展数据处理和成像。假设数据采集在*0y平面，沿多条平行的测量线，完成对*一区域的测量，这种扫描方式为C-scan。C-scan数据是由多条测量线的Bscan数据所构成，可对媒质部形成立体全面的扫描8。23 探地雷达的显示方式231 Wiggle曲线的显示方式在分析探地雷达信号过程中，Wiggle曲线图能从整体上、直观的解释不同深度的接收到回波的变化情况，是解释地质情况的重要分析方法。在一个探地雷达的二维成像

26、图中，如果能够通过曲线描绘出各个通道振幅变化情况，则这样的曲线就是Wiggle曲线。实际上，Wiggle曲线是并不是由一条曲线而成的，每通道的曲线都是由许多小线段组成。线段的横坐标表示通道数，纵坐标表示时间深度，幅值则由数据值的大小来表示。Wiggle图能较为直观的看出通道波形的变化以及同相轴的位置2。Wiggle曲线的算法可以表述如下。探地雷达的数据值，表示幅值的大小，正负表示方向。并且在二维图像的位置不会由于数值大小的改变而发生位置改变的现象。根据这个特点，可以将每通道振幅归一到-1 ,1的围。归一化的方法是，遍历各个通道数据，找出各通道数据中绝对值最大的数据，与其对应的通道数据都以这个最

27、大值作为基准值进展相应的计算。己知，*采样点的值为，则归一化后的值为。显示Wiggle曲线，需要计算各道曲线的间隔。设定各个根本参量和其对应的计算表达式如下表2.1所示。表2.1 wiggle曲线坐标设置根本参数计算表达式设置显示的*轴长度Ntraces设置显示的Y轴长度Nscale横轴各道间隔Trace = ntraces/tracescount纵轴各采样点间隔Avsample = nscale/samplecount各点*轴坐标Trace*itrance + trace*Dshow各点Y轴坐标Avsample*jscale将各坐标点依次连接，在MATLAB中，主要是通过plot指令完成曲线

28、的绘制。图2.7 探地雷达wiggle图显示232 灰度图显示式在MATALB中，inshow (Data)用于显示导入的雷达采样数据图。在没指定颜况下，默认颜色映像是jet。colormap)函数用于设定颜色映像，colormap (gray)用于设定常用的灰度图显示方式。灰度图显示方式是以各个色度的大小来区分各个扫描道的雷达波的振幅值大小，同时，各个扫描道按照一定的位置以及时间的间隔进展顺次排列，显示出渐变且均匀的灰度图。图2.8 探地雷达灰度图显示233 彩色显示方式彩色显示方式是利用不同的色彩代表不同的扫描道的雷达波的振幅值。与灰度图显示方式一样，各个扫描道按照一定的位置和时间的间隔进

29、展顺次排列，显示出渐变、均匀并且无断层现象的彩色显示方式。为了能够更加直观表示雷达图像中的目标区域，在绘图中，用两种颜色分别表示雷达数据中的最大值和最小值，定义渐变的变化模式。彩色显示也是本研究中主要的显示方式。显示效果如下列图图2.9 探地雷达数据彩色显示探地雷达数据的预处理在探地雷达成像过程中，数据的预处理和成像算法是密不可分的两局部。本文的主要流程也是按照这两局部展开。由下列图为本文的主要流程图图3.1 探地雷达主要流程图从图中可以看出预处理是成像处理的根底。一般来说，信号预处理的目的是减少杂波对目标回波的干扰，尽可能的恢复目标回波的数据。在原始探地雷达信号中，由于干扰噪声的存在，导致目

30、标的反射回波无法得到有效的显示。因此，使得对于地下目标的探测变得非常困难，给探地雷达工作者带来了不必要的问题。针对探地雷达的噪声分析，发现各类的干扰噪声信号，都具有其自身的特点。可以通过相应的方法去除或抑制那些非目标的反射波信号，使目标回波得到增强，使得目标的得到有效的探测出来。31 探地雷达直达波的去除根据去除方法的不同，可以将直达波去除方法分为两类。一类是依据直达波与目标回波的不同特性，在不同的域直接滤除的方法，主要有时间窗去除法、二维滤波法、小波变换去除法、子空间投影去除法、主元分析去除法等。另一类是根据接收信号先估计出直达波的参考波形，然后在接收信号中减去参考波形，以到达抵消直达波的目

31、的，可统称为抵消法，这一类主要有参考波抵消法、平均抵消法、参数模型抵消法、信号模型去除法、自适应抵消法等。本文主要研究的是主元分析去除法作为直达波的去除方法8。如果根据直达波和目标回波目标回波在二维扫描图像上的差异，将其投影到不同的子空间中，在去除直达波所在子空间后，重构信号可以实现直达波的去除，这种方法即为子空间投影去除法。利用奇异值分解SVD、主元分析PCA、独立元分析ICA等方法，可对接收信号矩阵进展分解，实现子空间投影。本文所研究的是基于奇异值分解的主元分解的主元分析去除法9。主元分析方法PCA又被称为主成分分析方法，是一种建立在最小均方误差根底上的线性变换处理方法。作为一种多元统计方

32、法，主元分析可以将一组变量的原有特征进展组合，抽取所谓的主元作为新的特征，并保证熵到达极小值。实际应用中，探地雷达接收到的B-scan信号可以表示为矩阵，其中M为*测线的测量点数目，N为时间上的采样点数。对进展奇异值分解：3.1其中，为正交矩阵，其中列向量分别为对称矩阵和的特征向量。为对角阵，对角元素即为其奇异值，按从大到小排列，可以表示为。令，i=1,2,3N，即为主元，对应时间采样信息。接收信号矩阵就可以认为是主元与相应特征图像的的加权和，即。由于为正交矩阵，因此。接收信号矩阵可以表示为,其中与分别为目标回波信号矩阵和直达波信号矩阵。由于在同一测线的各个测量点直达波信号变化不是很大，因此在

33、理想情况下可以近似认为个测量点直达波一样，即直达波矩阵秩为1。而对于地下目标的反射回波信号会变化，在B-scan图上一般多呈双曲线形状，因此目标反射回波矩阵的秩大于1。由 3.2可知，当=1时，即为矩阵的最近秩1矩阵，由于直达波信号能量远远大于目标回波信号，矩阵在中占主导地位，其秩为1，及包含的全部或者绝大局部信息，仅含有少许的目标回波信号的信息，通过去除可以有效去除直达波，并较好的保存目标回波信号的信息，通过去除可以有效去除直达波，并较好地保存目标回波的信息。图3.2探地雷达去直达波后图像32 探地雷达背景噪声的去除背景噪声的产生，是由于多种因素引起的，包括地下介质，空气中电磁波以及人为使用

34、设备不当等因素。由于背景噪声的存在掩盖了回波中的重要信息，加强了获得地下信息的难度。其中，对浅层回波信号影响最大、干扰性最强的背景噪声是雷达信号在地面的屡次强反射回波信号10。由于地表回波中背景噪声的延时情况是一致的，因此在各道同一时间的采样数据值十分接近，而真实的回波在各道不同列的延时数据是有所区别的。通过以上分析，从时间的角度可以将背景噪声和真实的回波信息区分开。背景噪声被认为是比拟稳定噪声，所以背景滤波的思想可以认为是，将各道回波信号减去这个稳定的背景值。探地雷达背景噪声去除的算法描述为: 首先，对探地雷达数据进展简单的均值法处理，其目的是为了去除地表的屡次强反射。在此，均值法主要是将B

35、扫描的二维矩阵的各行中每个元素减去该行的均值，用数学方式表示为:3.3其中，仍表示原始的回波信号值，仍表示利用均值法去除背景噪声后的回波信号值，其中m为总道数，n为每到的采样点数。针对中的每通道回波进展离散傅里叶变换，其表达式为其中k=1,2，n。 3.4对求取的每个A扫描的频域信号的幅度进展中值滤波，取其区域长度为q，滤波之后为Y(n),q可以根据我们的需要取值。滤波后再进展逆离散傅里叶变换，得到。则由Z构成的回波数据，背景噪声根本得到去除，是我们想要得到的结果。该算法是在传统均值滤波方法去除背景噪声的改良。首先，事先，需要是利用均值滤波去除地外表的屡次强反射。将均值滤波后的数据，在再利用离

36、散傅里叶变换方法，在频域采用区域中值滤波的方式，作进一步的滤波处理。该方法，能够很好的去除地表的强回波信号，更好的抑制了反射波附近的杂波信号，到达的去除背景噪声的目的2。图3.3 探地雷达经滤波后的成像图像33 探地雷达回波中的增益处理电磁波在传播的过程中，会发生衰减现象。同样，探地雷达波在地下传播的过程中会发生衰减，传播的深度越深，能量衰减越大。影响衰减的主要因素包括地下介质对电磁波的吸收和传播过程中发生的损耗11。同时，电磁波在不同介质中的传播，能量衰减是不一样的，根本情况下是呈现指数衰减的，当探地雷达波到达深度时，其信号己经远远不及于探地雷达发射时的强度，因此，这对于通过解释回波的信息来

37、研究地下构造和异常是非常不利的。而且对于较深部的异常判断也是有难度的。因此通过对每道探地雷达信号进展增益调节，对信号进展加强，从而提高整体或局部信号增益，使信号突出，到达突出异常，更好的判断异常区域的目的。数据处理的各个阶段都可以对探地雷达数据进展增益调节，这是因为增益调节并不会改变信号本身的特征量。因此，需要对回波信号恢复幅值，做相应的增强处理。这样不仅使探测的深度加深，又能从图像解释中获得更为详细的深度构造信息。331 振幅恢复振幅恢复，即恢复信号真实的幅值。目的是尽量恢复接收回波的真实信息，该方法是主要是利用介质的吸收系数与振幅值之间的关系。假设所探测的介质是均匀介质，则与发射天线距离为

38、Z的发射波振幅A为:3.5式中，Ao表示理想的振幅值即反射波真实振幅值;R表示该介质的吸收系数;1/Z表示发射波的扩散因子;t表示接收到该反射回波需要的时间。则由上式求得的回波信号的真实振幅值A。为:3.6根据回波时间与速度，可以计算得到回波所走的实际距离z=Vt V表示回波在该介质中传播的平均速度，则真实振幅值:3.7由以上分析，根据实际测得的信息，可得到探地雷达回波的振幅值。332 道振幅平衡采集到的回波信号，随着深度增加信号会有衰减现象，随着深度的增加，回波能量减弱。为了到达使探地雷达剖面信息相对清晰成像的效果，可以采用道平衡的方法。道平衡是将采样道振幅值大的按比例压缩，振幅值小的按比例

39、放大。这样各个道的振幅控制在一个合理围之。为了到达这个效果，可以将该采样道记录的振幅值在不同的阶段与不同权值相乘即可。设表示处理后的振幅值，表示原始回波振幅值，表示权系数，了表示该道第个采样点，则用公式表述为:3.8公式中相对于构成后的波形是具有相似性的。取值是随着时间和该振幅值的变化而变化的，与振幅的取值成反比。可以看出，是一个平均振幅值的比例的序列，并且变化相对缓慢。根据以上分析，设*通道信号采样点数为N，假设将其平均划分为K段，则每段长度为，即:3.9其中，符号表示采样间隔，符号表示取整。令表示分割后各段的振幅和其表达式为:3.10其中，m表示分割的第m段。表示每段的平均振幅，则则可表示

40、为加权系数，则均衡处理后的振幅为:3.11由上面公式可以得到：3.12由以上分析，最终得到道振幅均衡处理结果为:3.13式中c表示是个常数，代表道平均系数。由使用者确定，可以整体调整振幅值。333 道间振幅均衡道间均衡与道均衡的原理一致，不同之处在于，前者是*通道回波的信号的均衡，后者是道与道之间的均衡问题。道间均衡使各道的能量能够强弱均衡，保持在一定的围，为道号，为采样点，用公式表示为：3.14式中，表示道间均衡后的*道振幅值；表示未做处理前原始振幅值；表示道间均衡的加权系数。假设共有M个采样道，整个采样道的平均能量记为，针对各个采样道的平均能量记作，则3.15和的值，可以通过以下公式求得：

41、3.163.17道间均衡处理后的结果为3.18其中，S表示道间平衡系数，用于调节道间均衡后的振幅值。本研究中基于简单程度和实用性能及处理效果综合考虑选用的道增幅增益为自动时变增益算法。处理的结果图如下：图3.4 探地雷达时变增益wiggle图显示4 探地雷达成像算法的研究41 探地雷达扫描方式雷达成像的根本在于雷达和目标的相对转动，对探地雷达而言即探地雷达在地面多个方位向对目标进展照射，多方位的回波信号就包括了目标的空域信息。可采用基于单个收发天线对的合成孔径技术或者是基于阵列构造的实孔径技术进展地下目标探测和成像13。411 合成孔径技术下列图所示，为合成孔径扫描示意图。图4.1 合成孔径测

42、量合成孔径情况下，一个收发天线和一个接收天线组成一个收发天线对。这一天线对在空间不同位置处地下区域进展照射，从而完成合成孔径扫描14。在每一个孔径点处，发射天线发射电磁波，电磁波遇到目标体后产生外表感应电流，感应电流作为二次鼓励源向外辐射电磁波，回波信号仅被此孔径点处的接收天线接收，从而完成一次探测15。412 实孔径技术如下列图所示，为实孔径扫描示意图。图4.2 实孔径测量实孔径情况下，数据采集方式和地震波反射测量方式类似：多个接收天线形成一个阵列，多个发射天线沿同一条直线组成阵列构造。当发射天线阵元向探测区域发射电磁波时，地下目标散射场向全空间传播，接收天线阵列各阵元同时进展接收。为了进展

43、精细成像，需要尽可能多地记录目标的散射场信号。在实际应用中，为了简化操作，通常采用单个收发天线对的合成孔径技术。探地雷达系统多采用收发共用的天线。42 探地雷达波前成像算法波前重建最早作为地震波迁徙技术出现在地震学中。基于地震波与电磁波传播的相似性，Cafforio和Soumekh提出了两个根本一样的SAR模型。在系统模型中Cafforio和Soumekh运用一阶波恩近似，假定目标电磁散射可以模型化为目标各个离散散射中心相应的矢量和。设有一个稳态目标区含有多个散射点，各点的散射强度为，位于空域处。则目标特征函数可定义为 4.1雷达天线位于处照射地下目标区域，则天线与第n个目标之间的距离为。当天

44、线运动一段距离后，一个合成孔径就产生了，探地雷达合成孔径如下列图，表示成像区域中心到合成孔径的垂直距离。图4.3 探地雷达成像几何模型设雷达发射冲激脉冲信号为，天线在方位向累积角围增益不变，则回波信号可以表示为 (4.2)上式中表示媒质中的波速，近场区波传播因子合并于上式中的。对上式在时间域变量t作傅里叶变换，得(4.3)式中表示波数。上式中的域回波信号可以表示为目标的特征函数与球面相位函数在空间域的一维卷积，即 (4.4)(4.5)式子中的定义为u域的卷积。相位函数为SAR系统在给定距离是u域的转移函数。SAR成像就是对在u域的解卷积以去除转移函数的影响。这一逆运算过程就是一种目标重建过程，

45、在此称为波前重建8-11。对，在空间u域进展傅里叶变换，得到(4.6)而目标特征函数的二维傅里叶变换得到(4.7)因此可以得出(4.8)为了运用快速傅里叶变换实现算法中的目标重建，需要得到回波信号在平面的等间隔分布。将极坐标中等间隔分布的数据转化为直角坐标中等间隔分布数据的处理过程称为Stolt插值，插值运算的根本关系可表示为(4.9)常用的插值算法有最邻近点近似、临近四点近似、线性插值、Cubic插值和三次样条插值法。这类成像方法中，傅里叶变换可以用到快速傅里叶变换实现，运算量最要集中在二维非线性插值12。波前成像处理过程如下：图4.4 波前成像算法流程图探地雷达波前成像重建可以在谱域实现，

46、也可以在时域实现。在重建质量上，二者没有区别。但是在工程实现上各有利弊。谱域实现主要涉及傅里叶变换和插值，时域主要涉及求和和插值。从算法实现上，由于快速傅里叶变换的采用和目前VLSI、VHSI技术水平的不断提高，谱域的快速实现成为了可能。对时域实现，逐点计算想干点位置仍然无快速方法，但是可以在成像中考虑波传播扩散和衰减的影响，并且在做补偿处理时跟容易。图4.5 波前成像算法成像结果42 探地雷达投影层析成像高频情况下，波传播的衍射效应可以忽略，接收到散射场可以表示为沿波阵面上目标各强散射点回波值的叠加。此时目标可以视为多个强散射中心的组合，目标反演就变成了各散射中心位置和散射强度的反演11。与

47、CT类似，通过接收空间个方位的散射场就可以通过层析技术实现目标的反演9。目标函数可以用逆拉顿变换得到。建立空域和谱域几何关系如下列图。ab图4.4 实际空间和傅里叶空间的目标用表示基于散射中心模型的目标散射强度函数，则单次雷达照射是产生(4.10)其中表示扫描坐标系。表示沿直线的散射点的积分值，记目标散射强度函数的空域二维傅里叶变换为(4.11)从而有(4.12)经过一次投影后，可以由上式中的的一维傅里叶变换导出的二维傅里叶变换的一个切片，即(4.13)这就是投影切片定理。为了得到目标谱域空间的另一切片，只需要在另一个角度下重复进展一次照射。当从到对目标进展照射时，得到的散射场回波的傅里叶变换

48、就可以填充整个傅里叶空间，从可以通过傅里叶反变换重构目标散射强度函数11。不同角度照射下目标散射场数据的谱域空间表示。图4.6 投影层析成像算法实现结果5 结论本文针对探地雷达探测原理及一些信号处理和成像方法，利用MATLAB强大的矩阵运算能力，进展基于MATLAB的探地雷达成像算法研究。在次研究中，本文运用了奇异值分解法去除探地雷达探测信号中的直达波成分，运用一种改良的中值滤波方法去除探地雷达探测信号中的背景噪声，运用的道增幅增益来减少电磁波在传播中的衰减。在成像算法中，基于合成孔径扫描方式的成像算法研究为重点，并从众多的成像算法中选取其中两种用MATLAB实现。由于本人知识水平和个人能力有

49、限，探地雷达成像成像的效果不尽人意。如在数据采集的过程中由于很多原因会产生数据坏道，在本研究中没有去出去坏道的数据。在自动时变增益中，最好的选择方法为道间增幅增益，但是经研究发现道间增幅增益实现算法过于复杂因此选用道增幅增益。所以增益的效果不是很好。在成像算法实现中，由于前面数据的预处理局部做的不够充分，导致成像效果还存在着一些瑕疵。附录clear all;close all;%探地雷达数据读入RadarFile,RadarPath = uigetfile(*.LTE,选择原始的雷达数据);fid = fopen(RadarPath RadarFile,r);if fid = -1 disp

50、(选择数据错误，请从新选择); return;else endyuanshishuju = fread(fid,inf,unsigned char);fseek(fid,4,-1);%获得采样点数samples = fread(fid,1,unsigned short);disp(采样点数=);disp(samples);%数据位数num_bit = fread(fid,1,unsigned short);disp(数据位数=);shujuzijie = num_bit/8;disp(num_bit);%数据零偏num_zero = fread(fid,1,unsigned short);di

51、sp(数据零偏=);disp(num_zero);%每秒扫描道数daoshu = fread(fid,1,float);disp(每秒扫描道数=);disp(daoshu);%测量文件总的数据量frewind(fid);fseek(fid,0,1);sum_num = ftell(fid);if sum_num = -1 fseek(fid,0,1);sum_num = ftell(fid);else disp(总数据量=);disp(sum_num);end%计算总道数c*shujuliang = (sum_num - 1024)/shujuzijie; %1024 是数据头局部num_tr

52、aces = c*shujuliang/samples;disp(数据总道数=);disp(num_traces);%过滤数据头局部frewind(fid);fseek(fid,1024,-1);%读取每道数据for j = 1:num_traces A(:,j) = fread(fid,samples,bit16);endfclose(fid);%显示原始探地雷达图像%自定义色图模式mdsj = 1/daoshu;mdiansj = mdsj/samples;fs = 1/mdiansj; %采样频率t = 0:1/fs:1;rk(1:32,1)=(0:1/31:1);gk(1:32,1)=

53、(0:1/31:1);bk(1:32,1)=(0.4:(1-0.4)/31:1);rk(33:64,1)=(1:-(1-0.4)/31:0.4);gk(33:64,1)=(1:-1/31:0);bk(33:64,1)=(1:-1/31:0);map = rk gk bk;figure(1);%imshow(A);imagesc(A);colormap(map);title(探地雷达原始数据图像彩色显示);figure(2);imagesc(A);colormap(gray);title(探地雷达原始数据图像灰度显示);figure(3); * = 1:1:size(A,2); t = 1:1:

54、size(A,1); showma*=100; wscale = 1; dma*=ma*(nanma*(abs(A); d* = *(2) - *(1); ntraces=length(*); d=ntraces/showma*; if d ntraces, break end *t=A(:,i)./dma*; plot(*(i) + *t, t ,k-,linewidth,.1); if i=1 hold on; end enda*is(min(*) ma*(*) 0 ma*(t);title(探地雷达wiggle图显示);figure(4);plot(A(:,1000);%去除直达波，采用

55、主元分析去除法%先对原始矩阵A进展奇异值分解U,D,V = svd(A);%对特征值矩阵按从大到小排序p = zeros(1,samples);for k = 1:samples p(k) = D(k,k);endP = sort(p,2,descend);for k = 1:samples D(k,k) = P(k);endV1 = V ;Uni = inv(U);Y = Uni*A;y1 = D(1,1)*V1(1,:);A1 = U(:,1)*Y(1,:);CZDBA = A - A1;figure(8);imagesc(CZDBA);colormap(map);title(探地雷达去直

56、达波数据图像);jz = zeros(samples,num_traces);ffzzlb1 = zeros(samples,num_traces);for k = 1:samples j = 1; B = 0; while j= num_traces B = B + CZDBA(k,j); j = j + 1; end b = 1; while b = num_traces jz(k,b) = B / num_traces; b = b + 1; endendjzlb = CZDBA - jz;ffjzlb = zeros(samples,num_traces);for j = 1:num_

57、traces ffjzlb(:,j) = fft(jzlb(:,j);endfuzhipu = abs(ffjzlb);*iangpinpu = angle(ffjzlb);ffzzlb = medfilt2(fuzhipu);for j = 1:num_traces k = 1; while k =samples ffzzlb1(k,j) = ffzzlb(k,j)*cos(*iangpinpu(k,j)+i*ffzzlb(k,j)*sin(*iangpinpu(k,j); k = k + 1; endendzzlb = zeros(samples,num_traces);for j =1:

58、num_traces zzlb(:,j) = ifft(ffzzlb1(:,j);endK1 = abs(zzlb); figure(9);imagesc(K1);colormap(map);title(探地雷达经滤波器后数据图像);%自动时变增益aver = zeros(samples,num_traces);K2 = zeros(samples,num_traces);N = 5; %时窗大小M = 1;sum1 = 0;for k = 1:num_traces M = 1; while M samples for j = 0:1:N-1 if j+M = samples sum1 = s

59、um1 + abs(K1(M+j,k); end end if sum1 = 0 sum1 = 1; end for j = 0:1:N-1 if j+M = samples aver(M+j,k) = sum1; end end M = M + N; endend*ishu = 1000;j = 1;for k = 1:num_traces j = 1; while j samples K2(j,k) = *ishu*K1(j,k)*(2*N+1)/aver(j,k); j = j + 1; endendfigure(3); * = 1:1:size(K2,2); t = 1:1:size(

60、K2,1); showma*=100; wscale = 1; dma*=ma*(nanma*(abs(K2); d* = *(2) - *(1); ntraces=length(*); d=ntraces/showma*; if d ntraces, break end *t=K2(:,i)./dma*;% if strcmp(style,vararea)=1,% *t1=*t;% *t1(find(*t1 0) = 0;% % fill(*(1) + (i-1)*d*) + *t,fliplr(*t1),t,fliplr(t),0 0 0);% fill(*(i) + *t,fliplr(