本文采用美国德克萨斯大学空间研究中心所提供的2003-01—2014-12的GRACE RL06时变地球重力场模型数据，分辨率为1个月，在其数据处理过程中，已扣除高频海洋大气质量变化及各种潮汐的影响，所以GRACE时变重力场模型在陆地区域所反映的主要为陆地水储量变化信息。利用GRACE时变重力场模型计算等效水高的公式为^[19]：

式中，$a$为地球平均半径，$a=$6 378.137 km；${\overline{P}}_{lm}(\mathrm{c}\mathrm{o}\mathrm{s}\theta )$是完全规则化的缔合勒让德函数，其中，m和l为球谐位系数对应的次数和阶数；${\rho }_{\mathrm{a}\mathrm{v}\mathrm{e}}$表示地球平均密度，${\rho }_{\mathrm{a}\mathrm{v}\mathrm{e}}$=5 517 kg/m³；${\rho }_w$表示水的平均密度，${\rho }_w$=1 000 kg/m³；$k_l$为顾及地球弹性形变的负荷勒夫数；$(\theta ,\lambda )$表示地球表面某点的纬度和经度；$W_l$为滤波函数。由于高阶次噪声的影响，对球谐位系数截断到60阶。由于GRACE数据2阶项噪声较大，利用卫星激光测距（satellite laser ranging，SLR）测得的$C_{\mathrm{20}}$项代替GRACE的$C_{\mathrm{20}}$项^[20]。采用半径为300 km的高斯滤波来降低球谐位系数的噪声，并结合Chen等^[21]提出的P3M6去相关滤波来消除南北条带误差的影响。采用单一尺度因子方法对GRACE数据由于截断和平滑滤波等引起的泄露误差进行改正，具体步骤为：将研究区域内信号设为1，区域外设为0，对该信号球谐展开，截断至与GRACE球谐位系数相同的阶次，再进行与GRACE数据相同的滤波处理，得到研究区域内的剩余信号，此时的尺度因子为1/剩余信号。

水文模型数据采用由美国宇航局和美国海洋与大气局国家环境预报中心共建的GLDAS数据，该模型是利用先进的地表建模和数据同化技术生成的最佳地表状态和通量场，被广泛应用于水文气象研究中^[22]。时间段选取与GRACE数据相同的2003-01‍—2014-12，滤波方法采用与GRACE数据相同的方法进行处理。数据空间分辨率为1°×1°，时间分辨率为1个月。

引入由全球降水气候中心（Global Precipitation Climatology Center，GPCC）收集全球众多气象站点降雨资料合成的全球降雨量数据集^[23]，采用与GRACE数据相同的滤波方法进行处理。数据空间分辨率为0.5°×0.5°，时间分辨率为1个月。该数据产品由美国海洋大气管理局的地球系统研究实验室提供。

利用WaterGAP全球水文模型（waterGAP global hydrology model，WGHM）数据来计算地下水变化，该模型数据对土壤水、雪水及地下水进行了严密估计，并多次被国内外学者运用到水储量变化的评估研究中^[24-26]，由于只收集到1960—2009年的数据，因此，本文只利用了2003-01—2008-12的数据进行地下水计算，并采用与GRACE数据相同的滤波方法进行处理。数据空间分辨率为0.5°×0.5°，时间分辨率为1‍个月。

设利用GRACE重力场模型计算$m$个月的水储量变化信息，其中，每个月含$n$个观测点，则$m$个月的所有观测值可用矩阵$\mathit{B}$表示：

式中，${\mathit{B}}_i$（i∈[1，n]）为$m\times \mathrm{1}$维列向量，表示第$i$个观测点对应不同月份的观测值。矩阵$\mathit{B}$的自协方差矩阵可表示为：

矩阵$\mathit{B}$的PCA分解可表示为：

式中，$\mathit{E}$由$\mathit{B}$的特征向量组成，代表PCA分解所得空间分布，在求解过程中，要对特征向量进行归一化处理，将特征值按照降序排列，并使得特征值和特征向量相对应；$P$由主成分组成，代表PCA分解所得时间序列，分解得到的时间序列无实际物理意义，仅表示数值大小。

ICA方法可以视为PCA方法的一种扩展方法，其目的是为了将PCA分解的时间序列或者空间特征尽可能分解为相互独立的成分^[27]。利用ICA算法对PCA所得时间序列或者空间特征进行旋转，得：

式中，$\mathit{U}$为源信号矩阵；$\mathit{S}$为混合矩阵；$r$表示所进行旋转的时间序列或者空间特征的数量；$\mathit{R}$表示正交旋转矩阵（$\mathit{R}{\mathit{R}}^{\mathrm{T}}=\mathit{I}$），该矩阵一般是利用特征矩阵联合近似对角化方法计算得到。本文是利用ICA算法对PCA分解的空间特征进行旋转，使研究区域的信号在空间特征上尽可能表现为独立成分，对ICA分解所得结果进行分析时，要将空间特征和其所对应时间序列相结合进行分析。

本文研究区域选定为非洲南部（10°E~40°E，24°S~3°N），该区域除奥卡万戈三角洲地区外，还有维多利亚湖、刚果河流域南部区域、坦噶尼喀湖、马拉维湖等地区，其概略位置关系见图 1。

图  1  研究区域概况图

Figure 1.  Survey Map of the Study Area

利用多项式拟合法拟合出研究区域2003-01—2014-12水储量变化的长期趋势如图 2所示。该区域内水储量变化信号主要在奥卡万戈三角洲地区和维多利亚湖地区，因此，本文直接分析了奥卡万戈三角洲地区和维多利亚湖地区的水储量变化，其时间序列如图 3所示。

图  2  研究区域水储量变化长期趋势

Figure 2.  Trend Map of Mass Variations of the Study Area

图  3  由GRACE数据计算研究区域特征性水储量变化信息

Figure 3.  Characteristic Change Information of Water Variations in the Study Area Calculated by GRACE Data

由图 3中奥卡万戈三角洲地区水储量的时间序列可以看出，该地区水储量变化有很强的年周期性，并且在2006—2012年呈现上升趋势，这与该地区实际情况相符合，实际资料显示，在2004‍—2012年，奥卡万戈三角洲地区年降水量与前25年相比，增加了约15%。利用GPCC降雨数据计算的1994-01—2003-12和2004-01—2011-12的月平均降雨量如图 4所示。

图  4  2004年前10年和后8年月平均降雨量

Figure 4.  Monthly Average Rainfall in the First 10 Years and Last 8 Years of 2004

经计算，奥卡万戈河流域1994-01—2003-12的月平均降雨量为72.25 mm，2004-01—2011-12的月平均降雨量为82.31 mm，即后8年较前10年月平均降雨量增加了13.9%。从图 3维多利亚湖地区水储量变化的时间序列可以看出，维多利亚湖地区的水储量在2007年前呈现出下降趋势，这是因为2003—2006年东非发生过干旱，且在2006年达到最小值，这与2006年东非严重干旱相吻合；在2007—2011年相对平稳，在2011—2014年呈现出上升趋势。

结合图 3中两个地区的时间序列变化特点，以不同时间节点将整个研究时段分开拟合出两个地区在不同时间段内的变化趋势如表 1所示。

将图 1所示的整个研究区域按照经纬度分割为1°×1°的格网等效为观测点，即该区域共有868个观测点，将平铺在格网上的等效水高值作为观测值来构造矩阵。在对整个研究区域解释率达到95%的情况下，提取出了11个主成分，特征值及累计方差贡献率如图 5所示。

图  5  PCA分解的特征值及方差贡献率

Figure 5.  Eigen Value and Variance Contribution of PCA

图 5中，前3个主成分的累计方差贡献率达到了90.2%，其中，PC1、PC2、PC3的方差贡献率分别达到了70.5%、11.6%和8.1%。由PCA方法的实际意义可知，前3个成分已经包含了该变量场中的主要信息，即前3个主成分可以解释该区域内的主要水储量变化情况。前3个主成分的空间特征分布及相对应的时间序列见图 6和图 7。

图  6  PCA分解前3个成分的空间特征分布

Figure 6.  Spatial Feature Distribution of the First Three Components in PCA Decomposition

图  7  PCA分解前3个成分的时间序列

Figure 7.  Time Series of the First Three Components in PCA Decomposition

从图 6（a）中PC1的空间特征分布可以看出，PC1中强信号主要集中在奥卡万戈河流域东北部与赞比亚境内，且为正值。从图 7（a）中PC1的时间序列可以看出，该地区水储量变化呈现出很强的季节周期性，结合空间特征分布和时间序列可以看出，从2006年以后呈现出缓慢增加的趋势。

从图 6（b）中PC2的空间特征分布看出，在奥卡万戈河流入奥卡万戈沼泽的汇集地区域出现了较强的信号；从图 7（b）可以看出，PC2的时间序列是周年信号和半周年信号及长期趋势的混合，且呈现出增加的长期趋势，PC2的空间特征分布在该地区的信号显示为正值。因此，可以分析在奥卡万戈三角洲地区水储量呈现增加趋势，这与实际资料显示该地区水储量增加的信息相符合。

图 6（c）中PC3的空间特征分布表示出的最强信号在维多利亚湖地区。图 7（c）中PC3的时间序列也表现出的是周期信号和长期趋势的混合。可以看出，2006年附近水储量明显小于其他年份，这与2006年东非地区干旱相符合。从PC3时间序列还可看出，在2011—2014年水储量表现为缓慢增加，且PC3的空间特征分布在维多利亚湖区域显示为正值。两者相结合分析可得，2011年后维多利亚湖地区的水储量呈缓慢增加趋势。

在PC3的空间特征分布中，虽然最强信号出现在维多利亚湖区域，但是在维多利亚湖区域西侧的刚果河流域和奥卡万戈三角洲区域都有很明显的信号，这导致图 7中PC3的时间序列不能代表维多利亚湖的主要变化，还受到其他较强信号的影响。为了改善PCA的分解结果，将区域性特征信号更好地分解出来，把上述PCA分解出的空间特征分布进行旋转，进行ICA分解，得到分解后的空间分布特征和时间序列分别见图 8和图 ‍9。

图  8  ICA分解前3个成分的空间特征分布

Figure 8.  Spatial Feature Distribution of the First Three Components in ICA Decomposition

图  9  ICA分解前3个成分的时间序列

Figure 9.  Time Series of the First Three Components in ICA Decomposition

从图 8中IC1和IC2的空间特征分布来看，ICA分解的前两个成分和PCA分解的前两个成分的信号大致相同；IC3与PC3的空间特征分布图相比，很明显强信号更加集中在维多利亚湖区域，而PC3出现的刚果河流域和奥卡万戈三角洲的较强信号在IC3的空间特征分布很小，不会对维多利亚湖地区的较强信号造成影响。

结合图 9中的IC1和IC2的时间序列来看，IC1和PC1的时间序列基本相同；IC2和PC2的时间序列变化趋势相同，只在个别年份有所不同。IC3时间序列在2006年的水储量明显减小到邻近年份的最低值，与图 3中维多利亚湖地区水储量时间序列中情况吻合，比图 7中PC3的时间序列的特征信号的表现更加明显。与PCA相比，ICA不仅在空间上将信号尽可能分离成独立信号，而且在时间尺度上提取出特定信号。与多项式拟合方法相比，ICA不需要建立拟合模型和任何先验信息，且不需要单独针对特定区域或者特定点，可以同时提取较大区域变量场中质量变化的时空特征信号，减少了模型误差及选取区域不合适等问题带来的困扰。

采用与处理GRACE数据相同尺度的ICA方法处理GLDAS水文模型数据，在累计方差贡献率达到90%的情况下，提取出了前4个成分。ICA方法分解GLDAS数据前4个成分的空间特征分布和对应的时间序列如图 10和图 11所示。

图  10  ICA分解前4个成分的空间特征分布

Figure 10.  Spatial Feature Distribution of the First Four Components in ICA Decomposition

图  11  ICA分解前4个成分的时间序列

Figure 11.  Time Series of the First Four Components in ICA Decomposition

从图 10中IC1的空间特征分布可以看出，GLDAS的分解结果与GRACE分解结果在第‍1‍成分中都显示了奥卡万戈河流域东北部与赞比亚境内的信号。结合图 11中IC1的时间序列可以看出，GLDAS水文数据在该地区也表现出很强的周期性变化。

由图 10中IC2的空间特征分布可以看出，信号主要出现在奥卡万戈三角洲北部的刚果河流域东南部地区。图 11中IC2的时间序列是长期趋势和周期性信号的混合。总体来看，其长期趋势呈现出缓慢下降的特征，在IC2的空间特征分布中，该地区的信号呈现负值，其时间序列和空间特征相乘表示水储量变化重构值，因此，该地区水储量呈现出缓慢增加的趋势。

由图 10中IC3的空间特征分布可以看出，信号主要出现在奥卡万戈三角洲地区，该信号与GRACE-IC2信号位置相同，说明GRACE计算结果与GLDAS水文数据的计算结果在该地区符合得很好。图 11中IC3的时间序列呈现出先上升后下降趋势，且在上升期内与GRACE-IC2的时间序列趋势相同。GLDAS水文数据的IC4的信号也出现在奥卡万戈三角洲的西南部地区，时间序列显示了一定的周期性信号。而对于GRACE-IC3中出现的维多利亚湖地区的信号，在GLDAS结果中却没有显示，这是由于GLDAS数据所表示的水储量是陆地表层水，不包含地下水和内陆湖泊水，而GRACE所表示的水储量是表层水和地下水的一个综合信息的反映，因此，GLDAS数据在维多利亚湖区域信号不明显。

为了更直观地对比GRACE数据处理结果与GLDAS数据处理结果是否具有相似性，以便验证GRACE数据所得结果的正确性，本文对两者空间特征分布信号位置相同的成分所对应的时间序列进行了相关性分析，结果如图 12所示。

图  12  GRACE和GLDAS相似成分结果对比

Figure 12.  Comparison of Similar Compositions Between GRACE and GLDAS

GRACE-IC1和GLDAS-IC1的空间特征分布中的信号出现在相同位置，因此，图 12（a）将GRACE-IC1和GLDAS-IC1进行了相关性分析。可以看出，两者的结果在此区域都显示出很强的周期性，说明第1成分所提取出的是该地区的1‍个周期性信号，相关系数为0.85，表示两者结果为极强相关。GRACE-IC2和GLDAS-IC3空间特征分布的信号出现在相同位置，两者的相似信号没有出现在同一成分，是由于GRACE数据和GLDAS数据的模型及对变量场中所含成分占比不同等因素导致的。从图 12（b）中可以看出，两者在上升阶段，整体趋势基本相同，为了更好地对比分析两者在长期趋势上的相似性，对两者的时间序列进行了13点滑动平均处理，可扣除季节性波动的影响，处理结果如图 13所示。

图  13  GRACE-IC2和GLDAS-IC3时间序列13点滑动平均结果对比

Figure 13.  Comparison of 13-Point Moving Average Results of GRACE-IC2 and GLDAS-IC3 Time Series

经过13点滑动平均处理后，可以明显看出，在2011-10（图 13中绿线）前，GRACE和GLDAS的数据结果均呈现出上升趋势，相关系数为0.81，为极强相关，说明该地区水储量在2011-10前呈增加趋势。在2004-08—2005-07，GLDAS的数据结果呈上升趋势，而GRACE的数据结果表现为下降趋势，原因是该时间段水储量变化主要受到地下水影响。利用WGHM水文模型计算出该地区2003-01—2008-12的地下水变化，结果如图‍14（a）所示。对该结果进行13点滑动平均处理，再将GRACE的结果和GLDAS的结果进行对比，如图 14（b）所示。

图  14  WGHM水文模型计算地下水变化结果

Figure 14.  WGHM Hydrological Model Calculates Groundwater Variation Results

从图 14可以看出，在2004-08—2005-07，地下水变化与GRACE的结果均表现为下降趋势，证明了水储量变化确实受地下水变化主导。在2011-10—2013年底，该地区水储量呈下降趋势，2011-10前一年（2010-11—2011-10）和后一年（2011-11—2012-10）平均降雨量对比见图 15。

图  15  2011-10前后1年平均降雨量

Figure 15.  Monthly Average Annual Rainfall of One Year Before and After October 2011

从图 15可以看出，2010-11—2011-10，该地区月平均降雨量为85.78 mm，2011-11—2012-10，该地区月平均降雨量为78.34 mm，同比减少了9.5%，因此，该时间段水储量下降是由于2011-10后一年（2011-11—2012-10）降雨量比前一年（2010-11—2011-10）减少，导致地表水迅速减少。

对于GLDAS数据的下降趋势要比GRACE数据的下降趋势更快，这是由于前面所述的GRACE和GLDAS数据本身的差异导致的，GLDAS数据只含有陆地表层水，因此，对于地表水的减少反映更加灵敏，而GRACE数据反映的是地表水及地下水的一个综合信息反映。

本文利用GRACE时变重力场模型数据分析了非洲南部奥卡万戈三角洲及周边地区的水储量变化，利用PCA方法和ICA方法对结果进行分解，并对结果进行了分析。将GLDAS水文模型数据与GRACE数据进行相同的滤波处理，利用与GRACE数据相同尺度的ICA方法进行分解，结合GPCC降雨数据和WGHM水文模型数据对该地区水储量变化特点进行分析，所得结论如下：

1）非洲奥卡万戈河流域地区东北部水储量变化显示出很强的周期性信号，且GRACE数据和GLDAS数据计算结果符合程度较高，相关系数达到0.85，表现为极强相关，验证了GRACE数据分析水储量结果的正确性。

2）在奥卡万戈三角洲区域水储量长期变化趋势方面，该地区GRACE数据计算结果和GLDAS数据计算结果在2003-01—2011-10的整个水储量上升期内，去除季节性波动影响后，两者结果符合很好，相关系数达到0.81，说明了该地区水储量在2003-01—2011-10之间呈增加趋势，这与该地区实际情况相符，2004—2012年该地区降雨量增加导致水储量增加。而在2010-10后，降雨量明显减少导致地表水减少，因此，GLDAS计算结果迅速下降，而GRACE相对下降较为缓慢，这是由于两种数据含量不同导致的结果，GLDAS数据只含有陆地表层水，而GRACE所表示的水储量信息是地表水及地下水的综合反映。

3）在不需要先验信息的情况下，ICA方法可以将变量场的信息尽可能地分解成相互独立的成分，并且能够得到每个信号的空间特征分布和相对应的时间序列，这有助于分离一个地区不同特征的水文信号，并进行单独的定性分析。且在空间和时间上，ICA方法对信号的分解能力都要强于PCA方法。因此，ICA方法更加适用于对区域内部复杂信号的分离提取分析。