基于自適應(yīng)稀疏分解的探地雷達(dá)水平分層介質(zhì)時(shí)延估計(jì)

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

（南昌大學(xué)，江西 南昌　330031）

基于自適應(yīng)稀疏分解的探地雷達(dá)水平分層介質(zhì)時(shí)延估計(jì)

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

（南昌大學(xué)，江西 南昌330031）

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

自適應(yīng)稀疏分解；探地雷達(dá)；水平分層介質(zhì)；過完備字典；時(shí)延估計(jì)

0　引　言

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

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

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

1　基于自適應(yīng)稀疏分解的探地雷達(dá)時(shí)延估計(jì)

一般來說，雷達(dá)發(fā)射信號(hào)是先驗(yàn)的，為了簡(jiǎn)化分析，假定時(shí)延后的發(fā)射信號(hào)s（t）經(jīng)幅度擴(kuò)展或壓縮構(gòu)成目標(biāo)回波信號(hào)。若雷達(dá)探測(cè)區(qū)域內(nèi)有K個(gè)目標(biāo)，則雷達(dá)接收的回波信號(hào)可表示為：

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

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

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

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

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

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

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

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

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

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

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

當(dāng)?shù)螖?shù)k超過目標(biāo)空間稀疏度K時(shí)，終止迭代過程。當(dāng)?shù)螖?shù)小于K時(shí)，無法得到每個(gè)目標(biāo)的時(shí)延；當(dāng)?shù)螖?shù)大于K時(shí)，雖然重構(gòu)精度會(huì)得到提高，但對(duì)時(shí)延估計(jì)精度并無影響。自適應(yīng)稀疏分解算法流程如下：

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

3.重復(fù)：

5.Φk←φj*，構(gòu)成子字典

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

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

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

9.?k←k+1

10.若k≥K則停止迭代，否則轉(zhuǎn)步驟4。

2　實(shí)驗(yàn)結(jié)果

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

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

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

圖1　BΔτ=1時(shí)ASD估計(jì)結(jié)果

為了對(duì)MUSIC方法和ASD方法進(jìn)行時(shí)延估計(jì)精度對(duì)比，MUSIC估計(jì)結(jié)果如圖2所示。兩種方法皆能精確計(jì)算探地雷達(dá)信號(hào)時(shí)延，但ASD精度更高。為了驗(yàn)證本文方法在時(shí)延估計(jì)方面的有效性，表1列出了MUSIC和ASD針對(duì)不同處理時(shí)間分辨的BΔτ估計(jì)結(jié)果，圖3為MUSIC和ASD在不同處理時(shí)間分辨下的時(shí)延估計(jì)相對(duì)誤差。從表1和圖3可以看出，兩種方法的估計(jì)精度皆與理論BΔτ值成正比，并且ASD估計(jì)可靠性高于MUSIC。

表1　不同BΔτ條件下的MUSIC，ASD時(shí)延估計(jì)結(jié)果

圖2　BΔτ=1時(shí)MUSIC估計(jì)結(jié)果

圖3　不同BΔτ下的時(shí)延估計(jì)誤差對(duì)比

表2和圖4列出了當(dāng)BΔτ=1時(shí)，不同信噪比條件下的MUSIC，ASD時(shí)延估計(jì)結(jié)果。兩種算法估計(jì)精度與SNR成正比，當(dāng)信噪比過低時(shí)MUSIC無法估計(jì)時(shí)延，而ASD仍能粗略估計(jì)時(shí)延。當(dāng)信噪比提高，MUSIC和ASD均能粗略估計(jì)時(shí)延，并且ASD精度更高。當(dāng)信噪比高于一定值時(shí)，兩種方法的估計(jì)精度都趨于飽和。

表2　不同信噪比條件下的M USIC，ASD時(shí)延估計(jì)結(jié)果

2.2實(shí)測(cè)數(shù)據(jù)

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

圖4　不同SNR下的時(shí)延估計(jì)誤差對(duì)比

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

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

3　結(jié)　語

本文提出一種基于自適應(yīng)稀疏分解的探地雷達(dá)水平分層介質(zhì)時(shí)延估計(jì)算法，假設(shè)場(chǎng)景是一個(gè)介電常數(shù)固定的空?qǐng)鼍?，以合適的分辨率劃分場(chǎng)景，則總能找到空?qǐng)鼍爸信c實(shí)際目標(biāo)時(shí)延相對(duì)應(yīng)的某一點(diǎn)，保證了算法能夠自適應(yīng)選取與目標(biāo)對(duì)應(yīng)的時(shí)延原子。該算法在時(shí)域中進(jìn)行，與MUSIC相比，計(jì)算量小，估計(jì)精度高，提高了雷達(dá)的處理時(shí)間分辨率，并且不需已知層介質(zhì)的介電常數(shù)、厚度等，適用于工程探測(cè)。仿真和實(shí)驗(yàn)結(jié)果證明了該算法的優(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—），女，福建寧德人，碩士。研究方向?yàn)橥ㄐ偶夹g(shù)。