黃忠來，張建中

1廈門大學(xué)信息科學(xué)與技術(shù)學(xué)院，廈門 361005

2中國海洋大學(xué)海洋地球科學(xué)學(xué)院，青島 266100

1 引 言

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

2 頻譜反演方法原理

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

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

圖1所示的是包含4個反射系數(shù)的5層介質(zhì)模型，虛線表示地面，其z軸坐標(biāo)定為0，地下第i個反射層的深度用di表示，從上到下各層的介電常數(shù)分別為ε1、ε2、ε3和ε4，各層電導(dǎo)率分別為σ1、σ2、σ3和σ4，空氣的介電常數(shù)為ε0，電導(dǎo)率為σ0.w（t）為雷達(dá)發(fā)射波，w0（t）、w1（t）、w2（t）、w3（t）和w4（t）分別是在地表和地下各反射面的一次反射回波.這里，我們假設(shè)介質(zhì)為非磁性介質(zhì)，并僅考慮雷達(dá)波入射角接近0°，即單天線形式（Monostatic mode）雷達(dá)或發(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ù)是電磁波在介質(zhì)中傳播距離以及介質(zhì)電性參數(shù)的函數(shù).

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

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

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

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

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

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

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

其實(shí)部和虛部分別為

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

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

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

3 反演問題的解法

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

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

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

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

3.2 分步反演方法

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

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

其幅度譜為

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

解之得

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

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

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

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

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

其中，

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

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

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

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

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

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

4 實(shí)驗(yàn)例子

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

4.1 理論模型實(shí)驗(yàn)

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

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

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

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

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

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

4.2 實(shí)際資料應(yīng)用

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

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

為了驗(yàn)證本文反演算法對于薄層的效果，我們收集了某瀝青公路的GPR數(shù)據(jù).雷達(dá)天線發(fā)射頻率900MHz，每道采樣時間長度為15ns，采集512個數(shù)據(jù).雷達(dá)波在瀝青中的波長大約為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，平均電導(dǎo)率分別為6×10－6和0.005.

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

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

5 結(jié) 論

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

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

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

（References）

［1］ 馮德山，張彬，戴前偉等.基于速度估計的改進(jìn)型線性變換有限差分偏移在探地雷達(dá)中的應(yīng)用.地球物理學(xué)報，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ù)?全極化探地雷達(dá)正演模擬.地球物理學(xué)報，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］ 田鋼，林金鑫，王幫兵等.探地雷達(dá)地面以上物體反射干擾特征模擬和分析.地球物理學(xué)報，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ā).用地質(zhì)雷達(dá)數(shù)據(jù)資料反演二維地下介質(zhì)的方法.地球物理學(xué)報，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］ 丁亮，韓波，劉潤澤等.基于探地雷達(dá)的混凝土無損檢測反演成像方法.地球物理學(xué)報，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］ 張蓓.路面結(jié)構(gòu)層材料介電特性及其厚度反演分析的系統(tǒng)識別方法：路面雷達(dá)關(guān)鍵技術(shù)研究［博士論文］.重慶：重慶大學(xué)，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］ 李大心.探地雷達(dá)方法與應(yīng)用.北京：地質(zhì)出版社，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ā)，劉四新，王者江等.探地雷達(dá)方法原理及應(yīng)用.北京：科學(xué)出版社，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］ 張霖斌，紀(jì)晨，姚振興.疊后地震道反演的隨機(jī)爬山法.石油地球物理勘探，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.