考慮激電效應和電磁效應的復電阻率法 二維數值模擬研究

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

(1.核工業北京地質研究院 中核集團鈾資源勘查與評價技術重點實驗室，北京 100029；2.中國地質科學院礦產資源研究所國土資源部成礦作用與資源評價重點實驗室，北京 100037；3.中國地質大學(北京)地球物理與信息技術學院，北京 100083；4.核工業二四三大隊，內蒙古 赤峰 024000)

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

復電阻率法；有限單元；Cole-Cole模型；數值模擬；電磁效應

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

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

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

1　復電祖率法二維正演

1.1　二維數值模擬

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

圖1　地電模型示意圖

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

(1)

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

(2)

式中：^為波數域的值；ky為波數值。

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

(8)

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

(10)

(12)

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

(14)

1.2　有限單元法分析

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

(16)

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

(17)

k1(i,j)=

(18)

k2(i,j)=

(19)

k3(i,j)=

(20)

k4(i,j)=

(21)

s1(i)=

(23)

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

(24)

1.3　稀疏矩陣的存儲

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

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

表1　CSR存儲格式表

D17826953475ID14691012JD13524234424

1.4　Cole-Cole模型

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

(25)

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

1.5　正演基本流程

圖2為正演流程圖。

2　模型計算

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

2.1　層狀模型

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

圖2　正演流程圖

圖3　層狀模型

2.2　二維棱柱體計算結果對比

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

圖4　層狀模型計算結果對比圖

圖5　二維模型

Kerry Key的2DCSEM和本文開發的程序對該模型進行計算， 其計算結果如圖6所示。2DCSEM計算結果用圓圈表示，本文程序計算結果用點表示，電場分量Ex實部曲線、虛部曲線基本吻合，說明了本文開發的有限單元數值模擬算法正確。

3　結　論

本文基于有限單元法中四節點的等參單元完成了既考慮電磁效應又考慮激電效應的復電阻率法二維數值模擬。基于二次場滿足的麥克斯韋方程組，通過傅里葉變換得到麥克斯韋方程組波數域形式，結合有限單元法，推導出正演矩陣方程，為考慮激電效應，引入Cole-Cole模型，且應用第一類邊界條件，該方程組由超松弛迭代雙共軛梯度法來求解，將該波數域結果由反傅里葉變換變換到空間域。正演矩陣方程中的大型稀疏矩陣采用按行壓縮存儲，大大降低了內存占用量。最后，通過采用Kerry Key公開的程序與本文開發的程序對理論模型進行試算，驗證了本文正演程序正確、可靠。

圖6　Ex實虛部對比曲線

[1]楊進.環境與工程地球物理[M].北京:地質出版社,2011.

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

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

[4]徐自生,張麗,唐偉,等.復電阻率法(CR)在內蒙古那吉河鉛鋅礦探測中的應用[J].礦產勘查,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]許傳建,徐自生,楊志成,等.復電阻率(CR)法探測油氣藏的應用效果[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維復電阻率電磁場正反演研究[D].長春:吉林大學,2007.

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

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

[13]趙廣茂.帶地形的復電阻率2.5維電磁場正反演研究[D].長春:吉林大學,2009.

[14]尚曉.起伏地形條件下2.5維CR有限元數值模擬研究[D].北京:中國地質科學院,2012.

[15]尹成芳,柯式鎮,張雷潔.復電極型復電阻率掃頻系統響應數值模擬[J].測井技術,2014,38(3):273-278.

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

[17]范翠松,李桐林,嚴加永.2.5維復電阻率反演及其應用試驗[J].地球物理學報,2012,55(12):4044-4050.

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

[19]徐世浙.地球物理中的有限單元法[M].北京:科學出版社,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.