2018, Vol. 38
目前，在GRACE水文学研究中，大多是直接利用谐波拟合等平稳时间序列分析GRACE得到的水储量时间序列变化特点，得到的主要是水储量变化的周年或半周年特性。由于陆地水储量变化具有很强的周年特性，所以利用多项式拟合可以得到其周年项。但是多项式拟合模型在拟合数据时，需知道拟合数据所具有的周期特性，导致无法拟合得出数据里周期特性并不清楚的分量。针对这种情况，可利用时频域分析方法将多种周期信号分离出来。常用的时频域分析方法包括小波分析、快速傅里叶变换(fast Fourier transformation, FFT)、经验模态分解(empirical mode decomposition，EMD)[1]等。其中，小波分析需要利用基函数，分解层数需要有先验条件；FFT需要利用得到的频谱图提出具体拟合方法；而EMD则具有自适应性、能对非平稳时间序列处理的优点，使其在时间序列处理方面比小波分析和FFT更加简单有效。Loomis等[2]利用EEMD/CEEMD方法对GRACE得到的等效水高序列进行分析，得到除周年和半周年外的水循环周期，并选取不同的周期分量信号进行叠加，得到水循环的年际变化和长期变化。但EMD存在模态混叠效应问题，即同一个本征模态(intrinsic mode function，IMF)中出现不同频率或尺度的信号。Wu等[3]研究白噪声信号的统计特性，提出总体平均经验模态分解(ensemble empirical mode decomposition，EEMD)，通过对原始信号多次加入不同的白噪声进行EMD分解，将结果取平均，得到最终IMF。Yeh等[4]提出一种补充的总体平均经验模态分解(complementary ensemble empirical mode decomposition，CEEMD)，其在EEMD的基础上添加两个相反的白噪声信号。郑进德等[5]在研究上述几个EMD模型后，提出改进的EEMD算法(modified EEMD，MEEMD)，利用排列熵来研究分解信号的随机性，避免了EEMD/CEEMD方法中不必要的集成平均，使分解结果中出现虚假信号的几率更小。然而在计算排列熵时，需要确定嵌入维数m和时间延迟λ两个参数，如果m过小，重构的向量中包含太少的状态，算法失去意义；如果m取值过大，相空间的重构将会均匀化时间序列，此时不仅计算比较耗时，而且无法反映序列的细微变化。而对于λ的选取，更没有一个明确的说法。MEEMD中采用经验选取得到m与λ，使得其不太灵活。
本文结合频谱分析和KS检验，对MEEMD进行改进，并将其应用于GRACE水循环周期的提取。利用2005~2010年的GRACE重力卫星数据研究长江、黄河、淮河和珠江流域的陆地水储量变化，并与水文模型GLDAS-NOAH[6]结果对比；利用改进的EEMD算法，提取使用GRACE得到的四大流域陆地水储量变化的周年项、年际变化和趋势项，并结合实际的陆地水储量变化进行分析。
1 南北状条带误差消除与流域水储量计算方法比较由于GRACE卫星轨道设计和观测精度的制约，原始的高阶球谐系数中存在较大的系统误差，导致直接利用球谐系数反演得的全球陆地水分布结果中存在南北状条带误差[7]。本文采用CSR提供的2005-01~2010-12共72个月的GRACE Level-2(Release-05)产品(已将非潮汐大气和高频海洋信号以及各种潮汐的影响扣除)，该数据产品最终以球谐系数的形式给出[8]。为减小南北条带误差的影响，采用如下处理程序[7]：
1) 球谐系数去均值。用卫星激光测距(SLR)获得的C20项替换GRACE数据的C20项[9]，并向球谐系数中加入Swenson等[10]计算的地心改正项，得到时变信号ΔClm、ΔSlm。
2) 使用去相关滤波(P3M8)去除同一次的奇偶项球谐系数中包含的系统误差，减小条带误差的影响，结合350 km高斯滤波，削弱球谐系数的高阶项噪声。
根据上述步骤，按照Swenson等[7]的方法，利用式(1)得到1°×1°分辨率的流域等效水高分布：
| $ \begin{array}{l} \Delta {\sigma _{{\rm{region}}}} = \frac{{a{\rho _{{\rm{ave}}}}}}{{{\rho _w}{\mathit{\Omega} _{{\rm{region}}}}}} \cdot \\ \sum\limits_{l = 0}^{60} {\sum\limits_{m = 0}^l {\frac{{2l + 1}}{{1 + {k_l}}}{W_l}\left( {\vartheta _{lm}^c\Delta {C_{lm}} + \vartheta _{lm}^s\Delta {S_{lm}}} \right)} } \end{array} $ | (1) |
式中，a为地球半径(6 371 km)，ρave为地球平均密度(5 517 kg/m3), ρw为水的密度(1 000 kg/m3)，Ωregion为流域对应的角面积，由dΩregion=cosθdθd
流域平滑核函数可由式(2)计算得到[7]：
| $ \vartheta \left( {\theta, \mathit{\emptyset} } \right) = \left\{ \begin{array}{l} 0, 在流域外\\ 1, 在流域内 \end{array} \right. $ | (2) |
截取至60阶，式(2)变为[7]：
| $ \vartheta \left( {\theta, \mathit{\emptyset} } \right) = \frac{1}{{4\pi }}\sum\limits_{l = 0}^{60} {\sum\limits_{m = 0}^l {{P_{lm}}\left( {\cos \theta } \right)\left\{ {\vartheta _{lm}^c\cos m\mathit{\emptyset} + \vartheta _{lm}^s\sin m\mathit{\emptyset} } \right\}} } $ | (3) |
| $ \left\{ \begin{array}{l} \vartheta _{lm}^c\\ \vartheta _{lm}^s \end{array} \right\} = \int {\vartheta \left( {\theta, \mathit{\emptyset} } \right){P_{lm}}\left( {\cos \theta } \right)\left\{ \begin{array}{l} \cos m\mathit{\emptyset} \\ \sin m\mathit{\emptyset} \end{array} \right\}{\rm{d}}\mathit{\Omega} } $ | (4) |
计算流域水储量还可利用流域纬度加权方法。为比较流域纬度加权方法与流域平滑核函数得到的结果，利用流域纬度加权方法计算流域水储量变化[11]：
| $ \Delta {\sigma _{{\rm{region}}}} = \frac{{\sum\limits_{\theta =-\frac{\pi }{2}}^{\frac{\pi }{2}} {\sum\limits_{\mathit{\emptyset} = 0}^{2\pi } {\Delta H\left( {\theta, \mathit{\emptyset} } \right) \cdot \cos \theta } } }}{{\sum\limits_{\theta =-\frac{\pi }{2}}^{\frac{\pi }{2}} {\sum\limits_{\mathit{\emptyset} = 0}^{2\pi } {\cos \theta } } }} $ | (5) |
以长江流域为例，利用GLDAS-NOAH模型计算流域平滑核函数与流域纬度加权，得到流域水储量变化(图 1)。由图 1可见，两种方法得到的结果非常接近，故后文只分析利用流域平滑核函数得到的流域水储量变化结果。
|
图 1 GLDAS-NOAH数据计算的长江流域水储量变化 Fig. 1 GLDAS-NOAH derived TWS changes of Yangtze river basin |
GLDAS-NOAH模型输入数据包括降水分析结果、太阳辐射、地表气压、湿度和NCEP(national centers of environmental prediction)提供的地表水平风速数据，输出结果包括土壤温度以及4层(从地表到2 m深)土壤水含量[6]。为与GRACE结果进行对比，本文采用空间分辨率为1°×1°的GLDAS-NOAH数据，该数据只包含模型输出的土壤水分和地表积雪变化，地下水变化信息不在其内。
为合理地比较GRACE得到的陆地水储量变化结果与水文模型的结果，对GLDAS-NOAH数据进行与GRACE类似的处理，即将GLDAS-NOAH格网数据球谐展开，截取至60阶。利用P3M8去相关滤波和高斯滤波，对展开得到的球谐系数进行处理，得到由GLDAS-NOAH计算的滤波后陆地水变化，再结合流域平滑核函数得到由GLDAS-NOAH计算的滤波后流域水储量变化σfiltered。
在GRACE数据后处理中会引入一些误差，如球谐系数的截断误差、高斯滤波和去相关滤波所引起的信号泄露误差。利用尺度因子可减小由GRACE数据后处理引入的误差。首先利用GLDAS-NOAH原始格网数据，得到流域陆地水变化序列σoriginal，然后利用尺度因子(scaling factor)k，使σoriginal与kσfiltered的残差平方和最小，最后将GRAEC得到的流域水储量乘以k，即可削弱数据的后处理误差。计算公式为[12]：
| $ \min = \sum\limits_{i = 1}^n {{{\left( {{\sigma _{{\rm{original}}}}-k{\sigma _{{\rm{filtered}}}}} \right)}^2}} $ | (6) |
计算得到四大流域的尺度因子分别为：长江0.986、黄河1.036、淮河1.776、珠江1.328，具体数值参见表 1。对比可知，流域面积越大，尺度因子越接近1。
|
表 1 四大流域水储量变化拟合结果 Tab. 1 Fitting results of the four basins of TWS |
图 2是GRACE和GLD AS-NOAH流域水储量变化序列的对比(GRACE的值是乘以尺度因子后的结果)。可以看出，GRACE与GLDAS-NOAH计算得到的等效水高序列呈现出很强的相关性，4个流域的相关系数分别为0.84、0.61、0.87、0.89。
|
图 2 四大流域GRACE计算陆地水储量和GLDAS-NOAH结果对比 Fig. 2 The comparison between GRACE derived TWS of the four basins and GLDAS-NOAH results |
2006年黄河流域出现严重旱情，从GRACE和GLDAS-NOAH序列上可反映出来，利用GRACE计算2006年黄河流域水储量较2005年减少了79.51%。
长江流域GRACE、GLDAS-NOAH时变序列呈现很强的相关特性。GRACE结果比GLDAS-NOAH振幅要大，因为GLDAS-NOAH只提取了地表土壤水和雪水变化结果。
淮河流域面积2.7×105 km2，GRACE和GLDAS-NOAH得到的流域水储量变化序列周年特性没有其他流域明显，其序列更加符合流域水储量的整体变化趋势。
2005年和2008年夏季，珠江流域降水量明显高于多年均值，多处出现洪水。从GRACE流域水储量序列也可看出，2008年珠江流域水储量高于多年均值。
3.2 2010年四大流域水储量变化2010年受厄尔尼诺现象的影响，我国南方出现严重的冬旱夏涝现象。图 3(b)是利用GRACE数据，结合350 km高斯滤波与去相关滤波(P3M8)进行处理，并利用尺度因子恢复泄露信号得到的2010年中国陆地水储量变化。可以看出，2010-01~04广西、云南、贵州、重庆等地陆地水储量明显偏低。据水利部门的资料显示，2009和2010年广西、云南、贵州和重庆都处于枯水年，2010-01~05，云南、贵州和四川出现几十年一遇的大旱，直到2010-05下旬旱情才有所缓解，这与李琼等[13]的结果一致。
|
图 3 2010年四大流域水储量变化 Fig. 3 TWS changes of the four basins in 2010 ①~④分别为长江流域、黄河流域、淮河流域和珠江流域 |
2010-06起，长江及珠江中下游湖北、安徽、江浙一带、两广地区等地陆地水储量急剧增加，2010-07~09珠江下游、广东等地、长江中下游、湖北、湖南等地出现特大洪水。利用尺度因子进行泄露信号恢复后，灾情严重情况更加明显，特别是淮河流域陆地水增加量明显增大。
可以看出，经过尺度因子的恢复，各大流域的水储量变化得到部分恢复，由于尺度因子随着流域面积的增大而趋于1，所以长江和黄河流域变化不是很明显。利用恢复后的流域平均水处理变化与没有利用尺度因子恢复的信号之差，乘以流域面积，得到由滤波处理导致的水储量泄露量。长江、黄河、淮河和珠江流域恢复的泄露量分别约为-541.1 m3、-142.1 m3、3.31×103 m3、1.03×103 m3。
3.3 流域水储量时间序列分析 3.3.1 多项式拟合对GRACE计算得到的等效水高时序进行拟合，得到周年项、半周年项和趋势项序列：
| $ S = {a_0} + {a_1}t + {a_2}\cos \left( {{\omega _1}t + {\varphi _1}} \right) + {a_3}\cos \left( {{\omega _2} + {\varphi _2}} \right) $ | (6) |
式中，a1为趋势项，a2为周年项系数，a3为半周年项，ω1=2π，ω2=4π。拟合得到的结果见表 1，表中相位延迟是GRACE与GPCP降水数据间的相位延迟，通过使两个时序相关系数最大得到。可以看出，长江、淮河和珠江流域的水储量变化呈增加趋势，而黄河流域呈减少趋势。
3.3.2 改进的EEMD算法利用随机噪声分布的功率谱密度也具有随机性，结合双样本KS检验验证分解得到的序列是否具有随机性[14]。在利用CEEMD分解信号时，需要对每一个生成的IMF进行KS2检验。具体步骤为：1)利用CEEMD分解方法，得到IMF1；2)利用FFT得到IMF1的功率谱密度PSDIMF1；3)利用双样本KS检验，检验PSDIMF1与随机分布噪声功率谱是否来自同一分布；4)如果检验通过，则说明该IMF1具有随机分布特性，继续利用CEEMD算法提取下一个IMF，重复步骤1)~4)，直到该IMF检验不通过，则直接采取EMD算法提取自该IMF以下的IMF。
以长江流域为例，利用改进的EEMD方法对GRACE和GLDAS陆地水储量序列进行分解，结果见图 4。其中，实线代表GRACE水储量变化及其分解结果，虚线代表GLDAS水储量变化及其分解结果，GLDAS水储量变化分解结果只包含4个IMF，没有出现同GRACE分解结果中IMF3特性的IMF，原因有待进一步研究。图 4中，IMF1具有随机特性，将其作为噪声剔除。可以看出，IMF2~IMF4(GLDAS对应IMF2~IMF3)具有不同的周期特性。一般来说，由于水负荷导致的形变具有的是周年特性，故本文将具有周年特性的IMF2作为水储量周年变化序列，IMF3~IMF4(GLDAS对应IMF3)为年际变化，IMF5(GLDAS对应IMF4)作为流域水储量变化趋势。
|
图 4 利用改进的EEMD算法对长江流域TWS的分解结果 Fig. 4 Modified EEMD in decomposition TWS changes of Yangtze river basin |
从图 5可以看出，黄河流域在2005年为丰水年，2006~2007年初进入枯水年，流域水储量较2005年减少约79.51%。期间，黄河整体水储量持续减少，这与水文模型的观测结果一致。通过观察图 4中的年际变化，可以得到流域的丰水和枯水循环周期：长江为3.05 a、黄河为2.37 a、淮河为2.37 a、珠江为2.13 a，这与实际的水循环周期基本一致[15]。对于流域水储量的趋势，长江和淮河流域保持着水储量缓慢增加的趋势，黄河流域为先快后慢的变化趋势，珠江则在2005~2008初呈现增加趋势，而后呈现减少趋势。
|
图 5 改进的EEMD算法在四大流域GRACE数据得到的TWS序列中的应用 Fig. 5 Modified EEMD applied in the GRACE derived TWS time series analysis for four basins |
图 6是利用改进的EEMD算法得到的四大流域GRACE与GLDAS陆地水储量周年和半周年项与年际变化，其中，红线和蓝线分别是GRACE、GLDAS陆地水储量周年和半周年项，绿线和黑线分别是GRACE、GLDAS陆地水储量年际变化(丰水、枯水循环周期)。可以看出，GRACE与GLDAS中都存在接近的水循环周期，GLDAS得到的4个流域丰水枯水循环周期分别为：长江3 a、黄河2.1 a、淮河2.2 a、珠江1.98 a，这与GRACE得到的丰水枯水循环周期相似。
|
图 6 利用改进的EEMD方法分别提取GRACE与GLDAS得到的4个流域TWS序列的周年和半年项及其年际变化 Fig. 6 Modified EEMD applied in extracting annual, semiannual and interannual components from GRACE and GLDAS derived TWS time series for four basins |
1) 利用GRACE得到的长江、黄河、淮河和珠江流域的水储量变化与全球水文模型(GLDAS-NOAH)得到的结果具有很强的相关性，相关系数分别为0.84、0.67、0.81、0.89。
2) 利用改进的EEMD算法对长江、黄河、淮河和珠江流域的GRACE水储量变化时序进行分析，得到4个流域的周年变化、年际变化和趋势项。结果表明，4个流域均具有丰水、枯水循环周期，长江为3.05 a、黄河为2.37 a、淮河为2.37 a、珠江为2.13 a。改进的EEMD在数据处理上具备自适应性特性，使其可以用于陆地水储量变化研究，如果对其分解得到的本征模态一一解释，会捕捉到其他地球物理信号。
3) 利用尺度因子可以恢复部分GRACE数据后处理泄露的信号，使得水储量反演结果更加接近真实值。尺度因子主要受流域大小影响，流域越小，得到的尺度因子也越大，信号恢复也就越真实。
| [1] |
Huang N E, Sheng Z, Long S R, et al. The Empirical Mode Decomposition and the Hillbert Spectrum for Nonlinear and Non-Stationary Time Series Analysis[J]. Proc Roy Soc London A: Mathematical, Physical and Engineering Sciences, 1998, 454(1971): 903-995 DOI:10.1098/rspa.1998.0193
(  0) |
| [2] |
Loomis B D, Luthcke S B. Optimized Signal Denoising and Adaptive Estimation of Seasonal Timing and Mass Balance from Simulated GRACE-Like Regional Mass Variations[J]. Advances in Adaptive Data Analysis, 2014, 6(1)
( 0) |
| [3] |
Wu Z H, Huang N E. Ensemble Empirical Mode Decomposition: A Noise Assisted Data Analysis Method[J]. Advance in Adaptive Data Analysis, 2009, 1(1): 1-41 DOI:10.1142/S1793536909000047
( 0) |
| [4] |
Yeh J R, Shieh J S. Complementary Ensemble Empirical Mode Decomposition:A Noise Enhanced Data Analysis Method[J]. Advance in Adaptive Data Analysis, 2010, 2(2): 135-156 DOI:10.1142/S1793536910000422
( 0) |
| [5] |
郑近德, 程军圣, 杨宇. 改进的EEMD算法及其应用研究[J]. 振动与冲击, 2013, 32(21): 21-26 (Zheng Jinde, Cheng Junsheng, Yang Yu. Modified EEMD Algorithm and Its Applications[J]. Journal of Vibration and Shock, 2013, 32(21): 21-26 DOI:10.3969/j.issn.1000-3835.2013.21.004)
( 0) |
| [6] |
Rodell M, Houser P R, Jambor U, et al. The Global Land Data Assimilation System[J]. Bull Amer Meteor Soc, 2004, 85(3): 381-394 DOI:10.1175/BAMS-85-3-381
( 0) |
| [7] |
Swenson S, Wahr J. Methods for Inferring Regional Surface-Mass Anomalies from Gravity Recovery and Climate Experiment(GRACE)[J]. J Geophys Res, 2002, 107(B9)
( 0) |
| [8] |
Bettadpur S. UTCSR Level-2 Gravity Field Product User Handbook, GRACE 327-734[Z]. Center for Space Research, the University of Texas Austin, Austin, 2007
( 0) |
| [9] |
Cheng M K, Tapley B D. Variations in the Earth's Oblateness during the Past 28 Years[J]. J Geophys Res, 2004, 109(B9)
( 0) |
| [10] |
Swenson S, Chambers D, Wahr J. Estimating EatimatingGeocenter Variations from a Combination of GRACE and Ocean Model Output[J]. J Geophys Res, 2008, 113(B8)
( 0) |
| [11] |
Chen J L, Wilson C R, Famiglietti J S, et al. Attenuation Effect on Seasonal Basin-Scale Water Storage Changes from GRACE Time-Variable Gravity[J]. Journal of Geodesy, 2007, 81(4): 237-245 DOI:10.1007/s00190-006-0104-2
( 0) |
| [12] |
Landerer F W, Swenson S C. Accuracy of Scaled GRACE Terrestrial Water Storage Estimates[J]. Water Resour Res, 2012, 48(4)
( 0) |
| [13] |
李琼, 罗志才, 钟波, 等. 利用GRACE时变重力场探测2010年中国西南干旱陆地水储量变化[J]. 地球物理学报, 2013, 56(6): 1843-1849 (Li Qiong, Luo Zhicai, Zhong Bo, et al. Terrestrial Water Storage Changes of the 2010 Southwest China Drought Detected by GRACE Temporal Gravity Field[J]. Chinese Journal of Geophysics, 2013, 56(6): 1843-1849 DOI:10.6038/cjg20130606)
( 0) |
| [14] |
Wouters B, Schrama E J O. Improved Accuracy of GRACE Gravity Solutions through Empirical Orthogonal Function Filtering of Spherical Harmonics[J]. Geophys Res Lett, 2007, 34(23): 229-241
( 0) |
| [15] |
王金花, 刘红梅, 康玲玲, 等. 黄河中游干旱的变化及区间遭遇分析[J]. 干旱区资源与环境, 2006, 20(6): 109-113 (Wang Jinhua, Liu Hongmei, Kang Lingling, et al. Analysis on the Changes of Drought in the Middle Reaches of Yellow River[J]. Journal of Arid Land Resources and Enviroment, 2006, 20(6): 109-113)
( 0) |