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

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

關鍵詞： 大地電磁法，數值模擬，Chebyshev 多項式，并行計算，直接求解器

中圖分類號：P631 文獻標識碼：A DOI：10. 13810/j. cnki. issn. 1000‐7210. 2024. 06. 020

0 引言

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

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

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

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