数字高程模型的资料（信息工程大学单建晨）

数字高程模型的资料（信息工程大学单建晨）摘要：针对10 800阶地形球谐系数模型构建的计算精度、计算稳定性及超大规模运算量等问题 本文首先通过比较分析矩形离散积分方法、Gauss-Legendre积分方法和Driscoll/Healy积分方法 验证了Driscoll/Healy积分方法具有更高的计算精度。然后 给出了改进的Belikov公式 将完全正常化缔和勒让德函数在全球范围内以优于10-12精度高效、稳定递推至10 800阶; 提出了联合FFT和基于OpenMP的多核并行方法求解超高阶球谐系数的优化策略 显著提高了球谐系数模型构建的计算效率。最后 利用Earth2014_TBI与全球海底地形模型STO_IEU2020格网数据建立了全球10 800阶地形球谐系数模型sph.10 800_IEU 总体精度与Earth2014_TBI2014.shc相当 在试验海域的相对精度略优于模型Earth2014_TBI20

本文内容来源于《测绘学报》2023年第

超高阶全球地形球谐系数模型的构建方法与实现

单建晨^{1 2}

李姗姗¹ 范雕¹ 李新星¹ 黄炎¹ 邢志斌³

1. 信息工程大学 河南 郑州 450001;

2. 31201部队 广东 广州 510080;

3. 航天工程大学 北京 102200

基金项目：国家自然科学基金(42174007; 42204009)

摘要：针对10 800阶地形球谐系数模型构建的计算精度、计算稳定性及超大规模运算量等问题 本文首先通过比较分析矩形离散积分方法、Gauss-Legendre积分方法和Driscoll/Healy积分方法 验证了Driscoll/Healy积分方法具有更高的计算精度。然后 给出了改进的Belikov公式 将完全正常化缔和勒让德函数在全球范围内以优于10^-12精度高效、稳定递推至10 800阶; 提出了联合FFT和基于OpenMP的多核并行方法求解超高阶球谐系数的优化策略 显著提高了球谐系数模型构建的计算效率。最后 利用Earth2014_TBI与全球海底地形模型STO_IEU2020格网数据建立了全球10 800阶地形球谐系数模型sph.10 800_IEU 总体精度与Earth2014_TBI2014.shc相当 在试验海域的相对精度略优于模型Earth2014_TBI2014.shc。

关键词：地形球谐系数模型 Gauss-Legendre积分法 Driscoll/Healy积分法 矩形离散积分法 FFT OpenMP

引文格式：单建晨 李姗姗 范雕 等. 超高阶全球地形球谐系数模型的构建方法与实现[J]. 测绘学报，2023，52(5)：748-759. DOI: 10.11947/j.AGCS.2023.20220003

SHAN Jianchen LI Shanshan FAN Diao et al. Design and implementation of ultra-high-degree spherical harmonic model of earth topography[J]. Acta Geodaetica et Cartographica Sinica 2023 52(5): 748-759. DOI: 10.11947/j.AGCS.2023.20220003

阅读全文：http://xb.chinasmp.com/article/2023/1001-1595/20230506.htm

引 言

地形球谐系数模型可以有效弥补数值地形模型表达频谱信息的不足，并且能够转换运用到地形重力位、扰动位模型的构建，可应用于可视化地球、大规模地球物理地质等地球科学研究中。当前国际公开发布的全球地形球谐系数模型有美国地理空间情报局(National Geospatial-Intelligence Agency)的2160阶模型DTM2006.0和澳大利亚科廷大学(Curtin University)2160阶模型Earth2012、10 800阶模型Earth2014。模型Earth2012模型、Earth2014均包含了多个版本，以模型Earth2014为例，它囊括了SUR(the physical surface)、BED(bedrock)、TBI(topography bedrock and ice)、ICE(ice sheet thickness)、REF(rock-equivalent topography)等多个球谐模型，满足了不同研究应用中对全球地形表面形态的需求。根据发布球谐模型的频谱特性，文献[1—2]利用模型DTM2006.0构建高分辨率剩余地形模型(residual terrain model，RTM)，实现了崎岖地形中遗漏误差的估算与EGM2008高程异常模型精度的提高。文献[3]依据模型Earth2014岩石地形模型推导了地形重力位的球谐系数，并与GOCE卫星重力测量飞行任务观测到的重力数据进行比对，该模型被欧洲空间局(European Space Agency，ESA)成功地用于更新GOCE用户工具箱。

球谐分析可将球面观测物理量表达成球谐函数的级数形式，是地球物理学和大地测量学中重要的球面数据处理工具，在重力场^[4]、磁场^[5]、地形^[6]和电离层^[7-8]等模型的构建与研究中有着广泛的应用。常用的球谐分析方法有数值积分方法和最小二乘估计方法^[9](least square estimation，LSE)。数值积分方法的本质是关于格网的加权求和，运算总量为o(N³)阶(N为球谐分析最大阶)且计算易于实现，数值积分方法的缺陷主要在于不能顾及输入数据的观测误差，且数值逼近过程中离散勒让德函数求和的非正交性严重影响计算效果。LSE方法可以提供计算结果的参数方差、协方差等信息，但存在大规模矩阵乘法和矩阵求逆的问题，总运算量为o((N 1)⁴)阶，对于超高阶模型构建需要大量的计算资源支撑。总体而言，在模型构建计算精度方面，LSE方法相较于数值积分方法更高；在计算速度方面，数值积分方法则效率更高。综合权衡计算精度、计算效率、建模阶次及地形测量数据本身精度，本文选择数值积分方法构建10 800阶地球地形球谐系数模型。为提高数值积分方法的计算精度，文献[10—11]分别探讨了球谐分析中的Gauss-Legendre积分方法和Driscoll/Healy积分方法，较好地保证了从连续函数过渡到离散函数中勒让德函数正交性^[10-11]。文献[12]利用B样条插值技术构造出轮胎面上的连续函数，解决了离散误差的平滑因子问题。文献[13]给出了基于面积平均值的ba-DH、ba-GLQ算法，进一步提高了数值积分方法的计算精度。基于以上讨论的数值积分法，文献[14]构建了系列全球重力场模型MOD99，最高阶达到1800阶。文献[15]构建了2160阶重力场模型UGM08。文献[16]利用Mark A. Wieczorek团队研发的SHTOOLS package v2.8程序包中Gauss-Legendre积分法构建了10 800阶地形球谐系数模型Earth2014，为保证计算精度，该模型采用了X-数法计算完全正常化缔和勒让德函数(fully normalized associated Legendre function fnALF)^[17]，但文献[18]指出改进的Belikov方法与X-数法均可以将fnALF稳定递推到超高阶，但改进的Belikov方法的计算速度更优势，因此超高阶球谐系数模型构建效率仍存在提升空间。

由于超高阶勒让德函数的计算精度与稳定性、超大运算量及计算效率等问题，目前国内涉及构建10 800阶地形球谐模型的研究较少。为实现10 800阶地形球谐系数模型的高效构建，本文深入研究了Gauss-Legendre和Driscoll/Healy两种不同积分加权求和方法，并依托Earth2014地形模型设计试验，对比验证了两种方法模型构建的计算精度。针对超高阶球谐系数解算的效率低、内存需求大等问题，采用了更为高效的Belikov方法计算fnALF，同时给出了详细的计算优化策略，显著提高了计算速率，降低了对计算机性能的要求。最终利用精度更高的Driscoll/Healy积分方法，结合最新海底地形模型STO_IEU2020数据，构建了全球10 800阶球谐系数模型sph.10 800_IEU，并选取试验海域的实测水深数据进行了模型精度检核。

1 原理与方法

1.1 数值积分法中离散勒让德函数非正交性

任意一个有限、单值和连续的球面函数都可以展开成面球谐函数的级数形式，因此地形球谐系数模型h(θ λ)可以表示为

(1)

式中，θ为余纬，范围为(0 π)；λ为球心经度；

为fnALF；

、

为球谐系数。

通过重新调整式(1)中的求和排序^[19-20]，在纬度方向上引入A_m(θ)、B_m(θ)，根据球谐函数正交性可得离散积分形式为

(2)

(3)

式中，h(θ_i λ_j)为等间隔的离散数值模型；n_max为球谐模型的截断阶数；N_lat、N_lon分别为h(θ_i λ_j)纬度圈、经度圈数量；δ为克罗内克尔符号。通常情况下，所采用的离散模型格网(N_lat×N_lon=n_max×2n_max)纬度圈与经度圈均等间隔分布全球。文献[21]证明了在经度圈上的按照等间隔分布采样点，可以有效保证由积分式向求和公式转换过程中三角函数的正交性。

截断阶数为n_max时，纬度方向的关系式表达为矩阵形式^[21]

(4)

式中

(5)

按照最小二乘法计算式(4)中的球谐系数矩阵，即等式两边同时乘以Y^T，则有

(6)

式中，Y^TY可表达为勒让德函数之间的内积形式

，但由于离散勒让德函数的非正交性，Y^TY的非对角元素不为零，因此不能够构成对角矩阵。

然而，球谐系数的积分式计算是建立在严密的勒让德函数正交性之上，所以按照式(3)中离散化处理，必会降低计算精度。为了减少数值积分法中计算精度的损失，文献[11 21]分别提供了满足式(7)的两组纬度相关先验权重w_GL、w_DH，较好地保证了离散勒让德函数求和对正交性的逼近

(7)

式中，N_lat^GL、N_lat^DH分别为两种方法所需格网的纬度圈数量。

1.2 Gauss-Legendre积分法

文献[21]引入数值分析中的Gauss-Legendre方法，根据截断阶数n_max建立高斯格网(N_lat^GL×N_lon^GL=n_max 1×2n_max 1)，要求Gauss-Legendre积分法的格网纬度(余纬θ)满足勒让德函数零点的条件，即

，且按照式(8)对各格网纬度圈数据赋权。需要特别指出的是，式(8)分母中P′n_max 1(cosθ_i)为勒让德多项式对cosθ_i的一阶导数，并非对θ_i的一阶导数

(8)

因此，利用Gauss-Legendre积分法分步计算球谐系数的公式可以表达为式(9)、式(10)

(9)

(10)

1.3 Driscoll/Healy积分法

针对离散勒让德非正交性问题，Driscoll/Healy积分法依据截断阶数n_max建立等间距格网(N_lat^DH×N_lon^DH=2n_max 2×2n_max 2)，按照式(11)为各纬度圈格网点赋权^[22]

(11)

因此，利用Driscoll/Healy积分法分步计算球谐系数的公式可以表达为式(12)、式(13)

(12)

(13)

1.4 完全正常化缔和勒让德函数

改进的Belikov方法计算公式为

(14)

式中

(15)

取θ=[0° 90°]，间隔为1°，分别由Belikov方法和改进的Belikov方法计算

至10 800阶次，利用式(16)进行精度统计，取以10为底的对数后精度情况如图 1、图 2所示

(16)

图 1 Belikov法计算精度

Fig. 1 Calculation accuracy of Belikov method

图选项

图 2 改进的Belikov法计算精度

Fig. 2 Calculation accuracy of improved Belikov method

图选项

由图 1、图 2可知，Belikov法在θ=7°时出现了溢出现象，而改进的Belikov法的计算精度在全球纬度区域均优于10^-12，且可以稳定递推至10 800阶，未出现溢出现象。

2 数值积分法构建地形球谐系数模型

2.1 数据准备

模型Earth2014_TBI2014.shc为科廷大学发布的10 800阶地形球谐系数系列模型之一^[17]，为避免输入数据的精度不一致对不同数值积分法还原球谐系数效果的影响，试验采用全球地形球谐系数模型Earth2014_TBI2014.shc(Topography-Bedrock-Ice layer)作为唯一的输入源，球谐综合构建各方法所需格网。

2.2 模型构建

试验利用Gauss-Legendre积分法和Driscoll/Healy积分法构建10 800阶地形球谐系数模型，为了更好地对比两种方法求解模型的精度，试验引入矩形离散积分法，即根据截断阶数n_max建立等间距格网(N_lat^Rect×N_lon^Rect=n_max×2n_max)，在考虑球面格网面积的情况下利用格网中点值代替格网值累加计算。模型构建具体步骤如下：

(1) 构建截断阶数n_max为10 800阶时的矩形离散积分等距格网、Gauss-Legendre格网和Driscoll/Healy等距格网，分别记3种方法的格网为Rect-grid、GL-grid和DH-grid。

(2) 以模型Earth2014_TBI2014.shc(10 800阶) 的地形球谐系数模型为输入，利用球谐综合计算完成各格网数据采样。

(3) 使用传统矩形离散积分法、Gauss-Legendre积分法和Driscoll/Healy积分法，构建全球10 800阶地形球谐系数模型sph.Rect10 800、sph.GL10 800和sph.DH10 800。

(4) 依据各10 800阶地形球谐系数模型恢复全球地形数值模型SHS-Rect10 800、SHS-GL10 800和SHS-DH10 800。

2.3 计算策略

通过串行累加计算得出球谐系数的方式在低阶次尚可完成任务，但随着计算阶次的升高，若不采用任何优化方法，计算耗时将难以忍受，而且普通计算机无法满足计算内存需求。因此试验根据分步法中各步骤函数计算特点制定计算优化方案，以达到减少计算机内存需要、提高计算速率的效果。

2.3.1 快速傅立叶变换技术

对于一组均匀采样序列x(n)，n=0 1 … N－1，其傅立叶变换满足式(17)

(17)

式中，X(k)为第k个FFT序列；

。

在分步法实现球谐分析构建球谐系数模型的第一个步骤中，即计算式(2)过程量Am(θ)、Bm(θ)时，各格网纬度圈上的点符合均匀采样的要求，本文按照式(18)使用快速傅立叶变换技术实现各纬度圈上的快速求解^[23]

(18)

式中 Re、Im分别代表变换所得实部与虚部。

2.3.2 基于OpenMP的多核并行技术

随着待求球谐系数阶数的增加，利用串行算法计算第二步骤中球谐系数

所需的时间呈指数级增长，且普通计算机硬件内存难以满足计算需要。OpenMP(open multi-processing)是一种多线程并行方法，可用于共享程序编程，旨在不增加硬件内存的前提下充分调度CPU现有线程以利用计算资源。因此在试验根据过程量A_m(θ)、B_m(θ)计算球谐系数

处引入基于OpenMP的多核并行技术，可以有效解决计算耗时过长和内存不足问题。

考虑式(3)中各阶次待求球谐系数仅仅与变量余纬θ相关，因而试验中以余纬θ为单元设计并行方案。这样的设计使得一个格网纬度圈上的所有格网点只需计算一次勒让德函数，将勒让德函数计算次数降到最低，减少了计算冗余。

在以θ为单元的并行过程中，前后计算相互关联，存在各并行线程同时访问读写同一变量或数组的现象，势必会因为内存位置读写冲突而导致程序中断。为保证各并行线程的畅通运行，试验利用pravite(·)子句将简单变量声明为各线程私有变量，将勒让德函数

、球谐系数等原有二维数组升维成三维数组

、

、

，在循环中通过对v的控制，实现各线程中对这些数组的独立访问与运算^[24]。

通过以上数组升维方法解决了多线程并行运算中内存位置读写冲突问题，但也造成了待开辟数组所需内存的急剧增加，例如在使用Gauss-Legendre格网数据(2161×4321)计算2161阶球谐系数模型时，2161个格网纬度圈意味着需要开辟2161个(2161×2160×2160)的三维数组，而一个双精度型数据占用8个字节，则仅勒让德函数就需要约75.12 GB内存，普通计算机内存难以达到要求，更不可能满足计算10 800阶球谐系数的需求。为解决这一问题，试验中将所有并行单元(N_lat) 平均分组计算，组数即为开辟线程数(num_threads)，若不能整除，余数为q，则将这q个单元分配至前q个线程参与计算。仍以Gauss-Legendre法求解2161阶球谐系数模型为例，假设开辟线程数为40时，则第1个线程分配到55个单元，其余39个线程分配到54个单元，此时开辟勒让德函数三维数组所需要内存仅为1.39 GB，适应了普通计算机的硬件水平。

基于上述讨论，试验提出联合快速傅立叶变换和基于OpenMP的多核并行技术的计算策略，流程图如下。为检验其计算效率，以Gauss-Legendre积分法为例，分别统计了使用串行累计方法与本试验策略计算不同截断阶数球谐系数的时间(表 1)。在表 2所示的计算环境条件下，计算360阶球谐系数加速比达到了279倍，较好地解决了串行编译耗时多、效率低的问题。快速傅立叶变换技术实现如下：

表 1 Gauss-Legendre积分法计算球谐系数时间

Tab. 1 Gauss Legendre quadrature for calculating the time of spherical harmonic coefficients

s

表选项

表 2 计算环境参数

Tab. 2 Computing environment Parameters

表选项

2.4 试验结果与分析

2.4.1 谱域分析

阶误差描述了各阶位场的功率谱分量的差异，可用于评估恢复模型精度，阶误差RMS为^[25]

(19)

式中，

为参考模型球谐系数，本文试验中为模型Earth2014_TBI2014.shc系数；

为待评估模型球谐系数。阶误差可以反映待评估模型与参考模型之间的差异，差值越小，说明两个模型符合程度越高。

图 3为各方法解算模型的阶误差RMS对比情况，为直观反映各阶误差量级情况，图中还绘制了模型Earth2014_TBI2014.shc的阶方差曲线作为参考。

图 3 模型球谐系数阶误差(10 800阶)

Fig. 3 Model spherical harmonic coefficient order error(10 800 degree)

图选项

如图 3所示，试验所构建的全球10 800阶地形球谐系数模型sph.Rect10 800、sph.GL10 800和sph.DH10 800的阶误差在各阶段都保持相对稳定，其阶误差水平分别在10^-5、10^-11和10^-13左右。由此看出，3种方法均可以有效恢复球谐系数模型，但恢复效果存在较大差异，相较于传统矩形离散积分法，Gauss-Legendre积分法和Driscoll/Healy积分法明显更具有优势，可在将阶误差控制在10^-10量级以内。

为了进一步分析模型解算精度，给出了各模型球谐系数误差谱，如图 4所示。由图 4可知，相较于各自模型的高次项，3个模型的各阶低次项处的误差均有不同程度增大。其中，模型sph.Rect10 800的部分高阶低次项误差在10^-4量级以上，模型sph.GL10 800全阶低次项在10^-10~10^-12量级之间，模型sph.DH10 800全阶低次项在10^-13量级左右，波动较小。整体而言，模型sph.DH10 800的误差谱更为平滑，误差更小，与初始模型更为接近。

图 4 模型球谐系数误差谱

Fig. 4 Model spherical harmonic coefficient error spectrum

图选项

2.4.2 空域分析

分别将全球地形数值模型SHS-Rect10 800、SHS-GL10 800和SHS-DH10 800与原始模型Earth2014_TBI2014.shc所恢复的数值模型作差比较，记为DIF-SHS-Rect、DIF-SHS-GL、DIF-SHS-DH，各数值模型差值统计信息见表 3。

表 3 数值模型差值统计

Tab. 3 Numerical model difference statistics

m

表选项

由表 3可知，差值模型DIF-SHS-Rect、DIF-SHS-GL、DIF-SHS-DH的标准差分别为11.801 3、1.326×10^-6和6.221×10^-9m，这表明所对应地形球谐系数模型分别能够以10、10^-6、10^-9m精度恢复出地形原貌。

统计不同纬度圈各模型差值均方根，其量级随纬度变化的情况如图 5所示，可以发现，尽管差值模型DIF-SHS-Rect绝大多数点可以维持在接近10^-2m的量级，但其在高纬度区域的均方根达到百米以上量级，与中低纬度精度不相匹配，下降约5个量级，是导致该模型精度较低的重要原因。由图 5中的DIF-SHS-GL曲线可以看出，虽然地形模型SHS-GL10 800能达到优于10^-4m的精度，但也存在高纬度地区精度相对较差的现象，下降约3个量级。而DIF-SHS-DH曲线较为平稳，在全球各纬度均能保持在10^-8~10^-10m区间，未出现高纬度区域精度衰退的现象。归纳总结3种数值积分方法的构建模型要求及效果情况见表 4。

图 5 模型差值均方根随纬度变化

Fig. 5 Model RMS variation with latitude

图选项

表 4 各方法建模要求及效果

Tab. 4 Modeling requirements and effects of each method

表选项

由表 4可知，高斯格网点数与矩形离散等距格网相当，仅多出约一个纬圈的数据量，但因Gauss-Legendre权重的干预，改善了离散勒让德函数的非正交性，使得构建效果显著提升。相比之下，Driscoll/Healy积分法既增加了约1倍格网点数，又考虑了非正交性问题进而赋予了合适的权重，它构建的模型精度最高。

比较矩形离散积分法和Driscoll/Healy积分法，两者最大的相同之处为均采用了等距格网，而两者的不同点在于Driscoll/Healy积分法格网点数约为矩形离散积分法的两倍，以及所赋予各纬度的权重。但正因为这些不同点，使得两种方法所求算模型有着极大的差距。为了进一步比较两种方法的特性，利用Earth2014的前360阶模型系数，构建了矩形离散等距格网(360×720)和Driscoll/Healy等距格网(722×722)，分别用于矩形离散积分法计算360阶球谐系数模型，记为sph.Rect360、sph.Rect360.DHgrid；将Driscoll/Healy积分法恢复出的Driscoll/Healy等距格网(722×722)720阶球谐表达记为sph.DH720。绘制3个模型的阶误差如图 6所示。可以发现sph.Rect360、sph.Rect360.DHgrid的阶误差相差大约1个量级，分别为10^-2和10^-3；sph.DH720的前360阶阶误差保持在10^-14量级水平，360阶之后的阶误差量级陡然升至10⁰左右，不再有效。因此可以得出以下结论：采用更密集格网数据可以改善矩形离散法的数值积分精度，但并不能达到Driscoll/Healy积分法的效果；同时，尽管Driscoll/Healy积分法格网具有两倍数量的采样点，但仅仅是为了符合Driscoll/Healy积分法要求而采取的物理分辨率上的加密，它并不含有360阶以上更高频率的频谱信息，因此所还原的更高阶次球谐系数并不能保持原有的精度。

图 6 模型球谐系数阶误差(360阶)

Fig. 6 Model spherical harmonic coefficient order error(360 degree)

图选项

3 全球10 800阶地形球谐系数模型构建

3.1 数据准备

Earth2014_TBI为地形球谐系数模型Earth2014_TBI2014.shc的1′×1′数值地形模型，它采用的地形数据为全球陆地地形、基岩、冰盖数据，即可以看作没有水的地球地形，模型陆地地形数据由SRTM V4.1、SRTM_30 PLUS提供，海洋地形数据来源于SRTM30_PLUS V9，南极洲和格陵兰岛的冰盖和基岩地形数据分别由模型Bedmap2和GBT V3提供。

STO_IEU2020为全球1′×1′海底地形模型，范围为(180°W—180°N，60°S—64°N)，平均相对精度优于6%。模型利用高精度全球重力数据SIO V29.1，并融合了最新的单波束水深测量数据、多波束水深测量数据、雷达传感器探测水深数据等多源水深测量结果构建而成^[26]。相较于2014年发布的模型Earth2014_TBI，模型STO_IEU2020所采用的构建数据更具时效性，因此试验利用模型STO_IEU2020数据取代模型Earth2014_TBI中相应海域海底地形数据，构建全球地形数值模型。由于两个模型的垂直基准均为平均海水面，所用数据均基于WGS-84椭球，且海底地形平均精度仍在百米量级，远远大于两个模型基准即平均海水面之间的差异，因此进行数据融合拼接时可忽略其影响。

3.2 模型系数求解

对比3种数值积分方法可知，Driscoll/Healy积分法构建球谐系数模型精度最高。因此试验根据所构建全球地形数值模型，利用双三次插值法采样获得Driscoll/Healy等距格网(21 602×21 602)，并采用联合快速傅立叶变换和基于OpenMP的多核并行技术的计算策略快速计算模型系数，将所构模型记为sph.10 800_IEU。综合以上模型构建策略，整理模型Earth2014_TBI2014.shc与sph.10 800_IEU的构建基本信息见表 5。

表 5 Earth2014_TBI2014.shc与sph.10 800_IEU信息对比

Tab. 5 Comparison of earth2014_TBI2014. SHC and sph.10 800_IEU

表选项

3.3 模型精度评估

如图 7所示，模型sph.10 800_IEU的阶方差与模型Earth2014_TBI2014.shc高度吻合，阶误差曲线平滑，说明试验构建模型的可靠性。

图 7 模型sph.10 800_IEU阶误差RMS

Fig. 7 Error degree RMS of model sph.10 800_IEU

图选项

为评估模型sph.10 800_IEU精度，分别选取全球4片海域的实测水深数据进行外部检核，记为海域1、海域2、海域3、海域4，检核点数量分别为2485、9065、6967、6459个，海域检核点分布如图 8所示。模型sph.10 800_IEU、模型Earth2014_ TBI2014.shc球谐综合计算与检核点海深差值在各海域的统计结果见表 6。其中，“一次检核”代表海域所有检核点统计结果；“二次检核”代表经过3σ准则剔除后剩余检核点统计结果；“数据保留率”为剩余检核点数量的占比；“相对精度”为差值标准差与各海域平均海深的比值。

图 8 海域检核点分布

Fig. 8 Distribution of inspection points in each sea area

图选项

表 6 球谐模型检核结果

Tab. 6 Spherical harmonic model check result statistics

表选项

由表 6可知，在试验选取的各个海域中，二次检核的数据保留率均达到97%以上，标准差均在百米量级，说明模型计算结果真实可信；模型sph.10 800_IEU的球谐综合结果略优于模型Earth2014_TBI2014.shc，两者在海域2和海域3的检核结果标准差差异分别为5.43、3.12 m，在其余两个海域检核结果标准差差异在1 m左右，模型sph.10 800_IEU在4片海域相对精度分别提升了0.02%、0.17%、0.09%、0.07%。

综合以上分析，全球10 800阶球谐系数模型sph.10 800_IEU精度与模型Earth2014_TBI2014.shc相当，试验建模结果质量较高。

4 结论

本文实现了基于矩形离散积分法、Gauss-Legendre积分法和Driscoll/Healy积分法的球谐系数求解，采用模型Earth2014_TBI2014.shc球谐综合数据解算了10 800阶全球地形球谐系数模型sph.Rect10 800、sph.GL10 800和sph.DH10 800，通过制定计算策略提高了计算效率，并从谱域和空域分别比较分析了各种数值积分方法的模型构建精度，最终结合最新海底地形模型STO_IEU2020格网数据，构建了10 800阶全球地形球谐系数模型sph.10 800_IEU，得到如下结论：

(1) 分析所构建模型的谱域信息，传统矩形离散积分方法构建模型阶误差在10^-5量级，而Gauss-Legendre积分方法和Driscoll/Healy积分方法构建模型阶误差可分别达到10^-11和10^-13量级，高精度还原了原模型的频谱信息，但3种方法构建模型均存在各阶低次项处误差偏大的现象。

(2) 分析所构建模型的空域信息，矩形离散积分方法构建模型精度为10 m，而Gauss-Legendre积分方法和Driscoll/Healy积分方法具有10^-6、10^-9m精度恢复出地形原貌的能力；矩形离散积分方法和Gauss-Legendre积分方法描绘高纬度区域精度较差，相比中低纬度各有约5个、3个量级的下降，而Driscoll/Healy积分方法可以稳定且精细表达全球各个纬度地形信息。

(3) 改进的Belikov法可以在全球纬度范围内将fnALF高效稳定递推至10 800阶，且精度均优于10^-12。

(4) 通过快速傅立叶变换技术优化经度方向计算，以及在并行条件下采用数组升维和分组计算技术处理纬度方向计算，可以充分发挥现有计算设备计算能力，提高球谐系数模型构建效率，总体加速比可达到279倍。

(5) 利用Driscoll/Healy积分方法，联合Earth2014_TBI与最新海底地形模型STO_IEU2020数据所构建的全球10 800阶球谐系数模型sph.10 800_IEU在试验海域的相对精度略高于模型Earth2014_TBI2014.shc，两者总体精度相当，说明模型sph.10 800_IEU质量可达到行业前列水平。

作者简介

第一作者简介：单建晨(1997—) 男 硕士 研究方向为物理大地测量。E-mail: 1170633968@qq.com

初审：张艳玲

复审：宋启凡

终审：金 君

资讯