基于梯度法和L-BFGS算法的探地雷達時間域全波形反演

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

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

0 引言

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

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

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

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

1 全波形反演方法原理

1.1 探地雷達正演

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

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

(1)

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

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

1.2 探地雷達波形反演

1.2.1 目標函數(shù)

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

(2)

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

mk+1=mk+αkpk

(3)

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

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

(4)

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

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

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

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

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

(5)

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

(6)

(7)

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

8)

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

1.2.3 搜索步長

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

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

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

(9)

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

(10)

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

(11)

(12)

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

(13)

(14)

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

1.2.4 反演策略

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

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

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

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

2 數(shù)值試算

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

1)介電常數(shù)單獨反演，介電常數(shù)真實模型為圖2(a)，電導(dǎo)率為0.05 S·m-1，觀測數(shù)據(jù)是由真實模型正演得到的，初始模型為相對介電常數(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)能夠非常準確地刻畫出來，并且在數(shù)值上，反演結(jié)果和真實模型也非常接近。圖4(a)為目標函數(shù)隨迭代次數(shù)變化的收斂曲線，從圖4(a)中可以看出，采用L-BFGS算法，目標函數(shù)收斂速度很快，迭代35次即可收斂。

2)電導(dǎo)率單獨反演，電導(dǎo)率真實模型為圖2(b)，相對介電常數(shù)為6，觀測數(shù)據(jù)是由真實模型正演得到的，初始模型為電導(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é)果都接近于真實模型。圖4(b)為目標函數(shù)隨迭代次數(shù)變化的收斂曲線，從圖中可以看出，采用L-BFGS算法，目標函數(shù)收斂速度很快，迭代40次即可收斂。

圖2 梯度法單參數(shù)反演結(jié)果Fig.2 Radient method of single parameter of the full waveform inversion results(a)真實相對介電常數(shù)；(b)真實電導(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 目標函數(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)率同時反演，真實模型分別為圖2(a)、圖2(b)，初始模型為相對介電常數(shù)為6、電導(dǎo)率為0.05 S·m-1的均勻半空間，分別采用梯度法和L-BFGS算法進行反演，圖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算法更加接近于真實模型，并且L-BFGS算法的收斂速度和精度都高于梯度法，這是因為在多參數(shù)全波形反演中，Hessian矩陣的非對角塊元素反應(yīng)不同參數(shù)之間的耦合效應(yīng)，而L-BFGS算法是一種擬牛頓算法，即構(gòu)造出一個不含二階導(dǎo)數(shù)的矩陣來近似Hessian矩陣，可以有效地解決多參數(shù)反演中不同參數(shù)之間的影響。圖7為目標函數(shù)隨迭代次數(shù)變化的收斂曲線，從圖7中可以看出，L-BFGS算法比梯度法的目標函數(shù)收斂速度更快，目標函數(shù)下降的更小。

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

3 結(jié)論與展望

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