时变重力场反映地球系统内物质的分布、运动和变化，是研究地球系统内物质运动状态及其动力学机制的重要约束信息，可为解决人类面临的资源、环境和灾害等紧迫问题提供重要的基础地球物理信息。近十几年来，GRACE(Gravity Recovery And Climate Experiment)卫星已成功获取全球均匀覆盖的时变重力观测数据，使地球时变重力场模型在分辨率和精度上实现了重大提升。GRACE月时变重力场在大空间尺度上定量揭示了全球环境变化，例如：陆地水储量变化、冰川消融、海洋变化、冰后回弹以及大地震等导致的地球表层质量重新分布。

GRACE重力卫星具有可重复、大范围的观测优势，但仅对重力场信号的中长波长部分较为敏感，对短波长部分敏感性较弱，不利于恢复较高阶的时变地球重力场模型。GRACE的实际空间分辨率约400 km，其数据存在空间域的南北向条带等噪声需做滤波处理，会进一步降低空间分辨率，难以满足环境保护、资源监测以及灾害预测等研究对高精度高分辨率时变重力场的需求。

地面重力时变观测在地表进行，通过定期同点位复测来获得重力场随时间的变化。与GRACE重力卫星相比，地面重力测量的优点在于重力信号衰减幅度小，观测结果的短波长部分更加精确，空间分辨率也更高，适于分析区域性、近地表的物质运动和迁移现象。但地面重力观测易受地理和环境等因素制约：一方面，很多区域（如南极、沙漠和冰川等）难以施测，存在数据空白区，无法像GRACE卫星观测那样实现数据全球覆盖；另一方面，在可施测区域内，测点空间位置也难以实现均匀分布。直接利用地面重力观测数据确定区域重力场模型时，由于数据仅区域覆盖且空间不均匀分布，无法满足球谐分析这样的经典重力场恢复理论对数据的要求（基函数在区域内部不满足正交性）。因此，需借助局部重力场恢复理论和方法。

目前应用较为广泛的局部重力场恢复方法包括最小二乘配置法、球冠谐分析法、Mascon法和特殊基函数法等。Slepian局部谱分析方法即为一种特殊的基函数法，最早由D. Slepian提出，旨在解决一维连续情况下的时域和频域能量集中问题，随后被不断发展并引入到重力场研究领域中(Simons and Dahlen, 2006; Simons et al., 2006; Simons, 2010)。Slepian方法通过构建在研究区域内部的正交基函数，可使信号能量集中在研究区域内部，适用于表达区域尺度的重力场变化，并且Slepian展开与相同阶次的经典球谐展开具有同样的意义(Simons, 2010; 陈石等, 2017)。

本文对利用Slepian局部谱分析恢复局部重力场的方法进行了研究，并将Slepian方法应用到我国华北地区，确定了华北2011–2013年尺度的120阶局部重力场时变模型，并与研究区域内GRACE卫星重力估计结果进行对比分析。

Slepian方法的基本原理及恢复局部重力场模型的流程详见韩建成等(2021)。Slepian方法是一种具有稀疏(sparse)特性的时空局部谱方法。图1直观反映了华北地区Slepian函数分别展开到90、120和180阶时聚集因子的稀疏分布情况。对于120阶Slepian展开，聚集因子总数量为14641个，其中大于0.001的共41个，仅占总数量的0.3%，剩余均接近0，分布具有明显的稀疏性。

图1 华北地区对应90、120和180阶Slepian展开的聚集因子分布

图2进一步给出了华北地区Slepian展开到90、120和180阶时基函数的空间分布情况，其中洋红色曲线标记重力数据覆盖范围的边界。由图2可以看出，对于每组展开，当聚集因子较接近1时，所对应基函数基本都集中在研究区域内部，能量很少泄漏到区域外部；当聚集因子很小时（如图2第3行中接近0.001的情况），所对应基函数基本集中在研究区域外部，在研究区域内部几乎没有贡献。

图2 华北地区3组Slepian展开的基函数空间分布，分别对应90阶（左边2列）、120阶（中间2列）和180阶（右边2列）展开. 对于每组展开，第1和第2行为聚集因子最大的前4个基函数，第3行为聚集因子接近0.001的基函数. 洋红色线为研究区域边界，黑色线为陆海边界

本文采用了华北地区重力测网2011–2013年间的实际重力观测成果，一共6期观测数据（每年上下半年各测1期）。重力网平差采用基于贝叶斯平差模型的pyBACGS软件完成(Chen et al., 2019)，该方法可以较好地解决多台相对重力仪器联合平差的权系数优化、非线性漂移估计和仪器格值的最优化估计等问题，6期观测平差结果的点值精度均优于10×10-8 m/s2. 图3给出了6期数据覆盖的重复测点位置的空间分布，共有364个。

图3 本文研究区域及重复测点位置分布. 洋红色线为研究区域边界，绿色十字代表重复测点位置，灰色线为省界，黑色线为陆海边界

在研究区域内部364个重复测点处，将相邻两年对应测期内的重力观测做差（后一年观测减去前一年）获得1年尺度重力变化。所有比较都限定在相对应的测期内，以减弱水文变化等季节性信号的影响。最终得到4组1年尺度重力变化，即基于第1期观测的2013.04–2012.03、2012.03–2011.03（“–”号表示前后两组数据做差），以及基于第2期观测的2013.08–2012.08和2012.08–2011.08. 获得研究区域重复测点处的重力变化后，接下来应用Slepian方法恢复华北地区120阶时变重力场模型（最佳截断数取26），空间变化结果（格网间距为0.1°×0.1°）如图4所示。最大展开阶数以及Slepian最佳截断数的确定，详见韩建成等(2021)。

除1年尺度时变重力场模型外，本文还解算了2组华北地区120阶的2年尺度重力变化模型，第1组利用基于上半年观测(2013.04–2011.03)的重力变化数据解算，第2组利用基于下半年观测(2013.08–2011.08)的重力变化数据解算，解算结果如图5所示。

图 4 华北地区由Slepian方法确定的120阶1年尺度重力场变化 (2011–2013) 空间分布：(a)-(c) 基于相邻两年第1期观测差值; (b)-(d) 基于相邻两年第2期观测差值. 洋红色线为研究区域边界，灰色线为省界，黑色线为陆海边界. 单位为10-8 m/s2

图 5 华北地区由Slepian方法确定的120阶2年尺度重力场变化空间分布：(a) 基于2013.04和2011.03的观测差值; (b) 基于2013.08和2011.08的观测差值. 洋红色线为研究区域边界，灰色线为省界，黑色线为陆海边界. 单位为10-8 m/s2

本文还利用Slepian方法计算了华北地区2011–2013年间2年尺度的GRACE重力变化（2013.04–2011.03和2013.07–2011.08），如图6所示。

图6华北地区由Slepian方法得到的GRACE重力变化空间分布，其中(a)、(b)分别为2013.04–2011.03和2013.07–2011.08结果. 洋红色线为研究区域边界，灰色线为省界，黑色线为陆海边界. 单位为10-8 m/s2

对于2年尺度的重力变化，地面结果（图5）与GRACE结果（图6）表现出较为一致的“南减北增”空间分布特征。进一步对比图5和图6，可发现2年尺度的地面与GRACE卫星结果差异最大区域出现在石家庄–衡水连线一带，2011–2013年区域内地面结果为增大趋势，GRACE结果为减少趋势。该区域长期存在严重的浅层、深层地下水超量开采以及地表沉降问题，这是影响区域内重力场长期变化的两个主要因素。对地面测量而言，地下水亏损会令重力观测减小，地表沉降会导致重力观测增加。根据Shen等(2015)的地面观测结果，2009–2013年在河北省南部（包含石家庄–衡水连线一带）由地表沉降导致的重力增加效应要大于地下水亏损带来的减小效应；由于2009–2013年间石家庄–衡水连线一带地表沉降速度基本维持不变(约-0.05 m/a)，我们推断2011–2013年区域内地表沉降的影响仍然大于地下水亏损的影响，因此地面重力结果呈增大趋势。对GRACE观测而言，在卫星轨道处地面沉降导致的重力变化已大幅衰减，此时地下水亏损导致的重力减小效应占主导地位，因此GRACE结果整体为减少趋势。

本文的创新及取得的成果如下：(1)              首次尝试利用Slepian局部谱分析方法和多期地面重力观测确定更高空间分辨率的时变重力场模型，讨论给出了相应构建流程。(2)              提出了Slepian区域重力场建模参数的优化方法。(3)              获得了2011–2013年华北地区120阶年尺度时变重力场模型，并验证了模型的可靠性。

本文研究结果可为区域重力场建模提供新的参考，也可为华北地区水资源变化监测、构造活动分析以及地震风险性评估等研究提供高分辨率的时变重力场模型支撑。

研究成果发表于“地球物理学报”。研究受国家自然科学基金(41974095)，中国博士后科学基金(2018M641424)，中国地震局地球物理研究所基本科研业务费专项(DQJB19B03,DQJB20X09)联合资助。

【相关文献列表】

[1].      Chen S, Zhuang J, Li X et al. 2019. Bayesian approach for network adjustment for gravity survey campaign: methodology and model test. Journal of Geodesy, 93:681-700, doi:10.1007/s00190-018-1190-7.

[2].      Shen C, Xuan S, Zou Z et al. 2015. Trends in gravity changes from 2009 to 2013 derived from ground-based gravimetry and GRACE data in North China. Geodesy and Geodynamics, 6(6):423-428, doi:10.1016/j.geog.2015.08.001.

[3].      Simons F J. 2010. Slepian functions and their use in signal estimation and spectral analysis, in: Freeden W, Nashed M Z, Sonar T (Eds.), Handbook of Geomathematics. Springer, Heidelberg, Germany, pp. 891-923.

[4].      Simons F J, Dahlen F A. 2006. Spherical Slepian functions and the polar gap in geodesy. Geophysical Journal International, 166(3):1039-1061, doi:10.1111/j.1365-246X.2006.03065.x.

[5].      Simons F J, Dahlen F A, Wieczorek M A. 2006. Spatiospectral Concentration on a Sphere. SIAM Review, 48(3):504-536, doi:10.1137/s0036144504445765.

[6].      陈石, 徐伟民, 王谦身. 2017. 应用Slepian局部谱方法解算中国大陆重力场球谐模型. 测绘学报, 46(8):952-960, doi:10.11947/j.AGCS.2017.20150542.

[7].      韩建成, 陈石, 卢红艳, 徐伟民. 2021. 基于Slepian 方法和地面重力观测确定时变重力场模型：以2011–2013 年华北地区数据为例. 地球物理学报, 64(5):1542-1557, doi: 10.6038/cjg2021O0240.

【作者简介】

韩建成2012年毕业于武汉大学测绘学院，获固体地球物理学博士学位。2014年至2017年于台湾交通大学从事博士后研究，2018年进入中国地震局地球物理研究所从事博士后研究，现为第四研究室助理研究员。目前主要研究方向为重力场恢复、时变重力场应用及高性能计算等。