CN104198991A - 基于改进Sinc插值的小范围高精度定位方法 - Google Patents

Abstract

本发明提供了一种基于改进Sinc插值的小范围高精度定位方法，包括：L个固定的参考节点接收目标节点发出锯齿调频波v(t)信号；对接收的信号进行幅度均衡并对锯齿调频波以时间间隔T进行采样,得到样值函数x_i(n)；对样值函数进行重构,得到重构函数y(k)；对发送的锯齿调频波v(t)以时间间隔T/(M+1)采样得到样值函数u(k),把y(k)和u(k)进行互相关运算得相关峰值位置；利用相关峰值位置的差分别计算信号到达各参考节点和到达第1个参考节点的时间差，得到TDOA测量值；将TDOA测量值和参考节点的坐标代入Chan氏算法计算目标节点的位置；求目标节点的定位计算结果和真实位置的均方根误差得到定位精度RMSE。本发明通过降低采样频率，保证了定位的高精度。

Description

基于改进Sinc插值的小范围高精度定位方法

技术领域

本发明涉及无线通信领域，尤其涉及一种基于改进Sinc插值的小范围高精度定位方法。

背景技术

全球定位系统(GPS)的使用大幅度的提高了无线定位精度，使定位精度可达到米级，但是在小范围环境内GPS定位效果并不好，米级的定位精度远远不足。在一些精度要求比较高的应用系统中，现有定位系统往往不能满足用户的要求。如机场大厅、超市、图书馆、停车场等复杂的小范围环境，由于受到室内环境等条件的影响，目前的定位技术在精度上往往达不到需求。如何在小范围定位系统中提高定位精度是各国学者研究的热点。在小范围条件下对于复杂环境而言，室内布置，材料结构，建筑物尺度的不同导致了信号的路径损耗很大，与此同时，小范围复杂的环境会引起信号的反射，绕射，折射和散射，形成多径现象，使得接收信号的幅度，相位和到达时间发生变化，造成信号的损失，定位的难度大。相比于GPS系统，小范围的定位系统具有比较鲜明的特点与优势。小范围的定位系统其精度可达到厘米级，而且对于复杂的环境具有较强的适应性，应用范围更加的广阔。而且随着嵌入式技术、通信技术、信息处理技术的发展和成熟，大量适用于复杂环境的小型无线设备的出现也给小范围无线定位的研究和应用带来了方便。因此，小范围定位技术在近几年逐渐成为人们关注的一个研究方向，在将来也会继续成为研究的热点。

小范围无线定位技术通过对接收到的无线电波的一些参数进行测量，根据待定的算法判断出被测物体的位置，其主要机理是根据接收到的定位信号的强度(SOA：Strength of Arrival)、信号到达的角度(AOA：Angle of Arrival)、信号到达的时间(TOA：Time of Arrival)和信号到达的时间差(TDOA：TimeDifference of Arrival)以及这些参数的不同组合来进行移动终端的定位过程。到达时间差的定位技术TDOA是目前无线定位系统中应用最广泛的一项技术。与TOA定位技术相比，它不要求参考节点和目标节点之间的精确同步，与AOA定位技术相比，它不需要特殊的天线阵列，因而更加经济。由于SOA技术受周围环境影响很大，因此定位精度不及TODA高。在TODA技术中，主要有基于最小二乘算法(LS算法)，Chan氏算法，泰勒级数展开法等定位算法。当对定位精度要求不是很高的时候，需要对目标节点位置做初始估计时，采用TDOA最小二乘算法。Chan氏算法是对LS算法的一种改进，先采用加权最小二乘算法(WLS)得到一个初始解，通过利用第一次得到的估计位置坐标及附加变量等已知的约束条件来进行第二次WLS估计能得到改进的估计位置。

TDOA作为一种无线定位技术，是基于时间测量的定位方法。它实现的原理可分为两个基本步骤：1、利用时延估计技术将不同参考节点所接收到的目标节点发出的信号到达时间差进行精确估计，同时利用距离和时间的关系构造非线性双曲线方程；2、将估计出的参数带入方程，通过适当的算法求解方程，得到移动台位置参数。TDOA技术主要是用于构建定位的基础数学模型。TDOA是通过检测信号到达两个基站的时间差，而不是到达的绝对时间来确定移动台的位置，降低了时间同步要求。采用三个不同的基站可以测到两个TDOA，移动站位于两个TDOA决定的双曲线的交点上。

传统的基于最小二乘算法的TDOA定位技术实现流程如附图1所示。其最小二乘算法为：假设待定位的目标节点坐标为(x,y)，(x_i,y_i)为第i个参考节点的已知位置，共L个参考节点，如果要实现对移动台的二维坐标的定位，至少需要三个基站参与TDOA的测量，也即获得两个TDOA测量值，从而构成非线性方程组：

$\left[\begin{array}{c}{r}_{\mathrm{2,1}}\\ r_{\mathrm{3,1}}\\ \·\\ \·\\ \·\\ r_{i,1}\\ \·\\ \·\\ \·\\ r_{L,1}\end{array}\right]=c\left[\begin{array}{c}{t}_{\mathrm{2,1}}\\ t_{\mathrm{3,1}}\\ \·\\ \·\\ \·\\ t_{i,1}\\ \·\\ \·\\ \·\\ t_{L,1}\end{array}\right]=\left[\begin{array}{c}{r}_2-r_1\\ r_3-r_1\\ \·\\ \·\\ \·\\ r_i-r_1\\ \·\\ \·\\ \·\\ r_L-r_1\end{array}\right]=\left[\begin{array}{c}\sqrt{{(x_2-x)}^2+{(y_2-y)}^2}-\sqrt{{(x_1-x)}^2+{(y_1-y)}^2}\\ \sqrt{{(x_3-x)}^2+{(y_3-y)}^2}-\sqrt{{(x_1-x)}^2+{(y_1-y)}^2}\\ \·\\ \·\\ \·\\ \sqrt{{(x_i-x)}^2+{(y_i-y)}^2}-\sqrt{{(x_1-x)}^2+{(y_1-y)}^2}\\ \·\\ \·\\ \·\\ \sqrt{{(x_L-x)}^2+{(y_L-y)}^2}-\sqrt{{(x_1-x)}^2+{(y_1-y)}^2}\end{array}\right]---\left(1\right)$

为了求解该非线性方程组，可以先对其进行线性化处理。

$r_{i,1}^2+{2r}_{i,1}{r}_1=K_i-{2x}_{i,1}x-{2y}_{i,1}y-K_1---\left(2\right)$

式中x_i,1＝x_i-x₁,y_i,1＝y_i-y₁,(其中，i取2,3…L)。x,y,r₁是未知数，则上式为线性方程组。利用最小二乘算法求解方程组。

利用 $r_{i,1}^2+{2r}_{i,1}{r}_1=K_i-{2x}_{i,1}x-{2y}_{i,1}y-K_1$ 建立方程组：

此方程可写作：

AX＝b   (3)

(3)式中：

$A=\left[\begin{array}{cc}{x}_{\mathrm{2,1}}& y_{\mathrm{2,1}}\\ x_{\mathrm{3,1}}& y_{\mathrm{3,1}}\\ \·& \·\\ \·& \·\\ \·& \·\\ x_{L,1}& y_{L,1}\end{array}\right],$ $X=\left[\begin{array}{c}x\\ y\end{array}\right],$ $b=-\{\left[\begin{array}{c}{r}_{\mathrm{2,1}}\\ r_{\mathrm{3,1}}\\ \·\\ \·\\ \·\\ r_{L,1}\end{array}\right]r_1+\frac{1}{2}\left[\begin{array}{c}{r}_{\mathrm{2,1}}^2-K_2+K_1\\ r_{\mathrm{3,1}}^2-K_3+K_1\\ \·\\ \·\\ \·\\ r_{L,1}^2-K_L+K_1\end{array}\right]\};$

现在需要求解X，使得残差r＝AX-b的平方和最小，即

f(X)＝(AX-b)²＝(AX-b)^T(AX-b)   (4)

对(4)式求导并令其为零，得

$\frac{\mathrm{df}\left(X\right)}{\mathrm{dX}}={2A}^T\mathrm{AX}-{2A}^{T}b=0$

若(A^TA)非奇异，则得到方程AX＝b的最小二乘解：

X＝(A^TA)^-1A^Tb

实际应用中通常采用加权最小二乘法，即根据每个测量值的精度，在最小二乘中

采用不同的权值以提高定位精度。残差加权平方和函数为：

f(X)＝(AX-b)^TW(AX-b)

此时得到的加权最小二乘解为:

X＝(A^TWA)^-1A^TWb

式中W为加权矩阵，通常取测量误差的方差矩阵的逆矩阵可使估计误差的方差最小。

上述定位技术在实际应用中直接求解非线性双曲线方程比较困难，，且运用TDOA技术对目标节点定位要想得到10厘米以内的精度需要接收端的采样频率达到10GHz，这在实际应用中很难实现。

发明内容

本发明目的是提供一种基于改进Sinc插值的小范围高精度定位方法，以解决上述问题。

本发明提供了一种基于改进Sinc插值的小范围高精度定位方法，包括如下步骤：

步骤一：L个固定的参考节点接收目标节点发送的锯齿调频波v(t)的信号，其中L为正整数；

步骤二：对接收的信号进行幅度均衡，并对M个连续的锯齿调频波以时间间隔T进行采样，得到样值函数x_i(n)，其中i＝1，2，…，M；n＝0，1，…，N-1，N为每个锯齿调频波的样值点个数，其中M为正整数；

步骤三：利用Sinc插值算法对样值函数进行重构,得到重构函数y(k)，其中k＝0，1，…，(N-1)(M+1)+M。

步骤四：对发送的锯齿调频波v(t)直接以间隔时间T/(M+1)采样，得到样值函数u(k),其中k＝0,1,…,(N-1)(M+1)+M，把重构函数y(k)和u(k)进行互相关运算得到相关峰值位置A_i，其中i取1,2,3…L；

步骤五：利用相关峰值位置的差分别计算信号到达第2,3,4,5,…,L个参考节点和到达第1个参考节点的时间差t_2,1,t_3,1,…,t_i,1,…t_L,1，得到TDOA测量值，其中t_i,1＝(A_i-A₁)*T,T为y(k)的样值点时间间隔；

步骤六：将TDOA测量值和参考节点的坐标(x₁,y₁),(x₂,y₂),…,(x_i,y_i),…,(x_L,y_L)带入Chan氏算法计算目标节点的位置；

步骤七：求目标节点的定位计算结果和真实位置(x,y)的均方根误差得到定位精度RMSE。

进一步，步骤二中时间间隔T的选择满足：1/T>2f_v，其中f_v为v(t)的信号带宽，f_v的计算公式如下：

f_v＝2(m_v+1)F_m

式中m_v为v(t)的调制指数,F_m为的锯齿调制信号的最大频率分量，F_m＝7F，其中F为锯齿波峰出现的频率。

进一步，步骤三具体包括：

A：用Sinc插值算法对步骤二得到的样值函数x_i(n)进行重构，得到重构函数y_i(h)，i＝1，2，…，M；h＝0，1，…，(N-1)(M+1)+M，其中N为正整数；

x_i(n)，(i＝1，2，…，M；n＝0，2，…，N-1)为以T为间隔经过采样得到的已知序列，在两个样值点之间插入M个点，构成长度为M的序列c_i(m,n),(i,m＝1，2，…，M；n＝0，1，…，N-1),则可得：

$c_i(m,n)=\underset{l=0}{\overset{N-1}{\mathrm{\Σ}}}{x}_i\left(l\right)\frac{\sin \left[(n+m/(M+1)-l)\mathrm{\π}\right]}{(n+m/(M+1)-l)\mathrm{\π}}---\left(5\right)$

把样值函数x_i(n)和样值点之间的M个插值点结合，得到重构函数y_i(h)；

y_i(0)＝x_i(0),y_i(M+1)＝x_i(1),…,y_i((N-1)(M+1))＝x_i(N-1),

y_i(1)＝c₁(i,0),y_i(M+2)＝c₁(i,1),…,y_i((N-1)(M+1)+1)＝c₁(i,N-1),

                              .

                              .

                              .

y_i(M)＝c_M(i,M),y_i(2M+1)＝c_M(i,M),…,y_i((N-1)(M+1)+M)＝c_M(i,N-1)；

即：

$y_i\left(h\right)=\left\{\begin{array}{cc}{x}_i\left(n\right),& h=n(M+1)\\ c_i(m,n),& h=n(M+1)+i\end{array}\right.$

i,m＝1,2,…,M n＝0,1,2,…,N-1 h＝0,1,2,…,(N-1)(M+1)+M

B：用y₁(h),y₂(h),y₃(h),…,y_M(h)再次重构得到重构函数y(k)为

y(0)＝y₁(0),y(M+1)＝y₁(M+1),…,y((N-1)(M+1))＝y₁((N-1)(M+1)),

y(1)＝y₁(1),y(M+2)＝y₁(M+2),…,y((N-1)(M+1)+1)＝y₁((N-1)(M+1)+1),

y(2)＝y₂(2),y(M+3)＝y₂(M+3),…,y((N-1)(M+1)+2)＝y₂((N-1)(M+1)+2),

                             .

                             .

                             .

y(i)＝y_i(i),y(M+1+i)＝y_i(M+1+i),…,y((N-1)(M+1)+i)＝y_i((N-1)(M+1)+i),

                                .

                                .

                                .

y(M)＝y_M(M),y(2M+1)＝y_M(2M+1),…,y((N-1)(M+1)+M)＝y_M((N-1)(M+1)+M)；

即：

$y\left(k\right)=\left\{\begin{array}{cc}{y}_1\left(n(M+1)\right),& k=n(M+1)\\ y_i(n(M+1)+i),& k=n(M+1)+i\end{array}\right.$

其中n＝0,1,2…N-1,i＝1,2,3,…,M。

进一步，步骤六具体包括：

A:TDOA构造构成非线性方程组：

假设待定位的目标节点坐标为(x,y)，(x_i,y_i)为第i个参考节点的已知位置，共L个参考节点，获得两个TDOA测量值，构成非线性方程组：

$\left[\begin{array}{c}{r}_{\mathrm{2,1}}\\ r_{\mathrm{3,1}}\\ \·\\ \·\\ \·\\ r_{i,1}\\ \·\\ \·\\ \·\\ r_{L,1}\end{array}\right]=c\left[\begin{array}{c}{t}_{\mathrm{2,1}}\\ t_{\mathrm{3,1}}\\ \·\\ \·\\ \·\\ t_{i,1}\\ \·\\ \·\\ \·\\ t_{L,1}\end{array}\right]=\left[\begin{array}{c}{r}_2-r_1\\ r_3-r_1\\ \·\\ \·\\ \·\\ r_i-r_1\\ \·\\ \·\\ \·\\ r_L-r_1\end{array}\right]=\left[\begin{array}{c}\sqrt{{(x_2-x)}^2+{(y_2-y)}^2}-\sqrt{{(x_1-x)}^2+{(y_1-y)}^2}\\ \sqrt{{(x_3-x)}^2+{(y_3-y)}^2}-\sqrt{{(x_1-x)}^2+{(y_1-y)}^2}\\ \·\\ \·\\ \·\\ \sqrt{{(x_i-x)}^2+{(y_i-y)}^2}-\sqrt{{(x_1-x)}^2+{(y_1-y)}^2}\\ \·\\ \·\\ \·\\ \sqrt{{(x_L-x)}^2+{(y_L-y)}^2}-\sqrt{{(x_1-x)}^2+{(y_1-y)}^2}\end{array}\right]---\left(6\right)$

进行线性化处理：

$r_{i,1}^2+{2r}_{i,1}{r}_1=K_i-{2x}_{i,1}x-{2y}_{i,1}y-K_1---\left(7\right)$

(6)式中x_i,1＝x_i-x₁,y_i,1＝y_i-y₁,且i取2,3…L；x,y,r₁是未知数，(7)式为线性方程组；

利用Chan氏算法求解目标节点的坐标；

B：Chan氏算法求解第一次加权平均：

令为未知矢量，其中z_p＝[x,y]^T，假定z_a的元素间相互独立，则z_a的ML估计为： $z_a={\left(G_a^T{\mathrm{\Ψ}}^{-1}{G}_a\right)}^{-1}{G}_a^T{\mathrm{\Ψ}}^{-1}h;$

其中 $G_a=-\left[\begin{array}{ccc}{x}_{\mathrm{2,1}}& y_{\mathrm{2,1}}& r_{\mathrm{2,1}}\\ x_{\mathrm{3,1}}& y_{\mathrm{3,1}}& r_{\mathrm{3,1}}\\ \·& \·& \·\\ \·& \·& \·\\ \·& \·& \·\\ x_{L,1}& y_{L,1}& r_{L,1}\end{array}\right],$ $h=\frac{1}{2}\left[\begin{array}{c}{r}_{\mathrm{2,1}}^2-x_2^2-y_2^2+x_1^2+y_1^2\\ r_{\mathrm{3,1}}^2-x_3^2-y_2^2+x_1^2+y_1^2\\ \·\\ \·\\ \·\\ r_{L,1}^2-x_L^2-y_L^2+x_1^2+y_1^2\end{array}\right],$ Ψ＝c²BQB， $B=\left[\begin{array}{ccccc}{r}_2^0& 0& 0& ...& 0\\ 0& r_3^0& 0& ...& 0\\ \·& \·& \·& & \\ \·& \·& \·& \·\·\·& 0\\ \·& \·& \·& & \\ 0& 0& 0& ...& r_L^0\end{array}\right],$ $Q=\left[\begin{array}{cccc}1& 0& ...& 0\\ 0& 1& ...& 0\\ \·& \·& & \·\\ \·& \·& \·\·\·& \·\\ \·& \·& & \·\\ 0& 0& ...& 1\end{array}\right];$

进一步的近似式为：

$z_a\≈{\left(G_a^T{Q}^{-1}{G}_a\right)}^{-1}{G}_a^T{Q}^{-1}h$

利用近似式得到初始解，利用所得到的初始解计算B矩阵，再利用 $z_a={\left(G_a^T{\mathrm{\Ψ}}^{-1}{G}_a\right)}^{-1}{G}_a^T{\mathrm{\Ψ}}^{-1}h$ 求得第一次的WLS计算结果；

C：Chan氏算法求解第二次加权平均求出定位结果：

进行第二次WLS计算，首先计算估计位置z_a的协方差矩阵；其协方差矩阵为：矢量z_a是均值为实际值、协方差矩阵由 $\mathrm{cov}\left(z_a\right)={\left(G_a^T{\mathrm{\Ψ}}^{-1}{G}_a\right)}^{-1}$ 确定的随机矢量；

求得方程组：ψ′＝h′-G′_az′_a，

其中 $h^{\′}=\left[\begin{array}{c}{(z_{a,1}-x_1)}^2\\ {(z_{a,2}-y_1)}^2\\ z_{a,3}^2\end{array}\right],$ $G_a^{\′}=\left[\begin{array}{cc}1& 0\\ 0& 1\\ 1& 1\end{array}\right],$ $z_a^{\′}=\left[\begin{array}{c}{(x-x_1)}^2\\ {(y-y_1)}^2\end{array}\right];$

定义ψ′为z_a的误差矢量；ψ′的协方差矩阵为Ψ′＝[ψ′ψ′^T]＝4B′cov(z_a)B′，其中， $B^{\′}=\left[\begin{array}{ccc}{x}^0-x_1& 0& 0\\ 0& y^0-y_1& 0\\ 0& 0& r_1^0\end{array}\right],$ B′中的x⁰,y⁰,即为第一次的WLS计算结果；

Ψ′为高斯分布，z′_a的最大似然估计为：求得目标节点的定位计算结果为： $z_p=\sqrt{z_a^{\′}}+\left[\begin{array}{c}{x}_1\\ y_1\end{array}\right];$

求得目标节点和第i参考节点与第1个参考节点的距离差r_i,1，其公式为:r_i,1＝t_i,1*c，其中c为电磁波的传播速度；将r_i,1、L和(x_i,y_i)代入到Chan氏算法中，求解非线性方程组式(6)，获得目标节点的坐标(x,y)，其中(x_i,y_i)为第i参考节点的坐标。

本发明提供了一种基于改进Sinc插值的小范围高精度定位方法，通过改进Sinc插值算法，利用多个样值函数进行样值估算，得到多个样值点后插入其中一个样值函数的两个样值点之间，实现了对样值函数的重构，其与现有技术相比，通过降低采样频率，保证了定位的高精度。

附图说明

图1为现有技术中基于最小二乘算法的TDOA定位技术的流程图；

图2为本发明基于改进Sinc插值的小范围高精度定位方法Sinc插值法的示意图；

图3为本发明基于改进Sinc插值的小范围高精度定位方法的流程图；

图4为本发明v(n)，y(k)两个函数TDOA估计广义互相关方法的流程图：

图5为本发明v(n)，y(k)两个函数广义互相关之后得到相关峰的示意图；

图6为本发明的算法与现有技术中的算法下定位精度比较的示意图。

具体实施方式

本实施例的思路是：参考节点的接收端在满足抽样定理的条件下抽样，能够获得信号的全部信息，用这些数据和插值算法对信号进行重构，以达到降低参考节点接收端采样频率的目的。由于Sinc函数在频域中对应有限带宽的矩形低通滤波器，只要其频带比信号的频带宽那么用Sinc函数理论上可以无损的还原信号。普通Sinc插值算法直接通过一个样值函数估算并在两个样值点之间插入多个估算值实现对样值函数的重构。为了提高定位精度，本发明提出了改进了Sinc插值算法，该算法利用多个样值函数进行样值估算，得到多个样值点后插入其中一个样值函数的两个样值点之间实现对样值函数的重构。

下面通过具体的实施例子并结合附图对本发明做进一步的详细描述。

参图2及图3所示，图2为本发明基于改进Sinc插值的小范围高精度定位方法Sinc插值法的示意图；图3为本发明基于改进Sinc插值的小范围高精度定位方法的流程图。

在本实施例中，一种基于改进Sinc插值的小范围高精度定位方法，包括：

步骤S101，L个固定的参考节点接收目标节点发出的锯齿调频波v(n)信号，经过信道到达接收端；其中，L为正整数。

步骤S102，对接收的信号进行幅度均衡并对M个连续的锯齿调频波以时间间隔T进行采样，得到样值函数x_i(n)，其中i＝1，2，…，M；n＝0，1，…，N-1，N为每个锯齿调频波的样值点个数，其中，M为正整数。

在本实施例中，时间间隔T可按下述方法选择：

1/T必须大于2倍v(t)的信号带宽f_v，即1/T>2f_v；v(t)的信号带宽用中等质量通信系统方法估算，即f_v＝2(m_v+1)F_m，式中f_v为v(t)的信号带宽，m_v为v(t)的调制指数,F_m为的锯齿调制信号的最大频率分量，F_m＝7F，其中F为锯齿波峰出现的频率；

步骤S103，利用改进Sinc插值算法对样值函数进行重构,得到重构函数y(k)，其中k＝0，1，…，(N-1)(M+1)+M；

步骤S103的具体包括如下子步骤：

A：用Sinc插值算法对步骤二得到的样值函数x_i(n)进行重构，得到重构函数y_i(h)，i＝1，2，…，M；h＝0，1，…，(N-1)(M+1)+M；

x_i(n)，(i＝1，2，…，M；n＝0，2，…，N-1)为以T为间隔经过采样得到的已知序列，在两个样值点之间插入M个点，构成长度为M的序列c_i(m,n),(i,m＝1，2，…，M；n＝0，1，…，N-1),则可得：

$c_i(m,n)=\underset{l=0}{\overset{N-1}{\mathrm{\Σ}}}{x}_i\left(l\right)\frac{\sin \left[(n+m/(M+1)-l)\mathrm{\π}\right]}{(n+m/(M+1)-l)\mathrm{\π}}---\left(5\right)$

把样值函数x_i(n)和样值点之间的M个插值点结合起来，得到了重构函数y_i(h)；

y_i(0)＝x_i(0),y_i(M+1)＝x_i(1),…,y_i((N-1)(M+1))＝x_i(N-1),

y_i(1)＝c₁(i,0),y_i(M+2)＝c₁(i,1),…,y_i((N-1)(M+1)+1)＝c₁(i,N-1),

                                .

                                .

                                .

y_i(M)＝c_M(i,M),y_i(2M+1)＝c_M(i,M),…,y_i((N-1)(M+1)+M)＝c_M(i,N-1)；

即：

$y_i\left(h\right)=\left\{\begin{array}{cc}{x}_i\left(n\right),& h=n(M+1)\\ c_i(m,n),& h=n(M+1)+i\end{array}\right.$

i,m＝1,2,…,M n＝0,1,2,…,N-1 h＝0,1,2,…,(N-1)(M+1)+M

B：用y₁(h),y₂(h),y₃(h),…,y_M(h)再次重构得到重构函数y(k)为

y(0)＝y₁(0),y(M+1)＝y₁(M+1),…,y((N-1)(M+1))＝y₁((N-1)(M+1)),

y(1)＝y₁(1),y(M+2)＝y₁(M+2),…,y((N-1)(M+1)+1)＝y₁((N-1)(M+1)+1),

y(2)＝y₂(2),y(M+3)＝y₂(M+3),…,y((N-1)(M+1)+2)＝y₂((N-1)(M+1)+2),

                              .

                              .

                              .

y(i)＝y_i(i),y(M+1+i)＝y_i(M+1+i),…,y((N-1)(M+1)+i)＝y_i((N-1)(M+1)+i),

                                .

                                .

                                .

y(M)＝y_M(M),y(2M+1)＝y_M(2M+1),…,y((N-1)(M+1)+M)＝y_M((N-1)(M+1)+M)；

即：

$y\left(k\right)=\left\{\begin{array}{cc}{y}_1\left(n(M+1)\right),& k=n(M+1)\\ y_i(n(M+1)+i),& k=n(M+1)+i\end{array}\right.$ (其中n＝0,1,2…N-1,i＝1,2,3,…,M)。

步骤S104，对发送的锯齿调频波v(t)直接以间隔时间T/(M+1)采样得到样值函数u(k),其中k＝0,1,…,(N-1)(M+1)+M，把重构函数y(k)和u(k)进行互相关运算得到相关峰值位置A_i(其中i取1,2,3…L)；

参图4及图5所示，图4为本发明v(n)，y(k)两个函数TDOA估计广义互相关方法的流程图；图5为本发明v(n)，y(k)两个函数广义互相关之后得到相关峰的示意图。

图5中A_i为每个相关峰值对应的样值点所在的位置，图中只示出了4个峰值，实际峰值个数等于参考节点个数L。

步骤S105，利用相关峰值位置的差分别计算信号到达第2,3,4,5,…,L个参考节点和到达第1个参考节点的时间差t_2,1,t_3,1,…,t_i,1,…t_L,1，其中t_i,1＝(A_i-A₁)*T,T为y(k)的样值点时间间隔。

步骤S106，将TDOA测量值，即t_2,1,t_3,1,…,t_i,1,…t_L,1和参考节点的坐标(x₁,y₁),(x₂,y₂),…,(x_i,y_i),…,(x_L,y_L)带入Chan氏算法计算目标节点的位置。

在本实施例中，步骤S106具体包括如下子步骤：

A:TDOA构造构成非线性方程组：

假设待定位的目标节点坐标为(x,y)，(x_i,y_i)为第i个参考节点的已知位置，共L个参考节点，如果要实现对移动台的二维坐标的定位，至少需要三个基站参与TDOA的测量，也即获得两个TDOA测量值，从而构成非线性方程组：

$\left[\begin{array}{c}{r}_{\mathrm{2,1}}\\ r_{\mathrm{3,1}}\\ \·\\ \·\\ \·\\ r_{i,1}\\ \·\\ \·\\ \·\\ r_{L,1}\end{array}\right]=c\left[\begin{array}{c}{t}_{\mathrm{2,1}}\\ t_{\mathrm{3,1}}\\ \·\\ \·\\ \·\\ t_{i,1}\\ \·\\ \·\\ \·\\ t_{L,1}\end{array}\right]=\left[\begin{array}{c}{r}_2-r_1\\ r_3-r_1\\ \·\\ \·\\ \·\\ r_i-r_1\\ \·\\ \·\\ \·\\ r_L-r_1\end{array}\right]=\left[\begin{array}{c}\sqrt{{(x_2-x)}^2+{(y_2-y)}^2}-\sqrt{{(x_1-x)}^2+{(y_1-y)}^2}\\ \sqrt{{(x_3-x)}^2+{(y_3-y)}^2}-\sqrt{{(x_1-x)}^2+{(y_1-y)}^2}\\ \·\\ \·\\ \·\\ \sqrt{{(x_i-x)}^2+{(y_i-y)}^2}-\sqrt{{(x_1-x)}^2+{(y_1-y)}^2}\\ \·\\ \·\\ \·\\ \sqrt{{(x_L-x)}^2+{(y_L-y)}^2}-\sqrt{{(x_1-x)}^2+{(y_1-y)}^2}\end{array}\right]---\left(6\right)$

进行线性化处理：

$r_{i,1}^2+{2r}_{i,1}{r}_1=K_i-{2x}_{i,1}x-{2y}_{i,1}y-K_1---\left(7\right)$

(6)式中x_i,1＝x_i-x₁,y_i,1＝y_i-y₁,且i取2,3…L；x,y,r₁是未知数，(7)式为线性方程组，利用Chan氏算法就可以求解目标节点的坐标；

B：Chan氏算法求解第一次加权平均：

令为未知矢量，其中z_p＝[x,y]^T，如果假定z_a的元素间相互独立，则z_a的ML估计为： $z_a={\left(G_a^T{\mathrm{\Ψ}}^{-1}{G}_a\right)}^{-1}{G}_a^T{\mathrm{\Ψ}}^{-1}h;$

其中 $G_a=-\left[\begin{array}{ccc}{x}_{\mathrm{2,1}}& y_{\mathrm{2,1}}& r_{\mathrm{2,1}}\\ x_{\mathrm{3,1}}& y_{\mathrm{3,1}}& r_{\mathrm{3,1}}\\ \·& \·& \·\\ \·& \·& \·\\ \·& \·& \·\\ x_{L,1}& y_{L,1}& r_{L,1}\end{array}\right],$ $h=\frac{1}{2}\left[\begin{array}{c}{r}_{\mathrm{2,1}}^2-x_2^2-y_2^2+x_1^2+y_1^2\\ r_{\mathrm{3,1}}^2-x_3^2-y_2^2+x_1^2+y_1^2\\ \·\\ \·\\ \·\\ r_{L,1}^2-x_L^2-y_L^2+x_1^2+y_1^2\end{array}\right],$ Ψ＝c²BQB， $B=\left[\begin{array}{ccccc}{r}_2^0& 0& 0& ...& 0\\ 0& r_3^0& 0& ...& 0\\ \·& \·& \·& & \\ \·& \·& \·& \·\·\·& 0\\ \·& \·& \·& & \\ 0& 0& 0& ...& r_L^0\end{array}\right],$ $Q=\left[\begin{array}{cccc}1& 0& ...& 0\\ 0& 1& ...& 0\\ \·& \·& & \·\\ \·& \·& \·\·\·& \·\\ \·& \·& & \·\\ 0& 0& ...& 1\end{array}\right];$

进一步的近似为：

$z_a\≈{\left(G_a^T{Q}^{-1}{G}_a\right)}^{-1}{G}_a^T{Q}^{-1}h$

利用该近似式可以得到一个初始解，利用这个初始解计算B矩阵，再利用 $z_a={\left(G_a^T{\mathrm{\Ψ}}^{-1}{G}_a\right)}^{-1}{G}_a^T{\mathrm{\Ψ}}^{-1}h$ 求得第一次的WLS计算结果；

C：Chan氏算法求解第二次加权平均求出定位结果：

进行第二次WLS计算，首先计算估计位置z_a的协方差矩阵；其协方差矩阵为：矢量z_a是一个均值为实际值、协方差矩阵由确定的随机矢量，因此我们可以得到另一个方程组：ψ′＝h′-G′_az′_a，

其中 $h^{\′}=\left[\begin{array}{c}{(z_{a,1}-x_1)}^2\\ {(z_{a,2}-y_1)}^2\\ z_{a,3}^2\end{array}\right],$ $G_a^{\′}=\left[\begin{array}{cc}1& 0\\ 0& 1\\ 1& 1\end{array}\right],$ $z_a^{\′}=\left[\begin{array}{c}{(x-x_1)}^2\\ {(y-y_1)}^2\end{array}\right];$

定义ψ′为z_a的误差矢量；ψ′的协方差矩阵为Ψ′＝[ψ′ψ′^T]＝4B′cov(z_a)B′其中 $B^{\′}=\left[\begin{array}{ccc}{x}^0-x_1& 0& 0\\ 0& y^0-y_1& 0\\ 0& 0& r_1^0\end{array}\right],$ B′中的x⁰,y⁰,即为第一次的WLS计算结果；

由于Ψ为高斯分布，故Ψ′也为高斯分布，z′_a的最大似然估计为： $z_a^{\′}={\left(G_a^T{\mathrm{\Ψ}}^{\′-1}{G}_a^{\′}\right)}^{-1}{G}_a^T{\mathrm{\Ψ}}^{\′-1}{h}^{\′}$ 则目标节点的定位计算结果为： $z_p=\sqrt{z_a^{\′}}+\left[\begin{array}{c}{x}_1\\ y_1\end{array}\right];$

由建立的非线性方程组(6)可得看出方程组的参数有参考节点数目L，目标节点和第i参考节点与第1个参考节点的距离差r_i,1，第i参考节点的坐标(x_i,y_i)以及要求的未知参数目标节点的坐标(x,y)；通过TODA值可求出r_i,1；r_i,1＝t_i,1*c，c为电磁波的传播速度；L和(x_i,y_i)为已知参数；把三个参数代入到Chan氏算法中，求解非线性方程组式(6)，即可求得目标节点的坐标(x,y)；

步骤S107，求目标节点的定位计算结果和真实位置(x,y)的均方根误差得到定位精度RMSE。

参图6所示，为本发明的算法与现有技术中的算法下定位精度比较的示意图。由图6可以看出在没有插值算法的情况下，当采用频率从10GHz下降到1GHz时，定位精度从小于0.01m下降到3.80m左右，定位精度下降明显。当在采用频率1GHz采用普通Sinc算法在每两个样值点之间插入9个点之后，定位精度和在采用频率为1GHz没用插值算法相比没用明显提升，定位精度只是从3.80m左右下降到3.4m左右。当在采用频率1GHz采用改进Sinc算法在每两个样值点之间插入9个点之后，定位精度会有明显的提升达到0.01m左右，其定位精度接近采样频率为10GHz的定位精度。

本发明通过提供了一种基于改进Sinc插值的小范围高精度定位方法，采用改进Sinc插值算法，利用多个样值函数进行样值估算，得到多个样值点后插入其中一个样值函数的两个样值点之间，实现了对样值函数的重构，降低采样频率，保证了定位的高精度。

Claims (4)

1.一种基于改进Sinc插值的小范围高精度定位方法，其特征在于，包括：

步骤一：L个固定的参考节点接收目标节点发送的锯齿调频波v(t)的信号，其中L为正整数；

步骤二：对接收的信号进行幅度均衡，并对M个连续的锯齿调频波以时间间隔T进行采样，得到样值函数x_i(n)，其中i＝1，2，…，M；n＝0，1，…，N-1，N为每个锯齿调频波的样值点个数，其中M为正整数；

步骤三：利用Sinc插值算法对样值函数进行重构,得到重构函数y(k)，其中k＝0，1，…，(N-1)(M+1)+M。

步骤四：对发送的锯齿调频波v(t)直接以间隔时间T/(M+1)采样，得到样值函数u(k),其中k＝0,1,…,(N-1)(M+1)+M，把重构函数y(k)和u(k)进行互相关运算得到相关峰值位置A_i，其中i取1,2,3…L；

步骤五：利用相关峰值位置的差分别计算信号到达第2,3,4,5,…,L个参考节点和到达第1个参考节点的时间差t_2,1,t_3,1,…,t_i,1,…t_L,1，得到TDOA测量值，其中t_i,1＝(A_i-A₁)*T,T为y(k)的样值点时间间隔；

步骤六：将TDOA测量值，和参考节点的坐标(x₁,y₁),(x₂,y₂),…,(x_i,y_i),…,(x_L,y_L)带入Chan氏算法计算目标节点的位置；

步骤七：求目标节点的定位计算结果和真实位置(x,y)的均方根误差得到定位精度RMSE。

2.如权利要求1所述的基于改进Sinc插值的小范围高精度定位方法，其特征在于，所述时间间隔T的选择满足：1/T>2f_v，其中f_v为v(t)的信号带宽，f_v的计算公式如下：

f_v＝2(m_v+1)F_m

式中m_v为v(t)的调制指数,F_m为的锯齿调制信号的最大频率分量，F_m＝7F，其中F为锯齿波峰出现的频率。

3.如权利要求2所述的基于改进Sinc插值的小范围高精度定位方法，其特征在于，步骤三具体包括：

A：用Sinc插值算法对步骤二得到的样值函数x_i(n)进行重构，得到重构函数y_i(h)，i＝1，2，…，M；h＝0，1，…，(N-1)(M+1)+M，其中N为正整数；

x_i(n)，(i＝1，2，…，M；n＝0，2，…，N-1)为以T为间隔经过采样得到的已知序列，在两个样值点之间插入M个点，构成长度为M的序列c_i(m,n),(i,m＝1，2，…，M；n＝0，1，…，N-1),则可得：

$c_i(m,n)=\underset{l=0}{\overset{N-1}{\mathrm{\Σ}}}{x}_i\left(l\right)\frac{\sin \left[(n+m/(M+1)-l)\mathrm{\π}\right]}{(n+m/(M+1)-l)\mathrm{\π}}---\left(5\right)$

把样值函数x_i(n)和样值点之间的M个插值点结合，得到重构函数y_i(h)；

y_i(0)＝x_i(0),y_i(M+1)＝x_i(1),…,y_i((N-1)(M+1))＝x_i(N-1),

y_i(1)＝c₁(i,0),y_i(M+2)＝c₁(i,1),…,y_i((N-1)(M+1)+1)＝c₁(i,N-1),

                                .

                                .

                                .

y_i(M)＝c_M(i,M),y_i(2M+1)＝c_M(i,M),…,y_i((N-1)(M+1)+M)＝c_M(i,N-1)；

即：

$y_i\left(h\right)=\left\{\begin{array}{cc}{x}_i\left(n\right),& h=n(M+1)\\ c_i(m,n),& h=n(M+1)+i\end{array}\right.$

i,m＝1,2,…,M n＝0,1,2,…,N-1 h＝0,1,2,…,(N-1)(M+1)+M

B：用y₁(h),y₂(h),y₃(h),…,y_M(h)再次重构得到重构函数y(k)为

y(0)＝y₁(0),y(M+1)＝y₁(M+1),…,y((N-1)(M+1))＝y₁((N-1)(M+1)),

y(1)＝y₁(1),y(M+2)＝y₁(M+2),…,y((N-1)(M+1)+1)＝y₁((N-1)(M+1)+1),

y(2)＝y₂(2),y(M+3)＝y₂(M+3),…,y((N-1)(M+1)+2)＝y₂((N-1)(M+1)+2),

                             .

                             .

                             .

y(i)＝y_i(i),y(M+1+i)＝y_i(M+1+i),…,y((N-1)(M+1)+i)＝y_i((N-1)(M+1)+i),

                             .

                             .

                             .

y(M)＝y_M(M),y(2M+1)＝y_M(2M+1),…,y((N-1)(M+1)+M)＝y_M((N-1)(M+1)+M)；

即：

$y\left(k\right)=\left\{\begin{array}{cc}{y}_1\left(n(M+1)\right),& k=n(M+1)\\ y_i(n(M+1)+i),& k=n(M+1)+i\end{array}\right.$

其中n＝0,1,2…N-1,i＝1,2,3,…,M。

4.如权利要求3所述的基于改进Sinc插值的小范围高精度定位方法，其特征在于，步骤六具体包括：

A:TDOA构造构成非线性方程组：

假设待定位的目标节点坐标为(x,y)，(x_i,y_i)为第i个参考节点的已知位置，共L个参考节点，获得两个TDOA测量值，构成非线性方程组：

$\left[\begin{array}{c}{r}_{\mathrm{2,1}}\\ r_{\mathrm{3,1}}\\ \·\\ \·\\ \·\\ r_{i,1}\\ \·\\ \·\\ \·\\ r_{L,1}\end{array}\right]=c\left[\begin{array}{c}{t}_{\mathrm{2,1}}\\ t_{\mathrm{3,1}}\\ \·\\ \·\\ \·\\ t_{i,1}\\ \·\\ \·\\ \·\\ t_{L,1}\end{array}\right]=\left[\begin{array}{c}{r}_2-r_1\\ r_3-r_1\\ \·\\ \·\\ \·\\ r_i-r_1\\ \·\\ \·\\ \·\\ r_L-r_1\end{array}\right]=\left[\begin{array}{c}\sqrt{{(x_2-x)}^2+{(y_2-y)}^2}-\sqrt{{(x_1-x)}^2+{(y_1-y)}^2}\\ \sqrt{{(x_3-x)}^2+{(y_3-y)}^2}-\sqrt{{(x_1-x)}^2+{(y_1-y)}^2}\\ \·\\ \·\\ \·\\ \sqrt{{(x_i-x)}^2+{(y_i-y)}^2}-\sqrt{{(x_1-x)}^2+{(y_1-y)}^2}\\ \·\\ \·\\ \·\\ \sqrt{{(x_L-x)}^2+{(y_L-y)}^2}-\sqrt{{(x_1-x)}^2+{(y_1-y)}^2}\end{array}\right]---\left(6\right)$

进行线性化处理：

$r_{i,1}^2+{2r}_{i,1}{r}_1=K_i-{2x}_{i,1}x-{2y}_{i,1}y-K_1---\left(7\right)$

(6)式中x_i,1＝x_i-x₁,y_i,1＝y_i-y₁,且i取2,3…L；x,y,r₁是未知数，(7)式为线性方程组；

利用Chan氏算法求解目标节点的坐标；

B：Chan氏算法求解第一次加权平均：

令为未知矢量，其中z_p＝[x,y]^T，假定z_a的元素间相互独立，则z_a的ML估计为： $z_a={\left(G_a^T{\mathrm{\Ψ}}^{-1}{G}_a\right)}^{-1}{G}_a^T{\mathrm{\Ψ}}^{-1}h;$

其中 $G_a=-\left[\begin{array}{ccc}{x}_{\mathrm{2,1}}& y_{\mathrm{2,1}}& r_{\mathrm{2,1}}\\ x_{\mathrm{3,1}}& y_{\mathrm{3,1}}& r_{\mathrm{3,1}}\\ \·& \·& \·\\ \·& \·& \·\\ \·& \·& \·\\ x_{L,1}& y_{L,1}& r_{L,1}\end{array}\right],$ $h=\frac{1}{2}\left[\begin{array}{c}{r}_{\mathrm{2,1}}^2-x_2^2-y_2^2+x_1^2+y_1^2\\ r_{\mathrm{3,1}}^2-x_3^2-y_2^2+x_1^2+y_1^2\\ \·\\ \·\\ \·\\ r_{L,1}^2-x_L^2-y_L^2+x_1^2+y_1^2\end{array}\right],$ Ψ＝c²BQB， $B=\left[\begin{array}{ccccc}{r}_2^0& 0& 0& ...& 0\\ 0& r_3^0& 0& ...& 0\\ \·& \·& \·& & \\ \·& \·& \·& \·\·\·& 0\\ \·& \·& \·& & \\ 0& 0& 0& ...& r_L^0\end{array}\right],$ $Q=\left[\begin{array}{cccc}1& 0& ...& 0\\ 0& 1& ...& 0\\ \·& \·& & \·\\ \·& \·& \·\·\·& \·\\ \·& \·& & \·\\ 0& 0& ...& 1\end{array}\right];$

进一步的近似式为：

$z_a\≈{\left(G_a^T{Q}^{-1}{G}_a\right)}^{-1}{G}_a^T{Q}^{-1}h$

利用所述近似式得到初始解，利用所得到的初始解计算B矩阵，再利用 $z_a={\left(G_a^T{\mathrm{\Ψ}}^{-1}{G}_a\right)}^{-1}{G}_a^T{\mathrm{\Ψ}}^{-1}h$ 求得第一次的WLS计算结果；

C：Chan氏算法求解第二次加权平均求出定位结果：

进行第二次WLS计算，首先计算估计位置z_a的协方差矩阵；其协方差矩阵为：矢量z_a是均值为实际值、协方差矩阵由 $\mathrm{cov}\left(z_a\right)={\left(G_a^T{\mathrm{\Ψ}}^{-1}{G}_a\right)}^{-1}$ 确定的随机矢量；

求得方程组：ψ′＝h′-G′_az′_a，

其中 $h^{\′}=\left[\begin{array}{c}{(z_{a,1}-x_1)}^2\\ {(z_{a,2}-y_1)}^2\\ z_{a,3}^2\end{array}\right],$ $G_a^{\′}=\left[\begin{array}{cc}1& 0\\ 0& 1\\ 1& 1\end{array}\right],$ $z_a^{\′}=\left[\begin{array}{c}{(x-x_1)}^2\\ {(y-y_1)}^2\end{array}\right];$

定义ψ′为z_a的误差矢量；ψ′的协方差矩阵为Ψ′＝[ψ′ψ′^T]＝4B′cov(z_a)B′，

其中， $B^{\′}=\left[\begin{array}{ccc}{x}^0-x_1& 0& 0\\ 0& y^0-y_1& 0\\ 0& 0& r_1^0\end{array}\right],$ B′中的x⁰,y⁰,即为第一次的WLS计算结果；

Ψ′为高斯分布，z′_a的最大似然估计为：求得目标节点的定位计算结果为： $z_p=\sqrt{z_a^{\′}}+\left[\begin{array}{c}{x}_1\\ y_1\end{array}\right];$

求得目标节点和第i参考节点与第1个参考节点的距离差r_i,1，其公式为:r_i,1＝t_i,1*c，其中c为电磁波的传播速度；将r_i,1、L和(x_i,y_i)代入到Chan氏算法中，求解非线性方程组式(6)，获得目标节点的坐标(x,y)，其中(x_i,y_i)为第i参考节点的坐标。