馮淑萍，高延民

（1.西安測繪總站，陜西 西安 710054）

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

馮淑萍1，高延民1

（1.西安測繪總站，陜西 西安 710054）

根據(jù)觀測面和延拓面測量數(shù)據(jù)的Poisson積分平面近似關(guān)系，結(jié)合快速傅立葉變換算法，將向下延拓轉(zhuǎn)換到頻率域進行計算，并采用迭代Tikhonov正則化方法，克服計算的不穩(wěn)定性，提高計算結(jié)果的精度，實現(xiàn)了航空重力測量數(shù)據(jù)的向下延拓。最后采用模擬航空重力測量數(shù)據(jù)驗證了該算法的有效性，取得了較好的延拓結(jié)果。

航空重力；向下延拓；正則化參數(shù)；迭代Tikhonov正則化法；快速傅立葉變換

在實際的工作中，航空重力測量常常是在起伏的航線上進行的，然而重力資料的定量解釋方法要求測量數(shù)據(jù)分布在一個平面上，因此，需要將實測資料向下延拓到一個平面上。而向下延拓是一個典型的不適定問題[1-2]，主要表現(xiàn)為計算的不穩(wěn)定性。隨著向下延拓深度的增大，會對重力測量數(shù)據(jù)中的高頻干擾信號起顯著的放大作用，從而不能分辨有效信號。為了解決這個問題，本文采用迭代Tikhonov正則化方法來實現(xiàn)航空重力測量數(shù)據(jù)的向下延拓。

1 向下延拓的基礎(chǔ)模型

航空重力測量數(shù)據(jù)向下延拓的基本原理如圖1所示。

圖1 向下延拓示意圖

圖中，Δgh（ξ，η）表示觀測面上的航空重力測量數(shù)據(jù)；Δg0（x，y）表示延拓面上的航空重力測量數(shù)據(jù)。計算觀測面以下至場源以上z=0這個平面的航空重力測量數(shù)據(jù)稱為航空重力測量數(shù)據(jù)的向下延拓。

根據(jù)航空重力測量數(shù)據(jù)向上延拓公式，可得觀測面航空重力測量數(shù)據(jù)與延拓面重力數(shù)據(jù)之間的近似關(guān)系公式[3-6]：

式中，h是向下延拓深度；r是延拓面上點（x，y，0）與觀測面上點（ξ，η，h）之間的距離。式（1）是第一類Fredholm積分方程，具有實對稱核，且可以表示為二維卷積：

其中：

將式（2）轉(zhuǎn)換到譜域里計算，其譜表達式為：

因為：

其中，f=（u2+v2）1/2；u、v分別表示空域變量x、y對應(yīng)的頻域變量。將式（5）代入式（4），可以得到：

通過變換可以得到：

對式（7）兩端進行Fourier逆變換，得到延拓面重力數(shù)據(jù)：

從式（8）可以看出，向下延拓的基礎(chǔ)數(shù)學(xué)模型原理比較簡單，但是由于其向下延拓算子 的不穩(wěn)定性對航空重力測量數(shù)據(jù)中的高頻噪聲有著顯著的放大作用，在延拓深度較大時會導(dǎo)致延拓結(jié)果精度不高，要解決這種不適定問題，可引入正則化因子。因此，本文研究了迭代Tikhonov正則化法數(shù)學(xué)模型。

2 迭代Tikhonov正則化方法

式（2）可以修改為：

式中，K表示第一類Fredholm積分算子。對于求解不適定問題采用Tikhonov正則化是一個常用的方法，它指的是求解一個極小化的正則化泛函，即

式中，α為正則化參數(shù)，用于平衡不穩(wěn)定性和光滑性。上式等價于如下的Euler方程：

K*為算子K的伴隨算子，由于航空重力測量數(shù)據(jù)向下延拓的Fredholm算子K為對稱的線性緊算子，有：

因此，式（11）可以表示為：

即

對式（14）兩邊同時作傅立葉變換，得到：

經(jīng)過調(diào)整，可以得到：

考慮到式（5）、（6），式（16）可以變形為：

在式（17）兩端進行Fourier逆變換，即可得到航空重力測量數(shù)據(jù)向下延拓的Tikhonov正則化法公式為：

由于Tikhonov正則化的飽和效應(yīng)，使得正則化解與準確解的誤差估計不能達到階數(shù)最優(yōu)，迭代的Tikhonov正則化對此進行了改進，其迭代形式如下：

將上式變形為：

根據(jù)數(shù)學(xué)歸納法，式（22）又可以寫為：

根據(jù)式（20），式（23）可以改化為：

再根據(jù)式（5）、（6），經(jīng)過化簡，可以得到：

在式（24）兩端進行Fourier逆變換，可得到航空重力測量數(shù)據(jù)向下延拓的迭代Tikhonov正則化法公式為[7～13]：

3 數(shù)值實驗與結(jié)果分析

3.1 航空重力測量數(shù)據(jù)仿真

采用2 160階的EGM2008地球重力場模型，對中國某地區(qū)的航空重力測量數(shù)據(jù)和地面重力數(shù)據(jù)進行仿真。將航空重力測量飛機的飛行高度設(shè)為3 000 m，分辨率設(shè)為5'×5'。圖2和圖3分別表示觀測面理論航空重力測量數(shù)據(jù)和延拓面理論重力數(shù)據(jù)。為了驗證本課題向下延拓算法的有效性，在觀測面理論航空重力測量數(shù)據(jù)中加入均值為0、標準差為3 mGal的高斯白噪聲，其等值線圖如圖4所示。

3.2 迭代Tikhonov正則化法實驗結(jié)果

在對迭代Tikhonov正則化方法的延拓精度進行測試時，首先需要對延拓誤差與正則化因子的關(guān)系進行驗證。表1表示迭代Tikhonov正則化法在取得最小延拓誤差時對應(yīng)的正則化因子和迭代次數(shù)。因此，表1的統(tǒng)計結(jié)果進一步證實了所得出的結(jié)論。

表1 最小延拓誤差對應(yīng)的正則化因子和迭代次數(shù)/mGal

圖2 觀測面理論航空重力測量數(shù)據(jù)等值線/mGal

圖3 延拓面理論重力數(shù)據(jù)等值線/mGal

圖4 含有高斯白噪聲的觀測面航空重力測量數(shù)據(jù)等值線/mGal

圖5給出了迭代Tikhonov正則化法在取得最小延拓誤差，即α=8.251時對應(yīng)的延拓結(jié)果的等值線圖。

從圖3和圖5的比較來看，迭代Tikhonov正則化方法延拓結(jié)果對延拓面上理論重力數(shù)據(jù)的逼近效果較好。

確定合適的正則化參數(shù)是迭代Tikhonov正則化法獲取最優(yōu)解的關(guān)鍵，如果α遠大于1，將導(dǎo)致逼近問題對原問題過于光滑，使計算結(jié)果與原問題的解相差甚遠；相反，如果α遠小于1，則逼近問題未能很好地改良算子的譜，正則化迭代法的效果不明顯。

圖5 迭代Tikhonov正則化法延拓結(jié)果等值線/mGal

[1] 欒文貴.場位解析延拓的穩(wěn)定化算法[J].地球物理學(xué)報,1983,26(3):263-274

[2] 梁錦文.位場向下延拓的正則化方法[J].地球物理學(xué)報,1989,32(5):600-608

[3] 王興濤,夏哲仁,石磐,等.航空重力測量數(shù)據(jù)向下延拓方法比較[J].地球物理學(xué)報,2004,47(6):1 017-1 022

[4] 羅志才,寧津生,晁定波.衛(wèi)星重力梯度向下延拓的譜方法[J].測繪學(xué)報,1997,26(2):168-175

[5] 孫中苗.航空重力測量理論方法及應(yīng)用研究[D].鄭州:信息工程大學(xué),2004

[6] 劉曉剛,李珊珊,吳星.衛(wèi)星重力梯度數(shù)據(jù)的向下延拓[J].大地測量與地球動力學(xué),2011,31(1):132-137

[7] TIKHONOV A N, ARSENIN V Y. Solution of Certain Integral Equations of the First Kind[J]. J Assoc Comput Mach, 1962 (9):84-97

[8] 成怡,郝燕玲,劉繁明.航空重力測量數(shù)據(jù)向下延拓及其影響因素分析[J].系統(tǒng)仿真學(xué)報,2008,20(8):2 190-2 194

[9] 郝燕玲,成怡,孫楓,等.Tikhonov正則化向下延拓算法仿真實驗研究[J].儀器儀表學(xué)報,2008,29(3):605-609

[10] 鄧凱亮,暴景陽,黃謨濤,等.航空重力數(shù)據(jù)向下延拓的Tikhonov正則化法仿真研究[J].武漢大學(xué)學(xué)報(信息科學(xué)版),2010,35(12):1 414-1 417

[11] 曾小牛,李夕海,韓紹卿,等.位場向下延拓三種迭代方法之比較[J].地球物理學(xué)進展,2011,26(3): 908-915

[12] 鄧凱亮,黃謨濤,暴景陽,等.向下延拓航空重力數(shù)據(jù)的Tikhonov雙參數(shù)正則化法[J].測繪學(xué)報,2011,40(6):690-696

[13] 蔣濤,李建成,王正濤,等.航空重力向下延拓病態(tài)問題的求解[J].測繪學(xué)報,2011,40(6):112-122

P223

B

1672-4623（2016）10-0053-03

10.3969/j.issn.1672-4623.2016.10.015

馮淑萍，工程師，主要從事航空重力研究。

2015-07-13。

項目來源：國家自然科學(xué)基金資助項目（41304022）；國家重點基礎(chǔ)研究發(fā)展計劃資助項目（61322201，2013CB733303）；高分專項青年創(chuàng)新基金資助項目（GFZX04060103-5-12）。