郭家松， 秦 策， 乃國(guó)茹, 謝卓良， 馮 凱

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

0 引言

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

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

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

大地電磁具有趨膚效應(yīng)，高頻電磁波主要受淺部介質(zhì)影響，低頻主要受深部介質(zhì)影響[9]。因而大地電磁正演問(wèn)題過(guò)程中不同頻率對(duì)網(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)格的自動(dòng)加密，優(yōu)化網(wǎng)格，以此提高大地電磁正演計(jì)算精度。

1 正演算法

1.1 邊值問(wèn)題

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

▽×E=iωμH

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

(1)

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

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

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

(2)

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

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

圖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)有限單元法

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

(4)

求解上述變分問(wèn)題，需要對(duì)研究區(qū)域進(jìn)行剖分。本次研究中，借助三維有限元網(wǎng)格生成器-Gmsh將二維研究區(qū)域剖分為非結(jié)構(gòu)化的四邊形網(wǎng)格(圖2)，節(jié)點(diǎn)和四邊形單元遍及整個(gè)研究區(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)場(chǎng)值組成的向量；K1e、K2e以及K3e是單元節(jié)點(diǎn)處的形函數(shù)組成矩陣；u是研究域所有節(jié)點(diǎn)場(chǎng)值形成的向量；Ke是將所有單元上求得的K1e、K2e以及K3e擴(kuò)展為由全體節(jié)點(diǎn)所組成的矩陣。

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

Κu=P

(6)

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

TE模式：

(7)

TM模式：

(8)

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

(9)

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

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。采用本文算法對(duì)該模型進(jìn)行模擬計(jì)算，得到數(shù)值解與一維解析解的結(jié)果如圖3所示。

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

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

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

圖4 時(shí)不同細(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時(shí)的網(wǎng)格自適應(yīng)加密結(jié)果。圖7顯示在斷層破碎帶以及層位分界處網(wǎng)格被明顯加密。

圖8(a)、圖8(b)分別是TE、TM兩種極化模式下斷層破碎帶模型的視電阻率-頻率擬斷面圖。從圖8中可以看出，K型地電模型的三層電性層可以從結(jié)果圖中清晰分辨，且電阻率與實(shí)際地層電阻基本一致，計(jì)算精度高。高頻部分?jǐn)鄬悠扑閹幊霈F(xiàn)等值線錯(cuò)斷。上盤(pán)的表層厚度反應(yīng)明顯比下盤(pán)厚度薄，反應(yīng)上盤(pán)地層有相對(duì)下盤(pán)地層而言有抬升。低頻時(shí)，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 對(duì)稱褶皺模型示意圖Fig.10 Schematic diagram of symmetry fold model

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

2.3 對(duì)稱褶皺模型

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

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

圖13(a)、圖13(b)分別顯示的是TE、TM兩個(gè)極化模式對(duì)稱褶皺模型視電阻率-頻率擬斷面圖。 圖13顯示高頻時(shí)， 可清晰地反應(yīng)地層起伏。頻率在100 Hz，視電阻率值等值線在褶皺處斷裂，褶皺處電阻率值較兩側(cè)地層低。低頻時(shí)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模式擬斷面圖在褶皺波谷處顯示為低阻圈閉異常。低頻時(shí)阻抗相位在波谷處顯示為地層的隆起。TE極化模式低頻時(shí)并未表現(xiàn)起伏變化。

3 結(jié)論

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