王威，李景富，劉潔

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1.3 預(yù)條件共軛梯度算法

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

對(duì)于方程組：

作變換：

如果矩陣C滿足：

那么易知：

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

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

2 模型分析

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

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

圖4為分別采用Hex和Tet網(wǎng)格求解的視電阻率圖。方框和十字分別代表采用Hex網(wǎng)格的ICCG解和SORCG解，菱形和圓圈分別代表采用Tet網(wǎng)格的ICCG解和SORCG解，黑色實(shí)線表示解析解?？梢钥闯霾还懿捎媚姆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是收斂分析。橫坐標(biāo)為迭代次數(shù)，縱坐標(biāo)為相對(duì)殘差。可以看出，達(dá)到預(yù)設(shè)精度時(shí)SORCG的收斂次數(shù)要少于ICCG，效率相對(duì)較高，但I(xiàn)CCG的迭代次數(shù)也沒(méi)有超過(guò)200次。需要說(shuō)明的是，SORCG方法需要測(cè)試松弛因子ω，具有一定的經(jīng)驗(yàn)性[7]，而ICCG則無(wú)需人工干預(yù)其求解過(guò)程。總的來(lái)看，兩者在采用均勻網(wǎng)格剖分時(shí)都是較優(yōu)的求解方法。

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

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

圖6 采用4種剖分密度的均勻Hex、Tet網(wǎng)格時(shí)ICCG和SORCG的求解時(shí)間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)格的不均勻剖分會(huì)嚴(yán)重影響有限元系數(shù)矩陣的性質(zhì)，有可能引起迭代解法的失敗。同樣基于上述物理模型參數(shù)，此處采用不均勻的網(wǎng)格剖分，六面體單元數(shù)為54×54×29，在此基礎(chǔ)上進(jìn)行二次剖分為四面體，注意兩者的單元結(jié)點(diǎn)數(shù)一致，均為55×55×30。具體剖分如下：x，y方向的剖分坐標(biāo)均為（-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方向（指向地下空間為正方向）坐標(biāo)為（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。

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

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

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

圖8 采用不均勻網(wǎng)格時(shí)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）對(duì)比發(fā)現(xiàn)，采用均勻網(wǎng)格時(shí)的迭代解法收斂率普遍較快（200次左右），而采用非均勻網(wǎng)格時(shí)，收斂率較慢（達(dá)500此左右），說(shuō)明網(wǎng)格的均勻性極大影響了收斂速率；

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

3 結(jié) 論

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

本文提出新的網(wǎng)格剖分方案，即在六面體網(wǎng)格基礎(chǔ)上進(jìn)行二次剖分得到四面體網(wǎng)格。模型分析發(fā)現(xiàn)這種四面體網(wǎng)格即使剖分不均勻，也可滿足ICCG構(gòu)建預(yù)條件因子的要求，從而成功求解。不僅如此，這種新的四面體網(wǎng)格相對(duì)于六面體網(wǎng)格還可以減少約50%的內(nèi)存空間，降低對(duì)計(jì)算設(shè)備的要求。SORCG也可高效求解直流電阻率有限元方程，但在構(gòu)建預(yù)條件因子時(shí)需要額外測(cè)試選擇松弛因子ω，需要一定的經(jīng)驗(yàn)性[7]。由于三維直流電阻率模擬的區(qū)域范圍較大，非均勻網(wǎng)格剖分具有更大的現(xiàn)實(shí)意義。本文提出的基于四面體網(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]翁愛(ài)華，祝嗣安.天然氣水合物的近海底直流電測(cè)深響應(yīng)[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]強(qiáng)建科，阮百堯，周俊杰，等.煤礦巷道直流三極法超前探測(cè)的可行性[J].地球物理學(xué)進(jìn)展，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]劉洋，吳小平.巷道超前探測(cè)的并行Monte Carlo方法及電阻率各向異性影響[J].地球物理學(xué)報(bào)，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].地球物理學(xué)進(jìn)展，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].地球物理學(xué)報(bào)，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]張錢(qián)江，戴世坤，陳龍偉，等.多源條件下直流電阻率法有限元三維數(shù)值模擬中一種近似邊界條件[J].地球物理學(xué)學(xué)報(bào)，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]李勇，林品榮，徐寶利，等.復(fù)雜地形三維直流電阻率有限元數(shù)值模擬[J].地球物理學(xué)進(jìn)展，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.