CN112986907B - 一种时钟偏差和时钟漂移条件下的运动目标定位方法 - Google Patents

Abstract

本发明公开了一种时钟偏差和时钟漂移条件下的运动目标定位方法，其构建包含源节点的位置和速度、时钟偏差和时钟漂移四个未知量的加权最小二乘WLS问题，然后采用半正定松弛方法将非凸的加权最小二乘WLS问题转化为凸的半正定规划问题，再利用内点法求解凸的半正定规划问题，以得到全局最优解；优点是其在时钟偏差和时钟漂移同时存在且未知的情况下，从加权最小二乘估计模型估计理论和凸优化理论的角度出发，在实现源节点的位置和速度估计的同时，实现时钟偏差和时钟漂移的联合估计，且定位精度高、计算复杂度低，适用于复杂部署环境。

Description

一种时钟偏差和时钟漂移条件下的运动目标定位方法

技术领域

本发明涉及一种目标定位技术，尤其是涉及一种时钟偏差和时钟漂移条件下的运动目标定位方法。

背景技术

近年来，随着无线传感器网络技术的快速发展，其目标定位技术在军事、工业和民用领域得到了广泛的应用。现有的目标定位技术中，根据信号测量方式的不同，可以分为到达时间(time-of-arrival，TOA)、到达频率(frequency-of-arrival，FOA)、到达时间差(time-difference-of-arrival，TDOA)、到达频率差(frequency-difference-of-arrival，FDOA)、到达角度(angle-of-arrival，AOA)和接收信号强度(received-signal-strength，RSS)等。一般来看，基于TOA的目标定位方法和基于TDOA的目标定位方法能够达到较高的定位精度。如果定位目标与传感器节点之间存在相对运动，则可以通过引入FOA或者FDOA的观测信号，能够在估计定位目标运动速度的基础上，进一步提升定位目标的定位精度。

在无线传感器网络中，源节点(即定位目标)和传感器节点之间时钟的精确同步是进行基于TOA的目标定位的必要条件，这也是基于TOA的目标定位实用化过程中遇到的挑战之一。在使用无线电波作为载体进行定位时，10纳秒的时钟偏差将会带来3米的误差(无线电波在空气中的传播速度为3×10⁸m/s)。在文献“Second-order cone relaxation forTOA-based source localization with unknown start transmission time(未知起始传输时间条件下的二阶锥松弛TOA定位方法)”(IEEE Trans.Veh.Technol.,vol.63,no.6,pp.2973-2977,2014)中，Wang等作者提出了一种联合估计时钟偏差和源节点位置的二阶锥松弛方法，实验结果表明，该方法比现有的最小-最大算法(min-max algorithm，MMA)具有更好的估计性能，然而，该方法仅针对源节点和传感器节点是静态的应用场景。对于定位目标与传感器节点之间存在相对运动的定位场景，在文献“Passive source localizationusing importance sampling based on TOA and FOA measurements(基于TOA和FOA观测信号的重要性采样无源定位方法)”(Frontiers of Information Technology&ElectronicEngineering,vol.18,no.8,pp.1167-1179,2017)中，Liu等作者提出了一种基于TOA和FOA观测信号的蒙特卡罗重要性采样(Monte Carlo Importance Sampling)定位方法，该方法能够同时估计源节点的位置和速度，同时在低噪声区间达到克拉美-罗界(Cramer-Raolower bound，CRLB)，但是该方法并未考虑时钟的同步问题。在已公告的中国发明专利“抑制传感器位置速度先验误差的加权多维标度TOA和FOA多源协同定位方法”(专利号：ZL202010335968.5)中，王鼎等人提出了一种抑制传感器位置速度先验误差的加权多维标度TOA和FOA多源协同定位方法，该方法利用线性最小二乘(Linear Least Squares，LLS)估计优化模型获得各个辐射源的位置向量和速度向量的估计值，虽然该方法考虑了多个源节点的定位问题，但是同样存在未考虑时钟同步的缺陷。综合现有的文献来看，利用TOA和FOA观测信号进行源节点的位置和速度估计时，对于同时存在时钟偏差和时钟漂移的场景，尚未有可行的解决方案。

发明内容

本发明所要解决的技术问题是提供一种时钟偏差和时钟漂移条件下的运动目标定位方法，其在时钟偏差(Clock Bias)和时钟漂移(Clock Drift)同时存在且未知的情况下，从加权最小二乘(Weighted Least Squares，WLS)估计模型估计理论和凸优化理论的角度出发，在实现源节点的位置和速度估计的同时，实现时钟偏差和时钟漂移的联合估计，且定位精度高、计算复杂度低，适用于复杂部署环境。

本发明解决上述技术问题所采用的技术方案为：一种时钟偏差和时钟漂移条件下的运动目标定位方法，其特征在于包括以下步骤：

步骤1：在一个平面空间或者立体空间中部署无线传感器网络，该无线传感器网络中存在1个位置和速度未知的用于发射无线信号的源节点、N个位置和速度已知的用于接收无线信号的传感器节点、1个用于估计源节点的位置和速度及源节点和传感器节点之间的时钟偏差和时钟漂移的中心节点，将源节点的位置记为x^o，将源节点的速度记为v^o，将N个传感器节点的位置分别记为s₁,s₂,…,s_i,…,s_N，将N个传感器节点的速度分别记为

其中，N为正整数，若无线传感器网络部署在平面空间中则取N≥5，若无线传感器网络部署在立体空间中则取N≥6，s₁表示第1个传感器节点的位置，s₂表示第2个传感器节点的位置，s_i表示第i个传感器节点的位置，s_N表示第N个传感器节点的位置,

表示第1个传感器节点的速度，

表示第2个传感器节点的速度，

表示第i个传感器节点的速度，

表示第N个传感器节点的速度，i为正整数，1≤i≤N；

步骤2：根据到达时间TOA观测模型，获取每个传感器节点接收到的无线信号的到达时间，将第i个传感器节点接收到的无线信号的到达时间记为τ_i，τ_i描述为：

并根据到达频率FOA观测模型，获取每个传感器节点接收到的无线信号的到达频率，将第i个传感器节点接收到的无线信号的到达频率记为f_i，f_i描述为：

其中，符号“||||”为求欧氏距离符号，c表示无线信号的传播速度，

表示未知的时钟偏差，f₀表示载波频率，上标符号“T”表示向量或矩阵的转置，

表示未知的时钟漂移，Δτ_i为TOA测量噪声，Δf_i为FOA测量噪声；

步骤3：每个传感器节点将其接收到的无线信号的到达时间和到达频率发送给中心节点；

步骤4：中心节点将每个传感器节点接收到的无线信号的到达时间转变为距离观测量，将第i个传感器节点接收到的无线信号的距离观测量记为d_i，d_i描述为：

同样，中心节点将每个传感器节点接收到的无线信号的到达频率转变为距离变化率观测量，将第i个传感器节点接收到的无线信号的距离变化率观测量记为b_i，b_i描述为：

其中，

n_i＝Δτ_ic，

然后令d表示距离观测向量，令b表示距离变化率观测向量，令d^o表示距离真实向量，令b^o表示距离变化率真实向量，d＝[d₁,…,d_i,…,d_N]^T，b＝[b₁,…,b_i,…,b_N]^T，

其中，d₁表示第1个传感器节点接收到的无线信号的距离观测量，d_N表示第N个传感器节点接收到的无线信号的距离观测量，b₁表示第1个传感器节点接收到的无线信号的距离变化率观测量，b_N表示第N个传感器节点接收到的无线信号的距离变化率观测量，

和

根据

计算得到，

和

根据

计算得到；并引入两个噪声向量，分别记为n和m，n＝[n₁,…,n_i,…,n_N]^T，m＝[m₁,…,m_i,…,m_N]^T，其中，n和m相互独立，n服从均值为0的高斯分布，且协方差矩阵Q_n＝E[nn^T]，m服从均值为0的高斯分布，且协方差矩阵Q_m＝E[mm^T]，E[]为求数学期望，n₁和n_N根据n_i＝Δτ_ic计算得到，m₁和m_N根据

计算得到；

步骤5：对

进行移项，变换为

然后对

两边进行平方操作，并忽略噪声平方项

得到

并将

代入

中，得到

再展开

后进行移项，并忽略二阶噪声项n_im_i，得到

步骤6：引入辅助向量y^o，y^o定义为

然后将

改写成矩阵形式，描述为：A₁y^o-h₁≈B₁n；并将

改写成矩阵形式，描述为：A₂y^o-h₂≈B₂n+D₂m；接着联立A₁y^o-h₁≈B₁n和A₂y^o-h₂≈B₂n+D₂m得到Ay^o-h≈Bn+Dm；再根据Ay^o-h≈Bn+Dm构建加权最小二乘WLS问题，描述为：

其中，

A₁的维数为N×(2k+6)，0_1×k表示维数为1×k且元素全为0的列向量，若无线传感器网络部署在平面空间中则k＝2，若无线传感器网络部署在立体空间中则k＝3，

h₁的维数为N×1，

表示将

作为对角元素组成对角矩阵，

A₂的维数为N×(2k+6)，

h₂的维数为N×1，

表示将

作为对角元素组成对角矩阵，

表示将

作为对角元素组成对角矩阵，

A的维数为2N×(2k+6)，

h的维数为2N×1，

B的维数为2N×N，

D的维数为2N×N，0_N×N表示维数为N×N且元素全为0的矩阵，

表示求使得(Ay-h)^TR^-1(Ay-h)的值最小时y的值，y为y^o对应的优化变量，R表示权重矩阵，R＝BQ_nB^T+DQ_mD^T＝[BD]Q[B D]^T，R^-1表示R的逆，[B D]为维数为2N×2N的矩阵；

步骤7：通过发掘y^o对应的优化变量y中的元素之间隐含的关系，将加权最小二乘WLS问题进一步描述为：

其中，“s.t.”表示“受约束于……”，y(1:k)表示y中的第1个元素至第k个元素所组成的向量，y(k+1:2k)表示y中的第(k+1)个元素至第2k个元素所组成的向量，y(2k+1)表示y中的第(2k+1)个元素，y(2k+2)表示y中的第(2k+2)个元素，y(2k+3)表示y中的第(2k+3)个元素，y(2k+4)表示y中的第(2k+4)个元素，y(2k+5)表示y中的第(2k+5)个元素，y(2k+6)表示y中的第(2k+6)个元素；

步骤8：引入矩阵Y＝yy^T，将加权最小二乘WLS问题

松弛为

其中，

表示求使得

的值最小时Y和y的值，tr()表示求矩阵的迹，

表示

是半正定的，rank()表示求矩阵的秩，

和

源于如下等价公式：

步骤9：将半正定松弛后的问题

中的约束条件||y(1:k)||²＝y(2k+3)重写为tr(Y(1:k,1:k))＝y(2k+3)、y^T(1:k)y(k+1:2k)＝y(2k+4)重写为tr(Y(1:k,k+1:2k))＝y(2k+4)、y²(2k+1)＝y(2k+5)重写为Y(2k+1,2k+1)＝y(2k+5)、y(2k+1)y(2k+2)＝y(2k+6)重写为Y(2k+1,2k+2)＝y(2k+6)、y(2k+1)y(2k+6)＝y(2k+2)y(2k+5)重写为Y(2k+1,2k+6)＝y(2k+2,2k+5)；其中，Y(1:k,1:k)表示Y中的第1行至第k行、第1列至第k列的元素所组成的矩阵，Y(1:k,k+1:2k)表示Y中的第1行至第k行、第(k+1)列至第2k列的元素所组成的矩阵，Y(2k+1,2k+1)表示Y中的第2k+1行第2k+1列的元素，Y(2k+1,2k+2)表示Y中的第2k+1行第2k+2列的元素，Y(2k+1,2k+6)表示Y中的第2k+1行第2k+6列的元素；

步骤10：丢弃

中的

并结合tr(Y(1:k,1:k))＝y(2k+3)、tr(Y(1:k,k+1:2k))＝y(2k+4)、Y(2k+1,2k+1)＝y(2k+5)、Y(2k+1,2k+2)＝y(2k+6)、Y(2k+1,2k+6)＝y(2k+2,2k+5)，得到半正定规划问题，描述为：

步骤11：令R的初始值为维数为2N×2N的单位矩阵I_2N；然后根据

计算得到F的值；接着利用内点法求解半正定规划问题

得到Y的初步全局最优解和y的初步全局最优解，对应记为Y^*1和y^*1；再根据y的定义及y^*1，得到

各自的初步估计值，对应记为

其中，y^*1(2k+1)表示y^*1中的第(2k+1)个元素，y^*1(2k+2)表示y^*1中的第(2k+2)个元素；

步骤12：根据

并结合R＝BQ_nB^T+DQ_mD^T＝[B D]Q[B D]^T，更新R的值，进而更新F的值；接着再次利用内点法求解半正定规划问题

得到Y的最终全局最优解和y的最终全局最优解，对应记为Y^*2和y^*2；再根据y的定义及y^*2，得到x^o、v^o、

各自的最终估计值，对应记为x^*2、v^*2、

x^*2＝y^*2(1:k)，v^*2＝y^*2(k+1:2k)，

再根据

和

和

得到

和

的最终估计值，对应记为

和

其中，y^*2(1:k)表示y^*2中的第1个元素至第k个元素所组成的向量，y^*2(k+1:2k)表示y^*2中的第(k+1)个元素至第2k个元素所组成的向量，y^*2(2k+1)表示y^*2中的第(2k+1)个元素，y^*2(2k+2)表示y^*2中的第(2k+2)个元素。

与现有技术相比，本发明的优点在于：

1)本发明方法利用半正定松弛技术将非凸的加权最小二乘WLS问题转化为凸的半正定规划问题，再利用内点法求解凸的半正定规划问题，这样能够获得全局最优解，从而克服了传统的最大似然ML估计方法陷入局部最优点和定位精度较低的缺点。

2)本发明方法通过引入包含源节点的位置、速度、时钟偏差和时钟漂移的复合向量，并从简洁度、平衡性等多个角度考虑，设计其变量的顺序和组合，使得松弛后的半正定规划问题尽量的紧，从而省去了因半正定规划问题不够紧而带来的后续操作，如使用高斯随机化、本地搜索等方法带来的复杂度提升。

3)本发明方法无需将时钟偏差和时钟漂移作为已知条件，直接通过接收到的无线信号就可以恢复源节点的位置、速度、时钟偏差和时钟漂移，且计算复杂度低。

4)本发明方法在时钟偏差和时钟漂移同时存在且未知的情况下，在实现源节点的位置和速度估计的同时，能够实现时钟偏差和时钟漂移的联合估计，适用于复杂部署环境。

附图说明

图1为本发明方法的总体实现框图；

图2为传感器节点的数量为5个，

的变化范围从10^-2m²变化至10³m²时，本发明方法(简称SDP)、ML方法、Raw WLS方法的位置估计的RMSE性能随

变化的曲线；

图3为传感器节点的数量为5个，

的变化范围从10^-2m²变化至10³m²时，本发明方法(简称SDP)、ML方法、Raw WLS方法的速度估计的RMSE性能随

变化的曲线；

图4为传感器节点的数量为5个，

的变化范围从10^-2m²变化至10³m²时，本发明方法(简称SDP)、ML方法的时钟偏差估计的RMSE性能随

变化的曲线；

图5为传感器节点的数量为5个，

的变化范围从10^-2m²变化至10³m²时，本发明方法(简称SDP)、ML方法的时钟漂移估计的RMSE性能随

变化的曲线；

图6为

固定为10m²，传感器的数量从4个增加至9个时，本发明方法(简称SDP)、ML方法、Raw WLS方法的位置估计的RMSE性能随传感器节点的数量变化的曲线；

图7为

固定为10m²，传感器的数量从4个增加至9个时，本发明方法(简称SDP)、ML方法、Raw WLS方法的速度估计的RMSE性能随传感器节点的数量变化的曲线；

图8为传感器节点的数量为5个，

固定为10m²，FOA测量噪声系数η从10^-3变化至10^-1时，本发明方法、TOA-SOCP方法的位置估计的RMSE性能随FOA测量噪声系数变化的曲线。

具体实施方式

以下结合附图实施例对本发明作进一步详细描述。

本发明提出的一种时钟偏差和时钟漂移条件下的运动目标定位方法，其总体实现框图如图1所示，其包括以下步骤：

步骤1：在一个平面空间或者立体空间中部署无线传感器网络，该无线传感器网络中存在1个位置和速度未知的用于发射无线信号的源节点、N个位置和速度已知的用于接收无线信号的传感器节点、1个用于估计源节点的位置和速度及源节点和传感器节点之间的时钟偏差和时钟漂移的中心节点，将源节点的位置记为x^o，将源节点的速度记为v^o，将N个传感器节点的位置分别记为s₁,s₂,…,s_i,…,s_N，将N个传感器节点的速度分别记为

其中，N为正整数，若无线传感器网络部署在平面空间中则取N≥5，若无线传感器网络部署在立体空间中则取N≥6，s₁表示第1个传感器节点的位置，s₂表示第2个传感器节点的位置，s_i表示第i个传感器节点的位置，s_N表示第N个传感器节点的位置,

表示第1个传感器节点的速度，

表示第2个传感器节点的速度，

表示第i个传感器节点的速度，

表示第N个传感器节点的速度，i为正整数，1≤i≤N。

步骤2：根据到达时间TOA观测模型，获取每个传感器节点接收到的无线信号的到达时间，将第i个传感器节点接收到的无线信号的到达时间记为τ_i，τ_i描述为：

并根据到达频率FOA观测模型，获取每个传感器节点接收到的无线信号的到达频率，将第i个传感器节点接收到的无线信号的到达频率记为f_i，f_i描述为：

其中，符号“||||”为求欧氏距离符号，c表示无线信号的传播速度，

表示未知的时钟偏差，f₀表示载波频率，上标符号“T”表示向量或矩阵的转置，

表示未知的时钟漂移，Δτ_i为TOA测量噪声，Δf_i为FOA测量噪声。

步骤3：每个传感器节点将其接收到的无线信号的到达时间和到达频率发送给中心节点。

步骤4：中心节点将每个传感器节点接收到的无线信号的到达时间转变为距离观测量ROA(Range of Arrival)，将第i个传感器节点接收到的无线信号的距离观测量记为d_i，d_i描述为：

同样，中心节点将每个传感器节点接收到的无线信号的到达频率转变为距离变化率观测量RROA(Range Rate of Arrival)，将第i个传感器节点接收到的无线信号的距离变化率观测量记为b_i，b_i描述为：

其中，

n_i＝Δτ_ic，

然后令d表示距离观测向量，令b表示距离变化率观测向量，令d^o表示距离真实向量，令b^o表示距离变化率真实向量，d＝[d₁,…,d_i,…,d_N]^T，b＝[b₁,…,b_i,…,b_N]^T，

其中，d₁表示第1个传感器节点接收到的无线信号的距离观测量，d_N表示第N个传感器节点接收到的无线信号的距离观测量，b₁表示第1个传感器节点接收到的无线信号的距离变化率观测量，b_N表示第N个传感器节点接收到的无线信号的距离变化率观测量，

和

根据

计算得到，

和

根据

计算得到；并引入两个噪声向量，分别记为n和m，n＝[n₁,…,n_i,…,n_N]^T，m＝[m₁,…,m_i,…,m_N]^T，其中，n和m相互独立，n服从均值为0的高斯分布，且协方差矩阵Q_n＝E[nn^T]，m服从均值为0的高斯分布，且协方差矩阵Q_m＝E[mm^T]，E[]为求数学期望，n₁和n_N根据n_i＝Δτ_ic计算得到，m₁和m_N根据

计算得到。

如果接下来采用最大似然估计(ML)模型进行求解，那么接下来的过程为：将所有的未知量构成一个向量，用θ^o表示，

其中，θ^o的维度为(2k+2)，若无线传感器网络部署在平面空间中则k＝2，若无线传感器网络部署在立体空间中则k＝3；再将d和b构成一个向量，用p表示，p＝[d^T,b^T]^T；并将d^o和b^o构成一个向量，用p^o表示，p^o＝[d^oT,b^oT]^T；最后构建最大似然估计(ML)模型，描述为：

其中，θ为θ^o对应的优化变量，θ＝[x^T,v^T,d₀,b₀]^T，x为x^o对应的优化变量，v为v^o对应的优化变量，d₀为

对应的优化变量，b₀为

对应的优化变量，

p(θ)为将p^o中的θ^o用θ代替得到，Q表示权重矩阵，Q＝blkdiag(Q_n,Q_m)，blkdiag()表示进行块对角矩阵操作，Q^-1表示Q的逆矩阵，

表示求使得

的值最小时θ的值，

表示求使得(p-p(θ))^TQ^-1(p-p(θ))的值最小时θ的值。

步骤5：由于对最大似然估计(ML)模型求解容易陷入局部最优点，影响定位精度，因此本发明对

进行移项，变换为

然后对

两边进行平方操作，并忽略噪声平方项

得到

并将

代入

中，得到

再展开

后进行移项，并忽略二阶噪声项n_im_i，得到

步骤6：引入辅助向量y^o，y^o定义为

然后将

改写成矩阵形式，描述为：A₁y^o-h₁≈B₁n；并将

改写成矩阵形式，描述为：A₂y^o-h₂≈B₂n+D₂m；接着联立A₁y^o-h₁≈B₁n和A₂y^o-h₂≈B₂n+D₂m得到Ay^o-h≈Bn+Dm；再根据Ay^o-h≈Bn+Dm构建加权最小二乘WLS问题，描述为：

其中，

A₁的维数为N×(2k+6)，0_1×k表示维数为1×k且元素全为0的列向量，若无线传感器网络部署在平面空间中则k＝2，若无线传感器网络部署在立体空间中则k＝3，

h₁的维数为N×1，

表示将

作为对角元素组成对角矩阵，

A₂的维数为N×(2k+6)，

h₂的维数为N×1，

表示将

作为对角元素组成对角矩阵，

表示将

作为对角元素组成对角矩阵，

A的维数为2N×(2k+6)，

h的维数为2N×1，

B的维数为2N×N，

D的维数为2N×N，0_N×N表示维数为N×N且元素全为0的矩阵，

表示求使得(Ay-h)^TR^-1(Ay-h)的值最小时y的值，y为y^o对应的优化变量，R表示权重矩阵，R＝BQ_nB^T+DQ_mD^T＝[BD]Q[B D]^T，R^-1表示R的逆，[B D]为维数为2N×2N的矩阵。

步骤7：由于上述加权最小二乘WLS问题仍旧是一个非凸的优化问题，求解比较困难，因此本发明通过发掘y^o对应的优化变量y中的元素之间隐含的关系，将加权最小二乘WLS问题重新描述为：

其中，“s.t.”表示“受约束于……”，y(1:k)表示y中的第1个元素至第k个元素所组成的向量，y(k+1:2k)表示y中的第(k+1)个元素至第2k个元素所组成的向量，y(2k+1)表示y中的第(2k+1)个元素，y(2k+2)表示y中的第(2k+2)个元素，y(2k+3)表示y中的第(2k+3)个元素，y(2k+4)表示y中的第(2k+4)个元素，y(2k+5)表示y中的第(2k+5)个元素，y(2k+6)表示y中的第(2k+6)个元素。

步骤8：引入矩阵Y＝yy^T，将加权最小二乘WLS问题

松弛为

其中，

表示求使得

的值最小时Y和y的值，tr()表示求矩阵的迹，

表示

是半正定的，rank()表示求矩阵的秩，

和

源于如下等价公式：

步骤9：将半正定松弛后的问题

中的约束条件||y(1:k)||²＝y(2k+3)重写为tr(Y(1:k,1:k))＝y(2k+3)、y^T(1:k)y(k+1:2k)＝y(2k+4)重写为tr(Y(1:k,k+1:2k))＝y(2k+4)、y²(2k+1)＝y(2k+5)重写为Y(2k+1,2k+1)＝y(2k+5)、y(2k+1)y(2k+2)＝y(2k+6)重写为Y(2k+1,2k+2)＝y(2k+6)、y(2k+1)y(2k+6)＝y(2k+2)y(2k+5)重写为Y(2k+1,2k+6)＝y(2k+2,2k+5)；其中，Y(1:k,1:k)表示Y中的第1行至第k行、第1列至第k列的元素所组成的矩阵，Y(1:k,k+1:2k)表示Y中的第1行至第k行、第(k+1)列至第2k列的元素所组成的矩阵，Y(2k+1,2k+1)表示Y中的第2k+1行第2k+1列的元素，Y(2k+1,2k+2)表示Y中的第2k+1行第2k+2列的元素，Y(2k+1,2k+6)表示Y中的第2k+1行第2k+6列的元素。

步骤10：丢弃

中的

并结合tr(Y(1:k,1:k))＝y(2k+3)、tr(Y(1:k,k+1:2k))＝y(2k+4)、Y(2k+1,2k+1)＝y(2k+5)、Y(2k+1,2k+2)＝y(2k+6)、Y(2k+1,2k+6)＝y(2k+2,2k+5)，得到半正定规划问题，描述为：

步骤11：令R的初始值为维数为2N×2N的单位矩阵I_2N；然后根据

计算得到F的值；接着利用内点法求解半正定规划问题

得到Y的初步全局最优解和y的初步全局最优解，对应记为Y^*1和y^*1；再根据y的定义及y^*1，得到

各自的初步估计值，对应记为

其中，y^*1(2k+1)表示y^*1中的第(2k+1)个元素，y^*1(2k+2)表示y^*1中的第(2k+2)个元素。

步骤12：根据

并结合R＝BQ_nB^T+DQ_mD^T＝[B D]Q[B D]^T，更新R的值，进而更新F的值；接着再次利用内点法求解半正定规划问题

得到Y的最终全局最优解和y的最终全局最优解，对应记为Y^*2和y^*2；再根据y的定义及y^*2，得到x^o、v^o、

各自的最终估计值，对应记为x^*2、v^*2、

x^*2＝y^*2(1:k)，v^*2＝y^*2(k+1:2k)，

再根据

和

和

得到

和

的最终估计值，对应记为

和

其中，y^*2(1:k)表示y^*2中的第1个元素至第k个元素所组成的向量，y^*2(k+1:2k)表示y^*2中的第(k+1)个元素至第2k个元素所组成的向量，y^*2(2k+1)表示y^*2中的第(2k+1)个元素，y^*2(2k+2)表示y^*2中的第(2k+2)个元素。

通过仿真实验，可以验证本发明方法的有效性和可行性。

在二维平面空间上部署9个(N＝9)传感器节点，其位置和速度如表1所示，二维笛卡尔坐标为(X，Y)。

表1二维场景下的9个传感器节点的位置和速度

源节点的位置在一个圆环区域内随机产生，该圆环的中心在原点，外半径和内半径分别为1000和200，单位为米；源节点的速度在一个圆内随机产生，该圆的中心在原点，半径为20，单位为米/秒。因此，源节点有可能落在传感器节点组成的凸包(Convex Hull)内部，也有可能落在传感器节点组成的凸包(Convex Hull)外部。时钟偏差引起的

从一个均匀分布U(0,150)中随机产生(单位为米)，时钟漂移造成的

从一个均匀分布U(0,15)中随机产生(单位为米/秒)。

每个传感器节点接收到的无线信号的距离观测量和距离变化率观测量用步骤4产生。测量噪声的协方差矩阵分别为

和

表示TOA测量噪声的方差，

表示FOA测量噪声的方差，I_N表示维数为N×N的单位矩阵。一般来看，FOA测量噪声要小于TOA测量噪声，因此，可以用

来描述两者的关系，且FOA测量噪声系数η的取值范围一般为[0.001,1]，本发明的仿真中，前2个场景(即场景1和场景2)使用的FOA测量噪声系数η取0.1。

目标定位方法的性能可通过均方误差(RMSE)表示，源节点的位置估计性能RMSE定义为：

其中，M₁表示源节点的个数，M₂表示每个源节点做蒙特卡洛仿真的次数，1≤q≤M₁，1≤j≤M₂，x_q表示第q个源节点在第j次蒙特卡洛仿真中的位置真实值，

表示第q个源节点在第j次蒙特卡洛仿真中的位置估计值。源节点的速度估计及时钟偏差和时钟漂移的估计性能RMSE可以用相同的方法来定义，且在此取M₁＝50(50个源节点)，M₂＝500(每个源节点做500次蒙特卡洛仿真)。

为了体现本发明方法的RMSE性能，引入ML方法(最大似然估计方法)进行比较。ML方法是参数估计中常用的方法，在具体实现过程中，该方法需要一个初始值，然后进行迭代计算，在实验中将本发明方法的解作为初始值代入ML方法进行计算。

为了体现处理时钟偏差和时钟漂移的重要性，还引入了Raw WLS方法进行比较。Raw WLS方法把时钟偏差和时钟漂移作为噪声处理，直接忽略其影响，且Raw WLS方法只估计源节点的位置和速度。Raw WLS方法为在本发明方法的基础上修改得到，其具体步骤如下：

在步骤4中的ROA观测向量中，若忽略时钟偏差

则变为d_i＝||x^o-s_i||+n_i,i＝1,…,N，在RROA观测向量中，若忽略时钟漂移

则变为

参照步骤5，对ROA和RROA进行展开并忽略二阶噪声项，分别得到

和

之后，参照步骤6的方法，引入辅助向量

定义为

然后将

改写成矩阵形式，描述为：

并将

改写成矩阵形式，描述为：

接着联立

和

得到

再根据

构建加权最小二乘WLS问题，描述为：

其中

A₃的维数为N×(2k+2)，0_1×k表示维数为1×k且元素全为0的列向量，若无线传感器网络部署在平面空间中则k＝2，若无线传感器网络部署在立体空间中则k＝3，

h₃的维数为N×1，B₃＝diag(-2d₁,…,-2d_N)，B₃＝diag(-2d₁,…,-2d_N)表示将-2d₁,…,-2d_N作为对角元素组成对角矩阵，

A₄的维数为N×(2k+2)，

h₄的维数为N×1，B₄＝diag(-b₁,…,-b_N)，diag(-b₁,…,-b_N)表示将-b₁,…,-b_N作为对角元素组成对角矩阵，D₄＝diag(-d₁,…,-d_N)，diag(-d₁,…,-d_N)表示将-d₁,…,-d_N作为对角元素组成对角矩阵，

A_RWLS的维数为2N×(2k+2)，

h_RWLS的维数为2N×1，

B_RWLS的维数为2N×N，

D_RWLS的维数为2N×N，0_N×N表示维数为N×N且元素全为0的矩阵，

表示求使得

的值最小时y_RWLS的值，y_RWLS为

对应的优化变量，R_RWLS表示权重矩阵，

表示R_RWLS的逆，[B_RWLS D_RWLS]为维数为2N×2N的矩阵。

场景1：传感器节点取表1中的前5个，

的变化范围从10^-2m²变化至10³m²。

图2给出了本发明方法(简称SDP)、ML方法、Raw WLS方法的位置估计的RMSE性能随

变化的曲线，图3给出了本发明方法(简称SDP)、ML方法、Raw WLS方法的速度估计的RMSE性能随

变化的曲线。从图2和图3中可以看出，本发明方法与ML方法的RMSE性能非常接近，在小噪声和中等噪声范围内，本发明方法可以达到CRLB(Cramer-RaoLower Bound)界。由于Raw WLS方法没有考虑时钟偏差和时钟漂移，因此Raw WLS方法的RMSE性能在小噪声和中等噪声区间内受时钟偏差和时钟漂移的控制，远差于本发明方法。在噪声非常大的情况下

时钟偏差和时钟漂移的影响要小于测量噪声的影响，Raw WLS方法的RMSE性能要优于本发明方法和ML方法，原因在于Raw WLS方法估计的参量更少。

图4给出了本发明方法(简称SDP)、ML方法的时钟偏差估计的RMSE性能随

变化的曲线，图5给出了本发明方法(简称SDP)、ML方法的时钟漂移估计的RMSE性能随

变化的曲线。从图4和图5中可以看出，本发明方法在小噪声和中等噪声区间内同样能够达到CRLB界，表现出了较好的RMSE性能。为了评估本发明方法松弛的程度，表3给出了每一个

对应的噪声条件下，M₁＝50、M₂＝500即50个源节点进行的25000次蒙特卡洛仿真实验中本发明方法的Y的最终全局最优解的秩为1的数量(Rank-1数量)。从表2中可以看出，大部分情况下，本发明方法是紧的，即使在

为30m²的时候，Y的最终全局最优解的秩为1的数量(Rank-1数量)的比例都超过99.74％。

表2二维场景下SDP解的Rank-1数量

场景2：

固定为10m²，传感器节点的数量从表1中的前4个增加至9个。

图6给出了本发明方法(简称SDP)、ML方法、Raw WLS方法的位置估计的RMSE性能随传感器节点的数量变化的曲线，图7给出了本发明方法(简称SDP)、ML方法、Raw WLS方法的速度估计的RMSE性能随传感器节点的数量变化的曲线。从图6和图7中可以看出，本发明方法除了在传感器节点的数量为4的这个点以外，其他都能够达到CRLB界。另外，在该测量噪声条件下，无论是位置估计还是速度估计，Raw WLS方法的RMSE性能都比本发明方法和ML方法的RMSE性能要差。这同时体现了处理时钟偏差和时钟漂移的重要性。

场景3：为了比较同时使用TOA和FOA观测信号进行定位与仅使用TOA观测信号进行定位的RMSE性能差别，将本发明方法与文献“Second-order cone relaxation for TOA-based source localization with unknown start transmission time(未知起始传输时间下的基于TOA的二阶锥松弛定位)”(G.Wang,S.Cai,Y.Li,and M.Jin，IEEETrans.Veh.Technol.,vol.63,no.6,pp.2973-2977,2014.)的方法(简称TOA-SOCP)进行对比。

在此测试了二维情况下的RMSE性能，该场景使用了表1中的前5个传感器节点；测量噪声

固定为10m²。

图8给出了FOA测量噪声系数η从10^-3变化至10^-1时，本发明方法、TOA-SOCP方法的位置估计的RMSE性能随FOA测量噪声系数变化的曲线。从图8可以看出，由于FOA测量噪声系数的变化不影响TOA观测信号，因此，TOA-SOCP方法的定位性能无变化。另外，同时使用TOA和FOA观测信号进行定位，与仅使用TOA观测信号进行定位相比，能够显著提升定位性能，且FOA观测信号的噪声越小，对于整体定位性能的提升效果越好。同时，本发明方法与TOA-SOCP方法相比，更加靠近CRLB，验证了本发明方法的性能更好。图8中，TOA-CRLB表示仅使用TOA观测信号时的CRLB界，TOA/FOA-CRLB表示同时使用TOA和FOA观测信号的CRLB界。

Claims (1)

1.一种时钟偏差和时钟漂移条件下的运动目标定位方法，其特征在于包括以下步骤：

步骤1：在一个平面空间或者立体空间中部署无线传感器网络，该无线传感器网络中存在1个位置和速度未知的用于发射无线信号的源节点、N个位置和速度已知的用于接收无线信号的传感器节点、1个用于估计源节点的位置和速度及源节点和传感器节点之间的时钟偏差和时钟漂移的中心节点，将源节点的位置记为x^o，将源节点的速度记为v^o，将N个传感器节点的位置分别记为s₁,s₂,…,s_i,…,s_N，将N个传感器节点的速度分别记为

其中，N为正整数，若无线传感器网络部署在平面空间中则取N≥5，若无线传感器网络部署在立体空间中则取N≥6，s₁表示第1个传感器节点的位置，s₂表示第2个传感器节点的位置，s_i表示第i个传感器节点的位置，s_N表示第N个传感器节点的位置,

表示第1个传感器节点的速度，

表示第2个传感器节点的速度，

表示第i个传感器节点的速度，

表示第N个传感器节点的速度，i为正整数，1≤i≤N；

步骤2：根据到达时间TOA观测模型，获取每个传感器节点接收到的无线信号的到达时间，将第i个传感器节点接收到的无线信号的到达时间记为τ_i，τ_i描述为：

并根据到达频率FOA观测模型，获取每个传感器节点接收到的无线信号的到达频率，将第i个传感器节点接收到的无线信号的到达频率记为f_i，f_i描述为：

其中，符号“|| ||”为求欧氏距离符号，c表示无线信号的传播速度，

表示未知的时钟偏差，f₀表示载波频率，上标符号“T”表示向量或矩阵的转置，

表示未知的时钟漂移，Δτ_i为TOA测量噪声，Δf_i为FOA测量噪声；

步骤3：每个传感器节点将其接收到的无线信号的到达时间和到达频率发送给中心节点；

步骤4：中心节点将每个传感器节点接收到的无线信号的到达时间转变为距离观测量，将第i个传感器节点接收到的无线信号的距离观测量记为d_i，d_i描述为：

同样，中心节点将每个传感器节点接收到的无线信号的到达频率转变为距离变化率观测量，将第i个传感器节点接收到的无线信号的距离变化率观测量记为b_i，b_i描述为：

其中，

n_i＝Δτ_ic，

然后令d表示距离观测向量，令b表示距离变化率观测向量，令d^o表示距离真实向量，令b^o表示距离变化率真实向量，d＝[d₁,…,d_i,…,d_N]^T，b＝[b₁,…,b_i,…,b_N]^T，

其中，d₁表示第1个传感器节点接收到的无线信号的距离观测量，d_N表示第N个传感器节点接收到的无线信号的距离观测量，b₁表示第1个传感器节点接收到的无线信号的距离变化率观测量，b_N表示第N个传感器节点接收到的无线信号的距离变化率观测量，

和

根据

计算得到，

和

根据

计算得到；并引入两个噪声向量，分别记为n和m，n＝[n₁,…,n_i,…,n_N]^T，m＝[m₁,…,m_i,…,m_N]^T，其中，n和m相互独立，n服从均值为0的高斯分布，且协方差矩阵Q_n＝E[nn^T]，m服从均值为0的高斯分布，且协方差矩阵Q_m＝E[mm^T]，E[]为求数学期望，n₁和n_N根据n_i＝Δτ_ic计算得到，m₁和m_N根据

计算得到；

步骤5：对

进行移项，变换为

然后对

两边进行平方操作，并忽略噪声平方项

得到

并将

代入

中，得到

再展开

后进行移项，并忽略二阶噪声项n_im_i，得到

步骤6：引入辅助向量y^o，y^o定义为

然后将

改写成矩阵形式，描述为：A₁y^o-h₁≈B₁n；并将

改写成矩阵形式，描述为：A₂y^o-h₂≈B₂n+D₂m；接着联立A₁y^o-h₁≈B₁n和A₂y^o-h₂≈B₂n+D₂m得到Ay^o-h≈Bn+Dm；再根据Ay^o-h≈Bn+Dm构建加权最小二乘WLS问题，描述为：

其中，

A₁的维数为N×(2k+6)，0_1×k表示维数为1×k且元素全为0的列向量，若无线传感器网络部署在平面空间中则k＝2，若无线传感器网络部署在立体空间中则k＝3，

h₁的维数为N×1，

表示将

作为对角元素组成对角矩阵，

A₂的维数为N×(2k+6)，

h₂的维数为N×1，

表示将

作为对角元素组成对角矩阵，

表示将

作为对角元素组成对角矩阵，

A的维数为2N×(2k+6)，

h的维数为2N×1，

B的维数为2N×N，

D的维数为2N×N，0_N×N表示维数为N×N且元素全为0的矩阵，

表示求使得(Ay-h)^TR^-1(Ay-h)的值最小时y的值，y为y^o对应的优化变量，R表示权重矩阵，R＝BQ_nB^T+DQ_mD^T＝[BD]Q[B D]^T，R^-1表示R的逆，[B D]为维数为2N×2N的矩阵；

步骤7：通过发掘y^o对应的优化变量y中的元素之间隐含的关系，将加权最小二乘WLS问题进一步描述为：

其中，“s.t.”表示“受约束于……”，y(1:k)表示y中的第1个元素至第k个元素所组成的向量，y(k+1:2k)表示y中的第(k+1)个元素至第2k个元素所组成的向量，y(2k+1)表示y中的第(2k+1)个元素，y(2k+2)表示y中的第(2k+2)个元素，y(2k+3)表示y中的第(2k+3)个元素，y(2k+4)表示y中的第(2k+4)个元素，y(2k+5)表示y中的第(2k+5)个元素，y(2k+6)表示y中的第(2k+6)个元素；

步骤8：引入矩阵Y＝yy^T，将加权最小二乘WLS问题

松弛为

其中，

表示求使得

的值最小时Y和y的值，tr()表示求矩阵的迹，

表示

是半正定的，rank()表示求矩阵的秩，

和

源于如下等价公式：

步骤9：将半正定松弛后的问题描述

中的约束条件||y(1:k)||²＝y(2k+3)重写为tr(Y(1:k,1:k))＝y(2k+3)、y^T(1:k)y(k+1:2k)＝y(2k+4)重写为tr(Y(1:k,k+1:2k))＝y(2k+4)、y²(2k+1)＝y(2k+5)重写为Y(2k+1,2k+1)＝y(2k+5)、y(2k+1)y(2k+2)＝y(2k+6)重写为Y(2k+1,2k+2)＝y(2k+6)、y(2k+1)y(2k+6)＝y(2k+2)y(2k+5)重写为Y(2k+1,2k+6)＝y(2k+2,2k+5)；其中，Y(1:k,1:k)表示Y中的第1行至第k行、第1列至第k列的元素所组成的矩阵，Y(1:k,k+1:2k)表示Y中的第1行至第k行、第(k+1)列至第2k列的元素所组成的矩阵，Y(2k+1,2k+1)表示Y中的第2k+1行第2k+1列的元素，Y(2k+1,2k+2)表示Y中的第2k+1行第2k+2列的元素，Y(2k+1,2k+6)表示Y中的第2k+1行第2k+6列的元素；

步骤10：丢弃

中的

并结合tr(Y(1:k,1:k))＝y(2k+3)、tr(Y(1:k,k+1:2k))＝y(2k+4)、Y(2k+1,2k+1)＝y(2k+5)、Y(2k+1,2k+2)＝y(2k+6)、Y(2k+1,2k+6)＝y(2k+2,2k+5)，得到半正定规划问题，描述为：

步骤11：令R的初始值为维数为2N×2N的单位矩阵I_2N；然后根据

计算得到F的值；接着利用内点法求解半正定规划问题

得到Y的初步全局最优解和y的初步全局最优解，对应记为Y^*1和y^*1；再根据y的定义及y^*1，得到

各自的初步估计值，对应记为

其中，y^*1(2k+1)表示y^*1中的第(2k+1)个元素，y^*1(2k+2)表示y^*1中的第(2k+2)个元素；

步骤12：根据

并结合R＝BQ_nB^T+DQ_mD^T＝[B D]Q[B D]^T，更新R的值，进而更新F的值；接着再次利用内点法求解半正定规划问题

得到Y的最终全局最优解和y的最终全局最优解，对应记为Y^*2和y^*2；再根据y的定义及y^*2，得到x^o、v^o、

各自的最终估计值，对应记为x^*2、v^*2、

v^*2＝y^*2(k+1:2k)，

再根据

和

和

得到

和

的最终估计值，对应记为

和

其中，y^*2(1:k)表示y^*2中的第1个元素至第k个元素所组成的向量，y^*2(k+1:2k)表示y^*2中的第(k+1)个元素至第2k个元素所组成的向量，y^*2(2k+1)表示y^*2中的第(2k+1)个元素，y^*2(2k+2)表示y^*2中的第(2k+2)个元素。

Images