李建平，尚通曉, 關(guān)藝曉

(1.山東科技大學(xué) 山東省沉積成礦作用與沉積礦產(chǎn)重點(diǎn)實(shí)驗(yàn)室，青島 266590；2. 國土資源部地裂縫地質(zhì)災(zāi)害重點(diǎn)實(shí)驗(yàn)室，南京 210018 3. 江蘇省地質(zhì)調(diào)查研究院，南京 210018)

0 引言

積分方程法是計(jì)算三維異常體電磁場的一種有效方法[1-2]，與有限元、有限差分等方法需要全空間進(jìn)行剖分不同[3～5]，該方法只在異常區(qū)域進(jìn)行剖分，計(jì)算量相對較小，尤其適合計(jì)算三維電磁場正反演。

在野外勘探中，常用的激發(fā)源有接地的電偶源和不接地的回線源，它們均為規(guī)則形狀，但受實(shí)際條件的限制，很多條件下激發(fā)源不能布設(shè)成規(guī)則形狀，此外受激發(fā)源強(qiáng)度和接收裝置的影響，每次觀測采集的數(shù)據(jù)個(gè)數(shù)都是一定的，但各次觀測中的激發(fā)點(diǎn)坐標(biāo)，接收點(diǎn)坐標(biāo)等參數(shù)卻不一樣，所以不能把不同參數(shù)的觀測數(shù)據(jù)直接用于反演。本次研究針對此問題，首先利用文獻(xiàn)[6]給出的方法，求出任意形狀源產(chǎn)生的一次場，再將激發(fā)點(diǎn)和接收點(diǎn)不同的多組電磁場數(shù)據(jù)統(tǒng)一考慮，實(shí)現(xiàn)了此種條件下三維電磁場的反演。

1 任意形狀源一次場的計(jì)算

2 積分方程基本理論

根據(jù)麥克斯韋方程組、積分方程理論及電場和磁場的張量格林函數(shù)，可得均勻大地中三維異常體的積分方程為

(1)

為方便計(jì)算，將三維異常體剖分成N個(gè)小立方單元，并假設(shè)電阻率在每個(gè)剖分單元內(nèi)是均勻分布的，且每個(gè)小單元內(nèi)的電磁場值近似等于其中心點(diǎn)的值，則各剖分單元內(nèi)電場可用式(2)表示。

(2)

解方程(2)求得異常體內(nèi)各剖分單元總電場E(rm)后，再利用式(1)的剖分形式

(3)

可求得到空間內(nèi)任意一點(diǎn)的電場或磁場值。

3 反演方法概述

目前三維電磁場常用的反演方法有共軛梯度法[7-11]、最小二乘法[12-15]、α中心方法[16]、神經(jīng)網(wǎng)絡(luò)法[17]等，本研究采用文獻(xiàn)[15]給出的方法，把靈敏度矩陣分為線性項(xiàng)和非線性的偏微分項(xiàng)，減少了計(jì)算量，提高了精度，再通過經(jīng)典的阻尼最小二乘法 ，用正演模擬數(shù)據(jù)擬合實(shí)測數(shù)據(jù)，并逐步修改地電模型參數(shù)值 ，最終達(dá)到最優(yōu)擬合，反演得到地下三維異常體電阻率的分布。

反演過程可簡單描述如下：

(4)

(5)

上式中實(shí)測場和理論場值都有x,y,z三個(gè)分量，且下標(biāo)i=1、2、…、m表示每個(gè)觀測位置點(diǎn)或工作頻率點(diǎn)。

(6)

(7)

此時(shí)擬合誤差被表示成為電導(dǎo)率模型修改量Δx1,Δx2,…,Δxn的多元函數(shù)，其取極小值的條件為：

(8)

進(jìn)一步推導(dǎo)得

(9)

這樣分別取j=1、2、…、n，可以推導(dǎo)出如下求解模型修改量的線性方程組:

(PTP+λ)·ΔX=S

(10)

4 源多次激發(fā)條件下雅克比矩陣和右端矢量的求取

假設(shè)共進(jìn)行了N次測量，每次測量中測點(diǎn)個(gè)數(shù)均為m個(gè)，但場源和測點(diǎn)的位置各不相同。先按照前文所述方法求出各次測量時(shí)的雅克比矩陣P1、P2、…、PN，它們均為m行N列的矩陣。

N次測量后的總雅克比矩陣可用各次測量的雅克比矩陣表示為:

(11)

式(11)為N·m行n列矩陣。

多次測量后的右端總矢量為T=(t1,t2…,tn)，其中元素為：

(12)

將Q、T代入式(10)中并替換PTP和S，即：

(QTQ+λ)·ΔX=T

(13)

利用式(13)代替式(10)求取模型參數(shù)的修改量ΔX，從而將多次測量數(shù)據(jù)應(yīng)用到一次反演中。式(13)的雅克比矩陣元素為：

(14)

將式(4)代入式(14)得：

(15)

又由式(3)得：

(16)

(17)

求解式(17)需要首先求得異常體各單元中心點(diǎn)處電場對各單元電導(dǎo)率的偏導(dǎo)數(shù)，然后再疊加求出地表總場對各單元電導(dǎo)率的偏導(dǎo)數(shù)。前文中的Fi、pin具有x、y、z三個(gè)分量，展開得：

(18)

(19)

(20)

仿照式(2)的求取方法，異常體各單元中心點(diǎn)處電場對各單元電導(dǎo)率偏導(dǎo)數(shù)的矩陣方程可表示為：

(21)

5 三維電阻率異常體反演算例

針對任意形狀源激發(fā)的三維異常體反演，設(shè)計(jì)了以下兩個(gè)模型算例，第一個(gè)算例相對簡單，電阻率在異常體內(nèi)是不變的，剖分的小塊數(shù)目較少；第二個(gè)算例相對復(fù)雜，異常體內(nèi)部的電阻率是變化的，剖分的小塊數(shù)目較多，而且兩個(gè)算例激發(fā)源的形狀不同，從而說明算法的有效性和通用性。

5.1 模型1

模型1如圖1所示，均勻大地中存在一個(gè)電阻率異常體，其中心點(diǎn)坐標(biāo)為(0，0，100)，大小為200 m×200 m×50 m ；激發(fā)源為一梯形大回線(圖2)，中心坐標(biāo)為(0，0，0)，回線內(nèi)部有兩條測線，共25個(gè)測點(diǎn)；激發(fā)源工作頻率選擇為0.01 Hz、0.1 Hz、1 Hz、10 Hz；圍巖電阻率為ρ0=100 Ω·m，異常體電阻率真值為ρ=10 Ω·m，將異常體剖分為4×4×1個(gè)小塊，設(shè)各小塊內(nèi)部電阻率分布均勻且等于其中心點(diǎn)電阻率值。

圖1 模型1三維示意圖Fig．1 The 3D sketch map of model 1

圖2 梯形大回線示意圖Fig．2 The sketch map of the trapezium loop

先利用一次激發(fā)下25個(gè)觀測點(diǎn)的數(shù)據(jù)進(jìn)行反演，反演模型的初始電阻率值為ρ=50 Ω·m，阻尼因子為“1”，阻尼因子縮放系數(shù)為10。利用本文算法反演計(jì)算各小塊的電阻率，經(jīng)過13次迭代后，反演結(jié)果如圖3所示，此時(shí)擬合差為4.32×10-6，大于給定的閥值1×10-7，且擬合差不再隨迭代次數(shù)的增加而減小，各小塊電阻率反演值與真值相差較大，反演失敗。

經(jīng)分析可知，反演失敗的原因在于參加反演的數(shù)據(jù)量不足。如果用兩次觀測的數(shù)據(jù)進(jìn)行反演，設(shè)激發(fā)源中心點(diǎn)分別為(0，0，0)和(300，0，0)，每次觀測均為25個(gè)測點(diǎn)，共50個(gè)測點(diǎn)，其余參數(shù)與前文相同。經(jīng)過9次迭代后，反演結(jié)果如圖4所示，此時(shí)擬合差為1.06×10-8，小于給定的閥值1×10-7，反演結(jié)束，且各小塊電阻率反演值收斂于真值，說明算法是正確的。

圖3 單次觀測異常體電阻率參數(shù)反演結(jié)果Fig．3 Resistivity inversion results of abnormal body in single measurement

圖4 兩次觀測異常體電阻率參數(shù)反演結(jié)果Fig．4 Resistivity inversion results of abnormal body in two measurements

5.2 模型2

模型2如圖5所示，均勻大地中存在一個(gè)電阻率異常體，其中心點(diǎn)坐標(biāo)為(0，0，200)，大小為300 m×300 m×150 m；激發(fā)源為一等邊三角形大回線，邊長1 000 m(圖6)，回線內(nèi)部有5條測線，42個(gè)測點(diǎn)，利用三次觀測的數(shù)據(jù)進(jìn)行反演，中心點(diǎn)分別為(-500，0，0)、(0，0，0)和(500，0，0)；激發(fā)源工作頻率選擇為0.01 Hz、0.1 Hz、1 Hz、10 Hz、100 Hz；圍巖電阻率為ρ0=100 Ω·m，異常體電阻率見下圖模型真值，將異常體剖分為 個(gè)小塊，設(shè)各小塊內(nèi)部電阻率分布均勻且等于其中心點(diǎn)電阻率值。反演中模型的初始電阻率值為ρ=300 Ω·m，阻尼因子為“1”，阻尼因子縮放系數(shù)為10，經(jīng)過11次迭代后，此時(shí)擬合差為9.21×10-8，小于給定的閥值1×10-7，反演結(jié)束，反演結(jié)果如圖7～圖9所示。

圖9 第三層電阻率的真值和反演值Fig．9 The true and invert resistivity of the third layer(a) 第三層電阻率真值;(b) 第三層電阻率反演值

圖7、圖8、圖9為模型2各層的反演結(jié)果。模型2在z方向上剖分為三層 ，圖7(b)為第一層剖分小塊在150 m處xoy面上的反演結(jié)果，圖8(b)為第二層剖分小塊在200 m處xoy面上的反演結(jié)果，圖9(b)為第三層剖分小塊在250m處xoy面上的反演結(jié)果。從圖中可看出，盡管電阻率在模型體內(nèi)是變化的，而且模型體剖分后的小塊數(shù)目也比較多，但把三次觀測的數(shù)據(jù)同時(shí)用于計(jì)算，反演的結(jié)果仍收斂于真值，說明本文算法不僅適用于異常體電阻率不變的簡單情況，也適用于異常體電阻率變化的復(fù)雜情況。

6 結(jié)論

本次研究提出并實(shí)現(xiàn)了均勻大地中任意形狀源多個(gè)位置激發(fā)下的三維電阻率反演算法，通過模型試算得出如下結(jié)論：

1)積分方程只針對異常體單元反演，反演算法的效率較高，計(jì)算量相對較小，為三維電阻率反演應(yīng)用奠定了堅(jiān)實(shí)的基礎(chǔ)。

2)反演中激發(fā)源可以為任意形狀，更加符合野外實(shí)際情況。

3)將激發(fā)點(diǎn)和接收點(diǎn)不同的多次觀測的電磁場數(shù)據(jù)用于一次反演中，有效解決了反演數(shù)據(jù)不足導(dǎo)致的參數(shù)不收斂問題。

參考文獻(xiàn)：

[1] HOHMANN G W.Three dimensional induced polarization and electromagnetic modeling [J].Geophysics 1975, 40(2):309-324.

[2] WANNAMAKER P E, HOHMANN G W， SANFILIPO W A. Electromagnetic modeling of three dimensional bodies in layered earths using integral equations[J].Geophysics 1984,49(1):60-74.

[3] 朱伯芳.有限單元法原理與應(yīng)用[M].北京:水利出版社,1979.

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

[5] 周熙襄. 電法勘探數(shù)值模擬技術(shù)[M]. 成都:四川科學(xué)技術(shù)出版社,1986.

[6] 李建平,李桐林,張亞東. 層狀介質(zhì)任意形狀回線源瞬變電磁場正反演研究[J].物探與化探，2012,36(2):256-259.

[7] 林昌洪，譚捍東，舒晴，等.可控源音頻大地電磁三維共軛梯度反演研究[J].地球物理學(xué)報(bào)，2012，55(11)：3829-3839.

[8] MACKIE R L,MADDEN T R.Three-dimensional magnetotelluric inversion using conjugate gradients[J]. Geophys,J.Int,1993,115:215-229.

[9] NEWMAN G A,ALUMBAUGH D L.Three-dimensional magnetotelluric inversion using non-linear conjugate gradients[J]. Geophys,J.Int,2000,140:410-424.

[10] 胡祖志，胡祥云，何展翔.大地電磁非線性共軛梯度擬三維反演[J].地球物理學(xué)報(bào)，2006，49(4)：1226-1234.

[11] LIN C H,TAN H D,TONG T.Three-dimension conjugate gradient inversion of magnetotelluric full information data[J]. Applied Geophysics,2011,8(1):1-10.

[12] 劉斌，李術(shù)才，李樹忱，等.基于不等式約束的最小二乘法三維電阻率反演及其算法優(yōu)化[J].地球物理學(xué)報(bào)，2012，55(1)：260-268.

[13] 宛新林，席道瑛，高爾根，等.用改進(jìn)的光滑約束最小二乘正交分解法實(shí)現(xiàn)電阻率三維反演[J].地球物理學(xué)報(bào)，2005，48(2):439-444.

[14] SASAKI Y.3-D resistivity inversion using the finite-element method [J].Geophysics,1994,59(11):1839-1848.

[15] EATON P A. 3D electromagnetic inversion using integral equation[J]. Geophysical prospecting, 1989, 37:407-426.

[16] PETRICK W R JR, SILL W R, WARD S H. Three-dimensional resistivity inversion using alpha centers[J]. Geophysics,1981,46(8):1148-1163.

[17] EL-QADY G,USHIJIMA K. Inversion of DC resistivity data using neural networks[J].Geophysical Prospecting,2001,49(4):417-430.