CN113986906B - A Trajectory Reconstruction Method Based on Situational Targets - Google Patents

Abstract

The invention discloses a trajectory reconstruction method based on a situational target, and aims to solve the problem that the existing method cannot handle target trajectory data with uneven distribution of trajectory data points. The technical scheme is: preprocessing the original data of the situational target trajectory, that is, cleaning, converting the format and extracting the data in the same period of the target trajectory data points from different means; then setting the sampling frequency to perform a given sampling rate on the original trajectory data. Trajectory reconstruction is performed to obtain a resampling data set at the sampling rate f; finally, the trajectory data velocity vector correction is performed on the resampling data set based on equal time intervals to obtain a trajectory reconstruction data set. By adopting the invention, a trajectory reconstruction data set with uniform distribution of trajectory data points under a custom sampling rate can be obtained, the noise of the original data of the target trajectory is reduced, the fitting error of the sparse target trajectory data can be reduced, and the target trajectory data can be divided into three parts. The trajectory lines on the digital earth are clear and ensure high accuracy of the velocity vector.

Description

A Trajectory Reconstruction Method Based on Situational Goals

technical field

The invention relates to the field of geographic information situation, in particular to a method for reconstructing target track data.

Background technique

There are many monitoring methods and data sources for situational target trajectory data, resulting in a large amount of data and a high degree of complexity. Situational target trajectory data includes target state data monitored by different means at the same time, including three-dimensional coordinates in geographic space, target action speed, etc. Although visualization tools such as 3D digital earth can effectively display multi-dimensional space-time target trajectory data, due to the source of trajectory data Due to different methods and different acquisition frequencies, the visualized trajectory data of the 3D Digital Earth is chaotic, and it is impossible to observe the trajectory rules. In order to enhance the observation effect of the target's action trajectory and enhance the utilization value of the target trajectory data, it is necessary to clean and reconstruct the target trajectory data.

At present, the processing of trajectory data obtained by multiple methods and multiple acquisition frequencies is relatively rough. Generally, the original data is directly displayed on the 3D digital earth, and the dense points formed by the staff observing the original data are connected in time series to estimate the movement of the target. Behavior, this method is limited by the intensity of data and the experience of observers: when the original trajectory point data is unevenly distributed (referring to a certain period of time or a certain moment there are a large number of trajectory data points, a certain period of time a certain moment the number of data points is too large less or missing), it is more difficult for the staff to predict the target trajectory movement behavior and the accuracy of the prediction is difficult to control.

Existing trajectory data cleaning and reconstruction methods mainly include data noise reduction and curve fitting. Due to the fact that these methods process data with a single dimension and require uniform data distribution, the processing of data obtained by various means and multiple acquisition frequencies is difficult. The results of a large amount of trajectory data are poor, especially the problem of uneven distribution of trajectory data points: such as the commonly used filter fitting method, when the distribution of data points is uneven, a large amount of sawtooth noise will be generated in the dense data, and the fitting error in the sparse data Too large, which exacerbates the sparse distribution of some data, causing a large number of overlaps of some data. Therefore, conventional data cleaning and reconstruction methods to deal with unevenly distributed data often make the reconstructed data unable to represent the original data information, which is not conducive to other later algorithms. Analyze, mine and use data.

How to solve the problem that existing trajectory data cleaning and reconstruction methods cannot handle target trajectory data with uneven distribution of trajectory data points is a technical issue that is of great concern to those skilled in the art.

Contents of the invention

The technical problem to be solved by the present invention is to provide a trajectory data reconstruction method based on situational targets, which solves the problem that the existing methods cannot handle target trajectory data with uneven distribution of trajectory data points, reduces the noise in data-intensive places, and reduces the Fitting error where data is sparse.

In order to solve the above-mentioned technical problems, the technical solution of the present invention is: preprocess the original data of the target trajectory of the situation, that is, clean the target trajectory data points from different means, convert the format and extract the data in the same period; then set the sampling frequency to use The curve interpolation algorithm reconstructs the trajectory of the original trajectory data with a given sampling rate, and obtains the resampled data set at the sampling rate f; finally, based on the equal time interval, the trajectory data velocity vector correction is performed on the resampled data set, and the trajectory reconstruction data is obtained Set D.

The present invention comprises following four steps:

The first step is to build a situational target trajectory data reconstruction system. The situational target trajectory data reconstruction system is composed of a data preprocessor, a trajectory resampler, and a speed corrector.

The data preprocessor is connected with the trajectory resampler. The data preprocessor reads the original data set of the situational target trajectory Data-original from the original data set of the target trajectory, performs data cleaning, data format conversion, and data extraction at the same time period on the Data-original. Get the preprocessed data set Data-precision, and send the Data-precision to the trajectory resampler.

The trajectory resampler is connected with the data preprocessor and the speed corrector, receives the preprocessed data set Data-precision from the data preprocessor, and uses the cubic spline curve difference algorithm to reconstruct the trajectory with a given sampling rate for the Data-precision , get the resampled data set Data-sample, and send the Data-sample to the speed corrector.

The velocity corrector is connected with the trajectory resampler, receives the Data-sample from the trajectory resampler, performs trajectory data velocity vector correction on the resampled data points in the Data-sample, and obtains the trajectory reconstruction data set D.

The second step, the data preprocessor preprocesses the original data set of the situational target trajectory Data-original, including cleaning the situational target trajectory data, format conversion, and data extraction at the same time period, to obtain the target sampling period [T ₁ , T ₂ ] The trajectory data set Data-precision in the method is as follows:

2.1: The data preprocessor reads the original data set Data-original of the target trajectory from the original data set of the target trajectory,

Data-original＝{(t _i ,long _i ,lat _i ,alt _i ,v _i ,vθ _i )|i=1,...,I} (1)

Data-original is the target in the sampling period [T ₁ ,T ₂ ](t _i ∈[T ₁ ,T ₂ ], T ₁ is the starting point of the sampling period set by the user, T ₂ is the end point of the sampling period set by the user, and sampling refers to The process of sampling the behavioral state of the situational target, the sampling period refers to the time interval for sampling the behavioral state of the situational target), a collection of I data points collected by various means, I is a positive integer, and the i-th The sampling point data contains 6 dimensions, and each dimension is:

Time t _i (data format: "year-month-day hour: minute: second");

The longitude coordinate long _i of the data point, (data format: "(E/W) degrees: minutes: seconds")

The latitude coordinate lat _i of the data point (data format: "(N/S) degree: minute: second");

Altitude alt _i , indicating the altitude above sea level, in meters;

Velocity v _i represents the velocity of the target at a certain moment, and the unit is km/h;

Velocity direction vθ _i indicates the velocity direction of the target at a certain moment (data format: "degree:minute:second"), the direction of true north is 0, and the direction of clockwise is positive);

2.2: The data preprocessor performs data cleaning on Data-original, removes missing data in longitude, latitude or height fields, and obtains the cleaned data set Data-noNAN, the method is;

2.2.1: Let i=1;

2.2.2: Determine whether the i-th data point in the Data-original has missing longitude, latitude, or altitude fields. If any of these three items is missing, delete the i-th data point from the Data-original. If If none of these three items is missing, save the i-th data point to Data-noNAN;

2.2.3: let i=i+1; if i≤I, go to 2.2.2; if i>I, it means that the elimination is completed, and the cleaned data set Data-noNAN is obtained;

Data-noNAN＝{(t _in ,long _in ,lat _in ,alt _in ,v _in ,vθ _in )|in=1,...,IN} (2)

IN is the number of data points in Data-noNAN, 1≤IN≤I.

2.3: The data preprocessor converts the format of the data points in Data-noNAN: converts the time to a timestamp, and converts longitude and latitude to double precision in radians (north latitude is positive, south latitude is negative; east longitude is positive , west longitude is negative), convert the altitude to a double in meters (sea level is 0), convert the velocity magnitude to a double in m/s, and convert the velocity direction to a double in radians ( North is 0, clockwise is positive), the method is:

2.3.1: let in=1;

2.3.2: Use time conversion function ft, longitude conversion function flong, latitude conversion function flat, height conversion function fal, speed size conversion function fv, speed direction conversion function fvθ to the in-th data point in Data-noNAN (t _in , long _in ,lat _in ,alt _in ,v _in ,vθ _in ) are converted separately:

Make the converted moment t _in ′=ft(t _in ); ft(t _in ) means to use the java time tool class to convert the format of t _in : "year-month-day hour: minute: second" into a timestamp ( Timestamp (that is, the total number of seconds from GMT to the current time));

Make the converted longitude coordinate long _in ′=flong(long _in ), and flong(long _in ) means that the longitude data format of long _in : "(E/W) degree: minute: second" is converted into a double-precision floating-point number, and the unit is In radians (accuracy: retain three decimal places, that is, convert the angle unit from degrees, minutes, and seconds to the international unified unit of radians, (E/W) means east longitude and west longitude), the method is: take the prime meridian as 0, and the east longitude as Positive, west meridian is negative, determine whether long _in is positive or negative, if negative, add "-" before long _in , after converting to radian system, keep the accuracy of long _in to three decimal places;

Make the converted latitude coordinate lat _in '=flat(lat _in ), flat(lat _in ) means to convert the latitude data format of lat _in : "(N/S) degree: minute: second" into a double-precision floating-point number, unit In radians (accuracy: keep three decimal places, that is, convert the angle unit from degrees, minutes and seconds to the international unified unit of radians, (N/S) means north latitude and south latitude), the method is: take the equator as 0, and the north latitude as positive, If the south latitude is negative, determine whether lat _in is positive or negative. If it is negative, add "-" before lat _in , convert to radian system and keep the accuracy of lat _in to three decimal places;

Make the converted height alt _in '=fal(alt _in ), fal(alt _in ) means that the height alt _in is converted into the double precision with the unit of meter, the method is: take sea level as 0, above sea level as positive, sea level Below the plane is negative, determine whether alt _in is positive or negative, if negative, add "-" before alt _in , and keep the accuracy of alt _in to three decimal places;

Make the converted speed v _in '=fv(v _in ), fv(v _in ) means that the speed v _in (in km/h) is converted into a double-precision floating-point number (in m/s, precision: reserved three decimal places);

Make the converted velocity direction vθ _in ′=fvθ(vθ _in ), fvθ(vθ _in ) means converting the velocity direction vθ _in (data format: "degree:minute:second") into a double-precision floating-point number (that is, the angle unit is degrees, minutes and seconds into the international unified unit radian system), the method is: the direction of true north is 0, clockwise is positive, counterclockwise is negative, determine whether vθ _in is positive or negative, if negative, add "-" before vθ _in , and Keep the precision of alt _in to four decimal places;

Save the converted in-th data point (ie (t _in ′, long _in ′, lat _in ′, alt _in ′, v _in ′, vθ _in ′) to the converted data Data-format;

2.3.3: let i=i+1; if i≤IN, go to 2.3.2; if i>IN, it means that the conversion is completed, and the data set Data-format after format conversion is obtained;

Data-format={((t _in ′,long _in ′,lat _in ′,alt _in ′,v _in ′,vθ _in ′))|in=1,…,IN} (3)

2.4: Extract the data of the same time period within [T ₁ , T ₂ ] for Data-format. Because the data is collected by multiple methods, there are often a large number of data points at the same time or in a short time range, and the accuracy of these data points is different (for example: time data, some are accurate to the second, and some are accurate to the minute, resulting in The number of decimal places is different), so it is necessary to extract data from Data-format in the same time period (the same time period is the same time or a shorter time range) within [T ₁ , T ₂ ], the method is:

Divide the sampling period [T ₁ , T ₂ ] into P time periods dt (dt generally takes 60 seconds), obtain data points in each time period from Data-format, compare the accuracy of each data point, and take the precision dimension PN (The number of decimal places satisfying the number of dimensions of decimal places after format conversion in step 2.3 is the number of precision dimensions), and the data point with the largest PN is used as the feature point in the dt. The specific method is as follows:

2.4.1: Let the variable p=1;

2.4.1.1: Let there be L+1 data points in the pth time period dt _p of Data-format, and let these L+1 data points form the data point set DD _{p in dt p :}

DD _p ＝{D _p0 ,D _p1 ,…,D _pL |st:max(t _p0 ,t _p1 ,…,t _pL )-min(t _p0 ,t _p1 ,…,t _pL )≤dt _p } (4 )

D _pl is the l+1th data point, and the data points placed in DD _p satisfy the maximum value max(t _p0 ,t _p1 ,…,t _pL ) of the L+1 data points in DD _p and L in DD _p The difference between the minimum value min(t _p0 ,t _p1 ,…,t _pL ) of +1 data point is within the range of dt _p .

2.4.1.2: Calculate the number of precision dimensions PN of D _p0 , D _p1 ,..., D _pL :

(Because the time dimension in (t _in ′, long _in ′, lat _in ′, alt _in ′, v _in ′, vθ _in ′) is an integer. It is not necessary to calculate the effective digits of precision, so the maximum PN is 5), only calculate Long _in ′, lat _in ′, alt _in ′, v _in ′, vθ _in ′, the precision dimensions of these 5 dimensions, the method is:

2.4.1.2.1 Let the variable l=0;

2.4.1.2.2 Check the number of effective digits of data precision of long _in ′, lat _in ′, alt _in ′, v _in ′, vθ _in ′ of D _pl , if there are w dimensions that satisfy the decimal point after format conversion in step 2.3 digits, then let the precision dimension of D _pl be

转2.4.1.3；2.4.1.2.3 Let l=l+1; if l≤L, go to 2.4.1.2.2; if l>L, the calculation is completed, and the number of data precision dimensions of D _p0 , D _p1 ,...,D _pl is

Go to 2.4.1.3;

中的最大值，令该最大值的数据点为D_pmax：2.4.1.3: Found

The maximum value in , let the data point of the maximum value be D _p max:

D _p max＝(dpt _p ,dplong _p ,dplat _p ,dpalt _p ,dpv _p ,dpvθ _p ) (5)

Among them, dpt p , dplong _p , dplat _p , dpalt _p , dpv _p , and dpvθ _p are the _values of the six dimensions of the data points corresponding to the maximum number of precision dimensions in the p-th time period dt _p .

2.4.2: Put D _p max in the trajectory data set Data-precision;

2.4.3: Let p=p+1; if p≤P, go to 2.4.1.1; if p>P, it means that the data extraction in the same period is completed, and the trajectory data set of the target within the time period [T ₁ , T ₂ ] is obtained Data-precision:

Where P is the number of time segments divided by the sampling period [T ₁ , T ₂ ], that is, the number of trajectory data in Data-precision.

In the third step, the trajectory resampler reconstructs the trajectory with a given sampling rate for Data-precision based on the cubic spline curve difference algorithm, and obtains the resampled data set Data-sample at the sampling rate f. The method is:

3.1: Set the sampling rate f for [T ₁ , T ₂ ] (the sampling rate is generally taken as 1, that is, data is collected every 1 second) (let the collection time interval t _f =1/f), generate IS-1, etc. At the time of the interval time point, IS is a positive integer, which is the sample capacity of the time set T _f . Put these IS-1 equal interval time points and T ₁ , T ₂ into the time set T _f of the sampling rate f, and the last Take the minimum value of T ₁ + (IS-1)×t _f , T ₂ at all times:

T _f ＝{T ₁ ,T ₁ +1×t _f ,T ₁ +2×t _f ,…,min(T ₁ +(IS-1)×t _f ,T ₂ )} (7)

3.2: Use cubic spline curve difference algorithm (cubic spline interpolation algorithm see: [1] Song Yang. Research on ACARS track reconstruction algorithm based on adaptive cubic spline interpolation. Civil Aviation University of China, 2017.) to calculate T _f The longitude, latitude, height, velocity, and velocity direction corresponding to each moment in , form a resampling data set under the sampling rate f, the specific method is:

3.2.1 Select data from Data-precision and combine them into two-dimensional data (X, Y), where X={dpt _p |p=1,...,P} is time data, Y={Y _k |k=1 ,...,5} contains other 5 dimension data;

3.2.2: let k=1;

3.2.3: For the data set (X,Y _k ), set P data points as (x ₀ ,y ₀ ),(x ₁ ,y ₁ ),…,(x _P-1 ,y _P-1 ), Use the cubic spline difference algorithm to solve the cubic spline coefficient vector (a _q ,b _q ,c _q ,d _q ) of each curve segment, so that in each interval x _q ≤x≤x _q+1 , q= 0,1,…,P-1, the spline function for data fitting is:

f _q (x)＝(a _q ,b _q ,c _q ,d _q )×(1,(xx _q ),(xx _q ) ² ,(xx _q ) ³ ) ^T (9)

3.2.4: Bring the sampling time into formula (9), calculate the fitting result, and store the fitting result in the resampling data set Data-sample, the method is:

3.2.4.1: let is=1,

3.2.4.2 Take the is-th data of time set T _f , let it be t _is , t _is =T ₁ +(is-1)×t _f ;

3.2.4.3: Calculate the k-th dimension data at time t _is

保存于拟合后矩阵

的k行第is个数据上；3.2.4.4: Will

Save in the fitted matrix

On the is-th data of k row;

3.2.4.5: Let is=is+1:

If is≤IS, the k-th dimension data has not been fitted and calculated, go to 3.2.4.2,

If is>IS, the k-th dimension data fitting calculation is completed, go to 3.2.4.6;

3.2.4.6: Let k=k+1:

If k≤5, go to 3.2.3;

为：If k>5, all dimension data fitting calculations are completed, and a fitting matrix is generated

for:

得到采样率f下的重采样数据集Data-sample：Fitted data based on time data T _f and other 5 dimensions

Get the resampled data set Data-sample at the sampling rate f:

Data-sample＝{(dst _is ,dslong _is ,dslat _is ,dsalt _is ,dsv _is ,dsvθ _is )|t _is ∈T _f ,is=1,...,IS}(12)

Moment dst _is =T _{f, is} (T _{f, is} represents the is data of T _f );

dslong _is ,dslat _is ,dsalt _is ,dsv _is ,dsvθ _is are the results of data fitting for each dimension, namely

In the fourth step, the speed corrector corrects the trajectory data velocity vector of the Data-sample based on equal time intervals to obtain the trajectory reconstruction data set D. The method is:

The data set Data-sample is the data at equal time intervals, and the velocity vector is the resampling of the monitoring data. In practice, the velocity vector can be calculated according to time, longitude, latitude, and height. In order not to introduce calculation errors, the existing The method generally observes the velocity vector through a sensor, but the velocity vector collected by the sensor has large fluctuations and cannot be directly used for data analysis. Perform comparison and correction to obtain better velocity vector data, the specific steps are as follows:

4.1: Calculate the space velocity of each data point according to the Data-sample, and obtain the space velocity resampling data set, the method is:

4.1.1: let is=1;

4.1.2 According to the is data point in Data-sample

与空间速度方向

方法如下：(dst _is ,dslong _is ,dslat _is ,dsalt _is ,dsv _is ,dsvθ _is ) time, longitude, latitude and height value calculation to get the space velocity

and space velocity direction

Methods as below:

4.1.2.1 Use the Haversine formula to calculate the is-th data point and js-th data point according to latitude and longitude (dst _js , dslong _js , dslat _js , dsalt _js , dsv _js , dsvθ _js ) (js=is+λ, λ is the step size , the default is 1) the distance on the horizontal plane, the Haversin formula is:

Among them: haversin(θ)=sin ² (θ/2)=(1-cos(θ))/2, R is the radius of the earth, which is 6371km;

The distance d _is,js on the horizontal plane between the data point of the is item and the data point of the js item is calculated as:

与空间速度方向

方法为：4.1.2.2 Calculation of space velocity

and space velocity direction

The method is:

a, b are the parameters of the calculation process, and the parameter a is cos(lat _is+1 )sin(long _is+1 -long _is );

The parameter b is

放到空间速度重采样数据集

4.1.3 Will

into the space velocity resampling dataset

4.1.4 Let is=is+1; if is≤IS, go to 4.1.2; if is>IS, it means that the calculation is completed, and the speed resampling data set is obtained:

4.2: Correct the velocity vector according to the set error rate (η), the method is:

Based on the calculated space velocity vector as the standard, correct the velocity vector collected in the Data-sample. If the collected velocity vector is within the error rate range, keep the collected velocity vector; otherwise, use the calculated space velocity vector to replace the data in the Data-sample. Acquired velocity vector.

4.2.1 The error rate h for receiving user input from the keyboard is generally 0.05;

4.2.2: let is=1;

中第is项数据点的空间速度大小

与空间速度方向

计算dsv_is和

的差距、dsvθ_is和

的差距，根据差距与h的关系，确定修正后的速度矢量的大小v_is*和速度方向vθ_is*，方法如公式(18)所示：4.2.3 Take the velocity magnitude dsv _is and the velocity direction dsvθ _is of the data point of the is item of Data-sample, and take

The space velocity of the is item data point in

and space velocity direction

Calculate dsv _is and

The gap, dsvθ _is and

According to the relationship between the gap and h, determine the magnitude of the velocity vector v _is * and the velocity direction vθ _is * after correction, the method is shown in formula (18):

4.2.4 Use (v _is *, vθ _is *) to replace (dsv _is , dsvθ _is ) in the data set Data-sample, and replace (dst _is , dslong _is , dslat _is , dsalt _is , v _is *, vθ _is * ) into the trajectory reconstruction data set D.

4.2.5 Let is=is+1; if is≤IS, go to 4.2.2; if is>IS, it means that the replacement is completed, and the trajectory reconstruction data with the sampling rate f set in the target time period [T ₁ , T ₂ ] is obtained , including trajectory reconstruction datasets D and T _f :

By adopting the present invention to process the original trajectory data, the trajectory reconstruction data set under the set sampling rate can be obtained. Adopt the present invention can reach following technical effect:

1. By adopting the present invention, a trajectory reconstruction data set with uniform distribution of trajectory data points under a self-defined sampling rate can be obtained, which solves the problem that existing methods cannot process target trajectory data with uneven distribution of trajectory data points.

2. In the second step of the present invention, the data preprocessor adopts the data extraction method of the same time period to filter a large number of data points of the same time period, thereby reducing the noise at places where the original data of the target trajectory is dense.

3. The trajectory resampler in the third step of the present invention uses a trajectory reconstruction method with a given sampling rate to reconstruct the trajectory, and uses the set sampling rate to make the target trajectory data points evenly distributed in the time dimension, reducing the sparseness of the target trajectory data The fitting error at , makes the trajectory line of the fitted target trajectory data clear on the three-digit Earth.

4. The speed corrector in the fourth step of the present invention compares and corrects the calculated space velocity vector and the velocity vector observed and collected by the sensor based on equal time intervals to obtain better velocity vector data and ensure high velocity vector accuracy.

The trajectory reconstruction data set obtained by the present invention can help researchers to better use trajectory data for data analysis, help front-line staff to better observe the clear process of target activities, enhance the observation effect of the target's action trajectory, and enhance Utilization value of target trajectory data.

Description of drawings

Fig. 1 is the overall flow chart of the present invention.

Fig. 2 is a logical structure diagram of the situation target trajectory reconstruction system constructed in the first step of the present invention.

Detailed ways

The present invention will be further described in detail below in conjunction with the accompanying drawings and specific embodiments.

Fig. 1 is the general flowchart of the present invention, including four steps.

The first step is to build a situational target trajectory data reconstruction system. The situational target trajectory data reconstruction system is shown in Figure 2, which is composed of a data preprocessor, a trajectory resampler, and a speed corrector.

The data preprocessor is connected with the trajectory resampler. The data preprocessor reads the original data set of the situational target trajectory Data-original from the original data set of the target trajectory, performs data cleaning, data format conversion, and data extraction at the same time period on the Data-original. Get the preprocessed data set Data-precision, and send the Data-precision to the trajectory resampler.

The trajectory resampler is connected with the data preprocessor and the speed corrector, receives the preprocessed data set Data-precision from the data preprocessor, and uses the cubic spline curve difference algorithm to reconstruct the trajectory with a given sampling rate for the Data-precision , get the resampled data set Data-sample, and send the Data-sample to the speed corrector.

The velocity corrector is connected with the trajectory resampler, receives the Data-sample from the trajectory resampler, performs trajectory data velocity vector correction on the resampled data points in the Data-sample, and obtains the trajectory reconstruction data set D.

The second step, the data preprocessor preprocesses the original data set of the situational target trajectory Data-original, including cleaning the situational target trajectory data, format conversion, and data extraction at the same time period, to obtain the target sampling period [T ₁ , T ₂ ] The trajectory data set Data-precision in the method is as follows:

2.1: The data preprocessor reads the original data set Data-original of the target trajectory from the original data set of the target trajectory,

Data-original＝{(t _i ,longi,lat _i ,alt _i ,v _i ,vθ _i )|i=1,...,I} (2)

Data-original is the target in the sampling period [T ₁ ,T ₂ ](t _i ∈[T ₁ ,T ₂ ], T ₁ is the starting point of the sampling period set by the user, T ₂ is the end point of the sampling period set by the user, and sampling refers to The process of sampling the behavioral state of the situational target, the sampling period refers to the time interval for sampling the behavioral state of the situational target), a collection of I data points collected by various means, I is a positive integer, and the i-th The sampling point data contains 6 dimensions, and each dimension is:

Time t _i (data format: "year-month-day hour: minute: second", such as 2021-08-15 12:20:35);

The longitude coordinate long _i of the data point, (data format: "(E/W) degrees: minutes: seconds")

The latitude coordinate lat _i of the data point (data format: "(N/S) degree: minute: second", such as N23:35:15 north latitude 23 degrees 35 minutes and 15 seconds);

Altitude alt _i , indicating the altitude above sea level, in meters;

Velocity v _i represents the velocity of the target at a certain moment, and the unit is km/h;

Velocity direction vθ _i indicates the velocity direction of the target at a certain moment (data format: "degree:minute:second"), the true north direction is 0, clockwise is positive, for example, 23:35:15 is 23 degrees 35 minutes east of north 15 seconds);

2.2: The data preprocessor performs data cleaning on Data-original, removes missing data in longitude, latitude or height fields, and obtains the cleaned data set Data-noNAN, the method is;

2.2.1: Let i=1;

2.2.2: Determine whether the i-th data point in the Data-original has missing longitude, latitude, or altitude fields. If any of these three items is missing, delete the i-th data point from the Data-original. If If none of these three items is missing, save the i-th data point to Data-noNAN;

2.2.3: let i=i+1; if i≤I, go to 2.2.2; if i>I, it means that the elimination is completed, and the cleaned data set Data-noNAN is obtained;

Data-noNAN＝{(t _in ,long _in ,lat _in ,alt _in ,v _in ,vθ _in )|in=1,...,IN} (2)

IN is the number of data points in Data-noNAN, 1≤IN≤I.

2.3: The data preprocessor converts the format of the data points in Data-noNAN: converts the time to a timestamp, and converts longitude and latitude to double precision in radians (north latitude is positive, south latitude is negative; east longitude is positive , west longitude is negative), convert the altitude to a double in meters (sea level is 0), convert the velocity magnitude to a double in m/s, and convert the velocity direction to a double in radians ( North is 0, clockwise is positive), the method is:

2.3.1: let in=1;

2.3.2: Use time conversion function ft, longitude conversion function flong, latitude conversion function flat, height conversion function fal, speed size conversion function fv, speed direction conversion function fvθ to the in-th data point in Data-noNAN (t _in , long _in ,lat _in ,alt _in ,v _in ,vθ _in ) are converted separately:

Make the converted moment t _in ′=ft(t _in ); ft(t _in ) means to use the java time tool class to convert the format of t _in : "year-month-day hour: minute: second" into a timestamp ( Timestamp (the total number of seconds from Greenwich Mean Time to the current time), for example: "2021-10-01 12:12:12" is converted to a timestamp: 1633061532000;

Make the converted longitude coordinate long _in ′=flong(long _in ), and flong(long _in ) means that the longitude data format of long _in : "(E/W) degree: minute: second" is converted into a double-precision floating-point number, and the unit is In radians (accuracy: retain three decimal places, that is, convert the angle unit from degrees, minutes, and seconds to the international unified unit of radians, (E/W) means east longitude and west longitude), the method is: take the prime meridian as 0, and the east longitude as Positive, west longitude is negative, determine whether long _in is positive or negative, if negative, add "-" before long _in , and after converting to radian system, keep the accuracy of long _in to three decimal places, such as the longitude of Tiananmen Square in Beijing (East longitude: 116°23′17″) converted to (2.030rad);

Make the converted latitude coordinate lat _in '=flat(lat _in ), flat(lat _in ) means to convert the latitude data format of lat _in : "(N/S) degree: minute: second" into a double-precision floating-point number, unit In radians (accuracy: keep three decimal places, that is, convert the angle unit from degrees, minutes and seconds to the international unified unit of radians, (N/S) means north latitude and south latitude), the method is: take the equator as 0, and the north latitude as positive, South latitude is negative, determine whether lat _in is positive or negative, if it is negative, add "-" before lat _in , convert to radian system and keep the accuracy of lat _in to three decimal places, such as the latitude of Tiananmen Square in Beijing (North latitude: 39° 54′27″) converted to (0.696rad);

Make the converted height alt _in '=fal(alt _in ), fal(alt _in ) means that the height alt _in is converted into the double precision with the unit of meter, the method is: take sea level as 0, above sea level as positive, sea level Below the plane is negative, determine whether alt _in is positive or negative, if negative, add "-" before alt _in , and keep the accuracy of alt _in to three decimal places, for example, the altitude of 122.455m is converted to 122.455m above sea level, and the depth of 122.455m is converted to Altitude -122.455m;

Make the converted speed v _in '=fv(v _in ), fv(v _in ) means that the speed v _in (in km/h) is converted into a double-precision floating-point number (in m/s, precision: reserved Three digits after the decimal point), for example, 360.000km/h is converted to 100.000m/s;

Make the converted velocity direction vθ _in ′=fvθ(vθ _in ), fvθ(vθ _in ) means converting the velocity direction vθ _in (data format: "degree:minute:second") into a double-precision floating-point number (that is, the angle unit is degrees, minutes and seconds into the international unified unit radian system), the method is: the direction of true north is 0, clockwise is positive, counterclockwise is negative, determine whether vθ _in is positive or negative, if negative, add "-" before vθ _in , and Keep the accuracy of alt _in to four decimal places, for example, convert 23 degrees 35 minutes and 15 seconds north to east into 0.4115rad;

Save the converted in-th data point (ie (t _in ′, long _in ′, lat _in ′, alt _in ′, v _in ′, vθ _in ′) to the converted data Data-format;

2.3.3: let i=i+1; if i≤IN, go to 2.3.2; if i>IN, it means that the conversion is completed, and the data set Data-format after format conversion is obtained;

Data-format={((t _in ′,long _in ′,lat _in ′,alt _in ′,v _in ′,vθ _in ′))|in=1,…,IN} (3)

2.4: Extract the data of the same time period within [T ₁ , T ₂ ] for Data-format. Because the data is collected by multiple methods, there are often a large number of data points at the same time or in a short time range, and the accuracy of these data points is different (for example: time data, some are accurate to the second, and some are accurate to the minute, resulting in The number of decimal places is different), so it is necessary to extract data from Data-format in the same time period (the same time period is the same time or a shorter time range) within [T ₁ , T ₂ ], the method is:

Divide the sampling period [T ₁ , T ₂ ] into P time periods dt (dt generally takes 60 seconds), obtain data points in each time period from Data-format, compare the accuracy of each data point, and take the precision dimension PN (The dimension whose number of decimal points satisfies the agreed number of valid digits is the number of precision dimensions, such as the data of timestamp 1633061532000 (2.030, 0.696, 122.455, 100.000, 0.4115), the number of decimal places of the data in the five dimensions all meet the 2.3 step format The number of significant digits after conversion, the number of precision dimensions PN is 5); the data point with the largest PN is used as the feature point in the dt, the specific method is as follows:

2.4.1: Let the variable p=1;

2.4.1.1: Let there be L+1 data points in the pth time period dt _p of Data-format, and let these L+1 data points form the data point set DD _{p in dt p :}

DD _p ＝{D _p0 ,D _p1 ,…,D _pL |st:max(t _p0 ,t _p1 ,…,t _pL )-min(t _p0 ,t _p1 ,…,t _pL )≤dt _p } (4 )

D _pl is the l+1th data point, and the data points placed in DD _p satisfy the maximum value max(t _p0 ,t _p1 ,…,t _pL ) of the L+1 data points in DD _p and L in DD _p The difference between the minimum value min(t _p0 ,t _p1 ,…,t _pL ) of +1 data point is within the range of dt _p .

2.4.1.2: Calculate the number of precision dimensions PN of D _p0 , D _p1 ,..., D _pL :

(Because the time dimension in (t _in ′, long _in ′, lat _in ′, alt _in ′, v _in ′, vθ _in ′) is an integer. It is not necessary to calculate the effective digits of precision, so the maximum PN is 5), only calculate Long _in ′, lat _in ′, alt _in ′, v _in ′, vθ _in ′, the precision dimensions of these 5 dimensions, the method is:

2.4.1.2.1 Let the variable l=0;

2.4.1.2.2 Check the number of effective digits of data precision of long _in ′, lat _in ′, alt _in ′, v _in ′, vθ _in ′ of D _pl , if there are w dimensions that satisfy the decimal point after format conversion in step 2.3 number of digits, then let the number of precision dimensions of D _pl be

转2.4.1.3；2.4.1.2.3 Let l=l+1; if l≤L, go to 2.4.1.2.2; if l>L, the calculation is completed, and the number of data precision dimensions of D _p0 , D _p1 ,...,D _pl is

Go to 2.4.1.3;

中的最大值，令该最大值的数据点为D_pmax：2.4.1.3: Found

The maximum value in , let the data point of the maximum value be D _p max:

D _p max＝(dpt _p ,dplong _p ,dplat _p ,dpalt _p ,dpv _p ,dpvθ _p ) (5)

Among them, dpt p , dplong _p , dplat _p , dpalt _p , dpv _p , and dpvθ _p are the _values of the six dimensions of the data points corresponding to the maximum number of precision dimensions in the p-th time period dt _p .

2.4.2: Put D _p max in the trajectory data set Data-precision;

2.4.3: Let p=p+1; if p≤P, go to 2.4.1.1; if p>P, it means that the data extraction in the same period is completed, and the trajectory data set of the target within the time period [T ₁ , T ₂ ] is obtained Data-precision:

Where P is the number of time segments divided by the sampling period [T ₁ , T ₂ ], that is, the number of trajectory data in Data-precision.

In the third step, the trajectory resampler reconstructs the trajectory with a given sampling rate for Data-precision based on the cubic spline curve difference algorithm, and obtains the resampled data set Data-sample at the sampling rate f. The method is:

3.1: Set the sampling rate f for [T ₁ , T ₂ ] (the sampling rate is generally taken as 1, that is, data is collected every 1 second) (let the collection time interval t _f =1/f), generate IS-1, etc. At the time of the interval time point, IS is a positive integer, which is the sample capacity of the time set T _f . Put these IS-1 equal interval time points and T ₁ , T ₂ into the time set T _f of the sampling rate f, and the last Take the minimum value of T ₁ +(IS-1)×t _f , T ₂ at all times:

T _f ＝{T ₁ ,T ₁ +1×t _f ,T ₁ +2×t _f ,…,min(T ₁ +(IS-1)×t _f ,T ₂ )} (7)

3.2: Use cubic spline curve difference algorithm (cubic spline interpolation algorithm see: [1] Song Yang. Research on ACARS track reconstruction algorithm based on adaptive cubic spline interpolation. Civil Aviation University of China, 2017.) to calculate T _f The longitude, latitude, height, velocity, and velocity direction corresponding to each moment in , form a resampling data set under the sampling rate f, the specific method is:

3.2.1 Select data from Data-precision and combine them into two-dimensional data (X, Y), where X={dpt _p |p=1,...,P} is time data, Y={Y _k |k=1 ,...,5} contains other 5 dimension data;

3.2.2: let k=1;

3.2.3: For the data set (X,Y _k ), set P data points as (x ₀ ,y ₀ ),(x ₁ ,y ₁ ),…,(x _P-1 ,y _P-1 ), Use the cubic spline difference algorithm to solve the cubic spline coefficient vector (a _q ,b _q ,c _q ,d _q ) of each curve segment, so that in each interval x _q ≤x≤x _q+1 , q= 0,1,…,P-1, the spline function for data fitting is:

f _q (x)＝(a _q ,b _q ,c _q ,d _q )×(1,(xx _q ),(xx _q ) ² ,(xx _q ) ³ ) ^T (9)

3.2.4: Bring the sampling time into formula (9), calculate the fitting result, and store the fitting result in the resampling data set Data-sample, the method is:

3.2.4.1: let is=1,

3.2.4.2 Take the is-th data of time set T _f , let it be t _is , t _is =T ₁ +(is-1)×t _f ;

3.2.4.3: Calculate the k-th dimension data at time t _is

保存于拟合后矩阵

的k行第is个数据上；3.2.4.4: Will

Save in the fitted matrix

On the is-th data of k row;

3.2.4.5: let is=is+1:

If is≤IS, the k-th dimension data has not been fitted and calculated, go to 3.2.4.2,

If is>IS, the k-th dimension data fitting calculation is completed, go to 3.2.4.6;

3.2.4.6: Let k=k+1:

If k≤5, go to 3.2.3;

为：If k>5, all dimension data fitting calculations are completed, and a fitting matrix is generated

for:

得到采样率f下的重采样数据集Data-sample：Fitted data based on time data T _f and other 5 dimensions

Get the resampled data set Data-sample at the sampling rate f:

Data-sample＝{(dst _is ,dslong _is ,dslat _is ,dsalt _is ,dsv _is ,dsvθ _is )|t _is ∈T _f ,is＝1,…,IS} (12)

Moment dst _is =T _{f, is} (T _{f, is} represents the is data of T _f );

dslong _is ,dslat _is ,dsalt _is ,dsv _is ,dsvθ _is are the results of data fitting for each dimension, namely

In the fourth step, the speed corrector corrects the trajectory data velocity vector of the Data-sample based on equal time intervals to obtain the trajectory reconstruction data set D. The method is:

The data set Data-sample is the data at equal time intervals, and the velocity vector is the resampling of the monitoring data. In practice, the velocity vector can be calculated according to time, longitude, latitude, and height. In order not to introduce calculation errors, the existing The method generally observes the velocity vector through a sensor, but the velocity vector collected by the sensor has large fluctuations and cannot be directly used for data analysis. Perform comparison and correction to obtain better velocity vector data, the specific steps are as follows:

4.1: Calculate the space velocity of each data point according to the Data-sample, and obtain the space velocity resampling data set, the method is:

4.1.1: let is=1;

4.1.2 According to the is data point in Data-sample

与空间速度方向

方法如下：(dst _is ,dslong _is ,dslat _is ,dsalt _is ,dsv _is ,dsvθ _is ) time, longitude, latitude and height value calculation to get the space velocity

and space velocity direction

Methods as below:

4.1.2.1 Use the Haversine formula to calculate the is-th data point and js-th data point according to latitude and longitude (dst _js , dslong _js , dslat _js , dsalt _js , dsv _js , dsvθ _js ) (js=is+λ, λ is the step size , the default is 1) the distance on the horizontal plane, the Haversin formula is:

Among them: haversin(θ)=sin ² (θ/2)=(1-cos(θ))/2, R is the radius of the earth, which is 6371km;

The distance d _is,js on the horizontal plane between the data point of the is item and the data point of the js item is calculated as:

与空间速度方向

方法为：4.1.2.2 Calculation of space velocity

and space velocity direction

The method is:

a, b are the parameters of the calculation process, and the parameter a is cos(lat _is+1 )sin(long _is+1 -long _is );

The parameter b is

放到空间速度重采样数据集

4.1.3 Will

into the space velocity resampling dataset

4.1.4 Let is=is+1; if is≤IS, go to 4.1.2; if is>IS, it means that the calculation is completed, and the speed resampling data set is obtained:

4.2: Correct the velocity vector according to the set error rate (h), the method is:

Based on the calculated space velocity vector as the standard, correct the velocity vector collected in the Data-sample. If the collected velocity vector is within the error rate range, keep the collected velocity vector; otherwise, use the calculated space velocity vector to replace the data in the Data-sample. Acquired velocity vector.

4.2.1 The error rate h for receiving user input from the keyboard is generally 0.05;

4.2.2: let is=1;

中第is项数据点的空间速度大小

与空间速度方向

计算dsv_is和

的差距、dsvθ_is和

的差距，根据差距与h的关系，确定修正后的速度矢量的大小v_is*和速度方向vθ_is*，方法如公式(18)所示：4.2.3 Take the velocity magnitude dsv _is and the velocity direction dsvθ _is of the data point of the is item of Data-sample, and take

The space velocity of the data point of the is item in

and space velocity direction

Calculate dsv _is and

The gap, dsvθ _is and

According to the relationship between the gap and h, determine the magnitude of the velocity vector v _is * and the velocity direction vθ _is * after correction, the method is shown in formula (18):

4.2.4 Use (v _is *, vθ _is *) to replace (dsv _is , dsvθ _is ) in the data set Data-sample, and replace (dst _is , dslong _is , dslat _is , dsalt _is , v _is *, vθ _is * ) into the trajectory reconstruction data set D.

4.2.5 Let is=is+1; if is≤IS, go to 4.2.2; if is>IS, it means that the replacement is completed, and the trajectory reconstruction data with the sampling rate f set in the target time period [T ₁ , T ₂ ] is obtained , including trajectory reconstruction datasets D and T _f :

Claims (7)

1. A trajectory reconstruction method based on situational target, is characterized in that comprising the following steps: The first step is to build a situational target trajectory data reconstruction system. The situational target trajectory data reconstruction system is composed of a data preprocessor, a trajectory resampler, and a speed corrector; The data preprocessor is connected with the trajectory resampler. The data preprocessor reads the original data set of the situational target trajectory Data-original from the original data set of the target trajectory, performs data cleaning, data format conversion, and data extraction at the same time period on the Data-original. Get the preprocessed data set Data-precision, and send the Data-precision to the trajectory resampler; The trajectory resampler is connected with the data preprocessor and the speed corrector, receives the preprocessed data set Data-precision from the data preprocessor, reconstructs the trajectory with a given sampling rate for the Data-precision, and obtains the resampled data set Data- sample, send the Data-sample to the speed modifier; The velocity corrector is connected to the trajectory resampler, receives the Data-sample from the trajectory resampler, and corrects the trajectory data velocity vector for the resampled data points in the Data-sample to obtain the trajectory reconstruction data set D; The second step, the data preprocessor preprocesses the original data set of the situational target trajectory Data-original, including cleaning the situational target trajectory data, format conversion, and data extraction at the same time period, to obtain the target sampling period [T ₁ , T ₂ ] The trajectory data set Data-precision in the method is as follows: 2.1: The data preprocessor reads the original data set Data-original of the target trajectory from the original data set of the target trajectory, Data-original＝{(t _i ,long _i ,lat _i ,alt _i ,v _i ,vθ _i )|i=1,...,I} (1) Data-original is a collection of I data points collected by various means within the sampling period [T ₁ , T ₂ ], I is a positive integer, t _i ∈ [T ₁ , T ₂ ], T ₁ is the starting point of the sampling period set by the user, T ₂ is the end point of the sampling period set by the user, sampling refers to the process of sampling the action state of the situational target, and the sampling period refers to the time interval for sampling the action state of the situational target; the i-th sampling The point data contains 6 dimensions, and each dimension is respectively: time t _i , longitude coordinate long _i of the data point, latitude coordinate lat _i of the data point, height alt _i , velocity magnitude v _i ; velocity direction vθ _i ; 2.2: The data preprocessor performs data cleaning on Data-original, removes missing data in longitude, latitude or altitude fields, and obtains the cleaned data set Data-noNAN, Data-noNAN＝{(t _in ,long _in ,lat _in ,alt _in ,v _in ,vθ _in )|in=1,...,IN} (2) IN is the number of data points in Data-noNAN, 1≤IN≤I; 2.3: The data preprocessor performs format conversion on the data points in Data-noNAN: convert time to timestamp, convert longitude and latitude to double precision in radians, north latitude is positive, south latitude is negative; east longitude is positive , the west longitude is negative; convert the altitude into a double precision with the unit of meter, and calculate the altitude with sea level as 0; convert the velocity magnitude into a double precision with the unit of m/s; convert the velocity direction into a double with the unit of radian Accuracy, true north is 0, clockwise is positive; get the data set Data-format after format conversion; Data-format={((t _in ′,long _in ′,lat _in ′,alt _in ′,v _in ′,vθ _in ′))|in=1,…,IN} (3) t _in ′ is the moment after t _in conversion; long _in ′ is the longitude coordinate after long _in conversion; lat _in ′ is the latitude coordinate after lat _in conversion; alt _in ′ is the altitude after alt _in conversion; v _in ′ is The magnitude of the velocity after the conversion of v _in ; vθ _in ′ is the direction of the velocity after the conversion of vθ _in ; 2.4: Extract the data of the same period of time in [T ₁ , T ₂ ] for Data-format, the method is: Divide the sampling period [T ₁ , T ₂ ] into P time periods dt, obtain the data points in each time period from Data-format, compare the accuracy of each data point, take the precision dimension PN, and divide the data point with the largest PN As a feature point in the dt, the number of precision dimensions refers to the number of dimensions whose number of decimal places meets the agreed number of effective digits. The specific method is as follows: 2.4.1: Let the variable p=1; 2.4.1.1: Let there be L+1 data points in the pth time period dt _p of Data-format, and let these L+1 data points form the data point set DD _{p in dt p :} DD _p ＝{D _p0 ,D _p1 ,…,D _pL |st:max(t _p0 ,t _p1 ,…,t _pL )-min(t _p0 ,t _p1 ,…,t _pL )≤dt _p } (4 ) D _pl is the l+1th data point, and the data points placed in DD _p satisfy the maximum value max(t _p0 ,t _p1 ,…,t _pL ) of the L+1 data points in DD _p and L in DD _p The difference between the minimum value min(t _p0 ,t _p1 ,…,t _pL ) of +1 data point is within the range of dt _p ; 2.4.1.2: Calculate the precision dimension PN of D _p0 , D _p1 ,...,D _pL , and obtain the data precision dimension of D _p0 , D _p1 ,...,D _pl as

2.4.1.3: Found

The maximum value in , let the data point of the maximum value be D _p max: D _p max＝(dpt _p ,dplong _p ,dplat _p ,dpalt _p ,dpv _p ,dpvθ _p ) (5) Among them, dpt p , dplong _p , dplat _p , dpalt _p , dpv _p , and dpvθ _p are the _values of the six dimensions of the data points corresponding to the maximum number of precision dimensions in the p-th time period dt _p ; 2.4.2: Put D _p max in the trajectory data set Data-precision; 2.4.3: Let p=p+1; if p≤P, go to 2.4.1.1; if p>P, it means that the data extraction in the same period is completed, and the trajectory data set of the target within the time period [T ₁ , T ₂ ] is obtained Data-precision:

Where P is the number of time periods divided by the sampling period [T ₁ , T ₂ ], that is, the number of trajectory data in Data-precision; In the third step, the trajectory resampler reconstructs the trajectory with a given sampling rate for Data-precision based on the cubic spline curve difference algorithm, and obtains the resampled data set Data-sample at the sampling rate f. The method is: 3.1: Set the sampling rate f for [T ₁ , T ₂ ], set the sampling time interval t _f =1/f, and generate IS-1 time points at equal intervals, IS is a positive integer, which is the time set T _f Sample size, put these IS-1 equally spaced time points and T ₁ , T ₂ into the time set T _f of sampling rate f, and take T ₁ + (IS-1)×t _f , T ₂ for the last time Minimum value: T _f ＝{T ₁ ,T ₁ +1×t _f ,T ₁ +2×t _f ,…,min(T ₁ +(IS-1)×t _f ,T ₂ )} (7) 3.2: Use the cubic spline curve difference algorithm to calculate the longitude, latitude, height, velocity, and velocity direction corresponding to each moment in T _f to form a resampling data set at the sampling rate f. The specific method is: 3.2.1 Select data from Data-precision and combine them into two-dimensional data (X, Y), where X={dpt _p |p=1,...,P} is time data, Y={Y _k |k=1 ,...,5} contains other 5 dimension data;

3.2.2: let k=1; 3.2.3: For the data set (X,Y _k ), set P data points as (x ₀ ,y ₀ ),(x ₁ ,y ₁ ),…,(x _P-1 ,y _P-1 ), Use the cubic spline difference algorithm to solve the cubic spline coefficient vector (a _q ,b _q ,c _q ,d _q ) of each curve segment, so that in each interval x _q ≤x≤x _q+1 , q= 0,1,…,P-1, the spline function for data fitting is: f _q (x)＝(a _q ,b _q ,c _q ,d _q )×(1,(xx _q ),(xx _q ) ² ,(xx _q ) ³ ) ^T (9) 3.2.4: Bring the sampling time into formula (9), calculate the fitting result, and store the fitting result in the resampling data set Data-sample, the method is: 3.2.4.1: let is=1, 3.2.4.2 Take the is-th data of time set T _f , let it be t _is , t _is =T ₁ +(is-1)×t _f ; 3.2.4.3: Calculate the k-th dimension data at time t _is

3.2.4.4: Will

Save in the fitted matrix

On the is-th data of k row; 3.2.4.5: Let is=is+1: If is≤IS, the k-th dimension data has not been fitted and calculated, go to 3.2.4.2, If is>IS, the k-th dimension data fitting calculation is completed, go to 3.2.4.6; 3.2.4.6: Let k=k+1: If k≤5, go to 3.2.3; If k>5, all dimension data fitting calculations are completed, and a fitting matrix is generated

for:

Fitted data based on time data T _f and other 5 dimensions

Get the resampled data set Data-sample at the sampling rate f: Data-sample＝{(dst _is ,dslong _is ,dslat _is ,dsalt _is ,dsv _is ,dsvθ _is )|t _is ∈T _f ,is=1,…,IS} (12) Time dst _is = T _{f, is} , T _{f, is} represents the is data of T _f ; dslong _is ,dslat _is ,dsalt _is ,dsv _is ,dsvθ _is are the results of data fitting for each dimension, namely

In the fourth step, the speed corrector corrects the trajectory data velocity vector of the Data-sample based on equal time intervals to obtain the trajectory reconstruction data set D. The method is: 4.1: Calculate the space velocity of each data point according to the Data-sample, and obtain the space velocity resampling data set, the method is: 4.1.1: let is=1; 4.1.2 According to the is data point in Data-sample (dst _is ,dslong _is ,dslat _is ,dsalt _is ,dsv _is ,dsvθ _is ) time, longitude, latitude and height value calculation to get the space velocity

and space velocity direction

4.1.3 Will

into the space velocity resampling dataset

4.1.4 Let is=is+1; if is≤IS, go to 4.1.2; if is>IS, it means that the calculation is completed and the speed resampling data set is obtained

4.2: Correct the velocity vector according to the set error rate h, that is, use the calculated space velocity vector as a standard to correct the velocity vector collected in the Data-sample, if the collected velocity vector is within the range of the error rate, keep the collected velocity vector , otherwise replace the velocity vector collected in Data-sample with the calculated space velocity vector: 4.2.1 Receive the error rate h of user input from the keyboard; 4.2.2: let is=1; 4.2.3 Take the velocity magnitude dsv _is and the velocity direction dsvθ _is of the data point of the is item of Data-sample, and take

The space velocity of the data point of the is item in

and space velocity direction

Calculate dsv _is and

The gap, dsvθ _is and

According to the relationship between the gap and h, determine the magnitude of the velocity vector v _is ^* and the velocity direction vθ _is ^* after correction, the method is shown in formula (18):

4.2.4 Use (v _is ^* , vθ _is ^* ) to replace (dsv _is , dsvθ _is ) in the data set Data-sample, and replace (dst _is , dslong _is , dslat _is , dsalt _is , v _is ^* , vθ _is ^* ) into the trajectory reconstruction data set D; 4.2.5 Let is=is+1; if is≤IS, go to 4.2.2; if is>IS, it means that the replacement is completed, and the trajectory reconstruction data with the sampling rate f set in the target time period [T ₁ , T ₂ ] is obtained , including trajectory reconstruction datasets D and T _f :

2. a kind of trajectory reconstruction method based on situation target as claimed in claim 1, it is characterized in that the data format of t of the _ith sampling point data in Data-original is " year-month-day hour: minute: second "; long _i data format is "(E/W) degrees: minutes: seconds", (E/W) means east longitude and west longitude; lat _i data format is "(N/S) degrees: minutes: seconds" (N /S) indicates the north latitude and south latitude; alt _i is the altitude, the unit is meter; v _i indicates the speed of the target at a certain moment, the unit is km/h; vθ _i indicates the speed direction of the target at a certain moment, and the data format is "degrees :minute:second", true north is 0, and clockwise is positive. 3. a kind of trajectory reconstruction method based on situation target as claimed in claim 1, it is characterized in that the data preprocessor of step 2.2 carries out the method of data cleaning to Data-original is; 2.2.1: Let i=1; 2.2.2: Determine whether the i-th data point in the Data-original has missing longitude, latitude, or altitude fields. If any of these three items is missing, delete the i-th data point from the Data-original. If If none of these three items is missing, save the i-th data point to Data-noNAN; 2.2.3: Let i=i+1; if i≤I, go to 2.2.2; if i>I, it means that the elimination is completed, and the cleaned data set Data-noNAN is obtained. 4. a kind of trajectory reconstruction method based on situation target as claimed in claim 1, it is characterized in that the method that 2.3 step data preprocessors in 2.2 steps carry out format conversion to the data point in Data-noNAN is: 2.3.1: let in=1; 2.3.2: Use time conversion function ft, longitude conversion function flong, latitude conversion function flat, height conversion function fal, speed size conversion function fv, speed direction conversion function fvθ to the in-th data point in Data-noNAN (t _in , long _in ,lat _in ,alt _in ,v _in ,vθ _in ) are converted separately: Make the converted moment t _in '=ft (t _in ); ft (t _in ) means to use the java time tool class, the format of t _in : "year-month-day hour: minute: second" is converted to a timestamp; Make the converted longitude coordinate long _in ′=flong (long _in ), and flong (long _in ) means that the long _in data format of long in: "(E/W) degree: minute: second" is converted into a double-precision floating-point number, ( E/W) indicates east longitude and west longitude, the method is: take the prime meridian as 0, east longitude as positive, and west longitude as negative, determine whether long _in is positive or negative, if negative, add "-" before long _in to convert to radians After making, keep the precision of long _in to three decimal places; Make the converted latitude coordinate lat _in '=flat(lat _in ), flat(lat _in ) means to convert the latitude data format of lat _in : "(N/S) degree: minute: second" into a double-precision floating-point number, ( N/S) means north latitude and south latitude. The method is: take the equator as 0, north latitude as positive, and south latitude as negative. Determine whether lat _in is positive or negative. If it is negative, add "-" before lat _in and convert it to radians. And keep the accuracy of lat _in to three decimal places; Make the converted height alt _in '=fal(alt _in ), fal(alt _in ) means that the height alt _in is converted into the double precision with the unit of meter, the method is: take sea level as 0, above sea level as positive, sea level Below the plane is negative, determine whether alt _in is positive or negative, if negative, add "-" before alt _in , and keep the accuracy of alt _in to three decimal places; Make the converted velocity v _in ′=fv(v _in ), fv(v _in ) means converting the velocity v _in with the unit of km/h into a double-precision floating-point number, the unit is m/s, precision: keep the decimal point last three digits; Let the converted velocity direction vθ _in ′=fvθ(vθ _in ), fvθ(vθ _in ) means converting the "degree:minute:second" of the velocity direction _vθin into a double-precision floating-point number, the method is: the true north direction is 0 , clockwise is positive, counterclockwise is negative, determine whether vθ _in is positive or negative, if negative, add "-" before vθ _in , and keep the accuracy of alt _in to four decimal places. Save the converted in-th data point (ie (t _in ′, long _in ′, lat _in ′, alt _in ′, v _in ′, vθ _in ′) to the converted data Data-format; 2.3.3: Let i=i+1; if i≤IN, go to 2.3.2; if i>IN, it means the conversion is completed and Data-format is obtained. 5. a kind of trajectory reconstruction method based on situation target as claimed in claim 1, it is characterized in that described dt takes 60 seconds; Described sampling rate f takes 1, promptly collects once data every interval 1 second; Described lambda defaults is 1; the error rate h is 0.05. 6. a kind of trajectory reconstruction method based on situation target as claimed in claim 1, it is characterized in that the described calculation D _p0 of 2.4.1.2 step, D _p1 ,..., the method for the precision dimension number PN of D _pL is: 2.4.1.2.1 Let the variable l=0; 2.4.1.2.2 Check the number of effective digits of data precision of long _in ′, lat _in ′, alt _in ′, v _in ′, vθ _in ′ of D _pl , if there are w dimensions that satisfy the decimal point after format conversion in step 2.3 digits, then let the precision dimension of D _pl be

2.4.1.2.3 Let l=l+1; if l≤L, go to 2.4.1.2.2; if l>L, the calculation is completed, and the number of data precision dimensions of D _p0 , D _p1 ,...,D _pl is

7. a kind of trajectory reconstruction method based on situation target as claimed in claim 1, it is characterized in that the described in 4.1.2 step calculates the size of space velocity according to the is data point in Data-sample

and space velocity direction

The method is as follows: 4.1.2.1 Use the Haversine formula to calculate the distance between the is-th data point and the js-th data point (dst _js , dslong _js , dslat _js , dsalt _js , dsv _js , dsvθ _js ) on the horizontal plane according to the latitude and longitude, js=is+λ , λ is the step size, and the Haversin formula is:

Among them: haversin(θ)=sin ² (θ/2)=(1-cos(θ))/2, R is the radius of the earth, which is 6371km; The distance d _is,js on the horizontal plane between the data point of the is item and the data point of the js item is calculated as:

4.1.2.2 Calculation of space velocity

and space velocity direction

The method is:

a, b are the parameters of the calculation process, and the parameter a is cos(lat _is+1 )sin(long _is+1 -long _is ); The parameter b is

Images