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

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

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

0 引言

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

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

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

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

1 全波形反演方法原理

1.1 探地雷達正演

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

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

(1)

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

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

1.2 探地雷達波形反演

1.2.1 目標函數

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

(2)

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

mk+1=mk+αkpk

(3)

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

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

(4)

1.2.2 局部優化算法

局部優化算法分為梯度類算法和牛頓類算法，梯度法具有工作量少，所需存儲變量少等優點，但是收斂速度慢，計算效率低，易陷入局部極小點，而牛頓法具有精度高，收斂速度快等優點，但是所需存儲變量多，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)構造出近似矩陣，m的取值一般為3～20，本文取為5。

1.2.3 搜索步長

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

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

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

(9)

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

(10)

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

(11)

(12)

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

(13)

(14)

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

1.2.4 反演策略

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

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

1.2.5 時間復雜度分析

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

2 數值試算

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

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

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

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

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

圖4 目標函數Fig.4 Objective function(a)相對介電常數；(b)電導率

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

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

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

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

圖7 目標函數Fig.7 Objective function

3 結論與展望

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