基于探地雷達的路基地下異常體全波形反演

胡磊磊，陳 康，黃德軍，楊 榮

(1.廣西壯族自治區地球物理勘察院，廣西 柳州 545005；2.廣西壯族自治區第七地質隊，廣西 柳州 545005)

探地雷達(Ground Penetrating Radar,GPR)根據地下介質的電磁性質差異來查找目標異常，是一種高分辨率、無損的淺部地球物理探測方法，在工程質量檢測[1]、公路質量檢測[2]、地質災害研究[3]等領域均有廣泛應用。

在道路病害研究中，識別并提取GPR 剖面中的道路病害特征一直是GPR 能否成功應用的關鍵。在實際工作中，往往存在不同程度的環境干擾，導致檢測精度低、測量效果變差。對此，差值檢測是一種有效提高數據質量的方法[4]。建立地下病害屬性劃分標準，則主要依據回波波組形態、振幅和相位特性、吸收衰減特性等。在理論研究方面，主要通過正演研究電磁信號的傳播特性，模擬得到典型病害體的特征圖譜來對實際情況進行指導[5-6]；在工程上，研究探測的重點區域及其測線布置原則，采集參數設置及數據采集質量的評判方法。同時，對正常道路、典型干擾源和典型道路隱性病害進行標記，并對應其地球物理特征和GPR 波組特征，可以很好地對相似特征下的病害進行解釋，也可為路面養護維修提供決策支持[7-9]。除病害外，道路中瀝青層覆蓋的混凝土裂縫也是GPR 應用于道路無損檢測的重點方向[10]。同時，對道路采集得到的GPR 數據進行特定的處理，能夠有效突出病害特征，再將正演或標記的病害特征與處理后的數據進行綜合對比分析[11]，提高路基中隱藏病害的探明精度和后續解釋的準確度[12]。

反演是一種對所研究區域進行成像分析的有效手段，可以更加直觀地展示病害類型，為病害解釋提供理論依據。在眾多GPR 反演技術中，基于非線性理論基礎的全波形反演(Full Waveform Inversion，FWI)是目前最常用的一種方法。該方法有效利用了全波場信息，通過波動方程約束條件來解決非線性優化問題，經過多次反演迭代得到貼合地下實際情況的物性參數分布，如幾何信息特征和介質的介電特性[13]；在FWI 的應用中，激勵源子波對反演的結果具有很大影響，使用不同的源子波進行FWI，可以驗證不同飽和度的土壤-含水層系統[14]；稀疏盲反卷積技術是一種優化子波的可靠技術手段，建立地下反射率序列的稀疏表示，提高了反演結果的準確率[15]；通過分析測量和模型軌跡的擬合以及最終模型的剩余梯度來對獨立測量的測井實測數據進行驗證分析得出，子波估計主要通過影響初始速度模型進而對FWI 結果造成影響，初始速度模型與目標函數聯系緊密[16]；在實際工程實踐中儀器所產生的激勵源子波估計往往并不準確，若能在建立目標函數的過程中，消除子波影響，可以有效提高FWI 的準確性和實用性[17-18]。

傳統目標函數FWI 算法都是基于激勵源子波已知的假設條件，然而在實際勘探中非常難以獲取準確的激勵源子波形態，激勵源子波估計不準確會降低反演結果的可靠性和準確性。因此，如何降低或消除子波對反演結果的影響至關重要。本文在前人的基礎上，將褶積型目標函數的FWI 表達形式應用于GPR 反演中，在建立目標函數過程中，消除激勵源子波估計不準確的影響，提高反演精度，以期為實際工程應用提供理論依據和數據基礎。

1 方法原理

1.1 探地雷達正演

二維橫磁(Transverse Magnetic,TM)波Maxwell方程如下：

式中：H為磁場強度，A/m；E為電場強度，V/m；ε為介質的介電常數，F/m；μ為介質的磁導率，H/m；σ為電導率，S/m；t為時間，s；J為激勵源；x，y，z為方向。本文選用時域有限差分法(Finite Difference Time Domain,FDTD)算法對GPR 數據進行正演模擬[19]，選擇卷積完全匹配層(Convolutional Perfectly Matched Layer,CPML)吸收邊界用來控制有限模擬區域[20]。

1.2 傳統目標函數全波形反演

FWI 就是建立一個目標函數，通過最優化策略找到某個參數值，使得目標函數最小。對于合成數據，由于子波的形態是已知的，因此，可以假設實際采集到的探地雷達數據為Ez分量，對其建立如下的極小化目標泛函[21]：

式中：Ez為模擬的探地雷達數據(采用猜測參數的正演模擬數據)；tT為記錄總時長；上標obs 為實際的觀測數據(實際測量數據、合成剖面數據)；模擬波動方程必須滿足Maxwell 方程，采用拉格朗日乘子法將無約束優化問題轉化為約束優化問題，目標泛函變為：

式中：[ ψx,ψy,?]為 拉格朗日乘子函數；[e1,e2,e3]為伴隨方程，其表達式為：

利用分部積分法，注意到初始條件、終止條件為零，且滿足自由邊界條件。因此，得到：

式中：δ為變量的增量。

式中：xr為x方向源點位置。

令式(5)、式(6)、式(7)等于0，得到伴隨方程表達式為：

伴隨方程式(8)的形式與式(4)不同。最后得到相應的梯度公式為：

由于本文使用擬牛頓算法，因此，不需要直接求取Hessian 矩陣。

1.3 褶積型目標函數全波形反演

本文采用建立褶積型目標函數的方法，以消除激勵源子波的影響，與傳統目標函數不同，褶積型目標泛函為：

式中：下標ref 為相應的參考道，使用Green 函數表示方法，波動方程可以簡寫為：

式中：d為波動方程表達式，G與f為拆分的矩陣表達式。

將式(11)代入式(10)，得：

式中：上標cal 為相應的計算值。

從化簡后的式(12)中可以分析得到，構成目標函數殘差項的2 個數據項可以看作由相同激勵源信號激發產生的，消除了子波估計不準確對于反演的影響。對式(10)求取梯度，有：

式中：?為相關運算。此外，得到的梯度與式(9)相同。

2 合成數據實驗

2.1 路面塌陷異常褶積型目標函數與傳統目標函數FWI 對比

為了對比2 種不同目標函數對于FWI 結果的影響并說明本文算法的正確性，建立如圖1 所示的地面塌陷模型。塌陷模型分為三層，分別為空氣層(εr=1)，混凝土層(εr=9)及基巖層(εr=6)，路面塌陷分別為無填充(εr=1)及沙土填充(εr=3)，并且在地下具有一個異常空洞體(εr=1)。

圖1 路面塌陷異常相對介電常數分布Fig.1 Relative dielectric constant distribution in case of road collapse exceptions

實驗所用激勵源為400 MHz 的雷克子波，采用屏蔽天線B-Scan 的觀測方式，得到如圖2 所示的合成數據正演剖面。

圖2 路面塌陷異常正演剖面Fig.2 Forward section of road collapse exceptions

由圖2 可以看出，路面的2 個塌陷異常反映在直達波上，無填充的塌陷部分在直達波部分有明顯的下移，并且不存在反射波，但在角點存在繞射異常；沙土填充層的直達波連續，但是由于反射系數的改變導致直達波能量降低，并且具有很明顯的反射回波，且角點也存在繞射波；路面內部的空洞位置出現了一個明顯的雙曲線反射波。約2 ns 處的路面模型內部具有十分明顯的分層特征，同時路面內部結構的塌陷更為清晰。然而對于約10 ns 處空洞下方的層卻出現了同相軸的錯斷，原因是使用地面觀測系統時，由于該處上方具有空洞，導致波傳播到層界面的時間與其他位置不同。

為了說明激勵源子波估計不準確對于FWI 結果的影響，將激勵源調整為450 MHz 的雷克子波。采用混凝土層介電常數構成的均一介質當作反演的初始模型，傳統目標函數FWI 的結果如圖3a 所示。改為褶積型目標函數使用同樣的參數對該模型進行反演，得到如圖3b 所示結果。

對比圖3 可知，當激勵源子波估計不準確時，使用傳統目標函數進行反演得到的結果僅僅可以看出模型的大體趨勢，直達波連續無錯斷，也沒有扭曲變形，無法對塌陷和空洞的位置和形態信息進行準確反映，同時，0.6 ns 處的分層也并不清晰，且反演的具體介電常數數值是不準確的，極大地影響了對于下地介質分布的人工判讀，容易造成錯誤的解釋。而使用褶積型目標函數進行反演得到的結果，準確性得到了很大的提升，2 個地面塌陷、空洞和層位信息都得到了較好的恢復，且介電常數數值對應良好，能夠從介電常數的分布上清楚地確定該處的異常體類型。

圖3 路面塌陷異常傳統目標函數和褶積型目標函數FWI結果對比Fig.3 Comparison between results from conventional objective function and convolution type objective function

特別需要注意的是圖3b 中對于路面內部的空洞異常體，反演得到的異常體大小并不對應，并且在空洞異常體上方出現了虛假的高介電常數異常。該現象的原因是地面觀測這種觀測方式導致了對于深層和靠近模型邊界的異常得到的回波信號較少，導致該位置在剖面上的波形并不完整，且反演的多解性也在一定程度上深化了這一現象。若想規避這一現象則需要獲取更多的地下介質回波信息的約束，如鉆孔雷達、聲波結果等。

為了能夠量化說明FWI 算法的精確性，建立如下式所示的重構誤差函數:

式中：p為評價的參數，這里為相對介電常數；p0為初始模型參數；ptrue為真實模型參數；pk為第k次迭代后的模型參數。計算 Δ (p)得到如圖4 所示的模型重構誤差曲線。

由圖4 可知，隨著反演迭代次數的增加，模型重構誤差逐漸減小，最后趨于穩定，同時由于探地雷達觀測方式的限制，信息量不足以將地下介質信息完全恢復，因此，重構誤差會逐漸穩定且不再下降。

圖4 路面塌陷異常FWI 重構誤差曲線Fig.4 Error curve of ground surface collapse model reconstruction

在實際工作中通常還需要考慮算法的效率問題，本文基于Ubuntu 20.04 操作系統，使用Python3.8 編程語言進行全波形反演實驗，CPU 配置為Intel Xeon(R)CPU E5-2690 3.00 GHz × 20，反演過程中使用20 個進程并行的策略，得到FWI 效率對比(表1)。

由表1 可知，由于本文算法需要進行卷積運算，因此，單次迭代耗時略大于傳統算法，FWI 總耗時與單次迭代耗時并不是線性關系，所以傳統目標函數FWI 算法的總耗時反而大于本文提出的算法，這是由于錯誤的激勵源子波估計，造成了整體算法的不穩定，需要多次進行迭代步長的選取，而本文算法由于消除了激勵源子波估計對于反演的影響，因此，整體算法更為穩定，從而獲得了更高的算法效率。

表1 合成數據全波形反演效率對比Table 1 Comparison of full-waveform inversion efficiency for synthetic data

2.2 地下空洞異常褶積型目標函數FWI 實驗

為了說明基于褶積型目標函數FWI 對于地下復雜介質的反演效果和算法的正確性，建立如圖5 所示的地下模型。該模型分為兩層，分別為地表層(εr=4)及基巖層εr=9)，并且在地表層設置了6 個異常體，分別用于模擬地下空洞(εr=1)和地下管線(εr=6,εr=8,εr=9,εr=10)。

圖5 地下空洞異常相對介電常數分布Fig.5 Relative dielectric constant distribution in case of underground void exceptions

采用400 MHz 的雷克子波作為激勵源，采用屏蔽天線B Scan 的觀測方式，得到如圖6 所示的合成數據正演剖面。

圖6 路面塌陷異常正演剖面Fig.6 Forward section of road collapse exceptions

由圖6 可知，模型正演剖面十分復雜，由于異常體的相對介電常數各不相同，因此，剖面中的反射波能量也有所差異，具體表現為：相對介電常數相差越大，得到的反射波能量也越強。同時需要注意的是，由于管線存在垂直方向的重疊，正演剖面中的反射回波雙曲線存在耦合。

為了說明激勵源子波估計不準確對于FWI 結果的影響，采用600 MHz 的雷克子波作為反演所使用的激勵源，采用地表層介電常數構成的均一介質當作全波形反演的初始模型，使用褶積型目標函數，得到如圖7 所示的反演結果。

圖7 褶積型目標函數FWI 結果Fig.7 Convolution type objective function FWI result

對比圖5 和圖7 可知，使用褶積型目標函數FWI得到的結果由于消除了激勵源子波估計不準確帶來的影響，得到了較好的反演結果，較為明顯地區分了地層，6 個異常體位置及數值也得到了較好的對應。同時需要注意的是圖7 中反演得到的異常體大小并不對應，同時空洞異常體上方出現了虛假的高介電常數異常，造成這個現象的原因是觀測系統采用地面觀測，這種觀測方式導致了對于深層和靠近模型邊界的異常得到的回波信號較少，導致該位置在剖面上的波形并不完整，并且由于全波形算法本身存在的多解性，從而導致了這一現象。同時，對于距離較近的不同相對介電常數的異常體也能夠通過反演算法進行良好的對應，進而說明了基于褶積型目標函數的FWI 算法對于復雜地下介質的正確性。

為了能夠量化說明算法的精確性和效率，使用式(16)所示的重構誤差函數，得到如圖8 所示的模型重構誤差曲線，并且通過記錄得到反演效率(表2)。

圖8 地下復雜介質模型重構誤差曲線Fig.8 Error curve of reconstruction for underground complex medium models

表2 合成數據全波形反演效率Table 2 Full-waveform inversion efficiency for synthetic data

由圖8 可知，隨著反演迭代次數的增加，模型重構誤差逐漸減小，且下降幅度逐漸變緩，同時由于探地雷達觀測方式的限制，信息量不足以將地下介質信息完全恢復。

由表2 可知，由于使用了更細的模型網格，更深的測量深度增加了記錄時長，因此，增加了單次迭代耗時，同時算法的總耗時隨之增加。通過分析得出，對于更為精細的反演模型需要采用更細的網格，但同時造成了反演算法效率的下降，因此，需要根據實際的反演目標選擇合適的網格大小。

3 實測數據全波形反演試驗

3.1 地面塌陷實測數據

為了驗證本文算法的實用性，采用實測的GPR 剖面數據進行算法試驗。本次試驗路面數據采集自某城市路面塌陷項目，儀器選用MALA 探地雷達儀器，測量區域表面為平整的混凝土路面，周圍道路暫時封閉，無其他人為干擾。圖9 為現場實際工作記錄圖，其中紅色方框為實際的未填充的塌陷區域，黑色有向線段表示選取測線的實際走向。此次探測的測量區域深度為2.00 m，測線長度為1.00 m，使用400 MHz 的屏蔽天線，時窗長度為40 ns，采用距離測量模式獲取589道數據，每道的采樣點數設置為512。

圖9 實測路面現場工作記錄Fig.9 Field work record for measured pavement

由圖10 可知，直達波同相軸具有明顯的斷裂，深部存在一個能量不強的異常反射。可以判斷該段道路數據存在多個異常，且路面存在塌陷部分。剖面的末端淺部存在大量的異常波形，但由于同相軸沒有發生斷裂現象，可以推測異常不在地表，可能為地下淺埋管線。

圖10 實測路面的GPR 剖面Fig.10 GPR profile of measured pavement

使用本文的褶積型目標函數FWI 算法對該數據進行反演計算，采用均一模型作為反演的初始模型，激勵源選用400 MHz 雷克子波，得到的反演結果如圖11所示。

圖11 實測數據全波形反演介電常數分布Fig.11 Dielectric constant distribution in full-waveform inversion for measured data

由圖11 可知，傳統目標函數的反演結果中沒有任何有效信息，這是由于實際工作中的儀器激勵源子波未知，在反演的過程中所使用的子波估計所迭代的結果中的錯誤累計會逐漸增大，最終導致這一結果。而褶積型目標函數的反演結果反映了該路面不僅存在路面塌陷，在路面塌陷下方存在明顯的充填不密實，并且在路面深部同樣存在高介電常數的方形異常體，根據相對介電常數的分布可推測為規則排列的金屬管線。同樣，對于剖面末端的淺層異常，根據相對介電常數的分布，可推測為淺埋的金屬電纜，并且通過查閱相關資料，最終確定了該處存在2 處不同埋深的管線分布。

3.2 地下異常實測數據

選擇MALA 探地雷達儀器進行GPR 實際探測，測量區域表面為干燥的沙土路面，采用測線模式采集數據。已知地下埋藏未知深度的管線，此次探測的測量區域深度為4.00 m，測線長度為8.00 m，使用600 MHz 的屏蔽天線，時窗長度為60 ns，采用距離測量模式獲取682 道數據，每道的采樣點數設置為512。采集到的實測數據如圖12 所示。

圖12 實測路面GPR 剖面Fig.12 GPR profile of measured pavement

由圖12 可知，該數據表面存在層異常，并且能夠看到3 個反射異常。使用本文的褶積型目標函數FWI 算法對該數據進行反演計算，采用均一模型作為反演的初始模型，采用600 MHz 雷克子波作為反演的激勵源，得到的反演結果如圖13 所示

圖13 實測數據全波形反演介電常數分布Fig.13 Dielectric constant distribution in full-waveform inversion for measured data

由圖13 可知，顯示了3 個管線異常的相對位置及埋深情況，且層位信息得到了較好的恢復，并且從反演得到的相對介電常數分布，可以推測出該管線為非金屬管線，初步判斷為PVC 材質的小管徑管線。

同時，為了說明本文提出的算法針對實測數據的反演效率，記錄了模型參數和反演效率(表3)。

表3 實測數據全波形反演效率Table 3 Full-waveform inversion efficiency for measured data

通過對實測數據的FWI 試驗，證明了本文算法對于實測數據的適用性，同時，由于無法對實測數據激勵源子波進行良好的估計，因此，采用傳統目標函數FWI 將對反演結果造成較大的影響，但是本文提出的褶積型目標函數則能夠有效地避免這一問題，并且記錄了相應的算法耗時，結果表明該算法效率雖然不能做到實時處理，但針對復雜且重要的未知地下異常分布的分析依然能夠進行指導，為生產實際提供有效的依據。

4 結 論

a.針對城市道路塌陷隱患探測問題，提出探地雷達適用的褶積型目標函數的全波形反演算法。

b.該算法能夠克服常規反演目標函數由于激勵源子波估計不準確對全波形反演結果準確性造成的影響，取得良好的反演結果；在計算效率和反演準確度上，均優于常規的全波形反演。

c.將該算法應用于GPR 實測數據反演中，不僅較好識別異常體位置及層位信息，且可從反演的相對介電常數中推測出異常體材質，驗證了本文方法的實用性。

d.該算法可以得到較為準確的相對介電常數分布，與 GPR 剖面相互驗證，提高反演精度，為實際工程應用提供理論依據和數據基礎。