劉曉剛，吳曉平，王 凱

1 信息工程大學測繪學院，鄭州 450052

2 西安測繪研究所，西安 710054

擾動重力梯度張量單分量和組合分量最小二乘配置法模型的建立

劉曉剛1，2，吳曉平1，王 凱1

1 信息工程大學測繪學院，鄭州 450052

2 西安測繪研究所，西安 710054

建立了利用擾動重力梯度張量Tzz分量和Txx＋Tyy、Tzz－Txx－Tyy組合分量確定地球重力場的調和分析法模型，進一步推導了擾動重力梯度張量對角線三分量的自協(xié)方差和互協(xié)方差函數的級數展開式，推導了單分量、組合分量與重力位系數之間協(xié)方差函數的實用計算公式，給出了利用單分量和組合分量解算地球重力場模型的最小二乘配置法基本原理公式.結果表明，最小二乘配置法具有一定的抗差能力，隨著觀測數據誤差的不斷增大，其恢復的重力場模型有效階次不斷降低，精度也不斷下降；Tzz－Txx－Tyy組合分量解算重力場模型的精度最高，其次為Tzz分量，Txx＋Tyy組合分量最差.

地球重力場模型，調和分析法，最小二乘配置法，衛(wèi)星重力梯度，GOCE

1 引 言

最小二乘配置法（LSC：Least Squares Collocation）是近代地球重力場逼近中非常重要的一項研究成果.特別是隨著GOCE衛(wèi)星測量計劃的提出和實施，該方法作為確定全球或局部區(qū)域地球重力場模型的方法之一也備受關注，其理論也得到了進一步的發(fā)展和完善.

Schwarz和Krynski［1］研究了聯合衛(wèi)星重力梯度（SGG：Satellite Gravity Gradient）數據和陸地重力測量數據精化局部大地水準面的LSC方法；Tscherning［2］研究了最小二乘法和配置法在從空間測量數據中提取重力場信息以及將其與地面重力測量數據進行融合中的應用，并指出兩種方法在觀測量和參數個數相等時是等價的；羅志才［3］提出的頻域LSC法可以有效提高計算速度；Li和Sideris［4］聯合采用衛(wèi)星測高獲取大地水準面高數據和船測重力異常數據，對LSC法、頻域最小二乘法以及多輸入／單輸出系統(tǒng)理論恢復局部地球重力場的精度進行了比較；張傳定［5］完善了LSC理論，提出了最小二乘復配置理論，得到了各類單定邊值問題的解析解、調和分析解、最小二乘解和最小二乘復配置解，建立了超定邊值問題的最小方差解、最小二乘解和最小二乘復配置解；Tscherning［6］研究了利用LSC法估計全球重力場位系數及其誤差估計的方法；Tscherning等［7］對僅采用擾動引力數據時由FFT（Fast Fourier Transform）方法確定的局部大地水準面，與聯合擾動引力和GPS／水準數據由LSC方法確定的局部大地水準面精度進行了比較；Sansò和Tscherning［8］提出了快速球諧配置法，采用LSC法聯合Tzz數據和地面重力異常數據進行了模擬試算；李迎春［9］利用引力位的協(xié)方差函數導出了位系數與引力位二階導數的協(xié)方差函數模型，并推導出LSC法的實用解算公式；張嗥、陳瓊［10］通過空間擾動位協(xié)方差函數特性，得出SGG數據與擾動位系數的相關協(xié)方差函數，利用LSC法，推導出由SGG數據直接解算引力位系數的函數表達式；汪海洪等［11］對不同分辨率的數據融合問題進行了初步研究，詳細闡述了多分辯率LSC法的基本原理，推導了該方法的具體公式；Kotsakis［12］研究了在盡量不改變精度的前提下，通過消除LSC法內在的平滑效應來增強其恢復地球重力場能力的方法；章傳銀等［13］在分析局部重力場LSC法技術特點的基礎上，推導出一種能綜合多種類型、不同高度重力場元經驗協(xié)方差函數的通用表達方法；吳星等［14］推導了擾動引力梯度數據與球諧系數之間的協(xié)方差陣和自協(xié)方差陣，并有效利用FFT技術將其降階，得到了引力梯度徑向分量LSC法的完整計算公式.

本文首先建立了利用擾動重力梯度張量Tzz分量和Txx＋Tyy、Tzz－Txx－Tyy組合分量解算地球重力場的調和分析法模型，進一步推導了擾動重力梯度張量對角線三分量之間的自協(xié)方差和互協(xié)方差函數的級數展開式，推導了單分量、組合分量與重力位系數之間協(xié)方差函數的計算公式，給出了利用單分量和組合分量解算地球重力場模型的LSC法基本原理公式，最后對本文建立的LSC法模型進行了數值實驗驗證.

2 單分量和組合分量的調和分析法模型

擾動重力梯度張量對角線三分量的地心球坐標表達式為［15］其中，fM為引力常數；r為衛(wèi)星軌道上任意一點的地心向徑；θ和λ分別為它的地心余緯和地心經度；R為地球平均半徑為完全正?；厍驍_動引力位系數；n和m分別為球諧系數的階和次；nm（cosθ）為完全正?；喓侠兆尩潞瘮?

將（2）式和（3）式相加，可以得到Txx＋Tyy組合分量的表達式：

利用（1）式減去（4）式，可以得到Tzz－Txx－Tyy組合分量的表達式：

同理，可得利用Txx＋Tyy、Tzz－Txx－Tyy組合分量解算位系數和的調和分析法公式分別為

上述Tzz分量以及Txx＋Tyy、Tzz－Txx－Tyy組合分量的調和分析法模型采用的都是格網中點值，如果采用格網平均值，則需要加入平滑因子，關于平滑因子的計算，請參閱文獻［17］.

3 單分量和組合分量的最小二乘配置法模型

在研究地球重力場的有關問題時，擾動位的協(xié)方差是最基本的，其它重力場元（重力異常、高程異常、垂線偏差和重力梯度張量等）之間的協(xié)方差都可以由它求得.空間擾動位協(xié)方差函數的級數展開式為［5，18］

其中

式中Kn在一定階次N以內的值，通常是利用參考重力場模型系數來計算，高于N以上階次的值一般利用重力異常的階方差模型求得，其關系為

如果采用Moritz提出的兩分量模型，其形式為［19］

其中，α1，α2，S1，S2為實常數；A，B為整常數.取

α1，α2的單位是（mGal）2，其它量是無量綱的數.（Δg）不能由（12）式計算，一般采用（Δg）＝7.5（mGal）2.

根據擾動重力梯度張量對角線三分量與擾動位T之間的關系式（1）—（3），并考慮到擾動位協(xié)方差函數的級數展開式（9），可以推導出擾動重力梯度張量對角線三分量的自協(xié)方差和互協(xié)方差函數級數展開式為

由擾動重力梯度張量Tzz分量解算位系數的調和分析法公式（6），考慮到Tzz分量的協(xié)方差函數級數展開式（14），調換求和與積分次序，并根據球諧函數的正交性，可以得到Tzz分量與位系數的協(xié)方差公式為

根據Txx＋Tyy組合分量解算位系數的調和分析法模型（7）式，考慮到Txx和Tyy分量的自協(xié)方差以及互協(xié)方差函數級數展開式（15）、（16）、（19），調換求和與積分次序，并根據球諧函數的正交性，可以得到Txx＋Tyy組合分量與位系數的協(xié)方差公式為

根據Tzz－Txx－Tyy組合分量解算位系數的調和分析法模型（8）式，考慮到Txx、Tyy和Tzz分量的自協(xié)方差以及互協(xié)方差函數級數展開式（14）—（19），調換求和與積分次序，并根據球諧函數的正交性，可以得到Tzz－Txx－Tyy組合分量與位系數的協(xié)方差公式為

下面給出利用擾動重力梯度張量Tzz分量以及Txx＋Tyy、Tzz－Txx－Tyy組合分量解算地球重力場模型的LSC法的基本原理公式［3，5，9，14］：

式中，J表示位系數的列向量.其它幾個量表述如下：

（1）單分量Tzz

對于單分量Tzz，則（23）式中表示位系數與Tzz分量的協(xié)方差矩陣，由（20）式進行計算；表示Tzz分量之間的協(xié)方差矩陣，由（14）式進行計算；Tij表示全球規(guī)則平均擾動重力梯度張量Tzz分量形成的列向量.

（2）組合分量Txx＋Tyy

對于組合分量Txx＋Tyy，則（23）式中表示位系數與Txx＋Tyy組合分量的協(xié)方差矩陣，由（21）式進行計算；表示Txx＋Tyy組合分量之間的協(xié)方差矩陣，由（15）、（16）、（19）式進行計算；Tij表示全球規(guī)則平均擾動重力梯度張量Txx＋Tyy組合分量形成的列向量.

（3）組合分量Tzz－Txx－Tyy

對于組合分量Tzz－Txx－Tyy，則（23）式中表示位系數與Tzz－Txx－Tyy組合分量的協(xié)方差矩陣，由（22）式進行計算，而表示Tzz－Txx－Tyy組合分量之間的協(xié)方差矩陣，由（14）—（19）式進行計算；Tij表示全球規(guī)則平均擾動重力梯度張量Tzz－Txx－Tyy組合分量形成的列向量.

4 數值實驗與結果分析

在具體數值實驗中，首先利用EGM2008重力場模型生成了GOCE衛(wèi)星軌道上的重力梯度數據［20］，通過延拓［21］和格網化方法［22］將重力梯度數據歸算為250km平均軌道高度處45′×45′的格網平均值，并對兩極數據空白區(qū)進行填充，得到全球范圍內的重力梯度格網平均值數據；向45′×45′的全球擾動重力梯度格網平均值數據中分別加入標準差為0、0.1、1×10－12s－2和10×10－12s－2的零均值白噪聲，通過LSC法，采用FFT技術解算地球重力場模型，來測試不同的數據噪聲對重力場模型解算結果的影響，其對應的重力場模型階誤差RMS、累計大地水準面和重力異常誤差如圖1所示，表1表示加入不同誤差的Tzz數據恢復120階重力場模型對應的精度統(tǒng)計結果.

其次，利用含有0.5×10－12s－2零均值白噪聲的45′×45′全球擾動重力梯度格網平均值數據，對本文推導出的單分量Tzz和組合分量Txx＋Tyy、Tzz－Txx－Tyy的LSC法模型的有效性進行測試，其對應的重力場模型階誤差RMS、累計大地水準面和重力異常誤差如圖2所示，表2表示單分量Tzz和組合分量Txx＋Tyy、Tzz－Txx－Tyy數據恢復120階重力場模型對應的精度統(tǒng)計結果.

從圖1和表1可以看出，當觀測數據含有不同的誤差時對解算結果的精度影響較大，隨著加入誤差的不斷增大，其恢復的重力場模型有效階次不斷降低，精度也不斷下降.在Tzz數據中分別加入0.1×10－12s－2、1×10－12s－2和10×10－12s－2的零均值白噪聲時，其恢復的重力場模型有效階數分別約為140、135和110階.含有10×10－12s－2誤差的數據對最小二乘配置法解算結果的精度影響較大，而含有0.1×10－12s－2和1×10－12s－2誤差的數據，其解算結果在2～240階次內與不含誤差的數據結果比較一致，這說明最小二乘配置法在解算含有誤差的數據時具有一定的優(yōu)勢.其中，加入0.1×10－12s－2的誤差數據恢復的120階重力場模型，其對應的累計大地水準面誤差達到5.938cm、累計重力異常誤差達到1.095mGal的精度.

從圖2和表2可以看出，Tzz－Txx－Tyy組合分量解算結果的精度最高，其次為Tzz分量和Txx＋Tyy組合分量，這是因為擾動重力梯度張量的對角線三分量中，Tzz分量數據的量級最大，精度也最高，Txx和Tyy分量次之，在模型解算時，Tzz－Txx－Tyy組合分量不僅顧及了對角線三分量Tzz、Txx、Tyy的整體影響，而且各分量之間求差也在一定程度上削弱了數據中所含噪聲對解算結果的影響；Tzz分量的解算結果優(yōu)于Txx＋Tyy組合分量的原因是，在對角線三分量中，Tzz分量對模型解算的貢獻最大，另外，Txx＋Tyy的組合在一定程度上放大了數據中所含噪聲對解算結果的影響.單分量Tzz和組合分量Txx＋Tyy、Tzz－Txx－Tyy恢復的120階重力場模型對應的累計大地水準面誤差分別達到5.994cm、6.111cm、5.944cm，累計重力異常誤差分別達到1.105mGal、1.127mGal、1.096mGal的精度.

表1 加入不同誤差的Tzz數據恢復120階重力場模型的精度統(tǒng)計結果Table 1 The precision statistical results of the EGM in 120degree recovered by Tzzdata with different errors

表2 擾動重力梯度張量單分量和組合分量數據恢復120階重力場模型的精度統(tǒng)計結果Table 2 The precision statistical results of the EGM in 120degree recovered by single component and combined components of disturbing gravity gradients

5 結 語

LSC法需要將重力梯度數據進行歸算和格網化處理，得到平均軌道球面上的規(guī)則格網點或格網平均值數據，如果能預先確定格網的大小，則格網點的個數也就可以確定，從而法方程的維數也是一定的，并且由于該方法建立的法方程都是呈塊狀對角形式，因此方程的求逆等計算將非常簡便、快捷；另外，該方法也可以綜合各類觀測數據確定地球重力場模型，在數據處理時能顧及觀測量的誤差，并且能夠給出待估量的誤差協(xié)方差信息.它的缺點是在將觀測數據進行歸算和格網化處理時，不可避免地引入了一些歸算和格網化誤差，并且要求有可靠的先驗協(xié)方差模型，在恢復高階次地球重力場模型時，需要解算超大型協(xié)方差矩陣的逆矩陣，由于大型矩陣求逆存在的病態(tài)性問題，從而影響了數值解的穩(wěn)定性.

本文推導得到了利用擾動重力梯度張量Tzz分量以及Txx＋Tyy、Tzz－Txx－Tyy組合分量解算地球重力場的最小二乘配置法模型，進一步充實了利用重力梯度數據解算地球重力場模型的方法，從而使得最小二乘配置法模型的內容更加豐富.

（References）

［1］ Schwarz K P，Krynski J.Improvement of the geoid in local areas by satellite gradiometry.Bulletin Géodésique，1977，51：163－176.

［2］ Tscherning C C.Collocation and least squares methods as a tool for handling gravity field dependent data obtained through space research techniques.Bulletin Géodésique，1978，52：199－212.

［3］ 羅志才.利用衛(wèi)星重力梯度數據確定地球重力場的理論與方法［博士論文］.武漢：武漢測繪科技大學，1996.Luo Z C.Theories and method of the determination of the earth′s gravity field using satellite gravity gradients［Ph.D.thesis］（in Chinese）.Wuhan：Wuhan Technical University of Surveying and Mapping，1996.

［4］ Li J，Sideris M G.Marine gravity and geoid determination by optimal combination of satellite altimetry and shipborne gravimetry data.Journal of Geodesy，1997，71：209－216.

［5］ 張傳定.衛(wèi)星重力測量——基礎、模型化方法與數據處理算法［博士論文］.鄭州：信息工程大學測繪學院，2000.Zhang C D.Satellite gravimetry：foundation，modeling methods，and data processing algorithms［Ph.D.thesis］（in Chinese）.Zhengzhou：Institute of Surveying and Mapping，2000.

［6］ Tscherning C C.Computation of spherical harmonic coefficients and their error estimate using least－squares collocation.Journal of Geodesy，2001，75：12－18.

［7］ Tscherning C C，Radwan A，Tealeb A A，et al.Local geoid determination combining gravity disturbances and GPS－levelling—a case study in the Lake Nasser area，Aswan，Egypt.Journal of Geodesy，2001，75：343－348.

［8］ SansòF，Tscherning C C.Fast spherical collocation：theory and examples.Journal of Geodesy，2003，77：101－112.

［9］ 李迎春.利用衛(wèi)星重力梯度測量數據恢復地球重力場的理論與方法［碩士論文］.鄭州：信息工程大學測繪學院，2004.Li Y C.Theories and method of the recovery of the earth′s gravity field using satellite gravity gradients［Master′s thesis］（in Chinese）.Zhengzhou：Institute of Surveying and Mapping，2004.

［10］ 張皞，陳瓊.衛(wèi)星重力梯度數據解算位系數的最小二乘配置法.測繪與空間地理信息，2005，28（6）：34－35.Zhang H，Chen Q.Least squares collocation in calculating potential coefficient with satellite gravity gradiometry data.Geomatics and Spatial Information Technology（in Chinese），2005，28（6）：34－35.

［11］ 汪海洪，羅志才，羅佳，等.多分辨率最小二乘配置法初探.大地測量與地球動力學，2006，26（1）：115－118.Wang H H，Luo Z C，Luo J，et al.Preliminary research on multi－resolution least－square collocation.Journal of Geodesy and Geodynamics（in Chinese），2006，26（1）：115－118.

［12］ Kotsakis C.Least－squares collocation with covariancematching constrains.Journal of Geodesy，2007，81：661－677.

［13］ 章傳銀，丁劍，晁定波.局部重力場最小二乘配置通用表示技術.武漢大學學報（信息科學版），2007，32（5）：431－434.Zhang C Y，Ding J，Chao D B.General expression of least squares collocation in local gravity field.Geomatics and Information Science of Wuhan University（in Chinese），2007，32（5）：431－434.

［14］ 吳星，張傳定，劉曉剛.衛(wèi)星重力徑向梯度數據的最小二乘配置調和分析.測繪學報，2010，39（5）：471－477.Wu X，Zhang C D，Liu X G.Least squares collocation harmonic analysis of the radial satellite gravity gradients.Geodaetica et Cartographica Sinica（in Chinese），2010，39（5）：471－477.

［15］ Reed G B.Application of Kinematical Geodesy for Determining the Short Wave Length Components of the Gravity Field by Satellite Gradiometry.Ohio State University，Dept.of Geod.Sciences，Rep.No.201，Columbus，Ohio，1973.

［16］ 陸仲連.球諧函數.鄭州：測繪學院出版社，1984.Lu Z L.Spherical Harmonic Functions（in Chinese）.Zhengzhou：Institute of Surveying and Mapping Publishing House，1984.

［17］ 陸仲連.地球重力場理論與方法.北京：解放軍出版社，1996.Lu Z L.Theory and method of the earth′s gravity field（in Chinese）.Beijing：PLA Publishing House，1996.

［18］ 吳星.衛(wèi)星重力梯度數據處理理論與方法［博士論文］.鄭州：信息工程大學測繪學院，2009.Wu X.Research of Methods of Spherical Harmonic Analysis of the Earth’s Gravity Field［Ph.D.thesis］（in Chinese）.Zhengzhou：Institute of Surveying and Mapping，2009.

［19］ Moritz H.Advanced physical geodesy.England：Abacus Press，1980.

［20］ 劉曉剛，吳曉平，趙東明，等.擾動重力梯度的非奇異表示.測繪學報，2010，39（5）：450－457.Liu X G，Wu X P，Zhao D M，et al.Non－singular expression of the disturbing gravity gradients.Geodaetica et Cartographica Sinica（in Chinese），2010，39（5）：450－457.

［21］ 劉曉剛，李珊珊，吳星.衛(wèi)星重力梯度數據的向下延拓.大地測量與地球動力學，2011，31（1）：132－137.Liu X G，Li S S，Wu X.Downward continuation of the satellite gravity gradient data.Journal of Geodesy and Geodynamics（in Chinese），2011，31（1）：132－137.

［22］ 劉曉剛，吳曉平，王寶軍，等.衛(wèi)星重力梯度數據的網格化方法研究.大地測量與地球動力學，2010，30（6）：60－65.Liu X G，Wu X P，Wang B J，et al.Study on gridding methods of the satellite gravity gradient data.Journal of Geodesy and Geodynamics（in Chinese），2010，30（6）：60－65.

Construction of the least squares collocation models for single component and composite components of disturbed gravity gradients

LIU Xiao－Gang1，2，WU Xiao－Ping1，WANG Kai1
1 Institute of Surveying and Mapping，Information and Engineering University，Zhengzhou450052，China
2 Xi′an Research Institute of Surveying and Mapping，Xi′an 710054，China

This work constructs the spherical harmonic analysis models for single component Tzzand composite components Txx＋Tyy，Tzz－Txx－Tyyof disturbed gravity gradients to compute earth′s gravitational field model（EGM）.Besides，the series expansion formulae of variance and covariance functions of the diagonal components of disturbed gravity gradients are deduced.Then，the practical computational formulae of covariance functions between single component，composite components and the gravity potential coefficients are deduced.Finally，this work presents the fundamental formula of the least squares collocation（LSC）to compute EGM using the single component and composite components.The results show that the LSC is robust to a certain extent and the effective degree and precision of EGM decay gradually with the increase of errors.The precision of the EGM computed by the composite components Tzz－Txx－Tyyis the best，while that of the composite components Txx＋Tyyis the least.

Earth′s gravitational field model（EGM），Spherical harmonic analysis method，Least squares collocation（LSC），Satellite gravity gradient（SGG），GOCE

10.6038／j.issn.0001－5733.2012.05.015

P223

2011－09－25，2012－05－10收修定稿

國家自然科學基金（41174026，41104047，41174017），中國科學院研究生院地球科學學院博士后基金（2010046），信息工程大學博士學位論文創(chuàng)優(yōu)基金資助.

劉曉剛，男，1983年生，博士，助理研究員，主要從事衛(wèi)星重力測量研究.E－mail：liuxiaogang＿1949＠163.com

劉曉剛，吳曉平，王凱.擾動重力梯度張量單分量和組合分量最小二乘配置法模型的建立.地球物理學報，2012，55（5）：1572－1580，

10.6038／j.issn.0001－5733.2012.05.015.

Liu X G，Wu X P，Wang K.Construction of the least squares collocation models for single component and composite components of disturbed gravity gradients.Chinese J.Geophys.（in Chinese），2012，55（5）：1572－1580，doi：10.6038／j.issn.0001－5733.2012.05.015.

（本文編輯 汪海英）