戚志鵬，李 貅＊，吳 瓊，孫懷鳳，楊增林

1 長安大學地質(zhì)工程與測繪學院，西安 710054

2 山東大學巖土與結構工程研究中心，濟南 250061

1 引 言

瞬變電磁擬地震成像是實現(xiàn)瞬變電磁三維反演的有效手段，其中關鍵問題之一是實現(xiàn)由瞬變電磁擴散場到虛擬波場的轉換.只有進行了波場逆變換，才有可能實現(xiàn)基于波動方程的電磁法擬地震解釋.

所謂的電磁法擬地震解釋主要有三種，第一種是根據(jù)電磁波在導電介質(zhì)中和地震波在彈性介質(zhì)中傳播的相似性以及反射函數(shù)表達式之間的一致性把地震勘探中的一些基本方法應用于電磁測深資料解釋.主要成果有：Levy（1988）根據(jù)彈性波與大地電磁場之間的類比性，獲得了與反射地震類似的擬反射函數(shù)的脈沖響應時間斷面圖，將大地電磁資料用線性規(guī)劃技術對一維和簡單的二維地電斷面進行了成像計算，并推導了直流問題的擬脈沖響應函數(shù)［1］；王家映根據(jù)電磁波在導電介質(zhì)和地震波在彈性介質(zhì)中傳播的相似性及其反射函數(shù)表達式之間的一致性，在國內(nèi)率先提出了對大地電磁資料進行擬地震解釋的方法和思路［2－5］；郭文波、薛國強等基于中心回線裝置觀測成果與MT成果的一致性，借助大地電磁擬地震解釋，進行了TEM擬地震成像解釋［6－7］.第二種是借助偏移成像的廣義概念在擴散場條件下實現(xiàn)電磁偏移.主要成果有：前蘇聯(lián)學者Kjahos吸取了“偏移成像”的廣義概念，在電磁法中確立了正則偏移和解析延拓偏移兩種方法［8－10］；Zhdanov等［11－13］借鑒了地震勘探中的逆時偏移概念，對時間域的瞬變電磁場和頻率域的大地電磁場進行了逆時偏移成像的深入研究，提出了偏移電磁場的概念，并且建立在逆時偏移電磁場基礎上對二、三維反演問題也開展系列研究.第三種方法是基于Lavrent′ev（1980）的研究即瞬變電磁擴散場與虛擬波場之間存在著準確的數(shù)學變換關系，實現(xiàn)擴散場到波場的變換，在虛擬波場的條件下引入地震勘探處理和解釋的一些基本理論進行電磁成像［14］.主要成果有陳本池、李金銘等利用奇異值分解實現(xiàn)瞬變電磁波場變換，利用有限差分進行二維瞬變電磁擬地震成像解釋［15－16］；李貅等對瞬變電磁時間信號進行了分段處理，結合正則化算法分段實現(xiàn)了瞬變電磁波場變換，利用Kirchhoff積分實現(xiàn)了瞬變電磁三維曲面延拓成像［17－21］，并基于此開展了提高虛擬波場分辨率的系列研究［22－24］.但是，目前這方面的研究在算法和方法適用性方面仍存在不少問題.在對復雜介質(zhì)進行大地電磁擬地震解釋時遇到一些困難，如電磁波速度求取問題；基于擴散方程的逆時偏移成果僅限于單個或者多個孤立異常體的成像；利用擴散場到波場的變換關系可以實現(xiàn)復雜介質(zhì)的成像，但是由于波場反變換方程是第一類Fredholm積分，屬于典型的不適定問題［25－27］，李貅等早期研究為了降低問題的病態(tài)程度，采取了時間分段手段，有效降低了矩陣的病態(tài)程度，但這也引入了段與段的銜接問題［18，20，28］.

本文在李貅教授原有的研究基礎上進行瞬變電磁的波場變換算法，采用預條件正則化共軛梯度法實現(xiàn)全時域的波場變換算法，既有效避免了分段，又實現(xiàn)了由擴散場到波動場的變換.利用虛擬波場可以在波動方程下進行虛擬波場偏移成像.

2 波場變換的基本原理

早在1972年，Kunetz［29］在對大地電磁解釋和反演問題的研究中發(fā)現(xiàn)，電磁場滿足的擴散方程與波動方程之間具有相似關系.隨后，Lavrent′ev于1980 年，在 《Ill－posed problems of mathematical physics and analysis》一書中對擴散方程與波動方程之間的對應關系進行了深入研究，并給出了兩個方程之間準確的數(shù)學變換關系式［14］.Lee在1987年與1989年的研究成果中又進一步證明了這種數(shù)學變換適用于任意的場分量［30－31］.

在導電介質(zhì)中，忽略位移電流，瞬變電磁場滿足擴散方程.為了不失一般性，取f（x，y，z，t）為瞬變電磁場的電場或磁場分量函數(shù)，其滿足的擴散方程如式（1）所示：

其中μ為介質(zhì)磁導率，σ為電導率.引入函數(shù)U（x，y，z，τ），其滿足波動方程如式（2）所示：

根據(jù)文獻［14，16，21，29］可知，由擴散方程到波動方程轉換的對應關系表達式為

（3）式為第一類Fredholm型積分方程，由擴散場求波場是典型的不適定問題.將其離散后得到的線性代數(shù)方程組是病態(tài)的，且隨著階數(shù)的增加，矩陣條件數(shù)急劇增大，病態(tài)性更加嚴重，因此必須選用可靠的離散方式和穩(wěn)定的數(shù)值方法.本文采用預條件正則化共軛梯度法實現(xiàn)波場反變換的計算［32－33］.

3 波場變換的數(shù)值離散

計算式（4）的關鍵是如何選取一組τj（j＝1，2，…n）以及對應的一組wj（j＝1，2，…n）使其最好的滿足式（4）.對于虛擬時間與積分系數(shù)的選擇，文獻［18，20］通過給定的特殊積分，根據(jù)兩步最優(yōu)化選定一組最優(yōu)的虛擬時間和一組最優(yōu)的積分系數(shù).但是

將式（3）進行離散，寫成數(shù)值積分形式文中矩陣A的條件數(shù)較大，矩陣方程嚴重病態(tài).雖然如此，我們?nèi)钥梢越柚厥夥e分選取合適的離散方式，提取最優(yōu)的虛擬時間τ和最優(yōu)的積分系數(shù)w.

我們對比了一些常規(guī)的離散方式，如等間距離離散，對數(shù)等間距離離散，等面積離散和高度等間距離離散四種方式，得到不同離散方式的虛擬時間分布示意圖（圖2）.根據(jù)不同離散方式選取虛擬時間，采用數(shù)值積分方法進行積分，比較各數(shù)值積分精度和離散后各矩陣條件數(shù)，結果如表1所示.通過比較選取均方誤差最小，形成系數(shù)矩陣條件數(shù)最好的等間距離散.

表1 不同離散方式比較Table 1 The comparison of different discrete method

4 預條件正則化共軛梯度法實現(xiàn)波場反變換

考慮到文中方程組系數(shù)矩陣較大，單純的正則化共軛梯度法（RCG）或者預條件共軛梯度法（PCG），都不能很好地解決問題.因此，將兩種方法相結合形成預條件正則化共軛梯度法（PRCG）.由于在實際應用當中主要采集磁場的垂直分量，因此文中以磁場垂直分量為例進行分析.

將式（4）寫成矩陣形式為

其中A＝ ［aijwj］m×n，系數(shù)矩陣A包含積分系數(shù)wj，U＝［uj］n×1為虛擬子波，F(xiàn)＝h［］im×1為接收的瞬變場時間信號.

圖1 不同時間道的核函數(shù)展布圖（a）0.000001～0.0001s時間范圍；（b）0.1～10s時間范圍.Fig.1 Kernel functions at different times

圖2 不同離散方式虛擬時間分布示意圖（a）等面積離散；（b）對數(shù)等間距離散；（c）高度等間距離散；（d）線性等間距離散.Fig.2 The time display with different discrete way（a）Equal area discrete；（b）Logarithmic equal space discrete；（c）Equal altitude discrete；（d）Linear equal space discrete.

為了利用共軛梯度迭代，將式（5）轉化為

只要A是列滿秩矩陣，ATA就是對稱正定矩陣，因此可以利用共軛梯度法.但是ATA的條件數(shù)較A的條件數(shù)更大，使方程的病態(tài)更加嚴重.

為了進一步降低矩陣的條件數(shù)，改善方程的病態(tài)程度，在進行正則化共軛梯度之前，對系數(shù)矩陣進行預條件.對于預條件矩陣的構造采用超松弛預條件法，這是因為對稱超松弛預條件是一種較為有效的預條件方法，不僅預條件子容易求得，而且有效降低矩陣的條件數(shù)［34－37］.數(shù)學上已經(jīng)證明經(jīng)過超松弛預條件后，矩陣條件數(shù)降為原來的平方根［36］.

設矩陣可分解為

其中

K，Cl和Cu分別為S的對角元、下三角元和上三角元，ω為（0，2）內(nèi)的參數(shù).于是我們可以選擇預條件矩陣P如式（7）所示：

假設根據(jù)矩陣A（v）構造的預條件子為M（v），構造新的方程如式（8），矩陣M（v）－1A（v）接近單位陣，因此迭代很快收斂.具體計算過程如圖3所示，其中kmax為外層循環(huán)最大迭代次數(shù)，ε為正則化共軛梯度法迭代終止條件，lmax為內(nèi)層循環(huán)最大的迭代次數(shù)，ξ為內(nèi)層共軛梯度迭代終止條件，v為選定的正則化參數(shù).

其中，xk為第k次迭代的x值；x的初值x（0）選為單位向量；A（v）＝vI＋ATA.

圖3 預條件正則化共軛梯度法計算框圖Fig.3 Computer block diagram of PRCG method

正則化參數(shù)v的選擇非常重要，正則化參數(shù)v（δ）使U在近似性與穩(wěn)定性之間進行優(yōu)化選擇.Zhdanov等［38－39］提出正則化因子不斷遞減的自適應算法.并與L曲線法做了比較，認為采用自適應算法所得反演結果不遜于L曲線法所得結果，由于L曲線法需要通過多次反演計算來確定最佳的正則化因子，計算量將成倍增長，而自適應算法只需在每次迭代前確定參與當次反演的正則化因子，因此可以大大減少計算時間.在文獻［40］中，王彥飛依據(jù)偏差原理提出了重開始共軛梯度法（RSCG）計算最佳正則化因子，并與CGNR和CGER方法進行比較，認為RSCG方法更加穩(wěn)定［40］.正則化方法在大的電磁反演中已有廣泛的應用，如文獻［41－44］等均利用正則化方法實現(xiàn)了大地電磁反演，并給出了相應的正則化因子的確定方法.但是，瞬變電磁場較大地電磁場時間范圍與動態(tài)范圍均更廣，不適定性更加嚴重，不能簡單地引用大地電磁反演中使用的正則化方法.

經(jīng)分析，文獻［40］中所述的方法更適合本文情況.因此，文中依照文獻［40］重開始共軛梯度（RSCG）計算方法，根據(jù)部分先驗信息來選擇最優(yōu)的正則化參數(shù)值.正則化因子v為一漸變的量，其初始值v0為數(shù)據(jù)擬合泛函與穩(wěn)定泛函的比值，在前期大量模擬計算的基礎上基本可以確定正則化參數(shù)的初值，根據(jù)經(jīng)驗可以取v0等于0.00005為分段最大的正則化因子.在此后的迭代過程中，如果數(shù)據(jù)擬合殘差隨迭代次數(shù)逐漸變小，正則化因子可保持不變，否則按照（9）式進行選擇：

上式中ξ為經(jīng)驗系數(shù)，有ξ＞1，k表示預條件正則化共軛梯度（PRCG）的迭代過程中第k次迭代.

在完成積分離散以及正則化參數(shù)選擇后，對方程（4）式運用預條件正則化共軛梯度法進行求解計算.取波場理論值為U（x，y，z，τ）＝1，這時方程為一簡單積分，求得函數(shù)值H（x，y，z，t）＝表2所示為適用時間區(qū)間［0.000078，0.006280］的正則化參數(shù).利用（PRCG）法求得反變換結果如圖4a所示，最大誤差小于2%，可知文中所述方法滿足要求.

表2 正則化參數(shù)估計結果Table 2 Regularization parameter estimation result

圖4 利用PRCG與RCG方法實現(xiàn)的波場反變換結果（a）利用PRCG方法獲得的波場反變換結果；（b）利用RCG方法求得的波場反變換結果.Fig.4 Comparison wave－field inverse transform result with PRCG method and RCG method（a）The wave－field inverse transform result of synthetic wave records with PRCG method；（b）The wave－field inverse transform result of synthetic wave records with RCG method.

為了便于比較，給出采用正則化共軛梯度法（RCG）獲得的反變換結果，如圖4b所示.比較圖4a與圖4b可知，采用預條件處理后，正則化參數(shù)可以選取較小（v＝0.0011285）就能達到很好的收斂效果；不采用預條件處理，正則化參數(shù)選擇較大時（v＝0.008），其收斂效果尚不如采用預條件技術的結果.由此可知預條件處理是必要的.

眾所周知，正則化參數(shù)需要在穩(wěn)定性與分辨率之間平衡，正則化參數(shù)越大，分辨率越低.為了進一步說明問題，現(xiàn)給出同一模型兩種方法的波場變換結果，且兩方法選擇的正則化參數(shù)一致，為v＝0.003，如圖5所示，圖5a為采用PRCG方法獲得的波場變換結果，圖5b為采用RCG方法獲得波場變換結果.由圖可知模型具有兩個電性分界面，且兩個界面之間為一薄層，采用預條件正則化共軛梯度方法的結果可以分辨出薄層以及兩個電性界面，而正則化共軛梯度方法獲得的波場結果不能區(qū)別出兩個界面.由此可以說明，文中所述的PRCG反變換方法比RCG方法的分辨率高.

文獻［45］中，張軍、李貅（2009）等對不同時刻進行了分段，為了便于比較，截取相同時段計算結果與其進行對比，如圖6所示.對比可知，在相同時段內(nèi)，本文方法計算精度高于分段正則化算法精度.

5 瞬變電磁虛擬波場分析

5.1 不同地電模型波場變換響應

上文根據(jù)特殊子波對方法進行了驗證，說明方法的正確性.為了驗證方法對復雜介質(zhì)的適用性，利用文中波場變換方法分別對兩層模型和三層模型響應進行波場變換分析.瞬變電磁響應反映地下介質(zhì)的電阻率分布情況，當介質(zhì)為高電阻率時信號衰減較快，當介質(zhì)為低電阻率時信號衰減較慢.根據(jù)瞬變電磁衰減信號的特征，我們可以判斷層狀介質(zhì)的電阻率相對變化情況.以高斯脈沖為子波給出虛擬波場值，通過式（3）可以獲得與之對應的瞬變電磁衰減信號，并判斷地下介質(zhì)電阻率分布情況.根據(jù)瞬變電磁衰減信號可知，含有一個子波的為兩層模型，含有兩個子波的為三層模型，且子波為正表示下層電阻率較上一層電阻率低，子波為負表示下層電阻率較上一層電阻率高.

根據(jù)特殊子波U（x，y，z，τ）＝1可以確定正則化參數(shù)的合理范圍.選擇正則化參數(shù)后，對給定的地電模型進行模擬.對兩層模型的反變換結果如圖7、8所示，由圖可知，反變換所得波場與已知虛擬波場基本吻合.

三層模型的波場反變換結果如圖9、10、11、12所示，由圖可知，反變換結果與已知波場基本吻合，說明方法適合于復雜介質(zhì).圖中第二波峰較第一個波峰振幅有很大衰減，子波寬度略有增加.這與波在有耗介質(zhì)中的傳播規(guī)律相符，從而證明了算法的正確性.另一方面，模型的反演結果充分說明了瞬變電磁方法對淺部異常分辨率較高.隨著深度的增加，電磁波高頻部分損失嚴重，響應主要為低頻響應，因此對深部異常的分辨率降低.

5.2 含噪信號的波場變換分析

野外測量數(shù)據(jù)受各種因素影響難免會有干擾，因此有必要研究含噪信號的波場變換.文中以單脈沖響應的波場變換為例進行分析，加性噪聲分別為1%、2%、3%、5%.如圖13、14、15、16所示為不同信噪比的D型模型時域響應信號和與之對應的波場變換結果；圖17、18、19、20為不同信噪比的G型地電模型時間域響應和與之對應的波場變換結果.

由圖16b與20b可以明顯看出噪聲使得時間域響應曲線變得很不平滑.由波場反變換結果來看，當信噪比為小于5%時，反演值與期望值吻合得較好，當信噪比超過5%時，反演值波動劇烈，與預期結果差距較大.因此對于野外信號采集必須控制噪音為5%以下，并且對于局部不光滑數(shù)據(jù)不能直接進行處理，需要對數(shù)據(jù)進行平滑去噪后才能進行波場變換.因此本文方法完全適用于實際數(shù)據(jù)的處理.當然為提高波場變換的適用性也可以采取其他手段，如合成孔徑方法來提高波場的抗噪能力.

5.3 實測數(shù)據(jù)分析

在前文中已經(jīng)利用模型對方法進行了驗證，現(xiàn)對實測數(shù)據(jù)進行處理，進一步對算法進行驗證.以某煤田采空區(qū)實測數(shù)據(jù)為例進行分析.儀器選用GDP－32電法工作站，用中心回線方式進行測量，測線為1320、1360、1400，測點300—900.

測區(qū)視電阻率剖面如圖21所示.由圖可知，區(qū)域視電阻率變化較大，變化范圍從幾十到幾百歐姆米，總體表現(xiàn)為地表黃土覆蓋視電阻率高，基底視電阻率較高.在點號300—500深度120～350m左右范圍內(nèi)，有一低阻異常，結合工區(qū)資料，推斷為富水采空區(qū).根據(jù)視電阻率斷面圖可大概推斷采空區(qū)分布在點號（300，500）區(qū)間，由于深層電阻率變化緩慢，不容易估計采空區(qū)底界面，從視電阻率斷面圖可初步估計底界面在300～400m范圍內(nèi).

圖21 測區(qū)視電阻率斷面圖（a）1320線視電阻率深度斷面圖；（b）1360線視電阻率深度斷面圖；（c）1400線視電阻率深度斷面圖.Fig.21 Apperant resistivity section of the survey data（a）Apperant resistivity section of line 1320；（b）Apperant resistivity section of line 1360；（c）Apperant resistivity section of line 1400.

為了驗證方法效果，給出測線1320、1360、1400，點號為300到500間11個測點數(shù)據(jù)的全域波場變換結果.如圖22所示，其橫軸為點數(shù)，縱坐標為虛擬時間，單位為，在虛擬時間為 100處有一明顯的電性分界面，在虛擬時間為 400左右也有一電性分界面分布.根據(jù)圖21可知地下介質(zhì)類似于三層模型，與圖22波場變換結果相符.但是波場變換結果顯示的界面與視電阻率界面并不一致，這是由于圖22縱軸為虛擬時間，要想得到真實的地下介質(zhì)起伏形態(tài)還需要進行速度分析與延拓成像，虛擬波場速度是一個與地下介質(zhì)電阻率有關的量，經(jīng)過波場延拓后，得到深度剖面，其顯示結果應與視電阻率起伏形態(tài)相似.綜上所述，本文方法正確有效，能夠用于實際數(shù)據(jù)處理，并且對推斷多層采空具有明顯的優(yōu)勢.

6 結 論

文中主要研究由擴散方程到波動方程轉換的數(shù)值計算方法.首先，對前人推導的波場表達式進行數(shù)值離散，分析了多種離散方式的優(yōu)劣，選用效果較好的等間距剖分進行離散.然后，在前人工作基礎上，采用預條件正則化共軛梯度法，實現(xiàn)了全時域波場反變換.正則化共軛梯度法適合病態(tài)方程的求解，采用SSOR預條件子進行預條件，將系數(shù)矩陣條件數(shù)降為原條件數(shù)的平方根，有效改善方程的病態(tài)性，這樣選擇較小的正則化參數(shù)就可以得到較高的精度，有利于運用正則化共軛梯度實現(xiàn)方程的迭代計算.

實現(xiàn)波場反變換后，對算法進行了驗證，并對瞬變電磁虛擬波場的性質(zhì)進行了詳細討論.利用均勻模型（即子波為常數(shù)時），對反變換計算方法進行驗證，得到的子波與理論子波誤差小于2%，說明PRCG算法是有效的.對兩層（子波為單脈沖）和三層（子波為雙脈沖）模型進行反變換，得到的子波與理論子波相似，由此可以證明該方法適合于復雜介質(zhì).對兩層模型和三層模型進行了詳細的分析，可知虛擬波場對淺層目標分辨較好，隨著深度的增加，虛擬波場幅值迅速衰減，子波展寬嚴重，分辨率逐漸降低.

圖22 測區(qū)全域波場變換結果（a）1320線全域波場變換結果；（b）1360線全域波場變換結果；（c）1400線全域波場變換結果.Fig.22 The full－time wave－field transformation of the survey data（a）The full－time wave－field transformation of line 1320；（b）The full－time wave－field transformation of line 1360；（c）The full－time wave－field transformation of line 1400.

利用文中所述方法和選擇的正則化參數(shù)，在選定時段對于含有噪聲的瞬變電磁信號進行波場變換，當干擾在5%以下時波場變換結果與理論值吻合較好.但是當干擾大于5%時，波場變換結果與理論值相差較大.因此瞬變電磁信號的前處理，即信號平滑與去噪是必要的.當然我們也可以增大正則化參數(shù)，但這樣做勢必影響波場的分辨率.

經(jīng)過波場反變換后，虛擬波場不僅滿足波動方程，而且與地震子波類似，具有波的傳播特性.這為后續(xù)利用波動方程偏移成像方法提供了良好的基礎.

（

）

［1］ Levv S，Oldenburg D，Wang J.Subsurface imaging using magnetotelluric data.Geophysics，1988，53（1）：104－117.

［2］ 王家映.我國大地電磁測深研究新進展.地球物理學報，1997，40（S1）：206－216.Wang J Y.New development of magnetotelluric sounding in China.ActaGeophysicaSinica（ChineseJ.Geophys.）（in Chinese），1997，40（S1）：206－216.

［3］ 王家映，Oldenburg D，Levy S.大地電磁測深的擬地震解釋法.石油地球物理勘探，1985，20（1）：66－79.Wang J Y，Oldenburg D，Levy S.The magnetotelluric interpretation simuiating seismic method.OilGeophysical Prospecting（in Chinese），1985，20（1）：66－79.

［4］ 左海燕，王家映，方勝.關于大地電磁的平均速度問題.石油物探，1986，25（1）：84－97.Zuo H Y，Wang J Y，F(xiàn)ang S.On magnetotelluric average velocity.GeophysicalProsectingforPetroleum（in Chinese），1986，25（1）：84－97.

［5］ 王家映.大地電磁擬地震解釋法.北京：石油工業(yè)出版社，1995.Wang J Y.Magnetotelluric Sounding Interpretation Method（in Chinese）.Beijing：Petroleum Industry Press，1995.

［6］ 薛國強，李貅，宋建平等.回線源瞬變電磁成像的理論分析及數(shù)值計算.地球物理學報，2004，47（2）：338－343.Xue G Q，Li X，Song J P，et al.Theoretical analysis and numerical calculation of loop－source transient electromagnetic imaging.ChineseJ.Geophys.（in Chinese），2004，47（2）：338－343.

［7］ 郭文波.瞬變電磁擬地震解釋法研究［博士論文］.西安：長安大學，2000.Guo W B.Transient electromagnetic pseudo－seismic interpretation method research［Ph.D.thesis］（in Chinese）.Xi′an：Chang′an University，2000.

［8］ 牛之璉.時間域電磁法原理.長沙：中南大學出版社，2007.Niu Z L.Principle of Time Domain Electromagnetic Method（in Chinese）.Changsha：Central South University Press，2007.

［9］ 薛國強，李貅，底青云.瞬變電磁法正反演問題研究進展.地球物理學進展，2008，23（4）：1165－1172.Xue G Q，Li X，Di Q Y.Research progress in TEM forward modeling and inversion calculation.ProgressinGeophysics（in Chinese），2008，23（4）：1165－1172.

［10］ 李貅.瞬變電磁測深的理論與應用.西安：陜西科學技術出版社，2002.Li X.Transient Electromagnetic Sounding Theory and Application （in Chinese）.Xi′an.Shaanxi Science and Technology Press，2002.

［11］ Zhdanov M S，Dmitriev V I，F(xiàn)ang S，et al.Quasi－analytical approximations and series in electromagnetic modeling.Geophysics，2000，65（6）：1746－1757.

［12］ Zhdanov M S，Portniaguine O.Time－domain electromagnetic migration in the solution of inverse problems.Geophysics，1997，131（2）：293－309.

［13］ Zhdanov M S，Traynint P，Booker J R.Underground imaging by frequency－domain electromagnetic migration.Geophysics，1996，61（3）：666－682.

［14］ Lavrent′ev M M，Rornanov V G，Shishatskii S P.Ill－posed problems of mathematical physics and analysis（in Russian）.Providence RI：American Mathematical Society，1986.

［15］ 陳本池，李金銘，周鳳桐.瞬變電磁場擬波動方程偏移成像.石油地球物理勘探，1999，34（5）：546－554.Chen B C，Li J M，Zhou F T.Quasi wave equation migration of transient electromagnetic field.OGP（in Chinese），1999，34（5）：546－554.

［16］ 陳本池，周鳳桐，李金銘.瞬變電磁場的波場變換研究.物探與化探，1999，23（3）：195－201.Chen B C，Zhou F T，Li J M.The wavefield transformation study of transient electromagnetic field.GeophysicalandGeochemicalExploration（in Chinese），1999，23（3）：195－201.

［17］ 郭文波，李貅，薛國強等.瞬變電磁快速成像解釋系統(tǒng)研究.地球物理學報，2005，48（6）：187－192.Guo W B，Li X，Xue G Q，et al.A study of the interpretation system for TEM tomography.ChineseJ.Geophys.（in Chinese），2005，48（6）：187－192.

［18］ 李貅.瞬變電磁虛擬波場的三維曲面延拓成像研究［博士論文］.西安：西安交通大學，2005.Li X. The study about 3－D surface extension imaging technique in transient electromagnetic fictitious wave－field［Ph.D.thesis］ （in Chinese）.Xi′an：Xi′an Jiaotong University，2005.

［19］ 李貅，郭文波，胡建平.瞬變電磁測深快速擬地震解釋方法及應用效果.西安工程學院學報，2001，23（3）：42－45.LI X，Guo W B，Hu J P.The method and application effects of pseudo－seismic interpretation of TEM.JournalofXi′an EngineeringUniversity（in Chinese），2001，23（3）：42－45.

［20］ 李貅，戚志鵬，薛國強等.瞬變電磁虛擬波場的三維曲面延拓成像.地球物理學報，2010，53（12）：3005－3011.Li X，Qi Z P，Xue G Q，et al.Three dimensional curved surface continuation image based on tem pseudo wave－field.ChineseJ.Geophys.（in Chinese），2010，53（12）：3005－3011.

［21］ 李貅，薛國強，宋建平等.從瞬變電磁場到波場的優(yōu)化算法.地球物理學報，2005，48（5）：1185－1190.Li X，Xue G Q，Song J P，et al.An optimize method for transient electromagnetic field－wave field conversion.Chinese J.Geophys.（in Chinese），2005，48（5）：1185－1190.

［22］ 朱宏偉，李貅，張軍等.瞬變電磁法三維擬地震成像信息提取技術.地球物理學進展，2010，25（5）：1648－1656.Zhu H W，Li X，Zhang J，et al.Information collecting technology in 3－D pseudo－seismic imaging of transient electromagnetics.ProgressinGeophysics（in Chinese），2010，25（5）：1648－1656.

［23］ 劉銀愛.合成孔徑瞬變電磁偏移成像技術研究［博士論文］.西安：長安大學地球探測與信息技術，2010.Liu Y A.A research on TEM imaging method based on synthetic－aperture technology ［Ph.D.thesis］（in Chinese）.Xi′an：Chang′an University，2010.

［24］ Xue G Q，Yan Y J，Li X.Control of the waveform dispersion effect and applications in a TEM imaging technique for identifying underground objects.JournalofGeophysicsand Engineering，2011，8（2）：195－201，doi：10.1088／1742－2132／8／2／007.

［25］ 劉繼軍.不適定問題的正則化方法及應用.北京：科學出版社，2005.Liu J J. Regularization Methods and Applications （in Chinese）.Beijing：Science Press，2005

［26］ 王彥飛.反演問題的計算方法及其應用.北京：高等教育出版社，2007.Wang Y F.Computational Methods for Inverse Problems and Their Applications（in Chinese）.Beijing：Higher Education Press，2007.

［27］ 王彥飛，斯捷潘諾娃I E，提塔連科V N等.地球物理數(shù)值反演問題.北京：高等教育出版社，2011.Wang Y F，Stefan Nova I E，Titalianke V N，et al.Inverse Problems in Geophysics and Solution Methods（in Chinese）.Beijing：Advanced Education Press，2011.

［28］ 沈梅芳.瞬變電磁的虛擬波場偏移成像研究［博士論文］.西安：長安大學，2006.Shen M F.The fictitious wavefield migration imaging of transient electromagnetic method［Ph.D.thesis］（in Chinese）.Xi′an：Chang′an University，2006.

［29］ Kunetz G.Processing and interpretation of Magnetotelluric soundings.Geophysics，1972，37（6）：1005－1021.

［30］ Lee K H，Liu G，Morrison H F.A new approach to modeling the electromagnetic response of conductive media.Geophysics，1989，54（9）：1180－1192.

［31］ Lee S，McMechan G A，Aiken C L V.Phase－field imaging：The electromagnetic equivalent of seismic migration.Geophysics，1987，52（5）：678－693.

［32］ Bai Z，Zhang S.A regularized conjugate gradient method for symmetric positive definite system of linear equations.JournalofComputationalMathematics，2002，20（4）：437－448.

［33］ 周翠榮.改進的正則化共軛梯度法［博士論文］.杭州：電子科技大學，2010.Zhou C R.The improved regularized conjugate gradient method（in Chinese）［Ph.D.thesis］.Hangzhou：University of Electronic Science and Technology，2010.

［34］ 何小祥，劉梅林.SSOR預處理技術在二維電磁特性TDFEM分析中的應用.南京航空航天大學學報，2006，38（6）：670－673.He X X，Liu M L.Application of SSOR preconditioning technique in TDFEM for 2－D electromagnetic analysis.JournalofNanjingUniversityofAeronautics＆Astronautics（in Chinese），2006，38（6）：670－673.

［35］ 韋志輝，周榮富.SSOR方法的數(shù)值穩(wěn)定性.東南大學學報（自然科學版），1990，20（3）：108－113.Wei Z H，Zhou F R.Numerical stability of SSOR method.JournalofSoutheastUniversity（in Chinese），1990，20（3）：108－113.

［36］ 胡家贛.線性代數(shù)方程組的迭代解法.北京：科學出版社，1991.Hu J G.Linear Algebraic Equations Iterative Method （in Chinese）.Beijing：Science Press，1991.

［37］ 趙俊華.改進的SAOR預條件共軛梯度法［博士論文］.太原：太原理工大學，2005.Zhao J H.Modified SAOR preconditioned conjugate Gradient method ［Ph.D.thesis］（in Chinese）.Taiyuan：Taiyuan University of Technology，2005.

［38］ Zhdanov M S，Ellisz R，Mukherjee S.Three－dimensional regularized focusing inversion of gravity gradient tensor component data.Geophysics，2004，69（4）：925－937.

［39］ Zhdanov M S，Tolstaya E.A novel approach to the model appraisal and resolution analysis of regularized geophysical inversion.Geophysics，2006，71（6）：R79－R90.

［40］ Wang Y F.A restarted conjugate gradient method for illposed problems.ActaMathematicaeApplicataeSinica（EnglishSeries），2003，19（1）：31－40.

［41］ 陳曉斌，趙國澤，湯吉等.大地電磁自適應正則化反演算法.地球物理學報，2005，48（4）：937－946.Chen X B，Zhao G Z，Tang J，et al.An adaptive regularized inversion algorithm for magnetotelluric data.ChineseJ.Geophys.（in Chinese），2005，48（4）：937－946.

［42］ Tong X Z，Liu J X，Xu L H.Damped gauss－newton optimization algorithm for tow－dimensional magnetotelluric regularization inversion.ICIEC，2009，12

［43］ 劉小軍，王家林，吳健生.二維大地電磁正則化共軛梯度法反演算法.上海地質(zhì)，2007，（1）：71－74.Liu X J，Wang J L，Wu J S.Inversion algorithm of 2－D magnetotelluric data using regularized conjugate gradient method.ShanghaiGeological（in Chinese），2007，（1）：71－74.

［44］ 戴亦軍，童孝忠，張連偉等.利用一維正則化反演進行大地電磁測深數(shù)據(jù)擬二維反演解釋.物化探計算技術，2012，34（1）：33－38.Dai Y J，Tong X Z，Zhang L W，et al.Pseudo－2Dinversion interpretation for magnetotelluric data using 1Dregularization inversion method.ComputingTechniquesforGeophysical andGeochemicalExploration（in Chinese），2012，34（1）：33－38.

［45］ 張軍，李貅，趙瑩等.瞬變電磁虛擬波場高分辨成像技術研究.地球物理學進展，2011，26（3）：1077－1084.Zhang J，Li X，Zhao Y，et al.A technology research of highresolation imaging for the transient electromagnetic pseudo wave field.ProgressinGeophysics（in Chinese），2011，26（3）：1077－1084.