基于頻域有限差分的三維大地電磁的低感應數(shù)預條件求解方法研究

蔣濤 （油氣資源與勘探技術(shù)教育部重點實驗室 （長江大學），湖北 武漢430100 長江大學工程技術(shù)學院，湖北 荊州434020）

胡文寶 （油氣資源與勘探技術(shù)教育部重點實驗室 （長江大學），湖北 武漢430100）

正演數(shù)值計算的算法研究是地球物理勘探的基礎之一，只有充分了解和認識了已知場源所產(chǎn)生的地球物理異常，才可能從獲得的數(shù)據(jù)中反演引起異常的異常體電磁學性質(zhì)。常用的數(shù)值計算方法主要有基于積分方程的積分方程法、邊界單元法和基于微分方程的有限差分法、有限元方法。隨著計算機技術(shù)的發(fā)展和新算法技術(shù)的引進，近十年來大地電磁三維正反演取得了豐碩的成果［1，2］。有限差分法因其物理意義明確、簡單直觀、計算量小、計算速度快被很多學者選用［3～5］①Newman G A，Alumbaugh D L．Three-dimensional electromagnetic modeling on massively parallel computers．Sandia National Laboratories Report，SAND96-0582，1996．。譚捍東等［6］采用雙共軛梯度穩(wěn)定解法求解三維交錯采樣有限差分離散后得到的大型系數(shù)矩陣方程組，從迭代收斂穩(wěn)定、計算精度高、速度快等方面，通過2個理論模型檢驗了該算法的正確性和計算精度。肖騎彬等［7］從網(wǎng)格精度、形成系統(tǒng)方程的方式、邊界條件以及預條件線性算子等，對大地電磁有限差分數(shù)值解作了對比，結(jié)果表明，合適的預條件再配合好的線性算子，不僅可以加速收斂，而且可以降低迭代次數(shù)。李焱等［8］為了更有效地提高大地電磁三維正演的計算速度，根據(jù)各頻率之間求取電磁場值的過程是相互獨立的特點，引入并行處理技術(shù)，對2個理論模型的大地電磁響應采用三維交錯網(wǎng)格有限差分數(shù)值模擬，高效地實現(xiàn)了正演計算。在上述研究中，如何快速、準確地求解經(jīng)過有限網(wǎng)格離散后的大型稀疏線性方程組成為關(guān)鍵。

該次研究采用Yee［9］提出的交錯網(wǎng)格有限差分法導出了頻域Maxwell方程的離散化關(guān)系式，將三維電磁響應的模擬問題轉(zhuǎn)化成求解大型線性代數(shù)方程組的問題。筆者將總場分離成一次場 （背景場）和二次場 （散射場），背景場用解析方法求取，二次場用有限差分實現(xiàn)。導出的線性代數(shù)方程組采用預條件方法的Krylov子空間類迭代方法求解。

為了解決頻域低于1MHz時采用頻域有限差分方法得到的線性方程組在迭代求解過程中收斂緩慢的問題，Newman等［3，4］提出了一種能顯著減少三維頻域有限差分數(shù)值計算時間的方法，即基于電場的Helmholtz分解定理的低感應數(shù)預條件方法。Maxwell微分方程經(jīng)過頻域有限差分網(wǎng)格離散后得到的線性方程組的系數(shù)矩陣具有大型復數(shù)稀疏對稱的性質(zhì)，廣泛采用QMR （quasi-minimal residual，準最小殘差法）進行迭代求解。而低感應數(shù)預條件能顯著加快線性方程組迭代求解過程中的收斂速度，比簡單的Jacobi預條件加速時間快1／10，可以將迭代次數(shù)從未做預條件處理的上千次降低到幾十次。筆者在此基礎上，研究了在利用QMR迭代求解線性方程組的過程中，為構(gòu)建低感應數(shù)預條件矩陣，2種不同算法 （共軛梯度法和最小殘差法）對正演計算迭代次數(shù)和時間的差別，分析了原因并給出了結(jié)論。

1 頻域有限差分法

對于諧變電磁場，采用時諧因子eiωt（其中，i為虛數(shù)；ω為角頻域，rad／s；t為時間，s）在各向同性介質(zhì)中散射電場Es分布的控制方程可以從Maxwell方程組獲得：

式中：△為哈密頓算子，后面接 “×”表示旋度；ES為散射電場，V／m；μ為介質(zhì)的磁導率，H／m；σ為介質(zhì)的電導率，S／m；σ0為背景介質(zhì)的電導率，S／m；EP為背景介質(zhì)中的一次場，V／m。按照Newman等［3］的方法，式 （1）在交錯笛卡爾網(wǎng)格中離散，ES位于元胞邊的中心，利用中心差分代替微分，可得到大型病態(tài)線性方程組：

式中：A為經(jīng)過頻域有限差分離散后得到的大型復數(shù)對稱稀疏矩陣，每一行最多含有13個非零元素；E為散射電場矢量；b為與源或者邊界條件相關(guān)的矢量。線性方程組的左邊第m行可以寫成二矢量之積αTmum（其中，αm為系數(shù)矩陣A的第m行13個非零元素組成的向量；u為待求電場，V／m）。在磁導率變化不大的條件下 （μ＝μb）（其中，μb為背景介質(zhì)的磁導率，H／m），X方向的散射電場的線性方程可寫為αTmum：

式中：i、j、k分別為直角坐標X、Y、Z方向上的網(wǎng)格節(jié)點序號；Δ為系數(shù)矩陣A的對角元素；δ為相鄰網(wǎng)格節(jié)點之間的距離，m；δ＊為相鄰網(wǎng)格節(jié)點中心之間的距離，m；Ex、Ey、Ez分別為待求的散射電場三分量，V／m；x （i）、y （j）、z （k）分別為X、Y、Z方向上網(wǎng)格節(jié)點坐標，m。

系數(shù)矩陣A的對角元素為：

式中：σ（i＋1／2，j，k）為沿著元胞棱邊的4個相鄰單元按體積權(quán)重所取的電導率，S／m。

式 （4）～ （6）部分地給出了式 （2）的復數(shù)對稱系數(shù)矩陣，容易通過無矩陣內(nèi)積方法實現(xiàn)。在利用QMR的每一步迭代過程中，需要計算一次矩陣－矢量積。近似解xi收斂到實際解的速度依賴于條件數(shù)（其中，下標P為LP范數(shù)），這里系數(shù)矩陣條件數(shù)很大，所以收斂過程緩慢。

2 預條件技術(shù)

對中高頻率的大地電磁波，采用簡單的Jacobi預處理的QMR能夠有效地求解式 （2）。但是，當場源頻率低到接近似穩(wěn)場時，用頻域有限差分離散得到的方程組條件數(shù)很大，若用Krylov子空間類迭代方法求解式 （2），則因為方程的病態(tài)程度嚴重，使得迭代過程的收斂速度很緩慢。Druskin等［10］利用與頻域有限差分近似的SLDM （spectral Lanczos decomposition method）報道了求解過程中的問題。為了加快迭代收斂的速度，必須要使用預條件技術(shù) （又稱預優(yōu)技術(shù)）。

Newman等［4］和Weiss等［12］根據(jù)Helmholtz原理把電場 （E）分解成無散場 （F）與無旋場 （△f）之和，即：

將式 （7）代入式 （1），“△×△×”算符可化為矢性拉普拉斯算符：

式中：▽后不接符號表示梯度。

而式 （9）的解并不滿足電流的連續(xù)性方程，為了滿足該條件，對式 （1）取散度然后代入式 （7）可得：

式中：▽后接 “·”表示散度。

這樣，當頻率很低即低感應數(shù)條件滿足時，經(jīng)過交錯網(wǎng)格頻域有限差分離散，可采用共軛梯度類算法先后求解式 （9）、（10），可以使式 （8）的近似解達到限制的誤差范圍內(nèi)。

改善收斂緩慢的方法是，定義一預條件矩陣M，把式 （2）的求解可以有效轉(zhuǎn)換為求解（AM－1）x＝b（其中，x＝ME為解矢量）。若M為A的主對角線，則M稱為Jacobi預條件矩陣。如果A＝A1＋A2，當A1比A2占優(yōu)勢時，可取M ＝A1為預條件矩陣。在這里經(jīng)過預條件后的線性方程組變?yōu)椋?/p>

其中：A1E＝x。

這種預條件可以被應用到頻率較低的準靜態(tài)場的計算中，因為此時電場的旋度部分比起散度部分占優(yōu)勢。Newman等［4］提出的低感應數(shù)預條件方法正是基于此，將無散場部分F的擴散方程的解與無旋場▽f的泊松方程的解之和作為預條件矩陣。低感應數(shù)法對QMR迭代收斂速度的加速效果可以使迭代次數(shù)從未做預條件處理的上千次降低到幾十次，優(yōu)于不完全LU分解和簡單Jacobi預條件方法。

3 低感應數(shù)預條件矩陣的求取方法

由于線性方程組的系數(shù)矩陣A是一個大型稀疏復對稱的非埃爾米特矩陣，Krylov子空間類迭代方法是求解該類線性方程最為有效的方法。常用的求解稀疏復線性方程組的Krylov子空間類迭代方法主要是QMR，因為其穩(wěn)定性好且收斂速度也較快。在利用QMR求解線性方程組時所需的低感應數(shù)預條件矩陣的求取方法主要有2種：共軛梯度法 （conjugate gradient）和最小殘差法 （minimum residual）。該次研究對采用2種方法得到的低感應數(shù)預條件矩陣對大地電磁三維正演結(jié)果的影響進行了對比研究（2種方法的詳細內(nèi)容見文獻 ［5］）。

共軛梯度法和最小殘差法均屬于共軛方向算法的一部分，等同于使誤差函數(shù)在每次迭代過程的子空間中最小。考慮二次函數(shù)f（x）＝xTHx－h(huán)Tx（其中，H為對稱正定矩陣；h是線性方程組Hx＝h的右端向量）的最小值出現(xiàn)在＝H－1h （其中，是函數(shù)f（x）取最小值時的向量）并且值為f（）＝－h(huán)TH－1h （其中，f（＾x）是函數(shù)f（x）的最小值），因此找到二次函數(shù)的最小值相當于找到方程Hx＝h的解。

3.1 共軛梯度法

若A為對稱正定，對于共軛梯度法有H＝A和h＝b。因此在最小化f （x）的每一步驟中，解矢量x和殘差矢量r＝b－Ax根據(jù)下列方程更新：

式中：αi是沿著pi方向進行第i次迭代搜索時的步長；pi為第i次搜索方向。

選擇搜索方向pi，使得函數(shù)f （x）對于沿著搜索方向的距離保持不變。當rTi＋1pi＝0，給定，可以得到搜索方向的更新方程：

式中：βi是共軛梯度方向更新因子。

選擇βi使得Api＝0（搜索方向是線性獨立且與A是共軛的），同時也得到：

通過選擇p0＝r0，最小值出現(xiàn)在由 ｛r0，Ar0，A2r0，… ｝分隔的Krylov子空間中。

3.2 最小殘差法

而最小殘差法是選擇H＝ATA和h＝ATb來最小化函數(shù)f （x）。在該方法中，對于x、r、p的更新與共軛梯度法一致，α和β需要重新定義，以使得Api＝0和 （Api＋1）TApi＝0，這也就給出了：

最小殘差法與共軛梯度法一樣，也是選擇p0＝r0，使最小值出現(xiàn)在由 ｛r0，Ar0，A2r0，… ｝ 分隔的Krylov子空間中。

4 模型研究

為了比較共軛梯度法和最小殘差法2種求解算法在構(gòu)建低感應數(shù)預條件矩陣后對正演結(jié)果的影響，筆者設計了一種理論模型進行試算。設計的三維棱柱體模型如圖1所示，幾何尺寸為4km×4km×0．1km，電阻率為1Ω·m，頂面埋深為1000m，圍巖電阻率為100Ω·m。

圖1 三維棱拄體模型

研究了大地電磁頻率為10Hz時 （因大地電磁的TE（正交極化或水平極化）、TM （平行極化或垂直極化）模式對該次研究沒有影響，故未區(qū)分），在QMR迭代計算地表電場的過程中，分別采用最小殘差法和共軛梯度法獲取的低感應數(shù)預條件矩陣作為預條件，研究相對誤差和時間隨每一次迭代過程的變化關(guān)系。從圖2、3可以看出，在每一次線性方程組的迭代求解過程中，采用最小殘差法獲得的相對誤差及需要的時間都要小于共軛梯度法，并且線性方程組的迭代次數(shù)比共軛梯度法所需要的62次少1次。究其原因在于每一次迭代過程中需要通過先后計算矢量勢 （無散場部分F）和標量勢 （無旋場部分▽f），從而來構(gòu)建低感應數(shù)預條件矩陣，而利用最小殘差法需要的迭代次數(shù)均要小于共軛梯度法，如圖4、5所示。

圖2 2種方法求低感應數(shù)預條件矩陣時相對誤差隨QMR迭代次數(shù)變化的曲線對比

圖3 2種方法求低感應數(shù)預條件矩陣時求解時間隨QMR迭代次數(shù)變化的曲線對比

圖4 2種方法分別求標量勢時的迭代次數(shù)與QMR迭代次數(shù)的關(guān)系

圖5 2種方法分別求矢量勢時的迭代次數(shù)與QMR迭代次數(shù)的關(guān)系

5 結(jié)語

在對研究區(qū)域經(jīng)過長方體網(wǎng)格的頻域有限差分離散后，三維大地電磁的正演計算，時間主要耗費在線性方程組的迭代求解過程中，以此來達到所要求的計算精度。為了加快迭代收斂的速度，必須使用預條件技術(shù)。當場源頻率低到接近似穩(wěn)場時，低感應數(shù)預條件矩陣具有明顯的加速效果。采用最小殘差法構(gòu)建低感應數(shù)預條件矩陣后，線性方程組的迭代次數(shù)和和求解時間都較同一類型的共軛梯度法要少，具有一定的速度優(yōu)勢；同樣，在構(gòu)建低感應數(shù)預條件矩陣的過程中，最小殘差法需要先行計算的標量勢和矢量勢的迭代次數(shù)均比采用共軛梯度法所需的迭代次數(shù)要少。

［1］Siripunvaraporn W，Egbert G，Lencury Y．Three-dimensional magnetotelluric inversion：data-space method ［J］ ．Physics of the Earth and Planetary Interiors，2005，150：3～14．

［2］Siripunvaraporn W，Sarakorn W．An efficient data space conjugate gradient Occam’s method for three-dimensional magnetotelluric inversion ［J］．Geophys J Int，2011，186：567～579．

［3］Newman G A，Alumbaugh D L．Frequency-domain modeling of airborne electromagnetic responces using staggered finite differences［J］．Geophys Prosp，1995，43：1021～1042．

［4］Newman G A，Alumbaugh D L．Three-dimensional induction logging problems，Part 2：A finite-difference solution ［J］ ．Geophysics，2002，67：484～491．

［5］Mackie R L，Smith J T，Madden T R．Three-dimensional electromagnetic modeling using finite difference equations：The magnetotelluric example［J］．Radio Sci，1994，29：923～935．

［6］譚捍東，余欽犯，Booker John，等 ．大地電磁法三維交錯采樣有限差分數(shù)值模擬 ［J］．地球物理學報，2003，46（5）：705～711．

［7］肖騎彬，趙國澤 ．大地電磁有限差分數(shù)值解對比 ［J］．地球物理學報，2010，53（3）：622～630．

［8］李焱，胡祥云，楊文采，等 ．大地電磁三維交錯網(wǎng)格有限差分數(shù)值模擬的并行計算研究 ［J］．地球物理學報，2012，55（12）：4036～4043．

［9］Yee K S．Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media ［J］．IEEE Trans Ant Prop，1996，14：302～307．

［10］Druskin V L，Knizhnerman L A，Lee P．New spectral Lanczos decompostion method for induction modeling in arbitrary 3-D geometry［J］．Geophysics，1999，64：701～706．

［11］Weiss C J，Newman G A．Electromagnetic induction in a generalized 3Danisotropic earth，Part 2：The LIN preconditioner ［J］．Geophysics，2003，68：922～930．