Abstract
The four-dimensional variational data assimilation (4D-Var) method has been widely employed as an operational scheme in mainstream numerical weather prediction (NWP) centers. In addition to the ensemble data assimilation method, the randomization technique is still used to diagnose the standard deviations of background error in variational data assimilation (VAR) systems; however, such randomization techniques induce sampling noise, which may contaminate the quality of the standard deviations. First, this paper studies the properties of the sampling noise induced by the randomization technique. The results show that the sampling noise is on a small scale displaying high-frequency oscillations around the estimate compared with the estimate and this difference motivates the use of filtering techniques to eliminate the sampling noise effects. The characteristics of the standard deviation field of the control variables are also investigated, and the standard deviation fields of different model parameters have different scales and vary with the vertical model levels. To eliminate such sampling noise, the spectral filtering method used widely in the operational system and a modified spatial averaging approach are investigated. Although both methods have splendid performance in eliminating sampling noise, the spatial averaging approach is more efficient and easier to implement in operational systems. In addition, the optimal filtered results from the spatial averaging approach are dependent on model parameters and vertical levels, which is consistent with the variation in the standard deviation field. Finally, the spatial averaging approach is tested on the operational system at the global scale based on the YH4DVAR and the global NWP system, and the results indicate that the spatial averaging approach has positive effects on both analysis and forecast quality.
1. Introduction
A data assimilation system optimizes the combination of observations and a short-term forecast of the state variables to provide the best estimate of the initial state of the atmosphere [1]. Reports from the famous numerical weather prediction (NWP) Centre, European Centre for Medium-Range Weather Forecasts (ECMWF), show that the assimilated observations contribute 15% of the influence in any one analysis, and the complementary 85% is from the background which contains information from earlier assimilated observations [2]. The information from observations and a short-term forecast are combined and weighted by their corresponding error covariance matrix, the observation error covariance matrix (R matrix), and the background error covariance matrix (B matrix). The four-dimensional variational data assimilation (4D-Var) method is the mainstream scheme used by many leading NWP centers to estimate the initial state of the atmosphere for the NWP model [3]. 4D-Var can assimilate observations of nontraditional types, in particular from satellites or radar, and these observations are generally considered to be a key factor in the continuous improvement of the quality of NWP [4–6]. In the 4D-Var system, the B matrix has a profound impact on the analysis because it determines the weight of the prior state, propagates the information from observations, and imposes a balanced relationship between the model control variables [7, 8].
The specification of the B matrix is a key component of the 4D-Var system, but it remains a great challenge. The dimensions of the model state vector are of the order of in a current NWP system; thus, it is far too large for the B matrix to be explicitly represented. The lack of knowledge about the statistical properties of background error is also a limitation. To overcome these difficulties, significant efforts have been made to approximate the B matrix. In the operational variational data assimilation (VAR) system, generally, the B matrix has to be modeled by using control variable transforms (CVT) and represented as a product of operator form, , where C is the background error correlation operator and is the diagonal matrix of background error standard deviations [9–14]. Considering the constraints between multiple variables, the standard deviations need to be diagnosed in the VAR system. Fisher and Courtier [15] suggested that the standard deviations of background error can be estimated by a randomization technique if the B matrix is in the form . The randomization approximation of the standard deviations is widely implemented in the current operational centers [16].
The YH4DVAR system is an operational deterministic 4D-Var system that uses the global spectral model as a constraint to impose a dynamic balance on the assimilation [17]. Considering that the sampling noise introduced by the randomization technique has a negative effect on the estimate, a spectral filtering method truncated at a fixed wavenumber in an operational setting solves the issue of sampling noise in the YH4DVAR system [18]. The spectral filtering method enables the removal of most of the sampling noise while preserving the standard deviations of interest. However, the local variation characteristics of the standard deviations may be destroyed through fixed wavenumber truncation. In this paper, we will introduce the modified spatial averaging approach to solve the issue of sampling noise.
The article is organized as follows. First, the current operational version of the background error covariance matrix model in the YH4DVAR system is introduced in Section 2, together with the spatial structure of sampling noise and the characteristics of the local variation in the standard deviations for further investigations. In Section 3, filtering approaches are implemented on the operational system and the filtered results are obtained. The results of applying the filtering approaches to the operational system are presented in Section 4. Finally, conclusions follow in Section 5.
2. Motivations
In a 4D-Var system, the cost function measures the distance between the model trajectory and the observations over an assimilation window. The cost function [19] can be written in terms of increment bywhere subscript t is the time index, is the increment, and is the increment evolved by the tangent linear model from the initial time to time index t. and are the covariance matrices of background errors and observation errors at time index t, respectively. is the innovation vector at each time step, , where is the observation vector, is the background at time index t, and is the observation operator with a linear approximation .
To simplify the background term in , most leading NWP centers are modeling the B matrix in the CVT framework. A general CVT form of a covariance matrix can be written as an eigenvector decomposition [20]. In the VAR system, the B matrix can be modeled as , which is not computed explicitly in the assimilation but implied by the operator . The construction of operator is a challenge in background error covariance modeling, and it is currently decomposed into the parameter transform , horizontal transform , and vertical transform in operational settings [21]:where introduces balance constraints to the assimilations and the parameters of are thought to be nearly uncorrelated with each other [22–24]. Therefore, the background error covariance matrix can be modeled by , where the matrix is the covariance of .
2.1. Randomization Technique and Sampling Noise
In the YH4DVAR operational system, the control variables contain the vorticity , the divergence , the temperature and surface pressure , and the specific humidity . The matrix transforms into , and the parameters of correspond to the unbalanced components of the control variables; thus, is assumed to be a block-diagonal matrix with no correlation between the parameters. Therefore, the background error covariance matrix can be written explicitly by the autocovariance matrix and parameter transform , , with the background autocovariance matrix for each variable. In the VAR system, it is crucial to input the standard deviations of the background error of control variables into the autocovariance matrix to estimate the B matrix, and the standard deviations of other variables are implicitly given by the parameter transform.
Fisher and Courtier suggested that the randomization technique can be used to diagnose the standard deviations of background error in the assimilation system, thereby alleviating the impact of the imposition of balance conditions on the correlations of background error [15]. As the B matrix is of the form , the standard deviations of background error can be estimated by the sample standard deviations of model vectors, that is, applying to a set of N (N = 10) random vectors drawn from a Gaussian distribution with . Figure 1 shows the vorticity standard deviation field at approximately 1000 hPa obtained by the randomization technique. Figures 1(a)–1(f) correspond to the raw estimates obtained with 10, 50, 100, 200, 1000, and 10000 samples, respectively, which are abbreviated as , , , , , . The sampling noise tends to be characterized by a relatively small scale in such that the estimate strongly deteriorated, and the sampling noise is significantly reduced when the sample size increases to 50 that the main large-scale features of the vorticity standard deviation field are well captured. The values of the vorticity estimate are smaller in the tropics and data-dense areas and higher in the high latitude areas, especially in the southern hemisphere in . As the number of samples increases, the sampling noise is dramatically removed, and the scale of the standard deviation field is visible. When the number of samples reaches 1000, the basic shape of the estimate ceases to meaningfully change; for example, is similar to .
(a)
(b)
(c)
(d)
(e)
(f)
The randomization technique is widely used to diagnose the standard deviations of background error in VAR systems in many operational centers, and the sampling noise is inevitably introduced into the estimate when undersampling. Although increasing the sample number can alleviate this issue, it is usually unrealistic to employ a sample size over 100 for estimation in the operational system. It needs to be a compromise between the accuracy of the estimates and operational time constraints that considers maintaining affordable computational costs in the operational system. Therefore, it is important to remove the sampling noise efficiently in the case of undersampling.
To separate the sampling noise induced by the randomization technique from the estimate, it is necessary to investigate the difference between the sampling noise and the standard deviation field. The spatial structure of the sampling noise in the assimilation system will be explored in the operational version of the B matrix. The true state of background error standard deviations for vorticity is never exactly known, but it is known from Figure 1 that when the number of samples reaches 1000, the estimate of the vorticity field is sufficiently accurate. Thus, the reference for the vorticity standard deviation field is generated from 1000 samples. The sampling noise included in the estimates is obtained by subtracting the reference from the estimated values of 50, 100, and 200 samples. The sampling noise is defined by
The energy spectrum of the signal is an effective means to examine the spatial structure of the signal. The wavenumber in the spectral space corresponds to the scale of the spatial structure in the grid-point space that a large-scale spatial structure means that the energy of the small wavenumber in the spectral space is much larger than that of the large wavenumber, whereas it is a small scale. Figure 2 shows the sampling noise and the corresponding energy spectrum of a specific transect on the spatial fields from the 50 samples, 100 samples, and 200 samples. There are two striking features worth noting in Figure 2. First, the sampling noise is clearly presented in the vorticity standard deviation field and displays high-frequency oscillations around the globe on a small scale. From the corresponding energy spectrum in the right column, the energy spectrum of the noise over the entire band is more consistent and is close to the white noise. Another striking feature is that the structure distribution of sampling noise from different samples is similar, but as the number of samples increases, the absolute extreme value of the noise decreases, which is consistent with the fact that more samples generate a more accurate estimate. The information in Figures 1 and 2 indicates that the spatial structure of the sampling noise is different from that of the estimates, and this difference motivates the reduction in the sampling noise by filtering techniques.
(a)
(b)
(c)
(d)
(e)
(f)
2.2. Characteristics of the Standard Deviation Field
The characteristics of the standard deviations of background error obtained by the randomization technique in the assimilation system will be explored in the operational version of the B matrix.
The characteristics of the standard deviation estimate vary with the model parameters, as illustrated in Figure 3. An overview of Figure 3 shows that low values of vorticity, divergence, temperature, and U-wind estimates are found in the tropics, and higher values are found in the high latitude areas, especially in the southern hemisphere. It is also shown that the standard deviation fields of the four model parameters have different distributions: low values of temperature and divergence are more concentrated at low latitudes, and higher values are found in higher latitude areas; however, there are relatively large vorticity fields along the intertropical convergence zone, and U-wind has a larger value in the data-spare and Antarctic regions and a smaller value in the densely observed land area.
(a)
(b)
(c)
(d)
The temporal evolution of the standard deviation field is examined in Figure 4 at 1200 UTC and 2400 UTC. It appears that the evolution of the standard deviations over 12 hours is not quite significant, and the distribution of the temperature standard deviation field at 500 hPa (1000 hPa), stretching along the west-east axis, is quite similar over the period. In the VAR system, the covariances matrix is estimated by the climatic information; thus, the characteristics of the standard deviation field obtained by the randomization technique have a certain similarity. In addition, the distribution of the standard deviation field is connected with the dynamics of the synoptic situation. There is another striking feature in Figure 4 in the differences in the distribution of the standard deviation field at the vertical levels of the model. The values near the surface are larger than those at 500 hPa, which can be explained by the complicated physical process near the surface directly affecting the relative error statistics. The corresponding energy spectra of the standard deviation at 500 hPa and 1000 hPa are plotted in Figure 5. Both of the standard deviations at different model levels have relatively larger energy at large scales, and the energy spectrum of the standard deviation at 1000 hPa has a larger amplitude at medium scales than that at 500 hPa.
(a)
(b)
(c)
(d)
(a)
(b)
The length-scale of the estimate of the standard deviations varies with the model parameters and the vertical levels, and these variations have a profound impact on the procedure for eliminating sampling noise. Therefore, the method of removing the sampling noise should not only preserve the standard deviations of interest but also consider the variation in the standard deviation field with model parameters and levels.
3. Filtering Approaches
The standard deviations of background error estimated by the randomization technique are contaminated by sampling noise, and this noise has a negative effect on the estimate. The previous paragraph shows that the sampling noise differs in the spatial structure based on the standard deviations and that the standard deviations are model parameter-dependent and vertical level-dependent. In this section, the spectral filtering method used in the operational system and the concept of a spatial averaging approach will be introduced, with a discussion of the implementation of these approaches in the YH4DVAR system. The filtered results of spatial averaging on the estimate will also be discussed in detail.
3.1. Filtering Approaches and Implementation
In the YH4DVAR operational system, a spectral filtering method, truncated at a fixed wavenumber , is employed to remove the sampling noise by multiplying a coefficient with each wavenumber. The coefficient is defined bywhere is the wavenumber in the spectral space. The spectral filtering method achieves the elimination of sampling noise by setting a reasonable truncated wavenumber , based upon empirical settings, and the efficiency of the spectral filtering method is more significant when the standard deviation field is isotropic and varies slowly.
For the filtered results in grid-point space, the spatial averaging [25] of the estimate is defined bywhere , is the length of spatial averaging, and is the number of grid points in domain . The double sum in this equation indicates that spatial averaging has a potential effect of increasing the sampling realizations. This means that, in each model level, the total sampling realizations will be dramatically increased by , where is the sampling number. Thus, the averaged results have less sampling noise when the independent sampling realizations are increased.
It is important to investigate the implementation and computational cost of filtering approaches when considering their application to operational systems. The implementation of the spectral filtering and spatial averaging approach on the operational system is discussed, and the flowchart is shown in Figure 6. The practical implementation of the spectral filtering method using fixed wavenumber truncation in the operational system can be summarized by the following steps:(1)For a set of N (N = 50) random vectors drawn from a Gaussian distribution with , apply U to each vector to obtain a set of model vectors in spectral space, in which (2)Compute the standard deviations after transforming the model parameters into grid-point space(3)Transform the standard deviations into spectral space, then multiply the coefficient (4)Transform the filtered standard deviations into grid-point space
Spatial averaging can be performed directly in the grid-point space to obtain the filtered results after the standard deviations are obtained in the second step above; therefore, the spatial averaging approach does not need to transform standard deviations into spectral space and has a higher degree of parallelism in the grid-point space. There is another striking feature that once the optimal spatial average length is determined, there is no need to relocate the grid points satisfying the distance because the sample realizations in the field are invariable. The calculation of the optimal spatial averaging length L can be performed offline and updated with the frequency of the background error covariance statistics.
3.2. Filtering Effects
Figure 7(a) is the reference for the vorticity standard deviation field generated from 1000 samples, and the raw estimate computed from 50 samples is represented in Figure 7(b). It is possible to capture the main large-scale features with only 50 samples, and the sampling noise tends to be characterized by a relatively small scale, such that the details of the estimate strongly deteriorated. Note that in the operational assimilation system, it is usually unrealistic to employ a sample size over 100 for estimating because of the computational cost. The filtered result of the 50 samples is illustrated in Figures 7(c) and 7(d) and are abbreviated as and , respectively.
(a)
(b)
(c)
(d)
Compared to Figure 7(b), and show that the sampling noise is weakened in the averaged results, and the scale of the standard deviation is visible. The signal-to-noise ratios of and are 23.74 and 25.41, respectively, both of which are larger than that of , and a larger signal-to-noise ratio of strongly supports the fact that the spatial averaging is able to remove most of the detrimental noise. still has obvious noise, which is caused by the fixed truncation of the wavenumber. Figure 8 shows the coefficient of the fixed truncation of the wavenumber of the spectral filtering method. The used here is a conservative spectral filtering function. It sets most of the wavenumber with a coefficient equal to 1, which mainly cuts off some low-frequency noise and does not consider the local variation of the standard deviation field. The physical structure captured by is very detailed and is even as accurate as the estimate in . The averaged results have less sampling noise as the independent sampling realizations are increased, , which is coherent with the fact that the total sample size of the filtered 50 samples is close to the sample size of the raw 1000 samples. This attractive result shows the potential efficiency of the spatial averaging in that the averaging estimates can reach the accuracy of large sample estimates by significantly eliminating the noise.
Figure 9 shows the relative error of the filtering result of vorticity , temperature , and U-wind at 1000 hpa:where is the standard deviation filtered by the spatial averaging approach and is the reference standard deviation field by 1000 realizations. It can be seen that, with the increase in the filtering length, the relative error of the filtering result can be expressed as a U shape. There is a significant effect on reducing the relative error as the averaging length increases, and there is an optimal length at which the error of the filtering result is minimized. When the averaging length becomes excessive, the relative error increases again with a greater length. This can be explained by spatial filtering increasing the accuracy of the filtered results by increasing the independent sample realizations; when the length becomes excessive, the relevant points are also included in the realizations that result in a reduction in the filtering result. The optimal filtering length for vorticity is 250 km, for temperature is 750 km, and for U-wind is 500 km. The difference in the optimal length of these model parameters shows the difference in the standard deviation field, which is consistent with the characteristics of the standard deviation field mentioned in 2.2. It is absolutely necessary to consider the effects of the model parameters on the spatial averaging length in the procedure for removing sampling noise since the optimal filtering length is dependent on model parameters.
The efficiency of spatial averaging at the different vertical levels for different model parameters is shown in Figure 10. Figure 10(a) shows the optimal averaging length is different at the different model levels; for example, the optimal length of temperature around the surface is 750 km, which is larger than 500 km at 500 hpa. In addition, the trend in the optimal length varying with the vertical levels is different from that for model parameters. The optimal length of the U-wind increases obviously with the model levels, while the optimal length of the temperature decreases gradually from 1000 hPa to 500 hPa and then increases with the vertical model levels. In general, the optimal length of the variable increases with the vertical levels. This can be explained by the fact that the correlation length-scale of the variable field at higher levels is greater than that at the surface, and a greater length-scale corresponds to a greater optimal length for more independent sampling realizations. Figure 10(b) shows the relative error of U-wind raw estimates and optimal filtering results at different vertical levels. The relative error of raw estimates remains around 0.08 while there is a significant reduction in the relative error with the optimal length. The essence of the spatial averaging to reduce sampling noise is to increase the independent sampling realizations by optimal length. A larger optimal length may result in the fact that the sampling realizations are not completely independent, and this is coherent with the relative error associated with optimal length being increasing above level 64.
(a)
(b)
4. Effects on Analysis and Forecast Quality
The previous results show that the large-scale features of background error standard deviation can be captured by a randomization technique, but this technique may be contaminated by the sampling noise because of the finite sample sizes. Furthermore, the fixed wavenumber truncation in the spectral filtering method has an inferior impact on the local variation characteristics of the estimates, but the spatial averaging approach offers some advantages in addressing this issue. In this section, the effect of the spatial averaging approach on assimilation and forecast quality will be examined in the operational NWP system consisting of YH4DVAR and the global spectral model.
Applications for assimilation diagnostics and forecasting will be considered with the operational NWP system on a low-resolution grid with 91 levels (abbreviated as T399 L91). There are two groups. In one, the number N of random vectors is taken as 50 in the first group, which is numbered CN, using the spectral filtering method truncated with a fixed wavenumber to eliminate the sampling noise. In the other experimental group, the number N of random vectors is taken as 10 with spatial averaging, numbered SP. The results conducted with the operational system on a high-resolution grid with 91 levels (abbreviated as T799 L91) are used as a reference. The operational system run started at 00 UTC on 26 August 2015, ended at 00 UTC on 26 September 2015, and was updated every 12 hours. The statistical properties of the B matrix are generated using the popular National Meteorological Centre (NMC) method, an approach that uses the differences between forecasts of different lengths to evaluate forecasts error [9], and the optimal averaging length of the model parameters is calculated at 00 UTC on 26 August 2015 on each vertical level.
Figure 11 displays the root mean square error (RMSE) of the geopotential height as a function of time for two schemes: CN and SP. The Northern and southern hemisphere extratropics correspond to the regions of 30°N–90°N and 30°S–90°S, respectively. The RMSE in the southern hemisphere is generally higher than that in the northern hemisphere and tropics, and RMSE in the tropics is slightly smaller than that in the northern hemisphere. The RMSE of the SP analysis in all three regions was less than that of the control, which leads to the same conclusion as that for 850 hPa in Figure 11(b). It can be observed that the RMSE in different regions is significantly reduced by the use of the spatial averaging approach in Table 1, and this tends to be consistent with the result that the spatial filtering has a positive impact on improving the accuracy of the estimate.
(a)
(b)
The difference between the SP and the CN analysis fields at a level near 500 hPa and 850 hPa is shown in Figure 12. A visible common feature is that the differences in low latitudes were significantly smaller than those in high latitudes. This can be interpreted to mean that the equatorial waves coupled to convection explain most of the standard deviation in the tropics, especially in the lower troposphere; thus, the filtering effect of the spatial averaging approach is similar to that of the spectral method. Another interesting point is that the larger values are mainly located in the vicinity of ridges and troughs, which is in agreement with previous studies.
(a)
(b)
Applications for prediction have been achieved by examining the forecast quality of filtered results on a global spectral model. Figure 13 shows the RMSE of the forecast for temperature (500 hPa) at different times compared with the references at the corresponding time of T799. The RMSE has a small value in the tropics and northern hemispheres, whereas much higher values are found in the southern hemispheres in the data-poor areas. It is also observed that the RMSE in the three regions increases with a long forecast time, and the growth of RMSE in the tropics is slower than that in the nontropical area. The filtered estimates have a positive impact on the forecast quality, and its forecast improvement is generally weaker than that in the analysis. This may be due to the filtering effect of the numerical model, and continued research is important to improve the skill of the forecast.
(a)
(b)
(c)
(d)
(e)
(f)
5. Discussion and Conclusions
In a VAR system, the randomization technique is widely implemented in the current operational centers to diagnose the standard deviations of background error for modeling the B matrix. In this case, the standard deviations are estimated from a finite sample and will inevitably be affected by sampling noise because of the finite size of the samples. Therefore, some filtering techniques must be implemented to eliminate this sampling noise.
The spatial structure of the sampling noise introduced by the randomization technique is explored in the operational version of the B matrix. The estimated standard deviations from a finite sample are deteriorated by sampling noise, and the noise tends to be characterized by a relatively small scale and evenly distributed across the global map compared to the estimated reference field. The absolute value of noise dramatically decreases as the sample size increases, and the large-scale of standard deviations is more visible in the estimated field while the global distribution of sampling noise remains even. Such a difference in the structure between sampling noise and the estimates is the key to designing filtering techniques to filter out the small scale noise and preserve the signal of interest. There is a spectral filtering method truncated at a fixed wavenumber in an operational setting to resolve the issue of sampling noise in the YH4DVAR system. The application of this truncated spectral method is a kind of low-pass filter in spectral space that enables the small scale noise to be filtered out, preserves the large-scale information, and works well for eliminating sampling noise in the operational system.
Our investigation has also shown that the length-scale of the estimate of the standard deviations varies with the model parameters. The low values of the divergent standard deviations are more concentrated at low latitudes, and higher values are found in higher latitude areas, but there are relatively large vorticity fields along the intertropical convergence zone. In addition, the distribution of the standard deviation field is different between model vertical levels. The length-scale of the estimates from the randomization technique varies with the model parameters and vertical levels; thus, the fixed wavenumber truncation in the spectral filtering method has an inferior impact on the local variation characteristics of the standard deviations. The modified spatial averaging approach has some advantages in addressing this issue. In an operational context, the physical structure captured by the averaged results of 50 samples is fine, as the reference with the noise is dramatically eliminated, the spatial averaging approach has the potential to increase the independent sample realizations by averaging length, and the averaged results with the optimal length maximize the total sampling numbers. The optimal averaging length is connected with the model parameters and vertical level, which is consistent with the variation in the estimated standard deviations field; thus, the spatial averaging approach has better performance in eliminating the noise and preserving the characteristics of local variation in the estimates. In addition, the spatial averaging approach is easy to implement in the operational system and has a higher degree of parallelism in the grid-point space. The RMSE of the analysis in different regions is significantly reduced by the use of the spatial averaging approach, and relatively small values can be observed over data-spare regions such as the southern hemisphere. In addition, the results also show that the spatial averaging approach has a positive effect on improving forecast skill scores.
In the next step, we need to make further improvements to the spatial averaging approach. The optimal averaging length is computed offline, and the averaging length could be suboptimal with the variations in flow-dependent information in the assimilation period. It would be interesting to consider updating the optimal averaging length adapted with the variations. In addition, some qualitative analysis conclusions given in this paper need to be examined in detail in future work. The homogenous and uniform filtering type of spatial averaging approach ignores the anisotropic characteristics of the signal; thus, some important characteristics may be smoothed or weakened in averaging, and it would be worth the effort to introduce an adaptive averaging algorithm in the future.
Data Availability
The observations used in this paper are obtained from the NCEP and NOAA website (ftp://ftpprd.ncep.noaa.gov/pub/data/nccf/com/gfs/prod). The background data are not freely available because it is a short-term forecast by a global spectral model NWP system at the National University of Defense Technology that is subject to commercial confidentiality.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (Grant nos. 41675097, 41375113, and 41605070) and the National Key R&D Program of China (2018YFC1506704).