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

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

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

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

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

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

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

1 方法原理

1.1 探地雷達正演

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

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

1.2 傳統(tǒng)目標函數(shù)全波形反演

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

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

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

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

式中：δ為變量的增量。

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

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

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

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

1.3 褶積型目標函數(shù)全波形反演

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

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

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

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

式中：上標cal 為相應(yīng)的計算值。

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

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

2 合成數(shù)據(jù)實驗

2.1 路面塌陷異常褶積型目標函數(shù)與傳統(tǒng)目標函數(shù)FWI 對比

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

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

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

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

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

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

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

圖3 路面塌陷異常傳統(tǒng)目標函數(shù)和褶積型目標函數(shù)FWI結(jié)果對比Fig.3 Comparison between results from conventional objective function and convolution type objective function

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖7 褶積型目標函數(shù)FWI 結(jié)果Fig.7 Convolution type objective function FWI result

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

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

圖8 地下復(fù)雜介質(zhì)模型重構(gòu)誤差曲線Fig.8 Error curve of reconstruction for underground complex medium models

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

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

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

3 實測數(shù)據(jù)全波形反演試驗

3.1 地面塌陷實測數(shù)據(jù)

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

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

由圖10 可知，直達波同相軸具有明顯的斷裂，深部存在一個能量不強的異常反射?？梢耘袛嘣摱蔚缆窋?shù)據(jù)存在多個異常，且路面存在塌陷部分。剖面的末端淺部存在大量的異常波形，但由于同相軸沒有發(fā)生斷裂現(xiàn)象，可以推測異常不在地表，可能為地下淺埋管線。

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

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

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

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

3.2 地下異常實測數(shù)據(jù)

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

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

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

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

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

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

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

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

4 結(jié) 論

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

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

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

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