假定某时刻有M个传感器对一移动辐射源开展联合协同无源定位，定位场景如图1所示。

图  1  多平台联合定位示意图

记移动辐射源在该时刻的三维位置向量为${{ {u}}^o} = {[{x^o},{y^o},{z^o}]^{\rm T} }$，三维速度向量为${{\dot{ u}}^o} = {[{\dot x^o},{\dot y^o},{\dot z^o}]^{\rm T} }$。记该时刻第$i{\rm{(}}i = 1,2, \cdots ,M{\rm{)}}$个传感器的所在真实位置为${ {s}}_i^o = {[x_i^o,y_i^o,z_i^o]^{\rm T} }$，真实速度为${\dot{ s}}_i^o = {[\dot x_i^o,\dot y_i^o,\dot z_i^o]^{\rm T} }$。一般地，定位算法并不能获取上述传感器的真实坐标与真实速度，而只能获取传感器所在平台的导航设备提供的位置向量与速度向量。记该时刻导航设备给出的第$i$个传感器的所在位置为${{ {s}}_i} = {[{x_i},{y_i},{z_i}]^{\rm T} }$，速度为${{\dot{ s}}_i} = {[{\dot x_i},{\dot y_i},{\dot z_i}]^{\rm T} }$。

为简便行文，将各传感器的位置向量与速度向量收集在一起，并记为${ {\beta }} = {[{{ {s}}^{\rm T} },{{\dot{ s}}^{\rm T} }]^{\rm T} }$，其中${ {s}} = {[{ {s}}_1^{\rm T} , { {s}}_2^{\rm T} , \cdots ,}$${{ {s}}_M^{\rm T} ]^{\rm T} }$；${\dot{ s}} = {[{\dot{ s}}_1^{\rm T} , {\dot{ s}}_2^{\rm T} , \cdots ,{\dot{ s}}_M^{\rm T} ]^{\rm T} }$。若记传感器位置误差向量为$\Delta { {\beta }}$，那么$\Delta { {\beta }}\! =\! { {\beta }} - {{ {\beta }}^o} \!=\! {{\rm{[}}\Delta {{ {s}}^{\rm{T}}},\Delta {{\dot{ s}}^{\rm{{\rm T}}}}{\rm{]}}^{\rm{{\rm T}}}}$，其中，$\Delta { {s}}= $${[\Delta { {s}}_1^{\rm T} ,\Delta { {s}}_2^{\rm T} , \cdots ,\Delta { {s}}_M^{\rm T} ]^{\rm T} }$；$\Delta {\dot{ s}} = {[\Delta {\dot{ s}}_1^{\rm T} ,\Delta {\dot{ s}}_2^{\rm T} , \cdots ,\Delta {\dot{ s}}_M^{\rm T} ]^{\rm T} }; \Delta {{ {s}}_i} = {{ {s}}_i} - $${ {s}}_i^o$；$\Delta {{\dot{ s}}_i} = {{\dot{ s}}_i} - {\dot{ s}}_i^o$。一般地，$\Delta { {\beta }}$通常认为是零均值的高斯向量，其协方差矩阵为${{ {Q}}_{ {\beta }}}$。

记$\theta _i^o$与$\phi _i^o$分别表示该时刻辐射源相对于第$i{\rm{(}}i = 1,2, \cdots ,M{\rm{)}}$个传感器的方位角与高低角的真实值，表示为：

式中，$\theta _i^o \in \left( { - {\text π} ,{\text π} } \right)$；$\phi _i^o \in \left( {{{ - {\text π} } / {2,}}{{\text π} / 2}} \right)$。在实际中，所观测到的方位角与高低角均带有测量噪声，即有${ {\theta }} = {{ {\theta }}^o} + {{ {n}}_\theta }$与${ {\phi }} = {{ {\phi }}^{\rm{o}}} + {{ {n}}_\phi }$。记上式中的${ {\theta }} = {[{\theta _1}, {\theta _2}, }$$ \cdots ,{\theta _M}]^{\rm T}$与${ {\phi}} = {{\rm{[}}{\phi _1}, {\phi _2}, \cdots ,{\phi _M}{\rm{]}}^{\rm{T}}}$分别表示该时刻方位角与高低角的测量向量，而${{ {\theta }}^o} = {[\theta _1^o, \theta _2^o,\cdots ,\theta _M^o]^{\rm T} }$与${{ {\phi}} ^{\rm{o}}} = {{\rm{[}}\phi _1^o, \phi _2^o, \cdots ,\phi _M^o{\rm{]}}^{\rm{T}}}$分别表示该时刻方位角与高低角的真实值向量，以及${{ {n}}_\theta } = {[{n_{\theta 1}}, {n_{\theta 2}},\cdots ,{n_{\theta M}}]^{\rm T} }$与${{ {n}}_\phi } = {[{n_{\phi 1}}, {n_{\phi 2}}, \cdots ,{n_{\phi M}}]^{\rm T} }$为对应的零均值高斯噪声向量，其协方差阵分为记为${{ {Q}}_\theta }$、${{ {Q}}_\phi }$。并且假设测量噪声${{ {n}}_\theta }$与${{ {n}}_\phi }$跟传感器位置误差向量$\Delta { {\beta }}$统计独立。

设1号传感器为参考传感器，那么该时刻则有M−1个TDOA观测量作为无源定位算法的输入。${r_{i1}}$为第$i{\rm{(}}i = 2,3, \cdots ,M{\rm{)}}$个传感器与第1个传感器之间的TDOA观测量。一般地，${r_{i1}} = r_{i1}^o + {n_{i1}}$，其中${n_{i1}}$为测量噪声，$r_{i1}^o$为第i个传感器与第1个传感器之间的TDOA真值，并且有$r_{i1}^o = r_i^o - r_1^o$，其中$r_k^o = $$ \left\| {{{ {u}}^o} - { {s}}_k^o} \right\|$代表该时刻移动辐射源与第$k{\rm{(}}k = 1,2, \cdots ,M{\rm{)}}$个传感器之间的欧氏距离。记在该时刻M−1个TDOA观测向量为${ {r}} = {{\rm{[}}{r_{21}}, {r_{31}},\cdots ,{r_{M1}}{\rm{]}}^{\rm{T}}}$，则${ {r}} = {{ {r}}^o} + {{ {n}}_t}$，其中${{ {r}}^o} = {\rm{[}}r_{21}^o, r_{31}^o,\cdots ,$ $r_{M1}^o{{\rm{]}}^{\rm{T}}}$为TDOA真实值向量，${{ {n}}_t} = {\rm{[}}{n_{21}}, {n_{31}},\cdots , $${n_{M1}}{\rm{]}}^{\rm{T}}$为对应的观测噪声向量。假设噪声向量${{ {n}}_t}$为零均值、协方差矩阵为${{ {Q}}_t}$的高斯分布向量，并且与传感器位置误差向量$\Delta { {\beta }}$、方位角噪声向量${{ {n}}_\theta }$以及高低角噪声向量${{ {n}}_\phi }$统计独立。

同样地，设1号传感器为参考传感器，那么该时刻则有M−1个FDOA观测量作为无源定位算法的输入。${{\dot r}_{i1}}$为第$i{\rm{(}}i = 2,3, \cdots ,M{\rm{)}}$个传感器与第1个传感器之间的FDOA观测量。一般地，${{\dot r}_{i1}} = {\dot r}_{i1}^o + {\dot n_{i1}} = {\dot r}_i^o - {\dot r}_1^o + {\dot n_{i1}}$，其中${\dot n_{i1}}$为测量噪声，${\dot r}_{i1}^o$为第i个传感器与第1个传感器之间的FDOA真值，并且有${\dot r}_k^o = {{{{{\rm{(}}{{ {u}}^o} - { {s}}_k^o{\rm{)}}}^{\rm{T}}}{\rm{(}}{{{\dot{ u}}}^o} - {\dot{ s}}_k^o{\rm{)}}} / {r_k^o}}$，$k = 1,2, \cdots ,M $。记在该时刻M−1个FDOA观测向量为${\dot{ r}} = {{\rm{[}}{{\dot r}_{21}}, {{\dot r}_{31}},\cdots ,{{\dot r}_{M1}}{\rm{]}}^{\rm{T}}}$，则${\dot{ r}} = {{\dot{ r}}^o} + {{ {n}}_\nu }$，${{\dot{ r}}^o} = {{\rm{[}}{\dot r}_{21}^o, {\dot r}_{31}^o,\cdots ,{\dot r}_{M1}^o{\rm{]}}^{\rm{T}}}$为FDOA真实值向量，${{ {n}}_\nu } = {{\rm{[}}{\dot n_{21}}, {\dot n_{31}},\cdots ,{\dot n_{M1}}{\rm{]}}^{\rm T}}$为对应的观测噪声向量。假设噪声向量${{ {n}}_\nu }$为零均值、协方差矩阵为${{ {Q}}_\nu }$的高斯分布向量，并且与传感器位置误差向量$\Delta { {\beta }}$、方位角噪声向量${{ {n}}_\theta }$、高低角噪声向量${{ {n}}_\phi }$以及TDOA噪声向量${{ {n}}_t}$统计独立。

需要说明的是，许多文献(如文献[12])在定义TDOA表达式时将光速c也包含在内，在定义FDOA表达式时同时将光速c和辐射源信号载频${f_o}$也包含在内。而这两类定义在已知辐射源信号载频${f_o}$时没有本质区别。因此，本文不特意区分上述两类定义，其含义可由上下文得到。

综上，为简化算法推导，将上述3类观测量汇集起来。记${ {\alpha }} = {[{{ {\theta }}^{\rm{T}}},{{ {\phi}} ^{\rm{T}}},{{ {r}}^{\rm{T}}},{{\dot{ r}}^{\rm{T}}}]^{\rm{T}}}$为该时刻的AOA、TDOA及FDOA测量值向量，$\Delta { {\alpha }} = {{\rm{[}}{ {n}}_\theta ^{\rm{T}},n_\phi ^{\rm{T}},{ {n}}_t^{\rm{T}},{ {n}}_\nu ^{\rm{T}}{\rm{]}}^{\rm{T}}}$为对应的观测噪声向量。一般地，观测噪声向量$\Delta { {\alpha }}$被建模为零均值、协方差矩阵${{ {Q}}_\alpha }$为高斯向量，其中${{ {Q}}_\alpha }$的定义为：

那么，移动辐射源AOA-TDOA-FDOA联合无源定位算法所涉问题为：已知某一时刻的AOA、TDOA及FDOA测量向量${ {\alpha }}$与传感器位置向量${ {\beta }}$，以及二者对应的噪声/误差向量的协方差矩阵${{ {Q}}_\alpha }$与${{ {Q}}_{ {\beta }}}$，尽可能准确地估计该时刻的移动辐射源位置${{ {u}}^o}$。

移动辐射源AOA-TDOA-FDOA联合定位算法采用两阶段估计流程。第一阶段，首先建立起移动辐射源位置参数${{ {u}}^o}$、${{\dot{ u}}^o}$与观测向量${ {\alpha }}$及传感器位置向量${ {\beta }}$的伪线性关系，进而得到${{ {u}}^o}$与${{\dot{ u}}^o}$粗略估计值；第二阶段，修正第一阶段的粗估计值，从而精确估计出某时刻移动辐射源的位置。

方位角$\theta _i^o$、高低角$\phi _i^o$与移动辐射源位置参数${{ {u}}^o}$的伪线性关系可以通过以下代数变换推导得出。由式(1)可知，对观测真值$\theta _i^o$取正切运算，有：

式中，$\theta _i^o$、$x_i^o$以$y_i^o$及均为真实值，实际中并不可得，用带噪值或带误差值代替上述各量，则式(3)变为：

式中，${\eta _{\theta i}}$为测量噪声${n_{\theta i}}$以及传感器位置误差$\Delta {{ {s}}_i}$带来的近似误差。记${{ {g}}_{\theta i}} = {{\rm{[}} - \sin {\theta _i},\cos {\theta _i},0{\rm{]}}^{\rm{T}}}$，那么式(4)可以等价表达为：

当测量噪声${n_{\theta i}}$很弱时，$\cos \left( {{n_{\theta i}}} \right) \approx 1$及$\sin \left( {{n_{\theta i}}} \right) \approx $$ {n_{\theta i}}$，$\sin {\theta _i} \approx \;\sin \theta _i^o + {n_{\theta i}}\cos \theta _i^o$及$\cos {\theta _i} \approx \cos \theta _i^o -$$ {n_{\theta i}}\sin \theta _i^o$。将上述两近似等式以及${x_i} = x_i^o + \Delta {x_i},{y_i} = y_i^o + $$ \Delta {y_i}$带入式(4)，有${\eta _{\theta i}} = \cos \theta _i^o\Delta {y_i}- \sin \theta _i^o\Delta {x_i} + {\rm{\{ (}}{x^o} - x_i^o{\rm{)}}\cos \theta _i^o{\rm{ + }}$$ ({y^o} - y_i^o{\rm{)}}\sin \theta _i^o{\rm{\} }}{n_{\theta i}}$。再利用图1中的一个几何关系式${\rm{(}}{x^o} \!-\! x_i^o{\rm{)}}\cos \theta _i^o \!+\! {\rm{(}}{y^o} \!-\! y_i^o{\rm{)}}\sin \theta _i^o = l_i^o{\cos ^2}\theta _i^o + l_i^o{\sin ^2}\theta _i^o =$$ l_i^o = r_i^o\cos \phi _i^o$，并定义${ {g}}_{\theta i}^o = {{\rm{[}} - \sin \theta _i^o,\cos \theta _i^o,0{\rm{]}}^{\rm{T}}}$，那么有：

式(5)与式(6)分别建立了向量${{ {u}}^o}$与方位角${\theta _i}$之间的近似等式关系，以及近似误差${\eta _{\theta i}}$与测量噪声${n_{\theta i}}$及传感器位置误差$\Delta {{ {s}}_i}$之间的函数关系。

同样地，对式(1)中的观测真值$\phi _i^o$取正切运算，并用${\theta _i},{\phi _i},{x_i},{y_i},{z_i}$分别代替其真实值，有：

式中，${{ {g}}_{\phi i}} = {{\rm{[}} - \sin {\phi _i}\cos {\theta _i}, - \sin {\phi _i}\sin {\theta _i},\cos {\phi _i}{\rm{]}}^{\rm{T}}}$；${\eta _{\phi i}}$为测量噪声${n_{\phi i}}$以及传感器位置误差$\Delta {{ {s}}_i}$带来的近似误差。当测量噪声${n_{\phi i}}$很弱时，$\sin {\phi _i} = \sin \phi _i^o + {n_{\phi i}}\cos \phi _i^o$及$\cos {\phi _i} = \cos\phi _i^o - {n_{\phi i}}\sin \phi _i^o$近似成立，并将上述两近似关系及${{ {s}}_i} = { {s}}_i^o + \Delta {{ {s}}_i}$代入式(7)，忽略二阶误差项，有：

式中，${ {g}}_{\phi i}^o = {{\rm{[}} - \sin \phi _i^o\cos \theta _i^o, - \sin \phi _i^o\sin \theta _i^o,\cos \phi _i^o{\rm{]}}^{\rm{T}}}$。上式推导过程中使用了图1中两个几何关系$0 = {\rm{(}}{x^o} - x_i^o{\rm{)}}\times$$\sin \theta _i^o - {\rm{ (}}{y^o} - y_i^o{\rm{)}}\cos \theta _i^o$ 及${\rm{(}}{x^o} - \;x_i^o{\rm{)}}\cos \theta _i^o\cos \phi _i^o+ {\rm{(}}{y^o} - y_i^o{\rm{)}}\times$$ \sin \theta _i^o\cos \phi _i^o + {\rm{(}}{z^o} - z_i^o{\rm{)}}\sin \phi _i^o = r_i^o$。式(7)与式(8)分别建立了向量${{ {u}}^o}$与高低角${\phi _i}$之间的近似等式关系，以及近似误差${\eta _{\phi i}}$与测量噪声${n_{\phi i}}$及传感器位置误差$\Delta {{ {s}}_i}$之间的函数关系。

事实上，${ {s}}_i^o - {{ {u}}^o}$为三维空间中的一个向量，而向量${ {g}}_{\theta i}^o$与向量${ {g}}_{\phi i}^o$是此三维空间中既彼此正交又同时与向量${ {s}}_i^o - {{ {u}}^o}$成正交关系。当存在测量噪声和传感器位置误差时，定位算法只能得到${{ {g}}_{\theta i}}$与${{ {g}}_{\phi i}}$。向量${{ {g}}_{\theta i}}$与向量${{ {g}}_{\phi i}}$虽然彼此正交，但却不再与向量${ {s}}_i^o - {{ {u}}^o}$正交，这也就是残差${\eta _{\theta i}}$与${\eta _{\phi i}}$的产生原因。

如文献[10-12]及文献[7]所示，TDOA观测量${r_{i1}},$ $i{\rm{(}}i = 2,3, \cdots ,M{\rm{)}}$可以通过如下方式进行线性化处理。首先，对等式$r_{i1}^o + r_1^o = r_i^o$的两边同时取平方运算，有：

将上式中的$r_{i1}^o$用${r_{i1}} - {n_{i1}}$代替，${ {s}}_i^o$用${{ {s}}_i} - \Delta {{ {s}}_i}$代替，忽略二阶误差项整理可得：

式中，${\varepsilon _{t,i}}$代表近似误差，其定义为：

借助于辅助变量$r_1^o$，式(10)表征了TDOA观测量${r_{i1}}$与辐射源位置向量${{ {u}}^o}$之间的近似线性关系，同时式(11)也给出了近似误差${\varepsilon _{t,i}}$与测量噪声${n_{i1}}$及传感器位置误差$\Delta {{ {s}}_i}$之间的函数表达关系。

对于FDOA观测量${{\dot r}_{i1}},i{\rm{(}}i = 2,3, \cdots ,M{\rm{)}}$的线性化处理方式为同时对式(10)及式(11)的右端对时间求导数^[10-12]，即：

借助于辅助变量$r_1^o$与${\dot r}_1^o$，式(12)不仅给出了FDOA观测量${{\dot r}_{i1}}$与移动辐射源位置参数${{ {u}}^o}$与${{\dot{ u}}^o}$之间的近似线性关系，同时也给出了近似误差${\varepsilon _{\nu ,i}}$、测量噪声${n_{i1}}$、${\dot n_{i1}}$、传感器位置误差$\Delta {{ {s}}_i}$和速度误差$\Delta {{\dot{ s}}_i}$之间的函数表达关系。

1)第一阶段：粗估计。

从上面的分析中可以知道，借助于辅助变量$r_1^o$与${\dot r}_1^o$，移动辐射源位置参数${{ {u}}^o}$、${{\dot{ u}}^o}$与无源定位输入观测向量${ {\alpha }}$及传感器位置向量${ {\beta }}$之间能建立起近似线性关系，并且近似误差也有显示表达式。因此，在第一阶段可以选取${ {\varphi }}_1^o = {[{{{ {u}}^{o{\rm T}}}}, r_1^o,{{{\dot{ {u}}}^{o{\rm T}}} },{{\dot {r}}_1^o]^{\rm T}} }$为移动辐射源位置参数的估计参数向量。

记${ {\eta }} \;=\; {{\rm{[}}{\eta _{\theta 1}}, {\eta _{\theta 2}}, \cdots \;,{\eta _{\theta M}},\;{\eta _{\phi 1}},{\eta _{\phi 2}},\; \cdots \;,{\eta _{\phi M}}{\rm{]}}^{\rm{T}}},\;{{ {\varepsilon }}_t} = $${\rm{[}}{\varepsilon _{t,2}},{\varepsilon _{t,3}},\;\cdots ,{\varepsilon _{t,M}}{\rm{]}}^{\rm{T}}$，以及${{ {\varepsilon }}_\nu } = {[{\varepsilon _{\nu ,2}}, {\varepsilon _{\nu ,3}}, \cdots ,{\varepsilon _{\nu ,M}}]^{\rm T} }$，并定义${{ {\varepsilon }}_1} = {\rm{[}}{{ {\eta }}^{\rm{T}}},{ {\varepsilon }}_t^{\rm{T}},$${ {\varepsilon }}_\nu ^{\rm{T}}{\rm{]}}^{\rm{T}}$，那么联合式(5)~式(8)、式(10)~式(12)，有：

式中，

式(14)中的矩阵${{ O}_{i \times j}}$代表$i \times j$维的全零矩阵，矩阵${{ {h}}_{\theta \phi }},{{ {G}}_{\theta \phi }},{{ {h}}_t},{{ {h}}_\nu },{{ {G}}_t},{{ {G}}_\nu }$分别定义为：

式中，$ \odot $表示向量Schur乘法；${{{\textit{1}}}_{m \times 1}}$与${{{\textit{0}}}_{m \times 1}}$分别表征m维取值全1的列向量与全0的列向量；${{ {G}}_\theta } = {{\rm{[}}{{ {g}}_{\theta 1}}, {{ {g}}_{\theta 2}},\cdots ,{{ {g}}_{\theta M}}{\rm{]}}^{\rm{T}}},{{ {G}}_\phi } = {{\rm{[}}{{ {g}}_{\phi 1}}, {{ {g}}_{\phi 2}},\cdots ,{{ {g}}_{\phi M}}{\rm{]}}^{\rm{T}}} $；S=${\rm{[}}{{ {s}}_1}, {{ {s}}_2}, \cdots ,{{ {s}}_M}{\rm{]}}^{\rm{T}}$。

式(14)中的矩阵${ {B}}_{\theta \phi }^o,{ {D}}_{\theta \phi }^o$的定义分别为：

式中，$ { {B}}_{\theta s}^o\! = \!{\rm{diag}}{\rm{(}}{ {G}}_\theta ^o{\rm{); }}{ {B}}_{\phi s}^o \!=\! {\rm{diag}}{\rm{(}}{ {G}}_\phi ^o{\rm{); }}{ {B}}_{\phi \omega }^o \!\!=\! {\rm{diag}}{\rm{([}}r_1^o, r_2^o, \cdots $$ r_M^o{\rm{]); }}{ {B}}_{\theta n}^o = {\rm{diag}}{\rm{([}}r_1^o\cos \phi _1^o, r_2^o\cos \phi _2^o, \cdots ,r_M^o\cos \phi _M^o{\rm{])}} $。

式(14)中的矩阵${{ {B}}^o},{{\dot{ B}}^o}$分别定义为：

式(14)中的矩阵${{ {D}}^o},{{\dot{ D}}^o}$分别按如下方式定义其每一行，即矩阵${{ {D}}^o},{{\dot{ D}}^o}$的第$i{\rm{(}}i = 1,2, \cdots ,M - 1{\rm{)}}$行定义为：

式中，$ {{ {p}}_k} = ({{ {u}}^o} - { {s}}_k^o);{{\dot{ p}}_k} = ({{\dot{ u}}^o} - {\dot{ s}}_k^o),k = 1,2, \cdots ,M $。

基于式(13)，可以通过加权最小二乘来求解参数${ {\varphi }}_1^o$。权矩阵通常选为${ {W}} = {{\rm{(}}{ {B}}_1^o{{ {Q}}_\alpha }{ {B}}_1^{o{\rm{T}}} + { {D}}_1^o{{ {Q}}_\beta }{ {D}}_1^{o{\rm{T}}}{\rm{)}}^{ - 1}}$。但这种构造方式需要获取各传感器位置和速度真值、各传感器相对辐射源的方位角/高低角真值、各传感器到辐射源距离与距离变化率真值，并且还需已知辐射源位置与速度等真值。然而这些量在实际工程实践中不可得。

解决上述矛盾的一种折中的处理方式是采取两步加权最小二乘计算${ {\varphi }}_1^o$，即：

① 令权矩阵${ {W}} = { {Q}}_\beta ^{ - 1}$，则有：

那么，$ {{\hat{ u}}^o} \!=\! {{\hat{ {\varphi} }}_1^o}\left( {1:3} \right),\;\;{{{\hat {\dot { u}}}^o}} = {\hat{ {\varphi} }}_1^o\left( {5:7} \right),\;\;{{\hat r}_i^o} = \left\| {{{{\hat{ u}}}^o} - {{ {s}}_i}} \right\|,\;\; $${{\hat {\dot r}}_i^o} = {{{({{{\hat{ u}}}^o} - {{{ {s}}_i})}^{\rm T}} }({{{{\hat {\dot { u}}}}^o} - {{{\dot{ s}}}_i})} / {{\hat r}_i^o}}$。将${{\hat{ u}}^o}$与${{ {s}}_i}$代入式(1)可得${\hat \theta} _i^o$与$\hat \phi _i^o$。将上述各真值的替代值代入式(16)、式(17)与式(18)构造矩阵${\hat{ B}}_{\theta \phi }^o,{\hat{ D}}_{\theta \phi }^o,{{\hat{ B}}^o},{{\hat {\dot { B}}}^o},{{\hat{ D}}^o},{{\hat {\dot { D}}}^o}$，进而可构造出矩阵${\hat{ B}}_1^o$与${\hat{ D}}_1^o$。

② 选取权${{ {W}}_1} = {{\rm{(}}{\hat{ B}}_1^o{{ {Q}}_\alpha }{\hat{ B}}_1^{o{\rm{T}}} + {\hat{ D}}_1^o{{ {Q}}_\beta }{\hat{ D}}_1^{o{\rm{T}}}{\rm{)}}^{ - 1}}$，则${ {\varphi }}_1^o$的估计值${{ {\varphi }}_1}$为：

根据最小二乘理论，可知第一阶段参数的估计误差协方差阵为${\rm{ cov}} {\rm{(}}{{ {\varphi }}_1}{\rm{)}} = {{\rm{(}}{ {G}}_1^{\rm{T}}{{ {W}}_1}{{ {G}}_1}{\rm{)}}^{ - 1}}$。

2)第二阶段：精估计。

对于第一阶段的参数估计值${{ {\varphi }}_1}$，其前三维分量${{ {\varphi }}_1}(1:3)$是移动辐射源在某时刻的位置向量${{ {u}}^o}$的估计值，其第五至第七维分量${{ {\varphi }}_1}(5:7)$是该辐射源在该时刻的瞬时速度向量${{\dot{ u}}^o}$的估计值。但根据最小二乘理论，由式(20)估计得出${{ {\varphi }}_1}$的前提条件是$(r_1^o,{\dot r}_1^o)$与$({{ {u}}^o},{{\dot{ u}}^o})$无关，但这又与第一阶段估计参数${ {\varphi }}_1^o$的定义矛盾，因而按照这种方式估计出的移动辐射源的位置向量与速度向量不精确。因此，在第二阶段，本文使用第一阶段的参数估计值${{ {\varphi }}_1}(4)$与${{ {\varphi }}_1}(8)$来构建另一个最小二乘模型，以提高${{ {u}}^o}$与${{\dot{ u}}^o}$估计精度，从而实现对该时刻移动辐射源的无源定位。

记向量${{ {\varphi }}_1}$对应的估计误差向量为$\Delta {{ {\varphi }}_1}$，由于${{ {\varphi }}_1}(1:3) = {{ {u}}^o} + \Delta {{ {\varphi }}_1}(1:3)$，则有：

此外，从${ {\varphi }}_1^o$的定义可以得知${{ {\varphi }}_1}(4) = r_1^o + \Delta {{ {\varphi }}_1}(4) =$$ \left\| {{{ {u}}^o} - { {s}}_1^o} \right\| + \Delta {{ {\varphi }}_1}(4)$，将${ {s}}_1^o = {{ {s}}_1} - \Delta {{ {s}}_1}$代入上式，并忽略二阶误差项，则有：

类似地，有：

以及：

基于上面的讨论，第二阶段参数向量可选择为：

联合式(21)～式(24)，有：

式中，

另一方面，针对等式误差${{ {\varepsilon }}_2}$，有：

式中，矩阵${{ {B}}_2},{{ {D}}_2}$的定义如下所示：

式(28)决定了第二阶段加权最小二乘的权矩阵${{ {W}}_2}$的选取形式^[10]，即${ {W}}_2^{ - 1} = \{ {{ {B}}_2}{\rm{cov}}({{ {\varphi }}_1}){ {B}}_2^{\rm T} + {{ {D}}_2}{{ {Q}}_\beta }{ {D}}_2^{\rm T} + $${{ {B}}_2}{\rm{cov}}({{ {\varphi }}_1}){ {G}}_1^{\rm T} {{ {W}}_1}{{ {D}}_1}{{ {Q}}_\beta }{ {D}}_2^{\rm T} + {{ {D}}_2}{{ {Q}}_\beta }{ {D}}_1^{\rm T} {{ {W}}_1}{{ {G}}_1}{\rm{cov}}({{ {\varphi }}_1}){ {B}}_2^{\rm T} \} $。那么第二阶段估计参数${ {\varphi }}_2^o$的加权最小二乘估计值为：

且${ {\varphi }}_2^o$的估计误差为：

对应的协方差矩阵为${\rm{cov}} {\rm{(}}{{ {\varphi }}_2}{\rm{)}} = {{\rm{(}}{ {G}}_2^{\rm{T}}{{ {W}}_2}{{ {G}}_2}{\rm{)}}^{ - 1}}$。

一旦得到${{ {\varphi }}_2}$，那么位置向量${{ {u}}^o}$及速度向量${{\dot{ u}}^o}$的最终估计值可由以下方式计算得到，即：

式中，${ {\varPi }} = {\rm{diag}}[{{ {\varphi }}_1}(1:3) - {{ {s}}_1}]$；“./”表示对应元素相除。矩阵${ {\varPi }}$的存在避免了对${{ {\varphi }}_2}$元素开平方时的正负符号模糊问题。

式(32)最终给出了某时刻移动辐射源AOA-TDOA-FDOA的联合定位闭合解。下面给出定位解${ {u}}$与${\dot{ u}}$对应的协方差矩阵，若记算法给出的定位解向量${ {\varTheta }} = {{\rm{[}}{{ {u}}^{\rm{T}}},{{\dot{ u}}^{\rm{T}}}{\rm{]}}^{\rm{T}}}$，那么对式(25)求一阶全微分有：

式中，${ {\varXi }} = \left[ {\begin{array}{*{20}{c}} {2{\rm{diag}}[{{ {u}}^o} - {{ {s}}_1}]}&{{{ {{ O}}}_{3 \times 3}}} \\ {2{\rm{diag}}[{{{\dot{ u}}}^o} - {{{\dot{ s}}}_1}]}&{2{\rm{diag}}[{{ {u}}^o} - {{ {s}}_1}]} \end{array}} \right]$。

在第二阶段的估计过程中，${{ {s}}_1}$与${{\dot{ s}}_1}$被视作常量，因而$\Delta {{ {s}}_1}$与$\Delta {{\dot{ s}}_1}$未出现在式(33)里。式(28)表明当噪声强度不大时，${{ {\varepsilon }}_2}$近似满足零均值，因此式(31)所示的$\Delta {{ {\varphi }}_2}$同样近似为零均值。所以，本文所提的AOA-TDOA-FDOA联合定位给出的闭合解同样近似满足零均值特性。那么，由式(33)可得定位解向量${ {\varTheta }}$的理论协方差矩阵为：

需要指出的是，文献[10-11]分析了基于TDOA与FDOA观测量的联合定位闭合解算法，其第一阶段估计参数向量与本文选取的向量${ {\varphi }}_1^o$一致。${ {\varphi }}_1^o$是8维向量，因此文献[10-11]所提算法至少需要5个协同平台参与定位任务，才能按照式(19)及式(20)估计${ {\varphi }}_1^o$。而在实际的电子对抗作战过程中难以满足这一条件，本文在继承文献[10-11]所提算法框架的前提下，将AOA观测量引入定位算法，这就使得实际作战所需的协同平台个数降至3个，既节约了定位成本，也提高了无源定位闭合形式算法的战场生存能力。

上节详细讨论了某时刻移动辐射源AOA-TDOA-FDOA的联合定位闭合解，本节将对上述算法进行理论分析，并给出其估计误差的理论误差下界，即CRLB界。

AOA、TDOA以及FDOA测量值向量${ {\alpha}} $与传感器位置和速度向量${ {\beta }}$均服从高斯分布，并且两者独立，因此数据向量${ {\mu }} = {{\rm{[}}{{ {\alpha }}^{\rm{T}}},{{ {\beta }}^{\rm{T}}}{\rm{]}}^{\rm{T}}}$与参数向量${ {\vartheta}} = {[{{ {u}}^{\rm{T}}},{{\dot{ u}}^{\rm{T}}},{{ {s}}^{\rm{T}}},{{\dot{ s}}^{\rm{T}}}]^{\rm{T}}}$的概率密度函数为：

式中，${k_1} = - {{\ln \{ {{(2{ {{\text π}}} )}^{4M - 2}}\left| {{{ {Q}}_a}} \right|\} } / 2};{k_2} = - {{\ln \{ {{(2{ {{\text π}}} )}^{6M}}\left| {{{ {Q}}_{{\beta }}}} \right|\} } / 2},$都是常数。那么根据定义，参数${ {\theta}} $的CRLB界为：

式中，

由向量求导定义可知，式(37)中的偏导${\partial {{ {\alpha}} ^o}} / $${\partial { {\varTheta}} }$与${{\partial {{ {\alpha}} ^o}} / {\partial { {\beta}} }}$实质是向量${{ {\alpha}} ^o}$对向量${ {u}}$、向量${\dot{ u}}$、向量${ {s}}$及向量${\dot{ s}}$求偏导，即有：

以及：

下面逐一给出式(38)、(39)中每个偏导矩阵的具体表达式。根据定义有${{\partial {{ {\theta}} ^o}} / {\partial {\dot{ u}}}} = {{ O}_{M \times 3}}, {{\partial {{ {\phi}} ^o}} / {\partial {\dot{ u}}}} =$ $ {{ O}_{M \times 3}}$以及${{\partial {{ {\theta}} ^o}}/ {\partial {\dot{ s}}}} = {{ O}_{M \times 3M}},{{\partial {{ {\phi}} ^o}}/ {\partial {\dot{ s}}}} = {{ O}_{M \times 3M}}$。记${{ {\lambda }}_i}, {{ {\kappa }}_i}, $ ${\rm{ (}}i = 1,2, \cdots ,\;M{\rm{)}}$，具有如下表达形式：

那么，$ {{\partial {{ {\theta}} ^o}} / {\partial { {u}}}} = {{\rm{[}}{{ {\lambda }}_1},{{ {\lambda }}_2}, \cdots ,{{ {\lambda }}_M}{\rm{]}}^{\rm{T}}}\dot = {{ {A}}_{u1}},{{\partial {{ {\phi}} ^o}} / {\partial { {u}}}} = $${\rm{[}}{{ {\kappa }}_1}, {{ {\kappa }}_2}, \cdots , {{ {\kappa }}_M}{\rm{]}}^{\rm{T}}\;\dot =\;{{ {A}}_{u2}},{{\partial {{ {\theta}} ^o}} / {\partial { {s}}}} = {\rm{diag}}{\rm{(}}{{ {A}}_{u1}}{\rm{)}}\;\dot =\; {{ {A}}_{s1}},{{\partial {{ {\phi}} ^o}} / {\partial { {s}}}} =$${\rm{diag}}{\rm{(}}{{ {A}}_{u2}}{\rm{)}}\dot = {{ {A}}_{s2}} $。

记${{ {E}}_{3M \times 3}} = {\left[ {{{ {I}}_{3 \times 3}}, \cdots ,{{ {I}}_{3 \times 3}}} \right]^{\rm{T}}}$，并定义两个(M−1)×3M的矩阵${{ {C}}_1}$与矩阵${{ {C}}_2}$，其第$i{\rm{(}}i = 1,2, \cdots ,M - 1{\rm{)}}$行的定义为：

式中，$ {{ {a}}_i} \!=\! ({{ {u}}^o} - { {s}}_i^o)/r_i^o;{{ {b}}_i} \!=\! ({{\dot{ u}}^o} - {\dot{ s}}_i^o)/r_i^o - {\dot r}_i^o({{ {u}}^o} - { {s}}_i^o)/r_i^{o2}$。那么，则有：

综合上述各个偏导矩阵，那么式(38)与式(39)所表征的偏导矩阵的具体表达形式为：

将式(43)代入式(37)，即可得6M+6维方阵${\rm{CRLB}}\left( { {\vartheta }} \right)$的具体形式。

由矩阵求逆引理可得移动辐射源位置及速度参数${ {\varTheta }}$的CRLB界，即：

式中，${{{X}}^{ - 1}}$表示当测角与时差/频差传感器不存在位置误差及速度误差时的参数${ {\varTheta }}$的CRLB界，因此式(44)中的第二项表征了由测角与时差/频差传感器位置及速度误差引起的参数${ {\varTheta }}$的CRLB界的增量部分。式(44)的迹表征了任一种多平台AOA-TDOA-FDOA联合定位算法的估计误差均方根的下界。

文献[10]已证明，当TDOA与FDOA测量误差与传感器位置与速度误差都服从高斯分布并且噪声不大，并且传感器距离辐射源较远时，定位解向量${ {\varTheta }}$的理论协方差矩阵${\rm{cov}}({ {\varTheta }})$能达到其理论误差下界。基于类似的证明过程，当无源定位算法加入AOA观测量后，本文所提算法所给出定位解向量${ {\varTheta }}$的理论协方差矩阵${\rm{cov}}({ {\varTheta }})$也能达到对应的CRLB界。由于篇幅限制，在此省略这一结论的推导过程，下面将给出仿真实验，验证本文提出算法的均方根误差在噪声不大时同样能达到CRLB界。

为验证本文提出的移动辐射源AOA-TDOA- FDOA联合定位闭合解算法的性能，本节首先给出一组3平台协同定位一个移动辐射源位置的仿真实例。移动辐射源初始位置为[57 891, −12 536, $ 125{\rm{]}}^{\rm{T}}$(单位为m，下同)，初始速度为${{\rm{[}}12,30,0{\rm{]}}^{\rm{T}}}$(单位为m/s，下同)；协同平台主站初始位置为[5 000，$1\;500, \,3\;500{\rm{]}}^{\rm{T}}$，初始速度为${{\rm{[}}120, - 20,0{\rm{]}}^{\rm{T}}}$；协同平台辅站1的初始位置为${{\rm{[}}21\;500,\;15\;000,5\;800{\rm{]}}^{\rm{T}}}$，初始速度为${{\rm{[}}130,0,0{\rm{]}}^{\rm{T}}}$；协同平台辅站2的初始位置为${\rm{[}}20\;500, $$- 45\;000,5\;800{\rm{]}}^{\rm{T}}$，初始速度为${{\rm{[}}110,50,0{\rm{]}}^{\rm{T}}}$。辐射源载波频率为$1.125 \times {10^9}$ Hz。为评价无源定位算法的定位精度，本文采用均方根误差为评价标准，其定义为：

式中，${\rm{MSE}}$为某时刻无源定位所估计的辐射源位置的均方根误差；${{ {u}}^o}$为该时刻移动辐射源真实位置；${{ {u}}_i}$为第i次独立仿真实验所估计的辐射源位置；N为独立仿真实验次数。需要指出的是，在N次独立仿真实验过程中，每次仿真中辐射源移动轨迹以及每个平台的真实运动轨迹都一致。本节后续部分将首先以此为基础给出几组仿真实例，验证本文所提算法在各种噪声强度下的定位能力；然后再给出两组对比实验，说明本文所提算法在定位精度与定位时间等方面较其他已有定位算法所具备的优势。

为验证本文所提算法对平台位置误差的适应能力，在这组仿真中将3架空中平台的图1中x、y、z向上的位置误差标准差变化范围同时设置为10⁻²～10⁵ m，所有平台图1中x、y、z向上的速度误差标准差固定为0.2 m/s，并且将所有AOA测量噪声标准差均固定为3°，所有TDOA测量噪声标准差均固定为100 ns，以及所有FDOA测量噪声标准差均固定为100 Hz。5 000次独立仿真形成的均方根定位误差曲线如图2所示。为直观说明本文所提算法对平台位置误差的适应能力，每种设置下的定位误差理论下界(即式(44)所示的${\rm{CRLB}}\left( { {\varTheta }} \right)$矩阵的前3个对角元素和的平方根)也同时绘制在图2中。

图  2  平台位置误差与定位精度关系曲线

图2表明当所有协同平台在图1中x、y、z这3个方向的位置误差标准差均不超过5×10³ m时，本文所提算法在3平台协同定位场合下的均方根定位误差能达到对应的CRLB界。

为验证本文所提算法对测角噪声强度的适应能力，在这组仿真中，将每个平台的方位角测量噪声标准差变化范围同时设置为10⁻²°～10²°，而所有平台在图1中x、y、z向的位置误差标准差都固定为5 m、速度误差标准差都固定为0.2 m/s，并且将所有高低角测量噪声标准差均固定为3°，所有TDOA测量噪声标准差均固定为100 ns，所有FDOA测量噪声标准差均固定为100 Hz。5 000次独立仿真形成的均方根定位误差曲线如图3所示。同时为直观说明本文所提算法对测角噪声强度的适应能力，每种设置下的CRLB界也同时绘制在图3中。

图  3  方位角误差与定位精度关系曲线

图3表明当所有协同平台所测得的方位角的测量误差标准差不大于20°时，本文所提算法在3平台协同定位场合下的均方根定位误差能达到对应的CRLB界。

为验证本文所提算法对TDOA噪声强度的适应能力，在这组仿真中，将每个TDOA测量噪声标准差变化范围同时设置为10⁻²～10⁶ ns，所有平台在x、y、z向的位置误差标准差都固定为5 m，速度误差标准差都固定为0.2 m/s，并且将所有方位/高低角测量噪声标准差均固定为3°，所有FDOA测量噪声标准差均固定为100 Hz。5 000次独立仿真形成的均方根定位误差曲线如图4所示。同时为直观说明本文所提算法对TDOA测量噪声强度的适应能力，每种设置下的CRLB界也同时绘制在图4中。

图  4  TDOA误差与定位精度关系曲线

图4表明当所有TDOA观测量的测量噪声标准差不大于10⁵ ns时，本文所提算法在3平台协同定位场合下的均方根定位误差能达到对应的CRLB界。

为验证本文所提算法对FDOA噪声强度的适应能力，在这组仿真中，将每个FDOA测量噪声标准差变化范围同时设置为10⁻²～10⁴ Hz，所有平台在x、y、z向的位置误差标准差都固定为5 m，速度误差标准差都固定为0.2 m/s，并且将所有方位/高低角测量噪声标准差均固定为3°，所有TDOA测量噪声标准差均固定为100 ns。5 000次独立仿真形成的均方根定位误差曲线如图5所示。同时为直观说明本文所提算法对FDOA测量噪声强度的适应能力，每种设置下的CRLB界也同时绘制在图5中。

图  5  FDOA误差与定位精度关系曲线

图5表明当所有FDOA观测量的测量噪声标准差不大于200 Hz时，本文所提算法在3平台协同定位场合下的均方根定位误差能达到对应的CRLB界。

针对移动辐射源TDOA-FDOA联合无源定位闭合解算法需要至少5个协同平台参与一次定位任务的问题，本文将AOA观测量引入，提出一种新的基于AOA-TDOA-FDOA的联合定位闭合解算法。本文所提算法在兼顾平台位置不确定性以及保持代数闭合解框架不变的前提下，能将参与定位的平台个数有效地减少至3个。同时，理论分析与仿真实验结果都表明所提算法在面对工程上常见的噪声强度时，定位性能可以达到理论误差下界。这表明，本文所提算法既能拓宽算法的战场应用范围，也能提高算法的战场生存能力。