A ground radiation source multi-satellite TDOA / FDOA joint 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 radiating source multi-star TDOA / FDOA joint positioning method based on the constraint of the earth ellipsoid. In order to avoid multi-dimensional parameter joint optimization, two binary polynomial equations about the distance and distance rate of change between the ground-based radiating source and the main star are constructed, and a monomial seventh-degree polynomial equation about the distance between the ground-based radiating source and the main star is constructed based on sequential linear convolution operation. A Newton iteration method for explicit optimization of only a single Lagrange multiplier is proposed based on polynomial root operation and implicit function derivation theory. Finally, the optimal value of the Lagrange multiplier and the optimal value of the two auxiliary variables are used to obtain the position parameters of the ground-based radiating source. The simulation experiment results show that, compared with the existing multi-star TDOA / FDOA joint positioning method, the present application can greatly improve the positioning accuracy.

Description

Technical Field

[0001] This invention relates to the field of wireless signal positioning technology, and in particular to a multi-satellite TDOA / FDOA joint positioning method for ground-based radiation sources based on Earth ellipsoid constraints. Background Technology

[0002] Wireless signal positioning technology is widely used in many industrial and technological fields such as smart cities, radio astronomy, seismic surveying, autonomous driving, and emergency search and rescue. It is also an indispensable supporting technology in military fields such as battlefield environment monitoring, tactical command and coordination, swarm control, and missile guidance. Undoubtedly, the ability to provide accurate, real-time, and reliable location information services has become an important indicator of a nation's comprehensive strength and a crucial strategic support for national political, economic, and defense security initiatives.

[0003] Location of ground-based radiation sources based on communication satellite platforms is a common positioning system, including multi-satellite TDOA positioning and multi-satellite TDOA / FDOA joint positioning. Compared with multi-satellite TDOA positioning, multi-satellite TDOA / FDOA joint positioning utilizes frequency domain observations, thus achieving higher positioning accuracy. However, the multi-satellite TDOA / FDOA joint positioning problem is essentially a multi-dimensional parameter joint optimization problem, which is prone to problems such as local convergence and iterative divergence [Guo FC, Ho K C. Aquadratic 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.]. Summary of the Invention

[0004] To address the aforementioned issues, this invention proposes a multi-satellite TDOA / FDOA joint positioning method for Earth surface radiation sources based on Earth ellipsoidal constraints. Utilizing polynomial root-finding and implicit function differentiation theory, a Newton-style iterative method is proposed that explicitly optimizes only a single Lagrange multiplier. This effectively integrates Earth elliptic constraints, transforming the multidimensional parameter optimization problem into an optimization problem for a single Lagrange multiplier, thus improving robustness and global convergence in multidimensional parameter optimization and reducing computational complexity. Furthermore, the new method utilizes weighted multidimensional scaling analysis theory to construct pseudo-linear observation equations, thereby improving robustness to large observation errors. Finally, the combination of Earth ellipsoidal constraints significantly enhances the accuracy of multi-satellite TDOA / FDOA joint positioning.

[0005] To achieve the above objectives, the present invention adopts the following technical solution:

[0006] This invention proposes a multi-satellite TDOA / FDOA joint positioning method for ground-based radiation sources based on Earth ellipsoidal constraints. First, TDOA / FDOA observations of the ground-based radiation source are obtained using one primary satellite and multiple secondary satellites. These observations are then converted into range difference observations and range difference rate of change observations using known signal propagation speeds and transmission frequencies, respectively. Next, two scalar product matrices and two pseudo-inverse matrices are constructed using these two types of observations, and a pseudo-linear observation matrix is ​​built by introducing two auxiliary variables. Then, a first-order error analysis method is used to calculate four first-order perturbation matrices and an observation error coefficient matrix to quantitatively describe the impact of range difference observation errors and range difference rate of change observation errors on the pseudo-linear observation equations. Singular value decomposition of the matrix is ​​used to avoid the singularity problem of the pseudo-linear observation error covariance matrix, thereby calculating the optimal weighting matrix. Finally, using the Earth ellipsoidal constraints obeyed by the ground-based radiation source and the algebraic properties of the two auxiliary variables, a positioning optimization model with triple quadratic equality constraints is constructed. To avoid joint optimization of multidimensional parameters, two bivariate polynomial equations are constructed simultaneously regarding the distance between the ground-based radiation source and the host star, and the rate of change of that distance. A univariate seventh-degree polynomial equation regarding the distance between the ground-based radiation source and the host star is then constructed based on sequential linear convolution. A Newton-style iterative method is proposed, based on polynomial root-finding and implicit function differentiation theory, to explicitly optimize only a single Lagrange multiplier. Finally, the optimal value of this Lagrange multiplier and the optimal values ​​of the two auxiliary variables are used to obtain the position parameters of the ground-based radiation source. This method specifically includes:

[0007] Step 1: Select M communication satellites to locate the ground radiation source, and use these M communication satellites to obtain the TDOA observations of the ground radiation source signal reaching the m-th satellite and reaching the 1-th satellite. and FDOA observations Furthermore, distance difference observations were obtained using TDOA and FDOA observations. Observations on the rate of change of distance difference

[0008] Step 2: Utilize Construct an M×M scalar product matrix use and Construct an M×M scalar product matrix

[0009] Step 3: Utilize Construct an M×4 order observation matrix use Construct an M×4 order observation matrix

[0010] Step 4: Utilize Construct an M×5 pseudo-inverse matrix use Construct an M×5 pseudo-inverse matrix

[0011] Step 5: Utilize and Construct a 2M×6 pseudo-linear observation matrix

[0012] Step 6: Based on and One or more parameters are used to calculate four first-order perturbation matrices in sequence. as well as

[0013] Step 7: Let the iteration index k:=0, set the iteration threshold value δ, and based on... Calculate the initial value of the iteration This yields the initial values ​​for the iterative iteration of the ground radiation source location vector.

[0014] Step 8: Based on as well as Calculate the 2M×2(M-1) order observation error coefficient matrix

[0015] Step 9: For the matrix Perform singular value decomposition to determine a 2M×2(M-1) order column orthogonal matrix. and utilize Calculate the 2(M-1)×2(M-1) order pseudolinear observation error covariance matrix.

[0016] Step 10: Based on and Calculate the optimal weighting matrix of order 2M×2M.

[0017] Step 11: Based on and Calculate the 3×3 matrix sequentially 3 3×1 vectors as well as 6 scalars as well as And set the initial value of the Lagrange multiplier λ3.

[0018] Step 12: Based on and One or more parameters in the vector are used to compute the Lagrange multiplier vector. Bivariate polynomial coefficients as well as And calculate as well as Regarding the estimates of the Lagrange multipliers first derivative as well as

[0019] Step 13: Based on and One or more parameters are used to calculate the bivariate polynomial coefficient vector with respect to the ground radiation source location vector u. as well as And calculate as well as about first derivative as well as

[0020] Step 14: Construct two simultaneous parameters related to the distance β between the ground radiation source and the first satellite, and the rate of change of that distance. The bivariate polynomial equation, based on

[0021] and The coefficients of the above bivariate polynomial equation are calculated using one or more parameters. and And calculate and about first derivative and

[0022] Step 15: Construct a univariate seventh-degree polynomial equation for the distance β between the ground radiation source and the first satellite, based on... and The coefficients of the univariate seventh-degree polynomial equation are calculated using the linear convolution operation of sequences. The positive real roots of the univariate seventh-degree polynomial equation are then found using the power method.

[0023] Step 16: Based on and Iterative values ​​for calculating the rate of change β of the distance between the ground radiation source and the first satellite.

[0024] Step 17: Calculate the vector about first derivative

[0025] Step 18: Based on and Iterative values ​​for calculating the location vector of ground radiation sources And calculate about first derivative

[0026] Step 19: [Regarding...] Perform Newton iterations to obtain the updated value. If update volume Then proceed to step 20; otherwise, let Proceed to step 12;

[0027] Step 20: Let as well as like If the calculation stops, then stop; otherwise, update the iteration index k:=k+1 and go to step 8.

[0028] Further, in step 1, the distance difference observation is calculated in the following manner.

[0029]

[0030] in This represents the observable difference in distance between the m-th satellite and the satellite reaching the 1-th satellite. This represents the TDOA observations of the ground radiation source signal reaching the m-th satellite and reaching the 1-th satellite, where c is the signal propagation speed, u is the ground radiation source position vector, and s... m Let s1 be the position vector of the m-th satellite, 2≤m≤M, and s1 be the position vector of the 1-th satellite. Δr m1 This represents the observation error of the distance difference between the m-th satellite and the 1st satellite;

[0031] The observation of the rate of change of distance difference is calculated as follows:

[0032]

[0033] In the formula This represents the observed rate of change of the distance difference between the m-th satellite and the 1st satellite. This represents the FDOA observations of the ground-based radiation source signal reaching the m-th satellite and reaching the 1-th satellite, where f0 represents the signal transmission frequency. Let m be the velocity vector of the m-th satellite. Let this be the velocity vector of the first satellite. This represents the observation error of the rate of change of the distance difference between the m-th satellite and the 1st satellite.

[0034] Furthermore, in step 2, the structure is constructed as follows: and

[0035]

[0036] In the formula Let m represent the distance between the m1-th satellite and the m2-th satellite, and This represents the inner product of the position vector difference and the velocity vector difference between the m1-th satellite and the m2-th satellite;

[0037] Furthermore, in step 3, the structure is constructed as follows: and

[0038]

[0039] Furthermore, in step 4, the structure is constructed as follows: and

[0040]

[0041] In the formula 1 M Represents an M×1 column vector of all 1s; 0 M This represents an M×1 column vector consisting entirely of zeros.

[0042] Furthermore, in step 5, the structure is constructed as follows:

[0043]

[0044] In the formula This represents the 5th column vector in the 5×5 identity matrix I5; Representation matrix The first column vector in; Represented by matrix The matrix formed by columns 2 to 4 in the matrix; Representation matrix The 5th column vector in; Representation matrix The 6th column vector in.

[0045] Further, in step 6, the calculation is performed as follows: as well as

[0046]

[0047]

[0048] In the formula I M-1 and I M These represent (M-1)×(M-1) and M×M identity matrices, respectively; the diag function is used to create diagonal matrices; 0 M-1 Represents an (M-1)×1 order all-zero column vector; O (M-1)×(M-1) and O (3M+1)×(M-1) Let represent (M-1)×(M-1) and (3M+1)×(M-1) matrices, respectively, all zeros.

[0049] and Let Π represent the distance difference observation vector and the distance difference change rate observation vector, respectively; Π represents the permutation matrix, which satisfies the following condition: The vec function is used to convert a matrix into a column vector in column order.

[0050] I4 represents a 4×4 identity matrix, and 04 represents a 4×1 column vector of all zeros.

[0051] Furthermore, in step 7, the initial values ​​for the iteration are calculated as follows:

[0052]

[0053] in express The column vector consisting of the first 4 elements, express The last element;

[0054] This leads to the initial values ​​of the ground radiation source location vector for iteration. I3 represents a 3×3 identity matrix; 03 represents a 3×1 column vector of all zeros.

[0055] Further, in step 8, the calculation is performed as follows:

[0056]

[0057] In the formula and O M×5M Representing M×M 2 The order of M×5M is an all-zero matrix.

[0058] Furthermore, in step 9, the matrix... Singular value decomposition yields

[0059]

[0060] In the formula Represents a 2M×2M left singular matrix; Denotes a 2(M-1)×2(M-1) order right singular matrix; Represents a 2(M-1)×2(M-1) order singular value diagonal matrix; O 2×2(M-1) Represents a 2×2(M-1) order all-zero matrix;

[0061] Then use Calculate the 2(M-1)×2(M-1) order pseudolinear observation error covariance matrix.

[0062]

[0063] In the formula, E represents the 2(M-1)×2(M-1) order observation error covariance matrix.

[0064] Furthermore, in step 10, the optimal weighting matrix of order 2M×2M is calculated as follows:

[0065]

[0066] Further, in step 11, the calculations are performed sequentially as follows: as well as

[0067]

[0068] Further, in step 12, the calculation is performed as follows: as well as

[0069]

[0070]

[0071] In the formula e = 0.081819790992113 represents the first eccentricity of the Earth;

[0072]

[0073] Then calculate as follows:

[0074] In the formula

[0075]

[0076] Further, in step 13, the calculation is performed as follows: as well as

[0077]

[0078]

[0079] Then calculate as follows: as well as

[0080] Furthermore, in step 14, two parameters are constructed simultaneously regarding the distance β between the ground radiation source and the first satellite, and the rate of change of that distance. Bivariate polynomial equations:

[0081]

[0082] The coefficients of the above bivariate polynomial equation are calculated as follows:

[0083]

[0084]

[0085] In the formula, the symbols <·>1 and <·>2 represent the first and second elements in the vector, respectively;

[0086] Then calculate as follows: and

[0087]

[0088] Furthermore, in step 15, a univariate seventh-degree polynomial equation concerning the distance β between the ground radiation source and the first satellite is constructed as follows:

[0089]

[0090] In the formula Obtained through linear convolution of sequences:

[0091]

[0092] In the formula, * represents the linear convolution operation of sequences.

[0093] Further, in step 16, the calculation is performed as follows:

[0094] Further, in step 17, the calculation is performed as follows:

[0095]

[0096] In the formula

[0097]

[0098] Further, in step 18, the calculation is performed as follows:

[0099] Then calculate as follows:

[0100]

[0101] Furthermore, in step 19, the following methods are used to... Perform Newton iteration updates:

[0102]

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

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

[0105] This invention proposes a multi-satellite TDOA / FDOA joint positioning method for ground-based radiation sources based on Earth ellipsoid constraints. This method utilizes polynomial root-finding operations and implicit function differentiation theory to transform the multidimensional parameter optimization problem into an optimization problem for a single Lagrange multiplier, effectively improving the robustness and global convergence of the optimization algorithm and reducing computational complexity. Furthermore, this method, combined with Earth ellipsoid constraints and weighted multidimensional scaling analysis theory, significantly improves the accuracy of multi-satellite TDOA / FDOA joint positioning. Attached Figure Description

[0106] Figure 1 This is a flowchart of a multi-satellite TDOA / FDOA joint positioning method for ground-based radiation sources based on Earth ellipsoid constraints, according to an embodiment of the present invention.

[0107] 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.

[0108] 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.

[0109] Figure 4 This is a curve showing the variation of the root mean square error of the ground radiation source location vector estimation as a function of the observation error standard deviation σ, provided in this embodiment of the invention.

[0110] Figure 5 This is a curve showing the variation of the root mean square error of the ground radiation source location vector estimation with the longitude of the ground radiation source, as provided in the embodiments of the present invention. Detailed Implementation

[0111] The present invention will be further explained below with reference to the accompanying drawings and specific embodiments:

[0112] like Figure 1 As shown, a multi-satellite TDOA / FDOA joint localization method for ground-based radiation sources based on Earth ellipsoid constraints includes:

[0113] Step 1: Select M communication satellites to locate the ground radiation source, and use them to obtain the TDOA observations of the signal from the ground radiation source reaching the m-th (2≤m≤M) satellite (called the auxiliary satellite) and the 1st satellite (called the primary satellite). and FDOA observations Furthermore, distance difference observations were obtained using TDOA and FDOA observations. Observations on the rate of change of distance difference

[0114] Step 2: Utilize distance difference observation Observations on the rate of change of distance difference Construct an M×M scalar product matrix and M×M scalar product matrix

[0115] Step 3: Utilize distance difference observation Observations on the rate of change of distance difference Construct an M×4 order observation matrix and M×4 order observation matrix

[0116] Step 4: Utilize the observation matrix and observation matrix Construct an M×5 pseudo-inverse matrix and M×5 pseudo-inverse matrix

[0117] Step 5: Use the scalar product matrix scalar product matrix and pseudo-inverse matrix and pseudo-inverse matrix Construct a 2M×6 pseudo-linear observation matrix

[0118] Step 6: Calculate the four first-order perturbation matrices sequentially. as well as

[0119] Step 7: Let the iteration index k:=0, set the iteration threshold value δ, and based on... Calculate the initial value of the iteration This yields the initial values ​​for the iterative iteration of the ground radiation source location vector.

[0120] Step 8: Calculate the 2M×2(M-1) order observation error coefficient matrix

[0121] Step 9: For the matrix Perform singular value decomposition to determine a 2M×2(M-1) order column orthogonal matrix. and utilize Calculate the 2(M-1)×2(M-1) order pseudolinear observation error covariance matrix.

[0122] Step 10: Based on and Calculate the optimal weighting matrix of order 2M×2M.

[0123] Step 11: Calculate the 3×3 matrix sequentially. 3 3×1 vectors as well as 6 scalars as well as And set the initial value of the Lagrange multiplier λ3.

[0124] Step 12: Calculate the Lagrange multiplier vector Bivariate polynomial coefficients as well as And calculate their estimates with respect to the Lagrange multipliers. first derivative as well as

[0125] Step 13: Calculate the bivariate polynomial coefficient vector with respect to the ground radiation source location vector u. as well as And calculate their estimates with respect to the Lagrange multipliers. first derivative as well as

[0126] Step 14: Construct two simultaneous parameters related to the distance β between the ground-based radiation source and the host star, and the rate of change of that distance. Given a bivariate polynomial equation, calculate the coefficients of the above bivariate polynomial equation. and And calculate and Regarding the estimates of the Lagrange multipliers first derivative and

[0127] Step 15: Construct a univariate seventh-degree polynomial equation for the distance β between the ground-based radiation source and the host star, based on... and The coefficients of the univariate seventh-degree polynomial equation are calculated using the linear convolution operation of sequences. The positive real roots of the univariate seventh-degree polynomial equation are then found using the power method (denoted as ). ).

[0128] Step 16: Calculate the rate of change of distance between the ground-based radiation source and the primary star. The iteration value (denoted as) ).

[0129] Step 17: Calculate the vector Regarding the estimates of the Lagrange multipliers first derivative

[0130] Step 18: Calculate the iterative value of the ground radiation source location vector. And calculate Regarding the estimates of the Lagrange multipliers first derivative

[0131] Step 19: Estimate the Lagrange multipliers Perform Newton iterations to obtain the updated value. If update volume Then proceed to step 20; otherwise, let Proceed to step 12.

[0132] Step 20: Let as well as like If the calculation stops, then stop; otherwise, update the iteration index k:=k+1 and go to step 8.

[0133] Further, in step 1, M communication satellites are selected to locate the ground radiation source, where the first satellite is the primary satellite and the rest are secondary satellites. The position vector of the ground radiation source is denoted as u, and the position and velocity vectors of the m-th satellite are denoted as s, respectively. m and They used this to obtain TDOA observations of the ground-based radiation source signal reaching the m-th (2≤m≤M) satellite and reaching the host star. and FDOA observations TDOA observations Multiplying by the signal propagation speed c yields the distance differential measurement. The corresponding expression is

[0134]

[0135] In the formula, s1 is the position vector of the first satellite, and Δr m1 This indicates the distance difference observation error;

[0136] Then the FDOA observations Multiplying by the signal propagation speed c and then dividing by the signal transmission frequency f0 yields the observable rate of change of distance difference. The corresponding expression is

[0137]

[0138] In the formula Let this be the velocity vector of the first satellite. This represents the observation error of the rate of change of distance difference.

[0139] Furthermore, in step 2, distance difference observation is used. Observations on the rate of change of distance difference Construct an M×M scalar product matrix and M×M scalar product matrix The elements in these two matrices are respectively

[0140]

[0141] In the formula This represents the distance between the m1-th satellite and the m2-th satellite; This represents the inner product of the position vector difference and the velocity vector difference between the m1-th satellite and the m2-th satellite;

[0142] Furthermore, in step 3, distance difference observation is used. Observations on the rate of change of distance difference Construct an M×4 order observation matrix and M×4 order observation matrix The elements in these two matrices are respectively

[0143]

[0144] Furthermore, in step 4, the observation matrix is ​​used. and observation matrix Construct an M×5 pseudo-inverse matrix and M×5 pseudo-inverse matrix The formulas for calculating these two matrices are as follows:

[0145]

[0146] In the formula 1 M Represents an M×1 column vector of all 1s; 0 M This represents an M×1 column vector consisting entirely of zeros.

[0147] Furthermore, in step 5, the scalar product matrix is ​​used. scalar product matrix and pseudo-inverse matrix and pseudo-inverse matrix Construct a 2M×6 pseudo-linear observation matrix The formula for calculating this matrix is:

[0148]

[0149] In the formula This represents the 5th column vector in the 5×5 identity matrix I5; Representation matrix The first column vector in; Represented by matrix The matrix formed by columns 2 to 4 in the matrix; Representation matrix The 5th column vector in; Representation matrix The 6th column vector in.

[0150] Furthermore, in step 6, four first-order perturbation matrices are calculated sequentially. as well as The formulas for calculating these four matrices are as follows:

[0151]

[0152] In the formula I M-1 and I M These represent (M-1)×(M-1) and M×M identity matrices, respectively; the diag function is used to create diagonal matrices; 0 M-1 Represents an (M-1)×1 order all-zero column vector; O (M-1)×(M-1) and O (3M+1)×(M-1) Let represent (M-1)×(M-1) and (3M+1)×(M-1) matrices, respectively, all zeros.

[0153] and Let Π represent the distance difference observation vector and the distance difference change rate observation vector, respectively; Π represents the permutation matrix, which satisfies the following condition: The vec function is used to convert a matrix into a column vector in column order.

[0154] I4 represents a 4×4 identity matrix, and 04 represents a 4×1 column vector of all zeros.

[0155] Furthermore, in step 7, the iteration index k:=0 is set, the iteration threshold value δ is set, and the initial iteration value is calculated. The corresponding calculation formula is:

[0156]

[0157] in express The column vector consisting of the first 4 elements, express The last element.

[0158] This yields the initial values ​​for the iterative iteration of the ground radiation source location vector. Where I3 represents a 3×3 identity matrix; 03 represents a 3×1 column vector of all zeros.

[0159] Furthermore, in step 8, a 2M×2(M-1) order observation error coefficient matrix is ​​calculated. The corresponding calculation formula is:

[0160]

[0161] In the formula and O M×5M Representing M×M 2 The order of M×5M is an all-zero matrix.

[0162] Furthermore, in step 9, the matrix... Singular value decomposition yields

[0163]

[0164] In the formula Represents a 2M×2M left singular matrix; Denotes a 2(M-1)×2(M-1) order right singular matrix; Represents a 2(M-1)×2(M-1) order singular value diagonal matrix; O 2×2(M-1) This represents a 2×2(M-1) order all-zero matrix.

[0165] Then, using a 2M×2(M-1) order column orthogonal matrix Calculate the 2(M-1)×2(M-1) order pseudolinear observation error covariance matrix. The corresponding calculation formula is:

[0166]

[0167] In the formula, E represents the 2(M-1)×2(M-1) order observation error covariance matrix.

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

[0169]

[0170] Furthermore, in step 11, a 3×3 matrix is ​​calculated sequentially. 3 3×1 vectors as well as 6 scalars as well as The corresponding calculation formulas are as follows:

[0171]

[0172] Then set the initial value of the Lagrange multiplier λ3.

[0173] Furthermore, in step 12, the calculation of the Lagrange multiplier vectors is performed. Bivariate polynomial coefficients as well as The corresponding calculation formula is:

[0174]

[0175]

[0176] In the formula Where e = 0.081819790992113 represents the first eccentricity of the Earth, and the diag function is used to create a diagonal matrix;

[0177] Then they calculated their estimates of the Lagrange multipliers. first derivative as well as The corresponding calculation formula is:

[0178]

[0179] In the formula

[0180]

[0181] Further, in step 13, the bivariate polynomial coefficient vector with respect to the ground radiation source location vector u is calculated. as well as The corresponding calculation formula is:

[0182]

[0183] Then they calculated their estimates of the Lagrange multipliers. first derivative as well as The corresponding calculation formula is:

[0184]

[0185] Furthermore, in step 14, two parameters are constructed simultaneously regarding the distance β between the ground-based radiation source and the host star, and the rate of change of that distance. Given a bivariate polynomial equation, calculate the coefficients of the two polynomials. and The corresponding calculation formula is:

[0186]

[0187]

[0188] In the formula

[0189]

[0190] In the formula, the symbols <·>1 and <·>2 represent the first and second elements in the vector, respectively.

[0191] Then calculate separately. and Regarding the estimates of the Lagrange multipliers first derivative and The corresponding calculation formula is:

[0192]

[0193] Furthermore, in step 15, a univariate seventh-degree polynomial equation concerning the distance β between the ground radiation source and the primary star is constructed, and its polynomial coefficients are calculated. The corresponding expression is

[0194]

[0195] polynomial coefficients It can be obtained through linear convolution of sequences, and the corresponding calculation formula is:

[0196]

[0197] In the formula, * represents the linear convolution operation of sequences.

[0198] Then, the positive real roots of the equation are found using the power method, and denoted as .

[0199] Furthermore, in step 16, the rate of change of distance between the ground radiation source and the primary star is calculated. The iteration value, and denoted as The corresponding calculation formula is:

[0200]

[0201] Furthermore, in step 17, the vector is calculated. Regarding the estimates of the Lagrange multipliers first derivative The corresponding calculation formula is:

[0202]

[0203] In the formula

[0204]

[0205]

[0206] Furthermore, in step 18, the iterative value of the ground radiation source location vector is calculated. The corresponding calculation formula is:

[0207]

[0208] Then calculate Regarding the estimates of the Lagrange multipliers first derivative The corresponding calculation formula is:

[0209]

[0210] Furthermore, in step 19, the estimated values ​​of the Lagrange multipliers are... Perform Newton iteration update, the corresponding calculation formula is:

[0211]

[0212] In the formula r e =6378.137km represents the Earth's equatorial radius. (If updated...) Then proceed to step 20; otherwise, let Proceed to step 12.

[0213] Furthermore, in step 20, let as well as like If the calculation stops, then stop; otherwise, update the iteration index k:=k+1 and go to step 8.

[0214] As a concrete example, consider locating a radiation source on the Earth's surface with longitude of 133.14° and latitude of 32.20°. Its position vector in the ECEF coordinate system is u = [-3694.03942.03379.2]. T (km), there are 7 existing communication satellites performing TDOA / FDOA joint positioning of the radiation source. The latitude, longitude and altitude of these 7 communication satellites are shown in Table 1, and the movement speed of these 7 communication satellites is shown in Table 2.

[0215] Table 1. Latitude, longitude, and altitude of the 7 communication satellites

[0216]

[0217] Table 2. Movement speeds of 7 communication satellites

[0218]

[0219] The observation error covariance matrix is ​​set as Where I6 represents a 6×6 identity matrix, 16×6 Represents a 6×6 matrix of all 1s, O 6×6 Let σ represent a 6×6 matrix of all zeros, and let σ represent the standard deviation of the observation error.

[0220] First, let the standard deviation of the observation error σ be set to σ = 0.5. Figure 2 The scatter plot of the ground radiation source location results and the elliptic curve of the location error (XY plane coordinates in the ECEF coordinate system) are presented. Figure 3 The scatter plot of the ground radiation source location results and the elliptic curve of the location error (YZ plane coordinates in the ECEF coordinate system) are presented. From Figure 2 and Figure 3 As can be seen from the diagram, the shape of the scatter plot of the positioning result of the positioning method disclosed in this patent is consistent with the shape of the positioning error ellipse, and it corresponds to a large-area ellipse with a high probability and a small-area ellipse with a low probability, thus verifying the effectiveness of the new method.

[0221] Next, change the value of the standard deviation σ of the observation error. Figure 4 The curve of the root mean square error of the ground radiation source location vector estimation as a function of the standard deviation σ of the observation error is given; then, setting the standard deviation σ of the observation error to σ = 0.5, and changing the longitude of the ground radiation source, Figure 5 The root mean square error of the ground radiation source location vector estimation is given as a function of the longitude of the ground radiation source.

[0222] from Figure 4 and Figure 5 It can be seen from the following: (1) The multi-satellite TDOA / FDOA joint positioning method for ground radiation sources based on Earth ellipsoid constraints disclosed in this patent can asymptotically approximate the corresponding Cramer-Rao boundary for estimating the position vector of the ground radiation source, thus verifying the asymptotic statistical optimality of the new method; (2) Compared with the multi-satellite TDOA / FDOA joint positioning method that does not utilize Earth ellipsoid constraints, the positioning accuracy of the new method is significantly improved, and the performance gain obtained will increase with the increase of the standard deviation of the observation error and the longitude of the ground radiation source.

[0223] The above description is only a preferred embodiment of the present invention. It should be noted that those skilled in the art can make several improvements and modifications without departing from the principle of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A method for TDOA / FDOA joint positioning of ground radiation source multi-satellites based on constraints of the Earth ellipsoid, characterized in that, comprising: Step 1: select M communication satellites to locate the ground radiation source, and obtain TDOA observation of the ground radiation source signal arriving at the mth satellite and arriving at the first satellite using the M communication satellites and FDOA observation and further obtain range difference observation and range difference rate observation using the TDOA observation and the FDOA observation and range difference rate observation Step 2: Utilize Constructing the M x M order scalar product matrix Utilizing and Constructing the M x M order scalar product matrix Step 3: Utilizing Constructing an M x 4 order observation matrix Utilizing Constructing an M x 4 order observation matrix Step 4: Utilizing Constructing the M x 5 order pseudo-inverse matrix Utilizing Constructing the M x 5 order pseudo-inverse matrix Step 5: Utilizing and Constructing a 2M x 6 order pseudo-linear observation matrix Step 6: Based on one or more parameters in and the four first order perturbation matrices are computed sequentially and Step 7: Let iteration index k: = 0, set iteration threshold value δ, and based on Calculate iteration initial value Thus, the iteration initial value of the ground radiation source position vector is obtained Step 8: Based on and Computing a 2M x 2(M-1) order observation error coefficient matrix Step 9: On the matrix perform singular value decomposition to determine a 2M x 2(M-1) order column orthogonal matrix and use to calculate a 2(M-1) x 2(M-1) order pseudo-linear observation error covariance matrix Step 10: Based on and Computing the optimal weighting matrix of order 2M x 2M Step 11: Based on and 3x3 matrix 3 3x1 vectors and 6 scalars and Set initial value of Lagrange multiplier λ3 Step 12: Based on and One or more parameters in the vector are used to compute the Lagrange multiplier vector. Bivariate polynomial coefficients as well as And calculate as well as Regarding the estimates of the Lagrange multipliers first derivative as well as Step 13: Compute the bivariate polynomial coefficient vector and based on one or more parameters in and and compute and the first derivative of and ​ Step 14: Constructing two bivariate polynomial equations simultaneously on the distance β and the rate of change of the distance between the ground radiating source and the 1st satellite based on one or more parameters in and ​​​​​​​ Step 15: Construct a monic seventh-degree polynomial equation in terms of the distance β between the ground radiating source and the first satellite, based on and Calculate the coefficients of the monic seventh-degree polynomial equation using sequential linear convolution operations and solve for the positive real root of the monic seventh-degree polynomial equation using the power method Step 16: Based on and Calculate the rate of change of the distance between the ground radiating source and the 1st satellite the iteration value of Step 17: Compute the vector With respect to The first derivative of Step 18: Based on and Computing an iterative value of the ground radiance source position vector and computing the first derivative of with respect to Step 19: Perform Newton iteration update on to get updated value If the update is go to Step 20; otherwise let and go to Step 12; Step 20: Let and If then stop the computation; otherwise update the iteration index k := k + 1 and go to Step 8.

2. The TDOA / FDOA joint positioning method based on the constraint of the Earth ellipsoid according to claim 1, characterized in that, In step 1, the distance difference observation is computed as follows wherein denotes the range difference observation of the mth satellite to the 1st satellite, denotes the TDOA observation of the ground radiating source signal arriving at the mth satellite to the 1st satellite, c is the signal propagation speed, u is the ground radiating source position vector, s m is the position vector of the mth satellite, 2≤m≤M, s1 is the position vector of the 1st satellite, Δr m1 denotes the range difference observation error of the mth satellite to the 1st satellite; The distance difference change rate observation is calculated in the following way wherein denotes the rate of change of the range difference between the mth satellite and the 1st satellite, denotes the FDOA observation of the ground radiated source signal arriving at the mth satellite and arriving at the 1st satellite, f0 denotes the signal transmitting frequency, is the velocity vector of the mth satellite, is the velocity vector of the 1st satellite, denotes the observation error of the rate of change of the range difference between the mth satellite and the 1st satellite.

3. The method of claim 1, wherein the method is a TDOA / FDOA joint positioning method based on the constraint of the Earth ellipsoid, and the method comprises the following steps: In step 2, the M x M order scalar product matrix is constructed in the following way and the M x M order scalar product matrix In the formula denotes the distance between the m1th satellite and the m2th satellite, and 1≤m1; m2≤M; denotes the inner product of the position vector difference and the velocity vector difference between the m1th satellite and the m2th satellite; 4. The TDOA / FDOA joint positioning method based on the constraint of the Earth ellipsoid according to claim 1, characterized in that, In step 3, the M x 4 order observation matrix is constructed in the following manner and the M x 4 order observation matrix where s1is the position vector of the first satellite, s m is the position vector of the mth satellite, 2≤m≤M, is the velocity vector of the first satellite, is the velocity vector of the mth satellite.

5. The TDOA / FDOA joint positioning method based on the constraint of the Earth ellipsoid according to claim 1, characterized in that, In step 4, the M x 5 order pseudo-inverse matrix is constructed in the following manner and the M x 5 order pseudo-inverse matrix where 1 M denotes an M x 1 all-one column vector; 0 M denotes an M x 1 all-zero column vector.

6. The TDOA / FDOA joint positioning method based on the constraint of the Earth ellipsoid according to claim 1, characterized in that, In step 5, the 2M x 6 pseudo-linear observation matrix is constructed in the following manner In the formula Represents the 5th column vector in the 5×5 identity matrix I5; 0 M Represents an M×1 column vector consisting entirely of zeros; Representation matrix The first column vector in; Represented by matrix The matrix formed by columns 2 to 4 in the matrix; Representation matrix The 5th column vector in; Representation matrix The 6th column vector in.

7. The TDOA / FDOA joint positioning method based on the constraint of the Earth ellipsoid according to claim 1, characterized in that, In step 6, the four first order perturbation matrices are calculated as follows and wherein I M-1 and I M respectively represent (M-1) x (M-1) and M x M identity matrices; diag function is used to create a diagonal matrix; 0 M-1 denotes an (M-1) x 1 order all-zero column vector; O (M-1)×(M-1) and O (3M+1)×(M-1) denotes an (M-1) x (M-1) order and (3M+1) x (M-1) order all-zero matrix, respectively; and denote the range difference observation vector and the range difference rate observation vector, respectively; Π denotes a permutation matrix satisfying The vec function is used to convert into a column vector in the order of matrix rows; 1 M denotes an M x 1 order all-one column vector; I4represents a 4 x 4 order identity matrix, and 04represents a 4 x 1 order all-zero column vector.

8. The TDOA / FDOA joint positioning method based on the constraint of the Earth ellipsoid according to claim 6, characterized in that, In step 7, the iteration initial value is calculated as follows wherein represents a column vector consisting of the first 4 elements of represents the last 1 element of Further, the iterative initial value of the ground radiation source position vector is obtained I3 represents a 3x3 order unit matrix; 03 represents a 3x1 order all-0 column vector.

9. The TDOA / FDOA joint positioning method based on the constraint of the Earth ellipsoid according to claim 8, characterized in that, In step 8, the 2M x 2(M-1) order observation error coefficient matrix is calculated in the following manner wherein and O M×5M denote M x M 2 and M x 5M order all-zero matrices, respectively; denotes the 5th column vector in the 5 x 5 order identity matrix I5.

10. The TDOA / FDOA joint positioning method based on the constraint of the Earth ellipsoid for a ground radiation source multi-satellite according to claim 6, characterized in that, In step 9, singular value decomposition is performed on the matrix to obtain wherein denotes a 2M x 2M left-singular matrix; denotes a 2(M-1) x 2(M-1) right-singular matrix; denotes a 2(M-1) x 2(M-1) singular-value diagonal matrix; O 2×2(M-1) denotes a 2 x 2(M-1) all-zero matrix; Then, the following is used Computing a 2(M-1) x 2(M-1) order pseudo-linear observation error covariance matrix where E represents a 2(M-1) x 2(M-1) order observation error covariance matrix; In step 10, the 2M x 2M optimal weight matrix is calculated in the following manner In step 11, the 3x3 matrix is calculated sequentially in the following way 3 3x1 vectors and 6 scalars and In step 12, the following is calculated and wherein s1is the position vector of the first satellite, is the velocity vector of the first satellite; e = 0.081819790992113 represents the first eccentricity of the Earth, and the diag function is used to create a diagonal matrix. Then the following is calculated and where In step 13, the following is calculated and Then the following is calculated and In step 14, a bivariate polynomial equation is constructed in terms of the distance β between the ground-based radiating source and the first satellite and the rate of change of the distance β as follows: β + β' = a + bβ + cβ2+ dβ3+ eβ4+ fβ5+ gβ6+ hβ7+ iβ8+ jβ9+ kβ10+ lβ11+ mβ12+ nβ13+ oβ The coefficients of the above binary polynomial equation are calculated in the following manner: where the symbols <·>1 and <·>2 represent the 1st and 2nd elements in the vector, respectively; Then the following are calculated separately and with respect to the first derivative and In the step 15, a monomial 7th order polynomial equation with respect to the distance β between the ground radiation source and the 1st satellite is constructed in the following manner: In the formula Obtained by a sequence linear convolution operation: where * represents a sequential linear convolution operation; The rate of change of the distance between the ground radiating source and the first satellite is calculated in step 16 in the following manner of the iteration value In said step 17, the vector is calculated in the following way With respect to the first derivative of where In step 18, the iterative value of the ground radiating source position vector is calculated in the following manner Then the following is calculated with respect to the Lagrange multiplier estimate of the first derivative In step 19, the Lagrange multiplier estimate is updated in the following manner Newton iteration update is performed: where r e = 6378.137 km represents the Earth equatorial radius.

Images