基于MLS形函數(shù)處理邊界的電阻率法有限單元正演模擬

麻昌英, 趙文學(xué), 湯文武1,*, 柳建新, 閆玲玲, 周聰1,,秦臻1,, 鐘煒城, 程流燕

1 東華理工大學(xué)江西省防震減災(zāi)與工程地質(zhì)災(zāi)害探測(cè)工程研究中心, 南昌 330013 2 中南大學(xué)有色金屬成礦預(yù)測(cè)與地質(zhì)環(huán)境監(jiān)測(cè)教育部重點(diǎn)實(shí)驗(yàn)室, 長(zhǎng)沙 410083 3 東華理工大學(xué)地球物理與測(cè)控技術(shù)學(xué)院, 南昌 330013

0 引言

有限單元法(Finite element method, FEM)在求解區(qū)域內(nèi)根據(jù)節(jié)點(diǎn)拓?fù)潢P(guān)系剖分單元，基于單元節(jié)點(diǎn)在單元內(nèi)進(jìn)行線性或高次插值模擬場(chǎng)函數(shù)，具有模擬精度高，計(jì)算效率高，適應(yīng)性強(qiáng)等特點(diǎn).20世紀(jì)70年代初，Coggon (1971)首次將其引入到電法勘探正演模擬領(lǐng)域，國(guó)內(nèi)學(xué)者徐世浙(1994)將有限單元法應(yīng)用于電法勘探中，對(duì)有限單元法進(jìn)行了系統(tǒng)的闡述.此后，有限單元法成為目前國(guó)內(nèi)外電阻率法正演模擬中主要的數(shù)值模擬方法(e.g., Sasaki, 1994; 強(qiáng)建科和羅延鐘, 2007; 柯敢攀和黃清華, 2009; Huang and Lin, 2010; Xiong et al., 2018；麻昌英等, 2014; Hu et al., 2021; 戴世坤等, 2022).近年來(lái)，針對(duì)復(fù)雜起伏地形和地下復(fù)雜目標(biāo)體模擬，基于非結(jié)構(gòu)化單元和自適應(yīng)加密策略的有限單元法，在電阻率法、大地電磁及可控源電磁法等正演模擬中備受關(guān)注(e.g., Bing and Greenhalgh, 2001; 任政勇和湯井田, 2009; 湯井田等, 2010; Ren and Tang, 2010; Tang et al., 2010, 2015; Vachiratienchai and Siripunvaraporn, 2013; 韓騎等, 2015; 蔡紅柱等, 2015; 吳小平等, 2015; Yan et al., 2016; Chou et al., 2016; 李勇等, 2017; 陳漢波等, 2018; 曹曉月等, 2018; Ren et al., 2018; 趙寧等, 2019； Guo et al., 2022)，借助于圖形幾何學(xué)相關(guān)單元剖分器實(shí)現(xiàn)任意節(jié)點(diǎn)分布的非結(jié)構(gòu)化網(wǎng)格剖分，可適應(yīng)起伏地形和地下目標(biāo)體的復(fù)雜形態(tài)變化，并利用后驗(yàn)誤差策略指導(dǎo)單元自適應(yīng)加密，僅需建立初始單元模型，正演過(guò)程中可在需要的局部區(qū)域自適應(yīng)加密單元提高模擬精度，為推進(jìn)高靈活性、高精度的正演模擬研究發(fā)展做出了卓越貢獻(xiàn).

有限單元法采用單元節(jié)點(diǎn)在單元內(nèi)構(gòu)造線性或高次形函數(shù)對(duì)場(chǎng)函數(shù)進(jìn)行分片插值，因此有限單元法中網(wǎng)格單元剖分不可避免，單元之間模擬的場(chǎng)函數(shù)連續(xù)性欠佳.Lancaster和Salkauskas(1981)在研究曲面和曲線的擬合時(shí)，基于傳統(tǒng)的最小二乘近似理論，提出了移動(dòng)最小二乘法(Moving Least Square, MLS)構(gòu)造形函數(shù)(MLS形函數(shù))，相比于有限單元法形函數(shù)，MLS形函數(shù)具有良好的連續(xù)性，在有限單元法后處理中也得到了應(yīng)用(Tabbara et al., 1994; Tang et al., 2021).基于MLS形函數(shù)，Belytschko等(1994)在擴(kuò)散元法及對(duì)近似導(dǎo)數(shù)確定，強(qiáng)加本質(zhì)邊界條件及數(shù)值積分處理等問(wèn)題改進(jìn)的基礎(chǔ)上提出了無(wú)單元法(Element free method, EFM).無(wú)單元法不再需要基于節(jié)點(diǎn)拓?fù)潢P(guān)系的網(wǎng)格單元剖分，自提出之后在多個(gè)正演模擬領(lǐng)域得到了廣泛應(yīng)用.在地球物理領(lǐng)域方面，無(wú)單元法中的無(wú)單元Galerkin法(Element free Galerkin method, EFM)在地震波場(chǎng)(Jia and Hu, 2006; 王月英, 2007; Jie et al., 2008)、探地雷達(dá)(戴前偉和王洪華, 2013; 馮德山等, 2013; Feng et al., 2015)和大地電磁(蘇洲等, 2012; 嚴(yán)家斌和李俊杰, 2014; 嵇艷鞠等, 2016)等正演模擬研究中已取得了良好的研究成果，展示了無(wú)單元Galerkin法在地球物理勘探領(lǐng)域的有效性.在電阻率法正演模擬中，麻昌英等(麻昌英等, 2017, 2018; Ma et al., 2019)分別利用徑向基點(diǎn)插值(RPIM)形函數(shù)和移動(dòng)最小二乘(MLS)形函數(shù)研究了2.5維電阻率法無(wú)單元Galerkin法正演，實(shí)現(xiàn)了任意復(fù)雜形體、地形起伏、局部節(jié)點(diǎn)加密的電阻率法正演模擬.以上研究成果表明，無(wú)單元Galerkin法突破了傳統(tǒng)數(shù)值模擬方法受單元約束的限制，與傳統(tǒng)數(shù)值模擬方法相輔相成，具有高靈活性、高適應(yīng)性、高精度等特點(diǎn)，特別適用于任意復(fù)雜地電模型正演，在地球物理領(lǐng)域具有良好的發(fā)展前景.然而，無(wú)單元法通常采用高斯積分求積，由于無(wú)單元法的近似函數(shù)一般為有理函數(shù)而不是多項(xiàng)式，為了保證模擬精度須使用較多的高斯點(diǎn)和較高階的高斯積分點(diǎn)進(jìn)行計(jì)算(Dolbow and Belytschko, 1999)，同時(shí)需要為每一個(gè)高斯積分點(diǎn)構(gòu)造形函數(shù)，而無(wú)單元法計(jì)算量與積分點(diǎn)數(shù)量成正比，這將導(dǎo)致無(wú)單元法計(jì)算量大，計(jì)算效率低(麻昌英等, 2017)，在實(shí)際應(yīng)用中受到了一定的限制.

在電阻率法正演模擬中，現(xiàn)有正演模擬方法中采用第一類邊界條件時(shí)，認(rèn)為截?cái)噙吔缟蠄?chǎng)值為零，需要足夠大的計(jì)算域以滿足第一類邊界條件.采用第三類邊界條件(即混合邊界條件)時(shí)，為降低邊界對(duì)模擬精度的影響，通常仍要求截?cái)噙吔缗c供電點(diǎn)距離較大，通過(guò)擴(kuò)大計(jì)算域來(lái)降低截?cái)噙吔缬绊?阮百堯等, 2001; 黃俊革等, 2003; 麻昌英等, 2018)，這使得建模難度和計(jì)算量增加，期望進(jìn)一步探索建模更方便且計(jì)算量更少的方法.Blome等(2009)采用無(wú)限元法處理邊界問(wèn)題，與Dirichlet和混合邊界條件相比縮小了計(jì)算范圍要求.湯井田和公勁喆(2010)在有限元-無(wú)限元法三維電阻率法正演中，提出了一種最優(yōu)的無(wú)限元形函數(shù)，并其與非結(jié)構(gòu)化四面體有限元相結(jié)合，取代了傳統(tǒng)混合邊界條件，使得電位在無(wú)限域內(nèi)連續(xù)并在無(wú)限遠(yuǎn)處衰減為零，該方法可以在近似測(cè)區(qū)大小的計(jì)算范圍內(nèi)得到與混合邊界條件相當(dāng)?shù)挠?jì)算精度，優(yōu)于相同計(jì)算范圍下齊次邊界條件的解，有利于減少計(jì)算節(jié)點(diǎn)數(shù).Ren和Tang (2014)在三維電阻率法有限單元正演中，多電極情況下將處于所有供電電極近似中間位置作為中心點(diǎn)，在計(jì)算域足夠大時(shí)將該中心點(diǎn)計(jì)算邊界條件近似代替所有電極的邊界條件，建立統(tǒng)一的多電極系統(tǒng)邊值問(wèn)題，使多電極情況下具有共同的系統(tǒng)矩陣.Yuan等(2016)在電阻率法有限單元2.5維正演中，采用一種新型的無(wú)限元形函數(shù)來(lái)代替無(wú)窮遠(yuǎn)邊界上的傳統(tǒng)混合邊界條件，有效減少了因人為施加邊界條件不適當(dāng)而引起的數(shù)值誤差，并且最終系統(tǒng)矩陣是稀疏的、對(duì)稱的，并且與源電極無(wú)關(guān).張錢(qián)江等(2016)在三維電阻率法有限單元正演中提出了一種近似邊界條件方法，該方法將與場(chǎng)源位置相關(guān)的邊界系數(shù)矩陣從剛度矩陣中分離出來(lái)，使得分離后的剛度矩陣與場(chǎng)源位置無(wú)關(guān)，并將邊界系數(shù)矩陣與邊界處一次電場(chǎng)向量的乘積移到線性方程組右端項(xiàng)中，當(dāng)場(chǎng)源位置發(fā)生改變時(shí)，只需要重新計(jì)算右端項(xiàng)就可以實(shí)現(xiàn)快速回代求解，數(shù)值精度優(yōu)于混合邊界條件.

基于戴前偉和王洪華(2013)、Feng等(2015)、Ma等(2019)提出的探地雷達(dá)和電阻率無(wú)單元Galerkin法正演中的移動(dòng)最小二乘(MLS)形函數(shù)，本文將MLS形函數(shù)應(yīng)用于電阻率法有限單元法正演的第三類邊界條件處理，以期望進(jìn)一步降低正演模擬對(duì)計(jì)算域范圍要求.在此基礎(chǔ)之上，提出電阻率法有限單元-移動(dòng)最小二乘(FEM-MLS)耦合正演算法.通過(guò)不同模擬方法對(duì)層狀模型測(cè)深和異常體模型進(jìn)行模擬對(duì)比，驗(yàn)證了本文算法的有效性，結(jié)果表明采用第三類邊界條件時(shí)，本文提出的FEM-MLS耦合法相比于有限單元法，在滿足精度要求情況下可進(jìn)一步縮小計(jì)算域，在小范圍計(jì)算域內(nèi)便可獲得較好的正演模擬精度，并提高了計(jì)算效率，具有更好的模擬性能.

1 基于三角形單元的電阻率法有限單元法正演

當(dāng)采用第三類邊界條件時(shí)，2.5維電阻率法波數(shù)域總電位U(x,λ)的邊值問(wèn)題及其對(duì)應(yīng)的變分問(wèn)題分別表示為(1)和(2)式(徐世浙，1994)：

(1)

(2)

(1)和(2)式中σ為介質(zhì)電導(dǎo)率，λ為波數(shù)，I0為電流，δ(A)為狄拉克函數(shù)，A為場(chǎng)源點(diǎn)，K0、K1分別為第二類零階、一階修正貝塞爾函數(shù)，Ω為問(wèn)題域(全局域)，Γs為地表邊界，Γ∞為截?cái)噙吔?，rA為場(chǎng)源點(diǎn)到截?cái)噙吔绲狞c(diǎn)的距離，n為邊界單位外法向，cos(rA,n)為n與矢徑rA構(gòu)成的夾角余弦.為了更好模擬復(fù)雜幾何地質(zhì)體和起伏地形，采用非結(jié)構(gòu)化三角形單元剖分計(jì)算域.有限單元法求解上述變分問(wèn)題相關(guān)理論已非常成熟，通過(guò)線性插值和對(duì)(1)式第一式求變分(Zienkiewicz and Taylor, 1989; 徐世浙, 1994)，在每個(gè)三角單元Ωe內(nèi)(1)式等價(jià)為(3)式.

KeUe=Fe,

(3)

式中Ke=Ke1+Ke2+Ke3為有限單元法單元子剛度系數(shù)矩陣，其中Ke1和Ke2對(duì)應(yīng)(2)式中的體積分，Ke3對(duì)應(yīng)截?cái)噙吔绶e分，F(xiàn)e為右端項(xiàng).本文對(duì)有限單元法中截?cái)噙吔邕M(jìn)行研究，Ke3表達(dá)式為

(4)

Ke3=(Ke3ij),i,j=1,2,3

(5)

式中的各個(gè)元素表達(dá)式寫(xiě)為

(6)

式中

(7)

(8)

(9)

KFEUFE=F,

(10)

式中KFE為由全部三角單元子域Ωe的Ke相加組成的有限單元法總系數(shù)矩陣，UFE為有限單元法所有三角單元節(jié)點(diǎn)上的波數(shù)域電位子向量Ue組合成的列向量，F(xiàn)為有限單元法方程右端項(xiàng)列向量.

2 移動(dòng)最小二乘近似

利用移動(dòng)最小二乘(Moving Least Square, MLS)近似構(gòu)造的形函數(shù)已在無(wú)單元法中被廣泛應(yīng)用.無(wú)單元法正演模擬中，需要為計(jì)算域Ω中任意點(diǎn)xT=(x,z)構(gòu)造局部支持域Ωq，假設(shè)Ωq中包含有n個(gè)節(jié)點(diǎn){x1,x2,…,xn}，它們對(duì)應(yīng)的波數(shù)域電位值為{U1,U2,…,Un}.假設(shè)計(jì)算域Ω中波數(shù)域電位函數(shù)U(x)，則計(jì)算域內(nèi)任意一點(diǎn)x處U(x)的MLS近似表達(dá)式為

(11)

式中，p(x)為坐標(biāo)xT=(x,z)的基函數(shù)向量，本文采用線性基函數(shù)pT(x)=(1xz).(11)式中a(x)為近似插值的系數(shù)向量，具體表示為

aT(x)=[a1(x)a2(x)a3(x)].

(12)

利用最小二乘法對(duì)(11)式進(jìn)行推導(dǎo)(Liu and Gu, 2005; 馮德山等, 2013; Ma et al., 2019)，最后得到

ΦT(x)=[φ1(x)φ2(x) …φn(x)]

=pT(x)A-1(x)B(x),

(13)

式中Φ(x)為計(jì)算點(diǎn)x的支持域中n個(gè)節(jié)點(diǎn)的MLS形函數(shù)向量.A(x)為基函數(shù)加權(quán)后形成的系數(shù)矩陣，定義為

(14)

圖1 移動(dòng)最小二乘近似示意圖Fig.1 Sketch showing the moving least squares approximation

B(x)=[w(r1)p(x1)w(r2)p(x2)…w(rn)p(xn)].

(15)

在移動(dòng)最小二乘形函數(shù)中，其系數(shù)a(x)為坐標(biāo)x的函數(shù)，這使得該加權(quán)最小二乘近似可連續(xù)移動(dòng)，故只要合理選擇連續(xù)的權(quán)函數(shù)，則移動(dòng)最小二乘形函數(shù)在全局計(jì)算域內(nèi)是連續(xù)的，其全局連續(xù)特性對(duì)于連續(xù)場(chǎng)的正演模擬是有利的.在無(wú)單元法正演模擬應(yīng)用中，指數(shù)函數(shù)和樣條函數(shù)常被作為權(quán)函數(shù)，本文中采用形式如(16)式的三次樣條函數(shù)作為權(quán)函數(shù)(Liu and Gu, 2005; Ma et al., 2019)，它具有二階連續(xù)性.

(16)

(17)

式中d=|x-xi|為計(jì)算點(diǎn)x與支持域內(nèi)節(jié)點(diǎn)xi的距離，rw為權(quán)函數(shù)域尺寸.

3 MLS形函數(shù)處理邊界

(18)

式中，Φ為MLS形函數(shù)向量，(18)式展開(kāi)得到

(19)

(20)

對(duì)于單個(gè)高斯點(diǎn)，子邊界系數(shù)矩陣Ke3的元素表達(dá)式為

(21)

式中，i,j=1,2,…,n表示支持域Ωq中的節(jié)點(diǎn)局部編號(hào).高斯點(diǎn)處電導(dǎo)率采用MLS形函數(shù)插值(22)式計(jì)算獲得.

(22)

4 模型算例

4.1 層狀模型

采用具有解析解的三層電阻率模型驗(yàn)證本文算法的有效性.建立水平方向?qū)?80 m(X：-190～190 m)，垂直方向深108 m(Z：0～108 m)的矩形計(jì)算域Ω1.如圖2所示，在地表設(shè)置52根電極作為供電和觀測(cè)電極，在地表X=0處進(jìn)行對(duì)稱四極裝置測(cè)深觀測(cè).計(jì)算域Ω1及其節(jié)點(diǎn)分布如圖3a所示，采用不規(guī)則分布的節(jié)點(diǎn)將計(jì)算域離散，節(jié)點(diǎn)在靠近電極和地層界面區(qū)域采用密集的節(jié)點(diǎn)分布，其他區(qū)域采用稀疏的節(jié)點(diǎn)分布，總節(jié)點(diǎn)數(shù)為14280.首先基于非結(jié)構(gòu)化三角形單元的有限單元法(FEM)對(duì)22組不同收發(fā)距(AB/2)對(duì)稱四極裝置測(cè)深觀測(cè)進(jìn)行正演模擬.本文利用具有高精度的漢克爾變換計(jì)算三層電阻率模型對(duì)稱四極裝置測(cè)深的解析解(蔡盛, 2014).

圖2 三層電阻率模型和對(duì)稱四極裝置測(cè)深觀測(cè)示意圖Fig.2 Sketch showing the three-layer resistivity model and symmetrical quadrupole device

圖3 三層電阻率模型計(jì)算域Ω1(a)和Ω2(b)及其節(jié)點(diǎn)分布示意圖Fig.3 Sketch showing the nodal distribution of three-layer resistivity model and computational domains Ω1 (a) and Ω2 (b)

如圖3b在計(jì)算域Ω1的基礎(chǔ)上，采用快速延伸的節(jié)點(diǎn)擴(kuò)展計(jì)算域，形成計(jì)算域Ω2，水平方向?qū)挸^(guò)2000 m，垂直方向高超過(guò)超過(guò)1000 m，共16800個(gè)節(jié)點(diǎn).對(duì)計(jì)算域Ω1分別采用有限單元法和FEM-MLS耦合法進(jìn)行電阻率法對(duì)稱四極測(cè)深模擬，對(duì)計(jì)算域Ω2采用有限單元法模擬.采用FEM-MLS耦合法時(shí)，邊界積分計(jì)算中在子邊界上使用兩個(gè)高斯點(diǎn)計(jì)算，使用三角形單元的三個(gè)節(jié)點(diǎn)作為高斯點(diǎn)支持域.

圖4和圖5分別為中間層為高阻和低阻時(shí)，不同方法和不同計(jì)算域的正演視電阻率和相對(duì)誤差曲線.由圖4和圖5的正演視電阻率和相對(duì)誤差曲線分析可知，無(wú)論中間層為高阻或低阻，有限單元法采用計(jì)算域Ω1時(shí)，正演視電阻率曲線尾部數(shù)值偏低，相對(duì)誤差較大，F(xiàn)EM-MLS耦合法采用計(jì)算域Ω1和有限單元法采用計(jì)算域Ω2時(shí)，正演視電阻率曲線尾部模擬精度明顯得到改善，模擬結(jié)果與解析解吻合很好.表1為不同方法和不同計(jì)算域正演視電阻率平均相對(duì)誤差，相比于有限單元法采用計(jì)算域Ω1，F(xiàn)EM-MLS耦合法采用計(jì)算域Ω1和有限單元法采用計(jì)算域Ω2時(shí)均可有效提高模擬精度，它們的平均相對(duì)誤差均小于1%，且相差很小.上述分析表

圖4 中間層為高阻時(shí)三層電阻率模型模擬視電阻率(a)和相對(duì)誤差(b)曲線Fig.4 The apparent resistivity and relative errors of three-layer resistivity model with the resistance middle layer(a) Apparent resistivity; (b) Relative errors.

圖5 中間層為低阻時(shí)三層電阻率模型模擬視電阻率(a)和相對(duì)誤差(b)曲線Fig.5 The apparent resistivity and relative errors of three-layer resistivity model with the conductive middle layer(a) Apparent resistivity; (b) Relative errors.

表1 不同方法和不同計(jì)算域三層模型測(cè)深視電阻率平均相對(duì)誤差(%)Table 1 The relative errors (%) of three-layer resistivity model by using different methods and computational domains

明，通過(guò)擴(kuò)邊減小邊界的影響可有效提高有限單元法模擬精度，F(xiàn)EM-MLS耦合法可在較小范圍計(jì)算域內(nèi)獲得良好的模擬精度，與有限單元法采用擴(kuò)邊后的計(jì)算域正演模擬效果相當(dāng)，相比于采用相同較小范圍計(jì)算域的有限單元法平均精度(平均相對(duì)誤差)提高了約一倍，驗(yàn)證了FEM-MLS耦合法的有效性.因此FEM-MLS耦合法可使用較小的計(jì)算域便可獲得良好的模擬精度，相比于有限單元法不需要擴(kuò)邊處理，可使得正演建模更加方便，使用更小的計(jì)算域也使得計(jì)算量更小，提高正演模擬效率.

為了研究FEM-MLS耦合法中MLS形函數(shù)計(jì)算邊界積分時(shí)參數(shù)選擇與模擬精度的關(guān)系，首先在子邊界上分別采用高斯點(diǎn)數(shù)G=2，4，6對(duì)三層模型對(duì)稱四極測(cè)深模擬，使用三角形單元的三個(gè)節(jié)點(diǎn)作為高斯點(diǎn)支持域，采用計(jì)算域Ω1計(jì)算.圖6和圖7分別為中間層為高阻和低阻時(shí)，采用不同高斯點(diǎn)數(shù)的正演視電阻率和相對(duì)誤差曲線.由圖6和圖7分析可知，無(wú)論中間層為高阻或低阻，采用不同高斯點(diǎn)數(shù)獲得正演視電阻率和相對(duì)誤差基本是重合的，表明采用不同高斯點(diǎn)數(shù)模擬精度基本一致，同等條件下高斯點(diǎn)數(shù)對(duì)邊界積分影響不大.因此，為減少計(jì)算量在子邊界采用兩個(gè)高斯點(diǎn)計(jì)算邊界積分.

圖6 中間層為高阻時(shí)不同高斯點(diǎn)數(shù)三層電阻率模型模擬視電阻率(a)和相對(duì)誤差(b)曲線(a) 視電阻率; (b) 相對(duì)誤差.Fig.6 The apparent resistivity and relative errors of three-layer resistivity model with the resistance middle layer by using different number of Gauss points(a) Apparent resistivity; (b) Relative errors.

圖7 中間層為低阻時(shí)不同高斯點(diǎn)數(shù)三層電阻率模型模擬視電阻率(a)和相對(duì)誤差(b)曲線.Fig.7 The apparent resistivity and relative errors of three-layer resistivity model with the conductive middle layer by using different number of Gauss points(a) Apparent resistivity; (b) Relative errors.

分別采用支持域節(jié)點(diǎn)數(shù)D=2，3，4，其中當(dāng)D=2時(shí)為采用子邊界端點(diǎn)對(duì)應(yīng)的2個(gè)節(jié)點(diǎn)，即支持域兩個(gè)節(jié)點(diǎn)全為子邊界端點(diǎn)，當(dāng)D=3時(shí)為使用三角形單元的三個(gè)節(jié)點(diǎn)，當(dāng)D=4時(shí)為處三角形單元的三個(gè)節(jié)點(diǎn)外，加上距離高斯點(diǎn)最近的一個(gè)節(jié)點(diǎn).圖8和圖9分別為中間層為高阻和低阻時(shí)，支持域采用不同節(jié)點(diǎn)數(shù)的正演視電阻率和相對(duì)誤差曲線.由圖8和圖9分析可知，無(wú)論中間層為高阻或低阻，支持域采用不同節(jié)點(diǎn)數(shù)獲得正演視電阻率和相對(duì)誤差基本是重合的，支持域采用不同節(jié)點(diǎn)數(shù)模擬精度基本一致，表明同等條件下FEM-MLS耦合法邊界積分中支持域采用不同節(jié)點(diǎn)數(shù)對(duì)邊界積分影響不大.支持域內(nèi)在包含單元三個(gè)節(jié)點(diǎn)基礎(chǔ)上若再增加節(jié)點(diǎn)數(shù)時(shí)(D>3)，模擬結(jié)果與解析解也基本吻合(基本與上述模擬結(jié)果重合，因此未圖示模擬結(jié)果)，支持域節(jié)點(diǎn)數(shù)增加意味著計(jì)算量增加.

圖8 中間層為高阻時(shí)支持域不同節(jié)點(diǎn)數(shù)三層電阻率模型模擬視電阻率(a)和相對(duì)誤差(b)曲線Fig.8 The apparent resistivity and relative errors of three-layer resistivity model with the resistance middle layer by using different number of nodes in the support domain(a) Apparent resistivity; (b) Relative errors.

圖9 中間層為低阻時(shí)支持域不同節(jié)點(diǎn)數(shù)三層電阻率模型模擬視電阻率(a)和相對(duì)誤差(b)曲線Fig.9 The apparent resistivity and relative errors of three-layer resistivity model with the conductive middle layer by using different number of nodes in the support domain(a) Apparent resistivity; (b) Relative errors.

為了研究FEM-MLS耦合法中MLS形函數(shù)計(jì)算邊界積分時(shí)支持域節(jié)點(diǎn)位置與模擬精度的關(guān)系，當(dāng)D=2時(shí)，采用子邊界端點(diǎn)其中對(duì)應(yīng)的一個(gè)節(jié)點(diǎn)和三角形單元非子邊界上的一個(gè)節(jié)點(diǎn)作為支持域，即支持域兩個(gè)節(jié)點(diǎn)非全為子邊界端點(diǎn).圖10和圖11分別為中間層為高阻和低阻時(shí)正演視電阻率和相對(duì)誤差曲線.由圖10和圖11分析可知，無(wú)論中間層為高阻或低阻，當(dāng)支持域兩個(gè)節(jié)點(diǎn)全為子邊界端點(diǎn)時(shí)模擬結(jié)果與解析解吻合很好，相對(duì)誤差較小，而當(dāng)支持域兩個(gè)節(jié)點(diǎn)全非為子邊界端點(diǎn)時(shí)正演視電阻率曲線尾支與解析解出現(xiàn)一定的偏離，相對(duì)誤差曲線尾支逐漸增大，模擬精度下降.

通過(guò)上述模擬結(jié)果分析表明，當(dāng)采用MLS形函數(shù)計(jì)算邊界積分時(shí)，在同等情況下，F(xiàn)EM-MLS耦合法中子邊界上使用不同的高斯點(diǎn)數(shù)對(duì)模擬精度影響不大，采用不同高斯點(diǎn)數(shù)時(shí)模擬結(jié)果基本一致；高斯點(diǎn)支持域內(nèi)采用不同節(jié)點(diǎn)數(shù)及其組合形式時(shí)對(duì)模擬精度影響較大，其中當(dāng)子邊界端點(diǎn)對(duì)應(yīng)的兩個(gè)節(jié)點(diǎn)不全被包含于支持域內(nèi)時(shí)，可明顯降低收發(fā)距相對(duì)較大時(shí)的模擬精度.因此，在采用FEM-MLS耦合法模擬時(shí)，需確保子邊界兩個(gè)端點(diǎn)對(duì)應(yīng)的節(jié)點(diǎn)均包含于支持域內(nèi).此外，由于僅采用子邊界端點(diǎn)對(duì)應(yīng)的兩個(gè)節(jié)點(diǎn)作為支持域時(shí)，兩個(gè)節(jié)點(diǎn)在同一直線上，有可能使得MLS形函數(shù)構(gòu)造不穩(wěn)定，而若采用較多支持域節(jié)點(diǎn)時(shí)則需要更多的節(jié)點(diǎn)搜索.直接采用三角形單元的三個(gè)節(jié)點(diǎn)作為支持域無(wú)需額外構(gòu)造支持域，同時(shí)只要三角形單元形狀非特別尖銳情況下，通?？杀ＷCMLS形函數(shù)構(gòu)造穩(wěn)定性.基于本小節(jié)對(duì)三層地電模型模擬結(jié)果分析，F(xiàn)EM-MLS耦合法可采用三角形單元的三個(gè)節(jié)點(diǎn)作為支持域構(gòu)造MLS形函數(shù)，子邊界上采用兩個(gè)高斯點(diǎn)計(jì)算邊界積分，在確保模擬精度前提下以減少計(jì)算量和保證計(jì)算穩(wěn)定性.

圖10 中間層為高阻時(shí)支持域采用2個(gè)節(jié)點(diǎn)三層電阻率模型模擬視電阻率(a)和相對(duì)誤差(b)曲線Fig.10 The apparent resistivity and relative errors of three-layer resistivity model with the resistance middle layer by using two nodes in the support domain(a) Apparent resistivity; (b) Relative errors.

圖11 中間層為低阻時(shí)支持域采用2個(gè)節(jié)點(diǎn)三層電阻率模型模擬視電阻率(a)和相對(duì)誤差(b)曲線Fig.11 The apparent resistivity and relative errors of three-layer resistivity model with the conductive middle layer by using two nodes in the support domain(a) Apparent resistivity; (b) Relative errors.

4.2 水平圓柱狀異常體

采用水平圓柱狀異常體模型進(jìn)一步測(cè)試FEM-MLS耦合法的模擬性能.建立水平方向?qū)?00 m(X：-100～100 m)，垂直方向高100 m(Z：0～100 m)的計(jì)算域Ω1，背景介質(zhì)電阻率ρ1=100 Ωm，在模型中添加一個(gè)水平圓柱狀異常體，其電阻率為ρ2=1000 Ωm或ρ2=10 Ωm，水平圓柱狀異常體圓心位于(x,z)=(0 m, 15 m)，半徑為5 m，水平圓柱狀異常體模型和節(jié)點(diǎn)分布如圖12a所示，共12038個(gè)節(jié)點(diǎn).如圖12b所示，在計(jì)算域Ω1的基礎(chǔ)上，采用快速延伸的節(jié)點(diǎn)擴(kuò)展計(jì)算域，形成計(jì)算域Ω2，水平方向?qū)挸^(guò)2000 m，垂直方向高超過(guò)1000 m，共13872個(gè)節(jié)點(diǎn).對(duì)計(jì)算域Ω1分別采用有限單元法(FEM)和FEM-MLS耦合法進(jìn)行電阻率法對(duì)稱四極測(cè)深模擬，對(duì)計(jì)算域Ω2采用有限單元法模擬.

圖12 水平圓柱狀異常體模型不同計(jì)算域Ω1～Ω5及其節(jié)點(diǎn)分布示意圖(a)—(e) 分別為計(jì)算域Ω1～Ω5.Fig.12 Sketch showing different computational domains Ω3～Ω5 and the nodal distribution of the horizontal cylindrical anomaly model(a)—(e) is computational domains Ω3～Ω5, respectively.

圖13為有限單元法采用計(jì)算域Ω1正演的模擬結(jié)果，圖14為有限單元法采用計(jì)算域Ω2正演的模擬結(jié)果，圖15為本文FEM-MLS耦合法采用計(jì)算域Ω1正演的模擬結(jié)果.基于4.1節(jié)層狀模型分析結(jié)論，將有限單元法采用計(jì)算域Ω2正演的模擬結(jié)果圖14作為參考模擬結(jié)果.對(duì)比分析圖13和圖 14可知，有限單元法兩種計(jì)算域的模擬結(jié)果有明顯的差異，采用計(jì)算域Ω1正演的模擬結(jié)果與層狀模型類似，在收發(fā)距(AB/2)相對(duì)較大時(shí)正演視電阻率偏低，模擬精度較差.對(duì)比分析圖14和圖15可知，F(xiàn)EM-MLS耦合法采用計(jì)算域Ω1正演的模擬結(jié)果與參考結(jié)果基本相同.以參考結(jié)果(圖14)為基礎(chǔ)，計(jì)算FEM-MLS耦合法采用計(jì)算域Ω1正演的模擬結(jié)果相對(duì)誤差，相對(duì)誤差分布如圖16所示.由圖16分析，無(wú)論異常體為高阻還是低阻，相對(duì)誤差均很小，異常體為高阻時(shí)最大相對(duì)誤差不超過(guò)0.06%，異常體為低阻時(shí)最大相對(duì)誤差不超過(guò)0.05%，表明FEM-MLS耦合法在較小范圍計(jì)算域條件下可獲得良好的模擬精度，進(jìn)一步驗(yàn)證了FEM-MLS耦合法的有效性.

圖13 不同電阻率異常體有限單元法采用計(jì)算域Ω1獲得的視電阻率擬斷面圖(a) 高阻水平圓柱狀異常體; (b) 低阻水平圓柱狀異常體.Fig.13 The apparent resistivity of forward modeling results by using the FEM method with computational domain Ω1(a) Resistance anomaly; (b) Conductive anomaly.

圖14 不同電阻率異常體有限單元法采用計(jì)算域Ω2獲得的視電阻率擬斷面圖(a) 高阻水平圓柱狀異常體; (b) 低阻水平圓柱狀異常體.Fig.14 The apparent resistivity of forward modeling results by using the FEM method with computational domain Ω2(a) Resistance anomaly; (b) Conductive anomaly.

圖16 不同電阻率異常體FEM-MLS耦合法采用計(jì)算域Ω1模擬結(jié)果相對(duì)誤差分布(a) 高阻水平圓柱狀異常體; (b) 低阻水平圓柱狀異常體.Fig.16 The relative errors of forward modeling results by using the FEM-MLS coupling method with computational domain Ω1(a) Resistance anomaly; (b) Conductive anomaly.

為了進(jìn)一步研究FEM-MLS耦合法的模擬性能，采用范圍更小的計(jì)算域?qū)λ綀A柱狀異常體模型進(jìn)行模擬.在計(jì)算域Ω1及其節(jié)點(diǎn)分布基礎(chǔ)上，將計(jì)算域范圍縮小為不同范圍的計(jì)算域Ω3～Ω5(計(jì)算域按序號(hào)逐漸減小)，計(jì)算域Ω3～Ω5分別如圖12(c—e)所示.計(jì)算域Ω3：水平方向?qū)?60 m(X：-80～80 m)，垂直方向高60 m(Z：0～60 m)的，共11149個(gè)節(jié)點(diǎn)；計(jì)算域Ω4：水平方向?qū)?40 m(X：-70～70 m)，垂直方向高50 m(Z：0～50 m)的，共8519個(gè)節(jié)點(diǎn)；計(jì)算域Ω5：水平方向?qū)?30 m(X：-65～65 m)，垂直方向高40 m(Z：0～40 m)的，共6699個(gè)節(jié)點(diǎn).利用FEM-MLS耦合法分別采用計(jì)算域Ω3～Ω5進(jìn)行正演模擬，模擬結(jié)果如圖17所示.對(duì)比分析圖17和圖14可知，F(xiàn)EM-MLS耦合法采用計(jì)算域Ω3～Ω5的模擬結(jié)果與參考結(jié)果基本相同.以參考結(jié)果(圖14)為基礎(chǔ)，計(jì)算FEM-MLS耦合法采用計(jì)算域Ω3～Ω5的模擬結(jié)果相對(duì)誤差，相對(duì)誤差分布如圖18所示.

圖17 FEM-MLS耦合法采用計(jì)算域Ω3～Ω5獲得的視電阻率擬斷面圖(a) 高阻異常體; (b) 低阻異常體.Fig.17 The apparent resistivity of forward modeling results by using the FEM-MLS coupling method with different computational domains Ω3～Ω5(a) Resistance anomaly; (b) Conductive anomaly.

圖18 FEM-MLS耦合法采用計(jì)算域Ω3～Ω5的模擬結(jié)果相對(duì)誤差分布(a) 高阻異常體; (b) 低阻異常體.Fig.18 The relative errors of forward modeling results by using the FEM-MLS coupling method with different computational domains Ω3～Ω5(a) Resistance anomaly; (b) Conductive anomaly.

對(duì)比圖18和圖16可知，采用計(jì)算域Ω3～Ω5模擬時(shí)，水平圓柱狀異常體為高阻和低阻的模擬結(jié)果相對(duì)誤差均有所增加，并隨著計(jì)算域范圍的減小相對(duì)誤差逐漸增大，其中異常體為高阻時(shí)增加較低阻異常體明顯.表2列出了FEM-MLS耦合法采用計(jì)算域Ω3～Ω5時(shí)模擬結(jié)果的最大相對(duì)誤差和平均相對(duì)誤差，異常體為高阻或低阻時(shí)平均相對(duì)誤差均很小，異常體為高阻時(shí)最大相對(duì)誤差不超過(guò)0.4%，異常體為低阻時(shí)最大相對(duì)誤差不超過(guò)0.2%，表明當(dāng)計(jì)算域范圍進(jìn)一步減小時(shí)，F(xiàn)EM-MLS耦合法仍然能獲得良好的模擬精度.以上分析進(jìn)一步表明，在保證模擬精度前提下，有限單元法中采用無(wú)單元法計(jì)算邊界積分，可有效地縮小計(jì)算域范圍要求，相比于有限單元法，采用更小的計(jì)算域便可得到良好的模擬結(jié)果.

表2 FEM-MLS耦合法采用計(jì)算域Ω1，Ω3～Ω5模擬結(jié)果的相對(duì)誤差Table 2 The relative errors (%) of simulation results by using the FEM-MLS coupling method with different computational domains: Ω1，Ω3～Ω5

縮小計(jì)算域范圍要求不僅使得建模更加方便，同時(shí)也可減少節(jié)點(diǎn)和單元數(shù)量，提高計(jì)算效率.表3為不同方法和不同計(jì)算域單元數(shù)、節(jié)點(diǎn)數(shù)和計(jì)算時(shí)間(為對(duì)52根電極進(jìn)行了電位正演模擬耗時(shí)).有限單元法采用計(jì)算域Ω2時(shí)，由于計(jì)算域Ω2在計(jì)算域Ω1基礎(chǔ)上進(jìn)行了擴(kuò)邊，使用了更多的節(jié)點(diǎn)和單元，因此增加了計(jì)算時(shí)間.FEM-MLS耦合法采用計(jì)算域Ω1計(jì)算不需要擴(kuò)邊處理，計(jì)算效率與采用計(jì)算域Ω1的有限單元法一致.FEM-MLS耦合法采用計(jì)算域Ω3～Ω5計(jì)算時(shí)，由于計(jì)算域范圍更小，相比于計(jì)算域Ω1和Ω2，節(jié)點(diǎn)和單元數(shù)更少，因此計(jì)算量更少，有效提高了計(jì)算效率.結(jié)果表明，在計(jì)算精度(平均相對(duì)誤差)基本一致的情況下，但本文算法采用較小范圍的計(jì)算域，建模方便，使用更少的單元和節(jié)點(diǎn)，相比于采用較大計(jì)算域滿足邊界條件的有限單元法計(jì)算效率提高了約一倍.

表3 不同方法和不同計(jì)算域節(jié)點(diǎn)數(shù)、單元數(shù)和計(jì)算時(shí)間Table 3 Number of nodes, number of cells and calculational cost of different methods using different computational domains

5 結(jié)論

本文提出一種具有高精度、高效的電阻率法有限單元-移動(dòng)最小二乘(FEM-MLS)耦合正演算法，該算法將無(wú)單元法中的移動(dòng)最小二乘(MLS)形函數(shù)應(yīng)用于電阻率法有限單元正演的第三類邊界條件處理.本文算例分析表明，采用第三類邊界條件時(shí)，F(xiàn)EM-MLS耦合法可在較小范圍計(jì)算域內(nèi)獲得良好的模擬精度，與有限單元法采用較大范圍計(jì)算域的正演模擬效果相當(dāng)，驗(yàn)證了本文算法的有效性.本文算法采用較小范圍的計(jì)算域，建模方便，使用更少的單元和節(jié)點(diǎn)，相比于采用較大計(jì)算域滿足邊界條件的有限單元法計(jì)算效率提高了約一倍，相比于采用相同較小范圍計(jì)算域的有限單元法平均精度提高了約一倍.通過(guò)數(shù)值模擬分析了FEM-MLS耦合法的參數(shù)選擇對(duì)正演模擬的影響.分析表明，F(xiàn)EM-MLS耦合法中子邊界上使用不同的高斯點(diǎn)數(shù)對(duì)模擬精度影響不大，不同高斯點(diǎn)支持域內(nèi)節(jié)點(diǎn)數(shù)及其組合形式對(duì)模擬精度影響較大，需確保子邊界兩個(gè)端點(diǎn)對(duì)應(yīng)的節(jié)點(diǎn)均包含于支持域內(nèi).基于對(duì)三層地電模型模擬結(jié)果分析，F(xiàn)EM-MLS耦合法在處理邊界條件時(shí)，建議子邊界上采用三角形單元的三個(gè)節(jié)點(diǎn)作為支持域構(gòu)造MLS形函數(shù)，采用兩個(gè)高斯點(diǎn)計(jì)算邊界積分.

致謝感謝審稿人對(duì)完善本文提出的寶貴的具有建設(shè)性的意見(jiàn)和建議.在本文研究過(guò)程中，中南大學(xué)的劉海飛副教授、郭榮文研究員、崔益安教授、孫婭副教授及劉嶸特聘副教授給予了很多幫助和寶貴意見(jiàn)，在此表示感謝.