![]()
本文内容来源于《测绘学报》2024年第9期（审图号GS京(2024)1896号）
缔合Legendre函数的快速插值计算及其在六边形网格模型重力异常计算中的应用李新星![]()
,1,2,3, 范昊鹏![]()
,1, 万宏发1, 范雕1, 冯进凯1
1.信息工程大学地理空间信息学院，河南 郑州 450001
2.地理信息工程国家重点实验室，陕西 西安 710054
3.智慧地球重点实验室，北京 100020
摘要
鉴于局部区域非等纬度分布点的超高阶次球谐综合计算效率较低，本文深入研究了缔合Legendre函数的插值算法，结合模型重力异常求解的展开式，提出了超高阶模型空间重力异常插值快速计算方法，为了验证本文方法在东西狭长分布的点集更具有应用优势，采用球谐旋转变换技术进一步提升了计算效率。试验结果表明，利用2160阶次EGM2008模型计算日本地区同一高度的30 303个非等纬度分布的六边形网格点处模型重力异常，插值快速计算方法相比逐点计算，在误差不超±0.005 mGal水平下，计算耗时从3 669.41 s缩减到98.05 s，同时经过球谐旋转变换，将南北狭长的日本区域点集旋转为东西狭长分布，使得上述相应计算内容的耗时进一步由98.05 s缩减到19.06 s，计算效率相比最初方法提速比达到近200倍，有效解决了非等纬度分布点的超高阶模型重力异常快速解算的效率难题，且验证了该方法在东西狭长分布的情况下具有更高的提速比。
关键词
地球重力场模型; 球谐综合; 勒让德函数; 六边形网格; 球谐旋转变换; 插值方法国家自然科学基金(U23A2028)(41404020)(42174007)(42074013); 地理信息工程国家重点实验室基金(SKLGIE2023-Z-1-1); 中国地质大学（武汉）中央高校基本科研业务费(2024XLB13)
第一作者：李新星（1988—），男，博士，讲师，主要研究方向为物理大地测量。E-mail：minibad@126.com
通讯作者: 范昊鹏 E-mail：fanhaopeng2008@163.com
本文引用格式
李新星, 范昊鹏, 万宏发, 范雕, 冯进凯. 缔合Legendre函数的快速插值计算及其在六边形网格模型重力异常计算中的应用[J]. 测绘学报, 2024, 53(9): 1737-1747 doi:10.11947/j.AGCS.2024.20230110
LI Xinxing, FAN Haopeng, WAN Hongfa, FAN Diao, FENG Jinkai. A fast method for interpolation of the associated Legendre functions and its application to the calculation of local hexagonal grid point gravity anomalies using an ultra-high-degree gravity field model[J]. Acta Geodaetica et Cartographica Sinica, 2024, 53(9): 1737-1747 doi:10.11947/j.AGCS.2024.20230110
地球重力场模型可用于计算地球及空间任意一点的重力异常、垂线偏差、高程异常等重力场元，伴随着观测资料和计算资源的极大丰富，模型截断阶次已经由EGM96模型的360阶次[1]发展到EGM2008模型的2190阶2159次[2]，直至最新模型XGM2019e的5540阶次[3]。模型阶次的升高在提升重力场分辨率和精度的同时，也对使用者的计算机性能要求越来越高，考虑到超高阶模型重力场元的计算在利用“移去-恢复”方法构建厘米级大地水准面[4-5]、重力数据的精化处理、利用更高阶次的地形模型系数求解地形引力效应值等方面具有重要应用[6]，其快速高精度计算问题一直受国内外学者关注。球谐综合即利用某场元的球谐级数展开式和相应的球谐系数求解空间该场元值的过程，对于利用N阶球谐模型计算全球相应分辨率下2N2个场元值的球谐综合过程，计算复杂度为O(N4)，其中高阶次的缔合勒让德函数(associated Legendre functions, ALFs)的计算最为耗时。目前针对ALFs高效快速计算，以递推方法为主，包括行推、列推、跨阶次递推和Belikov递推等[7-10]，虽然各种递推方法的计算效率存在差异，但差异相对较小[11]，因此尽可能地避免计算ALFs成为提升球谐综合过程效率的重要途径。如，文献[12-13]求和方法不需要计算全部ALFs的值，只需要计算其中扇谐项的值即可完成模型值求解，该方法有提速效果但是并不明显[14]；一维、二维快速傅里叶变换(1D-FFT、2D-FFT)方法利用ALFs与三角函数之间的关系规避ALFs的递推，实现了球谐综合效率的飞速提升，计算复杂度由O(N4)分别减小至O(N3log2N)和O(N2log2N)，但该系列方法通常用于一定宽度的纬度带或全球范围的计算[15-16]；文献[17]通过交换球谐综合过程中阶次的求和顺序减小ALFs的重复计算，达到提升解算效率的目的，但仅限于按地理网格等纬度分布且高度一致的区域规则点；文献[18]提出非全次ALFs计算方法，通过省去无效阶次ALFs的递推过程来增加球谐综合的解算效率，但该方法仅适用于高纬度区域。为将上述各类快速方法扩展应用至计算非等纬度分布点、非等高度分布点的情况，非等间距傅里叶变换(nonequispaced FFT, NFFT)和泰勒级数展开等手段被使用，分别解决了全球随机离散点、局部非等高度离散点下的计算问题[19-23]。然而，对于局部区域大量的非等纬度分布点来说，如何提升其球谐综合的效率仍值得进一步研究[24-25]。文献[26]给出了ALFs的埃米尔特插值算法，经文献[27]验证，可引入该算法解决局部区域非等纬度分布点的计算问题。本文深入研究ALFs的插值方法，通过与球谐综合求解过程融合，实现局部区域大量非等纬度分布点的超高阶模型重力异常的快速高精度计算，效率提升显著。1 缔合Legendre函数的插值计算
1.1 缔合Legendre递推求解
![]()
(1)
![]()
(2)
式中，n和m分别对应ALFs的阶和次，且满足n≥m；球坐标系下x通常用cosθ代替，θ表示计算点位处的地心余纬。
ALFs可使用超几何函数、递推等多种方法求解，目前，递推方法依然是更加高效的选择，如行推、列推、跨阶次递推和Belikov递推等[11]。X-数法[28]、正则化的Belikov方法[8]等均是针对高阶次函数值溢出的问题，对上述4种基本递推方法的改进，但在解决溢出问题的同时，必然增加了计算耗时。鉴于跨阶次递推在20 000阶以内的超稳定特性[10]，这里只简述ALFs的跨阶次递推方法，递推公式为![]()
(3)
![]()
(4)
![]()
(5)
当m>2时，递推系数αnm、βnm、γnm均小于1，所以跨阶次递推方法是稳定的，递推过程如图1所示。图1
![]()
图1 ALFs的跨阶次递推
Fig.1 The cross-degree-order recursion sequences to compute ALFs
其他3种递推方法与跨界次方法类似，且在计算效率方面差异不大，具体后文将进行论述。1.2 缔合Legendre函数插值方法
ALFs的计算耗时较大，导致超高阶次球谐综合过程的效率低下。为了加快ALFs的计算效率，可以考虑利用已经计算好的![]()
(cosθ0)值来“间接”快速计算附近余纬θ处的ALFs值。根据式(1)定义，![]()
(cosθ)在余纬方向的泰勒级数展开公式为![]()
(6)
因此，当q0处的![]()
、![]()
、![]()
已知，且(θ-θ0)足够小的时候，可利用式(6)快速近似计算![]()
。然而经过分析，对于2160阶次，纬差(θ-θ0)为1'且采用式(6)中泰勒二阶展开进行计算时，![]()
的相对误差可达0.01%，因此对于高精度的解算来说是完全不够的。根据定义式(1)、式(2)可知，Pn(x)与![]()
为x的最高次幂不超过n的多项式，且![]()
在全球范围内有除了南北两个极点处函数值为0外，还有n-m个零点纬度，即满足![]()
的余纬值有n-m个；因为m≥0，零点纬度最多为n个，所以在Δθ≈π/n范围内，![]()
的曲线通常不会包含超过2个极值点(零点纬度在极区附近分布更密集，可能会存在超过3个极值点的情况)。因此在Δθ空间内，可使用一个关于余纬θ的多项式f(θ)来拟合![]()
在Δθ区间内的曲线，如图2所示。图2
![]()
图2 ALFs的内插
Fig.2 The interpolation process for calculating ALFs
因此，当计算任意余纬θ处的ALFs值时，可先行利用递推方法计算θi、θi+1、θi+2 3个等间隔分布余纬处ALFs及其导数的值并拟合出多项式f(θ)，然后利用f(θ)内插求出该范围中间任意点θ处的![]()
值。步长Δθ最大不能超过由奈奎斯特(Nyquist)采样定律确定的间隔π/N，其中N是所求![]()
中阶n的最大值，或者计算重力场元所使用的重力场模型的最高阶次。由图2可知，n4>n1,n-m越小即被拟合的函数![]()
越光滑、步长Δθ越小、多项式f(θ)的最高幂次越高，拟合的精度越好。为保证一定的插值精度，经反复试验，选择θi、θi+1、θi+2 3个连续等间隔分布余纬处的![]()
、![]()
、![]()
共9个值作为已知值，拟合一个8次多项式f(θ)，表达式为![]()
(7)
根据多项式f(θ)在θi、θi+1、θi+2 3个余纬处的值及其一阶、二阶导数值为递推方法求得的![]()
、![]()
、![]()
共9组已知值，以此作为约束，并将相应值进行如下简洁表示![]()
(8)
![]()
(9)
将(9)式回代入式(7)，得到![]()
在区间θi,θi+2内的近似表达![]()
(10)
式中，θi≤θ≤θi+2;b1—b9是只与所求点余纬有关而与n、m无关的系数，其表达式为![]()
(11)
其中
![]()
(12)
上述方法就是ALFs的插值快速计算方法，该方法对于余纬区间[θi,θi+2]内分布的大量任意点，只需要采用常规递推方法计算θi、θi+1、θi+2 3个余纬处的ALFs，及其一阶、二阶导数值，然后利用上述插值方法即可快速地求解其余大量点的ALFs值，而无须再进行效率低下的逐点递推，从而达到提升效率的目的。上述求解过程中ALFs的一阶和二阶导数依然采用相应的递推公式进行求解，具体公式此处从略，详见文献[29]。1.3 试验分析
为对比ALFs递推方法和插值方法的解算效率，采用VS2019(C++语言)在处理器为Intel(R)Core(TM)i7-9750H CPU@2.60 GHz，内存为32 Gb的便携式计算机上，使用传统4种递推方法和插值方法分别计算余纬[45°,45°10']范围内间隔5″共121个纬度处的2160阶次![]()
和余纬[45°,45°4']范围内间隔2″共121个纬度处的5400阶次![]()
的值，计算耗时统计见表1。表1 计算用时比较
Tab.1 The comparison of run-time
| 最高阶次 | 行推方法 | 列推方法 | 跨阶次递推 | Belikov递推 | 插值方法 |
|---|
| 2160 | 23.8 | 17.5 | 18.5 | 6.6 | 3.5 |
| 5400 | 135.3 | 115.0 | 115.2 | 68.7 | 21.6 |
新窗口打开| 下载CSV
由表1可知，利用Nyquist准则确定间距的插值方法计算的ALFs值，其计算效率相比跨阶次方法提升了近5倍，相比当前最快的Belikov方法也提升了近2倍的效率，且随着阶次越高，提升效果越明显。至于计算精度，将从后文试验中予以验证，此处不再赘述。至此，采用ALFs快速插值方法计算任意纬度处的模型重力场元点值，效率提升2~5倍，但对于上万点的计算来说，依然还是一件非常耗时的工程。基于此，本文将继续研究ALFs的快速插值方法如何在超高阶模型重力场元计算中发挥更好的提速作用。2 超高阶模型重力异常快速计算
2.1 模型重力异常计算公式
利用球谐系数求解空间一点(r,θ,λ)处的模型重力异常严密公式为[30]![]()
(13)
式中，r、θ、λ表示空间计算点的地心向径、地心余纬和经度；GM表示地心引力常数；R表示所使用的球谐系数对应的参考球半径，在利用EGM2008模型计算时该值取为6 378 136.3 m。
在实际工程应用中，由于ALFs的解算速度太慢耗时太长，并不会直接采用式(13)逐点计算模型重力异常，通常忽略高程的影响，直接采用球近似下的重力异常计算公式![]()
(14)
显然，式(14)在计算地心纬度相等的点集模型重力异常时，只需要计算1组相同的gE和gF，即可计算该纬度上所有点的重力异常值，从而提升计算效率。但对于地形起伏较大的区域，忽略地形的影响会造成较大的计算误差，针对该问题，文献[21-23]利用径向方向的泰勒级数展开方法予以解决。因此，本文后续将着重研究等高程面上，非等纬度分布点的模型重力异常快速计算问题。2.2 模型重力异常计算的插值方法
在计算纬度相等的点集的模型重力异常时，按照式(14)只需要计算一组该纬度的gE和gF，即可应用于所有点来求解模型重力异常，但对于局部区域大量纬度不相等的点，尚没有很好的高效计算方法。缔合Legendre函数的插值方法给了该问题高效求解的新思路，确定ALFs计算节点的示意如图3所示。图3
![]()
图3 确定ALFs计算节点
Fig.3 Determining the nodes for calculating ALFs
计算图3中随机分布点的模型重力异常，首先需确定该系列待计算点的余纬区间[θN,θS]，所有的计算点均位于该区间内，并根据计算模型重力异常所采用的模型阶次N和精度要求，确定该范围内ALFs递推节点的间隔Δθ(Δθ≤π/N)及对应节点。图3中确定出了覆盖区间[θN,θS]的最小节点数(θ1,θ2,…,θn)。然后基于跨阶次递推方法求解每一个节点纬度处的ALFs及其一阶、二阶导数值，那么所有待计算点处的ALFs值均可由上下相邻的3个节点的ALFs及其导数值内插求得。将式(10)代入式(14)，可得![]()
(15)
其中
![]()
(16)
因此，只需要递推求解θi、θi+1、θi+2 3个余纬处的gEm、gFm、gE'm、gF'm、gE″m及gF″m共18组参数，则可利用式(15)求解区间[θi,θi+2]内所有点的模型重力异常值，其中插值系数b1—b9由所计算点余纬θ与节点余纬θi的差异和步长Δθ决定，详见式(11)和式(12)。2.3 试验分析
为验证ALFs的插值方法在超高阶模型重力异常快速计算中的应用效果，本文在2.2节的计算环境下，利用全球均匀分布的等面积ISEA4H六边形结构[31-32]，如图4所示，使用Level=11层次剖分，生成日本地区30 303个非等纬度分布的等面积六边形网格点，单个网格平均实地面积约12.16 km2，网格平均间距约3.77 km。选择日本地区是因为其南北跨度相对较大，可充分验证算法的优势，也为下一节的改进效果提供参考。图4
![]()
图4 ISEA4H网格分布
注：L=4, Num=2562。
Fig.4 Voronoi diagram of ISEA4H grid
利用2160阶次EGM2008重力场模型系数，基于传统方法的式(14)计算该区域内30 303个点的模型重力异常值，其计算结果和耗时统计如表2和图5所示。表2 常规方法计算六边形网格重力异常点值结果统计
Tab.2 Results for calculating point gravity anomaly values on hexagonal grid using conventional method
| 重力异常点数 | CPU耗时/s | 重力异常值/mGal |
|---|
| 最大值 | 最小值 | 平均值 | 标准差 |
|---|
| 30 303 | 3 669.41 | 177.48 | －52.84 | 62.41 | 35.81 |
新窗口打开| 下载CSV
图5
![]()
图5 日本范围内ISEA4H六边形网格重力异常点值
Fig.5 Point gravity anomaly values on ISEA4H hexagonal grid within Japan
为分析不同节点的步长Δθ对计算结果精度以及计算效率的影响，分别采用5'、2'30″、2'和1'的步长利用插值快速方法，使用式(15)计算相同点位处的2160阶次模型重力异常，结果见表3。表3 使用ALFs插值算法计算的模型空间重力异常与传统方法之间的差异
Tab.3 The difference between the computed point free-air gravity anomalies using the ALFs interpolation algorithm and those by a conventional method
| 步长Δθ | 节点数 | 耗时/s | 模型重力异常计算误差/mGal |
|---|
| 最大值 | 最小值 | 平均值 | 标准差 |
|---|
| 10′ | 117 | 47.94 | 7.71×10－1 | －5.42×10－1 | 1.65×10－3 | 7.69×10－2 |
| 5′ | 233 | 98.05 | 3.34×10－3 | －2.90×10－3 | －3.06×10－6 | 3.14×10－4 |
| 2.5′ | 465 | 180.85 | 8.45×10－6 | －7.32×10－6 | 9.87×10－10 | 7.41×10－7 |
| 1′ | 1163 | 439.56 | 2.20×10－9 | －1.88×10－9 | 6.84×10－13 | 2.03×10－10 |
| 30″ | 2327 | 884.37 | 1.27×10－10 | －1.32×10－10 | 1.79×10－14 | 9.28×10－12 |
新窗口打开| 下载CSV
由表3可知，对于2160阶次的模型重力异常求解，节点间的步长Δθ应小于由Nyquist采样定律确定的5'，否则插值将造成较大的近似误差。如，当步长取10'的时候，最大模型重力异常误差为0.77 mGal，不满足严密计算需求，而5'步长的最大误差只有不足0.005 mGal，能够满足当前绝大多数应用需求。另一方面，步长减小必然造成节点数量的增多，从而增加了节点处ALFs值及其一、二阶导数的递推计算耗时；加之节点间内插点数相应变少，解算效率提升越发不明显，甚至更糟；当步长小到节点间只有1~2个内插点时，插值方法优势尽失，且耗时将增至传统逐点计算方法的3倍甚至更高，这是因为传统方法不需要计算ALFs的一、二阶导数值。因此，本文快速插值方法使用的原则为在满足精度要求的情况下，尽量减小节点的数量。此外，研究过程发现，本文方法在东西狭长点集分布区域具有更大的改善效果。下文采用旋转变换的方法对其说明。3 球谐旋转变换下的快速计算
3.1 球谐系数的旋转变换
球谐旋转变量变换(spherical harmonic rotation, SHR)可将南北向分布纬度跨度较大的区域和球谐系数同时进行坐标系旋转[18,33]，使得旋转后的区域呈“东西向”分布，减小纬度跨度；从而利用前述快速方法时，可在不减小计算点数量和密度的情况下，减小节点数量，进一步提升计算效率。已知坐标系O-xyz下的球谐系数为![]()
和![]()
，坐标系沿着z轴旋转α角，记为Rz(α)；则在新的坐标系O-x'y'z'下，可得到![]()
和![]()
旋转变换后的球谐系数![]()
和![]()
。该组系数结合旋转后的新的坐标，可求得与原坐标系中对应点一致的调和值，即保持地球空间的调和函数值不变。当坐标系沿z轴旋转α角时，相应的球谐系数的变换公式为![]()
(17)
当O-xyz沿着y轴旋转β角时，记为Ry(β)，得到新的坐标系O-x″y″z″，相应的球谐系数的变换公式为![]()
(18)
式中，变量![]()
的递推求解公式从略，详见文献[18]。总之，对于坐标系的任意旋转，球谐系数总可以进行相应的旋转变换，确保球谐模型计算同一点的场元值保持不变。
3.2 区域点集的旋转参数确定
对于确定的计算点集，可尝试通过坐标系旋转，使得该区域在纬度方向上的跨度最小，且纬度范围近似对称于纬度0°，以利用ALFs函数的对称特性节省一半的计算耗时。实现上述目的，首先，确定计算点集的凸包[34]，并采用旋转卡壳法获取凸包的最小卡壳距离和对应的方位角；然后，以最小卡壳距离的方位旋转为0°、最小卡壳宽度线Hmin的中点位于纬度0°为约束，确定按照z-y-z轴旋转的旋转角，进而将坐标和球谐系数分别进行旋转。旋转卡壳法是由文献[35]提出用于寻找凸多边形外接最长或最短宽度的算法。图6为旋转卡壳法，凸多边形对应的最小卡壳宽度一定是点线式的，因此只需考虑点线式的旋转卡壳宽度即可，对构成凸多边形的n条边进行循环，查找所有点与凸包各边之间的最远距离，再确定上述所有最远距离中最小的那一组，即最小卡壳宽度。平行线卡壳的对应点称为对踵点，所有的卡壳宽度线均起始于对踵点。图6
![]()
图6 旋转卡壳法
Fig.6 Rotating calipers method
卡壳宽度线的求解，利用叉乘计算对踵点与对应边的三角形面积，再计算点到直线的距离，图6中三角形ABC的面积可以表示为![]()
(19)
以BC边对踵点A为例，可以计算出向量AC和AB和边的长度![]()
，继而得到顶点A至直线BC的距离H，即为线BC和对踵点A对应的卡壳宽度，依次循环凸包上每一条边对应的对踵点及相应卡壳宽度，最后确定所有卡壳宽度中的最小值Hmin。旋转卡壳算法流程为
| 算法：旋转卡壳法 |
| 输入：凸多边形P={p1,…,pn} |
| 输出：最短宽度、最短宽度的方向 |
| initial width[n] |
| for i=1,2,…,n do |
|
| j=i+1 |
|
| while Area(pi,pi+1,pj)<Area(pi,pi+1,pj+1)do |
|
|
| j=j+1 |
|
|
| if j>n then |
|
|
|
| j=1 |
|
|
| end if |
|
| end while |
|
| width[i]=Width And Direction Compute(pi,pi+1,pj) |
| end for |
| return Width, Direction=Minimum (width) |
新窗口打开| 下载CSV
当求得最小卡壳宽度线的宽度值Hmin和方位角α后，为了使得计算区域纬度跨度足够小且对称于赤道，需要对原点位分布进行三维坐标旋转得到，而相应的旋转矩阵表达式为![]()
(20)
式中，θ、λ为该区域最小卡壳宽度线中点的地心余纬和经度；Rx、Ry、Rz分别表示绕x、y、z轴旋转的旋转矩阵，具体表达式本文不再详细给出。而任意的一个空间旋转又可以分解为只绕两个固定轴的旋转，例如z-y-z轴顺序的旋转，由旋转矩阵R元素求解对应轴旋转角φ的计算公式为
![]()
(21)
由式(21)即可求得，当区域内点集分别进行![]()
旋转，相应球谐系数也进行相应旋转变换。最终该区域的最小卡壳宽度线沿着子午线方向，使得所有点纬度跨度范围最小，且关于赤道对称，达到了旋转的目的。3.3 试验分析
将日本地区30 303个六边形网格点分别对坐标和球谐系数进行旋转，然后再按照插值快速计算方法计算该点集的模型重力异常值，并与不进行旋转的方法进行精度和效率的对比。试验流程如图7所示。图7
![]()
图7 球谐旋转变换的快速计算流程
Fig.7 Fast calculation process flowchart for spherical harmonic rotation transformation
根据原始点位分布，30 303个点集的纬度跨度达到19.40°。如图8(a)所示，利用卡壳方法，确定了最小卡壳宽度线的宽度值及其方位值，经过式(20)约束求解，确定了对30 303个点集沿z-y-z轴的旋转角为160.925 44°、58.172 21°和136.281 58°，旋转后使得点集分布沿纬度分布最小且纬度区间关于赤道对称，如图8(b)所示，纬度跨度只有6.54°，而对于节点的选取，利用ALFs的对称性，只需在纬度0°~3.27°之间选取。同时，利用EGM2008的2160阶次球谐系数，经过上述3个角的旋转，得到了新的球谐系数EGM2008R，该转换过程总耗时604.12 s。考虑到该过程与所计算点的数量无关，且可根据计算区域提前转换，故本文暂对此过程不作讨论。图8
![]()
图8 日本地区点集的旋转
Fig.8 Illustrating the rotation of point sets in the Japanese region
利用EGM2008R和旋转后的坐标，采用不同的节点步长计算2160阶次模型重力异常值，结果见表4。表4 坐标旋转后使用ALFs插值算法计算的模型空间重力异常与传统方法之间的差异
Tab.4 The difference between the computed point free-air gravity anomalies using the ALFs interpolation algorithm after coordinate rotation and those by a conventional method
| 步长Δθ | 节点数 | 耗时/s | 模型重力异常计算误差/mGal |
|---|
| 最大值 | 最小值 | 平均值 | 标准差 |
|---|
| 10′ | 21 | 11.13 | 3.57×10－1 | －1.62×10－1 | 1.96×10－3 | 1.88×10－2 |
| 5′ | 41 | 19.06 | 3.61×10－3 | －4.10×10－3 | －6.53×10－6 | 3.41×10－4 |
| 2.5′ | 81 | 35.62 | 1.01×10－5 | －1.11×10－5 | －9.77×10－10 | 8.20×10－7 |
| 1′ | 199 | 85.74 | 2.85×10－9 | －2.44×10－9 | 1.50×10－12 | 2.26×10－10 |
| 30″ | 395 | 167.13 | 1.33×10－10 | －1.40×10－10 | 8.62×10－13 | 1.08×10－11 |
新窗口打开| 下载CSV
由表4可知，日本地区南北狭长分布的点集经过坐标系旋转后，显著减小了纬度范围，从而减小插值快速计算的节点数量，解算效率进一步提高；且球谐系数的旋转对结果精度的影响非常小，可忽略不计。结果验证了本文所提快速计算方法在东西狭长分布的区域具有更显著的优势。4 结论
本文深入研究了ALFs的插值方法，并将其应用于球谐综合求解过程，显著提升了非等纬度分布点集的模型重力异常的计算效率，解决了局部区域大量非等纬度分布点球谐综合效率低下的难题，同时证明了该方法在东西向狭长地区的独特优势。本文以EGM2008模型计算日本地区30 303个球面非等纬度六边形网格点2160阶次模型重力异常值为算例，得到如下结论。(1)ALFs的插值相比经典的4种ALFs递推方法，计算耗时最小，效率是Belikov方法的2倍左右。ALFs的插值解算精度依赖插值过程中节点步长Δθ的大小选择，具体可按Nyquist采样定理设置。
(2)相比传统逐点求解，本文方法的计算耗时从3 669.41 s减小至98.05 s，在同一计算环境下效率提升近40倍；将日本地区旋转为东西狭长区域后，本文方法的应用优势更加凸显：节点步长大小为5'时，计算耗时由98.05 s减小至19.06 s，且重力异常值结果的最大误差不足0.005 mGal，最终相比逐点计算方法的3 669.41 s，效率提高近200倍。这也说明，本文方法在计算东西方向狭长分布的点集和密集程度大的点集模型值时，具有显著优势。(1)所计算点集位于同一高程面，针对实际任意高度分布的点，可进一步采用高程方向泰勒级数展开方法进行求解。(2)高阶次球谐系数的旋转变换过程有待进一步提升效率。(3)对于全球点位任意分布的稀疏离散点集或者沿纬度方向跨度较大的点集，本文方法可能效果有限，此时建议采用NFFT方法或者逐点计算方法。(4)该方法目前只适用于计算模型扰动位、重力异常、高程异常等无须计算ALFs一阶、二阶导数的重力场元，而对于计算模型垂线偏差、扰动重力、梯度等重力场元，该方法需要对ALFs的一阶、二阶导数值进行插值，其优势有待进一步验证。
![]()
![]()
