基于有限元超收斂的三維節(jié)點型有限元后處理技術(shù)

湯井田 廖濤山 陳 煌* 皇祥宇 周 峰 張林成

(①中南大學有色金屬成礦預(yù)測與地質(zhì)環(huán)境監(jiān)測教育部重點實驗室，湖南長沙 410083； ②中南大學地球科學與信息物理學院，湖南長沙 410083； ③中國能源建設(shè)集團湖南省電力設(shè)計院有限公司，湖南長沙 410083； ④東華理工大學地球物理與測控技術(shù)學院，江西南昌 330013； ⑤湖南城市學院信息與電子工程學院，湖南益陽 413002)

0 引言

在地球物理電磁法中，正演是反演的基礎(chǔ)。正演方法主要包括邊界單元法、有限差分法、積分方程法、有限元法等，其中有限元法因其對任意復雜地形及復雜地質(zhì)構(gòu)造的適應(yīng)能力強、計算精度較高等優(yōu)點被廣泛應(yīng)用[1]。無論是關(guān)于位的還是關(guān)于場的方程，都需要對有限元的數(shù)值解進行數(shù)值微分后處理才能進行阻抗計算，進而獲得視電阻率和相位響應(yīng)。直接從有限元解計算的梯度值在單元邊界上不連續(xù)且整體精度不高，雖然加密剖分網(wǎng)格或者增加單元形函數(shù)插值次數(shù)一定程度上可以提高精度，然而隨著網(wǎng)格的加密和插值次數(shù)的提高，有限元方法產(chǎn)生的線性方程組的未知數(shù)將按幾何級數(shù)急劇增加，這極大地增加了對計算機內(nèi)存的需求，計算時間急劇增加。因此，對有限元數(shù)值解進行恰當?shù)暮筇幚硪蕴岣哂邢拊馓荻戎档木?，是有限元?shù)值模擬中的一項重要內(nèi)容。

對于有限單元法的后處理，最早是直接對單元形函數(shù)微分(Shape-function differentiation，SFD)求取輔助場，這種方法在單元邊界上誤差較大[2]；Rodi[3]針對大地電磁法(Magnetotellurics，MT)矩形單元雙線性插值的二維正演使用矩量法(the Method of Moment，MOM)，其實質(zhì)是將輔助場的計算函數(shù)定義為與主場形函數(shù)相同的形式，用主場值重新定義函數(shù)中的列矢量，目的是讓輔助場自動滿足MT的邊界條件，但實現(xiàn)過程非常復雜且繁瑣；陳樂壽等[4]進行MT二維正演時，為了適應(yīng)傾斜界面，整體采用矩形網(wǎng)格剖分，而在分界面上將矩形單元按對角線分為兩個直角三角形，求主場時直接在三角形單元上采用二次插值，然后基于形函數(shù)微分提高輔助場的精度，取得了不錯的效果，但這種方法依然存在計算量和對內(nèi)存需求過大的問題；陳小斌等[5-7]提出利用有限元直接迭代法實現(xiàn)MT二維正演模擬及線源頻率電磁響應(yīng)的計算，由有限元直接迭代格式得到地表三角形單元上各邊中點的主場值，從而在單元上構(gòu)造二次插值的形函數(shù)，然后基于形函數(shù)微分計算輔助場，該方法得到的輔助場精度有所提高，但并不適于以有限元數(shù)值模擬為前提的后處理。

對于結(jié)構(gòu)化網(wǎng)格的節(jié)點有限元正演，拉格朗日插值(LI)法是常用的后處理方法?；谠摲椒?，徐世浙[8]實現(xiàn)了MT的正演模擬，張繼峰[9]實現(xiàn)了可控源電磁法(Controlled Source Electromagnetic Method，CSEM)的三維有限元正演，張林成等[10]實現(xiàn)了CSEM的三維有限元—無限元耦合正演。LI法盡管理論簡單且計算量非常小，但其精度嚴重依賴網(wǎng)格大小，且要求測點附近的網(wǎng)格被等距剖分。在非結(jié)構(gòu)化網(wǎng)格有限元的后處理中，移動最小二乘插值法(Moving least-square interpolation，MLSI)應(yīng)用最為廣泛。Badea等[11]在開展基于Coulomb位的有限元井中電磁正演模擬時，采用該算法進行后處理；Vladimir等[12]、蔡紅柱等[13]及陳漢波等[14]在開展基于Coulomb位的完全非結(jié)構(gòu)化網(wǎng)格有限元正演時，都沿用了MLSI算法，僅在權(quán)函數(shù)的使用上有所差別。盡管MLSI法可以恢復任意點的梯度值，但精度依然不高。

為獲得高精度且穩(wěn)定的后處理數(shù)值結(jié)果，本文引入一種基于超收斂單元片梯度值恢復(Super-convergence Patch Recovery，SPR)的后處理技術(shù)，并以可控源電磁法節(jié)點有限元正演作為應(yīng)用基礎(chǔ)，對該算法進行了詳細的說明和驗證。

Zienkiewicz等[15-20]發(fā)現(xiàn)有限元解的導數(shù)在某些特殊點上具有特別高的精度，并將這種有限元的超收斂性質(zhì)用于漸進準確的Z-Z后驗誤差估計方法，系統(tǒng)地提出了SPR技術(shù)。由于SPR技術(shù)具有計算簡潔、容易理解、效果顯著及現(xiàn)有的有限元接口方便等優(yōu)勢，在計算流體力學、結(jié)構(gòu)力學等領(lǐng)域已經(jīng)得到廣泛應(yīng)用。在地球物理領(lǐng)域，Ren等[21]基于該方法實現(xiàn)了直流電阻率法的三維非結(jié)構(gòu)化網(wǎng)格自適應(yīng)有限元正演。

本文以可控源電磁法為應(yīng)用基礎(chǔ)，基于電場二次場的雙旋度方程，采用結(jié)構(gòu)化六面體單元的伽遼金有限元法對電場雙旋度方程進行離散化；然后以單元片上高斯點的電場分量梯度值進行最小二乘曲面擬合，恢復單元片上地表測線節(jié)點上的電場梯度；最后根據(jù)電場分量的旋度方程計算地表測線上的磁場分量。地電模型分析表明，在保證主場有限元數(shù)值模擬程序可靠的前提下，與SFD、LI和MLSI方法相比，SPR后處理技術(shù)能夠以極小的計算量顯著提高磁場分量的計算精度，且相對誤差較穩(wěn)定。

1 SPR方法原理

1.1 CSEM滿足的邊值問題

對CSEM的電場二次場的雙旋度方程引入散度校正項[9]

=iωμ0(σ-σp)Ep

(1)

式中：ω為角頻率；μ0為真空中磁導率；ε0為真空中介電常數(shù)；σ為電導率；σp為背景電導率；Ep和Es分別是電場一次場和二次場。

可控源電磁法數(shù)值模擬中，為保證控制微分方程(式(1))解的唯一性，常加入無窮遠邊界條件，令外邊界Γ上的二次場Es衰減為零

Es=0 ∈Γ

(2)

本文采用伽遼金節(jié)點有限元法求解邊值問題(式(1)和式(2))，剖分網(wǎng)格采用結(jié)構(gòu)化六面體單元，單元形函數(shù)采用三線性插值。整個研究區(qū)域劃分為目標區(qū)域和網(wǎng)格邊界區(qū)域：測線和源所在的區(qū)域是目標區(qū)域，采用均勻網(wǎng)格剖分；向外延伸區(qū)域是網(wǎng)格邊界區(qū)域，采用逐漸向外擴大的網(wǎng)格剖分,以降低邊界的影響，提高數(shù)值模擬的精度。關(guān)于有限元求解的細節(jié)參考文獻[9]。

1.2 有限元后處理算法

本文節(jié)點有限元數(shù)值模擬的解僅針對電場。為了獲得卡尼亞視電阻率、阻抗和相位等電磁響應(yīng)，還需要利用有效的數(shù)值微分后處理算法，由電場的旋度方程恢復出高精度的磁場，其計算公式為

(3)

下面首先闡述本文引入SPR技術(shù)的原理。為了與常規(guī)的有限元后處理算法進行對比，分別介紹、對比SFD、LI和MLSI三種有限元后處理算法。

1.2.1 SPR后處理技術(shù)

根據(jù)Herrmann理論，有限元解及其梯度在某些特殊的點上具有至少高1階的精度，這種現(xiàn)象稱為超收斂，這些特殊點稱為超收斂點[15]。有限元解的超收斂點在節(jié)點上，而有限元解梯度的超收斂點一般與高斯點對應(yīng)。

SPR技術(shù)的原理簡述如下：由圍繞任一公共節(jié)點的所有相關(guān)單元組成一單元片，一般將單元片上高斯點的有限元解梯度作為采樣對象，設(shè)近似多項式進行最小二乘曲面擬合，以恢復單元片上節(jié)點或其他特定點的梯度值，恢復的梯度值同樣具有超收斂性。

單元片一般是由以待恢復點作為公共節(jié)點的所有單元組合而成，如圖1所示。若待恢復點位于介質(zhì)分界面或者介質(zhì)邊界時，選擇界面任意一側(cè)的鄰近節(jié)點作為中心組成單元片，如圖2所示。

圖1 二維線性插值單元片

圖2 介質(zhì)邊界及分界面上的二維單元片

不是所有單元高斯點上梯度的收斂階都與最佳收斂階相對應(yīng)，如二次插值的三角形單元，由于高斯點上的梯度值仍然具有很高的精度，對于梯度恢復仍然具有重要意義。單元形狀和有限元插值階次的不同都會造成高斯點上梯度的收斂階的差異，因而恢復結(jié)果的精度也存在差異。在近似多項式的設(shè)置上，一般與有限元單元形函數(shù)的形式保持一致，可以達到最佳恢復效果。在這一過程中，若待恢復點位于介質(zhì)邊界或者不同介質(zhì)分界面時，則該點可同時隸屬于兩側(cè)不同介質(zhì)的單元片，因而根據(jù)其所屬的單元片恢復結(jié)果取均值。

本文有限元正演采用六面體單元及三線性插值的形函數(shù)。由于測線一般都位于地空分界面(即地面)上，以地下一側(cè)的單元合成單元片，則測點位于單元片的邊界上，如圖3所示。具體后處理過程如下。

圖3 六面體單元組成的單元片(三線性插值)示意圖

首先，由單元形函數(shù)計算單元片上各采樣點上的電場的梯度值?Ei，其中i表示采樣點編號。如圖3所示六面體單元，以位于單元幾何中心的高斯點作為采樣點，每個單元片上共計n=8個采樣點。

在單元片上一般以單元片集中點作為原點建立坐標系(圖3)。在單元片上設(shè)擬合多項式的形式與三線性插值保持一致

?E*=Pm

(4)

式中：?E*表示電場梯度擬合式；多項式P=[1，x,y,z,xy,xz,yz,xyz]； 系數(shù)向量m=[m1,m2,m3,m4,m5,m6,m7,m8]T。

在每個單元片上，對于n個采樣點的電場梯度值，極小化最小二乘泛函

(5)

然后，分別對系數(shù)向量m的各分量求偏導，整理后得到

m=A-1k

(6)

由上述方程組求解得到多項式的系數(shù)向量m后，便可根據(jù)式(4)恢復單元片上的地表測線上節(jié)點的電場分量梯度，并由式(3)計算磁場分量。

1.2.2 單元形函數(shù)微分法(SFD)

在剖分區(qū)域內(nèi)采用矩形六面體網(wǎng)格進行離散，母單元與子單元坐標對應(yīng)關(guān)系為

(7)

式中：(ξ，η，γ)是母單元坐標； (x,y,z)是子單元坐標； (x0,y0,z0)是子單元的中心坐標；a、b、c分別是子單元x、y、z三個方向的棱長。

三線性插值的單元形函數(shù)統(tǒng)一表示為

(8)

式中(ξi，ηi，γi)表示母單元上第i個節(jié)點的坐標。

六面體單元上任意點的電場分量可表示為

(9)

根據(jù)上式可得空間各方向上的電場分量梯度為

(10)

最后，根據(jù)式(3)可計算磁場分量。

1.2.3 拉格朗日插值(LI)法

LI法是已知節(jié)點的電場分量，令插值基函數(shù)lj(x)在點xj處取值為1，在其他點xi(i≠j)處取值為0。

所選測線處于均勻剖分的網(wǎng)格區(qū)域，在x和y方向采用三點中心差分、在z方向地表以下采用四點向前差分求電場梯度，可避免高階插值出現(xiàn)振蕩現(xiàn)象。計算公式為

(11)

1.2.4 移動最小二乘(MLSI)法

MLSI法被廣泛應(yīng)用于非結(jié)構(gòu)化網(wǎng)格有限元正演的后處理。與SPR后處理技術(shù)不同，它的采樣對象是節(jié)點上的電場分量，同時引入權(quán)函數(shù)，使距離中心點越遠的點被賦予越小的權(quán)重，恢復局部域中任意位置的電場梯度[22]。

對于任意點K恢復其梯度值，以該點為原點、以半徑R確定一個球體范圍局部域，以域中節(jié)點上的電場分量f*作為采樣對象。

設(shè)局部域上電場分量的線性插值函數(shù)為

f=Uα

(12)

式中：多項式U=[1,x,y,z]；系數(shù)向量α=[α1,α2,α3,α4]T。

設(shè)加權(quán)函數(shù)隨與點K的距離增大而單調(diào)降低

ω(j)=e-β2[(xj/xm)2+(yj/ym)2+(zj/zm)2]

(13)

式中：j表示采樣點編號；xm=Max(|xj|)；ym=Max(|yj|)； 系數(shù)β控制權(quán)函數(shù)逐漸遞減到零的速度，本文取β=2。

對于局部域上的N個采樣點求極小化問題

(14)

與前文泛函求極值同理，對系數(shù)向量α求偏導，可得簡化的方程組

α=B-1v

(15)

通過式(15)求解出方程組中的系數(shù)向量α后，便可獲得點K處電場分量在x、y和z三個方向上的梯度值

(16)

最后，結(jié)合式(16)和式(3)便可計算出點K處磁場的各個分量。

2 數(shù)值對比與分析

2.1 主程序驗證

為了確保后處理算法對比的可靠性，首先需要驗證有限元數(shù)值模擬主程序的可靠性。設(shè)計如圖4所示的含一低阻(10Ω·m )薄層異常體的層狀模型。水平電偶極源的中心位于坐標軸原點，沿x軸布設(shè)，長度為1m；供電電流I=1A，計算時考慮位移電流；發(fā)射頻率為128Hz。將二次場的解析解作為參考解。

圖4 一維低阻層狀模型

采用矩形六面體網(wǎng)格對計算區(qū)域進行剖分；對目標區(qū)域進行均勻剖分，邊界區(qū)域剖分網(wǎng)格逐漸變大。剖分的具體方案為：x方向[-20000m,20000m]，剖分為105層；y方向[-10000m,10000m]，剖分為42層；z方向[-10000m,5000m]，剖分為34層，其中空氣部分被劃分為11層。

圖5 二次電場分量128Hz時(左)和(右)的幅值相對誤差

2.2 后處理算法數(shù)值對比分析

2.2.1 低阻層狀異常體模型的二次場分析

采用圖4所示低阻層狀異常體模型，發(fā)射頻率取128Hz，測線為y=1500m。

針對此低阻層狀模型分別利用四種后處理算法，對二次磁場分量實部與虛部的相對誤差進行分析，結(jié)果見圖6。表1為單條測線的精度、計算時間和計算內(nèi)存對比。由圖6和表1可知，SFD法和MLSI法的相對誤差較大，最高達26%，整體相對誤差波動非常大。相對這兩種方法，LI法的相對誤差較小，但是在中垂線附近較大，最大相對誤差約為14%。經(jīng)SPR處理的計算結(jié)果相對誤差非常小且很穩(wěn)定，最大相對誤差僅為2.67%，各項平均相對誤差均在1.8%以內(nèi)，基本沒有放大二次電場的誤差，恢復所得二次磁場在中垂線附近效果良好。從平均相對誤差來看，二次磁場的結(jié)果優(yōu)于二次電場。從表1可知，SPR法在計算時間和計算內(nèi)存的需求上都高于其他三種方法，但是對于幾分鐘的計算時間和幾十GB內(nèi)存消耗的有限元正演模擬，這樣的代價是可接受的。

表1 低阻層狀模型四種后處理方法單測線相關(guān)參數(shù)統(tǒng)計

圖6 低阻層狀模型不同后處理算法二次磁場分量相對誤差曲線(y=1500m，f=128Hz)

另外，SFD法是直接由單元形函數(shù)求取節(jié)點上的電場梯度，而SPR后處理技術(shù)是對單元形函數(shù)在高斯點上的電場梯度進行最小二乘擬合，恢復所得節(jié)點電場梯度與高斯點上具有相同的精度級別，因而SPR后處理技術(shù)恢復所得磁場的精度遠高于SFD方法，說明高斯點上的電場梯度值的確比節(jié)點上值精度高，從這個角度也印證超收斂現(xiàn)象。

2.2.2 高阻層狀異常體模型二次場分析

本節(jié)采用的高阻層狀異常體模型與圖4模型類似，區(qū)別是異常層的電阻率為高阻(1000Ω·m)，發(fā)射頻率為256Hz，測線仍取y=1500m。

表2 高阻層狀模型四種后處理方法二次磁場分量平均相對誤差統(tǒng)計

圖7 高阻層狀模型不同后處理算法二次磁場分量相對誤差(y=1500m，f=256Hz)

2.2.3 均勻半空間模型的磁場總場分析

本節(jié)采用與圖4類似的均勻半空間模型，發(fā)射頻率f=256Hz，測線取y=1300m。

圖8為均勻半空間模型四種后處理算法得到的磁場總場分量的實部與虛部相對誤差，表3為單條測線對應(yīng)參數(shù)的統(tǒng)計數(shù)據(jù)。由圖8和表3可知，經(jīng)SFD法處理后的Hx和Hy實部的相對誤差分別為7.25%和9.90%，對應(yīng)的虛部相對誤差小于3%，單點的最大相對誤差達到19%。經(jīng)MLSI后處理的各項平均相對誤差均在10%以內(nèi)，最大相對誤差約為16%。經(jīng)LI法處理的數(shù)據(jù)相對誤差較小，各項相對誤差均低于2.3%。與LI法相比，經(jīng)SPR處理的數(shù)據(jù)相對誤差更小，平均相對誤差最大值為1.55%，且誤差非常穩(wěn)定。

圖8 高阻層狀模型不同后處理算法磁場總場分量相對誤差圖(y=1500m，f=256Hz)

表3 均勻半空間模型四種后處理方法磁場總場分量平均相對誤差統(tǒng)計

2.2.4 三維高阻體模型的二次場分析

為了說明本文引入的SPR后處理對三維模型模擬的有效性，建立了圖9所示高阻體(1000Ω·m)三維模型。背景模型參數(shù)同圖4所示模型。測線為y=1400m，發(fā)射頻率f=256Hz。

圖9 三維高阻體模型

對三維體空間網(wǎng)格進行加密，使得到的數(shù)值解無限趨近于真實解，將收斂的數(shù)值解作為近似解析解(AAS)。這里僅分析磁場，以LI法計算所得磁場數(shù)據(jù)作為磁場的近似解析解。本文在網(wǎng)格相對稀疏的情況下進行后處理方法的對比分析。

將異常體剖分為312×202×186個網(wǎng)格單元，以LI法計算的磁場二次場作為近似解析解。對于稀疏的剖分網(wǎng)格(196×156×59)，不同方法后處理結(jié)果如圖10所示。可以看出，與其他三種后處理方法結(jié)果相比，SPR后處理所得磁場二次場曲線光滑性更好，不存在任何突變點，且總體上與近似解析解一致性更高。

圖10 三維高阻體模型不同后處理方法磁場二次場分量(左)和(右)幅值曲線(y=1400m，f=256Hz)

3 結(jié)論

本文引入基于超收斂單元片梯度值恢復(SPR)的后處理技術(shù)，該算法充分利用有限元解的梯度值在高斯點上具有至少高一階精度的超收斂現(xiàn)象，使用高斯點上的有限元解的梯度值進行最小二乘曲面擬合，恢復得到高精度的磁場。以結(jié)構(gòu)化網(wǎng)格的三維可控源電磁法節(jié)點有限元正演為應(yīng)用基礎(chǔ)，與SFD、MLSI法和LI法相比，SPR后處理技術(shù)在幾乎不增加正演數(shù)值模擬時間和內(nèi)存的情況下，顯著提高了輔助場的精度，且穩(wěn)定性強。

超收斂現(xiàn)象同樣存在于非結(jié)構(gòu)化網(wǎng)格節(jié)點有限元，因而SPR后處理技術(shù)可進一步推廣應(yīng)用到非結(jié)構(gòu)化網(wǎng)格節(jié)點有限元正演，恢復結(jié)果的精度與采樣點上有限元解梯度值的收斂階數(shù)、網(wǎng)格形狀和形函數(shù)階數(shù)等有關(guān)。由于超收斂性是基于節(jié)點有限元的，因而理論上SPR技術(shù)不適用于矢量有限元正演后處理。