陈青玉，周辉林，王玉皞，段荣行

（南昌大学，江西 南昌　330031）

基于自适应稀疏分解的探地雷达水平分层介质时延估计

陈青玉，周辉林，王玉皞，段荣行

（南昌大学，江西 南昌330031）

传统的时延估计方法受瑞利限限制，一些基于子空间的方法可以达到高分辨或超分辨，然而，子空间方法需要处理的数据量大，计算复杂度高。提出一种新的基于自适应稀疏分解的探地雷达水平分层介质时延估计方法，自适应稀疏分解将传统的参数估计问题转化为字典学习问题。自适应稀疏分解从过完备字典中选取少量原子，其对应的时间参数即为所需估计的时延。与传统的基于子空间方法相比，该算法在低信噪比时具有更高的估计正确率，并且直接在时域进行，减少了计算复杂度。仿真及实测结果表明，该算法与MUSIC算法对比在时延估计方面具有更大的优势。

自适应稀疏分解；探地雷达；水平分层介质；过完备字典；时延估计

0　引　言

探地雷达（Ground Penetrating Radar，GPR）又称为地质雷达，是一种利用超高频（106～109Hz）脉冲电磁波进行地下浅层目标（诸如管道、电缆、地雷以及路面隐患）探测的设备，具有探测速度快、探测过程连续、操作方便灵活、分辨率高、不损坏被探测目标等特点［1］，广泛应用于考古［2］、矿产资源勘探［3］、路面分析［4］和地雷探测［5］等众多领域。

探地雷达回波的幅度和时延包含目标位置及电磁散射特性等信息，因此为了对后续数据进行有效处理以及对目标精准定位，需要对回波幅度和时延进行精确估计。传统的时延估计技术包括基于观测矢量的逆傅里叶变换和互相关函数方法［6］。这两种方法容易实现，但分辨率受瑞利限限制，并且，当存在多个回波时，估计精度下降。文献［7-9］提出多种高分辨率正弦信号参数估计方法，如多重信号分类（Multiple Signal Classification，MUSIC）、线性预测、最大似然法等，然而这些方法必须应用于复幅度信号，并且计算量大。

本文假设发射信号已知或已测量得到结果，提出一种基于自适应稀疏分解（Adaptive Sparse Decomposition，ASD）的时延估计方法，能够大大提高时延估计精度，最小化地下介质衰减属性和雷达系统参数的影响。冗余字典的构造是稀疏分解过程最关键的步骤，常用的方法是基于模型和基于学习的字典［10］。传统的字典，如Fourier字典和小波字典构造简单，能够较好地处理一维信号，但无法处理更加复杂的或高维的信号。Wenbin Shao提出一种基于稀疏分解过程中再定义的自适应过完备Gabor字典的稀疏表示方法［11］，但并未将该字典应用于时延估计。本文提出的过完备字典由场景中不同位置的反射信号组成，该方法不仅可应用于稀疏重构，还可应用于时延估计。

1　基于自适应稀疏分解的探地雷达时延估计

一般来说，雷达发射信号是先验的，为了简化分析，假定时延后的发射信号s（t）经幅度扩展或压缩构成目标回波信号。若雷达探测区域内有K个目标，则雷达接收的回波信号可表示为：

式中：K为回波个数；αq为目标q的反射系数；τq为发射天线经目标q反射到达接收天线的时延；n（t）为独立同分布的高斯白噪声。

冗余字典的类型有多种，本文采用的是匹配字典，其根据与信号最匹配原则构造原子，减少了搜索匹配原子过程的计算量。首先，将成像区域离散化，对于二维成像区域，在方位向和距离向将其分为NT=×个像素，这样便可建立一个目标空间位置矢量B=｛π1，π2，…，πNT｝，其中πj是一个关于目标位置的二维矢量。定义一个指示目标空间属性的矢量α=［α1，α2，…，αNT］T，如果网格πj位置存在目标，则相应αj为非零元素，否则为0。

忽略天线直达波和杂波时，根据式（1），第i个收发位置接收天线的回波信号与指示函数的关系可用一个矩阵Ψi联系起来：

式中τi（πj）为第i个收发天线的双程时延到第 j个网格的时延。由式（3）可知，在不需要已知幅度因子的条件下，只需计算成像场景中各网格点到第i个收发天线的双程时延即可构造GPR位于第i个收发位置时相对各目标位置的字典。

假设探地雷达信号具有N个时间采样点，探测区域离散化成NT个网格点。假设探测区域是一个介电常数固定的空场景，则可构造时延字典（Ψi=［φ1，φ2，…，φNT］），稀疏系数αj在原子的迭代选择过程中优化。由以上分析可知，字典中的每个原子是发射波形经幅度伸缩和时间延迟的信号，保证了字典对探地雷达信号自适应。

一般来说背景和目标的介电常数差值有限，因此只要保证成像区域足够大，则目标的时延总能对应于场景内某一点的时延，这保证了算法能够选择与目标反射回波的最佳匹配时延原子。然后寻找第k次迭代过程中与残差信号互相关值最大的前k个原子φj*：

式中：rkj（τ）为时延原子φj（t）与残差信号s~k-1的互相关函数。这k个原子组成一个子字典Φk。通过下式更新系数矢量：

Φk由k次迭代中的子字典组成。最后，更新残差：

可增大自适应稀疏分解算法初始迭代阈值提高重构精度，但时延估计精度不受迭代次数影响，每次迭代过程的子字典Φk中的第一个原子的波峰或波谷对应的时间参数即为各目标的时延。对应给定的带宽B和两回波之间的最小时延Δτ可定义处理时间分辨率BΔτ描述时间分辨力。

为了确定各目标对应的时延原子，设置一个多层水平介质场景，当电磁波遇到存在电性差异的分界面时发生反射，反射系数可表示为：

式中：r为反射系数；εl为第l层介质的相对介电常数；εl+1为第l+1层介质的相对介电常数。从式（7）可以看出，反射系数幅度依赖于上下层介质的电性差，且反射信号的幅度正比于反射系数。因此，迭代过程中子字典Φk中的第一个原子对应的时间参数对应于电性差异值由大到小的水平分界面反射时延。

当迭代次数k超过目标空间稀疏度K时，终止迭代过程。当迭代次数小于K时，无法得到每个目标的时延；当迭代次数大于K时，虽然重构精度会得到提高，但对时延估计精度并无影响。自适应稀疏分解算法流程如下：

2.计算Ψi，即第i道数据的时延字典

3.重复：

5.Φk←φj*，构成子字典

6.Φk=［Φk-1，Φk］，填充矩阵

7. αk=Φk）-1s，更新系数矢量

8. s~k=s-Φkαk，更新残差

9. k←k+1

10.若k≥K则停止迭代，否则转步骤4。

2　实验结果

2.1仿真数据

使用GprMax建立水平层状介质探测模型，发射雷克子波中心频率为1 GHz，层介质厚度以步进间隔0.01 m 从0.04～0.08 m不断增加，介电常数设为6.0，天线距离层介质上表面0.01 m。

本实验雷达回波包含直达波、上下表面反射波以及噪声，首先移除杂波，则根据ASD原理，当迭代次数等于目标数时即可估计每个目标对应的时延。本实验中，目标反射即为介质层上下表面的反射，因此当迭代次数等于2时，子字典Φ1和Φ2的第一个原子分别对应上下表面的时延，由此得到电磁波在介质层内部经历的时延。当d=0.06m，即BΔτ=1时，时延估计结果如图1所示。目标在介质层经历的时延为1.007 9 ns，十分接近理论值0.983 1 ns，该算法有效。该结果还验证了字典原子皆是经幅度扩展和时间延迟后的发射波。

图1　BΔτ=1时ASD估计结果

为了对MUSIC方法和ASD方法进行时延估计精度对比，MUSIC估计结果如图2所示。两种方法皆能精确计算探地雷达信号时延，但ASD精度更高。为了验证本文方法在时延估计方面的有效性，表1列出了MUSIC和ASD针对不同处理时间分辨的BΔτ估计结果，图3为MUSIC和ASD在不同处理时间分辨下的时延估计相对误差。从表1和图3可以看出，两种方法的估计精度皆与理论BΔτ值成正比，并且ASD估计可靠性高于MUSIC。

表1　不同BΔτ条件下的MUSIC，ASD时延估计结果

图2　BΔτ=1时MUSIC估计结果

图3　不同BΔτ下的时延估计误差对比

表2和图4列出了当BΔτ=1时，不同信噪比条件下的MUSIC，ASD时延估计结果。两种算法估计精度与SNR成正比，当信噪比过低时MUSIC无法估计时延，而ASD仍能粗略估计时延。当信噪比提高，MUSIC和ASD均能粗略估计时延，并且ASD精度更高。当信噪比高于一定值时，两种方法的估计精度都趋于饱和。

表2　不同信噪比条件下的M USIC，ASD时延估计结果

2.2实测数据

本实验数据来源于2012年12月17日第三方厚度检测单位对吉安至莲花项目的检测报告。收发天线间隔为0.06 m，中心频率为1.6 GHz。天线每隔18.18 cm采集一道数据，本实验采集300道数据进行时延估计。截取时间窗大小为6.428 9 ns，采样间隔为0.032 31 ns，利用GPR数据处理软件“Reflexw”读取原始数据，结果如图5（a）所示。图中，第一道波为直达波，出现在1 ns左右，实际处理中需要抑制直达波。第二、三道波分别为上下层界面的反射波，从图中可读取出电磁波在此面层中的速度为0.122 47 m/ns，由此可得此面层相对介电常数为6.0，此面层厚度大约为5 cm。

图4　不同SNR下的时延估计误差对比

为了得到雷达发射波形，可以在雷达正下方放置一块金属板，由此得到发射波形的全反射波形，进而得到发射波形。选取每一道数据进行时延估计，MUSIC 和ASD的时延反演结果如图5（b）所示，对比图5（a）和图5（b），ASD反演结果更加准确，MUSIC抑制扰动能力相对ASD较差。

图5　原始数据与时延反演结果

3　结　语

本文提出一种基于自适应稀疏分解的探地雷达水平分层介质时延估计算法，假设场景是一个介电常数固定的空场景，以合适的分辨率划分场景，则总能找到空场景中与实际目标时延相对应的某一点，保证了算法能够自适应选取与目标对应的时延原子。该算法在时域中进行，与MUSIC相比，计算量小，估计精度高，提高了雷达的处理时间分辨率，并且不需已知层介质的介电常数、厚度等，适用于工程探测。仿真和实验结果证明了该算法的优越性。

［1］ZHANG C G，SONG M Y，WANG L G.Application of GPR in hole survey of the subway project in Beijing，China［J］.Energy procedia，2012，16：702-706.

［2］ZHAOW，FORTE E，PIPAN M，et al.Ground penetrating radar（GPR）attribute analysis for archaeological prospection［J］. Journal of applied geophysics，2013，97（10）：107-117.

［3］FRANCKE J.Applications of GPR in mineral resource evaluations［C］//Proceedings of 2010 13th International Conference on Ground Penetrating Radar（GPR）.Lecce：IEEE，2010：1-5.

［4］DENIS A，HUNEAU F，HERL S，et al.GPR data processing for fractures and flakes detection in sandstone［J］.Journal of applied geophysics，2009，68（2）：282-288.

［5］YAROVOY A.Landmine and unexploded ordnance detection and classification with ground penetrating radar［M］//Anon. Ground penetrating radar theory and applications.Amsterdam：Elsevier，2008：445-478.

［6］SAKHNOV K，VERTELETSKAYA E，SIMAK B.Echo delay estimation using algorithms based on cross-correlation［J］.Journal of convergence information technology，2011，6（4）：1-11.

［7］SHRESTHA S M，ARAI I.Signal processing of ground penetrating radar using spectral estimation techniques to estimate the position of buried targets［J］.EURASIP journal on advances in signal processing，2003（12）：1-12.

［8］BASTARD C L，BALTAZART V，WANG Y，et al.Thin-pavement thickness estimation using GPR with high-resolution and superresolution methods［J］.IEEE transactions on geoscience and remote sensing，2007，45（8）：2511-2519.

［9］LI J，WU R.An efficient algorithm for time delay estimation［J］. IEEE transactions on signal processing，1998，46（8）：2231-2235.

［10］RUBINSTEIN R，BRUCKSTEIN A M，ELAD M.Dictionaries for sparse representation modeling［J］.Proceedings of the IEEE，2010，98（6）：1045-1057.

［11］SHAOW，BOUZERDOUM A，PHUNG S L.Sparse representation of GPR traces with application to signal classification［J］. IEEE transactions on geoscience and remote sensing，2013，51（7）：3922-3930.

GPR horizontally-layered m edium tim e delay estim ation based on adap tive sparse decom position

CHEN Qingyu，ZHOU Huilin，WANG Yuhao，DUAN Rongxing
（Nanchang University，Nanchang 330031，China）

Since the traditional time delay estimationmethod is restricted by the Rayleigh limit，some subspace-basedmethods can reach high resolution or super-resolution.However，the subspace-based methods have large data size for processing and high computation complexity.A new horizontally-layered medium time delay estimation method of using ground penetrating radar （GPR）is proposed，which is based on adaptive sparse decomposition（ASD），and in which the traditional parameter estimation is turned to the dictionary learning by ASD.A few atoms are selected from the over-complete dictionary by means of ASD，and their corresponding time parameter is the estimating time delay.In comparison with the traditional subspace-based methods，this algorithm has higher estimation accuracy at low signal-to-noise ratio（SNR），and can directly run in time domain，so as to reduce computational complexity.The simulation and actualmeasured results show that the proposed algorithm has better advantage than the MUSIC（multiple signal classification）algorithm in the aspect of time delay estimation.

adaptive sparse decomposition；ground penetrating radar；horizontally-layered medium；over-complete dictionary；time delay estimation

TN95-34

A

1004-373X（2016）11-0009-04

10.16652/j.issn.1004-373x.2016.11.003

2015-09-02

陈青玉（1990—），女，福建宁德人，硕士。研究方向为通信技术。