CN106162865A - 终端的定位方法及装置 - Google Patents

Abstract

本发明提供了一种终端的定位方法及装置，其中，该方法包括：获取终端的物理无线访问接入点AP的位置信息和虚拟AP的位置信息，其中，虚拟AP的位置信息为物理AP的镜像位置信息；根据物理AP的位置信息和虚拟AP的位置信息按照预设算法对终端进行定位。通过本发明解决了在非视距环境下，相关技术中运用无线定位算法对终端进行定位会导致严重误差的问题，进而达到了提高终端进定位精度的效果。

Description

终端的定位方法及装置

技术领域

本发明涉及通信领域，具体而言，涉及终端的定位方法及装置。

背景技术

目前无线定位技术引起了人们越来越多的研究，其应用也越来越广泛，如雷达、声纳、通信、传感器网络。传统无线定位算法主要是通过测量直射径的信息来实现对目标的定位。但是在没有直射径或者是非直射径占主导地位的场景中，这种测量直射径的定位方法会导致严重的非视距(Non Line Of Sight，简称为NLOS)误差，尤其在复杂的室内环境中这一影响会更加突出。但是非直射径中包含着丰富的定位信息，所以很多定位算法利用多径定位来减少NLOS误差。

基于多径的定位的算法大致可以分为两类：几何定位法和统计定位法。统计定位法通常是将非直射径的测量建模成直射径的测量加上统计误差的方式来减小NLOS的误差，但是这类方法易受模型不匹配的影响。几何定位法通过确定非直射径与移动台位置的额外几何关系来提高定位精度。

常见的几何定位法如三角定位法利用已知的多个参考节点的坐标位置来对未知节点进行定位。参考节点能够发射射频信号，未知节点根据接收到的参考节点的射频信号来获得自身节点到参考节点的距离，未知节点只需测得与至少三个参考节点的距离就可以在二维坐标上实现对未知节点的定位。三角定位法的基本原理如图12所示。

A₁,A₂,A₃为已知坐标的参考节点，其坐标分别为(x₁,y₁),(x₂,y₂),(x₃,y₃)。B为未知节点，坐标设为(x,y)。B为三个圆的交点，A₁,A₂,A₃到节点B的距离为d₁,d₂,d₃，则可以得到

$\left\{\begin{array}{c}{(x-x_1)}^2+{\left({y-y}_1\right)}^2={d_1}^2\\ {\left({x-x}_2\right)}^2+{\left({y-y}_2\right)}^2={d_2}^2\\ {\left({x-x}_3\right)}^2+{\left({y-y}_3\right)}^2={d_3}^2\end{array}\right.$

求解上式可得：

$\left[\begin{array}{c}x\\ y\end{array}\right]=\frac{1}{2}\left[\begin{array}{cc}(x_1-x_3)& (y_1-y_3)\\ (x_2-x_3)& (y_2-y_3)\end{array}\right]\left[\begin{array}{c}{x_1}^2-{x_3}^2+{y_1}^2-{y_3}^2+{d_3}^2-{d_1}^2\\ {x_2}^2-{x_3}^2+{y_2}^2-{y_3}^2+{d_3}^2-{d_2}^2\end{array}\right]$

由于复杂的室内环境，这些方法需要严格的假设才能获得非直射径的测量信息。比如一些文献提出的方法是通过严格假设反射参数，使用直射径和非直射径确定移动终端的位置。而在没有直射径场景中，通过非直射径的反射参数准静态假设，采用扩展卡尔曼滤波器卡尔曼滤波(Extended Kalman Filter，简称为EKF)和概率数据关联滤波(Probabilistic Data AssociationFilter，简称为PDA)方法来跟踪移动终端，但是由于受非直射径参数的时变特性影响，从而导致累计定位误差的增加。

无线定位算法主要是通过测量直射径的信息来确定。但是在没有直射径或者是非直射径占主导地位的场景中，这种定位方法会导致严重的NLOS误差，尤其在复杂的室内环境中这一影响会更加突出。

针对相关技术中，运用无线定位算法对终端进行定位会导致严重的NLOS误差的问题，还未提出有效的解决方案。

发明内容

本发明提供了一种终端的定位方法及装置，以解决相关技术中运用无线定位算法对终端进行定位会导致严重的NLOS误差的问题。

根据本发明的一个方面，提供了一种终端的定位方法，包括：获取终端的物理无线访问接入点AP的位置信息和虚拟AP的位置信息，其中，所述虚拟AP的位置信息为所述物理AP的镜像位置信息；根据所述物理AP的位置信息和所述虚拟AP的位置信息按照预设算法对所述终端进行定位。

进一步地，获取所述终端的虚拟AP的位置信息包括：通过最大期望运算和迭代算法根据所述终端的估计位置信息、所述终端与所述物理AP和/或所述虚拟AP的距离信息获取所述虚拟AP的位置信息。

进一步地，获取所述虚拟AP的位置信息包括：通过以下公式确定所述虚拟AP的位置信息的最大化似然估计值，根据该最大化似然估计值确定所述虚拟AP的位置信息： ${\hat{A}}_{N+1:M}=\underset{A_{N+1:M}}{\arg \max }\left\{\underset{U(1:K)}{\mathrm{\Σ}}p(Z(1:K),U(1:K)|A_{N+1:M};A_{1:N})\right\},$ 其中，上标^表示估计运算，A_1:N是N个所述AP的位置坐标，A_N+1:M是M-N个所述虚拟AP的位置坐标，U(1:K)为所述终端的移动路径上从第1时刻到第K时刻各时间点的位置信息，Z(1:K)为所述终端的移动路径上从第1时刻到第K时刻各时间点到所述物理AP或者所述虚拟AP的距离信息，p(Z(1:K),U(1:K)|A_N+1:M；A_1:N)表示在给定A_1:N和A_N+1:M的条件下，发生Z(1:K)和U(1:K)的概率，N、M为自然数。

进一步地，根据所述最大化似然估计值确定所述虚拟AP的位置信息包括：对所述最大化似然估计值按照以下公式进行期望运算和最大化运算，得到所述虚拟AP的位置信息：所述期望运算通过以下公式实现：所述U(1:K)的概率分布q^t(U(1:K))为： $q^t\left(\hat{U}(1:K)\right)=p\left(\hat{U}(1:K)\right|Z(1:K),{\hat{A}}_{N+1:M}^{t-1};A_{1:N});$ 其中，上标t和t-1分别表示第t次和第t-1次迭代步骤；对于所述U(1:K)的最大似然函数的均值为： $Q(A_{N+1:M},q^t\left(\hat{U}(1:K)\right))=E_{q^t\left(\hat{U}(1:K)\right)}p(Z(1:K),{\hat{U}}^t(1:K)|A_{N+1:M};A_{1:N});$ 通过所述迭代运算，如下公式所示，得到使得最大化虚拟AP位置的估计值： $\begin{array}{c}{\hat{A}}_{N+1:M}=\underset{A_{N+1:M}}{\arg \max }\left\{Q(A_{N+1:M},q^t\left(\hat{U}(1:K)\right))\right\}\\ =\underset{A_{N+1:M}}{\arg \max }\left\{\underset{m=N+1}{\overset{M}{\mathrm{\Σ}}}\log E_{q^t\left(\hat{U}(1:K)\right)}p(Z{(1:K)}_m,{\hat{U}}^t(1:K)|A_m)\right\}\end{array},$ 所述虚拟AP位置的估计值表示为 ${\hat{A}}_m=\underset{A_m}{\arg \max }\left\{E_{q^t\left(\hat{U}(1:K)\right)}p(Z{(1:K)}_m,{\hat{U}}^t(1:K)|A_m)\right\},m=N+1:M.$

进一步地，所述通过轨迹粒子U^(i,t)(1:K)确定：U^(i,t)(1:K)＝[U^(j,t)(1:K-1),U^(i,t)(K)],i＝1,2.....,I，其中，为所述终端作为第i个粒子在第k个时刻、第t次迭代过程中的位置，U^(i,t)(1:K)为从第1时刻到第K时刻所有的集合，I为自然数，为在粒子滤波算法中生成的粒子总数。

进一步地，通过如下公式获取所述虚拟AP的位置信息： ${\hat{A}}_m=\underset{A_m}{\arg \max }\left\{\underset{i=1}{\overset{I}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}p(Z_m\left(k\right),{\hat{U}}^{(i,t)}\left(k\right)|A_m)\right\};$ 其中，Z_m(k)为所述终端在第k个时刻到第m个所述物理AP或者所述虚拟AP的直射径测量距离，Z(1:K)为所述终端的移动路径上从第1时刻到第K时刻各时间点到所述物理AP或者所述虚拟AP的距离信息的集合。

根据本发明的另一个方面，还提供了一种终端的定位装置，包括：获取模块，用于获取终端的物理无线访问接入点AP的位置信息和虚拟AP的位置信息，其中，所述虚拟AP的位置信息为所述物理AP的镜像位置信息；定位模块，用于根据所述物理AP的位置信息和所述虚拟AP的位置信息按照预设算法对所述终端进行定位。

进一步地，所述获取模块还用于通过最大期望运算和迭代算法根据所述终端的估计位置信息、所述终端与所述物理AP和/或所述虚拟AP的距离信息获取所述虚拟AP的位置信息。

进一步地，所述获取模块还用于通过以下公式确定所述虚拟AP的位置信息的最大化似然估计值，根据该最大化似然估计值确定所述虚拟AP的位置信息： ${\hat{A}}_{N+1:M}=\underset{A_{N+1:M}}{\arg \max }\left\{\underset{U(1:K)}{\mathrm{\Σ}}p(Z(1:K),U(1:K)|A_{N+1:M};A_{1:N})\right\},$ 其中，上标^表示估计运算，A_1:N是N个所述AP的位置坐标，A_N+1:M是M-N个所述虚拟AP的位置坐标，U(1:K)为所述终端的移动路径上从第1时刻到第K时刻各时间点的位置信息，Z(1:K)为所述终端的移动路径上从第1时刻到第K时刻各时间点到所述物理AP或者所述虚拟AP的距离信息，p(Z(1:K),U(1:K)|A_N+1:M；A_1:N)表示在给定A_1:N和A_N+1:M的条件下，发生Z(1:K)和U(1:K)的概率，N、M为自然数。

进一步地，所述获取模块还用于对所述最大化似然估计值按照以下公式进行期望运算和最大化运算，得到所述虚拟AP的位置信息：所述期望运算通过以下公式实现：所述U(1:K)的概率分布q^t(U(1:K))为： $q^t\left(\hat{U}(1:K)\right)=p\left(\hat{U}(1:K)\right|Z(1:K),{\hat{A}}_{N+1:M}^{t-1};A_{1:N});$ 其中，上标t和t-1分别表示第t次和第t-1次迭代步骤；对于所述U(1:K)的最大似然函数的均值为：

$Q(A_{N+1:M},q^t\left(\hat{U}(1:K)\right))=E_{q^t\left(\hat{U}(1:K)\right)}p(Z(1:K),{\hat{U}}^t(1:K)|A_{N+1:M};A_{1:N});$ 通过所述迭代运算，如下公式所示，得到使得最大化虚拟AP位置的估计值：

$\begin{array}{c}{\hat{A}}_{N+1:M}=\underset{A_{N+1:M}}{\arg \max }\left\{Q(A_{N+1:M},q^t\left(\hat{U}(1:K)\right))\right\}\\ =\underset{A_{N+1:M}}{\arg \max }\left\{\underset{m=N+1}{\overset{M}{\mathrm{\Σ}}}\log E_{q^t\left(\hat{U}(1:K)\right)}p(Z{(1:K)}_m,{\hat{U}}^t(1:K)|A_m)\right\}\end{array},$ 所述虚拟AP位置的估计值表示为 ${\hat{A}}_m=\underset{A_m}{\arg \max }\left\{E_{q^t\left(\hat{U}(1:K)\right)}p(Z{(1:K)}_m,{\hat{U}}^t(1:K)|A_m)\right\},m=N+1:M.$

进一步地，所述通过轨迹粒子U^(i,t)(1:K)确定：U^(i,t)(1:K)＝[U^(j,t)(1:K-1),U^(i,t)(K)],i＝1,2.....,I，其中，为所述终端作为第i个粒子在第k个时刻、第t次迭代过程中的位置，U^(i,t)(1:K)为从第1时刻到第K时刻所有的集合，I为自然数，为在粒子滤波算法中生成的粒子总数。

进一步地，通过如下公式获取所述虚拟AP的位置信息： ${\hat{A}}_m=\underset{A_m}{\arg \max }\left\{\underset{i=1}{\overset{I}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}p(Z_m\left(k\right),{\hat{U}}^{(i,t)}\left(k\right)|A_m)\right\};$ 其中，Z_m(k)为所述终端在第k个时刻到第m个所述物理AP或者所述虚拟AP的直射径测量距离，Z(1:K)为所述终端的移动路径上从第1时刻到第K时刻各时间点到所述物理AP或者所述虚拟AP的距离信息的集合。

通过本发明，采用获取终端的物理无线访问接入点AP的位置信息和虚拟AP的位置信息，其中，虚拟AP的位置信息为物理AP的镜像位置信息；根据物理AP的位置信息和虚拟AP的位置信息按照预设算法对终端进行定位。解决了相关技术中运用无线定位算法对终端进行定位会导致严重的NLOS误差的问题，进而达到了对终端进行准确定位的效果。

附图说明

此处所说明的附图用来提供对本发明的进一步理解，构成本申请的一部分，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：

图1是根据本发明实施例的终端的定位方法的流程图；

图2是根据本发明实施例的终端的定位装置的结构框图；

图3是根据本发明实施例的AP的坐标和虚拟AP的坐标示意图；

图4是根据本发明实施例的真实AP和虚拟AP坐标以及多径示意图；

图5是根据本发明实施例的MSPE的仿真结果示意图；

图6是根据本发明实施例的状态预测粒子生成示意图；

图7是根据本发明实施例的状态预测粒子权重计算示意图；

图8是根据本发明实施例的状态粒子重采样示意图；

图9是根据本发明实施例的轨迹粒子生成示意图；

图10是根据本发明实施例的EM算法第t次迭代中粒子滤波示意图；

图11是根据本发明实施例的基于EM和粒子滤波的虚拟AP位置估计方案流程图；

图12是三角定位法示意图。

具体实施方式

下文中将参考附图并结合实施例来详细说明本发明。需要说明的是，在不冲突的情况下，本申请中的实施例及实施例中的特征可以相互组合。

在本实施例中提供了一种终端的定位方法，图1是根据本发明实施例的终端的定位方法的流程图，如图1所示，该流程包括如下步骤：

步骤S102，获取终端的物理无线访问接入点AP的位置信息和虚拟AP的位置信息，其中，虚拟AP的位置信息为物理AP的镜像位置信息；

步骤S104，根据物理AP的位置信息和虚拟AP的位置信息按照预设算法对终端进行定位。

通过上述步骤，运用终端的AP的位置信息和虚拟AP的位置信息对终端进行定位，相比于现有技术中，在没有直射径或者是非直射径占主导地位的场景中，无线定位算法会导致严重的NLOS误差，上述步骤解决了相关技术中运用无线定位算法对终端进行定位会导致严重的NLOS误差的问题，进而达到了对终端进行准确定位的效果。

上述步骤S102中涉及到获取虚拟AP的位置信息，在一个可选实施例中，通过最大期望运算和迭代算法根据终端的估计位置信息、终端与物理AP和/或虚拟AP的距离信息获取虚拟AP的位置信息。

在一个可选实施例中，通过以下公式确定虚拟AP的位置信息的最大化似然估计值，根据最大化似然估计值确定该虚拟AP的位置信息： ${\hat{A}}_{N+1:M}=\underset{A_{N+1:M}}{\arg \max }\left\{\underset{U(1:K)}{\mathrm{\Σ}}p(Z(1:K),U(1:K)|A_{N+1:M};A_{1:N})\right\},$ 其中，上标^表示估计运算，A_1:N是N个该AP的位置坐标，A_N+1:M是M-N个该虚拟AP的位置坐标，U(1:K)为终端的移动路径上从第1时刻到第K时刻各时间点的位置信息，Z(1:K)为终端的移动路径上从第1时刻到第K时刻各时间点到该物理AP或者该虚拟AP的距离信息，p(Z(1:K),U(1:K)|A_N+1:M；A_1:N)表示在给定A_1:N和A_N+1:M的条件下，发生Z(1:K)和U(1:K)的概率，N、M为自然数，进而可以获取虚拟AP的位置信息。

在一个可选实施例中，对该最大化似然估计值按照以下公式进行期望运算和最大化运算，得到该虚拟AP的位置信息：该期望运算通过以下公式实现：该U(1:K)的概率分布q^t(U(1:K))为： $q^t\left(\hat{U}(1:K)\right)=p\left(\hat{U}(1:K)\right|Z(1:K),{\hat{A}}_{N+1:M}^{t-1};A_{1:N});$ 其中，上标t和t-1分别表示第t次和第t-1次迭代步骤；对于该U(1:K)的最大似然函数的均值为： $Q(A_{N+1:M},q^t\left(\hat{U}(1:K)\right))=E_{q^t\left(\hat{U}(1:K)\right)}p(Z(1:K),{\hat{U}}^t(1:K)|A_{N+1:M};A_{1:N});$ 通过该迭代运算，如下公式所示，得到使得最大化虚拟AP位置的估计值： $\begin{array}{c}{\hat{A}}_{N+1:M}=\underset{A_{N+1:M}}{\arg \max }\left\{Q(A_{N+1:M},q^t\left(\hat{U}(1:K)\right))\right\}\\ =\underset{A_{N+1:M}}{\arg \max }\left\{\underset{m=N+1}{\overset{M}{\mathrm{\Σ}}}\log E_{q^t\left(\hat{U}(1:K)\right)}p(Z{(1:K)}_m,{\hat{U}}^t(1:K)|A_m)\right\}\end{array},$ 该虚拟AP位置的估计值表示为 ${\hat{A}}_m=\underset{A_m}{\arg \max }\left\{E_{q^t\left(\hat{U}(1:K)\right)}p(Z{(1:K)}_m,{\hat{U}}^t(1:K)|A_m)\right\},m=N+1:M,$ 从而可以根据最大化似然估计值确定虚拟AP的位置信息。

在一个可选实施例中，通过轨迹粒子U^(i,t)(1:K)确定：U^(i,t)(1:K)＝[U^(j,t)(1:K-1),U^(i,t)(K)],i＝1,2.....,I，其中，为该终端作为第i个粒子在第k个时刻、第t次迭代过程中的位置，U^(i,t)(1:K)为从第1时刻到第K时刻所有的集合，I为自然数，为在粒子滤波算法中生成的粒子总数。

在一个可选实施例中，通过如下公式获取该虚拟AP的位置信息： ${\hat{A}}_m=\underset{A_m}{\arg \max }\left\{\underset{i=1}{\overset{I}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}p(Z_m\left(k\right),{\hat{U}}^{(i,t)}\left(k\right)|A_m)\right\};$ 其中，Z_m(k)为该终端在第k个时刻到第m个该物理AP或者该虚拟AP的直射径测量距离，Z(1:K)为该终端的移动路径上从第1时刻到第K时刻各时间点到该物理AP或者该虚拟AP的距离信息的集合。

在本实施例中还提供了一种终端的定位装置，该装置用于实现上述实施例及优选实施方式，已经进行过说明的不再赘述。如以下所使用的，术语“模块”可以实现预定功能的软件和/或硬件的组合。尽管以下实施例所描述的装置较佳地以软件来实现，但是硬件，或者软件和硬件的组合的实现也是可能并被构想的。

图2是根据本发明实施例的终端的定位装置的结构框图，如图2所示，该装置包括：获取模块22，用于获取终端的物理无线访问接入点AP的位置信息和虚拟AP的位置信息，其中，该虚拟AP的位置信息为该物理AP的镜像位置信息；定位模块24，用于根据该物理AP的位置信息和该虚拟AP的位置信息按照预设算法对该终端进行定位。

进一步地，获取模块22还用于通过最大期望运算和迭代算法根据该终端的估计位置信息、该终端与该物理AP和/或该虚拟AP的距离信息获取该虚拟AP的位置信息。

进一步地，获取模块22还用于通过以下公式确定该虚拟AP的位置信息的最大化似然估计值，根据该最大化似然估计值确定该虚拟AP的位置信息： ${\hat{A}}_{N+1:M}=\underset{A_{N+1:M}}{\arg \max }\left\{\underset{U(1:K)}{\mathrm{\Σ}}p(Z(1:K),U(1:K)|A_{N+1:M};A_{1:N})\right\},$ 其中，上标^表示估计运算，A_1:N是N个该AP的位置坐标，A_N+1:M是M-N个该虚拟AP的位置坐标，U(1:K)为该终端的移动路径上从第1时刻到第K时刻各时间点的位置信息，Z(1:K)为该终端的移动路径上从第1时刻到第K时刻各时间点到该物理AP或者该虚拟AP的距离信息，p(Z(1:K),U(1:K)|A_N+1:M；A_1:N)表示在给定A_1:N和A_N+1:M的条件下，发生Z(1:K)和U(1:K)的概率，N、M为自然数。

进一步地，获取模块22还用于对该最大化似然估计值按照以下公式进行期望运算和最大化运算，得到该虚拟AP的位置信息：该期望运算通过以下公式实现：该U(1:K)的概率分布q^t(U(1:K))为： $q^t\left(\hat{U}(1:K)\right)=p\left(\hat{U}(1:K)\right|Z(1:K),{\hat{A}}_{N+1:M}^{t-1};A_{1:N});$ 其中，上标t和t-1分别表示第t次和第t-1次迭代步骤；对于该U(1:K)的最大似然函数的均值为：

$Q(A_{N+1:M},q^t\left(\hat{U}(1:K)\right))=E_{q^t\left(\hat{U}(1:K)\right)}p(Z(1:K),{\hat{U}}^t(1:K)|A_{N+1:M};A_{1:N});$ 通过该迭代运算，如下公式所示，得到使得最大化虚拟AP位置的估计值：

$\begin{array}{c}{\hat{A}}_{N+1:M}=\underset{A_{N+1:M}}{\arg \max }\left\{Q(A_{N+1:M},q^t\left(\hat{U}(1:K)\right))\right\}\\ =\underset{A_{N+1:M}}{\arg \max }\left\{\underset{m=N+1}{\overset{M}{\mathrm{\Σ}}}\log E_{q^t\left(\hat{U}(1:K)\right)}p(Z{(1:K)}_m,{\hat{U}}^t(1:K)|A_m)\right\}\end{array},$ 该虚拟AP位置的估计值表示为 ${\hat{A}}_m=\underset{A_m}{\arg \max }\left\{E_{q^t\left(\hat{U}(1:K)\right)}p(Z{(1:K)}_m,{\hat{U}}^t(1:K)|A_m)\right\},m=N+1:M.$

进一步地，该通过轨迹粒子U^(i,t)(1:K)确定：U^(i,t)(1:K)＝[U^(j,t)(1:K-1),U^(i,t)(K)],i＝1,2.....,I，其中，为该终端作为第i个粒子在第k个时刻、第t次迭代过程中的位置，U^(i,t)(1:K)为从第1时刻到第K时刻所有的集合，I为自然数，为在粒子滤波算法中生成的粒子总数。

进一步地，通过如下公式获取该虚拟AP的位置信息： ${\hat{A}}_m=\underset{A_m}{\arg \max }\left\{\underset{i=1}{\overset{I}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}p(Z_m\left(k\right),{\hat{U}}^{(i,t)}\left(k\right)|A_m)\right\};$ 其中，Z_m(k)为该终端在第k个时刻到第m个该物理AP或者该虚拟AP的直射径测量距离，Z(1:K)为该终端的移动路径上从第1时刻到第K时刻各时间点到该物理AP或者该虚拟AP的距离信息的集合。

需要说明的是，上述各个模块是可以通过软件或硬件来实现的，对于后者，可以通过以下方式实现，但不限于此：上述各个模块均位于同一处理器中；或者，上述各个模块分别位于第一处理器和第二处理器中。

针对相关技术中存在的上述问题，下面结合可选实施例进行说明，在本可选实施例中结合了上述可选实施例及其可选实施方式。

本可选实施例针对室内定位系统，设计了一种可有效地利用多径信息来解决NLOS误差的定位算法。在本可选实施例中，为了描述简单，将基站称为AP。

由于基站和移动终端都可以视为“无线信号发射和接收装置”，因此本可选实施例的方法也可以用于在若干个已知位置的“无线信号发射和接收装置”与一个未知位置的“无线信号发射和接收装置”构成的网络中，通过已知位置的“无线信号发射和接收装置”与未知位置的“无线信号发射和接收装置”之间的无线信号多径信息来提高未知位置“无线信号发射和接收装置”定位的精度。

定义A_1:N是N个真实AP的坐标，A_N+1:M是M-N个虚拟AP的坐标(虚拟AP是指根据反射AP发射的无线信号的镜面反射原理得到的真实AP的镜像点，如图3所示)。

定义状态向量 $X\left(k\right)={\left[\begin{array}{cc}U{\left(k\right)}^T& \tilde{U}{\left(k\right)}^T\end{array}\right]}^T,$ 其中U(k)表示移动台在第k时刻2维位置坐标，上标～表示求导操作，即表示移动台2维速度向量，观测向量Z(k)＝[Z₁(k),…Z_m(k),…Z_M(k)]^T，其中Z_m(k)表示移动台到第m个(虚拟或真实)AP的直射径测量距离，该距离可以通过测量移动台与AP见多个传播路径的传播时延乘以光速得到。可以得到如下状态转移和测量方程：

系统状态转移方程为

X(k)＝F·X(k-1)+V(k)， (1)

V(k)表示均值为0方差为C_V的高斯噪声， $C_V=\left[\begin{array}{cc}0& 0\\ 0& {{\mathrm{\Δt\σ}}_V}^{2}I\end{array}\right],$ I表示2×2维的单位阵，σ_v ²为移动终端速度的方差，Δt为采样间隔，状态矩阵：

$F=\left[\begin{array}{cccc}1& 0& \mathrm{\Δt}& 0\\ 0& 1& 0& \mathrm{\Δt}\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]$

系统的观测方程为：

Z(k)＝g(X(k))+W(k) (2)

其中g(X(k))＝[||U(k)-A₁||,...,||U(k)-A_M||]^T，||·||代表欧式距离。W(k)是均值为0方差为 $\mathrm{\σ}=\mathrm{diag}[{\mathrm{\σ}}_1^2,..{\mathrm{\σ}}_m^2...{\mathrm{\σ}}_M^2]$ 的高斯噪声。

本可选实施例提出的基于最大期望算法(Expectation Maximization Algorithm，简称为EM)和粒子滤波的射频多径定位的主要思路是：定义基于到达时间(Time of Arrival，TOA)时间序列来测量移动台移动路径上不同点距离AP的Z(1:K)，用来确定移动终端轨迹。移动台的定位精度的提升可视为如下式表达的最大化似然估计。下式右侧的概率越大，表示对待定位设备定位的精度越高。

${\hat{A}}_{N+1:M},\hat{U}(1:K)=\underset{A_{N+1:M},U(1:K)}{\arg \max }p\left(Z(1:K)\right|A_{N+1:M},U(1:K);A_{1:N})---\left(3\right)$

其中，上标^表示估计操作，概率函数p(a|b,c；e)表示在参数b,c,e条件下a发生的概率，由于参数b,c和参数e分别为待估计参数和已知参数，中间采用分号分隔。p(Z(1:K)|A_N+1:M,U(1:K)；A_1:N)具体表示在已知的N个AP位置和在假定的M-N个AP的位置A_N+1:M和用户位置U(1:K)的条件下，移动台在移动路径上距离Z(1:K)发生的概率。

首先利用在已知多径测量距离的条件下虚拟AP位置与移动台位置的依赖性，即已知多径测量距离Z(1:K)的条件下，给定虚拟AP坐标A_N+1:M，可以得到用户位置U(1:K)，反过来，如果给定用户位置U(1:K)，可以得到虚拟AP坐标A_N+1:M，此时可以将U(1:K)看作是A_N+1:M的潜变量，进而可以采用潜变量模型将上述最大似然估计等效为

$\begin{array}{c}{\hat{A}}_{N+1:M}=\underset{A_{N+1:M}}{\arg \max }\left\{p\left(Z(1:K)\right|A_{N+1:M};A_{1:N})\right\}\\ =\underset{A_{N+1:M}}{\arg \max }\left\{\underset{U(1:K)}{\mathrm{\Σ}}p(Z(1:K),U(1:K)|A_{N+1:M};A_{1:N})\right\}\end{array}.---\left(4\right)$

这样上述最大化问题可以通过基于潜变量的最大期望算法(Expectation MaximizationAlgorithm，简称为EM)来解决，可以分为两个具体的步骤：

第一步通过EM算法来估计虚拟AP。本步骤重点在于虚拟AP位置确定。考虑到EM算法为迭代算法，有若干相似的迭代步骤组成，下面以其中第t次迭代步骤为例，进行介绍。

(a)估计虚拟AP位置的第t步迭代过程：

期望步骤(E-Step)：假设EM算法迭代到第t-1步时，虚拟AP的估计坐标为那么本次迭代步骤中(即第t次迭代步骤中)潜变量U(1:K)的概率分布q^t(U(1:K))为：

$q^t\left(\hat{U}(1:K)\right)=p\left(\hat{U}(1:K)\right|Z(1:K),{\hat{A}}_{N+1:M}^{t-1};A_{1:N}),---\left(5\right)$

其中上标t和t-1分别表示第t次和第t-1次迭代步骤。

那么对于U(1:K)的最大似然函数的均值为：

$Q(A_{N+1:M},q^t\left(\hat{U}(1:K)\right))=E_{q^t\left(\hat{U}(1:K)\right)}p(Z(1:K),{\hat{U}}^t(1:K)|A_{N+1:M};A_{1:N})---\left(6\right)$

最大化步骤(M-Step)：在第t步，虚拟AP的坐标可以通过最大化似然函数的期望值获得，也就有

$\begin{array}{c}{\hat{A}}_{N+1:M}=\underset{A_{N+1:M}}{\arg \max }\left\{Q(A_{N+1:M},q^t\left(\hat{U}(1:K)\right))\right\}\\ =\underset{A_{N+1:M}}{\arg \max }\left\{\underset{m=N+1}{\overset{M}{\mathrm{\Σ}}}\log E_{q^t\left(\hat{U}(1:K)\right)}p(Z{(1:K)}_m,{\hat{U}}^t(1:K)|A_m)\right\}\end{array}---\left(7\right)$

因为不同AP之间的测量误差是独立的，所以(7)中第二个式子可以等价成M-N个独立的最大化问题：

${\hat{A}}_m=\underset{A_m}{\arg \max }\left\{E_{q^t\left(\hat{U}(1:K)\right)}p(Z{(1:K)}_m,{\hat{U}}^t(1:K)|A_m)\right\},m=N+1:M---\left(8\right)$

由于得到一般是很困难的，这样导致了上述EM算法不可行性。但是蒙特卡洛算法通常可以作为一个替代的方法来描述后验概率密度函数[9]，那么基于粒子滤波的序贯重要性采样(Sequential Importance Sampling，SIS)可以用来估计EM算法中的下面采用SIS粒子滤波方法来近似移动台轨迹的概率分布

(b)利用粒子滤波方法近似第t次EM迭代步骤中移动台轨迹概率分布假设移动终端的初始状态 $X\left(0\right)={\left[\begin{array}{cc}U{\left(0\right)}^T& \tilde{U}{\left(0\right)}^T\end{array}\right]}^T.$ 考虑到粒子滤波方法为迭代算法，以粒子滤波迭代步骤中的第k步(对应第k时刻)迭代过程为例介绍。

首先，第k时刻按照重要性密度函数θ(X^(i,t)(k)|X^(i,t)(k-1),Z(k))生成粒子为其中上标(i，t)表示在第t次EM迭代步骤(以区别不同的迭代步骤)中生成的第i个粒子，I表示生成的粒子总数，在SIS粒子滤波方法中该重要性密度函数按照下式选取：

θ(X^(i,t)(k)|X^(i,t)(k-1),Z(k))＝p(X^(i,t)(k)|X^(i,t)(k-1)) (9)

其中条件概率p(X^(i,t)(k)|X^(i,t)(k-1))表示第k-1时刻状态X^(i,t)(k-1)转移到第k时刻状态X^(i,t)(k)的概率，该转移概率可以从系统状态转移方程(1)中获得：

p(X^(i,t)(k|k-1)|X^(i,t)(k-1))＝p(X^(i,t)(k|k-1)-F·X^(i,t)(k-1))。

由状态转移方程(1)可得X^(i,t)(k|k-1)-F·X^(i,t)(k-1)服从N(0,C_V)表示的正态分布。

对应粒子的权值可以表示为：

${\mathrm{\ω}}^{(i,t)}\left(k\right)=p\left(Z\left(k\right)\right|X^{(i,t)}\left(k\right|k-1),{\hat{A}}_{N+1:M}^{t-1}),i=\mathrm{1,2}....,I,---\left(10\right)$

其中可以根据系统测量方程(2)获得：

$\begin{array}{c}{\mathrm{\ω}}^{(i,t)}\left(k\right)=p\left(Z\left(k\right)\right|X^{(i,t)}\left(k\right|k-1),{\hat{A}}_{N+1:M}^{t-1}),i=\mathrm{1,2}...,I\\ =\underset{m=1}{\overset{M}{\mathrm{\Π}}}p\left(z_m\left(k\right)\right|X_m^{(i,t)}\left(k\right|k-1),{\hat{A}}_m^{t-1})\\ =\underset{m=1}{\overset{M}{\mathrm{\Π}}}p\left(\right||{\hat{U}}^{(i,t)}\left(k\right|k-1)-{\hat{A}}_m^{t-1}||-Z_m\left(k\right))\end{array}---\left(11\right)$

其中表示在第t次EM迭代中，第i个用户轨迹样本中第k时刻用户估计位置与第m个AP位置之间的距离与测量距离之间的误差，由系统测量方程(2)可得该误差满足已知概率分布，在本实现样例中服从正态分布；位置U(k|k-1)为预测状态向量X(k|k-1)的一部分，因为 $X\left(k\right)={\left[\begin{array}{cc}U{\left(k\right)}^T& \stackrel{\·}{U}{\left(k\right)}^T\end{array}\right]}^T.$ 得到归一化粒子权值： ${\stackrel{\‾}{\mathrm{\ω}}}^{(i,t)}\left(k\right)=\frac{{\mathrm{\ω}}^{(i,t)}\left(k\right)}{\underset{i=1}{\overset{I}{\mathrm{\Σ}}}{\mathrm{\ω}}^{(i,t)}\left(k\right)},i=\mathrm{1,2}.....,I$

为了防止粒子退化，在上述迭代过程中按照下式对预测粒子重采样:

$P_r(X^{(i,t)}\left(k\right)=X^{(j,t)}\left(k\right|k-1))={\stackrel{\‾}{\mathrm{\ω}}}^{(j,t)}\left(k\right)---\left(12\right)$

该表达式表示：按照粒子的权重，进行重采样，各生成的粒子以自己的权重大小为概率，出现在更新后的粒子中。于是经过重采样，权重较大的粒子，会在更新后的粒子中多次出现，而权重较小的粒子，会在更新后的粒子中出现频次较少，甚至完全消失。其中上标(j，t)表示第t次EM迭代步骤中第j个粒子，j表示更新前的粒子X^(j,t)(k|k-1)的粒子序号，i表示经过重采样后的更新粒子X^(i,t)(k)的粒子序号，则经过重采样后的更新粒子权重都相同，那么移动终端的第i次轨迹样本实现(或称作第i个轨迹粒子)可以表示为

U^(i,t)(1:k)为：U^(i,t)(1:k)＝[U^(j,t)(1:k-1),U^(i,t)(k)],i＝1,2.....,I (13)

由于移动终端在不同时间序列上得到的更新后的状态粒子是对应该时刻状态概率的独立同分布采样，相应地，由相应时间序列的状态粒子组成的移动终端轨迹样本实现(或轨迹粒子)也可以看成是服从轨迹概率分布的独立同分布采样实现，所以轨迹粒子可以看成是独立等概率的轨迹粒子，且这些轨迹粒子模拟了轨迹概率分布。

(c)将上述的EM算法和粒子滤波算法结合起来，也就是本可选实施例基于虚拟AP的EM和粒子滤波定位算法。将获得的轨迹粒子U^(i,t)(1:K)带入到EM算法中，那么调整后的EM算法为：

E-Step：

$\begin{array}{c}Q(A_{N+1:M},q^t\left(\hat{U}(1:K)\right))\text{=}{E}_{q^t\left(\hat{U}(1:K)\right)}p(Z(1:K),{\hat{U}}^t(1:K)|A_{N+1:M};A_{1:N})\\ =\underset{i=1}{\overset{I}{\mathrm{\Σ}}}p(Z(1:K),{\hat{U}}^{(i,t)}(1:K)|A_{N+1:M};A_{1N})\end{array}---\left(14\right)$

M-Step：

$\begin{array}{c}{\hat{A}}_m=\underset{A_m}{\arg \max }\left\{E_{q^t\left(\hat{U}(1:K)\right)}p(Z_m(1:K),\hat{U}(1:K)|A_m)\right\}\\ =\underset{A_m}{\arg \max }\left\{\underset{i=1}{\overset{I}{\mathrm{\Σ}}}p(Z_m(1:K),{\hat{U}}^{(i,t)}(1:K)|A_m)\right\}\\ =\underset{A_m}{\arg \max }\left\{\underset{i=1}{\overset{I}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}p(Z_m\left(k\right),{\hat{U}}^{(i,t)}\left(k\right)|A_m)\right\}\end{array}---\left(15\right)$

最终通过(15)完成了虚拟AP的位置确定。

第二步通过三角定位法，利用真实AP和虚拟AP以及终端与真实AP之间的距离来估计移动终端的轨迹。本步骤可以使用现有的最小二乘(Least Square，简称为LS)或者约束最小二乘(Constrained Least Square，简称为CLS)或者EKF算法实现。

为验证本可选实施例算法的性能，并与其他一些定位算法比较，进行了matlab仿真实验。图4是根据本发明实施例的真实AP和虚拟AP坐标以及多径示意图，仿真平台基于二维平面场景，如图4所示。在本次仿真中，目标移动台个数为1，在30m*30m的区域内运动(平面墙体)，三个固定射频AP位于相邻的角落，其坐标分别为(0,15),(30,15),(30,-15)。为说明简单，仿真只考虑一次反射多径(虚拟AP为一次镜像点)，这样仿真环境中就产生了6个虚拟AP，相应的理论坐标为(60,15),(0,-45),(-30,15),(0,-45),(30,45),(-30，-15)，虚拟AP的坐标我们可以使用本发明算法进行估计。3个AP和6个虚拟AP如图2所示，其中黑色实线表示直射径，黑色虚线表示非直射径。假定移动终端的初始位置为(0,0)，初始速度为[1，0]/s，移动速度扰动是均值为0方差为0.2m/s²的随机噪声，共计测量50次直射径和非直射径的信息，且测量时间间隔为0.5秒。

均方定位误差(MSPE)是衡量定位算法性能优劣标准之一，实验是在独立仿真100次后，最终得出的平均数据。图3横坐标代表测量距离的平均噪声纵坐标代表均方定位误差(MSPE)。本次实验对比了四种不同的定位算法，包括基于本发明EM和粒子滤波的约束最小二乘(CLS)定位算法和扩展卡尔曼(EKF)定位算法，以及传统的最小二乘(CLS)定位算法和扩展卡尔曼(EKF)定位算法。

图5是根据本发明实施例的MSPE的仿真结果示意图，如图5所示，基于EM和粒子滤波的约束最小二乘(CLS)定位算法和扩展卡尔曼(EKF)定位算法的定位精度和性能明显要好于传统的最小二乘(CLS)定位算法和扩展卡尔曼(EKF)定位算法。因为可选实施例的定位算法不仅利用到了直射径的定位信息，而且也利用到了非直射径的定位信息，克服了传统方法利用非直射径需要严格的假设的问题，减少了非直射径参数的时变特性影响，所以其定位精度要比传统的算法要高。

下面结合一个可选实施例进行详细说明。

首先给定本方案实施的初始条件：给定N个真实AP的位置A_1:N，移动终端的初始状态 $X\left(0\right)={\left[\begin{array}{cc}U{\left(0\right)}^T& \tilde{U}{\left(0\right)}^T\end{array}\right]}^T$ (初始位置和初始速度)，共计获得K个时刻点的多径测量距离Z(1:K)，测量时间间隔△t，同时已知下列统计参数：测量误差矩阵σ，和用户移动速度扰动矩阵C_V。

根据系统的状态转移方程(1)和测量方程(2)，采用本发明中基于EM和粒子滤波的位置估计算法首先实现对虚拟AP位置估计。这里主要介绍算法的第t步EM迭代过程，主要包括下面两步。

(a)基于粒子滤波的E(Estimation)步骤

假设在t-1次迭代步骤中的虚拟AP位置估计，在E(Estimation)步骤中，采用粒子滤波方法生成能够模拟移动终端轨迹概率分布q^t(U(1:K))的轨迹粒子U^(i,t)(1:K)，令粒子数为I＝100，过程如下：

(a.1)状态预测粒子生成

同样考虑到粒子滤波方法也是迭代算法，以粒子滤波迭代步骤中的第k步(对应第k时刻)迭代过程介绍为例：令第k-1次迭代得到的状态粒子为X^(i,t)(k-1)，则在第k时刻按照重要性密度函数θ(X^(i,t)(k)|X^(i,t)(k-1),Z(k))生成I＝100个该时刻预测状态粒子，记为在SIS粒子滤波方法中该重要性密度函数按照下式选取：

θ(X^(i,t)(k)|X^(i,t)(k-1),Z(k))＝p(X^(i,t)(k)|X^(i,t)(k-1))

＝p(X^(i,t)(k)-F·X^(i,t)(k-1))

其中由系统状态转移方程(1)可知，变量X^(i,t)(k)-F·X^(i,t)(k-1)服从N(0,C_V)表示的正态分布。

图6是根据本发明实施例的状态预测粒子生成示意图，为简单示意状态粒子生成过程，图6中给出第t次EM迭代步骤中的k＝2时刻的位置预测粒子生成示意图，为简单描述，其中粒子数取I＝3。

(a.2)计算状态预测粒子权重

则各粒子对应的权重值可以表示为：

$\begin{array}{c}{\mathrm{\ω}}^{(i,t)}\left(k\right)=p\left(Z\left(k\right)\right|X^{(i,t)}\left(k\right|k-1),{\hat{A}}_{N+1:M}^{t-1}),i=\mathrm{1,2}...,I\\ =\underset{m=1}{\overset{M}{\mathrm{\Π}}}p\left(z_m\left(k\right)\right|X_m^{(i,t)}\left(k\right|k-1),{\hat{A}}_m^{t-1})\\ =\underset{m=1}{\overset{M}{\mathrm{\Π}}}p\left(\right||{\hat{U}}^{(i,t)}\left(k\right|k-1)-{\hat{A}}_m^{t-1}||-Z_m\left(k\right))\end{array}$

其中变量表示在第t次EM迭代中，第i个用户轨迹样本中第k时刻用户估计位置与第m个AP位置之间的距离与测量距离之间的误差，由系统状态测量方程(2)可得该误差满足已知概率分布，在本实现样例中服从正态分布；

得到归一化粒子权值： ${\stackrel{\‾}{\mathrm{\ω}}}^{(i,t)}\left(k\right)=\frac{{\mathrm{\ω}}^{(i,t)}\left(k\right)}{\underset{i=1}{\overset{I}{\mathrm{\Σ}}}{\mathrm{\ω}}^{(i,t)}\left(t\right)},i=\mathrm{1,2}.....,I$

图7是根据本发明实施例的状态预测粒子权重计算示意图，为简单示意计算状态预测粒子全重，图7中给出第t次EM迭代步骤中的k＝2时刻的位置预测粒子全重计算示意图，为简单描述，其中粒子数取I＝3。

(a.3)重采样更新状态粒子和生成轨迹粒子

为了防止粒子退化，按照下式进行粒子重采样得到第k时刻的更新粒子

$p(X^{(i,t)}\left(k\right)=X^{(j,t)}\left(k\right|k-1))={\stackrel{\‾}{\mathrm{\ω}}}^{(j,t)}\left(k\right)$

其中j表示第i个更新粒子X^(i,t)(k)重采样到的预测粒子X^(j,t)(k|k-1)对应的粒子序号，并且重采样的粒子权重都相同，那么移动终端的对应时刻1,2,…,k的位置轨迹U^(i,t)(1:k)为：

U^(i,t)(1:k)＝[U^(j,t)(1:k-1),U^(i,t)(k)],i＝1,2.....,I，

其中位置U(k)为状态向量X(k)的一部分，因为 $X\left(k\right)={\left[\begin{array}{cc}U{\left(k\right)}^T& \tilde{U}{\left(k\right)}^T\end{array}\right]}^T,$ j表示重采样得到第i个更新粒子X^(i,t)(k)取到预测粒子X^(j,t)(k|k-1)的粒子序号。

对应每次迭代步骤中经过重采样得到的更新粒子X^(i,t)(k)实际上是对应式(6)概率分布的独立同分布采样，每个状态粒子X^(i,t)(k)具有相同的概率。因此，经过粒子滤波得到的轨迹粒子同样可以看成是独立同分布实现，该分布近似等于式(1)中潜变量U(1:K)的概率分布q^t(U(1:K))。

图8是根据本发明实施例的状态粒子重采样示意图，如图8所示，示意了重采样更新状态粒子过程，该图描述第t次EM迭代步骤中的k＝2时刻的位置粒子更新过程，粒子数取I＝3。

图9是根据本发明实施例的轨迹粒子生成示意图，如图9所示，示意了轨迹粒子生成，该图描述第t次EM迭代步骤中的轨迹粒子U^(i,t)(1:2)生成过程，粒子数取I＝3。

图10是根据本发明实施例的EM算法第t次迭代中粒子滤波示意图，为简单示意粒子滤波的迭代过程，图10给出了第t次EM迭代步骤中的粒子滤波(K＝3，粒子数I＝3)的过程。

(a.4)基于轨迹粒子修正EM算法E步骤

然后将粒子滤波得到的移动终端运动轨迹粒子代入到EM算法E步骤中

修正E-Step

$\begin{array}{c}Q(A_{N+1:M},q^t\left(\hat{U}(1:K)\right))=E_{q^t\left(\hat{U}(1:K)\right)}p(Z(1:K),{\hat{U}}^t(1:K)|A_{N+1:M};A_{1:N})\\ =\underset{i=1}{\overset{I}{\mathrm{\Σ}}}p(Z(1:K),{\hat{U}}^{(i,t)}(1:K)|A_{N+1:M};A_{1:N})\end{array}---\left(13\right)$

(b)采用粒子滤波修正EM算法M(Maximization)步骤

进一步，将粒子滤波得到的移动终端运动轨迹粒子代入到EM算法M步骤中：

修正M-Step：

$\begin{array}{c}{\hat{A}}_m=\underset{A_m}{\arg \max }\left\{E_{q^t\left(\hat{U}(1:K)\right)}p(Z{(1:K)}_m,\hat{U}(1:K)|A_m)\right\}\\ =\underset{A_m}{\arg \max }\left\{\underset{i=1}{\overset{I}{\mathrm{\Σ}}}p(Z{(1:K)}_m,{\hat{U}}^{(i,t)}(1:K)|A_m)\right\}\end{array}---\left(14\right)$

由于各虚拟AP位置可以通过EM算法独立求解，那么将上式进行离散搜索

$\begin{array}{c}{\hat{A}}_m=\underset{A_m}{\arg \max }\left\{\underset{i=1}{\overset{I}{\mathrm{\Σ}}}p(Z{(1:K)}_m,{\hat{U}}^{(i,t)}(1:K)|A_m)\right\}\\ =\underset{A_m}{\arg \max }\left\{\underset{i=1}{\overset{I}{\mathrm{\Σ}}}p\left(Z{(1:K)}_m\right|A_m,{\hat{U}}^{(i,t)}(1:K))\right\}\\ =\underset{A_m}{\arg \max }\left\{\underset{k}{\mathrm{\Σ}}\underset{i=1}{\overset{I}{\mathrm{\Σ}}}\log \right[p\left(Z_m\left(k\right)\right|A_m,{\hat{U}}^{(i,t)}\left(k\right))\left]\right\}\\ =\underset{{\left\{A_m^l\right\}}_{l=1}^L}{\arg \max }\left\{\underset{k}{\mathrm{\Σ}}\underset{i=1}{\overset{I}{\mathrm{\Σ}}}\log \left({W_m}^{\′}(k,i,L)\right)\right\}\end{array}---\left(15\right)$

其中表示在第t次EM迭代中，第i个用户轨迹样本中第k时刻用户估计位置与第m个AP位置之间的距离与测量距离之间的误差，该误差满足已知概率分布，在本实现样例中服从正态分布；表示虚拟AP的L个可能的位置，这些可能的位置可以通过对虚拟AP可能范围进行离散化得到，比如在极坐标系统中对虚拟AP位置的可能半径范围和角度范围分别进行离散化后得到。

图11是根据本发明实施例的基于EM和粒子滤波的虚拟AP位置估计方案流程图，上述基于EM和粒子滤波的虚拟AP位置估计方案流程图如图11所示，如图11所示，运用了期望运算(即E步骤)和最大化运算(即M步骤)。

(c)基于位置已知AP的传统定位方法

根据测量的距离，结合真实AP和虚拟AP的坐标，采用传统的CLS定位算法和EKF定位算法来获得移动终端的运动轨迹，即可完成定位过程。

综上所述，本可选实施例采用基于EM和粒子滤波的射频多径定位方案，引入了虚拟AP的技术方案。首次将非直射径转化为移动终端与相应虚拟AP的直射径，根据移动终端位置与虚拟AP位置之间的独立性，使用EM算法和粒子滤波相结合的思想来估计虚拟AP的位置，最后根据真实AP和虚拟AP来实现对移动终端的定位。有效地利用了非直射径的定位信息，克服了传统方法利用非直射径需要严格的假设的问题，减少了非直射径参数的时变特性影响，显著提高了移动终端定位精度和系统性能。

在另外一个实施例中，还提供了一种软件，该软件用于执行上述实施例及优选实施方式中描述的技术方案。

在另外一个实施例中，还提供了一种存储介质，该存储介质中存储有上述软件，该存储介质包括但不限于：光盘、软盘、硬盘、可擦写存储器等。

显然，本领域的技术人员应该明白，上述的本发明的各模块或各步骤可以用通用的计算装置来实现，它们可以集中在单个的计算装置上，或者分布在多个计算装置所组成的网络上，可选地，它们可以用计算装置可执行的程序代码来实现，从而，可以将它们存储在存储装置中由计算装置来执行，并且在某些情况下，可以以不同于此处的顺序执行所示出或描述的步骤，或者将它们分别制作成各个集成电路模块，或者将它们中的多个模块或步骤制作成单个集成电路模块来实现。这样，本发明不限制于任何特定的硬件和软件结合。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (12)

1.一种终端的定位方法，其特征在于，包括：

获取终端的物理无线访问接入点AP的位置信息和虚拟AP的位置信息，其中，所述虚拟AP的位置信息为所述物理AP的镜像位置信息；

根据所述物理AP的位置信息和所述虚拟AP的位置信息按照预设算法对所述终端进行定位。

2.根据权利要求1所述的方法，其特征在于，获取所述终端的虚拟AP的位置信息包括：

通过最大期望运算和迭代算法根据所述终端的估计位置信息、所述终端与所述物理AP和/或所述虚拟AP的距离信息获取所述虚拟AP的位置信息。

3.根据权利要求2所述的方法，其特征在于，获取所述虚拟AP的位置信息包括：通过以下公式确定所述虚拟AP的位置信息的最大化似然估计值，根据该最大化似然估计值确定所述虚拟AP的位置信息：

${\hat{A}}_{N+1:M}=\underset{A_{N+1:M}}{\arg \max }\left\{\underset{U(1:K)}{\mathrm{\Σ}}p(Z(1:K),U(1:K)|A_{N+1:M};A_{1:N})\right\}$

其中，上标^表示估计运算，A_1:N是N个所述AP的位置坐标，A_N+1:M是M-N个所述虚拟AP的位置坐标，U(1:K)为所述终端的移动路径上从第1时刻到第K时刻各时间点的位置信息，Z(1:K)为所述终端的移动路径上从第1时刻到第K时刻各时间点到所述物理AP或者所述虚拟AP的距离信息，p(Z(1:K),U(1:K)|A_N+1:M；A_1:N)表示在给定A_1:N和A_N+1:M的条件下，发生Z(1:K)和U(1:K)的概率，N、M为自然数。

4.根据权利要求3所述的方法，其特征在于，根据所述最大化似然估计值确定所述虚拟AP的位置信息包括：对所述最大化似然估计值按照以下公式进行期望运算和最大化运算，得到所述虚拟AP的位置信息：

所述期望运算通过以下公式实现：

所述U(1:K)的概率分布q^t(U(1:K))为： $q^t\left(\hat{U}(1:K)\right)=p\left(\hat{U}(1:K)\right|Z(1:K),{\hat{A}}_{N+1:M}^{t-1};A_{1:N}),$

其中，上标t和t-1分别表示第t次和第t-1次迭代步骤；

对于所述U(1:K)的最大似然函数的均值为：

$Q(A_{N+1:M},q^t\left(\hat{U}(1:K)\right))=E_{q^t\left(\hat{U}(1:K)\right)}p(Z(1:K),{\hat{U}}^t(1:K)|A_{N+1:M};A_{1:N});$

通过迭代运算，如下公式所示，得到使得最大化虚拟AP位置的估计值：

$\begin{array}{c}{\hat{A}}_{N+1:M}=\underset{A_{N+1:M}}{\arg \max }\left\{Q(A_{N+1:M},q^t\left(\hat{U}(1:K)\right))\right\}\\ =\underset{A_{N+1:M}}{\arg \max }\left\{\underset{m=N+1}{\overset{M}{\mathrm{\Σ}}}\log E_{q^t\left(\hat{U}(1:K)\right)}p(Z{(1:K)}_m,{\hat{U}}^t(1:K)|A_m)\right\}\end{array};$

所述虚拟AP位置的估计值表示为

${\hat{A}}_m=\underset{A_m}{\arg \max }\left\{E_{q^t\left(\hat{U}(1:K)\right)}p(Z{(1:K)}_m,{\hat{U}}^t(1:K)|A_m)\right\}m=N+1:M.$

5.根据权利要求4所述的方法，其特征在于，所述通过轨迹粒子U^(i,t)(1:K)确定：

U^(i,t)(1:K)＝[U^(j,t)(1:K-1),U^(i,t)(K)],i＝1,2.....,I，其中，为所述终端作为第i个粒子在第k个时刻、第t次迭代过程中的位置，U^(i,t)(1:K)为从第1时刻到第K时刻所有的集合，I为自然数，为在粒子滤波算法中生成的粒子总数。

6.根据权利要求5所述的方法，其特征在于，通过如下公式获取所述虚拟AP的位置信息：

${\hat{A}}_m=\underset{A_m}{\arg \max }\left\{\underset{i=1}{\overset{I}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}p(Z_m\left(k\right),{\hat{U}}^{(i,t)}\left(k\right)|A_m)\right\},$

其中，Z_m(k)为所述终端在第k个时刻到第m个所述物理AP或者所述虚拟AP的直射径测量距离，Z(1:K)为所述终端的移动路径上从第1时刻到第K时刻各时间点到所述物理AP或者所述虚拟AP的距离信息的集合。

7.一种终端的定位装置，其特征在于，包括：

获取模块，用于获取终端的物理无线访问接入点AP的位置信息和虚拟AP的位置信息，其中，所述虚拟AP的位置信息为所述物理AP的镜像位置信息；

定位模块，用于根据所述物理AP的位置信息和所述虚拟AP的位置信息按照预设算法对所述终端进行定位。

8.根据权利要求7所述的装置，其特征在于，所述获取模块还用于通过最大期望运算和迭代算法根据所述终端的估计位置信息、所述终端与所述物理AP和/或所述虚拟AP的距离信息获取所述虚拟AP的位置信息。

9.根据权利要求8所述的装置，其特征在于，所述获取模块还用于通过以下公式确定所述虚拟AP的位置信息的最大化似然估计值，根据该最大化似然估计值确定所述虚拟AP的位置信息：

${\hat{A}}_{N+1:M}=\underset{A_{N+1:M}}{\arg \max }\left\{\underset{U(1:K)}{\mathrm{\Σ}}p(Z(1:K),U(1:K)|A_{N+1:M};A_{1:N})\right\}$

其中，上标^表示估计运算，A_1:N是N个所述AP的位置坐标，A_N+1:M是M-N个所述虚拟AP的位置坐标，U(1:K)为所述终端的移动路径上从第1时刻到第K时刻各时间点的位置信息，Z(1:K)为所述终端的移动路径上从第1时刻到第K时刻各时间点到所述物理AP或者所述虚拟AP的距离信息，p(Z(1:K),U(1:K)|A_N+1:M；A_1:N)表示在给定A_1:N和A_N+1:M的条件下，发生Z(1:K)和U(1:K)的概率，N、M为自然数。

10.根据权利要求9所述的装置，其特征在于，所述获取模块还用于对所述最大化似然估计值按照以下公式进行期望运算和最大化运算，得到所述虚拟AP的位置信息：

所述期望运算通过以下公式实现：

所述U(1:K)的概率分布q^t(U(1:K))为： $q^t\left(\hat{U}(1:K)\right)=p\left(\hat{U}(1:K)\right|Z(1:K),{\hat{A}}_{N+1:M}^{t-1};A_{1:N}),$

其中，上标t和t-1分别表示第t次和第t-1次迭代步骤；

对于所述U(1:K)的最大似然函数的均值为：

$Q(A_{N+1:M},q^t\left(\hat{U}(1:K)\right))=E_{q^t\left(\hat{U}(1:K)\right)}p(Z(1:K),{\hat{U}}^t(1:K)|A_{N+1:M};A_{1:N});$

通过迭代运算，如下公式所示，得到使得最大化虚拟AP位置的估计值：

$\begin{array}{c}{\hat{A}}_{N+1:M}=\underset{A_{N+1:M}}{\arg \max }\left\{Q(A_{N+1:M},q^t\left(\hat{U}(1:K)\right))\right\}\\ =\underset{A_{N+1:M}}{\arg \max }\left\{\underset{m=N+1}{\overset{M}{\mathrm{\Σ}}}\log E_{q^t\left(\hat{U}(1:K)\right)}p(Z{(1:K)}_m,{\hat{U}}^t(1:K)|A_m)\right\}\end{array},$

所述虚拟AP位置的估计值表示为

${\hat{A}}_m=\underset{A_m}{\arg \max }\left\{E_{q^t\left(\hat{U}(1:K)\right)}p(Z{(1:K)}_m,{\hat{U}}^t(1:K)|A_m)\right\}m=N+1:M.$

11.根据权利要求10所述的装置，其特征在于，所述通过轨迹粒子U^(i,t)(1:K)确定：

U^(i,t)(1:K)＝[U^(j,t)(1:K-1),U^(i,t)(K)],i＝1,2.....,I，其中，为所述终端作为第i个粒子在第k个时刻、第t次迭代过程中的位置，U^(i,t)(1:K)为从第1时刻到第K时刻所有的集合，I为自然数，为在粒子滤波算法中生成的粒子总数。

12.根据权利要求11所述的装置，其特征在于，通过如下公式获取所述虚拟AP的位置信息：

${\hat{A}}_m=\underset{A_m}{\arg \max }\left\{\underset{i=1}{\overset{I}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}p(Z_m\left(k\right),{\hat{U}}^{(i,t)}\left(k\right)|A_m)\right\},$

其中，Z_m(k)为所述终端在第k个时刻到第m个所述物理AP或者所述虚拟AP的直射径测量距离，Z(1:K)为所述终端的移动路径上从第1时刻到第K时刻各时间点到所述物理AP或者所述虚拟AP的距离信息的集合。