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 PDF

Info

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
Application number
CN202110208858.7A
Other languages
Chinese (zh)
Other versions
CN112986907A (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.)
Zhejiang Wanli University
Original Assignee
Zhejiang Wanli University
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 Zhejiang Wanli University filed Critical Zhejiang Wanli University
Priority to CN202110208858.7A priority Critical patent/CN112986907B/en
Publication of CN112986907A publication Critical patent/CN112986907A/en
Application granted granted Critical
Publication of CN112986907B publication Critical patent/CN112986907B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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

Moving target positioning method under clock deviation and clock drift conditions
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 as
Figure BDA0002951693760000031
Wherein 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,
Figure BDA0002951693760000032
representing the speed of the 1 st sensor node,
Figure BDA0002951693760000033
representing the speed of the 2 nd sensor node,
Figure BDA0002951693760000034
representing the speed of the ith sensor node,
Figure BDA0002951693760000035
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:
Figure BDA0002951693760000036
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:
Figure BDA0002951693760000037
wherein the symbol "| | |" is a euclidean distance-solving symbol, c represents the propagation speed of the wireless signal,
Figure BDA0002951693760000038
representing unknown clock deviations, f0Representing the carrier frequency, the superscript symbol "T" representing the transpose of a vector or matrix,
Figure BDA0002951693760000039
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:
Figure BDA00029516937600000310
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:
Figure BDA0002951693760000041
Figure BDA0002951693760000042
wherein the content of the first and second substances,
Figure BDA0002951693760000043
ni＝Δτic，
Figure BDA0002951693760000044
Figure BDA0002951693760000045
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，
Figure BDA0002951693760000046
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,
Figure BDA00029516937600000417
and
Figure BDA00029516937600000418
according to
Figure BDA0002951693760000047
The result of the calculation is that,
Figure BDA0002951693760000048
and
Figure BDA0002951693760000049
according to
Figure BDA00029516937600000410
Calculating 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 to
Figure BDA00029516937600000411
Calculating to obtain;
and 5: to pair
Figure BDA00029516937600000412
Make item shifting and conversion into
Figure BDA00029516937600000413
Then to
Figure BDA00029516937600000414
Squaring both sides and ignoring the noise squared term
Figure BDA00029516937600000415
To obtain
Figure BDA00029516937600000416
And will be
Figure BDA0002951693760000051
Substitution into
Figure BDA0002951693760000052
In (1) obtaining
Figure BDA0002951693760000053
Re-spread
Figure BDA0002951693760000054
Then, the term shift is carried out, and the second-order noise term n is ignoredimiTo obtain
Figure BDA0002951693760000055
Step 6: introducing an auxiliary vector yo，yoIs defined as
Figure BDA0002951693760000056
Then will be
Figure BDA0002951693760000057
Rewritten in matrix form, described as: a. the1yo-h1≈B1n; and will be
Figure BDA0002951693760000058
Rewritten 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:
Figure BDA0002951693760000059
wherein the content of the first and second substances,
Figure BDA00029516937600000510
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,
Figure BDA00029516937600000511
h1has the dimension of N x 1,
Figure BDA00029516937600000512
Figure BDA00029516937600000513
show that
Figure BDA00029516937600000514
The diagonal matrix is composed as diagonal elements,
Figure BDA00029516937600000515
A2has a dimension of N x (2k +6),
Figure BDA0002951693760000061
h2has the dimension of N x 1,
Figure BDA0002951693760000062
Figure BDA0002951693760000063
show that
Figure BDA0002951693760000064
The diagonal matrix is composed as diagonal elements,
Figure BDA0002951693760000065
show that
Figure BDA0002951693760000066
The diagonal matrix is composed as diagonal elements,
Figure BDA0002951693760000067
the dimension of A is 2N x (2k +6),
Figure BDA0002951693760000068
the dimension of h is 2N x 1,
Figure BDA0002951693760000069
the dimension of B is 2N x N,
Figure BDA00029516937600000610
dimension of D is 2 NxN, 0N×NRepresenting a matrix of dimension N x N and elements all 0,
Figure BDA00029516937600000611
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:
Figure BDA00029516937600000612
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) problem
Figure BDA0002951693760000071
Is relaxed to
Figure BDA0002951693760000072
Wherein the content of the first and second substances,
Figure BDA0002951693760000073
express to make
Figure BDA0002951693760000074
Y and Y, tr () represents the trace of the matrix,
Figure BDA0002951693760000075
to represent
Figure BDA0002951693760000076
Is semi-positive, rank () represents the rank of the matrix,
Figure BDA0002951693760000077
and
Figure BDA0002951693760000078
derived from the equivalent equation:
Figure BDA0002951693760000079
and step 9: problems after loosening of the semipositive fixation
Figure BDA0002951693760000081
The 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 the
Figure BDA0002951693760000091
In (1)
Figure BDA0002951693760000092
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:
Figure BDA0002951693760000093
step 11: let R initial value be an identity matrix I with dimension of 2 Nx 2N2N(ii) a Then according to
Figure BDA0002951693760000094
Calculating to obtain the value of F; then solving the semi-definite programming problem by using an interior point method
Figure BDA0002951693760000101
Obtaining 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 obtain
Figure BDA0002951693760000102
The respective preliminary estimate values are correspondingly recorded as
Figure BDA0002951693760000103
Wherein, 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 to
Figure BDA0002951693760000104
And 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 again
Figure BDA0002951693760000105
Obtaining 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、
Figure BDA0002951693760000106
The respective final estimated values, corresponding to x*2、v*2、
Figure BDA0002951693760000107
x*2＝y*2(1:k)，v*2＝y*2(k+1:2k)，
Figure BDA0002951693760000108
Figure BDA0002951693760000109
Then according to
Figure BDA00029516937600001010
And
Figure BDA00029516937600001011
and
Figure BDA00029516937600001012
to obtain
Figure BDA00029516937600001013
And
Figure BDA00029516937600001014
is expressed as
Figure BDA00029516937600001015
And
Figure BDA00029516937600001016
Figure BDA00029516937600001017
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,
Figure BDA0002951693760000111
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 on
Figure BDA0002951693760000112
A profile of change;
figure 3 shows that the number of sensor nodes is 5,
Figure BDA0002951693760000113
from 10-2m2Change to 103m2In the method of the inventionRMSE Performance of velocity estimation by (SDP for short), ML method, Raw WLS method
Figure BDA0002951693760000114
A profile of change;
figure 4 shows that the number of sensor nodes is 5,
Figure BDA0002951693760000121
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 followed
Figure BDA0002951693760000122
A profile of change;
figure 5 shows that the number of sensor nodes is 5,
Figure BDA0002951693760000123
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 invention
Figure BDA0002951693760000124
A profile of change;
FIG. 6 is a drawing showing
Figure BDA0002951693760000125
Fixed 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 showing
Figure BDA0002951693760000126
Fixed 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;
FIG. 8 is a diagram of sensingThe number of the node is 5 and,
Figure BDA0002951693760000127
fixed at 10m2FOA measurement noise coefficient η 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.
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 as
Figure BDA0002951693760000128
Wherein 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,
Figure BDA0002951693760000131
representing the speed of the 1 st sensor node,
Figure BDA0002951693760000132
representing the speed of the 2 nd sensor node,
Figure BDA0002951693760000133
representing the speed of the ith sensor node,
Figure BDA0002951693760000134
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:
Figure BDA0002951693760000135
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:
Figure BDA0002951693760000136
wherein the symbol "| | |" is a euclidean distance-solving symbol, c represents the propagation speed of the wireless signal,
Figure BDA0002951693760000137
representing unknown clock deviations, f0Representing the carrier frequency, the superscript symbol "T" representing the transpose of a vector or matrix,
Figure BDA0002951693760000138
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:
Figure BDA0002951693760000139
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:
Figure BDA0002951693760000141
Figure BDA0002951693760000142
wherein the content of the first and second substances,
Figure BDA0002951693760000143
ni＝Δτic，
Figure BDA0002951693760000144
Figure BDA0002951693760000145
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，
Figure BDA0002951693760000146
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,
Figure BDA0002951693760000147
and
Figure BDA0002951693760000148
according to
Figure BDA0002951693760000149
The calculation results in that,
Figure BDA00029516937600001410
and
Figure BDA00029516937600001411
according to
Figure BDA00029516937600001412
Calculating 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 to
Figure BDA00029516937600001413
And (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,
Figure BDA00029516937600001414
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:
Figure BDA0002951693760000151
wherein theta is thetaoCorresponding optimization variable, θ ═ xT,vT,d0,b0]TX is xoCorresponding optimization variable, v is voCorresponding optimization variable, d0Is composed of
Figure BDA0002951693760000152
Corresponding optimization variables, b0Is composed of
Figure BDA0002951693760000153
The corresponding optimization variables are set to be the same as,
Figure BDA0002951693760000154
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,
Figure BDA0002951693760000155
express to make
Figure BDA0002951693760000156
The value of theta at which the value of (b) is minimum,
Figure BDA0002951693760000157
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 right
Figure BDA0002951693760000158
Make item shifting and transformation into
Figure BDA0002951693760000159
Then to
Figure BDA00029516937600001510
Squaring both sides and ignoring the noise squared term
Figure BDA00029516937600001511
To obtain
Figure BDA00029516937600001512
And will be
Figure BDA00029516937600001513
Substitution into
Figure BDA00029516937600001514
In (1) obtaining
Figure BDA00029516937600001515
Re-spread
Figure BDA00029516937600001516
Then, the term shift is carried out, and the second-order noise term n is ignoredimiTo obtain
Figure BDA00029516937600001517
And 6: introducing an auxiliary vector yo，yoIs defined as
Figure BDA00029516937600001518
Then will be
Figure BDA00029516937600001519
Rewritten in matrix form, described as: a. the1yo-h1≈B1n; and will be
Figure BDA00029516937600001520
Rewritten 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:
Figure BDA0002951693760000161
wherein the content of the first and second substances,
Figure BDA0002951693760000162
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,
Figure BDA0002951693760000163
h1has the dimension of (a) of N x 1,
Figure BDA0002951693760000164
Figure BDA0002951693760000165
show that
Figure BDA0002951693760000166
A diagonal matrix is composed as diagonal elements,
Figure BDA0002951693760000167
A2has a dimension of N x (2k +6),
Figure BDA0002951693760000168
h2has the dimension of N x 1,
Figure BDA0002951693760000169
Figure BDA00029516937600001610
show that
Figure BDA00029516937600001611
The diagonal matrix is composed as diagonal elements,
Figure BDA00029516937600001612
show that
Figure BDA00029516937600001613
The diagonal matrix is composed as diagonal elements,
Figure BDA00029516937600001614
the dimension of A is 2N x (2k +6),
Figure BDA00029516937600001615
the dimension of h is 2N x 1,
Figure BDA00029516937600001616
the dimension of B is 2N x N,
Figure BDA00029516937600001617
dimension of D is 2 NxN, 0N×NRepresenting a matrix of dimension N x N and elements all 0,
Figure BDA00029516937600001618
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:
Figure BDA0002951693760000171
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 problem
Figure BDA0002951693760000172
Is relaxed to
Figure BDA0002951693760000181
Wherein the content of the first and second substances,
Figure BDA0002951693760000182
express to make
Figure BDA0002951693760000183
Y and Y, tr () represents the trace of the matrix,
Figure BDA0002951693760000184
to represent
Figure BDA0002951693760000185
Is semi-positive, rank () represents the rank of the matrix,
Figure BDA0002951693760000186
and
Figure BDA0002951693760000187
derived from the equivalent equation:
Figure BDA0002951693760000188
and step 9: problems after loosening of the semipositive fixation
Figure BDA0002951693760000189
The 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 the
Figure BDA0002951693760000191
In (1)
Figure BDA0002951693760000192
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:
Figure BDA0002951693760000201
step 11: let R initial value be an identity matrix I with dimension of 2 Nx 2N2N(ii) a Then according to
Figure BDA0002951693760000202
Calculating to obtain the value of F; then solving the semi-definite programming problem by using an interior point method
Figure BDA0002951693760000203
Obtaining 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 obtain
Figure BDA0002951693760000204
The respective preliminary estimate values are correspondingly recorded as
Figure BDA0002951693760000205
Wherein, 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 to
Figure BDA0002951693760000206
And 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 again
Figure BDA0002951693760000211
Obtaining 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、
Figure BDA0002951693760000212
The respective final estimated values, corresponding to x*2、v*2、
Figure BDA0002951693760000213
x*2＝y*2(1:k)，v*2＝y*2(k+1:2k)，
Figure BDA0002951693760000214
Figure BDA0002951693760000215
Then according to
Figure BDA0002951693760000216
And
Figure BDA0002951693760000217
and
Figure BDA0002951693760000218
to obtain
Figure BDA0002951693760000219
And
Figure BDA00029516937600002110
is expressed as
Figure BDA00029516937600002111
And
Figure BDA00029516937600002112
Figure BDA00029516937600002113
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
Figure BDA00029516937600002114
Figure BDA0002951693760000221
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
Figure BDA0002951693760000222
Figure BDA0002951693760000223
Randomly generated (in meters) from a uniformly distributed U (0,150), caused by clock drift
Figure BDA0002951693760000224
Figure BDA0002951693760000225
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 respectively
Figure BDA0002951693760000226
And
Figure BDA0002951693760000227
Figure BDA0002951693760000228
represents the variance of the noise of the TOA measurement,
Figure BDA0002951693760000229
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 used
Figure BDA00029516937600002210
The 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:
Figure BDA00029516937600002211
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,
Figure BDA0002951693760000231
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 ignored
Figure BDA0002951693760000232
Then become di＝||xo-si||+niI 1, …, N, in the RROA observation vector, if clock drift is neglected
Figure BDA0002951693760000233
Then become into
Figure BDA0002951693760000234
Referring to step 5, the ROA and RROA are expanded and the second order noise terms are ignored to obtain
Figure BDA0002951693760000235
And
Figure BDA0002951693760000236
thereafter, referring to the method of step 6, an auxiliary vector is introduced
Figure BDA0002951693760000237
Is defined as
Figure BDA0002951693760000238
Then will be
Figure BDA0002951693760000239
Rewritten in matrix form, described as:
Figure BDA00029516937600002310
and will be
Figure BDA00029516937600002311
Rewritten in matrix form, described as:
Figure BDA00029516937600002312
then are connected together
Figure BDA00029516937600002313
And
Figure BDA00029516937600002314
to obtain
Figure BDA0002951693760000241
Then according to
Figure BDA0002951693760000242
Constructing a weighted least squares WLS problem, described as:
Figure BDA0002951693760000243
wherein
Figure BDA0002951693760000244
A3Has 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,
Figure BDA0002951693760000245
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,
Figure BDA0002951693760000246
A4has a dimension of N x (2k +2),
Figure BDA0002951693760000247
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,
Figure BDA0002951693760000248
ARWLShas a dimension of 2N x (2k +2),
Figure BDA0002951693760000249
hRWLShas a dimension of 2N x 1,
Figure BDA00029516937600002410
BRWLShas a dimension of 2N x N,
Figure BDA00029516937600002411
DRWLShas a dimension of 2 NxN, 0N×NRepresenting a matrix of dimension N x N and elements all 0,
Figure BDA00029516937600002412
express to make
Figure BDA00029516937600002413
Y is the smallest value ofRWLSValue of (a), yRWLSIs composed of
Figure BDA00029516937600002414
Corresponding optimization variable, RRWLSA matrix of weights is represented by a matrix of weights,
Figure BDA00029516937600002415
Figure BDA00029516937600002416
represents RRWLSInverse of (B)RWLS DRWLS]A matrix with dimensions 2N × 2N.
Scene 1: the sensor nodes take the first 5 in table 1,
Figure BDA0002951693760000251
from 10-2m2Change to 103m2。
FIG. 2 shows the RMSE performance of the location estimation of the present invention method (SDP for short), ML method, and Raw WLS method
Figure BDA0002951693760000252
The 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 of
Figure BDA0002951693760000253
The 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 noise
Figure BDA0002951693760000254
The 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 method
Figure BDA0002951693760000255
The 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 method
Figure BDA0002951693760000256
The 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 the
Figure BDA0002951693760000257
Corresponding 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 field
Figure BDA0002951693760000258
Is 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
Figure BDA0002951693760000259
Figure BDA0002951693760000261
Scene 2:
Figure BDA0002951693760000262
fixed at 10m2The number of sensor nodes increases from the first 4 to 9 in table 1.
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 noise
Figure BDA0002951693760000271
Fixed 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 as
Figure FDA0002951693750000011
Wherein 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,
Figure FDA0002951693750000012
representing the speed of the 1 st sensor node,
Figure FDA0002951693750000013
representing the speed of the 2 nd sensor node,
Figure FDA0002951693750000014
representing the speed of the ith sensor node,
Figure FDA0002951693750000015
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:
Figure FDA0002951693750000016
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:
Figure FDA0002951693750000017
wherein the symbol "| | |" is a euclidean distance-solving symbol, c represents the propagation speed of the wireless signal,
Figure FDA0002951693750000018
representing unknown clock deviations, f0Representing the carrier frequency, the superscript symbol "T" representing the transpose of a vector or matrix,
Figure FDA0002951693750000019
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:
Figure FDA0002951693750000021
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:
Figure FDA0002951693750000022
Figure FDA0002951693750000023
wherein the content of the first and second substances,
Figure FDA0002951693750000024
ni＝Δτic，
Figure FDA0002951693750000025
Figure FDA0002951693750000026
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，
Figure FDA0002951693750000027
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,
Figure FDA0002951693750000028
and
Figure FDA0002951693750000029
according to
Figure FDA00029516937500000210
The calculation results in that,
Figure FDA00029516937500000211
and
Figure FDA00029516937500000212
according to
Figure FDA00029516937500000213
Calculating 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 to
Figure FDA0002951693750000031
Calculating to obtain;
and 5: to pair
Figure FDA0002951693750000032
Make item shifting and conversion into
Figure FDA0002951693750000033
Then to
Figure FDA0002951693750000034
Squaring both sides and ignoring the noise squared term
Figure FDA0002951693750000035
To obtain
Figure FDA0002951693750000036
And will be
Figure FDA0002951693750000037
Substitution into
Figure FDA0002951693750000038
In (1) obtaining
Figure FDA0002951693750000039
Re-spread
Figure FDA00029516937500000310
Then, the term shift is carried out, and the second-order noise term n is ignoredimiTo obtain
Figure FDA00029516937500000311
Step 6: introducing an auxiliary vector yo，yoIs defined as
Figure FDA00029516937500000312
Then will be
Figure FDA00029516937500000313
Rewritten in matrix form, described as: a. the1yo-h1≈B1n; and will be
Figure FDA00029516937500000314
Rewritten 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:
Figure FDA00029516937500000315
wherein the content of the first and second substances,
Figure FDA00029516937500000316
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，
Figure FDA0002951693750000041
h1Has the dimension of (a) of N x 1,
Figure FDA0002951693750000042
Figure FDA0002951693750000043
show that
Figure FDA0002951693750000044
The diagonal matrix is composed as diagonal elements,
Figure FDA0002951693750000045
A2has a dimension of N x (2k +6),
Figure FDA0002951693750000046
h2has the dimension of N x 1,
Figure FDA0002951693750000047
Figure FDA0002951693750000048
show that
Figure FDA0002951693750000049
The diagonal matrix is composed as diagonal elements,
Figure FDA00029516937500000410
Figure FDA00029516937500000411
show that
Figure FDA00029516937500000412
The diagonal matrix is composed as diagonal elements,
Figure FDA00029516937500000413
the dimension of A is 2N x (2k +6),
Figure FDA00029516937500000414
the dimension of h is 2N x 1,
Figure FDA00029516937500000415
the dimension of B is 2N x N,
Figure FDA00029516937500000416
dimension of D is 2 NxN, 0N×NRepresenting a matrix of dimension N x N and elements all 0,
Figure FDA00029516937500000417
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:
Figure FDA0002951693750000051
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 problem
Figure FDA0002951693750000052
Is relaxed to
Figure FDA0002951693750000053
Wherein the content of the first and second substances,
Figure FDA0002951693750000054
express to make
Figure FDA0002951693750000061
Y and Y, tr () represents the trace of the matrix,
Figure FDA0002951693750000062
Figure FDA0002951693750000063
to represent
Figure FDA0002951693750000064
Is semi-positive, rank () represents the rank of the matrix,
Figure FDA0002951693750000065
and
Figure FDA0002951693750000066
derived from the equivalent equation:
Figure FDA0002951693750000067
and step 9: the problem after semi-positive relaxation is described
Figure FDA0002951693750000068
The 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 the
Figure FDA0002951693750000071
In
Figure FDA0002951693750000072
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:
Figure FDA0002951693750000073
step 11: let R initial value be an identity matrix I with dimension of 2 Nx 2N2N(ii) a Then according to
Figure FDA0002951693750000074
Calculating to obtain the value of F; then solving the semi-definite programming problem by using an interior point method
Figure FDA0002951693750000081
Obtaining 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 obtain
Figure FDA0002951693750000082
Preliminary estimate of each, pairShould be noted as
Figure FDA0002951693750000083
Wherein, 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 to
Figure FDA0002951693750000084
And 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 again
Figure FDA0002951693750000085
Obtaining 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、
Figure FDA0002951693750000086
The respective final estimated values, corresponding to x*2、v*2、
Figure FDA0002951693750000087
v*2＝y*2(k+1:2k)，
Figure FDA0002951693750000088
Figure FDA0002951693750000089
Then according to
Figure FDA00029516937500000810
And
Figure FDA00029516937500000811
and
Figure FDA00029516937500000812
to obtain
Figure FDA00029516937500000813
And
Figure FDA00029516937500000814
is expressed as
Figure FDA00029516937500000815
And
Figure FDA00029516937500000816
Figure FDA00029516937500000817
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).
CN202110208858.7A 2021-02-25 2021-02-25 Moving target positioning method under clock deviation and clock drift conditions Active CN112986907B (en)

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 (3)

* Cited by examiner, † Cited by third party
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
CN114910864B (en) * 2022-06-14 2023-08-15 中国人民解放军战略支援部队信息工程大学 Multi-platform Doppler positioning method with unknown signal propagation speed and signal frequency drift

Citations (5)

* Cited by examiner, † Cited by third party
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)

* Cited by examiner, † Cited by third party
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

Patent Citations (5)

* Cited by examiner, † Cited by third party
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)

* Cited by examiner, † Cited by third party
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
CN110658490B (en) RSS (really simple syndication) and AOA (automatic optical inspection) based three-dimensional wireless sensor network non-cooperative positioning method
CN111929640B (en) Sensor network positioning method under unknown transmission power condition
CN110662163A (en) RSS (really simple syndication) and AOA (automatic optical inspection) based three-dimensional wireless sensor network cooperative positioning method
Guo et al. Multi-source localization using time of arrival self-clustering method in wireless sensor 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
Chen et al. Improved two-step weighted least squares algorithm for TDOA-based source localization
CN111698695A (en) LTE fingerprint type positioning method based on neural network
Wang et al. A novel estimator for TDOA and FDOA positioning of multiple disjoint sources in the presence of calibration emitters
Cheng et al. Direction-of-arrival estimation with virtual antenna array: Observability analysis, local oscillator frequency offset compensation, and experimental results
WO2016095694A1 (en) Improved source localization algorithm in which sensor error is present
Ferreira et al. Optimal positioning of autonomous marine vehicles for underwater acoustic source localization using TOA measurements
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
CN107566981B (en) Indoor high-precision positioning method, device and system based on optimal path
CN114325581A (en) Elliptical target positioning method with clock synchronization error
Brückner et al. Phase difference based precise indoor tracking of common mobile devices using an iterative holographic extended Kalman filter
Tao et al. Guaranteed stability of sparse recovery in distributed compressive sensing MIMO radar
CN110673088B (en) Target positioning method based on arrival time in mixed line-of-sight and non-line-of-sight environment
Fadakar et al. Deep learning aided multi-source passive 3D AOA wireless positioning using a moving receiver: A low complexity approach
CN113923590A (en) TOA positioning method under condition of uncertain anchor node position
CN110988800A (en) Semi-positive relaxation positioning method based on acoustic energy
CN112986913A (en) Underwater target positioning method based on differential Doppler and arrival time delay difference

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