陳青玉，周輝林，王玉皞，段榮行

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

基于自適應稀疏分解的探地雷達水平分層介質時延估計

陳青玉，周輝林，王玉皞，段榮行

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

傳統(tǒng)的時延估計方法受瑞利限限制，一些基于子空間的方法可以達到高分辨或超分辨，然而，子空間方法需要處理的數(shù)據(jù)量大，計算復雜度高。提出一種新的基于自適應稀疏分解的探地雷達水平分層介質時延估計方法，自適應稀疏分解將傳統(tǒng)的參數(shù)估計問題轉化為字典學習問題。自適應稀疏分解從過完備字典中選取少量原子，其對應的時間參數(shù)即為所需估計的時延。與傳統(tǒng)的基于子空間方法相比，該算法在低信噪比時具有更高的估計正確率，并且直接在時域進行，減少了計算復雜度。仿真及實測結果表明，該算法與MUSIC算法對比在時延估計方面具有更大的優(yōu)勢。

自適應稀疏分解；探地雷達；水平分層介質；過完備字典；時延估計

0　引　言

探地雷達（Ground Penetrating Radar，GPR）又稱為地質雷達，是一種利用超高頻（106～109Hz）脈沖電磁波進行地下淺層目標（諸如管道、電纜、地雷以及路面隱患）探測的設備，具有探測速度快、探測過程連續(xù)、操作方便靈活、分辨率高、不損壞被探測目標等特點［1］，廣泛應用于考古［2］、礦產(chǎn)資源勘探［3］、路面分析［4］和地雷探測［5］等眾多領域。

探地雷達回波的幅度和時延包含目標位置及電磁散射特性等信息，因此為了對后續(xù)數(shù)據(jù)進行有效處理以及對目標精準定位，需要對回波幅度和時延進行精確估計。傳統(tǒng)的時延估計技術包括基于觀測矢量的逆傅里葉變換和互相關函數(shù)方法［6］。這兩種方法容易實現(xiàn)，但分辨率受瑞利限限制，并且，當存在多個回波時，估計精度下降。文獻［7-9］提出多種高分辨率正弦信號參數(shù)估計方法，如多重信號分類（Multiple Signal Classification，MUSIC）、線性預測、最大似然法等，然而這些方法必須應用于復幅度信號，并且計算量大。

本文假設發(fā)射信號已知或已測量得到結果，提出一種基于自適應稀疏分解（Adaptive Sparse Decomposition，ASD）的時延估計方法，能夠大大提高時延估計精度，最小化地下介質衰減屬性和雷達系統(tǒng)參數(shù)的影響。冗余字典的構造是稀疏分解過程最關鍵的步驟，常用的方法是基于模型和基于學習的字典［10］。傳統(tǒng)的字典，如Fourier字典和小波字典構造簡單，能夠較好地處理一維信號，但無法處理更加復雜的或高維的信號。Wenbin Shao提出一種基于稀疏分解過程中再定義的自適應過完備Gabor字典的稀疏表示方法［11］，但并未將該字典應用于時延估計。本文提出的過完備字典由場景中不同位置的反射信號組成，該方法不僅可應用于稀疏重構，還可應用于時延估計。

1　基于自適應稀疏分解的探地雷達時延估計

一般來說，雷達發(fā)射信號是先驗的，為了簡化分析，假定時延后的發(fā)射信號s（t）經(jīng)幅度擴展或壓縮構成目標回波信號。若雷達探測區(qū)域內有K個目標，則雷達接收的回波信號可表示為：

式中：K為回波個數(shù)；αq為目標q的反射系數(shù)；τq為發(fā)射天線經(jīng)目標q反射到達接收天線的時延；n（t）為獨立同分布的高斯白噪聲。

冗余字典的類型有多種，本文采用的是匹配字典，其根據(jù)與信號最匹配原則構造原子，減少了搜索匹配原子過程的計算量。首先，將成像區(qū)域離散化，對于二維成像區(qū)域，在方位向和距離向將其分為NT=×個像素，這樣便可建立一個目標空間位置矢量B=｛π1，π2，…，πNT｝，其中πj是一個關于目標位置的二維矢量。定義一個指示目標空間屬性的矢量α=［α1，α2，…，αNT］T，如果網(wǎng)格πj位置存在目標，則相應αj為非零元素，否則為0。

忽略天線直達波和雜波時，根據(jù)式（1），第i個收發(fā)位置接收天線的回波信號與指示函數(shù)的關系可用一個矩陣Ψi聯(lián)系起來：

式中τi（πj）為第i個收發(fā)天線的雙程時延到第 j個網(wǎng)格的時延。由式（3）可知，在不需要已知幅度因子的條件下，只需計算成像場景中各網(wǎng)格點到第i個收發(fā)天線的雙程時延即可構造GPR位于第i個收發(fā)位置時相對各目標位置的字典。

假設探地雷達信號具有N個時間采樣點，探測區(qū)域離散化成NT個網(wǎng)格點。假設探測區(qū)域是一個介電常數(shù)固定的空場景，則可構造時延字典（Ψi=［φ1，φ2，…，φNT］），稀疏系數(shù)αj在原子的迭代選擇過程中優(yōu)化。由以上分析可知，字典中的每個原子是發(fā)射波形經(jīng)幅度伸縮和時間延遲的信號，保證了字典對探地雷達信號自適應。

一般來說背景和目標的介電常數(shù)差值有限，因此只要保證成像區(qū)域足夠大，則目標的時延總能對應于場景內某一點的時延，這保證了算法能夠選擇與目標反射回波的最佳匹配時延原子。然后尋找第k次迭代過程中與殘差信號互相關值最大的前k個原子φj*：

式中：rkj（τ）為時延原子φj（t）與殘差信號s~k-1的互相關函數(shù)。這k個原子組成一個子字典Φk。通過下式更新系數(shù)矢量：

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

可增大自適應稀疏分解算法初始迭代閾值提高重構精度，但時延估計精度不受迭代次數(shù)影響，每次迭代過程的子字典Φk中的第一個原子的波峰或波谷對應的時間參數(shù)即為各目標的時延。對應給定的帶寬B和兩回波之間的最小時延Δτ可定義處理時間分辨率BΔτ描述時間分辨力。

為了確定各目標對應的時延原子，設置一個多層水平介質場景，當電磁波遇到存在電性差異的分界面時發(fā)生反射，反射系數(shù)可表示為：

式中：r為反射系數(shù)；εl為第l層介質的相對介電常數(shù)；εl+1為第l+1層介質的相對介電常數(shù)。從式（7）可以看出，反射系數(shù)幅度依賴于上下層介質的電性差，且反射信號的幅度正比于反射系數(shù)。因此，迭代過程中子字典Φk中的第一個原子對應的時間參數(shù)對應于電性差異值由大到小的水平分界面反射時延。

當?shù)螖?shù)k超過目標空間稀疏度K時，終止迭代過程。當?shù)螖?shù)小于K時，無法得到每個目標的時延；當?shù)螖?shù)大于K時，雖然重構精度會得到提高，但對時延估計精度并無影響。自適應稀疏分解算法流程如下：

2.計算Ψi，即第i道數(shù)據(jù)的時延字典

3.重復：

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

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

7.?αk=Φk）-1s，更新系數(shù)矢量

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

9.?k←k+1

10.若k≥K則停止迭代，否則轉步驟4。

2　實驗結果

2.1仿真數(shù)據(jù)

使用GprMax建立水平層狀介質探測模型，發(fā)射雷克子波中心頻率為1 GHz，層介質厚度以步進間隔0.01 m 從0.04～0.08 m不斷增加，介電常數(shù)設為6.0，天線距離層介質上表面0.01 m。

本實驗雷達回波包含直達波、上下表面反射波以及噪聲，首先移除雜波，則根據(jù)ASD原理，當?shù)螖?shù)等于目標數(shù)時即可估計每個目標對應的時延。本實驗中，目標反射即為介質層上下表面的反射，因此當?shù)螖?shù)等于2時，子字典Φ1和Φ2的第一個原子分別對應上下表面的時延，由此得到電磁波在介質層內部經(jīng)歷的時延。當d=0.06m，即BΔτ=1時，時延估計結果如圖1所示。目標在介質層經(jīng)歷的時延為1.007 9 ns，十分接近理論值0.983 1 ns，該算法有效。該結果還驗證了字典原子皆是經(jīng)幅度擴展和時間延遲后的發(fā)射波。

圖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實測數(shù)據(jù)

本實驗數(shù)據(jù)來源于2012年12月17日第三方厚度檢測單位對吉安至蓮花項目的檢測報告。收發(fā)天線間隔為0.06 m，中心頻率為1.6 GHz。天線每隔18.18 cm采集一道數(shù)據(jù)，本實驗采集300道數(shù)據(jù)進行時延估計。截取時間窗大小為6.428 9 ns，采樣間隔為0.032 31 ns，利用GPR數(shù)據(jù)處理軟件“Reflexw”讀取原始數(shù)據(jù)，結果如圖5（a）所示。圖中，第一道波為直達波，出現(xiàn)在1 ns左右，實際處理中需要抑制直達波。第二、三道波分別為上下層界面的反射波，從圖中可讀取出電磁波在此面層中的速度為0.122 47 m/ns，由此可得此面層相對介電常數(shù)為6.0，此面層厚度大約為5 cm。

圖4　不同SNR下的時延估計誤差對比

為了得到雷達發(fā)射波形，可以在雷達正下方放置一塊金屬板，由此得到發(fā)射波形的全反射波形，進而得到發(fā)射波形。選取每一道數(shù)據(jù)進行時延估計，MUSIC 和ASD的時延反演結果如圖5（b）所示，對比圖5（a）和圖5（b），ASD反演結果更加準確，MUSIC抑制擾動能力相對ASD較差。

圖5　原始數(shù)據(jù)與時延反演結果

3　結　語

本文提出一種基于自適應稀疏分解的探地雷達水平分層介質時延估計算法，假設場景是一個介電常數(shù)固定的空場景，以合適的分辨率劃分場景，則總能找到空場景中與實際目標時延相對應的某一點，保證了算法能夠自適應選取與目標對應的時延原子。該算法在時域中進行，與MUSIC相比，計算量小，估計精度高，提高了雷達的處理時間分辨率，并且不需已知層介質的介電常數(shù)、厚度等，適用于工程探測。仿真和實驗結果證明了該算法的優(yōu)越性。

［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，F(xiàn)ORTE 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—），女，福建寧德人，碩士。研究方向為通信技術。