多輻射源地空瞬變電磁響應三維數值模擬研究

李貅， 胡偉明， 薛國強

1 長安大學地質工程與測繪學院， 西安 710054 2 中國科學院礦產資源研究重點實驗室， 中國科學院地質與地球物理研究所， 北京 100029 3 中國科學院大學地球與行星科學學院， 北京 100049 4 中國科學院地球科學研究院， 北京 100029 5 西北有色地質礦業集團有限公司， 西安 710054

0 引言

電磁法通過觀測天然或人工源激發的電磁響應獲取地下電性信息(何繼善和薛國強，2018;底青云等，2019).傳統地面電磁法觀測方式在沙漠、戈壁、山地和濕地等地區無法高效地開展工作，航空電磁法可發揮替代作用(圖1).航空電磁法根據發射源布設方式可分為全航空電磁法和地空電磁法，其中全航空電磁法采集效率高，但由于受到載重限制，發射功率有限.地空電磁法將發射系統放置于地面，接收系統布設在飛行平臺上(圖2)，充分地結合了傳統地面電磁法和全航空電磁法的優點，可實現探測深度和工作效率之間的平衡.

圖1 航空電磁探測系統示意圖Fig.1 Schematic diagram of airborne electromagnetic detection system

圖2 多輻射源地空電磁探測系統示意圖Fig.2 Schematic diagram of multi-radiation source semi-airborne electromagnetic detection system

地空瞬變電磁探測系統通常采用單一接地線源向地下發射電磁場信號，并根據地面長偏移瞬變電磁法的模式進行數據處理和解釋(Mogi et al.，1998，2009；嵇艷鞠等，2013；李肅義等，2013).嵇艷鞠等(2013)采用吉林大學研發的地空瞬變電磁系統在江蘇省如東縣開展了海水侵蝕探測，并在內蒙古錫林郭勒盟巴彥寶利格盆地開展了地下水資源探測.Wang等(2013)引入基于多分辨分析的小波分析方法，克服飛行中的基線漂移問題.Ji等(2016)利用小波閾值法在最大限度上剔除數據中的白噪聲，實現非平穩噪聲的識別與剔除.李貅等(2015)利用TEM三維擬地震解釋技術，為地空電磁成像解釋打下基礎.Wu等(2019)提出一種基于小波神經網絡的噪聲處理方法.

國內學者對多輻射場源下地空電磁法的數據處理與反演開展了探索性研究.張瑩瑩等(2015)提出一種基于多源的地空瞬變電磁方法，并給出多分量全域視電阻率定義.張瑩瑩等(2016)推導了多源地空瞬變電磁的視縱向電導和視深度計算公式.Zhou等(2016a,b)解決了多源發射系統的野外裝置布設問題.Zhou等(2018)提出一種張量探測方法，以提升基于一對正交電性源的地空電磁系統對3D目標體的探測精度.本文采用三維矢量有限元法對兩個不同地質體進行多輻射源地空瞬變電磁響應開展模擬研究，分析了多輻射源在不同輻射方向、不同飛行高度的電磁響應分布特征，從而推進多輻射源地空瞬變電磁理論方法和技術體系的建立.

1 瞬變電磁場三維矢量有限元法原理

1.1 邊值問題

時間域麥克斯韋方程組為(李賀，2016)：

(1)

(2)

其中，E為電場強度，H為磁場強度，σ為介質電導率，μ0為介質的磁導率，Js為源電流密度.對(1)式兩端取旋度：

(3)

將(2)式代入(3)式中得：

(4)

整理可得：

(5)

根據電磁場理論，在無源區的兩種導電介質界面，電磁場滿足如下四個邊界條件：

(6)

式中，n為兩種介質界面處的單位法向分量，B為磁感應強度，D為電位移矢量.

對于無窮遠邊界，電場和磁場在無窮遠邊界上的切向分量為零，即：

(7)

1.2 矢量變分方程

根據加權余量法，電場控制方程的加權余量方程為

(8)

其中，f為矢量基函數.根據矢量分析恒等式、高斯公式以及邊界條件，可將(8)式寫為

=0，

(9)

(9)式即為矢量有限元法的矢量變分方程.

1.3 矢量有限元法

圖3 單元中節點與棱邊的關系(引自：李賀，2016)Fig.3 Relationship between nodes and edges in an element (from Li H, 2016)

采用六面體單元進行網格剖分，在單元中棱邊與節點的關系如圖3所示.在六面體單元中，將電場切向分量賦予各個單元的棱邊上.采用Whitney型插值基函數(李賀，2016)，使得插值基函數的散度為0，而旋度不等于0.因此，在計算區域存在介質不連續性時，電場強度的切向分量連續，法向分量不連續，可以有效地避免“偽解”的現象.

在(9)式中，含有電場對時間的偏導數項，采用向后差分的方式進行離散，即：

(10)

將(10)式代入(9)式，通過單元分析，可將(9)式寫成矩陣形式：

(11)

其中：

(12)

·NidV，

(13)

(14)

即可求得電磁響應的垂直分量.

1.4 源的加載

圖4 多輻射源加載示意圖(參考：孫懷鳳，2013)Fig.4 Schematic diagram of multi-radiation source loading(reference to Sun H F, 2013)

將電流密度直接施加到與電場的水平分量重合的單元棱邊上，如圖4所示.一般瞬變電磁場的正演均采用矩形脈沖進行激發，實際工作中，發射機不可能做到矩形脈沖發射，為了符合實際情況，在正演計算中，本文采用梯形波作為激發源.

2 多輻射源地空瞬變電磁法響應模擬

為了對比多輻射場源和單輻射場源的電磁響應，設計兩個不同的三維地質模型，并利用自主開發的瞬變電磁三維模擬軟件對其單輻射源和多輻射源下的電磁響應特征進行了模擬分析.地空電磁系統所采集的為垂直磁場分量或其感生電動勢.

2.1 模型一模擬結果

模型一：半空間電阻率為100 Ωm；異常體尺寸為200 m×200 m×30 m，電阻率1 Ωm，頂板埋深100 m；測線位于y=0 m處，飛行高度分別為10 m、30 m、50 m；分別計算在單源和雙源激發情況下的電磁響應.

(1)單一地質體、單一輻射源電磁響應特征

單輻射場源與單一三維模型立體圖如圖5a所示，其中源1長300 m，位于(-150，-150，0)至(150，-150，0)，供電電流15 A.沿測線(y=0 m)在不同飛行高度10 m、30 m、50 m的感生電動勢多測道圖分別如圖5b、c、d所示，由圖可知垂直分量多測道圖呈單峰分布，10 m飛行高度的感生電動勢幅值較大，而50 m飛行高度的感生電動勢幅值較小.

(2)單一地質體、多輻射源磁場響應特征

當兩個輻射源的電流方向相反時的地電模型如圖6a所示，源1與以上模型一致，源2位于(150，150，0)至(-150，150，0)，供電電流均為15 A.圖6b、c、d分別為沿測線(y=0)在不同飛行高度10 m、30 m、50 m的感生電動勢多測道圖，由圖可知感生電動勢多測道圖呈單峰分布，10 m飛行高度的感生電動勢幅值較大，而50 m飛行高度的感生電動勢幅值較小.同時，通過與圖5計算結果比較，可以得知：當兩個輻射源以相反方向的電流進行發射時，所激發的能量在空間具有一定的疊加作用，地下目標體的響應特征雖然彼此相似，但是兩個電流方向相反且相互平行的輻射源電磁響應幅值明顯大于單一發射源的情況，從而大大加強了探測信號強度.

2.2 模型二模擬結果

模型二：半空間電阻率為100 Ωm，異常體為兩個直立板目標體，尺寸均為300 m×200 m×50 m，頂板埋深均為100 m，電阻率均為1 Ωm，兩板間水平距離500 m，關于y軸對稱放置(如圖7a所示)；電源長均為500 m，電源距O點500 m，供電電流15 A.取飛行高度20 m，分別計算在單源、雙源和四源激發情況下的電磁響應.

(1)兩個地質體單輻射源垂直磁場分量響應特征

單輻射源激發模型如圖7a所示，激發源1長500 m，由(-250，-500，0)至(250，-500，0)，供電電流15 A，飛行高度20 m.圖7b所示的垂直磁場平面分布圖，說明由單源激發時，由垂直分量響應不能區分這兩個目標體.

(2)兩個地質體、兩個輻射源磁場垂直分量響應特征

兩個目標體水平距離500 m，模型參數不變，如圖8a所示，激發源1與圖7模型一致，激發源2長500 m，位于(250，500，0)至(-250，500，0)，供電電流15 A.當利用雙源同時激發時，取飛行高度30 m，在圖8b所示的垂直磁場平面圖中，出現了兩個異常體的外邊界輪廓，但是兩個異常體之間的邊界仍然顯得比較模糊.

(3) 兩個地質體、四個輻射源磁場垂直分量響應特征

兩個目標體模型參數不變，如圖9a所示，利用四個激發源同時進行激發(源長均為500 m，源1、源2與以上模型一致，源3位于(500，-250，0)至(500，250，0)，源4位于(-500，250，0)至(-500，250，0)，供電電流均為15A，飛行高度為30 m)，可以明顯看到兩個目標體分異明顯.當利用四個輻射源激發時，由圖9b所示的垂直磁場平面分布圖可知，兩個目標體的異常對稱，異常幅值較高，其內邊界的清晰度提高.通過與圖8所示的計算結果相比，可以認為四個發射源的垂直磁場響應比兩個發射源的響應效果更好，所反映的地質目標更真實.

3 結論

地空電磁法兼具探測深度大和探測效率高的優點，對于深部資源調查具有重要意義.在目前的實際應用中，地空電磁法一般單源激發，導致所接收到的信號較弱，難以獲得高精度解釋結果.本文基于自主研發的三維矢量有限元法正演模擬軟件對多輻射源地空瞬變電磁響應進行三維模擬，當采用單源激發時，由于僅在一個方向激發，場的幅值和分辨率受到限制，獲取的地質體的信息不全面，或者目標體并不能得到有效的識別.當采用多輻射場源作為地空電磁法的發射源時，能夠獲得不同角度的電磁場的輻射信息，從而獲得比單源激發更高的分辨率.特別是四個源同時進行激發時，對于兩個目標體有很好的識別能力.利用源的排列及電流方向等因素對信號影響的差異，合理布設電性源，可以達到加大勘探深度，提高多個目標體的分辨能力的目的.

圖5 不同飛行高度感生電動勢多測道圖(a) 模型三維立體圖； (b) 飛行高度10 m感生電動勢多測道圖； (c) 飛行高度30 m感生電動勢多測道圖； (d) 飛行高度50 m感生電動勢多測道圖.Fig.5 Multi-channel diagram of the dBz/dt with different flight heights(a) 3D view of model ； (b) Multi-channel diagram of dBz/dt with flight altitude 10 m； (c) Multi-channel diagram of dBz/dt with flight altitude 30 m; (d) Multi-channel diagram of dBz/dt with flight altitude 50 m.

圖6 不同飛行高度感生電動勢多測道圖(a) 模型三維立體圖； (b) 飛行高度10 m感生電動勢多測道圖； (c) 飛行高度30 m感生電動勢多測道圖； (d) 飛行高度50 m感生電動勢多測道圖.Fig.6 Multi-channel diagram of the dBz/dt with different flight heights(a) Model overhead view; (b) Multi-channel diagram of dBz/dt with flying height 10m; (c) Multi-channel diagram of dBz/dt with flying height 30m; (d) Multi-channel diagram of dBz/dt with flying height 50 m.

圖7 雙地質體、單一輻射源模型和垂直磁場分量平面圖(a) 模型三維立體圖； (b) 飛行高度20 m垂直磁場分量平面圖.Fig.7 Plane diagram of the Bz of two geological bodies with a single radiation source(a) 3D view of model； (b) Plane diagram of Bz with flight altitude 20 m.

圖8 雙地質體、雙輻射源模型和垂直磁場分量平面圖(a) 模型三維立體圖； (b) 飛行高度30 m垂直磁場分量平面圖.Fig.8 Plane diagram of the Bz of two geological bodies with two radiation sources(a) 3D view of model； (b) Plane diagram of Bz with flight altitude 20 m.

圖9 雙地質體、四個輻射源模型和垂直磁場分量平面圖(a) 模型三維立體圖； (b) 飛行高度20 m垂直磁場分量平面圖.Fig.9 Plane diagram of the Bz of two geological bodies with four radiation sources(a) 3D view of model； (b) Plane diagram of Bz with flight altitude 30 m.