Abstract

Landslides endanger regular industrial production and human safety. Displacement trend analysis gives us an explicit way to observe and forecast landslides. Although satellite-borne remote sensing methods such as synthetic aperture radar have gradually replaced manual measurement in detecting deformation trends, they fail to observe displacement in a north-south direction. Wireless low-cost GPS sensors have been developed to assist remote sensing methods in north-south deformation monitoring because of their high temporal resolution and wide usage. In our paper, a DLM-LSTM framework is developed to extract and predict north-south land deformation trends from meter accuracy GPS receivers. A dynamic linear model is introduced to model the relation between measurement and the state vector, including the trend, periodic variation, and autoregressive factors in a discontinuous low-cost latitude time series. The deformation trend with submeter-level accuracy is extracted by a Kalman filter and smoother. With validated input as in previous work, the power of an LSTM network is also shown in its ability to predict deformation trends in submeter-level accuracy. A submeter-level deformation trend is detected from wireless low-cost GPS sensors with meter-level navigation error. The framework will have broad application prospects in geological disaster monitoring.

1. Introduction

Due to long-term mining, landslides threaten the safety of coal mine workers and nearby residents [1]. Continuous deformation on the surface of a coal mine is the main cause of landslides. The surface deforms slowly or relatively steadily most of the time [2]. Once the displacement accumulates and exceeds the maximum tolerance of the land surface, large slope failure becomes unstoppable [3, 4]. Constant deformation monitoring captures necessary information for landslide research. After summarizing the laws of long time land surface motion, early warnings can be communicated to local people [5].

Coal mining in the Fushun Western Open-Pit Coal Mine (FWOCM), which is located at 41°5038.33 north latitude and 123°5321.77 east longitude, started in the early 20th century. Production of coal used to be a pillar of the Fushun economy. Decades of mining has caused frequent surface deformation and landslides that threaten the city, Fushun. An open pit with 6600 m east-to-west length and 2200 m north-to-south width stands among the residential districts [6, 7]. Therefore, a landslide is a serious threat to miners and local people. According to our field investigation, deformation in the north-south direction is much worse than in the east-west direction. The velocity of displacement on the north and south sides is less than 3 m/year. A fault zone was spotted between the coal mine’s conveyor belts and residential buildings. To protect the buildings and people nearby, an effective approach is required to detect north-south land deformation.

For large-scale slopes, synthetic aperture radar (SAR) and differential interferometric synthetic aperture radar (DInSAR) are typical area surveying techniques [8, 9] with submeter- or centimeter-level accuracy. SAR systems provide full-time and full-weather observations for land motion. However, effective mine surface deformation monitoring requires high spatial resolution remote sensing images, which are expensive and hard to obtain. Direct temporal prediction of landslides is, however, scarce and constrained [10, 11]. Moreover, if the slope does not correspond with the ascending or descending passes of SAR satellites, it is concentrated in a few pixels of the satellite image. Therefore, traditional SAR or DInSAR can detect submeter or centimeter displacements in up-down and east-west directions (Figure 4 in [12, 13]), but have difficulty in observing north-south ground motion. Although a number of studies have engaged in modifying conventional SAR image processing algorithms [14] to detect north-south deformation, SAR and DInSAR solutions are still too sophisticated for practice uses. Therefore, north-south deformation time series should be observed directly.

Time series constructed from global positioning system (GPS) observation provide indispensable information about north-south land surface deformation. Compared with large-scale observation, such as SAR, GPS receivers have the advantage of high time resolution and north-south displacement observation. In past cases, GPS latitude time series of one location were obtained by reading the value periodically. These GPS series supplement remote sensing approaches in north-south deformation detection. With the rapid development of information and communication technology (ICT), interconnected GPS sensors form wireless sensor networks that can substitute manual measurement in long-term, high temporal resolution land deformation monitoring. Spatially distributed autonomous sensor systems generate massive data streams while saving plenty of labor costs.

However, time series obtained from wireless GPS sensors suffer in terms of both accuracy and continuity. Current wireless GPS sensor networks use low-cost GPS receivers that are widely found in smartphones, vehicles, and unmanned aerial vehicles (UAV). They receive satellite C/A (coarse acquisition) codes on the L1 frequency (1.57542 GHz). Under standard position services (SPS), the accuracy of a stand-alone C/A code receiver that uses only satellites is within 10 meters at the metric of 95% all-in-view horizontal, which exhibits significant errors in positioning. On the other hand, because most wireless GPS sensors are deployed in rural areas, unstable power supply becomes inevitable. Limited energy leads to GPS sensor dysfunction, which generates discontinuous points in latitude time series. At the current stage, problems in both sensor materials and power supply cannot be solved immediately. Therefore, an advanced processing paradigm is required to exploit a low-cost, discontinuous GPS latitude stream.

Submeter displacement monitoring is required at FWOCM. Real-time differential GPS techniques, such as real-time kinematic (RTK) [15], improve positional accuracy through sending corrections from a base station to an RTK terminal in real time. Network RTK (NRTK) uses several widely spaced permanent stations, named continuously operating reference stations (CORS), to provide correction information to the RTK user terminals. NRTK improves the accuracy, reliability, and availability of RTK. Centimeter or millimeter accuracy of displacement measurements can be obtained with NRTK. However, the cost of establishing RTK and NRTK and their operation is high, which hinders the application of RTK and NRTK. Combining only low-cost GPS receivers has demonstrated another solution [16, 17]. Averaging, linear filters and maximum likelihood estimation have been used to process the homogeneous GPS receivers with the aim of improving accuracy. The paper [18] suggests that using a dense array of low-cost GPS receivers may achieve high-precision displacement measurements. It also suggests that removing the fluctuation may improve the displacement measurement of low-cost GPS receivers.

Like most time series, the 21-month-long time series of the low-cost GPS sensors deployed at FWOCM can be assumed to consist of three components: trend, seasonal variation, and unsystematic irregular fluctuations. For noisy fragmented time series, the trend expresses the dynamic movement explicitly, and thus, we put trend analysis as one of our overriding tasks in land deformation monitoring using low-cost GPS sensors. Recently, methods proposed to trend analysis can be summarized as two-step procedures, which are trend extraction and prediction. The unprocessed, raw time series is firstly smoothed to obtain validated trends. Then the extracted trends are used as training samples to estimate future trends.

There is a wealth of literature concerning the trend extraction and prediction of low-cost GPS time series [18]. However, few comprehensive trend analysis frameworks have been developed. For trend extraction, researchers firstly utilized Kalman filters [19]. Although the Kalman filter has good performance in denoising Gaussian noise, the simplicity of Kalman state space function limits its power in filtering out complex periodical variations. To address this issue, empirical mode decomposition (EMD) is introduced in Gao et al.’s paper [20, 21] to enhance the Kalman filter in nonlinear and non-Gaussian fluctuation removal. However, the unexplained intrinsic mode functions (IMF) generated during EMD sifting makes it less reliable. In the mean time, experiments in early warning of abrupt displacement change at the Yemaomian Landslide of the Three Gorge Region, China [22], shows that displacement fluctuations can be detected by an autoregressive model (AR), which is seldom considered in previous Kalman or EMD approaches. Thus, a unified model that combines spline smoothing, periodical variation detection, and AR process estimations should be proposed.

For trends prediction, recurrent neural networks (RNNs), particularly with long short-term memory (LSTM) hidden units, are widely acknowledged to be powerful state-of-art prediction models for sensor time series [23]. However, LSTM has rarely been shown in recent GPS land deformation trend analysis papers. Moreover, most LSTM networks for land deformation trend prediction received unprocessed noisy raw data as input, which confuses the network with wrong prior knowledge. Since the Kalman filter has been proved to enhance LSTM networks’ performance [24], it is natural for us to think about applying LSTM networks with a unified trend extraction model as a preprocessed step.

This paper studies the feasibility and practical implementation of a unified comprehensive framework that can satisfy both land deformation trend extraction and prediction from low-cost GPS sensors. We model slowly moving displacement, which means annual movement of less than 2 meters [2], as variations of mean level in daily measurements, and treat collected time series as samples from the distribution of each measurement. In this case, all of the discontinuity within a day can be avoided. Then, by refining the dynamic linear model (DLM) proposed in [25] and combining it with LSTM networks, we develop an end-to-end DLM-LSTM framework to solve trend analysis in a discontinuous low-cost GPS sensor time series. The major contribution of our framework lies as follows: (1)The introduced DLM model expands both the observation matrix and the evolution matrix of a state space function to model periodical variation and autocorrelation in residuals. These expansions allow DLM to extract displacement trends with higher precision from noisy, discontinuous, periodically variable, and autoregressive low-cost GPS time series.(2)For sequence prediction, the introduced LSTM networks receive training-preprocessed samples. Those 5training samples provide validated prior knowledge of land deformation trends, which enhances LSTM networks’ power in both linear and nonlinear time series prediction. LSTM networks will receive a lower mean square error, which means that prediction sequences will be closer to real land deformation trends.(3)The fusion of DLM and LSTM provides an end-to-end north-south land deformation trend analysis workflow, which can undertake both trend extraction and prediction jobs. Our framework is able to model discontinuous low-cost GPS time series without other operations. The raw sensor data received by wireless sensor networks can be used directly as input.

To test our DLM-LSTM framework, a spatial-temporal distributed low-cost latitude time series (LCLTS) dataset is constructed using deployed wireless GPS sensors on both the moving southern slope of the Fushun Western Open-Pit Coal Mine (FWOCM) and the static comparative field in Baoxie, Wuhan, China. All GPS sensors that provide latitude time series are of the same type. Experiments have shown that our framework extracts and predicts submeter-level north-south deformation trends in FWOCM from February to March in 2016, which coincides with manual measurement over the same period. However, the trend extracted in the Baoxie field from November to December 2016 is static at the submeter level. The two results validate our framework and prove its time and space independence. Like [26], our sensor network helps SAR in high-resolution north-south deformation detection. Steady observations reveal local features of 2-dimensional mine surface motion, which are valuable resources for researchers and authorities.

The paper is organized as follows. In Section 2, the DLM-LSTM framework is presented. The processing results for the LCLTS dataset are shown in Section 3. The conclusion is given in Section 4.

2. DLM-LSTM Framework for Deformation Trend Extraction and Forecast

2.1. Description of DLM-LSTM Framework

By comparing the characteristics of GPS latitude time series and field observations, the land deformation shows a stepwise character, which means the average daily latitude changed over adjacent days, but the latitudes within each day remain stationary. Long-term measurement implies that land displacement contains periodic components, which triggers fluctuations that mask the monotonic trend. Based on these observations, we propose a framework to filter out these disturbances and obtain monotonic deformation trends. Then the land deformation can be predicted. Figure 1 shows the flow chart of our DLM-LSTM framework. Firstly, the latitudes received every day are modeled as samples from Gaussian distributions. Mean level and standard deviation are regarded as the input variables to state space functions of the dynamic linear model. Then, we use a recursive Kalman filter to get the optimal states in the DLM model equations. Finally, the land deformation trend is extracted by a recursive Kalman smoother. We apply a simulation smoother in the last step of DLM to describe the full joint distribution of all states and the characteristic of extracted trends. After DLM, the deformation trend is transferred to LSTM networks for training. When network parameters are estimated, LSTM networks are able to predict the displacement using prior trend knowledge. In the end, extracted and predicted deformation trends are compared with observed ground truth displacement to evaluate the performance of our framework.

2.2. Dynamic Linear Model

The dynamic linear model was first defined and introduced in analyzing trends in stratosphere ozone time series [25]. From the definition, “dynamic” was explained as the evolution of the regression coefficient over time and “linear” means that the state space equations are all linear. From our understanding, DLM is the general form of Kalman filter functions. By adding additional operators to the observation matrix and the evolution matrix, simple Kalman Filters can denoise complex time-varying time series to obtain validated trends. In the following subsections, we will introduce the state space function of DLM designed for low-cost GPS time series and a Kalman recursive method for parameter estimation. The land deformation trend will be extracted after the parameters are fully estimated.

2.2.1. State Space Functions

According to [27], the general Kalman state space model is written as where are the observations at time , . The unobserved states, in our case the deformation trends, are defined as vector , which evolves in time with noise according to evolution matrix in (2). The relationship between observed and unobserved states is described in observation (1), where the result of transformed by observation matrix combining noise leads to . Both noises and lie in the Gaussian distribution with covariance and .

The DLM model first expands vectors in (1) and (2) in the following form: where is the mean level of latitude each day and is the change in the level from day to day . Correspondingly, is decomposed to two components: the Gaussian noise distribution with standard deviation and . This form of DLM model is able to detect smooth variation in the mean level and infer changes that happened in the land surface, but it does not consider periodic fluctuations.

To detect periodic fluctuations, the DLM models these fluctuations using harmonic functions. The corresponding operators in the observation matrix is defined as where means the harmonic detected in daily level series within 30 days. We adjust number to reach the best performance during the experiment.

Because paper [22] reveals an autoregressive characteristic in the time series of slowly moving surface displacement, we use a first-order autoregressive model for residual components, , with . In DLM form, we have

To extract land deformation trends from daily latitude time series, the DLM model must combine all the operators above, which leads to a larger-scale evolution and observation matrix as shown below: where and are the 5th and 10th harmonies of 30-day circulation. These two harmonies are found to outperform the denoising abilities of other combinations in periodic variations. Other operators remain the same as above. Respectively, we have a state vector as shown below:

These DLM state space equations contain unknown parameters, including covariance matrix , AR coefficient , and so on In the following step, these parameters will be fully estimated to be able to apply DLM in real situations.

2.2.2. Model Parameter Estimation

In DLM, the unknowns are divided into two major categories: model state and auxiliary parameters defined in state space equations. In this paper, we select the Kalman formula method proposed in paper [25] to obtain the unknown trends and parameters. The essence of the Kalman formula is the recursive prediction and comparison. By continuously updating the covariance matrix in the state space equations, the Kalman formula constructs the connections between measurement and hidden real states. Firstly, we perform Kalman filter forward recursion to predict states , and we assume the conditional probability of predicted states follow normal distributions as below: where represents all the unknown parameters in vector form. Based on this assumption, the daily latitude mean level and standard deviation can be applied as prior input knowledge. Then, the Kalman filter is used as shown below:

In (9), we first calculate the residual between the predicted value and our prediction. Then, we obtain the Kalman gain in 12 via prior and posterior covariance calculated in 10 and 11. With the Kalman gain, we are able to get the next state’s prior mean and covariance of . Then, we apply the Kalman smoother backward to get the smoothed land deformation trend where is the covariance of the observation noise. is the smoothed latitude trend, and is the smoothed state covariance.

2.3. LSTM Recurrent Networks

Time series prediction is about taking some previous input terms, putting them through some hidden units, and predicting the next term. In the past, computer scientists have favored the hidden Markov model (HMM). As the sequence grows larger and larger, the bits of information HMM needs to remember increase largely. Therefore, recurrent neural networks (RNN) have been developed to substitute HMM by storing more prior knowledge in a distributed way. With several different inputs remembered at once, RNN is able to predict more complicated dynamics.

In our paper, we developed an LSTM RNN for land deformation trend prediction. LSTM, so-called long-short time memory units, is a special type of RNN neural. As Figure 2 shows, in each LSTM unit, three gates are designed. The gates serve to help LSTM interact between the memory cell itself and its environment. The input gate can allow an incoming signal to update the state of the memory cell or block it. At the other end, the output gate can allow the state of the memory cell to have an effect on other neurons or prevent it. Finally, the forget gate can modulate the memory cell’s self-recurrent connection, allowing the cell to remember or forget its previous state, as needed. The whole network contains four LSTM memory cells. After the time sequence passes these units, the outputs will be pooled by mean pooling units to get the predicted value.

3. Case Studies

With the aim at illustrating the effectiveness of the proposed DLM-LSTM framework, one land surface deformation case and one static case are studied in the following subsections.

3.1. Case 1: Fushun Western Open-Pit Coal Mine
3.1.1. Dataset

We deployed a wireless sensor network (WSN) on the south slope near the fault zone. The yellow mark in Figure 3 indicates the position of our WSN, which is 41°4954.36 north latitude and 123°5200.31 east longitude. An inexpensive system including hardware and software was designed for monitoring the southern slope surface displacement of the FWOCM. The system consists of two parts, a WSN and a remote server. The sensor nodes of the WSN transmit the data to the remote server as soon as the sensing data are received. The received data in the server are stored, processed, and published on the web. No processing is running on the sensor nodes due to the lack of electricity in the field.

The WSN with a star topology was composed of two kinds of nodes, sensor nodes and a coordinator. All the sensor nodes send their data to the coordinator, and the coordinator forwards the data to the server through a connected General Packet Radio Service (GPRS) module. Both the sensor node and the coordinator were developed using the CC2530 ZigBee development board. For sensor nodes, the sensor was connected to the board through peripherals, such as GPIO or UART. The sensors included a 3-dimensional digital compass (DCM308), GPS (Ublox NEO-6M), and low-cost inertial measurement unit (IMU), soil humidity, soil temperature, atmospheric pressure, atmospheric temperature, atmospheric humidity, and a rain gauge. The communication protocol for data transmission between multiple sensor nodes and coordinators was developed based on Z-stack (TI’s API for low-cost ZigBee communication). The protocol adapts to the solar-powered network. The sensor node can find the coordinator and establish a communication link automatically. Figure 3 shows the position of the WSN. The WSN has been working well since December 1, 2014.

Figure 4 shows the star topology on the left and the layout of three GPS + IMU sensor nodes on the right. The three nodes were located near N41°50, E123°52. GPS sensors were positioned in a line on the north side, while another one was placed on the south side. The horizontal position accuracy of a stand-alone GPS receiver (L1 frequency, C/A code) is 2.5 m. As the sensor data accumulate, long-term land motion calculations based on current and historical data become practical and feasible. In our paper, we selected the north-side GPS sensor as our data source and collected latitude data from February 1 to March 2, 2016. The latitudes are numbered according to the interval between the date of the latitude and February 1, which means the latitudes on March 1 are numbered as day 30.

Figure 5 shows the latitude samples collected each day. As shown in the figure, there is a significant discontinuity in the latitude time series due to the solar power. The latitudes from February 1 to February 6 are all missing, while latitudes from February 6 to February 20 are lacking in quantity. Under the assumption that the latitude changes daily, the time series were composed of daily mean latitude. Our DLM-LSTM framework will be proved effective to extract land deformation trends in submeter accuracy from the time series.

3.1.2. Results

We first calculated the mean level and standard deviation each day to produce a mean latitude level time series. Then, we chose time series from day 6 to day 31 as the input to our DLM model. The mean and variance of latitude within a day are used as input of DLM, which tackles the discontinuous latitude measurements in a day. For days that measurements are totally unavailable, we use linear interpolation of near days as input. Figure 6 shows the result of the DLM model, seasonal-trend decomposition procedure based on Loess (STL), and traditional Kalman filter in land deformation trend extraction. As shown in Figure 6, the DLM (full line) extracts a linearized trend with lesser oscillation than STL, while overfitting may occur in the Kalman filter. In addition, we constructed an estimated trend from the daily mean latitude points, and the trend of latitudes falls in the 95% confidential interval. The 95% probability envelope contains all the mean latitude-level points in the figure, which shows that DLM has a good understanding of the changes in latitude by deformation. Figure 7 shows the result of diagnostic analysis on the residuals by plotting an estimated autocorrelation function and normal probability. In the autocorrelation function (ACF) line, a dashed horizontal line represents the approximative region where the coefficient does not significantly deviate from zero. Figures 6 and 7 validate the extracted trend via estimating autocorrelation functions of DLM residuals. After latitude degrees are transferred to submeters, we are able to compare the extracted displacement with ground truth, which is the professional measurement validated by local authorities.

The difference between the ground truth and the trends extracted by DLM, STL, and Kalman filter, respectively, is given in Table 1 by RMSE, MAPE, and . RMSE and are widely used in statistics. MAPE is the most common measure of forecast error. MAPE functions best when there are no extremes to the data (including zeros). According to the table, the computation results show consistent results with Figure 6 that DLM has better performance in characterizing land deformation trends, especially with lower error rates. Figure 8 shows the comparison of the extracted deformation trend with the ground truth. The ground truth comes from daily professional manual measurement of the displacement in the same period at the south slope of FWOCM. From Figure 8, we can see that our extracted trend shows a consistency with the manual measurement that was recognized by local authorities.

The results in Table 1 clearly indicate that the extraction performance of DLM is encouraging. The error between extracted trend and ground truth remains at the submeter level. As our GPS latitude sensors have a 2.5 m horizontal navigation error, the result of DLM is a big step in accuracy improvement. Figure 9 shows the trend results extracted by the DLM model from April 2 to 30 2016. From Figures 69, DLM is proved to be effective for trend extraction from noisy and discontinuous time series of low-cost GPS sensors.

Because DLM’s performance has been validated, the trend can be used as input for our four-unit LSTM networks. We selected the trend from day 1 to day 26 as our training sample and tried to predict the following 13 days’ trend by the LSTM network. The start point of the break period was predicted by LSTM, and then this prediction was used as input until the validated input was available. The LSTM network parameters are shown in Table 2. The length of 1D mean pooling is 4 according to Figure 2. We use stochastic gradient descent with mini batches to optimize the cost function of the LSTM network. Mini-batch gradient descent is a variation of the gradient descent algorithm that splits the training dataset into small batches that are used to calculate model loss and update model coefficients. In our paper, we set the three training samples in a batch. In addition, we use the Kalman filter to get the dynamic function of the previous 26 days and to predict the next 13 days by this function. Figure 10 shows the prediction result of LSTM and the Kalman filter compared to ground truth measurement. Our LSTM network predicts 13 days’ deformation trend with 7.7194 mm RMSE, as shown in Table 3. Compared to Kalman filter’s 36.7102 mm, this proves that the LSTM network has an accuracy improvement effect in trend prediction. Future works may focus on improving the LSTM network’s performance.

3.2. Case 2: Baoxie Sensor Web Experimental Field
3.2.1. Dataset

To prove that our DLM has a locational invariant propriety, we need to find another place to conduct the validation experiment. The Baoxie sensor web experiment field, located in the east of the central Wuhan area, is the main Sensor Web research and experiment center for Wuhan University’s State Key Laboratory for Information Engineering in Surveying, Mapping and Remote Sensing (LIESMARS). It contains four stations: Baoxie Landslide Monitoring Station, Baoxie Meteorological Experimental Station, Baoxie Edaphic and Meteorological Monitoring Station, and Baoxie Soil Temperature and Moisture Monitoring Station. These stations, in total, are about 800 m2 with more than 70 sensors. All of the sensors are powered by photovoltaic solar energy. Previously, LIESMARS successfully built a cyber-physical geographical information paradigm for in situ land surface monitoring [28], which provided a solid basis for our experiment. As Figure 11 implies, we selected a GPS sensor in a wireless sensor node deployed since 2013 as our experimental data source. The position of this GPS sensor is N30°280, E114°310. The horizontal position accuracy of the GPS receiver (L1 frequency, C/A code) is the same as the GPS receiver at FWOCM, which is 2.5 m.

We collected 14 days of latitude data from December 1 to December 14. Figure 12 shows the latitude samples received each day. It is clear that the GPS sensor in Baoxie is more stable than FWOCM, which has larger deviation in sample numbers each day. However, the Baoxie sensor suffers discontinuity, as does the FWOCM. Also, in situ station measurement shows that the location N30°28, E114°31 remains static with no submeter-level deformation in the same period, which provides the ground truth.

3.2.2. Results

Figures 13 and 14 show the extracted trend and residual analysis of GPS sensor from the Baoxie sensor web experimental field. From Figures 13 and 14, we can see that an obvious static trend is obtained by DLM. For further statistical tests, because the ground truth is recognized as zero, only RMSE can be calculated in Baoxie. The RMSE shown in Table 4 is 2.3542 mm, which obviously proves DLM’s capability in extracting a trend from low-cost GPS sensors.

4. Conclusions

A DLM-LSTM framework is proposed and applied for north-south deformation trend analysis from low-cost GPS time series. A dynamic linear model is introduced to model trend, noise, periodic variation, and autoregressive components hidden in a latitude sequence. Then, the land deformation trend is extracted using a recursive Kalman filter and smoother. LSTM recurrent neural networks are introduced to predict the next states in deformation trend. By analyzing two land deformation cases, the performance of the proposed method is evaluated quantitatively. Calculated displacement in the extracted trend from DLM is excellently consistent with manual measurement from the latitude series of FWOCM. The submeter-level accuracy can be obtained in a north-south deformation trend, which is difficult to obtain by satellite remote sensing approaches. Meanwhile, submeter-level accuracy is also achieved in trend prediction by LSTM networks. The results of this work indicate that the DLM-LSTM framework is a powerful and accurate end-to-end method to analyze land surface deformation trends.

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.

Acknowledgments

This work was supported by the National Key Research and Development Program of China under Grant no. 2016YFB0502601. It was also supported by the Fundamental Research Funds for the Central Universities under Grant no. 2042017kf0211.