基于非結構化三角網格的大地電磁二維h-型自適應有限元正演模擬

乃國茹, 王緒本, 秦 策, 謝卓良, 陳先潔

(1.成都理工大學 地球物理學院，成都 610059；2.河南理工大學，焦作 454002)

0 引言

經過幾代學者不懈地努力，大地電磁測深法得到了蓬勃發展，因其具有勘探深度大、施工方便、成本低、對低阻體反應靈敏等優點，在地球物理勘探、地質勘查和工程探測等眾多領域都發揮著重要作用。大地電磁測深法的最終目的，是得到與實際實際地質情況最為吻合的地電模型。該過程要求對實測資料進行正反演計算，使模型響應與實際資料達到最佳擬合。正演是反演的基礎，因此，對復雜地質模型進行高精度快速正演計算是MT數值模擬研究的重點[1-3]。

二維MT的正演問題已相對成熟，如有限差分、積分方程、有限單元等方法均有穩定實現[4-7]。這里著重探討有限單元法求解大地電磁正演問題。Rodi[8]、Wannamaker等[9]使用規則化矩形網格剖分對大地電磁二維正演模擬進行了研究；徐世浙等[10-11]將原來簡單的網格剖分發展為三角單元剖分和三角單元-矩形單元剖分，很大程度提高了大地電磁場場值求解的精度和計算速度。但復雜的地電模型中，傳統結構化網格為了保證內節點拓撲結構的一致性，往往需要更多的網格。Kerry Key等[12]、柳建新等[13]利用非結構化自適應有限元網格對二維大地電磁進行研究；Ren等[14]提出了面向目標的自適應矢量有限元法，大大提高了計算精度和適應性；羅天涯等[15]、楊振武等[16]相繼利用非結構化網格進行了二維MT正演研究。筆者認為，將非結構化三角網格與h-型自適應算法相結合，一定程度上可提升復雜構造二維正演的模擬精度與計算時間。

筆者主要研究基于非結構化三角網格劃分的自適應有限元解決大地電磁的二維正演問題。有限元法首先通過計算將研究域分成許多小單元結構的每個節點電磁場值，接著每個小單元的場值可利用二次插值函數得到，最后可計算出整個研究域的電磁場值。基于非結構化三角形網格剖分的自適應有限元算法，對不規則與復雜結構模型的邊界有了更好地模擬，為了提升數值模擬計算精度，采用了基于后處理技術的后驗誤差估計方法(Dual Error Estimate Weighting, DEW)[12]，對局部網格實現了自動加密。

1 基本理論

1.1 大地電磁的邊值問題

考慮電性參數沿走向不變的二維地電模型，則在二維地電結構中，大地電磁場滿足偏微分方程：

▽·(τ▽u)+λu=0

(1)

在大地電磁二維正演中，式(1)中電磁波的性質為定態時諧波，記e-iwt為時諧因子。TE模式下，u=Ex,τ=1/iμε,λ=σ-iωε，由于上邊界與地面之間距離足夠遠，取u=1；TM模式下，u=Hx,τ=1/(σ-iωε),λ=iωμ，因空氣中的磁場分量不受地下介質電性分布的影響，且為常量，可取地面上的u=1；對于下邊界與左右兩邊界，TE、TM兩種極化模式具有相同的邊界條件。

1.2 加權余量法

記φ為權重函數，引入加權余量法，式(1)兩邊同時乘以φ并進行積分，可得式(2)。

(2)

針對式(2)左邊的第一項，利用矢量運算式與積分變換式，并將相應邊界條件代入可得式(3)。

(3)

將式(3)代入式(2)可得式(4)。

(4)

1.3 二次插值有限元分析

如圖1所示，假設O點為三角形的重心，取三角單元內二次插值形函數向量為：

圖1 三角單元節點編碼及面積坐標Fig.1 Triangular element node coding and area coordinates

(5)

其中形函數Ni與面積坐標Li的關系為

N1=(2L1-1)L1

N2=(2L2-1)L2

N3=(2L3-1)L3

N4=4L1L2

N5=4L2L3

N6=4L1L3

假設u是三角單元內任意點的電磁場值，則

(6)

將式(6)代入式(4)，消去變分項后得：

(K1-K2+K3)ue=0

(7)

其中，

根據有限單元法的基本法則，進行單元矩陣的擴展，可得總體剛度矩陣：

(8)

代入邊界條件后，可得到有限元方程：

Ku=P

(9)

式中：Κ為對稱矩陣，是由單元上計算的系數矩陣擴展而來的總體系數矩陣，為了節約計算內存，利用定帶寬存儲方式，只將矩陣的上三角或下三角元素進行存儲，并只對非零元素進行存儲。在有限元方程求解過程中，首先采用乘數法對有限元方程添加邊界條件，該方法得優勢在于對方程的改動較小，程序容易實現，且效率較高[17]；然后利用直接稀疏LU分解法對系數矩陣組成的線性方程組進行求解，得到研究區上的節點場值u，有限元求解完成。

1.4 后驗誤差估計

在有限元中，當全局網格滿足一定的條件后，局部網格的密度對于有限元解的精度有著很大影響，這表明相對于全局網格的細化，更好的細化方案是根據需求有選擇性的進行局部區域網格地細化。利用3D-Gmsh進行研究域的粗剖分，將生成的粗網格作為初始網格，然后求解出每個單元節點處的有限元解uh，進而計算出解的梯度值▽uh，得到每個單元內的局部誤差為式(10)。

re=‖T▽uh-▽uh‖2(e)=‖(T-I)▽uh‖2(e)

(10)

式中：T為恢復因子；I為單位算子。其中，恢復因子T的超收斂特性對后驗誤差估計的有效性存在重要影響，故較精確的誤差估計表達式為式(11)。

(11)

其中：u為微分方程的真實解；h為每個三角單元外接圓的直徑。該方法屬于一種全局網格細化的基本誤差估計方法，為了實現對局部網格的加密，需要求解一個對偶問題，對原后驗誤差估計值使用對偶問題的解的后驗誤差估計值進行加權，即對局部誤差re增加一個加權項，多位學者對該問題的定義都已做了詳細的論述[14,18]。加權后驗誤差估計公式為式(12)。

Re=‖k(T-I)▽uh‖2(e)‖(T-I)▽wh‖2(e)

(12)

式中：TE模式時k=1；TM模式時，k=1/σ，wh為對偶問題的解。

1.5 自適應網格剖分

這里采用開源軟件Gmsh建立非結構化三角網格，該網格的靈活度較高，可以在場變化劇烈的區域及電性突變界面處的加密網格，對于復雜地電模型有較好的模擬效果。為了更好地符合邊界條件，在研究區域以外設定網格稀疏的擴展區域，如圖2所示。

圖2 非結構化自適應網格Fig.2 Unstructured adaptive grid

自適應網格細化方式有：h-型、p-型和hp-型三類[19-20]。h-型是通過逐步加密網格，而形函數的階次保持不變的情況下滿足精度要求；p-型是通過改變形函數的階次，但單元網格的大小維持不變的情況下逼近精確解；hp-型是結合h-型和p-型兩種細化方式，雖然細化效果優于前兩種方法，但是實現較為困難。筆者采用h-型自適應網格細化技術，利用初始粗網格得到的局部誤差為指導進行局部區域的網格細化，對細化后的網格再次進行誤差計算，直至細化后的網格誤差滿足精度要求。

2 算例

2.1 一維層狀介質

如圖3所示，構建了一個H型水平層狀地電模型。其中ρ1=200 Ω·m，h1=2 000 m；ρ2= 100Ω·m，h2= 2 000 m；ρ3=300 Ω·m，層厚無限延伸。測線距離L= 20 km，測點個數N= 101個，點距r= 200 m； 頻率f為 10-4Hz～104Hz，取對數(log10)后，以0.2為采樣間隔，頻點n= 41個，自適應加密后網格數目為8 652。對H型地電模型采用本文算法數值模擬，將模擬結果與解析解、2D矩形網格剖分有限元解成圖(圖4)分析。由圖4可知，視電阻率曲線與解析解視電阻率曲線基本吻合，在高頻段，阻抗相位曲線與解析解阻抗相位曲線存在誤差。

圖3 一維層狀介質地電模型Fig.3 One-dimensional layered dielectric geoelectric model

圖4 一維層狀介質響應曲線Fig.4 Response curve of one-dimensional layered media(a)視電阻率曲線;(b)相位曲線

表1 誤差統計表

使用電阻率相對誤差百分比和相位誤差絕對值來評估精度，其計算公式如下：

(13)

(14)

式中：ρs為解析解視電阻率值；ρa為有限元解的視電阻率值；φs為解析解阻抗相位；φa為有限元解的阻抗相位；n為頻點個數。分別選擇矩形網格與非結構化三角網格進行數值模擬，計算兩種模式下有限元解與數值解的誤差(表1)。

在網格節點數量大致相同的情況下，通過表1可知，2D非結構化三角網格FEM數值解的視電阻率的誤差比矩形網格小，且在計算速度上得到了一定的提升。

2.2 COMMEMI-2D1模型

對COMMEMI-2D1模型(圖5)進行數值模擬，旨在驗證該算法對二維構造的有效性。如圖5所示，在電阻率ρ2=100 Ω·m的地下均勻半空間存在一個電阻率ρ1= 0.5 Ω·m的低阻異常體，其寬為1 km，高為2 km，頂面距離地表0.25 km；設置長度L=40 km的測線，測點數N=81，測點間距r=0.5 km。為了與Zhdanov等[21]發布的結果進行對比驗證，使用本文算法選取10 Hz頻率進行數值模擬(圖6)。

圖5 COMMEMI-2D1模型Fig.5 COMMEMI-2D1 model

由圖6可知，本文算法在TE和TM兩種模式的模擬結果與Zhdanov等[21]提供的數據基本吻合，說明了計算方法的正確性。根據響應曲線可知，在低阻異常體的邊界處，TE模式視電阻率產生明顯的突跳增大，而TM模式相對變化平緩，且TE模式視電阻率受低阻體的影響范圍較TM模式更廣，曲線體現出向上的開口更大。

圖6 COMMEMI-2D1模型Fig.6 COMMEMI-2D1 model(a) TE模式;(b) TM模式

2.3 起伏地電模型

2.3.1 地壘模型

如圖7所示，構建一個地壘模型。模型上方電阻率ρAi r= 108Ω·m，下方電阻率ρ= 100 Ω·m，起伏地表水平范圍為-800 m～800 m，高度為500 m，凸起部分范圍為-200 m～200 m，研究區范圍為-3 000 m～3 000 m，在地表設置86個測點。采用非結構化三角形單元進行網格粗剖分，經過10次細化迭代，自適應加密之后的網格單元數為7 640，網格節點為3 467(圖8)。利用自適應有限元法對該模型進行正演計算，選取0.01 Hz、1 Hz、100 Hz 3個頻點，可得兩個模式下的正演響應曲線(圖9)。

圖7 地壘模型示意圖Fig.7 Schematic diagram of the land barrier model

圖8 地壘模型自適應網格Fig.8 Adaptive grid for basement model

地壘模型對圖9(a)、圖9(b)、圖9(c)與圖9(d)中的曲線都產生了不同程度的影響。其中，對于視電阻率曲線而言，兩種模式下的曲線是反向的， TE模式下曲線隨頻率的波動范圍要明顯低于TM模式，且圖9(c)曲線在地形變化處(-800 m、800 m、-200 m和200 m)的突變程度明顯高于圖9(a)曲線； 而對于相位曲線而言，地壘模型對二者的影響很接近，其與頻率的變化大致成正相關，且在地形變化處(-800 m、800 m、-200 m和200 m)曲線的突跳大致相同。

2.3.2 地塹模型

如圖10所示，構建一個地塹模型。模型參數與地壘模型參數相同，采用非結構化三角形單元進行網格粗剖分，經過12細化迭代，自適應加密之后的網格單元數為9 860，網格節點為4 617(圖11)。采用自適應有限元法對該模型進行正演計算，選取0.01 Hz、1 Hz、100 Hz 3個頻點，可得兩個模式下的正演響應曲線(圖12)。

圖10 地塹模型示意圖Fig.10 Schematic diagram of the mantle model

圖11 地塹模型自適應網格Fig.11 Adaptive grid of graben model

由圖12可知，地嵌模型對兩種極化模式下的視電阻率曲線和阻抗相位曲線都產生了一定程度的影響。對于視電阻率曲線，兩種極化模式下的曲線是反向的， TE模式下曲線受地形起伏影響比TM模式更加明顯，且圖12(c)曲線在地形變化處(-800 m、800 m、-200 m和200 m)的突變程度明顯高于圖12(a)曲線； 而對于相位曲線，地塹模型對二者的影響很接近，其與頻率的變化大致成正相關，且在地形變化處(-800 m、800 m、-200 m和200 m)隨著頻率的增大曲線突跳明顯。

圖12 地塹模型自適應有限元二維正演響應曲線Fig.12 Adaptive finite element 2D forward response curve of graben model(a)TE模式視電阻率曲線;(b)TE模式阻抗相位曲線;(c)TM模式視電阻率曲線;(d)TM模式阻抗相位曲線

2.4 復雜構造模型

2.4.1 組合異常體模型

構建如圖13所示的地電模型。在電阻率為100 Ω·m的均勻半空間，存在電阻率為10 Ω·m的低阻塊體和電阻率為1 000 Ω·m的高阻塊體。兩塊體大小相同，距離地面埋深均為1 000 m，兩塊體之間相距4 000 m。

圖13 組合異常體模型示意圖Fig.13 Schematic diagram of combined abnormal body model

采用非結構化三角形單元進行網格粗剖分，經過3次細化迭代，自適應加密之后的網格單元數為6 100，網格節點為2 788(圖14)。采用自適應有限元法對該模型進行正演計算，得到二維正演響應剖由圖15可知，兩種極化模式下的正演結果都很好地刻畫了兩個異常體的基本特征，TM極化模式對異常體的正演響應顯得更加靈敏，而TE極化模式對異常體的分辨率稍差。再者，正演響應結果對低阻異常體的分辨率比較高，而對高阻異常體的分辨率較差。

圖14 組合異常體模型自適應網格Fig.14 Combined anomaly model adaptive grid

圖15 組合異常體模型自適應有限元二維正演響應剖面圖Fig.15 Adaptive finite element two-dimensional forward response profile of combined abnormal body model(a)TE模式視電阻率圖(b)TE模式阻抗相位圖;(c)TM模式視電阻率圖;(d)TM模式阻抗相位圖

2.4.2 斷層模型

構建如圖16所示逆斷層模型，斷層上下盤均為層狀均勻介質，地電結構為KH型，電性參數如圖16所示。斷層傾角為45°，斷層上下盤斷距為500 m。采用非結構化三角形單元進行網格粗剖分，經過3次細化迭代，自適應加密之后的網格單元數為3 497，網格節點為1 430(圖17)。

圖16 斷層模型示意圖Fig.16 Schematic diagram of fault model

圖17 斷層模型自適應網格Fig.17 Fault model adaptive grid

圖18是TE、TM兩種極化模式下斷層模型的二維正演響應剖面圖。從圖18(a)、圖18(c)可得出，對斷層上盤的K型地電模型的三層地層有較為清楚地分辨。高頻部分斷層破碎帶處出現等值線錯斷。上盤的表層厚度反應明顯比下盤厚度薄，反應上盤地層有相對下盤地層而言有抬升；低頻時， TM極化模式下的視電阻率圖中等值線出現錯動。從圖18(b)、圖18(d)可得出，高頻時，對淺表及中部地層的抬升分辨較高，底界面埋深差別明顯分辨，左高右低；低頻時，TE極化模式下的等值線變化不明顯，TM極化模式下的等值線在斷裂處出現扭動。

圖18 斷層模型自適應有限元二維正演響應剖面圖Fig.18 Fault model adaptive finite element 2D forward response profile(a) TE模式視電阻率圖;(b)TE模式阻抗相位圖;(c)TM模式視電阻率圖;(d)TM模式阻抗相位圖

3 結論

通過Gmsh剖分的非結構化三角網格作為初始網格，應用對偶誤差估計加權法對初始網格進行自適應細化，構建了一維層狀模型、COMMEMI-2D1模型與起伏地電模型，通過解析解、2D 非結構化三角網格FEM解與2D矩形網格FEM數值解的試算對比分析，得到以下結論：

1)通過一維層狀模型，驗證了本算法的準確性。

2) 通過與2D矩形網格FEM數值解地對比，驗證了本文算法有較高精度。

3)通過復雜模型與起伏地電模型，計算并分析了兩種模型的二維MT響應特征，可作為實測MT數據質量評價、靜態效應研究的參考資料，進一步表明本文算法存在使用價值。

4)自適應網格加密相對全局網格加密能夠節省計算量，非結構化網格相比矩形網格能夠更好地適應起伏地形和復雜異常體。