基于對(duì)數(shù)目標(biāo)函數(shù)的跨孔雷達(dá)頻域波形反演

孟旭， 劉四新， 傅磊， 王憲楠， 劉新彤， 王文天， 蔡佳琪

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

基于對(duì)數(shù)目標(biāo)函數(shù)的跨孔雷達(dá)頻域波形反演

孟旭， 劉四新*， 傅磊， 王憲楠， 劉新彤， 王文天， 蔡佳琪

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

摘要波形反演在探地雷達(dá)領(lǐng)域的應(yīng)用已有十余年歷史，但絕大部分算例屬于時(shí)間域波形反演.頻率域波形反演由于能夠靈活地選擇迭代頻率并可以使用不同類(lèi)型的目標(biāo)函數(shù)，因而更加多樣化.本文的頻率域波形反演基于時(shí)間域有限差分(FDTD)法，采用對(duì)數(shù)目標(biāo)函數(shù)，可在每一次迭代過(guò)程中同時(shí)或者單獨(dú)反演介電常數(shù)和電導(dǎo)率.文中詳細(xì)推導(dǎo)了頻率域波形反演的理論公式，給出對(duì)數(shù)目標(biāo)函數(shù)下的梯度表達(dá)式，并使用離散傅氏變換(DFT)實(shí)現(xiàn)數(shù)據(jù)的時(shí)頻變換，能夠有效地減少大模型反演的內(nèi)存需求.在后向殘場(chǎng)源的時(shí)頻域轉(zhuǎn)換過(guò)程中，提出僅使用以當(dāng)前頻點(diǎn)為中心的一個(gè)窄帶數(shù)據(jù)，可以消除高頻無(wú)用信號(hào)的干擾，獲得可靠的反演結(jié)果.為加速收斂，采用每迭代十次則反演頻率跳躍一定頻帶寬度的反演策略.實(shí)驗(yàn)證明適當(dāng)?shù)念l率跳躍能夠在不降低分辨率的基礎(chǔ)上有效地提高反演效率.通過(guò)兩組不同情形下合成數(shù)據(jù)反演的分析對(duì)比，證明基于對(duì)數(shù)目標(biāo)函數(shù)的波形反演結(jié)果準(zhǔn)確可靠.最后，將該方法應(yīng)用到一組實(shí)際數(shù)據(jù)，得到較好的反演結(jié)果.

關(guān)鍵詞波形反演； 對(duì)數(shù)目標(biāo)函數(shù)； 介電常數(shù)； 電導(dǎo)率； 鉆孔雷達(dá)

In this paper, waveform inversion is implemented by means of a finite-difference time-domain solution of Maxwell′s equations and a logarithmic objective function is applied. Permittivity and conductivity can be updated simultaneously or separately at each iterative step. The derivation process of the formulas is described in detail and we show the specific expression of gradient under the logarithmic objective function. It is important to note that the gradient formula of the frequency domain waveform inversion is different from the gradient formula of the time domain waveform inversion. The reason is that the cost function is essentially different. Discrete Fourier transform (DFT) is applied to transform time domain data into the frequency domain, which only increases a few calculations in the inversion. The method can greatly reduce memory requirement when the inverted model is on a large scale. When transforming the back-residual source, we present that only a narrow-band data whose center is the current frequency is used. The method can effectively reduce the interference of high frequency information, so reliable inversion results can be obtained.

In order to accelerate the convergence, a special frequency strategy is applied. The inversion frequency skips a number of bandwidth after every ten times iterative step. Results turn out that the strategy can efficiently improve the inversion efficiency and does not influence the resolution. To verify the effect of our algorithm, we test it on two different synthetic data. Then we compare with the results of conventional waveform inversion: (1) Two small cylindrical bodies of permittivity in a layered medium. (2) A layered medium with multiple embedded cylindrical inclusions of permittivity and conductivity. Finally, we apply the logarithmic algorithm to field data.

The results on synthetic data show that the logarithmic waveform inversion is better than the conventional waveform inversion. When permittivity and conductivity are inverted at the same time, the results of logarithmic waveform inversion can accurately reconstruct both shape and location of the inclusions. This is because that the amplitude of the residual is normalized by the amplitude of the modeled wavefield. The inversion results of the field data also prove the practical value of logarithmic waveform inversion.

1引言

探地雷達(dá)層析成像在地質(zhì)學(xué)、工程學(xué)和水文地質(zhì)學(xué)中應(yīng)用廣泛(Pettinelli et al.,2009; Liu and Sato, 2002; Liu et al., 2011,2014；傅磊等，2014)，但傳統(tǒng)的層析成像技術(shù)(比如走時(shí)成像和衰減成像)由于只能利用一部分的信號(hào)信息，所提供的結(jié)果分辨率有限(Williamson and Worthington,1993).波形反演技術(shù)能夠提供亞波長(zhǎng)級(jí)分辨率的地下介質(zhì)分布圖像，因此在近些年深受矚目.Kuroda等(2005)使用全波形反演的方法對(duì)跨孔雷達(dá)數(shù)據(jù)成像，但只考慮介電常數(shù).Ernst等(2007a)的波形反演方法能夠通過(guò)級(jí)聯(lián)的方式既考慮介電常數(shù)又考慮電導(dǎo)率.Meles等(2010)提出了矢量全波形反演技術(shù)，實(shí)現(xiàn)了介電常數(shù)和電導(dǎo)率的同步迭代.在實(shí)際數(shù)據(jù)處理方面，Cai等(1996)使用程函方程對(duì)實(shí)際數(shù)據(jù)進(jìn)行反演；Zhou等(2007)的時(shí)間域重建法屬于波形反演方法，使用實(shí)際數(shù)據(jù)反演了介電常數(shù)和電導(dǎo)率；Ernst等(2007b)將全波形反演應(yīng)用到跨孔雷達(dá)實(shí)際數(shù)據(jù)處理中，對(duì)源子波估計(jì)和實(shí)際數(shù)據(jù)的三維校正等問(wèn)題作了一定分析.國(guó)內(nèi)方面，王兆磊等(2007)使用雷達(dá)數(shù)據(jù)反演了地下介質(zhì)的物性參數(shù)；周輝等(2014)提出了一種不需要提取激發(fā)脈沖的探地雷達(dá)波形反演方法.

目前探地雷達(dá)波形反演方法大部分屬于時(shí)間域波形反演，但在地震波形反演中，頻率域波形反演占據(jù)重要地位.自從Pratt等(1998)將后向傳播技術(shù)應(yīng)用在大尺度地質(zhì)模型的頻率域波形反演后，該技術(shù)被廣泛使用在頻率域地震波形反演中.相對(duì)于時(shí)間域波形反演，頻率域波形反演具有多變的頻率選擇策略并能夠選擇多種目標(biāo)函數(shù).Song等(1995)的串行反演策略是將低頻數(shù)據(jù)波形反演結(jié)果作為相鄰高頻數(shù)據(jù)波形反演結(jié)果的初始模型，該方法因具有較高的穩(wěn)定性被廣泛應(yīng)用于波形反演中.Pratt(1999), Pratt和Shipp (1999)將有效頻段內(nèi)的頻率域數(shù)據(jù)按從低頻到高頻分組，較低頻組的反演結(jié)果作為相鄰較高頻組的初始模型，該策略為組間串行反演策略，有利于壓制噪聲和反演假象.Sirgue和Pratt(2004)對(duì)近層狀模型反射地震數(shù)據(jù)多頻反演的頻率間隔選擇策略做了研究.目標(biāo)函數(shù)方面，即使在地震波形反演中，很長(zhǎng)一段時(shí)間內(nèi)反演目標(biāo)函數(shù)一直在使用理論波場(chǎng)與觀測(cè)波場(chǎng)殘差的l2范數(shù)(Tarantola, 1984，1986; Pratt, 1990; Pratt and Worthington,1990; Pratt et al., 1998).Shin和Min(2006)提出了基于自然對(duì)數(shù)波場(chǎng)的反演目標(biāo)函數(shù)，即對(duì)數(shù)目標(biāo)函數(shù)，并取得了較好的效果.Chung等(2007)使用混合l1/l2范數(shù)的Huber函數(shù)作為頻率域波形反演的目標(biāo)函數(shù).

本文在Maxwell方程的頻率域表達(dá)式基礎(chǔ)上，詳細(xì)推導(dǎo)了基于對(duì)數(shù)目標(biāo)函數(shù)的頻率域波形反演理論公式，并指出其梯度公式與時(shí)間域波形反演梯度公式之間不是簡(jiǎn)單對(duì)應(yīng)關(guān)系.盡管這兩者在表達(dá)式上是對(duì)應(yīng)的，但不能僅通過(guò)對(duì)公式的時(shí)頻轉(zhuǎn)換得到.反演中的正演過(guò)程使用時(shí)間域有限差分(FDTD)法，之后利用離散傅里葉變換(DFT)將時(shí)間域數(shù)據(jù)轉(zhuǎn)換為頻率域數(shù)據(jù)，并在頻率域計(jì)算梯度.相對(duì)于時(shí)間域波形反演，使用DFT技術(shù)不需要保存所有時(shí)間步長(zhǎng)的電場(chǎng)值，能夠極大地減少大尺度模型反演過(guò)程中的計(jì)算機(jī)內(nèi)存需求.在計(jì)算后向傳播場(chǎng)時(shí)，殘場(chǎng)源是在頻率域計(jì)算的，因此需要將殘場(chǎng)源轉(zhuǎn)換到時(shí)間域.為了避免高頻無(wú)用信號(hào)的干擾，我們提出在殘場(chǎng)源的時(shí)頻域轉(zhuǎn)換過(guò)程中，只使用以當(dāng)前反演頻率為中心的一個(gè)窄帶內(nèi)頻域數(shù)據(jù)，有利于提高穩(wěn)定性.文中反演頻率的選擇采用串行反演策略，為了加速收斂，每迭代若干次后將反演頻率增加10 MHz.由于在跨孔模式下每組發(fā)射天線對(duì)應(yīng)的接收數(shù)據(jù)是相互獨(dú)立的，我們采用MATLAB分布式計(jì)算引擎(Matlab Distributed Computing Engine，MDCE)，通過(guò)建立一個(gè)局域網(wǎng)機(jī)群，實(shí)現(xiàn)正演的并行運(yùn)算，能夠明顯地提高反演效率.

通過(guò)兩組合成數(shù)據(jù)來(lái)檢驗(yàn)本文波形反演方法的效果.第一組合成數(shù)據(jù)只反演介電常數(shù)，用來(lái)驗(yàn)證文中提出的反演策略；第二組模型同時(shí)反演介電常數(shù)和電導(dǎo)率，結(jié)果證明基于對(duì)數(shù)目標(biāo)函數(shù)的頻率域波形反演比基于傳統(tǒng)目標(biāo)函數(shù)的頻率域波形反演更加準(zhǔn)確.在實(shí)際數(shù)據(jù)處理方面，仍需要考慮數(shù)據(jù)的三維校正和源子波估計(jì)等問(wèn)題，在此基礎(chǔ)上給出了實(shí)際數(shù)據(jù)的反演結(jié)果.

2理論

2.1正演問(wèn)題的頻域表達(dá)形式

在探地雷達(dá)地球物理問(wèn)題中，使用介電常數(shù)和電導(dǎo)率來(lái)描述地下介質(zhì)的分布情況，并假設(shè)磁導(dǎo)率是均勻不變.Maxwell方程可以表示為如下矩陣形式：

(1)

(2)

(3)

(4)

2.2對(duì)數(shù)目標(biāo)函數(shù)及梯度計(jì)算

不同于Shin和Min(2006)通過(guò)對(duì)目標(biāo)函數(shù)和聲波方程求導(dǎo)來(lái)推導(dǎo)梯度的辦法，本文的理論推導(dǎo)基于Meles等(2010)的波形反演理論(吳俊軍等，2014)與對(duì)數(shù)目標(biāo)函數(shù)的聯(lián)合.(5)式為頻率域波形反演中使用最廣泛的傳統(tǒng)目標(biāo)函數(shù)，它是正演模擬數(shù)據(jù)與實(shí)際數(shù)據(jù)差的l2范數(shù)：

(5)

本文使用的目標(biāo)函數(shù)為基于自然對(duì)數(shù)的正演模擬數(shù)據(jù)與實(shí)際數(shù)據(jù)差的l2范數(shù)：

(7)

+iθj(ω)+2πnj],

(8)

(9)

nm，nj和nf為計(jì)算得到的整數(shù)值.則

(10)

正如(10)式中表明的，需要對(duì)相位進(jìn)行反纏繞來(lái)獲得有意義的相位信息，但這個(gè)過(guò)程通常難以實(shí)現(xiàn).假定nm+nj=nf.該假定隱含了正演模型必須與真實(shí)模型接近這個(gè)條件，從而反纏繞過(guò)后的實(shí)際數(shù)據(jù)的相位與正演數(shù)據(jù)的相位具有相同的度數(shù)(ShinandMin，2006).則目標(biāo)函數(shù)寫(xiě)成：

(11)

從(11)式發(fā)現(xiàn)，目標(biāo)函數(shù)為相位誤差及對(duì)數(shù)振幅誤差的l2范數(shù)形式，這也是對(duì)數(shù)目標(biāo)函數(shù)的特點(diǎn)，能夠?qū)崿F(xiàn)振幅和相位的自然分離.全波形反演層析成像的目的是尋找目標(biāo)函數(shù)S(ε,σ)最小時(shí)對(duì)應(yīng)的介電常數(shù)ε和電導(dǎo)率σ空間分布.本文使用最速下降法來(lái)尋找目標(biāo)函數(shù)最小.為了得到目標(biāo)函數(shù)的梯度，寫(xiě)出擾動(dòng)系統(tǒng)下的Maxwell方程組：

(12)

(12)式減去(1)式得：

(13)

(14)

在空間域，這是一個(gè)關(guān)于擾動(dòng)正演電場(chǎng)的方程.將(14)式代入(3)式中：

(15)

(16)

(17)

則

(18)

目標(biāo)函數(shù)梯度的推導(dǎo)使用擾動(dòng)目標(biāo)函數(shù)的一階近似：

+O(δ ε2,δ σ2),

(19)

擾動(dòng)目標(biāo)函數(shù)為：

(21)

其中，Rn(ε,σ)代表高階項(xiàng).將(21)式代入(20)式中展開(kāi)化簡(jiǎn)得到：

(22)

(23)其中

(24)

(25)

(26)

(25)式為最終的梯度表達(dá)式.需要注意的是，梯度公式中的Re表示取復(fù)數(shù)的實(shí)部.這是因?yàn)樵谕茖?dǎo)梯度過(guò)程中，(22)式與梯度相關(guān)的項(xiàng)在簡(jiǎn)化時(shí)虛部互相抵消了.盡管推導(dǎo)的方式不同，式(25)與頻率域地震波形反演表達(dá)類(lèi)似，但如果僅根據(jù)時(shí)間域波形反演的梯度公式做傅里葉變換得到頻率域的梯度公式，其表達(dá)是不正確的.本文(23)式的推導(dǎo)過(guò)程中做了一些省略，具體情形參考Meles等(2010)及吳俊軍等(2014)的論文.

2.3后向傳播殘場(chǎng)源

由于殘場(chǎng)源首先在頻率域計(jì)算，之后經(jīng)傅里葉變換至?xí)r間域作為后向傳播FDTD中的源項(xiàng).與基于傳統(tǒng)目標(biāo)函數(shù)的頻率域波形反演不同，不能直接使用正演數(shù)據(jù)和實(shí)際數(shù)據(jù)的頻域場(chǎng)值來(lái)計(jì)算后向殘場(chǎng)源，這是因?yàn)?24)式中包含的除法會(huì)在計(jì)算過(guò)程放大原來(lái)無(wú)用的高頻噪聲.本文的正演在時(shí)間域進(jìn)行，反演在頻率域進(jìn)行，因此后向傳播過(guò)程中需要經(jīng)歷二次傅里葉變換和一次反傅里葉變換(后向殘場(chǎng)源的正反變換和整個(gè)模型空間電場(chǎng)值的正變換).

圖1a顯示了一道實(shí)際數(shù)據(jù)及與它對(duì)應(yīng)的正演數(shù)據(jù)的時(shí)間域波形，圖1b為兩者的頻譜.從圖1b中可以看出，信號(hào)的有效頻率范圍大致在250 MHz以下，也就是說(shuō)高于250 MHz的信息為無(wú)用的噪聲.使用圖1a中的正演波形和實(shí)際波形對(duì)應(yīng)的全部頻域數(shù)據(jù)，經(jīng)(24)式計(jì)算并做反傅里葉變換，得到圖1c中后向殘場(chǎng)源的時(shí)間域波形.圖1d中實(shí)線為圖1c對(duì)應(yīng)的頻譜，虛線為(24)式直接計(jì)算的結(jié)果.對(duì)比發(fā)現(xiàn)使用所有頻域數(shù)據(jù)的結(jié)果完全喪失了本應(yīng)具備(直接計(jì)算結(jié)果)的特征，這是因?yàn)?24)式包含的除法，放大了250 MHz之后的噪聲信息.這些無(wú)意義的噪聲在頻譜中占據(jù)了主導(dǎo)地位，經(jīng)過(guò)反傅里葉變換得到圖1c中的無(wú)效波形.圖1f為本文使用的有限點(diǎn)頻譜轉(zhuǎn)換方法：選擇以當(dāng)前反演頻率為中心的，帶寬在20 MHz左右范圍內(nèi)的頻率數(shù)據(jù)，且指定其他頻率數(shù)據(jù)極小(10-8+10-8j).之后反傅里葉變換得到時(shí)間域的后向殘場(chǎng)源.圖1e為以80 MHz為中心的后向殘場(chǎng)源.圖1f中實(shí)線為使用頻率域數(shù)據(jù)直接計(jì)算的頻域后向殘場(chǎng)源；點(diǎn)虛線、短虛線和長(zhǎng)虛線分別為以80 MHz，130 MHz和180 MHz為中心計(jì)算的時(shí)間域后向殘場(chǎng)源的頻譜.從圖1f中可以看出，我們的方法能夠有效保證后向傳播殘場(chǎng)源在對(duì)應(yīng)頻率范圍內(nèi)的正確性.

圖1為了顯示方便，我們最多給出了300 MHz內(nèi)的頻譜.需要注意的是：圖1b中300 MHz之后的頻域數(shù)據(jù)極??；圖1d和圖1f中直接計(jì)算得到的殘場(chǎng)源頻譜隨著頻率的增加急劇增大，之后的頻域數(shù)值比我們關(guān)注的300 MHz以?xún)?nèi)的頻域數(shù)值高出好幾個(gè)數(shù)量級(jí)(對(duì)比圖1d中虛線和圖1f中實(shí)線，兩者為顯示范圍不同下的同一數(shù)據(jù))，這使得傳統(tǒng)帶通濾波方法壓制高頻噪聲的效率很低.因此本文提出利用有限頻點(diǎn)來(lái)計(jì)算后向殘場(chǎng)源，之后的合成數(shù)據(jù)以及實(shí)際數(shù)據(jù)結(jié)果也證明了該方法的有效性.

圖1　后向傳播殘場(chǎng)源的變換(a)實(shí)線為正演數(shù)據(jù)，虛線為實(shí)際數(shù)據(jù)； (b) 實(shí)線為正演數(shù)據(jù)的頻譜，虛線為實(shí)際數(shù)據(jù)的頻譜； (c) 使用所有頻點(diǎn)計(jì)算得到的后向殘場(chǎng)源； (d) 實(shí)線為(c)對(duì)應(yīng)的頻譜，虛線為直接由式(24)計(jì)算的頻譜.左側(cè)y軸對(duì)應(yīng)實(shí)線振幅，右側(cè)y軸為虛線振幅； (e) 以80 MHz為中心的窄帶數(shù)據(jù)計(jì)算的時(shí)間域后向殘場(chǎng)源； (f) 中心頻率為80 MHz、130 MHz、180 MHz及直接使用公式(24)計(jì)算的后向殘場(chǎng)源頻譜.Fig.1　Transformation of the back-residual source(a) Solid line shows model data, and dashed line shows field data; (b) Frequency spectrum of (a); (c) Time domain back-residual source computed using all frequency points; (d) Solid line is frequency spectrum of (c), and dashed line is frequency spectrum directly computed from formula (24). Left ‘y’ axis corresponds to amplitude of solid line and right ‘y’ axis corresponds to amplitude of dashed line; (e) Time domain back-residual source under the center frequency of 80 MHz; (f) Frequency spectrum of back-residual source under the center frequency of 80 MHz, 130 MHz, 180 MHz and by using formula (24), respectively.

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

(27)

(28)

(29)

(30)

這里需要注意的是，κε、κσ為不同的小穩(wěn)定因子，必須小心地為其選擇適當(dāng)?shù)闹挡⑶以诜囱葸^(guò)程中隨著迭代更新.

在本文的頻率域波形反演中，正演使用時(shí)間域有限差分(FDTD，O(2,4)).梯度的計(jì)算是在頻率域完成的，也就是說(shuō)每次正演之后都需要將時(shí)間域的電場(chǎng)值轉(zhuǎn)換到頻率域.我們使用離散傅里葉變換(DFT)來(lái)進(jìn)行時(shí)頻轉(zhuǎn)換，這樣無(wú)需保存所有時(shí)間的場(chǎng)值(能夠極大地減少內(nèi)存需求)，僅需在FDTD的每一次時(shí)間迭代中計(jì)算一次然后疊加即可：

(31)

整個(gè)反演過(guò)程如圖2所示，正演都在時(shí)間域進(jìn)行，而梯度及殘場(chǎng)源的計(jì)算在頻率域進(jìn)行.需要注意的是，后向傳播場(chǎng)的殘場(chǎng)源，是在頻率域計(jì)算后經(jīng)反傅里葉變換得到的.在本文中，為了避免高頻無(wú)用信號(hào)的干擾，提出只使用以當(dāng)前反演頻率為中心的一個(gè)窄頻帶內(nèi)的數(shù)據(jù)進(jìn)行時(shí)頻轉(zhuǎn)換.在計(jì)算完梯度，迭代步長(zhǎng)是在時(shí)間域計(jì)算的，之后對(duì)模型更新直至收斂完成.整個(gè)波形反演中各個(gè)步驟之間的計(jì)算順序是按照?qǐng)D2自上而下進(jìn)行的.

圖2　頻率域波形反演流程圖虛線代表數(shù)據(jù)的時(shí)頻域轉(zhuǎn)換，黑線代表波形反演流程.Fig.2　Flow chart of frequency domain waveform inversionDashed lines stand for Fourier transform and black lines represent waveform inversion process.

3算例

3.1復(fù)雜目標(biāo)體介電常數(shù)成像及頻率選擇策略

為了驗(yàn)證文中頻率域?qū)?shù)波形反演層析成像的能力，建立如圖3a所示的復(fù)雜模型.模型大小為6 m×6 m，包含3個(gè)地層：其中第1層與第3層地層參數(shù)相同，相對(duì)介電常數(shù)都為5；中間層的相對(duì)介電常數(shù)為5.5，并埋藏有兩根異常體管線；異常體位于深度3 m處，間隔1 m，相對(duì)介電常數(shù)為7；認(rèn)定所有地層與異常體的電導(dǎo)率相同且已知，均為0.001 S·m-1.在本次模擬中共13組發(fā)射(圓圈表示)，每組發(fā)射對(duì)應(yīng)13組接收(叉號(hào)表示).發(fā)射器與接收器位置均為垂直方向0 m到6 m之間，并以0.5 m等間隔排列.在本次反演中，我們認(rèn)為電導(dǎo)率是已知的，只對(duì)介電常數(shù)成像.初始的介電常數(shù)模型采用均勻地層，相對(duì)介電常數(shù)值為5.5，與中間層相同.

圖3　不同頻率策略下的相對(duì)介電常數(shù)反演結(jié)果(a) 真實(shí)模型； (b) fS=0 MHz； (c) fS=10 MHz； (d) fS=20 MHz； (e) fS=30 MHz； (f) fS=40 MHz.Fig.3　Relative permittivity results of different frequency strategies(a) Real model; (b) fS=0 MHz; (c) fS=10 MHz; (d) fS=20 MHz; (e) fS=30 MHz; (f) fS=40 MHz.

反演的起始頻率為50 MHz，每經(jīng)過(guò)一次反演則反演頻率增加2 MHz.為了加速收斂，每十次反演之后，反演頻率進(jìn)行一次跳躍.例如：當(dāng)前反演次數(shù)為n，則反演頻率為50+fS·int(n/10)+2(n-1) MHz，其中int表示取正整數(shù)，fS為跳躍的頻率寬度.為了驗(yàn)證該頻率策略，以圖3a為真實(shí)模型進(jìn)行反演：圖3b為無(wú)頻率跳躍時(shí)的對(duì)數(shù)波形反演結(jié)果；圖3c、3d、3e、3f分別為頻率跳躍寬度為10 MHz、20 MHz、30 MHz、40 MHz時(shí)的對(duì)數(shù)波形反演結(jié)果.為了使不同結(jié)果之間更具對(duì)比性，在每個(gè)頻點(diǎn)處只反演一次，且最終的反演頻率都為228 MHz.

從結(jié)果來(lái)看，5種不同模式下的對(duì)數(shù)波形反演結(jié)果都能清晰地重建上下兩層界面及中間層中埋藏的異常體管線，但圖3f中兩個(gè)異常體管線的成像效果較弱.圖4a為5種模式波形反演對(duì)應(yīng)的目標(biāo)函數(shù)曲線：青綠色曲線為無(wú)頻率跳躍時(shí)的目標(biāo)函數(shù)曲線，紅色、綠色、藍(lán)色和洋紅色曲線分別對(duì)應(yīng)跳躍頻率寬度為10 MHz、20 MHz、30 MHz和40 MHz時(shí)的目標(biāo)函數(shù)曲線.由于前10次反演的頻率完全相同，故在圖4a中所有曲線在10次之前重合.無(wú)頻率跳躍的曲線收斂較慢，這是由于該模式下低頻部分的反演次數(shù)較多造成的.圖4a同時(shí)也反映了不同模式的反演效率：無(wú)頻率跳躍模式下共反演90次，頻率跳躍寬度為10 MHz時(shí)共反演60次，隨著頻率跳躍寬度的增加反演次數(shù)分別為50、45和30.圖4b為沿AA′方向的介電常數(shù)切線，黑線為真實(shí)模型，可以看出只有當(dāng)頻率跳躍寬度為40 MHz時(shí)反演效果出現(xiàn)明顯的降低，而其他結(jié)果無(wú)較大差別.

圖3證明在頻率域波形反演中，適當(dāng)?shù)念l率跳躍能夠有效地加速反演過(guò)程的收斂且不影響效果.需要注意的是：隨著頻率跳躍寬度的增大，計(jì)算效率的提升是遞減的；如果頻率跳躍寬度過(guò)大，則反演頻率的不足會(huì)影響反演結(jié)果(圖3f).因此，要在計(jì)算效率和反演效果之間取舍，選擇恰當(dāng)?shù)念l率跳躍寬度.另外，頻率跳躍的可行性建立在第一次頻率跳躍之前反演過(guò)程已經(jīng)取得較好收斂效果的基礎(chǔ)上.這隱含了一個(gè)條件，即：在頻率跳躍之前，目標(biāo)函數(shù)已經(jīng)收斂到全局最小附近.在本文之后的算例中，選擇較為保守的10 MHz為頻率跳躍寬度.

3.2介電常數(shù)和電導(dǎo)率同時(shí)成像

以上算例證明了基于對(duì)數(shù)目標(biāo)函數(shù)的頻率域波形反演在單介電常數(shù)成像(認(rèn)定電導(dǎo)率均勻且已知)的能力和文中反演策略的有效性.在實(shí)際情況中，介電常數(shù)和電導(dǎo)率是同時(shí)變化的，本例對(duì)介電常數(shù)和電導(dǎo)率同時(shí)反演并與基于傳統(tǒng)目標(biāo)函數(shù)的頻率域波形反演結(jié)果對(duì)比.

真實(shí)模型如圖5a和圖5d所示，發(fā)射天線(用圓圈表示)位于水平距離0 m，深度從0 m到20 m，共41個(gè)發(fā)射，間距0.5 m；接收天線(用叉號(hào)表示)位于水平距離10 m處，深度同樣從0 m到20 m，間距0.5 m共41個(gè).三根直徑為0.5 m的管線異常體位于深度6.5 m處，互相之間間隔2 m，相對(duì)介電常數(shù)為7，電導(dǎo)率為0.008 S·m-1.另一個(gè)較大的管線異常體位于深度12 m，直徑為1 m，相對(duì)介電常數(shù)為6.5，電導(dǎo)率為0.005 S·m-1.整個(gè)模型分為三層，所有的異常體都埋藏在中間層，頂層和底層的電性參數(shù)相同(吳俊軍等，2014).將反演的起始頻率定在30 MHz，采用與之前算例相同的頻率策略.圖5b、5e、5c、5f為迭代60次之后的反演結(jié)果.

圖5c和圖5f分別為對(duì)數(shù)目標(biāo)函數(shù)下介電常數(shù)和電導(dǎo)率反演結(jié)果，圖5b和圖5e分別為傳統(tǒng)目標(biāo)函數(shù)下介電常數(shù)和電導(dǎo)率結(jié)果.能夠很明顯地看出對(duì)數(shù)目標(biāo)函數(shù)的結(jié)果優(yōu)于傳統(tǒng)目標(biāo)函數(shù)的結(jié)果：深度6.5 m處的三根管線異常體，在圖5c和圖5f中表現(xiàn)得十分清晰，且互相之間分隔很好；而圖5b和圖5e中異常連到一起，無(wú)法辨認(rèn)出異常體的數(shù)量和形狀.對(duì)于深度12 m處的大異常體，傳統(tǒng)目標(biāo)函數(shù)的結(jié)果表現(xiàn)出明顯的畸變(水平方向拉長(zhǎng))，而對(duì)數(shù)目標(biāo)函數(shù)的結(jié)果與真實(shí)模型十分接近，表現(xiàn)為正圓而非橢圓.

圖6a為沿AA′方向的介電常數(shù)切線，圖6b為沿BB′方向的電導(dǎo)率切線.對(duì)比真實(shí)模型、對(duì)數(shù)目標(biāo)函數(shù)結(jié)果和傳統(tǒng)目標(biāo)函數(shù)結(jié)果在相同位置切線的形態(tài)和數(shù)值，對(duì)數(shù)目標(biāo)函數(shù)結(jié)果更接近真實(shí)模型，也更能表現(xiàn)出異常體的形態(tài).從圖6b發(fā)現(xiàn)，對(duì)數(shù)目標(biāo)函數(shù)結(jié)果與真實(shí)模型在數(shù)值上的一致性相當(dāng)好，這說(shuō)明了在同時(shí)反演的情況下，電導(dǎo)率成像效果好于介電常數(shù)成像(吳俊軍等，2014).在傳統(tǒng)的射線理論反演中，介電常數(shù)反演要比電導(dǎo)率反演容易，這是因?yàn)樽邥r(shí)成像在理論上比衰減成像完善.而在波形反演中，電導(dǎo)率成像效果好于介電常數(shù)成像效果是由于波形中振幅信息對(duì)電導(dǎo)率的變化更加敏感(相對(duì)于相位信息對(duì)介電常數(shù)的變化).

圖7為第21個(gè)源(Depth：10 m)發(fā)射，60次迭代后的模型的模擬接收波形與實(shí)際數(shù)據(jù)波形的對(duì)比(接收器波形為41道，為顯示方便，每隔四道選擇一道).盡管圖5中傳統(tǒng)目標(biāo)函數(shù)波形反演成像結(jié)果與真實(shí)模型差別較大，但在圖7a中，藍(lán)線代表的反演數(shù)據(jù)波形與實(shí)際數(shù)據(jù)波形仍符合很好，這說(shuō)明傳統(tǒng)目標(biāo)函數(shù)波形反演結(jié)果是收斂的.圖7c中的藍(lán)線和紅線分別為圖7a中的藍(lán)線與黑線的差和圖7b中紅線與黑線的差.在圖7c第25道至41道，能夠明顯發(fā)現(xiàn)藍(lán)線有異常，說(shuō)明圖7a中傳統(tǒng)目標(biāo)函數(shù)結(jié)果與真實(shí)模型之間的差要大于圖7b中對(duì)數(shù)目標(biāo)函數(shù)結(jié)果與真實(shí)模型的差.這也證明了對(duì)數(shù)目標(biāo)函數(shù)波形反演優(yōu)于傳統(tǒng)目標(biāo)函數(shù)波形反演.

圖4　(a)目標(biāo)函數(shù)收斂曲線； (b) 深度3 m介電常數(shù)切線(沿AA′方向)Fig.4　(a) Objective function convergence curves; (b) Permittivity in the depth of 3 m (along AA′ line)

圖5　介電常數(shù)和電導(dǎo)率同時(shí)成像(a)介電常數(shù)真實(shí)模型； (b) 傳統(tǒng)目標(biāo)函數(shù)下介電常數(shù)反演結(jié)果； (c) 對(duì)數(shù)目標(biāo)函數(shù)下介電常數(shù)反演結(jié)果； (d) 電導(dǎo)率真實(shí)模型； (e) 傳統(tǒng)目標(biāo)函數(shù)下電導(dǎo)率反演結(jié)果； (f) 對(duì)數(shù)目標(biāo)函數(shù)下電導(dǎo)率反演結(jié)果.Fig.5　Simultaneous tomography of permittivity and conductivity(a) Real model of permittivity; (b) Permittivity results of conventional object function; (c) Permittivity results of logarithmic object function; (d) Real model of conductivity; (e) Conductivity results of conventional object function; (f) Conductivity results of logarithmic object function.

圖6　(a) 介電常數(shù)AA′切線； (b) 電導(dǎo)率BB′切線Fig.6　(a) Permittivity along AA′ line; (b) Conductivity along BB′ line

圖7　第21個(gè)源(depth：10 m)發(fā)射，接收器接收到的波形對(duì)比(a)60次迭代后傳統(tǒng)目標(biāo)函數(shù)波形反演數(shù)據(jù)(藍(lán)線)與實(shí)際數(shù)據(jù)(黑線)對(duì)比； (b) 60次迭代后對(duì)數(shù)目標(biāo)函數(shù)波形反演數(shù)據(jù)(紅線)與實(shí)際數(shù)據(jù)(黑線)對(duì)比； (c) 藍(lán)線為(a)中藍(lán)線與黑線的差，紅線為(b)中紅線與黑線的差.Fig.7　Transmitter gathers generated by the source placed at a 10 m depth in Fig.5(a) Blue and black lines express the differences of radar trace generated from the results of conventional waveform inversion and field data after 60 iterations; (b) Red and black lines show the differences of radar trace generated from the results of logarithmic waveform inversion and field data after 60 iterations; (c) Blue line represents the difference between blue and black lines in (a), red line shows the difference of red and black lines in (b).

3.3實(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í)成像(王飛等，2013)和重心頻率下降法(Liu et al., 1998)得到.

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

圖8a和圖8c分別為基于走時(shí)反演結(jié)果計(jì)算的相對(duì)介電常數(shù)和重心頻率下降法得到的電導(dǎo)率分布圖像，圖8b和圖8d為對(duì)數(shù)目標(biāo)函數(shù)波形反演得到的相應(yīng)結(jié)果.圖9給出了第10個(gè)源(Depth：19 m)發(fā)射，基于射線結(jié)果正演和波形反演結(jié)果的正演波形對(duì)比.對(duì)比圖8a和圖8b發(fā)現(xiàn)，對(duì)數(shù)波形反演的相對(duì)介電常數(shù)結(jié)果相對(duì)于走時(shí)反演結(jié)果變化不大，這主要是因?yàn)樽邥r(shí)反演得到的相對(duì)介電常數(shù)結(jié)果已經(jīng)較為準(zhǔn)確.但需要注意的是，對(duì)比圖9a和圖9b發(fā)現(xiàn)，波形反演得到的波形相對(duì)于基于射線結(jié)果正演得到的波形，與實(shí)際數(shù)據(jù)波形之間的相位符合更好.這說(shuō)明，波形反演在走時(shí)反演的基礎(chǔ)之上仍有提高，從而得到更加準(zhǔn)確的相對(duì)介電常數(shù)分布.對(duì)比圖8c和圖8d，波形反演得到的電導(dǎo)率結(jié)果較重心頻率下降法得到的電導(dǎo)率結(jié)果有明顯的提高.無(wú)論是波形反演結(jié)果還是射線結(jié)果，都能明顯觀察到深度28 m位置處的地下異常體，而波形反演得到的電導(dǎo)率結(jié)果呈現(xiàn)了更多的細(xì)節(jié).另外，波形反演得到的相對(duì)介電常數(shù)在異常體內(nèi)部表現(xiàn)出了一個(gè)條帶狀的較低介電常數(shù)區(qū)域，與波形反演得到的電導(dǎo)率結(jié)果在相同位置處具有同樣的特征.由于缺乏對(duì)應(yīng)的地質(zhì)信息，無(wú)法判斷異常體的具體類(lèi)型.需要注意的是，盡管對(duì)實(shí)際數(shù)據(jù)做了三維校正，但反演結(jié)果仍受一些因素的影響(如反演的多解性，噪聲以及天線等)，因此反演結(jié)果中的一些細(xì)節(jié)變化并不一定可靠，需通過(guò)波形對(duì)比來(lái)判斷反演結(jié)果整體的好壞.從圖9c也能直觀地看出，波形反演的結(jié)果無(wú)論是在相位還是振幅上都比射線結(jié)果更加貼近實(shí)際數(shù)據(jù).因此，波形反演得到的結(jié)果更加準(zhǔn)確，而且相對(duì)介電常數(shù)和電導(dǎo)率之間的符合較好.因此，對(duì)數(shù)波形反演方法在實(shí)際數(shù)據(jù)處理中的應(yīng)用結(jié)果表明，該方法優(yōu)于傳統(tǒng)的射線類(lèi)反演方法.

4總結(jié)與討論

本文詳細(xì)推導(dǎo)了對(duì)數(shù)目標(biāo)函數(shù)下頻率域探地雷達(dá)波形反演層析成像方法，與地震波形反演通過(guò)對(duì)目標(biāo)函數(shù)和波動(dòng)方程求導(dǎo)來(lái)得到梯度表達(dá)式不同，本文借鑒了Meles等的擾動(dòng)理論推導(dǎo)梯度，并指出頻率域梯度表達(dá)式中Re的意義是時(shí)間域波形反演梯度公式中體現(xiàn)不出的.使用離散傅里葉變換(DFT)將每次正演的時(shí)間域數(shù)據(jù)轉(zhuǎn)換到頻率域的方法只增加極少量的計(jì)算量，避免了巨量時(shí)間域數(shù)據(jù)的存儲(chǔ)，有效減少計(jì)算機(jī)內(nèi)存需求.在計(jì)算殘場(chǎng)源時(shí)，為了避免高頻無(wú)用信號(hào)的干擾，提出只使用以當(dāng)前反演頻率為中心的一個(gè)窄頻帶數(shù)據(jù)即可.在頻率迭代策略上，文中采用串行頻率策略，并且為了加速收斂每迭代10次則反演頻率增加10 MHz.不同于地震波形反演每個(gè)頻點(diǎn)都反演多次，本文每個(gè)頻點(diǎn)只進(jìn)行一次反演，這有兩點(diǎn)原因：一是文中每次反演頻率只增加很少(2 MHz)，反映在波長(zhǎng)上變化微弱；二是每次反演過(guò)程中都計(jì)算迭代步長(zhǎng)，保證反演更有效地收斂.文中三個(gè)算例也證明本文的頻率策略是可行的.

圖8　貴州實(shí)際數(shù)據(jù)反演結(jié)果(a) 走時(shí)反演得到的相對(duì)介電常數(shù)； (b) 對(duì)數(shù)波形反演得到的相對(duì)介電常數(shù)； (c) 重心頻率下降法得到的電導(dǎo)率； (d) 對(duì)數(shù)波形反演得到的電導(dǎo)率Fig.8　Inversion results of field data from Guizhou(a) Permittivity results from first-arrival times method; (b) Permittivity results of logarithmic waveform inversion; (c) Conductivity results from applying centroid frequency downshift method; (d) Conductivity results of logarithmic waveform inversion.

圖9　貴州實(shí)際數(shù)據(jù)反演結(jié)果的波形對(duì)比(a) 黑線和藍(lán)線分別為實(shí)際數(shù)據(jù)波形和基于射線模型的正演波形； (b) 黑線和紅線分別為實(shí)際數(shù)據(jù)波形和最終全波形反演結(jié)果對(duì)應(yīng)的波形； (c) 藍(lán)線和紅線分別為(a)中黑線與藍(lán)線的差和(b)中黑線和紅線的差.所有數(shù)據(jù)都根據(jù)實(shí)際數(shù)據(jù)的最大振幅歸一化.Fig.9　Transmitter gathers generated by the source placed at a 19 m depth in Fig.8(a) Blue and black lines show radar traces generated from results of ray tomograms and field data; (b) Red and black lines show radar traces generated from results of logarithmic waveform inversion and field data; (c) Blue and red lines show difference between blue and black lines in (a) and between red and black lines in (b), respectively. Amplitudes in all panels are normalized with respect to maximum amplitude of field data.

合成數(shù)據(jù)的反演結(jié)果證明對(duì)數(shù)目標(biāo)函數(shù)下波形反演優(yōu)于傳統(tǒng)目標(biāo)函數(shù)下波形反演，原因是使用對(duì)數(shù)目標(biāo)函數(shù)時(shí)，波場(chǎng)殘差受正演波場(chǎng)的自然比例調(diào)節(jié)作用.盡管如此，由于在推導(dǎo)過(guò)程中相位存在假設(shè)，使用對(duì)數(shù)目標(biāo)函數(shù)的探地雷達(dá)波形反演對(duì)介電常數(shù)初始模型的要求較高.在實(shí)際數(shù)據(jù)處理中，數(shù)據(jù)的三維校正和源子波估計(jì)仍是重點(diǎn)，另外需要盡可能準(zhǔn)確地選擇初始模型來(lái)保證反演過(guò)程收斂.

未來(lái)我們需要關(guān)注探地雷達(dá)三維波形反演，其中FDTD仍將發(fā)揮重要作用.另外初始模型的獲取可用Laplace域波形反演及Laplace-Fourier域波形反演替代傳統(tǒng)的射線理論反演，從而獲得更可靠的初始模型.在頻率域探地雷達(dá)波形反演中波場(chǎng)正演方法、頻率選擇策略、目標(biāo)函數(shù)設(shè)置方式、源子波估計(jì)等問(wèn)題依然需要進(jìn)一步的研究.

致謝感謝康涅狄格大學(xué)土木與環(huán)境工程學(xué)院劉瀾波老師提供重心頻率下降法的程序.

References

Cai W Y, Qin F H, Schuster G T. 1996. Electromagnetic velocity inversion using 2-D Maxwell′s equations.Geophysics, 61(4): 1007-1021.

Chung W, Ha T, Ha W, et al. 2007. Robust seismic waveform inversion using back-propagation algorithm.∥ 77th Annual International Meeting, Expanded Abstracts: 1780-1784.

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

Fu L, Liu S X, Liu L B, et al. 2014. Airborne ground penetrating radar numerical simulation and reverse time migration.ChineseJ.Geophys. (in Chinese), 57(5): 1636-1646, doi: 10.6038/cjg20140526.

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-1179. 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.IEEEGeoscienceandRemoteSensingLetters, 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.

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 R G. 1990. Inverse theory applied to multi-source cross-hole tomography, Part 2: Elastic wave-equation method.GeophysicalProspecting, 38(3): 311-329. Pratt R G, Worthington M H. 1990. Inverse theory applied to multi-source cross-hole tomography, Part 1: Acoustic wave-equation method.GeophysicalProspecting, 38(3): 287-310.

Pratt G, Shin C, Hicks G J. 1998. Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion.GeophysicalJournalInternational, 133(2): 341-362.

Pratt R G. 1999. Seismic waveform inversion in the frequency domain, Part 1: Theory and verification in a physical scale

model.Geophysics, 64(3): 888-901.

Pratt R G, Shipp R M. 1999. Seismic waveform inversion in the frequency domain, Part 2: Fault delineation in sediments using crosshole data.Geophysics, 64(3): 902-914.

Shin C, Min D J. 2006. Waveform inversion using a logarithmic wavefield.Geophysics, 71(3): R31-R42.

Sirgue L, Pratt R G. 2004. Efficient waveform inversion and imaging: A strategy for selecting temporal frequencies.Geophysics, 69(1): 231-248. Song Z M, Willianson P R, Pratt R G. 1995. Frequency-domain acoustic-wave modeling and inversion of crosshole data: Part Ⅱ: Inversion method, synthetic experiments and real-data results.Geophysics, 60(3): 796-809.

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

Tarantola A. 1986. A strategy for nonlinear elastic inversion of seismic reflection data.Geophysics, 51(10): 1893-1903.

Wang F, Liu S X, Qu X X, et al. 2013. Crosshole radar traveltime tomography without ray tracing using the high accuracy fast marching method.ChineseJ.Geophys. (in Chinese), 56(11): 3896-3907, doi: 10.6038/cjg20131131.

Wang Z L, Zhou H, Li G F. 2007. Inversion of ground-penetrating radar data for 2D electric parameters.ChineseJ.Geophys. (in Chinese), 50(3): 897-904.

Williamson P R, Worthington M H. 1993. Resolution limits in ray tomography due to wave behavior: Numerical experiments.Geophysics, 58(5): 727-735.

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. 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. Zhou H, Sato M, Takenaka T, et al. 2007. Reconstruction from antenna-transformed radar data using a time-domain reconstruction method.IEEETransactionsonGeoscienceandRemoteSensing, 45(3): 689-696.

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

傅磊, 劉四新, 劉瀾波等. 2014. 機(jī)載探地雷達(dá)數(shù)值模擬及逆時(shí)偏移成像. 地球物理學(xué)報(bào), 57(5): 1636-1646, doi: 10.6038/cjg20140526.

王飛, 劉四新, 曲昕馨等. 2013. 基于HAFMM的無(wú)射線追蹤跨孔雷達(dá)走時(shí)層析成像. 地球物理學(xué)報(bào), 56(11): 3896-3907, doi: 10.6038/cjg20131131.

王兆磊, 周輝, 李國(guó)發(fā). 2007. 用地質(zhì)雷達(dá)數(shù)據(jù)資料反演二維地下介質(zhì)的方法. 地球物理學(xué)報(bào), 50(3): 897-904.

吳俊軍, 劉四新, 李彥鵬等. 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.

(本文編輯何燕)

Frequency domain waveform inversion of cross-hole GPR data based on a logarithmic objective function

MENG Xu, LIU Si-Xin*，F(xiàn)U Lei, WANG Xian-Nan, LIU Xin-Tong, WANG Wen-Tian, CAI Jia-Qi

CollegeofGeoExplorationScienceandTechnology,JilinUniversity,Changchun130026,China

AbstractGround penetrating radar (GPR) tomography plays an important role in geology, hydrogeology and engineering investigations. Traditional tomography (i.e., the first-arrival times and maximum first-cycle amplitudes) based on ray theory cannot provide high-resolution images because only a fraction of the information contained in the radar data is used in the inversion. In recent years, waveform inversion is one of the biggest concerns because it can provide sub-wavelength resolution images. Waveform inversion has been applied in GPR over ten years, but most of the results are computed in the time domain. In the frequency domain, the choice of inverted frequencies is flexible and different types of objective functions can be used, therefore the results of frequency-domain waveform inversion are more diversified than that of the time domain.

KeywordsWaveform inversion; Logarithmic objective function; Permittivity; Conductivity; Cross-hole radar

基金項(xiàng)目國(guó)家高技術(shù)研究發(fā)展計(jì)劃(2013AA064603),國(guó)家自然科學(xué)基金(40874073，41074076)資助.

作者簡(jiǎn)介孟旭，男，1987年生，博士研究生，主要從事探地雷達(dá)數(shù)據(jù)處理與解釋、電磁波數(shù)值模擬及全波形反演的研究.E-mail：mengxu519@126.com*通訊作者劉四新，男，1966年生，山西太谷人，日本東北大學(xué)工學(xué)博士，博士生導(dǎo)師，主要從事探地雷達(dá)、鉆孔雷達(dá)及電磁波測(cè)井等的方法理論和應(yīng)用方面的研究.E-mail：liusixin@jlu.edu.cn

doi:10.6038/cjg20160530 中圖分類(lèi)號(hào)P631

收稿日期2015-05-26，2016-03-04收修定稿

孟旭， 劉四新， 傅磊等.2016. 基于對(duì)數(shù)目標(biāo)函數(shù)的跨孔雷達(dá)頻域波形反演.地球物理學(xué)報(bào)，59(5):1875-1887,doi:10.6038/cjg20160530.

Meng X, Liu S X, Fu L, et al. 2016. Frequency domain waveform inversion of cross-hole GPR data based on a logarithmic objective function.ChineseJ.Geophys. (in Chinese),59(5):1875-1887,doi:10.6038/cjg20160530.