基于預(yù)處理共軛梯度的大地電磁快速正演

張繼鋒 ，湯井田，王燁，肖曉

(1. 長(zhǎng)安大學(xué) 地質(zhì)工程與測(cè)繪學(xué)院，陜西 西安，710054；2. 中南大學(xué) 信息物理工程學(xué)院，湖南 長(zhǎng)沙，410083)

大地電磁方法因具有成本低、工作方便、不受高阻層屏蔽等優(yōu)點(diǎn)，已廣泛應(yīng)用于深部地球物理勘探以及石油天然氣等能源的勘查。由于有限元方法具有模擬復(fù)雜地形的能力，適應(yīng)性比較強(qiáng)，其應(yīng)用越來(lái)越廣泛。Coggon[1]利用能量最小化原理推導(dǎo)了電磁變分方程，采用低階次三角形單元對(duì)半空間下二維地電結(jié)構(gòu)進(jìn)行了有限元數(shù)值模擬；Wannamaker等[2]研究了二維帶地形的大地電磁正演問(wèn)題；Key等[3]采用非結(jié)構(gòu)化網(wǎng)格的自適應(yīng)有限元對(duì)二維地電模型進(jìn)行了模擬；陳小斌等[4-9]采用有限元法研究了大地電磁的二維正演問(wèn)題。這些研究主要從方法或計(jì)算精度上進(jìn)行了探討，對(duì)于如何壓縮存儲(chǔ)和快速求解等問(wèn)題，張永杰等[10-14]研究了系數(shù)矩陣的存貯和ICCG方法(不完全Cholesky分解預(yù)處理共軛梯度法)的求解，但僅限于實(shí)數(shù)范圍；Smith等[15-16]研究了雙共軛梯度法的求解，謝德馨等[17-18]研究了電磁場(chǎng)分析中大型復(fù)系數(shù)方程組的求解，但大都是針對(duì)非地球物理電磁場(chǎng)進(jìn)行研究。在此，本文作者從地球物理電磁場(chǎng)中大地電磁二維正演出發(fā)，采用全稀疏按行壓縮存儲(chǔ)方式，把不完全LDLT(改進(jìn)的 Cholesky分解)應(yīng)用于矩陣的預(yù)處理，結(jié)合共軛梯度迭代法快速求解大型稀疏復(fù)系數(shù)方程組，以實(shí)現(xiàn)二維大地電磁的快速正演模擬。

1 二維大地電磁有限元正演

假定地下電性結(jié)構(gòu)是二維的，取走向?yàn)?z軸，x軸與 z軸垂直，保持水平，y軸垂直向上。當(dāng)平面電磁波以任意角度入射地面時(shí)，地下介質(zhì)中的電磁波總以平面波形式幾乎垂直向下傳播。將麥克斯韋方程組中的2個(gè)旋度方程按分量展開，并考慮到0/=??z，得2個(gè)獨(dú)立的方程組，分別命名為TE模式和TM模式，可統(tǒng)一表示為式中：Ez為電場(chǎng)的z軸分量；Hz為磁場(chǎng)的z軸分量；σ為地下介質(zhì)的電導(dǎo)率；ω為角頻率；μ為大地的磁導(dǎo)率；ε為介質(zhì)的介電常數(shù)。式(4)即為有限元正演所依賴的變分公式。其中：是上邊界；是下邊界；側(cè)邊界為齊次邊界條件，不予以考慮。

2 壓縮存儲(chǔ)

傳統(tǒng)的一維變帶寬或二維定帶寬存儲(chǔ)方式只存儲(chǔ)帶寬內(nèi)的所有元素，節(jié)省了許多存儲(chǔ)空間，但即使在帶寬內(nèi)，也包含著大量的元素 0，若節(jié)點(diǎn)的編號(hào)很不規(guī)則，則會(huì)形成1個(gè)非帶狀的矩陣，不能采用該存儲(chǔ)方式。按行壓縮存儲(chǔ)方案，根據(jù)有限元?jiǎng)偠染仃嚧笮拖∈璧奶匦裕瑑H存儲(chǔ)剛度矩陣中的非零元素，不需要考慮帶寬和網(wǎng)格的劃分規(guī)則，大大節(jié)省了存儲(chǔ)空間，而且在使用迭代法求解時(shí)，元素0不參與運(yùn)算，節(jié)省了計(jì)算時(shí)間，特別適合于各種高效快速的迭代方法求解。

按行壓縮是把剛度矩陣中的非零元素按行依次存儲(chǔ)在1個(gè)一維數(shù)組中，然后，通過(guò)2個(gè)索引數(shù)組以確定該一維數(shù)組中的元素與原矩陣中元素的對(duì)應(yīng)關(guān)系。假設(shè)有1個(gè)稀疏系數(shù)矩陣M(如式(5)所示)，其中大部分元素為0，通過(guò)3個(gè)數(shù)組a，ia和ja存儲(chǔ)該矩陣中的非零元素。首先，把該矩陣中的非零元素按行依次存儲(chǔ)在數(shù)組a中；然后，通過(guò)索引數(shù)組ia和ja確定數(shù)組a中的元素在矩陣M中具體位置： ja表示列索引，即a(k)=M(i, j)，那么 ja(k)=j；ia表示行索引，若 M(i, j)是第i行第一個(gè)非零元素，并且a(k) = M(i, j)，則ia(i)= k。這樣，數(shù)組a就和矩陣M建立了一一對(duì)應(yīng)關(guān)系，為各個(gè)元素參與共軛梯度法迭代求解提供了基礎(chǔ)。稀疏矩陣按行壓縮存儲(chǔ)表見表1。

表1 稀疏矩陣按行壓縮存儲(chǔ)表Table 1 Row compressed store for sparse matrix

上例直觀地給出了矩陣的元素，可以很清楚地看到哪些是非零元素。但是，實(shí)際有限元程序的編制中直觀的剛度矩陣并不存在，所以，必須了解非零元素在剛度矩陣中的位置，以便真正實(shí)現(xiàn)按行壓縮存儲(chǔ)方案，其基本步驟歸納如下：

表2 二維不同存儲(chǔ)方法的存儲(chǔ)量比較Table 2 Comparison of memory capacitance for two dimensional different store methods

(1) 尋找包含第 i (i=1, 2, …, m)個(gè)節(jié)點(diǎn)的所有單元。

(2) 找出與節(jié)點(diǎn)i相關(guān)的所有節(jié)點(diǎn)，然后對(duì)其進(jìn)行排序，相同的節(jié)點(diǎn)只出現(xiàn)1次。

(3) 確定節(jié)點(diǎn)的稀疏關(guān)系。

(4) 把所有單元?jiǎng)偠染仃嚱M裝加到全局剛度矩陣中。

為了更直觀地看到壓縮存儲(chǔ)對(duì)于節(jié)省內(nèi)存空間的優(yōu)越性，將一維變帶寬存儲(chǔ)對(duì)規(guī)則的四邊形網(wǎng)格劃分形成的系數(shù)矩陣的存儲(chǔ)量進(jìn)行比較，如表2所示。

從表2可以看到：隨著網(wǎng)格節(jié)點(diǎn)數(shù)增多，壓縮比例越來(lái)越大，也就是說(shuō)，剛度系數(shù)矩陣越大，越利于壓縮存儲(chǔ)。所以，在進(jìn)行大型的有限元計(jì)算時(shí)，必須進(jìn)行壓縮存儲(chǔ)，才能夠有效地利用存儲(chǔ)空間；其次，按行壓縮的存儲(chǔ)量隨著網(wǎng)格節(jié)點(diǎn)的增加，其壓縮比例急劇增大，遠(yuǎn)遠(yuǎn)優(yōu)于一維變帶寬存儲(chǔ)。這是因?yàn)榘葱袎嚎s只存儲(chǔ)非零元素，方程組的系數(shù)矩陣不斷增大，而每一行的非零元素變化不大，所以，它的壓縮量最大，充分顯示了有限元系數(shù)矩陣的稀疏性。

3 預(yù)處理共軛梯度

3.1 不完全LDLT預(yù)處理

有限元方法中最后形成的大型系數(shù)矩陣常常具有很大的條件數(shù)，是一個(gè)非常病態(tài)的矩陣。對(duì)于這種矩陣，若直接采用共軛梯度法，則收斂速度非常慢，甚至不收斂，所以，必須采取合適的預(yù)處理方法，通過(guò)改變矩陣的一些關(guān)鍵元素值，降低矩陣的條件數(shù)，使共軛梯度法能夠順利迭代，提高其收斂速度。這里采用不完全 LDLT預(yù)處理，避免開方運(yùn)算，可用于非正定方程組的求解。

首先，把LDLT具體的算法簡(jiǎn)要陳述為：

其中：dii為分解后對(duì)角陣元素；lij為分解后下三角陣元素；aij為原始矩陣元素。

但是，在LDLT預(yù)處理過(guò)程中，產(chǎn)生的L矩陣是1個(gè)稠密矩陣，原來(lái)的系數(shù)矩陣稀疏性被破壞，這與壓縮存儲(chǔ)的思想不符，所以，必須采用不完全 LDLT分解法，使存儲(chǔ)量和改善系數(shù)矩陣性態(tài)之間達(dá)到平衡。

3.2 基于LDLT預(yù)處理共軛梯度法

給定1個(gè)初始猜想的x0，有

對(duì)i=0, 1, 2, …，進(jìn)行如下計(jì)算：

時(shí)迭代停止。其中：iα，1+ix，1+ir，iβ和1+ip表示第i次迭代的值；ε為迭代誤差的判斷值；〉〈βα,表示α與β的內(nèi)積。

3.3 逆矩陣與向量乘積的快速計(jì)算

以上算法式中，存在(LDLT)-1與任一向量的乘積，運(yùn)算時(shí)間長(zhǎng)。如何快速高效地求解(LDLT)-1r是預(yù)處理共軛梯度法速度的關(guān)鍵。為了解決此問(wèn)題，設(shè)

可見：把(LDLT)-1r轉(zhuǎn)換為求解2個(gè)線性方程組，由于L是一個(gè)下三角陣，且僅存非零元素，因而，以上順代和回代僅涉及少量的非零元素的乘積和加法運(yùn)算，求解非常迅速。

3.4 數(shù)值算例

針對(duì)預(yù)處理共軛梯度算法，編制相應(yīng)的程序，并且在計(jì)算機(jī)上進(jìn)行驗(yàn)證。計(jì)算機(jī)配置內(nèi)存為2G，主頻為2.26 Hz，硬盤容量150 Gbyte。圖1和圖2所示是不同自由度下共軛梯度法的收斂速度，可以看出：ILDLT預(yù)處理共軛梯度法收斂速度非?？?，但在收斂過(guò)程中，迭代誤差會(huì)上下波動(dòng)，最終收斂到準(zhǔn)確解。計(jì)算時(shí)間和迭代次數(shù)見表3。從表3可以看出：采用該方法計(jì)算上萬(wàn)個(gè)自由度的大型復(fù)系數(shù)方程組，時(shí)間在1 s之內(nèi)，這表明本文算法非常高效。

圖1 2 925個(gè)自由度的ILDLTCG收斂率Fig.1 ILDLTCG convergence rate of 2 925 degree of freedom

圖2 13 020個(gè)自由度的ILDLTCG收斂率Fig.2 ILDLTCG convergence rate of 13 020 degree of freedom

表3 自由度不同時(shí)的計(jì)算時(shí)間及迭代次數(shù)Table 3 Computation time and iteration number for different degrees of freedom

4 地電模型正演

為了檢驗(yàn)程序的正確性，建立了1個(gè)層狀模型，分別對(duì)二維大地電磁的2種極化模式進(jìn)行計(jì)算，并與解析解進(jìn)行對(duì)比。模型層電阻率分別為：ρ1=10 Ω·m，ρ2=100 Ω·m，ρ3=10 Ω·m，ρ4=100 Ω·m，ρ5=10 Ω·m；層厚度分別為：h1=300 m，h2=700 m，h3=1 500 m，h4=3 100 m。采用如下40個(gè)頻點(diǎn)：320，240，160，120，80，60，40，30，20，15，10，7.5，6.0，4.5，3.0，2.25，1.5，1.125，0.75，0.562 5，0.375，0.281 25，0.187 5，0.140 625，0.093 75，0.070 312 5，0.046 875，0.035 156 25，0.023 437 5，0.017 578 125，0.011 718 75，0.008 789 062 5，0.005 859 375，0.004 394 531 25，0.002 929 687 5， 0.002 197 265 625，0.001 464 843 5，0.001 098 632 5，0.000 732 421 75，0.000 549 316 25。地電模型的視電阻率曲線和相應(yīng)曲線分別見圖 3和圖4。

圖3 地電模型的視電阻率曲線Fig.3 Curve of apparent resistivity for geo-electrical model

圖4 地電模型的視相位曲線Fig.4 Curve of apparent phase for geo-electrical model

經(jīng)計(jì)算，對(duì)于TE模式和TM模式，視電阻率ρs平均相對(duì)誤差分別為 0.384%和 0.374%，最大相對(duì)誤差分別為 1.299%和 2.029%；視相位φ平均相對(duì)誤差分別為0.108%和0.122%，最大相對(duì)誤差分別為0.320%和0.663%，而且最大相對(duì)誤差都出現(xiàn)在極低頻段。

5 結(jié)論

(1) 按行壓縮存儲(chǔ)1萬(wàn)個(gè)自由度矩陣時(shí)，壓縮率達(dá)到 99.9%以上，比常規(guī)的變帶寬存儲(chǔ)方法節(jié)省了大量?jī)?nèi)存空間，而且矩陣越大，壓縮比例越大，尤其有利于存儲(chǔ)超大型稀疏矩陣。

(2) 不完全 LDLT預(yù)處理極大地提高了共軛梯度法的收斂速度，使得求解1萬(wàn)個(gè)自由度的方程組時(shí)所用時(shí)間在1 s以內(nèi)。

(3) 建立了大地電磁層狀模型，分別對(duì)TE模式和TM模式2種極化方式進(jìn)行了數(shù)值模擬，并與解析解進(jìn)行比較，結(jié)果表明：采用大地電磁層狀模型計(jì)算精度高，平均相對(duì)誤差在0.5%以內(nèi)，為快速精確反演提供了保證。

[1] Coggon J H. Electromagneic and electrical modeling by the finite element method[J]. Geophysics, 1971, 36(1): 132-155.

[2] Wannamaker P E, Stodt J A, Rijo L. A stable finite element solution for two-dimensional magnetotelluric modeling[J].Geophysical Journal International, 1987, 88(1): 277-296.

[3] Key K, Weiss C. Adaptive finite-element modeling using unstructured grids: The 2D magnetotelluric example[J].Geophysics, 2006, 71(6): 291-299.

[4] 陳小斌. 有限元直接迭代算法[J]. 物探化探計(jì)算技術(shù), 1999,21(2): 165-171.CHEN Xiao-bin. Direct iterative algorithm of finite element[J].Computing Techniques for Geophysical and Geochemical Exploration, 1999, 21(2): 165-171.

[5] 徐世浙. 地球物理中的有限單元法[M]. 北京: 科學(xué)出版社,1994: 229-247.XU Shi-zhe. Finite element method in geophysics[M]. Beijing:Science Press, 1994: 229-247.

[6] 陳小斌, 張翔, 胡文寶. 有限元直接迭代算法在 MT二維正演計(jì)算中的應(yīng)用[J]. 石油地球物理勘探, 2000, 35(4): 487-496.CHEN Xiao-bin, ZHANG Xiang, HU Wen-bao. Application for direct iterative algorithm in MT 2D forward computation[J]. Oil Geophysical Prospecting, 2000, 35(4): 487-496.

[7] 史明娟, 徐世浙, 劉斌. 大地電磁二次函數(shù)插值的有限元正演模擬[J]. 地球物理學(xué)報(bào), 1997, 40(3): 421-430.SHI Ming-juan, XU Shi-zhe, LIU Bin. Finite element method using quadratic element in MT forward modeling[J]. Chinese J Geophys, 1997, 40(3): 421-430.

[8] 馬為, 陳小斌, 趙國(guó)澤. 大地電磁測(cè)深二維正演中輔助場(chǎng)的新算法[J]. 地震地質(zhì), 2008, 30(2): 525-533.MA Wei, CHEN Xiao-bin, ZHAO Guo-ze. A new algorithm for the calculation of auxiliary field in MT 2D forward modeling[J].Seismology and Geology, 2008, 30(2): 525-533.

[9] 劉小軍, 王家林, 于鵬. 基于二次場(chǎng)的二維大地電磁有限元數(shù)值模擬[J]. 同濟(jì)大學(xué)學(xué)報(bào), 2007, 35(8): 1113-1117.LIU Xiao-jun, WANG Jia-lin, YU Peng. Secondary field based on two dimensional magnetotelluric numerical modeling by finite element method[J]. Journal of Tongji University, 2007,35(8): 1113-1117.

[10] 張永杰, 孫秦, 李江海. 大型稀疏線性方程組的改進(jìn)ICCG方法[J]. 計(jì)算物理, 2007, 24(5): 581-584.ZHANG Yong-jie, SUN Qin, LI Jiang-hai. An improved ICCG method for large scale sparse linear equations[J]. Chinese Journal of Computational Physics, 2007, 24(5): 581-584.

[11] 張永杰, 孫秦. 大型稀疏線性方程組新的ICCG方法[J]. 數(shù)值計(jì)算與計(jì)算機(jī)應(yīng)用, 2007, 43(36): 19-20.ZHANG Yong-jie, SUN Qin. Preconditioned biconjugate gradient method of large-scale complex linear equations[J].Computer Engineering and Applications, 2007, 43(36): 19-20.

[12] 張永杰, 孫秦. 大型稀疏線性方程組的全稀疏存貯策略[J].陜西理工學(xué)院學(xué)報(bào): 自然科版, 2005, 21(4): 67-68.ZHANG Yong-jie, SUN Qin. Fully sparse strategy of large scale linear equations[J]. Journal of Shaanxi University of Technology:Natural Science Edition, 2005, 21(4): 67-68.

[13] 吳小平, 徐果明, 李時(shí)燦. 解大型稀疏方程組的ICCG方法及其計(jì)算機(jī)實(shí)現(xiàn)[J]. 煤田地質(zhì)與勘探, 1999, 27(6): 54-56.WU Xiao-ping, XU Guo-ming, LI Shi-can. ICCG Method for solving large sparse equations and its computer programming[J].Coal Geology & Exploration, 1999, 27(6): 54-56.

[14] 陳志, 高旅端. 求解大規(guī)模稀疏線性方程組的算法[J]. 北京工業(yè)大學(xué)學(xué)報(bào), 2001, 27(3): 262-265.CHEN Zhi, GAO Lü-duan. An algorithm for solving large-scale sparse group of linear equations[J]. Journal of Beijing Polytechnic University, 2001, 27(3): 262-265.

[15] Smith J T. Conservative modeling of 3-D electromagnetic fields,Part II bi-conjugate gradient solution and an accelerator[J].Geophysics, 1996, 61(5): 1319-1324.

[16] 吳海容. 復(fù)雙共軛梯度法的結(jié)構(gòu)[J]. 電機(jī)與控制學(xué)報(bào), 1996,19(2): 133-141.WU Hai-rong. Structure of the complex bi-conjugate gradient method[J]. Electric Machines and Control, 1996, 19(2):133-141.

[17] 謝德馨, 姚纓英, 白保東. 電磁場(chǎng)分析中大型稀疏對(duì)稱線性方程組預(yù)處理法的改進(jìn)[J]. 電機(jī)與控制學(xué)報(bào), 1997, 1(2):98-101.XIE De-xin,YAO Ying-ying, BAI Bao-dong. The modified preconditioning approach for large sparse symmetric linear systems in electromagnetic field analysis[J]. Electric Machines and Control, 1997, 1(2): 98-101.

[18] 文代剛, 姜可薰, 黃鍵. 求解有限元復(fù)代數(shù)方程組的實(shí)型ICCG法[J]. 電工技術(shù)學(xué)報(bào), 1996, 11(4): 62-65.WEN Dai-gang, JIANG Ke-xun, HUANG Jian. Real type of ICCG method for complex algebraic equations of FEM[J].Transaction of China Electrotechnical Society, 1996, 11(4):62-65.