周牧，龙玥辛，蒲巧林，王勇，何维

（1.重庆邮电大学通信与信息工程学院，重庆 400065；2.重庆邮电大学移动通信技术重庆市重点实验室，重庆 400065）

1 引言

随着移动互联网时代向物联网时代的逐步发展，基于位置的服务（LBS,location-based service）得到了越来越广泛的应用。全球定位系统（GPS,global positioning system）[1]和蜂窝定位系统[2]是目前最成熟的2 种室外定位系统，它们可为室外用户提供准确的位置信息，但由于室内环境的复杂性以及人员走动对信号传播的影响，用户在室内难以稳定地接收来自卫星和蜂窝基站的信号。为此，众多学者展开研究并提出了射频识别定位系统、ZigBee定位系统、蓝牙定位系统、Wi-Fi 定位系统等多种室内定位系统。其中，由于Wi-Fi 网络具有部署成本低、环境适应性强、覆盖范围广等优势，Wi-Fi室内定位系统[3]得到了广泛应用。

Wi-Fi 室内定位系统通常使用2 种测量特征：接收信号强度（RSS,received signal strength）和信道状态信息（CSI,channel state information）。RSS指纹特征维数低且空间分辨率低，导致定位精度低。相比于RSS，CSI 包含了信号传输过程中更细粒度和多样化的物理层信息（如多通道子载波相位和幅度信息），可用于刻画信号散射、衰落和功率随距离衰减的综合效应，从而被广泛应用于Wi-Fi高精度定位系统。近年来，国内外学者提出了基于几何、基于指纹以及基于数据融合等的多种CSI 室内定位技术，其中，基于CSI 指纹的室内定位技术因其处理效率高、抗干扰性强、定位精度高等优点得到了广泛应用，但仍然面临诸多挑战。

首先，由于手机采集的CSI 相较于RSS 来说，更容易受到环境的影响，稳定性较差，因此需要进行降噪处理。早期的局部平滑滤波算法，如高斯滤波[4]和均值滤波[5]，通过平滑处理来去除噪声，但这类方法难以处理图像的纹理信息等非平滑部分，所以不能在降噪的同时有效保留纹理信息。空域降噪算法[6]以随机噪声的零和特点为理论基础进行降噪处理，但其主要针对随机噪声的降噪，而未去除因设备缺陷导致的一些脉冲噪声。考虑到矩阵的表示形式容易破坏数据的原始空间结构，文献[7]提出了一种基于非局部相似的非负Tucker 分解方法，实现了空间域的非局部相似性和光谱域的全局相似性。该方法利用非局部相似性，在合适的窗口大小下，将高光谱图像中提取的三维全波段斑块进行分组并形成三阶张量。实验结果表明，该方法在减少计算量的同时提高了高光谱图像的质量。

其次，降噪后的CSI 中稀疏分布着有用的结构和成分，为了提取其隐含的特征，需要进行特征提取。文献[8-9]提出一种主体组合成分分析和线性决策分析的方法来提取最相关的特征向量以减少维数。文献[10]提出了一种基于二维离散小波变换的跨域尺度变换方法，在对序列图像进行二维小波变换获得主要特征点的同时，实现了数据的降维和噪声抑制；构造了一种新的归一化相关代价函数，并利用黄金分割算法高效地找到最优旋转角速度。但由于二维小波变换仅着眼于数据的频域相关特性，忽略了时域−频域之间的潜在关系和重要的判别信息，因此其无法适用于三维CSI 的高低频信号变化信息的特征提取。

另外，基于CSI 指纹的定位技术在室内定位领域受到了广泛关注，如何利用CSI 实现高精度的室内定位成为众多学者研究的方向。文献[11]提出了一种基于CSI 幅值指纹的窄带物联网定位算法，其采用多维标度分析方法计算目标点与参考点之间的欧氏距离和时间反转共振强度，然后采用K 近邻（KNN,K nearest neighbor）算法进行位置估计。多维标度分析方法通过直观的空间图再现研究对象之间的关系，它的缺陷在于认为各维度对目标的贡献相同。文献[12]提出了一种基于位置不确定性约束的机器人室内定位方法，并利用贪婪算法来求解优化问题。但贪婪算法依赖于当前已经做出的选择，所以其无法保证得到最优解。文献[13]提出利用矢量格式的CSI 数据作为位置指纹，再将三层半连接神经网络和一层全连接神经网络相组合用于位置估计。然而，当面对较复杂的非线性问题时，神经网络容易产生收敛速度慢、网络不稳定和陷入局部最优等一系列问题。

针对上述问题，本文提出了一种基于CSI 张量分解的室内Wi-Fi 指纹定位方法，方法流程如图1 所示。

2 CSI 预处理

鉴于三阶张量可以描述CSI 的信息和结构，本文采用张量的形式来表达复杂的CSI。其目的在于保留数据原有形式的同时，最大限度地保留图像内在的结构信息，通过建立一种以多维度分析为主线的张量处理框架，实现张量的降噪处理和特征提取，提升数据处理分析的能力。

2.1 CSI 降噪处理

一般地，智能手机采集到的原始CSI 可表示为一幅X轴为子载波、Y轴为数据包且Z轴为CSI 幅值的三维图像[14]，如图2 所示，子载波、数据包、CSI 幅值的数量分别为L1、L2、L3。由于子载波的CSI 幅值在不同时刻的波动较大且容易受到环境的影响，导致原始CSI 的三维图像中出现了较多的奇异值，进而影响指纹库的构建甚至定位的性能，因此需对原始CSI 进行降噪处理。

本文将原始CSI 三维图像表示为一个二元三阶张量（其元素仅包含0 和1）。具体而言，令张量的形状为[L1,L2,L3]，即第i(i=1,…,3)维（或称“阶”）有Li个元素，将该张量记为O∈RL1×L2×L3，并通过秩一张量的线性组合将其表示为

其中，μr(∈RL1)、υr(∈RL2)和ωr(∈RL3)分别为第r(r=1,…,M)个秩一张量在含噪张量O的3 个维度上分解得到的单位向量；符号°表示向量的外积运算；λ为第rr个秩一张量的组分奇异值，其刻画了第r个秩一张量的组分在整体中的比重；M为用于重构O的秩一张量的个数。

为了利用CSI 的时域−频域互补信息，在保持幅值连续性的同时实现高保真降噪的目的，本文提出基于平行因子（PARAFAC,parallel factor）分析模型的张量分解算法（如图3 所示）来估计用于重构O中无噪张量S（即S=O−N）的秩一张量的个数（即分解级数）k，其中，N为张量O中的噪声。由于该算法未对各维度进行平滑处理，因此没有引入额外的信息，保持了图像时域−频域的一致性，其相较于其他张量分解算法，可以取得更好的降噪效果以提升图像质量。此外，k的表达式为

其中，SNR 为张量O的信噪比（即信号和噪声的强度比），Li和Ki分别为张量S中第i维度的维度数和张量秩。

由于SNR 等价于信号和噪声的方差比，因此本文对含噪张量O分块进行信号和噪声的方差估计。先将张量O分成A个体积相等的子张量，再将每个子张量分为B个体积相等的局部张量。将第f(f=1,…,A)个子张量的方差表示为，其中，为第f个子张量中第e个局部张量的灰度值，为第f个子张量的平均灰度值，则噪声方差的估计值可表示为

其中，为第f个子张量的方差，为A个子张量的方差的平均值。此时，信号方差的估计值可表示为

于是，可得

为了得到无噪张量S在第i维度上的张量秩Ki，本文根据赤池信息准则（AIC,Akaike information criterion）[15-16]进行估计，其估计式（具体推导过程见附录1）为

张量S在3 个维度上的展开过程如图4 所示，为了使重构的无噪张量S尽可能逼近理想无噪张量，本文令两者之间的均方根误差最小，即。同时，考虑到张量模型中各维度的秩一张量可在求解中合并为因子矩阵，本文采用交替最小二乘（ALS,alternate least squares）迭代算法[17]来求解上述最优化问题。为此，令迭代次数为t，第t次迭代得到的各维度因子矩阵分别为和，且第t次迭代得到的权矩阵为。ALS 迭代算法的基本思路为固定2 个矩阵来求解剩余矩阵，本文选择固定矩阵V和W来求解矩阵U。由于3 个因子矩阵均未知，现初始化因子矩阵为全1 矩阵，权矩阵Λ0为单位矩阵。将第t−1 次迭代结果的第 1 维度展开矩阵表示为，则第t次迭代得到的加权因子矩阵分别为

利用更新后的因子矩阵可得重构第t次迭代结果为

2.2 CSI 特征提取

降噪后的张量中稀疏分布着有用的结构和成分，为了找到其隐含的结构和成分，需对张量进行特征提取，以找到具有一定物理意义的数据表示。然而，现有的特征提取算法主要着眼于二维数据的频域相关特性，忽略了时域−频域存在的潜在关系，损失了重要的判别信息。因此，本文引入张量小波分解算法来获取子载波、数据包以及CSI 幅值这3 个维度上的高低频信号变化信息来实现CSI 特征提取。张量分解后得到的高低频分量的关系为

其中，符号⊕和⊗分别表示张量的直和运算和克罗内克积运算，L和H分别表示作用于X、Y和Z这3 个维度上的离散小波低通和带通滤波器。此时，将得到的八组小波子成分分别记为LLL、LLH、LHL、LHH、HLL、HLH、HHL和HHH。

其中，Pm(i,j,k)为小波子成分中位于(i,j,k)处的元素值。由于张量小波子成分所包含的多方向和多频率信息已被文献[18]证明能够提供类别判别信息，因此本文将其应用于构造CSI 位置指纹。具体而言，构造第n(n=1,…,Nf)个参考点处的CSI 位置指纹为，其中，Nf为参考点的个数，为第n个参考点处第m组小波子成分的小波系数。CSI 降噪处理及特征提取过程如算法1 所示。

算法1CSI 预处理算法

3 基于PLSR 的建模与定位

CSI 指纹定位方法包括2 个阶段：离线阶段和在线阶段。离线阶段设置参考点并进行CSI 的采集和预处理以得到标签样本（含参考点的CSI 位置指纹和位置坐标），然后建立定位模型；在线阶段设置测试点并进行CSI 的采集和预处理以得到待定位样本（含测试点的CSI 位置指纹），然后将其代入定位模型以得到测试点的位置坐标估计。

为了研究CSI 位置指纹和位置坐标之间多重相关变量的相互依赖关系，本文以参考点的CSI 位置指纹为自变量且位置坐标(x n,yn)为因变量，利用偏最小二乘回归（PLSR,partial least squares regression）算法进行多因变量对多自变量的线性回归建模。PLSR 算法不仅可以解决典型相关分析中无法获得自变量到因变量直接映射的问题，又可避免多元线性回归分析中因为自变量之间存在相关性而导致的过拟合问题。

为了求解w1和c1，引入拉格朗日乘子L为

分别对w1和c1求偏导并令其为0，得到

经推导，可得w1和c1分别为对称矩阵XTYYTX和YTXXTY的最大特征值所对应的单位特征向量。将求解后的w1和c1分别代入t1=Xw1和u1=Yc1，可得X和Y的第一对主成分t1和u1。

为了解决从X到Y的映射问题，根据主成分回归思想将X和Y分别对它们的主成分t1和u1进行回归建模，得到，其中，E1和G1为残差矩阵。将X中主成分t1不能解释的残差矩阵E1作为新的X，且Y中主成分u1不能解释的残差矩阵G1作为新的Y，不断提取新的主成分，直到主成分数量达到上限（即X的秩），算法结束。令算法结束时共得到a对主成分，则可将初始X和Y分别表示为

4 实验结果

4.1 实验场景

本文实验分别在走廊（即场景1）和实验室（即场景2）中进行，如图5 所示。场景1 和场景2 的面积分别为42.08 m×3.12 m和22.72 m×8.04 m。图5 中，圆点表示参考点（RP,reference point），三角形表示测试点（TP,test point）。场景1 和场景2中各设置了40 个RP，并随机设置了6 个和8 个TP。本文实验选择商用D-Link 作为无线接入点（AP,access point）以及装有Nexmon 测试平台的Google Nexus 6 智能手机作为接收端。在每个位置处，接收端采集300 个包含64 个子载波（带宽为20 MHz）的CSI 数据包，并将采集到的CSI 转换为一个二元三阶张量以进行降噪处理和特征提取。

4.2 降噪结果分析

在不同场景中采集到的CSI 受到不同程度的噪声干扰，导致某些子载波在不同时刻采集到的CSI幅值出现较大偏差，且由于手机采集到的CSI 本身具有波动幅度大、稳定性差等问题，因此本文实验需对原始CSI 数据进行子载波筛选。原始CSI 分布如图6(a)和图6(b)所示，图6 中“+”表示在采样过程中CSI 幅值出现极值的情况，“箱线”表示各子载波的CSI 幅值分布，其从上到下依次表示CSI幅值的最大值、上四分位值、中位值、下四分位值和最小值。由于场景1 相对空旷，人员流动少，因此相比于场景2，在场景1 中采集到的数据受到噪声的干扰更小，极值也更少。为了在降低计算开销的同时剔除幅值偏差较大的子载波，本文筛选出采样过程中幅值方差的均值最小的30 个最优子载波，如图6(c)和图6(d)所示。由图6 可知，经过筛选后各子载波的CSI 幅值分布更加集中，极值也更少。

本文引入小波降噪算法作为对比算法，选择的小波基及其设置的分解层数N分别为：Haar（N=5）、Coiflet3（N=1）、Symlet8（N=8）。另外，在实验过程中，这3 种算法阈值选取都采用固定阈值估计，阈值函数都选用硬阈值，并规定了阈值处理不随噪声水平变化。将筛选后的CSI 分别采用本文提出的基于PARAFAC 分析模型的张量分解算法和传统的小波降噪算法（如Haar 小波降噪、Coiflets 小波降噪和Symlets 小波降噪）进行降噪处理，得到图7 和图8所示的降噪效果。可以看出，经本文算法处理得到的图像更加平滑，CSI 幅值分布更加集中，偏差也更小，由此说明本文算法的降噪效果优于其他3 种算法。为了更准确地比较和分析不同算法的降噪性能，本文采用SNR、均方误差（MSE,mean square error）以及峰值信噪比（PSNR,peak signal noise ratio）3 个参数对2 种场景下采用不同算法进行降噪处理的性能进行评估，如图9 所示。

由图9 可知，分别在场景1 和场景2 中采集到的数据经本文算法处理后的CSI 图像的SNR 分别为64.7 dB 和61.8 dB，相比于原始CSI 图像的SNR，提升了15 dB 以上；MSE 分别为6.21 和6.20，较传统小波降噪算法，经本文算法处理后的图像更接近于理想无噪图像；PSNR分别为40.20 dB 和40.21 dB，进一步证实了本文算法具有更好的降噪性能。分别对比同一场景下不同算法的性能参数可知，针对场景1中采集到的数据，选用本文算法或Symlets 小波降噪算法[20]进行降噪处理的性能更优；而针对场景2中采集到的数据，选用本文算法或Haar 小波降噪算法[21]进行降噪处理的性能更优。由此进一步验证了本文算法在2 种场景下都具有更优的降噪性能。

4.3 特征提取结果分析

图10 对比了各小波子成分在不同数据包个数下所对应的ASM 值分布，可见经张量小波分解得到的不同位置处同一小波子成分的ASM 值分布是相对集中的，且由于小波子成分HHH作为细节分量包含了分解前CSI 的绝大部分信息，因此当选取不同个数的数据包时，小波子成分HHH的ASM 值为8 个小波子成分中最大的，并随着数据包个数的增加，相同数据包个数下不同小波子成分ASM 值的差值增大。

4.4 定位性能分析

4.4.1 定位误差分析

在保留本文所提其余算法的情况下，本节依次更换了降噪处理、特征提取和定位阶段的算法来分析本文中3 种主要算法（即PARAFAC、张量小波分解和PLSR）对定位性能的影响，如图11 所示。

基于PARAFAC分析模型的张量分解算法在去除噪声分量的同时，最大限度地保留了数据内在的结构信息。图11(a)比较了在RP 取不同个数的情况下，分别经本文算法降噪处理和无降噪处理时的定位误差累积分布函数（CDF,cumulative distribution function）。当RP 个数相同时，经本文算法进行降噪处理后的定位性能明显优于无降噪处理时的定位性能，且随着RP 数的增加，算法的定位性能也越好。归其原因，在于当RP 数越多时，经降噪处理和特征提取后得到的标签样本也越多，于是基于更多的样本进行回归建模所得到的模型将更接近于真实模型，进而最终的定位结果也会更加准确。图11(b)比较了选用不同降噪算法时的定位CDF。观察可知，本文算法在定位误差为4 m 时的置信概率（即CDF 的值）为94.88%，高于其他3 种小波降噪算法的置信概率，分别为91.52%、85.88%和80.21%。

由图11(c)可知，相较于方向梯度直方图（HOG,histogram of oriented gradient）特征算法[22]（定位误差为4 m 时置信概率为84.84%）和局部二值模式（LBP,local binary pattern）特征算法[23]（定位误差为4 m 时置信概率为83.08%），采用本文算法进行特征提取将对定位性能的提升提供更大帮助。此外，由图11(d)可知，本文算法在定位误差为4 m 时的置信概率（即94.88%）优于加权K 近邻（WKNN,weighted K nearest neighbor）算法[24]的 84.9%（Np=3）、88.54%（Np=5）、82.72%（Np=7）和79.97%（Np=9），近似最佳三角形内点测试（APIT,approximate perfect point-in-triangulation test）算法[25]的81.97%以及贝叶斯（Bayes）算法[26]的73.14%，其中，Np 表示邻近点个数。

4.4.2 时间开销分析

除了定位误差以外，时间开销也是评估定位性能的重要指标。不同RP 数下4 种定位算法的运行时间如表1 所示，用于定位的RP 数越多，对应的运行时间也越长。但无论有多少个RP，本文所提PLSR 算法的运行时间都是最短的。以RP=40 为例，此时WKNN、APIT 和Bayes 算法的运行时间分别为0.83 s、0.75 s 和0.74 s，均大于本文算法的运行时间0.67 s，这表明本文算法在计算复杂度方面比WKNN、APIT 和Bayes 算法更有优势。

表1 不同RP 数下4 种定位算法的运行时间

5 结束语

本文提出采用张量的形式来表达CSI 的优势在于，保留数据原有存在形式的同时最大限度地保留了图像内在的结构信息。本文首先研究了基于PARAFAC 分析模型的张量分解算法和ALS 迭代算法相结合用于降噪处理的可行性；然后，利用张量小波分解算法在CSI 的3 个维度上进行小波分解实现特征提取，有效降低了CSI 维数并得到CSI 位置指纹；最后，基于PLSR 算法建立定位模型，对位置坐标进行预测，进而实现定位。实验结果表明，本文算法在定位误差3.5 m、4 m 和4.5 m 内的置信概率分别为89.81%、94.88%和98.05%，均明显优于其他现有算法，由此验证了本文提出的基于CSI张量分解的室内Wi-Fi 指纹定位方法在提升数据处理分析能力和拟合CSI 位置指纹和位置坐标关系的同时，还具有更优的定位性能。

基于Tucker分析模型的CSI降噪处理是一个值得研究的问题，下一步工作将对该问题进行深入探讨；此外，CSI 图像的灰度共生矩阵是进行特征提取的有力工具，但它不能直接提供进行类别判断的特性，所以作者还将在灰度共生矩阵的基础上研究用于定量描述图像特征的统计属性，并由这些统计属性来构造用于室内Wi-Fi 指纹定位的CSI位置指纹。

附录1 式(6)推导过程

将式(6)改写为

为了求解张量秩K，需从可供选择的模型中选择AIC值最小的模型。AIC 的构造如下

假设噪声是与信号无关的均值为0 的高斯随机过程，则噪声的协方差矩阵可表示为σ2I。为了确定r(r≤min (L,N))的值，构造模型x=As+n，其中，矩阵A∈RL×r，向量s∈Rr×1。于是，可得模型的协方差矩阵R=Ψ+σ2I，矩阵Ψ=ASAH，信号的协方差矩阵S=E(ssH)。假设矩阵A为列满秩，则A(Φj)线性无关，并假设矩阵S为满秩，则Ψ秩为r，那么Ψ的L−r个最小特征值等于0。矩阵R∈RL×L的各特征值之间的关系为λ1≥λ2≥…≥λL，其中，最小的L−r个特征值等于σ2，即λr+1=λr+2=…=λ L=σ2，前r个特征值对应的特征向量表示为Vj。

r的值可由矩阵R的最小特征值的多重性（即个数）推测得到，但实际上协方差矩阵R是未知的，所以基于线性代数中的谱表示定理，将R( r)=Ψ(r)+σ2I中的R(r)表示为，模型的参数向量表示为。利用最大似然估计算法来计算已知样本X={x1,…,xN}所对应的参数向量。因此，计算样本的联合概率密度为

步骤 2求解Num。由于Θ(r)中参数个数为(r+1 +2Lr)，则此时的自由度为(r+2Lr)，又由于Θ(r)中自由参数的个数等于由Θ(r)所张成的空间的自由度，因此需对特征向量进行正交化和标准化处理，可得

观察可知，式(31)中的未知量仅为r。于是，可通过调整r∈{0,…,L−1}的取值来改变AIC 的值，AIC 取得最小值时所对应的r即张量秩K。