李志旋 岳明鑫? 周官群2)

1) (中國(guó)科學(xué)技術(shù)大學(xué)地球和空間科學(xué)學(xué)院, 合肥 230026)

2) (安徽惠洲地質(zhì)安全研究院股份有限公司, 合肥 231202)

(2018 年 8 月 21 日收到; 2018 年 12 月 3 日收到修改稿)

采用矢量有限元法實(shí)現(xiàn)了三維電磁擴(kuò)散場(chǎng)數(shù)值模擬, 并成功將其應(yīng)用在大地電磁的正演研究中. 為靈活精確地?cái)M合起伏地形和地下不規(guī)則構(gòu)造, 采用由不規(guī)則四面體單元組成的非結(jié)構(gòu)化網(wǎng)格, 可根據(jù)模型設(shè)計(jì)的需要調(diào)整網(wǎng)格的大小. 引入了基于二次場(chǎng)理論, 將解析的一次場(chǎng)從總場(chǎng)中扣除, 直接計(jì)算二次場(chǎng), 使得誤差僅局限于相對(duì)較小的二次場(chǎng), 以提高總場(chǎng)計(jì)算精度. 常規(guī)的節(jié)點(diǎn)有限元法不滿足電性分界面上法向電場(chǎng)不連續(xù)和無(wú)源區(qū)單元內(nèi)電流密度無(wú)散, 違反麥克斯韋方程組. 為克服節(jié)點(diǎn)有限元法的弊端, 使用矢量有限元法求解基于二次電場(chǎng)的偏微分方程. 另外, 在算法設(shè)計(jì)中, 考慮了磁導(dǎo)率參數(shù)的變化, 可以模擬磁導(dǎo)率不均勻的模型.通過(guò)與COMMEMI模型已發(fā)表的結(jié)果對(duì)比, 證明了本文算法的正確性和精確性. 為突顯非結(jié)構(gòu)網(wǎng)格優(yōu)勢(shì), 計(jì)算了橢球異常體模型和任意地形模型的MT響應(yīng), 并詳細(xì)討論了地形和磁化效應(yīng)對(duì)三維數(shù)值模擬結(jié)果的影響.

1 引 言

電磁擴(kuò)散場(chǎng)數(shù)值模擬方法中, 積分方程法、有限差分法和有限元法是三種經(jīng)典并廣泛使用的方法. 積分方程法[1?3]受并矢格林函數(shù)的復(fù)雜性限制,僅適用于簡(jiǎn)單規(guī)則模型的模擬. 有限差分法[4?6]的優(yōu)點(diǎn)是易于實(shí)現(xiàn), 其方程的離散形式?jīng)Q定了該方法只能模擬可四邊形或六面體剖分的計(jì)算域, 對(duì)于復(fù)雜模型, 如起伏地形和不規(guī)則地下結(jié)構(gòu), 采取臺(tái)階狀近似誤差大, 且剖分難實(shí)現(xiàn). 相對(duì)而言, 采用非結(jié)構(gòu)化網(wǎng)格時(shí), 有限元法則能夠靈活精確地處理各種復(fù)雜模型[7].

1971年, 有限元法最早由Coggon[8]引入地球物理電磁場(chǎng)數(shù)值模擬. 早期的有限元法使用規(guī)則四邊形和六面體剖分計(jì)算域. 為處理傾斜的地形,Fox等[9]和Wannamaker等[10]將四邊形單元分成四個(gè)三角單元, 該方法本質(zhì)上仍是基于結(jié)構(gòu)化網(wǎng)格剖分, 對(duì)于復(fù)雜的模型仍難以剖分. 非結(jié)構(gòu)化網(wǎng)格的采用, 使得這一問(wèn)題迎刃而解, 拓展了有限元模擬的應(yīng)用范圍[11?13]. 另外, 直接計(jì)算二次場(chǎng)以提高計(jì)算精度, 基于二次場(chǎng)的算法其一次場(chǎng)可解析地求得, 不存在誤差, 雖然二次場(chǎng)由數(shù)值法解得, 存在誤差, 但二次場(chǎng)相對(duì)于一次場(chǎng)較小, 誤差比重不大,從而整體上提高了總場(chǎng)的精度.

可從兩個(gè)角度來(lái)求解三維電磁場(chǎng)問(wèn)題: 其一,基于勢(shì)(場(chǎng)的矢勢(shì)和標(biāo)勢(shì))求解, 然后對(duì)勢(shì)微分求場(chǎng), 微分過(guò)程會(huì)影響解的精度[14]; 其二, 直接基于場(chǎng)求解. 本文采用后者, 從二次電場(chǎng)的偏微分方程來(lái)求解. 節(jié)點(diǎn)有限元法在處理三維矢量電磁場(chǎng)時(shí),存在一些問(wèn)題: 第一, 不滿足電性分界面上法向電場(chǎng)不連續(xù); 第二, 不滿足無(wú)源區(qū)單元內(nèi)電流密度散度為零. 上述問(wèn)題直接違反麥克斯韋方程組, 造成偽解. Jin[15]和柳建新等[16]引入散度校正來(lái)壓制偽解, 但不能完全消除. 本文采用矢量有限元法, 使用Nédélec型矢量形函數(shù), 很好地克服了節(jié)點(diǎn)有限元法的弊端. 另外, 在算法設(shè)計(jì)中, 考慮了磁導(dǎo)率參數(shù)的變化, 可以模擬磁導(dǎo)率不均勻地對(duì)地球電磁擴(kuò)散產(chǎn)生的影響.

作為研究地球內(nèi)部電性構(gòu)造的一種重要的地球物理手段, 大地電磁(MT) 法廣泛應(yīng)用于多個(gè)領(lǐng)域. 深部低頻段, 地殼和上地幔的MT電阻率結(jié)構(gòu)成像, 對(duì)巖石圈深部和地球動(dòng)力學(xué)的演化過(guò)程的研究起著重要作用[17,18]; 淺部高頻段, 廣泛應(yīng)用于金屬礦產(chǎn)資源勘探、地?zé)峥碧?、工程調(diào)查和地下水檢測(cè)等方面[19?21]. 本文首先在前人的基礎(chǔ)上, 采用更為靈活的非結(jié)構(gòu)化網(wǎng)格剖分方式, 導(dǎo)出了三維MT基于二次電場(chǎng)的矢量有限元公式及阻抗求取公式, 對(duì)復(fù)雜不規(guī)則模型進(jìn)行了模擬, 取得了理想的結(jié)果; 然后, 討論了三維地形對(duì)模擬結(jié)果的影響,通過(guò)與二維的模型正演結(jié)果的對(duì)比, 表明了對(duì)三維數(shù)據(jù)進(jìn)行二維處理的不合理; 最后, 詳盡分析了不均勻磁導(dǎo)率對(duì)模擬結(jié)果的影響.

2 正演理論分析

2.1 方程推導(dǎo)

式 中E和H為 分 別 電 場(chǎng) 和 磁 場(chǎng) ;σ為 電 導(dǎo) 率 ;μ=μrμ0為磁導(dǎo)率,μ0為真空中磁導(dǎo)率,μr為相對(duì)磁導(dǎo)率;ε為真空中介電常數(shù);ω為角頻率. 由方程 (1)可導(dǎo)出

將總場(chǎng)分成一次場(chǎng)和二次場(chǎng), 即 E =Ep+Es , 其中下標(biāo)p和s分別表示一次場(chǎng)(primary field)和二次場(chǎng)(secondary field). 在背景模型中(見(jiàn)圖1),背景電導(dǎo)率和磁導(dǎo)率分別為 σ0 和 μ0 :

由 (2)和(3)式聯(lián)合, 考慮磁導(dǎo)率的差異, 可得

圖1 背景模型示意圖Fig.1. Schematic diagram for background model.

2.2 邊界條件

得到微分方程后, 還要有邊界條件才能組成邊值問(wèn)題. 當(dāng)背景模型和實(shí)際一維模型不同時(shí), 比如實(shí)際模型為三層模型(含異常體), 將背景模型設(shè)為兩層模型, 假設(shè)外邊界距離異常體足夠遠(yuǎn), 異常體對(duì)邊界的影響可忽略, 邊界上二次場(chǎng)為實(shí)際一維模型和背景模型一次場(chǎng)的差:

式中下標(biāo)1和2分別表示背景模型和實(shí)際一維模型. 當(dāng)背景模型和實(shí)際一維模型相同時(shí), 外邊界上僅有一次場(chǎng), 二次場(chǎng)為零, (5)式變成為齊次第一類邊界條件; 當(dāng)?shù)叵履Ｐ蜑槎S模型, 其中包含一個(gè)三維異常體時(shí), 仍將一維模型設(shè)為背景模型, 用二維有限元法計(jì)算出邊界上二維模型的數(shù)值解, 代入(5)式中的本文利用走向方向延伸的長(zhǎng)度遠(yuǎn)大于剖面(與走向垂直)上的尺寸的三維模型來(lái)近似二維模型, 故背景模型均為一維層狀介質(zhì)模型.

2.3 一次場(chǎng)的計(jì)算

背景模型(見(jiàn)圖1), 其一次場(chǎng)可解析地求出[16].為避免由電場(chǎng)和磁場(chǎng)中指數(shù)增加項(xiàng)引起的計(jì)算機(jī)數(shù)據(jù)類型的越界和計(jì)算機(jī)對(duì)較大數(shù)據(jù)的截?cái)嗾`差,假設(shè)源為x方向, 將一次電場(chǎng)和磁場(chǎng)的形式解設(shè)為

令

空氣層中, 在模型頂面z=z1上:

可解得

由(7)式邊界條件, 可得

將a1和b1代入(9)式依次求出求各層的系數(shù)ai和bi, 將ai和bi代入 (6) 式, 可計(jì)算模型中任意點(diǎn)的一次場(chǎng).

2.4 加權(quán)余量法

通過(guò)加權(quán)余量法將微分方程過(guò)渡到有限元方程, 進(jìn)而采取有限元法來(lái)求解. 本文使用加權(quán)余量法,取權(quán)重為 δ Es , dv表示單元體積分, 由 (4)式可得

利用矢量運(yùn)算公式?·(A×B)=?×A·B-?×B·A, 邊界上滿足第一類邊界條件, 且單元內(nèi) 物性參數(shù)為常數(shù), (10)式可化為

2.5 矢量有限元法求解

節(jié)點(diǎn)有限元法能有效地模擬直流問(wèn)題, 但處理交變的矢量電磁場(chǎng)時(shí), 存在一些問(wèn)題: 第一, 不滿足無(wú)源的區(qū)域單元內(nèi)電流密度的散度為零; 第二,在電性分界面上, 不滿足法向電流密度連續(xù), 即法向電場(chǎng)不連續(xù). 這些問(wèn)題明顯違反麥克斯韋方程組, 是節(jié)點(diǎn)有限元法理論上的缺陷. 雖然引入一些校正手段, 如散度校正, 可以在一定程度上壓制偽解, 但并不能從根本上消除. 矢量有限元法由于將場(chǎng)賦到棱邊上, 自然滿足切向場(chǎng)的連續(xù), 對(duì)法向場(chǎng)沒(méi)有強(qiáng)制性的約束, 允許法向場(chǎng)的不連續(xù), 為法向電流密度連續(xù)提供了條件. 本文采用非結(jié)構(gòu)四面體單元網(wǎng)格進(jìn)行空間離散(Gmsh開(kāi)源軟件產(chǎn)生), 選用Nédélec型矢量形函數(shù)[15]:

式中 li 為第 i條邊的長(zhǎng)度; L i1 和 L i2 分別為第 i 條邊的 i1 和 i2 節(jié)點(diǎn)的局部坐標(biāo)(體積坐標(biāo)).

由(13)式可知, 矢量有限元法自然滿足無(wú)源區(qū)單元內(nèi)電流密度的無(wú)散.

三維四面體單元中棱邊的編號(hào)規(guī)則如圖2所示, 單元內(nèi)任一點(diǎn)的場(chǎng)可由形函數(shù)(12)式和插值公式(14)得到:

圖2 四面體單元及其棱邊編號(hào)示意圖Fig.2. Schematic diagram for tetrahedral element and edge numbers.

將(14)式代入有限元公式(11)可得

式 中Es和Ep為 6×1 的 標(biāo) 量 列 向 量 ;K1 和K2 為6×6的矩陣,

方程(15)的右端項(xiàng)為源項(xiàng), 其中Ep為已知的一次場(chǎng)棱邊值. 在計(jì)算一次場(chǎng)時(shí), 只能計(jì)算具體某節(jié)點(diǎn)的值, 而不能直接得到棱邊值, 需要將節(jié)點(diǎn)的一次場(chǎng)轉(zhuǎn)化到棱邊上. 對(duì)于四面體單元, 已知四個(gè)節(jié)點(diǎn)的一次場(chǎng), 如何求六個(gè)棱邊的值, 有兩種處理方法: 方法一, 解方程法, 四個(gè)節(jié)點(diǎn), 每個(gè)節(jié)點(diǎn)一次場(chǎng)分為x,y和z三個(gè)分量, 即總共 12 個(gè)方程, 求6個(gè)未知量(棱邊值), 顯然是一個(gè)矛盾方程, 可在方程的兩邊同時(shí)左乘系數(shù)矩陣的轉(zhuǎn)置, 得到一個(gè)近似解; 方法二, 投影法, 棱邊中點(diǎn)的一次場(chǎng)向棱邊方向的投影值近似為棱邊值. 節(jié)點(diǎn)往棱邊上轉(zhuǎn)化時(shí), 不是嚴(yán)格的一一對(duì)應(yīng), 上述兩種方法都要近似,本文采用易于操作實(shí)現(xiàn)的方法二.

將局部單元?jiǎng)偠染仃嚁U(kuò)展成總體剛度矩陣, 其

有限元方程:

邊界條件(5)式的加載, 對(duì)方程(16)式系數(shù)矩陣和右端項(xiàng)做部分修改[22]. (16)式方程條件數(shù)較大, 采用迭代法求解收斂慢[23], 特別是低頻部分.本文利用PARDISO求解器(MKL庫(kù)函數(shù))進(jìn)行求解, 系數(shù)矩陣采用CSR格式壓縮存儲(chǔ), 算法速度快, 穩(wěn)定性強(qiáng).

3 MT 響應(yīng)的計(jì)算

由(16)式得到所有棱邊值后, 通過(guò)(14)式的插值公式可計(jì)算計(jì)算域內(nèi)任一點(diǎn)二次電場(chǎng), 將二次電場(chǎng)代入(1)式麥克斯韋方程, 可計(jì)算出二次磁場(chǎng)(輔助場(chǎng)). 將二次場(chǎng)和解析的一次場(chǎng)相加得到總場(chǎng),通過(guò)(17)式計(jì)算阻抗Z:

目前國(guó)內(nèi)學(xué)者仍多沿用處理二維的方式處理三維問(wèn)題, 默認(rèn)Zxx=0 和Zyy=0 , 則Zxy=Ex/Hy和Zyx=Ey/Hx. 在三維情況下Zxx和Zyy都存在,不等于零的, 上述方法必然會(huì)引入誤差. 本文采用兩個(gè)正交的源, 即x和y方向, 對(duì)模型兩次模擬得到 兩組數(shù)據(jù), 由(18)式計(jì)算阻抗張量:

式中電場(chǎng)和磁場(chǎng)的下標(biāo)1和2分別表示x和y方向不同的源模擬所得數(shù)據(jù). 得到阻抗張量后, 可通過(guò)(19)式計(jì)算視電阻率和相位:

4 算法驗(yàn)證

20世紀(jì)80年代以來(lái), 隨著計(jì)算機(jī)計(jì)算能力的提升, 三維電磁場(chǎng)數(shù)值模擬方法隨之快速發(fā)展. 對(duì)于MT問(wèn)題, 一維層狀介質(zhì)和個(gè)別簡(jiǎn)單二維模型具有解析解, 三維模型則不存在解析解. 對(duì)于三維模型, 當(dāng)不同的數(shù)值模擬方法得到不同的結(jié)果時(shí),孰優(yōu)孰劣, 沒(méi)有一個(gè)客觀的比對(duì)標(biāo)準(zhǔn). COMMEMI(Comparison of Modelling Methods for Electromagnetic Induction)是一個(gè)眾多學(xué)者參與的國(guó)際合作項(xiàng)目[23], 不同的學(xué)者采用不同數(shù)值模擬方法計(jì)算相同模型的MT響應(yīng), 然后計(jì)算均值和標(biāo)準(zhǔn)差,為新的算法提供一個(gè)比對(duì)標(biāo)準(zhǔn), 當(dāng)新的算法的模擬結(jié)果落在參考結(jié)果的標(biāo)準(zhǔn)差以內(nèi), 認(rèn)為該算法模擬結(jié)果正確.

4.1 一維層狀模型

一維模型如圖3所示. 第一層電阻率為100 ?·m ,層厚度為2 km; 第二層為基底層延伸到到無(wú)窮遠(yuǎn),電阻率為400 ?·m; 該模型屬于G型地電模型.該模型剖分為31154個(gè)節(jié)點(diǎn)182806個(gè)單元, 單頻的計(jì)算時(shí)間為8.84 s, 測(cè)點(diǎn)位于模型中心, 離邊界距離為20 km. 圖 4為基于總場(chǎng)和二次場(chǎng)的二次插值的模擬結(jié)果與解析解的對(duì)比. 由圖4可看出, 總

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

場(chǎng)和二次場(chǎng)算法的模擬結(jié)果都與解析解吻合, 進(jìn)一步研究發(fā)現(xiàn), 相對(duì)于總場(chǎng)算法, 二次場(chǎng)算法的模擬結(jié)果更為精確, 其中相位曲線的差異較為明顯. 這是由于基于二次場(chǎng)的算法, 其一次場(chǎng)可解析地求得, 不存在誤差, 雖然二次場(chǎng)由數(shù)值法解得存在誤差, 但二次場(chǎng)相對(duì)于一次場(chǎng)較小, 誤差比重不大,從而提高了總場(chǎng)的精度. 因此相比于總場(chǎng)算法, 二次場(chǎng)算法的誤差較小, 結(jié)果也較為理想.

圖4 基于總場(chǎng)和二次場(chǎng)的二次插值的模擬結(jié)果與解析解的對(duì)比T為周期) (a) 視電阻率 ρs; (b) 相位 ?Fig.4. Comparisons of modelling results of quadratic interpolation and analytical solutions based on total fields and secondary fields: (a) Apparent resistivity ρs; (b) phase ? .

圖5 COMMEMI-3D1 模型, 其中 (a)圖中虛線為 x 和 y 測(cè)線Fig.5. COMMEMI-3D1 model. The dashed lines in (a) are the x and y survey lines.

4.2 三維模型COMMEMI-3D1

對(duì)于三維地電模型, 本文采用COMMEMI-3D1 (見(jiàn)圖 5, 其中ρ0,ρ分別代表背景和異常體的電阻率)模型驗(yàn)證算法的正確性, 生成的四面體網(wǎng)格82992個(gè), 由圖6可知(三角符號(hào)為平均值, 豎線為標(biāo)準(zhǔn)差), 對(duì)于 0.1 Hz的剖面, XY 模式和YX 模式,x測(cè)線和y測(cè)線, 模擬的結(jié)果均在參考結(jié)果的標(biāo)準(zhǔn)差以內(nèi), 且較為靠近均值; 對(duì)于10 Hz的剖面(見(jiàn)圖7), YX模式勉強(qiáng)落在標(biāo)準(zhǔn)差以內(nèi), 對(duì)于XY模式,x測(cè)線上?500—500 m和y測(cè)線上?1000—1000 m(即異常所在的區(qū)域), 模擬的結(jié)果明顯比參考結(jié)果偏低. Mitsuhata 和 Uchida[25]基于勢(shì) (電流密度的矢勢(shì)和標(biāo)勢(shì))的算法, 采用結(jié)構(gòu)化六面體網(wǎng)格, 混合使用矢量有限元和節(jié)點(diǎn)有限元法, 模擬得到的XY模式10 Hz剖面的結(jié)果相對(duì)于參考結(jié)果偏低, 為了驗(yàn)證其結(jié)果, Mitsuhata 和 Uchida[25]還用交錯(cuò)網(wǎng)格有限差分算法進(jìn)行了模擬, 結(jié)果也偏低, 本文結(jié)果與其結(jié)論相吻合. 綜合以上討論, 驗(yàn)證了本文算法的正確性.

4.3 三維模型COMMEMI-3D2

為了進(jìn)一步驗(yàn)證本文算法的精度, 本文計(jì)算了COMMEMI-3D2 (見(jiàn)圖8)并將結(jié)果和基于電場(chǎng)標(biāo)量勢(shì)和磁場(chǎng)矢量勢(shì)的T-?算法進(jìn)行對(duì)比[25]. 該模型生成的四面體網(wǎng)格384992個(gè), 計(jì)算頻率0.001 Hz,耗時(shí) 23.86 s, 由圖 9 可知, 對(duì)于x方向測(cè)線, XY模式下本文結(jié)果和T-?算法吻合良好.

圖6 (a), (b) x 和 (c), (d) y 測(cè)線模擬結(jié)果和 Zhdanov 等[24]的結(jié)果對(duì)比 (剖面為 0.1 Hz)Fig.6. Comparisons of modelling results and Zhdanov et al[24] of (a), (b) x survey line and and (c), (d) y survey line (profile is 0.1 Hz).

圖7 (a), (b) x 和 (c), (d) y 測(cè)線模擬結(jié)果和 Zhdanov 等[24]的結(jié)果對(duì)比 (剖面為 10 Hz)Fig.7. Comparisons of modelling results and Zhdanov et al[24] of (a), (b) x survey line and and (c), (d) y survey line (profile is 10 Hz).

圖8 COMMEMI-3D2 模型Fig.8. COMMEMI-3D2 model.

圖9 模擬結(jié)果和 T -? 算法的結(jié)果[25]對(duì)比 (剖面為 0.001 Hz)(a) 視電阻率; (b)相位Fig.9. Comparisons of modelling results and T -? algorithm results[25] (profile is 10 Hz): (a) Apparent resistivity;(b) phase.

5 算 例

5.1 橢球體模型

相比于規(guī)則六面體網(wǎng)格, 非結(jié)構(gòu)化網(wǎng)格的長(zhǎng)處在于模擬不規(guī)則地下結(jié)構(gòu), 現(xiàn)設(shè)計(jì)六面體網(wǎng)格難以實(shí)現(xiàn)的復(fù)雜模型以突顯本文算法的優(yōu)點(diǎn). 橢球體模型如圖10所示, 橢球體模型模擬結(jié)果切片圖如圖 11 所示, 頻率代表似深度. 由圖 11 可看出, 在異常體中心處及其附近的切片, XY和YX模式的視電阻率和相位剖面對(duì)異常體均有明顯的反映, 且隨著切片遠(yuǎn)離異常體, 異常逐漸減小.

圖10 橢球體模型 (a) xy 剖面; (b) xz 剖面Fig.10. Ellipsoidal model: (a) xy profile; (b) xz profile.

圖11 橢球體模型模擬結(jié)果切片圖 (a) XY 模式視電阻率; (b) YX 模式視電阻率; (c) XY 模式相位; (d) YX 模式相位Fig.11. Slices of modelling results for ellipse model: (a) XY-mode apparent resistivities; (b) YX-mode apparent resistivities;(c) XY-mode phase; (d) YX-mode phase.

5.2 起伏地形模型

討論三維地形模型之前, 先計(jì)算二維純地形模型(見(jiàn)圖12), 通過(guò)對(duì)比三維和二維算法的模擬結(jié)果, 來(lái)驗(yàn)證本文算法對(duì)帶地形模型的正確性. 三維數(shù)值模擬時(shí), 用一個(gè)xoz剖面上同圖12和y方向(走向方向)延伸20 km的三維模型來(lái)近似二維模型, 二維數(shù)值模擬采用有限元法. 由圖13可知, 三維和二維算法模擬結(jié)果吻合得很好, 且二維地形所引起的異常特征和Franke等[14]的結(jié)論一致, 從而證實(shí)了本文算法對(duì)帶地形模型的正確性.

為對(duì)比三維地形對(duì)MT響應(yīng)的影響, 仍沿用圖10模型, 僅將其水平地表改成起伏地表(見(jiàn)圖14),山峰底部長(zhǎng)寬均為 3 km, 頂部長(zhǎng)寬均為 1 km, 高

圖12 二維純地形模型Fig.12. Two-dimensional pure topographical model.

圖13 三維和二維模擬結(jié)果對(duì)比Fig.13. Comparisions of three-dimensional and two-dimensional modelling results.

圖14 山峰模型及其非結(jié)構(gòu)化網(wǎng)格剖分(左圖點(diǎn)為測(cè)線)Fig.14. Peak model and discretization using unstructured grids (the points in left diagram are survey line).

為 0.5 km, 以原點(diǎn)為中心對(duì)稱分布, 類同 Ren 等[26]的帶地形模型, 測(cè)線沿著x軸分布, 生成有限元網(wǎng)格20355個(gè). 對(duì)于圖14帶地形地電模型, 其背景模型為兩層模型, 水平地表為空氣和地下介質(zhì)的分界面, 山峰看作空氣層中的異常體來(lái)處理. 由圖14可看出, 在異常體邊界及地表附近場(chǎng)變化劇烈的區(qū)域?qū)W(wǎng)格進(jìn)行了加密, 另外將黃色計(jì)算域進(jìn)行了適當(dāng)?shù)南蛲庋由? 使邊界條件更為貼近實(shí)際.

在實(shí)際生產(chǎn)中, 考慮到數(shù)據(jù)三維處理的復(fù)雜性和效率, 往往將三維數(shù)據(jù)逐個(gè)剖面二維處理, 然后拼接在一起. 現(xiàn)將三維地形模型簡(jiǎn)化成測(cè)線所在剖面 (y= 0 剖面)二維地形模型 (見(jiàn)圖 12), 對(duì)比二維算法與三維算法正演模擬的結(jié)果差異, 為全面討論, 另加一個(gè)同圖14一樣尺寸的山谷模型. 為突顯地形影響, 背景電阻率和橢球電阻率均設(shè)為100 ?·m , 僅含地形因素導(dǎo)致視電阻率異常. 由圖 15(a)和圖 15(b)可知, 對(duì)于山峰模型, 三維的XY模式和二維的TM模式的模擬結(jié)果趨勢(shì)大致相同, 但具體量值有區(qū)別; 三維的YX模式和二維的TE模式幾乎完全不同, 隨著頻率的降低差別越來(lái)越大, 電阻率異常完全相反, 剖面為 0.1 Hz時(shí),二維顯示輕微的正異常, 三維則顯示為負(fù)異常, 這是由于二維模型認(rèn)為走向方向上無(wú)限延伸, 不存在地形起伏, 而三維模型包含地形在走向上的變化,深部低頻段電阻率負(fù)異常恰恰是由于走向方向的地形起伏造成的. 由圖13(c)和圖13(d)可看出,對(duì)于山谷模型, 結(jié)論和山峰模型一致, 但三維的XY模式和二維的TM模式差異相比山峰模型更大. 相比于二維地形, TE模式高頻段和TM模式全頻段受地形影響, 顯然, 三維地形對(duì)模擬結(jié)果的影響更為嚴(yán)重和復(fù)雜.

圖15 三維和二維模擬結(jié)果對(duì)比 (a), (b)山峰模型; (c), (d)山谷模型Fig.15. Comparisions of three-dimensional and two-dimensional modelling results: (a), (b) Peak model; (c), (d) valley model.

5.3 磁導(dǎo)率異常模型

以往的三維數(shù)值模擬將磁導(dǎo)率假設(shè)為真空中的磁導(dǎo)率, 且認(rèn)為在整個(gè)計(jì)算域均勻分布, 在磁導(dǎo)率異常較大的礦區(qū)該假設(shè)明顯不合理. Mukherjee和Everett[27]和Ren等[26]指出磁導(dǎo)率對(duì)數(shù)值模擬的結(jié)果有重要的影響. 模型仍沿用圖10中的橢球體模型, 將背景電阻率設(shè)為 100 ?·m, 背景磁導(dǎo)率設(shè)為真空磁導(dǎo)率, 測(cè)點(diǎn)分布在x軸上, 該模型生成的四面體單元個(gè)數(shù)為15499個(gè). 對(duì)于圖16(a)和圖16(b), 為突出磁導(dǎo)率因素的影響, 將橢球體電阻率也設(shè)為100 ?·m, 該模型地層即為電性均勻半空間, 僅含磁導(dǎo)率不均勻, 模擬結(jié)果表現(xiàn)為視電阻率正異常, 且隨著相對(duì)磁導(dǎo)率的增加, XY和YX模式的視電阻率異常也隨著增加; 對(duì)于圖16(c)和圖16(d), 討論電性參數(shù)和磁性參數(shù)均存在異常的模型響應(yīng), 將橢球體電阻率設(shè)為 0.5 ?·m, 橢球體電阻率引起的視電阻率負(fù)異常被磁導(dǎo)率引起的正異常抵消, 且負(fù)異常隨著相對(duì)磁導(dǎo)率的增加逐漸減小. 觀察有限元方程(15)的右端源項(xiàng), 可分為兩部分: 一部分為電阻率異常引起, 另一部分為磁導(dǎo)率異常引起, 且電阻率異常引起的源項(xiàng)隨頻率的降低而減小, 而磁導(dǎo)率異常引起的源項(xiàng)不隨頻率變化. 因此隨著頻率的降低, 磁導(dǎo)率異常對(duì)模擬結(jié)果影響越來(lái)越大, 在磁導(dǎo)率異常較大的區(qū)域有可能占主導(dǎo)地位. 以上討論表明, 在磁導(dǎo)率存在異常的區(qū)域, 忽略磁導(dǎo)率因素, 僅考慮電阻率參數(shù)的MT數(shù)值模擬會(huì)引入很大的誤差, 磁導(dǎo)率必須作為一個(gè)獨(dú)立的參數(shù)來(lái)對(duì)待.

圖16 不同磁導(dǎo)率模擬結(jié)果對(duì)比 (a), (b) 無(wú)電性異常模型; (c), (d)含電性異常模型Fig.16. Comparisons of modelling results for different magnetic permeability: (a), (b) Model without electrical anomaly;(c), (d) model with electrical anomaly.

6 結(jié) 論

本文采用非結(jié)構(gòu)化網(wǎng)格來(lái)離散計(jì)算域, 能更靈活精確地?cái)M合起伏地形和地下不規(guī)則異常體, 網(wǎng)格的大小可根據(jù)模型設(shè)計(jì)的需求調(diào)整, 如高頻段加密網(wǎng)格、低頻段增大網(wǎng)格和電性分界面加密網(wǎng)格等,提高精度的同時(shí)盡可能地減小計(jì)算量. 同時(shí), 將解析的一次場(chǎng)從總場(chǎng)中分離, 直接計(jì)算二次場(chǎng), 使得數(shù)值誤差僅局限于相對(duì)一次場(chǎng)較小的二次場(chǎng), 從而提高了總場(chǎng)的精度. 矢量有限元法滿足電性分界面上法向電場(chǎng)不連續(xù)和無(wú)源區(qū)單元內(nèi)電流密度散度為零, 克服了節(jié)點(diǎn)有限元法的弊端. 另外, 在推導(dǎo)基于二次電場(chǎng)的雙旋度矢量方程過(guò)程中考慮了介質(zhì)電導(dǎo)率和磁導(dǎo)率的不均勻性, 可以模擬磁導(dǎo)率不均勻模型的電磁場(chǎng). 通過(guò)和COMMEMI模型的已發(fā)表結(jié)果對(duì)比, 證明了本文算法的正確性. 另外,對(duì)不規(guī)則異常體和任意地形的復(fù)雜模型的模擬, 均取得了理想的結(jié)果. 三維地形影響比二維更為嚴(yán)重和復(fù)雜, 用二維算法處理三維模型或數(shù)據(jù)會(huì)引入很大的誤差. 在磁導(dǎo)率存在異常的區(qū)域, 磁導(dǎo)率對(duì)數(shù)值模擬的結(jié)果有著重要的影響, 磁導(dǎo)率必須作為一個(gè)獨(dú)立的參數(shù)來(lái)對(duì)待.