基于TDOA/GROA的多站连续定位模型的制作方法

本发明属于多传感器无源定位技术，涉及利用到达增益比等量测信息对辐射源目标的定位问题，提供了一种基于到达增益比的多站连续定位模型。
背景技术：
：无源定位技术一直是信号处理领域的重要研究课题，具有隐蔽性好、探测距离较远、适用性强等优点，在雷达，无线传感器网络等领域受到国内外学者的广泛关注。目前主要使用的量测信息有方位角、到达时差(TDOA)、到达频差以及到达增益比(GROA)等，或以上量测组合得到的多种联合定位系统。传统的基于方位角的辐射源定位方法，需要在不同的角度对同一信号进行检测，并且容易因为坐标系不同和观测站姿态引入新的误差项。而基于到达时差和到达增益比量测的无源定位技术无需测量方位角、俯仰角等信息，因此不会受到坐标系转换和传感器载体姿态误差的影响，具有稳定性强、精度高等优点。目前的时差定位算法中，多通过引入额外的中间变量(辐射源到参考观测站的径向距离)建立线性时差方程。在此定位模型基础上，利用两步加权最小二乘(TS-WLS)算法、约束总体最小二乘(CTLS)算法、束加权最小二乘(CWLS)算法等进一步提高了定位精度。另一方面，由于辐射源信号会随着传播而产生能量衰减，将达到增益比量测信息加入到TDOA体制中，构造TDOA/GROA联合定位体制，并利用GROA量测对中间变量的修正，改善对辐射源的定位精度。以上都是使用的基于中间变量的定位模型的定位算法，对于多运动站的情况下，由于观测站不断运动，中间变量即辐射源到参考观测站的径向距离也在不断变化，使得上述定位模型及其算法无法直接推广到多运动站连续定位，需要建立新的定位模型。技术实现要素：本发明的目的在于针对固定辐射源目标，提供一种基于TDOA/GROA的多站连续定位模型，并利用CWLS算法求解辐射源位置信息。考虑到中间变量造成的局限性，从TDOA和GROA的量测方程出发，得到了只与目标状态向量有关的定位模型，从而避免了引入中间变量。然后根据量测模型推导了定位模型的误差项，得到了一种基于CWLS的TDOA/GROA联合定位算法。该算法适用于多运动站场景下的固定辐射源目标连续定位，其方法流程如图1所示，包括以下步骤：1问题描述假设定位场景中有M个运动观测站接收固定辐射源的信号，通过测量两个不同站之间的达到时间差和到达增益比来估计辐射源位置。N维空间时，由文献(15]可知得到闭式解的条件是M≥N+2。本文以2维空间为研究背景，所采用的方法可推广到3维空间中。记k时刻未知辐射源的位置向量为s0＝[x0(k)y0(k)]T，第i个观测站位置向量为si＝[xi(k)yi(k)]T(i＝1,2,…,M)。以第1个观测站为参考站，则真实时差量测为t‾i1=r‾i1c,i=1,2,...,M---(1)]]>式中：为真实距离差量测，ri＝|si-s0|(i＝1,2,…,M)表示第i个观测站到辐射源目标的径向距离，c为电磁波传播速度。典型定位场景示意图如图2所示。根据信号传播理论，信号的传播损耗与信号源和传感器之间距离的n次方成正比，在自由空间中，不考虑多径效应时，可将损耗因子n设定为常数2。同样以第1个观测站为参考站，则真实增益比量测为g‾i1=(rir1)2=|si-s0|2|s1-s0|2,i=2,3,...,M---(2)]]>2TDOA/GROA定位模型(1)TODA定位模型基于TDOA量测方程，为表示方便，省略时间标签，该模型的向量形式为1t21(|s2|2-|s1|2)-1t(i+2)1(|s(i+2)|2-|s1|2)+c2t(i+2)2=[2t21(s2T-s1T)-2t(i+2)1(s(i+2)T-s1T)]s0(i=1,2,...,M-2)---(3)]]>对式(3)做变换，可得无测量误差时的TDOA线性方程为A‾t(k)x=b‾t(k)---(4)]]>式中：A‾t(k)=[αt,1T(k),αt,2T(k),...,αt,M-2T(k)]Tαt,i(k)=[2r‾(i+2)1(k)x21(k)-2r21(k)x(i+2)1(k)2r‾(i+2)1(k)y21(k)-2r‾21(k)y(i+2)1(k)]b‾t(k)=[βt,1(k),βt,2(k),...,βt,(M-2)(k)]Tβt,i(k)=r‾(i+2)1(k)(d22(k)-d12(k))-r‾21(k)(di+22(k)-d12(k))+(r‾(i+2)1(k))2r‾21(k)-(r‾21(k))2r‾(i+2)1(k)xij(k)=xi(k)-xj(k),yij(k)=yi(k)-yj(k),(i,j=1,2,...,M;i≠j)---(5)]]>实际到达距离差(RDOA)测量值中含有测量噪声，则实际RDOA量测模型为ri1(k)=r‾i1(k)+δri1(k),i=2,3,...,M---(6)]]>式中：δgi1为噪声引起的RDOA测量误差，服从零均值的高斯分布。记到RDOA误差向量为δr(k)＝[δr21(k),δr31(k),…,δrM1(k)]T，其量测噪声的协方差矩阵为E[δr(k)δr(k)T]=Rrr=c2▿rr-1---(7)]]>▿rr=T3πω3ρ21+Mρ(MIM-1-1(M-1)×(M-1))---(8)]]>式中：T为观测时间，ω为信号带宽，ρ＝Psignal/Pnoise为信号的信噪比，Psignal和Pnoise分别为信号和噪声的功率谱密度函数，IM-1为M-1阶单位矩阵，1(M-1)×(M-1)为M-1阶全1矩阵。(2)GRDA定位模型注意到对式(2)进行等式变换，省略时间标签可得gi1|s1|2-|si|2-2(gi1s1T-siT)s0gi1-1=-|s0|2,i=2,3,...,M---(9)]]>令i＝2，上式可写作g21|s1|2-|s2|2-2(g21s1T-s2T)s0g21-1=-|s0|2---(10)]]>式(9)和式(10)中均含有未知项-|s0|2，两两相减消去，整理后得到基于GROA定位模型的向量形式为(gi1-1)(g21|s1|2-|s2|2)-(g21-1)(gi1|s1|2-|si|2)=[2(gi1-1)(g21s1T-s2T)-2(g21-1)(gi1s1T-siT)]s0(i=3,4,...,M)---(11)]]>这样我们就得到了无测量误差时关于辐射源目标状态的GROA线性方程A‾g(k)x=b‾g(k)---(12)]]>式中：A‾g(k)=[α‾g,1T(k),α‾g,2T(k),...,α‾g,(M-2)T(k)]Tα‾g,i(k)=2(g‾(i+2)1(k)-1)(g‾21(k)x1(k)-x2(k))2(g‾(i+2)1(k)-1)(g‾21(k)y1(k)-y2(k))Tb‾g(k)=[β‾g,1T(k),β‾g,2T(k),...,β‾g,(M-2)T(k)]Tβ‾g,i(k)=(g‾(i+2)1(k)-1)(g‾21(k)r12(k)-r22(k))-(g‾21(k)-1)(g‾(i+2)1(k)r12(k)-ri+22(k))---(13)]]>实际到达增益比测量值中含有测量噪声δgi1(k)，则实际GROA量测模型为gi1(k)=g‾i1(k)+δgi1(k),i=2,3,...,M---(14)]]>式中：δgi1为噪声引起的到达增益比测量误差，服从零均值的高斯分布。记到达增益比误差向量为δg(k)＝[δg21(k),δg31(k),…,δgM1(k)]T，其协方差矩阵为E[δg(k)δg(k)T]=Rgg=▿gg-1---(15)]]>▿gg=Tπωρ21+Mρ(MIM-1+1-Mρ1+Mρ1(M-1)×(M-1))---(16)]]>原有需要中间变量的GROA定位模型通过联合TDOA和GROA量测得到ri1＝(gi1-1)r1，可知原定位模型通过加入GROA量测信息提高对中间变量r1的估计精度，进而提高算法性能，并没有通过GROA量测直接改善目标定位精度。由式(3)和式(11)可知所提无需中间变量的GROA定位模型能够直接估计目标状态，并避免引入中间变量，适用于多运动站的无源定位。综上所述，通过引入GROA定位模型，构建TDOA/GROA联合定位模型，可以充分利用GROA中包含的目标位置信息，提高算法性能，并且保持了无需中间变量的优点，适用于多运动站对固定目标的连续定位。其联合定位模型为A‾(k)x=b‾(k)---(17)]]>式中：3基于CWLS的定位算法将式(6)和式(14)带入到和中，忽略高阶误差项后，整理可得A(k)=A‾(k)+ΔA(k)b(k)=b‾(k)+Δb(k)---(18)]]>式中：A(k)和b(k)分别是带入实际测量值的观测矩阵和观测向量，且△A(k)＝[-B1(k)v(k)-B2(k)v(k)]，△b(k)＝B3(k)v(k)，忽略时间标签，误差项推导结果如式(19)所示。B1=blkdiag([b1-2x21IM-2],[b1-2x21IM-2]),b1=[2x31,2x41,...,2xM1]TB2=blkdiag([b2-2y21IM-2],[-b12y21IM-2]),b2=[2y31,2y41,...,2yM1]TB3=blkdiag([bt1,bt2],[bg,(r12-r22)IM-2])bt1=d12-d32+r‾312-2r‾21r‾31d12-d42+r‾412-2r‾21r‾41...d12-dM2+r‾M12-2r‾21r‾M1Tbt2=blkdiag([d22-d12-r‾212+2r31r21d22-d12-r‾212+2r41r21...d22-d12-r‾212+2rM1r21])bg=d32-d12d42-d12...dM2-d12Tv=δr21δr31...δrM1δg21δg31...δgM1T---(19)]]>式中：blkdaig(·)为分块对角矩阵。将式(15)带入式(13)，整理后可得b(k)-A(k)x＝△b(k)-△A(k)x＝C(k)v(k)(20)式中：C(k)＝B1(k)x0+B2(k)y0+B3(k)定义加权矩阵W(k)＝E(C(k)v(k)vT(k)CT(k))＝C(k)RCT(k)(21)式中协方差矩阵R＝E(v(k)vT(k))＝blkdiag(Rrr,Rgg)。定义增广矩阵Au＝[-Ab](22)和增广解向量X＝[x1]T(23)式中：A＝[AT(1)AT(2)…AT(k)]T，b的定义与A类似。易得Au=A‾u+ΔAu---(24)]]>ΔAu=[B1vB2vB3v]v=v1Tv2T...vkTTBi=blkdiag(Bi1,Bi2,...,Bik),i=1,2,3---(25)]]>从而得到CWLS算法的最小化代价函数为ϵ=XTAuTW-1AuX---(26)]]>式中：ε的定义与A类似；W＝blkdiag(W(1),W(2),…W(k))。带入式(24)并对ε求期望，可得平均代价函数为E(ϵ)=XTA‾uTW-1A‾uX+XTE(ΔAuTW-1Au)X---(27)]]>等式右边两项均为非负且与X相关。第一项为理想代价函数，当时，可取得最小值0。正是第二项的存在，当最小化E(ε)时，使解偏离真值从而产生偏差。因此，当约束第二项为常数时，E(ε)将在X取得真值时达到最小值。于是可构造CWLS问题为minXXTAuTW-1AuXs.t.XTΩX=const---(28)]]>式中：const为任意非负常数。利用拉格朗日乘子法可以解上述约束最小化问题。记λmin为(AuTW-1Au,Ω)]]>最小的广义特征值，λmin对应的广义特征值向量为，则辐射源的状态估计为x^=X^(1:2)/X^(3)---(29)]]>下面求解约束矩阵Ω，带入△Au＝[B1vB2vB3v]整理后可得Ω=Ω11Ω12Ω13Ω21Ω22Ω23Ω31Ω32Ω33---(30)]]>式中：注意到随着时间的推进，矩阵维数快速增大，导致计算量剧增，可利用式(31)和式(32)避免高维矩阵直接相乘，从而减小计算量，也可以对该算法加窗处理，提高算法的实时性。AuTW-1Au=Σi=1kAuT(i)W-1(i)Au(i)---(31)]]>Ωmn=Σi=1ktr(BmT(i)W-1(i)Bn(i)R)---(32)]]>注意到加权矩阵W与未知的辐射源的状态有关，计算时首先将W设为单位矩阵，得到辐射源的状态估计后，重新计算W，进而得到最终的目标状态估计。计算时，用实际测量值代替W和约束矩阵中的真实距离差和真实达到增益比。附图说明图1：技术方案流程图；图2：典型定位场景。具体实施方式以下结合说明书附图对本发明作进一步详细描述。参照说明书附图1，本发明中基于TDOA/GROA的多站连续定位方法分为以下几个步骤：步骤1：令k＝1，初始化协方差矩阵RR＝blkdiag(Rrr,Rgg)式中Rrr和Rgg分别为RDOA量测和GROA量测的协方差矩阵。步骤2：根据TDOA和GORA量测信息以及站址信息计算定位模型Au(k)Au(k)＝[-A(k)b(k)]A(k)=AtT(k)AgT(k)Tb(k)=btT(k)bgT(k)T]]>At(k)=[αt,1T(k),αt,2T(k),...,αt,M-2T(k)]Tαt,i(k)=[2r(i+2)1(k)x21(k)-2r21(k)x(i+2)1(k)2r(i+2)1(k)y21(k)-2r21(k)y(i+2)1(k)]bt(k)=[βt,1(k),βt,2(k),...,βt,(M-2)(k)]Tβt,i(k)=r(i+2)1(k)(d22(k)-d12(k))-r21(k)(di+22(k)-d12(k))+(r(i+2)1(k))2r21(k)-(r21(k))2r(i+2)1(k)xij(k)=xi(k)-xj(k),yij(k)=yi(k)-yj(k),(i,j=1,2,...,M;i≠j)]]>Ag(k)=[αg,1T(k),αg,2T(k),...,αg,(M-2)T(k)]Tαg,i(k)=2(g(i+2)1(k)-1)(g21(k)x1(k)-x2(k))2(g(i+2)1(k)-1)(g21(k)y1(k)-y2(k))Tbg(k)=[βg,1T(k),βg,2T(k),...,βg,(M-2)T(k)]Tβg,i(k)=(g(i+2)1(k)-1)(g21(k)r12(k)-r22(k))-(g21(k)-1)(g(i+2)1(k)r12(k)-ri+22(k))]]>其中ri1(k)和gi1(k)分别为k时刻的RODA和GROA的实际测量值且i＝2,3,…,M，xi和yi为k时刻观测站位置且i＝1,2,…,M；省略时间标签k，误差项B1、B2和B3为B1=blkdiag([b1-2x21IM-2],[b1-2x21IM-2]),b1=[2x31,2x41,...,2xM1]TB2=blkdiag([b2-2y21IM-2],[-b22y21IM-2]),b2=[2y31,2y41,...,2yM1]TB3=blkdiag([bt1,bt2],[bg,(d12-d22)IM-2])bt1=d12-d32+r‾312-2r‾21r‾31d12-d42+r‾412-2r‾21r‾41...d12-dM2+r‾M12-2r‾21r‾M1Tbt2=blkdiag([d22-d12-r‾212+2r31r21d22-d12-r‾212+2r41r21...d22-d12-r‾212+2rM1r21])bg=d32-d12d42-d12...dM2-d12T]]>其中：blkdaig(·)为分块对角矩阵，IM-2为M-2阶单位矩阵，且i＝1,2,…,M。步骤3：若已经得到目标状态估计则加权矩阵W(k)为W(k)=E(C(k)v(k)vT(k)CT(k))=C(k)RCT(k)C(k)=B1(k)x^0+B2(k)y^0+B3(k)]]>若目标状态估计未知，令加权矩阵W(k)为单位阵；步骤4：计算和ΩAuTW-1Au=Σi=1kAuT(i)W-1(i)Au(i)]]>Ω=Ω11Ω12Ω13Ω21Ω22Ω23Ω31Ω32Ω33]]>Ωmn=Σi=1ktr(BmT(i)W-1(i)Bn(i)R)]]>并通过对其广义特征值分解，得到最小广义特征值对应的广义特征值向量；步骤5：更新目标状态估计x^=X^(1:2)/X^(3)]]>令k＝k+1，并从步骤2继续运行。当前第1页1 2 3