TMCnet News
Reduced-Rank Sigma-Point Kalman Filter and Its Application in ENSO Model [Journal of Atmospheric and Oceanic Technology](Journal of Atmospheric and Oceanic Technology Via Acquire Media NewsEdge) ABSTRACT The huge computational expense has been a main challenge while applying the sigma-point unscented Kalman filter (SPUKF) to a high-dimensional system. This study focuses on this issue and presents two methods to construct a reduced-rank sigma-point unscented Kalman filter (RRSPUKF). Both techniques employ the truncated singular value decomposition (TSVD) to factorize the covariance matrix and reduce its rank through truncation. The reduced-rank square root matrix is used to select the most important sigma points that can retain the main statistical features of the original sigma points. In the first technique, TSVD is applied on the covariance matrix constructed in the data space [RRSPUKF(D)], whereas in the second technique TSVD is applied on the covariance matrix constructed in the ensemble space [RRSPUKF(E)]. The two methods are applied to a realistic El Niño-Southern Oscillation (ENSO) prediction model [Lamont- Doherty Earth Observatory model, version 5 (LDEO5)] to assimilate the sea surface temperature (SST) anomalies. The results show that both the methods are more computationally efficient than the full-rank SPUKF, in spite of losing some estimation accuracy. When the truncation reaches a trade-off between cost expense and estimation accuracy, both methods are able to analyze the phase and intensity of all major ENSO events from 1971 to 2001 with comparable estimation accuracy. Furthermore, the RRSPUKF is compared against ensemble square root filter (EnSRF), showing that the overall analysis skill of RRSPUKF and EnSRF are comparable to each other, but the former is more robust than the latter. (ProQuest: ... denotes formulae omitted.) 1. Introduction The purpose of data assimilation is to optimally estimate the state of a physical system using imperfect numerical models and noisy observations. The Kalman filter (KF) and its derivatives have been widely applied in atmospheric and oceanic sciences in recent years (e.g., Nerger 2003; Robert et al. 2006; Evensen 2003; Deng and Tang 2009). The standard KF is a linear analysis based on the Bayes's theorem, whereas the extended Kalman filter (EKF) is an expansion of KF, in which the linearization is employed for nonlinear systems. The main drawback of EKF is the requirement of Jacobian or tangent linear model (TLM) for the linearization of the nonlinear models, which is extremely difficult to implement in a large-dimensional system (Ambadan and Tang 2009; Julier and Uhlmann 1997). In addition, the EKF only retains the first- and second-order statistical moments in evaluating error statistics, often leading to underestimation in error covariance (e.g., Evensen 1992; Miller et al. 1994; Gauthier et al. 1993). Instead, the ensemble Kalman filter (EnKF), proposed by Evensen (1994), estimates the error statistics using ensembles generated by integrating the nonlinear models (Evensen 1994; Houtekamer and Mitchell 1998). Unlike EKF, EnKF neither requires linearization of the dynamical model nor neglects the contribution from higher-order terms to error statistics (Evensen 1997). EnKF has become a powerful and widely used method within the data assimilation community because of its algorithmic simplicity, relative ease of implementation, and affordable computational cost (Evensen 2003; Li and Xiu 2008). An extensive review of the literature on EnKF can be found in Evensen (2003, 2009). To avoid the underestimation of forecast error covariance due to the decrease of ensemble spread, the observation is often perturbed in the standard EnKF, also known as stochastic EnKF (Houtekamer and Mitchell 1998; Burgers et al. 1998). However, Whitaker and Hamill (2002) found that the perturbed observations can act as an additional source of sampling error, making the estimate of analysis error covariance less accurate. They developed a deterministic EnKF known as the ensemble square root filter (EnSRF) to avoid the perturbation of observation, which uses the traditional Kalman gain to update the ensemble mean and the reformulated Kalman gain to update the deviation from the mean. For the same ensemble size, EnSRF can yield better results than the EnKF without any additional computational expense (e.g., Whitaker and Hamill 2002; Tippett et al. 2003). Amain challenge inEnKFandEnSRFis howto choose an appropriate ensemble size.Asmall ensemble size often results in inaccurate representation of the error covariance matrix, whereas a large ensemble size is not computationally feasible for high-dimensional systems (Lorenc 2003; Anderson 2012). The true prediction error covariance is usually unknown in EnKF and EnSRF. Therefore, some specific questions are likely to arise:How to generate, with a greater degree of accuracy, the finite ensembles for the optimal estimate of prediction error covariance? How accurate is the estimated prediction error covariance for a given ensemble size? The methods currently used for generating ensemble members in the standard EnKF seem unable to answer these questions (Tang et al. 2014b). Recently, Ambadan and Tang (2009), as well as Luo and Moroz (2009), have introduced a new ensemblebased filter to the field of earth science using the simple Lorenz system that is called the sigma-point Kalman filter (SPKF). SPKF is a derivativeless Kalman filter optimal for nonlinear state-space estimation (Julier et al. 1995; Julier and Uhlmann 2004; Julier 2002; Van der Merwe et al. 2004). In SPKF, the error statistics and model state are analyzed using the transformation of deterministically chosen minimal set of weighted sample points called sigma points, which captures the true mean and covariance of the prior random variable completely and the posterior mean and covariance to the second and third order for Gaussian inputs (Julier et al. 1995; Julier and Uhlmann 2002; Van der Merwe et al. 2004). The SPKFs family includes the sigma-point unscentedKalman filter (SPUKF) based on scaled unscented transformation (Julier et al. 1995; Julier 2002), the sigmapoint central differenceKalman filter (SPCDKF) based on Stirling's interpolation formulas (Nørgaard et al. 2000a,b; Ito and Xiong 2000), and their square root variants (Van der Merwe 2004; Van der Merwe and Wan 2001). This paper confines to only the SPUKF for simplicity. A remarking challenge in implementing SPKF for a realistic high-dimensional system is the large computational expense (Fan and You 2009). For a high-dimensional system like an ocean or atmosphere GCM, it is computationally not feasible to satisfy the requirement that the number of sigma points should be greater than twice the number of augmented state1 (Ambadan and Tang 2009; Luo and Moroz 2009; Han and Li 2008; Chandrasekhar et al. 2008). In the case of additive noise, the state augmentation can be avoided using iterated and trapezoidal approaches and the number of sigma points can be reduced to 2LM11 for an LM-dimensional system (Fan and You 2009; Van der Merwe and Wan 2001). There are some reduced-order sigma-point schemes using simplex sigma points (e.g., Julier and Uhlmann 2002; Julier 2003), which can decrease sigma points to LM11 for an LM-dimensional system. This number is still computationally infeasible for a realistic climate model. Recently, Ambadan and Tang (2009) suggested a subspace projection approach for the dimensional compression by the application of truncated principal component analysis (TPCA) on the multidimensional sigma-point space. In this technique the sigma points in the principal component space, which will retain the main features of the original sigma-point space, are used to approximate the error propagation. Luo and Moroz (2009) suggested that the truncated singular value decomposition (TSVD) can be applied on the error covariance matrix to reduce the sigma points. Another method to reduce the rank of SPUKF is throughCholesky decomposition (e.g., Stewart 1998; Chandrasekar et al. 2008). This study focuses on the investigation of the possibility of applying a reduced-rank SPUKF to a realistic climate model. The important properties and the estimation accuracy of the reduced-rank SPUKF are explored. Emphasis is placed on the implementation of the reduced-rank SPUKF and on its comparison against a full-rank SPUKF using a realistic climate model, which was not reported before. The model used is the Lamont- Doherty Earth Observatory, version 5 (LDEO5), version of the Zebiak-Cane El Niño prediction model, which is an intermediate complex model. The paper is organized as follows. Section 2 describes the SPUKF and its practical algorithm including the dimensional reduction methods. In section 3, there is a brief description for the Zebiak-CaneENSOmodel (ZCmodel). Section 4 describes the experiments and results. Section 5 contains the discussions and conclusions. 2. Methods a. EnKF, EnSRF, and SPUKF Suppose Xt is the state of a dynamical system of dimension L at time t and Yt is the observation at time t. The state-space equations of the dynamical system can be represented as below: ... (1) ... (2) where f ( ) is the nonlinear model that takes state Xt21 to Xt , and qt21 is the random model error following a Gaussian distribution with zero mean and covariance Qt21, h( ) is the observation operator, and rt is the observation noise, which is also Gaussian with mean zero and covariance Rt . 1) THE FRAMEWORK OF ENKF With the EnKF introduced by Evensen (Evensen 1994; Houtekamer and Mitchell 1998), ensemble members are integrated forward and updated whenever new observations are available, namely, ... (3) ... (4) ... (5) ... (6) where Xatis the analyzed state and H is the linearized observation operator. The subscript b indicates the model forecast state, Kt is the Kalman gain, Pbt is the forecast error covariance, E[.] represents the expectation value, the overbar denotes the ensemble mean,and Pat is the analysis error covariance; I is an identity matrix. 2) FORMULATION OF ENSRF In the classic EnKF schemes, randomly perturbed observations are usually used to avoid the underestimation of the analysis error covariance (Houtekamer and Mitchell 1998; Burgers et al. 1998). However, perturbing the observation will act as another source of sampling error and decrease the accuracy of the analysis error covariance matrix (Whitaker and Hamill 2002). Thus, an optimal approach that makes use of unperturbed observations but avoids systematical undermining of analysis error covariance is developed by Whitaker and Hamill (2002) and called EnSRF. For convenience of discussion, Eq. (3) can be rewritten as the sum of mean analysis state Xat and the deviation from the mean state Xa0t , namely, ... (7) ... (8) where Kt is the traditional Kalman gain and fKt is the gain used to update the deviation from the mean. For the traditional EnKF (Kt 5fKt) and Y0t is the randomly perturbed observation errors (Whitaker and Hamill 2002; Burgers et al. 1998). In EnSRF, Eq. (8) will be written as ... (9) Here fKt is chosen such that it will satisfy the analysis error covariance Eq. (6), namely, ... (10) In EnSRF, observations can be also processed one at a time, in that case HPb t HT and R reduced to scalars. Rewriting Eq. (9), if fKt 5aKt , where a is a constant, on solving for a, the solution obtained is (Whitaker and Hamill 2002) ... (11) 3) SPUKF In this section the algorithm of SPUKF is briefly introduced. In EnKF, the ensemble members are usually generated randomly using some well-designed random perturbation schemes (e.g., Evensen 2003). In EnKF and EnSRF, the ensemble size is usually determined considering the limitations dictated by complexity of the assimilation system, or calculation capacity, or in some cases by sensitivity studies (Mitchell et al. 2002; Lorenc 2003; Curn et al. 2012; Liang 2007), which is not deterministic in a rigorous statistical sense. In contrast, SPUKF uses a deterministically chosen scheme to generate ensemble members (i.e., sigma points) by the scaled unscented transformation formula (Julier 2002): ... (12) There are a set of weights associated with these sigma points ... The weights associated with the covariance calculation are ... or i 5 0, and ... (here the superscripts m and c correspond to mean and covariance, respectively). Equation (12) generates 2L11 ensemble members, where L is the dimension of augmented states concatenating the model states, process noise, and observation noise, and ... denotes the ith column of the weighted square root of the analysis error covariance matrix. Here ... is a scaling parameter, a varies between 0 and 1, 0 # a # 1, and determines the spread of the sigma points around themean state ... it should be ideally small tominimize the higher-order effects when the nonlinearities are strong. The term k is a control parameter that guarantees positive semidefiniteness of the covariance matrix and is a very important issue in implementing SPUKF. The value of k also depends on the value of a, and k makes sure that l is positive, otherwise the weight of the zeroth sigma point w(m) 0 ,0 and the nonpositive weights canmake the covariance nonpositive semidefinite (Van der Merwe 2004; Julier 2002). The term b is a nonnegative weighting term introduced to reduce the approximation error. If the state distribution is Gaussian, b52 is an optimal choice (Van der Merwe 2004; Julier 2002). In this paper the values of a51, k50, and b52 are used; their rationales and the detailed derivations can be found in Van der Merwe (2004) and references therein. SPUKF makes use of the following reformulated equation to find out the Kalman gain Kt and analysis error covariance matrix ... ... (13) ... (14) where Pxtyt is the cross covariance matrix of the state and observation prediction error, and Pytyt is the covariance matrix of the observation prediction error. For a detailed derivation see Van der Merwe (2004) and Ambadan and Tang (2009). When the new observation arrives, states are updated using the following equation and the analysis error covariance are updated using Eq. (14): ... (15) At this step one can generate new sigma points, given by ... (16) The sigma-point vector is propagated through Eqs. (1) and (2). The mean and covariance are calculated using the following formula: ... (17a) ... (17b) The cross covariance Pxtyt and projection covariance Pytyt are evaluated by the following equations: ... (18) Using the SPUKF, the second-order accuracy can be preserved in the mean and covariance. Higher-order information of the system in the covariance can be included by using the scaling parameters a and b. (Julier 2002; Julier and Uhlmann 2004). b. Reduced-rank sigma-point unscented Kalman filter (RRSPUKF) As discussed above, Eq. (12) produces 2L11 sigma points used for Eqs. (13) and (14) to estimate prediction error covariance and Kalman gain. This is computationally infeasible for a realistic model. To significantly decrease the number of sigma points, some methods have been proposed as mentioned in the introduction. 1) TSVD IN THE DATA SPACE, RRSPUKF(D) In this method, TSVD is applied on the covariance matrix constructed in the data space, so the filter will be denoted as RRSPUKF (Data) [RRSPUKF(D)]. As can be seen in Eq. (12), the numbers of sigma points are determined by the dimension of error covariance matrix. Thus, the central point is to effectively reduce the rank of the error covariance matrix, which will be achieved by the TSVD method. The analysis error covariance Pat in Eq. (14) is symmetric and can be decomposed as ... (19) where ... ) is a diagonal matrix consisting of the eigenvalues s2t ,i, sorted in descending order; that is ... is the matrix consisting of the corresponding eigenvectors ... To reduce the number of sigma points, a reduced-rank error covariance matrix, denoted by cPa t , can be approximately obtained by the truncation of Eq. (19), namely, ... (20) where l is the truncation number. The new sigma points are generated as follows: ... (21) where s is the square root of the eigenvalue and e is the eigenvector. The new sigma-point vector now becomes ... (22) The number of sigma points by Eq. (22) is 2l 11. 2) TSVD IN THE ENSEMBLE SPACE, RRSPUKF(E) In this technique, TSVD is applied on the covariance matrix constructed in the ensemble space (Ensemble) or in short [RRSPUKF(E)]. The RRSPUKF(D) can effectively save the computational time through reducing the sigma points from 2L 11 to 2l11. The saving is especially significant when the error covariance Pat can be approximated well by a few leading modes so that a small l can have a good truncation accuracy. However, a large challenge in RRSPUKF(D) is the computation of Pat itself. When the state dimension is large, it is computationally expensive, even unavailable in storage, to compute such a big matrix. Another approach to seek the most important sigma points, which can avoid the aforementioned challenge, is through TSVD in the ensemble space (Tang et al. 2014b). When the state dimension LM is much greater than the ensemble dimension n, that is, LM n, it is possible to compute the SVD on the (n 3 n) covariance matrix (Von Storch and Hannöschock 1984; Wilks 2011), which is computationally affordable for a system with large dimensions. In this method all the ensembles are updated according to Eq. (3), and then from these members we find the analysis error covariance matrix PaE t in the ensemble subspace according to ... (23) The SVD is performed on the (n 3 n) covariance matrix: ... (24) where Fat is the eigenvector and Gat is a diagonal matrix of the eigenvalues. The eigenvectors of PaE t and Pat are different, but the leading n eigenvectors of Pat can be computed from the eigenvector Fat of PaE t using their relationship (Von Storch and Hannoschöck 1984; Wilks 2011), namely, ... (25) where l51, . . . , n. The role of the denominator is tomake sure that the resulting Ea t has unit length. The maximum value of l is (n 2 1)/2, so that the total number of ensemble members is n for every time step. The new sigma points are generated in the same way as in Eq. (21). c. Covariance localization To avoid the possible impact of spurious correlation on the analysis, covariance localization is applied in this study (Hamill et al. 2001; Houtekamer and Mitchell 2001). The background error covariance matrix Pb t is (Schur product) multiplied by a distance-dependent localization correlation matrix r, and their product will be also a covariance matrix (Horn 1990). The modified covariance matrix will be roPbt : ... (26) And the members of the correlation matrix ... (27) where r is a localization function that has a value of 1 at the analysis grid and its value decreases as the distance from the analysis grid increases, d is the decorrelation length, and r is the distance between the grid i and grid j. The localization radius is fixed to be 458. 3. The Zebiak-Cane model The model used is the LDEO5 of the Zebiak-Cane model (Chen et al. 2004), which is a classic El Niño prediction model (Zebiak and Cane 1987). The model computes anomalies of atmospheric and oceanic fields, relative to a specified monthly mean climatology. The atmosphere dynamics follows Gill (1980), using a steady-state linear shallow water equation on an equatorial beta plane. The circulation is forced by a heating anomaly that depends on the sea surface temperature (SST) anomaly and moisture convergence. The ocean dynamics uses the reduced gravity model (Cane 1984). The surface currents are generated by spinning up the model with monthly mean climatological wind. The thermodynamics describes the SST anomaly and heat flux change. The model time step is 10 days. The model domain is confined in the tropical Pacific Ocean (298S-298N, 101.258E-73.1258W). The grid for ocean dynamics is 28 longitude 3 0.58 latitude and 5.6258 longitude 3 28 latitude for the atmosphere model and SST physics, respectively. The reconstructed sea surface temperature anomaly (SSTA) data by Kaplan et al. (1998) are used for the true state and are also located on the model grids. This dataset has been used by Chen et al. (2004) for the successful long-term retrospective prediction of ENSO over 100 years (1856- 2003). 4. Experiments and results a. ENSO state estimation with RRSPUKF(D) In the state estimation experiments, Kaplan SSTA reanalysis dataset (27.58N-27.58S, 122.58E-92.58W) at the resolution of 58 3 58 in both longitude and latitude (Kaplan et al. 1998) is assimilated. Assimilation experiments are performed every month for the period from January 1971 to December 2000. All the experiments assume the observation errors and model errors to be uncorrelated in space and time. To solve the underestimation problem of the background error covariance, a multiplicative inflation scheme is considered as in other assimilation systems (e.g., Hamill and Whitaker 2005). Namely, ... (28) where ... is raw background error covariance and v is the inflation coefficient usually obtained by tuning experiments (also called tunable parameter). In this study, the value v is set to 0.15, for the ensemble size of 41, which is found to result in the best analysis. This direct inflation of background error covariance is equivalent to the inflation of model states if the model errors are assumed to be proportional to the background error (Li et al. 2009) but can avoid the divergence and inbreeding problem of the filter of finite ensemble size (Furrer and Bengtsson 2007; Karspeck and Anderson 2007; Hamill and Whitaker 2005). The total number of spatial grid points within the assimilation domain is 540. When SPUKF is implemented, the state vector is refined as the concatenation of the model states, model errors and measurement errors, thus the augmented state dimension will become 3 3 540. Consequently, the number of a full SPUKF sigma points is 23(33540)11. To decrease the sigma points, the RRSPUKF(D) is applied to the analysis error covariance matrix Pa t . The total number of sigma points 2L11 is now reduced to 2l11. The estimate accuracy of truncated covariance matrix depends on the choice of l. The selection of l should be chosen cautiously in such way that it is not too small or too large. If l is too small, some of the information from the Pa t will be lost, and a large l will create computational burden. b. Sensitivity experiments In this section the performance of RRSPUKF(D), as a function of the number of truncated modes of the SVD, is examined. A comparison among different truncated modes is shown in terms of root-mean-square error (RMSE) between the analysis values and the observed values as shown in Fig. 1 (solid line). Different numbers of sigma points ranging from 5 to 101 are tested. As seen in Fig. 1 (solid line), when the number of sigma points increases, the RMSE of Niño-3.4 SSTA analysis against the observation decreases. However, when the number of sigma points used are beyond 40 sigma points (i.e., l 5 20), there is little improvement in the analysis. One of the reasons for higher RMSE when the ensemble number is small is the underestimation of the covariance matrix by the finite truncated modes. The RMSE and correlation skill for the period 1971- 2000 over the entire basin are shown in Figs. 2 and 3. Like Fig. 1, the analysis has a large error when ensemble size is small, as indicated by the higher RMSE and lower correlation coefficient. This is a case of filter divergence, which usually occurs when the variance of the prior distribution is small compared to the measurement uncertainties and the filter assigns small weight to the measurement (Karspeck and Anderson 2007; Furrer and Bengtsson 2007). When the ensemble number increases, RMSE decreases; but beyond 40 ensembles, there is little improvement in the analysis, as evidenced in these figures. This suggests that an ensemble size of 40 is sufficient to approximate the mean and spread of the prior distribution as well as to characterize the coherent structure between the measurement and the prior state. In Figs. 2 and 3, we can also observe that different ensemble sizes from 31 and 41 to 51 can lead to similar behavior in both RMSE and correlation, but below 21 ensembles the skill is very poor, indicating the sensitivity of assimilation performance to the truncated modes is not high if the basic characteristics can be captured by these modes. c. Comparison of RRSPUKF(D) and SPUKF Based on sensitivity experiments of the analysis error to the truncation modes l, it is found that 41 sigma points (i.e., l 5 20) are enough to fairly estimate the mean and covariance. The resulting analysis is very good and comparable with a full SPUKF as shown in Figs. 4 and 5. The variance explained by the 41 sigma points is more than 85% of the total variance (also see Fig. 10). Figures 4b and 4a show the RMSE of the Niño-3.4 index (the averageSSTAin 58N-58S, 1208-1708W)between the analysis and the observation for both the full-rank SPUKF and RRSPUKF(D) with 41 sigma points, respectively. As seen in Fig. 4, the RMSE of RRSPUKF(D) is similar to that of the full-rank SPUKF in most of the assimilation window. This can also be seen in Fig. 5 where the RMSE and correlation between the analysis and observation for both the schemes are compared in the entire basin. The bootstrap experiments are conducted to perform the statistical significance test for RMSE and correlation skill, as follows. 1) For a given grid point, the analysis and observed sample of the same time are combined to form a pair sample during the period of 1971-2001; 2) 95% of the pair samples are randomly chosen to calculate the correlation coefficient and RMSE between analysis and observation; 3) process 2 is repeated 10 000 times to produce 10 000 correlations and RMSEs, from which the standard deviation is obtained; and 4) the 95% confidence interval is used as the threshold values of sampling errors. Figures 5e and 5f show the difference of RMSE and correlation between SPUKF and RRSPUKF(D), with the statistical significance area shaded at the confidence level of 95%. As seen in Fig. 5f, the correlation skill of RRSPUKF(D) has no significant difference from that of a full-rank SPUKF for most domains of the assimilation. However, the difference of RMSE between RRSPUKF(D) and SPUKF is significant for almost the entire assimilation domain as shown in Fig. 5e. This is mainly because the RMSE itself is very small in both schemes, leading to very small criteria values. Thus, it can be argued that the RRSPUKF could fairly well simulate the performance of the full-rank SPUKF in terms of the measure of correlation skill, and it can also generate a small RMSE. Figure 6 shows the one-step forecast (one-month lead forecast) of Niño-3.4 SSTA index [i.e., Xb t by Eq.(1)] from the RRSPUKF(D) with 41 sigma points during the time period 1971-2000. The solid line is the model state before assimilation (i.e., the prediction initialized from the analysis of last step analysis) and the plus sign (1) denotes the observations. As seen in Fig. 6, the RRSPUKF(D) filter has a good capability for estimating the phase and intensity of all major ENSO events during the entire period. The large error occurs at few time steps when there is high-frequency variability in SSTA, which is reasonable given that the observation used for assimilation can also contain uncertainties. Another possible explanation is due to the decadal variability of ENSO, which may be not well captured by the ZC model (Karspeck and Anderson 2007). d. Forecast skill To further explore the impact of different assimilation schemes on ENSO prediction, the SSTA hindcasts of the duration of 12 months, initialized from the first day of each calendar month, were performed during the period 1971-2000. The hindcast skill is evaluated in terms of the RMSE and correlation skill of predicted SSTA against the observed counterpart. Figure 7 shows the RMSE and correlation coefficient skill of Niño-3.4 SSTA, as a function of lead time for different ensemble sizes. As seen in the figure, both the RMSE and correlation skills decrease with lead time, as characterized in a nonlinear or stochastic dynamical system like ENSO (Jin et al. 1994; Chen et al. 2004; Penland and Sardeshmukh 1995; Kirtman and Schopf 1998). On the other hand, the RMSE decreases and correlation skill increases as the ensemble size increases. The prediction skills start to saturate and have little improvement when the number of ensemble size is approximately above 20. This result is consistent with those obtained in the above sensitivity experiments. e. Comparison between RRSPUKF(D) and RRSPUKF(E) approach An alternate method for finding the dominant sigma points is through TSVD in the ensemble subspace as introduced in section 2b. This method is computationally feasible even for a high-dimensional model. Figure 1 compares the RMSE skill of analyzed Niño-3.4 SSTA index obtained by RRSPUKF(D) (solid line) and RRSPUKF(E) (dotted line) methods for the period from 1971 to 2000 against the observation. As seen in Fig. 1, the RRSPUKF(D) is better than RRSPUKF(E) when the ensemble size is small, indicating that the latter is more inefficient in using small ensembles to capture the main features of the full covariance matrix. This is because the RRSPUKF(D) directly discomposes the full covariance, whereas the RRSPUKF(E) uses the ensemble to approximate the covariance matrix. However, as the ensemble size increases, the difference between the two schemes becomes very small. Figure 8 shows the RMSE (Figs. 8a,c) and the correlation (Figs. 8b,d) skill of analyzed SSTA using the two schemes for the period 1971-2000 when the ensemble size is 41. Figures 8e and 8f show the differences in RMSE and correlation between the two schemes. Similar to Fig. 5 as discussed in section 4c, the bootstrap method is used to perform the statistical significance test, as shaded in the figure. The confidence level is 95%. Apparently, there is no significant difference of correlation and RMSE skill between the two schemes for almost the entire domain of assimilation, indicating that the skill of RRSPUKE(E) is completely comparable to the RRSPUKF(D). Figure 9 shows the explained variance as a function of ensemble size for both RRSPUKF(E) and the RRSPUKF(D), where the variance explained by truncated modes is calculated as below, with each variable having the same denotation as that in section 2b: ... (29) ... (30) As seen in Fig. 9, the explained variance is larger in RRSPUKF(E) than in the RRSPUKF(D). This seems interesting since the RRSPUKF(D) usually has better performance than the RRSPUKF(E) for the same ensemble size as shown in Fig. 1. One possible reason is that both schemes have different total variances (denominators), making them incomparable. It should be noted that the RRSPUKF(E) only uses half of the ensemble (i.e., throwing out the other half) to calculate the covariance matrix to keep the same ensemble size of 2l11, given a prescribed number of l. Furthermore, since the ensemble size is not large enough, there is an inaccurate representation or underestimation for the prediction covariance. The SVD that is performed on such an inaccurate covariance matrix can lead to suboptimal sigma points, suggesting that explained variance cannot be taken as the single entity to decide the number of sigma points to be used. f. Comparison between the RRSPUKF(D) and EnSRF Assimilation experiments are also performed using EnSRF algorithm with the same experimental setup as those in RRSPUKF(D). In all the experiments both the RRSPUKF(D) and EnSRF start from the same initial conditions. Shown in Fig. 10a is the RMSE of Niño-3.4 index analysis by RRSPUKF(D) and EnSRF (Fig. 10b), both with ensemble size of 41 and against the observation counterpart. The RMSE and correlation over the entire basin is also compared in Fig. 11, showing that the RRSPUKF(D) results in the smaller RMSE and slightly higher correlation skill than EnSRF, although their differences are not statistically significant for most of domain as shown in Figs. 11e and 11f. The better performance of RRSPUKF over the EnSRF is probably because of its superior mathematical properties explained in sections 1 and 2. Namely, the RRSPUKF can better estimate the first- and secondorder statistical moments using SVD technique to approximate the known true statistical moments in contrast to an estimate of unknown true moments in EnSRF. In addition, in RRSPUKF(D) the truncation error can be quantified but in the case of EnSRF the selection of number of ensemble is arbitrary as explained in section 2. Like Fig. 7, Fig. 12 shows the hindcast experiments of Niño-3.4 SSTA initialized by the RRSPUKF(D), EnSRF, and RRSPUKF(E). The hindcast settings are the same as those in Fig. 7. Bootstrap experiments are performed to test the statistical significance for both RMSE and correlation skill. The experiment setting was similar to that described in the section 4c. The error bars indicate the sampling errors at 95% confidence interval. As shown in Fig. 12, all the three schemes have similar skills at lead times shorter than 4-5 months. As the lead time increases the RMSE of RRSPUKF(D) and RRSPUKF(E) schemes are significantly lower than that of the EnSRF scheme but correlation skills are still comparable for all three schemes, indicating the RMSE skill is more sensitive to the assimilation scheme. 5. Discussion and conclusions The ensemble-based Kalman filters like the ensemble square root filter (EnSRF) have successfully demonstrated their ability in oceanic and atmospheric data assimilation. The main disadvantages of these methods are the uncertainty in constructing the ensembles and the requirement of a greater number of ensembles for the accurate approximation of the mean and variance of the prior and posterior random variables (Ambadan and Tang 2009; Kim et al. 2007; Chandrasekhar et al. 2008). The sigma-point unscented Kalman filter (SPUKF), which is relatively new in atmospheric and oceanographic data assimilation, deterministically chooses ensemble members and thereby reduces uncertainty in generating ensembles (Julier 2003). Also, the SPUKF makes use of a derivative-less optimal estimation technique that uses a minimal number of statistically weighted samples called sigma points to calculate the error statistics. When transformed through the nonlinear model they capture the true mean and covariance accurately to the second order in the Taylor series expansion (Julier 2002). The SPUKF uses reformulated equations to find the Kalman gain and error covariance so that there is no need to calculate the Jacobian or tangent linear model (TLM) as in the case of extended Kalman filter (EKF), and there is no need to linearize the nonlinear measurement function as in the ensemble Kalman filter (EnKF) family of filters (Tang and Ambadan 2009; Tang et al. 2014a). Despite all the merits of SPUKF, the need for twice the number of sigma points more than the system states is a major computational concern for the application of SPUKF onto realistic ocean or atmosphere models. Given the computational constraints of full-order error covariance estimation, the data assimilation community is very interested in the development of a reduced-order scheme for the state estimation. In this study, the reduced-rank sigma-point unscented Kalman filter (Data) [RRSPUKF(D)] and reduced-rank sigma-point unscented Kalman filter (Ensemble) [RRSPUKF(E)] are introduced and experimentally implemented to a realistic climatemodel. The challenge of the computational burden for SPUKF is successfully overcome by employing truncated SVD technique for dimensional compression or rank reduction of the error covariance matrix. The sensitivity experiments are carried out using different numbers of truncated modes after the application of SVD on the error covariance matrix. The results showed that the number of sigma points, even as low as 41, are enough to get an analysis very close to the fullorder SPUKF analysis. The analysis obtained using 41 sigma points is successfully able to estimate the phase and intensity of all majorENSOevents during the study period of 1971-2000. The estimation skill of RRSPUKF(D&E) with EnSRF in terms of RMSE and correlation coefficients is compared. When the same number of ensembles are used, the performance of RRSPUKF is slightly better than the EnSRF, which suggests that the RRSPUKF scheme can obtain more accurate estimates because of its better algorithmic properties. The advantage of RRSPUKF(E) over RRSPUKF(D) is in its affordable computational cost and lower memory usage. The main challenge of RRSPUKF(D) is to represent the full analysis covariance matrix explicitly, which is still computationally not affordable for highdimensional models. In the case of RRSPUKF(E) the ensembles are used to approximate the covariancematrix and it is calculated in the reduced dimension (ensemble space).Our results showed thatwhen the ensemble size is as small as 41, the performance of both the schemes RRSPUKF(D) and RRSPUKF(E) are comparable. This study explores the application of SPUKF and RRSPUKF(D&E) for a climate model of intermediate complexity. The results showed that RRSPUKF with 41 ensembles can have estimation accuracy close to that of a full-rank SPUKF, making it a potential candidate for data assimilation in large-dimensional atmospheric or ocean GCM. We are currently working on the application of RRSPUKF for data assimilation in an atmospheric GCM. Acknowledgments. This work is supported by the NSERC Discovery program and National Science Foundation of China (41276029). Yanjie Cheng is supported by National Science foundation of China (41175074). 1 The augmented state used by the SPKF concatenate the model states, process noise, and observation noise. Namely, if numerical model is an LM-dimensional system, the number of augmented state is 3LM. REFERENCES Ambadan, J. T., and Y. Tang, 2009: Sigma-point Kalman filter data assimilation methods for strongly nonlinear systems. J. Atmos. Sci., 66, 261-285, doi:10.1175/2008JAS2681.1. Anderson, J. L., 2012: Localization and sampling error correction in ensemble Kalman filter data assimilation. Mon. Wea. Rev., 140, 2359-2371, doi:10.1175/MWR-D-11-00013.1. Burgers, G., P. J. Van Leeuwen, and G. Evensen, 1998: Analysis scheme in the ensemble. Kalman filter. Mon. Wea. Rev., 126, 1719-1724, doi:10.1175/1520-0493(1998)126,1719: ASITEK.2.0.CO;2. Cane, M. A., 1984: Modeling sea level during El Niño. J. Phys. Oceanogr., 14, 1864-1874, doi:10.1175/1520-0485(1984)014,1864: MSLDEN.2.0.CO;2. Chandrasekar, J., I. S. Kim, D. S. Benstein, and A. J. Ridley, 2008: Reduced-rank unscented Kalman filtering using Cholesky-based decomposition. Int. J. Control, 81, 1779- 1792, doi:10.1080/00207170801891724. Chen, D., M. A. Cane, A. Kaplan, S. E. Zebiak, and D. Huang, 2004: Predictability of El Niño over the past 148 years. Nature, 428, 733-736, doi:10.1038/nature02439. Curn, J., D. Marinescu, G. Lacey, and V. Cahill, 2012: Estimation with non-white Gaussian observation noise using a generalised ensemble Kalman filter. Proc. IEEE Int. Symp. on Robotic and Sensors Environments (ROSE), Magdeburg, Germany, IEEE, 85-90, doi:10.1109/ROSE.2012.6402618. Deng, Z., and Y. Tang, 2009: Reconstructing the past wind stresses over the tropical Pacific Ocean from 1875 to 1947. J. Appl. Meteor.Climatol., 48, 1181-1198, doi:10.1175/2008JAMC2049.1. Evensen, G., 1992: Using the extended Kalman filter with a multilayer quasi-geostrophic ocean model. J. Geophys. Res., 97, 17 905-17 924, doi:10.1029/92JC01972. -, 1994: Sequential data assimilation with a nonlinear quasigeostrophic model using Monte Carlo methods to forecast error statistics. J. Geophys. Res., 99, 10 143-10 162, doi:10.1029/94JC00572. -, 1997: Advanced data assimilation for strongly nonlinear dynamics. Mon. Wea. Rev., 125, 1342-1354, doi:10.1175/ 1520-0493(1997)125,1342:ADAFSN.2.0.CO;2. -, 2003: The ensemble Kalman filter: Theoretical formulation and practical implementation. Ocean Dyn., 53, 343-367, doi:10.1007/s10236-003-0036-9. -, 2009: Data Assimilation: The Ensemble Kalman Filter. 2nd ed. Springer, 307 pp. Fan, C., and Z. You, 2009: Highly efficient sigma points filter for spacecraftattitude and rate estimation. Math. Probl. Eng., 2009, 507370, doi:10.1155/2009/507370. Furrer, R., and T. Bengtsson, 2007: Estimation of high dimensional prior and posterior covariance matrices in Kalman filter variants. J. Multivar. Anal., 98, 227-255, doi:10.1016/ j.jmva.2006.08.003. Gauthier, P., P. Courtier, and P. Moll, 1993: Assimilation of simulated wind lidar data with a Kalman filter. Mon. Wea. Rev., 121, 1803-1820, doi:10.1175/1520-0493(1993)121,1803: AOSWLD.2.0.CO;2. Gill, A. E., 1980: Some simple solutions for heat-induced tropical circulation. Quart. J. Roy. Meteor. Soc., 106, 447-462, doi:10.1002/qj.49710644905. Hamill, T. M., and J. S. Whitaker, 2005: Accounting for the error due to unresolved scales in ensemble data assimilation: A comparison of different approaches. Mon. Wea. Rev., 133, 3132-3147, doi:10.1175/MWR3020.1. -, -, and C. Snyder, 2001: Distance-dependent filtering of background error covariance estimates in an ensemble Kalman filter. Mon. Wea. Rev., 129, 2776-2790, doi:10.1175/ 1520-0493(2001)129,2776:DDFOBE.2.0.CO;2. Han, X., and X. Li, 2008: An evaluation of the nonlinear/ non-Gaussian filters for the sequential data assimilation. Remote Sens. Environ., 112, 1434-1449, doi:10.1016/ j.rse.2007.07.008. Horn, R., 1990: The Hadamard product. Matrix Theory and Applications, C. R. Johnson, Ed., Proceedings of Symposia in Applied Mathematics, Vol. 40, American Mathematical Society, 87-169. Houtekamer, P. L., and H. L. Mitchell, 1998: Data assimilation using an ensemble Kalman filter technique. Mon. Wea. Rev., 126, 796-811, doi:10.1175/1520-0493(1998)126,0796: DAUAEK.2.0.CO;2. -, and -, 2001: A sequential ensemble Kalman filter for data assimilation. Mon. Wea. Rev., 129, 123-137, doi:10.1175/ 1520-0493(2001)129,0123:ASEKFF.2.0.CO;2. Ito, K., and K. Xiong, 2000: Gaussian filters for nonlinear filtering problems. IEEETrans. Autom. Control, 45, 910-927, doi:10.1109/ 9.855552. Jin, F. F., J. D. Neelin, and M. Ghil, 1994: El Niño on the devil's staircase-Annual subharmonic steps to chaos. Science, 264, 70-72, doi:10.1126/science.264.5155.70. Julier, S. J., 2002: The scaled unscented transformation. Proceedings of the 2002 American Control Conference, Vol. 6, IEEE, 4555-4559, doi:10.1109/ACC.2002.1025369. -, 2003: The spherical simplex unscented transformation. Proceedings of the 2003 American Control Conference, Vol. 3, IEEE, 2430-2434, doi:10.1109/ACC.2003.1243439. -, and J. K. Uhlmann, 1997: A new extension of the Kalman filter to nonlinear systems. Signal Processing, Sensor Fusion, and Target Recognition VI, I. Kadar, Ed., International Society for Optical Engineering (SPIE Proceedings 3068), 182- 193, doi:10.1117/12.280797. -, and -, 2002: Reduced sigma point filters for the propagation of means and covariances through nonlinear transformation. Proceedings of the 2002 American Control Conference, Vol. 2, IEEE, 887-892, doi:10.1109/ ACC.2002.1023128. -, and -, 2004: Unscented filtering and nonlinear estimation. Proc. IEEE, 92, 401-422, doi:10.1109/JPROC.2003.823141. -, -, and F. Durrant-Whyte, 1995: A new approach for filtering nonlinear systems. Proceedings of the 1995 American Control Conference, Vol. 3, IEEE, 1628-1632, doi:10.1109/ ACC.1995.529783. Kaplan,A.,M.A.Cane,Y.Kushnir,A. C. Clement,M.B. Blumenthal, and B. Rajagopalan, 1998: Analysis of global sea surface temperature 1856-1991. J. Geophys. Res., 103, 18 567-18 589, doi:10.1029/97JC01736. Karspeck, A. R., and J. L. Anderson, 2007: Experimental implementation of an ensemble adjustment filter for an intermediate ENSO model. J. Climate, 20, 4638-4658, doi:10.1175/ JCLI4245.1. Kim, I. S., J. Chandrasekar, H. J. Palanthandalam-Madapusi, A. Ridley, and D. S. Bernstein, 2007: State estimation for large-scale systems based on reduced-order error-covariance propagation. Proc. American Control Conf., New York, NY, IEEE, 5700-5705, doi:10.1109/ACC.2007.4282477. Kirtman, B. P., and P. S. Schopf, 1998: Decadal variability in ENSO predictability and prediction. J. Climate, 11, 2804-2822, doi:10.1175/1520-0442(1998)011,2804:DVIEPA.2.0.CO;2. Li, H., E. Kalnay, T. Miyoshi, and C. M. Danforth, 2009: Accounting for model errors in ensemble data assimilation. Mon. Wea. Rev., 137, 3407-3419, doi:10.1175/2009MWR2766.1. Li, J., and D. Xiu, 2008: On numerical properties of the ensemble Kalman filter for data assimilation. Comput. Methods Appl. Mech. Eng., 197, 3574-3583, doi:10.1016/j.cma.2008.03.022. Liang, B., 2007: An ensemble Kalman filter module for automatic history matching. Ph.D. thesis, The University of Texas at Austin, 335 pp. Lorenc, A. C., 2003: The potential of the ensemble Kalman filter for NWP-Acomparison with 4D-Var. Quart. J. Roy. Meteor. Soc., 129, 3183-3203, doi:10.1256/qj.02.132. Luo, X., and I. M. Moroz, 2009: Ensemble Kalman filter with the unscented transform. Physica D, 238, 549-562, doi:10.1016/ j.physd.2008.12.003. Miller, R., M. Ghil, and F. Gauthiez, 1994: Advanced data assimilation in strongly nonlinear dynamical systems. J. Atmos. Sci., 51, 1037-1056, doi:10.1175/1520-0469(1994)051,1037: ADAISN.2.0.CO;2. Mitchell, H. L., P. L. Houtekamer, and G. Pellerin, 2002: Ensemble size, balance, and model-error representation in an ensemble Kalman filter. Mon. Wea. Rev., 130, 2791-2808, doi:10.1175/ 1520-0493(2002)130,2791:ESBAME.2.0.CO;2. Nerger, L., 2003: Parallel filter algorithms for data assimilation in oceanography. Ph.D. thesis, University of Bremen, Polar and Marine Research Rep. 487, 174 pp. Nørgaard, M., N. K. Poulsen, and O. Ravn, 2000a: Advances in derivative-free state estimation for nonlinear systems. Dept. of Mathematical Modeling, Technical University of Denmark Tech. Rep. IMM-REP-1998-15, 33 pp. -,-, and -, 2000b: New developments in state estimation of nonlinear systems. Automatica, 36, 1627-1638, doi:10.1016/ S0005-1098(00)00089-3. Penland, C., and P. D. Sardeshmukh, 1995: The optimal growth of tropical sea surface temperature anomalies. J. Climate, 8, 1999-2024, doi:10.1175/1520-0442(1995)008,1999: TOGOTS.2.0.CO;2. Robert, C., E. Blayo, and J. Verron, 2006: Comparison of reducedorder, sequential and variational data assimilation methods in the tropical PacificOcean. OceanDyn., 56, 624-633, doi:10.1007/ s10236-006-0079-9. Stewart, G. W., 1998: Basic Decompositions. Matrix Algorithms, Vol. 1, SIAM, 184 pp. Tang, Y., and J. Ambadan, 2009: Reply. J. Atmos. Sci., 66, 3501- 3503, doi:10.1175/2009JAS3273.1. -, -, and D. Chen, 2014a: Nonlinear measurement function in the ensemble Kalman filter. Adv. Atmos. Sci., 31, 551-558, doi:10.1007/s00376-013-3117-9. -, Z. Deng, K. K. Manoj, and D. Chen, 2014b: A practical scheme of the sigma-point Kalman filter for high-dimensional systems. J. Adv. Model. Earth Syst., 6, 21-37, doi:10.1002/ 2013MS000255. Tippett,M. K., J. L. Anderson, C. H. Bishop, T. M. Hamill, and J. S. Whitaker, 2003: Ensemble square root filters. Mon.Wea. Rev., 131, 1485-1490, doi:10.1175/1520-0493(2003)131,1485: ESRF.2.0.CO;2. Van der Merwe, R., 2004: Sigma-point Kalman filters for probabilistic inference in dynamic state-space models. Ph.D. thesis, Oregon Health and Science University, 377 pp. -, andE.A.Wan, 2001: The square-root unscentedKalman filter for state and parameter estimation. Proceedings of the IEEE International Conference on Acoustics, Speech, and Signal Processing, 2001 (ICASSP '01), Vol. 6, IEEE, 3461-3464. -, -, and S. I. Julier, 2004: Sigma-point Kalman filters for nonlinear estimation and sensor fusion: Applications to integrated navigation. AIAA Guidance, Navigation and Control Conf., Providence, RI, American Institute of Aeronautics and Astronautics, 5120-5122. Von Storch, H., and G. Hannoschöck, 1984: Comments on ''Empirical orthogonal function- analysis of wind vectors over the tropical Pacific region.'' Bull. Amer. Meteor. Soc., 65, 162. Whitaker, J. S., and T. Hamill, 2002: Ensemble data assimilation without perturbed observations. Mon. Wea. Rev., 130, 1913-1924, doi:10.1175/1520-0493(2002)130,1913: EDAWPO.2.0.CO;2. Wilks, D. S., 2011: Statistical Methods in the Atmospheric Sciences. Academic Press, 676 pp. Zebiak, S. E., and M. A. Cane, 1987: A model El Niño-Southern Oscillation. Mon. Wea. Rev., 115, 2262-2278, doi:10.1175/ 1520-0493(1987)115,2262:AMENO.2.0.CO;2. K. K. MANOJ Environmental Science and Engineering, University of Northern British Columbia, Prince George, British Columbia, Canada YOUMIN TANG Environmental Science and Engineering, University of Northern British Columbia, Prince George, British Columbia, Canada, and State Key Laboratory of Satellite Ocean Environment Dynamics, Hangzhou, China ZIWANG DENG Environmental Science and Engineering, University of Northern British Columbia, Prince George, British Columbia, and Department of Statistics, York University, Toronto, Ontario, Canada DAKE CHEN State Key Laboratory of Satellite Ocean Environment Dynamics, Hangzhou, China YANJIE CHENG National Climate Center, China Meteorological Administration, Beijing, China (Manuscript received 27 August 2013, in final form 14 March 2014) Corresponding author address: Dr. Youmin Tang, Environmental Science and Engineering, University of Northern British Columbia, 3333 University Way, Prince George BC V2N 4Z9, Canada. E-mail: [email protected] (c) 2014 American Meteorological Society |