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

1 長(zhǎng)安大學(xué)地質(zhì)工程與測(cè)繪學(xué)院，西安 710054

2 山東大學(xué)巖土與結(jié)構(gòu)工程研究中心，濟(jì)南 250061

1 引 言

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

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

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

2 波場(chǎng)變換的基本原理

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

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

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

根據(jù)文獻(xiàn)［14，16，21，29］可知，由擴(kuò)散方程到波動(dòng)方程轉(zhuǎn)換的對(duì)應(yīng)關(guān)系表達(dá)式為

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

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

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

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

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

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

4 預(yù)條件正則化共軛梯度法實(shí)現(xiàn)波場(chǎng)反變換

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

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

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

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

圖2 不同離散方式虛擬時(shí)間分布示意圖（a）等面積離散；（b）對(duì)數(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）轉(zhuǎn)化為

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

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

設(shè)矩陣可分解為

其中

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

假設(shè)根據(jù)矩陣A（v）構(gòu)造的預(yù)條件子為M（v），構(gòu)造新的方程如式（8），矩陣M（v）－1A（v）接近單位陣，因此迭代很快收斂.具體計(jì)算過(guò)程如圖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 預(yù)條件正則化共軛梯度法計(jì)算框圖Fig.3 Computer block diagram of PRCG method

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

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

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

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

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

圖4 利用PRCG與RCG方法實(shí)現(xiàn)的波場(chǎng)反變換結(jié)果（a）利用PRCG方法獲得的波場(chǎng)反變換結(jié)果；（b）利用RCG方法求得的波場(chǎng)反變換結(jié)果.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）獲得的反變換結(jié)果，如圖4b所示.比較圖4a與圖4b可知，采用預(yù)條件處理后，正則化參數(shù)可以選取較?。╲＝0.0011285）就能達(dá)到很好的收斂效果；不采用預(yù)條件處理，正則化參數(shù)選擇較大時(shí)（v＝0.008），其收斂效果尚不如采用預(yù)條件技術(shù)的結(jié)果.由此可知預(yù)條件處理是必要的.

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

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

5 瞬變電磁虛擬波場(chǎng)分析

5.1 不同地電模型波場(chǎng)變換響應(yīng)

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

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

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

5.2 含噪信號(hào)的波場(chǎng)變換分析

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

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

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

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

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

圖21 測(cè)區(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.

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

6 結(jié) 論

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

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

圖22 測(cè)區(qū)全域波場(chǎng)變換結(jié)果（a）1320線全域波場(chǎng)變換結(jié)果；（b）1360線全域波場(chǎng)變換結(jié)果；（c）1400線全域波場(chǎng)變換結(jié)果.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ù)，在選定時(shí)段對(duì)于含有噪聲的瞬變電磁信號(hào)進(jìn)行波場(chǎng)變換，當(dāng)干擾在5%以下時(shí)波場(chǎng)變換結(jié)果與理論值吻合較好.但是當(dāng)干擾大于5%時(shí)，波場(chǎng)變換結(jié)果與理論值相差較大.因此瞬變電磁信號(hào)的前處理，即信號(hào)平滑與去噪是必要的.當(dāng)然我們也可以增大正則化參數(shù)，但這樣做勢(shì)必影響波場(chǎng)的分辨率.

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

（

）

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

［2］ 王家映.我國(guó)大地電磁測(cè)深研究新進(jìn)展.地球物理學(xué)報(bào)，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.大地電磁測(cè)深的擬地震解釋法.石油地球物理勘探，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］ 左海燕，王家映，方勝.關(guān)于大地電磁的平均速度問(wèn)題.石油物探，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］ 薛國(guó)強(qiáng)，李貅，宋建平等.回線源瞬變電磁成像的理論分析及數(shù)值計(jì)算.地球物理學(xué)報(bào)，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］ 郭文波.瞬變電磁擬地震解釋法研究［博士論文］.西安：長(zhǎng)安大學(xué)，2000.Guo W B.Transient electromagnetic pseudo－seismic interpretation method research［Ph.D.thesis］（in Chinese）.Xi′an：Chang′an University，2000.

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

［9］ 薛國(guó)強(qiáng)，李貅，底青云.瞬變電磁法正反演問(wèn)題研究進(jìn)展.地球物理學(xué)進(jìn)展，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］ 李貅.瞬變電磁測(cè)深的理論與應(yīng)用.西安：陜西科學(xué)技術(shù)出版社，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］ 陳本池，李金銘，周鳳桐.瞬變電磁場(chǎng)擬波動(dòng)方程偏移成像.石油地球物理勘探，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］ 陳本池，周鳳桐，李金銘.瞬變電磁場(chǎng)的波場(chǎng)變換研究.物探與化探，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］ 郭文波，李貅，薛國(guó)強(qiáng)等.瞬變電磁快速成像解釋系統(tǒng)研究.地球物理學(xué)報(bào)，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］ 李貅.瞬變電磁虛擬波場(chǎng)的三維曲面延拓成像研究［博士論文］.西安：西安交通大學(xué)，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］ 李貅，郭文波，胡建平.瞬變電磁測(cè)深快速擬地震解釋方法及應(yīng)用效果.西安工程學(xué)院學(xué)報(bào)，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］ 李貅，戚志鵬，薛國(guó)強(qiáng)等.瞬變電磁虛擬波場(chǎng)的三維曲面延拓成像.地球物理學(xué)報(bào)，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］ 李貅，薛國(guó)強(qiáng)，宋建平等.從瞬變電磁場(chǎng)到波場(chǎng)的優(yōu)化算法.地球物理學(xué)報(bào)，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］ 朱宏偉，李貅，張軍等.瞬變電磁法三維擬地震成像信息提取技術(shù).地球物理學(xué)進(jìn)展，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］ 劉銀愛.合成孔徑瞬變電磁偏移成像技術(shù)研究［博士論文］.西安：長(zhǎng)安大學(xué)地球探測(cè)與信息技術(shù)，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］ 劉繼軍.不適定問(wèn)題的正則化方法及應(yīng)用.北京：科學(xué)出版社，2005.Liu J J. Regularization Methods and Applications （in Chinese）.Beijing：Science Press，2005

［26］ 王彥飛.反演問(wèn)題的計(jì)算方法及其應(yīng)用.北京：高等教育出版社，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ù)值反演問(wèn)題.北京：高等教育出版社，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］ 沈梅芳.瞬變電磁的虛擬波場(chǎng)偏移成像研究［博士論文］.西安：長(zhǎng)安大學(xué)，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］ 周翠榮.改進(jìn)的正則化共軛梯度法［博士論文］.杭州：電子科技大學(xué)，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預(yù)處理技術(shù)在二維電磁特性TDFEM分析中的應(yīng)用.南京航空航天大學(xué)學(xué)報(bào)，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)定性.東南大學(xué)學(xué)報(bào)（自然科學(xué)版），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ù)方程組的迭代解法.北京：科學(xué)出版社，1991.Hu J G.Linear Algebraic Equations Iterative Method （in Chinese）.Beijing：Science Press，1991.

［37］ 趙俊華.改進(jìn)的SAOR預(yù)條件共軛梯度法［博士論文］.太原：太原理工大學(xué)，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］ 陳曉斌，趙國(guó)澤，湯吉等.大地電磁自適應(yīng)正則化反演算法.地球物理學(xué)報(bào)，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］ 戴亦軍，童孝忠，張連偉等.利用一維正則化反演進(jìn)行大地電磁測(cè)深數(shù)據(jù)擬二維反演解釋.物化探計(jì)算技術(shù)，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］ 張軍，李貅，趙瑩等.瞬變電磁虛擬波場(chǎng)高分辨成像技術(shù)研究.地球物理學(xué)進(jìn)展，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.