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

Abstract

The invention provides a three-stage arrival time difference positioning method based on multidimensional scaling sub space analysis in the field of a wireless network and a mobile calculation field. The method comprises three stages which comprise a first stage of obtaining a sensor position coordinate on a plane, an arrival time difference, the variance of measurement error of the arrival time difference and a signal spread speed, calculating a corrected measurement scalar product matrix, through characteristic value decomposition, calculating a first estimation point and a distance estimation value, and calculating a direction vector, a second stage of calculating the parameter corresponding to a second estimation point on the basis of the first estimation point by using the distance estimation value and the direction vector, and a third stage of taking the parameter corresponding to the second estimation point as an initial value, through a dichotomy root process, obtaining the parameter corresponding to the third estimation point, further calculating the third estimation point, obtaining the final esitmation value of the coordinate, and thus determining the position of a signal source. The three-stage arrival time difference positioning method has advantages of wide range of application, high positioning precision and high accuracy.

Description

Three time of arrival in stage based on multidimensional scaling subspace analysis differ from localization method

Technical field

What the present invention relates to is wireless network and mobile computing field, specifically differs from localization method a kind of three time of arrival in stage based on multidimensional scaling subspace analysis.

Background technology

In the applications such as radar, sonar, mobile communication, multimedia, wireless sensor network, usually face a major issue, namely according to difference information time of arrival, a signal source is positioned.What is called difference time of arrival refers to: send signal by signal source, by distributing in space, position sensor that is known and time phase mutually synchronization receives this signal, and measuring-signal arrives the time of other sensors, calculate time that signal that signal source sends arrives each sensor and the difference of time arriving first sensor thus.

Through finding the retrieval of prior art, He ?WenWei, RongPeng, QunWan, Zhang ?XinChen, the paper " Multidimensionalscalinganalysisforpassivemovingtargetloc alizationwithTDOAandFDOAmeasurements " that & Shang ?FuYe delivered in volume 3 phase in March in 2010 the 58th of periodical IEEETransactionsonSignalProcessing, proposes and differs from the time of arrival utilizing multidimensional scaling to analyze passive sensor array and the poor method obtaining source location of arrival rate.But in the method when sensor array becomes wire or is similar to wire, required for do inversion operation matrix will become rank deficient matrix or height ill-condition matrix, by force it is inverted and will produce huge positioning error.

Chinese patent literature CN103648164A, open (bulletin) day 2014.03.19, disclose a kind of based on differing from time of arrival and the wireless-sensor network distribution type localization method of Gossip algorithm.Be characterized in, one, anchor node obtain self position coordinates; Two, Distributed Time is realized synchronous; Three, anchor node wakes monitoring unknown node at random up; Four, wake anchor node up and preserve Received signal strength moment and local coordinate system; Five, all anchor nodes whether preserve by settling signal monitoring and data; Six, anchor node j receives the data of its all M adjacent anchor node; Seven, the initial estimate of anchor node j for unknown node the unknown is obtained; Eight, all anchor nodes obtain unknown node position initial estimate; Nine, Gossip algorithm Stochastic choice adjacent anchor node switching locator data is run; Ten, algorithm stops.But the unknown node location estimation value that in the method, different anchor node collects is likely not quite identical, the location estimation value that how to calculate maximum likelihood from obtained positional information is not mentioned in application.

Summary of the invention

The present invention is directed to prior art above shortcomings, propose poor localization method a kind of three time of arrival in stage based on multidimensional scaling subspace analysis, by making subspace analysis to measurement scalar product matrix, in three stages source location is estimated, obtain the more accurate position of signal source.

The present invention is achieved by the following technical solutions,

The present invention includes following steps:

Step 1, according to plane upper sensor position coordinates, time of arrival poor (time that the signal that namely signal source sends arrives each sensor and the difference of time arriving first sensor) obtaining measuring scalar product matrix, and by described time of arrival difference the variance of measuring error revise measuring scalar product matrix; By carrying out Eigenvalues Decomposition process to revised measurement scalar product matrix, calculate the first estimation point, and calculate a direction vector, and a range estimation, specifically comprise the 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 wherein: represent measuring error variance, c represents signal velocity;

Described measurement scalar product matrix is wherein: matrix the i-th row, jth column element be

${\[\hat{B}\]}_{i,j}=\frac{1}{2}{({\hat{d}}_{i1}-{\hat{d}}_{j1})}^2-\frac{1}{2}\[{(x_i-x_j)}^2+{(y_i-y_j)}^2\];$

Step 1.2), measurement scalar product matrix is revised: obtain revised measurement scalar product matrix 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 range estimation, specifically comprise 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, eigenvalue matrix diag (s ₁..., s _m) expression diagonal element is s ₁..., s _mdiagonal matrix, to eigenwert s ₁..., s _mdescending sort is carried out according to absolute value, namely | s ₁|>=...>=| s _m|;

Step 1.3.2), obtain coefficient vector by the linear combination of proper vector: v ₁..., 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 of this estimation 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), calculated direction vector p: when then otherwise $p=\frac{{\[-(y_2-y_1),x_2-x_1\]}^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;

Step 2, 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{\Δ}=a_1^2-a_0;$ Then as Δ >0, define during n=1 define during n=2 otherwise, define during n=1 during n=2 without definition;

Step 3, with parameter corresponding to the second estimation point for initial value, by dichotomy rooting process, obtain the parameter that the 3rd estimation point is corresponding, calculate the 3rd estimation point further, obtain the final estimated value of coordinate, determine the position of signal source thus, concrete steps comprise:

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

Step 3.2), dichotomy rooting g (t)=0, find out the 3rd estimation point u ⁽³⁾corresponding parametric t ₃, the component of u (t) is designated as x _u, y _u, the component of p is designated as l _x, l _y, that is: 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, d _m=|| w _m||, ${\stackrel{\·}{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{\·}{Z}=\left[\begin{array}{ccc}-l_x& ...& -l_x\\ -l_y& ...& -l_y\\ {\stackrel{\·}{d}}_1& ...& {\stackrel{\·}{d}}_M\end{array}\right],\stackrel{\·}{B}={\stackrel{\·}{Z}}^T\·diag(1,1,-1)\·Z+Z^T\·diag(1,1,-1)\·\stackrel{\·}{Z},$ F (t)=tr ((B-B ₁) ²), wherein: the mark of tr () representing matrix;

Step 3.2.1), calculate g (t ₂): u (t ₂)=u ⁽¹⁾+ t ₂p, as g (t ₂) >0, then 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, as g (t ₂-4 σ)≤0, then get 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 σ), as g (t ₂-8 σ)≤0, then get 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, as g (t ₂-12 σ)≤0, then get 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, as g (t ₂-16 σ)≤0, then get t _left=t ₂-16 σ, t _right=t ₂-12 σ, go to step 3.2.10, otherwise get t ₃=t ₂-16 σ, go to step 3.3;

Step 3.2.6), calculate g (t ₂+ 4 σ): u (t ₂+ 4 σ)=u ⁽¹⁾+ (t ₂+ 4 σ) p, as g (t ₂+ 4 σ) >0, then get 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, as g (t ₂+ 8 σ) >0, then get 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, as g (t ₂+ 12 σ) >0, then get 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, as g (t ₂+ 16 σ) >0, then get t _left=t ₂+ 12 σ, t _right=t ₂+ 16 σ, go to step 3.2.10, otherwise get t ₃=t ₂+ 16 σ, go to step 3.3;

Step 3.2.10), get calculate g (t _middle): u (t _middle)=u ⁽¹⁾+ t _middlep, when | g (t _middle) |≤10 ^-3, then t is got ₃=t _middle, go to step 3.3, otherwise go to step 3.2.11;

Step 3.2.11), as g (t _middle) >0, then get t _right=t _middle, otherwise get t _left=t _middle, work as t _rightt-t _left≤10 ^-2, then get go to step 3.3, otherwise go to step 3.2.10;

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), to the second estimation point u ⁽²⁾corresponding parametric t ₂set new initial value: get n=n+1, when there is definition, then get go to step 3.2.1; 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 the final estimated value of source location coordinate.

Technique effect

Compared with prior art, the present invention is calculated by difference time of arrival under multidimensional scaling framework and arrives range difference, be applicable to the various sensor formation such as nominal sensor formation, approximate wire sensor array type, line sensor formation, and when sensor image data is also not exclusively compatible, still the source location estimated value that square error is minimum be can calculate, positioning precision and accuracy improve.

Accompanying drawing explanation

Fig. 1 is process flow diagram of the present invention.

Embodiment

Elaborate to embodiments of the invention below, the present embodiment is implemented under premised on technical solution of the present invention, give detailed embodiment and concrete operating process, but protection scope of the present invention is not limited to following embodiment.

Embodiment 1

8 sensor devices in a known plane, the position coordinates of these 8 sensor devices is respectively:

$u_1=\left[\begin{array}{c}11.8998\\ 50.5957\end{array}\right],u_2=\left[\begin{array}{c}49.8364\\ 69.9077\end{array}\right],u_3=\left[\begin{array}{c}95.9744\\ 89.0903\end{array}\right],u_4=\left[\begin{array}{c}34.0386\\ 95.9291\end{array}\right],u_5=\left[\begin{array}{c}58.5268\\ 54.7216\end{array}\right],$

$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 actual position of signal source is $u_0=\left[\begin{array}{c}84.0717\\ 25.4282\end{array}\right].$

Three time of arrival in stage based on multidimensional scaling subspace analysis, poor localization method comprised following three phases:

The described first stage comprises the following steps:

Step 1.1), generate and measure scalar product matrix:

Known signal velocity of propagation is c, is normalized to c=1.Measure signal source u ₀to each sensor u _mtime of arrival and signal source u ₀to the 1st sensor u ₁the difference of time of arrival be: ${\hat{t}}_{51}=-37.7748,{\hat{t}}_{61}=-13.4115,{\hat{t}}_{71}=-62.2653,{\hat{t}}_{81}=-\mathrm{18.3493.}$ Then signal source u ₀to each sensor u _marrival distance and signal source u ₀to the 1st sensor u ₁arrival distance difference be: ${\hat{d}}_{31}=c\·{\hat{t}}_{31}=-11.4547,{\hat{d}}_{41}=c\·{\hat{t}}_{41}=-10.5034,{\hat{d}}_{51}=c\·{\hat{t}}_{51}=-37.7748,$ ${\hat{d}}_{61}=c\·{\hat{t}}_{61}=-13.4115,{\hat{d}}_{71}=c\·{\hat{t}}_{71}=-62.2653,{\hat{d}}_{81}=c\·{\hat{t}}_{81}=-18.3493;$

Knownly work as m=2 ..., when 8 the variance of measuring error be all then error variance 2 σ ²=2 × 0.3 ²;

Generate and measure scalar product matrix obtain: $\hat{B}={10}^3\×\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), measurement scalar product matrix is revised: $B_1={10}^3\×\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), calculate the first estimation point and a range estimation comprises the following steps:

Step 1.3.1), to B ₁-σ ²i _mobtain as Eigenvalues Decomposition: $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 proper vector: $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), using coefficient vector v as weighting coefficient, to position coordinates matrix $\left[\begin{array}{ccc}{x}_1& ...& x_8\\ y_1& ...& y_8\end{array}\right]$ Column vector do linear combination, calculate the first estimation point thus $u^{\left(1\right)}=\left[\begin{array}{c}83.6662\\ 25.5762\end{array}\right]$ With a range estimation

Step 1.4), calculated direction vector p: because $\left|1_M^T{v}_5\right|=0.5155>{10}^{-10},$ Therefore $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]:$

Described subordinate phase utilization orientation vector p and range estimation calculate the second estimation point u ⁽²⁾corresponding parametric t ₂: a ₁=-64.7006, a ₀=50.7965, Δ=4.1354 × 10 ³; Because of Δ > 0, therefore

The described phase III comprises:

Step 3.1), with the second estimation point u ⁽²⁾corresponding parametric t ₂for initial value: get f _min=∞, n=1, then $t_2=t_2^{\left(1\right)}=129.0074;$

Step 3.2), according to parametric t ₂initial value, utilize dichotomy rooting g (t)=0, find out the 3rd estimation point u ⁽³⁾corresponding parametric t ₃:

Step 3.2.1), calculate $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 ⁵, because g is (t ₂)≤0, therefore go to step 3.2.6;

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

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

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

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

Step 3.3), calculate the 3rd 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 ⁷, due to f (t ₃)≤f _min, then get $\hat{u}=u^{\left(3\right)}=\left[\begin{array}{c}-0.7827\\ 129.3684\end{array}\right],$ And get f _min=2.9857 × 10 ⁷, go to step 3.4;

Step 3.4), to the second estimation point u ⁽²⁾corresponding parametric t ₂set new initial value: get n=n+1=2, $t_2=t_2^{\left(2\right)}=0.3937,$ Go to step 3.2.1;

Step 3.2.1), calculate $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 ⁴, because g is (t ₂) >0, therefore go to step 3.2.2;

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

Step 3.2.10), calculate | 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, therefore go to step 3.2.11;

Step 3.2.11), calculate t _rightt-t _left: because g (t _middle)≤0, therefore get t _left=-0.2063, because t _rightt-t _left>10 ^-2, therefore go to step 3.2.10;

Step 3.2.10), calculate | 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, therefore go to step 3.2.11;

Step 3.2.11), calculate t _rightt-t _left: because g (t _middle) >0, therefore get t _right=0.0937, because t _rightt-t _left>10 ^-2, therefore 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, therefore go to step 3.2.11;

Step 3.2.11), calculate t _rightt-t _left: because g (t _middle)≤0, therefore get t _left=-0.0563, because t _rightt-t _left>10 ^-2, therefore 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, therefore go to step 3.2.11;

Step 3.2.11), calculate t _rightt-t _left: because g (t _middle)≤0, therefore get t _left=0.0187, because t _rightt-t _left>10 ^-2, therefore 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, therefore go to step 3.2.11;

Step 3.2.11), calculate t _rightt-t _left: because g (t _middle) >0, therefore get t _right=0.0562, because t _rightt-t _left>10 ^-2, therefore go to step 3.2.10;

Step 3.2.10), calculate | 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, therefore go to step 3.2.11;

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

Step 3.2.10), calculate | 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, therefore go to step 3.2.11;

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

Step 3.3), calculate the 3rd 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 is (t ₃)≤f _min, therefore get $\hat{u}=u^{\left(3\right)}=\left[\begin{array}{c}83.6395\\ 25.6089\end{array}\right],$ And get f _min=1.0242 × 10 ⁴, go to step 3.4;

Step 3.4), with the second estimation point u ⁽²⁾corresponding parametric t ₂for initial value: get n=n+1=3, because nothing definition, therefore retain the f in step 3.3 _min=1.0242 × 10 ⁴and correspondence $\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 source location coordinate: f _min=1.0242 × 10 ⁴time, $\hat{u}=\left[\begin{array}{c}83.6395\\ 25.6089\end{array}\right]$ For the final estimated value of source location coordinate.

The final estimated value of the source location coordinate of acquisition is calculated through the present embodiment mode $\hat{u}=\left[\begin{array}{c}83.6395\\ 25.6089\end{array}\right]$ With source location true coordinate $u_0=\left[\begin{array}{c}84.0717\\ 25.4282\end{array}\right]$ Substantially conform to, error is minimum.

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 ${\[\hat{B}\]}_{i,j}=\frac{1}{2}{({\hat{d}}_{i1}-{\hat{d}}_{j1})}^2-\frac{1}{2}\[{(x_i-x_j)}^2+{(y_i-y_j)}^2\];$

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{{\[-(y_2-y_1),x_2-x_1\]}^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{\Δ}=a_1^2-a_0;$ Then as Δ >0, define during n=1 $t_2^{\left(1\right)}=-a_1+\sqrt{\mathrm{\Δ}},$ Define during n=2 $t_2^{\left(2\right)}=-a_1-\sqrt{\mathrm{\Δ}};$ 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{\·}{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{\·}{Z}=\left[\begin{array}{ccc}-l_x& ...& -l_x\\ -l_y& ...& -l_y\\ {\stackrel{\·}{d}}_1& ...& {\stackrel{\·}{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{\·}{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.