基于譜元法的廣域電磁法三維正演模擬

徐錦通， 湯井田,2,3,4*

1 中南大學地球科學與信息物理學院， 長沙 410083 2 中南大學有色金屬成礦預測與地質環境監測教育部重點實驗室， 長沙 410083 3 有色資源與地質災害探查湖南省重點實驗室， 長沙 410083 4 自然資源部覆蓋區深部資源勘查工程技術創新中心， 合肥 230001

0 引言

可控源音頻大地電磁法(CSAMT)是一種重要的地球物理勘探方法.具有觀測范圍廣、工作效率高、測量精度高、適應性強等特點，被廣泛應用于礦產資源勘查、地下水及工程勘查等領域.但也存在一定的缺點，可控源音頻大地電磁法沿用了大地電磁法的觀測方式以及卡尼亞視電阻率計算公式，為滿足遠區測量要求，通常需要在距離發射源數千米至數十千米的區域內接收正交的水平電磁場信號來獲得卡尼亞視電阻率和相位信息，背離了采用人工源使信號增強的初衷，并且在野外實際施工中，由于測區場地條件的限制，有時難以滿足遠區測量的要求.采用的卡尼亞視電阻率計算公式是在假設遠區條件下，將電磁場精確表達式中的高次項近似舍去，然后基于阻抗定義推導得到的， 引入了一定的人為誤差(何繼善，2010).

廣域電磁法是一類人工源頻率域電磁勘探方法(湯井田和何繼善，1994，2005；何繼善，2010)，廣域電磁法繼承了CSAMT 使用人工場源克服場源隨機性的優點，采用較短的收發距進行測量，摒棄了CSAMT遠區信號微弱的劣勢，用適合于全域的不進行簡化的公式進行視電阻率計算，擴展了觀測適用的范圍.同時，廣域電磁法發射源一次發送的是含多個主頻成分偽隨機電流信號，而非變頻發送，并且在測量時只需要觀測電磁場的一個分量，而不是觀測一對正交的電磁場分量，大大提高了觀測速度、精度和野外工作效率.

廣域電磁法解釋目前主要采用的是一維和二維反演技術.然而在地形和地質條件復雜地區，地下目標體實際是典型的三維異常體.此時，若仍然采用一維和二維反演，很難得到可靠的的反演結果，必須發展三維電磁反演技術.三維反演的一個重要難點是如何實現快速高精度的三維正演，目前主流的三維電磁數值模擬方法主要包括積分方程法(Xiong et al.,1992; Kuvshinov et al.2002; Zhdanov et al.,2006)、有限體積法(Jahandari and Farquharson, 2015;Lu and Farquharson, 2020)、有限差分法(Haber et al.,2000; Weiss and Newman,2003；盧杰和李予國, 2019)、有限元法(Ren et al., 2013, 2014; Grayver,2015,秦策等, 2020).廣域電磁法三維正演方面，李帝銓等(2013)采用積分方程法實現了E-Ex廣域電磁法三維正演；彭榮華等(2018)實現了基于二次耦合勢的有限體積法廣域電磁法三維正演計算；周印明等(2021)實現了基于矢量位和標量位的廣域電磁法三維有限元數值模擬.

譜元法是將有限元和譜方法結合發展而來的一種高精度數值模擬方法.同時兼具有限元的靈活性和譜方法的高精度.與傳統高階有限元方法相比，由于采用了高斯正交基函數，插值節點為高斯分布型節點，可避免高階等距插值的Runge現象.Patera(1984)最早提出基于Gauss-Lobatto-Chebyshev 基函數的譜元法，將其應用流體力學數值模擬.R?nquist和Patera(1987)提出了另一種基函數Gauss-Lobatto-Legendre基函數的譜元法.基于GLL基函數的譜元法結合降階的高斯積分可得到完全對角的質量矩陣，基于GLC多項式的譜元法可以解析計算單元矩陣，理論上可獲得更高的精度.在地球物理領域里，譜元法最早被應用于地震勘探中(Seriani et al.,1997;Komatitsch and Tromp, 1999，2002a, 2002b).在電磁勘探領域，Lee等(2006)將譜元法應用于電阻抗層析成像的正演求解器開發，證明了譜元法的高精度及快速收斂特性.Lee等(2006)將基于GLL多項式的混合階譜元法用于求解三維矢量電磁波動方程和三維瞬變電磁問題中，與高階有限元相比計算時間大幅減小.Zhou等(2016)將譜元法與區域分解結合應用于大規模航空瞬變電磁正演計算，有效的節省了計算時間，克服了一定程度的較低頻率下由于計算機精度有限導致的系統矩陣嚴重病態的低頻崩潰問題.Zhou和Shi等(2017)提出一種混合譜元法并且將散度約束條件加入到譜元法中，來克服低頻崩潰問題，在井下電磁勘探正演中取得了理想的效果.Huang 等(2017，2019)實現了基于GLL基函數譜元法頻率域航空電磁法正演及各向異性正演.劉玲等(2018)實現了基于GLC基函數譜元法的海洋電磁正演模擬.陳漢波等(2019)實現了基于混合階矢量基函數的海洋電磁法各項異性正演.Yin等(2019)等實現了基于GLC多項式的航空電磁正演模擬.

譜元法當前主要應用于航空、海洋電磁法，在陸地可控源電磁法中的應用尚未見到，本文在前人工作的基礎上將譜元法引入陸地可控源電磁法的分支——廣域電磁法的正演計算.目前針對廣域電磁法三維正演研究主要是基于位場方程，需要通過差分間接得到電場，會帶來精度上的損失，本文采用電場雙旋度方程作為控制方程直接獲得高精度的電場解，利用一維模型和三維驗證算法的正確性、精度及適用性，設計典型三維地電模型研究廣域電磁法響應特征及其探測能力.

1 基本理論

1.1 廣域電磁法邊值問題

從微分形式的麥克斯韋方程組出發，取時諧因子為e-iω t，可得電性源廣域電磁法滿足的頻率域麥克斯韋方程組

(1)

(2)

(3)

(4)

式中：D、E為電位移矢量和電場強度；B、H為磁感應強度和磁場強度；σ為電導率；ω為角頻率；μ為磁導率；ε為介電常數；Js為外部激發的電性源.

人工源電磁法正演需考慮源的加載及源奇異性的影響，為了避免場源的奇異性，采用二次場算法，將總場分解成背景場和二次場，背景場由解析解計算得到，二次場滿足的公式(1)、(2)的形式為

(5)

(6)

其中

(7)

式中，σp為背景電導率；Ep為一次電場.

將(6)式代入(5)式，消去磁場可得二次電場滿足的矢量赫姆霍茲方程：

(8)

其中，k2=iωμσ+ω2με，為波數.

為保證電場的唯一性，需要考慮邊界條件的加載，本文在正演時對計算區域進行擴邊處理，當擴邊區域足夠大時可以近似認為二次電場已衰減為0，因此本文在外邊界處選取 Dirichlet 邊界條件如下

Es|?Ω=0,

(9)

(8)式和(9)式共同構成廣域電磁法邊值問題.

1.2 邊值問題的譜元求解

采用譜元法求解上述邊值問題，首先將求解區域劃分成若干個六面體單元，在每個單元內二次電場可用關于x，y，z方向的高階矢量插值基函數展開

(10)

利用伽遼金加權殘差法處理電場矢量方程(8)，令單元內的加權殘差積分為0，并應用第一矢量格林定理與Dirichlet邊界條件可得

(11)

式中，Ne為單元總數，Nb為單元棱邊總數.

(11)式可寫成矩陣方程的形式：

(A-k2M)Es=B,

(12)

式中，A為剛度矩陣，M為質量矩陣，Es為待求二次場列向量，B為右端源向量.

譜元法中可采用基于Gauss-Lobatto-Legendre(GLL)多項式的基函數和Gauss-Lobatto-Chebyshev(GLC)多項式基函數，本文選取的是GLL基函數，GLL基函數是定義在標準參考區間[-1，1]上的，在三維參考單元中ξ,η,ζ∈[-1,1]×[-1,1]×[-1,1]中，N階矢量基函數可表示為一維基函數的混合階張量積(Lee et al.,2006)：

(13)

其中一維GLL基函數表達式為

(14)

式中，LN(ξ)、L′N(ξ)為N階Legendre正交多項式及其導數.ξi為GLL插值節點，即方程(1-ξ2)L′N(ξ)=0的零點.

物理單元的插值基函數與參考單元插值基函數的映射關系如下

(15)

(16)

其中J為坐標映射雅克比矩陣

(17)

由此可將(11)式中的物理坐標系下單元矩陣轉換為標準參考坐標系下的表達形式

(18)

(19)

進而通過在參考域中求解線性方程組(12)可以獲得各棱邊上的場值，并通過式(10)插值獲得任意點的場值.方程組的求解本文采用直接解法開源求解器PARDISO進行求解.

1.3 E-Ex廣域視電阻率

廣域電磁法根據場源形式以及觀測方式可以分為多種，考慮到野外實際情況，目前為止采用水平電流源發射信號測量電場的x單分量的E-Ex廣域電磁法應用最為廣泛，本文研究主要針對這類廣域電磁法.

這里以電場水平分量Ex來闡述E-Ex廣域電磁法和廣域視電阻率的概念.均勻大地表面水平電流源的電場x分量的計算公式為

(20)

式中，I為電流;dL為電偶極源的長度；k為均勻半空間的波數；r為收發距，即觀測點距偶極子中心的距離；φ為電偶極源方向和源的中點到接收點矢徑之間的夾角.

根據式(21)可以定義廣域意義上的視電阻率：

(21)

式中，

ΔVMN=Ex·MN,

(22)

(23)

式中，MN為接收電極距.不難看出，公式(21)為自變量為復數的指數復變函數，采用一般的代數方法無法簡單的提取出其中含有的的電阻率信息，廣域電磁法提出的方法是采用計算機迭代的方法進行提取.

2 數值試驗

為驗證本文算法的正確性和精度，設計了典型地電模型進行正演計算與對比分析.計算平臺為Linux Fedora 21小型AMAX服務器，內存250 GB，主頻為2.3 GHz.

2.1 一維層狀模型精度驗證

為驗證本文算法的正確性及精度，設計如圖1所示的三層模型.第一層電阻率為10 Ωm ，厚度為500 m，第二層電阻率為100 Ωm，厚度為200 m，第三層電阻率為500 Ωm.發射源為位于坐標原點的電偶極子，電流為1 A，長度為1 m，接收點位置坐標為(6000 m，0，0)，發射頻率為1～8192 Hz對數等間隔的25個頻點.整個計算區域為100 km×100 km×100 km.

圖1 一維層狀模型示意圖

分別采用3階、4階基函數的譜元法進行正演計算，并與有限單元法及解析解對比.圖2給出了不同階數基函數的Ex分量幅值隨頻率變化曲線及與解析解的相對誤差.從圖中可以看出，有限元法、譜元法Ex響應均與解析解擬合較好，高頻部分誤差相對于低頻要小，精度更高，這是因為高頻電磁波衰減較快，更加符合擴邊區域外邊界二次電場衰減為零的邊界條件假設. 對于譜元法，隨著基函數階數的提升，相對誤差顯著減小，表明譜元法可以通過提高基函數的階數達到提高精度的目的.表1統計了3階譜元法與有限元兩種方法的計算參數，從表中可以看出，當基函數的階數為3階時，兩者的精度相當，但譜元法的自由度相對于有限元更少，相應的計算時間也更少，說明高階譜元法可以在稀疏網格下以更少的自由度獲得更高精度的解，對于網格依賴性較低，從而可以節省計算資源消耗.

圖2 一維層狀模型Ex響應及相對誤差曲線

表1 譜元法與有限元法計算時間統計表

2.2 三維模型對比測試

為進一步驗證本文算法的精度及對三維模型的有效性，設計如圖3所示的三維低阻體模型.背景電阻率為50 Ωm，低阻異常體電阻率為5 Ωm，頂部埋深為100 m，幾何尺寸為120 m×200 m×400 m，發射源為沿x方向的電偶極子源，長度為100 m，中心坐標為(50 m，0 m，0 m)，距離異常體中心1 km，發射電流為1 A，發射頻率為3 Hz.采用背景為50 Ωm均勻半空間模型進行一次場計算.在異常體正上方沿x方向布置一條測線，測線范圍為400～1400 m，點距為20 m.整個模擬區域為80 km×80 km×80 km,最小網格尺寸為50 m×50 m×50 m，計算時間為153 s.

圖3 三維低阻體模型示意圖

圖4為本文算法結果與積分方程法(湯井田等，2018)，有限元法(Ansari and Farquharson，2014)等計算結果Ex分量實部與虛部對比曲線.從圖中結果可以看出，本文的結果與前人積分方程、有限元法結果均吻合較好，驗證了本文算法對三維模型良好的適用性及精度，可以用來進行廣域電磁法三維正演分析.

圖4 低阻體模型Ex響應曲線

2.3 廣域電磁法響應特征分析

2.3.1 低阻異常體

為分析廣域電磁法響應特征及探測能力，設計如圖5所示的三維低阻異常體模型，電阻率為1 Ωm的低阻體賦存在電阻率為100 Ωm的背景介質中，異常體大小為500 m×500 m×500 m，頂部埋深為300 m.發射采用沿x方向的電偶極子源，長度為1 m，發射電流為1 A，頻率為0.5～8192 Hz對數等間隔的15個頻點.源的坐標為(0 m，-5000 m，0 m)，取異常體中心在地面投影為坐標原點.在異常體正上方沿x坐標軸布置一條測線，長度為2 km，點距為50 m.整個模擬區域為80 km×80 km×80 km,最小網格尺寸為100 m×100 m×50 m，計算時間為127 s.

圖5 三維低阻體模型示意圖

利用本文算法對該模型進行正演，并同時計算卡尼亞視電阻率和廣域視電阻率.圖6為該模型卡尼亞視電阻率與廣域視電阻率頻率切片圖.從圖中可以看出無論是卡尼亞視電阻率還是廣域視電阻率均能較好的反映出低阻異常體的存在，表現出明顯的低阻響應.但在頻率低于10 Hz時，卡尼亞視電阻率整體都出現了明顯的升高，大多數觀測值都高于背景電阻率，這是因為隨著頻率的降低，電磁波在地下介質傳播的趨膚深度變大，此時測線所在位置已經不滿足遠區(五倍以上趨膚深度)的條件，出現了由場源引起的非平面波效應，卡尼亞視電阻率公式不再適用，產生了畸變.而廣域視電阻率在低頻段仍能較好的反映出異常體存在及地下電阻率信息，不受遠區測量的限制，表明廣域電磁法相對于傳統可控源電磁法具有更大的觀測范圍.

圖6 低阻體模型卡尼亞視電阻率與廣域視電阻率切片圖

在異常體上方布置一2 km×2 km的測網，線距為100 m，點距為100 m.發射頻率為8 Hz，其他源參數保持不變.圖7為卡尼亞視電阻率與廣域視電阻率測區平面圖.源的位置為(0 m，-5000 m,0 m),圖中從上到下觀測位置離源的距離依次減小.圖中可以看出，卡尼亞視電阻率與廣域視電阻率均對低阻異常體有響應，隨著測線與源的距離的減小，卡尼亞視電阻率出現了急劇的增大，無法真實客觀的反映出地下電阻率信息.而廣域視電阻率在整個測區內均能較好的反應出異常體的位置及地下電阻率信息.表明了廣域電磁法在較小收發距下也能獲得地下真實地電信息，即在同等收發距下相對于CSAMT 具有更深的探測能力.

圖7 低阻體模型卡尼亞視電阻率與廣域視電阻率平面分布圖(8 Hz)

2.3.2 組合異常體

為進一步探究廣域電磁法對復雜異常體的分辨能力，設計如圖8 所示的組合異常體模型.均勻半空間背景電阻率為100 Ωm，低阻異常體電阻率為1 Ωm，高阻異常體電阻率為1000 Ωm，頂部埋深均為200 m，大小均為400 m×400 m×400 m.在異常體上方布置一2 km×2 km的測網，測線沿x方向布置，線距100 m，點距50 m.發射源的位置為(0 m，-5000 m，0 m)，采用x方向電偶源，長度為1 m，發射電流為1 A，發射頻率為8 Hz、32 Hz.整個模擬區域為80 km×80 km×80 km,最小網格尺寸為100 m×100 m×50 m，計算時間為162 s.

圖8 組合異常體模型示意圖

圖9為電場Ex分量和磁場Hy分量平面等值線圖.由圖可知，電場、磁場均表現出隨收發距增大而衰減的規律.電場較好的反映出了地下異常體的存在，低阻體上方呈現低電場值的特征，等值線向源一側彎曲，高阻體上方呈現高電場值特征，等值線向背離源一側彎曲.而磁場未出現明顯的高低阻異常特征.圖10分別為8 Hz、32 Hz卡尼亞視電阻率與廣域視電阻率平面分布圖.從圖中可以看出，在32 Hz時，所有測點基本滿足遠區條件，卡尼亞視電阻率與廣域視電阻率表現出類似的特征，未出現明顯的畸變.頻率為8 Hz時，部分測點已不滿足遠區條件，卡尼亞視電阻率受非平面波效應的影響，在靠近源一側的部分測點視電阻率值出現明顯的增大，與地下真實電阻率相差較大，并且對高阻異常體響應有一定的覆蓋.而廣域視電阻率仍能較好的反映出兩個異常體的幾何邊界位置和地下電阻率信息，體現了廣域電磁法對復雜地下異常體良好的分辨能力.此外，我們可以看出，無論是卡尼亞視電阻率還是廣域視電阻率對于低阻體的響應相對于高阻體都要更強.這是因為低阻異常體吸引電流，使得異常體上方電場等值線變密，而高阻體排斥電流，使得異常體上方電場等值線變疏，低阻異常體電場響應更強.而磁場對高阻體，低阻體均不敏感，因此無論是基于電場和磁場的比值提取到的卡尼亞視電阻率還是僅用電場值提取到的廣域視電阻率對于低阻體都更加敏感.

圖9 組合異常體模型Ex與Hy平面等值線圖

圖10 組合異常體卡尼亞視電阻率與廣域視電阻率平面分布

3 結論

(1)本文將基于GLL基函數的譜元法引入廣域電磁法三維正演，為避免場源奇異性，采用二次場算法計算.通過與解析解、有限元法對比，驗證了算法的高精度與高效率，設計三維模型驗證了其良好的三維適用性和有效性.

(2)設計典型三維地電模型進行廣域電磁法正演，分析了廣域電磁法響應特征及探測能力，結果表明廣域電磁法可以克服傳統可控源的電磁法的非平面波效應的限制，廣域視電阻率在非遠區仍能較好的反映地下真實的地電信息，拓展了人工源電磁法的觀測范圍，對于提高探測深度及野外施工效率具有十分重要的意義.

(3)本文研發的廣域電磁法譜元法正演求解器具有高精度計算復雜電磁響應的能力，能夠為廣域電磁數據三維反演提供核心支撐，預期能夠大幅度提高廣域電磁數據三維反演水平.

致謝感謝國家自然科學基金重點項目對本文研究的資金支持.感謝兩位匿名審稿專家及編輯對本文提出的寶貴意見和建議.