孫 文，吳曉平，王慶賓，劉曉剛，王 凱

1．信息工程大學(xué)地理空間信息學(xué)院，河南 鄭州 450001；2．地理信息工程國家重點實驗室，陜西 西安 710054；3．西安測繪研究所，陜西 西安 710054；4．武漢大學(xué)地球空間環(huán)境與大地測量教育部重點實驗室，湖北 武漢 430079；5．總裝備部工程設(shè)計研究總院，北京 100028

航空重力數(shù)據(jù)向下延拓的波數(shù)域迭代Tikhonov正則化方法

孫 文1，吳曉平1，王慶賓1，劉曉剛2，3，4，王 凱5

1．信息工程大學(xué)地理空間信息學(xué)院，河南 鄭州 450001；2．地理信息工程國家重點實驗室，陜西 西安 710054；3．西安測繪研究所，陜西 西安 710054；4．武漢大學(xué)地球空間環(huán)境與大地測量教育部重點實驗室，湖北 武漢 430079；5．總裝備部工程設(shè)計研究總院，北京 100028

探討了Tikhonov正則化方法以及迭代Tikhonov正則化方法在航空重力數(shù)據(jù)向下延拓中的應(yīng)用，指出Tikhonov正則化方法在利用GCV法或L曲線法選取正則化參數(shù)時存在的不足以及空間域迭代法迭代不收斂的問題。引入波數(shù)域Tikhonov迭代正則化法，有效解決了上述方法存在的問題；針對迭代次數(shù)的確定問題，提出了基于重力異常差異熵的迭代次數(shù)確定方法。算例結(jié)果表明，波數(shù)域迭代法迭代過程收斂，且計算精度高、速度快，值得廣泛應(yīng)用于航空重力數(shù)據(jù)的向下延拓。

航空重力數(shù)據(jù)；向下延拓；迭代法；波數(shù)域；正則化參數(shù)；Tikhonov正則化；重力異常差異熵

1 引 言

隨著科技的進步與發(fā)展，重力數(shù)據(jù)的獲取不再局限于傳統(tǒng)的地面重力測量，利用快速高效的航空重力或衛(wèi)星重力測量方法測量載體飛行高度面的重力異常，然后通過數(shù)學(xué)手段求得地面重力異常，成為獲取大面積區(qū)域尤其是地面重力觀測難以進行區(qū)域重力數(shù)據(jù)的重要方法，由此產(chǎn)生了將空中重力異常向下延拓至地面高度的需求。

眾所周知，向下延拓屬于地球物理反問題，一般難以獲得穩(wěn)定解。對此，國內(nèi)外學(xué)者展開了深入研究，獲得了眾多具有理論意義或?qū)嶋H應(yīng)用價值的成果。目前，航空重力數(shù)據(jù)向下延拓主要有以下幾種方法：正則化法、配置法、迭代法、“倒錐法”、直接代表法、點質(zhì)量法，等。嚴(yán)格來說，迭代法屬于正則化法的一種，但本文將其單獨作為一種方法，因為迭代法與正則化法在正則化參數(shù)的選擇方面有較大不同。

文獻［1—2］探討了正則化方法在向下延拓中的應(yīng)用，比較分析了正則化參數(shù)選取的方法以及效果。文獻［3］將最小二乘配置法應(yīng)用于重力數(shù)據(jù)的向下延拓。文獻［4］提出了向下延拓重力數(shù)據(jù)的“倒錐法”，能夠完成高精度的向下延拓計算。文獻［5—6］分析了地面與航空重力數(shù)據(jù)的頻域特征，指出隨著高度的增加，空中重力異常的頻譜衰減劇烈，并據(jù)此提出了直接代表法，以空中重力異常作為地面重力異常的低頻項，利用DEM數(shù)據(jù)進行向下延拓的高頻改正，得到了較為滿意的結(jié)果。文獻［7］比較分析了多種向下延拓方法的精度，說明了正則化方法的優(yōu)勢所在。文獻［8］詳細(xì)研究了“移去—恢復(fù)”方法在向下延拓正則化方法中的應(yīng)用，指出該方法能夠提高正則化方法向下延拓的精度。文獻［9］提出基于多尺度邊緣約束的重力場信號向下延拓方法，進一步豐富了向下延拓的理論與方法。文獻［10］比較了Tikhonov正則化、嶺估計和廣義嶺估計應(yīng)用于向下延拓的效果，指出廣義嶺估計相比其余兩種延拓方法在精度、穩(wěn)定性和抗差性方面均具有更好效果。文獻［11］提出基于信噪比的正則化方法并將其用于航空重力向下延拓。文獻［12］詳細(xì)分析了地形對向下延拓的影響，文獻［13］從Poisson積分公式出發(fā)，研究了向下延拓的不同機制，指出點延拓與面延拓應(yīng)采用不同的平滑因子，從而減小離散誤差的影響。文獻［14］比較了3種迭代正則化方法，并以球體重力場模型為基礎(chǔ)對實際延拓效果進行了分析。

這些方法的提出與應(yīng)用，一定程度上解決了重力場向下延拓的精度問題，但仍然存在一些不足。以正則化方法為例，該方法需要準(zhǔn)確選取正則化參數(shù)并需要進行大型矩陣的奇異值分解計算，計算速度非常慢，無法應(yīng)用于大面積區(qū)域；配置法由于協(xié)方差函數(shù)確定的巨大困難導(dǎo)致該方法延拓精度有限；“倒錐法”需要在延拓區(qū)域有地面實測重力數(shù)據(jù)的支持，通常難以滿足。

本文引入迭代Tikhonov正則化方法，利用其空間域和波數(shù)域兩類迭代求解航空重力異常的向下延拓，并與Tikhonov正則化方法作比較，得出一些有益的結(jié)論，為后續(xù)研究提供一定的參考。

2 向下延拓的數(shù)學(xué)模型

重力異常的向上延拓是關(guān)于球外部Dirichlet

問題的求解過程，其結(jié)果由球近似下的Poisson

積分公式給出［15］式中，R為地球平均半徑；Δgh、Δg分別為空中與地面重力異常；r＝R＋h；h為向上延拓高度；（φ，λ）和（φ′，λ′）分別為空中計算點與地面流動點的坐標(biāo)；為兩點間距離，φ為兩點間球面角距，可由式（2）計算

在不影響精度的條件下，式（1）可平面近似為

對于式（1），已知Δgh求解Δg的問題即為重力異常的向下延拓。通常，反問題都是不適定問題。此時式（1）為第一類Fredholm積分方程，當(dāng)Δgh存在誤差時，直接求逆結(jié)果極不穩(wěn)定，從而限制了向下延拓重力異常的精度。

3 向下延拓的正則化方法

3．1 Tikhonov正則化方法及其譜分解

在求解地球物理反問題的方法中，正則化方法作為一種重要理論在其中起著不可替代的作用，而其中尤以Tikhonov提出的正則化方法［16］應(yīng)用最為廣泛，其基本思想為［1］：求解參數(shù)x的估值，使得

式中，α＞0為正則化參數(shù)。由式（4）的約束條件，可以得到x的正則化解表達(dá)式為

對系數(shù)陣A進行奇異值分解得

式中，U、V分別是A的左右酉矩陣，U＝［ui］，V＝［vi］，i＝1，2，．．．，n；Λ為對角陣，其對角線元素由A矩陣的奇異值λi按照遞減的次序排列而成。則式（5）的譜分解形式為

由式（7）可以看出，正則化參數(shù)α調(diào)整奇異值單調(diào)趨于0的趨勢，改善了最終解的穩(wěn)定性。利用Tikhonov正則化方法求解不適定問題的關(guān)鍵在于正則化參數(shù)α的選取，主要方法有廣義交叉驗證（GCV）法［17］、L－曲線法［18］、圖解－數(shù)值迭代法［1］等。

GCV法基本原理是：選取參數(shù)α，使得GCV函數(shù)取得最小值，其中tr（）表示矩陣的跡，Qα＝A（αI＋ATA）－1AT。顯然，GCV法求取正則化參數(shù)的過程是迭代計算的過程。

L－曲線法的原理是：選取適當(dāng)?shù)膮?shù)α，以平衡lg Axα－L2與lg xα2的值。實際操作時，選取L－曲線（lg Axα－L2，lg xα2）的拐點所對應(yīng)的參數(shù)，此α值即為所求最佳正則化參數(shù)。

3．2 空間域Tikhonov迭代法

空間域Tikhonov迭代法（以下簡稱空間域迭代法）可表示為［19－20］

式（10）即為空間域迭代法的實用公式。給定最大迭代次數(shù)N或很小的正數(shù)ε，當(dāng)滿足＜ε或n＝N 時停止迭代?？梢宰C明［19］，對于給定的正數(shù)α，當(dāng)n→∞時，xnα→x。

3．3 波數(shù)域Tikhonov迭代法

由于向上延拓過程穩(wěn)定可靠，故可直接采用快速傅里葉算法（FFT）計算。式（3）右端是卷積形式，據(jù)此，文獻［21］導(dǎo)出了向上延拓算子為Φ（ω）＝e－ωh，稱為原始延拓算子，則式（3）的波數(shù)域表達(dá)式為［22］

式（11）即為向下延拓的FFT法。

波數(shù)域Tikhonov迭代法（以下簡稱波數(shù)域迭代法）即在空間域迭代過程中引入FFT方法，將其轉(zhuǎn)化到波數(shù)域進行迭代求解，然后通過傅里葉逆變換求得最終的延拓效果。對式（9）兩邊作傅里葉變換可得根據(jù)文獻［14］，由數(shù)學(xué)歸納法可以最終推得波數(shù)域迭代法的向下延拓算子為，n為迭代次數(shù)，α為給定正數(shù)。

則波數(shù)域迭代法的最終表達(dá)式為

式中，F(xiàn)、F－1分別表示傅里葉變換和逆變換過程。

4 算例與分析

4．1 數(shù)據(jù)

本文將覆蓋范圍為1．5°×1．5°、分辨率為2′× 2′的某區(qū)域?qū)崪y地面格網(wǎng)平均重力異常用式（1）延拓至3500 m高度，積分半徑為30′，以此作為向下延拓的初始數(shù)據(jù)，以原始地面重力異常作為向下延拓結(jié)果的檢核數(shù)據(jù)。地面（h＝0 m）重力異常以及h＝3500 m高度處航空重力異常數(shù)據(jù)統(tǒng)計情況列于表1。

表1 地面與航空重力異常數(shù)據(jù)統(tǒng)計Tab．1 Statistics of ground and airborne gravity anomaly m Gal

為避免計算過程中邊界效應(yīng)對最終結(jié)果的影響，取30′為積分半徑，計算全部區(qū)域向下延拓的結(jié)果，但在檢核延拓效果時取位于該地區(qū)中央的1°×1°區(qū)域，因該區(qū)域數(shù)據(jù)受邊界效應(yīng)的影響較??；在延拓過程中，“移去—恢復(fù)”理論作為向下延拓的重要方法也同時加以應(yīng)用，采用的參考模型為EGM2008全球重力場模型［23］的前360階次。

4．2 結(jié)果與分析

4．2．1 病態(tài)性分析及最小二乘解

如上所述，向下延拓問題的法矩陣一般呈現(xiàn)出病態(tài)性特征，其條件數(shù)往往遠(yuǎn)超病態(tài)臨界值1000。經(jīng)計算，由該區(qū)域數(shù)據(jù)組成的法矩陣條件數(shù)高達(dá)13 328，屬高度病態(tài)矩陣，若直接采用最小二乘方法（LS）求解，將造成結(jié)果的極不穩(wěn)定。系數(shù)陣的奇異值分布情況如圖1所示。

圖1 系數(shù)陣奇異值分布圖Fig．1 Singular values of coefficient matrix

由圖1可以看出，奇異值單調(diào)趨于0，這是造成最小二乘直接求逆計算結(jié)果不穩(wěn)定的主要原因。直接求逆的結(jié)果與地面實際值的差值如圖2所示（經(jīng)緯度坐標(biāo)已經(jīng)過處理，非實際坐標(biāo)，下同）。

圖2 最小二乘法差值分布Fig．2 Error distribution of LS method

從圖2可以看出，由于法矩陣嚴(yán)重病態(tài)，最小二乘直接求逆的向下延拓結(jié)果極不穩(wěn)定，無法在實際數(shù)據(jù)處理過程中采用。

4．2．2 Tikhonov正則化法

Tikhonov正則化方法能夠調(diào)整奇異值趨于0的程度，從而改善最小二乘計算結(jié)果。基于GCV法和L－曲線法選取正則化參數(shù)，選取范圍定為系數(shù)陣最大奇異值與最小奇異值之間［2］。GCV函數(shù)值隨正則化參數(shù)α變化趨勢以及L曲線如圖3（a）和圖3（b）所示。

圖3 GCV函數(shù)及L－曲線隨α變化趨勢示意圖Fig．3 GCV function values and L－curve with differentα

由圖3可以看出，該區(qū)域GCV函數(shù)并不收斂，且L－曲線不存在拐點，只能選擇曲率半徑最小處作為L－曲線的解。這是Tikhonov正則化方法存在的主要問題。取最小GCV值對應(yīng)的α參數(shù)值（L－曲線曲率半徑最小處的α值與其相同），則α＝0．011 7，對應(yīng)系數(shù)陣最小奇異值。利用該α值進行Tikhonov正則化向下延拓，其最終結(jié)果與實際地面重力異常的差值分布如圖4所示。

比較圖2和圖4可以看出，Tikhonov正則化方法能夠在一定程度上消除法矩陣病態(tài)的影響，但是在局部區(qū)域的延拓結(jié)果仍然存在較大誤差，其原因可能是由于GCV函數(shù)不收斂，L－曲線拐點不存在，導(dǎo)致無法確定最優(yōu)正則化參數(shù)，從而影響了最終的延拓效果。

圖4 Tikhonov正則化法差值分布Fig．4 Error distribution of Tikhonov regularization method

4．2．3 空間域迭代法

為了分析正則化參數(shù)對空間域迭代法延拓精度的影響，取不同的正則化參數(shù)，分別利用式（10）進行解算，其均方根值隨迭代次數(shù)及正則化參數(shù)α的變化趨勢如圖5所示。

圖5 空間域迭代Tikhonov正則化法精度變化趨勢Fig．5 Space domain iterative Tikhonov regularization precision change with iterated times andα

由圖5可以看出，空間域迭代法的延拓精度與迭代次數(shù)相關(guān)，迭代次數(shù)越多，精度越高；另外，隨著α的減小，前幾次迭代精度的變化越為明顯；無論α的取值如何，迭代10次基本都能達(dá)到該方法所能達(dá)到的最高精度。

但圖5也說明了該方法的重要問題所在：即無論α如何取值，其精度變化曲線不收斂，這說明雖然理論證明式（10）收斂，仍然需要考慮其適用性問題。圖6是α＝1、迭代次數(shù)為20時延拓結(jié)果與實際值的誤差分布圖，與圖4比較可以看出，兩種方法的結(jié)果較為接近，僅有細(xì)微的差別。這說明，若不考慮收斂性問題，空間域迭代法由于對正則化參數(shù)不敏感，從而在應(yīng)用時無需精確選取正則化參數(shù)，這是空間域迭代法相比Tikhonov正則化法的優(yōu)勢所在，此時起正則化作用的是迭代次數(shù)。

圖6 空間域迭代法差值分布Fig．6 Error distribution of space domain iterative method

4．2．4 波數(shù)域迭代法

首先分析延拓算子在不同的迭代次數(shù)和正則化參數(shù)下的頻率響應(yīng)。為便于說明，頻率響應(yīng)分析按照一維情況處理，二維情況與一維相似，注意此處h＝3500 m。圖7（a）、圖7（b）分別是迭代次數(shù)為5次和50次時，不同的正則化參數(shù)所對應(yīng)延拓算子的頻率響應(yīng)大小；圖8（a）、圖8（b）分別是正則化參數(shù)取0．05和0．5時，不同的迭代次數(shù)所對應(yīng)延拓算子的頻率響應(yīng)大小。

圖7 不同迭代次數(shù)的頻率響應(yīng)Fig．7 Frequency response with different iterated times

圖8 不同正則化參數(shù)的頻率響應(yīng)Fig．8 Frequency response with different regularization parameters

比較圖7、圖8可以發(fā)現(xiàn)，波數(shù)域迭代法的延拓算子是帶通濾波器，其最大頻率響應(yīng)隨著迭代次數(shù)的增加而增大，隨著正則化參數(shù)的增大而減小。當(dāng)?shù)螖?shù)過多時，其延拓精度會隨著迭代的進行而變得極不穩(wěn)定。從上述分析也可以看出，波數(shù)域迭代法與空間域迭代法類似，迭代次數(shù)在計算過程中起著正則化作用。圖9是波數(shù)域迭代法的向下延拓精度隨迭代次數(shù)及α的變化趨勢圖，從中可以看出，無論α如何取值，隨著迭代的進行，其結(jié)果先收斂后發(fā)散，且隨著α的增大，收斂所需的迭代次數(shù)也在增加，例如α分別取0．05、0．5和1時，收斂所需的迭代次數(shù)分別為3次、38次和69次。

圖9 波數(shù)域迭代法精度隨迭代次數(shù)及α的變化趨勢Fig．9 Wave number domain iterative method precision change with iterated times andα

取α＝1、迭代69次作為波數(shù)域迭代法的計算參數(shù)，其向下延拓結(jié)果與實際值的差值分布如圖10所示。

圖10 波數(shù)域迭代法差值分布Fig．10 Error distribution of wave number domain iterative method

4．2．5 比較與分析

表2列出了4種方法精度及運算速度的比較統(tǒng)計情況，所用計算機為：Windows XP SP3操作系統(tǒng)，Matlab R2011b試驗平臺，主頻3．2 GHz，1 GB內(nèi)存。其中，統(tǒng)計消耗時間時，沒有考慮系數(shù)矩陣構(gòu)建所消耗的時間；兩種迭代法的迭代次數(shù)均為100次，α均取為1，最終結(jié)果所對應(yīng)的迭代次數(shù)分別為20次與69次；Tikhonov正則化方法所消耗時間包括系數(shù)矩陣奇異值分解、GCV函數(shù)計算所消耗時間，α取值0．011 7。

表2 4種方法計算精度與速度比較Tab．2 Comparation of four methods

比較圖2、圖4、圖6、圖10以及表2可以明顯看出，在向下延拓精度方面，除去最小二乘法，波數(shù)域迭代法的精度最高，空間域迭代法次之，Tikhonov正則化法最低，其原因主要是由于后兩種方法存在不收斂的現(xiàn)象，即對于Tikhonov正則化法，GCV函數(shù)不收斂，對于空間域迭代法，迭代過程不收斂。故兩種方法雖然能夠完成向下延拓的解算，但其精度受到較大限制。波數(shù)域迭代法的精度隨著迭代次數(shù)的增加明顯存在先收斂后發(fā)散的過程，其最終結(jié)果相比另外兩種方法更加可靠。在消耗時間方面，波數(shù)域迭代法也有較大優(yōu)勢，空間域迭代法次之，Tikhonov正則化法最慢。其原因是由于迭代法需要進行矩陣求逆運算，而Tikhonov正則化法需要計算矩陣的奇異值分解，兩種運算都非常消耗時間。

值得說明的是，若考慮計算機內(nèi)存的問題，Tikhonov正則化法由于需要計算矩陣的奇異值分解，對內(nèi)存的需求量巨大，例如上述試驗所用計算機在進行5000階以上矩陣的奇異值分解運算時會出現(xiàn)內(nèi)存不夠的問題，而迭代法消耗內(nèi)存量相比Tikhonov正則化法較少，這一點在進行較大區(qū)域延拓計算時應(yīng)加以考慮。

4．2．6 基于差異熵確定迭代次數(shù)

由圖9可知，迭代法的關(guān)鍵在于確定準(zhǔn)確的迭代終止條件，一般取廣義交叉驗證方法，但上文的研究已經(jīng)表明，該方法不收斂，故此處無法應(yīng)用。當(dāng)?shù)孛嬷亓Ξ惓Ｎ粗獣r，每次迭代的結(jié)果精度無法評估，更加難以選擇合適的迭代次數(shù)。為了準(zhǔn)確確定不同的正則化參數(shù)對應(yīng)的迭代次數(shù)，受信息熵概念的啟發(fā)，提出一種基于差異熵的迭代次數(shù)確定方法。

地球物理學(xué)中，區(qū)域重力異常熵定義為［24］式中，H（Δg）表示重力異常熵；p（Δgi，j）表示Δgi，j的概率。p（Δgi，j）理論上應(yīng)使用大量數(shù)據(jù)經(jīng)統(tǒng)計得出，但在區(qū)域范圍內(nèi)應(yīng)用時也可直接計算

差異熵更能反映區(qū)域內(nèi)重力異常的變化［24］，其定義與式（15）形式相同，所不同的是p（Δgi，j）。以表示該區(qū)域重力異常的均值，定義Ci，j＝

利用差異熵確定迭代次數(shù)的基本思想為：熵值體現(xiàn)了系統(tǒng)的不確定程度，較大的向下延拓誤差體現(xiàn)為較大的熵值，故在每次迭代過程中，計算每次延拓結(jié)果的差異熵，最小的差異熵對應(yīng)最高精度的延拓結(jié)果。該方法的主要誤差源是p（Δgi，j）的計算，包含兩個方面的因素：①無論采取何種手段，向下延拓的結(jié)果誤差不可避免，導(dǎo)致p（Δgi，j）的計算誤差；②利用部分?jǐn)?shù)據(jù)樣本代替總體概率統(tǒng)計所造成的誤差，這應(yīng)是主要的誤差源所在。

為了驗證所提方法的效果，在上述計算過程中加入差異熵的計算，根據(jù)最小差異熵獲取不同的正則化參數(shù)所對應(yīng)的迭代次數(shù)，并與圖9中實際最高精度對應(yīng)的迭代次數(shù)作比較，其差值結(jié)果如圖11所示。

圖11 基于重力異常差異熵確定迭代次數(shù)的精度Fig．11 Precision of iterated times determination based on gravity anomaly variance entropy

由圖11可以看出，基于差異熵的迭代次數(shù)確定方法是合理且有效的，當(dāng)正則化參數(shù)小于1時，該方法確定的迭代次數(shù)與實際所需迭代次數(shù)差值在3次以內(nèi)；若進一步縮小正則化參數(shù)的取值，則確定迭代次數(shù)的準(zhǔn)確率更高，這充分驗證了該方法的有效性。同時，該方法的準(zhǔn)確性與正則化參數(shù)的大小存在一定關(guān)系，正則化參數(shù)選取較大時，其準(zhǔn)確率有所降低，原因尚不清楚，故在實際應(yīng)用時，建議選取較小的正則化參數(shù)（小于0．5）。

5 結(jié) 論

本文針對重力異常向下延拓的病態(tài)性問題，引入迭代Tikhonov正則化方法，分別介紹了空間域和波數(shù)域兩種迭代法的理論計算方法，并進行了3種方法的實際解算和比較，指出Tikhonov正則化方法和空間域迭代法存在的問題。通過以上計算和分析，可以得到如下幾點結(jié)論：

（1）在進行航空重力數(shù)據(jù)向下延拓的求解過程中，最小二乘方法所得結(jié)果極不穩(wěn)定，采用不同的正則化方法能夠不同程度地改善向下延拓的精度和穩(wěn)定性。

（2）對于Tikhonov正則化方法，使用GCV方法求解正則化參數(shù)時，GCV函數(shù)會不收斂；利用L－曲線法時會出現(xiàn)L曲線拐點無法確定的情況，此時得到的正則化參數(shù)雖然能夠改善求解精度，但是限制了精度的進一步提高。

（3）對于空間域迭代法，雖然理論證明該方法最終收斂，但是在實際應(yīng)用時，對于固定的正則化參數(shù)，其均方根曲線存在不收斂的情況，其原因仍待進一步研究，在使用時應(yīng)加以注意。

（4）波數(shù)域迭代法將向下延拓問題轉(zhuǎn)化為帶通濾波器，從而實現(xiàn)快速、穩(wěn)定的重力異常向下延拓，該方法計算速度快，穩(wěn)定性好，精度高，通過迭代次數(shù)實現(xiàn)正則化，對正則化參數(shù)不敏感，故無需進行正則化參數(shù)的選取，且不會出現(xiàn)類似空間域迭代法或者Tikhonov正則化方法存在的不收斂問題，值得廣泛用于航空重力異常向下延拓的解算。

（5）針對迭代法迭代次數(shù)的確定問題，提出了基于重力異常差異熵的迭代次數(shù)確定方法，實際計算結(jié)果表明，該方法能夠準(zhǔn)確確定迭代次數(shù)，特別是當(dāng)正則化參數(shù)取值較小時，效果更佳。

［1］ WANG Xingtao，SHI Pan，ZHU Feizhou．Regularization Methods and Spectral Decomposition for the Downward Continuation of Airborne Gravity Data［J］．Acta Geodaetica et Cartographica Sinica，2004，33（1）：33－38．（王興濤，石磐，朱非洲．航空重力測量數(shù)據(jù)向下延拓的正則化算法及其譜分解［J］．測繪學(xué)報，2004，33（1）：33－38．）

［2］ DENG Kailiang，HUANG Motao，BAO Jingyang，et al．Tikhonov Two－parameter Regulation Algorithm in Downward Continuation of Airborne Gravity Data［J］．Acta Geodaetica et Cartographica Sinica，2011，40（6）：690－696．（鄧凱亮，黃謨濤，暴景陽，等．向下延拓航空重力數(shù)據(jù)的Tikhonov雙參數(shù)正則化法［J］．測繪學(xué)報，2011，40（6）：690－696．）

［3］ KELLER W，HIRSCH M．Downward Continuation versus Free Air Reduction in Airborne Gravimetry［C］∥Geodesy and Physics of the Earth：Geodetic Contributions to Geodynamics．Berlin：Springer－Verlag，1992：266－272．

［4］ YANG Qiangwen，WU Xiaoping．“Qusai－cone Method”，a New Way of Downward Continuation of Airborne Gravimetry Data［J］．Journal of the PLA Institute of Surveying and Mapping，1998，15（3）：169－172．（楊強文，吳曉平．航空重力數(shù)據(jù)向下延拓的倒錐法［J］．解放軍測繪學(xué)院學(xué)報，1998，15（3）：169－172．）

［5］ SHI Pan，WANG Xingtao．Determination of the Terrain Surface Gravity Field Using Airborne Gravimetry and DEM［J］．Acta Geodaetica et Cartographica Sinica，1997，26（2）：117－121．（石磐，王興濤．利用航空重力測量和DEM確定地面重力場［J］．測繪學(xué)報，1997，26（2）：117－121．）

［6］ SHI Pan，WANG Xingtao．The Frequency Domain Analysis of the Determination of Terrestrial Mean Gravity Anomaly Using Airborne Gravimetry［J］．Acta Geodaetica et Cartographica Sinica，1995，24（4）：301－308．（石磐，王興濤．空中測量地面平均重力異常的頻域分析［J］．測繪學(xué)報，1995，24（4）：301－308．）

［7］ WANG Xingtao，XIA Zheren，SHI Pan，et al．A Comparison of Different Downward Continuation Methods for Airborne Gravity Data［J］．Chinese Journal of Geophysics，2004，47（6）：1017－1021．（王興濤，夏哲仁，石磐，等．航空重力測量數(shù)據(jù)向下延拓方法比較［J］．地球物理學(xué)報，2004，47（6）：1017－1021．）

［8］ HAO Yanling，CHENG Yi，SUN Feng，et a1．Simulation Research on Tikhonov Regularization Algorithm in Downward Continuation［J］．Chinese Journal of Scientific Instrument，2008，29（8）：2190－2194．（郝燕玲，成怡，孫楓，等．Tikhonov正則化向下延拓算法仿真實驗研究［J］．儀器儀表學(xué)報，2008，29（8）：2190－2194．）

［9］ NING Jinsheng，WANG Haihong，LUO Zhicai．Downward Continuation of Gravity Signals Based on the Multiscale Edge Constraint［J］．Chinese Journal of Geophysics，2005，48（1）：63－68．（寧津生，汪海洪，羅志才．基于多尺度邊緣約束的重力場信號的向下延拓［J］．地球物理學(xué)報，2005，48（1）：63－68．）

［10］ JIANG Tao，LI Jiancheng，WANG Zhengtao，et al．Solution of Ill－posed Problem in Downward Continuation of Airborne Gravity［J］．Acta Geodaetica et Cartographica Sinica，2011，40（6）：684－689．（蔣濤，李建成，王正濤，等．航空重力向下延拓病態(tài)問題的求解［J］．測繪學(xué)報，2011，40（6）：684－689．）

［11］ GU Yongwei，GUI Qingming．Study of Regularization Based on Signal－to－noise Index in Airborne Gravity Downward to the Earth Surface［J］．Acta Geodaetica et Cartographica Sinica，2010，39（5）：458－464．（顧勇為，歸慶明．航空重力測量數(shù)據(jù)向下延拓基于信噪比的正則化方法的研究［J］．測繪學(xué)報，2010，39（5）：458－464．）

［12］ WANG Y M．The Effect of Topography on the Determination of the Geoid Using Analytical Downward Continuation［J］．Bulletin Géodésique，1990，64：231－246．

［13］ GOLI M，NAJAFI－ALAMDARI M，VANICEK P．Numerical Behaviour of the Downward Continuation of GravityAnomalies［J］．Studia Geophysica et Geodaetica，2011，55：191－202．

［14］ ZENG Xiaoniu，LI Xihai，H AN Shaoqing，et al．A Comparison of Three Iteration Methods for Downward Continuation of Potential Fields［J］．Progress in Geophysics，2011，26（3）：908－915．（曾小牛，李夕海，韓紹卿，等．位場向下延拓三種迭代方法之比較［J］．地球物理學(xué)進展，2011，26（3）：908－915．）

［15］ HEISKANEN W A，MORITZ H．Physical Geodesy［M］．Beijing：Surveying and Mapping Press，1979．（海斯卡涅W A，莫里斯H．物理大地測量學(xué)［M］．北京：測繪出版社，1979．）

［16］ TIKHONOV A N．Solution of Ill－posed Problem［M］．New York：John Wiley ＆Arsenin，1977．

［17］ KILMER M E，O’LEARY D P．Choosing Regularization Parameters in Iterative Methods for Ill－posed Problems［J］．SIAM Journal on Matrix Analysis and Applications，2001，22（4）：1204－1221．

［18］ HANSEN P C，O’LEARY D P．The Use of the L－curve in the Regularization of Discrete Ill－posed Problems［J］．SIAM Journal on Scientific Computing，1993，14（6）：1487－1503．

［19］ ZHANG Luyin，ZHANG Yuhai，QIAN Kunming．On the Iterated Tikhonov Regularization for Ill－posed Problems［J］．Journal of Shandong University：Natural Science，2011，46（4）：29－33．（張路寅，張玉海，錢坤明．關(guān)于不適定問題的迭代Tikhonov正則化方法［J］．山東大學(xué)學(xué)報：理學(xué)版，2011，46（4）：29－33．）

［20］ KALTENBACHER B，NEUBAUER A，SCHERZER O．Iterative Regularization Methods for Nonlinear Ill－posed Problems［M］．New York：Walter de Gruyter，2008．

［21］ DEAN W C．Frequency Analysis for Gravity and Magnetic Interpretation［J］．Geophysics，1958，23：97－127．

［22］ HOU Zhongchu，LIU Xiufang．Gravity and Magnetic Field Conversion Using Two Dimension Fourier Transform［J］．Journal of Beijing Normal University，1978，14（2）：54－69．（侯重初，劉秀芳．用二維富氏變換進行重磁位場的轉(zhuǎn)換［J］．北京師范大學(xué)學(xué)報，1978，14（2）：54－69．）

［23］ ZHANG Chuanyin，GUO Chunxi，CHEN Junyong，et al．EGM2008 and Its Application Analysis in Chinese Mainland［J］．Acta Geodaetica et Cartographica Sinica，2009，38（4）：283－289．（章傳銀，郭春喜，陳俊勇，等．EGM2008地球重力場模型在中國大陸適用性分析［J］．測繪學(xué)報，2009，38（4）：283－289．）

［24］ XU Xiaosu，ZHANG Yiqun．Application of Modified Terrain Entropy Algorithm in Terrain Aided Navigation［J］．Journal of Chinese Inertial technology，2008，16（5）：595－598．（徐曉蘇，張逸群．改進的地形熵算法在地形輔助導(dǎo)航中的應(yīng)用［J］．中國慣性技術(shù)學(xué)報，2008，16（5）：595－598．）

（責(zé)任編輯：陳品馨）

Wave Number Domain Iterative Tikhonov Regularization Method for Downward Continuation of Airborne Gravity Data

SUN Wen1，WU Xiaoping1，WANG Qingbin1，LIU Xiaogang2，3，4，WANG Kai5
1．Institute of Geospatial Information，Information Engineering University，Zhengzhou 450001，China；2．Xi’an Research Institute of Surveying and Mapping，Xi’an 710054，China；3．National Key Laboratory of Geo－Information Engineering，Xi’an 710054，China；4．Key laboratory of Geo－space Environment and Geodesy of Ministry of Education，Wuhan University，Wuhan 430079，China；5．Center for Engineering Design and Research under the Headquarters of General Equipment，Beijing 100028，China

The application of Tikhonov regularization method and space domain iterative Tikhonov method in the downward continuation was studied in this paper．There were some problems in these two methods．For the former，both the GCV and the L－curve methods showed the unconvergence property during the regularization parameter determination process，while for the latter，the RMS curve was unconvergent．The wave number domain iterative method was introduced，which showed better behavior in the convergence property than the above two methods．For the determination of iterated times，a method based on gravity anomaly variance entropy was proposed．An experiment was carried out using four methods and result showed that，compared to the other three methods，the wave number domain iterative method was better in precision and efficiency，which was worth to be wildly used in the downward continuation of airborne gravity data．

airborne gravity data；downward continuation；iterative method；wave number domain；regularization parameter；Tikhonov regularization；gravity anomaly variance entropy

SUN Wen（1987—），male，PhD candidate，majors in physical geodesy．

WANG Kai

P223

A

1001－1595（2014）06－0566－09

國家自然科學(xué)基金（41104047；41274029；41304022）；國家高技術(shù)研究發(fā)展專項（2013AA122502）；地球空間環(huán)境和大地測量教育部重點實驗室開放基金（11－01－03）

2013－01－11

孫文（1987—），男，博士生，研究方向為物理大地測量。

王凱

SUN Wen，WU Xiaoping，WANG Qingbin，et al．Wave Number Domain Iterative Tikhonov Regularization Method for Downward Continuation of Airborne Gravity Data［J］．Acta Geodaetica et Cartographica Sinica，2014，43（6）：566－574．（孫文，吳曉平，王慶賓，等．航空重力數(shù)據(jù)向下延拓的波數(shù)域迭代Tikhonov正則化方法［J］．測繪學(xué)報，2014，43（6）：566－574．）

10．13485／j．cnki．11－2089．2014．0080

修回日期：2013－10－03

E－mail：gravitysunwen＠gmail．com

E－mail：wsdbld＠aliyun．com