預條件共軛梯度法求解三維地電場有限元方程的網(wǎng)格分析*

王威，李景富，劉潔

（1.中山大學地球科學與工程學院∥廣東省地質(zhì)過程與礦產(chǎn)資源探查重點實驗室，廣東廣州510275；2.廣東省有色地質(zhì)環(huán)境中心，廣東廣州510060）

直流電阻率法作為經(jīng)濟高效的地球物理勘探方法，在礦產(chǎn)資源勘探，水文地質(zhì)調(diào)查，工程勘察領(lǐng)域獲得了廣泛的應用[1-6]。正演數(shù)值模擬不僅是分析特定地電結(jié)構(gòu)響應特征的有效手段，而且作為反演的必要組成步驟，其求解效率直接關(guān)系到反演是否具有實用性。常用的數(shù)值模擬方法有積分方程法、有限差分法和有限元法，其中有限元法因其堅實的理論背景和靈活的模擬策略，已成為電法勘探二、三維數(shù)值模擬中廣泛采用的方法[7-11]。有限元法模擬直流電阻率問題時需要求解大型線性方程組。直接解法由于矩陣分解需要占用大量的內(nèi)存，一般的個人計算機難以承受。迭代算法中的預條件共軛梯度（Pre-conditioned Conjugate Gradient）法無需分解系數(shù)矩陣，可有效減少內(nèi)存需求和提升計算速度，主要有不完全Cholesky共軛梯度（ICCG）算法和超松弛預條件共軛梯度（SORCG）算法。

網(wǎng)格剖分是有限元模擬的重要課題。規(guī)則六面體網(wǎng)格由于剖分簡單而在直流電法有限元模擬中被廣泛采用。以往研究認為ICCG算法只適合求解三維有限差分法形成的系統(tǒng)方程，而不適合于三維有限元方程的求解[12]。然而，這一結(jié)論的給出是基于六面體網(wǎng)格的。本文提出了一種新的四面體剖分方案，經(jīng)過系統(tǒng)的網(wǎng)格對比分析，結(jié)果表明我們提出的網(wǎng)格剖分方案不僅可以使ICCG法穩(wěn)定求解有限元方程，而且存儲量只需采用常規(guī)六面體網(wǎng)格時的50%。本研究推動了三維直流電阻率法有限元模擬的發(fā)展，為其實際應用奠定了基礎(chǔ)。

1 直流電阻率法的有限元求解

1.1 有限元方程的構(gòu)建

關(guān)于二次場電位的點源三維地電場邊值問題可表示為[13]：

其中σ為電導率，它與電阻率ρ互為倒數(shù)關(guān)系；Up是由背景電導率σp引起的背景電位，Us由電導率異常體σs引起的二次場電位；r是源點到測點的向量；r是源點到測點的距離；n是求解域邊界的法方向；Γs、?！薹謩e是地表和計算區(qū)域剩余部分的邊界。

根據(jù)變分原理，該邊值問題可轉(zhuǎn)換為求解泛函的極值問題：

將求解域離散為若干單元，假設(shè)每個單元內(nèi)電位符合線性分布。令泛函變分為零，可得有限元方程：（具體過程可參見文獻[7]）

其中有限元系數(shù)K一個稀疏對稱的矩陣，與地下結(jié)構(gòu)的電阻率分布有關(guān)，右端項F是與激發(fā)源分布有關(guān)的向量，Us為需要求解的各個網(wǎng)格結(jié)點上的二次電位值。解出Us之后加上解析方法求出的一次場電位Up即可得最終需要的總場電位U。

1.2 有限元網(wǎng)格分析

網(wǎng)格剖分是關(guān)系預條件共軛梯度法能否穩(wěn)定求解的重要因素。規(guī)則六面體網(wǎng)格（以下簡稱Hex）由于其剖分簡單而獲得廣泛應用。規(guī)則六面體網(wǎng)格的內(nèi)部結(jié)點只與26個結(jié)點相鄰（圖1）。由有限元理論可知最終形成的系數(shù)矩陣K每行的非零元素最多只有27個，又由于K具有對稱性，故每行只需存儲14個非零元素。

圖1 六面體網(wǎng)格的結(jié)點連接關(guān)系Fig.1 The nodes'connection of hexahedral mesh

本文提出在六面體網(wǎng)格基礎(chǔ)上進行二次剖分得到四面體網(wǎng)格。具體方案為將每個六面體單元剖分為5個四面體單元，注意對鄰接六面體單元進行的是兩種不同的劃分（圖2）。圖2（a）中單元的結(jié)點編號分別為（1，5，4，2）、（8，4，5，7）、（3，7，2，4）、（6，2，7，5）、（4，5，7，2），圖2（b）中單元的結(jié)點編號分別為（1，6，3，2）、（8，3，6，7）、（4，1，8，3）、（5，8，1，6）、（1，3，6，8）。如圖 3所示，這種剖分方法得到的四面體網(wǎng)格一半的內(nèi)部結(jié)點具有18個鄰接結(jié)點，另一半的內(nèi)部結(jié)點具有6個鄰接結(jié)點。同樣，由有限元理論可知，系數(shù)矩陣平均每行只有7個非零元素需要存儲（由于系數(shù)矩陣的對稱性，只需存儲下三角矩陣），故其存儲需求約為采用六面體網(wǎng)格的50%。

圖2 在六面體網(wǎng)格基礎(chǔ)上二次剖分得到的四面體網(wǎng)格Fig.2 The tetrahedral mesh derived from hexahedral mesh

圖3 四面體網(wǎng)格的結(jié)點連接關(guān)系Fig.3 The nodes'relationship of tetrahedral mesh

1.3 預條件共軛梯度算法

直流電阻率法的有限元模擬需要求解大型線性方程組Ax＝b。預條件共軛梯度法是求解該問題的高效算法，其基本思想如下：

對于方程組：

作變換：

如果矩陣C滿足：

那么易知：

即變換后的線性系統(tǒng)（5）的系數(shù)矩陣近似于單位矩陣I，求解將快速收斂。定義M＝CCT，M即稱為預條件矩陣。問題轉(zhuǎn)化為如何尋找到合適的預條件矩陣M。ICCG和SORCG的預條件矩陣的獲取可參考文獻[7，12]，在此不再贅述。

在以下的模型計算中，我們采用殘差的L2范數(shù)作為預條件共軛梯度法的收斂準則，即＜10-8，其中x0為初始解，xi為第i步的解。

2 模型分析

2.1 均勻網(wǎng)格剖分

層狀介質(zhì)在地球物理勘探中是十分常見的。此處計算一個兩層的地電模型，第一層電阻率為100 Ω·m，厚度為12 m，第二層為半無限空間，電阻率為10Ω·m。實際模擬計算都是在有限的區(qū)域中進行，設(shè)定計算區(qū)域在x，y，z方向上的尺度分別為（-204～204）m，（-204～204）m，（0～204）m。點電源在坐標原點激發(fā)，測量裝置為單極-單極（pole-pole）。將計算區(qū)域剖分為102×102×51個單元的均勻六面體（Hex）網(wǎng)格（單元邊長為4 m），并在此基礎(chǔ)上進行二次剖分，將單元劃分為四面體（Tet）網(wǎng)格。注意這兩種剖分方法并不會帶來單元結(jié)點數(shù)的差異，即兩者結(jié)點數(shù)均為103×103×52。以下分別采用ICCG和SORCG求解，來分析它們的求解精度和收斂效率。

圖4為分別采用Hex和Tet網(wǎng)格求解的視電阻率圖。方框和十字分別代表采用Hex網(wǎng)格的ICCG解和SORCG解，菱形和圓圈分別代表采用Tet網(wǎng)格的ICCG解和SORCG解，黑色實線表示解析解?？梢钥闯霾还懿捎媚姆N網(wǎng)格，ICCG和SORCG都可以獲得與解析解非常一致的結(jié)果，精度均符合要求。

圖4 采用均勻Hex、Tet網(wǎng)格的ICCG和SORCG解Fig.4 Solutions from ICCG and SORCG by adopting uniform hexahedral and tetrahedral meshes

圖5是收斂分析。橫坐標為迭代次數(shù)，縱坐標為相對殘差。可以看出，達到預設(shè)精度時SORCG的收斂次數(shù)要少于ICCG，效率相對較高，但ICCG的迭代次數(shù)也沒有超過200次。需要說明的是，SORCG方法需要測試松弛因子ω，具有一定的經(jīng)驗性[7]，而ICCG則無需人工干預其求解過程?？偟膩砜?，兩者在采用均勻網(wǎng)格剖分時都是較優(yōu)的求解方法。

為了詳細比較不同網(wǎng)格剖分密度時ICCG和SORCG的求解效率，我們設(shè)計了4種密度的網(wǎng)格：68×68×34，102×102×51，136×136×68，204×204×102，其單元邊長分別為6、4、3和2 m。我們統(tǒng)計了不同方案的計算時間（圖6），可以看出：不論是用ICCG還是SORCG求解，只要采用均勻網(wǎng)格，都可以獲得較高的效率。即使高達400萬的網(wǎng)格結(jié)點，在單機（DELL T7810工作站，CPU E5-2609，內(nèi)存32G）上求解時間也均不超過20 min。

圖5 采用均勻Hex、Tet網(wǎng)格時ICCG和SORCG的收斂情況Fig.5 The convergence analysis of ICCG and SORCG by employing uniform hexahedral and tetrahedral meshes

圖6 采用4種剖分密度的均勻Hex、Tet網(wǎng)格時ICCG和SORCG的求解時間Fig.6 The ICCG and SORCG solving time when employing uniform hexahedral and tetrahedral meshes of four different mesh sizes

2.2 非均勻網(wǎng)格剖分

網(wǎng)格的不均勻剖分會嚴重影響有限元系數(shù)矩陣的性質(zhì)，有可能引起迭代解法的失敗。同樣基于上述物理模型參數(shù)，此處采用不均勻的網(wǎng)格剖分，六面體單元數(shù)為54×54×29，在此基礎(chǔ)上進行二次剖分為四面體，注意兩者的單元結(jié)點數(shù)一致，均為55×55×30。具體剖分如下：x，y方向的剖分坐標均為（-500，-375，-275，-200，-150，-120，-90，-70，-55，-42.5，-30，-20，-15，-12，-9，-7，-5.5，-4.25，-3.25，-2.5，-2，-1.5，-1.2，-1，-0.8，-0.5，-0.2，0，0.2，0.5，0.8，1，1.2，1.5，2，2.5，3.25，4.25，5.5，7，9，12，15，20，30，42.5，55，70，90，120，150，200，275，375，500），z方向（指向地下空間為正方向）坐標為（0，0.2，0.5，0.75，1，1.25，1.5，2，2.5，3，4，5，7，9，12，15，20，25，32.5，42.5，55，70，90，120，150，200，250，325，400，500），單位：m。

計算結(jié)果如圖7所示。我們發(fā)現(xiàn)，正如已有的研究結(jié)果[12]，采用不均勻Hex網(wǎng)格時ICCG的確會求解失敗。究其原因，是由于ICCG在構(gòu)建預條件因子時，對角矩陣D的元素出現(xiàn)了大量負數(shù)。對于此例，統(tǒng)計發(fā)現(xiàn)13.98%的對角元素為負，直接導致了ICCG求解失敗。但是，采用本文提出的在Hex基礎(chǔ)上繼續(xù)二次剖分得到的Tet網(wǎng)格形成的有限元方程，其系數(shù)矩陣滿足ICCG的要求，即對角矩陣D的元素全部為正數(shù)，求解成功。

比較采用Hex網(wǎng)格的SORCG（圖7黑色圓圈）解、采用Tet網(wǎng)格的ICCG解（菱形）和SORCG解（十字），都與解析解（實線）吻合較好，并且符合理論預期，即小極距視電阻率趨近于第一層實際電阻率100Ω·m，大極距視電阻率趨近于第二層（半空間）的實際電阻率10Ω·m。分析不同方案的計算結(jié)果（圖8）可以得出：

圖7 采用不均勻網(wǎng)格的兩層模型的視電阻率Fig.7 The apparent resistivities of a two-layer model by employing non-uniform meshes

圖8 采用不均勻網(wǎng)格時ICCG和SORCG的收斂情況Fig.8 The convergence analysis of ICCG and SORCG by employing non-uniform meshes

（1）采用Hex不均勻網(wǎng)格的ICCG求解失敗，但是采用本文提出的Tet不均勻網(wǎng)格的ICCG求解獲得成功；

（2）與均勻網(wǎng)格結(jié)果（圖5）對比發(fā)現(xiàn)，采用均勻網(wǎng)格時的迭代解法收斂率普遍較快（200次左右），而采用非均勻網(wǎng)格時，收斂率較慢（達500此左右），說明網(wǎng)格的均勻性極大影響了收斂速率；

（3）采用Tet網(wǎng)格可以獲得和Hex網(wǎng)格相當?shù)那蠼饩群托?，但正如前文網(wǎng)格分析中所述，本文提出的Tet網(wǎng)格的存儲要求只需Hex網(wǎng)格的1/2，由此降低了對計算設(shè)備的要求。

3 結(jié) 論

ICCG是求解直流電阻率有限差分方程的主流方法，但在求解有限元方程有可能會失敗。本文通過詳細的網(wǎng)格分析后認為，基于六面體網(wǎng)格剖分的有限元系數(shù)矩陣會嚴重依賴于網(wǎng)格剖分是否均勻，即網(wǎng)格均勻，ICCG可以求解成功，但網(wǎng)格不均勻，會引起構(gòu)建預條件因子時對角矩陣元素出現(xiàn)大量負數(shù)而求解失敗。

本文提出新的網(wǎng)格剖分方案，即在六面體網(wǎng)格基礎(chǔ)上進行二次剖分得到四面體網(wǎng)格。模型分析發(fā)現(xiàn)這種四面體網(wǎng)格即使剖分不均勻，也可滿足ICCG構(gòu)建預條件因子的要求，從而成功求解。不僅如此，這種新的四面體網(wǎng)格相對于六面體網(wǎng)格還可以減少約50%的內(nèi)存空間，降低對計算設(shè)備的要求。SORCG也可高效求解直流電阻率有限元方程，但在構(gòu)建預條件因子時需要額外測試選擇松弛因子ω，需要一定的經(jīng)驗性[7]。由于三維直流電阻率模擬的區(qū)域范圍較大，非均勻網(wǎng)格剖分具有更大的現(xiàn)實意義。本文提出的基于四面體網(wǎng)格剖分的ICCG、SORCG算法可成功求解三維直流電阻率有限元方程。

[1]RUCKER D F，LOKE M H，LEVITT M T，et al.Electrical resistivity characterization of an industrial site using long electrodes[J].Geophysics，2010，75（4）：WA95-WA104.

[2]翁愛華，祝嗣安.天然氣水合物的近海底直流電測深響應[J].石油地球物理勘探，2010，45（3）：458-461.WENG A H，ZHU S A.Near sea-bottom DC sounding response for gas hydrate[J].Oil Geophysical Prospecting，2010，45（3）：458-461.

[3]強建科，阮百堯，周俊杰，等.煤礦巷道直流三極法超前探測的可行性[J].地球物理學進展，2011，26（1）：320-326.QIANG J K，RUAN B Y，ZHOU JJ，et al.The feasibility of advanced detection using DCthree-electrode method in coal-mine tunnel[J].Progress in Geophysics，2011，26（1）：320-326.

[4]LOKE M H，CHAMBERS J E，RUCKER D F，et al.Recent developments in the direct-current geoelectrical imaging method[J].Journal of Applied Geophysics，2013，95：135-156.

[5]劉洋，吳小平.巷道超前探測的并行Monte Carlo方法及電阻率各向異性影響[J].地球物理學報，2016，59（11）：4297-4309.LIU Y，WU X P.Parallel Monte Carlo method for advanced detection in tunnel incorporating anisotropic resistivity effect[J].Chinese Journal of Geophysics，2016，59（11）：4297-4309.

[6]MITCHELL M A，OLDENBURG D W.Data quality control methodologies for large，non-conventional DCresistivity datasets[J].Journal of Applied Geophysics，2016，135：163-182.

[7]王威，吳小平.電阻率任意各向異性三維有限元快速正演[J].地球物理學進展，2010，25（4）：1365-1371.WANG W，WU X P.Rapid finite element resistivity forward modeling for 3D arbitrary anisotropic structures[J].Progress in Geophysics，2010，25（4）：1365-1371.

[8]REN Z Y，TANG J T.A goal-oriented adaptive finite-element approach for multi-electrode resistivity system[J].Geophysical Journal International，2014，199（1）：136-145.

[9]吳小平，劉洋，王威.基于非結(jié)構(gòu)網(wǎng)格的電阻率三維帶地形反演[J].地球物理學報，2015，58（8）：2706-2717.WU X P，LIU Y，WANG W.3D resistivity inversion incorporating topography based on unstructured meshes[J].Chinese Journal of Geophysics，2015，58（8）：2706-2717.

[10]張錢江，戴世坤，陳龍偉，等.多源條件下直流電阻率法有限元三維數(shù)值模擬中一種近似邊界條件[J].地球物理學學報，2016，59（9）：3448-3458.ZHANG Q J，DAISK，CHEN L W，et al.An approximation boundary condition for FEM-based 3-D numerical simulation with multi-source direct current resistivity method[J].Chinese Journal of Geophysics，2016，59（9）：3448-3458.

[11]李勇，林品榮，徐寶利，等.復雜地形三維直流電阻率有限元數(shù)值模擬[J].地球物理學進展，2009，24（3）：1039-1046.LI Y，LIN PR，XU B L，et al.FEM numerical modeling of 3-D DC resistivity under complicated terrain[J].Progress in Geophysics，2009，24（3）：1039-1046.

[12]WU X P.A 3-Dfinite-element algorithm for DCresistivity modelling using the shifted incomplete Cholesky conjugate gradient method[J].Geophysical Journal International，2003，154：947-956.

[13]WANG W，SPITZER K，WU X P.Three-dimensional DC anisotropic resistivity modelling using finite elements on unstructured grids[J].Geophysical Journal International，2013，193（2）：734-746.