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

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

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

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

1 塊對角最小二乘的基本原理

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

(1)

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

(2)

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

(3)

(4)

當重力數(shù)據(jù)分布于整個旋轉橢球面、數(shù)據(jù)在經(jīng)度方向分辨率一致、數(shù)據(jù)的權與經(jīng)度無關且關于赤道對稱時，權矩陣P是個對角陣，因此法矩陣N=ATPA可表示為

(5)

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

(6)

(7)

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

2 塊對角化最小二乘方法的改進

2.1 法矩陣約化求解

由上述分析可知，利用分塊矩陣求解法矩陣，可簡化為對一個塊的求解

(8)

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

(9)

(10)

由上述公式可看出，矩陣At具有極強的規(guī)律性，為簡化矩陣計算，設

(11)

(12)

則可得

(13)

則法矩陣Nt的計算如下

(14)

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

(15)

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

(16)

則

(17)

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

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

Tab.1Astatisticaltableoftheimprovementofthestorageobtainedbythemethodofblockdiagonalandmatrixreduction

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

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

2.2 Legendre函數(shù)計算存儲策略

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

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

2.3 Legendre函數(shù)對稱性

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

(18)

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

3 試驗及分析

3.1 精度改進

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

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

(19)

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

3.2 效率改進

為驗證本文所提改進方法在構建全球重力場模型的優(yōu)勢，使用3.2 GHz主頻、2 GB內存PC計算機，基于vs2008試驗平臺，進行模擬計算。

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

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

表2最小二乘方法恢復重力場模型時間

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

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

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

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

(20)

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

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

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

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

4 結 論

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

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

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

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

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

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