Research on Leakage Location of Spacecraft in Orbit Based on Frequency Weighting Matrix Beamforming Algorithm by Lamb Waves

Featured Application: The Frequency Weighting Matrix Beamforming Algorithm proposed in this paper has the potential to locate leakages in large structures such as spacecraft, and the algorithm can adapt the characteristics of leaks. Abstract: Clashes between space debris and spacecraft in orbit may cause air leakages, which pose a substantial danger to the crew and the spacecraft. Lamb wave dispersion in spacecraft structures and the randomness of leak holes are the difficulties in leak location. To solve these problems, a frequency weighting matrix beamforming algorithm is proposed in this paper. The elastic Lamb waves that are caused by leakages are acquired by an ‘L’ shaped sensor array consisting of eight acoustic emission sensors. The angle of a leak can be obtained through the superposition of different time delays, and the intersection of two angles can be used to find the location of the leak. Traditional beamforming is improved by matching the wave speeds in different frequency bands and weightings according to the energy distribution. Narrowband filtering is used to delay overlay different signal speeds with different frequency bands via a dispersion curve. The weighting method is used to compensate the frequency band response of different leak holes. The detailed location algorithm process is introduced and verified by experiments. For 1.5 and 2 mm leak holes, location direction accuracies of 1.33° and 1.93° for one sensor array were obtained, respectively.


Introduction
Orbital objects, such as meteoroids and space debris, are among the serious and inevitable threats to spacecraft [1]. Technically, collisions caused by debris larger than 10 cm can be avoided by using databases that are compiled by debris-tracking systems [2]. Meanwhile, strikes by objects that are smaller than 10 cm lack an efficient detecting approach. A collision may damage the bulkhead of the spacecraft, causing internal pressure to leak, which seriously threatens the flight and the safety of astronauts and leads to severe consequences [3]. Thus, the detection and location of leakages would be of significant use for spacecraft in orbit.
For the problem of leak source location, some related findings have been reported. According to different principles, the currently representative methods include the optical imaging method [4], the optical fiber method [5], the helium mass spectrometry method [6], and the acoustic detection method [7]. Compared with other methods, the acoustic detection method can achieve the non-destructive detection and real-time location of leakages in a larger range with a high location accuracy and stability. Moreover, acoustic sensors are easy to integrate with spacecraft structures, thus eliminating the need for complex processes such as laying optical fibers or using noble gases. Therefore, the acoustic detection method is a potential leakage location technology for spacecraft in orbit [8].
For leakage location based on acoustic methods, the key difficulty lies in the continuity of the leakage acoustic signal. For a collision with a time point of occurrence, location can be performed by calculating the time difference of arrival. To solve the above problems, Holland et al. [9,10] carried out research on cross-correlation and two-dimensional Fourier transform techniques to determine air leakages on the International Space Station (ISS). In the study, two 16 × 16 sensor arrays were applied to locate a leak. All possible cross-correlations between the 256 sensor positions in each array were recorded. A large number of sensors and millions of data samples led to a time-and financeconsuming task. Al-Jumaili et al. [11] proposed a T-mapping method. This method obtains timedomain eigenvalues from a prior meshing experiment and uses a clustering algorithm to identify the position of the signal source. This scheme has good adaptability to complex structures, but the initial workload is huge and it has no application potential for huge spacecraft.
In order to reduce the number of sensors and the algorithm complexity of the leak location scheme, the beamforming technique was firstly applied to the method of leakage detections with acoustic emission sensors by McLaskey et al. [12]. The beamforming algorithm has been primarily used in radar, sonar and exploratory seismology [13]. However, the specimen tested in their experiment was a steel-reinforced concrete bridge ramp instead of a thin plate. The work was based on Rayleigh wave theory, which is not applicable to the thin plates used in spacecraft. Tian et al. [14] tested a thin plate that was equipped with near-field beamforming analysis, which was suitable for Lamb wave theory. Burst acoustic emission signals that were created by breaking a mechanical pencil lead on the surface of the specimens were applied in their tests, and the propagation speeds of the signals were obtained by time difference of arrival (TDOA) techniques. This method is not ideal for a leak-generated continuous signal. However, the difficulties of using beamforming techniques for acoustic emission signal processing lie in the confirmation of wave speeds due to their dispersion. Zhang et al. [15] considered the effect of dispersion, but their scheme did not consider the diversity of the frequency bands that are generated by leakage in practical applications. Thus, it can be seen that the difficulties in the leakage location of thin plates such as those used spacecraft are the effects of dispersion on wave speed and the band adaptability of practical applications.
In order to solve the above problems and to improve the accuracy of the beamforming algorithm for leak location, a frequency weighting matrix beamforming (FWMB) algorithm is proposed in this paper. In this method, narrowband filtering is used to match the dispersion curve to achieve the matrix superposition of wave speeds at different frequencies. An energy evaluation on the frequency band of the signal is performed to determine the selected frequency band and its weighting coefficient matrix, and the weighting method is used to compensate for the frequency band response that is caused by different leaks to achieve a frequency band adaptation. The FWMB algorithm was validated by laboratory tests that used an L-shaped acoustic emission array, and its performance was evaluated. The method in this paper has the potential to be used in spacecraft in orbit due to its light weight and easy integration.
The remainder of this paper is organized as follows. In Section 2, the FWMB process is proposed. In Section 3, the experiment platform in the laboratory is introduced. In Section 4, the conducted experiments are described in order to demonstrate the performance of the FWMB algorithm . The conclusions are presented in Section 5.

Delay-and-Sum Beamforming Algorithm
Beamforming is a mature digital signal processing technology that is commonly used in the fields of radar, sonar, and seismic detection. The basic idea of the beamforming algorithm is to calculate the energy of a wave in a certain hypothetical direction by summing the delays of the signal of each array element to give an estimate of the wave arrival direction. The beamforming method uses a sensor array with a fixed spatial position to measure the spatial sound field, and it processes the signals that are measured by each sensor to obtain detailed sound source information. In the leak detection method that is covered in this paper, the direction of a leak can be obtained by a set of sensors. The final location point can be obtained by the intersection of two sets of sensor arrays.
In beamforming algorithms, the choice of sensor array shape is critical. According to the characteristics of spacecraft structures, a planar array can be used for detection. The shape of a planar array is usually circular, square, triangular, cross, and the like. For spacecraft, the load needs to be as small and light as possible. Research by Cui et al. [16] showed that L-type arrays can obtain the best experimental results with the smallest number of sensors. Therefore, in the following research, an Lshaped array with 8 sensors was used as an example. Moreover, the method proposed in this paper is applicable to other array formats. The demonstration of the beamforming algorithm is shown in  As shown in Figure 1, the sensor array in the figure is labeled n, including #0-7. Among them, the #0 sensor is for reference. When a leak occurs, sound waves propagate through the plate and are acquired by the sensor. f(t,n) is defined as the signal that is acquired by each sensor at the time of t. The relative angle θ between the leak source and the array is defined as the actual direction of the leak source, and θ' represents the hypothetical direction.
When the distance between the sensor array and the leakage source is far greater than the size of the array, the signal waves are considered to propagate along parallel paths. Under this assumption, the signal sampled by the reference sensor #0 is the replica of time-delayed signals that are recorded by other array elements. The time delays between array elements are determined by both by the relative distance from the signal source to array element #n compared to the reference sensor #0 under the direction of arrival and by the wave speeds. Usually, the direction of arrival is unknown; a hypothetical direction of arrival θ' is set. According to the geometric relationship of Figure 1, the time delay of each signal for #0 can be expressed as: where dn represents the relative distance as shown in Figure 1 and v represents the wave speed. dn can be defined as: where d represents the distance between the sensors; this distance is the same as the diameter of the sensor and is 8 mm.
If θ' is consistent with the true angle θ, the signal will be concentrated through the superposition of time delay. Therefore, the signals of each sensor are delayed and superimposed to obtain the following delayed sum signal g(t,θ'): The energy of the signal can be obtained by squaring and integrating the superimposed signal in the time domain. By scanning calculations for different hypothetical angles, an energy function B (θ') related to the angle can be obtained: After obtaining the angle of the peak, it can be used as the result. The location of the leak source can be obtained by the intersection of the two arrays' straight lines.

Dispersion of Lamb Wave
Considering that the bulkhead structure of an orbiting spacecraft is a thin metal plate, due to its large curvature, it is approximated as a metal flat plate in this paper to simplify the research. In thin metal plates, ultrasonic waves propagate in the form of Lamb waves. Lamb waves have less attenuation than body waves, which is beneficial to the location of leaks. However, the propagation velocity of Lamb wave varies with wave frequency, which is called the dispersion characteristic of Lamb waves. Frequency dispersion complicates the characteristics of ultrasonic waves, and it is difficult to directly determine the speed of waves to complete the beamforming algorithm. As a consequence, studying the dispersion characteristics of Lamb waves is critical. There are two modes of waves-symmetric and anti-symmetric-that independently propagate in plates. By equipping a numerical method, equations known as the Rayleigh-Lamb frequency relations for both symmetric and anti-symmetric waves expressed as Equations (5) and (6) can be solved [17].
Anti-symmetric mode: where k represents the wave number, ω is the angular frequency, and d is the thickness of the plate. When the thickness of a plate and the frequency of the signals are given, a set of numerical solutions, denoted by ka and ks, are acquired by solving the Rayleigh-Lamb equations. ka and ks represent the wave numbers of anti-symmetric and symmetric waves, respectively. This phenomenon proves that multiple Lamb wave modes exist in the plate, and each mode has its own phase velocities and group velocities. The number of the numerical solutions increases with the frequency of the signals and the thickness of the plates, which means that the number of wave modes increases. Usually, symmetric modes are represented by S0, S1, …, Sn, and anti-symmetric modes are represented by A0, A1, …, An. In each mode, the corresponding phase velocity also varies with the frequency, and the phenomenon is called the frequency dispersion. Dispersion curves consist of a set of curves can be drawn to depict these characteristics. From the dispersion curves, the theoretical values of phase speeds at any particular frequency-thickness product of any specific mode can be obtained.

Frequency Weighting Matrix Beamforming Algorithm
The traditional beamforming algorithm and the dispersion phenomenon in the plate are introduced above. From the analysis, it can be understood that the traditional beamforming method has the following two difficulties for the leakage location of thin plates such as spacecraft shells: (1) Lamb wave velocity is difficult to determine, and there is an error in the delay matrix.
(2) For practical applications, it is difficult to determine the characteristics of frequency bands by different leaks.
The FWMB algorithm that is proposed in this paper solves the above problems through the following two angles: (1) Narrowband filtering is used to delay different overlaid signal speeds with different frequency bands via the dispersion curve.
(2) A weighting method is used to compensate the frequency band response of different leak holes in real applications. The weighting matrix is determined by the different energy ratios of different frequency bands to adapt to the frequency band changes.
A flowchart comparison between the traditional beamforming [18] and the FWMB algorithm is shown in Figure 2. In the figure, the main differences between the two algorithms are marked by a red border. The following describes the FWMB algorithm based on the traditional beamforming algorithm with some improvements in detail. First, the frequency band of the signal needs to be determined. The signal from sensor #0 is processed by the fast Fourier transform and meshed at 10 kHz intervals (for example, 100-110 kHz). Since the number of points in each frequency band is consistent, all the values in each frequency band are summed and normalized as the energy value of the frequency band. The frequency bands that reach more than 50% of the peak value are reserved, and the bands are defined by the internal frequency as f1, f2, ..., fk (for example, 100-110 kHz is 105 kHz). The total number of the reserved frequency bands is k. The weight of the peak frequency band is 1, and the corresponding weights determine the frequency weighting matrix mf1, mf2, ..., mfk of the different frequency bands.
The signals of sensors #0-7 are narrow-band filtered on the above frequency bands to obtain a processed signal matrix. The processed signal is represented as fk(t,n). For leakage excitation in the thin plate, a large out-of-plane displacement is generated [19]. For the out-of-plane motion, compared to the A0 mode, the energy of the S0 mode in the Lamb wave is relatively small [17]. Because the beamforming method mainly focuses on energy, not the arrival time in the TDOA, the wave velocity of the A0 mode is selected as the wave velocity for each bands. The wave velocity matrixes sf1, sf2, ..., sfk at the above frequencies are determined by the dispersion curve. Therefore, the signals of each sensor are delayed and superimposed to obtain the following delayed sum signal h(t,θ') (corresponds to g(t,θ') in Section 2.1.): Through the above changes, the FWMB algorithm can improve the traditional beamforming algorithm's ability to locate thin plate leaks and its adaptability to different leaks.

Experimental Setup
The performance of the proposed method was evaluated experimentally. The experimental setup is first introduced in this section. The schematic diagram and photo of the experimental setup are shown in Figures 3 and 4, respectively.  A 5A06 aluminum alloy plate, which is consistent with the material of spacecraft bulkheads, was chosen. The test plate had a diameter of 100 × 100 cm and a thickness of 2 mm. A 1 mm diameter hole was drilled in the center of the board, and this was used to simulate the damage on the spacecraft. Grids with dimensions of 5 × 5 cm were drawn onto the plate to indicate the position of the sensor array.
The vacuum pump was used to provide a stable leakage pressure differential to simulate actual leakage conditions. After putting the vacuum-to-plate adapter under the leak hole of the plate, the pumped was turned on until the adapter was attached tightly to the plate. After adjusting the vacuum pump control valve, the subsequent experimental research was able to be started when the measured value of the vacuum degree was less than 10,000 Pa.
The signal acquisition part of the experiment setup was made up by an acoustic emission sensor array, pre-amplifiers, and an acoustic emission instrument. The sensor used in the experiment was a Nano30 sensor (Physical Acoustics Corp., New Jersey, USA) that was 8 mm in diameter and 7.5 mm in height. The Nano30 sensor has a stable frequency response from 100 to 750 KHz. Eight sensors, numbered from 0 to 7, were fixed by special fixtures as 'L,' as shown in Figure 5. The piezoelectric acoustic emission sensors were mounted perpendicular to the plate. The elastic Lamb waves that were caused by the leakage were acquired by the sensor array, and the out-of-plane displacements were mainly measured [20]. The acoustic emission instrument was a DS2-16A instrument (Soft Island Corp., Beijing, China), which can realize multi-channel synchronous acquisition with a 3 MHz sampling rate. The sensors were connected to the pre-amplifiers (Soft Island Corp., Beijing, China). The signals were amplified by 40 dB through the preamplifiers and collected by the acoustic emission instrument.

Experimental Process
A hole located at (0, 0) was drilled into the plate, and it is marked as red circular in Figure 6. Once the set-up was equipped with the air exhaust system, the air evacuated from the hole to the vacuum-to-plate adapter could simulate the leakage. Leak-generated ultrasounds spreading through the plate as guided Lamb waves were acquired by the sensors. The sensor array was placed at position A, of which the actual azimuthal direction of the leak was 9°, and the signal sampling was started as depicted in Section 3.1. The collection of simultaneous leak signals of 8 sensors was repeated 5 times; each data series lasted 0.5 s and was split into 4 segments. Thus, a total of up to 20 segments of data at each position were gained at position A. A series of experiments was done with our array at 9 positions. The positions are marked with green marks in Figure 6. The corresponding azimuthal direction of the leak at positions B, C, D, E, F, G, H and I was 13°, 20°, 24°, 31°, 45°, 57°, 73°, and 79°, respectively. The distance from the sensor array to the leak source was maintained at R = 30 cm during the experiment. Thus, 180 segments of data were acquired in total for one leakage hole. In order to test the leakage of different holes to verify the adaptability of the method in this paper, two plates with 1.5 and 2.0 mm holes were used. From the FWMB theory, the location of signal source can be determined by applying two sensor arrays. Thus, in the rest part of the paper, the orientation accuracy of the leak source is discussed instead of the positioning accuracy.

Results and Discussions
In this section, the above experimental data are processed and described. First, by taking the experimental data as an example, the calculation method of the related matrix is demonstrated to introduce the method proposed in this paper. Subsequently, the results of the experiment are counted to evaluate the performance of the FWMB algorithm.

The Frequency Weighting Factors with Different Holes
One advantage of the FWMB algorithm is that it can be adaptively optimized for different leakage holes. In practical applications, the causes of leaks are abundant, so the shapes and sizes of leaks are various. For different leaks, the spectral characteristics may vary greatly. In order to verify this problem, 1.5 and 2.0 mm standard circular leaks were selected for testing.
Because signals below 100 kHz have strong background noises and signals above 500 kHz are very weak, signals between 100 and 500 kHz were chosen for signal processing. The time domain plot, frequency spectrum, and frequency weighting factors for the signals that were sampled at position A of sensor #0 are illustrated in Figure 7. The amplitudes of the acoustic emission signals that were generated by the two leaks in the time domain were close, and from the waveform point of view, they were both difficult to obtain effective information from and were confusing. After performing FFT on the signal, the spectrum energy of different frequency bands was integrated at 10 kHz intervals and normalized to obtain its weighting coefficient matrix. In the figure, the center frequency of the band is used as the label (for example, 100-110 kHz is represented as 105 kHz). It can be seen from the results that the acoustic emission signal that was generated by the 1.5 mm leak was in a higher frequency band. This is consistent with the physical understanding that smaller gaps result in sharper and higher frequency sounds.
In order to further optimize the algorithm, a band energy threshold was set. The bands with energy smaller than 50% of the peak value were discarded in subsequent processing. It can be seen from the figure that the 1.5 mm leak reserved nine frequency bands of 100-190 kHz, while the 2.0 mm leak reserved eight frequency bands of 100-180 kHz. For the above two leaks, although the frequency bands were similar, the weighting coefficients were significantly different. For holes with more specific shapes, larger differences will occasionally be encountered. The above results illustrate the adaptive ability of the method in this paper to different leaks.

The Speed Matrix
In the previous section, the frequency weighting matrix was obtained. In the FWMB algorithm, another key parameter is the wave velocity matrix. Its role is to accurately match the wave speed in different frequency bands to solve the wave speed error that is caused by the dispersion effect in the flat plate. From the perspective of beamforming theory, wave speed is a critical parameter. An inaccurate wave speed results in the wrong delay, which leads to the wrong location result.
Due to the dispersion characteristics of Lamb waves, signals with different frequency components have different phase velocities. It was assumed herein that the phase velocities were approximately the same in a narrow band of 10 kHz. The numerical solution of the Rayleigh-Lamb equation of the test plate is shown in Figure 8.  The data were processed into eight or nine frequency bands by narrow band-pass filtering: 100-110, 110-120, ..., 180-190 kHz. The corresponding phase velocities of different frequency components in the dispersion curve were recorded as a velocity matrix. For the low frequency thickness product, between 0.3 and 0.6 MHz • mm, there were only two wave modes (the A0 and S0 modes). Since the A0 mode is dominant in the energy of the wave [17], the A0 mode signal is discussed in the rest of this article. Table 1 lists the corresponding phase speeds cp and the weighting factors in different frequency bands. Once the phase speeds and the geometry of a sensor array are given, the time delay matrix can be obtained via Equations (1) and (2). With the process in Figure 2, the FWMB algorithm can be produced. When the power gets its maximum, the corresponding assumed angle is the best estimate of the actual direction of arrival.

Location Results
This section gives a preliminary introduction and discussion of the location results. The 20° leakage at sensor array position C for 1.5 mm is taken as an example. Figure 9 shows the results of the traditional beamforming algorithm and the FWMB algorithm proposed in this paper.

Normalized Energy
For the traditional beamforming, a separate frequency band is usually used to approximate the wave speed information. In order to have more contrast with the FWMB algorithm, a narrow band similar to that in the FWMB algorithm was selected. Figure 9a shows the results in the case of the 120-130 kHz filtering. It can be seen that a maximum value could be obtained near 20°, and a good result was obtained. However, the obtained result was 19°, which, although very close to 20°, still had some errors. In Figure 9b, another frequency band is used as an example. In this example, two peaks appeared, although there was also a significant peak at 20°. However, a higher peak appeared at around 70° and misjudgment occurred. The reason for the above phenomenon is probably that the traditional beamforming only used one single narrow frequency band, and the energy carried by this frequency band was limited and not necessarily the largest. Therefore, there is a certain possibility that the superposition energy obtained from multiple angles was approximated, thereby obtaining multiple peaks. At the same time, the reflection was also a main interference factor for location. Though the energy was reduced for reflection, it was still a source of interference for positioning. The FWMB method can remove the effect of false peaks by superimposing more frequency band results. The above results show that the traditional beamforming algorithm is prone to errors because it uses less frequency band information and is less robust.
To overcome the above problems, under the conditions of this paper, 64 sets of data overlay were used in the FWMB algorithm to improve the eight sets of data overlay in the traditional beamforming. Therefore, it can be seen in Figure 9c that the result graph of the FWMB algorithm has more fluctuations compared to traditional beamforming. Though the data are not smooth, more accurate results were obtained due to the weighted superposition of the multiple sets of frequency band data. In view of the above preliminary conclusions, we show the data of several other points as further displays.
The results for the leakages in positions I and B are shown in Figure 10. The graphics use a circular coordinate system for easier display and understanding. The red curve and angle value represent the FWMB result, and the blue curve and angle value represent the traditional beamforming result. The frequency band in the label indicates that it is used in traditional beamforming. For the FWMB algorithm, the frequency band is adaptive and does not require manual selection. It can be further seen from the results that the FWMB algorithm optimized the results from two perspectives for traditional beamforming. First, due to the matrix-type multi-data superposition, the discrimination of the results was further enhanced and small errors were resolved, such as within 2° in Figure 10a,d,e,f. In addition, through the selection of the energy band and the calculation of the weighting coefficient, the situations where multiple peaks appeared and misjudgment occurred, such as in Figure 10b,c, were improved.
It can also be seen that the FWMB algorithm results also had obvious fluctuations. For example, there is a multiple peak at about 13° in Figure 10a-c. This shows that the FWMB algorithm may have misjudged the location, especially when the leakage was near the edge, which was mainly related to the reflection of the boundary. However, due to the large number of FWMB overlays, this problem was compensated for by averaging, and correct results were almost obtained.

Results Statistics and Evaluation
This section summarizes the results and evaluates the performance of the FWMB algorithm. The statistics of the leakage angle location results for two different holes are shown in Figure 11. The upper and lower ends of the error bar in the figure indicate the maximum and minimum values in the measurement data. As can be seen from the figure, the FWMB algorithm was able to get a close location result for different angles. Judging from the average error, reading errors of 1.33° and 1.93° were obtained in the holes of 1.5 and 2.0 mm, respectively, and the comprehensive error was 1.63°. The above results prove that the FWMB method proposed in this paper is adaptive to different leaks. The above advantages are obtained by weighting the energy ratios of different frequency bands in the FWMB method. Due to the use of narrow frequency band, the dispersion problem of Lamb waves is solved, and the wave velocity parameter in the algorithm is more accurate.

Discussion and Prospection
The above experiments confirm the improvement of traditional beamforming by the FWMB algorithm. However, the FWMB algorithm still has some issues that need to be addressed in future research. The leak hole that was used in this paper was a circular hole. In subsequent research, irregular holes that are caused by real collisions will be tried to more comprehensively evaluate the effectiveness of the method. • This paper evaluated the positioning performance at different angles, and further experimental and theoretical analyses of the effective distance of the method and the key influencing factors will be made in subsequent research centers.

•
Boundary reflection is also an important factor that causes mispositioning, and no relevant experiments were done in this paper. For large spacecraft structures, the attenuation of waves reduces the effects of reflections. However, the above issues are still worthy of discussion and research. In subsequent experiments, a non-reflective case by adding sound-absorbing cement can be tried.

•
In an actual application environment, the structure may have complex structures such as curved surfaces or stiffeners. The effectiveness or compensation method of the FWMB algorithm is also worthy of attention and research.

•
The ability to locate leaks at different locations needs to be further evaluated, especially for leaks that are close to the boundary. Though the FWMB algorithm's multi-layer overlay scheme can obtain more accurate results, there will still be multi-peak situations, and the algorithm needs to be improved in the follow-up.
In summary, for continuous leakage, the FWMB algorithm provides significant improvements. In subsequent developments, optimization and further testing are needed based on actual applications.

Conclusions
This paper proposes a frequency weighting matrix beamforming algorithm to solve the problem of locating continuous leaking sound sources in plate structures. Compared with the traditional beamforming algorithm, this algorithm uses a narrow frequency band to match the wave speed to solve the influence of the Lamb wave dispersion of the thin plate on the signal positioning. Through the energy-weighted method, the algorithm can adapt to different leaks, which has obvious advantages in application. The above scheme was verified through experiments, and 1.5 and 2 mm leaks were tested. The results showed that the frequency weighting matrix beamforming algorithm that is proposed in this paper can obtain a location accuracy of 1.63°. The above results can be used for the leak detection of on-orbit spacecraft, because the shape of the leak holes that are generated by space debris is inconsistent, so the adaptability of this method can solve the above problems.
Author Contributions: L.Q., X.R. and L.S. established the theoretical method; X.L. and L.W. designed the experiments; G.Y. and T.L. performed the experiment; L.Q. and Y.Z. analyzed the data. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.