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

劉廣東

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

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

劉廣東

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

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

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

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

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

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

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

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

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

1 問題描述

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

1.1 正問題

正問題：已知激勵源和媒質的電磁參數，求時域電磁場。

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

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

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

1.2 逆問題

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

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

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

其次，定義復值的目標函數為：

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

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

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

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

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

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

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

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

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

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

2 逆問題解法

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

為了抑制逆問題的病態特性，引入一階吉洪諾夫（Tikhonov）正則化技術［14］。基于最小二乘準則，定義最小化問題的泛函為：

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

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

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

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

3 應用算例

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

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

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

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

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

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

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

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

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

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

4 結束語

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

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

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

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

［1］鄭圣談,曾昭發,劉四新,等.寬帶高頻電磁場數據反演方法研究［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］劉廣東,張業榮.一種處理分層有耗色散介質的時域逆散射方法［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，Fuyang Teachers College，Fuyang 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