基于三維有限元的航空大地電磁傾子響應(yīng)特征

劉彥濤，彭莉紅，孫棟華，張偉盟，王海紅

(1.核工業(yè)航測遙感中心，河北 石家莊 050002； 2.中核集團(tuán)鈾資源地球物理勘查技術(shù)(重點實驗室)，河北 石家莊 050002； 3.河北省航空探測與遙感技術(shù)重點實驗室，河北 石家莊 050002)

0 引言

大地電磁法利用天然電磁場作為探測場源，具有探測深度大、設(shè)備輕便等優(yōu)點[1]。為對高山、荒漠、沼澤地等勘探人員無法到達(dá)區(qū)域進(jìn)行地球物理勘探，可選方法之一為運用航空器搭載測量儀器對這些區(qū)域進(jìn)行探測，也因此航空大地電磁近些年得到了快速發(fā)展[2-5]。與傳統(tǒng)大地電磁方法原理相同，航空大地電磁利用天然電磁場作為場源，接收由于地下地質(zhì)體橫向不均勻激勵的垂直磁場分量，結(jié)合地面接收到的磁場水平分量，換算傾子參數(shù)，反演地下地質(zhì)體電阻率分布，從而達(dá)到探測目的。遠(yuǎn)方天然電磁場可看作平面電磁波垂直入射到地面，只有當(dāng)?shù)刭|(zhì)體存在橫向電性變化時，才會激勵起垂直磁場[6]。不同于傳統(tǒng)大地電磁法，由于航空大地電磁空中探測無法接地，只能采集磁場分量，通過磁場分量可換算得到傾子參數(shù)，再利用傾子參數(shù)對地下地質(zhì)體進(jìn)行地球物理解釋。

傾子是對地質(zhì)體電性橫向不均勻性的一種反映[7]，由于磁場垂直分量強(qiáng)度小且受干擾大，早期大地電場法勘探中沒有得到成功應(yīng)用，1989年Labson提出遠(yuǎn)參考技術(shù)對原始數(shù)據(jù)進(jìn)行降噪處理，并通過優(yōu)化磁場接收裝置，成功使用傾子參數(shù)進(jìn)行了地球物理解釋[8]。近些年對大地電磁傾子響應(yīng)的研究較多[9-12]，尤其是利用數(shù)值模擬算法模擬傾子對典型地電模型的響應(yīng)特征，分析其響應(yīng)規(guī)律[13-17]，研究覆蓋層、橫向電阻率差異大小、地形、測線高度等因素對傾子的影響[18-20]。前人數(shù)值研究中所選測線大多垂直巖性分界面，使異常體傾子響應(yīng)理想化，而在實際探測中，不能保證設(shè)計測線與所有橫向地電斷面垂直。本文首先簡單介紹了航空大地電磁勘探原理以及航空大地電磁ZTEM探測系統(tǒng)，利用三維有限元方法對測區(qū)磁場分量分布特征進(jìn)行模擬分析，之后建立典型地電模型，研究了當(dāng)測線與橫向地電斷面斜交時的傾子響應(yīng)特征。

1 方法原理

天然電磁場以平面波的形式垂直入射到地下，入射的平面電磁波中不存在磁場垂直分量，當(dāng)?shù)叵碌刭|(zhì)體存在橫向電性變化時，將會激勵垂直磁場，這個垂直磁場與兩個水平磁場存在復(fù)線性關(guān)系：

Hz(ω)=Tx(ω)Hx(ω)+Ty(ω)Hy(ω) ，

(1)

式中：Hx(ω)、Hy(ω)、Hz(ω)分別為不同頻率x、y、z磁場分量，Tx(ω)、Ty(ω)為傾子系數(shù)。傾子是復(fù)數(shù)值且隨頻率變化。傾子參數(shù)與地下地質(zhì)體電性分布相關(guān)而不與激發(fā)電磁波極化方向相關(guān)，所以傾子可以作為解釋地下電性分布的參數(shù)。

在式(1)兩邊分別乘以Hx(ω)、Hy(ω)的共軛復(fù)數(shù)值[6]，構(gòu)建方程組：

(2)

通過方程組(2)可求取傾子值Tx、Ty：

(3)

2 ZTEM系統(tǒng)

ZTEM系統(tǒng)是目前全球應(yīng)用較為廣泛的商用航空大地電磁系統(tǒng)，系統(tǒng)主要由兩部分組成，如圖1。

圖1 ZTEM系統(tǒng)Fig.1 System of ZTEM

空中接收線圈用于接收磁場Hz分量，線圈直徑7.4 m，考慮到直升機(jī)對接收信號的影響，同時保證飛行員的可操控性，線圈吊掛于直升機(jī)下方85 m處。線圈上安置3個GPS天線，實時監(jiān)測線圈姿態(tài)，用于消除水平和垂直電磁場交叉耦合的效應(yīng)。數(shù)據(jù)采集過程中控制線圈離地高度約100 m，飛機(jī)飛行速度80 km/h。

地面接收線圈主要是由3個正交的邊長為 3.05 m的方形接收線圈組成，安置在工作區(qū)25 km范圍內(nèi)且人文干擾相對較小的固定地點，用來觀測水平磁場分量Hx和Hy?；? km范圍內(nèi)最好無高壓線通過，并且沒有無線電基站、信號發(fā)射站等發(fā)射無線信號的裝置。

圖2 頻率-趨膚深度曲線Fig.2 Curve of frequency-skin depth

從圖中可以看出，ZTEM系統(tǒng)6個頻率的趨膚深度線性增加，即ZTEM所選用的6個探測頻率使得6個頻率的趨膚深度大致呈均勻等間距分布，這種分布形式利于后期數(shù)據(jù)成圖與反演計算。

3 三維有限元算法

本文使用商用軟件COMSOL Multiphysics進(jìn)行正演模擬。該軟件是基于三維有限元算法的多物理場模擬軟件，使用AC/DC模塊建立典型三維地質(zhì)模型，對大地電磁響應(yīng)進(jìn)行模擬，分析傾子響應(yīng)特征。

(4)

綜合使用磁絕緣和理想磁導(dǎo)體邊界條件，設(shè)置垂直場源的邊界為理想磁導(dǎo)體邊界，平行場源的邊界為磁絕緣邊界，從而達(dá)到模擬無限域的效果。使用自由四面體對模型進(jìn)行網(wǎng)格剖分，對異常體網(wǎng)格加密，提高求解精度，模型剖分效果如圖3所示。

圖3 模型網(wǎng)格剖分示例Fig.3 Example of model meshing

依據(jù)ZTEM系統(tǒng)所用頻率，數(shù)值模擬中使用頻率為25、37、75、150、300、600 Hz，模型背景電阻率100 Ω·m，充分考慮模擬邊界效應(yīng)以及最低頻率最大探測深度應(yīng)小于模型高度，設(shè)置模型長3 000 m，寬3 000 m，高1 400 m。為突出傾子對異常響應(yīng)特點，異常體上方不設(shè)置覆蓋層。

4 大地電磁三維正演模擬

4.1 全區(qū)域磁場三分量分布特征

航空大地電磁ZTEM系統(tǒng)數(shù)據(jù)采集過程中，基站固定于測區(qū)某個人文干擾低的地方采集磁場Hx、Hy分量，而不是沿測線進(jìn)行數(shù)據(jù)采集。這種采集方式基于磁場Hx、Hy分量全區(qū)變化可忽略的假設(shè)，為驗證上述假設(shè)，建立模型對大地電磁磁場三分量進(jìn)行正演模擬，分析其分布特征。

如圖4所示建立三維模型，模型長3 000 m，寬3 000 m，高1 400 m，背景電阻率為100 Ω·m，設(shè)置2個相同尺寸不同電阻率的異常體，單個異常體長1 500 m，寬600 m，高1 000 m，異常體電阻率值分別為10 Ω·m與1 000 Ω·m，選取沿紅線方向的磁場三分量數(shù)據(jù)，進(jìn)行成圖分析。

圖4 復(fù)雜地質(zhì)模型示意Fig.4 Schematic diagram of complex geological modle

所選測線既有沿巖性分界面的情況，又有垂直巖性分界面的情況，相對具有代表性。選取37 Hz與300 Hz兩個頻率磁場三分量數(shù)據(jù)進(jìn)行成圖(圖5)，兩幅圖反映的規(guī)律類似：所選測線方向磁場水平分量Hx、Hy變化很小，基本可以忽略不計，磁場垂直分量Hz變化非常明顯，在測線垂直相交的巖性分界面上出現(xiàn)極值點；由此驗證了在實際數(shù)據(jù)采集過程中，測區(qū)固定地點采集Hx、Hy分量、沿測線采集Hz分量的采集方式的可行性。同時，從圖中還可以看出2個頻率出現(xiàn)極值點的位置相同，但極值絕對值大小不同，高頻磁場Hz分量極值絕對值大于低頻磁場Hz分量極值絕對值。

圖5 2個典型頻率的磁場三分量變化趨勢Fig.5 Variation trend of three components of magnetic field at two typical frequencies

4.2 低阻體傾子響應(yīng)特征

利用數(shù)值模擬對傾子響應(yīng)特征進(jìn)行分析的文章較多，所研究的響應(yīng)特征皆為測線垂直巖性分界面的情況，但實際應(yīng)用中測線設(shè)計不能每次都做到完全垂直巖性分界面。下面利用有限元法建立三維模型，研究測線與巖性分界面斜交時的傾子響應(yīng)特征。

如圖6所示，三維模型長3 000 m，寬3 000 m，高1 400 m，背景電阻率設(shè)為100 Ω·m；異常體長1 500 m，寬800 m，高1 000 m，異常體電阻率設(shè)為10 Ω·m，坐標(biāo)原點為模型上表面中心點。利用三維有限元方法對大地電磁法進(jìn)行數(shù)值模擬，研究傾子沿巖性分界面斜交測線低阻異常體的響應(yīng)特征。

圖6 三維模型示意Fig.6 3D model sketch

選擇3條具有代表性的測線，如圖7所示(圖中為模型俯視圖)：測線a與x方向巖性分界面斜交；測線b過x方向和y方向巖性分界面交點，沿異常體對角線方向；測線c與y方向巖性分界面斜交。沿3條測線分別繪制6個頻率的傾子Tx與傾子Ty曲線(圖8～圖10，圖中虛線代表測線與巖性分界面相交位置)，分析傾子響應(yīng)特征。數(shù)值模擬和實際應(yīng)用中傾子實部響應(yīng)規(guī)律較傾子虛部響應(yīng)規(guī)律好，文中所提傾子皆指傾子實部。

圖7 測線走向示意Fig.7 Sketch map of survey line

圖8為沿測線a傾子Tx與傾子Ty曲線，可以看出傾子Ty響應(yīng)幅值高出傾子Tx響應(yīng)幅值一個數(shù)量級，且傾子Ty響應(yīng)極值點位置對巖性分界面位置顯示較為準(zhǔn)確。傾子Tx響應(yīng)幅值較低,高頻曲線對巖性分界面無反映，低頻曲線對巖性分界面的響應(yīng)特征與傾子Ty的響應(yīng)特征類似，但曲線極值點與巖性分界面位置存在偏差，傾子響應(yīng)極值點寬于巖性分界點。測線a沿y軸方向存在巖性變化，理論上表現(xiàn)為傾子Ty響應(yīng)，由于測線非完全垂直于巖性分界面，存在微弱傾子Tx響應(yīng)，該響應(yīng)對巖性分界面的表現(xiàn)較差。

圖8 測線a傾子低阻響應(yīng)特征Fig.8 Low-resistance response characteristics of tipper on survey line a

圖9顯示，沿測線b傾子Tx響應(yīng)幅值略低于傾子Ty響應(yīng)幅值，各頻率對于巖性分界面皆有響應(yīng)，傾子Tx曲線極值點對巖性分界面的位置反映較為準(zhǔn)確，但Ty曲線極值點對巖性分界面的位置反映存在偏差，曲線極值點范圍窄于實際巖性分界面位置。

圖9 測線b傾子低阻響應(yīng)特征Fig.9 Low-resistance response characteristics of tipper on survey line b

圖10顯示，沿測線c傾子Tx響應(yīng)幅值高于傾子Ty響應(yīng)幅值，且前者曲線平滑度優(yōu)于后者，傾子Tx對巖性分界面的位置反映較為準(zhǔn)確。傾子Ty響應(yīng)幅值較低，高頻曲線對巖性分界面反映微弱，低頻曲線極值點對巖性分界面位置反映存在偏差，曲線極值點反映的異常體寬度窄與實際異常體寬度。

圖10 測線c傾子低阻響應(yīng)特征Fig.10 Low-resistance response characteristics of tipper on survey line c

縱向?qū)Ρ葓D8、圖9、圖10，傾子Tx響應(yīng)曲線反映x軸方向巖性變化，隨著測線偏離x軸方向越大，傾子Tx響應(yīng)幅值越低，對于巖性分界面的反映能力越差，曲線極值點過渡越平滑；而傾子Ty響應(yīng)曲線反映y軸方向巖性變化，隨著測線偏離y軸方向越大，傾子Ty響應(yīng)幅值越低，對于巖性分界面的反映能力越差，曲線極值點過渡越平滑。

4.3 高阻體傾子響應(yīng)特征

將圖6模型異常體電阻率設(shè)置為1 000 Ω·m，模型背景圍巖電阻率仍為100 Ω·m，坐標(biāo)原點為模型上表面中心點。利用三維有限元方法對大地電磁法進(jìn)行數(shù)值模擬，繪制沿圖7測線傾子響應(yīng)曲線，研究傾子沿巖性分界面斜交測線高阻異常體響應(yīng)特征。

圖11～圖13分別對應(yīng)3條測線傾子響應(yīng)曲線。從圖11可以看出，沿測線a傾子Tx響應(yīng)幅值低于傾子Ty響應(yīng)幅值，兩者各個頻率曲線極值對異常體邊界皆有反映，傾子Tx極值點曲線過渡相對平滑，極值點位置與巖性分界點位置略微存在偏差，但偏差不大。圖12表明，沿測線b傾子Tx響應(yīng)與傾子Ty響應(yīng)曲線類似且幅值相近，傾子Tx與傾子Ty曲線極值點對巖性分界面的位置反映皆較為準(zhǔn)確。圖13顯示，沿測線c傾子Tx響應(yīng)幅值高于傾子Ty響應(yīng)幅值，傾子Tx與傾子Ty各頻率曲線極值點對巖性分界面的位置反映皆較為準(zhǔn)確。

對比圖11、圖12、圖13可以看出，相同地質(zhì)模型下，隨測線偏離x軸方向越大，測線方向越接近y軸，傾子Tx響應(yīng)幅值越低，傾子Ty響應(yīng)幅值越高。對比低阻體響應(yīng)曲線可以看出，高阻體傾子響應(yīng)曲線較低阻體響應(yīng)曲線光滑，響應(yīng)幅值較低的測線上，傾子高頻曲線不存在對異常體分界面無反映的狀態(tài)。同時，對于低阻異常體，傾子曲線在巖性分界面形成極值，在異常體中段無巖性變化位置，傾子曲線存在值為零的水平段，如圖10中傾子Tx響應(yīng)，而高阻體模型中，傾子響應(yīng)曲線在兩極值間不存在傾子為零的水平段。實際上，在高阻異常體數(shù)值模擬中發(fā)現(xiàn)，只有當(dāng)兩邊巖性分界面最夠遠(yuǎn)時，傾子響應(yīng)曲線才會在異常體中心位置出現(xiàn)傾子為零的水平段。

圖11 測線a傾子高阻響應(yīng)特征Fig.11 High-resistance response characteristics of tipper on survey line a

圖12 測線b傾子高阻響應(yīng)特征Fig.12 High-resistance response characteristics of tipper on survey line b

5 結(jié)論與討論

1)航空大地電磁只采集磁場分量，在測區(qū)范圍內(nèi)磁場水平分量變化很小，只有地質(zhì)體存在橫向不均勻性時才能激發(fā)磁場垂直分量，所以航空大地電磁數(shù)據(jù)采集過程中，只需磁場垂直分量沿測線進(jìn)行采集，水平分量固定地點測量即可。

2)傾子作為航空大地電磁數(shù)據(jù)解釋中重要參數(shù)，在巖性橫向分界面上表現(xiàn)為極值異常。當(dāng)測線與巖性分界面斜交時，傾子Tx與傾子Ty對異常分界面皆有反映。隨測線偏離x軸方向越大，測線方向越接近y軸，傾子Tx響應(yīng)幅值越低，傾子Ty響應(yīng)幅值越高，反之亦然。

3)傾子Tx對x軸方向巖性變化反映最好，當(dāng)測線與x軸方向巖性變化分界面斜交時，傾子Ty對巖性分界面存在微弱反映，且當(dāng)異常體為低阻體時，傾子Ty對巖性分界面位置反映存在偏差；同樣，傾子Ty對y軸方向巖性變化反映最好，當(dāng)測線與y軸方向巖性變化分界面斜交時，傾子Tx對巖性分界面存在微弱反映，且當(dāng)異常體為低阻體時，傾子Tx對巖性分界面位置反映存在偏差。

4)相同地質(zhì)模型，不同異常電阻率條件下，當(dāng)測線與y軸方向巖性變化分界面斜交時，高阻異常體傾子Tx曲線較低阻異常體傾子Tx曲線平滑，且低阻異常體傾子Tx低頻曲線對巖性分界面反映較差。

值得注意的是數(shù)值模擬為無噪聲影響理想條件下的大地電磁響應(yīng)特征，傾子曲線中高頻響應(yīng)幅值高于低頻響應(yīng)幅值，實際野外數(shù)據(jù)采集過程中，頻率較高時，由于采樣窗口窄，數(shù)據(jù)信噪較低。本文只討論了傾子實部響應(yīng)特征，可進(jìn)一步對傾子虛部與相位響應(yīng)特征進(jìn)行研究，使用傾子多個參數(shù)對地質(zhì)異常進(jìn)行解釋，提高航空大地電磁探測準(zhǔn)確度。