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

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

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

0 引言

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

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

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

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

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

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

1 SPR方法原理

1.1 CSEM滿足的邊值問題

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

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

(1)

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

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

Es=0 ∈Γ

(2)

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

1.2 有限元后處理算法

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

(3)

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

1.2.1 SPR后處理技術

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

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

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

圖1 二維線性插值單元片

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

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

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

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

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

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

?E*=Pm

(4)

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

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

(5)

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

m=A-1k

(6)

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

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

在剖分區(qū)域內采用矩形六面體網格進行離散，母單元與子單元坐標對應關系為

(7)

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

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

(8)

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

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

(9)

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

(10)

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

1.2.3 拉格朗日插值(LI)法

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

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

(11)

1.2.4 移動最小二乘(MLSI)法

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

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

設局部域上電場分量的線性插值函數為

f=Uα

(12)

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

設加權函數隨與點K的距離增大而單調降低

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

(13)

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

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

(14)

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

α=B-1v

(15)

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

(16)

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

2 數值對比與分析

2.1 主程序驗證

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

圖4 一維低阻層狀模型

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

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

2.2 后處理算法數值對比分析

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

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

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

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

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

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

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為單條測線對應參數的統(tǒng)計數據。由圖8和表3可知，經SFD法處理后的Hx和Hy實部的相對誤差分別為7.25%和9.90%，對應的虛部相對誤差小于3%，單點的最大相對誤差達到19%。經MLSI后處理的各項平均相對誤差均在10%以內，最大相對誤差約為16%。經LI法處理的數據相對誤差較小，各項相對誤差均低于2.3%。與LI法相比，經SPR處理的數據相對誤差更小，平均相對誤差最大值為1.55%，且誤差非常穩(wěn)定。

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

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

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

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

圖9 三維高阻體模型

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

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

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

3 結論

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

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