基于Curvelet變換的隧道裂隙水GPR數(shù)據(jù)處理研究

朱自強(qiáng)，朱 賀*，魯光銀，王凡，譚 潔

（1．中南大學(xué)地球科學(xué)與信息物理學(xué)院，長沙 410083；2．中國水電顧問集團(tuán)貴陽勘測設(shè)計(jì)研究院，貴陽 550081）

基于Curvelet變換的隧道裂隙水GPR數(shù)據(jù)處理研究

朱自強(qiáng)1，朱 賀1*，魯光銀1，王凡2，譚 潔1

（1．中南大學(xué)地球科學(xué)與信息物理學(xué)院，長沙 410083；2．中國水電顧問集團(tuán)貴陽勘測設(shè)計(jì)研究院，貴陽 550081）

應(yīng)用地質(zhì)雷達(dá)探測隧道內(nèi)部裂隙水的分布狀況時(shí)，受隧道鋼筋網(wǎng)的影響以及噪聲干擾，導(dǎo)致探測信號(hào)中弱有效信號(hào)難以識(shí)別。Curvelet變換可以對(duì)二維信號(hào)從頻率、角度和空間位置實(shí)現(xiàn)有效反射波和干擾波的分離及降噪處理，應(yīng)用Curvelet變換對(duì)含鋼筋網(wǎng)干擾的裂隙水模型的探地雷達(dá)模擬數(shù)據(jù)進(jìn)行波場分離和降噪處理，在壓制直達(dá)波、去除鋼筋層反射信號(hào)及多次反射的基礎(chǔ)上，可以有效地提取出含裂隙水所引起的反射波信號(hào)。將該方法應(yīng)用到汝郴高速某隧道實(shí)測探地雷達(dá)數(shù)據(jù)中，較好地去除了鋼筋層反射信號(hào)的影響，準(zhǔn)確地判斷出了裂隙水的空間分布位置，證明了所提出方法的有效性。

Curvelet變換；探地雷達(dá)；隧道裂隙水；波場分離；降噪

0 引言

探地雷達(dá)法［1－2］是近些年來迅速發(fā)展起來的一種高分辨率、高效率的勘探方法，目前在隧道勘測中，特別是探測巖溶裂隙水方面應(yīng)用較多［3－6］。由于隧道周圍裂隙水被密集的鋼筋層覆蓋，用探地雷達(dá)對(duì)隧道裂隙水探測時(shí)會(huì)受到鋼筋層反射的干擾，因此，需要對(duì)雷達(dá)數(shù)據(jù)進(jìn)行處理，去除干擾和噪聲，從而提取有效的異常信號(hào)。劉斌等［7］將復(fù)信號(hào)技術(shù)應(yīng)用到探地雷達(dá)對(duì)地下裂隙水的探測中，對(duì)雷達(dá)數(shù)據(jù)中的異常分別從振幅、頻率和相位等各方面進(jìn)行綜合分析，提高了預(yù)報(bào)精度和準(zhǔn)確性，但沒能從根本上克服探地雷達(dá)法抗干擾能力弱的缺點(diǎn)。劉四新等［8］用鉆孔探地雷達(dá)對(duì)巖溶裂隙進(jìn)行勘探，勘探效果較好，但具有探測上方位不確定性和成本較高的缺點(diǎn)。鄒海林等［9］基于小波變換對(duì)探地雷達(dá)信號(hào)進(jìn)行處理，李才明［10］基于小波能譜分析對(duì)巖溶區(qū)探地雷達(dá)進(jìn)行目標(biāo)識(shí)別。小波變換很好地解決了時(shí)頻同步問題，而且對(duì)一維分段平穩(wěn)信號(hào)的表達(dá)具有最佳性能，但是對(duì)于描述二維圖像信息卻存在諸多缺點(diǎn)［11－13］。

Candes等人［14－15］在1999年提出了第一代Curvelet變換，并于2002年提出了實(shí)現(xiàn)更簡單、運(yùn)行速度更快、更便于直觀理解的第二代Curvelet變換算法。在此之后，他們又于2005年提出了兩種基于第二代Curvelet 變換理論的快速離散實(shí)現(xiàn)方法［16］，即USFFT（unequally-spaced fast Fourier transforms）和Wrapping兩種快速離散算法。2008年以來，該方法在地球物理資料處理尤其是地震信號(hào)處理中得到了迅速地研究和應(yīng)用［17－20］。在國內(nèi)，仝中飛等人［21］提出了迭代閾值法壓制地震信號(hào)中的隨機(jī)噪聲，取得了很好的效果。潘雪輝［22］提出了一種用Curvelet變換壓制地震資料中線性噪聲的方法，為地震信號(hào)中線性噪聲的壓制問題提供了較好地解決方法。

Curvelet變換建立在小波變換的基礎(chǔ)上，增加了一個(gè)方位參數(shù)，解決了小波變化在處理二維信號(hào)時(shí)的不足，而探地雷達(dá)檢測數(shù)據(jù)主要為描述一定深度范圍內(nèi)剖面的二維信號(hào)，這讓作者萌生思路，利用Curvelet變換對(duì)探地雷達(dá)數(shù)據(jù)進(jìn)行處理研究。在對(duì)Curvelet變換方法研究基礎(chǔ)上，將其應(yīng)用到隧道裂隙水的探地雷達(dá)檢測資料處理中。本研究首先模擬了在含有覆蓋鋼筋層和裂隙水的混凝土中的探地雷達(dá)波場，并通過Curvelet變換對(duì)模擬波場進(jìn)行去除直達(dá)波、提取水層信息、降噪等一系列數(shù)據(jù)處理流程，提取出有效的裂隙水異常信息；隨后將該方法應(yīng)用到汝郴高速某隧道實(shí)測地質(zhì)雷達(dá)數(shù)據(jù)，獲取了隧道內(nèi)裂隙水異常的空間分布，通過鉆孔注漿驗(yàn)證，證明了此方法的有效性。

1 Curvelet變換基本原理

小波變換在信號(hào)處理中的應(yīng)用得到了很大的發(fā)展，遺憾的是，由一維小波所生成的可分離小波只具有有限的方向，因此在處理圖像邊緣時(shí)的效果不是很好。Curvelet變換除了與小波變換一樣具有尺度和位移參數(shù)，還增加了一個(gè)方位參數(shù)，因此具有更好的方位識(shí)別能力，可以精確表達(dá)出圖像中邊緣的方向信息。

1．1 連續(xù)Curvelet變換

Curvelet變換是利用基函數(shù)與信號(hào)（或函數(shù)）的內(nèi)積形式來實(shí)現(xiàn)信號(hào)（或函數(shù)）的稀疏表示，表示為：

其中：φj，k，l表示Curvelet函數(shù)；j、k、l分別表示尺度、位置和方向參數(shù)。

Curvelet變換在頻域內(nèi)采用窗函數(shù)U來實(shí)現(xiàn)的。分別定義徑向窗函數(shù)W（r）、r∈（1／2，2）和角度窗函數(shù)V（t），t∈［－1，1］，二者滿足式（2）與式（3）。

對(duì)于每一個(gè)j≥j0，在頻域中定義窗函數(shù)Uj如式（4）。

綜上所述，我們可以得到Curvelet變換公式：

圖1 Curvelet變換示意圖Fig．1 Curvelet transform schematic diagram

1．2 離散Curvelet變換

在直角坐標(biāo)系下f［t1，t2］，0≤t1，t2＜n為輸入，Curvelet變換離散形式為：

1．3 離散Curvelet變換快速實(shí)現(xiàn)方法

本研究采用基于USFFT算法的快速離散Cur-velet變換方法，實(shí)現(xiàn)過程如下：

1）對(duì)于給定的直角坐標(biāo)系下二維函數(shù)f［t1，t2］，0≤t1，t2＜ω進(jìn)行2DFFT，得到其二維頻率域表達(dá)式：

2）在頻率域，對(duì)每一對(duì)（i，j），對(duì)函數(shù)＾f［n1，n2］重新采樣，得到函數(shù)采樣值：

其中：Pj＝｛（n1，n2）：n1，0≤n1≤n1，0＋L1，j，n2，0≤n2≤n2，0＋L2，j｝；L1，j表示窗函數(shù)支撐區(qū)間的長度；L2，j表示支撐區(qū)間的寬度。

1．4 Curvelet變換的非線性逼近能力

Curvelet變換之所以比小波變換更適合表示二維信號(hào)，其主要原因在于Curvelet變換比小波在二維平面具有更好的稀疏性。假設(shè)f∈L2［0，1］為一條曲線的光滑部分，分別用小波變換和Curvelet變換對(duì)其分解，再分別選擇最佳的M個(gè)系數(shù)對(duì)函數(shù)進(jìn)行重構(gòu)，得到三種變換下的的最佳逼近它們的逼近誤差分別為：

從重構(gòu)誤差式中可以看出，Curvelet變換的逼近率要優(yōu)于小波變換。Curvelet變換之所在表示二維曲線時(shí)有更好的逼近率，主要在于它在小波變換的基礎(chǔ)上增加了方位參數(shù)，使得它具有多方向性、多尺度性、波動(dòng)性和各向異性的特點(diǎn)，這也使得用更少的Curvelet系數(shù)就能表示二維信號(hào)的主要特征。

圖2 小波和Curvelet對(duì)二維曲線逼近示意圖Fig．2 Schematic of two－dimensional curve approximation with wavelet and Curvelet

2 模擬數(shù)據(jù)處理分析

針對(duì)表層存在鋼筋干擾下的含裂隙水混凝土模型進(jìn)行數(shù)值模擬，模型大小為5．2 m×1．28 m，混凝土介電常數(shù)ε0為9，電導(dǎo)率σ0為0．01 s／m。距離表層0．3 m處設(shè)置一排鋼筋，半徑0．01 m。鋼筋層下方設(shè)置一個(gè)水平和一個(gè)傾斜含裂隙水，橫截面大小均為0．8 m×0．08 m，介電常數(shù)ε1為81，電導(dǎo)率σ1為0．05 s／m，模擬天線頻率為400 MHz。模型剖面和雷達(dá)模擬剖面如圖3所示。

圖3 雷達(dá)數(shù)值模擬Fig．3 GPR numerical simulation

在圖3中，表層直達(dá)波強(qiáng)度較高，弱化了下面鋼筋層和裂隙水層的反射信號(hào)。將總波場變換到Curvelet域后，然后在角度窗函數(shù)下選擇變換域中θ＝0（θ為方向因子）附近的Curvelet系數(shù)進(jìn)行重構(gòu)，得到直達(dá)波的波場信號(hào)，進(jìn)而分離得到的直達(dá)波波場信號(hào)和剩余波場信號(hào)（圖4）。

圖4 波場分離Fig．4 Wave field separation

當(dāng)直達(dá)波去除后，下面鋼筋層和裂隙水層的反射信息加強(qiáng)，但是二次反射和干擾信息的存在干擾了異常信號(hào)，通過Curvelet變換的降噪功能可以去掉二次反射及其他干擾信號(hào)。基于Curvelet變換的特點(diǎn)，前人提出了各種關(guān)于降噪的閾值方法［23－24］，本研究采取一種自適應(yīng)閾值方法［25］，即隨著尺度的改變，每一層閾值可以相應(yīng)地改變。閾值函數(shù)表達(dá)式如下：其中 sigma為噪聲的方差估計(jì)值；N為Curvelet域系數(shù)矩陣大?。籹為Curvelet域尺度分解下的層數(shù)。

圖5 波場降噪Fig．5 Wavefield denoising

圖5為sigma系數(shù)設(shè)為700時(shí)降噪的效果圖，降噪后，鋼筋和裂隙水的異常信息已非常明顯。原始二維信號(hào)在經(jīng)過Curvelet變換后會(huì)得到的一系列Curvelet系數(shù)矩陣，將這些系數(shù)矩陣在MATLAB中生成時(shí)間－頻率圖像，通過對(duì)比原始圖像異常位置可以識(shí)別出包含有效信號(hào)和干擾信號(hào)的矩陣。如果我們只需要裂隙水的異常信息，只需對(duì)這些離散化的系數(shù)矩陣進(jìn)行人工識(shí)別，并挑選出含裂隙水信息的系數(shù)矩陣，進(jìn)行Curvelet反變換，得到只含裂隙水異常信息的二維數(shù)據(jù)。圖6為應(yīng)用Curvelet變換波場分離方法得到的裂隙水異常信息數(shù)據(jù)。

圖6 提取裂隙水波場信息Fig．6 Wave signal extraction of fissure water

處理結(jié)果顯示，用Curvelet變換對(duì)雷達(dá)模擬數(shù)據(jù)的處理效果良好，波場分離和降噪都能夠基本達(dá)到目的。兩個(gè)方向有差異的裂隙水的異常信息在處理過程中表現(xiàn)出不同的處理效果，這也說明了Curvelet變換在處理圖像信息時(shí)帶有很強(qiáng)的方向選擇性。由于水平的裂隙水在方向上與上面鋼筋層保持平行，用Curvelet變換離散化圖像信息時(shí)，水平裂隙水與鋼筋層的信息在角度窗函數(shù)選擇下沒能夠完全分離，只能夠通過尺度來分離部分信息。這樣傾斜的裂隙水異常信息才會(huì)顯得比水平裂隙更明顯。

3 實(shí)測數(shù)據(jù)處理分析

在Curvelet變換方法研究和模擬處理的基礎(chǔ)上，對(duì)汝郴高速公路某隧道內(nèi)采集所得探地雷達(dá)數(shù)據(jù)進(jìn)行處理。探測目的是找出隧道邊墻和底部存在的裂隙水，確定這些裂隙水的位置及相關(guān)信息，為工程治理工作提供較好的參考信息。探測中使用的天線為400 MHz，采樣點(diǎn)數(shù)設(shè)為512，參考設(shè)計(jì)方案，鋼筋間距約為50 cm。

圖7（a）為截取采集數(shù)據(jù)中某一異常位置信息，直達(dá)波和鋼筋產(chǎn)生的干擾信息強(qiáng)度較強(qiáng)，同時(shí)數(shù)據(jù)中含有一些噪聲干擾。首先我們要將有效異常信息和鋼筋層干擾信息及直達(dá)波信息分離開，提取含有有效異常信息的數(shù)據(jù)。圖7（b）為分離出的鋼筋層信息及一些干擾信息，圖7（c）為提取得到的含有裂隙水層異常信息的波場。由于波場分離更多在于角度窗下的選擇，圖7（b）中鋼筋層干擾和直達(dá)波與圖7（c）中有效信息被基本分離開，但仍然有一些因?yàn)榕c裂隙水異常信息方向一致而未被分離的遺留信息，且存在一些噪聲。遺留下的干擾信息強(qiáng)度較弱，可以通過Curvelet降噪功能，應(yīng)用自適應(yīng)尺度閾值法，通過尺度窗的選擇進(jìn)行降噪，取sigma系數(shù)為48 000得到降噪后的異常圖7（d）。圖7（d）異常信息清晰，進(jìn)而可以判斷裂隙水的空間位置。

根據(jù)探測和處理結(jié)果，對(duì)發(fā)現(xiàn)的異常位置進(jìn)行鉆孔驗(yàn)證如圖8所示。驗(yàn)證結(jié)果支持了該方法對(duì)于裂隙水空間位置的判斷，并在該異常位置進(jìn)行了注漿治理，總注漿量為10 T左右水泥和38 T左右化學(xué)漿。對(duì)于本次探測中發(fā)現(xiàn)的其他幾個(gè)異常位置，用文中所述方法進(jìn)行分析并打鉆注漿，在判斷得到的異常較大位置都能注入不同量的漿液，驗(yàn)證了該處理方法的準(zhǔn)確性與可行性。

4 結(jié)論

1）基于Curvelet變換，從波場分離的角度來提高探地雷達(dá)異常信號(hào)的辨識(shí)度，較好地解決了直達(dá)波和鋼筋層強(qiáng)反射對(duì)下層異常信息的干擾問題。

2）角度參數(shù)的增加使得Curvelet變換在二維信號(hào)的分離上有了很強(qiáng)的優(yōu)勢，但是在干擾信息與有效信息方向性一致的情況下，分離效果并不理想。

圖7 工程應(yīng)用實(shí)例Fig．7 Examples of application engineering

圖8 現(xiàn)場鉆孔照片F(xiàn)ig．8 Photos of drilling

3）閾值函數(shù)的選取直接關(guān)系到Curvelet降噪的效果，采取適應(yīng)不同尺度的自適應(yīng)閾值函數(shù)，基本能夠達(dá)到消除噪聲的目的，但閾值函數(shù)thresh沒有從角度去適應(yīng)系數(shù)矩陣，找到同時(shí)適應(yīng)尺度和角度的閾值函數(shù)將是下一步急需解決的問題。

［1］ 李大心．探地雷達(dá)方法與應(yīng)用［M］．北京：地質(zhì)出版社，1994．

［2］ 楊天春，呂紹林，王齊仁．探地雷達(dá)檢測道路厚度結(jié)構(gòu)的應(yīng)用現(xiàn)狀及進(jìn)展［J］．物探與化探，2003，27（1）：79－82．

［3］ 楊龍江，閻滿存．探地雷達(dá)在隱伏巖溶區(qū)工程地質(zhì)勘探中的應(yīng)用［J］．物探化探計(jì)算技術(shù)，1999，21（2）：172－176．

［4］ 葛雙成，邵長云．巖溶勘察中的探地雷達(dá)技術(shù)及應(yīng)用［J］．地球物理學(xué)進(jìn)展，2005，20（2）：476－481．

［5］ 李孟娟，李川．探地雷達(dá)檢測隧道襯砌厚度的研究［J］．物探化探計(jì)算技術(shù)，2008，30（3）：231－234．

［6］ 李仁海，楊磊，許新剛，等．地質(zhì)雷達(dá)探測技術(shù)在巖溶地形勘察中的應(yīng)用［J］．物探化探計(jì)算技術(shù)，2009，31（5）：442－446．

［7］ 劉斌，李術(shù)才，李樹忱，等．復(fù)信號(hào)分析技術(shù)在地質(zhì)雷達(dá)預(yù)報(bào)巖溶裂隙水中的應(yīng)用研究［J］．巖土力學(xué)，2009，30（7）：2191－2196．

［8］ 劉四新，曾昭發(fā)，徐波．利用鉆孔雷達(dá)探測地下含水裂縫［J］．地球物理學(xué)進(jìn)展，2006，21（2）：620－624．

［9］ 鄒海林，寧書年，林捷．小波理論在探地雷達(dá)信號(hào)處理中的應(yīng)用［J］．地球物理學(xué)進(jìn)展，2004，19（2）：268－275．

［10］ 李才明，王良書，徐鳴潔，等．基于小波能譜分析的巖溶區(qū)探地雷達(dá)目標(biāo)識(shí)別［J］．地球物理學(xué)報(bào)，2006，49（5）：1499－1504．

［11］WALKER J S，CHEN Y．Image denoising using treebased wavelet subband correlations and shrinkage［J］．Optical Engineering，2000，39（11）：2900－2908．

［12］楊鳳娟，羅省賢．小波變換在探地雷達(dá)檢測鋼筋中的應(yīng)用［J］．物探化探計(jì)算技術(shù)，2009，31（4）：354－360．

［13］黃敏，朱德兵，郭政學(xué)，等．連續(xù)小波變換在探地雷達(dá)信號(hào)分析中的應(yīng)用研究［J］．物探化探計(jì)算技術(shù)，2012，34（5）：593－598．

［14］CANDES E J，DONOHO D L．Curvelets：A surprisingly effective nonadaptive representation for objects with edges［R］．DTIC Document，2000．

［15］STARCK J，CANDèS E J，DONOHO D L．The curvelet transform for image denoising［J］．Image Processing，IEEE Transactions on．2002，11（6）：670－684．

［16］CANDES E，DEMANET L，DONOHO D，et al．Fast discrete curvelet transforms［J］．Multiscale Modeling＆Simulation，2006，5（3）：861－899．

［17］Neelamani R，Baumstein A I，Gillard D G，et al．Coherent and random noise attenuation using the curvelet transform［J］．The Leading Edge．2008，27（2）：240－248．

［18］姜宇東，楊勤勇，何柯，等．基于曲波變換的地面微地震資料去噪方法研究［J］．石油物探，2012，51（6）：620－624．

［19］吳愛弟，趙秀玲．基于曲波變換的地震信號(hào)去噪方法［J］．中國石油大學(xué)學(xué)報(bào)：自然科學(xué)版，2010，34（003）：30－33．

［20］張?jiān)茝?qiáng)，張培林，王國德，等．基于曲波變換和色度模型的彩色圖像去噪［J］．中國圖象圖形學(xué)報(bào)，2012，17（012）：1472－1477．

［21］仝中飛，王德利，劉冰．基于Curvelet變換閾值法的地震數(shù)據(jù)去噪方法［J］．吉林大學(xué)學(xué)報(bào)：地球科學(xué)版，2008：48－52．

［22］潘雪輝．曲波變換在地震信號(hào)去噪中的應(yīng)用［D］．西安：西安理工大學(xué)，2009．

［23］HERRMANN F J，HENNENFENT G．Non-parametric seismic data recovery with curvelet frames［J］．Geophysical Journal International，2008，173（1）：233－248．

［24］HERRMANN F J，MOGHADDAM P，STOLK C C．Sparsity-and continuity-promoting seismic image recovery with curvelet frames［J］．Applied and Computational Harmonic Analysis，2008，24（2）：150－173．

［25］劉靜．地震圖像噪聲壓制在Curvelet域的閾值選取與噪聲定位［D］．成都：西南交通大學(xué)，2011．

Processing of GPR data in tunnel fissure water based on Curvelet transform

ZHU Zi-qiang1，ZHU He1*，LU Guang-yin1，WANG Fan2，TAN Jie1
（1．School of Geosciences and Info-Physics，Central South University，Changsha 410083，China；2．HYDROChina Guiyang Engineering Corporation，Guiyang 550081，China）

When we detect the distribution of fissure water in tunnel with ground penetrating radar，it is hard to identify the weak effective signal because of the influence of mesh reinforcement and noise．Curvelet transform can achieve effective wave field separation and denoising of two－dimensional signal from the frequency，angle and spatial location．The Curvelet transform is applied to wave field separation and signal denoise processing by using of GPR simulated data from fissured water which is disturbed seriously by the existence of mesh reinforcement．After the direct signal suppression，removal of the reflected signal and multiple reflections from mesh reinforcement，we can effectively extract the reflected wave signal which is caused by the water bearing structure crack．Our method is applied to real data measured in a tunnel of the Ruchen highway．The results show that after removal of the reflected signal from mesh reinforcement we can accurately determine the location of fissure water，which proves the effectiveness of the proposed method．

Curvelet transform；ground-penetrating radar（GPR）；fissure water in tunnel；wavefield separation；denoising

P 631．4

A

10．3969／j．issn．1001-1749．2014．05．10

1001-1749（2014）05-0571-06

2014-03-33 改回日期：2014-06-22

國家自然科學(xué)基金項(xiàng)目（41174061）；中南大學(xué)自由探索計(jì)劃（2011QNZT011）

朱自強(qiáng)（1964-），男，教授，博士生導(dǎo)師，主要從事地質(zhì)災(zāi)害探測與監(jiān)測工作，E-mail：13507319431＠139．com。

*通訊作者：朱賀（1988-），男，碩士，主要從事地球物理信號(hào)處理研究，E-mail：zhuhe24＠sina．com。