向下延拓航空重力數據的Tikhonov雙參數正則化法

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

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

向下延拓航空重力數據的Tikhonov雙參數正則化法

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

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

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

向下延拓；航空重力數據；Tikhonov雙參數正則化法；GCV法

1 引 言

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

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

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

2 向下延拓原理

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

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

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

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

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

取

則式（3）可寫為

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

積分區域內大地水準面上的重力異常的最小二乘解為

3 Tikhonov雙參數正則化法

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

3.1 向下延拓模型的譜分解

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

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

3.2 Tikhonov雙參數正則化法

Tikhonov正則化法的實質是尋找適當的濾波因子，用相鄰的適定解去逼近原問題的解［14］。對于觀測方程式（4）而言，其正則化函數為

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

Tikhonov正則化解的濾波因子為

將式（10）進行譜分解，得到

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

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

Tikhonov雙參數正則化法的濾波因子為

式中，k為截斷參數。

Tikhonov雙參數正則化法的譜分解形式為

Tikhonov雙參數正則化解的均方誤差為

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

3.3 正則化參數的選取

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

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

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

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

式中

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

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

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

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

圖2 正則化參數分布圖Fig.2 Characteristic regularization parameters

4 仿真試驗及分析

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

4.1 數據準備

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

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

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

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

表1 仿真區域各項重力異常特性的統計Tab.1 Statistics of gravity anomalies at simulation area mGal

4.2 試驗步驟

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

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

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

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

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

4.3 結果比較與分析

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

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

基于三種計算方法獲得的大地水準面重力異常與原始海面重力異常比較的差值統計結果見表2。

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

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

由表2可以看出：

（1）方法1的法矩陣的條件數達到了14 886.62，遠遠大于1 000的病態界點，表明此時的法矩陣是嚴重病態的。

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

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

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

5 結束語

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

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

Tikhonov雙參數正則化法的關鍵是獲取最優的截斷參數和正則化參數。采用GCV函數作為判斷標準，GCV函數隨著截斷參數和正則化參數的同時引進而作了相應的改進。通過基于EGM2008重力場模型仿真的航空重力數據進行向下延拓試驗，驗證了Tikhonov雙參數正則化法要優于僅有正則化參數的Tikhonov正則化法。

［1］ SUN Zhongmiao.Theory，Methods and Applications of Airborne Gravimetry［D］.Zhengzhou：Information Engineering University，2004.（孫中苗.航空重力測量理論、方法及應用研究［D］.鄭州：信息工程大學，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.（楊強文，吳曉平.航空重力數據向下延拓的倒錐法［J］.解放軍測繪學院學報，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］.測繪學報，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.（王興濤，夏哲仁，石磐，等.航空重力測量數據向下延拓方法比較［J］.地球物理學報，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.（王興濤，石磐，朱非洲.航空重力測量數據向下延拓的正則化算法及其譜分解［J］.測繪學報，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］.儀器儀表學報，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.（成怡，郝燕玲，劉繁明.航空重力測量數據向下延拓及其影響因素分析［J］.系統仿真學報，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］.武漢：中國科學院測量與地球物理研究所，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.（鄧凱亮.海域多源重力數據的處理、融合及應用研究［D］.大連：海軍大連艦艇學院，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］.大地測量學與地球動力學，2002，33（3）：11－14.）

［14］ LIU Jijun.Regularization Method and Application of the Ill－posed Equation［M］.Beijing：Science Press，2005.（劉繼軍.不適定問題的正則化方法及應用［M］.北京：科學出版社，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和正則化算法的低軌衛星軌道平滑［J］.地球物理學報，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］.測繪學報，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）；國家自然科學基金（41074002；61071006）；國家重大海洋勘測專項資助（4200301）

叢樹平）

2010－10－11

2011－03－17

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

E－mail：dengkailiang036＠163.com