Azimuth observation data real-time processing method and system with backtracking restoration function

Abstract

The invention provides an azimuth observation data real-time processing method and system with a backtracking and repairing function, and belongs to the field of data processing. The problems that phase delay exists during real-time smoothing and outlier processing of underwater target observation data, and spot type outlier detection and correction processing are improper are solved. According to azimuth and course observation data received in real time, the azimuth observation data are smoothed in real time by adopting a local polynomial regression method based on a Gaussian kernel, and outliers existing in the azimuth data and the course data are identified in real time by adopting an innovation chi-square detection method based on Kalman filtering and are eliminated. And according to the existing non-outlier effective data prediction, correcting the azimuth data at the outlier, and carrying out backtracking restoration on the interval of the outlier data. The method is small in phase delay and good in data real-time performance; the Gaussian kernel property is combined with the past and future data secondary correction of the outlier interval, the accuracy of target position calculation can be improved, the accumulative error of the calculation algorithm is reduced, and the method has good theoretical and engineering application values.

Description

Technical Field

[0001] The present invention relates to the field of data processing technology, and in particular to a method and system for real-time processing of azimuth observation data with a backtracking repair function. Background Art

[0002] In the process of azimuth estimation, due to the influence of various factors such as environmental noise interference, sensor inaccuracy, and improper operation of experimenters, the received azimuth observation data contains high-frequency noise and abnormal values ​​that seriously deviate from the normal azimuth change law, namely, wild values.

[0003] High-frequency noise caused by environmental noise and sensor factors can cause small fluctuations in the azimuth data compared to the true value, which in turn leads to accumulated errors when using the observed data to solve the target motion elements. Common methods for removing noise include: median filtering, Gaussian filtering, sliding average filtering, exponential moving average filtering, Savitzky-Golay filtering, etc. In practical applications, it is hoped that the process of smoothing the observed data can be carried out in real time. Although most filtering methods can suppress high-frequency noise well, there is a large phase delay when applied to online estimation, and the smoothing effect is not good at the beginning and end of the data.

[0004] A large number of studies on methods for dealing with outliers have been carried out in China. Outliers are generally divided into isolated outliers and spot-type outliers. In experiments, if the operation is improper or affected by the environment, there will also be "step-type outliers" (cliff-type outliers). General methods for detecting outliers include: Wright criterion method, detection method based on Kalman innovation change, etc. Most outlier detection methods do not produce ideal results when detecting the start and end points of spot-type outliers. And if outliers are directly eliminated, it will often lead to data loss, and the data of the outlier part, especially the data interval of spot-type outliers, needs to be repaired. In "Polynomial Moving Smoothing Algorithm and Its Application in Flight Test Data Preprocessing", Yang Yajun et al. used a polynomial fitting method to estimate the data within the interval of the sliding window, and detected the outliers by judging whether the observed data at the next moment differed by 3 to 5 times the standard deviation, and used a polynomial to correct the outliers for the predicted value at the next moment; in "A Study on Outlier Removal and Smoothing Method of Flight Parameter Data", An Li et al. used k times the standard deviation to judge the outliers, considered the processing method of step-type outliers, and finally used Newton interpolation to correct the outlier data; in "Research on Key Technologies of Multi-Sensor Information Fusion", Kang Jian used a method for detecting outliers based on new information changes, and compensated the data of outlier points through a weighted function; in "A Real-time Outlier Removal Method for Azimuth Observation Data", Qian Jin et al. compared the residuals of the heading data at different times with 3 times the standard deviation σ for uniformly moving targets, recorded the data greater than 3σ as outliers, and used the characteristics of uniform motion of the target to correct the outlier data.

[0005] Based on the analysis of the above domestic and foreign literature, the current research on real-time smoothing and outlier removal of target observation data in underwater environments has problems such as phase delay in smoothing methods, improper detection of spot-type outliers, and improper correction processing. These problems will eventually lead to large deviations from the actual data. Summary of the invention

[0006] The technical problems to be solved by the present invention are:

[0007] In order to solve the problems of phase delay, improper detection and correction of spot-type outliers in existing underwater target observation data during real-time smoothing and outlier processing.

[0008] The present invention adopts the following technical solutions to solve the above technical problems:

[0009] The present invention provides a real-time processing method for azimuth observation data with a backtracking repair function, comprising the following steps:

[0010] S100, real-time reading of the current t n Heading data y HD 、Direction data y BR If the current shot is the first 5 shots, skip the data processing process and convert the current shot's heading data y HD 、Direction data y BR Store in the corresponding array, if the current beat is not the previous 5 beats of data, proceed to the next step;

[0011] S200, reading an identifier indicating whether the current shot data is valid; if the current shot data is a valid identifier, proceeding to the next step; if the current shot data is an invalid identifier, using a local polynomial fitting based on a Gaussian kernel according to the stored data to predict the data of the current shot;

[0012] S300, assuming that the current noise is Gaussian noise, and determining whether the position and heading data of the shot are outliers through chi-square detection;

[0013] S400, if it is determined to be an outlier, then the outlier type is determined in combination with the existing data, and the outlier types include spot-type outliers, step-type outliers and azimuthal spot-type outliers, and the determined outliers are repaired using the local fitting polynomial based on the Gaussian kernel in step S200;

[0014] S500, according to step S200 and step S300, combining the existing data n Smoothing the data;

[0015] S600, judging whether there are spot-type wild values ​​that can be retroactively repaired based on the existing data, and if so, performing retroactive repair;

[0016] S700, record the current shot's azimuth and heading data and the processed azimuth and heading data, continue to receive the next shot's data, return to step S100, and start the data processing flow for a new shot.

[0017] Further, in step S200, the fitting polynomial includes:

[0018] Assume the window width is N, and the data at each position within the window width is recorded as t i , the current shot data is recorded as t n , the time series within the window width is represented by t, t i It is expressed as:

[0019] t i ∈t={t i |0<t n-1 -t i ≤N} (1)

[0020] The Gaussian kernel is expressed as:

[0021]

[0022] Among them, y i Represents the observed azimuth and heading data;

[0023] Assume sample set Set the basis function Let polynomial for:

[0024]

[0025] Where k represents the degree of each factor in the polynomial, and M represents the highest degree of the polynomial;

[0026] The extraction coefficient matrix is ​​θ = [a 0 a 1 …a M ] T , the basis function matrix is The matrix form is:

[0027]

[0028] Let Φ(t)=[φ(t 0 )φ(t 1 )…φ(t n )] T , assuming that the regression function under ideal conditions is F(x), and selecting the second norm l 2 As a weighted loss function, let J be:

[0029]

[0030] Among them, w i represents weight, F(x)=[f(x 0 )f(x 1 )…f(x n )] T , the weight matrix is ​​calculated by the Gaussian kernel:

[0031]

[0032] Derivative of J(θ):

[0033]

[0034] θ=(Φ T WΦ) -1 Φ T WF (9)

[0035] In formula (3), when M=1, the model is a linear regression model; the polynomial fitted based on the existing valid data is used to predict the data of the current shot, and the predicted value is:

[0036] φ(t n )=θ T φ(t n ) (10)

[0037] Use the predicted value φ(t n ) replaces the invalid position and heading data of the current shot.

[0038] Furthermore, in step S300, it specifically includes:

[0039] Current shot data t n By n The data before the data is taken Y(t n-1 )The polynomial prediction value obtained by fitting is The actual observed data is y i ,

[0040] Calculate the Kalman filter innovation, the tth n The new value of the data is υ(t n )as follows:

[0041]

[0042] Among them, y(t n ) represents the tth n The actual observation data of the shooting data;

[0043] Assume that receiving n When taking data, the sequence of the first N-1 new messages is as follows:

[0044] I(tn-1 )={υ(t n-N ),υ(t n-N+1 ),...,υ(t n-1 )} (12)

[0045] Its mean is μ I (t 1:n-1 ), with variance σ I (t 1:n-1 ), the chi-square value is calculated as:

[0046]

[0047] The probability density function of the chi-square distribution is:

[0048]

[0049] The cumulative distribution function of the chi-square distribution is:

[0050]

[0051] Among them, Γ(·) is the gamma function, γ(·) is the incomplete gamma function, and d is the degree of freedom;

[0052] The significance level is α = 0.05, and the critical value is 3.84. 2 >3.84, then the t n The data is an outlier.

[0053] Furthermore, in step S400, the method for determining the outlier type includes:

[0054] If there are continuous outliers not less than l SPK If the value is σ, it is considered to be a spot type outlier: For the azimuth observation data, find the first time the data deviates from 1 standard deviation σ SPK Position t SPK , take the position as the starting point of the azimuth speckle type wild value, use the fitting polynomial of the position to predict the data in the speckle type wild value interval shot by shot, and replace the observed data; after that, the observed data is detected to be no more than 1 times the standard deviation σ SPK The position of the spot type wild value is taken as the actual end position; for the heading observation data, look forward to find the first time that the data deviates from 1 standard deviation σ SPK Position t SPK , taking this position as the starting point of the heading maneuver, and then detecting that the observed data is no greater than 1 times the standard deviation σ SPK The position of is taken as the end position of heading maneuver;

[0055] When the n The data is an outlier, if |υ(t n )-μI (t 1:n-1 )|>4σ(t 1:n-1 ) and duration Δt u Less than l speck beat, then it is considered that the t n The data is a step-type outlier; look forward to find the first time the data deviates from 1 times the standard deviation σ speck1 Position t speck1 , take this position as the starting point of the azimuth step wild value, use the polynomial fitted before this position to predict the data in the step wild value interval shot by shot, and replace the observed data; after that, the observed data is detected to be no more than 1 times the standard deviation σ speck1 The position of is taken as the actual end position of the step wild value.

[0056] Furthermore, in step S500, it specifically includes:

[0057] Assume the window width is N, and the data at each position within the window width is recorded as t, which can be expressed as:

[0058] t i ∈t={t i |0<t n -t i ≤N} (16)

[0059] The Gaussian kernel is expressed as:

[0060]

[0061] Other formulas are shown in (3)-(10). n The result of smoothing the data is:

[0062] Furthermore, in step S600, it specifically includes:

[0063] Let the current time be the tth n Shot data, assuming the last shot of the spot to be traced back is Let the set be Other outliers in the range belong to the set S OTH , to satisfy Backtracking of wild values;

[0064] Assume that the center of the spot type outlier is The length of the spot type field value is 2r, and the start of the spot type field value is Shoot, end The wild value interval of the spot type is:

[0065]

[0066] The length of the wild-value interval is 2r, the total window width is 2N+2r, and the interval π is selected in the complete azimuth data set. REV :

[0067]

[0068] According to formula (17), the Gaussian kernel W is obtained:

[0069]

[0070] In the entire orientation data set, excluding S SPK The outlier interval is used to perform local polynomial regression on the data before and after the outlier interval.

[0071] The local polynomial regression method is used to back-repair each value in the wild value interval of length 2r; suppose that the sth value in the interval i If the data is shot for backtracking repair, the corresponding Gaussian kernel interval is as follows:

[0072]

[0073] The Gaussian kernel is expressed as:

[0074]

[0075] Get a series of weight values:

[0076]

[0077] The Gaussian kernel weight matrix is ​​expressed as:

[0078]

[0079] The sth value in the outlier interval is obtained by formula (8)-(10) i The backtracking repair value of the shot data.

[0080] The present invention provides a real-time processing system for azimuth observation data with a backtracking and repairing function. The system has a program module corresponding to the above steps, and executes the steps in the real-time processing method for azimuth observation data with a backtracking and repairing function when running.

[0081] The present invention provides a computer-readable storage medium, wherein the computer-readable storage medium stores a computer program, and the computer program is configured to implement the steps of a real-time processing method for azimuth observation data with a backtracking repair function when called by a processor.

[0082] Compared with the prior art, the present invention has the following beneficial effects:

[0083] The present invention provides a real-time processing method and system for azimuth observation data with a backtracking repair function. The experimental data can be smoothed and processed in real time by setting relevant parameters according to the duration of the wild values ​​in the statistical data. The present invention uses an improved Gaussian kernel for local polynomial regression, with a small phase delay and good data real-time performance. The present invention can detect isolated wild values, spot wild values, and step wild values, and can obtain the starting and ending points of spot wild values ​​and step wild values, and can identify and correct wild value data in real time. The present invention provides a "backtracking repair" data processing method, which uses the properties of the Gaussian kernel and combines the past data and future data of the wild value interval for secondary correction. This method can improve the accuracy of target position solution, reduce the cumulative error of the solution algorithm, and has good theoretical and engineering application value. BRIEF DESCRIPTION OF THE DRAWINGS

[0084] Figure 1 It is a flow chart of a method for real-time processing of azimuth observation data with a backtracking repair function in an embodiment of the present invention;

[0085] Figure 2 The curve diagram of the heading angle and azimuth angle observation data containing isolated, spot-type and step-type wild values ​​in the embodiment of the present invention;

[0086] Figure 3 This is a diagram of outlier detection results in an embodiment of the present invention;

[0087] Figure 4 The distribution diagram of isolated, spot, and step wild values ​​in the embodiment of the present invention;

[0088] Figure 5 This is a data distribution diagram before real-time smoothing of the position data in an embodiment of the present invention;

[0089] Figure 6 This is a data distribution diagram after backtracking and repairing the step-type wild value interval in an embodiment of the present invention;

[0090] Figure 7 This is a real-time smoothing effect diagram of each shot of data in an embodiment of the present invention;

[0091] Figure 8 This is a diagram showing the effect of backtracking repair of spot-type wild values ​​in an embodiment of the present invention;

[0092] Fig. 9 This is a comparison diagram of the effects of the Savitzky-golay filtering method in the embodiment of the present invention and the method of the present invention after real-time smoothing of the azimuth data;

[0093] Fig.10 This is a comparison diagram of the effects of the sliding average filtering method in the embodiment of the present invention and the method of the present invention after real-time smoothing of the azimuth data;

[0094] Fig.11 This is a diagram showing the real-time processing result of data when it is assumed that the invalid interval is the 224th to 404th beats in an embodiment of the present invention. DETAILED DESCRIPTION

[0095] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, specific embodiments of the present invention are described in detail below with reference to the accompanying drawings.

[0096] Since the moving target's heading maneuver and environmental interference will cause outliers in the observed data, it is necessary to identify and correct isolated, spot, and step outliers. Figure 4 As shown in , when the effective observation data cannot be received due to environmental interference, Fig.11 As shown, the orientation data can be estimated online; then, the orientation observation data is smoothed online in real time according to the orientation and heading observation data received in real time, and finally, the data in the wild value interval is backtracked and repaired.

[0097] Specific implementation plan 1: Combine Figures 1 to 8 As shown, the present invention provides a real-time processing method for azimuth observation data with a backtracking repair function, comprising the following steps:

[0098] S100, read the current t n Heading data y HD 、Direction data y BR . In order to ensure the accuracy of the outlier recognition effect, it is necessary to calculate the standard deviation of the normal state observation data based on the experimental data in advance. The present invention uses the data of the first 5 shots of each experiment as the training data for calculating the standard deviation of the normal state observation data. If the current shot is the data of the first 5 shots of the experiment, the data processing process is skipped directly, and the heading data and azimuth data of the current shot are stored in the corresponding array. If the current shot is not the data of the first 5 shots of the experiment, proceed to the next step;

[0099] S200, read the flag indicating whether the current beat data is valid. If there is packet loss in the observed data, the data of this beat is invalid. The external interface should give a representation of invalid data in this beat. If the current beat data is a valid flag, proceed to the next step; if the current beat data is an invalid flag, fit a polynomial based on the stored data to predict the data of this beat;

[0100] The data fitting polynomial includes: assuming the window width is N, the data at each position within the window width is recorded as t i , the time series within the window width is expressed as t and is expressed as:

[0101]

[0102] The Gaussian kernel is expressed as:

[0103]

[0104] Assume sample set y i Represents the observed azimuth and heading data, and sets the basis function Let polynomial for:

[0105]

[0106] Where k represents the degree of each factor in the polynomial, and M represents the highest degree of the polynomial;

[0107] The extraction coefficient matrix is ​​θ = [a 0 a 1 …a M ] T , the basis function matrix is The matrix form is:

[0108]

[0109] Let Φ(t)=[φ(t 0 )φ(t 1 )…φ(t n )] T , assuming that the regression function under ideal conditions is F(x), and selecting the second norm l 2 (square error) is used as the weighted loss function and is set to J:

[0110]

[0111] Among them, w i represents weight, F(x)=[f(x 0 )f(x 1 )…f(x n )] T , the weight matrix is ​​calculated by the Gaussian kernel:

[0112]

[0113] To make F(x) and Φ(x) as close as possible within the window width, the weighted loss function J(θ) must be minimized, and the derivative of J(θ) is:

[0114]

[0115] θ=(Φ T WΦ) -1 Φ T WF (9)

[0116] In formula (3), when M=1, the model is a linear regression model, which conforms to the angle variation law of the present invention;

[0117] The polynomial fitted based on the existing valid data is used to predict the current beat data. The current beat is t n , the predicted value is:

[0118]

[0119] Use the predicted value φ(t n ) Replace the invalid position and heading data of the current shot;

[0120] S300, assuming that the noise in the experiment is Gaussian noise, and judging whether the position and heading data of the shot are outliers by chi-square test; assuming that the current shot is t n , the polynomial obtained by fitting the existing valid data is used to predict the data of this shot, as shown in formula (10), and the predicted value is

[0121] The current shot data is t n , t n The data is taken by t n The observed data Y(t n-1 )The polynomial prediction value obtained by fitting is The actual observed data is y i ,

[0122] Calculate the Kalman filter innovation, the tth n The new value of the data is υ(t n ),as follows:

[0123]

[0124] Among them, y(t n ) represents the tth n The actual observation data of the shooting data;

[0125] Assume that receiving n When taking data, the sequence of the first N-1 new messages is as follows:

[0126] I(t n-1 )={υ(t n-N ),υ(t n-N+1 ),...,υ(t n-1 )} (12)

[0127] Its mean is μ I (t 1:i-1 ), with variance σ I (t 1:i-1 ), the chi-square value is calculated as:

[0128]

[0129] The probability density function of the chi-square distribution is:

[0130]

[0131] The cumulative distribution function of the chi-square distribution is:

[0132]

[0133] Where Γ(·) is the gamma function and γ(·) is the incomplete gamma function. The degree of freedom d is 1, and the significance level is α = 0.05. At this time, the critical value is 3.84. If there is χ 2 >3.84, then the t n The data is an outlier;

[0134] S400, judging the wild value type based on the existing data and performing data repair, the judging method is:

[0135] If there are continuous outliers not less than l SPK If the data is not taken, it is considered to be a spot type outlier. For azimuth observation data, it is necessary to find the first time that the data deviates from 1 times the standard deviation σ SPK Position t SPK , take this position as the starting point of the azimuth speckle type wild value, use the polynomial fitted before this position to predict the data in the speckle type wild value interval shot by shot, and replace the observed data; after that, the observed data is detected to be no more than 1 times the standard deviation σ SPK The position of the spot type wild value is taken as the actual end position; for the heading observation data, the first deviation of the data from 1 standard deviation σ is also found forward SPK Position t SPK , taking this position as the starting point of the heading maneuver, and then detecting that the observed data is no greater than 1 times the standard deviation σ SPK The position of is taken as the end position of heading maneuver;

[0136] If the x-beat data detected continuously after the first beat is not an outlier, if the number of consecutive beats exceeds the minimum length l allowed for interruption of the spot-type outlier DLY , then the spot type wild value is considered to end, and the data change from this beat to the last beat is less than 1 times the standard deviation σ SPK The position of the spot type outlier is regarded as the end position of the spot type outlier. If the total length of the outlier is less than the minimum length l for identifying as a spot type outlier SPK , it is not treated as a spot-type outlier, but only as an isolated outlier, and is corrected using the fitted polynomial extrapolation result;

[0137] When the n The data is an outlier, if |υ(t n )-μ I (t 1:n-1 )|>4σ(t1:n-1 ) and duration Δt u Less than l speck beat, then it is considered that the t n The data is a step-type wild value; the characteristics of a step-type wild value are short duration, large jump amplitude, and will not return to the original data change trend; find the first time the data deviates from 1 times the standard deviation σ speck1 Position t speck1 , take this position as the starting point of the step-type wild value, use the polynomial fitted before this position to predict the data in the step-type wild value interval shot by shot, and replace the observed data; after that, the observed data is detected to be no greater than 1 times the standard deviation σ speck1 The position of is taken as the actual end position of the step-type wild value;

[0138] S500, according to step S200 and step S300, combining the existing data n Smoothing the data, assuming the window width is N, the data at each position within the window width is recorded as t, expressed as:

[0139] t i ∈t={t i |0<t n -t i ≤N} (16)

[0140] The Gaussian kernel is expressed as:

[0141]

[0142] Other formulas are shown in (3)-(10). n The result of smoothing the data is:

[0143] S600: Determine whether there are spot-type wild values ​​that can be retroactively repaired based on the existing data. If so, perform retroactive repair, specifically including:

[0144] Let the current time be the tth n Shot data, assuming the last shot of the spot to be traced back is Set Other outliers in the range belong to the set S OTH , only satisfied The condition can be met to backtrack the wild value;

[0145] Assume that the center of the spot type outlier is The length of the spot type field value is 2r, and the start of the spot type field value is Shoot, end The wild value interval of the spot type is:

[0146]

[0147] The length of the wild-value interval is 2r, the total window width is 2N+2r, and the interval π is selected in the complete azimuth data set. REV :

[0148]

[0149] According to formula (17), the Gaussian kernel W is obtained:

[0150]

[0151] In the entire orientation data set, excluding S SPK The outlier interval is used to perform local polynomial regression on the data before and after the outlier interval.

[0152] The local polynomial regression method is used to backtrack and repair each value in the wild value interval of length 2r; suppose that the i-th shot data in the interval is backtracked and repaired, then the corresponding Gaussian kernel interval is as follows:

[0153]

[0154] The Gaussian kernel is expressed as:

[0155]

[0156] Get a series of weight values:

[0157]

[0158] The Gaussian kernel weight matrix is ​​expressed as:

[0159]

[0160] The backtracking repair value of the i-th shot data in the wild value interval is obtained by formulas (8)-(10);

[0161] S700, record the current shot's azimuth and heading data and the processed azimuth and heading data, continue to receive the next shot's data, return to step S100, and start the data processing flow for a new shot.

[0162] Specific implementation scheme 2: The present invention provides a real-time processing system for azimuth observation data with a backtracking and repairing function. The system has a program module corresponding to the above steps, and executes the steps in the above-mentioned real-time processing method for azimuth observation data with a backtracking and repairing function during operation.

[0163] The other combinations and connection relationships of this embodiment are the same as those of the first embodiment.

[0164] Specific implementation scheme three: The present invention provides a computer-readable storage medium, which stores a computer program, and the computer program is configured to implement the steps of a real-time processing method for azimuth observation data with a backtracking repair function when called by a processor.

[0165] The other combinations and connection relationships of this embodiment are the same as those of the first embodiment.

[0166] Simulation experiment

[0167] Receive the target’s azimuth and heading observation data in real time, recorded as y BR and HD , record the real-time smoothed data as y RCO , the repaired output data is recorded as y RCV The test data used in this embodiment is as follows Figure 2 As shown, the target has undergone heading maneuvers at 95-180 beats and 712-800 beats, and the azimuth angle in this range is reflected as a spot-type wild value, the azimuth angle data at 272-313 beats is a spot-type wild value, and 419-432 is a step-type wild value.

[0168] The present invention uses a local polynomial regression method based on a Gaussian kernel to smooth the observed data. The window width is 80, and the current time is t n If the current shot is the data of the first 5 shots of the experiment, the data processing process is skipped directly, and the heading data and azimuth data of the current shot are stored in the corresponding array.

[0169] Let the data at each position within the window width be t i , the window width is N, expressed as:

[0170] t i ∈t={t i |||t n -t i ||≤N} (28)

[0171] The Gaussian kernel is expressed as:

[0172]

[0173] In order to achieve online smoothing and avoid large phase delays, such as Figure 5 As shown, only t i Conforms to the following formula:

[0174] t i ∈t={t i |0≤(t n -t i )≤N} (30)

[0175] Assume sample set Set the basis function Let polynomial as follows:

[0176]

[0177] The extraction coefficient matrix is ​​θ = [a 0 a 1 …a n ] T , the basis function matrix is The matrix form is:

[0178]

[0179] Φ(t)=[φ(t 0 )φ(t 1 )…φ(t n )] T .

[0180] Assume that the regression function under ideal conditions is F(x), and select the second norm l 2 (square error) is used as the weighted loss function and is set to J:

[0181]

[0182] where F(x) = [f(x 0 )f(x 1 )…f(x n )] T , the weight matrix is ​​calculated by the Gaussian kernel:

[0183]

[0184] To make F(x) and Φ(x) as close as possible within the window width, the weighted loss function J(θ) must be minimized, and the derivative of J(θ) is:

[0185]

[0186] θ=(Φ T WΦ) -1 Φ T WF (37)

[0187] In formula (31), when M=1, the model is a linear regression model, which conforms to the angle variation law of the present invention. The present invention uses a first-order polynomial to smooth the data, and the smoothing value is:

[0188]

[0189] The present invention adopts chi-square test to determine whether the position and heading data of the shot are outliers. Assume that the current shot is t n, the polynomial fitted based on the existing valid data is used to predict the current data, as shown in (10), and the predicted value is Assume that the data per shot within the window width is t i , t i The data is taken by t i Previous data Y(t i-1 )The polynomial prediction value obtained by fitting is The actual observed data is y i

[0190] Calculate the Kalman filter innovation, the tth n The new value of the data is υ(t n ),as follows:

[0191]

[0192] Assume that receiving i There is a front l when taking data win -1 beat new information sequence is as follows:

[0193] I(t n-1 )={υ(t n-N ),υ(t n-N+1 ),...,υ(t n-1 )} (40)

[0194] Its mean is μ I (t 1:n-1 ), with variance σ I (t 1:n-1 ), the chi-square value is calculated as:

[0195]

[0196] The degree of freedom is 1, and the significance level is α = 0.05. At this time, the critical value is 3.84. If there is a χ 2 >3.84, then the t n The data is an outlier.

[0197] The polynomial fitted at the previous moment is used to predict the time t i The prediction data. Let the coefficient matrix θ corresponding to the basis function matrix T , the predicted value is:

[0198]

[0199] In the following test, if {t 1 ,t 2 ,...,t 10} n If there is l near the data SPKIn order to find the start and end range of the spot type wild value more accurately, we can find the first time that the data change exceeds 1 times the standard deviation σ from the first beat of the spot type wild value. SPK This beat is used as the actual starting position of the spot-type wild value, and the predicted value is obtained through fitting polynomial prediction as shown in formula (44), which replaces the orientation data.

[0200] If the x-beat data detected continuously after the first beat is not an outlier, if the number of consecutive beats exceeds the minimum length l allowed for interruption of the spot-type outlier DLY , then the spot type wild value is considered to end, and the data change from this beat to the last beat is less than 1 times the standard deviation σ SPK The position of the spot type outlier is regarded as the end position of the spot type outlier. If the total length of the outlier is less than the minimum length l for identifying as a spot type outlier SPK , it is not treated as a spot-type outlier, but only as an isolated outlier, and is corrected using the fitted polynomial extrapolation result.

[0201] If the first-shot field value has |υ(t n )-μ I (t 1:n-1 )|>4σ(t 1:n-1 ), it is regarded as a step-type wild value, and the step of looking forward for the starting point of the spot-type wild value is performed, and the azimuth data is corrected in real time as shown in formula (44).

[0202] At the same time, the heading data is read to determine in real time whether a heading maneuver has occurred. The determination method is as follows (39)-(43). If a spot-type wild value is detected, it is considered that the target has performed a heading maneuver within the wild value interval, and the azimuth data of this section is regarded as a spot-type wild value. DLY After capturing valid data with non-wild values, backtracking repair is performed.

[0203] If the detection data is invalid at a certain moment, the data predicted in (44) is used as the orientation data at that moment.

[0204] Let the current time be the tth n Shot data, taking the step-type wild value near the 426th shot data in this embodiment as an example, assuming that the last shot of the spot to be traced back is Set Other outliers in the range belong to the set S OTH , only full The condition can be met to trace back the outlier. Let the center of the step-type outlier be The length of the step-type wild value is 2r, and the step-type wild value starts at the Shoot, end The step-type wild value interval is:

[0205]

[0206] The length of the wild-value interval is 2r, the total window width is 2N+2r, and the interval π is selected in the complete azimuth data set. REV :

[0207]

[0208] According to formula (17), the Gaussian kernel W is obtained:

[0209]

[0210] In the entire orientation data set, excluding S SPK The outlier interval is used to perform local polynomial regression on the data before and after the outlier interval.

[0211] Using the local polynomial regression method, each value in the wild value interval of length 2r is backtracked and repaired. Assuming that the i-th shot data in the interval is backtracked and repaired, the corresponding Gaussian kernel interval is as follows:

[0212]

[0213] The Gaussian kernel is expressed as:

[0214]

[0215] Get a series of weight values:

[0216]

[0217] The Gaussian kernel weight matrix is ​​expressed as:

[0218]

[0219] The backtracking repair value of the i-th beat data in the wild value interval is obtained by (8)-(10). The schematic diagram of the step-type wild value backtracking repair near 426 beats in this embodiment is shown in Figure 6 It can be seen that the data after backtracking repair has a better effect in smoothing the step-type wild values.

[0220] In this embodiment, when the “backtracking repair” method is not used to correct the data, the data obtained is as follows: Figure 7 As shown in the figure, after using the "backtracking repair" method to correct the data, the repaired data obtained is as follows Figure 8 As shown in the figure, it can be seen that the curve after backtracking repair has no burrs, is relatively smooth overall, and retains certain characteristic information of the original data. At the same time, the Savitzky-golay smoothing filter and recursive moving average filter method are designed to perform real-time smoothing on the test data used in this embodiment. The results are shown in the figure. Fig. 9 , Fig.10 As shown, it can be seen that the effect of smoothing using other methods is not good, especially not suitable for the processing of spot-type wild values ​​and economical wild values ​​near 272 and 426 beats. The root mean square error RMSE, signal-to-noise ratio SNR, smoothness R and other indicators are selected to evaluate the effect of real-time smoothing of experimental data by different methods. The root mean square error RMSE is calculated according to the following formula:

[0221]

[0222] Where: f(i) is the original signal, is the denoised signal, and n is the signal length.

[0223] The root mean square error reflects the difference between the original signal and the denoised signal. In actual use, the smaller the root mean square error, the better the denoising effect.

[0224] The signal-to-noise ratio SNR is calculated as follows:

[0225] SNR = 10 × log 10 (power signal / power noise ) (59)

[0226] in: power signal is the power of the original signal, power noise is the power of the noise, f(i) is the original signal, is the signal after denoising, and n is the signal length. It is generally believed that the higher the signal-to-noise ratio, the better the denoising effect.

[0227] The smoothness R is calculated according to the following formula:

[0228]

[0229] Where: f(i) is the original signal, is the denoised signal, and n is the signal length.

[0230] Smoothness is an important indicator for judging the effect of abnormal data processing results. It is generally believed that the smoother the signal, the smaller the value of the smoothness index, and the better the denoising effect.

[0231] The indicators of real-time smoothing data of different algorithms are shown in the following table:

[0232] Table 1 Indicators of different real-time smoothing methods

[0233]

[0234] It can be seen from Table 1 that the effect of real-time data smoothing of this method is better than the other two methods.

[0235] In this embodiment, a test experiment for invalid data is conducted. Assuming that the invalid data interval is 224 to 404 beats, the processing of invalid data in the present invention is based on the predicted value obtained by extrapolating the polynomial fitted by the least squares method as the observed data of the beat. The test results are shown in Fig.11 It can be seen that the prediction result given by the prediction method involved in the present invention is relatively accurate and does not deviate from the changing trend of the original data.

[0236] Although the present invention is disclosed as above, the protection scope of the present invention is not limited thereto. Those skilled in the art may make various changes and modifications without departing from the spirit and scope of the present invention, and these changes and modifications will fall within the protection scope of the present invention.

Claims

1. A real-time processing method for azimuth observation data with a backtracking repair function, characterized in that: The following steps are involved: S100, real-time reading of the current t n Heading data y HD 、Direction data y BR If the current shot is the first 5 shots, skip the data processing process and convert the current shot's heading data y HD 、Direction data y BR Store in the corresponding array, if the current beat is not the previous 5 beats of data, proceed to the next step; S200, reading an identifier indicating whether the current shot data is valid; if the current shot data is a valid identifier, proceeding to the next step; if the current shot data is an invalid identifier, using a local polynomial fitting based on a Gaussian kernel according to the stored data to predict the data of the current shot; S300, assuming that the current noise is Gaussian noise, and determining whether the position and heading data of the shot are outliers through chi-square detection; S400, if it is determined to be an outlier, then the outlier type is determined in combination with the existing data, and the outlier types include spot-type outliers, step-type outliers and azimuthal spot-type outliers, and the determined outliers are repaired using the local fitting polynomial based on the Gaussian kernel in step S200; S500, according to step S200 and step S300, combining the existing data n Smoothing the data; S600, judging whether there are spot-type wild values ​​that can be retroactively repaired based on the existing data, and if so, performing retroactive repair; S700, record the current shot's azimuth and heading data and the processed azimuth and heading data, continue to receive the next shot's data, return to step S100, and start the data processing flow for a new shot.

2. The method for real-time processing of azimuth observation data with a backtracking repair function according to claim 1, characterized in that: In step S200, the fitting polynomial includes: Assume the window width is N, and the data at each position within the window width is recorded as t i , the current shot data is recorded as t n , the time series within the window width is represented by t, t i It is expressed as: t i ∈t={t i |0<t n-1 -t i ≤N} (1) The Gaussian kernel is expressed as: Among them, y i Represents the observed azimuth and heading data; Assume sample set Set the basis function Let polynomial for: Where k represents the degree of each factor in the polynomial, and M represents the highest degree of the polynomial; The extraction coefficient matrix is ​​θ=[a0 a1 … a M ] T , the basis function matrix is The matrix form is: Let Φ(t)=[φ(t0) φ(t1) … φ(t n )] T , let the regression function under ideal state be F(x), select the two-norm l2 as the weighted loss function, set it as J: Among them, w i represents the weight, F(x)=[f(x0) f(x1) … f(x n )] T , the weight matrix is ​​calculated by the Gaussian kernel: Derivative of J(θ): θ=(Φ T (WΦ) -1 F T WF (9) In formula (3), when M=1, the model is a linear regression model; the polynomial fitted based on the existing valid data is used to predict the data of the current shot, and the predicted value is: Use the predicted value φ(t n ) replaces the invalid position and heading data of the current shot.

3. The method for real-time processing of azimuth observation data with a backtracking repair function according to claim 2, characterized in that: In step S300, specifically including: Current shot data t n By n The data before the data is taken Y(t n-1 )The polynomial prediction value obtained by fitting is The actual observed data is y i , Calculate the Kalman filter innovation, the tth n The new value of the data is υ(t n )as follows: Among them, y(t n ) represents the tth n The actual observation data of the shooting data; Assume that receiving n When taking data, the sequence of the first N-1 new messages is as follows: I(t n-1 )={υ(t n-N ),υ(t n-N+1 ),...,υ(t n-1 )} (12) Its mean is μ I (t 1:n-1 ), with variance σ I (t 1:n-1 ), the chi-square value is calculated as: The probability density function of the chi-square distribution is: The cumulative distribution function of the chi-square distribution is: Among them, Γ(·) is the gamma function, γ(·) is the incomplete gamma function, and d is the degree of freedom; The significance level is α = 0.05, and the critical value is 3.

84. 2 >3.84, then the t n The data is an outlier.

4. The method for real-time processing of azimuth observation data with a backtracking repair function according to claim 3 is characterized in that: In step S400, the method for determining the outlier type includes: If there are continuous outliers not less than l SPK If the value is σ, it is considered to be a spot type outlier: For the azimuth observation data, find the first time the data deviates from 1 standard deviation σ SPK Position t SPK , take the position as the starting point of the azimuth speckle type wild value, use the fitting polynomial of the position to predict the data in the speckle type wild value interval shot by shot, and replace the observed data; after that, the observed data is detected to be no more than 1 times the standard deviation σ SPK The position of the spot type wild value is taken as the actual end position; for the heading observation data, look forward to find the first time that the data deviates from 1 standard deviation σ SPK Position t SPK , taking this position as the starting point of the heading maneuver, and then detecting that the observed data is no greater than 1 times the standard deviation σ SPK The position of is taken as the end position of heading maneuver; When the n The data is an outlier, if |υ(t n )-μ I (t 1:n-1 )|>4σ(t 1:n-1 ) and duration Δt u Less than l speck beat, then it is considered that the t n The data is a step-type outlier; look forward to find the first time the data deviates from 1 times the standard deviation σ speck1 Position t speck1 , take this position as the starting point of the azimuth step wild value, use the polynomial fitted before this position to predict the data in the step wild value interval shot by shot, and replace the observed data; after that, the observed data is detected to be no more than 1 times the standard deviation σ speck1 The position of is taken as the actual end position of the step wild value.

5. The method for real-time processing of azimuth observation data with a backtracking and repairing function according to claim 4 is characterized in that: In step S500, specifically including: Assume the window width is N, and the data at each position within the window width is recorded as t, which can be expressed as: t i ∈t={t i |0<t n -t i ≤N} (16) The Gaussian kernel is expressed as: Other formulas are shown in (3)-(10). n The result of smoothing the data is:

6. The method for real-time processing of azimuth observation data with a backtracking and repairing function according to claim 5, characterized in that: In step S600, specifically including: Let the current time be the tth n Shot data, assuming the last shot of the spot to be traced back is Let the set be Other outliers in the range belong to the set S OTH , to satisfy Backtracking of wild values; Assume that the center of the spot type outlier is The length of the spot type field value is 2r, and the start of the spot type field value is Shoot, end The wild value interval of the spot type is: The length of the wild-value interval is 2r, the total window width is 2N+2r, and the interval π is selected in the complete azimuth data set. REV : According to formula (17), the Gaussian kernel W is obtained: In the entire orientation data set, excluding S SPK The outlier interval is used to perform local polynomial regression on the data before and after the outlier interval. The local polynomial regression method is used to back-repair each value in the wild value interval of length 2r; suppose that the sth value in the interval i If the data is shot for backtracking repair, the corresponding Gaussian kernel interval is as follows: The Gaussian kernel is expressed as: Get a series of weight values: The Gaussian kernel weight matrix is ​​expressed as: The sth value in the outlier interval is obtained by formula (8)-(10): i The backtracking repair value of the shot data.

7. A real-time processing system for azimuth observation data with a backtracking repair function, characterized in that: The system has a program module corresponding to the steps of any one of claims 1 to 6, and executes the steps of the above-mentioned real-time processing method of azimuth observation data with backtracking repair function when running.

8. A computer-readable storage medium, characterized in that: The computer-readable storage medium stores a computer program, and the computer program is configured to implement the steps of the real-time processing method for azimuth observation data with a retrospective repair function described in any one of claims 1 to 6 when called by a processor.

Images