自適應稀疏分解的分層介質參數反演

,,,

(南昌大學信息工程學院,江西南昌330031)

0 引言

目前關于探地雷達(Ground Penetrating Radar,GPR)電磁波在層狀體系中的傳播模型和介電特性模型反演缺乏一套嚴謹的理論,這嚴重制約了GPR無損技術的發展,使得GPR信號分析一直依賴簡化公式或根據經驗人工調試參數。現有研究成果大多數只能確保結構層單層厚度的精確檢測,對結構層多層厚度及結構層參數反演的研究和其他道路指標的檢測結果不盡人意。因此,開展高效準確的反演優化策略研究,探索層狀結構GPR檢測反演理論及其工程應用具有重要意義。

現有分層介質參數反演的手段主要分為時域和頻域[1]兩個方面。時域反演中傳統的參數反演算法包括通過估計的回波信號幅度和時延等信息進行介質參數反演的剝層(Layer Stripping)反演法和基于電磁模型參數優化過程理論的高斯牛頓(Gauss Newton)反演法[2],Caorsi和Stasolla提出一種基于能量比的層剝離技術(Energy-Based Layer Stripping)能夠在不同工作環境下重構分層介質的幾何和電性能參數,但要求相鄰的兩個回波信號之間的時間間隔比發射脈沖持續時間長[3]。頻域反演在薄層介質更具優勢,秦瑤等利用電磁波在層狀體系中的傳輸函數反演出地下層結構的厚度和廣義反射系數[4],但這種方法無法解析出地下層的電性參數。

有鑒于此,本文結合時間延遲完備字典和稀疏自適應匹配追蹤算法(Sparse Adaptive Matching Pursuit,SAMP)[5],以及基于能量比的層剝離技術提出了一種新的分層介質參數反演算法。該算法能夠對探地雷達信號進行稀疏分解,估計信號的反射系數和時間延遲,并在此基礎上實現分層介質的厚度和介電常數的估計。克服了層剝離估計精度較低和反卷積算法執行時間長等缺點,具有較短的執行時間和較高的估計精度。

1 解決方案

1.1 稀疏分解模型

GPR接收信號x(t)可建模為分層介質的廣義反射系數矢量a(t)與雷達發射信號s(t)的卷積加上測量噪聲及模型誤差n(t),寫成矩陣形式:

式中,X=[x(t1),x(t2),…,x(t N)]T,S=[s(t-τ1),s(t-τ2),…,s(t-τNT)]是一個經過Ψ描述后的與延時相關的托普利茲矩陣,α=[α1,α2,…,αNT]T,N=[n(t1),n(t2),…,n(t N)]T,廣義反射系數可以通過去卷積計算α=S-1X[6]。不過這種方法很費時。利用稀疏表示(Sparse Representation,SR)[7]可以很有效且快速地解決這個問題,SR是利用一個過完備字典S∈R N×NT,N＜N T中的少量原子線性的表示一個原始信號,即

式中,P0是一個NP-Hard[8]問題。實際應用中,一般采用貪婪算法或凸松弛算法解決此問題。

1.2 構造過完備字典

SR的關鍵步驟是構造過完備字典,過完備字典的類型有多種,本文采用匹配字典,以與信號最匹配為原則構造原子,減少了搜索匹配原子過程中的計算量。

為了求解式(1),假設探測場景的每個網格都是均勻的,則每一個網格點的目標回波信號時延是可以計算的,則接收信號可以表示為

式中,X=[x(t1),x(t2),…,x(t N)],t1為起始測量時間,t N=t1+(N-1)Δt,Δt為采樣時間間隔,N為時間采樣總數。則N×N T維的矩陣Ψ的第j列可表示為

式中,τ(πj)為收發器到第j個網格點目標的雙倍程時延,在不需要已知發射信號幅度的條件下,只需計算成像場景中各網格點到收發器的雙程時延即可構造需要的字典。

1.3 基于稀疏約束的廣義反射系數和時延估計

根據GPR信號傳播理論,當雷達波在地下傳播時,只有在不連續介質交界面處才會產生反射,因此廣義反射系數矢量α(t)是稀疏的,從而可以利用該先驗信息,將傳統的盲估計問題轉化為稀疏約束下的優化問題,即

式(5)可以被基于過完備字典的稀疏分解方法求出,SAMP算法的具體求解過程如下:

輸入:N×N T的時延過完備字典Ψ,N×1維觀測向量X,步長S。

輸出:信號的稀疏表示系數估計,含有各目標時延參數的原子集合Ψt1。

1)初始化r0=X,Λ0=?,L=S,t=1;

2)計算u=abs[ΨTr t-1](即計算〈r t-1,φj〉,1≤j≤N T),選擇u中L個最大值,將這些值對應的Ψ的列序號j構成集合S k(列集合序號);

3)令Ck=Λt-1∪Sk,Ψt={φj}(for allj∈Ck);

4)求X=Ψtαt的 最 小 二 乘 解:

5)從^αt中選出絕對值最大的L項記為,對應的Ψt中的L列記為ΨtL,對應的Λ的列 序號記為ΛtL,更新集合F=ΛtL;

7)如果殘差r tnew=0,則停止迭代進入第8)步;如果 ‖r tnew‖ ≥ ‖r t-1‖2,更新步長L=L+S,返回第2)步繼續迭代;前面兩個條件依次都不滿足,則Λt=F,r t=r tnew,t=t+1,如果t≤M則停止迭代進入第8)步,否則返回第2)步繼續迭代;

1.4 基于層剝離技術的分層介質參數反演

探地雷達收發系統距離分層介質上表面高度為d0,分層介質為K層水平分層的均勻非磁性介質。假設發射信號為s(t),x1(t)為第i-1層和第i層交界面的反射信號,εi和d i分別為第i層的介電常數和厚度,在TM極化方式下,第i-1層和第i層交界面處的反射系數和雙向透射系數分別為

式中,ρm-1,m為第m-1層和第m層交界面處的雙向透射系數,k m,k j為電磁波在第m層和第j層介質中的波數,d m,d j為第m層和第j層介質厚度,τk為第k層界面反射信號的雙倍程時延。估計出信號的時間延遲和廣義反射系數之后,就可以利用層剝離算法(Layer Stripping)反演出每層介質的介電常數和厚度,其詳細的反演過程如下:

1)由反卷積算法或SAMP算法估計反射系數矢量和各層界面時延;

5)知道τ1,τ2和ε1后可以估算出d1和k1,再由公式求出R1,2;

7)整個參數反演過程就是以上幾個步驟的重復迭代過程。

2 實驗結果分析

2.1 回波信號有重疊時參數反演結果

探測場景各層的參數為L0(ε0=1.0,d0=0.40 m),L1(ε1=6.0,d1=0.05 m),L2(ε2=10.0,d2=0.19 m),信號幅度和時延估計結果如圖1所示。

圖1 3種算法在回波信號有重疊情況下的信號重構和時延估計結果

表1給出了3種方法的參數反演結果和所需要的執行時間。

表1 回波信號有重疊時3種方法反演結果和所需要的執行時間

根據表中結果可知,當回波信號有重疊的情況下,反卷積方法對介質參數的估計結果是最精確的但很費時,Energy-Based Layer Stripping方法執行時間最短,但估計結果不精確,而SAMP方法估計結果及所需的執行時間都是比較合適的。因此,與反卷積和Energy-Based Layer Stripping方法相比SAMP方法有著明顯的優勢。

2.2 時間分辨率對處理結果的影響

定義BΔτ=B(2dε/c)為信號的時間處理分辨率,B為信號帶寬,Δτ為能夠區分兩個反射信號的最小時間延遲。在實驗中,給定帶寬時,d1和ε1決定了回波信號是否有混疊。圖2給出了當ε1=6.0固定不變,改變層厚度d1時回波信號幅度和時延估計結果,相應的各層介質參數反演結果如表2所示。

圖2 ε1=6.0固定不變,d1變化時回波信號幅度和時延估計結果

表2 ε1=6.0不變,d1變化時薄層介質參數反演結果

根據表中的參數反演結果,圖3給出了參數估計誤差與d1的關系曲線圖,從圖3可以看出當ε1不變時,分層介質參數估計誤差總的趨勢是隨d1增大,即BΔτ增大而減小。

圖3 參數估計誤差與d1的關系曲線

相同地,圖4給出了當d1=0.08 m固定不變,ε1變化時回波信號時延和幅度估計結果,相應的參數反演結果如表3所示,圖5給出了參數估計誤差與ε1的關系曲線圖。

圖4 d1=0.08 m固定不變,ε1變化時回波信號幅度和時延估計結果

表3 d1=0.08 m不變,ε1變化時薄層介質參數反演結果

圖5 參數估計誤差與ε1的關系曲線

由圖5可知,d1固定不變時,隨著ε1增大即BΔτ的增大,ε1,ε2的估計誤差減小。而厚度的估計誤差與ε1的變化不成比例關系,但總的趨勢是誤差隨BΔτ增大而減小。

3 結束語

本文主要研究GPR信號在分層介質中的傳播模型和各層介質的參數反演問題。研究結果表明,SAMP使用過完備字典中少量的原子就能估算出回波信號時延和幅度,大大減少了運算時間。在此基礎上利用Layer Stripping方法反演各層介質的參數。與傳統的反卷積方法和單純的Energy-Based Layer Stripping方法相比,該方法具有執行時間短和估計精度高等優點。

[1]HUANG Zhonglai,ZHANG Jianzhong.Determination of Parameters of Subsurface Layers Using GPR Spectral Inversion Method[J].IEEE Trans on Geoscience and Remote Sensing,2014,52(12):7527-7533.

[2]LIU Guangdong,ZHANG Kaiyin.A Time-Domain Gauss-Newton Inversion Algorithm for Solving Two-Dimensional Electromagnetic Inverse Scattering Problems[J].Acta Physica Sinica,2014,63(3):1-15.

[3]CAORSI S,STASOLLA M.A Layer Stripping Approach for EM Reconstruction of Stratified Media[J].IEEE Trans on Geoscience and Remote Sensing,2014,52(9):5855-5869.

[4]QIN Y,WANG Q.Using GPR Spectrum Inversion Method to Improve Recognition Ability of Thin-Layer[J].International Journal of Digital Content Technology and Its Applications,2012,6(10):136-143.

[5]康乃馨,何明浩,王冰切,等.基于壓縮感知的多徑LFM信號參數估計[J].雷達科學與技術,2016,14(3):291-296.KANG Naixin,HE Minghao,WANG Bingqie,et al.Parameter Estimation of Multipath LFM Signal Based on Compressive Sensing[J].Radar Science and Technology,2016,14(3):291-296.(in Chinese)

[6]ZHAO S,SHANGGUAN P,AL-QADI I L.Application of Regularized Deconvolution Technique for Predicting Pavement Thin Layer Thicknesses from Ground Penetrating Radar Data[J].NDT&E International,2015,73(3):1-7.

[7]鄭純丹,周代英.稀疏分解在雷達一維距離像中的應用[J].雷達科學與技術,2013,11(1):55-58.ZHENG Chundan,ZHOU Daiying.Research on High Resolution Range Profile Based on Sparse Decomposition in Radar[J].Radar Science and Technology,2013,11(1):55-58.(in Chinese)

[8]SHAO W,BOUZERDOUM A,PHUNG S L.Sparse Representation of GPR Traces with Application to Signal Classification[J].IEEE Trans on Geoscience and Remote Sensing,2013,51(7):3922-3930.