摘要：GPR管道探測數值模擬主要採用基於Maxwell方程組的時域有限差分（FDTD）法與基於幾何光學原理的射線追蹤法，本文對兩種方法進行了闡述。並對地下管道模型GPR數值模擬與物理模擬特徵進行了探討。研究結果表明：管道周圍介質導電率與對雷達波的吸收衰減成正比，同相軸幅度隨著來自管道頂、底回波能量的減弱而降低。
關鍵詞： GPR管道探測模擬 射線追蹤法 探地雷達 時域有限差分模擬
PR（GPR）是利用天線向介質中發射高頻電磁脈衝波來探測介質中目標體的一種物理探測技術；GPR探測的物理基礎是被探測目標與周圍介質之間存在電磁性（導電率、介電係數及磁導率）差異[1-2]。GPR通常採用反射剖面法進行探測，GPR剖面記錄來自地下介質中的電磁反（散）射回波的走時及幅度資訊，根據GPR剖面上來自目標體回波的運動學與動力學特徵來推演地下目標的位置、埋深及屬性等相關資訊。
通常條件下，由各種不同材質構成的地下管道與其周圍土壤介質之間的導電率與介電係數的差異幾乎是普遍存在的；另外，各種地下管道一般都具有較平整、光滑的外表面，因而對雷達波而言是良好的散射截面。因此，GPR探測地下管道具有優良的地球物理前提條件。與傳統的基於電磁感應原理的管線探測技術相比，由於GPR具有解析度高、影象直觀，而且既能探測金屬管道也能探測非金屬管道等顯著優點[3-4]，因此，GPR目前已成為探測地下管道的主要方法之一。
1、GPR管道探測數值模擬方法介紹
目前，GPR管道探測數值模擬主要採用基於Maxwell方程組的時域有限差分（FDTD）法與基於幾何光學原理的射線追蹤法。
1.1 時域有限差分（FDTD）法簡介
時域有限差分（FDTD）法由K.S.Yee[12]於1966年提出，因為其計算過程直觀、通用性好而得到廣泛應用。時域有限差分法是利用離散的二階中心差分方程近似替代連續的Maxwell微分方程。因此，構建合理的差分格式是保證微分方程解準確性的關鍵步驟。同時，又由於離散化取樣的過程會導致電磁波的頻散現象的發生，在FDTD計算過程中，時間步長Δt與空間步長Δl的取值必須滿足一定的限定條件，否則差分方程解將會不穩定。當電磁波在介質中傳播時，電磁場的邊界總被看作是開放的，電磁場分佈空間被看作是無限的。而時域有限差分是透過計算機對電磁場的傳播及空間分佈進行數值模擬，由於計算機記憶體是有限的，所以其差分網格必然會被截斷。為了使電磁波不在差分網格截斷處發生明顯反射，就必須設定吸收邊界條件以確保計算機模擬的有限空間與自然界無限大空間之間的差異儘量的小[13]。
1.2 射線追蹤法簡介
射線追蹤法已廣泛應用於波傳播的正演模擬研究。射線追蹤方法計算速度快，結果直觀，在層析成像、偏移、反演及模型試算中都佔有很重要的地位。射線追蹤方法首先是在地震勘探和研究中發展很快，並取得了實際的進展。由於GPR理論和地震勘探理論的相似性，因而近年來也有學者將射線追蹤方法應用於GPR波的傳播和正演研究，如Cai（1995）應用射線追蹤方法進行了二維介質中雷達波的傳播與模擬研究[9]。射線追蹤即對射線路徑進行準確計算和描述，對於GPR正演模擬，射線路徑是指雷達波的傳播軌跡。射線追蹤法的理論基礎是在高頻近似條件下電磁波的主能量沿射線軌跡傳播。試射法是一種基於斯奈爾定律的射線追蹤法，其基本思路是：射線以某一給定的初始發射角從發射點射出，在介質中波的傳播遵循斯奈爾定律，記錄射線的傳播軌跡，直至到達指定的接收點附近滿足一定精度的區域內；若射線終止位置不滿足該精度，則調整初始出射角，重新計算直至射線滿足精度為止。
基於射線追蹤法的GPR數值模擬的核心是射線追蹤，射線路徑越準確，模擬出的結果才更接近於真實情況。基於試射法的GPR數值模擬實現主要有以下幾個步驟[15]：
（1）建立管道模型；（2）用一組射線對模型進行試射，根據斯奈爾定律求出各條射線的終點，找出滿足一定精度要求的射線的入射角的大致範圍；（3）在區間內應用二分法進行迭代，使射線滿足一定的精度要求；（4）根據射線與各個介面的交點和各個層位的物性引數，求出射線的旅行時間和各個介面的反射波幅度；（5）用已求出的反射波幅度的時間序列與脈衝子波進行褶積，即得出單道雷達記錄；將各道雷達記錄按順序排列，便得出模擬雷達剖面。因此，從理論上講，圓柱狀管道GPR時—距曲線的張開度（離心率）的大小僅與管道周圍介質雷達波速度有關，即波速值越大，雙曲線的離心率越大，雙曲線張開度越大。由於GPR管道探測總是以地面位置為參考來分析和判斷來自管道回波同相軸的變化，因此就有可能把由於管道直徑的變化所引起的回波同相軸張開度的變化誤判為管道周圍介質中雷達速度變化引起的，或者做出相反的判斷。
2、地下管道模型GPR數值模擬與物理模擬特徵分析
為了對地下管道探測GPR響應特徵進行模擬研究，必須建立包含管道的各種地電模型，這必然涉及到各種管道的材質、管道周圍介質以及管道內充填介質等。
2.1 不同材質與管內充填不同性質流體的管道GPR響應模擬分析
圖1a、b、c和d為採用EKKO100型GPR系統採用主頻為225ZHz天線對金屬與非金屬管道進行物理模擬探測的GPR記錄（原始資料由加拿大Sensor&.Software公司提供）。管道埋設在電性不均勻的砂質土壤中，圖1a為直徑為15.5cm充滿空氣的金屬管道；圖2b、c、d分別為充滿自來水（電阻率為34.5Ω·m）、空氣以及鹹水（電阻率為0.212Ω·m）直徑為16cm的塑膠管道。所有管道頂面埋深均為1.25m。圖2a表明，由於金屬材質的高反射率與強吸收衰減作用，GPR剖面只能看到來自管道頂面的回波同相軸，管頂回波幅度強，雙曲線形同相軸對稱性好。圖1b顯示充滿自來水的塑膠管道的GPR剖面上不但能看到管道頂、底的清晰回波，還能隱約看到來自管道頂、底之間的多次波同相軸。圖1c顯示充滿空氣的塑膠管道由於來自管道頂、底回波時差很小（由於電磁波在空氣中的速度約為水中速度的9倍，因此在圖1c中的GPR記錄中不但看不到來自管底的強反射回波，而且由於來自管底回波與管頂回波部分重疊而產生相消干涉而導致GPR記錄上管頂回波幅度明顯減弱了。圖1d顯示充滿高導電率的鹹水管道由於管內鹹水的強吸收衰減作用GPR記錄上只能看到來自管頂的回波同相軸。
圖1採用225MHz天線0.5m天線距對金屬與非金屬管道進行物理模擬的記錄。（a）、充滿自來水的塑膠管道（b）、充滿空氣的塑膠管道（c）及充滿鹹水的塑膠管道（d）進行物理模擬的GPR記錄。圖（b）中T為管頂反射，B為管底反射，M為管內頂、底之間的多次波。
2.2 周圍介質屬性對管道GPR響應特徵的影響模擬分析
管道一般埋設在土壤等介質中，管道周圍介質屬性不同，GPR管道探測的響應特徵也會不同。為了考察周圍介質屬性變化對管道GPR探測響應特徵的影響，建立以下管道模型：管道材質為混凝土，管徑1.2m，管道周圍介質分別為淤泥、粉質黏土和自來水；管頂埋深1.0m，管壁厚度0.05m；雷達激勵源子波中心頻率為200MHz。採用射線追蹤法進行模擬。
當管道埋設在粉質粘土中時，由於粉質黏土導電率相對較低，模擬記錄上管道頂、底回波同相軸幅度強；當管道埋設在在淤泥中時，由於淤泥的導電率較高，對雷達波具有強吸收衰減，管道頂、底回波能量弱，同相軸幅度低（圖8b）；當管道位於在自來水中時，由於自來水的導電率介於粉質黏土與淤泥之間，加上水的介電常數值最大，雷達波在其中的傳播速度最低，因此，管頂回波同相軸幅度較淤泥介質中的明顯，但比粉質黏土中的明顯偏低。另外，管道周圍介質雷達波速度越低，管頂回波同相軸形狀越偏狹，張開度越小。
3、討論與結論
基於單道反射的GPR剖面上圓柱狀管道的回波時—距曲線為雙曲線，該曲線的離心率及其漸近線的斜率與管道半徑大小無關，只與管道周圍介質中雷達波速度有關，漸近線斜率與管道周圍介質雷達波速度成反比，即管道周圍介質中雷達波速度越高，管道回波時—距曲線張開度越大。因此，管道周圍介質的電磁波速度是影響管道回波同相軸形態的主要因素。管徑越大，管道回波同相軸雙曲線張開度越大，管頂和管底回波時差越大，回波同相軸分離越開；反之，當管道半徑足夠小時，管頂和管底反射同相軸趨於重合。透過對GPR記錄上來自圓柱狀管道回波同相軸的擬合可以確定其時—距曲線即雙曲線方程，並由此可以反演管道周圍介質的平均雷達波速度及管道半徑等數，根據波速可以確定管頂埋深。
對不同材質與管內充填不同性質流體的管道GPR響應物理模擬與數值模擬結果表明：由於金屬材質的高反射率與強吸收衰減作用，金屬管道GPR剖面只有來自管頂的回波同相軸，回波幅度強，雙曲線形同相軸對稱性好；充滿自來水的塑膠管道的GPR剖面上不但有來自管道頂、底的清晰回波存在，還能隱約看到來自管道頂、底之間的多次波同相軸；充滿空氣的塑膠管道由於來自管道頂、底回波時差很小，因此GPR記錄中不但看不到來自管底的強反射回波，而且由於來自管底回波與管頂回波部分重疊而產生相消干涉而導致GPR記錄上管頂回波幅度明顯減弱；充滿高導電率的鹹水管道由於管內鹹水的強吸收衰減作用GPR記錄上只能看到來自管頂的回波同相軸。在潮溼黏土介質中隨管頂埋深的增加，雷達波受吸收衰減及球面擴散的影響，管道回波能量逐漸減弱；反射回波同相軸雙曲線張開度也隨管頂埋深加大而增大。管道周圍介質的電性不但影響管道回波同相軸的形態，也會明顯影響其幅度；管道周圍介質導電率越高，對雷達波的吸收衰減越強，來自管道頂、底回波能量約弱，同相軸幅度越低。
