砂巖體裂隙滲流區(qū)的探地雷達信號動態(tài)響應機制

關鍵詞：探地雷達；裂隙滲流；水文模擬；瞬時屬性

石窟寺作為石質(zhì)文物的典型代表，具有極高的歷史、藝術和科學價值[1]。近年來，我國對石窟類文物保護的工作力度逐步加大。巖體裂隙滲水病害是我國石窟寺最普遍、危害最嚴重的病害類型之一，具有滲水通道微小，滲流結構復雜，滲流現(xiàn)象微弱等特點，對它的探測和治理是石窟寺保護工作的重點也是難點。

由于文物的不可再生性，文物裂隙的探測要盡可能地使用無損檢測的手段，探地雷達（GPR）是研究和保護具有較高文化和歷史價值古代建筑的理想無損探測方法[2]。裂隙通常是被嵌在均勻巖層中的一層[3]。探地雷達對裂隙的檢測來源于薄層反射理論，裂縫和圍巖之間存在電性差異，電磁波在裂縫壁內(nèi)產(chǎn)生多次反射，使得亞米級波長的電磁波對毫米級孔徑的裂縫產(chǎn)生響應（薄層響應）[4?5]。探地雷達方法已經(jīng)探討了對巖體內(nèi)裂縫的檢測[6?7]。Vickers等[8]將探地雷達技術使用在墨西哥州的科峽谷印第安遺跡勘探過程中，這是探地雷達在考古與文物保護領域中的首次應用，此后，探地雷達技術頻繁地使用在文物保護過程中[4，9?10]。裂縫控制著巖石中流體的流動[11]。Giertzuch[12]把時移差分探地雷達探測法應用于監(jiān)測生理鹽水示蹤劑在亞毫米裂隙中的流動。Dong等[13]利用一種新的時移全極化探地雷達和各向異性分析方法檢測毫米尺度裂縫。Eskandari[14]使用SVD算法分析探地雷達圖像以減少雜波，使得目標信號明顯增強，該技術能夠檢測到寬度大于1.3mm的裂縫。YulianitaDina[15]使用天線頻率為60MHz的探地雷達檢測湖體滲水，探測出湖體滲水的深度。Guo等[16]使用時移探地雷達數(shù)據(jù)繪制土壤中的地下優(yōu)先流動。研究表明，探地雷達多數(shù)應用于工程探測和水文地質(zhì)領域的裂隙探測，在石窟巖體微裂隙滲流探測方面的研究還應加強。

因此，有必要針對石窟寺裂隙的特點，進行流體滲流過程中電磁波場響應及其機理研究。文中首先構建了砂巖裂隙滲流水文模型，對裂隙水的滲流過程和不同時間點的滲流狀態(tài)進行了瞬態(tài)模擬，隨后進行雷達波場數(shù)值模擬，分析其波場特征，并利用希爾伯特變換對正演模擬的結果進行瞬時振幅的提取與分析，總結出不同時刻裂隙滲流模型探地雷達的波場與瞬時屬性特征。在實際砂巖裂隙滲流的探測中，利用實際的波場特征判斷滲流裂隙的位置，對比不同滲流裂隙的含水量大小，為探地雷達實現(xiàn)超精度探測砂巖體微小滲流裂隙提供重要的理論依據(jù)。

1 針對裂隙滲流的探地雷達數(shù)值模擬方法

將砂巖裂隙水動力學模擬的模型與探地雷達波場數(shù)值模擬相結合，分析了裂隙滲流和雷達波響應的機制關系，流程如圖1所示。首先，基于Richards方程[17]和VanGenuchten持水模型[18]建立砂巖裂隙滲流水文模型，再將水含量轉(zhuǎn)換為雷達波場模擬所需的介電常數(shù)和電導率值，通過水含量與瞬時振幅的關系得到雷達波地下傳播特征對砂巖體裂隙滲流的動態(tài)響應機制。

1.1 砂巖裂隙滲流模擬

通過求解滲流場滿足的偏微分方程式和邊界條件式實現(xiàn)滲流有限元的計算分析。Richards方程[17]是由達西定律[19]演變而來的，描述了可變飽和多孔介質(zhì)中的流動，闡述流體在填充某些孔隙并從其他孔隙排出時水力屬性的變化，表明非飽和水流也符合達西定律，即非飽和水流的滲流速度也與總土水勢能成正比，且與土體中孔隙通道的幾何性質(zhì)有關[19]。

通過求解Richards方程和VanGenuchten模型有效模擬了砂巖中裂隙水的滲流過程。在公式（1）中，設置液體密度ρ為1000kg/m3，流體動力黏度μ為10-3Pa·s，表示水的流體屬性。根據(jù)砂巖性質(zhì)，在VanGenuchten模型中，n的值被設置為2。最后，利用Comsol軟件對式（1）（2）求解，得到巖體裂隙模型在不同時刻的含水飽和度。

1.2 砂巖裂隙滲流電性模型建立

利用探地雷達測定土壤含水率的Topp公式[20]，將砂土的含水飽和度代入公式（3）轉(zhuǎn)化為相對介電常數(shù)：

式中：εr表示相對介電常數(shù)；θ表示含水飽和度。砂巖模型電導率可利用砂土含水量與電導率的經(jīng)驗公式進行計算[21]：

水文模擬結果可以進一步應用于探地雷達數(shù)值模擬。利用中公式（3）和（4）可以計算出相對介電常數(shù)εr和電導率σ，通過數(shù)值模擬得到雷達波場響應振幅。這一過程完成了從水文參數(shù)到電性參數(shù)，再到雷達波場響應振幅的轉(zhuǎn)換。

1.3 砂巖裂隙滲流探地雷達數(shù)值模擬

為了準確地模擬探地雷達探測砂巖體裂隙滲流模型時的電磁波在介質(zhì)中的傳播過程，使用時域有限差分數(shù)值模擬進行砂巖體裂隙滲流模型的電磁波場響應計算，完成從電性參數(shù)到探地雷達響應振幅的轉(zhuǎn)換。時域有限差分法是一種常用的直接求解偏微分方程邊值問題的數(shù)值方法[22]。Yee[22]提出了Yee網(wǎng)格空間離散方式，將目標空間進行網(wǎng)格離散，電場和磁場分量分布于網(wǎng)格之中，在交織的網(wǎng)格空間中電磁場進行交替迭代計算，當數(shù)值計算滿足穩(wěn)定性條件并且收斂時，完成計算過程。FDTD計算方法有很多優(yōu)點，例如，簡單直觀易于理解，計算量小，可行性高，被廣泛應用于工程領域。

時域有限差分法由麥克斯韋旋度方程的微分形式出發(fā)。Maxwell方程的旋度方程組為

2 砂巖體中陡傾角單裂隙滲流的探地雷達數(shù)值模擬

2.1 裂隙水動力學模擬

石窟寺及摩崖造像在我國石質(zhì)文物中處于主體地位，其內(nèi)容豐富、規(guī)模龐大、數(shù)量眾多、分布廣泛，具有地質(zhì)體和人工營造建筑及藝術品的雙重特征[23]。石窟寺文物病害定義的范疇，既要考慮石窟寺文物本體的結構及雕塑藝術品的特點，又要考慮依托自然山體地質(zhì)環(huán)境的特點，同時，還要結合我國石窟寺及石刻保護工作的特點和保護實踐的需求[24]。在石窟寺保護領域，石窟寺文物載體的概念是對石窟寺文物本體起支撐作用，并與文物本體結構安全密切相關的構造體及地質(zhì)環(huán)境，包括石窟寺洞窟構筑物與雕刻藝術品所依托的山體，以及相關聯(lián)的地質(zhì)體[24]。文中的研究對象位于重慶大足的北山石刻實驗場，與大足石刻摩崖造像相連，是文物本體相關聯(lián)的地質(zhì)體。據(jù)資料記載，北山造像巖體主要分上下2個巖性層，含多條泥質(zhì)弱層，上部棕褐色砂巖，對應崖面位置處于石窟造像以上，含3條明顯的弱層；下部位于造像砂巖以下，主要為一套紫灰、紫紅色厚層塊狀細粒長石石英砂巖和紫紅色泥巖、粉砂質(zhì)泥巖組成，兩部分整合接觸[25]。

根據(jù)大足北山石刻區(qū)域地質(zhì)條件，建立了砂巖質(zhì)裂隙滲流模型，模擬了砂巖中裂隙水的滲流過程，結合瞬態(tài)計算探究水在巖體內(nèi)部的滲流規(guī)律。建立大小為0.8m×0.8m的砂巖體裂隙滲流模型，如圖2所示。裂隙起始位置為（0.5m，0.8m），末端位置為（0.3m，0.3m），裂隙厚度df設置為0.05m，背景介質(zhì)為砂巖。裂隙滲流由裂隙頂部端點給定的壓力所產(chǎn)生，由公式（1）和（2）計算的裂隙滲流0、1、20、40、80s時刻的體積含水量模型分別如圖3所示。由圖可知，當裂隙中有水向外滲出時，形成了裂隙周邊的含水區(qū)域，以及外圍還未入滲的不含水區(qū)域。

2.2 裂隙滲流電性模型和探地雷達數(shù)值模擬

對上文6個時刻的滲流模型進行探地雷達數(shù)值模擬。正演模擬參數(shù)如表1所示。其中，激勵源采用雷克子波，天線中心頻率為1.6GHz。天線距離模型5cm，采用發(fā)射天線和接收天線同時移動的方式，每移動2.5cm采集一道數(shù)據(jù)，每個模型共采集了26道數(shù)據(jù)。由含水飽和度與模型介電常數(shù)的關系公式（3）計算得到數(shù)值模擬0、1、20、40、60、80s時刻的滲流介電模型如圖4所示。

圖5分別對應圖4滲流時刻的模型模擬結果，從雷達剖面可以看出，由于滲流點與周邊圍巖介電常數(shù)的差異，雷達的反射波在剖面中呈現(xiàn)雙曲線的形態(tài)。但隨著水不斷在介質(zhì)中擴散，含水的高介電常數(shù)介質(zhì)的面積不斷擴大，雙曲線特征逐漸消失。

3 分析

3.1 瞬時振幅屬性提取

使用時域有限差分法進行探地雷達數(shù)值模擬得到時間域內(nèi)的信號Ez（t），使用Hilbert變換實現(xiàn)瞬時信號的提取。Hilbert變換是復信號分析中一項重要的工具，Hilbert變換又稱為90。相移濾波，實質(zhì)是將信號的相位譜做90。相移，而保持振幅譜不變[26]。在時間域內(nèi)連續(xù)信號Ez（t）的Hilbert變換定義為

由式（12）可知Hilbert變換就是在時間域內(nèi)將一個原始信號與Hilbert變換因子h（t）=1πt進行卷積運算。瞬時屬性的概念來源于地震勘探。雖然，地震勘探與探地雷達所測參數(shù)和源不同，但有著相同的運動學特征，因此，地震勘探中的研究手段也可以用于探地雷達的采集和處理[27]。Young等[28]首次把地震屬性分析技術用于探地雷達數(shù)據(jù)分析；Senecha等[29]用地震處理和解釋軟件進行探地雷達屬性的提取和分析，探地雷達屬性分析技術取得了較大進展。瞬時振幅可以衡量探地雷達反射波的強度，反映能量因介質(zhì)介電常數(shù)差異、介質(zhì)吸收和傳播距離等因素變化而發(fā)生的衰減趨勢。從數(shù)學計算公式角度分析，瞬時振幅等于探地雷達反射信號的實部和虛部總能量的平方根，可有效地突出特殊巖層的變化情況[30]。探地雷達接收天線中接收的電磁波信號的三角函數(shù)表達式為式中，A（t）即為瞬時振幅。

探地雷達屬性分析能夠提取雷達反射波中記錄的信息。對預處理后的原始信號Ez（t）使用公式（12）進行Hilbert變換得到瞬時振幅A（t），并做出圖6不同時刻滲流模型的瞬時振幅剖面。瞬時振幅正比于該時刻雷達反射信號總能量的平方根，體現(xiàn)了因介質(zhì)介電常數(shù)差異、介質(zhì)吸收和傳播距離等變化而發(fā)生的能量衰減趨勢[31]。由瞬時振幅屬性剖面可以看出，隨著滲流時間的增加滲流面積不斷擴大，導致電磁波傳播穿過介電常數(shù)較高的區(qū)域時發(fā)生強反射，能量較強的亮帶因此不斷變寬。

3.2 含水飽和度與瞬時振幅關系

為了更清晰地分析含水飽和度與雷達響應特征之間的關系，選取了20、40、60、80s時刻部分單道信號的瞬時振幅最大值，提取其對應的含水飽和度，得到含水飽和度與瞬時振幅交會圖，如圖7所示。通過函數(shù)擬合，得到雷達反射波瞬時振幅隨飽和度變化的曲線。可以看出，瞬時振幅與含水飽和度成正比關系，表明隨著滲流時間的增加，含水飽和度不斷增大，雷達反射波的振幅也隨著增大。

4 結論

文中提出了一種基于水文模型和探地雷達數(shù)值模擬分析石窟巖體裂隙滲流區(qū)探地雷達瞬時屬性響應特征的方法。水文建模模擬了水在裂隙中的流動以及在圍巖中的擴散現(xiàn)象，基于結果建立了不同滲流時刻的介電常數(shù)模型，進行探地雷達正演模擬，分析了含水飽和度與雷達反射波瞬時振幅之間的關系。得出以下結論，在0s時刻，滲流處于初始值時刻，滲流點相當于類表面散射體，雷達響應呈雙曲線特征；在1s時刻，滲流區(qū)只有距裂隙周圍很近的地方，滲流體積很小，滲流區(qū)和含水裂隙與所使用電磁波的波長相比可近似為一個長度無限長寬度較小的偶極子散射體，雷達剖面呈現(xiàn)雙曲線特征[32]；在20、40、60、80s時刻，滲流區(qū)面積不斷擴大，滲流模型相當于平面散射體，雙曲線特征消失。在瞬時振幅剖面的結果中可以看出，水隨著時間在裂隙巖體中擴散，導致電磁波傳播穿過介電常數(shù)較高的區(qū)域時發(fā)生強反射而形成的不斷變寬的亮帶。通過函數(shù)擬合，得到當含水飽和度增大時，雷達波反射的瞬時振幅也增大的趨勢。結果表明，探地雷達具有探測滲水裂隙的能力，能夠?qū)柡投瓤焖僮兓臐B水裂隙區(qū)做出不同的響應，這在石質(zhì)文物的保護中發(fā)揮著重要的作用。研究結果為探地雷達實現(xiàn)超精度探測砂巖體微小滲流裂隙提供了重要的理論支撐。