CN112986907B - Moving target positioning method under clock deviation and clock drift conditions - Google Patents
Moving target positioning method under clock deviation and clock drift conditions Download PDFInfo
- Publication number
- CN112986907B CN112986907B CN202110208858.7A CN202110208858A CN112986907B CN 112986907 B CN112986907 B CN 112986907B CN 202110208858 A CN202110208858 A CN 202110208858A CN 112986907 B CN112986907 B CN 112986907B
- Authority
- CN
- China
- Prior art keywords
- sensor node
- matrix
- representing
- node
- dimension
- 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
Links
- 238000000034 method Methods 0.000 title claims abstract description 121
- 238000005457 optimization Methods 0.000 claims abstract description 16
- 238000004364 calculation method Methods 0.000 claims abstract description 9
- 239000011159 matrix material Substances 0.000 claims description 64
- 239000013598 vector Substances 0.000 claims description 43
- 230000008859 change Effects 0.000 claims description 16
- 238000009826 distribution Methods 0.000 claims description 6
- 241000287196 Asthenes Species 0.000 claims description 4
- 238000006467 substitution reaction Methods 0.000 claims description 3
- 238000006243 chemical reaction Methods 0.000 claims description 2
- 238000005259 measurement Methods 0.000 description 13
- 238000007476 Maximum Likelihood Methods 0.000 description 5
- 238000000342 Monte Carlo simulation Methods 0.000 description 5
- 230000008569 process Effects 0.000 description 5
- 238000010586 diagram Methods 0.000 description 3
- 230000000694 effects Effects 0.000 description 3
- 230000004807 localization Effects 0.000 description 3
- 238000005070 sampling Methods 0.000 description 3
- 241001189642 Theroa Species 0.000 description 2
- 230000005540 biological transmission Effects 0.000 description 2
- 229960001948 caffeine Drugs 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 238000002474 experimental method Methods 0.000 description 2
- 230000002401 inhibitory effect Effects 0.000 description 2
- 238000004088 simulation Methods 0.000 description 2
- RYYVLZVUVIJVGH-UHFFFAOYSA-N trimethylxanthine Natural products CN1C(=O)N(C)C(=O)C2=C1N=CN2C RYYVLZVUVIJVGH-UHFFFAOYSA-N 0.000 description 2
- 238000004422 calculation algorithm Methods 0.000 description 1
- 239000002131 composite material Substances 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000000691 measurement method Methods 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 230000005855 radiation Effects 0.000 description 1
- 230000008707 rearrangement Effects 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
- 238000009827 uniform distribution Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S5/00—Position-fixing by co-ordinating two or more direction or position line determinations; Position-fixing by co-ordinating two or more distance determinations
- G01S5/02—Position-fixing by co-ordinating two or more direction or position line determinations; Position-fixing by co-ordinating two or more distance determinations using radio waves
- G01S5/0257—Hybrid positioning
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S5/00—Position-fixing by co-ordinating two or more direction or position line determinations; Position-fixing by co-ordinating two or more distance determinations
- G01S5/02—Position-fixing by co-ordinating two or more direction or position line determinations; Position-fixing by co-ordinating two or more distance determinations using radio waves
- G01S5/06—Position of source determined by co-ordinating a plurality of position lines defined by path-difference measurements
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02D—CLIMATE CHANGE MITIGATION TECHNOLOGIES IN INFORMATION AND COMMUNICATION TECHNOLOGIES [ICT], I.E. INFORMATION AND COMMUNICATION TECHNOLOGIES AIMING AT THE REDUCTION OF THEIR OWN ENERGY USE
- Y02D30/00—Reducing energy consumption in communication networks
- Y02D30/70—Reducing energy consumption in communication networks in wireless communication networks
Landscapes
- Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Mobile Radio Communication Systems (AREA)
Abstract
The invention discloses a moving target positioning method under the conditions of clock deviation and clock drift, which constructs a weighted least square WLS problem comprising four unknowns of the position and the speed of a source node, the clock deviation and the clock drift, then converts a non-convex weighted least square WLS problem into a convex semi-positive planning problem by adopting a semi-positive relaxation method, and then solves the convex semi-positive planning problem by utilizing an interior point method to obtain a global optimal solution; the method has the advantages that under the condition that the clock deviation and the clock drift exist and are unknown at the same time, the method realizes the joint estimation of the clock deviation and the clock drift while realizing the position and speed estimation of the source node from the perspective of a weighted least square estimation model estimation theory and a convex optimization theory, has high positioning precision and low calculation complexity, and is suitable for a complex deployment environment.
Description
Technical Field
The invention relates to a target positioning technology, in particular to a moving target positioning method under the conditions of clock deviation and clock drift.
Background
In recent years, with the rapid development of wireless sensor network technology, the target positioning technology is widely applied in the fields of military, industry and civilian use. In the existing target positioning technology, the time-of-arrival (TOA), frequency-of-arrival (FOA), time-difference-of-arrival (TDOA), frequency-difference-of-arrival (FDOA), angle-of-arrival (AOA), and Received Signal Strength (RSS) may be classified according to the signal measurement method. Generally, the TOA-based target location method and the TDOA-based target location method can achieve higher location accuracy. If relative motion exists between the positioning target and the sensor node, the positioning accuracy of the positioning target can be further improved on the basis of estimating the motion speed of the positioning target by introducing an observation signal of FOA or FDOA.
In a wireless sensor network, precise synchronization of clocks between a source node (i.e., positioning target) and a sensor node is a necessary condition for TOA-based target positioning, which is one of the challenges encountered in the TOA-based target positioning practical process. When using radio waves as a carrier for positioning, a clock skew of 10 nanoseconds will introduce an error of 3 meters (the propagation velocity of radio waves in air is 3 × 10)8m/s). In the document "Second-order con re-arrangement for TOA-based source localization with unknown start transmission time" (IEEE trans. veh. technol., vol.63, No.6, pp.2973-2977,2014), Wang et al propose a Second-order cone relaxation method that jointly estimates clock skew and source node position, and experimental results show that the method has better estimation performance than the existing min-max algorithm (MMA), however, the method is only for application scenarios where the source node and the sensor node are static. For a positioning scene with relative motion between a positioning target and a sensor node, a Passive positioning method based on importance sampling of TOA and FOA observation signals is disclosed in the literature "Passive positioning method using localization and imaging based on TOA and FOA observation signals" (Frontiers of Information Technology)&Electronic Engineering, vol.18, No.8, pp.1167-1179,2017), Liu et al propose a Monte Carlo Importance Sampling (Monte Carlo Importance Sampling) positioning method based on TOA and FOA observation signals, which can estimate the position and speed of the source node at the same time, and reach the Cramer-Rao lower bound (CRLB) in the low noise region, but the method is not limited to the above-mentioned methodThe method does not consider the problem of clock synchronization. In the published Chinese invention patent "weighted multidimensional scaling TOA and FOA multi-source co-location method for inhibiting sensor position and speed prior errors" (patent number: ZL 202010335968.5), Wang et al propose a weighted multidimensional scaling TOA and FOA multi-source co-location method for inhibiting sensor position and speed prior errors, the method utilizes a Linear Least Squares (LLS) estimation optimization model to obtain position vectors and estimated values of speed vectors of all radiation sources, although the method considers the problem of location of a plurality of source nodes, the method also has the defect of not considering clock synchronization. In view of the existing literature, when the TOA and FOA observation signals are used to estimate the position and the speed of the source node, no feasible solution exists for the scenario in which clock deviation and clock drift exist simultaneously.
Disclosure of Invention
The invention aims to solve the technical problem of providing a moving target positioning method under the conditions of Clock deviation and Clock Drift, which starts from the angles of Weighted Least Square (WLS) estimation model estimation theory and convex optimization theory under the condition that the Clock deviation (Clock Bias) and the Clock Drift (Clock Drift) exist and are unknown at the same time, realizes the joint estimation of the Clock deviation and the Clock Drift while realizing the position and speed estimation of a source node, has high positioning precision and low calculation complexity, and is suitable for complex deployment environments.
The technical scheme adopted by the invention for solving the technical problems is as follows: a method for positioning a moving target under the conditions of clock deviation and clock drift is characterized by comprising the following steps:
step 1: deploying a wireless sensor network in a planar space or a stereo space, wherein 1 source node with unknown position and speed for transmitting wireless signals, N sensor nodes with known position and speed for receiving wireless signals, and 1 central node for estimating the position and speed of the source node and clock offset and clock drift between the source node and the sensor nodes exist in the wireless sensor network, and the position of the source node is marked as xoThe velocity of the source node is recorded asvoThe positions of the N sensor nodes are respectively marked as s1,s2,…,si,…,sNThe velocities of the N sensor nodes are respectively recorded asWherein N is a positive integer, N is more than or equal to 5 if the wireless sensor network is deployed in a planar space, N is more than or equal to 6 if the wireless sensor network is deployed in a three-dimensional space, and s is greater than or equal to1Indicating the position of the 1 st sensor node, s2Indicating the position of the 2 nd sensor node, siIndicating the location of the ith sensor node, sNIndicating the location of the nth sensor node,representing the speed of the 1 st sensor node,representing the speed of the 2 nd sensor node,representing the speed of the ith sensor node,the speed of the Nth sensor node is represented, i is a positive integer, and i is more than or equal to 1 and less than or equal to N;
step 2: acquiring the arrival time of the wireless signal received by each sensor node according to the time of arrival (TOA) observation model, and recording the arrival time of the wireless signal received by the ith sensor node as taui,τiThe description is as follows:and acquiring the arrival frequency of the wireless signal received by each sensor node according to the arrival frequency FOA observation model, and recording the arrival frequency of the wireless signal received by the ith sensor node as fi,fiThe description is as follows:wherein the symbol "| | |" is a euclidean distance-solving symbol, c represents the propagation speed of the wireless signal,representing unknown clock deviations, f0Representing the carrier frequency, the superscript symbol "T" representing the transpose of a vector or matrix,representing unknown clock drift, Δ τiNoise, Δ f, measured for TOAiMeasuring noise for the FOA;
and step 3: each sensor node sends the arrival time and the arrival frequency of the wireless signals received by the sensor node to a central node;
and 4, step 4: the central node converts the arrival time of the wireless signal received by each sensor node into a distance observed quantity, and records the distance observed quantity of the wireless signal received by the ith sensor node as di,diThe description is as follows:similarly, the central node converts the arrival frequency of the wireless signal received by each sensor node into a distance change rate observed quantity, and records the distance change rate observed quantity of the wireless signal received by the ith sensor node as bi,biThe description is as follows: wherein,ni=Δτic, then let d denote the distance observation vector, let b denote the distance change rate observation vector, let d denoteoRepresenting the distance true vector, let boRepresenting the true vector of the rate of change of distance, d ═ d1,…,di,…,dN]T,b=[b1,…,bi,…,bN]T,Wherein d is1Representing a distance observation of a wireless signal received by the 1 st sensor node, dNRepresenting a distance observation of a wireless signal received by the Nth sensor node, b1Representing a range rate observation of a wireless signal received by the 1 st sensor node, bNA range rate observation representing a wireless signal received by the nth sensor node,andaccording toThe result of the calculation is that,andaccording toCalculating to obtain; and two noise vectors are introduced, which are respectively marked as n and m, and n is equal to [ n ═ n1,…,ni,…,nN]T,m=[m1,…,mi,…,mN]TWhere n and m are independent of each other, n follows a Gaussian distribution with a mean of 0, and a covariance matrix Qn=E[nnT]M obeys a Gaussian distribution with mean 0 and covariance matrix Qm=E[mmT],E[]To find the mathematical expectation, n1And nNAccording to ni=Δτic is calculated to obtain m1And mNAccording toCalculating to obtain;
and 5: to pairMake item shifting and conversion intoThen toSquaring both sides and ignoring the noise squared termTo obtainAnd will beSubstitution intoIn (1) obtainingRe-spreadThen, the term shift is carried out, and the second-order noise term n is ignoredimiTo obtain
Step 6: introducing an auxiliary vector yo,yoIs defined asThen will beRewritten in matrix form, described as: a. the1yo-h1≈B1n; and will beRewritten in matrix form, described as: a. the2yo-h2≈B2n+D2m; then combine A1yo-h1≈B1n and A2yo-h2≈B2n+D2m to obtain Ayo-h ≈ Bn + Dm; according to Ayo-h ≈ Bn + Dm constructs a weighted least squares WLS problem, described as:wherein,A1has a dimension of N × (2k +6), 01×kA column vector having dimensions of 1 × k and elements of 0 all is represented, where k is 2 if the wireless sensor network is deployed in a planar space, k is 3 if the wireless sensor network is deployed in a stereo space,h1has the dimension of N x 1, show thatThe diagonal matrix is composed as diagonal elements,A2has a dimension of N x (2k +6),h2has the dimension of N x 1, show thatThe diagonal matrix is composed as diagonal elements,show thatThe diagonal matrix is composed as diagonal elements,the dimension of A is 2N x (2k +6),the dimension of h is 2N x 1,the dimension of B is 2N x N,dimension of D is 2 NxN, 0N×NRepresenting a matrix of dimension N x N and elements all 0,expression of the result (Ay-h)TR-1The value of y when (Ay-h) is the smallest, y being yoCorresponding optimization variables, R represents a weight matrix, R ═ BQnBT+DQmDT=[B D]Q[B D]T,R-1Represents the inverse of R, [ B D ]]Is a matrix with dimension of 2 Nx 2N;
and 7: by digging yoThe implicit relationship between elements in the corresponding optimization variable y, further describes the weighted least squares WLS problem as:wherein "s.t." means "constrained to … …", y (1: k) denotes a vector of the 1 st to the kth elements in y, y (k +1:2k) denotes a vector of the (k +1) th to the 2 k-th elements in y, y (2k +1) denotes a (2k +1) th element in y, y (2k +2) denotes a (2k +2) th element in y, y (2k +3) denotes a (2k +3) th element in y, y (2k +4) denotes a (2k +4) th element in y, y (2k +5) denotes a (2k +5) th element in y, y (2k +6) denotes a (2k +6) th element in y;
and 8: the lead-in matrix Y is yyTWeighted Least Squares (WLS) problemIs relaxed toWherein,express to makeY and Y, tr () represents the trace of the matrix,to representIs semi-positive, rank () represents the rank of the matrix,andderived from the equivalent equation:
and step 9: problems after loosening of the semipositive fixationThe constraint of (1: k) | | y (1: k) | non-calculation2Y (2k +3) is rewritten as tr (Y (1: k,1: k)) ═ Y (2k +3), YT(1: k) Y (k +1:2k) ═ Y (2k +4) is rewritten as tr (Y (1: k, k +1:2k)) ═ Y (2k +4), Y2(2k +1) ═ Y (2k +5) rewrite to Y (2k +1) ═ Y (2k +5), Y (2k +1) Y (2k +2) ═ Y (2k +6) rewrite to Y (2k +1,2k +2) ═ Y (2k +6), Y (2k +1) Y (2k +6) ═ Y (2k +2) Y (2k +5) rewrite to Y (2k +1,2k +6) ═ Y (2k +2,2k + 5); wherein Y (1: k ) represents a matrix composed of elements of rows 1 to k and columns 1 to k in Y, Y (1: k, k +1:2k) represents a matrix composed of elements of rows 1 to k and columns (k +1) to 2k in Y, Y (2k +1) represents an element of rows 2k +1 and columns 2k +1 in Y, Y (2k +1,2k +2) represents an element of rows 2k +2 in Y, and Y (2k +1,2k +6) represents an element of rows 2k +1 and columns 2k +6 in Y;
step 10: discard theIn (1)And combining tr (Y (1: k,1: k)) ═ Y (2k +3), tr (Y (1: k, k +1:2k)) ═ Y (2k +4), Y (2k +1) ═ Y (2k +5), Y (2k +1,2k +2) ═ Y (2k +6), Y (2k +1,2k +6) ═ Y (2k +2,2k +5), we get a semi-positive programming problem, described as:
step 11: let R initial value be an identity matrix I with dimension of 2 Nx 2N2N(ii) a Then according toCalculating to obtain the value of F; then solving the semi-definite programming problem by using an interior point methodObtaining a primary global optimal solution of Y and a primary global optimal solution of Y, and correspondingly marking as Y*1And y*1(ii) a According to the definition of y and y*1To obtainThe respective preliminary estimate values are correspondingly recorded asWherein, y*1(2k +1) represents y*1The (2k +1) th element, y*1(2k +2) represents y*1The (2k +2) th element in (a);
step 12: according toAnd in combination with R ═ BQnBT+DQmDT=[B D]Q[B D]TUpdating the value of R and further updating the value of F; then, solving the semi-definite programming problem by utilizing the interior point method againObtaining the final global optimal solution of Y and the final global optimal solution of Y, and correspondingly recording as Y*2And y*2(ii) a According to the definition of y and y*2To obtain xo、vo、The respective final estimated values, corresponding to x*2、v*2、x*2=y*2(1:k),v*2=y*2(k+1:2k), Then according toAndandto obtainAndis expressed asAnd wherein, y*2(1: k) represents y*2The vector y consisting of the 1 st element to the k th element of (A)*2(k +1:2k) represents y*2The (k +1) th element to the 2 k-th element of (a), y*2(2k +1) represents y*2The (2k +1) th element, y*2(2k +2) represents y*2The (2k +2) th element in (b).
Compared with the prior art, the invention has the advantages that:
1) the method converts the non-convex weighted least square WLS problem into the convex semi-positive definite programming problem by using a semi-positive definite relaxation technology, and solves the convex semi-positive definite programming problem by using an interior point method, so that a global optimal solution can be obtained, and the defects that the traditional maximum likelihood ML estimation method falls into a local optimal point and the positioning precision is low are overcome.
2) According to the method, the sequence and the combination of variables of the composite vector containing the position, the speed, the clock deviation and the clock drift of the source node are designed from multiple angles of simplicity, smoothness, balance and the like, so that the relaxed semi-definite programming problem is as tight as possible, and subsequent operations caused by the fact that the semi-definite programming problem is not tight enough are omitted, such as complexity improvement caused by the methods of Gaussian randomization, local search and the like.
3) According to the method, the position, the speed, the clock deviation and the clock drift of the source node can be recovered directly through the received wireless signal without taking the clock deviation and the clock drift as known conditions, and the calculation complexity is low.
4) Under the condition that the clock deviation and the clock drift exist and are unknown at the same time, the method can realize the joint estimation of the clock deviation and the clock drift while realizing the position and the speed estimation of the source node, and is suitable for complex deployment environments.
Drawings
FIG. 1 is a block diagram of an overall implementation of the method of the present invention;
figure 2 shows that the number of sensor nodes is 5,from 10-2m2Change to 103m2In the time, the RMSE performance of the position estimation of the method (SDP for short), the ML method and the Raw WLS method of the invention is dependent onA profile of change;
figure 3 shows that the number of sensor nodes is 5,from 10-2m2Change to 103m2In the method of the inventionRMSE Performance of velocity estimation by (SDP for short), ML method, Raw WLS methodA profile of change;
figure 4 shows that the number of sensor nodes is 5,from 10-2m2Change to 103m2In time, the RMSE performance of the clock bias estimation of the method (SDP for short) and the ML method of the invention is followedA profile of change;
figure 5 shows that the number of sensor nodes is 5,from 10-2m2Change to 103m2In time, the RMSE performance of the clock drift estimation of the method (SDP for short) and the ML method of the inventionA profile of change;
FIG. 6 is a drawing showingFixed at 10m2When the number of the sensors is increased from 4 to 9, the RMSE performance of the position estimation of the method (abbreviated as SDP), the ML method and the Raw WLS method changes along with the number of the sensor nodes;
FIG. 7 is a drawing showingFixed at 10m2When the number of the sensors is increased from 4 to 9, the RMSE performance of the speed estimation of the method (called SDP for short), the ML method and the Raw WLS method changes along with the number of the sensor nodes;
Detailed Description
The invention is described in further detail below with reference to the accompanying examples.
The general implementation block diagram of the moving target positioning method under the clock deviation and clock drift conditions provided by the invention is shown in fig. 1, and the method comprises the following steps:
step 1: deploying a wireless sensor network in a planar space or a stereo space, wherein 1 source node with unknown position and speed for transmitting wireless signals, N sensor nodes with known position and speed for receiving wireless signals, and 1 central node for estimating the position and speed of the source node and clock offset and clock drift between the source node and the sensor nodes exist in the wireless sensor network, and the position of the source node is marked as xoLet the velocity of the source node be voThe positions of the N sensor nodes are respectively marked as s1,s2,…,si,…,sNThe velocities of the N sensor nodes are respectively recorded asWherein N is a positive integer, N is more than or equal to 5 if the wireless sensor network is deployed in a planar space, N is more than or equal to 6 if the wireless sensor network is deployed in a three-dimensional space, and s is greater than or equal to1Indicating the position of the 1 st sensor node, s2Indicating the position of the 2 nd sensor node, siIndicating the location of the ith sensor node, sNIndicating the location of the nth sensor node,representing the speed of the 1 st sensor node,representing the speed of the 2 nd sensor node,representing the speed of the ith sensor node,and the speed of the Nth sensor node is shown, i is a positive integer, and i is more than or equal to 1 and less than or equal to N.
Step 2: acquiring the arrival time of the wireless signal received by each sensor node according to the time of arrival (TOA) observation model, and recording the arrival time of the wireless signal received by the ith sensor node as taui,τiThe description is as follows:and acquiring the arrival frequency of the wireless signal received by each sensor node according to the arrival frequency FOA observation model, and recording the arrival frequency of the wireless signal received by the ith sensor node as fi,fiThe description is as follows:wherein the symbol "| | |" is a euclidean distance-solving symbol, c represents the propagation speed of the wireless signal,representing unknown clock deviations, f0Representing the carrier frequency, the superscript symbol "T" representing the transpose of a vector or matrix,representing unknown clock drift, Δ τiNoise, Δ f, measured for TOAiNoise was measured for FOA.
And step 3: each sensor node transmits the time of arrival and frequency of arrival of the wireless signals it receives to the central node.
And 4, step 4: the central node transmits each signalThe arrival time of the wireless signal received by the sensor node is converted into a distance observation amount roa (range of arrival), and the distance observation amount of the wireless signal received by the ith sensor node is recorded as di,diThe description is as follows:similarly, the central node converts the arrival frequency of the wireless signal received by each sensor node into a distance Rate observed quantity rroa (range Rate of arrival), and records the distance Rate observed quantity of the wireless signal received by the i-th sensor node as bi,biThe description is as follows: wherein,ni=Δτic, then let d denote the distance observation vector, let b denote the distance change rate observation vector, let d denoteoRepresenting the distance true vector, let boRepresenting the true vector of the rate of change of distance, d ═ d1,…,di,…,dN]T,b=[b1,…,bi,…,bN]T,Wherein d is1Representing a distance observation of a wireless signal received by the 1 st sensor node, dNRepresenting a distance observation of a wireless signal received by the Nth sensor node, b1Representing a range rate observation of a wireless signal received by the 1 st sensor node, bNA range rate observation representing a wireless signal received by the nth sensor node,andaccording toThe calculation results in that,andaccording toCalculating to obtain; and two noise vectors are introduced, which are respectively marked as n and m, and n is equal to [ n ═ n1,…,ni,…,nN]T,m=[m1,…,mi,…,mN]TWhere n and m are independent of each other, n follows a Gaussian distribution with a mean of 0, and a covariance matrix Qn=E[nnT]M obeys a Gaussian distribution with mean 0 and covariance matrix Qm=E[mmT],E[]To find the mathematical expectation, n1And nNAccording to ni=Δτic is calculated to obtain m1And mNAccording toAnd (4) calculating.
If the maximum likelihood estimation (ML) model is then used to solve, then the following process is: forming all unknowns into a vector by thetaoIt is shown that,wherein, thetaoHas a dimension of (2k +2), k if the wireless sensor network is deployed in a planar space2, if the wireless sensor network is deployed in a stereo space, k is 3; d and b form a vector, which is expressed by p, and p is ═ dT,bT]T(ii) a And d isoAnd boForm a vector, with poIs represented by the formula po=[doT,boT]T(ii) a Finally, a maximum likelihood estimation (ML) model is constructed, which is described as:wherein theta is thetaoCorresponding optimization variable, θ ═ xT,vT,d0,b0]TX is xoCorresponding optimization variable, v is voCorresponding optimization variable, d0Is composed ofCorresponding optimization variables, b0Is composed ofThe corresponding optimization variables are set to be the same as,p (theta) isoθ inoIs substituted by θ, Q represents a weight matrix, Q ═ blkdiag (Q)n,Qm) Blkdiag () denotes performing a block diagonal matrix operation, Q-1The inverse matrix of Q is represented by,express to makeThe value of theta at which the value of (b) is minimum,expression of equation (p-p (theta))TQ-1(p-p (θ)) is the value of θ at the minimum.
And 5: since the solution to the maximum likelihood estimation (ML) model is easy to fall into a local optimum point, the positioning accuracy is influenced,thus the invention is rightMake item shifting and transformation intoThen toSquaring both sides and ignoring the noise squared termTo obtainAnd will beSubstitution intoIn (1) obtainingRe-spreadThen, the term shift is carried out, and the second-order noise term n is ignoredimiTo obtain
And 6: introducing an auxiliary vector yo,yoIs defined asThen will beRewritten in matrix form, described as: a. the1yo-h1≈B1n; and will beRewritten in matrix form, described as: a. the2yo-h2≈B2n+D2m; then combine A1yo-h1≈B1n and A2yo-h2≈B2n+D2m to obtain Ayo-h ≈ Bn + Dm; then according to Ayo-h ≈ Bn + Dm constructs a weighted least squares WLS problem, described as:wherein,A1has a dimension of N × (2k +6), 01×kA column vector having dimensions of 1 × k and elements of 0 all is represented, where k is 2 if the wireless sensor network is deployed in a planar space, k is 3 if the wireless sensor network is deployed in a stereo space,h1has the dimension of (a) of N x 1, show thatA diagonal matrix is composed as diagonal elements,A2has a dimension of N x (2k +6),h2has the dimension of N x 1, show thatThe diagonal matrix is composed as diagonal elements,show thatThe diagonal matrix is composed as diagonal elements,the dimension of A is 2N x (2k +6),the dimension of h is 2N x 1,the dimension of B is 2N x N,dimension of D is 2 NxN, 0N×NRepresenting a matrix of dimension N x N and elements all 0,expression of the result (Ay-h)TR-1The value of y when (Ay-h) is the smallest, y being yoCorresponding optimization variables, R represents a weight matrix, R ═ BQnBT+DQmDT=[B D]Q[B D]T,R-1Represents the inverse of R, [ B D ]]A matrix with dimensions 2N × 2N.
And 7: since the weighted least squares WLS problem is still a non-convex optimization problem that is difficult to solve, the present invention exploits yoThe implicit relationship between the elements in the corresponding optimization variable y, re-describes the weighted least squares WLS problem as:where "s.t." means "constrained to … …", y (1: k) denotes a vector composed of the 1 st to the kth elements in y, y (k +1:2k) denotes a vector composed of the (k +1) th to the 2 k-th elements in y, y (2k +1) denotes the (2k +1) th element in y, y (2k +2) denotes the (2k +2) th element in y, y (2k +3) denotes the (2k +3) th element in y, y (2k +4) denotes the (2k +4) th element in y, y (2k +5) denotes the (2k +5) th element in y, and y (2k +6) denotes the (2k +6) th element in y.
And 8: the lead-in matrix Y is yyTWeighted least squares WLS problemIs relaxed toWherein,express to makeY and Y, tr () represents the trace of the matrix,to representIs semi-positive, rank () represents the rank of the matrix,andderived from the equivalent equation:
and step 9: problems after loosening of the semipositive fixationThe constraint of (1: k) | | y (1: k) | non-calculation2Y (2k +3) is rewritten as tr (Y (1: k,1: k)) ═ Y (2k +3), YT(1: k) Y (k +1:2k) ═ Y (2k +4) is rewritten as tr (Y (1: k, k +1:2k)) ═ Y (2k +4), Y2(2k +1) ═ Y (2k +5) rewrite to Y (2k +1) ═ Y (2k +5), Y (2k +1) Y (2k +2) ═ Y (2k +6) rewrite to Y (2k +1,2k +2) ═ Y (2k +6), Y (2k +1) Y (2k +6) ═ Y (2k +2) Y (2k +5) rewrite to Y (2k +1,2k +6) ═ Y (2k +2,2k + 5); wherein Y (1: k ) represents a matrix composed of elements of rows 1 to k and columns 1 to k in Y, Y (1: k, k +1:2k) represents a matrix composed of elements of rows 1 to k and columns (k +1) to 2k in Y, Y (2k +1) represents elements of rows 2k +1 and columns 2k +1 in Y, Y (2k +1,2k +2) represents elements of rows 2k +1 and columns 2k +2 in Y, and Y (2k +1,2k +6) represents elements of rows 2k +1 and columns 2k +6 in Y.
Step 10: discard theIn (1)And combining tr (Y (1: k,1: k)) ═ Y (2k +3), tr (Y (1: k, k +1:2k)) ═ Y (2k +4), Y (2k +1) ═ Y (2k +5), Y (2k +1,2k +2) ═ Y (2k +6), Y (2k +1,2k +6) ═ Y (2k +2,2k +5), we get a semi-positive programming problem, described as:
step 11: let R initial value be an identity matrix I with dimension of 2 Nx 2N2N(ii) a Then according toCalculating to obtain the value of F; then solving the semi-definite programming problem by using an interior point methodObtaining a primary global optimal solution of Y and a primary global optimal solution of Y, and correspondingly marking as Y*1And y*1(ii) a According to the definition of y and y*1To obtainThe respective preliminary estimate values are correspondingly recorded asWherein, y*1(2k +1) represents y*1The (2k +1) th element, y*1(2k +2) represents y*1The (2k +2) th element in (b).
Step 12: according toAnd in combination with R ═ BQnBT+DQmDT=[B D]Q[B D]TUpdating the value of R and further updating the value of F; then, solving the semi-definite programming problem by utilizing the interior point method againObtaining the final global optimal solution of Y and the final global optimal solution of Y, and correspondingly recording as Y*2And y*2(ii) a According to the definition of y and y*2To obtain xo、vo、The respective final estimated values, corresponding to x*2、v*2、x*2=y*2(1:k),v*2=y*2(k+1:2k), Then according toAndandto obtainAndis expressed asAnd wherein, y*2(1: k) represents y*2The vector y consisting of the 1 st element to the k th element of (A)*2(k +1:2k) represents y*2The (k +1) th element to the 2 k-th element of (c), y*2(2k +1) represents y*2The (2k +1) th element, y*2(2k +2) represents y*2The (2k +2) th element in (b).
The effectiveness and feasibility of the method can be verified through simulation experiments.
The method is characterized in that 9 (N ═ 9) sensor nodes are deployed on a two-dimensional plane space, the positions and the speeds of the sensor nodes are shown in table 1, and the two-dimensional Cartesian coordinates are (X, Y).
TABLE 1 position and velocity of 9 sensor nodes in a two-dimensional scene
The position of the source node is randomly generated in a circular ring area, the center of the circular ring is at the origin, the outer radius and the inner radius are respectively 1000 and 200, and the unit is meter; the velocity of the source node is randomly generated within a circle centered at the origin with a radius of 20 m/s. Therefore, the source node may fall inside or outside the Convex Hull (Convex Hull) formed by the sensor nodes. Caused by clock skew Randomly generated (in meters) from a uniformly distributed U (0,150), caused by clock drift Randomly generated (in meters per second) from a uniform distribution U (0, 15).
A range observation and a range rate observation of the wireless signal received by each sensor node are generated with step 4. The covariance matrices of the measured noise are respectivelyAnd represents the variance of the noise of the TOA measurement,representing FOA measurement noiseVariance of sound, INRepresenting an identity matrix of dimension N x N. In general, the FOA measurement noise is less than the TOA measurement noise, and therefore, can be usedThe relationship between the two is described, and the value range of the FOA measurement noise coefficient eta is generally [0.001, 1%]In the simulation of the present invention, the FOA measurement noise coefficient η used in the first 2 scenes (i.e., scene 1 and scene 2) is 0.1.
The performance of the target positioning method can be expressed by mean square error (RMSE), which is defined as:wherein M is1Indicates the number of source nodes, M2Representing the times of Monte Carlo simulation of each source node, q is more than or equal to 1 and less than or equal to M1,1≤j≤M2,xqRepresenting the true position value of the qth source node in the jth monte carlo simulation,the position estimate of the qth source node in the jth monte carlo simulation is shown. The estimated performance RMSE of the speed estimation and clock bias and clock drift of the source node can be defined in the same way, and M is taken here150(50 source nodes), M2500 (500 monte carlo simulations per source node).
To embody the RMSE performance of the method of the invention, an ML method (maximum likelihood estimation method) was introduced for comparison. The ML method is a commonly used method in parameter estimation, and in the specific implementation process, the method needs an initial value, then iterative computation is carried out, and the solution of the method is taken as the initial value to be substituted into the ML method for computation in the experiment.
In order to embody the importance of processing clock deviation and clock drift, a Raw WLS method is also introduced for comparison. The Raw WLS method treats clock bias and clock drift as noise, directly ignores its influence, and estimates only the position and velocity of the source node. The Raw WLS method is obtained by modifying the method of the invention, and comprises the following specific steps:
in the ROA observation vector in step 4, if the clock deviation is ignoredThen become di=||xo-si||+niI 1, …, N, in the RROA observation vector, if clock drift is neglectedThen become intoReferring to step 5, the ROA and RROA are expanded and the second order noise terms are ignored to obtainAndthereafter, referring to the method of step 6, an auxiliary vector is introducedIs defined asThen will beRewritten in matrix form, described as:and will beRewritten in matrix form, described as:then are connected togetherAndto obtainThen according toConstructing a weighted least squares WLS problem, described as:whereinA3Has a dimension of N × (2k +2), 01×kA column vector having dimensions of 1 × k and elements of 0 all is represented, where k is 2 if the wireless sensor network is deployed in a planar space, k is 3 if the wireless sensor network is deployed in a stereo space,h3has dimension of Nx 1, B3=diag(-2d1,…,-2dN),B3=diag(-2d1,…,-2dN) Represents that 2d is1,…,-2dNThe diagonal matrix is composed as diagonal elements,A4has a dimension of N x (2k +2),h4has dimension of Nx 1, B4=diag(-b1,…,-bN),diag(-b1,…,-bN) Is represented by1,…,-bNForming a diagonal matrix as diagonal elements, D4=diag(-d1,…,-dN),diag(-d1,…,-dN) Is represented by1,…,-dNThe diagonal matrix is composed as diagonal elements,ARWLShas a dimension of 2N x (2k +2),hRWLShas a dimension of 2N x 1,BRWLShas a dimension of 2N x N,DRWLShas a dimension of 2 NxN, 0N×NRepresenting a matrix of dimension N x N and elements all 0,express to makeY is the smallest value ofRWLSValue of (a), yRWLSIs composed ofCorresponding optimization variable, RRWLSA matrix of weights is represented by a matrix of weights, represents RRWLSInverse of (B)RWLS DRWLS]A matrix with dimensions 2N × 2N.
FIG. 2 shows the RMSE performance of the location estimation of the present invention method (SDP for short), ML method, and Raw WLS methodThe curves of the variation, and FIG. 3 shows the RMSE performance of the velocity estimation of the method of the present invention (SDP for short), ML method, and Raw WLS method as a function ofThe curve of the change. As can be seen from FIGS. 2 and 3, the RMSE performance of the method of the present invention is very close to that of the ML method, and the method of the present invention can reach the CRLB (Cramer-Rao Lower Bound) boundary in the range of small noise and medium noise. Because the Raw WLS method does not take clock bias and clock drift into account, the RMSE performance of the Raw WLS method is controlled by clock bias and clock drift in the small and medium noise intervals, much worse than the method of the present invention. In case of very loud noiseThe effects of clock bias and clock drift are less than the effects of measurement noise and the RMSE performance of the Raw WLS method is superior to the inventive method and ML method because the Raw WLS method estimates fewer parameters.
FIG. 4 shows the RMSE performance of the clock bias estimation of the present invention method (SDP for short), ML methodThe curves of the variation, FIG. 5 shows the RMSE performance of the clock drift estimation of the method of the invention (SDP for short), ML methodThe curve of the change. As can be seen from FIGS. 4 and 5, the method of the present invention can reach the CRLB boundary in the small noise and medium noise interval, and shows better RMSE performance. To assess the degree of relaxation of the process of the invention, Table 3 gives each of theCorresponding to noise conditions, M1=50、M2The Rank of the final global optimal solution for Y of the inventive method in 25000 monte carlo simulation experiments performed on 500, i.e. 50 source nodes, is the number of 1 (Rank-1 number). As can be seen from Table 2, in most cases the process of the invention is tight, even when used in the fieldIs 30m2At that time, the ratio of the number of ranks of the final global optimal solution of Y being 1 (Rank-1 number) exceeds 99.74%.
TABLE 2 Rank-1 number of SDP solutions in two-dimensional scenarios
Fig. 6 shows the curves of RMSE performance of position estimation with the number of sensor nodes in the methods of the present invention (SDP for short), ML method, and Raw WLS method, and fig. 7 shows the curves of RMSE performance of speed estimation with the number of sensor nodes in the methods of the present invention (SDP for short), ML method, and Raw WLS method. As can be seen from fig. 6 and 7, the method of the present invention can reach the CRLB boundary except at the point where the number of sensor nodes is 4. In addition, under this measurement noise condition, both the position estimate and the velocity estimate, the RMSE performance of the Raw WLS method is worse than that of the present invention method and the ML method. This represents an importance for dealing with both clock skew and clock drift.
Scene 3: in order to compare the RMSE performance differences between positioning using TOA and FOA observation signals simultaneously and positioning using TOA observation signals only, the method of the present invention was compared with the method of the document "Second-order con relaxation positioning based on TOA at unknown initial transmission time" (g.wang, s.cai, y.li, and m.jin, IEEE trans.veh.technol., vol.63, No.6, pp.2973-2977,2014.) (abbreviated TOA-SOCP).
RMSE performance in the two-dimensional case was tested here, using the first 5 sensor nodes in table 1; measuring noiseFixed at 10m2。
FIG. 8 shows the FOA measurement noise figure η from 10-3Change to 10-1In time, the RMSE performance of the position estimation of the method and the TOA-SOCP method is changed along with the FOA measurement noise coefficient. As can be seen from fig. 8, since the variation of the noise figure of the FOA measurement does not affect the TOA observation signal, the localization performance of the TOA-SOCP method is unchanged. In addition, use TOA and FOA observation signal to fix a position simultaneously, compare with only using TOA observation signal to fix a position, can show and promote the positioning performance, and the noise of FOA observation signal is less, and is better to the promotion effect of whole positioning performance. Meanwhile, compared with the TOA-SOCP method, the method is closer to the CRLB, and the method has better performance. In FIG. 8, TOA-CRLB indicates the CRLB boundary when the signal is observed using only TOA, and TOA/FOA-CRLB indicates the CRLB boundary when both the signal is observed using TOA and FOA.
Claims (1)
1. A method for positioning a moving target under the conditions of clock deviation and clock drift is characterized by comprising the following steps:
step 1: deploying a wireless sensor network in a planar space or a stereo space, wherein 1 source node with unknown position and speed for transmitting wireless signals, N sensor nodes with known position and speed for receiving wireless signals, and 1 source node for estimating the position and the speed of the wireless sensor network existThe position and the speed of the point, the clock deviation and the clock drift between the source node and the sensor node are taken as the central node, and the position of the source node is recorded as xoLet the velocity of the source node be voThe positions of the N sensor nodes are respectively marked as s1,s2,…,si,…,sNThe velocities of the N sensor nodes are respectively recorded asWherein N is a positive integer, N is more than or equal to 5 if the wireless sensor network is deployed in a planar space, N is more than or equal to 6 if the wireless sensor network is deployed in a three-dimensional space, and s is greater than or equal to1Indicating the position of the 1 st sensor node, s2Indicating the position of the 2 nd sensor node, siIndicating the location of the ith sensor node, sNIndicating the location of the nth sensor node,representing the speed of the 1 st sensor node,representing the speed of the 2 nd sensor node,representing the speed of the ith sensor node,the speed of the Nth sensor node is represented, i is a positive integer, and i is more than or equal to 1 and less than or equal to N;
step 2: acquiring the arrival time of the wireless signal received by each sensor node according to the time of arrival (TOA) observation model, and recording the arrival time of the wireless signal received by the ith sensor node as taui,τiThe description is as follows:and according to the frequency of arrival FOA observationThe model is used for acquiring the arrival frequency of the wireless signal received by each sensor node and recording the arrival frequency of the wireless signal received by the ith sensor node as fi,fiThe description is as follows:wherein the symbol "| | |" is a euclidean distance-solving symbol, c represents the propagation speed of the wireless signal,representing unknown clock deviations, f0Representing the carrier frequency, the superscript symbol "T" representing the transpose of a vector or matrix,representing unknown clock drift, Δ τiNoise, Δ f, measured for TOAiMeasuring noise for the FOA;
and step 3: each sensor node sends the arrival time and the arrival frequency of the wireless signals received by the sensor node to a central node;
and 4, step 4: the central node converts the arrival time of the wireless signal received by each sensor node into a distance observed quantity, and records the distance observed quantity of the wireless signal received by the ith sensor node as di,diThe description is as follows:similarly, the central node converts the arrival frequency of the wireless signal received by each sensor node into a distance change rate observed quantity, and records the distance change rate observed quantity of the wireless signal received by the ith sensor node as bi,biThe description is as follows: wherein,ni=Δτic, then let d denote the distance observation vector, let b denote the distance change rate observation vector, let d denoteoRepresenting the distance true vector, let boRepresenting the true vector of the rate of change of distance, d ═ d1,…,di,…,dN]T,b=[b1,…,bi,…,bN]T,Wherein d is1Representing a distance observation of a wireless signal received by the 1 st sensor node, dNRepresenting a distance observation of a wireless signal received by the Nth sensor node, b1Representing a range rate observation of a wireless signal received by the 1 st sensor node, bNA range rate observation representing a wireless signal received by the nth sensor node,andaccording toThe calculation results in that,andaccording toCalculating to obtain; and two noise vectors are introduced, which are respectively marked as n and m, and n is equal to [ n ═ n1,…,ni,…,nN]T,m=[m1,…,mi,…,mN]TWhere n and m are independent of each other, n follows a Gaussian distribution with a mean of 0, and a covariance matrix Qn=E[nnT]M obeys a Gaussian distribution with mean 0 and covariance matrix Qm=E[mmT],E[]To find the mathematical expectation, n1And nNAccording to ni=Δτic is calculated to obtain m1And mNAccording toCalculating to obtain;
and 5: to pairMake item shifting and conversion intoThen toSquaring both sides and ignoring the noise squared termTo obtainAnd will beSubstitution intoIn (1) obtainingRe-spreadThen, the term shift is carried out, and the second-order noise term n is ignoredimiTo obtain
Step 6: introducing an auxiliary vector yo,yoIs defined asThen will beRewritten in matrix form, described as: a. the1yo-h1≈B1n; and will beRewritten in matrix form, described as: a. the2yo-h2≈B2n+D2m; then combine A1yo-h1≈B1n and A2yo-h2≈B2n+D2m to obtain Ayo-h ≈ Bn + Dm; according to Ayo-h ≈ Bn + Dm constructs a weighted least squares WLS problem, described as:wherein,A1has a dimension of N × (2k +6), 01×kA column vector having dimensions of 1 × k and all elements of 0 is represented, where k is 2 if the wireless sensor network is deployed in a planar space and k is 2 if the wireless sensor network is deployed in a stereoscopic spacek=3,h1Has the dimension of (a) of N x 1, show thatThe diagonal matrix is composed as diagonal elements,A2has a dimension of N x (2k +6),h2has the dimension of N x 1, show thatThe diagonal matrix is composed as diagonal elements, show thatThe diagonal matrix is composed as diagonal elements,the dimension of A is 2N x (2k +6),the dimension of h is 2N x 1,the dimension of B is 2N x N,dimension of D is 2 NxN, 0N×NRepresenting a matrix of dimension N x N and elements all 0,expression of the result (Ay-h)TR-1The value of y when the value of (Ay-h) is minimum, y beingoCorresponding optimization variables, R represents a weight matrix, R ═ BQnBT+DQmDT=[B D]Q[B D]T,R-1Represents the inverse of R, [ B D ]]Is a matrix with dimension of 2 Nx 2N;
and 7: by digging yoThe implicit relationship between elements in the corresponding optimization variable y, further describes the weighted least squares WLS problem as:wherein "s.t." means "constrained to … …", y (1: k) denotes a vector of the 1 st to the kth elements in y, y (k +1:2k) denotes a vector of the (k +1) th to the 2 k-th elements in y, y (2k +1) denotes a (2k +1) th element in y, y (2k +2) denotes a (2k +2) th element in y, y (2k +3) denotes a (2k +3) th element in y, y (2k +4) denotes a (2k +4) th element in y, y (2k +5) denotes a (2k +5) th element in y, y (2k +6) denotes a (2k +6) th element in y;
and 8: the lead-in matrix Y is yyTWeighted least squares WLS problemIs relaxed toWherein,express to makeY and Y, tr () represents the trace of the matrix, to representIs semi-positive, rank () represents the rank of the matrix,andderived from the equivalent equation:
and step 9: the problem after semi-positive relaxation is describedThe constraint of (1: k) | | y (1: k) | non-calculation2Y (2k +3) is rewritten as tr (Y (1: k,1: k)) ═ Y (2k +3), YT(1: k) Y (k +1:2k) ═ Y (2k +4) is rewritten as tr (Y (1: k, k +1:2k)) ═ Y (2k +4), Y2(2k +1) ═ Y (2k +5) is rewritten to Y (2k +1) ═ Y (2k +5), Y (2k +5)+1) Y (2k +2) ═ Y (2k +6) rewrite to Y (2k +1,2k +2) ═ Y (2k +6), Y (2k +1) Y (2k +6) ═ Y (2k +2) Y (2k +5) rewrite to Y (2k +1,2k +6) ═ Y (2k +2,2k + 5); wherein Y (1: k ) represents a matrix composed of elements of rows 1 to k and columns 1 to k in Y, Y (1: k, k +1:2k) represents a matrix composed of elements of rows 1 to k and columns (k +1) to 2k in Y, Y (2k +1) represents an element of rows 2k +1 and columns 2k +1 in Y, Y (2k +1,2k +2) represents an element of rows 2k +2 in Y, and Y (2k +1,2k +6) represents an element of rows 2k +1 and columns 2k +6 in Y;
step 10: discard theInAnd combining tr (Y (1: k,1: k)) ═ Y (2k +3), tr (Y (1: k, k +1:2k)) ═ Y (2k +4), Y (2k +1) ═ Y (2k +5), Y (2k +1,2k +2) ═ Y (2k +6), Y (2k +1,2k +6) ═ Y (2k +2,2k +5), we get a semi-positive programming problem, described as:
step 11: let R initial value be an identity matrix I with dimension of 2 Nx 2N2N(ii) a Then according toCalculating to obtain the value of F; then solving the semi-definite programming problem by using an interior point methodObtaining a primary global optimal solution of Y and a primary global optimal solution of Y, and correspondingly marking as Y*1And y*1(ii) a According to the definition of y and y*1To obtainPreliminary estimate of each, pairShould be noted asWherein, y*1(2k +1) represents y*1The (2k +1) th element, y*1(2k +2) represents y*1The (2k +2) th element in (a);
step 12: according toAnd in combination with R ═ BQnBT+DQmDT=[B D]Q[B D]TUpdating the value of R and further updating the value of F; then, solving the semi-definite programming problem by utilizing the interior point method againObtaining the final global optimal solution of Y and the final global optimal solution of Y, and correspondingly recording as Y*2And y*2(ii) a According to the definition of y and y*2To obtain xo、vo、The respective final estimated values, corresponding to x*2、v*2、v*2=y*2(k+1:2k), Then according toAndandto obtainAndis expressed asAnd wherein, y*2(1: k) represents y*2The vector y consisting of the 1 st element to the k th element of (A)*2(k +1:2k) represents y*2The (k +1) th element to the 2 k-th element of (a), y*2(2k +1) represents y*2The (2k +1) th element, y*2(2k +2) represents y*2The (2k +2) th element in (b).
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110208858.7A CN112986907B (en) | 2021-02-25 | 2021-02-25 | Moving target positioning method under clock deviation and clock drift conditions |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110208858.7A CN112986907B (en) | 2021-02-25 | 2021-02-25 | Moving target positioning method under clock deviation and clock drift conditions |
Publications (2)
Publication Number | Publication Date |
---|---|
CN112986907A CN112986907A (en) | 2021-06-18 |
CN112986907B true CN112986907B (en) | 2022-05-17 |
Family
ID=76350349
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202110208858.7A Active CN112986907B (en) | 2021-02-25 | 2021-02-25 | Moving target positioning method under clock deviation and clock drift conditions |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN112986907B (en) |
Families Citing this family (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN113835064B (en) * | 2021-08-13 | 2022-09-23 | 中国人民解放军战略支援部队信息工程大学 | Weighted multi-dimensional scale TDOA (time difference of arrival) positioning method for cooperative correction source observation information |
CN113835061B (en) * | 2021-08-13 | 2022-07-05 | 中国人民解放军战略支援部队信息工程大学 | Single-platform Doppler two-stage closed positioning method in presence of signal carrier frequency prior error |
CN114035181A (en) * | 2021-09-30 | 2022-02-11 | 宁波大学 | Moving target elliptic positioning method for unknown transmitter position and speed |
CN114910864B (en) * | 2022-06-14 | 2023-08-15 | 中国人民解放军战略支援部队信息工程大学 | Multi-platform Doppler positioning method with unknown signal propagation speed and signal frequency drift |
Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103293512A (en) * | 2012-03-02 | 2013-09-11 | 瑞士优北罗股份有限公司 | Positioning using a local wave-propagation model |
CN107765216A (en) * | 2017-08-29 | 2018-03-06 | 宁波大学 | Target location and timing parameter combined estimation method in unsynchronized wireless networks |
KR20180107964A (en) * | 2017-03-23 | 2018-10-04 | 국방과학연구소 | Method and apparatus for positioning using combination of tdoa/fdoa |
CN110988800A (en) * | 2020-02-28 | 2020-04-10 | 浙江万里学院 | Semi-positive relaxation positioning method based on acoustic energy |
CN111929640A (en) * | 2020-06-19 | 2020-11-13 | 浙江万里学院 | Sensor network positioning method under condition of unknown transmitting power |
Family Cites Families (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US7663547B2 (en) * | 2007-04-13 | 2010-02-16 | Glowlink Communications Technology, Inc. | Determining a geolocation solution of an emitter on earth based on weighted least-squares estimation |
-
2021
- 2021-02-25 CN CN202110208858.7A patent/CN112986907B/en active Active
Patent Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103293512A (en) * | 2012-03-02 | 2013-09-11 | 瑞士优北罗股份有限公司 | Positioning using a local wave-propagation model |
KR20180107964A (en) * | 2017-03-23 | 2018-10-04 | 국방과학연구소 | Method and apparatus for positioning using combination of tdoa/fdoa |
CN107765216A (en) * | 2017-08-29 | 2018-03-06 | 宁波大学 | Target location and timing parameter combined estimation method in unsynchronized wireless networks |
CN110988800A (en) * | 2020-02-28 | 2020-04-10 | 浙江万里学院 | Semi-positive relaxation positioning method based on acoustic energy |
CN111929640A (en) * | 2020-06-19 | 2020-11-13 | 浙江万里学院 | Sensor network positioning method under condition of unknown transmitting power |
Non-Patent Citations (5)
Title |
---|
Moving Source Localization Using TOA and FOA Measurements With Imperfect Synchronization;Jiong Shi等;《Signal Processing》;20210406;第1-15页 * |
TOA-based joint synchronization and source localization with random errors in sensor positions and sensor clock biases;Yinggui Wang等;《Ad Hoc Networks》;20141220;第99-111页 * |
基于到达时间与频率估计的高精度多星无源定位方法研究;刘梦竹;《中国优秀博硕士学位论文全文数据库(硕士)工程科技Ⅱ辑》;20200215(第2期);第49-52页 * |
基于到达时间的稳健非视距定位研究;张圣金;《中国优秀博硕士学位论文全文数据库(硕士)社会科学Ⅰ辑》;20180215(第2期);第14-16页 * |
存在站址误差下的时频差稳健定位算法;高向颖 等;《雷达学报》;20201031;第9卷(第5期);第916-924页 * |
Also Published As
Publication number | Publication date |
---|---|
CN112986907A (en) | 2021-06-18 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN112986907B (en) | Moving target positioning method under clock deviation and clock drift conditions | |
CN106842128B (en) | The acoustics tracking and device of moving target | |
Salari et al. | Mobility-aided wireless sensor network localization via semidefinite programming | |
CN111929640B (en) | Sensor network positioning method under unknown transmission power condition | |
CN110658490B (en) | RSS (really simple syndication) and AOA (automatic optical inspection) based three-dimensional wireless sensor network non-cooperative positioning method | |
CN111157943B (en) | TOA-based sensor position error suppression method in asynchronous network | |
CN109581281A (en) | Moving objects location method based on reaching time-difference and arrival rate difference | |
CN110662163A (en) | RSS (really simple syndication) and AOA (automatic optical inspection) based three-dimensional wireless sensor network cooperative positioning method | |
CN103323815A (en) | Underwater acoustic locating method based on equivalent sound velocity | |
CN107765216B (en) | Target position and timing parameter combined estimation method in unsynchronized wireless networks | |
Wang et al. | Iterative constrained weighted least squares estimator for TDOA and FDOA positioning of multiple disjoint sources in the presence of sensor position and velocity uncertainties | |
Wang et al. | A novel estimator for TDOA and FDOA positioning of multiple disjoint sources in the presence of calibration emitters | |
Chen et al. | Improved two-step weighted least squares algorithm for TDOA-based source localization | |
Zhao et al. | Closed-form two-way TOA localization and synchronization for user devices with motion and clock drift | |
Cheng et al. | Direction-of-arrival estimation with virtual antenna array: Observability analysis, local oscillator frequency offset compensation, and experimental results | |
CN114325581A (en) | Elliptical target positioning method with clock synchronization error | |
CN107566981B (en) | Indoor high-precision positioning method, device and system based on optimal path | |
Fadakar et al. | Deep learning aided multi-source passive 3D AOA wireless positioning using a moving receiver: A low complexity approach | |
CN110568406B (en) | Positioning method based on acoustic energy under condition of unknown energy attenuation factor | |
Li et al. | Robust kernel-based machine learning localization using NLOS TOAs or TDOAs | |
CN114051209B (en) | Fingerprint positioning method based on intelligent reflecting surface and scene geometric model | |
Brückner et al. | Phase difference based precise indoor tracking of common mobile devices using an iterative holographic extended Kalman filter | |
CN108415005A (en) | A kind of passive location delay time estimation method and device | |
CN110673088B (en) | Target positioning method based on arrival time in mixed line-of-sight and non-line-of-sight environment | |
CN113923590A (en) | TOA positioning method under condition of uncertain anchor node position |
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 |