基于调和分析法的全球地形球谐系数模型构建

单建晨，李姗姗，范 雕，李新星，黄 炎，邢志斌

（1.信息工程大学，河南 郑州 450051；2.航天工程大学，北京 102200）

地形球谐系数模型是指描述地球地形地貌的一组基本参数的集合，是对全球地形数据的解析表达。相较于数值模型，地形球谐系数模型具有可以表达全球地形频谱特性的优点，并且依托地球物理公式关系便于实现与重力位、扰动位[1]和地球密度的转换运用，同时避免了在获取特定区域、点位物理量时计算积分、内插等复杂过程。目前地形球谐系数模型在地形正向建模、布格异常或重力扰动计算、球谐地形引力位模型构建和未测区域地形数据预测等方面拥有实际运用及广阔前景[2,3]。

通常球谐系数模型的求解方法可分为调和分析法（Spherical Harmonic Analysis，SHA）和最小二乘法（Least-squares，LSQ）。最小二乘法能够提供求算球谐系数的参数方差、协方差信息，但解算甚（超）高阶球谐系数时，需要大量的计算资源对超大型非对角法方程矩阵进行求逆，对计算机性能提出了非常高的要求。调和分析法的出发点是对积分计算的数值评估，计算易于实现，但因为离散勒让德函数的非正交性使得精度欠佳，为了提高调和分析法求解模型精度国内外学者做了大量研究[4,5]。

本文研究推导了现有调和分析方法，以国际发布地形球谐系数模型Earth2012 球谐综合计算的地形值作为输入数据，分别利用Gauss-Legendre 积分法和矩形离散积分法完成了2160 阶地形球谐系数模型的构建，探讨分析了各个方法的优劣，着重分析了不同积分法恢复的模型在高纬度区域的表达精度，并提出了联合快速傅里叶变换（FFT）技术和基于OpenMP 多核并行技术的模型系数计算策略，大大提高了计算效率。

1 原理与方法

1.1 矩形离散积分法

根据展开理论，一个任意的球面上的函数可以展开为面球谐函数，因此全球地形可由球谐系数表达成关于球心经纬度的完全正规谐函数形式：

其中f(θ,λ)为全球地形函数，θ为以北极方向为起始方向的极距，范围为(0,π)，λ为球心经度，为归一化缔和Legendre 函数，、为全球地形完全正规谐函数系数。

利用正交关系，可以依据式(1)得到、的表达式：

通常现实情况下，f(θ,)λ表达为等间隔的离散数值模型，与式(2)中数值积分计算的连续函数要求相悖。调和分析法利用f(θ,)λ离散值代入求和逼近数值积分，从而避免了获取连续全球地形函数的难题。通过调和分析法所得面球谐系数存在最大（截断）阶次nmax，nmax与输入的全球地形数据分辨率R存在关系nmax=180°/R。

一种计算处理式(2)的方法是将关于球心经纬度的二重积分转换成全球区域等角格网数据的离散求和计算：

其中，n=0,1…nmax，m=0,1…n；、为格网中心点的极距和球心经度；Δθ=π/nmax，Δλ=2π/2nmax。

1.2 Gauss-Legendre 积分法

由于离散勒让德函数的非（弱）正交性，按照矩形离散积分方法处理积分式将会导致较大的精度损失，相比之下，Gauss-Legendre 积分法通过设置Gauss-Legendre 权重保证了从连续函数过渡到离散函数中勒让德函数正交性。该方法引入数值分析中Gauss-Legendre 积分的概念对数据进行重新采样、赋权并求和。区别于上文方法中的等间距网格（nmax×2nmax）的构建，Gauss-Legendre 积分法依据最大（截断）阶数nmax建立高斯格网（（nmax+1）×（2nmax+1）），该格网的纬度圈为满足方程Pnmax+1(cosθi)=0（θi为极矩）的Gauss-Legendre 节点，格网子午圈则等间距分布；同时为了便于编程实现对离散求和的公式进行了改进，将式(2)中的面球谐系数展开至经纬度方向分步拆解计算，对于经度方向：

其中m=0,1,2…,∞。应用离散求和法计算式(4)中定积分近似值：

其中m=0,1,2…,nmax。对于纬度方向：

其中n=m,m+1,m+2…,∞。运用数值分析中Gauss-Legendre 方法对式(6)离散求积：

其中，n=m,m+1,m+2…∞；式(5)(7)中δm0为克罗内克尔符号；式(7)中wi为Gauss-Legendre 权，定义为：

其中i=0,1,2…nmax。需要特别指出的是式(8)分母中(cosθ)为勒让德多项式对cosθ的一阶导数，并非对θ的一阶导数。

2 试验结果与分析

2.1 数据准备和处理

本文选取科廷大学发布的Earth2012.topo_bathy（Earth2012 系列模型之一）2160 阶地形球谐系数模型作为试验数据，该球谐模型是基于SRTM V4.1、SRTM30_PLUS、Etopo1 等数据源所构建的5 ′×5 ′数值模型并做球谐分析得到。为保证用于Gauss-Legendre 和矩形离散积分法两种球谐分析的输入数据具有一致性，试验利用Earth2012.topo_bathy模型球谐综合完成所构建格网的数据采样，并将Earth2012.topo_bathy 模型作为真值，以检验两种方法所得地形球谐系数的还原效果。试验采用改进的Belikov 方法计算Legendre 函数，可以较好地保证函数计算精度与稳定性[6]。在编程实现球谐综合计算时，使用OpenMP 并行算法和FFT 技术，实现了快速运算。

2.2 Gauss-Legendre 格网纬度圈计算

根据Gauss-Legendre 格网构建要求，格网纬度圈取值为函数(x)在区间x∈[-1,1]的零点（nmax+1个）。试验使用搜索迭代法（scan-iteration method,SIM）求算方程解[7]。因为函数(x)关于原点对称，为减少计算量只需要考虑求算函数在区间x∈(0,1]的零点。SIM 计算步骤如下：

（1）设定区间宽度w，根据给定区间宽度，通过区间端点函数值的异号性判断是否含有零点，并标记含零点区间；

（2）不断缩小搜索区间宽度，直至 Int[nmax/2]个零点独立分布在标记区间；

（3）在每一个标记区间使用割线迭代法计算零点值。

2.3 计算效率

综合上述式(5)(7)，Gauss-Legendre 积分法计算球谐系数的实质运算公式为：

其中，n=0,1…nmax、m=0,1…n。若直接按照式(9)累加的方式串行编译实现球谐系数的求解，计算效率低下，计算耗时约为121.82h，难以适应应用需求，实际试验中分别采用FFT 技术和OpenMP 并行算法对式(5)(7)中的Am(θ)、Bm(θ)与Cnm、Snm的计算求解。

式(5)为卷积表达形式，可以按照式(10)对Am(θ)、Bm(θ)进行快速傅里叶变换处理，其中Re、Im 分别代表变换所得实部与虚部[8,9]。

其中m=0,1,2…,nmax。考虑式(7)中各阶次待求球谐系数仅仅与唯一变量极矩θ相关，因而考虑以极矩θ为单元并行计算，此种并行方式还可以将勒让德函数计算次数降到最低，即一个采样纬度圈上的所有采样点只需计算一次勒让德函数。为避免内存位置读写冲突，利用pravite()子句对简单变量实现线程私有化处理，对勒让德函数[n][m]、s球谐系数Cnm[n][m]、Snm[n][m]等数组升维成[v][n][m]、Cnm[v][n][m]、Snm[v][n][m]，通过对v的控制，实现各线程中对这些数组的独立访问与运算[10]。

试验所用计算器处理器为Intel®Core™ i7-9750H CPU @ 2.60 GHz，内存32.00 GB，采用VS2013 编程平台10 个线程并行计算，总耗时466.083 s，相比于串行累加算法耗时，加速比约为983，较好地解决了串行编译耗时多、效率低的问题。

2.4 精度评价

按照上文原理算法，建立矩形等角格网和Gauss-Legendre 格网，并利用格网数据恢复地形球谐模型。其中Gauss-Legendre 积分法恢复模型记为Sph.GL 模型，矩形离散积分法恢复模型记为Sph.Rect模型，试验以原模型Earth2012 为标准，引入模型系数阶方差、系数阶误差等标准检验模型精度。

以Earth2012 作为标准模型，计算两个恢复模型的阶误差RMSn，公式如下：

其中，n=0,1,2…nmax；、为Earth2012模型球谐系数；、为恢复模型球谐系数。阶误差反应了恢复模型与标准模型之间的差异，差值越小，说明模型符合程度越高，恢复效果更好。

如图1所示用阶误差RMSn对模型进行对比，为便于直观比较，同时绘制了Earth2012、Sph.Rect、Sph.GL 模型的阶方差曲线。

图1 模型球谐系数阶误差Fig.1 Order error of spherical harmonic models

图1(a)-(d)显示，各模型阶方差与阶误差曲线均连续且光滑，所恢复的球谐系数模型Sph.Rect、Sph.GL的阶方差曲线在各频段与原模型Earth2012 近乎吻合，可以认为球谐分析结果正确。

比较图1(a)-(d)可以看出，模型Sph.GL 阶误差在各频段都小于模型Sph.Rect，模型Sph.GL 比模型Sph.Rect 更接近于模型Earth2012。其中图1(a)中，尽管各模型零阶球谐系数量级达到 103的情况下，计算数据显示还原模型Sph.GL 与标准模型零阶系数的阶误差达到 10-6量级。图1(a)-(c)显示模型Sph.GL 在这些频段的阶误差维持在 10-1以下，而模型Sph.Rect 的阶误差在 10-1以上。值得关注的是，在图1(d)中，约2150 阶以后，模型Sph.Rect 的阶方差和阶误差曲线均出现了急剧下降的现象。经过反复研究比较模型数据，发现直至2159 阶次，模型Sph.Rect 球谐系数的量级与模型Earth2012 保持一致，然而在2160 阶的高次项处模型Sph.Rect的系数量级出现了偏移：、量级为10-27，远小于认为不再有效，因此按照试验中矩形离散积分法所求算模型严格意义上只精确到了2159 阶。相比之下Sph.GL 模型球谐系数的量级与Earth2012 模型在全阶次保持一致，可以认为精确到了2160 阶，达到模型恢复阶数要求。

为了进一步比较分析两种方法恢复模型的精度，给出了球谐系数的误差谱，如图2所示。比较图2(a)和图2(b)，在中、低阶次部分，模型Sph.GL 球谐系数误差小于模型Sph.Rect，两者在高阶次部分差异较小。

图2 模型Sph.Rect 和模型Sph.GL 球谐系数误差谱Fig2 Coefficient error spectral of models

综合以上分析，模型Sph.GL 的阶误差较小且中、低阶次球谐系数误差明显小于模型Sph.Rect，因此相较于矩形离散积分法，Gauss-Legendre 积分法具有更好的球谐系数模型还原效果。

分别利用模型Earth2012、还原模型Sph.Rect 及Sph.GL 的球谐系数进行球谐综合计算，获得分辨率为5 ′×5 ′的全球数值模型SHS-Earth2012、SHS-Rect 和SHS-GL。以SHS-Earth2012 为标准，分别记SHS-Rect和 SHS-GL 与其差值为“ DIF-SHS-Rect”、“DIF-SHS-GL”，统计信息示于表1。

表1 球谐综合结果差值统计（单位：米）Tab.1 Difference statistics of spherical harmonic synthesis results (Unit:m)

表1显示，球谐模型GL 恢复的数值模型与Earth2012 更加接近，两者差值绝对值的均值为1.857×10-6m，远小于DIF-SHS-Rect 的4.605 m；两者的标准差差异较大，DIF-SHS-GL 差值大小分布更为集中；DIF-SHS-GL 极大、小值均较小。就所对比的各项指标而言，可以认为模型SHS-GL 均优于模型SHS-Rect，更为接近标准模型。

统计DIF-SHS-Rect、DIF-SHS-GL 两者在±500 m之间的点数量，绘制柱状分布统计图如图3所示。

图3 全球球谐模型误差分布图Fig.3 Error distribution of global spherical harmonic models

图3显示，DIF-SHS-GL 数据更集中地分布在0附近。为了更好地比较两个模型差值的分布情况，分别统计差值在±1 m、±10 m、±100 m、±500 m 区间的点数量和占比情况，如表2所示。

表2 全球球谐模型误差数值分布统计Tab.2 Numerical distribution statistics of global spherical harmonic model error

表2显示，两组差值在|δh|＜500m范围内的点数量都远超过90%。定义差值在|δh|＜1m的点为高精度还原点，比较两表，DIF-SHS-GL 的高精度点数超过DIF-SHS-Rect 5 倍，达到总点数的12.160%。同时DIF-SHS-GL 数值在|δh|＜100m范围内的点数占比超过了 90%，占到了试验点的绝大多数。认为|δh|＞ 100m 的点还原效果较差，定义为异常点，其在两个模型占比分别约为23%、10%。绘制DIF-SHS-Rect、DIF-SHS-GL 两模型中异常点位置分布图，示于图4、图5。

图4、图5中大陆架轮廓清晰可见，两模型中异常点多数沿着海陆交界处分布，笔者认为可能的原因是陆地与海洋交界处地形变化剧烈，但纬度圈上的格网点是等角距分布，所包含信息不够丰富，不能细致地描绘这些区域；DIF-SHS-GL 中异常点的数量、分布密度都明显小于DIF-SHS-Rect；在南北方向高纬度区域SHS-GL 与SHS-Earth2012 吻合地很好，对比与矩形离散积分法，Gauss-Legendre 积分方法在这些区域具有独特优势。

图4 DIF-SHS-Rect 异常点位分布图Fig.4 Distribution of anomalous points in DIF-SHS-Rect

图5 DIF-SHS-GL 中异常点位分布图Fig.5 Distribution of anomalous points in DIF-SHS-GL

为进一步比较两种方法所恢复模型在高纬度区域的精度，分别研究了模型DIF-SHS-Rect、DIF-SHS-GL在区域60°N～90°N、60°S～90°S的表现效果，北半球区域记为DIF-SHS-Rect-N、DIF-SHS-GL-N，南半球区域记为DIF-SHS-Rect-S、DIF-SHS-GL-S，两个模型高纬度区域差值分布如图6所示。比较图6(a)与图6(b)、图6(c)与图6(d)，DIF-SHS-GL 在南、北高纬度区域的值更加集中于 0m 附近，明显优于DIF-SHS-Rect。

图6 高纬度区域差值分布图Fig.6 Difference distribution diagram at high latitudes

表3为两组球谐综合结果在高纬度的差值统计情况，结合两种方法的原理，由于Gauss-Legendre 积分法赋予采样点Gauss-Legendre 权重，从而很好地改善了积分法在高纬度区域表达精度问题。如表3所示，DIF-SHS-GL 在南、北高纬度区域的标准差约为DIF-SHS-Rect 的1/6。对比表1和表3，还可以发现研究区域DIF-SHS-GL-N、DIF-SHS-GL-S 的标准差远小于全球区域的标准差20.920 m，约为其1/2，表达精度优于全球平均标准；而 DIF-SHS-Rect-N、DIF-SHS-Rect-S 的标准差比表1中全球区域标准差大，低于全球平均标准。

表3 高纬度区域球谐综合结果差值统计（单位：米）Tab.3 Difference statistics of spherical harmonic synthesis results at high latitudes (Unit:m)

表4为球谐模型误差在南、北高纬度区域数值分布统计信息，比较发现DIF-SHS-GL 在各数值区间的数量占比均明显大于DIF-SHS-Rect。对比表4中高纬度误差与表2中全球误差分布，各模型在高纬度区域数量占比均有提高，其中：DIF-SHS-Rect 和DIF-SHS-GL 在|δh|＜1m区间数量占比分别平均提高1.678%和16.444%；在|δh|＜10m 区间数量占比分别平均提高12.041%和30.474%；在|δh|＜100m 区间数量占比分别平均提高 13.037%和 8.596%。除|δh|＜100m 区间外，DIF-SHS-GL 提高更显著，但DIF-SHS-GL 在南、北高纬度区域|δh|＜100m 区间的差值数量占比均超过98%，远大于DIF-SHS-Rect。分析球谐模型表达高纬度与全球精度存在差异原因为：球面采样中高纬度采样点更密集，富含更多地形特征信息，因而恢复模型描绘高纬度区域能力更强，在各区间数量占比也会高于全球平均水平。

表4 高纬度区域差值模型数值分布统计Tab.4 Numerical distribution statistics of difference models at high latitudes

3 结 论

本文首先建立矩形等角格网和Gauss-Legendre 格网；然后依据Earth2012 模型球谐综合计算离散格网点地形数据，并分别用矩形离散积分法、Gauss-Legendre 积分法构建全球2160 阶地形球谐系数模型Sph.Rect、Sph.GL，并还原出全球5′×5 ′地形数值模型SHS-GL、SHS-Rect；最后讨论了两种方法构建地形球谐模型的精度。可得到以下结论：

（1）比对 Sph.GL 模型和 Sph.Rect 模型与Earth2012 模型的阶误差表明，在输入数据近似相等的情况下，Gauss-Legendre 积分法还原效果更好，与标准模型更接近；

（2）比对数值模型SHS-GL 和SHS-Rect 表明，在海陆交界等地形变化剧烈区域，矩形离散积分法和Gauss-Legendre 积分法描绘地形细节能力均显不足；

（3）相比矩形离散积分法，Gauss-Legendre 积分方法表达高纬度区域地形更具优势，且该方法描绘高纬度地区地形精度高于描绘全球的平均水平，能较好解决两极附近存在采样点较密、数据不均匀的问题；

（4）通过数组拓展升维等方法可以有效解决并行中内存读写冲突问题，合理的并行方案及快速傅里叶变换技术可以充分利用计算机计算资源，进而大大提高调和分析计算中耗时长、效率低的问题。可为高阶地形球谐模型构建方法选择提供参考。