用有限單元法模擬各向異性介質(zhì)中二維電性異常體的大地電磁響應(yīng)*

熊 彬，羅天涯，李長(zhǎng)偉，王有學(xué)，張 智，覃小鋒，丁彥禮

(桂林理工大學(xué)地球科學(xué)學(xué)院，廣西桂林 541004）

各向異性對(duì)大地電磁探測(cè)數(shù)據(jù)的處理解釋備受國(guó)內(nèi)外地球物理學(xué)家的關(guān)注，徐世浙等［1］曾給出二維各向異性地電斷面大地電磁場(chǎng)的有限元解法并指出各向異性體的異常幅值約為各向同性體的5.3倍;Osella和Martinelli［2］研究了二維各向異性大地電磁響應(yīng)，他們發(fā)現(xiàn)TE模式對(duì)于二維各向異性地電結(jié)構(gòu)的響應(yīng)不靈敏，它不存在各向異性現(xiàn)象，而TM模式則與TE模式的情況相反;他們關(guān)于TE、TM模式的定義與文獻(xiàn) ［6］的定義不一致，本文沿用文獻(xiàn) ［6］的定義方式;Li［3］研究二維電性各向異性電磁響應(yīng)的有限元解法，假設(shè)走向沿x軸，y垂直于x軸，保持水平，z軸垂直向下，他研究了不同層面傾角下各向異性介質(zhì)視電阻率的變化特征，并得出ρxy不依賴于傾角，這說明ρxy不受各向異性的影響，而ρxy曲線在異常體處的極值點(diǎn)偏向反傾向一側(cè)，且偏移程度隨傾角增大而增大;沈金松等［4］對(duì)各向異性層狀介質(zhì)中視電阻率與磁場(chǎng)響應(yīng)進(jìn)行研究，他們認(rèn)為對(duì)于各向異性地層，沿地層走向的測(cè)線視電阻率最大，沿地層傾向的視電阻率最小，相比之下，磁電阻率對(duì)各向異性參數(shù)的分辨率較低。綜上所述，前人工作首先指出介質(zhì)的各向異性不可忽略，然后對(duì)TE、TM模式的電磁響應(yīng)進(jìn)行比較以及對(duì)傾斜角度等各向異性參數(shù)進(jìn)行研究。雖然前人對(duì)于各向異性介質(zhì)中的大地電磁場(chǎng)做了很多研究，但仍有一些各向異性參數(shù)值得討論?；诖?，本文首先對(duì)比TEz模式中各向同性與各向異性介質(zhì)的大地電磁場(chǎng)特征并與徐世浙［5］的結(jié)論對(duì)比，接著研究?jī)A角對(duì)各向異性的影響并與Li［3］的研究對(duì)比，再研究ρ‖和ρ┸等參數(shù)的變化對(duì)各向異性介質(zhì)中大地電磁場(chǎng)響應(yīng)的影響。

1 各向異性介質(zhì)中大地電磁場(chǎng)二維變分問題求解

假設(shè)地下電性結(jié)構(gòu)是二維的，取其走向?yàn)閦軸，x軸與z軸垂直，保持水平，y軸垂直向上，并設(shè)平行層面的電導(dǎo)率為δ‖，垂直層面的電導(dǎo)率為δ┸。根據(jù)場(chǎng)矢量的組成，可分成TMz和TEz兩種模式［6］。注意到各向異性介質(zhì)中TMz模式與電導(dǎo)率為δ‖的各向同性介質(zhì)中TMz模式滿足同樣的微分方程，文中直接針對(duì)TEz模式進(jìn)行討論。首先根據(jù)各向異性介質(zhì)中TEz模式下大地電磁場(chǎng)滿足的變分問題，利用有限單元法得到線性代數(shù)方程組，在解之前，代入第一類邊界條件，即可得到各節(jié)點(diǎn)的磁場(chǎng)強(qiáng)度，計(jì)算輔助場(chǎng)，從而可獲得地表的視電阻率［5］。文中采用的區(qū)域剖分方式為，通過引矩形單元的兩條對(duì)角線，將其細(xì)分為4個(gè)小三角形［7－8］，最終消去中間虛設(shè)的節(jié)點(diǎn)，這樣既便于模擬傾斜物性界面，又不會(huì)增加節(jié)點(diǎn)數(shù)。

2 數(shù)值模擬

2.1 向斜良導(dǎo)體模型

建立向斜模型［5］，取各向同性良導(dǎo)體的電阻率等于各向異性良導(dǎo)體的平均電阻率 (ρ1==1 Ω·m)，這是為了說明各向同性體與各向異性體中大地電磁場(chǎng)異常幅值的差異是由介質(zhì)本身電性結(jié)構(gòu)不同而造成的，而不是電阻率差異所致。從圖1中可知，即使各向同性良導(dǎo)體的電阻率等于各向異性良導(dǎo)體的平均電阻率，二者呈現(xiàn)出不同的異常幅值，例如在f=10－3Hz時(shí)，各向同性體異常幅值約為24倍，而各向異性體異常幅值約為84倍，這與徐世浙［5］得到的結(jié)果是一致的?？梢?，考慮各向異性對(duì)大地電磁數(shù)據(jù)處理與解釋很有必要。同時(shí)，可以看到隨著頻率的增大，各向異性體對(duì)大地電磁場(chǎng)的影響減弱，曲線逐步反映淺部各向同性圍巖的視電阻率。

圖1 各向同性與各向異性介質(zhì)大地電磁剖面比較Fig.1 Comparison of the magnetotelluric profile between isotropic and anisotropic media

2.2 不同傾角模型

如圖2所示各向異性低阻體，以逆時(shí)針為正，傾角α分別取0°，－45°，－90°，－135°得到如圖3所示TEz模式下視電阻率變化規(guī)律。從圖3中可以直觀地看到，當(dāng)角度為－45°和－135°時(shí)，曲線的谷值偏離模型中心，偏離方向與傾向相反。這與Li［3］得到的結(jié)果是一致的。

2.3 各向異性半空間模型

建立各向異性半空間 ρ‖=0.4 Ω·m，ρ┸=2.5 Ω·m，視電阻率變化規(guī)律如圖4所示。從圖4中可以清晰的看出，頻率從102Hz開始，視電阻率隨頻率增大而減小，而頻率在102～10－4Hz這個(gè)頻段視電阻率基本不改變，這說明在各向異性介質(zhì)中大地電磁勘探所采用的頻率應(yīng)小于102Hz這個(gè)數(shù)量級(jí)。

2.4 高阻異常體與低阻異常體視電阻率曲線對(duì)比

圖4 各向異性半空間中視電阻率隨頻率的變化Fig.4 Apparent resistivity varying with frequency in the anisotropic half space

采用2.1向斜良導(dǎo)體模型的幾何形態(tài)，但圍巖電阻率取ρ2=100 Ω·m，異常體分兩種情況討論:

1)各向異性高阻異常體(ρ1)‖=400 Ω·m，(ρ1)┷=2500 Ω·m，各向同性高阻異常體 ρ1=1000 Ω·m;計(jì)算結(jié)果如圖5所示。

2)各向異性低阻異常體，(ρ1)‖=4 Ω·m，(ρ1)┷=25 Ω·m，各向同性低阻異常體 ρ1=10 Ω·m，計(jì)算結(jié)果如圖6所示。這兩種下，均取各向同性體的電阻率等于各向異性體的平均電阻率，即ρ1=。從圖5和6可以得出3點(diǎn)結(jié)論:①高阻異常體的情況下，大地電磁場(chǎng)對(duì)各向異性體與各向同性的響應(yīng)區(qū)別很小，例如在f=10－3Hz時(shí)各向同性體的異常幅值約為1.81倍，各向異性體的異常幅值也僅僅約為1.77倍，從這里看出，各向同性體的異常幅值甚至比各向異性體的異常幅值還大;相比之下，異常體為低阻時(shí)，同樣在f=10－3Hz時(shí)，各向同性體的異常幅值約為8.5倍，各向異性體的異常幅值卻高達(dá)28倍。從這一對(duì)比可以得出，高阻體時(shí)，各向異性體異常幅值約為各向同性體的0.97倍，低阻體時(shí)，各向異性體異常幅值約為各向同性體的3.31倍，這說明各向異性低阻體對(duì)大地電磁場(chǎng)的影響大于各向異性高阻體的影響。②不討論各向異性的情況下，低阻體的異常幅值也大于高阻體的異常幅值。③高阻體時(shí)，視電阻率隨頻率變化幅度很小。

2.5 矩形模型

如圖7所示的矩形異常體，電阻率的變化分兩種情況進(jìn)行討論:①(ρ1)‖=1 Ω·m，(ρ1)┷取(10、102、103、104、105、106)Ω·m 等6個(gè)值，計(jì)算結(jié)果如圖8所示。② (ρ1)┷=1 Ω·m，(ρ1)‖取 (10、102、103、104、105、106)Ω·m 等6個(gè)值，計(jì)算結(jié)果如圖9所示。

圖7 矩形異常體Fig.7 Rectangle anomalous body

從圖8中可以得出，當(dāng)異常體平行于層面的電阻率(ρ1)‖保持不變，垂直于層面的電阻率(ρ1)┷不斷增加，乃至(ρ1)┷大于圍巖電阻率，斷面中一直呈現(xiàn)出低阻異常體;而從圖9中可以看出，當(dāng)(ρ1)┷保持不變，(ρ1)‖在低于圍巖電阻率范圍內(nèi)不斷增加，斷面中出現(xiàn)低阻異常體，但當(dāng)(ρ1)‖大于或者等于圍巖的電阻率時(shí)，異常體的視電阻率分布規(guī)律發(fā)生改變;(ρ1)‖大于圍巖電阻率時(shí)，斷面中出現(xiàn)高阻異常體，但(ρ1)‖等于圍巖電阻率時(shí)，無法分辨出各向異性異常體，只表現(xiàn)出圍巖的視電阻率。從圖8和圖9中可以得出，平行于層面的電阻率對(duì)各向異性介質(zhì)平均電阻率的影響遠(yuǎn)遠(yuǎn)大于垂直層面的電阻率的影響。

3 結(jié) 論

文中首先對(duì)比各向同性與各向異性介質(zhì)中大地電磁場(chǎng)的異常特征;其次，對(duì)各向異性介質(zhì)參數(shù)變化對(duì)各向異性介質(zhì)的影響進(jìn)行討論，從算例中得出以下結(jié)論:

1)介質(zhì)各向異性不可忽略。同一頻率下 (f=10－3Hz)，即使各向同性體電阻率與各向異性體平均電阻率相同，各向異性體視電阻率的異常幅值約為各向同性體的異常幅值的3.5倍。因此在大地電磁測(cè)深工作中，介質(zhì)的各向異性是不可忽視的因素。

2)層面傾角對(duì)視電阻率的影響。層面傾角的存在使得視電阻率曲線的極值點(diǎn)偏離模型中心，偏離方向與地質(zhì)體傾向相反。

3)頻率對(duì)視電阻率的影響。高頻對(duì)視電阻率影響比較明顯，在解釋大地電磁數(shù)據(jù)時(shí)應(yīng)注意考慮頻率的影響。

4)高、低阻體的異常特征。高阻體的各向異性特征不明顯，例如頻率為10－3Hz時(shí)各向異性體的異常幅值僅約為各向同性體的0.97倍;而低阻體情況下，各向異性異常幅值約為各向同性體的3.31倍;且高阻體隨頻率的變化不如低阻體的那么明顯。5)ρ‖與ρ┸對(duì)各向異性介質(zhì)的影響。平行于層面的電阻率對(duì)各向異性介質(zhì)平均電阻率的影響遠(yuǎn)遠(yuǎn)大于垂直層面的電阻率的影響，因而它對(duì)大地電磁場(chǎng)變化規(guī)律影響更大。

［1］徐世浙，趙生凱.二維各向異性地電斷面大地電磁場(chǎng)的有限元法解法［J］.地震學(xué)報(bào)，1985，7(1):80 －90.

［2］OSELLA A M，MARTINELLI P.Magnetotelluric response of anisotropic 2-D structures［J］.Geophys J Int，1993，115:819－828.

［3］LI Y.A finite-element algorithm for electromagnetic induction in two-dimensional anisotropic conductivity structures［J］.Geophys J Int，2002，148:389 －401.

［4］沈金松，郭乃川.各向異性層狀介質(zhì)中視電阻率與磁場(chǎng)響應(yīng)研究［J］.地球物理學(xué)報(bào)，2008，51(5):1608－1619.

［5］徐世浙.地球物理中的有限單元法［M］.北京:科學(xué)出版社，1994.

［6］TAFLOVE A，HAGNESS S C.Computational electrodynamics:The finite-difference time-domain method［M］.Boston，London:Artech House Publishers，2000:70 －72.

［7］羅延鐘，張桂清.電子計(jì)算機(jī)在電法勘探中的應(yīng)用［M］.武漢:武漢地質(zhì)學(xué)院出版社，1987.

［8］熊彬，羅延鐘.電導(dǎo)率分塊均勻的瞬變電磁2.5維有限元模擬［J］.地球物理學(xué)報(bào)，2006，49(2):590－597.