CN105259533A - Three-stage arrival time difference positioning method based on multidimensional scaling sub space analysis - Google Patents

Abstract

A three-stage time difference of arrival positioning method based on multidimensional scale subspace analysis in the field of wireless networks and mobile computing, including three stages, the first stage: obtaining the coordinates of the sensor position on the plane, the time difference of arrival, and the variance of the measurement error of the time difference of arrival and signal propagation speed, calculate the corrected measurement scalar product matrix, through eigenvalue decomposition, calculate the first estimated point and a distance estimated value, and calculate a direction vector; the second stage: on the basis of the first estimated point, Using the estimated distance value and direction vector, calculate the parameters corresponding to the second estimated point; the third stage: use the parameters corresponding to the second estimated point as the initial value, and obtain the corresponding parameters, and further calculate the third estimated point to obtain the final estimated value of the coordinates, thereby determining the position of the signal source. The invention has the characteristics of wide application range and high positioning precision and accuracy.

Description

Three-stage Time Difference of Arrival Localization Method Based on Subspace Analysis of Multidimensional Scaling Method

technical field

The invention relates to the field of wireless network and mobile computing, in particular to a three-stage time difference of arrival positioning method based on multi-dimensional scale subspace analysis.

Background technique

In the application fields of radar, sonar, mobile communication, multimedia, wireless sensor network, etc., an important problem is often faced, that is, to locate a signal source according to the time difference of arrival information. The so-called time difference of arrival refers to: a signal is sent by a signal source, and the signal is received by a sensor distributed in space, whose position is known and time-synchronized with each other, and the time when the signal arrives at other sensors is measured, and the signal source is calculated from this. The difference between the time a signal arrives at each sensor and the time it arrives at the first sensor.

After searching the existing technologies, it was found that He‐WenWei, RongPeng, QunWan, Zhang‐XinChen, & Shang‐FuYe published the paper "Multidimensional scaling analysis for passive moving target localization with TDOAandFDOAmeasurements" in the journal IEEE Transactions on Signal Processing, Volume 58, Issue 3, March 2010. A method of analyzing the arrival time difference and arrival frequency difference of the passive sensor array to obtain the position of the signal source. However, in this method, when the sensor array becomes linear or nearly linear, the matrix that needs to be inverted will become a rank deficient matrix or a highly ill-conditioned matrix, and forcibly inverting it will cause a huge positioning error.

Chinese Patent Document No. CN103648164A, publication (announcement) date 2014.03.19, discloses a wireless sensor network distributed positioning method based on time difference of arrival and Gossip algorithm. Its characteristics are: 1. The anchor node obtains its own position coordinates; 2. Realizes distributed time synchronization; 3. The anchor node randomly wakes up and monitors unknown nodes; 4. Wakes up the anchor node to save the received signal time and local coordinates; Complete signal monitoring and data storage; 6. Anchor node j receives the data of all M adjacent anchor nodes; 7. Obtain the initial estimated value of anchor node j’s unknown unknown node; 8. All anchor nodes obtain the location of unknown nodes Initial estimated value; 9. Run the Gossip algorithm to randomly select adjacent anchor nodes to exchange positioning data; 10. Terminate the algorithm. However, in this method, the unknown node position estimates collected by different anchor nodes may not be completely consistent, and how to calculate the most likely position estimate from the obtained position information is not mentioned in the application.

Contents of the invention

Aiming at the above-mentioned deficiencies in the prior art, the present invention proposes a three-stage time difference of arrival positioning method based on multi-dimensional scale subspace analysis, and estimates the position of the signal source in three stages by performing subspace analysis on the measurement scalar product matrix , to obtain a more accurate position of the signal source.

The present invention is achieved through the following technical solutions,

The present invention comprises the following steps:

Step 1. Obtain the measurement scalar product matrix according to the position coordinates of the sensors on the plane and the time difference of arrival (i.e. the time difference between the time when the signal sent by the signal source reaches each sensor and the time when the signal arrives at the first sensor), and obtain the measured scalar product matrix through the time difference of arrival time The variance of the measurement error corrects the measurement scalar product matrix; by performing eigenvalue decomposition processing on the corrected measurement scalar product matrix, the first estimated point is calculated, and a direction vector and a distance estimate are calculated, specifically including the following step:

Step 1.1), generate measurement scalar product matrix The sensor position coordinates on the plane are u _m =[x _m ,y _m ] ^T , m=1,...,M, where: M represents the number of sensors and the number is greater than or equal to 5, and u _m represents the mth The position coordinates of the sensor, x _m represents the x-axis coordinate of the mth sensor, and y _m represents the y-axis coordinate of the mth sensor;

According to the arrival time difference and the signal propagation speed, the arrival distance difference is obtained: when m=1 When m=2,...,M in: Indicates the difference between the measured arrival time from signal source u ₀ to each sensor u _m and the arrival time from signal source u ₀ to the first sensor u ₁ , c indicates the signal propagation speed;

When m=2,...,M, according to The variance of the measurement error and the signal propagation speed get error variance of in: express The measurement error variance of , c represents the signal propagation speed;

The measured scalar product matrix is where: matrix The elements of row i and column j of are

${<span\; class="notranslate">\; \&lsqb;</span>\stackrel{<span\; class="notranslate">\; ^</span>}{\mathrm{<span\; class="notranslate">\; B</span>}}<span\; class="notranslate">\; \&rsqb;</span>}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; j</span>}}<span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; 2</span>}}{<span\; class="notranslate">\; (</span>{\stackrel{<span\; class="notranslate">\; ^</span>}{\mathrm{<span\; class="notranslate">\; d</span>}}}_{\mathrm{<span\; class="notranslate">\; i</span>}\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; -</span>{\stackrel{<span\; class="notranslate">\; ^</span>}{\mathrm{<span\; class="notranslate">\; d</span>}}}_{\mathrm{<span\; class="notranslate">\; j</span>}\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; -</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; \&lsqb;</span>{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; j</span>}}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; +</span>{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; j</span>}}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; \&rsqb;</span><span\; class="notranslate">\; ;</span>$

Step 1.2), correcting the measured scalar product matrix: obtain the corrected measured scalar product matrix Where: I _M represents the M×M unit matrix, and 1 _M represents the M-dimensional column vector whose elements are all 1;

Step 1.3), calculating the first estimated point and the estimated distance, specifically include the following steps:

Step 1.3.1), perform eigenvalue decomposition on B ₁ -σ ² I _M :

B ₁ -σ ² I _M ＝[v ₁ ,...,v _M ]diag(s ₁ ,...,s _M )[v ₁ ,...,v _M ] ^T , where: v ₁ ,. ..,v _M ∈ R ^M is a pairwise orthogonal vector with a modulus of 1. The eigenvalue matrix diag(s ₁ ,...,s _M ) indicates that the diagonal elements are s ₁ ,...,s The diagonal matrix of _M , the eigenvalues s ₁ ,...,s _M are arranged in descending order of absolute value, that is, |s ₁ |≥...≥|s _M |;

Step 1.3.2), the coefficient vector is obtained by the linear combination of the eigenvectors: v ₁ ,...,v _M is a set of orthonormal basis of ^RM , so the coefficient vector v in ^RM can be expressed as v ₁ , ..., the linear combination of v _M v=k ₁ v ₁ +...+k ₅ v ₅ , when m=1,...,5, the combination coefficient of base v _m And when m≥6, the combination coefficients k ₆ ,...,k _M of base v _m are all 0;

Step 1.3.3), with the coefficient vector v as the weighting coefficient, the position coordinate matrix $\left[\begin{array}{ccc}{x}_1& ...& x_m\\ {\mathrm{the\; y}}_1& ...& {\mathrm{the\; y}}_m\end{array}\right]$ The column vectors of are linearly combined to calculate the first estimated point $u^{\left(1\right)}=\frac{\left[\begin{array}{ccc}{x}_1& ...& x_m\\ {\mathrm{the\; y}}_1& ...& {\mathrm{the\; y}}_m\end{array}\right]v}{1_m^{T}v}$ and an estimate of the distance from the estimated point $d_1^{\left(1\right)}=-\frac{\left[\begin{array}{ccc}{\hat{d}}_{11}& ...& {\hat{d}}_{m1}\end{array}\right]v}{1_m^{T}v};$

Step 1.4), calculate the direction vector p: when but otherwise $p=\frac{{\&lsqb;-({\mathrm{the\; y}}_2-{\mathrm{the\; y}}_1),x_2-x_1\&rsqb;}^T}{\left|\right|u_2-u_1\left|\right|},$ in: $u_f=\frac{\left[\begin{array}{ccc}{x}_1& ...& x_m\\ {\mathrm{the\; y}}_1& ...& {\mathrm{the\; y}}_m\end{array}\right]v_5}{1_m^T{v}_5},$ ||·|| represents the Euclidean norm of the vector;

Step 2, on the basis of the first estimated point, using the estimated distance and the direction vector p, calculate the intermediate parameter t ₂ corresponding to the second estimated point:

definition $t_2=t_2^{\left(\mathrm{no}\right)},a_0=\left|\right|u^{\left(1\right)}-u_1|{|}^2-{\left(d_1^{\left(1\right)}\right)}^2,$ a ₁ =p ^T (u ⁽¹⁾ -u ₁ ), $\mathrm{\&Delta;}=a_1^2-a_0;$ Then when Δ>0, n=1, define Define when n=2 Otherwise, define when n=1 When n=2 undefined;

Step 3. Using the parameter corresponding to the second estimated point as the initial value, obtain the parameter corresponding to the third estimated point through the dichotomy root-finding process, further calculate the third estimated point, and obtain the final estimated value of the coordinate, thereby determining the signal source location, the specific steps include:

Step 3.1), setting the parameter t2 corresponding to the second estimation point u ⁽²⁾ _as the initial value: define f _min =∞, n=1,

Step 3.2), find the root g(t)=0 by dichotomy, find out the parameter t ₃ corresponding to the third estimation point u ⁽³⁾ , record the components of u(t) as x _u , y _u , and the components of p Recorded as l _x ,l _y , namely: u(t)=[x _u ,y _u ] ^T , p=[l _x ,l _y ] ^T , for m=1,...,M, set w _m = [x _u -x _m ,y _u -y _m ] ^T ，d _m ＝||w _m ||， ${\stackrel{\&CenterDot;}{d}}_m=p^T{w}_m/\left|\right|w_m\left|\right|,$ further order $Z=\left[\begin{array}{ccc}{x}_1-x_u& ...& x_m-x_u\\ {\mathrm{the\; y}}_1-{\mathrm{the\; y}}_u& ...& {\mathrm{the\; y}}_m-{\mathrm{the\; y}}_u\\ d_1& ...& d_m\end{array}\right],$ B=Z ^T diag(1,1,-1) Z, $\stackrel{\&CenterDot;}{Z}=\left[\begin{array}{ccc}-l_x& ...& -l_x\\ -l_{\mathrm{the\; y}}& ...& -l_{\mathrm{the\; y}}\\ {\stackrel{\&Center\; Dot;}{d}}_1& ...& {\stackrel{\&Center\; Dot;}{d}}_m\end{array}\right],\stackrel{\&Center\; Dot;}{B}={\stackrel{\&Center\; Dot;}{Z}}^T\&Center\; Dot;diag(1,1,-1)\&Center\; Dot;Z+Z^T\&Center\; Dot;diag(1,1,-1)\&CenterDot;\stackrel{\&Center\; Dot;}{Z},$ f(t)=tr((BB ₁ ) ² ), Among them: tr( ) represents the trace of the matrix;

Step 3.2.1), calculate g(t ₂ ): u(t ₂ )=u ⁽¹⁾ +t ₂ p, when g(t ₂ )>0, go to step 3.2.2, otherwise go to step 3.2.6 ;

Step 3.2.2), calculate g(t ₂ -4σ): u(t ₂ -4σ)=u ⁽¹⁾ +(t ₂ -4σ)p, when g(t ₂ -4σ)≤0, then take t _left ＝t ₂ -4σ, t _right ＝t ₂ , go to step 3.2.10, otherwise go to step 3.2.3;

Step 3.2.3), calculate g(t ₂ -8σ): u(t ₂ -8σ)=u ⁽¹⁾ +(t ₂ -8σ)p, calculate g(t ₂ -8σ), when g(t ₂ -8σ)≤0, then take t _left =t ₂ -8σ, t _right =t ₂ -4σ, go to step 3.2.10, otherwise go to step 3.2.4;

Step 3.2.4), calculate g(t ₂ -12σ): u(t ₂ -12σ)=u ⁽¹⁾ +(t ₂ -12σ)p, when g(t ₂ -12σ)≤0, take t _left ＝t ₂ -12σ, t _right ＝t ₂ -8σ, go to step 3.2.10, otherwise go to step 3.2.5;

Step 3.2.5), calculate g(t ₂ -16σ): u(t ₂ -16σ)=u ⁽¹⁾ +(t ₂ -16σ)p, when g(t ₂ -16σ)≤0, take t _left =t ₂ -16σ, t _right =t ₂ -12σ, go to step 3.2.10, otherwise take t ₃ =t ₂ -16σ, go to step 3.3;

Step 3.2.6), calculate g(t ₂ +4σ): u(t ₂ +4σ)=u ⁽¹⁾ +(t ₂ +4σ)p, when g(t ₂ +4σ)>0, then take t _left ＝t ₂ , t _right ＝t ₂ +4σ, go to step 3.2.10, otherwise go to step 3.2.7;

Step 3.2.7), calculate g(t ₂ +8σ): u(t ₂ +8σ)=u ⁽¹⁾ +(t ₂ +8σ)p, when g(t ₂ +8σ)>0, then take t _left ＝t ₂ +4σ, t _right ＝t ₂ +8σ, go to step 3.2.10, otherwise go to step 3.2.8;

Step 3.2.8), calculate g(t ₂ +12σ): u(t ₂ +12σ)=u ⁽¹⁾ +(t ₂ +12σ)p, when g(t ₂ +12σ)>0, then take t _left ＝t ₂ +8σ, t _right ＝t ₂ +12σ, go to step 3.2.10, otherwise go to step 3.2.9;

Step 3.2.9), calculate g(t ₂ +16σ): u(t ₂ +16σ)=u ⁽¹⁾ +(t ₂ +16σ)p, when g(t ₂ +16σ)>0, then take t _left ＝t ₂ +12σ, t _right ＝t ₂ +16σ, go to step 3.2.10, otherwise take t ₃ ＝t ₂ +16σ, go to step 3.3;

Step 3.2.10), take Calculate g(t _middle ): u(t _middle )＝u ⁽¹⁾ +t _middle ·p, when |g(t _middle )|≤10 ^-3 , then take t ₃ =t _middle , go to step 3.3, otherwise go to Step 3.2.11;

Step 3.2.11), when g(t _middle )>0, take t _right ＝t _middle , otherwise take t _left ＝t _middle , when t _rightt -t left _t≤10 ^-2 , take Go to step 3.3, otherwise go to step 3.2.10;

Step 3.3), calculate f(t ₃ ) and the third estimated point u ⁽³⁾ : u ⁽³⁾ =u ⁽¹⁾ +t ₃ p, when f(t ₃ )≤f _min , then take f _min = f(t ₃ ), go to step 3.4;

Step 3.4), set a new initial value for the parameter _t2 corresponding to the second estimated point u ⁽²⁾ : take n=n+1, when is defined, take Go to step 3.2.1; otherwise keep f _min and its corresponding f in step 3.3 Go to step 3.5;

Step 3.5), determine the final estimated value of the position coordinates of the signal source: keep the step 3.4 as the final estimate of the signal source location coordinates.

technical effect

Compared with the prior art, the present invention calculates the arrival distance difference through the time difference of arrival under the framework of the multidimensional scaling method, and is applicable to various sensor formations such as normal sensor formations, approximate linear sensor formations, and linear sensor formations. When the collected data are not completely compatible, the estimated value of the signal source position with the smallest mean square error can still be calculated, which improves the positioning accuracy and accuracy.

Description of drawings

Fig. 1 is the flow chart of the present invention.

detailed description

The embodiments of the present invention are described in detail below. This embodiment is implemented on the premise of the technical solution of the present invention, and detailed implementation methods and specific operating procedures are provided, but the protection scope of the present invention is not limited to the following implementation example.

Example 1

It is known that there are 8 sensor devices on a plane, and the position coordinates of these 8 sensor devices are:

${\mathrm{<span\; class="notranslate">\; u</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; =</span>\left[\begin{array}{c}\mathrm{<span\; class="notranslate">\; 11.8998</span>}\\ \mathrm{<span\; class="notranslate">\; 50.5957</span>}\end{array}\right]<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; u</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; =</span>\left[\begin{array}{c}\mathrm{<span\; class="notranslate">\; 49.8364</span>}\\ \mathrm{<span\; class="notranslate">\; 69.9077</span>}\end{array}\right]<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; u</span>}}_{\mathrm{<span\; class="notranslate">\; 3</span>}}<span\; class="notranslate">\; =</span>\left[\begin{array}{c}\mathrm{<span\; class="notranslate">\; 95.9744</span>}\\ \mathrm{<span\; class="notranslate">\; 89.0903</span>}\end{array}\right]<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; u</span>}}_{\mathrm{<span\; class="notranslate">\; 4</span>}}<span\; class="notranslate">\; =</span>\left[\begin{array}{c}\mathrm{<span\; class="notranslate">\; 34.0386</span>}\\ \mathrm{<span\; class="notranslate">\; 95.9291</span>}\end{array}\right]<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; u</span>}}_{\mathrm{<span\; class="notranslate">\; 5</span>}}<span\; class="notranslate">\; =</span>\left[\begin{array}{c}\mathrm{<span\; class="notranslate">\; 58.5268</span>}\\ \mathrm{<span\; class="notranslate">\; 54.7216</span>}\end{array}\right]<span\; class="notranslate">\; ,</span>$

$u_6=\left[\begin{array}{c}22.3812\\ 13.8624\end{array}\right],u_7=\left[\begin{array}{c}75.1267\\ 14.9294\end{array}\right],u_8=\left[\begin{array}{c}25.5095\\ 25.7508\end{array}\right];$ The real position of the signal source is $u_0=\left[\begin{array}{c}84.0717\\ 25.4282\end{array}\right].$

The three-stage TDOA localization method based on multidimensional scaled subspace analysis includes the following three stages:

The first phase described includes the following steps:

Step 1.1), generate measurement scalar product matrix:

It is known that the signal propagation speed is c, and the normalization is c=1. The difference between the measured arrival time from signal source u ₀ to each sensor u _m and the arrival time from signal source u ₀ to the first sensor u ₁ is: ${\hat{t}}_{51}=-37.7748,{\hat{t}}_{61}=-13.4115,{\hat{t}}_{71}=-62.2653,{\hat{t}}_{81}=-\mathrm{18.3493.}$ Then the difference between the arrival distance from the signal source u ₀ to each sensor u _m and the arrival distance from the signal source u ₀ to the first sensor u ₁ is: ${\hat{d}}_{31}=c\&CenterDot;{\hat{t}}_{31}=-11.4547,{\hat{d}}_{41}=c\&CenterDot;{\hat{t}}_{41}=-10.5034,{\hat{d}}_{51}=c\&CenterDot;{\hat{t}}_{51}=-37.7748,$ ${\hat{d}}_{61}=c\&CenterDot;{\hat{t}}_{61}=-13.4115,{\hat{d}}_{71}=c\&CenterDot;{\hat{t}}_{71}=-62.2653,{\hat{d}}_{81}=c\&CenterDot;{\hat{t}}_{81}=-18.3493;$

It is known that when m=2,...,8 The variance of the measurement error is the same as but The error variance 2σ ² =2×0.3 ² ;

generate matrix of scalar products of measurements get: $\hat{B}={10}^3\&times;\left[\begin{array}{cccccccc}0& -0.6917& -4.0296& -1.2175& -0.3821& -0.6397& -0.6964& -0.2329\\ -0.6917& 0& -1.2056& 0.0237& -0.0074& -1.9208& -0.9675& -1.2680\\ -4.2096& -1.2056& 0& -1.7003& -0.9454& -5.5357& -1.6764& -4.4648\\ -1.2175& 0.0237& -1.7003& 0& 0.0165& -3.1495& -1.4770& -2.0826\\ -0.3821& 0.0074& -0.9454& 0.0165& 0& -1.1912& -0.6296& -0.7760\\ -0.6397& -1.9208& -5.5357& -3.1495& -1.1912& 0& -0.1983& -0.0634\\ -0.6964& -0.9675& -1.6764& -1.4770& -0.6296& -0.1983& 0& -0.3252\\ -0.2329& -1.2680& -4.4648& -2.0826& -0.7760& -0.0634& 0.3252& 0\end{array}\right];$

Step 1.2), correcting the measurement scalar product matrix: $B_1={10}^3\&times;\left[\begin{array}{cccccccc}0& -0.6918& -4.0297& -1.2176& -0.3822& -0.6398& -0.6965& -0.2330\\ -0.6917& 0& -1.2056& 0.0236& -0.0075& -1.9209& -0.9676& -1.2681\\ -4.2097& -1.2056& 0& -1.7004& -0.9455& -5.5358& -1.6764& -4.4649\\ -1.2176& 0.0236& -1.7004& 0& 0.0164& -3.1495& -1.4770& -2.0827\\ -0.3822& 0.0075& -0.9455& 0.0164& 0& -1.1913& -0.6297& -0.7761\\ -0.6398& -1.9209& -5.5358& -3.1495& -1.1913& 0& -0.1984& -0.0635\\ -0.6965& -0.9676& -1.6766& -1.4770& -0.6297& -0.1984& 0& -0.3253\\ -0.2330& -1.2681& -4.4649& -2.0827& -0.7761& -0.0635& 0.3253& 0\end{array}\right];$

Step 1.3), calculating the first estimated point and a distance estimate comprises the following steps:

Step 1.3.1), perform eigenvalue decomposition on B ₁ -σ ² I _M to get: $v_1=\left[\begin{array}{c}0.3192\\ 0.2147\\ 0.5923\\ 0.3385\\ 0.1424\\ 0.4581\\ 0.1955\\ 0.3440\end{array}\right],v_2=\left[\begin{array}{c}0.3014\\ -0.1219\\ -0.6419\\ -0.2215\\ -0.0590\\ 0.5053\\ 0.1312\\ 0.3965\end{array}\right],$ $v_3=\left[\begin{array}{c}0.3658\\ 0.3804\\ -0.4021\\ 0.5600\\ 0.2731\\ -0.2697\\ -0.3141\\ -0.0110\end{array}\right],v_4=\left[\begin{array}{c}-0.0861\\ 0.4891\\ -0.0201\\ -0.4562\\ 0.4610\\ -0.3272\\ 0.3761\\ 0.2895\end{array}\right],v_5=\left[\begin{array}{c}-0.1346\\ -0.1949\\ -0.2087\\ 0.3640\\ 0.2960\\ 0.1066\\ 0.7038\\ -0.4167\end{array}\right]v_6=\left[\begin{array}{c}-0.1195\\ -0.6444\\ 0.0439\\ 0.2027\\ 0.3445\\ -0.3160\\ -0.0632\\ 0.5521\end{array}\right],v_7=\left[\begin{array}{c}-0.2098\\ -0.0480\\ 0.0203\\ -0.1848\\ 0.6725\\ 0.4493\\ -0.4503\\ -0.2492\end{array}\right],$ $v_8=\left[\begin{array}{c}0.7679\\ -0.3154\\ 0.1709\\ -0.3269\\ 0.1749\\ -0.2047\\ 0.0495\\ -0.3161\end{array}\right];$ s ₁ =-1.1334×10 ⁴ , s ₂ =8.5322×10 ³ , s ₃ =2.8357×10 ³ , s ₄ =-66.9765, s ₅ =31.8773, s ₆ =5.7825×10 ^-13 , s ₇ =-5.1324 ×10 ^-13 , s ₈ =-1.1186×10 ^-13 ; k ₁ =5.3656×10 ^-5 , k ₂ =1.0543×10 ^-5 , k ₃ =1.9163×10 ^-5 , k ₄ =0.4293, k ₅ = 1.3348;

Step 1.3.2), obtain a coefficient vector by the linear combination of feature vectors: $v=\left[\begin{array}{c}-0.2165\\ -0.0501\\ -0.2872\\ 0.2901\\ 0.5931\\ 0.0017\\ 1.1008\\ -0.4319\end{array}\right];$

Step 1.3.3), with the coefficient vector v as the weighting coefficient, the position coordinate matrix $\left[\begin{array}{ccc}{x}_1& ...& x_8\\ {\mathrm{the\; y}}_1& ...& {\mathrm{the\; y}}_8\end{array}\right]$ The column vectors of are linearly combined to calculate the first estimated point $u^{\left(1\right)}=\left[\begin{array}{c}83.6662\\ 25.5762\end{array}\right]$ and a distance estimate

Step 1.4), calculate the direction vector p: because $\left|1_m^T{v}_5\right|=0.5155>{10}^{-10},$ so $u_f=\left[\begin{array}{c}83.4118\\ 25.8888\end{array}\right],$ $p=\left[\begin{array}{c}-0.6311\\ 0.7757\end{array}\right]:$

The second stage uses the direction vector p and the distance estimate Calculate the parameter t ₂ corresponding to the second estimated point u ⁽²⁾ : a ₁ =-64.7006, a ₀ =50.7965, Δ=4.1354×10 ³ ; since Δ>0, so

The third phase described includes:

Step 3.1), taking the parameter t ₂ corresponding to the second estimated point u ⁽²⁾ as the initial value: take f _min =∞, n=1, then $t_2=t_2^{\left(1\right)}=129.0074;$

Step 3.2), according to the initial value of the parameter _t2 , use the dichotomy to find the root g(t)=0, and find out the corresponding parameter t3 of the third estimated point u ⁽³⁾ _:

Step 3.2.1), calculation $g\left(t_2\right):u\left(t_2\right)=\left[\begin{array}{c}2.2467\\ 125.6451\end{array}\right],$ g(t ₂ )=-2.0785×10 ⁵ , since g(t ₂ )≤0, go to step 3.2.6;

Step 3.2.6), calculation $g(t_2+4\mathrm{\&sigma;}):u(t_2+4\mathrm{\&sigma;})=\left[\begin{array}{c}1.4893\\ 126.5759\end{array}\right],$ g(t ₂ +4σ)＝-2.0240×10 ⁵ , since g(t ₂ +4σ)≤0, go to step 3.2.7;

Step 3.2.7), calculation $g(t_2+8\mathrm{\&sigma;}):u(t_2+8\mathrm{\&sigma;})=\left[\begin{array}{c}0.7320\\ 127.5067\end{array}\right],$ g(t ₂ +8σ)＝-1.9712×10 ⁵ , since g(t ₂ +8σ)≤0, go to step 3.2.8;

Step 3.2.8), calculation $g(t_2+12\mathrm{\&sigma;}):u(t_2+12\mathrm{\&sigma;})=\left[\begin{array}{c}-0.0254\\ 128.4376\end{array}\right],$ g(t ₂ +12σ)＝-1.9201×10 ⁵ , since g(t ₂ +12σ)≤0, go to step 3.2.9;

Step 3.2.9), calculation $g(t_2+16\mathrm{\&sigma;}):u(t_2+16\mathrm{\&sigma;})=\left[\begin{array}{c}-0.7827\\ 129.3684\end{array}\right],$ g(t ₂ +16σ)=-1.8707×10 ⁵ , because g(t ₂ +16σ)≤0, so take t ₃ =t ₂ +16σ=133.8074, go to step 3.3;

Step 3.3), calculate the third estimation point $u^{\left(3\right)}:u^{\left(3\right)}=u^{\left(1\right)}+t_{3}p=\left[\begin{array}{c}-0.7827\\ 129.3684\end{array}\right],$ f(t ₃ )＝2.9857×10 ⁷ , since f(t ₃ )≤f _min , then take $\hat{u}=u^{\left(3\right)}=\left[\begin{array}{c}-0.7827\\ 129.3684\end{array}\right],$ And take f _min =2.9857×10 ⁷ , go to step 3.4;

Step 3.4), setting a new initial value for the parameter _t2 corresponding to the second estimation point u ⁽²⁾ : getting n=n+1=2, $t_2=t_2^{\left(2\right)}=0.3937,$ Go to step 3.2.1;

Step 3.2.1), calculation $g\left(t_2\right):u\left(t_2\right)=\left[\begin{array}{c}83.4177\\ 25.8816\end{array}\right],$ g(t ₂ )=2.9567×10 ⁴ , since g(t ₂ )>0, go to step 3.2.2;

Step 3.2.2), calculation $g(t_2-4\mathrm{\&sigma;}):u(t_2-4\mathrm{\&sigma;})=\left[\begin{array}{c}84.1750\\ 24.9508\end{array}\right],$ g(t ₂ -4σ)=-6.6714×10 ⁴ , because g(t ₂ -4σ)≤0, so take t _left ＝-0.8063, t _right ＝0.3937, go to step 3.2.10;

Step 3.2.10), calculating |g(t _middle )|: t _middle =-0.2063, $u\left(t_{middle}\right)=\left[\begin{array}{c}83.7963\\ 25.4162\end{array}\right],$ g(t _middle )＝-2.0526×10 ⁴ , because |g(t _middle )|>10 ^-3 , so turn to step 3.2.11;

Step 3.2.11), calculate t _rightt -t _left : since g(t _middle )≤0, take t left t= _-0.2063 , and since t _rightt -t _left >10 ^-2 , go to step 3.2.10;

Step 3.2.10), calculating |g(t _middle )|: t _middle =0.0937, $u\left(t_{middle}\right)=\left[\begin{array}{c}83.6070\\ 25.6489\end{array}\right],$ g(t _middle )＝4.0541×10 ³ , because |g(t _middle )|>10 ^-3 , so turn to step 3.2.11;

Step 3.2.11), calculate t _right -t _left : since g(t _middle )>0, take t _right = _0.0937 , and because t right -t _left >10 ^-2 , go to step 3.2.10;

Step 3.2.10), calculate t _middle =-0.0563, $u\left(t_{middle}\right)=\left[\begin{array}{c}83.7017\\ 25.5326\end{array}\right],$ g(t _middle )＝-8.3555×10 ³ , because |g(t _middle )|>10 ^-3 , so turn to step 3.2.11;

Step 3.2.11), calculate t _rightt -t _left : because g(t _middle )≤0, take t left t= _-0.0563 , and because t _rightt -t _left >10 ^-2 , go to step 3.2.10;

Step 3.2.10), calculate |g(t _middle )|: t _middle =0.0187, $u\left(t_{middle}\right)=\left[\begin{array}{c}83.6543\\ 25.5907\end{array}\right],$ g(t _middle )＝-2.1802×10 ³ , because |g(t _middle )|>10 ^-3 , so turn to step 3.2.11;

Step 3.2.11), calculate t _rightt -t _left : since g(t _middle )≤0, take t _left =0.0187, and since t _rightt -t _left >10 ^-2 , go to step 3.2.10;

Step 3.2.10), calculate |g(t _middle )|: t _middle =0.0562, $u\left(t_{middle}\right)=\left[\begin{array}{c}83.6307\\ 25.6198\end{array}\right],$ g(t _middle )＝929.6178, because |g(t _middle )|>10 ^-3 , so turn to step 3.2.11;

Step 3.2.11), calculate t _right -t _left : since g(t _middle )>0, take t _right = _0.0562 , and because t right -t _left >10 ^-2 , go to step 3.2.10;

Step 3.2.10), calculating |g(t _middle )|: t _middle =0.0375, $u\left(t_{middle}\right)=\left[\begin{array}{c}83.6425\\ 25.6053\end{array}\right],$ g(t _middle )＝-627.1421, because |g(t _middle )|>10 ^-3 , so turn to step 3.2.11;

Step 3.2.11), calculate t _rightt -t _left : since g(t _middle )≤0, take t _left =0.0375, then t _rightt -t _left >10 ^-2 , so go to step 3.2.10;

Step 3.2.10), calculating |g(t _middle )|: t _middle =0.0469, $u\left(t_{middle}\right)=\left[\begin{array}{c}83.6366\\ 25.6126\end{array}\right],$ g(t _middle )＝150.7786, because |g(t _middle )|>10 ^-3 , so turn to step 3.2.11;

Step 3.2.11), calculate t _rightt -t _left : since g(t _middle )>0, take t _right =0.0469, then t _rightt -t _left ≤10 ^-2 , so take t ₃ =0.0422, go to step 3.3 ;

Step 3.3), calculate the third estimation point $u^{\left(3\right)}:u^{\left(3\right)}=\left[\begin{array}{c}83.6395\\ 25.6089\end{array}\right],$ Calculate f(t ₃ )=1.0242×10 ⁴ , because f(t ₃ )≤f _min , so take $\hat{u}=u^{\left(3\right)}=\left[\begin{array}{c}83.6395\\ 25.6089\end{array}\right],$ And take f _min =1.0242×10 ⁴ , go to step 3.4;

Step 3.4), take the parameter t2 corresponding to the second estimated point u ⁽²⁾ _as the initial value: get n=n+1=3, because there is no The definition of , so keep f _min =1.0242×10 ⁴ and its corresponding value in step 3.3 $\hat{u}=\left[\begin{array}{c}83.6395\\ 25.6089\end{array}\right],$ Go to step 3.5;

Step 3.5), determine the final estimated value of the position coordinates of the signal source: when f _min =1.0242×10 ⁴ , $\hat{u}=\left[\begin{array}{c}83.6395\\ 25.6089\end{array}\right]$ is the final estimated value of the position coordinates of the signal source.

The final estimated value of the position coordinates of the signal source calculated through the method of this embodiment $\hat{u}=\left[\begin{array}{c}83.6395\\ 25.6089\end{array}\right]$ and the real coordinates of the signal source position $u_0=\left[\begin{array}{c}84.0717\\ 25.4282\end{array}\right]$ Basically the same, the error is very small.

Claims (4)

1., based on a poor localization method three time of arrival in stage for multidimensional scaling subspace analysis, it is characterized in that, comprise following three phases:

First stage: know plane upper sensor position coordinates, time of arrival is poor, the variance of the measuring error that time of arrival is poor and signal velocity, calculate the measurement scalar product matrix revised, pass through Eigenvalues Decomposition, calculate the first estimation point and a range estimation, and calculate a direction vector;

Subordinate phase: on the basis of the first estimation point, utilizes range estimation and direction vector, calculates the parameter that the second estimation point is corresponding;

Phase III: with the parameter corresponding to the second estimation point for initial value, by dichotomy rooting process, obtain the parameter corresponding to the 3rd estimation point, calculate the 3rd estimation point further, obtain the final estimated value of coordinate, determine the position of signal source thus.

2. differ from localization method three time of arrival in stage based on multidimensional scaling subspace analysis according to claim 1, it is characterized in that, the described first stage comprises following steps:

Step 1.1) generate and measure scalar product matrix: described plane upper sensor position coordinates is u _m=[x _m, y _m] ^t, m=1 ..., M, wherein: M represents number of sensors and quantity is more than or equal to 5, u _mrepresent the position coordinates of m sensor, x _mrepresent the x-axis coordinate of m sensor, y _mrepresent the y-axis coordinate of m sensor;

Obtain arriving range difference with signal velocity according to difference time of arrival: as m=1 work as m=2 ..., during M wherein: represent the signal source u measured ₀to each sensor u _mtime of arrival and signal source u ₀to the 1st sensor u ₁the difference of time of arrival, c represents signal velocity;

Work as m=2 ..., during M, according to measuring error variance and signal velocity obtain error variance be wherein: represent measuring error variance, c represents signal velocity;

Described measurement scalar product matrix is wherein: matrix capable, the jth n column element of the i-th m be ${\&lsqb;\hat{B}\&rsqb;}_{i,j}=\frac{1}{2}{({\hat{d}}_{i1}-{\hat{d}}_{j1})}^2-\frac{1}{2}\&lsqb;{(x_i-x_j)}^2+{(y_i-y_j)}^2\&rsqb;;$

Step 1.2) measurement scalar product matrix is revised: wherein: I _mrepresent M × M unit matrix, 1 _mrepresent that element is all the M dimensional vector of 1;

Step 1.3) calculate the first estimation point and a range estimation comprises the following steps:

Step 1.3.1) to B ₁-σ ²i _mmake Eigenvalues Decomposition: B ₁-σ ²i _m=[v ₁..., v _m] diag (s ₁..., s _m) [v ₁..., v _m] ^t, wherein: v ₁..., v _m∈ R ^mthe vector that pairwise orthogonal and mould are 1, diag (s ₁..., s _m) expression diagonal element is s ₁..., s _mdiagonal matrix, [v ₁..., v _m] ^tthe transposition of representing matrix, to s ₁..., s _mdescending sort is carried out according to absolute value, namely | s ₁|>=...>=| s _m|;

Step 1.3.2) obtain coefficient vector v:v by the linear combination of proper vector ₁..., v _mfor R ^mone group of orthonormal basis, therefore R ^min coefficient vector v can be expressed as v ₁..., v _mlinear combination v=k ₁v ₁+ ...+k ₅v ₅, work as m=1 ..., when 5, base v _mcombination coefficient and when m>=6, base v _mcombination coefficient k ₆..., k _mbe 0;

Step 1.3.3) using coefficient vector v as weighting coefficient, to position coordinates matrix $\left[\begin{array}{ccc}{x}_1& ...& x_M\\ y_1& ...& y_M\end{array}\right]$ Column vector do linear combination, calculate the first estimation point thus $u^{\left(1\right)}=\frac{\left[\begin{array}{ccc}{x}_1& ...& x_M\\ y_1& ...& y_M\end{array}\right]v}{1_M^{T}v}$ With a range estimation $d_1^{\left(1\right)}=-\frac{\left[\begin{array}{ccc}{\hat{d}}_{11}& ...& {\hat{d}}_{M1}\end{array}\right]v}{1_M^{T}v};$

Step 1.4) calculated direction vector p: when then otherwise $p=\frac{{\&lsqb;-(y_2-y_1),x_2-x_1\&rsqb;}^T}{\left|\right|u_2-u_1\left|\right|},$ Wherein: $u_f=\frac{\left[\begin{array}{ccc}{x}_1& ...& x_M\\ y_1& ...& y_M\end{array}\right]v_5}{1_M^T{v}_5},$ || || represent the euclideam norm of vector.

3. according to claim 1 based on multidimensional scaling subspace analysis three time of arrival in stage difference localization method, it is characterized in that, described subordinate phase on the basis of the first estimation point, the range estimation described in utilization with direction vector p, calculate the intermediate parameters t that the second estimation point is corresponding ₂:

Definition $t_2=t_2^{\left(n\right)},a_0=\left|\right|u^{\left(1\right)}-u_1|{|}^2-{\left(d_1^{\left(1\right)}\right)}^2,$ a ₁＝p ^T(u ⁽¹⁾-u ₁)， $\mathrm{\&Delta;}=a_1^2-a_0;$ Then as Δ >0, define during n=1 $t_2^{\left(1\right)}=-a_1+\sqrt{\mathrm{\&Delta;}},$ Define during n=2 $t_2^{\left(2\right)}=-a_1-\sqrt{\mathrm{\&Delta;}};$ Otherwise, define during n=1 $t_2^{\left(1\right)}=-a_1,$ During n=2 without definition.

4. differ from localization method three time of arrival in stage based on multidimensional scaling subspace analysis according to claim 1, it is characterized in that, the described phase III comprises following steps:

Step 3.1) with the second estimation point u ⁽²⁾corresponding parametric t ₂for initial value: definition f _min=∞, n=1,

Step 3.2) process of dichotomy rooting g (t)=0, find out the 3rd estimation point u ⁽³⁾corresponding parametric t ₃, wherein: note u (t)=[x _u, y _u] ^t, p=[l _x, l _y] ^t, to m=1 ..., M, makes w _m=[x _u-x _m, y _u-y _m] ^t, ${\stackrel{\&CenterDot;}{d}}_m=p^T{w}_m/\left|\right|w_m\left|\right|,$ Further order $Z=\left[\begin{array}{ccc}{x}_1-x_u& ...& x_M-x_u\\ y_1-y_u& ...& y_M-y_u\\ d_1& ...& d_M\end{array}\right],$ B＝Z ^T·diag(1,1,-1)·Z， $\stackrel{\&CenterDot;}{Z}=\left[\begin{array}{ccc}-l_x& ...& -l_x\\ -l_y& ...& -l_y\\ {\stackrel{\&CenterDot;}{d}}_1& ...& {\stackrel{\&CenterDot;}{d}}_M\end{array}\right],$ B＝Z ^T·diag(1,1,-1)·Z，f(t)＝tr((B-B ₁) ²)， $g\left(t\right)=2tr\left(\right(B-B_1\left)\stackrel{\&CenterDot;}{B}\right),$ Wherein: the mark of tr () representing matrix;

Step 3.3) calculate f (t ₃) and the 3rd estimation point u ⁽³⁾: u ⁽³⁾=u ⁽¹⁾+ t ₃p, as f (t ₃)≤f _min, then get f _min=f (t ₃), go to step 3.4;

Step 3.4) with the second estimation point u ⁽²⁾corresponding parametric t ₂for initial value: get n=n+1, when there is definition, then get go to step 3.2; Otherwise the f retained in step 3.3 _minand correspondence go to step 3.5;

Step 3.5) determine the final estimated value of source location coordinate: step 3.4 is retained as signal source u ₀actual coordinate.