馬博文 黃清華,?

北京大學(xué)學(xué)報(bào)(自然科學(xué)版) 第59卷 第3期 2023年5月

Acta Scientiarum Naturalium Universitatis Pekinensis, Vol. 59, No. 3 (May 2023)

10.13209/j.0479-8023.2023.015

國家自然科學(xué)基金(42021003)資助

2022–04–30;

2022–05–27

探地雷達(dá)疊前逆時(shí)偏移成像影響因素研究

馬博文1,2黃清華1,2,?

1.北京大學(xué)地球與空間科學(xué)學(xué)院地球物理學(xué)系, 北京 100871; 2.河北紅山巨厚沉積與地震災(zāi)害國家野外科學(xué)觀測站, 邢臺(tái) 054000; ?通信作者, E-mail: huangq@pku.edu.cn

為了分析探地雷達(dá)數(shù)據(jù)處理過程中各種可能的誤差對(duì)逆時(shí)偏移算法成像效果的影響, 通過數(shù)值模擬測試, 對(duì)模型存在速度誤差和數(shù)據(jù)存在噪聲的情況進(jìn)行分析。對(duì)于模型存在速度誤差的情況, 選用不同模型速度、不同異常體深度以及簡化速度模型的情況進(jìn)行成像效果測試。結(jié)果表明, 隨著速度誤差加大, 異常體位置偏離也變大, 且異常體深度越大, 位置偏離越大。對(duì)于各層速度相差不大的層狀結(jié)構(gòu), 使用簡化的均勻半空間模型進(jìn)行成像, 依然可以較好地分辨地下的異常體。對(duì)于存在噪聲的情況, 選用加入不同噪聲水平的數(shù)據(jù)進(jìn)行成像測試, 結(jié)果表明逆時(shí)偏移(RTM)算法對(duì)噪聲有較好的壓制效果, 即使接收的數(shù)據(jù)中有效信號(hào)的信噪比在?15dB 左右, 該算法仍然能將對(duì)應(yīng)的結(jié)構(gòu)成像出來。

探地雷達(dá); 三維逆時(shí)偏移; 模型速度誤差; 數(shù)據(jù)噪聲

探地雷達(dá)(ground penetrating radar, GPR)使用天線向地下發(fā)出高頻電磁波(0.01～1GHz), 通過接受電性異常體的反射波來推斷地下結(jié)構(gòu), 是探測地下淺層結(jié)構(gòu)的重要手段[1], 具有簡單快捷、無損探測和分辨率高等優(yōu)勢[2], 應(yīng)用領(lǐng)域非常廣泛[3], 已在地質(zhì)勘探[4]、沉積學(xué)研究[5]、考古[6]、地質(zhì)災(zāi)害評(píng)估[7]、城市工程[8]以及地下水監(jiān)測等多個(gè)領(lǐng)域發(fā)揮重要作用。

GPR 采集的原始數(shù)據(jù)是時(shí)間記錄, 反映出來的是一個(gè)假的未聚焦圖像[9], 相對(duì)于真實(shí)情況, 異常體的結(jié)構(gòu)和位置可能存在較大的偏差。偏移處理可以將原始圖像從時(shí)間域轉(zhuǎn)換到深度域, 并將圖像中的異常體還原為原本的形狀和位置。考慮到 GPR使用的電磁波工作頻率高, 與彈性波遵循相同形式的波動(dòng)方程, 因此勘探地震中的許多偏移處理方法也被引入 GPR 的數(shù)據(jù)處理過程[10], 逆時(shí)偏移(rever-se time migration, RTM)便是其中一種處理方法。多數(shù)偏移方法對(duì)波場的模擬計(jì)算, 求解的是波動(dòng)方程的近似形式——單程波動(dòng)方程, RTM 則是通過直接求解原始波動(dòng)方程來模擬波場傳播, 其原理比較直觀, 并且能夠應(yīng)用于處理具有大傾角反射體或多次波的復(fù)雜結(jié)構(gòu), 與其他偏移方法相比, 具有很高的精度[11?12]。

Baysal 等[13]和 McMechan[14]分別提出彈性波的疊前和疊后 RTM 成像算法。之后, Fisher 等[9]將二維 RTM 成像算法引入 GPR 數(shù)據(jù)處理中, 并指出二維 RTM 在大角度的反射體、產(chǎn)生衍射的結(jié)構(gòu)界面和空間不均勻結(jié)構(gòu)等情況下可以帶來顯著的改善。朱尉強(qiáng)[15]發(fā)展了 GPR 的并行三維 RTM 算法, 張崇民[16]利用 GPR 三維 RTM 成像方法研究工程隱蔽病害, 焦翠翠[17]將 RTM 應(yīng)用在礦井富水區(qū)的 GPR 三維成像研究中。

對(duì)于實(shí)際數(shù)據(jù), RTM 成像質(zhì)量可能受數(shù)據(jù)質(zhì)量和波場模擬電性模型的影響。楊仁虎等[18?20]對(duì)地震波 RTM 算法中震源子波形態(tài)、直達(dá)波、速度結(jié)構(gòu)和噪聲等因素對(duì)成像精度的影響進(jìn)行過一系列研究, 對(duì) RTM 在地震成像實(shí)踐中的應(yīng)用起到一定的參考作用。本文針對(duì) GPR 三維 RTM 成像算法進(jìn)行研究, 探討各個(gè)影響因素對(duì) GPR 的 RTM 成像效果的影響。

1 探地雷達(dá)逆時(shí)偏移成像算法

RTM 成像方法通過直接求解雙向波動(dòng)方程實(shí)現(xiàn)波場外推, 比其他偏移算法精度更高。圖 1 為疊前 RTM 算法流程示意圖。對(duì)于共發(fā)射點(diǎn)的 GPR 數(shù)據(jù), 首先根據(jù)實(shí)際情況建立電性結(jié)構(gòu)模型, 在發(fā)射點(diǎn)位置激發(fā)實(shí)際發(fā)射波形, 計(jì)算介質(zhì)中的正傳波場(,,,); 然后在接收點(diǎn)將接收波形進(jìn)行反向時(shí)間軸上的激發(fā), 計(jì)算介質(zhì)中的反傳波場(,,,); 最后對(duì)模擬得到的正傳波場和反傳波場應(yīng)用互相關(guān)成像條件, 得到成像結(jié)果(,,,)[21]:

其中, T 為模擬中波場傳播的總時(shí)間, τ 為波場傳播的某一時(shí)刻。

對(duì)于一個(gè)區(qū)域的 RTM 成像, 需要將發(fā)射源布置在不同的位置來得到反射數(shù)據(jù), 然后對(duì)每個(gè)位置的共發(fā)射點(diǎn)數(shù)據(jù)分別進(jìn)行疊前 RTM 成像, 最后將成像結(jié)果疊加, 得到最終的地下結(jié)構(gòu)信息。

在對(duì)波場和的模擬中, 采用空間 4 階時(shí)間 2階的時(shí)域有限差分法[22]求解波動(dòng)方程。算法的穩(wěn)定性條件[23]要求庫朗數(shù)小于 6/7, 即波在一個(gè)時(shí)間步長內(nèi)走過的距離小于約 0.85 個(gè)空間步長; 算法的頻散關(guān)系要求一個(gè)波長包含至少 10 個(gè)空間網(wǎng)格[24]。

2 逆時(shí)偏移成像效果影響因素分析

模擬數(shù)據(jù)與實(shí)際數(shù)據(jù)存在很大的差別。數(shù)值模擬試驗(yàn)中, 速度結(jié)構(gòu)為準(zhǔn)確已知, 通過共中心點(diǎn)剖面(common midpoint profile, CMP)或鉆孔等方法得到的地下速度模型, 可能與實(shí)際情況存在一定程度的偏差, 不能準(zhǔn)確地反映地下實(shí)際速度。另外, 實(shí)際數(shù)據(jù)中可能出現(xiàn)背景噪音和雜波等不需要的信號(hào), 對(duì)成像造成干擾。針對(duì)這些情況, 我們通過數(shù)值模擬實(shí)驗(yàn), 討論速度模型差異以及噪聲對(duì) RTM成像效果的影響。

2.1 模型設(shè)置及理想狀況下的成像結(jié)果

數(shù)值模擬采用的模型設(shè)置如下: 模型長度為 5m, 寬度為 2m, 深度為 4m, 背景相對(duì)介電常數(shù)為5, 相對(duì)磁導(dǎo)率為 1, 電導(dǎo)率為零。如圖 2 所示, 在軸中間深度為 2m 處, 沿軸依次放置 3 個(gè)半徑為0.1m 的點(diǎn)狀空腔。發(fā)射源采用主頻為 400MHz 的雷克子波。進(jìn)行數(shù)值計(jì)算時(shí), 空間網(wǎng)格間距取為0.02m×0.025m×0.02m, 模型區(qū)域共 250×80×200 個(gè)網(wǎng)格點(diǎn), 一個(gè)波長范圍內(nèi)包含大約 16 個(gè)網(wǎng)格, 避免空間頻散。時(shí)間步長為 0.0158ns, 庫朗數(shù)為 0.4, 滿足穩(wěn)定性條件。

圖 3 給出理想狀態(tài)下 RTM 成像結(jié)果方向的豎直切片。經(jīng)過逆時(shí)偏移成像, 3 個(gè)異常體的位置都得到很好的還原。

2.2 模型速度誤差

對(duì)于實(shí)際數(shù)據(jù), 需要通過各種方法(如 CMP 和鉆孔等)獲得介質(zhì)速度, 建立 RTM 初始模型。但是, 獲得的速度可能存在一定的偏差, 不能準(zhǔn)確地代表實(shí)際的地下介質(zhì)速度結(jié)構(gòu)。本文通過模擬實(shí)驗(yàn), 測試模型速度結(jié)構(gòu)差異對(duì)成像結(jié)果的影響。

大圓點(diǎn)示意 3 個(gè)點(diǎn)狀空腔的尺度和位置; 三角形示意發(fā)射源位置, 其坐標(biāo)為(2.5m, 1m, 0); 點(diǎn)線示意接收點(diǎn)的 7 條測線, 從y=0.25m到 y=1.75m, 每隔 0.25m 布置一條測線, 每條測線上的測點(diǎn)(點(diǎn)線上的小圓點(diǎn))間距為 0.1 m

設(shè)置 RTM 初始模型的介質(zhì)速度分別比實(shí)際情況快或慢 2%, 5%和 10%, 分別在這些初始模型上進(jìn)行 RTM 成像。取各個(gè)結(jié)果的方向豎直切片進(jìn)行對(duì)比觀察, 結(jié)果如圖 4 所示。當(dāng)介質(zhì)速度誤差為 2%時(shí), 異常體的形態(tài)與理想狀態(tài)差別不大, 介質(zhì)速度的誤差幾乎沒有影響(圖 4(a))。當(dāng)介質(zhì)速度誤差達(dá)到 5%時(shí), 對(duì)異常體的聚焦效果減弱, 導(dǎo)致異常體的成像形態(tài)有所展寬, 成像位置也有所偏差(圖 4(b))。當(dāng)介質(zhì)速度誤差達(dá)到 10%時(shí), 異常體的波形已無法聚焦到一點(diǎn), 呈現(xiàn)類似球形波陣面的形狀(圖 4(c)和(d)); 成像位置的偏離也非常大, 介質(zhì)速度比實(shí)際情況小的時(shí)候成像位置偏上(圖 4(c)), 介質(zhì)速度比實(shí)際情況大的時(shí)候成像位置偏下(圖 4(d))。

在介質(zhì)速度誤差相同的情況下, 如果異常體的深度不同, 對(duì)成像結(jié)果的影響也不同。圖 5 給出介質(zhì)速度誤差為 5%的情況下, 分別在異常體深度為 1m 和 3m 的成像結(jié)果。結(jié)合圖 4(b)可以看到, 對(duì)于相同的速度誤差, 異常體深度越大, 成像結(jié)果的誤差越大。

對(duì)于水平層狀結(jié)構(gòu), 使用均一介質(zhì)模型測試RTM 的成像效果。如圖 6(a)所示, 在原始模型的基礎(chǔ)上, 分別在深度 0.8～1.2m 和 1.2～1.6m 之間加入速度與背景相差 10%的高速層和低速層, 對(duì)應(yīng)的相對(duì)介電常數(shù)為 4.1 和 6.2, RTM 成像所用的介質(zhì)模型不變, 仍然為均一速度模型。從成像結(jié)果(圖 6(b)) 可以看到, 界面和異常體的位置均能成像出來。雖然由于界面反射信號(hào)太強(qiáng)而使得異常體的信號(hào)較弱, 但是其位置沒有很大的偏差, 說明對(duì)于水平層狀結(jié)構(gòu), 可以取各層速度的大致平均值作為均一速度結(jié)構(gòu)進(jìn)行 RTM 成像, 也能較好地恢復(fù)地下異常體的位置。

紅色虛線表示異常體的實(shí)際位置。成像結(jié)果為波場互相關(guān)的值, 表明地下結(jié)構(gòu)信息, 該數(shù)值本身沒有實(shí)際意義, 故未給出色標(biāo), 下同

2.3 噪聲穩(wěn)定性測試

對(duì)于實(shí)際的 GPR 數(shù)據(jù), 噪聲往往是一個(gè)不可忽略的問題, 各種成像算法也需要對(duì)數(shù)據(jù)中的噪聲有一定的魯棒性。本研究的數(shù)值模擬測試中, 對(duì)接收數(shù)據(jù)分別加入相對(duì)于有效信號(hào)具有不同噪聲水平的高斯白噪聲, 然后進(jìn)行 RTM 成像, 檢測噪聲干擾對(duì)RTM 成像效果的影響。

圖 7(a1)～(a3)展示部分加入不同噪聲之后的接收數(shù)據(jù)。噪聲水平與有效信號(hào)相當(dāng)(信噪比為 0dB) 時(shí), 有效信號(hào)較為清晰(圖 7(a1)); 信噪比為?15dB時(shí), 有效信號(hào)勉強(qiáng)可見(圖 7(a2)); 信噪比為?30dB時(shí), 有效信號(hào)肉眼不可見(圖 7(a3))。

(a)速度比實(shí)際情況小 2%; (b)速度比實(shí)際情況小 5%; (c)速度比實(shí)際情況小 10%; (d)速度比實(shí)際情況大 10%

異常體深度: (a)為1 m; (b)為3 m

(a)水平層狀結(jié)構(gòu)下生成數(shù)據(jù)的相對(duì)介電常數(shù)模型; (b)均一介質(zhì)初始模型的RTM成像結(jié)果。虛線示意分層界面和異常體的實(shí)際位置

圖 7(b1)～(b3)給出對(duì)應(yīng)噪聲水平下的 RTM 成像結(jié)果。在 0dB 信噪比的情況下, 噪聲對(duì)成像結(jié)果基本上沒有影響(圖 7(b1)); 在?15dB 信噪比的情況下, 噪聲對(duì)成像結(jié)果有較明顯的影響, 背景出現(xiàn)很多雜亂的條紋, 但依然能夠較好地識(shí)別異常體(圖 7(b2)); 在?30dB 信噪比的情況下, 只能勉強(qiáng)分辨出異常體(圖 7(b3)); 當(dāng)信噪比低于?30dB 時(shí), 異常體則完全無法成像。結(jié)合前面對(duì)含噪聲原始信號(hào)的分析可以得出, 只要原始數(shù)據(jù)信噪比大于?15dB 左右, RTM成像也能較好地給出異常體信息(圖 7(c)和(d))。該數(shù)值測試結(jié)果表明, RTM 算法對(duì)噪聲有較好的壓制效果。

3 結(jié)論

本文通過對(duì) GPR 的 RTM 成像算法各種影響因素的研究, 分析 GPR 三維 RTM 算法在速度結(jié)構(gòu)誤差和噪聲影響下的魯棒性, 研究結(jié)果可為實(shí)際數(shù)據(jù)采集中數(shù)據(jù)質(zhì)量的控制提供參考。

對(duì)速度結(jié)構(gòu)的分析表明, 在模型介質(zhì)速度與實(shí)際速度結(jié)構(gòu)相差不大的情況下, 成像結(jié)果幾乎沒有偏差; 隨著速度誤差增大, 成像結(jié)果中異常體位置的偏離也變大, 并且對(duì)于相同的速度誤差, 異常體深度越大, 成像結(jié)果偏離越大。所以, 進(jìn)行 RTM成像要盡可能獲得準(zhǔn)確的速度結(jié)構(gòu)。另外, 對(duì)于速度誤差不是特別大的層狀結(jié)構(gòu), 使用簡化的均勻半空間速度模型進(jìn)行 RTM 成像, 可以較好地分辨地下的異常體和層狀結(jié)構(gòu)。

加入高斯白噪聲的信噪比: (a1)和(b1)為 0dB; (a2)和(b2)為?15dB; (a3)和(b3)為?30 dB

對(duì)含噪聲數(shù)據(jù)的分析結(jié)果表明, RTM 算法對(duì)噪聲有較好的壓制效果, 即使接收數(shù)據(jù)中的有效信號(hào)僅能勉強(qiáng)分辨(例如信噪比大于?15dB 時(shí)), RTM 算法依然能較有效地壓制噪聲, 并將對(duì)應(yīng)的結(jié)構(gòu)成像出來。

[1] 曾昭發(fā). 探地雷達(dá)方法原理及應(yīng)用. 北京: 科學(xué)出版社, 2006

[2] Mellett J S. Ground penetrating radar applications in engineering, environmental management, and geo-logy. Journal of Applied Geophysics, 1995, 33: 157– 166

[3] 趙衛(wèi)平, 潘和平, 李清松, 等. 井中雷達(dá)應(yīng)用進(jìn)展. 工程地球物理學(xué)報(bào), 2005, 2(4): 297–303

[4] Zhang Di, Li Jiacun, Liu Shaotang, et al. Multi-frequencies GPR measurements for delineating the shallow subsurface features of the Yushu strike slip fault. Acta Geophys, 2019, 67(2): 501–515

[5] Bakker M A J, Van Der Meer J J M. Structure of a pleistocene push moraine revealed by GPR: the eas-tern Veluwe Ridge, The Netherlands. Geological So-ciety London Special Publications, 2003, 211(1): 143– 151

[6] Trinks I, Johansson B, Gustafsson J, et al. Efficient, large-scale archaeological prospection using a true three-dimensional ground-penetrating radar array sys-tem. Archaeological Prospection, 2010, 17(3): 175–186

[7] Benson A K. Applications of ground penetrating radar in assessing some geological hazards: examples of groundwater contamination, faults, cavities. Journal of Applied Geophysics, 1995, 33: 177–193

[8] Busetti A, Calligaris C, Forte E, et al. Non-invasive methodological approach to detect and characterize high-risk sinkholes in urban cover evaporite karst: integrated reflection seismics, PS-InSAR, Leveling, 3D-GPR and ancillary data. A NE Italian case study. Remote Sensing, 2020, 12(22): 3814

[9] Fisher E, McMechan G A, Annan A P, et al. Examples of reverse-time migration of single-channel, ground- penetrating radar profiles. Geophysics, 1992, 57(4): 577–586

[10] Berkhout A J. Wavefield extrapolation techniques in seismic migration, a tutorial. Geophysics, 1981, 46 (12): 1638–1656

[11] Zhu Jinming, Lines L R. Comparison of Kirchhoff and reverse-time migration methods with applications to prestack depth imaging of complex structures. Geo-physics, 1998, 63(4): 1166–1176

[12] Yoon K, Shin C, Sangyong S U H, et al. 3D reverse-time migration using the acoustic wave equation: an experience with the SEG/EAGE data set. Leading Edge, 2003, 22(1): 38–41

[13] Baysal E, Kosloff D D, Sherwood J W C. Reverse time migration. Geophysics, 1983, 48(11): 1514–1524

[14] McMechan G A. Migration by extrapolation of time-dependent boundary values. Geophysical Prospecting, 1983, 31(3): 413–420

[15] 朱尉強(qiáng). 逆時(shí)偏移成像技術(shù)在探地雷達(dá)數(shù)據(jù)處理中的應(yīng)用[D]. 北京: 北京大學(xué), 2016

[16] 張崇民. 工程隱蔽病害地質(zhì)雷達(dá)三維逆時(shí)偏移成像方法[D]. 濟(jì)南: 山東大學(xué), 2019

[17] 焦翠翠. 基于RTM的礦井富水區(qū)成像研究[D]. 西安: 西安科技大學(xué), 2020

[18] 楊仁虎, 常旭, 劉伊克. 疊前逆時(shí)偏移影響因素分析. 地球物理學(xué)報(bào), 2010, 53(8): 1902–1913

[19] 楊仁虎. 逆時(shí)偏移隨機(jī)噪聲影響分析. 地球物理學(xué)進(jìn)展, 2021, 36(6): 2618–2627

[20] 楊仁虎. 逆時(shí)偏移速度誤差影響分析. 地球物理學(xué)進(jìn)展, 2021, 36(4): 1630–1639

[21] Claerbout J F. Toward a unified theory of reflector mapping. Geophysics, 1971, 36(3): 467–481

[22] Georgakopoulos S V, Birtcher C R, Balanis C A, et al. Higher-order finite-difference schemes for electro-magnetic radiation, scattering, and penetration, part I: theory. IEEE Antennas Propag Mag, 2002, 44(1): 134–142

[23] Fang Jiayuan. Time domain finite difference computa-tion for Maxwell’s equations [D]. Berkeley: Univer-sity of California, 1989

[24] Liu Qinghuo. The PSTD algorithm: a time-domain method requiring only two cells per wavelength. Mic-rowave and Optical Technology Letters, 1997, 15(3): 158–165

Study on Influencing Factors of Pre-Stack Reverse Time Migration Imaging of Ground Penetrating Radar

MA Bowen1,2, HUANG Qinghua1,2,?

1. Department of Geophysics, School of Earth and Space Sciences, Peking University, Beijing 100871; 2. Hebei Hongshan National Observatory on Thick Sediments and Seismic Hazards, Xingtai 054000; ? Corresponding author, E-mail: huangq@pku.edu.cn

To analyze the influence of various possible errors in the actual ground-penetrating radar data processing on the imaging effect of the reverse time migration algorithm, numerical simulation tests are used to analyze the effect of the velocity errors in model and noises in data. For velocity error, models with different velocities and abnormal body depths, and simplified velocity models are used to test the subsequent imaging effect. Results show that with the increase of the velocity error and the depth of the abnormal body, the deviation of the position of the abnormal body also increases. For the layered structure with slight differences in velocity between each layer, imaging using a simplified uniform half-space model can still effectively distinguish the anomalies. Imaging tests are also performed on data with different noise levels. The results show that the reverse time migration (RTM) algorithm is effective on noise suppression. Even if the signal-to-noise ratio of the valid signal in the received data is around ?15 dB, the RTM algorithm can still image the corresponding structure.

ground penetrating radar; 3D reverse time migration; velocity errors; data noise