田家磊,李新星,吳曉平,邢志斌
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ǔ)、高效率的解算。
重力異常與位系數(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ǔ)難題。
由上述分析可知,利用分塊矩陣求解法矩陣,可簡(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)單高效。
對(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]
通過(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í)間。
采用上述改進(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í)。
為驗(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ù)求解,并且效率得到了大的提升。
為了說(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ì)。
數(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)方法最大的局限性。