0 概述

位置信息在基于位置的服务和应用中起到关键作用。由于到达时差(Time Difference of Arrival, TDOA)定位方式不需要在发送的信号中加入时间戳, 具有硬件复杂度低、定位精度高的优势^[1], 因此其在室内定位中得到广泛应用。

TDOA定位实质上是对双曲线方程组的求解, 然而传统算法多存在不足, 例如:Taylor级数展开法^[2]需要选取一个与实际位置较为接近的展开点以保证算法的收敛; Chan算法^[3]在测量值误差服从高斯分布时能达到理论最优值, 但当测量值存在非视距误差时性能会显著下降; 扩展卡尔曼滤波^[4-5]是一种对非线性系统进行线性近似处理的目标跟踪算法, 但跟踪精度容易受信道条件的影响; 加权最小二乘法^[6]和约束总体最小二乘法^[7]等基于最小二乘的改进算法在运算过程中存在矩阵求逆, 当测量值较少或测量误差较大时易产生矩阵奇异的情况。

除了上述传统解析算法以外, 粒子群优化算法^[8]、遗传算法^[9]和樽海鞘群算法^[10]等智能优化算法也被应用于TDOA定位。此类算法在搜索范围内初始化大量随机点, 通过建立随机点与测量值之间的适应度函数来评价其优劣程度, 再通过算法的寻优机制迭代更新随机点的位置, 最终收敛于最优的目标估计位置。智能优化算法省略了复杂的求解过程, 不存在解析算法不可导或无解的问题。由NFL (No-Free-Lunch)理论^[11]可知, 少有适用于所有优化问题的智能优化算法, 例如在应用过程中可能存在后期收敛速度慢、陷入局部最优的问题。对此, 研究者相继提出解决方案:文献[12]提出自适应权重策略以平衡算法的全局搜索和局部开发能力, 提高收敛速度; 文献[13]通过引入变异机制, 使粒子在优化过程中向不同方向运动, 从而保证算法可以跳出局部最优; 文献[14]提出将人工蜂群算法与模糊C均值聚类算法相结合的改进算法, 提高了寻优精度。

2019年, HEIDARI等人提出哈里斯鹰优化(Harris Hawk Optimization, HHO)算法^[15], 该算法有较强的全局搜索能力, 并且需要调节的参数较少。文献[16]将长期记忆引入HHO算法, 使个体参考过去的经验进行运动, 增强了种群多样性。文献[17]采用动态控制参数减小HHO算法陷入局部最优的概率, 并通过变异算子进一步提高全局搜索效率。

本文提出一种改进的哈里斯鹰优化(Improved HHO, IHHO)算法用于TDOA定位, 通过对适应度函数和初始种群进行优化, 提高算法的寻优精度和收敛速度。

1 TDOA定位模型

TDOA定位几何模型如图 1所示, 其中, 标签周期性地向基站发送定位信号。

假设基站数量为N, 则定位信号到达基站i的时间为:

$ {t_i} = {\dot t_i} + {n_{{t_i}}} $

(1)

其中, i=1, 2, …, N, ${{\dot t}_i}$为信号到达时间真实值, n_ti为到达时间噪声误差, 服从零均值高斯分布。TDOA定位在测量出标签到两个不同基站的到达时间差值后, 将其乘以信号传播速度即可得到距离差测量值。定义r_i为标签到基站i的距离, r_{i, j}为标签到基站i与到基站j的距离差测量值, 表达式为:

$ \begin{array}{*{20}{l}} {{r_{i,j}} = {\rm{c}}({t_i} - {t_j}) = {\rm{c}}({{\dot t}_i} + {n_{{t_i}}} - {{\dot t}_j} - {n_{{t_j}}}) = }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {r_i} + {n_{{r_i}}} - {r_j} - {n_{{r_j}}}} \end{array} $

(2)

$ {{r_i} = {\rm{c}}{{\dot t}_i} = \sqrt {{{(x - {x_i})}^2} + {{(y - {y_i})}^2}} } $

(3)

$ {{n_{{r_i}}} = {\rm{c}}{n_{{t_i}}}} $

(4)

其中, i, j=1, 2, …, N, i≠j, c为电磁波的传播速度, (x, y)和(x_i, y_i)分别为标签和基站i的坐标, n_ri为距离噪声误差, 服从均值为0、方差为σ_r²的高斯分布。

根据标签到两个不同基站的距离差能建立唯一的一条双曲线, 然后凭借多基站建立双曲线方程组, 在此基础上, 利用数学方法求解即可得到未知标签的估计位置。

2 改进的哈里斯鹰优化算法 2.1 改进的适应度函数

定义基站1为主基站, 其余基站为从基站。采用最大似然法估计标签位置$(\hat x, \hat y)$, r_{i, 1}服从均值为(r_i－r₁)、方差为σ_d²=2σ_r²的高斯分布, 且各个测量值相互独立, 则似然函数为:

$ {\kern 1pt} \begin{array}{*{20}{l}} {L = \prod\limits_{i = 2}^N {\left\{ {\frac{1}{{\sqrt {2{\rm{ \mathsf{ π} }}} {\sigma _{\rm{d}}}}}\exp \left[ { - \frac{{{{({r_{i,1}} - {r_i} + {r_1})}^2}}}{{2\sigma _{\rm{d}}^2}}} \right]} \right\}} = }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {{\left( {\frac{1}{{\sqrt {2{\rm{ \mathsf{ π} }}} {\sigma _{\rm{d}}}}}} \right)}^{N - 1}}\exp \left[ { - \frac{{\sum\limits_{i = 2}^N {{{({r_{i,1}} - {r_i} + {r_1})}^2}} }}{{2\sigma _{\rm{d}}^2}}} \right]} \end{array} $

(5)

求解使似然函数最大的坐标值, 等价于求解下式:

$ (\hat x,\hat y) = {\rm{argmin}} [\sum\limits_{i = 2}^N {{{({r_{i,1}} - {r_i} + {r_1})}^2}} ] $

(6)

当该非线性函数方程值最小时, 得到标签坐标的估计值, 用解析法求解比较困难。因此, 本文采用哈里斯鹰优化算法求解。利用式(7)推导出基站1为主基站时的适应度函数:

$ {f_1}(\mathit{\boldsymbol{X}}) = \sum\limits_{i = 2}^N {{{({r_{i,1}} - {r_i} + {r_1})}^2}} $

(7)

其中, X 为哈里斯鹰个体的位置矢量。适应度函数的选取对算法的寻优精度有直接影响。结合式(2)和式(7)可以发现, 主基站的信号到达时间t₁会影响适应度函数中的所有平方项, 而从基站的信号到达时间t_i只会影响其中第(i－1)个平方项, 这使得适应度函数对主基站的信号到达时间较为敏感, 尤其是当主基站测量误差较大时, 适应度值会明显增大, 此时适应度函数会失真, 不能很好地反映解的优劣程度。

定义基站i作为主基站时的适应度函数为:

$ {f_i}(\mathit{\boldsymbol{X}}) = \sum\limits_{j = 1}^N {{{({r_{j,i}} - {r_j} + {r_i})}^2}} $

(8)

其中, j≠i, i∈{1, 2, …, N}。首先将每个基站作为主基站, 按照最大似然法计算一个适应度函数值, 然后将其中的最小值作为改进的适应度函数f(X)代入非线性函数方程进行求解, 表达式为:

$ {f(\mathit{\boldsymbol{X}}) = \min [{f_1}(\mathit{\boldsymbol{X}}),{f_2}(\mathit{\boldsymbol{X}}), \cdots ,{f_N}(\mathit{\boldsymbol{X}})]} $

(9)

$ {(\hat x,\hat y) = {\rm{argmin}} [f(\mathit{\boldsymbol{X}})]} $

(10)

2.2 改进的初始种群

定义种群中个体的位置表示为 X =[x₁, x₂], 搜索空间的上界为 u^b=[u₁^b, u₂^b], 下界为 l^b=[l₁^b, l₂^b]。随机生成初始种群中个体的位置:

$ \mathit{\boldsymbol{X}} = {\rm{rand }}(1,2) \times ({\mathit{\boldsymbol{u}}^{\rm{b}}} - {\mathit{\boldsymbol{l}}^{\rm{b}}}) + {\mathit{\boldsymbol{l}}^{\rm{b}}} $

(11)

其中, rand(1, 2)表示生成2维随机向量, 其中元素为[0, 1]之间的随机数。

初始种群的分布对算法的收敛速度有很大的影响。由于没有任何先验知识, 因此大部分智能优化算法的初始位置都随机生成。利用混沌映射产生初始种群对提高收敛速度有一定帮助^[18], 但对于TDOA定位问题不是最佳解决方案。本文通过计算复杂度较低的Chan算法^[3]计算出一个初始解, 然后从生成的初始种群中随机挑选出一个个体用初始解代替, 从而得到改进的初始种群。由于初始解距离最终优化结果较近, 因此该算法可以减少不必要的全局搜索, 在不影响种群多样性的前提下, 达到快速收敛的目的。

2.3 哈里斯鹰优化算法

哈里斯鹰优化算法^[15]是一种模拟哈里斯鹰捕食行为的智能优化算法, 主要由搜索阶段、搜索与开发的转换和开发阶段三部分组成。

2.3.1 搜索阶段

哈里斯鹰随机栖息在某个地方, 通过2种策略找到猎物, 可表示为:

$ \begin{array}{*{20}{l}} {X(\tau + 1) = }\\ {\left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{X}}_{{\rm{ rand }}}}(\tau ) - {r_1}|{\mathit{\boldsymbol{X}}_{{\rm{ rand }}}}(\tau ) - 2{r_2}\mathit{\boldsymbol{X}}(\tau )|,q \ge 0.5}\\ {[{\mathit{\boldsymbol{X}}_{{\rm{ rabbit }}}}(\tau ) - {\mathit{\boldsymbol{X}}_{\rm{m}}}(\tau )] - {r_3}[{\mathit{\boldsymbol{l}}^{\rm{b}}} + {r_4}({\mathit{\boldsymbol{u}}^{\rm{b}}} - {\mathit{\boldsymbol{l}}^{\rm{b}}})],q < 0.5} \end{array}} \right.} \end{array} $

(12)

在式(12)中, X(τ)和 X(τ+1)分别为当前和下一次迭代时个体的位置, τ为迭代次数, X_rand(τ)为随机选出的个体位置, X_rabbit(τ)为猎物位置, 即拥有最优适应度的个体位置, r₁~r₄和q都是[0, 1]之间的随机数, q用于随机选择要采用的策略, X_m(τ)为个体平均位置, 表达式为:

$ {\mathit{\boldsymbol{X}}_{\rm{m}}}(\tau ) = \frac{1}{M}\sum\limits_{k = 1}^M {{\mathit{\boldsymbol{X}}_k}} (\tau ) $

(13)

其中, X_k(τ)为种群中第k个个体的位置, M表示种群规模。

2.3.2 搜索与开发的转换

HHO算法根据猎物的逃逸能量在搜索和不同的开发行为之间转换, 逃逸能量定义为:

$ E = 2{E_0}\left( {1 - \frac{\tau }{T}} \right) $

(14)

其中, E₀是猎物的初始能量, 为[－1, 1]之间的随机数, 每次迭代时自动更新, τ为迭代次数, T为最大迭代次数。当| E |≥1时, 进入搜索阶段; 当| E | < 1时, 进入开发阶段。

2.3.3 开发阶段

定义r为[0, 1]之间的随机数, 用于选择不同的开发策略。

当0.5≤ |E| < 1且r≥0.5时, 采取软围攻策略进行位置更新:

$ \mathit{\boldsymbol{X}}(\tau + 1) = \Delta \mathit{\boldsymbol{X}}(\tau ) - E|J{\mathit{\boldsymbol{X}}_{{\rm{rabbit}}}}(\tau ) - \mathit{\boldsymbol{X}}(\tau )| $

(15)

其中, Δ X(τ)= X_rabbit(τ)－ X(τ)表示猎物位置与个体当前位置的差值, J为[0, 2]之间的随机数。

当|E| < 0.5且r≥0.5时, 采取硬围攻策略进行位置更新:

$ \mathit{\boldsymbol{X}}(\tau + 1) = {\mathit{\boldsymbol{X}}_{{\rm{rabbit}}}}(\tau ) - E|\Delta \mathit{\boldsymbol{X}}(\tau )| $

(16)

当0.5≤ |E| < 1且r < 0.5时, 采取渐近式快速俯冲的软包围策略进行位置更新:

$ {\mathit{\boldsymbol{X}}(\tau + 1) = \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{Y}},f(\mathit{\boldsymbol{Y}}) < f(\mathit{\boldsymbol{X}}(\tau ))}\\ {\mathit{\boldsymbol{Z}},f(\mathit{\boldsymbol{Z}}) < f(\mathit{\boldsymbol{X}}(\tau ))} \end{array}} \right.} $

(17)

$ {\mathit{\boldsymbol{Y}} = {\mathit{\boldsymbol{X}}_{{\rm{rabbit}}}}(\tau ) - E|J{\mathit{\boldsymbol{X}}_{{\rm{rabbit}}}}(\tau ) - \mathit{\boldsymbol{X}}(\tau )|} $

(18)

$ {\mathit{\boldsymbol{Z}} = \mathit{\boldsymbol{Y}} + \mathit{\boldsymbol{S}} \times {\rm{LF}}(2)} $

(19)

其中, f(·)为适应度函数, S 为二维随机向量, 其中元素为[0, 1]之间的随机数, LF(·)是莱维飞行的数学表达式。

当| E | < 0.5且r < 0.5时, 采取渐近式快速俯冲的硬包围策略进行位置更新:

$ {\mathit{\boldsymbol{X}}(\tau + 1) = \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{Y}},f(\mathit{\boldsymbol{Y}}) < f(\mathit{\boldsymbol{X}}(\tau ))}\\ {\mathit{\boldsymbol{Z}},f(\mathit{\boldsymbol{Z}}) < f(\mathit{\boldsymbol{X}}(\tau ))} \end{array}} \right.} $

(20)

$ {\mathit{\boldsymbol{Y}} = {\mathit{\boldsymbol{X}}_{{\rm{rabbit}}}}(\tau ) - E|J{\mathit{\boldsymbol{X}}_{{\rm{rabbit}}}}(\tau ) - {\mathit{\boldsymbol{X}}_{\rm{m}}}(\tau )|} $

(21)

$ {\mathit{\boldsymbol{Z}} = \mathit{\boldsymbol{Y}} + \mathit{\boldsymbol{S}} \times {\rm{LF}}(2)} $

(22)

2.4 基于IHHO算法的TDOA定位

使用IHHO算法进行TDOA定位的具体步骤如下:

步骤1  种群初始化。根据搜索空间每一维的上界和下界, 利用式(11)初始化每个个体, 然后利用Chan算法计算一个初始解随机代替种群中的一个个体。

步骤2  计算初始适应度。利用式(9)计算所有个体的适应度值, 将适应度最优的个体位置设为当前猎物位置。

步骤3  位置更新。通过式(14)更新猎物逃逸能量, 然后根据逃逸能量和生成的随机数执行搜索或开发行为中对应的位置更新策略。

步骤4  计算适应度。通过式(9)计算位置更新后的个体适应度, 并与猎物适应度值进行比较, 若位置更新后的个体适应度值优于猎物, 则以适应度值更优的个体位置作为新的猎物位置。

重复步骤3和步骤4, 当算法迭代次数达到最大迭代次数时, 输出当前猎物位置作为目标的估计位置。

3 仿真与结果分析

本文在Matlab 2018b环境下对IHHO算法进行定位性能仿真测试, 并与带有变异机制的动态哈里斯鹰优化算法(Dynamic Harris Hawks Optimization algorithm with Mutation Mechanism, DHHO/M)^[17]、增强鲸鱼优化算法(Enhanced Whale Optimization Algorithm, EWOA)^[19]、通过竞争选择改进的蚁狮优化算法(Improved AntLion Optimizer algorithm via Tournament Selection, IALOT)^[20]以及混沌樽海鞘群算法(Chaotic Salp Swarm Algorithm, CSSA)^[21]进行比较。

仿真场景和参数设置如下:在20 m×20 m的范围内部署8个基站, 基站坐标分别为(0, 0)、(0, 10)、(0, 20)、(10, 20)、(20, 20)、(20, 10)、(20, 0)、(10, 0), 1 000个测试点均匀分布在定位场景内, 搜索空间上界 u^b=[20, 20], 下界 l^b=[0, 0], 最大迭代次数T=60。选取均方根误差(Root Mean Square Error, RMSE)作为定位精度评价指标:

$ {\rm{RMSE}} = \sqrt {\frac{1}{P}\sum\limits_{p = 1}^P {{{({x_p} - {{\hat x}_p})}^2}} + {{({y_p} - {{\hat y}_p})}^2}} $

(23)

其中, (x_p, y_p)和$({{\hat x}_p}, {{\hat y}_p})$分别为第p个测试点的真实位置和估计位置, P为测试点总数。通常根据均方根误差是否接近克拉美罗下界(Cramer-Rao Lower Bound, CRLB)来评估定位算法的性能。克拉美罗下界规定了无偏估计量所能达到的标准差下限, 代表理论最佳估计精度, 计算公式为:

$ {\rm{CRLB}} = \sqrt { {\rm{trace}} [{{({\mathit{\boldsymbol{G}}^{\rm{T}}}{\mathit{\boldsymbol{Q}}^{ - 1}}\mathit{\boldsymbol{G}})}^{ - 1}}]} $

(24)

$ \mathit{\boldsymbol{G}} = \left[ {\begin{array}{*{20}{c}} {\frac{{x - {x_2}}}{{{r_2}}} - \frac{{x - {x_1}}}{{{r_1}}}}&{\frac{{y - {y_2}}}{{{r_2}}} - \frac{{y - {y_1}}}{{{r_1}}}}\\ {\frac{{x - {x_3}}}{{{r_3}}} - \frac{{x - {x_1}}}{{{r_1}}}}&{\frac{{y - {y_3}}}{{{r_3}}} - \frac{{y - {y_1}}}{{{r_1}}}}\\ \vdots & \vdots \\ {\frac{{x - {x_N}}}{{{r_N}}} - \frac{{x - {x_1}}}{{{r_1}}}}&{\frac{{y - {y_N}}}{{{r_N}}} - \frac{{y - {y_1}}}{{{r_1}}}} \end{array}} \right] $

(25)

$ \mathit{\boldsymbol{Q}} = {\left[ {\begin{array}{*{20}{c}} {2\sigma _{\rm{r}}^2}&{\sigma _{\rm{r}}^2}& \cdots &{\sigma _{\rm{r}}^2}\\ {\sigma _{\rm{r}}^2}&{2\sigma _{\rm{r}}^2}& \cdots &{\sigma _{\rm{r}}^2}\\ \vdots & \vdots &{}& \vdots \\ {\sigma _{\rm{r}}^2}&{\sigma _{\rm{r}}^2}& \cdots &{2\sigma _{\rm{r}}^2} \end{array}} \right]_{(N - 1) \times (N - 1)}} $

(26)

其中, trace(·)函数用于求矩阵的迹, σ_r为距离噪声标准差, (x, y)和(x_i, y_i)分别为标签和基站i的坐标, r_i为标签到基站i的距离, i=1, 2, …, N, N为基站数量。

3.1 定位精度分析

将种群规模M设为30, 基站数量N设为8, 考虑实际室内定位系统的测距精度, 将距离噪声标准差σ_r的值设定在0.1 m~1.0 m范围内进行实验, 得到不同距离噪声标准差下的定位精度, 如表 1所示。可以看出, 随着距离噪声标准差的增大, 不同算法的RMSE和CRLB都逐渐增大, 而本文提出的IHHO算法在不同的距离噪声标准差环境下都取得了更低的定位误差, RMSE相比其他算法更接近CRLB, 具有较高的定位精度。

在距离噪声标准差σ_r为0.5 m的条件下, 分析基站数量对定位效果的影响, 分别在4个~8个基站参与定位的情况下进行实验, 得到不同基站数量的定位精度, 如表 2所示。可以看出, 不同算法的RMSE和CRLB随着基站数量增加而减小, 在不同基站数量情况下, IHHO算法都达到了较低的RMSE, 相比其他算法更接近CRLB。当基站数量由4个逐渐增加到8个时, IHHO算法的RMSE相比DHHO/M算法分别减小了2.79 %、10.07 %、10.56 %、13.92 %和16.87 %, RMSE减小的比例随着基站数量增加而增大, 与EWOA、IALOT和CSSA算法对比也可以得出同样的结论, 这说明IHHO算法可以更充分地利用测量值的冗余来减小定位误差, 在基站数量较多的高精度定位场景中具有更好的定位性能。

3.2 收敛性能分析

智能优化算法需要一定的种群规模和迭代次数来完成全局最优解的搜索。收敛能力好的智能优化算法可以在较小种群规模或较少迭代次数的条件下完成收敛, 降低计算量。在N=8、σ_r=0.5 m的条件下, 改变种群规模M进行实验, 得到RMSE与种群规模的关系, 如图 2所示。可以看出, EWOA算法的RMSE在种群规模达到25以上时不再显著下降, DHHO/M和IALOT算法的RMSE在种群规模达到15以上时趋于稳定, CSSA和IHHO算法可以在M=5的条件下完成位置解算, 定位精度对种群规模不敏感。IHHO算法由于引入了初始解, 使种群个体可以迅速集中到全局最优解附近进行搜索, 因此提高了寻优效率。

将种群规模M设为30, 在迭代次数不断增加的情况下, 比较不同算法适应度值和RMSE的收敛情况, 各算法适应度值和RMSE与迭代次数的关系如图 3和图 4所示。

可以看出, 适应度值和RMSE都随着迭代次数的增加呈降低趋势, 然后逐渐趋于稳定。DHHO/M、EWOA、IALOT和CSSA算法由于都采用基于最大似然估计的适应度函数, 在达到最大迭代次数后得到的定位效果比较接近。IHHO算法由于采用了改进的适应度函数, 可以对个体位置优劣做出更准确的评价, 使算法可以收敛到更小的适应度值, 并计算出RMSE更小的定位结果。CSSA算法在每次迭代中都保有一部分个体不随机运动, 降低了种群总体的随机性, 导致收敛速度较慢。DHHO/M、EWOA和IALOT算法由于前期要在整个二维空间做充分的全局搜索, 因此需要较多次迭代才能完成收敛。IHHO算法通过引入初始解简化全局搜索的步骤, 使适应度值和RMSE有一个较小的初始值, 可以通过较少的迭代次数计算出目标的位置, 克服了智能优化算法需要多次迭代的缺点。

4 结束语

本文设计一种改进的哈里斯鹰优化算法用于TDOA定位的非线性问题求解, 并对初始种群和适应度函数进行优化。与其他智能优化算法的定位性能比较结果表明, 改进的哈里斯鹰优化算法具有定位精度高、收敛速度快的优势。下一步将把该算法应用于到达时差-到达角度联合定位, 进一步提升定位精度。