球諧旋轉(zhuǎn)變換結(jié)合非全次Legendre方法的局部六邊形網(wǎng)格重力場(chǎng)球諧綜合

李新星, 李建成, 劉曉剛, 范昊鵬, 靳超

1 武漢大學(xué)測(cè)繪學(xué)院， 武漢 430079 2 信息工程大學(xué)地理空間信息學(xué)院， 鄭州 450001 3 西安測(cè)繪研究所， 西安 710054 4 61365部隊(duì)， 天津 300000

0 引言

傳統(tǒng)基于等角地理網(wǎng)格的地球重力場(chǎng)描述面臨著許多無法克服的局限性，包括網(wǎng)格面積不等、重力場(chǎng)平滑因子復(fù)雜、觀測(cè)數(shù)據(jù)冗余、積分離散化誤差大等，而六邊形的“離散球面網(wǎng)格系統(tǒng)”(Discrete Global Grid System,DGGS)因其良好的特性，在遙感、氣象、水文等領(lǐng)域取得了廣泛應(yīng)用(Sahr et al.,2003). 作為地理網(wǎng)格的延伸，美國著名學(xué)者富勒(Fuller)最早進(jìn)行了球面離散化的相關(guān)研究并且發(fā)明了“網(wǎng)格球頂”(Geodesic Dome)結(jié)構(gòu)，之后許多學(xué)者關(guān)于球面網(wǎng)格系統(tǒng)的研究都直接或間接地受到了他的影響(張永生等,2007). 近年來，隨著DGGS網(wǎng)格系統(tǒng)不斷完善，其已經(jīng)在諸如矢量數(shù)據(jù)結(jié)構(gòu)(Dutton，1999)、空間數(shù)據(jù)索引(Dutton,1996a)、制圖綜合(Dutton,1996b)、全球動(dòng)態(tài)數(shù)據(jù)結(jié)構(gòu)(Goodchild and Shiren,1992)、全球土地監(jiān)測(cè)(Suess et al.,2004)、全球海洋分析(Kidd et al.,2003)、全球氣候模型(Randall et al.,2002)、海洋路徑規(guī)劃、城市交通規(guī)劃(Brodsky,2018)等方面得到了應(yīng)用.

基于六邊形的DGGS網(wǎng)格本應(yīng)在重、磁場(chǎng)這一類具有各向同性特征的物理場(chǎng)中發(fā)揮更大作用，但是目前相關(guān)研究很少. 類六邊形網(wǎng)格結(jié)構(gòu)在重磁場(chǎng)應(yīng)用方面，Vestine等(1963a,1963b)較早研究了基于球面均勻分布點(diǎn)的地磁場(chǎng)的Poisson積分應(yīng)用，Eicker(2008)將多種點(diǎn)位均勻分布方式引入到徑向基函數(shù)中進(jìn)而求解衛(wèi)星重力場(chǎng)取得了較好結(jié)果，Slobbe等(2012)將Fibonacci均勻分布引入球形斯列普基函數(shù)求解平均海水面和海面地形，Ran等(2018)將Fibonacci格網(wǎng)用于Mascon方法來反演冰蓋質(zhì)量變化，以上均是在數(shù)值積分求解方面引入均勻分布結(jié)構(gòu)，而鮮有在球諧分析與綜合中開展研究，其中最主要的原因在于，基于六邊形網(wǎng)格系統(tǒng)的點(diǎn)位并不像地理網(wǎng)格那樣是等緯度分布的，如圖1所示，從而需要大量締合Legendre函數(shù)(associated Legendre functions,ALFs)的求解，這對(duì)計(jì)算能力提出了更高要求(Pavlis,1988; 李新星等,2014; 田家磊等,2018).

圖1 基于全球六邊形網(wǎng)格的360階全球重力異常分布圖Fig.1 360-degree gravity anomaly values based on global hexagonal grids

同時(shí),觀測(cè)資料的精度與數(shù)量的增多使得對(duì)應(yīng)球諧展開的階次也不斷升高,對(duì)計(jì)算能力同樣提出了挑戰(zhàn). 例如目前有720階次的NGDC-720地殼磁異常模型(Maus,2010),有經(jīng)典的2160階次地球重力場(chǎng)模型EGM2008(Pavlis et al.,2012),還有10800階次的地形、布格異常、均衡異常模型(Balmino et al.,2012; Hirt and Rexer,2015). 受限于球諧展開所需龐大的計(jì)算量,現(xiàn)有網(wǎng)格數(shù)據(jù)的分辨率已經(jīng)遠(yuǎn)遠(yuǎn)走在現(xiàn)有模型階次對(duì)應(yīng)分辨率的前列. 例如地磁場(chǎng)的磁異常網(wǎng)格數(shù)據(jù)EMAG2v3分辨率達(dá)到2′,根據(jù)Nyquest采樣定律,因磁場(chǎng)為偶極子場(chǎng),其對(duì)應(yīng)10800階次的球諧展開; sandwell和DTU15版本的海洋衛(wèi)星測(cè)高重力異常分辨率達(dá)到1′(Sandwell and Smith,2008; Andersen and Knudsen,2016), 因重力場(chǎng)為標(biāo)量場(chǎng),所以其對(duì)應(yīng)10800階次; GEBCO2014和GEBCO2019海底地形網(wǎng)格數(shù)據(jù)的分辨率分別達(dá)到了30″和15″,對(duì)應(yīng)21600和43200階次; SRTM3、ARTER和ALOS全球DEM數(shù)據(jù)分辨率分別達(dá)到90 m、30 m和12 m,對(duì)應(yīng)的球諧展開階次更高,因此,球諧展開中的嚴(yán)重耗時(shí)問題也是急需解決的.

針對(duì)ALFs及其相應(yīng)積分、導(dǎo)數(shù)等函數(shù)的計(jì)算問題,許多學(xué)者從遞推方法、計(jì)算精度、計(jì)算效率和溢出問題等方面做了大量的工作. 對(duì)于ALFs的遞推方法,主要包括行推法、列推法、Belikov和跨階次遞推等四種方法,其中Belikov方法的計(jì)算效率最佳(李新星,2013). 為了解決遞推過程中的溢出問題,Fukushima (2012)提出了X數(shù)法實(shí)現(xiàn)了任意階次ALFs的遞推計(jì)算,代價(jià)是損失了一定的計(jì)算效率; Xing等(2020)改進(jìn)了Belikov方法,解決了Belikov原遞推公式在高階次出現(xiàn)溢出的情況,例如緯度為75°的ALFs當(dāng)階次達(dá)到4800時(shí)，ALFs的計(jì)算出現(xiàn)溢出; 張傳定等(2004)首次提出ALFs的跨階次遞推方法,并經(jīng)驗(yàn)證在20000階次以內(nèi)其計(jì)算準(zhǔn)確且不會(huì)發(fā)生溢出(于錦海等,2015).

關(guān)于ALFs遞推計(jì)算效率的進(jìn)一步提升,目前還沒有較為顯著的方法. 學(xué)者大多繞過ALFs的計(jì)算,通過盡量減小ALFs的計(jì)算量來尋求提升球諧綜合(Spherical Harmonic Synthesis,SHS)和球諧分析(Spherical Harmonic Analysis,SHA)計(jì)算效率的方法,例如在SHS和SHA中要求相應(yīng)的計(jì)算點(diǎn)等緯度網(wǎng)格分布,從而計(jì)算同一緯圈上所有點(diǎn)的模型值只需要計(jì)算一次該緯度的ALFs. Clenshaw(1957)提出Clenshaw求和方法避免大部分Legendre函數(shù)值的求解來提高SHS的計(jì)算效率,同時(shí)快速傅里葉變換(Fast Fourier Transform,FFT)技術(shù)的使用大大加速了球諧分析與綜合的計(jì)算,例如Colombo(1981)在球諧分析中使用了一維FFT,Sneeuw和Bun (1996)利用球面到輪胎面映射實(shí)現(xiàn)了二維FFT的球諧分析與綜合. 盡管FFT技術(shù)加快了SHS的計(jì)算速度,但其僅限于在計(jì)算規(guī)則網(wǎng)格分布下具有相同高度的點(diǎn)集,針對(duì)非等緯度、不同高度分布點(diǎn)的SHS和SHA,有學(xué)者采用了泰勒級(jí)數(shù)展開或等間隔內(nèi)插的思想來提升效率(Kunis and Potts,2003; Moazezi et al.,2016; Seif et al.,2018).

實(shí)際應(yīng)用中,更多是針對(duì)局部重力場(chǎng)或磁場(chǎng)數(shù)據(jù)的處理,為了解決因球諧方法計(jì)算量龐大而不便于局部使用的問題,學(xué)者們避開球諧模型,使用局部擬合、等效源、矩諧分析、球冠諧分析等方法來進(jìn)行局部物理場(chǎng)的建模(徐文耀和朱崗昆,1984; 高金田等,2005; Younis,2013),而上述方法存在的邊界效應(yīng)、基函數(shù)不正交、非整階Legendre遞推等問題還需要深入研究,如果球諧模型變得更加高效,也必然能夠在局部重磁場(chǎng)模型構(gòu)建中發(fā)揮重要作用.

本文從ALFs求解方面入手,尋求提升其計(jì)算效率的方法. Jekeli等(2007)總結(jié)了不同緯度下ALFs函數(shù)值的量級(jí)隨階次升高的衰減關(guān)系,通過利用圖2a所示的逆向行推方法,求解ALFs值的有效部分,對(duì)于絕對(duì)值較小(例如<10-15)的值默認(rèn)為0而不進(jìn)行計(jì)算,從而節(jié)約了遞推過程中的計(jì)算量,在全球SHS過程中,該方法理論上能夠節(jié)約36%的計(jì)算量,但是由于其逆向行推方法計(jì)算效率偏低,并沒有給出最終的計(jì)算耗時(shí)比較.

圖2 ALFs的逆向行推法(a)和跨階次遞推(b)方法Fig.2 The increasing order (a) and cross-degree-order (b) recursion formula for ALFs

本文利用該思想,結(jié)合跨階次遞推方法,提出了非全次Legendre(Non-full-order Legendre,NFOL)計(jì)算方法,實(shí)現(xiàn)了高緯度區(qū)域任意點(diǎn)的快速球諧綜合. NFOL方法是在利用跨階次遞推方法計(jì)算ALFs時(shí),當(dāng)遞推到達(dá)ALFs有效值(例如>10-15)邊界即可停止遞推,從而減小計(jì)算量,提升效率,由于有效值分布在全部“階”中“次”偏小的區(qū)域,所以該遞推只需要遞推全部階的部分“次”.

對(duì)于中低緯度的ALFs計(jì)算,由于ALFs有效值區(qū)域幾乎遍布整個(gè)階次,所以NFOL方法在低緯度區(qū)域并不適用,針對(duì)該問題,本文引入球諧旋轉(zhuǎn)變量變換(Spherical Harmonic Rotation,SHR)方法,該方法實(shí)現(xiàn)了坐標(biāo)系旋轉(zhuǎn)下,對(duì)球諧展開系數(shù)的變換. 利用該方法將低緯度待求解區(qū)域的中心通過坐標(biāo)系旋轉(zhuǎn)移至北極點(diǎn),從而使得旋轉(zhuǎn)后的點(diǎn)位均位于“高緯度”,使得與之對(duì)應(yīng)的ALFs有效值區(qū)域集中于低“次”部分,從而可采用NFOL方法提升SHS效率.

SHR方法在多領(lǐng)域具有重要作用. 地磁場(chǎng)模型中,該方法對(duì)于不同坐標(biāo)系下磁場(chǎng)模型球諧系數(shù)間的比較、日變磁場(chǎng)Sq的內(nèi)/外源場(chǎng)分離等具有重要作用(De Santis et al.,1996). 該方法求解過程中,Wigner D矩陣起到了重要的作用,許多學(xué)者研究了D矩陣及其遞推求解方法(Aubert,2013; Fukushima,2017).國內(nèi),許厚澤和蔣福珍(1964)較早研究了SHR原理并推導(dǎo)了公式,其在飛行器軌道擾動(dòng)引力快速賦值中得到應(yīng)用(任萱,1985; 鄭偉等,2011; 王建強(qiáng)等,2013),但是這些應(yīng)用中采用的球諧系數(shù)階次較低,且其變換精度均沒有經(jīng)過嚴(yán)格的精度評(píng)估. Fukushima (2017)利用X數(shù)法實(shí)現(xiàn)了2160階次球諧系數(shù)的90°旋轉(zhuǎn).

本研究基于NFOL和SHR方法,分別實(shí)現(xiàn)了高緯度南極地區(qū)低分辨率的六邊形網(wǎng)格點(diǎn)重力異常、低緯度加里曼丹島地區(qū)的高分辨率六邊形網(wǎng)格點(diǎn)重力異常的構(gòu)建,通過比較分析,驗(yàn)證了方法的精度與效率,得到了有益結(jié)果,成果值得推廣.

1 非全次Legendre函數(shù)的遞推

1.1 球諧展開基本公式

對(duì)于球面上的平方可積函數(shù),通過SHA能夠得到對(duì)應(yīng)的球諧展開系數(shù),這里以重力場(chǎng)中重力異常的球諧展開為例,其與重力場(chǎng)位系數(shù)之間的關(guān)系如下：

(1)

(2)

(3)

(4)

n≥2,n≠m

(5)

其中

(6)

遞推初始值為

(7)

當(dāng)m≥2時(shí),采用以下遞推公式(于錦海等,2015)：

(8)

其中

(9)

(10)

(11)

當(dāng)m>2時(shí),遞推系數(shù)cn m,dn m,hn m的值均小于1,因此跨階次遞推是穩(wěn)定可靠的,遞推過程如圖2b所示.

1.2 非全次Legendre方法

圖3 m=100時(shí),θ=1°、45°和89°的ALFs的大小Fig.3 The value of ALFs with θ=1°,45° and 89° when m=100

m≈nsinθ.

(12)

圖4 ALFs值中穩(wěn)定振蕩區(qū)與快速衰減區(qū)的分界線,圖中綠色區(qū)域值較小,近似為0Fig.4 The boundary between the stable oscillation zone and the fast decay zone of the value of ALFs,the green area in the figure illustrates that the value is extremely small,approximately 0

由于ALFs的值在該分界線周邊有一個(gè)衰減過程,如圖5所示,因此想要保證整個(gè)ALFs的計(jì)算精度,不能直接以式(12)所給出的直線作為分界線,而應(yīng)該適當(dāng)?shù)叵騧增大和n減小的方向移動(dòng)該直線,得到新的分界線

m≈nsinθ+t,

(13)

圖5 分界線周邊ALFs值的衰減過程,圖中藍(lán)色區(qū)域值小于10-20,可近似為0Fig.5 The decay process of the values of ALFs near the boundary,the blue area in the figure represents that the values are less than 10-20,approximately 0

其中參數(shù)t的大小決定了ALFs遞推值的精度和計(jì)算效率,t越大,精度越有保證,但ALFs需要遞推的區(qū)域會(huì)變大(見圖6),相應(yīng)計(jì)算效率會(huì)降低,關(guān)于參數(shù)t的確定,會(huì)在后續(xù)實(shí)驗(yàn)部分給出.

圖6 參數(shù)t的含義Fig.6 The meaning of parameter t

因此,上述方法確定了ALFs的遞推范圍,我們開創(chuàng)性地避開Jekeli等(2007)所采用的逆向行推方法,采用跨階次遞推方法計(jì)算了非全次遞推范圍內(nèi)的ALFs值,解決了逆向行推不穩(wěn)定且計(jì)算效率低下的問題.由圖5可知,對(duì)于非高緯度位置處的ALFs,顯然上述方法并不能有效提升計(jì)算效率,因此引入球諧系數(shù)的旋轉(zhuǎn)變量變換,通過換極將計(jì)算區(qū)域移至“高緯”區(qū)域,從而提升計(jì)算效率.

2 球諧旋轉(zhuǎn)變量變換(SHR)

2.1 沿z軸的旋轉(zhuǎn)變換

已知球面上一點(diǎn)Q(r,θ,λ),當(dāng)坐標(biāo)系O-xyz沿著z軸旋轉(zhuǎn)α角,記為Rz(α),得到新的坐標(biāo)系O-x′y′z,則Q點(diǎn)在新坐標(biāo)系下的球坐標(biāo)為(r,θ,λ1),其中λ1=λ-α,即λ=λ1+α，如圖7所示.

圖7 沿z軸的坐標(biāo)系旋轉(zhuǎn)Fig.7 Rotation of the coordinate system along axis z

因此Q點(diǎn)的重力異常在新的球坐標(biāo)系O-x′y′z下可展開為新的球諧系數(shù)

(14)

根據(jù)λ與λ1的以下三角關(guān)系式：

cosmλ=cosm(λ1+α)=cosmλ1cosmα

-sinmλ1sinmα，

sinmλ=sinm(λ1+α)=sinmλ1cosmα

+cosmλ1sinmα，

(15)

代入式(1)并與式(14)進(jìn)行比較,可得兩個(gè)坐標(biāo)系下球諧系數(shù)之間的關(guān)系

(16)

式(16)即為坐標(biāo)系沿z軸旋轉(zhuǎn)α角時(shí),對(duì)應(yīng)的球諧系數(shù)的變換公式.

2.2 沿y軸的旋轉(zhuǎn)變換

已知球面上一點(diǎn)Q(r,θ,λ),當(dāng)坐標(biāo)系O-xyz沿著y軸旋轉(zhuǎn)β角,記為Ry(β),得到新的坐標(biāo)系O-x″yz″.

如圖8,Q點(diǎn)在新坐標(biāo)系下的球坐標(biāo)為(r,Θ,Λ),其中Θ為新坐標(biāo)系下的球心余緯,Λ為新坐標(biāo)系下的球心經(jīng)度,Λ1與Λ互補(bǔ),Q點(diǎn)的重力異常在新的球坐標(biāo)系O-x″yz″下可展開為新的球諧系數(shù),表達(dá)式與式(1)類似，

圖8 沿y軸的坐標(biāo)系旋轉(zhuǎn)Fig.8 Rotation of the coordinate system along axis y

(17)

(18)

(19)

(20)

(21)

(22)

其中

(23)

(24)

(25)

(26)

(26)式可利用D函數(shù)的對(duì)稱性及其與該兩組變量之間的對(duì)應(yīng)關(guān)系得證,具體見Aubert(2013),即圖9d中k=4與圖9c中m=4對(duì)應(yīng)的元素值相等.

圖10 β=90°的30階值的切片圖(a) 30階球諧轉(zhuǎn)換矩陣的值; (b) 與比較說明式(29)對(duì)稱結(jié)構(gòu).Fig.10 Slice graph of the value of with β=90° and the n up to 30(a) Values of transformation matrix of in Eq.(29).

2.3 沿y軸90°的旋轉(zhuǎn)變換

(27)

其中

(28)

(29)

(30)

(31)

(32)

其中

(33)

(34)

(35)

上述遞推方法的計(jì)算過程中,Hn,m隨著n的增大快速減小,當(dāng)其絕對(duì)值小于2.225×10-308的時(shí)候,數(shù)值發(fā)生溢出而強(qiáng)制為0,會(huì)影響后續(xù)遞推工作,因此在該遞推過程中應(yīng)使用X-數(shù)法進(jìn)行求解.

2.4 將(θ,λ)視為北極點(diǎn)的旋轉(zhuǎn)變換

Zotter (2009)文中指出,Ry(θ)可分解為以下旋轉(zhuǎn)的組合：

Ry(θ)=Rz(π/2)Ry(π/2)Rz(θ)Ry(-π/2)Rz(-π/2)

(36)

還可以寫為

Ry(θ)=Rz(π/2)Ry(π/2)Rz(θ-π)Ry(π/2)Rz(π/2)

(37)

因此,只需要計(jì)算沿y軸旋轉(zhuǎn)90°的球諧系數(shù)變換以及簡(jiǎn)單的z軸變換,就可以實(shí)現(xiàn)任意角度的旋轉(zhuǎn)變換,即

Ry(θ)Rz(λ)=Rz(π/2)Ry(π/2)Rz(θ-π)Ry(π/2)

×Rz(π/2)Rz(λ)，

(38)

其中系數(shù)的變換采用公式(16)和(18).

3 數(shù)值實(shí)驗(yàn)及分析

為了驗(yàn)證NFOL和SHR方法在局部六邊形網(wǎng)格重力場(chǎng)快速球諧綜合方面的應(yīng)用,共進(jìn)行三個(gè)實(shí)驗(yàn),計(jì)算環(huán)境如下：2.30GHz主頻的Intel Core i5-6200U CPU和8 GB主內(nèi)存的PC機(jī),64位Windows XP OS操作系統(tǒng)下VS2017編譯環(huán)境.

3.1 參數(shù)t的確定

第1節(jié)中指出,使用NFOL方法進(jìn)行球諧綜合,需要確定公式(13)中的參數(shù)t,t取值越大,計(jì)算精度越高而對(duì)應(yīng)效率越低,反之亦然. 因此,需要在滿足精度要求的情況下,合理的選擇t的大小來提高球諧綜合的計(jì)算效率.

利用2160階2159次的EGM2008模型位系數(shù),計(jì)算地球表面點(diǎn)的模型重力異常值,利用不同緯度、不同t取值的NFOL方法計(jì)算重力異常,并與常規(guī)全階次方法計(jì)算的結(jié)果作為真值進(jìn)行比較,結(jié)果統(tǒng)計(jì)見表1.

表1 不同緯度、不同t的2160階NFOL方法的球諧綜合結(jié)果統(tǒng)計(jì)

通過以上試驗(yàn),采用2160階次模型進(jìn)行計(jì)算時(shí),如果t取0,計(jì)算結(jié)果是錯(cuò)誤的,說明t參數(shù)不能省略,取t=110可保證地面任意點(diǎn)位10-19m·s-2量級(jí)的計(jì)算精度,因此后續(xù)實(shí)驗(yàn)采用t=110. 同時(shí)發(fā)現(xiàn),當(dāng)緯度低于45°的時(shí)候,利用NFOL方法并沒有明顯的優(yōu)勢(shì),緯度為60°～70°時(shí)利用NFOL方法效率提升近1倍.

3.2 南極地區(qū)(高緯度)六邊形網(wǎng)格重力場(chǎng)計(jì)算

本文基于張永生等(2007)所采用的基于二十面體的Synder等積投影(Synder,1992)及相應(yīng)剖分方法生成ISEA4H(Icosahedral Snyder Equal Area aperture 4 Hexagon)網(wǎng)格系統(tǒng),其中南極區(qū)域內(nèi)包含17732個(gè)非等緯度分布的六邊形網(wǎng)格點(diǎn),單個(gè)網(wǎng)格平均實(shí)地面積780 km2,網(wǎng)格平均間距約15.76 km. 這里利用2160階次EGM2008模型,采用全次和NFOL兩種方法計(jì)算網(wǎng)格點(diǎn)處的重力異常值,驗(yàn)證NFOL方法在高緯度區(qū)域球諧綜合中的計(jì)算效能和精度,其計(jì)算結(jié)果和精度統(tǒng)計(jì)見表2和圖11、圖12.

表2 NFOL方法求解南極地區(qū)17732個(gè)點(diǎn)重力異常的性能統(tǒng)計(jì)(單位：10-5m·s-2)Table 2 Performance of the non-full order Legendre method for calculating 17732 point gravity anomalies in Antarctica (Unit: 10-5m·s-2)

圖11 南極地區(qū)17732個(gè)六邊形網(wǎng)格點(diǎn)重力異常分布Fig.11 Distribution of gravity anomalies at 17732 hexagonal grid points in Antarctica

圖12 南極地區(qū)六邊形網(wǎng)格點(diǎn)重力異常NFOL解的精度Fig.12 Accuracy of gravity anomalies calculated by the non-full order Legendre method on hexagonal grids in Antarctica

由以上圖表可見,NFOL方法能夠保證平均10-19m·s-2量級(jí)精度水平的重力異常球諧綜合,計(jì)算耗時(shí)不到常規(guī)全階次求解耗時(shí)的一半,具有明顯的性能提升. 需要注意的是,當(dāng)利用NFOL方法計(jì)算緯度小于45°的區(qū)域重力異常時(shí),將失去其計(jì)算優(yōu)勢(shì). 為了將NFOL方法應(yīng)用于中低緯度區(qū)域,我們結(jié)合SHR方法,構(gòu)建位于赤道附近的加里曼丹島區(qū)域的六邊形高分辨率網(wǎng)格重力場(chǎng).

3.3 加里曼丹島區(qū)域(低緯度)六邊形網(wǎng)格重力場(chǎng)計(jì)算

利用高分辨率的ISEA4H網(wǎng)格結(jié)構(gòu)剖分加里曼丹主島區(qū)域,生成15125個(gè)非等緯度分布的六邊形網(wǎng)格點(diǎn),單個(gè)網(wǎng)格平均實(shí)地面積48.75 km2,網(wǎng)格平均間距約3.94 km.

首先確定加里曼丹島區(qū)域中心的地心坐標(biāo)(θ=89.141°,λ=114.2°),然后根據(jù)2.4節(jié)所述方法,將北極點(diǎn)旋轉(zhuǎn)至該中心點(diǎn),利用2160階次EGM2008模型位系數(shù),通過(38)式的旋轉(zhuǎn)變換,得到新的球諧系數(shù)EGM2008R,該轉(zhuǎn)換過程總耗時(shí)不足500 s,通?？筛鶕?jù)計(jì)算區(qū)域的位置提前計(jì)算.在旋轉(zhuǎn)變換后的坐標(biāo)系中,加里曼丹島的六邊形網(wǎng)格點(diǎn)均位于新北極點(diǎn)附近,故均屬于“高緯”,因此使用NFOL方法,以EGM2008R系數(shù)和旋轉(zhuǎn)后的坐標(biāo)作為輸入,計(jì)算15125個(gè)點(diǎn)的重力異常值,并與利用EGM2008模型和旋轉(zhuǎn)前的坐標(biāo)、基于常規(guī)全階次方法的計(jì)算結(jié)果進(jìn)行比較,結(jié)果見表3以及圖13和圖14.

表3 SHR+NFOL方法求解加里曼丹島地區(qū)15125個(gè)點(diǎn)重力異常的性能統(tǒng)計(jì)(單位：10-5m·s-2)Table 3 Performance of the NFOL combined with SHR for calculating 15125 point gravity anomalies in the Kalimantan area (Unit: 10-5m·s-2)

圖13 加里曼丹島地區(qū)15125個(gè)六邊形網(wǎng)格點(diǎn)基于SHR+NFOL方法計(jì)算的重力異常結(jié)果Fig.13 The results of gravity anomalies calculated by the non-full order Legendre method combined with SHR at 15125 hexagonal grid points in Kalimantan Island area

圖14 加里曼丹島地區(qū)15125個(gè)六邊形網(wǎng)格點(diǎn)基于SHR+NFOL方法計(jì)算重力異常的精度情況Fig.14 The accuracy of gravity anomalies calculated by the non-full order Legendre method combined with SHR at 15125 hexagonal grid points in Kalimantan Island area

本次計(jì)算的重力異常的精度水平為10-16m·s-2,比3.2節(jié)中直接使用NFOL方法求解精度水平10-19m·s-2低,說明SHR由于采用坐標(biāo)變換和球諧系數(shù)的變換,引入了計(jì)算誤差,該誤差導(dǎo)致重力異常的計(jì)算誤差約10-16m·s-2; 同時(shí)可以發(fā)現(xiàn),此次實(shí)驗(yàn)中，NFOL方法計(jì)算效率是全階次解的近5倍,加速效果更加明顯,這說明,相比3.2節(jié)中的大區(qū)域、低分辨率六邊形網(wǎng)格點(diǎn)重力異常的球諧綜合,NFOL方法在小區(qū)域、高分辨率離散點(diǎn)的球諧綜合中更具優(yōu)勢(shì).

4 結(jié)論

為實(shí)現(xiàn)局部六邊形網(wǎng)格重力場(chǎng)元模型值的快速計(jì)算,以南極洲和加里曼丹島區(qū)域的六邊形網(wǎng)格點(diǎn)重力異常計(jì)算為例,研究了NFOL和SHR方法在球諧綜合快速計(jì)算中的應(yīng)用.

首先,對(duì)ALFs遞推過程及其數(shù)值特征進(jìn)行了詳細(xì)分析,給出了其值從穩(wěn)定振蕩區(qū)到快速衰減區(qū)的分界線的理論表達(dá)式,利用分界線表達(dá)式對(duì)“次”進(jìn)行有效截?cái)?聯(lián)合跨階次遞推,提出了基于NFOL方法的高緯度點(diǎn)快速球諧綜合. 在低緯度區(qū)域六邊形網(wǎng)格重力異常計(jì)算方面,首次將SHR與NFOL方法結(jié)合,通過坐標(biāo)系旋轉(zhuǎn),將低緯度區(qū)域的點(diǎn)變換到“高緯度”區(qū)域,從而利用NFOL方法實(shí)現(xiàn)低緯度區(qū)域的快速球諧綜合. 通過變量變換理論推導(dǎo),借助X-數(shù)法,實(shí)現(xiàn)了2160階次球諧系數(shù)的任意角度旋轉(zhuǎn)下的變換.

通過數(shù)值實(shí)驗(yàn)及分析,我們得到以下結(jié)論：對(duì)于2160階次的Legendre遞推,若保證全球任意點(diǎn)處計(jì)算重力異常的精度滿足10-19m·s-2,ALFs從穩(wěn)定振蕩區(qū)到快速衰減區(qū)的分界線理論表達(dá)式中t的最佳取值為110; 在高緯度區(qū)域(≥60°),NFOL方法能夠平均提升ALFs求解和球諧綜合的效率近1倍,緯度越高,效率提升越明顯,最高可提速近10倍,且在雙精度數(shù)值運(yùn)算中,能夠保證重力異常10-19m·s-2的平均精度水平,基于此構(gòu)建了南極洲地區(qū)的六邊形網(wǎng)格重力異常; 在中低緯度區(qū)域(<60°),NFOL方法借助SHR,同樣能夠?qū)崿F(xiàn)加速效果,其重力異常計(jì)算精度在10-16m·s-2水平,精度足夠滿足應(yīng)用需求，且本文提出的SHR+NFOL方法在進(jìn)行小區(qū)域高分辨率大量離散點(diǎn)的球諧綜合應(yīng)用中具有明顯優(yōu)勢(shì).

需要強(qiáng)調(diào)說明,NFOL方法是對(duì)ALFs計(jì)算效率的提升,而對(duì)于非等緯度點(diǎn)球諧綜合的計(jì)算效率提升方法還有很多,例如采用Seif等(2018)提出的Legendre內(nèi)插方法,能夠更快地實(shí)現(xiàn)六邊形網(wǎng)格重力異常的計(jì)算,而本文提出的NFOL方法與其他提升球諧綜合效率的方法并不沖突(FFT方法除外,因?yàn)槠洳恍枰?jì)算ALFs),加速效果可疊加.

致謝感謝戰(zhàn)略支援部隊(duì)信息工程大學(xué)李?yuàn)檴櫧淌趯?duì)文章提出的指導(dǎo)性修改意見,賁進(jìn)、童曉沖教授在生成六邊形網(wǎng)格方面給予的幫助和指導(dǎo).