曹創(chuàng)華，童孝忠，柳建新

(1.中南大學(xué) 地球科學(xué)與信息物理學(xué)院，長(zhǎng)沙 410083；2.中南大學(xué) 有色資源與地質(zhì)災(zāi)害探查湖南省重點(diǎn)實(shí)驗(yàn)室，長(zhǎng)沙 410083)

頻率域人工源電磁測(cè)深法勘探自從推廣以來(lái)已經(jīng)得到廣泛的應(yīng)用[1]。近年來(lái)不論是正演模擬，還是反演計(jì)算，地球物理工作者都嘗試了很多，也相應(yīng)地取得了一定的成功[2]。盡管三維建模和反演算法現(xiàn)已取得很多進(jìn)步[3]，野外試驗(yàn)也逐漸展開，由于它們計(jì)算耗時(shí)和操作麻煩，成本較大，解釋技術(shù)還不成熟，還沒有廣泛應(yīng)用。

目前，二維反演分析仍然在實(shí)際應(yīng)用中扮演主要角色，但在國(guó)內(nèi)外常用的還是以二維的TM極化模式進(jìn)行反演，成功與失敗案例平分秋色。本文作者在科研和野外工作中發(fā)現(xiàn)，可以嘗試線源可控源音頻大地電流法(以下皆稱作CSAET法)剖面的二維反演。國(guó)內(nèi)幾乎沒有單位線源二維CSAET法反演進(jìn)行野外試驗(yàn)嘗試，數(shù)值模擬和反演計(jì)算更鮮有文獻(xiàn)論述。因此，嘗試優(yōu)化的線源二維CSAET反演算法顯得尤為必要。

眾所周知，有限單元法(FEM)和有限差分法(FDM)是對(duì)頻率域大地電磁法響應(yīng)進(jìn)行數(shù)值模擬的主要工具[4-5]，三維頻率域電磁法正演響應(yīng)也經(jīng)常用積分方程法[6]。對(duì)于一維層狀介質(zhì)，其一般都有解析解，較為簡(jiǎn)單；但當(dāng)面對(duì)二維復(fù)雜地形和復(fù)雜介質(zhì)模擬時(shí)，電磁場(chǎng)的正演響應(yīng)沒有解析解，且有限差分法在模擬起伏地形條件下的頻率域電磁響應(yīng)時(shí)沒有有限單元法精度高[7]，因此，本文作者利用有限單元法對(duì)線源二維CSAET法進(jìn)行數(shù)值模擬計(jì)算。

近年來(lái)，國(guó)內(nèi)外學(xué)者在電磁數(shù)據(jù)反演方面也取得了新的進(jìn)展，比如TORRES-VERDIN等[8]實(shí)現(xiàn)了起伏地形下的Born近似層析電磁數(shù)據(jù)反演成像問題；SAKASHITA等[9]利用這一原理從鉆孔電磁數(shù)據(jù)中提取了鉆孔的電阻率和磁化率特征曲線；ZHDANOV等[10]發(fā)展了一種把非線性方程組線性化的擬線性策略。本文作者在其基礎(chǔ)上，提出了一種平滑約束技術(shù)，通過(guò)反演模型參數(shù)轉(zhuǎn)換，來(lái)更好的解決這個(gè)問題；通過(guò)這種方法，反演迭代中的電導(dǎo)率更好的得到了限制，提高了反演的準(zhǔn)確度。

本研究中，利用平滑度約束反演方法反演了由線電流源產(chǎn)生的二維頻域電磁響應(yīng)數(shù)據(jù)。首先，采用有限單元法進(jìn)行了線源二維CSAET數(shù)值模擬得到了正演響應(yīng)數(shù)據(jù)。為了檢查算法的有效性，在文中把解析解和數(shù)值結(jié)果做了對(duì)比分析。其次，開發(fā)了一個(gè)平滑約束最小二乘反演方法，并利用伴隨方程法計(jì)算了反演過(guò)程中的靈敏度矩陣。最后，通過(guò)合成模型驗(yàn)證了本文算法的正確性。

1 正演計(jì)算

1.1 變分方程

本研究中假設(shè)的地電模型如圖1所示，y方向?yàn)樽呦蚍较?，沿y方向電阻率不變，電阻率的變化在xz平面，假設(shè)電磁波按照時(shí)間域諧波e-iωt變化，由電流源產(chǎn)生的電磁波在二維情況下，圖示測(cè)量模式為TE模式[11]，如下式所示：

式中：Εy表示y方向的電場(chǎng)強(qiáng)度，I為電流源的電流強(qiáng)度，δ為狄拉克函數(shù)，ω是角頻率，σ為模擬空間的電導(dǎo)率，μ(=4π×10-7H/m)為磁導(dǎo)率。

圖1 二維正演模型三維坐標(biāo)系示意圖Fig.1 3D Coordinate system for 2D forward modeling

將水平電偶極沿地層走向延伸到無(wú)限遠(yuǎn)，電偶極源變成了線源，那么對(duì)于二維地電模型來(lái)說(shuō)，場(chǎng)也就是二維，對(duì)于外邊界條件，采用無(wú)窮遠(yuǎn)邊界條件，即認(rèn)為電場(chǎng)在外邊界已經(jīng)衰減到非常小，接近于零，對(duì)結(jié)果的影響不大，可以忽略。式(2)變分式與邊值式(1)等價(jià)，如下：

式中：?為計(jì)算空間，即為圖1中的Site1～Site3之間的地表以下垂直于異常體的斷面。

1.2 有限單元法

為了解式(2)微分方程的解，本研究利用有限單元法，把計(jì)算域剖分成能夠進(jìn)行模擬的四邊形，如圖2所示。x方向上的網(wǎng)格在源附近劃分最密，在測(cè)區(qū)部位相對(duì)稀疏一些，其他部分劃分稀疏；y方向上以地表為界，空氣層劃分稀疏，地層下網(wǎng)格剖面在探測(cè)區(qū)域劃分較為密集，深部相對(duì)稀疏，本研究中采用矩形單元將整個(gè)區(qū)域剖分，矩形單元采用雙二次插值，在每個(gè)單元上取8個(gè)點(diǎn)(4個(gè)端點(diǎn)和4條邊的中點(diǎn))。每個(gè)矩形網(wǎng)格內(nèi)的電場(chǎng)場(chǎng)值由若干個(gè)節(jié)點(diǎn)值相加而成，用表示，如下式：

圖2 有限元剖分模型示意圖Fig.2 Schematic model of finite element subdivision

利用上述方法及式(4)插值公式對(duì)式(2)進(jìn)行了離散插值，由方程式(2)得到：

式中：P為與電流源相關(guān)的參數(shù)。假設(shè)每個(gè)單元的計(jì)算方程在滿足邊界條件下都可以得到解，然后將它們擴(kuò)展矩陣相加，得到y(tǒng)?=K ΕP，其中K為總體擴(kuò)展矩陣，Εye代表單元各節(jié)點(diǎn)y方向上的電位場(chǎng)值的列向量。

1.3 正演算例驗(yàn)證

為了驗(yàn)證有限單元法正演模擬的正確性，設(shè)置地電模型為均勻半空間，把利用有限單元法求出的數(shù)值解和均勻半空間條件下的解析解進(jìn)行了對(duì)比。均勻半空間的解析解的表達(dá)式已經(jīng)由WARD等[11]給出，如式(8)所示：

式中：K1為第二類1階的修正貝塞爾函數(shù)；k0和k1都為電磁波在模擬介質(zhì)中的波數(shù)。

假設(shè)均勻半空間的電阻率為100 Ω·m，電流源的方向沿y軸中心坐標(biāo)在地表的x=0 m處。在頻率100 Hz沿著x軸的有限單元數(shù)值解的場(chǎng)值計(jì)算結(jié)果如圖3所示，其結(jié)果用紅色小圓圈表示；細(xì)藍(lán)色實(shí)線表示同樣條件下由式(8)計(jì)算出的解析解，發(fā)現(xiàn)兩者對(duì)應(yīng)較好。通過(guò)此例說(shuō)明正演算法是有效正確的。

圖3 電阻率100 Ω?m的均勻半空間解析解和有限單元法數(shù)值解對(duì)比圖Fig.3 Comparison of finite element method response and analytical(solid line)solution on surface of 100 Ω ?m half-space(Symbols represent the FEM solution and solid line represents the analytical solution)

2 二維反演

2.1 平滑約束最小二乘反演

頻率域電磁數(shù)據(jù)的反演問題是病態(tài)的、非線性的，可以通過(guò)正則化方法對(duì)其進(jìn)行求解[13]。如下式所表示：

式中：()φ m為反演目標(biāo)函數(shù)，()S m為穩(wěn)定因子函數(shù)，β為正則化參數(shù)，正則化參數(shù)在求解 ()P m 最小值時(shí)起權(quán)衡作用。其中()φ m和()S m可用下兩式表示：

式中：F(m)為正演響應(yīng)函數(shù)，m表示模型空間向量，d表示數(shù)據(jù)空間向量，C表示模型參數(shù)的權(quán)重矩陣。

目前解此類方程有很多種數(shù)學(xué)方法[14-16]，此處利用平滑約束最小二乘反演方法求解正則化反演問題。使得式(9)線性化，如下式表示：

式中：Δd表示觀測(cè)的視參數(shù)和正演響應(yīng)數(shù)據(jù)的誤差向量，Δm表示每次反演計(jì)算迭代過(guò)程中的模型誤差向量，J為反演過(guò)程中的靈敏度矩陣或者由正演響應(yīng)函數(shù)F產(chǎn)生的雅可比矩陣，L為拉普拉斯(二階)平滑因子。

2.2 靈敏度計(jì)算

理論上，靈敏度的計(jì)算有如下3種方式：1)直接計(jì)算方法；2)靈敏度方程法；3)伴隨方程法[17]。這3種方法的計(jì)算時(shí)間分別與 N ×Mf、N×Mf和Mo×Mf的正演計(jì)算相關(guān)(N是正演參數(shù)的數(shù)量，Mf為參與計(jì)算的頻率的數(shù)量的數(shù)量，Mo為測(cè)量點(diǎn)的數(shù)量)。因?yàn)樵诜囱萦?jì)算過(guò)程中，模型參數(shù)的數(shù)量遠(yuǎn)遠(yuǎn)大于觀測(cè)數(shù)據(jù)的量，所以伴隨方程法計(jì)算靈敏度矩陣最為有效。

利用伴隨方程法計(jì)算了靈敏度矩陣，用此方法計(jì)算了電磁數(shù)據(jù)的反演問題。為了達(dá)到反演的目的，電導(dǎo)率被表示為基函數(shù)的有限線性組合τ( r)[18]：

式中：jψ表示基函數(shù)，jτ表示權(quán)系數(shù)；r為搜索半徑。這樣靈敏度矩陣的計(jì)算就可以由下式得出：

式中： ′Ε表示在電磁場(chǎng)空間Ω里面的輔助場(chǎng)場(chǎng)值。

2.3 正則化參數(shù)的選擇

正則化參數(shù)平衡著最小的數(shù)據(jù)誤差和模型粗糙度之間的關(guān)系，是評(píng)判反演過(guò)程穩(wěn)定性和分辨率的重要參數(shù)，使得它們兩者達(dá)到平衡。在此過(guò)程中利用自適應(yīng)正則化參數(shù)計(jì)算方法，計(jì)算式如下：

式中：k表示第k次反演擬合迭代的次數(shù)，正則化參數(shù)選擇與第k次迭代關(guān)系如表達(dá)式：

式中：0.01≤c≤0.5，βk為正則化因子最大值；＊β由式(15)自適應(yīng)算法的極小值得到。

圖4 模型一正演初始模型和反演結(jié)果Fig.4 Performance test of 2D smoothness-constrained inversion for simple model:(a)2D true model used for resistivity inversion showing both buried two rectangular bodies located starting at 500 m depth from surface;(b)Inversion result for frequency-domain electromagnetic data

3 算例

3.1 模型一

如圖4(a)所示，設(shè)計(jì)的模型在x和z軸組成的平面上，設(shè)計(jì)剖面的長(zhǎng)度為6000 m，深度為3000 m，其中x軸上的坐標(biāo)范圍為(-3000 m，3000 m)，z軸上的坐標(biāo)范圍為(-3000 m，0 m)，圍巖電阻率設(shè)置為500 Ω·m；異常體設(shè)置 2個(gè)，位置在 x軸(-1000 m，-500 m)、z軸(-1000 m，-500 m)的正方形異常設(shè)置為1000 Ω·m，位置在x軸(500 m，1000 m)、z軸(-1000 m，-500 m)的正方形異常設(shè)置為250 Ω·m。

反演的初始網(wǎng)格由52列和26行組成，共計(jì)1352個(gè)剖分單元，選擇的背景場(chǎng)電阻率值為500 Ω·m，測(cè)量點(diǎn)為30個(gè)，點(diǎn)距200 m，選擇的頻點(diǎn)為分別為0.125、0.25、0.5、1、2、4、8、16、32、64 和 128 Hz，共計(jì)11個(gè)頻點(diǎn)。圖4(b)所示即為反演結(jié)果，從圖4(b)中可以看出：反演后最大的電阻率變?yōu)?60 Ω·m，最小電阻率為320 Ω·m，從反演結(jié)果的模型電阻率等值線斷面可以清楚的識(shí)別低阻區(qū)域和高阻區(qū)域。這從側(cè)面說(shuō)明平滑度約束反演算法處理CSAET頻率域電磁數(shù)據(jù)時(shí)能很好地反演簡(jiǎn)單的組合模型，但高阻異常與初始模型相比，電阻率變低；低阻異常與初始模型相比，電阻率變高。

圖5 模型二正演初始模型(a)和反演結(jié)果(b)Fig.5 Performance test of 2D smoothness-constrained inversion for model as same as one tested by SASAKI[19]:(a)2D true model;(b)Inversion result for frequency-domain electromagnetic data

3.2 模型二

圖5(a)所示為SASAKI提出的電磁法經(jīng)典模型[19]。剖面線總長(zhǎng)為17500 m，深度4000 m，其背景電阻率值為50 Ω·m，剖面深部鑲嵌兩個(gè)5 Ω·m的低阻異常，在地表附近設(shè)置了一個(gè)100 Ω·m的相對(duì)高阻異常和一個(gè)10 Ω·m的相對(duì)中等電阻率異常。此剖面設(shè)置為35個(gè)測(cè)點(diǎn)，點(diǎn)距為500 m，選擇的頻點(diǎn)仍分別為0.125、0.25、0.5、1、2、4、8、16、32、64 和 128 Hz共11個(gè)頻點(diǎn)。

此模型的正演響應(yīng)仍然利用上文中提到的有限單元法，其初始網(wǎng)格由40列、和28行組成，共計(jì)1120個(gè)剖分單元，利用本研究中提到的方法進(jìn)行反演。圖5(b)所示為此模型的反演結(jié)果，可見反演結(jié)果分層較為明顯，電阻率范圍與模型一類似，具有高于背景電阻率的異常反演后電阻率變低而低于背景電阻率的異常反演后電阻率值變高的現(xiàn)象。

通過(guò)上述兩個(gè)模型的計(jì)算，表明二維平滑度約束方法可以有效地對(duì)CSAET頻率域電磁數(shù)據(jù)進(jìn)行反演。

4 結(jié)論

1)推導(dǎo)了垂直頻率域線源頻率域CSAET法二維電場(chǎng)微分方程，并推導(dǎo)了其變分方程，利用有限單元法進(jìn)行程序編制，實(shí)現(xiàn)了其正演響應(yīng)的求解。

2)編制了二維平滑度約束方法反演算法，實(shí)現(xiàn)了簡(jiǎn)單模型和SASAKI復(fù)雜模型的反演計(jì)算，對(duì)設(shè)置地層和異常體反演結(jié)果較好。

3)反演過(guò)程中采用自適應(yīng)正則化參數(shù)技術(shù)實(shí)現(xiàn)了自動(dòng)調(diào)節(jié)反演過(guò)程中的正則化因子。

4)從理論上分析了伴隨方程近似法計(jì)算靈敏度的耗時(shí)較小，并成功的進(jìn)行了反演計(jì)算。

[1] 劉春明,佟鐵鋼,何繼善.多種電磁法在某金礦的野外勘探應(yīng)用[J].中國(guó)有色金屬學(xué)報(bào),2013,23(9):2422-2429.LIU Chun-ming,TONG Tie-gang,HE Ji-shan.Exploration of various electromagnetic methods in some gold mine[J].The Chinese Journal of Nonferrous Metals,2013,23(9):2422-2429.

[2] UNSWORTH M J,TRAVIS B J,CHAVE A D.Electromagnetic induction by a finite electric dipole source over a 2-D earth[J].Geophysics,1993,58,198-214.

[3] COMMER M,NEWMAN G A.New advancesin three-dimensional controlled-source electromagnetic inversion[J].Geophysical Journal International,2008,172:513-535.

[4] COGGON J.Electromagnetic and electrical modeling by finite element method[J].Geophysics,1971,36:132-155.

[5] 閻 述,陳明生.線源頻率電磁測(cè)深二維正演(一)[J].煤田地質(zhì)與勘探,1999,27(5):60-62.YAN Shu,CHEN Ming-sheng.Two-dimensionalforward modeling of line source frequency electromagnetic sounding(Ⅰ)[J].Coal Geology&Exploration,1999,27(5):60-62.

[6] 李帝銓,謝 維,程黨性.E-Ex廣域電磁法三維數(shù)值模擬[J].中國(guó)有色金屬學(xué)報(bào),2013,23(9):2459-2470.LI Di-quan,XIE Wei,CHEN Dang-xing.Three-dimensional modeling for E-Ex wide field electromagnetic methods[J].The Chinese Journal of Nonferrous Metals,2013,23(9):2459-2470.

[7] 柳建新,童孝忠,郭榮文,李愛勇,楊 生.大地電磁測(cè)深法勘探—資料處理、反演與解釋[M].北京:科學(xué)出版社,2012:58-98.LIU Jian-xin,TONG Xiao-zhong,GUO Rong-wen,LI Ai-yong,YANG Sheng.The magnetotelluric sounding prospecting—Data processing,inversion and interpretation[M].Beijing:Science Press,2012:58-98.

[8] TORRES-VERDIN C,HABASHY T M.Rapid 2.5-D forward modeling and inversion via a new nonlinearscattering approximation[J].Radio Science,1994,29:1051-1079.

[9] SAKASHITA S,FUKUOKA K.An imaging technique for magnetic susceptibility and resistivity by electromagnetic tomography-Numerical experiments: Buturi-Tansa[J].Geophysics Exploration,1998,51:493-506.

[10]ZHDANOV M S,FANG S.Three-dimensional quasi-linear electromagnetic inversion[J].Radio Science,2012,31:741-754.

[11]WARDSH,HOHMANNGW,NABIGHIANMN.Electromagnetic Methods in applied geophysics[C]//Electromagnetic Theory for Geophysical Applications.Tulsa,USA:Society for Exploration in Geophysics,1988:131-312.

[12]徐世浙.地球物理中的有限單元法[M].北京:科學(xué)出版社,1994:5-28.XU Shi-zhe.Finite element method in geophysics[M].Beijing:Science Press,1994:5-28.

[13]TIKHONOV A N,ARSENIN V Y.Solution of Ill-Posed Problems[M].New York:Wiley,1977:3-22.

[14]ZHDANOV M S.Geophysical inverse theory and regularization problems[M].Amsterdam:Elsevier,2002:5-29.

[15]MITSUHATA Y,UCHIDA T,AMANO H.2.5-D inversion of frequency-domain electromagnetic data generated by a grounded-wire source[J].Geophysics,2002,74:1753-1768.

[16]CANDANSAYAR M E.Two-dimensionalinversion of magnetotelluric data with consecutive use of conjugate gradient and least-squares solution with singular value decomposition algorithms[J].Geophysical Prospecting,2008,56:141-157.

[17]MCGILLIVRAY P R,OLDENBURG D W.Methods for calculating Fréchetderivatives and sensitivities for the non-linear inverse problem:A comparative study[J].Geophys Prospect,1990,38:499-524.

[18]FARQUHARSON C G,OLDENBURG D W.Approximate sensitivities for the electromagnetic inverse problem[J].Geophysical Journal International,1996,126:235-252.

[19]SASAKI Y.Two-dimensional joint inversion of magnetotelluric and dipole-dipole resistivity data[J].Geophysics,1989,54:254-262.