鐵軌干擾?？h地震臺(tái)地電阻率觀測(cè)的有限元分析

路中慧 李志濤 張達(dá) 郭少峰

（中國(guó)河南 456250 河南省浚縣地震臺(tái)）

0 引言

經(jīng)過50 多年的發(fā)展，地電阻率定點(diǎn)連續(xù)觀測(cè)已成為地震地球物理觀測(cè)的重要測(cè)項(xiàng)之一。在多次中強(qiáng)地震前，特別是在1976 年唐山7.8 級(jí)、2008 年汶川8.0 級(jí)地震前，均記錄到顯著的地電阻率異常。但是，目前地電阻率臺(tái)站受到的干擾愈來愈多，且受干擾情況較復(fù)雜，特別是地鐵運(yùn)營(yíng)、城際高鐵、高壓輸電線路、管道輸送等大型工程建設(shè)產(chǎn)生的干擾，嚴(yán)重影響了地電臺(tái)站觀測(cè)數(shù)據(jù)產(chǎn)出質(zhì)量；另外，還有水渠、大型鐵質(zhì)用具、金屬管道布設(shè)、鐵絲網(wǎng)等金屬管網(wǎng)及工業(yè)、農(nóng)業(yè)游散電流干擾，這些干擾均會(huì)引起地電阻率測(cè)值的變化。

地電阻率觀測(cè)數(shù)據(jù)用于地震預(yù)測(cè)時(shí)，首先要解決干擾問題。解滔等（2015）、劉素珍等（2017）基于影響系數(shù)理論，采用數(shù)值模擬方法對(duì)地震電阻率觀測(cè)受局部介質(zhì)電阻率變化的影響機(jī)理進(jìn)行了研究；解滔等（2015，2016）依據(jù)三維影響系數(shù)在地表的分布情況定性判斷干擾源對(duì)地電阻率觀測(cè)的影響形態(tài)，并結(jié)合不同測(cè)道測(cè)值排除或鎖定某些干擾源；解滔等（2015）、王同利等（2017）采用有限元數(shù)值計(jì)算方法建立三維模型，分析了地電阻率測(cè)區(qū)位于地表的金屬導(dǎo)線和局部異常體對(duì)觀測(cè)產(chǎn)生的干擾形態(tài)和幅度隨時(shí)間的變化特征。本文利用?？h地震臺(tái)測(cè)區(qū)鐵質(zhì)異常體相關(guān)信息及地質(zhì)剖面、巖層電性、電測(cè)深等資料，建立三維有限元模型，經(jīng)模型檢驗(yàn)后，結(jié)合實(shí)際觀測(cè)定量地分析局部鐵質(zhì)異常體對(duì)觀測(cè)的影響幅度和形態(tài)，以期為進(jìn)一步跟蹤?？h地震臺(tái)地電阻率觀測(cè)數(shù)據(jù)異常變化和地震前兆異常判定提供量化依據(jù)。

1 ?？h地震臺(tái)簡(jiǎn)介

?？h地震臺(tái)位于位于河南省浚縣白寺鄉(xiāng)中國(guó)人民解放軍73569 部隊(duì)農(nóng)場(chǎng)東部。該臺(tái)地處魯西隆起、太行山隆起及開封斷塊之間，構(gòu)造單元以塊狀結(jié)構(gòu)為主，凹凸相間。該區(qū)域的濮陽(yáng)斷塊為NE 向菱形塊體，其北界為磁縣—大名斷裂，南界為新鄉(xiāng)—商丘斷裂，西有太行山山前斷裂與太行山隆起接壤，東為聊蘭斷裂與魯西隆起分界。區(qū)域地層受青羊口斷裂、湯西斷裂及EW 向的安陽(yáng)南斷裂控制（圖1），浚縣地震臺(tái)距活動(dòng)構(gòu)造湯東斷裂僅25 km 左右。臺(tái)站海拔58 m，測(cè)區(qū)地貌平坦開闊，坡度高差極小，周圍主要為農(nóng)田，局部為林地，植被發(fā)育，環(huán)境干擾源較少。巖性表層360 m 以上為粘土層，以下為新近系泥灰?guī)r。地表大部分為全新統(tǒng)、上更新統(tǒng)地層所覆蓋，而且第四紀(jì)覆蓋層僅185 m，觀測(cè)環(huán)境較有利。

圖1 ?？h地震臺(tái)位置及周邊斷裂分布Fig.1 Location of Xunxian Seismic Station and distribution of peripheral faults

?？h地震臺(tái)地電阻率觀測(cè)系統(tǒng)于2013 年11 月30 日建成，2014 年初投入運(yùn)行，觀測(cè)數(shù)據(jù)質(zhì)量符合規(guī)范要求，年變化形態(tài)正常。臺(tái)站地電阻率觀測(cè)布極方式采用對(duì)稱四極法（圖2），共布設(shè)NS、EW 兩個(gè)測(cè)道，觀測(cè)室偏離布極區(qū)中心的距離為180 m。供電電極極距均為1 km，測(cè)量電極級(jí)距均為0.3 km，電極均采用900 mm×900 mm×8 mm 的鉛板電極，外線路采用地埋方式，電極埋設(shè)深度3.0 m，觀測(cè)儀器為ZD8M 型地電儀。

圖2 ?？h地震臺(tái)地電阻率觀測(cè)布極Fig.2 The distribution of observation poles of resistivity at Xunxian Seismic Station

在日常數(shù)據(jù)監(jiān)控中發(fā)現(xiàn)，自2016 年5 月13 日起EW 測(cè)道地電阻率數(shù)據(jù)出現(xiàn)明顯的下降階變，之后趨于穩(wěn)定。經(jīng)環(huán)境調(diào)查發(fā)現(xiàn)，2016 年5 月12 日，浚縣地震臺(tái)地電阻率EW 測(cè)線東供電極東端，新建1 條三角鐵軌，該鐵軌NS 向長(zhǎng)約300 m，寬約0.05 m，距東供電極約9.2 m。由于異常出現(xiàn)時(shí)間與鐵軌鋪設(shè)時(shí)間較接近，因此，EW 測(cè)道的數(shù)據(jù)變化是否與鐵軌有關(guān)，以及鐵軌對(duì)地電阻率觀測(cè)的影響需要進(jìn)一步驗(yàn)證。

2 有限元數(shù)值分析

2.1 穩(wěn)恒電流場(chǎng)有限元法

地電阻率定點(diǎn)觀測(cè)中采用對(duì)稱四極裝置，觀測(cè)過程中，在供電電極A、B 輸入直流電流，在測(cè)量電極M、N 測(cè)量電勢(shì)差。此問題可視為穩(wěn)恒電流場(chǎng)計(jì)算，滿足Maxwell 方程組和電荷守恒定律，因此穩(wěn)恒電流場(chǎng)問題可表示為如下泊松方程

式中，V為由電流源產(chǎn)生的電位；σ為介質(zhì)電導(dǎo)率；δ(x,y,z)為Dirac delta 函數(shù)。有限介質(zhì)空間的全部邊界為Г。其中，一部分邊界沒有電流流出（如地表），滿足Neumann 邊界條件，記為ГS；其余邊界記為ГV，滿足Dirichlet 邊界條件。因此，式（1）滿足如下邊界條件

應(yīng)用虛功原理可得到穩(wěn)恒電流場(chǎng)泊松方程的有限元弱解形式

式中，Ω為計(jì)算區(qū)域；φ為任意的虛位移函數(shù)，在滿足Dirichlet 邊界條件的邊界上，虛位移函數(shù)φ=0。

地電阻率觀測(cè)在地表自然地滿足Neumann 邊界條件，在水平方向和垂直方向 （深度）可視為無窮遠(yuǎn)邊界，可以施加Dirichlet 邊界條件（V=0），也可以施加Neumann 邊界條件（Coggon，1971）。但是，建立的模型在水平、垂直方向上的尺度不可能是無限的，對(duì)于一固定尺寸的模型，在供電極距AB 大于一定的值后，對(duì)無窮遠(yuǎn)邊界施加Diriehlet 邊界條件時(shí)計(jì)算得到的地電阻率值將小于實(shí)際值；而對(duì)無窮遠(yuǎn)邊界施加Neumann 邊界條件時(shí)計(jì)算得到的地電阻率值將大于實(shí)際值（Dey et al，1979；Li et al，2005）。對(duì)固定的供電極距AB，模型尺寸越大，邊界對(duì)計(jì)算結(jié)果的影響就越小。但是模型越大，計(jì)算量也就越大，因此需要合理地選擇模型水平方向的尺寸和最底層厚度。

模型經(jīng)單元離散化和施加電流源、邊界條件后可對(duì)單元節(jié)點(diǎn)上的自由度（電位）進(jìn)行數(shù)值求解，求解出電位分布后可以獲得測(cè)量電極間的電位差，進(jìn)而依據(jù)對(duì)稱四極裝置系數(shù)計(jì)算地電阻率（Huang et al，2010）。

2.2 模型建立

根據(jù)?？h地震臺(tái)測(cè)區(qū)內(nèi)電測(cè)深和電性剖面的測(cè)量資料，該臺(tái)測(cè)區(qū)場(chǎng)地電性結(jié)構(gòu)橫向分布較均勻，具有明顯的水平層狀結(jié)構(gòu)的特點(diǎn)。該區(qū)域第四系沉積層主要由粉土和泥灰?guī)r組成，粉土層電阻率為10—40 Ω·m，泥灰?guī)r層電阻率為9—57 Ω·m，總體電阻率較低。

依據(jù)?？h地震臺(tái)EW 向電測(cè)深曲線，反演得到測(cè)區(qū)電測(cè)深曲線（圖3），將其作為三維有限元模型的基礎(chǔ)數(shù)據(jù)，建立水平層狀模型并假定淺層3 層介質(zhì)電阻率均勻，參考地電阻率觀測(cè)資料得到測(cè)區(qū)地下電性結(jié)構(gòu)，具體數(shù)值見表1。一些研究者（解滔等，2013；王同利等，2017）通過對(duì)多個(gè)臺(tái)站地電阻率觀測(cè)資料進(jìn)行有限元數(shù)值分析認(rèn)為，模型水平尺寸越大，地層越厚，計(jì)算結(jié)果就越接近理論值，同時(shí)計(jì)算量也越大。當(dāng)模型水平尺寸D＞6 倍AB、模型厚度H＞2 倍AB 時(shí)，計(jì)算結(jié)果越接近理論值。為了保證計(jì)算結(jié)果與理論值間具有較好的一致性，建立一定水平尺寸、不同厚度的模型。經(jīng)過多次試驗(yàn)，最終確定模型尺寸為長(zhǎng)6 000 m，寬6 000 m，厚3 000 m，模型水平面中心始終位于臺(tái)站觀測(cè)裝置布極中心，此時(shí)邊界效應(yīng)對(duì)計(jì)算的影響已小于儀器的觀測(cè)精度，可滿足模型計(jì)算的要求。

圖3 ?？h地震臺(tái)EW 向電測(cè)深曲線Fig.3 Resistivity sounding curve of Xunxian Seismic Station in EW direction

表1 浚縣地震臺(tái)EW 向電測(cè)深曲線反演電性結(jié)構(gòu)Table 1 The resistivity structure inverted from sounding curve of Xunxian Seismic Station in EW direction

設(shè)測(cè)區(qū)內(nèi)鐵軌為長(zhǎng)方形金屬鐵板，長(zhǎng)為300 m，寬為0.05 m，計(jì)算鐵板的面積，在模型中使用導(dǎo)線單元對(duì)鐵板劃分網(wǎng)格，從里到外依次采用10×10、20×20、50×50 進(jìn)行有限元網(wǎng)格劃分，圖4 分別為鐵軌所在位置有限元模型水平剖面和鐵軌所在水平面中心布極區(qū)網(wǎng)格劃分剖面。

圖4 鐵軌模型網(wǎng)格劃分（a）鐵軌所在位置有限元模型水平剖面；（b）鐵軌所在水平面中心布極區(qū)網(wǎng)格劃分剖面Fig.4 The meshing of rail model

2.3 計(jì)算結(jié)果

表2 為鐵軌鋪設(shè)前后?？h地震臺(tái)地電阻率EW 向觀測(cè)值與計(jì)算值。由表2 可見，鐵軌鋪設(shè)前后地電阻率觀測(cè)值變化為-0.04 Ω·m，有限元模型計(jì)算變化值為-0.08 Ω·m，此結(jié)果與實(shí)際觀測(cè)值的變化趨勢(shì)基本相同，但變化量較實(shí)際觀測(cè)值大。鐵軌對(duì)觀測(cè)值的影響趨勢(shì)與實(shí)際觀測(cè)值的下降同步，這主要是由局部低阻體的影響所致，與解滔等（2016）低阻干擾源影響所表現(xiàn)出的電阻率下降的研究結(jié)果相符。所以，2016 年5 月?？h地震臺(tái)地電阻率EW 測(cè)向測(cè)值下降變化與鋪設(shè)的鐵軌有關(guān)。

表2 鐵軌鋪設(shè)前后地電阻率觀測(cè)值、計(jì)算值（單位：Ω·m）Table 2 The observed and calculated results of apparent resistivity before and after the laying of rail

3 結(jié)論

采用三維有限元模型對(duì)?？h地震臺(tái)地電阻率測(cè)區(qū)內(nèi)鐵軌對(duì)觀測(cè)數(shù)據(jù)的影響進(jìn)行了數(shù)值模擬分析，得到以下結(jié)論。

（1）三維有限元初始模型的建立較關(guān)鍵。在建立模型時(shí)發(fā)現(xiàn)，當(dāng)對(duì)鐵質(zhì)異常體分別采用導(dǎo)線單元和鐵板進(jìn)行網(wǎng)格劃分時(shí)，計(jì)算所得結(jié)果有所不同，鐵板模擬比導(dǎo)線單元模擬計(jì)算結(jié)果明顯偏低，更接近實(shí)際觀測(cè)值，但不同模型的計(jì)算結(jié)果變化形態(tài)和趨勢(shì)一致。計(jì)算表明，鐵質(zhì)異常體產(chǎn)生的干擾與其相對(duì)電極的位置（王同利等，2016）和自身長(zhǎng)度間也有關(guān)系，當(dāng)?shù)妥梵w位于供電極附近時(shí)，干擾形態(tài)為下降階變，并且距電極越近，干擾幅度越大；當(dāng)干擾源長(zhǎng)度不同時(shí)，產(chǎn)生干擾的幅度也不同。

（2）由有限元數(shù)值計(jì)算結(jié)果發(fā)現(xiàn)，?？h地震臺(tái)東供電極附近的鐵軌會(huì)導(dǎo)致地電阻率觀測(cè)值的下降，數(shù)值分析結(jié)果比實(shí)際觀測(cè)數(shù)據(jù)變化偏大，這可能是數(shù)值分析模型與測(cè)區(qū)地下實(shí)際電性結(jié)構(gòu)之間存在一定差異，而干擾幅值和形態(tài)分析結(jié)果與實(shí)際觀測(cè)結(jié)果趨勢(shì)一致，則主要是由局部低阻異常體的影響所致。因此認(rèn)為，?？h地震臺(tái)地電阻率EW 測(cè)向測(cè)值出現(xiàn)的下降階變是由環(huán)境干擾引起的。