基于GLC 多項(xiàng)式譜元法的大地電磁場二維正演模擬

摘要： 為了提高大地電磁場數(shù)值模擬的精度和效率，提出基于Gauss‐Lobatto‐Chebyshev（GLC）基函數(shù)譜元法的大地電磁場正演模擬方法。該方法首先從麥克斯韋方程組出發(fā)，推導(dǎo)了二維大地電磁場邊值問題；然后基于Galerkin 加權(quán)余量法將微分形式的邊值問題轉(zhuǎn)換成積分弱形式； 最后采用GLC 正交多項(xiàng)式插值基函數(shù)對全局問題進(jìn)行離散化，利用Pardiso 求解器求解大型稀疏線性方程組得到大地電磁場，實(shí)現(xiàn)了二維大地電磁場數(shù)值模擬。為了提高數(shù)值模擬效率，文中算法采用變密度規(guī)則網(wǎng)格剖分技術(shù)，即在電性復(fù)雜區(qū)域使用細(xì)網(wǎng)格，在電性均勻區(qū)域使用粗網(wǎng)格，采用OpenMP 編程模式實(shí)現(xiàn)了多個(gè)頻點(diǎn)的并行計(jì)算，以此達(dá)到縮短計(jì)算時(shí)間的目的。通過一維層狀介質(zhì)模型數(shù)值模擬結(jié)果驗(yàn)證了文中算法的正確性和精度，相比于Gauss‐Lobatto‐Legendre（GLL）多項(xiàng)式譜元法的數(shù)值模擬結(jié)果，GLC 多項(xiàng)式譜元法的模擬結(jié)果精度更高。針對國際標(biāo)準(zhǔn)模型COMMEMI 2D‐1 和帶地形模型，進(jìn)行了基于GLC 多項(xiàng)式譜元法、有限差分法和基于三角網(wǎng)格有限元法的正演模擬，三者數(shù)值結(jié)果對比表明，GLC 多項(xiàng)式譜元法計(jì)算精度更高、網(wǎng)格依賴性更低。

關(guān)鍵詞： 大地電磁法，數(shù)值模擬，Chebyshev 多項(xiàng)式，并行計(jì)算，直接求解器

中圖分類號：P631 文獻(xiàn)標(biāo)識碼：A DOI：10. 13810/j. cnki. issn. 1000‐7210. 2024. 06. 020

0 引言

大地電磁測深是蘇聯(lián)學(xué)者Tikhonov[1]和法國學(xué)者Cagniard[2]于20 世紀(jì)50 年代提出來的。電磁法利用天然電磁場作為場源，探測地下電性的垂直變化。該方法在構(gòu)造單元劃分、基底起伏形態(tài)研究、沉積盆地電性分層、油氣田勘探、地?zé)嵴{(diào)查和地震預(yù)測等方面應(yīng)用廣泛。大地電磁場正反演方法是地球物理學(xué)研究的重要領(lǐng)域。正演是反演的基礎(chǔ)，研究高效、高精度的正演算法至關(guān)重要[3]。目前電磁場數(shù)值模擬方法主要包括：積分方程法（IE）[4]、有限差分法（FDM）[5‐6]、有限體積法（FVM）[7‐8] 和有限元法（FEM）[9‐11]。IE 只需對異常區(qū)域進(jìn)行離散，計(jì)算速度快[12]，但它并不適合復(fù)雜異常體的數(shù)值計(jì)算。FDM 具有實(shí)現(xiàn)簡單、計(jì)算速度快等優(yōu)點(diǎn)，然而，由于必須將計(jì)算域剖分為規(guī)則單元，因此在處理復(fù)雜模型時(shí)，其適用性受到限制[13]。FEM 可以采用規(guī)則網(wǎng)格和非規(guī)則網(wǎng)格對研究區(qū)域進(jìn)行剖分，模擬復(fù)雜地形和地質(zhì)體的能力較強(qiáng)，目前該方法已經(jīng)成為一種高效的電磁場正演方法。傳統(tǒng)的FEM 采用一階線性插值，需要精細(xì)離散網(wǎng)格以滿足精度要求，因此需要大量的計(jì)算資源。近年來，很多學(xué)者開展了高階FEM 研究[14]，其特點(diǎn)是在單元內(nèi)采用高階插值基函數(shù)擬合單元內(nèi)的場值變化，該方法精度高但是實(shí)現(xiàn)難度較大，且傳統(tǒng)高階插值會產(chǎn)生Runge 現(xiàn)象[15]，影響數(shù)值模擬結(jié)果的準(zhǔn)確性。

譜元法（Spectral Element Method ，SEM ）是近年來發(fā)展起來的一種基于有限元法和譜方法的高精度數(shù)值模擬方法，它本質(zhì)是一種非線性插值高階有限元法，即P 型有限元法。20 世紀(jì)80 年代，Patera[16]首次提出Gauss‐Lobatto‐Chebyshev（GLC）基函數(shù)譜元法，并將其應(yīng)用到流體力學(xué)數(shù)值模擬。該方法具有精度高和快速收斂的優(yōu)點(diǎn)，許多學(xué)者對其進(jìn)行了大量研究。Maday 等[17]擴(kuò)展了插值基函數(shù)的形式，將Gauss‐ Lobatto‐ Legendre（GLL）多項(xiàng)式作為譜元法的基函數(shù)求解偏微分方程。20 世紀(jì)90 年代，地球物理學(xué)者開始將該方法引入地球物理場數(shù)值模擬。Seriani 等[18]利用GLC 多項(xiàng)式譜元法模擬二維聲波場，取得了很好的結(jié)果； 隨后，Seriani[19]和Maday 等[20]相繼開展了一系列相關(guān)研究，使得譜元法在地震領(lǐng)域得到快速發(fā)展；20 世紀(jì)末，Komatitsch團(tuán)隊(duì)[21‐26]發(fā)表了一系列關(guān)于地震波譜元法正演模擬的文章，并開發(fā)了譜元法開源軟件SPECFEM3DCartesian，在譜元法的基礎(chǔ)上逐步加入混合網(wǎng)格技術(shù)、并行技術(shù)和區(qū)域分解技術(shù)，并取得了很好的應(yīng)用效果。劉有山等[27]、李琳等[28]分別將基于Fekete和Cohen 節(jié)點(diǎn)的三角形網(wǎng)格譜元法應(yīng)用于地震正演模擬，證明了三角形網(wǎng)格譜元法對于復(fù)雜模型也有很好的模擬效果。孟雪莉等[29]為了解決傳統(tǒng)譜元法數(shù)值積分精度不足的問題，提出了一種優(yōu)化積分算法，提高了利用譜元法模擬地震波的數(shù)值精度。21 世紀(jì)初，杜克大學(xué)Liu 等[30]首次將譜元法引入計(jì)算電磁學(xué)，提出了一種Chebyshev 偽譜法，求解矢量波動方程和薛定諤方程。經(jīng)過多年的發(fā)展，譜元法已經(jīng)廣泛應(yīng)用于多種電磁場正演問題的求解。Lee 等[31]利用混合階譜元法求解三維矢量電磁波動方程，相比于有限元法大大縮短了計(jì)算時(shí)間；劉玲等[32]實(shí)現(xiàn)了基于譜元法的頻率域三維海洋可控源電磁場正演，模擬結(jié)果表明，譜元法對網(wǎng)格的依賴程度較之于傳統(tǒng)方法更低； 朱姣等[33]將基于Warp 和Blend 點(diǎn)集的非結(jié)構(gòu)網(wǎng)格譜元法應(yīng)用于任意各向異性介質(zhì)三維直流電阻率正演模擬，解決了傳統(tǒng)譜元法無法模擬帶復(fù)雜地形和異常體的任意各向異性模型電磁響應(yīng)的問題； 方小姣等[34]和周翔等[35]將GLL 多項(xiàng)式譜元法應(yīng)用于大地電磁正演模擬，結(jié)果表明譜元法是大地電磁正演模擬的有效方法。

通常情況下，譜元法選取GLL 多項(xiàng)式為插值基函數(shù)。采用GLL 插值多項(xiàng)式及GLC 積分可以形成對角質(zhì)量矩陣[36]，通常被稱為Mass‐ Lumping 技術(shù)[37]。這種方法顯著減少了系數(shù)矩陣中非零元素的數(shù)量，增強(qiáng)了矩陣的稀疏性，便于求解大型稀疏線性方程。GLL 基函數(shù)的特點(diǎn)在于擁有重合的插值點(diǎn)與積分點(diǎn)，對時(shí)間域問題的迭代求解有顯著優(yōu)勢。但是，為了滿足這種一致性需要，采用降階積分運(yùn)算會帶來一定的精度損失[38]。相較于GLL 多項(xiàng)式，GLC 多項(xiàng)式的優(yōu)勢在于有解析的積分形式，可以保證更高的精度。此外，因?yàn)镃hebyshev 級數(shù)在逼近函數(shù)時(shí)是一致收斂的，理論上GLC 多項(xiàng)式作為基函數(shù)能更好地?cái)M合電磁場變化[39]，可以避免傳統(tǒng)高階插值產(chǎn)生的Runge 現(xiàn)象[15]。

綜上所述，本文將基于GLC 多項(xiàng)式的譜元法應(yīng)用于大地電磁場正演，以此提高大地電磁數(shù)值模擬精度。首先介紹基于GLC 多項(xiàng)式譜元法的大地電磁場正演方法和實(shí)施步驟； 然后利用一維層狀模型的數(shù)值模擬結(jié)果驗(yàn)證本文算法的正確性，并與基于GLL 多項(xiàng)式譜元法的模擬精度進(jìn)行對比； 通過二維模型以及帶地形模型算例，對比分析本文算法與三角網(wǎng)格有限元法以及有限差分法在數(shù)值模擬精度和效率上的差異； 最后，分析并探討并行計(jì)算在提高本文算法效率方面的作用。