嚴(yán)家斌， 胡 濤， 趙 玄， 皇祥宇

(1. 中南大學(xué) 地球科學(xué)與信息物理學(xué)院，長沙 410083；2. 有色資源與地質(zhì)災(zāi)害探查湖南省重點實驗室，長沙 410083 )

0 引言

電磁法勘探從卡尼亞電阻率的出現(xiàn)到如今已發(fā)展近70年，基于在地表測量正交分量求取阻抗消除場源因素的理論已經(jīng)十分成熟。雖已有諸多學(xué)者對可控源電磁法的垂直場分量Ez、Hz進行了理論推導(dǎo)和研究[1]，同時電磁場地表的Hz也對構(gòu)造的水平不均勻性有較好的反映。但對地下空間場分量的分布特征淺層構(gòu)造的研究相對較少，這是因為地下垂直電場Ez的信號十分微弱，數(shù)量級遠小于水平電場，增大了觀測時的誤差。雖然在地表電磁勘探中并無Ez分量，但若同時測量鉆井中的多道電場，分析電磁場中多道電場的變化的空間關(guān)系，對于油氣藏監(jiān)測等有實際的借鑒性意義。

通過前人的研究可以總結(jié)出測量Ez的方法，包括直接測量法和間接測量法兩種。直接測量法就是利用鉆井技術(shù)的井中(或海洋)電極測量Ez，早期Jang等[2]研究了海洋沉積物中埋藏的深阻層的垂直電磁響應(yīng)，結(jié)果表明雖然垂直電場的信號振幅比水平場小得多，但垂直發(fā)射機產(chǎn)生的垂直電流對電阻層是極其敏感的；Sharqawy等[3]利用地井鉆孔技術(shù)進行電磁測量，說明測量垂直電場在技術(shù)上是可行的；王智等[4]對頻率域地井垂直電磁法響應(yīng)特征進行了數(shù)值模擬，說明利用鉆井技術(shù)測量的垂直電場能夠反映地電結(jié)構(gòu)垂向變化。而間接測量法就是需要通過磁場間接測量電場，Rocroi[5]、Takacs等[6]首先提出探測可控源電磁法中的垂直電場，提高了對深度目標(biāo)體的分辨率。GFZ(德國赫姆霍茲地學(xué)研究中心)的Tietze K等[7]推導(dǎo)了可控源電磁法的垂直電場，并討論其對地電模型的電阻率橫向與縱向變化敏感度。L?seth等[8]在海上電磁測量中采用水平電偶極源和水平磁偶極源觀測垂直電磁場分量，對推導(dǎo)出的Ez分量響應(yīng)進行分析，結(jié)果表明利用Ez分量來監(jiān)測儲油層的變化，能獲得較好的效果。Schaller[9]采集了荷蘭Schoonebeek陸上油田模型的合成電磁數(shù)據(jù)證實了雖然Ez的振幅小，但它比水平分量對儲層的反映更敏感。一系列研究表明無論是通過直接或間接測量得出Ez，在技術(shù)和理論上都是可行的。但利用可控源鉆孔技術(shù)對Ez分量進行監(jiān)測會遇到不可避免的問題，那就是場源難以做到一直發(fā)射，未來可能會出現(xiàn)利用天然場源監(jiān)測多道Ez和Hz的空間變化異常分布。

筆者在皇祥宇博士編寫的三維矢量有限元法程序的基礎(chǔ)上，重點研究了不均勻性引起的z分量的空間變化規(guī)律，對不同模型的垂直電場變化規(guī)律進行了數(shù)值模擬，量化電磁場垂直電場的響應(yīng)特征。

1 矢量有限元理論基礎(chǔ)

隨著計算機技術(shù)的不斷發(fā)展和學(xué)者對三維正反演算法的大量研究，對頻率域電磁法進行三維正演已不再困難。三維電磁模擬時，基于節(jié)點的標(biāo)量有限元由于自由度賦予節(jié)點難以模擬突變和間斷，在電性界面突變處難以獲得高精度的場分量?；谑噶縿莸娜S有限元方法通過進一步求取微分來獲取電場和磁場，需要采取較細(xì)的網(wǎng)格獲得高精度Ez。因此為獲得高精度的垂直分量Ez、Hz，筆者采用矢量有限元來進行三維正演。

1.1 變分問題

根據(jù)電磁場理論，忽略位移電流時的頻率域麥克斯韋方程組為

▽×E=iωμH

▽×H=σE

(1)

式中:ω為角頻率；μ為磁導(dǎo)率；ε為介電常數(shù)；σ為電導(dǎo)率。

在不考慮磁導(dǎo)率變化的情況下[10]，將方程組(1)求旋度消除H，得到三維雙旋度方程：

▽×▽×E-k2E=0

(2)

式中：k為波數(shù)，k2=iωμσ。三維情況下異常體界面有場突變，不能計算邊界處的電場散度，電場各分量相互耦合，無法單獨分離出來，這時電場可表示成三分量形式：

E=Exex+Eyey+Ezez

(3)

式中：ex、ey、ez分別表示三個方向的單位矢量。設(shè)三維研究區(qū)域為Ω，區(qū)域邊界為Γ，根據(jù)加權(quán)余量法計算得到加權(quán)余量方程式(4)。

圖1 矢量六面體單元中節(jié)點與棱邊編號Fig.1 Node and edge number in the vector hexahedral cell

(4)

式中：δE為E的變分。邊界切向分量采用一維層狀模型計算電場值作為狄利克雷邊界條件，代入可形成最后的方程[11](S表示計算區(qū)域的底界面)：

(5)

1.2 網(wǎng)格離散和矢量單元分析

矢量有限元由于其基函數(shù)的無散性和允許場法向分量的突變，直接計算電場具有較高的精度[12]。為將研究區(qū)域離散化，對模擬區(qū)域進行矩形網(wǎng)格剖分，每個單元節(jié)點和棱邊編號如圖1所示。

(6)

(7)

(8)

組裝所有單元的系數(shù)矩陣Ke，得到有限元的線性方程組

KE=b

(9)

式中：K表示總體系數(shù)矩陣，為大型稀疏復(fù)對稱矩陣；E=(E1,E2,…,En)；n為所有棱邊總數(shù)，表示要求解電場；b由強加邊界條件決定。通過賦大數(shù)法[13]將計算區(qū)域頂界面邊界條件代入式(9)，求解線性方程組，這樣就能夠得到網(wǎng)格單元棱邊上的Ex、Ey、Ez的值。

三維正演速度的關(guān)鍵在于大型稀疏線性方程組的求解[14]。由于矢量有限元離散后得到的剛度矩陣為病態(tài)矩陣，求解方程組時必須采用預(yù)處理技術(shù)以提高收斂速度。經(jīng)過數(shù)值模擬對比發(fā)現(xiàn)，三維下隨著網(wǎng)格剖分不均勻及頻率降低，需選擇更穩(wěn)定的預(yù)條件，因此采用穩(wěn)定的SSOR預(yù)處理技術(shù)[15]。為了進一步提高收斂速度，采用一維層狀電場作為迭代解法的初始解，可有效提高收斂速度，同時本文的模擬頻率較高，避免了低感應(yīng)數(shù)(低頻難以收斂)等問題的出現(xiàn)[16]。

2 模型分析

2.1 單個異常體模型

電阻率為100 Ω·m的均勻半空間中含有一個電阻率為1 Ω·m(或10 000 Ω·m)的異常體，異常體幾何尺寸為30 m×30 m×45 m，頂面埋深40 m(如圖2)。研究區(qū)域大小為3 000 m× 3 000 m×4 500 m，采用非均勻矩形網(wǎng)格剖分，x、y、z方向網(wǎng)格數(shù)為33×33×33(z方向含5層空氣層，電阻率108Ω·m) ，模擬測點數(shù)為34個，記錄頻點數(shù)為25個，異常體中心在地表投影坐標(biāo)為(0，0，0)。各方向網(wǎng)格剖分如下(單位m)：

x=[200 180 160 140 120 100 80 60 40 20 10 10 10 5 5 5 5 5 5 5 10 10 10 20 40 60 80 100 120 140 160 180 200]

y=[200 180 160 140 120 100 80 60 40 20 10 10 10 5 5 5 5 5 5 5 10 10 10 20 40 60 80 100 120 140 160 180 200]

z=[100 80 60 40 20 5 5 5 5 5 5 5 5 5 5 5 5 5 10 10 10 10 10 10 20 30 40 50 60 70 80 90 100]

2.1.1 低阻異常體

設(shè)定觀測頻率為3 500 Hz，模擬中心測點(0 m，0 m，0 m)的Hz、Ez沿z軸的變化規(guī)律(圖3)。Hz幅值變化規(guī)律體現(xiàn)了該頻率探測深度范圍內(nèi)介質(zhì)電性橫向變化的電磁響應(yīng)。分析圖3(a)可以看出，Hz分量的曲線光滑流暢，在穿越異常體模型界面時不會發(fā)生很明顯的變化。觀察圖3(b)能發(fā)現(xiàn)，地表(地-空界面)上流向空氣的電流密度為零，故Ez分量在z=0 m分界面處很小接近于零。沿z軸不斷向下，異常體的影響逐漸顯現(xiàn)，Ez分量的幅值不斷增加，并在臨近異常體模型的上界面(z=40 m)處達到極值，在穿過異常體上界面幅值發(fā)生先升高后降低的突變，進入異常體后迅速減小，然后趨于平緩，當(dāng)?shù)竭_異常體的下界面(z=85 m)處Ez分量的幅值又發(fā)生突變，達到極值，且上界面處極值數(shù)值大于下界面。在低阻體中，Ez分量均較小，這可由歐姆微分定理可知，電流的法向分量連續(xù)，必然造成z分量的突變。

因此利用Ez分量幅值在電性界面發(fā)生突變的性質(zhì)我們可以判定異常體的上下分界面及厚度。對比Hz、Ez分量，不難發(fā)現(xiàn)Ez分量在分辨異常體的上下界面上優(yōu)于Hz分量。

為探討Ez分量對地下異常的影響范圍，分析地表多個測點上的Ez振幅測深變化。測點分布如圖4所示，分別為中心點A(0 m，0 m，0 m)、內(nèi)部點B(-5 m，0 m，0 m)、棱邊點C(-15 m，0 m，0 m)和外部點D(-22.5 m，0 m，0 m)。

從圖5可以發(fā)現(xiàn)，位于異常體內(nèi)部的測點曲線大致相似，在異常體界面位置，棱邊點Ez振幅變化較大異常體中心處較小。當(dāng)測點為外部點D時，Ez曲線變化平緩，在上下界面處不發(fā)生突變。

圖2 單個異常體模型網(wǎng)格剖分及示意圖Fig.2 Grid division and schematic diagram of a single anomalous body model(a)單個異常體模型網(wǎng)絡(luò)剖分圖；(b)單個異常體模型示意圖

圖3 中心點A處Hz、Ez幅值Fig.3 The Hz & Ez amplitude of center measurement point A(a)Hz分量；(b)Ez分量

圖4 地表測點分布圖Fig.4 Measuring point distribution

圖5 多個測點Ez振幅Fig.5 The Ez amplitude of measurement points

圖6 Ez分量振幅平面圖Fig.6 The Ez component amplitude plan

圖7 歸一化處理測點分布Fig.7 Normalization processing point

圖8 歸一化處理結(jié)果圖Fig.8 Normalization processing results

為了更加突出Ez分量對異常體模型上下界面的反映作用，圖6給出了頻率為3 500 Hz時Ez分量振幅在剖面y=0 m上的平面圖(虛線框表示異常體的位置)。從圖中可看出，Ez分量幅值在異常體模型上界面(z=40 m)處和下界面(z=85 m)處呈現(xiàn)高值現(xiàn)象，上界面高值大于下界面，符合上述Ez分量在垂直方向上的變化規(guī)律，驗證了垂直分量Ez對異常體上下界面的有效識別作用。

為驗證Ez分量計算結(jié)果的合理性，我們采用類似可控源大地電磁法中的處理方法(圖7)，將上下兩個測點的Ez1、Ez2數(shù)據(jù)相除，用它們的比值作為歸一化處理結(jié)果，獲得穩(wěn)定的傳輸函數(shù)。通過將Ez1和Ez2相比發(fā)現(xiàn)(圖8)，曲線基本上在數(shù)值1附近波動，波動最劇烈處也在誤差允許范圍內(nèi)，說明我們對Ez數(shù)據(jù)的求取是可行的。

圖9 中心點A處Hz、Ez幅值Fig.9 The Hz & Ez amplitude of center measurement point A(a)Hz分量；(b)Ez分量

圖10 多個測點Ez振幅Fig.10 The Ez amplitude of measurement points

圖11 Ez分量振幅平面圖Fig.11 The Ez component amplitude plan

圖12 組合異常體模型示意圖Fig.12 Schematic diagram of composite abnormal body model

2.1.2 高阻異常體

圖9為中心測點A的Hz、Ez分量沿z軸的變化規(guī)律圖。分析圖9(a)能得出Hz分量不能分辨異常體上下界面，但在高阻體中心附近(z=62.5 m)存在極值。從圖9(b)可知，與低阻模型一樣，Ez分量在穿過高阻異常體上下界面時，也會發(fā)生突變異常，在界面處達到極大值。但與低阻體相比，垂直電場振幅曲線在穿越異常體界面時會有先降低后升高的小幅度變化，說明觀測總場變化與低阻情況時相反。

圖3(b)與圖9(b)為對比低、高異常體模型計算結(jié)果，在低阻異常體內(nèi)Ez幅值相對于背景圍巖明顯且穩(wěn)定，呈“U”變化，而在高阻異常體內(nèi)Ez幅值則呈“V”變化，形成這種差異的原因顯然與電磁場在高、低阻體內(nèi)的感應(yīng)特性有關(guān)。通過上述分析還可以發(fā)現(xiàn)低阻異常體引起的突變要比高阻異常體引起的突變更為嚴(yán)重，這可能是因為電磁場在低阻體內(nèi)的趨膚深度要小，因而使得異常體的幾何尺寸大小與電磁場的趨膚深度相當(dāng)，此時可能產(chǎn)生感應(yīng)效應(yīng)。

圖10為多個測點的Ez垂直響應(yīng)，由圖10可以發(fā)現(xiàn)，高阻體內(nèi)的測點特征相似，但相較于低阻體，Ez的幅值小，推測可能是由于高阻體阻礙了電流的流動，棱邊周邊離背景低阻體更近，電荷難以積累。當(dāng)測點位于異常體外時，Ez分量則不受影響。

圖11是Ez分量振幅平面圖。分析圖不難發(fā)現(xiàn)，同低阻模型一樣，Ez分量幅值在高阻體上下界面(z=40 m，z=85 m)存在高值，上界面處大于下界面，但對比低阻體，高值的數(shù)值較小。

2.2 組合異常體模型

電阻率為100 Ω·m的均勻半空間中含有一個電阻率為10 000 Ω·m的高阻異常體(上)和一個電阻率為1 Ω·m的低阻異常體(下)，兩異常體幾何尺寸均為30 m×30 m×30 m，高阻體頂面埋深30 m，低阻體頂面埋深95 m，兩者在z方向相距35 m，其余參數(shù)同單個異常體模型，如圖12所示。

圖13是組合模型中心測點A上垂直分量曲線圖。由圖13(a)可知，Hz分量在兩異常體中心位置達到極值，但不能分辨模型的上下界面。Ez分量曲線(圖13(b))大致吻合單個異常體時的響應(yīng)特征，在高、低阻異常體的上下界面處都有突變，到達極值，位置分別為z=30 m、z=60 m、z=95 m和z=125 m附近。但與單個高阻異常體時不同，組合模型中上方的高阻體上界面極值小于下界面，其垂直響應(yīng)受到模型下方低阻體的影響較大。說明了組合模型垂直響應(yīng)既保留了單個異常體模型時的特征，又反映了模型高、低阻異常體之間的相互作用。

圖13 中心點A處Hz、Ez幅值Fig.13 The Hz & Ez amplitude of center measurement point A(a)Hz分量；(b)Ez分量

圖14 多個測點Ez振幅Fig.14 The Ez amplitude of measurement points

圖15 Ez分量振幅平面圖Fig.15 The Ez component amplitude plan

對比圖14中的多個測點響應(yīng)特征，發(fā)現(xiàn)Ez分量對模型內(nèi)部測點的響應(yīng)異常明顯，對模型外測點基本上沒有影響，且在異常體的頂點位置反應(yīng)最明顯，同樣說明了頂點處電荷聚集多。

圖16 低阻體電荷分布示意圖Fig.16 The schematic diagram of the low resistivity charge distribution

分析Ez分量在y=0 m處的振幅平面圖(圖15)，發(fā)現(xiàn)Ez分量在組合模型上下四個界面處均呈現(xiàn)高值現(xiàn)象，與測點的測深曲線圖結(jié)果一致。且受低阻體影響，高阻體的下界面高值現(xiàn)象更強。

3 垂直分量變化物理機制分析

電磁波向下傳播過程中，當(dāng)遇到電性分界面時會有感應(yīng)電荷積累，其運動產(chǎn)生的二次電場會與一次場疊加[17]，如圖16所示。

圖16中的實線箭頭表示一次場，虛線箭頭表示二次場。由圖可知，當(dāng)異常體為良導(dǎo)時，在異常體上方二次電場與一次電場的方向相反，觀測總場是減小的。在異常體的兩端邊界附近，二次電場與一次電場是同向的，觀測總場是增大的。同理可知異常體為高阻時相反[18]。

根據(jù)電磁場理論[19]，設(shè)異常體界面兩端電場法向分量為En1和En2，則邊界條件為

En2-En1=ρs/ε0

σ1En1=σ2En2

(10)

式中:ρs為面電荷;ε0為真空介電常數(shù);σ為電導(dǎo)率。

由式(10)的第一式可知面電荷的產(chǎn)生會使面的法向分量產(chǎn)生間斷現(xiàn)象，而面電荷并不會產(chǎn)生切向分量的異常突變，因為面電荷產(chǎn)生的切向分量在界面上下并不產(chǎn)生方向的突變，而z分量則會產(chǎn)生方向的突變，比如在正電荷上方會產(chǎn)生向上的z分量，而在下方則會產(chǎn)生向下的z分量，因此面電荷對z分量的影響至關(guān)重要。由于地下介質(zhì)的橫向不均勻性，當(dāng)存在非水平面分布且互不等量的積累電荷時，就會有沿水平方向傳播的電場，使得總場的方向不再垂直向下，出現(xiàn)Ez異常響應(yīng)。

聯(lián)合上式可推導(dǎo)出電荷分布的表達式

(11)

根據(jù)公式(11)可知，地下垂直電場受積累面電荷影響，而異常體的物性特征如電導(dǎo)率等決定了面電荷的積累量[20]。

4 結(jié)論

通過對高頻情況下三維模型垂直響應(yīng)進行研究，分析了不同異常模型下垂直分量響應(yīng)的變化規(guī)律及特征。模型結(jié)果中的Ez分量對電性變化異常靈敏，且易于在鉆井中采集，因此或可以采用多道電場同時采集，將相鄰兩道相除實現(xiàn)歸一化，可以獲得相對穩(wěn)定的傳輸函數(shù)。通過對z分量的模擬得出以下結(jié)論：

1)當(dāng)電磁場穿越異常體上下界面時，垂直電場Ez會發(fā)生突變；垂直磁場Hz不會發(fā)生突變，但在異常體中心處出現(xiàn)極值。電磁場垂直分量Ez能夠很好地反映地下構(gòu)造的上下界面，在區(qū)分異常體上下界面時Ez分量響應(yīng)優(yōu)于Hz分量。

2)界面處Ez的突變主要是由界面上的面電荷導(dǎo)致，界面產(chǎn)生的面電荷主要聚集在異常體棱邊上，由式(10)可知，面電荷的產(chǎn)生只會對面上法向分量的間斷產(chǎn)生影響，并未產(chǎn)生切向分量的間斷。