基于重力卫星数据探测核幔边界动态过程信号

摘要

使用GRACE（Gravity Recovery and Climate Experiment）重力场数据，联合多源数据模型剔除地表过程的重力信号，以分离得到核幔边界处包含液态外核流体运动特征的重力信号，同时根据CHAOS-7地磁场模型重建核幔边界与地球表面的地核地磁场变化，以此为基础使用主成分分析（PCA）与独立主成分分析（ICA）方法对重力信号与磁场信号进行统计分解，得到2种信号的时变模态与空间模态，对比分析2种信号分解结果的相关性，研究表明，重力信号第3模式的时变模态与磁场信号第2模式的时变模态相关性较好，同时，其空间模态反映了太平洋低纬度地区较为明显的东西向“偶极子”形式的变化特征，与地磁场二阶导数在太平洋地区的空间变化特征相对应，在2014年以后重力信号的时变模态与磁场信号二阶导数的时变模态均表现出相同的线性快速变化特征，与现有的研究成果相符合。

关键词

大地测量学; GRACE任务; 时变重力场; 核幔边界; 质量迁移

直到20年前，虽然人们对地球重力场的“静态”部分了解相对准确，但对其时间变化的了解仅限于较低的球谐阶次。2002年，随着GRACE（Gravity Recovery and Climate Experiment）任务的启动，情况得到了大幅改善；2018年，其后续任务GRACE-FO（Gravity Recovery and Climate Experiment-Follow On）成功实施，重力卫星任务使人们能够以前所未有的时空分辨率和精度测量地球重力随时间的变化。

在年际和年代际时间尺度上，陆地水储量和海洋质量的变化似乎主导了观测到的行星尺度重力变化。然而对于季节项的周期变化，观测值和预测值拟合得并不理想。其原因是这些地表过程的模型仍然不完善，尤其是与陆地蓄水量有关的模型。同时，地球液核中的核心流在年际到10年的时间尺度上也能影响地球重力场变化，部分无法解释的信号也可能与地核流动有关^［1］。地球流体核心的动力学过程导致的质量再分配主要以下述几种不同的方式引起重力变化，如密度异常导致流体外核表层对流、分层密度层中的垂直运动产生密度迁移、由随时间变化的核心流引起的压力变化推动和拉动核幔边界（core-mantle boundary， CMB）^［2-3］。

最近的一些研究结果表明，观测到的特定大地测量信号可能与地核动力学过程有关。首先，Mandea等^［4］证实，在对已知地表贡献进行校正后，卫星记录的时变重力场主要时空模态与磁场垂直分量的二阶主要时空模态密切相关。研究涵盖了2002年年中至2010年年中之间的有限时间跨度，在几个地理位置，特别是在太平洋区域，剩余重力信号与磁场径向分量的长期加速度（径向分量的二阶导数）之间存在良好的相关性。此外，还对全球重力场和全球磁场径向长期加速度时间序列进行了奇异值分解（singular value decomposition， SVD）。虽然每个时间序列的主导模式具有不同的空间结构，但它们的年际变化具有良好的相关性^［5］。同时后续研究^［4-5］表明，良好的相关性一直持续到2015年底。文献［6-7］研究探测到了全球定位系统（GPS）网络记录的地表位移中行星级尺度约为1mm的振幅信号，每个周期约为6年^［6-7］。重力观测和极移观测^［8］中也观察到周期为6年的大地测量信号。在文献［6-9］中，研究者认为，地核中的动力学过程是观测信号的来源，主要是因为核心流携带的轴向角动量与6年期间的日长变化之间存在联系，这说明通过大地测量观测可以捕捉地核内的动力学过程^［10］。

本文基于反演得到的GRACE重力场数据，利用多源数据模型在扣除各种非地核信号后，得到包含有地核信号的重力残项序列，同时根据CHAOS-7地磁场模型，构建出地磁场径向分量、分量一阶导数、分量二阶导数的全球模型，并使用统计分解方法对重力场信号与磁场信号进行处理，得到重力信号与磁场信号的空间模态与时变模态，发现重力信号第3模式的时变特征与磁场信号第2模式的时变信号相关性较好，在2014年以后，两者均存在线性快速变化的特征，同时，其空间模态反映了太平洋低纬度地区较为明显东西向“偶极子”形式的空间特征，与地磁场二阶导数在太平洋地区的空间变化特征相对应。

1 实验数据

1.1　GRACE重力场数据

GRACE重力场数据可以公开获得，主要由德克萨斯大学空间研究中心（CSR）、德国地学中心（GFZ）、美国喷气推进实验室（JPL）三大机构发布。本文使用的是CSR GRACE/GRACE-FO RL06 Solutions数据，由CSR官网解算提供，其中，潮汐影响包括海洋潮汐、固体潮和极潮，都已在数据处理中去除，非潮汐的大气和海洋的影响也同时在数据处理中剔除。选用CSR机构提供的阶数为60的模型，将官网发布的GRACE球谐系数解转换为网格形式的质量异常（详细的计算方法及有关处理见GRACE CSR用户手册）。应用的后处理包括：①添加对地球中心运动和地球扁率的独立估计（C20改正）；②过滤以南北条纹为特征的相关误差；③与地心有关的一阶项改正；④采用ICE6G-D模型进行冰川均衡调整（GIA）的更正。

1.2　GLDAS模型

全球陆地数据同化系统（GLDAS）模型由美国宇航局协同美国国家环境预报中心进行发布。该模型使用的是卫星观测数据和地面观测数据，通过先进的陆地表面模型和数据同化技术可以得到地表流体的变化，并输出陆地地表的各项物理参数。这一强大的新陆地表面建模系统集成了来自先进观测系统的数据，用以支持改进水文气象的研究与调查^［11］。

根据GLDAS-NOAH模型，将总体陆地水（土壤湿度、根系水、植被水、径流、冠层水、地表积雪等）储量变化统一归化成与GRACE数据相同的等效水高形式，并用到GRACE数据的观测结果中，以扣除GRACE信号中陆地水储量的成分，该模型时间分辨率为1个月，并以NetCDF4的格式按1°$\times$1°的单元网格形式发布。

1.3　ECCO模型

海洋环流和气候估算（ECCO）联合会成立于1999年，主要对海洋环流及其气候领域进行详尽的研究。如今，在NASA物理海洋学、建模构建冰冻圈项目的持续支持下，ECCO正在提高海洋状态的估计范围和估计真实度。与通常的数据同化（如，在数值气候预报中使用的数据同化）不同，ECCO提供了海洋随时间演化的物理性特征一致性估计。ECCO的长期目标是综合考虑海洋环流、海冰、生物化学、大气等因素提供全球高分辨率的海洋状态估计模型^［12］。

使用ECCO Version 4、Release 4产品数据，ECCO V4R4综合了总环流模型（MITgcm）和大多数可用的卫星数据及现场实测数据，以产生一个物理性一致的水循环闭合海洋模型，选取ECCO产品中海表面高度数据（SSH）以用于扣除海洋质量变化信号。

1.4　WOA模型

世界海洋地图集2018（WOA 2018）由美国国家海洋和大气管理局气候与全球变化计划的拨款促成，该计划由国家海洋数据中心（现为国家环境信息中心）统筹与管理，其目的是编制研究海洋学数据库和基于这些数据库进行诊断研究。该地图集所包含的数据收集于2018年世界海洋数据库（WOD 2018）中，统计和客观分析海洋的实地测量数据文件由国家环境信息中心（NCEI）在线分发。

实验采用的《2018年世界海洋地图集》（WOA 2018）于2019年7月发布，其对早期版本的《世界海洋地图集》进行了数据更新，包括自上一次发布以来向世界海洋数据库（WOD）添加的约300万个新海洋学模型以及更新了的质量约束条件。发布的海洋物理参量有：温度、盐度、含氧量、硝酸盐含量、磷酸盐含量、硅酸盐含量、电导率、密度。

1.5　ICE6G-D模型

目前国际上GIA全球模型有：ICE-3G、ICE-4G（VM2）、ICE-5G（VM2）以及ICE-6G_C（VM5a）等ICE模型和Paulson07模型和Geruo13模型；GIA区域模型主要针对南极地区，例如IJ05模型、IJ05_R2模型、W12和W12a模型。

多伦多大学根据ICE-NG（VMX）序列发布了冰川均衡调整模型，称为ICE-6G_C（VM5a），最初以球谐系数的形式发布，阶数为256，并已被证明可为来自北美、欧亚大陆和南极洲的GPS观测数据、GRACE卫星反演得到的时变重力场数据以及放射性碳测定的相对海平面历史数据提供精确的数据拟合，Purcell等^［13］认为ICE-6G_C模型存在改进空间，提出了ICE-6G_D模型，在一定程度上优化了模型结果，但仍然存在一定物理层面上的差异。

本文选用JPL发布的ICE-6G_D模型进行GIA改正，以球谐系数的形式发布，阶数为256阶。

1.6　CHAOS-7模型

丹麦大学国家空间研究所基于低地球轨道卫星Swarm、CryoSat-2、CHAMP、SAC-C和Örsted收集的磁场观测以及地面观测站月平均测量值的年度差异发布了1999年至2020年期间与时间相关的近地地磁场CHAOS-7模型。CHAOS-7模型包括一个球谐度为20的时变地核场、一个与LCS-1岩石圈场模型合并的25阶以上的静态内部场、磁层场及其对应感应场的模型，描述卫星矢量磁强计对准的欧拉角估计数据以及CryoSat-2的磁强计校准参数。只有来自磁静区域的数据满足严格的地磁平静时间标准才用于磁场估计，使用的模型参数解算方法为迭代重加权正则化最小二乘法，与以前版本的CHAOS模型相比，在高阶系数中，时间相关的地核内源场的正则化条件更为宽松。

CHAOS-7是过去15年DTU开发的一系列时间相关地磁场模型中的最新模型^［14-16］。在过去的6年中，Swarm卫星三星组提供的卫星数据达到了全球完整和同质的数据覆盖。CHAOS-7的数据来源除了基于磁测卫星上收集的磁强计数据，还有经修订的月平均地面观测站长期变化数据系列的更新版^［15］。本次实验所提及的地磁场均为地核地磁场。

2 理论原理

2.1　数据处理流程

主要思路是利用多源数据模型对预处理完成后的GRACE重力信号进行地表浅层物质迁移的改正，以此来突出核幔边界的重力信号。在对基于GRACE网格的重力时间序列进行了所有影响（陆地蓄水、陆地冰和海洋的影响）的校正后，将PCA-ICA方法使用于剩余GRACE重力场以及地磁场模型，以研究其时空特异性。其主要的年际重力信号与磁场信号相关，含有发生在核幔边界（CMB）的过程特征。流程示意图如图1。

图1  实验流程示意

Fig. 1  Flowchart of experiment

2.2　GRACE时变重力场反演

基于不同的反演数据，主要有两大类反演方法：第一类，基于GRACE的Level-2月地球重力场的球谐系数反演，这种方法首次由Wahr等^［17］提出，是目前较为常用的方法，称为位理论法；第二类，根据GRACE的Level-1B数据进行反演，如质量瘤法（mascon）法，该方法在一定程度上可以提高反演地表质量变化的时空分辨率，但方法的技术原理较为复杂，处理流程难度较大。

在上述2种方法的基础上，发展出新的方法，比如球谐系数拟合法、点质量模型法和最优先验模型法等。本文采用最为常用的位理论法。

$$\mathrm{\Delta }h\left(\theta ,\lambda \right)=\frac{a{\rho }_{\mathrm{a}\mathrm{v}\mathrm{e}}}{\mathrm{3}{\rho }_{\mathrm{w}}}{\displaystyle \sum _{l=\mathrm{0}}^{\mathrm{\infty }}}{\displaystyle \sum _{m=\mathrm{0}}^l}{\overline{P}}_{lm}(\mathrm{c}\mathrm{o}\mathrm{s}\theta )\frac{\mathrm{2}l+\mathrm{1}}{\mathrm{1}+k_l}·$$

$$(\mathrm{\Delta }{\overline{C}}_{lm}\mathrm{c}\mathrm{o}\mathrm{s}(m\lambda )+\mathrm{\Delta }{\overline{S}}_{lm}\mathrm{s}\mathrm{i}\mathrm{n}(m\lambda ))$$

（1）

式中：$\mathrm{\Delta }h\left(\theta ,\lambda \right)$为等效水高变化，反映地球地表水的质量变化，其中$\theta \mathrm{为余}\mathrm{纬},\lambda \mathrm{为经}\mathrm{度}$；$a$为地球平均参考半径；${\rho }_{\mathrm{a}\mathrm{v}\mathrm{e}}$为地球平均密度；${\rho }_{\mathrm{w}}$为水的密度；${\overline{P}}_{lm}(\mathrm{c}\mathrm{o}\mathrm{s}\theta )$为完全规格化的缔合勒让德函数，其中l为阶，m为次；$k_l\mathrm{为}l$阶的负荷勒夫数；$\mathrm{\Delta }{\overline{C}}_{lm}$和$\mathrm{\Delta }{\overline{S}}_{lm}$为完全规格化的球谐系数变化，分别是$\mathrm{c}\mathrm{o}\mathrm{s}(m\lambda )$的系数、$\mathrm{s}\mathrm{i}\mathrm{n}(m\lambda )$的系数。式（1）即为用GRACE数据反演水储量变化的基本公式。

2.3　地磁场模型数据

CHAOS-7模型的基本参数化与CHAOS系列早期版本的参数化相同，假设测量是在没有电流的区域中收集的，因此，在电磁学的准静态近似下，矢量磁场$\mathit{B}$可以用标量势表示，使得$\mathit{B}=-\nabla V$。磁标量势$V=V^{\mathrm{i}\mathrm{n}\mathrm{t}}+V^{\mathrm{e}\mathrm{x}\mathrm{t}}$，由内部（主要是核心和岩石圈）源的势$V^{\mathrm{i}\mathrm{n}\mathrm{t}}$、外部（此处假设为磁层）源的势$V^{\mathrm{e}\mathrm{x}\mathrm{t}}$组成。内部和外部都以球面谐波（SH）形式展开。

对于内源场，在地心地固系下（ECEF），对其进行球谐展开^［16］，有

$$V^{\mathrm{i}\mathrm{n}\mathrm{t}}(R,\theta ,\lambda )=a{\displaystyle \sum _{l=\mathrm{1}}^{N_{\mathrm{i}\mathrm{n}\mathrm{t}}}}{\displaystyle \sum _{m=\mathrm{0}}^l}\left(g_{lm}\mathrm{c}\mathrm{o}\mathrm{s}m\lambda +h_{lm}\mathrm{s}\mathrm{i}\mathrm{n}m\lambda \right)$$$$·$$$${\left(\frac{a}{R}\right)}^{l+\mathrm{1}}{P}_{lm}\left(\mathrm{c}\mathrm{o}\mathrm{s}\theta \right)$$

（2）

式中：$(R,\theta ,\lambda )$为球面坐标；$N_{\mathrm{i}\mathrm{n}\mathrm{t}}\mathrm{为描}\mathrm{述内}\mathrm{源场}\mathrm{的最}\mathrm{大阶}\mathrm{数};g_{lm}$与$h_{lm}$为描述内源场的高斯系数，在CHOAS-7模型中，1~20阶的内源场高斯系数$g_{lm}、h_{lm}$与时间有关，21阶及以上的高斯系数是静态的；$P_{lm}$为施密特半归一化缔合勒让德函数。

2.4　基于统计信息的信号分解

2.4.1　主成分分析（PCA）

将经验正交函数分析方法（EOF）应用于分解重力残差项，用来研究改正后的GRACE重力场剩余信号和时空特性，该方法也称主成分分析（PCA），是一种分析矩阵数据中结构特征、提取主要数据特征量的统计分解方法^［18］。将一组GRACE单元格网点时间序列视为一组观测值时间序列数据矩阵，对于整个格网的第$i$个格网点进行$m$次观测，则对应该格网点的观测向量${\mathit{X}}_i$可表示为

$${\mathit{X}}_i=\left(\begin{array}{c}{\mathit{x}}_{\mathrm{1},i}\\ {\mathit{x}}_{\mathrm{2},i}\\ \begin{array}{l}⋮\\ {\mathit{x}}_{m,i}\end{array}\end{array}\right)$$

（3）

由于实验地区有$n$个格网点，且对每个格网点都进行过$m$次观测，所以所有格网点上的距平观测值时间序列按列排列可以用一个维度为$s\times n$的矩阵$\mathit{Y}$表示，为

$$\mathit{Y}=({\mathit{X}}_{\mathrm{1}},{\mathit{X}}_{\mathrm{2}},\cdots ,{\mathit{X}}_n)=\left(\begin{array}{cc}\begin{array}{cc}{x}_{\mathrm{1,1}}& x_{\mathrm{1,2}}\\ x_{\mathrm{2,1}}& x_{\mathrm{2,2}}\end{array}& \begin{array}{cc}\cdots & x_{\mathrm{1},n}\\ \cdots & x_{\mathrm{2},n}\end{array}\\ \begin{array}{cc}⋮& ⋮\\ x_{s,\mathrm{1}}& x_{s,\mathrm{2}}\end{array}& \begin{array}{cc}& ⋮\\ \dots & x_{s,n}\end{array}\end{array}\right)$$

（4）

$\mathrm{式中}:s$代表观测时间序列的长度。为了凸显变化信息，式（4）对矩阵$\mathit{Y}$通常进行距平处理，主成分分解的结果可以写为

$$\mathit{Y}=\mathit{P}{\mathit{E}}^{\mathrm{T}}$$

（5）

式中：P为时变模态矩阵；E为空间模态矩阵。矩阵$\mathit{E}$中的$s$列数据代表观测值$\mathit{Y}$的$s$个空间模态，即经验正交函数，观测值矩阵$\mathit{Y}$在第$i$个对应空间模态上的投影${\mathit{p}}_i=\mathit{Y}{\mathit{e}}_i$就是该空间模态对应的时变模态，也称为主成分，${\mathit{p}}_i$是矩阵P的第i列，${\mathit{e}}_i$是矩阵$\mathit{E}$的第i列。

式（5）右侧的矩阵可以采用奇异值分解（SVD）方法进行矩阵的统计信号分解，以研究在应用所有校正后重力残项仍然存在的信号及其时空变异性。

基于SVD分解可以将式（5）分解成为

$$\mathit{Y}=\widehat{\mathit{P}}\mathit{\Lambda }{\mathit{E}}^{\mathrm{T}}$$

（6）

$\mathrm{式中}:\widehat{\mathit{P}}$包含了正交化后的时变模态，即${\widehat{\mathit{P}}}^{\mathrm{T}}\widehat{\mathit{P}}=\mathit{I};\mathit{\Lambda }$为包含按照从大到小顺序排列的奇异值对角矩阵。

对于分解后保留的主要的$r<\mathrm{m}\mathrm{i}\mathrm{n}(s,n)$个空间模态和对应的$r$个时变模态，观测值矩阵可以通过式（7）进行近似恢复：

$$\mathit{Y}\approx {\mathit{Y}}_r={\widehat{\mathit{P}}}_r{\mathit{\Lambda }}_r{\mathit{E}}_r^{\mathrm{T}}$$

（7）

式中：${\mathit{Y}}_r$为在保留前r个时变模态的情况下矩阵Y的近似矩阵；${\widehat{\mathit{P}}}_r$、${\mathit{\Lambda }}_r$和${\mathit{E}}_r$的维度分别为$s\times r$、$r\times r$和$n\times r$，且${\widehat{\mathit{P}}}_r^{\mathrm{T}}{\widehat{\mathit{P}}}_r={\mathit{I}}_r$，${\mathit{E}}_r^{\mathrm{T}}{\mathit{E}}_r={\mathit{I}}_r$，每一个${\widehat{\mathit{p}}}_j{\lambda }_j{\mathit{e}}_j^{\mathrm{T}}(j=\mathrm{1,2},\mathrm{3},\cdots ,r$）被称为观测值变量矩阵的“模式”，其中，${\widehat{\mathit{p}}}_j$代表${\widehat{\mathit{P}}}_j$中的第$j$列数据，${\lambda }_j$代表${\mathit{\Lambda }}_i$对角线上第$j$个值，${\mathit{e}}_j$代表${\mathit{E}}_r$中第$j$列数据。

2.4.2　独立主成分分析（ICA）

PCA方法是一种基于二阶统计信息的分解方法，在分解过程中，不管是基于时间序列构造的自协方差矩阵还是相关性矩阵，都只包含观测信号中的二阶统计信息，而忽略了观测值概率密度函数的高阶矩中所含的信号，且经过PCA分解得到的不同模式仅表现为彼此不相关，而并非相互独立。

基于上述讨论，需要对传统的PCA方法进行改进扩展，为了从基于观测数据的概率密度函数中获取更多的信息，就需要在分解过程中顾及概率密度函数的高阶统计信息作为影响因子，由此有关学者便提出独立主成分分析（ICA）方法。ICA方法的基本思想是利用高阶统计信息，使得分解得到的成分尽可能彼此独立，旨在在没有先验信息的条件下，利用高阶统计信息，将观测信号分解成多个彼此相互独立的成分，可以视为是基于二阶统计信息分解方法的一种拓展延伸。

将统计信息的信号分解问题一般化，转化为寻找将观测值矩阵$\mathit{Y}$转化为源信号矩阵$\mathit{S}$的函数$f\mathrm{。}$

$$\mathit{S}=f(\mathit{Y})$$

（8）

在多数应用场景下，式（8）转化关系为了简化计算，可以简单定义为线性转换，因此，式（8）的函数变换可以用矩阵乘积的形式表示，为

$$\mathit{S}=\mathit{W}\mathit{Y}$$

（9）

$\mathrm{式中}:\mathit{W}$为待求的解混矩阵。

目前常用的确定解混矩阵$\mathit{W}$的实用计算方法也是SVD分解方法，如式（10）：

$$\mathit{Y}=\widehat{\mathit{P}}\mathit{\Lambda }{\mathit{E}}^{\mathrm{T}}$$

（10）

ICA分解可以看作是2.4.1节中所讲到的PCA分解的旋转扩展应用，通过上述分析，基于正交模式旋转的ICA分解可表示为

$$\mathit{Y}\approx {\mathit{Y}}_r={\widehat{\mathit{P}}}_r{\mathit{\Lambda }}_r{\mathit{R}}_r{\mathit{R}}_r^{\mathrm{T}}{\mathit{E}}_r^{\mathrm{T}}$$

（11）

式中：${\mathit{R}}_r$为维度是$r\times r$的矩阵，为了得到ICA分解的模态，便需要确定此旋转矩阵，该矩阵对${\mathit{E}}_r$和${\widehat{\mathit{P}}}_r$进行矩阵旋转时使它们的列之间尽可能在统计上保持独立，并且其旋转矩阵本身也具有正交性，即${\mathit{R}}_r{\mathit{R}}_r^{\mathrm{T}}={\mathit{I}}_r$。那么，根据式（11）得到统计独立源信号可表示为

$${\overline{\mathit{S}}}_{r\times n}={\mathit{R}}_r^{\mathrm{T}}{\mathit{\Lambda }}_r^{-\mathrm{1}}{\widehat{\mathit{P}}}_r^{\mathrm{T}}{\mathit{Y}}_r=\mathit{W}{\mathit{Y}}_r$$

（12）

基于式（12）得到的源信号${\overline{\mathit{S}}}_{r\times n}$的行在空间上尽可能达到独立。以上分解称为ICA分解。

3 实验结果

3.1　重力场模型数据处理

为统一各个模型的时间跨度，分析2002年4月至2015年12月GRACE重力场的月时变序列，并对重力场模型进行相应的校正，提取与地核物理过程有关的深部地球重力信号，剔除的成分包括大气负荷、陆地流域水文变化和海洋环流变化引起的海洋质量变化以及GIA。

3.1.1　GRACE时变重力场

使用CSR机构提供的GRACE重力场模型，该数据可以公开获得，并根据2.2节提及的相关原理，将其数据转化为表征网格质量异常的等效水高形式，由于GRACE的卫星轨道形状对重力场的低阶项（$C_{\mathrm{20}}$项）不敏感，其低阶带谐项$C_{\mathrm{20}}$的精度相对较低，所以以卫星激光测距技术（SLR）测得的$C_{\mathrm{20}}$项代替GRACE的$C_{\mathrm{20}}$项。得到的结果如图2所示。

图2  GRACE全球时变重力场等效水高（2002年4月）

Fig. 2  GRACE global time-varying gravity field equivalent water height(April 2022)

3.1.2　大气负荷与海洋压力

根据CSR发布的GRACE用户手册可知，对于全球大气和海洋负载质量变化所引起的重力场变化，在Level-1B阶段，根据AOD1B（Atmosphere and Ocean De-Aliasing Level-1B）Release 6数据进行校正。

3.1.3　陆地水文成分

选用GLDAS水文模型来剔除GRACE重力场中的陆地水文成分，综合总体陆地水储量变化异常（土壤湿度、根系水、植被水、径流、冠层水、地表积雪等），将其统一归化成与GRACE数据相同的等效水高形式，调整空间分辨率并选取对应的时间跨度用到GRACE数据的观测结果中。

3.1.4　海洋质量变化

实验使用ECCO Version 4、Release 4产品中发布的海洋表面高（SSH）数据，得到全球测高的总体海平面变化，联合WOA18数据集的盐度分层模型与温度分层模型计算得到比容海平面变化，同时顾及冰后回弹对海面高的影响，计算得到由于全球海水质量变化引起的等效水高^［19］。

3.1.5　GIA影响

实验选用JPL发布的ICE-6G_D模型进行GIA修正，将球谐系数统一转化为等效水高的形式。

3.2　地磁场模型数据构建

丹麦大学国家空间研究所使用CHAOS-7模型研究了地磁场的时间变化特征。本文基于CHAOS-7模型构建地球表面以及核幔边界处地磁场径向分量强度和一阶导数、二阶导数模型。图3—5中，F、F '、F″分别表示地磁场径向强度、地磁场径向强度关于时间的一阶导数、地磁场径向强度关于时间的二阶导数。

图3  全球表面地磁场径向分量强度（2019年1月，截断阶数20阶）

Fig. 3  Radial magnetic field on surface of the earth (Jan. 2019, up to SH degree 20)

图4  全球表面地磁场径向分量强度一阶导数（2019年1月，截断阶数20阶）

Fig. 4  First time derivative of radial magnetic field on surface of the earth (Jan. 2019, up to SH degree 20)

图5  全球表面地磁场径向分量强度二阶导数（2019年1月，截断阶数20阶）

Fig. 5  Second time derivative of radial magnetic field on surface of the earth (Jan. 2019, up to SH degree 20)

图3—5显示了时间历元为2019年1月的CHAOS-7模型在地表的主磁场径向分量、其一阶导数（长期变化，SV）和二阶导数（持续加速度，SA）的全球分布图，截断的球谐阶数为20阶。

地球表面径向磁场分量具有一些显著特征，比如南大西洋的高纬度斑块和赤道区域较弱的径向磁场。径向场一阶导数显示，最大的径向场增加区域位于南美洲的东北角，同时有一条从南部非洲向西南延伸的径向磁场带。地球表面的径向场二阶导数显示了太平洋中强烈的局部偶极结构（在图5的中部区域可见），东太平洋（包括夏威夷附近）出现负加速度、澳大利亚东北部与西太平洋出现正加速度^［20-21］。

图6  核幔边界地磁场径向分量强度（2019年1月，截断阶数13阶）

Fig. 6  Radial magnetic field at core-mantle boundary (Jan. 2019, up to SH degree 13)

图7  核幔边界地磁场径向分量强度一阶导数（2019年1月，截断阶数17阶）

Fig. 7  First time derivative of radial magnetic field at core-mantle boundary (Jan. 2019, up to SH degree 17)

图6—8 显示了时间历元为2019 年1 月的CHAOS-7 模型在核幔边界的主磁场径向分量、其一阶导数（长期变化，SV）和二阶导数（持续加速度，SA）的全球分布图。针对核幔边界的研究区域，相关研究者指出在几年或更长的时间尺度内，地幔中电磁感应在与时间有关的内源场模型中起次要作用。同时，为了确保核幔边界的地磁场强度（包括一阶导数和二阶导数）图的图像平滑与稳定，实验选取的截断阶数分别为13、17、15^［16］。对于核幔边界的信号而言，在非洲与南美洲之间的赤道地区，存在最强的地磁场径向分量一阶导数信号，该特征是区域径向场向西运动的结果^［15-16］；对于核幔边界的二阶导数图，可以看出，位于东经$\mathrm{60}°~\mathrm{90}°$的太平洋地带存在大规模的加速现象，这与地表特征相对应，但位于西经$\mathrm{60}°~\mathrm{90}°$地区的大规模减速现象相对于地表特征来说并不明显。

图8  核幔边界地磁场径向分量强度二阶导数（2019年1月，截断阶数15阶）

Fig. 8  Second time derivative of radial magnetic field at the core-mantle boundary (Jan. 2019, up to SH degree 15)

3.3　重力场数据统计分解

以纬度$-\mathrm{60}°~\mathrm{60}°$之间的地区作为实验区域，采用PCA-ICA方法将实验区域内的GRACE重力残项序列分解为相应的空间模态和时变模态。图9展示了采用PCA-ICA分解后的不同数量GRACE模式的能量占比，图中横坐标表示GRACE的模式数，纵坐标表示模式能量之和占总能量的比例，例如横坐标为10时，其对应点的纵坐标表示前10个GRACE重力残项模式的能量之和占所有GRACE模式能量之和的分数。

图9  不同数量GRACE模式的总能量占比

Fig. 9  Percentage of total energy in GRACE modes

根据图9的统计结果，结合有关学者的研究成果，选取PCA-ICA分解的前6个模式作为主要模式进行研究。图10展示了全球实验区域6个主要模式的空间模态。对应实验区域6个主要模式的时变模态见图11所示。

图10  PCA-ICA方法分解出的GRACE空间模态

Fig. 10  Spatial mode of GRACE using PCA-ICA

图11  PCA-ICA方法分解出的前6个GRACE模式的时变模态

Fig. 11  Time-varying mode of GRACE modes using PCA-ICA

由10图可知，前2个模式存在明显的季节性特征，后续4个模式则反映了长周期变化，其中除了第3模式，另外3个模式除了长周期信号特征还反映了实验区域的季节性信号。前6个模式的能量占比分别为21.60%、13.23%、9.27%、5.45%、3.29%、2.72%。

3.4　地磁场数据统计分解

根据已有研究成果可知，重力残项与磁场径向分量的二阶导数存在一定相关性，本实验重点针对地磁场模型径向分量的二阶导数数据，采用相同的PCA-ICA分解方法对实验区域（纬度$-\mathrm{60}°~\mathrm{60}°$）内的磁场模型进行统计分解。

3.4.1　地球表面地磁场

对于地磁场径向分量的统计分解实验，与重力场数据相统一，选取PCA-ICA分解的前6个模式作为主要模式。不同数量地球表面地磁场模式的总能量占比见图12。图13展示全球实验区域6个主要模式的空间模态，前6个模式的信号能量占比分别为39.42%、17.73%、16.05%、8.53%、5.18%、4.09%。其对应的6个模式的时变模态结果如图14。

图12  不同数量地球表面地磁场模式的总能量占比

Fig. 12  Percentage of total energy in geomagnetic field modes on surface of the earth

图13  PCA-ICA方法分解出的地球表面地磁场（二阶导数）空间模态

Fig. 13  Spatial mode of second time derivative of radial magnetic field on surface of the earth using PCA-ICA

图14  PCA-ICA方法分解出的前6个地球表面地磁场（二阶导数）时变模态

Fig. 14  Time-varying mode of second time derivative of radial magnetic field on surface of the earth using PCA-ICA

3.4.2　核幔边界地磁场

对于核幔边界地磁场径向分量的统计分解实验，选取PCA-ICA分解的前6个模式作为主要模式。不同数量核幔边界地磁场模式的总能量占比见图15。图16展示全球实验区域6个主要模式的空间模态，其信号能量占比分别为53.00%、17.70%、15.38%、6.72%、2.68%、2.03%。其对应的6个模式的时变模态结果如图17。

图15  不同数量核幔边界地磁场模式的总能量占比

Fig. 15  Percentage of total energy in geomagnetic field modes at the core-mantle boundary

图16  PCA-ICA方法分解出的核幔边界地磁场（二阶导数）空间模态

Fig. 16  Spatial mode of second time derivative of radial magnetic field at the core-mantle boundary using PCA-ICA

图17  PCA-ICA方法分解出的前6个地球表面地磁场（二阶导数）模式时变模态

Fig. 17  Time-varying mode of second time derivative of radial magnetic field at the core-mantle boundary using PCA-ICA

3.5　重力信号与磁场信号的对比分析

根据图11中6个模式对应的时变模态信号特征分析可得，第1、4模式存在明显的季节性周期特征，应该是由于在进行陆地水文信号、海洋质量变化信号等非地核信号的扣除时并未扣除干净，导致重力残项中存在海洋洋流或是陆地水储量的季节性周期信号残留，对于第2、3、5、6模式，其对应的时变模态主要是长周期信号，信号的能量强度分别为13.23%、9.27%、3.29%、2.72%，其中可能包含有核幔边界的特征信号。将重力残项的时变模态与地磁场信号6个模式的时变模态进行分析对比（图14和图17），结合有关学者的研究成果，认为第3模式主要表征地核信号的变化特征^［15-16］。

核幔边界和地球表面地磁场第2模式的时变模态均与重力残项第3模式的时变模态存在高度相关性，三者时变模态的信号特征如图18所示。

图18  重力残项第3时变模态与地磁场（地球表面与核幔边界）第2时变模态信号

Fig. 18  Third time-varying mode of gravity field and second time-varying mode of geomagnetic field (surface and core mantle boundary of the earth)

根据图18可以看出，考虑其时间模态，对于重力残项的第3模式与地磁场信号的第2模式而言，均存在较高的相关性；考虑其空间模态，结合图5和图8，地磁场径向分量的二阶导数在太平洋地区也存在较为明显的变化信号，该变化特征在经过PCA-ICA分解后得到凸显（见图13），尤其是地球表面的地磁场信号，其第1、2、4、6模式在太平洋地区均表现出较强的变化信号，与重力残项的第3模式空间模态特征相对应。

自2014年以来，太平洋地区地磁场径向分量的二阶导数存在着线性迅速变化的现象，局部地区的加速度变化在地球表面呈现东西向排列偶极子形式的简单结构；对于重力信号的第3模式，在2014年以后，重力信号也呈现出线性迅速变化的现象，与磁场信号的时变特征存在着较高的相关性，但其第3空间模态在太平洋低纬度地区并未表现出明显的“偶极子”形式的变化特征，其地区的空间模态存在较多的短波变化。

同时对比分析CHAOS-7模型重构地球表面和核幔边界的地磁场结果可知，目前将径向分量二阶导数从地球表面深入到核幔边界仍存在一定难度，因为其功率谱在核幔边界处随着球谐阶数的增加而增加，意味着占主导地位的短波长成分并未受到良好的约束，此外，下地幔的电导率也存在不确定性，增加了向核幔边界反演的难度，尽管如此，在太平洋中部低纬度地区地球表面的地磁场二阶导数信号依然能在一定程度上反映核幔边界相应地区的一致性变化。

结合地磁场径向分量二阶导数信号（见图5、图8）及其第2模式的时变模态与空间模态（见图18）可知，本实验成功探测出自2014年以后太平洋地区地磁场二阶导数信号快速线性变化的特征，同时发现，重力残项第3模式的时变模态的信号变化特征与地磁场径向分量二阶导数的时变模态的信号变化特征存在较好的相关性（见图18）。

4 结语

基于2002年4月至2015年12月的GRACE卫星数据，使用公开发布的陆地水文模型（GLDAS模型）、海洋模型（ECCO模型）以及冰川均衡调整模型（ICE6G-D模型）对重力信号进行改正，以消除陆地水文、陆地冰川、海洋质量变化等非地核信号，从而分离出有效的地核信号。对于地磁场模型，采用CHAOS-7模型，反演出地球表面和核幔边界的地磁地核场时变模型，时间跨度覆盖2002年1月至2015年12月。采用PCA-ICA方法，对重力残项信号与地磁场信号进行统计分解，最后基于分解出的空间模态与时间模态对比分析两者的相关性及一致性。

研究发现，重力残项的第1、4模式为尚未扣除干净的陆地水文信号和海洋洋流信号，存在明显的周期性特征，其余模式存在长周期特征，其中重力信号第3模式的时变特征与磁场信号第2模式的时变特征在同期内存在相同的变化趋势，两者表现出较强的相关性，尤其在2014年以后，重力信号与磁场信号均呈现线性快速变化，与文献［4-6］研究结果符合。磁场信号的空间模态反映了太平洋低纬度地区较为明显的“偶极子”形式的变化特征。由于地磁地核场信号起源于地核，且实验所得重力残项第3模态与地磁场信号表现出较强相关性，所以认为本实验得到的重力信号中包含来自地核动力学过程的特征信息。综上所述，卫星重力信号在经过地表信号的剔除后，其包含的地核信号强度有所凸显，也能在一定程度上反映核幔边界处的动力学过程。

作者贡献声明

王正涛：提出研究主题、构建研究框架、提出研究思路。

张勇刚：完善实验方法、收集实验数据、进行数据实验、模型推演与论文撰写。