Maneuvering target tracking method and system based on orthogonal decomposition and parameter adaptation

Abstract

The invention discloses a maneuvering target tracking method and system based on orthogonal decomposition and parameter adaptation. The method comprises the following steps: acquiring an initial state of a target; a self-adaptive interactive multi-model filtering framework is provided, complex three-dimensional space motion of a target is orthogonally decomposed into a horizontal plane and two perpendicular planes which are perpendicular to each other, and complete description and tracking are carried out on any space maneuver; and designing an online parameter adaptive mechanism, intelligently adjusting the smooth strength of the quasi-uniform B spline according to the real-time maneuvering mode probability output by the IMM filter, robustly identifying the turning angular velocity in each plane through a multi-section geometric average strategy based on the optimal smooth trajectory, and dynamically feeding back the identification result to the corresponding collaborative turning model. According to the method, a complex high-dimensional nonlinear tracking problem is converted into a plurality of parallel and adaptive low-dimensional problems for solving, and finally, high-precision, fast-convergence and strong-robustness tracking of the high-maneuvering target is realized.

Description

Technical Field

[0001] This invention relates to the field of maneuvering target tracking technology, and in particular to a maneuvering target tracking method and system based on orthogonal decomposition and parameter adaptation. Background Technology

[0002] In fields such as UAV swarm cooperative positioning and air traffic control, accurate and real-time tracking of highly dynamic and maneuverable flying targets is a crucial core technology. Traditional target tracking methods, such as Kalman filters using single motion models like uniform velocity (CV) or uniform acceleration (CA), can achieve good tracking results when the target is in stable, non-maneuvering flight.

[0003] However, in real-world scenarios, targets often perform complex maneuvers, including but not limited to high-speed turns, S-shaped maneuvers, spiral climbs, or dive turns. In such cases, the target's true motion pattern exhibits a significant mismatch with any single fixed model. Continuing to use a single model filter will lead to several serious drawbacks:

[0004] (1) Tracking accuracy drops sharply: When the target starts to maneuver, the filter based on the non-maneuvering model cannot respond to its state changes in time, resulting in a rapid increase in tracking error.

[0005] (2) Risk of filter divergence: The continuous accumulation of model mismatch may lead to errors in the calculation of the covariance matrix of the filter, which will eventually cause the estimation result to deviate completely from the true value, i.e. filter divergence, resulting in tracking loss.

[0006] To address the limitations of single-model tracking of maneuvering targets, the Interactive Multiple Model (IMM) algorithm was proposed and has been widely applied. The IMM algorithm runs multiple different motion models in parallel, dynamically adjusting their probability weights based on the degree of matching between each model and actual observation data, and finally weighted and fused the state estimates of all models. This method can adapt to target maneuvering changes to a certain extent.

[0007] Nevertheless, traditional IMM algorithms still face the following challenges and shortcomings in engineering practice:

[0008] (1) The contradiction between the completeness of the model set and the computational burden. The performance of the traditional IMM algorithm is highly dependent on the pre-defined model set. In order to cover all possible maneuver modes of the target, a large and complete model library is required. However, this will greatly increase the computational burden, which is unacceptable for systems that require high real-time performance. Conversely, if the model set is too small or poorly designed, the performance of the traditional IMM algorithm will still be limited and cannot reach the optimal level when the actual maneuver of the target (e.g., turning angular velocity) is not included in any of the pre-defined models.

[0009] (2) Fixed model parameters and uncertainty of target maneuvering. The maneuvering model in traditional IMM algorithms (such as the cooperative turning CT model) usually has fixed parameters, such as a few preset typical turning angular velocities. However, the actual turning angular velocity of the target is continuously changing and unpredictable. Fixed model parameters mean that the algorithm can only "approximate" but not "match" the target's real maneuvering, which will still introduce a significant tracking error in strong maneuvering scenarios.

[0010] (3) Model switching delay affects instantaneous accuracy. In the traditional IMM algorithm, the update of model probability depends on the cumulative effect of information (residual). When the target motion pattern changes abruptly, the transfer of model probability often has a certain delay. This lag causes the filter to still assign higher weights to the old, mismatched model in the initial stage of target maneuver switching, thus affecting the instantaneous response speed and accuracy of tracking.

[0011] Therefore, how to design a method that can adjust its internal model parameters in real time and adaptively under limited computing resources to accurately match the unknown and time-varying maneuvering modes of the target, thereby achieving fast, robust and accurate tracking of highly dynamic targets, is a technical problem that urgently needs to be solved in the field of target tracking. Summary of the Invention

[0012] The purpose of this invention is to provide a maneuvering target tracking method and system based on orthogonal decomposition and parameter adaptation, so as to achieve accurate estimation of the position and motion state of highly maneuvering flying targets.

[0013] The technical solution for achieving the objective of this invention is as follows: Firstly, this invention provides a maneuvering target tracking method based on orthogonal decomposition and parameter adaptation, comprising the following steps:

[0014] S1. Initialization: The target azimuth and elevation angle measurements are obtained by at least two passive observation platforms. An initial three-dimensional position sequence is obtained based on multi-frame triangulation and smoothing. The initial velocity is obtained by least squares linear fitting, forming the initial state variables.

[0015] S2. Adaptive Interactive Multi-Model Filtering by Orthogonal Decomposition: Construct a model set consisting of uniform velocity, uniform acceleration, and cooperative turning models corresponding to a horizontal plane and two mutually perpendicular vertical planes, respectively, to accurately describe arbitrary spatial maneuvers. Then, execute the filtering process of model interactive mixing, parallel extended Kalman filtering, model probability update, and state fusion.

[0016] S3. Online parameter adaptation: Based on the fused trajectory in step S2, the trajectory is smoothed using an adaptive quasi-uniform B-spline within a fixed-length sliding window, and the real-time turning angular velocity of each plane is calculated using multi-segment three-point geometric curvature in three orthogonal planes; the smoothing intensity of the adaptive quasi-uniform B-spline is adaptively adjusted by the model probability output in step S2, and the identified angular velocity is dynamically written back to the state transition matrix of the corresponding cooperative turning model for use in the next filtering cycle.

[0017] Furthermore, the initialization in step S1 specifically includes:

[0018] S11: For the two passive observation platforms in the... The positions of the first and second observation platforms are as follows: and The measured target azimuth and elevation angles are ( , )and( , A unit direction vector pointing towards the target is constructed from measurements:

[0019] ;

[0020] in, Construct a system of linear equations:

[0021] ;

[0022] in, It is in the The coefficient matrix of the linear equation system constructed in the frame consists of two unit direction vectors. It is the vector of constant terms in the linear equation system, representing the position difference vector between the two observation platforms. The distance parameter vector to be solved is as follows:

[0023] ;

[0024] in and This represents the scalar distance from the first and second observation platforms to the target point; the distance parameters are solved using the least squares method. :

[0025] ;

[0026] Get the first Initial position estimate of the frame Triangulation is performed on all the set initial processing frames to obtain the initial position sequence;

[0027] S12: Perform a moving average on the initial position sequence obtained in step S11 to obtain the processed position sequence, and take the nearest value. The arithmetic mean of the frames is used as the initial position of the filter. ;

[0028] S13: Use the position sequence after moving average to determine the target's three coordinates Perform a first-order least-squares linear fit over time. ,in It describes the target position over time. A changing function, and It is the slope and intercept of the first-order minimum fitted line, and the fitted slope is taken. As the initial velocity and with physical limiting;

[0029] S14: Set the initial triaxial accelerations to 0, i.e., the initial acceleration vector is... , Given the initial velocity obtained in step 3, we obtain the initial state. And establish the initial state covariance in block diagonal form. :

[0030] ;

[0031] in , , These are the prior variance matrices for position, velocity, and acceleration, respectively.

[0032] Furthermore, the uniform velocity CV, uniform acceleration CA, and cooperative turning CT model set in step S2 specifically include:

[0033] The state transition matrices for uniform velocity CV and uniform acceleration CA are as follows:

[0034] ;

[0035] ;

[0036] in It is the sampling period. It is a 3x3 identity matrix;

[0037] For any CT plane Its turning angular velocity Instead of preset fixed parameters, the parameters are identified and updated in real time during each filtering cycle based on the online parameter adaptive process in step S3; assuming the turning angular velocity of each plane is... ,make Its velocity rotation matrix and position-velocity coupling submatrix for:

[0038] The position-velocity coupling submatrix is:

[0039] , ;

[0040] , ;

[0041] , ;

[0042] Velocity rotation matrix Described in a plane Inside, after time Then, how does the velocity vector rotate; the position-velocity coupling submatrix. Described in a plane Internally, the effect of velocity changes on position; Acting on the corresponding planar velocity components, Acting on the corresponding planar position-velocity coupling, embedded Matrix obtained :

[0043] .

[0044] Furthermore, the measurement model in the parallel extended Kalman filter in step S2 specifically includes:

[0045] In the implementation of the two observation platforms, for the first One observation platform, Let the target relative quantity be , , Horizontal distance Total distance Take measurement vector azimuth With pitch angle definition:

[0046] ;

[0047] Its Jacobian matrix The non-zero elements are only located in the first three columns corresponding to the positional components, satisfying:

[0048] ;

[0049] ;

[0050] For angular residuals Perform surround normalization:

[0051] ;

[0052] in, It is the angle residual normalization function, where r is the angle value to be normalized. This indicates the function to round down; then EKF update is performed.

[0053] Furthermore, the adaptive interactive multi-model filtering in step S2 specifically includes:

[0054] Based on the model probability of the previous time step Total number of models and the Markov transition probability matrix between models elements in Calculate the probability of the prediction model and mixture probability :

[0055] ;

[0056] ;

[0057] Based on the posterior state estimates of each model at the previous time step covariance For each model Calculate its mixed initial state and covariance. ;

[0058] For each model Given a mixed initial state, a state transition model is used for prediction, and EKF updates are performed in parallel to obtain the posterior state of each model. Covariance residual and new information covariance ;

[0059] According to each model residual New information covariance and measurement vector dimension Calculate the likelihood function value :

[0060] ;

[0061] Given the total number of models, the likelihood values ​​are weighted and normalized based on the predicted model probabilities to obtain the posterior model probability at the current time step. :

[0062] ;

[0063] The model probability is used both for the final state fusion in the next step and for the control basis for adaptively adjusting the smoothing intensity of the adaptive QUBS in step S3.

[0064] Posterior state estimation of all filter outputs Based on its latest model probability We perform weighted fusion to obtain the final target state estimate. Covariance :

[0065] ;

[0066] .

[0067] Furthermore, the online parameter adaptation in step S3 includes:

[0068] S31: Set the sampling period Segment duration ,make:

[0069] ;

[0070] Window length :

[0071] ;

[0072] The length is maintained in a first-in, first-out manner. position buffer ;

[0073] S32: Construct the objective function :

[0074] ;

[0075] in For the basis function matrix, It is a second-order difference matrix. These are the control point vectors of the B-spline, and the unknowns to be solved. They are obtained by solving... Determine a smooth curve. It is the adaptive smoothing intensity;

[0076] S33: Probability from CV model Calculate smoothing strength :

[0077] ;

[0078] The upper and lower limits are defined as follows: and exponential smoothing Suppressing oscillations, smoothing coefficient To obtain the adaptive smoothing intensity ;

[0079] S34: By integer offset Construct multiple sets of three points in a loop. For each set of three points, select the three-point index as... , Using the current time marker, the instantaneous velocity at the endpoint is estimated using second-order backward difference. The signed discrete curvature on each orthogonal plane is calculated using the three-point geometric relationship. Thus, the initial angular velocity is obtained. ;

[0080] S35: Calculate the arithmetic mean of all the angular velocity values ​​calculated in each group and subject to physical amplitude constraints to obtain the estimated turning angular velocity in the three planes at the current moment. It is then dynamically written into the state transition matrix of the corresponding CT model for use in the next filtering cycle.

[0081] Secondly, the present invention provides a maneuvering target tracking system based on orthogonal decomposition and parameter adaptation, for implementing the method described in the first aspect, the system comprising:

[0082] The first module is used to obtain the target azimuth and elevation angle measurements through at least two passive observation platforms, obtain an initial three-dimensional position sequence based on multi-frame triangulation and smoothing, and obtain the initial velocity by least squares linear fitting to form the initial state quantity.

[0083] The second module is used to construct a model set consisting of uniform speed, uniform acceleration, and cooperative turning models corresponding to a horizontal plane and two mutually perpendicular vertical planes, respectively, to accurately describe arbitrary spatial maneuvers. Then, it executes a filtering process of model interaction and mixing, parallel extended Kalman filtering, model probability update, and state fusion.

[0084] The third module is used to smooth the trajectory within a fixed-length sliding window using adaptive quasi-uniform B-splines based on the fusion trajectory of the second module, and to calculate the real-time turning angular velocity of each plane in three orthogonal planes using multi-segment three-point geometric curvature. The smoothing intensity of the adaptive quasi-uniform B-splines is adaptively adjusted by the model probability output in step S2, and the identified angular velocities are dynamically written back to the state transition matrix of the corresponding cooperative turning model for use in the next filtering cycle.

[0085] Thirdly, the present invention provides a computer device including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the program to implement the steps of the method described in the first aspect.

[0086] Fourthly, the present invention provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the steps of the method described in the first aspect.

[0087] Fifthly, the present invention provides a computer program product, including a computer program that, when executed by a processor, implements the steps of the method described in the first aspect.

[0088] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0089] (1) Significantly improves tracking accuracy: By estimating and updating the angular velocity parameters of the CT model online and in real time, the model of this invention can "match" rather than "approximate" the actual maneuver of the target. This accurate model matching greatly reduces the model mismatch error, especially when the target is performing strong maneuvers or long turns, the tracking accuracy is much higher than that of the traditional IMM algorithm using a fixed parameter model.

[0090] (2) Enhanced adaptability to complex three-dimensional maneuvers: The present invention constructs a set of turning models covering a horizontal plane and two mutually perpendicular vertical planes, and through the idea of ​​orthogonal decomposition, it can effectively handle various complex combined maneuvers of the target in three-dimensional space, such as spiral climb or dive turn, and has stronger robustness and wider applicability.

[0091] (3) Improved rapid response capability to maneuver changes: The closed-loop feedback mechanism of this invention enables the model parameters to be adjusted closely with the changes in the target's actual motion state. At the same time, the innovative adaptive QUBS smoothing mechanism can perform strong smoothing to suppress noise when the target is not maneuvering, and weak smoothing to preserve maneuver details when the target is maneuvering, ensuring the accuracy and real-time performance of parameter identification. This compensates for the lag in model probability updates in the traditional IMM algorithm, enabling the system to converge faster when the target's maneuver mode changes abruptly. Attached Figure Description

[0092] To more clearly illustrate the specific embodiments of the present invention or the technical solutions in the prior art, the drawings used in the description of the specific embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained from these drawings without creative effort.

[0093] Figure 1This is a schematic diagram of the aircraft detection method based on the adaptive interactive multi-model algorithm of the present invention;

[0094] Figure 2 This is an example diagram illustrating a scenario of the present invention;

[0095] Figure 3 This is a flowchart of the initialization module of the present invention;

[0096] Figure 4 This is a flowchart of the CT model parameter adaptive module of the present invention;

[0097] Figure 5 This is a comparison chart of the angular velocity identification results of the present invention. Detailed Implementation

[0098] This invention proposes a passive cooperative tracking method for aerial 3D maneuvering targets based on orthogonal decomposition and parameter adaptation, aiming to solve the problems of low tracking accuracy and slow response caused by model fixation, incomplete description, and slow startup in existing technologies. First, this method employs an initialization algorithm to quickly obtain the target's high-precision initial state, enabling rapid tracking convergence. Second, the core of this invention lies in an adaptive interactive multi-model (IMM) filtering framework. This framework innovatively orthogonally decomposes the complex 3D spatial motion of the target into a horizontal plane and two mutually perpendicular vertical planes, thereby providing a complete description and tracking of arbitrary spatial maneuvers in a "decomposition-parallel filtering-fusion" mode. This invention further designs an online parameter adaptation mechanism. The core of this mechanism lies in the deep coupling of trajectory smoothing and parameter identification: based on the real-time maneuver mode probability output by the IMM filter, the mechanism intelligently adjusts the smoothing intensity of the quasi-uniform B-spline (QUBS), and based on this optimal smoothed trajectory, robustly identifies the turning angular velocity in each plane through a multi-segment geometric averaging strategy. Finally, the identification results are dynamically fed back to the corresponding cooperative turning (CT) model. This invention, through the organic combination of initialization, decomposition tracking and parameter adaptation, creatively transforms a complex high-dimensional nonlinear tracking problem into multiple parallel, adaptive low-dimensional problems for solution, ultimately achieving high-precision, fast-convergence and robust tracking of highly maneuverable targets.

[0099] The technical solution of the present invention will now be described in detail with reference to the accompanying drawings.

[0100] This invention addresses the complexity and uncertainty of the three-dimensional spatial motion of highly maneuverable aerial targets, as well as the low tracking accuracy and slow response of traditional tracking methods due to fixed model parameters and incomplete descriptions of maneuver patterns. It proposes a high-precision tracking method for maneuvering targets based on orthogonal decomposition and parameter adaptation. This method mainly consists of three parts: an initialization module, an innovative orthogonal decomposition adaptive IMM filtering module, and a core online parameter adaptation module. The relationships between these parts are shown in the attached figure. Figure 1 As shown, a closed-loop feedback system of "filtering-identification-update" is formed, which can be divided into three main steps:

[0101] (1) Initialization: The target azimuth and elevation angle measurements are obtained by at least two passive observation platforms. An initial three-dimensional position sequence is obtained based on multi-frame triangulation and smoothing. The initial velocity is obtained by least squares linear fitting, forming the initial state quantity.

[0102] (2) Orthogonal decomposition adaptive interactive multi-model (IMM) filtering: Construct a model set consisting of uniform velocity (CV), uniform acceleration (CA) and cooperative turning (CT) models corresponding to a horizontal plane and two mutually perpendicular vertical planes respectively, to accurately describe arbitrary spatial maneuvers. Then, execute the filtering process of model interactive mixing, parallel extended Kalman filter (EKF), model probability update and state fusion.

[0103] (3) Online parameter adaptation: Based on the fused trajectory in step S2, adaptive quasi-uniform B-spline (QUBS) is used for trajectory smoothing within a fixed-length sliding window, and the real-time turning angular velocity of each plane is calculated using multi-segment three-point geometric curvature in three orthogonal planes; the smoothing intensity of the adaptive QUBS is... The model probability output from step S2 is adaptively adjusted, and the identified angular velocity is dynamically written back to the state transition matrix of the corresponding CT model for use in the next filtering cycle.

[0104] The specific implementation methods for each of the above steps are explained in detail below:

[0105] Step 1: Initialization process

[0106] This embodiment uses two observation platforms as an example for illustration, but the invention is also applicable to scenarios with more than two observation platforms. A "Local Fixed Northeast (ENU)" coordinate system is defined as a globally unified reference coordinate system. The origin is set at the initial position of the first observation platform, with its X-axis pointing due east, Y-axis pointing due north, and Z-axis perpendicular to the local horizontal plane pointing upwards (zenith). The positions and target states of all observation platforms are represented in this coordinate system. A specific scenario is as follows: Figure 2 As shown.

[0107] The initialization process used in this invention is as follows: Figure 3 As shown, using the pure angle information from the initial sampling periods, a total of The first frame estimates the initial three-dimensional position sequence of the target using triangulation. In the second frame... The frames and the positions of the first and second observation platforms are respectively and The measured target azimuth and elevation angles are respectively ( , )and( , Construct a unit direction vector pointing towards the target:

[0108]

[0109] in Ideally, the two lines of sight should intersect at a single point, namely the target location. In reality, due to measurement errors, the two observation lines usually do not intersect precisely. Therefore, the problem is transformed into solving an overdetermined system of equations. The least squares problem, where, It is in the The coefficient matrix of the linear equation system constructed in the frame consists of two unit direction vectors. It is the vector of constant terms in the linear equation system, representing the position difference vector between the two observation platforms. Let be the distance parameter vector to be solved.

[0110]

[0111] in and This represents the scalar distance from the first and second observation platforms to the target point; the distance parameters are solved using the least squares method. :

[0112]

[0113] The first distance parameter obtained The initial position estimate of the j-th frame is obtained. Repeat this process for all initialization frames to obtain the initial position estimation sequence.

[0114] The obtained initial position estimation sequence is smoothed using a moving average filter to suppress noise and outliers. Specifically, a length of [length missing] is used. The window, for the first The first frame 3D coordinates Perform smoothing calculations:

[0115]

[0116] Will Combining according to the x, y, z dimensions in order This yields a smoothed complete sequence of initial positions. Then, to further improve the stability of the initial positions, the last position of this smoothed trajectory is taken. The arithmetic mean of the frames is calculated. As the final initial position for filter activation:

[0117]

[0118] Based on the position sequence obtained above The initial velocity of the target is estimated by fitting a first-order polynomial (linear line) to each coordinate dimension. For each dimension, the linear model to be solved is:

[0119]

[0120] in yes Values ​​in a vector It is the velocity (slope) to be determined. It is the first The timestamp corresponding to the frame data point. It is the intercept of the linear model to be solved. For error terms; by constructing the design matrix and observation vector This linear regression problem can be solved using the least squares method with a closed-form solution:

[0121]

[0122] Take the slope As the initial velocity for this dimension, an upper and lower bound constraint is imposed on the estimated velocity value:

[0123]

[0124] Among them, the upper and lower limits of speed and These are preset parameters. This represents the final determined initial velocity value in the d-dimensional plane. Different constraints can be imposed on the velocities in the three directions, and these constraints can be combined to obtain the initial velocity vector. .

[0125] After completing the above positioning and fitting, this invention uses a 9-dimensional state vector. A complete description of the target's three-dimensional motion state:

[0126]

[0127] in: For location; For speed; It is acceleration.

[0128] Assuming the initial acceleration is zero, construct the initial state vector of the filter. :

[0129]

[0130] In this invention, the initial covariance matrix It is constructed as a 9x9 block diagonal matrix, reflecting that the estimation errors of initial position, velocity, and acceleration are independent of each other, with the following structure:

[0131]

[0132] in:

[0133] It is a 3x3 diagonal matrix representing the covariance of the initial position. Its diagonal elements are the variances of the uncertainties of the initial position in the three directions; It is a 3x3 diagonal matrix representing the covariance of the initial velocity. Its diagonal elements are the variances of the uncertainty of the initial velocity in the three directions; It is a 3x3 diagonal matrix representing the covariance of the initial acceleration.

[0134] Process noise covariance matrix This invention describes the differences between the motion model and real motion. For each model... Define a process noise covariance matrix respectively. For the first A model is defined with its process noise standard deviation set as . Through a 9x3 noise transfer matrix Mapping continuous acceleration noise to a discrete state vector:

[0135]

[0136] in, It is the sampling period. It is a 3x3 identity matrix.

[0137] Finally, the first The process noise covariance matrix of each model is:

[0138]

[0139] Measurement noise covariance matrix This describes the degree of uncertainty in the sensor measurements. It is assumed that the four angle measurement errors are independent and follow a zero-mean Gaussian distribution with a standard deviation of... Then the measurement noise covariance matrix can be set as follows:

[0140]

[0141] in, It is a 4x4 identity matrix.

[0142] The above steps provide a complete and reliable set of initialization parameters for subsequent adaptive IMM filtering, including initial state, state covariance, process noise, and measurement noise.

[0143] Step 2: Adaptive Interactive Multiple Model (IMM) Filtering Process Based on Orthogonal Decomposition

[0144] This invention employs five different motion models to track the target state in parallel, with each model corresponding to a 9-dimensional state transition matrix to comprehensively describe the target's various maneuvers. These five models specifically include:

[0145] CV Model: This model is suitable for targets performing non-maneuvering or weakly maneuvering flight phases. Its core physical assumption is that the target's acceleration is zero and its velocity remains constant. State transition matrix. It can be represented as:

[0146]

[0147] CA Model: This model is applicable to phases of sustained linear acceleration or deceleration of a target. Its core physical assumption is that the target's acceleration remains constant. State transition matrix. It can be represented as:

[0148]

[0149] CT model: for any CT plane Its turning angular velocity These are not preset, fixed parameters, but are identified and updated in real time during each filtering cycle; Its velocity rotation matrix With position-velocity coupling matrix for:

[0150] ,

[0151] ,

[0152] ,

[0153] Will Acting on the corresponding planar velocity components, Acting on the corresponding planar position-velocity coupling, embedded Matrix obtained :

[0154]

[0155]

[0156]

[0157] As can be seen from the above, the three CT models describe turning maneuvers in three mutually orthogonal planes in space (horizontal plane XOY, vertical plane YOZ, and vertical plane XOZ). When the target maneuvers with a complex spatial trajectory (e.g., having turning components in both the horizontal and vertical directions simultaneously), different CT models can each capture its turning components in the corresponding plane, thus achieving an accurate description of maneuvers in any direction in space through the combination of model sets.

[0158] Based on the above model set, the interactive multi-model filter will strictly execute the following four core steps in each sampling period:

[0159] (1) Interactive hybrid: based on the posterior state estimates of each model at the previous time step covariance probability and the Markov transition probability matrix between models For each model Calculate its mixed initial state. First, calculate the model. Predicted probability :

[0160]

[0161] in This represents the total number of models. Next, we calculate the number of models... To model Mixed probability :

[0162]

[0163] Then, using these mixture probabilities as weights, the states of each model at the previous time step are mixed proportionally to obtain the model. Mixed initial state estimation Covariance :

[0164]

[0165]

[0166] (2) Parallel filtering: For each model, starting with the mixed initial state and covariance, perform an extended Kalman filter (EKF) process once:

[0167] Based on the posterior state estimate and covariance at time k-1, the prior state estimate is obtained. :

[0168]

[0169] The covariance prediction equation yields the prior covariance. :

[0170]

[0171] in For the model The state transition matrix (defined by the model set). This is the process noise covariance matrix of the corresponding model.

[0172] Measurement prediction: based on prior state estimation Through nonlinear measurement functions Calculate theoretical measurement prediction values .

[0173] Jacobian matrix calculation: Calculating measurement functions Jacobian matrix of the state vector The matrix is ​​a 4x9 matrix, with its non-zero elements concentrated in the first three columns. Let the target predicted position be the same as the observation position. The relative position vector is , , Let the horizontal distance The square of the total distance . All non-zero elements in the matrix are as follows:

[0174] For the azimuth angle of observation machine 1 Partial derivative:

[0175]

[0176] For the azimuth angle of the observation machine 2 Partial derivative:

[0177]

[0178] Pitch angle of observation instrument 1 Partial derivative:

[0179]

[0180] Pitch angle of observation instrument 2 Partial derivative:

[0181]

[0182] Calculate actual measurements With predictive measurement residuals between :

[0183]

[0184] Since the measurement includes angles, the angular component in the residual needs to be normalized to ensure its value is within the acceptable range. Within the interval:

[0185]

[0186] Where r is the angle value to be normalized. This represents the floor function.

[0187] Calculate the covariance matrix of the residuals :

[0188]

[0189] Kalman gain calculation :

[0190]

[0191] State and covariance update: obtaining the model posterior estimation Posterior covariance The update uses the Joseph form, which has better numerical stability:

[0192]

[0193]

[0194] in It is an identity matrix.

[0195] (3) Model probability update: based on the residuals and innovation covariance of each filter and the dimension of the measurement vector. Calculate its likelihood function value This value follows a Gaussian distribution:

[0196]

[0197] Combining prediction model probabilities We obtain the posterior model probability at the current time by weighting and normalizing the likelihood values. :

[0198]

[0199] (4) State fusion: The posterior state estimates of all filter outputs are weighted and fused according to their latest model probabilities to obtain the final target state estimate. Covariance :

[0200]

[0201]

[0202] Step 3: Online parameter adaptation process

[0203] This step is crucial for achieving adaptive high-precision tracking. Its core objective is to identify the target's current turning angular velocity in various orthogonal planes in real time along the time axis based on the historical trajectory fused from the filter, and then use the identification results to dynamically update model parameters, forming a closed loop. The online parameter adaptation process is as follows: Figure 4 As shown.

[0204] Trajectory data caching: The system will store the latest fused output of the target's 3D position. Store a fixed length First-in-first-out sliding window buffer To balance response delay and statistical stability in maneuver identification, the window length... Selection and segment duration Associated. Based on sampling interval Calculate the frame interval of the segment. :

[0205]

[0206] Total length of sliding window Set as:

[0207]

[0208] This design ensures that there are always enough data points in the buffer to support the subsequent construction of three intervals. The calculation is performed on the three-point group of the frame.

[0209] Adaptive smoothing preprocessing: When identifying maneuver parameters from historical trajectories, directly calculating curvature on noisy discrete trajectory data greatly amplifies the noise, leading to unstable identification results. To address this issue, this invention introduces a quasi-uniform B-spline (QUBS) to buffer the trajectory before calculating curvature. The data points within the range are preprocessed for smoothing. This method achieves smoothing by solving a regularized least squares problem with a smoothing penalty term. Its core is defined by the following objective function. At the same time, it takes into account both data fidelity and curve smoothness:

[0210]

[0211] This objective function consists of two key parts, through a smoothing factor. To weigh:

[0212] (1) Fitting error term : Requires the B-spline basis function matrix and the vector of control point coefficients to be determined Defined smoothed curve , compared with the original noisy data points The sum of squared errors between these parameters is minimized. This ensures that the final smooth curve does not deviate excessively from the true trend of the original trajectory.

[0213] (2) Smoothing penalty term :in, It is a second-order difference matrix, and its effect on the control point coefficients is approximately equivalent to calculating the second derivative of a B-spline curve. Smoothing factor. It is an adjustable hyperparameter used to balance data fidelity and curve smoothness.

[0214] In the objective function, construct the B-spline basis function matrix. Choose cubic splines, i.e., spline order. Based on the number of window frames and sampling interval Construct the parameter sequence:

[0215]

[0216] Constructing the junction vector using a quasi-uniform method Its form is that there are at both ends Repeated nodes are used to ensure that the curve passes through both the first and last data points:

[0217]

[0218] B-spline basis functions Defined by the Cox-de Boor recurrence relation. For :

[0219]

[0220] in, For the i-th node, This is the (i+1)th node. For higher orders... :

[0221]

[0222] basis matrix The number of control points is basis function matrix Each element That is, at the parameter point The above is the first indivual The results of evaluating the first-order basis functions: .

[0223] Construct a sparse matrix Its structure is as follows, used to approximate the second derivative of a curve:

[0224]

[0225] To minimize the objective function Regarding its control point coefficients Find the gradient and set it to zero to obtain the normal equation for solving the optimal control point:

[0226]

[0227] The solution is:

[0228]

[0229] in, These are the control point coefficients obtained by solving for the three-dimensional coordinates. Finally, multiplying the obtained optimal control point coefficients by the basis function matrix yields the smoothed trajectory.

[0230]

[0231] This smoothed trajectory This will serve as input for subsequent curvature and angular velocity calculations, providing a stable and reliable data foundation for the accurate identification of maneuver parameters.

[0232] To balance noise immunity and maneuver detail preservation, this embodiment introduces a model probability-driven adaptive adjustment mechanism for smoothing intensity in the QUBS smoothing process. Specifically, in each filtering cycle... At the end, the IMM filter gives the posterior model probability vector. The elements inside represent the probabilities of the five models at this moment:

[0233]

[0234] Define the range of values ​​for the smoothing factor And the probability of the current CV model. Mapped to the current smoothing factor:

[0235]

[0236] when (This indicates that the target is in a straight, constant-speed segment or a non-maneuvering segment.) The smoothing process will place greater emphasis on noise reduction and smoothing to filter out measurement noise; when (This indicates that the target is currently in a maneuvering turning phase.) The smoothing process weakens the constraints to preserve the curvature details of the trajectory as much as possible. This is to suppress model probability jitter that causes... Rapid fluctuations can be further smoothed using exponential smoothing methods:

[0237]

[0238] in: This is the smoothing coefficient. Ultimately, it will... Substituting these values ​​into the QUBS solution process above, we obtain the smoothed trajectory for the current period:

[0239]

[0240]

[0241]

[0242] After obtaining the smooth trajectory, this invention employs a multi-segment three-point group averaging strategy to identify the target's turning angular velocity on three orthogonal planes (xoy, yoz, xoz). This strategy aims to improve the robustness and noise resistance of the identification by fusing calculation results under multiple geometric configurations. Specifically, the system uses an integer offset... (From 0 to ) Loop construction Three-point groups, among which This is the preset frame interval. For each... From smooth trajectory The three points selected are indexes. ,in This is the index of the latest point in the buffer at the current moment. Let these three points be... , , superscript of variables Both indicate that the physical quantity is generated by the first... Intermediate results obtained from group calculations, including:

[0243]

[0244] The terminal point is estimated using the three-point second-order backward difference formula. The instantaneous velocity. With Axial velocity components For example, the calculation formula is as follows:

[0245]

[0246] in , , There are three points respectively , , of The coordinate axes can be calculated similarly. Axial velocity components and Axial velocity components Then, the velocity vector is projected onto three orthogonal planes, and the magnitude of the velocity in each plane is calculated. In the xoy plane:

[0247]

[0248] The same applies to the other planes.

[0249] For each set of three points, project its three-dimensional coordinates onto three orthogonal planes xoy, yoz, and xoz, respectively, and calculate its signed discrete curvature on these planes using the geometric relationships between the three points. Side length These are the pairwise Euclidean distances between three points on the projection plane:

[0250]

[0251] Directed area :

[0252]

[0253] in, These are the projected two-dimensional coordinates.

[0254] radius of curvature :

[0255]

[0256] Signed curvature :

[0257]

[0258] Based on the basic principles of turning kinematics Calculate the first The initial angular velocities of the three-point group in each plane :

[0259]

[0260] The initial angular velocity is physically constrained to ensure it does not exceed the preset maximum angular velocity. The angular velocity value after physical constraint processing is obtained. :

[0261]

[0262] The arithmetic mean of the valid and constrained angular velocity values ​​calculated from all groups is used to obtain the most reliable estimate of the turning angular velocity in the three planes at the current moment. :

[0263]

[0264] in This is the estimated turning angular velocity on a certain plane. It is a group of all calculated effective angular velocity values. The set, It is the number of elements in the set.

[0265] Successfully identified the turning angular velocity on three orthogonal planes. Subsequently, the system will use these real-time maneuvering parameters to update the state transition matrix of the CT models in the model set, replacing the turning rate parameters in Model 3, Model 4, and Model 5 respectively, for use in the next stage of adaptive IMM filtering. In this way, the model parameters can be adaptively adjusted according to the changes in the target's actual maneuvering state, completing a closed-loop feedback process that combines online parameter identification and model updating.

[0266] Example

[0267] This embodiment simulates a scenario where two observation aircraft perform cooperative passive tracking of a single highly maneuverable target. The total simulation duration is 80 seconds, and the data sampling period is... The initial filter values ​​are given by using the initial 15 frames.

[0268] Observation aircraft 1 starts at the origin of the coordinate system [0,0,0] (m) and flies along the Y-axis (north direction) at a constant speed of [0,300,0] (m / s). Observation aircraft 2 starts at [10000,0,1500] (m) and also flies along the Y-axis at a constant speed of [0,300,0] (m / s), maintaining a parallel flight path to observation aircraft 1.

[0269] The tracked target's trajectory consists of a horizontal S-shaped maneuver. In the initial phase (0~16.67s), the target moves at a constant speed of 300m / s from its initial position [5000, 50000, 800] (m) along the negative Y-axis. During the maneuver phase (16.67s~72.67s), the target performs two consecutive horizontal S-shaped maneuvers, each with a turning angular velocity of 5.3° / s. In the final phase (72.67s~80s), the target resumes its constant speed linear motion.

[0270] The measurement data in the simulation are not ideal values, but rather simulate the errors of real sensors. First, based on the actual positions of the observation platform and the target, the true geometric angles (azimuth and elevation angles) at each moment are calculated. Then, a zero-mean Gaussian white noise with a standard deviation of 1 mrad is superimposed on each angle value to generate the noisy measurement sequence used as the final filter input.

[0271] When preprocessing the initial triangulation trajectory, the length of the moving average smoothing window is... Set to 41; used to calculate the final initial position confirmation frame count. Set to 6; apply physical constraints to the initial velocities fitted by the least squares method, with the x and y axis velocity components limited to ±350 m / s and the z axis velocity component limited to ±70 m / s; initial state covariance matrix. Composed of three diagonal matrices The location covariance velocity covariance The initial acceleration covariance for each model standard deviation of diagonal elements respectively .

[0272] Assuming the standard deviation of angle measurement error If the value is 1 mrad, then the measurement noise covariance matrix is ​​set as follows:

[0273]

[0274] The standard deviation of process noise is set for each model. They are respectively These correspond to CV and CA, respectively. , , Model.

[0275] Set the initial model probability vector Markov transition probability matrix The probability that the model remains unchanged Set to 0.98, the probability of transitioning to any other model. Set it to 0.005.

[0276] Used to construct the segmented duration of a three-point group Set to 0.8 seconds; QUBS smoothing factor The value range is from 0.001 to 0.1; the physical constraint limit for the online identified turning angular velocity is set to ±12° / s.

[0277] A traditional IMM algorithm is set as a benchmark for comparison. This benchmark algorithm uses the exact same filter framework, state vector definition, and initialization process as the present invention. However, its key difference lies in that the three CT models in its model set use fixed, preset turning angular velocities, and it lacks the online parameter identification and model update capabilities of the present invention; it uses QUBS with a fixed smoothing factor. The value is fixed at 0.005; the root-mean square error (RMSE) of each parameter is used as the evaluation index. The calculation method is as follows:

[0278]

[0279] In the formula Indicates the total number of frames. and They represent in The actual value and the estimated value at time t.

[0280] To quantitatively evaluate the performance of the method of this invention, it was compared with the traditional IMM algorithm in the same simulation scenario. The tracking error analysis results are shown in Table 1. The data shows that the method of this invention significantly outperforms the traditional IMM algorithm in both position and velocity tracking accuracy. In position tracking, the root mean square error of the total distance of the method of this invention is improved by approximately 11% compared to the traditional IMM algorithm. Furthermore, the errors of the method of this invention are also lower than those of the traditional algorithm across all three axes: X, Y, and Z. In velocity tracking, the advantages of this invention are even more pronounced. Especially in the Y direction, the accuracy improvement is as high as 26.9%. In the X and Z directions, the velocity accuracy is also improved by approximately 6.7% and 7.2%, respectively.

[0281] Table 1

[0282]

[0283] The horizontal turning angular velocity estimated online in this embodiment is compared with the true value of the target maneuver set in the simulation scenario, and the result is as follows: Figure 5 As shown. The results demonstrate that the angular velocity curve estimated by this invention can track the changes in the actual angular velocity very accurately, proving the effectiveness of the online maneuver parameter identification module.

[0284] In summary, through simulation comparison with the traditional IMM algorithm, the method proposed in this invention demonstrates significant advantages in tracking accuracy, adaptability to complex maneuvers, and response speed to maneuver changes.

Claims

1. A maneuvering target tracking method based on orthogonal decomposition and parameter adaptation, characterized in that, Includes the following steps: S1. Initialization: The target azimuth and elevation angle measurements are obtained by at least two passive observation platforms. An initial three-dimensional position sequence is obtained based on multi-frame triangulation and smoothing. The initial velocity is obtained by least squares linear fitting, forming the initial state variables. S2. Adaptive Interactive Multi-Model Filtering by Orthogonal Decomposition: Construct a model set consisting of uniform velocity, uniform acceleration, and cooperative turning models corresponding to a horizontal plane and two mutually perpendicular vertical planes, respectively, to describe arbitrary spatial maneuvers, and execute a filtering process of model interactive mixing, parallel extended Kalman filtering, model probability update, and state fusion. S3. Online parameter adaptation: Based on the fused trajectory in step S2, the trajectory is smoothed using an adaptive quasi-uniform B-spline within a fixed-length sliding window, and the real-time turning angular velocity of each plane is calculated using multi-segment three-point geometric curvature in three orthogonal planes; the smoothing intensity of the adaptive quasi-uniform B-spline is adaptively adjusted by the model probability output in step S2, and the identified angular velocity is dynamically written back to the state transition matrix of the corresponding cooperative turning model for use in the next filtering cycle.

2. The method according to claim 1, characterized in that, The initialization in step S1 specifically includes: S11: For the two passive observation platforms in the... The positions of the first and second observation platforms are as follows: and The measured target azimuth and elevation angles are ( , )and( , A unit direction vector pointing towards the target is constructed from measurements: ； in, Construct a system of linear equations: ； in, It is in the The coefficient matrix of the linear equation system constructed in the frame consists of two unit direction vectors. It is the vector of constant terms in the linear equation system, representing the position difference vector between the two observation platforms. The distance parameter vector to be solved is as follows: ； in and This represents the scalar distance from the first and second observation platforms to the target point; the distance parameters are solved using the least squares method. : ； Get the first Initial position estimate of the frame Triangulation is performed on all the set initial processing frames to obtain the initial position sequence; S12: Perform a moving average on the initial position sequence obtained in step S11 to obtain the processed position sequence, and take the nearest value. The arithmetic mean of the frames is used as the initial position of the filter. ; S13: Use the position sequence after moving average to determine the target's three coordinates Perform a first-order least-squares linear fit over time. ,in It describes the target position over time. A changing function, and It is the slope and intercept of the first-order minimum fitted line, and the fitted slope is taken. As the initial velocity and with physical limiting; S14: Set the initial triaxial accelerations to 0, i.e., the initial acceleration vector is... , Given the initial velocity obtained in step 3, we obtain the initial state. And establish the initial state covariance in block diagonal form. : ； in , , These are the prior variance matrices for position, velocity, and acceleration, respectively.

3. The method according to claim 2, characterized in that, The uniform velocity CV, uniform acceleration CA, and cooperative turning CT model set in step S2 specifically include: The state transition matrices for uniform velocity CV and uniform acceleration CA are as follows: ； ； in It is the sampling period. It is a 3x3 identity matrix; For any CT plane Its turning angular velocity Instead of preset fixed parameters, the parameters are identified and updated in real time during each filtering cycle based on the online parameter adaptive process in step S3; assuming the turning angular velocity of each plane is... ,make Its velocity rotation matrix and position-velocity coupling submatrix for: The position-velocity coupling submatrix is: ， ； ， ； ， ； Velocity rotation matrix Described in a plane Inside, after time Then, how does the velocity vector rotate; the position-velocity coupling submatrix. Described in a plane Internally, the effect of velocity changes on position; Acting on the corresponding planar velocity components, Acting on the corresponding planar position-velocity coupling, embedded Matrix obtained : 。 4. The method according to claim 3, characterized in that, The measurement model in the parallel extended Kalman filter in step S2 specifically includes: In the implementation of the two observation platforms, for the first One observation platform, Let the target relative quantity be , , Horizontal distance Total distance Take measurement vector azimuth With pitch angle definition: ； Its Jacobian matrix The non-zero elements are only located in the first three columns corresponding to the positional components, satisfying: ； ； For angular residuals Perform surround normalization: ； in, It is the angle residual normalization function, where r is the angle value to be normalized. This indicates the function to round down; then EKF update is performed.

5. The method according to claim 4, characterized in that, Step S2, the adaptive interactive multi-model filtering, specifically includes: Based on the model probability of the previous time step Total number of models and the Markov transition probability matrix between models elements in Calculate the probability of the prediction model and mixture probability : ； ； Based on the posterior state estimates of each model at the previous time step covariance For each model Calculate its mixed initial state and covariance. ; For each model Given a mixed initial state, a state transition model is used for prediction, and EKF updates are performed in parallel to obtain the posterior state of each model. Covariance residual and new information covariance ; According to each model residual New information covariance and measurement vector dimension Calculate the likelihood function value : ； Given the total number of models, the likelihood values ​​are weighted and normalized based on the predicted model probabilities to obtain the posterior model probability at the current time step. : ； The model probability is used both for the final state fusion in the next step and for the control basis for adaptively adjusting the smoothing intensity of the adaptive QUBS in step S3. Posterior state estimation of all filter outputs Based on its latest model probability We perform weighted fusion to obtain the final target state estimate. Covariance : ； 。 6. The method according to claim 5, characterized in that, The online parameter adaptation in step S3 includes: S31: Set the sampling period Segment duration ,make: ； Window length : ； The length is maintained in a first-in, first-out manner. position buffer ; S32: Construct the objective function : ； in For the basis function matrix, It is a second-order difference matrix. These are the control point vectors of the B-spline, and the unknowns to be solved. They are obtained by solving... Determine a smooth curve. It is the adaptive smoothing intensity; S33: Probability from CV model Calculate smoothing strength : ； The upper and lower limits are defined as follows: and exponential smoothing Suppressing oscillations, smoothing coefficient To obtain the adaptive smoothing intensity ; S34: By integer offset Construct multiple sets of three points in a loop. For each set of three points, select the three-point index as... , Using the current time marker, the instantaneous velocity at the endpoint is estimated using second-order backward difference. The signed discrete curvature on each orthogonal plane is calculated using the three-point geometric relationship. Thus, the initial angular velocity is obtained. ; S35: Calculate the arithmetic mean of all the angular velocity values ​​calculated in each group and subject to physical amplitude limiting constraints to obtain the estimated turning angular velocity values ​​in the three planes at the current moment. It is then dynamically written into the state transition matrix of the corresponding CT model for use in the next filtering cycle.

7. A maneuvering target tracking system based on orthogonal decomposition and parameter adaptation, characterized in that, The system for implementing the method according to any one of claims 1 to 6 comprises: The first module is used to obtain the target azimuth and elevation angle measurements through at least two passive observation platforms, obtain an initial three-dimensional position sequence based on multi-frame triangulation and smoothing, and obtain the initial velocity by least squares linear fitting to form the initial state quantity. The second module is used to construct a model set consisting of uniform speed, uniform acceleration, and cooperative turning models corresponding to a horizontal plane and two mutually perpendicular vertical planes, respectively, to describe arbitrary spatial maneuvers and to execute filtering processes such as model interaction and mixing, parallel extended Kalman filtering, model probability update, and state fusion. The third module is used to smooth the trajectory within a fixed-length sliding window using adaptive quasi-uniform B-splines based on the fusion trajectory of the second module, and to calculate the real-time turning angular velocity of each plane in three orthogonal planes using multi-segment three-point geometric curvature. The smoothing intensity of the adaptive quasi-uniform B-splines is adaptively adjusted by the model probability output in step S2, and the identified angular velocities are dynamically written back to the state transition matrix of the corresponding cooperative turning model for use in the next filtering cycle.

8. A computer device, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, When the processor executes the program, it implements the steps of any of the methods described in claims 1-6.

9. A computer-readable storage medium having a computer program stored thereon, characterized in that, When executed by a processor, the program implements the steps of the method described in any one of claims 1-6.

10. A computer program product, comprising a computer program, characterized in that, When executed by a processor, the computer program implements the steps of the method described in any one of claims 1-6.

Images