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

鄧凱亮，黃謨濤，暴景陽，歐陽永忠，陸秀平，劉傳勇

1.大連艦艇學(xué)院 海測工程系，遼寧大連116018；2.海軍海洋測繪研究所，天津300061

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

鄧凱亮1，2，黃謨濤2，暴景陽1，歐陽永忠2，陸秀平2，劉傳勇2

1.大連艦艇學(xué)院 海測工程系，遼寧大連116018；2.海軍海洋測繪研究所，天津300061

為避免正則化參數(shù)對向下延拓過程可靠成分的修正影響，提出Tikhonov雙參數(shù)正則化法。引進(jìn)截斷參數(shù)，將法矩陣的奇異值分為相對較大的奇異值（可靠部分）和相對較小的奇異值（不可靠部分）；引進(jìn)正則化參數(shù)，只對法矩陣的小奇異值進(jìn)行修正，以抑制高頻誤差對向下延拓解的影響。采用廣義交互確認(rèn)法（GCV）確定截斷參數(shù)和正則化參數(shù)?；贓GM2008重力場模型仿真一組航空重力數(shù)據(jù)，驗證該方法對航空重力數(shù)據(jù)向下延拓過程的有效性。

向下延拓；航空重力數(shù)據(jù)；Tikhonov雙參數(shù)正則化法；GCV法

1 引 言

航空重力測量不僅能快速經(jīng)濟(jì)地獲得中高頻重力場信息，而且能遠(yuǎn)距離測量，有效地突破測量范圍的局限性，是衛(wèi)星重力測量和海面重力測量的重要補充［1－2］。

在大地水準(zhǔn)面等應(yīng)用上，需要將航空重力數(shù)據(jù)向下延拓到大地水準(zhǔn)面上，并采用球近似。向下延拓的基本方法是逆Poisson方法，這是一個信號放大的非平穩(wěn)過程，很小的觀測噪聲往往引起待估參數(shù)較大的誤差，屬于不適定問題［2］。對此，國內(nèi)外學(xué)者作了許多有益的研究［3－10］，提出迭代法［3］、配置法［4］、倒錐法［5］、直接代表法［6］、Tikhonov正則化法［7－10］等；文獻(xiàn)［7—8］對上述方法在航空重力測量數(shù)據(jù)向下延拓中的應(yīng)用進(jìn)行了比較分析，表明造成向下延拓結(jié)果不適定的主要原因是觀測誤差的高頻部分被法矩陣的小奇異值嚴(yán)重放大。Tikhonov正則化方法通過正則化參數(shù)對法矩陣小奇異值的修正，能有效抑制高頻誤差對參數(shù)估值的影響，同時驗證了該方法的優(yōu)越性。但是正則化參數(shù)不僅修正了較小的奇異值，而且對相對較大的奇異值也進(jìn)行了修正，這樣使模型中的可靠成分也受到了畸變的影響［11－12］。

為了避免正則化參數(shù)對向下延拓過程可靠成分的修正影響，筆者提出Tikhonov雙參數(shù)正則化法?；舅悸肥牵涸谙蛳卵油刈V分解的基礎(chǔ)上，引進(jìn)截斷參數(shù)，將法矩陣的奇異值分為相對較大的奇異值（可靠部分）和相對較小的奇異值（不可靠部分）；引進(jìn)正則化參數(shù)，只對法矩陣的小奇異值進(jìn)行修正，以抑制高頻誤差對向下延拓解的影響。參數(shù)的選擇采用廣義交互確認(rèn)法（generalized cross－validation，GCV）?；贓GM2008重力場模型仿真了一組航空重力數(shù)據(jù)，以驗證該方法對航空重力數(shù)據(jù)向下延拓過程的有效性。

2 向下延拓原理

根據(jù)位理論和Dirichlet邊值理論，重力異常的向上延拓可由Poisson積分公式［3］計算，公式為

式中，R為地球半徑；h是延拓的高度；r為延拓高度處的球半徑（r＝R＋h）；ω為以R為半徑的球面；K（r，ψ，R）是積分核函數(shù)，其公式［3］為

式中，ψ是計算點與球面流動點之間的角距，其表達(dá)式為

一般情況下，航空重力異常和大地水準(zhǔn)面上的重力異常均是離散值，故將式（1）離散化，表達(dá)式可化為

式中，ωc是有效積分范圍；N為積分范圍ωc內(nèi)的測量點數(shù)；Δgω－ωc（r，φ，λ）為遠(yuǎn)區(qū)影響，引起的誤差為截斷誤差，可由重力場模型近似計算。

取

則式（3）可寫為

當(dāng)航空重力異常Δgh（r）向下延拓到大地水準(zhǔn)面時，Δgh（r）是已知觀測量，大地水準(zhǔn)面上的重力異常Δg（R，θj，λj）是待求量，這是Poisson積分的逆問題。

積分區(qū)域內(nèi)大地水準(zhǔn)面上的重力異常的最小二乘解為

3 Tikhonov雙參數(shù)正則化法

航空重力測量數(shù)據(jù)的向下延拓是信號放大的不適定過程，很小的觀測噪聲往往能引起待估參數(shù)較大的誤差，這正是正則化法所要解決的問題。在眾多的正則化方法中，以Tikhonov正則化法應(yīng)用最為廣泛。

3.1 向下延拓模型的譜分解

在航空測量數(shù)據(jù)的向下延拓中，觀測方程系數(shù)矩陣的奇異值單調(diào)遞減，最小二乘解被放大的觀測噪聲污染。對系數(shù)矩陣A作奇異值分解（SVD）

式中，U、V分別是A的左右奇異向量，U＝［u1u2… un］，V＝［v1v2… vn］，且滿足是對角矩陣，其對角線上的元素為A的遞減的奇異值λi。將式（6）代入式（5）中得最小二乘解的譜分解形式

3.2 Tikhonov雙參數(shù)正則化法

Tikhonov正則化法的實質(zhì)是尋找適當(dāng)?shù)臑V波因子，用相鄰的適定解去逼近原問題的解［14］。對于觀測方程式（4）而言，其正則化函數(shù)為

式中，α＞0是正則化參數(shù)；‖x‖表示x的范數(shù)。根據(jù)式（9）的約束條件，可得Tikhonov正則化解

Tikhonov正則化解的濾波因子為

將式（10）進(jìn)行譜分解，得到

由式（12）、式（13）可以看出，正則化參數(shù)不僅對小的奇異值進(jìn)行了修正，也對相對較大的奇異值即模型的可靠部分進(jìn)行了修正，使得最小二乘解發(fā)生了畸變。

對此，設(shè)計了Tikhonov雙參數(shù)正則化法。主要思路是：引進(jìn)Tikhonov截斷參數(shù)，將奇異值分為相對較大的奇異值和相對較小的奇異值，即可靠部分和不可靠部分；沿用正則化參數(shù)，使其只對模型的不可靠部分進(jìn)行修正。這樣既抑制了小奇異值對觀測噪聲高頻段的放大影響，又避免模型的可靠部分受到修正影響。

Tikhonov雙參數(shù)正則化法的濾波因子為

式中，k為截斷參數(shù)。

Tikhonov雙參數(shù)正則化法的譜分解形式為

Tikhonov雙參數(shù)正則化解的均方誤差為

由式（16）可以看出，Tikhonov雙參數(shù)正則化法的估值誤差包括兩部分：前者是測量誤差引起的估值誤差，隨著截斷參數(shù)k的增大而減小，隨著正則化參數(shù)α的增大而減??；后者是正則化引起的估值誤差，隨著截斷參數(shù)k的增大而減小，隨著正則化參數(shù)α的增大而增大。因此，存在最佳的截斷參數(shù)k和正則化參數(shù)α，使得由測量誤差與正則化誤差引起的估值誤差達(dá)到最佳的平衡。

3.3 正則化參數(shù)的選取

近年來，統(tǒng)計學(xué)界和大地測量學(xué)界提出了多種方法選擇正則化參數(shù)。廣義交互確認(rèn)法（generalized cross－validation，GCV）得到了廣泛的應(yīng)用［15－16］，該法是基于“遺漏一個”（leave－out－one）的思想。該思想轉(zhuǎn)化為一個便于操作的最小化問題時，就是選擇合適正則化參數(shù)GCV函數(shù)最小。該函數(shù)定義為

Qα是所謂的影響矩陣，由Axα＝QαL 定義。Wahba給出了求得GCV函數(shù)的另一種形式

式中，n為觀測值個數(shù)；tr為矩陣的跡。最佳的截斷參數(shù)k和正則化參數(shù)α對應(yīng)于GCV函數(shù)的最小值。

對于Tikhonov雙參數(shù)正則化法，基于Qα的定義，容易得到

式中

式中，VN、UN由Poisson方程的法矩陣的譜分解式N＝UNdiag（λN）VTN得到；λN（k，α）由式（20）得到

在數(shù)據(jù)量大的情況下執(zhí)行GCV，應(yīng)考慮正則化參數(shù)的收斂速度［17］：在得到法矩陣的奇異值λN后，假設(shè)正則化參數(shù)為q個，則取αq＝λn，ratio＝，修正參數(shù)為

在本文的仿真試驗中，法矩陣的奇異值λN分布圖見圖1，所考慮的正則化參數(shù)α的分布圖見圖2。

圖1 法矩陣的奇異值分布圖Fig.1 Singular values characteristic of the Poisson normal matrix

圖2 正則化參數(shù)分布圖Fig.2 Characteristic regularization parameters

4 仿真試驗及分析

為了驗證基于Tikhonov雙參數(shù)正則化法對航空重力數(shù)據(jù)進(jìn)行向下延拓結(jié)果的有效性，設(shè)計了航行高度為5km的仿真試驗。

4.1 數(shù)據(jù)準(zhǔn)備

基于EGM2008重力場模型計算的重力異常作為仿真試驗的數(shù)據(jù)基礎(chǔ)。EGM2008重力場模型［18］是由NGA（national geospatial－intelligence agency）釋放的全球超高階地球重力場模型，由衛(wèi)星重力測量、衛(wèi)星測高和大地水準(zhǔn)面重力觀測等資料聯(lián)合解算得到，模型階數(shù)達(dá)到2 160。EGM2008重力場模型導(dǎo)出的重力異常在我國大陸的總體精度為10.5mGal［19］（1Gal＝0.01m／s2）。

基于EGM2008重力場模型計算某區(qū)域2 160階的大地水準(zhǔn)面上的重力異常Δg0（圖3）和5km高空面的重力異常Δg5（圖4）。它的范圍為1°× 1°，格網(wǎng)間距為2′×2′。仿效移去－恢復(fù)技術(shù)的應(yīng)用，取EGM2008重力場模型的360階作為參考模型，得到參考重力異常Δg′0和參考航空重力異常Δg′5。各項重力異常的統(tǒng)計特性見表1。

圖3 仿真的大地水準(zhǔn)面重力異常Δg0Fig.3 Simulative ground gravity dataΔg0

圖4 仿真的航空重力異常Δg5Fig.4 Simulative airborne gravity dataΔg5

表1 仿真區(qū)域各項重力異常特性的統(tǒng)計Tab.1 Statistics of gravity anomalies at simulation area mGal

4.2 試驗步驟

為了模擬航空重力測量數(shù)據(jù)的處理結(jié)果的誤差，對仿真的航空重力異常Δg5引入三種觀測誤差分別是零均值的白噪聲：e1（σ＝±1mGal），e2（σ＝±3mGal），e3（σ＝±5mGal）。試驗步驟如下：

（1）模擬產(chǎn)生3個量級的誤差并加入Δg5中，得到含誤差的仿真航空重力異常Δg5e1、Δg5e2、Δg5e3，并從中移去參考航空重力異常Δg′5，得到殘差航空重力異常δΔg5e1、δΔg5e2、δΔg5e3。

（2）基于Tikhonov雙參數(shù)正則化法，利用GCV法，求得各殘差航空重力異常δΔg5e1、δΔg5e2、δΔg5e3向下延拓時的截斷參數(shù)k和正則化參數(shù)α（圖5）。為了驗證雙參數(shù)法的效果，設(shè)計了Tikhonov法的兩種方案，方案1是只有正則化參數(shù)的Tikhonov法，方案2是Tikhonov雙參數(shù)正則化法。

（3）按式（12）對各殘差航空重力異常δΔg5e1、 δΔg5e2、δΔg5e3向下延拓，得到延拓后的殘差大地水準(zhǔn)面重力異常δΔg0e1、δΔg0e2、δΔg0e3，并恢復(fù)參考重力異常Δg′0，得到大地水準(zhǔn)面重力異常Δg0e1、Δg0e2、Δg0e3。

圖5 各誤差條件下兩種方案確定的正則化參數(shù)示意圖Fig.5 Regularization parameters given by two projects each under the conditions of three error levels

4.3 結(jié)果比較與分析

為了分析Tikhonov雙參數(shù)正則化法對航空重力數(shù)據(jù)基于逆Poisson積分方程的向下延拓模型的有效性，將向下延拓得到的大Δg0e1、Δg0e2、Δg0e3與仿真的大地水準(zhǔn)面重力異常Δg0進(jìn)行比較，以檢驗向下延拓后數(shù)據(jù)的精度。

為了驗證本方法的效果，設(shè)計了三種計算方法：方法1，直接逆Poisson積分方程的向下延拓法；方法2，基于GCV法的只有正則化參數(shù)的Tikhonov法的向下延拓法；方法3，基于GCV法的Tikhonov雙參數(shù)正則化法的向下延拓法。

基于三種計算方法獲得的大地水準(zhǔn)面重力異常與原始海面重力異常比較的差值統(tǒng)計結(jié)果見表2。

表2 與模擬的重力異常真值Δg0比較的差值統(tǒng)計Tab.2 Statistics of the differences between calculated values（Δg0e1，Δg0e2，Δg0e3）and the simulative true valuesΔg0in the comparison

基于三種計算方法獲得的大地水準(zhǔn)面上的重力異常見圖6。

由表2可以看出：

（1）方法1的法矩陣的條件數(shù)達(dá)到了14 886.62，遠(yuǎn)遠(yuǎn)大于1 000的病態(tài)界點，表明此時的法矩陣是嚴(yán)重病態(tài)的。

（2）在設(shè)計的三種誤差條件下，方法2使得法矩陣的條件數(shù)控制在70以內(nèi)，表明此時的向下延拓結(jié)果具有良好的穩(wěn)定性；方法3將法矩陣的條件數(shù)控制在30以內(nèi)，表明該方法能有效抑制微小誤差被小奇異值放大程度，在保證向下延拓結(jié)果的穩(wěn)定性方面具有優(yōu)越性。

（3）在設(shè)計的三種誤差條件下，采用方法2，標(biāo)準(zhǔn)差分別由原來的20.28mGal、54.53mGal、83.94mGal降低到4.36mGal、6.45mGal、 8.18mGal，而方法3，能進(jìn)一步提高成果的精度，分別降到了3.27mGal、4.29mGal和6.08mGal。表明方法2和方法3都可有效地抑制觀測誤差對向下延拓結(jié)果的影響，提高成果的精度，后者的改善效果更加明顯。

圖6 各誤差條件下三種方法的向下延拓成果Fig.6 Downward continuation results given by three methods each under the conditions of three error levels

5 結(jié)束語

利用直接逆Poisson積分方程進(jìn)行航空重力數(shù)據(jù)向下延拓，延拓結(jié)果極不穩(wěn)定（用條件數(shù)表示）。在Tikhonov正則化處理過程中，基于GCV法得到與誤差環(huán)境相匹配的正則化參數(shù)，能抑制小奇異值對觀測誤差的放大影響，改進(jìn)延拓結(jié)果的穩(wěn)定性和精度。但是此時正則化參數(shù)對法矩陣的可靠部分即相對大的奇異值也進(jìn)行了改正，這樣將導(dǎo)致模型的可靠部分發(fā)生畸變。

在對向下延拓過程法矩陣進(jìn)行譜分解的基礎(chǔ)上，為避免正則化參數(shù)對法矩陣奇異值可靠成分的修正影響，提出Tikhonov雙參數(shù)正則化法。即在正則化參數(shù)的基礎(chǔ)上，引進(jìn)截斷參數(shù)，將法矩陣的奇異值分為相對較大的奇異值（可靠部分）和相對較小的奇異值（不可靠部分），而正則化參數(shù)只對法矩陣的不可靠部分進(jìn)行修正，以抑制高頻誤差對向下延拓解的影響。

Tikhonov雙參數(shù)正則化法的關(guān)鍵是獲取最優(yōu)的截斷參數(shù)和正則化參數(shù)。采用GCV函數(shù)作為判斷標(biāo)準(zhǔn)，GCV函數(shù)隨著截斷參數(shù)和正則化參數(shù)的同時引進(jìn)而作了相應(yīng)的改進(jìn)。通過基于EGM2008重力場模型仿真的航空重力數(shù)據(jù)進(jìn)行向下延拓試驗，驗證了Tikhonov雙參數(shù)正則化法要優(yōu)于僅有正則化參數(shù)的Tikhonov正則化法。

［1］ SUN Zhongmiao.Theory，Methods and Applications of Airborne Gravimetry［D］.Zhengzhou：Information Engineering University，2004.（孫中苗.航空重力測量理論、方法及應(yīng)用研究［D］.鄭州：信息工程大學(xué)，2004.）

［2］ KERN M.An Analysis of the Combination and Downward Continuation of Satellite，Airborne and Terrestrial Gravity Data［D］.Calgary：University of Calgary，2003.

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

［4］ 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.

［5］ 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.）

［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］ 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.）

［9］ HAO Yanling，CHENG Yi，SUN Feng，et al.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.）

［10］ CHENG Yi，HAO Yanling，LIU Fanming.Analysis of Effect and Downward Continuation of Airborne Gravity Data［J］.Journal of System Simulation，2008，20（3）：605－609.（成怡，郝燕玲，劉繁明.航空重力測量數(shù)據(jù)向下延拓及其影響因素分析［J］.系統(tǒng)仿真學(xué)報，2008，20（3）：605－609.）

［11］ WANG Zhenjie.Research on the Regularization Solutions of Ill－Posed Problems in Geodesy［D］.Wuhan：Institute of Geodesy and Geophysics，Chinese Academy of Sciences，2003.（王振杰.大地測量中不適定問題的正則化解法研究［D］.武漢：中國科學(xué)院測量與地球物理研究所，2003.）

［12］ DENG Kailiang.Research on the Procession，Combination and Application of the Muti－source Gravity Data on the Sea［D］.DaLian：Dalian Naval Academy，2011.（鄧凱亮.海域多源重力數(shù)據(jù)的處理、融合及應(yīng)用研究［D］.大連：海軍大連艦艇學(xué)院，2011.）

［13］ SHEN Yunzhong，XU Houze.Spectral Decomposition Formula of Regularization Solution for Ill－posed Equation［J］.Journal of Geodesy and Geodynamics，2002，33（3）：11－14.（沈云中，許厚澤.不適定方程正則化算法的譜分解式［J］.大地測量學(xué)與地球動力學(xué)，2002，33（3）：11－14.）

［14］ LIU Jijun.Regularization Method and Application of the Ill－posed Equation［M］.Beijing：Science Press，2005.（劉繼軍.不適定問題的正則化方法及應(yīng)用［M］.北京：科學(xué)出版社，2005.）

［15］ KUSCHE J，KLEES R.Regularization of Gravity Field Estimation from Satellite Gravity Gradients［J］.Journal of Geodesy，2002，76（6－7）：359－368.

［16］ LIU Xianglin，DITMAR P.Smoothing a Satellite Orbit on the Basis of B－spline and Regularization［J］.Chinese Journal of Geophysics，2006，49（1）：99－105.（柳響林，DITMAR P.基于B－spline和正則化算法的低軌衛(wèi)星軌道平滑［J］.地球物理學(xué)報，2006，49（1）：99－105.）

［17］ GOLUB G H，VON MATT U.Generalized Cross－validation for Large－scale Problems［J］.Journal of Computational and Graphical Statistics，1997，6（1）：1－34.

［18］ PAVLIS N K，HOLMES S A，KENYON S C，et al.An Earth Gravitational Model to Degree 2160：EGM2008［J］.Geophysical Research Abstracts，2008，10：13－18.

［19］ ZHANG Chuanyin，GUO Chunxi，CHEN Junyong，et al.EGM2008and 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.）

Tikhonov Two－parameter Regulation Algorithm in Downward Continuation of Airborne Gravity Data

DENG Kailiang1，2，HUANG Motao2，BAO Jingyang1，OUYANG Yongzhong2，LU Xiuping2，LIU Chuanyong2
1.Department of Hydrography and Cartography，Dalian Naval Academy，Dalian116018，China；2.Naval Institute of Hydrographic Surveying and Charting，Tianjin 300061，China

In order to avoid the correct effect provided by regulation parameter on the downward continuation credible portion，Tikhonov two－parameter regulation algorithm is presented.Using the truncation parameter，the singular values of the normal matrix are divided into the relative large singular values（the credible portion）and the relative small singular values（the incredible portion）.Using regulation parameter，only the relative small singular values are corrected，and the effect induced by the high frequency error is suppressed.The truncation parameter and the regulation parameter are determined by GCV method.Using the simulative airborne gravity data based on the EGM2008，as true values of the gravity field，the results of numerical calculation examples have clearly demonstrated that this method for downward continuation of airborne gravity data is fairly efficient，and shown significant advantages.

downward continuation；airborne gravity data；Tikhonov two－parameter regulation algorithm；GCV method

DENG Kailiang（1983—），male，PhD candidate，majors in marine geodetic survey.

1001－1595（2011）06－0690－07

P223

A

國家863計劃（2006AA06A202；2009AA121405）；國家自然科學(xué)基金（41074002；61071006）；國家重大海洋勘測專項資助（4200301）

叢樹平）

2010－10－11

2011－03－17

鄧凱亮（1983—），男，博士生，研究方向為海洋大地測量。

E－mail：dengkailiang036＠163.com