A Fast and Accurate LEO Satellite-Based Direct Position Determination Assisted by TDOA Measurements

Shijie Li,Qianyun Zhang,*,Boyu Deng,Biyi Wu,Yue Gao

1 School of Cyber Science and Technology,BeiHang University,Beijing 100191,China

2 Tsinghua University,China,Beijing 100084,China

3Beijing Institute of Technology,Beijing 100081,China

4 School of Computer Science Fudan University,Shanghai 200433,China

Abstract: This paper focuses on the trusted vessel position acquisition using passive localization based on the booming low-earth-orbit (LEO) satellites.As the high signal-to-noise ratio (SNR) reception cannot always be guaranteed at LEO satellites, the recently developed direct position determination(DPD)is adopted.For LEO satellite-based passive localization systems, an efficient DPD is challenging due to the excessive exhaustive search range leading from broad satellite coverage.In order to reduce the computational complexity,we propose a time difference of arrival-assisted DPD (TA-DPD) which minimizes the searching area by the time difference of arrival measurements and their variances.In this way, the size of the searching area is determined by both geometrical constraints and qualities of received signals,and signals with higher SNRs can be positioned more efficiently as their searching areas are generally smaller.Both two-dimensional and three-dimensional passive localization simulations using the proposed TA-DPD are provided to demonstrate its efficiency and validity.The superior accuracy performance of the proposed method, especially at low SNRs conditions, is also verified through the comparison to conventional two-step methods.Providing a larger margin in link budget for satellite-based vessel location acquisition,the TA-DPD can be a competitive candidate for trusted marine location service.

Keywords:direct position determination; LEO satellite constellation;passive localization;TDOA

I.INTRODUCTION

With the rapid development of global trade, worldwide shipping has been significantly increased.For many maritime services, like smart ocean, logistics,and navigation, the position information is crucial to provide customized and efficient location-based services [1].Besides, trusted vessel location acquisition is of great importance to fight against forbidden marine activities like illegal, unreported and unregulated fishing(IUU),hazardous cargo transports, smuggling of goods and humans [2, 3].The automatic identification system (AIS) is currently used to facilitate the tracking and monitoring of vessel location and movement[4].However,due to the lack of security considerations,the position information can be easily modified illegally or deliberately by spoofing[5],tampering[6],false data injections[7],and et al.

Solely relying on the propagation characteristics of wireless signals,the passive localization technique determines the emitter position based on the received signal waveform by one or multiple sensors.Therefore, it is immune to malicious attacks in network and application layers.Two-step passive localization techniques, like angle of arrival (AOA), time difference of arrival (TDOA), and frequency difference of arrival (FDOA), have been well developed over the past decades.Determining emitter locations based on estimated intermediate parameters, the localization accuracy of these approaches deteriorates significantly as the signal-to-noise ratio (SNR) gets lower[8-10].An alternative localization technique-onestep method, also known as direct positioning determination(DPD)-was proposed and studied in recent years [8, 11-13].In one-step positioning, instead of extracting intermediate parameters, all received signals are proceeded simultaneously at a data fusion center for position determination.Avoiding the process of correlating measured parameters with the signal source,the DPD has been demonstrated to exhibit higher accuracy than the two-step methods,especially at low SNR[12,13].

Owing to the broad coverage range and few geographical restrictions, passive positioning using the low earth orbit(LEO)satellites has been of great significance in marine localization.However, although a higher positioning accuracy is able to be achieved[14], the use of DPD in a satellite-based localization system can be computationally prohibitive because of the wide observation area.Besides,as the relationship between the cost function and position of the emitter is highly non-linear, it is quite difficult to provide a closed-form expression for the emitter position.

Recently, various studies have been made to improve the efficiency of the one-step localization process.The most widely used approach employs the coarse-to-fine strategy, which first searches on coarse grids over a large area and then follows a fine grid search around peaks on the coarse grids [15-17].To further improve the searching efficiency, the finegrained grid search can be replaced by the parabolic interpolation as well[18].Another strategy to accelerate the peak searching in DPD is achieved by simplifying the highdimensional exhaustion search into multiple one-dimensional searches.The implementation of this strategy includes the expectation-maximization[19] and alternating projection [20].Moreover, the peak detection of DPD cost functions can also be efficiently determined by intelligent optimization techniques.In[21],the cost function peak is first roughly estimated by the niche particle swarm optimization and clustering, and then precisely determined by the gradient search.However, due to the wide satellite field-of-view(FoV),the footprint of a LEO satellite on Earth’s surface is up to thousands of kilometers,and a coarse-grained grid search with kilometer resolution is computationally prohibitive.Decreasing the coarsegrid resolution helps to reduce the computational burden but might causes positioning failure, as the cost function peak might be completely ignored due to the large grid size.Furthermore,an overly large searching area will also introduce unnecessary errors like false alarms with inherent phase ambiguity [22].Therefore,applying DPD to efficient vessel localization using LEO satellites is still quite challenging.

To realize efficient LEO satellite-based trusted vessel location acquisition,we propose a TDOA-assisted DPD(TA-DPD).The TA-DPD minimizes the searching area by considering both the geometric constraints and the satellite TDOA measurements.Based on satellites FoVs and line-of-sight (LoS) radio frequency(RF) propagation assumptions, the multi-satellites common-view region is first obtained using geometric rules.Another region that contains the vessel position is constructed by using TDOAs and their variances.From the statistical theory,the true TDOA value falls into the 3-Sigma limits of the measured TDOA mean with 99.7% probability, and the lower and upper 3-Sigma limits of TDOA between two satellites project as a belt on the Earth’s surface.Therefore, the intersection of such two belts generates an area where the vessel locates.The overlap of these two regions produces the minimized searching area for the DPD,and then the emitter position can be determined by finding the cost function peak in this searching area.The peak finding process can be performed using either the conventional exhaustive approach or iterative approaches mentioned above.

The rest of this paper is organized as follows.In Section II,the LEO-satellite passive localization problem to be addressed is formulated.The localization algorithm and the searching area determination process are illustrated in detail in Section III.Then, the position estimation error lower bound i.e.CRLB for this specific localization scenario is analytically derived in Section IV.In Section V, numerical simulations and comparisons are presented to evaluate the accuracy and performance of the proposed approach.Finally,Section VI draws the conclusion.

II.PROBLEM STATEMENT

As shown in Figure 1, signals from an emitter on the Earth’s surface are captured byLsatellites deployed on the LEO.The received signals from all satellites are available for processing at a data fusion center.Coordinates of all LEO satellites are known and denoted as p1,p2,....,pL,and the location of emitter puis the unknown to be determined.Dimensions of all position vectors is 2 in two-dimension cases and 3 in three-dimension cases,i.e.pl,pu ∈R2,3.

Figure 1. An illustration of passive localization using LEO satellites for an emitter on the Earth.

Distances between the emitter on the Earth’s surface and thelthLEO satellite are expressed by

where||·||is the Euclidean norm.Moreover, since the LEO satellites have a velocity around 7km/s and the signal carrier frequencies of most emitters are normally MHz/GHz (106Hz-109Hz), the Doppler effect between the emitter and the satellite receiver is not negligible.

2.1 Signal Model

As in most cases the carrier frequencyfcof a signal is much larger than the bandwidthBof the signal waveforms(t), the transmitted signal can be expressed bys(t)e-j2πfct.Assuming the LEO satellites intercept signal transmitted by the emitter on the Earth, pland vldenote the coordinates vector and the velocity vector of thelthsatellite at the interception interval.The complex signal observed by thelthsatellite at the interception time interval,after frequency downconversion by the carrier frequencyfc,can be expressed as

in whicht ∈[-T/2,T/2]is the observation time interval ands(t)is the signal waveform during the interception interval.τl=dl/c0is the signal propagation time andblrepresents the complex signal attenuation.wlis the wide-sense stationary, zero-mean, and complex Gaussian noise.The variableflis the observed Doppler frequency shift given by

where vlis the velocity vector of thelthsatellite.To ensure the complete signal interception,τl,k ≪Tis required.Furthermore, the observation time length must be short enough so that changes of the satellite positions do not exceed the desired positioning error during the observation duration, i.e.τl,kcan be seen as constant.

2.2 Direct Position Determination

In this subsection, we briefly summarize the basic principle of the DPD and then point out the problem to be addressed in the aforementioned localization scenario using the LEO satellite cluster.Assuming each receiver collectsNtime samples in a given time interval,with sampling intervalTs,of the frequency down coverted signal(2).Discrete signal arrays are defined as follows

in whichtn=nTs,n=1,2,...,N.Write the received signal in a compact form,we have

where Flrepresents a down shifting operation which shifts s byindices.Without loss of generality,we assume the transmitted signal satisfying||s||= 1 and the covariance of wlisσ2I.Unlike conventional two-step passive localization approaches,the DPD estimates the emitter position puthrough the received signal sequencesirectly.

Using the maximum-likelihood estimator,the emitter position corresponds to the minimum value of the cost function that measures the difference between received signals and estimated received signals by (5),i.e.

The minimization (6) requires the path attenuation scalars satisfying

Substituting(7)into(6),the cost function equals to

As rlis independent of parameter p in (8), the minimization of(6)and(8)is equivalent to the maximization of the following cost function

Here,the defined Hermitian matrix Q∈ZN×Nis

For signals with known waveform,DPD estimates the emitter position by

If the signal waveform is unknown, the maximization of the cost functionC2(p)can be replaced by an equivalent one given by

whereλmaxis the largest eigenvalue andV as fully explained by[23].

Note thatflandτlused to evaluate the cost functions are all functions of the assumed emitter position p,and the position determination(11)can be achieved by grid searching.Therefore,the complexity of DPD is determined by two factors,i.e.the evaluation complexity of cost function(9)or(12)for a given p and the number of points in the area to be searched.However,for DPD of emitters using LEO satellite clusters, the grid searching process can be quite time-consuming due to the large size of satellite coverage.

To summarize,the problem discussed here is briefly stated as follow: how to find a sufficiently small area on the Earth’s surface where the emitter locates for DPD execution.

III.SEARCHING AREA DETERMINATION

This section presents a simple yet robust approach to determine the searching area on the Earth’s surface that contains the emitter.To ensure this area is as small as possible,the proposed method jointly adopts satellite common-view region and TDOA measurements.

3.1 Multi-Satellite Common-View Region

Generally speaking, the coverage footprint of a satellite on the Earth’s surface is a spherical cap, and the size of the cap is determined by the satellite orbit altitude and the FoV angle (i.e.the aperture of the receiver antenna on the satellite).As depicted in Figure 2a, the coverage footprint can be denoted by a spherical capC(D,θ) where D is the satellite subpoint, andθis the angle between vector D and the vector from the center of the Earth to a point on the coverage footprint rim.Apparently, the value ofθis determined by the radius of the Earth, the altitude of the satellite orbit, and the largest slant range of the satellite able to receive emitter signals.The intersection of two spherical caps on the Earth’s surface as shown in Figure 2b is quite complicated, and the analytical expression describing the common-view area of two satellites varies with cases as explained in[24].Let alone the analytical formula for the common-view region of more than three satellites.

Figure 2. (a)The FoV of a satellite;(b)FoVs of two satellites and their intersection.

In this work,a numerical approach is used to determine the common-view area of multiple satellites.Assuming the Earth’s surface is a perfectly spherical surface,and the it is firstly sampled evenly.The equidistributed points on the sphere surface can be achieved by choosing circles of latitude at constant intervalsdϑand on these circles points are separated bydφ,such thatdϑ ≈dφ.The density of sample points can be controlled bydϑordφwhereϑ ∈[0,π] andφ ∈[0,2π]represents the polar angle and azimuth angle,respectively.

Based on the simple geometric relation, it is straightforward to judge whether the satellite receiver can receive a RF signal transmitted from one point on the Earth’s surface.Therefore, the multi-satellite common-view area is filled with point clouds that can be observed by all satellites.To obtain the commonview boundary with the required accuracy or resolution, the boundary can be refined by generating a set of equal distribution points with higher resolution near it.Then repeating the wireless link availability judgment for a new point cloud, the refined boundary of the common view can be decided as depicted in Figure 3.For more than three satellites,the common-view determination process is similar and omitted here for brevity.

Figure 3. The boundary refining process for the multisatellite common-view area determination.

3.2 Area Estimation Based on TDOA

The area on the Earth’s surface that may contain the emitter can be established using TDOA measurements and their 3-Sigma limits.It has been shown that only two TDOA measurements from three satellites are sufficient to perform the geolocation of a target on the Earth’s surface[25], so we consider the passive positioning with three satellites.This approach converts the geolocation process to an algebraic solution of finding the positive root of a polynomial,i.e.the solution of puin terms ofd1.

wheredi,1=di -d1andR0=||pu||.In (13),d2,1andd3,1are evaluated by the speed of lightc0and the TDOA measurements Δτ2,1and Δτ3,1.Finally, the one dimensional Newton’s search is adopted to efficiently find the positive root of this fourth order polynomial and to determine pu.

In the signal processing community,there are a variety of approaches to evaluate the TDOA of a signal from two sensors.Since received signals are all available at the fusion center, the TDOA can be approximated by the cross-correlation technique owing to its accuracy and robustness.Using the popular generalized cross-correlation method with phase transformation(GCC-PHAT)[26],the TDOA Δˆτ1,2between two received signals r1and r2is

whereG1,2(ω) is the cross-spectrum of the received signal r1and r2.The performance of the crosscorrelation processing is determined by the signal bandwidth, and a signal with a wider bandwidth usually has a narrower cross-correlation peak.Theoretically, the variance of TDOA using cross-correlation[27]is

whereξeffis the effective SNR given by

andξ1andξ2are SNRs of received signals,BandTare the signal bandwidth and pulse duration, respectively.The root-mean-square(RMS)bandwidth of the transmitted signals(f)is

In this work, we use the variance of TDOA measures in an alternative way.Under rather general assumptions, the TDOA between two sensors has a multivariate normal distribution [28], thence the estimated TDOA is normally distributed with its mean equal to the true value of TDOA and the covariance matrix.Therefore, the true TDOA Δτ1,2locates inside the 3-Sigma limits of the measured TDOA,i.e.Pr(Δ-3σTDOA<Δτ1,2<Δ+3σTDOA)>99.7%.Using the lower and upper limits of Δτ1,2and Δτ1,3of three satellites and the localization procedure(13),the area on the Earth’s surface containing the emitter position with large probability can be established and we call it the“3-Sigma TDOA area”.Obviously,a larger TDOA varianceσTDOAnormally produces a larger 3-Sigma TDOA area on earth’s surface.If the signal has a higher SNR or a larger (RMS) bandwidth, the 3-Sigma TDOA area will be smaller.

3.3 Searching Area Refinement

The localization process in DPD is performed by searching the multi-satellite common-view region on the Earth’s surface.However, this searching area is still too large for efficient localization using DPD.By overlapping the area obtained from the second step with the multi-satellite common-view region using boolean operations,the searching field can be narrowed.

Furthermore,by increasing the signal duration time or equivalently enlarging the signal recording time at the sensor end, the searching area for DPD will also narrow down.For the localization scenario withL>3 satellites, theoretically()3-Sigma TDOA areas can be established in total.By finding the common region of these areas obtained from TDOA measurements and their variances,the searching area for DPD can be further refined.

Consisting of four steps, the workflow of the TADPD is illustrated in Figure 4.The common region covered by multiple LEO satellites is firstly found.The second step determines an area based on the TDOA measurements and their 3-Sigma limits.Thirdly, a minimized searching area is obtained by overlapping the region obtained in the second step with the multi-satellite common-view area.Finally,points are evenly picked over the intersection area,and the DPD is carried out.

Figure 4. Flowchart of the proposed TA-DPD algorithm.

IV.LOCALIZATION ERROR ANALYSES

The CRLB is the lower bound of the covariance for an unbiased estimator, and it is often used as a criterion to evaluate the performance of an estimation approach.In the following of this section,we derive the CRLB of passive localization using LEO satellites and estimate the theoretical accuracy that the TA-DPD can achieve.In the problem of localization,the CRLB relates to the unknown parameter vector pu

where J is the Fisher Information Matrix (FIM).Theijthentry of the FIM is given by

in whichφiis theithelement of the unknown parameter vector pu.The emitter position puon the Earth’s surface is denoted by(xu,yu,zu)in the earth-centered earth-fixed (ECEF) coordinates or by (R,θ,φ) in a spherical coordiate system as shown in Figure 5.Therefore, the FIM is a 2×2 matrix as the Earth radiusRcan be seen as a constant.The vector m is the combination of the noise-free received signals from allLreceivers given as follow

σ2in(19)is the additive noise variance level.For the sake of simplicity, we introduce a local East-North-Up (ENU) coordinate system centered at the emitter.Thus,the derivatives can be calculated by

The east and north vectors of the local ENU coordinate are given by= (-sinφ,cosφ,0) and=(-cosφcosθ,-sinφcosθ,sinθ)in the ECEF coordinates,respectively.

Using the chain rule transformation formula, the derivative of the expected receiving signal vector mlfrom thelthsatellite in m with respect to variablex′andy′are given by

and the derivative of mlwith respect to the Doppler frequencyfland time delayτlis

with

and

The partial derivative of the Doppler frequency and the time delay with respective tox′andy′are

whereare angles between the line connecting the transmitter and thelthsatellite and unit vectorrespectively,as shown in Figure 5.Here,we assume the emitter is stationary or quasi-stationary during the signal interception interval.

Figure 5. Geometries of the ECEF coordinate system and the local ENU coordinate system.

The partial derivatives of the Doppler frequency are

whereφlis the angle between the satellite velocity vector and the line connecting thelthsatellite and the transmitter.are velocities of thelthsatellite in thex′andy′axis shown in Figure 5.In summary,entries of the FIM can be written as follow

whereα1=|Rbl|2,α2=|Rsinθbl|2, andα3=sinθ|Rbl|2are signal attenuation scalars.

V.NUMERICAL SIMULATIONS

In this section,we first demonstrate the passive localization process of the proposed TA-DPD for the 2D case,and the general 3D scenario with multiple satellites is then presented.

5.1 2D Case

Consider a 2D localization problem as shown in Figure 6, in which the emitter locates at{0,15}[km].Moving sensors with RF receivers are at{20,30}[km],{-20,30}[km], and{0,0}[km], and their velocities are{2600,7030}[m/s],{-7350,-1500}[m/s], and{4780,-5780}[m/s].FoV angles of these three sensors are 19.3°,36.4°and 53.1°,respectively.This assumption holds in many applications since the receiving antennas have directional radiation patterns.We also assume that all sensors only receive signals from the LoS channel.In this way,by intersecting FoVs of three sensors, the multi-sensor common-view region is determined and highlighted as a shadow area in Figure 6.

Figure 6.The illustration of a 2D localization problem with three sensors and one emitter.The three sensors all have directional receiving patterns,and their common-view area is shown in shadow.

In simulation, the emitter to be located transmits signal with the carrier frequency of 162MHz and the bandwidth of 25KHz, and the signal is modulated by the minimum frequency shift keying (MSK).Reception SNRs at all sensors are -10dB, and the number of signal interceptionKequals to one.TDOAs Δτ1,2and Δ1,3are measured through the generalized crosscorrelation approach with received signals, and their variances are evaluated using (15) based on the receiving SNR.The 3-Sigma upper and lower limits for Δτ1,2and Δτ1,3are illustrated by the solid lines and dash lines at the inset of Figure 7.In this 2D case,the TDOA line of position is hyperbolic,and thus four hyperbolic curves can be plotted by the 3-Sigma TDOA upper and lower limits.The intersection of these hyperbolic curves also produces a shadow region that contains the emitter with large possibility as depicted in Figure 7.

Figure 7. The intersection region(in shadow)from the hyperbolic lines of 3-Sigma limits of TDOA measurements.Inset:3-Sigma limits of TDOA measurements.

The overlapping area of shadow regions in Figure 6 and Figure 7 is presented in Figure 8.Subdividing this overlapping area into grids,and evaluating the cost function on each node over the searching grid.A typical result of the cost function on the magenta shadow area is presented by the inset of Figure 8,in which the peak of the cost function is the estimated emitter location.To gather sufficient statistics, simulation results based on 100 Monte Carlo runs are presented in Figure 9 where the emitter position estimated by the conventional TDOA method and the proposed TA-DPD with-10dB SNR are marked by red dots and blue cross in this figure.It is noticed from Figure 8 and Figure 9,due to the existence of noise,some TDOA estimations are out of the multi-sensor common-view region.To evaluate the confidence level of these two methods,the elliptical error probable (EEP), an ellipse containing estimated positions by a probability of 95%, is also considered here.As we can see from Figure 9, the EEP of TDOA is 5 to 6 times larger than the counterpart of the proposed method in this example.

Figure 8. The TA-DPD searching area obtained by overlapping sensor common-view regions and the intersection region of TDOA lines of position.Inset: Cost function to determine the emitter location over the searching area.

Figure 9. Localization results of the TDOA method and TADPD under 100 Monte Carlo simulations.The EEP 95 of TDOA estimates is shown as red-dash line,and that of TADPD results is in blue dash-dot line.

5.2 3D Simulation

By demonstrating a general 3D positioning scene,this subsection explains the localization process of TADPD and evaluates its performance on multiple LEO satellites equipped with omnidirectional antennas.Assume locations and velocities of the three LEO satellites are quasi-static during the signal interception interval, and their specific values are given in Table 1.The emitter locates at [-3079.2,5532.2,700][km] in ECEF coordinate, i.e.a position with 6.3°north latitude and 119.1°east longitude in the pacific ocean.The carrier frequency, bandwidth, and modulation of the transmitted signal are identical to those in the previous example and represents a typical signal frame of AIS in nautical applications.

Table 1. Satellite positions and velocities in the ECEF coordinate.

With the proposed method, searching areas vary with SNRs of received signals.The multi-satellite common-view region on the Earth’s surface based on the LoS assumption is bounded by the outermost blue line shown in Figure 10.This figure also illustrates boundaries of the searching area with signal SNR level of -15dB, -12dB, -9dB, and -6dB.It is noticed as the received signal quality improves, the searching area gradually shrinks to 99.8%, 59.9%, 19.2%, and 5.8%of the common-view region.For the SNR higher than-9dB,in this case,the searching area is completely determined by the 3-Sigma TDOA.

Figure 10. Searching area boundaries from the threesatellite common-view region and generated by the proposed method at different SNRs.

Then, through comparing the proposed TA-DPD with TDOA method and traditional DPD, its localization accuracy and efficiency are discussed.With the signal SNR of -5dB, emitter positions estimated by TDOA approach and the proposed method are presented in Figure 11 under 100 Monte Carlo runs.From this figure, the root-means-square error (RMSE) of TDOA estimates is around 2.5 times larger than TADPD at the studied SNR level.Similar to the previous section, the contours represent the EEP curve.Taking the grid spacing as 2.5km, running times of both the TA-DPD and traditional DPD are presented in Table 2.As can be seen from the table, the higher the SNR is, the more efficient the TA-DPD exhibits over traditional DPD.It is noticed whenSNR=-15dB,the traditional DPD is slightly faster than the TA-DPD.This is because at this SNR level, the searching area refined by the proposed method has the same size as the common-view region, which can be observed in Figure 10.Yet, extra time is spent on TDOA estimation for the TA-DPD.

Figure 11. Localization errors of the TDOA method and TA-DPD for a general 3D case under 100 Monte Carlo simulations.EEP 95 of TDOA estimates is shown as red-dash line,and that of TA-DPD results is in blue dash-dot line.

Table 2. The comparison of geolocation efficiency between the proposed TA-DPD and traditional DPD.

Finally, we examine the RMSE performance of the proposed TA-DPD under different SNRs and results are given at Figure 12.RMSEs of the TDOA method are also presented here for comparison,and the CRLB derived in the previous section is presented as a reference as well.Because the proposed method is based on the DPD, its error diverges much slower than the TDOA localization method as demonstrated in the figure.

Figure 12. Localization RMSEs under different SNR levels for different localization methods.

VI.CONCLUSION

With the rapid development of global shipping and increasing illegal marine activities, the passive location technology independent of mainstream Global Navigation Satellite Systems(GNSSs)is of significant importance from the cybersecurity aspect.In this paper,the implementation of the recently developed DPD approach based on the LEO satellites is investigated.As the efficiency of DPD is determined by searching area size and because of the vast region covered by satellite,a searching area minimization approach is developed.Besides the multi-satellite common-view, another region containing the emitter based on the TDOA and its variance is determined on the Earth’s surface, and the searching area on the Earth’s surface is obtained by overlapping these two regions.As a result, the searching area size is determined by the quality of the received signal, i.e, the signal with larger SNR or wider bandwidth can be positioned more efficiently.Numerical simulations and comparisons demonstrate the validity and accuracy of the proposed TA-DPD in both 2D and realistic 3D scenarios.We believe the TA-DPD provides a legitimate tool for trusted and robust localization-based maritime service in the era of emerging LEO constellations.

VII.ACKNOWLEDGEMENT

This work was supported in part by the National Key Research and Development Program of China under Grant No.2019YFB1803200,and the National Natural Science Foundation of China(NSFC)under Grant No.61901020,and the Civil Aviation Administration of China.