基于改進殘差法的定向鉆孔雷達三維成像算法

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

吉林大學地球探測科學與技術學院，長春 130026

0 引言

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

傳統的鉆孔雷達發射天線和接收天都是全方位的[10]，因此，其只能探測目標地質體深度以及其距井孔的距離，并不能對其方位進行探測。目前，對地下介質進行三維立體成像仍然是一件非常困難的事情。為了解決這個問題，發展了定向鉆孔雷達技術。

定向鉆孔雷達不僅可以探測目標的深度和距井孔距離，而且可以探測其方位。定向鉆孔雷達的實現一般有兩種方式：其一是用陣列天線接收的方式進行測量，主要使用均勻圓陣，根據不同的波至時間提取方位信息；其二是采用定向發射天線的方式，天線向井周定向發射電磁波，目標體的方位位于電磁波的發射方向[11]。

對于以上兩種定向鉆孔雷達系統，人們分別開發了相應的數據處理方法。對于第一種定向鉆孔雷達系統，即用陣列天線接收的定向鉆孔雷達系統，采用分析各個接收天線所接收到波形的先后次序來分析反射波的方位信息；對于第二種定向鉆孔雷達系統，天線以一定的角度發射電磁波，如果該角度下接收天線可以接收到信號，說明該方位存在目標地質體，否則，該方位不存在目標地質體[12]。

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

本文改進殘差法算法，采用設計窗口逐點移動的辦法，直接計算信號由信源到達圓陣中心時間，僅對方位角進行殘差法擬合。這樣，一方面大幅提高了算法的運算速度，另一方面，可以對任意多個信源信號進行處理。同時，本文提出并實現了利用方位角數據合成橫向切片和縱向切片的算法，再利用所得的縱橫切片合成三維數組,再利用所得的三維數組對地質體進行三維成像的算法，獲得三維成像結果。

1 改進的殘差法算法

Ebihara等[15]提出了使用殘差法對入射信號進行方位估計的算法，主要是針對在一個鉆孔中同軸線饋電的偶極天線陣列，如圖2所示。算法包括如下假設：只有一個入射信號入射到天線陣元中；各個天線陣元接收到的信號不發生畸變；各個陣元的噪聲水平是相同的。在選取方位角時，取正北方向為0°，順時針方向為正方向，因此，第i個天線陣元接收到信號的時間ti可以表示為

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

a.各方位接收天線接收信號；b.北向接收天線接收信號;c.改變θ和t北向接收天線接收信號；d.南向接收天線接收信號。圖2 殘差法波形示意圖Fig.2 Waveform of residual algorithm

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

(1)

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

(2)

最小。式中：tij表示第i個接收天線在第j個采樣點的實測時間數據；fij表示每一個v和θ的改變值在第i個接收天線在第j個采樣點對應的擬合信號。最小殘差所對應的t和θ組合即為所求結果?？梢钥闯?，在一個深度，只有一個方位。

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

φk=φrk。

(3)

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

(4)

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

先假設入射角為1°，2°，…，360°，分別計算各個天線的接收數據；然后用實測數據與理論數據的差計算絕對值，殘差最小所對應的方位角即為目標方位角。

a.各方位接收天線接收信號；b.窗口內北向接收天線接收信號;c.改變θ之后的北向接收天線接收信號；d.南向接收天線接收信號。圖3 改進的殘差法示意圖Fig.3 Theory of improved residual algorithm

殘差法的具體算法如下:

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

2) 選取一定大小的窗口，對接收到的信號進行處理。數據處理前要設定一個閾值，當信號大于閾值才算作有用信號，否則視為噪聲。

3) 如圖4所示，選取北向接收天線作為理論計算的標準。假設入射角為θ，陣列半徑為r，那么北向接收天線與東向接收天線的時間差為

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

(5)

時間差所對應的采樣點數為ΔnE=ΔtE/dt。

圖4 殘差法計算示意圖Fig.4 Theory of residual algorithm

于是,東向接收天線的計算數據為

TEc=TNm(t-ΔnEdt)。

(6)

北向接收天線與南向接收天線的時間差為

ΔtS=2rcos(θ)/v，

(7)

時間差所對應的采樣點數為ΔnS=ΔtS/dt，于是，南向接收天線的計算數據為

TSc=TNm(t-ΔnSdt)。

(8)

西向接收天線與東向接收天線的時間差為

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

(9)

時間差所對應的采樣點數:ΔnW=ΔtW/dt，于是，西向接收天線的計算數據為

TWc=TNm(t-ΔnWdt)。

(10)

4)最后計算殘差：

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

(11)

最小殘差所對應的角度即為方位角。

殘差法算法進行方位識別的主要原理是利用4個接收天線所接收信號的相對位置計算方位角。一旦信噪比過低，導致4個接收天線接收到信號的相對位置發生較大的變化，必然影響方位識別結果。如果信號中存在少量噪聲，對于4個接收天線接收到信號的相對位置并無太大影響，那么殘差法算法依然可以得到非常好的計算結果。

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

2 合成數據算例

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

表1 合成數據采用的材料性質

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

表2 合成數據采用的參數

圖5 用于鉆孔雷達合成數據的模型Fig.5 Model for the borehole radar synthetic data

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

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

1)設計一個三維數組，讓道數、方位角數和采樣點數為三維數組的3個維數。

2)根據各道的方位角數據和接收到的電場數據對三維數組進行賦值。

4)對于該三維數組中的數據進行顯示，得到該模型的三維成像結果圖。

圖6 各天線雷達信號波形Fig.6 Signals for four antennas

圖7 深度為11.75 m的方位識別結果圖Fig.7 Azimuth angle trace at the depth of 11.75 m

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

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

圖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

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

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

3 結論

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

2)本文實現了對殘差法的算法改進，成功實現了殘差法對于多目標的方位識別，并在此基礎上提出并實現了鉆孔雷達縱橫切片以及三維成像方法。

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

[1] 曾昭發，李文奔，習建軍，等. 基于DOA估計的陣列式探地雷達逆向投影目標成像方法[J].吉林大學學報(地球科學版)，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] 朱自強，彭凌星，魯光銀，等. 基于互相關函數對鉆孔雷達層析成像的改進[J].吉林大學學報(地球科學版)，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] 張理輕,馬曄,楊宇,等.鉆孔雷達數據處理技術及分析[J].地震工程學報,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] 陳建勝,陳從新.鉆孔雷達技術的發展和現狀[J].地球物理學進展,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] 曾昭發,劉四新,王者江,等.探地雷達方法原理及應用[M].北京：科學出版社,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.