CN105554881A - 一种非直达波环境下的gsm-r网络干扰源定位方法 - Google Patents

Abstract

本发明属于铁路无线通信网干扰源定位技术领域，涉及一种非直达波环境下的GSM-R网络干扰源定位方法。本发明主要包括：先从GSM-R网络中获取n个到达时间测量值和n个到达角度测量值，构建干扰源与移动站之间的几何关系模型；根据几何关系模型，在忽略非直达波误差的条件下，获取干扰源位置坐标估计值和误差矢量；根据获得的干扰源位置坐标估计值和误差矢量，构建用于抑制非直达波误差的非线性加权代价函数，然后获取抑制非直达波误差后的干扰源位置坐标。本发明能带的有益效果为，消除了非直达波误差对干扰源定位结果的影响，进一步提高了干扰源在非直达波环境下的定位精度，弥补了现有的GSM-R网络中非直达波环境下干扰源定位精度不高的缺陷。

Description

一种非直达波环境下的GSM-R网络干扰源定位方法

技术领域

本发明属于铁路无线通信网干扰源定位技术领域，涉及一种非直达波环境下的GSM-R网络干扰源定位方法。

背景技术

GSM-R(GlobalSystemforMobileCommunications-Railway)是专门为铁路通信设计的综合专用铁路无线通信系统。由国际铁路联盟UIC与欧洲电信标准组织ETSI为欧洲各国新一代铁路无线通信开发的技术标准。GSM-R在GSM的标准基础上增加了调度通信功能，并且适合高速环境下使用，为铁路安全高效运营提供了一种可以定制附加功能的专用数字综合移动通信平台，使得用户可以在这个信息平台上轻松开发各种各样的铁路应用。

GSM-R系统是我国铁路通信的发展趋势，它能够实现铁路各种移动信息资源采集、传输，为现代化的铁路调度、指挥、控制提供通信平台，同时提高了铁路的运转效率，减小了运行成本。2003年，铁道部选定GSM-R作为我国铁路未来综合数字调度网的技术体制。

随着我国电子工业和无线通信系统的迅猛发展和普及，大量的电子设备产品应用于社会中，在造福百姓的同时也会给一些专用频段的通信系统带来频率干扰的问题。各种电子产品和通信产品会发出各种频率不同、信号量不同的电磁波并辐射到周围环境中，使空间中的电磁环境变得越来越复杂，往往两系统共用一个频段或者两系统的频段非常接近，所以系统之间的电磁干扰问题已经变得日益突出。我国GSM-R的工作频段主要工作在885-889/930-934MHz范围内，共4MHz带宽，它虽然独享这段带宽，但随时都会受到如移动、联通、电信和无线电台由于设备老化或故障造成发射的无线电信号产生频率漂移、杂波等产生的邻频干扰或互调干扰，还有非法不授权的基站或天线造成的同频干扰。这些无线干扰信号对GSM-R正常通信的干扰，将会造成系统通信掉话、通话质量差、信道拥塞等问题，这些问题会给铁路的安全运行以及旅客的生命安全和国家的发展带来隐忧。因此，有必要发展高效的GSM-R干扰源定位算，消除干扰源对GSM-R系统造成的影响，确保中国铁路专用无线通信系统的安全运行。

目前GSM-R网络干扰源的定位算法主要分为测距定位算法和非测距定位算法两类，非测距算法主要包括质心算法、Dv-Hop算法和APIT算法等。由于非测距算法只能提供粗精度的定位服务，因此发展了基于测距的定位算法。基于测距的定位算法是通过测量网络中节点间的距离或方向角度以及信号强度来计算定位节点坐标的。主要测量技术有AOA(arrivalofangle到达角度)、RSS(receivedsignalstrength接收信号强度)，TOA(timeofarrival到达时间)、TDOA(timedifferenceofarrival到达时差)等。定位算法主要包括基于TOA和TDOA的定位算法以及基于TDOA/AOA，TDOA/RSS，TOA/AOA和TOA/RSS等的混合定位算法。在GSM-R网络中，由于我国铁路覆盖面积大，铁路线长，沿线地形复杂，因此对于基于测距的定位算法往往存在非直达波环境。目前虽然已经提出了很多抑制非直达波误差的定位算法，但都存在定位精度不高、运算量大等诸多缺陷。

发明内容

本发明目的是针对现有的GSM-R网络干扰源定位技术中存在的缺点与不足，在TOA/AOA混合定位算法中采用非线性加权代价函数消除非直达波误差的影响，发展了一种非直达波环境下的GSM-R网络干扰源定位方法。

本发明的技术方案是：一种非直达波环境下的GSM-R网络干扰源定位方法，其特征在于，包括以下步骤：

a.采用基于信号到达时间的定位算法和基于信号到达角度的定位算法，从GSM-R网络中获取n个到达时间测量值和n个到达角度测量值；其中，n为GSM-R网络中移动站的数目；构建干扰源与移动站之间的几何关系模型；

b.根据步骤a中获得的几何关系模型，在忽略非直达波误差的条件下，获取干扰源位置坐标估计值和误差矢量；

c.根据步骤b中获得的干扰源位置坐标估计值和误差矢量，构建用于抑制非直达波误差的非线性加权代价函数，然后获取抑制非直达波误差后的干扰源位置坐标。

本发明的详细方法包括：

步骤1：通过获得的n个TOA和n个AOA测量值构建几何模型。

假设GSM-R网络中的n个移动站的位置坐标为(x_i,y_i),i＝1,2,...,n，干扰源的位置坐标为(x,y)，因此可构建如下几何方程组。

$\left\{\begin{array}{c}{r}_i^2={(x-x_i)}^2+{(y-y_i)}^2\\ {\mathrm{tan\θ}}_i=\frac{{\mathrm{sin\θ}}_i}{{\mathrm{cos\θ}}_i}=\frac{y-y_i}{x-x_i}\end{array}\right.,i=1,2,\mathrm{...},n---\left(1\right)$

其中r_i为干扰源到第i个移动站的距离值，θ_i为干扰源到第i个移动站的弧度。

将式(1)展开化简可得：

$\left\{\begin{array}{c}2x_{i}x+2y_{i}y-k=k_i-r_i^2\\ {\mathrm{sin\θ}}_{i}x-{\mathrm{cos\θ}}_{i}y=x_i{\mathrm{sin\θ}}_i-y_i{\mathrm{cos\θ}}_i\end{array}\right.,i=1,\mathrm{...},n---\left(2\right)$

其中，k＝x²+y²，

将式(2)写成矩阵形式为：

Y＝GZ(3)

其中， $G=\left\{\begin{array}{ccc}2x_1& 2y_1& -1\\ .& .& .\\ .& .& .\\ .& .& .\\ 2x_n& 2y_n& -1\\ {\mathrm{sin\θ}}_1& -{\mathrm{cos\θ}}_1& 0\\ .& .& .\\ .& .& .\\ .& .& .\\ {\mathrm{sin\θ}}_n& -{\mathrm{cos\θ}}_n& 0\end{array}\right\},Y=\left\{\begin{array}{c}{x}_1^2+y_1^2-r_1^2\\ .\\ .\\ .\\ x_n^2+y_n^2-r_n^2\\ x_1{\mathrm{sin\θ}}_1-y_1{\mathrm{cos\θ}}_1\\ .\\ .\\ .\\ x_n{\mathrm{sin\θ}}_n-y_n{\mathrm{cos\θ}}_n\end{array}\right\},Z=\left\{\begin{array}{c}x\\ y\\ k\end{array}\right\}.$

步骤2：忽略非直达波误差，对干扰源位置坐标进行初步估计。

步骤1所建立的矩阵方程中，r_i和θ_i为实测的距离值和弧度值，由于GSM-R通信环境比较复杂，存在非直达波环境，因此在实际的测量值中存在非直达波误差，如下所示：

$\{\begin{array}{c}{r}_i={\hat{r}}_i+n_{r_i}+n_{rNLOS}\\ {\mathrm{\θ}}_i={\widehat{\mathrm{\θ}}}_i+n_{{\mathrm{\θ}}_i}+n_{\mathrm{\θ}NLOS}\end{array},i=1,2,\mathrm{...},n---\left(4\right)$

其中，和为真实的距离值和弧度值。和分别为测量的距离误差和弧度误差，服从均值为零的高斯分布，方差分别为和n_rNLOS和n_θNLOS为测量距离和弧度的非直达波误差，服从指数分布。

将式(4)代入式(3)中可得残差为：

e＝Y-GZ(5)

式(5)的最大似然估计值为：

Z＝(G^Tcov(e)^-1G)^-1G^Tcov(e)^-1Y(6)

其中，cov(e)为误差矢量的协方差矩阵：

cov(e)＝E(ee^T)(7)

将式(5)展开可得：

$\begin{array}{c}{e}_i=x_i^2+y_i^2-r_i^2-2x_{i}x-2y_{i}y+k\\ =x_i^2+y_i^2-{({\hat{r}}_i+n_{r_i})}^2-2x_{i}x-2y_{i}y+k\\ =x_i^2+y_i^2-{\hat{r}}_i^2-2{\hat{r}}_i{n}_{r_i}-n_{r_i}^2-2x_{i}x-2y_{i}y+k\\ e_{n+i}=x_i{\mathrm{sin\θ}}_i-y_i{\mathrm{cos\θ}}_i-x{\mathrm{sin\θ}}_i+y{\mathrm{cos\θ}}_i\\ =(x_i-x)\sin ({\widehat{\mathrm{\θ}}}_i+n_{{\mathrm{\θ}}_i})-(y_i-y)\cos ({\widehat{\mathrm{\θ}}}_i+n_{{\mathrm{\θ}}_i})\\ =(x_i-x)(\sin {\widehat{\mathrm{\θ}}}_i\cos n_{{\mathrm{\θ}}_i}+\sin n_{{\mathrm{\θ}}_i}\cos {\widehat{\mathrm{\θ}}}_i)-(y_i-y)(\cos {\widehat{\mathrm{\θ}}}_i\cos n_{{\mathrm{\θ}}_i}-\sin {\widehat{\mathrm{\θ}}}_i\sin n_{{\mathrm{\θ}}_i})\\ =\[(x_i-x)\sin {\widehat{\mathrm{\θ}}}_i-(y_i-y)\cos {\widehat{\mathrm{\θ}}}_i\]\cos n_{{\mathrm{\θ}}_i}+\[(x_i-x)\cos {\widehat{\mathrm{\θ}}}_i+(y_i-y)\sin {\widehat{\mathrm{\θ}}}_i\]\sin n_{{\mathrm{\θ}}_i}\\ i=1,2,\mathrm{...},n\end{array}---\left(8\right)$

由式(3)可知 $x_i^2+y_i^2-{\hat{r}}_i^2=2x_{i}x+2y_{i}y-k,(x_i-x)\sin {\widehat{\mathrm{\θ}}}_i=(y_i-y)\cos {\widehat{\mathrm{\θ}}}_i,$ 因此对式(8)化简并忽略高次项可得：

$\begin{array}{c}{e}_i=-2{\hat{r}}_i{n}_{r_i}\\ e_{n+i}=\[(x_i-x)\cos {\widehat{\mathrm{\θ}}}_i+(y_i-y)\sin {\widehat{\mathrm{\θ}}}_i\](n_{{\mathrm{\θ}}_i}-\frac{n_{{\mathrm{\θ}}_i}^3}{3!}+\mathrm{...}+{\left(-1\right)}^n\frac{n_{{\mathrm{\θ}}_i}^{2n+1}}{(2n+1)!})\\ =\[(x_i-x)\cos {\widehat{\mathrm{\θ}}}_i+(y_i-y)\sin {\widehat{\mathrm{\θ}}}_i\]n_{{\mathrm{\θ}}_i}\\ i=1,2,\mathrm{...},n\end{array}---\left(9\right)$

由式(9)可得e_ie_j，e_n+ie_n+j，和的期望值为：

$E\left(e_i{e}_j\right)=E\[(-2{\hat{r}}_i{n}_{r_i})(-2{\hat{r}}_j{n}_{r_j})\]=0,i,j=1,2,...,n,i\≠j---\left(10\right)$

$\begin{array}{c}E\left(e_{n+i}{e}_{n+j}\right)=E\{\[(x_i-x)\cos {\widehat{\mathrm{\θ}}}_i+(y_i-y)\sin {\widehat{\mathrm{\θ}}}_i\]n_{{\mathrm{\θ}}_i}\[(x_j-x)\cos {\widehat{\mathrm{\θ}}}_j+(y_j-y)\sin {\widehat{\mathrm{\θ}}}_j\]n_{{\mathrm{\θ}}_j}\}\\ =0,i,j=1,2,\mathrm{...},n,i\≠j\end{array}---\left(11\right)$

$E\left(e_i^2\right)=E\[{(-2{\hat{r}}_i{n}_{r_i})}^2\]=4{\hat{r}}_i^{2}E\left(n_{r_i}^2\right)=4{\hat{r}}_i^2{\mathrm{\σ}}_r^2,i=1,2,...,n---\left(12\right)$

$\begin{array}{c}E\left(e_{n+i}^2\right)=E\left\{{\[\left(x_i-x\right)\cos {\widehat{\mathrm{\θ}}}_i+(y_i-y)\sin {\widehat{\mathrm{\θ}}}_i\]}^2{n}_{{\mathrm{\θ}}_i}^2\right\}\\ ={\[(x_i-x)\cos {\widehat{\mathrm{\θ}}}_i+(y_i-y)\sin {\widehat{\mathrm{\θ}}}_i\]}^{2}E\left(n_{{\mathrm{\θ}}_i}^2\right)\\ ={\[(x_i-x)\cos {\widehat{\mathrm{\θ}}}_i+(y_i-y)\sin {\widehat{\mathrm{\θ}}}_i\]}^2{\mathrm{\σ}}_{\mathrm{\θ}}^2\\ i=1,2,\mathrm{...},n\end{array}---\left(13\right)$

将式(10),(11),(12),(13)代入式(7)可得误差矢量的协方差矩阵为：

$\mathrm{cov}\left(e\right)=\left\{\begin{array}{cc}{B}_r{Q}_r{B}_r& \\ & B_{\mathrm{\θ}}{Q}_{\mathrm{\θ}}{B}_{\mathrm{\θ}}\end{array}\right\}---\left(14\right)$

其中，B_r＝2diag{[r₁…r_n]}， $Q_r=diag\left\{\left[\begin{array}{ccc}{\mathrm{\σ}}_{r1}^2& ...& {\mathrm{\σ}}_{rn}^2\end{array}\right]\right\},$

$Q_{\mathrm{\θ}}=diag\left\{\left[\begin{array}{ccc}{\mathrm{\σ}}_{\mathrm{\θ}1}^2& ...& {\mathrm{\σ}}_{\mathrm{\θ}n}^2\end{array}\right]\right\},$

$B_{\mathrm{\θ}}=diag\left\{\left[\begin{array}{ccc}(x_1-x)cos{\widehat{\mathrm{\θ}}}_1+(y_1-y)sin{\widehat{\mathrm{\θ}}}_1& ...& (x_n-x)cos{\widehat{\mathrm{\θ}}}_n+(y_n-y)sin{\widehat{\mathrm{\θ}}}_n\end{array}\right]\right\},$ B_θ中含有未知量x和y，可先用TOA进行初步估计，再将估计值代入B_θ，即：

$Z_{toa}={\left(G_{toa}^T{B_r}^{-1}{G}_{toa}\right)}^{-1}{G}_{toa}^T{B_r}^{-1}{Y}_{toa}---\left(15\right)$

其中， $G_{toa}=\left\{\begin{array}{ccc}2x_1& 2y_1& -1\\ .& .& .\\ .& .& .\\ .& .& .\\ 2x_n& 2y_n& -1\end{array}\right\},Y_{toa}=\left\{\begin{array}{c}{x}_1^2+y_1^2-r_1^2\\ .\\ .\\ .\\ x_n^2+y_n^2-r_n^2\end{array}\right\},Z_{toa}=\left\{\begin{array}{c}x\\ y\\ k\end{array}\right\}.$

步骤3：通过干扰源位置坐标初步估计值和误差矢量，构建非线性加权代价函数，抑制非直达波误差。

常规的定位方法通常是基于均方代价函数或加权均方代价函数，而这些定位方法在只存在高斯噪声的直达波环境中是最优的，而在非直达波环境中，由于非直达波误差分布往往偏离高斯分布，这些常规方法失去了最优性。而且均方代价函数随着误差的增长呈抛物线大幅度增长，对均方代价函数的最小化并不能抑制大误差，因此有必要对其改进，所以可以采用非线性加权代价函数抑制非直达波误差提高精确定位。

首先，构建非线性加权代价函数：

$J=\underset{i=1}{\overset{n}{\Σ}}{w}_i\mathrm{\ρ}\left(e_i\right)---\left(16\right)$

其中，ρ(e_i)为抑制非直达波误差的非线性加权代价函数，w_i为权值。

在非直达波环境中对干扰源进行定位需要一种对大误差不敏感的代价函数。Talvar估计器就是这类函数，具有良好的非线性，能够通过权值抑制大误差。因此本发明采用Talvar函数，该函数定义为：

其中，β为门限值，可通过误差e_i的绝对中位差求取：

β＝med{|e_i-med{e_i}|}(18)

将步骤2中式(6)所求得的干扰源位置坐标代入式(5)可求得误差矢量，再将其代入式(18)中，即可求得非线性代价函数的门限值。

定义影响函数：

定义加权函数：

将式(5)代入式(16)中，并对其求导使其最小化，可得：

$\underset{i=1}{\overset{n}{\Σ}}{w}_i{G}_i^T\mathrm{\ψ}\left(e_i\right)=0---\left(21\right)$

将式(20)代入式(21)，可得：

$\underset{i=1}{\overset{n}{\Σ}}{w}_i{G}_i^T(Y_i-G_{i}Z)q\left(e_i\right)=0---\left(22\right)$

化简式(22)，可得：

$\underset{i=1}{\overset{n}{\Σ}}{w}_i{G}_i^{T}q\left(e_i\right)Y_i=\underset{i=1}{\overset{n}{\Σ}}{G}_i^T{G}_i{w}_{i}q\left(e_i\right)Z---\left(23\right)$

将式(23)写成矩阵形式，可得：

G^TWQY＝G^TWQGZ(24)

其中，Q＝diag{[q(e₁)q(e₂)…q(e_n)]}，加权矩阵W为误差矢量的协方差的逆：

W＝cov(e)^-1(25)

由式(24)通过变换可得到抑制非直达波误差的干扰源位置坐标估计值：

Z＝(G^TWQG)^-1G^TWQY(26)

步骤4：精确定位干扰源。

通过式(26)求得的干扰源位置坐标未考虑x,y,k之间的相关性，由于k＝x²+y²，因此可以进一步优化干扰源位置坐标估计值。

首先，假设x,y,k的估计误差为μ₁,μ₂,μ₃，因此式(26)求得的估计值可表示为：

[z]₁＝x+μ₁，[z]₂＝y+μ₂，[z]₃＝k+μ₃(27)

根据式(27)可构造如下误差矢量：

e'＝Y'-G'Z'(28)

其中 $\begin{array}{cccc}{Y}^{\′}=\left[\begin{array}{c}{\left({\[Z\]}_1^2\right)}^2\\ {\left({\[Z\]}_2^2\right)}^2\\ {\[Z\]}_3\end{array}\right],& G^{\′}=\left[\begin{array}{cc}1& 0\\ 0& 1\\ 1& 1\end{array}\right],& Z^{\′}=\left[\begin{array}{c}{x}^2\\ y^2\end{array}\right]& e^{\′}=\left[\begin{array}{c}{e_1}^{\′}\\ {e_2}^{\′}\\ {e_3}^{\′}\end{array}\right]\end{array}.$

将式(27)代入到式(28)，并忽略二次项可得：

$\begin{array}{c}{e^{\′}}_1=2{\mathrm{x\μ}}_1+{\mathrm{\μ}}_1^2\≈2{\mathrm{x\μ}}_1\\ {e^{\′}}_2=2{\mathrm{y\μ}}_2+{\mathrm{\μ}}_2^2\≈2{\mathrm{y\μ}}_2\\ {e^{\′}}_3={\mathrm{\μ}}_3\end{array}---\left(29\right)$

由式(29)可得e'的协方差矩阵为：

cov(e')＝E(e'e'^T)＝D{cov(Z)}D(30)

其中，D＝diag{[2x2y1]}≈diag{[2[z]₁2[z]₂1]}，cov(Z)为式(26)求得的Z的协方差矩阵，可通过扰动分析方法求解。定义Δ为误差扰动分量，则式(26)可写成：

$\left(\right({\hat{G}}^T+{\mathrm{\ΔG}}^T\left)WQ\right(\hat{G}+\mathrm{\Δ}G\left)\right)(\hat{Z}+\mathrm{\Δ}Z)=({\hat{G}}^T+{\mathrm{\ΔG}}^T)WQ(\hat{Y}+\mathrm{\Δ}Y)---\left(31\right)$

由于因此式(31)展开并忽略高次项可得：

${\hat{G}}^{T}WQ\hat{G}\mathrm{\Δ}Z+{\hat{G}}^{T}WQ\mathrm{\Δ}G\hat{Z}={\hat{G}}^{T}WQ\mathrm{\Δ}Y---\left(32\right)$

对式(32)化简移项得：

$\mathrm{\Δ}Z={\left({\hat{G}}^{T}WQ\hat{G}\right)}^{-1}{\hat{G}}^{T}WQ(\mathrm{\Δ}Y-\mathrm{\Δ}G\hat{Z})={\left({\hat{G}}^{T}W\hat{G}\right)}^{-1}{\hat{G}}^{T}W(\mathrm{\Δ}Y-\mathrm{\Δ}G\hat{Z})---\left(33\right)$

由于所以由式(33)可得：

$\mathrm{\Δ}Z={\left({\hat{G}}^{T}W\hat{G}\right)}^{-1}{\hat{G}}^{T}We---\left(34\right)$

由式(34)可得cov(Z)为：

$\begin{array}{c}\mathrm{cov}\left(Z\right)=E\[{\mathrm{\ΔZ\ΔZ}}^T\]\\ ={\left({\hat{G}}^{T}W\hat{G}\right)}^{-1}{\hat{G}}^{T}WE\left({\mathrm{ee}}^T\right)W\hat{G}{\left({\hat{G}}^{T}W\hat{G}\right)}^{-1}\\ ={\left({\hat{G}}^{T}W\hat{G}\right)}^{-1}\end{array}---\left(35\right)$

对式(28)求最大似然估计可得：

$\begin{array}{c}{Z}^{\′}=\arg \min \left\{{(Y^{\′}-G^{\′}{Z}^{\′})}^T\mathrm{cov}{\left(e^{\′}\right)}^{-1}(Y^{\′}-G^{\′}{Z}^{\′})\right\}\\ ={\left(G^{\′T}\mathrm{cov}{\left(e^{\′}\right)}^{-1}{G}^{\′}\right)}^{-1}{G}^{\′T}\mathrm{cov}{\left(e^{\′}\right)}^{-1}{Y}^{\′}\end{array}---\left(36\right)$

由式(36)最终可求得优化后的干扰源位置坐标估计值为：

$Z^{\′\′}=\left[\begin{array}{c}x\\ y\end{array}\right]=\left[\begin{array}{c}sign\left({\[Z\]}_1\right)\sqrt{{\[Z^{\′}\]}_1}\\ sign\left({\[Z\]}_2\right)\sqrt{{\[Z^{\′}\]}_2}\end{array}\right]---\left(37\right)$

本发明能带的有益效果为，本发明利用TOA/AOA混合定位算法获得干扰源位置坐标的初步估计值和误差矢量，构建了非线性加权代价函数，消除了非直达波误差对干扰源定位结果的影响，进一步提高了干扰源在非直达波环境下的定位精度，弥补了现有的GSM-R网络中非直达波环境下干扰源定位精度不高的缺陷。

附图说明

图1本发明的干扰源定位流程图；

图2本发明的移动站和干扰源坐标随机部署图；

图3本发明的干扰源位置坐标定位效果图；

图4本发明的干扰源位置坐标定位误差图。

具体实施方式

下面结合附图和实施例对本发明进行详细的描述

如图1所示，本发明的流程为：首先对从GSM-R网络中获得的TOA和AOA测量值进行几何建模，之后忽略非直达波误差对干扰源定位结果的影响，进行初步估计，然后通过获得的初步估计值和定位误差，构建非线性加权代价函数，最后对干扰源进行精确定位。

实施例

如图2所示，GSM-R网络二维平面上，部署了8个移动站用于测量干扰源的TOA和AOA值，坐标分别为(0,0)、(0,5000)、(1250,2500)、(2500,5000)、(2500,0)、(3750,2500)、(5000,5000)、(5000,0)，以(x_i,y_i),i＝1,2,...,8表示。干扰源随机分布其中，本次实施例中干扰源位置坐标为(3254.8839,2035.6002)，以(x,y)表示。

步骤1：通过获得的n个TOA和n个AOA测量值构建几何模型。

首先对移动站和干扰源进行几何建模，通过获得的8个TOA和8个AOA测量值，构建如下几何方程组：

$\{\begin{array}{c}{r}_i^2={(x-x_i)}^2+{(y-y_i)}^2\\ {\mathrm{tan\θ}}_i=\frac{{\mathrm{sin\θ}}_i}{{\mathrm{cos\θ}}_i}=\frac{y-y_i}{x-x_i}\end{array},i=1,2,\mathrm{...},8---\left(38\right)$

其中r_i为干扰源到第i个移动站的距离值，θ_i为干扰源到第i个移动站的弧度。

对式(38)展开化简可得：

$\{\begin{array}{c}2x_{i}x+2y_{i}y-k=k_i-{r_i}^2\\ {\mathrm{sin\θ}}_{i}x-{\mathrm{cos\θ}}_{i}y=x_i{\mathrm{sin\θ}}_i-y_i{\mathrm{cos\θ}}_i\end{array},i=1,\mathrm{...},8---\left(39\right)$

其中，k＝x²+y²，

将式(39)写成矩阵形式为：

Y＝GZ(40)

其中， $G=\left\{\begin{array}{ccc}2x_1& 2y_1& -1\\ .& .& .\\ .& .& .\\ .& .& .\\ 2x_8& 2y_8& -1\\ {\mathrm{sin\θ}}_1& -{\mathrm{cos\θ}}_1& 0\\ .& .& .\\ .& .& .\\ .& .& .\\ {\mathrm{sin\θ}}_8& -{\mathrm{cos\θ}}_8& 0\end{array}\right\},Y=\left\{\begin{array}{c}{x}_1^2+y_1^2-r_1^2\\ .\\ .\\ .\\ x_8^2+y_8^2-r_8^2\\ x_1{\mathrm{sin\θ}}_1-y_1{\mathrm{cos\θ}}_1\\ .\\ .\\ .\\ x_8{\mathrm{sin\θ}}_8-y_8{\mathrm{cos\θ}}_8\end{array}\right\},Z=\left\{\begin{array}{c}x\\ y\\ k\end{array}\right\}.$

步骤2：忽略非直达波误差，对干扰源位置坐标进行初步估计。

步骤1所建立的矩阵方程中，r_i和θ_i为实测的距离值和弧度值，由于GSM-R通信环境比较复杂，存在非直达波环境，因此在实际的测量值中存在非直达波误差，如下所示：

$\{\begin{array}{c}{r}_i={\hat{r}}_i+n_{r_i}+n_{rNLOS}\\ {\mathrm{\θ}}_i={\widehat{\mathrm{\θ}}}_i+n_{{\mathrm{\θ}}_i}+n_{\mathrm{\θ}NLOS}\end{array},i=1,2,\mathrm{...},8---\left(41\right)$

其中，和为真实的距离值和弧度值。和分别为测量的距离误差和弧度误差，服从均值为零的高斯分布，标准差分别为σ_r＝10和σ_θ＝0.02。n_rNLOS和n_θNLOS为测量距离和弧度的非直达波误差，服从参数为50和0.2的指数分布。

将式(41)代入式(39)中可得残差为：

e＝Y-GZ(42)

式(42)的协方差矩阵为：

$\mathrm{cov}\left(e\right)=\left\{\begin{array}{cc}{B}_r{Q}_r{B}_r& \\ & B_{\mathrm{\θ}}{Q}_{\mathrm{\θ}}{B}_{\mathrm{\θ}}\end{array}\right\}---\left(43\right)$

其中，B_r＝2diag{[r₁…r₈]}， $Q_r=diag\left\{\left[\begin{array}{ccc}{\mathrm{\σ}}_{r1}^2& ...& {\mathrm{\σ}}_{r8}^2\end{array}\right]\right\},Q_{\mathrm{\θ}}=diag\left\{\left[\begin{array}{ccc}{\mathrm{\σ}}_{\mathrm{\θ}1}^2& ...& {\mathrm{\σ}}_{\mathrm{\θ}8}^2\end{array}\right]\right\},$ $B_{\mathrm{\θ}}=diag\left\{\left[\begin{array}{ccc}(x_1-x)cos{\widehat{\mathrm{\θ}}}_1+(y_1-y)sin{\widehat{\mathrm{\θ}}}_1& ...& (x_8-x)cos{\widehat{\mathrm{\θ}}}_8+(y_8-y)sin{\widehat{\mathrm{\θ}}}_8\end{array}\right]\right\},$ B_θ中含有未知量x和y，可先用TOA进行初步估计，再将估计值代入B_θ，即：

$Z_{toa}={\left(G_{toa}^T{B_r}^{-1}{G}_{toa}\right)}^{-1}{G}_{toa}^T{B_r}^{-1}{Y}_{toa}=\left[\begin{array}{c}3.2665\×{10}^3\\ 3.0442\×{10}^3\\ 1.4870\×{10}^7\end{array}\right]---\left(44\right)$

其中， $G_{toa}=\left\{\begin{array}{ccc}2x_1& 2y_1& -1\\ .& .& .\\ .& .& .\\ .& .& .\\ 2x_8& 2y_8& -1\end{array}\right\},Y_{toa}=\left\{\begin{array}{c}{x}_1^2+y_1^2-r_1^2\\ .\\ .\\ .\\ x_8^2+y_8^2-r_8^2\end{array}\right\},Z_{toa}=\left\{\begin{array}{c}x\\ y\\ k\end{array}\right\}.$

由式(42)和式(43)可得忽略非直达波误差的干扰源位置坐标粗略估计值为：

$Z={\left(G^T\mathrm{cov}{\left(e\right)}^{-1}G\right)}^{-1}{G}^T\mathrm{cov}{\left(e\right)}^{-1}Y=\left[\begin{array}{c}3.2666\×{10}^3\\ 2.0353\×{10}^3\\ 1.4828\×{10}^7\end{array}\right]---\left(45\right)$

将式(45)求得的干扰源位置坐标粗略估计值代入到式(42)中，得到定位误差矢量：

$e=Y-GZ=\left[\begin{array}{c}-2.0068\×{10}^5\\ 0.6878\×{10}^5\\ 0.4965\×{10}^5\\ -0.5280\×{10}^5\\ 0.8343\×{10}^5\\ 0.0321\×{10}^5\\ 0.0491\×{10}^5\\ -0.8894\×{10}^5\\ 0.0030\×{10}^5\\ -0.0047\×{10}^5\\ 0.0001\×{10}^5\\ -0.0039\×{10}^5\\ -0.0017\×{10}^5\\ -0.0004\×{10}^5\\ -0.0070\×{10}^5\\ -0.0033\×{10}^5\end{array}\right]---\left(46\right)$

步骤3：通过干扰源位置坐标初步估计值和误差矢量，构建非线性加权代价函数，抑制非直达波误差。

构建非线性加权代价函数：

$J=\underset{i=1}{\overset{8}{\Σ}}{w}_i\mathrm{\ρ}\left(e_i\right)---\left(47\right)$

其中，ρ(e_i)为抑制非直达波误差的非线性代价函数，e_i为步骤2求得的定位误差，w_i为权值。

在非直达波环境中对干扰源进行定位需要一种对大误差不敏感的代价函数。Talvar估计器就是这类函数，具有良好的非线性，能够通过权值抑制大误差。

由式(46)求得Talvar函数所需的门限值为：

$\mathrm{\β}=med\left\{\right|e_i-med\left\{e_i\right\}\left|\right\}=\{\begin{array}{c}5.6792\×{10}^4,i=1,2,\mathrm{...},8\\ 559.6836,i=9,10,\mathrm{...},16\end{array}---\left(48\right)$

由式(48)可得Talvar函数为：

有式(49)可得影响函数为：

由式(50)可得加权函数为：

将式(42)代入式(47)中，并对其求导使其最小化，可得：

$\underset{i=1}{\overset{8}{\Σ}}{w}_i{G}_i^T\mathrm{\ψ}\left(e_i\right)=0---\left(52\right)$

将式(51)代入到式(52)，可得：

$\underset{i=1}{\overset{8}{\Σ}}{w}_i{G}_i^T(Y_i-G_{i}Z)q\left(e_i\right)=0---\left(53\right)$

化简式(53)，可得：

$\underset{i=1}{\overset{8}{\Σ}}{w}_i{G}_i^{T}q\left(e_i\right)Y_i=\underset{i=1}{\overset{8}{\Σ}}{G}_i^T{G}_i{w}_{i}q\left(e_i\right)Z---\left(54\right)$

将式(54)写成矩阵形式，可得：

G^TWQY＝G^TWQGZ(55)

其中，Q＝diag{[q(e₁)q(e₂)…q(e_n)]}，W＝cov(e)^-1。

由式(55)通过变换可得抑制非直达波误差的干扰源位置坐标估计值为：

$Z={\left(G^{T}WQG\right)}^{-1}{G}^{T}WQY=\left[\begin{array}{c}3.2582\×{10}^3\\ 2.0299\×{10}^3\\ 1.4740\×{10}^7\end{array}\right]---\left(56\right)$

式(56)的协方差矩阵为：

$\mathrm{cov}\left(Z\right)={\left({\hat{G}}^{T}W\hat{G}\right)}^{-1}---\left(57\right)$

步骤4：精确定位干扰源。

通过式(56)求得的干扰源位置坐标未考虑x,y,k之间的相关性，由于k＝x²+y²，因此可以进一步优化干扰源位置坐标估计值。

假设x,y,k的估计误差为μ₁,μ₂,μ₃，因此式(56)求得的估计值可表示为：

[z]₁＝x+μ₁，[z]₂＝y+μ₂，[z]₃＝k+μ₃(58)

根据式(58)可构造如下误差矢量：

e'＝Y'-G'Z'(59)

其中 $\begin{array}{cccc}{Y}^{\′}=\left[\begin{array}{c}{\left({\[Z\]}_1^2\right)}^2\\ {\left({\[Z\]}_2^2\right)}^2\\ {\[Z\]}_3\end{array}\right],& G^{\′}=\left[\begin{array}{cc}1& 0\\ 0& 1\\ 1& 1\end{array}\right],& Z^{\′}=\left[\begin{array}{c}{x}^2\\ y^2\end{array}\right]& e^{\′}=\left[\begin{array}{c}{e_1}^{\′}\\ {e_2}^{\′}\\ {e_3}^{\′}\end{array}\right]\end{array},$

则式(59)的协方差矩阵为：

cov(e')＝E(e'e'^T)＝D{cov(Z)}D(60)

其中，D＝diag{[2x2y1]}≈diag{[2[z]₁2[z]₂1]}。

对式(59)求最大似然估计可得：

$Z^{\′}={\left(G^{\′T}\mathrm{cov}{\left(e^{\′}\right)}^{-1}{G}^{\′}\right)}^{-1}{G}^{\′T}\mathrm{cov}{\left(e^{\′}\right)}^{-1}{Y}^{\′}=\left[\begin{array}{c}1.0614\×{10}^7\\ 0.4118\×{10}^7\end{array}\right]---\left(61\right)$

由式(61)最终可求得优化后的干扰源的位置坐标估计值为：

$Z^{\′\′}=\left[\begin{array}{c}x\\ y\end{array}\right]=\left[\begin{array}{c}3.2579\×{10}^3\\ 2.0294\×{10}^3\end{array}\right]---\left(62\right)$

式(62)即为最终的抑制非直达波误差的干扰源定位结果，定位效果图如图3所示，从图3可以看出本发明的干扰源定位结果非常接近干扰源的真实坐标值，有效地抑制了非直达波误差对定位结果的影响。

为了进一步说明本发明方法的定位效果，在图4中给出了测量距离在不同的非直达波误差环境下的定位误差效果。图4中测量距离的非直达波误差服从参数分别为50、120、200、280、350的指数分布，从图4中可以看出本发明方法在非直达波环境下存在较小的定位误差。

从上述验证结果可以看出，采用本发明方法可以有效地抑制非直达波误差估计出干扰源位置坐标，定位效果较好。

Claims (5)

1.一种非直达波环境下的GSM-R网络干扰源定位方法，其特征在于，包括以下步骤：

a.采用基于信号到达时间的定位算法和基于信号到达角度的定位算法，从GSM-R网络中获取n个到达时间测量值和n个到达角度测量值；其中，n为GSM-R网络中移动站的数目；构建干扰源与移动站之间的几何关系模型；

b.根据步骤a中获得的几何关系模型，在忽略非直达波误差的条件下，获取干扰源位置坐标估计值和误差矢量；

c.根据步骤b中获得的干扰源位置坐标估计值和误差矢量，构建用于抑制非直达波误差的非线性加权代价函数，然后获取抑制非直达波误差后的干扰源位置坐标。

2.根据权利要求1所述的一种非直达波环境下的GSM-R网络干扰源定位方法，其特征在于，步骤a中所述几何关系模型的具体构建方法为：

假设GSM-R网络中的n个移动站的位置坐标为(x_i,y_i),i＝1,2,...,n，干扰源的位置坐标为(x,y)，其中，下标i指代移动站的序号；构建如下几何方程组：

$\left\{\begin{array}{c}{r}_i^2={(x-x_i)}^2+{(y-y_i)}^2\\ {\mathrm{tan\θ}}_i=\frac{{\mathrm{sin\θ}}_i}{{\mathrm{cos\θ}}_i}=\frac{y-y_i}{x-x_i}\end{array}\right.$

其中r_i为干扰源到第i个移动站的距离值，θ_i为干扰源到第i个移动站的弧度；

该几何方程组可展开化简为：

$\left\{\begin{array}{c}2x_{i}x+2y_{i}y-k=k_i-r_i^2\\ {\mathrm{sin\θ}}_{i}x-{\mathrm{cos\θ}}_{i}y=x_i{\mathrm{sin\θ}}_i-y_i{\mathrm{cos\θ}}_i\end{array}\right.$

其中，k＝x²+y²，

写成矩阵形式为：

Y＝GZ；

其中， $G=\left\{\begin{array}{ccc}2x_1& 2y_1& -1\\ .& .& .\\ .& .& .\\ .& .& .\\ 2x_n& 2y_n& -1\\ {\mathrm{sin\θ}}_1& -{\mathrm{cos\θ}}_1& 0\\ .& .& .\\ .& .& .\\ .& .& .\\ {\mathrm{sin\θ}}_n& -{\mathrm{cos\θ}}_n& 0\end{array}\right\},$ $Y=\left\{\begin{array}{c}{x}_1^2+y_1^2-r_1^2\\ .\\ .\\ .\\ x_n^2+y_n^2-r_n^2\\ x_1{\mathrm{sin\θ}}_1-y_1{\mathrm{cos\θ}}_1\\ .\\ .\\ .\\ x_n{\mathrm{sin\θ}}_n-y_n{\mathrm{cos\θ}}_n\end{array}\right\},$ $Z=\left\{\begin{array}{c}x\\ y\\ k\end{array}\right\}.$

3.根据权利要求2所述的一种非直达波环境下的GSM-R网络干扰源定位方法，其特征在于，步骤b中所述获取干扰源位置坐标估计值和误差矢量的具体方法为：

假设实测的距离值和弧度值与真实的距离值和弧度值之间的关系为：

$\left\{\begin{array}{c}{r}_i={\hat{r}}_i+n_{r_i}+n_{rNLOS}\\ {\mathrm{\θ}}_i={\widehat{\mathrm{\θ}}}_i+n_{{\mathrm{\θ}}_i}+n_{\mathrm{\θ}NLOS}\end{array}\right.$

其中，和为真实的距离值和弧度值，和分别为测量的距离误差和弧度误差，服从均值为零的高斯分布，方差分别为和n_rNLOS和n_θNLOS为测量距离和弧度的非直达波误差，服从指数分布；

结合矩阵Y＝GZ可得误差矢量e为e＝Y-GZ；

忽略非直达波误差后，干扰源位置坐标估计值为e＝Y-GZ的最大似然估计值，可得：Z＝(G^Tcov(e)^-1G)^-1G^Tcov(e)^-1Y，其中cov(e)为误差矢量的协方差矩阵。

4.根据权利要求3所述的一种非直达波环境下的GSM-R网络干扰源定位方法，其特征在于，步骤c中所述构建用于抑制非直达波误差的非线性加权代价函数的具体方法为：

构建非线性加权代价函数：

$J=\underset{i=1}{\overset{n}{\Σ}}{w}_i\mathrm{\ρ}\left(e_i\right)$

其中，ρ(e_i)为抑制非直达波误差的非线性加权代价函数，w_i为权值，采用具有良好的非线性且能够通过权值抑制大误差的Talvar函数作为ρ(e_i)，e_i为误差矢量，定义为：

其中，β为门限值，可通过误差e_i的绝对中位差求取：

β＝med{|e_i-med{e_i}|}

则，定义影响函数ψ(e_i)为：

加权函数为：

将e＝Y-GZ代入构建的加权代价函数可得：

$\underset{i=1}{\overset{n}{\Σ}}{w}_i{G}_i^T\mathrm{\ψ}\left(e_i\right)=0$

结合加权函数可得：

$\underset{i=1}{\overset{n}{\Σ}}{w}_i{G}_i^T(Y_i-G_{i}Z)q\left(e_i\right)=0$

化简后可得：

$\underset{i=1}{\overset{n}{\Σ}}{w}_i{G}_i^{T}q\left(e_i\right)Y_i=\underset{i=1}{\overset{n}{\Σ}}{G}_i^T{G}_i{w}_{i}q\left(e_i\right)Z$

写成矩阵形式为：G^TWQY＝G^TWQGZ；其中，Q＝diag{[q(e₁)q(e₂)…q(e_n)]}，加权矩阵W为误差矢量的协方差的逆：W＝cov(e)^-1；

可得到抑制非直达波误差后的干扰源位置坐标估计值：Z＝(G^TWQG)^-1G^TWQY。

5.根据权利要求4所述的一种非直达波环境下的GSM-R网络干扰源定位方法，其特征在于，还包括步骤：

d.根据获得的抑制非直达波误差后的干扰源位置坐标估计值，消除误差获得更精确的干扰源位置坐标估计值，具体方法为：

假设x,y,k的估计误差为μ₁,μ₂,μ₃，则Z＝(G^TWQG)^-1G^TWQY可表示为：

[z]₁＝x+μ₁，[z]₂＝y+μ₂，[z]₃＝k+μ₃

构成误差矢量e'为：e'＝Y'-G'Z'；

其中， $Y^{\′}=\left[\begin{array}{c}{\left({\[Z\]}_1^2\right)}^2\\ {\left({\[Z\]}_2^2\right)}^2\\ {\[Z\]}_3\end{array}\right],$ $G^{\′}=\left[\begin{array}{cc}1& 0\\ 0& 1\\ 1& 1\end{array}\right],$ $Z^{\′}=\left[\begin{array}{c}{x}^2\\ y^2\end{array}\right]$ $e^{\′}=\left[\begin{array}{c}{e_1}^{\′}\\ {e_2}^{\′}\\ {e_3}^{\′}\end{array}\right];$

将[z]₁＝x+μ₁，[z]₂＝y+μ₂，[z]₃＝k+μ₃分别代入误差矢量e'可得：

${e^{\′}}_1=2{\mathrm{x\μ}}_1+{\mathrm{\μ}}_1^2\≈2{\mathrm{x\μ}}_1$

${e^{\′}}_2=2{\mathrm{y\μ}}_2+{\mathrm{\μ}}_2^2\≈2{\mathrm{y\μ}}_2$

e'₃＝μ₃

可得e'的协方差矩阵为：cov(e')＝E(e'e'^T)＝D{cov(Z)}D；其中，D＝diag{[2x2y1]}≈diag{[2[z]₁2[z]₂1]}，cov(Z)为Z＝(G^TWQG)^-1G^TWQY的协方差矩阵；定义Δ为误差扰动分量，则Z＝(G^TWQG)^-1G^TWQY可写成：

$\left(\right({\hat{G}}^T+{\mathrm{\ΔG}}^T\left)WQ\right(\hat{G}+\mathrm{\Δ}G\left)\right)(\hat{Z}+\mathrm{\Δ}Z)=({\hat{G}}^T+{\mathrm{\ΔG}}^T)WQ(\hat{Y}+\mathrm{\Δ}Y)$

由于将其展开并忽略高次项可得：

${\hat{G}}^{T}WQ\hat{G}\mathrm{\Δ}Z+{\hat{G}}^{T}WQ\mathrm{\Δ}G\hat{Z}={\hat{G}}^{T}WQ\mathrm{\Δ}Y$

化简移项后得：

$\mathrm{\Δ}Z{\left({\hat{G}}^{T}WQ\hat{G}\right)}^{-1}{\hat{G}}^{T}WQ(\mathrm{\Δ}Y-\mathrm{\Δ}G\hat{Z})={\left({\hat{G}}^{T}W\hat{G}\right)}^{-1}{\hat{G}}^{T}W(\mathrm{\Δ}Y-\mathrm{\Δ}G\hat{Z})$

由于 $e=\mathrm{\Δ}Y-\mathrm{\Δ}G\hat{Z},$ 可得：

$\mathrm{\Δ}Z={\left({\hat{G}}^{T}W\hat{G}\right)}^{-1}{\hat{G}}^{T}We$

可得cov(Z)为：

$\begin{array}{c}\mathrm{cov}\left(Z\right)=E\[{\mathrm{\ΔZ\ΔZ}}^T\]\\ ={\left({\hat{G}}^{T}W\hat{G}\right)}^{-1}{\hat{G}}^{T}WE\left({\mathrm{ee}}^T\right)W\hat{G}{\left({\hat{G}}^{T}W\hat{G}\right)}^{-1}\\ ={\left({\hat{G}}^{T}W\hat{G}\right)}^{-1}\end{array}$

对e'＝Y'-G'Z'求最大似然估计可得：

Z'＝argmin{(Y'-G'Z')^Tcov(e')^-1(Y'-G'Z')}

＝(G'^Tcov(e')^-1G')^-1G'^Tcov(e')^-1Y'

可得优化后的干扰源位置坐标估计值为：

$Z^{\′\′}=\left[\begin{array}{c}x\\ y\end{array}\right]=\left[\begin{array}{c}sign\left({\[Z\]}_1\right)\sqrt{{\[Z^{\′}\]}_1}\\ sign\left({\[Z\]}_2\right)\sqrt{{\[Z^{\′}\]}_2}\end{array}\right].$