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

吉林大學地球探測科學與技術學院， 長春　130026

基于對數(shù)目標函數(shù)的跨孔雷達頻域波形反演

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

吉林大學地球探測科學與技術學院， 長春130026

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

關鍵詞波形反演； 對數(shù)目標函數(shù)； 介電常數(shù)； 電導率； 鉆孔雷達

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引言

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

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

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

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

2理論

2.1正演問題的頻域表達形式

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

(1)

(2)

(3)

(4)

2.2對數(shù)目標函數(shù)及梯度計算

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

(5)

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

(7)

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

(8)

(9)

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

(10)

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

(11)

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

(12)

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

(13)

(14)

在空間域，這是一個關于擾動正演電場的方程.將(14)式代入(3)式中：

(15)

(16)

(17)

則

(18)

目標函數(shù)梯度的推導使用擾動目標函數(shù)的一階近似：

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

(19)

擾動目標函數(shù)為：

(21)

其中，Rn(ε,σ)代表高階項.將(21)式代入(20)式中展開化簡得到：

(22)

(23)其中

(24)

(25)

(26)

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

2.3后向傳播殘場源

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

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

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

圖1　后向傳播殘場源的變換(a)實線為正演數(shù)據(jù)，虛線為實際數(shù)據(jù)； (b) 實線為正演數(shù)據(jù)的頻譜，虛線為實際數(shù)據(jù)的頻譜； (c) 使用所有頻點計算得到的后向殘場源； (d) 實線為(c)對應的頻譜，虛線為直接由式(24)計算的頻譜.左側(cè)y軸對應實線振幅，右側(cè)y軸為虛線振幅； (e) 以80 MHz為中心的窄帶數(shù)據(jù)計算的時間域后向殘場源； (f) 中心頻率為80 MHz、130 MHz、180 MHz及直接使用公式(24)計算的后向殘場源頻譜.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迭代步長的求解

(27)

(28)

(29)

(30)

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

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

(31)

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

圖2　頻率域波形反演流程圖虛線代表數(shù)據(jù)的時頻域轉(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復雜目標體介電常數(shù)成像及頻率選擇策略

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

圖3　不同頻率策略下的相對介電常數(shù)反演結(jié)果(a) 真實模型； (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)過一次反演則反演頻率增加2 MHz.為了加速收斂，每十次反演之后，反演頻率進行一次跳躍.例如：當前反演次數(shù)為n，則反演頻率為50+fS·int(n/10)+2(n-1) MHz，其中int表示取正整數(shù)，fS為跳躍的頻率寬度.為了驗證該頻率策略，以圖3a為真實模型進行反演：圖3b為無頻率跳躍時的對數(shù)波形反演結(jié)果；圖3c、3d、3e、3f分別為頻率跳躍寬度為10 MHz、20 MHz、30 MHz、40 MHz時的對數(shù)波形反演結(jié)果.為了使不同結(jié)果之間更具對比性，在每個頻點處只反演一次，且最終的反演頻率都為228 MHz.

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

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

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

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

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

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

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

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

圖4　(a)目標函數(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ù)和電導率同時成像(a)介電常數(shù)真實模型； (b) 傳統(tǒng)目標函數(shù)下介電常數(shù)反演結(jié)果； (c) 對數(shù)目標函數(shù)下介電常數(shù)反演結(jié)果； (d) 電導率真實模型； (e) 傳統(tǒng)目標函數(shù)下電導率反演結(jié)果； (f) 對數(shù)目標函數(shù)下電導率反演結(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) 電導率BB′切線Fig.6　(a) Permittivity along AA′ line; (b) Conductivity along BB′ line

圖7　第21個源(depth：10 m)發(fā)射，接收器接收到的波形對比(a)60次迭代后傳統(tǒng)目標函數(shù)波形反演數(shù)據(jù)(藍線)與實際數(shù)據(jù)(黑線)對比； (b) 60次迭代后對數(shù)目標函數(shù)波形反演數(shù)據(jù)(紅線)與實際數(shù)據(jù)(黑線)對比； (c) 藍線為(a)中藍線與黑線的差，紅線為(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ù)據(jù)

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

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

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

4總結(jié)與討論

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

圖8　貴州實際數(shù)據(jù)反演結(jié)果(a) 走時反演得到的相對介電常數(shù)； (b) 對數(shù)波形反演得到的相對介電常數(shù)； (c) 重心頻率下降法得到的電導率； (d) 對數(shù)波形反演得到的電導率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ù)據(jù)反演結(jié)果的波形對比(a) 黑線和藍線分別為實際數(shù)據(jù)波形和基于射線模型的正演波形； (b) 黑線和紅線分別為實際數(shù)據(jù)波形和最終全波形反演結(jié)果對應的波形； (c) 藍線和紅線分別為(a)中黑線與藍線的差和(b)中黑線和紅線的差.所有數(shù)據(jù)都根據(jù)實際數(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é)果證明對數(shù)目標函數(shù)下波形反演優(yōu)于傳統(tǒng)目標函數(shù)下波形反演，原因是使用對數(shù)目標函數(shù)時，波場殘差受正演波場的自然比例調(diào)節(jié)作用.盡管如此，由于在推導過程中相位存在假設，使用對數(shù)目標函數(shù)的探地雷達波形反演對介電常數(shù)初始模型的要求較高.在實際數(shù)據(jù)處理中，數(shù)據(jù)的三維校正和源子波估計仍是重點，另外需要盡可能準確地選擇初始模型來保證反演過程收斂.

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

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

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.

附中文參考文獻

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

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

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

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

周輝, 陳漢明, 李卿卿等. 2014. 不需提取激發(fā)脈沖的探地雷達波形反演方法. 地球物理學報, 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

基金項目國家高技術研究發(fā)展計劃(2013AA064603),國家自然科學基金(40874073，41074076)資助.

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

doi:10.6038/cjg20160530 中圖分類號P631

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

孟旭， 劉四新， 傅磊等.2016. 基于對數(shù)目標函數(shù)的跨孔雷達頻域波形反演.地球物理學報，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.