CN113986906B - Track reconstruction method based on situation target - Google Patents
Track reconstruction method based on situation target Download PDFInfo
- Publication number
- CN113986906B CN113986906B CN202111422153.1A CN202111422153A CN113986906B CN 113986906 B CN113986906 B CN 113986906B CN 202111422153 A CN202111422153 A CN 202111422153A CN 113986906 B CN113986906 B CN 113986906B
- Authority
- CN
- China
- Prior art keywords
- data
- precision
- track
- point
- time
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
- G06F16/20—Information retrieval; Database structures therefor; File system structures therefor of structured data, e.g. relational data
- G06F16/21—Design, administration or maintenance of databases
- G06F16/215—Improving data quality; Data cleansing, e.g. de-duplication, removing invalid entries or correcting typographical errors
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
- G06F16/20—Information retrieval; Database structures therefor; File system structures therefor of structured data, e.g. relational data
- G06F16/25—Integrating or interfacing systems involving database management systems
- G06F16/254—Extract, transform and load [ETL] procedures, e.g. ETL data flows in data warehouses
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
- G06F16/20—Information retrieval; Database structures therefor; File system structures therefor of structured data, e.g. relational data
- G06F16/25—Integrating or interfacing systems involving database management systems
- G06F16/258—Data format conversion from or to a database
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
- G06F16/20—Information retrieval; Database structures therefor; File system structures therefor of structured data, e.g. relational data
- G06F16/29—Geographical information databases
Landscapes
- Engineering & Computer Science (AREA)
- Databases & Information Systems (AREA)
- Theoretical Computer Science (AREA)
- Data Mining & Analysis (AREA)
- Physics & Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Quality & Reliability (AREA)
- Remote Sensing (AREA)
- Complex Calculations (AREA)
Abstract
The invention discloses a situation target-based track reconstruction method, and aims to solve the problem that the existing method cannot process target track data with non-uniform track data point distribution. The technical scheme is as follows: preprocessing the situation target track original data, namely cleaning target track data points from different means, converting formats and extracting data in the same time period; then setting a sampling frequency to carry out trajectory reconstruction with a given sampling rate on the original trajectory data to obtain a resampling data set under the sampling rate f; and finally, carrying out track data speed vector correction on the resampled data set based on the equal time interval to obtain a track reconstruction data set. By adopting the method, the track reconstruction data set with uniformly distributed track data points under the user-defined sampling rate can be obtained, the noise at the dense part of the original data of the target track is reduced, the fitting error at the sparse part of the target track data is reduced, the track line of the target track data on the three-digit earth is clear, and the high precision of the velocity vector is ensured.
Description
Technical Field
The invention relates to the field of geographic information situation, in particular to a method for reconstructing target track data.
Background
The situation target track data monitoring means are numerous, the data sources are wide, and the data volume is large and the data complexity degree is high. The situation target trajectory data comprises target state data of the same moment monitored by different means, wherein the target state data comprises geographic space three-dimensional coordinates, target action speed and the like, and the target trajectory data of multidimensional space and time can be effectively displayed through visualization tools such as a three-dimensional digital earth and the like. In order to enhance the observation effect of the action track of the target and enhance the utilization value of the target track data, the target track data needs to be cleaned and reconstructed.
At present, the processing of track data acquired by multiple means and multiple acquisition frequencies is rough, the original data is generally directly displayed on a three-dimensional digital earth, and the movement behavior of a target is estimated by performing time series connection on dense points formed by observing the original data by a worker, wherein the method is limited by the density degree of the data and the experience of the observer: when the original trajectory point data is unevenly distributed (which means that a large number of trajectory data points exist in a certain period of time or a certain moment, and the number of data points in a certain period of time is too small or missing), the difficulty of estimating the motion behavior of the target trajectory by a worker is increased, and the estimation accuracy is difficult to control.
The existing track data cleaning and reconstructing methods mainly comprise data noise reduction and curve fitting, and due to the fact that the methods are single in data processing dimension and uniform in data distribution requirements, the results of processing a large amount of track data acquired by various means and various acquisition frequencies are poor, and particularly the problem that track data points are not uniformly distributed cannot be processed: for example, in a common filtering fitting method, when data points are unevenly distributed, a large amount of sawtooth noise is generated at a data dense part, and fitting errors at a data sparse part are large, which aggravates the sparseness of partial data distribution and enables a large amount of partial data to be overlapped, so that the reconstructed data cannot represent original data information when the conventional data cleaning and reconstructing method processes unevenly distributed data, and is not beneficial to analysis, mining and use of the data by other algorithms in a later period.
How to solve the problem that the existing track data cleaning and reconstructing method can not process the target track data with non-uniform track data point distribution is a technical problem which is of great concern to technicians in the field.
Disclosure of Invention
The invention aims to solve the technical problem of providing a situation target-based track data reconstruction method, solve the problem that the existing method cannot process target track data with non-uniform track data point distribution, reduce noise at the position of dense data and reduce fitting errors at the position of sparse data.
In order to solve the technical problems, the technical scheme of the invention is as follows: preprocessing situation target track original data, namely cleaning target track data points from different means, converting formats and extracting data in the same time period; then setting sampling frequency and carrying out track reconstruction with a given sampling rate on the original track data by using a curve interpolation algorithm to obtain a resampling data set under the sampling rate f; and finally, carrying out track data velocity vector correction on the resampled data set based on the equal time interval to obtain a track reconstruction data set D.
The invention comprises the following four steps:
the method comprises the following steps of firstly, constructing a situation target track data reconstruction system, wherein the situation target track data reconstruction system is composed of a data preprocessor, a track resampler and a speed corrector.
The Data preprocessor is connected with the track resampler, reads a situation target track original Data set Data-original from a target track original Data set, performs Data cleaning, data format conversion and Data extraction in the same time period on the Data-original to obtain a preprocessed Data set Data-precision, and sends the Data-precision to the track 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, performs trajectory reconstruction with a given sampling rate on the Data-precision by adopting a cubic spline curve difference algorithm to obtain a resampled Data set Data-sample, and sends the Data-sample to the speed corrector.
And the speed corrector is connected with the track resampler, receives the Data-sample from the track resampler, and performs track Data speed vector correction on the resampled Data point in the Data-sample to obtain a track reconstruction Data set D.
Secondly, the Data preprocessor preprocesses the Data-original Data set of the situation target track, including cleaning, format conversion and Data extraction in the same time period of the situation target track Data to obtain a target sampling period [ T [ ] 1 ,T 2 ]Data-precision of an inner trajectory Data set, the method is as follows:
2.1: the Data preprocessor reads a situation target track original Data set Data-original from a target track original Data set,
Data-original＝{(t i ,long i ,lat i ,alt i ,v i ,vθ i )|i＝1,…,I} (1)
targeting the Data-origin at the sampling period T 1 ,T 2 ](t i ∈[T 1 ,T 2 ]，T 1 Is the start of the sampling period, T, set by the user 2 Is the end point of a sampling period set by a user, the sampling refers to the process of sampling the action state of a situation target, and the sampling period refers to the time interval of sampling the action state of the situation target), a set consisting of I data points acquired by various means, wherein I is a positive integer, the ith sampling point data comprises 6 dimensions, and each dimension is respectively:
time t i (data format: year-month-day: minutes: seconds);
longitudinal coordinate long of data point i And (data format: "(E/W) degree: minute: second")
Latitude coordinate lat of data point i (data format: "(N/S) degree: minutes: seconds");
height alt i Altitude, in meters;
magnitude of velocity v i The unit of the speed of the target at a certain moment is km/h;
velocity direction v θ i The speed direction of the target at a certain moment is shown (data format: degree: minute: second), the north direction is 0, and the clockwise direction is positive);
2.2: the Data preprocessor cleans Data-original, eliminates Data with missing longitude, latitude or height field, and obtains cleaned Data set Data-nann;
2.2.1: let i =1;
2.2.2: judging whether the ith Data point in the Data-original has the deletion of a longitude field, a latitude field or an altitude field, if any of the 3 items is absent, deleting the ith Data point from the Data-original, and if none of the 3 items is absent, storing the ith Data point into a Data-nann;
2.2.3: let i = i +1; if I is less than or equal to I, rotating to 2.2.2; if I is larger than I, the elimination is finished, and a cleaned Data set Data-nann 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 the Data-nans, and IN is more than or equal to 1 and less than or equal to I.
2.3: the Data preprocessor performs format conversion on Data points in the Data-nana: converting time to a time stamp, converting longitude and latitude to double precision (north latitude is positive, south latitude is negative; east longitude is positive, west longitude is negative) in radian, converting altitude to double precision (sea level is 0) in meter, converting speed to double precision in m/s, and converting speed direction to double precision (north direction is 0, clockwise is positive) in radian, the method comprises the following steps:
2.3.1: let in =1;
2.3.2: using time conversion function ft, longitude conversion function floating, latitude conversion function flat, altitude conversion function fal, speed magnitude conversion function fv, speed direction conversion function fv theta to the in-th Data point (t) in Data-nann in ,long in ,lat in ,alt in ,v in ,vθ in ) Each item in (a) is converted separately:
let the converted time t in ′＝ft(t in )；ft(t in ) Indicating Using java time tool class, let t in The format of (c): "year-month-day time: minutes: seconds" is converted into a time stamp (i.e., the total number of seconds from greenwich time to current time));
let the converted longitude coordinate Long in ′＝flong(long in )，flong(long in ) Indicate will long in Longitude data format of (1): the method comprises the following steps of (E/W) degree, minute and second are converted into double-precision floating point numbers with radian as a unit (precision, three digits after decimal point is reserved, namely, the angle unit is minute and second is converted into international unified unit radian measure, and (E/W) represents east longitude and west longitude), wherein the method comprises the following steps: let the meridian of the first meridian be 0, the east meridian be positive, and the west meridian be negative, determine long in Whether it is positive or negative, and if it is negative at long in Adding "-" to the front, and turning to long after making radian in The three bits behind the decimal point are reserved for the precision of (2);
let the converted latitude coordinate lat in ′＝flat(lat in )，flat(lat in ) Indicate will lat in Latitude data format of (1): the unit of radian is converted into double-precision floating point number (precision: three bits after decimal point is reserved, namely angle unit is minute second and is converted into international unified unit radian system, and the method is as follows: let the equator be 0, north latitude be positive and south latitude be negative, determine lat in Whether it is positive or negative, if it is negative at lat in Adding "-" to the front, converting into radian system, and then adding lat in The three bits behind the decimal point are reserved for the precision of (2);
let the converted height alt in ′＝fal(alt in )，fal(alt in ) Indicates the height alt in Converting into double precision with the unit of meter, the method is as follows: let the sea level be 0, positive above sea level and negative below sea level, determine alt in Whether it is positive or negative, if it is negative at alt in Add "-" and let alt in The precision of (2) is reserved three bits after decimal point;
let the converted speed magnitude v in ′＝fv(v in )，fv(v in ) Indicating the magnitude v of the velocity in (unit is km/h) into double-precision floating point number (unit is m/s, precision: three bits after decimal point is reserved);
let the converted speed direction v theta in ′＝fvθ(vθ in )，fvθ(vθ in ) Indicates the velocity direction v theta in (data format: "degree: minutes: seconds") to a double-precision floating-point number (i.e. to international when the angle unit is minutes and seconds)Unified unit arc system) by the following method: the positive north direction is 0, the clockwise direction is positive, the anticlockwise direction is negative, and v theta is determined in Whether it is positive or negative, if it is negative at v θ in Add "-" and let alt in The four digits after the decimal point are reserved for the precision of (2);
the converted in data point (i.e., (t) in ′,long in ′,lat in ′,alt in ′,v in ′,vθ in ') save to converted Data-format;
2.3.3: let i = i +1; if i is less than or equal to IN, rotating to 2.3.2; if i is greater than IN, the conversion is finished, and a Data-format of the Data set 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: performing [ T ] on Data-format 1 ,T 2 ]And extracting data in the same time period. Because Data is acquired by multiple means, a large number of Data points are often acquired at the same time or in a short time range, and the Data points are different in precision (for example, the precision of time Data is different to seconds, and the precision of time Data is different to minutes, so that the decimal number after conversion is different), the Data-format is required to be in T 1 ,T 2 ]The data extraction is carried out in the same time period (the same time period is the same time or a shorter time range), and the method comprises the following steps:
will sample period [ T 1 ,T 2 ]Dividing into P time periods dt (dt generally takes 60 seconds), acquiring Data points in each time period from Data-format, comparing the precision of each Data point, taking precision dimension number PN (the dimension number of the decimal point number after the decimal point number meets 2.3 steps of format conversion is the precision dimension number), and taking the Data point with the maximum PN as a characteristic point in dt, wherein the specific method comprises the following steps:
2.4.1: let variable p =1;
2.4.1.1: let the p-th time period dt of Data-format p There are L +1 data points, let these L +1 data points constitute dt p Inner set of data points DD 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 For the l +1 th data point, the DD is put p Satisfies DD p Maximum value max (t) at time of middle L +1 data points p0 ,t p1 ,…,t pL ) And DD p Time minimum min (t) of middle L +1 data points p0 ,t p1 ,…,t pL ) Difference of (d) in dt p Within the range.
2.4.1.2: calculating D p0 ,D p1 ,…,D pL Precision dimension number PN:
(due to (t) in ′,long in ′,lat in ′,alt in ′,v in ′,vθ in ') the time dimension is an integer. No significant number of precision bits need to be calculated, so PN is at most 5), only long is calculated in ′,lat in ′,alt in ′,v in ′,vθ in ' precision dimension number of these 5 dimensions, the method is:
2.4.1.2.1 let variable l =0;
2.4.1.2.2 examination D pl Long of (2) in ′,lat in ′,alt in ′,v in ′,vθ in If w dimensions satisfy the decimal point number after 2.3 step format conversion, let D pl Degree of precision dimension of
2.4.1.2.3 let l = l +1; if L is less than or equal to L, rotating to 2.4.1.2.2; if l>L, obtaining D after the calculation is finished p0 ,D p1 ,…,D pl Has a data precision dimension of
Rotating for 2.4.1.3;
2.4.1.3: find out
Let the data point of the maximum be D p max：
D p max＝(dpt p ,dplong p ,dplat p ,dpalt p ,dpv p ,dpvθ p ) (5)
Wherein dpt p ,dplong p ,dplat p ,dpalt p ,dpv p ,dpvθ p Respectively the maximum value of the precision dimension number in the p time period dt p The values of the 6 dimensions of the inner corresponding data point.
2.4.2: will D p max is put into the trajectory Data set Data-precision;
2.4.3: let p = p +1; if P is less than or equal to P, rotating to 2.4.1.1; if p is>P, indicating that the data extraction is finished in the same time interval, and obtaining the target at [ T ] 1 ,T 2 ]Trajectory Data-precision within a time period:
where P is the sampling period [ T ] 1 ,T 2 ]The number of divided time periods, that is, the number of trace Data in Data-precision.
Thirdly, the track resampler carries out track reconstruction with a given sampling rate on Data-precision based on a cubic spline difference algorithm to obtain a resampled Data set Data-sample under the sampling rate f, and the method comprises the following steps:
3.1: to [ T ] 1 ,T 2 ]Setting the sampling rate f (the sampling rate is typically 1, i.e. data is acquired every 1 second) (let the acquisition time interval t be f = 1/f), the time of generating IS-1 equally spaced time points, IS being a positive integer and a time set T f IS-1 equally spaced time points and T 1 ,T 2 Set of moments T put at a sampling rate f f In, the last moment is T 1 +(IS-1)×t f ,T 2 Minimum value of (c):
T f ＝{T 1 ,T 1 +1×t f ,T 1 +2×t f ,…,min(T 1 +(IS-1)×t f ,T 2 )} (7)
3.2: the cubic spline curve difference algorithm is used (cubic spline interpolation algorithm see: [1 ]]Songyang, ACARS track reconstruction algorithm research based on adaptive cubic spline interpolation, china university of civil aviation, 2017.) f The longitude, latitude, altitude, speed and speed direction corresponding to each moment in the process form a resampling data set under the sampling rate f, and the specific method is as follows:
3.2.1 selecting Data from Data-precision, combining into two-dimensional Data (X, Y), where X = { dpt) p I P =1, \ 8230, P is time data, Y = { (Y) } k L k =1, \8230 |, 5} contains other 5 dimensional data;
3.2.2: let k =1;
3.2.3: for data set (X, Y) k ) Let P data points be (x) 0 ,y 0 ),(x 1 ,y 1 ),…,(x P-1 ,y P-1 ) Solving the coefficient vector (a) of cubic spline curve of each curve segment by using cubic spline curve difference algorithm q ,b q ,c q ,d q ) So that in each interval x q ≤x≤x q+1 Q =0,1, \ 8230;, P-1, the spline function to which the data is fitted is:
f q (x)＝(a q ,b q ,c q ,d q )×(1,(x-x q ),(x-x q ) 2 ,(x-x q ) 3 ) T (9)
3.2.4: substituting the sampling time into a formula (9), calculating a fitting result, and storing the fitting result into a resampled Data set Data-sample by the following method:
3.2.4.1: let is =1 and the ratio of the average power,
3.2.4.2 time set T f The ith data of (1), let as t is ，t is ＝T 1 +(is-1)×t f ；
3.2.4.3: calculating the time t is Of the k-th dimensionAccording to
3.2.4.4: will be provided with
Saved in the fitted matrix
The k lines of (1) are to the ith data;
3.2.4.5: let is = is +1:
if IS less than or equal to IS, the kth dimension data IS not fit and calculated, then 3.2.4.2 IS converted,
if IS greater than IS, the fitting calculation of the kth dimension data IS completed, and then the operation IS switched to 3.2.4.6;
3.2.4.6: let k = k +1:
if k is less than or equal to 5, rotating to 3.2.3;
if k is>5, fitting and calculating all dimension data to generate a fitting matrix
Comprises the following steps:
based on time data T f Data fitted with other 5 dimensions
Obtaining a resampling Data-sample under 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 T f Is data);
dslong is ,dslat is ,dsalt is ,dsv is ,dsvθ is the result of fitting for each dimension of data, i.e.
Fourthly, the velocity corrector carries out track Data velocity vector correction on the Data-sample based on the equal time interval to obtain a track reconstruction Data set D, and the method comprises the following steps:
the Data-sample is Data under equal time interval, the velocity vector is resampling of monitoring Data, in practice, the velocity vector can be calculated according to time, longitude, latitude and altitude, in order to not introduce calculation error, the existing method generally observes the velocity vector through the sensor, but the velocity vector acquired by the sensor has large volatility, can not be directly used for Data analysis, in order to ensure higher precision of the velocity vector, the invention uses the calculated space velocity vector to compare and correct with the velocity vector acquired by the sensor observation, and obtains better velocity vector Data, and the specific steps are as follows:
4.1: calculating the space velocity of each Data point according to the Data-sample to obtain a space velocity resampling Data set, wherein the method comprises the following steps:
4.1.1: let is =1;
4.1.2 Data-sample based on the is Data Point
(dst is ,dslong is ,dslat is ,dsalt is ,dsv is ,dsvθ is ) Calculating the time, longitude, latitude and height value to obtain the space velocity
And direction of space velocity
The method comprises the following steps:
4.1.2.1 Haverine's formula is used to calculate the is data point and the js data point (dst) according to longitude and latitude js ,dslong js ,dslat js ,dsalt js ,dsv js ,dsvθ js ) (js = is + λ, λ being the step size, default is 1) distance on the horizontal plane, the Haversin formula being:
wherein: haversin (theta) = sin 2 (theta/2) = (1-cos (theta))/2, R is the radius of the earth, and 6371km is taken;
calculating the distance d between the horizontal plane of the data point of the is item and the horizontal plane of the data point of the js item is,js Comprises the following steps:
4.1.2.2 calculation of spatial velocity magnitude
And direction of space velocity
The method comprises the following steps:
a and b are parameters of the calculation process, and the parameter a is cos (lat) is+1 )sin(long is+1 -long is )；
Parameter b is
4.1.3 will
Put into space velocity resample data set
4.1.4 order is = is +1; if IS less than or equal to IS, turning to 4.1.2; if IS > IS, the calculation IS completed, and a speed resampling data set IS obtained:
4.2: correcting the velocity vector according to a set error rate (η) by:
correcting the speed vector collected in the Data-sample by taking the calculated space speed vector as a standard, if the collected speed vector is in an error rate range, keeping the collected speed vector, and if not, replacing the speed vector collected in the Data-sample by the calculated space speed vector.
4.2.1 receiving the error rate h input by the user from the keyboard, generally 0.05;
4.2.2: let is =1;
4.2.3 speed size dsv of Data-sample Data point of item is is Direction of velocity dsv θ is Get it
Space velocity size of the middle is item data point
And direction of space velocity
Calculating dsv is And
difference, dsv θ of is And
according to the relation between the difference and h, the magnitude v of the corrected velocity vector is determined is * And a velocity direction v θ is * The method is shown as formula (18):
4.2.4 use (v) is *,vθ is * ) Replacement of (dsv) in Data-sample of Data set is ,dsvθ is ) Will (dst) is ,dslong is ,dslat is ,dsalt is ,v is *,vθ is * ) Into the trajectory reconstruction data set D.
4.2.5 order is = is +1; if IS less than or equal to IS, rotating to 4.2.2; if is>IS, indicating the completion of the replacement, obtaining a target time period [ T 1 ,T 2 ]Trajectory reconstruction data for setting a sampling rate f, comprising trajectory reconstruction data sets D and T f ：
The method is adopted to process the original track data, and a track reconstruction data set under a set sampling rate can be obtained. The invention can achieve the following technical effects:
1. by adopting the method and the device, the trajectory reconstruction data set with uniformly distributed trajectory data points under the user-defined sampling rate can be obtained, and the problem that the existing method cannot process target trajectory data with non-uniformly distributed trajectory data points is solved.
2. The second step data preprocessor of the invention adopts the same time interval data extraction method to filter a large amount of data points in the same time interval, thereby reducing the noise of the dense position of the original data of the target track.
3. In the third step, the trajectory resampler performs trajectory reconstruction by adopting a trajectory reconstruction method with a given sampling rate, target trajectory data points are uniformly distributed in a time dimension by using the set sampling rate, fitting errors at sparse positions of the target trajectory data are reduced, and the trajectory of the fitted target trajectory data on the three-position digital earth is clear.
4. The fourth step speed corrector of the invention uses the calculated space speed vector and the speed vector observed and collected by the sensor to carry out comparison correction based on equal time intervals, obtains better speed vector data and can ensure higher precision of the speed vector.
The track reconstruction data set obtained by the method can help researchers to better utilize the track data to analyze data, can help front-line workers to better observe the clear process of target activities, enhances the observation effect on the action track of the target and enhances the utilization value of the target track data.
Drawings
FIG. 1 is a general flow diagram of the present invention.
Fig. 2 is a logical structure diagram of a situation target trajectory reconstruction system constructed in the first step of the present invention.
Detailed Description
The invention is described in further detail below with reference to the drawings and the detailed description.
FIG. 1 is a general flow diagram of the present invention, comprising four steps.
The method comprises the following steps of firstly, constructing a situation target track data reconstruction system, wherein the situation target track data reconstruction system is composed of a data preprocessor, a track resampler and a speed corrector, and is shown in figure 2.
The Data preprocessor is connected with the track resampler, reads a situation target track original Data set Data-original from a target track original Data set, performs Data cleaning, data format conversion and Data extraction in the same time period on the Data-original to obtain a preprocessed Data set Data-precision, and sends the Data-precision to the track 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, performs trajectory reconstruction with a given sampling rate on the Data-precision by adopting a cubic spline curve difference algorithm to obtain a resampled Data set Data-sample, and sends the Data-sample to the speed corrector.
And the speed corrector is connected with the track resampler, receives the Data-sample from the track resampler, and performs track Data speed vector correction on the resampled Data point in the Data-sample to obtain a track reconstruction Data set D.
Secondly, preprocessing a Data-original Data set of the situation target track by a Data preprocessor, wherein the preprocessing comprises the steps of cleaning, format conversion and Data extraction in the same time period of the situation target track Data to obtain a target sampling period [ T [ ] 1 ,T 2 ]Data-precision of an inner trajectory Data set, the method is as follows:
2.1: the Data preprocessor reads a situation target track original Data set Data-original from the target track original Data set,
Data-original＝{(t i ,longi,lat i ,alt i ,v i ,vθ i )|i＝1,…,I} (2)
targeting the Data-origin at the sampling period T 1 ,T 2 ](t i ∈[T 1 ,T 2 ]，T 1 Is the start of the sampling period, T, set by the user 2 Is the end point of a sampling period set by a user, the sampling refers to the process of sampling the action state of a situation target, the sampling period refers to the time interval of sampling the action state of the situation target), a set consisting of I data points acquired by various means, I is a positive integer, the ith sampling point data comprises 6 dimensions, and each dimension respectively:
time t i (data format: year-month-day: minute: second, such as 2021-08-15;
longitude coordinate long of data point i And (data format: "(E/W) degree: minute: second")
Latitude coordinate lat of data point i (data format: "(N/S) degree: minutes: seconds", such as N23: 35;
height alt i Altitude, in meters;
magnitude of velocity v i The unit of the speed of the target at a certain moment is km/h;
direction of velocity v θ i Indicating the speed direction at a certain moment of the target (data format: degree: minute: second), the north-positive direction is 0, the clockwise direction is positive, such as 23;
2.2: the Data preprocessor cleans Data of Data-original, eliminates Data with missing longitude, latitude or height field, and obtains a cleaned Data set Data-nann;
2.2.1: let i =1;
2.2.2: judging whether the ith Data point in the Data-original has the deletion of a longitude field, a latitude field or an altitude field, if any of the 3 items is absent, deleting the ith Data point from the Data-original, and if none of the 3 items is absent, storing the ith Data point into a Data-nann;
2.2.3: let i = i +1; if I is less than or equal to I, rotating to 2.2.2; if I is larger than I, the elimination is finished, and a cleaned Data set Data-nann 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 the Data-nanan, and IN is more than or equal to 1 and less than or equal to I.
2.3: the Data preprocessor performs format conversion on Data points in the Data-nans: converting time to a time stamp, converting longitude and latitude to double precision (north latitude is positive, south latitude is negative; east longitude is positive, west longitude is negative) in radian, converting altitude to double precision (sea level is 0) in meter, converting speed to double precision in m/s, and converting speed direction to double precision (north direction is 0, clockwise is positive) in radian, the method comprises the following steps:
2.3.1: let in =1;
2.3.2: using time conversion function ft, longitude conversion function floating, latitude conversion function flat, altitude conversion function fal, speed magnitude conversion function fv, speed direction conversion function fv theta to the in-th Data point (t) in Data-nann in ,long in ,lat in ,alt in ,v in ,vθ in ) Each item in (a) is converted separately:
let the converted time t in ′＝ft(t in )；ft(t in ) Indicating Using the java time tool class, will t in The format of (1): "year-month-day: minute: second" is converted into a time stamp (total number of seconds from greenwich mean time to current time)), for example: "2021-10-01-12: 1633061532000;
let the converted longitude coordinate long in ′＝flong(long in )，flong(long in ) Indicate will long in Longitude data format of (1): the unit of radian is converted into double-precision floating point number (precision: three bits after decimal point is reserved, namely angle unit is minute second is converted into international unified unit radian system, (E/W) represents east longitude and west longitude), and the method is as follows: let the meridian of the first meridian be 0, the east meridian be positive, and the west meridian be negative, determine long in Whether it is positive or negative, and if it is negative at long in Adding "-" to the front, and turning to long after making radian in The accuracy of (D) is retained three decimal places, for example, longitude (east longitude: 116 deg. 23 '17') of Beijing Tiananmen square is converted to (2.030 rad);
let the converted latitude coordinate lat in ′＝flat(lat in )，flat(lat in ) Indicate will lat in Latitude data format of (1): the unit of radian is converted into double-precision floating point number (precision: three bits after decimal point is reserved, namely angle unit is minute second and is converted into international unified unit radian system, and the method is as follows: let the equator be 0, north latitude be positive and south latitude be negative, determine lat in Whether it is positive or negative, if it is negative at lat in Adding "-" to the front, converting into radian system, and then adding lat in The precision of (2) is reserved for three decimal places, for example, after conversion to (0.696 rad) from latitude (north latitude: 39 DEG 54 '27') of Tiananmen square of Beijing;
let the converted height alt in ′＝fal(alt in )，fal(alt in ) Indicates the height alt in Converting into double precision with the unit of meter, the method is as follows: the sea level is 0 and abovePositive, negative below sea level, determine alt in Whether it is positive or negative, if it is negative at alt in Add "-" and let alt in The precision of (d) is preserved three decimal places, e.g., height 122.455m translates to altitude 122.455m, depth 122.455m translates to altitude-122.455 m;
let the converted speed magnitude v in ′＝fv(v in )，fv(v in ) Indicating the magnitude v of the velocity in (in km/h) to a double precision floating point number (in m/s, precision: three bits after decimal point hold), for example 100.000m/s after 360.000km/h conversion;
let the converted speed direction v theta in ′＝fvθ(vθ in )，fvθ(vθ in ) Indicates the velocity direction v theta in (data format: degree: minutes: second) is converted into a double-precision floating point number (namely, the angle unit is divided into minutes and seconds and is converted into an international unified unit radian system), and the method comprises the following steps: the positive north direction is 0, the clockwise direction is positive, the anticlockwise direction is negative, and v theta is determined in Whether it is positive or negative, if it is negative at v θ in Add "-" and let alt in The precision of (c) is maintained four digits after the decimal point, for example, 23 degrees 35 minutes 15 seconds north to east is converted to 0.4115rad;
the converted in data point (i.e., (t) in ′,long in ′,lat in ′,alt in ′,v in ′,vθ in ') save to converted Data-format;
2.3.3: let i = i +1; if i is less than or equal to IN, rotating to 2.3.2; if i is greater than IN, the conversion is finished, and a Data-format of the Data set 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: performing [ T ] on Data-format 1 ,T 2 ]And extracting data in the same time period. Because the data is acquired by multiple means, a large number of data points exist at the same time or in a shorter time range, and the data points have different accuracies (for example, the accuracy of the time data is different from the second accuracy to the second accuracy, and the accuracy of the time data is different from the minute accuracy to the minute accuracy, which causes the decimal place after conversion to be different), the data are acquired by multiple meansThis is required for Data-format at [ T ] 1 ,T 2 ]The data extraction is carried out in the same time period (the same time period is the same time or a shorter time range), and the method comprises the following steps:
will sample period [ T 1 ,T 2 ]Dividing into P time periods dt (dt is 60 seconds generally), acquiring Data points in each time period from Data-format, comparing the precision of each Data point, and acquiring precision dimension number PN (the decimal point number satisfies the dimension of appointed effective digital number as the precision dimension number, such as Data (2.030, 0.696, 122.455, 100.000, 0.4115) of a timestamp 1633061532000, wherein the decimal point numbers of the Data with 5 dimensions all satisfy the effective digital number after 2.3-step format conversion, and the precision dimension number PN is 5); taking the data point with the maximum PN as the characteristic point in dt, the specific method is as follows:
2.4.1: let variable p =1;
2.4.1.1: let the p-th time period dt of Data-format p There are L +1 data points, let these L +1 data points constitute dt p Inner set of data points DD 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 For the l +1 th data point, the DD is put p Satisfies DD p Maximum value max (t) at time of middle L +1 data points p0 ,t p1 ,…,t pL ) And DD p Time minimum min (t) of middle L +1 data points p0 ,t p1 ,…,t pL ) Difference of (d) in dt p Within the range.
2.4.1.2: calculating D p0 ,D p1 ,…,D pL Precision dimension number PN:
(due to (t) in ′,long in ′,lat in ′,alt in ′,v in ′,vθ in ') the time dimension is an integer. No significant number of precision bits need to be calculated, so PN is at most 5), only long is calculated in ′,lat in ′,alt in ′,v in ′,vθ in ' precision dimension number of these 5 dimensions, the method is:
2.4.1.2.1 let variable l =0;
2.4.1.2.2 examination D pl Long of (2) in ′,lat in ′,alt in ′,v in ′,vθ in ' if there are w dimensions satisfying the decimal place number after 2.3 step format conversion, let D pl Degree of precision dimension of
2.4.1.2.3 let l = l +1; if L is less than or equal to L, rotating to 2.4.1.2.2; if l>L, after the calculation, D is obtained p0 ,D p1 ,…,D pl Has a data precision dimension of
Rotating for 2.4.1.3;
2.4.1.3: find out
Let the data point of the maximum be D p max：
D p max＝(dpt p ,dplong p ,dplat p ,dpalt p ,dpv p ,dpvθ p ) (5)
Wherein dpt p ,dplong p ,dplat p ,dpalt p ,dpv p ,dpvθ p Respectively the maximum value of the precision dimension number in the p time period dt p Values of 6 dimensions of the inner corresponding data point.
2.4.2: will D p max is put into the trajectory Data set Data-precision;
2.4.3: let p = p +1; if P is less than or equal to P, rotating to 2.4.1.1; if p is>P, indicating that the data extraction is finished in the same time interval, and obtaining the target at [ T ] 1 ,T 2 ]Trajectory Data set Data-precision over time period:
where P is the sampling period [ T ] 1 ,T 2 ]The number of divided time periods, that is, the number of trace Data in Data-precision.
Thirdly, the track resampler performs track reconstruction with a given sampling rate on Data-precision based on a cubic spline difference algorithm to obtain a resampled Data set Data-sample under a sampling rate f, and the method comprises the following steps:
3.1: to [ T ] 1 ,T 2 ]Setting the sampling rate f (the sampling rate is typically 1, i.e. data is acquired every 1 second) (let the acquisition time interval t be f = 1/f), the time of generating IS-1 equally spaced time points, IS being a positive integer and a time set T f IS-1 equally spaced time points and T 1 ,T 2 Set of moments T put at a sampling rate f f In, the last moment is T 1 +(IS-1)×t f ,T 2 Minimum value of (c):
T f ＝{T 1 ,T 1 +1×t f ,T 1 +2×t f ,…,min(T 1 +(IS-1)×t f ,T 2 )} (7)
3.2: the cubic spline curve difference algorithm is used (cubic spline interpolation algorithm see: [1 ]]Songyang, ACARS track reconstruction algorithm research based on adaptive cubic spline interpolation, china university of civil aviation, 2017.) f The longitude, latitude, altitude, speed and speed direction corresponding to each moment form a resampling data set under the sampling rate f, and the specific method comprises the following steps:
3.2.1 selecting Data from Data-precision, combining into two-dimensional Data (X, Y), where X = { dpt) p I P =1, \ 8230, P is time data, Y = { (Y) } k I k =1, \ 8230 |, 5} contains the other 5 dimensional data;
3.2.2: let k =1;
3.2.3: for data set (X, Y) k ) Let P data points be(x 0 ,y 0 ),(x 1 ,y 1 ),…,(x P-1 ,y P-1 ) Solving cubic spline coefficient vector (a) of each curve segment by using cubic spline curve difference algorithm q ,b q ,c q ,d q ) So that in each interval x q ≤x≤x q+1 Q =0,1, \8230;, P-1, the spline function to which the data was fitted is:
f q (x)＝(a q ,b q ,c q ,d q )×(1,(x-x q ),(x-x q ) 2 ,(x-x q ) 3 ) T (9)
3.2.4: substituting the sampling time into a formula (9), calculating a fitting result, and storing the fitting result into a resampled Data set Data-sample by the following method:
3.2.4.1: let is =1 and the ratio of the average power,
3.2.4.2 time set T f The is-th data of (1), let as t is ，t is ＝T 1 +(is-1)×t f ；
3.2.4.3: calculating the time t is K-th dimension data of (2)
3.2.4.4: will be provided with
Saved in the fitted matrix
K lines of (1) is th data;
3.2.4.5: let is = is +1:
if IS less than or equal to IS, the k-dimension data IS not fit for calculation, then 3.2.4.2 IS carried out,
if IS > IS, fitting and calculating the kth dimension data, and turning to 3.2.4.6;
3.2.4.6: let k = k +1:
if k is less than or equal to 5, rotating to 3.2.3;
if k is>5, fitting and calculating all dimension data to generate a fitting matrix
Comprises the following steps:
based on time data T f Data fitted with other 5 dimensions
Obtaining a resampling Data-sample under 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 T f Is data);
dslong is ,dslat is ,dsalt is ,dsv is ,dsvθ is the result of fitting for each dimension of data, i.e.
Fourthly, the velocity corrector carries out track Data velocity vector correction on the Data-sample based on the equal time interval to obtain a track reconstruction Data set D, and the method comprises the following steps:
the Data-sample is Data under equal time interval, the velocity vector is resampling of monitoring Data, in practice, the velocity vector can be calculated according to time, longitude, latitude and altitude, in order to not introduce calculation error, the existing method generally observes the velocity vector through the sensor, but the velocity vector acquired by the sensor has large volatility, can not be directly used for Data analysis, in order to ensure higher precision of the velocity vector, the invention uses the calculated space velocity vector to compare and correct with the velocity vector acquired by the sensor observation, and obtains better velocity vector Data, and the specific steps are as follows:
4.1: calculating the spatial velocity of each Data point according to the Data-sample to obtain a spatial velocity resampling Data set, wherein the method comprises the following steps:
4.1.1: let is =1;
4.1.2 Data-sample based on the is Data Point
(dst is ,dslong is ,dslat is ,dsalt is ,dsv is ,dsvθ is ) Calculating time, longitude, latitude and height values of the space to obtain space velocity
And direction of space velocity
The method comprises the following steps:
4.1.2.1 Haverine formula is used for calculating the data point of the is th item and the data point of the js th item (dst) according to the longitude and latitude js ,dslong js ,dslat js ,dsalt js ,dsv js ,dsvθ js ) (js = is + λ, λ is the step size, default is 1) the distance on the horizontal plane, the Haversin formula:
wherein: haversin (theta) = sin 2 (theta/2) = (1-cos (theta))/2, R is the radius of the earth, and 6371km is taken;
calculating the distance d between the horizontal plane of the data point of the is item and the horizontal plane of the data point of the js item is,js Comprises the following steps:
4.1.2.2 calculating the magnitude of the space velocity
And direction of space velocity
The method comprises the following steps:
a, b are parameters of the calculation process, and a is cos (lat) is+1 )sin(long is+1 -long is )；
The parameter b is
4.1.3 will
Put into space velocity resample data set
4.1.4 order is = is +1; if IS less than or equal to IS, rotating to 4.1.2; if IS > IS, the calculation IS finished, and a speed resampling data set IS obtained:
4.2: correcting the velocity vector according to a set error rate (h) by:
correcting the speed vector collected in the Data-sample by taking the calculated space speed vector as a standard, if the collected speed vector is in an error rate range, keeping the collected speed vector, and if not, replacing the speed vector collected in the Data-sample by the calculated space speed vector.
4.2.1 receiving the error rate h input by the user from the keyboard, generally 0.05;
4.2.2: let is =1;
4.2.3 get the velocity magnitude dsv of the Data-sample item Data point of the is is Direction of velocity dsv θ is Get it
Space velocity magnitude of the data point of the middle is item
And direction of space velocity
Calculating dsv is And
difference of (dsv θ) is And
according to the relation between the difference and h, the magnitude v of the corrected velocity vector is determined is * And a velocity direction v θ is * The method is shown in formula (18):
4.2.4 use (v) is *,vθ is * ) Replacement of (dsv) in Data-sample of Data set is ,dsvθ is ) Will (dst) is ,dslong is ,dslat is ,dsalt is ,v is *,vθ is * ) Is placed in the trajectory reconstruction data set D.
4.2.5 order is = is +1; if IS less than or equal to IS, rotating to 4.2.2; if is>IS, which indicates the completion of the replacement, obtains a target time period [ T ] 1 ,T 2 ]Number of track reconstructions for setting sampling rate fAccording to, including trajectory reconstruction data sets D and T f ：
Claims (7)
1. A track reconstruction method based on situation targets is characterized by comprising the following steps:
the method comprises the following steps of firstly, constructing a situation target track data reconstruction system, wherein the situation target track data reconstruction system consists of a data preprocessor, a track resampler and a speed corrector;
the Data preprocessor is connected with the track resampler, reads a situation target track original Data set Data-original from a target track original Data set, performs Data cleaning, data format conversion and Data extraction in the same time period on the Data-original to obtain a preprocessed Data set Data-precision, and sends the Data-precision to the track resampler;
the trajectory resampler is connected with the Data preprocessor and the speed corrector, receives the Data-precision of the preprocessed Data set from the Data preprocessor, carries out trajectory reconstruction with a given sampling rate on the Data-precision to obtain a resampled Data set Data-sample, and sends the Data-sample to the speed corrector;
the speed corrector is connected with the track resampler, receives the Data-sample from the track resampler, and performs track Data speed vector correction on the resampled Data point in the Data-sample to obtain a track reconstruction Data set D;
secondly, preprocessing a Data-original Data set of the situation target track by a Data preprocessor, wherein the preprocessing comprises the steps of cleaning, format conversion and Data extraction in the same time period of the situation target track Data to obtain a target sampling period [ T [ ] 1 ,T 2 ]The inner trajectory Data-precision, method is as follows:
2.1: the Data preprocessor reads a situation target track original Data set Data-original from a target track original Data set,
Data-original＝{(t i ,long i ,lat i ,alt i ,v i ,vθ i )|i＝1,…,I} (1)
targeting the Data-original at the sampling period [ T ] 1 ,T 2 ]In the method, a set consisting of I data points acquired by multiple means, wherein I is a positive integer, and t i ∈[T 1 ,T 2 ]，T 1 Is the start of the sampling period, T, set by the user 2 The method comprises the following steps that a sampling period is set by a user, sampling refers to the process of sampling the situation target action state, and the sampling period refers to the time interval of sampling the situation target action state; the ith sample point data contains 6 dimensions, each dimension is: time t i Longitude coordinate of data point long i Latitude coordinate of data point lat i Height alt i Magnitude of velocity v i (ii) a Direction of velocity v θ i ；
2.2: the Data preprocessor cleans the Data original, eliminates the Data with missing longitude, latitude or height field to obtain the cleaned Data set Data-nann,
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 the Data-nans, and IN is more than or equal to 1 and less than or equal to I;
2.3: the Data preprocessor performs format conversion on Data points in the Data-nana: converting the time into a timestamp, converting the longitude and the latitude into double precision with the unit of radian, wherein north latitude is positive and south latitude is negative; the east meridian is positive, and the west meridian is negative; converting the altitude into double precision with the unit of meter, and calculating the altitude by taking the sea level as 0; converting the speed into double precision with the unit of m/s; converting the speed direction into double precision with radian as a unit, wherein the north direction is 0, and the clockwise direction is positive; obtaining a Data-format of the Data set after format conversion;
Data-format＝{((t in ′,long in ′,lat in ′,alt in ′,v in ′,vθ in ′))|in＝1,…,IN} (3)
t in is' t in The time after the conversion; long in Is long in The converted longitude coordinates; lat in Is' lat in The latitude coordinate after conversion; alt in ' is alt in A converted height; v. of in ' is v in The converted speed; v θ in Is' v θ in The converted speed direction;
2.4: performing [ T ] on Data-format 1 ,T 2 ]Data extraction in the same time period is carried out, and the method comprises the following steps:
will sample period [ T 1 ,T 2 ]Dividing into P time periods dt, acquiring Data points in each time period from Data-format, comparing the precision of each Data point, taking a precision dimension number PN, taking the Data point with the maximum PN as a characteristic point in the dt, wherein the precision dimension number refers to the number of dimensions of which the decimal point number meets the agreed effective digital digit number, and the specific method comprises the following steps:
2.4.1: let variable p =1;
2.4.1.1: let the p-th time period dt of Data-format p There are L +1 data points, let these L +1 data points constitute dt p Inner set of data points DD 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 For the l +1 th data point, the DD is put p Satisfies DD p Maximum value max (t) at time of middle L +1 data points p0 ,t p1 ,…,t pL ) And DD p Time minimum min (t) of middle L +1 data points p0 ,t p1 ,…,t pL ) Difference of (d) in dt p Within the range;
2.4.1.2: calculating D p0 ,D p1 ,…,D pL To the precision dimension number PN of D to obtain D p0 ,D p1 ,…,D pl Has a data precision dimension of
2.4.1.3: find out
Let the data point of the maximum be D p max：
D p max＝(dpt p ,dplong p ,dplat p ,dpalt p ,dpv p ,dpvθ p ) (5)
Wherein dpt p ,dplong p ,dplat p ,dpalt p ,dpv p ,dpvθ p Respectively the maximum value of the precision dimension number in the p time period dt p Values of 6 dimensions of the inner corresponding data point;
2.4.2: will D p max is put into the trajectory Data set Data-precision;
2.4.3: let p = p +1; if P is less than or equal to P, rotating to 2.4.1.1; if p is>P, indicating that the data extraction is finished in the same time interval, and obtaining the target at [ T ] 1 ,T 2 ]Trajectory Data set Data-precision over time period:
where P is the sampling period [ T ] 1 ,T 2 ]The number of divided time periods, namely the number of the trace Data in the Data-precision;
thirdly, the track resampler carries out track reconstruction with a given sampling rate on Data-precision based on a cubic spline difference algorithm to obtain a resampled Data set Data-sample under the sampling rate f, and the method comprises the following steps:
3.1: to [ T ] 1 ,T 2 ]Setting a sampling rate f to make the collection time interval t f =1/f, the time of IS-1 equally spaced time points IS generated, IS IS a positive integer and IS a time set T f IS-1 equally spaced time points and T 1 ,T 2 Set of moments T put at the sampling rate f f In, the last moment is T 1 +(IS-1)×t f ,T 2 Minimum value of (d):
T f ＝{T 1 ,T 1 +1×t f ,T 1 +2×t f ,…,min(T 1 +(IS-1)×t f ,T 2 )} (7)
3.2: computing T using cubic spline curve difference algorithm f The longitude, latitude, altitude, speed and speed direction corresponding to each moment in the process form a resampling data set under the sampling rate f, and the specific method is as follows:
3.2.1 selecting Data from Data-precision, combining into two-dimensional Data (X, Y), where X = { dpt) p I P =1, \ 8230, P is time data, Y = { (Y) } k I k =1, \ 8230 |, 5} contains the other 5 dimensional data;
3.2.2: let k =1;
3.2.3: for the data set (X, Y) k ) Let P data points be (x) 0 ,y 0 ),(x 1 ,y 1 ),…,(x P-1 ,y P-1 ) Solving the coefficient vector (a) of cubic spline curve of each curve segment by using cubic spline curve difference algorithm q ,b q ,c q ,d q ) So that in each interval x q ≤x≤x q+1 Q =0,1, \8230;, P-1, the spline function to which the data was fitted is:
f q (x)＝(a q ,b q ,c q ,d q )×(1,(x-x q ),(x-x q ) 2 ,(x-x q ) 3 ) T (9)
3.2.4: substituting the sampling time into a formula (9), calculating a fitting result, and storing the fitting result into a resampled Data set Data-sample by the following method:
3.2.4.1: let is =1 and the ratio of the average power,
3.2.4.2 time set T f The ith data of (1), let as t is ，t is ＝T 1 +(is-1)×t f ；
3.2.4.3: calculating the time t is Data of k-th dimension
3.2.4.4: will be provided with
Saved in the fitted matrix
K lines of (1) is th data;
3.2.4.5: let is = is +1:
if IS less than or equal to IS, the k-dimension data IS not fit for calculation, then 3.2.4.2 IS carried out,
if IS greater than IS, the fitting calculation of the kth dimension data IS completed, and then the operation IS switched to 3.2.4.6;
3.2.4.6: let k = k +1:
if k is less than or equal to 5, rotating to 3.2.3;
if k is>5, fitting and calculating all dimension data to generate a fitting matrix
Comprises the following steps:
based on time data T f Data fitted with other 5 dimensions
Obtaining a Data-sample of a resampling Data set under a 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 T f The is-th data of (1);
dslong is ,dslat is ,dsalt is ,dsv is ,dsvθ is the result of fitting for each dimension of data, i.e.
Fourthly, the speed corrector carries out track Data speed vector correction on the Data-sample based on equal time intervals to obtain a track reconstruction Data set D, and the method comprises the following steps:
4.1: calculating the space velocity of each Data point according to the Data-sample to obtain a space velocity resampling Data set, wherein the method comprises the following steps:
4.1.1: let is =1;
4.1.2 Data-sample based on the is Data Point
(dst is ,dslong is ,dslat is ,dsalt is ,dsv is ,dsvθ is ) Calculating time, longitude, latitude and height values of the space to obtain space velocity
And direction of space velocity
4.1.3 will
Put into space velocity resample data set
4.1.4 order is = is +1; if IS less than or equal to IS, rotating to 4.1.2; if is>IS, descriptionAfter the calculation is finished, a speed resampling data set is obtained
4.2: correcting the speed vector according to the set error rate h, namely correcting the speed vector collected in the Data-sample by taking the calculated space speed vector as a standard, if the collected speed vector is in the error rate range, keeping the collected speed vector, and otherwise, replacing the speed vector collected in the Data-sample by the calculated space speed vector:
4.2.1 receiving an error rate h input by a user from a keyboard;
4.2.2: let is =1;
4.2.3 get the velocity magnitude dsv of the Data-sample item Data point of the is is Direction of velocity dsv θ is Get it
Space velocity magnitude of the data point of the middle is item
And direction of space velocity
Calculating dsv is And
difference of (dsv θ) is And
according to the relation between the difference and h, the magnitude v of the corrected velocity vector is determined is * And a velocity direction v θ is * The method is shown in formula (18):
4.2.4 use (v) is * ,vθ is * ) Replacement of (dsv) in Data-sample of Data set is ,dsvθ is ) Will (dst) is ,dslong is ,dslat is ,dsalt is ,v is * ,vθ is * ) Putting the data into a track reconstruction data set D;
4.2.5 order is = is +1; if IS less than or equal to IS, rotating to 4.2.2; if is>IS, indicating the completion of the replacement, obtaining a target time period [ T 1 ,T 2 ]Trajectory reconstruction data set for setting a sampling rate f, comprising trajectory reconstruction data sets D and T f ：
2. The method as claimed in claim 1, wherein t is the Data of the ith sampling point in the Data-original i The data format of the data is 'year-month-day: minute: second'; long i The data format is "(E/W) degree: minutes: seconds", (E/W) represents east longitude and west longitude; lat i The data format is "(N/S) degree: minutes: seconds" (N/S) represents north latitude and south latitude; alt i Altitude, in meters; v. of i The unit of the speed of the target at a certain moment is km/h; v θ i The speed direction of the target at a certain moment is shown, the data format is 'degree: minutes: seconds', the north direction is 0, and the clockwise direction is positive.
3. The posture target-based track reconstruction method as claimed in claim 1, wherein the Data-original is cleaned by the Data preprocessor in 2.2 steps;
2.2.1: let i =1;
2.2.2: judging whether the ith Data point in the Data-original has the deletion of a longitude field, a latitude field or an altitude field, if any of the 3 items is absent, deleting the ith Data point from the Data-original, and if none of the 3 items is absent, storing the ith Data point into a Data-nann;
2.2.3: let i = i +1; if I is less than or equal to I, rotating to 2.2.2; and if I is larger than I, the elimination is finished, and a cleaned Data set Data-nann is obtained.
4. The method for reconstructing a trajectory based on a situational target according to claim 1, wherein the 2.2-step 2.3-step Data preprocessor performs format conversion on Data points in Data-nans by:
2.3.1: let in =1;
2.3.2: using time conversion function ft, longitude conversion function floating, latitude conversion function flat, altitude conversion function fal, speed magnitude conversion function fv, speed direction conversion function fv theta to the in-th Data point (t) in Data-nann in ,long in ,lat in ,alt in ,v in ,vθ in ) Each of which is converted separately:
let the converted time t in ′＝ft(t in )；ft(t in ) Indicating Using the java time tool class, will t in The format of (c): converting the time stamp into a year-month-day time minute: second;
let the converted longitude coordinate long in ′＝flong(long in )，flong(long in ) Indicate will long in Longitude data format of (1): "(E/W) degree: minutes: seconds" into a double precision floating point number, (E/W) representing east and west meridians by: let the meridian of the first meridian be 0, the east meridian be positive, and the west meridian be negative, determine long in Whether it is positive or negative, and if it is negative at long in Adding minus before, converting into arc system, and then long in The precision of (2) is reserved three bits after decimal point;
let the converted latitude coordinate lat in ′＝flat(lat in )，flat(lat in ) Indicate will lat in Latitude data format of (1): the degree (N/S) is minute: second, and the system is converted into double-precision floating point number, and the degree (N/S) represents north latitude and south latitude by the following method: at the equator0, positive north latitude and negative south latitude, and determining lat in Whether it is positive or negative, if it is negative at lat in Adding minus before, converting into radian system, and adding lat in The precision of (2) is reserved three bits after decimal point;
let the converted height alt in ′＝fal(alt in )，fal(alt in ) Indicates the height alt in Converting into double precision with the unit of meter, and the method comprises the following steps: let sea level be 0, above sea level be positive, below sea level be negative, determine alt in Whether it is positive or negative, if it is negative at alt in Add "-" and let alt in The precision of (2) is reserved three bits after decimal point;
let the converted speed magnitude v in ′＝fv(v in )，fv(v in ) Representing the magnitude v of the velocity in km/h in Conversion to double precision floating point number, unit is m/s, precision: three bits after the decimal point are reserved;
let the converted speed direction v theta in ′＝fvθ(vθ in )，fvθ(vθ in ) Indicates the velocity direction v theta in Minutes to seconds is converted into a double-precision floating point number, and the method comprises the following steps: the positive north direction is 0, the clockwise direction is positive, the anticlockwise direction is negative, and v theta is determined in Whether it is positive or negative, if it is negative at v θ in Add "-" and let alt in The four decimal places are retained.
The converted in data point (i.e., (t) in ′,long in ′,lat in ′,alt in ′,v in ′,vθ in ') save to converted Data-format;
2.3.3: let i = i +1; if i is less than or equal to IN, rotating to 2.3.2; if i is greater than IN, the conversion is completed, and Data-format is obtained.
5. The method for reconstructing a trajectory based on situational targets according to claim 1, wherein dt is 60 seconds; the sampling rate f is 1, namely data are collected every 1 second; the lambda is 1 by default; the error rate h is taken to be 0.05.
6. The method of claim 1A track reconstruction method based on situation targets is characterized in that the step D is calculated in 2.4.1.2 p0 ,D p1 ,…,D pL The method for measuring the degree PN by the precision comprises the following steps:
2.4.1.2.1 let variable l =0;
2.4.1.2.2 examination D pl Long of in ′,lat in ′,alt in ′,v in ′,vθ in ' if there are w dimensions satisfying the decimal place number after 2.3 step format conversion, let D pl Degree of precision dimension of
2.4.1.2.3 let l = l +1; if L is less than or equal to L, rotating to 2.4.1.2.2; if l>L, after the calculation, D is obtained p0 ,D p1 ,…,D pl Has a data precision dimension of
7. The method as claimed in claim 1, wherein the step 4.1.2 of calculating the space velocity according to the is Data point in the Data-sample
And direction of space velocity
The method comprises the following steps:
4.1.2.1 Haverine's formula is used to calculate the is data point and the js data point (dst) according to longitude and latitude js ,dslong js ,dslat js ,dsalt js ,dsv js ,dsvθ js ) Distance on the horizontal plane, js = is + λ, λ is the step size, and the Haversin formula is:
wherein: haversin (theta) = sin 2 (theta/2) = (1-cos (theta))/2, R is the radius of the earth, and 6371km is taken;
calculating the distance d between the horizontal plane of the data point of the is item and the horizontal plane of the data point of the js item is,js Comprises the following steps:
4.1.2.2 calculating the magnitude of the space velocity
And direction of space velocity
The method comprises the following steps:
a and b are parameters of the calculation process, and the parameter a is cos (lat) is+1 )sin(long is+1 -long is )；
Parameter b is
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111422153.1A CN113986906B (en) | 2021-11-26 | 2021-11-26 | Track reconstruction method based on situation target |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111422153.1A CN113986906B (en) | 2021-11-26 | 2021-11-26 | Track reconstruction method based on situation target |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113986906A CN113986906A (en) | 2022-01-28 |
| CN113986906B true CN113986906B (en) | 2022-11-04 |
Family
ID=79732117
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202111422153.1A Active CN113986906B (en) | 2021-11-26 | 2021-11-26 | Track reconstruction method based on situation target |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113986906B (en) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117576294B (en) * | 2024-01-16 | 2024-03-26 | 太原理工大学 | Data transmission parameter optimization method of digital twin system of mining equipment |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107450472B (en) * | 2017-08-31 | 2019-06-07 | 华中科技大学 | The method for realizing cutter path parameter arc length based on three bezier curve interpolation |
| CN107450471B (en) * | 2017-08-31 | 2019-06-11 | 华中科技大学 | A method of cutter path parameter arc length is realized based on cubic PH curve interpolation |
-
2021
- 2021-11-26 CN CN202111422153.1A patent/CN113986906B/en active Active
Also Published As
| Publication number | Publication date |
|---|---|
| CN113986906A (en) | 2022-01-28 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109816664B (en) | Three-dimensional point cloud segmentation method and device | |
| CN109472802B (en) | Surface mesh model construction method based on edge feature self-constraint | |
| CN113986906B (en) | Track reconstruction method based on situation target | |
| CN111650579B (en) | InSAR mining area three-dimensional deformation estimation method and device for rock migration parameter adaptive acquisition and medium | |
| CN103077559B (en) | Based on the fruit ear three-dimensional rebuilding method of sequence image | |
| Van Niekerk et al. | Input variable selection for interpolating high-resolution climate surfaces for the Western Cape | |
| CN110038291B (en) | Motion information acquisition system and motion posture recovery method for rowing athletes | |
| CN115859116A (en) | Marine environment field reconstruction method based on radial basis function regression interpolation method | |
| CN115423696B (en) | Remote sensing orthographic image parallel generation method of self-adaptive thread parameters | |
| CN113156082B (en) | Method and system for identifying depth of active layer of crop root system | |
| Neyman et al. | Areal spread of the effect of cloud seeding at the Whitetop experiment | |
| CN115953085A (en) | Method for evaluating influence of composite dry heat event on vegetation growth | |
| CN116152512A (en) | Height measuring and calculating method based on building shadow restoration | |
| CN116126980A (en) | Object-oriented multi-scale landform type dividing method | |
| CN112380984B (en) | Remote sensing-based salt-biogas vegetation slow-flow capacity space evaluation method | |
| CN114418911A (en) | Method for carrying out statistical regression downscaling on remote sensing image and improving water body definition | |
| CN117649611B (en) | DEM data production processing method based on twice orientation | |
| CN117315002B (en) | Method, system and device for measuring three-dimensional phenotype of object to be measured | |
| CN114882100B (en) | Land coverage area estimation method and system based on sub-pixel mapping | |
| CN118364663B (en) | Tidal flat terrain remote sensing inversion method and system based on tide level complementary cumulative distribution function | |
| CN112700227B (en) | Method and system for calculating standard working hours of aerial survey interior work image amplitude and land class processing and application | |
| CN117670802A (en) | Ground extraction and correction method based on seed points | |
| CN114626184A (en) | Distributed photovoltaic total radiation estimation method and system based on satellite remote sensing | |
| CN118332026A (en) | Visualization method and device for ocean primary productivity distribution and computer equipment | |
| CN118363007A (en) | Tidal flat terrain time sequence monitoring method and system based on tidal flat remote sensing inundation frequency correction |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |