Uniform-speed moving target positioning method under asynchronous transmitter and unknown position

Abstract

The invention discloses a uniform-speed moving target positioning method under asynchronous transmitters and unknown positions, which aims at a multi-base positioning scene, utilizes time delay information of continuous moments to replace a traditional time delay and Doppler information mixing method, and reduces system hardware requirements. Aiming at a distance measurement model, eliminating part of coupling variables by using numerical characteristics between measurement values, then separating unknown variables, and converting a positioning problem into a weighted least square problem with constraints; performing first-order Taylor expansion on a target function of the problem to obtain a weighted least square problem with constraint and bias reduction capability; converting the problem into a convex positive semidefinite programming problem; and finally obtaining the optimal estimation values of the initial coordinates and the moving speed of the target. The method has the advantages that under the condition that the transmitter is asynchronous and the position is unknown, the target position and speed are estimated through distance measurement values observed at continuous time points, the equipment cost is low, the calculation complexity is low, and the obvious deviation reduction capacity is achieved.

Description

Technical Field

[0001] The invention relates to a target positioning technology, in particular to a method for positioning a uniformly moving target when a transmitter is asynchronous and the position is unknown. Background Art

[0002] In the field of target positioning, the positioning of mobile targets is more challenging than that of stationary targets, especially in application scenarios such as target tracking, sensor networks, search and rescue, etc. Traditional positioning technology mainly relies on measurements such as time delay (TD), angle of arrival (AOA), and received signal strength (RSS). In order to estimate the speed of mobile targets more accurately, modern positioning systems must also use Doppler frequency shift (DFS) information.

[0003] In recent years, mobile target positioning methods that combine time and frequency information have been widely studied and developed. These methods mainly include the joint application of Time Different of Arrival (TDOA) and Frequency Different of Arrival (FDOA).

[0004] Based on these studies, a multi-base station positioning model suitable for mobile targets was proposed. These models usually use passive coherent positioning technology to estimate the position and speed of mobile targets by analyzing the time delay and Doppler frequency shift information in the signals transmitted from base stations such as TV towers and radio towers and the echo signals reflected by the targets. Since the transmission sources such as TV towers and radio towers are not designed for positioning, the system usually cannot directly obtain their exact location. However, through the timing function of signal transmission, the corresponding signal transmission time can be obtained, thereby indirectly obtaining the time delay information. The main advantage of multi-base station positioning technology is that it can use existing infrastructure without the need to deploy additional dedicated transmitting and receiving equipment, thereby reducing the cost of the system.

[0005] Multi-base station positioning technology not only performs well in common application scenarios, but also has great potential in scenarios requiring high-precision and wide-coverage positioning. The multi-base positioning system achieves accurate positioning of mobile targets through the collaborative work of multiple transmitters and receivers, using a variety of measurement methods such as time difference of arrival (TDOA) and Doppler frequency shift (DFS).

[0006] Although the introduction of Doppler frequency shift information can significantly improve the reliability of the estimation of the position and velocity of the moving target, this method is mainly applicable to narrowband signals. For broadband signals, the accuracy of frequency information is often not enough to meet the requirements of precise positioning. In order to overcome this limitation, some studies have proposed a model that assumes that the moving target moves at a constant speed in a short period of time, and determines the initial position and velocity of the moving target by collecting corresponding distance or angle measurements at multiple consecutive time points, thereby depicting the target's motion trajectory; however, with the introduction of multi-time measurement values, the complexity of existing positioning methods has also increased.

[0007] At present, although there have been some studies on multi-base mobile target positioning, there has been no research on estimating the position and speed of mobile targets based on the multi-base positioning model using distance measurements observed at consecutive time points. This new positioning problem has challenged existing positioning technology, but also provided an important opportunity for the development of new adaptive positioning algorithms, which has high research significance. Summary of the invention

[0008] The technical problem to be solved by the present invention is to provide a method for locating a uniformly moving target when the transmitter is asynchronous and the position is unknown. When the transmitter is asynchronous and the position is unknown, the distance measurement values ​​observed at continuous time points are used to estimate the target position and speed, which solves the problem that target speed estimation cannot be achieved due to equipment cost limitations. The method has low computational complexity and the ability to reduce deviations.

[0009] The technical solution adopted by the present invention to solve the above technical problems is: a method for locating a uniformly moving target with an asynchronous transmitter and an unknown position, characterized by comprising the following steps:

[0010] Step 1: In a multi-base positioning scenario, assume that there are several sensors and a target. Some sensors are receivers with known coordinate positions. Different receivers have the same reference clock, that is, the clocks of different receivers are synchronized. The remaining sensors are transmitters with unknown coordinate positions. The transmitter and receiver are asynchronous in time, that is, the reference clocks are different. There is a clock synchronization error between the transmitter and the receiver. The initial coordinate position of the target is unknown and it moves at an unknown constant speed.

[0011] Step 2: At different time points, the signal transmitted by the transmitter reaches the receiver via a direct path or an indirect path reflected by the target. The receiver can obtain the direct path distance measurement value and the indirect path distance measurement value observed at the continuous time points, and accordingly construct a distance measurement model, including a direct path distance measurement part and an indirect path distance measurement part; wherein, the signal transmitted by the transmitter carries the clock information when the signal is transmitted;

[0012] Step 3: For the indirect path distance measurement part of the distance measurement model, the indirect path distance measurement values ​​obtained by the remaining receivers except the first receiver are subtracted from the indirect path distance measurement value obtained by the first receiver, and the unknown variables are separated through mathematical processing; for the direct path distance measurement part of the distance measurement model, the unknown variables are separated through mathematical processing; then, all unknown variables are used to transform the distance measurement model after the unknown variables are separated into the original optimization problem, which is a weighted least squares problem with constraints;

[0013] Step 4: Perform a first-order Taylor expansion on the objective function in the original optimization problem; then obtain a constraint on the deviation term based on the first-order Taylor expansion; then combine the original optimization problem and the constraint on the deviation term to obtain a constrained weighted least squares problem with the ability to reduce deviation;

[0014] Step 5: Use the semi-positive definite relaxation technique to relax the constrained weighted least squares problem with bias reduction capability into a convex semi-positive definite programming problem that is easy to solve;

[0015] Step 6: Directly solve the convex semi-positive programming problem to obtain the optimal estimate of the initial coordinate position of the target, the optimal estimate of the target's moving speed, the optimal estimate of the coordinate position of each transmitter, and the optimal estimate of the clock synchronization error of each transmitter.

[0016] In step 1, the multi-base positioning scenario is two-dimensional or three-dimensional, the number of receivers is N, and the coordinate position of the nth receiver is s n , n=1,2,…,N, the number of transmitters is M, and the coordinate position of the mth transmitter is The clock synchronization error of the mth transmitter is The initial coordinate position of the target is p o , the target's moving speed is v o .

[0017] In step 2, the constructed distance measurement model is described as: Among them, K represents the total number of time points, r mn,k represents the indirect path distance measurement value obtained when the signal transmitted by the mth transmitter at time k reaches the nth receiver through the indirect path reflected by the target, d mn,k represents the direct path distance measurement value obtained when the signal transmitted by the mth transmitter at time k reaches the nth receiver through the direct path, ε rmn,k Represents r mn,k The corresponding measurement noise, ε dmn,k Indicates d mn,k The corresponding measurement noise, p o +kv orepresents the coordinate position of the target at time k, and “||||” is the Euclidean norm operator.

[0018] In step 3, the specific process of separating the unknown variables in the indirect path distance measurement part of the distance measurement model is as follows:

[0019] Step 3.1a: Add the indirect path distance measurement part of the distance measurement model

[0020]

[0021] Split into two types, namely

[0022] and

[0023]

[0024] Step 3.2a: Subtract the two equations obtained in step 3.1a to obtain

[0025]

[0026] Step 3.3a: For r in step 3.1a m1,k In the expression Perform mathematical processing of transposing the terms, and then perform mathematical processing of square expansion on the result after transposing the terms, ignoring the second-order noise term, and obtain: The unknown variables are separated, where The superscript "T" indicates the transpose operation of a matrix or vector; for r in step 3.2a mn,k -r m1,k ||p in the expression o +kv o -s 1 ||Perform mathematical processing of transposing terms, and then perform mathematical processing of square expansion on the result after transposing terms to obtain: The unknown variables are separated, where

[0027] In step 3, the specific process of separating the unknown variables in the direct path distance measurement part of the distance measurement model is as follows:

[0028] Direct path distance measurement part of the distance measurement model In Perform mathematical processing of transposing the terms, and then perform mathematical processing of square expansion on the result after transposing the terms, ignoring the second-order noise term, and obtain: The unknown variables are separated, where

[0029] In step 3, the process of obtaining the original optimization problem is as follows:

[0030] Step 3.1b: Combine equation 1, equation 2, and equation 3 to obtain the distance measurement model after separating the unknown variables, which is described as: Then all unknown variables are expressed as vectors, denoted as x o , in, All are coupled variables;

[0031] Step 3.2b: Utilize x o The distance measurement model after separating the unknown variables is transformed into the original optimization problem. The original optimization problem is a weighted least squares problem with constraints, which is described as: Among them, x is the optimization variable in the original optimization problem, x and x o Corresponding, (g-Hx) T W(g-Hx) is the objective function in the original optimization problem, L represents the dimension of the multi-base positioning scenario, L is 2 or 3, and x is defined i Represents the i-th element of x, and defines x (i:j) represents the column vector consisting of the i-th element to the j-th element of x, B=blkdiag(B 1 ,...,B K ),B k =blkdiag(B r1,k ,B r2,k ,B d ), Definition I i represents the identity matrix of dimension i×i, Definition 1 i represents a vector of all 1s of length i, and defines 0 i×j represents an all-zero matrix of dimension i×j, Diag() represents the element diagonalization operation function, that is, the vectors in it are arranged diagonally to form a matrix, blkdiag() represents the matrix diagonalization operation function, that is, each matrix in it is arranged diagonally to form a matrix, Q k Represents ε k The corresponding covariance matrix, ε k represents the composite noise vector at time k, ε k Subject to the mean of 0 and the covariance matrix is ​​Q k Gaussian distribution, rn,k =[r 1n,k ,...,r Mn,k ] T , d n,k =[d 1n,k ,...,d Mn,k ] T , Represented by the matrix The submatrix consisting of the M+1th row to the MNth row and all columns, is the Kronecker product operator.

[0032] The specific process of step 4 is as follows:

[0033] Step 4.1: Define optimization variables in, α is an intermediate variable; then the objective function (g-Hx) in the original optimization problem is T W(g-Hx) is transformed into y T R T WRy, where R represents the coefficient matrix carrying noise, R = [-H, g]; then the objective function y T R T WRy is expanded by the first order Taylor to obtain the first order Taylor expansion, which is described as: in, It means to find the mathematical expectation, R o Represents the coefficient matrix without noise, r mn,k The corresponding indirect path distance measurement without noise d mn,k The corresponding direct path distance measurement without noise Substitute r in R = [-H, g] mn,k d mn,k Get R o , “⊙” is the Hadamard product operator, ΔR=[-ΔH,Δg], A k,r1 =[0 M×(M+2)L ,I M ,0 M×(k-1) ,1 M ,0 M×(K-k) ,0 M×M(K+1) ], A k,r2 =[0 M(N-1)×((M+2)L+(M+k-1)) ,-1 M(N-1) ,0 M(N-1)×(K-k) ,0 M(N-1)×M(K+1) ], represents the composite noise vector after deformation at K moments, represents the composite noise vector after deformation at time k,

[0034] Step 4.2: Based on the first-order Taylor expansion, introduce Form a constraint on the deviation term, described as: in,

[0035] Step 4.3: Combining the original optimization problem and the constraint on the bias term, we get a constrained weighted least squares problem with bias reduction capability, which is described as: Among them, y end Represents the last element of y, defining y (i:j) represents the column vector consisting of the i-th element to the j-th element of y, and defines y i Represents the i-th element of y.

[0036] In step 5, the convex semi-positive definite programming problem is described as: Where Y represents the newly introduced optimization variable, and Y=yy T ,tr{R T WRY} is the objective function in the convex semi-positive definite programming problem, tr{·} means taking the diagonal elements of the matrix and summing them, Y is the optimization variable of the convex semi-positive definite programming problem, and is defined as It means from the i-th row to the j-th row and the l-th column to the The submatrix consisting of the elements corresponding to the columns is defined as Y (i,j) represents the element corresponding to the i-th row and j-th column of Y, Y (end,end) is the index of the last row or column of Y. Y ≥ 0 means that the matrix Y is positive semidefinite.

[0037] In step 6, the initial coordinate position p of the target is o The optimal estimate of * , Set the target's moving speed v o The optimal estimate of * , The coordinate position of the mth transmitter The optimal estimate of The clock synchronization error of the mth transmitter is The optimal estimate of in, Y * represents the optimal estimate of Y, and defines Represents the submatrix consisting of the elements corresponding to the i-th row to the j-th row and the last column of Y, and defines Represents the element corresponding to the i-th row and j-th column of Y.

[0038] Compared with the prior art, the advantages of the present invention are:

[0039] 1) The method of the present invention considers a new positioning problem, in which the distance measurement values ​​of continuous time points are used instead of the traditional use of distance and Doppler mixed measurement values ​​to simultaneously estimate the target position and moving speed; it avoids the inability to estimate the moving target speed due to the performance limitation of the signal acquisition equipment, and at the same time reduces the requirements for signal parameters compared with the distance and Doppler mixed method.

[0040] 2) The method of the present invention utilizes the numerical relationship between the distance measurement values ​​to eliminate some coupled variables, thereby constructing a weighted least squares problem with constraints with a smaller number of coupled unknown variables; the reduction of the dimension of unknown variables can effectively reduce the search dimension of the problem, thereby reducing the computational complexity of the problem and facilitating the tightening of the optimization problem.

[0041] 3) The MSE analysis conducted by the present invention shows that the estimation of unknown variables by the constrained weighted least squares problem with bias reduction capability before relaxation under Gaussian noise can be close to the performance of CRLB; it shows that the method of the present invention can effectively estimate the target position and target velocity in this positioning scenario, and at the same time, the elimination of coupling variables and reduction of estimation bias processing applied in forming the optimization problem optimizes the results without reducing the accuracy of the estimation results; in addition, a theoretical deviation expression is derived to evaluate the size of the deviation. BRIEF DESCRIPTION OF THE DRAWINGS

[0042] Figure 1 It is a general implementation block diagram of the method of the present invention;

[0043] Figure 2 It is a schematic diagram of the system model of the present invention;

[0044] Figure 3 A comparison diagram of the logarithmic mean square error (MSE) of the target position and target velocity estimation by the method of the present invention (SDP, BR-SDP) and the Cramer-Rao lower bound (CRLB) as the measurement noise power changes;

[0045] Figure 4 It is a comparison diagram of the logarithmic deviation (Bias) of the target position and target speed estimation by the method of the present invention (SDP, BR-SDP) and the theoretical deviation (Theoretical) as the measurement noise power changes. DETAILED DESCRIPTION

[0046] The present invention is further described in detail below with reference to the accompanying drawings.

[0047] The present invention proposes a method for locating a uniformly moving target with an asynchronous transmitter and unknown position. The overall implementation block diagram is as follows: Figure 1 As shown, it includes the following steps:

[0048] Step 1: Figure 2 As shown in the figure, in a multi-base positioning scenario, there are several sensors and one target. Some sensors are used as receivers and their coordinate positions are known. Different receivers have the same reference clock, that is, the clocks of different receivers are synchronized. The remaining sensors are used as transmitters and their coordinate positions are unknown. The time between the transmitter and the receiver is asynchronous, that is, the reference clocks are different. There is a clock synchronization error between the transmitter and the receiver. The initial coordinate position of the target is unknown, and it moves at an unknown constant speed.

[0049] Specifically, in step 1, the multi-base positioning scenario is two-dimensional or three-dimensional, the number of receivers is N, and in this embodiment, N=6, and the coordinate position of the nth receiver is s n , n = 1, 2, ..., N, the number of transmitters is M, in this embodiment, M = 2, the coordinate position of the mth transmitter is The clock synchronization error of the mth transmitter is The initial coordinate position of the target is p o , the target's moving speed is v o .

[0050] Step 2: At different consecutive time points, the signal transmitted by the transmitter reaches the receiver via a direct path or an indirect path reflected by the target. The receiver can obtain the direct path distance measurement values ​​and the indirect path distance measurement values ​​observed at consecutive time points, and construct a distance measurement model based on this, including a direct path distance measurement part and an indirect path distance measurement part; wherein, the signal transmitted by the transmitter carries the clock information when the signal is transmitted.

[0051] It is worth noting that this positioning scenario is different from the existing positioning scenario. When the target is moving, it considers collecting only the delay measurement value. The present invention considers a new positioning scenario for the problem of mobile target positioning. Specifically: in traditional positioning methods, if it is necessary to achieve simultaneous estimation of the position and speed of a moving target, it is usually necessary to use delay information and Doppler information at the same time. The use of Doppler information requires that the signal acquisition equipment can achieve high-precision time-frequency measurement, and has the ability to sample and process signals at high speed. The method of the present invention innovatively uses only the delay information (i.e., the distance information measurement value r mn,k , d mn,k ), replacing the traditional method of mixing time delay and Doppler information. This significantly reduces the demand and complexity of system hardware, reduces dependence on precision equipment, and improves the feasibility of the system.

[0052] Further, the constructed distance measurement model is described as: Among them, K represents the total number of time points, r mn,k represents the indirect path distance measurement value obtained when the signal transmitted by the mth transmitter at time k reaches the nth receiver through the indirect path reflected by the target, d mn,k It represents the direct path distance measurement value obtained when the signal transmitted by the m-th transmitter at time k reaches the n-th receiver via the direct path. Specifically, after the signal transmitted by the m-th transmitter at time k reaches the n-th receiver via the direct path or the indirect path reflected by the target, the n-th receiver determines the signal flight time by comparing its own clock information when the signal arrives with the clock information carried by the signal, thereby obtaining the direct path distance measurement value or the indirect path distance measurement value. Represents r mn,k The corresponding measurement noise is, Indicates d mn,k The corresponding measurement noise, p o +kv o represents the coordinate position of the target k at time k, "||||" is the Euclidean norm operator, used to represent the distance, ||p o +kv o -s n || represents the distance between the kth target and the nth receiver at time instant, represents the distance between the target k and the mth transmitter at time instant, Represents the distance between the mth transmitter and the nth receiver.

[0053] Step 3: For the indirect path distance measurement part of the distance measurement model, the indirect path distance measurement values ​​obtained by the remaining receivers except the first receiver are subtracted from the indirect path distance measurement value obtained by the first receiver, and the unknown variables are separated through mathematical processing; for the direct path distance measurement part of the distance measurement model, the unknown variables are separated through mathematical processing; then all the unknown variables are used to transform the distance measurement model after the separation of the unknown variables into the original optimization problem, which is a weighted least squares problem with constraints.

[0054] Specifically, in step 3, the specific process of separating the unknown variables in the indirect path distance measurement part of the distance measurement model is:

[0055] Step 3.1a: Add the indirect path distance measurement part of the distance measurement model

[0056]

[0057] Split into two types, namely

[0058] and

[0059]

[0060] Step 3.2a: Subtract the two equations obtained in step 3.1a to obtain This difference step effectively reduces the number of coupled variables, achieving the purpose of tightening the final optimization problem. The present invention cleverly utilizes the numerical relationship between the target reflection path distance information measurement values ​​between different receivers, and eliminates the common norm term before forming the optimization problem. This process has the following advantages: (1) By eliminating the common norm term, the number of coupled variables in the final optimization problem is effectively reduced, thereby reducing the complexity of the calculation. The reduction in the number of variables means that when solving the optimization problem, the calculation time and resource consumption can be shortened, and the solution efficiency can be improved. (2) Removing the common norm term can make the optimization problem more compact and the constraints more stringent. This has the effect of "tightening" the optimization problem to a certain extent, making the final solution more accurate and stable. Such a tightening effect helps to improve the convergence of the optimization problem and reduce the feasible domain of the solution, thereby enhancing the robustness and accuracy of the algorithm.

[0061] Step 3.3a: For r in step 3.1a m1,k In the expression Perform mathematical processing of transposing terms, and then perform mathematical processing of square expansion on the result after transposing terms, ignoring the second-order noise term get: The unknown variables are separated, where The superscript "T" indicates the transpose operation of a matrix or vector. 1 、r m1,k and Respectively represent s n 、r mn,k and The result when n = 1; for r in step 3.2a mn,k -r m1,k ||p in the expression o +kv o -s 1 ||Perform mathematical processing of transposing terms, and then perform mathematical processing of square expansion on the result after transposing terms to obtain: The unknown variables are separated, where

[0062] Specifically, in step 3, the specific process of separating the unknown variables in the direct path distance measurement part of the distance measurement model is as follows:

[0063] Direct path distance measurement part of the distance measurement model In Perform mathematical processing of transposing the terms, and then perform mathematical processing of square expansion on the result after transposing the terms, ignoring the second-order noise term, and obtain: The unknown variables are separated, where

[0064] Specifically, in step 3, the process of obtaining the original optimization problem is as follows:

[0065] Step 3.1b: Combine equation 1, equation 2, and equation 3 to obtain the distance measurement model after separating the unknown variables, which is described as: Then all unknown variables are expressed as vectors, denoted as x o , in, are coupled variables.

[0066] Step 3.2b: Utilize x o The distance measurement model after separating the unknown variables is transformed into the original optimization problem. The original optimization problem is a weighted least squares problem with constraints, which is described as: Among them, min means to minimize the objective function, and st means "subject to the constraints...". Consider x o The constraints generated by the numerical relationship between the unknown variables in the optimization problem are: x is the optimization variable in the original optimization problem, and x and x o correspond, p represents the target position variable, v represents the target velocity variable, t 1 represents the first transmitter position variable, t M represents the Mth transmitter position variable, τ 1 represents the clock synchronization error variable of the first transmitter, τ M represents the clock synchronization error variable of the Mth transmitter, ξ 1,1 ,...,ξ 1,K , They represent the coupled variables formed by different combinations of the previous independent variables, (g-Hx) T W(g-Hx) is the objective function in the original optimization problem, L represents the dimension of the multi-base positioning scenario, L is 2 or 3, and x is defined i represents the i-th element of x, x ((M+2)L+M+k) Represents x i The result when i=(M+2)L+M+k, x ((M+2)L+K+kM+m) Represents x i The result when i=(M+2)L+K+kM+m, x ((M+2)L+m)Represents x i The result when i=(M+2)L+m, x ((M+2)L+(M+1)K+M+m) Represents x i The result when i=(M+2)L+(M+1)K+M+m, define x (i:j) represents the column vector consisting of the i-th element to the j-th element of x, x (1:L) Represents x (i:j) The result when i=1, j=L, x (L+1：2L) Represents x (i:j) The result when i=L+1,j=2L, x ((m+1)L+1:(m+2)L) Represents x (i:j) The result when i=(m+1)L+1,j=(m+2)L is: g, g r1,k , g r2,k and d,k are the introduced coefficient vectors, s 2 Indicates n The result when n=2, s N Indicates n The result when n = N, r 11,k Represents r mn,k The result when n=1, m=1, r M1,k Represents r mn,k The result when n=1,m=M, r 12,k Represents r mn,k The result when n=2, m=1, r M2,k Represents r mn,k The result when n=2, m=M, r 1N,k Represents r mn,k The result when n=N,m=1, r MN,k Represents r mn,k The result when n=N,m=M, d 11,k Indicates d mn,k The result when n=1, m=1, d M1,k Indicates d mn,k The result when n=1,m=M, d 1N,k Indicates d mn,k The result when n=N,m=1, d MN,k Indicates d mn,k The result when n=N,m=M is: B=blkdiag(B 1 ,...,B K ),B k =blkdiag(B r1,k ,B r2,k ,B d ), B.B k , B r1,k , B r2,k , B d are the introduced coefficient matrices, and I is defined i represents the identity matrix of dimension i×i, I MN Indicates I i The result when i = MN, I M Indicates I i The result when i=M, Definition 1 i represents a vector of all 1s of length i, 1 N-1 Representation 1 i The result when i=N-1, define 0 i×j represents a matrix of all zeros with dimension i×j, 0 M×MN Represents 0 i×j The result when i=M, j=MN is 0 M(N-1)×M(N-1) Represents 0 i×j The result when i=M(N-1), j=M(N-1), Diag() represents the element diagonalization operation function, that is, the vectors are arranged diagonally according to each element to form a matrix, blkdiag() represents the matrix diagonalization operation function, that is, each matrix is ​​arranged diagonally to form a matrix, represents the distance between the target and the first transmitter at time k, represents the distance between the target k and the Mth transmitter at time instant, Indicates the distance between the first transmitter and the first receiver, represents the distance between the Mth transmitter and the first receiver, represents the distance between the 1st transmitter and the Nth receiver, represents the distance between the Mth transmitter and the Nth receiver, represents the distance between the target and the second receiver at time k, represents the distance between the target k and the Nth receiver at time k, Q k Represents ε k The corresponding covariance matrix, ε k represents the composite noise vector at time k, Respectively represent r 11,k 、r M1,k 、r 1N,k 、r MN,k d 11,k d M1,k d 1N,k d MN,k The corresponding measurement noise, ε kSubject to the mean of 0 and the covariance matrix is ​​Q k Gaussian distribution, H.H k , H r1,k , H r1,k1 , H r1,k2 , H r1,k3 , H r2,k , H r2,k1 , H r2,k2 , H d,k , H d,k1 , H d,k2 , H d,k3 are the introduced coefficient matrices, 1 M Representation 1 i The result when i=M is ​​1 N Representation 1 i The result when i=N is 0 M×ML , 0 M×M , 0 M×(k-1) , 0 M×(K-k) , 0 M×(k-1)M , 0 M×(K-k)M , 0 M(N-1)×M(L+1) , 0 M(N-1)×M(K+1) , 0 M(N-1)×(k-1) , 0 M(N-1)×(K-k) , 0 MN×2L and 0 MN×K(M+1) Respectively represent 0 i×j When i=M,j=ML,i=M,j=M,i=M,j=k-1,i=M,j=Kk,i=M,j=M(k-1),i=M,j=M(Kk),i=M(N-1),j=M(L+ 1), the result when i=M(N-1),j=M(K+1),i=M(N-1),j=k-1,i=M(N-1),j=Kk,i=MN,j=2L,i=MN,j=K(M+1), r n,k =[r 1n,k ,...,r Mn,k ] T , r 1n,k Represents r mn,k The result when m=1, r Mn,k Represents r mn,k The result when m=M, d n,k =[d 1n,k ,...,d Mn,k ] T , d 1n,k Indicates d mn,k The result when m=1, d Mn,k Indicates d mn,k The result when m=M is, Represented by the matrix The submatrix consisting of the M+1th row to the MNth row and all columns, is the Kronecker product operator.

[0067] Step 4: Perform a first-order Taylor expansion on the objective function in the original optimization problem; then obtain a constraint on the bias term based on the first-order Taylor expansion; then combine the original optimization problem and the constraint on the bias term to obtain a constrained weighted least squares problem with bias reduction capability.

[0068] Specifically, the specific process of step 4 is:

[0069] Step 4.1: Define optimization variables in, α is an intermediate variable; then the objective function (g-Hx) in the original optimization problem is T W(g-Hx) is transformed into y T R T WRy, where R represents the coefficient matrix carrying noise, R = [-H, g]; considering that the measurement values ​​in the coefficient matrix R carry measurement errors, the objective function y T R T WRy is expanded by the first order Taylor to obtain the first order Taylor expansion, which is described as: in, It means to find the mathematical expectation, R o Represents the coefficient matrix without noise, r mn,k The corresponding indirect path distance measurement without noise d mn,k The corresponding direct path distance measurement without noise Substitute r in R = [-H, g] mn,k d mn,k Get R o , “⊙” is the Hadamard product operator, ΔR=[-ΔH,Δg], A k,r1 =[0 M×(M+2)L ,I M ,0 M×(k-1) ,1 M ,0 M×(K-k) ,0 M×M(K+1) ], A k,r2 =[0 M(N-1)×((M+2)L+(M+k-1)) ,-1 M(N-1) ,0 M(N-1)×(K-k) ,0 M(N-1)×M(K+1) ], represents the composite noise vector after deformation at K moments, represents the composite noise vector after deformation at time k, express The result when k=1 is, express The result when k = K is, express The result when m=1 is, express The result when m=M is, express The result when m=1 is, express Results when m=M.

[0070] Step 4.2: Based on the first-order Taylor expansion, introduce Form a constraint on the deviation term, described as: in,

[0071] Step 4.3: Combining the original optimization problem and the constraint on the bias term, we get a constrained weighted least squares problem with bias reduction capability, which is described as: Among them, y end Represents the last element of y, defining y (i:j) represents the column vector consisting of the i-th element to the j-th element of y, y (1:L) Represents y (i:j) The result when i=1,j=L, y ((L+1):2L ) represents y (i:j) The result when i=L+1,j=2L, y (((m+1)L+1):(m+2)L ) represents y (i:j) The result when i=(m+1)L+1,j=(m+2)L is defined as y i represents the i-th element of y, y ((M+2)L+M+k) Represents y i The result when i=(M+2)L+M+k is y ((M+2)L+K+kM+m) Represents y i The result when i=(M+2)L+K+kM+m, y ((M+2)L+m) Represents y i The result when i=(M+2)L+m, y ((M+2)L+(M+1)K+M+m) Represents y i The result when i=(M+2)L+(M+1)K+M+m.

[0072] Step 5: Use the semi-positive definite relaxation technique to relax the constrained weighted least squares problem with bias reduction capability into a convex semi-positive definite programming problem that is easy to solve.

[0073] Specifically, in step 5, the convex semi-positive definite programming problem is described as: Where Y represents the newly introduced optimization variable, and Y=yy T ,tr{R T WRY} is the objective function in the convex semi-positive definite programming problem, tr{·} means taking the diagonal elements of the matrix and summing them, Y is the optimization variable of the convex semi-positive definite programming problem, and is defined as It means from the i-th row to the j-th row and the l-th column to the The submatrix Y (1:L,1:L) express When i=1, j=L, l=1, The result, Y ((L+1):2L,(L+1):2L) express When i=L+1,j=2L,l=L+1, The result, Y (1:L,(L+1):2L) express When i=1, j=L, l=L+1, The result, Y (1:L,end) express When i=1, j=L, The result, Y ((L+1):2L,end) express When i=L+1,j=2L, The result, Y (((m+1)L+1):(m+2)L,((m+1)L+1):(m+2)L )express When i=(m+1)L+1, j=(m+2)L, l=(m+1)L+1, The result, Y (1:L,((m+1)L+1):(m+2)L) express When i=1, j=L, l=(m+1)L+1, The result, Y ((L+1):2L,((m+1)L+1):(m+2)L) express When i=L+1,j=2L,l=(m+1)L+1, The result when Y (i,j) represents the element corresponding to the i-th row and j-th column of Y, Y ((M+2)L+m,(M+2)L+m) Represents Y (i,j) The result when i=(M+2)L+m, j=(M+2)L+m, Y ((M+2)L+M+k,(M+2)L+M+k) Represents Y (i,j) The result when i=(M+2)L+M+k, j=(M+2)L+M+k is Y (end,end) Represents Y (i,j) The result when i=j=end, Y ((M+2)L+K+kM+m,end) Represents Y (i,j) The result when i=(M+2)L+K+kM+m,j=end, Y ((M+2)L+m,(M+2)L+M+k) Represents Y (i,j) The result when i=(M+2)L+m, j=(M+2)L+M+k is Y((M+2)L+(M+1)K+M+m,end) Represents Y (i,j) The result when i=(M+2)L+(M+1)K+M+m,j=end, Y (end,end) is the index of the last row or column of Y. Y ≥ 0 means that the matrix Y is positive semidefinite.

[0074] Step 6: Use existing toolboxes such as the interior point method to directly solve the convex semi-positive programming problem to obtain the optimal estimate of the initial coordinate position of the target, the optimal estimate of the target's moving speed, the optimal estimate of the coordinate position of each transmitter, and the optimal estimate of the clock synchronization error of each transmitter.

[0075] Specifically, in step 6, the initial coordinate position p of the target is o The optimal estimate of * , Set the target's moving speed v o The optimal estimate of * , The coordinate position of the mth transmitter The optimal estimate of The clock synchronization error of the mth transmitter is The optimal estimate of in, Y * represents the optimal estimate of Y, and defines represents the submatrix consisting of the elements corresponding to the i-th row to the j-th row and the last column of Y, express The result when i=1,j=L is: express The result when i=L+1,j=2L is: express The result when i=(m+1)L+1,j=(m+2)L is defined as represents the element corresponding to the i-th row and j-th column of Y, express The result when i=(M+2)L+m,j=end, express The result when i=j=end.

[0076] In order to verify the feasibility and effectiveness of the method of the present invention, a simulation test was carried out on the method of the present invention.

[0077] Assume that there are 2 transmitters (i.e., M=2) and 6 receivers (i.e., N=6) in a three-dimensional space, i.e., a three-dimensional coordinate system (i.e., L=3). Their real Cartesian coordinate positions in the three-dimensional space are randomly generated, and the real Cartesian coordinate position (initial coordinate position) and real moving speed of the target are also randomly generated. The range is shown in Table 1:

[0078] Table 1 Scenario range

[0079]

[0080] At the same time, the scene and the measurement values ​​are scaled by 10 times (to avoid calculation problems in the toolbox due to excessively large measurement values ​​during the problem solving process). The measurement noise is assumed to have a mean of 0 and a covariance matrix Q at time k. k =σ 2 I 2MN Gaussian noise, where σ 2 represents the variance of the measurement noise, I 2MN Represents the identity matrix of dimension 2MN×2MN.

[0081] Based on the above parameter settings, 10 test scenes were randomly generated in the simulation, and the optimal estimated value of the initial coordinate position of the target, the optimal estimated value of the moving speed of the target and the logarithmic mean square error (MSE) and logarithmic bias (Bias) of the method of the present invention were tested in each test scene. To verify the stability of the method, each test scene was subjected to 1000 Monte Carlos, and the final result was the average of 1000 Monte Carlos and 10 test scenes. To verify the performance of the method, the Cramer-Rao bound (CRLB) and theoretical bias (Theoretical) were introduced for comparison.

[0082] Figure 3 A comparison diagram of the logarithmic mean square error (MSE) of the target position and target velocity estimation by the method of the present invention (SDP, BR-SDP) and the Cramer-Rao lower bound (CRLB) as the measurement noise power changes is shown. The upper figure is a schematic diagram of the logarithmic mean square error (MSE) of the target position estimation as the measurement noise power changes, and the lower figure is a schematic diagram of the logarithmic mean square error (MSE) of the target velocity estimation as the measurement noise power changes; Figure 4 A comparison chart of the logarithmic deviation (Bias) of target position and target speed estimation using the method of the present invention (SDP, BR-SDP) and the theoretical deviation (Theoretical) is shown. The upper figure is a schematic diagram of the logarithmic deviation (Bias) of target position estimation changing with the measured noise power, and the lower figure is a schematic diagram of the logarithmic deviation (Bias) of target speed estimation changing with the measured noise power. Figure 3 and Figure 4In the above, SDP refers to a semi-definite relaxation method, i.e., a method obtained by the method of the present invention without the process of reducing the deviation capability, and BR-SDP refers to a semi-definite relaxation method with the ability to reduce the deviation. Figure 3 As can be seen, when the measured noise power increases by 10log 10 σ 2 When changes, the method of the present invention can approach the Cramer-Rao lower bound under the entire noise condition, indicating that the method of the present invention can accurately estimate the target position and target speed using only the time delay information. At the same time, the elimination of the coupling variables does not affect the performance of the final estimation result and also reduces the computational complexity of the problem. Figure 4 It shows that the method of the present invention can effectively reduce the logarithmic deviation of the final estimation result by reducing the deviation, which has an advantage of about 15dB compared with the result without deviation reduction, and can be close to the theoretical deviation when the noise is medium or below medium.

[0083] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principle of the present invention should be included in the protection scope of the present invention.

Claims

1. A method for locating a uniformly moving target with an asynchronous transmitter and unknown position, characterized in that The following steps are involved: Step 1: In a multi-base positioning scenario, assume that there are several sensors and a target. Some sensors are receivers with known coordinate positions. Different receivers have the same reference clock, that is, the clocks of different receivers are synchronized. The remaining sensors are transmitters with unknown coordinate positions. The transmitter and receiver are asynchronous in time, that is, the reference clocks are different. There is a clock synchronization error between the transmitter and the receiver. The initial coordinate position of the target is unknown and it moves at an unknown constant speed. Step 2: At different time points, the signal transmitted by the transmitter reaches the receiver via a direct path or an indirect path reflected by the target. The receiver can obtain the direct path distance measurement value and the indirect path distance measurement value observed at the continuous time points, and accordingly construct a distance measurement model, including a direct path distance measurement part and an indirect path distance measurement part; wherein, the signal transmitted by the transmitter carries the clock information when the signal is transmitted; Step 3: For the indirect path distance measurement part of the distance measurement model, the indirect path distance measurement values ​​obtained by the remaining receivers except the first receiver are subtracted from the indirect path distance measurement value obtained by the first receiver, and the unknown variables are separated through mathematical processing; for the direct path distance measurement part of the distance measurement model, the unknown variables are separated through mathematical processing; then, all unknown variables are used to transform the distance measurement model after the unknown variables are separated into the original optimization problem, which is a weighted least squares problem with constraints; Step 4: Perform a first-order Taylor expansion on the objective function in the original optimization problem; then obtain a constraint on the deviation term based on the first-order Taylor expansion; then combine the original optimization problem and the constraint on the deviation term to obtain a constrained weighted least squares problem with the ability to reduce deviation; Step 5: Use the semi-positive definite relaxation technique to relax the constrained weighted least squares problem with bias reduction capability into a convex semi-positive definite programming problem that is easy to solve; Step 6: Directly solve the convex semi-positive programming problem to obtain the optimal estimate of the initial coordinate position of the target, the optimal estimate of the target's moving speed, the optimal estimate of the coordinate position of each transmitter, and the optimal estimate of the clock synchronization error of each transmitter.

2. The method for locating a uniformly moving target with an asynchronous transmitter and unknown position according to claim 1 is characterized in that In step 1, the multi-base positioning scenario is two-dimensional or three-dimensional, the number of receivers is N, and the coordinate position of the nth receiver is s n , n=1,2,…,N, the number of transmitters is M, and the coordinate position of the mth transmitter is The clock synchronization error of the mth transmitter is The initial coordinate position of the target is p o , the target's moving speed is v o .

3. The method for locating a uniformly moving target with an asynchronous transmitter and unknown position according to claim 2 is characterized in that In step 2, the constructed distance measurement model is described as: ， Among them, K represents the total number of time points, r mn,k represents the indirect path distance measurement value obtained when the signal transmitted by the mth transmitter at time k reaches the nth receiver through the indirect path reflected by the target, d mn,k represents the direct path distance measurement value obtained when the signal transmitted by the mth transmitter at time k reaches the nth receiver through the direct path, Represents r mn,k The corresponding measurement noise is, Indicates d mn,k The corresponding measurement noise, p o +kv o represents the coordinate position of the target at time k, and "||| |" is the Euclidean norm operator.

4. The method for locating a uniformly moving target with an asynchronous transmitter and unknown position according to claim 3 is characterized in that In step 3, the specific process of separating the unknown variables in the indirect path distance measurement part of the distance measurement model is as follows: Step 3.1a: Add the indirect path distance measurement part of the distance measurement model Split into two types, namely and Step 3.2a: Subtract the two equations obtained in step 3.1a to obtain Step 3.3a: For r in step 3.1a m1,k In the expression Perform mathematical processing of transposing the terms, and then perform mathematical processing of square expansion on the result after transposing the terms, ignoring the second-order noise term, and obtain: (Formula 1) realizes the separation of unknown variables, where The superscript "T" indicates the transpose operation of a matrix or vector; for r in step 3.2a mn,k -r m1,k ||p in the expression o +kv o -s1|| performs mathematical processing of transposing the term, and then performs mathematical processing of square expansion on the result after transposing the term, and obtains: (Formula 2) realizes the separation of unknown variables, where In step 3, the specific process of separating the unknown variables in the direct path distance measurement part of the distance measurement model is as follows: Direct path distance measurement part of the distance measurement model In Perform mathematical processing of transposing the terms, and then perform mathematical processing of square expansion on the result after transposing the terms, ignoring the second-order noise term, and obtain: (Equation 3) realizes the separation of unknown variables, where 5. The method for locating a uniformly moving target with an asynchronous transmitter and unknown position according to claim 4 is characterized in that In step 3, the process of obtaining the original optimization problem is as follows: Step 3.1b: Combine equation 1, equation 2, and equation 3 to obtain the distance measurement model after separating the unknown variables, which is described as: ； Then all unknown variables are expressed as vectors, denoted as x o , in, All are coupled variables; Step 3.2b: Utilize x o The distance measurement model after separating the unknown variables is transformed into the original optimization problem. The original optimization problem is a weighted least squares problem with constraints, which is described as: Among them, x is the optimization variable in the original optimization problem, x and x o Corresponding, (g-Hx) T W(g-Hx) is the objective function in the original optimization problem, L represents the dimension of the multi-base positioning scenario, L is 2 or 3, and x is defined i Represents the i-th element of x, and defines x (i:j) represents the column vector consisting of the i-th element to the j-th element of x, B=blkdiag(B1,...,B K ),B k =blkdiag(B r1,k ,B r2,k ,B d ), ， Definition I i represents the identity matrix of dimension i×i, Definition 1 i represents a vector of all 1s of length i, and defines 0 i×j represents an all-zero matrix of dimension i×j, Diag( ) represents the element diagonalization operation function, that is, the vectors are arranged diagonally to form a matrix, blkdiag( ) represents the matrix diagonalization operation function, that is, each matrix is ​​arranged diagonally to form a matrix, Q k Represents ε k The corresponding covariance matrix, ε k represents the composite noise vector at time k, ε k Subject to the mean of 0 and the covariance matrix is ​​Q k Gaussian distribution, r n,k =[r 1n,k ,...,r Mn,k ] T , d n,k =[d 1n,k ,...,d Mn,k ] T , Represented by moment Array The submatrix consisting of the M+1th row to the MNth row and all columns, is the Kronecker product operator.

6. The method for locating a uniformly moving target with an asynchronous transmitter and unknown position according to claim 5 is characterized in that The specific process of step 4 is as follows: Step 4.1: Define optimization variables in, α is an intermediate variable; then the objective function (g-Hx) in the original optimization problem is T W(g-Hx) is transformed into y T R T WRy, where R represents the coefficient matrix carrying noise, R = [-H, g]; then the objective function y T R T WRy is expanded by the first order Taylor to obtain the first order Taylor expansion, which is described as: ， in, represents the mathematical expectation, R o Represents the coefficient matrix without noise, r mn,k The corresponding indirect path distance measurement without noise d mn,k The corresponding direct path distance measurement without noise is Substitute r in R = [-H, g] mn,k ,d mn,k Get R o , "⊙" is the Hadamard product operator, ΔR = [-ΔH, Δg], A k,r1 =[0 M×(M+2)L ,I M ,0 M×(k-1) ,1 M ,0 M×(K-k) ,0 M×M(K+1) ], A k,r2 ＝[0 M(N-1)×((M+2)L+(M+k-1)) ,-1 M(N-1) ,0 M(N-1)×(K-k) ,0 M(N-1)×M(K+1) ]， represents the composite noise vector after deformation at K moments, represents the composite noise vector after deformation at time k, Step 4.2: Based on the first-order Taylor expansion, introduce Form a constraint on the deviation term, described as: in, Step 4.3: Combining the original optimization problem and the constraint on the bias term, we get a constrained weighted least squares problem with bias reduction capability, which is described as: Among them, y end Represents the last element of y, defining y (i:j) represents the column vector consisting of the i-th element to the j-th element of y, and defines y i Represents the i-th element of y.

7. The method for locating a uniformly moving target with an asynchronous transmitter and unknown position according to claim 6 is characterized in that In step 5, the convex semi-positive definite programming problem is described as: AND ((M+2)L+K+kM+m,end) =-tr{Y (((m+1)L+1):(m+2)L,((m+1)L+1):(m+2)L) }+2tr{And (1:L,((m+1)L+1):(m+2)L) }+2ktr{And ((L+1):2L,((m+1)L+1):(m+2)L) }+Y ((M+2)L+m,(M+2)L+m) +2Y ((M+2)L+m,(M+2)L+M+k) AND ((M+2)L+(M+1)K+M+m,end) =And ((M+2)L+m,(M+2)L+m) -tr{And (((m+1)L+1):(m+2)L,((m+1)L+1):(m+2)L) } ， Where Y represents the newly introduced optimization variable, and Y=yy T ,tr{R T WRY} is the objective function in the convex semi-positive definite programming problem, tr{·} means taking the diagonal elements of the matrix and summing them, Y is the optimization variable of the convex semi-positive definite programming problem, and is defined as It means from the i-th row to the j-th row and the l-th column to the The submatrix consisting of the elements corresponding to the columns is defined as Y (i,j) represents the element corresponding to the i-th row and j-th column of Y, Y (end,end) is the index of the last row or column of Y. Y ≥ 0 means that the matrix Y is positive semidefinite.

8. The method for locating a uniformly moving target with an asynchronous transmitter and unknown position according to claim 7 is characterized in that In step 6, the initial coordinate position p of the target is o The optimal estimate of * , Set the target's moving speed v o The optimal estimate of * , The coordinate position of the mth transmitter The optimal estimate of The clock synchronization error of the mth transmitter is The optimal estimate of in, Y * represents the optimal estimate of Y, and defines Represents the submatrix consisting of the elements corresponding to the i-th row to the j-th row and the last column of Y, and defines Represents the element corresponding to the i-th row and j-th column of Y.

Images