image
Home Journals Progress in Geophysics
Progress in Geophysics

Abbreviation (ISO4): Prog Geophy      Editor in chief:

About  /  Aim & scope  /  Editorial board  /  Indexed  /  Contact  / 
Online first  /  Current issue  /  All issues  /  Top access  /  RSS

Progress in Geophysics >

2025 , Vol. 40 >Issue 1: 11 - 24

DOI: https://doi.org/10.6038/pg2025II0047

Feasibility analysis of constructing a regional gravitational model using the method of spherical harmonics analysis

  • XinXing LI , 1, 2, 3 ,
  • JinKai FENG , 2, * ,
  • HaoPeng FAN 2 ,
  • XiaoGang LIU 4 ,
  • Diao FAN 2 ,
  • ShanShan LI 2
Expand
  • 1 State Key Laboratory of Geographic Information Engineering, Xi'an 710054, China
  • 2 Institude of Geospatial Information, Information Engineering University, Zhengzhou 450001, China
  • 3 Key Laboratory of Smart Earth, Beijing 100020, China
  • 4 Xi'an Research Institute of Surveying and Mapping, Xi'an 710054, China

Received date: 2024-03-29

  Online published: 2025-03-13

Copyright

Copyright ©2025 Progress in Geophysics. All rights reserved.
Fold

Abstract

Technical methods related to the regional gravity fields, include precise geoid determination, downward continuation of gravity anomalies, and construction of regional models, which are greatly faced with challenges, such as low accuracy or efficiency, computation approximations, instability, and lack of criterion. This study introduces the fast implementation technique of ultra-high degree spherical harmonics analysis into the construction of regional gravity field, combined with satellite gravity field models of low to medium frequency spectrum and GNSS/leveling data. By employing a 5400-degree spherical harmonics analysis (SHA), it is possible to achieve a regional gravity field reconstruction with a resolution of 2 arc-minutes and an error less than 1 mGal for gravity anomalies, less than 1 cm for (quasi) geoid, and less than 0.4 arc-seconds for vertical deflection. Meanwhile, SHA can also be used for stable downward continuation of airborne gravity anomalies, which is verified to be more simple, accurate, and stable. Therefore, SHA is worth promoting and applying in the modeling of regional gravity fields and related data processing problems.

Key words： Spherical harmonic analysis; Spherical harmonic synthesis; Downward continuation; Regional gravity field; Geoid; Normal ellipsoid

Cite this article

XinXing LI , JinKai FENG , HaoPeng FAN , XiaoGang LIU , Diao FAN , ShanShan LI . Feasibility analysis of constructing a regional gravitational model using the method of spherical harmonics analysis[J]. Progress in Geophysics, 2025 , 40(1) : 11 -24 . DOI: 10.6038/pg2025II0047

0 引言

利用大量重力场观测资料构建超高精度、超高空间分辨率和时间分辨率的地球重力场模型并生成其附属产品，一直是大地测量和地球物理学科领域的重要任务(李新星, 2013; Pavlis et al., 2012; 党亚民等, 2023).但受硬件、环境、政策、成本等因素的限制，使得全球高精度网格重力数据的获取并非易事.虽然目前已有学者发布了分辨率1′的重力异常网格数据(Sandwell et al., 2014; Zhu et al., 2020)，但海洋部分的真实分辨率不足2′，陆地大量重力空白区则是采用高分辨率地形数据均衡填补得到(Hirt et al., 2013; Ince et al., 2020; 王凯等, 2014; Balmino et al., 2012)，其精度水平有待评估.所以在EGM2008模型发布之后，虽然有更高阶次的超高阶全球重力场模型公布，但均并没有显著的精度提升(Ince et al., 2019).
与之相比，得力于卫星重磁的发展以及局部数据的整合，区域重磁场的建模方法层出不穷，模型精度不断提升(徐文耀等, 2011; Haines, 1985; Thébault et al., 2006; 吴晓平, 1984; 翟振和孙中苗, 2009).以大地水准面为例，许多实验山区的大地水准面精度已经达到厘米甚至毫米级，我国大量省市也已经建成了厘米级局部大地水准面模型(宁津生等, 2004; 李建成, 2012).目前局部重力场建模方法依然面临着诸多难题.
第一，局部建模方法多种多样，但各有利弊.例如纯数学拟合的方法，包括幂函数级数拟合、样条函数拟合、双三角函数拟合等，这些方法简单易行，局部拟合效果很好，但是缺少物理内涵，外推预报效果不佳；等效源、径向基函数、Slepian基函数等方法具有一定的物理含义，但结果的好坏依赖等效源和基函数的数量、位置、强度等，有时还会与实际不符而出现解释错误(Haines, 1985)；球冠谐分析(Spherical Cap Harmonic Analysis, SCHA)、矩谐分析(Rectangular Harmonic Analysis, RHA)等方法也有物理含义，但RHA方法适用于范围和高度较小的区域，不能很好的实现场元的向上向下延拓(Alldredge, 1981; Jiang et al., 2014)，SCHA方法在小区域建模中级数变长且收敛变慢甚至不准确，且不能很好的仿真场元随径向的变化，延拓效果也较差(Haines, 1985)；修正SCHA (Revised SCHA, R-SCHA) 相比SCHA虽然能够较好的改善场元径向的变化，但是其引入的Mehler圆锥函数会使得余纬方向的边界效应增大(Thébault et al., 2006).
第二，大部分局部重力场模型只能描述单一重力场元要素，不能很好的解决多种重力场元之间的相互转换，且三维空间数据的延拓不严密.例如局部大地水准面模型与局部重力异常数据之间的关系依赖于Stokes积分，航空重力数据的向下延拓主要基于Poisson积分进行等；同时局部重力场模型构建和数据处理过程中面临诸多细节性难题，例如前述Stokes积分和Poisson积分面临着存在积分离散化误差、非全球积分以及数据空白带来的远区效应、积分过程中积分半径大小的选择、中央网格积分奇异性、正则化过程中因子参数的确定等问题，均与所构建区域地形、重力数据复杂程度和模型构建过程紧密相关，处理方法和参数没有统一标准(欧阳永忠等, 2011).以航空重力数据向下延拓这一问题为例，学者们提出了最小二乘法、Tikhonov正则化法、直接代表法、逆Poisson积分迭代法、解析延拓法、虚拟点质量法等诸多方法(王兴涛等, 2004)；为改进上述方法的精度和稳定性，学者们又陆续作了大量工作，例如解决航空数据的边界效应、对Poisson积分核的改化和带限处理、分离剩余地形产生的高频信号、确定正则化参数、顾及延拓算子非线性项、引入FFT(Fast Fourier Transform)方法加快计算等，但一定程度上也使处理流程和计算公式更加复杂(蒋涛等, 2018; 黄谟涛等, 2018, 2022; 刘晓刚等, 2018; 马健等, 2023).
球谐分析方法及全球重力场模型的构建在整个20世纪得到了较快的发展，随着硬件和算力的提升，主流方法由早期效率较高的数值积分或调和分析方法(Rapp and Pavlis, 1990)逐渐发展为理论更加严密的最小二乘方法(Pavlis et al., 2012)，但是随着地形数据的极大丰富，对球谐分析的阶次提出了更高要求，所以在FFT方法的提速下，数值积分方法成为10800阶以上模型构建的核心方法(Hirt and Rexer, 2015)，田家磊等(2018)通过改进块对角最小二方法，使得目前利用块对角最小二乘方法实现2160甚至更高阶次的球谐分析(Spherical Harmonic Analysis, SHA)已经不再是难事.进一步结合高精度的区域重力资料，利用SHA方法解决局部地区重力场建模问题也成为一种选择.重力场模型“剪裁法”利用了局部重力资料精化地球重力场模型从而更适用于该局部区域，本质不属于局部建模而是全球建模问题(陆洋和许厚泽, 1994; 许厚泽等, 2002)；梁伟等(2023)使用SHA联合航空重力数据构建了适用于局部区域的超高阶重力场模型，虽然其利用SHA方法改善了局部区域重力场，但本质上还是全球重力场建模问题；吴文京(1984)中曾指出，SHA方法不能应用于局部地磁场建模，该结论的得出是基于其仅仅利用局部观测数据资料来建模型，并没有对局部区域以外的测点分布进行约束，所以无法实现球谐系数的稳定求解；丁子凡等(2022)使用SHA方法较好的解决了重力异常径向方向各阶次导数的计算，并通过迭代实现重力异常的向下延拓，这里SHA只是发挥了辅助作用.
为了验证SHA在局部重力场模型构建中的作用，本文从Laplace方程的球谐级数展开出发，引入SHA的最小二乘快速实现方法，并完成对局部区域重力场多场元要素的逼近、航空重力异常数据的向下延拓等内容，以验证SHA在局部重力场模型构建中应用的可行性及其优势，给出有益结论.

1 球谐综合与分析

1.1 球谐综合

球谐函数是在球坐标系统下表示的Laplace方程的正交解集，Laplace方程表达式为(Hofmann-Wellenhof and Moritz, 2006)：
$\Delta F=0, $
其中“Δ”表示Laplace算子，F是空间的一个标量函数，满足式(1)的F又称为调和函数.地球重力场理论中，地球引力位V在地球外部空间是调和的，可展开为如下的球谐级数形式：
$\begin{aligned}& V(r, \theta, \lambda)= \frac{\mathrm{GM}}{r} \sum\limits_{n=0}^{\infty}\left(\frac{R}{r}\right)^n \sum\limits_{m=0}^n\left(\bar{C}_{n m} \cos m \lambda+\right. \\&\left.\bar{S}_{n m} \sin m \lambda\right) \bar{P}_{n m}(\cos \theta), \end{aligned}$
式中(r, θ, λ)表示所计算点在球坐标下的的向径、余纬和经度；R表示参考半径，其大小与计算过程中选择的参考椭球没有关系，只和使用的地球重力场模型有关，以EGM2008模型和XGM2019e模型为例，R的大小均为6378136.3 m；GM为地心引力常数；n和m分别表示球谐展开的阶(degree)和次(order)；Pnm(cosθ)表示完全归一化缔合勒让德函数；(Cnm, Snm)为对应球谐展开的高斯系数，又称为球谐系数.
进一步，引入正常椭球后根据重力场元与扰动位间的泛函关系，可得到球近似下的空间重力异常Δg、高程异常ζ、垂线偏差(ξ, η)、扰动重力δg的计算公式如下：
$\left\{\begin{array}{l}\Delta g=\frac{\mathrm{GM}}{R^2} \sum\limits_{n=1}^{\infty}(n-1)\left(\frac{R}{r}\right)^{n+2} \sum\limits_{m=0}^n\left(\bar{C}_{n m}^* \cos m \lambda+\bar{S}_{n m} \sin m \lambda\right) \bar{P}_{n m}(\cos \theta) \\\zeta=\frac{\mathrm{GM}}{r \gamma} \sum\limits_{n=1}^{\infty}\left(\frac{R}{r}\right)^n \sum\limits_{m=0}^n\left(\bar{C}_{n m}^* \cos m \lambda+\bar{S}_{n m} \sin m \lambda\right) \bar{P}_{n m}(\cos \theta) \\\xi=\frac{\mathrm{GM}}{r R \gamma} \sum\limits_{n=1}^{\infty}\left(\frac{R}{r}\right)^n \sum\limits_{m=0}^n\left(\bar{C}_{n m}^* \cos m \lambda+\bar{S}_{n m} \sin m \lambda\right) \frac{\mathrm{d} \bar{P}_{n m}(\cos \theta)}{\mathrm{d} \theta} \text{ }\text{ }\text{ }\text{ }\text{ }\text{ }\text{ }\text{ }\text{ }\text{ }\text{ }\text{ }, \\\eta=\frac{\mathrm{GM}}{r R \gamma \sin \theta} \sum\limits_{n=1}^{\infty}\left(\frac{R}{r}\right)^n \sum\limits_{m=0}^n m\left(\bar{C}_{n m}^* \sin m \lambda-\bar{S}_{n m} \cos m \lambda\right) \bar{P}_{n m}(\cos \theta) \\\delta g=\frac{\mathrm{GM}}{r^2} \sum\limits_{n=1}^{\infty}(n+1)\left(\frac{R}{r}\right)^n \sum\limits_{m=0}^n\left(\bar{C}_{n m}^* \cos m \lambda+\bar{S}_{n m} \sin m \lambda\right) \bar{P}_{n m}(\cos \theta)\end{array}\right.$
其中，γ表示计算点处的正常重力值，可利用椭球参数的点位坐标计算得到；ξ、η分别表示垂线偏差的子午分量和卯酉分量；φ表示地心纬度，φ=90°-θ. 在球近似下，计算球面上的高程异常与大地水准面高的表达式是一致的，此处不再区分.以上计算过程均可称为球谐综合(Spherical Harmonic Synthesis, SHS).

1.2 球谐分析

地球重力场模型的构建是利用地球空间重力场元观测资料求解球谐系数的过程，称为球谐分析(SHA)，是球谐综合的逆过程.该过程求解包括两种成熟方法，一是数值积分，又称为调和分析方法；另一种是最小二乘方法.
(1) 数值积分方法
目前国际上5400阶次以上的模型大部分采用数值积分方法，该方法的特点是计算效率相对较高，但计算精度略低.学者分别引入了Gauss网格和D-H网格的离散化方法提升解算精度，同时为提升计算效率，引入了1D-FFT、2D-FFT方法(Wieczorek and Meschede, 2018).
已知覆盖半径R球面的重力异常信息，即公式中r=R，则球谐系数可逐个利用式(4)求得：
$\begin{array}{l}\left.\begin{array}{l}\bar{C}_{n m} \\\bar{S}_{n m}\end{array}\right\}=\frac{R^2\left(1+\delta_0^m\right)}{4 \mathtt{π} G M(n-1)} \iint_\mathtt{σ} \Delta g(\theta, \lambda) \bar{P}_{n m}(\cos \theta)\left\{\begin{array}{l}\cos m \lambda \\\sin m \lambda\end{array}\right\} \mathrm{d} \sigma .\\\begin{array}{c}\approx & \frac{R^2\left(1+\delta_0^m\right)}{4 \mathtt{π} G M(n-1)} \sum\limits_{i=1}^{\text {row }} \sum\limits_{j=1}^{\text {column }} \iint_\mathtt{σ} \Delta g\left(\theta_i, \lambda_j\right) \\& \bar{P}_{n m}\left(\cos \theta_i\right)\left\{\begin{array}{l}\cos m \lambda_j \\\sin m \lambda_j\end{array}\right\} \sin \theta_i \Delta \theta \Delta \lambda,\end{array}\end{array}$
式中，δ0m函数表示当m=0时值为1，m取其他值时均为零；dσ=sinθdθdλ，row表示重力异常网格数据的总纬度行数，column表示总经度列数，θi表示第i行网格数据的余纬值，λi表示第j列网格数据的余纬值，Δθ、Δλ表示数据在经纬度方向的分辨率.
(2) 最小二乘方法
相比数值积分方法，最小二乘方法顾及了球谐系数中“次m”相同的系数间的相关性，因而求解精度较高.但该方法计算效率低，即使使用了“块对角”稀疏化，计算量依然庞大.对此，Sneeuw (1994)通过分离球谐因子在经度方向和纬度方向的正交性，改化了系数矩阵的结构从而提升块对角最小二乘方法的解算效率.
求解过程中，重力异常数据通常按地理网格全球分布，因面球谐因子的正交性，法矩阵N呈稀疏状态，且Cnm与Snm之间、Cnm中“次m”不相等的系数之间、Snm中“次m”不相等的系数之间、“阶n”为偶数与“阶n”为奇数的系数之间均是独立不相关的，所以对于最高阶次Nmax模型的SHA，法方程可以拆分为4Nmax个法方程并分别求解：
$\left\{\begin{array}{ll}\boldsymbol{N}_{m, c 1} \hat{\boldsymbol{X}}_{m, c 1}+\boldsymbol{U}_m=0 & n( { odd }), \left(m=0, 1, 2, \cdots, N_{\max }-1\right) \\\boldsymbol{N}_{m, c 2} \hat{\boldsymbol{X}}_{m, c 2}+\boldsymbol{U}_m=0 & n( { exen}), \left(m=0, 1, 2, \cdots, N_{\max }-1\right) \\\boldsymbol{N}_{m, s 1} \boldsymbol{X}_{m, s 1}+\boldsymbol{U}_m=0 & n( { odd }), \left(m=1, 2, \cdots, N_{\max }\right) \\\boldsymbol{N}_{m, s 2} \hat{\boldsymbol{X}}_{m, s 2}+\boldsymbol{U}_m=0 & n( { exen}), \left(m=1, 2, \cdots, N_{\max }\right)\end{array}, \right.$
其中Nm, c1表示与Cnm中n为奇数，m为固定值的所有球谐系数对应的法矩阵，$ \hat{\boldsymbol{X}}_{m, c1}$表示球谐系数待求向量，其余与之类似，Um表示法方程的自由项向量，相关计算和矩阵存储方式详见(田家磊等., 2018).

2 SHA在局部重力场建模中的应用

目前普遍认为，只有局部重力资料而没有全球分布的相当精度和分辨率的重力场资料是无法构建全球重力场模型的，所以SHA一直被认为是一种全球重力场建模方法(吴文京, 1984)，更多的是作为精化局部重力场的辅助方法，例如构建地球重力场模型的“剪裁法”等(陆洋和许厚泽, 1994; 梁伟等, 2023).
下面，在验证SHA算法严密性的基础上，将SHA应用于局部重力场建模，分析其建模精度及其优劣，给出相应的使用建议，并在航空重力数据延拓方面给出有益的应用体验.

2.1 SHA算法精度评估

为验证本文所提SHA方法的正确性，使用目前国际上发布的阶次最高的5400阶次重力场模型XGM2019e(Zingerle et al., 2020)为参考值，利用公式计算了全球(R=6378136.3) 球面2′分辨率地理网格中点的重力异常值，如图 1所示.
图1 XGM2019e模型计算的球近似下的全球空间重力异常

Fig 1 Global spherical approximated free-air gravity anomaly calculated from the XGM2009e model

在表 1所示计算环境下，利用块对角最小二乘方法，反算5400阶次的球谐系数XGM2019e_R，耗时约6.0 h，计算过程中进程内存最大1.3 GB，对于当前5400阶次的SHA计算量来说，个人计算机性能完全能够满足，移植到服务器上并加入并行计算后，实现目前国际上SHA达到的最高阶次43200也不是难事(Rexer and Hirt, 2015)，而且本文使用的最小二乘方法有区别于国际上使用的软件SHTOOLS和数值积分方法，在计算精度、精度估计，以及卫星模型联合求解方面更具有优势(Wieczorek and Meschede, 2018).
表1 计算环境参数

Table 1 Computing environment parameters

计算环境 参数
操作系统 Windows10 64bit
CPU Intel(R) Core(TM) i7-9750H CPU @ 2.60 GHz
内存/GB 32
线程数 1
软件环境 VS 2019/C++/Eigen/Intel MKL
以XGM2019e模型为真值，对重建的XGM2019e_R进行比较分析，恢复系数与原系数真值差异如图 2所示.同时，对重建模型系数的误差的“阶RMS”进行统计，结果如图 3所示.计算公式为：
$\sigma_n=\sqrt{\frac{1}{2 n+1} \sum\limits_{m=0}^n\left(\bar{C}_{n m}^R-\bar{C}_{n m}^{\mathrm{XGM}}\right)^2+\left(\bar{S}_{n m}^R-\bar{S}_{n m}^{\mathrm{XGM}}\right)^2}, $
图2 SHA的球谐系数误差分布图

Fig 2 Distribution of spherical harmonic coefficient errors of SHA

图3 SHA的球谐系数误差的阶RMS图

Fig 3 Order Root-mean-square plot of spherical harmonic coefficient errors of SHA

其中，CnmXGM, SnmXGM是XGM2019e模型系数，CnmR, SnmR是重建的模型系数.
以上实验结果表明，5400甚至更高阶次SHA的实现和精度方面均已经达到了较高水平，也为用于构建局部重力场模型提供了可能.

2.2 局部球谐模型的建立

同样使用XGM2019e模型作为“真值”，选取我国西南部山区8°×10°范围(27°N~35°N, 90°E~100°E)开展仿真实验，其地形如图 4所示.
图4 实验区地形

Fig 4 The terrain of experimental area

首先计算该区域内球面上2′分辨率地理网格中点处的重力异常Δg0作为“观测值”，同时计算区域内高程异常ζ0、垂线偏差(ξ0, η0)值以及3000 m、8000 m高度处的重力异常值Δgair3000、Δgair8000作为局部重力场精化的目标真值，数据大小均为240行×300列，其仿真结果如图 5所示.
图5 局部重力场各场元仿真结果

(a)地面重力异常；(b)3000 m重力异常；(c)8000 m重力异常；(d)地面高程异常；(e)地面垂线偏差子午分量；(f)地面垂线偏差卯酉分量.

Fig 5 The simulation results of various field elements of the local gravitational field

(a) Ground gravity anomaly; (b) 3000 m gravity anomaly; (c) 8000 m gravity anomaly; (d) Abnormal ground elevation; (e) Ground vertical deviation meridian component; (f) The vertical deviation of the ground is the Maoyou component.

对区域重力异常观测量Δg0之外的所有区域，用2′分辨率网格的“0”值重力异常填补，以此构成一套全球2′网格重力异常数据ΔgModel10；随后对其进行SHA生成表征该局部区域重力场的5400阶次球谐系数模型Model1.需要注意该模型球谐系数已不再具有严格的物理含义，只能说是对该区域空间重力场的最佳一致逼近.
利用Model1重建该区域相应重力场元Δg1、ζ1、(ξ1, η1)、Δgair13000和Δgair18000，与原仿真真值比较，统计结果及差异值分布情况分别如表 2和图 6所示.
表2 Model1模型重建局部重力场误差统计表

Table 2 Error statistics of the reconstructed local gravity field using Model1

统计量/单位 数据量 Max Min Mean Std Rms
Δg1-Δg0(0 m)/mGal 240×300 19.217 -19.731 -0.090 0.567 0.575
Δgair1-Δgair(3000 m)/mGal 240×300 58.108 -63.672 -0.012 2.419 2.419
Δgair1-Δgair(8000 m)/mGal 240×300 72.672 -62.945 0.105 3.690 3.692
ζ1-ζ0/m 240×300 53.832 34.768 41.580 2.630 41.663
ξ1-ξ0/arc sec(″) 240×300 56.610 -10.017 1.141 4.056 4.214
η1-η0/arc sec(″) 240×300 22.902 -3.795 5.350 1.370 5.523
图6 Model1模型重建局部重力场误差分布图

(a)地面重力异常；(b)3000 m重力异常；(c)8000 m重力异常；(d)地面高程异常；(e)地面垂线偏差子午分量；(f)地面垂线偏差卯酉分量.

Fig 6 Distribution of reconstructed local gravity field errors using Model1

(a) Ground gravity anomaly; (b) 3000 m gravity anomaly; (c) 8000 m gravity anomaly; (d) Abnormal ground elevation; (e) Ground vertical deviation meridian component; (f) The vertical deviation of the ground is the Maoyou component.

由表 2统计结果可见，由于区域以外全球重力资料的缺失，导致SHA构建的模型相对原局部重力场重建的精度较差，不能够满足高精度应用的需求；但是通过图 6的误差分布图可见，地面及空中重力异常中，误差较大的数据统一分布在边界区域，这是由边界效应产生的，因此局部重力场模型只能适用于扣除一定边界的中间区域.
分别对上述场元扣除重建区域0.5°、1.0°、2°的边界，对剩余中间区域进行误差的统计，结果见表 3.SHA方法在扣除2°的边界效应后能够以1.0 mGal以内的误差实现山区空间重力异常的重建.
表3 Model1模型重建局部重力场扣除边界区域后的误差统计表

Table 3 Error statistics of the reconstructed local gravity field using Model1 after deducting a certain boundary

统计量 扣除边界大小 数据量 Max Min Mean Std Rms
Δg1-Δg0(0 m)/mGal 0.5° 210×270 2.869 -2.916 -0.090 0.258 0.273
1.0° 180×240 1.503 -1.906 -0.089 0.232 0.248
2.0° 120×180 0.781 -1.020 -0.089 0.204 0.223
Δgair1-Δgair(3000 m)/mGal 0.5° 210×270 1.389 -1.117 -0.041 0.170 0.175
1.0° 180×240 0.411 -0.349 -0.045 0.077 0.089
2.0° 120×180 0.066 -0.163 -0.048 0.030 0.057
Δgair1-Δgair(8000 m)/mGal 0.5° 210×270 1.137 -2.346 -0.022 0.266 0.267
1.0° 180×240 1.137 -0.549 0.030 0.199 0.202
2.0° 120×180 0.252 -0.159 0.021 0.068 0.071
ζ1-ζ0/m 0.5° 210×270 48.705 36.554 41.464 2.027 41.514
1.0° 180×240 46.479 37.532 41.379 1.599 41.410
2.0° 120×180 43.859 38.925 41.275 0.987 41.286
ξ1-ξ0/arc sec(″) 0.5° 210×270 11.010 -3.964 0.757 2.527 2.638
1.0° 180×240 6.894 -2.829 0.643 1.976 2.078
2.0° 120×180 3.697 -1.682 0.559 1.218 1.340
η1-η0/arc sec(″) 0.5° 210×270 11.496 1.619 5.375 1.080 5.482
1.0° 180×240 9.390 3.154 5.356 0.885 5.428
2.0° 120×180 6.722 3.858 5.307 0.610 5.342
对(似)大地水准面和垂线偏差来说，扣除大的边界并不能有效降低相应误差水平，说明除了边界效应以外，该误差还包含了大的系统性差异，分析其原因是区域以外没有数据导致重力场中低频信号缺失造成的.
考虑到目前纯卫星重力场模型阶次已高至300阶次(赵永奇等, 2023; Göttl et al., 2018; Flechtner et al., 2010)，借鉴“移去-恢复”技术，分别使用XGM2019e模型前36阶和300阶次系数，对实验区域以外的2′网格的重力异常值进行填补，生成全球2′网格重力异常数据集ΔgModel236和ΔgModel3300；然后使用SHA构建该区域局部重力场球谐模型Model2和Model3.利用这两个模型重建区域内相应重力场元并再次与参考值比较，其误差结果统计见表 4.Model3重建的局部重力场误差扣除2°边界后的分布如图 7所示.
表4 Model2、Model3模型重建局部重力场误差统计表

Table 4 Error statistics of the reconstructed local gravity field using Model2 and Model3

统计量 扣除边界大小 数据量 Max Min Mean Std Rms
Δg236-Δg0(0 m)/mGal 0° 240×300 18.405 -18.864 -0.022 0.543 0.543
2.0° 120×180 0.885 -0.945 -0.021 0.203 0.205
Δg3300-Δg0(0 m)/mGal 0° 240×300 11.439 -12.522 -0.001 0.463 0.463
2.0° 120×180 0.872 -0.945 -0.001 0.200 0.200
ζ236-ζ0/m 0° 240×300 8.396 -2.463 1.030 1.178 1.564
2.0° 120×180 2.038 0.281 0.986 0.418 1.070
ζ3300-ζ0/m 0° 240×300 1.166 -0.664 0.064 0.066 0.092
2.0° 120×180 0.103 0.013 0.056 0.014 0.058
ξ236-ξ0/arc sec(″) 0° 240×300 48.253 -9.891 1.111 2.468 2.706
2.0° 120×180 1.606 -0.019 0.492 0.266 0.559
ξ3300-ξ0/arc sec(″) 0° 240×300 12.323 -16.475 -0.008 0.523 0.523
2.0° 120×180 0.159 -0.150 0.003 0.034 0.035
η236-η0/arc sec(″) 0° 240×300 19.116 -6.782 0.368 1.506 1.551
2.0° 120×180 1.793 -1.367 0.090 0.559 0.566
η3300-η0/arc sec(″) 0° 240×300 7.516 -7.600 0.005 0.323 0.323
2.0° 120×180 0.235 -0.142 0.010 0.037 0.038
Δgair236-Δgair(3000 m)/mGal 0° 240×300 53.611 -58.275 0.107 2.154 2.156
2.0° 120×180 0.099 -0.089 -0.007 0.024 0.025
Δgair3300-Δgair(3000 m)/mGal 0° 240×300 31.599 -21.087 0.007 0.810 0.810
2.0° 120×180 0.070 -0.074 -0.0005 0.016 0.016
图7 Model3模型重建局部重力场误差分布图

(a)重力异常；(b)高程异常；(c)垂线偏差子午分量.

Fig 7 Distribution of reconstructed local gravity field errors using Model3

(a) Gravity anomaly; (b) Abnormal elevation; (c) The meridian component of vertical deviation.

比较表 3与表 4结果可见，外部区域36阶和300阶模型的填充并没有显著提升地面重力异常的重建精度，说明地面重力异常重建的误差是高频信号混频引入的，该结论与图 7a显示结果一致，本文分析该误差主要来自于区域边界数据的突变；而300阶模型对局部以外区域的填补对于小于10 cm误差的高程异常和小于0.4″精度垂线偏差的重建是不可或缺的；垂线偏差的重建精度已经完全满足要求，而且其误差中也主要是高频信号，由图 7c可见；而对于似大地水准面的重建，显然10 cm误差对于现代局部高程基准来说是非常大的，不能不予考虑.由图 7b可见，高程异常重建的误差具有明显的趋势项，说明本文方法能够很好的重建似大地水准面高频信号，但由于区域以外300阶以上重力场大量信号的缺失，导致对该区域内高程异常造成长波项误差.传统局部(似)大地水准面构建均会引入GNSS/水准数据做控制，其主要控制的就是该长波项内容.
假设在局部区域内均匀分布了9个GNSS/水准点，能够将该项长波误差准确的标定出来(这里直接取9个点位处的重建误差值)，进行曲面拟合后对该区域重建的高程异常误差进行Δζ修正，其修正后的精度统计情况如表 5和图 8所示.
表5 Model3模型重建局部高程异常误差统计表(单位：m)

Table 5 Error statistics of the reconstructed local height anomaly using Model3 (unit: m)

统计量 数据量 Max Min Mean Std Rms
ζ3300-ζ0+Δζ 120×180 0.0058 -0.0069 -0.0005 0.0013 0.0014
图8 9个点控制修正后Model3模型重建高程异常的误差分布图

Fig 8 Distribution of reconstructed local gravity field errors using Model3 after 9-point control corrections

由上述结果可知，利用SHA重建局部高程异常的精度在GNSS/水准数据控制的情况下误差不超过1 cm，且最终精度取决于GNSS/水准的精度，由图 8可见，剩余误差依然呈趋势变化，可知进一步提升高程异常精度的方法主要是增加外部区域重力场的信号或者增加GNSS/水准数据量，而不是增加重力测量数据的分辨率.
至此，本文联合区域2′网格重力异常数据，300阶次卫星重力场模型和区域内9个GNSS/水准数据，利用SHA方法成功构建了该区域重力场模型，使得扣除2°边界区域后的中央区域，单因方法造成的重力异常误差 < 1 mGal、(似)大地水准面误差 < 1 cm、垂线偏差误差 < 0.4″，充分说明了SHA在局部重力场建模中具有良好的表现，需要注意的是上述实验均为理论验证实验，数据中不含任何测量误差.

2.3 航空重力数据向下延拓中的应用

在局部重力场精化中，航空重力数据延拓也是重点关注内容之一.下面通过实验，观察SHA构建局部重力场模型在航空重力数据向下延拓问题中的表现.
依然选取27°N~35°N, 90°E~100°E(8°×10°)范围开展仿真实验，利用XGM2019e模型计算球面以上3000 m处分辨率2′的重力异常值Δgair3000作为观测资料，球面0 m高度处重力异常Δg0作为向下延拓的“真值”.首先利用XGM2019e模型前300阶次系数计算填补区域外侧3000 m高度球面空间2′的重力异常生成全球3000 m高度处2′网格重力异常数据ΔgModel43000，然后利用SHA方法生成该局部区域重力场球谐系数模型Model4，观测方程计算过程中取r=R+3000.Model4模型可直接计算该范围内球面高度0 m处的2′网格重力异常ΔgModel40，其与“真值”Δg0的差异即可视为延拓误差，因边界效应，延拓结果依然是在扣除一定边界后的中央区域正确，统计结果见表 6，中央区延拓的误差分布如图 9所示.
表6 SHA方法延拓无噪声局部航空重力数据误差统计表

Table 6 Error statistics table for downward continuation of noise-free local airborne gravity anomalies by SHA

统计量 扣除边界大小 数据量 Max Min Mean Std Rms
ΔgModel40-Δg0/mGal 0° 240×300 116.462 -25.488 0.0007 4.147 4.147
0.5° 210×270 8.308 -8.135 -0.003 0.704 0.704
1.0° 180×240 4.031 -3.907 -0.001 0.538 0.538
2.0° 120×180 1.793 -1.893 -0.002 0.413 0.413
图9 SHA方法延拓局部航空重力异常的误差分布图

Fig 9 Distribution of errors in downward continuation of local airborne gravity anomalies by SHA

由统计和图 9所示结果可见，延拓结果中有2 mGal以内的延拓误差，均方根为0.413 mGal，图中可见这些误差是高频噪声，由局部区域外侧300阶次以上重力场信号缺失造成，所以该部分误差只能通过增加全球重力场信号的方法予以削弱，例如增加外部空间填补重力异常所使用的重力场模型的阶次的方式.
在实际航空重力数据中往往是包含测量误差的，因此在向下延拓过程中必然会将测量误差放大，为了分析噪声在向下延拓中的影响，在上述Δgair3000重力异常中加入3 mGal的白噪声，并采用同样的填补方式得到ΔgModel53000，采用SHA方法生成5400阶次的球谐模型Model5，利用5400阶次Model5可直接计算球面0 m的重力异常ΔgModel50，通过与目标真值Δg0比较，结果见表 7，发现噪声被放大非常严重，噪声高达17.5 mGal，显然已经远远超出可使用范围.
表7 SHA方法向下延拓具有3 mGal噪声的局部航空重力数据误差统计表

Table 7 Error statistics table for downward continuation of local airborne gravity anomalies with 3 mGal noise by SHA

统计量 截断阶次 Max Min Mean Std Rms
ΔgModel50-Δg0/mGal 5400 68.259 -68.791 -0.028 17.520 17.519
ΔgModel50(4320)-Δg0/mGal 4320 51.176 -53.418 -0.027 10.581 10.580
ΔgModel50(3600)-Δg0/mGal 3600 67.196 -72.251 -0.029 9.747 9.747
ΔgModel50(3240)-Δg0/mGal 3240 76.301 -88.343 -0.030 10.112 10.112
为了抑制航空重力数据对地面高频噪声的污染，最有效的方式就是将噪声超过信号本身的高频频段予以舍弃，以损失小的真实信号为代价换取大的噪声抑制效果.为实现该目标，对Model5模型分别进行不同阶次的截断，即计算地面2′重力异常值时最高阶次分别只取到4320、3600和3240，得到球面0 m的重力异常ΔgModel50(4320)、ΔgModel50(3600)和ΔgModel50(3240).经过比较分析可以发现，通过对Model5模型高阶系数进行截断，可以有效抑制向下延拓过程中对航空重力数据高频噪声的放大，但是损失的是延拓到地面后重力异常数据的空间分辨率，以3600阶次截断为例，相当于2′的航空重力数据延拓到地面后分辨率降为3′，这是抑制噪声所必须付出的代价.从表 7结果看，延拓结果即使因边界效应扣除2°的边界，误差依然高达±70 mGal，结果看似不正确，与许多航空重力向下延拓结果精度有出入，但实际情况就是如此.
为了说明该问题，下面增加一组实验.利用XGM2019e模型的前2160阶次模型仿真该区域3000 m高度处和地面0 m高度处2′分辨率的重力异常ΔgX21603000和ΔgX21600，在ΔgX21603000中加入3 mGal白噪声后作为航空重力数据观测值，使用XGM2019e前300阶次系数模型填补测区外2′的重力异常值，得到全球3000 m高度重力异常ΔgModel63000，然后使用SHA方法生成5400阶次局部重力场球谐模型Model6，利用截断至2160阶次的Model6计算地面0 m高程2′的重力异常值ΔgModel60(2160)完成延拓过程.
根据上一节中的结论，外部区域使用36阶和300阶模型的填充不会显著的提升地面重力异常的重建精度，所以在只有航空重力数据向下延拓的任务中，也无需对外部区域进行重力异常填补，即补“0”值即可.上述实验中ΔgX21603000中加入3 mGal白噪声后作为航空重力数据观测值，区域外全部直接补充2′的“0”值，得到全球重力异常ΔgModel73000，依次重复Model6相同的过程，得到延拓结果ΔgModel70(2160)，两种方式延拓统计结果见表 8.由结果可见，在采用SHA方法进行向下延拓时，对航空重力数据以外补“0”值也能够实现较高精度的向下延拓而无需依赖于中低阶卫星重力场模型对空白区进行填补.
表8 截断至2160阶的SHA方法向下延拓具有3 mGal噪声的航空重力数据的误差统计表

Table 8 Error statistics table for downward continuation of local airborne gravity anomalies with 3 mGal noise by SHA truncated to 2160 degree

统计量 Max Min Mean Std Rms
ΔgModel60(2160)-Δg0/mGal 7.815 -7.795 -0.027 2.061 2.061
ΔgModel70(2160)-Δg0/mGal 7.603 -8.349 -0.152 2.117 2.122
结果说明，在3000 m高度2′重力异常加入3 mGal噪声的情况下，向下延拓后与地面“真值”比较最大差异在10 mGal以内，且差值的均方根小于3 mGal，该结果与目前最新向下延拓仿真实验成果一致且更优(黄谟涛等, 2022)，同时更加稳定.
对比分析Model5、Model6两个实验可知，Model5中有效信号是5400阶次，而噪声添加在2′网格航空重力数据中，相当于在第5400阶这一频段加入了大的噪声，SHA过程会强制将噪声与相应信号混淆并通过混频影响到其他阶次信号，而在向下延拓过程中又不可避免的会将高频信号中的误差进行放大，结果是不可接受的.但是在Model6中，虽然加入的噪声依然是5400阶这一频段，但是有效信号只有2160阶，所以在SHA过程中会强制将航空重力数据中的大部分噪声吸收到2160~5400阶次系数中，而Model6延拓结果只使用了0~2160阶次信号而没有使用包含噪声及其放大信号的2160~5400的系数，所以有较好的延拓效果.
综上，有效提升航空重力数据向下延拓精度并抑制噪声的手段，是将有用信号与噪声信号的频带进行显著区分，目前许多延拓方法中顾及地形效应(刘敏等, 2016)，先移去地形引力信号，延拓后再恢复地形引力信号，就是将有用信号中的高频部分移除，使得被延拓的有用信号与高频噪声在频带上予以了有效区分.

3 结论

为解决局部重力场模型构建精度难题，本文尝试将超高阶次的SHA方法应用到局部重力场建模中，在个人PC机上实现了5400阶次球谐模型的高精度构建，并将其应用于局部重力场建模和航空重力异常向下延拓中，均取得了很好的应用效果，并得到了相关有益经验.通过仿真实验，得出以下结论：
(1) 5400阶次的高精度SHA能够在PC计算机上实现，内存要求1.3 GB，耗时约6CPUh.
(2) SHA应用于2′分辨率的局部重力场模型构建，需要联合300阶次的卫星重力场模型对局部外侧空白区域进行数据的填补，重建的局部球谐模型无法解决边界效应难题，扣除2°边界后的中央区重建精度较好，各个高度空间重力异常重建误差小于1 mGal，垂线偏差重建误差小于0.4″，(似)大地水准面重建误差小于10 cm且其误差主要来自外部区域300阶以上信号的缺失，该误差表现为系统性中长波误差，配合GNSS/水准做控制后可达到mm级精度.
(3) SHA方法能够应用于航空重力异常向下延拓，可有效抑制高频噪声的放大且能够充分顾及向下延拓过程的非线性特征，仅从方法而言是目前积分延拓方法顾及远区效应、核改化、核带限等多种因素后的极限表现，且与迭代方法相比，也是一种简单且稳健的方法，值得在后续应用中推广.
(4) 对于具有噪声的航空重力数据的向下延拓，必须以牺牲数据的高频有效信号为代价换取对相应高频噪声的抑制，这一矛盾是不可调和的，所以在实际应用中通常利用高分辨率地形数据将高频信号进行提取，实现延拓的信号与高频噪声在频段上有效分离，然后向下延拓过程抑制高频信号和噪声，最后再补充地形部分的高频信号实现高精度的向下延拓.
本文实验过程还存在如下局限性：
(1) 整个过程是在球近似下展开的，而更严密的应该采用椭球谐分析，其中差异有待进一步研究分析；
(2) 局部重力场模型构建在球面上进行，而没有考虑起地形的起伏与其质量效应.
综上所述，SHA方法结合卫星重力场模型可实现高精度的局部重力场重建，方法简单且稳定性和通用性强，其缺点在于计算量大，浪费计算资源.

感谢审稿专家提出的修改意见和编辑部的大力支持！

References

Publishing order | Descend order by publishing year | Descend order by cited within
Alldredge L R. Rectangular harmonic analysis applied to the geomagnetic field. Journal of Geophysical Research: Solid Earth, 1981, 86(B4): 3021- 3026.

Balmino G, Vales N, Bonvalot S, et al. Spherical harmonic modelling to ultra-high degree of Bouguer and isostatic anomalies. Journal of Geodesy, 2012, 86(7): 499- 520.

DOI

Dang Y M, Jiang T, Yang Y X, et al. Research progress of geodesy in China (2019-2023). Acta Geodaetica et Cartographica Sinica, 2023, 52(9): 1419- 1436.

Ding Z F, Li J C, Xu X Y. Analytical downward continuation of gravity anomaly based on global spherical harmonic expansion iteration. Science of Surveying and Mapping, 2022, 47(10): 8- 14.

Flechtner F, Dahle C, Neumayer K H, et al. 2010. The release 04 CHAMP and GRACE EIGEN gravity field models. //Flechtner F M, Gruber T, Güntner A, et al eds. System Earth via Geodetic-Geophysical Space Techniques. Berlin, Heidelberg: Springer, 41 -58.

Göttl F, Schmidt M, Seitz F. Mass-related excitation of polar motion: an assessment of the new RL06 GRACE gravity field models. Earth, Planets and Space, 2018, 70(1): 195

DOI

Haines G V. Spherical cap harmonic analysis. Journal of Geophysical Research: Solid Earth, 1985, 90(B3): 2583- 2591.

DOI

Hirt C, Claessens S, Fecher T, et al. New ultrahigh-resolution picture of Earth's gravity field. Geophysical Research Letters, 2013, 40(16): 4279- 4283.

DOI

Hirt C, Rexer M. Earth2014: 1 arc-min shape, topography, bedrock and ice-sheet models-available as gridded data and degree-10, 800 spherical harmonics. International Journal of Applied Earth Observation and Geoinformation, 2015, 39: 103- 112.

DOI

Hofmann-Wellenhof B, Moritz H. 2006. Physical Geodesy. 2nd ed. Vienna: Springer.

Huang M T, Liu M, Deng K L, et al. Analytical solution of downward continuation for airborne gravimetry based on upward continuation method. Chinese Journal of Geophysics, 2018, 61(12): 4746- 4757.

DOI

Huang M T, Deng K L, Wu T Q, et al. Rigorous modification model of upward continuation and its applications on the downward continuation of gravity anomaly. Acta Geodaetica et Cartographica Sinica, 2022, 51(1): 41- 52.

Ince E S, Barthelmes F, Reißland S, et al. ICGEM-15 years of successful collection and distribution of global gravitational models, associated services, and future plans. Earth System Science Data, 2019, 11(2): 647- 674.

DOI

Ince E S, Abrykosov O, Förste C, et al. Forward gravity modelling to augment high-resolution combined gravity field models. Surveys in Geophysics, 2020, 41(4): 767- 804.

DOI

Jiang T, Li J C, Dang Y M, et al. Regional gravity field modeling based on rectangular harmonic analysis. Science China Earth Sciences, 2014, 57(7): 1637- 1644.

DOI

Jiang T, Xiao X N, Dang Y M, et al. Evaluation of the external accord accuracy of airborne gravity data with upward continuation. Acta Geodaetica et Cartographica Sinica, 2018, 47(5): 567- 574.

Li J C. The recent Chinese terrestrial digital height datum model: gravimetric quasi-geoid CNGG2011. Acta Geodaetica et Cartographica Sinica, 2012, 41(5): 651- 660. 651-660, 669

Li X X. 2013. Building of an ultra-high-degree geopotential model [Master's thesis](in Chinese). Zhengzhou: PLA Information Engineering University.

Liang W, Xu X Y, Li J C, et al. 2023. Research on determination of gravity field models in regional area using airborne gravity data and reference model. Geomatics and Information Science of Wuhan University (in Chinese). https://doi.org/10.13203/j.whugis20230002.

Liu M, Huang M T, Ouyang Y Z, et al. Test and analysis of downward continuation models for airborne gravity data with regard to the effect of topographic height. Acta Geodaetica et Cartographica Sinica, 2016, 45(5): 521- 530.

Liu X G, Sun Z M, Guan B, et al. Downward continuation of airborne gravimetry data based on improved poisson integral iteration method. Acta Geodaetica et Cartographica Sinica, 2018, 47(9): 1188- 1195.

Lu Y, Xu H Z. The regional model of high-order gravitational field and the potential coefficient model of the Qinghai-Tibet region. Chinese Journal of Geophysics, 1994, 37(4): 487- 498.

Ma J, Zhal Z H, Feng C Q, et al. Optimization algorithm for analytical downward continuation based onimproved radial derivative. Acta Geodaetica et Cartographica Sinica, 2023, 52(10): 1650- 1660.

Ning J S, Luo Z C, Li J C. Present situations and technical modes for refining provincial (municipal) geoid in China. Journal of Geodesy and Geodynamics, 2004, 24(1): 4- 8.

Ouyang Y Z, Deng K L, Huang M T, et al. The choice of integral radius and its far-zone effects in the continuation of the gravity anomaly. Hydrographic Surveying and Charting, 2011, 31(4): 5- 7. 5-7, 12

Pavlis N K, Holmes S A, Kenyon S C, et al. The development and evaluation of the Earth Gravitational Model 2008 (EGM2008). Journal of Geophysical Research: Solid Earth, 2012, 117(B4): B04406

DOI

Rapp R H, Pavlis N K. The development and analysis of geopotential coefficient models to spherical harmonic degree 360. Journal of Geophysical Research: Solid Earth, 1990, 95(B13): 21885- 21911.

DOI

Rexer M, Hirt C. Ultra-high-degree surface spherical harmonic analysis using the gauss-legendre and the driscoll/healy quadrature theorem and application to planetary topography models of earth, mars and moon. Surveys in Geophysics, 2015, 36(6): 803- 830.

DOI

Sandwell D T, Müller R D, Smith W H F, et al. New global marine gravity model from CryoSat-2 and Jason-1 reveals buried tectonic structure. Science, 2014, 346(6205): 65- 67.

DOI

Sneeuw N. Global spherical harmonic analysis by least-squares and numerical quadrature methods in historical perspective. Geophysical Journal International, 1994, 118(3): 707- 716.

DOI

Thébault E, Schott J J, Mandea M. Revised spherical cap harmonic analysis (R-SCHA): validation and properties. Journal of Geophysical Research: Solid Earth, 2006, 111(B1): B01102

Tian J L, Li X X, Wu X P, et al. Fast realization of ultra-high-degree geopotential model by improved least-squares method. Acta Geodaetica et Cartographica Sinica, 2018, 47(11): 1437- 1445.

Wang K, Zhang C D, Wu X, et al. Construction and quas-stable analysis of topographic/isostatic potential coefficients of gravity field. Chinese Journal of Geophysics, 2014, 57(3): 760- 769.

DOI

Wang X T, Xia Z R, Shi P, et al. A comparison of different downward continuation methods for airborne gravity data. Chinese Journal of Geophysics, 2004, 47(6): 1017- 1022.

Wieczorek M A, Meschede M. SHTools: tools for working with spherical harmonics. Geochemistry, Geophysics, Geosystems, 2018, 19(8): 2574- 2592.

DOI

Wu W J. On the application of geomagnetic spherical harmonic analysis method in local areas. Chinese Science Bulletin, 1984, 29(16): 999- 1003.

DOI

Wu X P. Point-mass model of local gravity field. Acta Geodaetica et Cartographica Sinica, 1984, 13(4): 249- 258.

Xu H Z, Lu Y, Zhang K F. A tailored method for altimetry-gravimetry boundary value problems. Acta Geodaetica et Cartographica Sinica, 2002, 31(S1): 12- 15.

Xu W Y, Qu J M, Du A M. Geomagnetic field modelling for the globe and a limited region. Progress in Geophysics, 2011, 26(2): 398- 415.

DOI

Zhai Z H, Sun Z M. Continuation model construction and application analysis of local gravity field based on least square collocation. Chinese Journal of Geophysics, 2009, 52(7): 1700- 1706.

DOI

Zhao Y Q, Li J C, Xu X Y, et al. Determination of static gravity field model by using satellite data of GOCE andGRACE. Chinese Journal of Geophysics, 2023, 66(6): 2322- 2336.

DOI

Zhu C, Guo J Y, Gao J Y, et al. Marine gravity determined from multi-satellite GM/ERM altimeter data over the South China Sea: SCSGA V1. 0. Journal of Geodesy, 2020, 94(5): 50

DOI

Zingerle P, Pail R, Gruber T, et al. The combined global gravity field model XGM2019e. Journal of Geodesy, 2020, 94(7): 66

DOI

党 亚民, 蒋 涛, 杨 元喜, 等. 中国大地测量研究进展(2019—2023). 测绘学报, 2023, 52(9): 1419- 1436.

丁 子凡, 李 建成, 徐 新禹. 一种重力异常向下延拓的球谐分析迭代方法研究. 测绘科学, 2022, 47(10): 8- 14.

黄 谟涛, 刘 敏, 邓 凯亮, 等. 基于向上延拓的航空重力向下解析延拓解. 地球物理学报, 2018, 61(12): 4746- 4757.

DOI

黄 谟涛, 邓 凯亮, 吴 太旗, 等. 重力异常向上延拓严密改化模型及向下延拓应用. 测绘学报, 2022, 51(1): 41- 52.

蒋 涛, 肖 学年, 党 亚民, 等. 航空重力向上延拓的外部精度评估. 测绘学报, 2018, 47(5): 567- 574.

李 建成. 最新中国陆地数字高程基准模型: 重力似大地水准面CNGG2011. 测绘学报, 2012, 41(5): 651- 660. 651-660, 669

李新星. 2013. 超高阶地球重力场模型的构建[硕士论文]. 郑州: 解放军信息工程大学.

梁伟, 徐新禹, 李建成, 等. 2023. 联合先验模型和航空重力数据构建适用于局部区域的超高阶重力场模型. 武汉大学学报(信息科学版), https://doi.org/10.13203/j.whugis20230002.

刘 敏, 黄 谟涛, 欧阳 永忠, 等. 顾及地形效应的重力向下延拓模型分析与检验. 测绘学报, 2016, 45(5): 521- 530.

刘 晓刚, 孙 中苗, 管 斌, 等. 航空重力测量数据向下延拓的改进Poisson积分迭代法. 测绘学报, 2018, 47(9): 1188- 1195.

陆 洋, 许 厚泽. 区域高阶重力场模型与青藏地区局部位系数模型. 地球物理学报, 1994, 37(4): 487- 498.

马 健, 翟 振和, 冯 长强, 等. 基于改进径向导数的解析向下延拓优化算法. 测绘学报, 2023, 52(10): 1650- 1660.

宁 津生, 罗 志才, 李 建成. 我国省市级大地水准面精化的现状及技术模式. 大地测量与地球动力学, 2004, 24(1): 4- 8.

欧阳 永忠, 邓 凯亮, 黄 谟涛, 等. 重力异常延拓计算积分半径的选择及远区效应. 海洋测绘, 2011, 31(4): 5- 7. 5-7, 12

田 家磊, 李 新星, 吴 晓平, 等. 超高阶重力场模型最小二乘快速实现. 测绘学报, 2018, 47(11): 1437- 1445.

王 凯, 张 传定, 吴 星, 等. 地形/均衡重力场模型构制及其拟稳分析. 地球物理学报, 2014, 57(3): 760- 769.

DOI

王 兴涛, 夏 哲仁, 石 磐, 等. 航空重力测量数据向下延拓方法比较. 地球物理学报, 2004, 47(6): 1017- 1022.

吴 文京. 关于地磁球谐分析方法应用于局部地区的问题. 科学通报, 1984, 29(16): 999- 1003.

吴 晓平. 局部重力场的点质量模型. 测绘学报, 1984, 13(4): 249- 258.

许 厚泽, 陆 洋, 张 克非. 测高-重力边值问题的局部剪裁解. 测绘学报, 2002, 2002(S1): 12- 15.

徐 文耀, 区 加明, 杜 爱民. 地磁场全球建模和局域建模. 地球物理学进展, 2011, 26(2): 398- 415.

DOI

翟 振和, 孙 中苗. 基于配置法的局部重力场延拓模型构建与应用分析. 地球物理学报, 2009, 52(7): 1700- 1706.

DOI

赵 永奇, 李 建成, 徐 新禹, 等. 利用GOCE和GRACE卫星观测数据确定静态重力场模型. 地球物理学报, 2023, 66(6): 2322- 2336.

DOI

Options
Outlines

/

〈 〉