乃國(guó)茹, 王緒本, 秦 策, 謝卓良, 陳先潔

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

0 引言

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

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

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

1 基本理論

1.1 大地電磁的邊值問題

考慮電性參數(shù)沿走向不變的二維地電模型，則在二維地電結(jié)構(gòu)中，大地電磁場(chǎng)滿足偏微分方程：

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

(1)

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

1.2 加權(quán)余量法

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

(2)

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

(3)

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

(4)

1.3 二次插值有限元分析

如圖1所示，假設(shè)O點(diǎn)為三角形的重心，取三角單元內(nèi)二次插值形函數(shù)向量為：

圖1 三角單元節(jié)點(diǎn)編碼及面積坐標(biāo)Fig.1 Triangular element node coding and area coordinates

(5)

其中形函數(shù)Ni與面積坐標(biāo)Li的關(guān)系為

N1=(2L1-1)L1

N2=(2L2-1)L2

N3=(2L3-1)L3

N4=4L1L2

N5=4L2L3

N6=4L1L3

假設(shè)u是三角單元內(nèi)任意點(diǎn)的電磁場(chǎng)值，則

(6)

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

(K1-K2+K3)ue=0

(7)

其中，

根據(jù)有限單元法的基本法則，進(jìn)行單元矩陣的擴(kuò)展，可得總體剛度矩陣：

(8)

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

Ku=P

(9)

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

1.4 后驗(yàn)誤差估計(jì)

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

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

(10)

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

(11)

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

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

(12)

式中：TE模式時(shí)k=1；TM模式時(shí)，k=1/σ，wh為對(duì)偶問題的解。

1.5 自適應(yīng)網(wǎng)格剖分

這里采用開源軟件Gmsh建立非結(jié)構(gòu)化三角網(wǎng)格，該網(wǎng)格的靈活度較高，可以在場(chǎng)變化劇烈的區(qū)域及電性突變界面處的加密網(wǎng)格，對(duì)于復(fù)雜地電模型有較好的模擬效果。為了更好地符合邊界條件，在研究區(qū)域以外設(shè)定網(wǎng)格稀疏的擴(kuò)展區(qū)域，如圖2所示。

圖2 非結(jié)構(gòu)化自適應(yīng)網(wǎng)格Fig.2 Unstructured adaptive grid

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

2 算例

2.1 一維層狀介質(zhì)

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

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

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

表1 誤差統(tǒng)計(jì)表

使用電阻率相對(duì)誤差百分比和相位誤差絕對(duì)值來評(píng)估精度，其計(jì)算公式如下：

(13)

(14)

式中：ρs為解析解視電阻率值；ρa(bǔ)為有限元解的視電阻率值；φs為解析解阻抗相位；φa為有限元解的阻抗相位；n為頻點(diǎn)個(gè)數(shù)。分別選擇矩形網(wǎng)格與非結(jié)構(gòu)化三角網(wǎng)格進(jìn)行數(shù)值模擬，計(jì)算兩種模式下有限元解與數(shù)值解的誤差(表1)。

在網(wǎng)格節(jié)點(diǎn)數(shù)量大致相同的情況下，通過表1可知，2D非結(jié)構(gòu)化三角網(wǎng)格FEM數(shù)值解的視電阻率的誤差比矩形網(wǎng)格小，且在計(jì)算速度上得到了一定的提升。

2.2 COMMEMI-2D1模型

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

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

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

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

2.3 起伏地電模型

2.3.1 地壘模型

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

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

圖8 地壘模型自適應(yīng)網(wǎng)格Fig.8 Adaptive grid for basement model

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

2.3.2 地塹模型

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

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

圖11 地塹模型自適應(yīng)網(wǎng)格Fig.11 Adaptive grid of graben model

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

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

2.4 復(fù)雜構(gòu)造模型

2.4.1 組合異常體模型

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

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

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

圖14 組合異常體模型自適應(yīng)網(wǎng)格Fig.14 Combined anomaly model adaptive grid

圖15 組合異常體模型自適應(yīng)有限元二維正演響應(yīng)剖面圖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 斷層模型

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

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

圖17 斷層模型自適應(yīng)網(wǎng)格Fig.17 Fault model adaptive grid

圖18是TE、TM兩種極化模式下斷層模型的二維正演響應(yīng)剖面圖。從圖18(a)、圖18(c)可得出，對(duì)斷層上盤的K型地電模型的三層地層有較為清楚地分辨。高頻部分?jǐn)鄬悠扑閹幊霈F(xiàn)等值線錯(cuò)斷。上盤的表層厚度反應(yīng)明顯比下盤厚度薄，反應(yīng)上盤地層有相對(duì)下盤地層而言有抬升；低頻時(shí)， TM極化模式下的視電阻率圖中等值線出現(xiàn)錯(cuò)動(dòng)。從圖18(b)、圖18(d)可得出，高頻時(shí)，對(duì)淺表及中部地層的抬升分辨較高，底界面埋深差別明顯分辨，左高右低；低頻時(shí)，TE極化模式下的等值線變化不明顯，TM極化模式下的等值線在斷裂處出現(xiàn)扭動(dòng)。

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

3 結(jié)論

通過Gmsh剖分的非結(jié)構(gòu)化三角網(wǎng)格作為初始網(wǎng)格，應(yīng)用對(duì)偶誤差估計(jì)加權(quán)法對(duì)初始網(wǎng)格進(jìn)行自適應(yīng)細(xì)化，構(gòu)建了一維層狀模型、COMMEMI-2D1模型與起伏地電模型，通過解析解、2D 非結(jié)構(gòu)化三角網(wǎng)格FEM解與2D矩形網(wǎng)格FEM數(shù)值解的試算對(duì)比分析，得到以下結(jié)論：

1)通過一維層狀模型，驗(yàn)證了本算法的準(zhǔn)確性。

2) 通過與2D矩形網(wǎng)格FEM數(shù)值解地對(duì)比，驗(yàn)證了本文算法有較高精度。

3)通過復(fù)雜模型與起伏地電模型，計(jì)算并分析了兩種模型的二維MT響應(yīng)特征，可作為實(shí)測(cè)MT數(shù)據(jù)質(zhì)量評(píng)價(jià)、靜態(tài)效應(yīng)研究的參考資料，進(jìn)一步表明本文算法存在使用價(jià)值。

4)自適應(yīng)網(wǎng)格加密相對(duì)全局網(wǎng)格加密能夠節(jié)省計(jì)算量，非結(jié)構(gòu)化網(wǎng)格相比矩形網(wǎng)格能夠更好地適應(yīng)起伏地形和復(fù)雜異常體。