1 引言

近年来，相比传统体制的雷达，MIMO雷达以其在检测，跟踪和参数估计更优的性能获得了大量关注。MIMO雷达分为共置天线 ^[1] 和分布式两种体制 ^[2] 。本文主要研究的是分布式MIMO雷达多目标跟踪背景下的发射站，接收站优化选择。分布式MIMO雷达属于无线传感网络的一种特殊形式。在传感器网络中，由于频谱资源，通信资源，能量资源的限制，需要通过对传感器资源进行合理的优化配置来达到更高的效率。分布式MIMO雷达也提出了对收发传感器进行优化选择的需求。

传感器选择问题作为传感器网络管理的一个重要子课题，涌现出了大量有价值的成果。文献[3– 5]提出了一系列关于一般传感器优化选择的方法。文献[3]最早提出了一个通用的凸优化选择模型来进行传感器选择。在该文中，传感器选择被建立为一个尽可能小化估计误差椭圆体积情况下的布尔变量的选择问题。通过该模型可以在满足固定传感器数量需求的情况下实现参数估计性能的性能最优选择。文献[4]通过对卡尔曼滤波器的广义推导，利用半正定松弛技术，实现了在各个传感器之间存在相关噪声下的优化选择，但是其优化模型忽略了当前时刻费歇尔信息矩阵的先验信息，所以其性能选择准则是非最优的。文献[5]用扩展卡尔曼滤波器中的卡尔曼增益的列稀疏性来表征传感器选择的有无，实现了跟踪性能和被选传感器数量的折中。

针对分布式MIMO雷达体制下的发射站接收站选择，Hana Godrich根据现有的传感器选择方法较早的关注了在分布式MIMO雷达体制下的收发站选择问题。她在文献[6]中提出了一种贪婪迭代算法，即每一次迭代中选择性能最好的传感器，在较低计算量的前提下获得了发射站和接收站的近似最优选择。然而这种方法更多的是基于对分布式MIMO雷达的性能分析下的一种直觉选择。在文献[7]中，在MIMO雷达网络在定位的应用背景下，实现了一种基于2阶椎规划(Second Order Cone Programming, SOCP)的分两步走发射站选择和功率分配方法，即在先选择发射站后再对发射功率进行分配。该方法虽然寻求联合资源分配，但是一方面优化求解过程不是联合进行，另一方面没有考虑接收站选择在MIMO雷达网络中的影响。在MIMO雷达目标跟踪应用背景下，利用后验克拉美罗下界(Posterior Cramer-Rao Lower Bound, PCRLB)作为性能指标，文献[8]中把优化问题设计为尽可能少的获取选择发射站的数量同时尽可能高的提高跟踪性能。该模型比较符合实际任务需求，但是由于其目标函数的非凸性，只能采取启发式方法进行求解，其有效性无法总能得到保证，两种代价函数之间的优化系数也没有给出一种较好的求解方案。文献[9]考虑了MIMO雷达网络在分散发射聚焦接收工作模式下的多目标接收站选择。虽然文中给出了基于凸优化的选择方案，但是文中假设发射站雷达的全向发射模式现代雷达还难以达到。

在现有的传感器选择方法中。对于一般通用的传感器选择算法，在很大程度上因为具体模型的不同，无法直接应用于分布式MIMO的传感器选择。另外一方面，现有的分布式MIMO雷达选择算法要么采用的是一种启发式的算法，要么只能选择雷达网络中发射站和接收站中的其中一种，无法保证其最优性。本文将优化多目标中跟踪性能最差的目标后验克拉美罗下界为目标函数，构建关于收发站的代价函数，利用分块坐标迭代，实现了多目标跟踪背景下分布式MIMO雷达收发站的联合选择。仿真实验，验证了本方法相比其它算法获得了更优的选择性能，并且在拥有更低计算量的情况下逼近全局搜索的最优选择。

2 系统信号模型

考虑这样一个拥有 M个发射站， N个接收站非相参分布式MIMO雷达系统。在雷达监视范围内一共存在 Q个目标。在 k时刻，第 q个目标的状态向量可以表示为 ${{x}}_k^q = {\left[ {x_k^q,\dot x_k^q,y_k^q,\dot y_k^q} \right]^{\rm{T}}}$ , $q \in {\cal Q}$ 。其中 ${\cal Q}$ 定义为目标集合， ${\left( \cdot \right)^{\rm{T}}}$ 表示求转置运算。第 m个发射站和第 n个接收站的位置分别定义为 $\left( {{x_m},{y_m}} \right)$ ， $m \in {\cal M}$ 和 $\left( {{x_n},{y_n}} \right)$ ， $n \in {\cal N}$ 。其中 ${\cal M}$ 是发射站集合， ${\cal N}$ 是接收站集合。对于第 q个目标，第 $\left( {m,n} \right)$ 条路径所对应的延时和多普勒频移可以分别表示为：

$\tau _{m,n,k}^q = \frac{{R_{m,k}^q + R_{n,k}^q}}{\rm c} \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad $

(1)

$\begin{array}{l}f_{m,n,k}^q = \frac{{{f_{\rm{c}}}}}{{\rm{c}}}\left[ {\frac{{\dot x_k^q\left( {{x_{m,k}} - x_k^q} \right) + \dot y_k^q\left( {{y_{m,k}} - y_k^q} \right)}}{{R_{m,k}^q}}} \right.\\ \quad \quad \quad \quad \left. { + \frac{{\dot x_k^q\left( {{x_{n,k}} - x_k^q} \right) + \dot y_k^q\left( {{y_{n,k}} - y_k^q} \right)}}{{R_{n,k}^q}}} \right]\end{array}$

(2)

其中， $R_{m,k}^q = \sqrt {{{\left( {x_k^q - {x_{m,k}}} \right)}^2} + {{\left( {y_k^q - {y_{m,k}}} \right)}^2}} $ 为第 q个目标到第 m个发射站的距离，而 $R_{n,k}^q \!=\! $ $ \sqrt {{{\left( {x_k^q - {x_{n,k}}} \right)}^2} + {{\left( {y_k^q - {y_{n,k}}} \right)}^2}} $ 表示其到第 n个发射站的距离。c表示光速。

发射站发射一系列能量归一化的正交波形满足 $\int_{ - \infty }^{ + \infty } {{{\left| {{{\tilde s}_m}\left( t \right)} \right|}^2}{\rm d}t = 1} $ ，并且对于 $\forall m \ne n$ ， $ \displaystyle\int_{ - \infty }^{ + \infty } {{{\tilde s}_m}\left( t \right)\tilde s_n^*\left( {t - \tau } \right){{\rm e}^{ -{\rm j}2{{π}} ft}}{\rm d}t = 0\;\;} $ 。第 q个目标被第 m个发射站照射并且被第 n个接收站所接收到的基带信号表示为：

$\tilde r_{m,n}^q\left( t \right) \!=\! \tilde a_{m,n}^qh_{m,n,k}^q{\tilde s_m}\left( \!{t \!-\! \tau _{m,n,k}^q} \!\right){{\rm e}^{{\rm j}2{{π}} f_{m,n,k}^qt}} \!+\! {\tilde n_{m,n}}\left( t \right)$

(3)

其中，目标的RCS起伏被建模为一个复高斯随机变量 $\tilde a_{m,n}^q \sim {\cal C}{\cal N}\left( {0,2\sigma _{m,n}^2} \right)$ 。双站之间的路径衰落和发射站信号发射功率用 $h_{m,n}^q \propto {\left[ {{P_m}/{{\left( {R_{m,k}^qR_{n,k}^q} \right)}^2}} \right]^{1/2}}$ 来表征。 ${\tilde n_{m,n}}\left( t \right)$ 是一个循环对称零均值方差为 N ₀的复高斯随机过程。定义 ${\phi _{m,k}} = {\tan ^{ - 1}}\left( {\displaystyle\frac{{y_k^q - {y_{m,k}}}}{{x_k^q - {x_{m,k}}}}} \right)$ 为目标 q相对发射站 m的方位角， ${\varphi _{n,k}} \!= \! {\tan ^{ - 1}} \! \!\left( \! {\displaystyle\frac{{y_k^q \!- \! {y_{n,k}}}}{{x_k^q \!- \! {x_{n,k}}}}} \!\right)$ 为目标 q相对于接收站 n的方位角。

3 目标动态模型和量测模型

假设目标运动模型被建模为随机加速度线性运动模型，那么对于第 q个目标可以表示为：

${{x}}_k^q = {{Fx}}_{k - 1}^q + {{v}}_{k - 1}^q$

(4)

其中， ${{x}}_k^q = {\left[ {x_k^q,\dot x_k^q,y_k^q,\dot y_k^q} \right]^{\rm{T}}}$ 是关于目标状态的4维向量。 ${{v}}_k^q$ 是目标运动的过程噪声，满足0均值高斯分布，其对应的协方差矩阵为：

${{Q}} = \xi {{{I}}_2} \otimes \left[ {\begin{array}{*{20}{c}}{\frac{{\Delta {T\;^3}}}{3}} & {\frac{{\Delta {T\;^2}}}{2}}\\{\frac{{\Delta {T\;^2}}}{2}} & {\Delta T}\end{array}} \right]$

(5)

$\Delta T \,$ 是目标运动的采样间隔， $\xi $ 表示过程噪声的强度系数，而 ${{{I}}_2}$ 是一个维度为2的单位矩阵。 $ \otimes $ 表示kronecker积。目标运动的状态转移矩阵可以表示为：

${{F}} = {{{I}}_2} \otimes \left[ {\begin{array}{*{20}{c}}1 & {\Delta T}\\0 & 1\end{array}} \right]$

(6)

在一般的文献中，通常会将每个接收站获得的全部采样数据传输至数据融合中心来统一进行信号级数据融合以获得目标的位置速度最优估计。但是为了降低数据传输率和时频不同步等其它问题，本文采用的是分布式非直接估计方法 ^[10] ，即每个接收站都对目标对应 M条反射波路径的延时和多普勒频移进行最大似然估计然后再传输至数据融合中心，而不必传输全部采样数据。在非直接参数估计中，本文定义一组关于目标的未知参数 $ {{θ}}\! \left( {{{x}}_k^q} \right) \!=\!\left[\! { {{\left( {{{θ}} _{1,1,k}^q} \right)}^{\rm{T}}}, }\right.\cdots ,{\left. {{{\left( {{{θ}} _{M,N,k}^q} \right)}^{\rm{T}}}} \right]^{\rm{T}}} $ 其中 ${{θ}} _{m,n,k}^q = {\left[ {\tau _{m,n,k}^q,\omega _{m,n,k}^q} \right]^{\rm{T}}}$ 表示第 q个目标在第 mn条路径下的待估延时和多普勒频。它可以在第 n个接收站通过最大似然估计得到

$\hat {{θ}} _{m,n,k}^q = \arg \mathop {\max }\limits_{\hat \tau _{m,n,k}^q,\hat f_{m,n,k}^q} \sum\limits_{m = 1}^M \int {\Big| {\tilde r_n^q\left( t \right)\tilde s_m^*\left( {t - \hat \tau _{m,n,k}^q} \right)}\Big. }\\ \quad \quad \quad \ \ \, {\Big. { {{\rm{e}}^{ - {\rm{j}}2{{π}}\hat f_{m,n,k}^qt}}} \Big|^2}{\rm{dt}}$

(7)

其中， $\tilde r_n^q\left( t \right) = \sum\limits_{m = 1}^M {\tilde r_{m,n}^q\left( t \right)} $ 是第 n个发射站的回波信号之和。于是非线性量测方程为：

${{y}}_k^q = {{{h}}_k}\left( {{{x}}_k^q} \right) + {{w}}_k^q$

(8)

其中， ${{{h}}_k}\left( \cdot \right)$ 表示目标的非线性量测方程， ${{w}}_k^q$ 为量测噪声。 ${{y}}_k^q$ 即为各条路径下对目标延时多普勒频率的量测值。在大多数应用背景下，由于在高信噪比条件下，对目标未知参数 ${{θ}} _k^q$ 的最大似然估计总是能到达其对应的克拉美罗下界(Cramer-Rao Lower Bound, CRLB)，所以对未知参数 ${{θ}} _k^q$ 的量测噪声矩阵可以表示为其费歇尔信息矩阵的逆

${{J}}_{\rm{Z}}^{ - 1}\!\left( {{{θ}} _k^q} \right) \!\!=\!\! {\left[ {{\rm{blkdiag}}\left(\! {{{{J}}_{\rm{Z}}}\left( {{{θ}} _{1,1,k}^q} \right),\! \cdots\! ,{{{J}}_{\rm{Z}}}\left( {{{θ}} _{M,N,k}^q} \right)} \right)} \right]^{ - 1}}$

(9)

其中， ${\rm{blkdiag}}\left( \cdot \right)$ 表示矩阵分块对角。由各条反射波路径的独立性可知，关于未知参数 ${{θ}} _k^q$ 的费歇尔信息矩阵可以表示为矩阵分块对角形式。第 q个目标在第 mn条路径下未知参数 ${{θ}} _{m,n,k}^q$ 的费歇尔信息矩阵在文献[11]中被给出

${{{J}}_{\rm{Z}}}\left( {{{θ}} _{m,n,k}^q} \right) = \frac{{8{{{π}} ^2}{{\left| {a_{m,n}^q} \right|}^2}h_{m,n,k}^q}}{{{N_0}}}\left[ {\begin{array}{*{20}{c}}{\beta _m^2} & {{\xi _m}}\\{{\xi _m}} & {T_m^2}\end{array}} \right]$

(10)

其中， ${\beta _m} = {\left[ {\frac{\displaystyle{\int {{f^2}{{\left| {{S_m}\left( f \right)} \right|}^2}{\rm d}\,f} }}{\displaystyle{\int {{{\left| {{S_m}\left( f \right)} \right|}^2}{\rm d}\,f} }}} \right]^{1/2}}$ 为信号等效带宽， ${T_m} = {\left[ {\frac{\displaystyle{\int {{t^2}{{\left| {{s_m}\left( t \right)} \right|}^2}{\rm d}t} }}{\displaystyle{\int {{{\left| {{s_m}\left( t \right)} \right|}^2}{\rm d}t} }}} \right]^{1/2}}$ 为信号等效时宽，并定义 ${\xi _m}$ 为 ${\xi _m} = {\mathop{\rm Im}\nolimits} \left[ {\displaystyle\int {t{s_m}\left( t \right)\frac{{\partial s_m^*\left( t \right)}}{{\partial t}}{\rm d}t} } \right]$ 。

4 后验克拉罗美罗下界下的跟踪性能

后验克拉罗美罗下界(PCLRB)是常用的用来衡量目标跟踪性能的指标 ^[12] 。对于一个任意未知参数向量 a ，通过观测向量 b ，它的后验克拉罗美罗下界(PCLRB)可以表达为：

${{\mathbb{E}}_{{{a,b}}}}\left\{ {\left[ {{\overset{\frown} {{a}}} \left( {{b}} \right) - {{a}}} \right] {{\left[ {{\overset{\frown} {{a}}} \left( {{b}} \right) - {{a}}} \right]}^{\rm{T}}}} \right\} \succeq {{J}}_{\rm{B}}^{ - 1}$

(11)

其中， ${{{J}}_{\rm{B}}}$ 是贝叶斯信息矩阵(BFIM)， $ \succeq $ 表示广义不等式。 ${\mathbb{E}}\left( \cdot \right)$ 表示期望求取符号。BFIM可以通过联合概率密度函数 ${p_{{{a,b}}}}\left( {{{a,b}}} \right)$ 来计算

$ {{{J}}_{\rm{B}}} = - {\mathbb{E}}\left[ {\Delta _{{a}}^{{a}}\log {p_{{{a}}, {{b}}}}\left( {{a}, {{b }}} \right)} \right] $

(12)

其中， $\Delta _{{b}}^{{a}} = {\nabla _{{a}}}\left[ {\nabla _{{b}}^{\rm{T}}} \right]$ ，而 $$\nabla $ 是梯度算子。根据 ${p_{{{a}}, {{b}}}}\left( {{{a,b}}} \right) = {p_{{{b|a}}}}\left( {{{b|a}}} \right){p_{{a}}}\left( {{a}} \right)$ 可以得到

${{{J}}_{\rm{B}}} = {{{J}}_{\rm{Z}}} + {{{J}}_{\rm{P}}}$

(13)

其中， ${{{J}}_{\rm P}}$ 是和先验信息有关的费歇尔信息矩阵，而 ${{{J}}_{\rm Z}}$ 是关于目标状态和量测数据有关的费歇尔信息矩阵分别定义为：

$\left. \begin{array}{l}{{{J}}_{\rm{Z}}} = - {{\mathbb{E}}_{{a}}}\left[ {\Delta _{{a}}^{{a}}{p_{{{b|a}}}}\left( {{{b|a}}} \right)} \right]\\{{{J}}_{\rm{P}}} = - {{\mathbb{E}}_{{a}}}\left[ {\Delta _{{a}}^{{a}}{p_{{a}}}\left( {{a}} \right)} \right] \end{array} \right\}$

(14)

在本文的目标跟踪模型中， ${\hat {{x}}}_k^q\left( {{{Y}}_{1:k}^q} \right)$ 代表了根据从1到 k时刻观测值所获得观测值对目标状态 $${{x}}_k^q$ 的估计。那么其估计误差协方差矩阵被PCRLB最小界定为：

${{P}}_k^q \!=\! {\mathbb{E}}\!\left\{ \!\!{\left[ {{\overset{\frown} {{x}}}_k^q\left( {{{Y}}_{1:k}^q} \right) \!-\! {{x}}_k^q} \right]\!\!{{\left[ {{\overset{\frown} {{x}}}_k^q\!\left( {{{Y}}_{1:k}^q} \right) \!-\! {{x}}_k^q} \right]}^{\rm{T}}}} \!\right\} \!\!\succeq\!\! {\left( {{{J}}_k^q} \right)^{ - 1}}$

(15)

要利用所有时刻的观测数据来计算PCRLB非常困难，文献[12]给出了BFIM ${{J}}_k^q$ 根据观测数据更新迭代计算形式

${{J}}_{k + 1}^q = {{D}}_k^{22} - {{D}}_k^{21}{\left( {{{J}}_k^q + {{D}}_k^{11}} \right)^{ - 1}}{{D}}_k^{12}$

(16)

这里为了符号简洁省略了BFIM的下标B来简化书写。式中各项可以分别表示为：

$\left. {\begin{array}{*{20}{l}} \ \;{{{D}}_k^{11} = {{\mathbb{E}}_{{{x}}_k^q,{{x}}_{k + 1}^q}}\left[ { - \Delta _{{{x}}_k^q}^{{{x}}_k^q}\ln p\left( {{{x}}_{k + 1}^q|{{x}}_k^q} \right)} \right]}\\ \ \;{{{D}}_k^{12} = {{D}}_k^{21} = {{\mathbb{E}}_{{{x}}_k^q,{{x}}_{k + 1}^q}}\left[ { - \Delta _{{{x}}_k^q}^{{{{x}}_{k + 1}}}\ln p\left( {{{x}}_{k + 1}^q|{{x}}_k^q} \right)} \right]}\\\begin{array}{l}{{D}}_k^{22} = {{D}}_k^{22}\left( 1 \right) + {{D}}_k^{22}\left( 2 \right) \\\ \;\ \;\ \ = {{\mathbb{E}}_{{{x}}_k^q,{{x}}_{k + 1}^q}}\left[ { - \Delta _{{{x}}_{k + 1}^q}^{{{x}}_{k + 1}^q}\ln p\left( {{{x}}_{k + 1}^q|{{x}}_k^q} \right)} \right]\\\ \;\ \; \ \ \ \ \ + {{\mathbb{E}}_{{{x}}_{k + 1}^q,{{y}}_{k + 1}^q}}\left[ { - \Delta _{{{x}}_{k + 1}^q}^{{{x}}_{k + 1}^q}\ln p\left( {{{y}}_{k + 1}^q|{{x}}_k^q} \right)} \right]\end{array}\end{array}} \!\!\!\! \right\}$

(17)

由于本文假设目标满足线性运动模型，通过矩阵求逆引理以及一系列的变换和化简，先验费歇尔信息矩阵可以表示为：

${{{J}}_{\rm{P}}}\left( {{{x}}_k^q} \right) = {\left[ {{{{Q}}^q} + {{F}}{{\left( {{{J}}_{k - 1}^q} \right)}^{ - 1}}{{{F}}^{\rm{T}}}} \right]^{ - 1}}$

(18)

而关于量测信息的费歇尔信息矩阵为：

$\begin{array}{l}{{{J}}_{\rm{Z}}}\left( {{{x}}_k^q} \right) = {{\mathbb{E}}_{{{x}}_k^q,{{y}}_k^q}}\left[ { - \Delta _{{{{x}}_k}}^{{{{x}}_k}}\ln p\left( {{{y}}_k^q|{{x}}_k^q} \right)} \right]\\\;\;\;\;\;\;\;\;\;\;\;\;\; = {{\mathbb{E}}_{{{x}}_k^q}}\left\{ {{{\mathbb{E}}_{{{y}}_k^q|{{x}}_k^q}}\left[ { - \Delta _{{{x}}_k^q}^{{{x}}_k^q}\ln p\left( {{{y}}_k^q|{{x}}_k^q} \right)} \right]} \right\}\\\;\;\;\;\;\;\;\;\;\;\;\;\; = {{\mathbb{E}}_{{{x}}_k^q}}\left[ {{{\left( {{{H}}_k^q} \right)}^{\rm{T}}}{{{J}}_{\rm{Z}}}\left( {{{θ}} _k^q} \right){{H}}_k^q} \right]\end{array}$

(19)

其中 ${{H}}_k^q \!=\!\! \left[ \!{{\nabla _{{{x}}_k^q}}{{\left( {{{θ}} _k^q} \right)}^{\rm{T}}}}\! \right] \!\!=\!\! {\left[ {{{\left( {{{H}}_{1,1,k}^q} \right)}^{\rm{T}}},\! \cdots \!,{{\left( {{{H}}_{M,N,k}^q} \right)}^{\rm{T}}}} \right]^{\rm{T}}}$ 是目标量测参数 ${{θ}} _k^q$ 关于目标状态 ${{x}}_k^q$ 的雅克比矩阵。 ${\left( {{{H}}_{m,n,k}^q} \right)^{\rm{T}}} \!\!=\!\! \left[\! {{\nabla _{{{x}}_k^q}}{{\left( {{{θ}} _{m,n,k}^q} \right)}^{\rm{T}}}} \!\right] \!\!=\!\! \left[\!\! {\begin{array}{*{20}{c}} {{\nabla _{{{x}}_k^q}}\tau _{m,n,k}^q} & {{\nabla _{{{x}}_k^q}}f_{m,n,k}^q} \end{array}} \!\!\right],$ $m \in {\cal M},n \in {\cal N}$ 是对应 mn条路径量测参数 ${{θ}} _{m,n,k}^q$ 关于目标状态 ${{x}}_k^q$ 的雅克比矩阵。 ${\nabla _{{{x}}_k^q}}\tau _{m,n,k}^q$ 和 ${\nabla _{{{x}}_k^q}}f_{m,n,k}^q$ 又可以分别表示为：

${\nabla _{{{x}}_k^q}}\tau _{m,n,k}^q = \frac{1}{\rm c}\left[ {\begin{array}{*{20}{c}}{\cos \phi _{m,k}^q + \cos \varphi _{n,k}^q}\\{\sin \phi _{m,k}^q + \sin \varphi _{n,k}^q}\\0\\0\end{array}} \right]$

(20)

和

${\nabla _{{{x}}_k^q}}f_{m,n,k}^q = - \frac{{2{{π}} }}{\lambda }\left[ {\begin{array}{*{20}{c}}\begin{array}{l}\dot x_k^q\left( {\frac{1}{{R_{m,k}^q}} + \frac{1}{{R_{n,k}^q}}} \right)\\ - \frac{{\cos \phi _{m,k}^q\left( {\dot x_k^q\cos \phi _{m,k}^q + \dot y_k^q\sin \phi _{m,k}^q} \right)}}{{{{\left( {R_{m,k}^q} \right)}^2}}} - \frac{{\cos \varphi _{n,k}^q\left( {\dot x_k^q\cos \varphi _{n,k}^q + \dot y_k^q\sin \varphi _{n,k}^q} \right)}}{{{{\left( {R_{n,k}^q} \right)}^2}}}\end{array}\\\begin{array}{l}\dot y_k^q\left( {\frac{1}{{R_{m,k}^q}} + \frac{1}{{R_{n,k}^q}}} \right)\\ - \frac{{\sin \phi _{m,k}^q\left[ {\dot x_k^q\cos \phi _{m,k}^q + \dot y_k^q\sin \phi _{m,k}^q} \right]}}{{{{\left( {R_{m,k}^q} \right)}^2}}} - \frac{{\sin \varphi _{n,k}^q\left[ {\dot x_k^q\cos \varphi _{n,k}^q + \dot y_k^q\sin \varphi _{n,k}^q} \right]}}{{{{\left( {R_{n,k}^q} \right)}^2}}}\end{array}\\{\cos \phi _{m,k}^q + \cos \varphi _{n,k}^q}\\{\sin \phi _{m,k}^q + \sin \varphi _{n,k}^q}\end{array}} \right]$

(21)

通常来说 ${{J}}_{\rm{Z}}^q$ 中期望的计算需要通过蒙特卡洛方法来进行计算，也就是说从概率密度函数 $p\left( {{{{x}}_k}|{{{x}}_{k - 1}}} \right)$ 中采样 N _s 个样本 $\left\{ {{{x}}_{k|k - 1}^{\left( j \right)}} \right\}_{j = 1}^{N_s}$ 以相同的权重1/ N _s 加权求和。另一方面，如果目标运动的过程噪声很小，则可以去掉期望符号。在这种条件下由于采用预先规划的方式来进行发射站接收站选择，所以可以利用预测值 ${{\overset{\frown} {{x}}}_{k|k - 1}}$ 来进行计算。值得注意的是，本文假设线性运动模型是为了降低计算PCRLB的计算量，但是对于问题模型的建立和分析并无影响。如果是非线性运动模型，那么在计算贝叶斯费歇尔信息矩阵的时候，先验费歇尔信息矩阵就不再是 ${{{J}}_{{\rm{P}},k}} = {\left[ {{{{Q}}^q} + {{F}}{{\left( {{{{J}}_{k - 1}}} \right)}^{ - 1}}{{{F}}^{\rm{T}}}} \right]^{ - 1}}$ ，而需要根据利用蒙特卡洛积分对 ${{{D}}^{11}}$ , ${{{D}}^{12}}$ , ${{{D}}^{21}}$ 和 ${{{D}}^{22}}\left( 1 \right)$ 的期望进行计算。对于机动目标运动模型，因为不满足计算PCRLB所需的观测值和参数向量的联合概率密度的对数是2次可微的条件，需要寻求其它性能准则。

5 发射站接收站优选模型及求解

在本节中，分布式MIMO雷达多目标发射站接收站联合选择模型被建立。首先本文引入两组选择向量 ${{{t}}^{\rm{T}}} \!\!= \!\! \left[ {{{\left( {{{{t}}^1}} \right)}^{\rm{T}}}, \cdots ,{{\left( {{{{t}}^Q}} \right)}^{\rm{T}}}} \right]$ 和 ${{{r}}^{\rm{T}}} \!\!=\! \! \left[ {{{\left( {{{{r}}^1}} \right)}^{\rm{T}}}, \cdots ,{{\left( {{{{r}}^Q}} \right)}^{\rm{T}}}} \right]$ 。其中， ${\left( {{{{t}}^q}} \right)^{\rm{T}}} = \left[ {t_1^q, \cdots ,t_M^q} \right],t_m^q \in \left\{ {0,1} \right\}$ 和 ${\left( {{{{r}}^q}} \right)^{\rm{T}}} = $ $ \left[ {r_1^q, \cdots ,r_N^q} \right],r_n^q \in \left\{ {0,1} \right\},q \in {\cal Q}$ 分别表示第 q个目标的选择向量。如果 $t_m^q = 1$ 表示针对第 q个目标，第 m个发射站被选中，反之如果 $t_m^q = 0$ 则没有。类似的， $r_n^q = 1$ 和 $r_n^q = 0$ 则分别表示针对第 q个目标，第 n个接收站是否被选中。因为数据费歇尔信息矩阵 ${{{J}}_{\rm{Z}}}\left( {{{x}}_k^q} \right)$ 表现为各条路径对应求和的形式，所以发射站和接收站的选择向量通过影响量测费歇尔信息矩阵来影响最终的跟踪精度

${{{J}}_{\rm{Z}}}\!\left( {{{x}}_k^q,{{t}},{{r}}} \right) \!\!=\!\! \sum\limits_{m = 1}^M\! {\sum\limits_{n = 1}^N {t_m^qr_n^q \!{{\left(\! {{{H}}_{m,n,k}^q} \right)}^{\rm{T}}}\!\!\!{{{J}}_{\rm{Z}}}\!\left( {{{θ}} _{m,n,k}^q} \right)\!\!{{H}}_{m,n,k}^q} } $

(22)

费歇尔信息矩阵逆矩阵的迹是一个常用的用来进行优化选择的准则 ^{[6, 7]} 。由于雷达网络资源的限制，针对每一个目标只能有有限数量的发射站和接收站对其进行观察。本文通过选择发射站和接收站来获得更低的PCRLB，通常情况下也能获得更低的MSE。同时由于多目标的存在，需要平衡每个目标之间的跟踪性能。在本文的框架中以优化各个目标中跟踪性能最差的目标达到最优来实现系统整体性能的优化。于是本文给出以下优化模型

$\left. {\begin{array}{*{20}{l}}{\mathop {\min }\limits_{{{t}},{{r}}} \;\mathop {\max }\limits_q \left[ {{\mathop{\rm tr}\nolimits} {{\left( {{{{J}}_k}} \right)}^{ - 1}}} \right]}\\\begin{array}{l}{\rm{subject\; to}}\;K_{\min }^q \le {\left( {{{u}}_1^q} \right)^{\rm{T}}}{{t}} \le K_{\max }^q,q \in {\cal Q}\\\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;L_{\min }^q \le {\left( {{{u}}_2^q} \right)^{\rm{T}}}{{r}} \le L_{\max }^q,q \in {\cal Q}\\\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\left( {{{u}}_3^m} \right)^{\rm{T}}}{{t}} \le 1, \ m \in {\cal M}\\\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\left( {{{u}}_4^n} \right)^{\rm{T}}}{{r}} \le 1, \ n \in {\cal N}\\\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;t_m^q \in \left\{ {0,1} \right\}, \ m \in {\cal M},q \in {\cal Q}\\\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;r_n^q \in \left\{ {0,1} \right\}, \ n \in {\cal N},q \in {\cal Q}\end{array}\end{array}} \!\!\!\!\!\!\!\right\}$

(23)

其中， ${{u}}_1^q$ 和 ${{u}}_2^q$ 为两组向量分别除了第 $\left( {q - 1} \right)M + 1$ 到 qM和 $\left( {q - 1} \right)N + 1$ 到 qN为1，其余元素都为0的向量。 ${{u}}_3^m$ 和 ${{u}}_4^n$ 为两组向量其中分别第 $\left( {q - 1} \right)M + m$ 和第 $\left( {q - 1} \right)N + n$ 个元素为1其余元素为0。该优化问题中，第1个和第2个约束限制了可被选择的发射站和接收站数量的最大最小数目。最小数量约束用来保证最大似然估计的有效性，而最大数量约束保证所选发射接收站数量不超过系统所提供的数量限制。值得注意的是，最优选择数量的确定不在本文的讨论范围之内。第3个和第4个约束是基于假设每一个发射站和接收站在某一时刻只能照射和观察一个目标。

因为最后两个非凸约束的存在导致该优化问题不可解。所以为了有效的求解该优化问题，本文把最后两个非凸约束松弛为其对应的凸包 $t_m^q \in \left[ {0,1} \right],m \in {\cal M},q \in {\cal Q}$ 和 $r_n^q \!\in\! \left[ {0,1} \right],n \!\in\! {\cal N},q \!\in\! {\cal Q}$ 。可以证明对于 $\forall {{X}} \in {{\mathbb{S}}^ + }$ , ${\mathop{\rm tr}\nolimits} \left( {{{{X}}^{ - 1}}} \right)$ 是 X 的凸函数 ^[13] 。又因为BFIM是选择向量 t 和 r 的仿射函数，所以该问题关于选择向量 t 和 r 分别为一个凸优化问题。引入一组辅助向量 ${{v}} = {\left[ {{v_1}, \cdots ,{v_L}} \right]^{\rm{T}}}$ 满足

${v_l}\! \ge\! {{e}}_l^{\rm{T}}\!{\left[ {{{J}}_{{\rm P},k}^q \!+\! {{J}}_{\rm Z}^q\left( {{{x}}_k^q,{{t}},{{r}}} \right)} \right]^{ - 1}}\!{{{e}}_l}, \ l \!\in\! {\cal L},q \!\in\! {\cal Q}$

(24)

利用舒尔补定理，可以把优化问题转化为一个半正定规划(SemiDefinite Programming, SDP)问题 1 0

$\left. {\begin{array}{*{20}{l}}{\mathop {\min }\limits_{{{t}},{{r}},{{v}}} \;{{{1}}^{\rm{T}}}{{v}}}\\\begin{array}{l}{\rm{subject\; to}}\!\left[\!\! {\begin{array}{*{20}{c}}{{{{J}}_{{\rm{P}},k}} + {{{J}}_{\rm{Z}}}\left( {{{x}}_k^q,{{t}},{{r}}} \right)} & {{{{e}}_l}}\\{{{e}}_l^{\rm{T}}} & {{v_l}}\end{array}} \!\!\right]\\ \quad \quad \quad \quad \ \!\succeq\! {{{0}}_{L + 1}}, \ l \in {\cal L},q \in {\cal Q}\\\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\left( {{{u}}_1^q} \right)^{\rm{T}}}{{t}} = {K_{\max }}, \ q \in {\cal Q}\\\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\left( {{{u}}_2^q} \right)^{\rm{T}}}{{r}} = {L_{\max }}, \ q \in {\cal Q}\\\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\left( {{{u}}_3^m} \right)^{\rm{T}}}{{t}} \le 1, \ m \in {\cal M}\\\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\left( {{{u}}_4^n} \right)^{\rm{T}}}{{r}} \le 1, \ n \in {\cal N}\\\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;t_m^q \in \left[ {0,1} \right],m \in {\cal M}, \ q \in {\cal Q}\\\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;r_n^q \in \left[ {0,1} \right],n \in {\cal N}, \ q \in {\cal Q}\end{array}\end{array}} \right\}$

(25)

其中， e _l 是 L维单位矩阵的第 l列向量。由于该转化后的优化问题关于选择 t 和 r 分别为一个凸优化问题， 所以本文直觉的选择坐标下降法(Bolck Coordinate Descent, BCD)来进行优化。坐标下降法可以描述为在每一次迭代中，本文最小化 t 或者 r 中的一个，同时保持另一个选择向量固定。通过此种方法，发射站和接收站联合选择可以总结步骤为：

1. 利用目标状态预测估计状态 ${\hat {{x}}}_{k|k - 1}^q$ 计算各条路径的先验费歇尔信息矩阵 ${{{J}}_{{\rm P},k}}$ 以及量测费歇尔信息矩阵 ${{{J}}_{\rm{Z}}}\left[ {{{{θ}} _{m,n}}\left( {{\hat {{x}}}_{k|k - 1}^q} \right)} \right],m \in {\cal M},n \in {\cal N}$ 。

2. 给定一个迭代初始可行点 ${{\hat {{r}}}_0} \in {{\cal F}_2}$ 和收敛容限 $\delta $ 。

3. 针对式(25)中的问题，对 $s = 1,2,3 \cdots ,S$ 进行 S次迭代

　(1) $\left( {{{{\hat {{t}}}}_s},{{v}}} \right) = \arg \mathop {\min }\limits_{{{t}} \in {{\cal F}_1},\Omega } \;f\left( {{{t}},{{{\hat {{r}}}}_{s - 1}},{{v}}} \right)$ ；

　(2) $\left( {{{{\hat {{r}}}}_s},{{v}}} \right) = \arg \mathop {\min }\limits_{{{r}} \in {{\cal F}_2},\Omega } \;f\left( {{{{\hat {{t}}}}_s},{{r}},{{v}}} \right)$ ；

　(3) 如果 $\left| {f({{{t}}_s},{{{r}}_s},{{v}}) - f\left( {{{{t}}_{s - 1}},{{{r}}_{s - 1}},{{v}}} \right)} \right| < \delta $ 或达到 S停止，得到最优解 t ^*和 r ^*。

4. 对 $q = 1,2, \cdots ,Q$

　(a) 令 ${\left( {{{{t}}^q}} \right)^*} = {{{t}}^*}\left( {\left( {q - 1} \right)M:qM} \right)$ 和 ${\left( {{{{r}}^q}} \right)^*} = $ $ {{{r}}^*}\left( {\left( {q - 1} \right)N:qN} \right)$

　(b) 对 $m \in {\cal M}$ ，定义示性函数 ${\left( {t_m^q} \right)^ * } =$ $\left\{ {\begin{array}{*{20}{c}} 1, & {{t_{mq}} \ge \min \left\{ {{{\left[ {t_m^q} \right]}_{K_{\min }^q}},1/2} \right\}}\\ 0, & {{\rm{otherwise}}} \quad \quad \quad \quad \quad \quad \ \,\end{array}} \right.$

　(c) 对 $n \in {\cal N}$ ，定义示性函数 ${\left( {r_n^q} \right)^*} = $ $ \left\{ {\begin{array}{*{20}{c}} 1, & {{r_{nq}} \ge \min \left\{ {{{\left[ {r_n^q} \right]}_{L_{\min }^q}},1/2} \right\}}\\ 0, & {{\rm{otherwise}}} \quad \quad \quad \quad \quad \quad \ \, \end{array}} \right.$

5. 得到发射站和接收站的选择集合 ${{\cal T}_q}$ 和 ${{\cal R}_q}$ 对 $q \in {\cal Q}$ 。

其中， ${{\cal F}_1}$ 和 ${{\cal F}_2}$ 分别是 t 和 r 的可行集。 ${\left[ {t_m^q} \right]_{K_{\min }^q}}$ 和 ${\left[ {r_n^q} \right]_{L_{\min }^q}}$ 代表 ${{{t}}^q}$ 和 ${{{r}}^q}$ 中最大的第 $K_{\min }^q$ 个和最大的第 $L_{\min }^q$ 个元素。 f是优化问题2的目标函数。采样重要性重采样(Sample Importance Resample)粒子滤波 ^[14] 算法被应用于本文的目标跟踪并且假设数据关联算法已经被有效采用。在每次迭代中，可以通过内点法 ^[15] 进行有效的求解半正定规划问题(SDP)。文献[16]指出，只要在每一次迭代中，目标函数是凸函数或者可行解唯一，分块坐标下降的收敛性就可以保证。尽管在理论上无法保证全局最优性，但是仿真实验说明本文提出的方法是确实有效的。

6 仿真实验

在本节中，本文进行了计算机仿真实验来验证所提出的优化模型的有效性。仿真实验结果表明，本文所提出的动态规划方法使取得了比固定选择更优的性能，并且能在计算量更小的前提下获得和全局穷举搜索近似的性能。

在第1个仿真场景中，为了呈现发射站和接收站的动态选择过程，首先考虑一个规模相对较大的分布式MIMO雷达网络。12个发射站和12个接收站分别随机布置在 $\left[ {100\;{\rm{km}} \times 100\;{\rm{km}}} \right]$ 大小的区域中。总共2个目标的初始状态分别被设定为 $\left[ {\begin{array}{*{20}{c}} {10\;\;{\rm{km}}} \!&\! {10\;\;{\rm{km}}} \!&\! {150\;\;{\rm{m/s}}} \!&\! {150\;\;{\rm{m/s}}} \end{array}} \right]$ 和 $ \Big[ 90\;{\rm{km}} {90\;{\rm{km}} \ - 150\;{\rm{m}}/{\rm{s}} \ - 150{\rm{m}}/{\rm{s}}} \Big] $ 。假设跟踪每次采样间隔为1 s，即 $\Delta T = 1\;{\rm{s}}$ 。允许的被选择的发射站和接收站最小和最大的数目分别为3和4。仿真中随机的选择4个发射站和4个接收站作为固定选择情况下的对比实验组。

在 图1中，本文给出了在时刻50 s, 200 s, 300 s下针对目标1和目标2的发射站接收站动态实时选择情况，其中“ $\Delta $ ”符号和“ $\nabla $ ”符号分别代表了发射站和接收站的位置，红色“ $ \Box $ ”框和“ $\diamondsuit $ ”框分别代表了被目标1和目标2选中的发射站和接收站。可以看到，通过所提出的方法，离目标最近的发射站和接收站总是被优先选择来提高接收信号的信噪比。另一方面，被选中发射站接收站尽量在空间上尽量以多个角度观察目标来取得更好的空间多样性。当目标1和目标2靠近后，选择方案对两者的性能需求进行了折中。值得说明的是，在式中假设每一个发射站和接收站在某一时刻只能照射和观察一个目标是为了仿真结果的有效展示而不至于混淆，各个雷达具备同时观察更多个目标的能力并不影响优化模型的建立和求解，该方法依然是有效的。

本文定义每个时刻所跟踪目标中性能最差的RMSE为：

${\rm{RMS}}{{\rm{E}}_k} = \mathop {\max }\limits_q \!\!{\left[\! {\sqrt {\frac{1}{{{\rm{MC}}}}\sum\limits_i^{{\rm{MC}}} {\left\| {{{x}}_k^q\left( i \right) \!\!-\!\!{\hat {{x}}}_{k|k}^q\left( i \right)} \right\|_2^2} } } \right]^{1/2}} \ \ \ \ $

(26)

其中，MC是蒙特卡洛仿真实验的次数。 ${{{x}}_k}$ 是目标状态向量而 ${{\hat {{x}}}_{k|k}}$ 是目标状态估计向量。在 图2中，为了对各种方法进行对比，对某次实验的选择方案进行重复的1000次蒙特卡洛实验进行性能对比。本文只给出了每个时刻两种方法下跟踪性能最差的目标的RMSE曲线来简化图形。可以看到，本文所提出的动态优化选择方法相比于固定的选择方法，在不同时刻都能以更优的性能完成对最差性能目标的跟踪，实现了系统的整体性能优化。也就是说，本文提出的方法可以利用每一时刻的先验信息动态来调整针对各个目标的发射站接收站选择结果，从而获得更小的跟踪误差。

当发射站和接收站数量较大的时候穷举搜索法因为计算量太大而变得不可行。为了验证所提方法的全局最优性，本文考虑将发射站和接收站的总数分别减少为8个和8个，各个目标被允许的发射站和接收站数目保持4个不变。同 图2一样，在给出一次完成选择方案后进行1000次蒙特卡洛仿真实验求得两种方法的RMSE。可以看到， 图3中，本文所提方法与穷举搜索有非常相近的跟踪性能。

为了验证算法的稳定性，在同一场景下100次重复试验下文中方法的平均PCRLB对比和由穷举搜索获得的PCRLB。 图4表明本文基于分块坐标下降方法的优化选择方法在保证稳健性的基础上取得了近乎穷举搜索的性能。

给出本文算法准确的计算量比较困难，相比于穷举搜索指数级增长的计算量 ${\cal O} \sim \left( {{2^{QM + QN}}} \right)$ ，本文基于半正定规划的计算量约为 ${\cal O}\! \sim \!\left[ \! {S{{\left( {QM} \right)}^3} \!\!+\!\! S{{\left( {QN} \right)}^3}} \!\right]$ ^[15] 。为了直观对比两者计算量的差距， 图5给出了在同样的计算机下不同场景下两者的执行时间对比图。其中Tx表示发射站，Rx表示接收站。可见在8发8收的场景下进行选择，穷举搜索法的计算时间略小于本文的方法。但是当处于9发9收的情况，穷举搜索法计算量出现组合爆炸现象，基本已经失去了实用意义。

7 结论

本文研究了在分布式MIMO雷达网络中，多目标跟踪背景下的发射站接收站联合选择问题。一种基于以后验卡拉美罗界为性能指标，旨在优化多目标跟踪下跟踪性能最差的目标的发射站接收站优化问题被提出。通过分块坐标下降该优化问题得以求解，实现了发射站和接收站的动态闭环选择。通过仿真实验，在每一个时刻，本文提出的方法可以利用前面时刻的先验信息形成预先规划针对每一个目标的发射站接收站选择。相比于非动态规划下的固定选择，本文的优化选择方法取得了更好的性能。于此同时，对比全局最优的穷举搜索在更小的计算量的前提下获得了近似的性能。