考慮激電效應(yīng)和電磁效應(yīng)的復(fù)電阻率法 二維數(shù)值模擬研究

張志勇，劉固望，譚捍東，腰善叢，黃　笑，付振興

(1.核工業(yè)北京地質(zhì)研究院 中核集團(tuán)鈾資源勘查與評(píng)價(jià)技術(shù)重點(diǎn)實(shí)驗(yàn)室，北京 100029；2.中國(guó)地質(zhì)科學(xué)院礦產(chǎn)資源研究所國(guó)土資源部成礦作用與資源評(píng)價(jià)重點(diǎn)實(shí)驗(yàn)室，北京 100037；3.中國(guó)地質(zhì)大學(xué)(北京)地球物理與信息技術(shù)學(xué)院，北京 100083；4.核工業(yè)二四三大隊(duì)，內(nèi)蒙古 赤峰 024000)

復(fù)電祖率實(shí)測(cè)數(shù)據(jù)是包含激電效應(yīng)和電磁效應(yīng)的影響，而傳統(tǒng)的復(fù)電祖率法二維數(shù)值模擬只考慮激電效應(yīng)會(huì)影響后續(xù)資料處理結(jié)果的精度，因此本文開展同時(shí)考慮激電效應(yīng)和電磁效應(yīng)的復(fù)電祖率法二維正演研究。本文正演是基于二次場(chǎng)滿足的麥克斯韋方程組，通過傅里葉變換將其變換到波數(shù)域，結(jié)合有限單元法，引入Cole-Cole模型，形成正演矩陣方程，采用超松弛迭代雙共軛梯度法求解該方程組可得波數(shù)域結(jié)果，再由反傅里葉變換將其轉(zhuǎn)換至空間域，完成了正演研究。然后通過幾個(gè)模型算例的計(jì)算驗(yàn)證了該算法的正確性，也為復(fù)電阻率法二維反演奠定了基礎(chǔ)。

復(fù)電阻率法；有限單元；Cole-Cole模型；數(shù)值模擬；電磁效應(yīng)

復(fù)電阻率法的物質(zhì)基礎(chǔ)是巖石、礦石的電化學(xué)性質(zhì)差異，通過在地面觀測(cè)大地的頻譜，從而達(dá)到尋找電性異常體的目的[1]。相對(duì)其他電法分支方法，復(fù)電阻率法有以下優(yōu)點(diǎn)：①該方法能得到的電性參數(shù)較多，結(jié)合這些參數(shù)的對(duì)比解釋能提供更豐富的地電信息；②采用輕便的采集裝備，方便其在地形復(fù)雜地區(qū)開展工作。目前，復(fù)電阻率法已廣泛應(yīng)用在固體礦產(chǎn)勘探[2-4]、水文地質(zhì)[5]、油氣資源[6]、監(jiān)測(cè)環(huán)境[7]、工程地質(zhì)[8]等方面。因此，研究復(fù)電祖率法正演有一定的實(shí)際應(yīng)用價(jià)值。

國(guó)內(nèi)外許多學(xué)者做了關(guān)于復(fù)電祖率法數(shù)值模擬方面的工作。Loke等完成了復(fù)電阻率法二維正反演，其未考慮電磁效應(yīng)而僅僅考慮了激電效應(yīng)，其反演的電阻率和極化率初始模型是由一種近似反演方法獲取[9]。國(guó)內(nèi)，徐凱軍[10]、蔡軍濤等[11]、楊曉弘等[12]、趙廣茂[13]、尚曉[14]、尹成芳等[15]、張志勇等[16]都是基于泊松方程，完成了電阻率分塊均勻、網(wǎng)格剖分為三角形或四邊形的復(fù)電阻率法二維有限元數(shù)值模擬。范翠松等[17]、歐鷗[18]分別研究了復(fù)電阻率法二維正反演。

本文復(fù)電阻率法二維正演是基于二次場(chǎng)滿足的麥克斯韋方程組，經(jīng)過傅里葉變換將其變換到波數(shù)域，結(jié)合有限單元法，形成正演矩陣方程，同時(shí)考慮激電效應(yīng)和電磁效應(yīng)，將大地的實(shí)際電阻率由Cole-Cole模型定義得到的復(fù)電阻率替換，該方程采用超松弛迭代雙共軛梯度法求解，再由反傅里葉變換得空域結(jié)果，二維正演完成。

1　復(fù)電祖率法二維正演

1.1　二維數(shù)值模擬

本文采用的地電模型如圖1所示，其中y軸為構(gòu)造走向方向，假定其電性參數(shù)僅僅在x-z是變化的，而在走向方向是恒定無變化的。

圖1　地電模型示意圖

假設(shè)時(shí)間因子是eiωt，忽略位移電流的影響，二次場(chǎng)滿足的麥克斯韋方程組，見式(1)。

(1)

傅里葉變換，見式(2)。

(2)

式中：^為波數(shù)域的值；ky為波數(shù)值。

對(duì)式(1)沿著走向方向進(jìn)行傅里葉變換，可將其從空間域變換到波數(shù)域，整理后得式(3)～(8)。

(8)

求解式(5)和式(6)，得式(9)和式(10)。求解式(3)和式(8)，得式(11)和式(12)。

(10)

(12)

將式(9)和式(11)代入式(4)，可得式(13)。將式(10)和式(12)代入式(7)，可得式(14)。

(14)

1.2　有限單元法分析

結(jié)合有限單元法[19]中四節(jié)點(diǎn)的等參單元法推導(dǎo)，對(duì)式(13)和式(14)應(yīng)用伽里金法、格林公式，且正演中邊界條件為第一類邊界條件(也即是電磁場(chǎng)分量在邊界網(wǎng)格上的值都為零)，得到二次場(chǎng)滿足的電磁場(chǎng)有限單元法計(jì)算公式，見式(15)和式(16)。

(16)

經(jīng)過單元分析，將式(15)和式(16)寫為式(17)。

(17)

k1(i,j)=

(18)

k2(i,j)=

(19)

k3(i,j)=

(20)

k4(i,j)=

(21)

s1(i)=

(23)

為求解正演矩陣方程(式(17))，要先得到方程右端項(xiàng)S中的波數(shù)域一次電場(chǎng)分量。而其可以參考電磁法理論[20]進(jìn)行推導(dǎo)，得到式(24)[21]。

(24)

1.3　稀疏矩陣的存儲(chǔ)

式(17)中系數(shù)矩陣K是稀疏且對(duì)稱的，全部存儲(chǔ)K矩陣需要的內(nèi)存較大。由于矩陣K每行的非零元素個(gè)數(shù)有限，因此可采用按行壓縮(CSR)方式將其進(jìn)行存儲(chǔ)(表1)，這樣不僅減少內(nèi)存占用量，而且便于使用共軛梯度法來求解方程組。

假設(shè)K為N×N的稀疏且對(duì)稱矩陣，采用按行壓縮的方式存儲(chǔ)，矩陣K的存儲(chǔ)可用一維數(shù)組P、IP、JP進(jìn)行。K中每行非零元素值依次存放在數(shù)組P中；數(shù)組IP共有N+1個(gè)元素，每行首個(gè)非零元素在數(shù)組K中的位置存放在前N個(gè)數(shù)中，非零元素總數(shù)加1等于第N+1個(gè)數(shù)；K中每行非零元素所在列號(hào)存放在數(shù)組JP中。假設(shè)稀疏矩陣M如下所示。

表1　CSR存儲(chǔ)格式表

D17826953475ID14691012JD13524234424

1.4　Cole-Cole模型

為同時(shí)考慮激電效應(yīng)和電磁效應(yīng)，正演數(shù)值模擬時(shí)大地的實(shí)際電阻率由Cole-Cole模型定義的復(fù)電阻率來替換。1978年，Pelton等通過研究總結(jié)出可以用Cole-Cole模型來表征均勻巖石、礦石的復(fù)電阻率頻譜，見式(25)[22]。

(25)

式中：ρ(iω)為復(fù)電阻率；c為頻率相關(guān)系數(shù)；m為極化率；ρ0為零頻電阻率；τ為時(shí)間常數(shù)。

1.5　正演基本流程

圖2為正演流程圖。

2　模型計(jì)算

通過分別設(shè)計(jì)層狀模型、二維模型來驗(yàn)證本文二維數(shù)值模擬程序的正確性，將其計(jì)算結(jié)果分別與Kerry Key的1DCSEM程序[23]、Kerry Key的2DCSEM程序[24]計(jì)算結(jié)果進(jìn)行對(duì)比。下面兩個(gè)算例中，沿x軸水平方向放置發(fā)射源，且位于x=y=z=0處，從x=50 m到x=1 000 m放置25個(gè)接收點(diǎn)，頻率f為8 Hz、16Hz、32 Hz。

2.1　層狀模型

如圖3所示，在背景電阻率為100 Ω·m的地下介質(zhì)中存在電阻率為110 Ω·m，厚度為100 m的低阻層，低阻層頂界面埋深為120 m。將計(jì)算結(jié)果與Kerry Key的1DCSEM計(jì)算結(jié)果進(jìn)行對(duì)比驗(yàn)證。圖4中，1DCSEM計(jì)算結(jié)果用圓圈表示，而本文程序計(jì)算結(jié)果用點(diǎn)表示，電場(chǎng)分量Ex實(shí)部曲線、虛部曲線基本吻合，說明本文有限單元數(shù)值模擬算法的正演結(jié)果正確。

圖2　正演流程圖

圖3　層狀模型

2.2　二維棱柱體計(jì)算結(jié)果對(duì)比

如圖5所示，在背景電阻率為100 Ω·m的地下介質(zhì)中存在電阻率為10 Ω·m的二維棱柱體，其頂界面埋深為120 m， 長(zhǎng)200 m，厚100 m。分別采用

圖4　層狀模型計(jì)算結(jié)果對(duì)比圖

圖5　二維模型

Kerry Key的2DCSEM和本文開發(fā)的程序?qū)υ撃Ｐ瓦M(jìn)行計(jì)算， 其計(jì)算結(jié)果如圖6所示。2DCSEM計(jì)算結(jié)果用圓圈表示，本文程序計(jì)算結(jié)果用點(diǎn)表示，電場(chǎng)分量Ex實(shí)部曲線、虛部曲線基本吻合，說明了本文開發(fā)的有限單元數(shù)值模擬算法正確。

3　結(jié)　論

本文基于有限單元法中四節(jié)點(diǎn)的等參單元完成了既考慮電磁效應(yīng)又考慮激電效應(yīng)的復(fù)電阻率法二維數(shù)值模擬。基于二次場(chǎng)滿足的麥克斯韋方程組，通過傅里葉變換得到麥克斯韋方程組波數(shù)域形式，結(jié)合有限單元法，推導(dǎo)出正演矩陣方程，為考慮激電效應(yīng)，引入Cole-Cole模型，且應(yīng)用第一類邊界條件，該方程組由超松弛迭代雙共軛梯度法來求解，將該波數(shù)域結(jié)果由反傅里葉變換變換到空間域。正演矩陣方程中的大型稀疏矩陣采用按行壓縮存儲(chǔ)，大大降低了內(nèi)存占用量。最后，通過采用Kerry Key公開的程序與本文開發(fā)的程序?qū)碚撃Ｐ瓦M(jìn)行試算，驗(yàn)證了本文正演程序正確、可靠。

圖6　Ex實(shí)虛部對(duì)比曲線

[1]楊進(jìn).環(huán)境與工程地球物理[M].北京:地質(zhì)出版社,2011.

[2]楊振威,嚴(yán)加永,陳向斌.頻譜激電法在安徽沙溪斑巖銅礦中的應(yīng)用[J].地球物理學(xué)進(jìn)展,2013,28(4):2014-2023.

[3]曹蔚杰,單明良,高勇,等.復(fù)電阻率法(CR)在銅鉬礦勘查中的應(yīng)用效果[J].工程地球物理學(xué)報(bào),2014,11(2):203-207.

[4]徐自生,張麗,唐偉,等.復(fù)電阻率法(CR)在內(nèi)蒙古那吉河鉛鋅礦探測(cè)中的應(yīng)用[J].礦產(chǎn)勘查,2015(2):165-170.

[5]REVIL A,KARAOULIS M,JOHNSON T,et al.Review:Some low-frequency electrical methods for subsurface characterization and monitor in hydrogeology[J].Hydrogeology Journal,2012,20(4):617-658.

[6]許傳建,徐自生,楊志成,等.復(fù)電阻率(CR)法探測(cè)油氣藏的應(yīng)用效果[J].石油地球物理勘探,2004(zk):31-35.

[7]FLORES Orozco A,KEMNA A,OBERD?RSTER C,et al.Delineation of subsurface hydrocarbon contamination at a former hydrogenation plant using spectral induced polarization imaging[J].Journal of Contaminant Hydrology,2012,136-137:131-144.

[8]BREEDE K,KEMNA A,ESSER O,et al.Spectral induced polarization measurements on variably saturated sand-clay mixtures[J].Near surface geophysics,2012,10(6):479-489.

[9]LOKE M H,CHAMBERS J E，OGILVY R D.Inversion of 2D spectral induced polarization imaging data[J].Geophysical Prospecting,2006,54(3):287-301.

[10]徐凱軍.2.5維復(fù)電阻率電磁場(chǎng)正反演研究[D].長(zhǎng)春:吉林大學(xué),2007.

[11]蔡軍濤,阮百堯,趙國(guó)澤,等.復(fù)電阻率法二維有限元數(shù)值模擬[J].地球物理學(xué)報(bào),2007,50(6):1869-1876.

[12]楊曉弘,何繼善,童孝忠.頻率域激電有限元數(shù)值模擬[J].地球物理學(xué)進(jìn)展,2008,23(4):1186-1189.

[13]趙廣茂.帶地形的復(fù)電阻率2.5維電磁場(chǎng)正反演研究[D].長(zhǎng)春:吉林大學(xué),2009.

[14]尚曉.起伏地形條件下2.5維CR有限元數(shù)值模擬研究[D].北京:中國(guó)地質(zhì)科學(xué)院,2012.

[15]尹成芳,柯式鎮(zhèn),張雷潔.復(fù)電極型復(fù)電阻率掃頻系統(tǒng)響應(yīng)數(shù)值模擬[J].測(cè)井技術(shù),2014,38(3):273-278.

[16]張志勇,周峰.復(fù)電阻率2.5D有限單元法正演[J].地球物理學(xué)進(jìn)展,2014,29(5):2326-2331.

[17]范翠松,李桐林,嚴(yán)加永.2.5維復(fù)電阻率反演及其應(yīng)用試驗(yàn)[J].地球物理學(xué)報(bào),2012,55(12):4044-4050.

[18]歐鷗.起伏地形條件下2.5維復(fù)電祖率法數(shù)值模擬與反演成像研究[D].成都:成都理工大學(xué),2015.

[19]徐世浙.地球物理中的有限單元法[M].北京:科學(xué)出版社,1994.

[20]NABIGHIAN M.N.Electromagnetic Methods in Applied Geophysics[R].vol.1:Theory.Society of Exploration Geophysicists,1988.

[21]LU X Y.Inversion of Controlled-source Audio-Frequency Magnetotelluric data[D].Seattle：University of Washington,1999.

[22]PELTON W H,WARD S H,HALLOF P G,et al.Mineral discrimination and removal of inductive coupling with multifrequency IP[J].Geophysics,1978,43(3):588-609.

[23]KEY K.1D inversion of multicomponent,multifrequency marine CSEM data:Methododology and synthetic studies for resolving thin resistive layers[J].Geophysics,2009,74(2):F9-F20.

[24]KEY K，OVALL J.A parallel goal-oriented adaptive finite element method for 2.5-D electromagnetic modeling[J].Geophysical Journal International,2011,186(1):137-154.