技术领域
[0001] 本
发明涉及定位测量领域，尤其涉及一种高精度室内无线定位方法。
背景技术
[0002] 随着无线通信和移动网络的发展，人们对定位的需求日益增加，但是由于
信号遮挡等因素的存在，使得如北斗、GPS等通用的定位技术在室内环境下难以得到部署。
[0003] 现有针对室内的无线定位技术主要有超宽带(UWB)、WiFi和Beacon定位技术等，每种技术均有特定的适用场景。上述技术的基本思路是通过计算“锚
节点”(已知布置坐标的定位节点)的信号到达时间或接收功率，获得被测节点和锚节点的距离，然后通过三
角定位法、质心
算法或最小二乘估计算法确定被测节点的
位置信息。
[0004] 其中，三点定位法是利用不在同一条直线上的三个锚节点到未知节点的距离，通过解方程得到未知节点的坐标。质心算法是选出距离被测节点较近的多个锚节点，将其坐标依据到被测节点的距离进行加权平均。最小二乘估计算法则是将被测节点到多个锚节点的距离聚合成向量，将对应的锚节点坐标聚合成矩阵，通过求矩阵的伪逆得到待测节点的坐标估计。但是，上述三种方法均存在
缺陷：
[0005] 1、三点定位算法需要解方程组，但由于测量误差的存在，经常无法得到唯一解，需要进行近似，影响定位精度。
[0006] 2、质心算法在锚节点数量很少的情况下定位精度较差。
[0007] 3、最小二乘估计需要计算矩阵的伪逆，运算复杂度高，并且在矩阵不满秩的情况下无法求解。
[0008] 4、在实际环境中待测物体具有通常处于连续的移动状态，此时刻的坐标通常与前一时刻的位置信息存在较大的相关性，而上述三种方法均未有效利用此相关性。
发明内容
[0009] 本发明的目的是提供一种高精度室内无线定位方法，能够克服上述三种室内无线定位方法中的弊端，高效准确的获得室内待测节点的高精度位置信息。
[0010] 本发明采用下述技术方案：
[0011] 一种高精度室内无线定位方法，依次包括以下步骤：
[0012] A：设在t时刻被测节点可利用的锚节点个数为N，N≥3，根据每个锚节点的序列号得到对应的锚节点的坐标(xn,yn,zn)，并随机将N个锚节点编号为[1:N]；根据到达时间或接收功率得到被测节点在t时刻到第n个锚节点的距离 1≤n ≤N；
[0013] B：设t时刻被测节点的真实坐标为(x(t),y(t),z(t))，并定义为向量形式α(t)；
[0014] (x(t) y(t) z(t))T@α(t) (1)；
[0015] 将步骤A中在t时刻获得的所有N个锚节点到被测节点的距离 及锚节点的坐标(xn,yn,zn)，按(2)式计算并
整理为向量r(t)；
[0016]
[0017] 其中， 为向量r(t)中第n个元素；
[0018] 将t时刻获得的所有N个锚节点坐标，按(3)式计算并整理为矩阵Φ(t)；
[0019]
[0020] 其中， 为矩阵Φ(t)中第n行向量， 为矩阵Φ(t)中第n行第1列元素； 为矩阵Φ(t)中第n行第2列元素； 为矩阵Φ(t)中第n行第3列元素；
[0021] 由距离和坐标的运算关系
[0022] 根据(1)式-(3)式，可得到
[0023] r(t)＝Φ(t)α(t)+ω(t) (4)；
[0024] 其中，ω(t)表示实际测量过程中存在的误差向量；设误差服从均值为零、方差为σI的高斯分布，记为N(ω(t)；0,σI)，其中矩阵σI为大小为N×N的协方差矩阵；
[0025] 通过卡尔曼滤波挖掘相邻时刻的空间相关性，将t时刻在x、y、z三个坐标上的卡尔曼
状态方程系数定义为 和 并将所有t＝{1,...,T} 时刻的上述向量和矩阵归纳为如下集合形式，{α(1),...,α(T)}@α，{r(1),...,r(T)}@r， {Φ(1),...,Φ(T)}@Φ，
[0026] C：根据(4)式所示的运算关系，利用全概率公式、变量间的隐
马尔科夫特性和卡尔曼状态方程，对全局概率分布P(r,α,Φ,A,B)进行因式分解：
[0027]
[0028] (5)式中表达式P(g)表示概率分布；似然函数 可分解为定义为函数
表示数学期望为 方差为σI的高斯分布；由于连续的两次观测之间具有马尔可夫特性，表示为 u为服从标准高斯分布的噪声，因此在(5) 式中：
[0029] 坐标x(t),x(t-1)之间的函数约束为 定义为函数
[0030] 坐标y(t),y(t-1)之间的函数约束为 定义为函数
[0031] 坐标z(t),z(t-1)之间的函数约束为 定义为函数
[0032] 参数A和B的先验均假设为均匀分布，为P(A)＝U(0,1)，P(B)＝U(0,1)；其中U(a,b) 表示区间为[a,b]的均匀分布；
[0033] D：利用步骤C中得到的因子分解进行因子图建模；
[0034] 根据步骤C所示的因式分解进行因子图建模的过程为，定义(5)式中的每个函数为函数节点；定义每个变量为变量节点；将所有函数节点和与之相关的变量节点通过线段相连，便构成了因子图模型；
[0035] 引入映射节点 和对应的映射函数 映射节点 和对应的映射函数 以及映射节点 和对应的映射函数 其中δ(g)
表示delta函数；
[0036] 将上述因子图模型分为晶格网络部分、卡尔曼网络部分和卡尔曼参数估计部分，其中晶格网络部分对应于函数节点 到变量节点 之间的网络，卡尔曼网络部分对应于函数节点 到变量节点 之间的网络，卡尔曼参
数估计部分对应于变量节点 之间互相连接的网络；
[0037] E：利用步骤D中所得的因子图模型，通过消息传递算法对晶格网络部分、卡尔曼网络部分和卡尔曼参数估计部分进行消息计算，得到在t时刻被测节点的高精度位置信息。
[0038] 所述的步骤E包括以下具体步骤：
[0039] E1：通过消息传递算法对晶格网络部分进行消息计算；
[0040] 根据消息传递算法中的平均场计算规则分别得到 到 的消息到 的消息 以及 到 的消息
[0041]
[0042]
[0043]
[0044] 其中，变量 和 分别为坐标 和 的数学期望，其值可以由(14)式和(15)式得到，σ为测距噪声方差；
[0045] 根据消息传递算法，变量 到函数 的消息 计算为：
[0046]
[0047] 则映射函数 到节点x(t)的消息为：
[0048]
[0049] 同理，映射函数 到节点 的消息计算为：
[0050]
[0051] 映射函数 到节点 的消息计算为：
[0052]
[0053] E2：通过消息传递算法对卡尔曼网络部分进行消息计算；
[0054] 设函数 到变量x(t-1)的消息为 其值在(12) 式中更新，则利用置信传播规则计算消息 为
[0055]
[0056] (10)式中变量 和 分别定义为
[0057]
[0058] 设卡尔曼状态方程系数 的数学期望分别为 和 状态方程系数的计算方法分别在(19)式和(21)式中，则利用平均场规则计算得到中间函数
进而利用置信传播规则计算消息 为
[0059]
[0060] 其中， 分别表示消息 的期望和方差；由于定位算法不能根据t时刻的数据修正t-1时刻的位置估计，称之为无追溯性，则变量x(t)的置信与 到x(t)的消息无关，即
[0061]
[0062] 其中，b(x(t))表示变量x(t)的后验概率估计，定义为置信； 分别表示变量x(t)置信的期望和方差；
[0063] 同理，得到坐标y、z的置信
[0064]
[0065]
[0066] 在消息传递算法中，变量置信的物理意义为该变量的联合后验概率估计，即 (13)式—(15)式中得到了在t时刻被测节点的高精度位置信息；
[0067] E3：通过消息传递算法对卡尔曼参数估计部分进行消息计算；
[0068] 根据t和t-1时刻坐标x的置信b(x(t))和b(x(t-1))，由平均场规则可以计算消息为
[0069]
[0070] 其中，表达式b(x)表示以x为自变量对b(x)和f的乘积求积分；上式中变量 和分别表示消息 的期望和方差，可以计算为
[0071]
[0072] 假定连续Q个时刻状态方程系数不变，则可以得到自适应参数 的置信[0073]
[0074] 其中，自适应参数 的期望和方差 和 分别为
[0075]
[0076] 利用
期望最大化算法，将自适应参数 的置信 近似为
[0077]
[0078] 对于自适应参数 的估计，假设中间变量 则函数 到 的消息 根据消息传递算法计算为
[0079]
[0080] 其中表达式b(x)表示以x为自变量对b(x)和f的乘积求积分；上式中变量 定义为
[0081] 利用连续Q个时刻可以得到参数 的估计值为 从而得到变量的估计值
[0082]
[0083] 本发明能够克服上述三种室内无线定位方法中的弊端，有效提升测量精度，获得显著的性能增益，高效准确的获得室内待测节点的高精度位置信息。
附图说明
[0084] 图1为本发明步骤A中被测节点随时间发生连续的位移的示意图；
[0085] 图2为本发明
实施例中锚节点部署位置和待测设备移动轨迹的俯视图；
[0086] 图3为本发明实施例中锚节点部署位置和待测设备移动轨迹的三维视图；
[0087] 图4为本发明实施例中效锚节点个数N＝4时定位精度随噪声功率变化曲线图；
[0088] 图5为本发明实施例中效锚节点个数N＝6时定位精度随噪声功率变化曲线图；
[0089] 图6为本发明实施例中效锚节点个数N＝8时定位精度随噪声功率变化曲线图；
[0090] 图7为本发明实施例中效锚节点个数N＝12时定位精度随噪声功率变化曲线图；
[0091] 图8为发明实施例中定位精度随有效锚节点个数变化曲线图。
具体实施方式
[0092] 以下结合附图和实施例对本发明作以详细的描述：
[0093] 如图1至图8所示，本发明所述的高精度室内无线定位方法，依次包括以下步骤：
[0094] A：设在t时刻被测节点可利用的锚节点个数为N，N≥3，根据每个锚节点的序列号得到对应的锚节点的坐标(xn,yn,zn)，并随机将N个锚节点编号为[1:N]；根据到达时间或接收功率得到被测节点在t时刻到第n个锚节点的距离 1≤n ≤N；
[0095] 实施过程如附图1所示，其中五角星代表固定位置的锚节点，圆形代表被测节点，假定被测节点随时间发生连续的位移，移动路线如曲线所示。
[0096] B：设t时刻被测节点的真实坐标为(x(t),y(t),z(t))，并定义为向量形式α(t)；
[0097] (x(t) y(t) z(t))T@α(t) (1)；
[0098] 将步骤A中在t时刻获得的所有N个锚节点到被测节点的距离 及锚节点(t)的坐标(xn,yn,zn)，按(2)式计算并整理为向量r ；
[0099]
[0100] 其中， 为向量r(t)中第n个元素；
[0101] 将t时刻获得的所有N个锚节点坐标，按(3)式计算并整理为矩阵Φ(t)；
[0102]
[0103] 其中， 为矩阵Φ(t)中第n行向量， 为矩阵Φ(t)中第n行第1列元素； 为矩阵Φ(t)中第n行第2列元素； 为矩阵Φ(t)中第n行第3列元素；
[0104] 由距离和坐标的运算关系 根据(1)式- (3)式，可得到
[0105] r(t)＝Φ(t)α(t)+ω(t) (4)；
[0106] 其中，ω(t)表示实际测量过程中存在的误差向量；设误差服从均值为零、方差为σI的高斯分布，记为N(ω(t)；0,σI)，其中矩阵σI为大小为N×N的协方差矩阵；
[0107] 通过卡尔曼滤波挖掘相邻时刻的空间相关性，将t时刻在x、y、z三个坐标上的卡尔曼状态方程系数定义为 和 并将所有t＝{1,...,T} 时刻的上述向量和矩阵归纳为如下集合形式，{α(1),...,α(T)}@α，{r(1),...,r(T)}@r， {Φ(1),...,Φ(T)}@Φ，
[0108] C：根据(4)式所示的运算关系，利用全概率公式、变量间的隐马尔科夫特性和卡尔曼状态方程，对全局概率分布P(r,α,Φ,A,B)进行因式分解：
[0109]
[0110] (5)式中表达式P(g)表示概率分布；似然函数 可分解为定义为函数
表示数学期望为 方差为σI的高斯分布；由于连续的两次观测之间具有马尔可夫特性，表示为 u为服从标准高斯分布的噪声，因此在(5) 式中：
[0111] 坐标x(t),x(t-1)之间的函数约束为 定义为函数
[0112] 坐标y(t),y(t-1)之间的函数约束为 定义为函数
[0113] 坐标z(t),z(t-1)之间的函数约束为 定义为函数
[0114] 参数A和B的先验均假设为均匀分布，为P(A)＝U(0,1)，P(B)＝U(0,1)；其中U(a,b) 表示区间为[a,b]的均匀分布；
[0115] D：利用步骤C中得到的因子分解进行因子图建模；
[0116] 根据步骤C所示的因式分解进行因子图建模的过程为，定义(5)式中的每个函数为函数节点；定义每个变量为变量节点；将所有函数节点和与之相关的变量节点通过线段相连，便构成了因子图模型；
[0117] 为了表达清楚，此处引入映射节点 和对应的映射函数 映射节点 和对应的映射函数 以及映射节点 和对应的映射函数
其中δ(g)表示delta函数；
[0118] 将上述因子图模型分为晶格网络部分、卡尔曼网络部分和卡尔曼参数估计部分，其中晶格网络部分对应于函数节点 到变量节点 之间的网络，卡尔曼网络部分对应于函数节点 到变量节点 之间的网络，卡尔曼参
数估计部分对应于变量节点 之间互相连接的网络；
[0119] E：利用步骤D中所得的因子图模型，通过消息传递算法对晶格网络部分、卡尔曼网络部分和卡尔曼参数估计部分进行消息计算；
[0120] 步骤E包括以下具体步骤：
[0121] E1：通过消息传递算法对晶格网络部分进行消息计算；
[0122] 根据消息传递算法中的平均场计算规则分别得到 到 的消息、 到 的消息以及 到 的消息：
[0123]
[0124]
[0125]
[0126] 其中，变量 和 分别为坐标 和 的数学期望，其值可以由 (14)式和(15)式得到，σ为测距噪声方差；
[0127] 根据消息传递算法，变量 到函数 的消息 计算为：
[0128]
[0129] 则映射函数 到节点x(t)的消息为：
[0130]
[0131] 同理，映射函数 到节点 的消息为：
[0132]
[0133] 映射函数 到节点 的消息为：
[0134]
[0135] E2：通过消息传递算法对卡尔曼网络部分进行消息计算；
[0136] 设函数 到变量x(t-1)的消息为 其值在(12) 式中更新，则利用置信传播规则计算消息 为
[0137]
[0138] 其中，变量 和 分别定义为
[0139]
[0140] 设状态方程系数 的数学期望分别为 和 状态方程系数 的计算方法分别在(19)式和(21)式中，则利用平均场规则计算得到中间函数
进而利用置信传播规则计算消息 为
[0141]
[0142] 其中 分别表示消息 的期望和方差。由于定位算法不能根据t 时刻的数据修正t-1时刻的位置估计，称之为无追溯性，则变量x(t)的置信与 到 x(t)的消息无关，即
[0143]
[0144] 其中，b(x(t))表示变量x(t)的后验概率估计，定义为置信； 分别表示变量x(t)置信的期望和方差；
[0145] 同理，可以得到坐标y、z的置信
[0146]
[0147]
[0148] 在消息传递算法中，变量置信的物理意义为该变量的联合后验概率估计，即(13)式-(15)式中得到了在t时刻被测节点的高精度位置信息。
[0149] 由于在所有时刻针对x,y,z坐标的消息计算具有完全相同的形式，简洁起见在步骤E2和步骤E3部分仅给出t时刻针对x坐标的消息更新公式，而对其他时刻和y,z坐标的更新可以类比得到，仅仅需要将对应的坐标x替换为y、z，时间t 替换为其他时刻即可。
[0150] E3：通过消息传递算法对卡尔曼参数估计部分进行消息计算；
[0151] 由于本发明能自适应调整卡尔曼方法中状态方程系数，可以称之为自适应卡尔曼方法，其中系数 和 也称之为自适应系数。根据t和t-1时刻坐标x的置信b(x(t))和b(x(t-1))，由平均场规则可以计算消息 为
[0152]
[0153] 其中表达式b(x)表示以x为自变量对b(x)和f的乘积求积分。上式中变量 和分别表示消息 的期望和方差，可以计算为
[0154]
[0155] 假定连续Q个时刻状态方程系数不变，则可以得到自适应参数 的置信[0156]
[0157] 其中，自适应参数 的期望和方差 和 分别为
[0158]
[0159] 利用期望最大化算法，将自适应参数 的置信 近似为
[0160]
[0161] 对于自适应参数 的估计，假设中间变量 则函数 到 的消息可以根据消息传递算法计算为
[0162]
[0163] 其中表达式b(x)表示以x为自变量对b(x)和f的乘积求积分，上式中变量 定义为
[0164] 利用连续Q个时刻可以得到参数 的估计值为 从而得到变量的估计值
[0165]
[0166] 根据本发明提出的定位算法，可以建立仿真环境，对所提算法以及文献中已有算法进行对比。本实施例中建立仿真环境对本文所提算法以及文献中已有算法进行数值仿真和对比。仿真环境设定为一个高1m、长宽均为5m的空间，每间隔 1m部署一个锚节点(锚节点部署时存在一定的装配误差)。待测设备在x轴和y轴方向上线性移动，在z轴方向上以正弦方式周期移动，并且移动轨迹并非光滑，而是有一定的随机抖动，如附图2和附图3所示。在附图2和附图3中，锚节点的部署位置由星号表示，设备的移动轨迹由实线表示，其中附图2和附图3分别为对应的俯视图和三维视图。
[0167] 假定待测设备在移动时检测到距离最近的N'个锚节点，但其中可能存在某个锚节点受到遮挡导致测距发生突变，经筛选后假设有效的锚节点个数为N＝4:12。与本发明所提算法进行对比的算法有
基础的三点质心算法、质心算法以及最小二乘估计方法。本文所提方法在估计自适应参数时设定连续Q＝5次观测中状态方程具有相同系数。仿真中通过在真实距离 上
叠加高斯噪声模拟测量误差，噪声功率为-4dB至-16dB，定位精度由
最小均方误差(MSE)表示，计算如下：
[0168]
[0169] 附图4至附图8给出了所有算法的定位精度随测量误差的变化曲线，其中附图4至附图7分别给出了接收机接收到锚节点个数N＝4、6、8和12时的仿真结果，从中可以看出发明所提算法能够有效提升测量精度。附图8展示了测量精度随有效锚节点个数的变化曲线，从中可以看出在锚节点个数较少时，本发明所提算法并不能表现出明显的性能增益，而随观测到的锚节点个数增加，本发明所提算法定位精度迅速提升。总之，由仿真结果表明，相对于已知算法本发明提出的动态定位算法能够获得显著的性能增益，具有很高的实用价值。
