基于多偏移距探地雷達數(shù)據(jù)的包絡-波形反演

槐 楠 曾昭發(fā) 李 靜 王 卓

①(吉林大學地球探測科學與技術學院 長春 130026)②(吉林建筑大學測繪與勘查工程學院 長春 130118)

1 引言

探地雷達(Ground Penetrating Radar, GPR)是一種用于研究淺地表電磁特性的無損檢測技術，主要基于電磁波傳播原理，通過激發(fā)并記錄高頻電磁波(20～1000 MHz)在地下介質中(地下障礙物或不連續(xù)地質體)的反射、繞射等傳播特性來探測電性結構分布。目前被廣泛應用于土木工程場地勘查[1]、地雷探測[2]、環(huán)境與水文地質研究[3,4]、農業(yè)評估[5,6]、考古調查[7,8]、極地考察[9]以及航空航天[10,11]等前沿工程與科學研究領域。

探地雷達的分辨率和探測深度極易受到脈沖源頻率和地下電導率分布的限制，其中，高頻電磁脈沖具有較高的分辨率，但其在地下介質中的衰減更快，因此穿透深度較小。根據(jù)測量方式的差異，探地雷達可進一步被劃分為共偏移距和多偏移距兩種類型。相較于共偏移距的測量方式，多偏移距測量可以靈活調整發(fā)射和接收天線之間的距離，捕獲到的繞射波信息較為可靠，有助于實現(xiàn)精確的速度分析以及對地下電性參數(shù)的定量估計。因此開展基于多偏移距探地雷達數(shù)據(jù)的研究是十分必要的。

目前，探地雷達數(shù)據(jù)的解釋工作已經從單純推斷目標體空間位置、形態(tài)、尺寸以及層界面等信息，上升到對介質電性參數(shù)(介電常數(shù)、電導率)的反演估計和對近地表結構的精細刻畫。現(xiàn)行的反演方法主要可分為基于射線理論的層析成像方法和基于全波場理論的全波形反演方法。其中，全波形反演(Full-Waveform Inversion, FWI)最初起源于時域地震成像，后被推廣到頻率域，其同時兼顧了波場的運動學和動力學特征(振幅和相位)，是一種高分辨率、高精度的成像方法。全波形反演已在地震勘探領域獲得廣泛應用，但由于地面雷達測量方式對地下介質的照明能力有限，基于地面雷達數(shù)據(jù)的全波形反演具有挑戰(zhàn)性。Lavoué等人[12]采用擬牛頓算法針對多偏移距地面雷達數(shù)據(jù)開展頻率域FWI，同時重構了介電常數(shù)和電導率的2維分布。Feng等人[13,14]將井間地震FWI和GPR FWI集成到一個聯(lián)合程序中，用于對不同地球物理參數(shù)(P波速度、介電常數(shù)和電導率)模型的定量成像；Feng等人[15]將多尺度反演策略和吉洪諾夫正則化(Tikhonov regularization, TV)約束方法結合應用于時域GPR FWI過程，實現(xiàn)了介電常數(shù)和電導率的同步反演；王珣等人[16]在時間域采用改進的全變差正則化策略，開展了GPR雙參數(shù)多尺度同步反演研究；Huai等人[17]在時間域提出了一種基于模型的層剝離FWI，有效增強了模型深部的重構質量和反演分辨率。

然而，當利用梯度類最優(yōu)化算法求解基于Born近似框架下的全波形反演問題時，一個嚴重偏離正解的初始模型會在波形匹配過程中產生“周波跳躍”問題，使得目標函數(shù)傾向陷入到局部極值[18]。因此，全波形反演高度依賴精確的初始模型來確保反演過程的正確收斂。低頻成分對于恢復長波背景速度結構進而構建大尺度初始模型至關重要，但是實際采集的數(shù)據(jù)中所攜帶的低頻信息十分有限，這就使得研究低頻缺失情況下的全波形反演方法顯得尤為重要。最初在地震勘探領域提出并發(fā)展起來的包絡反演是用于重構低頻成分比較知名的一類方法。Bozda?等人[19]最早構建了基于瞬時相位和包絡的全波形反演目標函數(shù)，該種方法旨在利用Hilbert變換來實現(xiàn)波形相位和振幅信息的有效分離，從而降低了反問題的非線性程度。Chi等人[20,21]和Wu等人[22]構建了基于包絡殘差的全波形反演目標函數(shù)，并利用伴隨狀態(tài)法求取了梯度算子，利用包絡反演成功重建了地下介質的大尺度背景信息。Wu等人[23]提出了調制-褶積信號模型并在研究中賦予了包絡數(shù)據(jù)以明確的物理意義。劉新彤等人[24]研究了低頻缺失下基于包絡目標函數(shù)的全波形反演，在有效還原低頻信息的同時提供了對地下介質的定量解釋，然而這一研究僅針對跨孔雷達的觀測模式。與跨孔雷達相比，地面雷達對地下介質的照明角度有限，這本身就增加了反問題的非線性和不確定性；當?shù)孛胬走_缺失低頻成分的情況下，也會極大地制約全波形反演的穩(wěn)定性和精確度。因此，研究基于地面雷達數(shù)據(jù)的包絡全波形反演具有重要意義。

此外，引用分頻處理的多尺度反演策略也有利于緩解全波形反演的“周波跳躍”問題，提高反演的穩(wěn)定性。頻率多尺度方法可按低頻到高頻依次反演多個離散頻率或頻率組，其中，利用低頻記錄對局部極值的不敏感特性來重建平滑的大尺度背景結構，隨后再利用高頻記錄刻畫局部精細結構[25]。Boonyasiriwat等人[26]最初在地震勘探領域通過選擇最優(yōu)頻帶和維納濾波來實現(xiàn)更高效的時間域多尺度全波形反演；隨后，這一策略被成功借鑒到探地雷達的全波形反演中。Meles等人[27]通過在模型更新過程中由低到高逐漸增加頻率成分來擴展數(shù)據(jù)的頻率帶寬，從而實現(xiàn)了對跨孔-井地數(shù)據(jù)的時-頻域聯(lián)合反演。Lavoué等人[12]開展了與頻率采樣有關的多參數(shù)全波形反演研究，通過數(shù)值實驗證實了寬頻帶數(shù)據(jù)的同時反演可以有效增強雙參數(shù)同步重建質量。Li等人[28]提出了一種基于跨孔雷達數(shù)據(jù)的Laplace域波形反演，并提出了頻率跳躍的多尺度波形反演策略，實現(xiàn)了反演過程的快速收斂。

本文在現(xiàn)有的探地雷達時間域全波形反演框架下，回顧了地震勘探中信號包絡提取和常規(guī)包絡反演方法的基本原理與實施流程，并將其成功推廣到探地雷達領域，提出了基于地面多偏移距雷達數(shù)據(jù)的包絡-波形反演方法，實現(xiàn)了對地下介電常數(shù)模型的有效重構。采用3層結構介電常數(shù)模型仿真模擬了數(shù)據(jù)低頻信息缺失的情況，驗證了包絡-波形反演在增強數(shù)據(jù)所攜帶的大尺度構造信息方面的能力與優(yōu)勢，有效克服了常規(guī)全波形反演在缺失低頻信息時無法正常工作的缺陷。此外，在所構建的包絡-波形反演框架下，又有針對性地加入了頻率多尺度策略，在重建大尺度宏觀背景構造的基礎上，也有利于實現(xiàn)對小尺度結構的精細刻畫。本文還進行了抗噪性測試并驗證了包絡-波形反演在噪聲條件下仍具有較強的反演能力。

2 包絡-波形反演

2.1 探地雷達信號包絡提取

因此，探地雷達信號的包絡也可以通過求取解析信號的振幅而獲得，即

圖1(a)和圖1(b)分別給出了主頻為80 MHz的源子波及其包絡的波形與頻譜數(shù)據(jù)；接下來，利用高通濾波器剔除了低于40 MHz的全部以及40～80 MHz的部分低頻成分，缺失低頻的子波及其包絡的波形與頻譜數(shù)據(jù)分別如圖1(c)和圖1(d)所示。從圖1(d)不難看出：盡管源子波缺失了低頻成分，但是其包絡數(shù)據(jù)仍包含豐富的低頻信息。

圖1 子波及其包絡的波形與頻譜信息

2.2 反演框架

常規(guī)包絡反演目標函數(shù)的定義方式借鑒了全波形反演定義目標函數(shù)的基本思想，即采用觀測數(shù)據(jù)包絡與模擬數(shù)據(jù)包絡的殘差的2范數(shù)作為目標函數(shù)[23]

2.3 頻率多尺度策略

在包絡-波形反演框架下，又進一步引入頻率多尺度策略，該策略的基本思想在于：借助濾波手段將觀測數(shù)據(jù)和子波分解為所需的頻率成分，在反演的初始階段先利用低頻成分構造大尺度背景速度構造，并將其作為后續(xù)反演的初始模型；接下來，再利用高頻成分精細刻畫地下介質中的小尺度弱擾動。這種多尺度策略從 “空間尺度”來細化地下介質模型，可以有效確保反演過程的正確收斂并提高反演結果的精度。

濾波是實現(xiàn)頻率多尺度的關鍵，這里借鑒Boonyasiriwat等人[26]的濾波思想，即通過維納濾波器將子波和觀測數(shù)據(jù)分解到不同的頻帶范圍。在包絡-波形反演框架下，借助子波包絡設計的維納濾波器公式為

其中，ω表示角頻率，fWi表示維納濾波器，Eor表示原始子波去低頻后取包絡，這里假定預先已知，Eta為目標頻率子波的包絡，ε表示控制數(shù)值溢出的穩(wěn)定因子。

2.4 包絡-波形反演實施步驟

本文所采用的包絡-波形反演可概括為一種分段式反演策略，即先通過包絡反演來重建介質的大尺度背景信息，并將其作為后續(xù)反演的初始模型，再利用常規(guī)全波形反演構造小尺度精細結構。值得注意的是：在低頻缺失的情況下，由于常規(guī)全波形反演無法構建高頻子波到低頻子波的維納濾波器，因此無法采用頻率多尺度策略。低頻成分缺失情況下包絡-波形反演的具體實施步驟可概括如下：

階段1 包絡反演：步驟1 輸入：初始模型、觀測記錄、源子波；步驟2 利用高通濾波器，構造缺失低頻的觀測記錄；

步驟3 選取頻率多尺度策略所需反演頻率，設定迭代次數(shù)；

步驟4 計算模擬記錄；

步驟5 利用式(5)計算缺失低頻觀測記錄的包絡和模擬記錄的包絡；

步驟6 確定當前迭代頻率，根據(jù)式(9)構造維納濾波器，作用于觀測記錄和模擬記錄包絡上，獲得指定頻率的觀測記錄和模擬記錄的包絡；

步驟7 根據(jù)梯度式(8)，計算整個模型的梯度；步驟8 利用不精確線搜索方法計算更新步長；步驟9 更新模型，迭代次數(shù)+1，并重復步驟4—步驟9，直到完成設定的全部迭代運算；

步驟10 輸出反演結果，并將其作為第2階段常規(guī)全波形反演的初始模型。

階段2 常規(guī)全波形反演：

步驟1 輸入：包絡反演的結果作為初始模型、觀測記錄、源子波；

步驟2 利用高通濾波器，構造缺失低頻的觀測記錄和源子波；

3 數(shù)值測試

3.1 完整頻率信息情況測試

為了確保反演過程的正確收斂及反演結果的置信度，本文采用了頻率多尺度策略，即按照2次函數(shù)的變化趨勢，從15～80 MHz中選擇了10個反演頻率，分別為15, 19, 24, 30, 37, 44, 52, 61, 70,80 MHz。包絡-波形反演的第1階段包絡反演和第2階段常規(guī)全波形反演，分別利用上述的10個頻率按由低到高的次序依次進行10次迭代，共計執(zhí)行200次迭代運算，包絡-波形反演結果如圖3(a)所示。為了便于進一步對比不同方法的反演效果，這里給出了僅采用常規(guī)全波形反演的結果(圖3(b))，反演過程中，上述10個反演頻率各進行20次迭代，以確保與包絡-波形反演的迭代運算次數(shù)一致。

圖3 3層結構介質模型的反演結果

兩種反演策略都能夠清晰地反映模型中異常的整體分布特征。其中，包絡-波形反演主要利用低頻成分重構大尺度背景信息，因此對于具有起伏底界面的中間層的刻畫要明顯優(yōu)于常規(guī)全波形反演。圖4給出了分別沿圖2(a)所示的真實模型和圖3所示的反演結果中不同位置(x=1.5 m, x=3.0 m, x=6.0 m)提取的參數(shù)縱向分布曲線；可以看出：包絡-波形反演結果中間層左下方介質的分布要更接近真實模型。盡管如此，由于常規(guī)全波形反演方法本身包含了全部波形信息，因此其在小尺度結構的刻畫精度和分辨率則更高。這里還分別計算了包絡-波形反演結果和常規(guī)全波形反演結果與真實模型的相關系數(shù)，分別為0.9868和0.9828，兩者與真實模型均極為相似，但包絡-波形反演重構能力略高。

圖2 3層結構介質模型

圖4 反演結果抽道對比

3.2 缺失低頻成分情況測試

在全波形反演中，低頻數(shù)據(jù)對保證反演的穩(wěn)定性至關重要。研究表明：采用低頻波場進行波形匹配時，發(fā)生“周波跳躍”的臨界值較大，因而低頻波場更有利于緩解反演的跳周問題。然而，在數(shù)據(jù)實際采集過程中，往往難以獲得數(shù)據(jù)有效的低頻分量，因此，研究低頻信息缺失情況下的反演方法則顯得十分重要。本節(jié)仍采用3層結構介質模型來測試低頻信息缺失下包絡-波形方法的反演能力。

首先，采用主頻為80 MHz的子波對如圖2(a)所示的介電常數(shù)真實模型進行正演計算，獲得探地雷達合成數(shù)據(jù)。接下來，利用高通濾波器濾除數(shù)據(jù)的低頻成分從而對低頻缺失情況進行仿真模擬。圖5給出了高通濾波前后第1個源對應的第1道雷達記錄及其相應的頻譜。從中可以看出原始雷達合成數(shù)據(jù)包含了完整的頻率信息，頻帶范圍為0～300 MHz(圖5(b))；而高通濾波器有效剔除了低于40 MHz的全部以及40～80 MHz的部分低頻成分(圖5(d))；此外，通過對比濾波前后的單道雷達記錄(圖5(a)和圖5(c))不難發(fā)現(xiàn)，缺失了低頻成分的雷達記錄的波形在振幅上也有所變化。

圖5 第1個源對應的第1道探地雷達合成記錄及其頻譜

基于2.4節(jié)給出的低頻成分缺失情況下包絡-波形反演的具體實施步驟，仍然采用3.1節(jié)所選取的10個反演頻率：在包絡反演階段，10個反演頻率先各進行10次迭代；由于低頻缺失條件下常規(guī)全波形反演無法應用頻率多尺度策略，因而常規(guī)全波形反演階段僅利用缺失低頻的源子波(原始主頻為80 MHz)進行100次迭代運算。由此，包絡-波形反演共計進行200次迭代更新，反演結果如圖6(a)所示。這里仍然給出了僅采用常規(guī)全波形反演的結果(圖6(b))，利用缺失低頻的源子波(原始主頻為80 MHz)進行了200次迭代運算。

圖7給出了分別沿圖2(a)所示的真實模型和圖6所示的反演結果中不同位置(x=1.5 m, x=3.0 m,x=6.0 m)提取的參數(shù)縱向分布曲線。通過兩種反演結果對比(圖6)以及抽道對比結果(圖7)不難發(fā)現(xiàn)：在低頻信息缺失的情況下，常規(guī)全波形反演恢復的介電常數(shù)模型圖像模糊，幾乎無法識別關于地下異常和地層分布的有效反射信號(圖6(b))；此外，與真實模型相比，重建結果的背景值被錯誤估計，推測反演結果可能陷入了局部極小值，由此可見，缺失低頻波場對全波形反演的重建質量會造成一定的影響。相比之下，包絡-波形反演很好地恢復了模型從淺到深的整體特征，對模型的還原度較高；反演圖像中，具有起伏底界面的層狀介質清晰可見，并且可以準確辨別模型深部嵌入的兩個夾雜體的位置信息，與3.1節(jié)具有完整頻帶信息的反演結果(圖3(a))相似度極高，這足以說明：包絡-波形反演可以在低頻成分缺失且常規(guī)全波形反演失效時發(fā)揮作用。此外，缺失低頻情況下，兩種反演結果與真實模型相關系數(shù)分別為0.9747和0.9351，相較于常規(guī)全波形反演，包絡-波形反演有效降低了反演誤差。

圖6 低頻缺失情況下模型的反演結果

圖7 反演結果抽道對比

3.3 抗噪性測試

在實際探地雷達勘探中，采集的雷達記錄中往往包含嚴重的干擾噪聲，這對后續(xù)的反演和成像都將造成極大的影響。因此，本節(jié)將對基于探地雷達數(shù)據(jù)的包絡-波形反演方法進行抗噪性測試。真實模型仍然采用如圖2(a)所示的3層結構介質模型。發(fā)射和接收天線的布設方式、反演參數(shù)、反演策略的選取均與3.1節(jié)無噪測試相同。在合成的觀測記錄中加入了經過均值濾波處理的白噪聲，這樣做是為了降低白噪聲的頻率以更好地匹配原始觀測記錄的主頻；噪聲總能量約為原始合成記錄總能量的0.5%。圖8(a)和圖8(b)分別給出了包絡-波形反演以及常規(guī)全波形反演的重建結果。兩種反演方法都能準確地恢復模型的宏觀構造信息，對具有起伏底界面的層狀介質的層位刻畫都比較準確，但是常規(guī)全波形反演圖像中，模型第2層左下角(4～6 m)介質的均勻性欠佳；圖9給出了分別沿圖2(a)所示的真實模型和圖8所示的反演結果中不同位置(x=1.5 m, x=3.0 m,x=6.0 m)提取的參數(shù)縱向分布曲線；可以看出：模型深度范圍在6 m以下的部分包絡-波形反演(藍色實線)要略優(yōu)于常規(guī)全波形反演(綠色實線)。通過分別計算兩種反演結果與真實模型的相似度也可有力證明上述結論(相關系數(shù)分別為0.9808和0.9784)。

圖8 添加噪聲后模型的反演結果

圖9 反演結果抽道對比

4 結束語

本文依據(jù)電磁波與地震波在動力學和運動學特征上的高度相似性，將地震勘探中用于解決低頻缺失情況的包絡反演方法遷延至探地雷達研究領域，構建了適用于地面多偏移距雷達采集模式的包絡-波形反演框架。經過研究分析，本文得到如下幾方面的結論：

(1)低頻數(shù)據(jù)的缺乏會影響大尺度擾動介質的恢復，從而導致常規(guī)全波形反演對初始模型具有高度依賴性。雷達記錄的包絡波動和衰減攜帶了有效的低頻信息，可用于估計大尺度(長波長)的背景結構。包絡-波形反演既保證了地下大尺度強擾動構造的有效重建，又兼顧了全波形反演可以精細刻畫小尺度弱擾動目標體的優(yōu)勢，因此在近地表電磁勘探領域具有很高的應用前景。

(2)本文通過對3層結構介電常數(shù)模型進行模擬研究發(fā)現(xiàn)：當雷達記錄包含了完整的頻帶信息時，常規(guī)全波形反演方法在細節(jié)信息的刻畫方面要優(yōu)于包絡-波形反演；但在低頻信息缺失的情況下，常規(guī)全波形反演則無法提供對地下物性參數(shù)的準確估計。相反，包絡-波形反演在低頻成分缺失時仍能給出穩(wěn)定且質量較高的反演結果。

(3)在包絡-波形反演框架下，加入了頻率多尺度的策略，實現(xiàn)了組間頻率由低到高串行反演，該策略可以有效壓制反演假象，提高反演精度。然而，低頻缺失的情況下，常規(guī)全波形反演方法無法構建高頻子波到低頻子波的維納濾波器，因而無法開展頻率多尺度策略。通過抗噪性測試，進一步驗證了包絡-波形反演方法對白噪聲具有不錯的抵抗性。