王文天，劉四新，鹿 琪，李宏卿，傅 磊

吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院，長(zhǎng)春 130026

0 引言

探地雷達(dá)是以地下介質(zhì)的電磁特性差異為基礎(chǔ)的地球物理探測(cè)技術(shù)，其通過發(fā)射天線發(fā)送脈沖波形式的高頻電磁波[1]。電磁波在地下介質(zhì)傳播過程中，遇到存在電磁特性差異的界面時(shí)便發(fā)生反射，返回的電磁波被接收后，通過對(duì)信號(hào)進(jìn)行處理，形成雷達(dá)時(shí)間剖面圖[2]；然后根據(jù)圖中信號(hào)強(qiáng)度、波形、雙程走時(shí)等參數(shù)分析探測(cè)目標(biāo)體的空間位置、電性狀況與幾何形態(tài)等[3]。按照工作方式的不同，探地雷達(dá)分為地面探地雷達(dá)和鉆孔雷達(dá)[4]。鉆孔雷達(dá)通過將雷達(dá)天線置于鉆孔中進(jìn)行探測(cè)，天線更接近地下目標(biāo)體，具有更大的探測(cè)范圍和分辨率[5]。相對(duì)于地面雷達(dá)，鉆孔雷達(dá)儀器在井中的探測(cè)范圍得到較大的擴(kuò)展,在相對(duì)導(dǎo)電的巖石中,探測(cè)范圍可達(dá)10～20 m[6]。鉆孔雷達(dá)的發(fā)展歷程大致如下：1978年，Rubin等[7]在輝綠巖中進(jìn)行了鉆孔雷達(dá)實(shí)驗(yàn)；1980年國(guó)際STRIPA(international stripa project)計(jì)劃開始實(shí)行,瑞典RAMAC鉆孔雷達(dá)系統(tǒng)則是在該計(jì)劃下開發(fā)完成的[8]；20世紀(jì)90年代后,日本東北大學(xué)Sato 等[9]研究開發(fā)了極化鉆孔雷達(dá)系統(tǒng),能夠進(jìn)行井下的全極化測(cè)量。

傳統(tǒng)的鉆孔雷達(dá)發(fā)射天線和接收天都是全方位的[10]，因此，其只能探測(cè)目標(biāo)地質(zhì)體深度以及其距井孔的距離，并不能對(duì)其方位進(jìn)行探測(cè)。目前，對(duì)地下介質(zhì)進(jìn)行三維立體成像仍然是一件非常困難的事情。為了解決這個(gè)問題，發(fā)展了定向鉆孔雷達(dá)技術(shù)。

定向鉆孔雷達(dá)不僅可以探測(cè)目標(biāo)的深度和距井孔距離，而且可以探測(cè)其方位。定向鉆孔雷達(dá)的實(shí)現(xiàn)一般有兩種方式：其一是用陣列天線接收的方式進(jìn)行測(cè)量，主要使用均勻圓陣，根據(jù)不同的波至?xí)r間提取方位信息；其二是采用定向發(fā)射天線的方式，天線向井周定向發(fā)射電磁波，目標(biāo)體的方位位于電磁波的發(fā)射方向[11]。

對(duì)于以上兩種定向鉆孔雷達(dá)系統(tǒng)，人們分別開發(fā)了相應(yīng)的數(shù)據(jù)處理方法。對(duì)于第一種定向鉆孔雷達(dá)系統(tǒng)，即用陣列天線接收的定向鉆孔雷達(dá)系統(tǒng)，采用分析各個(gè)接收天線所接收到波形的先后次序來分析反射波的方位信息；對(duì)于第二種定向鉆孔雷達(dá)系統(tǒng)，天線以一定的角度發(fā)射電磁波，如果該角度下接收天線可以接收到信號(hào)，說明該方位存在目標(biāo)地質(zhì)體，否則，該方位不存在目標(biāo)地質(zhì)體[12]。

本文針對(duì)的是第一種定向鉆孔雷達(dá)系統(tǒng)的實(shí)現(xiàn)方式，包括1個(gè)發(fā)射天線和4個(gè)接收天線，4個(gè)天線分布于一個(gè)垂直于井眼的圓環(huán)位置上，呈等角度分布(圖1)。為了分析方便，我們把4個(gè)天線指定為東、南、西、北4個(gè)方向，順時(shí)針旋轉(zhuǎn)。當(dāng)北天線方向和地理正北一致時(shí)，其余天線分別和東、西、南一致，否則會(huì)有偏差角，這可以通過測(cè)量北天線相對(duì)于地理正北的夾角來校正。針對(duì)這種系統(tǒng)，2000年，Ebihara等[13]引入MUSIC算法進(jìn)行數(shù)據(jù)處理，取得了較好的效果，同時(shí)確定了MUSIC算法在低信噪比情況下的局限性。2007年， Takayama等[14]利用反正切法進(jìn)行數(shù)據(jù)處理。2013年，Ebihara等[15]提出了使用計(jì)算殘差來確定目標(biāo)方位的算法，該算法對(duì)信號(hào)由信源到達(dá)圓陣中心的時(shí)間和入射角進(jìn)行殘差法擬合，來求取入射方位，但只能對(duì)一個(gè)信源信號(hào)進(jìn)行方位識(shí)別處理。

本文改進(jìn)殘差法算法，采用設(shè)計(jì)窗口逐點(diǎn)移動(dòng)的辦法，直接計(jì)算信號(hào)由信源到達(dá)圓陣中心時(shí)間，僅對(duì)方位角進(jìn)行殘差法擬合。這樣，一方面大幅提高了算法的運(yùn)算速度，另一方面，可以對(duì)任意多個(gè)信源信號(hào)進(jìn)行處理。同時(shí)，本文提出并實(shí)現(xiàn)了利用方位角數(shù)據(jù)合成橫向切片和縱向切片的算法，再利用所得的縱橫切片合成三維數(shù)組,再利用所得的三維數(shù)組對(duì)地質(zhì)體進(jìn)行三維成像的算法，獲得三維成像結(jié)果。

1 改進(jìn)的殘差法算法

Ebihara等[15]提出了使用殘差法對(duì)入射信號(hào)進(jìn)行方位估計(jì)的算法，主要是針對(duì)在一個(gè)鉆孔中同軸線饋電的偶極天線陣列，如圖2所示。算法包括如下假設(shè)：只有一個(gè)入射信號(hào)入射到天線陣元中；各個(gè)天線陣元接收到的信號(hào)不發(fā)生畸變；各個(gè)陣元的噪聲水平是相同的。在選取方位角時(shí)，取正北方向?yàn)?°，順時(shí)針方向?yàn)檎较颍虼?，第i個(gè)天線陣元接收到信號(hào)的時(shí)間ti可以表示為

Tx. 發(fā)射天線；Rx. 接收天線。圖1 單發(fā)-四收定向鉆孔雷達(dá)示意圖Fig.1 Single-transmitter four-receiver borehole rada

a.各方位接收天線接收信號(hào)；b.北向接收天線接收信號(hào);c.改變?chǔ)群蛅北向接收天線接收信號(hào)；d.南向接收天線接收信號(hào)。圖2 殘差法波形示意圖Fig.2 Waveform of residual algorithm

ti(θ,t)=t-τcos(φi-θ)。

(1)

式中：t為波傳到圓陣中心所用的時(shí)間；τ為電磁波在兩個(gè)天線陣元之間最大傳播時(shí)間的一半，即如果圓陣半徑為r，電磁波的傳播速度為v，那么τ=r/v；φi為第i個(gè)天線陣元所在的方位角；θ為入射波的方位角。將第i個(gè)天線陣元所接收到的實(shí)測(cè)數(shù)據(jù)定義為ti，每一個(gè)t和θ的改變值對(duì)應(yīng)數(shù)據(jù)信號(hào)為fi。假設(shè)共有N個(gè)采樣點(diǎn)，接收天線的數(shù)目為M，在本文中，由于存在4個(gè)接收天線，因此M為4。通過一定范圍內(nèi)對(duì)于t和θ的改變，找到使得測(cè)量值與計(jì)算值的殘差達(dá)到最小的情況，即使得殘差

(2)

最小。式中：tij表示第i個(gè)接收天線在第j個(gè)采樣點(diǎn)的實(shí)測(cè)時(shí)間數(shù)據(jù)；fij表示每一個(gè)v和θ的改變值在第i個(gè)接收天線在第j個(gè)采樣點(diǎn)對(duì)應(yīng)的擬合信號(hào)。最小殘差所對(duì)應(yīng)的t和θ組合即為所求結(jié)果。可以看出，在一個(gè)深度，只有一個(gè)方位。

本文改進(jìn)該算法，通過選取寬度為w的窗口(圖3)，只對(duì)窗口中的信號(hào)求取殘差，這樣，可以對(duì)任意多信源進(jìn)行處理，而之前的殘差法只能對(duì)一個(gè)信源信號(hào)進(jìn)行處理。將窗口內(nèi)求取的方位角賦值給φk，公式為

φk=φrk。

(3)

其中：φk為窗口中央點(diǎn)的方位角；k為該道窗口中央的采樣點(diǎn)號(hào)；φrk為在一定的窗口下利用殘差法求取方位角，即

(4)

本方法使得t可以通過窗口中心的走時(shí)直接計(jì)算得出，即t=kdt，其中dt為時(shí)間采樣間隔。而之前的殘差法需要對(duì)t、θ兩個(gè)變量進(jìn)行殘差擬合求出。

先假設(shè)入射角為1°，2°，…，360°，分別計(jì)算各個(gè)天線的接收數(shù)據(jù)；然后用實(shí)測(cè)數(shù)據(jù)與理論數(shù)據(jù)的差計(jì)算絕對(duì)值，殘差最小所對(duì)應(yīng)的方位角即為目標(biāo)方位角。

a.各方位接收天線接收信號(hào)；b.窗口內(nèi)北向接收天線接收信號(hào);c.改變?chǔ)戎蟮谋毕蚪邮仗炀€接收信號(hào)；d.南向接收天線接收信號(hào)。圖3 改進(jìn)的殘差法示意圖Fig.3 Theory of improved residual algorithm

殘差法的具體算法如下:

1) 輸入圓陣半徑r、采樣時(shí)間間隔dt等參數(shù)，東、南、西、北4個(gè)接收天線所接收到的信號(hào)為TEm、TSm、TWm、TNm。

2) 選取一定大小的窗口，對(duì)接收到的信號(hào)進(jìn)行處理。數(shù)據(jù)處理前要設(shè)定一個(gè)閾值，當(dāng)信號(hào)大于閾值才算作有用信號(hào)，否則視為噪聲。

3) 如圖4所示，選取北向接收天線作為理論計(jì)算的標(biāo)準(zhǔn)。假設(shè)入射角為θ，陣列半徑為r，那么北向接收天線與東向接收天線的時(shí)間差為

ΔtE=rcos(θ)/v-rsin(θ)/v。

(5)

時(shí)間差所對(duì)應(yīng)的采樣點(diǎn)數(shù)為ΔnE=ΔtE/dt。

圖4 殘差法計(jì)算示意圖Fig.4 Theory of residual algorithm

于是,東向接收天線的計(jì)算數(shù)據(jù)為

TEc=TNm(t-ΔnEdt)。

(6)

北向接收天線與南向接收天線的時(shí)間差為

ΔtS=2rcos(θ)/v，

(7)

時(shí)間差所對(duì)應(yīng)的采樣點(diǎn)數(shù)為ΔnS=ΔtS/dt，于是，南向接收天線的計(jì)算數(shù)據(jù)為

TSc=TNm(t-ΔnSdt)。

(8)

西向接收天線與東向接收天線的時(shí)間差為

ΔtW=rcos(θ)/c+rsin(θ)/v，

(9)

時(shí)間差所對(duì)應(yīng)的采樣點(diǎn)數(shù):ΔnW=ΔtW/dt，于是，西向接收天線的計(jì)算數(shù)據(jù)為

TWc=TNm(t-ΔnWdt)。

(10)

4)最后計(jì)算殘差：

V=|TEc-TEm|+|TSc-TSm|+|TWc-TWm|。

(11)

最小殘差所對(duì)應(yīng)的角度即為方位角。

殘差法算法進(jìn)行方位識(shí)別的主要原理是利用4個(gè)接收天線所接收信號(hào)的相對(duì)位置計(jì)算方位角。一旦信噪比過低，導(dǎo)致4個(gè)接收天線接收到信號(hào)的相對(duì)位置發(fā)生較大的變化，必然影響方位識(shí)別結(jié)果。如果信號(hào)中存在少量噪聲，對(duì)于4個(gè)接收天線接收到信號(hào)的相對(duì)位置并無太大影響，那么殘差法算法依然可以得到非常好的計(jì)算結(jié)果。

本文以分離立交橋梁施工過程中施工風(fēng)險(xiǎn)為研究對(duì)象，運(yùn)用WBS-RBS耦合矩陣分析方法對(duì)整個(gè)分離立交橋梁施工中的工作進(jìn)行逐層分解，同時(shí)將項(xiàng)目過程中可能存在的風(fēng)險(xiǎn)因素逐層向下延伸，最終形成耦合矩陣.而后運(yùn)用粗糙集確定各風(fēng)險(xiǎn)權(quán)重，確定分析出各風(fēng)險(xiǎn)的重要程度及主次要大小，為分離立交橋梁施工的安全生產(chǎn)管理提供參考.

2 合成數(shù)據(jù)算例

用GprMax合成數(shù)據(jù)[16]，其模型包括2個(gè)裂縫，一個(gè)裂縫在正東方向，另一個(gè)裂縫在正南方向，具體如圖5所示。本模型包含2個(gè)目標(biāo)地質(zhì)體，由于傳統(tǒng)的殘差法不符合“只有一個(gè)入射信號(hào)入射到天線陣元中”的假設(shè)，因此是無法進(jìn)行方位識(shí)別的。利用該模型驗(yàn)證改進(jìn)的殘差法的有效性，模型參數(shù)如表1、表2所示，圍巖為花崗巖，裂縫內(nèi)為水。作為探地雷達(dá)的目標(biāo)探測(cè)而言，最重要的2項(xiàng)參數(shù)是相對(duì)介電常數(shù)和電導(dǎo)率。

表1 合成數(shù)據(jù)采用的材料性質(zhì)

模型對(duì)應(yīng)的波形如圖6所示，發(fā)射天線的起始深度為14.00 m，接收天線的起始深度為11.00 m，發(fā)射天線與接收天線的中間位置的起始深度為12.75 m，接收天線的4個(gè)陣元分別位于東、南、西、北4個(gè)方位，發(fā)射天線和接收天線均不考慮形狀，為點(diǎn)發(fā)射點(diǎn)接收，探地雷達(dá)由井底向井口移動(dòng)。從圖6可以清楚看到2個(gè)反射信號(hào)，距離井眼近的信號(hào)為正南方向的裂縫引起的反射波，距離井眼較遠(yuǎn)的信號(hào)為正東方向的裂縫所引起，與模型相符合。數(shù)據(jù)處理前，首先要去除直達(dá)波，將深度為12.25 m到12.75 m的五道前270個(gè)采樣點(diǎn)的平均波形作為直達(dá)波，從記錄中減去便可實(shí)現(xiàn)直達(dá)波的消除。

表2 合成數(shù)據(jù)采用的參數(shù)

圖5 用于鉆孔雷達(dá)合成數(shù)據(jù)的模型Fig.5 Model for the borehole radar synthetic data

利用殘差法計(jì)算方位角，選取11.75 m深度處的道的方位識(shí)別結(jié)果進(jìn)行顯示，如圖7所示。從圖7可以清楚看到：在50～100 ns處有一個(gè)信號(hào)，方位角為180°，即正南方向；在100～150 ns處有一個(gè)信號(hào)，方位角為90°，為正東方向，與模型吻合。

利用計(jì)算的方位角合成縱橫切片數(shù)據(jù)，基本流程如下：

1)設(shè)計(jì)一個(gè)三維數(shù)組，讓道數(shù)、方位角數(shù)和采樣點(diǎn)數(shù)為三維數(shù)組的3個(gè)維數(shù)。

2)根據(jù)各道的方位角數(shù)據(jù)和接收到的電場(chǎng)數(shù)據(jù)對(duì)三維數(shù)組進(jìn)行賦值。

4)對(duì)于該三維數(shù)組中的數(shù)據(jù)進(jìn)行顯示，得到該模型的三維成像結(jié)果圖。

圖6 各天線雷達(dá)信號(hào)波形Fig.6 Signals for four antennas

圖7 深度為11.75 m的方位識(shí)別結(jié)果圖Fig.7 Azimuth angle trace at the depth of 11.75 m

選取第11.75 m處的道，其橫向切片如圖8所示。從該切片圖中可以明顯看出正東方向存在一個(gè)

異常，正南方向存在一個(gè)異常，與模型相符合。90°縱向切片如圖9a所示，可以很明顯看到一個(gè)異常，因此和模型正東方向的裂縫吻合良好；180°縱向切片如圖9b所示，從圖中很明顯看到一個(gè)異常，與模型正南方向吻合良好。說明殘差法取得了非常好的效果。

圖8 第10道切片圖(深度11.75 m)Fig.8 Horizontal cross-section at the depth of 11.75 m

a.90°縱向切片；b. 180°縱向切片。圖9 縱向切片示意圖Fig.9 Vertical cross-section

最后，利用所求的方位角進(jìn)行三維成像顯示，結(jié)果如圖10所示。從圖10可以清楚看出，在井的正東方向存在一個(gè)異常，與原模型正東方向的裂縫吻合良好；在井的正南方向同樣存在一個(gè)異常，與原模型正南方向的裂縫吻合良好?？紤]到定向鉆孔雷達(dá)是從井底向井口運(yùn)行的，殘差法算法只能識(shí)別裂縫在豎直方向上的規(guī)模，類似地面物探中的二維測(cè)量。如果有一排鉆孔，就能看出水平向變化。綜上所述，殘差法在該模型中取得了非常好的應(yīng)用。

圖10 模型的三維成像結(jié)果Fig.10 Three-dimensional imaging results

3 結(jié)論

1)本文研究的定向鉆孔雷達(dá)包含1個(gè)發(fā)射天線和4個(gè)接收天線，通過所接收到波形到達(dá)的先后次序，可以有效地確定目標(biāo)地質(zhì)體的精確方位。

2)本文實(shí)現(xiàn)了對(duì)殘差法的算法改進(jìn)，成功實(shí)現(xiàn)了殘差法對(duì)于多目標(biāo)的方位識(shí)別，并在此基礎(chǔ)上提出并實(shí)現(xiàn)了鉆孔雷達(dá)縱橫切片以及三維成像方法。

3)利用GprMax合成數(shù)據(jù)，檢驗(yàn)改進(jìn)算法的適用性，可以看出改進(jìn)后殘差法可以有效判斷目標(biāo)地質(zhì)體的精確方位，并且在三維成像中取得了非常好的效果。

[1] 曾昭發(fā)，李文奔，習(xí)建軍，等. 基于DOA估計(jì)的陣列式探地雷達(dá)逆向投影目標(biāo)成像方法[J].吉林大學(xué)學(xué)報(bào)(地球科學(xué)版)，2017,47(4)：1308-1318.

Zeng Zhaofa, Li Wenben, Xi Jianjun, et al. Inverse Direction Imaging Method of Array Type GPR Based on DOA Estimation[J]. Journal of Jilin University(Earth Science Edition), 2017,47(4)：1308-1318.

[2] Sato M, Takayama T. A Novel Directional Borehole Radar System Using Optical Electric Field Sensors[J]. IEEE Transactions on Geoscience and Remote Sensing, 2007, 45(8):2529-2535.

[3] Lytle R J, Laine E F. Design of Miniature Directional Antenna for Geophysical Probing from Borehole[J]. IEEE Transactions on Geoscience and Remote Sensing, 1978, 16(6): 304-307.

[4] 朱自強(qiáng)，彭凌星，魯光銀，等. 基于互相關(guān)函數(shù)對(duì)鉆孔雷達(dá)層析成像的改進(jìn)[J].吉林大學(xué)學(xué)報(bào)(地球科學(xué)版)，2014 , 44 (2):668-674.

Zhu Ziqiang, Peng Lingxing, Lu Guangyin, et al. Improved Borehole-GPR Tomography Based on Cross-Correlation[J].Journal of Jilin University( Earth Science Edition),2014,44(2):668-674.

[5] 張理輕,馬曄,楊宇,等.鉆孔雷達(dá)數(shù)據(jù)處理技術(shù)及分析[J].地震工程學(xué)報(bào),2014,36(4):1107-1112.

Zhang Liqing, Ma Ye,Yang Yu, et al. Study on Data Processing Techniques of Borehole Radar[J]. China Earthquake Engineering Journal, 2014, 36(4):1107-1112.

[6] Lytle R J, Laine E F. Design of Miniature Directional Antenna for Geophysical Probing from Borehole[J]. IEEE Transactions on Geoscience and Remote Sensing, 1978, 16(6): 304-307.

[7] Olsson O, Falk L, Forslund O, et al. Borehole Radar Applied to the Characterization of Hydraulic Conductive Fracture Zones in Crystalline Rock[J]. Geophysical Prospecting, 2010, 40(2):109-142.

[8] 陳建勝,陳從新.鉆孔雷達(dá)技術(shù)的發(fā)展和現(xiàn)狀[J].地球物理學(xué)進(jìn)展,2008,23(5):300-306.

Chen Jiansheng, Chen Congxin. The Review of Borehole Radar Technology[J]. Progress in Geophysics, 2008, 23(5):300-306.

[9] Sato M, Ohkubo T, Niitsuma H.Cross-Polarimetric Borehole Radar Measurement with a Slot Antenna[J]. Journal of Applied Geophysics,1995,33(1):53-61.

[10] Ebihara S, Sasakura A. Mode Effect on Direct Wave in Single-Hole Borehole Radar[J]. IEEE Transactions on Geoscience and Remote Sensing, 2011, 49(2):854-867.

[11] Falk L R, Olsson O L, Sandberg E V. Combined Interpretation of Fracture Zones in Crystalline Rock Using Single-Hole, Crosshole Tomography and Directional Borehole-Radar Data[J]. Society of Petrophysicists and Well-Log Analysts,1991, 32(2):12.

[12] 曾昭發(fā),劉四新,王者江,等.探地雷達(dá)方法原理及應(yīng)用[M].北京：科學(xué)出版社,2006.

Zeng Zhaofa, Liu Sixin, Wang Zhejiang, et al. Ground Penetrating Radar Method and Application[M]. Beijing:Beijing Science Press, 2006.

[13] Ebihara S. Super-Resolution of Coherent Targets by a Directional Borehole Radar[J]. IEEE Transactions on Geoscience & Remote Sensing, 2000, 38(4):1725-1732.

[14]Takayama T, Sato M. A Novel Direction-Finding Algorithm for Directional Borehole Radar[J]. IEEE Transactions on Geoscience & Remote Sensing, 2007, 45(8):2520-2528.

[15] Ebihara S, Kawai H, Wada K. Estimating the 3D Position and Inclination of a Planar Interface with a Directional Borehole Radar[J]. Near Surface Geophysics, 2013, 11(2):185-195.

[16] Giannopoulos A. Modelling Ground Penetrating Ra-dar by GprMax[J]. Construction & Building Materials, 2005, 19:755-762.