劉廣東

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

二維時域變形玻恩迭代法及其在地下勘探中的應用

劉廣東

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

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

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

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

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

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

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

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

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

1 問題描述

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

1.1 正問題

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

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

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

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

1.2 逆問題

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

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

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

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

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

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

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

其中，I（ω）表示激勵時域脈沖對應的頻域形式。

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

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

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

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

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

L×1復值列向量˙a定義為：

2 逆問題解法

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

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

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

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

定義第i步迭代時的相對剩余誤差為：

時域變形玻恩迭代法（DBIM）的計算流程如圖2所示。

3 應用算例

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

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

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

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

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

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

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

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

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

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

4 結(jié)束語

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

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

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

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

［1］鄭圣談,曾昭發(fā),劉四新,等.寬帶高頻電磁場數(shù)據(jù)反演方法研究［J］.地球物理學報,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ì)的時域逆散射方法［J］.電子學報,2011,39（12）：2856-2862.

［4］劉廣東.乳腺癌檢測的時域微波成像算法研究［D］.南京：南京郵電大學,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］張建軍,王芹.二值圖像恢復的一個非線性正則化方法［J］.應用科學學報,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］劉廣東,張開銀.二維電磁逆散射問題的時域高斯-牛頓反演算法［J］.物理學報,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

國家自然科學基金項目（61071022）；阜陽師范學院重點基金項目（2011FSKJ01ZD）

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