A ground radiation source TDOA positioning method based on earth ellipsoid constraint

Abstract

The present application belongs to the technical field of wireless signal positioning, and discloses a ground-based radiated source TDOA positioning method based on the constraint of the earth ellipsoid. First, the TDOA observation quantity about the ground-based radiated source is obtained by using multiple adjacent space observation platforms, and is converted into a distance difference observation quantity. Then, a scalar product matrix is constructed by using the distance difference observation quantity, and a pseudo-linear observation equation is constructed based on the same. Next, the influence of the distance difference observation error and the prior observation error of the observation platform position vector on the pseudo-linear observation equation is quantitatively analyzed by using a first-order error analysis method, and the singularity problem of the pseudo-linear observation error covariance matrix is avoided by using matrix QR decomposition, so that an optimal weighting matrix is obtained. Then, by using the earth ellipsoid constraint to which the ground-based radiated source is subjected and the algebraic characteristics of the auxiliary variable, a positioning optimization model containing double quadratic equality constraints is constructed. The present application utilizes the earth ellipsoid constraint to which the ground-based radiated source is subjected, and effectively improves the positioning precision.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of wireless signal positioning, and particularly relates to a ground radiation source TDOA positioning method based on an earth ellipsoid constraint. BACKGROUND

[0002] Wireless signal positioning technology is widely used in intelligent city, radio astronomy, seismic survey, automatic driving, emergency rescue and many other industrial technical fields, and is also an indispensable supporting technology in military fields such as battlefield environment monitoring, tactical command cooperation, cluster control and missile guidance. Whether accurate, real-time and reliable position information services can be provided has become an important symbol of the national comprehensive strength.

[0003] As known, the observation quantity for radiation source positioning involves multi-domain parameters such as space, time, frequency and energy, and TDOA (equivalent to distance difference) is a kind of observation quantity frequently used. TDOA positioning technology is to position through the TDOA of the radiation source signal collected by multiple observation platforms, and the time difference of the radiation source signal arriving at two different observation platforms can determine a hyperbolic surface (line), and the intersection of multiple hyperbolic surfaces (lines) can obtain a radiation source position vector. With the continuous development of modern communication technology and TDOA measurement technology, TDOA positioning technology has become one of the most mainstream radiation source positioning methods.

[0004] When the radiation source is located on the earth's surface (i.e., a ground radiation source), it can be positioned by TDOA using multiple near-space observation platforms. In order to improve the TDOA positioning accuracy, the constraint of the earth ellipsoid to which the position vector of the radiation source is subjected should be integrated into the positioning method. For positioning optimization models containing equality constraints, the conventional positioning method needs to perform multi-dimensional parameter joint optimization, which is prone to problems of iteration divergence and local convergence [Guo F C, Ho K C. A quadratic constraint solution method for TDOA and FDOA localization [A]. Proceedings of the IEEE International Conference on Acoustic, Speech and Signal Processing [C]. Prague, Czech: IEEE Press, May 2011: 2588-2591.][Li Q, Chen B X, Yang ML. Improved two-step constrained total least-squares TDOA localization algorithm based on the alternating direction method of multipliers [J]. IEEE Sensors Journal, 2020, 20(22): 13666-13673.]. SUMMARY

[0005] To solve the above problems, the present application uses the implicit function derivation theory to propose a ground radiation source TDOA positioning method based on the constraint of the earth ellipsoid, which only optimizes a single Lagrange multiplier explicitly, effectively improving the robustness and global convergence of the optimization algorithm, effectively avoiding multi-dimensional parameter joint optimization problems, and reducing the computational complexity. In addition, since the new method uses the constraint of the earth ellipsoid to which the ground radiation source is subjected, the positioning accuracy is effectively improved.

[0006] To achieve the above purpose, the present application adopts the following technical solutions:

[0007] The application provides a ground radiation source TDOA positioning method based on an earth ellipsoid constraint. First, TDOA observation quantities about the ground radiation source are obtained by using multiple adjacent space observation platforms, and the TDOA observation quantities are converted into range difference observation quantities. Then, a scalar product matrix is constructed by using the range difference observation quantities, and a pseudo-linear observation equation is constructed based on the scalar product matrix. Next, the influence of the range difference observation error and the prior observation error of the observation platform position vector on the pseudo-linear observation equation is quantitatively analyzed by using a first-order error analysis method, and the singularity problem of the pseudo-linear observation error covariance matrix is avoided by using matrix QR decomposition, so that an optimal weighting matrix is obtained. Then, a positioning optimization model containing double quadratic equality constraints is constructed by using the earth ellipsoid constraint to which the ground radiation source is subjected and the algebraic characteristics of the auxiliary variables. In order to avoid joint optimization of multi-dimensional parameters, a quadratic polynomial equation between two Lagrange multipliers and unknown parameters and a monomial cubic polynomial equation are constructed, and a new calculation method of only performing explicit optimization on a single Lagrange multiplier is proposed by using the implicit function derivation theory. Finally, the positioning result of the ground radiation source is obtained by using the optimal value of the Lagrange multiplier. The method specifically comprises:

[0008] Step 1: obtaining TDOA observation quantities of the ground radiation source signal arriving at the mth observation platform and arriving at the first observation platform by using M observation platforms arranged in the near space and further obtaining range difference observation quantities by using the TDOA observation quantities

[0009] Step 2: obtaining prior observation quantities of the position vectors of all observation platforms with the first observation platform as a reference

[0010] Step 3: constructing an MxM order scalar product matrix by using and

[0011] Step 4: constructing an Mx4 order observation matrix by using and calculating an Mx5 order pseudo-linear observation matrix based on and

[0012] Step 5: calculating an M 2 x3M order observation platform position matrix by using

[0013] Step 6: setting an iteration index k:=0, setting an iteration threshold value δ, and calculating an iteration initial value based on and further obtaining an iteration initial value of the position vector of the ground radiation source​​​​​

[0014] Step 7: Based on Compute 5 x 1 order vector

[0015] Step 8: Based on and Compute M x (M-1) order distance difference observation error coefficient matrix Based on and Compute M x 3(M-1) order observation platform position vector prior observation error coefficient matrix

[0016] Step 9: Perform QR decomposition on matrix to determine M x (M-1) order column orthogonal matrix and (M-1) x (M-1) order upper triangular matrix and based on and Compute (M-1) x (M-1) order pseudo-linear observation error covariance matrix

[0017] Step 10: Based on and Compute M x M order optimal weighting matrix

[0018] Step 11: Based on and Compute 3 x 3 order matrix and two 3 x 1 order vectors and and set the initial value of Lagrange multiplier λ2

[0019] Step 12: Based on and Compute quadratic polynomial coefficients about Lagrange multiplier λ1 and and compute and the first derivative about the Lagrange multiplier estimate value and

[0020] Step 13: Based on and Compute quadratic polynomial coefficient vector about ground radiation source position vector u and and compute and about ​the first derivative of and

[0021] Step 14: based on and calculate the coefficients in the cubic equation and and calculate and the first derivative of and

[0022] Step 15: based on and establish a cubic equation about the distance between the ground radiation source and the first observation platform ρ1 = ||u-s1||2, s1 is the position vector of the first observation platform, and solve the positive real root of the equation by using the power method

[0023] Step 16: calculate the first derivative of

[0024] Step 17: based on and calculate the intermediate estimate of the ground radiation source position vector and calculate the first derivative of

[0025] Step 18: perform Newton iteration update on to obtain the updated value if the update amount go to Step 19; otherwise, let and go to Step 12;

[0026] Step 19: let and if stop the calculation; otherwise, update the iteration index k: = k+1, and go to Step 7.

[0027] Further, in the step 1, the distance difference observation is calculated in the following way

[0028]

[0029] wherein denotes the distance difference observation between the mth observation platform and the first observation platform, c is the signal propagation speed, u is the ground radiation source position vector, s m ​​​is the position vector of the mth observation platform, sm is the position vector of the 1st observation platform, is the TDOA observation of the ground radiation source signal arriving at the mth observation platform and arriving at the 1st observation platform, Δr m1 denotes the distance difference observation error between the mth observation platform and the 1st observation platform.

[0030] Further, in the step 2, the prior observation of all observation platform position vectors is obtained in the following way

[0031]

[0032] wherein is the prior observation of the 1st observation platform position vector, is the prior observation of the mth observation platform position vector, 2≤m≤M, sm m is the position vector of the mth observation platform, Δs m denotes the prior observation error.

[0033] Further, in the step 3, the M×M order scalar product matrix is constructed in the following way

[0034]

[0035] wherein denotes the prior observation of the distance between the m1th observation platform and the m2th observation platform, 1≤m1;m2≤M;

[0036] Further, in the step 4, the M×4 order observation matrix is constructed in the following way

[0037]

[0038] The M×5 order pseudo-linear observation matrix is calculated in the following way

[0039]

[0040] wherein denotes the 1st column vector in the matrix ; denotes the matrix composed of the 2nd column to the 4th column in the matrix ; denotes the 5th column vector in the matrix ; 1 M denotes the M×1 order all-1 column vector.

[0041] Further, in the step 5, M is calculated in the following way 2 ×3M-order observation platform position matrix

[0042]

[0043] where the blkdiag function is used to create a block diagonal matrix.

[0044] Further, in the step 6, the iteration initial value is calculated in the following way

[0045]

[0046] Thus, the iteration initial value of the ground radiation source position vector is obtained I3 represents a 3 × 3-order unit matrix, and 03 represents a 3 × 1-order all-0 column vector.

[0047] Further, in the step 7, the 5 × 1-order vector is calculated in the following way

[0048]

[0049] where 1 M represents an M × 1-order all-1 column vector.

[0050] Further, in the step 8, the M × (M-1)-order distance difference observation error coefficient matrix is calculated in the following way

[0051]

[0052] The M × 3(M-1)-order observation platform position vector prior observation error coefficient matrix is calculated in the following way

[0053]

[0054] where I M-1 and I M respectively represent (M-1) × (M-1)-order and M × M-order unit matrices; I 3(M-1) represents a 3(M-1) × 3(M-1)-order unit matrix; I4 represents a 4 × 4-order unit matrix; 04 represents a 4 × 1-order all-0 column vector; 0 M-1 represents an (M-1) × 1-order all-0 column vector; O 3×3(M-1) represents a 3 × 3(M-1)-order all-0 matrix; 1 M represents an M × 1-order all-1 column vector; represents a distance difference observation vector; Π represents a permutation matrix, which satisfies vec function is used to convert into column vector in the order of matrix row; scalar denotes the 1st element in vector denotes the 2nd to 5th element in vector

[0055] Further, in step 9, QR decomposition is performed on matrix to obtain

[0056]

[0057] wherein denotes M×(M-1) order orthogonal matrix; denotes M×1 order vector; denotes (M-1)×(M-1) order upper triangular matrix; 0 M-1 denotes (M-1)×1 order all-0 column vector;

[0058] Then, based on and (M-1)×(M-1) order pseudo-linear observation error covariance matrix

[0059]

[0060] wherein P r denotes distance difference observation error covariance matrix; P s denotes observation platform position vector prior observation error covariance matrix.

[0061] Further, in step 10, M×M order optimal weighting matrix

[0062]

[0063] Further, in step 11, 3×3 order matrix and two 3×1 order vectors and

[0064]

[0065] Further, in step 12, quadratic polynomial coefficients about Lagrange multiplier λ1 and

[0066]

[0067] wherein​​ Where e = 0.081819790992113 represents the first eccentricity of the Earth;

[0068]

[0069] Calculated as follows and Regarding the estimates of the Lagrange multipliers first derivative as well as

[0070]

[0071] Further, in step 13, the quadratic polynomial coefficient vector with respect to the ground radiation source location vector u is calculated in the following manner. as well as

[0072]

[0073] Calculated as follows as well as about first derivative as well as

[0074]

[0075] Furthermore, in step 14, the coefficients of the cubic equation are calculated in the following manner. as well as

[0076]

[0077] Calculated as follows as well as about first derivative as well as

[0078]

[0079] Furthermore, in step 15, the established expression for the cubic equation in one variable is:

[0080]

[0081] Further, in step 16, the calculation is performed as follows: about First derivative:

[0082]

[0083] in This represents the value of ρ1 when the iteration index is k;

[0084] Further, in step 17, the intermediate estimate of the ground radiation source location vector is calculated in the following manner.

[0085]

[0086] Calculated as follows about First derivative:

[0087]

[0088] Furthermore, in step 18, the following formula is used to... Perform Newton iterations to obtain the updated value.

[0089]

[0090] In the formula r e =6378.137km represents the radius of the Earth's equator.

[0091] Compared with the prior art, the present invention has the following beneficial effects:

[0092] This invention proposes a TDOA (Terrestrial Radiation Source Location) method based on Earth ellipsoid constraints for Earth surface radiation sources. The new method only requires explicit optimization of a single Lagrange multiplier, which effectively improves the robustness and global convergence of the optimization algorithm and reduces computational complexity. In addition, the new method can effectively improve positioning accuracy by utilizing the Earth ellipsoid constraints that the ground radiation source obeys. Attached Figure Description

[0093] Figure 1 This is a flowchart of a ground radiation source TDOA positioning method based on Earth ellipsoid constraint according to an embodiment of the present invention;

[0094] Figure 2 This is a scatter plot of the ground radiation source location results and an elliptic curve of the location error (XY plane coordinates in the ECEF coordinate system) provided in the embodiments of the present invention.

[0095] Figure 3 This is a scatter plot of the ground radiation source location results and an elliptic curve of the location error (YZ plane coordinates in the ECEF coordinate system) provided in the embodiments of the present invention.

[0096] Figure 4 is a curve of the root mean square error of the ground radiation source position vector estimation provided by the embodiment of the application varying with the distance difference observation error standard deviation σ1;

[0097] Figure 5 is a curve of the root mean square error of the ground radiation source position vector estimation provided by the embodiment of the application varying with the observation platform position vector prior observation error standard deviation σ2. DETAILED DESCRIPTION

[0098] The application will be further explained in connection with the accompanying drawings and specific embodiments:

[0099] As shown in the figure, a ground radiation source TDOA positioning method based on the constraint of the earth ellipsoid comprises: Figure 1 Step 1: M observation platforms are placed in the near space, and TDOA observation quantities of a certain ground radiation source signal arriving at the mth (2≤m≤M) observation platform and arriving at the first observation platform (i.e. the main observation platform) are obtained by using them

[0100] and distance difference observation quantities are further obtained by using the TDOA observation quantities

[0101] Step 2: the prior observation quantities of the position vectors of all observation platforms are obtained with the main observation platform as a reference

[0102] Step 3: the prior observation quantities of the position vectors of the observation platforms and the distance difference observation quantities are used to construct an M×M order scalar product matrix

[0103] Step 4: the prior observation quantities of the position vectors of the observation platforms and the distance difference observation quantities are used to construct an M×4 order observation matrix and based on and the scalar product matrix an M×5 order pseudo-linear observation matrix is calculated

[0104] Step 5: the prior observation quantities of the position vectors of the observation platforms are used to calculate an M 2 ×3M order observation platform position matrix

[0105] Step 6: let the iteration index k:=0, set an iteration threshold value δ, and based on calculate an iteration initial value thereby obtaining an iteration initial value of the ground radiation source position vector ​

[0106] Step 7: Compute 5 x 1 order vector

[0107] Step 8: Compute M x (M-1) order distance difference observation error coefficient matrix

[0108] Step 9: Perform QR decomposition on matrix to determine M x (M-1) order column orthogonal matrix and (M-1) x (M-1) order upper triangular matrix

[0109] Step 10: Compute M x M order optimal weighting matrix based on

[0110] Step 11: Compute 3 x 3 order matrix and two 3 x 1 order vectors and set the initial value of Lagrange multiplier λ2

[0111] Step 12: Compute the quadratic polynomial coefficients about Lagrange multiplier λ1 based on and compute the first derivative about the Lagrange multiplier estimate

[0112] Step 13: Compute the quadratic polynomial coefficient vector about ground radiation source position vector u based on and ​​​​​​​​​​​​​​​​​​​​​​​​​the first derivative of the estimate of the Lagrange multiplier and

[0113] Step 14: based on and calculate the coefficients in the cubic equation and and and the first derivative of the estimate of the Lagrange multiplier and

[0114] Step 15: based on and establish a cubic equation about the distance between the ground radiation source and the main observation platform ρ1 = ||u-s1||2, and solve the positive real root of the equation by using the power method

[0115] Step 16: calculate the first derivative of the estimate of the Lagrange multiplier

[0116] Step 17: based on and calculate the intermediate estimate of the ground radiation source position vector and the first derivative of the estimate of the Lagrange multiplier

[0117] Step 18: perform Newton iteration update on the estimate of the Lagrange multiplier to obtain the updated value If the update amount go to Step 19; otherwise, let and go to Step 12.

[0118] Step 19: let and If stop the calculation; otherwise, update the iteration index k: = k + 1, and go to Step 7.

[0119] Further, in Step 1, M observation platforms are placed in the near space, and they are used to perform TDOA positioning on a ground radiation source. The position vector of the ground radiation source is denoted as u, and the position vector of the mth observation platform is denoted as sm m Using them, the TDOA observation quantity of the ground radiation source signal arriving at the mth (2≤m≤M) observation platform and arriving at the first observation platform (i.e., the main observation platform) can be obtained ​​​TDOA observations Multiplying by the signal propagation speed c gives the range difference observations The corresponding expression is

[0120]

[0121] where Δr m1 denotes the range difference observation error.

[0122] Further, in step 2, the prior observations of all the platform position vectors are obtained with reference to the master platform where the prior observation of the master platform position vector is equal to its true value, and the prior observations of the other platform position vectors are based on the true values with a random error superimposed, i.e.

[0123]

[0124] where Δs m denotes the prior observation error.

[0125] Further, in step 3, the prior observations of the platform position vectors and the range difference observations are used to construct an M x M order scalar product matrix The elements in the matrix are

[0126]

[0127] where denotes the prior observation of the distance between the m1th platform and the m2th platform;

[0128] Further, in step 4, the prior observations of the platform position vectors and the range difference observations are used to construct an M x 4 order observation matrix The elements in the matrix are

[0129]

[0130] Then, based on the observation matrix and the scalar product matrix an M x 5 order pseudo-linear observation matrix is calculated The corresponding calculation formula is

[0131]

[0132] where denotes the first column vector in the matrix ; and denotes the second column vector in the matrix the matrix formed by the 2nd to 4th columns in the matrix the 5th column vector in M denotes an M x 1 all-one column vector.

[0133] Further, in the step 5, the prior observation of the observation platform position vector is utilized to calculate M 2 the 3M order observation platform position matrix The corresponding calculation formula is

[0134]

[0135] The blkdiag function is used to create a block diagonal matrix.

[0136] Further, in the step 6, the iteration index k is set to 0, the iteration threshold value δ is set, and the iteration initial value is calculated The corresponding calculation formula is

[0137]

[0138] Thus, the iteration initial value of the ground radiation source position vector is obtained where I3 represents a 3 x 3 order unit matrix; 03 represents a 3 x 1 order all-zero column vector.

[0139] Further, in the step 7, the 5 x 1 order vector is calculated The corresponding calculation formula is

[0140]

[0141] Further, in the step 8, the M x (M-1) order distance difference observation error coefficient matrix and the M x 3(M-1) order observation platform position vector prior observation error coefficient matrix The corresponding calculation formula is

[0142]

[0143] In the formula, I M-1 and I M respectively represent (M-1) x (M-1) order and M x M order unit matrices; I 3(M-1) represents a 3(M-1) x 3(M-1) order unit matrix; I4 represents a 4 x 4 order unit matrix; 04 represents a 4 x 1 order all-zero column vector; 0 M-1 represents a (M-1) x 1 order all-zero column vector; O 3×3(M-1) represents a 3 x 3(M-1) order all-zero matrix; ​Π represents the distance difference observation vector; Π represents the permutation matrix, which satisfies The `vec` function is used to convert a matrix into a column vector in column order; scalar Representing vectors The first element in the vector Representing vectors The column vector formed by the 2nd to 5th elements in the array.

[0144] Furthermore, in step 9, the matrix... QR decomposition yields

[0145]

[0146] In the formula Denotes an M×(M-1) order orthogonal matrix; Represents an M×1 order vector; Let represent an (M-1)×(M-1) order upper triangular matrix. Then, based on these two matrices, calculate the (M-1)×(M-1) order pseudolinear observation error covariance matrix. The corresponding calculation formula is:

[0147]

[0148] In the formula P r P represents the distance difference observation error covariance matrix; s This represents the prior observation error covariance matrix of the observation platform's position vector.

[0149] Furthermore, in step 10, the optimal weighting matrix of order M×M is calculated. The corresponding calculation formula is:

[0150]

[0151] Furthermore, in step 11, a 3×3 matrix is ​​calculated. and two 3×1 vectors and The corresponding calculation formula is:

[0152]

[0153] Then set the initial value of the Lagrange multiplier λ2.

[0154] Furthermore, in step 12, the coefficients of the quadratic polynomial with respect to the Lagrange multiplier λ1 are calculated. as well as The corresponding calculation formula is:

[0155]

[0156] where where e = 0.081819790992113 represents the first eccentricity of the Earth;

[0157]

[0158] Then calculate and The first derivative of the Lagrange multiplier estimate value is calculated as and The corresponding calculation formula is

[0159]

[0160] Further, in step 13, the quadratic polynomial coefficient vector about the ground radiation source position vector u is calculated and The corresponding calculation formula is

[0161]

[0162] Then calculate and The first derivative of the Lagrange multiplier estimate value is calculated as and The corresponding calculation formula is

[0163]

[0164] Further, in step 14, the coefficient in the cubic equation of one variable is calculated and The corresponding calculation formula is

[0165]

[0166] Then calculate and The first derivative of the Lagrange multiplier estimate value is calculated as and The corresponding calculation formula is

[0167]

[0168] Further, in step 15, a cubic equation of one variable about the distance ρ1 = ||u-s1||2 between the ground radiation source and the subjective observation platform is established, and the corresponding equation expression is

[0169]

[0170] Then the positive real root of the equation is solved by using the power method, and is denoted as

[0171] Further, in the step 16, the value of is calculated as The first order derivative of the Lagrange multiplier estimate is calculated as

[0172]

[0173] where denotes the value of at iteration index k.

[0174] Further, in the step 17, the intermediate estimate of the ground radiation source position vector is calculated as

[0175]

[0176] Then the value of is calculated as The first order derivative of the Lagrange multiplier estimate is calculated as

[0177]

[0178] Further, in the step 18, the Newton iteration is performed on the Lagrange multiplier estimate to obtain the updated value The corresponding calculation formula is

[0179]

[0180] where r e = 6378.137 km denotes the Earth equatorial radius. If the update is larger than a predefined threshold, then the process goes to step 19; otherwise, let and goes to step 12.

[0181] Further, in the step 19, let and If then the process stops; otherwise, update the iteration index k: = k + 1, and goes to step 7.

[0182] As a specific example, consider the positioning of a ground radiation source with longitude 128.14° and latitude 32.20°, whose position vector in the ECEF coordinate system is u = [-3336.4424 48.9337 9.2] T(km), and the existing 7 adjacent space observation platforms are used to locate the radiation source by TDOA. The position coordinates of the 7 adjacent space observation platforms in the ECEF coordinate system are shown in Table 1.

[0183] Table 1 Position coordinates of 7 adjacent space observation platforms in the ECEF coordinate system (unit: km)

[0184]

[0185] The distance difference observation error covariance matrix is where I6 represents a 6x6 order unit matrix, 1 6×6 represents a 6x6 order all-1 matrix, and σ1 represents the distance difference observation error standard deviation. The observation platform position vector prior observation error covariance matrix is where σ2 represents the observation platform position vector prior observation error standard deviation.

[0186] First, the distance difference observation error standard deviation v1 is set to σ1 = 0.1 (km), and the observation platform position vector prior observation error standard deviation σ2 is set to σ2 = 0.1 (km), Figure 2 The ground radiation source positioning result scatter plot and the positioning error ellipse curve (X-Y plane coordinates in the ECEF coordinate system) are given, Figure 3 The ground radiation source positioning result scatter plot and the positioning error ellipse curve (Y-Z plane coordinates in the ECEF coordinate system) are given. From Figure 2 and Figure 3 It can be seen that the positioning result scatter plot shape of the positioning method disclosed in the patent is consistent with the positioning error ellipse shape, and a large probability corresponds to a large area ellipse and a small probability corresponds to a small area ellipse, thereby verifying the effectiveness of the new method.

[0187] Then, the observation platform position vector prior observation error standard deviation σ2 is set to σ2 = 0.1 (km), and the value of the distance difference observation error standard deviation σ1 is changed, Figure 4 The ground radiation source position vector estimation root mean square error curve with the change of the distance difference observation error standard deviation σ1 is given; and then the distance difference observation error standard deviation σ1 is set to σ1 = 0.1 (km), and the value of the observation platform position vector prior observation error standard deviation σ2 is changed, Figure 5 The ground radiation source position vector estimation root mean square error curve with the change of the observation platform position vector prior observation error standard deviation σ2 is given.

[0188] From Figure 4 and Figure 5It can be seen that: (1) the TDOA positioning method based on the constraint of the earth ellipsoid disclosed in the patent can asymptotically approach the corresponding Cramer-Rao bound for the root mean square error of the position vector estimation of the ground radiation source, thereby verifying the asymptotic statistical optimality of the new method; (2) compared with the TDOA positioning method without using the constraint of the earth ellipsoid, the positioning accuracy of the new method is significantly improved.

[0189] Finally, it should be pointed out that the new method disclosed in the patent avoids joint optimization of multi-dimensional parameters, and only explicitly optimizes a single Lagrange multiplier, effectively improving the robustness and global convergence of the optimization algorithm, and reducing the computational complexity.

[0190] The above only shows the preferred embodiments of the present application, and it should be noted that for those skilled in the art, without departing from the principles of the present application, a number of improvements and refinements can be made, and these improvements and refinements should also be considered within the scope of protection of the present application.

Claims

1. A ground-based radiated source TDOA positioning method based on the constraint of the Earth ellipsoid, characterized in that, Includes: Step 1: Obtain TDOA observations of the ground radiated source signal arriving at the mth observation platform and arriving at the first observation platform using M observation platforms placed in the near space And further obtain range difference observations using the TDOA observations Step 2: Obtain the prior observations of all the observation platform position vectors with the first observation platform as the reference Step 3: Utilize and Constructing the M x M order scalar product matrix Step 4: Utilize and construct an M x 4 order observation matrix and based on and compute an M x 5 order pseudo-linear observation matrix Step 5: Utilizing Compute M 2 x 3 M order observation platform position matrix Step 6: Let the iteration index k: = 0, set the iteration threshold value δ, and based on Calculate the iteration initial value Further, the iteration initial value of the ground radiation source position vector is obtained Step 7: Based on Computing 5 x 1 order vector Step 8: Based on and Compute M x (M-1) order distance difference observation error coefficient matrix Based on and Compute M x 3(M-1) order observation platform position vector prior observation error coefficient matrix Step 9: QR decomposition of the matrix determines an M x (M-1) order column-orthogonal matrix and an (M-1) x (M-1) order upper triangular matrix and based on and calculates an (M-1) x (M-1) order pseudo-linear observation error covariance matrix Step 10: Based on and Computing the M x M optimal weighting matrix Step 11: Based on and compute a 3x3 matrix of moments and two 3x1 vectors of moments and and set the initial value of the Lagrange multiplier λ2 Step 12: Based on and the quadratic polynomial coefficients about the Lagrange multiplier λ1 and and the calculation of and the first derivative about the Lagrange multiplier estimate value and ​ Step 13: Based on and calculating a quadratic polynomial coefficient vector and and calculating and the first derivative and ​ Step 14: Based on and Calculate the coefficients in the cubic equation and and calculate and The first derivative of and ​ Step 15: Based on and establish a cubic equation about the distance between the ground radiation source and the first observation platform, ρ1=||u-s1||2, u is the position vector of the ground radiation source, s1 is the position vector of the first observation platform, and solve the positive real root of the equation by using the power method Step 16: Calculate With respect to first derivative; Step 17: Based on and Computing an intermediate estimate of the ground radiance source position vector and computing the first derivative with respect to ; Step 18: Perform Newton iteration update on to get updated value If the update is go to Step 19; otherwise let and go to Step 12; Step 19: Let and If then stop the calculation; otherwise update the iteration index k := k + 1 and go to Step 7.

2. The method according to claim 1, wherein, In step 1, the distance difference observation is computed as follows wherein represents the distance difference observation between the mth observation platform and the 1st observation platform, c is the signal propagation speed, u is the ground radiation source position vector, s m is the position vector of the mth observation platform, s1 is the position vector of the 1st observation platform, is the TDOA observation of the ground radiation source signal arriving at the mth observation platform and arriving at the 1st observation platform, △r m1 represents the distance difference observation error between the mth observation platform and the 1st observation platform.

3. The method of claim 1, wherein the method is based on the constraint of the Earth ellipsoid. In step 2, the prior observation of all the observed platform position vectors is obtained in the following way wherein is a prior observation of the first observed platform position vector, is a prior observation of the mth observed platform position vector, 2≤m≤M, s m is a position vector of the mth observed platform, △s m denotes a prior observation error.

4. The method of claim 1, wherein the method is based on the constraint of the Earth ellipsoid. In step 3, the MxM order scalar product matrix is constructed in the following way In the formula represents the prior observation of the distance between the m1th observation platform and the m2th observation platform, 1≤m1;m2≤M; 5. The method of claim 1, wherein the method is based on the constraint of the Earth ellipsoid. In step 4, the M x 4 order observation matrix is constructed in the following manner The M x 5 order pseudo-linear observation matrix is calculated in the following way wherein denotes the first column vector in matrix denotes the matrix consisting of the second to fourth columns in matrix denotes the fifth column vector in matrix M denotes an M x 1 order all-one column vector.​​​ 6. The method of claim 1, wherein the method is based on the constraint of the Earth ellipsoid. In step 5, M is calculated as follows 2 x3M order observation platform position matrix where the blkdiag function is used to create a block diagonal matrix.

7. The method of claim 5, wherein the method is based on the constraint of the Earth ellipsoid. In step 6, the iteration initial value is calculated as follows Thus, the iterative initial value of the position vector of the ground radiation source is obtained I3 represents a 3x3 order unit matrix, and 03 represents a 3x1 order all-0 column vector.

8. The method of claim 1, wherein the method is based on the constraint of the Earth ellipsoid. In step 7, the 5x1 order vector is calculated as follows where 1 M denotes an M x 1 all-one column vector.

9. The method of claim 1, wherein the method is based on the constraint of the Earth ellipsoid. In step 8, the M x (M-1) order distance difference observation error coefficient matrix is calculated in the following manner The M x 3 (M-1) order observation platform position vector prior observation error coefficient matrix is calculated in the following manner In the formula I M-1 and I M Let I represent (M-1)×(M-1) and M×M identity matrices, respectively; 3(M-1) I represents a 3(M-1)×3(M-1) order identity matrix; I4 represents a 4×4 order identity matrix; 04 represents a 4×1 order all-zero column vector; 0 M-1 Represents an (M-1)×1 order all-zero column vector; O 3×3(M-1) Represents a 3×3(M-1) order all-zero matrix; 1 M Represents an M×1 order column vector of all 1s; Π represents the distance difference observation vector; Π represents the permutation matrix, which satisfies The `vec` function is used to convert a matrix into a column vector in column order; scalar Representing vectors The first element in the vector Representing vectors The column vector formed by the 2nd to 5th elements in the array.

10. The method of claim 5, wherein the method is based on the constraint of the Earth ellipsoid. In step 9, QR decomposition is performed on the matrix to obtain wherein denotes an M x (M-1) order orthogonal matrix; denotes an M x 1 order vector; denotes an (M-1) x (M-1) order upper triangular matrix; 0 M-1 denotes an (M-1) x 1 order all-zero column vector; Then based on and calculating an (M-1) x (M-1) order pseudo-linear observation error covariance matrix where P r denotes the range difference observation error covariance matrix; P s denotes the prior observation error covariance matrix of the observation platform position vector; In step 10, the M x M optimal weight matrix is calculated in the following manner In step 11, the 3x3 matrix is computed as follows and two 3x1 vectors and In said step 12, the quadratic polynomial coefficients with respect to the Lagrange multiplier λ1 are calculated in the following way and In the formula where e = 0.081819790992113 represents the first eccentricity of the Earth; The following is calculated and with respect to the Lagrange multiplier estimate of the first derivative and In step 13, the quadratic polynomial coefficient vector for the ground radiating source position vector u is computed in the following way and The first derivative of the function f(x) is calculated as follows and with respect to the first derivative of the function f(x) and In said step 14, the coefficients in the monic cubic equation are calculated in the following way and The first derivative of and with respect to the first derivative of and In step 15, the expression of the cubic equation is established as: In said step 16, the following is calculated With respect to the first derivative of wherein denotes the value of p1 for iteration index k; In step 17, the intermediate estimate of the ground radiance source position vector is computed in the following way The following is calculated With respect to The first derivative of In step 18, the Newton iteration update is performed according to ​ where r e = 6378.137 km represents the Earth equatorial radius.

Images