Abstract
Transient electromagnetic apparent resistivity imaging technology is one of the more promising methods for external inspection of metallic oil and gas pipelines. Through the research on the transient electromagnetic response and imaging technology of pipelines, it is found that the accuracy and real-time performance of the apparent resistivity calculation are the key to its application. To achieve fast imaging, a three-layer BP neural network is designed with the kernel function of the secondary field as the input and the transient parameter value as the output; the nonlinear equation of transient response is fitted by the neural network to solve the apparent resistivity, and inversion depth is calculated based on smoke ring theory. Aiming at the shortcomings of the traditional BP network, such as slow convergence rate and the ease of falling into local minima, the genetic algorithm is designed to optimize the initial weight and threshold of the network. In the model pipeline experiment, the measured data are brought into the trained GA-BP network, and calculation time is greatly shortened. The obtained sectional image can directly and accurately reflect the pipeline shape. The validity and practicability of the transient electromagnetic apparent resistivity imaging technology based on the GA-BP neural network are verified, which is expected to be a powerful tool for real-time evaluation of pipeline corrosion detection.
1. Introduction
As of the beginning of 2018, the total mileage of China’s long-distance oil and gas pipelines has accumulated to approximately 133,100 kilometers. This increased length of pipelines makes in-service corrosion inspection more and more important for the safe operation of oil and gas pipeline networks. At present, internal inspection devices and technologies based on magnetic flux leakage and ultrasonic waves are mostly adopted, but they are all affected by the pipeline structure, operating environment, and transmission medium, especially the buried pipelines, which are limited by high cost and difficulties in implementation and accurate positioning of defects. The noncontact pipeline external inspection technology with the advantages of trenchless and nonstop transmission will be the main development direction for breaking through the bottleneck of the pipeline testing industry in the future.
The transient electromagnetic method is a kind of controllable artificial-source time-domain electromagnetic testing method based on the electromagnetic induction principle: the primary magnetic field generated by the pulse current is transmitted by the transmitting coil above the pipeline, and the eddy current secondary field signal carrying the electromagnetic characteristics of the metal pipe is received by the receiving coil; apparent resistivity imaging technology is used for signal analysis and inversion interpretation; finally, the position and metal quantity change of the pipeline can be displayed visually to achieve noncontact external inspection of oil and gas pipelines.
Tartaras and Zhdanov first proposed the S-inversion method based on the “floating thin plate model,” which converted the time derivative matrix of the detected vertical magnetic-field component into the apparent conductivity-depth value to obtain a continuous image [1, 2]. Lee et al., drawing on the thought of seismic interpretation, used wave-field conversion to realize quasiwave equation migration imaging [3]. These methods are all approximate interpretations and are suitable for mineral exploration in deep formation. If used for inspection of buried pipelines, the resolution is low and the precision is not sufficient. According to the smoke ring theory proposed by Nabighian [4], the apparent resistivity-depth image of the pipeline can be obtained by time-depth conversion; the key is to accurately determine the all-time apparent resistivity in order to effectively reflect the position and shape of the pipeline. Ward and Hohmann derived the analytical expression of apparent resistivity under the central-loop device [5]. Raiche and Spies proposed an iterative algorithm for calculating the late-time apparent resistivity accurately [6]. Christensen obtained a single-valued all-time apparent resistivity by calculating the vertical magnetic field component [7]. More research is focused on using dichotomy or Newton’s iteration method to solve the apparent resistivity nonlinear equation, which has high accuracy and wide applicability, but they require high computing resources and time and cannot achieve real-time fast imaging when survey lines are long and measuring points are dense.
Therefore, the method of solving the transient electromagnetic apparent resistivity under the central-loop device using the BP neural network optimized by the genetic algorithm is proposed. The BP (backpropagation) neural network is a supervised multilayer feedforward network trained by the error backpropagation algorithm [8]. The self-learning ability of the neural network is used to fit the nonlinear equation, and the transient parameters are calculated according to the kernel function value based on the measured data. Then, apparent resistivity is obtained and fast imaging is realized. Finally, the feasibility and accuracy of the method are verified by pipeline model experiments.
2. Transient Electromagnetic Apparent Resistivity Imaging Technology
In the detection using the transient electromagnetic method, the central-loop device is optimally coupled with the oil and gas pipelines, and the reflected abnormal amplitude and lateral resolution are the largest. Therefore, this paper discusses the inversion imaging technique of transient electromagnetic response using the central-loop model under the step pulse current excitation.
The time-domain expression of the step pulse current excitation function iswhere I is the amplitude of the emission current (A) and t is the transient delay time after the emission current is turned off.
According to the Duhamel integral, the analytical expression of the transient response under the central-loop device is obtained as [4, 5, 9]where Bz is the vertical component of the secondary field magnetic induction intensity, is the partial derivative of Bz to time, μ0 is the homogeneous half-space magnetic permeability (approximates to 4π × 10−7 H/m), a is the radius of the transmitting coil, ρ is the homogeneous half-space apparent resistivity, u is the transient parameter defined by ρ, and is the error function. The computational expression of ρ is obtained from equation (4) as
The transient response expressions (2) and (3) are normalized, and and are defined as kernel functions of Bz and , respectively; then,where and are obtained by Bz and measured by means of a detecting instrument, and the remaining variables are determined by the operating parameters of the device. From equation (5), the key to solving the apparent resistivity ρ is to find the root u of the nonlinear equations (6) and (7).
According to Nabighian’s smoke ring theory [10], the step transient response in homogeneous half-space excited by the surface transmitting coil can be expressed by annular electric current, which diffuses downward and outward with time. Through the time-depth conversion, the time-varying transient electromagnetic response measured on the surface is converted into the curve of resistivity with depth to realize fast imaging interpretation.
Oil and gas pipelines are generally laid horizontally, which can be regarded as one layer of layered medium. The vertical depth d of annular current at a certain moment is
The diffusion velocity of annular current iswhere ρj and ρi are the all-time apparent resistivity corresponding to the sampling times tj and ti of the adjacent time channels.
Every layer’s apparent resistivity ρl can be expressed aswhere tij is the arithmetic mean of two adjacent sampling times.
The corresponding apparent depth Hl is
In summary, the apparent resistivity imaging process is shown in Figure 1.
3. BP Neural Network Optimized by Genetic Algorithm
The premise of realizing transient electromagnetic fast inversion imaging is to accurately calculate the apparent resistivity value of a pipeline, but there is a relatively complicated implicit function relationship between the transient response value of buried pipelines and the apparent resistivity. Conventional numerical methods (such as iterative methods) are complex and inefficient. Nevertheless, the artificial neural network (ANN) is particularly suitable for the processing of uncertain and unstructured data. With the help of previous data, the relationship between the kernel function and the transient parameter can be established directly, which makes fast imaging possible.
The BP neural network includes one input layer, several hidden layers, and one output layer. The learning algorithm is realized by signal forward propagation and error backpropagation: sample signals enter the hidden layer from the input layer after weighted cumulative processing, and hidden layers generate the output signal to the output layer through complex nonlinear transformation; if the result does not match the expected output, the error will be returned to the input layer through the hidden layer; each layer modifies the weight and threshold of every neuron according to the error signal obtained by itself so that the error decreases along the gradient direction; after updating the connection weights and thresholds, the forward propagation process of the signal is repeated; repeated training will continue until the output error meets the desired range or the predetermined learning times are reached. The structure of the BP neural network is shown in Figure 2.
Take a 3-layer BP network with one hidden layer as an example. The input layer has m neurons, and x = (x1, x2, ..., xm) is the input vector of the sample; the hidden layer has l neurons; the output layer has n neurons, and y = (y1, y2, ..., yn) is the output vector.
The BP network belongs to the multilayer feedforward backpropagation network. Each neuron in each layer can only connect with the neurons in the former and the latter layer: when the signals propagate forward, the neuron receives the signals from all the former layer’s neurons and transmits the output to all the neurons in the next layer. No information is transmitted between neurons in the same layer.
For a given sample , the input of the j-th neuron in the hidden layer iswhere is the connection weight of the i-th neuron in the input layer to the j-th neuron in the hidden layer, θj is the threshold of the j-th neuron in the hidden layer, f is the activation function, and xi is the input (i = 1, 2, …, m).
The input of the k-th neuron in the output layer iswhere is the connection weight of the j-th neuron in the hidden layer to the k-th neuron in the output layer and θk is the threshold of the k-th neuron in the output layer.
The output error Ep of the sample iswhere dk is the expected output of the k-th neuron in the output layer and yk is the actual output ().
The conventional BP network training algorithm adopts the gradient descent method or Gauss–Newton method to correct the connection weight and threshold of the hidden layer, input layer, and output layer according to the principle of error gradient descent.
In order to enhance the performance of the BP network, the Levenberg–Marquardt (L-M) algorithm is used to train the network, which is better for solving nonlinear problems.
The L-M algorithm is a combination of the gradient descent method and the Gauss–Newton method, which can achieve seamless movement between the two methods and greatly improve the convergence rate. The idea is that each iteration does not follow a single negative gradient direction but allows errors to be searched along the direction of deterioration and dynamically adjusts the convergence direction of iteration. By adaptively optimizing the weights between the gradient descent method and the Gauss–Newton method, the training can converge effectively, and the convergence rate and generalization ability of the network can be greatly improved. New weights and thresholds are updated as
The correction quantity ∆ iswhere I is the identity matrix, is the error vector, s is the learning step, and μ is the adaptive adjusted damping factor. When the solution is close, μ ⟶ 0, the weight adjustment approaches the Gauss–Newton method so that the algorithm converges quickly near the neighborhood of the solution; when the distance from the solution is long, μ ⟶ + ∞, the weight adjustment approaches the gradient descent method to ensure that the cost function is reduced when the operation is difficult. is the Jacobian matrix composed of partial derivatives of network errors to weights. Take t weights as an example,
After getting the new weights and thresholds according to the correction quantity ∆, the square sum of error Ep is calculated: if Ep decreases, after dividing μ by β, the network output is recalculated using the new weights and thresholds, and the weights and thresholds are continuously updated. If Ep increases, β is multiplied by μ without updating the weights and thresholds, and equations (16) and (17) are used to calculate new weights and thresholds; β takes a real number greater than 1.
For each given sample, the connection weights and thresholds of the network are updated layer by layer based on the error. When all given samples are trained, the network global error E is calculated as
If E satisfies the preset precision or the learning times reach the set value, the network meets the requirements and ends the algorithm; otherwise, the next sample is selected to continue the training.
The BP neural network has good performance in self-learning, self-adaptation, and generalization capabilities, but for nonlinear systems and complex objective functions, the use of the unoptimized gradient descent method will lead to the problem of long training time and local optima. The genetic algorithm (GA), a global optimization search algorithm, is selected to optimize the BP neural network to accelerate convergence and ensure the global optimum of training results.
The GA uses a heuristic global random search method to simulate the process of biological evolution in nature [11]. By means of a series of genetic operations such as selection, crossover, and mutation, the initial population containing the candidate solution set of the problem is gradually evolved to a state containing the approximate optimal solution.
The steps to optimize the BP network using the genetic algorithm (to get the GA-BP network) are as follows:(1)The transient electromagnetic method inspection data are normalized; the network topology, the number of layers, the number of neurons in each layer, etc. are determined; and the network is generated.(2)The initial weights and thresholds are encoded to form the initial population; each individual in the population represents a combination of initial weights and thresholds of the network.(3)The reciprocal of the objective function of the BP network is used as the fitness function to calculate the fitness of each individual, and the selection, crossover, and mutation genetic operators are applied to evolve the current population to generate a new population.(4)The above steps are repeated until the population reaches the maximum evolutional generation, and the genetic operation is ended; the individual with the greatest fitness in the current population is selected as the optimal individual.(5)The optimal individual is decoded, and the obtained weights and thresholds are substituted into the BP neural network as the initial value for training until the network performance meets the desired requirements.
The flowchart of the GA-BP network algorithm is shown in Figure 3.
4. Solution of Apparent Resistivity by GA-BP Neural Network
In view of the excellent performance of the GA-BP neural network in nonlinear curve fitting, according to the above-mentioned research, the neural network with kernel function as the input and transient parameter u as the output based on the nonlinear curve fitting method can be used to calculate the transient parameter u corresponding to the kernel function value. In practical application, the corresponding u is calculated by the trained neural network through the kernel function value based on the measured data so as to avoid specific complex electromagnetic field numerical calculations, and the apparent resistivity is calculated according to equation (5) to achieve the imaging.
4.1. Determination of Topological Structure of BP Neural Network
According to equations (6) and (7), the kernel function of time partial derivative of the secondary field is a dual-valued function, so u corresponding to is not unique; sectional treatment is required when establishing a network, which leads to complicated training and difficulty in interpreting the results. The kernel function of the quadratic vertical magnetic field Bz is a single-valued function, and and the transient parameter u are a one-to-one match, which can better reflect the electrical structure of the pipe-electric section. Therefore, the kernel function is taken as the input to create a BP neural network with u as the output.
Generalization ability of the BP network is determined by the complex degree of the problem, the structure and size of the network, and the number and representative of given samples. According to Kolmogorov’s theorem, a three-layer BP network can achieve the approximation of any nonlinear function. In this paper, we design a 3-layer BP network containing one hidden layer with as the input and u as the output.
There is no general method to determine the number of neurons in the hidden layer. If the number is too small, the network performance will be very poor and the law of samples cannot be reflected; if the number is too large, the “overfitting” phenomenon will appear, the training time will be longer, and it will easily fall into the local optimum solution. The errors of various neural network models with different numbers of hidden neurons are analyzed through simulation research, and the calculation results are shown in Figure 4. It can be seen from Figure 4 that the number of neurons in the hidden layer is not a simple linear relationship with the error. When the number is 13, the error is the smallest.
Based on the above analysis, the final designed BP neural network has a three-layer architecture, as shown in Figure 5. The network has one input layer with single neuron , one hidden layer with 13 neurons, and one output layer with single neuron u.
The BP network belongs to supervised learning, and the training samples are determined by numerical calculation according to equation (6). The target test object is a metal pipe, and the environmental media are mainly soil and air, so the resistivity of the inspection area is less than 104 Ω·m. Combining the acquisition parameters of transmitting and receiving devices, including the radius of coils and sampling interval, the value range of u determined by equation (5) is [10−3, 200]. In this range, u is equal-interval sampled to calculate the corresponding . A total of 20000 groups of data are obtained: 80% of them are taken as the training set and the remaining 20% as the validation set. The differentiable nonlinear sigmoid function is selected as the network activation function. The goal of the square sum of error is set at 0.001, and the maximum number of iterations is 1000.
It is known by the backpropagation algorithm that the larger the learning rate is, the greater the weight change and the faster the convergence speed will be. However, excessive learning rate will cause network oscillation, leading to low performance, and it is easy to cause algorithm instability when exceeding a certain extreme value. The smaller learning rate can ensure the stability and convergence of the system, so the learning rate in this paper is set at 0.05.
4.2. Design of Genetic Algorithm
After the network structure is determined, the weight set and threshold set will be generated randomly, and global optimization is carried out on them by the genetic algorithm. The obtained global optimal or approximate optimal weights and thresholds optimized by the genetic algorithm are used as initial weights and thresholds of the BP algorithm for further training the neural network.
When designing a genetic algorithm, the connection weights and thresholds should be coded first. They are transformed into chromosomes in the population with the representation of the gene string code to construct the search space which the algorithm can deal with. The input layer and output layer of the BP network have 1 neuron, and the hidden layer has 13 neurons, so the number of connection weights is 1 × 13 + 1 × 13 = 26. The hidden layer needs 13 thresholds, and the output layer has 1. Then, the total number of parameters that need to be optimized by the genetic algorithm is 26 + 13 + 1 = 40. Conventional binary coding would result in a too large solution space for training and low search efficiency. In this paper, real coding is used: the gene code of each chromosome is represented by a real number in a certain range, which facilitates the introduction of heuristic information related to the problem domain to increase the search ability of the algorithm.
After population initialization coding is completed, based on the simple and general principle, the fitness of each chromosome is evaluated by evaluation function derived from the reciprocal of the objective function of the neural network, which measures that the degree of individuals in the group can achieve or contribute to finding the optimal solution in the calculation. As one of the control parameters, the population size—the total number of chromosomes—has a significant impact on the efficiency of the algorithm: too small size is not conducive to evolution, while too large size will lead to an increase in computational complexity and an excessively long convergence time. Generally, the population size is between 10 and 200, and 50 is selected in this study.
After all chromosomes are evaluated, the selection operator is used to extract chromosomes from the current population to create an intermediate population. The selection operation uses a combination of optimal individual preservation and roulette: the best individual preservation strategy is adopted first; that is, the individuals with the highest fitness do not participate in the crossover and mutation operations but are directly copied to the next generation to ensure that the optimal solution of a certain generation will not be destroyed. In order to avoid falling into local optima, roulette selection is adopted for other individuals in the current generation; that is, the probability of an individual being selected is proportional to its fitness; the higher the fitness, the greater the probability of being selected, which not only ensures the diversity of individuals in the population but also keeps the fitness function value of individuals close to the optimal solution.
Crossover and mutation operations are carried out in the intermediate population generated by the selection operation to produce the next-generation population representing the new solution set. Crossover operation is the operation of replacing and recombining the partial structure of two parent individuals to generate new individuals, which is the main method to generate new individuals and determines the global search ability of the algorithm. In this paper, the arithmetic crossover method matching with real number coding is used to generate new individuals from the linear combination of two individuals.
Mutation operation is to replace the values at some loci in an individual’s coding string with other alleles at those loci to form a new individual, which is the auxiliary method to generate new individuals but determines the local search ability of the algorithm, and it is possible to repair important genes lost during the crossover process. We adopt nonuniform mutation to make random perturbations to the original gene values, and the results are used as the new gene values after the mutation. Each locus is mutated with the same probability, which is equivalent to a slight change in the whole solution space.
The key to control crossover and mutation operations is crossover probability and mutation probability . A larger can enhance the ability of the algorithm to open up new search areas, but the high-performance excellent schema is more likely to be destroyed. If is too small, more individuals will be copied directly to the next generation, and the search may stall. A larger can increase the diversity of antibody population, but it is possible to make the algorithm tend to pure random search; meanwhile, a smaller is not conducive to the generation of new individuals and inhibition of premature convergence.
and are generally selected within the optimum parameters suggested by predecessors. However, the adaptive genetic algorithm proposed by Srinvivas is used in this paper, which is closely related to fitness, and can ensure the convergence of the algorithm while maintaining the diversity of the population:where fmax is the maximum fitness value in the population, f is the larger fitness value of the two chromosomes participating in the crossover, favg is the average fitness value of the current population, fmin is the minimum fitness value in the population, and f′ is the fitness value of the individual to be mutated.
After the new population is generated by selection, crossover, and mutation operators, the evaluation and creation procedures of all chromosomes are repeated, and the fitness of the population is continuously improved until the stop criteria are met. Then, the chromosomes in the final optimal population are decoded as the initial weights and thresholds of the network; the backpropagation training is carried out using the L-M algorithm until the performance index meets the expected value, and the final optimal value of the BP network is obtained.
After the parameters of the BP neural network and genetic algorithm and samples are determined, the neural network is trained by classical BP and GA-BP, respectively, and the results are compared, as shown in Table 1. By comparison, it can be concluded that the GA-BP algorithm has achieved better results than the classical BP algorithm in terms of convergence speed, elapsed time, and accuracy. The results show that the GA algorithm can shorten the search space at a faster speed and that it is not easy to fall into local minima, while the BP algorithm has the characteristics of high local search efficiency. Combining the two methods can achieve better learning results than conventional algorithms.
Finally, the GA-BP neural network composed of the optimal weights and thresholds between each neuron is used as the tool to solve apparent resistivity: value obtained from the measured data is used as the input to directly map the corresponding transient parameter u, and then the apparent resistivity at each sampling point can be calculated.
5. Experimental Study
Through the model pipeline experiment, the constructed GA-BP network is tested and its effectiveness and practicability are analyzed. In order to verify whether the imaging technology established in this paper can accurately reflect the shape of the pipeline, the model pipeline is composed of three steel pipes of different diameters and wall thicknesses. The specifications are φ275 × 10 mm, φ246 × 13 mm, and φ271 × 5.5 mm, and the total length is 6700 mm. The entire measuring line is 8 m long. Measuring points are arranged along the axis of the pipeline and isometrically divided into 41 points with a distance of 200 mm. Points 5–37 cover the whole length of the model pipeline, while points 1–4 and 38–41 are located on the ground outside the pipe end for comparison and reference. The experimental model pipeline section and the layout of measuring points are shown in Figure 6.
The experimental testing device is the “oil and gas pipelines’ transient electromagnetic method external inspection equipment” developed by the laboratory. It is composed of an antenna device, a transmitter, a receiver, and an upper computer and has the characteristics of an integrated transceiver, remote operation, and so on. The host and working connection mode of the transient electromagnetic method external inspection equipment are shown in Figure 7.
The experimental parameters, device parameters, and working parameters are shown in Table 2.
Before the experiment, the ambient noise level and background field of the pipeline-free site are measured first, and the results show that the experimental environment is basically free from interference and the background field is relatively stable. Then, the pipe is laid flat on the ground and the testing is started; during the inspection, each measuring point should be observed twice at least, and the relative error of the two observations should not exceed 3%; if it exceeds, more observations should be made, and the average of 2-3 sets of data with the least deviation is taken as the response data of the measuring point. The experimental site is shown in Figure 8.
(a)
(b)
In order to ensure the consistency of the detection data and facilitate postprocessing and interpretation, the induced voltage value normalized by the emission current () is recorded in unit of μV/A by the equipment. The multichannel profile constituted by measured data, as shown in Figure 9, is the connection set of induction electromotive force at the same sampling time of all measuring points in one test, which also reflect the variation of measured voltage with time attenuation on the pipe section.
The measured data of the first few time windows reflect the trend of the secondary field changing with the shape of the pipeline, but further quantitative inversion interpretation is still needed.
Before calculating the apparent resistivity, the induced voltage is converted into the transient field value , and the kernel function obtained by equation (6) is used to enter the trained GA-BP neural network to calculate the transient parameter u; then, the apparent resistivity is obtained from equation (5), and the apparent resistivity-depth relationship is calculated according to the smoke ring inversion theory. Because the obtained apparent resistivity of each measuring point is discontinuous, the nearest neighbor interpolation method, which is more effective for filling a uniform data interval area and a nonvalue area, is used between each measuring point. Finally, the smooth pipeline section image is obtained by the Surfer software, as shown in Figure 10.
The distribution and depth of the pipe can be clearly shown from the image. The pipeline layout range starts from measuring point 4, and the apparent resistivity is gradually increased relative to the environmental medium; the abnormal apparent resistivity displayed at measuring points 5–37 reveals different outer diameters of the three sections of pipes; measuring points 17 and 26 correspond to the joint face of pipelines, which is in accordance with the actual model. Starting from measuring point 37, the testing coil is gradually removed from the model pipeline, and the apparent resistivity is decreased. The image perfectly matches with the actual pipe shape. The time required for calculating apparent resistivity by the conventional numerical algorithm and trained GA-BP neural network is shown in Table 3.
The results show that the GA-BP neural network with powerful parallel computing ability is obviously superior to the conventional numerical method in terms of computing speed; the computational time is less than 1 s, and the accuracy is also in line with the requirements. It is proved that the theoretical method and transient electromagnetic method inspection equipment studied can be applied to the detection of pipeline burial depth and metal loss, which can be used as a general survey tool for pipeline corrosion.
6. Conclusion
The key of transient electromagnetic apparent resistivity inversion imaging technology is to calculate the apparent resistivity quickly and accurately. In this paper, a BP neural network with the input of the kernel function of the secondary field intensity and the output of the transient parameters is constructed; the apparent resistivity under the central-loop device is solved based on the fitting of the nonlinear equation; the depth algorithm of smoke ring inversion is used to achieve the real-time imaging. In order to overcome the disadvantages of the conventional BP network, such as slow convergence, the genetic algorithm is used to optimize the initial weight and threshold of the network according to the objective function to improve the network performance. The model pipe experiment and imaging results show that the apparent resistivity inversion technology based on the GA-BP network can accurately describe the shape and distribution position of the pipe, and the trained network can be used as a fast and effective tool to solve the apparent resistivity under the central-loop device, which provides an intuitive and practical means to determine pipeline metal loss through the inversion results.
In the process of experiment, the pipeline model with obvious morphological differences is specially adopted to test the feasibility of the trenchless pipeline transient electromagnetic inspection method established in this paper. In the next step, the prefabricated buried pipeline with pre-prepared defects will be tested in real time, focusing on the imaging effect of minor defects and the influence of soil characteristics on the detection results, and the application range and precision of neural network inversion imaging technology will be studied further.
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 research was supported by the National Key Research and Development Program of China (Grant no. 2016YFC0802100).