黃忠來 張建中 黃吉林

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

探地雷達薄層信號的譜反演算法＊

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

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

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

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

1 引言

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

2 譜反演算法過程

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

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

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

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

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

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

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

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

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

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

以包含 4個反射系數(shù)的模型為例,圖 3中反射系數(shù) r1～r4分別為 0.1、0.3、-0.2、0.4,t1～t4分別是 4個反射系數(shù)的距離地面的時間位置,T1、T2是前后兩對反射系數(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ù)對及其奇偶分解Fig.3 Coefficients pairs and their even-odd decompositions

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

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

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

3 廣義反射系數(shù)

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

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

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

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

4 反演算法

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

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

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

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

5.1 楔形模型仿真

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

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

圖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è)計模型在地下有 4個反射界面,反射系數(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個層,上方層的厚度設(shè)計成電磁波在層中波長的 1/10,下方層的厚度設(shè)計成波長的 1/2,即上方是個相對薄層,而下方是相對厚層。

表1 實際值與反演結(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é)果與實際值比較Fig.10 Comparison between inversion results and true values

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

6 結(jié)語

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

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)識別方法——路面雷達關(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].石油地球物理勘探,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)

黃忠來,博士研究生,主要研究方向為探地雷達及隨機信號處理.E-mail：huangzhonglai@163.com

P631.2

A