High-precision satellite passive positioning and tracking method for moving target

Abstract

The invention provides a high-precision satellite passive positioning and tracking method for a moving target, and the method comprises the steps: receiving a target radiation source signal through at least two satellites, and calculating time difference parameters and frequency difference parameters corresponding to the target radiation source information of different satellites through a cross-ambiguity function; according to the time difference parameters, the target position is solved through a three-satellite time difference positioning system or a four-satellite time difference positioning system; according to the target position, the frequency difference parameter and the satellite speed information, constructing a linear equation to solve the target speed; and decoupling the three-dimensional motion state of the target into an XYZ three-axis independent model, and carrying out filtering optimization by using an extended Kalman filtering algorithm based on the target position and the target speed so as to obtain a target positioning and tracking result. Therefore, the target motion state with higher precision can be obtained while the calculation amount is reduced, and continuous and stable tracking of the target is realized.

Description

Technical Field

[0001] The present invention relates to the technical field of positioning, and in particular, to a high-precision satellite passive positioning and tracking method for moving targets. Background Art

[0002] Currently, there are two types of target positioning methods: active positioning and passive positioning. Active positioning is that the sensor actively emits signals outward and completes the positioning of the target by measuring the reflected signals. Passive positioning, the sensor node itself does not emit signals outward, but determines the position of the radiation source by receiving signals from unknown radiation sources. For moving targets, using satellites carrying sensors to continuously observe and receive radiation source signals over a large area and accurately position and track them has broad application prospects in both military and civilian fields, and is a very important research content in the field of passive positioning.

[0003] The current existing research mainly focuses on improving the single positioning accuracy. However, from the analysis of the entire process of target movement, the movement states of a target between consecutive moments are interrelated. If the historical information of observations is introduced, the accuracy of estimating the target movement state can be effectively improved, but there is relatively little relevant research in the field of passive positioning. Summary of the Invention

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

[0005] To achieve the above purpose, the present invention provides a high-precision satellite passive positioning and tracking method for moving targets, including the steps of:

[0006] Receiving target radiation source signals through at least two satellites, and calculating the corresponding time difference parameter and frequency difference parameter between the target radiation source signals of different satellites by using the cross ambiguity function;

[0007] According to the time difference parameter, solving the target position by using a three-satellite time difference positioning system or a four-satellite time difference positioning system;

[0008] According to the target position, the frequency difference parameter and the satellite speed information, constructing a linear equation to solve the target speed;

[0009] Decoupling the target three-dimensional motion state into an XYZ three-axis independent model, and based on the target position and the target speed, using the extended Kalman filter algorithm for filtering optimization to obtain the target positioning and tracking result.

[0010] Optionally, the calculation formula of the cross ambiguity function is:

[0011]

[0012] Among them, s1(t) and s2(t) are the target radiation source signals received by 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, according to the time difference parameter, calculating the target position by using a three-star time difference positioning system or a four-star time difference positioning system includes:

[0014] When the number of satellites receiving signals is 3, the two independent time difference equations of the satellites and the earth's spherical surface equation are combined based on the three-star time difference positioning system to obtain the first set of simultaneous equations;

[0015] When the number of satellites receiving signals is greater than or equal to 4, the three independent time difference equations of the satellites are combined based on the four-satellite time difference positioning system to obtain the second set of simultaneous equations;

[0016] Initial iteration parameters are configured, and the first set of simultaneous equations or the second set of simultaneous equations is iteratively solved by the Levenberg-Marquardt method to calculate the target position.

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

[0018]

[0019] The second set of simultaneous equations is:

[0020]

[0021] Among them, the position of each satellite is (x i ,y i , z i ), i is the satellite tag, i = 1, 2, 3, 4; the position of the radiation source is (x, y, z), c is the speed of light, R e is the radius of the Earth’s sphere, τ 12 is the time difference between the unknown radiation source signal reaching satellite 1 and satellite 2, τ 13 is the time difference between the unknown radiation source signal reaching satellite 1 and satellite 3, τ 14 is the time difference between the unknown radiation source signal reaching satellite 1 and satellite 4.

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

[0023]

[0024] The linear equation for solving the target speed based on the four-star time difference positioning system is:

[0025]

[0026] Wherein:

[0027]

[0028]

[0029] (v x , v y , v z ) is the velocity of the radiation source on the XYZ three axes, and f d12 is the frequency difference of the unknown radiation source signal arriving between satellite 1 and satellite 2, and f d13 is the frequency difference of the unknown radiation source signal arriving between satellite 1 and satellite 3, and f d14 is the frequency difference of the unknown radiation source signal arriving between satellite 1 and satellite 4; the velocities of each satellite on the XYZ three axes are (v xi , v yi , v zi ), i is the satellite label, i = 1, 2, 3, 4; the position 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] Wherein, Φ 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 is the adjacent observation time interval, t k is the current moment, and t k-1 is the previous observation moment; n x[k-1] is t k-1 The system disturbance noise of the target X-axis at time t, n y [k-1] is t k-1 The system disturbance noise of the target Y-axis at time t, n z [k-1] is t k-1 The system disturbance noise of the target Z-axis at time t.

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

[0038]

[0039] Optionally, the steps of filtering optimization using the extended Kalman filter algorithm include:

[0040] Update the filtering parameters;

[0041] Establish a state prediction equation, an error covariance prediction equation and a measurement equation according to the radiation source motion state change model;

[0042] Calculate the Kalman filter gain according to the state prediction equation, the error covariance prediction equation and the measurement equation;

[0043] Obtain the state filtering value at the current time according to the Kalman filter gain;

[0044] When a new target radiation source signal is received, calculate the new target position and new target speed through the steps of time difference frequency difference value calculation, target position calculation, and target speed calculation, and update the filtering parameters at the new time based on the calculated new target position and new target speed. Optionally, the state prediction equation is:

[0045]

[0046] Wherein, is the state prediction value at time t k , is the state filtering value at time t k-1 , Φ[k,k-1] is the state transition matrix corresponding to time t k-1 and time t k ;

[0047] The error covariance prediction equation is:

[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 one-step prediction mean square error matrix corresponding to time t k-1 and time t k The one-step prediction mean square error matrix at time t k-1 The estimated mean square error matrix at time t k-1 The system disturbance matrix at time t x / y / z and Q k-1 is the system noise at time t

[0050] The measurement equation at time t k is as follows

[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 one-step prediction mean square error matrix corresponding to time t k-1 and time t k The one-step prediction mean square error matrix at time t x / y / z and R k is the observation noise at time t

[0055] The state filtering value at the current time is calculated based on the following formula

[0056]

[0057] where is the state prediction value at time t k K[k] is the Kalman filter gain at time t k and O x / y / z [k] is the observation value at time t k

[0058] ​Through the extended Kalman filter algorithm, the present invention iteratively optimizes using the observation data at the previous and current moments. The obtained result has a higher accuracy compared to the single positioning result. At the same time, high-order motion state information such as velocity and acceleration can also be obtained, and the estimation effect of the target motion state is better. Based on the calculated target position, the present invention also introduces the frequency difference information to construct a linear equation system to solve the target velocity, with fast calculation speed, and provides high-dimensional input for the extended Kalman filter algorithm, making its effect better. At the same time, when performing Kalman filtering, the present invention innovatively decouples the XYZ three dimensions for separate filtering. On the one hand, since the three-dimensional motion intensities of the target must be different, the three-dimensional decoupling can set more adaptable parameters according to the actual situation. On the other hand, after decoupling, the order of the calculation matrix required in a single dimension is reduced, and at the same time, parallel calculation methods can be used for the three dimensions to calculate synchronously, greatly improving the overall calculation efficiency. BRIEF DESCRIPTION OF THE DRAWINGS

[0059] Figure 1 It is a flowchart of the steps of the high-precision satellite passive positioning and tracking method for moving targets provided by an embodiment of the present invention;

[0060] Figure 2 It is a schematic diagram of the specific execution process of the high-precision satellite passive positioning and tracking method for moving targets provided by an embodiment of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0061] In order to make the objectives, technical solutions and advantages of the present invention clearer, the present invention will be further described in detail below with reference to the drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention, and are not used to limit the present invention.

[0062] It should be noted that the references to "one embodiment", "embodiment", "example embodiment", etc. in this specification mean that the described embodiment may include specific features, structures or characteristics, but not every embodiment must include these specific features, structures or characteristics. In addition, such expressions do not refer to the same embodiment. Further, when combining embodiments to describe specific features, structures or characteristics, whether or not there is an explicit description, it has been shown that it is within the knowledge of those skilled in the art to combine such features, structures or characteristics into other embodiments.

[0063] In addition, in the specification and the subsequent claims, certain terms are used to refer to specific components or parts. Those with ordinary knowledge in the relevant field should understand that the manufacturer can use different nouns or terms to refer to the same component or part. The specification and the subsequent claims do not use the difference in names as a way to distinguish components or parts, but use the difference in the functions of components or parts as the criterion for distinction. The terms "comprising" and "including" mentioned throughout the specification and the subsequent claims are open-ended terms and should be interpreted as "including but not limited to". In addition, the term "connected" herein includes any direct and indirect electrical connection means. Indirect electrical connection means include connection through other devices.

[0064] Before describing the embodiments of the present application in detail, first briefly describe the technical concept of the present application: The present invention solves the problems of low single-time accuracy and discontinuous motion state estimation in traditional passive positioning through joint time difference and frequency difference measurement and iterative optimization of historical data. First, use the cross ambiguity function of satellite received signals to extract time difference and frequency difference parameters; based on the time difference information, solve the target position through a three-satellite or four-satellite positioning system; combine the frequency difference parameters with satellite speed information to construct a linear equation to directly solve the target speed; finally, decouple the target three-dimensional motion state into an independent model of the XYZ axes, and use the extended Kalman filter (EKF) to perform parallel filtering optimization on the position, speed, and acceleration. By fusing historical observation data and real-time calculation results, iteratively correct the target motion parameters, while reducing the computational complexity, significantly improving the positioning accuracy and tracking stability.

[0065] Next, combine specific embodiments to describe the specific principle of the high-precision satellite passive positioning and tracking method for moving targets of the present application.

[0066] Figures 1 to 2 A high-precision satellite passive positioning and tracking method for moving targets provided by an embodiment of the present invention is shown, including the steps:

[0067] S101: Receive the target radiation source signal through at least two satellites, and use the cross ambiguity function to calculate the corresponding time difference parameter and frequency difference parameter between the target radiation source signals of different satellites.

[0068] In specific implementation, at the current moment, use two satellites to receive the signals of the target radiation source, denoted as s1(t) and s2(t), and use the cross ambiguity function to measure the time difference and frequency difference values at this time.

[0069] The calculation formula of the cross ambiguity function is:

[0070]

[0071] Among them, 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; * represents conjugate; by adjusting different τ and f values ​​for 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: According to the time difference parameters, the target position is calculated by using a three-star time difference positioning system or a four-star time difference positioning system. Step S102 of this embodiment includes: when the number of satellites receiving signals is 3, two independent time difference equations of the satellites and the earth spherical equation are combined based on the three-star time difference positioning system to obtain a first set of simultaneous equations; when the number of satellites receiving signals is greater than or equal to 4, three independent time difference equations of the satellites are combined based on the four-star time difference positioning system to obtain a second set of simultaneous equations; iterative initial parameters are configured, and the first set of simultaneous equations or the second set of simultaneous equations are iteratively solved by the Levenberg-Marquardt method to solve the target position.

[0073] In the specific implementation, this embodiment specifically adopts the Levenberg-Marquardt method to solve the target position, which is specifically implemented through the following sub-steps:

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

[0075] After obtaining the time difference values ​​between multiple satellites at the current moment, the positioning system can be selected according to the number of satellites receiving signals to preliminarily calculate the target position at this time. When the number of satellites is less than 3, the moving target cannot be positioned; when the number of satellites is 3, the three-star time difference positioning system can be used to calculate the target position; when the number of satellites is 4, the four-star time difference positioning system can be used to calculate the target position; when the number of satellites is greater than 4, four of them can be selected to use the four-star time difference positioning system to calculate the target position.

[0076] The three-star time difference positioning system or the four-star time difference positioning system used will be explained below.

[0077] The three-star time difference positioning system uses two independent time difference equations for the radiation source signal to reach three satellites to solve the spherical equations determined by the earth's surface to solve the coordinates of the radiation source, that is, the first set of simultaneous equations is:

[0078]

[0079] If the number of satellites receiving signals is three, denoted as satellite 1, satellite 2 and satellite 3, it is only necessary to combine the independent time difference equations corresponding to satellite 1 and the other two satellites (satellite 2 and satellite 3) with the spherical equations determined by the earth's surface to solve the coordinates of the radiation source.

[0080] The four-satellite time difference positioning system uses three independent time difference equations of the radiation source signal reaching the four satellites to solve the coordinates of the radiation source, that is, the second set of simultaneous equations is:

[0081]

[0082] Among them, the position of each satellite is (x i ,y i , z i ), i is the satellite tag, i = 1, 2, 3, 4; the position of the radiation source is (x, y, z), c is the speed of light, R e is the radius of the Earth’s sphere, τ 12 is the time difference between the unknown radiation source signal reaching satellite 1 and satellite 2, τ 13 is the time difference between the unknown radiation source signal reaching satellite 1 and satellite 3, τ 14 is the time difference between the unknown radiation source signal reaching satellite 1 and satellite 4.

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

[0084] Among them, the position of each satellite is (x i ,y i , z i ), i is the satellite tag, i = 1, 2, 3, 4; the position of the radiation source is (x, y, z), c is the speed of light, R e is the radius of the Earth’s sphere, τ 12 is the time difference between the unknown radiation source signal reaching satellite 1 and satellite 2, τ 13 is the time difference between the unknown radiation source signal reaching satellite 1 and satellite 3, τ 14 is the time difference between the unknown radiation source signal reaching satellite 1 and satellite 4; τ 12 , τ 13 and τ 14 The solution is obtained through the above step S101.

[0085] Since the positioning equation is a nonlinear equation, it can be solved using common nonlinear equation solving methods, so the corresponding Jacobian matrix and objective function need to be clarified.

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

[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 as follows:

[0089]

[0090] Where:

[0091]

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

[0093] Denote the iteration number as j. At the same time, select the initial iteration point (x j , y j , z j ) according to engineering experience within the satellite coverage area, j = 0, and take the parameters ε, μ0, γ1, γ2, λ1, λ2; and make

[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 Jacobian matrix J(x j , y j , z j ) and objective function F(x j , y j , z j ) corresponding to this iteration (x j , y j , z j ).

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

[0101] Specifically, determine whether the following inequality holds. If it holds, end the iteration and output the current (x j , y j , z j ) as the calculated target position. Otherwise, proceed to the next step.

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

[0103] Subsequently, for simplicity of expression, denote J 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 the vector d j :

[0106]

[0107] where the matrix I is the identity matrix.

[0108] Calculate the value ρ j :

[0109]

[0110] where, 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] Make an interval judgment on the value ρ j . When ρ j < λ1, then μ j+1 = γ2μ j ; when λ1 ≤ ρ j < λ2, then μ j+1 = μ j ; when λ2 ≤ ρ j , then μ j+1 = γ1μ j ;

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

[0113] S103: Construct a linear equation to solve for the target velocity based on the target position, the frequency difference parameter, and the satellite velocity information. After calculating the target position at the current moment, combining it with the known satellite position and velocity information and the frequency difference parameter, the target velocity can be calculated.

[0114] After obtaining the target position using the three-satellite time difference positioning system, the target velocity can be characterized in the following matrix form:

[0115]

[0116] After obtaining the target position using the four-satellite time difference positioning system, the target velocity can be characterized in the following matrix form:

[0117]

[0118] Both of the above matrix representation forms are linear, so the target velocity can be quickly obtained by matrix inversion.

[0119] Where:

[0120]

[0121]

[0122] (v x , v y , v z ) is the velocity of the radiation source on the XYZ three axes, f d12 is the frequency difference between the signals from the unknown radiation source arriving at satellite 1 and satellite 2, f d13 is the frequency difference between the signals from the unknown radiation source arriving at satellite 1 and satellite 3, f d14 is the frequency difference between the signals from the unknown radiation source arriving at satellite 1 and satellite 4; f d12 , f d13 and f d14 are solved through the above step S101. The velocities of each satellite on the XYZ three axes are (v xi , v yi , v zi ), i is the satellite label, i = 1, 2, 3, 4; the position of the radiation source is (x, y, z), and c is the speed of light.

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

[0124] After obtaining the target position and target speed information at the current moment, higher-precision target position and speed information can be obtained through the tracking algorithm. In this embodiment, based on the extended Kalman filter algorithm framework, decoupling processing is performed on the XYZ three dimensions, and filtering processing is performed separately on each dimension. On the one hand, the algorithm scale is reduced, and on the other hand, parallel computing can be used for simultaneous processing, effectively improving the computing efficiency.

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

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

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

[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] Among them, Φ 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 is the adjacent observation time interval, t k is the current moment, t k-1 is the previous observation moment; n x [k-1] is the system perturbation noise of the target X axis at time t k-1 , n y[k - 1] is t k-1 The system disturbance noise of the target Y - axis at time t, n z [k - 1] is t k-1 The system disturbance noise of the target Z - axis at time t.

[0134] s x [k] is t k The state matrix of the target X - axis at time t, s x [k - 1] is t k-1 The state matrix of the target X - axis at time t, s y [k] is t k The state matrix of the target Y - axis at time t, s y [k - 1] is t k-1 The state matrix of the target Y - axis at time t, s z [k] is t k The state matrix of the target Z - axis at time t, 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 The time instant and t k The state transition matrix between time instants, Γ[k - 1] is t k-1 The system disturbance 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 ) is the position of the target on the X - axis at time t k , v x (t k ) is the velocity of the target on the X - axis at time t k , a x (t k ) is the acceleration of the target on the X - axis at time t k , b x (t k ) is the second - order acceleration of the target on the X - axis at time t k ; y(t k ) is the position of the target on the Y - axis at time t k , v y (t k ) is the velocity of the target on the Y - axis at time t k , a x (t k ) is the acceleration of the target on the Y - axis at time t k , b y (t k ) is the acceleration of the target on the Y - axis at time t kThe second-order acceleration of the target Y-axis at a moment; z(t k ) is the position of the target Z-axis at time t k , v z (t k ) is the velocity of the target Z-axis at time t k , a z (t k ) is the acceleration of the target Z-axis at time t k , b z (t k ) is the second-order acceleration of the target Z-axis at time t k .

[0146] Since the XYZ three axes are completely dual, for simplicity of expression, the XYZ three axes will be uniformly expressed below, that is, the state matrix of the XYZ three axes at time t k is represented as s x / y / z [k], the system disturbance noise is represented as n x / y / z [k], the state filtering value is represented as the state prediction value is represented as the observation value is represented as O x / y / z [k], the system noise is represented as O x / y / z , and the observation noise is represented as R x / y / z .

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

[0148] Record the time mark number k = 1. After obtaining the target position (x, y, z) and velocity (v x , v y , v z ) for the first time, the initial state filtering value of the extended Kalman filtering algorithm is set as:[[]]

[0149]

[0150] At the same time, set the initial estimated mean square error matrix P[1], the system noise coefficients Q x / Q y / Q z of the XYZ three axes, and the observation noise parameters R x / R y / R z .

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

[0152] After receiving the signal at a new moment, after the above time difference and frequency difference value calculation, target position calculation, and target speed calculation, new target position and speed information is obtained. Using a filtering algorithm on it can obtain a position and speed value with higher accuracy. In this embodiment, the Kalman filtering algorithm framework is used for processing, and the specific implementation steps include:

[0153] 1. Update the filtering parameters.

[0154] Increment the time stamp count k by one. At this time, the new signal reception time is denoted as t k , and the previous signal reception time is denoted as t k-1 time. Thus, the time difference between t k-1 time and t k time is obtained, and the state transition matrix Φ[k,k - 1] at t k-1 time and the system disturbance matrix Γ[k - 1] at t x , v y , v z ). Update the observation values O x / y / z [k] of the XYZ three axes based on the target position (x, y, z) and speed (v

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

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

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

[0158] 2. Establish a state prediction equation, an error covariance prediction equation, and a measurement equation according to the radiation source motion state change model.

[0159] Specifically, the state prediction equation is:

[0160]

[0161] where is the state prediction value at t k time, is the state filtering value at t k-1 time, and Φ[k,k - 1] is the state transition matrix corresponding to t k-1 time and t k time;

[0162] The error covariance prediction equation is:

[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 one - step prediction mean - square error matrix corresponding to time t k-1 and time t k , P[k - 1] is the estimated mean - square error matrix at time t k-1 , Γ[k - 1] is the system disturbance matrix at time t k-1 , Q x / y / z is the system noise at time t k-1 ;

[0165] The measurement equation at time t k is as follows:

[0166]

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

[0168] The Kalman filter gain at time t k 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 one - step prediction mean - square error matrix corresponding to time t k-1 and time t k , R x / y / z is the observation noise at time t k ;

[0171] 4. Obtain the state filtering value at the current time according to the Kalman filter gain.

[0172] The state filtering value at the current time t k is calculated based on the following formula:

[0173]

[0174] where, is the state prediction value at time t k , K[k] is the Kalman filter gain at time t k and Ox / y / z [k] is the observed value at time t k .

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

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

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

[0178] The first element of the state filtering value corresponds to the position of the target on the X-axis, the second element corresponds to the velocity of the target on the X-axis, and the third element corresponds to the acceleration of the target on the X-axis; the state filtering value The first element corresponds to the position of the target on the Y-axis, the second element corresponds to the velocity of the target on the Y-axis, and the third element corresponds to the acceleration of the target on the Y-axis; the state filtering value The first element corresponds to the position of the target on the Z-axis, the second element corresponds to the velocity of the target on the Z-axis, and the third element corresponds to the acceleration of the target 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, calculate the new target position and new target velocity through the steps of time difference and frequency difference value calculation, target position calculation, and target velocity calculation (i.e., step S101, step S102, step S103), and update the filtering parameters at the new time 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 time).

[0180] In summary, through the extended Kalman filter algorithm, the present invention iteratively optimizes using the observation data at the previous and current moments. The obtained result has higher accuracy compared to the single-positioning result. At the same time, high-order motion state information such as velocity and acceleration can also be obtained, and the estimation effect of the target motion state is better. Based on the calculated target position, the present invention also introduces frequency difference information to construct a linear equation system to solve the target velocity, with fast calculation speed, and provides high-dimensional input for the extended Kalman filter algorithm, making its effect better. At the same time, when performing Kalman filtering, the present invention innovatively decouples the XYZ three dimensions for separate filtering. On the one hand, since the three-dimensional motion intensities of the target must be different, the three-dimensional decoupling can set more adaptable parameters according to the actual situation. On the other hand, after decoupling, the order of the calculation matrix required in a single dimension is reduced, and at the same time, parallel calculation methods can be used for the three dimensions to calculate synchronously, greatly improving the overall calculation efficiency.

[0181] It should be noted that the present application can be implemented in software and / or a combination of software and hardware. For example, it can be implemented using an application-specific integrated circuit (ASIC), a general-purpose computer, or any other similar hardware device. In one embodiment, the software program of the present application can be executed by a processor to implement the above steps or functions. Similarly, the software program of the present application (including related data structures) can be stored in a computer-readable recording medium, such as a RAM memory, a magnetic or optical drive, or a floppy disk and similar devices. In addition, some steps or functions of the present application can be implemented using hardware, for example, as a circuit that cooperates with the processor to execute each step or function.

[0182] The method according to the present invention can be implemented on a computer as a computer-implemented method, or in dedicated hardware, or in a combination of both. The executable code or a part thereof for the method according to the present 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-temporary program code components stored on a computer-readable medium for executing the method according to the present invention when the program product is executed on a computer.

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

Claims

1. A high-precision satellite passive positioning and tracking method for moving targets, characterized in that: Includes steps: Receive target radiation source signals through at least two satellites, and calculate corresponding time difference parameters and frequency difference parameters between the target radiation source signals of different satellites using a mutual ambiguity function; According to the time difference parameters, the target position is calculated by a three-star time difference positioning system or a four-star time difference positioning system; Constructing a linear equation to solve the target speed according to the target position, the frequency difference parameter and the satellite speed information; The three-dimensional motion state of the target is decoupled into an XYZ three-axis independent model, and based on the target position and the target speed, an extended Kalman filter algorithm is used to perform filter optimization to obtain a target positioning and tracking result.

2. The high-precision satellite passive positioning and tracking method for moving targets according to claim 1, characterized in that: The calculation formula of the mutual fuzzy function is: Among them, s1(t) and s2(t) are the target radiation source signals received by 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.

3. The high-precision satellite passive positioning and tracking method for moving targets according to claim 1, characterized in that: According to the time difference parameters, the target position is calculated by using a three-star time difference positioning system or a four-star time difference positioning system, including: When the number of satellites receiving signals is 3, the two independent time difference equations of the satellites and the earth's spherical surface equation are combined based on the three-star time difference positioning system to obtain the first set of simultaneous equations; When the number of satellites receiving signals is greater than or equal to 4, the three independent time difference equations of the satellites are combined based on the four-satellite time difference positioning system to obtain the second set of simultaneous equations; Initial iteration parameters are configured, and the first set of simultaneous equations or the second set of simultaneous equations is iteratively solved by 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 set of simultaneous equations is: The second set of simultaneous equations is: Among them, the position of each satellite is (x i ,y i , z i ), i is the satellite tag, i = 1, 2, 3, 4; the position of the radiation source is (x, y, z), c is the speed of light, R e is the radius of the Earth’s sphere, τ 12 is the time difference between the unknown radiation source signal reaching satellite 1 and satellite 2, τ 13 is the time difference between the unknown radiation source signal reaching satellite 1 and satellite 3, τ 14 is the time difference between the unknown radiation source signal reaching 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 linear equation for solving the target speed based on the three-star time difference positioning system is: The linear equation for solving the target speed based on the four-star time difference positioning system is: in: (v x , v y , v z ) is the speed of the radiation source on the XYZ axis, f d12 is the frequency difference between the unknown radiation source signal reaching satellite 1 and satellite 2, f d13 is the frequency difference between the unknown radiation source signal reaching satellite 1 and satellite 3, f d14 is the frequency difference between the unknown radiation source signal reaching satellite 1 and satellite 4; the speed of each satellite on the XYZ axis is (v xi , v yi , v zi ), i is the satellite tag, i = 1, 2, 3, 4; the position of the radiation source is (x, y, z), and c is the speed of light.

6. The high-precision satellite passive positioning and tracking method for moving targets according to claim 1, characterized in that: The XYZ three-axis independent model is a motion state change model, and the expression of the motion state change model is: s x [k]=Φ[k,k-1]s x [k-1]+Γ[k-1]n x [k-1]; s y [k]=Φ[k,k-1]s y [k-1]+Γ[k-1]n y [k-1]; s z [k]=Φ[k,k-1]s z [k-1]+Γ[k-1]n z [k-1]; Wherein, Φ is the state transfer matrix of the extended Kalman filter, Γ is the disturbance matrix of the extended Kalman filter; Δt=t k -t k-1 is the time interval between adjacent observations, t k is the current moment, t k-1 is the last observation time; n x [k-1] is t k-1 The system disturbance noise of the target X axis at the moment, n y [k-1] is t k-1 The system disturbance noise of the target Y axis at the moment, n z [k-1] is t k-1 The system disturbance noise of the target Z axis at moment .

7. The high-precision satellite passive positioning and tracking method for moving targets according to claim 6, characterized in that: The initial state filter value of the extended Kalman filter algorithm is set as:

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

9. The high-precision satellite passive positioning and tracking method for moving targets according to claim 8, characterized in that: The state prediction equation is: in, t k The predicted value of the state at time t k-1 The state filter value at time t, Φ[k,k-1] is the state filter value at time t k-1 Time and t k The state transfer matrix at time; The error covariance prediction equation is: P'[k / k-1]=Φ[k,k-1]P[k-1]Φ T [k,k-1]+Γ[k-1]Q x / y / z C T [k-1]; Among them, 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, P[k-1] is k-1 The estimated mean square error matrix at time t, Γ[k-1] is k-1 The system perturbation matrix at time , Q x / y / z t k-1 System noise at each moment; t k The measurement equation at time is:

10. The high-precision satellite passive positioning and tracking method for moving targets according to claim 9, characterized in that: The Kalman filter gain is calculated based on the following formula: K[k]=P'[k / k-1]H T [k](H[k]P'[k / k-1]H T [k]+R x / y / z [k]) -1 ; Among them, 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 t k The observation noise at the moment; The state filter value at the current moment is calculated based on the following formula: in, t k The state prediction value at time t, K[k] is k The Kalman filter gain at time t, O x / y / z [k] is t k Observed value at time.

Images