CN110673196B - Time difference positioning method based on multidimensional calibration and polynomial root finding - Google Patents

Time difference positioning method based on multidimensional calibration and polynomial root finding Download PDF

Info

Publication number
CN110673196B
CN110673196B CN201910893455.3A CN201910893455A CN110673196B CN 110673196 B CN110673196 B CN 110673196B CN 201910893455 A CN201910893455 A CN 201910893455A CN 110673196 B CN110673196 B CN 110673196B
Authority
CN
China
Prior art keywords
matrix
vector
order
time difference
root
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201910893455.3A
Other languages
Chinese (zh)
Other versions
CN110673196A (en
Inventor
王鼎
唐涛
张莉
王成
杨宾
魏帅
杨泽宇
李建兵
张光伟
韩博臣
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Information Engineering University of PLA Strategic Support Force
Original Assignee
Information Engineering University of PLA Strategic Support Force
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Information Engineering University of PLA Strategic Support Force filed Critical Information Engineering University of PLA Strategic Support Force
Priority to CN201910893455.3A priority Critical patent/CN110673196B/en
Publication of CN110673196A publication Critical patent/CN110673196A/en
Application granted granted Critical
Publication of CN110673196B publication Critical patent/CN110673196B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01T—MEASUREMENT OF NUCLEAR OR X-RADIATION
    • G01T1/00—Measuring X-radiation, gamma radiation, corpuscular radiation, or cosmic radiation
    • G01T1/29—Measurement performed on radiation beams, e.g. position or section of the beam; Measurement of spatial distribution of radiation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Data Mining & Analysis (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Theoretical Computer Science (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Algebra (AREA)
  • Computing Systems (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • High Energy & Nuclear Physics (AREA)
  • Molecular Biology (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Radar Systems Or Details Thereof (AREA)
  • Position Fixing By Use Of Radio Waves (AREA)

Abstract

The invention belongs to the technical field of radiation source positioning, and discloses a time difference positioning method based on multi-dimensional calibration and polynomial root finding, which comprises the following steps: step 1: respectively counting the time difference of the radiation source signal reaching the mth sensor and the 1 st sensor, and multiplying the time difference by the signal propagation speed to obtain the observed distance difference
Figure DDA0002209502900000011
Wherein M is the number of sensors; step 2: observed quantity by distance difference
Figure DDA0002209502900000012
Constructing a scalar product matrix, and constructing a multi-dimensional calibration pseudo linear equation through the scalar product matrix; and step 3: and solving a multi-dimensional calibration pseudo linear equation by using a Newton method, selecting a real root, eliminating a virtual root, and obtaining a position vector of the radiation source through the real root. The invention improves the accuracy of multi-sensor time difference positioning.

Description

Time difference positioning method based on multidimensional calibration and polynomial root finding
Technical Field
The invention belongs to the technical field of radiation source positioning, and particularly relates to a time difference positioning method based on multi-dimensional calibration and polynomial root finding.
Background
As is well known, radiation source location technology plays an important role in many industrial and information fields, such as target monitoring, navigation telemetry, seismic surveying, radio astronomy, emergency assistance, safety management, etc. The basic process of radiation source location is to extract parameters (also called location observation) related to the position of the radiation source from the electromagnetic signal, and then to use the parameters to calculate the radiation source position information. Observations for radiation source localization relate to multi-domain parameters such as space, time, frequency, energy, etc., where time difference of arrival (which may be equivalently distance difference) is a type of observation that is applied more frequently. The multi-sensor time difference positioning technology is characterized in that positioning is carried out through the arrival time difference of radiation source signals acquired by a plurality of sensors, one hyperboloid (line) is determined according to the arrival time difference of a target radiation source to two different sensors, and the radiation source coordinates can be obtained by intersecting a plurality of hyperboloids (lines). With the continuous development of modern communication technology and time difference measurement technology, the time difference positioning technology has become one of the most mainstream radiation source positioning means.
Based on the algebraic characteristics of the equation of time difference observation, researchers at home and abroad propose many time difference positioning methods with excellent performance, including iterative methods (Yang K, An J P, Bu X Y, Sun G C. constrained total least squares-precision location estimation time-difference-of-arrival measurements [ J ]. IEEE Transactions on temporal Technology,2010,59(3):1558 × 1562.) (Qu X M, Xie L H, Tan W R. objective constrained weighted prediction results localization A and FDmeasure measurements [ J ]. IEEE Transactions on Signal Processing,2017,65(15): 3) and TDO Processing methods [ IEEE transaction ] for solving the equation of time difference estimation J. TDO. approximate, TDO. simulation, I, J.), 2009,57(12): 4598-. In the closed solution positioning method, a multidimensional calibration time difference positioning method is provided, and the method obtains a pseudo linear equation related to the position of a radiation source by constructing a scalar product matrix and obtains a closed solution of the position parameter of the radiation source based on the pseudo linear equation. However, this kind of method fails to utilize quadratic equation constraints to which the augmented vector obeys, and therefore its positioning accuracy has not yet reached asymptotically optimal solution.
Disclosure of Invention
The invention provides a time difference positioning method based on multi-dimensional calibration and polynomial root solving, aiming at the problem of poor time difference positioning accuracy of a multi-sensor, fully considering quadratic equation constraint obeyed by an augmented vector, converting the time difference positioning problem into the polynomial root solving problem, and obtaining a radiation source position vector by solving the polynomial root.
In order to achieve the purpose, the invention adopts the following technical scheme:
a time difference positioning method based on multi-dimensional calibration and polynomial root solving comprises the following steps:
step 1: respectively counting the time difference of the radiation source signal reaching the mth sensor and the 1 st sensor, and multiplying the time difference by the signal propagation speed to obtain the observed distance difference
Figure BDA0002209502880000021
Wherein M is the number of sensors;
step 2: observed quantity by distance difference
Figure BDA0002209502880000022
Constructing a scalar product matrix, and constructing a multi-dimensional calibration pseudo linear equation through the scalar product matrix;
and step 3: and solving a multi-dimensional calibration pseudo linear equation by using a Newton method, selecting a real root, eliminating a virtual root, and obtaining a position vector of the radiation source through the real root.
Further, the step 2 comprises:
step 2.1: using sensor position vector sm}1≤m≤MSum-distance difference observed quantity
Figure BDA0002209502880000023
Constructing an (M +1) × (M +1) order distance matrix D:
Figure BDA0002209502880000024
in the formula dmn＝||sm-sn||2(1≤m、n≤M)，dmnRepresents the distance between the m-th sensor and the n-th sensor;
step 2.2: calculating an (M +1) × (M +1) order scalar product matrix W using the distance matrix D:
Figure BDA0002209502880000031
in the formula
Figure BDA0002209502880000032
Is a transformation matrix, in which IM+1Represents an identity matrix of order (M +1) × (M +1), 1M+1Representing a full 1 vector of length M + 1;
step 2.3: by sensor position vector sm}1≤m≤MSum-distance difference observed quantity
Figure BDA0002209502880000033
Constructing a matrix G of (M +1) × 4 order, constructing a matrix G of (M +1) × 5 order
Figure BDA0002209502880000034
t1Representing the 1 st column vector in the matrix T, T2A matrix formed by vectors of 2 nd to 5 th columns in the matrix T is represented:
Figure BDA0002209502880000035
in the formula (I), the compound is shown in the specification,
Figure BDA0002209502880000036
in the formula
Figure BDA0002209502880000037
Is the distance difference observed for the 1 st sensor, which is constantly zero;
step 2.4: setting the iteration index k:equalto 0, setting an iteration threshold value delta, and calculating an iteration initial value through a scalar product matrix W and a matrix T
Figure BDA0002209502880000038
Step 2.5: observed quantity by matrix T and distance difference
Figure BDA0002209502880000039
Calculating an (M +1) × (M-1) order matrix
Figure BDA00022095028800000310
Calculating an (M +1) × (M-1) order matrix from the scalar product matrix W and the matrix G
Figure BDA00022095028800000311
Step 2.6: by means of a matrix
Figure BDA00022095028800000312
And
Figure BDA00022095028800000313
calculating an (M +1) × (M-1) order matrix
Figure BDA0002209502880000041
And to the matrix
Figure BDA0002209502880000042
Making singular valuesDecomposition of
Figure BDA0002209502880000043
Wherein H(k)Is an orthogonal matrix of (M +1) × (M-1) order, V(k)Is an orthogonal matrix of (M-1) × (M-1) order(k)Is a diagonal matrix of (M-1) × (M-1) order, the diagonal elements of which are matrices
Figure BDA0002209502880000044
The singular value of (a);
step 2.7: after decomposition by singular values
Figure BDA0002209502880000045
Calculating a weighting matrix omega(k)＝Σ(k)V(k)TEV(k)Σ(k)TWherein E is a covariance matrix of range difference observation errors;
step 2.8: by scalar product matrix W, matrix T, H(k)And omega(k)Calculating the matrix phi in turn(k)Sum vector
Figure BDA0002209502880000046
Figure BDA0002209502880000047
Step 2.9: through phi(k)Construction matrix
Figure BDA0002209502880000048
For matrix
Figure BDA0002209502880000049
And (3) carrying out characteristic value decomposition:
Figure BDA00022095028800000416
wherein, O3×1Representing a 3-dimensional all-zero column vector, O1×3Representing a 3-dimensional all-zero row vector, matrix P(k)Is a matrix of
Figure BDA00022095028800000410
Is determined by the feature vector of (a),
Figure BDA00022095028800000411
is a matrix
Figure BDA00022095028800000412
A characteristic value of (d);
step 2.10: through phi(k)、
Figure BDA00022095028800000415
P(k)And the 1 st sensor position vector s1The vectors are calculated sequentially as follows
Figure BDA00022095028800000413
And
Figure BDA00022095028800000414
Figure BDA0002209502880000051
step 2.11: by passing
Figure BDA0002209502880000052
And
Figure BDA0002209502880000053
sequentially calculating scalars
Figure BDA0002209502880000054
Figure BDA0002209502880000055
And
Figure BDA0002209502880000056
step 2.12: by passing
Figure BDA0002209502880000057
s1And
Figure BDA0002209502880000058
sequentially calculating scalars
Figure BDA0002209502880000059
And
Figure BDA00022095028800000510
step 2.13: by passing
Figure BDA00022095028800000511
And
Figure BDA00022095028800000512
calculating a scalar quantity
Figure BDA00022095028800000513
Step 2.14: by passing
Figure BDA00022095028800000514
And constructing a multi-dimensional calibration pseudo linear equation.
Further, the step 3 comprises:
step 3.1: solving a multi-dimensional calibration pseudo linear equation by using a Newton method, selecting a real root, and removing a pseudo root;
step 3.2: updating iteration index k: ═ k +1, and calculating iteration result
Figure BDA00022095028800000515
If it is
Figure BDA00022095028800000516
Stopping iteration, otherwise, turning to the step 2.5;
step 3.3: and determining the position vector of the radiation source by using the final iteration convergence result.
Further, the matrix
Figure BDA00022095028800000517
The calculation formula of (2) is as follows:
Figure BDA0002209502880000061
in the formula
Figure BDA0002209502880000062
Representing a vector
Figure BDA0002209502880000063
The m-th element of (1);
the matrix
Figure BDA0002209502880000064
The calculation formula of (2) is as follows:
Figure BDA0002209502880000065
in the formula
Figure BDA0002209502880000066
Representing a vector
Figure BDA0002209502880000067
The m-th element of (1);
Figure BDA0002209502880000068
representing a vector
Figure BDA0002209502880000069
The m-th element of (1); o is4×(M-1)Represents an all-zero matrix of 4 (M-1) th order.
Further, the
Figure BDA00022095028800000610
And
Figure BDA00022095028800000611
is calculated byThe formula is as follows:
Figure BDA0002209502880000071
wherein < >jRepresenting the jth component in the vector.
Further, the
Figure BDA0002209502880000072
And
Figure BDA0002209502880000073
the calculation formula of (2) is as follows:
Figure BDA0002209502880000074
further, the scalar quantity
Figure BDA0002209502880000075
The calculation formula of (2) is as follows:
Figure BDA0002209502880000081
further, the step 3.1 comprises:
the multidimensional calibration pseudo linear equation is as follows:
Figure BDA0002209502880000082
let the root of the equation be
Figure BDA0002209502880000083
Selecting a real root by
Figure BDA0002209502880000084
Figure BDA0002209502880000091
In the formula
Figure BDA0002209502880000092
For the distance-difference observation vector,
Figure BDA0002209502880000093
indicating the use of the jth eigenvalue
Figure BDA0002209502880000094
The obtained position vector of the radiation source,
Figure BDA0002209502880000095
further, the
Figure BDA0002209502880000096
The calculation formula of (2) is as follows:
Figure BDA0002209502880000097
further, the step 3.3 comprises: using final iterative convergence results
Figure BDA0002209502880000098
The first 3 components of which serve as the position vector of the radiation source.
Compared with the prior art, the invention has the following beneficial effects:
the method converts the time difference positioning problem into a polynomial root solving problem by fully utilizing quadratic equation constraint obeyed by the augmentation vector, and obtains the radiation source position vector by solving the polynomial root.
Drawings
Fig. 1 is a basic flowchart of a time difference positioning method based on multidimensional calibration and polynomial root finding according to an embodiment of the present invention;
FIG. 2 is a plot of root mean square error of radiation source position estimate versus standard deviation σ of range observation error;
FIG. 3 is a plot of the RMS error of the source position estimate as a function of the parameter c;
fig. 4 is a graph of root mean square error of radiation source position estimation (σ ═ 1) as a function of parameter c;
fig. 5 shows the variation of the root mean square error of the radiation source position estimate with the parameter c (σ ═ 2).
Detailed Description
The invention is further illustrated by the following examples in conjunction with the accompanying drawings:
as shown in fig. 1, a time difference positioning method based on multidimensional scaling and polynomial root finding includes:
step S11: respectively counting the time difference of the radiation source signal reaching the mth sensor and the 1 st sensor, and multiplying the time difference by the signal propagation speed to obtain the observed distance difference
Figure BDA0002209502880000101
Wherein M is the number of sensors;
step S12: observed quantity by distance difference
Figure BDA0002209502880000102
Constructing a scalar product matrix, and constructing a multi-dimensional calibration pseudo linear equation through the scalar product matrix;
step S13: and solving a multi-dimensional calibration pseudo linear equation by using a Newton method, selecting a real root, eliminating a virtual root, and obtaining a position vector of the radiation source through the real root.
Specifically, in step S11, M sensors are placed in space and used to time-difference locate the radiation source. The radiation source position vector is u, and the position vector of the mth sensor is sm＝[xm ym zm]T(1. ltoreq. m.ltoreq.M), with which it is possible to obtain the arrival of the radiation source signal at the M (2. ltoreq. m.ltoreq.M) th sensor and at the 1 st sensorTime difference
Figure BDA0002209502880000103
Will be time difference
Figure BDA0002209502880000104
Multiplying by the signal propagation velocity to obtain the observed distance difference
Figure BDA0002209502880000105
The expression is as follows:
Figure BDA0002209502880000106
in the formula ofmRepresenting the range difference observation error.
Specifically, the step S12 includes:
step S12.1: using sensor position vector sm}1≤m≤MSum-distance difference observed quantity
Figure BDA0002209502880000107
Constructing an (M +1) × (M +1) order distance matrix D:
Figure BDA0002209502880000108
in the formula dmn＝||sm-sn||2(1≤m、n≤M)，dmnRepresents the distance between the m-th sensor and the n-th sensor;
step S12.2: calculating an (M +1) × (M +1) order scalar product matrix W using the distance matrix D:
Figure BDA0002209502880000109
in the formula
Figure BDA0002209502880000111
LMIs a transformation matrix, in which IM+1Represents an identity matrix of order (M +1) × (M +1), 1M+1Is expressed as length MA full 1 vector of + 1;
step S12.3: by sensor position vector sm}1≤m≤MSum-distance difference observed quantity
Figure BDA0002209502880000112
Constructing a matrix G of (M +1) × 4 order, constructing a matrix G of (M +1) × 5 order
Figure BDA0002209502880000113
t1Representing the 1 st column vector in the matrix T, T2A matrix formed by vectors of 2 nd to 5 th columns in the matrix T is represented:
Figure BDA0002209502880000114
in the formula (I), the compound is shown in the specification,
Figure BDA0002209502880000115
in the formula
Figure BDA0002209502880000116
Is the distance difference observed for the 1 st sensor, which is constantly zero;
step S12.4: setting the iteration index k:equalto 0, setting an iteration threshold value delta, and calculating an iteration initial value through a scalar product matrix W and a matrix T
Figure BDA0002209502880000117
Step S12.5: observed quantity by matrix T and distance difference
Figure BDA0002209502880000118
Calculating an (M +1) × (M-1) order matrix
Figure BDA0002209502880000119
Calculating an (M +1) × (M-1) order matrix from the scalar product matrix W and the matrix G
Figure BDA00022095028800001110
Step S12.6: by means of a matrix
Figure BDA00022095028800001111
And
Figure BDA00022095028800001112
calculating an (M +1) × (M-1) order matrix
Figure BDA00022095028800001113
And to the matrix
Figure BDA00022095028800001114
Performing singular value decomposition
Figure BDA0002209502880000121
Wherein H(k)Is an orthogonal matrix of (M +1) × (M-1) order, V(k)Is an orthogonal matrix of (M-1) × (M-1) order(k)Is a diagonal matrix of (M-1) × (M-1) order, the diagonal elements of which are matrices
Figure BDA0002209502880000122
The singular value of (a);
step S12.7: after decomposition by singular values
Figure BDA0002209502880000123
Calculating a weighting matrix omega(k)＝Σ(k)V(k)TEV(k)Σ(k)TWherein E is a covariance matrix of range difference observation errors;
step S12.8: by scalar product matrix W, matrix T, H(k)And omega(k)Calculating the matrix phi in turn(k)Sum vector
Figure BDA0002209502880000124
Figure BDA0002209502880000125
Step S12.9: through phi(k)Construction matrix
Figure BDA0002209502880000126
For matrix
Figure BDA0002209502880000127
And (3) carrying out characteristic value decomposition:
Figure BDA0002209502880000128
wherein, O3×1Representing a 3-dimensional all-zero column vector, O1×3Representing a 3-dimensional all-zero row vector, matrix P(k)Is a matrix of
Figure BDA0002209502880000129
Is determined by the feature vector of (a),
Figure BDA00022095028800001210
is a matrix
Figure BDA00022095028800001211
A characteristic value of (d);
step S12.10: through phi(k)、
Figure BDA00022095028800001212
P(k)And the 1 st sensor position vector s1The vectors are calculated sequentially as follows
Figure BDA00022095028800001213
And
Figure BDA00022095028800001214
Figure BDA0002209502880000131
step S12.11: by passing
Figure BDA0002209502880000132
And
Figure BDA0002209502880000133
sequentially calculating scalars
Figure BDA0002209502880000134
Figure BDA0002209502880000135
And
Figure BDA0002209502880000136
step S12.12: by passing
Figure BDA0002209502880000137
s1And
Figure BDA0002209502880000138
sequentially calculating scalars
Figure BDA0002209502880000139
And
Figure BDA00022095028800001310
step S12.13: by passing
Figure BDA00022095028800001311
And
Figure BDA00022095028800001312
calculating a scalar quantity
Figure BDA00022095028800001313
Step S12.14: will be provided with
Figure BDA00022095028800001314
And constructing the multi-dimensional calibration pseudo linear equation as the polynomial coefficient of the multi-dimensional calibration pseudo linear equation which needs to be constructed finally.
Specifically, the step S13 includes:
step S13.1: solving a multi-dimensional calibration pseudo linear equation by using a Newton method, selecting a real root, and removing a pseudo root;
step S13.2: updating iteration index k: ═ k +1, and calculating iteration result
Figure BDA00022095028800001315
If it is
Figure BDA00022095028800001316
Stopping iteration, otherwise, turning to the step S12.5;
step S13.3: and determining the position vector of the radiation source by using the final iteration convergence result.
In particular, the matrix
Figure BDA00022095028800001317
The calculation formula of (2) is as follows:
Figure BDA0002209502880000141
in the formula
Figure BDA0002209502880000142
Representing a vector
Figure BDA0002209502880000143
The m-th element of (1);
the matrix
Figure BDA0002209502880000144
The calculation formula of (2) is as follows:
Figure BDA0002209502880000145
in the formula
Figure BDA0002209502880000146
Representing a vector
Figure BDA0002209502880000147
The m-th element of (1);
Figure BDA0002209502880000148
representing a vector
Figure BDA0002209502880000149
The m-th element of (1); o denotes an all-zero matrix, e.g. O4×(M-1)Represents an all-zero matrix of 4 (M-1) th order.
Specifically, the
Figure BDA00022095028800001410
And
Figure BDA00022095028800001411
the calculation formula of (2) is as follows:
Figure BDA0002209502880000151
wherein < >jRepresenting the jth component in the vector.
Specifically, the
Figure BDA0002209502880000152
And
Figure BDA0002209502880000153
the calculation formula of (2) is as follows:
Figure BDA0002209502880000154
in particular, the scalar quantity
Figure BDA0002209502880000155
The calculation formula of (2) is as follows:
Figure BDA0002209502880000161
in particular, said step S13.1 comprises:
the multidimensional calibration pseudo linear equation is as follows:
Figure BDA0002209502880000162
let the root of the equation be
Figure BDA0002209502880000163
Selecting a real root by
Figure BDA0002209502880000164
Figure BDA0002209502880000171
In the formula
Figure BDA0002209502880000172
For the distance-difference observation vector,
Figure BDA0002209502880000173
indicating the use of the jth eigenvalue
Figure BDA0002209502880000174
The obtained position vector of the radiation source,
Figure BDA0002209502880000175
specifically, the
Figure BDA0002209502880000176
The calculation formula of (2) is as follows:
Figure BDA0002209502880000177
in particular, the stepsStep S13.3 includes: using final iterative convergence results
Figure BDA0002209502880000178
The first 3 components of which serve as the position vector of the radiation source.
It is worth noting that in the process of constructing the multi-dimensional calibration pseudo-linear equation (polynomial) in the application, each matrix, each vector, each scalar, such as the matrix G, T, is constructed,
Figure BDA0002209502880000179
Φ(k)、
Figure BDA00022095028800001710
P(k)、
Figure BDA00022095028800001711
Equal, vector of
Figure BDA00022095028800001712
Etc. scalar quantity
Figure BDA00022095028800001713
And so on, are intermediate quantities set for simplifying mathematical expressions and calculation processes.
To verify the effect of the present invention, the following experiment was performed:
assuming a total of 7 observation stations (sensors) using TDOA (distance Difference information) to locate a radiation source, the position coordinates of the observation stations are shown in Table 1, the distance Difference observation error vector ε follows a mean value of zero and a covariance matrix of
Figure BDA00022095028800001714
A gaussian distribution of (a).
TABLE 1 Observation station 3-dimensional position coordinate (Unit: m)
Figure BDA00022095028800001715
Figure BDA0002209502880000181
The radiation source is first set to two cases: the 1 st type is a near field source with a position vector u ═ 220024002100]T(m); the 2 nd is a far-field source with a position vector u ═ 720084008900]T(m) of the reaction mixture. Figure 2 shows the variation of the root mean square error of the radiation source position estimate with the standard deviation sigma. The standard deviation σ is then set to two cases: 1 st is σ ═ 1; type 2 is σ ═ 2, and the radiation source position vector is set to u ═ 220022002200]T+[400 400 400]Tc (m). Fig. 3 shows the root mean square error of the radiation source position estimate as a function of the parameter c. As can be seen from fig. 2 and 3, the root mean square error of the positioning method disclosed in this patent for the position estimation of the radiation source can reach the cramer-perot limit (i.e., the lower theoretical limit). In addition, as can be seen from fig. 2 and 3, as the distance between the radiation source and the observation station increases, the positioning accuracy of the radiation source is gradually reduced, and the positioning accuracy of the radiation source for the near-field source is higher than that of the radiation source for the far-field source.
Comparing the positioning method disclosed in this patent with the existing multidimensional calibration method, it is noted that the existing multidimensional calibration method does not utilize quadratic equation constraints obeyed by the augmentation vector. The simulation parameters are the same as those in fig. 3, wherein σ is set to 1, and fig. 4 shows a variation curve of the root mean square error of the radiation source position estimation along with the parameter c; then, let σ be 2, fig. 5 shows the root mean square error of the radiation source position estimate as a function of the parameter c. As can be seen from fig. 4 and 5, since the existing multidimensional calibration method does not utilize the quadratic equation constraint obeyed by the augmented vector, the final positioning error is increased, and the negative effect on the positioning accuracy is related to the relative position between the radiation source and the observation station. As can be seen from fig. 4 and 5, the positioning method disclosed in this patent can indeed improve the positioning accuracy of the radiation source.
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 (9)

1. A time difference positioning method based on multi-dimensional calibration and polynomial root solving is characterized by comprising the following steps:
step 1: respectively counting the time difference of the radiation source signal reaching the mth sensor and the 1 st sensor, and multiplying the time difference by the signal propagation speed to obtain the observed distance difference
Figure FDA0002751656930000011
Wherein M is the number of sensors;
step 2: observed quantity by distance difference
Figure FDA0002751656930000012
Constructing a scalar product matrix, and constructing a multi-dimensional calibration pseudo linear equation through the scalar product matrix; the step 2 comprises the following steps:
step 2.1: using sensor position vector sm}1≤m≤MSum-distance difference observed quantity
Figure FDA0002751656930000013
Constructing an (M +1) × (M +1) order distance matrix D:
Figure FDA0002751656930000014
in the formula dmn＝||sm-sn||2(1≤m、n≤M)，dmnRepresents the distance between the m-th sensor and the n-th sensor;
step 2.2: calculating an (M +1) × (M +1) order scalar product matrix W using the distance matrix D:
Figure FDA0002751656930000015
in the formula
Figure FDA0002751656930000016
Is a transformation matrix, in which IM+1Represents an identity matrix of order (M +1) × (M +1), 1M+1Representing a full 1 vector of length M + 1;
step 2.3: by sensor position vector sm}1≤m≤MSum-distance difference observed quantity
Figure FDA0002751656930000017
Constructing a matrix G of (M +1) × 4 order, constructing a matrix G of (M +1) × 5 order
Figure FDA0002751656930000018
t1Representing the 1 st column vector in the matrix T, T2A matrix formed by vectors of 2 nd to 5 th columns in the matrix T is represented:
Figure FDA0002751656930000021
in the formula (I), the compound is shown in the specification,
Figure FDA0002751656930000022
in the formula
Figure FDA0002751656930000023
Is the distance difference observed for the 1 st sensor, which is constantly zero;
step 2.4: setting the iteration index k:equalto 0, setting an iteration threshold value delta, and calculating an iteration initial value through a scalar product matrix W and a matrix T
Figure FDA0002751656930000024
Step 2.5: observed quantity by matrix T and distance difference
Figure FDA0002751656930000025
Computing(M +1) × (M-1) order matrix
Figure FDA0002751656930000026
Calculating an (M +1) × (M-1) order matrix from the scalar product matrix W and the matrix G
Figure FDA0002751656930000027
Step 2.6: by means of a matrix
Figure FDA0002751656930000028
And
Figure FDA0002751656930000029
calculating an (M +1) × (M-1) order matrix
Figure FDA00027516569300000210
And to the matrix
Figure FDA00027516569300000211
Performing singular value decomposition
Figure FDA00027516569300000212
Wherein H(k)Is an orthogonal matrix of (M +1) × (M-1) order, V(k)Is an orthogonal matrix of (M-1) × (M-1) order(k)Is a diagonal matrix of (M-1) × (M-1) order, the diagonal elements of which are matrices
Figure FDA00027516569300000213
The singular value of (a);
step 2.7: after decomposition by singular values
Figure FDA00027516569300000214
Calculating a weighting matrix omega(k)＝Σ(k)V(k)TEV(k)Σ(k)TWherein E is a covariance matrix of range difference observation errors;
step 2.8: by scalar product matrix W, matrix T, H(k)And omega(k)Sequentially meterCalculation matrix phi(k)Sum vector
Figure FDA0002751656930000031
Figure FDA0002751656930000032
Step 2.9: through phi(k)Construction matrix
Figure FDA0002751656930000033
For matrix
Figure FDA0002751656930000034
And (3) carrying out characteristic value decomposition:
Figure FDA0002751656930000035
wherein, O3×1Representing a 3-dimensional all-zero column vector, O1×3Representing a 3-dimensional all-zero row vector, matrix P(k)Is a matrix of
Figure FDA0002751656930000036
Is determined by the feature vector of (a),
Figure FDA00027516569300000323
is a matrix
Figure FDA0002751656930000037
Characteristic value of (1)3Representing a 3 × 3 order identity matrix;
step 2.10: through phi(k)、
Figure FDA0002751656930000038
P(k)And the 1 st sensor position vector s1The vectors are calculated sequentially as follows
Figure FDA0002751656930000039
And
Figure FDA00027516569300000310
Figure FDA00027516569300000311
step 2.11: by passing
Figure FDA00027516569300000312
And
Figure FDA00027516569300000313
sequentially calculating scalars
Figure FDA00027516569300000314
Figure FDA00027516569300000315
And
Figure FDA00027516569300000316
step 2.12: by passing
Figure FDA00027516569300000317
s1And
Figure FDA00027516569300000318
sequentially calculating scalars
Figure FDA00027516569300000319
And
Figure FDA00027516569300000320
step 2.13: by passing
Figure FDA00027516569300000321
And
Figure FDA00027516569300000322
calculating a scalar quantity
Figure FDA0002751656930000041
Step 2.14: by passing
Figure FDA0002751656930000042
Constructing a multi-dimensional calibration pseudo linear equation;
and step 3: and solving a multi-dimensional calibration pseudo linear equation by using a Newton method, selecting a real root, eliminating a virtual root, and obtaining a position vector of the radiation source through the real root.
2. The time difference positioning method based on multi-dimensional calibration and polynomial root finding as claimed in claim 1, wherein said step 3 comprises:
step 3.1: solving a multi-dimensional calibration pseudo linear equation by using a Newton method, selecting a real root, and removing a pseudo root;
step 3.2: updating iteration index k: ═ k +1, and calculating iteration result
Figure FDA0002751656930000043
If it is
Figure FDA0002751656930000044
Stopping iteration, otherwise, turning to the step 2.5;
step 3.3: and determining the position vector of the radiation source by using the final iteration convergence result.
3. The time difference positioning method based on multi-dimensional calibration and polynomial root finding as claimed in claim 1, wherein said matrix
Figure FDA0002751656930000045
The calculation formula of (2) is as follows:
Figure FDA0002751656930000046
in the formula
Figure FDA0002751656930000047
Representing a vector
Figure FDA0002751656930000048
The m-th element of (1);
the matrix
Figure FDA0002751656930000049
The calculation formula of (2) is as follows:
Figure FDA0002751656930000051
in the formula
Figure FDA0002751656930000052
Representing a vector
Figure FDA0002751656930000053
The m-th element of (1);
Figure FDA0002751656930000054
representing a vector
Figure FDA0002751656930000055
The m-th element of (1); o is4×(M-1)Represents an all-zero matrix of 4 (M-1) th order, O3×(M-1)Represents an all-zero matrix of order 3 × (M-1), O1×(M-1)Represents an all-zero matrix of order 1 × (M-1), O(M+1)×1Denotes an (M +1) × 1 order all-zero matrix, O4×1Representing a 4 x 1 order all-zero matrix.
4. The time difference positioning method based on multi-dimensional calibration and polynomial root finding as claimed in claim 1, wherein said method comprises
Figure FDA0002751656930000056
And
Figure FDA0002751656930000057
the calculation formula of (2) is as follows:
Figure FDA0002751656930000058
wherein < >jRepresenting the jth component in the vector.
5. The time difference positioning method based on multi-dimensional calibration and polynomial root finding as claimed in claim 1, wherein said method comprises
Figure FDA0002751656930000059
And
Figure FDA00027516569300000510
the calculation formula of (2) is as follows:
Figure FDA0002751656930000061
6. the time difference positioning method based on multi-dimensional calibration and polynomial root finding as claimed in claim 1, wherein the scalar quantity is
Figure FDA0002751656930000062
The calculation formula of (2) is as follows:
Figure FDA0002751656930000063
7. the time difference positioning method based on multi-dimensional calibration and polynomial root finding as claimed in claim 2, wherein said step 3.1 comprises:
the multidimensional calibration pseudo linear equation is as follows:
Figure FDA0002751656930000071
let the root of the equation be
Figure FDA0002751656930000072
Selecting a real root by
Figure FDA0002751656930000073
Figure FDA0002751656930000074
In the formula
Figure FDA0002751656930000075
For the distance-difference observation vector,
Figure FDA0002751656930000076
indicating the use of the jth eigenvalue
Figure FDA0002751656930000077
The obtained position vector of the radiation source,
Figure FDA0002751656930000078
8. the time difference positioning method based on multi-dimensional calibration and polynomial root finding as claimed in claim 7, whereinCharacterized in that
Figure FDA0002751656930000079
The calculation formula of (2) is as follows:
Figure FDA00027516569300000710
9. the time difference positioning method based on multi-dimensional calibration and polynomial root finding as claimed in claim 8, wherein said step 3.3 comprises: using final iterative convergence results
Figure FDA00027516569300000711
The first 3 components of which serve as the position vector of the radiation source.
CN201910893455.3A 2019-09-20 2019-09-20 Time difference positioning method based on multidimensional calibration and polynomial root finding Active CN110673196B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910893455.3A CN110673196B (en) 2019-09-20 2019-09-20 Time difference positioning method based on multidimensional calibration and polynomial root finding

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910893455.3A CN110673196B (en) 2019-09-20 2019-09-20 Time difference positioning method based on multidimensional calibration and polynomial root finding

Publications (2)

Publication Number Publication Date
CN110673196A CN110673196A (en) 2020-01-10
CN110673196B true CN110673196B (en) 2021-01-22

Family

ID=69077053

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910893455.3A Active CN110673196B (en) 2019-09-20 2019-09-20 Time difference positioning method based on multidimensional calibration and polynomial root finding

Country Status (1)

Country Link
CN (1) CN110673196B (en)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111551897B (en) * 2020-04-25 2021-01-22 中国人民解放军战略支援部队信息工程大学 TDOA (time difference of arrival) positioning method based on weighted multidimensional scaling and polynomial root finding under sensor position error
CN111551896B (en) * 2020-04-25 2021-01-26 中国人民解放军战略支援部队信息工程大学 Weighted multidimensional scale TOA and FOA multi-source co-location method for inhibiting sensor position and speed prior errors
CN111551895B (en) * 2020-04-25 2021-01-26 中国人民解放军战略支援部队信息工程大学 Method for positioning TDOA and FDOA of motion source based on weighted multidimensional scale and Lagrange multiplier
CN112068099B (en) * 2020-07-31 2023-12-22 西安电子科技大学 Multi-radiation-source rapid positioning and speed measuring method and device based on error compensation

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102590786A (en) * 2012-02-03 2012-07-18 中国电子科技集团公司第三十八研究所 Multilateral positioning system based on distributed clock
EP2950119A1 (en) * 2014-05-26 2015-12-02 Ion Beam Applications S.A. System and method for verifying a particle beam
CN105353351A (en) * 2015-10-27 2016-02-24 杭州电子科技大学 Improved positioning method based on multi-beacon arrival time differences
CN105866735A (en) * 2016-04-06 2016-08-17 上海交通大学 Correction cost function time difference of arrival (TDOA) iteration positioning method based on multidimensional scaling (MDS) model
CN109239654A (en) * 2018-08-21 2019-01-18 中国人民解放军战略支援部队信息工程大学 Positioning using TDOA result method for correcting error neural network based
CN109298388A (en) * 2018-08-21 2019-02-01 中国人民解放军战略支援部队信息工程大学 Over-the-horizon target geographical coordinate direct method estimating based on azimuth information

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108414973B (en) * 2018-05-08 2020-03-13 中国人民解放军战略支援部队信息工程大学 Multi-target direct positioning method based on neural network calculation

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102590786A (en) * 2012-02-03 2012-07-18 中国电子科技集团公司第三十八研究所 Multilateral positioning system based on distributed clock
EP2950119A1 (en) * 2014-05-26 2015-12-02 Ion Beam Applications S.A. System and method for verifying a particle beam
CN105353351A (en) * 2015-10-27 2016-02-24 杭州电子科技大学 Improved positioning method based on multi-beacon arrival time differences
CN105866735A (en) * 2016-04-06 2016-08-17 上海交通大学 Correction cost function time difference of arrival (TDOA) iteration positioning method based on multidimensional scaling (MDS) model
CN109239654A (en) * 2018-08-21 2019-01-18 中国人民解放军战略支援部队信息工程大学 Positioning using TDOA result method for correcting error neural network based
CN109298388A (en) * 2018-08-21 2019-02-01 中国人民解放军战略支援部队信息工程大学 Over-the-horizon target geographical coordinate direct method estimating based on azimuth information

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
A Novel Weighted Multidimensional Scaling Analysis for Time-of-Arrival-Based Mobile Location;Wei H W;《IEEE Transactions on Signal Processing》;20081231;第3018页第1栏第1行-第3020页第2栏最后1行 *
基于最小二乘法的牛顿迭代信源定位算法;刘利军;《弹箭与制导学报》;20061231;第325页第1栏第2段-第328页第1栏第1段 *

Also Published As

Publication number Publication date
CN110673196A (en) 2020-01-10

Similar Documents

Publication Publication Date Title
CN110673196B (en) Time difference positioning method based on multidimensional calibration and polynomial root finding
CN111551895B (en) Method for positioning TDOA and FDOA of motion source based on weighted multidimensional scale and Lagrange multiplier
CN111551896B (en) Weighted multidimensional scale TOA and FOA multi-source co-location method for inhibiting sensor position and speed prior errors
CN108375752B (en) Amplitude-phase error single-radiation-source direction finding method based on all-angle search
CN111551897B (en) TDOA (time difference of arrival) positioning method based on weighted multidimensional scaling and polynomial root finding under sensor position error
CN104297718B (en) Interferometer array integrated correction method
CN105022026B (en) The two-dimentional angle estimation method of L-type array
CN106646453A (en) Doppler radar target tracking method based on predicted value measurement conversion
CN104793177B (en) Microphone array direction-finding method based on least square method
CN112379327A (en) Two-dimensional DOA estimation and cross coupling correction method based on rank loss estimation
CN109507635A (en) Utilize the array amplitude phase error evaluation method of two unknown orientation auxiliary sources
CN106353720A (en) Multi-station continuous positioning model based on TDOA/GROA (time different of arrival/gain ratio of arrival)
CN105629278A (en) GNSS pseudo-range single-point positioning-based high-precision mutual difference value median weighted positioning method
CN107797091B (en) Novel pure-direction target positioning method based on subspace
CN109212466B (en) Quantum dragonfly evolution mechanism-based broadband direction finding method
CN113835063A (en) Unmanned aerial vehicle array amplitude and phase error and signal DOA joint estimation method
CN111830465B (en) Two-dimensional Newton orthogonal matching pursuit compressed beam forming method
CN113835064B (en) Weighted multi-dimensional scale TDOA (time difference of arrival) positioning method for cooperative correction source observation information
CN112068099B (en) Multi-radiation-source rapid positioning and speed measuring method and device based on error compensation
CN112533284B (en) Near-far field unified positioning method based on arrival angle
Liang et al. Application of Taylor-Chan algorithm based on TDOA in sound source location
CN112540343A (en) Mobile target source positioning method based on mobile receiver cooperative analysis
CN104166795A (en) Combined sine-wave frequency estimation method based on multi-observation vector sparse representation
Li et al. An Underwater Source Localization Method Using Bearing Measurements
CN116883469B (en) Point cloud registration method based on EIV model description under plane feature constraint

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant