田家磊，李新星，吳曉平，邢志斌

1. 信息工程大學(xué)地理空間信息學(xué)院，河南 鄭州 450001； 2. 63850部隊(duì)，吉林 白城 137000； 3. 武漢大學(xué)測(cè)繪學(xué)院，湖北 武漢 430079

地球重力場(chǎng)模型是對(duì)地球重力場(chǎng)地面和空間觀測(cè)數(shù)據(jù)實(shí)施數(shù)學(xué)上的解析逼近的結(jié)果，在大地測(cè)量學(xué)、地球物理學(xué)、地震學(xué)、地質(zhì)學(xué)、海洋科學(xué)、國(guó)防科學(xué)、空間科學(xué)中的應(yīng)用非常廣泛[1-5]。利用重力場(chǎng)模型和相應(yīng)的球諧展開式可以計(jì)算多種擾動(dòng)重力場(chǎng)元(重力異常、擾動(dòng)重力、垂線偏差、大地水準(zhǔn)面高等)[6-8]。為了獲得高精度、高分辨率的地球重力場(chǎng)模型，國(guó)內(nèi)外眾多機(jī)構(gòu)投入了巨大的人力、物力與財(cái)力[9-10]，目前地球重力場(chǎng)模型已經(jīng)發(fā)展到2159完全階次[11]。2008年美國(guó)國(guó)家地理空間情報(bào)局(US National Geospatial-Intelligence Agency)發(fā)布了目前應(yīng)用范圍最廣、高精度超高階地球重力場(chǎng)模型EGM2008[12-13]。2014年德國(guó)地學(xué)研究中心(Geo Forschungs Zentrum)和法國(guó)空間大地測(cè)量組(Groupe Recherches Geodesie Spatiale)聯(lián)合發(fā)布了EIGEN(European Improved Gravity Model of the Earth by New Techniques)系列模型中的最新模型EIGEN-6C4[14-16]。這兩個(gè)模型是目前普遍采用的超高階重力場(chǎng)模型。

衛(wèi)星重力數(shù)據(jù)用來(lái)解算地球重力場(chǎng)模型的中低階部分，地面與海洋重力數(shù)據(jù)用來(lái)解算地球重力場(chǎng)模型的高階部分[17-19]。目前重力場(chǎng)模型位系數(shù)的求解方法主要有兩種：數(shù)值積分方法(numerical quadrature，NQ)又稱調(diào)和分析方法，還有最小二乘方法(least squares，LS)。聯(lián)合衛(wèi)星重力資料確定中低階位系數(shù)的方法也有多種，包括基于原始觀測(cè)數(shù)據(jù)的處理和基于衛(wèi)星重力場(chǎng)模型位系數(shù)的處理，其中利用最小二乘方法解算位系數(shù)的法矩陣疊加方法更具有優(yōu)勢(shì)，能夠更加詳細(xì)地顧及地面數(shù)據(jù)和衛(wèi)星觀測(cè)數(shù)據(jù)之間的權(quán)重關(guān)系[20]。采用最小二乘方法構(gòu)建超高階重力場(chǎng)模型過(guò)程中，觀測(cè)方程中的未知數(shù)(模型位系數(shù))和觀測(cè)值數(shù)量大，相應(yīng)的法方程規(guī)模巨大，大型法方程的構(gòu)建與直接求解難以實(shí)現(xiàn)[21-23]。目前國(guó)內(nèi)外利用最小二乘方法計(jì)算超高階重力場(chǎng)模型，均使用超算機(jī)或者采用集成計(jì)算機(jī)網(wǎng)絡(luò)，基于并行算法進(jìn)行的，計(jì)算成本高，耗時(shí)長(zhǎng)[24-25]?；谧钚《朔椒ㄔ谇蠼獬唠A模型問題上的優(yōu)勢(shì)和面臨的難題，本文探索該海量計(jì)算過(guò)程中難題的解決方案，力爭(zhēng)實(shí)現(xiàn)超高階重力場(chǎng)模型的小存儲(chǔ)、高效率的解算。

1 塊對(duì)角最小二乘的基本原理

重力異常與位系數(shù)之間關(guān)系的基本數(shù)學(xué)模型為[6]

(1)

最小二乘方法可得位模型系數(shù)的解及相應(yīng)的協(xié)方差陣[20]

(2)

根據(jù)式(1)，矩陣A的各個(gè)元素由式(3)和式(4)表示

(3)

(4)

當(dāng)重力數(shù)據(jù)分布于整個(gè)旋轉(zhuǎn)橢球面、數(shù)據(jù)在經(jīng)度方向分辨率一致、數(shù)據(jù)的權(quán)與經(jīng)度無(wú)關(guān)且關(guān)于赤道對(duì)稱時(shí)，權(quán)矩陣P是個(gè)對(duì)角陣，因此法矩陣N=ATPA可表示為

(5)

為方便研究，本文近似權(quán)矩陣P為單位陣，根據(jù)球諧函數(shù)的正交特性知，由以上觀測(cè)方程得到的法方程矩陣N是一稀疏矩陣，即

(6)

(7)

以2160階次重力場(chǎng)模型前2159階位系數(shù)的最小二乘方法求解為例，Nmax=2160，系數(shù)矩陣列分塊個(gè)數(shù)以及生成的法矩陣塊個(gè)數(shù)均為4×2159，法方程的求解可簡(jiǎn)化為對(duì)4×2159個(gè)法矩陣分別求解。這些法矩陣塊中最大的塊為int(Nmax/2)×int(Nmax/2)，即1080×1080，因此最大的系數(shù)列塊矩陣A1的大小為(2160×4320)×1080，存儲(chǔ)A1需要75.1 G的空間，因此對(duì)于普通性能的計(jì)算機(jī)，塊對(duì)角化后的超高階模型求解仍面臨存儲(chǔ)難題。

2 塊對(duì)角化最小二乘方法的改進(jìn)

2.1 法矩陣約化求解

由上述分析可知，利用分塊矩陣求解法矩陣，可簡(jiǎn)化為對(duì)一個(gè)塊的求解

(8)

利用式(8)求解法矩陣過(guò)程中，需要存儲(chǔ)系數(shù)列塊矩陣At，最大的At矩陣為A1，存儲(chǔ)該矩陣需要75.1 GB空間，較難實(shí)現(xiàn)，針對(duì)此問題，試分析矩陣At的特點(diǎn)，尋找簡(jiǎn)化方法。At可表示為式(9)、式(10)兩種形式之一

(9)

(10)

由上述公式可看出，矩陣At具有極強(qiáng)的規(guī)律性，為簡(jiǎn)化矩陣計(jì)算，設(shè)

(11)

(12)

則可得

(13)

則法矩陣Nt的計(jì)算如下

(14)

綜合式(11)、式(12)，可得Nt的表達(dá)式為

(15)

因?yàn)楦窬W(wǎng)重力異常分布是按照全球等角格網(wǎng)分布的，即經(jīng)度間隔是一致的，可得

(16)

則

(17)

對(duì)于求解2160階次的超高階重力場(chǎng)模型而言，使用5′×5′格網(wǎng)重力異常，共計(jì)2160個(gè)緯度值，因此所有Tm θi值的存儲(chǔ)量為2160×1080，即需要開辟17.8 MB的空間。綜上所述，經(jīng)過(guò)塊對(duì)角化以及本文提出的法矩陣快速計(jì)算方案，求解一法矩陣的過(guò)程，大大節(jié)約了存儲(chǔ)空間，提升效果見表1。

表1塊對(duì)角以及法矩陣約化求解方法對(duì)存儲(chǔ)量改進(jìn)的效果統(tǒng)計(jì)表

Tab.1Astatisticaltableoftheimprovementofthestorageobtainedbythemethodofblockdiagonalandmatrixreduction

矩陣類型觀測(cè)方程/TB對(duì)角化后最大“塊”/GB約化求解法存儲(chǔ)的Tm矩陣/MB系數(shù)矩陣316.875.117.8法矩陣158.48.98.9

由表1可知，經(jīng)過(guò)待求位系數(shù)的重新排序的塊對(duì)角化以及法矩陣約化求解改進(jìn)后，使得系數(shù)矩陣和法矩陣的必要存儲(chǔ)空間有了大幅度的縮減，使得計(jì)算變得簡(jiǎn)單高效。

2.2 Legendre函數(shù)計(jì)算存儲(chǔ)策略

對(duì)于2160階次的重力場(chǎng)模型求解，當(dāng)一次性存取所有的Legendre函數(shù)值，需要存儲(chǔ)變量個(gè)數(shù)為2160×(2159×2164)/2，對(duì)應(yīng)的存儲(chǔ)空間37.6 GB，這么大的數(shù)據(jù)量很難使用內(nèi)存存儲(chǔ)，如果寫入文件進(jìn)行讀取的話會(huì)大大降低計(jì)算效率。如果每求一個(gè)At矩陣便遍歷求解一次legendre值，那么需要循環(huán)計(jì)算次數(shù)為2160×4，這也將大大降低計(jì)算效率。

圖1 Legendre遞推示意圖[3]Fig.1 Legendre recursive sketch map[3]

2.3 Legendre函數(shù)對(duì)稱性

通過(guò)對(duì)系數(shù)矩陣At進(jìn)行分析可知，當(dāng)n-m為偶數(shù)時(shí)，At矩陣上半部分與下半部分是對(duì)稱的，當(dāng)n-m為奇數(shù)時(shí)，At矩陣上半部分與下半部分?jǐn)?shù)值相等，但符號(hào)相反。究其原因，是由于全球格網(wǎng)數(shù)據(jù)分布是關(guān)于赤道對(duì)稱的，滿足θi=π-θN-1-i，且締合Legendre函數(shù)滿足以下特性

(18)

因此在計(jì)算矩陣At及法矩陣時(shí)，只需要計(jì)算上半部分即可，最終求得的法矩陣N乘以2即為最終的法矩陣，該方法能夠節(jié)省一半的存儲(chǔ)空間及一半的計(jì)算時(shí)間。

3 試驗(yàn)及分析

3.1 精度改進(jìn)

采用上述改進(jìn)的計(jì)算方法和存儲(chǔ)方式，利用EGM2008模型2159階位系數(shù)模擬計(jì)算全球5′×5′格網(wǎng)重力異常點(diǎn)值數(shù)據(jù)，實(shí)現(xiàn)2159階次超高階重力場(chǎng)模型的最小二乘方法實(shí)現(xiàn)，并與EGM2008模型前2159階位系數(shù)進(jìn)行比較，作為模型恢復(fù)方法精度的衡量標(biāo)準(zhǔn)，并與輪胎調(diào)和分析方法恢復(fù)重力場(chǎng)模型的精度進(jìn)行對(duì)比分析[26]，結(jié)果見圖3。圖3分別統(tǒng)計(jì)了兩種方法恢復(fù)的模型系數(shù)與EGM2008系數(shù)差的絕對(duì)值的自然對(duì)數(shù)。

使用EGM2008作為標(biāo)準(zhǔn)模型，計(jì)算兩種方法恢復(fù)的2159階模型系數(shù)的誤差階RMS并進(jìn)行比較，比較結(jié)果見圖4，其計(jì)算公式為[3]

(19)

由圖4可看出，塊對(duì)角最小二乘方法恢復(fù)模型系數(shù)的誤差階RMS明顯小于使用數(shù)值積分方法恢復(fù)的模型系數(shù)，最小二乘方法恢復(fù)的模型整體精度比數(shù)值積分方法高出至少5個(gè)數(shù)量級(jí)。

3.2 效率改進(jìn)

為驗(yàn)證本文所提改進(jìn)方法在構(gòu)建全球重力場(chǎng)模型的優(yōu)勢(shì)，使用3.2 GHz主頻、2 GB內(nèi)存PC計(jì)算機(jī)，基于vs2008試驗(yàn)平臺(tái)，進(jìn)行模擬計(jì)算。

試驗(yàn)統(tǒng)計(jì)了分別使用3°×3°、1°×1°、30′×30′及5′×5′的格網(wǎng)重力異常數(shù)據(jù)，利用傳統(tǒng)塊對(duì)角最小二乘法恢復(fù)59階、179階、359階和2159階次模型的計(jì)算時(shí)間，并利用改進(jìn)后的最小二乘方法進(jìn)行同樣的計(jì)算，統(tǒng)計(jì)結(jié)果見表2。

如表2所示，采用塊對(duì)角最小二乘方法計(jì)算模型位系數(shù)，效率相對(duì)低下，計(jì)算359階位系數(shù)需要近2 h，對(duì)于2159階重力場(chǎng)模型，由于其系數(shù)矩陣At存儲(chǔ)量龐大，無(wú)法實(shí)現(xiàn)系數(shù)的求解。相比傳統(tǒng)塊對(duì)角最小二乘法，采用改進(jìn)后的方法計(jì)算效率有了明顯提升，與恢復(fù)359階模型使用2 h相比，改進(jìn)方法只需要20 s。更重要的是，改進(jìn)后的方法實(shí)現(xiàn)了2159超高階重力場(chǎng)模型的求解，耗時(shí)不足2 h。

表2最小二乘方法恢復(fù)重力場(chǎng)模型時(shí)間

Tab.2 Time of recovering geopotential modelby least-squares method ms

綜上，改進(jìn)后的最小二乘方法不僅實(shí)現(xiàn)了普通計(jì)算機(jī)上超高階模型的位系數(shù)求解，并且效率得到了大的提升。

3.3 數(shù)據(jù)的精度估計(jì)

為了說(shuō)明相比數(shù)值積分方法，最小二乘方法在精度評(píng)估方面的優(yōu)勢(shì)，本文對(duì)EGM2008模型前359階位系數(shù)模擬的30′×30′格網(wǎng)重力異常加入白噪聲，利用最小二乘方法恢復(fù)模型位系數(shù)過(guò)程中，記錄VTPV的估計(jì)值，利用式(20)可獲取觀測(cè)值的單位權(quán)中誤差

(20)

式中，n為觀測(cè)量個(gè)數(shù)；t為未知數(shù)個(gè)數(shù)。

通過(guò)對(duì)上述30′×30′重力異常數(shù)據(jù)分別加入5×10-5m/s2和10×10-5m/s2白噪聲后，進(jìn)行模型位系數(shù)的恢復(fù)，最終得到的統(tǒng)計(jì)結(jié)果見表3。

表3 精度估計(jì)情況統(tǒng)計(jì)

根據(jù)表3結(jié)果可知，最小二乘方法能夠在一定程度上反映出數(shù)據(jù)與模型的吻合程度，即一定程度上評(píng)估數(shù)據(jù)的精度情況，這是最小二乘相對(duì)于數(shù)值積分方法不可比擬的優(yōu)勢(shì)。

4 結(jié) 論

數(shù)值積分方法實(shí)質(zhì)上就是法矩陣為對(duì)角陣的最小二乘方法，即認(rèn)為模型系數(shù)之間均是不相關(guān)的，協(xié)方差均為0，而塊對(duì)角最小二乘則估計(jì)了同次不同階系數(shù)之間的相關(guān)性，其更符合實(shí)際情況，也能夠恢復(fù)更高精度的重力場(chǎng)模型。數(shù)值積分方法不能對(duì)重力異常數(shù)據(jù)進(jìn)行精度評(píng)估，塊對(duì)角最小二乘方法可以。相比數(shù)值積分方法，最小二乘方法無(wú)法求解Nmax階系數(shù)，例如5′×5′分辨率的格網(wǎng)重力異常只能求解2159階次的系數(shù)，第2160階的系數(shù)必需使用數(shù)值積分方法才能求得。另外，數(shù)值積分方法使用的重力異常數(shù)據(jù)通常為點(diǎn)值，即包含全頻段信號(hào)，所反演的系數(shù)才是實(shí)際模型位系數(shù)，而最小二乘方法所使用的格網(wǎng)重力異常數(shù)據(jù)應(yīng)只包含所求階次的能量，否則其信號(hào)在位系數(shù)之間會(huì)發(fā)生泄漏。

圖2 標(biāo)準(zhǔn)向前列推法及跨階次發(fā)方法的溢出情況比較[3]Fig.2 Comparison of overflow between standard forward method and cross order method

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

圖4 誤差階RMS比較Fig.4 Comparison of error degree RMS

本文基于Pavlis提出的塊對(duì)角最小二乘方法，通過(guò)對(duì)法矩陣求解方程進(jìn)行約化、對(duì)Legendre函數(shù)的計(jì)算和存儲(chǔ)方式進(jìn)行設(shè)計(jì)，結(jié)合締合Legendre函數(shù)關(guān)于赤道的對(duì)稱性，解決了大型矩陣存儲(chǔ)及計(jì)算效率低下的難題，實(shí)現(xiàn)了超高階重力場(chǎng)模型最小二乘方法的小存儲(chǔ)、高效率的解算。結(jié)果顯示，精度相比數(shù)值積分方法提高近5個(gè)量級(jí)，計(jì)算效率提高近300倍，并精確估計(jì)了重力異常數(shù)據(jù)的精度。

本文最小二乘改進(jìn)方法的局限性：采用上述高效的超高階模型位系數(shù)最小二乘方法，在式(8)中，有權(quán)矩陣P為單位權(quán)陣這一假設(shè)。對(duì)其進(jìn)行拓展，當(dāng)位于同一經(jīng)度的格網(wǎng)重力異常精度相同，或位于同一緯度的格網(wǎng)重力異常精度相同，這兩種情況下也可以使用上述改進(jìn)的快速方法。當(dāng)出現(xiàn)每一片區(qū)域都有自己不同的精度這一情況下，上述方法是不適用的，這也是該改進(jìn)方法最大的局限性。