基于FDTD的空間步進頻雷達煤層異質(zhì)區(qū)探測正演分析

賀 順，羅彩禎，楊志偉，賈遠航

（1.西安科技大學 通信與信息工程學院，陜西 西安 710054；2.西安電子科技大學 雷達信號處理重點實驗室，陜西 西安 710071）

0 引 言

現(xiàn)如今由煤層異質(zhì)區(qū)如采空區(qū)、富水區(qū)、瓦斯突出區(qū)等引發(fā)的礦井事故多有發(fā)生，了解掘進面前方煤體的結構對煤礦安全生產(chǎn)尤為重要。目前，多種地球物理勘探技術已被廣泛應用于煤礦井下地質(zhì)預測中，如礦井直流電法、瞬變電磁法、地震勘探法等。礦井直流法受地形的影響較大，需要花費大量的時間進行設備的鋪設；瞬變電磁法存在淺部探測盲區(qū)且容易受到鐵磁介質(zhì)的影響；在多種勘探技術中，探地雷達有抗干擾能力強、受地形的影響較小、操作方便、不需要進行龐雜的設備布置等優(yōu)點。許多學者對探地雷達應用于煤礦井下探測進行了研究。文獻[2?3]使用地質(zhì)雷達探測井下瓦斯富集區(qū)并對其特征進行研究；文獻[4]使用GprMax和Matlab數(shù)值軟件對井下二維空間磁場模型進行正演模擬，研究不同激勵源與中心頻率對模擬結果的影響。文獻[5]使用超寬帶探地雷達開展了煤層厚度的探測試驗。井下探測使用的探地雷達多為單發(fā)單收雷達系統(tǒng)，這種系統(tǒng)雖然操作簡單，探測方便，但需要多次移動進行多點測量，可獲得的信息量有限。

多輸入多輸出探地雷達探測一次可獲得多組數(shù)據(jù)，相比單站雷達能獲得更多的信息量，更加充分地利用了回波信號的空時信息，并提高了系統(tǒng)的空間分辨能力。文獻[6?7]先后對陣列式探地雷達目標成像問題以及陣列式探地雷達信號極化場特征進行了研究；文獻[8]對車載陣列探地雷達實時檢測地雷進行了研究。目前，多輸入多輸出探地雷達多應用于檢測地下藏雷、管線結構、路面地下空洞等方面，在煤礦井下探測的應用較少。

探地雷達通過增大信號帶寬和天線合成孔徑來提高雷達探測的距離分辨力和方位分辨力。硬件發(fā)展的水平限制了超寬帶信號的產(chǎn)生。多輸入多輸出體制下的空間步進頻雷達通過在相鄰的陣元上依次發(fā)射頻率步進的信號，在接收端合成一個寬帶信號，解決了探地雷達系統(tǒng)中目標的距離分辨力與雷達信號帶寬之間的矛盾，同時相比單站雷達方位分辨力又有所提高。本文主要研究多輸入多輸出雷達體制下的空間步進頻信號在煤礦井下環(huán)境的傳播，建立二維煤礦井下異質(zhì)區(qū)環(huán)境模型，使用時域有限差分方法（Finite?difference Time?domain，F(xiàn)DTD）進行空間步進頻信號探測異質(zhì)區(qū)的數(shù)值模擬，分析不同異質(zhì)區(qū)模型下的回波結果，得到對應模型下的傳播規(guī)律。

1 空間步進頻信號模型

空間步進頻雷達對煤礦井下異質(zhì)區(qū)進行探測時，將雷達天線陣列放置在井下巷道中。天線陣列同時平行于巷道地面與煤礦開采掘進面，如圖1所示。發(fā)射天線向掘進面發(fā)射電磁波信號，信號在電性參數(shù)不同的物質(zhì)分界面上發(fā)生反射和透射現(xiàn)象，接收天線接收反射信號，通過分析接收信號的時間、相位等信息對異質(zhì)區(qū)的位置等特征進行解釋。

圖1 空間步進頻雷達井下放置示意圖

電磁波信號在介質(zhì)中傳播的快慢，影響回波響應出現(xiàn)的時間，根據(jù)相鄰回波響應之間的時間差，可計算電磁波傳播的距離。雷達發(fā)射信號在介質(zhì)中的傳播速度可根據(jù)式（1）求得：

式中：為電磁波在自由空間傳播的速度；′為介質(zhì)的相對介電常數(shù)?？芍殡姵?shù)越大，信號在介質(zhì)中的傳播速度越慢。

空間步進頻雷達實現(xiàn)空間上的頻率步進信號，使用等間隔的均勻線陣組成線陣的每個陣元發(fā)射頻率步進的信號?？臻g步進頻雷達發(fā)射陣列結構如圖2所示，發(fā)射天線為個陣元的均勻線陣，其陣元間距為。第個陣元的發(fā)射信號x ()為：

圖2 空間步進頻雷達發(fā)射陣列結構

（1）空間位置管理需要采用矢量數(shù)據(jù)結構。在地理位置表述時，矢量數(shù)據(jù)具有精確度高、數(shù)據(jù)量小的特性，以方便全面系統(tǒng)的闡述空間位置關系，對于數(shù)據(jù)圖形、屬性及數(shù)據(jù)恢復均能在短時間內(nèi)實現(xiàn)，其具體格式如圖4。

2 FDTD時域有限差分原理

井下探地雷達探測異質(zhì)區(qū)屬于電磁類探測方法。FDTD十分適用于探地雷達建模，因為它能有效地模擬具有不同幾何形狀和介電特性的地下分層結構。FDTD方法將微分形式的時域麥克斯韋方程轉(zhuǎn)化為差分方程，使用迭代的方法求解有限空間、時間范圍內(nèi)的電場和磁場分量。一般情況下，對于非均勻、各項同性、線性的媒質(zhì)材料，其麥克斯韋電磁方程組微分形式如下所示：

式中：表示電場強度矢量，單位為V/m；表示磁場強度矢量，單位為A/m；表示電位移矢量，單位為C/m；表示磁感應強度，單位為T；表示電流密度矢量，單位為A/m；J 表示磁通量密度矢量，單位為V/m；表示電荷密度，單位為C/m。

同時有電磁場對材料介質(zhì)作用的關系式，即本構方程：

式中：表示物質(zhì)材料的介電常數(shù)，單位為F/m；表示物質(zhì)材料的電導率，單位為S/m；表示物質(zhì)材料的磁導率，單位為H/m；σ表示計算磁損耗的磁阻率，單位為Ω/m。

FDTD方法對計算空間進行網(wǎng)格的剖分，圖3為著名的Yee元胞，元胞在坐標軸方向上的大小稱為該方向上的空間步長，規(guī)定每一個（或）場分量周圍有4個（或）場分量環(huán)繞，應用這種離散方式將含時間變量的麥克斯韋旋度方程轉(zhuǎn)化為差分方程。

圖3 FDTD離散中的元胞

對于二維問題，所有物理量與無關，??=0。以TM波為例，其麥克斯韋旋度方程為式（8），式（9）~式（11）為推導所得任意節(jié)點處電場和磁場的FDTD差分公式，式中分別為（，+1 2），（+1 2，），（，）。

由于使用差分方法對麥克斯韋方程進行數(shù)值計算會產(chǎn)生色散現(xiàn)象，即波的傳播速度會隨波長而改變。這種色散會引起波形的畸變，且網(wǎng)格越粗色散越嚴重。為了使差分方程的計算結果能更近似于原方程的解，需要對時間步長和空間步長進行一定的限制，在二維問題中應該滿足式（12）。

3 正演模擬

3.1 實驗模型及參數(shù)設置

在煤礦開采的過程中，地質(zhì)條件復雜多變，將存在于煤體中的非煤區(qū)域稱為異質(zhì)區(qū)。煤層中常見的異質(zhì)區(qū)有陷落柱、斷層、采空區(qū)、富水區(qū)、瓦斯富集區(qū)等。本文對采空區(qū)、富水區(qū)以及巖層三種主要煤礦異質(zhì)區(qū)進行模擬。未考慮巷道中各種儀器、線纜等產(chǎn)生的干擾，并假設介質(zhì)為均勻同性，建立掘進工作方向煤層內(nèi)部的物理模型如圖4所示。模型設置異質(zhì)區(qū)上方煤層厚度為0.5 m，采空區(qū)厚度為0.5 m，富水區(qū)厚度為0.5 m，薄泥巖層厚度為0.5 m。

圖4 異質(zhì)區(qū)物理模型

本文使用收發(fā)分置的天線陣，收發(fā)天線均為5陣元的均勻線陣，陣元間隔為0.1 m。收發(fā)天線陣水平放置，與掘進方向煤層表面平行。收發(fā)天線位于同一水平面，發(fā)射天線陣放置距煤層表面0.8 m，接收天線陣放置距煤層表面0.6 m。發(fā)射天線發(fā)射時長為1 ns的正弦信號，選定線陣一端的陣元作為第一陣元，其所發(fā)射信號的中心頻率為200 MHz，其余陣元的中心頻率依次以200 MHz的頻率增量進行步進，這樣第5個陣元發(fā)射信號的中心頻率為1 GHz。

模擬中的最高頻率選擇為陣元中最高中心頻率的3倍，由最高頻率確定的最小波長為0.1 m，計算得到空間步長不應超過0.01 m，本文實驗中3個二維異質(zhì)區(qū)模型在方向與方向上的空間網(wǎng)格步長均設定為0.003 m。由式（12）可算得時間步長為7.076 0×10s。本文使用完美匹配層（Perfectly Matched Layer，PML）作為吸收邊界條件對空間進行截斷。異質(zhì)區(qū)模型中所用介質(zhì)的電性參數(shù)如表1所示。

表1 模型電性參數(shù)

3.2 傳播特性分析

不改變發(fā)射天線以及接收天線的位置，對上述三種異質(zhì)區(qū)進行數(shù)值計算，得到接收天線5個陣元處隨時間迭代次數(shù)變化的電場信號數(shù)值。將5個接收點結果進行繪圖顯示，得到三種異質(zhì)區(qū)的仿真結果如圖5~圖7所示。

圖5 采空區(qū)結果

圖7 薄泥巖層結果

圖5為采空區(qū)模擬結果，煤與空氣界面表現(xiàn)為一條明顯的軸線，出現(xiàn)在第1 157次時間迭代處，通過計算可知該點時間為8.187 ns。采空區(qū)上下表面回波對應的反射軸線分別出現(xiàn)在第2 337次與第2 789次時間迭代處，計算得到對應時間分別為16.537 ns和19.735 ns。某位置的時間等于該位置的時間迭代次數(shù)乘以時間步進間隔，時間步進間隔可由式（12）求得。

圖6為富水區(qū)模型的模擬結果，可以清楚地看到在煤與空氣界面形成的響應后僅有一條清楚的響應軸線，該反射軸線表示煤層與富水區(qū)上表面的回波。由于水的介電常數(shù)大，電磁波在水中的衰減速度快。當信號穿過煤層進入富水區(qū)之后很快被吸收，因此僅能看到煤與水上表面形成的回波響應。

圖6 富水區(qū)結果

圖7為薄泥巖層模型的模擬結果，與采空區(qū)模型相似，同樣可以看到煤與泥巖上下交界面處的兩個回波，但與采空區(qū)相比較為微弱。

實驗設置天線陣平行于掘進煤層表面，即發(fā)射信號為垂直入射，根據(jù)雙程走時′和速度可計算煤體到異質(zhì)區(qū)上表面厚度=′2，電磁波在煤中的傳播速度利用式（1）求得。表2為根據(jù)實驗數(shù)據(jù)對模型厚度進行計算的結果，表3為對模型內(nèi)介質(zhì)速度進行計算的結果。速度從大到小排序為空氣、煤、泥巖、水。根據(jù)仿真數(shù)據(jù)得到的計算結果與理論值存在一定的誤差，厚度的誤差不超過0.05 m，速度的誤差不超過8%。除仿真過程中產(chǎn)生的誤差外，還有結果讀取和計算時不精確造成的誤差。

表2 模型厚度分析 m

表3 模型介質(zhì)速度分析

異質(zhì)區(qū)的物理模型不變，改變異質(zhì)區(qū)上方煤層厚度為5 m，采空區(qū)、富水區(qū)、泥巖區(qū)三種異質(zhì)區(qū)的厚度均為5 m，和方向的空間步長調(diào)整為0.01 m，其余參數(shù)不變，得到的結果如圖8~圖10所示。同樣為采空區(qū)與泥巖區(qū)表現(xiàn)為兩條響應軸線，泥巖區(qū)的響應較弱，富水區(qū)表現(xiàn)為一條響應軸線，與厚度為0.5 m模型下的模擬結果一致，僅到達時間不同。

圖8 5 m厚采空區(qū)模擬結果

圖9 5 m厚富水區(qū)模擬結果

圖10 5 m厚泥巖區(qū)模擬結果

為了更好地觀察空間步進頻信號在整個仿真區(qū)域的傳播情況，計算了三種模型空間中所有網(wǎng)格的電場強度，繪制頻率為1 GHz情況下二維仿真模型平面上的電場分布圖，如圖11所示。

圖11 三種模型的電場分布圖

電場分布圖中的橫縱坐標分別為模型掘進方向與垂直掘進方向上的距離坐標。設置=0的界面為煤礦開采工作面，橫坐標的正半軸表示掘進面上方的巷道；橫坐標的負半軸表示掘進面下方的煤層內(nèi)部。

右側(cè)各發(fā)射陣元發(fā)射頻率步進的信號向左側(cè)煤層中進行傳播。圖11a）為采空區(qū)模型下的電場分布情況，電磁波信號由自由空間傳入煤層中，采空區(qū)部分波紋與煤層上表面自由空間處波紋一致；圖11b）為富水區(qū)的電場分布情況，由于能量衰減快在傳播中消耗殆盡，電磁波基本不能穿過水層；圖11c）為薄泥巖層的電場分布情況，電磁波從煤層到泥巖中波紋變密。電磁波信號傳播過程中波紋有密有疏，波紋的稀疏情況隨波長大小發(fā)生變化。在頻率不變的情況下，波速越大波長就越大，波紋越稀疏，根據(jù)電場分布圖得到模型中介質(zhì)速度由大到小排序分別為空氣、煤、泥巖、水，與理論一致。

4 結 語

本文基于時域有限差分方法對空間步進頻信號在三種不同井下異質(zhì)區(qū)模型中的傳播情況進行正演模擬，模擬結果表明：通過使用空間步進頻信號的時域波形進行數(shù)值計算，得到的煤層厚度與正演模型設置的參數(shù)基本吻合；在三種異質(zhì)區(qū)正演結果中，采空區(qū)模型與薄泥巖層模型下均可看到信號在異質(zhì)區(qū)上下表面產(chǎn)生的兩條反射軸線，但采空區(qū)的回波更為明顯，而富水區(qū)模型下僅可看到水煤上交界面的反射軸線；根據(jù)模擬的電場分布圖，通過波紋疏密程度可以對介質(zhì)的速度大小進行大致的排序；信號傳播時衰減的快慢不同，在水中的衰減最大，泥巖次之，空氣衰減最小。

由以上仿真結果可知：進行實際探測時，若探測數(shù)據(jù)出現(xiàn)明顯的回波異常區(qū)域，則可以判斷該處存在與煤體不同的異質(zhì)區(qū)，并可根據(jù)異?；夭ǔ霈F(xiàn)的時間點對異質(zhì)區(qū)域上方的煤層厚度進行估計；若信號在產(chǎn)生響應后急速衰減則可判斷該區(qū)域存在水，驗證了使用多發(fā)多收體制下的空間步進頻雷達對煤體異質(zhì)區(qū)探測的可行性。