三維高頻遠場平面波電磁勘探中的Ez研究

嚴家斌， 胡 濤， 趙 玄， 皇祥宇

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

0 引言

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

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

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

1 矢量有限元理論基礎

隨著計算機技術的不斷發展和學者對三維正反演算法的大量研究，對頻率域電磁法進行三維正演已不再困難。三維電磁模擬時，基于節點的標量有限元由于自由度賦予節點難以模擬突變和間斷，在電性界面突變處難以獲得高精度的場分量。基于矢量勢的三維有限元方法通過進一步求取微分來獲取電場和磁場，需要采取較細的網格獲得高精度Ez。因此為獲得高精度的垂直分量Ez、Hz，筆者采用矢量有限元來進行三維正演。

1.1 變分問題

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

▽×E=iωμH

▽×H=σE

(1)

式中:ω為角頻率；μ為磁導率；ε為介電常數；σ為電導率。

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

▽×▽×E-k2E=0

(2)

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

E=Exex+Eyey+Ezez

(3)

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

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

(4)

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

(5)

1.2 網格離散和矢量單元分析

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

(6)

(7)

(8)

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

KE=b

(9)

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

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

2 模型分析

2.1 單個異常體模型

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

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

因此利用Ez分量幅值在電性界面發生突變的性質我們可以判定異常體的上下分界面及厚度。對比Hz、Ez分量，不難發現Ez分量在分辨異常體的上下界面上優于Hz分量。

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

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

圖2 單個異常體模型網格剖分及示意圖Fig.2 Grid division and schematic diagram of a single anomalous body model(a)單個異常體模型網絡剖分圖；(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 歸一化處理結果圖Fig.8 Normalization processing results

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

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

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

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

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

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

2.2 組合異常體模型

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

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

圖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中的多個測點響應特征，發現Ez分量對模型內部測點的響應異常明顯，對模型外測點基本上沒有影響，且在異常體的頂點位置反應最明顯，同樣說明了頂點處電荷聚集多。

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

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

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

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

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

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

En2-En1=ρs/ε0

σ1En1=σ2En2

(10)

式中:ρs為面電荷;ε0為真空介電常數;σ為電導率。

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

聯合上式可推導出電荷分布的表達式

(11)

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

4 結論

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

1)當電磁場穿越異常體上下界面時，垂直電場Ez會發生突變；垂直磁場Hz不會發生突變，但在異常體中心處出現極值。電磁場垂直分量Ez能夠很好地反映地下構造的上下界面，在區分異常體上下界面時Ez分量響應優于Hz分量。

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