基于積分方程法的任意形狀源多次激發三維電磁場反演

李建平，尚通曉, 關藝曉

(1.山東科技大學 山東省沉積成礦作用與沉積礦產重點實驗室，青島 266590；2. 國土資源部地裂縫地質災害重點實驗室，南京 210018 3. 江蘇省地質調查研究院，南京 210018)

0 引言

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

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

1 任意形狀源一次場的計算

2 積分方程基本理論

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

(1)

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

(2)

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

(3)

可求得到空間內任意一點的電場或磁場值。

3 反演方法概述

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

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

(4)

(5)

上式中實測場和理論場值都有x,y,z三個分量，且下標i=1、2、…、m表示每個觀測位置點或工作頻率點。

(6)

(7)

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

(8)

進一步推導得

(9)

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

(PTP+λ)·ΔX=S

(10)

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

假設共進行了N次測量，每次測量中測點個數均為m個，但場源和測點的位置各不相同。先按照前文所述方法求出各次測量時的雅克比矩陣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)求取模型參數的修改量ΔX，從而將多次測量數據應用到一次反演中。式(13)的雅克比矩陣元素為：

(14)

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

(15)

又由式(3)得：

(16)

(17)

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

(18)

(19)

(20)

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

(21)

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

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

5.1 模型1

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

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

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

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

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

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

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

5.2 模型2

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

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

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

6 結論

本次研究提出并實現了均勻大地中任意形狀源多個位置激發下的三維電阻率反演算法，通過模型試算得出如下結論：

1)積分方程只針對異常體單元反演，反演算法的效率較高，計算量相對較小，為三維電阻率反演應用奠定了堅實的基礎。

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

3)將激發點和接收點不同的多次觀測的電磁場數據用于一次反演中，有效解決了反演數據不足導致的參數不收斂問題。

參考文獻：

[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] 朱伯芳.有限單元法原理與應用[M].北京:水利出版社,1979.

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

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

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

[7] 林昌洪，譚捍東，舒晴，等.可控源音頻大地電磁三維共軛梯度反演研究[J].地球物理學報，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].地球物理學報，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] 劉斌，李術才，李樹忱，等.基于不等式約束的最小二乘法三維電阻率反演及其算法優化[J].地球物理學報，2012，55(1)：260-268.

[13] 宛新林，席道瑛，高爾根，等.用改進的光滑約束最小二乘正交分解法實現電阻率三維反演[J].地球物理學報，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.