基于无线传感器网络(wireless sensor network, WSN)的目标定位技术因其广阔的应用前景，近年来广受关注和研究^[1-2]。大多数关于目标定位的研究都集中在有源定位^[3-4]上，即目标需携带无线设备。与有源定位不同，无源定位^[5-6]无需为目标配备任何无线收发设备，其通过分析目标对无线链路的阴影效应估计目标位置，可应用于入侵检测、移动医疗、边界保护、安防监控和智能感知等诸多场景中^[7]。基于无线传感器网络的无源定位技术主要可分为4类：1)基于几何的无源定位技术^[8]；2)基于指纹的无源定位技术^[9]；3)基于无线层析成像(radio tomographic imaging, RTI)的无源定位技术^[10]；4)基于压缩感知(compressive sensing, CS)的无源定位技术^[11]。其中，基于几何的方法通过寻找被遮挡无线链路的重叠覆盖区域来估计目标位置。该方法定位精度较低，且对目标间距有限制。为提高定位精度，在无源定位中引入指纹识别的方法。基于指纹的方法根据现有的指纹库和当前测量值估计目标位置，能有效提升定位性能。但是，指纹的采集和指纹库的建立需要进行费时费力的现场调查。基于无线层析成像的方法采用计算机断层扫描技术，可以在不了解目标数量和指纹数据库的前提下定位许多目标，且定位精度较高。但是该方法需要密集部署无线链路，这将导致较高的硬件成本和能源问题。为减少测量次数，压缩感知理论被应用于无源定位中。该理论将数据采样和压缩过程相结合，在采集信号的同时完成数据压缩，有效避免了采样后再压缩导致的资源浪费。该理论表明，通过对高维的稀疏信号进行欠采样，可得到低维的测量数据，再利用稀疏重构算法可以较大概率估计出原信号。基于压缩感知的方法将目标位置估计转化为位置向量的稀疏重构，所需部署的无线节点数量大大降低，且能在测量值较少的情况下，有效利用目标位置分布的空间稀疏性，获得较高的定位精度。

在基于压缩感知的无源定位中，无源字典的建立直接影响到系统的定位性能^[12]。其建立方法主要分为指纹训练和理论模型。由于采集指纹需要耗费大量人力和物力资源，所以通常采用基于理论模型的方法建立无源字典。传统的理论模型包括矩形模型^[13]、椭圆模型^[14]及衍射模型^[15]。但这些模型都存在拟合目标阴影效应不好的问题，进而导致定位精度不高。文献[16]提出更加精确的鞍面模型。鞍面模型在计算阴影效应的过程中考虑了目标位置与链路的垂直距离，更贴合实际定位环境中的情况，有效提高了定位精度。

传统的压缩感知无源定位方法都是基于接收信号强度(received signal strength, RSS)的^[17]。上述模型均使用RSS变化量拟合目标阴影效应。但RSS作为粗粒度指标，本身定位精度不高，且RSS测量值会受到反射、绕射与散射等现象的影响。相位是无线信号的一个重要特性，比信号强度更具细粒度，基于信号相位的定位方法理论上具有更高的定位精度。文献[18]在鞍面模型的基础上提出基于相位信息的理论模型，并首次将其运用于基于无线层析成像的无源定位中，取得了优于RSS模型的定位结果。

为有效利用接收信号的相位信息，本文在压缩感知框架下设计基于相位偏移(phase shift, PS)的无源定位方法。该方法将接收信号相位偏移值作为观测数据，结合变分贝叶斯推理，恢复目标位置稀疏向量，估计目标位置。同时进行网格裁剪，降低运算复杂度，减少计算时间。该方法克服了传统基于RSS的压缩感知无源定位方法定位精度不高、在变化环境中定位误差增大的问题，有效提升了系统定位性能。

1 系统模型

无源多目标定位场景如图 1所示。定位区域内随机分布有K个目标，为估计目标位置，在定位区域四周部署能够采集相位信息的传感器节点，测量各无线链路的接收信号相位，从而计算相位偏移值。为更好地覆盖定位区域，将传感器均匀部署在定位区域的边界上。由于压缩感知理论针对的是离散信号，为了将压缩感知中的稀疏恢复算法用到无源定位中，需对定位区域进行网格化处理。将定位区域划分为N个大小相同的网格，并按顺序编号，即1, 2, 3, …, n, …, N。用一个N维向量θ表示K个目标的位置分布

$ \begin{equation*} \boldsymbol{\theta}=\left[\theta_{1}, \theta_{2}, \cdots, \theta_{n}, \cdots, \theta_{N}\right]^{\mathrm{T}} . \end{equation*} $

(1)

其中：θ∈R^N×1为目标位置向量。当第n个格点上存在目标时，θ_n=1；否则θ_n=0。定位区域中目标的个数可表示为K=‖θ‖₀。此外，K≪N，θ为K稀疏向量。

设传感器数量为2M，每2个传感器(一发一收)组成一条无线链路，则共有M条无线链路。当目标位于定位区域中时，会对某些无线链路的信号传播路径产生遮挡，这时无线链路的接收信号相位会发生改变。而且，一个目标在不同的位置会对周围的无线信号产生不同的影响。因此，可根据目标对无线链路的阴影效应，估计目标位置。这里，使用接收信号的相位偏移值拟合目标阴影效应。将第m条无线链路的接收信号相位表示为φ_m，则其相位偏移值Δφ_m可计算为

$ \begin{equation*} \Delta \varphi_{m}=\left|\varphi_{m}-\varphi_{m}^{0}\right| \text {. } \end{equation*} $

(2)

其中：φ_m⁰为第m条无线链路空闲时的接收信号相位(m≤M)。

根据PRS权重模型^[18]，无线链路的空间影响区域如图 2所示。图中阴影区域即为无线链路的空间影响区域，只有在空间影响区域内的目标才会对该无线链路的接收信号相位产生影响。

目标阴影效应所导致的相位偏移可计算为

$ \begin{align*} \phi_{l, i}= & A_{l}\left(\frac{1-c_{l}}{\left(0.5 d_{l}\right)^{2}} U_{i}^{2}+c_{l}\left(1-\frac{V_{i}^{2}}{R_{l}^{2}}\right)\right) . \\ & \text { s.t. } \frac{U_{i}^{2}}{\left(0.5 d_{l}\right)^{2}}+\frac{V_{i}^{2}}{R_{l}^{2}} \leqslant 1 . \end{align*} $

(3)

其中：ϕ_{l, i}表示位于第i个网格的目标对链路l的阴影效应所导致的相位偏移，(U_i, V_i) 为第i个网格在U-V坐标系中的坐标，d_l为链路l的长度，R_l为空间影响区域的短轴长，A_l为链路l相位偏移的最大值，c_l为链路l中点的归一化相位偏移(A_l和c_l均为与环境相关的参数)。根据式(3)，可构建无源字典Φ

$ \boldsymbol{\varPhi}=\left[\begin{array}{cccc} \phi_{1, 1} & \phi_{1, 2} & \cdots & \phi_{1, N} \\ \phi_{2, 1} & \phi_{2, 2} & \cdots & \phi_{2, N} \\ \vdots & \vdots & \ddots & \vdots \\ \phi_{M, 1} & \phi_{M, 2} & \cdots & \phi_{M, N} \end{array}\right]. $

(4)

通过部署在定位区域四周的接收传感器所采集到的相位信息，根据式(2)，可以计算出接收信号的相位偏移值，并将其作为观测数据，由此可建立测量模型

$ \begin{equation*} {\boldsymbol{y}}={\boldsymbol{\varPhi}} {\boldsymbol{\theta}}+{\boldsymbol{\varepsilon}} . \end{equation*} $

(5)

其中：y∈R^M×1为观测数据，表示M条无线链路的相位偏移值，其第m个元素y_m=Δφ_m；ε表示噪声向量。

2 基于相位偏移的定位算法 2.1 分层先验模型

为诱导目标位置向量θ的稀疏性，引入一个2层高斯先验模型。在先验模型的第1层，θ被定义为隐藏变量，服从高斯分布。具体地，将θ的高斯先验分布定义为

$ \begin{align*} & p(\boldsymbol{\theta} \mid \boldsymbol{\alpha})=\prod\limits_{n=1}^{N} \mathrm{~N}\left(\theta_{n} \mid 0, \alpha_{n}^{-1}\right)= \\ & (2 {\rm{ \mathsf{ π}}})^{-\frac{N}{2}}|\boldsymbol{\varLambda}|^{\frac{1}{2}} \exp \left(-\frac{1}{2} \boldsymbol{\theta}^{\mathrm{T}} \boldsymbol{\varLambda} \boldsymbol{\theta}\right) . \end{align*} $

(6)

其中: α_n为θ_n的逆方差，α=[α₁, α₂, …, α_n, …α_N]^T，Λ=diag(α)。

在先验模型的第2层，α_n被视为隐藏变量，假设α的先验分布为Gamma分布，即

$ \begin{gather*} p(\boldsymbol{\alpha} ; c, \boldsymbol{d})=\prod\limits_{n=1}^{N} \operatorname{Gamma}\left(\alpha_{n} \mid c, d_{n}\right)= \\ \prod\limits_{n=1}^{N} \frac{1}{\Gamma(c)} d_{n}^{c} \alpha_{n}^{c-1} \cdot \exp \left(-d_{n} \alpha_{n}\right) . \end{gather*} $

(7)

其中: $ \Gamma(c)=\int_{0}^{\infty} x^{c-1} \mathrm{e}^{-x} \mathrm{~d} x, c$和d={d₁, d₂, …, d_n, …, d_N} 是Gamma分布的确定性参数。

ε为零均值噪声向量，服从高斯分布，假设其逆方差为β。基于式(5)中提出的测量模型，似然函数可表示为

$ p(\boldsymbol{y} \mid \boldsymbol{\theta}, \beta)=\left(2 {\rm{ \mathsf{ π}}} \beta^{-1}\right)^{-\frac{M}{2}} \exp \left(-\frac{\beta}{2}\|\boldsymbol{y}-\boldsymbol{\varPhi} \boldsymbol{\theta}\|_{2}^{2}\right). $

(8)

将β也视为一个隐藏变量，假设其遵循一个带有确定性参数a和b的Gamma先验分布

$ \begin{equation*} p(\beta ; a, b)=\operatorname{Gamma}(\beta \mid a, b)=\frac{1}{\Gamma(a)} b^{a} \beta^{a-1} \mathrm{e}^{-b \beta}. \end{equation*} $

(9)

2.2 基于变分贝叶斯推理的稀疏重构

在上面提出的分层先验模型中，y为观测数据，z={θ,α, β} 为隐藏变量，可以利用变分贝叶斯推理得到z的后验分布。此外，Ω={a, b, c,d} 为先验模型中的确定性参数，通常被赋给非常小的值，为z提供非信息先验。根据变分贝叶斯推理，z的对数后验近似为

$ \ln q(\boldsymbol{\theta})=\langle\ln p(\boldsymbol{y}, {\boldsymbol{z}} ; \boldsymbol{\Omega})\rangle_{q(\boldsymbol{\alpha}) q(\beta)}+\text { const. } $

(10)

$ \ln q(\boldsymbol{\alpha})=\langle\ln p(\boldsymbol{y}, {\boldsymbol{z}} ; \boldsymbol{\Omega})\rangle_{q(\boldsymbol{\theta}) q(\beta)}+\text { const. } $

(11)

$ \ln q(\boldsymbol{\beta})=\langle\ln p(\boldsymbol{y}, {\boldsymbol{z}} ; \boldsymbol{\Omega})\rangle_{q(\boldsymbol{\theta}) q(\boldsymbol{\alpha})}+\text { const. } $

(12)

其中: const表示常数，用以补全后验概率; 〈·〉表示求期望, p(y,z; Ω) 表示隐藏变量z和观测数据y的联合分布。根据概率链式法则，p(y,z; Ω) 可分解为

$ p(\boldsymbol{y}, \boldsymbol{z} ; \boldsymbol{\Omega})=p(\boldsymbol{y} \mid \boldsymbol{\theta}, \beta) p(\boldsymbol{\theta} \mid \boldsymbol{\alpha}) p(\beta ; a, b) p(\boldsymbol{\alpha} ; c, \boldsymbol{d}). $

(13)

根据α，β的分布，将式(13)代入式(10)，可得到θ的后验分布为

$ \begin{equation*} q(\boldsymbol{\theta})=\mathrm{N}(\boldsymbol{\theta} \mid \boldsymbol{\mu}, \boldsymbol{\varSigma}) . \end{equation*} $

(14)

由于θ服从高斯分布，根据式(14)，其后验分布的均值向量μ和协方差矩阵Σ可表示为

$ \boldsymbol{\mu}=\langle\beta\rangle \boldsymbol{\varSigma} \boldsymbol{\varPhi}^{\mathrm{T}} \boldsymbol{y} . $

(15)

$ \boldsymbol{\varSigma}=\left(\langle\beta\rangle \boldsymbol{\varPhi}^{\mathrm{T}} \boldsymbol{\varPhi}+\langle\boldsymbol{\varLambda}\rangle\right)^{-1} . $

(16)

根据θ的分布，将式(13)代入式(11)，可得到α的后验分布为

$ \begin{equation*} q(\boldsymbol{\alpha})=\prod\limits_{n=1}^{N} \operatorname{Gamma}\left(\alpha_{n} \mid \tilde{c}, \tilde{d}_{n}\right) . \end{equation*} $

(17)

由于α服从Gamma分布，根据式(17)，其后验分布的参数可表示为

$ \tilde{c}=c+\frac{1}{2} . $

(18)

$ \tilde{d}_{n}=d_{n}+\frac{1}{2}\left(\mu_{n}\right)^{2}+\frac{1}{2} \sum\nolimits_{n, n} . $

(19)

根据θ的分布，将式(13)代入式(12)，可得到β的后验分布为

$ \begin{equation*} q(\beta)=\operatorname{Gamma}(\beta \mid \tilde{a}, \tilde{b}) . \end{equation*} $

(20)

由于β服从Gamma分布，根据式(20)，其后验分布的参数可表示为

$ \tilde{a}=a+\frac{M}{2} . $

(21)

$ \tilde{b}=b+\frac{1}{2}\left(\operatorname{tr}\left(\boldsymbol{\varPhi} \boldsymbol{\varSigma} \boldsymbol{\varPhi}^{\mathrm{T}}\right)+\|\boldsymbol{y}-\boldsymbol{\varPhi} \boldsymbol{\mu}\|_{2}^{2}\right) . $

(22)

根据上述结果，隐藏变量α_n、β的期望可分别表示为$ \left\langle\alpha_{n}\right\rangle=\tilde{c} / \tilde{d}_{n}, \langle\beta\rangle=\tilde{a} / \tilde{b}$。设计基于相位偏移的定位算法，如表 1所示。

当迭代更新隐藏变量的近似后验时，α的许多分量通常在收敛时取非常大的值。当α_n足够大时，θ_n就会变成较小的值，这种情况下，网格n对观测数据y的影响可以忽略不计。这意味着，可以将该网格从网格集Π中移除，在迭代过程中进行实时网格裁剪，以减少计算量。我们主要根据α的当前近似后验分布来裁剪网格，设定阈值α_th，则每次迭代中，网格集Π被更新为

$ \begin{equation*} \Pi^{\text {new }}=\Pi^{\text {old }}-\left\{n \mid\left\langle\alpha_{n}\right\rangle>\alpha_{\text {th }}\right\} . \end{equation*} $

(23)

其中：Π^new表示裁剪后的网格集，Π^old表示当前网格集。根据裁剪后的网格集，可以对无源字典Φ进行相应的裁剪。从而，所提定位算法的计算量将随着网格数量的减少而降低。

在每一次迭代中，测量残差γ计算为

$ \begin{equation*} \gamma=\|\boldsymbol{y}-\boldsymbol{\varPhi} \boldsymbol{\mu}\|_{2} . \end{equation*} $

(24)

设定测量残差阈值γ_th，当γ＜γ_th或迭代次数τ达到最大迭代次数τ_max时，停止判据。设定稀疏阈值μ_th，滤除均值μ中较小的非零分量。最后，根据最终得到的μ恢复位置向量θ，估计目标位置。

3 仿真结果及分析

本节将进行一系列仿真实验来验证所提算法的性能。在仿真中，目标监测区域A被设定为一个6.5 m×6.5 m的正方形区域。将其等分为N=169个边长为0.5 m的正方形网格。K=4个目标随机分布在A中，M=26条无线链路沿A的四周均匀分布。使用M条无线链路的接收信号相位偏移值作为观测数据。为检验所提算法的可靠性和鲁棒性，在每个测量值上添加高斯白噪声。测量值的信噪比定义为SNR=10lg(‖Φθ‖₂²/(Mσ²))，其中σ²表示噪声向量ε的方差。在仿真中，设定信噪比的默认值为20 dB。为验证所提的基于相位偏移的压缩感知无源定位算法(device-free localization-phase shift, DFL-PS)的有效性和鲁棒性，将其与基于接收信号强度的压缩感知无源定位算法(device-free localization-received signal strength, DFL-RSS)进行了比较。

首先，评估网格裁剪过程的有效性，并研究裁剪阈值α_th对所提算法定位性能的影响。图 3展示了不同裁剪阈值(α_th=1, 3, 10, 20, ∞，其中α_th=∞为不进行网格裁剪)下，所提算法的平均定位误差(这里，平均定位误差定义为真实目标位置和估计目标位置之间的平均欧氏距离)和平均定位时间。如图 3(a)所示，当α_th=1时，平均定位误差最大。这是因为α_th的值越小，越容易滤除目标存在的网格，而错误的裁剪将导致定位误差的大幅上升。当α_th=10、α_th=20和α_th=∞时，定位误差较小，且基本相同。因此，α_th选择合适的值进行网格裁剪，不会影响系统的定位精度。如图 3(b)所示，随着最大迭代次数τ_max的增加，不进行网格裁剪将导致平均定位时间呈线性增加，而实施网格裁剪的方案的平均定位时间只呈现小幅上升。同时，α_th的值越小，平均定位时间越短。由此可见，网格裁剪过程可有效降低运算复杂度，节省计算成本。基于以上仿真结果，选择裁剪阈值α_th=10作为定位精度和计算成本之间的权衡。

接着，模拟一个变化环境，用来验证所提算法的鲁棒性。为模拟变化环境中的环境影响因素，在环境相关字典参数(A_l和c_l)中加入高斯白噪声。在仿真中，A_l和c_l的值可计算为

$ A_{l}\left(E_{\mathrm{d}}\right)=A_{0}+\sum\limits_{i=1}^{E_{\mathrm{d}}} \varepsilon_{A}, $

(25)

$ c_{l}\left(E_{\mathrm{d}}\right)=c_{0}+\sum\limits_{i=1}^{E_{\mathrm{d}}} \varepsilon_{c} . $

(26)

其中：A₀和c₀分别表示A_l和c_l的初值；ε_A和ε_c为加性高斯白噪声，其方差分别设为0.5和0.01；E_d表示环境影响因子等级；A_l(E_d)和c_l(E_d) 表示环境影响因子等级为E_d时，所对应的环境相关字典参数值。

根据以上参数，可以建立一个基于变化环境的无源字典，得到模拟变化环境的测量值。同时，将由克拉美罗下界(Cramer-Rao lower bound, CRLB)推导而来的定位误差作为各算法的对比参照，具体推导过程见附录。图 4绘制了DFL-RSS和所提算法DFL-PS的平均定位误差随环境影响因子等级E_d变化的曲线。从中可以看出，当E_d从1增加到9时，2种无源定位方法的平均定位误差都相应增加。DFL-PS的定位精度明显高于DFL-RSS，且随着环境的变化，DFL-PS的平均定位误差增长幅度小于DFL-RSS。由此可见，DFL-PS具有更好的鲁棒性。

然后比较不同目标数下2种算法的定位精度。图 5展示了不同目标数K下，2种算法估计目标位置和真实目标位置的差异。从中可以看出，随着K的增大，2种算法的定位精度都有所下降。这是因为K的增大使目标位置向量θ的稀疏性变差，导致稀疏恢复算法性能下降。但任意目标数下，DFL-PS的定位精度均高于DFL-RSS。

最后比较不同大小定位区域下，DFL-PS和DFL-RSS的定位性能。表 2展示了DFL-PS和DFL-RSS在6.5 m×6.5 m定位区域和14 m×14 m定位区域中的平均定位误差和平均定位时间(平均超过5 000次实验)。结果表明，DFL-PS在不同大小定位区域中均比DFL-RSS更精确，且DFL-PS采用的网格裁剪机制有效降低了定位时间。

4 结束语

针对传统基于RSS的压缩感知无源定位方法定位精度不高的问题，本文设计实现了基于相位偏移的压缩感知无源多目标定位方法DFL-PS。该方法将接收信号相位偏移值作为观测数据，结合变分贝叶斯推理，恢复目标位置稀疏向量，估计目标位置。同时进行网格裁剪，降低运算复杂度，节省计算成本。仿真实验结果表明，相较于DFL-RSS，该方法定位精度有较大提升(在6.5 m×6.5 m的监测区域中，定位精度提升超过1倍；在14 m×14 m的监测区域中，定位精度提升接近50%)，且平均定位时间仅为前者的1/3。同时，该方法在变化环境中，定位性能没有明显下降，具有较强的鲁棒性。