High-precision satellite passive location and tracking method for moving targets

Abstract

The application provides a high-precision satellite passive positioning and tracking method for a moving target, target radiation source signals are received by at least two satellites, time difference parameters and frequency difference parameters corresponding to target radiation source information of different satellites are calculated by using mutual ambiguity functions; a target position is solved by a three-satellite time difference positioning system or a four-satellite time difference positioning system according to the time difference parameters; a linear equation is constructed to solve a target speed according to the target position, the frequency difference parameters and satellite speed information; a target three-dimensional motion state is decoupled into XYZ three-axis independent models, and filtering optimization is performed by using an extended Kalman filtering algorithm based on the target position and the target speed, so that a target positioning and tracking result is obtained. In this way, the application can obtain a higher-precision target motion state while reducing the calculation amount, and realize continuous and stable tracking of the target.

Description

Technical Field

[0001] This invention relates to the field of positioning technology, and in particular to a high-precision passive satellite positioning and tracking method for moving targets. Background Technology

[0002] Currently, there are two main types of target localization methods: active localization and passive localization. Active localization involves sensors actively emitting signals and measuring the reflected signals to locate the target. Passive localization, on the other hand, involves sensor nodes that do not emit signals but instead determine the location of radiation sources by receiving signals from unknown sources. For moving targets, using satellites equipped with sensors to continuously observe and receive radiation source signals over a wide area, and then accurately locate and track them, has broad application prospects in both military and civilian fields, and is a very important research area within passive localization.

[0003] Current research mainly focuses on improving the accuracy of single-shot positioning. However, from the perspective of the entire process of target movement, the motion state of a target at different times is interconnected. Introducing historical observation information can effectively improve the accuracy of estimating the target's motion state, but there is relatively little research on this topic in the field of passive positioning. Summary of the Invention

[0004] The purpose of this invention is to provide a high-precision passive satellite positioning and tracking method for moving targets, which can obtain higher accuracy target motion status while reducing computational load, and achieve continuous and stable tracking of the target.

[0005] To achieve the above objectives, this invention provides a high-precision passive satellite positioning and tracking method for moving targets, comprising the following steps:

[0006] By receiving the target radiation source signal from at least two satellites, the time difference parameter and frequency difference parameter between the target radiation source signals from different satellites are calculated using a mutual ambiguity function;

[0007] Based on the time difference parameters, the target location is calculated using a three-star time difference positioning system or a four-star time difference positioning system.

[0008] Based on the target location, the frequency difference parameter, and the satellite velocity information, a linear equation is constructed to solve for the target velocity;

[0009] The target's three-dimensional motion state is decoupled into an independent XYZ three-axis model, and based on the target's position and velocity, the extended Kalman filter algorithm is used for filtering optimization to obtain the target positioning and tracking results.

[0010] Optionally, the formula for calculating the mutual ambiguity function is:

[0011]

[0012] Wherein, s1(t) and s2(t) are the target radiation source signals received by the two satellites respectively; τ and f are the time difference parameter and the frequency difference parameter corresponding to the maximum value determined by two-dimensional search.

[0013] Optionally, the target location can be calculated using a three-star or four-star time difference positioning system based on the time difference parameters, including:

[0014] When the number of satellites receiving signals is 3, the first set of simultaneous equations is obtained by combining the two independent time difference equations of the satellites and the equation of the Earth's sphere based on the three-satellite time difference positioning system.

[0015] When the number of satellites receiving signals is greater than or equal to 4, the second set of simultaneous equations is obtained by combining the three independent time difference equations of the four-satellite time difference positioning system.

[0016] Configure the initial parameters for iteration, and iteratively solve the first or second simultaneous equation system using the Levenberg-Marquardt method to calculate the target position.

[0017] Optionally, the first system of simultaneous equations is:

[0018]

[0019] The second system of simultaneous equations is:

[0020]

[0021] The positions of each satellite are (x i y i , z i ), i is the satellite marker, i = 1, 2, 3, 4; the location of the radiation source is (x, y, z), c is the speed of light, R e Let τ be the radius of the Earth's sphere. 12 τ is the time difference between the arrival of the signal from the unknown radiation source on satellite 1 and satellite 2. 13 τ is the time difference between the arrival of the signal from the unknown radiation source on satellite 1 and satellite 3. 14 The time difference between the arrival of the signal from the unknown radiation source on satellite 1 and satellite 4.

[0022] Optionally, the linear equation for the target velocity, based on the three-star time difference positioning system, is as follows:

[0023]

[0024] The linear equation for the target velocity, derived from the four-satellite time difference positioning system, is as follows:

[0025]

[0026] in:

[0027]

[0028]

[0029] (v x v y v z Let f be the velocity of the radiation source along the X, Y, and Z axes. d12 f is the frequency difference between the arrival signals from the unknown radiation source on satellite 1 and satellite 2. d13 f represents the frequency difference between the arrival signals from the unknown radiation source on satellite 1 and satellite 3. d14 The frequency difference between satellite 1 and satellite 4 is the signal from the unknown radiation source; the velocities of each satellite on the X, Y, and Z axes are (v... xi v yi v zi ), i is the satellite marker, i = 1, 2, 3, 4; the location of the radiation source is (x, y, z), and c is the speed of light.

[0030] Optionally, the XYZ three-axis independent model is a motion state change model, and the expression of the motion state change model is:

[0031] s x [k]=Φ[k,k-1]s x [k-1]+Γ[k-1]n x [k-1];

[0032] s y [k]=Φ[k,k-1]s y [k-1]+Γ[k-1]n y [k-1];

[0033] s z [k]=Φ[k,k-1]s z [k-1]+Γ[k-1]n z [k-1];

[0034] Where Φ is the state transition matrix of the extended Kalman filter, and Γ is the perturbation matrix of the extended Kalman filter;

[0035]

[0036] Δt=t k -t k-1 t represents the time interval between adjacent observations. k For the current time, t k-1 n is the time of the last observation. x[k-1] is t k-1 System disturbance noise along the target's X-axis at any given time, n y [k-1] is t k-1 System disturbance noise along the target Y-axis at any given time, n z [k-1] is t k-1 System disturbance noise along the Z-axis of the target at any given time.

[0037] Optionally, the initial state filter value of the extended Kalman filter algorithm is set as follows:

[0038]

[0039] Optionally, the step of optimizing the filter using the extended Kalman filter algorithm includes:

[0040] Update the filter parameters;

[0041] Establish state prediction equations, error covariance prediction equations, and measurement equations based on the radiation source motion state change model;

[0042] Calculate the Kalman filter gain based on the state prediction equation, error covariance prediction equation, and measurement equation.

[0043] Based on the Kalman filter gain, the current state filter value is obtained;

[0044] When a new target radiation source signal is received, the new target position and velocity are calculated through steps of time difference / frequency difference calculation, target position calculation, and target velocity calculation. The filter parameters for the new time period are then updated based on the calculated new target position and velocity. Optionally, the state prediction equation is:

[0045]

[0046] in, For t k The predicted state value at time 10:00. For t k-1 The state filter value at time t, Φ[k,k-1] is the value corresponding to t. k-1 Time and t k The state transition matrix at time t;

[0047] The error covariance prediction equation is as follows:

[0048] P'[k / k-1]=Φ[k,k-1]P[k-1]Φ T [k,k-1]+Γ[k-1]Q x / y / z Γ T [k-1];

[0049] Where P′[k / k-1] is the value corresponding to t k-1 Time and t k The one-step prediction mean square error matrix at time t, where P[k-1] represents the mean square error at time t. k-1 The mean square error matrix at time t, where Γ[k-1] represents the mean square error at time t. k-1 The system perturbation matrix at time Q x / y / z For t k-1 System noise at any given moment;

[0050] t k The measurement equation for time is:

[0051]

[0052] Optionally, the Kalman filter gain is calculated based on the following formula:

[0053] K[k]=P'[k / k-1]H T [k](H[k]P'[k / k-1]H T [k]+R x / y / z [k]) -1 ;

[0054] Where P′[k / k-1] is the value corresponding to t k-1 Time and t k The one-step prediction mean square error matrix at time t, R x / y / z For t k Observation noise at any given moment;

[0055] The current state filter value is calculated based on the following formula:

[0056]

[0057] in, For t k The predicted state value at time t, K[k] is the state prediction value at time t. k Kalman filter gain at time O x / y / z [k] is t k The observed value at that moment.

[0058] This invention utilizes an extended Kalman filter algorithm, iteratively optimizing the results using observation data from different time points. The resulting position is more accurate than a single positioning test, and it also obtains higher-order motion state information such as velocity and acceleration, leading to better estimation of the target's motion state. Furthermore, based on the calculated target position, this invention introduces frequency difference information to construct a linear equation system to solve for the target velocity. This results in faster calculation speed and provides a high-dimensional input for the extended Kalman filter algorithm, enhancing its performance. In addition, this invention innovatively decouples the XYZ dimensions for separate filtering during Kalman filtering. Firstly, since the target's three-dimensional motion intensity inevitably varies, three-dimensional decoupling allows for setting more suitable parameters based on the actual situation. Secondly, after decoupling, the required order of the computation matrix in each dimension is reduced, while the three dimensions can be calculated synchronously using parallel computing methods, significantly improving overall computational efficiency. Attached Figure Description

[0059] Figure 1 This is a flowchart illustrating the steps of a high-precision satellite passive positioning and tracking method for moving targets according to an embodiment of the present invention.

[0060] Figure 2 This is a schematic diagram illustrating the specific execution flow of the high-precision satellite passive positioning and tracking method for moving targets provided in an embodiment of the present invention. Detailed Implementation

[0061] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0062] It should be noted that references to "an embodiment," "embodiment," "example embodiment," etc., in this specification refer to the described embodiment including specific features, structures, or characteristics, but not every embodiment must include these specific features, structures, or characteristics. Furthermore, such expressions do not refer to the same embodiment. Moreover, when describing specific features, structures, or characteristics in conjunction with embodiments, whether or not explicitly described, it is indicated that incorporating such features, structures, or characteristics into other embodiments is within the knowledge of those skilled in the art.

[0063] Furthermore, certain terms are used in the specification and subsequent claims to refer to specific components or parts. Those skilled in the art will understand that manufacturers may use different names or terms to refer to the same component or part. This specification and subsequent claims do not distinguish components or parts by differences in name, but rather by differences in function. The terms "comprising" and "including" used throughout the specification and subsequent claims are open-ended and should be interpreted as "including but not limited to." Additionally, the term "connection" here includes any direct and indirect electrical connection means. Indirect electrical connection means include connections made through other means.

[0064] Before describing the embodiments of this application in detail, the technical concept of this application is briefly described first: This invention solves the problems of low single-shot accuracy and discontinuous motion state estimation in traditional passive positioning by using joint time difference and frequency difference measurement and iterative optimization of historical data. First, the time difference and frequency difference parameters are extracted using the mutual ambiguity function of the satellite received signals; based on the time difference information, the target position is calculated using a three-star or four-star positioning system; combining the frequency difference parameters and satellite velocity information, a linear equation is constructed to directly solve for the target velocity; finally, the target's three-dimensional motion state is decoupled into an independent XYZ axis model, and the position, velocity, and acceleration are optimized by parallel filtering using an extended Kalman filter (EKF). By fusing historical observation data and real-time calculation results, the target motion parameters are iteratively corrected, significantly improving positioning accuracy and tracking stability while reducing computational complexity.

[0065] The specific principles of the high-precision satellite passive positioning and tracking method for moving targets of this application will be described below with reference to specific embodiments.

[0066] Figures 1-2 This invention illustrates a high-precision satellite passive positioning and tracking method for moving targets according to an embodiment of the present invention, comprising the following steps:

[0067] S101: Receive the target radiation source signal through at least two satellites, and calculate the time difference parameter and frequency difference parameter between the target radiation source signals from different satellites using a mutual ambiguity function.

[0068] In practice, at the current moment, the signals from the target radiation source are received by two satellites and denoted as s1(t) and s2(t). The time difference and frequency difference at this moment are measured using a mutual ambiguity function.

[0069] The formula for calculating the mutual ambiguity function is as follows:

[0070]

[0071] Wherein, s1(t) and s2(t) are the target radiation source signals received by the two satellites respectively; τ and f are the time difference parameter and the frequency difference parameter corresponding to the maximum value determined by two-dimensional search; that is, τ and f are the time difference and frequency difference values ​​to be estimated; * indicates conjugate; by adjusting different values ​​of τ and f to perform two-dimensional search, the τ and f values ​​corresponding to the maximum value of the function are the time difference parameter and frequency difference parameter of the signal.

[0072] S102: Based on the time difference parameters, the target position is calculated using a three-satellite time difference positioning system or a four-satellite time difference positioning system. Step S102 in this embodiment includes: when the number of satellites receiving signals is 3, a first set of simultaneous equations is obtained by combining two independent time difference equations of the satellites and the Earth's sphere equation based on the three-satellite time difference positioning system; when the number of satellites receiving signals is greater than or equal to 4, a second set of simultaneous equations is obtained by combining three independent time difference equations of the satellites based on the four-satellite time difference positioning system; initial iteration parameters are configured, and the first or second set of simultaneous equations is iteratively solved using the Levenberg-Marquardt method to calculate the target position.

[0073] In specific implementation, this embodiment uses the Levenberg-Marquardt method to calculate the target position, which is achieved through the following sub-steps:

[0074] (2-1) Determine the positioning system and, based on this, clarify the nonlinear problem model.

[0075] After obtaining the time difference values ​​between multiple satellites at the current moment, the target position can be initially calculated based on the number of satellites receiving the signal and the selected positioning system. Specifically, when the number of satellites is less than 3, it is impossible to locate the moving target; when the number of satellites is 3, the target position can be calculated using a three-satellite time difference positioning system; when the number of satellites is 4, the target position can be calculated using a four-satellite time difference positioning system; and when the number of satellites is greater than 4, four satellites can be selected to use a four-satellite time difference positioning system to calculate the target position.

[0076] The following will explain the three-star or four-star time difference positioning system used.

[0077] The three-satellite time difference positioning system uses two independent time difference equations for the arrival of radiation source signals from three satellites, combined with the equations for the sphere determined by the Earth's surface, to solve for the coordinates of the radiation source. The first set of equations is as follows:

[0078]

[0079] If there are three satellites receiving the signal, denoted as Satellite 1, Satellite 2, and Satellite 3, then the coordinates of the radiation source can be solved by combining the independent time difference equations corresponding to Satellite 1 with the equations of the sphere determined by the Earth's surface.

[0080] The four-satellite time difference positioning system uses three independent time difference equations for the arrival of the radiation source signal at four satellites to solve for the coordinates of the radiation source. The second simultaneous equation system is as follows:

[0081]

[0082] The positions of each satellite are (x i y i , z i ), i is the satellite marker, i = 1, 2, 3, 4; the location of the radiation source is (x, y, z), c is the speed of light, R e Let τ be the radius of the Earth's sphere. 12 τ is the time difference between the arrival of the signal from the unknown radiation source on satellite 1 and satellite 2. 13 τ is the time difference between the arrival of the signal from the unknown radiation source on satellite 1 and satellite 3. 14 The time difference between the arrival of the signal from the unknown radiation source on satellite 1 and satellite 4.

[0083] If there are four or more satellites receiving signals, the four satellites with the best geometric distribution can be selected first for positioning calculation to reduce positioning error; these are denoted as satellite 1, satellite 2, satellite 3 and satellite 4. The independent time difference equations corresponding to satellite 1 and the other three satellites (satellite 2, satellite 3 and satellite 4) are then combined to solve for the coordinates of the radiation source.

[0084] The positions of each satellite are (x i y i , z i ), i is the satellite marker, i = 1, 2, 3, 4; the location of the radiation source is (x, y, z), c is the speed of light, R e Let τ be the radius of the Earth's sphere. 12 τ is the time difference between the arrival of the signal from the unknown radiation source on satellite 1 and satellite 2. 13 τ is the time difference between the arrival of the signal from the unknown radiation source on satellite 1 and satellite 3. 14 τ is the time difference between the arrival of the signal from the unknown radiation source on satellite 1 and satellite 4. 12 τ 13 and τ 14 Solve using the above steps S101.

[0085] Since the positioning equation is a nonlinear equation, it can be solved using common nonlinear equation solution methods. Therefore, it is necessary to define the corresponding Jacobian matrix and objective function.

[0086] For Samsung's time zone positioning system, its Jacobian matrix J(x, y, z) and objective function F(x, y, z) are:

[0087]

[0088] For the four-star time difference positioning system, its Jacobian matrix J(x, y, z) and objective function F(x, y, z) are:

[0089]

[0090] in:

[0091]

[0092] (2-2) Configure the initial parameters for iteration.

[0093] Let the iteration number be j, and select the initial iteration point (x) within the satellite coverage area based on engineering experience. j ,y j ,z j ), j = 0, and take parameters ε, μ0, γ1, γ2, λ1, λ2; and make such that

[0094] 0 < ε < 1

[0095] μ0>0

[0096] 0 < γ1 < 1 < γ2

[0097] 0 < λ1 < λ2 ≤ 1.

[0098] (2-3) Calculate the Jacobian matrix and objective function corresponding to this iteration.

[0099] That is, calculate the current iteration (x) j ,y j ,z j The corresponding Jacobian matrix J(x) j ,y j ,z j ) and objective function F(x) j ,y j ,z j ).

[0100] (2-4) Determine whether the iteration termination condition is met. If so, output the target position of this solution.

[0101] Specifically, determine whether the following inequalities are true. If they are true, end the iteration and output the value of (x) at this point. j ,y j ,z j If the calculated target position is not specified, proceed to the next step.

[0102] ||J T (x j ,y j ,z j )F(x j ,y j ,z j ||≤ε;

[0103] For the sake of simplicity, let J be the following. j =J(x) j ,y j ,z j ), F j =F(x) j ,y j ,z j ).

[0104] (2-5) Update the relevant iteration parameters.

[0105] Calculate vector d j :

[0106]

[0107] Here, matrix I is the identity matrix.

[0108] Calculate the numerical value ρ j :

[0109]

[0110] Among them, let (x j ',y j ',z j ') T =(x j ,y j ,z j ) T +d j F j '=F(x j ',y j ',z j ').

[0111] Logarithmic value ρ j Perform interval judgment, when ρ j <λ1, then μ j+1 =γ2μ j When λ1≤ρ j If λ < 2, then μ j+1 =μ j When λ2≤ρ j Then we have μ j+1 =γ1μ j ;

[0112] Take (x) j+1 ,y j+1 ,z j+1 )=(x j ',y j ',z j '), increment the iteration number j by 1, and return to sub-step 3 to calculate the Jacobian matrix and objective function corresponding to this iteration.

[0113] S103: Based on the target position, the frequency difference parameter, and the satellite velocity information, a linear equation is constructed to solve for the target velocity. Once the target position at the current moment is calculated, it is combined with the known satellite position and velocity information and the frequency difference parameter to calculate the target velocity.

[0114] Once the target position is obtained using the Samsung time difference positioning system, the target velocity can be represented as a matrix form as follows:

[0115]

[0116] Once the target's position is obtained using a four-satellite time difference positioning system, the target's velocity can be represented in the following matrix form:

[0117]

[0118] Both of the above matrix representations are linear, so the target speed can be quickly obtained by inverting the matrix.

[0119] in:

[0120]

[0121]

[0122] (v x v y v z Let f be the velocity of the radiation source along the X, Y, and Z axes. d12 f is the frequency difference between the arrival signals from the unknown radiation source on satellite 1 and satellite 2. d13 f represents the frequency difference between the arrival signals from the unknown radiation source on satellite 1 and satellite 3. d14 f represents the frequency difference between the arrival signals from the unknown radiation source on satellite 1 and satellite 4. d12 f d13 and f d14 The solution is obtained through step S101 above. The velocities of each satellite on the X, Y, and Z axes are (v... xi v yi v zi ), i is the satellite marker, i = 1, 2, 3, 4; the location of the radiation source is (x, y, z), and c is the speed of light.

[0123] S104: Decouple the target's three-dimensional motion state into an independent XYZ three-axis model, and based on the target's position and velocity, use the extended Kalman filter algorithm for filtering optimization to obtain the target positioning and tracking results.

[0124] After obtaining the target's current position and velocity information, a tracking algorithm can be used to obtain more accurate target position and velocity information. This embodiment decouples the XYZ dimensions based on the extended Kalman filter algorithm framework, performing filtering in each dimension separately. This reduces the algorithm's size and allows for simultaneous processing via parallel computing, effectively improving computational efficiency.

[0125] See Figure 2 The specific implementation steps of step S104 are as follows:

[0126] (4-1) Establish a model of the change in the motion state of the radiation source.

[0127] The XYZ three-axis independent model is a motion state change model. In this embodiment, the motion state change models are established for the XYZ axes of the target radiation source as follows:

[0128] s x [k]=Φ[k,k-1]s x [k-1]+Γ[k-1]n x [k-1];

[0129] s y [k]=Φ[k,k-1]s y [k-1]+Γ[k-1]n y [k-1];

[0130] s z [k]=Φ[k,k-1]s z [k-1]+Γ[k-1]n z [k-1];

[0131] Where Φ is the state transition matrix of the extended Kalman filter, and Γ is the perturbation matrix of the extended Kalman filter;

[0132]

[0133] Δt=t k -t k-1 t represents the time interval between adjacent observations. k For the current time, t k-1 n is the time of the last observation. x [k-1] is t k-1 System disturbance noise along the target's X-axis at any given time, n y[k-1] is t k-1 System disturbance noise along the target Y-axis at any given time, n z [k-1] is t k-1 System disturbance noise along the Z-axis of the target at any given time.

[0134] s x [k] is t k The state matrix of the target along the X-axis at time s x [k-1] is t k-1 The state matrix of the target along the X-axis at time s y [k] is t k The state matrix of the target along the Y-axis at time s y [k-1] is t k-1 The state matrix of the target along the Y-axis at time s z [k] is t k The state matrix of the target along the Z-axis at time s z [k-1] is t k-1 The state matrix of the target Z-axis at time t, Φ[k,k-1] is t k-1 Time and t k The state transition matrix between time steps, Γ[k-1] is the state transition matrix between time steps t. k-1 The system perturbation matrix at time t.

[0135] Specifically, s x [k]=[x(t k )v x (t k )a x (t k )] T ;

[0136] s x [k-1]=[x(t k-1 )v x (t k-1 )a x (t k-1 )] T ;

[0137] s y [k]=[y(t k )v y (t k )a y (t k )] T ;

[0138] s y [k-1]=[y(t k-1 )v y (t k-1 )a y (t k-1)] T ;

[0139] s z [k]=[z(t k )v z (t k )a z (t k )] T ;

[0140] s z [k-1]=[z(t k-1 )v z (t k-1 )a z (t k-1 )] T ;

[0141]

[0142] n x [k-1]=b x (t k-1 );

[0143] n y [k-1]=b y (t k-1 );

[0144] n z [k-1]=b z (t k-1 );

[0145] Where x(t) k ) for t k The position of the target on the X-axis at any given time, v x (t k ) for t k The velocity of the target along the X-axis at any given time, a x (t k ) for t k The acceleration of the target along the X-axis at any given time, b x (t k ) for t k The second-order acceleration of the target along the X-axis at time t; y(t) k ) for t k The position of the target on the Y-axis at any given time, v y (t k ) for t k The velocity of the target along the Y-axis at any given time, a x (t k ) for t k The target's acceleration along the Y-axis at any given time, b y (t k ) for t kThe second-order acceleration of the target along the Y-axis at time t; z(t) k ) for t k The position of the target on the Z-axis at any given time, v z (t k ) for t k The velocity of the target along the Z-axis at any given time, a z (t k ) for t k The target's acceleration along the Z-axis at any given time, b z (t k ) for t k The second-order acceleration of the target along the Z-axis at any given time.

[0146] Since the XYZ axes are completely dual, for the sake of simplicity, the XYZ axes will be described uniformly below, that is, the XYZ axes in t k The state matrix at time s is represented as s x / y / z [k], the system disturbance noise is characterized as n x / y / z [k], the state filter value is characterized as State prediction value is characterized as The observed value is represented as O x / y / z [k], system noise is characterized as O x / y / z The observation noise is characterized as R x / y / z .

[0147] (4-2) Set the initial state.

[0148] The time marker number is k=1. In the first calculation, the target position (x, y, z) and velocity (v) are obtained. x v y v z After that, the initial state filter value of the extended Kalman filter algorithm is set as follows:

[0149]

[0150] Meanwhile, based on engineering experience, the initial estimated mean square error matrix P[1] and the noise figure Q of the XYZ triaxial system are set. x / Q y / Q z XYZ triaxial observation noise parameter R x / R y / R z .

[0151] (4-3) Perform filtering at the new moment to obtain high-precision motion state parameters.

[0152] After receiving the signal at the new time, and after calculating the time difference / frequency difference, target position, and target velocity as described above, new target position and velocity information is obtained. A filtering algorithm can then be used to obtain more accurate position and velocity values. This embodiment uses a Kalman filter algorithm framework for processing, and the specific implementation steps include:

[0153] 1. Update the filter parameters.

[0154] Increment the time stamp number k by one, and record the new signal reception time as t. k The last time the signal was received is denoted as t. k-1 Time, thus obtaining t k-1 Time and t k The state transition matrix Φ[k,k-1] between time steps t k-1 The system perturbation matrix Γ[k-1] at time t, and the target position (x, y, z) and velocity (v) obtained in this calculation. x v y v z Update XYZ triaxial observations O x / y / z [k]:

[0155] O x [k]=[xv x 0] T ;

[0156] O y [k]=[yv y 0] T ;

[0157] O z [k]=[zv z 0] T .

[0158] 2. Establish state prediction equations, error covariance prediction equations, and measurement equations based on the radiation source motion state change model.

[0159] Specifically, the state prediction equation is:

[0160]

[0161] in, For t k The predicted state value at time 10:00. For t k-1 The state filter value at time t, Φ[k,k-1] is the value corresponding to t. k-1 Time and t k The state transition matrix at time t;

[0162] The error covariance prediction equation is as follows:

[0163] P'[k / k-1]=Φ[k,k-1]P[k-1]Φ T [k,k-1]+Γ[k-1]Q x / y / z Γ T [k-1];

[0164] Where P′[k / k-1] is the value corresponding to t k-1 Time and t k The one-step prediction mean square error matrix at time t, where P[k-1] represents the mean square error at time t. k-1 The mean square error matrix at time t, where Γ[k-1] represents the mean square error at time t. k-1 The system perturbation matrix at time Q x / y / z For t k-1 System noise at any given moment;

[0165] t k The measurement equation for time is:

[0166]

[0167] 3. Calculate the Kalman filter gain based on the state prediction equation, error covariance prediction equation, and measurement equation.

[0168] t k The Kalman filter gain at time t is calculated based on the following formula:

[0169] K[k]=P'[k / k-1]H T [k](H[k]P'[k / k-1]H T [k]+R x / y / z [k]) -1 ;

[0170] Where P′[k / k-1] is the value corresponding to t k-1 Time and t k The one-step prediction mean square error matrix at time t, R x / y / z For t k Observational noise at any given moment.

[0171] 4. Obtain the current state filter value based on the Kalman filter gain.

[0172] Current time t k The state filter value is calculated based on the following formula:

[0173]

[0174] in, For t k The predicted state value at time t, K[k] is the state prediction value at time t. k Kalman filter gain at time Ox / y / z [k] is t k The observed value at that moment.

[0175] Simultaneously calculate the estimated mean square error matrix at the current time:

[0176] P[k]=(IK[k]H[k])P'[k / k-1];

[0177] Where P[k] is t k The mean square error matrix for time t, where I is the identity matrix and K[k] is the mean square error matrix for time t. k The Kalman filter gain at time t, H[k] is the Kalman filter gain at time t. k The measurement matrix at time t, P′[k / k-1] is the matrix corresponding to time t. k-1 Time and t k The one-step prediction mean square error matrix at time t.

[0178] State filter value The first element corresponds to the target's position on the X-axis, the second element corresponds to the target's velocity on the X-axis, and the third element corresponds to the target's acceleration on the X-axis; state filter value. The first element corresponds to the target's position on the Y-axis, the second element corresponds to the target's velocity on the Y-axis, and the third element corresponds to the target's acceleration on the Y-axis; state filter value. The first element corresponds to the target's position on the Z-axis, the second element corresponds to the target's velocity on the Z-axis, and the third element corresponds to the target's acceleration on the Z-axis; by outputting the above parameters, high-precision motion state information of the target can be obtained.

[0179] 5. When a new target radiation source signal is received, the new target position and new target velocity are calculated through the steps of time difference frequency difference calculation, target position calculation, and target velocity calculation (i.e., steps S101, S102, and S103), and the filter parameters at the new moment are updated based on the calculated new target position and new target velocity (i.e., return to step 1 to update the motion state information at the new moment).

[0180] In summary, this invention utilizes the extended Kalman filter algorithm and iteratively optimizes the results using observation data from different time points, achieving higher accuracy compared to single-shot positioning. It also obtains higher-order motion state information such as velocity and acceleration, resulting in better estimation of the target's motion state. Furthermore, based on the calculated target position, this invention introduces frequency difference information to construct a linear equation system to solve for the target velocity, resulting in faster calculation speed and providing a high-dimensional input for the extended Kalman filter algorithm, further enhancing its performance. Moreover, this invention innovatively decouples the XYZ dimensions for separate filtering during Kalman filtering. This allows for more tailored parameter settings based on the specific circumstances, given the inherent differences in the target's three-dimensional motion intensity. Additionally, decoupling reduces the order of the computational matrix required in each dimension, while allowing for simultaneous parallel computation in the three dimensions, significantly improving overall computational efficiency.

[0181] It should be noted that this application can be implemented in software and / or a combination of software and hardware, for example, using an application-specific integrated circuit (ASIC), a general-purpose computer, or any other similar hardware device. In one embodiment, the software program of this application can be executed by a processor to implement the steps or functions described above. Similarly, the software program of this application (including related data structures) can be stored in a computer-readable recording medium, such as RAM memory, magnetic or optical drives, floppy disks, and similar devices. Furthermore, some steps or functions of this application can be implemented in hardware, for example, as circuitry that works with a processor to perform the various steps or functions.

[0182] The method according to the invention can be implemented on a computer as a computer-implemented method, or in dedicated hardware, or a combination of both. Executable code or portions thereof for the method according to the invention can be stored on a computer program product. Examples of computer program products include memory devices, optical storage devices, integrated circuits, servers, online software, etc. Preferably, the computer program product includes non-transitory program code components stored on a computer-readable medium so as to execute the method according to the invention when the program product is executed on a computer.

[0183] Of course, the present invention may have other various embodiments. Without departing from the spirit and essence of the present invention, those skilled in the art can make various corresponding changes and modifications according to the present invention, but these corresponding changes and modifications should all fall within the protection scope of the appended claims.

Claims

1. A high-precision satellite passive positioning and tracking method for moving targets, characterized in that, Including the following steps: By receiving the target radiation source signal from at least two satellites, the time difference parameter and frequency difference parameter between the target radiation source signals from different satellites are calculated using a mutual ambiguity function; Based on the time difference parameters, the target location is calculated using a three-star time difference positioning system or a four-star time difference positioning system. Based on the target location, the frequency difference parameter, and the satellite velocity information, a linear equation is constructed to solve for the target velocity; The target's three-dimensional motion state is decoupled into an independent XYZ three-axis model, and based on the target's position and velocity, the extended Kalman filter algorithm is used for filtering optimization to obtain the target positioning and tracking results. The linear equation for the target velocity, derived using the Samsung time difference positioning system, is as follows: ； The linear equation for the target velocity, derived from the four-satellite time difference positioning system, is as follows: ； in: ； ； ( , , (x) represents the velocity of the radiation source along the X, Y, and Z axes. The frequency difference between the arrival signals from the unknown radiation source on satellite 1 and satellite 2. The frequency difference between the arrival signals from the unknown radiation source on satellite 1 and satellite 3. The frequency difference between satellite 1 and satellite 4 is the signal from the unknown radiation source; the velocities of each satellite on the XYZ axes are ( ), i Mark the satellite. i =1, 2, 3, 4; the location of the radiation source is ( x , y , z ), c The speed of light; The XYZ three-axis independent model is a motion state change model, and the expression of the motion state change model is: ； in, To extend the state transition matrix of the Kalman filter, For corresponding Time and The state transition matrix at time t, for The state matrix of the target along the X-axis at any given time. for The state matrix of the target along the Y-axis at any given time. for The state matrix of the target along the Z-axis at any given time. for The system perturbation matrix at time t. The perturbation matrix for the extended Kalman filter; ； ； The time interval between adjacent observations For the current moment, This refers to the last observation time; for System disturbance noise along the target's X-axis at any given time. for System disturbance noise along the target Y-axis at any given time. for System disturbance noise along the Z-axis of the target at any given time; The initial state filter value of the extended Kalman filter algorithm is set as follows: ； in( , , ) represents the velocity of the radiation source along the XYZ axes.

2. The high-precision satellite passive positioning and tracking method for moving targets according to claim 1, characterized in that, The formula for calculating the mutual ambiguity function is as follows: ； in, and The signals from the target radiation source received by the two satellites at time t respectively; and To determine the time difference parameter and frequency difference parameter corresponding to the maximum value through a two-dimensional search; * represents conjugate.

3. The high-precision satellite passive positioning and tracking method for moving targets according to claim 1, characterized in that, Based on the time difference parameters, the target location is calculated using a three-star or four-star time difference positioning system, including: When the number of satellites receiving signals is 3, the first set of simultaneous equations is obtained by combining the two independent time difference equations of the satellites and the equation of the Earth's sphere based on the three-satellite time difference positioning system. When the number of satellites receiving signals is greater than or equal to 4, the second set of simultaneous equations is obtained by combining the three independent time difference equations of the four-satellite time difference positioning system. Configure the initial parameters for iteration, and iteratively solve the first or second simultaneous equation system using the Levenberg-Marquardt method to calculate the target position.

4. The high-precision satellite passive positioning and tracking method for moving targets according to claim 3, characterized in that, The first system of simultaneous equations is: ； The second system of simultaneous equations is: ； The positions of each satellite are ( ), i Mark the satellite. i =1, 2, 3, 4; the location of the radiation source is ( x , y , z ), c At the speed of light, The radius of the Earth's sphere. The time difference between the arrival of the signal from the unknown radiation source on satellite 1 and satellite 2. The time difference between the arrival of the signal from the unknown radiation source on satellite 1 and satellite 3. The time difference between the arrival of the signal from the unknown radiation source on satellite 1 and satellite 4.

5. The high-precision satellite passive positioning and tracking method for moving targets according to claim 1, characterized in that, The steps for filtering optimization using the extended Kalman filter algorithm include: Update the filter parameters; Establish state prediction equations, error covariance prediction equations, and measurement equations based on the radiation source motion state change model; Calculate the Kalman filter gain based on the state prediction equation, error covariance prediction equation, and measurement equation. Based on the Kalman filter gain, the current state filter value is obtained; When a new target radiation source signal is received, the new target position and new target velocity are calculated through the steps of time difference frequency difference calculation, target position calculation, and target velocity calculation, and the filter parameters at the new moment are updated based on the calculated new target position and new target velocity.

6. The high-precision satellite passive positioning and tracking method for moving targets according to claim 5, characterized in that, The state prediction equation is: ； in, for The predicted state value at time 10:

00. for The state filter value at time 10:

00. For corresponding Time and The state transition matrix at time t; The error covariance prediction equation is as follows: ； in, For corresponding Time and The one-step prediction mean square error matrix at time 1. for The mean square error matrix for time step 1. for The system perturbation matrix at time t. for System noise at any given moment; The measurement equation for time is: 。 7. The high-precision satellite passive positioning and tracking method for moving targets according to claim 6, characterized in that, The Kalman filter gain is calculated based on the following formula: ； in, For corresponding Time and The one-step prediction mean square error matrix at time 1. for Observation noise at any given moment; The current state filter value is calculated based on the following formula: ； in, for The predicted state value at time 10:

00. K [ k ]for Kalman filter gain at time t. for The observed value at that moment.

Images