Passive positioning method based on TDOA (time difference of arrival) and FDOA (frequency difference of arrival) joint estimation

Abstract

The invention relates to a passive positioning method based on TDOA and FDOA joint estimation. The method comprises the following steps: S0, predicting the initial position of a signal source; the method comprises the following steps: S1, calculating signal time delay and Doppler frequency shift of each sensing node signal according to a signal source and the current position of each sensing node; s2, setting a sensing node signal as a reference signal, calculating a signal time delay difference and a Doppler frequency shift difference of each sensing node signal relative to the reference signal, and compensating the corresponding sensing node signal; s3, according to a signal difference between the compensated sensing node signal and a reference signal, establishing a function representation of a signal residual error of each sensing node by taking a signal source position as a parameter; and S4, constructing a target function based on the sum of squares of signal residuals of all the sensing nodes except the reference node, and solving a signal source position enabling the target function to be minimum by using a numerical optimization algorithm. According to the method, the passive positioning resolving speed can be improved on the premise of ensuring the precision.

Description

Technical Field

[0001] The present invention relates to the technical field of passive positioning, and in particular to a passive positioning method for joint estimation of TDOA and FDOA. Background Art

[0002] Based on the classification of observation information, the commonly used ones are Direction of Arrival / DOA, Time difference of Arrival / TDOA, Frequency difference of Arrival / FDOA and Received signal Strength / RSS. In different application scenarios, the difficulty of obtaining each observation information, the credibility of the obtained information and the accuracy of the final positioning solution (estimation) are all different. For example, under the positioning requirements for non-cooperative units, the timestamp of signal transmission cannot be obtained, that is, TOA information cannot be obtained. At the same time, the wavelength of the non-cooperative signal is also unknown, and it is impossible to determine whether the baseline is much larger than half the wavelength of the signal, that is, it is impossible to confirm whether the AOA estimate avoids phase ambiguity.

[0003] The most classic direct positioning method is the grid search method, which assumes that the signal source is at each grid point. When the signal source is at a certain grid point, the distance between the assumed position and the sensing node can be calculated. When the distance is calculable, the TDOA information and FDOA information between different sensing nodes can be calculated. Then, the calculated signal is matched with the actual received signal. It is believed that the higher the matching degree, the closer the assumed position is to the unknown signal source. The cost function is used to measure the degree of signal matching. Generally speaking, the larger the value of the cost function, the higher the matching degree.

[0004] A major disadvantage of the grid search method is the large amount of calculation. Not only does it need to calculate the position of each grid point, but it also needs to calculate what form the signal received by each sensing node at that point should be. The degree of matching usually requires a large number of signal-related operations. If the cost function is a convex function (the convex function in optimization theory is actually a "concave" function), then the position of the signal source is the position corresponding to the minimum value of the function. The steepest descent method and Gauss-Newton iteration are both classic nonlinear programming algorithms. The steepest descent method is the simplest and can guarantee convergence, but the solution speed is very slow and has not been widely used; Gauss-Newton iteration has improved speed, but the current solution needs to be close to the correct solution, otherwise the algorithm result may be unstable and non-convergent. Summary of the invention

[0005] The technical problem to be solved by the present invention is to provide a passive positioning method for joint estimation of TDOA and FDOA, which can improve the passive positioning solution speed while ensuring accuracy.

[0006] The technical solution adopted by the present invention to solve the technical problem is: to provide a passive positioning method for joint estimation of TDOA and FDOA, which is applied to a system including a stationary signal source and a plurality of sensing nodes, and includes the following steps:

[0007] S0 predicts the initial position of the signal source;

[0008] S1 calculates the signal delay and Doppler frequency shift of each sensing node signal according to the signal source and the current location of each sensing node;

[0009] S2 sets a sensing node signal as a reference signal, calculates the signal delay difference and Doppler frequency shift difference of each sensing node signal relative to the reference signal, and uses the calculated signal delay difference and Doppler frequency shift difference to perform signal compensation on the corresponding sensing node signal;

[0010] S3 establishes a function representation of the residual error of each sensing node signal with the signal source position as a parameter according to the signal difference between the compensated sensing node signal and the reference signal;

[0011] S4 constructs an objective function based on the sum of squares of signal residuals of all sensing nodes except the reference node, and uses a numerical optimization algorithm to solve the signal source position that minimizes the objective function.

[0012] Furthermore, the method of using a numerical optimization algorithm to solve the signal source position that minimizes the objective function includes:

[0013] S400 sets the damping factor and maximum number of iterations;

[0014] S401 calculates the value of the objective function according to the signal source and the current position of each sensing node;

[0015] S402 calculates the derivative matrix of the objective function, and solves the update amount of the signal source position according to the calculated derivative matrix, thereby updating the signal source position;

[0016] S403 calculates an updated value of the objective function using the updated signal source position;

[0017] S404 compares the updated value of the objective function with its current value, and updates the damping factor and the value of the objective function according to the comparison result;

[0018] S405 determines whether the norm of the update amount of the signal source position is less than the set threshold. If it is less than the set threshold, the updated signal source position is output. Otherwise, it is determined whether the maximum number of iterations is reached. If the maximum number of iterations is not reached, it returns to step S401. Otherwise, if the maximum number of iterations is reached, the signal source position is output.

[0019] Further, step S404 includes:

[0020] If the updated value of the objective function is less than its current value, the damping factor is adjusted using the first coefficient and the updated value of the objective function is set as the value of the objective function; otherwise, the damping factor is adjusted using the second coefficient and the current value of the objective function is set as the value of the objective function.

[0021] Furthermore, before the step of outputting the signal source position in step S405, the method further includes:

[0022] Determine whether the norm of the update amount of the signal source position is less than the set threshold, and if it is less than the set threshold, output the signal source position.

[0023] Furthermore, the derivative matrix is ​​a Jacobian matrix, and the update amount of the signal source position is calculated by the following formula:

[0024]

[0025] Among them, Δp i is the update amount of the signal source position, is the objective function The Jacobian matrix of , λ is the damping factor, v is the residual parameter, and I is the identity matrix.

[0026] Furthermore, the residual parameter v is set to the maximum value of the residual of the signal of each sensing node calculated according to the signal source and the current position of the sensing node.

[0027] Furthermore, the predicted signal source initial position is obtained by performing a round of grid search algorithm.

[0028] Furthermore, the sensing node observes the signal source along the set moving trajectory, and its signal delay is calculated by the following formula:

[0029]

[0030] Among them, τ l,k is the signal delay, p l,k is the position of the lth sensing node at the kth observation, p is the position of the signal source, and c is the signal propagation speed.

[0031] Furthermore, the Doppler frequency shift of the sensing node is calculated by the following formula:

[0032]

[0033] Among them, f c is the frequency of the signal source, v l,k is the speed of the lth sensor node at the kth observation.

[0034] Furthermore, the objective function is the mean of the sum of squares of signal residuals of all sensing nodes except the reference node.

[0035] Beneficial Effects

[0036] Due to the adoption of the above technical scheme, the present invention has the following advantages and positive effects compared with the prior art: the present invention combines multiple information types to build a model, adopts the method of joint estimation of two observation information, TDOA and FDOA, to solve the positioning requirements of multiple mobile sensing nodes and a single stationary signal source scenario without prior information of the signal source; the patent of the present invention adopts a passive positioning method, for the case where the target itself is the signal source, only the sensing node is used to receive the signal, and the target is positioned by analyzing and processing the received signal, which is not easy to be detected and tracked by the enemy, and is more concealed and safe; the present invention adopts a direct positioning method to directly locate and solve the observed signal, and directly gives the position of the target by constructing the relationship between the signal model and the position estimation, reducing the error introduced by the intermediate measurement value estimation, especially in a low signal-to-noise ratio environment, it has higher positioning robustness; the present invention adopts an improved Levenberg-Marquardt nonlinear programming method, which can effectively solve the problems of "slow positioning" of the steepest descent method and "inaccurate positioning" of the Gauss-Newton iteration method under noisy conditions and when there is a large gap between the initial value of the passive positioning of the wireless signal and the actual position. BRIEF DESCRIPTION OF THE DRAWINGS

[0037] Figure 1 is a flow chart of an embodiment of the present invention;

[0038] Figure 2 is a schematic diagram of the coarse-grained results of the grid search method according to an embodiment of the present invention;

[0039] Figure 3 is a schematic diagram of fine-grained results of a grid search method according to an embodiment of the present invention;

[0040] Figure 4 is a flow chart of a positioning algorithm according to an embodiment of the present invention;

[0041] Figure 5 4 is a graph comparing the performance of algorithms according to the embodiments of the present invention. DETAILED DESCRIPTION

[0042] The present invention will be further described below in conjunction with specific embodiments. It should be understood that these embodiments are only used to illustrate the present invention and are not intended to limit the scope of the present invention. In addition, it should be understood that after reading the content taught by the present invention, those skilled in the art can make various changes or modifications to the present invention, and these equivalent forms fall within the scope limited by the appended claims of the application equally.

[0043] The embodiments of the present invention relate to a passive positioning method for joint estimation of TDOA and FDOA, such as Figure 1 As shown, it is applied to a system including a stationary signal source and a plurality of sensing nodes, and includes the following steps:

[0044] S0 predicts the initial position of the signal source;

[0045] S1 calculates the signal delay and Doppler frequency shift of each sensing node signal according to the signal source and the current location of each sensing node;

[0046] S2 sets a sensing node signal as a reference signal, calculates the signal delay difference and Doppler frequency shift difference of each sensing node signal relative to the reference signal, and uses the calculated signal delay difference and Doppler frequency shift difference to perform signal compensation on the corresponding sensing node signal;

[0047] S3 establishes a function representation of the residual error of each sensing node signal with the signal source position as a parameter according to the signal difference between the compensated sensing node signal and the reference signal;

[0048] S4 constructs an objective function based on the sum of squares of signal residuals of all sensing nodes except the reference node, and uses a numerical optimization algorithm to solve the signal source position that minimizes the objective function.

[0049] The initial position of the signal source is directly located by performing a grid search in a large area of ​​interest. Figure 2 and Figure 3 As shown. Simply put, the grid search method assumes that the signal source is at each grid point. When at a certain grid point, the distance between the assumed position of the signal source and the sensing node is calculable. When the distance is calculable, the TDOA information and FDOA information between different sensing nodes can be calculated. Then, the calculated signal is matched with the actual received signal. It is believed that the higher the matching degree, the closer the assumed position is to the unknown signal source. The cost function is used to measure the degree of signal matching. Generally speaking, the larger the value of the cost function, the higher the matching degree.

[0050] The condition for terminating the search of the grid search method is generally that the distance between the grid points is less than a certain threshold. When the number of iterations is very large and the distance between the grid points is close enough, the positioning accuracy can be extremely high. However, its major disadvantage is the large amount of calculation. Not only the position of each grid point must be calculated, but also the form of the signal received by each sensing node at that point must be calculated. The matching degree usually requires a large number of signal-related operations.

[0051] The application scenario of this embodiment is set as a stationary signal source and L mobile sensing nodes in the region of interest and frequency band of interest. The sensing nodes can be stationary or mobile. The sensing nodes observe the signal source K times in the specified moving trajectory. The position of the signal source is denoted as p = (x, y), and the position of the lth sensing node at the kth observation is p l,k , the instantaneous speed is vx ,k , where 1≤l≤L, 1≤k≤K. The frequency point where the signal source is located is f c , then the signal received by the lth sensor node at the kth observation is

[0052]

[0053] Among them, s k (t) is the signal envelope, τ l,k is the signal delay, f l,k is the Doppler shift, w l,k (t) is noise. Signal delay and Doppler frequency shift can be calculated by the following formula:

[0054]

[0055] In order to facilitate the calculation of the signal after delay and Doppler, the signal r received by the lth sensing node at the kth observation l.k (t) is subjected to discrete Fourier transform to obtain r l.k Frequency domain representation of (t)

[0056]

[0057] in, and They are k (t) and w l,k (t) is the discrete Fourier transform.

[0058] The goal of the nonlinear programming algorithm in the iterative process is to minimize the residual function, so the definition of the residual function affects the algorithm's computational complexity and performance to a certain extent. Since TDOA and FDOA information is used, and the "difference" is defined as the difference between the lth sensing node and the first sensing node, this implementation method takes the signal of the first sensing node as the reference signal, and the "residual" is the difference between the reference signal and the remaining L-1 sensing nodes after some compensation.

[0059] In the i-th iteration process, the current signal source position is solved as The delay and Doppler shift of each node can be calculated as follows:

[0060]

[0061] Therefore, the delay difference vector between the lth sensing node and the first sensing node at the kth observation is obtained: and the Doppler frequency shift difference vector

[0062]

[0063] The signal r actually received by the lth sensing node l,k (t) Compensation for time delay and Doppler shift difference The compensated signal Or frequency domain form

[0064]

[0065] Then, choose the residual function here for

[0066]

[0067] because It is related to the signal source position p, so the residual function can also be expressed as

[0068]

[0069] The residual is defined as

[0070]

[0071] The LM algorithm minimizes the following function in each iteration

[0072]

[0073] The derivative of f can be expressed by the Jacobian matrix J with respect to p(x,y). The items in the Jacobian matrix are

[0074]

[0075] This implementation uses the Levenberg-Marquardt algorithm to find the signal source position that minimizes the objective function f(p), such as Figure 4 As shown, the specific steps are as follows:

[0076] 1) Get the initial position of LM iteration

[0077] 2) Set the damping factor λ, the iteration termination threshold ∈, and the maximum number of iterations N max , adjustable coefficients ρ1, ρ2;

[0078] 3) Let the number of iterations i = 1;

[0079] 4) Calculate the residual e k , the objective function Let v = max(e k );

[0080] 5) Calculate the Jacobian matrix In digital signal processing, It can be calculated by numerical difference method;

[0081] 6) Find the update amount of the signal source position

[0082]

[0083] Among them, ν is the residual parameter, I is a diagonal matrix with all diagonal elements equal to 1, that is, the identity matrix;

[0084] 7) Calculate the new position of the signal source

[0085] 8) Calculate the new residual and get the updated value of the objective function

[0086] 9) If λ=λ×ρ1;

[0087] 10) If λ=λ×ρ2,

[0088] 11) If ||Δp i ||<∈andi <N max , execute step 12), otherwise i=i+1, repeat steps 4) to 10);

[0089] 12) Output the final signal source position estimate

[0090] The performance comparison of this implementation method with the grid search direct positioning algorithm and the two-step positioning algorithm based on TDOA and FDOA is shown in the figure below. Figure 5 As shown. Combined Figure 5 As can be seen from the table below, the LM fast positioning method proposed in this embodiment can achieve positioning error performance that is basically the same as the grid search direct positioning method, but the calculation time is much shorter than the grid search direct positioning method.

[0091] algorithm Average positioning time (seconds) Grid search direct positioning 1.125 Two-step positioning based on TDOA and FDOA 2.828 Levenberg-Marquardt localization 0.606

Claims

1. A passive positioning method based on joint estimation of TDOA and FDOA, characterized in that: The system is applied to a system including a stationary signal source and a plurality of sensing nodes, and comprises the following steps: S0 predicts the initial position of the signal source; S1 calculates the signal delay and Doppler frequency shift of each sensing node signal according to the signal source and the current location of each sensing node; S2 sets a sensing node signal as a reference signal, calculates the signal delay difference and Doppler frequency shift difference of each sensing node signal relative to the reference signal, and uses the calculated signal delay difference and Doppler frequency shift difference to perform signal compensation on the corresponding sensing node signal; S3 establishes a function representation of the residual error of each sensing node signal with the signal source position as a parameter according to the signal difference between the compensated sensing node signal and the reference signal; S4 constructs an objective function based on the sum of squares of signal residuals of all sensing nodes except the reference node, and uses a numerical optimization algorithm to solve the signal source position that minimizes the objective function.

2. The method according to claim 1, characterized in that The method of using a numerical optimization algorithm to solve the signal source position that minimizes the objective function includes: S400 sets the damping factor and maximum number of iterations; S401 calculates the value of the objective function according to the signal source and the current position of each sensing node; S402 calculates the derivative matrix of the objective function, and solves the update amount of the signal source position according to the calculated derivative matrix, thereby updating the signal source position; S403 calculates an updated value of the objective function using the updated signal source position; S404 compares the updated value of the objective function with its current value, and updates the damping factor and the value of the objective function according to the comparison result; S405 determines whether the norm of the update amount of the signal source position is less than the set threshold. If it is less than the set threshold, the updated signal source position is output. Otherwise, it is determined whether the maximum number of iterations is reached. If the maximum number of iterations is not reached, it returns to step S401. Otherwise, if the maximum number of iterations is reached, the signal source position is output.

3. The method according to claim 2, characterized in that Step S404 includes: If the updated value of the objective function is less than its current value, the damping factor is adjusted using the first coefficient and the updated value of the objective function is set as the value of the objective function; otherwise, the damping factor is adjusted using the second coefficient and the current value of the objective function is set as the value of the objective function.

4. The method according to claim 2, characterized in that: Before the step of outputting the signal source position in step S405, the method further includes: Determine whether the norm of the update amount of the signal source position is less than the set threshold, and if it is less than the set threshold, output the signal source position.

5. The method according to claim 2, characterized in that: The derivative matrix is ​​a Jacobian matrix, and the update amount of the signal source position is calculated by the following formula: Among them, Δp i is the update amount of the signal source position, is the objective function The Jacobian matrix of , λ is the damping factor, v is the residual parameter, and I is the identity matrix.

6. The method according to claim 5, characterized in that The residual parameter v is set to the maximum value of the signal residuals of each sensing node calculated based on the signal source and the current location of the sensing node.

7. The method according to claim 1, characterized in that The predicted signal source initial position is obtained by performing a round of grid search algorithm.

8. The method according to claim 1, characterized in that The sensing node observes the signal source along the set moving trajectory, and its signal delay is calculated by the following formula: Among them, τ l,k is the signal delay, p l,k is the position of the lth sensing node at the kth observation, p is the position of the signal source, and c is the signal propagation speed.

9. The method according to claim 8, characterized in that The Doppler frequency shift of the sensing node is calculated by the following formula: Among them, f c is the frequency of the signal source, v l,k is the speed of the lth sensor node at the kth observation.

10. The method according to claim 1, characterized in that The objective function is the mean of the sum of squares of the signal residuals of all sensing nodes except the reference node.

Images