確定大地水準(zhǔn)面的Tikhonov最小二乘配置法

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

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

大地水準(zhǔn)面是定義正高高程系統(tǒng)的高程基準(zhǔn)面，也是反映地球內(nèi)部結(jié)構(gòu)和密度分布特征的物理面。確定高精度高分辨率的大地水準(zhǔn)面，已成為21世紀(jì)大地測(cè)量學(xué)科發(fā)展全局性的戰(zhàn)略目標(biāo)［1－3］。國內(nèi)外學(xué)者就大地水準(zhǔn)面的確定做了許多有益的研究，提出了Stokes理論、Molodensky理論、Bjerhammar理論和最小二乘配置理論等［4－6］。其中以統(tǒng)計(jì)理論為基礎(chǔ)的最小二乘配置理論，由于它具有能對(duì)多種類型的重力觀測(cè)量進(jìn)行聯(lián)合處理的特性，在大地水準(zhǔn)面確定的應(yīng)用上受到廣泛關(guān)注［7－16］。利用最小二乘配置法確定大地水準(zhǔn)面的關(guān)鍵是協(xié)方差函數(shù)的確定［7－9］，但是即使擬合的協(xié)方差函數(shù)能充分表達(dá)研究區(qū)域范圍內(nèi)重力場(chǎng)特性，由于在大地水準(zhǔn)面的確定過程中需要對(duì)協(xié)方差矩陣進(jìn)行求逆，而協(xié)方差矩陣的求逆過程是信號(hào)放大的非平穩(wěn)過程，協(xié)方差矩陣的小奇異值將放大觀測(cè)誤差對(duì)配置結(jié)果的影響，導(dǎo)致配置結(jié)果不穩(wěn)定且精度偏低［17］。當(dāng)已知的重力觀測(cè)量存在觀測(cè)誤差時(shí)，最小二乘配置法難以得到穩(wěn)定精確的大地水準(zhǔn)面解。

本文在確定大地水準(zhǔn)面的最小二乘配置法中，引入Tikhinov正則化法，對(duì)協(xié)方差矩陣進(jìn)行正則化處理，以抑制協(xié)方差矩陣的小奇異值對(duì)觀測(cè)誤差的放大影響，得到穩(wěn)定且高精度的大地水準(zhǔn)面高?；贓GM2008重力場(chǎng)模型計(jì)算了一組重力異常數(shù)據(jù)，以該重力異常作為基礎(chǔ)數(shù)據(jù)，聯(lián)合最小二乘配置法和Tikhonov正則化法確定大地水準(zhǔn)面，以驗(yàn)證該方法的有效性。

1 最小二乘配置法

1.1 基本原理

重力場(chǎng)的所有重力場(chǎng)觀測(cè)量都可看成空間平穩(wěn)隨機(jī)場(chǎng)的隨機(jī)量，任意重力場(chǎng)觀測(cè)量l可表示為

式中，Li為由擾動(dòng)位T表示重力場(chǎng)觀測(cè)量l的線性泛函算子；e為觀測(cè)噪聲，其向量形式為

式中，t＝LT，是l的信號(hào)部分。

最小二乘配置公式［4－5］為

依據(jù)式（3），以重力異常Δg作為基礎(chǔ)數(shù)據(jù)，基于最小二乘配置法確定大地水準(zhǔn)面的公式為

式中，CNΔg為大地水準(zhǔn)面與重力異常的互協(xié)方差矩陣；CΔgΔg為重力異常的協(xié)方差矩陣；Cnn是觀測(cè)噪聲的協(xié)方差陣。

1.2 協(xié)方差函數(shù)的確定

擾動(dòng)位T的協(xié)方差函數(shù)可表示為［7］

依據(jù)位理論、邊值理論和協(xié)方差傳播定律，基于移去－恢復(fù)思想，重力異常與大地水準(zhǔn)面之間協(xié)方差函數(shù)可寫為

球面重力異常經(jīng)驗(yàn)協(xié)方差C（ψ）的公式［11，18］為

式中，θ為地心余緯；λ為地心經(jīng)度；α為方位角。若記球面上重力異常的格網(wǎng)值為f，面積為B，則式（7）的離散形式為

依據(jù)重力異常經(jīng)驗(yàn)協(xié)方差C（ψ）擬合式（6）中的第二項(xiàng)，求得常數(shù)a0、A和RB，進(jìn)而確定重力異常和大地水準(zhǔn)面之間的互協(xié)方差函數(shù)。

2 Tikhonov正則化法

依據(jù)式（4），在大地水準(zhǔn)面的確定過程中需要對(duì)協(xié)方差矩陣CΔgΔg進(jìn)行求逆，而協(xié)方差矩陣的條件數(shù)較大，其求逆過程是信號(hào)放大的非平穩(wěn)過程，小的觀測(cè)誤差往往會(huì)引起結(jié)果的較大誤差，屬于不適定問題［17］。正則化算法的實(shí)質(zhì)就是通過選擇合適的正則化參數(shù)來抑制觀測(cè)噪聲對(duì)參數(shù)估值的影響，以得到穩(wěn)定、精確的解。在眾多的正則化方法中，以Tikhonov正則化法應(yīng)用最為廣泛［19－21］。

2.1 Tikhonov正則化法

Tikhonov正則化法的實(shí)質(zhì)是用相鄰的適定解去逼近源問題的解［22］。取觀測(cè)值Δg的個(gè)數(shù)為q，令

則，x是一個(gè)q向量，式（4）可寫為

由式（10）可看出，是q個(gè)CNi的線性組合，由于CNi可由協(xié)方差函數(shù)計(jì)算得到，故只要得到穩(wěn)定精確的x，就能確定穩(wěn)定精確的。

將式（9）寫為

式中，A＝CΔgΔg＋Cnn，l＝Δg

則式（11）的正則化函數(shù)為

式中，α＞0是正則化參數(shù)；‖x‖表示x的范數(shù)。

根據(jù)式（12）的約束條件，可得Tikhonov正則化解

在正則化解的兩邊乘以CNΔg，則得到

所以求得穩(wěn)定精確的大地水準(zhǔn)面高的關(guān)鍵是得到穩(wěn)定精確的

式（11）在頻域內(nèi)的形式為

為了分析正則化參數(shù)的影響，引入均方誤差MSE

由式（16）可知，Tikhonov正則化法的估值誤差包括兩部分：前者是測(cè)量誤差引起的估值誤差，隨著正則化參數(shù)α的增大而減??；后者是正則化引起的估值誤差，隨著正則化參數(shù)α的增大而增大。如何選擇正則化參數(shù)α是整個(gè)正則化算法的關(guān)鍵。

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

近年來，統(tǒng)計(jì)學(xué)界和大地測(cè)量學(xué)界提出了多種方法選擇正則化參數(shù)，比如L曲線法［19－20］和廣義交互確認(rèn)法（generalized cross－validation，GCV）［21－26］等。這里選擇GCV法選擇正則化參數(shù)。

該函數(shù)定義為

式中，Qα是所謂的影響矩陣，由Axα＝QαL定義；n為觀測(cè)值個(gè)數(shù)；tr為矩陣的跡。

最佳的正則化參數(shù)α對(duì)應(yīng)于GCV函數(shù)的最小值。

本文的試驗(yàn)取q為100［19］，以選取的正則化參數(shù)的示意圖見圖1。

3 仿真試驗(yàn)及精度分析

為了驗(yàn)證聯(lián)合最小二乘配置法和Tikhonov正則化法確定大地水準(zhǔn)面的有效性，設(shè)計(jì)了以重力異常作為基礎(chǔ)數(shù)據(jù)確定大地水準(zhǔn)面的仿真試驗(yàn)。

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

基于EGM2008重力場(chǎng)模型計(jì)算的重力異常作為仿真試驗(yàn)的基礎(chǔ)數(shù)據(jù)。EGM2008重力場(chǎng)模型［27］是由 NGA（National Geospatial－intelligence Agency）釋放的全球超高階地球重力場(chǎng)模型，由衛(wèi)星重力測(cè)量、衛(wèi)星測(cè)高和地面重力觀測(cè)等資料聯(lián)合解算得到，模型階數(shù)達(dá)到2160。EGM2008重力場(chǎng)模型導(dǎo)出的重力異常在我國大陸的總體精度為10.5mGal（1mGal＝10－5m／s2）［28］。

基于EGM2008重力場(chǎng)模型分別計(jì)算山區(qū)、丘陵和海域2160的重力異常Δg山、Δg丘和Δg海（見圖1），大地水準(zhǔn)面高N山、N丘和N海（見圖2）。區(qū)域范圍都為1°×1°，格網(wǎng)間距都為2′×2′?？紤]移去－恢復(fù)技術(shù)的應(yīng)用，取EGM2008重力場(chǎng)模型的360階作為參考模型，得到參考重力異常Δg′山、Δg′丘和 Δg′海和參考大地水準(zhǔn)面高N′山、N′丘和N′海。重力異常和大地水準(zhǔn)面的統(tǒng)計(jì)特性見表1和表2。

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

表2 仿真區(qū)域大地水準(zhǔn)面特性的統(tǒng)計(jì)Tab.2 Characteristics of Geoids at simulation area mGal

3.2 試驗(yàn)步驟

為了模擬重力異常的觀測(cè)誤差，在仿真的重力異常Δg山、Δg丘和Δg海引入3種觀測(cè)誤差分別是零均值的白噪聲：e1（σ＝±1mGal）、e2（σ＝±3mGal）、e3（σ＝±5mGal）。以山區(qū)重力異常為例，試驗(yàn)步驟如下：

圖1 仿真的重力異常Fig.1 Simulative gravity anomalies

圖2 仿真的大地水準(zhǔn)面Fig.2 Simulative geoids

表3 各誤差條件下協(xié)方差函數(shù)的參數(shù)Tab.3 Parameters of the covariance function in three kinds of errors

（4）在已知協(xié)方差函數(shù)的系數(shù)A和RB的基礎(chǔ)上，依據(jù)式（6）中計(jì)算重力異常與大地水準(zhǔn)面的協(xié)方差矩陣CNΔg和階方差矩陣CΔgΔg。

（5）依據(jù)Tikhonov正則化原理，利用GCV法計(jì)算階方差矩陣CΔgΔg求逆時(shí)的正則化參數(shù)（見圖3）。

丘陵區(qū)域和海洋區(qū)域的仿真試驗(yàn)和山區(qū)區(qū)域的仿真試驗(yàn)類似。

3.3 結(jié)果比較和分析

圖3 各誤差條件下協(xié)方差矩陣CΔgΔg用GCV法確定正則化參數(shù)示意圖Fig.3 Regularization parameters chosen by the GCV method in three kinds of errors

為了驗(yàn)證本方法的效果，設(shè)計(jì)了3種計(jì)算方法。

方法1：Stokes法（積分半徑取1度）；

方法2：直接的最小二乘配置法；

方法3：聯(lián)合最小二乘配置法和Tikhonov正則化法的算法。

比較結(jié)果見表4。

由表4可以看出：

（1）在設(shè)計(jì)的3種誤差條件下，方法2計(jì)算的大地水準(zhǔn)面出現(xiàn)千米級(jí)誤差，表明選擇取全區(qū)域觀測(cè)值擬合協(xié)方差系數(shù)將增大協(xié)方差矩陣之間的相關(guān)性，使得矩陣嚴(yán)重病態(tài)。

（2）在設(shè)計(jì)的3種誤差條件下，方法3計(jì)算的大地水準(zhǔn)面的標(biāo)準(zhǔn)差，在山區(qū)區(qū)域分別為9.19cm、9.14cm和8.98cm；在丘陵區(qū)域分別為3.69cm、3.80cm和3.36cm；在海洋區(qū)域分別為2.14cm、1.98cm和1.96cm，遠(yuǎn)優(yōu)于方法2計(jì)算的大地水準(zhǔn)面，表明方法3能有效抑制觀測(cè)誤差對(duì)結(jié)果的影響，得到穩(wěn)定精確的大地水準(zhǔn)面高。

（3）方法1計(jì)算的大地水準(zhǔn)面的標(biāo)準(zhǔn)差，在山區(qū)區(qū)域分別為15.01cm、15.06cm和15.09cm；在丘陵區(qū)域分別為5.36cm、5.51cm和5.51cm；在海洋區(qū)域分別為2.93cm、3.14cm和3.19cm，與方法3比較，二者精度相當(dāng)。

（4）在誤差增大的情況下，方法3的正則化參數(shù)顯著增大，在山區(qū)區(qū)域分別為495.50、10 796.80和41 125.74；在丘陵區(qū)域分別為90.54、953.37和2 577.74；在海洋區(qū)域分別為1 867.72、7 219.84和38 015.51。對(duì)應(yīng)的CΔgΔg條件數(shù)依次減小，在山區(qū)區(qū)域分別為8 724.81、1 922.93和981.87；在丘陵區(qū)域分別為981.87、302.84和182.93；在海洋區(qū)域分別為302.84、154.63和66.75，同時(shí)大地水準(zhǔn)面的標(biāo)準(zhǔn)差也相應(yīng)減小，表明正則化參數(shù)嚴(yán)重影響方法3的穩(wěn)定性和精度。

表4 與大地水準(zhǔn)面N0的比較Tab.4 Comparison Results of Geoids Differences between calculated results and Simulative Geoid N0 cm

4 結(jié)束語

最小二乘配置法由于能融合不同種類重力觀測(cè)數(shù)據(jù)確定大地水準(zhǔn)面的特性而受到廣泛關(guān)注。但由于協(xié)方差矩陣存在病態(tài)性，微小的觀測(cè)誤差將被觀測(cè)數(shù)據(jù)協(xié)方差矩陣的小奇異值放大，導(dǎo)致計(jì)算的配置結(jié)果不穩(wěn)定且精確偏低。

在最小二乘配置法中引入Tikhonov正則化法，利用正則化參數(shù)修正協(xié)方差矩陣的小奇異值，能抑制其對(duì)觀測(cè)誤差的放大影響?；赥ikhonov－LSC法計(jì)算大地水準(zhǔn)面，能有效提高穩(wěn)定性和精度。通過以EGM2008重力場(chǎng)模型分別計(jì)算的山區(qū)、丘陵和海域重力異常作為基礎(chǔ)數(shù)據(jù)確定相應(yīng)區(qū)域大地水準(zhǔn)面的試驗(yàn)，驗(yàn)證了該方法的有效性。

Tikhonov正則化法能有效改進(jìn)基于重力異常利用最小二乘配置理論計(jì)算大地水準(zhǔn)面的精度和穩(wěn)定性。但Tikhonov正則化法是否適用于最小二乘配置理論的最優(yōu)正則化法，還有待進(jìn)一步研究。

［1］ CHAO Dingbo，SHEN Wenbin，WANG Zhengtao.Investigation of the Possibility and Method of Determining Global Centimeter－level Geoid［J］.Acta Geodaetica et Cartographica Sinica，2007，36（4）：370－376.（晁定波，申文斌，王正濤.確定全球厘米級(jí)大地水準(zhǔn)面的可能性與方法探討［J］.測(cè)繪學(xué)報(bào)，2007，36（4）：370－376.）

［2］ XU Xi，ZHU Jianjun.Relative Accuracy Estimation for Determining Regional Gravimetric Geoid［J］.Acta Geodaetica et Cartographica Sinica，2009，38（5）：383－390.（許曦，朱建軍.區(qū)域重力大地水準(zhǔn)面確定的相對(duì)精度估計(jì)［J］.測(cè)繪學(xué)報(bào)，2009，38（5）：383－390.）

［3］ LIU Zhenyu，GAO Binghao.Establishing the West Jilin Quasi－geoid Based on CQG2000［J］.Acta Geodaetica et Cartographica Sinica，2010，39（5）：441－443.（劉振宇，高炳浩.基于CQG2000的吉林省西部地區(qū)似大地水準(zhǔn)面的建立［J］.測(cè)繪學(xué)報(bào)，2010，39（5）：441－443.）

［4］ MORITZ H.Advanced Physical Geodesy［M］.Karlsruhe：Wichmann，1980.

［5］ HEISKANEN W A，MORITZ H.Physical Geodesy［M］.San Francisco：Freeman and Company，1967.

［6］ BJERHAMMAR A.A New Theory of Geodetic Gravity［R］.Stockholm：Freeman and Company，1967.

［7］ WU Xing，ZHANG Chuanding，LIU Xiaogang.Least－squares Collocation Harmonic Analysis of the Radial Satellite Gravity Gradients［J］.Acta Geodaetica et Cartographica Sinica，2010，39（5）：471－477.（吳星，張傳定，劉曉剛.衛(wèi)星重力徑向梯度數(shù)據(jù)的最小二乘配置調(diào)和分析［J］.測(cè)繪學(xué)報(bào)，2010，39（5）：471－477.）

［8］ CHAI Hongzhou，CUI Yue，MING Feng.The Determination of Chinese Mainland Crustal Movement Model Using Leastsquares Collocation［J］.Acta Geodaetica et Cartographica Sinica，2009，38（1）：61－65.（柴洪洲，崔岳，明鋒.最小二乘配置方法確定中國大陸主要塊體運(yùn)動(dòng)模型［J］.測(cè)繪學(xué)報(bào)，2009，38（1）：61－65.）

［9］ ZHANG Chuanyin，DING Jian，CHAO Dingbo.General Expression of Least Squares Collocation in Gravity Field［J］.Geomatics and Information Science of Wuhan University，2007，32（5）：431－434.（章傳銀，丁劍，晁定波.重力場(chǎng)最小二乘配置通用表示技術(shù)［J］.武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，2007，32（5）：431－434.）

［10］ ZOU Xiancai，LI Jiancheng.A Geoid Determination Using Least－Squares Collocation［J］.Geomatics and Information Science of Wuhan University，2004，29（3）：218－222.（鄒賢才，李建成.最小二乘配置方法確定大地水準(zhǔn)面的研究［J］.武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，2004，29（3）：218－222.）

［11］ TSCHERNING C C，RAPP R H.Closed Covariance Expressions for Gravity Anomalies，Geoid Undulations and Deflections of the Vertical Implied by Anomaly Degree Variance Model［R］.Ohio：Ohio State University，1974.

［12］ KNUDSEN P.Estimation and Modelling of the Local Empirical Covariance Function Using Gravity and Satellite Altimeter Data［J］.Bulletin Geodesique，1987，61：145－160.

［13］ WEN Hanjiang.The Estimation of Covariance Function in Least Squares Collocation［J］.Science of Surveying and Mapping，2000，25（3）：37－40.（文漢江.最小二乘配置法中協(xié)方差函數(shù)的計(jì)算［J］.測(cè)繪科學(xué)，2000，25（3）：37－40.）

［14］ AYHAN M E.Geoid Determination in Turkey （TG－91）［J］.Journal of Geodesy，1993，67（1）：10－22.

［15］ ARABELOS D N，TSCHERNING C C.Error－Covariance of the Estimates of Spherical Harmonic Coefficients Computed by LSC，Using Second－order Radial Derivative Functional Associated with Realistic GOCE Orbits［J］.Journal of Geodesy，2009，83（5）：419－430.

［16］ KOTSAKIS C.Least－Squares Collocation with Covariancematching Constraints［J］.Journal of Geodesy，2007，81（10）：661－677.

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

［18］ LU Zhonglian.Theory and Method of Earth’s Gravity Field［M］.Beijing：PLA Publishing House，1996.（陸仲連.地球重力場(chǎng)的理論與方法［M］.北京：解放軍出版社，1996.）

［19］ 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.（王振杰.大地測(cè)量中不適定問題的正則化解法研究［D］.武漢：中國科學(xué)院測(cè)量與地球物理研究所，2003.）

［20］ DENG Kailiang，BAO Jingyang，HUANG Motao，et al.Simulation of Tikhonov Regulation Algorithm in Downward Continuation of Airborne Gravity Data［J］.Geomatics and Information Science of Wuhan University，2010，35（12）：1414－1417.（鄧凱亮，暴景陽，黃謨濤，等.航空重力數(shù)據(jù)向下延拓的Tikhonov正則化法仿真研究［J］.武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，2010，35（12）：1414－1417.）

［21］ 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］.大地測(cè)量學(xué)與地球動(dòng)力學(xué)，2002，33（3）：11－14.）

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

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

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

［25］ PAVLIS N K，HOLMES S A，KENYON S C，et al.An Earth Gravitational Model to Degree 2160：EGM2008［C］∥Proceedings of the 2008General Assembly of the European Geosciences Union，Vienna：［s.n.］，2008.

［26］ ZHANG Chuanying，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地球重力場(chǎng)模型在中國大陸適用性分析［J］.測(cè)繪學(xué)報(bào)，2009，38（4）：283－289.）