二維時(shí)域變形玻恩迭代法及其在地下勘探中的應(yīng)用

劉廣東

（阜陽師范學(xué)院物理與電子科學(xué)學(xué)院，安徽阜陽236037）

二維時(shí)域變形玻恩迭代法及其在地下勘探中的應(yīng)用

劉廣東

（阜陽師范學(xué)院物理與電子科學(xué)學(xué)院，安徽阜陽236037）

為了方便直接利用時(shí)域測量數(shù)據(jù)，重構(gòu)地下目標(biāo)的電參數(shù)，將二維頻域變形玻恩迭代法（BDIM）發(fā)展到時(shí)域形式。迭代反演過程中，正問題用時(shí)域有限差分（FDTD）法求解；而逆問題的病態(tài)特性由吉洪諾夫（Tikhonov）正則化技術(shù)抑制。在地下勘探應(yīng)用算例中，噪聲污染被加入，數(shù)值結(jié)果初步證實(shí)了改進(jìn)方法的可行性。

變形玻恩迭代法；探地雷達(dá)；時(shí)域有限差分；吉洪諾夫正則化

地下勘探是指利用探地雷達(dá)（ground penetrating radar，GPR）發(fā)射高頻寬帶脈沖電磁波，對(duì)地下目標(biāo)進(jìn)行探測的地球物理方法，工程上可廣泛用于地下水調(diào)查、地下固體廢物與污染物探測、考古、土壤評(píng)價(jià)及城市管線探測等許多方面［1］。其物理基礎(chǔ)是目標(biāo)（object of interest，OI）和背景媒質(zhì)的電參數(shù)（介電常數(shù)和電導(dǎo)率）差異。因此，地下勘探屬于一類電磁逆散射問題［2］。

電磁逆散射（也常稱為反演、重建、重構(gòu)等）是一類富有挑戰(zhàn)的不適定問題，非線性和病態(tài)性是其中的兩個(gè)主要困難［2-4］。發(fā)展至今，尚未形成一套通用的解決方案。常用的做法是采樣線性化近似結(jié)合迭代法處理非線性問題［5］，正則化技術(shù)抑制病態(tài)特性［6-7］。

從方法論角度，處理電磁逆散射問題，可以從頻域進(jìn)行，相應(yīng)方法主要有高斯－牛頓法（Gauss－Newton，GN）［7］、玻恩迭代法（Born iterative method，BIM）［8］和變形玻恩迭代法（distorted Born iterative method，DBIM）［9］；也可以從時(shí)域進(jìn)行，相應(yīng)方法主要有正反時(shí)間步進(jìn)法（forward－backward time－stepping，F(xiàn)BTS）［10］、時(shí)域玻恩迭代法（BIM）［11］和時(shí)域高斯－牛頓法（GN）［12］。研究表明：頻域方法對(duì)窄帶信號(hào)而言可行、經(jīng)濟(jì)；而對(duì)高頻寬帶脈沖而言，時(shí)域方法則顯示了信息量大、重建分辨率高等優(yōu)勢［10］。

為了直接利用時(shí)域測量數(shù)據(jù)，文獻(xiàn)［11］將玻恩迭代法（BIM）從頻域發(fā)展到時(shí)域，重建無耗媒質(zhì)的介電常數(shù)分布圖像具有超分辨率。

然而，文獻(xiàn)［9］指出：變形玻恩迭代法（DBIM）的收斂速度快于玻恩迭代法（BIM），更適合雙電參數(shù)（介電常數(shù)和電導(dǎo)率）同時(shí)重建。

利用上述結(jié)論，并借鑒文獻(xiàn)［11］的思想，本文將二維（two－dimensional，2－D）變形玻恩迭代法（DBIM）從頻域改進(jìn)到時(shí)域，并應(yīng)用到地下勘探的數(shù)值算例中。迭代過程中，每步產(chǎn)生的正問題由時(shí)域有限差分法（finite－difference time－domain，F(xiàn)DTD）求解［13］。

1 問題描述

首先，假設(shè)所有媒質(zhì)是各向同性、線性和無磁的（即相對(duì)磁導(dǎo)率μr＝1）。其次，直角坐標(biāo)系下，設(shè)在平行于z軸的地下淺層包含若干個(gè)柱狀異常體，被置于地表（或者近地表）且用平行于z軸的線電流源激勵(lì)的發(fā)射天線T產(chǎn)生的橫磁波（transverse magnetic，TM）照射，產(chǎn)生的總電場由接收天線R接收，x－y橫截面如圖1所示，其中，S表示未知電參數(shù)的不均勻目標(biāo)（OI）區(qū)域，Ω表示已知電參數(shù)的均勻背景媒質(zhì)區(qū)域。記媒質(zhì)的相對(duì)介電常數(shù)和電導(dǎo)率分別為εr和σ。M根發(fā)射天線和N根接收天線分別位于和1， 2，．．．，N）位置。

1.1 正問題

正問題：已知激勵(lì)源和媒質(zhì)的電磁參數(shù)，求時(shí)域電磁場。

二維（2－D）橫磁波（TM）照射下，時(shí)域電磁場滿足麥克斯韋（Maxwell）方程組［13］

其中，位置矢量r＝（x，y），μ0和ε0分別表示真空磁導(dǎo)率和介電常數(shù)，Jz表示激勵(lì)電流源。

本文采用時(shí)域有限差分法（FDTD）離散媒質(zhì)空間，周圍采用卷積完全匹配層（Convolution perfectly matched layer，CPML）吸收邊界條件［13］。相關(guān)實(shí)現(xiàn)細(xì)節(jié)見文獻(xiàn)［13］。

1.2 逆問題

逆（或反）問題：已知測量（或仿真）的時(shí)域電磁場和均勻背景媒質(zhì)的電參數(shù)，求不均勻目標(biāo)的電參數(shù)。

首先，利用快速傅里葉變換（fast Fourier transform， FFT），獲得頻域復(fù)值的散射電場為：

其中，j表示虛數(shù)單位，角頻率ω＝2πf，f表示頻率。

其次，定義復(fù)值的目標(biāo)函數(shù)為：

其中，εr（r）和σ（r）分別表示目標(biāo)的實(shí)值相對(duì)介電常數(shù)和電導(dǎo)率，εbr和σb分別表示背景媒質(zhì)的實(shí)值相對(duì)介電常數(shù)和電導(dǎo)率。

類似于文獻(xiàn)［11］，由于不均勻散射體（目標(biāo)）的存在，可以導(dǎo)出產(chǎn)生的頻域復(fù)值散射場滿足：

其中，k0表示真空實(shí)值波數(shù)，利用對(duì)易原理，不均勻目標(biāo)區(qū)的格林（Green）函數(shù)滿足［9］：

其中，I（ω）表示激勵(lì)時(shí)域脈沖對(duì)應(yīng)的頻域形式。

離散計(jì)算空間為L＝LxLy個(gè)網(wǎng)格，其中Lx、Ly分別表示x和y方向剖分的網(wǎng)格數(shù)。采用脈沖基函數(shù)，展開目標(biāo)函數(shù)為：

其中，˙el和˙al為復(fù)值展開系數(shù)，脈沖基函數(shù)定義為［11］：

其中，sl表示第l個(gè)離散網(wǎng)格。這樣，積分方程（8）可轉(zhuǎn)化為復(fù)值的矩陣形式

其中，J×1復(fù)值列向量˙e定義為：

這里J＝MNW，W表示采樣頻點(diǎn)數(shù)，上角標(biāo)t表示轉(zhuǎn)置運(yùn)算（下文同）；J×L復(fù)值矩陣˙K的元素K·j，l為：

L×1復(fù)值列向量˙a定義為：

2 逆問題解法

為了處理逆問題的非線性，本文將矩陣方程（13）轉(zhuǎn)化為最小化問題迭代求解，迭代初始，采用玻恩近似（Born approximation，BA）線性化逆問題［11］。

為了抑制逆問題的病態(tài)特性，引入一階吉洪諾夫（Tikhonov）正則化技術(shù)［14］?；谧钚《藴?zhǔn)則，定義最小化問題的泛函為：

其中，上、下標(biāo)2分別表示平方、L2范數(shù)運(yùn)算，R表示L× L正則化矩陣，γ表示正的正則化參數(shù)，其選取辦法可以參考文獻(xiàn)［11］。最小化問題等價(jià)于法方程：

其中，上角標(biāo)H表示復(fù)共軛轉(zhuǎn)置運(yùn)算。該方程可用共軛梯度法（conjugate gradient，CG）迭代求解［15］。

定義第i步迭代時(shí)的相對(duì)剩余誤差為：

時(shí)域變形玻恩迭代法（DBIM）的計(jì)算流程如圖2所示。

3 應(yīng)用算例

為了檢驗(yàn)本文方法的性能，應(yīng)用它進(jìn)行地下勘探，即利用GPR測量的時(shí)域電場重建地下二維（x－y截面）目標(biāo)的相對(duì)介電常數(shù)和電導(dǎo)率分布。地質(zhì)雷達(dá)測量的幾何結(jié)構(gòu)如圖3所示。其中，粗實(shí)線正方形區(qū)域?yàn)殡妳?shù)未知的目標(biāo)區(qū)域（目標(biāo)區(qū)域中含有四個(gè)等大的正方形異常體S1、S2、S3和S4，邊長等于0.45 m），正方形邊長等于4.50 m，周圍為已知電參數(shù)的背景媒質(zhì)（空氣的相對(duì)介電常數(shù)和電導(dǎo)率分別為εr＝1，σ＝0；大地的相對(duì)介電常數(shù)和電導(dǎo)率分別為εr＝4，σ＝1.0×10－3S／m），邊界間距1.50 m。

圖3中，黑色圓點(diǎn)表示GPR的8個(gè)均勻分布的天線陣元（點(diǎn)模型），相鄰間距為0.45 m，置于水平地面上方0.75 m，各陣元被依次激活，當(dāng)其中某個(gè)陣元用作發(fā)射天線時(shí)，則其余7個(gè)用作接收天線，即M＝8，N＝7。激勵(lì)的時(shí)域信號(hào)選為

其中，fc1＝20 MHz，fc2＝60 MHz。即輻照時(shí)間T＝40 ns，頻點(diǎn)數(shù)W＝2。

迭代過程中，每步產(chǎn)生的正問題由時(shí)域有限差分法（FDTD）求解［13］，剖分網(wǎng)格尺寸為Δx＝Δy?Δ＝ 0.15 m，時(shí)間步長為，其中c0為真空光速。虛線區(qū)域周圍由10層卷積完全匹配層（CPML）吸收邊界截?cái)啵?3］。

本文利用FDTD法仿真時(shí)域電場代替測量數(shù)據(jù)。此時(shí)，目標(biāo)區(qū)域中各異常體的電參數(shù)見表1。為了讓仿真更接近真實(shí)測量，避免過于樂觀的重建，特作幾點(diǎn)嘗試：（1）剖分采用兩倍精細(xì)網(wǎng)格；（2）在仿真數(shù)據(jù)中加入均勻分布的隨機(jī)噪聲［3］，信噪比取20 dB；（3）大地中除了明顯的異常體以外，可能存在一些電參數(shù)差異較小的雜質(zhì)，為了考察它們的影響，取目標(biāo)區(qū)域的電特性參數(shù)在1±5%的范圍內(nèi)隨機(jī)變化。

逆問題求解過程中，選用Polak－Ribière－Polyak（PRP）共軛梯度法（CG）［15］求解方程（18），正則化矩陣R取為單位矩陣I，正則化參數(shù)取為γ＝0.001。

按照圖2所示的計(jì)算流程，迭代求解目標(biāo)函數(shù)，處理的部分結(jié)果分別如圖4～圖6所示。其中，圖4和圖5分別給出了相對(duì)介電常數(shù)εr和電導(dǎo)率σ的真實(shí)分布以及在迭代次數(shù)分別為i＝0、30和60時(shí)的重建結(jié)果，圖6是相對(duì)誤差隨迭代次數(shù)的變化關(guān)系。其中，第30和60次迭代的相對(duì)誤差分別為0.18%和0.03%。

因此，可以得出：（1）總體上，2D時(shí)域DBIM應(yīng)用到該具體算例是收斂的，收斂速度先快后慢；且隨著迭代次數(shù)的增加，相對(duì)誤差逐漸減小，重建精度逐漸提高（當(dāng)然，迭代次數(shù)越多，計(jì)算耗時(shí)也越長，實(shí)際的工程應(yīng)用中，可以選擇一個(gè)重建精度與實(shí)時(shí)重建的合理折中）。

（2）當(dāng)?shù)螖?shù)i＝60時(shí)，盡管四個(gè)主要散射體S1、S2、S3和S4的電特性參數(shù)高于或低于背景媒質(zhì)（大地），但均收到了較好的重建效果，較好地重現(xiàn)了它們的形狀、尺寸、位置信息（各參數(shù)和相應(yīng)的真實(shí)值尚有一定的差異，適當(dāng)增加仿真時(shí)間步數(shù)或者適當(dāng)加大迭代次數(shù)是可能的解決方法）。（3）就4個(gè)主要散射體的位置而言，S1和S2對(duì)比S3和S4，前者的重建效果略差于后者，可能是由于S1和S2位于地下較深位置，距離GPR較遠(yuǎn)，信號(hào)在大地中經(jīng)受衰減，來自S1和S2的散射信號(hào)較弱所致（解決方法是適當(dāng)降低激勵(lì)源的頻率。當(dāng)然，這也同時(shí)降低了重建的分辨率，在實(shí)際的工程應(yīng)用中，應(yīng)該選擇一個(gè)檢測深度和重建分辨率的合理折中）。

綜合起來，由于測量視角的受限和測量數(shù)據(jù)的有限，增加了逆問題的病態(tài)性質(zhì)，盡管如此，在噪聲環(huán)境下，在電參數(shù)假設(shè)為1±5%隨機(jī)變化的大地中，地下（異常體）目標(biāo)仍被成功地探測到。這可能得益于本文時(shí)域DBIM的魯棒性和正則化技術(shù)等因素（當(dāng)然，正則化方案已有多種［16］，正則化參數(shù)的選擇也是個(gè)難題，基本思路是：根據(jù)具體問題，通過數(shù)值測試確定最優(yōu)參數(shù)）。

4 結(jié)束語

利用GPR發(fā)射高頻寬帶脈沖進(jìn)行地下目標(biāo)勘探，并進(jìn)行時(shí)域反演是一種重要的地球物理方法，在水文地質(zhì)、工程地質(zhì)、地質(zhì)調(diào)查和軍事探測等領(lǐng)域具有良好的應(yīng)用前景。

本文將二維（2－D）變形玻恩迭代法（DBIM）從頻域改進(jìn)到時(shí)域。迭代過程中，正問題由時(shí)域有限差分法（FDTD）求解，經(jīng)正則化的逆問題由快速傅里葉變換（FFT）結(jié)合共軛梯度法（CG）求解。

地下勘探的數(shù)值結(jié)果顯示：目標(biāo)的介電常數(shù)和電導(dǎo)率均被成功地重建，再現(xiàn)了散射體的位置、尺寸和形狀等目標(biāo)信息。本文的雙參數(shù)反演方法豐富了寬帶高頻電磁脈沖在有耗媒質(zhì)中的傳播理論，拓展了電磁勘探方法，體現(xiàn)了一定的工程應(yīng)用價(jià)值。

當(dāng)然，諸多問題尚待進(jìn)一步研究：（1）本文方法能否應(yīng)用于其它領(lǐng)域（如微波醫(yī)學(xué)成像）；（2）其它正則化方案（如自適應(yīng)正則化）的性能；（3）推廣到三維情形。

［1］鄭圣談,曾昭發(fā),劉四新,等.寬帶高頻電磁場數(shù)據(jù)反演方法研究［J］.地球物理學(xué)報(bào),2008,51（1）：266-272.

［2］Zhang Peng,Zhang Xiaojuan,Fang Guangyou.Comparison of the imaging resolutions of time reversaland backprojection algorithms in EM inverse scattering［J］.IEEE Geosci.Remote Sens.Lett.,2013,10（2）：357-361.

［3］劉廣東,張業(yè)榮.一種處理分層有耗色散介質(zhì)的時(shí)域逆散射方法［J］.電子學(xué)報(bào),2011,39（12）：2856-2862.

［4］劉廣東.乳腺癌檢測的時(shí)域微波成像算法研究［D］.南京：南京郵電大學(xué),2011.

［5］Gennarelli G,Soldovieri F.A linear inverse scattering algorithm for radar imaging in multipath environments［J］.IEEE Geosci.Remote Sens.Lett.,2013,10（5）：1085-1089.

［6］張建軍,王芹.二值圖像恢復(fù)的一個(gè)非線性正則化方法［J］.應(yīng)用科學(xué)學(xué)報(bào),2009,27（1）：74-78.

［7］Mojabi P,LoVetri J.A prescaled multiplicative regularized Gauss-New ton inversion［J］.IEEE Trans.Antennas Propag.,2011,59（8）：2954-2963.

［8］Wang Y M,Chew W C.An iterative solution of twodimensional electromagnetic inverse scattering problem［J］.Int.J.Imaging Syst.Technol.,1989,1（1）：100-108.

［9］Chew W C,Wang Y M.Reconstruction of two-dimensional permittivity distribution using the distorted Born iterativemethod［J］.IEEE Trans.Med.Imaging,1990,9（2）：218-225.

［10］Takenaka T,Jia H,Tanaka T.M icrowave imaging of electrical property distributions by a forward-backward time-stepping method［J］.J.Electromagn.Waves Appl., 2000,14（12）：1609-1626.

［11］A liM A,Moghaddam M.3D nonlinear super-resolution microwave inversion technique using time-domain data［J］.IEEE Trans.Antennas Propag.,2010,58（7）：2327-2336.

［12］劉廣東,張開銀.二維電磁逆散射問題的時(shí)域高斯-牛頓反演算法［J］.物理學(xué)報(bào),2014,63（3）：116-130.

［13］Taflove A,Hagness S C.Computational electrodynamics：The finite-difference time-domain method（3rd ed.）［M］.A rtech：House press,2005.

［14］Tikhonov A N,A rsenin V Y.Solutions of ill-posed problems［J］.SIAM Rev,1979,21（2）：266-267.

［15］Dai Y H.A family of hybrid conjugate gradientmethods for unconstrained optimization［J］.Math.Comput., 2003,17（2）：1317-1328.

［16］Mojabi P,LoVetri J.Overview and classification of some regularization techniques for the Gauss-New ton inversion method applied to inverse scattering problems［J］.IEEE Trans.Antennas Propag.,2009,57（9）：22658-22665.

Two-Dimensional Time-Domain Distorted Born Iterative Method and Its App lication in Underground Exp loration

LIU Guangdong
（School of Physical and Electronics，F(xiàn)uyang Teachers College，F(xiàn)uyang 236037，China）

The two-dimensional（2-D）distorted Born iterativemethod（DBIM）is developed from frequency domain to time domain in order to reconstruct the dielectric parameters of the underground objectand directly use the time-domainmeasured data.In the iterative inversion process，the forward problem is solved by the finite-difference time-domain（FDTD）method，while the ill-peculiarity of the corresponding inverse problem is restrained by Tikhonov's regularization technique. Furthermore，the noise pollution is considered in a numerical example of underground exploration application，and the numerical results preliminarily demonstrate the feasibility ofmodified algorithm.

distorted Born iterativemethod（DBIM）；ground penetrating radar（GPR）；finite-difference time-domain（FDTD）；Tikhonov's regularization

TN95

A

1673-1549（2014）04-0056-05

10.11863／j.suse.2014.04.14

2014-03-17

國家自然科學(xué)基金項(xiàng)目（61071022）；阜陽師范學(xué)院重點(diǎn)基金項(xiàng)目（2011FSKJ01ZD）

劉廣東（1972-），男，江蘇灌云人，副教授，博士，主要從事電磁散射和逆散射方面的研究，（E-mail）liu＿guang＿dong＠126.com