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

Three-stage arrival time difference positioning method based on multidimensional scaling sub space analysis Download PDF

Info

Publication number
CN105259533A
CN105259533A CN201510710744.7A CN201510710744A CN105259533A CN 105259533 A CN105259533 A CN 105259533A CN 201510710744 A CN201510710744 A CN 201510710744A CN 105259533 A CN105259533 A CN 105259533A
Authority
CN
China
Prior art keywords
calculate
arrival
estimation point
vector
stage
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN201510710744.7A
Other languages
Chinese (zh)
Other versions
CN105259533B (en
Inventor
蒋武扬
徐昌庆
裴凌
郁文贤
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Shanghai Jiao Tong University
Original Assignee
Shanghai Jiao Tong University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Shanghai Jiao Tong University filed Critical Shanghai Jiao Tong University
Priority to CN201510710744.7A priority Critical patent/CN105259533B/en
Publication of CN105259533A publication Critical patent/CN105259533A/en
Application granted granted Critical
Publication of CN105259533B publication Critical patent/CN105259533B/en
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S5/00Position-fixing by co-ordinating two or more direction or position line determinations; Position-fixing by co-ordinating two or more distance determinations
    • G01S5/02Position-fixing by co-ordinating two or more direction or position line determinations; Position-fixing by co-ordinating two or more distance determinations using radio waves
    • G01S5/06Position of source determined by co-ordinating a plurality of position lines defined by path-difference measurements
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S5/00Position-fixing by co-ordinating two or more direction or position line determinations; Position-fixing by co-ordinating two or more distance determinations
    • G01S5/02Position-fixing by co-ordinating two or more direction or position line determinations; Position-fixing by co-ordinating two or more distance determinations using radio waves
    • G01S5/0252Radio frequency fingerprinting

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Position Fixing By Use Of Radio Waves (AREA)

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.

经过对现有技术的检索发现,He‐WenWei,RongPeng,QunWan,Zhang‐XinChen,&Shang‐FuYe在期刊IEEETransactionsonSignalProcessing的2010年3月第58卷3期发表的论文“MultidimensionalscalinganalysisforpassivemovingtargetlocalizationwithTDOAandFDOAmeasurements”,提出了利用多维标度法分析无源传感器阵列的到达时间差和到达频率差得到信号源位置的方法。但是该方法中当传感器阵列成为线状或近似于线状时,所需要做求逆运算的矩阵将成为秩亏矩阵或高度病态矩阵,强行对其求逆将产生巨大的定位误差。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.

中国专利文献号CN103648164A,公开(公告)日2014.03.19,公开了一种基于到达时间差和Gossip算法的无线传感器网络分布式定位方法。其特点是,一、锚节点获取自身位置坐标;二、实现分布式时间同步;三、锚节点随机唤醒监测未知节点;四、唤醒锚节点保存接收信号时刻和本地坐标;五、所有锚节点是否全部完成信号监测和数据保存;六、锚节点j收到其所有M个相邻锚节点的数据;七、获取锚节点j对于未知节点未知的初始估计值;八、所有锚节点获取未知节点位置初始估计值;九、运行Gossip算法随机选择相邻锚节点交换定位数据;十、算法终止。但是该方法中不同锚节点采集到的未知节点位置估计值有可能不完全一致,如何从所获得的位置信息中计算出最大可能性的位置估计值在申请中并未提及。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:

步骤1、根据平面上传感器位置坐标、到达时间差(即信号源所发出的信号到达各个传感器的时间与到达第一个传感器的时间之差)得到测量标量乘积矩阵,并通过所述的到达时间差的测量误差的方差对测量标量乘积矩阵进行修正;通过对修正后的测量标量乘积矩阵进行特征值分解处理,计算得到第一估计点,并计算出一个方向向量,和一个距离估计值,具体包括以下步骤: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:

步骤1.1)、生成测量标量乘积矩阵所述的平面上传感器位置坐标为um=[xm,ym]T,m=1,...,M,其中:M表示传感器数量且数量大于等于5个,um表示第m个传感器的位置坐标,xm表示第m个传感器的x轴坐标,ym表示第m个传感器的y轴坐标;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;

根据到达时间差与信号传播速度得到到达距离差:当m=1时当m=2,...,M时其中:表示已测量到的信号源u0到各个传感器um的到达时间与信号源u0到第1个传感器u1的到达时间之差,c表示信号传播速度;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 0 to each sensor u m and the arrival time from signal source u 0 to the first sensor u 1 , c indicates the signal propagation speed;

当m=2,...,M时,根据的测量误差方差与信号传播速度得到的误差方差其中:表示的测量误差方差,c表示信号传播速度;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;

所述的测量标量乘积矩阵为其中:矩阵的第i行、第j列元素为The measured scalar product matrix is where: matrix The elements of row i and column j of are

[[ BB ^^ ]] ii ,, jj == 11 22 (( dd ^^ ii 11 -- dd ^^ jj 11 )) 22 -- 11 22 [[ (( xx ii -- xx jj )) 22 ++ (( ythe y ii -- ythe y jj )) 22 ]] ;;

步骤1.2)、对测量标量乘积矩阵进行修正:得到修正后的测量标量乘积矩阵其中:IM表示M×M单位矩阵,1M表示元素全部为1的M维列向量;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;

步骤1.3)、计算第一估计点和距离估计值,具体包括以下步骤:Step 1.3), calculating the first estimated point and the estimated distance, specifically include the following steps:

步骤1.3.1)、对B12IM作特征值分解:Step 1.3.1), perform eigenvalue decomposition on B 12 I M :

B12IM=[v1,...,vM]diag(s1,...,sM)[v1,...,vM]T,其中:v1,...,vM∈RM是两两正交、而且模均为1的向量,特征值矩阵diag(s1,...,sM)表示对角元为s1,...,sM的对角矩阵,对特征值s1,...,sM依绝对值进行降序排列,即|s1|≥...≥|sM|;B 12 I M =[v 1 ,...,v M ]diag(s 1 ,...,s M )[v 1 ,...,v M ] T , where: v 1 ,. ..,v M ∈ R M is a pairwise orthogonal vector with a modulus of 1. The eigenvalue matrix diag(s 1 ,...,s M ) indicates that the diagonal elements are s 1 ,...,s The diagonal matrix of M , the eigenvalues s 1 ,...,s M are arranged in descending order of absolute value, that is, |s 1 |≥...≥|s M |;

步骤1.3.2)、由特征向量的线性组合得到系数向量:v1,...,vM为RM的一组标准正交基,故RM中的系数向量v可以表示为v1,...,vM的线性组合v=k1v1+...+k5v5,当m=1,...,5时,基vm的组合系数而当m≥6时,基vm的组合系数k6,...,kM均为0;Step 1.3.2), the coefficient vector is obtained by the linear combination of the eigenvectors: v 1 ,...,v M is a set of orthonormal basis of RM , so the coefficient vector v in RM can be expressed as v 1 , ..., the linear combination of v M v=k 1 v 1 +...+k 5 v 5 , when m=1,...,5, the combination coefficient of base v m And when m≥6, the combination coefficients k 6 ,...,k M of base v m are all 0;

步骤1.3.3)、以系数向量v作为加权系数,对位置坐标矩阵 x 1 ... x M y 1 ... y M 的列向量作线性组合,由此计算出第一估计点 u ( 1 ) = x 1 ... x M y 1 ... y M v 1 M T v 和该估计点的一个距离估计值 d 1 ( 1 ) = - d ^ 11 ... d ^ M 1 v 1 M T v ; Step 1.3.3), with the coefficient vector v as the weighting coefficient, the position coordinate matrix x 1 ... x m the y 1 ... the y m The column vectors of are linearly combined to calculate the first estimated point u ( 1 ) = x 1 ... x m the y 1 ... the y m v 1 m T v and an estimate of the distance from the estimated point d 1 ( 1 ) = - d ^ 11 ... d ^ m 1 v 1 m T v ;

步骤1.4)、计算方向向量p:当否则 p = [ - ( y 2 - y 1 ) , x 2 - x 1 ] T | | u 2 - u 1 | | , 其中: u f = x 1 ... x M y 1 ... y M v 5 1 M T v 5 , ||·||表示向量的欧几里德范数;Step 1.4), calculate the direction vector p: when but otherwise p = [ - ( the y 2 - the y 1 ) , x 2 - x 1 ] T | | u 2 - u 1 | | , in: u f = x 1 ... x m the y 1 ... the y m v 5 1 m T v 5 , ||·|| represents the Euclidean norm of the vector;

步骤2、在第一估计点的基础上,利用所述的距离估计值和方向向量p,计算得到第二估计点对应的中间参数t2Step 2, on the basis of the first estimated point, using the estimated distance and the direction vector p, calculate the intermediate parameter t 2 corresponding to the second estimated point:

定义 t 2 = t 2 ( n ) , a 0 = | | u ( 1 ) - u 1 | | 2 - ( d 1 ( 1 ) ) 2 , a1=pT(u(1)-u1), Δ = a 1 2 - a 0 ; 则当Δ>0,n=1时定义n=2时定义否则,n=1时定义n=2时无定义;definition t 2 = t 2 ( no ) , a 0 = | | u ( 1 ) - u 1 | | 2 - ( d 1 ( 1 ) ) 2 , a 1 =p T (u (1) -u 1 ), Δ = a 1 2 - a 0 ; Then when Δ>0, n=1, define Define when n=2 Otherwise, define when n=1 When n=2 undefined;

步骤3、以第二估计点对应的参数为初始值,通过二分法求根过程,得到第三估计点对应的参数,进一步计算第三估计点,得到坐标的最终估计值,由此确定信号源的位置,具体步骤包括: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:

步骤3.1)、设定第二估计点u(2)对应的参数t2为初始值:定义fmin=∞,n=1, Step 3.1), setting the parameter t2 corresponding to the second estimation point u (2) as the initial value: define f min =∞, n=1,

步骤3.2)、二分法求根g(t)=0,找出第三估计点u(3)对应的参数t3,将u(t)的分量记为xu,yu,将p的分量记为lx,ly,即:u(t)=[xu,yu]T,p=[lx,ly]T,对m=1,...,M,令wm=[xu-xm,yu-ym]T,dm=||wm||, d · m = p T w m / | | w m | | , 进一步令 Z = x 1 - x u ... x M - x u y 1 - y u ... y M - y u d 1 ... d M , B=ZT·diag(1,1,-1)·Z, Z · = - l x ... - l x - l y ... - l y d · 1 ... d · M , B · = Z · T · d i a g ( 1 , 1 , - 1 ) · Z + Z T · d i a g ( 1 , 1 , - 1 ) · Z · , f(t)=tr((B-B1)2),其中:tr(·)表示矩阵的迹;Step 3.2), find the root g(t)=0 by dichotomy, find out the parameter t 3 corresponding to the third estimation point u (3) , 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 ||, d · m = p T w m / | | w m | | , further order Z = x 1 - x u ... x m - x u the y 1 - the y u ... the y m - the y u d 1 ... d m , B=Z T diag(1,1,-1) Z, Z · = - l x ... - l x - l the y ... - l the y d &Center Dot; 1 ... d &Center Dot; m , B &Center Dot; = Z &Center Dot; T &Center Dot; d i a g ( 1 , 1 , - 1 ) &Center Dot; Z + Z T &Center Dot; d i a g ( 1 , 1 , - 1 ) · Z &Center Dot; , f(t)=tr((BB 1 ) 2 ), Among them: tr( ) represents the trace of the matrix;

步骤3.2.1)、计算g(t2):u(t2)=u(1)+t2p,当g(t2)>0,则转步骤3.2.2,否则转步骤3.2.6;Step 3.2.1), calculate g(t 2 ): u(t 2 )=u (1) +t 2 p, when g(t 2 )>0, go to step 3.2.2, otherwise go to step 3.2.6 ;

步骤3.2.2)、计算g(t2-4σ):u(t2-4σ)=u(1)+(t2-4σ)p,当g(t2-4σ)≤0,则取tleft=t2-4σ,tright=t2,转步骤3.2.10,否则转步骤3.2.3;Step 3.2.2), calculate g(t 2 -4σ): u(t 2 -4σ)=u (1) +(t 2 -4σ)p, when g(t 2 -4σ)≤0, then take t left =t 2 -4σ, t right =t 2 , go to step 3.2.10, otherwise go to step 3.2.3;

步骤3.2.3)、计算g(t2-8σ):u(t2-8σ)=u(1)+(t2-8σ)p,计算g(t2-8σ),当g(t2-8σ)≤0,则取tleft=t2-8σ,tright=t2-4σ,转步骤3.2.10,否则转步骤3.2.4;Step 3.2.3), calculate g(t 2 -8σ): u(t 2 -8σ)=u (1) +(t 2 -8σ)p, calculate g(t 2 -8σ), when g(t 2 -8σ)≤0, then take t left =t 2 -8σ, t right =t 2 -4σ, go to step 3.2.10, otherwise go to step 3.2.4;

步骤3.2.4)、计算g(t2-12σ):u(t2-12σ)=u(1)+(t2-12σ)p,当g(t2-12σ)≤0,则取tleft=t2-12σ,tright=t2-8σ,转步骤3.2.10,否则转步骤3.2.5;Step 3.2.4), calculate g(t 2 -12σ): u(t 2 -12σ)=u (1) +(t 2 -12σ)p, when g(t 2 -12σ)≤0, take t left =t 2 -12σ, t right =t 2 -8σ, go to step 3.2.10, otherwise go to step 3.2.5;

步骤3.2.5)、计算g(t2-16σ):u(t2-16σ)=u(1)+(t2-16σ)p,当g(t2-16σ)≤0,则取tleft=t2-16σ,tright=t2-12σ,转步骤3.2.10,否则取t3=t2-16σ,转步骤3.3;Step 3.2.5), calculate g(t 2 -16σ): u(t 2 -16σ)=u (1) +(t 2 -16σ)p, when g(t 2 -16σ)≤0, take t left =t 2 -16σ, t right =t 2 -12σ, go to step 3.2.10, otherwise take t 3 =t 2 -16σ, go to step 3.3;

步骤3.2.6)、计算g(t2+4σ):u(t2+4σ)=u(1)+(t2+4σ)p,当g(t2+4σ)>0,则取tleft=t2,tright=t2+4σ,转步骤3.2.10,否则转步骤3.2.7;Step 3.2.6), calculate g(t 2 +4σ): u(t 2 +4σ)=u (1) +(t 2 +4σ)p, when g(t 2 +4σ)>0, then take t left =t 2 , t right =t 2 +4σ, go to step 3.2.10, otherwise go to step 3.2.7;

步骤3.2.7)、计算g(t2+8σ):u(t2+8σ)=u(1)+(t2+8σ)p,当g(t2+8σ)>0,则取tleft=t2+4σ,tright=t2+8σ,转步骤3.2.10,否则转步骤3.2.8;Step 3.2.7), calculate g(t 2 +8σ): u(t 2 +8σ)=u (1) +(t 2 +8σ)p, when g(t 2 +8σ)>0, then take t left =t 2 +4σ, t right =t 2 +8σ, go to step 3.2.10, otherwise go to step 3.2.8;

步骤3.2.8)、计算g(t2+12σ):u(t2+12σ)=u(1)+(t2+12σ)p,当g(t2+12σ)>0,则取tleft=t2+8σ,tright=t2+12σ,转步骤3.2.10,否则转步骤3.2.9;Step 3.2.8), calculate g(t 2 +12σ): u(t 2 +12σ)=u (1) +(t 2 +12σ)p, when g(t 2 +12σ)>0, then take t left =t 2 +8σ, t right =t 2 +12σ, go to step 3.2.10, otherwise go to step 3.2.9;

步骤3.2.9)、计算g(t2+16σ):u(t2+16σ)=u(1)+(t2+16σ)p,当g(t2+16σ)>0,则取tleft=t2+12σ,tright=t2+16σ,转步骤3.2.10,否则取t3=t2+16σ,转步骤3.3;Step 3.2.9), calculate g(t 2 +16σ): u(t 2 +16σ)=u (1) +(t 2 +16σ)p, when g(t 2 +16σ)>0, then take t left =t 2 +12σ, t right =t 2 +16σ, go to step 3.2.10, otherwise take t 3 =t 2 +16σ, go to step 3.3;

步骤3.2.10)、取计算g(tmiddle):u(tmiddle)=u(1)+tmiddle·p,当|g(tmiddle)|≤10-3,则取t3=tmiddle,转步骤3.3,否则转步骤3.2.11;Step 3.2.10), take Calculate g(t middle ): u(t middle )=u (1) +t middle ·p, when |g(t middle )|≤10 -3 , then take t 3 =t middle , go to step 3.3, otherwise go to Step 3.2.11;

步骤3.2.11)、当g(tmiddle)>0,则取tright=tmiddle,否则取tleft=tmiddle,当trightt-tleft≤10-2,则取转步骤3.3,否则转步骤3.2.10;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;

步骤3.3)、计算f(t3)和第三估计点u(3):u(3)=u(1)+t3p,当f(t3)≤fmin,则取fmin=f(t3),转步骤3.4;Step 3.3), calculate f(t 3 ) and the third estimated point u (3) : u (3) =u (1) +t 3 p, when f(t 3 )≤f min , then take f min = f(t 3 ), go to step 3.4;

步骤3.4)、对第二估计点u(2)所对应的参数t2设定新的初始值:取n=n+1,当有定义,则取转步骤3.2.1;否则保留步骤3.3中的fmin及其对应的转至步骤3.5;Step 3.4), set a new initial value for the parameter t2 corresponding to the second estimated point u (2) : 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;

步骤3.5)、决定信号源位置坐标的最终估计值:将步骤3.4保留的作为信号源位置坐标的最终估计值。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

图1为本发明流程图。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.

实施例1Example 1

已知一个平面上8个传感器装置,这8个传感器装置的位置坐标分别为:It is known that there are 8 sensor devices on a plane, and the position coordinates of these 8 sensor devices are:

uu 11 == 11.899811.8998 50.595750.5957 ,, uu 22 == 49.836449.8364 69.907769.9077 ,, uu 33 == 95.974495.9744 89.090389.0903 ,, uu 44 == 34.038634.0386 95.929195.9291 ,, uu 55 == 58.526858.5268 54.721654.7216 ,,

u 6 = 22.3812 13.8624 , u 7 = 75.1267 14.9294 , u 8 = 25.5095 25.7508 ; 信号源的真实位置为 u 0 = 84.0717 25.4282 . u 6 = 22.3812 13.8624 , u 7 = 75.1267 14.9294 , u 8 = 25.5095 25.7508 ; The real position of the signal source is u 0 = 84.0717 25.4282 .

基于多维标度子空间分析的三阶段到达时间差定位方法包含以下三个阶段: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:

步骤1.1)、生成测量标量乘积矩阵:Step 1.1), generate measurement scalar product matrix:

已知信号传播速度为c,归一化为c=1。已测量到信号源u0到各个传感器um的到达时间与信号源u0到第1个传感器u1的到达时间之差为: t ^ 51 = - 37.7748 , t ^ 61 = - 13.4115 , t ^ 71 = - 62.2653 , t ^ 81 = - 18.3493. 则信号源u0到各个传感器um的到达距离与信号源u0到第1个传感器u1的到达距离之差为: d ^ 31 = c · t ^ 31 = - 11.4547 , d ^ 41 = c · t ^ 41 = - 10.5034 , d ^ 51 = c · t ^ 51 = - 37.7748 , d ^ 61 = c · t ^ 61 = - 13.4115 , d ^ 71 = c · t ^ 71 = - 62.2653 , d ^ 81 = c · t ^ 81 = - 18.3493 ; 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 0 to each sensor u m and the arrival time from signal source u 0 to the first sensor u 1 is: t ^ 51 = - 37.7748 , t ^ 61 = - 13.4115 , t ^ 71 = - 62.2653 , t ^ 81 = - 18.3493. Then the difference between the arrival distance from the signal source u 0 to each sensor u m and the arrival distance from the signal source u 0 to the first sensor u 1 is: d ^ 31 = c · t ^ 31 = - 11.4547 , d ^ 41 = c · t ^ 41 = - 10.5034 , d ^ 51 = c · t ^ 51 = - 37.7748 , d ^ 61 = c · t ^ 61 = - 13.4115 , d ^ 71 = c · t ^ 71 = - 62.2653 , d ^ 81 = c · t ^ 81 = - 18.3493 ;

已知当m=2,...,8时的测量误差的方差同为的误差方差2σ2=2×0.32It is known that when m=2,...,8 The variance of the measurement error is the same as but The error variance 2σ 2 =2×0.3 2 ;

生成测量标量乘积矩阵得到: B ^ = 10 3 × 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 ; generate matrix of scalar products of measurements get: B ^ = 10 3 × 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 ;

步骤1.2)、对测量标量乘积矩阵进行修正: B 1 = 10 3 × 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 ; Step 1.2), correcting the measurement scalar product matrix: B 1 = 10 3 × 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 ;

步骤1.3)、计算第一估计点和一个距离估计值包括以下步骤:Step 1.3), calculating the first estimated point and a distance estimate comprises the following steps:

步骤1.3.1)、对B12IM作特征值分解得到: v 1 = 0.3192 0.2147 0.5923 0.3385 0.1424 0.4581 0.1955 0.3440 , v 2 = 0.3014 - 0.1219 - 0.6419 - 0.2215 - 0.0590 0.5053 0.1312 0.3965 , v 3 = 0.3658 0.3804 - 0.4021 0.5600 0.2731 - 0.2697 - 0.3141 - 0.0110 , v 4 = - 0.0861 0.4891 - 0.0201 - 0.4562 0.4610 - 0.3272 0.3761 0.2895 , v 5 = - 0.1346 - 0.1949 - 0.2087 0.3640 0.2960 0.1066 0.7038 - 0.4167 v 6 = - 0.1195 - 0.6444 0.0439 0.2027 0.3445 - 0.3160 - 0.0632 0.5521 , v 7 = - 0.2098 - 0.0480 0.0203 - 0.1848 0.6725 0.4493 - 0.4503 - 0.2492 , v 8 = 0.7679 - 0.3154 0.1709 - 0.3269 0.1749 - 0.2047 0.0495 - 0.3161 ; s1=-1.1334×104,s2=8.5322×103,s3=2.8357×103,s4=-66.9765,s5=31.8773,s6=5.7825×10-13,s7=-5.1324×10-13,s8=-1.1186×10-13;k1=5.3656×10-5,k2=1.0543×10-5,k3=1.9163×10-5,k4=0.4293,k5=1.3348;Step 1.3.1), perform eigenvalue decomposition on B 12 I M to get: v 1 = 0.3192 0.2147 0.5923 0.3385 0.1424 0.4581 0.1955 0.3440 , v 2 = 0.3014 - 0.1219 - 0.6419 - 0.2215 - 0.0590 0.5053 0.1312 0.3965 , v 3 = 0.3658 0.3804 - 0.4021 0.5600 0.2731 - 0.2697 - 0.3141 - 0.0110 , v 4 = - 0.0861 0.4891 - 0.0201 - 0.4562 0.4610 - 0.3272 0.3761 0.2895 , v 5 = - 0.1346 - 0.1949 - 0.2087 0.3640 0.2960 0.1066 0.7038 - 0.4167 v 6 = - 0.1195 - 0.6444 0.0439 0.2027 0.3445 - 0.3160 - 0.0632 0.5521 , v 7 = - 0.2098 - 0.0480 0.0203 - 0.1848 0.6725 0.4493 - 0.4503 - 0.2492 , v 8 = 0.7679 - 0.3154 0.1709 - 0.3269 0.1749 - 0.2047 0.0495 - 0.3161 ; s 1 =-1.1334×10 4 , s 2 =8.5322×10 3 , s 3 =2.8357×10 3 , s 4 =-66.9765, s 5 =31.8773, s 6 =5.7825×10 -13 , s 7 =-5.1324 ×10 -13 , s 8 =-1.1186×10 -13 ; k 1 =5.3656×10 -5 , k 2 =1.0543×10 -5 , k 3 =1.9163×10 -5 , k 4 =0.4293, k 5 = 1.3348;

步骤1.3.2)、由特征向量的线性组合得到一个系数向量: v = - 0.2165 - 0.0501 - 0.2872 0.2901 0.5931 0.0017 1.1008 - 0.4319 ; Step 1.3.2), obtain a coefficient vector by the linear combination of feature vectors: v = - 0.2165 - 0.0501 - 0.2872 0.2901 0.5931 0.0017 1.1008 - 0.4319 ;

步骤1.3.3)、以系数向量v作为加权系数,对位置坐标矩阵 x 1 ... x 8 y 1 ... y 8 的列向量作线性组合,由此计算出第一估计点 u ( 1 ) = 83.6662 25.5762 和一个距离估计值 Step 1.3.3), with the coefficient vector v as the weighting coefficient, the position coordinate matrix x 1 ... x 8 the y 1 ... the y 8 The column vectors of are linearly combined to calculate the first estimated point u ( 1 ) = 83.6662 25.5762 and a distance estimate

步骤1.4)、计算方向向量p:因为 | 1 M T v 5 | = 0.5155 > 10 - 10 , u f = 83.4118 25.8888 , p = - 0.6311 0.7757 : Step 1.4), calculate the direction vector p: because | 1 m T v 5 | = 0.5155 > 10 - 10 , so u f = 83.4118 25.8888 , p = - 0.6311 0.7757 :

所述的第二阶段利用方向向量p和距离估计值计算出第二估计点u(2)所对应的参数t2:a1=-64.7006,a0=50.7965,Δ=4.1354×103;因Δ>0,故 The second stage uses the direction vector p and the distance estimate Calculate the parameter t 2 corresponding to the second estimated point u (2) : a 1 =-64.7006, a 0 =50.7965, Δ=4.1354×10 3 ; since Δ>0, so

所述的第三阶段包括:The third phase described includes:

步骤3.1)、以第二估计点u(2)所对应的参数t2为初始值:取fmin=∞,n=1,则 t 2 = t 2 ( 1 ) = 129.0074 ; Step 3.1), taking the parameter t 2 corresponding to the second estimated point u (2) as the initial value: take f min =∞, n=1, then t 2 = t 2 ( 1 ) = 129.0074 ;

步骤3.2)、根据参数t2的初始值,利用二分法求根g(t)=0,找出第三估计点u(3)所对应的参数t3Step 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 (3) :

步骤3.2.1)、计算 g ( t 2 ) : u ( t 2 ) = 2.2467 125.6451 , g(t2)=-2.0785×105,因为g(t2)≤0,故转步骤3.2.6;Step 3.2.1), calculation g ( t 2 ) : u ( t 2 ) = 2.2467 125.6451 , g(t 2 )=-2.0785×10 5 , since g(t 2 )≤0, go to step 3.2.6;

步骤3.2.6)、计算 g ( t 2 + 4 σ ) : u ( t 2 + 4 σ ) = 1.4893 126.5759 , g(t2+4σ)=-2.0240×105,因为g(t2+4σ)≤0,故转步骤3.2.7;Step 3.2.6), calculation g ( t 2 + 4 σ ) : u ( t 2 + 4 σ ) = 1.4893 126.5759 , g(t 2 +4σ)=-2.0240×10 5 , since g(t 2 +4σ)≤0, go to step 3.2.7;

步骤3.2.7)、计算 g ( t 2 + 8 σ ) : u ( t 2 + 8 σ ) = 0.7320 127.5067 , g(t2+8σ)=-1.9712×105,因为g(t2+8σ)≤0,故转步骤3.2.8;Step 3.2.7), calculation g ( t 2 + 8 σ ) : u ( t 2 + 8 σ ) = 0.7320 127.5067 , g(t 2 +8σ)=-1.9712×10 5 , since g(t 2 +8σ)≤0, go to step 3.2.8;

步骤3.2.8)、计算 g ( t 2 + 12 σ ) : u ( t 2 + 12 σ ) = - 0.0254 128.4376 , g(t2+12σ)=-1.9201×105,因为g(t2+12σ)≤0,故转步骤3.2.9;Step 3.2.8), calculation g ( t 2 + 12 σ ) : u ( t 2 + 12 σ ) = - 0.0254 128.4376 , g(t 2 +12σ)=-1.9201×10 5 , since g(t 2 +12σ)≤0, go to step 3.2.9;

步骤3.2.9)、计算 g ( t 2 + 16 σ ) : u ( t 2 + 16 σ ) = - 0.7827 129.3684 , g(t2+16σ)=-1.8707×105,因为g(t2+16σ)≤0,故取t3=t2+16σ=133.8074,转步骤3.3;Step 3.2.9), calculation g ( t 2 + 16 σ ) : u ( t 2 + 16 σ ) = - 0.7827 129.3684 , g(t 2 +16σ)=-1.8707×10 5 , because g(t 2 +16σ)≤0, so take t 3 =t 2 +16σ=133.8074, go to step 3.3;

步骤3.3)、计算第三估计点 u ( 3 ) : u ( 3 ) = u ( 1 ) + t 3 p = - 0.7827 129.3684 , f(t3)=2.9857×107,由于f(t3)≤fmin,则取 u ^ = u ( 3 ) = - 0.7827 129.3684 , 并取fmin=2.9857×107,转步骤3.4;Step 3.3), calculate the third estimation point u ( 3 ) : u ( 3 ) = u ( 1 ) + t 3 p = - 0.7827 129.3684 , f(t 3 )=2.9857×10 7 , since f(t 3 )≤f min , then take u ^ = u ( 3 ) = - 0.7827 129.3684 , And take f min =2.9857×10 7 , go to step 3.4;

步骤3.4)、对第二估计点u(2)所对应的参数t2设定新的初始值:取n=n+1=2, t 2 = t 2 ( 2 ) = 0.3937 , 转步骤3.2.1;Step 3.4), setting a new initial value for the parameter t2 corresponding to the second estimation point u (2) : getting n=n+1=2, t 2 = t 2 ( 2 ) = 0.3937 , Go to step 3.2.1;

步骤3.2.1)、计算 g ( t 2 ) : u ( t 2 ) = 83.4177 25.8816 , g(t2)=2.9567×104,因为g(t2)>0,故转步骤3.2.2;Step 3.2.1), calculation g ( t 2 ) : u ( t 2 ) = 83.4177 25.8816 , g(t 2 )=2.9567×10 4 , since g(t 2 )>0, go to step 3.2.2;

步骤3.2.2)、计算 g ( t 2 - 4 σ ) : u ( t 2 - 4 σ ) = 84.1750 24.9508 , g(t2-4σ)=-6.6714×104,因为g(t2-4σ)≤0,故取tleft=-0.8063,tright=0.3937,转步骤3.2.10;Step 3.2.2), calculation g ( t 2 - 4 σ ) : u ( t 2 - 4 σ ) = 84.1750 24.9508 , g(t 2 -4σ)=-6.6714×10 4 , because g(t 2 -4σ)≤0, so take t left =-0.8063, t right =0.3937, go to step 3.2.10;

步骤3.2.10)、计算|g(tmiddle)|:tmiddle=-0.2063, u ( t m i d d l e ) = 83.7963 25.4162 , g(tmiddle)=-2.0526×104,因为|g(tmiddle)|>10-3,故转步骤3.2.11;Step 3.2.10), calculating |g(t middle )|: t middle =-0.2063, u ( t m i d d l e ) = 83.7963 25.4162 , g(t middle )=-2.0526×10 4 , because |g(t middle )|>10 -3 , so turn to step 3.2.11;

步骤3.2.11)、计算trightt-tleft:因为g(tmiddle)≤0,故取tleft=-0.2063,因为trightt-tleft>10-2,故转步骤3.2.10;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;

步骤3.2.10)、计算|g(tmiddle)|:tmiddle=0.0937, u ( t m i d d l e ) = 83.6070 25.6489 , g(tmiddle)=4.0541×103,因为|g(tmiddle)|>10-3,故转步骤3.2.11;Step 3.2.10), calculating |g(t middle )|: t middle =0.0937, u ( t m i d d l e ) = 83.6070 25.6489 , g(t middle )=4.0541×10 3 , because |g(t middle )|>10 -3 , so turn to step 3.2.11;

步骤3.2.11)、计算trightt-tleft:因为g(tmiddle)>0,故取tright=0.0937,因为trightt-tleft>10-2,故转步骤3.2.10;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;

步骤3.2.10)、计算tmiddle=-0.0563, u ( t m i d d l e ) = 83.7017 25.5326 , g(tmiddle)=-8.3555×103,因为|g(tmiddle)|>10-3,故转步骤3.2.11;Step 3.2.10), calculate t middle =-0.0563, u ( t m i d d l e ) = 83.7017 25.5326 , g(t middle )=-8.3555×10 3 , because |g(t middle )|>10 -3 , so turn to step 3.2.11;

步骤3.2.11)、计算trightt-tleft:因为g(tmiddle)≤0,故取tleft=-0.0563,因为trightt-tleft>10-2,故转步骤3.2.10;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;

步骤3.2.10)、计算|g(tmiddle)|:tmiddle=0.0187, u ( t m i d d l e ) = 83.6543 25.5907 , g(tmiddle)=-2.1802×103,因为|g(tmiddle)|>10-3,故转步骤3.2.11;Step 3.2.10), calculate |g(t middle )|: t middle =0.0187, u ( t m i d d l e ) = 83.6543 25.5907 , g(t middle )=-2.1802×10 3 , because |g(t middle )|>10 -3 , so turn to step 3.2.11;

步骤3.2.11)、计算trightt-tleft:因为g(tmiddle)≤0,故取tleft=0.0187,因为trightt-tleft>10-2,故转步骤3.2.10;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;

步骤3.2.10)、计算|g(tmiddle)|:tmiddle=0.0562, u ( t m i d d l e ) = 83.6307 25.6198 , g(tmiddle)=929.6178,因为|g(tmiddle)|>10-3,故转步骤3.2.11;Step 3.2.10), calculate |g(t middle )|: t middle =0.0562, u ( t m i d d l e ) = 83.6307 25.6198 , g(t middle )=929.6178, because |g(t middle )|>10 -3 , so turn to step 3.2.11;

步骤3.2.11)、计算trightt-tleft:因为g(tmiddle)>0,故取tright=0.0562,因为trightt-tleft>10-2,故转步骤3.2.10;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;

步骤3.2.10)、计算|g(tmiddle)|:tmiddle=0.0375, u ( t m i d d l e ) = 83.6425 25.6053 , g(tmiddle)=-627.1421,因为|g(tmiddle)|>10-3,故转步骤3.2.11;Step 3.2.10), calculating |g(t middle )|: t middle =0.0375, u ( t m i d d l e ) = 83.6425 25.6053 , g(t middle )=-627.1421, because |g(t middle )|>10 -3 , so turn to step 3.2.11;

步骤3.2.11)、计算trightt-tleft:因为g(tmiddle)≤0,故取tleft=0.0375,则trightt-tleft>10-2,故转步骤3.2.10;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;

步骤3.2.10)、计算|g(tmiddle)|:tmiddle=0.0469, u ( t m i d d l e ) = 83.6366 25.6126 , g(tmiddle)=150.7786,因为|g(tmiddle)|>10-3,故转步骤3.2.11;Step 3.2.10), calculating |g(t middle )|: t middle =0.0469, u ( t m i d d l e ) = 83.6366 25.6126 , g(t middle )=150.7786, because |g(t middle )|>10 -3 , so turn to step 3.2.11;

步骤3.2.11)、计算trightt-tleft:因为g(tmiddle)>0,故取tright=0.0469,则trightt-tleft≤10-2,故取t3=0.0422,转步骤3.3;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 3 =0.0422, go to step 3.3 ;

步骤3.3)、计算第三估计点 u ( 3 ) : u ( 3 ) = 83.6395 25.6089 , 计算f(t3)=1.0242×104,因为f(t3)≤fmin,故取 u ^ = u ( 3 ) = 83.6395 25.6089 , 并取fmin=1.0242×104,转步骤3.4;Step 3.3), calculate the third estimation point u ( 3 ) : u ( 3 ) = 83.6395 25.6089 , Calculate f(t 3 )=1.0242×10 4 , because f(t 3 )≤f min , so take u ^ = u ( 3 ) = 83.6395 25.6089 , And take f min =1.0242×10 4 , go to step 3.4;

步骤3.4)、以第二估计点u(2)所对应的参数t2为初始值:取n=n+1=3,因为无的定义,故保留步骤3.3中的fmin=1.0242×104及其对应的 u ^ = 83.6395 25.6089 , 转步骤3.5;Step 3.4), take the parameter t2 corresponding to the second estimated point u (2) as the initial value: get n=n+1=3, because there is no The definition of , so keep f min =1.0242×10 4 and its corresponding value in step 3.3 u ^ = 83.6395 25.6089 , Go to step 3.5;

步骤3.5)、决定信号源位置坐标的最终估计值:fmin=1.0242×104时, u ^ = 83.6395 25.6089 为信号源位置坐标的最终估计值。Step 3.5), determine the final estimated value of the position coordinates of the signal source: when f min =1.0242×10 4 , u ^ = 83.6395 25.6089 is the final estimated value of the position coordinates of the signal source.

经本实施例方式计算获得的信号源位置坐标的最终估计值 u ^ = 83.6395 25.6089 与信号源位置真实坐标的 u 0 = 84.0717 25.4282 基本相符,误差极小。The final estimated value of the position coordinates of the signal source calculated through the method of this embodiment u ^ = 83.6395 25.6089 and the real coordinates of the signal source position u 0 = 84.0717 25.4282 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 0to each sensor u mtime of arrival and signal source u 0to the 1st sensor u 1the 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 [ B ^ ] i , j = 1 2 ( d ^ i 1 - d ^ j 1 ) 2 - 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 12i mmake Eigenvalues Decomposition: B 12i m=[v 1..., v m] diag (s 1..., s m) [v 1..., v m] t, wherein: v 1..., v m∈ R mthe vector that pairwise orthogonal and mould are 1, diag (s 1..., s m) expression diagonal element is s 1..., s mdiagonal matrix, [v 1..., v m] tthe transposition of representing matrix, to s 1..., s mdescending sort is carried out according to absolute value, namely | s 1|>=...>=| s m|;
Step 1.3.2) obtain coefficient vector v:v by the linear combination of proper vector 1..., v mfor R mone group of orthonormal basis, therefore R min coefficient vector v can be expressed as v 1..., v mlinear combination v=k 1v 1+ ...+k 5v 5, work as m=1 ..., when 5, base v mcombination coefficient and when m>=6, base v mcombination coefficient k 6..., k mbe 0;
Step 1.3.3) using coefficient vector v as weighting coefficient, to position coordinates matrix x 1 ... x M y 1 ... y M Column vector do linear combination, calculate the first estimation point thus u ( 1 ) = x 1 ... x M y 1 ... y M v 1 M T v With a range estimation d 1 ( 1 ) = - d ^ 11 ... d ^ M 1 v 1 M T v ;
Step 1.4) calculated direction vector p: when then otherwise p = [ - ( y 2 - y 1 ) , x 2 - x 1 ] T | | u 2 - u 1 | | , Wherein: u f = x 1 ... x M y 1 ... y M 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 2:
Definition t 2 = t 2 ( n ) , a 0 = | | u ( 1 ) - u 1 | | 2 - ( d 1 ( 1 ) ) 2 , a 1=p T(u (1)-u 1), Δ = a 1 2 - a 0 ; Then as Δ >0, define during n=1 t 2 ( 1 ) = - a 1 + Δ , Define during n=2 t 2 ( 2 ) = - a 1 - Δ ; Otherwise, define during n=1 t 2 ( 1 ) = - 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 (2)corresponding parametric t 2for initial value: definition f min=∞, n=1,
Step 3.2) process of dichotomy rooting g (t)=0, find out the 3rd estimation point u (3)corresponding parametric t 3, 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, d · m = p T w m / | | w m | | , Further order Z = x 1 - x u ... x M - x u y 1 - y u ... y M - y u d 1 ... d M , B=Z T·diag(1,1,-1)·Z, Z · = - l x ... - l x - l y ... - l y d · 1 ... d · M , B=Z T·diag(1,1,-1)·Z,f(t)=tr((B-B 1) 2), g ( t ) = 2 t r ( ( B - B 1 ) B · ) , Wherein: the mark of tr () representing matrix;
Step 3.3) calculate f (t 3) and the 3rd estimation point u (3): u (3)=u (1)+ t 3p, as f (t 3)≤f min, then get f min=f (t 3), go to step 3.4;
Step 3.4) with the second estimation point u (2)corresponding parametric t 2for 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 0actual coordinate.
CN201510710744.7A 2015-10-28 2015-10-28 The three stage reaching time-difference localization methods based on multidimensional scaling subspace analysis Expired - Fee Related CN105259533B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201510710744.7A CN105259533B (en) 2015-10-28 2015-10-28 The three stage reaching time-difference localization methods based on multidimensional scaling subspace analysis

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201510710744.7A CN105259533B (en) 2015-10-28 2015-10-28 The three stage reaching time-difference localization methods based on multidimensional scaling subspace analysis

Publications (2)

Publication Number Publication Date
CN105259533A true CN105259533A (en) 2016-01-20
CN105259533B CN105259533B (en) 2017-07-18

Family

ID=55099294

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201510710744.7A Expired - Fee Related CN105259533B (en) 2015-10-28 2015-10-28 The three stage reaching time-difference localization methods based on multidimensional scaling subspace analysis

Country Status (1)

Country Link
CN (1) CN105259533B (en)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105891776A (en) * 2016-04-06 2016-08-24 上海交通大学 Direct method time difference of arrival positioning method based on multidimensional scaling (MDS) model
CN108279411A (en) * 2018-02-01 2018-07-13 电子科技大学 A kind of passive MIMO time difference positioning methods based on MDS
CN111551897A (en) * 2020-04-25 2020-08-18 中国人民解放军战略支援部队信息工程大学 TDOA localization method based on weighted multidimensional scaling and polynomial root finding in the presence of sensor location prior observation errors

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2009010832A1 (en) * 2007-07-18 2009-01-22 Bang & Olufsen A/S Loudspeaker position estimation
CN101354435A (en) * 2008-09-05 2009-01-28 清华大学 Self-localization method of sensor network nodes based on distance order relationship
US20110140969A1 (en) * 2006-04-19 2011-06-16 Mustafa Ergen Method And System For Hybrid Positioning Using Partial Distance Information
CN104038901A (en) * 2014-05-30 2014-09-10 中南大学 Indoor positioning method for reducing fingerprint data acquisition workload

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20110140969A1 (en) * 2006-04-19 2011-06-16 Mustafa Ergen Method And System For Hybrid Positioning Using Partial Distance Information
WO2009010832A1 (en) * 2007-07-18 2009-01-22 Bang & Olufsen A/S Loudspeaker position estimation
CN101354435A (en) * 2008-09-05 2009-01-28 清华大学 Self-localization method of sensor network nodes based on distance order relationship
CN104038901A (en) * 2014-05-30 2014-09-10 中南大学 Indoor positioning method for reducing fingerprint data acquisition workload

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
HE-WEN WEI等: "multidimensional scaling analysis for passive moving target localization with TDOA and FDOA measurement", 《IEEE TRANSACTIONS ON SIGNLA PROCESSING》 *
彭鑫: "无线传感器网络中基于多维标度的节点定位算法", 《中国优秀硕士学位论文全文数据库 信息科技辑》 *
王琳等: "一种多维尺度分析到达时间差定位算法", 《导航定位学报》 *

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105891776A (en) * 2016-04-06 2016-08-24 上海交通大学 Direct method time difference of arrival positioning method based on multidimensional scaling (MDS) model
CN105891776B (en) * 2016-04-06 2018-06-12 上海交通大学 Direct method reaching time-difference localization method based on MDS models
CN108279411A (en) * 2018-02-01 2018-07-13 电子科技大学 A kind of passive MIMO time difference positioning methods based on MDS
CN108279411B (en) * 2018-02-01 2020-04-14 电子科技大学 A Passive MIMO Time Difference Location Method Based on MDS
CN111551897A (en) * 2020-04-25 2020-08-18 中国人民解放军战略支援部队信息工程大学 TDOA localization method based on weighted multidimensional scaling and polynomial root finding in the presence of sensor location prior observation errors
CN111551897B (en) * 2020-04-25 2021-01-22 中国人民解放军战略支援部队信息工程大学 TDOA localization method based on weighted multidimensional scaling and polynomial root finding under sensor position error

Also Published As

Publication number Publication date
CN105259533B (en) 2017-07-18

Similar Documents

Publication Publication Date Title
CN105866735B (en) The reaching time-difference iteration localization method of amendment cost function based on MDS models
CN105954712B (en) The direct localization method of the multiple target of associated wireless electric signal complex envelope and carrier phase information
CN106793087B (en) Array antenna indoor positioning method based on AOA and PDOA
CN104853435B (en) A kind of indoor orientation method based on probability and device
CN101221238B (en) Dynamic deviation estimation method based on gauss average value mobile registration
CN105158730B (en) TDOA localization methods based on MDS subspaces the 4th and the 5th characteristic vector
CN106405533B (en) Radar target combined synchronization and localization method based on constraint weighted least-squares
CN104619020A (en) WIFI Indoor Positioning Method Based on RSSI and TOA Ranging
CN104038901B (en) Indoor positioning method for reducing fingerprint data acquisition workload
CN103369466B (en) A kind of map match assists indoor orientation method
CN107124762B (en) Wireless positioning method for efficiently eliminating non-line-of-sight errors
CN109195110B (en) Indoor localization method based on hierarchical clustering technology and online extreme learning machine
CN104684081A (en) Node localization algorithm for selecting anchor nodes based on distance clustering in wireless sensor networks
CN108051779A (en) A kind of positioning node preferred method towards TDOA
CN106353720A (en) Multi-station continuous positioning model based on TDOA/GROA (time different of arrival/gain ratio of arrival)
CN111551897A (en) TDOA localization method based on weighted multidimensional scaling and polynomial root finding in the presence of sensor location prior observation errors
CN104581945B (en) The WLAN indoor orientation methods of semi-supervised APC clustering algorithms based on distance restraint
CN102200573A (en) Method for determining incoming wave direction of near-field target signal
CN105259533A (en) Three-stage arrival time difference positioning method based on multidimensional scaling sub space analysis
CN110673196A (en) Time difference positioning method based on multidimensional calibration and polynomial root finding
CN113923590B (en) A TOA positioning method under the condition of uncertain anchor node position
CN105891776B (en) Direct method reaching time-difference localization method based on MDS models
CN114325581A (en) Elliptical target positioning method with clock synchronization error
CN104635206B (en) A kind of method and device of wireless location
WO2017049914A1 (en) Terminal positioning method, apparatus, and system

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20170718

Termination date: 20191028