俞海龍， 馮 晅,2， 恩和得力海， 趙建宇， 孫成城

(1.吉林大學(xué) 地球探測科學(xué)與技術(shù)學(xué)院，長春 130026;2.地球信息探測儀器教育部重點(diǎn)實(shí)驗(yàn)室(吉林大學(xué))，長春 130026)

0 引言

探地雷達(dá)(Ground Penetrating Radar)，是一種利用高頻無線電波來確定地下物質(zhì)分布規(guī)律的地球物理方法[1]，它具有采樣率高、精度高、抗干擾性強(qiáng)、無損探測等優(yōu)點(diǎn)。由于探地雷達(dá)具有能夠?qū)\層地下物質(zhì)高分辨率成像的特點(diǎn)，因此，它在工程地球物理勘探，環(huán)境研究以及基礎(chǔ)設(shè)施研究中都被廣泛地應(yīng)用[2]。

探地雷達(dá)波形反演在最優(yōu)化理論的框架下，充分利用雷達(dá)記錄的波形、走時(shí)、相位等信息，通過擬合觀測數(shù)據(jù)與模擬數(shù)據(jù)，使殘差達(dá)到最小，來反演地下介質(zhì)的物性參數(shù)(如電導(dǎo)率、介電常數(shù)、磁導(dǎo)率等)的方法，而電導(dǎo)率σ和介電常數(shù)ε最能夠反映出地下介質(zhì)的電磁學(xué)性質(zhì)[3]，因此反演電導(dǎo)率和介電常數(shù)具有著重要的意義。

大量的全波形反演方法包括聲波、彈性波、粘彈性波、各向異性聲波方程等的有限差分或有限元法分別在時(shí)間域[4-6]和頻率域[7-8]中研究并且應(yīng)用于地震數(shù)據(jù)，但是到目前為止，探地雷達(dá)反演的發(fā)展仍具有一定的局限性，這是因?yàn)槔走_(dá)信號的振幅和相位依賴于天線的方向以及波場的傳播路徑，所以我們必須了解天線的輻射方向以及雷達(dá)數(shù)據(jù)的矢量特性。對于跨孔雷達(dá)，Ernst等[9]、Kuroda等[10]提出了基于二維麥克斯韋方程的時(shí)域有限差分全波形反演方法；Meles 等[11]考慮了全波形反演電參數(shù)的矢量特性，在反演中允許電導(dǎo)率與介電常數(shù)同時(shí)更新；Kalogeropoulos等[12]將全波形反演應(yīng)用于評價(jià)混凝土中的濕度和氯化物；王兆磊等[13]利用探地雷達(dá)全波形反演方法反演出二維介質(zhì)電導(dǎo)率和介電常數(shù)； Lambot等[14]、Klotzsche等[15]、Streich等[16]分別研究了探地雷達(dá)波形反演。

筆者首先對TM模式下的麥克斯韋方程進(jìn)行時(shí)域有限差分正演模擬，并采用單軸各向異性(UPML)吸收邊界條件以消除邊界反射的影響[17]，分別進(jìn)行電導(dǎo)率和介電常數(shù)的單參數(shù)反演，以及電導(dǎo)率和介電常數(shù)雙參數(shù)同時(shí)反演，在反演的過程中，由于不同參數(shù)之間的相互耦合、相互干涉的影響，繼而增加了反演的非線性程度，另外不同參數(shù)的物理量綱不同，致使進(jìn)一步增加了反演的病態(tài)性[18-19]，在多參數(shù)的全波形反演中，考慮到近似Hessian矩陣的非對角塊元素反映不同參數(shù)之間的耦合效應(yīng)，因此在反演中利用近似Hessian矩陣的逆來夠消除不同參數(shù)之間的影響[20-21]，進(jìn)而提高目標(biāo)函數(shù)的收斂速度，獲得準(zhǔn)確的反演結(jié)果，采用L-BFGS算法進(jìn)行探地雷達(dá)全波形反演，并和梯度法做了對比，結(jié)果顯示L-BFGS算法能夠更好的消除電導(dǎo)率和介電常數(shù)之間的影響。

1 全波形反演方法原理

1.1 探地雷達(dá)正演

眾所周知，正問題是相對于反問題而言的，正演模擬的精度嚴(yán)重影響著反演的結(jié)果，筆者采用TM模式下的麥克斯韋方程組，利用高精度空間10階進(jìn)行波場模擬。

TM模式下麥克斯韋旋度方程[17]為式(1)。

(1)

其中：ε為相對介電常數(shù)；μ為相對磁導(dǎo)率；σ為電導(dǎo)率；Ez為電場強(qiáng)度z分量；Hx、Hy分別為磁場強(qiáng)度x、y分量。

由于實(shí)際計(jì)算區(qū)域是有限的，當(dāng)波場傳到邊界時(shí)會發(fā)生邊界反射，為消除邊界反射的影響，采用吸收效果較好的單軸各向異性(UPML)吸收邊界條件。

1.2 探地雷達(dá)波形反演

1.2.1 目標(biāo)函數(shù)

全波形反演通過擬合正演數(shù)據(jù)和實(shí)測數(shù)據(jù)，使兩者之差的平方達(dá)到最小，目標(biāo)函數(shù)可以寫為式(2)。

(2)

全波形反演是大規(guī)模的無約束非線性反演問題，一般采用局部優(yōu)化方法來求解，模型參數(shù)的迭代更新公式可表示為式(3)。

mk+1=mk+αkpk

(3)

其中：m為待反演的參數(shù)，本文中為介電常數(shù)ε和電導(dǎo)率σ；αk為迭代步長；pk為迭代方向；k為迭代次數(shù)。

目標(biāo)函數(shù)關(guān)于模型參數(shù)的梯度[11]為式(4)。

(4)

1.2.2 局部優(yōu)化算法

局部優(yōu)化算法分為梯度類算法和牛頓類算法，梯度法具有工作量少，所需存儲變量少等優(yōu)點(diǎn)，但是收斂速度慢，計(jì)算效率低，易陷入局部極小點(diǎn)，而牛頓法具有精度高，收斂速度快等優(yōu)點(diǎn)，但是所需存儲變量多，Hessian矩陣以及Hessian矩陣的逆求取困難，為了避免Hessian矩陣的直接求取，進(jìn)而提出了擬牛頓算法，其中L-BFGS算法是一種公認(rèn)的最有效的擬牛頓類算法[22]。

梯度法中模型迭代公式為式(5)。

εk+1=εk-αkε▽Sε

σk+1=σk-αkσ▽Sσ

(5)

擬牛頓法中模型迭代公式為式(6)。

(6)

(7)

對式(7)進(jìn)行m次遞推可得式(8)。

8)

因此，只需記錄m個(gè)向量{pi,qi}，i=k-m、…、k-1，就可以按式(8)構(gòu)造出近似矩陣，m的取值一般為3～20，本文取為5。

1.2.3 搜索步長

搜索步長的選取對反演結(jié)果有著重要的影響，搜索步長太小，需要增加迭代次數(shù)以使得目標(biāo)函數(shù)收斂，進(jìn)而增加了計(jì)算量；搜索步長太大，可能會破壞迭代過程的穩(wěn)定性，且極其容易陷入局部極小值。當(dāng)前常用的步長搜索方法有線性搜索方法以及拋物線擬合法等，線性搜索方法雖然簡單，但是對于高維數(shù)問題的求解耗時(shí)較多；拋物線法在每一次迭代中都需進(jìn)行多次正演以得到最優(yōu)步長，而全波形反演需要多次迭代才能夠得到最終結(jié)果，因此拋物線法的計(jì)算量也是相當(dāng)大的。

筆者采用A. Pica*等[23]提出的步長求解方法，根據(jù)式(2)全波形反演目標(biāo)函數(shù)為式(9)。

S(σk+1,εk+1)=S(σk+αk▽Skσ,εk+αk▽Skε)=

(9)

對目標(biāo)函數(shù)S(σk+αk▽Skσ,εk+αk▽Skε)求導(dǎo)，令導(dǎo)數(shù)為零，以求取使得目標(biāo)函數(shù)達(dá)到最小時(shí)所對應(yīng)的步長為式(10)。

(10)

最終得到步長的表達(dá)式為式(11)。

(11)

(12)

式(11)中介電常數(shù)和電導(dǎo)率采用相同的步長，Meles等[11]通過試驗(yàn)說明式(11)求取的步長主要取決于介電常數(shù)的梯度，而電導(dǎo)率的梯度貢獻(xiàn)很小，幾乎不起作用，對此，分別求取電導(dǎo)率和介電常數(shù)的步長，電導(dǎo)率和介電常數(shù)的步長求取如下

(13)

(14)

其中：穩(wěn)定因子δσ、δε的選取需滿足式(12)，從式(13)、式(14)可以看出，用此種方法求取步長只需做兩次正演模擬即可，相比較拋物線法可以大大減少正演次數(shù)，進(jìn)而提高計(jì)算效率。

1.2.4 反演策略

反演主要分三步：①對介電常數(shù)單獨(dú)反演，即假定電導(dǎo)率已知，介電常數(shù)為未知參數(shù)；②對電導(dǎo)率單獨(dú)反演，即假定介電常數(shù)已知，電導(dǎo)率為未知參數(shù)；③對介電常數(shù)和電導(dǎo)率同時(shí)反演，即電導(dǎo)率和介電常數(shù)都為未知參數(shù)。在雙參數(shù)同時(shí)反演時(shí)，分別采用梯度法和L-BFGS法，并就兩者反演的精度和計(jì)算效率進(jìn)行了對比。圖1是全波形反演計(jì)算流程，該計(jì)算流程是對介電常數(shù)和電導(dǎo)率同時(shí)反演，對于單參數(shù)反演，只需計(jì)算相應(yīng)參數(shù)的梯度、步長以及模型更新量即可。

圖1 探地雷達(dá)全波形反演計(jì)算流程Fig.1 Flow chat of GPR full waveform inversion

1.2.5 時(shí)間復(fù)雜度分析

模型的空間采樣間隔為0.01 m，網(wǎng)格數(shù)為50*50，時(shí)間采樣間隔為0.02 ns，F(xiàn)DTD時(shí)間迭代次數(shù)為1 000次，共有9個(gè)場源位置，采用單精度(4 bit)保存，電場值保存需要內(nèi)存為50*50*1000*9/1024/1024=21.457 7 M，計(jì)算反傳波場同樣需要內(nèi)存21.457 7 M，共計(jì)42.915 4 M內(nèi)存。在一次迭代過程中每一個(gè)源都需要4次正演計(jì)算，第一次得到正演數(shù)據(jù)，第二次波場反傳計(jì)算梯度，另外兩次得到迭代步長，因?yàn)楣灿?個(gè)場源，因此每一次迭代共需要36次正演計(jì)算。

2 數(shù)值試算

筆者設(shè)計(jì)的試算模型為方格模型(圖(2))，圖2(a)中，背景模型的相對介電常數(shù)為6，左右兩個(gè)方形異常體的相對介電常數(shù)分別為5和7。圖2(b)中，背景模型的電導(dǎo)率為0.05 S·m-1，左右兩個(gè)方形異常體的電導(dǎo)率分別為0.01 S·m-1和0.1 S·m-1。模型大小為0.5 m×0.5 m，網(wǎng)格大小為50×50，網(wǎng)格間距為0.01 m，采用一點(diǎn)激發(fā)多點(diǎn)接收的采集方式，發(fā)射天線位于垂直位置為0 m處，發(fā)射天線初始點(diǎn)位于第5個(gè)網(wǎng)格點(diǎn)，發(fā)射天線之間間隔為5個(gè)網(wǎng)格點(diǎn)，共有9個(gè)發(fā)射天線，接收天線位于垂直位置為0 m以及水平位置為0 m和0.5 m處，每個(gè)位置均有50個(gè)接收點(diǎn)，時(shí)間采樣間隔為0.02 ns，記錄時(shí)間為12 ns，場源函數(shù)選用主頻為800 MHz的雷克子波。

1)介電常數(shù)單獨(dú)反演，介電常數(shù)真實(shí)模型為圖2(a)，電導(dǎo)率為0.05 S·m-1，觀測數(shù)據(jù)是由真實(shí)模型正演得到的，初始模型為相對介電常數(shù)為6，電導(dǎo)率為0.05 S·m-1的均勻半空間，反演采用L-BFGS算法，迭代次數(shù)為73次，圖2(c)是介電常數(shù)反演結(jié)果，圖3(a)是在深度為0.25 m處穿過兩異常體中心的橫向剖面，從圖2(c)和圖3(a)中可以看出，異常體的形態(tài)能夠非常準(zhǔn)確地刻畫出來，并且在數(shù)值上，反演結(jié)果和真實(shí)模型也非常接近。圖4(a)為目標(biāo)函數(shù)隨迭代次數(shù)變化的收斂曲線，從圖4(a)中可以看出，采用L-BFGS算法，目標(biāo)函數(shù)收斂速度很快，迭代35次即可收斂。

2)電導(dǎo)率單獨(dú)反演，電導(dǎo)率真實(shí)模型為圖2(b)，相對介電常數(shù)為6，觀測數(shù)據(jù)是由真實(shí)模型正演得到的，初始模型為電導(dǎo)率為0.05 S·m-1，相對介電常數(shù)為6的均勻半空間，反演采用L-BFGS算法，迭代次數(shù)為100次，圖2(d)是電導(dǎo)率反演結(jié)果，圖3(b)是在深度為0.25 m處穿過兩異常體中心的橫向剖面，從圖2(d)和圖3(b)可以看出，無論在異常體的形態(tài)上還是在數(shù)值上，反演結(jié)果都接近于真實(shí)模型。圖4(b)為目標(biāo)函數(shù)隨迭代次數(shù)變化的收斂曲線，從圖中可以看出，采用L-BFGS算法，目標(biāo)函數(shù)收斂速度很快，迭代40次即可收斂。

圖2 梯度法單參數(shù)反演結(jié)果Fig.2 Radient method of single parameter of the full waveform inversion results(a)真實(shí)相對介電常數(shù)；(b)真實(shí)電導(dǎo)率；(c)單參數(shù)介電常數(shù)反演結(jié)果；(d)單參數(shù)電導(dǎo)率反演結(jié)果

圖3 穿過兩異常體中心的單參數(shù)反演結(jié)果橫向剖面Fig.3 The transverse profile of through the two abnormal body center of single parameter inversion results(a)相對介電常數(shù)；(b)電導(dǎo)率

圖4 目標(biāo)函數(shù)Fig.4 Objective function(a)相對介電常數(shù)；(b)電導(dǎo)率

圖5 梯度法雙參數(shù)全波形反演結(jié)果Fig.5 Gradient method of double parameters of the full waveform inversion results(a)相對介電常數(shù)；(b)電導(dǎo)率,L-BFGS雙參數(shù)反演結(jié)果；(c)相對介電常數(shù)；(d)電導(dǎo)率

圖6 穿過兩異常體中心的雙參數(shù)反演結(jié)果橫向剖面Fig.6 Through the center of the two abnormal body double parameters inversion results transverse section(a)相對介電常數(shù)；(b)電導(dǎo)率

3)介電常數(shù)與電導(dǎo)率同時(shí)反演，真實(shí)模型分別為圖2(a)、圖2(b)，初始模型為相對介電常數(shù)為6、電導(dǎo)率為0.05 S·m-1的均勻半空間，分別采用梯度法和L-BFGS算法進(jìn)行反演，圖5(a)、圖5(b)分別是梯度法的介電常數(shù)和電導(dǎo)率的反演結(jié)果，圖5(c)、圖5(d)分別是L-BFGS算法的介電常數(shù)和電導(dǎo)率反演結(jié)果，圖6分別是介電常數(shù)(圖6(a))和電導(dǎo)率(圖6(b))在深度為0.25 m處穿過兩異常體中心的橫向剖面。

從圖6中可以看出，L-BFGS算法的反演結(jié)果比梯度法的反演的結(jié)果，在異常體的形態(tài)上相差不大，但是在數(shù)值上，L-BFGS算法更加接近于真實(shí)模型，并且L-BFGS算法的收斂速度和精度都高于梯度法，這是因?yàn)樵诙鄥?shù)全波形反演中，Hessian矩陣的非對角塊元素反應(yīng)不同參數(shù)之間的耦合效應(yīng)，而L-BFGS算法是一種擬牛頓算法，即構(gòu)造出一個(gè)不含二階導(dǎo)數(shù)的矩陣來近似Hessian矩陣，可以有效地解決多參數(shù)反演中不同參數(shù)之間的影響。圖7為目標(biāo)函數(shù)隨迭代次數(shù)變化的收斂曲線，從圖7中可以看出，L-BFGS算法比梯度法的目標(biāo)函數(shù)收斂速度更快，目標(biāo)函數(shù)下降的更小。

圖7 目標(biāo)函數(shù)Fig.7 Objective function

3 結(jié)論與展望

從TM模式下麥克斯韋方程組出發(fā)，采用時(shí)域有限差分的方法進(jìn)行雷達(dá)波場的數(shù)值模擬，并采用單軸各向異性(UPML)吸收邊界條件以消除邊界反射的影響，對于探地雷達(dá)波形反演，給出了電導(dǎo)率和介電常數(shù)梯度方向的計(jì)算方法，并以步長為自變量通過求取目標(biāo)函數(shù)為極值的方式來確定最優(yōu)步長，分別對介電常數(shù)、電導(dǎo)率進(jìn)行單參數(shù)反演，另外也對介電常數(shù)與電導(dǎo)率雙參數(shù)同時(shí)反演，并利用擬牛頓算法L-BFGS法來消除不同參數(shù)之間的耦合影響，從反演的結(jié)果可以看出，對于單參數(shù)反演，無論是電導(dǎo)率還是介電常數(shù)，反演結(jié)果都十分接近于真實(shí)模型，對于雙參數(shù)同時(shí)反演，反演結(jié)果的異常體形態(tài)接近于真實(shí)模型，但是由于電導(dǎo)率和介電常數(shù)之間的耦合影響，致使在數(shù)值上不能像單參數(shù)反演那樣得到準(zhǔn)確的反演結(jié)果，但是和真實(shí)模型相差不大。本文波形反演是在時(shí)間域進(jìn)行的，今后將采用并行算法以解決巨量的正演計(jì)算，進(jìn)而提高計(jì)算效率。