Abstract

The solar radio flux at 10.7cm (F10.7) is a direct monitor and an important indicator of solar variability, and F10.7 is commonly used in empirical atmospheric models, ionosphere models, etc. The source regions of F10.7 are mainly in the corona above the active regions, and the extreme ultraviolet (EUV) images reflect the coronal thermal structure. In this paper, an index is defined as based on the intensity values of solar EUV images to represent the coronal contribution to F10.7. The Spearman correlation coefficient between the observed values of F10.7 and is 0.85 in 304 Å EUV images. Based on the high correlation, an empirical model is constructed. Combining the EUV data of SDO/AIA and the twin STEREO/EUVI, solar full-disk EUV images can be generated, and the future 27-day values of can be calculated. Then, a realistic estimation of F10.7 from 1 to 27 days in advance can be provided by the empirical model. Compared to the predictive values of F10.7 by the 54th-order autoregressive models in 2012-2013, the error drop-rate of our model is 12.54%, and our method has significant advantages in the upcoming 3 to 27 days’ forecast.

1. Introduction

F10.7 is an index of solar radio flux at the frequency of 2800 MHz and the wavelength of 10.7 cm. There are primarily two kinds of radiation mechanisms: gyroresonance emission and bremsstrahlung emission [14]. Gyroresonance emission only occurs where the magnetic field is sufficiently strong. Schonfeld et al. [5] found that, on 2011 December 9, in the rising phase of Cycle 24, of the variable component of the F10.7 flux was associated with the gyroresonance emission mechanism, although this percentage varies significantly over the activity cycle. The strength of bremsstrahlung emission is related to the plasma density. Because the active regions and the flares are denser than the quiet Sun, the source regions of F10.7 are mainly in the coronal heights above the active regions [6]. Livingston et al. [7] found that there is no more the linear relation traditionally assumed to exist between F10.7 and the sunspot number during Cycle 24. Selhorst et al. [8] found that F10.7 and the number of active regions correlate well during the period of 1992-2013. It may be caused that F10.7 are more sensitive to weaker magnetic fields than those necessary to form sunspots, of the order of 1500 G. The observation of F10.7 started in 1947, and it has never been interrupted. At present, the commonly used value of F10.7 in the world is observed at local noontime (20:00 UT) by the Dominion Radio Astrophysical Observatory in Canada, and it is expressed in units of sfu (solar radiation flux density, 1 sfu = 10−22·W ·m−2·Hz−1) [9].

The solar extreme ultraviolet (EUV) irradiance heats and ionizes the Earth’s atmosphere. Because the EUV irradiance cannot be observed from the ground, F10.7 has been used as a proxy to represent EUV [5]. The upper atmospheric models have the computer code which also use F10.7 as the proxy for solar irradiances [10]. The model of upper atmosphere is commonly used to calculate the atmospheric density to the orbital prediction of low Earth orbit satellites, so it is necessary to input the future F10.7 in the case of orbit prediction [1113]. Therefore, the prediction of F10.7 plays a major role in the accuracy of orbit prediction.

The approximately 11-year long-period and average 27-day mid-period variations have already been discovered in solar activity. The 11-year periodicity is the solar variational cycle, and the 27-day periodicity is due to the period of the solar rotation [14]. Prediction of the solar radiation index should correspond with the periodic features of solar activity. Singular spectrum analysis (SSA) is a nonparametric spectral estimation method, which can be an aid in the decomposition of time series into a sum of components and forecast the time series by these components. SSA can be effectively used to analyze the time series with periodic oscillations [15]. In machine learning, an autoregressive model (AR-model) could learn from a series of timed steps and take the previous results as inputs for a regression model, in order to predict the value of the next time step [16]. Thus, Zhong et al. [17] and Liu et al. [18] used the SSA method and the 54th-order AR-model, respectively, to forecast 27-day F10.7 because both the SSA method and the AR-model are suited to periodic and cyclical statistics. The methods of these two articles both showed predictive precision higher than that achieved by the United States Air Force (USAF) during the solar minimum in the Cycle 23. However, when a large active region (AR) rapidly appears on or disappears from the solar disk, the predictive precision of the 54th-order AR-model is unsatisfactory. Additionally, Wang et al. [19] proposed a mathematical method to extend the mid-term prediction of F10.7 to as long as 54 days without the need for extra solar observation data.

F10.7 is generated at coronal heights and related to the presence of active regions and the occurrence of flares [8]. Some empirical prediction models based on major solar features are described below. Wen et al. [20] forecast F10.7 using the areas and positions of sunspots. Henney et al. [21] forecast F10.7 utilizing advanced predictions of the global solar magnetic field generated by a flux transport model. The two prediction models above were constructed only using data from the Earth-side solar hemisphere; many scholars have focused on the far-side solar disk. Quémerais and Bertaux [22] forecast future 14-day F10.7 using the interplanetary Lyman α background data obtained by the Solar Wind Anisotropies (SWAN) telescopes on board the Solar and Heliospheric Observatory (SOHO) [23]. Lindsey and Braun [24] used seismic maps of the far-side solar hemispheres for solar active regions (ARs) forecasting. Although the Lyman α background data and seismic maps of the far-side solar disk have the longer observed time span, the extreme ultraviolet (EUV) images of the twin Solar Terrestrial Relations Observatory (STEREO) spacecraft used in this study can show the far-side solar ARs more directly and describe their real-time variation characteristic more accurately.

The different wavelength EUV images are expected to observe the different regions of the solar atmosphere, such as the coronal hole (CH), quiet sun (QS), active region (AR), and flare plasma [25], and the values of F10.7 have a high correlativity with the magnetic active regions [26]. F10.7 can be visualized as viewing the chromosphere (electron temperature ~10,000K) through a lower corona of varying optical thickness due to the changing extent and density of the trapped coronal plasma overlying active regions and other active structures. For the discussion we formulate a new index, , which is a proxy for the low-coronal, free-free contribution to F10.7, derived from the EUV data of the Atmospheric Imaging Assembly (AIA) on board the Solar Dynamics Observatory (SDO) [27]. Balan et al. [28] demonstrated that the shorter (extreme ultraviolet and ultraviolet) and longer (10.7 cm) wavelength solar fluxes have a nonlinear relationship during entire solar cycles. This is because the decreasing presence of gyroresonant absorption compared with free-free emission at low activity levels and the contribution of the magnetic field strength below about 1500 Gauss to F10.7 is undercounted [29]. Therefore, we define a function of F10.7 in terms of . Furthermore, since February 2011, the entire far-side solar disk could be observed directly by the Extreme Ultraviolet Imager (EUVI) instruments, the telescopes in the two STEREO’s Sun Earth Connection Coronal and Heliospheric Investigation (SECCHI) imaging suite [30]. The solar full-disk EUV images can be obtained by combining data from STEREO/EUVI with that of SDO/AIA. Then, the future 27-day values of can be calculated to forecast F10.7. In this paper, an empirical method is constructed according to the ideas above to predict 27-day F10.7 based on the EUV images.

The aim of this paper is to predict solar F10.7 up to one Carrington rotation in advance utilizing full solar composite 304 Å EUV images from SDO/AIA and STEREO/EUVI. We first defined an index as based on the intensity values of solar EUV images to represent the low-coronal, free-free contribution to F10.7. This paper importantly demonstrates the tangible benefits that 360 degree solar observations provide for the prediction of solar activity. This prediction method is found to perform better in the 3-27 days period, particularly in cases where active regions emerge on the far side of the Sun.

After introducing the background of forecasting F10.7 in Section 1, Section 2 presents and analyzes the dataset and establishes the method. Section 3 contains the results and discussion, and Section 4 presents the conclusion.

2. Data Processing and Method Establishment

The two EUV channels (171 Å and 304 Å) of SDO/AIA are chosen, because they are also in STEREO/EUVI. The daily level-1 512 × 512 data FITS (Flexible Image Transport System) files of SDO/AIA from May 2010 to December 2015 are downloaded from the available database of the Joint Science Operations Center (JSOC) at Stanford University (http://jsoc.stanford.edu/), and the daily 512 × 512 data FITS files of STEREO/EUVI from January 2011 to December 2013 are downloaded from the available database of the website https://stereoftp.nascom.nasa.gov/data/beacon/ahead/secchi/img/euvi/.

The above data files are updated daily at approximately 20:00 UT, which corresponds to the observed time of the F10.7 index. The FITS files are eliminated if their quality (FITS header keyword) is not equal to 0. The data files of F10.7 can be downloaded from the websites of the National Oceanic and Atmospheric Administration (NOAA) (ftp://ftp.ngdc.noaa.gov/STP/space-weather/solar-data/solar-features/solar-radio/noontime-flux/penticton/penticton_observed/listings/listing_drao_noontime-flux-observed_daily.txt).

The SDO/AIA datasets from May 2010 to December 2015 are used to determine the form of the function between F10.7 and . The SDO/AIA datasets of 2012-2013 and STEREO/EUVI datasets are taken as the testing sets.

2.1. Coronal Contribution to F10.7

To match the EUVI data, the previous work should use the overlapping and interchangeable AIA and EUVI data to represent measurements from the same plasma [31]. This previous work has been developed as an automated procedure in the SSWIDL software routines (ssc_form_euvi_synoptic.pro). In this routine, to match the SECCHI/EUVI data, the SDO/AIA data has been processed by the following formulas:

For 171Å,

For 304Å,

The “” is the original SDO/AIA data; the “exptime” is the exposure time of SDO/AIA in seconds. The “” is date and time when observation of this image started. The parameter “1.1” is the rough conversion factor for AIA 171 Å images to EUVI 171 Å. The routine “ssc_get_aia_304_factor.pro” returns the rough conversion factor for AIA 304 Å images to EUVI 304 Å. Finally, the “” is the SDO/AIA data which matches the SECCHI/EUVI data. To calibrate the datasets of three EUV cameras, all SDO/AIA datasets have been processed by the above routine and all our results of SDO/AIA are based on the “”.

Vernazza et al. [32], Krista and Gallagher [33], and Pérez-Suárez et al. [34] found different solar regions (coronal holes, quiet sun, and active regions) corresponding to the different distributions in the intensity histogram of EUV images. Schonfeld et al. [5] showed that the EUV images collected by SDO/AIA could represent the bremsstrahlung component of F10.7, so a proxy is defined in our paper to represent the coronal contribution to F10.7 in EUV images by the following formula:where is the pixel intensity threshold of source regions (SRs). The pixel point belongs to an SR if is greater than or equal to . The parameter is the sum of pixel numbers on the Earth-side EUV images, so is the integral effect of all SRs in the Earth-side corona.

The Spearman Regression correlation coefficient () is defined as the Pearson correlation coefficient between the ranked variables, and assesses how well the relationship between two variables can be described using a monotonic function. is computed from where is the covariance of the rank variables. and are the standard deviations of the rank variables [35, 36].

There is a nonlinear relationship between and F10.7 [28], so the between them in different values is to define the values of in the EUV images of two channels. The results are shown in Figure 1. The max values are 0.75 in 171 Å and 0.85 in 304 Å.

The 171 Å passband mainly reflects the upper transition regions and the quiet coronal regions of the solar atmosphere [27]. The 304 Å passband shows the upper chromosphere/transition regions and less variability with solar activity on a short time scale [37]. F10.7 can be visualized as viewing the chromosphere (electron temperature ~10,000K) through a lower corona of varying optical thickness due to the changing extent and density of the trapped coronal plasma overlying active regions and other active structures [9]. The source regions of F10.7 are closer to the regions of solar 304 Å than 171 Å, so the agreement between the F10.7 and 304 Å is better than 171 Å. When is equal to 103.2 DN/s in 304 Å, the between and F10.7 is up to 0.85. Thus, the pixels where intensities are greater than 103.2 DN/s show some information about the F10.7 source regions.

Based on the above analysis, the index associated with the F10.7 source regions in 304 Å EUV images is defined as follows:

The green part in Figure 2(b) shows where the is on the EUV image. Most of this area covers the ARs, which is consistent with previous research: the source regions of F10.7 are mainly in the corona above the active regions [5, 6].

2.2. Establishing the Method

Figure 3 shows the scatter diagram between F10.7 and from May 2010 to December 2015. Considering the exponential relationship between and F10.7 in Figure 3, an empirical function of and F10.7 is defined as follows:where is a constant term that is expected to be the coronal base value. The parameters , , and are the undetermined coefficients.

To derive the best-fit coefficients for (4), all available daily and F10.7 from May 2012 to December 2015 are used for fitting by the nonlinear least square methods. The fitting function in Figure 3 is Y = 65 + 16.38 (lg X – 5.5) + 17.5 (lg X - 5.5)3.12. The parameter “5.5” is the logarithmic total intensities while there is no active region on solar surface. That means the minimum value of F10.7 is equal to 65 sfu while equals 5.5. The higher order term means that the rate of F10.7 increasing with growth is nonlinear. To reduce the error of fitting function, we add a linear term. The correlation coefficient between fitted Y and F10.7 is 0.86. In consideration of the degradation of EUV instruments, the parameters a1, a2, and a3 are calculated by sliding fitting the daily and F10.7 of previous 14 CRs (Carrington rotations, 1CR = 27 days) using the nonlinear least square methods. The correlation coefficient of sliding fitting is up to 0.92 (in Figures 10(a) and 10(d)): these are described in detail on the part three.

The STEREO/EUVI system provides a direct observation of the far side of the solar disk, showing areas which will rotate onto the side of the Sun visible from Earth in the next few days. There is a special IDL procedure for generating the full-disk EUV image (Figure 4) in the SolarSoft tree for STEREO ($SSW/stereo/ssc/idl/beacon/ ssc_form_euvi_synoptic.pro) that combines the 304 Å data from SDO/AIA with the data from STEREO/EUVI nearest to 20:00 UT. The “+” symbols under the numbers 0-27 in Figure 4 are the diurnal projection of the Earth on the solar surface along the Earth-Sun line from 27 March 2013 to 27 days later. The latitude and longitude of Earth’s projection can be calculated by the IDL procedure in the SolarSoft tree for STEREO ($SSW/stereo/gen/idl/spice/get_stereo_lonlat.pro). The green regions show where the pixel intensities are greater than 103.2 DN/s. The angle between the Earth-Sun line and the normal direction of pixel in the EUV image is defined as . While the range of belongs to [0°, 90°], the pixel is in the Earth-side disk. Then converting into the same coordinates with SDO/AIA images, the equivalent predictive of the next 27 days from EUV images is defined as in (5).where the parameter represents the day. The parameter is the intensity of pixel in the next days’ image.

The area and intensity of AR will change in the next 27 days, especially in the solar maximum. The features of an AR’s appearance, development, and disappearance are complex and unique. Therefore, the future values of are very difficult to forecast even though the far-side solar disk is observed before 13.5 days. So we compare the previous with on the same days to analyze the relationship between the and . For example, the value of on 27 March 2013 corresponds to the value of on 28 March 2013. The first scatter diagram at the top left corner of Figure 5 shows the linear relationship between the values of from 2 January 2011 to 28 February 2013 and the values of from 3 January 2011 to 1 March 2013. The last scatter diagram at the bottom right corner of Figure 5 shows the relationship between the values of from 2 January 2011 to 28 February 2013 and the values of from 29 January 2011 to 27 March 2013. Even though the points disperse and the linear correlation coefficients dwindle with the growing parameter t, we assumed there are 27 kinds of linear relationships between the values of and in the next 27 days, which is defined in (6). After adjusting for (6), the parameter is the corrected and predictive future 27-day values of .where and are the undetermined coefficients in the fitted linear equation (the green solid line in Figure 5), which are given in Table 1. The parameter t represents the day. Finally, the values of F10.7 can be predictive 27 days in advance through substitution of for in (4).

3. Results and Discussion

In the application, there is no data after the forecasting date. Only the daily and F10.7 before the forecasting date can be used to obtain the best-fit coefficients for (4). The correlation between Y and F10.7 is performed for a certain number of Carrington rotations (CRs) before the start date of the series. This date is then advanced one day at a time through 2012 and 2013, generating a new correlation on each day. The minimum, average, and maximum correlation from this series is then recovered for that fitting window. The fitting window is varied from 1 to 28 CRs to generate Figure 6. The minimal correlation coefficient () reaches the maximum in 14 CRs, and the three parameters , , and are steady after 14 CRs. Thus, the length of the fitting window is 14 CRs.

The 0-day predictive and observed F10.7 values are compared in Figure 7. The 0-day predictive values are consistent with the observed F10.7 values, especially from August 2012 to August 2013. The mean relative daily forecast error is defined as in (7). The 0-day forecast relative errors () of our method are 5.49% in 2012, 4.86% in 2013, and 5.18% in 2012-2013.where the parameter , which represents the number of testing samples, is equal to 611 in 2012-2013. The parameter represents the day, and and denote the observed and predictive values, respectively.

To show the precise of our prediction, Figure 8 shows a comparison of the 54th-order AR-model and our method with respect to the average daily predictive relative errors in advance of 1-27 days. The predictive F10.7 values of 54th-order AR-model are calculated with the method of Liu et al. [18].

The two kinds of increase with the growth of advance time. While the error of our model increases with a steady rate in 1-27 days, the other error grows faster in 1-9 days and approximates to horizontal lines in 10-27 days. From the contrast, it is discovered that there are two kinds of fitting errors in our model: one is in (4), and the other is in (6). in our method is the error of fitting in (4), which is up to 5.18%, and the of AR-model is below 4%. Thus, the F10.7 short-term predictive precision of our model is lower than that of AR-model. The reason is that the known preceding value has the tallest weight in the predicted value of AR-model, and the F10.7 varies slowly and autocorrelates very strongly. Therefore, the in our method is greater than that in the AR-model for 1-3 days. The growth rate of in our method is the error of fitting in (6). Even so, because of the addition of solar back disk data in our model, the accuracy of prediction is improved obviously for 4-27 days.

There are two classical evaluation indicators used to reflect the errors of the predictive results: the mean absolute error and the mean relative error. To assess the accuracy of every testing case, the errors are defined as follows:where and denote the observed and predictive values, respectively, and the parameter is the day.

To assess the accuracy of all testing cases, another two kinds of error are defined in (10) and (11).where the parameter represents the number of testing samples, represents the day, and and denote the observed and predictive values, respectively.

To assess the drop-rate between our model and AR-model, the parameters and are defined as follows:where the and represent and of AR-model, respectively, during 2012-2013. and represent the and of our model, respectively, during 2012-2013.

Figure 9 shows the observed values of F10.7 and the values of in the 54th-order AR-model and our method during 2012-2013. Due to the uncontrolled spin of satellites, instrument failure, or their inappropriate relative positions, SDO and the twin STEREO spacecraft cannot capture the full-disk EUV images clearly and completely, which should be omitted from the testing sample set (yellow vertical short line at the top in Figure 9). Finally, there are 611 samples for testing during 2012-2013, with 304 samples in 2012 and 307 samples in 2013.

The tendencies of in the two methods are approximated in Figure 9. The of our method is slightly less than that of AR-model, especially during June 2012 and October 2013. A threshold of is defined here to evaluate the accuracy of the forecasting results.

The in our method is less than 10% from July to November in 2012, which means that the prediction of our method is satisfied when the period of F10.7 is stable. Compared with the in the 54th-order AR-model, that in our method decreases from 11.54% to 10.09% with a 12.52% drop-rate during 2012-2013 (shown in Table 2). The drop-rate of is up to 10.35%. Table 2 shows that the accuracy of our method is desirable during testing. Additionally, two forecasting cases are chosen for further analysis in Figure 10.

The in our method is much less than that of the 54th-order AR-model in the case of 25 March 2013. Figure 10(a) shows a scatter diagram of observed and F10.7 in the previous 14 CRs and the future 27 days. The fitting line in Figure 10(a) is fitted by the data of the previous 14 CRs. The correlation coefficient (R) between fitted Y and F10.7 is as high as 0.92. Figure 10(b) shows a comparison of observed F10.7, , , and in the previous and future 27 days from the test date, 25 March 2013. Figure 10(c) shows a comparison of observed and predictive F10.7. The results of our method are closer to the observed F10.7 than that of the 54th-order AR-model. The ARs (NOAA 11711, 11715, 11716, 11717, and 11718) on the Earth-side SDO/AIA EUV images on 8 April 2013 (Figure 11(c)) were already captured 14 days in advance by the EUV images on 25 March 2013 (Figure 11(a)). Additionally, these ARs move towards the Earth without disappearance, and no new ARs appear. In summary, when the ARs on the far-side EUV images change slowly in the next Carrington rotation, the forecast accuracy of our method is preferable to that of AR-model. This is the greatest advantage of our model in comparison to the AR-model.

Another case on 28 February 2012 is shown in Figures 10(d), 10(e), and 10(f). There is a crest of from 1 to 14 March 2012 along with a trough of and (see Figure 10(e)). Contrasting the three images in Figure 12, new ARs (NOAA 11429 and 11430) appear, but they are too weak in the solar full-disk EUV image (Figure 12(a)). These ARs then strengthen while they rotate towards the Earth on 8 March 2012 (Figure 12(c)). of our method decreases from 15.65% on 28 February 2012 to 8.45% on 6 March 2012, as seen in Figure 13. Our model could recognize the approximate position of NOAA 11429 and 11430 on 4 March 2012 because the forecast values of F10.7 remain close to the observed values after 14 March 2012 (corresponding to 10, 9, and 8 on horizontal axis in Figures 13(d), 13(e), and 13(f), respectively). Thus, once the solar full-disk EUV image captures the new SRs, the forecast results of our model are adjusted rapidly, and the predictive precision is improved quickly. This is the second advantage of our model in comparison to the AR-model.

In comparison with AR-model in Figures 10(c) and 10(f), the accuracy of prediction is unsatisfactory in the coming 1-2 days. This finding agrees with the results in Figure 8. Thus, to gain a more satisfactory predictive precision, we can combine the upcoming 1-2 days predictive values of the 54th-order AR-model with the upcoming 3-27 days predictive values of our method in practice.

4. Conclusion

Through the error analysis and case studies on the medium-term forecast experiments of daily F10.7 during 2012-2013, we can obtain the following conclusions.

The full-disk EUV images can provide coronal information about the far-side solar disk, which is 13.5 days earlier than that of other models using only Earth-side information or F10.7 itself. Thus, the F10.7 medium-term forecast accuracy of our method is better than that of the 54th-order AR-model, especially for upcoming 3-27 days.

The sensitivity of our model is much higher than that of AR-model. The forecast results of our model can adjust rapidly, and the predictive precision is improved quickly. The drop-rate of in our method is 12.52% during 2012-2013.

This was the first attempt in which we found a proxy in solar EUV images to represent the coronal contribution to F10.7. It was also the first attempt to forecast the upcoming 27-day values of F10.7 based on the solar full-disk EUV images. Moreover, to gain a more satisfactory predictive precision, combining the upcoming 1-2 days predictive values of the 54th-order AR-model with the upcoming 3-27 days predictive values of our method in practice should be considered.

Although there are problems with receiving real-time data from satellites STEREO/EUVI and SDO/AIA, this paper importantly demonstrates the tangible benefits that 360 degree solar observations provide for the prediction of solar activity. The Lagrangian 5 (L5) point lies at the third corners of the equilateral triangles in the plane of orbit whose common base is the line between the centers of Sun and Earth. So the viewing angle in L5 point can reach -150 degrees in Stonyhurst heliographic coordinates and the L5 observing platform can provide the EUV images before about 11.25 day. If the current data could be afforded by the ability of the L5 observing platform, it would enable this technique to forecast about 11.25 days F10.7 in practice.

Data Availability

All data used in the manuscript can be downloaded from the available database of the websites. And below are these URL of websites. (1) The F10.7 index data used to support the findings of this study can be downloaded from the available database of the National Oceanic and Atmospheric Administration (NOAA) (ftp://ftp.ngdc.noaa.gov/STP/space-weather/solar-data/solar-features/solar-radio/noontime-flux/penticton/penticton_observed/listings/listing_drao_noontime-flux-observed_daily.txt). (2) The daily level-1 FITS (Flexible Image Transport System) files of SDO/AIA from May 2010 to December 2015 are downloaded from the available database of the Joint Science Operations Center (JSOC) at Stanford University (http://jsoc.stanford.edu/). (3) The daily FITS files of STEREO/EUVI from January 2011 to December 2013 are downloaded from the available database of STEREO Science Center (https://stereoftp.nascom.nasa.gov/data/beacon/ahead/secchi/img/euvi/ and https://stereoftp.nascom.nasa.gov/data/beacon/behind/secchi/img/euvi/).

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

The 10.7cm solar radio flux data are provided as a service by the National Research Council of Canada. We thank the SDO/AIA instrument team for providing the coronal observations and Stanford University’s JSOC for providing the available database. We also thank the STEREO/EUVI instrument team and NOAA of the USA for providing their data. This work was supported by the National Defense Science and Technology Innovation Special Zone, the National Natural Science Foundation of China (Grant no. Y75037A070 and Grant no. 41604149), the Beijing Municipal Science and Technology Project (project number Z181100002918004), and the Strategic Priority Research Program of Chinese Academy of Sciences (Grant no. XDA17010302).