頻率域線源CSAET二維反演數(shù)值計算

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

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

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

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

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

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

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

1 正演計算

1.1 變分方程

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

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

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

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

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

1.2 有限單元法

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

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

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

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

1.3 正演算例驗證

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

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

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

圖3 電阻率100 Ω?m的均勻半空間解析解和有限單元法數(shù)值解對比圖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)的、非線性的，可以通過正則化方法對其進行求解[13]。如下式所表示：

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

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

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

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

2.2 靈敏度計算

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

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

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

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

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

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

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

式中：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è)計的模型在x和z軸組成的平面上，設(shè)計剖面的長度為6000 m，深度為3000 m，其中x軸上的坐標范圍為(-3000 m，3000 m)，z軸上的坐標范圍為(-3000 m，0 m)，圍巖電阻率設(shè)置為500 Ω·m；異常體設(shè)置 2個，位置在 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行組成，共計1352個剖分單元，選擇的背景場電阻率值為500 Ω·m，測量點為30個，點距200 m，選擇的頻點為分別為0.125、0.25、0.5、1、2、4、8、16、32、64 和 128 Hz，共計11個頻點。圖4(b)所示即為反演結(jié)果，從圖4(b)中可以看出：反演后最大的電阻率變?yōu)?60 Ω·m，最小電阻率為320 Ω·m，從反演結(jié)果的模型電阻率等值線斷面可以清楚的識別低阻區(qū)域和高阻區(qū)域。這從側(cè)面說明平滑度約束反演算法處理CSAET頻率域電磁數(shù)據(jù)時能很好地反演簡單的組合模型，但高阻異常與初始模型相比，電阻率變低；低阻異常與初始模型相比，電阻率變高。

圖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]。剖面線總長為17500 m，深度4000 m，其背景電阻率值為50 Ω·m，剖面深部鑲嵌兩個5 Ω·m的低阻異常，在地表附近設(shè)置了一個100 Ω·m的相對高阻異常和一個10 Ω·m的相對中等電阻率異常。此剖面設(shè)置為35個測點，點距為500 m，選擇的頻點仍分別為0.125、0.25、0.5、1、2、4、8、16、32、64 和 128 Hz共11個頻點。

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

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

4 結(jié)論

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

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

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

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

[1] 劉春明,佟鐵鋼,何繼善.多種電磁法在某金礦的野外勘探應(yīng)用[J].中國有色金屬學報,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] 閻 述,陳明生.線源頻率電磁測深二維正演(一)[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].中國有色金屬學報,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] 柳建新,童孝忠,郭榮文,李愛勇,楊 生.大地電磁測深法勘探—資料處理、反演與解釋[M].北京:科學出版社,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].北京:科學出版社,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.