Abstract

Multipath is the main systematic error of the Global Navigation Satellite System (GNSS) short baseline positioning. Multipath cannot be eliminated by the double-differenced technique and is difficult to parameterize, which severely restrict the high-precision GNSS positioning application. Based on the spatiotemporal repeatability of multipath, the sidereal filtering in coordinate-domain (SF-CD), the sidereal filtering in observation-domain (SF-OD), and the multipath hemispherical map (MHM) can be used to mitigate the multipath in real-time. However, the multipath model with large matrix for multi-GNSS multipath mitigation is difficult to achieve lightweight calculation and the SF-CD cannot be applied to mitigate the multi-GNSS multipath. In this paper, we propose a new multipath mitigation strategy in the coordinate-domain that shakes off the formation mechanism of multipath, a CNN (convolutional neural network)-LSTM (long short-term memory) method is used to mine the deep multipath features in GNSS coordinate series. Furthermore, multipath will be mitigated in real-time by constantly predicting the value of the next epoch. The experimental results show that the CNN-LSTM effectively mitigates the multi-GNSS multipath. The method can reduce the average RMS (root-mean square) of multi-GNSS positioning errors in the east, north, and vertical directions by 62.3%, 70.8%, and 66.0%. Moreover, comparing with the SF-CD, SF-OD, and MHM, CNN-LSTM can more effectively mitigate the effects of the GPS multipath, and the ability of multipath mitigation is almost not affected over time.

1. Introduction

The Global Navigation Satellite System (GNSS), as a real-time and high-precision positioning technology, is widely used in many fields such as navigation, geodesy, deformation monitoring, and photogrammetry [14]. However, many errors limit the GNSS positioning accuracy [5], such as receiver clock error, satellite clock error, tropospheric delay error, ionospheric delay error, and multipath. The double-differenced technique can eliminate satellite and receiver clock errors and significantly weaken tropospheric and ionospheric delay errors in GNSS short baseline positioning; however, it cannot eliminate or weaken multipath. Therefore, the multipath is the primary errors source in GNSS short baseline positioning [6]. Research on multipath mitigation has become an important research issue, and many scholars have conducted lots of researches around this issue [710]. In addition to site selection, current researches can be divided into two categories: hardware-based and software-based approaches. The hardware-based approaches mainly mitigate multipath by improving antenna and receiver performance [1114]; however, that cannot eliminate multipath, theoretically [15].

There are many software-based approaches for multipath mitigation. Improved stochastic models based on the signal-to-noise ratio (SNR) and carrier-to-noise-power-density ratio are used to suppress the impact of multipath on baseline vector accuracy [16, 17]. Antenna array technique is an effective method to mitigate multipath. Ray et al. [10, 18] proposed to estimate the multipath based on the multiple close space antenna system, which significantly suppressed the influence of code and carrier multipath. However, the polarization of multipath signals is severe in complex and changeable environments [19]. By considering the multipath frequency-domain characteristics, signal processing methods are applied for mitigating the multipath, such as a Vondrak filter with cross-validation [9], empirical mode decomposition (EMD) [20], independent component analysis with reference signal [8], and wavelet analysis [21]. These signal processing methods provide ways to mine the characteristics of multipath and appear to mitigate multipath effectively.

Multipath is mainly affected by the satellite, receiver, and reflector positions. The operation of the satellites directly affects the multipath, when the receiver and reflector positions are relatively stable. The corresponding multipath also has a similar periodic variation considering that the satellite operating cycles are usually stable. According to this characteristic, the sidereal filtering (SF) was employed to mitigate the multipath [2225].

It is necessary to model different satellites separately since there are slight differences in the operating cycle of each satellite. The advanced sidereal filtering (ASF) was proposed to solve this problem. Based on the “zero mean” assumption [26], Zhong et al.[27] used the ASF to convert the double-difference residual series of the GPS carrier phase observations into a single-difference residual series, which was to reduce the multipath of different GPS satellite and to improve the coordinates accuracy. Similarly, Zhang et al.[28] and Ye et al.[29] successfully mitigate the multipath of three different types of satellites with the ASF method: GEO (geostationary orbits), IGSO (inclined geosynchronous orbits), and MEO (medium Earth orbits) in BeiDou Navigation Satellite System (BDS). Although the SF and ASF methods can effectively mitigate the impact of multipath on positioning accuracy, both methods have to calculate the operating cycles of the satellites in advance. However, the operating cycle of the satellite is not fixed, which changes slowly with time. Therefore, the degree of multipath mitigation in SF and ASF methods decreases continuously over time. Besides, when the multi-GNSS single-epoch solution is performed, the computational efficiency and reliability of the ASF method will be reduced due to the significant increase in the number of visible satellites.

Regarding the identified issues, Dong et al.[7] constructed a single-difference observation equation based on a single receiver with multiple antennas and proposed to establish a multipath hemispherical map (MHM) with satellite altitude and azimuth as independent variables. The MHM makes full utilization of the nonlinear relationships among the geometric relationship of the satellites, reflectors, stations, and multipath. In a stable environment, the MHM method appears to effectively mitigate multipath and directly displays the multipath of satellites at different positions. Wang et al.[30] proposed an improved MHM method, named trend surface analysis-based MHM (T-MHM), to improve the calibration accuracy of multipath in the lattice. Besides, MHM has been continuously being developed by scholars and being used for precise point positioning [31] and dynamic short-distance relative positioning [32]. The MHM and its improved methods overcome the dependence of the SF and its improved methods on the satellite operating cycle, and the operating mode is simpler and more effective than the SF method. However, when these methods are used to mitigate the multi-GNSS multipath, the multipath model will have a very large matrix, which makes it impossible to achieve lightweight calculation. Moreover, an increase in number of satellite and signal frequency types [33] by using multi-GNSS observation could lead the multipath to be more complicated [34].

Deep learning, as an effective tool for solving nonlinear problems [35, 36], has been successfully applied many fields in recent years, such as satellite clock error prediction [37], intelligent satellite selection [38], indoor positioning [39], GNSS reflectometry [40], and vehicle navigation [41]. The salient feature of deep learning is the data-driven with multiple layers of a perceptron. There is no need to establish an accurate physical model; only training and learning of historical data are required to obtain the optimal parameter characteristics of the data system, thereby completing data mining [42]. As mentioned above, the multipath has complex nonlinear relationships with satellites, receiver, and reflectors positions, especially in multi-GNSS observations. It is challenging to build a mathematical model of GNSS multipath due to this complex nonlinearity [43]. Therefore, we propose the use of deep learning methods for feature mining of multipath in GNSS coordinate series, real-time prediction, and mitigation of multipath. The specific deep learning method used in this paper is the combination of the convolutional neural network and long short-term memory (CNN-LSTM), which combined the advantages of local feature extraction of CNN and the prediction ability advantage of LSTM [44, 45].

The rest of this paper is organized as follows. Section 2 briefly introduces the principle of multipath and related processing methods. Section 3 designs different experiments to validate the performance of the proposed method. Finally, conclusions are drawn in Section 4.

2. Basic Principle and Algorithm

In this section, we firstly describe the mechanism of multipath generation, and then we introduce the basic principles of CNN and LSTM. Finally, we give the strategies and algorithms for real-time multipath mitigation using the CNN-LSTM method.

2.1. Basic Principle of Multipath

In GNSS measurement, the receivers not only receive direct signals transmitted by satellites but also inevitably receive signals reflected or diffracted by surrounding objects. The direct signals and reflected signals are superimposed to produce interference signals, which are treated as standard signals. Since the length of the propagation path of the reflected or the diffracted signals are usually longer than the direct signal, the interference signals will produce some deviations in the ranging results, which will severely damage the GNSS positioning accuracy [46].

For GNSS carrier phase observation, there is a phase delay between the reflected signal and the direct signal. At the same time, the signal energy will be attenuated after reflection. The degree of attenuation is related to the reflection factor of the reflector. When selecting a site, it is generally chosen to be far away from tall buildings. Therefore, one reflected signal SR can be expressed aswhere is the reflection factor, and , A, and are the phase of direct signal, the amplitude, and the phase offset of the reflected signal. The multipath Smth caused by a single reflected signal can be expressed aswhere is the wavelength of the satellite signal.

For a single GNSS satellite signal, it may be reflected by obstacles in different spatial locations. Meanwhile, the reflection factors of different obstacles are different, and the multipath errors caused are also different. When multi-GNSS combined observation, it can be seen from equation (2) that even the satellite signals of the same reflection path will have a difference in the value of multipath errors due to the signal wavelength. Especially GLONASS, the signals broadcast by each of its satellites have different frequencies.

Since the GPS satellite constellation in the space repeats nearly every sidereal day (23 h 56 m 04 s), the multipath errors will repeat approximately every sidereal day if the external environment and receiver location remain unchanged. Therefore, the periodic characteristic can be applied to model and mitigate multipath errors. However, for multi-GNSS combination or BDS with heterogeneous constellations, different navigation systems have different satellite orbital heights and satellite repetition periods. Therefore, multipath no longer has periodic characteristics and exhibit strong nonlinearity.

2.2. Convolutional Neural Network

CNN constructs multiple filters and extracts data structure features through convolution operations and pooling operations. The main types of layers to build a typical CNN architecture: input layer, convolutional layer, pooling layer, fully connected layer, and output layer. Among them, the convolutional layer comprises multiple feature maps, each of which is composed of multiple neurons, and each neuron is connected to a local area of the previous feature map through a convolution kernel. The pooling layer includes multiple feature maps, each of which uniquely corresponds to a feature map of the layer above it without changing the number of feature maps. There are usually several convolutional layers and pooling layers. The convolutional layers and pooling layers are alternately set. That is, one convolutional layer is connected to one pooling layer, and the pooling layer is connected to another convolutional layer.

In the convolutional layer, the feature map of the previous layer is convolved with the convolution kernel of the current layer. The result of the convolution operation is calculated by the activation function to form the feature map of this layer. The feature map of the output iswhere is the feature map corresponding to the j-th convolution kernel of the l-th convolution layer; is the activation function, usually using functions such as sigmoid and tanh; Mj is the acceptance domain of the current neuron; is the i-th weighting coefficient of the j-th convolution kernel of the l-th layer; and is the bias term.

The pooling layer uses local correlation to subsample the data to reduce the data dimension while retaining useful information. At the same time, the pooling operation is used to maintain the features, so that the features have a translation, rotation, and distortion invariance. The pooling operation can be expressed aswhere is the weighting coefficient of the down-sampling layer and is the down-sampling function.

2.3. Long Short-Term Memory

In 1986, Rumelhart et al [47] proposed a recurrent neural network (RNN), which uses a directed loop to establish relationships between input data at different times. Although RNN can effectively handle nonlinear time series, problems such as gradient vanishing and explosion still exist. Then, long short-term memory (LSTM) was proposed to effectively alleviate the above problems by introducing gate operations [48]. A typical LSTM architecture consists of a cell and three “regulators”, which form the information flow inside the LSTM unit: input gate, output gate, and forget gate. Cells are connected end to end and selectively memorize information to overcome the long-term dependence problem of nonlinear physical modeling, through the cooperative work of input gate, output gate, and forget gate.

The typical LSTM workflow consists of four parts: the decision to discard information, determining to update information, updating cell status, and outputting information.

2.3.1. Decision to Discard the Information

The forget gate determines how much information is discarded from the previous memory information stream. The output value of the hidden layer at time t-1 and the input value xt at time t are linearly combined, compressed using the sigmoid function to output a value between 0 and 1 (0 means completely discarded, 1 means completely reserved). The formula for discarding information ratio ft iswhere is the sigmoid function; is the weight of the forget gate; and is the bias of the forget gate.

2.3.2. Determining to Update Information

Determine what new information is stored in the cell state. There are two parts: the first part, the input gate determines the value that needs to be updated; in the second part, a tanh function creates a new candidate value vector and adds it to the cell state. The formulas for these two parts arewhere is the input gate; is the weight of the input gate; is the bias of update gate; tanh is the hyperbolic tangent function; is a candidate value; is the weight of update candidate; and is the bias of the update candidate.

2.3.3. Updating Cell Status

The old state of cell multiplies by discarding information ratio to determine what information needs to be discarded, and it is changed according to how much each status is updated as follows:where is the new status value and denotes elementwise multiplication.

2.3.4. Outputting Information

The output is based on the state of the cell, which is a filtered value. Run a sigmoid layer to determine which parts of the cell state should be output. Then, we place the cell state through tanh function and multiply it by the output of the sigmoid layer to output the determined parts. The formulas are:where is the updated output parts; is the weight of the updated output value; is the bias of the updated output value; and is the finalized output.

2.4. CNN-LSTM for Multipath Real-Time Mitigation

Based on the principle of the algorithms, CNN has a strong ability to mine local features of data, and LSTM can better obtain time-series features and has an excellent predictive ability. Therefore, we propose to use the CNN-LSTM integrated method to predict GNSS multipath (Figure 1). This work can be divided into two parts: one is to build a CNN-LSTM network model based on a certain amount of training data; the other is based on the established network model to perform real-time prediction and multipath mitigation. The prediction way of the method is shown as follows [49]:where denotes “predict,” is the coordinate value of the sliding window, and denotes the final prediction. Multipath is corrected in real time by constantly predicting the value of the next epoch.

The CNN-LSTM network model is trained on the training data, which is described as follows: firstly, we apply the convolution layers to perform convolution calculations on the data in each time series to extract local features and then apply the pooling layers to transform the output dimensions. ReLU function is applied in this part as the activation function for improving the generalization ability of the convolution calculations. Secondly, the local features are input into the LSTM network for deep mining and prediction of the data. To prevent overfitting, the dropout method and the fully connected layer (dense) are used to output the target vectors.

Based on the established CNN-LSTM network model, real-time prediction and multipath mitigation are described as follows: (1) set the currently obtained data as a simulation data stream by continuously adding new data to the training set; (2) set up a window to update the CNN-LSTM network model based on the previous network model and the new data and predict the multipath at the next moment; (3) according to the predicted multipath, the GNSS coordinate series at the corresponding time is corrected to mitigate the multipath in real-time.

Our hardware environment comprises a graphics processing unit (GPU) NVIDIA GTX 1060, a central processing unit (CPU) Intel Core i7-7700 K (processor frequency is 4.2 GHz, quad-core and eight-thread), and 16 GB of memory. Our algorithm is developed on Keras API (Application Programming Interface) framework based on Python 3.6 and Tensorflow framework.

When using the CNN-LSTM method for network training, the choice of hyperparameters will affect the final prediction result. Moreover, the data complexity determines the input size and output size of a typical time series prediction. To make the local characteristics of the prediction series more significant, we set the input to 100 and output to 1. Moreover, there is no deterministic principle about the optimal network structure, and the final hyperparameter selection is determined through multiple experiments [50].

Mean absolute percentage error (MAPE) was introduced as an evaluation indicator for hyperparameter selection. It can be expressed aswhere is the original series, is the predicted series, and is the length of the series. The value range of MAPE is [0, +∞), the smaller the value, the better the value selection of the hyperparameters.

We use the experimental data set in the next section to calculate and analyze MAPE. The results are shown in Table 1. The MAPE with a smaller value can be obtained with different schemes, and the difference between the MAPE values corresponding to the different schemes is not significant. It shows that using different hyperparameter schemes in Table 1 can obtain high-precision fitting results. Therefore, we choose the relatively optimal scheme 1 for subsequent data processing and analysis. Its parameters setting is as follows: the number of convolution filters is 3, the convolution kernel size is 12, the dropout is 0.1, and the number of hidden layers in the LSTM is 12.

3. Experiments and Results

3.1. Data Collection

The experimental data were collected on the roof of the School of Spatial Information and Geomatics Engineering, Anhui University of Science and Technology. Because of the short baseline distance, the small height difference, and the same receivers, we assume that the double-differenced technique can eliminate the influence of tropospheric delay and ionospheric delay in the coordinate results. Therefore, it is assumed that the coordinate series is only affected by noise and multipath.

The base and the rover stations simultaneously receive GPS, BDS, and GLONASS signals. The sampling interval is 1 s, and the satellite cutoff angle is 15°. To validate the feasibility of the CNN-LSTM method, we continuously observed between day of year (DOY) 300 and 308 in 2019 (9 days in total). We perform a single-frequency and single-epoch calculation with the baseline length constraint to ensure the ambiguity-fixed solution. Finally, we obtain a three-dimensional coordinate series for subsequent multipath mitigation. The general environment of the continuous operating observation station is shown in Figure 2. Moreover, the long-term static data are calculated as the true value to evaluate the solution result.

3.2. CNN-LSTM for Multi-GNSS Multipath Mitigation

To test the ability of the CNN-LSTM method to mitigate the multipath in multi-GNSS (GPS/BDS/GLONASS) series, we use the parameter model trained by DOY 300 to perform real-time multipath mitigation on DOY 301-308. Figure 3 shows the original series of DOY 301, the predicted series, and their residual series. The original series subtracts the true value to display their amplitude clearly. In their residuals (Figure 3(b)), the series can fluctuate at 0, indicating that the CNN-LSTM method can effectively mitigate the multipath. However, there is still lower frequency trend item interference in the residual. Therefore, it is necessary to analyze the multipath mitigation capability of the proposed method from the perspective of frequency-domain analysis.

The power spectral density (PSD) can reflect the relationship between signal power and frequency. The main interference frequency domain of multipath is generally considered to be lower than 0.02 Hz [51]. As shown in Figure 4, take the E direction series, for example, the DOY 301 original series and the real-time-corrected series are used for spectrum analysis. The low-frequency signal dominates the signal component of the original series. When it is mapped into the time domain, the multipath seriously affects the positioning accuracy. Obviously, in the mitigated multipath series, the low-frequency signal is significantly mitigated.

We count the results of the 8 consecutive days in the E, N, and U directions, respectively, and these statistical results are shown in Figure 5. The CNN-LSTM method effectively mitigates the multipath and obtains excellent positioning accuracy. Moreover, the proposed method is not subject to the degradation of positioning accuracy caused by the time passage. The 8-day improvements remained stable with high improvement results. The average accuracy of the original series is 1.32, 1.48, and 4.20 mm in the E, N, and U directions, respectively; the average accuracy after mitigating the multipath by the CNN-LSTM method is 0.50, 0.43, and 1.43 mm; the average improvement is 62.3%, 70.8%, and 66.0%.

3.3. Comparison of  CNN-LSTM  with Other Methods

In order to verify the superiority of the CNN-LSTM method, we conduct comparative experiments in the GPS system by different methods, such as the sidereal filtering in coordinate-domain (SF-CD), the sidereal filtering in observation-domain (SF-OD), and the multipath hemispherical map (MHM). It is worth mentioning that all multipath models are denoised by the EMD method. Moreover, the time delay of SF-CD is fixed at 236 seconds; the time delay of SF-OD is obtained by calculating ephemeris; the MHM method is based on the double-difference residual, because it is difficult for conventional receivers to remove the receiver clock error by single difference. Since the CNN-LSTM method uses the first day data for training and is used to predict subsequent data, for other methods, we also use the first day observation data to model the multipath.

Taking the MHM method as an example, we give the multipath model established by this method, as shown in Figure 6. Since the double-difference residual is used to map to the spatial position of satellites (elevation and azimuth angles), when the satellite with the highest altitude angle is used as the pivot satellite at some time, the satellite will not have the double-difference residual. The position of the pivot satellites is marked with a black frame in Figure 6. Obviously, the double-difference residual in the space domain cannot reflect the multipath spatial correlation, because the transformation of the pivot satellite causes the change of the double-difference residual, but that does not affect the MHM based on spatial repeatability to mitigate the multipath.

The following 8-day original series is corrected by the multipath model. Figure 7 shows the original series in the E direction and the corrected series of the four methods to analyze the degree of the multipath mitigation in the residual series. Among them, the SF-OD and MHM methods have an obvious mutation in the series because the pivot satellite transformation does not occur at a fixed time in each DOY, which leads to obvious “end effects” in the use of multipath models based on double-difference residuals by using the SF-OD method. That also results in obvious “distortion” in the double-difference residual value of the MHM method in selecting the nearest satellite position. This phenomenon can be avoided by using SF-CD and CNN-LSTM methods to mitigate the multipath in the coordinate-domain; however, the SF-CD method can only be used in navigation systems with the same orbit. The SF-CD and CNN-LSTM methods can significantly mitigate the multipath, but the mitigated results need to be analysed in the frequency-domain.

Figure 8 shows the power spectral density of the original series in the E direction and the mitigated series of the four methods. The CNN-LSTM method significantly mitigates the impact of multipath below 0.02 Hz, indicating that the CNN-LSTM has learned the deep characteristics of multipath in the network training stage. Only at the 2E-5 Hz, the ability of CNN-LSTM to mitigate multipath is slightly worse than that of SF-CD, which is the reason why there are a small amount of low-frequency multipath in the mitigated series in Figure 7. This does not affect the overall attenuation of the multipath because the magnitude of this effect can be ignored. Moreover, as shown in Figure 8, the mitigated results of the two methods based on the observation-domain are equivalent in the frequency-domain. The ability of MHM to handle low-frequency multipath is slightly better than SF-OD and slightly worse than SF-OD at high-frequency. The result obtained in the frequency-domain is consistent with the conclusion of the literature [7].

We mitigate the coordinate series for the following 8 days, and obtain the statistical results and improvement of the positioning accuracy, as shown in Figures 9 and 10. There is no significant difference in the RMS of the original 8-day series, but the RMS of the series mitigated by different methods has different rules. The correction results of SF-CD and SF-OD in DOY301-304 are basically the same, but the correction accuracy will decrease over time. The difference is that the accuracy attenuation of the series after SF-OD mitigation is mainly caused by the different time of the pivot satellite transformation. For the MHM method, the time lapse has little effect on the correction accuracy, and the accuracy decay is slower than that of the sidereal filtering methods. This shows that the spatial repeatability based on the multipath can effectively mitigate the multipath, even if only using the first day data for modeling. For the CNN-LSTM method, the corrected results are stable and will not decrease over time, and the three directions maintain excellent multipath mitigation capabilities, because the CNN-LSTM method can learn the deep characteristics of the multipath and can effectively predict the multipath.

Table 2 shows the average RMS of the original series and mitigated series from DOY 301 to 304. In the horizontal component, the correction accuracy of SF-CD, SF-OD, and MHM methods are basically the same. The CNN-LSTM method is better than these three methods, and the accuracy of the E and N directions of the original series has been increased from 2.61 mm and 3.35 mm to 0.80 mm and 1.04 mm. In the vertical component, the methods based on the coordinate-domain are significantly better than the methods based on the observation domain. The original positioning accuracy of the CNN-LSTM method is increased from 9.42 mm to 2.46 mm.

4. Conclusions

We consider the complex nonlinear change mechanism of multipath and the influencing factors, and we propose to use the deep learning method to conduct deep mining of multipath changing mechanism and establish a CNN-LSTM network model of multipath. Based on the established network model, a real-time mitigation strategy for multipath is proposed, and the feasibility and superiority of the new method are validated through different experiments.

The results show that the new method can dig deep into the change mechanism of multipath in coordinate series. Based on the trained model, the CNN-LSTM method accurately mitigates the impact of low-frequency multipath. The proposed method can effectively mitigate the multi-GNSS multipath in coordinate series.

Compared with the SF-CD, SF-OD, and MHM methods, the proposed method does not depend on the multipath repeatability that can be applied to the real-time multipath mitigation in the coordinate series. Also, the proposed method is more effective than the SF-CD, SF-OD, and MHM methods for multipath mitigation to a certain extent, and its improvement of multipath mitigation is not affected by the extension of data collecting time.

In general, the multipath mitigation based on deep learning proposed in this paper does not consider the formation mechanism of multipath but mines the inherent change mechanism of multipath in the coordinate domain, which has fewer constraints and stronger scalability. This paper provides a new solution for the multipath mitigation research.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Authors’ Contributions

Chao Liu and Yuan Tao were involved in the conceptualization; Chao Liu and Yuan Tao were responsible for the software; Yuan Tao and Tengfei Zhou were responsible for the methodology; Chao Liu and Yuan Tao validated the study; Chao Liu and Yuan Tao performed the formal analysis; Yuan Tao, Xingwang Zhao, Haojie Hu, and Haiqiang Xin investigated the study; Chao Liu and Yuan Tao were involved in collecting the resources; Yuan Tao and Chunyang Liu curated the data; Yuan Tao wrote the original draft; Chunyang Liu, Tianyang Chen, and Haiqiang Xin wrote, reviewed, and edited the draft; Xingwang Zhao, Haojie Hu, and Tengfei Zhou visualized the study; Tianyang Chen and Haiqiang Xin supervised the study; Chao Liu was responsible for the project administration; Chao Liu was involved in the funding acquisition. All the authors have read and agreed to the published version of the manuscript.

Acknowledgments

The work was supported by the National Key Research and Development Program of China (No. 2016YFC0803109), the National Natural Science Foundation of China (No. 41772130), and the Anhui University of Science and Technology Graduate Innovation Foundation (2019CX2077).