一种机载雷达网络中目标跟踪的机动资源分配方法

1.本发明涉及雷达资源分配领域，具体涉及一种机载雷达网络中目标跟踪的机动资源分配方法。


背景技术：

2.随着集成电路技术的发展，在机载平台上装备了越来越多的雷达，因而也建立了一种新型的雷达资源——机动资源。而在雷达工作时，雷达资源预算与跟踪性能之间存在矛盾，所以资源分配对于雷达网络系统中的目标跟踪至关重要，必须协同分配多个节点的资源，以确保雷达网络能够在有限的资源预算下完成跟踪任务。
3.为了实现这一点，提出了许多面向跟踪任务的雷达资源分配方法，其中将资源分配方案设计为优化问题，优化每个节点的机动资源，设计其轨迹，以获得更好的雷达空间观测多样性和目标跟踪性能，来最大限度地提高雷达网络的资源利用效率。但是大多数现有方法没有考虑雷达的机动性，并且现有的航迹规划工作也没有考虑到雷达探测背景。


技术实现要素：

4.针对现有技术中存在的问题，本发明的目的在于提供一种机载雷达网络中目标跟踪的机动资源分配方法，来实现对机动资源的充分利用，以获得更好的雷达目标跟踪性能。
5.为了达到上述目的，本发明采用以下技术方案予以实现。
6.一种机载雷达网络中目标跟踪的机动资源分配方法，包括以下步骤：
7.步骤1，建立一组具有i＝1，2，...，n个节点的机载雷达网络，n为大于0的整数；定义一个分配间隔t0；令k表示k时刻，k的初始值为1，k∈{1，2，
…
，k}，k 为预先设定的最大跟踪时刻，k为大于0的正整数；
8.步骤2，n个雷达节点将接收到的量测进行融合处理，融合中心根据所述回波数据构建目标量测模型，确定第k时刻目标状态和节点状态；
9.步骤3，根据第k时刻目标状态和节点状态，确定运动模型；
10.步骤4，根据贝叶斯克拉美罗界确定目标函数；根据目标存在环境的实际情况，确定目标的约束条件；
11.步骤5，消除非积极约束；
12.步骤6，对积极约束线性化；
13.步骤7，利用近端交替方向乘子法，求解优化问题，得到机载雷达网络中各个节点雷达在x方向和y方向上的速度；通过雷达在x方向和y方向上的速度，再得出各个节点雷达的平飞速度和转角速度
14.步骤8，当k未达到预先设定的最大跟踪时刻k时，令k加1，返回步骤2；当k达到预先设定的最大跟踪时刻k时，迭代重复执行过程停止，得到第1时刻各个节点雷达的平飞速度和转角速度至第k时刻各个节点雷达的平飞速度和转角速度并记为一种机载
雷达网络中目标跟踪的机动资源分配的结果。
15.与现有技术相比，本发明的有益效果为：在保证收敛的条件下解决了一个具有多重线性不等式/等式约束的非凸目标函数的问题，同时避免了应用梯度投影方法在多面体上进行投影时计算代价的高昂；能实现对机动资源的充分利用，能在更短的计算时间内获得良好的目标跟踪精度。
附图说明
16.下面结合附图和具体实施例对本发明做进一步详细说明。
17.图1为本发明机载雷达网络中目标跟踪的机动资源分配方法的流程示意图
18.图2为本发明的仿真实验设定的机载雷达网络中各个雷达、目标位置以及威胁区域的空间关系图
19.图3为本发明在每个分配时间单元k上的数学表达式中第三个约束处于非积极状态说明图；
20.图4为本发明在每个分配时间单元k上的数学表达式中第四个约束处于非积极状态说明图；
21.图5为扩展点的说明图；
22.图6为约束条件下扩展点的说明图；
23.图7为平飞速度和转角速度的限制的近似范围说明图；
24.图8为在本发明方法、对照方案1方法和对照方案2方法下，机动资源分配的结果对比图；
25.图9为在本发明方法、对照方案1方法和对照方案2方法下，两个相近雷达之间的距离对比图；
26.图10为在本发明方法、对照方案1方法和对照方案2方法下，目标跟踪的贝叶斯克拉美罗下界对比图；
27.图11为本发明方法与对照方案2方法下计算时间的对比图。
具体实施方式
28.下面将结合实施例对本发明的实施方案进行详细描述，但是本领域的技术人员将会理解，下列实施例仅用于说明本发明，而不应视为限制本发明的范围。
29.参考图1，为本发明机载雷达网络中目标跟踪的机动资源分配方法的流程示意图；该机载雷达网络中目标跟踪的机动资源分配方法，包括以下步骤：
30.步骤1，建立一组具有i＝1，2，...，n个节点的机载雷达网络，n为大于0的整数；定义一个分配间隔t0；令k表示k时刻，k的初始值为1，k∈{1，2，
…
，k}，k 为预先设定的最大跟踪时刻，k为大于0的正整数；
31.步骤2，n个雷达节点将接收到的量测进行融合处理，融合中心根据所述回波数据构建目标量测模型，确定第k时刻目标状态和节点状态。
32.在第k个分配区间内，雷达i可以为目标产生一系列的量测m
i，k
≥1。定义为雷达i在k
i，m
(表示时间单元)的第m个量测值。
33.所述目标量测模型，其表达式为式(1)：
[0034][0035]
式(1)中，表示k时刻第i个雷达的第m个量测所对应的状态向量，是一个复合状态向量，包含目标的状态向量和机载雷达的状态向量；其中目标k时刻的状态向量表示为其中下标t代表目标，x
t，k
为目标的水平位置，为目标的水平速度，y
t，k
为目标的垂直位置，为目标的垂直速度；机载雷达的状态向量表示为其中表示雷达i的空间方位角，r代表雷达，为雷达i的水平位置，为雷达i的垂直位置；将n个雷达状态向量链接到单个矢量中，得到
[0036]
式(1)中，链接了不同的量测值，其中表示第i个雷达在第k个时刻得到的第m个量测中目标的距离，表示第i个雷达在第k时刻得到的第m个量测中目标的方位角，和表示雷达i在这次量测中的位置；
[0037]
式(1)中，是雷达i的量测噪声，设为零均值的高斯噪声，其协方差为式(2)：
[0038][0039]
式(2)中，是雷达的平飞速度，是雷达转角速度；和是gps测得信息的对应的误差方差，与目标信噪比无关；和是距离和方位测量的误差方差，且与目标信噪比呈反线性，如式(3)所示：
[0040][0041]
式(3)中，信噪比是第i个雷达节点在第k个时间单元得到第m 个量测时的信噪比，它与每个节点直接相关的轨迹控制向量d
i，k
有关；d
i，k
是需要分配的机动雷达资源，它是机载雷达网络的轨迹控制向量；β
i2
和为雷达i 的传输信号带宽和3db接收波束宽度。
[0042]
步骤3，根据第k时刻目标状态和节点状态，确定运动模型；
[0043]
具体的，将目标的状态向量链接到中，可得到一个扩展向量并且这个状态向量的维数是3n 4；对于ξk，其状态转移方程如式(4)：
[0044]
ξk＝f(ξ
k-1
) ukꢀꢀꢀ
(4)
[0045]
式(4)中，其中f
t
(
·
)为已知的状态转移函数；uk是零均值高斯过程噪声，其协方差为qk，；是机载雷达网络的轨迹控制向量；fr(
·
)表示机载雷达运动函数，表示为式(5)：
[0046][0047]
式(5)中，
[0048]
步骤4，根据贝叶斯克拉美罗界确定目标函数；根据目标存在环境的实际情况，确定目标的约束条件；
[0049]
直观地说，向目标移动肯定会增加目标的信噪比，减少其测量的误差方差。然而，由于一些物理约束，多个节点可能不会直接向目标移动。从数学上讲，机动资源分配方案可以表述为一个最大化目标跟踪精度的问题，并受到一些机动资源的约束，所考虑的参数是多个机载雷达的机动资源变量
[0050]
具体的，通过使用复合量测，可以将克拉美罗下界矩阵近似为式(6)：
[0051][0052]
在式(6)中，f是状态转移函数的雅可比矩阵；是在ξ
k 1
下关于第m次量测估计的雅可比矩阵；j是fisher信息矩阵；
[0053]
为了使目标精度最大化，机动资源分配方案在每个分配时间单元k上的数学表示可以表示为式(7)：
[0054][0055]
其目标函数如式(8)所示，用来描述目标的跟踪精度。
[0056][0057]
式(8)中，λ是归一化矩阵，如式(9)：
[0058][0059]
式(9)中，表示克罗内克算子，i2表示2阶的单位矩阵。
[0060]
式(7)约束意味着平飞速度和转角速度受到最小值和最大值的限制，每一个机载雷达都不能通过多个威胁区，因为它需要避免被敌人的雷达探测到并且任何两个机载雷达都需要在每个分配时间单元上间隔一个最小距离d
min
。
[0061]
在模型中，对于威胁区中每个敌方雷达都被认为是静止的，用和来表示敌方第l个雷达的探测半径和中心。
[0062]
步骤5，消除非积极约束。
[0063]
鉴于机动资源分配方案在每个分配时间单元k上的数学表达式通常是一个困难的非线性优化问题，因为它有一个非凸的目标函数和太多的非凸约束。尽管目标函数具有非凸性，但有效的求解仍然面临两种困难。一方面，非线性和非凸约束带来了通过优化算法去生成可行性迭代的困难。另一方面，约束的数量很大，这将大大增加计算负担。为了解决这两个困难，设计了一种包含非积极约束消除和积极约束线性化的三步求解技术，以解决原始的非凸优化问题。
[0064]
在实践当中，因为雷达与威胁区域以及雷达与雷达之间的距离可能比较远，因此，许多约束可能是非积极的。消除非积极约束将大大减少计算量，因此，要进行非积极约束的消除。
[0065]
在以下命题中，给出了对于判断末位两个约束是否积极的充分条件。
[0066]
命题1：设r为雷达在分配间隔内可以移动的最大距离，为雷达i在第k 个分配时间单元上的已知位置。然后末位两个约束在以下条件下是不积极的：
[0067]
a)对于如果那么第三个约束条件是不积极的。
[0068]
b)对于如果那么第四个约束条件是不积极的。
[0069]
命题1的证明是显然的，因为只需要证明那些被消除的非积极约束对于给定的前两个约束总是成立的。对于命题1最直观的解释如图3所示。
[0070]
参考图3，黑色十字叉表示雷达i的位置，红色圆圈区域表示威胁区域，蓝色扇区区域表示雷达在下一个分配时间内所能达到的可能位置。如果雷达远离威胁区域，只要前两个约束条件保持不变，则第三个约束始终为非积极约束。参考图4，也可以得出类似的结论。第四个约束是避撞约束。如果雷达彼此远离，它们将保持安全距离以避撞。
[0071]
因此，机动资源分配方案在每个分配时间单元k上的数学表达式可以改写为式(10)。
[0072][0073]
式(10)中，φ
i，k
和μ
i，k
为两个集合，定义如下：
[0074]
对于节点i的威胁区域约束，如果满足则我们称该约束是积极的，且第i个雷达在第k个时间段内的全体威胁区域l的集合为φ
i，k
；对于节点i的碰撞约束，如果满足则我们称该约束是积极的，且第i个雷达在第k个时间段内的全体可能碰撞的节点j的集合为μ
i，k
。
[0075]
步骤6，对积极约束线性化。
[0076]
具体的，引入一组辅助变量：
[0077][0078][0079]
其中，和分别表示雷达在x方向和y方向上的速度分量。
[0080]
基于辅助向量，约束条件可以近似为线性约束。
[0081]
机动资源分配方案在每个分配时间单元k上的数学表达式共有四个约束条件。
[0082]
子步骤6.1，对第三个约束线性化。
[0083]
根据不等式x2 y2≥2xy，可以将第三个非线性约束放松为线性约束。具体来说，对于任何的有因此，第三个非线性约束的左侧允许一个线性形式的下界，如式(11)：
[0084][0085]
式(11)中，和是扩展点。
[0086]
根据式(11)，第三个约束可以作为一个线性约束来放松，如式(12)：
[0087][0088]
第三个约束条件的边界线可以指定为式(13)：
[0089][0090]
在数学上，对线性化过程有以下性质：
[0091]
命题2：对于任何的都要满足方程上式的边界线与第1个威胁区的边界相切，是一个相切点。
[0092]
证明：与边界线之间的距离可以表示为：
[0093][0094]
因此，边界线与威胁区相切。将带入边界线，可以验证它是边界线上的一个点。另外，位于威胁区的边界上，因此看到是一个切点。这就完成了证明。
[0095]
参考图5，可以看出，扩展点最好选择在弧的中间，该弧位于威胁区与可到达区域相交的边界(图5中的蓝色区域)；经过线性化处理，红色圆圈外的区域被近似为了黑线外的区域。
[0096]
将辅助变量g
i，k
和机载雷达的状态转移函数表达式代入已经放松后的第三个约束条件，得到式(13)：
[0097]
χ
i，l，k
·gi，k
≥ζ
i，l，k
ꢀꢀꢀ
(13)
[0098]
式(13)中，χ
i，l，k
和ζ
i，l，k
与预定义的参数相关，如扩展点雷达的位置敌方雷达的中心和探测半径等。
[0099]
子步骤6.2，对第四个约束进行线性化。
[0100]
为处理第四个非线性约束，也转向不等式x2 y2≥2xy。在此基础上，第四个约束可以放宽为式(14)：
[0101][0102]
从式(14)可以看出，第四个约束与呈线性关系，位于四维欧几里得空间中；当给出或时，式(14)的左侧在二维欧几里得空间中指定了两条平行线。
[0103]
因此，第四个约束可以看作是二维欧几里得空间中的两个平行线性约束，其中可行区域对应于两条平行线的外部，如图6所示。
[0104]
这两条平行的线可以写成式(15)：
[0105][0106]
通过正确选择和式(15)指定的线有以下属性。
[0107]
命题3：对于任意和满足下面两个条件的线，则上式中两条平行线之间的距离等于d
min
，两条线都垂直于由点和确定的线。
[0108][0109][0110]
证明：由和确定的线的斜率为同时，由于和满足第一个条件，两条平行线的斜率可以推导为
[0111]
因此，有：
[0112][0113]
这意味着这两条平行线垂直于由点和确定的线。此外，与命题 2的证明相似，可以计算出当和满足第二个条件时两条平行线之间距离为d
min
。这就完成了证明。
[0114]
命题3表明，如果扩展点和满足条件，则放宽后第四个约束(线性约束)可以保证和之间的距离大于或等于d
min
。
[0115]
将辅助变量g
i，k
和机载雷达的状态转移函数表达式代入已经放松后的第四个约束条件，有式(16)：
[0116][0117]
式(16)中，γ
i，j，k
和ε
i，j，k
与多个预定义参数相关，如扩展点和两个雷达和的位置、以及分配时间间隔t0等。
[0118]
子步骤6.3，对第一个和第二个约束进行线性化。
[0119]
用g
i，k
将第一个与第二个约束进行重新表达，且对重新表达的式子进行线性化。上述线性化过程可以有效地将非线性约束转换为线性约束。
[0120]
用g
i，k
代替和则有：
[0121][0122][0123]
线性平飞速度和转角速度约束被表示为非线性约束，这很难处理。
[0124]
然后，将前两个约束表示的非凸的扇形区域近似为一个凸的梯形区域，如图7所示。为了得到梯形区域的四个边，需要得到顶点a、b和o的坐标。
[0125][0126][0127][0128]
同时，图7中的线l1和l2可以写成：
[0129][0130][0131]
根据oa、ob、l1和l2，可以得到图7中的梯形面积。
[0132]
这个梯形面积用四个线性约束表示，可以缩写为式(19)：
[0133][0134]
式(19)中，a
i，k
和c
i，k
与和相关。
[0135]
子步骤6.4，结合上述线性化过程，机动资源分配方案在每个分配时间单元k上的数学表达式可以写为式(20)。
[0136][0137]
步骤7，利用近端交替方向乘子法，求解优化问题，得到机载雷达网络中各个节点雷达在x方向和y方向上的速度；通过雷达在x方向和y方向上的速度，再得出各个节点雷达的平飞速度和转角速度
[0138]
式(20)是一个具有多重线性不等式/等式约束的非凸目标函数。为了保证收敛性且解决这类问题，通常采用基于梯度投影的方法。然而，在多面体上的投影(对应于式(20)中的线性不等式/等式约束)会造成很大的计算量。
[0139]
具体的，在式(20)中，两个变量g
i，k
、g
j，k
是线性耦合的。
[0140]
引入解耦变量来解耦，要优化的问题可以写为式 (21)：
[0141][0142]
使用近端交替方向乘子法，将辅助变量ek，sk引入不等式约束，并将式(21) 重写为式(22)：
[0143][0144]
近端交替方向乘子法的一个显著特点是引入一个平滑迭代(即指数加权) 序列，并在每次迭代中包含以当前平滑原始迭代为中心的二次近端项的增广拉格朗日函数。
[0145]
令：
[0146][0147]
式(23)是增广拉格朗日函数，其中χk＝[...，χ
i，l，k
，...]和ζk＝[...，ζ
i，l，k
，...]表
示的是一组χ
i，l，k
和ζ
i，l，k
，b
kgk
＝zk表示的是一组约束y＝[ye，ys，yw]表示对偶变量；ρe，ρs，ρw＞0是原始的惩罚参数；p是一个预定义的正参数，可以确保l(gk，zk，ek，sk，bk；y)在gk处是强凸的；bk是引入的平滑迭代(即指数加权)序列，将该序列插入增广的拉格朗日函数中，就可以保证下一次的迭代不会太过偏离，保证不精确的增广拉格朗日乘子方法在非凸问题上的稳定。
[0148]
近端乘子交替方向算法可以总结如下：
[0149]
第一步，令ρe，ρs，ρw＞0，λe，λs，λw＞0，0＜β≤1，c1，c2＞0；
[0150]
第二步，初始化y
e，0
，y
s，0
，y
w，0
，其中和均大于等于0；
[0151]
第三步，确定φ
i，k
和μ
i，k
；
[0152]
第四步，计算χ
i，l，k
和γ
i，l，k
和
[0153]
第五步，对于t＝1，2，3，
…
，开始做循环，循环以下公式：
[0154][0155][0156][0157][0158]
且通过梯度投影法进行更新；
[0159]
其中，和分别通过以下三个公式计算得出：
[0160][0161][0162][0163]
其中，[
·
]

是操作符，将一个向量投影到非负象限，确保向量的所有元素都为正。
[0164]
第六步，结束循环；
[0165]
得到机载雷达网络中各个节点雷达在x方向和y方向上的速度；
[0166]
通过雷达在x方向和y方向上的速度，根据计算得出各个节点雷达的平飞速度和转角速度
[0167]
步骤8，当k未达到预先设定的最大跟踪时刻k时，令k加1，返回步骤2；当k达到预先设定的最大跟踪时刻k时，迭代重复执行过程停止，得到第1时刻各个节点雷达的平飞速度和转角速度至第k时刻各个节点雷达的平飞速度和转角并记为一种机载雷达网络中目标跟踪的机动资源分配的结果。
[0168]
下面通过仿真实验对本发明效果作进一步验证说明。
[0169]
(一)仿真条件：
[0170]
仿真条件是在一台2.9ghz英特尔(r)酷睿(tm)i7-4790cpu和安装了 matlab9.8.0(2020a)的计算机上进行的。
[0171]
(二)仿真内容与结果分析：
[0172]
参考图2，为本发明的仿真实验设定的机载雷达网络中各个雷达、目标位置以及威胁区域的空间关系图。
[0173]
机载雷达数为n＝5，最大转角速度w
max
＝0.03rad/s，最小平飞速度v
min
和最大平飞速度v
max
设为(25，100)m/s。
[0174]
假设目标的雷达横截面是非波动的，因此目标的信噪比只与雷达-目标的距离有关。当雷达-目标距离为120km时，我们将一个参考目标的信噪比设置为12db。假设目标以恒定的50m/s的速度移动。将连续分配时间单元之间的间隔设置为t0＝20s，并利用20个分配间隔来支持我们的仿真。为了避撞，d
min
应该大于2v
max
t0，因此我们设置d
min
＝5km。
[0175]
雷达的初始状态、目标和威胁区域的位置如图2所示。
[0176]
为了初始化近端交替方向乘子法的参数，设置了原始步长c1＝0.05和 c2＝103。设置惩罚参数为ρe＝2，ρs＝1
×
10-4
，ρw＝3
×
10-3
来平衡迭代过程中梯度的不同项。最大迭代次数设置为120次。此外，初始化了对偶和近端步长的值λe＝4，λs＝5
×
10-5
，λw＝1.5
×
10-3
和β＝0.9(0＜β≤1)。
[0177]
为了说明本发明的优越性，我们将资源分配结果、相应的克拉美罗下界和计算时间与两个对照方案进行了比较。
[0178]
对照方案1，提前预定每个节点资源的实际情况。在对照方案1中，假设机载雷达以恒定的速度80m/s移动，而其初始位置和方向角度如图2所示。
[0179]
对照方案2，使用遗传算法来解决问题。遗传算法的初始条件如下：种群规模为200，交叉率为0.8，突变率为0.2，最大代数为100。
[0180]
图8为三种方法的机动资源分配结果。可以看出，看到雷达3可以在我们的方案和对照方案2中，避开威胁区域，继续向目标移动，这说明了本发明的正确性。在对照方案1中，雷达3以恒定的轨迹移动，因此它可能不会移动到目标，更糟糕的是，如果其方向角度设置不正确，它将会通过威胁区域。
[0181]
这三种方法对应的雷达1与雷达2之间的距离如图9所示。我们的方案和对照方案2都可以确保距离大于d
min
。相反，在对照方案1中，由于雷达已经提前预定了轨迹来进行移动，两个雷达之间的距离不能得以保证。
[0182]
从图10中可以看出，本发明和对照方案2相比，对照方案1有最差的跟踪精度，因为它采用了已经预定好的机动参数。
[0183]
为了分析计算负荷，将本发明的计算时间与对照方案2进行了比较。超过10次蒙特卡罗试验的平均运行时间，相应的结果如图11所示。参考图11 可以看出，与对照方案2相比，本发明可以在消耗更少的计算时间的同时实现等效的目标跟踪精度。
[0184]
综上所述，仿真实验验证了本发明的正确性，有效性和可靠性。
[0185]
虽然，本说明书中已经用一般性说明及具体实施方案对本发明作了详尽的描述，但在本发明基础上，可以对之作一些修改或改进，这对本领域技术人员而言是显而易见的。
因此，在不偏离本发明精神的基础上所做的这些修改或改进，均属于本发明要求保护的范围。