李新星,吳曉平,李姍姍,劉曉剛,崔志偉

1．信息工程大學(xué)地理空間信息學(xué)院,河南鄭州 450001;2．西安測繪研究所,陜西西安 710054; 3．61206部隊,北京 100042

塊對角最小二乘方法在確定全球重力場模型中的應(yīng)用

李新星1,吳曉平1,李姍姍1,劉曉剛2,崔志偉3

1．信息工程大學(xué)地理空間信息學(xué)院,河南鄭州 450001;2．西安測繪研究所,陜西西安 710054; 3．61206部隊,北京 100042

分析研究最小二乘求解重力場模型的6種位系數(shù)排列方式下的塊對角形態(tài)及3種不同條件下的塊對角近似:BD-1、BD-2、BD-3,探討塊對角最小二乘方法在聯(lián)合早期衛(wèi)星重力場模型和最新GRACE-only模型中的應(yīng)用。試驗結(jié)果表明,塊對角最小二乘方法較之于積分方法,能更好地提高所恢復(fù)模型的精度,說明在衛(wèi)星重力飛速發(fā)展、地面重力數(shù)據(jù)不斷完善的今天,塊對角最小二乘法在超高階地球重力場模型構(gòu)建方面的優(yōu)勢逐漸突出。

塊對角矩陣;最小二乘;超高階地球重力場模型;EGM2008;正交性

1 引 言

自1937年首次推出6階地球重力位模型開始,經(jīng)過70多年的發(fā)展,地球重力場模型已經(jīng)發(fā)展到2160完全階次。在模型的構(gòu)建方法中,基于數(shù)值積分方法在對精度評估方面的缺陷,最小二乘方法已成為目前國際上構(gòu)建重力場模型的主要方法。2008年NGA(US National Geospatial-Intelligence Agency)發(fā)布的目前精度最好的超高階地球重力場模型EGM2008的解算采用最小二乘方法,在具有高質(zhì)量重力數(shù)據(jù)的地區(qū),EGM2008模型大地水準面與獨立的GPS/水準數(shù)據(jù)之間差異小于±10 cm,其確定的垂線偏差相對于獨立的天文大地垂線偏差的差異在±1．1″和±1．3″。我國由于全球數(shù)據(jù)的缺乏,目前所構(gòu)建的超高階地球重力場模型精度較差,且由于均采用數(shù)值積分方法,沒有給出合適的精度評價。另外,隨著衛(wèi)星重力的不斷發(fā)展,獲得的衛(wèi)星重力場模型不斷精化超高階地球重力場模型中低階部分,最小二乘方法在聯(lián)合地面重力數(shù)據(jù)和衛(wèi)星重力場方面也起到了至關(guān)重要的作用。所以,要構(gòu)建我國自己的高精度的超高階全球重力場模型,一方面要獲取全球高質(zhì)量的重力數(shù)據(jù),另一方面應(yīng)當發(fā)展當前普遍使用的最小二乘方法[1-6]。

重力場模型構(gòu)建過程中,由于觀測方程中的未知量(模型位系數(shù))太多,且觀測值數(shù)量巨大,所以在構(gòu)建法方程以及法方程求解過程中會遇到海量計算問題,直接求解是不可行的,而球諧函數(shù)的正交特性使得由觀測方程得到的法方程矩陣是一稀疏矩陣,那么,通過對未知數(shù)的重新排列,獲得具有塊對角形式的法矩陣,這對于法方程的解算具有重要作用[7-8]。本文重點討論了利用塊對角最小二乘方法構(gòu)建全球重力場模型以及聯(lián)合衛(wèi)星重力場模型的方法,并通過試驗計算,與數(shù)值積分方法進行了相關(guān)比較。

2 基本原理

全球格網(wǎng)平均重力異常采用球諧分析求解地球重力場模型的基本觀測方程為[9-17]

ˉrij為第i行第j列格網(wǎng)平均重力異常的地心向徑;λj表示第j列格網(wǎng)中點經(jīng)度;(θiN、θis)表示第i行格網(wǎng)的北邊界緯度與南邊界緯度;a為正常橢球的長半徑;Δˉgcij為經(jīng)過橢球改正、大氣改正、二階梯度項改正等歸算后滿足調(diào)和邊值問題的地面格網(wǎng)平均重力異常值;ˉC?nm和ˉSnm即為所求的一組重力場模型位系數(shù)。

將觀測方程式(1)表示為最小二乘方法中的通用表達式

根據(jù)上述公式能夠求得最終的位系數(shù),然而對于利用全球5′×5′格網(wǎng)重力異常估計2159階位系數(shù)(不包括C00),A矩陣大小為9 331 200× 4 665 596,法矩陣ATPA大小為4 665 596× 4 665 596,所以龐大的計算量目前是不可能完成的。顧及龐大的計算量,為了提高計算效率,有必要討論法方程在地面重力異常數(shù)據(jù)滿足以下條件時的一些特殊性質(zhì):

(1)數(shù)據(jù)分布于一旋轉(zhuǎn)曲面上(例如旋轉(zhuǎn)橢球面上)。

(2)格網(wǎng)數(shù)據(jù)覆蓋整個曲面,且經(jīng)度方向的分辨率一致。

(3)數(shù)據(jù)中誤差與經(jīng)度無關(guān),即數(shù)據(jù)的權(quán)重不依賴于經(jīng)度的大小。

(4)數(shù)據(jù)權(quán)重關(guān)于赤道對稱,即緯度φ與－φ的格網(wǎng)數(shù)據(jù)精度相同。

顯然,滿足以上條件的權(quán)矩陣是個對角陣。根據(jù)式(1),矩陣A的各個元素由式(12)、式(13)表示

這里要注意區(qū)分變量r和ˉrij,相應(yīng)的U＝ATPLb的元素為

理想情況下,假設(shè)N×2N個重力異常分布于半徑為a的球面上,且重力異常都有相同的標準差σ,那么權(quán)陣P＝σ2I,在不影響問題討論下令σ＝1以簡化推導(dǎo),γ＝GM/a2,則式(16)可簡化為

另外,全球格網(wǎng)數(shù)據(jù)分布是關(guān)于赤道對稱的,滿足θi＝π－θN－1－i,再根據(jù)締合Legendre函數(shù)的特性,有

法矩陣在計算數(shù)據(jù)滿足以上4個條件的情況下,呈稀疏矩陣,有助于提高計算效率。滿足式(26)條件下,通過對法方程中的未知系數(shù)進行重新排序,生成具有較好的結(jié)構(gòu)的塊對角矩陣,這種塊對角法矩陣稱為第1類塊對角近似BD-1 (block-diagonal-1)。文獻[11]中給出了6種排列方式,表1簡要描述了這6種排列方式。

表1 6種塊對角系數(shù)排序方式Tab．1 6 kinds of block-diagonal coefficients sorting patterns

上述排列方法中,對Ⅴ排列方式進行進一步細化,區(qū)分n－m為奇數(shù)偶數(shù)的情況,例如,對于N＝6的排列,排列順序記為C20、C40、C60、C30、C50、C21、C41、C61、C31、C51、S21、S41、S61、S31、S51、C22、C42、…。下面采用了全球8×16的格網(wǎng)平均重力異常計算N＝6情況下法矩陣N[]的情況,其中權(quán)矩陣P采用單位陣I,具體見圖1,圖中白色區(qū)域表示為0元素,黑色像素點為非零元素。

在實際計算過程中,數(shù)據(jù)經(jīng)過空白區(qū)的填補和數(shù)據(jù)解析延拓等處理后能近似滿足(1)、(2)條件,但是因為各個局部區(qū)域的測量精度都存在差異,所以一般情況下不滿足(3)、(4)條件。所以當上述4個條件只有(4)條件不滿足時,法矩陣

該塊對角形式稱為第2類近似塊對角(BD-2)結(jié)構(gòu),第3類近似塊對角結(jié)構(gòu)(BD-3)結(jié)構(gòu)的法矩陣

下面根據(jù)上述定義式計算Ⅴ排列下,3種近似塊對角結(jié)構(gòu)的法矩陣形式:

通過圖2能夠很清晰地看出,對于Ⅴ排列,重力異常數(shù)據(jù)滿足上述條件越多,其法矩陣中非零元素越少,“塊”越小,塊對角形式越簡單。說明在模型計算過程中,隨著全球重力異常數(shù)據(jù)全球覆蓋不斷完善,分辨率和精度不斷提高,法矩陣可近似為越來越簡單的形式,便于模型的快速求解。

當模型位系數(shù)最大階數(shù)為Nmax、采用Ⅴ排列情況下,3種BD結(jié)構(gòu)的特點統(tǒng)計如表2所示。

表2 3種塊對角近似的特點Tab．2 Characters of 3 kinds of block-diagonal approximate conditions

圖1 6種系數(shù)排列方式下的法矩陣稀疏狀態(tài)Fig．1 Sparse status of normal matrix in 6 kinds of sorting patterns

圖2 3種近似條件下BD-1、SD-2、BD-3的法矩陣特點Fig．2 Normal matrix characters of BD-1,BD-2 and BD-3 in 3 kinds of approximate conditions

根據(jù)上述理論,相比數(shù)值積分方法,塊對角最小二乘具有以下特點:①數(shù)值積分方法確定的每個系數(shù)彼此之間是相互獨立的,而對于塊對角最小二乘而言,即使使用的重力異常數(shù)據(jù)之間不相關(guān),但確定的系數(shù)之間也是相關(guān)的;②數(shù)值積分方法不能對重力異常數(shù)據(jù)進行精度評估,塊對角最小二乘方法可以;③使用30′×30′離散數(shù)據(jù),根據(jù)Nyquist法則,數(shù)值積分方法能夠恢復(fù)到360階,而塊對角最小二乘只能最大恢復(fù)到359階。

3 塊對角聯(lián)合衛(wèi)星重力場模型

大多數(shù)模型求解中均使用最小二乘方法聯(lián)合衛(wèi)星重力場模型得到最終的聯(lián)合解。隨著衛(wèi)星重力場模型以及地面重力觀測數(shù)據(jù)的不斷完善,塊對角方法在聯(lián)合衛(wèi)星重力場模型方面也在發(fā)生著變化。

根據(jù)衛(wèi)星重力場模型及其誤差協(xié)方差矩陣,能夠完全再現(xiàn)計算該模型的法方程和對應(yīng)的U,誤差協(xié)方差矩陣的逆對應(yīng)于法矩陣

同樣根據(jù)地面數(shù)據(jù),能夠得到地面的法方程

如果使用簡單的等權(quán)相加來聯(lián)合上面兩個法方程,即將衛(wèi)星的法矩陣疊加到地面法矩陣上,就能夠得到聯(lián)合后的法方程

上式求得的未知數(shù)^x就是概念上的聯(lián)合解,在實際解算過程中,必須考慮兩種法方程聯(lián)合中的權(quán)重分配問題,該問題本文暫不討論。

3．1 早期衛(wèi)星重力場模型的聯(lián)合

在CHAMP(challenging mini-satellite payload)、GRACE(gravity recovery and climate experiment)等重力衛(wèi)星發(fā)射之前,早期的衛(wèi)星重力場模型是綜合衛(wèi)星軌道跟蹤、衛(wèi)星激光測距等手段求得的。文獻[10]中計算分析了使用地面數(shù)據(jù)計算的法矩陣和衛(wèi)星重力場模型的協(xié)方差矩陣(即法矩陣)中非塊對角元素的大小及其影響,結(jié)果顯示,在不損失結(jié)果精度的情況下,很難將當時的衛(wèi)星重力場模型的法矩陣近似為上述任何3種塊對角(BD)型,因此在前期模型構(gòu)建過程中使用的衛(wèi)星重力場模型協(xié)方差陣,沒有稀疏特性,但因為其階數(shù)較低,計算能夠?qū)崿F(xiàn)。另外,由于地面數(shù)據(jù)的精度和分布較差,所以地面觀測數(shù)據(jù)的法矩陣采用塊對角近似條件最弱的BD-3。

聯(lián)合衛(wèi)星重力場模型的計算過程中,使用了文獻[10]中所謂的“falling kite”結(jié)構(gòu),見圖3(d)。本文利用全球8×16格網(wǎng)數(shù)據(jù)模擬計算地面階數(shù)Dmax＝6的法矩陣,并聯(lián)合階數(shù)Dsat＝4的衛(wèi)星法矩陣。圖3(a)是地面法矩陣BD-3結(jié)構(gòu),與之對應(yīng)的衛(wèi)星系數(shù)也按Ⅴ排列,則其占滿的法矩陣的重新分布情況見圖3(b)黑色部分,顯然這種聯(lián)合后的法矩陣會產(chǎn)生一個很大的“塊”,幾乎是整個法矩陣的大小,所以該排列形式不利于高效的計算。

圖3 最小二乘中法矩陣的聯(lián)合Fig．3 Combination of normal matrix in LS combination solution

為構(gòu)造最小的“塊”,將要求解的未知數(shù)Cnm、Snm分為3組:①n＞Dsat,m＞Dsat;②n＞Dsat, m≤Dsat;③n≤Dsat,m≤Dsat。每一組中的未知數(shù)都按之前的Ⅴ方式排列,得到的地面法矩陣的結(jié)構(gòu)見圖3(c),衛(wèi)星法矩陣與地面聯(lián)合的法矩陣見圖3(d)。

對于“falling kite”結(jié)構(gòu)的計算,采用分組方法,考慮到法矩陣為對稱矩陣,將其按照圖4分割為若干個小矩陣,其中G11和G22是“純地面”的塊對角陣,G23由非零元素的矩形塊構(gòu)成,與G22有相同的行分割,G33則是包含地面和衛(wèi)星信息的對稱矩陣。

圖4 聯(lián)合法方程矩陣求解示意圖Fig．4 Sketch map of combination with normal matrix

根據(jù)式(32),^X1可單獨利用塊對角G11求解而不影響其他系數(shù),剩下的部分為

上式系數(shù)矩陣G為對稱矩陣,可對其進行Cholesky分解來快速穩(wěn)定的計算。由于最終模型的協(xié)方差陣的計算非常復(fù)雜,故在實際解算過程中,只計算以下3個內(nèi)容:①1系數(shù)的塊對角協(xié)方差矩陣G;②2每個系數(shù)的方差;③3系數(shù)的滿協(xié)方差矩陣[9]。

3．2 GRACE-only衛(wèi)星重力場模型的聯(lián)合

早期模型在聯(lián)合過程中使用了“falling kite”結(jié)構(gòu)的主要原因是因為在沒有損失結(jié)果質(zhì)量情況下很難將衛(wèi)星模型系數(shù)的法方程近似為任何塊對角型。而在最新超高階重力場模型的聯(lián)合中,因為使用的全球數(shù)據(jù)分辨率及精度高,全球分布較完善,所以在其求解過程中使用的是條件最嚴格的BD-1塊對角近似,另外,隨著重力衛(wèi)星技術(shù)的不斷發(fā)展,高精度的GRACE-only模型也可以獲取。GRACE信息的誤差特性允許其協(xié)方差矩陣采用塊對角近似,因此,在聯(lián)合過程中將GRACE-only模型的協(xié)方差矩陣近似為BD-1,與地面法方程疊加,而舍棄之前的BD-3結(jié)構(gòu)和“falling kite”結(jié)構(gòu),在不損害結(jié)果質(zhì)量下使得聯(lián)合解算既簡潔又高效。此外,由于橢球諧分析比球諧分析精度好,所以模型構(gòu)建均采用橢球諧分析,而在球諧系數(shù)和橢球諧系數(shù)轉(zhuǎn)換過程中,并沒有改變法矩陣的塊對角形式,這一點對于能夠采用塊對角方法來計算模型是至關(guān)重要的[18-20]。

將衛(wèi)星重力場模型Dsat＝4的BD-1法矩陣結(jié)構(gòu)疊加到地面數(shù)據(jù)確定的Dmax＝6的BD-1法矩陣上,對應(yīng)U矩陣的疊加與N矩陣對應(yīng),見圖5。

圖5 GRACE-only法方程聯(lián)合情況Fig．5 Sketch map of GRACE-only normal equation combination

根據(jù)圖5,顯然可以一次性聯(lián)合解一個對角塊,這樣,對于2160階次的超高階地球重力場模型,所需要求逆的最大的對稱矩陣是1080×1080規(guī)模的,這對于目前的計算水平來說并不算什么,這種近似相比于早期模型聯(lián)合衛(wèi)星模型,大大簡化了平差難度。

4 試驗計算分析

為驗證塊對角最小二乘在構(gòu)建全球重力場模型的優(yōu)勢及其實用性,本文使用3．2 GHz主頻、2 GB內(nèi)存計算機,基于vs2008試驗平臺,進行模擬計算。使用塊對角方法計算模型過程中,矩陣A最大,計算ATPA最耗時間,當計算一個2160階的超高階地球重力場模型,計算法矩陣中最大的“塊”所使用的A的大小為(2160－1)(2160＋3)×1080,開辟該數(shù)組需要至少37．6 GB內(nèi)存,所以基于硬件限制,本文模擬計算360階模型。首先使用EGM2008模型截至360階的模型系數(shù)模擬生成全球30′×30′格網(wǎng)平均重力異常,因為對于30′×30′分辨率數(shù)據(jù),最小二乘方法只能恢復(fù)到359階,所以使用模擬生成的格網(wǎng)平均重力異常采用數(shù)值積分以及塊對角最小二乘兩種方法恢復(fù)EGM2008模型前359階位系數(shù),結(jié)果比對見圖6。

圖6 模型恢復(fù)誤差的比較Fig．6 Comparison of error of recovered model

圖6中通過統(tǒng)計恢復(fù)系數(shù)與EGM2008系數(shù)差的絕對值的自然對數(shù),獲得數(shù)值積分方法和塊對角最小二乘方法恢復(fù)模型位系數(shù)的精度情況。圖6中明顯看出塊對角最小二乘方法恢復(fù)的模型整體精度比數(shù)值積分方法高出10個數(shù)量級,唯獨在偶階帶諧項C2k,0有較大誤差。

由基本原理可知,數(shù)值積分方法求解的位系數(shù)之間是相互獨立的,而塊對角最小二乘方法求得的位系數(shù)是相關(guān)的,因此在使用塊對角最小二乘方法求解模型位系數(shù)過程中,必須根據(jù)Nyquist法則對重力異常數(shù)據(jù)的頻譜能量進行限制,即無論使用多大分辨率數(shù)據(jù)恢復(fù)n階位系數(shù),每一離散數(shù)據(jù)中必須只包含0～n階系數(shù)的貢獻。下面以30′×30′格網(wǎng)平均重力異常恢復(fù)240階位系數(shù)為例,進行分析計算。

首先使用EGM2008模型截至360階的模型系數(shù)模擬生成全球30′×30′格網(wǎng)平均重力異常Δ0～360,采用數(shù)值積分方法恢復(fù)0～360階模型位系數(shù),利用恢復(fù)系數(shù)中的241～360階計算包含241～360階頻譜能量的30′×30′格網(wǎng)平均重力計算包含0～240階頻譜能量的30′×30′格網(wǎng)平均重力異常,利用該30′×30′格網(wǎng)重力異常Δ0～240,采用塊對角最小二乘方法恢復(fù)0～240階模型系數(shù)。使用EGM2008作為標準模型,計算兩種方法恢復(fù)的360階和240階模型系數(shù)的誤差階RMS并進行比較,比較結(jié)果見圖7,其計算公式為

圖7 誤差階RMS比較Fig．7 Comparison of error degree RMS

由圖7可看出,無論恢復(fù)359階還是240階模型位系數(shù),塊對角最小二乘方法恢復(fù)模型系數(shù)的誤差階RMS明顯小于使用數(shù)值積分方法恢復(fù)的模型系數(shù)。其中塊對角恢復(fù)的359階位系數(shù)的誤差階RMS曲線分為兩條,上面的是偶數(shù)階,下面是奇數(shù)階,其之間的差異由偶階帶諧項C2k,0的計算誤差引起的,而恢復(fù)的240階系數(shù)由于奇階和偶階系數(shù)之間的相關(guān)性發(fā)生頻譜疊效應(yīng),所以誤差階RMS曲線為一條。通過圖7比較分析可知,塊對角最小二乘相比數(shù)值積分方法,能夠更好地恢復(fù)模型,再加上其聯(lián)合衛(wèi)星重力場模型簡潔高效的特點,越來越凸顯其在模型構(gòu)建中的優(yōu)勢。

5 結(jié) 論

相比于數(shù)值積分方法,基于塊對角最小二乘構(gòu)建的全球重力場模型的優(yōu)勢主要體現(xiàn)在:①聯(lián)合中低階衛(wèi)星重力場模型更加簡潔高效;②能夠給出所構(gòu)建的全球重力場模型每階位系數(shù)的精度評價,相比中低階部分采用最小二乘平差方法、高階部分采用數(shù)值積分方法構(gòu)建的模型,塊對角最小二乘方法避免了誤差譜曲線在計算方法變化的階數(shù)部分不連續(xù)的現(xiàn)象;③模型恢復(fù)的精度更高?；谝陨蟽?yōu)勢,加上目前的計算能力已經(jīng)能夠滿足使用塊對角最小二乘方法求解2160階超高階地球重力場模型,因此,有必要進一步加強塊對角最小二乘方法計算模型方面的研究,這對于我國構(gòu)建自主高精度超高階地球重力場模型,提高我國區(qū)域大地水準面的精度具有重要作用。

[1] PAVLIS N K,HOLMES S A,KENYON S C,et al．The Development and Evaluation of the Earth Gravitational Model 2008[J]．Journal of Geophysical Research,2012, 117(04406)．DOI:10．1029/2011JB008916．

[2] WANG Zhengtao,DANG Yamin,CHAO Dingbo,et al．Theory and Methodology of Ultra-high-degree Geopotential Model Determination[M]Beijing:Surveying and Mapping Press,2011．(王正濤,黨亞民,晁定波．超高階地球重力場模型確定的理論與方法[M]．北京:測繪出版社,2011．)

[3] ZHENG Wei,XU Houze,ZHONG Min,et al．Progress and Present Status of Research on Earth’s Gravitational Field Models[J]．Journal of Geodesy and Geodynamics, 2009,30(4):83-90．(鄭偉,許厚澤,鐘敏,等．地球重力場模型研究進展和現(xiàn)狀[J]．大地測量與地球動力學(xué)．2009, 30(4):83-90．)

[4] NING Jinsheng．The Progress in the Earth’s Gravity Field in China[J]．Northeast Surveying and Mapping,2002,25 (4):6-9．(寧津生．我國地球重力場研究的進展[J]．東北測繪,2002,25(4):6-9．)

[5] SHU Chanfang,LI Fei,HAO Weifeng．Evaluation of EGM 2008 and Its Application Analysis over a Particular Region of China[J]．Geomatics and Information Science of Wuhan University,2011,36(8):919-922．(束蟬方,李斐,郝衛(wèi)峰．EGM2008模型在中國某地區(qū)的檢核及適用性分析[J]．武漢大學(xué)學(xué)報:信息科學(xué)版,2011,36(8):919-922．)

[6] ZHANG Chuanyin,GUO Chunxi,CHEN Junyong,et al．EGM 2008 and 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．)

[7] NING Jinsheng,QIU Weigen,TAO Benzao．Theory of Geopotential Model Determination[M]．Wuhan:Wuhan University Press,1990．(寧津生,邱衛(wèi)根,陶本藻．地球重力場模型理論[M]．武漢:武漢大學(xué)出版社,1990．)

[8] LIU Xiaogang．Theory and Methods of the Eath’s Gravity Field Model Recovery from GOCE Data[D]．Zhengzhou: Information Engineering University,2010．(劉曉剛．GOCE衛(wèi)星測量恢復(fù)地球重力場模型的理論與方法[D]．鄭州:信息工程大學(xué),2010．)

[9] PAVLIS N K,The Block-diagonal Least-squares Approach [R]．Washington D C:NASA Goddard Space Flight Center,1998．

[10] LEMOINE F G．The Development of the Joint NASA GSFC and the National Imagery and Mapping Agency(NIMA)Geopotential Model EGM96[R]．Washington D C: NASA Goddard Space Flight Center,1998．

[11] PAVLIS N K．Modeling and Estimation of a Low Degree Geopotential Model from Terrestrial Gravity Data[R]．Columbus:Ohio State University,1988．

[12] RAPP R H,PAVLIS N K．The Development and Analysis of Geopotential Coefficient Models to Spherical Harmonic Degree 360[J]．Journal of Geophysical Research,1990, 95(B13):21885-21911．

[13] COLOMBO O L．Numberical Methods for Harmonic Analysis on the Sphere[R]．Columbus:Ohio State University,1981．

[14] LU Zhonglian．Theory and Method of Earth’s Gravity[M]．Beijing:The People’s Liberation Army Press,1996．(陸仲連．地球重力場理論與方法[M]．北京:解放軍出版社,1996．)

[15] HEISKANEN W,MORITZ H．Physical Geodesy[M]．Beijing:Surveying and Mapping Press,1979．(海斯卡涅W A,莫里斯H．物理大地測量學(xué)[M]．北京:測繪出版社,1979．)

[16] HOFMANN-WELLENHOF B,MORITZ H．Physical Geodesy[M]．2nd edition．Berlin:Springer,2005．

[17] LI Jiancheng,CHEN Junyong,NING Jinsheng,et al．Approximation Theory of the Earth Gravity Field and Determination of the Chinese Gravity Geoid Model 2000 [M]．Wuhan:Wuhan University Press,2006．(李建成,陳俊勇,寧津生,等．地球重力場逼近理論與中國2000似大地水準面的確定[M]．武漢:武漢大學(xué)出版社,2006．)

[18] PVLIS N K,SALEH J．Error Propagation with Geographic Specificity for Very High Degree Geopotential Models in Gravity[C]∥IAG Symposia of Geoid and Space Missions．Berlin:Springer,2004．

[19] JEKELI C．The Exact Transformation between Ellipsoidal and Spherical Harmonic Expansions[J]．Manuscripta Geodaetica,1987,13:106-113．

[20] GLEASON D M．Comparing Ellipsoidal Corrections to the Transformation between the Geopotential’s Spherical and Ellipsoidal Spectrums[J]．Manuscripta Geodaetica,1988, 13:114-129．

(責(zé)任編輯:陳品馨)

The Application of Block-diagonal Least-squares Methods in Geopotential Model Determination

LI Xinxing1,WU Xiaoping1,LI Shanshan1,LIU Xiaogang2,CUI Zhiwei3
1．Institute of Geospatial Information,Information Engineering University,Zhengzhou 450001,China;2．Xi’an Research Institute of Surveying and Mapping,Xi’an 710054,China;3．61206 Troops,Beijing 100042,China

Six kinds of permutations of spherical harmonic coefficients and their block diagonal form of corresponding normal matrixes were studied as well as three kinds of approximation BDs:BD-1,BD-2,BD-3．The application of block-diagonal least-squares methods in geopotential model determination was analyzed and its advantage to quadrature formulas was showed by a simulated calculation．It illustrates that,with the global terrestrial gravity data and the satellite gravity field model whose accuracy is being improved continuously,block-diagonal least-squares methods are playing an increasingly important role in the calculation of the ultra-high-degree gravity field model．

block-diagonal;least-squares;ultra-high-degree geopotential model;EGM2008;orthogonality

LI Xinxing(1988—),male,master,majors in physical geodesy．

P223

A

1001-1595(2014)08-0778-08

國家自然科學(xué)基金(41274029;41304022);國家高技術(shù)研究發(fā)展專項(2013AA122502)

2013-12-06

李新星(1988—),男,碩士,主要從事物理大地測量研究。

E-mail:minibad＠126．com

LI Xinxing,WU Xiaoping,LI Shanshan,et al．The Application of Block-diagonal Least-squares Methods in Geopotential Model Determination[J]．Acta Geodaetica et Cartographica Sinica,2014,43(8):778-785．(李新星,吳曉平,李姍姍,等．塊對角最小二乘方法在確定全球重力場模型中的應(yīng)用[J]．測繪學(xué)報,2014,43(8):778-785．)

10．13485/j．cnki．11-2089．2014．0110

修回日期:2014-04-21