姚長利，李宏偉，鄭元滿，孟小紅，張聿文

1地下信息探測(cè)技術(shù)與儀器教育部重點(diǎn)實(shí)驗(yàn)室和地質(zhì)過程與礦產(chǎn)資源國家重點(diǎn)實(shí)驗(yàn)室，北京 100083

2中國地質(zhì)大學(xué)（北京）地球物理與信息技術(shù)學(xué)院，北京 100083

重磁位場(chǎng)轉(zhuǎn)換計(jì)算中迭代法的綜合分析與研究

姚長利，李宏偉，鄭元滿，孟小紅，張聿文

1地下信息探測(cè)技術(shù)與儀器教育部重點(diǎn)實(shí)驗(yàn)室和地質(zhì)過程與礦產(chǎn)資源國家重點(diǎn)實(shí)驗(yàn)室，北京 100083

2中國地質(zhì)大學(xué)（北京）地球物理與信息技術(shù)學(xué)院，北京 100083

處理轉(zhuǎn)換計(jì)算在重磁資料解釋中發(fā)揮著重要的作用，但一些計(jì)算如向下延拓、化極等有時(shí)是很不穩(wěn)定的，在頻率域中則表現(xiàn)為其轉(zhuǎn)換因子具有明顯的放大作用，所以其FFT理論計(jì)算結(jié)果是不穩(wěn)定的.因此，很多研究工作都是圍繞增加計(jì)算的穩(wěn)定性、提高計(jì)算效果進(jìn)行的，其中迭代法是近來在研究中受到普遍重視的方法技術(shù)，并取得了較好的成果.但也存在對(duì)迭代法研究還不夠深入，對(duì)其存在的缺點(diǎn)認(rèn)識(shí)不夠充分、客觀等問題，例如，迭代法進(jìn)行延拓及化極等計(jì)算時(shí)，對(duì)一些具體應(yīng)用雖能在一定程度上獲得較好的計(jì)算結(jié)果，但卻存在計(jì)算結(jié)果并不會(huì)隨著迭代次數(shù)的增加而得到持續(xù)改善的問題，對(duì)于原本不穩(wěn)定的計(jì)算，迭代法在迭代次數(shù)比較大時(shí)，所得的結(jié)果依然是不穩(wěn)定的.為此，本文在對(duì)迭代法進(jìn)行分析研究的基礎(chǔ)上，進(jìn)一步推導(dǎo)了迭代法的通式，并分析了對(duì)迭代法收斂性影響的各種因素.分析結(jié)果表明：迭代法收斂到FFT理論直接計(jì)算結(jié)果的決定因素是計(jì)算過程中如何選擇原始數(shù)據(jù)到目標(biāo)數(shù)據(jù)的映射函數(shù)；在選擇了合適的映射函數(shù)的情況下，迭代次數(shù)不僅僅是決定計(jì)算成本，而是決定結(jié)果好壞的關(guān)鍵因素；增加迭代次數(shù)雖然能夠使計(jì)算收斂到FFT直接計(jì)算理論結(jié)果，但如果該理論結(jié)果本身就是不穩(wěn)定的，則迭代法計(jì)算如果收斂，也是收斂到一個(gè)不穩(wěn)定的結(jié)果.所以針對(duì)位場(chǎng)處理轉(zhuǎn)換中一些不穩(wěn)定計(jì)算采用迭代法，并沒有從根本上解決計(jì)算的不穩(wěn)定性問題.

位場(chǎng)轉(zhuǎn)換，迭代法，F(xiàn)FT直接轉(zhuǎn)換，收斂性，收斂條件

1 引 言

重磁資料的分析解釋通常離不開選擇多種處理轉(zhuǎn)換方法進(jìn)行計(jì)算與分析對(duì)比的過程，即處理轉(zhuǎn)換在重磁資料解釋中發(fā)揮著重要的作用.但由于一些轉(zhuǎn)換計(jì)算如向下延拓、化極等有時(shí)是很不穩(wěn)定的，所以，如何得到穩(wěn)定有效的轉(zhuǎn)換計(jì)算結(jié)果就成為重磁位場(chǎng)的重要研究?jī)?nèi)容之一.

重磁位場(chǎng)的向下延拓是典型的不穩(wěn)定計(jì)算，其本質(zhì)是數(shù)學(xué)物理的反問題，所以在計(jì)算時(shí)具有不穩(wěn)定性.在形式上，向下延拓的不穩(wěn)定性在頻率域換算中表現(xiàn)為轉(zhuǎn)換因子對(duì)數(shù)據(jù)中高頻成份的放大作用，因此，一些研究就是基于頻率域里附加濾波器或改造下延轉(zhuǎn)換因子的手段，對(duì)造成不穩(wěn)定的高頻放大作用進(jìn)行一定的壓制，盡管該方式已經(jīng)取得了良好的效果，但還難以達(dá)到高精度勘探階段人們更高的期望.針對(duì)上述問題，近年來我國研究者對(duì)迭代法進(jìn)行了深入研究，取得了一些比較深入的研究成果.迭代法首先是Strakhov等[1]在研究曲化平時(shí)提出來的，他為了實(shí)現(xiàn)弱地形（即場(chǎng)源位于起伏地形最低點(diǎn)之下）上異常向下延拓到最低點(diǎn)平面，采取并系統(tǒng)地建立了積分迭代法，給出了計(jì)算迭代步驟，指出了迭代的收斂性分析.只是受當(dāng)時(shí)計(jì)算條件的約束，迭代法計(jì)算量相當(dāng)大限制了其應(yīng)用.徐世浙等[2－4]在此基礎(chǔ)上也進(jìn)行了系統(tǒng)的研究，并采取頻率域方法提高迭代法的計(jì)算速度，充分體現(xiàn)了迭代法門檻低的突出優(yōu)點(diǎn)，從而使迭代法受到普遍重視.劉東甲等[5]基于頻率域從數(shù)學(xué)上證明了迭代法向下延拓的等價(jià)計(jì)算關(guān)系式，依此方式可以直接計(jì)算等效的迭代結(jié)果，避免了前人方法的實(shí)際迭代過程，因而明顯提高了計(jì)算速度.另外，證明了迭代向下延拓計(jì)算隨著迭代次數(shù)的增加，所得到的結(jié)果無限趨近于FFT直接向下延拓的理論結(jié)果，依此作為迭代法下延收斂的重要依據(jù).張輝等[6]也對(duì)迭代向下延拓的收斂性進(jìn)行了類似的分析證明.曾小牛等[7]基于正則化手段壓制干擾信號(hào)高頻影響的實(shí)際需要，比較了迭代Tikhonov正則化法、Landweber正則化迭代法和積分迭代法，對(duì)提高迭代法下延中減弱高頻干擾的效果進(jìn)行了分析研究.于波等[8]也對(duì)實(shí)際資料中噪聲干擾的存在對(duì)迭代下延的影響進(jìn)行了分析研究.王彥國等[9]提出了泰勒級(jí)數(shù)迭代法，該方法主要是加快了迭代法計(jì)算時(shí)迭代過程的收斂速度.駱遙[10]通過對(duì)迭代下延中場(chǎng)分離的過程的分析，進(jìn)一步對(duì)迭代過程的收斂性加以解釋說明.由于迭代下延計(jì)算方法簡(jiǎn)單、易于實(shí)現(xiàn)，并且在很多情況下能取得相當(dāng)好的效果，所以，近年來該方法在實(shí)際中得到了普遍應(yīng)用（徐世浙等[11－14]；楊金玉等[15]；于波等[16－17]；王小兵等[18]；肖鋒等[19]；李成立等[20]；余海龍等[21]）.由此可見，迭代法在重磁位場(chǎng)下延中應(yīng)用研究已受到廣泛的關(guān)注和重視.

迭代法在向下延拓中的積極研究也將其影響擴(kuò)展到重磁位場(chǎng)轉(zhuǎn)換的其它方法中，比如磁異常的低緯度化極.化極轉(zhuǎn)換方法也稱化向地磁極（簡(jiǎn)稱化極）是最基本的磁測(cè)資料處理方法之一，其主要作用是消除地磁場(chǎng)對(duì)場(chǎng)源斜磁化造成的磁異常的復(fù)雜性，與原始磁異常相比，化極計(jì)算后的磁異常在形態(tài)上更加簡(jiǎn)單，與場(chǎng)源的分布對(duì)應(yīng)關(guān)系上也更加直接，因而成為磁異常分析解釋的基礎(chǔ).但是，在低緯度地區(qū)測(cè)量獲得的磁異常，其化極計(jì)算是不穩(wěn)定的，在頻率域中表現(xiàn)為化極轉(zhuǎn)換因子具有明顯的放大作用，在磁赤道附近，化極計(jì)算甚至因過分放大導(dǎo)致無法進(jìn)行.針對(duì)低緯度化極問題的長期研究進(jìn)展，這里不再一一列出，雖然取得了一些比較有效的成果，但仍然存在各種問題或不足.這里我們只對(duì)與迭代法相關(guān)的研究進(jìn)展加以介紹.最早文獻(xiàn)初步涉及到迭代法的是Keating等[22]的工作，他們嘗試運(yùn)用維納濾波的化極結(jié)果反算到原斜磁化方向與原始數(shù)據(jù)進(jìn)行對(duì)比，對(duì)化極結(jié)果進(jìn)行修正，不過實(shí)驗(yàn)及分析表明，當(dāng)修正次數(shù)超過3次時(shí)，結(jié)果不會(huì)隨著迭代次數(shù)的增加而改善.Li[23]運(yùn)用能量平衡法將化極結(jié)果反算到斜磁化方向與原始數(shù)據(jù)對(duì)化極結(jié)果修正，相當(dāng)于應(yīng)用了迭代過程，但其也只進(jìn)行了3次迭代.駱遙等[24]在借鑒迭代法下延取得成果的基礎(chǔ)上進(jìn)行了赤道化極研究，其重要結(jié)論是認(rèn)識(shí)到化極計(jì)算不但改變了異常的振幅而且改變了相位，所以直接通過計(jì)算值與實(shí)測(cè)值的殘差進(jìn)行修正的迭代法會(huì)導(dǎo)致結(jié)果不收斂，為此，通過技術(shù)性處理，即結(jié)合磁異?！暗瓜唷碧幚矸椒ǎǚ接瓐虻萚25]），實(shí)現(xiàn)了迭代法磁赤道化極，取得了比較好的效果.對(duì)于重磁位場(chǎng)處理轉(zhuǎn)換中迭代法的應(yīng)用而言，針對(duì)性的“倒相”轉(zhuǎn)換映射措施是對(duì)迭代法應(yīng)用研究的進(jìn)一步深化.

盡管迭代法應(yīng)用于向下延拓及低緯度化極等轉(zhuǎn)換計(jì)算已取得了一定的效果，但隨著研究的逐步深入，實(shí)際上已經(jīng)隱約暴露出迭代法應(yīng)用于重磁位場(chǎng)處理轉(zhuǎn)換的不足和問題.鑒于前人對(duì)迭代法的研究都是針對(duì)具體問題進(jìn)行孤立分析的，通過對(duì)前人研究成果加以分析、對(duì)比，可以發(fā)現(xiàn)，對(duì)于不同處理轉(zhuǎn)換問題，迭代法的收斂性和收斂條件是不盡相同的.所以，進(jìn)一步深入研究迭代法在重磁位場(chǎng)計(jì)算中的收斂性以及收斂條件的一般性規(guī)律，對(duì)于更好地應(yīng)用迭代法以及客觀評(píng)價(jià)其特點(diǎn)、認(rèn)識(shí)其存在的問題，具有重要意義.為此，本文對(duì)迭代法在重磁處理中的應(yīng)用進(jìn)行了統(tǒng)一的分析與研究，推導(dǎo)了迭代法通式，再針對(duì)具體問題推導(dǎo)出具體的收斂條件，并對(duì)各種因素對(duì)迭代法收斂性的影響進(jìn)行了分析.我們的研究表明，向下延拓及化極等不同的重磁處理轉(zhuǎn)換方法中需要選擇不同的迭代計(jì)算方式，相當(dāng)于是選擇了不同的映射函數(shù)；在選擇了合適的映射函數(shù)的情況下，如果迭代次數(shù)很大，迭代計(jì)算雖然能夠收斂到理論計(jì)算結(jié)果（在頻率域即FFT理論直接計(jì)算結(jié)果），但是一些處理轉(zhuǎn)換其理論計(jì)算結(jié)果本身通常是不穩(wěn)定的；如果迭代次數(shù)過少，則可能對(duì)有用信息的壓制過大，造成計(jì)算結(jié)果仍然存在較大的誤差；頻率域迭代法等效濾波計(jì)算關(guān)系式的存在證明了迭代法位場(chǎng)轉(zhuǎn)換在本質(zhì)上與改造轉(zhuǎn)換因子或附加特定濾波器的作用是等效的，所以，針對(duì)向下延拓及低緯度化極這類不穩(wěn)定的轉(zhuǎn)換計(jì)算，迭代法并沒有從根本上解決計(jì)算精度與穩(wěn)定性之間的矛盾.下面我們將對(duì)此進(jìn)行詳細(xì)的分析.

2 方法原理

重磁異常轉(zhuǎn)換計(jì)算中的迭代法應(yīng)用是一種通過多次的穩(wěn)定計(jì)算來代替不穩(wěn)定直接計(jì)算的過程，這點(diǎn)上明顯不同于其它領(lǐng)域里常用的迭代方式.以向下延拓為例，就是原本直接的下延計(jì)算，轉(zhuǎn)化為特定的分步迭代計(jì)算過程，每步迭代計(jì)算的主要過程是向上延拓，而不是向下延拓，這是和其它領(lǐng)域常規(guī)迭代法的明顯區(qū)別.現(xiàn)有文獻(xiàn)迭代法下延的具體迭代步驟為：首先將原始觀測(cè)面上的重磁異常直接當(dāng)成是待下延位置平面上的結(jié)果重磁異常，顯然，這只是一個(gè)初始近似，其近似程度（即誤差）是由延拓的高度差及重磁異常復(fù)雜程度共同決定的.遺憾的是其具體誤差是未知的，即無法評(píng)價(jià)其近似程度.為此，第二步是為評(píng)價(jià)及修改服務(wù)的，采取的是將該延拓面上的重磁異常初始近似結(jié)果再次向上延拓到原始觀測(cè)面上，這一步向上延拓是需要實(shí)際計(jì)算的.由于向上延拓是穩(wěn)定轉(zhuǎn)換計(jì)算，且具有很高的精度，因此這個(gè)計(jì)算過程引入的誤差很小，比較而言可以忽略.這樣，原始觀測(cè)面上就獲得了新的由上延計(jì)算得到的結(jié)果，原始數(shù)據(jù)與其對(duì)比的差異，一定程度上可以看成是下延平面上未知的理想結(jié)果與近似結(jié)果的差異的反映，基于此間接達(dá)到了對(duì)下延結(jié)果進(jìn)行近似評(píng)價(jià)的目的.如果該差異（即誤差）足夠小，則認(rèn)為下延的結(jié)果已足夠逼近理想結(jié)果，否則，需要設(shè)法提高、改進(jìn).這里，迭代法的改進(jìn)方式就是采取將原始數(shù)據(jù)與上延結(jié)果的差直接附加于待延拓面上的初始異常，附加修改后獲得的新的近似異常，則需要重復(fù)采用上延的方法計(jì)算與評(píng)價(jià)，這樣就形成了新的迭代過程.由此可見，這種迭代法下延是希望通過分步的上延計(jì)算（穩(wěn)定計(jì)算）代替一步的直接下延計(jì)算（不穩(wěn)定計(jì)算），其突出優(yōu)點(diǎn)是原理清晰、步驟簡(jiǎn)單.眾多研究者也都強(qiáng)調(diào)該計(jì)算過程的穩(wěn)定性，并以最終逼近理論直接下延結(jié)果證明其收斂性.

為了客觀評(píng)價(jià)迭代法在重磁位場(chǎng)轉(zhuǎn)換中的作用，有必要進(jìn)一步進(jìn)行歸納分析.下面我們?cè)陬l率域?qū)Φǖ牧鞒踢M(jìn)行說明，以期得到一般性的迭代轉(zhuǎn)換關(guān)系式.

設(shè)U0（u，v）為觀測(cè)平面上原始數(shù)據(jù)u0（x，y）的Fourier變換結(jié)果，Ua（u，v）為待求的目標(biāo)數(shù)據(jù)ua（x，y）的Fourier變換結(jié)果，其中x，y為空間域坐標(biāo)變量，u，v為頻率域頻率變量，ψ（u，v）為相應(yīng)轉(zhuǎn)換計(jì)算的頻率域轉(zhuǎn)換因子，則二者在頻率域轉(zhuǎn)換關(guān)系可以表示為

對(duì)于化極、向下延拓等轉(zhuǎn)換問題，（1）式中ψ（u，v）的一些值會(huì)很大甚至為無窮，直接計(jì)算是不穩(wěn)定或無法計(jì)算，傳統(tǒng)的針對(duì)性方法多數(shù)是對(duì)該轉(zhuǎn)換因子進(jìn)行改造，提高其計(jì)算的穩(wěn)定性.對(duì)于采用迭代法計(jì)算，為了得到穩(wěn)定的計(jì)算結(jié)果，首先直接將原始數(shù)據(jù)當(dāng)成是目標(biāo)數(shù)據(jù)的初始值（對(duì)于向下延拓）或“倒相”后作為初始值（對(duì)于赤道化極），這里我們將其作一般性推廣，即定義為一個(gè)U0（u，v）到的映射為目標(biāo)數(shù)據(jù)頻譜Ua（u，v）的近似值，其形式為

在頻率域里，由于我們所要解決的問題是線性的，所以映射也選擇線性映射.我們稱φ（u，v）為映射函數(shù)，理論上，映射函數(shù)在選擇上應(yīng)該是頻率域轉(zhuǎn)換因子的一個(gè)近似.

圖1 迭代法流程圖Fig.1 The scheme of iteration method

迭代法流程如圖1所示.為了求解（1）式中的Ua（u，v），我們?cè)O(shè)（2）式的是上述方程的近似解，但其近似程度實(shí)際上是無法評(píng)價(jià)的，為此，這里采取的措施是將其反算回去與原始數(shù)據(jù)對(duì)比，試圖由此間接地評(píng)價(jià)其近似程度，顯然，其前提是反運(yùn)算的誤差不大.設(shè)Ur（u，v）是原始數(shù)據(jù)與反運(yùn)算結(jié)果的殘差，則有

這里先假設(shè)通過選取合適的φ（u，v）使得結(jié)果能夠收斂，具體步驟及過程在后面會(huì)詳細(xì)論述，則雖然仍然不是方程的精確解，但是的改進(jìn).可以把再次反算回與原始數(shù)據(jù)對(duì)比、評(píng)價(jià)，如此，則形成迭代過程，直到誤差足夠小為止.

根據(jù)上面的假設(shè)，我們可以總結(jié)出迭代法的遞推公式：

下面對(duì)上述第n次的迭代結(jié)果進(jìn)行歸納推導(dǎo)：第一次迭代：

第二次迭代：

根據(jù)（3）、（4）、（5）式以及數(shù)學(xué)歸納法不難得出：第n次迭代結(jié)果為

根據(jù)等比數(shù)列求和公式將其化簡(jiǎn)可得

從上式可以看出只要

即對(duì)于所有的u、v都有下式成立：

迭代結(jié)果就會(huì)收斂到ψ（u，v）U0（u，v）.

由此可見，映射函數(shù)φ（u，v）決定著迭代法是否收斂.在收斂的前提下，迭代無窮次的結(jié)果為FFT理論直接計(jì)算結(jié)果ψ（u，v）U0（u，v），又與我們選擇的映射函數(shù)無關(guān).而當(dāng)有限次迭代時(shí)，收斂的速度受所選擇的映射函數(shù)控制.下面我們將結(jié)合具體處理轉(zhuǎn)換方法，分析對(duì)應(yīng)的迭代計(jì)算中如何選擇映射函數(shù)以及評(píng)價(jià)具體的迭代收斂特性.

3 不同處理轉(zhuǎn)換問題的收斂性具體分析

上面推導(dǎo)的迭代法通式中，映射函數(shù)如何選擇及收斂性如何是關(guān)鍵因素，下面我們將結(jié)合不同的具體處理轉(zhuǎn)換加以分析.

3.1 迭代延拓

3.1.1 迭代向下延拓

在頻率域，向下延拓可以表示為

其中：U0（u，v）是觀測(cè)面上異常的頻譜，Uh（u，v）是向下延拓結(jié)果的頻譜，s為徑向波數(shù)，h為下延的高度差，ψ（u，v）＝esh為下延轉(zhuǎn)換因子，是典型的高頻放大濾波器，高度差較大時(shí)，計(jì)算受高頻放大影響嚴(yán)重，計(jì)算結(jié)果很不穩(wěn)定.

根據(jù)上面的推導(dǎo)可知，對(duì)于迭代向下延拓計(jì)算，其逆變換因子ψ－1（u，v）＝e－sh，即為向上延拓因子，帶入收斂條件式（8）可得：｜1－φ（u，v）e－sh｜＜1，

進(jìn)一步化簡(jiǎn)可得

式（11）即為迭代法向下延拓中，為保證迭代收斂，所選擇的映射函數(shù)必須滿足的條件.在滿足該條件的情況下，如何選擇映射函數(shù)可以基于不同的考慮，如收斂速度、計(jì)算過程簡(jiǎn)捷等.

分析前人現(xiàn)有的迭代法向下延拓研究成果可以發(fā)現(xiàn)，實(shí)際上大家均隱含著選擇了最簡(jiǎn)單的常數(shù)為1的映射函數(shù).這里為了計(jì)算簡(jiǎn)便，也不影響一般性，我們?cè)谇叭嘶A(chǔ)上稍作擴(kuò)展，取映射函數(shù)φ（u，v）＝m，m為實(shí)常數(shù)，我們稱其為迭代速度系數(shù)（后面可以看出它的這個(gè)特點(diǎn)），則收斂條件可以轉(zhuǎn)化為

這里記區(qū)間S＝（0，2），則S為迭代向下延拓速度系數(shù)的收斂區(qū)間.下面我們來分析速度系數(shù)m、迭代次數(shù)n以及向下延拓高度差h對(duì)濾波器的影響情況.

由（12）式可見，迭代法下延計(jì)算等價(jià)于對(duì)直接下延轉(zhuǎn)換因子乘上一個(gè)附加濾波器，所以，在本質(zhì)上與傳統(tǒng)的針對(duì)下延不穩(wěn)定性而采取的通過附加穩(wěn)定濾波器，從而達(dá)到改造下延轉(zhuǎn)換因子、實(shí)現(xiàn)穩(wěn)定下延的作用是一樣的.

下面我們可以進(jìn)一步分析評(píng)價(jià)迭代下延所附加的濾波器的頻譜特征.將（12）式寫成下面的形式：

針對(duì)（13）式這樣的迭代下延附加濾波器，顯然其特性受三個(gè)參數(shù)控制，即迭代次數(shù)n、下延高度差h和速度系數(shù)m.為簡(jiǎn)單直觀起見，下面用圖示方式對(duì)這三個(gè)參數(shù)的變化對(duì)濾波器的影響進(jìn)行表示.

圖2 迭代下延附加濾波器頻譜曲線Fig.2 Spectrum curve of the iterative downward continuation filter

圖2分別顯示了速度系數(shù)、迭代次數(shù)、下延高度對(duì)迭代向下延拓濾波器的影響，由圖我們可以得到以下認(rèn)識(shí)和結(jié)論：

（1）迭代下延附加濾波器為典型的低通濾波器，所以，迭代下延相當(dāng)于在直接下延計(jì)算中附加了低通濾波作用，使得結(jié)果受高頻影響減小.該濾波器只是在極限情況下即迭代次數(shù)無窮大時(shí)，其數(shù)值趨近于1，即沒有任何附加濾波作用，也就是不對(duì)轉(zhuǎn)換計(jì)算產(chǎn)生任何影響.

（2）在迭代次數(shù)與下延深度不變的情況下，速度系數(shù)越大，濾波器通帶范圍越大（即向高頻移動(dòng)），其作用是速度系數(shù)越大迭代下延結(jié)果與FFT直接下延結(jié)果越接近，故較大的系數(shù)可以使迭代下延用較少的迭代次數(shù)達(dá)到收斂，但是當(dāng)速度系數(shù)過大時(shí)，會(huì)對(duì)低頻數(shù)據(jù)造成一定的影響，分析如下（示意如圖3）：

圖3 迭代次數(shù)奇偶性對(duì)濾波器的影響Fig.3 The effects of iterative number′s parity on the filter

（13）式所表示的濾波器可以認(rèn)為是y＝1－（1－x）n的函數(shù)，當(dāng)m趨近于2時(shí)，對(duì)于低頻信號(hào)而言，x＝me－sh會(huì)趨近于2，若此時(shí)迭代次數(shù)n不是很大，則y會(huì)隨著迭代次數(shù)選擇奇數(shù)及偶數(shù)的交替變化而發(fā)生震蕩現(xiàn)象（如圖3所示），所以此時(shí)的濾波器在壓制高頻信號(hào)的同時(shí)也對(duì)低頻信號(hào)產(chǎn)生了一定的影響.將圖2（a）的低頻部分局部放大，會(huì)得到如圖4所示的速度系數(shù)值對(duì)低頻的具體影響.

圖4 不同速度系數(shù)m對(duì)低頻的影響（迭代次數(shù)為50）Fig.4 The effects of velocity factor on low－frequency（50iterations）

根據(jù)上面的分析可知，隨著迭代次數(shù)的增加，迭代向下延拓雖然收斂到FFT向下延拓的結(jié)果，但是當(dāng)?shù)螖?shù)不大時(shí)，計(jì)算結(jié)果隨迭代波動(dòng)，表現(xiàn)為誤差曲線呈鋸齒狀震蕩（后面有圖例）.

（3）在下延深度與速度系數(shù)不變的情況下，迭代次數(shù)越多，濾波器通帶高通范圍越大，濾波器越接近“不濾波”，即結(jié)果越趨近于理論下延結(jié)果.

（4）在迭代次數(shù)與速度系數(shù)不變的情況下，下延深度小，附加濾波器通帶高通范圍越大，反之下延深度大，則濾波器主要表現(xiàn)為低通，可見下延深度越大，迭代下延需要更多的迭代次數(shù)才能與FFT直接下延結(jié)果越接近.

由此可見，當(dāng)速度系數(shù)處于穩(wěn)定區(qū)間之內(nèi)時(shí)，迭代次數(shù)與速度系數(shù)作用相似，均可以控制低通濾波的通帶范圍.

從上面的認(rèn)識(shí)我們可以看出，迭代向下延拓是逼近FFT直接向下延拓的過程，其本質(zhì)上為對(duì)直接下延進(jìn)行了一次低通濾波，從而壓制通常會(huì)引起不穩(wěn)定的高頻信息，使得在較少的迭代次數(shù)下表現(xiàn)為良好的計(jì)算穩(wěn)定性.另一方面由于其與其它壓制下延因子濾波器的方法在本質(zhì)上是一樣的，所以，并沒有從根本上徹底解決向下延拓問題中精度與穩(wěn)定性的矛盾.

3.1.2 迭代向上延拓

眾所周知，位場(chǎng)向上延拓是滿足Laplace方程的Dirichlet問題，屬于數(shù)學(xué)物理領(lǐng)域中的正問題，具有穩(wěn)定性，只要位場(chǎng)異常數(shù)據(jù)采樣率足夠密集，異常形態(tài)足夠完整，則上延計(jì)算的精度將可達(dá)很高精度.在頻率域向上延拓計(jì)算中，因上延轉(zhuǎn)換因子隨頻率變化其值始終是穩(wěn)定的，可以通過直接計(jì)算得到穩(wěn)定的向上延拓結(jié)果，所以，不需要選擇迭代法進(jìn)行計(jì)算.這里我們僅僅為了研究迭代法的適應(yīng)能力以及探索其各方面的特性，所以嘗試用迭代法來解向上延拓問題，便于對(duì)迭代法在位場(chǎng)轉(zhuǎn)換中的應(yīng)用有一個(gè)更全面和深入的了解.

與迭代向下延拓類似，易得迭代向上延拓的收斂條件為

如果選擇映射函數(shù)為簡(jiǎn)單的實(shí)常數(shù)（即速度系數(shù)），則對(duì)應(yīng)的收斂條件需滿足

上式中的smax是原始數(shù)據(jù)徑向波數(shù)的最大值，h仍為高度差.

迭代上延的通式與迭代下延形式類似，其表達(dá)式為

上式中的U－h(huán)（u，v）為向上延拓理論結(jié)果的譜，為迭代向上延拓n次結(jié)果.對(duì)（16）式進(jìn)行分析表明，迭代向上延拓相當(dāng)于在理論向上延拓的基礎(chǔ)上進(jìn)行一次高通濾波，所以，其計(jì)算結(jié)果理論上容易變得不穩(wěn)定.并且，只有m為足夠小的正實(shí)數(shù)時(shí)，迭代向上延拓才是收斂的，當(dāng)?shù)螖?shù)趨近于無窮大時(shí)，才能收斂到理論的向上延拓結(jié)果.

這里假設(shè)按照前人迭代下延的方法進(jìn)行迭代上延計(jì)算，即進(jìn)一步選擇映射函數(shù)為常數(shù)1，也就是首先將原始異常數(shù)據(jù)當(dāng)成為上延平面上的初次結(jié)果，因?yàn)閷?duì)其無法進(jìn)行直接評(píng)價(jià)，需要再將其下延到原平面進(jìn)行對(duì)比評(píng)價(jià)，將評(píng)價(jià)誤差再修正初始結(jié)果，采取逐步迭代修正方式進(jìn)行.顯然，（15）式表明m＝1在大多數(shù)條件下極有可能是不收斂的（如h稍大，或smax稍大），這是迭代上延與迭代下延的明顯區(qū)別，后面的理論模型計(jì)算也證明了這一點(diǎn).

圖5展示的是上延高度差、迭代次數(shù)以及速度系數(shù)對(duì)迭代上延的效果影響的情況.

圖5 迭代上延附加濾波器頻譜曲線Fig.5 Spectrum curve of the iterative upward continuation filter

從上面迭代上延可知，迭代法也可以應(yīng)用于解決穩(wěn)定的位場(chǎng)計(jì)算問題，但是如果原本穩(wěn)定計(jì)算的逆運(yùn)算是不穩(wěn)定的，那么對(duì)映射函數(shù)的選擇則要苛刻得多，否則選擇迭代法就會(huì)將原本穩(wěn)定計(jì)算問題轉(zhuǎn)化為不穩(wěn)定計(jì)算問題.

迭代下延與迭代上延的特性分析，為我們更全面地了解迭代法位場(chǎng)轉(zhuǎn)換的性質(zhì)，掌握具體的應(yīng)用方法打下了理論基礎(chǔ).

3.1.3 迭代曲化平

曲化平是將起伏觀測(cè)面上獲得的重磁異常轉(zhuǎn)換到某一平面上等效異常的計(jì)算方法，目的是消除觀測(cè)面起伏的影響.可以看出，向下延拓和向上延拓即平面之間的延拓是其特例，所以理論上同樣可以運(yùn)用迭代法來解決曲化平計(jì)算問題.

Strakhov等[1]及徐世浙等[4]將迭代法應(yīng)用于弱地形的曲化平計(jì)算，簡(jiǎn)化了曲化平的復(fù)雜性并取得了比較好的曲化平計(jì)算結(jié)果.如果從算法角度評(píng)價(jià)該方法，通過我們上面迭代法下延及迭代法上延分析可以理解，該曲化平迭代過程相當(dāng)于選擇了等價(jià)的映射函數(shù)為常數(shù)1，遺憾的是由于曲化平在頻率域沒有直接的解析表達(dá)式，所以，無法進(jìn)一步評(píng)價(jià)其等效的濾波特性了.

雖然迭代法曲化平計(jì)算過程簡(jiǎn)單并在條件較好時(shí)取得了不錯(cuò)的效果，但作者也指出其受地形的起伏大小、下延到地形最低點(diǎn)的高差等影響較大（注：迭代曲化平中間過程最重要一步是將曲面上的數(shù)據(jù)下延到地形最低點(diǎn)所在的平面上，所以，必須是場(chǎng)源較深的弱地形條件）.從理論上分析其原因，可以認(rèn)為是地形起伏大以及高差大，決定了未知的場(chǎng)源在起伏大的地形上產(chǎn)生的異常數(shù)據(jù)與在高差大的待轉(zhuǎn)化平面上產(chǎn)生的異常數(shù)據(jù)比較，無論在數(shù)值大小還是異常形態(tài)上均相差比較大，如果分析其對(duì)應(yīng)頻譜，則曲面數(shù)據(jù)的譜可能與平面數(shù)據(jù)的譜在振幅和相位都存在明顯的差異，此時(shí)仍然選取映射為1顯然是不合理的，特別是因?yàn)橄辔徊顒e大的情況，這在下面化極中表現(xiàn)得更加明顯.

由于曲化平?jīng)]有簡(jiǎn)單的迭代法解析表達(dá)式，這里就不對(duì)其迭代法曲化平的映射函數(shù)的選取及收斂性等進(jìn)行量化分析了.

3.2 迭代化極

化極轉(zhuǎn)換是最基本的磁測(cè)資料處理方法之一，所以，低緯度化極的困難一直受到重點(diǎn)研究.下面我們對(duì)迭代法在化極中的應(yīng)用加以類似的分析研究，以期掌握其客觀內(nèi)在規(guī)律.

在頻率域中，化極處理可以表示為

其中：U0（u，v）是斜磁化條件下異常的頻譜，Up（u，v）是垂直磁化條件下磁異常垂直分量的頻譜，I和D分別是地磁場(chǎng)的傾角和偏角且，這里我們假設(shè)磁化方向與地磁場(chǎng)方向一致.u、v是x、y方向上的波數(shù)，θ＝arctan（u/v），i是虛數(shù)單位.圖6所示為化極因子頻譜曲線，在I＝0的水平磁化條件下，當(dāng)θ→D±時(shí)（此圖中D＝0°，則θ＝90°和270°），頻率域化極因子的幅度急速上升，頻譜被無限放大，化極因子本身出現(xiàn)奇異，無法獲得穩(wěn)定的化極結(jié)果.

圖6 常規(guī)化極因子不同磁傾角的頻譜曲線及其放大作用Fig.6 The effects of inclination on the RTP factor

根據(jù)前面的推導(dǎo)可知，采用迭代法化極計(jì)算，化極的逆過程我們稱其為“逆化極”，對(duì)應(yīng)的濾波器為“逆化極”因子，其形式為

帶入（9）式可得

在迭代法中，φ（θ）取實(shí)常數(shù)更方便（如前面的迭代延拓），也更能體現(xiàn)迭代法易簡(jiǎn)單實(shí)現(xiàn)的優(yōu)點(diǎn)，考慮到化極計(jì)算比延拓復(fù)雜，所以這里取其為實(shí)數(shù)m（θ）.則將上式展開可得

我們?cè)O(shè)

對(duì)（19）式求解，可得收斂條件為：

或者：

由于θ是連續(xù)變化的，即cos2（D－θ）可以取到[0，1]區(qū)間內(nèi)的任意值，下面我們討論不同磁傾角情況下的收斂條件：

（1）當(dāng)I＞45°時(shí)，對(duì)于所有的θ都有：

所以對(duì)于所有的θ都有C（θ）＜0，所以此時(shí)的收斂條件為

上式中的Cmin為C（θ）的最小值.

（2）當(dāng)0°＜I≤45°時(shí)，C（θ）不全為負(fù)值，還可能存在0值或者正值，顯然當(dāng)m（θ）為常數(shù)時(shí)無法滿足收斂條件.如果要在此情況下保持收斂，那么m（θ）必須為一個(gè)隨頻率變化、需要根據(jù)具體的磁傾角和磁偏角加以針對(duì)性設(shè)計(jì)的變量了，其細(xì)節(jié)不是我們這里討論的重點(diǎn)，因此，這里就不展開詳細(xì)討論了.

（3）當(dāng)I＝0時(shí)（即水平磁化環(huán)境，相當(dāng)于位于磁赤道地區(qū)），C（θ）＝2，所以易得此時(shí)的收斂條件為

這種情況下，我們可以進(jìn)一步令m（θ）＝m為實(shí)常數(shù)，則I＝0收斂條件如下：

通過該收斂條件可知，當(dāng)取m＝1時(shí)，磁赤道迭代化極是不收斂的，即直接將實(shí)測(cè)磁異常當(dāng)成化極結(jié)果初始值這樣逐步迭代其結(jié)果將不收斂.駱遙等（2010）在針對(duì)赤道化極研究中，注意到赤道化極前后異常相位的倒相關(guān)系，采取每次迭代都對(duì)殘差“倒相”的處理，實(shí)現(xiàn)了迭代赤道化極，其對(duì)應(yīng)的速度系數(shù)為這里的m＝－1，恰好在我們的收斂區(qū)間式（20）內(nèi)，所以，這里對(duì)駱遙等（2010）研究中的經(jīng)驗(yàn)認(rèn)識(shí)從理論上給予了證明.

下面我們簡(jiǎn)單分析迭代次數(shù)n以及速度系數(shù)m對(duì)濾波器的影響：設(shè)為第n次迭代的結(jié)果，Up是FFT理論直接化極結(jié)果，則迭代化極等價(jià)于在FFT直接化極的基礎(chǔ)上進(jìn)行附加濾波，濾波器形式為

圖7分別顯示了速度系數(shù)、迭代次數(shù)對(duì)迭代化極濾波器的影響，經(jīng)過分析，我們可以得到如下結(jié)論：

圖7 迭代化極附加濾波器頻譜曲線Fig.7 Spectrum curve of the iterative RTP filter

（1）在迭代次數(shù)和磁傾角不變的情況下，速度系數(shù)絕對(duì)值越大，濾波器通帶范圍越大，即增大速度系數(shù)的絕對(duì)值可以使結(jié)果更趨近于FFT直接化極結(jié)果，并且增加收斂速度；但是和迭代向下延拓存在的問題類似，速度系數(shù)同樣存在穩(wěn)定區(qū)間，速度系數(shù)在穩(wěn)定區(qū)間之外時(shí)，計(jì)算會(huì)不穩(wěn)定.

為簡(jiǎn)單起見，當(dāng)取磁傾角、磁偏角均為0時(shí)，濾波器函數(shù)記為y＝1－（1＋mcos2（θ））n，其中θ取0到2π，不同速度系數(shù)的濾波曲線如圖8所示.

圖8 速度系數(shù)迭代次數(shù)對(duì)計(jì)算穩(wěn)定性的影響Fig.8 The effects of velocity factor and iterative number on the calculation stability

從上面圖8中可以看出，當(dāng)速度系數(shù)在穩(wěn)定區(qū)間之內(nèi)時(shí)（圖8（a）、圖8（b）），迭代次數(shù)的奇偶變化對(duì)濾波曲線的影響不大；而當(dāng)速度系數(shù)接近收斂區(qū)域的邊界，其穩(wěn)定性明顯降低，表現(xiàn)為迭代次數(shù)的奇偶變化會(huì)導(dǎo)致濾波曲線發(fā)生明顯變化，從而將導(dǎo)致計(jì)算結(jié)果發(fā)生震蕩（圖8（c）、圖8（d））.

（2）在速度系數(shù)不變的情況下，迭代次數(shù)越大，濾波器通帶范圍越大，迭代化極結(jié)果與FFT直接化極結(jié)果越相近，且當(dāng)?shù)螖?shù)趨近于無窮時(shí)，迭代化極的結(jié)果也趨近于理論化極的結(jié)果.

因此，迭代化極的本質(zhì)同樣是在原有化極因子的基礎(chǔ)上加上一個(gè)帶通濾波器.直接化極因子主要在相位處具有明顯的放大作用，而迭代法化極的附加濾波器恰好在這些地方予以壓制，故能使得計(jì)算結(jié)果更加穩(wěn)定.當(dāng)速度系數(shù)處于穩(wěn)定區(qū)間內(nèi)時(shí)，迭代次數(shù)與速度系數(shù)的絕對(duì)值對(duì)收斂性的貢獻(xiàn)的作用相似.

雖然我們得到了迭代法化極的關(guān)系式，并且能避免真正的迭代計(jì)算步驟，但是迭代法化極并沒有從本質(zhì)上解決低緯度化極因子奇異的問題，式（21）表明：迭代次數(shù)非常大時(shí)，迭代化極回歸為直接化極，顯然，迭代法化極是對(duì)直接化極的一種濾波近似，主要是壓制化極中不穩(wěn)定的頻率成份，因此，和其它針對(duì)低緯度化極采取的改造化極因子或附加其它濾波器是沒有本質(zhì)區(qū)別的.這里還應(yīng)該指出，迭代法的迭代次數(shù)是不好選擇的一個(gè)參數(shù)，所以，相對(duì)于其它方法，迭代法并不具有優(yōu)勢(shì).

4 數(shù)據(jù)實(shí)驗(yàn)分析

為了對(duì)迭代法有一個(gè)更加直觀的認(rèn)識(shí)，我們選取迭代延拓和迭代化極方法進(jìn)行模型計(jì)算實(shí)驗(yàn)，對(duì)迭代法的收斂性等加以分析評(píng)價(jià)，以期為合理應(yīng)用迭代法打下基礎(chǔ).

4.1 迭代下延模型實(shí)驗(yàn)

這里我們選取一個(gè)斜磁化球體模型的計(jì)算實(shí)驗(yàn)來對(duì)向下延拓計(jì)算的收斂性以及計(jì)算精度進(jìn)行分析.所選球體的磁化強(qiáng)度為0.5A/m，磁傾角45°，磁偏角15°，球體中心點(diǎn)的x、y坐標(biāo)均為0，埋深為2km，實(shí)驗(yàn)矩形網(wǎng)格數(shù)據(jù)參數(shù)為：測(cè)線數(shù)M＝301，每條測(cè)線上的測(cè)點(diǎn)數(shù)N＝301，點(diǎn)距、線距均為50m.其總場(chǎng)磁異常ΔT等值線圖見圖9所示，圖9a為0m平面上的磁異常，圖9b為1000m平面上的磁異常，即位于圖9a的上方1000m處.

針對(duì)該理論模型數(shù)據(jù)，我們應(yīng)用迭代法計(jì)算下延結(jié)果，其中誤差大小用均方差衡量，其計(jì)算式子為：

上式中ut（m，n）為理論結(jié)果，uc（m，n）為迭代計(jì)算的結(jié)果.

不同迭代次數(shù)計(jì)算結(jié)果如圖10所示.下延深度為20倍點(diǎn)距，即由高度為1000m的圖9b異常數(shù)據(jù)，迭代下延得到高度為0m的異常數(shù)據(jù)，選擇映射函數(shù)為1，即采取與前人相同的速度系數(shù).由結(jié)果圖可見，該模型迭代向下延拓計(jì)算時(shí)，在迭代次數(shù)較少時(shí)，所得的結(jié)果是穩(wěn)定的；而迭代次數(shù)越大，數(shù)據(jù)的高頻成分放大越明顯，所得結(jié)果與FFT理論直接下延結(jié)果越接近.

圖10 不同迭代次數(shù)迭代向下延拓結(jié)果Fig.10 Iterative downward continuation results with different iterative numbers

對(duì)該例而言，分析計(jì)算時(shí)達(dá)到最小誤差的迭代次數(shù)以及最小誤差與速度系數(shù)、下延深度的關(guān)系，可得結(jié)論如下：

（1）在下延深度一定且速度系數(shù)處于穩(wěn)定區(qū)間的情況下，達(dá)到最小誤差的迭代次數(shù)隨著速度系數(shù)的增大而減小，最小誤差值總體而言呈下降趨勢(shì).

（2）當(dāng)速度系數(shù)在穩(wěn)定區(qū)間之外時(shí)，迭代下延濾波器對(duì)高頻壓制的同時(shí)也會(huì)對(duì)低頻產(chǎn)生影響，表現(xiàn)為誤差曲線有明顯的震蕩.如圖11a所示，速度系數(shù)為1.0時(shí)，三條誤差曲線都是光滑的，無震蕩現(xiàn)象；但當(dāng)速度系數(shù)達(dá)到1.99時(shí)（見圖11d），不同下延深度的誤差曲線都會(huì)隨著迭代的進(jìn)行發(fā)生較為明顯的震蕩.

圖11 速度系數(shù)對(duì)迭代下延誤差曲線的影響Fig.11 The affected error curves by velocity factor

（3）在速度系數(shù)一定且處于穩(wěn)定區(qū)間內(nèi)的情況下，隨著下延深度的增加，所需要的迭代次數(shù)增加，同時(shí)最小誤差值也隨之增大.

（4）在向下延拓深度與速度系數(shù)不變的情況下，迭代下延的結(jié)果與理論結(jié)果的誤差隨著迭代次數(shù)的增加是先減小后又增大，相當(dāng)于有最佳迭代次數(shù)（或迭代區(qū)間），但由于理論結(jié)果是未知的，實(shí)際誤差無法準(zhǔn)確獲得，所以迭代法中只能間接評(píng)價(jià)誤差，其可靠性降低，因此在實(shí)際迭代下延中通常難以獲得最佳迭代次數(shù).

4.2 迭代上延模型實(shí)驗(yàn)

正如前面分析，實(shí)際中因?yàn)樯涎愚D(zhuǎn)換是穩(wěn)定計(jì)算，具有很高的精度，故沒有人會(huì)真正使用迭代法，那樣將使原本簡(jiǎn)單問題復(fù)雜化.這里僅僅是為了評(píng)價(jià)迭代法在上延應(yīng)用中如何表現(xiàn)，為此進(jìn)行必要的嘗試與特性分析.

針對(duì)上面的模型數(shù)據(jù)，現(xiàn)在采用迭代法將圖9a中0m平面的數(shù)據(jù)上延到1000m處.首先我們擬采取與迭代下延相同的方式，根據(jù)（15）式采取常數(shù)映射，但卻得不到實(shí)質(zhì)結(jié)果.因?yàn)榇藭r(shí)上延20倍點(diǎn)距，即1000m，smax＝0.089（m－1），速度系數(shù)取實(shí)常數(shù)值的范圍為：0＜m＜2.1582×10－39，即非常小，造成原始數(shù)據(jù)乘以該映射值后接近于0值了，所以，迭代的變化量將極其小，需要極其巨大的迭代次數(shù)才能得到結(jié)果.所以，基于常數(shù)映射的迭代法向上延拓計(jì)算，延拓的高度差不能太大.此結(jié)論在較小延拓高度差下得到了驗(yàn)證，如：由0m上延到50m，速度系數(shù)的取實(shí)常數(shù)值的范圍為：0＜m＜0.0117，比如取速度系數(shù)為0.01，即可得到下圖的收斂結(jié)果.

圖12即是取迭代速度系數(shù)為0.01，由0m上延到50m時(shí)的計(jì)算結(jié)果，其中圖12a為迭代1000次的上延數(shù)據(jù)分布圖，其結(jié)果與理論結(jié)果（圖略）的誤差為0.0002nT，其迭代中誤差收斂情況如圖12b所示，收斂還是很快的.但是，如果取迭代速度系數(shù)為1時(shí)，由于其不在收斂條件0＜m＜0.0117內(nèi)，實(shí)際計(jì)算也證明了其迭代是不收斂的，僅迭代5次其結(jié)果誤差已非常大了（圖13a），圖13b的收斂曲線可以看出其是發(fā)散的.

圖12 理論模型迭代上延0m到50m結(jié)果—速度系數(shù)取0.01Fig.12 Iterative upward continuation result（velocity factor 0.01，height 50m）and error curve

圖13 理論模型迭代上延0m到50m結(jié)果—速度系數(shù)取1.0Fig.13 Iterative upward continuation result（velocity factor 1，height 50m）and error curve

如果由0m平面迭代上延到250m高度平面，實(shí)際計(jì)算的成本將大大增加.此時(shí)smax＝0.089（m－1）不變，常數(shù)m 的取值范圍：0＜m＜2.1554×10－10，計(jì)算中取值1×10－10實(shí)驗(yàn)，最終迭代雖然收斂（見圖14），但迭代次數(shù)極大.

圖14 理論模型迭代上延0m到250m結(jié)果—速度系數(shù)取1.0×10－10Fig.14 Iterative upward continuation result（velocity factor 1，height 250m）and error curve

正如上面指出，按照同樣的方式，計(jì)算由0m到1000m，成本將無法承受，這是上延這樣的原本穩(wěn)定的直接計(jì)算，轉(zhuǎn)換為迭代法后如果采取簡(jiǎn)單映射處理方式造成的后果.這時(shí)，我們就不能取簡(jiǎn)單的實(shí)常數(shù)映射了，但如果采取變量形式，應(yīng)該易于解決遇到的新問題.

針對(duì)簡(jiǎn)單處理無法實(shí)現(xiàn)的0m到1000m的迭代上延，我們根據(jù)前面的（14）式的取值范圍重新取值，比如取φ（u，v）＝0.25e－sh，則只迭代了41次，誤差即達(dá)到0.00059nT，已完全收斂到非常接近理論結(jié)果的程度了，結(jié)果見圖15所示.

通過上面迭代上延的模型計(jì)算，使我們認(rèn)識(shí)到：在迭代下延中可以選擇簡(jiǎn)單的常數(shù)映射，但對(duì)于迭代上延，必須采取變量映射才能有效收斂，所以，映射函數(shù)的合理選取是迭代法重磁處理轉(zhuǎn)換的核心.

圖15 理論模型迭代上延0m到1000m結(jié)果—速度系數(shù)取φ（u，v）＝0.25e－shFig.15 Iterative upward continuation result（velocity factor 1，height 1000m）and error curve

上述迭代上延的試驗(yàn)與分析，目的不是應(yīng)用于實(shí)際上延計(jì)算（因?yàn)橹苯由涎愚D(zhuǎn)換已足夠精確），而是深化對(duì)迭代法的認(rèn)識(shí).

4.3 迭代化極數(shù)據(jù)試驗(yàn)

為了對(duì)迭代法在化極轉(zhuǎn)換中的應(yīng)用效果進(jìn)行評(píng)價(jià)，我們選取最困難的水平磁化條件下的化極例子進(jìn)行分析.下面選擇前人文獻(xiàn)中同樣的長方體模型來進(jìn)行化極計(jì)算實(shí)驗(yàn)，模型的長和寬都是20m，厚2m，上頂埋深1m，磁化強(qiáng)度為1A/m.數(shù)據(jù)網(wǎng)格參數(shù)為：點(diǎn)數(shù)M＝64，線數(shù)N＝64，點(diǎn)距Δx＝1m，線距Δy＝1m.該模型產(chǎn)生的理論磁異常如圖16所示.

圖16 理論模型總場(chǎng)異常ΔT（單位：nT）Fig.16 Theoretical total field anomalies

迭代化極的結(jié)果如圖17所示，計(jì)算中選擇映射函數(shù)為滿足收斂條件的常數(shù)－1.0（與駱遙等[24]（2010）中的“倒相”對(duì)應(yīng)）.迭代次數(shù)及不同映射函數(shù)（即速度系數(shù)）對(duì)該模型化極誤差影響見圖18.由這些結(jié)果圖可以看出：化極結(jié)果隨著迭代次數(shù)的增大誤差先減小然后又逐漸增大，直至無窮大.實(shí)際上赤道化極的FFT直接計(jì)算理論結(jié)果就是趨近于無窮大的，因此，迭代法赤道化極只是在較少的迭代次數(shù)時(shí)可以得到較好的化極結(jié)果，是理論結(jié)果的一個(gè)近似，而無限次迭代將使結(jié)果接近原本的直接計(jì)算，無論其原本是否穩(wěn)定、是否收斂.這也是我們應(yīng)用迭代法需要關(guān)注的客觀因素.

圖17 迭代化極結(jié)果Fig.17 Results of iterative RTP

5 結(jié) 論

本文主要從理論上推導(dǎo)了位場(chǎng)計(jì)算迭代法的通式，認(rèn)為迭代法收斂與否取決于映射函數(shù)的選取，并給出了具體轉(zhuǎn)換問題迭代法的收斂條件，分析了各種影響迭代收斂性的因素，例如：如果原本不穩(wěn)定計(jì)算的逆運(yùn)算是穩(wěn)定的，如迭代下延和低緯度赤道化極中，可以選擇簡(jiǎn)單的常數(shù)映射；但對(duì)于逆運(yùn)算不穩(wěn)定的情況，如向上延拓等（其逆運(yùn)算為不穩(wěn)定的向下延拓），那么對(duì)映射函數(shù)的選擇則要苛刻得多，必須采取變量映射才能有效收斂，所以，映射函數(shù)的合理選取是迭代法重磁處理轉(zhuǎn)換的核心.這就為迭代法在其他位場(chǎng)計(jì)算中的應(yīng)用奠定了理論基礎(chǔ).

映射函數(shù)的選取不僅對(duì)收斂性以及收斂速度起著絕對(duì)作用，并且對(duì)計(jì)算結(jié)果的誤差也有一定的影響.通常情況下，對(duì)于頻譜相位不發(fā)生改變的轉(zhuǎn)換，如相對(duì)高差不大的向下延拓等，映射函數(shù)的值為正值才能保證收斂；而對(duì)于低緯度化極等轉(zhuǎn)換，轉(zhuǎn)換前后數(shù)據(jù)頻譜的振幅和相位都發(fā)生了變化，映射函數(shù)的值為負(fù)值才能保證收斂；而對(duì)于曲化平等問題，轉(zhuǎn)換前后的頻譜的相位是否發(fā)生變化是不確定的，所以要根據(jù)具體情況來判斷，并且可能需要針對(duì)性設(shè)計(jì)變化的映射函數(shù).

圖18 速度系數(shù)對(duì)迭代化極誤差的影響Fig.18 The effects of velocity factor on the iterative RTP error

迭代法主要是用于解決一些位場(chǎng)計(jì)算中的不穩(wěn)定問題的，但從上面的理論分析以及模型數(shù)據(jù)計(jì)算也可以看出，迭代法雖然能夠通過試驗(yàn)不同參數(shù)獲得比較好的效果，但其結(jié)果的精度并不會(huì)隨著迭代次數(shù)的增加而保持增加，其本質(zhì)上相當(dāng)于還是一種濾波器（或其變形），其作用是對(duì)原有轉(zhuǎn)換算子放大的區(qū)域予以一定的濾波壓制；如果迭代次數(shù)過大，迭代法的結(jié)果將無限趨近于FFT理論計(jì)算結(jié)果，結(jié)果依然是不穩(wěn)定的；如果迭代次數(shù)過小則會(huì)壓制原始數(shù)據(jù)中的有效信號(hào)，使得結(jié)果誤差過大.所以針對(duì)位場(chǎng)處理轉(zhuǎn)換中一些不穩(wěn)定計(jì)算采用迭代法，并沒有從根本上解決計(jì)算的精度與不穩(wěn)定性之間的問題，這是我們應(yīng)用迭代法時(shí)需要認(rèn)識(shí)清楚的.

（References）

[1] CTPAXOBΒΗ，ДEBИЦЬIH B M.О ПPИBEДEHИИ HAБ？ ПОДEHHЪIX ЗHAЧEHИЙ ПОTEHЦИA？IЬHЬIXПO(jiān)JIEЙ КОДHOMУУPOBHIO.AKADEMIIA NAUK SSSR，IZUESTIIA，F(xiàn)IZIKA ZEML 1965，4：60－72.Strakhov A V，Devitsyn V N.The reduction of observed values of a potential field to values at a constant level.Proceedings of the Academy of Sciences.Physics of the Earth （in Russian），1965，4：60－72.

[2] 徐世浙.位場(chǎng)延拓的積分－迭代法.地球物理學(xué)報(bào)，2006，49（4）：1176－1182.Xu S Z.The integral－iteration method for continuation of potential fields.Chinese J.Geophys.（in Chinese）.2006，49（4）：1176－1182.

[3] 徐世浙.迭代法與FFT法位場(chǎng)向下延拓效果的比較.地球物理學(xué)報(bào)，2007，50（1）：285－289.Xu S Z.A comparison of effects between the iteration method and FFT for downward continuation of potential fields.Chinese J.Geophys.（in Chinese），2007，50（1）：285－289.

[4] 徐世浙，余海龍.位場(chǎng)曲化平的插值一迭代法.地球物理學(xué)報(bào)，2007，50（6）：1811－1815.Xu S Z，Yu H L.The interpolation—iteration method for potential field continuation from undulating surface to plane.Chinese J.Geophys.（in Chinese），2007，50（6）：1811－1815.

[5] 劉東甲，洪天求，賈志海等.位場(chǎng)向下延拓的波數(shù)域迭代法及其收斂性.地球物理學(xué)報(bào)，2009，52（6）：1599－1605.Liu D J，Hong T Q，Jia Z H，et al.Wave number domain iteration method for downward continuation of potential fields and its convergence.Chinese J.Geophys.（in Chinese），2009，52（6）：1599－1605.

[6] 張輝，陳龍偉，任治新等.位場(chǎng)向下延拓迭代法收斂性分析及穩(wěn)健向下延拓方法研究.地球物理學(xué)報(bào)，2009，52（4）：1107－1113.Zhang H，Chen L W，Ren Z X，et al.Analysis on convergence of iteration method for potential fields downward continuation and research on robust downward continuation method.Chinese J.Geophys.（in Chinese），2009，52（4）：1107－1113.

[7] 曾小牛，李夕海，韓紹卿等.位場(chǎng)向下延拓三種迭代方法之比較.地球物理學(xué)進(jìn)展，2011，26（3）：908－915.Zeng X N，Li X H，Han S Q，et al.A comparison of three iteration methods for downward continuation of potential fields.Progress in Geophysics.（in Chinese），2011，26（3）：908－915.

[8] 于波，翟國君，劉燕春等.噪聲對(duì)磁場(chǎng)向下延拓迭代法的計(jì)算誤差影響分析.地球物理學(xué)報(bào)，2009，52（8）：2182－2188.Yu B，Zhai G J，Liu Y C，et al.Analysis of noise effect on the calculation error of downward continuation with iteration method.Chinese J.Geophys.（in Chinese），2009，52（8）：2182－2188.

[9] 王彥國，張鳳旭，王祝文等.位場(chǎng)向下延拓的泰勒級(jí)數(shù)迭代法.石油地球物理勘探，2011，45（4）：657－662.Wang Y G，Zhang F X，Wang Z W，et al.Iteration method of Taylor series for downward continuation of potential fields.Oil Geophysical Prospecting （in Chinese），2011，45（4）：657－662.

[10] 駱遙.位場(chǎng)迭代法向下延拓的地球物理含義—以可下延異常逐次分離過程為例.地球物理學(xué)進(jìn)展，2011，26（4）：1197－1200.Luo Y.Geophysical implications of the iteration method for downward continuation of potential field：a case by computing regional anomaly can be continued downward.Progress in Geophysics.（in Chinese），2011，26（4）：1197－1200.

[11] 徐世浙，曹洛華，姚敬金.重力異常三維反演－視密度成像方法技術(shù)的應(yīng)用.物探與化探，2007，31（1）：25－28.Xu S Z，Cao L.H，Yao J.J.3D inversion of gravity anomaly：an application of the apparent density magery technology.Geophysical ＆ Geochemical Exploration （in Chinese），2007，31（1）：25－28.

[12] 徐世浙，余海龍，姚敬金等.新疆色爾特能地區(qū)視密度和視磁性反演.物探化探計(jì)算技術(shù)，2007，29（S1）：12－16.Xu S Z，Yu H L.Yao J.J，et al..The inversion of apparent density and apparent susceptibility in the seerteneng area，xin jiang.Computing Techniques for Geophysical and Geochemical Exploration （in Chinese），2007，29（S1）：12－16.

[13] 徐世浙，王瑞，周堅(jiān)鑫等，從航磁資料延拓出海平面的磁場(chǎng).海洋學(xué)報(bào)，2007，29（6）：53－57.Xu S Z，Wang R，Zhou J X，et al.To continue aeromagnetic data to the sea level.Acta Oceanologica Sinica，2007，29（6）：53－57.

[14] 徐世浙，王華軍，余海龍等，普光氣田重力異常的視密度反演.地球物理學(xué)報(bào)，2009，52（9）：2357－2363.Xu S Z，Wang H J，Yu H L，et al.Apparent density inversion of gravity anomalies of Puguang gas field.Chinese J.Geophys.（in Chinese），2009，52（9）：2357－2363.

[15] 楊金玉，徐世浙，余海龍等.視密度反演在東海及鄰區(qū)重力異常解釋中的應(yīng)用.地球物理學(xué)報(bào)，2008，51（6）：1909－1916.Yang J Y，Xu S Z，Yu H L，et al.Application of apparent density inversion method in the East China Sea and its adjacent area.Chinese J.Geophys.（in Chinese），2008，51（6）：1909－1916.

[16] 于波，劉雁春，邊剛等，海洋磁力測(cè)量垂直空間改正研究.武漢大學(xué)學(xué)報(bào).2008，33（9）：926－929.Yu B，Liu Y C.Bian G，et al.Correction of Vertical Free－Air in Marine Magnetic Surveying.Geomatics and Information，2008，33（9）：926－929.

[17] 于波，翟國君，劉雁春等，利用航磁數(shù)據(jù)向下延拓得到海平面磁場(chǎng)的方法.測(cè)繪學(xué)報(bào).2009，38（3）：202－209.Yu B，Zhai G J，Liu Y C，et al.The downward continuation method of aeromagnetic data to the sea level.Acta Geodaetica et Cartographica Sinica.2009，38（3）：202－209.

[18] 王小兵，張金生，喬玉坤等.頻率域延拓在地磁匹配中的應(yīng)用研究了.現(xiàn)代防御技術(shù)，2009，37（6）：87－90.Wang X B，Zhang J S，Qiao Y K，et al.Application of the frequency domain continuation in the geomagnetism matching guidance system.Modern Defence Technology，2009，37（6）：87－90.

[19] 肖鋒，李海俠，余海龍等.確定水下磁性體位置的一個(gè)實(shí)例.物探與化探，2009，33（1）：105－108.Xiao F，Li H X，Yu H L，et al.The Determination of the position of an under water magnetic object：a case study.Geophysical ＆ Geochemical Exploration （in Chinese），2009，33（1）：105－108.

[20] 李成立，崔瑞華，劉益中等.盆地基底巖性的綜合地球物理預(yù)測(cè)方法—以松遼盆地濱北地區(qū)基底巖性預(yù)測(cè)為例.地球物理學(xué)報(bào)，2011，54（2）：491－498.Li C L，Cun R H，Liu Y Z，et al.Comprehensive geophysical prediction method of basement lithology－Example of Binbei area，Songliao basin.Chinese J.Geophys.（in Chinese），2011，54（2）：491－498.

[21] 余海龍，徐世浙，李海俠等，曲面上航磁異常與梯度分量的轉(zhuǎn)換方法.浙江大學(xué)學(xué)報(bào)（工學(xué)版），2011，45（2）：397－404.Yu H L，Xu S Z，Li H X，et al.Transformation of total magnetic field anomaly and gradients on irregular surface.Journal of Zhejiang University （Engineering Science），2011，45（2）：397－404.

[22] Keating P，Zerbo L.An improved technique for reduction to the pole at low latitudes.Geophysics，1996，61（1）：131－137.

[23] Li X. Magnetic reduction－to－the－pole at low latitudes：Observations and considerations.The Leading Edge，2008，27（8）：900－1002.

[24] 駱遙，薛典軍.磁赤道處化極方法.地球物理學(xué)報(bào)，2010，53（12）：2998－3004.Luo Y，Xue D J.Reduction to the pole at the geomagnetic equator.Chinese J.Geophys.（in Chinese），2010，53（12）：2998－3004.

[25] 方迎堯，張培琴，劉浩軍.低磁緯度地區(qū)ΔT異常解釋的途徑與方法.物探與化探，2006，30（1）：48－54.Fang Y Y，Zhang P Q，Liu H J.Approaches to the interpretation of magnetic Δ Tanomalies in the low magnetic latitude area.Geophysical ＆ Geochemical Exploration （in Chinese），2006，30（1）：48－54.

Research on iteration method using in potential field transformations

YAO Chang－Li，LI Hong－Wei，ZHENG Yuan－Man，MENG Xiao－Hong，ZHANG Yu－Wen
1 Key Laboratory of Geo－detection （China University of Geosciences，Beijing），Ministry of Education，Beijing100083，China State Key Laboratory of Geological Processes and Mineral Resources，Beijing100083，China
2 School of Geophysics and Information Technology，China University of Geosciences，Beijing100083，China

Data processing and transformations play an important role in the interpretations of gravity and magnetic data，but some calculations，such as downward continuation，reduction to the pole of magnetic anomalies and so on，are unstable at times，because their transformation factors generally have a significant amplifying effect，so that the results obtained directly from FFT are unstable.Many studies are focused on increasing the stability and the effectiveness of the calculations，one of which is iterative method which has achieved good results and has recently

wide attention.But the study of iterative problem is not deep enough，andobviously there exists a lack of awareness of its shortcomings fully.In this paper，the general iterative formula is deduced based on an analytical study of the iterative method，and convergence effect of various factors on the iteration are analyzed.The results show that the main factor which determinants whether the iteration methods will convergence to the FFT theoretical results is how to choose the mapping function from the raw data to the target data，and in the case of a suitable mapping function chosen，the number of iteration in calculation is not just a cost，but is a key factor in determining the outcome good or bad.But if the FFT theoretical calculation of the transformation is unstable，the iteration will converge to an unstable result.So the iteration method can＇t fundamentally resolve the instability of an unstable.

Potential field transformation，Iteration method，F(xiàn)FT transformation，Convergence，Condition of convergence

10.6038/j.issn.0001－5733.2012.06.028

P631

2011－11－16，2012－05－30收修定稿

國家公益性項(xiàng)目（201011039）、國家高技術(shù)研究發(fā)展計(jì)劃（2007AA06Z134）和高等學(xué)校學(xué)科創(chuàng)新引智計(jì)劃（B07011）聯(lián)合資助.

姚長利，男，1965年生，教授，主要從事重磁勘探方法技術(shù)研究.E－mail：clyao＠cugb.edu.cn

姚長利，李宏偉，鄭元滿等：重磁位場(chǎng)轉(zhuǎn)換計(jì)算中迭代法的綜合分析與研究.地球物理學(xué)報(bào)，2012，55（6）：2062－2078，

10.6038/j.issn.0001－5733.2012.06.028.

Yao C L，Li H W，Zhen Y M，et al.Research on iteration method using in potential field transformations.Chinese J.Geophys.（in Chinese），2012，55（6）：2062－2078，doi：10.6038/j.issn.0001－5733.2012.06.028.

（本文編輯 劉少華）