CN110632555B - 一种基于矩阵特征值扰动的tdoa直接定位方法 - Google Patents

Abstract

本发明公开一种基于矩阵特征值扰动的TDOA直接定位方法，包括：M个观测站同步接收目标信号源辐射的信号，并将接收的信号汇总到中心站；生成Sinc内插函数矩阵，通过Sinc内插函数矩阵建立关于TDOA参数的信号模型；利用关于TDOA参数的信号模型建立关于目标位置的最大似然的优化函数；基于矩阵特征值扰动方法得出所述最大似然的优化函数的梯度向量和Hessian矩阵；设置迭代条件，基于所述梯度向量和Hessian矩阵，利用Newton迭代方法寻优，得到最终的目标位置估计结果。本发明能够有效的提高低信噪比下的定位估计精度，具有较高的可靠性和实际应用价值。

Description

一种基于矩阵特征值扰动的TDOA直接定位方法

技术领域

本发明属于无线信号定位技术领域，涉及一种基于矩阵特征值扰动的TDOA直接定位方法，特别是针对基于TDOA的单目标无源定位场景。

背景技术

无线信号定位技术在民用和军事领域都有着十分重要的意义，因其在电子对抗中的巨大作用，是各个国家非常重视的研究项目。经过几十年的研究发展，无源定位方法从定位方式上被分为了两步定位(M.Wax,T.Kailath,"Decentralized processing in sensorarrays."IEEE Transactions on Acoustics Speech&Signal Processing.33.5(1985):1123–1129.)和直接定位(王云龙，吴瑛.联合时延与多普勒频率的直接定位改进算法[J].西安交通大学学报，2015，49(4)：123-129.)(Direct Position Determination,DPD)两种。传统的无源定位一般均采用两步定位方法，其原理是首先从接收信号中提取包含信号定位信息的相关参数(主要包含空域、时域、频域以及能量域等参数)，例如：到达角度(Angle ofArrival,AOA)、到达时间(Time of Arrival)、到达时间差(Time Difference of Arrival,TDOA)和多普勒频差(Difference of Doppler,DD)、接收信号强度(Received Signal ofStrength,RSS)以及多参数联合来估计目标的位置参数。但是，大量的研究实验表明，传统的两步定位方法存在难以达到渐进最优、多目标分辨能力较弱等数据关联问题。

近年来，直接定位技术受到了国内外学者的广泛研究。其原理是从接收信号数据中直接提取信号发射源的位置参数，这样做省去了传统两步定位技术中估计中间观测量的过程，避免了参数估计和目标位置解算分离而导致的定位信息损失，较于两步定位方式，能够获得更高的估计精度。

本发明基于TDOA设计了一种直接定位方法，并设计了相应的快速算法对目标位置进行解算，相比于传统两步定位方法，本方法有着更高的估计精度，且减少了算法的计算量，更易实现。

发明内容

本发明针对到达时间差的定位问题，提供了一种基于矩阵特征值扰动的TDOA直接定位方法，能够快速稳定高效的提供信号辐射源目标位置参数估计。

为了实现上述目的，本发明采用以下技术方案：

一种基于矩阵特征值扰动的TDOA直接定位方法，包括：

步骤1：M个观测站同步接收目标信号源辐射的信号，并将接收的信号汇总到中心站；

步骤2：生成Sinc内插函数矩阵，通过Sinc内插函数矩阵建立关于TDOA参数的信号模型；

步骤3：利用关于TDOA参数的信号模型建立关于目标位置的最大似然的优化函数；

步骤4：基于矩阵特征值扰动方法得出所述最大似然的优化函数的梯度向量和Hessian矩阵；

步骤5：设置迭代条件，基于所述梯度向量和Hessian矩阵，利用Newton迭代方法寻优，得到最终的目标位置估计结果。

进一步地，所述关于TDOA参数的信号模型为：

其中，

为第1到M观测站接收的信号原始数据构成的观测向量，

表示各观测站接收到的噪声向量，

为观测站接收到的衰减后的信号波形，s为目标信号源发射信号原始波形；A(θ)＝[(a₁(p))^T(b_2,1a₂(p))^T …(b_M,1a_M(p))^T]^T为系数矩阵，θ＝[p^T b^T]^T为待估计参数p和b的组合，p＝[x y]^T为待估计的目标位置坐标，b＝[b_2,1b_3,1...b_M,1]^T为衰减系数构成的向量，b_m,1＝b_m/b₁,m＝2,...,M为归一化后的衰减系数，b₁为参考站接收到的信号传输过程中的衰减系数，b_m为第m个观测站接收到的信号传输过程中的衰减系数，A(θ)的子矩阵表示为：

其中l_m(k)＝sinc(k+τ_m,1/T_s)为Sinc函数，a_m(p)为第m个观测站的Sinc内插函数矩阵，k∈[-K,K]，K为Sinc函数的阶数，T_s为信号采样周期，m∈[1,M]。

进一步地，所述目标位置的最大似然的优化函数为：

其中，λ_max为M×M阶的半正定Hermitian矩阵X_RR＝R_a ^TR_a的最大特征值，

为经过变换后包含目标位置参数p的复杂矩阵。

进一步地，所述最大似然的优化函数的梯度向量和Hessian矩阵分别为：

其中，h(p)为最大似然的优化函数的梯度向量，H(p)为最大似然的优化函数的Hessian矩阵，X_RR(p)为包含变量p的M×M阶半正定Hermitian矩阵，e_M为X_RR的第M个特征向量，H₁(p)和H₂(p)为中间变量，表达式为：

其中，e_m为X_RR(p)的第m个特征向量，m∈[1,M-1]，λ_m表示X_RR(p)的第m个特征值，λ_M表示X_RR(p)的第M个特征值。

进一步地，所述步骤5包括：

步骤5.1：利用两步最小二乘方法粗略估计目标位置参数，将目标位置的粗略估计结果作为最大似然的优化函数的迭代初值p⁽⁰⁾，设置迭代的步长因子μ⁽⁰⁾∈(0,1)和步长衰减因子η∈(0,1)，设置迭代的停止条件ε∈(0,1)和初始迭代次数i:＝0；

步骤5.2：判断是否达到迭代的停止条件，若否则对矩阵X_RR(p⁽ⁱ⁾)进行特征分解，获得特征向量

和对应的特征值

若是则终止迭代；

步骤5.3：若

则μ^(i-1):＝ημ^(i-1)，i:＝i-1并跳转至步骤5.5，否则进行步骤5.4；

步骤5.4：按照最大似然的优化函数的梯度向量和Hessian矩阵公式分别计算h(p⁽ⁱ⁾)和H(p⁽ⁱ⁾)；

步骤5.5：利用Newton迭代方法迭代计算p⁽ⁱ⁺¹⁾＝p⁽ⁱ⁾-μ⁽ⁱ⁾(H(p⁽ⁱ⁾))^-1h(p⁽ⁱ⁾)，迭代次数i:＝i+1，跳转至步骤5.2。

与现有技术相比，本发明具有的有益效果：

相比于传统的两步定位方法，本发明所提方法能够克服信号源与数据不匹配的问题，针对基于TDOA定位的现实场景，算法能够直接利用各观测站接收到的信号原始数据，基于最大似然准则，直接估计得到目标的精确位置参数，较两步定位方法，定位精度更高。同时，本发明利用基于Hermitian矩阵扰动理论，利用Newton迭代的方法求解关于目标位置的最大似然的优化函数，收敛速度更快，能够有效减少计算量，具有更好的时效性。

附图说明

图1为本发明实施例一种基于矩阵特征值扰动的TDOA直接定位方法的基本流程图；

图2为本发明实施例一种基于矩阵特征值扰动的TDOA直接定位方法的两种目标定位场景示意图；其中(2a)为近场定位场景示意图，(2b)为远场定位场景示意图；

图3为本发明实施例一种基于矩阵特征值扰动的TDOA直接定位方法的结果对比示意图；其中(3a)为近场定位场景中目标位置估计均方根误差随着信噪比的变化曲线对比示意图，(3b)为远场定位场景中目标位置估计均方根误差随着信噪比的变化曲线对比示意图。

具体实施方式

下面结合附图和具体的实施例对本发明做进一步的解释说明：

如图1所示，一种基于矩阵特征值扰动的TDOA直接定位方法，包括：

步骤S101：M个观测站同步接收目标信号源辐射的信号，并将接收的信号汇总到中心站；

步骤S102：生成Sinc内插函数矩阵，通过Sinc内插函数矩阵建立关于TDOA参数的信号模型；

步骤S103：利用关于TDOA参数的信号模型建立关于目标位置的最大似然的优化函数；

步骤S104：基于矩阵特征值扰动方法得出所述最大似然的优化函数的梯度向量和Hessian矩阵；

步骤S105：设置迭代条件，基于所述梯度向量和Hessian矩阵，利用Newton迭代方法寻优，得到最终的目标位置估计结果。

具体地，步骤S101中，各观测站接收到的信号模型可表示为：

r_m(t)＝b_ms(t-τ_m)+n_m(t),0<t<T,m＝1,2,...,M

其中τ_m＝||p-u_m||/c，p＝[x y]^T为待估计目标的位置坐标，u_m＝[x_m y_m]^T表示第m个观测站的位置坐标，b_m为第m个观测站接收到的信号传输过程中的衰减系数，n_m(t)表示第m个观测站在t时刻的噪声，s(t-τ_m)表示经过时延τ_m的信号波形。这里假设信号的观测时间段为(0,T)，在这个时间段足够长且各站能够在各自观测时段内捕捉到同一段发射信号。选定其中一个站为参考站(中心站)，那么第m个观测站和参考站的接收信号到达时间差(TDOA参数)可表示为：

τ_m,1＝τ_m-τ₁＝(||p-u_m||-||p-u₁||)/c,m＝1,2,...,M

其中，τ_m＝||p-u_m||/c为信号从发射源到第m个观测站的传播时间，u_m＝[x_m y_m]^T表示第m个观测站的位置坐标，c为信号传输速率。

将参考站(第1个观测站)接收到的信号作为参考信号，各个观测站接收到的信号可以重新表示为：

其中

为观测站接收到的传播τ_m衰减后的信号波形，同时，对衰减系数归一化b_m,1＝b_m/b₁，b_m和b₁分别为第m个参考站(第1个观测站)接收到的信号传输过程中的衰减系数。

具体地，步骤S102中，通过引入Sinc内插函数矩阵，接收信号模型可以转化为：

其中v表示不定积分的积分变量，sinc(v-τ_m,1)表示Sinc函数，

表示经过时延v的波形函数

为观测站接收到的随时间t衰减后的信号波形，s为目标信号源发射信号原始波形。

以T_s为信号采样周期，信号表示为

其中r_m(nT_s)为第m个观测站第n个采样时刻的接收信号，n_m(nT_s)表示第m个观测站第n个采样时刻的噪声。

由于Sinc函数的阶数有限，所以有

其中k∈[-K,K]，K为Sinc函数的阶数，T_s为信号采样周期，m∈[1,M]。

而后利用包含信号源位置参数的矩阵A(p,b)构建观测方程

其中，

为第1到M观测站接收的信号原始数据构成的观测向量，

表示各观测站接收到的噪声向量，

为观测站接收到的衰减后的信号波形，s为目标信号源发射信号原始波形；A(θ)＝[(a₁(p))^T(b_2,1a₂(p))^T…(b_M,1a_M(p))^T]^T为系数矩阵，θ＝[p^T b^T]^T为待估计参数p和b的组合，p＝[x y]^T为待估计的目标位置坐标，b＝[b_2,1b_3,1...b_M,1]^T为衰减系数构成的向量，b_m,1＝b_m/b₁,m＝2,...,M为归一化后的衰减系数，b₁为参考站接收到的信号传输过程中的衰减系数，b_m为第m个观测站接收到的信号传输过程中的衰减系数A(θ)的子矩阵表示为：

其中l_m(k)＝sinc(k+τ_m,1/T_s)为Sinc函数，a_m(p)为第m个观测站的Sinc内插函数矩阵，k∈[-K,K]，K为Sinc函数的阶数，T_s为信号采样周期，m∈[1,M]。

具体地，步骤S103中，基于最大似然准则建立关于目标位置p的最大似然函数的优化函数并化简得到：

其中，lnL(ρ)为对数似然函数化简的结果，L(ρ)表示最大似然目标函数，

表示噪声功率，Const表示常数，可以省略常数项，ρ为

p、b的集合。

首先利用最小二乘估计信号波形估计结果

再将信号波形估计结果

带入最大似然优化函数(令

)

I_MN表示MN×MN阶单位矩阵。

可得θ估计结果θ_opt：

为了降低矩阵维度，减少计算量，由广义特征值的相关性质，关于目标函数(目标位置的最大似然的优化函数)可以改写为

其中λ_max为M×M阶的半正定Hermitian矩阵X_RR＝R_a ^TR_a的最大特征值，

R_a为经过变换后包含目标位置参数p的复杂矩阵。

具体地，步骤S104中，为了避免传统网格搜索带来的巨大计算量，基于矩阵特征值扰动理论，推导目标函数的梯度向量和Hessian矩阵用于后续Newton迭代。

由目标函数的闭式解可得，X_RR＝R_a ^TR_a是一个M×M阶的半正定Hermitian矩阵，假设该矩阵的M个特征值和其对应的特征向量分别为λ₁≤λ₂≤…≤λ_M和e₁,e₂,…,e_M，令该矩阵受到扰动后得到的半正定Hermitian矩阵为

矩阵扰动量为

扰动矩阵

的M个特征值分别为

并有如下关系式

其中，

下面基于此，推导目标函数的梯度向量和Hessian矩阵，由上述的推导，关于扰动矩阵

的特征值可表示为

表示在P点处的二阶无穷小量；

梯度向量h(p)和Hessian矩阵H(p)分别为

其中，h(p)为最大似然的优化函数的梯度向量，H(p)为最大似然的优化函数的Hessian矩阵，X_RR(p)为包含变量p的M×M阶半正定Hermitian矩阵，e_M为X_RR的第M个特征向量，H₁(p)和H₂(p)为中间变量，表达式为：

具体地，步骤S105中，根据步骤S104中得到的目标函数的梯度向量和Hessian矩阵，利用Newton迭代对目标函数进行优化，得到目标位置的精确估计结果。

所设计出的Newton迭代方法的实现步骤为：

步骤S1051：利用两步最小二乘方法，粗略估计出目标位置参数，将目标位置的粗略估计结果作为迭代的初值p⁽⁰⁾，设置迭代的步长因子μ⁽⁰⁾∈(0,1)和步长衰减因子η∈(0,1)，再设置迭代的停止条件ε和迭代次数i的初始值0，即i:＝0，进行步骤S1052，ε∈(0,1)；作为一种可实施方式，ε＝10^-5；

步骤S1052：判断是否达到迭代的停止条件，若否则对矩阵X_RR(p⁽ⁱ⁾)进行特征分解，获得特征向量

和对应的特征值

进行步骤S1053；若是则终止迭代；

步骤S1053：若

则μ^(i-1):＝ημ^(i-1)，i:＝i-1并跳转至步骤S1055，否则进行步骤S1054；即如果迭代后结果没有更好，缩短迭代的步长因子，返回上一步再次迭代更新，计数因子i-1，用新的步长因子重新迭代；

步骤S1054：利用步骤S1055中推导的梯度向量和Hessian矩阵的代数表达式，分别计算第i次迭代对应的梯度向量h(p⁽ⁱ⁾)和Hessian矩阵H(p⁽ⁱ⁾)，进行步骤S1055；

步骤S1055：利用Newton迭代方法迭代计算p⁽ⁱ⁺¹⁾＝p⁽ⁱ⁾-μ⁽ⁱ⁾(H(p⁽ⁱ⁾))^-1h(p⁽ⁱ⁾)，迭代次数i:＝i+1，跳转至步骤S1052。

为验证本发明效果，进行如下实验：

如图2中(2a)和(2b)所示，假设有四个观测站对目标进行定位，观测站的位置坐标分别为(0m,0m)，(-4376m,2540m)，(4376m,2540m)，(0m,-4250m)，且每个观测站都可以在观测时段内接收到信号的直达波。针对远场和近场两种定位场景，分别设置了目标位置为(400m,800m)(近场)和(10000m,-10000m)(远场)。

下面将本专利公开的一种基于矩阵特征值扰动的TDOA直接定位方法同现有的基于TDOA的两步定位方法进行比较。比较不同方法的目标位置估计结果的均方根误差随着信噪比(dB)的变化曲线，如图3中(3a)和(3b)所示，分别给出了近场和远场两种情况下的估计结果，在低信噪比下，本专利公开的目标直接定位方法估计结果明显好于传统的两步定位方法，且在高信噪比条件下，定位精度能够进一步提升。

以上所示仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (1)

1.一种基于矩阵特征值扰动的TDOA直接定位方法，其特征在于，包括：

步骤1：M个观测站同步接收目标信号源辐射的信号，并将接收的信号汇总到中心站；

步骤2：生成Sinc内插函数矩阵，通过Sinc内插函数矩阵建立关于TDOA参数的信号模型；

所述关于TDOA参数的信号模型为：

其中，

为第1到M观测站接收的信号原始数据构成的观测向量，

表示各观测站接收到的噪声向量，

为观测站接收到的衰减后的信号波形，s为目标信号源发射信号原始波形；A(θ)＝[(a₁(p))^T (b_2,1a₂(p))^T…(b_M,1a_M(p))^T]^T为系数矩阵，θ＝[p^T b^T]^T为待估计参数p和b的组合，p＝[x y]^T为待估计的目标位置坐标，b＝[b_2,1 b_3,1…b_M,1]^T为衰减系数构成的向量，b_m,1＝b_m/b₁,m＝2,...,M为归一化后的衰减系数，b₁为参考站接收到的信号传输过程中的衰减系数，b_m为第m个观测站接收到的信号传输过程中的衰减系数，A(θ)的子矩阵表示为：

其中l_m(k)＝sinc(k+τ_m,1/T_s)为Sinc函数，a_m(p)为第m个观测站的Sinc内插函数矩阵，k∈[-K,K]，K为Sinc函数的阶数，T_s为信号采样周期，m∈[1,M]；

步骤3：利用关于TDOA参数的信号模型建立关于目标位置的最大似然的优化函数；

所述目标位置的最大似然的优化函数为：

其中，λ_max为M×M阶的半正定Hermitian矩阵X_RR＝R_a ^TR_a的最大特征值，

为经过变换后包含目标位置参数p的复杂矩阵；

步骤4：基于矩阵特征值扰动方法得出所述最大似然的优化函数的梯度向量和Hessian矩阵；

所述最大似然的优化函数的梯度向量和Hessian矩阵分别为：

其中，h(p)为最大似然的优化函数的梯度向量，H(p)为最大似然的优化函数的Hessian矩阵，X_RR(p)为包含变量p的M×M阶半正定Hermitian矩阵，e_M为X_RR的第M个特征向量，H₁(p)和H₂(p)为中间变量，表达式为：

其中，e_m为X_RR(p)的第m个特征向量，m∈[1,M-1]，λ_m表示X_RR(p)的第m个特征值，λ_M表示X_RR(p)的第M个特征值；

步骤5：设置迭代条件，基于所述梯度向量和Hessian矩阵，利用Newton迭代方法寻优，得到最终的目标位置估计结果；

所述步骤5包括：

步骤5.1：利用两步最小二乘方法粗略估计目标位置参数，将目标位置的粗略估计结果作为最大似然的优化函数的迭代初值p⁽⁰⁾，设置迭代的步长因子μ⁽⁰⁾∈(0,1)和步长衰减因子η∈(0,1)，设置迭代的停止条件ε∈(0,1)和初始迭代次数i:＝0；

步骤5.2：判断是否达到迭代的停止条件，若否则对矩阵X_RR(p⁽ⁱ⁾)进行特征分解，获得特征向量

和对应的特征值

若是则终止迭代；

步骤5.3：若

则μ^(i-1):＝ημ^(i-1)，i:＝i-1并跳转至步骤5.5，否则进行步骤5.4；

步骤5.4：按照最大似然的优化函数的梯度向量和Hessian矩阵公式分别计算h(p⁽ⁱ⁾)和H(p⁽ⁱ⁾)；

步骤5.5：利用Newton迭代方法迭代计算p⁽ⁱ⁺¹⁾＝p⁽ⁱ⁾-μ⁽ⁱ⁾(H(p⁽ⁱ⁾))^-1h(p⁽ⁱ⁾)，迭代次数i:＝i+1，跳转至步骤5.2。

Images