黃忠來 張建中 黃吉林

(1)廈門大學(xué)信息科學(xué)與技術(shù)學(xué)院,廈門 361005 2)中國海洋大學(xué)海洋地球科學(xué)學(xué)院,青島 266100)

探地雷達(dá)薄層信號(hào)的譜反演算法＊

黃忠來1)張建中2)黃吉林1)

(1)廈門大學(xué)信息科學(xué)與技術(shù)學(xué)院,廈門 361005 2)中國海洋大學(xué)海洋地球科學(xué)學(xué)院,青島 266100)

為了從多層介質(zhì)的探地雷達(dá)回波數(shù)據(jù)中提取薄層的位置和厚度信息,推導(dǎo)和驗(yàn)證了在頻率域內(nèi)實(shí)現(xiàn)的譜反演算法。該算法給出了各反射面的廣義反射系數(shù),并建立了用于聯(lián)系反射系數(shù)序列頻譜和薄層參數(shù)的代價(jià)函數(shù)并進(jìn)行反演。仿真結(jié)果表明,對(duì)于厚度小于調(diào)諧厚度的薄層,該算法也可以進(jìn)行準(zhǔn)確反演,從而提高了對(duì)于薄層的分辨能力。

探地雷達(dá);薄層;譜反演算法;廣義反射系數(shù);代價(jià)函數(shù)

1 引言

探地雷達(dá) (Ground Penetrating Radar,GPR)通過向地下發(fā)射電磁脈沖,并接收地下回波信號(hào)來探測地下目標(biāo),它具有分辨率高、無損、快速等優(yōu)點(diǎn)。本文研究的對(duì)象是含有薄層的層狀結(jié)構(gòu)地下介質(zhì),目的是通過反演算法得到地下各反射面的位置以及薄層的厚度,這里的薄層是指厚度小于電磁波在其中波長 1/4的層。對(duì)于層狀介質(zhì)參數(shù)的反演,文獻(xiàn)[1]針對(duì)地震勘探共中心點(diǎn)方法 (CMP)中的模型進(jìn)行改進(jìn),并在時(shí)域反演,以求取地下兩層介質(zhì)的厚度和介電常數(shù),但是反演時(shí)需要在時(shí)域得到反射面回波的到達(dá)時(shí)間。文獻(xiàn)[2]研究路面層狀結(jié)構(gòu),利用雷達(dá)波傳播規(guī)律進(jìn)行時(shí)域反演,得到層的復(fù)介電常數(shù)和厚度,但是同樣沒有考慮目標(biāo)層是薄層的情況。本文在文獻(xiàn)[3]提出的地震勘探譜反演算法的基礎(chǔ)上,結(jié)合電磁波在層狀介質(zhì)中傳播的特點(diǎn),提出用廣義反射系數(shù)代替原算法中的反射系數(shù),推導(dǎo)并驗(yàn)證了對(duì)地下薄層 GPR信號(hào)進(jìn)行反演的頻域算法,提高了對(duì)地下薄層的分辨力。

2 譜反演算法過程

設(shè) re=(r1+r2)/2,ro=(r1-r2)/2,則：

其中 Re[g(f)]=2recos(πfT),Im[g(f)]=2rosin(πfT)分別是實(shí)部和虛部。從定義式可以看出 re是函數(shù) f(r1,r2)=r1δ(t-t1)+r2δ(t-t2)做奇偶分解后的偶部分,ro則是奇部分(圖 2)。

令 G(f)表示 g(f)的幅度譜,則：

圖1 反射系數(shù)對(duì)Fig.1 Reflection coefficients pairs

圖2 反射系數(shù)對(duì)奇偶分解Fig.2 Even-odd decomposition of coefficients pairs

對(duì) G(f)取導(dǎo)數(shù)得：

將 (3)、(4)兩式相乘得：

找到使得代價(jià)函數(shù)值為零的參數(shù)和后,就可以根據(jù)

求出上、下反射系數(shù)的奇偶分量,進(jìn)而得到反射系數(shù)值。

由于事先并不知道反射系數(shù)譜 g(f),所以在計(jì)算 G(f)時(shí)要利用接收到的雷達(dá)回波和雷達(dá)發(fā)射子波的頻譜。此外,算法建立代價(jià)函數(shù)的前提是分析窗口時(shí)間零點(diǎn)位置位于兩反射系數(shù)的中點(diǎn)。下面,將在此算法的基礎(chǔ)上推導(dǎo)適用于多反射系數(shù)的代價(jià)函數(shù),本文中所推導(dǎo)出的函數(shù)表達(dá)式和文獻(xiàn)[3]中的有些不同。

以包含 4個(gè)反射系數(shù)的模型為例,圖 3中反射系數(shù) r1～r4分別為 0.1、0.3、-0.2、0.4,t1～t4分別是 4個(gè)反射系數(shù)的距離地面的時(shí)間位置,T1、T2是前后兩對(duì)反射系數(shù)的時(shí)間厚度,則反射系數(shù)序列g(shù)(t)及其頻譜 g(f)可以表示為 g(t)=r1(t)+r2(t) =[r1δ(t-t1)+r2δ(t-t1-T1)]+[r3δ(t-t3)+r4δ (t-t3-T2)],以及 g(f)=r1(f)+r2(f)。

圖3 反射系數(shù)對(duì)及其奇偶分解Fig.3 Coefficients pairs and their even-odd decompositions

由前面兩反射系數(shù)的推導(dǎo)可知：

將 4個(gè)反射系數(shù) (兩個(gè)反射系數(shù)對(duì))推廣到 n個(gè)反射系數(shù)對(duì),則代價(jià)函數(shù)可表示為：

此代價(jià)函數(shù)含有以下未知量：反射系數(shù)對(duì)中,上層反射系數(shù)的時(shí)間位置 t,兩反射系數(shù)之間時(shí)間厚度T,反射系數(shù)對(duì)的偶分量 re以及反射系數(shù)對(duì)的奇分量 ro。以上就是譜分解算法的大致過程。通過最小化代價(jià)函數(shù)的值,可以得到地下各層的時(shí)間位置、厚度以及反射系數(shù)。但是,算法中討論的是簡單的反射系數(shù)褶積模型,未討論波在地下傳播時(shí)的衰減。對(duì)于探地雷達(dá)波來說,反演出的反射系數(shù)并不等于反射面真正的反射系數(shù),因?yàn)榈叵陆橘|(zhì)實(shí)際上是非理想的有耗介質(zhì),導(dǎo)致電磁波在傳播時(shí)會(huì)逐漸衰減,影響其傳播的除了反射系數(shù)、透射系數(shù)外還有衰減常數(shù)。因此原代價(jià)函數(shù)里的反射系數(shù)奇偶分量需要重新定義。

3 廣義反射系數(shù)

同理,第二個(gè)反射界面的廣義反射系數(shù)為：

以此類推,第 i個(gè)廣義反射系數(shù)為：

式中 zi為波在第 i層中單向傳播的距離,c為光速,并認(rèn)為α0=0。

圖4 電磁波在地下傳播路徑示意圖Fig.4 Electromagnetic wave propagation path under ground

4 反演算法

對(duì)于包含多參數(shù)的代價(jià)函數(shù)的反演,本文采用在模擬退火算法上改進(jìn)得到的隨機(jī)爬山法[4]。其算法步驟為：1)產(chǎn)生初始參數(shù)向量 X={x1,x2,…, xi,…,xN},每個(gè)參數(shù)的值可以在規(guī)定的范圍內(nèi)隨機(jī)生成,也可以根據(jù)實(shí)際情況下得到的數(shù)據(jù)進(jìn)行設(shè)定。這時(shí),向量對(duì)應(yīng)的代價(jià)函數(shù)為Obj。

2)在{x1,x2,…,xi,…,xN}中隨機(jī)選取一個(gè)參數(shù)進(jìn)行修改,例如選取的是第 i個(gè)變量 xi,則修改后xi變?yōu)榍?xi+sign＊d＊ rang,sign是在 -1和 1兩者間隨機(jī)選取的符號(hào),d為修改的步長,用以控制反演的速度和精度,rand為 0～1之間的隨機(jī)數(shù)。

重復(fù) 2)到 4)的操作,直到最后計(jì)算出的代價(jià)函數(shù)值滿足收斂條件,或者反演的次數(shù)超過設(shè)定值為止。每次修改時(shí)的步長 d不是固定的,開始時(shí)為了加快收斂速度,步長可以取得較大,以后逐漸減小,從而提高每次修改的精度,盡量使算法收斂到全局最優(yōu)值。

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

5.1 楔形模型仿真

為討論目標(biāo)層厚度對(duì)算法結(jié)果的影響,利用GPRMax軟件[5]進(jìn)行正演仿真,產(chǎn)生作為實(shí)際參考值的數(shù)據(jù)。建立上、下層面反射系數(shù)同號(hào)和異號(hào)兩種情況下的楔形模型A和B(圖 5和圖 6分別是A和B的時(shí)域回波圖)。模型 A中,令楔形上界面反射系數(shù)為 -0.25,下界面為 0.25。模型 B中的楔形,令上界面反射系數(shù) -0.14,下界面為 -0.11。楔形厚度最厚處為發(fā)射波在楔形中波長的 1/2,雷達(dá)發(fā)射波中心頻率設(shè)定為 300 MHz(圖 5～8)。

圖7和圖 8分別是楔形模型 A和 B的反射系數(shù)反演結(jié)果。從圖7和圖8可以看到,在兩種反射系數(shù)組合下的模型中,1/4波長調(diào)諧厚度均位于記錄道第 25道,當(dāng)楔形厚度逐漸減小到 1/8波長時(shí),算法仍然可以正確反演出反射面的位置和厚度,這說明了算法對(duì)于薄層的分辨效果非常好。

圖6 模型B回波Fig.6 Received waves ofmodelB

圖7 模型A反射系數(shù)反演結(jié)果Fig.7 Inversion results of coefficients in modelA

圖8 模型B反射系數(shù)反演結(jié)果Fig.8 Inversion results of coefficients in modelB

5.2 多反射系數(shù)層狀模型

設(shè)計(jì)模型在地下有 4個(gè)反射界面,反射系數(shù)分別是 r1=-0.250 0、r2=0.111 1、r3=0.200 0、r4= 0.500 0。其相應(yīng)的廣義反射系數(shù)分別是 gr1= -0.113 1、gr2=0.043 7、gr3= -0.030 6、gr4= 0.069 9。地下包含了 2個(gè)層,上方層的厚度設(shè)計(jì)成電磁波在層中波長的 1/10,下方層的厚度設(shè)計(jì)成波長的 1/2,即上方是個(gè)相對(duì)薄層,而下方是相對(duì)厚層。

表1 實(shí)際值與反演結(jié)果比較Tab.1 Comparison between i nversion results and true values

圖9 回波、發(fā)射波及反射系數(shù)序列頻譜Fig.9 Spectrum of reflected wave,incident wave and reflection coefficients

圖10 反演結(jié)果與實(shí)際值比較Fig.10 Comparison between inversion results and true values

圖9是雷達(dá)回波、發(fā)射波的頻譜以及由前兩者計(jì)算出的模型反射系數(shù)序列頻譜,圖 10是由反演出的參數(shù)計(jì)算得到的頻譜和真實(shí)反射系數(shù)序列頻譜之間的對(duì)比,上方是兩個(gè)頻譜的實(shí)部,下方是虛部,其中紅色星劃線是真實(shí)值,藍(lán)色圈劃線是反演值。從表 1的多個(gè)反射面模型反演結(jié)果可以看到,薄層的位置以及厚度的反演結(jié)果準(zhǔn)確,表明當(dāng)?shù)叵聻槎鄬咏Y(jié)構(gòu)時(shí),算法仍然正確有效。其中薄層的廣義反射系數(shù)反演值與真實(shí)值存在一定誤差。

6 結(jié)語

通過研究電磁波在地下層狀介質(zhì)中的傳播規(guī)律,重新定義了反射系數(shù)表達(dá)式,推導(dǎo)并驗(yàn)證了適用于 GPR信號(hào)的譜反演算法。正是因?yàn)橛辛诵碌膹V義反射系數(shù)定義式,才可以根據(jù)實(shí)際模型參數(shù)來驗(yàn)證算法反演結(jié)果的正確性。利用譜反演算法可以得到厚度小于 1/8波長的薄層的反射面位置,層厚以及每個(gè)反射面廣義反射系數(shù)的大小。仿真結(jié)果表明該算法對(duì)于提高薄層的分辨有很好的效果,但是在優(yōu)化算法的結(jié)果上存在多解性問題,有進(jìn)一步改進(jìn)的余地。

1 Kao Chienping,et al.Measurement of layer thickness and permittivity using a new multilayer model from GPR data [J].IEEE Transactions on geoscience and remote sensing, 2007,45(8)：2 463-2 470.

2 張蓓.路面結(jié)構(gòu)層材料介電特性及其厚度反演分析的系統(tǒng)識(shí)別方法——路面雷達(dá)關(guān)鍵技術(shù)研究 [D].重慶大學(xué),2003. (Zhang Bei.System identification method for back-calculating the dielectric property and thickness of pavement structures_study on applied technology of ground penetrating radar[D].ChongqingUniversity,2003)

3 Charles I Puryear and John P Castagna.Layer-thickness deter mination and stratigraphic interpretation using spectral inversion：Theory and application[J].Geophysics,2008, 73(2)：37-48.

4 張霖斌,紀(jì)晨,姚振興.疊后地震道反演的隨機(jī)爬山法[J].石油地球物理勘探,1997,32(1)：75-80.(Zhang Linbin,JiChen and Yao Zhenxing.Stochastic hill-climbing algorithm of post-stack seismic inversion[J].Oil Geophysical Prospecting,1997,32(1)：75-80)

5 Giannopoulos A.Modelling ground penetrating radar by Gpr Max[J].Construction and Building Materials,2005, 19：755-762.

6 秦德文.基于譜反演的薄層預(yù)測與反演方法研究[D].中國石油大學(xué),2009.(Qin Dewen.Thin bed prediction using spectral inversion and inversion method analysis[D].China University of Petroleum,2009)

7 Liu C R,et al.New model for estimating the thickness and per mittivity of subsurface layers from GPR data[J]. IEE Proe.-Rodur SonorNaviz,2002,149(6)：315-319.

A SPECTRAL INVERSI ON ALGORITHM FOR GPR SIGNALS OF THIN LAYERS

Huang Zhonglai1),Zhang Jianzhong2)and Huang Jilin1)

(1)School of Infor m ation Science and Technology,X iam en University,X iam en 361005 2)College of M arine Geo-science,Ocean University of China,Q ingdao 266100)

In order to extract infor mation including thin layer’s position and thickness,from reflected wave data of Ground Penetrating Radar(GPR)in multi-layered situation,a spectral inversion algorithm realized in frequency domain is presented.Through analysis of electromagnetic wave propagation in multi-layermedia,generalized reflection coefficients are derived.After dividing the coefficients into pairs followed by decomposing the pairs into even and odd components,a cost function is established to connect reflection coefficients’spectrum and layer’sparameters before inversion is carried out.Some simulations show that accurate inversion results can still be expected when thickness of the layer is less than tuning thickness,which therefore i mproves resolution of thin layers.

Ground Penetrating Radar;thin layer(GPR);spectral inversion algorithm;general reflection coefficients;cost function

waves ofmodelA

1671-5942(2011)04-0154-06

2011-03-16

國家自然科學(xué)基金(40774065)

黃忠來,博士研究生,主要研究方向?yàn)樘降乩走_(dá)及隨機(jī)信號(hào)處理.E-mail：huangzhonglai@163.com

P631.2

A