砂巖體裂隙滲流區(qū)的探地雷達(dá)信號(hào)動(dòng)態(tài)響應(yīng)機(jī)制

關(guān)鍵詞：探地雷達(dá)；裂隙滲流；水文模擬；瞬時(shí)屬性

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

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

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

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

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

1.1 砂巖裂隙滲流模擬

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

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

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

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

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

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

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

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

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

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

2.1 裂隙水動(dòng)力學(xué)模擬

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

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

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

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

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

3 分析

3.1 瞬時(shí)振幅屬性提取

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

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

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

3.2 含水飽和度與瞬時(shí)振幅關(guān)系

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

4 結(jié)論

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