CN111551897B - TDOA (time difference of arrival) positioning method based on weighted multidimensional scaling and polynomial root finding under sensor position error - Google Patents

Abstract

The invention discloses a TDOA (time difference of arrival) positioning method based on weighted multidimensional scaling and polynomial rooting under the condition of the prior observation error of a sensor position, which comprises the steps of firstly obtaining TDOA observed quantity of a radiation source signal by utilizing a plurality of sensors, and constructing a scalar product matrix by utilizing distance difference observed quantity, thereby forming a multidimensional scaling pseudo-linear equation; then, the influence of TDOA observation errors and sensor position prior observation errors on a pseudo linear equation is quantitatively analyzed, so that an optimal weighting matrix is determined; secondly, constructing secondary equality constraint by utilizing algebraic features of the augmented unknown vector, and constructing a secondary equality constraint weighted least square optimization model by combining a multi-dimensional scale pseudo linear equation and a sensor position prior observed value; and finally, converting the optimization problem into a polynomial root solving problem by utilizing a Lagrange multiplier technology, and obtaining the joint estimation of the radiation source position vector and the sensor position vector by utilizing a Newton root solving method. The invention can improve the positioning precision of the radiation source and can also obtain more accurate sensor position information.

Description

TDOA (time difference of arrival) positioning method based on weighted multidimensional scaling and polynomial root finding under sensor position error

Technical Field

The invention belongs to the technical field of radiation source positioning, and particularly relates to a TDOA (time difference of arrival) positioning method based on weighted multidimensional scaling and polynomial root solving under the existence of a priori observation error of a sensor position.

Background

As is well known, radiation source positioning technology plays an important role in a variety of industrial and electronic information fields, such as target monitoring, navigation telemetry, seismic surveying, radio astronomy, emergency assistance, safety management, and the like. The basic process of radiation source positioning is to extract parameters (also called positioning observation) related to the position of the radiation source from electromagnetic signals, and then to use the parameters to calculate a radiation source position vector. Observations for radiation source localization involve multi-domain parameters such as space, time, frequency, energy, etc., with TDOA (which may be equivalently distance difference) being one type of observation that is applied more frequently. The TDOA positioning technique is to perform positioning by using the arrival time difference of radiation source signals collected by a plurality of sensors, wherein the arrival time difference of the radiation source signals at two different sensors can determine 1 hyperboloid (line), and the position coordinates of the radiation source can be obtained by intersecting a plurality of hyperboloids (lines). With the continuous development of modern communication technology and time difference measurement technology, TDOA positioning technology has become one of the most mainstream radiation source positioning means. According to the algebraic features of the TDOA observation equation, domestic and foreign scholars propose a plurality of TDOA positioning methods with excellent performance, wherein the TDOA positioning methods comprise an iteration method and a closed solution method, and the closed solution method does not need iteration operation and can effectively avoid the problems of local convergence, iteration divergence and the like, so that the TDOA positioning methods are widely favored by the scholars. In recent years, among analytic class positioning methods, researchers have proposed a TDOA positioning method based on weighted multidimensional scaling (Wei H W, Wan Q, Chen Z X, Ye S F.A novel weighted multidimensional scaling analysis for time-of-arrival-based mobile location [ J ]. IEEE Transactions on Signal Processing,2008,56(7): 3018-3022). this method has obtained a pseudo-linear equation about the radiation source position vector by constructing a scalar product matrix, and has given a closed solution of the radiation source position vector, which can achieve a better positioning effect. However, this method does not utilize quadratic constraints satisfied by the augmented unknown vector, and therefore its positioning accuracy is not asymptotically optimal. On the other hand, in an actual positioning scene, a priori observation error of a sensor position often exists, for example, the sensor is randomly arranged or the sensor is installed on an airborne platform, a ship-borne platform or the like. The findings in the literature (Ho K C, Lu X, Kovavisaruch L. Source localization using TDOA and FDOA measures in the presence of receiver location errors: analysis and solution [ J ]. IEEE Transactions on Signal Processing,2007,55(2):684-696.) indicate that a priori observation errors in the sensor location severely degrade the positioning accuracy and need to be taken into account in the positioning method. The weighted multi-dimensional scale TDOA positioning method proposed in the literature (Zhu Guaiui, Von Dazheng, Nee Weike, sensor position error, time difference positioning algorithm [ J ] electronics report, 2016,44(1):21-26.) considers the influence of the prior observation error of the sensor position, and improves the positioning accuracy, but the method does not utilize the quadratic equation constraint satisfied by the augmented unknown vector, so the positioning accuracy is not asymptotically optimal. Based on the current research situation, the invention discloses a TDOA (time difference of arrival) positioning method based on weighted multidimensional scaling and polynomial root solving under the condition of the prior observation error of the sensor position. Compared with the existing TDOA positioning method based on the weighted multidimensional scale, the new method not only can remarkably improve the positioning precision of the radiation source, but also can obtain more accurate sensor position information.

Disclosure of Invention

Aiming at the problem of poor positioning accuracy of the existing TDOA positioning method based on weighted multidimensional scaling, the invention provides the TDOA positioning method based on weighted multidimensional scaling and polynomial root finding under the condition of the prior observation error of the sensor position.

In order to achieve the purpose, the invention adopts the following technical scheme:

a TDOA positioning method based on weighted multidimensional scaling and polynomial rooting in the presence of prior observation errors of sensor positions comprises the following steps:

step 1: obtaining the observed quantity of TDOA of the radiation source signal reaching the M-th sensor and the 1 st sensor by using M sensors arranged in the space, and further obtaining the observed quantity of distance difference by using the observed quantity of TDOA

Step 2: apriori observations using sensor position

Sum-distance difference observed quantity

Constructing an (M +1) × (M +1) order distance matrix D;

and step 3: calculating an (M +1) × (M +1) order scalar product matrix W using the distance matrix D;

and 4, step 4: first, a priori observations are made using sensor locations

Sum-distance difference observed quantity

Constructing an (M +1) x 4-order matrix G, and then calculating an (M +1) x 5-order matrix T by using the matrix G;

and 5: let the iteration index k:equalto 0, set the iteration threshold value delta, according to

Iterative initial value g of W and T calculation⁽⁰⁾And s⁽⁰⁾；

Step 6: according to

g⁽⁰⁾And s⁽⁰⁾Sequentially calculating an (M +1) × (M-1) order matrix B_t1(g^(k),s^(k)) And B_t2(g^(k),s^(k)) And (M +1) × 3M order matrix B_s1(g^(k),s^(k)) And B_s2(g^(k),s^(k))；

And 7: according to B_t1(g^(k),s^(k)) And B_t2(g^(k),s^(k)) And B_s1(g^(k),s^(k)) And B_s2(g^(k),s^(k)) Calculating an (M +1) × (M-1) order matrix B_t(g^(k),s^(k))＝B_t1(g^(k),s^(k))+B_t2(g^(k),s^(k)) And (M +1) × 3M order matrix B_s(g^(k),s^(k))＝B_s1(g^(k),s^(k))+B_s2(g^(k),s^(k)) And to matrix B_t(g^(k),s^(k)) Performing singular value decomposition;

and 8: first according to B_s(g^(k),s^(k))、B_t(g^(k),s^(k)) And B after singular value decomposition_t(g^(k),s^(k)) Calculating a (4M-1) × (4M-1) order weighting matrix (Ω)^(k))^-1Then using a weighting matrix (omega)^(k))^-1Calculating a (3M +4) × (3M +4) order matrix Φ^(k)And (3M +4) × 1 order column vector

And step 9: for (3M +4) × (3M +4) order matrix

Carrying out eigenvalue decomposition;

step 10: firstly, the eigenvalue decomposition result in step 9 is used to calculate the (3M +4) x 1 order column vector

And

then according to

And

calculating a scalar quantity

Step 11: decomposing the result and scalar quantity by using the characteristic value in step 9

Calculating a scalar quantity

Step 12: using Newton's method to solve

Selecting real roots and eliminating false roots, wherein the roots are roots of a unitary 6-degree polynomial of coefficients;

step 13: computing an iteration result g using the root selected in step 12^(k+1)And s^(k+1)If g | | |^(k+1)-g^(k)||₂If the value is less than or equal to delta, the step 14 is carried out, otherwise, the iteration index k is updated to be k +1, and the step 6 is carried out;

step 14: will iterate the sequence g^(k)The first 3 components of the convergence values of the sequence are used as the estimate of the radiation source position vector, and the sequence of iterations s is iterated^(k)The converged value of } is used as an estimate of the sensor position vector.

Further, the step 1 comprises:

according to the radiation source position vector u and the position vector of the mth sensor

Obtaining the observed quantity of TDOA of the radiation source signal reaching the m-th sensor and reaching the 1 st sensor

M is more than or equal to 2 and less than or equal to M, and measuring TDOA

Multiplying by the signal propagation velocity to obtain the observed distance difference

The corresponding expressions are respectively

Wherein epsilon_m1Representing the range difference observation error.

Further, the step 2 comprises:

apriori observations using sensor position

Sum-distance difference observed quantity

Constructing an (M +1) × (M +1) order distance matrix D, with the corresponding calculation formula of

In the formula

Further, the step 3 comprises:

calculating an (M +1) × (M +1) order scalar product matrix W using an (M +1) × (M +1) order distance matrix D, the corresponding calculation formula being

In the formula

Wherein I_M+1Represents an identity matrix of order (M +1) × (M + 1); 1_(M+1)×(M+1)Represents an (M +1) × (M +1) order all 1 matrix.

Further, the step 4 comprises:

firstly, calculating an (M +1) multiplied by 4-order matrix G, wherein the corresponding calculation formula is

In the formula 1_(M+1)×1Represents a (M +1) × 1 order all 1-column vector; o is_1×3Representing all 0 row vectors of order 1 × 3;

then, the matrix G is used for calculating a (M +1) multiplied by 5-order matrix T, and the corresponding calculation formula is

In the formula

Further, the step 5 comprises:

setting the iteration index k:equalto 0, setting an iteration threshold value delta, and calculating an iteration initial value g⁽⁰⁾And s⁽⁰⁾The corresponding calculation formula is

Vector t in the formula₁Represents the 1 st column vector in the matrix T; momentMatrix T₂A matrix formed by the 2 nd to 5 th columns in the matrix T is represented.

Further, the step 6 comprises:

sequentially calculating an (M +1) × (M-1) order matrix B according to the following formula_t1(g^(k),s^(k)) And B_t2(g^(k),s^(k)) And (M +1) × 3M order matrix B_s1(g^(k),s^(k)) And B_s2(g^(k),s^(k))：

In the formula

Wherein O is_(M+1)×MRepresents an (M +1) x M order all 0 matrix; 1_M×1Representing an M × 1 order all-1 column vector; i is_MRepresenting an M × M order identity matrix; i is_M-1Representing an identity matrix of order (M-1) × (M-1); o is_1×MRepresenting all 0 row vectors of order 1 × M; o is_1×(M-1)Represents all 0 row vectors of 1 × (M-1) order;

is represented by a vector s^(k)3m-2, 3m-1 and 3m elements of the column vector of order 3 × 1;

O_3×(M+1)represents a 3 × (M +1) order all 0 matrix;

O_2×(M-1)represents a 2 (M-1) order all 0 matrix; o is_(M+1)×3MRepresents (M + 1). times.3MA rank all 0 matrix;

is represented by a vector α (g)^(k),s^(k)) The 3 × 1 order column vector of the 2 nd, 3 rd and 4 th elements; o is_1×3MRepresenting all 0 row vectors of 1 × 3M order; o is_M×1Representing an M × 1 order all 0 column vector; 1_M×(M+1)Represents an M × (M +1) order all 1 matrix; i is₃Representing a 3 x 3 order identity matrix.

Further, the step 7 includes:

first, calculate the (M +1) × (M-1) order matrix B_t(g^(k),s^(k))＝B_t1(g^(k),s^(k))+B_t2(g^(k),s^(k)) And (M +1) × 3M order matrix B_s(g^(k),s^(k))＝B_s1(g^(k),s^(k))+B_s2(g^(k),s^(k)) And to matrix B_t(g^(k),s^(k)) Singular value decomposition is carried out to obtain

B_t(g^(k),s^(k))＝H^(k)Σ^(k)V^(k)T

In the formula H^(k)Represents an orthogonal matrix of (M +1) × (M-1) order; v^(k)Represents an orthogonal matrix of order (M-1) × (M-1); sigma^(k)Representing an (M-1) × (M-1) order diagonal matrix having diagonal elements of matrix B_t(g^(k),s^(k)) The singular value of (a).

Further, the step 8 includes:

first, calculate the (4M-1) × (4M-1) order weighting matrix (Ω)^(k))^-1Wherein the matrix Ω^(k)Is calculated by the formula

In the formula E_tRepresenting a TDOA observation error covariance matrix; e_sRepresenting a covariance matrix of prior observation errors of the sensor locations;

then using a weighting matrix (omega)^(k))^-1Calculating a (3M +4) × (3M +4) order matrix Φ^(k)And (3M +4) × 1 order column vector

The corresponding calculation formula is

In the formula O_4×3MRepresenting a 4 × 3M-order all 0 matrix; o is_3M×4Represents a 3 mx 4 order all 0 matrix; o is_3M×(M-1)Represents a 3 Mx (M-1) order all 0 matrix; o is_(M-1)×3MRepresenting an (M-1) × 3M order all 0 matrix; i is_3MRepresenting a 3M × 3M order identity matrix; the other expressions are

Wherein O is_4×1Representing a 4 x 1 order all 0 column vector; o is_1×4Representing all 0 row vectors of 1 × 4 order; i is₄Representing a 4 x 4 order identity matrix.

Further, the step 9 includes:

for (3M +4) × (3M +4) order matrix

Is subjected to eigenvalue decomposition to obtain

In the formula O_1×3(M-1)Representing a 1 x 3(M-1) order all 0 row vector；O_3(M-1)×1Represents a 3(M-1) × 1 order all 0 column vector; o is_{3(M-1)×3(M-1)}Represents a 3(M-1) × 3(M-1) order all 0 matrix; p^(k)Is a matrix made up of eigenvectors;

wherein

The eigenvalues are represented and arranged in descending order of absolute value from large to small, only the first 4 eigenvalues are non-zero eigenvalues, and the rest are zero eigenvalues.

Further, the step 10 includes:

firstly, the eigenvalue decomposition result in step 9 is used to calculate the (3M +4) x 1 order column vector

And

the corresponding calculation formula is

Then calculates the scalar therefrom

The corresponding calculation formula is

In the formula

Representing a vector

The jth element in (a);

representing a vector

The jth element in (a).

Further, the step 11 includes:

using in step 9

The first 4 eigenvalues and scalars of

Calculating a scalar quantity

The corresponding calculation formula is

Further, the step 12 includes:

using Newton's method to solve

Being the root of a univariate 6 th order polynomial of a coefficient, the corresponding polynomial equation can be expressed as

Let the root of the polynomial equation be

Selecting a real root using the following criteria

In the formula

And

respectively indicate the utilization of the jth root

The obtained position vector of the m-th sensor and the radiation source position vector are calculated by the corresponding formula

Wherein O is_3×(3M+1)Represents a 3 × (3M +1) order all 0 matrix; o is_3×4Represents a 3 × 4 order all 0 matrix;

express identity matrix I_MThe m-th column vector of (1).

Further, the step 13 includes:

using the root selected in step 12

Calculating an iteration result g^(k+1)And s^(k+1)The corresponding calculation formula is

In the formula O_4×4Representing a 4 x 4 order all 0 matrix;

if g | | |^(k+1)-g^(k)||₂δ is not greater than δ, go to step 14, otherwise update iteration index k ═ k +1, and go to step 6.

Further, the step 14 includes:

will iterate the sequence g^(k)The first 3 components of the convergence values of the sequence are used as the estimate of the radiation source position vector, and the sequence of iterations s is iterated^(k)The converged value of } is used as an estimate of the sensor position vector.

Compared with the prior art, the invention has the following beneficial effects:

firstly, obtaining TDOA observed quantity of a radiation source signal by utilizing a plurality of sensors in a 3-dimensional space, and constructing 1 scalar product matrix by utilizing the distance difference observed quantity, thereby forming a multi-dimensional scale pseudo linear equation; then, the influence of TDOA observation errors and sensor position prior observation errors on a pseudo linear equation is quantitatively analyzed, so that an optimal weighting matrix is determined; secondly, constructing secondary equality constraint by utilizing algebraic features of the augmented unknown vector, and constructing 1 secondary equality constraint weighted least square optimization model by combining a multi-dimensional scale pseudo linear equation and a sensor position prior observed value, wherein unknown parameters in the model simultaneously comprise a radiation source position vector and a sensor position vector; and finally, converting the optimization problem into a polynomial root solving problem by utilizing a Lagrange multiplier technology, and obtaining the joint estimation of the radiation source position vector and the sensor position vector by utilizing a Newton root solving method. The method is based on a weighted multidimensional scaling principle, converts the TDOA positioning problem into a polynomial root solving problem on the basis of the quadratic equation constraint satisfied by the augmented unknown vector, and obtains the joint estimation of the radiation source position vector and the sensor position vector by using a Newton root solving method. Compared with the existing TDOA positioning method based on the weighted multidimensional scale, the new method has two advantages: the 1 st advantage is that the positioning accuracy of the radiation source can be improved by utilizing quadratic equation constraint satisfied by the augmented unknown vector; the 2 nd advantage is that the joint estimation of the radiation source position vector and the sensor position vector is realized, the influence of the prior observation error of the sensor position is inhibited, and more accurate sensor position information can be obtained.

Drawings

FIG. 1 is a basic flow diagram of a TDOA locating method based on weighted multidimensional scaling and polynomial rooting in the presence of a priori observation errors of sensor locations in accordance with an embodiment of the present invention;

FIG. 2 is a positioning result scatter plot versus positioning error elliptic curve (X-Y coordinate plane);

FIG. 3 is a positioning result scatter plot versus positioning error elliptic curve (Y-Z coordinate plane);

FIG. 4 is a graph of root mean square error of radiation source position estimates as a function of standard deviation σ_sThe variation curve of (d);

FIG. 5 is a graph of root mean square error of sensor position estimates as a function of standard deviation σ_sThe variation curve of (d);

FIG. 6 is a plot of root mean square error of radiation source position estimate as a function of standard deviation σ_tThe variation curve of (d);

FIG. 7 is a graph of root mean square error of sensor position estimates as a function of standard deviation σ_tThe change curve of (2). (ii) a

FIG. 8 is a plot of the root mean square error of the source position estimate against the standard deviation c;

FIG. 9 is a plot of root mean square error of sensor position estimates as a function of standard deviation c;

figure 10 is a plot of the root mean square error of the radiation source position estimate against the standard deviation c.

Detailed Description

The invention is further illustrated by the following examples in conjunction with the accompanying drawings:

as shown in FIG. 1, a TDOA locating method based on weighted multidimensional scaling and polynomial rooting in the presence of a priori observation errors of sensor positions comprises the following steps:

step 1: placing M sensors in space, using them to obtain TDOA observations of the radiation source signal arriving at the M-th (M is 2-M) and arriving at the 1 st sensor, and using the TDOA observations to further obtain range-difference observations

Step 2: apriori observations using sensor position

Sum-distance difference observed quantity

Constructing an (M +1) × (M +1) order distance matrix D;

and step 3: calculating an (M +1) × (M +1) order scalar product matrix W using the distance matrix D;

and 4, step 4: first, a priori observations are made using sensor locations

Sum-distance difference observed quantity

Calculating an (M +1) x 4 order matrix G, and then calculating an (M +1) x 5 order matrix T by using the matrix G;

and 5: let the iteration index k:equalto 0, set the iteration threshold value delta, according to

Iterative initial value g of W and T calculation⁽⁰⁾And s⁽⁰⁾；

Step 6: according to

g⁽⁰⁾And s⁽⁰⁾Sequentially calculating an (M +1) × (M-1) order matrix B_t1(g^(k),s^(k)) And B_t2(g^(k),s^(k)) And (M +1) × 3M order matrix B_s1(g^(k),s^(k)) And B_s2(g^(k),s^(k))；

And 7: first according to B_t1(g^(k),s^(k)) And B_t2(g^(k),s^(k)) And B_s1(g^(k),s^(k)) And B_s2(g^(k),s^(k)) Calculating an (M +1) × (M-1) order matrix B_t(g^(k),s^(k))＝B_t1(g^(k),s^(k))+B_t2(g^(k),s^(k)) And (M +1) × 3M order matrix B_s(g^(k),s^(k))＝B_s1(g^(k),s^(k))+B_s2(g^(k),s^(k)) And to matrix B_t(g^(k),s^(k)) Performing singular value decomposition;

and 8: first according to B_s(g^(k),s^(k))、B_t(g^(k),s^(k)) And B after singular value decomposition_t(g^(k),s^(k)) Calculating a (4M-1) × (4M-1) order weighting matrix (Ω)^(k))^-1Then using a weighting matrix (omega)^(k))^-1Calculating a (3M +4) × (3M +4) order matrix Φ^(k)And (3M +4) × 1 order column vector

And step 9: for (3M +4) × (3M +4) order matrix

Carrying out eigenvalue decomposition;

step 10: firstly, the eigenvalue decomposition result in step 9 is used to calculate the (3M +4) x 1 order column vector

And

then according to

And

calculating 4 scalars

Step 11: decomposing the result and scalar quantity by using the characteristic value in step 9

Calculating 7 scalars

Step 12: using Newton's method to solve

Selecting real roots and eliminating false roots, wherein the roots are roots of a unitary 6-degree polynomial of coefficients;

step 13: computing an iteration result g using the root selected in step 12^(k+1)And s^(k+1)If g | | |^(k+1)-g^(k)||₂If the value is less than or equal to delta, the step 14 is carried out, otherwise, the iteration index k is updated to be k +1, and the step 6 is carried out;

step 14: will iterate the sequence g^(k)The first 3 components of the convergence values of the sequence are used as the estimate of the radiation source position vector, and the sequence of iterations s is iterated^(k)The converged value of } is used as an estimate of the sensor position vector.

Further, in step 1, M sensors are placed in space, and are used to perform TDOA localization of the radiation source. The radiation source position vector is u, the m sensor position vector is

Wherein the content of the first and second substances,

respectively representing the coordinates of the mth sensor in the directions of an x axis, a y axis and a z axis; using them, it is possible to obtain the observed quantities of TDOA from the radiation source signal arriving at the M (2. ltoreq. M. ltoreq. M) th sensor to the 1 st sensor

Measuring TDOA

Multiplying by the signal propagation velocity to obtain the observed distance difference

The corresponding expressions are respectively

Wherein c is the signal propagation speed; epsilon_m1Representing the range difference observation error.

Further, in the step 2, the observed quantity is observed a priori by using the position of the sensor

Sum-distance difference observed quantity

Constructing an (M +1) × (M +1) order distance matrix D, with the corresponding calculation formula of

In the formula

It is worth mentioning that it is possible to show,

are obtained in advance but contain errors therein.

Further, in step 3, the (M +1) × (M +1) order product matrix W is calculated by using the (M +1) × (M +1) order distance matrix D, and the corresponding calculation formula is

In the formula

Wherein I_M+1Represents an identity matrix of order (M +1) × (M + 1); 1_(M+1)×(M+1)Represents an (M +1) × (M +1) order all 1 matrix.

Further, in step 4, firstly, a (M +1) × 4 th order matrix G (G has no specific physical meaning, and is only an intermediate matrix) is calculated, and the corresponding calculation formula is

In the formula 1_(M+1)×1Represents a (M +1) × 1 order all 1-column vector; o is_1×3Representing all 0 row vectors of order 1 × 3;

then, a matrix G is used for calculating a (M +1) multiplied by 5-order matrix T (T has no specific physical meaning and is only a middle matrix), and the corresponding calculation formula is

In the formula

Further, in step 5, let the iteration index k:equalto 0, set the iteration threshold value δ, and calculate the iteration initial value g⁽⁰⁾And s⁽⁰⁾The corresponding calculation formula is

g⁽⁰⁾＝-(T₂ ^TW^TWT₂)^-1T₂ ^TW^TWt₁,

Vector t in the formula₁Represents the 1 st column vector in the matrix T; matrix T₂To representThe matrix T is a matrix formed by the 2 nd to 5 th columns (i.e., T ═ T)₁T₂])。

Further, in the step 6, the (M +1) × (M-1) order matrix B is calculated in sequence_t1(g^(k),s^(k)) And B_t2(g^(k),s^(k)) And (M +1) × 3M order matrix B_s1(g^(k),s^(k)) And B_s2(g^(k),s^(k))(B_t1(g^(k),s^(k))、B_t2(g^(k),s^(k))、B_s1(g^(k),s^(k))、B_s2(g^(k),s^(k)) Have no specific physical meaning, only intermediate parameters):

in particular, matrix B_t1(g^(k),s^(k)) Is calculated by the formula

In the formula O_(M+1)×MRepresents an (M +1) x M order all 0 matrix; 1_M×1Representing an M × 1 order all-1 column vector; the other expressions are

Wherein I_MRepresenting an M × M order identity matrix; i is_M-1Representing an identity matrix of order (M-1) × (M-1); o is_1×MRepresenting all 0 row vectors of order 1 × M; o is_1×(M-1)Represents all 0 row vectors of 1 × (M-1) order; the other expressions are

Wherein

Is represented by a vector s^(k)3m-2, 3m-1 and 3m elements of (a) to a 3 x 1 order column vector.

In particular, matrix B_t2(g^(k),s^(k)) Is calculated by the formula

In the formula

Wherein

O_3×(M+1)Represents a 3 × (M +1) order all 0 matrix;

O_2×(M-1)representing an all 0 matrix of order 2 (M-1).

In particular, matrix B_s1(g^(k),s^(k)) Is calculated by the formula

In the formula

Wherein O is_(M+1)×3MRepresenting an (M +1) × 3M order all 0 matrix.

In particular, matrix B_s2(g^(k),s^(k)) Is calculated by the formula

In the formula

Wherein

Is represented by a vector α (g)^(k),s^(k)) The 3 × 1 order column vector of the 2 nd, 3 rd and 4 th elements; o is_1×3MRepresenting all 0 row vectors of 1 × 3M order; o is_M×1Representing an M × 1 order all 0 column vector; 1_M×(M+1)Represents an M × (M +1) order all 1 matrix; i is₃Representing a 3 x 3 order identity matrix.

Further, in step 7, first, a matrix B of (M +1) × (M-1) order is calculated_t(g^(k),s^(k))＝B_t1(g^(k),s^(k))+B_t2(g^(k),s^(k)) And (M +1) × 3M order matrix B_s(g^(k),s^(k))＝B_s1(g^(k),s^(k))+B_s2(g^(k),s^(k)) And to matrix B_t(g^(k),s^(k)) Singular value decomposition is carried out to obtain

B_t(g^(k),s^(k))＝H^(k)Σ^(k)V^(k)T

In the formula H^(k)Represents an orthogonal matrix of (M +1) × (M-1) order; v^(k)Represents an orthogonal matrix of order (M-1) × (M-1); sigma^(k)Representing an (M-1) × (M-1) order diagonal matrix having diagonal elements of matrix B_t(g^(k),s^(k)) The singular value of (a).

Further, in step 8, first, a (4M-1) × (4M-1) order weighting matrix (Ω) is calculated^(k))^-1Wherein the matrix Ω^(k)Is calculated by the formula

In the formula E_tRepresenting a TDOA observation error covariance matrix; e_sRepresenting a covariance matrix of prior observation errors of the sensor locations;

then using a weighting matrix (omega)^(k))^-1Calculating a (3M +4) × (3M +4) order matrix Φ^(k)And (3M +4) × 1 order column vector

The corresponding calculation formula is

In the formula O_4×3MRepresenting a 4 × 3M-order all 0 matrix; o is_3M×4Represents a 3 mx 4 order all 0 matrix; o is_3M×(M-1)Represents a 3 Mx (M-1) order all 0 matrix; o is_(M-1)×3MRepresenting an (M-1) × 3M order all 0 matrix; i is_3MRepresenting a 3M × 3M order identity matrix; the other expressions are

Wherein O is_4×1Representing a 4 x 1 order all 0 column vector; o is_1×4Representing all 0 row vectors of 1 × 4 order; i is₄Representing a 4 x 4 order identity matrix.

Further, in the step 9, the (3M +4) × (3M +4) order matrix is processed

Is subjected to eigenvalue decomposition to obtain

In the formula O_1×3(M-1)Represents a 1 × 3(M-1) order all 0 row vector; o is_3(M-1)×1Represents a 3(M-1) × 1 order all 0 column vector; o is_{3(M-1)×3(M-1)}Represents a 3(M-1) × 3(M-1) order all 0 matrix; p^(k)Is a matrix made up of eigenvectors;

wherein

The eigenvalues are represented and arranged in descending order of absolute value from large to small, only the first 4 eigenvalues are non-zero eigenvalues, and the rest are zero eigenvalues.

Further, in the step 10, firstly, the (3M +4) × 1 order column vector is calculated by using the eigenvalue decomposition result in the step 9

And

the corresponding calculation formula is

Then 4 scalars are calculated therefrom

The corresponding calculation formula is

In the formula

Representing a vector

The jth element in (a);

representing a vector

The jth element in (a).

Further, in the step 11, the characteristic value in the step 9 is used

And scalar quantity

Calculating 7 scalars

The corresponding calculation formula is

Further, in the step 12, the solution is solved by using the Newton method

Being the root of a univariate 6 th order polynomial of a coefficient, the corresponding polynomial equation can be expressed as

Let the root of the polynomial equation be

Selecting a real root using the following criteria

In the formula

And

respectively indicate the utilization of the jth root

The obtained position vector of the m-th sensor and the radiation source position vector are calculated by the corresponding formula

Wherein O is_3×(3M+1)Represents a 3 × (3M +1) order all 0 matrix; o is_3×4Represents a 3 × 4 order all 0 matrix;

presentation sheetBit matrix I_MThe m-th column vector of (1).

Further, in step 13, the root selected in step 12 is utilized

Calculating an iteration result g^(k+1)And s^(k+1)The corresponding calculation formula is

In the formula O_4×4Representing a 4 x 4 order all 0 matrix;

if g | | |^(k+1)-g^(k)||₂δ is not greater than δ, go to step 14, otherwise update iteration index k ═ k +1, and go to step 6.

Further, in step 14, the iteration sequence g^(k)The first 3 components of the convergence values of the sequence are used as the estimate of the radiation source position vector, and the sequence of iterations s is iterated^(k)The converged value of } is used as an estimate of the sensor position vector.

To verify the effect of the invention, the following simulation experiment was performed:

assuming that the source is located using TDOA information (i.e., range difference information) obtained from 6 sensors whose location coordinates are shown in Table 1, the range difference observation error vector obeys a mean of zero and a covariance matrix of

A gaussian distribution of (a). The position vector of the sensor can not be accurately obtained, only the prior observation value can be obtained, and the prior observation error obeys that the mean value is zero and the covariance matrix is

A gaussian distribution of (a). σ here_tAnd σ_sAre all standard deviations.

TABLE 1 sensor 3D position coordinate (unit: m)

Setting the position vector of the radiation source as u [ -3300-]^T(m) the standard deviation σ_tAnd σ_sAre respectively set to sigma_t0.5 and σ_sFig. 2 shows a positioning result scatter diagram and a positioning error elliptic curve (X-Y coordinate plane); fig. 3 shows a positioning result scatter plot versus positioning error elliptic curve (Y-Z coordinate plane).

The location method disclosed in this patent is compared below to a weighted multidimensional scaling location method that does not take into account prior observation errors of the sensor locations. First, the position vector of the radiation source is set as u [ -4300-]^T(m) the standard deviation σ_sIs set to sigma_s0.6, change in standard deviation σ_tFigure 4 gives the root mean square error of the radiation source position estimate as a function of the standard deviation sigma_tThe variation curve of (d); FIG. 5 shows the root mean square error of the sensor position estimate as a function of the standard deviation σ_tThe change curve of (2).

Then, the position vector of the radiation source is set as u [ -4300-]^T(m) the standard deviation σ_tIs set to sigma_t0.3, change in standard deviation σ_sFigure 6 shows the root mean square error of the radiation source position estimate as a function of the standard deviation sigma_sThe variation curve of (d); FIG. 7 shows the root mean square error of the sensor position estimate as a function of the standard deviation σ_sThe change curve of (2).

Finally, the standard deviation sigma is calculated_tAnd σ_sAre respectively set to sigma_t0.4 and σ_s0.8, the radiation source position vector is set as u-1800-]^T+[-250 -250 -250]^Tc (m), FIG. 8 shows the variation curve of the root mean square error of the radiation source position estimation with the parameter c; figure 9 shows the root mean square error of the sensor position estimate as a function of the parameter c.

From FIG. 4 to FIG. 9It can be seen that: (1) the positioning method disclosed by the patent has higher positioning precision than a weighted multidimensional scaling positioning method without considering the prior observation error of the sensor position, and the performance difference of the two methods is along with the standard deviation sigma_sThe larger the sensor position prior observation error is, the more obvious the advantages of the positioning method disclosed by the patent are; (2) the root mean square error of the radiation source position estimation by the positioning method disclosed by the patent can reach the Cramer-Rao bound (namely the lower theoretical bound); (3) the positioning method disclosed by the patent can improve the estimation precision of the position of the sensor (compared with the prior observation precision), and the root mean square error of the estimation of the position of the sensor can reach the Cramer-Root bound (namely the lower theoretical bound).

The positioning method disclosed in this patent is compared below to a weighted multidimensional scaling positioning method that does not utilize quadratic constraints. The simulation parameters are the same as those of fig. 8 and 9, and fig. 10 shows the variation curve of the root mean square error of the radiation source position estimation along with the parameter c.

As can be seen from fig. 4 to 10: because the positioning method disclosed by the patent utilizes quadratic equation constraint obeyed by the augmented unknown vector, the positioning accuracy of the positioning method is obviously higher than that of a weighted multidimensional scaling positioning method which does not utilize quadratic equation constraint.

The above shows only the preferred embodiments of the present invention, and it should be noted that it is obvious to those skilled in the art that various modifications and improvements can be made without departing from the principle of the present invention, and these modifications and improvements should also be considered as the protection scope of the present invention.

Claims (15)

1. A TDOA positioning method based on weighted multidimensional scaling and polynomial rooting in the presence of prior observation errors of sensor positions is characterized by comprising the following steps:

step 1: obtaining the observed quantity of TDOA of the radiation source signal reaching the M-th sensor and the 1 st sensor by using M sensors arranged in the space, and further obtaining the observed quantity of distance difference by using the observed quantity of TDOA

Step 2: apriori observations using sensor position

Sum-distance difference observed quantity

Constructing an (M +1) × (M +1) order distance matrix D;

and step 3: calculating an (M +1) × (M +1) order scalar product matrix W using the distance matrix D;

and 4, step 4: first, a priori observations are made using sensor locations

Sum-distance difference observed quantity

Constructing an (M +1) x 4-order matrix G, and then calculating an (M +1) x 5-order matrix T by using the matrix G;

and 5: let the iteration index k equal to 0, set the iteration threshold value delta, according to

Iterative initial value g of W and T calculation⁽⁰⁾And s⁽⁰⁾；

Step 6: according to

g⁽⁰⁾And s⁽⁰⁾Sequentially calculating an (M +1) × (M-1) order matrix B_t1(g^(k),s^(k)) And B_t2(g^(k),s^(k)) And (M +1) × 3M order matrix B_s1(g^(k),s^(k)) And B_s2(g^(k),s^(k))；

And 7: according to B_t1(g^(k),s^(k)) And B_t2(g^(k),s^(k)) And B_s1(g^(k),s^(k)) And B_s2(g^(k),s^(k)) Calculating an (M +1) × (M-1) order matrix B_t(g^(k),s^(k))＝B_t1(g^(k),s^(k))+B_t2(g^(k),s^(k)) And (M +1) × 3M order matrix B_s(g^(k),s^(k))＝B_s1(g^(k),s^(k))+B_s2(g^(k),s^(k)) And to matrix B_t(g^(k),s^(k)) Performing singular value decomposition;

and 8: first according to B_s(g^(k),s^(k))、B_t(g^(k),s^(k)) And B after singular value decomposition_t(g^(k),s^(k)) Calculating a (4M-1) × (4M-1) order weighting matrix (Ω)^(k))^-1Then using a weighting matrix (omega)^(k))^-1Calculating a (3M +4) × (3M +4) order matrix Φ^(k)And (3M +4) × 1 order column vector

And step 9: for (3M +4) × (3M +4) order matrix

Carrying out eigenvalue decomposition;

step 10: firstly, the eigenvalue decomposition result in step 9 is used to calculate the (3M +4) x 1 order column vector

And

then according to

And

calculating a scalar quantity

Step 11: decomposing the result and scalar quantity by using the characteristic value in step 9

Calculating a scalar quantity

Step 12: using Newton's method to solve

Selecting real roots and eliminating false roots, wherein the roots are roots of a unitary 6-degree polynomial of coefficients;

step 13: computing an iteration result g using the root selected in step 12^(k+1)And s^(k+1)If g | | |^(k+1)-g^(k)||₂If the value is less than or equal to δ, the step 14 is carried out, otherwise, the iteration index k is updated to k +1, and the step 6 is carried out;

step 14: will iterate the sequence g^(k)The first 3 components of the convergence values of the sequence are used as the estimate of the radiation source position vector, and the sequence of iterations s is iterated^(k)The converged value of } is used as an estimate of the sensor position vector.

2. A TDOA location method based on weighted multidimensional scaling and polynomial rooting in the presence of a priori observation errors of sensor locations as recited in claim 1, wherein said step 1 comprises:

according to the radiation source position vector u and the position vector of the mth sensor

Obtaining the observed quantity of TDOA of the radiation source signal reaching the m-th sensor and reaching the 1 st sensor

Measuring TDOA

Multiplying by the signal propagation velocity to obtain the observed distance difference

The corresponding expressions are respectively

Wherein epsilon_m1Representing the range difference observation error.

3. A TDOA location method based on weighted multidimensional scaling and polynomial rooting in the presence of a priori observation errors of sensor locations as recited in claim 1, wherein said step 2 comprises:

apriori observations using sensor position

Sum-distance difference observed quantity

Constructing an (M +1) × (M +1) order distance matrix D, with the corresponding calculation formula of

In the formula

4. A TDOA location method based on weighted multidimensional scaling and polynomial rooting in the presence of a priori observation errors of sensor locations as recited in claim 1, wherein said step 3 comprises:

calculating an (M +1) × (M +1) order scalar product matrix W using an (M +1) × (M +1) order distance matrix D, the corresponding calculation formula being

In the formula

Wherein I_M+1Represents an identity matrix of order (M +1) × (M + 1); 1_(M+1)×(M+1)Represents an (M +1) × (M +1) order all 1 matrix.

5. A TDOA location method based on weighted multidimensional scaling and polynomial rooting in the presence of a priori observation errors of sensor locations as recited in claim 1, wherein said step 4 comprises:

firstly, calculating an (M +1) multiplied by 4-order matrix G, wherein the corresponding calculation formula is

In the formula 1_(M+1)×1Represents a (M +1) × 1 order all 1-column vector; o is_1×3Representing all 0 row vectors of order 1 × 3;

then, the matrix G is used for calculating a (M +1) multiplied by 5-order matrix T, and the corresponding calculation formula is

In the formula

6. A TDOA location method based on weighted multidimensional scaling and polynomial rooting in the presence of a priori observation errors of sensor locations as recited in claim 1, wherein said step 5 comprises:

setting the iteration index k to 0, setting an iteration threshold value delta, and calculating an iteration initial value g⁽⁰⁾And s⁽⁰⁾The corresponding calculation formula is

Vector t in the formula₁Represents the 1 st column vector in the matrix T; matrix T₂A matrix formed by the 2 nd to 5 th columns in the matrix T is represented.

7. A TDOA location method based on weighted multidimensional scaling and polynomial rooting in the presence of a priori observation errors of sensor locations as recited in claim 4, wherein said step 6 comprises:

sequentially calculating an (M +1) × (M-1) order matrix B according to the following formula_t1(g^(k),s^(k)) And B_t2(g^(k),s^(k)) And (M +1) × 3M order matrix B_s1(g^(k),s^(k)) And B_s2(g^(k),s^(k))：

In the formula

Wherein O is_(M+1)×MRepresenting the (M +1) × M order all 0 momentsArraying; 1_M×1Representing an M × 1 order all-1 column vector; i is_MRepresenting an M × M order identity matrix; i is_M-1Representing an identity matrix of order (M-1) × (M-1); o is_1×MRepresenting all 0 row vectors of order 1 × M; o is_1×(M-1)Represents all 0 row vectors of 1 × (M-1) order;

is represented by a vector s^(k)3m-2, 3m-1 and 3m elements of the column vector of order 3 × 1;

and 1 is less than or equal to m₂≤M；O_3×(M+1)Represents a 3 × (M +1) order all 0 matrix;

O_2×(M-1)represents a 2 (M-1) order all 0 matrix; o is_(M+1)×3MRepresents an (M +1) × 3M order all 0 matrix;

is represented by a vector α (g)^(k),s^(k)) The 3 × 1 order column vector of the 2 nd, 3 rd and 4 th elements; o is_1×3MRepresenting all 0 row vectors of 1 × 3M order; o is_M×1Representing an M × 1 order all 0 column vector; 1_M×(M+1)Represents an M × (M +1) order all 1 matrix; i is₃Representing a 3 x 3 order identity matrix.

8. A TDOA location method based on weighted multidimensional scaling and polynomial rooting in the presence of a priori observation errors of sensor locations as recited in claim 7, wherein said step 7 comprises:

first, calculate the (M +1) × (M-1) order matrix B_t(g^(k),s^(k))＝B_t1(g^(k),s^(k))+B_t2(g^(k),s^(k)) And (M +1) × 3M order matrix B_s(g^(k),s^(k))＝B_s1(g^(k),s^(k))+B_s2(g^(k),s^(k)) And to matrix B_t(g^(k),s^(k)) Singular value decomposition is carried out to obtain

B_t(g^(k),s^(k))＝H^(k)Σ^(k)V^(k)T

In the formula H^(k)Represents an orthogonal matrix of (M +1) × (M-1) order; v^(k)Represents an orthogonal matrix of order (M-1) × (M-1); sigma^(k)Representing an (M-1) × (M-1) order diagonal matrix having diagonal elements of matrix B_t(g^(k),s^(k)) The singular value of (a).

9. A TDOA location method based on weighted multidimensional scaling and polynomial rooting in the presence of a priori observation errors of sensor locations as recited in claim 8, wherein said step 8 comprises:

first, calculate the (4M-1) × (4M-1) order weighting matrix (Ω)^(k))^-1Wherein the matrix Ω^(k)Is calculated by the formula

In the formula E_tRepresenting a TDOA observation error covariance matrix; e_sRepresenting a covariance matrix of prior observation errors of the sensor locations;

then using a weighting matrix (omega)^(k))^-1Calculating a (3M +4) × (3M +4) order matrix Φ^(k)And (3M +4) × 1 order column vector

The corresponding calculation formula is

In the formula O_4×3MRepresenting a 4 × 3M-order all 0 matrix; o is_3M×4Represents a 3 mx 4 order all 0 matrix; o is_3M×(M-1)Represents a 3 Mx (M-1) order all 0 matrix; o is_(M-1)×3MRepresenting an (M-1) × 3M order all 0 matrix; i is_3MRepresenting a 3M × 3M order identity matrix; the other expressions are

Wherein O is_4×1Representing a 4 x 1 order all 0 column vector; o is_1×4Representing all 0 row vectors of 1 × 4 order; i is₄Representing a 4 x 4 order identity matrix.

10. A TDOA location method based on weighted multidimensional scaling and polynomial rooting in the presence of a priori observation errors of sensor locations as recited in claim 9, wherein said step 9 comprises:

for (3M +4) × (3M +4) order matrix

Is subjected to eigenvalue decomposition to obtain

In the formula O_1×3(M-1)Represents a 1 × 3(M-1) order all 0 row vector; o is_3(M-1)×1Represents a 3(M-1) × 1 order all 0 column vector; o is_{3(M-1)×3(M-1)}Represents a 3(M-1) × 3(M-1) order all 0 matrix; p^(k)Is a matrix made up of eigenvectors;

wherein

Represents a characteristic value, andthe absolute values are arranged in descending order from large to small, only the first 4 eigenvalues are non-zero eigenvalues, and the rest are zero eigenvalues.

11. A TDOA location method based on weighted multidimensional scaling and polynomial rooting in the presence of a priori observation errors of sensor locations as recited in claim 10, wherein said step 10 comprises:

firstly, the eigenvalue decomposition result in step 9 is used to calculate the (3M +4) x 1 order column vector

And

the corresponding calculation formula is

Then calculates the scalar therefrom

The corresponding calculation formula is

In the formula

Representing a vector

The jth element in (a);

representing a vector

The jth element in (a).

12. A TDOA location method based on weighted multidimensional scaling and polynomial rooting in the presence of a priori observation errors of sensor locations as recited in claim 11, wherein said step 11 comprises:

using in step 9

The first 4 eigenvalues and scalars of

Calculating a scalar quantity

The corresponding calculation formula is

13. A TDOA location method based on weighted multidimensional scaling and polynomial rooting in the presence of a priori observation errors of sensor locations as recited in claim 1, wherein said step 12 comprises:

using Newton's method to solve

Being the root of a univariate 6 th order polynomial of a coefficient, the corresponding polynomial equation can be expressed as

Let the root of the polynomial equation be

Selecting a real root using the following criteria

In the formula

And

respectively indicate the utilization of the jth root

The obtained position vector of the m-th sensor and the radiation source position vector are calculated by the corresponding formula

Wherein O is_3×(3M+1)Represents a 3 × (3M +1) order all 0 matrix; o is_3×4Represents a 3 × 4 order all 0 matrix;

express identity matrix I_MThe m-th column vector of (1).

14. A TDOA location method based on weighted multidimensional scaling and polynomial rooting in the presence of a priori observation errors of sensor locations as recited in claim 1, wherein said step 13 comprises:

using step 12Root of (2)

Calculating an iteration result g^(k+1)And s^(k+1)The corresponding calculation formula is

In the formula O_4×4Representing a 4 x 4 order all 0 matrix;

if g | | |^(k+1)-g^(k)||₂δ is not greater, go to step 14, otherwise update iteration index k ═ k +1, and go to step 6.

15. A TDOA location method based on weighted multidimensional scaling and polynomial rooting in the presence of a priori observation errors of sensor locations as recited in claim 1, wherein said step 14 comprises:

will iterate the sequence g^(k)The first 3 components of the convergence values of the sequence are used as the estimate of the radiation source position vector, and the sequence of iterations s is iterated^(k)The converged value of } is used as an estimate of the sensor position vector.

Images