Документ взят из кэша поисковой машины. Адрес оригинального документа : http://lnfm1.sai.msu.ru/~tempus/english/science/articles/ZotovShum.pdf Дата изменения: Tue May 25 16:33:36 2010 Дата индексирования: Mon Oct 1 23:39:32 2012 Кодировка: Поисковые слова: universe

Multichannel singular spectrum analysis of the gravity field data from GRACE satellites.
Zotov L.V. and Shum C.K.


Sternberg Astronomical Institute, Moscow State University, Universitetsky pr. 13, Moscow, Russia 119992 (wolftempus@gmail.com). School of Earth Sciences, Ohio State University, 275 Mendenhall Laboratory 125 South Oval Mall, Columbus OH 43210-1308 USA.

Abstract. The NASA/DLR GRACE (Gravity Recovery and Climate Experiment) twin-satellite mission was launched on 17.03.2002. They provide data products including monthly gravitational field solutions, which reflect the mass redistribution within the Earth system. However, these data need filtering and reduction of the so-called stripes or correlated high-frequency errors due primarily to model errors and imperfect observability of the gravity field inversion solution. We applied the multichannel singular spectrum analysis (MSSA) method to the post-process monthly GRACE data in the spectral and spatial domains. This method allows an improved separation of signals represented by the principal components containing periodic and secular terms, which can be associated with seasonal or longer redistribution of water masses, ice melting, glacial isostatic adjustment, ocean mass variations, and transient events such as the Sumatra earthquake in 2004. Flexibility of MSSA, which allows to filter out stripes and high-frequency noise, makes the method useful for GRACE data post-processing. Keywords: GRACE, temporal gravity field of the Earth, singular spectrum analysis PACS: 91.10.Fc; 91.10.Op; 91.10.Xa

INTRODUCTION
The GRACE twin-satellite mission which is developed jointly by American and German scientists was launched from Plesetsk cosmodrome in 2002. They have a mean orbital altitude of about 500 km, one following another with a mean separation distance of 220 km, and with near-polar orbits. They are equipped with a dual-frequency K-band microwave inter-satellite ranging system, GPS receivers, accelerometers, USOs, star trackers and laser corner cube retro-reflectors. The fundamental measurement is the distance between the satellites or its changes under the influence of accelerations due to the attracting masses GRACE overflights and due to other causes. The measured range or range-rate is the so-called Level 1 data product [2]. These data are then used along with numerous forward models accounting for the "known" forces to solve for the Earth's gravity field at a monthly temporal resolution. This process involves precision orbit determination with use of on-board GPS, non-gravitational modeling using the accelerometers measurements, correction for ocean and earth tides, as well as for barotropic responses of the ocean and solid earth to the atmospheric pressure changes. The results of processing form the Level 2 data product [1]. This is the monthly set of Stokes coefficients representing the Earth's gravitational potential complete to spherical harmonic degree 60 or 120, depending on the origin of various data centers.


FIGURE 1. The map of the first principal component of the difference between the observed gravitational field and GGM01C model in EWH, obtained by the MSSA method. PC 1 contains most of the stripes energy. Initial and final day of the year of the measurements epoch is shown above.

Our work is directed at providing an alternative method to reduce the so-called stripes or high-frequency, geographically correlated errors in the GRACE equivalent water heights (EWH) monthly maps (Figure 1) and separating principal components (PCs).

INITIAL DATA AND PROCESSING METHOD
The static or mean part of Earth's gravitational field is available from the models such as GGM01C, GGM02C, EGM2008 (up to degree 2160) [9], constructed using the data from CHAMP, GRACE and other satellites, altimetry, terrestrial and marine gravimetry. The variable component can be obtained by subtraction the model (GGM01C was used) from the monthly GRACE Level 2 data, calculated by one of the processing centers (we used the CSR RL04 Level 2 data). These datasets include Stokes coefficients Cnm , Snm representing the exterior geopotential spherical harmonics expansion V ( , , r) = GM r

n=2 m=0





n

a r

n

m (Cnm cos m + Snm sin m )Pn (sin ),

(1)

m where n is the degree, m is the order of the spherical harmonic, Pn are the fully normalized Associated Legendre Polynomials, as well as a is the mean equatorial radius of the Earth, the arguments , , r are latitude, longitude and radius respectively [12]. We used the data complete to degree 60. The coefficients of zero and first degree are zero due to the choice of the coordinate system and that GRACE is insensitive to degree one coefficients (geocenter). We processed 79 months of data since 8/2002 through 2/2009, but some months were not complete (6/2003, 1/2004). Stokes coefficient changes in time were linearly interpolated. In 2002-2003 the satellites calibration/validation was performed, so the data is not as precise as obtained later in operational mode.


The difference Cnm (t ), Snm (t ) between the observed coefficients and the model is expressed in terms of the EWH anomalies: 2(n + 1) m Wn (Cnm (t ) cos m + Snm (t ) sin m )Pn (sin ), 1 + kn n=2 m=0 (2) here pave and pw are average densities of the Earth and sea water respectively, kn is the load Love number of degree n, Wn is a filter coefficient. Since our goal was to demonstrate that MSSA method allows to perform filtering of observations, we have taken Wn = 1. In spectral domain MSSA-processing was conducted before using the expression (2). However, different authors use a variety of filtering methods for elimination or reduction of noise and stripes in the GRACE monthly gravity field solutions. Among them ­ Gaussian filtering with symmetric and asymmetric Gaussian function [7], Wiener [10] and regularizing [11] filters, whose coefficients depend not only on degree, but on order, de-stripping/smoothing [15] filters, operating on removal of anomalously large values from the resonant orders of the Stokes' coefficients. MSSA can also be used with these types of filtering. MSSA method [3], [8], also known as an Extended Empirical Orthogonal Functions (EEOF) analysis, is a generalization of the singular spectrum analysis (SSA) over the multidimensional time series. SSA in its turn is a generalization of principal component analysis (PCA) and consists of four stages: (a) the formation of trajectory matrix, (b) its singular values decomposition (SVD), (c) principal components grouping and (d) their recovery through hankelization. SSA algorithm is described in details in [4],[5]. The main parameter of the algorithm is the caterpillar-SSA length L (lag), which determines the dimensionality of the time series embedding space and can be chosen up to onehalf of their data span. We have chosen after some experiments L = 18, which gives sufficiently good results in separation of components. Trajectory matrix can be build for one time series component (channel), selected from the Stokes coefficients matrix, let's say for Ci j (tk ), k = 0, . . . , N - 1 as follows Ci j (t0 ) Ci j (t1 ) . . . Ci j (tK -1 ) Ci j (t1 ) Ci j (t2 ) . . . Ci j (tK ) XCij = (3) . . . ... . . . . . . Ci j (tL-1 ) Ci j (tL ) . . . Ci j (tN -1 ) h( , , t ) = a pave 3 pw



60

n

here K = N - L + 1. Matrix X, composed of such blocks for every channel Ci j and Si X = [XC2,0 , XS2,0 ..., , XCi j , XSi j , ..., XC60,60 XS60,60 ]

j

(4)

then can be used to calculate the lagged covariance matrix A = XT X. PC's can be obtained from the solution of the eigenvalue problem for A. They can be also obtained directly from the SVD of the matrix X X = USVT , so that squares of it's singular values si , which stand on the diagonal of S, are eigenvalues of A, and left eigenvectors vi , included as columns in V, form empirical orthogonal


FIGURE 2.

Singular values, representing the signal energy distribution over MSSA components.

functions (EOFs). Then the component i corresponds to the matrix Xi = si ui vT , i and forms a PC (in general, several singular values could be grouped in one PC. We did it with use of -correlations [5]). In MSSA we should reconstruct the vectorial time series of PC from this matrix Xi . Hankelization allows to reconstruct PC in every channel from the corresponding blocks of matrix, organized as in 4. Suppose we need to reconstrust the Cl m channel of the PC i, let's call it g, then each k-th count can be obtained from the averaging along the side diagonals of the corresponding matrix block of Xi , let's call it i Y = XClm , as follows
1 k+1 1 L k+1

gk =

m=1 L

y

m,k-m+2

,

0 k < L - 1 L - 1 k < K ,


1 N -k

m=1 N -K +1

y

m,k-m+2

,

m=k-K +2



y

m,k-m+2

K k < N

where L = min(L, K ), K = max(L, K ). It can be seen, that the first and the last L elements are calculated from fewer number of matrix values, so the first and the last points of PCs can be less consistent. It is supposed, that elements on side diagonals of matrix Y are approximately equal, and it is almost Hankel. In case, when it doesn't hold strictly, some kind of border effect appears (see. Fig 2). When L = 1 SSA method becomes PCA and principal components are obtained in fact from the analysis of the usual non-lagged signal covariance matrix A. PCA by the name of EOF-analysis was applied to the GRACE data in works [13], [16],[17]. In [13] SSA was also tested after PCA on the PCs, containing annual signal. It showed good approximation abilities. In [6] under the same name a bit different and very promising method was applied ­ a non-linear modifications of EOF-analysis for non-stationary time series, where PC's are obtained by means of time series envelopes calculations and orthogonalisation. But the MSSA, despite its mathematical complexity, has significant benefits. In the matrix A mutual correlations of different channels present, MSSA's PCs contain correlations, which appear in all the channels simultaneously. MSSA can allocate time-frequency (in case of analysis of the Stokes coefficients matrixes in spectral


FIGURE 3. Approximation of the S21 Stokes coefficient behavior by the PCs after MSSA (L=18) processing in spectral domain.

FIGURE 4.

Second PC, representing annual hydrological cycle.

domain) or spatio-temporal (in case of analysis of digital longitude-latitude maps of the anomalies in the spatial domain) correlations, inherent in all the signal coordinate components (channels).

RESULTS OF PROCESSING
First we applied the MSSA in the spectral domain to the Stokes coefficients. After calculation the principal components were converted into spatial maps of EWH. The distribution of signal energy, depending on the singular values is represented by Fig. 2. Typical behavior of PCs and how they approximate initial signal on the example of S21 Stokes coefficient is shown on Fig. 3. First PC appeared to have almost unchanging amplitude and geometry in time (Fig. 1 and 3 left). It contains most of the stripes energy and probably some constant component of the useful signal. This PC can be used for reducing the energy of the stripes, but at the same time, if PC 1 is removed, some constant part of the signal will be also removed. Additional information is needed to extract useful signal from PC 1. Second PC, obtained by grouping of signals corresponding to the second and third singular values (Fig. 4 and 3 right), contains the signal of annual hydrological cycle almost without stripes. Seasonal changes of the amount of water can be clearly seen in Amazon, Congo, Indus basins, together with changes from summer to winter in Russia, Canada, USA, China, etc. Third PC contains secular changes (Fig. 5 and 3 right). The gravity decrease in southeastern Greenland and West Antarctica to the east of the Ross Sea, caused by


FIGURE 5.

Third PC, representing secular changes.

melting of ice is clearly visible. The melting of mountain glaciers in Alaska, Himalayas and other regions is also evident. Higher order PCs contain high frequency components of noise, including those related to the stripes and transient events, such as the coseismic deformation after the Sumatra earthquake in December 2004, etc. Then MSSA was applied to 1-degree resolution GRACE monthly EWH maps, calculated from the expression (2). Principal components obtained by calculations in spatial domain appear to be, up to the order of singular values, similar to those, obtained in spectral domain. They also reflect permanent components of signal and stripes, signal with annual cycle, secular changes and high frequency components separated from each other. We put the animation of the obtained PCs on the web-site http://lnfm1.sai. msu.ru/~tempus/GRACE/index.htm.

Discussion
Preliminary results show, that MSSA helps to distinguish physically meaningful components, separate theme, reduce the stripes influence. Exact physical interpretation requires comparison with hydrological, GIA, ocean circulation models and other observations, such as ice mass balance, which is a topic for future study. We do not think that we need the rotation of PCs here to increase their meaningfullness, though it can be done in case of regional studies, as it has been recommended in [14]. Among the remaining questions include: (i) what is the useful part of the signal in PC 1, (ii) how to reduce boundary effects for the first and the last L points of PCs and (iii) how to better separate secular change from annual and, probably, other periodic signals? As for the last question ­ possible answer could be that better separation will be achievable when measurement time span increase, which will give more possibilities for the choice of L. Comparison of MSSA (L = 18) with PCA (L = 1) showed the following advantages of MSSA: if the first MSSA PC is almost invariable in amplitude, the first PCA component is changing periodically, though they both have spatially stable footprints. The second PC, which reflects the annual cycle, in PCA illustrates only amplitude variations of EWH in certain areas, such as Amazon and Congo basins, while MSSA-maps show


seasonal movements of these evolving patterns in the regions of their location. Thus, MSSA results reflect the dynamic of space-time correlation patterns better and allows better filtering. In addition, MSSA method has greater flexibility and could be useful in the analysis of other satellite observations, such as satellite altimetry and precipitation data, which are among the topics of future studies.

ACKNOWLEDGMENTS
This work is supported by the President of Russia grant MK-4234.2009.5. The GRACE data products are from the University of Texas at Austin Center for Space Research via NASA/JPL's PODAAC. The Ohio State University component of the research is supported by NASA and the Ohio State University Climate, Water and Carbon Program. Computations were performed with use of Ohio State Supercomputing Center cluster GLENN.

RE FE RE NCE S
1. Bettadpur S., Level-2 Gravity Field Product User Handbook, 2007, ftp://podaac.jpl.nasa. gov/pub/grace/doc/L2- UserHandbook_v2.3.pdf 2. Case K., Kruizinga G, Sien-Chong Wu, GRACE Level 1B Data Product User Handbook, 2004, ftp://podaac.jpl.nasa.gov/pub/grace/doc/Handbook_1B_v1.2.pdf 3. Ghil M., R.M. Allen, M.D. Dettinger et al., Advanced spectral methods for climatic time series, Rev. Geophys. 40(1), 3.1-3.41, 2002. 4. Golyandina N., Method "Caterpillar-SSA": analysis of the time series. (In russian) SPB., BBM, 2004 5. Golyandina N.,Nekrutkin V., Zhigljavskyet A., Analysis of time series structure: SSA and related techniques, Chapman & Hall/CRC, 2001 6. Hamlington B., R. Leben, R.S Nerem. Using Empirical Mode Decomposition to Calculate Patterns of Global Sea Level Change. OSTST meeting 2008. http://www.aviso.oceanobs.com/ fileadmin/documents/OSTST/2008/hamlington.pdf 7. Han Shin-Chan, Shum C.K., Jekeli Ch. et al. Non-isotropic filtering of GRACE temporal gravity for geophysical signal enhancement. Geophys. J. Int., Vol. 163, Issue 1 2005, pp. 18­25 8. Jollife I.T., Principal Component Analysis, Springer, 2001 9. Kenyon S. et al. Toward the next Earth gravitational model. 2007. http://earth- info.nga. mil/GandG/wgs84/gravitymod/new_egm/EGM08_papers/EGM- 2007- final.pdf 10. Klees R., Revtova E., Gunter B. et al., The design of an optimal filter for monthly GRACE gravity models. Geophys. J. Int., Vol. 175, Issue 5768, 2008, pp. 417­432, 11. Kusche J., R. Schmidt, S. Petrovic, et al., Decorrelated GRACE time-variable gravity solutions by GFZ and their validation using a hydrological model. J. of Geodesy, N 83, 2009, pp. 903-913. 12. Panteleev V.L., Theory of figure of the Earth, (lectures in russian). 2000. http://lnfm1.sai. msu.ru/grav/russian/lecture/tfe/index.html 13. Rangelova E., Wal W., Braun A. et al., Analysis of GRACE time-variable mass redistribution signals over North America by means of principal components analysis. J. of Geophys. Res., Vol. 112, 2007. 14. Rangelova E., Sideris M., Contributions of terrestrial and GRACE data to the study of the secular ё geoid changes in North America. J. of Geodynamics, 46, 2008, pp. 131БAS143. 15. Swenson S., and J. Wahr, Post-processing removal of correlated errors in GRACE data, Geophys. Res. Letters, 33, 2006. 16. Schrama E., Wouters B., Lavallee D. Signal and noise in Gravity Recovery and Climate Experiment (GRACE) observed surface mass variations. J. of Geophys. Res., Vol. 112, 2007. 17. Wouters B., Schrama E., Improved accuracy of GRACE gravity solution through empirical orthogonal function filtering of spherical harmonics. Geophys. Res. Letters, Vol. 34, 2007.