吳俊軍，劉四新，李彥鵬，張宇生，張麗麗，傅磊，王飛，孟旭，雷林林

1中國石油集團東方地球物理勘探有限責(zé)任公司新興物探開發(fā)處，涿州 072751

2吉林大學(xué) 地球探測科學(xué)與技術(shù)學(xué)院，長春 130026

3沈陽航空航天大學(xué) 電子信息工程學(xué)院，沈陽 110034

1 引言

跨孔雷達全波形反演是一種使用全波形信息反演兩鉆孔之間地下信息的層析成像技術(shù).在地震勘探中，全波形反演能夠提供小于波長的分辨率，在理想狀態(tài)下，分辨率能夠達到1／2～1／3波長（Pratt，1999）.利用聲波的全波形反演方法比利用射線方法在分辨率上能提高一個量級.探地雷達全波形反演的發(fā)展得益于其與地震學(xué)的高度相似性，以及近年來電磁學(xué)的快速發(fā)展.例如將麥克斯韋方程代替聲波方程或彈性波方程，而解的形式與地震學(xué)非常相似，所以地震全波反演方法在過去三十年里所取得的發(fā)展，對跨孔雷達全波反演方法的研究起到了很大的促進作用.Johnson等（2005）以及 Cai等（1996）分別使用了菲涅爾量以及程函方程的方法進一步推動了波形反演技術(shù)的發(fā)展.許多探地雷達學(xué)者采用了修正的散射成像技術(shù)進行波形成像（Zhou and Liu，2000；Zhou et al.，2007；Cui et al.，2004）.Kuroda等（2005）和 Meles等（Meles et al.，2010；Meles et al.，2011；Meles et al.，2012）使用全波反演的方法對跨孔雷達數(shù)據(jù)成像，并取得較好的效果.Kuroda等的方法只考慮介電常數(shù)，而Ernst等（Ernst et al.，2007）通過級聯(lián)更新的方法既考慮介電常數(shù)，又考慮電導(dǎo)率.后者的基本想法是固定電導(dǎo)率反演介電常數(shù)，反之亦然.兩種方法的共同點在于其本質(zhì)上都是標(biāo)量的.Meles等提出的波形反演方法考慮電參數(shù)的矢量性，反演過程中允許介電常數(shù)和電導(dǎo)率同時更新.與常規(guī)射線追蹤方法相比，全波反演方法能夠更加有效的識別亞波長異常體.Belina等對不同隨機模型下的全波反演子波估計進行了深入的研究（Belina et al.，2012a；Belina et al.，2012b；Belina et al.，2012c）.

本文全面推導(dǎo)了全波形跨孔雷達層析成像反演方法.通過基于局域網(wǎng)的分布式并行算法，有效地解決了巨型數(shù)據(jù)正演計算問題，從而實現(xiàn)了全波形跨孔層析成像方法在普通計算機上全面有效的運行.本文中首先建立了基于單軸各向異性介質(zhì)完全匹配層（UPML）吸收邊界的TE模式下的時間域有限差分 （FDTD）二 維 正 演 算 法 （Taflove Hagness，2005）.由于跨孔雷達測量方式中發(fā)射天線與接收天線極化方式為垂直于地表的z方向，因此選用FDTD正演算法時需要使用有Ex、Ez、及Hy的TE模式.反演算法采用最速梯度法求解，通過應(yīng)用包括時間維度在內(nèi)的全波場信息與殘場逆向傳播的全波場信息乘積來計算梯度方向，通過求取以步長為自變量的目標(biāo)函數(shù)的極值確定步長公式.由于介電常數(shù)與電導(dǎo)率在量級上的區(qū)別，因此在確定迭代步長時需要對兩種電性參數(shù)進行分別計算以提高收斂速度，這里需要分別使用不同的穩(wěn)定因子.文中對介電常數(shù)及電導(dǎo)率分別采用對數(shù)化處理，并推導(dǎo)出新的梯度公式.采用參數(shù)對數(shù)化處理能夠很好的提高收斂速度及穩(wěn)定性，并且拓寬了異常與周邊介質(zhì)物性參數(shù)對比度.

在不降低分辨率的情況下，通過對反演模型采用抽稀的辦法，將數(shù)據(jù)減小至原內(nèi)存需求的1／6，以解決大型模型對計算機內(nèi)存的巨大需求，文中分別對不同類型的模型進行了全波反演成像，取得了較好的結(jié)果.本文提出將第一步單介電常數(shù)反演作為初始模型（吳俊軍，2012），有效的解決了同步反演收斂速度較慢的問題.

2 跨孔雷達全波反演理論推導(dǎo)

2.1 正演問題的積分表達形式

為了更好地全面表達波場情況，麥克斯韋方程可以表示成如下形式（Meles et al.，2010）：

其中，Js為電流密度源，Es和Hs為產(chǎn)生的電場和磁場，上角標(biāo)s表示指定的源，M（ε，σ）為一線性算子，在本文中由于未涉及到磁介質(zhì)，因此本文中假定磁導(dǎo)率為常數(shù)，其值為μ0.M（ε，σ）代表麥克斯韋方程組，使得在任意點任意時刻方程組可以表示如下：

其中，ε（x）、σ（x）分別為空間變量x的介電常數(shù)和電導(dǎo)率分布.場量值是關(guān)于空間和時間的矢量，而電性參數(shù)僅是位置變量x的標(biāo)量.向量Es和Hs不是特指某一時刻某一空間點的值，而是包括整個空間域和時間域的所有場值，即

由于在以下的討論中均只使用了電場值，因此省略Hs，得到一個簡單的方程：

式中，Es（x，t）為源Js在 （x，t）產(chǎn)生的電場，G（x，t，x′，t′）為源Js（x′，t′）的格林函數(shù)張量.對于三維介質(zhì)，可以用下角標(biāo)的方式寫出每個方向分量的電場：

其中，Esi（x，t）為電場各個方向的分量，Gik（x，t，x′，t′）為格林函數(shù)張量各個方向的元素，Jsk（x′，t′）為源對應(yīng)各個方向的分量.本文中使用向量的方式來表達方程，在使用中很容易就能轉(zhuǎn)變?yōu)閺埩啃问剑∕eles et al.，2010）.

2.2 目標(biāo)函數(shù)方程的梯度計算

地球物理反演中有多種多樣的反演問題（Greenhalgh et al.，2006）.全波形層析成像反演問題主要目的是尋找介電常數(shù)ε與電導(dǎo)率σ的空間分布.利用正演模擬數(shù)據(jù)道與實測數(shù)據(jù)道之差的平方最小準(zhǔn)則來求解目標(biāo)函數(shù)S（ε，σ）.目標(biāo)函數(shù)可以寫為

式中，Es（ε，σ）和分別為正演數(shù)據(jù)和觀測數(shù)據(jù)，T為轉(zhuǎn)置符號.（7）式分別在源s，接收點d以及觀測時間τ三個方向上求和.誤差方程的關(guān)鍵部分是電場求解，由于電場是關(guān)于介電常數(shù)和電導(dǎo)率的方程.要找到誤差方程的梯度，需要對關(guān)于這些參數(shù)的電場進行一階近似.與（1）式相似，可以寫出擾動系統(tǒng)下的麥克斯韋方程組.

（8）式減去（1）式得

（9）式包含了與（1）式相同的算子M，此時作為源的這一項表示如下：

在空間域，這是一個關(guān)于擾動正演電場的方程.將（10）式代入到（4）式得到：

將（11）式展開得到：

引入一個新核函數(shù)（或稱敏感因子）Ls可以分別表示為如下：

則其在任一點x，每個時刻t的展開式如下：

利用delta函數(shù)的積分性質(zhì)得到：

則（12）式可表示如下：

（11）式可表示為

直接計算核函數(shù)Ls需要將每一個位置x′作為源都進行一次正演計算.但在本文的全波反演過程中，并沒有直接使用敏感算子Ls，而是將它應(yīng)用在殘差的處理上，在下文中將會介紹，在殘差計算時僅需解決一次正演問題.由于誤差方程（7）式并沒有考慮全部空間的波場（只考慮了接收點位置的正演數(shù)據(jù)與觀測數(shù)據(jù)之差），在下文中，將使用算子Fs，該算子線性化描述所有接收點和觀察時間的電場：

其中，

對于算子Fs，并不能像Ls那樣簡單定義.在下文中繼續(xù)使用Ls代替Fs進行梯度計算.值得注意的是Fs區(qū)間是Ls區(qū)間的一個子集，因此，計算Fs區(qū)間上等同于et al.，2010）.

誤差方程梯度的推導(dǎo)使用的是擾動目標(biāo)函數(shù)的一階近似.即

將擾動目標(biāo)函數(shù)正常展開：

通過比較（23）式與（24）式，并對梯度使用與（21）式相同的線性化手段，很容易得到目標(biāo)函數(shù)的梯度為

這里引進如下的接收點位置正演的電場與實測場的殘場標(biāo)記方法：

如2.1節(jié)所述的Fs及Ls的性質(zhì)（即在接收點xd，接收時間為τ的Ls區(qū)間即為Fs的區(qū)間），得到：

梯度找到后，將Ls的轉(zhuǎn)置應(yīng)用到殘場.更確切的說，每一個元素（即每一個發(fā)射源、每一個接收位置、每一個接收時刻的全尺寸空間的介電常數(shù)或電導(dǎo)率）的梯度通過Ls與“殘場源”［ΔEs］d，τ的內(nèi)積獲得.由于這個原因，本文之前推導(dǎo)了Ls詳細的表達式.

總結(jié)（13）式及（14）式可以寫出如下梯度式：

由于敏感因子Ls的計算需要正演，將Es項放在格林函數(shù)的左邊，對“殘場源”使用格林函數(shù)的轉(zhuǎn)置.這樣避免了直接對敏感因子Ls進行求解.上述提及的表達式很具有概括性，允許各種形式的數(shù)據(jù)反演（即允許任意的源和接收器的方向組合，也就是跨孔和井地測量（VRP））.

經(jīng)過如上述變換，由（28）式可得

殘差（“殘場源”）之和Rs為：

在（29）式中，Es表示在空間介質(zhì)中麥克斯韋方程的電場解，而G＾TRs可以解釋成在相同介質(zhì)中后向傳播的矢量場.（29）式簡單的形式說明梯度是通過將所有源產(chǎn)生的所有時刻及空間三個維度（x，y，z）的場值與“殘場源”產(chǎn)生相應(yīng)時刻與空間場值的內(nèi)積之和，該內(nèi)積由前向傳播向量場與后向傳播殘向量場相乘構(gòu)成.該內(nèi)積涉及整個時空域的Es、G＾TRs.對應(yīng)的空間梯度分量由空間delta函數(shù)δ（x-x′）表示.

積分算子及其轉(zhuǎn)置是通過其核函數(shù)的性質(zhì)定義的，而線性算子及其轉(zhuǎn)置擁有相同的核函數(shù)，唯一的區(qū)別為變量是積分還是相加（Tarantola，2005）.

在2.1節(jié)中，定義了算子G ＾對普通源Js的算法，為了給出GT清晰的說明，需要探討G ＾的細節(jié).對于源Js，接收點d，觀測時間為τ，（6）式詳細的給出了求解得到的電場的每個分量（i=x，y，z k=x，y，z）：

此時應(yīng)當(dāng)非常注意各個變量的位置以及格林函數(shù)、源和場的解對應(yīng)的角標(biāo).算子G ＾=Gik（x，t，x′，t′）是全部空間V（x′）和時間域t′相加產(chǎn)生指定時空域電場Esi（x，t）的一個函數(shù).而轉(zhuǎn)置G ＾T=Gik（x′，t′，x，t）是指定時空域產(chǎn)生全部空間和時間域的函數(shù).現(xiàn)定義一個新的矢量函數(shù)T，作為對殘差（即“殘場源”）［ΔEs］d，τ的作用：

顯示各個維度分量為

（33）式的形式與（31）不同，除了沒有時間和空間的積分（我們考慮在時間域及空間域是一個delta函數(shù)，將這樣一個積分作用到（33）式即與（31）式積分形式一致），格林函數(shù)的角標(biāo)也與源（即“殘場源”）的角標(biāo)不同，時空的自變量也發(fā)生了變換.

為了更好的解釋（33）式，這里使用互易定理（Harrington et al.，2001）：

即格林函數(shù)下標(biāo)變換的同時，函數(shù)變量表示空間位置的兩元素變換位置，此時函數(shù)值不變.因此，重寫（33）式，可得：這樣就可以解釋殘差為什么可以作為一個源的傳播問題.空 間 分 量（x′，t′）以及Gki張 量 與 向 量（Es（xd，τ）-（xd，τ））i的關(guān)系，就如（31）式中與源這一項關(guān)系一樣.另外由于該問題具有時間的不變性，這意味著格林張量滿足如下關(guān)系：

于是（35）式改寫如下：

此時引入Ts，d，τ為三分量波場的矢量形式，可以將（29）改寫為

將其展開，得到

該等式可以理解為模型內(nèi)任意點在時間域正向與反向傳播矢量場零相位互相關(guān)，利用delta函數(shù)積分性質(zhì)可得：Δ

利用分布性質(zhì)，可以將d和τ的求和移動到（40）式的積分里面，最后得到梯度的表達式：

Δ

需要注意的是，此時“殘差源”Ts，d，τ為多點源，即每一個接收點位置都有一個源，在時間上為逆時傳播.

2.3 迭代步長的求解

一旦確定了梯度方向，下一步需要確定介電常數(shù)及電導(dǎo)率的迭代步長.理論上，介電常數(shù)和電導(dǎo)率迭代步長應(yīng)遵循公式

本文中使用最優(yōu)化的方法能夠找到迭代步長ζ（Pica et al.，1990），這里涉及到尋找一個沿負梯度方向的最小目標(biāo)函數(shù)，即：

方程（43）中，當(dāng)ε、σ、Sε、Sσ確定時只有唯一的獨立變量（也就是迭代步長ζ），最小值通過目標(biāo)函數(shù)一階導(dǎo)數(shù)等于零很容易求得：

解（44）式可得（Pica et al.，1990）：

這里，Es（ε，σ）是模型參數(shù) （ε，σ）的正演波場，κ是經(jīng)驗建立的擾動量（穩(wěn)定因子），Es（ε+κΔSε，σ+κΔSσ）是計算介電常數(shù)及電導(dǎo)率在其各自的負梯度方向的正演電場（這個過程需要一個額外的時間域有限差分模擬（FDTD）），介電常數(shù)和電導(dǎo)率敏感度的巨大差異會導(dǎo)致復(fù)雜模型下同時反演的失敗.Meles等（2010）通過建立一個特殊的模型，進行三次數(shù)值模擬實驗量化這些差異，在同一模型下三種不同的擾動方式如下：

1.介電常數(shù)與電導(dǎo)率同時進行微擾Es（ε+κΔSε，σ+κΔSσ）-Es（ε，σ）；

2.純電導(dǎo)率擾動Es（ε，σ+κΔSσ）-Es（ε，σ）；

3.純介電常數(shù)擾動Es（ε+κΔSε，σ）-Es（ε，σ）.

試驗說明僅使用電導(dǎo)率微擾的振幅非常小，與使用其他兩種微擾方式存在兩位數(shù)的差異，因此式（45）的步長求取方法主要取決于介電常數(shù)的梯度，而電導(dǎo)率的梯度基本不起作用.Ernst提供了一個級聯(lián)方式反演（Ernst et al.，2007），該方法是對介電常數(shù)進行迭代反演時，保持電導(dǎo)率不變；類似的，對電導(dǎo)率進行反演時，保持介電常數(shù)不變.在介電常數(shù)反演中為了最小化電導(dǎo)率異常的影響，將模型每一道數(shù)據(jù)除以其振幅能量，使得殘場數(shù)據(jù)歸一化，而在電導(dǎo)率反演中則采用正常的殘場數(shù)據(jù).

相比之下，Tarantola（2005）建議對于不同類型的模型參數(shù)采用不同的迭代步長.對于GPR數(shù)據(jù)，迭代步長采用如下模式：

因此另外一種反演方法為一個介電常數(shù)和電導(dǎo)率在每一次迭代中的準(zhǔn)同步反演，該方法通過各自沿ΔSε、ΔSσ方向?qū)ふ覙O值點，目標(biāo)函數(shù)為

得到：

其中κε、κσ為不同的小穩(wěn)定因子.對于小穩(wěn)定因子κε、κσ要很小心的取其上下界限.這些穩(wěn)定因子在反演中需要更新，在本文數(shù)值實驗中，κε=10-9，κσ=105，盡管該方法需要一個額外的兩次正演模型計算，但是同步反演過程的性質(zhì)減少了所需FDTD的總數(shù)，因為在一次迭代中介電常數(shù)和電導(dǎo)率可以同時更新.本文全波反演過程在每個源位置每次迭代需要計算四次正問題：第一次計算正演數(shù)據(jù)，第二次計算梯度方向殘差，另外兩次計算迭代步長.與Ernst相比，其在一次正演模型中，并不能將介電常數(shù)和電導(dǎo)率同時算出，該方法共需6次正演計算.此外，采用級聯(lián)方式，需要去選擇變換介電常數(shù)和電導(dǎo)率時介電常數(shù)／電導(dǎo)率的迭代次數(shù)，而同步迭代則不需要，所以同步迭代實現(xiàn)起來更簡單一些.流程圖如圖1所示，其中res-FDTD表示殘場FDTD計算，add-FDTDε、add-FDTDσ分別表示計算介電常數(shù)和電導(dǎo)率迭代步長時額外的兩次FDTD計算（Reiter and Rodi，1996）.

2.4 對數(shù)化物性參數(shù)

推導(dǎo)的反演算法是基于自然參數(shù)ε，σ和μ0的，但是不同的參數(shù)化形式讓計算（也就是在模型空間中變換坐標(biāo)）更容易實施.通常采用一種對數(shù)的表示方法：

其中，ε0是真空中介電常數(shù)，σ0為任意的電導(dǎo)率，這里采用1S·m-1，因為模型空間為線性空間，使用對數(shù)化表示參數(shù)在反演中很重要（Tarantola，2005）.此外，該參數(shù)化能夠確保介電常數(shù)和電導(dǎo)率為正值，且可以壓縮模型的數(shù)值，使其可以適應(yīng)更大的數(shù)值跨度（Ernst et al.，2007）.在反演中這樣表達的變化僅僅影響到梯度方向的表達式.為求解新的梯度方程，首先引入一個新的目標(biāo)函數(shù)表達式：

圖1 跨孔雷達全波反演計算流程Fig.1 Flow chart of cross-hole GPR full-wave inversion

將新梯度函數(shù)S′ε＾，S′σ＾表示成Sε，Sσ的函數(shù).表達式如下：

根據(jù)（52）、（53）式的定義，從（29）式可得：

Δ

2.5 全波反演計算效率及分布式算法

MATLAB Distributed Computing Engine（MATLAB分布式計算引擎，MDCE）及Parallel Computing Toolbox可以開發(fā)運行基于機群系統(tǒng)的分布式及并行MATLAB算法，且無需改變原有的MATLAB科學(xué)計算開發(fā)環(huán)境（MathWork Product Documentation，2011）.

在2.2及2.3中提出了完整的求解介電常數(shù)及電導(dǎo)率全波反演的梯度和迭代步長算法，能夠處理3D或者2D問題.然而，目前計算機在計算多個3D的FDTD正演時計算能力尚且嚴(yán)重不足，因此本文全波反演的反演算法是在直角坐標(biāo)系的2D介質(zhì)中.使用的是線源（在y方向上），即假設(shè)在y方向上介質(zhì)參數(shù)及源沒有變化，6個分量的EM場被分成兩對獨立的方程組（Taflove et al.，2005）.在跨孔雷達測量方式中僅需計算其中一組，即麥克斯韋方程組的TE模式（該模式電場分量始終與y向橫切，與偶極天線匹配），這里存在三個分量：Ex，Ez，Hy.

為了在計算中確保穩(wěn)定性和避免數(shù)值色散，需要考慮介質(zhì)參數(shù)和源的頻率組合.通常的要求是最短波長應(yīng)大于10個網(wǎng)格長度（Taflove and Hagness，2005）.在本文的模型中，采用中心頻率為100MHz的雷克子波，因此采用空間網(wǎng)格距離為5cm，這樣每個模型有105～106個網(wǎng)格節(jié)點.根據(jù)發(fā)射天線到接收點的距離，需要幾千個時間步長（為滿足時間采樣間隔的穩(wěn)定性）去獲得較長時間窗口的波場.由于采用UPML吸收邊界，寬度上匹配層共占20層元胞.在計算迭代方向時，所有發(fā)射點產(chǎn)生的所有網(wǎng)格點的Ex，Ez場都需要保存在內(nèi)存中，這需要產(chǎn)生大約40×NsrcGB內(nèi)存，其中Nsrc為發(fā)射點的個數(shù).由于模型空間的分辨率遠小于正演模型網(wǎng)格離散度，采用抽稀的方法可有效降低內(nèi)存需求.

Meles將3×3個正演元胞用一個反演元胞代替（Meles et al.，2010），在反演過程中并沒有很明顯的減小空間分辨率.本文中，對于不同的模型采用不同的抽稀模式，即對于跨孔測量，因其橫向分辨率不如縱向分辨率，可選擇在橫向每隔2個節(jié)點取一個數(shù)，在縱向上盡量少抽稀，可以每隔一個節(jié)點保留一個電場值，在實際情況中這降低了大概一個量級的內(nèi)存需求.在計算中應(yīng)當(dāng)盡量避免內(nèi)存交換過程，提高單個機群節(jié)點的內(nèi)存使用，保證計算機的正常運算.由于每個發(fā)射點位置的計算相互獨立，計算程序可以在分布式計算機群上高效率的計算，該機群可以將每個發(fā)射點放置到一個CPU節(jié)點的核中.由于分派數(shù)據(jù)額外的傳輸時間消耗僅相當(dāng)于一個正演過程的20%.因此，如果盡量避免數(shù)據(jù)通訊，完成一個完整的同步反演的計算時間Tcomp表示如下：

其中Tforward是單次正演計算所需要的時間，Niter為迭代的次數(shù).

3 算例實現(xiàn)

3.1 小尺寸規(guī)則目標(biāo)體成像

上文中已經(jīng)系統(tǒng)的介紹了全波形反演層析成像技術(shù)的原理，下面首先以單介電常數(shù)反演為例，實現(xiàn)全波反演過程.模型1如圖2所示為圓柱目標(biāo)體模型，其在2D界面上為一個圓，圍巖相對介電常數(shù)為εr=5.5、電導(dǎo)率為σ=0.0001S·m-1.目標(biāo)體直徑為1m，垂直方向位于深度3m的位置，水平方向位于3m的位置，目標(biāo)體的介電常數(shù)為εr=7，電導(dǎo)率與圍巖相同的，該目標(biāo)體位于兩個6m深的鉆孔中間，兩鉆孔間距為6m.左邊有13個發(fā)射裝置，其間距為0.5m，用圓圈符號表示，在右邊有13個接收裝置，使用叉號表示.

從圖2（b，c）可以看到已經(jīng)對目標(biāo)體區(qū)域進行了一個整體的勾勒，圖2d為經(jīng)過100次迭代后的全波形反演的圖像（此時目標(biāo)函數(shù)已遠小于初始模型的目標(biāo)函數(shù)（如圖3a所示））.通過對比圖2b以及圖2d可以發(fā)現(xiàn)，全波反演成像分辨率得到了顯著提高.圖3b為第7個源第7接收器接收的第100次后子波與初始模型及實測模型的對比，可以看到經(jīng)過迭代后的子波波形已經(jīng)很好地與實測數(shù)據(jù)匹配.圖3c為沿深度為3m的A-A（如圖2a所示）所作切線，圖3d為沿水平位置3m的B-B線（如圖2a所示）所作的切線.圖中，黑色曲線為理論曲線，藍色曲線為經(jīng)過一次迭代后曲線，紅色曲線為經(jīng)過100次迭代全波反演的波形圖，可以看到經(jīng)過反演后的結(jié)果與觀測的子波匹配的很好.從圖3（c，d）可以看出，在水平方向上，異常高值區(qū)域在圖中跨度較大，可以看出跨孔方法在垂直方向上的分辨率要優(yōu)于水平方向的分辨率.

3.2 大模型成像

模型2為如圖4a所示，發(fā)射井與接收井深度均為20.5m，兩井之間距離10m，目標(biāo)體為矩形，其中心點位于垂直距離10m、水平距離5m的位置，該矩形區(qū)域為正方形，邊長為2m，其介電常數(shù)為εr=7，電導(dǎo)率為σ=0.0001S·m-1.圍巖介電常數(shù)為εr=5.5，電導(dǎo)率為σ=0.0001S·m-1.發(fā)射天線水平位置為0m，垂直位置從1m開始，到20.5m結(jié)束，共40個位置，間隔0.5m；接收天線水平位置位于10m，垂直位置從1m開始，到20.5m結(jié)束，同樣共40個位置，間隔0.5m.

圖4b為其全波反演結(jié)果，該圖能夠精確的反演出矩形目標(biāo)體的形狀及電性參數(shù)大小，并且非常精確的反演出矩形目標(biāo)體的上下兩條邊，左右兩條邊也能夠較準(zhǔn)確的與模型對應(yīng).模型2使用了大尺寸區(qū)域參與計算，這對計算機硬件要求則提高到另外一個量級，以本模型為例，模型離散間隔為0.05m，網(wǎng)格數(shù)為240×480，F(xiàn)DTD時間迭代次數(shù)為1732次，共有40個源位置，同樣采用單精度（4bit）保存，電場值保存需要內(nèi)存240×480×1732×40×4／1024／1024／1024=29.7318G，計算殘場同樣需要內(nèi)存29.7318G，共計59.4635G內(nèi)存.使用基于MDCE的并行算法后，每個worker需要承擔(dān)的內(nèi)存為59.4635G／40=1.4866G，本機群中的兩臺搭載6核的PC，其每個核對應(yīng)一個worker，因此這兩臺PC每臺需要承擔(dān)的內(nèi)存為1.4866G×6=8.9196G，雖然這兩臺電腦物理內(nèi)存為16G，但8.9196G對其來說也是一個較大的載荷.這里使用了抽稀的方法降低內(nèi)存占用空間，由于刻畫該矩形異常體需要的分辨率數(shù)遠低于FDTD的網(wǎng)格數(shù)，因此可以將介電常數(shù)的網(wǎng)格數(shù)抽稀.由于跨孔測量方法在縱向分辨率上表現(xiàn)更為出色，因此在縱向上每隔1個點抽取一個相對介電常數(shù)值，而在橫向上每隔兩個點抽取一個相對介電常數(shù)值，如此網(wǎng)格數(shù)量上將減少到原來的1／6，即每臺PC占用內(nèi)存1.4866G，這將大大降低內(nèi)存的占用.

圖5給出了由第20個源（Depth：10.5m）發(fā)射，介電常數(shù)反演后所有接收器接收的波形與初始模型及實測模型的對比.從圖中可以明顯的看到兩者之間的差別.圖5b中能夠看出反演后結(jié)果與實測數(shù)據(jù)重合的非常好.為了進一步觀察殘差，圖5c中可以看到相比初始模型與實測數(shù)據(jù)的殘差，反演后結(jié)果與實測數(shù)據(jù)殘差幾乎為零.

3.3 復(fù)雜及隨機模型成像

為了更好的體現(xiàn)全波形層析成像的能力，建立如圖6a所示的復(fù)雜模型3，發(fā)射井與接收井深度為20.5m，兩井之間距離10m，異常體位置及形狀，深藍色表示介電常數(shù)為εr=5，淡藍色表示介電常數(shù)為εr=5.5，綠色表示介電常數(shù)為εr=6，深紅色表示介電常數(shù)為εr=7.所有異常體電導(dǎo)率均為σ=0.0001S·m-1.發(fā)射天線水平位置為0m，垂直位置從1m開始，到20.5m結(jié)束，共40個位置，間隔0.5m；接收天線水平位置位于10m，垂直位置從1m開始，到20.5m結(jié)束，同樣共40個位置，間隔0.5m.

圖6b為其反演結(jié)果，從上到下依次可得：首先垂直深度2m處可以清晰的分辨層狀分界面，垂直深度6m附近的帶狀起伏薄層能夠準(zhǔn)確地反演出起伏的微觀形態(tài)，能夠清晰的分辨出垂直深度10m的長方形高值異常，對于14m深度附近水平傾角不是很大的傾斜異常，同樣可以清晰的反演出結(jié)果.底部18m以下的垂直分界面，在垂直深度18m附近位置可以清晰的觀察到其分界線，但越往深處，分界面逐漸被均質(zhì)介質(zhì)代替，原因為觀測系統(tǒng)對該類別界面無法識別.

圖2 （a）圓柱目標(biāo)體示意圖；（b）初至?xí)r射線追蹤法反演結(jié)果；（c）第一次迭代后介電常數(shù)；（d）100次迭代后全波反演結(jié)果Fig.2 （a）Schematic diagram of cylinder body；（b）Inversion result of first-arrival ray method；（c）The permittivity after first iteration；（d）Inversion results after 100iterations

圖3 （a）目標(biāo)函數(shù)；（b）子波匹配情況；（c）深度3m介電常數(shù)切線（沿A-A所做切線）；（d）水平3m介電常數(shù)切線（沿B-B所做切線）Fig.3 （a）Objective function；（b）The wavelet matching conditions；（c）The permittivity in the depth of 3M （along A-A line）；（d）The permittivity in the distance of 3M （along B-B line）

圖4 矩形目標(biāo)體（a）及其反演結(jié)果（b）Fig.4 The model（a）and inversion results（b）of rectangular object

隨機模型（模型4）如圖6c所示，該隨機模型的水平相關(guān)長度為5m，垂直相關(guān)長度為1m.相對介電常數(shù)最小值4.8，最大值6.2.電導(dǎo)率值均為σ=0.0001S·m-1.發(fā)射天線水平位置為0m，垂直位置從1m開始，到20.5m結(jié)束，共40個位置，間隔0.5m；接收天線水平位置位于10m，垂直位置從1m開始，到20.5m結(jié)束，同樣共40個位置，間隔0.5m.

從圖6d反演結(jié)果中可以看出對于該尺寸相關(guān)長度的隨機介質(zhì)，本文的反演方法能夠較清晰的分辨出0.5m以上的薄層.高介電常數(shù)值能夠清晰看到的部分為（從上到下）：3.5m深度反演區(qū)域左側(cè)的高值薄層、5m深度附近薄層、8m深度附近薄層、14m深度靠近發(fā)射源位置的高值區(qū)域，18m深度右側(cè)部分.而14m右側(cè)部分兩個相連的更薄層則只能分辨出相對高值的一個.對于小于0.5m的微觀結(jié)構(gòu)，暫不能很好的區(qū)分開來.該反演結(jié)果能夠很準(zhǔn)確的完成層狀結(jié)構(gòu)劃分的問題.

3.4 介電常數(shù)與電導(dǎo)率分別成像

模型5物性參數(shù)分別如圖7（a，c）所示，發(fā)射天線（用圈形符號表示）位于水平距離0m位置，深度從0m到20.5m，共40個發(fā)射天線，間距為0.5m；接收天線（用叉形符號表示）位于水平距離10m位置，深度同樣從0m到20.5m，共40個接收天線，間距為0.5m.

圖5 第20個源（Depth：10.5m）發(fā)射，介電常數(shù)反演后所有接收器接收的波形與初始模型及實測模型的對比Fig.5 For（a）and（b），the transmitter located at position 10.5mdepth in Fig.4

圖6 （a）（b）復(fù)雜模型（模型3）及其反演結(jié)果；（c）（d）水平相關(guān)長度5m，垂直相關(guān)長度1m（模型4）及其反演結(jié)果Fig.6 （a）（b）Complex model（model 3）and its inversion results；（c）（d）Model of correlation lengths is 5mand 1m （model 4）in xand z directions and its inversion results

圖7（b，d）分別為模型5的介電常數(shù)和電導(dǎo)率反演結(jié)果（均為10次反演迭代），介電常數(shù)反演結(jié)果能夠清晰的分辨地層，以及大小尺寸的管線和空洞.而對于電導(dǎo)率反演從圖中可以清晰的看出位于4m位置的水平界面以及位于16m附近的起伏地表，圓柱管線位置未能有效的區(qū)分開來，與介電常數(shù)反演相比，12.5m深度的空洞反演形狀略有畸變.且在沒有射線穿過0m以上和20.5m以下區(qū)域，也能夠反演出結(jié)果.

3.5 介電常數(shù)與電導(dǎo)率同時成像

圖7 大尺寸復(fù)雜模型（模型5）及介電常數(shù)與電導(dǎo)率分別反演結(jié)果Fig.7 Large-scale complex model（model 5）and the inversion results of permittivity and conductivity

圖8 結(jié)合VRP測量方式考慮地表的模型（模型6）Fig.8 The model combined with VRP measurement and considered the surface（model 6）

本節(jié)中，將對介電常數(shù)及電導(dǎo)率同時進行反演，并且加入了井地測量（VRP）測量系統(tǒng).為了有效解決初始模型問題，本文提出采用第一次介電常數(shù)梯度作為同時反演成像的初始模型，能夠有效提高收斂速度.圖8中該模型共進行26次迭代將目標(biāo)函數(shù)計算至初始目標(biāo)函數(shù)的1%以下水平.從圖中的介電常數(shù)和電導(dǎo)率反演結(jié)果可以看出，電導(dǎo)率的反演結(jié)果對管線的成像更清晰；加入VRP測量模式后，橫向分辨率有所提高，但也造成了管線成像形狀的畸變，同樣觀察11.6m深度的空洞成像，可看出其橫向分辨率有較大的提高，但也發(fā)生了一定程度的畸變.在接近下底面的位置，電導(dǎo)率反演出現(xiàn)了很大的區(qū)域不收斂，相應(yīng)的介電常數(shù)反演模型中同樣存在這樣的問題，在接近底部位置出現(xiàn)了較大的假高值區(qū)域.同樣因為電導(dǎo)率的問題，造成了介電常數(shù)反演結(jié)果在底部出現(xiàn)了邊緣效應(yīng).

反演結(jié)果可以看出電導(dǎo)率對管線和空洞的成像效果要好于介電常數(shù)成像.

4 結(jié)論

本文全面推導(dǎo)了全波形跨孔雷達層析成像反演方法，通過基于局域網(wǎng)的分布式并行算法，有效地解決了巨量數(shù)據(jù)的正演計算，將目前僅能在超級計算機上進行計算的全波形跨孔層析成像方法在普通計算機上全面有效的實現(xiàn).正演中采用了基于UPML吸收邊界的TE模式的二維FDTD算法.反演中通過應(yīng)用包括時間維度的全波場信息與殘場逆向傳播的全波場信息乘積計算梯度方向，且通過求取以步長為自變量的目標(biāo)函數(shù)的極值確定步長公式.文中對介電常數(shù)及電導(dǎo)率分別采用對數(shù)化處理，并推導(dǎo)出新的梯度公式.在不降低分辨率的情況下，通過對反演模型采用抽稀的辦法，將數(shù)據(jù)減小至抽稀前內(nèi)存需求的1／6，解決了大型模型對計算機內(nèi)存的巨大需求.通過構(gòu)建矩形目標(biāo)體模型可以發(fā)現(xiàn)全波反演能夠精確的反演出矩形目標(biāo)體的棱邊位置、電性參數(shù)等.通過建立水平相關(guān)長度及垂直相關(guān)長度各異的隨機介質(zhì)模型可以發(fā)現(xiàn)，在100MHz天線頻率下，介質(zhì)介電常數(shù)3～8之間時，能夠清晰的分辨出0.5m以上的水平薄層異常體.在介電常數(shù)及電導(dǎo)率同步反演中，本文提出將第一步單介電常數(shù)反演作為初始模型，有效地解決了同步反演收斂速度較慢的問題.

Belina F A，Irving J，Ernst J，et al.2012a.Evaluation of the reconstruction limits of a frequency-independent crosshole georadar waveform inversion scheme in the presence of dispersion.Journal of Applied Geophysics，78：9-19.

Belina F A，Irving J，Ernst J，et al.2012b.Analysis of an iterative deconvolution approach for estimating the source wavelet during waveform inversion of crosshole georadar data.Journal of Applied Geophysics，78：20-30.

Belina F A，Irving J，Ernst J，et al.2012c.Waveform inversion of crosshole georadar data：influence of source wavelet variability and the suitability of a single wavelet assumption.IEEE Transactions on Geoscience and Remote Sensing，50（11）：4610-4625.

Cai W Y，Qin F H，Schuster G T.1996.Electromagnetic velocity inversion using 2-D Maxwell′s equations.Geophysics，61（4）：1007-1021.

Cui T J，Qin Y，Wang G L，et al.2004.Low-frequency detection of two-dimensional buried objects using high-order extended born approximations.Inverse Problem，20（6）：S41-S62.

Ernst J R，Maurer H，Green A G，et al.2007.Full-waveform inversion of crosshole radar data based on 2-D finite-difference time-domain solutions of maxwell′s equations.IEEE Transactions on Geoscience and Remote Sensing，45（9）：2807-2828.

Greenhalgh S A，Bing Z，Green A.2006.Solutions algorithms and inter-relations for local minimization search geophysical inversion.Journal of Geophysics and Engineering，3（2）：101-113.

Harrington R F.2001.Time-Harmonic Electromagnetic Fields.New York：Wiley，IEEE Press.

Johnson T C，Routh P S，Knoll M D.2005.Fresnel volume georadar attenuation-difference tomography.Geophys.J.Int.，162（1）：9-24.

Kuroda S，Takeuchi M，Kim H J.2005.Full waveform inversion algorithm for interpreting cross-borehole GPR data.∥Proc.SEG Int.Expo.And 75th Annu.Meeting.Houston，TX，1176-1179.

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.

MathWorks-Product Documentation.2011.Accelerating the pace of engineering and science.http：∥ www.mathworks.cn／help／techdoc／ref／f16-6011.html.

Meles G A，Greenhalgh S A，Green A G，et al.2012.GPR Full-Waveform sensitivity and resolution analysis using an FDTD adjoint method.IEEE Transactions on Geoscience and Remote Sensing，50（5）：1881-1896.

Meles G A，Greenhalgh S A，Kruk J V D，et al.2011.Taming the non-linearity problem in GPR full-waveform inversion for high contrast media.Journal of Applied Geophysics，73（2）：174-186.

Meles G A，Kruk J V D，Stewart 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.IEEETransactionson GeoscienceandRemoteSensing，48（9）：3391-3407.

Pica A，Diet J P，Tarantola A.1990.Nonlinear inversion of seismic reflection data in a laterally invariant medium.Geophysics，55（3）：284-292.

Reiter D T，Rodi W.1996.Nonlinear waveform tomography applied to Crosshole seismic data.Geophysics，61（3）：902-913.

Tarantola A.2005.Inverse Problem Theory and Methods for Model Parameter Estimation.Philadelphia，PA：Soc.Ind.and Appl.Math.

Taflove A，Hagness S C.2005.Computational Electrodynamics：The Finite-Difference Time-Domain Method.3rd ed.Norwood，MA：Artech House.

Wu J J.2012.Research on cross-hole radar tomography using fullwaveform inversion method［Ph.D.thesis］.Changchun：Jilin University.

Zhou C G，Liu L B.2000.Radar-diffraction tomography using the modified quasi-linear approximation.IEEE Transactions on Geoscience and Remote Sensing，38（1）：404-415.

Zhou H，Sato M，Takenaka T，et al.2007.Reconstruction from antenna transformed radar data using a time-domain reconstruction method.IEEE Transactions on Geoscience and Remote Sensing，45（3）：689-696.

附中文參考文獻

吳俊軍.2012.跨孔雷達全波形層析成像反演方法的研究［博士論文］.長春：吉林大學(xué).