隧道施工不良地質探地雷達超前探測全波形反演研究

張崇民， 張鳳凱， 李 堯

(山東大學巖土與結構工程研究中心， 山東 濟南 250061)

0 引言

隨著我國經濟建設的發展，大量隧道工程破土動工，隧道工程施工中需要面臨的地質狀況越來越復雜。溶洞、巖體破碎帶、巖體裂隙等不良地質極易引發突水突泥、坍方等重大災害，并造成嚴重的經濟損失和人員傷亡[1]。因此，提前了解隧道掌子面前方賦存的不良地質體狀況對于施工安全具有極為重要的意義。

探地雷達(ground penetrating radar，GPR)基于向地下地質體發射高頻電磁信號并接收異常地質體反射信號來進行勘查，具有無傷、高效、高分辨率和對水體敏感的優點[2]，在隧道短距離不良地質超前探測中具有較好的應用效果。近年來，許多學者針對探地雷達隧道探測開展了廣泛的研究，并通過實際工程資料分析[3]、正演模擬[4]等總結了典型隧道不良地質體的響應特征和辨識方法，對于提高隧道不良地質超前探測的準確性起到了關鍵性作用。隨著探地雷達數據處理方法的發展，對于探測資料的解釋也由傳統的雷達剖面圖分析向偏移[5-6]、反演[7]、復信號分析[8-9]等多種方法發展，進一步提高了探地雷達數據分析的精度，為隧道探地雷達資料處理提供了可靠的思路?？偟膩碚f，探地雷達在隧道不良地質超前探測中獲得了廣泛應用，基于實際工程和正演模擬總結得到的不良地質體的響應特征，實現了對不良地質體的簡單辨識。由于電磁波傳播過程中存在多次反射波和繞射波，導致探地雷達數據中的同相軸走向與實際不良地質體的形態差異較大[4,8,10]，特別是地質狀況復雜區域的雷達探測結果中同相軸交錯疊加，僅憑經驗對雷達剖面圖分析無法獲得異常體的真實形態信息。

基于上述問題，本文提出將全波形反演運用到探地雷達數據處理中進行不良地質體判斷的方法。全波形反演(full waveform inversion,FWI)是一種基于全波場信息實現地下介質參數反演的高分辨率成像方法，最早由A. Tarantola[11]提出并應用于地震勘探領域，后來逐漸在雷達信號處理領域中應用發展。S. Busch等[12]實現了地面雷達數據的相對介電常數和電導率的定量估計； 吳俊軍[13]、孟旭[14]對時間域全波形反演方法做了系統推導及擴展； 李堯[15]將跨孔雷達運用到隧道施工不良地質超前探測中，并將全波形反演和逆時偏移成像方法相結合起來?？傮w來說，全波形反演方法能夠較為準確地反演出地下介質的介電常數信息，在隧道掌子面前方不良地質體定位和識別方面顯示出了良好的應用前景。

本文將全波形反演方法引入探地雷達探測數據處理流程，針對典型不良地質體進行探地雷達探測正演模擬和全波形反演研究，獲得探地雷達探測典型不良地質體的響應特征和反演成像規律，以期為實際工程應用提供理論基礎。

1 全波形反演基本原理及數據處理流程

全波形反演是一種基于全波場信息反演地下介質的反演方法，反演結果可以指示目標探測區域的較準確相對介電常數分布，成像結果更為直觀。全波形反演的思想是尋找使目標函數最小時的地下介質物性參數分布(本文搜索的是相對介電常數和電導率)。筆者采用的目標函數為正演模擬數據與實際觀測數據之間的二范數[13-16]。

(1)

式中：E(ε,σ),Eobs分別為正演模擬數據和實際觀測數據；s為探地雷達探測位置；t為探測時間。

當目標函數的值最小時，正演模擬數據和實際觀測數據最接近，此時得到的相對介電常數分布也就越接近真實模型。本文使用共軛梯度法進行目標函數尋優，目標函數梯度可由式(2)和式(3)計算[13-16]。

(2)

(3)

式(2)和式(3)中：I=GT·R；R=E-Eobs(R為殘差，G為格林算子，I為逆時傳播波場數據，均由正演模擬獲得)。

介電常數和電導率更新的迭代公式如式(4)和式(5)[13-16]所示。

εk+1=εk-ζε,k·gε。

(4)

σk+1=σk-ζσ,k·gσ。

(5)

式(4)和式(5)中：k為迭代次數；ζε,k和ζσ,k為第k次迭代中相對介電常數和電導率的迭代步長。

ζε,k和ζσ,k可分別表示為[13-16]：

(6)

(7)

式(6)和式(7)中：κε和κσ為小穩定因子，反演過程中必須選擇適當的值，并且隨著迭代更新。

根據上述梯度和迭代步長不斷更新，直至達到終止條件，完成全波形反演。

探地雷達探測數據全波形反演處理過程，具體分為以下8個步驟[11,13]。

1) 建立初始模型；

2) 將源子波加載到發射天線，進行正演模擬，并記錄全時刻正傳波場和接收天線處的正演波形；

3) 將實際接收數據與正演數據相減得到殘差E-Eobs和目標函數S；

4) 將殘差E-Eobs在時間上倒置，并以此為源加載到接收天線上，計算逆時殘場；

5) 通過正傳波場和逆時殘場計算相對介電常數迭代梯度gε和電導率迭代梯度gσ；

6) 計算對應于梯度方向的相對介電常數迭代步長ζε和電導率迭代步長ζσ；

7) 更新模型參數；

8) 計算目標函數并判斷是否達到終止條件，如未達到終止條件則跳轉至步驟2)，如已達到終止條件則完成全波形反演過程。

本文在基于全波形反演方法處理探地雷達正演模擬數據時，首先假設初始模型為均勻模型，使用全波形反演方法得到不良地質區域的相對介電常數。全波形反演流程如圖1所示。

圖1 全波形反演流程

2 隧道施工典型不良地質體探地雷達探測數值算例

隧道圍巖地質情況復雜多變，各種不良地質造成的安全事故也屢見不鮮，隧道施工不良地質常見于溶洞、巖體裂隙和巖體破碎帶等[17]。基于上述不良地質類型，本文以空氣、灰巖、潮濕黏土、水4種介質設計了4個數值算例,進行正演模擬和全波形反演結果分析。介質電性參數設置如表1所示。

表1 介質電性參數設置[18]

本文探地雷達采用單發射單接收型平板天線，發射天線和接收天線收發距固定。發射電磁波中心頻率為100 MHz，采樣時間為365 ns。圖2為探地雷達隧道檢測示意圖，左側為探地雷達測線，平板天線在掌子面上橫向移動，間隔0.1 m進行1次探測，并得到1道數據。當電磁波在掌子面前方巖體內傳播至不良地質體界面時發生發射，反射信號被接收天線接收。本文數值算例中采用的源子波類型為雷克子波，運用Matlab編制的時域有限差分程序對寬10 m、長20 m的模型進行模擬，正演網格空間步長0.05 m，整個程序運算均在1臺CPU配置為Intel Xeon E5-2603的小型服務器上進行。通過Matlab內置的parfor函數，在1次正演過程中同時進行12線程的運算，耗時約為81 s。單次迭代過程共包括正傳波場計算、逆時殘場計算和2次迭代步長計算等4次正演過程。全波形反演過程中對于內存需求主要集中在幾千個時間步長的正傳波場場量存儲上。在程序計算中，模型網格數為200×400，單個網格存儲類型為single單精度類型，大小為4字節，并且需要同時存儲12個線程的場量，這就需要產生大約10 GB的內存。

圖2 探地雷達隧道檢測示意圖

2.1 溶洞模型數值算例

溶洞成因于巖石的巖溶作用，其物性參數與周圍巖體存在明顯差異。溶洞內部常常會充填碎黏土、水、空氣等物質，使溶洞與周圍巖體之間的物性差異更加明顯[19]。溶洞模型如圖3所示。設置溶洞模型Ⅰ、模型Ⅱ、模型Ⅲ如圖3(a)、3(d)、3(g)所示，背景介質為灰巖，在距離掌子面5，7，13 m處分別設置3個大小約為1 m×1 m，2 m×2 m，4 m×4 m的不規則溶洞，溶洞內部分別充填空氣、潮濕黏土和水。在圖3(b)、3(e)、3(h)所示的正演結果中，黑色虛線為溶洞位置。溶洞響應均呈現出連續的雙曲線同相軸，弧頂處為溶洞的近掌子面邊界，遠端界面及溶洞形狀均無法獲得。模型Ⅱ和模型Ⅲ正演結果中的反射信號出現了明顯的同相軸反轉。由于電磁波信號在巖層中衰減，溶洞③響應并不強烈，反射信號勉強能夠辨識出其近掌子面位置。在如圖3(c)所示的模型Ⅰ全波形反演結果中，溶洞①和溶洞②的反演結果較為理想，溶洞形狀基本被反演出來，結果中指示的溶洞位置分別位于掌子面前方5 m和7 m處，與實際位置一致，相對介電常數值為0～2，也與預設模型相符合。模型Ⅰ正演結果中雙曲線弧面較長，僅能判斷近掌子面界面，而無法分辨其橫向寬度。溶洞③的反演結果并不理想，但是仍能辨識出溶洞③的前界面在掌子面前方13 m處及大致橫向寬度為4 m。在如圖3(f)和圖3(i)所示的模型Ⅱ和模型Ⅲ的全波形反演結果中，溶洞①和溶洞②的形狀與模型吻合較好，模型Ⅱ反演的溶洞①和溶洞②的介電常數值與預設模型接近，但是模型Ⅲ中溶洞①和溶洞②的反演結果在12～15，與預設模型中水的相對介電常數81差距較大。模型Ⅱ和模型Ⅲ的全波形反演結果中，均沒有反映出溶洞③的介電常數范圍,只能對其近掌子面位置進行判斷,考慮原因為接收到的溶洞③反射信號微弱，同時介質電性參數差異過大,影響反演結果。

(a) 溶洞充填空氣模型(模型Ⅰ)

(b) 模型Ⅰ正演結果

(c) 模型Ⅰ全波形反演所得相對介電常數分布

(d) 溶洞充填潮濕黏土模型(模型Ⅱ)

(e) 模型Ⅱ正演結果

(g) 溶洞充填水模型(模型Ⅲ)

(i) 模型Ⅲ全波形反演所得相對介電常數分布

2.2 巖體裂隙模型數值算例

巖體裂隙主要成因于巖體內部應力作用或地殼運動，其往往會造成巖體應力失穩[20]。巖體裂隙模型如圖4所示。設置巖體裂隙模型Ⅳ、模型Ⅴ、模型Ⅵ如圖4(a)、4(d)、4(g)所示，背景介質為灰巖，在距離掌子面7～12 m處設置2個傾角(以掌子面為基準)分別為30°、40°的裂隙，其長度分別為6 m和5 m，內部分別填充空氣、潮濕黏土、水。在圖4(b)、4(e)、4(h)所示的模型正演模擬結果中，黑色虛線指示裂隙長度走向。從電磁波反射信號中可以清楚地觀察到2組反射非常強烈的同相軸，方向與設置模型中的裂隙走向呈交叉狀，裂隙尖端呈現拋物線狀繞射信號，圖4(e)和圖4(h)所示模型Ⅴ和模型Ⅵ的正演結果中反射信號同相軸反轉。由于同相軸走向與真實裂隙形態差異較大，單純基于正演剖面圖分析無法獲知前方巖體裂隙形態。圖4(c)和圖4(f)分別為模型Ⅳ和模型Ⅴ的全波形反演結果。預制模型中灰巖的相對介電常數為7，空氣的相對介電常數為1，潮濕黏土的相對介電常數為15，反演結果中指示的背景圍巖相對介電常數值為6～8，圖4(c)和圖4(f)中異常體相對介電常數分別為0～2和14～15，均與模型較為吻合，反演得到的裂隙形態走向、長度與模型相比也較為相符，基本實現了對模型的復現，能夠通過全波形反演結果簡單地辨識出掌子面前方巖體裂隙的實際形態。圖4(i)所示模型Ⅵ的反演結果中，裂隙形態反映較好，介電常數值在18～20，與設置模型參數相差較大，但是與圖3(i)所示的模型Ⅲ反演結果相比，模型Ⅵ反演得到的高介電常數分布更加理想，其可能原因為裂隙狹長形態使得在探測中能夠反射更多的電磁波信息。

2.3 巖體破碎帶模型數值算例

巖體破碎帶是完整巖體在構造運動作用下產生的破碎區域，主要表現為區域巖體分塊破碎，塊體空隙內部可能充填空氣或水等[21]。巖體破碎帶模型如圖5所示。設置巖體破碎帶模型Ⅶ、模型Ⅷ、模型Ⅸ如圖5(a)、5(d)、5(g)所示，背景介質為灰巖，在距離掌子面8～12 m處為破碎區域，破碎巖體空隙分別填充空氣、潮濕黏土和水。圖5(b)為模型Ⅶ正演模擬得到的雷達剖面圖，可以看出異常體區域的電磁波反射信號同相軸交叉錯亂，因此考慮異常體內部介質性質變化較大； 同時異常信號反射區域的前后邊界面較為清晰，指示的異常區域也與設置模型基本吻合。然而，僅依靠正演結果分析無法對異常信號區域介質性質進行判斷。圖5(c)為模型Ⅶ的全波形反演結果，相較于正演結果，破碎帶前界面走向在圖5(c)得到了良好的重現； 異常體區域的相對介電常數主要分布在1～2和6～8，考慮到空氣和灰巖的相對介電常數，可以推斷不良地質區域存在大量破碎，破碎區域內部充填空氣。反演得到的低介電常數位置，基本實現了對巖體內部小空隙的歸位，基于全波形反演結果，可以更加直觀地了解破碎區域內部的巖體碎塊和空隙分布狀態。圖5(e)和圖5(h)分別為模型Ⅷ和模型Ⅸ正演結果，可以看出電磁波反射信號主要集中在破碎帶淺層區域，同相軸呈現明顯反轉。與圖5(b)所示模型Ⅶ正演結果相比，電磁波在含水地層中衰減比在空氣中更加嚴重，破碎帶深部區域反射信號微弱，特別是模型Ⅷ設置的潮濕黏土中，深層電磁波信號幾乎完全衰減。圖5(f)和圖5(i)分別為模型Ⅷ和模型Ⅸ的全波形反演結果，可以看出掌子面前方8 m處出現有明顯的相對介電常數高值區域，其范圍分別為12～13和9～10，且高值區域走向和預設模型吻合。模型Ⅸ的反演結果較差，反演結果與實際模型參數差值在70左右，但是考慮其較高的介電常數，仍可判斷異常區域含水。由于電磁波信號在含水地層衰減，模型Ⅷ和模型Ⅸ中破碎帶深層區域反演結果均不理想。

(a) 裂隙充填空氣模型(模型Ⅳ)

(b) 模型Ⅳ正演結果

(c) 模型Ⅳ全波形反演所得相對介電常數分布

(d) 裂隙充填潮濕黏土模型(模型Ⅴ)

(e) 模型Ⅴ正演結果

(f) 模型Ⅴ全波形反演所得相對介電常數分布

(g) 裂隙充填水模型(模型Ⅵ)

(h) 模型Ⅵ正演結果

(i) 模型Ⅵ全波形反演所得相對介電常數分布

(a) 巖體破碎帶充填空氣模型(模型Ⅶ)

(b) 模型Ⅶ正演結果

(c) 模型Ⅶ全波形反演所得相對介電常數分布

(d) 巖體破碎帶充填潮濕黏土模型(模型Ⅷ)

(e) 模型Ⅷ正演結果

(f) 模型Ⅷ全波形反演所得相對介電常數分布

(g) 巖體破碎帶充填水模型(模型Ⅸ)

(h) 模型Ⅸ正演結果

(i) 模型Ⅸ全波形反演所得相對介電常數分布

2.4 不良地質組合模型數值算例

在實際工程中，地質狀況更加復雜，多類不良地質體常常相伴出現，因此，本文在上述不良地質類型的基礎上設計了組合地質模型。不良地質組合模型如圖6所示。設置模型Ⅹ如圖6(a)所示，背景介質為灰巖。在掌子面前方9 m處設置1個大小為2 m×4 m的橢圓形溶洞，溶洞兩端連接裂隙，內部充填水。在模型Ⅹ的正演結果(圖6(b))中，可以看到3組不同走向的同相軸，中間1組雙曲線同相軸為溶洞響應，弧頂處為溶洞的近掌子面位置； 兩側同相軸為裂隙響應，其位置與裂隙的實際位置存在偏差。在模型Ⅹ的全波形反演結果中，裂隙起于掌子面前方8 m處，終止于掌子面前方12 m處，呈120°夾角，與實際模型對應良好； 反演結果中指示的內部介質相對介電常數值在12～15，與設置模型相差較大，但仍可依據反演結果中裂隙區域較高介電常數值推測地層含水； 溶洞前界面的反演結果較好，其位置位于掌子面前方9.5 m處，未能反演出溶洞的后界面。設置模型Ⅺ如圖6(d)所示，背景介質為灰巖。在掌子面前方8～9 m處設置1條破碎帶，破碎帶后方距掌子面13 m處有1個大小約為3 m×3 m的圓形溶洞，內部充填空氣。在模型Ⅺ的正演結果(圖6(e))中，可以看到明顯的雜亂反射區域，反射信號同相軸交錯疊加； 在模型中后部能觀察到一近似平行于掌子面的同相軸，但是反射信號響應十分微弱。在模型Ⅺ的全波形反演結果中，破碎位置與設置模型吻合，反演得到的介電常數值在2～3，與設置模型相近。溶洞反演結果較差，在圖6(f)中能辨識到其近掌子面界面，且介電常數較低，在3～5。設置模型Ⅻ如圖6(g)所示，背景介質為灰巖。在掌子面右側前方8 m處存在4 m×4 m的巖體破碎區域，破碎區域一端連接裂隙，內部充填空氣。模型Ⅻ的正演結果如圖6(h)所示，能明顯辨識出一傾斜同相軸，為裂隙位置。破碎區域電磁波反射信號較為雜亂，且多條同相軸與裂隙反射信號交錯。通過正演結果僅能識別裂隙位置，無法分辨破碎區域。在模型Ⅻ的全波形反演結果中將裂隙的位置很好地呈現出來，指示的內部介質相對介電常數值在1～2，與設置模型一致。反演結果中巖體破碎位置和范圍與設置模型吻合較好，對于破碎區域的相對介電常數反演結果，近掌子面結果與模型較為一致，深部區域的介電常數反演較差。

(a)溶洞-裂隙充填水模型(模型Ⅹ)

(b) 模型Ⅹ正演結果

(c) 模型Ⅹ全波形反演所得相對介電常數分布

(d) 巖體破碎-溶洞充填空氣模型(模型Ⅺ)

(e) 模型Ⅺ正演結果

(f) 模型Ⅺ全波形反演所得相對介電常數分布

(g) 巖體破碎-裂隙充填空氣模型(模型Ⅻ)

(h) 模型Ⅻ正演結果

(i) 模型Ⅻ全波形反演所得相對介電常數分布

2.5 全波形反演適應性分析

影響全波形反演的主要因素有初始模型、源子波形態、反演目標體形態等。本文所采用的初始模型為均勻模型，得到了良好的全波形反演結果。若基于先驗地質信息建立初始模型，可以進一步對反演結果進行優化。在數值模擬中，源子波形態是已知的，但是在處理實際數據時，并不知道實際的發射波形，這就需要對源子波進行估計，并在每次迭代過程中，更新源子波信號。全波形反演對于狹長的、平行于掌子面的不良地質體有更好的反演效果，其原因為狹長的界面能夠反射更多的電磁波信號?；诒疚恼故镜牧严赌Ｐ头囱萁Y果，均清晰地顯示出了多條裂隙的走向和形態。在溶洞模型中，目標體相對于雷達測線的反射面較小，能夠獲得的反射信息較少。對于靠近掌子面的小型溶洞，反演結果較好，位置及相對介電常數均較為清晰，而深層溶洞僅能反演得到前界面位置； 特別是在設置的組合模型Ⅺ中，破碎帶后方的溶洞響應極其微弱，反演結果得到的前界面也非常微弱。在破碎帶模型的反演結果中，破碎區域能夠反射更多的電磁波信號，淺層破碎區域走向較為清晰。在實際工程中，由于并不知道不良地質體相對于掌子面的走向和形態，因此可以采用2條測線(1條水平、1條垂直)，對多個方向的地質雷達數據進行反演分析。當反演區域內部的相對介電常數變化過大時，全波形反演過程中對于區域介電常數的迭代更新往往不足，反演結果與實際模型存在較大差值，像充填水模型反演得到的介電常數值與實際差值在50以上，但是能觀察到介電常數高值區域。此時，應結合實際地質資料進行綜合評判。

3 結論與討論

針對傳統隧道探地雷達數據處理中僅能依靠經驗對不良地質體進行判斷的局限性，將全波形反演方法應用到探地雷達數據處理中，針對典型不良地質體建立數值模型，通過系統的正演模擬和全波形反演研究，取得了較好的效果。

1)電磁波在圍巖中的傳播特征較為復雜，雷達剖面圖中的同相軸與實際異常體形態往往存在較大差異。在巖體裂隙數值算例中，電磁波反射信號同相軸與實際裂隙走向完全不同，直接依據雷達剖面圖無法獲得異常體的形態特征。

2)采用全波形反演方法可以反演出掌子面前方異常體的電性參數和位置形態，4個數值算例的全波形反演結果均有助于對預設模型的判斷。例如： 含水地層在反演結果中表現出較高的相對介電常數值，通常在15以上； 空洞在反演結果中表現出較低的相對介電常數，在0～2。

3)當不良地質體介電常數差距懸殊時，全波形反演結果不理想； 同時，電磁波在傳播過程中的衰減會影響全波形反演的結果。例如，溶洞算例中后界面反射信息較少、巖體破碎帶算例中深部破碎反射信號微弱，都影響了相關區域的全波形反演結果，但是綜合先驗地質資料分析，仍然可以對前方地質狀況進行大致判斷。

4)由于實際探測地質雷達數據信噪比較低，全波形反演在隧道超前探測中的應用效果并不理想。下一步將重點開展實際數據濾波降噪算法，以提高實際數據全波形反演效果。