华北地区是中国七大地理区划之一，包括山西省、河北省、北京市和天津市，其中京津冀地区经济发达、人口密集，在中国政治、经济以及文化等方面发挥着重要作用^[1]。然而，由于各种自然和人为因素，如过度开采地下水灌溉农作物，华北地区面临着严重的地表垂直形变，社会的可持续发展也因此受到了较大影响^[1-3]。因此，必须对该地区的地表位移进行及时且有效的监测，充分了解形变的机理，从而为灾害的准确预测提供数据保障和理论支撑。

全球导航定位系统（global navigation satellite system，GNSS）是监测地表垂直形变最有效的空间大地测量技术之一，为不同时空尺度的地球物理现象研究提供了宝贵的基础数据^[4]。利用GNSS数据解算的坐标时间序列定义为一组按照时间顺序排列的基准站坐标^[4]，其中包含了丰富的构造及非构造信号，后者主要表现为地表质量重分布、基岩热膨胀效应、共模误差等因素导致的季节性变化^[2，5]，准确地估计GNSS季节性位移并识别相应的源信号有助于深入理解固体地球形变与外层质量变化之间的相互作用。通常情况下，首先假定季节性信号的振幅为常数，然后通过最小二乘方法进行估计，然而恒定的振幅与真实季节性激发源呈现的非均匀变化相悖。因此，以经验正交函数（empirical orthogonal function，EOF）^[6]、主成分分析（principal component analysis，PCA）^[7-8]和奇异谱分析（singular spectrum analysis，SSA）^[9]为主的统计分析方法受到了广泛关注，它们可在无任何先验信息的条件下将GNSS时间序列分解为若干个正交分量，并提取其中振幅随时间变化的季节性分量。然而，上述方法仅考虑了时间序列的方差和协方差二阶统计量，未顾及到更高阶统计信息，由此获得的季节性分量很大程度上是多个物理源信号的混合^[5，10]。

若将GNSS时间序列视为来自不同物理源的独立信号的线性混合，由于混合模式未知，GNSS源信号提取实际是一个盲源分离问题。独立分量分析（independent component analysis，ICA）是处理该问题的有效方法^[10-11]，它在PCA的基础上引入高阶统计量，使分解的成分尽可能相互独立，从而有效提高了潜在源信号的可识别性。近年来，ICA已成功应用于分析各类地球物理信号，如陆地水储量^[12-13]、GNSS共模误差^[10，14-15]、地震信号^[16]等。鉴于ICA方法在信号分析中的独特能力，本文以华北地区为研究对象，利用ICA方法对20个GNSS基准站的垂直位移时间序列进行分析，探讨各独立分量的时空变化特性，评估与环境负载信号的相关性，由此确定驱动华北地区季节性垂直位移的主要源信号。

1 数据获取与处理

1.1   GNSS时间序列

本文通过中国大陆构造环境监测网络（crustal movement observation network of China，CMONOC）获取了20个GNSS基准站的垂直位移时间序列（http://www.cgps.ac.cn），时间跨度为2011-01-01—2019-11-30，站点位置分布如图 1所示。

图  1  研究区域GNSS站点分布

Figure  1.  Distribution of GNSS Stations in Study Area

下载: 全尺寸图片 幻灯片

数据预处理包括：根据台站日志文件修正设备故障和更换等引起的非地球物理偏移，同时通过地震统计信息（http://www.csi.ac.cn/）改正由同震及震后引发的偏移，改正量为偏移前、后各7天平均值的差值；利用ICE6G-D模型^[1]，通过球谐系数法^[17]估计冰川均衡调整对地壳垂直运动产生的效应，并将其从各站点时间序列中扣除，该效应在华北地区表现为长期变化趋势，范围为0.27~0.31 mm/a；基于四分位间距统计量检测异常值并标记为缺失数据，经计算得到各时间序列缺失比例的平均值为6.20%。

1.2   环境负载数据

环境负载主要指大气负载（atmospheric loading，ATML）、水文负载（hydrologic loading，HYDL）以及非潮汐海洋负载（non-tide ocean loading，NTOL）^[4，18]。根据以往研究^[2，5，19]，非潮汐海洋负载在华北地区诱发的垂直形变非常小，且时间序列较为复杂，因此本文暂不考虑。大气负载垂直位移取自德国地学研究中心利用非潮汐大气表面压力数据研制的格网数据产品^[20]（http://rz-vm115.gfz-potsdam.de:8080/repository），时间分辨率为3 h，空间分辨率为0.5° × 0.5°，对该格网数据产品进行克里金空间插值得到各站点位移时间序列，扣除其中的线性趋势后通过滑动平均将序列的时间分辨率转换为12 h。

重力场恢复与气候实验（gravity recovery and climate experiment，GRACE）卫星及后续计划GFO（GRACE follow-on）自发射以来，利用其数据解算的时变重力场模型已实现了全球或区域地表物质迁移的高精度、高时效监测^[18]。本文采用德克萨斯大学空间研究中心提供的2011-01—2019-11期间RL06月重力场数据表征水文负载，数据格式为96阶次球谐系数，且已扣除了海潮、固体潮、极潮、非潮汐大气和非潮汐海洋等影响，因此可反映陆地水的质量变化。相关后处理包括：因为GRACE/GFO无法监测地心运动，且对二阶项系数不敏感，所以在球谐系数中回加了一阶项系数，并替换了其中的C₂₀及C₃₀系数^[21]；结合P4M6去相关滤波和高斯平滑技术（半径取300 km）移除南北条带并抑制高频噪声^[22]；利用ICE6G-D模型移除冰川均衡调整效应，以保证与GNSS的可比性；通过球谐系数法计算水文负载在各站点引起的垂直形变并扣除形变序列中的线性趋势^[17-18]。

为进一步验证GRACE/GFO模型估计水文负载位移的可靠性，分别利用地表流量模型（land surface discharge model，LSDM）^[20]和全球陆面数据同化系统的流域陆面模型（catchment land surface model，CLSM）^[23]计算水文负载造成的垂直位移。首先将两种水文模型的格网数据分别转化为球谐系数，截断至96阶次，保持与GRACE/GFO一致，然后利用P4M6和高斯平滑技术进行滤波，最后通过球谐系数法分别得到基于LSDM和CLSM模型的水文负载位移。GRACE/GFO、LSDM及CLSM三种模型获取的各站点平均垂直位移的时间序列如图 2所示，其中阴影部分表示GRACE和GFO之间的数据缺失时段。由图 2可以看出，GRACE/GFO与两种水文模型获取的平均位移时间序列较为一致，GRACE/GFO与CLSM和LSDM的相关系数分别为0.66、0.60，具有较好的一致性，表明了GRACE/GFO模型能够准确地估计水文负载位移。

图  2  三种模型获取各站点平均垂直位移的时间序列

Figure  2.  Time Series of Average Vertical Displacement Derived from Three Models

下载: 全尺寸图片 幻灯片

2 ICA方法

ICA方法是一种盲源分离技术，可通过最大化高阶统计量将混合的GNSS时间序列分解为若干个统计独立的源信号^[24]。设${\mathit{X}}_{m\times n}$为扣除趋势后的GNSS坐标时间序列矩阵，m为观测历元数，n为台站数。首先，通过PCA进行矩阵分解，计算式为：

式中，P和E的列向量分别称为主成分分量（principal component，PC）和特征向量，两者均按特征值由大到小对应排列。

若前$r\left(r<n\right)$个方差占比较大的PC代表GNSS中的信号^[8]，则式（1）可变换为：

P矩阵中各列向量虽然互不相关，但并不独立。因此，利用酉矩阵W对P矩阵进一步旋转，式（2）扩展为：

式中，$\mathit{S}=\mathit{P}\mathit{W}$，其每一列向量表示一个独立的时间分量；$\mathit{A}=\mathit{E}\mathit{W}$，其列向量通常称为空间响应^[10]。每个时间分量乘以相应的空间响应即表示一个独立分量（independent component，IC）。本文采用特征矩阵联合近似对角化算法估计酉矩阵W，该算法具有较好的鲁棒性，其详细过程参见文献[12]。

考虑到数据缺失的干扰，特别是缺失比例较大的站点（如SXKL站缺失率达26.58%），传统的迭代或插值方法可能会引入虚假信息，导致提取的信号产生偏差。因此，本文采用文献[8，24]提出的改进方法求解式（1）中的主成分，该方法无需预先插值，且效果优于传统插值方法。

3 GNSS季节性垂直位移时空分解

首先，通过最小二乘方法拟合GNSS垂直位移时间序列，获取各序列的长期趋势。拟合函数的计算式为^[5]：

式中，y₀是偏移量；b是长期趋势；t为历元数；$c_f$、$d_f$是谐波分量系数；${\displaystyle \sum _{f=\mathrm{1}}^{\mathrm{2}}}\left[c_f\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{2}\mathrm{\pi }ft)+d_f\mathrm{c}\mathrm{o}\mathrm{s}(\mathrm{2}\mathrm{\pi }ft)\right]$表示季节项，包括周年（f=1）和半年（f=2）变化，其振幅为$\sqrt[]{c_f^{\mathrm{2}}+d_f^{\mathrm{2}}}$；ε是残差。

剔除长期趋势后，仅利用式（4）的季节项函数对去趋势后的GNSS序列拟合周年和半年信号，结果如图 3所示。由图 3可以看出，所有时间序列的周年振幅远大于半年振幅，范围分别为1.67~4.46 mm和0.18~1.27 mm。

图  3  各台站季节性垂直时间序列及拟合的周年和半年信号

Figure  3.  Seasonal Vertical Time Series and Fitted Annual and Semi-annual Signals

下载: 全尺寸图片 幻灯片

然后，利用ICA方法对GNSS季节性时间序列进行时空分解，根据特征值大小将主成分降序排列后，计算其累计方差占比。GNSS、大气和水文负载位移分解时主成分的累计方差占比如图 4所示。由图 4可知，前6个GNSS-PC累计方差占比达86%，在很大程度上已能表示所需的季节性特征，因此对前6个GNSS-PC进一步分解得到6个GNSS-IC。

图  4  GNSS、大气和水文负载位移的主成分累计方差占比

Figure  4.  Cumulative Variance Percentages of PCs in Decomposition of GNSS, ATML and HYDL Displacement

下载: 全尺寸图片 幻灯片

鉴于IC排列顺序具有随机性且方差均为1，采用文献[5，19]定义的贡献率对6个GNSS-IC重新排序，获得的时间分量和空间响应如图 5所示，其中红色箭头的长度和方向分别表示空间响应的量级和符号，箭头向上为正，向下为负，绿色虚线框表示具有空间集聚特征的局部区域。为使结果更加直观，将每个GNSS-IC的空间响应值除以绝对值最大的响应值进行归一化，而时间分量乘以该归一化因子实现标准化，因此所有的空间响应值均介于-100%~100%之间。此外，为了准确判定GNSS-IC的周期特性，计算其时间分量的功率谱，如图 6所示。

图  5  前6个GNSS-IC的时间分量和空间响应

Figure  5.  Temporal Components and Spatial Response of the First Six GNSS-ICs

下载: 全尺寸图片 幻灯片

图  6  前6个GNSS-IC时间分量的功率谱

Figure  6.  Power Spectra of Temporal Components of the First Six GNSS-ICs

下载: 全尺寸图片 幻灯片

由图 5可以看出，GNSS-IC1是显著的周年信号，其空间响应方向完全一致，且响应强度在整体上相对均匀，由北至南有缓慢递增的趋势。

由图 6可以看出，GNSS-IC2和GNSS-IC3均呈现出明显的周年变化，且还包含了一些年内信号，尤其是GNSS-IC2，存在明显的半年周期。尽管GNSS-IC2和GNSS-IC3的空间响应在整个区域内表现出符号一致性，但其量级并非完全系统或完全随机，反而具有局部集聚特征（空间集聚特征判定依据为局部区域的空间响应具有相同方向且量级差异小于10%^[10]），见图 5中绿色虚线框内站点。该分布特征与华北地区的水资源分布状况较类似^[3，25-26]，由于降水的时空分布极不均匀，各子区域的人类活动也不尽相同，如山西省大量的煤矿开采、河北省南部广泛的农业种植、京津地区快速的工业发展等，导致华北各子区域水资源损耗差异较大^[26]，由此推断GNSS-IC2和GNSS-IC3可能是针对不同类型水文质量变化的形变响应。GNSS-IC4~GNSS-IC6的季节性特征较弱，其空间响应分布也更为复杂，应是GNSS站点坐标中的共模误差^[7，10]，需要进一步分析。总体而言，利用ICA分解华北地区季节性垂直位移能够有效提取隐含的季节性周期变化，并反映不同季节性信号的空间分布迥异性。

4 GNSS季节性垂直位移源信号分析

4.1   环境负载位移时空滤波

采集的地球物理数据必定会受到噪声污染，因此大气和水文负载位移也需要进行相应的ICA滤波^[10，14-15]，求解的ICs分别如图 7、图 8所示，其中红色箭头的含义与图 5一致，图 8中时间分量的阴影部分表示GRACE和GFO之间的数据缺失时段。

图  7  前3个ATML-IC的时间分量和空间响应

Figure  7.  Temporal Components and Spatial Response of the First Three ATML-ICs

下载: 全尺寸图片 幻灯片

图  8  前3个HYDL-IC的时间分量和空间响应

Figure  8.  Temporal Components and Spatial Response of the First Three HYDL-ICs

下载: 全尺寸图片 幻灯片

对于大气负载位移，前2个ATML-IC累计方差占比达99.71%（图 4）。从图 7可以看出，ATML-IC1是纯粹的周年信号，且空间响应在整个区域中呈现出高度一致性；ATML-IC2与ATML-IC1类似，但时间分量中还包含一些高频信号；ATML-IC3的空间响应呈现阶梯状分布，由北向南逐渐减小，然而其时间分量明显由高频噪声主导。因此，本文仅利用前2个ATML-IC重构大气负载位移时间序列。

对于水文负载位移，前3个HYDL-IC累计方差占比约97.36%（图 4）。从图 8可以看出，HYDL-IC1的时间分量较平稳且振幅最强，具有显著的周年变化，其空间响应也呈现出较好的一致性；反观HYDL-IC2和HYDL-IC3，两者的时间分量在整个时段内波动较大（如HYDL-IC2中2017年，HYDL-IC3中2011—2013年），年周期性较差但年内变化较为丰富。HYDL-IC2的强空间响应主要聚集于太行山沿线及北京市周边，表明该区域存在致使水文变化明显的自然或人类活动；HYDL-IC3的空间响应则表现出南北相反的特点，反映了某水文分量在华北地区南北部分别处于补给和损耗状态。综上，本文选择前3个HYDL-IC重构水文负载位移时间序列。

4.2   GNSS-IC与环境负载垂直位移的比较

GRACE/GFO数据仅能估计总水文质量变化（即地表水、雪水、土壤水、地下水等水文分量的质量变化之和）引起的形变，而GNSS-IC2和GNSS-IC3可能是对不同水文分量变化的形变响应，因此将GNSS-IC2和GNSS-IC3合并为GNSS-IC2/3，再将其与水文负载位移时间序列进行比较，同时将GNSS-IC1与大气负载位移时间序列进行比较，计算各站点Pearson相关系数，结果如图 9所示。

图  9  GNSS-IC与环境负载位移的相关系数

Figure  9.  Correlation Coefficients Between GNSS-ICs and Environmental Loading Displacements

下载: 全尺寸图片 幻灯片

由图 9（a）可知，GNSS-IC1与大气负载位移在整体上较为吻合，区域内各站点相关系数比较接近，范围为0.48~0.63，平均值为0.60，表明大气负载是华北地区季节性垂直位移最重要的驱动源。由图 9（b）可知，GNSS-IC2/3与水文负载位移存在较好的关联性，相关系数平均值为0.50，说明水文负载是华北地区季节性垂直位移另一个重要的驱动源。然而，相关系数分布反映出明显的南北差异，较大值主要位于南部，这与华北地区不同的地下水开采方式密切相关^[1，27]。南部地下水主要开采自潜水层，虽然损耗量较大，但易受到降水的补给，因此干扰地表位移季节性变化的程度有限。反之，中北部地下水主要存储于深水层，长期的过度开采致使地下水损耗极为严重且难以补给，形成了多个地下水沉降漏斗，从而产生了更为复杂的地表形变。

站点HECD、HETS、HECX的两组相关系数均显著小于其他站点，可能是受当地效应影响较大所致。此外，本文尝试将GNSS-IC（IC2、IC3）与HYDL-IC（IC1、IC2、IC3）以及利用地下水、土壤水计算的垂直位移进行交叉比较，但并未发现较明显的相关性。因此，ICA获取的分量能够反映何种程度的真实物理源信号还待商榷。

5 结语

本文基于ICA方法分析了2011—2019年华北地区20个GNSS基准站的季节性垂直位移时间序列的时空变化特性，并确定了其主要的源信号，得出结论如下：

1）利用ICA分解季节性垂直位移能够有效探测其中隐含的周期信号，并且反映各信号迥异的空间分布特性。其中GNSS-IC1为显著的周年信号，空间响应呈现相对均匀的分布特征；GNSS-IC2和GNSS-IC3均包含了周年和年内周期信号，但空间响应反映出局部集聚特征，与该地区不均匀的水资源分布状况有关。

2）利用ICA对大气负载位移进行时空分解，有效分离了其中的周年信号及高频噪声，表明了ICA可作为一种可靠的信号滤波方法。此外，在很大程度上分离了水文负载位移的周年和年内变化，且空间响应分布反映了不同水文分量变化造成的地表垂直位移。

3）GNSS-IC1与大气负载位移时间序列吻合较好，各站点的平均相关系数为0.60，说明大气负载是华北季节性垂直位移最重要的源信号；GNSS-IC2/3与水文负载位移时间序列在各站点的平均相关系数为0.50，说明水文负载是重要程度仅次于大气负载的源信号。

上述结果体现了ICA方法在GNSS时间序列信号提取中的良好表现。然而，本文尚未分离出其他因素诱发的季节性位移，如非潮汐海洋负载、基岩及观测墩热膨胀效应等，后续需进一步分析，使结果更加全面、客观。