0 概述

对信号源进行探测、定位和跟踪在信号处理应用中具有重要价值。相对于有源定位系统, 无源定位系统具有不主动向外发射信号、隐蔽性强、生存能力好、探测距离远等优点, 因此, 越来越受到国内外相关学者和工程技术人员的广泛关注^[1-3]。到达时间差(Time Difference of Arrival, TDOA)定位技术作为无源定位的一种关键技术^[4], 引起了研究人员的极大重视, 并已成功应用于无线蜂窝网络的移动用户定位问题中^[5-6]。

TDOA测向交叉是常用的定位方法, 其实质是一种解析的方法, 即通过空间中的几何关系计算出目标位置。该方法虽然计算简单, 但是其定位精度并不高^[7]。最小二乘(Least Squares, LS)是一种常见的定位算法^[8]。该方法应用范围较广, 可用于线性阵列、平面阵列和三维空间中的均匀分布阵列等场合, 但是经典的LS法只考虑观测向量的误差, 假设系数矩阵没有误差或不考虑系数矩阵的误差。文献[9]提出的泰勒级数展开法也是一种常用的定位算法。该算法定位精度较高、鲁棒性强, 但是泰勒级数展开法对初始值的要求极高, 如果初始值位置偏离真实值, 则迭代可能不收敛以至于无法得到结果^[10-11]。文献[12]利用最小二乘估计值作为Taylor级数展开的初值, 提高了算法的测向精度。

在野外起伏地形部署传感器网络会导致TDOA出现偏测量问题^[13]。文献[14]认为, 即使传感器位置的随机误差很小, 也会大幅降低目标的测向精度。在这种情况下, LS估计值将不再适合作为Taylor级数展开算法的初始估计值, 因为Taylor级数展开算法对初始值的依赖性非常强。

为弥补LS法不考虑系统误差的理论缺陷, 本文提出总体最小二乘(Total Least Squares, TLS)波达方向估计方法。基于测向模型建立测向方程, 并将非线性方程转化为伪线性方程。同时构造一个增广矩阵, 对其进行奇异值分解后计算出目标的位置, 其无需获得信号源位置的粗略估计和进行迭代计算^[15]。在此基础上, 对TLS方法的克拉美罗下界(Cramer-Rao Lower Bound, CRLB)^[16]进行分析。

1 TLS波达方向估计方法 1.1 最小二乘测向

图 1为基于TDOA的短基线传感器网络定向示意图, 假设传感器网络由N个传感器构成, 在图 1中用△表示, S_i为第i个传感器, 其坐标为(x_i, y_i, z_i), 接收到的时间为t_i, 各传感器之间的距离大约为30 m, 高度差约为5 m。取距离信号源最近的传感器作为基准传感器, 即x₀=y₀=z₀=0。M为瞬时声源的波达方向上的一点, M到基准传感器的距离为2 km~3 km, M在XOY平面上的投影记为M′。

将基准传感器放在坐标原点, 令S=[(x₁, y₁, z₁); (x₂, y₂, z₂); …; (x_N-1, y_N-1, z_N-1)], S为除基准传感器外, 其余传感器与基准传感器之间的距离向量构成的矩阵^[17]。时间差向量τ=[τ₁, τ₂, …, τ_N-1]^T, 其中τ_i=t_i-t₀。S₀M与S₀M′的夹角为俯仰角θ, S₀M′与X轴正向的夹角为方位角φ, 令γ=[φ, θ]^T。

由图 1可知, S₀P的长度cτ_i等于S_i在S₀M上的投影, 即有如下TDOA方程:

$ c{\tau _i} = {\mathit{\boldsymbol{S}}_i}\mathit{\boldsymbol{K}} $

(1)

其中, K为声源信号的单位方向向量。

$ \mathit{\boldsymbol{K}} = \left[ {\begin{array}{*{20}{c}} {{k_x}}\\ {{k_y}}\\ {{k_z}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\cos \theta \cos \varphi }\\ {\cos \theta \sin \varphi }\\ {\sin \theta } \end{array}} \right] $

(2)

整理可得:

$ c{\tau _i} = {x_i}\cos \theta \cos \varphi + {y_i}\cos \theta \sin \varphi + {z_i}\sin \theta $

(3)

将N－1个TDOA方程写成如下形式:

$ c\mathit{\boldsymbol{\tau }} = \mathit{\boldsymbol{SK}} $

(4)

由于真实TDOA值是无法获得的, 因此得到的测量值都是带噪声的。如果n=[n₁   n₂  …   n_N－1]∈ $\mathbb{R}$^(N-1)×1表示TDOA测量噪声矢量, 则有如下公式:

$ \mathit{\boldsymbol{\hat \tau }} = \mathit{\boldsymbol{\tau }} + \mathit{\boldsymbol{n}} = \frac{{\mathit{\boldsymbol{SK}}}}{c} + \mathit{\boldsymbol{n}} $

(5)

通常TDOA的测量误差是独立不相关的, 其协方差矩阵如下:

$ \mathit{\boldsymbol{Cov}}\left( {\mathit{\boldsymbol{\hat \tau }}} \right) = E\left( {\mathit{\boldsymbol{\hat \tau }}{{\mathit{\boldsymbol{\hat \tau }}}^{\rm{T}}}} \right) - E\left( {\mathit{\boldsymbol{\hat \tau }}} \right)E\left( {{{\mathit{\boldsymbol{\hat \tau }}}^{\rm{T}}}} \right) $

(6)

其中, E(·)表示求期望值。

当传感器数目多于3个时, 方向矢量估计K逼近理论值的优化问题是超定的, 综合以上分析, 可以采用最小二乘方法进行求解, 其优化模型如下:

$ \mathit{\boldsymbol{\hat K}} = \mathop {\arg \min }\limits_\mathit{\boldsymbol{k}} \left\{ {{{\left[ {\mathit{\boldsymbol{\hat \tau }} - \frac{{\mathit{\boldsymbol{SK}}}}{c}} \right]}^{\rm{T}}}\mathit{\boldsymbol{Cov}}{{\left( {\mathit{\boldsymbol{\hat \tau }}} \right)}^{ - 1}}\left[ {\mathit{\boldsymbol{\hat \tau }} - \frac{{\mathit{\boldsymbol{SK}}}}{c}} \right]} \right\} $

(7)

根据最小二乘准则, 得到如下估计结果:

$ \mathit{\boldsymbol{\hat K}} = c{\left( {{\mathit{\boldsymbol{S}}^{\rm{T}}}\mathit{\boldsymbol{S}}} \right)^{ - 1}}{\mathit{\boldsymbol{S}}^{\rm{T}}}\mathit{\boldsymbol{\hat \tau }} $

(8)

1.2 波达方向估计

在测向问题中, 经典LS法只考虑观测向量的误差, 假设系数矩阵没有误差或者干脆不考虑误差。然而, 在数据的获取过程中, 由于物理、机械、技术、仪器或作业人员的观测条件的限制, 测量误差不可避免, 观测向量、系数矩阵都有可能存在误差。因此, 在这种情况下, LS法的测向精度会大幅降低。为了弥补LS法的理论缺陷, 本文提出TLS乘波达方向估计方法。

令传感器坐标误差为e, 则有:

$ \mathit{\boldsymbol{S}} = \mathit{\boldsymbol{\hat S}} + \mathit{\boldsymbol{e}} $

(9)

其中, e=[e₁  e₂  …  e_N－1]∈$\mathbb{R}$^(N-1)×1。

TLS的基本思想是用校正向量n去干扰观测向量τ, 同时用校正矩阵e去干扰系数矩阵S, 以便对S和τ中的误差进行联合补偿, 从而抑制观测误差和系统误差对测向精度的影响。具体表达式如下:

$ c\left( {\mathit{\boldsymbol{\hat \tau }} + \mathit{\boldsymbol{n}}} \right) = \left( {\mathit{\boldsymbol{\hat S}} + \mathit{\boldsymbol{e}}} \right)\mathit{\boldsymbol{K}} $

(10)

显然, 校正向量和校正矩阵都要尽可能小。因此, TLS问题可用下列约束优化问题来描述:

$ \begin{array}{l} \mathop {\min }\limits_{\mathit{\boldsymbol{e}},\mathit{\boldsymbol{n}},K} \left\| {\mathit{\boldsymbol{e}},\mathit{\boldsymbol{n}}} \right\|_2^2 = \left\| \mathit{\boldsymbol{e}} \right\|_2^2 + \left\| \mathit{\boldsymbol{n}} \right\|_2^2\\ {\rm{subject}}\;{\rm{to}}\;c\left( {\mathit{\boldsymbol{\hat \tau }} + \mathit{\boldsymbol{n}}} \right) = \left( {\mathit{\boldsymbol{\hat S}} + \mathit{\boldsymbol{e}}} \right)\mathit{\boldsymbol{K}} \end{array} $

(11)

由式(10)可知, 原矩阵方程cτ=SK可改写为:

$ \left( {\left[ {\mathit{\boldsymbol{\hat S}},c\mathit{\boldsymbol{\hat \tau }}} \right] + \left[ {\mathit{\boldsymbol{e}},c\mathit{\boldsymbol{n}}} \right]} \right)\left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{K}}\\ { - 1} \end{array}} \right] = {\bf{0}} $

(12)

其等价为式(13)。

$ \left( {\mathit{\boldsymbol{D}} + \Delta \mathit{\boldsymbol{D}}} \right)\left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{K}}\\ { - 1} \end{array}} \right] = {\bf{0}} $

(13)

在式(13)中, 增广矩阵D=[$\mathit{\boldsymbol{\hat S}}$, c $\mathit{\boldsymbol{\hat \tau }}$]和增广校正矩阵ΔD=[e, cn]均为(N-1)×4维矩阵, 而$\mathit{\boldsymbol{p}}{\rm{ = }}\left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{K}}\\ { - 1} \end{array}} \right]$为4×1维向量。

由式(12)可知, 解向量p为增广矩阵D=[$\mathit{\boldsymbol{\hat S}}$, c$\mathit{\boldsymbol{\hat \tau }}$]的最小奇异值σ_min所对应的右奇异向量, 也是[$\mathit{\boldsymbol{\hat S}}$, c$\mathit{\boldsymbol{\hat \tau }}$]^T[$\mathit{\boldsymbol{\hat S}}$, c$\mathit{\boldsymbol{\hat \tau }}$]的最小特征值λ_min所对应的特征向量。换言之, 解向量p是下列Rayleigh商最小化的无约束优化问题的解:

$ \begin{array}{*{20}{c}} {\mathop {\min }\limits_\mathit{\boldsymbol{K}} J\left( \mathit{\boldsymbol{K}} \right) = \frac{{{{\left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{K}}\\ { - 1} \end{array}} \right]}^{\rm{T}}}{{\left[ {\mathit{\boldsymbol{\hat S}},c\mathit{\boldsymbol{\hat \tau }}} \right]}^{\rm{T}}}\left[ {\mathit{\boldsymbol{\hat S}},c\mathit{\boldsymbol{\hat \tau }}} \right]\left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{K}}\\ { - 1} \end{array}} \right]}}{{{{\left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{K}}\\ { - 1} \end{array}} \right]}^{\rm{T}}}\left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{K}}\\ { - 1} \end{array}} \right]}} = }\\ {\frac{{\left\| {\mathit{\boldsymbol{\hat SK}} - c\mathit{\boldsymbol{\hat \tau }}} \right\|_2^2}}{{\left\| \mathit{\boldsymbol{K}} \right\|_2^2 + 1}}} \end{array} $

(14)

由于K^TK=k_x²+k_y²+k_z²=1, 因此式(11)可等价为一个含约束的标准LS问题, 解向量p可简化为下列Rayleigh商最小化的无约束优化问题的解:

$ \begin{array}{*{20}{c}} {\mathop {\min }\limits_\mathit{\boldsymbol{K}} J\left( \mathit{\boldsymbol{K}} \right) = \frac{{{{\left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{K}}\\ { - 1} \end{array}} \right]}^{\rm{T}}}{{\left[ {\mathit{\boldsymbol{\hat S}},c\mathit{\boldsymbol{\hat \tau }}} \right]}^{\rm{T}}}\left[ {\mathit{\boldsymbol{\hat S}},c\mathit{\boldsymbol{\hat \tau }}} \right]\left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{K}}\\ { - 1} \end{array}} \right]}}{{{{\left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{K}}\\ { - 1} \end{array}} \right]}^{\rm{T}}}\left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{K}}\\ { - 1} \end{array}} \right]}} = }\\ {\left\| {\mathit{\boldsymbol{\hat SK}} - c\mathit{\boldsymbol{\hat \tau }}} \right\|_2^2} \end{array} $

(15)

subject to K^TK=1

用Lagrange乘数法求解上式, 定义如下目标函数:

$ \begin{array}{*{20}{c}} {J\left( \mathit{\boldsymbol{K}} \right) = \left\| {\mathit{\boldsymbol{\hat{S} K}} - c\mathit{\boldsymbol{\hat \tau }}} \right\|_2^2 + \lambda \left( {{\mathit{\boldsymbol{K}}^{\rm{T}}}\mathit{\boldsymbol{K}} - 1} \right) = }\\ {\left\| {\mathit{\boldsymbol{Dp}}} \right\|_2^2 + \lambda \left( {{\mathit{\boldsymbol{K}}^{\rm{T}}}\mathit{\boldsymbol{K}} - 1} \right)} \end{array} $

(16)

其中, λ为Lagrange乘数。由于‖Dp‖₂²=p^TD^TDp, 因此由$\frac{{\partial \mathit{\boldsymbol{J}}\left( \mathit{\boldsymbol{K}} \right)}}{{\partial \mathit{\boldsymbol{K}}}}=0$, 得到如下等式:

$ {\mathit{\boldsymbol{D}}^{\rm{T}}}\mathit{\boldsymbol{Dp}} = \mathit{\boldsymbol{\lambda p}} $

(17)

式(17)表明, 应该选择矩阵D^TD的最小特征值, 即增广矩阵D的最小奇异值的平方, 作为Lagrange乘数。TLS的解向量p是最小奇异值${\sigma _{\min }}{\rm{ = }}\sqrt {{\lambda _{\min }}} $所对应的右奇异向量。

令(N-1)×4维增广矩阵D的奇异值分解为:

$ \mathit{\boldsymbol{D}} = \mathit{\boldsymbol{U \boldsymbol{\varSigma} }}{\mathit{\boldsymbol{H}}^{\rm{T}}} $

(18)

其奇异值按照顺序σ₁≥σ₂≥σ₃≥σ₄排列, 与这些奇异值对应的右奇异向量分别为h₁、h₂、h₃、h₄。综合以上分析可知, 总体最小二乘解为p=h₄, 原矩阵方程cτ=SK的最小二乘解由式(19)给出。

$ \mathit{\boldsymbol{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over K} }} = \left[ {\begin{array}{*{20}{c}} {\frac{{{\mathit{\boldsymbol{h}}_4}\left( 1 \right)}}{{{\mathit{\boldsymbol{h}}_4}\left( 4 \right)}}}\\ {\frac{{{\mathit{\boldsymbol{h}}_4}\left( 2 \right)}}{{{\mathit{\boldsymbol{h}}_4}\left( 4 \right)}}}\\ {\frac{{{\mathit{\boldsymbol{h}}_4}\left( 3 \right)}}{{{\mathit{\boldsymbol{h}}_4}\left( 4 \right)}}} \end{array}} \right] $

(19)

根据波达方向定位在三维空间中的三角函数关系, 得到方位角和俯仰角的估计值, 其计算过程如下:

$ \begin{array}{l} \varphi = \arctan \frac{{{k_y}}}{{{k_x}}}\\ \theta = \arctan \frac{{{k_z}}}{{\sqrt {k_x^2 + k_y^2} }} \end{array} $

(20)

2 TLS克拉美罗下界

CRLB是无源定位系统中一个重要的误差评估方法, 为待估计参数的无偏估计提供了误差方差的下界, 可用来分析各种定位方法所能达到的最优理论性能。

从组合矢量V=[τ^T, S^T]^T的概率密度函数开始推导CRLB, 由于n和e均服从零均值高斯分布且相互独立, 则矢量$\Delta \mathit{\boldsymbol{V}} = \mathit{\boldsymbol{V}} - \mathit{\boldsymbol{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over V} }}$是协方差矩阵为块对角矩阵的高斯随机矢量, 分块对角矩阵$\mathit{\boldsymbol{Q}} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{\tau }}}}&\bf{0}\\ \bf{0}&{{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{S}}}} \end{array}} \right]$, 块对角协方差矩阵的对角块如下:

$ {\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{S}}} = \left[ {\begin{array}{*{20}{c}} {\sigma _{{e_1}}^2}&{}&{}&{}\\ {}&{\sigma _{{e_2}}^2}&{}&{}\\ {}&{}& \ddots &{}\\ {}&{}&{}&{\sigma _{{e_{N - 1}}}^2} \end{array}} \right] $

$ {\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{\tau }}} = \left[ {\begin{array}{*{20}{c}} {\sigma _{{n_1}}^2}&{}&{}&{}\\ {}&{\sigma _{{n_2}}^2}&{}&{}\\ {}&{}& \ddots &{}\\ {}&{}&{}&{\sigma _{{n_{N - 1}}}^2} \end{array}} \right] $

令Ψ=[γ, S^T]^T, 将概率密度函数相对于Ψ进行参数化处理可以得到CRLB(γ), 数据矢量V=[τ^T, S^T]^T的条件概率密度函数的对数表达式如下:

$ \begin{array}{*{20}{c}} {\ln f\left( {\mathit{\boldsymbol{V}}\left| \mathit{\boldsymbol{ \boldsymbol{\varPsi} }} \right.} \right) = {k_1} - \frac{1}{2}{{\left( {\mathit{\boldsymbol{\hat \tau }} - \mathit{\boldsymbol{\tau }}} \right)}^{\rm{T}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{\tau }}^{ - 1}\left( {\mathit{\boldsymbol{\hat \tau }} - \mathit{\boldsymbol{\tau }}} \right) + {k_2} - }\\ {\frac{1}{2}{{\left( {\mathit{\boldsymbol{\hat S}} - \mathit{\boldsymbol{S}}} \right)}^{\rm{T}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{S}}^{ - 1}\left( {\mathit{\boldsymbol{\hat S}} - \mathit{\boldsymbol{S}}} \right)} \end{array} $

(21)

其中${k_1} = - \frac{1}{2}\ln ({(2{\rm{ {\rm{ \mathsf{ π} }} }} )^{N - 1}}\left| {{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{\tau }}}} \right|), {k_2} = - \frac{1}{2}\ln ({(2{\rm{{\rm{ \mathsf{ π} }} }} )^{N - 1}}\left| {{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{S}}}} \right|)$。Fisher信息矩阵如下:

$ \mathit{\boldsymbol{F}}\left( \mathit{\boldsymbol{ \boldsymbol{\varPsi} }} \right) = - E\left[ {\frac{{{\partial ^2}\ln f\left( {\mathit{\boldsymbol{V}}\left| \mathit{\boldsymbol{ \boldsymbol{\varPsi} }} \right.} \right)}}{{\partial \mathit{\boldsymbol{ \boldsymbol{\varPsi} }}\partial {\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}^{\rm{T}}}}}} \right] $

(22)

因此, 有如下表达式:

$ {\mathit{\boldsymbol{F}}^{ - 1}}\left( \mathit{\boldsymbol{ \boldsymbol{\varPsi} }} \right) = - E{\left[ {\frac{{{\partial ^2}\ln f\left( {\mathit{\boldsymbol{V}}\left| \mathit{\boldsymbol{ \boldsymbol{\varPsi} }} \right.} \right)}}{{\partial \mathit{\boldsymbol{ \boldsymbol{\varPsi} }}\partial {\mathit{\boldsymbol{ \boldsymbol{\varPsi} }}^{\rm{T}}}}}} \right]^{ - 1}} = {\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{F}}_1}}&{{\mathit{\boldsymbol{F}}_2}}\\ {{\mathit{\boldsymbol{F}}_3}}&{{\mathit{\boldsymbol{F}}_4}} \end{array}} \right]^{ - 1}} $

(23)

其中:

$ {\mathit{\boldsymbol{F}}_1} = {\left( {\frac{{\partial \mathit{\boldsymbol{\tau }}}}{{\partial \mathit{\boldsymbol{\gamma }}}}} \right)^{\rm{T}}}\mathit{\boldsymbol{Q}}_\tau ^{ - 1}\frac{{\partial \mathit{\boldsymbol{\tau }}}}{{\partial \mathit{\boldsymbol{\gamma }}}} + {\left( {\frac{{\partial \mathit{\boldsymbol{S}}}}{{\partial \mathit{\boldsymbol{\gamma }}}}} \right)^{\rm{T}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{S}}^{ - 1}\frac{{\partial \mathit{\boldsymbol{S}}}}{{\partial \mathit{\boldsymbol{\gamma }}}} $

(24)

$ {\mathit{\boldsymbol{F}}_2} = \mathit{\boldsymbol{F}}_3^{\rm{T}} = {\left( {\frac{{\partial \mathit{\boldsymbol{\tau }}}}{{\partial \mathit{\boldsymbol{\gamma }}}}} \right)^{\rm{T}}}\mathit{\boldsymbol{Q}}_\tau ^{ - 1}\frac{{\partial \mathit{\boldsymbol{\tau }}}}{{\partial {\mathit{\boldsymbol{S}}^{\rm{T}}}}} + {\left( {\frac{{\partial \mathit{\boldsymbol{S}}}}{{\partial \mathit{\boldsymbol{\gamma }}}}} \right)^{\rm{T}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{S}}^{ - 1}\frac{{\partial \mathit{\boldsymbol{S}}}}{{\partial \mathit{\boldsymbol{S}}}} $

(25)

$ {\mathit{\boldsymbol{F}}_4} = {\left( {\frac{{\partial \mathit{\boldsymbol{\tau }}}}{{\partial \mathit{\boldsymbol{S}}}}} \right)^{\rm{T}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{\tau }}^{ - 1}\frac{{\partial \mathit{\boldsymbol{\tau }}}}{{\partial {\mathit{\boldsymbol{S}}^{\rm{T}}}}} + {\left( {\frac{{\partial \mathit{\boldsymbol{S}}}}{{\partial \mathit{\boldsymbol{S}}}}} \right)^{\rm{T}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{S}}^{ - 1}\frac{{\partial \mathit{\boldsymbol{S}}}}{{\partial \mathit{\boldsymbol{S}}}} $

(26)

$\frac{{\partial \mathit{\boldsymbol{\tau }}}}{{\partial \mathit{\boldsymbol{\gamma }}}}$为(N－1)×2维矩阵, 如下:

$ \frac{{\partial \mathit{\boldsymbol{\tau }}}}{{\partial \mathit{\boldsymbol{\gamma }}}} = \left[ {\frac{{\partial {\mathit{\boldsymbol{\tau }}_1}}}{{\partial \mathit{\boldsymbol{\gamma }}}},\frac{{\partial {\mathit{\boldsymbol{\tau }}_2}}}{{\partial \mathit{\boldsymbol{\gamma }}}}, \cdots ,\frac{{\partial {\mathit{\boldsymbol{\tau }}_{N - 1}}}}{{\partial \mathit{\boldsymbol{\gamma }}}}} \right] $

(27)

其中, $\frac{{\partial {\tau _i}}}{{\partial \mathit{\boldsymbol{\gamma }}}} = \left[ {\frac{{\partial {\tau _i}}}{{\partial \theta }}, \frac{{\partial {\tau _i}}}{{\partial \varphi }}} \right]$。

$\frac{{\partial \mathit{\boldsymbol{\tau }}}}{{\partial \mathit{\boldsymbol{S}}}}$为(N－1)×3(N－1)维矩阵, 如下:

$ \frac{{\partial \mathit{\boldsymbol{\tau }}}}{{\partial \mathit{\boldsymbol{S}}}} = \left[ {\frac{{\partial {\tau _1}}}{{\partial {S_1}}},\frac{{\partial {\tau _2}}}{{\partial {S_2}}}, \cdots ,\frac{{\partial {\mathit{\boldsymbol{\tau }}_{N - 1}}}}{{\partial {\mathit{\boldsymbol{S}}_{N - 1}}}}} \right] $

(28)

其中, $\frac{{\partial {\tau _i}}}{{\partial {S_i}}} = \left[ {\frac{{\partial {\tau _i}}}{{\partial {x_i}}}, \frac{{\partial {\tau _i}}}{{\partial {y_i}}}, \cdots , \frac{{\partial {\tau _i}}}{{\partial {z_i}}}} \right]$。

$\frac{{\partial \mathit{\boldsymbol{S}}}}{{\partial \mathit{\boldsymbol{\gamma }}}}$为3(N－1)×2维矩阵, 如下:

$ \frac{{\partial \mathit{\boldsymbol{S}}}}{{\partial \mathit{\boldsymbol{\gamma }}}} = {\mathit{\boldsymbol{O}}_{3\left( {N - 1} \right) \times 2}} $

(29)

$\frac{{\partial \mathit{\boldsymbol{S}}}}{{\partial \mathit{\boldsymbol{S}}}}$为3(N－1)×3(N－1)维矩阵, 如下:

$ \frac{{\partial \mathit{\boldsymbol{S}}}}{{\partial \mathit{\boldsymbol{S}}}} = {\mathit{\boldsymbol{I}}_{3\left( {N - 1} \right) \times 3\left( {N - 1} \right)}} $

(30)

令${A_i} = \frac{{\partial {\tau _i}}}{{\partial \theta }}、{B_i} = \frac{{\partial {\tau _i}}}{{\partial \varphi }}$, 则有如下公式:

$ {\mathit{\boldsymbol{F}}_1} = {\left( {\frac{{\partial \mathit{\boldsymbol{\tau }}}}{{\partial \mathit{\boldsymbol{\gamma }}}}} \right)^{\rm{T}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{\tau }}^{ - 1}\frac{{\partial \mathit{\boldsymbol{\tau }}}}{{\partial \mathit{\boldsymbol{\gamma }}}} = \left[ {\begin{array}{*{20}{c}} {\sum\limits_{i = 1}^{N - 1} {\frac{1}{{\sigma _{{n_i}}^2}}A_i^2} }&{\sum\limits_{i = 1}^{N - 1} {\frac{1}{{\sigma _{{n_i}}^2}}{A_i}{B_i}} }\\ {\sum\limits_{i = 1}^{N - 1} {\frac{1}{{\sigma _{{n_i}}^2}}{B_i}{A_i}} }&{\sum\limits_{i = 1}^{N - 1} {\frac{1}{{\sigma _{{n_i}}^2}}B_i^2} } \end{array}} \right] $

(31)

因此, F^－1(Ψ)是[3(N－1)+2]×[3(N－1)+2]维矩阵, 其左上角的2×2维块矩阵为未知声源方向矢量F^－1(γ), 通过分块矩阵求逆公式, 得到如下公式:

$ \begin{array}{l} CRLB\left( \mathit{\boldsymbol{\gamma }} \right) = {\left( {{\mathit{\boldsymbol{F}}_1} - {\mathit{\boldsymbol{F}}_2}\mathit{\boldsymbol{F}}_4^{ - 1}\mathit{\boldsymbol{F}}_2^{\rm{T}}} \right)^{ - 1}} = \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\mathit{\boldsymbol{F}}_1^{ - 1} + \mathit{\boldsymbol{F}}_1^{ - 1}{\mathit{\boldsymbol{F}}_2}{\left( {{\mathit{\boldsymbol{F}}_4} - {\mathit{\boldsymbol{F}}_2}} \right)^{ - 1}}\mathit{\boldsymbol{F}}_2^{\rm{T}}\mathit{\boldsymbol{F}}_1^{ - 1} \end{array} $

(32)

分析式(32)可知, F₁^－1中只包含带估计目标Ψ=[γ^T, S^T]^T中γ的信息, 各传感器误差引入的信息则包含在F₂与F₄中, 因此F₁^－1为传感器位置无误差时的克拉美罗下界。式(32)中的第2项为传感器位置误差所引起的损失差ΔCRLB(γ), CRLB(γ)的迹为到达时间差测向任何无偏估计MSE的下界。

3 仿真结果与分析

在确定了总体最小二乘算法后, 对该测向方法进行仿真。整个仿真过程分为2个步骤。首先, 模拟产生8个传感器的坐标, 生成一个随机的信号源。然后, 结合以上信息对该方法进行计算机仿真, 检验本文方法的数学模型是否正确、各个环节是否能够满足设计要求, 并将本文算法与其他算法进行比较, 找出最佳测向算法。

模拟实验以声源测向为背景, 具体参数设置如下:

1) 采用Matlab工具进行仿真, 每轮蒙特卡洛模拟次数均为10 000次。

2) 设置8个传感器, 假设传感器位置的X轴和Y轴坐标在半径30 m的圆形内随机选取, Z轴坐标在[－5, 5]的范围内随机确定, 单位为m。

3) 随机产生声源方向, 方位角设置在[0°, 90°]之间, 俯仰角设置在[－30°, 30°]之间。

4) TDOA测量误差服从高斯分布, 标准偏差在[0, 20]范围内随机选取, 单位为ms。

5) 系统误差服从高斯分布, 标准偏差在[0, 10]范围内随机选取, 单位为m。

3.1 系统误差对测向精度的影响

累积分布函数(Cumulative Distribution Function, CDF)定义为误差小于给定的定位误差容限的概率, 图 2、图 3给出10 000次随机测试下的测向结果CDF图。图 2是传感器在有误差和无误差2种情况下, 通过LS算法得到的俯仰角的CDF图, 从图 2可以看出, 相对于无系统误差的情况, 含系统误差时所得的测向精度较低。图 3是传感器在有误差和无误差2种情况下, 通过LS算法所得到的方位角的CDF图像, 由图 3可知, 传感器位置误差对测向精度影响较大。综合以上分析, 系统误差将会降低LS法的测向精度。

3.2 TLS方法测向精度

假设已知传感器坐标的方差σ_e1²=5 m, 分别改变观测信号的噪声误差值, 调整噪声方差σ_ni², 使其在0.005 6到0.017 8之间变化, 即调整10 lg σ_ni²从－45 dB到－35 dB之间变化, 分别通过LS方法、TLS方法和LS-Taylor方法^[12]进行仿真, 得到俯仰角和方位角的均方根误差(Root Mean Square Error, RMSE), 如图 4、图 5所示。从图 4、图 5可以看出, 随着信号噪声的增加, RMSE值也随之增大, 相比于LS闭式解和LS-Taylor迭代方法, TLS方法的定位误差更小, 即本文方法抗传感器位置误差噪声的效果更好, 且无需迭代计算。

4 结束语

本文针对经典最小二乘法在传感器位置存在误差时测向精度低的问题, 提出总体最小二乘波达方向估计方法TLS, 同时考虑系统误差和测量误差, 弥补LS法不考虑系统误差的理论缺陷。仿真结果表明, 与经典LS方法和改进的LS-Taylor算法相比, 该方法的定位性能有较大提升。另外, 本文TLS方法不需要进行迭代计算或者获得声源位置的粗略估计, 保证了定位结果的全局收敛性。下一步将研究基于TLS方法与Taylor级数展开的混合定位方法, 以提高定位结果的抗误差能力。