基于非結(jié)構(gòu)化網(wǎng)格的大地電磁2D自適應(yīng)有限元正演模擬

郭家松， 秦 策， 乃國茹, 謝卓良， 馮 凱

(成都理工大學(xué) 地球物理學(xué)院，成都 610059; 2.河南理工大學(xué) 物理與電子信息學(xué)院，焦作 454000)

0 引言

20世紀(jì)50年代，蘇聯(lián)學(xué)者吉洪諾夫(Tikhonov)[1]以及法國學(xué)者卡尼爾(L.Cagniard)[2]提出了大地電磁測深法的理論基礎(chǔ)。該方法以天然交變電磁場為場源，基于電磁感應(yīng)原理，研究地下介質(zhì)在地表的電磁場分布規(guī)律，以此解決相關(guān)的地質(zhì)問題。經(jīng)過幾十年的發(fā)展，大地電磁測深法因其成本低廉，對低阻反應(yīng)敏感、勘探深度大等優(yōu)點(diǎn)被廣泛用于深部構(gòu)造研究，油氣勘探、地下水勘探與環(huán)境保護(hù)等方面[3-5]。

大地電磁測深法的最終目的是得到與實(shí)際地下介質(zhì)分布情況最接近的地電模型。這個過程需對研究區(qū)的實(shí)測資料進(jìn)行一系列正反演計(jì)算，使模型響應(yīng)與實(shí)際資料達(dá)到最佳擬合。正演是反演的基礎(chǔ)，大地電磁測深法的正演建立在Maxwell方程組之上。目前用于大地電磁測深法正演數(shù)值模擬的方法主要是有限差分法、有限單元法和積分方程法[6-7]。

這里著重研究有限單元法求解大地電磁正演問題。有限單元法求解變分問題的核心思想是將計(jì)算域劃分為有限個互不重疊的單元從而在單元上進(jìn)行插值求解[8]，求解域的剖分結(jié)果對計(jì)算結(jié)果有很大影響。在大地電磁測深法的正演問題中，采用規(guī)則化網(wǎng)格進(jìn)行研究區(qū)域解的剖分局限性很大，剖分結(jié)構(gòu)對復(fù)雜構(gòu)造的地電模型擬合程度差，易產(chǎn)生較大的幾何離散誤差。筆者采用非結(jié)構(gòu)化的四邊形網(wǎng)格進(jìn)行剖分，從而使剖分結(jié)果與實(shí)際模型能達(dá)到最大程度的擬合。

大地電磁具有趨膚效應(yīng)，高頻電磁波主要受淺部介質(zhì)影響，低頻主要受深部介質(zhì)影響[9]。因而大地電磁正演問題過程中不同頻率對網(wǎng)格的稀疏程度的要求不同，利用同一套網(wǎng)格求解不同頻率下的大地電磁響應(yīng)必然有誤差。因此，在非結(jié)構(gòu)化四邊形網(wǎng)格剖分基礎(chǔ)上引入自適應(yīng)算法，從求解域的粗網(wǎng)格出發(fā)，在不同頻點(diǎn)下利用網(wǎng)格單元的后驗(yàn)誤差估計(jì)指導(dǎo)網(wǎng)格的自動加密，優(yōu)化網(wǎng)格，以此提高大地電磁正演計(jì)算精度。

1 正演算法

1.1 邊值問題

Maxwell方程組描述了電場、磁場以及電流密度、電荷密度之間的關(guān)系，是電磁理論的核心。大地電磁測深正演理論基于Maxwell方程組，從方程組出發(fā)，有諧頻率為ω的定態(tài)電磁場方程為：

▽×E=iωμH

▽×H=(σ-iωε)E

(1)

式中:μ和σ分別表示介質(zhì)的磁導(dǎo)率電導(dǎo)率；ε為介電常數(shù)。

這里主要針對二維地電結(jié)構(gòu)，即介質(zhì)的電性參數(shù)在有明顯走向的傾斜巖層、背向斜等地質(zhì)構(gòu)造中只沿傾向和深度方向變化。因此，取x軸平行構(gòu)造走向，y軸垂直于x軸，z軸垂直水平面建立笛卡爾坐標(biāo)系，將電場和磁場沿坐標(biāo)軸進(jìn)行分解。根據(jù)上述定態(tài)電磁場方程，可得二維介質(zhì)電場分量與磁場分量所滿足的偏微分方程：

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

(2)

其中，對于TE極化模式，u=Ex，τ=1/iωμ，λ=σ-iωμ；對于TM極化模式，u=Hx，τ=1/iωμ，λ=σ-iωμ。

與一般求解偏微分方程一樣，上述偏微分方程需要在特定的邊界條件下求解。考慮如圖1所示的研究區(qū)域，其中TE極化模式下的研究區(qū)域?yàn)锳BCD，TM極化模式下的研究區(qū)域?yàn)锳′B′CD(2)式引入邊界條件后可形成如下邊值問題：

圖1 研究區(qū)域示意圖Fig.1 Schematic diagram of the study area

(3)

圖2 非結(jié)構(gòu)化四邊形網(wǎng)格剖分示意圖Fig.2 Schematic diagram of unstructured quadrilateral mesh

1.2 自適應(yīng)有限單元法

有限元是將邊值問題轉(zhuǎn)化為等價的變分問題進(jìn)而求解。變分問題如下：

(4)

求解上述變分問題，需要對研究區(qū)域進(jìn)行剖分。本次研究中，借助三維有限元網(wǎng)格生成器-Gmsh將二維研究區(qū)域剖分為非結(jié)構(gòu)化的四邊形網(wǎng)格(圖2)，節(jié)點(diǎn)和四邊形單元遍及整個研究區(qū)域。在單元內(nèi)采用二次插值，將連續(xù)函數(shù)的求解轉(zhuǎn)變?yōu)楣?jié)點(diǎn)上的離散函數(shù)值的求解，進(jìn)而將求解域離散為所有單元的線性組合，并利用網(wǎng)格單元內(nèi)的二次插值函數(shù)可使變分表達(dá)式轉(zhuǎn)變?yōu)橄蛄颗c矩陣相乘的形式為式(2)。

(5)

其中:ue是網(wǎng)格單元節(jié)點(diǎn)場值組成的向量；K1e、K2e以及K3e是單元節(jié)點(diǎn)處的形函數(shù)組成矩陣；u是研究域所有節(jié)點(diǎn)場值形成的向量；Ke是將所有單元上求得的K1e、K2e以及K3e擴(kuò)展為由全體節(jié)點(diǎn)所組成的矩陣。

對式(5)求變分并令其等于零，并加入變分問題的邊界條件，可得如下線性方程組

Κu=P

(6)

求解上述線性方程采用LU分解法求得上述線性方程組，從而得到各節(jié)點(diǎn)的電磁場值，利用數(shù)值算法求得場值沿垂向的偏導(dǎo)數(shù)?u/?z得到輔助場，將輔助場代入式(7)、式(8)求得視電阻率和阻抗相位響應(yīng)。

TE模式：

(7)

TM模式：

(8)

為滿足網(wǎng)格在不同頻率下的計(jì)算精度，在有限單元法求解大地電磁正演問題的過程中引入自適應(yīng)原則。在每個頻點(diǎn)下通過單元網(wǎng)格的后驗(yàn)誤差估計(jì)值去指導(dǎo)網(wǎng)格的加密。本文的誤差計(jì)算方法基于Kelly等[10]提出的有限單元后驗(yàn)誤差估計(jì)方法，通過計(jì)算單元邊界上場值梯度跳變值沿著單元邊界的積分來近似作為每個單元的誤差得到如下誤差計(jì)算式：

(9)

式中:cF=hf/2，hf表示是單元K中最長的對角線；〔·〕表示單元K邊界上場值梯度的跳變值。

2 模型試算

2.1 層狀模型

為檢驗(yàn)本文方法的有效性，構(gòu)建有三層層狀模型，模型的地電結(jié)構(gòu)為K型，參數(shù)如下：ρ1=10 Ω·m，h1=1 000 m，ρ2=100 Ω·m，h2=1 200 m，ρ3=10 Ω·m。采用本文算法對該模型進(jìn)行模擬計(jì)算，得到數(shù)值解與一維解析解的結(jié)果如圖3所示。

表1 解析解計(jì)算誤差Tab.1 Calculation error of analytical solution

從表1計(jì)算的解析解誤差值中可以看出，網(wǎng)格剖分的稀疏程度影響著大地電磁有限元正演模擬的精度，細(xì)化次數(shù)越多，誤差越小，有限元正演模擬的精度越高。圖3為TE極化模式下求得的數(shù)值解和一維解析解的對比，結(jié)果顯示網(wǎng)格細(xì)化程度越高，解析解和數(shù)值解逐漸吻合。但網(wǎng)格剖分較粗糙時，高頻段計(jì)算出的大地電磁響應(yīng)與解析解有很大誤差，說明粗網(wǎng)格不能滿足高頻時的計(jì)算精度要求。大地電磁波具有趨膚效應(yīng)，對研究域剖分結(jié)果的縱向網(wǎng)格的尺寸與電磁波的最小趨膚深度有關(guān)，因此在高頻時需要采用較小網(wǎng)格進(jìn)行計(jì)算[11-12]。圖4即為100 Hz時的網(wǎng)格細(xì)化結(jié)果。結(jié)果顯示在該頻點(diǎn)下，網(wǎng)格細(xì)化主要集中在淺部，隨著細(xì)化次數(shù)的增加，靠近地表的網(wǎng)格越密，尺寸越小，計(jì)算精確度有了明顯提高。通過以上分析，利用本文的自適應(yīng)有限元算法對網(wǎng)格進(jìn)行局部加密，可以滿足不同頻率下的網(wǎng)格的疏密程度，且在對初始網(wǎng)格細(xì)化兩次，數(shù)值解和解析解基本吻合，且誤差較低，證實(shí)了本文算法的有效性。

圖3 解析解和數(shù)值解的對比Fig.3 Comparison between analytic solution and numerical solution(a)視電阻率；(b)阻抗相位

圖4 時不同細(xì)化次數(shù)下的網(wǎng)格細(xì)化結(jié)果Fig.4 Mesh refinement results for different refinement times at f=100 Hz(a)N=0；(b)N=1；(c)N=2

圖5 斷層破碎帶模型示意圖Fig.5 Schematic diagram of the fault fracture zone

圖6 初始網(wǎng)格剖分Fig.6 Initial mesh generation

圖7 f=100 Hz網(wǎng)格自適應(yīng)結(jié)果Fig.7 Meshes adaptive result at 100 Hz

2.2 隱伏斷層破碎帶模型

構(gòu)建圖5所示逆斷層破碎帶模型，電性參數(shù)如圖中所示。在水平方向3 000 m～4 000 m范圍內(nèi)有一電阻率為10 Ω.m的低阻破碎帶，破碎帶寬度為150 m，傾角為60°。

圖6和圖7分別顯示的模型的初始網(wǎng)格剖分結(jié)果以及頻率f=100 Hz時的網(wǎng)格自適應(yīng)加密結(jié)果。圖7顯示在斷層破碎帶以及層位分界處網(wǎng)格被明顯加密。

圖8(a)、圖8(b)分別是TE、TM兩種極化模式下斷層破碎帶模型的視電阻率-頻率擬斷面圖。從圖8中可以看出，K型地電模型的三層電性層可以從結(jié)果圖中清晰分辨，且電阻率與實(shí)際地層電阻基本一致，計(jì)算精度高。高頻部分?jǐn)鄬悠扑閹幊霈F(xiàn)等值線錯斷。上盤的表層厚度反應(yīng)明顯比下盤厚度薄，反應(yīng)上盤地層有相對下盤地層而言有抬升。低頻時，TE模式破碎帶下方地層顯示為凹陷形態(tài)。TM極化模式斷面圖在破碎帶下方電阻率值較兩側(cè)低，形成縱向的帶狀低阻體，基底在破碎帶下方被低阻體切斷，形成斷裂。

圖8 視電阻率-頻率擬斷面圖Fig.8 Apparent resistivity-frequency pseudosection(a)TE; (b)TM

圖9 阻抗相位-頻率擬斷面圖Fig.9 Impedance phase-frequency pseudosections(a)TE; (b)TM

圖10 對稱褶皺模型示意圖Fig.10 Schematic diagram of symmetry fold model

圖9(a)、圖9(b)分別是TE、TM兩種極化模式斷層破碎帶模型的阻抗相位-頻率擬斷面圖。從圖9中可以看出，高頻時 可明顯看出淺表地層底界面埋深有差別，左高右低。在頻率為范圍內(nèi)，相位剖面顯示為低異常，TM極化模式的等值線出現(xiàn)斷裂。低頻時，破碎帶底部相位等值線呈現(xiàn)凹陷形態(tài)，TM模式為凸起，但反應(yīng)不明顯。

2.3 對稱褶皺模型

建立圖10所示的對稱褶皺模型，模型劃分有四層電性層，為HK型。電性層層厚及電阻率參數(shù)及褶皺的大小如圖10所示，其中對稱褶皺存在于水平方向3 000 m～6 000 m的范圍內(nèi)。以波長和波幅定義對稱褶皺構(gòu)造的大小，其中褶皺波長為2 000 m，波峰的幅度為500 m。

圖11和圖12分別顯示的對稱褶皺模型的初始網(wǎng)格剖分結(jié)果以及當(dāng)頻率f=100 Hz時，對初始網(wǎng)格細(xì)化兩次之后的自適應(yīng)加密結(jié)果。通過細(xì)化前后的網(wǎng)格剖分結(jié)果可以看出，細(xì)化兩次之后在淺部，網(wǎng)格被明顯細(xì)化加密。

圖13(a)、圖13(b)分別顯示的是TE、TM兩個極化模式對稱褶皺模型視電阻率-頻率擬斷面圖。 圖13顯示高頻時， 可清晰地反應(yīng)地層起伏。頻率在100 Hz，視電阻率值等值線在褶皺處斷裂，褶皺處電阻率值較兩側(cè)地層低。低頻時TE模式剖面顯示地層并未有明顯的起伏變化，而TM模式剖面圖清晰顯示存在地層起伏。且隨著頻率降低，褶皺的波谷處凹陷形態(tài)反應(yīng)越加明顯，且波谷波幅隨頻率變低增大。

圖11 模型初始網(wǎng)格剖分Fig.11 Initial mesh generation

圖12 f=100 Hz網(wǎng)格自適應(yīng)結(jié)果Fig.12 Meshes adaptive result at 100 Hz

圖13 視電阻率-頻率擬斷面圖Fig.13 Apparent resistivity-frequency pseudosection(a)TE; (b)TM

圖14 阻抗相位-頻率擬斷面圖Fig.14 Impedance phase-frequency pseudosections(a)TE; (b)TM

圖14(a)、圖14(b)分別顯示的是阻抗相位-頻率擬斷面圖。高頻部分，地層有明顯起伏，形態(tài)與模型褶皺形態(tài)吻合。但頻率越低，起伏幅度越小。在頻率范圍內(nèi)，同一頻點(diǎn)相位值出現(xiàn)斷裂。頻率范圍為1 Hz～10 Hz，TM模式擬斷面圖在褶皺波谷處顯示為低阻圈閉異常。低頻時阻抗相位在波谷處顯示為地層的隆起。TE極化模式低頻時并未表現(xiàn)起伏變化。

3 結(jié)論

采用了三維網(wǎng)格剖分器-Gmsh對研究模型進(jìn)行非結(jié)構(gòu)化四邊形網(wǎng)格剖分，可以使剖分結(jié)果與實(shí)際模型達(dá)到最大程度上的擬合，能夠真實(shí)的模擬地形起伏，斷層構(gòu)造，褶皺等復(fù)雜地質(zhì)結(jié)構(gòu)模型。采用自適應(yīng)的有限元方法求解大地電磁正演問題，在每個頻點(diǎn)下利用網(wǎng)格單元的后驗(yàn)誤差估計(jì)指導(dǎo)網(wǎng)格的局部加密，滿足了網(wǎng)格對不同頻點(diǎn)計(jì)算精度的要求，進(jìn)一步優(yōu)化了網(wǎng)格質(zhì)量和數(shù)量，避免了人為網(wǎng)格剖分所造成的影響，從而提高了二維正演計(jì)算精度。