基于數(shù)字圖像分割法的跨孔雷達走時層析成像

曲昕馨，李桐林，王 飛

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

基于數(shù)字圖像分割法的跨孔雷達走時層析成像

曲昕馨，李桐林，王 飛

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

跨孔雷達走時層析成像主要利用雷達波的走時進行反演，走時提取的正確與否將直接影響到層析成像的效果。數(shù)字圖像分割法基于凸集投影（POCS）方法，使用能量比彩色圖像分割技術準確提取走時。數(shù)字圖像分割法提取走時首先應用在折射地震波的數(shù)據處理中。筆者首次將數(shù)字圖像分割法提取走時的方法應用到跨孔雷達走時層析成像中，使用迭代線性反演算法重建了雷達波速度場。反演過程中，使用最小二乘QR分解法（LSQR）求解線性方程組，利用彎曲射線追蹤技術構建雅可比矩陣，走時的計算值則由多模板快速推進算法（MSFM）得到。為了驗證數(shù)字圖像分割法在走時層析成像中的效果，使用一組合成數(shù)據及一組實測數(shù)據分別對比了基于數(shù)字圖像分割法和基于傳統(tǒng)能量比法獲得的層析成像反演結果。對比結果表明，使用數(shù)字圖像分割法得到的層析成像結果更為精確，誤差更小，能為判斷地下雷達波速度場提供更為有力的幫助。

跨孔雷達；走時提??；走時層析成像；數(shù)字圖像分割法；凸集投影

0 引言

跨孔雷達走時層析成像是近年來非常流行的一種雷達數(shù)據反演方法，被廣泛地應用在水文［1］、環(huán)境［2］、巖土工程［3］和金屬礦探測［4］等領域。影響跨孔雷達走時層析成像的因素有很多，包括射線覆蓋角度、正則化因子、測量裝置排布以及走時提取等［5］?？缈桌走_走時層析成像主要是利用雷達波的走時進行反演，因此，走時提取的正確與否將直接影響到層析成像的效果。

直達波走時的提取首先是用于地震記錄的編輯與處理。傳統(tǒng)的走時拾取方法主要分為兩大類：一類是基于地震記錄瞬時特征的方法，其中最具代表性的是Coppens［6］于1985年提出的能量比法，這類方法對噪聲比較敏感，當?shù)卣鹩涗浀脑肼曒^嚴重時，走時拾取的精度會受到影響；另一類方法是基于地震記錄整體特征的方法，其中最常見的是相關法［7］，這類算法雖然對噪聲有較好的抑制作用，但受到地震道之間相似性等因素的影響，對于復雜地震記錄，難以準確拾取走時。因此，研究者們又相繼提出了一些改進算法。如Al-Ghamdi［8］通過使用自適應閾值改進了Coppens的方法，Sabbione等［9］使用熵和分維方法改進了Coppens的方法等。

在鉆孔雷達走時層析成像中，走時提取技術的研究發(fā)展得比較晚。Messinger［10］將互相關方法與STA/LTA（信號短時平均值/信號長時平均值）技術相結合來自動拾取走時；Irving等［11］將互相關方法應用到共射線角度集記錄的直達波走時的提取中；Tronicke［12］提出了一種基于統(tǒng)計標準的走時提取方法；Bernard等［13］提出了一種使用赤池信息準則和連續(xù)小波變換進行走時提取的方法。

筆者將一種新的走時提取方法——數(shù)字圖像分割法應用到跨孔雷達走時層析成像中。數(shù)字圖像分割法由Mousa提出，是Coppens方法的一種改進算法，最初應用在折射地震波數(shù)據的走時提取中。該方法使用基于凸集投影（POCS）［14］的能量比彩色圖像分割技術準確提取走時［15］。為了驗證該方法在走時層析成像中的效果，筆者使用一組合成數(shù)據及一組實測數(shù)據分別對比基于數(shù)字圖像分割法和基于傳統(tǒng)能量比法獲得的層析成像反演結果。在反演算法上，選取常用的迭代線性反演算法［16］，并使用最小二乘QR分解法（LSQR）［17］求解線性方程組。反演過程中雅可比矩陣的構造由彎曲射線追蹤技術［18］實現(xiàn)，走時的計算值由多模板快速推進算法（MSFM）［19］得到。

1 數(shù)字圖像分割法提取走時

數(shù)字圖像分割法的分割是在能量比圖像上進行的。由于走時的能量比通常比后續(xù)波的能量比高，可以在索引顏色圖上為它們分配不同的顏色。一旦獲得包含走時的分割后的圖像（走時被分配一個參考顏色），可以很容易地將可能走時位置的索引從能量比圖像上映射到它們對應的雷達記錄上。按照Coppens［6］的方法，計算能量比時，首先需要分別計算一個滑動窗口以及一個累積窗口內的能量，然后將滑動窗口內的能量與累積窗口內的能量相比?；瑒哟翱诘拈L度為L，累積窗口的長度為從初始采樣點開始一直到當前滑動窗口位置。能量比的公式如下：

式中：ERj是從第j個采樣點開始的能量比的值；Ai是第i個采樣點的信號振幅。

數(shù)字圖像分割法使用POCS技術進行圖像分割。POCS技術的基本思路是：給定一個參考顏色向量r，希望尋找到所有的顏色向量h∈I?H（I是給定希爾伯特空間H的一個子集的圖像），受限于一定的限制條件Ci∈H。對于m個已知的屬性，則有m個限制條件。h可以由聯(lián)立的POCS公式迭代求解獲得：

可以施加如下定義的限制條件：

式中：ε是h和r之間的最大距離；ρ是一個相關系數(shù)，滿足0≤｜ρ｜≤1；C1和C2在H中分別表示一個球體和一個圓錐體，因此它們都是凸集。在數(shù)字圖像分割法中，所需要的解的集合為這2個凸集的交集。基于公式（3）和（4）的拉格朗日公式，使用緊鄰原則，可以推導出公式（3）和（4）對應的相關投影算子公式：

下面給出基于POCS技術的數(shù)字圖像分割法的具體實現(xiàn)步驟：

1）對于一個給定的雷達數(shù)據集，首先使用公式（1）計算其能量比振幅。

2）將每個振幅分配一個索引顏色。

3）將這些能量比索引顏色轉換到一個3D像素顏色域（例如RGB域），可以得到一個彩色的能量比數(shù)字圖像。對于低信噪比的數(shù)據，還可以對能量比圖像使用中值濾波進行預處理以增強圖像質量。

4）選擇w1、w2、ε和ρ的值。

5）在能量比圖像中選取走時所對應的參考顏色r。

6）使用公式（2）、（5）和（6）計算滿足公式（3）和（4）的向量h。

7）如果hk＋1和hk之間的誤差小于一個給定的閾值，則停止；否則，重復步驟4）至6）。

8）分割后走時區(qū)域的下緣的索引即為走時的位置。

2 基于彎曲射線追蹤的走時層析成像算法

筆者使用Tikhonov正則化函數(shù)［20］作為走時層析成像反演的總目標函數(shù)：

式中：tmes和tc（s）分別為走時的觀測值（即提取的走時）和計算值；s和s0分別為待求的慢度向量和初始模型的慢度向量；Wd和Wm分別為數(shù)據加權矩陣和模型加權矩陣［21］；λ是正則化參數(shù)。最小化φ（s），可以推導出：

式中，J＝?t∕?s是雅可比矩陣，它由每個慢度元胞內的走時對慢度的偏導數(shù)組成。走時函數(shù)tc（s）取決于未知的慢度s，所以tc（s）是非線性的。為了獲得tc（s）的近似線性形式，對tc（s）進行泰勒級數(shù)展開，并且只保留一階項，tc（s）可以近似表示成tc（s）≈tc（s0）＋Js0（s－s0）。此時，式（8）變成：

為了提高反演的精度，使用迭代線性反演公式將式（9）寫成如下的迭代形式：

求解式（10）需要進行矩陣與矩陣的相乘，非常耗時間。一個更好的算法是尋找下式的最小二乘解［22］：

式（11）的解比（10）的解精確得多，并且具有更高的計算效率。筆者使用LSQR算法［17］求解式（11）。走時的觀測值tmes通過數(shù)字圖像分割法進行走時提取得到，走時的計算值tc（sk）則使用高精度的MSFM 算法［19］計算求得。彎曲射線追蹤算法［18］則被用來精確地追蹤射線路徑和構建雅可比矩陣Jk。

3 合成數(shù)據實例

為了驗證數(shù)字圖像分割法在跨孔雷達走時層析成像中的反演效果，首先使用數(shù)字圖像分割法對合成數(shù)據模型的數(shù)值模擬結果進行走時提取，并使用提取的走時進行層析成像反演。如圖1所示，合成數(shù)據模型為高11m、寬6m的隨機介質模型，采用二維時域有限差分（FDTD）算法對該模型進行數(shù)值模擬。發(fā)射天線和接收天線分別布設在距模型左邊界0.5m和右邊界0.5m的2個鉆孔中。發(fā)射天線從0.5m深度處開始移動，每次移動0.25m，共移動41次；對于每個發(fā)射天線位置，接收天線都從0.5 m深度處開始布設，每隔0.25m布設一個接收天線，每個發(fā)射天線對應41個接收天線。這樣，整個模型區(qū)域共有1 681條射線。

圖1 用來生成合成數(shù)據的隨機介質模型Fig.1 The random medium model used to generate synthetic data set

圖2 合成數(shù)據提取的走時Fig.2 The extracted traveltimes of synthetic data set

圖2為由數(shù)值模擬得到的合成數(shù)據的某幾道分別使用能量比法和數(shù)字圖像分割法提取的各條射線的走時對比。從圖2中可以看出，2條曲線的整體趨勢比較相似，但部分位置還是存在一定差異，尤其在各發(fā)射天線位置附近差異最大。這說明這2種方法提取的走時會使得兩者的反演結果在發(fā)射天線位置（即左側鉆孔處）產生較大的不同。至于這2種提取方法哪種更準確，由于隨機介質模型的理論走時無法獲得，很難進行對比。但是，可以將這2種方法提取的走時應用到走時層析成像反演中，通過對比反演結果的誤差而從一個側面來判定這2種方法提取走時的精確性。

層析成像反演時，初始模型選取為一個速度為0.1m/ns的均勻介質模型，反演網格距為0.05m，Wd和Wm分別選取單位矩陣和拉普拉斯算子。正則化參數(shù)λ的值選為25。為了測量反演后的速度模型接收天線處走時的計算值tc（sk）與提取的合成數(shù)據模擬結果的走時tmes之間的誤差，筆者使用了均方根誤差（RMS）：

式中，n為總射線數(shù)。

圖3a為基于能量比法提取的走時進行的層析成像反演結果。可以看出，反演的結果較好，但是在靠近層析成像圖的左側鉆孔的位置，反演出了很多原模型中并不存在的虛假異常。如在0.5、5、6和10.5m深度處均出現(xiàn)了比較明顯的虛假異常，這與圖2中2條走時曲線的最大差異處吻合得很好。圖3b為基于數(shù)字圖像分割法提取的走時進行的層析成像反演結果。圖3b與圖3a大體上很相似，但是圖3a中出現(xiàn)的虛假異常在圖3b中均沒有出現(xiàn)，并且圖3b的均方根誤差（0.06ns）明顯小于圖3a（0.13ns）。這說明，相比于能量比方法，數(shù)字圖像分割法提取的走時更加準確，并且基于數(shù)字圖像分割法所得到的層析成像反演結果精度更高，更加接近真實模型。

圖3 合成數(shù)據模型反演結果Fig.3 The reconstructed results of synthetic data set

4 實測數(shù)據實例

使用一組實測數(shù)據來評價基于數(shù)字圖像分割法的跨孔雷達走時層析成像方法的反演效果。該實測數(shù)據采自貴州。貴州位于我國西南地區(qū)，地下廣泛分布著可溶性的碳酸鹽巖（灰?guī)r、白云巖等），介電常數(shù)約為9（電磁波速度約為0.1m/ns）。由于該地區(qū)巖溶作用明顯，其地下多分布裂隙和溶洞。部分裂隙和溶洞中由于水的存在，其介電常數(shù)明顯增大，平均介電常數(shù)為14～18（電磁波速度為0.07～0.08 m/ns）；部分裂隙和溶洞中充滿空氣，使得其介電常數(shù)相對周圍圍巖（碳酸鹽巖）明顯減小（電磁波速度相應增大）。因此，可以通過觀察跨孔雷達走時層析成像的結果來判斷該地區(qū)地下溶洞和裂隙的分布情況及性質（含水或是含空氣）。

在數(shù)據采集過程中，2個鉆孔之間的水平距離為18m，垂直方向上的測量區(qū)域為5.6～48.0m。發(fā)射天線在左側鉆孔9～48m深度之間移動，每次移動1m。接收天線在右側鉆孔5.6～47.2m深度范圍內，每隔0.2m布設一個。反演中網格距為0.2m，其他參數(shù)與合成數(shù)據實例中的參數(shù)一致。

圖4為分別使用能量比法和數(shù)字圖像分割法對實測數(shù)據提取的發(fā)射天線位于某些深度處的走時。可以看出，兩者的差異還是比較大的。從整體上來看，各深度能量比法提取的走時幾乎都大于數(shù)字圖像分割法提取的走時。并且圖4b中還能觀察到能量比法提取的走時在某些位置上還存在著較大的突變；這是因為實測數(shù)據的某些道信噪比不高，而能量比法在處理低信噪比數(shù)據時會產生很大誤差。但是這些突變在數(shù)字圖像分割法提取的走時曲線中并沒有出現(xiàn)，說明數(shù)字圖像分割法在數(shù)據信噪比較低的情況下仍能獲得較高的走時提取精度。

圖5a和圖5b分別為基于能量比法和數(shù)字圖像分割法所提取的走時進行層析成像的反演結果。

圖4 實測數(shù)據走時Fig.4 The extracted traveltimes of field data set

圖5 實測數(shù)據模型反演結果Fig.5 The reconstructed results of field data set

可以看出：圖中最主要的異常為25～30m深度處的一個水平低速層（圖中用藍色表示），該低速層的電磁波速度為0.07～0.08m/ns，與含水溶洞對應得很好，因此可判斷該低速層為一含水溶洞；此外，圖5a和圖5b的右上角和左下角均分布零星的低速體（藍色區(qū)域），可以推斷出這2個位置也分布著一些小型的含水裂隙；而在兩圖的左上角和右下角，則分布著一些條帶狀的高速體，結合當?shù)氐牡刭|背景，可以推斷其為一些含空氣的裂隙。圖5中在11～40m處的那些較小的異常為虛假異常，是由于發(fā)射天線移動步長較大而引起的射線覆蓋較稀疏造成的。

圖5a和圖5b的區(qū)別主要體現(xiàn)在2個地方：一是圖5a在5～15m范圍內有一條傾斜的藍色低速帶，而圖5b中則沒有；二是圖5a的右下角在2條高速帶（紅色）之間出現(xiàn)了一條低速帶（藍色），而圖5b中的相應位置同樣沒有?？梢哉J為這些異常為虛假異常，推測是圖4b中所示的那些走時曲線突變造成的，這已被后期的鉆探結果證實。此外，通過對比圖5a和圖5b可以看出，圖5b對背景場的重建更加光滑，虛假異常更少一些。在反演誤差上，基于數(shù)字圖像分割法反演結果的誤差（0.95ns）要遠遠小于基于能量比法反演結果的誤差（3.62ns）。這說明，對于實測數(shù)據，基于數(shù)字圖像分割法所提取的走時的反演結果也好于基于能量比法所提取的走時的重建結果。

5 結論

1）合成數(shù)據模擬結果表明，基于數(shù)字圖像分割法提取的走時精度高于基于能量比法所提取的走時精度，數(shù)字圖像分割法具有更高的抑制干擾的能力。

2）通過對模型數(shù)據與實測數(shù)據反演結果的分析，可以發(fā)現(xiàn)基于數(shù)字圖像分割法能獲得遠遠好于基于能量比法的層析成像反演結果。它能夠減少發(fā)射天線附近以及信噪比較低情況下所產生的虛假異常，從而獲得更好的速度場重建結果。

總之，無論對于合成數(shù)據還是實測數(shù)據，使用數(shù)字圖像分割法得到的層析成像結果更為精確，誤差更小，能為判斷地下雷達波速度場提供更為有力的幫助。

（References）：

［1］Tronicke J，Holliger K，Barrash W，et al.Multivari-ate Analysis of Cross-Hole Georadar Velocity and Attenuation Tomograms for Aquifer Zonation［J］.Water Resources Research，2004，40（1）：W01519.

［2］Gloaguen E，Marcotte D，Giroux B，et al.Aubertin，Stochastic Borehole Radar Velocity and Attenuation Tomographies Using Cokriging and Cosimulation［J］.Journal of Applied Geophysics，2007，62（2）：141-157.

［3］Johnson T C，Routh P S，Barrash W，et al.A Field Comparison of Fresnel Zone and Ray-Based GPR Attenuation-Difference Tomography for Time-Lapse Imaging of Electrically Anomalous Tracer or Contaminant Plumes［J］.Geophysics，2007，72（2）：21-29.

［4］劉四新，周俊峰，吳俊軍，等.金屬礦鉆孔雷達探測的數(shù)值模擬［J］.吉林大學學報：地球科學版 ，2010，40（6）：1479-1484.

Liu Sixin，Zhou Junfeng，Wu Junjun，et al.Numerical Simulations of Borehole Radar for Metal Ore Detection［J］.Journal of Jilin University：Earth Science Edition，2010，40（6）：1479-1484.

［5］冉利民，劉四新，李玉喜，等.影響跨孔雷達層析成像效果的幾個因素［J］.吉林大學學報：地球科學版，2013，43（5）：1672-1680.

Ran Limin，Liu Sixin，Li Yuxi，et al.Several Factors Affecting Cross-Hole Radar Tomography［J］.Journal of Jilin University：Earth Science Edition，2013，43（5）：1672-1680.

［6］Coppens F.First Arrival Picking on Common-Offset Trace Collections for Automatic Estimation of Static Corrections［J］.Geophysical Prospecting，1985，33（8）：1212-1231.

［7］Molyneux J B，Schmitt D R.First-Break Timing：Arrival Onset Times by Direct Correlation［J］.Geophysics，1999，64（5）：1492-1501.

［8］Al-Ghamdi A.Automatic First Arrival Picking Using Energy Ratios［D］.Saudi Arabia：King Fahd University of Petroleum ＆ Minerals，2007.

［9］Sabbione J I， Velis D. Automatic First-Breaks Picking：New Strategies and Algorithms［J］.Geophysics，2010，75（4）：67-76.

［10］Messinger J.Effective Automatic Picking of Traveltime Data with High Precision［C］//Slob E，Yarovoy A，Rhebergen J.Proceedings of the 10th International Conference on GPR.Delft：Delft University of Technology，2004：91-94.

［11］Irving J D，Knoll M D，Knight R J.Improving Crosshole Radar Velocity Tomograms：A New Ap－proach to Incorporating High-Angle Traveltime Data［J］.Geophysics，2007，72（4）：31-41.

［12］Tronicke J.The Influence of High Frequency Uncorrelated Noise on First-Break Arrival Times and Crosshole Traveltime Tomography［J］.Journal of Environmental and Engineering Geophysics，2007，12（2）：173-184.

［13］Bernard G，Abderrezak B，Michal C.Assisted Travel-Time Picking of Crosshole GPR Data ［J］.Geophysics，2009，72（4）：35-48.

［14］Gubin L G，Polyak B T，Raik E V.The Method of Projections for Finding the Common Point in Convex Sets［J］. USSR Computational Mathematics and Mathematical Physics，1967，7（6）：1-24.

［15］Mousa W A，Al-Shuhail A A，Al-Lehyani A.A New Technique for First-Arrival Picking of Refracted Seismic Data Based on Digital Image Segmentation［J］.Geophysics，2011，76（5）：79-89.

［16］Pelton W H，Rijo L，Swift C M.Inversion of Two-Dimensional Resistivity and Induced-Polarization Data［J］.Geophysics，1978，43（4）：788-803.

［17］Paige C C，Saunders M A.LSQR：An Algorithm for Sparse Linear Equations and Sparse Least Squares［J］.ACM Transactions on Mathematical Software，1982，8（1）：43-71.

［18］Aldridge D F，Oldenburg D W.Two-Dimensional Tomographic Inversion with Finite-Difference Traveltimes［J］.Journal of Seismic Exploration，1993，2（3）：257-274.

［19］Hassouna M S，F(xiàn)arag A A.Multistencils Fast Marching Methods：A Highly Accurate Solution to the Eikonal Equation on Cartesian Domains［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，2007，29（9）：1-12.

［20］Tikhonov A N，Arsenin V Y.Solution of Ill-Posed Problems［M］.Washington：Winston ＆Sons，1977.

［21］Greenhalgh S A，Bing Z，Green A.Solutions，Algorithms and Inter-Relations for Local Minimization Search Geophysical Inversion［J］.Journal of Geophysics and Engineering，2006，3（2）：101-113.

［22］Lines L R，Treitel S.Tutorial：A Review of Least-Squares Inversion and Its Application to Geophysical Problems［J］.Geophysical Prospecting，1984，32（2）：159-186.

Cross-Hole Radar Travel-Time Tomography Based on Digital Image Segmentation

Qu Xinxin，Li Tonglin，Wang Fei

CollegeofGeoExplorationScienceandTechnology，JilinUniversity，Changchun130026，China

The effectiveness of cross-hole radar tomography depends mainly on the quality of the extracted first arrival-time.Digital image segmentation method，based on the projection onto convex sets（POCS）technique，is used to extract the first arrival-time by segmenting the color image of the energy ratio，and had previously been applied in refracted seismic first arrival-time extraction.We first applied digital image segmentation method into cross-hole radar travel-time tomography to reconstruct the velocity field using an iteratively linearized inversion approach.During the inversion，LSQR algorithm was employed to solve the system of linear equations，Jacobian matrix was constructed by the curved ray tracing technique，and the travel-time was calculated using Multistencils Fast Marching Method（MSFM）.We employed a synthetic data set and a field data set to test the effectiveness of the digital image segmentation method in travel-time tomography.For comparison，a traditional energy ratio method is also used for first arrival-time extraction.The result reflected that the tomography based on digital image segmentation method is more accurate with smaller residuals，and can provide more effective judgment for the underground velocity field.

cross-hole radar；first arrival-time extraction；travel-time tomography；digital image segmentation；projection onto convex sets

10.13278/j.cnki.jjuese.201404302

P631.2

A

曲昕馨，李桐林，王飛.基于數(shù)字圖像分割法的跨孔雷達走時層析成像.吉林大學學報：地球科學版，2014，44（4）：1340-1347.

10.13278/j.cnki.jjuese.201404302.

Qu Xinxin，Li Tonglin，Wang Fei.Cross-Hole Radar Travel-Time Tomography Based on Digital Image Segmentation.Journal of Jilin University：Earth Science Edition，2014，44（4）：1340-1347.doi：10.13278/j.cnki.jjuese.201404302.

2013-11-11

國家重大科學儀器設備開發(fā)專項（2011YQ05006009）；中國地質調查局吉黑東部綜合找礦方法研究項目（12120113098400）

曲昕馨（1984—，女，博士研究生，主要從事電磁法正反演研究，E-mail：xinxin2191＠sina.com

李桐林（1962—，男，教授，博士生導師，主要從事電磁法理論及其應用的研究，E-mail：lilaoshizh＠163.com。