曲昕馨，李桐林，王 飛

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

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

曲昕馨，李桐林，王 飛

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

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

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

0 引言

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

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

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

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

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

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

式中：ERj是從第j個(gè)采樣點(diǎn)開始的能量比的值；Ai是第i個(gè)采樣點(diǎn)的信號振幅。

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

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

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

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

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

2）將每個(gè)振幅分配一個(gè)索引顏色。

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

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

5）在能量比圖像中選取走時(shí)所對應(yīng)的參考顏色r。

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

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

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

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

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

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

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

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

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

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

3 合成數(shù)據(jù)實(shí)例

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

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

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

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

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

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

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

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

4 實(shí)測數(shù)據(jù)實(shí)例

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

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

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

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

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

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

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

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

5 結(jié)論

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

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

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

（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］劉四新，周俊峰，吳俊軍，等.金屬礦鉆孔雷達(dá)探測的數(shù)值模擬［J］.吉林大學(xué)學(xué)報(bào)：地球科學(xué)版 ，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］冉利民，劉四新，李玉喜，等.影響跨孔雷達(dá)層析成像效果的幾個(gè)因素［J］.吉林大學(xué)學(xué)報(bào)：地球科學(xué)版，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ù)字圖像分割法的跨孔雷達(dá)走時(shí)層析成像.吉林大學(xué)學(xué)報(bào)：地球科學(xué)版，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

國家重大科學(xué)儀器設(shè)備開發(fā)專項(xiàng)（2011YQ05006009）；中國地質(zhì)調(diào)查局吉黑東部綜合找礦方法研究項(xiàng)目（12120113098400）

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

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