不依賴(lài)源子波的跨孔雷達(dá)時(shí)間域波形反演

劉四新， 孟旭， 傅磊

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

?

不依賴(lài)源子波的跨孔雷達(dá)時(shí)間域波形反演

劉四新， 孟旭， 傅磊*

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

波形反演是近年來(lái)較熱門(mén)的反演方法，其分辨率可以達(dá)到亞波長(zhǎng)級(jí)別.在波形反演的實(shí)際應(yīng)用中，源子波的估計(jì)十分重要.傳統(tǒng)方法使用反褶積來(lái)估計(jì)源子波并隨著反演過(guò)程更新，該方法在合成數(shù)據(jù)波形反演中效果較好，但在實(shí)際數(shù)據(jù)反演過(guò)程中存在一系列的問(wèn)題.由于實(shí)際數(shù)據(jù)信噪比較低，在源子波估計(jì)過(guò)程中需要大量的人為干涉，且結(jié)果并不一定可靠.本文使用一種基于褶積波場(chǎng)的新型目標(biāo)函數(shù)，令反演過(guò)程不再依賴(lài)源子波.詳細(xì)推導(dǎo)了針對(duì)跨孔雷達(dá)波形反演的梯度及步長(zhǎng)公式，實(shí)現(xiàn)介電常數(shù)和電導(dǎo)率的同步反演.針對(duì)一個(gè)合成數(shù)據(jù)模型同時(shí)反演介電常數(shù)和電導(dǎo)率，結(jié)果表明該方法能夠反演出亞波長(zhǎng)尺寸異常體的形狀和位置.接著，將該方法應(yīng)用到兩組實(shí)際數(shù)據(jù)中，并與基于估計(jì)源子波的時(shí)間域波形反演結(jié)果進(jìn)行比較.結(jié)果表明不依賴(lài)源子波的時(shí)間域波形反演結(jié)果分辨率更高，也更準(zhǔn)確.

波形反演； 鉆孔雷達(dá)； 跨孔； 不依賴(lài)源子波； 介電常數(shù)； 電導(dǎo)率

1 引言

跨孔雷達(dá)是一種高效的淺層地球物理探測(cè)技術(shù)，能夠提供兩個(gè)鉆孔之間的高分辨率地質(zhì)剖面圖像(Pettinelli et al.,2009; Liu and Sato,2002；Liu et al., 2011, 2014).波形反演可以得到精準(zhǔn)的地下介質(zhì)物性參數(shù)，其分辨率達(dá)到亞波長(zhǎng)級(jí)別.結(jié)合跨孔雷達(dá)和波形反演，能得到準(zhǔn)確且高分辨率的地下介質(zhì)物性信息.探地雷達(dá)波形反演起步較晚，早期Moghaddam等(1991)在Tarantol(1984)的全波形反演技術(shù)上對(duì)小目標(biāo)體介電常數(shù)成像做了一些嘗試，但沒(méi)有取得滿(mǎn)意的效果，這主要是因?yàn)槿狈ο闰?yàn)信息來(lái)保證收斂.之后Kuroda等(2005)和Ernst等(2007b)分別使用全波形反演的方法對(duì)跨孔雷達(dá)數(shù)據(jù)成像，并取得了較好的效果.不同的是，Kuroda等的方法只考慮了介電常數(shù)，而Ernst等采用級(jí)聯(lián)更新的辦法同時(shí)考慮介電常數(shù)和電導(dǎo)率.Meles等(2010)提出了矢量全波形反演技術(shù)，并且在反演過(guò)程中同時(shí)迭代介電常數(shù)和電導(dǎo)率.Belina等(2012)對(duì)全波形反演中源的子波估計(jì)部分進(jìn)行了一些分析.國(guó)內(nèi)方面，吳俊軍等(2014)對(duì)全波形反演的理論基礎(chǔ)進(jìn)行了詳細(xì)的推導(dǎo)，使用最速下降法來(lái)更新模型參數(shù)；周輝等(2014)提出了一種不需要提取源子波的波形反演方法，但只考慮介電常數(shù).

在實(shí)際數(shù)據(jù)的波形反演中，由于源子波未知，通常的辦法是在反演過(guò)程中加入源子波這一變量，并隨著反演的進(jìn)行更新.當(dāng)反演結(jié)果與真實(shí)模型相同時(shí)，則估計(jì)得到的源子波與真實(shí)源子波相同.該方法在處理合成數(shù)據(jù)時(shí)十分有效，但在實(shí)際數(shù)據(jù)反演中表現(xiàn)不佳.這可能是因?yàn)閷?shí)際數(shù)據(jù)對(duì)應(yīng)的地下介質(zhì)復(fù)雜度高，且數(shù)據(jù)本身的信噪比較低.在地震波形反演中，已有許多不依賴(lài)于源子波的波形反演方法.原理大致可以分成兩類(lèi)：在頻率域，利用一道參考波形對(duì)數(shù)據(jù)進(jìn)行歸一化處理，在此過(guò)程中自然消除源子波的影響(Lee and Kim，2003；Xu et al.，2006)；另一種為對(duì)實(shí)際數(shù)據(jù)和正演數(shù)據(jù)同時(shí)分別卷乘參考道的實(shí)際和正演數(shù)據(jù)，這樣新數(shù)據(jù)中同時(shí)包含實(shí)際子波和理論子波，從而擺脫源子波的影響(Choi et al.，2005).在探地雷達(dá)波形反演中，周輝等(2014)借助已知源子波、地質(zhì)模型和采集系統(tǒng)的正演數(shù)據(jù)，將實(shí)際數(shù)據(jù)轉(zhuǎn)換成給定源子波的信號(hào)，以避免激發(fā)脈沖的求取.

本文實(shí)現(xiàn)了一種不依賴(lài)于源子波的跨孔雷達(dá)時(shí)間域波形反演.我們首先對(duì)實(shí)際數(shù)據(jù)和一道來(lái)自正演數(shù)據(jù)的參考道做褶積，同樣的，對(duì)正演數(shù)據(jù)和一道來(lái)自實(shí)際數(shù)據(jù)的參考道做褶積.新的目標(biāo)函數(shù)建立在褶積波場(chǎng)之上，理論上，不論正演過(guò)程使用哪種激發(fā)脈沖，新的褶積波場(chǎng)都共用相同的源子波.新目標(biāo)函數(shù)的另一個(gè)重要特征是正演波場(chǎng)使用的源子波對(duì)褶積波場(chǎng)起到低通濾波作用，因此不依賴(lài)源子波的時(shí)間域波形反演很容易實(shí)現(xiàn)多尺度策略.文中使用求導(dǎo)的方法詳細(xì)推導(dǎo)了不依賴(lài)源子波的時(shí)間域波形反演的梯度和步長(zhǎng)公式，實(shí)現(xiàn)了介電常數(shù)和電導(dǎo)率的同步反演.不同于Meles等(2010)根據(jù)擾動(dòng)來(lái)推導(dǎo)梯度，求導(dǎo)法更加簡(jiǎn)潔易懂.之后，我們將該方法應(yīng)用到三組合成數(shù)據(jù)和兩組實(shí)際數(shù)據(jù)中.

2 理論

2.1 目標(biāo)函數(shù)

通常情況下，目標(biāo)函數(shù)為實(shí)際數(shù)據(jù)和正演數(shù)據(jù)之間的l2范數(shù)：

(1)

(2)

g和s分別為格林函數(shù)和源子波，上角標(biāo)for和obs分別代表正演數(shù)據(jù)和實(shí)際數(shù)據(jù)；*號(hào)表示褶積過(guò)程.波形反演即通過(guò)更新模型參數(shù)使(2)式所代表的誤差最小，但在實(shí)際數(shù)據(jù)處理過(guò)程中，源子波sobs必須在反演過(guò)程中估計(jì).

為了避免在反演過(guò)程中對(duì)源子波的估計(jì)，本文實(shí)現(xiàn)一種不依賴(lài)源子波的波形反演.對(duì)于同一個(gè)源下的所有接收數(shù)據(jù)，正演數(shù)據(jù)與一道取自實(shí)際數(shù)據(jù)的參考道褶積；同樣的，實(shí)際數(shù)據(jù)也與一道取自正演數(shù)據(jù)的參考道褶積.新的目標(biāo)函數(shù)為：

(3)

(4)在(4)式中，等號(hào)右端的兩項(xiàng)共用了相同的新源子波sfor*sobs.基于褶積過(guò)程的自然線性特征，源子波的影響被消去了.因此，建立在(4)式基礎(chǔ)上的波形反演不再依賴(lài)于源子波.

2.2 目標(biāo)函數(shù)對(duì)應(yīng)的梯度

(3)式的l2范數(shù)寫(xiě)成點(diǎn)乘形式為：

(5)

我們將(6)式右端第一項(xiàng)寫(xiě)成積分形式：

(7)

令t-τ=ξ，則dτ=-dξ.(7)式可以重新寫(xiě)成：

(8)

根據(jù)附錄A，(8)式重寫(xiě)為：

(9)

(10)

(6)式的第二項(xiàng)采用相同的形式推導(dǎo)：

(11)

其中

t.

(12)

在(11)式和(12)式中，需要注意的是第二個(gè)相關(guān)殘場(chǎng)r″i,j只應(yīng)在第k個(gè)參考道位置后向傳播.

因此，不依賴(lài)源子波的時(shí)間域波形反演對(duì)應(yīng)的目標(biāo)函數(shù)梯度的后向傳播殘場(chǎng)源的第一項(xiàng)為：

(13)

對(duì)應(yīng)在每個(gè)接收器位置j.第二項(xiàng)后向傳播殘場(chǎng)源為：

(14)

只在第k個(gè)參考道位置后向傳播.由于第一項(xiàng)和第二項(xiàng)后向傳播是在同一時(shí)間進(jìn)行的，因此只需要一次正演即可完成.在(13)和(14)式中，?代表互相關(guān).對(duì)于(13)式，只有實(shí)際數(shù)據(jù)的參考道與殘場(chǎng)的每一項(xiàng)進(jìn)行相關(guān)；而在(14)式中，為實(shí)際數(shù)據(jù)和殘場(chǎng)之間逐道之間的相關(guān).在實(shí)際操作中，(14)式對(duì)應(yīng)的第二項(xiàng)后向傳播殘場(chǎng)源被累加在一起作為一道信號(hào)在參考道位置k處反向傳播.

最后，給出最終的梯度表達(dá)式：

(15)

其中

Ts,d,τ(x′,t′)=G(x′,τ-t′,xj,0)(r′+r″),

(16)

T為后向傳播場(chǎng)，(16)式中的格林函數(shù)G表明了其在時(shí)間上的反向傳播.

2.3 迭代步長(zhǎng)的求解

在本文的波形反演中，分別沿著介電常數(shù)和電導(dǎo)率的負(fù)梯度方向?qū)ふ易罴巡介L(zhǎng)，實(shí)現(xiàn)兩者的同步反演.第n次反演對(duì)應(yīng)的目標(biāo)函數(shù)為：

(17)

這里，p代表模型參數(shù)(ε和σ).則第n+1次目標(biāo)函數(shù)為：

(19)

F代表電場(chǎng)E在模型參數(shù)p處的線性因子.由于本文按照負(fù)梯度方向搜索步長(zhǎng)，因此：

(20)

將(20)式代入(18)式：

當(dāng)S(pn+1)最小時(shí)，滿(mǎn)足

(23)

這時(shí)，αn為最佳步長(zhǎng)：

(24)

(25)

(26)

迭代公式為：

(27)

(28)

3 算例

3.1 復(fù)雜模型下介電常數(shù)和電導(dǎo)率同步反演

本節(jié)給出復(fù)雜模型下，同步反演介電常數(shù)和電導(dǎo)率的結(jié)果.模型如圖1所示，圖1a為介電常數(shù)模型，圖1c為電導(dǎo)率模型.該模型分為三層，其中第一層和第三層介質(zhì)相同，相對(duì)介電常數(shù)為5，電導(dǎo)率為0.001 S·m-1.中間層的相對(duì)介電常數(shù)為5.5，電導(dǎo)率為0.0028 S·m-1.中間層埋藏了兩個(gè)管狀異常體，它們的直徑為0.5 m，間隔2 m，相對(duì)介電常數(shù)為7，電導(dǎo)率為0.008 S·m-1.發(fā)射和接收器的位置分別用圓圈和叉表示.初始模型為均勻地層，其相對(duì)介電常數(shù)為5.5，電導(dǎo)率為0.0028 S·m-1，與中間層相同.

圖1b和圖1d分別為同時(shí)反演的介電常數(shù)和電導(dǎo)率結(jié)果，圖中上下兩個(gè)界面清晰，也能夠分辨出兩根管狀異常體的大小和位置.從整體上看，電導(dǎo)率的反演結(jié)果要優(yōu)于介電常數(shù)，表現(xiàn)在兩個(gè)管狀異常體的電導(dǎo)率反演結(jié)果在形狀和數(shù)值上都與真實(shí)模型更接近.但反演結(jié)果中出現(xiàn)了較多的假異常，這是由于本文的波形反演在梯度計(jì)算過(guò)程中包含了數(shù)據(jù)之間的褶積和互相關(guān)處理，提高了反演過(guò)程的非線性程度.加上本節(jié)算例中模型的復(fù)雜度較高，且在反演過(guò)程中介電常數(shù)和電導(dǎo)率之間互相作用，導(dǎo)致源的近場(chǎng)干擾更加強(qiáng)烈.

3.2 實(shí)際數(shù)據(jù)——貴州

數(shù)據(jù)采集地位于中國(guó)西南的貴州省，地貌屬于高原山地，并且境內(nèi)喀斯特地貌發(fā)育豐富，巖溶分布范圍廣泛.采集過(guò)程使用中心頻率為100 MHz的MALA鉆孔雷達(dá)系統(tǒng)，探測(cè)的兩個(gè)鉆孔之間的水平距離為18 m.反演中使用的數(shù)據(jù)包含40組發(fā)射，其中發(fā)射源以1 m等間隔分布在9 m至48 m深度范圍內(nèi).接收器的總覆蓋范圍為0 m至47 m，但每個(gè)源對(duì)應(yīng)的接收器范圍和數(shù)目都不相同.由于缺乏對(duì)應(yīng)地質(zhì)資料，無(wú)法獲得鉆孔所在地詳細(xì)的地質(zhì)信息.初始模型分別由走時(shí)成像和重心頻率下降法(Liu et al., 1998)得到.

圖1 (a)介電常數(shù)真實(shí)模型；(b)介電常數(shù)反演結(jié)果；(c)電導(dǎo)率真實(shí)模型；(d)電導(dǎo)率反演結(jié)果Fig.1 (a) True model of permittivity; (b) Permittivity from source-independent waveform inversion; (c) True model of conductivity; (d) Conductivity from source-independent waveform inversion

圖2 貴州實(shí)際數(shù)據(jù)(a) 介電常數(shù)初始模型； (b) 使用估計(jì)源子波的介電常數(shù)反演結(jié)果； (c) 不依賴(lài)源子波的介電常數(shù)反演結(jié)果； (d) 電導(dǎo)率初始模型； (e) 使用估計(jì)源子波的電導(dǎo)率反演結(jié)果； (f) 不依賴(lài)源子波的電導(dǎo)率反演結(jié)果.Fig.2 Field data of Guizhou(a) Initial model of permittivity; (b) Permittivity from inversion using estimated source wavelets; (c) Permittivity from source-independent inversion; (d) Initial model of conductivity; (e) Conductivity from inversion using source wavelets; (f) Conductivity from source-independent inversion.

在對(duì)實(shí)際數(shù)據(jù)進(jìn)行波形反演之前，首先對(duì)實(shí)際數(shù)據(jù)進(jìn)行了一系列的預(yù)處理，包括提取數(shù)據(jù)、濾波.之后，采用Ernst等(2007a)的辦法對(duì)實(shí)際數(shù)據(jù)做三維校正，另外在實(shí)際數(shù)據(jù)的波形反演過(guò)程中，共進(jìn)行了三次源子波估計(jì).

圖2a和圖2b為反演中使用的介電常數(shù)和電導(dǎo)率初始模型，分別由走時(shí)成像和重心頻率下降法得到.走時(shí)成像得到的介電常數(shù)效果較好，而重心頻率下降法得到的電導(dǎo)率只能反映出地層中異常體的相對(duì)變化，其數(shù)值是不準(zhǔn)確的.圖2b和圖2e為使用反褶積法估計(jì)源子波的時(shí)間域波形反演結(jié)果，圖2c和圖2f為不依賴(lài)源子波的時(shí)間域波形反演結(jié)果.從整體來(lái)看，由波形反演得到的電導(dǎo)率結(jié)果分辨率更高，反映了地下介質(zhì)的精細(xì)結(jié)構(gòu)，而且與介電常數(shù)的反演結(jié)果吻合較好，都表現(xiàn)出左下和右上區(qū)域在數(shù)值上較高的特征.對(duì)比圖2c和圖2f，發(fā)現(xiàn)兩者相對(duì)于圖2a變化較小，這主要是因?yàn)樽邥r(shí)成像技術(shù)十分成熟，其結(jié)果已經(jīng)較為準(zhǔn)確.盡管圖2b和圖2c中均體現(xiàn)了更多的細(xì)節(jié)，這些細(xì)節(jié)難以解釋.我們認(rèn)為圖2c的結(jié)果更為準(zhǔn)確，比如圖2c在深度43 m處的異常與電導(dǎo)率一致，而圖2b中無(wú)這種特征.對(duì)比圖2e和圖2f，兩者都反映了地層特征，如深度18 m處、深度35 m處及深度45 m處.但圖2f的結(jié)果更加清晰，分辨率更高.除此之外，圖2f中在深度13 m位置有一水平薄層，而圖2e中無(wú)體現(xiàn).我們認(rèn)為不依賴(lài)源子波的時(shí)間域波形反演結(jié)果更加準(zhǔn)確可靠，優(yōu)于使用估計(jì)源子波的時(shí)間域波形反演.

3.3 實(shí)際數(shù)據(jù)——岫巖

岫巖位于遼東半島北部，是我國(guó)著名的玉石產(chǎn)地之一.1995年，于岫巖哈達(dá)碑鎮(zhèn)發(fā)現(xiàn)世界最大玉體——‘玉皇’.‘玉皇’高25 m，最大直徑30 m，總體積達(dá)2.4×104m3，總重量約6×104t.本次數(shù)據(jù)采集位置在巨型玉石‘玉皇’腳下，目的是探測(cè)地基穩(wěn)定狀況、采空區(qū)層數(shù)、充填情況及穩(wěn)定狀況.采集過(guò)程使用基于矢量網(wǎng)絡(luò)分析儀的自制鉆孔雷達(dá)天線系統(tǒng)，中心頻率為100 MHz.兩個(gè)鉆孔之間的距離為16 m，如圖3所示.共有32個(gè)發(fā)射源，范圍從29.5 m到67.5 m等間隔分布，每組發(fā)射對(duì)應(yīng)28個(gè)接收器.接收器序列1 m等間隔排列，且隨源的每次下移整體向下移動(dòng)1 m.

圖3 ‘玉皇’現(xiàn)場(chǎng)照片及鉆孔位置示意圖Fig.3 Scene photo of ‘Yuhuang’ and schematic diagram of borehole location

圖3中，開(kāi)裂的巖體為‘玉皇’，人站立的位置為巖石掉塊.根據(jù)已有的地質(zhì)資料，我們得知：‘玉皇’巨型玉體下部0～100 m深度范圍內(nèi)存在五層以上的巷道，分屬不同的礦山.這些巷道中，有的已經(jīng)坍塌，有的尚在使用中.

圖4a和圖4d分別為走時(shí)成像和重心頻率下降法得到的介電常數(shù)和電導(dǎo)率分布圖像.圖4b和圖4e為使用估計(jì)源子波的時(shí)間域波形反演結(jié)果，圖4c和圖4f為不依賴(lài)源子波的時(shí)間域波形反演結(jié)果.由于該組實(shí)際數(shù)據(jù)的信噪比不高，在使用反褶積估計(jì)源子波時(shí)僅使用了部分?jǐn)?shù)據(jù)(深度62.5～67.5 m).在使用估計(jì)源子波的整個(gè)反演過(guò)程中共進(jìn)行三次源子波估計(jì)，但后兩次僅用來(lái)修正估計(jì)源子波的振幅.由于只使用了部分?jǐn)?shù)據(jù)估計(jì)源子波，在反演過(guò)程中缺乏對(duì)源子波相位的修正，圖4b中深度67.5 m處的異?？赡苁怯晒烙?jì)源子波和真實(shí)源子波之間的相位差異造成的.圖4e的電導(dǎo)率反演結(jié)果的最大值超過(guò)了1 S·m-1，最小值低于0.01 mS·m-1，這種不合理反映了估計(jì)源子波在振幅上是不準(zhǔn)確的.這種不準(zhǔn)確表現(xiàn)在兩個(gè)互相矛盾的方面：以深度61.5 m處為分界線，上半部分表明估計(jì)源子波的振幅較小，因而整體電導(dǎo)率結(jié)果較大；下半部分表明估計(jì)源子波的振幅較大，所以電導(dǎo)率的結(jié)果較小.造成這種現(xiàn)象的原因在于數(shù)據(jù)質(zhì)量不高，反褶積法無(wú)法得到可靠的源子波.圖4f與圖4e在整體上表現(xiàn)一致，但圖4f在數(shù)值上更為合理，這是因?yàn)橄嗽吹挠绊?圖4f在深度61.5至67.5 m范圍內(nèi)表現(xiàn)出電導(dǎo)率較高的層狀結(jié)構(gòu)，圖4e在相同位置也有所體現(xiàn).因此，在實(shí)際數(shù)據(jù)質(zhì)量不高的情況下，使用估計(jì)源子波的波形反演難以得到好的反演結(jié)果，不依賴(lài)源子波的波形反演具有更大的優(yōu)勢(shì).

4 總結(jié)與討論

本文實(shí)現(xiàn)了一種不依賴(lài)源子波的時(shí)間域跨孔雷達(dá)波形反演，其原理是建立一個(gè)新的基于褶積波場(chǎng)的目標(biāo)函數(shù).在該目標(biāo)函數(shù)中，通過(guò)褶積使實(shí)際數(shù)據(jù)和正演數(shù)據(jù)共用相同的新源項(xiàng)，從而消去了源子波的影響.通過(guò)求導(dǎo)的方式詳細(xì)推導(dǎo)了目標(biāo)函數(shù)對(duì)應(yīng)的梯度和步長(zhǎng)公式，并引入地震波形反演中的虛擬源概念.

為了驗(yàn)證文中不依賴(lài)源子波的波形反演效果，同時(shí)對(duì)一個(gè)合成數(shù)據(jù)模型反演介電常數(shù)和電導(dǎo)率，結(jié)果表明該方法能夠重建亞波長(zhǎng)尺寸的異常體和地層結(jié)構(gòu).最后，我們將該方法應(yīng)用到兩組實(shí)際數(shù)據(jù)中.其中貴州實(shí)際數(shù)據(jù)質(zhì)量較好，使用估計(jì)源子波和不依賴(lài)源子波兩種波形反演均能得到可靠的反演結(jié)果，但通過(guò)對(duì)比我們認(rèn)為不依賴(lài)源子波的波形反演效果更佳，這可能是估計(jì)得到的源子波與實(shí)際源子波仍存在差異造成的.而岫巖實(shí)際數(shù)據(jù)的信號(hào)質(zhì)量較差，使用反褶積法估計(jì)源子波具有很大的難度，最終的反演結(jié)果較差.在這種情況下，更能體現(xiàn)出不依賴(lài)源子波波形反演的價(jià)值.

盡管不依賴(lài)源子波的波形反演在理論上更加先進(jìn)，這種方法也帶來(lái)一些問(wèn)題.最主要的是該方法在梯度計(jì)算過(guò)程中需要進(jìn)行褶積和互相關(guān)運(yùn)算，較大地提高了反演的非線性，因此需要更好的初始模型和反演策略來(lái)保證收斂.由于在頻率域更易于實(shí)現(xiàn)多尺度策略，未來(lái)我們考慮研究頻率域不依賴(lài)源子波的波形反演.

附錄A 虛擬源原理

對(duì)于Maxwell方程：

(A1)

在時(shí)間域，我們可以簡(jiǎn)寫(xiě)成：

(A2)

(A2)式兩端分別對(duì)模型參數(shù)p求導(dǎo)，并忽略磁場(chǎng)項(xiàng)：

(A3)

整理可得：

(A4)

其中

(A5)

v為虛擬源(Pratt et al., 1998; Shin et al., 2001)，對(duì)于不同的模型參數(shù)(ε，σ)具有不同的表達(dá)式；M-1等同于格林函數(shù)g.對(duì)于時(shí)間域Maxwell方程而言，有

圖4 岫巖實(shí)際數(shù)據(jù)(a) 介電常數(shù)初始模型； (b) 使用估計(jì)源子波的介電常數(shù)反演結(jié)果； (c) 不依賴(lài)源子波的介電常數(shù)反演結(jié)果； (d) 電導(dǎo)率初始模型； (e) 使用估計(jì)源子波的電導(dǎo)率反演結(jié)果； (f) 不依賴(lài)源子波的電導(dǎo)率反演結(jié)果.Fig.4 Field data of Xiuyan(a) Initial model of permittivity; (b) Permittivity from inversion using source wavelets; (c) Permittivity from source-independent inversion; (d) Initial model of conductivity; (e) Conductivity from inversion using source wavelets; (f) Conductivity from source-independent inversion.

vε=?tE,

(A6)

vσ=E.

(A7)

對(duì)于(A2)式，如果忽略磁場(chǎng)項(xiàng)，寫(xiě)成包含格林函數(shù)的形式：

E=g*J,

(A8)

則對(duì)比(A2)(A4)(A5)(A8)可以得出：

(A9)

(A10)

Belina F A, Irving J, Ernst J R, et al. 2012. Waveform inversion of crosshole georadar data: influence of source wavelet variability and the suitability of a single wavelet assumption.IEEETransactionsonGeoscienceandRemoteSensing, 50(11): 4610-4625. Choi Y, Shin C, Min D J, et al. 2005. Efficient calculation of the steepest descent direction for source-independent seismic waveform inversion: An amplitude approach.JournalofComputationalPhysics, 208(2): 455-468.

Ernst J R, Green A G, Maurer H, et al. 2007a. Application of a new 2D time-domain full-waveform inversion scheme to crosshole radar data.Geophysics, 72(5): J53-J64.

Ernst J R, Maurer H, Green A G, et al. 2007b. Full-waveform inversion of crosshole radar data based on 2-D finite-difference time-domain solutions of Maxwell′s equations.IEEETransactionsonGeoscienceandRemoteSensing, 45(9): 2807-2828. Kuroda S, Takeuchi M, Kim H J. 2005. Full waveform inversion algorithm for interpreting cross-borehole GPR data.∥ 75th Annual International Meeting, SEG, Expanded Abstracts: 1176-1180.

Lee K H, Kim H J. 2003. Source-independent full-waveform inversion of seismic data.Geophysics, 68(6): 2010-2015. Liu L B, Lane J W, Quan Y L. 1998. Radar attenuation tomography using the centroid frequency downshift method.JournalofAppliedGeophysics, 40(1-3): 105-116.

Liu S X, Sato M. 2002. Electromagnetic logging technique based on borehole radar.IEEETransactionsonGeoscienceandRemoteSensing, 40(9): 2083-2092.

Liu S X, Wu J J, Zhou J F, et al. 2011. Numerical simulations of borehole radar detection for metal ore.IEEETransactionsonGeoscienceandRemoteSensing, 8(2): 308-312.

Liu S X, Lei L L, Fu L, et al. 2014. Application of pre-stack reverse time migration based on FWI velocity estimation to ground penetrating radar data.JournalofAppliedGeophysics, 107: 1-7. Meles G A, van der Kruk J, Greenhalgh S A, et al. 2010. A new vector waveform inversion algorithm for simultaneous updating of conductivity and permittivity parameters from combination crosshole/borehole-to-surface GPR data.IEEETransactionsonGeoscienceandRemoteSensing, 48(9): 3391-3407.

Moghaddam M, Chew W C, Oristaglio M. 1991. Comparison of the born iterative method and Tarantola′s method for an electromagnetic time-domain inverse problem.InternationalJournalofImagingSystemsandTechnology, 3(4): 318-333.

Pettinelli E, Di Matteo A, Mattei E, et al. 2009. GPR response from buried pipes: Measurement on field site and tomographic reconstructions.IEEETransactionsonGeoscienceandRemoteSensing, 47(8): 2639-2645.

Pratt G, Shin C, Hicks G L. 1998. Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion.GeophysicalJournalInternational, 133(2): 341-362. Shin C, Jang S, Min D J. 2001. Improved amplitude preservation for prestack depth migration by inverse scattering theory.GeophysicalProspecting, 49(5): 592-606.

Tarantola A. 1984. Inversion of seismic reflection data in the acoustic approximation.Geophysics, 49(8): 1259-1266.

Wu J J, Liu S X, Li Y P, et al. 2014. Study of cross-hole radar tomography using full-waveform inversion.ChineseJ.Geophys. (in Chinese), 57(5): 1623-1635, doi: 10.6038/cjg20140525.

Xu K, Greenhalgh S A, Wang M Y. 2006. Comparison of source-independent methods of elastic waveform inversion.Geophysics, 71(6): R91-R100.

Zhou H, Chen H M, Li Q Q, et al. 2014. Waveform inversion of ground-penetrating radar data without needing to extract exciting pulse.ChineseJ.Geophys.(in Chinese), 57(6): 1968-1976, doi: 10.6038/cjg20140627.

附中文參考文獻(xiàn)

吳俊軍, 劉四新, 李彥鵬等. 2014. 跨孔雷達(dá)全波形反演成像方法的研究. 地球物理學(xué)報(bào), 57(5): 1623-1635, doi: 10.6038/cjg20140525.

周輝, 陳漢明, 李卿卿等. 2014. 不需提取激發(fā)脈沖的探地雷達(dá)波形反演方法. 地球物理學(xué)報(bào), 57(6): 1968-1976, doi: 10.6038/cjg20140627.

(本文編輯 何燕)

Source-independent time-domain waveform inversion of cross-hole GPR data

LIU Si-Xin, MENG Xu, FU Lei*

CollegeofGeo-ExplorationScienceandTechnology,JilinUniversity,Changchun130026,China

As a highly efficient near-surface geophysical technology, cross-hole GPR can provide high-resolution images between two boreholes. Waveform inversion, of which the resolution reaches sub-wavelength, can reconstruct accurate subsurface dielectric parameters such as permittivity and conductivity. We can get accurate and high-resolution subsurface dielectric information by combining waveform inversion and cross-hole GPR data. In the waveform inversion of field data, the source wavelet is unknown. Normally, the source wavelet can be estimated by using a deconvolution method. In this method, the source wavelet is a new unknown parameter added in the inversion and updated with iteration. When the results of inversion are same as the true models, the estimated source wavelet is same as the true source wavelet. This method is useful in the inversion of synthetic data, but it does not perform well in the field data. Lots of intervention are needed to choose the best source wavelet after the estimation. The reason may be that the subsurface medium of field data is highly complex and the signal-to-noise ratio of the field data is low.In this paper, we realize a source-independent time-domain waveform inversion. Firstly, the observed wavefields are convolved with a reference trace of the modeled wavefield, then the modeled wavefields are convolved with a reference trace of observed wavefield. A new object function is based on the convolved wavefields. In theory, no matter which source wavelet is used in the forward modeling. The source wavelet of the observed and the modeled wavefields are equally convolved with both terms in the object function, so that the effect of the source wavelet is removed. Another important feature of this object function is that the modeled wavefields act as a low-pass filter of the wavefield based on the frequency range of the source used for modeling. Though it is very easy to employ a frequency-selection strategy, we don′t discuss the strategy here. For each shot gather, we always choose the trace which is nearest to the antenna as the reference trace in this study.To check the effect of our algorithm, we simulate three different synthetic data: 1) Single cylindrical body of permittivity. 2) Single cylindrical body of conductivity. And 3) layered media with multiple embedded cylindrical inclusions of permittivity and conductivity. Finally, we apply the source-independent waveform inversion to the field data of Guizhou and Xiuyan.The first two results of synthetic data show that our algorithm can provide high-resolution images. When permittivity and conductivity are simultaneously updated, the results of synthetic data are good. The layers are well reconstructed and the two pipes can be clearly distinguished. But there are many false features which are caused by the artifacts of sources. We have to do convolution and crosscorrelation to compute the gradients. These convolution and crosscorrelation operations increase the nonlinearity of the inversion. The artifacts of the sources become stronger after the convolution and crosscorrelation is hard to remove when the complex models are of a small scale. Then we apply the source-independent waveform inversion to the field data of Guizhou and Xiuyan. The signal-to-noise ratio of Guizhou field data is high, so both the source-independent and source-estimated waveform inversion can obtain accurate results. We think the results of source-independent waveform inversion are better. The reason may be the difference between the estimated source wavelet and the true source wavelet, though they are very similar. However, it is difficult to get the ideal estimated source wavelet of Xiuyan field data because its signal-to-noise ratio is low. So the results of the Xiuyan field data can better reflect the value of our source-independent waveform inversion.

Waveform inversion; Borehole radar; Cross-hole; Source-independent; Dielectric constant; Conductivity

10.6038/cjg20161210.

國(guó)家重點(diǎn)研發(fā)計(jì)劃(2016YFC0600505),國(guó)家高技術(shù)研究發(fā)展計(jì)劃(2013AA064603), 國(guó)家自然科學(xué)基金(3A515AV44423), 博士后科學(xué)基金(801151010423)共同資助.

劉四新，男，1966年生，日本東北大學(xué)工學(xué)博士，博士生導(dǎo)師，主要從事探地雷達(dá)、鉆孔雷達(dá)及電磁波測(cè)井等的方法理論和應(yīng)用方面的研究. E-mail：liusixin@jlu.edu.cn

*通訊作者 傅磊，男，1985年生，講師，主要從事地球物理正反演及成像研究.E-mail: leifu@jlu.edu.cn

10.6038/cjg20161210

P631

2016-02-01，2016-11-02收修定稿

劉四新， 孟旭， 傅磊. 2016. 不依賴(lài)源子波的跨孔雷達(dá)時(shí)間域波形反演. 地球物理學(xué)報(bào)，59(12):4473-4482,

Liu S X, Meng X, Fu L. 2016. Source-independent time-domain waveform inversion of cross-hole GPR data.ChineseJ.Geophys. (in Chinese),59(12):4473-4482,doi:10.6038/cjg20161210.