2018, Vol. 38
2. 高速铁路运营安全空间信息技术国家地方联合工程实验室，成都市高新区西部园区，611756;
3. 四川测绘地理信息局测绘技术服务中心，成都市人民北路198号，610081
GRACE重力卫星由美国NASA和德国DLR合作发射，用来探测地球重力场的时空变化。Wahr等[1]提出了利用GRACE时变重力场模型估计地表质量变化的算法。GRACE观测数据能够每月解算一个空间分辨率达到400 km、精度达到1 cm等效水高的时变重力场模型[2]，即GRACE可以以一定的精度来研究和分析大尺度水文信号变化，但是重力场模型的中高阶项会产生大量的条带噪声，对其处理有不结合和结合球谐位系数先验精度信息两种方法。结合先验精度信息方法建模较为复杂，本文不作讨论；不结合先验精度信息方法一般是对球谐位系数直接进行频谱滤波，主要包括各向同性高斯滤波、各向异性扇形滤波、去相关滤波等，这些滤波在剔除部分条带误差的同时会抑制真实的信号而造成信号的泄漏[3]。信号提取方法主要是多项式拟合、主成分分析和独立成分分析等[4]。
主成分分析(principal component analysis，PCA)又称为经验正交函数方法(empirical orthogonal function, EOF)，是在20世纪初由Pearson和Hotlling两位学者各自研究并独立推广而来，当时其主要任务是用来求解空间和时间权重系数[5]。近年来，主成分分析被用于GRACE时变地球重力场信号的提取，主要包括3种策略[6-8]：1)直接对球谐位系数作主成分分析，能够最大程度地提取时变信号，但其提取的信号和实验前预期的地球物理信号相关性很小；2)先对不同月份的球谐位系数滤波，再计算等效水高并作主成分分析，进而提取全球水文变化的长期趋势，但由于在主成分分析前对所有阶次的球谐位系数使用了相同的滤波参数，导致噪声剔除和信号保留不能兼顾，达不到理想的效果；3)对不同机构同一时间段的等效水高作主成分分析，提取同一时间段各机构时变模型的共同信号，但不能用于长期趋势的研究。这表明主成分分析能有效地应用于时变重力场信号的提取，达到信号和噪声分离的效果。
Wouters等[4]首次提出了对球谐位系数进行主成分分析，相比其他的滤波方法，该方法提取的信号没有衰减趋势。Frederic等[6]利用GRACE重力场模型基于EOF方法提取了南美洲2003~2010年等效水高在空间域和时间域的变化，研究近些年干旱、洪涝等极端自然灾害的变化规律；柯宝贵等[7]对2002~2010年GRACE重力场模型反演的等效水高作主成分分析，提取了重力场变化的长期趋势；穆大鹏等[8]也对GRACE重力场模型反演的等效水高进行了主成分分析。在上述研究中对提取的各个主成分都是选取统一的滤波参数，没有考虑各个模态所含噪声不同的特点，故会造成噪声剔除不彻底和泄漏误差明显的结果。本文通过控制不同主成分球谐位系数反演全球等效水高的截断阶数，找出各主成分噪声特点，再根据这些特点选择最佳滤波参数，这样在剔除噪声的同时又最大程度地保留了真实信号。
1 基本原理 1.1 等效水高计算根据GRACE球谐位系数计算等效水高的计算式为[1]：
| $ \begin{array}{l} {\rm{EWH}}\left( {\theta, \lambda } \right) = \frac{{a{\rho _{{\rm{ave}}}}{\rm{}}}}{{3{\rho _{\rm{w}}}}}\sum\limits_{n = 0}^\infty {\frac{{2n + 1}}{{1 + {k_n}}}} {W_n}\sum\limits_{m = 0}^n {} \\ (\Delta {{\bar C}_{nm}}{\rm{cos}}m\lambda + \Delta {{\bar S}_{nm}}{\rm{sin}}m\lambda ){{\bar P}_{nm}}({\rm{cos}}\theta ) \end{array} $ | (1) |
式中，a、
本文时变重力场模型球谐位系数的截断阶数均取60，由于GRACE提供的球谐位系数是假定地球质心与坐标系原点重合的，没有考虑地球质心的移动，所以不考虑零阶项和一阶项；利用SLR测得的C2, 0项代替GRACE的C2, 0项[9]。将一个月的球谐位系数按照阶数和次数的顺序展开为一个3 717 ×1的列向量，再将108个月的球谐位系数构造为一个3 717×108的矩阵。根据主成分原理[6]，将108个月的球谐位系数作为样本X，表示如下：
| $ \mathit{\boldsymbol{X }} = \left\{ {\begin{array}{*{20}{c}} {\Delta {C_{2, 0(t = 1)}}}&{\Delta {C_{2, 0(t = 2)}}}&{ \ldots }&{\Delta {C_{2, 0(t = 108)}}}\\ {\Delta {C_{2, 1(t = 1)}}}&{\Delta {C_{2, 1(t = 2)}}}&{ \ldots }&{\Delta {C_{2, 1(t = 108)}}}\\ \vdots & \vdots & \vdots & \vdots \\ {\Delta {C_{60, 60(t = 1)}}}&{\Delta {C_{60, 60(t = 2)}}}& \ldots &{\Delta {C_{60, 60(t = 108)}}}\\ {\Delta {S_{2, 0(t = 1)}}}&{\Delta {S_{2, 0(t = 2)}}}&{ \ldots }&{\Delta {S_{2, 0(t = 108)}}}\\ {\Delta {S_{2, 1(t = 1)}}}&{\Delta {S_{2, 1(t = 2)}}}&{ \ldots }&{\Delta {S_{2, 1(t = 108)}}}\\ \vdots & \vdots & \vdots & \vdots \\ {\Delta {S_{60, 60(t = 1)}}}&{\Delta {S_{60, 60(t = 2)}}}&{ \ldots }&{\Delta {S_{60, 60(t = 108)}}} \end{array}} \right\} $ | (2) |
式中，t表示月份，将矩阵 X标准化后再计算其相关系数矩阵，最后计算出特征值λ和对应的特征向量 E，通过主成分分析得到新变量：
| $ \mathit{\boldsymbol{P}}{_N} = ({\rm{ }}\mathit{\boldsymbol{P}}{_1}, {\rm{ }}\mathit{\boldsymbol{P}}{_2}, \ldots, {\rm{ }}\mathit{\boldsymbol{P}}{_n}) = \mathit{\boldsymbol{ XE}} $ | (3) |
通过选取较大的k个特征值所对应的新变量PK=(P1, P2, …, Pk)重构原来的重力场模型，可以达到信息提取的目的[10]：
| $ \mathit{\boldsymbol{X}}{_{{\rm{PCA}}}} = {\rm{ }}\mathit{\boldsymbol{P}}{_K}\mathit{\boldsymbol{E}}{^{{\rm{T}}}} $ | (4) |
| $ \alpha = \frac{{\sum\limits_{i = 1}^k {{\rm{ }}{\lambda _k}} }}{{\sum\limits_{{\rm{ }}i = 1}^{108} {{\lambda _i}} }} $ | (5) |
式中，重构重力场模型变量XPCA所包含的信息可以由贡献率α来衡量。
根据GRACE提供的每月重力场模型通过式(1)计算全球格网(1°×1°)等效水高数据，按照经纬度顺序排列成具有360×180=64 800个元素的列向量(相当于把360×180的矩阵拉直成为一个列向量)，那么每个月的等效水高数据能形成一个空间场分量，108个月的等效水高数据为64 800×108的矩阵，对此矩阵进行主成分分析就是基于等效水高的主成分分析。
目前主要的时变重力场主成分分析法步骤如表 1所示。方法a和b原理一样，只是滤波和主成分分析顺序不同，信号的提取效果一样。由式(1)计算的等效水高是球谐位系数通过线性变换得到，故方法a和d本质上也一样。对上述两种数据分别进行球谐位系数和等效水高主成分分析，算出两种方法各模态的贡献率和累计方差贡献率，如图 1所示。可以看出，两者的各模态贡献率和累计方差贡献率基本相同，说明两种方法结果基本一致。本文主要探讨改进方法c相对于其他3种方法的优点，由于其他3种信号提取结果基本一致，故只要c优于任意一种方法即可。
|
表 1 主成分分析步骤 Tab. 1 The steps of principal component analysis |
|
图 1 两种方法各模态贡献率和累计方差贡献率 Fig. 1 Two methods' each mode contribution rate and cumulative variance contribution rate |
由图 1可知，前三模态累计贡献率接近80%，体现了绝大部分信号，故只讨论前三模态。利用方法c对数据进行主成分分析，得到前三模态的球谐位系数，再根据式(1)计算等效水高，如图 2所示。
|
图 2 球谐系数主成分前三空间和时间模态 Fig. 2 First three spatial and temporal modes of PCA based on spherical harmonic coefficients |
由图 2可以看出，前三空间模态都存在大量的南北条带噪声。第一模态中，强信号主要出现在巴西北部亚马逊流域、非洲中部刚果河流域、亚洲中西部湄公河流域；第二模态中，强信号主要出现在格陵兰岛、南极、北美东北部地区以及阿拉斯加冰川区域；第三模态中，强信号主要出现在巴西北部亚马逊流域。第一和第三时间模态均呈年周期的变化，第一时间模态在整个时间段内信号变化趋于稳定，振幅较大；第二时间模态在2005~2013年有逐年下降的趋势；第三时间模态信号同样没有明显的突变趋势，年周期变化明显。
时间模态的周年性变化与全球的季节性气候变化密切相关。亚马逊流域、刚果河流域、亚洲南部等处于赤道附近，旱季和雨季分明，降水量受季节影响非常大，地表水量变化会引起重力场变化，所以有较强的信号。由此可以总结出，空间模态表现出明显的区域特征，这些地区大多是降雨量受季节变化、冰川融化等影响，结合时变重力场产生的原因与这些地方有强信号也是相符的[11]。
2.2 球谐位系数主成分分析的改进由图 2发现，随着模态数的递增，所含的信号越来越少，噪声越来越多，且各模态所包含的南北条带噪声也不相同，第一模态噪声主要集中在(30°N~50°N、30°S~50°S)，第二模态噪声分布比较均匀，第三模态噪声较大，分布主要集中在(50°S~50°N)。对前三模态作P5M10去相关滤波和250 km高斯滤波处理，反演等效水高并把结果作为真值；计算前三模态的信噪比分别为2.669、2.164、0.721。由于各种模态噪声含量和分布均不相同，如果对不同模态仍采用相同的滤波方式，显然是不合理的，对此提出改进：针对不同模态所含噪声的特点选取不同的高斯滤波参数。
对第一模态，分别以5、10、15、…、60为截断阶数来计算等效水高，如图 3所示。从图 3可以看出，在25阶和30阶时表现为强信号，主要噪声从35阶开始；故对第一模态35阶以上的项进行高斯滤波，35阶以下的项不处理。最佳滤波半径的确定：对第一模态(仅处理35阶以上的项)采用半径为200 km、250 km、300 km、…、1 000 km的高斯滤波处理；结合上文中第一模态等效水高真值计算出不同半径高斯滤波等效水高的信噪比，如图 4中(a)所示。可见，当半径为350 km时信噪比最高为9.674，故对第一模态，从35阶开始采用半径为350 km的高斯滤波。
|
图 3 第一模态不同阶数时的等效水高 Fig. 3 Equivalent water pattern of the first modal with different orders |
|
图 4 三种模态不同滤波半径的信噪比 Fig. 4 Signal to noise ratio of different filter radius in three modes |
对第二和第三模态同样以5、10、15、…、60为截断阶数反演等效水高，来分析噪声特点。由于篇幅问题，不再给出等效水高图，第二模态和第三模态主要噪声分别从30阶、15阶开始。故对第二模态30阶以上的项进行高斯滤波，30阶以下的项不处理；对第三模态15阶以上的项进行高斯滤波，15阶以下的项不处理。同上，结合上文中第二模态、第三模态等效水高真值，计算第二模态(仅处理30阶以上的项)、第三模态(仅处理15阶以上的项)不同半径高斯滤波等效水高的信噪比，如图 4中(b)、(c)所示。可见，当滤波半径为300 km时第二模态信噪比最高为7.415，半径为500 km时第三模态信噪比最高为3.480。故对第二模态，从30阶开始采用半径为300 km的高斯滤波；对第三模态从15阶开始，采用半径为500 km的高斯滤波。
为了进行对比，引入两种传统的主成分分析：2005-01~2013-12(108个月)的球谐位系数分别采用半径为250 km和350 km的高斯滤波处理，计算等效水高，再进行主成分分析(在后文中简称方法A和方法B)。将3种方法提取的前三模态进行比较，如图 5所示。
|
图 5 3种方法前三主成分比较 Fig. 5 Compare the three main components of the three methods |
图 5中3种方法前三模态包含的强信号几乎一样，对比发现图 5(a)、(d)、(g)所含的南北条带噪声明显较多，可知方法A对噪声的剔除效果不明显。比较图 5(b)、(c)、(h)、(i)未发现明显差异，图 5(e)、(f)在格陵兰岛部分地区和南极差距较大，这表明方法B和改进法提取的第一模态、第三模态几乎相同，第二模态在格陵兰岛部分地区和南极差距较大。为了进一步探讨这种差距把方法A、B得到的前三模态与改进法得到的前三模态作差(图 6)。
|
图 6 方法A、B与改进法前三模态之差 Fig. 6 First three modes' differences between methods A, B and improved method |
图 6(a)、(c)、(e)几乎没有强信号，说明改进法对真实信号的保留效果与方法A相当，即改进法提取的前三模态信号没有大量泄漏。图 6(a)、(c)、(e)存在大量南北条带噪声，再次证明了改进法对南北条带噪声的剔除效果优于方法A。图 6(b)、(d)、(f)均未体现出南北条带噪声，这说明改进法对南北条带噪声的剔除效果与方法B相当。图 6(d)强信号主要集中在格陵兰岛和南极，结合图 5(d)、(e)、(f)以及图 6(c)得到如下结论：方法B的第二模态在格陵兰岛和南极泄漏现象明显，而改进法的第二模态在格陵兰岛和南极信号保留较好，泄漏误差较小。总的来说，改进法对南北条带噪声的剔除效果与方法B相当，对真实信号的保留效果与方法A相当。
实验结果表明，未改进的主成分分析法在信号保留和噪声剔除两方面不能兼顾：
1) 大半径的高斯滤波主成分分析法(如方法B)对南北条带噪声剔除效果明显，但高纬度地区的信号会被抑制从而造成信号泄漏。
2) 小半径的高斯滤波主成分分析法(如方法A)在高纬度地区泄露误差较小，但对南北条带噪声的剔除效果不明显。
改进后的方法考虑了各个模态噪声的特点，对不同模态采用不同的滤波参数，一定程度上解决了这个问题。相对于传统方法，改进后的方法不仅对中低纬度地区南北条带噪声的剔除效果明显，而且格陵兰岛和南极地区信号得到最大程度的保留，泄漏误差大大减小。
2.3 对比分析根据图 1知基于球谐位系数主成分分析和基于等效水高主成分的前7模态累计方差贡献率分别为0.861 3，0.859 7。以85%为重构阈值(前7模态)根据式(4)用改进法(第4至第7模态最佳滤波参数见表 2)和方法B，对2006-06和2006-12的等效水高数据进行重构，如图 7(c)、(d)、(e)、(f)所示。将2006-06和2006-12的GRACE重力场模型进行高斯250 km滤波再计算等效水高，如图 7(a)和(b)所示，其在高纬度地区的等效水高可以近似看作真值。为了对比两种方法重构结果的差异，把方法B与改进法的重构结果作差，如图 8所示。
|
图 7 原始数据小半径高斯滤波、方法B重构、改进法重构比较 Fig. 7 Original data small radius Gaussian filtering, method B reconstruction and improved method reconstruction comparison |
|
|
表 2 第4至第7模态改进法高斯滤波参数 Tab. 2 The fourth to seventh modal improvement method Gaussian filter parameters |
|
图 8 方法B与改进法重构之差 Fig. 8 Differences between method B reconstruction and improved reconstruction |
图 7(a)、(b)在高纬度部分地区的信号保留较好，但是在50°N~50°S之间条带噪声十分明显；图 7(c)、(e)，(d)、(f)中低纬度地区没有条带噪声，可以看出方法B和改进法在中低纬度地区对南北条带噪声的剔除效果都比较明显。图 8(a)、(b)中低纬度地区无强信号，值均接近于零，这说明两种方法对噪声的剔除效果相当；图 8中的强信号主要集中在格陵兰岛部分区域和南极。为了进一步比较两种方法重构结果在格陵兰岛地区的差异，计算格陵兰岛2006-06和2006-12原始数据未经滤波处理的等效水高，方法B、改进法提取的等效水高如图 9所示。
|
图 9 3种方法格林兰岛的等效水高 Fig. 9 Equivalent water high of three methods in Greenland |
从图 9可以看出，在格陵兰岛中部、西南部、东南部的部分地区原始数据计算等效水高为负值，而方法B重构结果变为正值造成信号失真。改进后的方法在一定程度上改善了这种状况，信号保留较好，泄漏误差大大减小。改进法重构结果更加接近于原始数据反演的等效水高，泄漏误差明显小于方法B。
综合上述结果可以得出，改进后的方法在中低纬度地区的噪声剔除效果和方法B几乎一致，在高纬度地区泄漏误差大大减小，故改进后的方法优于传统的主成分分析法。
3 结语本文基于球谐位系数进行主成分分析，根据前三个时间模态和空间模态分析出部分区域等效水高呈现出周年性变化，这与全球的季节性气候变化密切相关。
基于本文提出的改进思想，通过控制各主成分球谐位系数的截断阶数来分析各个主成分含噪声的特点，并结合信噪比选取最佳滤波半径，将改进法的前三模态与两种传统主成分分析法的前三模态进行比较，发现改进后的方法不仅在低纬度南北条带噪声的剔除效果明显，而且在格陵兰岛和南极真实信号保留较好，泄漏误差大大减小。
通过改进前后两种方法重构2006-06和2006-12的等效水高，并与原始球谐位系数通过小半径高斯滤波反演的等效水高相比较，结果表明改进后的主成分分析法不仅在低纬度地区的噪声剔除效果明显，而且在格陵兰岛地区泄漏误差明显小于改进前的主成分分析，说明改进后的主成分分析优于传统的主成分分析。
| [1] |
Wahr J, Molenaar M, Bryan F. Time Variability of the Earth's Gravity Field: Hydrological and Oceanic Effects and this Possible Detection Using GRACE[J]. Journal of Geophysical Reseach, 1998, 103(B12): 30 205-30 229 DOI:10.1029/98JB02844
(  0) |
| [2] |
Wahr J, Swenson S, Zlotnicki V, et al. Time-Variable Gravity from GRACE: First Results[J]. Geophysical Research Letters, 2004, 31(11)
( 0) |
| [3] |
许才军, 龚正. GRACE时变重力数据的后处理方法研究进展[J]. 武汉大学学报:信息科学版, 2016, 41(4): 503-511 (Xu Caijun, Gong Zheng. Review of the Post-Treatment Methods on GRACE Time Varyed Gravity Data[J]. Geomatics and Information Science of Wuhan University, 2016, 41(4): 503-511)
( 0) |
| [4] |
Wouters B, Schrama E J O. Improved Accuracy of GRACE Gravity Solutions through Empirical Orthogonal Function Filter of Spherical Harmonic[J]. Geophysical Research Letters, 2007, 34(23)
( 0) |
| [5] |
Frappart F, Seoane I, Ramillen G. Validation of GRACE-Derived Terrestrial Water Storage from a Regional Approach over South America[J]. Remote Sensing of Environment, 2013, 137: 69-83 DOI:10.1016/j.rse.2013.06.008
( 0) |
| [6] |
夏非, 张永战, 吴蔚. EOF分析在海岸地貌与沉积学研究中的应用进展[J]. 地理科学进展, 2009, 28(2): 174-186 (Xia Fei, Zhang Yongzhan, Wu Wei. Progress in Application of EOF Analysis in Reseach of Coastal Geomorphology and Sedimentology[J]. Progress in Geography, 2009, 28(2): 174-186 DOI:10.11820/dlkxjz.2009.02.003)
( 0) |
| [7] |
柯宝贵, 章传银, 张利明, 等. 应用经验正交函数分析方法研究[J]. 科学技术与工程, 2012, 12(2): 275-279 (Ke Baogui, Zhang Chuanyin, Zhang Liming, et al. Analysis on the Earth Gravity Temporal Variation Based on the Orthogonal Function Method[J]. Science Technology and Engineering, 2012, 12(2): 275-279)
( 0) |
| [8] |
穆大鹏, 郭金运, 孙中昶, 等. 基于主成分分析的GRACE重力场模型等效水高[J]. 地球物理学进展, 2014, 29(4): 1 512-1 517 (Mu Dapeng, Guo Jinyun, Sun Zhongchang, et al. Equivalent Water Height form GRACE Gravity Model Based on Principal Component Analysis[J]. Progress in Geophysics, 2014, 29(4): 1 512-1 517)
( 0) |
| [9] |
Cheng M K, Ries J C, Tapley B D. Geocenter Variations from Analysis of SLR Data[C]. IAG Symposia 2010, Reference Frames for Application in Geosciences(REFAG2010), Marne-La-Vallee, 2010
( 0) |
| [10] |
于红娟, 郭金运, 李九龙, 等. 基于主成分分析的相对重力观测零漂和固体潮提取[J]. 大地测量与地球动力学, 2015, 35(5): 906-910 (Yu Hongjuan, Guo Jinyun, Li Jiulong, et al. Zero Drift and Solid Earth Tide Extracted from Relative Gravimetric Data with Principal Component Analysis[J]. Journal of Geodesy and Geodynamics, 2015, 35(5): 906-910)
( 0) |
| [11] |
李琼, 罗志才, 钟波, 等. 利用GRACE时变重力场探测2010年中国西南干旱陆地水储量变化[J]. 地球物理学报, 2013, 56(6): 1 843-1 849 (Li Qiong, Luo Zhicai, Zhong Bo, et al. Terrestrial Water Storage Changes of the 2010 Southwest China Drought by GRACE Temporal Gravity Field[J]. Chinese Journal of Geophysics, 2013, 56(6): 1 843-1 849)
( 0) |
2. State-Province Joint Engineering Laboratory of Spatial Information Technology for High-Speed Railway Safety, West High-Tec Zone, Chengdu 611756, China;
3. Surveying and Mapping Technical Service Center of Sichuan Bureau of Surveying Mapping and Geoinformation, 198 North-Renmin Road, Chengdu 610081, China