基于非线性约束最小二乘的单站定位方法与流程

本发明涉及定位技术领域，具体涉及到一种宏蜂窝nlos环境下的单站定位技术。

背景技术：

1996年美国联邦通信委员会公布了e911定位需求，其要求在2001年10月前，各种无线蜂窝网络必须能对发出e911紧急呼叫的移动台提供精度在125m内的定位服务，且满足此定位精度的概率应不低于67％；并要求在2001年以后，提供更高的定位精度及三维位置信息。e911技术规范的公布激发了人们对无线定位技术的研究热情，并且随着通信技术的发展和移动设备越来越高的普及，人们对自身位置的定位服务需求日益激增，基于移动通信网络的定位技术得到飞速发展，无论是在人们的日常出行、交通工具定位导航，还是无线资源分配和整合、国家信息安全等方面都有着举足轻重的地位。当前主流的定位技术包括：广泛应用于室外的卫星定位技术、蜂窝网定位技术、声波定位技术、雷达定位技术等，以及基于wifi、rfid、蓝牙、毫米波等室内定位技术，其中蜂窝定位技术，以其定位速度快、成本低(不需要移动终端上添加额外的硬件)、耗电少等优势，作为一种轻量级的定位方法，越来越常用。

非视距传播是影响蜂窝网定位精度的主要因素，目前对蜂窝网非视距的定位研究主要包括：非直达波鉴别算法；基于散射体信息的非直达波定位算法；基于不等式约束和统计的非直达波定位方法；学习型定位方法和非直达波跟踪算法。其中基于散射体信息的定位是近年来发展起来抑制非直达波误差的新思路，主要包括基于散射体信道模型和基于散射体几何位置关系两个分支。基于散射体信道模型根据散射模型(如圆环，高斯等模型)对测量参数aoa，toa等进行重构，然后利用传统的los定位算法实现目标定位；基于散射体几何位置关系这类算法往往通过先获取散射体的相关几何位置信息，以此来构建出基站、散射体、目标三者之间的几何关系，构造出对应的线性定位方程，一般可同时定位出散射体和目标的位置。

但目前蜂窝网定位的主要存在以下两个问题：nlos传播问题和多基站定位问题。蜂窝网在城区或者郊区等环境下几乎不存在los路径，而nlos传播引起定位参数的测量误差，降低定位精度，同时现有的定位技术多偏向多站定位，其对时间和数据同步要求严格，系统的复杂度高，且存在可测性问题。

因此，有必要开发一种基于非线性约束最小二乘的单站定位方法。

技术实现要素：

本发明的目的是提出一种基于非线性约束最小二乘的单站定位方法，它能解决传统定位算法对基站个数的依赖性问题，利用单基站进行定位，不需要时钟同步，降低算法复杂度；且定位精度高。

本发明所述的基于非线性约束最小二乘的单站定位方法，包括以下步骤：

步骤一：根据宏蜂窝信号传播特性构造宏蜂窝单站的nlos(non-line-of-sight，非视距)信号传播模型；

步骤二：计算基站端每条反射径信号的到达角(angleofarrival,aoa)θi和到达时间(timeofarrival,toa)τi，其中i＝1,...,m，m表示反射径条数；

步骤三：利用多径信号参数aoa,toa及基站、目标和散射体间的位置关系建立定位方程组，此定位方程组为非线性欠定方程组；

步骤四：引入宏蜂窝的圆环模型，假设散射体均匀分布在以移动台(mobilestation,ms)为圆心、半径为r的圆上，移动台和基站(basestation,bs)间距离为d；

步骤五：根据圆环模型增加关于目标位置和散射体位置的非线性方程组，将非线性欠定方程组转化为超定方程组；

步骤六：实际测量时，考虑测量误差，将定位问题转化为非线性约束最小二乘问题；

步骤七：利用lm(levenberg-marquardt，列文伯格-马夸尔特法)算法求解非线性约束最小二乘问题；

步骤八：当残差较大时，引入拟牛顿bfgs(以其发明者broyden,fletcher,goldfarb和shanno的姓氏首字母命名)算法，对hesse(海赛)矩阵进行逼近，再利用lm算法求解，本发明称此改进算法为“lm+bfgs”算法；；

步骤九：算法迭代结束，输出估计的目标位置

进一步，所述步骤三具体为：

设信号经第i个散射体反射后到达基站，经反射后的路径长度为ri，则散射体坐标(xi,yi)为：

建立关于目标位置(xms,yms)和ri的非线性欠定方程组，如公式(2)所示：

其中：c表示光速，方程组中共有m个方程，m+2个未知变量；

将公式(2)简写为：

进一步，所述步骤五具体为：

圆半径r表示为：因此有：

cτ1-r1＝cτj+1-rj+1,j＝1,...,m-1(4)；

联合公式(3)，得到宏蜂窝的圆环模型下求解目标位置和散射体位置的超定方程组，如下式所示：

上式方程组共2m-1个方程，m+2个未知变量，当m>3时，式(5)即为超定方程组。

进一步，所述步骤六具体为：

由于aoa和toa测量误差，公式(5)中等号不总是成立，因此将公式(5)写为：

其中：x＝(xms,yms,ri)^t，χ＝(r1,rj+1)^t，εi(x)和ξj(χ)为残差；

故对目标位置的估计，通过最小化如下目标函数得到：

其中：

同时待估的目标位置(xms,yms)应处于以最小传播距离lmin＝cτmin为半径，bs为圆心的圆内，同时最大最小到达角将ms限制在圆心角为α的扇形区域，扇形的两条半径用向量表示为和长为lmin，α1为与bs和ms所在直线的夹角，α2为与bs和ms所在直线的夹角，α1,α2可以用方向矢量表示为：

其中：

其中：α1(xms,yms)表示α1是关于xms,yms的函数；α2(xms,yms)表示α2是关于xms,yms的函数；

因此，(xms,yms)应处在下列不等式的交集内，即可行域z内：

其中：xbs为基站的横坐标，ybs为基站的纵坐标，同时ri≤cτmax；

综上，式(5)的位置解算转换为求解如下的非线性约束最小二乘问题：

进一步，所述步骤八具体为：

步骤八-一：选取可行域内初始点x∈intz，给定lm算法初始化参数：尺度因子μ>0，增长因子β>0；给定bfgs算法参数：初始对称正定阵β1＝i，其中，ι为单位矩阵；给定终止误差0≤ε＜＜1，令k＝1；

步骤八-二：计算公式(7)的一阶导将j(xk),f(xk)简写为jk,fk，jk为雅克比矩阵，fk＝(f1,f2,...,f2m-1)^t，jk表示为如下形式：

步骤八-三：计算搜索方向dk＝-(bk+μi)-¹gk，其中，bk为对f(x)进行泰勒展开的二阶项的近似矩阵，其初始值为β1；

步骤八-四：由armijio搜索求步长λk，令xk+1＝xk+λkdk；

步骤八-五：计算下一时刻一阶导数计算pk＝xk+1-xk和qk＝gk+1-gk；

步骤八-六：根据bfgs算法更新βk+1，βk+1的更新表达式如下：

其中：βk+1为bk更新后的值；

步骤八-七：计算f(xk+1)，如果f(xk+1)<f(xk)，转入步骤八-八，否则转入步骤八-九；

步骤八-八：如果||dk||2≤ε，得到停止迭代，否则令μ:＝μ/β，同时k:＝k+1，并转入步骤八-二；

步骤八-九：如果||dk||2≤ε，得到停止迭代，否则令μ:＝μβ，同时k:＝k+1，并转入步骤八-二。

进一步，所述步骤八-四具体为：

步骤八-四-一：给定armijio算法的参数ρ∈(0,1),σ∈(0,0.5)，最大迭代次数mmax，令m＝0；

步骤八-四-二：计算f(xk+ρ^mdk)和f(xk)+σρ^mgk^tdk，若f(xk+ρ^mdk)≤f(xk)+σρ^mgk^tdk，转入步骤八-四-四，否则转入步骤八-四-三；

步骤八-四-三：令m:＝m+1，并判断m是否达到最大迭代次数mmax，如否，则转入步骤八-四-二，若是，则停止迭代，并执行步骤八-四-四；

步骤八-四-四：输出搜索步长λk＝ρ^m。

本发明的有益效果：它解决了传统定位算法对基站个数的依赖性问题，利用单基站进行定位，不需要时钟同步，降低了算法的复杂度；同时利用nlos环境下，多径信号作为定位路径，而非抑制多径，解决了nlos环境下定位误差大的问题；利用散射模型增加变量约束条件，提高了最小二乘算法的有效性，实现了宏蜂窝单站定位20-30米的定位精度。

附图说明

图1为本发明整体流程图；

图2为本发明宏蜂窝nlos环境下的散射圆环模型；

图3为目标可能出现的区域(扇形区域)，即迭代算法的可行域；

图4为lm算法流程图；

图5为“lm+bfgs”算法流程图；

图6为armijio算法流程图。

具体实施方式

下面结合具体实施例及附图对本发明作进一步详细说明：

如图1所示，本发明所述的基于非线性约束最小二乘的单站定位方法，包括以下步骤：

步骤一：根据宏蜂窝信号传播特性构造宏蜂窝单站的nlos信号传播模型，如图2所示。假设信号经历单次反射从移动台到达基站，基站位于坐标原点，ms位于以bs为坐标原点的x轴上，散射体s分布在ms周围，宏蜂窝环境下基站天线一般处于较高位置，因此周围不存在散射体。本发明仅考虑los(line-of-sight，视距)不存在的情况。

步骤二：计算基站端每条反射径信号的到达角(aoa)θi和到达时间(toa)τi，其中i＝1,...,m，m表示反射径条数，本发明中m≥3。

步骤三：利用多径信号参数aoa,toa及基站、目标和散射体间的位置关系建立定位方程组，此定位方程组为非线性欠定方程组。

设信号经第i个散射体反射后到达基站，经反射后的路径(即散射体和基站间的反射路径)长度为ri，则散射体坐标(xi,yi)为：

建立关于目标位置(xms,yms)和ri的非线性欠定方程组，如公式(2)所示：

其中：c表示光速，取c＝3×10⁸m/s。方程组中共有m个方程，m+2个变量。

公式(2)可简写为：

步骤四：引入宏蜂窝的圆环模型。假设散射体均匀分布在以ms为圆心的圆上，圆半径为r，ms、bs距离为d，如图2所示。本发明中取r＝100米，d＝1000米。

步骤五：根据圆环模型增加关于目标位置和散射体位置的非线性方程，将欠定方程组转化为超定方程组。

圆半径r可表示为：因此有：

cτ1-r1＝cτj+1-rj+1,j＝1,...,m-1(4)；

联合公式(3)，可以得到宏蜂窝的圆环模型下求解目标位置和散射体位置的超定方程组，如下式所示：

方程组共2m-1个方程，m+2个变量，当m>3时，式(5)即为超定方程组。

步骤六：实际测量时，由于aoa,toa测试误差，公式(5)中等号不总是成立，因此公式(5)可以写为：

其中x＝(xms,yms,ri)^t，χ＝(r1,rj+1)^t，εi(x)和ξj(χ)为残差。

故对目标的位置估计可以通过最小化如下目标函数得到：

其中

同时待估目标位置(xms,yms)应处于以最小传播距离lmin＝cτmin为半径，bs为圆心的圆内，同时最大最小到达角将ms限制在圆心角为α的扇形区域(即图3所示的扇形区域

内)，α为扇形圆心角，扇形的两条半径用向量表示为和长为lmin＝cτmin，α1为与bs和ms所在直线的夹角，α2为与bs和ms所在直线的夹角，α1,α2可以用方向矢量表示为：

其中：

α1(xms,yms)表示α1是关于xms,yms的函数，α2(xms,yms)表示α2是关于xms,yms的函数。

因此，(xms,yms)应处在下列不等式的交集内，即可行域z内，(xms,yms)所在扇形区域可以表示为下列不等式的交集：

其中xbs为基站的横坐标，ybs为基站的纵坐标，同时ri≤cτmax。

综上，式(5)的位置解算可以转换为求解如下的非线性约束最小二乘问题

步骤七：本发明选择经典的lm算法求解上述非线性约束最小二乘问题。公式(11)中的约束条件所述区域即是变量的可行域，记为z，如下式所示：

lm算法流程如图4所示：

步骤七-一：选取可行域内初始点x∈intz，给定初始的尺度因子参数μ>0，增长因子β>0，终止误差0≤ε＜＜1，令k＝1。

步骤七-二：计算雅克比矩阵j(xk)和hesse矩阵hk＝jk^tjk(这里为方便表达将h(xk)，j(xk)简写为hk，jk)；其中jk可以表示为如下形式：

步骤七-三：计算搜索方向dk＝-(hk+μi)^-1jk^tfk，其中ι为单位矩阵，fk＝(f1,f2,...,f2m-1)^t。

步骤七-四：由armijio线搜索算法求步长λk，使之满足

并令xk+1＝xk+λkdk得到下一时刻变量值。

步骤七-五：计算f(xk+1)，如果f(xk+1)<f(xk)，转入步骤七-六，否则转步入骤七-七。

步骤七-六：如果||dk||2≤ε，得到停止迭代，否则令μ:＝μ/β，同时k:＝k+1，重复步骤七-二至步骤七-五。

步骤七-七：如果||dk||2≤ε，得到停止迭代，否则令μ:＝μβ，同时k:＝k+1，重复步骤七-二至步骤七-五。

至此，完成了步骤七的全部操作。

步骤八：当残差较大时，最小二乘解算结果误差较大，主要原因是令hesse矩阵h＝j^tj时，忽略了其二阶项但是二阶偏导的计算比较复杂，因此引入拟牛顿bfgs算法，用不包含二阶偏导的矩阵对hesse矩阵进行更准确的逼近，再利用lm算法求解，此算法可以有效的减小算法误差，提高定位精度。本发明称此改进算法为“lm+bfgs”算法，其流程如图5所示。

步骤八-一：选取可行域内初始点x∈intz，给定lm算法参数：尺度因子μ>0，增长因子β>0；给定bfgs算法参数：初始对称正定阵β1＝i，ι为单位矩阵；给定终止误差0≤ε＜＜1，令k＝1，本发明中μ＝0.01,β＝10,ε＝10^-3。

步骤八-二：计算公式(7)的一阶导将j(xk),f(xk)简写为jk,fk，jk为雅克比矩阵，fk＝(f1,f2,...,f2m-1)^t，jk表示为如下形式：

步骤八-三：计算搜索方向dk＝-(bk+μi)^-1gk，其中，bk为对f(x)进行泰勒展开的二阶项的近似矩阵，其初始值为β1。

步骤八-四：由armijio搜索求步长λk，令xk+1＝xk+λkdk，armijio搜索求步长流程如图6所示。

步骤八-四-一：给定armijio算法的参数ρ∈(0,1),σ∈(0,0.5)，最大迭代次数mmax，令m＝0，本发明中ρ＝0.5,σ＝0.24,mmax＝20。

步骤八-四-二：计算f(xk+ρ^mdk)和f(xk)+σρ^mgk^tdk，若f(xk+ρ^mdk)≤f(xk)+σρ^mgk^tdk，转步骤八-四-四，否则转步骤八-四-三。

步骤八-四-三：令m:＝m+1，并判断m是否达到最大迭代次数mmax，如否，则转入步骤八-四-二，若是，则停止迭代，并执行步骤八-四-四。

步骤八-四-四：输出搜索步长λk＝ρ^m。

至此，完成了步骤八-四的全部操作。

步骤八-五：计算下一时刻一阶导数计算pk＝xk+1-xk和qk＝gk+1-gk。

步骤八-六：根据bfgs算法更新βk+1，βk+1的更新表达式如下：

其中：bk为对f(x)进行泰勒展开的二阶项的近似矩阵，βk+1为bk更新后的值。

步骤八-七：计算f(xk+1)，如果f(xk+1)<f(xk)，转入步骤八-八，否则转入步骤八-九。

步骤八-八：如果||dk||2≤ε，得到停止迭代，否则令μ:＝μ/β，同时k:＝k+1，重复步骤八-二至步骤八-七。

步骤八-九：如果||dk||2≤ε，得到停止迭代，否则令μ:＝μβ，同时k:＝k+1，重复步骤八-二至步骤八-七。

至此，完成了步骤八的全部操作。

步骤九：算法迭代结束后，得到估计参数最后输出估计的目标位置