李貅， 胡偉明， 薛國(guó)強(qiáng)

1 長(zhǎng)安大學(xué)地質(zhì)工程與測(cè)繪學(xué)院， 西安 710054 2 中國(guó)科學(xué)院礦產(chǎn)資源研究重點(diǎn)實(shí)驗(yàn)室， 中國(guó)科學(xué)院地質(zhì)與地球物理研究所， 北京 100029 3 中國(guó)科學(xué)院大學(xué)地球與行星科學(xué)學(xué)院， 北京 100049 4 中國(guó)科學(xué)院地球科學(xué)研究院， 北京 100029 5 西北有色地質(zhì)礦業(yè)集團(tuán)有限公司， 西安 710054

0 引言

電磁法通過(guò)觀測(cè)天然或人工源激發(fā)的電磁響應(yīng)獲取地下電性信息(何繼善和薛國(guó)強(qiáng)，2018;底青云等，2019).傳統(tǒng)地面電磁法觀測(cè)方式在沙漠、戈壁、山地和濕地等地區(qū)無(wú)法高效地開(kāi)展工作，航空電磁法可發(fā)揮替代作用(圖1).航空電磁法根據(jù)發(fā)射源布設(shè)方式可分為全航空電磁法和地空電磁法，其中全航空電磁法采集效率高，但由于受到載重限制，發(fā)射功率有限.地空電磁法將發(fā)射系統(tǒng)放置于地面，接收系統(tǒng)布設(shè)在飛行平臺(tái)上(圖2)，充分地結(jié)合了傳統(tǒng)地面電磁法和全航空電磁法的優(yōu)點(diǎn)，可實(shí)現(xiàn)探測(cè)深度和工作效率之間的平衡.

圖1 航空電磁探測(cè)系統(tǒng)示意圖Fig.1 Schematic diagram of airborne electromagnetic detection system

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

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

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

1 瞬變電磁場(chǎng)三維矢量有限元法原理

1.1 邊值問(wèn)題

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

(1)

(2)

其中，E為電場(chǎng)強(qiáng)度，H為磁場(chǎng)強(qiáng)度，σ為介質(zhì)電導(dǎo)率，μ0為介質(zhì)的磁導(dǎo)率，Js為源電流密度.對(duì)(1)式兩端取旋度：

(3)

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

(4)

整理可得：

(5)

根據(jù)電磁場(chǎng)理論，在無(wú)源區(qū)的兩種導(dǎo)電介質(zhì)界面，電磁場(chǎng)滿(mǎn)足如下四個(gè)邊界條件：

(6)

式中，n為兩種介質(zhì)界面處的單位法向分量，B為磁感應(yīng)強(qiáng)度，D為電位移矢量.

對(duì)于無(wú)窮遠(yuǎn)邊界，電場(chǎng)和磁場(chǎng)在無(wú)窮遠(yuǎn)邊界上的切向分量為零，即：

(7)

1.2 矢量變分方程

根據(jù)加權(quán)余量法，電場(chǎng)控制方程的加權(quán)余量方程為

(8)

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

=0，

(9)

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

1.3 矢量有限元法

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

采用六面體單元進(jìn)行網(wǎng)格剖分，在單元中棱邊與節(jié)點(diǎn)的關(guān)系如圖3所示.在六面體單元中，將電場(chǎng)切向分量賦予各個(gè)單元的棱邊上.采用Whitney型插值基函數(shù)(李賀，2016)，使得插值基函數(shù)的散度為0，而旋度不等于0.因此，在計(jì)算區(qū)域存在介質(zhì)不連續(xù)性時(shí)，電場(chǎng)強(qiáng)度的切向分量連續(xù)，法向分量不連續(xù)，可以有效地避免“偽解”的現(xiàn)象.

在(9)式中，含有電場(chǎng)對(duì)時(shí)間的偏導(dǎo)數(shù)項(xiàng)，采用向后差分的方式進(jìn)行離散，即：

(10)

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

(11)

其中：

(12)

·NidV，

(13)

(14)

即可求得電磁響應(yīng)的垂直分量.

1.4 源的加載

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

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

2 多輻射源地空瞬變電磁法響應(yīng)模擬

為了對(duì)比多輻射場(chǎng)源和單輻射場(chǎng)源的電磁響應(yīng)，設(shè)計(jì)兩個(gè)不同的三維地質(zhì)模型，并利用自主開(kāi)發(fā)的瞬變電磁三維模擬軟件對(duì)其單輻射源和多輻射源下的電磁響應(yīng)特征進(jìn)行了模擬分析.地空電磁系統(tǒng)所采集的為垂直磁場(chǎng)分量或其感生電動(dòng)勢(shì).

2.1 模型一模擬結(jié)果

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

(1)單一地質(zhì)體、單一輻射源電磁響應(yīng)特征

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

(2)單一地質(zhì)體、多輻射源磁場(chǎng)響應(yīng)特征

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

2.2 模型二模擬結(jié)果

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

(1)兩個(gè)地質(zhì)體單輻射源垂直磁場(chǎng)分量響應(yīng)特征

單輻射源激發(fā)模型如圖7a所示，激發(fā)源1長(zhǎng)500 m，由(-250，-500，0)至(250，-500，0)，供電電流15 A，飛行高度20 m.圖7b所示的垂直磁場(chǎng)平面分布圖，說(shuō)明由單源激發(fā)時(shí)，由垂直分量響應(yīng)不能區(qū)分這兩個(gè)目標(biāo)體.

(2)兩個(gè)地質(zhì)體、兩個(gè)輻射源磁場(chǎng)垂直分量響應(yīng)特征

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

(3) 兩個(gè)地質(zhì)體、四個(gè)輻射源磁場(chǎng)垂直分量響應(yīng)特征

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

3 結(jié)論

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

圖5 不同飛行高度感生電動(dòng)勢(shì)多測(cè)道圖(a) 模型三維立體圖； (b) 飛行高度10 m感生電動(dòng)勢(shì)多測(cè)道圖； (c) 飛行高度30 m感生電動(dòng)勢(shì)多測(cè)道圖； (d) 飛行高度50 m感生電動(dòng)勢(shì)多測(cè)道圖.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 不同飛行高度感生電動(dòng)勢(shì)多測(cè)道圖(a) 模型三維立體圖； (b) 飛行高度10 m感生電動(dòng)勢(shì)多測(cè)道圖； (c) 飛行高度30 m感生電動(dòng)勢(shì)多測(cè)道圖； (d) 飛行高度50 m感生電動(dòng)勢(shì)多測(cè)道圖.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 雙地質(zhì)體、單一輻射源模型和垂直磁場(chǎng)分量平面圖(a) 模型三維立體圖； (b) 飛行高度20 m垂直磁場(chǎng)分量平面圖.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 雙地質(zhì)體、雙輻射源模型和垂直磁場(chǎng)分量平面圖(a) 模型三維立體圖； (b) 飛行高度30 m垂直磁場(chǎng)分量平面圖.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 雙地質(zhì)體、四個(gè)輻射源模型和垂直磁場(chǎng)分量平面圖(a) 模型三維立體圖； (b) 飛行高度20 m垂直磁場(chǎng)分量平面圖.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.