Abstract

The measurement of radiated emission (RE) in an anechoic chamber becomes very challenging at high frequencies, up to 60 GHz, because the scanning plane of the receiver is in measurement standard deviation from the actual wavefront. As a result, the RE intensity of the devices may be underestimated, resulting in electromagnetic interference. The deviation between the electric field at the far-field vertical scanning point and the actual wavefront is researched. Then, in an anechoic chamber, a hybrid deep learning amendment model of convolutional neural network (CNN) and transformer is proposed to correct the RE measurement at a 3 m distance. The results indicate that the correction is reliable, with an average error of 6.35% for a 3 m distance in a semianechoic chamber and less than 4.83% for other test scenarios. The proposed method provides a promising solution for RE measurement at a millimeter wave band in an anechoic chamber.

1. Introduction

Higher frequency bands, such as 24–28 GHz, 37–40 GHz, and even 67–71 GHz in 5G NR communication, will inevitably become more used with the exponential growth of wireless devices and ever-increasing bandwidth usage [1]. The radiation emission (RE) of electronic equipment is an important determinant in electromagnetic compatibility (EMC) tests. The highest fundamental frequency of the equipment under test (EUT) determines the measurement frequency ranges [2]. Several standards for measuring RE at frequencies up to 6 GHz [3] and 18 GHz [48] are described, but only a few standards for millimeter wave (mmW) bands are presented. However, with the emergence of an increasing number of electronic devices for 5G NR in the market, RE measurement in the mmW band has become crucial and has drawn the attention of EMC Standard Committees [9]. The measured distances between the vertical scanning plane and the EUT are specified in the RE test standard as 1 m, 3 m, and 10 m. Because of its greater accuracy and lower cost, an anechoic chamber with a distance of 3 m is regularly used at high frequency. When the measured frequencies are not very high, this measurement distance approximates far-field conditions (measurement distance greater than or equal to , where is the maximum size of the antenna and is the wavelength of interest) [3]. To capture the maximum signal of the EUT radiated field under the limited 3 dB coverage, the receive antenna must be moved from 1 m to 4 m in height along a vertical axis on the scanning surface as the electrical size of the EUT increases. As a result, as the frequency increases, the typical measurement distances for EUT may no longer meet the far-field condition, resulting in the scanning surface no longer overlapping the wavefront. As the operating frequency of EUT increases to the mmW band, the scanning position rises and the scanning surface moves away from the wavefront, consequently making the actual measurement distance exceed the specified one. These inconsistencies lead to the underestimation of the actual RE level of the EUT, which may pose safety risks during equipment use [10].

The aperture radiation source is replaced by an equivalent vertical dipole radiation model in this paper, which is based on the theory of line antenna theory. To solve the inverse radiation problem, we investigate the noncorrespondence between the equiphase surface and the vertical scanning surfaces during the RE measurement in an anechoic chamber, i.e., how to obtain the maximum radiated intensity on the wavefront using the far-field scanning technique. Despite significant advances in full-wave electromagnetic simulation methodologies, the electromagnetic analysis tools currently available are inadequate to address the inverse radiation problem. To address this limitation, this study proposes a hybrid deep learning model of the convolutional neural network (CNN) and transformer to correct the RE measurement. The transformer model encodes input sequences to solve sequential data problems and consists of an encoder and a decoder based on the self-attention mechanism [11, 12]. The transformer is used in computer vision for image classification and object detection due to its performance in sequence data processing [1315]and because of its ability to accurately and efficiently capture long-range dependencies between inputs and labels, it is also used for long-term serial data trend prediction. Despite its advantages over traditional neural network models, there have been few reports on its use to determine microwave radiation emission. As a result, the proposed hybrid deep learning method can be used to radiate emission data over continuous frequencies while predicting trends and correcting deviations.

The remainder of this study is organized as follows. We first review the literature on radiated emissions calculation and deviation correction using the deep learning method in Section 2, followed by the proposed radiated emission calculation model in Section 3, detailing the algorithm and field data deviation analysis in radiated emission measurement. The RE measurement correction and analysis were performed by a transformer-based algorithm presented in Section 4. Finally, we conclude our work in Section 5, along with future work.

The analysis of electromagnetic radiation is an essential work in the problem of electromagnetic compatibility with electronic devices. Almost no equipment can function without connection accessories such as cables and electronic connectors, as well as holes and slots in the enclosure that serve as paths for electromagnetic power leakage. Worse, some slits may become radiation antennas, resulting in electromagnetic radiation overload. Because different electromagnetic sources produce various radiation fields, a simple method of calculating the radiation field is to replace the radiated emissions with a set of tiny dipoles with the same radiated fields. An equivalent dipole model was used in the literature [16] to simulate the radiated emissions within an aperture enclosure. The authors proposed a method for representing electromagnetic emissions from a printed circuit board using a near-field scanning equivalent dipole model [17]. The authors used transmission-line theory and a dipole antenna model to present a novel method for estimating radiated emissions from electrically long printed circuit board traces [18]. Some other studies on RE measurement include discussions of the effects of test sites on the measurement results concerning site voltage standing wave ratio [19], absorber reflectivity and chamber size [20], measurement uncertainty, and the effects of site variability and measurement equipment [21, 22].

All these studies performed well in terms of calculating radiated emissions. More research studies on deviation correction using numerical simulations or measured radiation emission data combined with neural network algorithms are still needed. The researchers presented a deep learning approach based on a convolutional neural network combined with far-field wave data generated from a near-field resonant metal body at microwave frequencies for subwavelength imaging in the far-field [23]. In other studies, researchers used near-field scanning microscopy or an equivalent set of elemental dipoles methods associated with genetic algorithms [24], convolutional neural networks [25, 26], hierarchical attention-based deep neural networks [27], extreme gradient boosting method [28], or strategies based on artificial neural networks and optimizer algorithm [29, 30]. We use the transformer model to capture the relationship on feature maps to establish the correlations between two multivariate data series. The transformer is an encoder model based on the self-attention tool that breaks the limitation that the recurrent neural network cannot be calculated in parallel. As a result, it has got much attention. In [31], the authors proposed a Transformer-based long-term sequence prediction model that reduces the time and space complexity through the sparse self-attention mechanism. In [32], Wu et al. proposed a new decomposition structure with an autocorrelation mechanism, which embeds the sequence decomposition techniques to improve prediction accuracy. The studies mentioned earlier have demonstrated the great potential of the transformer model in the field of sequence data analysis. However, there are limited studies on error correction in the existing EMC domains for the radiated emission measurement. Therefore, this paper provides an in-depth discussion of the application of the CNN-transformer model in multivariate sequential data prediction and correction.

3. Field Deviation in RE Measurement

3.1. RE Calculation Model

Figure 1 briefly describes the RE measurement method using CISPR 16-2-3. The EUT was mounted on a turntable platform 0.8 m above the ground to measure, the radiated electric field receiving antenna was scanned from 1 m to 4 m and set for horizontal and vertical polarization to determine the maximum electric field value for EUTs with a size less than 3 dB (the beam width of the receiving antenna), and a receiving antenna center height equal to the EUT. The receiving antenna was scanned in the vertical/horizontal direction for EUTs with vertical/horizontal sizes greater than 3 dB. The field measured by the receiving antenna on the scanning surface (at point ) is not the actual field (at point Q) on the same wavefront, as shown in Figure 2. The fields measured by the receiver scanning along the vertical plane differ from those at the wavefront, as shown by points and Q. Because the propagation path (R) changes after reaching the scanning surface, the total electric field on the equiphase surface was attenuated. The actual propagation path causes more significant attenuation as the scanning height increases.

For simplicity, the radiated electric fields of apertures and cables of the EUT were assumed to be equivalent to those of vertical dipole antennas. The calculation model of a vertical dipole above the perfect electric conductor (PEC) is shown in Figure 2. The difference in electric fields between the actual wavefront and the scanning surface was calculated in the fully anechoic chamber (FAC) and the semianechoic chamber (SAC), respectively. When RE is measured in the FAC, the corresponding model only considers elemental dipole A. However, when measuring in the SAC, the ground reflection must be regarded as a PEC condition. Therefore, the radiated electric field is produced by superimposing the electric fields generated by dipole A and its mirror image elemental dipole B.

The total field of the dipoles was calculated by discretizing current elements with length at each point in space, and the electric field was the superposition of the electric fields generated by the current elements at that point [33, 34]. Electric fields and of the current elements of the dipoles A and B, respectively, are described in the following equation.where is the current of the current element, is the number of waves in the free spaces, is the angular frequency, is the distance between the short electric dipole and the field point, and is the angle between the line connecting the element to the field point and the z-axis. The variable can take the value of 1 or 2. Then, the total electric field radiated by the current element at the two symmetrical positions on dipoles A and B is shown in the following equation.where and represent the electric fields radiated by current elements on dipole A and and are those by current elements on dipole B. The electric field expression was converted from the spherical to the Cartesian coordinate system to obtain the field component received by the vertically polarized antenna, as shown in the following equation.

The distance from a single current element on dipole A to the measurement point and the distance from the corresponding current element on the mirror dipole B to the measurement point are depicted in the following equation.where are the coordinates of the short current element, are the coordinates of the field point, is the distance between the origin and the field point, and is the angle between the line connecting the origin to the point and the -axis. It is assumed in the calculation that the current distribution on the entire vertical dipole antenna satisfies the following equation.where is the peak value of the current on the dipole, is the wavelength, (=0.8 m) is the height of the dipole off the ground, and is its length. Thus, the total radiated electric field of the vertical dipole antenna at the field is given by the following equation.

The accuracy of the simulations based on the above equations was numerically verified. First, the electric field on the scanning surface of the dipole antenna (length = 1.6 m; operating frequency = 300 MHz) at a measurement distance of 3 m was calculated and compared with the numerical electromagnetic code (NEC) [35]. The dipole antenna feedback point was used as the coordinate origin and the axis of the antenna as the z-axis of a Cartesian coordinate system. Figure 3 indicates that the NEC calculation results correspond well with the calculation program results, with only a slight difference in the extrema near the feedback point.

3.2. RE Calculation Results and Analysis

The electric fields on the equiphase and scanning surfaces were calculated for the frequency band 1–60 GHz, with antenna lengths of 2 U and 6 U (1 U = 44.45 mm) in the FAC and SAC test sites, respectively, to analyze the effect of the 3 m and 10 m measurement distances on the RE measurement. Because the size of the current communication devices is in the 2 U–6 U range, antenna lengths of 2 U and 6 U were chosen. The vertical polarization receiving antenna was used as an example because the linear polarization antenna could measure either the vertical or horizontal electric field component. The maximum total field strength on the equiphase surface and the maximum component on the scanning surface were calculated at the same frequency point (coordinates shown in Figure 2). Then, equation (7) was used to calculate the relative difference. The scanning surface heights corresponding to the two maximum electric fields were recorded, as required in the RE test. The scanning height corresponding to the maximum total electric field on the equiphase surface is the height of the point extending along the radial line to the scanning surface, as the scanning height of point Q and the height of point P shown in Figure 2.where represent the maximum component of the electric field on the scanning surface and represent the total electric field on the equiphase surface.

Figure 4 depicts the difference in maximum electric fields between the scanning and equiphase surfaces and the corresponding scanning heights for the 3 m measurement distance and antenna length of 2 U in the FAC. For frequencies, less than 5 GHz, the maximum electric fields on the equiphase and scanning surfaces are roughly equal, as are the scanning heights. However, there is a deviation in both the maximum position and the maximum radiation value for frequencies greater than 5 GHz. As the frequency increases, so does the divergence in the maximum electric field on the scanning surface and the corresponding scanning height. The far-field condition of 3 m is satisfied for an antenna with a length of 2 U if the frequency is less than 56 GHz. However, as the frequency increases, the electrical length of the antenna gradually increases, as does the number of lobes. The maximum radiation direction is also constantly changing, resulting in the sawtooth distribution of scanning height shown in Figure 4(a). Meanwhile, as the frequency increases, the electrical length of the antenna increases, resulting in the maximum radiation deviation increasing as the path attenuation. As a result, the height of the maximum electric field on the scanning surface changes and becomes smaller than on the equiphase surface, as shown in Figure 4(b). Therefore, at more considerable scanning heights, the equiphase surface still behaves approximately as a plane but with a more significant error. Figure 4 indicates that the deviation is more than 15% when the frequency is higher than 5 GHz. When the frequency exceeds 25 GHz, the variation is greater than 20%, with a maximum relative difference of 42.89%. Thus, the deviation between the scanning surface and the equiphase surface for the RE measurement must be considered.

Figure 5 shows the results as the electrical size of the EUT increases and the equivalent radiating antenna length becomes 6 U. It can be seen from Figure 5 that as the size of the antenna increases, the electrical size and the number of antenna lobes increases. The maximum electric fields on the scanning and equiphase surfaces and the corresponding scanning heights are consistent only when the frequency is lower than 1.6 GHz. Compared with the antenna length 2 U in Figure 4, when the frequency is greater than 6 GHz, the measurement distance of 3 m no longer satisfies the far-field conditions. Consequently, the component at the field point cannot be ignored, resulting in a significant deviation in the electric fields on the two surfaces and changes in corresponding positions, as shown in Figure 5(a).

When the far-field condition is unsatisfied, the relative deviation in the maximum electric fields on the two surfaces can be as high as 43% (near 26 GHz) or as low as 2% (near 45 GHz), as shown in Figure 5(b). The reason is that, in the near-field region, the projections of and components on the z-axis can be either in the same or opposite directions depending on the changes in their electrical dimensions. When the frequency is above 44 GHz, the relative deviation in the maximum electric fields on the two surfaces is at most 10% at some frequency points because of the influence of the component. However, with the electrical size of the antenna getting longer, its maximum radiation direction exceeds the scanning height of 4 m. Therefore, the relative difference between the maximum electric fields on the two surfaces within the scanning height range becomes smaller. If the measurement site is the SAC, the positions of the maximum electric fields on the two surfaces change more drastically than those in the FAC, owing to ground reflection. As shown in Figure 6, the relative difference between the maximum electric fields of the two surfaces is more significant than 20%, with a maximum relative deviation of 47.89% in the frequency range close to 90%.

The results of the antenna with a length of 6 U in the FAC with a measurement distance increased from 3 m to 10 m are shown in Figure 7, where the frequency bands below 21 GHz satisfy the far-field condition. Figure 7(a) shows that the deterioration in the far-field condition is less than that when the measurement distance is 3 m, and the scanning surface is closer to the equiphase surface. The variations in the scanning height corresponding to the maximum electric fields on the two surfaces resemble the pattern in Figure 4. The scanning heights of the maximum electric fields on the two surfaces coincide in the frequency range of 1–7.5 GHz. Therefore, the path and polarization attenuation decreased compared with the distance of 3 m, and the relative difference between the maximum electric fields on the two surfaces was reduced. As Figure 7(b) indicates, the relative deviation is lower than 8.6%.

The results are summarized in Table 1 to show the difference between the maximum electric fields on the scanning and equiphase surfaces and the corresponding heights for a given measurement distance. It also demonstrates that the uncertainty caused by the mismatch of the equiphase and scanning surfaces cannot be ignored when measuring the RE in the FAC. The deviations grow as the electrical size of the EUT increases. They are worsened by the deterioration of the far-field conditions caused by the increase in frequency, which is more dramatic at a measurement distance of 3 m.

4. Transformer Algorithm-Based Correction of RE Measurement

4.1. Transformer Algorithm

The transformer model is based on the encoder-decoder architecture, including a multihead attention module and the feed-forward neural network layers. The residual connection and layer normalization are used in the encoder and decoder sections to prevent model degradation, as shown in Figure 8. The encoder generates the vector corresponding to the input sequence, while the decoder generates the target sequence according to the output of the encoder. The encoder maps the input sequence with N features to the linear embedding sequence data and then fed to a multihead attention layer. The decoder then generates the output sequence . To extract the local features of the input data, three layers of one-dimensional convolution operations are performed on the original input data before inputting the embedding layer, each layer with 128 neurons and a convolution kernel size of 11. . The convolution of the input data sequence captures its local features using local receptive fields and shared weights, so it achieves a certain degree of displacement, scale, and distortion invariance. On the other hand, the hierarchical structure of convolutional kernels only learns from simple low-level textures to higher-order semantic patterns. However, compared to the transformer model, they cannot capture long-term dependencies and require deeper networks with several layers to increase their receptive fields. Combining the efficiency and inductive priors of CNN with the ability of attention to capture long-range information can create a better architecture for sequential data applications.

Because information about the positions of the input data sequence is valuable, some information about the relative position of the tokens in the sequence must be injected to make the multihead attention layer aware of the sequence order. Positional encoding with fixed sine and cosine functions is used to identify position information, and position encoding was implemented in the following equation.where is the position in the embedding tokens, is the embedding depth index that takes a value in , and the is the embedding depth. The values generated by sine and cosine functions are concatenated pairwise and added to the embedding of the input sequence.

The multihead attention layer is the main layer where attention scores are calculated, and the attention scores of the input sequence are modeled using a self-attention mechanism based on three main concepts: the query vector , the key vector , and the value vector . A single sequence query searches potential relationships by finding similarities in the sequence through keys. The comparison between the query and key pairs gives weight to the value. The interaction between the attention weights and values determines how much focus to place on other parts of the sequence while representing the current sequential data. The query, key, and value matrices are calculated by multiplying the input convolutional sequence with three different weight matrices: , and .

In equation (9), represent query, key, and value, respectively. The multihead attention calculations are implemented by scaled dot-product attention, as shown in the left half of Figure 8. Scaled dot-product attention can be denoted as equation (10). The weight of value is calculated through query and key, which is used to determine a weighted sum of values.where represents the dimension of the key vector sequence; Figure 9 presents the calculation of multihead attention. In the multihead attention mechanism, the calculation of the ith attention head can be represented as equation (11). The multihead attention is the concatenation of each attention head.

In equation (11), , , and denote the linear transformation of of the ith attention head, respectively; the parametric matrix , , , and are the linear projection parameters; refers to the number of heads in the multihead attention; represents the dimension of key and query, respectively; Concat represents concatenate operation, and represents the linear transformation of concatenated output. The multihead attention model can focus on the information of different representation spaces and simulate more levels of detailed information.

The feed-forward layer comprises two linear transformations and a rectified linear unit (ReLU) activation function. ReLU is used for mitigating gradient vanishing and gradient explosion. The computation of the feed-forward network is position-wise, and the weight of the linear transformation of each step is identical. The formula of the feed-forward network can be given as shown in the following equation.

The last two layers in the decoder are an affine transformation layer, which consists of a learnable scaling factor and a bias factor, and a multilayer perceptron, which consists of two linear transformations. The affine transformation layer and the multilayer perceptron are sequentially added for the trend prediction of the sequence data. The corresponding formula is expressed aswhere , , and , are linearly learnable weights and biases, respectively.

4.2. Transformer Algorithm-Based Correction of RE Measurement

The proposed CNN-transformer algorithm is implemented with Python 3.7 and PyTorch 1.9.1 and runs on a Tesla P100 GPU server. It is necessary to determine the input and output data when using a deep learning algorithm for data error correction. According to the analysis in Section 3, the electric length of the dipole, the maximum electric fields on the scanning plane, and the location information strongly correlate with the maximum value of electric field intensity on the equiphase surface and its location information. Therefore, input data include the electric length of the dipole, the maximum electric field strength, and the scanning height. The output data contain the maximum electric field strength and height on the equiphase surface. Because the radiated emissions from the dipole with a length of 6 U include the results of the corresponding dipole length of 2 U, all the two 591 data sequences are divided into a training set, a validation set, and a test set. To construct the sequential dataset, one striding is used to generate sequential data with a length of 96 data points, yielding a total of 1054 sequential data, and the final data dimension is (1054, 128, 5). The first 900 sequences of data are used as the training set, the middle 25 data sequences are used as the validation set, and the last sequence is used as the test set. The sequence length between the validation and the test set is 128, ensuring that it does not participate in the training and testing process. To improve the correction accuracy of the measured data, a convolution kernel of size 11 was operated on the original data sequentially to extract the local features [3639]. Due to the measurement site differences (FAC and SAC) and the measurement distances (3 m and 10 m), after combining the 2 U and 6 U vertical dipole test data, there are 4 cases: (1) the measurement site is FAC, and the measurement distance is 3 m; (2) the site is SAC, and the distance is 3 m; (3) the site is FAC, and the distance is 10 m; (4) the site is SAC, and the distance is 10 m.

We investigated the prediction results of several deep learning models to evaluate the prediction accuracy of the proposed models for radiated emissions. First, we conducted experiments on several models, which include CNN-GRU, CNN-LSTM, and CNN-transformer [40]. Then, to evaluate the performance of these models, one sequential data including 128 samples were selected from the test data for prediction. The predicted results are shown in Figures 912, wherein the red, brown, and light green curves are the predicted values by the three different models based on the test data. As observed from the figure, deep learning models can predict future trends. Moreover, the predicted value can be used to solve the problem of deviations. The proposed model for the performance of prediction is assessed on bias, mean absolute error (MAE), and mean relative error (MRE) metrics where the mathematical equations are shown in equation (14). Basically, the MAE metrics calculate the absolute difference between predicted values and the actual and predicted values from the presented deep learning models.where and are the actual and predicted values, respectively.

The performance of the CNN-transformer model is quite good when compared with the other two models in terms of the curve morphological characteristics of the predicted values and the calculated bias, MAE, and MRE values. The MRE of the CNN-transformer is 4.83%, 6.35%, 2.61%, and 6.32% in different test sites, respectively, as demonstrated in Table 2. From the comparison figure of the three different models and Table 2, it can be seen that the CNN-transformer model has the highest degree of fitting between the true value and the predicted value. After comparing these three models horizontally and vertically, we find that the prediction accuracies of CNN-transformer models are higher than that of CNN-LSTM and CNN-GRU models.

Figures 9 and 10 reveal that the prediction results of the proposed model in the 3 m FAC are better than the prediction results of the SAC, and the two correspond well to local minimum values, but the deviation is near the maximum value. The overall deviation of the prediction results in the SAC test is greater, the maximum deviation from the true value is 20.79%, with an average deviation of 6.35%, while the maximum deviation relative to the real value in the FAC test is approximately 12.57%, and the average deviation is 4.83%. Although there is a small deviation between the predicted value and the real value for the local maximum value, the overall trend is consistent, indicating that the model has captured the relationship between the maximum electric field strength on the equiphase and scanning surfaces. Figures indicate that the model can be used to correct the maximum electric field on the equiphase surface based on the field value on the scanning surface with a distance of 3 m in SAC or FAC sites.

Figures 11 and 12 show that the prediction results of the presented model in the 10 m anechoic chambers correspond well to the real values and that the two are consistent between the local minimum and maximum values. The prediction results for the FAC test site indicate that the maximum deviation from the true value is 8.29%, with an average deviation of 2.61%; the results for the SAC test site only show a significant deviation at the last local maximum value, and the predicted deviation relative to the real value is about 15.69% and the average deviation is 4.40%. The results of the tests in the anechoic chambers are nearly identical. These figures indicate that the model can correct the maximum electric field on the equiphase surface based on the value on the scanning surface in the 10 m method full-wave and half-wave anechoic chambers.

From the visual analysis of the bias between the predicted data and the actual data, it can be found that the model can accurately predict radiated emissions from a vertical antenna for different test sites. In a predicted data series, the predicted data change with the actual data, and the local maximum and minimum values can be perfectly fitted. In contrast, it demonstrates the CNN-transformer model’s capability of modeling the long-term dependencies in sequential data. Meanwhile, the CNN-GRU and CNN-LSTM are slightly inadequate in fitting the long-term patterns of the sequential data.

5. Conclusions

The electric fields on the equiphase and scanning surfaces at different measurement distances in the frequency band 1–60 GHz were calculated and analyzed using the symmetric dipole antenna theory. By comparing the maximum electric fields on the scanning and equiphase surfaces as well as the corresponding scanning heights, it was found that when the measurement distance is 3 m, the noncoincidence between the equiphase and scanning surfaces and the deterioration of the far-field conditions caused by the frequency increase brought large uncertainties in the RE measurement. Using the hybrid deep learning algorithm of CNN and the transformer model to learn the relationship between the maximum electric field strength of the scanning surface and the values of the equiphase surface enables the rectification of discrepancies in the maximum electric field on the scanning surface. The error was about 6.35%, and the average error of the other scenarios was within 4.83%. The correction accuracy was relatively reliable. While the existing test site and method remain unchanged, the actual radiation of the EUT at high frequencies can be obtained through deep learning algorithm correction. Future research could focus on the RE measurement errors caused in annular aperture slits or multiaperture arrays and explore the use of other deep learning algorithms for deviation correction.

Data Availability

The radiated emissions data used to support the findings of this study are deposited in the Gitee repository (https://gitee.com/mikefengshi/conv-transformer).

Conflicts of Interest

The authors declare that there are no conflicts of interest.

Acknowledgments

This work was supported by the National Natural Science Foundation of China under Grant 61877041.