黃忠來，張建中

1廈門大學信息科學與技術學院，廈門 361005

2中國海洋大學海洋地球科學學院，青島 266100

1 引 言

探地雷達（GPR）具有分辨率高、無損、快速等優(yōu)點，是探測近地表介質的主要技術之一.利用GPR，不僅可以得到地層的埋深和厚度，也可以獲取地下反射面的反射系數(shù)以及地層的介電常數(shù)和電導率等參數(shù)，為此，人們在正演［1－4］和反演方面［5－9］做了很多研究.目前層狀介質參數(shù)的反演方法大體可分為時域反演和頻域反演兩類.在時域反演算法中，Chien等［10］利用地震勘探中的共中心點方法（CMP），求取地下兩層介質的厚度和介電常數(shù)，張蓓［11］根據(jù)路面層狀結構，在時域對探地雷達波進行反演，得到層的復介電常數(shù)和厚度，但是這兩種方法都要求在時域能夠分辨出不同反射層的回波；另一些算法則要求有先驗信息作為基礎，如Huang等［12］采用參數(shù)校驗方法，通過在介電常數(shù)計算式中引入校驗參數(shù)來求取水平路基的層厚和介電常數(shù)，但需要事先利用鉆孔得到地下介質的有關參數(shù).在頻域反演方面，通過建立頻域正演模型，不斷改變模型參數(shù)，使模型計算的回波頻譜與實際數(shù)據(jù)頻譜一致，獲取地層的厚度和電性參數(shù).Qin等［13］通過求取電磁波在層狀介質中的傳輸函數(shù)，反演地層的厚度和廣義反射系數(shù)，但不能求得地層的電性參數(shù)；Minet以及 Lambot等［14－16］利用 Green函數(shù)建立目標函數(shù)，反演厚度和電性參數(shù)，但參數(shù)過多時，只有給出目標層的先驗信息后才能取得較好的結果.譜反演算法［17］是最近提出的利用地震記錄的頻譜，反演地層深度、厚度和反射系數(shù)的新方法.與時域反演方法相比，譜反演算法對薄層（如油氣儲層）有較強的適應性.本文在譜反演算法［17］的基礎上，結合電磁波在層狀介質中的傳播規(guī)律和GPR信號頻譜特征，提出了一種利用地下層狀介質GPR回波信號頻譜，同時估算介質幾何參數(shù)和電性參數(shù)的頻域反演算法；通過分析模型參數(shù)對GPR頻譜的影響規(guī)律，提出了分別單獨估計一個或幾個參數(shù)的分步反演方法，并以此結果作為進一步同時反演所有參數(shù)的初始值.最后，用理論模型的模擬數(shù)據(jù)和實測數(shù)據(jù)對方法進行了測試.

2 頻譜反演方法原理

GPR 測 量 目 標 往 往 處 于 天 線 場 遠 區(qū)［18－19］，且GPR多發(fā)射高頻電磁脈沖.利用傅里葉變換，電磁脈沖可以分解為一系列不同頻率的諧波，而這些諧波的傳播都可以近似為平面波的傳播形式.可見GPR的理論基礎是平面波在介質中的傳播規(guī)律［20］.因此，本文將雷達波作為平面波處理.對于GPR發(fā)射的電磁波來說，地下介質實際上是非理想的有耗介質，電磁波在傳播時會逐漸衰減，影響傳播的因素除了反射系數(shù)、透射系數(shù)外還有衰減常數(shù).當介質為非理想介質時，電導率σ≠0，等效介電常數(shù)為復數(shù)，此時波數(shù)傳播常量的定義為［21］

圖1 探地雷達波在地下傳播路徑示意圖Fig.1 Propagation path of GPR wave in layered media

圖1所示的是包含4個反射系數(shù)的5層介質模型，虛線表示地面，其z軸坐標定為0，地下第i個反射層的深度用di表示，從上到下各層的介電常數(shù)分別為ε1、ε2、ε3和ε4，各層電導率分別為σ1、σ2、σ3和σ4，空氣的介電常數(shù)為ε0，電導率為σ0.w（t）為雷達發(fā)射波，w0（t）、w1（t）、w2（t）、w3（t）和w4（t）分別是在地表和地下各反射面的一次反射回波.這里，我們假設介質為非磁性介質，并僅考慮雷達波入射角接近0°，即單天線形式（Monostatic mode）雷達或發(fā)射天線和接收天線之間的距離很小的情況.為清晰起見，圖中將入射波角度放大，并將回波路徑分開表示.第i層和i＋1層之間反射面的反射系數(shù)ri，i＋1和透射系數(shù)τi，i＋1分別為［22－23］：

且τi，i＋1＝1＋ri，i＋1.

式中，zi為波在第i層單向傳播的距離，α0＝0.可以看到，廣義反射系數(shù)是電磁波在介質中傳播距離以及介質電性參數(shù)的函數(shù).

定義了廣義反射系數(shù)后，我們構造用于反演目標層深度、厚度以及電性參數(shù)的代價函數(shù).以下考慮地下含有兩個水平反射界面的三層介質情況，如圖2a所示.若上、下廣義反射系數(shù)分別為r1和r2，將地面作為時間零點，上界面到地面的時間距離為t0，下界面到地面的時間距離為t1，中間層的時間厚度

圖2 三層介質反射系數(shù)模型Fig.2 Reflection coefficients of a three－layer model

為T，則反射系數(shù)序列可以表示為

如果將分析窗口的時間零點放在中間層的中點位置，那么反射系數(shù)序列的表達式變?yōu)?/p>

對上式進行傅里葉變換，得到該反射系數(shù)序列的頻譜：

將它表示成三角函數(shù)形式：

令re＝（r1＋r2）/2，ro＝ （r1－r2）/2，則［17］：

其實部和虛部分別為

若雷達發(fā)射子波頻譜為w（f），接收到的反射回波頻譜為s（f），則地下介質的反射系數(shù)序列的頻譜為

用雷達記錄反射系數(shù)序列實部和虛部分別與褶積模型反射系數(shù)序列實部和虛部之差的絕對值之和，定義的代價函數(shù)為

其中，abs表示絕對值；fH和fL分別是所用頻譜的頻率上限和下限，將根據(jù)實際數(shù)據(jù)中所含噪聲的分布情況確定，即盡量選取噪聲較少的頻帶.該代價函數(shù)的未知量包括目標層上反射面的時間位置t0，目標層的時間厚度T，廣義反射系數(shù)對的偶分量re和奇分量ro.由于廣義反射系數(shù)是電性參數(shù)的函數(shù)，把代價函數(shù)中的廣義反射系數(shù)奇、偶分量用電性參數(shù)替換，此時代價函數(shù)可表示成O（ε1，ε2，…，εn，σ1，σ2，…σn，t0，T）.通過求解代價函數(shù)的最小值問題，就可以得到地下各層的時間位置、厚度和電性參數(shù)，從而實現(xiàn)地下介質幾何參數(shù)和電性參數(shù)的同時反演.

3 反演問題的解法

3.1 全局優(yōu)化算法

對于求解像（17）式的多參數(shù)代價函數(shù)的極小值問題，我們采用在模擬退火算法上改進的隨機爬山法［24］.算法步驟為

（1）產生初始未知參數(shù)向量X＝｛x1，x2，…，xi，…，xN｝.每個參數(shù)的值可以在規(guī)定的范圍內隨機生成，也可以根據(jù)先驗信息給定.記初始參數(shù)向量對應的代價函數(shù)為O.

（2）在｛x1，x2，…，xi，…，xN｝中隨機選取一個參數(shù)進行修改.例如選取的是第i個變量xi，則修改后xi變?yōu)?，且＝xi＋sign·d·rand.sign是隨機選取的正負符號，－1或1.rand為0～1之間的隨機數(shù)，d為修改步長，用以控制修改參數(shù)的速度和精度.對于本文的幾何和電性參數(shù)，需給定不同的d值.

3.2 分步反演方法

通過分析有關參數(shù)與GPR回波信號頻譜之間的變化特征，提出了分別確定未知參數(shù)的分步反演方法.

考察廣義反射系數(shù)序列的頻譜：

其幅度譜為

對上式求導并令導數(shù)為0，即，

解之得

f表示幅度譜極大值點和極小值點的頻率，這些極大值點和極小值點分別是周期性分布的，相鄰極小值點頻率之差，即幅度譜凹陷周期為

（22）式表明，幅度譜凹陷周期由時間厚度T決定，而與其它參數(shù)無關.例如當T分別取T1＝6.67ns和T2＝3.33ns，其它參數(shù)保持不變時，幅度譜凹陷周期分別為Δf1＝150MHz和Δf2＝300MHz，如圖3所示.利用這一點，我們就可以在不考慮其它參數(shù)的情況下，僅利用反射系序列幅度譜凹陷周期來確定T.當層的厚度很小時，反射系數(shù)幅度譜凹陷周期將很大，若超出所設定的反演頻帶，就不能使用（22）式估算T.這時，可以找出幅度譜的第一個極大值點頻率，峰值頻率fp，再利用（21）式估算T，在式（21）中令n＝2，即得到T＝1/（2fp）.

從（21）和（22）式可以看出，隨著層厚度的減小，反射系數(shù)幅度譜的第一個峰值頻率和幅度譜凹陷周期將不斷增大.即使當層厚度小于調諧厚度時，反射系數(shù)幅度譜也有這樣明顯的變化規(guī)律，從而可以利用頻譜分辨或反演小于調諧厚度的薄層.

圖3 幅度譜凹陷周期與層時間厚度T的關系Fig.3 Relationship between amplitude spectral notches period and layer′s temporal thickness T

將相位譜分為兩個部分，第一部分不包含參數(shù)t，第二部分不包含反射系數(shù)的奇偶分量：

其中，

由于ro/re與頻率f無關，所以相位譜∠r（f）隨頻率f變化的快慢由t0和T決定，且θ1隨頻率f的變化速度比θ2慢2πft0.當T按上述方法被確定后，就可以利用相位譜的變化快慢估算t0.為說明ro/re對于相位譜的影響，我們隨機給定一組參數(shù)向量，并使其它參數(shù)保持不變，只改變ro/re的值，相應的相位譜如圖4所示.Δf是相位發(fā)生跳變前持續(xù)的頻率間隔，基本不隨ro/re變化，這說明ro/re比值主要影響相位譜的幅值，而對相位譜變化快慢的影響很小.這樣，在確定T參數(shù)后，就可以單獨估算t0.

圖4 相位譜隨ro/re值的變化Fig.4 Variation of phase spectrum with the value of ro/re

考察相位譜中不含參數(shù)t0的項θ1.由式（3）可以得，當?shù)叵掠腥龑咏橘|時，

由于r1和r2只與介電常數(shù)有關，α2由第i＋2層介質的電導率以及介電常數(shù)決定，z2由T以及第2層的介電常數(shù)決定.即當T和t0被確定以后，影響相位譜的只有3層的介電常數(shù)和第2層的電導率.這樣，僅利用相位譜就可以反演三個介電常數(shù)和中間層的電導率.

通過上面分析可知，不同參數(shù)對反射系數(shù)序列頻譜的影響是不同的.T單獨決定了幅度譜凹陷周期大小，t0和T共同影響相位譜的變化快慢，t0和各層的介電常數(shù)以及中間層電導率則完全確定了相位譜.因此，可以這樣分別確定這些參數(shù)：首先根據(jù)幅度譜的凹陷周期確定厚度參數(shù)T，再根據(jù)相位譜的變化速度反演時間位置參數(shù)t0，然后根據(jù)相位譜反演3層的介電常數(shù)和中間層的電導率，最后反演上層的電導率.采用隨機爬山法分別反演t0、三個介電常數(shù)以及中間層電導率.

這里用理論數(shù)據(jù)來說明分步求取參數(shù)的有效性.圖5和圖6中的實線是利用正確參數(shù)向量生成的反射系數(shù)序列的頻譜；點線是利用估算的參數(shù)向量生成的反射系數(shù)序列的頻譜.其中，圖5（a，b）中的點線分別是由隨機生成的一組參數(shù)向量計算出的反射系數(shù)序列頻譜的幅度譜和相位譜，與理論頻譜相差較大；圖5c中點線是利用幅度譜凹陷周期估計T后計算的振幅譜.可以看到在其它參數(shù)保持為原隨機生成的參數(shù)值的情況下，估算的T使幅度譜凹陷周期與理論一致，當然相位譜仍然相差很大，如圖5d所示.圖6（a，b）中的點線是在圖5估算的T的基礎上，利用相位譜估算t0，且其它參數(shù)不變時計算的幅度譜和相位譜.可以看出，估算的t0值使相位譜更接近理論值，且t0的改變并不影響幅度譜的凹陷周期.圖6（c，d）中的點線是由前面確定的T和t0值以及繼續(xù)反演三層的介電常數(shù)以及中間層的電導率后計算的幅度譜和相位譜，這時計算的頻譜與理論頻譜基本一致，其小的差別是因為上層的電導率沒有估算而隨機選取的緣故.可以看出，這樣求出的介質參數(shù)值與理論參數(shù)值比較接近，計算的頻譜也與理論頻譜比較接近.這種分步求取參數(shù)的好處是可以極大地減少每次反演的未知量個數(shù)，降低收斂到局部最小值的幾率，提高運算效率和反演結果的可靠性.特別是能為在沒有先驗信息情況下的多參數(shù)全局優(yōu)化算法提供接近真解的初值.

4 實驗例子

分別利用理論模型合成數(shù)據(jù)和實測資料對本文方法進行了測試和應用.理論模型采用楔形目標模型，以測試反演方法對不同厚度目標層的有效性.實測資料是收集到的公路路面GPR檢測數(shù)據(jù).

4.1 理論模型實驗

使用GPRMax［25］建立了楔形目標層模型并合成了理論GPR數(shù)據(jù).發(fā)射雷克子波中心頻率為300MHz；上中下三層的介電常數(shù)ε1、ε2和ε3分別為9、16和9，電導率σ1、σ2和σ3分別為0.02、0.01和0.02；楔形層上界面的反射系數(shù)為－0.25，下界面反射系數(shù)為0.25；楔形層上界面距離地表0.38m，楔形層的最大厚度為0.125m，是雷達發(fā)射波在楔形中波長的1/2，1/4波長的調諧厚度位于第25道處.模型及楔形層上下反射面位置、楔形厚度以及各層電性參數(shù)的反演結果如圖7和圖8所示.可以看到，當目標厚度小于調諧厚度時，依然可以得到正確的反演結果.

現(xiàn)在以模型第35道為例，說明采用分步反演結果作為初始值的效果.圖9給出了分別采用不同初始值時，代價函數(shù)隨迭代次數(shù)的變化.實線為初始值采用分步反演結果，當?shù)螖?shù)為1200次時就已收斂，對應的代價函數(shù)值是0.2834；點劃線的初值是隨機生成的，當?shù)螖?shù)超過2500時才收斂，對應的代價函數(shù)是0.4269.顯然，用分步反演結果作為初值時，反演收斂得更快且收斂于更小的代價函數(shù)值.表1列出了參數(shù)理論值、用作初值的分步反演結果和同時反演結果，圖10是用該反演結果計算的幅度譜（a）及相位譜（b）與理論頻譜之間的對比.

從上述反演結果可以看到，當楔形層的厚度大于1/8波長（即圖7中第37道位置）時，本文方法可以給出較準確的上反射面時間位置和楔形層的時間厚度，電性參數(shù)也基本都在理論值附近；當楔形層厚度小于1/8波長后，反演結果誤差較大.總的來說，深度和厚度參數(shù)的反演結果優(yōu)于電性參數(shù)的結果.對目標層的時間厚度T的反演較準確，這是因為幅度譜凹陷周期大小主要與它相關，而與其它參數(shù)的聯(lián)系不密切，分步反演給出了非常接近真值的初始值.同樣，使用與反射面位置參數(shù)t0相關而與其它參數(shù)不相關的相位譜屬性，使t0的反演結果也比較準確.介電常數(shù)既影響相位譜，又影響幅度譜.當需要反演的介電常數(shù)個數(shù)較多時，反演結果就容易陷入局部最優(yōu)，這也是介電常數(shù)的反演結果沒有T和t0結果準確的主要原因.相比之下，電導率反演結果不夠準確，因為它決定雷達波的衰減程度，主要體現(xiàn)在雷達信號的幅值變化上，但幅值不僅與電導率有關，還與反射面上、下介電常數(shù)的差值以及層厚度有關系.

表1 采用分步反演數(shù)據(jù)為初始值的參數(shù)反演結果Table 1 Parameter results of inversion starting with initial values given by the staged inversion method

為了進一步說明本文分步反演方法的性能，將本文方法與參考文獻［11］中的頻域反演方法做了比較.文獻［11］中的方法是針對地下三層介質，在地下第一層介質厚度、介電常數(shù)和電導率已知的情況下，反演第二層厚度和兩個反射界面的廣義反射系數(shù).因此，設計了中間有一薄層的三層模型.第一層厚度為40cm，中間薄層厚度為3cm，約為電磁波在層中波長的1/7.從上到下三層的相對介電常數(shù)分別是9，25，9，電導率為0.02，0.01，0.02.發(fā)射天線中心頻率為300MHz.分二種情況對這兩種方法進行比較：方案一是按照文獻［11］方法已知第一層的厚度和電性參數(shù)，隨機給定其余參數(shù)的初值進行反演；方案二是隨機給定包括第一層介質參數(shù)的所有參數(shù)的初值進行反演.兩種方法反演過程的代價函數(shù)隨迭代次數(shù)的變化曲線如圖11所示，反演出的中間薄層的厚度見表2.可以看出，當?shù)谝粚拥暮穸群碗娦詤?shù)已知時，兩種反演方法都能很快收斂，且反演結果與理論值誤差較小.當所有參數(shù)初始值隨機選取時，兩種方法達到收斂所需要的迭代次數(shù)增加，本文方法大約在迭代1200次時收斂，這時反演的薄層厚度與理論值的誤差為0.26cm，文獻［11］方法則大約在迭代2300次時才收斂，且代價函數(shù)仍然比前者大很多，反演的薄層厚度與理論值的誤差為0.35cm.可見，在沒有先驗信息獲得較好初值的情況下，本文的分步反演方法可以提高反演精度和效率.

表2 用兩種方法反演薄層厚度（單位：cm）Table 2 Thin－layer thickness inverted using the two algorithms

4.2 實際資料應用

圖11 兩種反演方法代價函數(shù)隨迭代次數(shù)的變化Fig.11 Cost functions versus iteration number for the two algorithms

第一個實測GPR數(shù)據(jù)來自廣東某公路的混凝土路面檢測資料.公路設計為水泥穩(wěn)定碎石基層加水泥混凝土路面，采用C30混凝土，水泥層設計厚度為26cm，基層設計厚度為22cm.通常水泥層介電常數(shù)為5.5，電導率約為0.01～0.05；碎石介電常數(shù)約為5～15，電導率為3～20.實際鉆芯位置在K0＋000－K1＋000處，鉆探顯示水泥的平均厚度為23cm.GPR天線發(fā)射頻率為900MHz，每道采樣時間長度為15ns，采集512個數(shù)據(jù).圖12是GPR實測數(shù)據(jù)剖面及反演結果，其中介電常數(shù)分別為4、6.5和15，電導率分別為0.01和0.025，把時間單位換算成長度單位，水泥層的平均厚度為24cm，基層平均厚度為21cm，反演結果與鉆探結果基本一致，但水泥層厚度卻略小于原設計厚度，也說明施工的水泥層厚度略微不夠.

為了驗證本文反演算法對于薄層的效果，我們收集了某瀝青公路的GPR數(shù)據(jù).雷達天線發(fā)射頻率900MHz，每道采樣時間長度為15ns，采集512個數(shù)據(jù).雷達波在瀝青中的波長大約為0.167m，瀝青層上面層的厚度為0.04m，略小于1/4波長，因此可以被視為薄層.在時域剖面圖上，瀝青上面層的下底面回波與地表回波疊加在一起，很難對兩者做出區(qū)分，而本文的頻譜反演算法能給出瀝青上層面的厚度.圖13為第600道數(shù)據(jù)的上面層回波反射系數(shù)幅度譜，由于該層厚度小，難以確定幅度譜凹陷周期，但可以獲取幅度譜的第一個峰值頻率fp為940MHz，根據(jù)式（21）得到T 的初值為0.53ns，約0.041m.通過全局優(yōu)化反演后得到該段瀝青路面的上面層和下面層的深度和厚度如圖14所示，其中反演的從上至下各層的平均相對介電常數(shù)的分別為0.83、5.3和3.8，平均電導率分別為6×10－6和0.005.

圖12 某混凝土公路檢測的GPR剖面（a）及反演結果（b）Fig.12 GPR profile from a highway detection and corresponding inversion result

圖13 實際反射系數(shù)幅度譜和根據(jù)反演模型計算的幅度譜Fig.13 GPR amplitude spectra obtained from measured data and calculated using inversion model，respectively

5 結 論

通過對電磁波在層狀介質中的傳播規(guī)律和單天線模式或發(fā)射天線與接收天線間距很小時的GPR信號頻譜特征的研究，基于地下層狀介質GPR回波信號頻譜及全局優(yōu)化算法，提出了一種同時估算介質幾何參數(shù)和電性參數(shù)的頻域反演方法.理論模型和實際資料測試結果表明，該方法效果良好，具有較高的分辨率，對厚度為1/8發(fā)射脈沖波長左右的薄層介質仍有較好結果，而且對幾何參數(shù)的反演精度優(yōu)于對電性參數(shù)的反演.

地下層狀介質的不同參數(shù)對于反射系數(shù)序列頻譜屬性的影響不同.如目標層時間厚度決定了幅度譜凹陷周期的大小，目標層時間厚度和上頂界面時間位置決定相位譜的變化快慢，上頂界面時間位置和各層的介電常數(shù)以及中間層電導率共同影響相位譜.據(jù)此，我們提出了利用不同的頻譜屬性分別估算不同未知參量的分步反演方法.以該分步反演結果作為全局優(yōu)化反演的初值，可以大大提高反演的收斂速度和反演結果的可靠性.

圖14 某瀝青路面的GPR剖面（a）及反演結果（b）Fig.14 GPR profile of asphalt pavement and corresponding inversion result

（References）

［1］ 馮德山，張彬，戴前偉等.基于速度估計的改進型線性變換有限差分偏移在探地雷達中的應用.地球物理學報，2011，54（5）：1340－1347.Feng D S，Zhang B，Dai Q W，et al.The application of the improved linear transformation of finite difference migration based on the velocity estimation in the GPR date processing.Chinese J.Geophys.（in Chinese），2011，54（5）：1340－1347.

［2］ 馮晅，鄒立龍，劉財?shù)?全極化探地雷達正演模擬.地球物理學報，2011，54（2）：349－357.Feng X，Zou L L，Liu C，et al.Forward modeling for fullpolarimetric ground penetrating radar.Chinese J.Geophys.（in Chinese），2011，54（2）：349－357.

［3］ 田鋼，林金鑫，王幫兵等.探地雷達地面以上物體反射干擾特征模擬和分析.地球物理學報，2011，54（10）：2639－2651.Tian G，Lin J X，Wang B B，et al.Simulation and analysis reflections interference from above surface objects of ground penetrating radar.Chinese J.Geophys.（in Chinese），2011，54（10）：2639－2651.

［4］ Huang Y Q，Zhang J Z，Liu Q H.Three－dimensional GPR ray tracing based on wavefront expansion with irregular cells.IEEE Transactions on Geoscience and Remote Sensing，2011，49（2）：679－687.

［5］ Huang Y Q，Liu Y H，Liu Q H，et al.Improved 3－D GPR detection by NUFFT combined with MPD method.Progress in Electromagnetics Research，2010，103：185－199.

［6］ Liu C R，Li J，Gan X，et al.New model for estimating the thickness and permittivity of subsurface layers from GPR data.IEE Proc.－Radar Sonar Navig，2002，149（6）：315－319.

［7］ Meles G A，Van der Kruk J，Greenhalgh S A，et al.A new vector waveform inversion algorithm for simultaneous updating of conductivity and permittivity parameters from combination crosshole/borehole－to－surface GPR Data.IEEE Transactionson Geoscienceand Remote Sensing，2010，48（9）：3391－3407.

［8］ 王兆磊，周輝，李國發(fā).用地質雷達數(shù)據(jù)資料反演二維地下介質的方法.地球物理學報，2007，50（3）：897－904.Wang Z L，Zhou H，Li G F.Inversion of ground－penetrating radar data for 2D electric parameters.Chinese J.Geophys.（in Chinese），2007，50（3）：897－904.

［9］ 丁亮，韓波，劉潤澤等.基于探地雷達的混凝土無損檢測反演成像方法.地球物理學報，2012，55（1）：317－326.Ding L，Han B，Liu R Z，et al.Inversion imaging method for concrete non－destructive testing based on GPR.Chinese J.Geophys.（in Chinese），2012，55（1）：317－326.

［10］ Kao C P，Li J，Wang Y Q，et al.Measurement of layer thickness and permittivity using a new multilayer model from GPR data.IEEE Transactions on Geoscience and Remote Sensing，2007，45（8）：2463－2470.

［11］ 張蓓.路面結構層材料介電特性及其厚度反演分析的系統(tǒng)識別方法：路面雷達關鍵技術研究［博士論文］.重慶：重慶大學，2003.Zhang B.System identification method for backcalculating the dielectric property and thickness of pavement structures：study on applied technology of ground penetrating radar（in Chinese）［Doctor′s thesis］.Chongqing：Chongqing University，2003.

［12］ Huang C L，Su Y.A new GPR calibration method for high accuracy thickness and permittivity measurement of multilayered pavement.∥Proceedings of the Tenth International Conference on Ground Penetrating Radar.Delft，The Netherlands：IEEE，2004：627－630.

［13］ Qin Y，Chen J，F(xiàn)ang G Y，et al.Research on thin－layer recognition technique based on the spectrum inversion method of ground penetrating radar.Journal of Electronics ＆Information Technology，2010，32（11）：2760－2763.

［14］ Minet J，Patriarca C，Slob E C，et al.Characterization of layered media using full－waveform inversion of proximal GPR data.∥2010URSL International Symposium on Electromagnetic Theory.Berlin：IEEE，2010：1004－1007.

［15］ Minet J，Lambot S，Delaide G，et al.A generalized frequency domain reflectometry modeling technique for soil electrical properties determination.Vadose Zone Journal，2010，9（4）：1063－1072.

［16］ Patriarca C， Lambot S， Mahmoudzadeh M R，et al.Reconstruction of sub－wavelength fractures and physical properties of masonry media using full－waveform inversion of proximal penetrating radar.Journal of Applied Geophysics，2011，74（1）：26－37.

［17］ Puryear C.I，Castagna J P.Layer－thickness determination and stratigraphic interpretation using spectral inversion：Theory and application.Geophysics，2008，73（2）：R37－R48.

［18］ 李大心.探地雷達方法與應用.北京：地質出版社，1994.Li D X.Methods and Applications of Ground Penetrating Radar（in Chinese）.Beijing：Geological Publishing House，1994.

［19］ Balanis C A.Antenna Theory：Analysis and Design.New Jersey：John Wiley，2005.

［20］ 曾昭發(fā)，劉四新，王者江等.探地雷達方法原理及應用.北京：科學出版社，2006.Zeng Z F，Liu S X，Wang Z J，et al.Principles and Applications of GPR Method（in Chinese）.Beijing：Science Press，2006.

［21］ Guru B S，H?z?rolu H R.Electromagnetic Field Theory Fundamentals.Cambridge： Cambridge University Press，2004.

［22］ Schmelzbach C，Scherbaum F，Tronicke J，et al.Bayesian frequency－domain blind deconvolution of ground－penetrating radar data.Journal of Applied Geophysics，2011，75（4）：615－630.

［23］ Jol H M.Ground Penetrating Radar：Theory and Applications.Singapore：Elsevier，2009.

［24］ 張霖斌，紀晨，姚振興.疊后地震道反演的隨機爬山法.石油地球物理勘探，1997，32（1）：75－80.Zhang L B，Ji C，Yao Z X.Stochastic hill－climbing algorithm for poststack seismic trace inversion.Oil Geophysical Prospecting （in Chinese），1997，32（1）：75－80.

［25］ Giannopoulos A. Modelling ground penetrating radar by GprMax.Construction and Building Materials，2005，19（10）：755－762.