Abstract
With the development of meteorological observation technology, satellite data have found increasingly wide use in the numerical weather prediction field. However, there are various observational biases in satellite data, including a random bias brought about by complex weather systems and a systematic bias caused by the instrument itself, which greatly influence the quality of satellite data. A gradient information assimilation method is proposed in this paper to eliminate systematic bias. This method uses a gradient operator for gradient transformation between the model variable and observation variable and reaches the objective of eliminating systematic bias. An ideal experiment of variational data assimilation is conducted using the Community Radiative Transfer Model (CRTM) and Advanced Microwave Sounding Unit-A (AMSU-A) data, indicating that only assimilating gradient information can eliminate the smooth systematic bias in observation data. Then, a numerical simulation of tropical cyclone (TC) Megi and data assimilation experiment are conducted using the Weather Research Forecast (WRF) and WRF Data Assimilation (WRFDA) model as well as the Atmospheric Infrared Sounder (AIRS) data. The results show that the method of gradient information assimilation can improve the accuracy of TC tracks forecast and is also applicable for dealing with unreliable satellite data.
1. Introduction
With the rapid development of space technology and remote sensing technology, satellite data have been widely used in data assimilation, weather analysis, and forecasting by virtue of their high temporal and spatial resolution. In particular, satellite data play a key role in TC research because of a lack of observations over the ocean. For example, Kidder et al. [1] studied the warm core structure of TCs and gradient wind inversion with the data obtained by Advanced Microwave Sounding Units (AMSU). The study shows that the AMSU data provide improvement to track forecasts of TC [2]. Wang et al. [3] used BDA Method combined with AMSU-A Data Assimilation Method to reconstruct the Mesoscale Information of TC. McNally et al. [4] added AIRS data to the assimilation system of the European Centre for Medium-Range Weather Forecast (ECMWF) and found that assimilating AIRS data improved the quality of analysis and forecasting. The distribution of clouds, precipitation, and thermal structure in the TC’s evolution were described in detail using Cloudsat data [5]. Although satellite data have been successfully applied to many studies, the bias of the observation instruments on satellites can lead to most of the analytical errors in the model of the initial field and cannot be ignored [6]. Due to the bias of production accuracy as well as a variety of calibration, positioning, and other aspects of observatory instruments carried out by satellites, the final observation results also contain bias. At present, bias avoidance mainly includes quality control and bias correction. Quality control can eliminate obvious bias by deleting data that are too disperse or of low reliability but cannot eliminate bias from the source and can only reduce the use of unreliable satellite data. In addition, systemic bias is a problem of global data rather than individual data that cannot be eliminated by quality control. The purpose of bias correction is to eliminate or reduce systematic bias, the sources of which include bias in radiative transfer modes and instrument calibration, and systematic bias forms the response characteristics of the sensor over time. At present, systematic bias cannot be fundamentally removed by bias correction. Therefore, satellite observation data need to be revised before utilization. Wang et al. [7] proposed that, based on the GRAPES model variational system and assuming a known bias correction, forcing the initial state of the model mapping bias to be consistent with satellite observations of the air mass bias can achieve the effect of bias correction. However, in practice, the mechanism of generating bias is unknown, so it is difficult to eliminate the impact of systematic bias. The ECMWF initially adopted global scan correction and linear air mass correction to reduce satellite bias. However, the change of radiation angle and air mass modification were not taken into the consideration [8] in these methods, so Eyre [9] cited cloud radiation to adjust the scheme. Harris and Kelly [10] first considered the latitude dependency of scan correction and divided the globe into 18 latitudes. Then, air mass bias correction was performed on the result of the scan line scan correction using the air mass bias prediction factor from the mode background field. Liu [11], with reference to the previous TIROS Operational Vertical Sounder (TOVS) bias correction scheme of ECMWF combined with the radiation characteristics of the Chinese Advanced TIROS Operational Vertical Sounder (ATOVS), performed the correction of the scanning line bias and air mass bias, and the former was obtained by multiple linear regression modelling based on the selected air mass forecasting factor. The experimental results show that the error distribution tends to a Gaussian distribution before and after the bias correction and that the bias standard deviation is also reduced accordingly. Bao et al. [12] developed the FY3A-Microwave Temperature Sounder (FY3A-MWTS) satellite data bias correction system based on the classical satellite data bias correction method, which can effectively eliminate the air mass deviation. Furthermore, the deviation of the corrected image after averaging the observation residuals is greatly improved, approaching 0. At present, the bias correction method can reduce but not eliminate systematic bias. Wang et al. [13] proposed a method to eliminate systematic bias by using a gradient operator for gradient transformation between the model variable and observation variable based on a shallow water wave mode. This method makes the assimilation effect better able to absorb the spatial and temporal distribution tendency rather than the data itself and is also applicable to low confidence data. However, this method only studies shallow water equations and other ideal experiments and has not been used in practical applications.
With the popularization and rapid development of satellite observation data, the quality control of satellite data and elimination of systematic bias of observation data urgently need to be addressed. In this paper, the gradient information of AMSU-A observed TBB data is introduced into the model as new assimilation data using the CRTM model and the feasibility of the gradient information assimilation method for satellite observation data is verified by comparing assimilation results between gradient information assimilation method and original variational assimilation method. On this basis, the WRFDA module was used to assimilate AIRS observations with higher horizontal resolution. Then, the TC “Megi” was simulated with the WRF model. Finally, the results were analyzed by comparing the assimilation effect under different schemes. Although AMSU-A data and AIRS data are different in objections, channels, and resolutions, gradient information extraction and model variation replacement are only in the horizontal direction. Therefore, characteristics of these data do not interfere with the assimilation results. The AMSU-A observational data are used to verify the feasibility of the gradient information assimilation method. The advantage of assimilating AMSU-A observational data is that it avoids the complex channel selection problem and reduces the computational complexity of the assimilation algorithm.
2. Assimilation Principle of Satellite Data
2.1. Cost Function
In the conventional assimilation method, the cost function is defined aswhere is the difference between the mode control variable and the background variable and is the difference between the simulated and observed AIRS data. The detailed variable settings are as follows:
where represents background error covariance matrix, the superscript represents matrix transpose, represents weight coefficient matrix which reflects the quality of the observed data credibility, represents forward radiative transfer operator, and represents the brightness temperature of the observed data.
2.2. CRTM Mode
CRTM, a rapid radiation transmission model developed by the US Joint Center for Satellite Data Assimilation (JCSDA), is one of the software programs in the numerical weather forecasting data assimilation system. In this paper, the CRTM model is used to transform the existing atmospheric profile information (including barometric pressure, temperature, humidity, and ozone level) into the required observation brightness temperature information using the clear sky module (ClearSky), which acts as the observation operator in the ideal experiment. In addition, the descent algorithm adopted in the process of assimilation is the finite inner quasi-Newton method (Nocedal [14], 1980; Liu and Nocedal [15], 1989).
2.3. WRFDA
WRFDA is an assimilation module of the WRF mode that assimilates the observed data into a background field obtained by numerical prediction and uses the variational assimilation technique to obtain the required analysis field. The current WRF-3DVar system can assimilate conventional observations and unconventional data like radar, occultation data, satellite radiometric data, and so on.
3. Variational Assimilation Scheme of Gradient Information
3.1. Construction of Horizontal Gradient Information
Suppose there exists the following truth vector of satellite bright temperature observation data:
where is the number of data in the data sequence and superscript represents the transpose. Assuming there exists systematic bias and random bias vector of satellite data,and the bias has an unbiased stochastic distribution. The mode status variable is . The observation operator is . It is equivalent to use the CRTM mode to convert the profile containing the temperature information into the brightness temperature.
The observed data vector with bias can be expressed as :
To eliminate , it is necessary to rebuild the model variables and observation vectors. First, the model variable is transformed into gradient information. Then, consider the mapping of the model variable to the observation space vector.
Then, performing the gradient information transformation,and is the new pattern variable after the gradient information transformation. is the spatial distance between the two data locations. The observation vector is rewritten as :
is the gradient information after the transformation of the observed vector, which eliminates the impact of , and is the residual expectation of the new variables after transformation:
3.2. Gradient Information Assimilation
The gradient information is introduced into the assimilation of satellite data. When only the gradient information is assimilated, the cost function is defined as
is the difference between the new control variable and new observation variable after the gradient information is transformed.
represents the gradient information transformation operator.
4. Ideal Test
The standard profile used in this paper is from the CRTM model. According to the range of pressure in the actual model, it extracts the 63 layer temperature profile (51.5–1010 hPa) that contains the atmospheric temperature profiles, humidity profiles, ozone profiles, and other information. In this paper, the temperature variation in profile is chosen as the assimilation test. According to the distribution of the pressure layer and channel weight function map of the AMSU-A satellite data, it selects 4–9 channel as a channel assimilation which is also the mode of the commonly used AMSU-A satellite data assimilation channel.
Table 1 shows the specific test settings of the ideal scheme. The purpose of gradient information assimilation is to eliminate the systematic bias of the instrument. To highlight the impact of systematic bias and not consider the background error covariance matrix, a fixed random bias is set in the range of ±1 K. The systematic bias is set to 0, 0.5, 1, 1.5, or 2.0 K. According to the test settings, , represent the true value of the brightness temperature of two standard profiles. , represent the brightness temperature of the two standard profiles that are a superposition of random bias and systematic bias. , represent the brightness temperature of two profiles obtained by variational assimilation, which set , as the initial profile, and set , as the true value of the brightness temperature. Then, the RMSE of , , and , is calculated as an average of the 4–9 channel. The results are shown in Figure 1.
Figure 1 shows the RMSE of different assimilation schemes under different systematic bias. The RMSE of the AMSU-A_GRD test does not change when the systematic bias increases. The change of systematic bias has no effect on the assimilation effect under the condition that only the observed brightness temperature gradient information is assimilated, and it can eliminate the systematic bias of the instrument. In addition, the RMSE of the AMSU-A_GRD_C test was small compared to that of the AMSU-A test and can improve the AMSU-A test when the gradient information is used as a constraint condition.
Figure 2 shows the decline curve of the cost function of the three schemes when the random bias is constant and systematic bias is fixed at 1 K. The figure shows that the cost function of the three experimental schemes can ultimately be reduced to a stable value. The AMSU-A, AMSU-A_GRD, and AMSU-A_GRD-C test reach a stable value after 3, 5, and 10 iterations, respectively. All three schemes can meet the standard cost function decline curve, and the number of iterations reaches a steady value in a few iterations. The iterative efficiency of the gradient information assimilation method is close to that of conventional assimilation method, which shows that the method of gradient information assimilation is feasible.
Figure 3 shows the gradient of the decline curve of the three schemes, when the random bias is constant and systematic bias is fixed at 1 K. The gradient descending curve of the AMSU-A test achieves a stable value after the third iteration, that of the AMSU-A_GRD test has a large fluctuation in the 10th iteration and afterwards is stable, and that of the AMSU-A_GRD_C test is also nearly stable after the 10th iteration. Based on the above simulation results, it is found that the gradient descent curve of the gradient information assimilation method is close to that of the conventional assimilation method, and the feasibility of the gradient information assimilation method is verified.
5. Numerical Simulation of TC
5.1. TC Simulation and Assimilation Scheme
The number 13 TC “Megi” in 2010 was generated on the surface of the northwest Pacific Ocean at 20:00 on October 13, intensified to a strong tropical storm on the evening of October 14, and strengthened to a TC at 05:00 on October 15. It landed on China’s Fujian Province at Zhangpu County at 12 o’clock on October 23. The center of the TC had a maximum wind level of 13, a wind speed of 38 m/s, and a minimum pressure 970 hPa upon landing. In this paper, we simulate the initial moment at UTC October 17, 2010 06:00, and forward forecast 72 hours to October 20, 2010 06:00. The model background field is generated from National Centers for Environmental Prediction (NCEP) reanalysis data of 1° × 1°. The model uses a triple move nested grid, with grid sizes of 54 km, 18 km, and 6 km, and the number of horizontal grid dimensions is 100 × 70, 136 × 94, and 100 × 94, respectively. The first grid is fixed, and the second and third grids are vorticity following. The physical scheme adopts the Thompson graupel cloud microphysical scheme, RRTM longwave radiation scheme, Dudhia shortwave radiation scheme, Yonsei University (YSU) boundary layer scheme, and the Kain-Fritsch (new Eta) cumulus parameterization scheme. The assimilation time window is set to , and there are 88 channels into assimilation mode. The assimilation test is performed in the outer layer region. The initial conditions of the inner layer region are obtained from an outer region interpolation. The intensity, center position, and moving speed of the tropical cyclone are taken from the TC Yearbook compiled by the China Meteorological Administration. The specific test program set up is shown in Table 2.
5.2. TCs Simulation Path and Improvement Ratio
Figure 4 shows the comparison between the TC path simulated in each group and the actual conditions. Figure 4(a) shows that the path simulated by the AIRS_GRD test is the closest to the actual path in the movement trajectory trend and distance bias, and the AIRS test differs greatly from observation. The numerical simulation results of the TC are worse after assimilating AIRS data by the conventional assimilation method, but they are obviously improved when assimilating the same AIRS data using the gradient information assimilation method. Figure 4(b) shows that the TC paths simulated by the three experiments except the AIRS test are close to the actual path. Comparing the AIRS_SYB and AIRS tests indicates that the results of the numerical simulation are obviously improved by using the conventional assimilation method after adding a systematic bias to the AIRS data.
(a)
(b)
Figure 5 shows track errors of the TC from experiments during the 72 h hindcast period. Figure 5(a) shows that, compared to the CTRL test and AIRS test, after the assimilation of the AIRS data, the TC has a larger path error from 12 to 24 h and 18 to 48 h than the CTRL test without assimilation. Figure 5(b) shows that, in addition to the AIRS test, the other three types of path simulation bias are relatively close. In the four schemes, the path deviation of the AIRS test is large, especially from 18 to 48 h and 48 to 72 h, and the errors are obviously increased compared with the other schemes.
(a)
(b)
Figure 6 shows the improvement of the TC track. From Figures 6(a) and 6(b), it can be observed that the TC path simulated by the AIRS_GRD test is improved based on the CTRL test and AIRS test and that the amplitude and times of the improvement are obviously more than those of the worsening. In addition, the TC numerical simulation improvement of the AIRS_GRD test is particularly evident in the second half of the time period. As seen from Figure 6(c), comparing the AIRS_SYB test with the AIRS test, the AIRS_SYB test was significantly improved in the later numerical simulations after subtracting a systematic bias. Comparing the AIRS_SYB_GRD test with the AIRS_GRD test in Figure 6(d), the times and magnitude of the improvement and nonimprovement in the simulation results of the path are similar and can be ignored when compared to the improvement ratio between the AIRS_SYB test and AIRS test in Figure 6(c). Thus, the conventional assimilation method is improved after adding a specific systematic bias, while, for the gradient information assimilation method, the effect of the conventional assimilation result is negligible after adding a specific systematic bias that compares to conventional assimilation results. It can be seen from the above analysis that, for the conventional assimilation method, if there is a systematic bias in the assimilated observation data, then the numerical simulation results will have a greater impact. In addition, after removing this part of the systematic bias, the simulation results will obviously improve. For the gradient information assimilation method, even if there is systematic bias in the observed data, its influence on the numerical simulation is also negligible.
(a)
(b)
(c)
(d)
5.3. Improvement of Initial Field
Figure 7 shows the initial 500 hPa wind field deviation. Figures 7(a) and 7(b) mainly show the effect of gradient information assimilation method on wind field. As seen from Figure 7(a), after the assimilation of the AIRS data by the conventional assimilation method, the TC region has a westerly wind field deviation based on the CTRL experiment. The TC track is in the east-west direction. It affects the progress of the TC, so it makes the simulated TC track slower in the AIRS test than in the CTRL test, resulting in a large deviation of the path trajectory. As seen from Figure 7(b), near the TC, the AIRS_GRD test has an easterly wind field deviation compared with the AIRS test and the TC wind direction contributes to the TC moving speed. Thus, the conventional AIRS data assimilation method reduces the accuracy of TC simulation for TC “Megi” and the gradient information of assimilating AIRS data can improve the path simulation of TCs. Figures 7(c) and 7(d) mainly show the effect of a systematic bias in observed data on the wind field. As seen in Figure 7(c), near the TC the AIRS_SYB test has an easterly wind field deviation compared to the AIRS test, and combined with the TC moving path, this wind field deviation increases the speed of the TC. The simulated results of the AIRS_SYB test are better than those of the AIRS test in the path simulation of TC “Megi,” and the reason for the error in the simulation of the AIRS test path is that the early moving speed is slower than that of the other test. As seen in Figure 7(d), the wind field variation between the AIRS_SYB_GRD test and AIRS_GRD test is negligible compared to the wind field deviation between the AIRS_SYB test and AIRS test. Thus, adding a systematic bias to the observed brightness temperature when using the conventional assimilation method to assimilate the AIRS data significantly improves the initial wind field, while, in the gradient information assimilation method, adding a systematic deviation to the influence of the initial wind field is almost negligible.
(a)
(b)
(c)
(d)
Figure 8 is the initial bias of the geopotential height in the vertical section at 18.5°N. Figure 8(a) shows that, after the assimilation of AIRS data in the AIRS test, a negative geopotential height deviation occurs in the upper 250 hPa region, indicating that after the assimilation of the AIRS data, the geopotential height decreases, which makes the TC intensity of the upper layer increase. The geopotential height deviation is basically positive in the region below 600 hPa, which indicates that, after the assimilation of AIRS data in this region, the geopotential height increases and the corresponding TC intensity decreases. As seen from Figures 8(b) and 8(c), both the AIRS_GRD test and AIRS_SYB test show a positive deviation of the geopotential height above 600 hPa. The results show that compared with the conventional data assimilation method, the gradient data assimilation and conventional data assimilation method with certain systematic deviation increase the geopotential height in this region, and the TC intensity of the corresponding region decreases. In the lower 600 hPa region the deviation is negative, indicating that the geopotential height of the region has decreased, and the corresponding strength of the TC is enhanced. Based on the above analysis, the gradient information assimilation method can effectively reduce the influence of the systematic bias of observed data on the initial field of TCs and improve the accuracy of the numerical simulation results after assimilation.
(a)
(b)
(c)
5.4. Assimilation Diagnostic Analysis
Figure 9 is the difference of the number of assimilated field-of-view between the AIRS test and the AIRS_GRD test, in other words, the number of field-of-view of AIRS test entering the assimilation mode more than the AIRS_GRD test for each channel. Figure 9 shows that, in addition to a few differences between the channel being relatively large, most of the differences between the channels are located in the vicinity of 10. The gradient information assimilation method removed a small part of the visual field point further more based on the conventional assimilation method. In a sense, in the theory of gradient information assimilation method, some of the dead pixels of the field-of-view that do not meet the new variables after the gradient transformation are eliminated, and it played a role in quality control once again. For several channels with the large amplitude difference, at this stage the test cannot give a reasonable explanation, and it will be further studied in the next step.
Figure 10 is the OMA and OMB’s AVE from each channel of AIRS test and AIRS_GRD test. OMA is the difference between the observed field brightness temperature and the analysis field brightness temperature value, and OMB is the difference between the observed field brightness temperature value and the background field brightness temperature value. Figure 10(a) shows that, in the AIRS test, the mean values of OMA and OMB are distributed around 0. Because of the brightness temperature value of a single channel on a single field-of-view, under the conditions of systematic bias and random bias, the difference of brightness temperature value of each field-of-view between the observed field and the background field has both positive and negative values, the same as the difference between the observed field and the analysis field. In general, the OMA values in each channel are closer to 0 than the OMB values in the AIRS test, indicating that, after assimilation, the analysis field approaches the observation field on the basis of the background field, and it achieves the purpose of data assimilation. Figure 10(b) shows that, in the AIRS_GRD test, the AVE of OMA and OMB are all positive, indicating that the gradient information of brightness temperature in the observation field is larger than both of the analysis field’s and of the background field’s. The reason for this phenomenon may be that the new variable information obtained by the gradient information transformation in the gradient information assimilation method contains a certain degree of gradient information of the horizontal direction physical variable field, but the concrete reason needs to be further studied in the following work. In general, the OMA values in each channel are closer to 0 than the OMB values in the AIRS_GRD test, indicating that, after assimilation, the analysis field approaches the observation field on the basis of the background field, and it achieves the purpose of data assimilation.
(a)
(b)
Figure 11 is the OMA and OMB’s RMS from each channel of AIRS test and AIRS_GRD test. After assimilation, all OMA values of each channel are smaller than the OMB values, which indicates that the assimilation effect is achieved. In addition, the OMA of the conventional data assimilation method underwent an obvious decrease when compared to the value of OMB (Figure 11(a)), while the decrease of gradient information assimilation is relatively small (Figure 11(b)). This is because gradient information assimilation is mainly the trend of the assimilating observation field, not the physical quantity itself, so the decline of the OMA value is not obvious. Although the decline of the OMA value in the conventional assimilation method is more obvious, the value of the observed field also has great systematic bias and the reliability cannot be guaranteed, so the numerical simulation results in the latter part may not be a good improvement.
(a)
(b)
6. Comparison of Gradient Based Bias Correction with WRFDA VarBC Method
The VarBC module is a bias correction module which the WRFDA mode comes with, and it introduces the bias correction parameter . and the control variable construct a new control variable . The cost function is defined as represents background error covariance matrix, the superscript represents matrix transpose, represents weight coefficient matrix which reflects the quality of the observed data credibility, and represents forward radiative transfer operator. represents the brightness temperature of the observed data.
Assuming that the control variable vector background field is independent of the bias correction parameter background field, the cost function is defined as
Finally, the cost function is minimized by conjugate gradient method, and the bias of the observed data is corrected.
All of the tests mentioned above were not assimilated by the VarBC method. In order to compare the bias correction effect between the gradient information assimilation method and the VarBC method, a numerical simulation contrast test was added up. The test settings are shown in the Table 3.
Figure 12 shows the comparison between the TC path simulated in each group and the actual conditions. It shows that the path simulated by the AIRS_GRD test is the closest to the actual path in the movement trajectory trend and distance bias. By contrast, the AIRS test differs greatly from observation and AIRS_VarBC test result is in between. When comparing the AIRS_VarBC test and AIRS_GRD test results, path simulation results are obviously improved by replacing VarBC method with gradient information assimilation method.
Figure 13 shows track errors of the TC from experiments during the 72 h hindcast period. Compared to the AIRS_VarBC test and AIRS test, after adding VarBC method, the TC has a smaller path error from 12 to 48 h and 66 to 72 h than AIRS test. Nevertheless, the track errors of the gradient information assimilation method are still the least which is better than that of the VarBC method.
To be more intuitive, Figure 14 gives the improvement of the TC track. Comparing the AIRS_GRD test with the AIRS_VarBC test, the AIRS_GRD test was significantly improved in the middle of numerical simulation after replacing VarBC method with gradient information assimilation method. In summary, the effect of the gradient information assimilation method is better than that of the VarBC method that comes from WRFDA model.
The VarBC method introduces the bias correction parameter , and bias is corrected by the variational method. It can only minimize the value of all the observed data bias according to the correlation between the background field and the observation field. But for a single observation point, the effect of bias corrections is not the best. By contrast, the gradient information assimilation method proposed in this paper can eliminate the systematic bias of each observation point effectively.
7. Conclusion
In this paper, the ideal experiment of AMSU-A data is performed with the CRTM model to verify the feasibility of gradient information assimilation. The WRF model and WRFDA assimilation AIRS data are used to simulate the TC “Megi,” and the following conclusions are obtained:(1)In the CRTM ideal experiment, gradient information assimilation can eliminate the influence of systematic bias, and the descending curve of the objective function and gradient satisfies the basic requirement of assimilation, which proves the feasibility of gradient information assimilation of satellite data.(2)In the WRF model, it is found that the gradient information assimilation method can improve the TC wind field and bias of geopotential height in the vertical section and can improve TC path simulation.(3)In the assimilation experiment of AIRS data in the WRFDA module, the systematic bias of the observed brightness temperature has little influence on the gradient information assimilation method and is negligible compared with the conventional assimilation method. From this point of view, the method of gradient information assimilation has a positive effect in eliminating the systematic bias in the observed data.(4)In the WRFDA assimilation diagnostic analysis, the OMA (AVE and RMS) of the gradient information assimilation method did not undergo an obvious decrease compared to the conventional assimilation method. However, the emphasis of the gradient information assimilation method is on the assimilation of observed field trends and is still applicable to observed data with low reliability, which greatly increases the usefulness of satellite data.(5)TC path simulation results are obviously improved by replacing VarBC method with gradient information assimilation method.
In this paper, we study the positive impact of satellite data gradient information assimilation on TC numerical simulation. The next step is to extract and assimilate gradient information from the vertical gradient and time dimension.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work was financially funded by National Natural Science Foundation of China (Grant nos. 41375106, 11271195, and 41230421).