李新星, 李建成, 刘晓刚, 范昊鹏, 靳超

1 武汉大学测绘学院， 武汉 430079 2 信息工程大学地理空间信息学院， 郑州 450001 3 西安测绘研究所， 西安 710054 4 61365部队， 天津 300000

0 引言

传统基于等角地理网格的地球重力场描述面临着许多无法克服的局限性，包括网格面积不等、重力场平滑因子复杂、观测数据冗余、积分离散化误差大等，而六边形的“离散球面网格系统”(Discrete Global Grid System,DGGS)因其良好的特性，在遥感、气象、水文等领域取得了广泛应用(Sahr et al.,2003). 作为地理网格的延伸，美国著名学者富勒(Fuller)最早进行了球面离散化的相关研究并且发明了“网格球顶”(Geodesic Dome)结构，之后许多学者关于球面网格系统的研究都直接或间接地受到了他的影响(张永生等,2007). 近年来，随着DGGS网格系统不断完善，其已经在诸如矢量数据结构(Dutton，1999)、空间数据索引(Dutton,1996a)、制图综合(Dutton,1996b)、全球动态数据结构(Goodchild and Shiren,1992)、全球土地监测(Suess et al.,2004)、全球海洋分析(Kidd et al.,2003)、全球气候模型(Randall et al.,2002)、海洋路径规划、城市交通规划(Brodsky,2018)等方面得到了应用.

基于六边形的DGGS网格本应在重、磁场这一类具有各向同性特征的物理场中发挥更大作用，但是目前相关研究很少. 类六边形网格结构在重磁场应用方面，Vestine等(1963a,1963b)较早研究了基于球面均匀分布点的地磁场的Poisson积分应用，Eicker(2008)将多种点位均匀分布方式引入到径向基函数中进而求解卫星重力场取得了较好结果，Slobbe等(2012)将Fibonacci均匀分布引入球形斯列普基函数求解平均海水面和海面地形，Ran等(2018)将Fibonacci格网用于Mascon方法来反演冰盖质量变化，以上均是在数值积分求解方面引入均匀分布结构，而鲜有在球谐分析与综合中开展研究，其中最主要的原因在于，基于六边形网格系统的点位并不像地理网格那样是等纬度分布的，如图1所示，从而需要大量缔合Legendre函数(associated Legendre functions,ALFs)的求解，这对计算能力提出了更高要求(Pavlis,1988; 李新星等,2014; 田家磊等,2018).

图1 基于全球六边形网格的360阶全球重力异常分布图Fig.1 360-degree gravity anomaly values based on global hexagonal grids

同时,观测资料的精度与数量的增多使得对应球谐展开的阶次也不断升高,对计算能力同样提出了挑战. 例如目前有720阶次的NGDC-720地壳磁异常模型(Maus,2010),有经典的2160阶次地球重力场模型EGM2008(Pavlis et al.,2012),还有10800阶次的地形、布格异常、均衡异常模型(Balmino et al.,2012; Hirt and Rexer,2015). 受限于球谐展开所需庞大的计算量,现有网格数据的分辨率已经远远走在现有模型阶次对应分辨率的前列. 例如地磁场的磁异常网格数据EMAG2v3分辨率达到2′,根据Nyquest采样定律,因磁场为偶极子场,其对应10800阶次的球谐展开; sandwell和DTU15版本的海洋卫星测高重力异常分辨率达到1′(Sandwell and Smith,2008; Andersen and Knudsen,2016), 因重力场为标量场,所以其对应10800阶次; GEBCO2014和GEBCO2019海底地形网格数据的分辨率分别达到了30″和15″,对应21600和43200阶次; SRTM3、ARTER和ALOS全球DEM数据分辨率分别达到90 m、30 m和12 m,对应的球谐展开阶次更高,因此,球谐展开中的严重耗时问题也是急需解决的.

针对ALFs及其相应积分、导数等函数的计算问题,许多学者从递推方法、计算精度、计算效率和溢出问题等方面做了大量的工作. 对于ALFs的递推方法,主要包括行推法、列推法、Belikov和跨阶次递推等四种方法,其中Belikov方法的计算效率最佳(李新星,2013). 为了解决递推过程中的溢出问题,Fukushima (2012)提出了X数法实现了任意阶次ALFs的递推计算,代价是损失了一定的计算效率; Xing等(2020)改进了Belikov方法,解决了Belikov原递推公式在高阶次出现溢出的情况,例如纬度为75°的ALFs当阶次达到4800时，ALFs的计算出现溢出; 张传定等(2004)首次提出ALFs的跨阶次递推方法,并经验证在20000阶次以内其计算准确且不会发生溢出(于锦海等,2015).

关于ALFs递推计算效率的进一步提升,目前还没有较为显著的方法. 学者大多绕过ALFs的计算,通过尽量减小ALFs的计算量来寻求提升球谐综合(Spherical Harmonic Synthesis,SHS)和球谐分析(Spherical Harmonic Analysis,SHA)计算效率的方法,例如在SHS和SHA中要求相应的计算点等纬度网格分布,从而计算同一纬圈上所有点的模型值只需要计算一次该纬度的ALFs. Clenshaw(1957)提出Clenshaw求和方法避免大部分Legendre函数值的求解来提高SHS的计算效率,同时快速傅里叶变换(Fast Fourier Transform,FFT)技术的使用大大加速了球谐分析与综合的计算,例如Colombo(1981)在球谐分析中使用了一维FFT,Sneeuw和Bun (1996)利用球面到轮胎面映射实现了二维FFT的球谐分析与综合. 尽管FFT技术加快了SHS的计算速度,但其仅限于在计算规则网格分布下具有相同高度的点集,针对非等纬度、不同高度分布点的SHS和SHA,有学者采用了泰勒级数展开或等间隔内插的思想来提升效率(Kunis and Potts,2003; Moazezi et al.,2016; Seif et al.,2018).

实际应用中,更多是针对局部重力场或磁场数据的处理,为了解决因球谐方法计算量庞大而不便于局部使用的问题,学者们避开球谐模型,使用局部拟合、等效源、矩谐分析、球冠谐分析等方法来进行局部物理场的建模(徐文耀和朱岗昆,1984; 高金田等,2005; Younis,2013),而上述方法存在的边界效应、基函数不正交、非整阶Legendre递推等问题还需要深入研究,如果球谐模型变得更加高效,也必然能够在局部重磁场模型构建中发挥重要作用.

本文从ALFs求解方面入手,寻求提升其计算效率的方法. Jekeli等(2007)总结了不同纬度下ALFs函数值的量级随阶次升高的衰减关系,通过利用图2a所示的逆向行推方法,求解ALFs值的有效部分,对于绝对值较小(例如<10-15)的值默认为0而不进行计算,从而节约了递推过程中的计算量,在全球SHS过程中,该方法理论上能够节约36%的计算量,但是由于其逆向行推方法计算效率偏低,并没有给出最终的计算耗时比较.

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

本文利用该思想,结合跨阶次递推方法,提出了非全次Legendre(Non-full-order Legendre,NFOL)计算方法,实现了高纬度区域任意点的快速球谐综合. NFOL方法是在利用跨阶次递推方法计算ALFs时,当递推到达ALFs有效值(例如>10-15)边界即可停止递推,从而减小计算量,提升效率,由于有效值分布在全部“阶”中“次”偏小的区域,所以该递推只需要递推全部阶的部分“次”.

对于中低纬度的ALFs计算,由于ALFs有效值区域几乎遍布整个阶次,所以NFOL方法在低纬度区域并不适用,针对该问题,本文引入球谐旋转变量变换(Spherical Harmonic Rotation,SHR)方法,该方法实现了坐标系旋转下,对球谐展开系数的变换. 利用该方法将低纬度待求解区域的中心通过坐标系旋转移至北极点,从而使得旋转后的点位均位于“高纬度”,使得与之对应的ALFs有效值区域集中于低“次”部分,从而可采用NFOL方法提升SHS效率.

SHR方法在多领域具有重要作用. 地磁场模型中,该方法对于不同坐标系下磁场模型球谐系数间的比较、日变磁场Sq的内/外源场分离等具有重要作用(De Santis et al.,1996). 该方法求解过程中,Wigner D矩阵起到了重要的作用,许多学者研究了D矩阵及其递推求解方法(Aubert,2013; Fukushima,2017).国内,许厚泽和蒋福珍(1964)较早研究了SHR原理并推导了公式,其在飞行器轨道扰动引力快速赋值中得到应用(任萱,1985; 郑伟等,2011; 王建强等,2013),但是这些应用中采用的球谐系数阶次较低,且其变换精度均没有经过严格的精度评估. Fukushima (2017)利用X数法实现了2160阶次球谐系数的90°旋转.

本研究基于NFOL和SHR方法,分别实现了高纬度南极地区低分辨率的六边形网格点重力异常、低纬度加里曼丹岛地区的高分辨率六边形网格点重力异常的构建,通过比较分析,验证了方法的精度与效率,得到了有益结果,成果值得推广.

1 非全次Legendre函数的递推

1.1 球谐展开基本公式

对于球面上的平方可积函数,通过SHA能够得到对应的球谐展开系数,这里以重力场中重力异常的球谐展开为例,其与重力场位系数之间的关系如下：

(1)

(2)

(3)

(4)

n≥2,n≠m

(5)

其中

(6)

递推初始值为

(7)

当m≥2时,采用以下递推公式(于锦海等,2015)：

(8)

其中

(9)

(10)

(11)

当m>2时,递推系数cn m,dn m,hn m的值均小于1,因此跨阶次递推是稳定可靠的,递推过程如图2b所示.

1.2 非全次Legendre方法

图3 m=100时,θ=1°、45°和89°的ALFs的大小Fig.3 The value of ALFs with θ=1°,45° and 89° when m=100

m≈nsinθ.

(12)

图4 ALFs值中稳定振荡区与快速衰减区的分界线,图中绿色区域值较小,近似为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的值在该分界线周边有一个衰减过程,如图5所示,因此想要保证整个ALFs的计算精度,不能直接以式(12)所给出的直线作为分界线,而应该适当地向m增大和n减小的方向移动该直线,得到新的分界线

m≈nsinθ+t,

(13)

图5 分界线周边ALFs值的衰减过程,图中蓝色区域值小于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

其中参数t的大小决定了ALFs递推值的精度和计算效率,t越大,精度越有保证,但ALFs需要递推的区域会变大(见图6),相应计算效率会降低,关于参数t的确定,会在后续实验部分给出.

图6 参数t的含义Fig.6 The meaning of parameter t

因此,上述方法确定了ALFs的递推范围,我们开创性地避开Jekeli等(2007)所采用的逆向行推方法,采用跨阶次递推方法计算了非全次递推范围内的ALFs值,解决了逆向行推不稳定且计算效率低下的问题.由图5可知,对于非高纬度位置处的ALFs,显然上述方法并不能有效提升计算效率,因此引入球谐系数的旋转变量变换,通过换极将计算区域移至“高纬”区域,从而提升计算效率.

2 球谐旋转变量变换(SHR)

2.1 沿z轴的旋转变换

已知球面上一点Q(r,θ,λ),当坐标系O-xyz沿着z轴旋转α角,记为Rz(α),得到新的坐标系O-x′y′z,则Q点在新坐标系下的球坐标为(r,θ,λ1),其中λ1=λ-α,即λ=λ1+α，如图7所示.

图7 沿z轴的坐标系旋转Fig.7 Rotation of the coordinate system along axis z

因此Q点的重力异常在新的球坐标系O-x′y′z下可展开为新的球谐系数

(14)

根据λ与λ1的以下三角关系式：

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

-sinmλ1sinmα，

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

+cosmλ1sinmα，

(15)

代入式(1)并与式(14)进行比较,可得两个坐标系下球谐系数之间的关系

(16)

式(16)即为坐标系沿z轴旋转α角时,对应的球谐系数的变换公式.

2.2 沿y轴的旋转变换

已知球面上一点Q(r,θ,λ),当坐标系O-xyz沿着y轴旋转β角,记为Ry(β),得到新的坐标系O-x″yz″.

如图8,Q点在新坐标系下的球坐标为(r,Θ,Λ),其中Θ为新坐标系下的球心余纬,Λ为新坐标系下的球心经度,Λ1与Λ互补,Q点的重力异常在新的球坐标系O-x″yz″下可展开为新的球谐系数,表达式与式(1)类似，

图8 沿y轴的坐标系旋转Fig.8 Rotation of the coordinate system along axis y

(17)

(18)

(19)

(20)

(21)

(22)

其中

(23)

(24)

(25)

(26)

(26)式可利用D函数的对称性及其与该两组变量之间的对应关系得证,具体见Aubert(2013),即图9d中k=4与图9c中m=4对应的元素值相等.

图10 β=90°的30阶值的切片图(a) 30阶球谐转换矩阵的值; (b) 与比较说明式(29)对称结构.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°的旋转变换

(27)

其中

(28)

(29)

(30)

(31)

(32)

其中

(33)

(34)

(35)

上述递推方法的计算过程中,Hn,m随着n的增大快速减小,当其绝对值小于2.225×10-308的时候,数值发生溢出而强制为0,会影响后续递推工作,因此在该递推过程中应使用X-数法进行求解.

2.4 将(θ,λ)视为北极点的旋转变换

Zotter (2009)文中指出,Ry(θ)可分解为以下旋转的组合：

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

(36)

还可以写为

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

(37)

因此,只需要计算沿y轴旋转90°的球谐系数变换以及简单的z轴变换,就可以实现任意角度的旋转变换,即

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

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

(38)

其中系数的变换采用公式(16)和(18).

3 数值实验及分析

为了验证NFOL和SHR方法在局部六边形网格重力场快速球谐综合方面的应用,共进行三个实验,计算环境如下：2.30GHz主频的Intel Core i5-6200U CPU和8 GB主内存的PC机,64位Windows XP OS操作系统下VS2017编译环境.

3.1 参数t的确定

第1节中指出,使用NFOL方法进行球谐综合,需要确定公式(13)中的参数t,t取值越大,计算精度越高而对应效率越低,反之亦然. 因此,需要在满足精度要求的情况下,合理的选择t的大小来提高球谐综合的计算效率.

利用2160阶2159次的EGM2008模型位系数,计算地球表面点的模型重力异常值,利用不同纬度、不同t取值的NFOL方法计算重力异常,并与常规全阶次方法计算的结果作为真值进行比较,结果统计见表1.

表1 不同纬度、不同t的2160阶NFOL方法的球谐综合结果统计

通过以上试验,采用2160阶次模型进行计算时,如果t取0,计算结果是错误的,说明t参数不能省略,取t=110可保证地面任意点位10-19m·s-2量级的计算精度,因此后续实验采用t=110. 同时发现,当纬度低于45°的时候,利用NFOL方法并没有明显的优势,纬度为60°～70°时利用NFOL方法效率提升近1倍.

3.2 南极地区(高纬度)六边形网格重力场计算

本文基于张永生等(2007)所采用的基于二十面体的Synder等积投影(Synder,1992)及相应剖分方法生成ISEA4H(Icosahedral Snyder Equal Area aperture 4 Hexagon)网格系统,其中南极区域内包含17732个非等纬度分布的六边形网格点,单个网格平均实地面积780 km2,网格平均间距约15.76 km. 这里利用2160阶次EGM2008模型,采用全次和NFOL两种方法计算网格点处的重力异常值,验证NFOL方法在高纬度区域球谐综合中的计算效能和精度,其计算结果和精度统计见表2和图11、图12.

表2 NFOL方法求解南极地区17732个点重力异常的性能统计(单位：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 南极地区17732个六边形网格点重力异常分布Fig.11 Distribution of gravity anomalies at 17732 hexagonal grid points in Antarctica

图12 南极地区六边形网格点重力异常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量级精度水平的重力异常球谐综合,计算耗时不到常规全阶次求解耗时的一半,具有明显的性能提升. 需要注意的是,当利用NFOL方法计算纬度小于45°的区域重力异常时,将失去其计算优势. 为了将NFOL方法应用于中低纬度区域,我们结合SHR方法,构建位于赤道附近的加里曼丹岛区域的六边形高分辨率网格重力场.

3.3 加里曼丹岛区域(低纬度)六边形网格重力场计算

利用高分辨率的ISEA4H网格结构剖分加里曼丹主岛区域,生成15125个非等纬度分布的六边形网格点,单个网格平均实地面积48.75 km2,网格平均间距约3.94 km.

首先确定加里曼丹岛区域中心的地心坐标(θ=89.141°,λ=114.2°),然后根据2.4节所述方法,将北极点旋转至该中心点,利用2160阶次EGM2008模型位系数,通过(38)式的旋转变换,得到新的球谐系数EGM2008R,该转换过程总耗时不足500 s,通常可根据计算区域的位置提前计算.在旋转变换后的坐标系中,加里曼丹岛的六边形网格点均位于新北极点附近,故均属于“高纬”,因此使用NFOL方法,以EGM2008R系数和旋转后的坐标作为输入,计算15125个点的重力异常值,并与利用EGM2008模型和旋转前的坐标、基于常规全阶次方法的计算结果进行比较,结果见表3以及图13和图14.

表3 SHR+NFOL方法求解加里曼丹岛地区15125个点重力异常的性能统计(单位：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 加里曼丹岛地区15125个六边形网格点基于SHR+NFOL方法计算的重力异常结果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 加里曼丹岛地区15125个六边形网格点基于SHR+NFOL方法计算重力异常的精度情况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

本次计算的重力异常的精度水平为10-16m·s-2,比3.2节中直接使用NFOL方法求解精度水平10-19m·s-2低,说明SHR由于采用坐标变换和球谐系数的变换,引入了计算误差,该误差导致重力异常的计算误差约10-16m·s-2; 同时可以发现,此次实验中，NFOL方法计算效率是全阶次解的近5倍,加速效果更加明显,这说明,相比3.2节中的大区域、低分辨率六边形网格点重力异常的球谐综合,NFOL方法在小区域、高分辨率离散点的球谐综合中更具优势.

4 结论

为实现局部六边形网格重力场元模型值的快速计算,以南极洲和加里曼丹岛区域的六边形网格点重力异常计算为例,研究了NFOL和SHR方法在球谐综合快速计算中的应用.

首先,对ALFs递推过程及其数值特征进行了详细分析,给出了其值从稳定振荡区到快速衰减区的分界线的理论表达式,利用分界线表达式对“次”进行有效截断,联合跨阶次递推,提出了基于NFOL方法的高纬度点快速球谐综合. 在低纬度区域六边形网格重力异常计算方面,首次将SHR与NFOL方法结合,通过坐标系旋转,将低纬度区域的点变换到“高纬度”区域,从而利用NFOL方法实现低纬度区域的快速球谐综合. 通过变量变换理论推导,借助X-数法,实现了2160阶次球谐系数的任意角度旋转下的变换.

通过数值实验及分析,我们得到以下结论：对于2160阶次的Legendre递推,若保证全球任意点处计算重力异常的精度满足10-19m·s-2,ALFs从稳定振荡区到快速衰减区的分界线理论表达式中t的最佳取值为110; 在高纬度区域(≥60°),NFOL方法能够平均提升ALFs求解和球谐综合的效率近1倍,纬度越高,效率提升越明显,最高可提速近10倍,且在双精度数值运算中,能够保证重力异常10-19m·s-2的平均精度水平,基于此构建了南极洲地区的六边形网格重力异常; 在中低纬度区域(<60°),NFOL方法借助SHR,同样能够实现加速效果,其重力异常计算精度在10-16m·s-2水平,精度足够满足应用需求，且本文提出的SHR+NFOL方法在进行小区域高分辨率大量离散点的球谐综合应用中具有明显优势.

需要强调说明,NFOL方法是对ALFs计算效率的提升,而对于非等纬度点球谐综合的计算效率提升方法还有很多,例如采用Seif等(2018)提出的Legendre内插方法,能够更快地实现六边形网格重力异常的计算,而本文提出的NFOL方法与其他提升球谐综合效率的方法并不冲突(FFT方法除外,因为其不需要计算ALFs),加速效果可叠加.

致谢感谢战略支援部队信息工程大学李姗姗教授对文章提出的指导性修改意见,贲进、童晓冲教授在生成六边形网格方面给予的帮助和指导.