CN108919253A - 一种基于帝国竞争算法的天波雷达多路径数据关联方法 - Google Patents

Abstract

一种基于帝国竞争算法的天波雷达多路径数据关联方法，本发明涉及一种基于帝国竞争算法的天波雷达多路径数据关联方法。本发明的目的是为了解决现有MPDA算法存在计算量大、算法执行时间过长的缺点。过程为：一，建立天波雷达跟踪***模型；二，初始化，计算滤波预测值；三，随机生成初始国家并计算权力；四，初始帝国；五，殖民地同化作用；六，殖民地革命；七，帝国增强；八，殖民竞争；九，若迭代结束或只剩一个帝国，输出结果，否则，返回五；十，计算目标状态估计值，十一，判断目标跟踪是否结束，如果没有，使k＝k+1，执行二～十。本发明用于雷达数据处理领域。

Description

一种基于帝国竞争算法的天波雷达多路径数据关联方法

技术领域

本发明涉及一种基于帝国竞争算法的天波雷达多路径数据关联方法，属于雷达数据处理领域。

背景技术

数据关联模块是雷达数据处理***的关键组成部分，数据关联的作用是在杂波中将量测数据与真实目标航迹进行关联匹配，从而实现对目标的精准跟踪。数据关联算法性能的优劣直接影响到跟踪***整体的性能。

与常规雷达不同，天波雷达面临的环境更为复杂，且存在着低检测概率、低数据率、低测量精度、高虚警率和多径效应等问题，使其数据处理变得非常困难。其中，多径效应的存在降低了单一传播模式下的检测概率，提高了航迹失跟率，是天波雷达迫切需要解决的问题。由此可见，传统的数据关联算法已经不能满足天波雷达迅速建立航迹形成早期预警的要求，研究、改进和发展新的天波超视距雷达数据关联算法是十分迫切和必要的。

随着电离层探测技术的进步，学者开始探讨利用电离层信息解决多径传播问题，提出了多路径概率数据关联算法(Multipath Probability Data Association，MPDA)等算法。其中，MPDA算法将各种传播模式的回波作为备选的竞争量测，在最大似然意义下选择最优的回波，进行传播模式与目标的关联以及状态更新，是应用最为广泛的一种算法，但由于该算法需要穷举所有的关联事件，计算量会随着杂波密度的增加而指数增加，使算法成为NP难题，存在计算量大、算法执行时间过长的缺点。

发明内容

本发明的目的是为了解决现有MPDA算法存在计算量大、算法执行时间过长的缺点，不适用于实际工程，从而提出一种基于帝国竞争算法的天波雷达多路径数据关联方法。

一种基于帝国竞争算法的天波雷达多路径数据关联方法过程为：

步骤一，建立天波雷达运动目标的运动模型和多路径量测模型，将天波雷达运动目标地面坐标系下的运动状态标记为L维的运动状态向量x，天波雷达运动目标雷达坐标系下的量测值标记为M维的量测向量y，根据运动模型获得天波雷达状态误差协方差矩阵P和天波雷达过程噪声协方差矩阵Q，根据量测模型获得量测协方差矩阵R；

步骤二，根据步骤一获得的第k时刻天波雷达运动目标在地面坐标系下的运动状态向量x(k)、天波雷达状态误差协方差矩阵P(k)，计算天波雷达运动目标第k+1时刻的运动状态预测值和状态预测协方差矩阵P(k+1|k)，并由此得到天波雷达运动目标第k+1时刻传播模式i＝1,...,n_k+1在雷达坐标系下的量测预测值和量测预测协方差S_i(k+1)；n_k是时刻k传播模式的个数；

步骤三，传感器获取当前时刻天波雷达运动目标在雷达坐标系下的量测值，随机生成N个初始量测-模式关联对，即初始国家，并计算每个国家的权力c_t，t＝1，2，3，...，N，初始化当前迭代次数g＝1，预设算法迭代次数G；N取值为正整数；

步骤四，依据国家权力大小从高到低排序，选择前N_imp个作为殖民国家，其余N_col个作为殖民地，按殖民国家的权力大小将殖民地分配给初始殖民国家，形成初始帝国；

步骤五，殖民地同化作用：殖民国家以自身量测-模式关联对(步骤三得到的)为基础，依据一定同化率随机地替换和改变殖民地路径中的α个量测-模式关联对；

步骤六，殖民地革命：殖民地按革命概率P_m随机地选取殖民地两个位置的关联量测进行交换，然后计算殖民地革命前后的权力，保留权力大的殖民国家，并重新评估殖民国家和殖民地的权力，若后者大于前者，则取代前者成为新的帝国；

步骤七，帝国增强，根据每个殖民国家权力，对权力最小的国家进行增强，增强方式为：随机移除该殖民国家中β个关联量测，重新建立量测-模式关联对，重新计算殖民国家增强前后的权力，保留权力大的殖民国家；

步骤八，殖民竞争，以帝国为单位，计算帝国内殖民国家和所有殖民地的标准化权力TC_n；通过评估所有帝国权力，将最弱帝国的最弱殖民地分配给最强帝国；当某个帝国内没有任何殖民地时，帝国消失，该帝国的殖民国家也被分配给最强帝国作为殖民地；

步骤九，帝国竞争算法不断迭代运算，若达到预设算法的迭代次数G或只剩一个帝国，算法结束，输出最终的量测-模式关联对作为数据关联的结果，否则，令g＝g+1，返回步骤五；

步骤十，根据步骤九输出的关联结果计算当前时刻地面坐标系下的目标状态估计值；

步骤十一，判断目标跟踪是否结束，如果没有结束，使k＝k+1，执行步骤二～步骤十，如果结束则完成数据关联方法。

本发明的有益效果为：

本发明将天波雷达的多路径数据关联问题建模为组合优化问题，并将帝国竞争算法引入其中，为该领域提供了一个基于社会群智能算法的新思路，在保证算法跟踪精度的前提下大大缩短了算法执行时间和计算量；如图8到图11给出50次蒙特卡洛仿真后目标各维度目标各维度的均方根误差，从图中可以看出，随着跟踪的增加，目标径向距离、径向距离变化率、方位角、方位角变化率的均方根误差均逐渐减少并趋近于0，跟踪精度提高10倍左右，同样证明本发明可以有效解决天波雷达多路径数据关联问题，达到了收敛的效果。

附图说明

图1为基于帝国竞争算法的天波雷达多路径数据关联方法流程图；

图2为天波雷达运动目标在雷达坐标系下的径向距量测图；

图3为天波雷达运动目标在雷达坐标系下的多普勒量测图；

图4为天波雷达运动目标在雷达坐标系下的方位角量测图；

图5为天波雷达运动目标在地面坐标系下的目标跟踪航迹图；

图6为天波雷达运动目标在地面坐标系下的目标径向距离变化率跟踪图；

图7为天波雷达运动目标在地面坐标系下的目标方位角变化率跟踪图；

图8为天波雷达运动目标的径向距离均方根误差随时间变化曲线图；

图9为天波雷达运动目标的径向距离变化率均方根误差随时间变化曲线图；

图10为天波雷达运动目标的方位角均方根误差随时间变化曲线图；

图11为天波雷达运动目标的方位角变化率均方根误差随时间变化曲线图。

具体实施方式

具体实施方式一：结合图1说明本实施方式，本实施方式的一种基于帝国竞争算法的天波雷达多路径数据关联方法具体过程为：

步骤一，建立天波雷达运动目标的运动模型和多路径量测模型(多路径量测模型就是多径效应下天波雷达的***量测模型，包括量测向量和量测方程y(k))，将天波雷达运动目标地面坐标系下的运动状态标记为L维的运动状态向量x，天波雷达运动目标雷达坐标系下的量测值标记为M维的量测向量y，根据运动模型(包括运动状态向量和运动状态方程)获得天波雷达状态误差协方差矩阵P和天波雷达过程噪声协方差矩阵Q，根据量测模型获得量测协方差矩阵R；

步骤二，根据步骤一获得的第k时刻天波雷达运动目标在地面坐标系下的运动状态向量x(k)、天波雷达状态误差协方差矩阵P(k)，计算天波雷达运动目标第k+1时刻的运动状态预测值和状态预测协方差矩阵P(k+1|k)，并由此得到天波雷达运动目标第k+1时刻传播模式i＝1,...,n_k+1在雷达坐标系下的量测预测值和量测预测协方差S_i(k+1)；n_k是时刻k传播模式的个数；

步骤三，传感器获取当前时刻天波雷达运动目标在雷达坐标系下的量测值，随机生成N个初始量测-模式关联对，即初始国家，并计算每个国家的权力c_t，t＝1，2，3，...，N，初始化当前迭代次数g＝1，预设算法迭代次数G；N取值为正整数；

步骤四，依据国家权力大小从高到低排序，选择前N_imp个作为殖民国家，其余N_col个作为殖民地，按殖民国家的权力大小将殖民地分配给初始殖民国家，形成初始帝国；

步骤五，殖民地同化作用：殖民国家以自身量测-模式关联对(步骤三得到的)为基础，依据一定同化率随机地替换和改变殖民地路径中的α个量测-模式关联对；

步骤六，殖民地革命：殖民地按革命概率P_m随机地选取殖民地两个位置的关联量测进行交换，然后计算殖民地革命前后的权力，保留权力大的殖民国家，并重新评估殖民国家和殖民地的权力，若后者大于前者，则取代前者成为新的帝国；

步骤七，帝国增强，根据每个殖民国家权力，对权力最小的国家进行增强，增强方式为：随机移除该殖民国家中β个关联量测，重新建立量测-模式关联对，重新计算殖民国家增强前后的权力，保留权力大的殖民国家；

步骤八，殖民竞争，以帝国为单位，计算帝国内殖民国家和所有殖民地的标准化权力TC_n；通过评估所有帝国权力，将最弱帝国的最弱殖民地分配给最强帝国；当某个帝国内没有任何殖民地时，帝国消失，该帝国的殖民国家也被分配给最强帝国作为殖民地；

步骤九，帝国竞争算法不断迭代运算，若达到预设算法的迭代次数G或只剩一个帝国，算法结束，输出最终的量测-模式关联对作为数据关联的结果，否则，令g＝g+1，返回步骤五；

步骤十，根据步骤九输出的关联结果计算当前时刻地面坐标系下的目标状态估计值；

步骤十一，判断目标跟踪是否结束，如果没有结束，使k＝k+1，执行步骤二～步骤十，如果结束则完成数据关联方法。

具体实施方式二：本实施方式与具体实施方式一不同的是：所述步骤一中建立天波雷达运动目标的运动模型和多路径量测模型(多路径量测模型就是多径效应下天波雷达的***量测模型，包括量测向量和量测方程y(k))，将天波雷达运动目标地面坐标系下的运动状态标记为L维的运动状态向量x，天波雷达运动目标雷达坐标系下的量测值标记为M维的量测向量y，根据运动模型(包括运动状态向量和运动状态方程)获得天波雷达状态误差协方差矩阵P和天波雷达过程噪声协方差矩阵Q，根据量测模型获得量测协方差矩阵R；具体过程为：

在采样时刻k，一个匀速天波雷达运动目标在地面坐标系下的运动状态向量表示为：

其中，ρ(k)表示天波雷达运动目标在地面坐标系下的径向距离，表示天波雷达运动目标在地面坐标系下的径向距离变化率(径向速度)，b(k)表示天波雷达运动目标在地面坐标系下的方位角，表示天波雷达运动目标在地面坐标系下的方位角变化率，L＝4，分别为ρ(k)、b(k)、为求的转置；

则天波雷达运动目标运动状态方程表示为：

x(k+1)＝F(k)x(k)+v(k)

其中，F(k)是L×L阶的状态转移矩阵，L维过程噪声v(k)是一个零均值的高斯白噪声，v(k)协方差矩阵为Q(k)；

在采样时刻k，天波雷达运动目标一次量测值在雷达坐标系下的量测向量表示为：

y(k)＝[R_g(k),R_r(k),A_z(k)]′

式中，R_g为天波雷达运动目标在在雷达坐标系下的径向距，R_r为天波雷达运动目标在在雷达坐标系下的多普勒，A_z为天波雷达运动目标在在雷达坐标系下的面方位角，M＝3，分别为R_g、R_r、A_z；[R_g(k),R_r(k),A_z(k)]′为求[R_g(k),R_r(k),A_z(k)]的转置；

则考虑多径效应(多径效应指天波雷达的多路径传播效应,多径效应导致了有多个传播模式，两者在概念上等价)下天波雷达的***量测模型表示为：

y(k)＝H_i(x(k))+w_i(k),i＝1,...,n_k

其中，H_i(x)是非线性量测方程；w_i(k)是零均值的高斯白噪声，量测协方差矩阵为R_i(k)；n_k是时刻k传播模式的个数；

天波雷达过程噪声协方差矩阵Q为：

Q(k)＝E[v(k)v(k)']

式中，E[v(k)v(k)']为求v(k)v(k)'的期望；

量测协方差矩阵R为：

R_i(k)＝E[w_i(k)w_i(k)']

其它步骤及参数与具体实施方式一相同。

具体实施方式三：本实施方式与具体实施方式一或二不同的是：所述步骤二中根据步骤一获得的第k时刻天波雷达运动目标在地面坐标系下的运动状态向量x(k)、天波雷达状态误差协方差矩阵P(k)，计算天波雷达运动目标第k+1时刻的运动状态预测值和状态预测协方差矩阵P(k+1|k)，并由此得到天波雷达运动目标第k+1时刻传播模式i＝1,...,n_k+1在雷达坐标系下的量测预测值和量测预测协方差S_i(k+1)；具体为：

当k＝0时，采用传感器获得天波雷达运动目标在雷达坐标系下的测量值向量，初始化天波雷达运动目标在地面坐标系下的运动状态向量x(0)和天波雷达状态误差协方差矩阵P(0)；

当k＞0时，用k时刻的滤波估计值和滤波协方差P(k|k)初始化天波雷达运动目标在地面坐标系下的运动状态向量x(k)和天波雷达状态误差协方差矩阵P(k)，P(k)＝P(k|k)；计算第k+1时刻的运动状态预测值和状态预测协方差矩阵P(k+1|k)，并由此得到天波雷达运动目标第k+1时刻传播模式i＝1,...,n_k+1的量测预测值和量测预测协方差S_i(k+1)；计算方法具体为：

P(k+1|k)＝F(k)P(k|k)F(k)

S_i(k+1)＝J_i(k+1)P(k+1|k)J_i(k+1)'+R_i(k),i＝1,...,n_k+1

其中，J_i(k+1)为量测函数H_i(x)在当前估计值下的雅可比矩阵，J_i(k+1)'为J_i(k+1)的转置，n_k+1是时刻k+1传播模式的个数。

其它步骤及参数与具体实施方式一或二相同。

具体实施方式四：本实施方式与具体实施方式一至三之一不同的是：所述步骤三中传感器获取当前时刻天波雷达运动目标在雷达坐标系下的量测值，随机生成N个初始量测-模式关联对，即初始国家，并计算每个国家的权力c_t，t＝1，2，3，...，N，初始化当前迭代次数g＝1，预设算法迭代次数G；N取值为正整数；

具体过程为：

随机形成N个初始国家C＝{C₁,C₂,...,C_N}，并计算每个国家的权力c＝{c₁,c₂,...,c_N}，将国家权力c与似然函数d建立相关性c_n＝1/D_n；

其中c_n，D_n分别表示第n个国家的权力和对应的似然函数，n＝1，2，3，...，N；

似然函数的计算公式为

式中，d_zi表示第z个量测与第i个传播模式关联的似然函数，y_z(k)表示第k时刻第z个量测，表示第k时刻第i个传播模式的量测预测值，S_i(k)表示第k时刻第i个传播模式的新息协方差，为的转置。

在本发明中，每个国家代表一组量测-传播模式关联对。假设有三条传播模式A、B、C，另有三个量测a、b、c，则初始化形成的量测-传播模式关联对表示为(A,c)-(B,a)-(C,b)，则该国家的权力c＝1/D，其中，D＝d_cA+d_aB+d_bC；

d_cA表示第c个量测与第A个传播模式关联的似然函数，d_aB表示第a个量测与第B个传播模式关联的似然函数，d_bC表示第b个量测与第C个传播模式关联的似然函数。

其它步骤及参数与具体实施方式一至三之一相同。

具体实施方式五：本实施方式与具体实施方式一至四之一不同的是：所述步骤四中依据国家权力大小从高到低排序，选择前N_imp个作为殖民国家，其余N_col个作为殖民地，按殖民国家的权力大小将殖民地分配给初始殖民国家，形成初始帝国；具体为：

根据国家权力大小从高到低排序，选择前N_imp个作为殖民国家，殖民国家的权力定义为c_imp_f，f＝1,2,…,N_imp；

其余N_col个作为殖民地，殖民地的权力定义为c_col_j，j＝1,2,…,N_col；

然后对殖民国家权力进行标准化：

式中，c_imp_h为第h个殖民国家的权力，h＝1,2,…,N_imp；

根据公式N_imp_f＝round{p_imp_f×N_col}将所有殖民地分配给殖民国家，第f个殖民国家分得N_imp_f个殖民地；

其中，round{p_imp_f×N_col}为对p_imp_f×N_col取整；N_imp_f为殖民地个数。

其它步骤及参数与具体实施方式一至四之一相同。

具体实施方式六：本实施方式与具体实施方式一至五之一不同的是：所述步骤五中殖民地同化作用：殖民国家以自身量测-模式关联对(步骤三得到的)为基础，依据一定同化率随机地替换和改变殖民地路径中的α个量测-模式关联对；具体过程为：

步骤五一、定义ρ为同化率，并满足0＜ρ＜1，在殖民地每个传播模式的位置(已知的，例如(A,c)-(B,a)-(C,b)分别是(0.1,0.6,0.5)，(A,c)-(B,a)-(C,b)分别是[0,1]之间的随机数，分别执行步骤五二、步骤五三和步骤五四；

步骤五二、随机数大于等于ρ的殖民地传播模式的位置处发生同化；

同化方式为：用殖民国家中与随机数大于等于ρ的殖民地传播模式的位置处殖民地传播模式关联的量测作为同化后殖民地相同传播模式的关联量测；

步骤五三、随机数小于ρ的传播模式的位置处，若殖民地与随机数小于ρ的传播模式的关联量测在同化时未出现过，则视作被保持；

保持方式为：殖民地在随机数小于ρ的传播模式的位置处的量测-传播模式关联对直接作为同化后殖民地的量测-传播模式关联对；

步骤五四、随机数小于ρ的传播模式的位置处，若殖民地与随机数小于ρ的传播模式的关联量测已在同化时出现过，则将殖民国家在随机数小于ρ的传播模式的位置处的关联量测替换此位置的关联量测；

步骤五五、将同化后的殖民地替换原殖民地。

其它步骤及参数与具体实施方式一至五之一相同。

具体实施方式七：本实施方式与具体实施方式一至六之一不同的是：所述步骤六中殖民地革命：殖民地按革命概率P_m随机地选取殖民地两个位置的关联量测进行交换，然后计算殖民地革命前后的权力，保留权力大的殖民国家，并重新评估殖民国家和殖民地的权力，若后者大于前者，则取代前者成为新的帝国；具体为：

步骤六一、定义p_m为革命概率：

其中，p_imp_f为殖民地所属殖民国家的标准化权力，p_imp_max为所有殖民国家中最大的标准化权力，g为当前运算迭代数，G为预设算法的迭代次数，p₀为初始革命概率，0＜p₀＜1；

步骤六二、每个殖民地生成一个[0,1]之间的随机数，若该随机数大于p_m，则该殖民地发生革命；

革命方式为：随机将殖民地的两个传播模式的关联量测交换，直至每个殖民地的量测交换完，即完成革命操作；

步骤六三、每个殖民地完成革命操作后，重新评估帝国内的殖民国家和殖民地权力并进行排序，若殖民地的权力大于殖民国家，则该殖民地和殖民国家交换身份。

其它步骤及参数与具体实施方式一至六之一相同。

具体实施方式八：本实施方式与具体实施方式一至七之一不同的是：所述步骤七中帝国增强，根据每个殖民国家权力，对权力最小的国家进行增强，增强方式为：随机移除该殖民国家中β个关联量测，重新建立量测-模式关联对，重新计算殖民国家增强前后的权力，保留权力大的殖民国家；具体过程为：

根据每个殖民国家权力，对权力中权力最小的殖民国家进行增强操作，增强方式为：随机移除β个关联量测，从剩余量测集中选取β个量测***移除的位置，重新计算该殖民国家的权力，若替换后的殖民国家的权力大于原来的殖民国家的权力，则增强成功，否则，保留原殖民国家。

其它步骤及参数与具体实施方式一至七之一相同。

具体实施方式九：本实施方式与具体实施方式一至八之一不同的是：所述步骤八中殖民竞争，以帝国为单位，计算帝国内殖民国家和所有殖民地的标准化权力TC_n；通过评估所有帝国权力，将最弱帝国的最弱殖民地分配给最强帝国；当某个帝国内没有任何殖民地时，帝国消失，该帝国的殖民国家也被分配给最强帝国作为殖民地；具体为：

步骤八一、计算帝国权力，帝国的权力由组成该帝国的殖民国家的权力与殖民地平均权力加权的和构成：

其中，TC_f为第f个帝国的权力，c_col_j为第f个帝国第j个殖民地的权力，N_imp_f为当前该帝国殖民地个数，δ为权系数；c_imp_f为第f个殖民国家的权力，f＝1,2,…,N_imp；

步骤八二、对帝国权力进行标准化，第h个帝国的标准化权力TP_h表示为：

且满足建立帝国权力向量TP：

式中，NTC_h、NTC_f为中间变量(NTC_f和NTC_h是同一个东西，只不过避免混淆下标不一样)；TC_h为第h个帝国的权力，h＝1,2,…,N_imp；

步骤八三、建立一个维数与帝国数量相同的，由均匀分布随机数r_h组成的向量：

h＝1,2,…,N_imp；

其中，NI_h表示第h个帝国内拥有国家的数量，U为均匀分布，并定义概率向量:

假设d_h为D中所有元素的最大值，h≤N_imp，则把最弱帝国的最弱殖民地分配给第h个帝国，更新最弱和最强帝国的殖民地数量；

第h个帝国即是最强帝国；

步骤八四、判断是否存在无殖民地的帝国，若是，则该帝国的殖民国家作为殖民地分配给最强帝国，更新帝国即殖民国家的数量N_imp＝N_imp-1。

其它步骤及参数与具体实施方式一至八之一相同。

具体实施方式十：本实施方式与具体实施方式一至九之一不同的是：所述步骤十中根据步骤九输出的关联结果计算当前时刻地面坐标系下的目标状态估计值；具体为：

按照步骤九输出的量测-模式关联对更新目标的状态：

式中，表示第z个量测与第i个传播模式相关联；为第k+1时刻天波雷达运动目标在地面坐标系下的滤波估计值，P(k+1|k+1)为第k+1时刻天波雷达运动目标在地面坐标系下的滤波协方差。

其它步骤及参数与具体实施方式一至九之一相同。

采用以下实施例验证本发明的有益效果：

下面结合附图和实施例详细说明本发明的实施方式。附图和实施方式中的描述仅仅是示例性的，不是对本发明的保护范围的限制。

实施例

为展示本发明的实施方式的效果，设置一组仿真环境进行验证：

假设电离层是暂时平稳且不倾斜的，具有典型的两层分层结构——E层与F层。设E层电离层高度h_E＝100km，F层高度h_F＝220km。雷达发射站与接收站距离d＝100km。目标的初始状态设为：

x(0)＝[1100km,0.15km/s,0.10472rad,8,72665e-5rad/s]′

目标朝向雷达做径向运动。采样时间为T＝20s，采样点数为50个。

由于电离层分为两层，故存在四种传播模式：E层发射E层接收(EE模式)，E层发射F层接收(EF模式)，F层发射E层接收(FE模式)，F层发射F层接收(FF模式)，即n_k＝4。仿真环境按照四种传播模式分别仿真，每种传播模式的检测概率P_d＝0.4，将量测噪声标准差设置为

如图1所示，首先，在步骤一，建立运动目标的状态模型和量测模型。由于仿真目标在地面作匀速直线运动，因此用匀速运动(CV)模型对其建模，则目标的运动模型为：

x(k+1)＝F(k)x(k)+v(k)

式中，目标在地面坐标系下的运动状态向量：

其中，ρ(k)表示径向距离，表示径向距离变化率(径向速度)，b(k)表示方位角，表示方位角变化率，即L＝4，

***的状态转移矩阵：

v(k)是L维过程噪声，具有方差：

***量测模型为：

y(k)＝H_i(x(k))+w_i(k),i＝1,...,n_k

即

式中，雷达坐标系下的量测向量：

y(k)＝[R_g(k),R_r(k),A_z(k)]′

其中，R_g(k)表示径向距，R_r(k)表示多普勒，A_z(k)表示面方位角，即M＝3，

非线性量测方程：

w_i(k)是零均值的高斯白噪声，具有方差

接着，在步骤二，第k时刻的运动状态向量x(k)、状态误差协方差P(k)，

当k＝0时

当k＞0时

P(k)＝P(k|k)

计算第k+1时刻的运动状态预测值和状态预测协方差P(k+1|k)

P(k+1|k)＝F(k)P(k|k)F(k)

第k+1时刻传播模式i＝1,2,3,4的量测预测值和量测预测协方差S_i(k+1)

y(k)＝H_i(x(k))+w_i(k),i＝1,2,3,4

S_i(k+1)＝J_i(k+1)P(k+1|k)J_i(k+1)'+R_i(k),i＝1,2,3,4

式中，J_i(k+1)量测函数H_i(x)在当前估计值下的雅可比矩阵

忽略时间k可以表示为：

式中，

η＝ρ-dsin(b)

然后，在步骤三，产生初始帝国，随机形成N个初始国家C＝{C₁,C₂,...,C_N}，在本实施例中，在本发明中，每个国家代表一组量测-传播模式关联对。假设有四个量测a、b、c，d，则初始化形成的量测-传播模式关联对可表示为(EE,c)-(EF,a)-(FE,b)-(FF,d)，即一个国家。其中，(EE,c)表示量测c与传播模式EE相关联，计算每个国家的权力c＝{c₁,c₂,...,c_N}，c_n＝1/D_n，其中D_n表示第n个国家对应的似然函数之和：

式中，d_jq表示第j个量测与第q个传播模式关联的似然函数，初始化当前迭代次数g，预设算法迭代次数G。

之后，在步骤四，依据国家权力大小从高到低排序，选择N_imp个作为殖民国家，它们的权力定义为c_imp_i，i＝1,2,…,N_imp其余N_col个作为殖民地，它们的权力定义为c_col_j，j＝1,2,…,N_col；然后对殖民国家权力进行标准化：

根据公式N_imp_i＝round{p_imp_i×N_col}将所有殖民地分配给殖民国家，第n个殖民国家分得N_imp_i个殖民地；

其中，round{p_imp_i×N_col}为对p_imp_i×N_col取整；N_imp_i为殖民地个数；

然后，在步骤五，发生殖民地同化作用，定义同化率ρ＝0.5，并在殖民地每个传播模式的位置生成一个[0,1]之间的随机数，若随机数大于等于ρ，则发生同化，殖民国家中与该传播模式关联的量测作为同化后殖民地相同传播模式的关联量测，若随机数小于ρ，且在殖民地与这些传播模式的关联量测在同化时未出现过，则视作被保持，保持方式为殖民地在这些位置处的量测-传播模式关联对直接作为同化后殖民地的量测-传播模式关联对，若随机数小于ρ，但殖民地与这些传播模式的关联量测已在同化时出现过，则此位置也发生同化，将殖民国家在这些位置处的关联量测替换此位置的关联量测，最后，将同化后的殖民地替换原殖民地。

然后，在步骤六，发生殖民地革命，定义p_m为革命概率：

其中，p_imp_i为殖民地所属殖民国家的标准化权力，p_imp_max所有殖民国家中最大的标准化权力，g为当前运算迭代数，G为总共迭代次数，p₀为初始革命概率，定义p₀＝0.3，每个殖民地生成一个[0,1]之间的随机数，若该随机数大于p_m，则该殖民地发生革命，革命方式为随机将殖民地的两个传播模式的关联量测交换，每个殖民地完成革命操作后，重新评估帝国内的殖民国家和殖民地权力并进行排序，若殖民地的最高权力大于殖民国家，则该殖民地和殖民国家交换身份。

接着，在步骤七，进行帝国增强，根据殖民国家权力表，对权力表中权力最小的殖民国家进行增强操作，增强方式为随机移除m个关联量测，从剩余量测集中选取相同数量的量测***移除的位置，重新计算该殖民国家的权力，若替换后的权力大于原来的权力，则增强成功，否则，保留原殖民国家。

然后，在步骤八，发生帝国竞争，竞争过程为：

(1)计算帝国权力，帝国的权力由组成该帝国的殖民国家的权力与殖民地平均权力加权和构成：

其中，TC_f为第f个帝国的权力，c_col_j为第f个帝国第j个殖民地的权力，N_imp_f为当前该帝国殖民地个数，δ为权系数；c_imp_f为第f个殖民国家的权力，f＝1,2,…,N_imp；

(2)对帝国权力进行标准化，第h个帝国的标准化权力TP_h表示为：

且满足建立帝国权力向量TP：

式中，NTC_h、NTC_f为中间变量；TC_h为第h个帝国的权力，h＝1,2,…,N_imp；

(3)建立一个维数与帝国数量相同的，由均匀分布随机数r_h组成的向量：

h＝1,2,…,N_imp；

其中，NI_h表示第h个帝国内拥有国家的数量，U为均匀分布，并定义概率向量:

假设d_h为D中所有元素的最大值，h≤N_imp，则把最弱帝国的最弱殖民地分配给第h个帝国，更新最弱和最强帝国的殖民地数量；

第h个帝国即是最强帝国；

(4)判断是否存在无殖民地的帝国，若是，则该帝国的殖民国家作为殖民地分配给最强帝国，更新帝国和殖民国家的数量N_imp＝N_imp-1。

然后，在步骤九，判断当前帝国数量是否为1或当前算法迭代次数g≥G，若是，算法迭代结束，输出最终的量测-模式关联对作为数据关联的结果，若否，返回步骤五，并令g＝g+1。

然后，在步骤十，根据步骤九输出的关联结果计算当前时刻地面坐标系下的目标状态估计值，计算公式为：

式中，表示第z个量测与第i个传播模式相关联；为第k+1时刻天波雷达运动目标在地面坐标系下的滤波估计值，P(k+1|k+1)为第k+1时刻天波雷达运动目标在地面坐标系下的滤波协方差。

最后，在步骤十一，返回执行所述步骤二，并使k＝k+1，直到目标跟踪结束。

选取均方根误差(RMSE)作为性能的量度，进行100次蒙特卡罗仿真实验。

仿真结果如图2到图11所示

图2到图4给出单次仿真的径向距、多普勒和方向角典型量测图，由图可见，对于同一目标，其四种传播模式的量测非常接近，并存在相互干扰的情况，而且每条传播模式的检测概率都很低，这也尽一步增大了数据关联的难度。

图5给出地面坐标系下目标跟踪航迹图，由图3可见，随着跟踪时间的增加，跟踪结果向着真实值逐渐趋于稳定，消除了量测噪声和杂波的影响，达到了收敛的效果，证明了本发明的有效性。

图6图7分别给出地面坐标系下目标径向距离变化率和方位角变化率的跟踪结果图，由图中可见，同样地，随着跟踪时间的增加，跟踪结果均向着真实值逐渐趋于稳定，消除了量测噪声和杂波的影响，达到了收敛的效果，证明了本发明的有效性。

图8到图11给出50次蒙特卡洛仿真后目标各维度目标各维度的均方根误差，从图中可以看出，随着跟踪的增加，目标径向距离、径向距离变化率、方位角、方位角变化率的均方根误差均逐渐减少并趋近于0，跟踪精度提高10倍左右，同样证明本发明可以有效解决天波雷达多路径数据关联问题，达到了收敛的效果。

本发明还可有其它多种实施例，在不背离本发明精神及其实质的情况下，本领域技术人员当可根据本发明作出各种相应的改变和变形，但这些相应的改变和变形都应属于本发明所附的权利要求的保护范围。

Claims (10)

1.一种基于帝国竞争算法的天波雷达多路径数据关联方法，其特征在于：

步骤一，建立天波雷达运动目标的运动模型和多路径量测模型，将天波雷达运动目标地面坐标系下的运动状态标记为L维的运动状态向量x，天波雷达运动目标雷达坐标系下的量测值标记为M维的量测向量y，根据运动模型获得天波雷达状态误差协方差矩阵P和天波雷达过程噪声协方差矩阵Q，根据量测模型获得量测协方差矩阵R；

步骤二，根据步骤一获得的第k时刻天波雷达运动目标在地面坐标系下的运动状态向量x(k)、天波雷达状态误差协方差矩阵P(k)，计算天波雷达运动目标第k+1时刻的运动状态预测值和状态预测协方差矩阵P(k+1|k)，并由此得到天波雷达运动目标第k+1时刻传播模式i＝1,...,n_k+1在雷达坐标系下的量测预测值和量测预测协方差S_i(k+1)；n_k是时刻k传播模式的个数；

步骤三，传感器获取当前时刻天波雷达运动目标在雷达坐标系下的量测值，随机生成N个初始量测-模式关联对，即初始国家，并计算每个国家的权力c_t，t＝1，2，3，...，N，初始化当前迭代次数g＝1，预设算法迭代次数G；N取值为正整数；

步骤四，依据国家权力大小从高到低排序，选择前N_imp个作为殖民国家，其余N_col个作为殖民地，按殖民国家的权力大小将殖民地分配给初始殖民国家，形成初始帝国；

步骤五，殖民地同化作用：殖民国家以自身量测-模式关联对为基础，依据一定同化率随机地替换和改变殖民地路径中的α个量测-模式关联对；

步骤六，殖民地革命：殖民地按革命概率P_m随机地选取殖民地两个位置的关联量测进行交换，然后计算殖民地革命前后的权力，保留权力大的殖民国家，并重新评估殖民国家和殖民地的权力，若后者大于前者，则取代前者成为新的帝国；

步骤七，帝国增强，根据每个殖民国家权力，对权力最小的国家进行增强，增强方式为：随机移除该殖民国家中β个关联量测，重新建立量测-模式关联对，重新计算殖民国家增强前后的权力，保留权力大的殖民国家；

步骤八，殖民竞争，以帝国为单位，计算帝国内殖民国家和所有殖民地的标准化权力TC_n；通过评估所有帝国权力，将最弱帝国的最弱殖民地分配给最强帝国；当某个帝国内没有任何殖民地时，帝国消失，该帝国的殖民国家也被分配给最强帝国作为殖民地；

步骤九，若达到预设算法的迭代次数G或只剩一个帝国，算法结束，输出最终的量测-模式关联对作为数据关联的结果，否则，令g＝g+1，返回步骤五；

步骤十，根据步骤九输出的关联结果计算当前时刻地面坐标系下的目标状态估计值；

步骤十一，判断目标跟踪是否结束，如果没有结束，使k＝k+1，执行步骤二～步骤十，如果结束则完成数据关联方法。

2.根据权利要求1所述一种基于帝国竞争算法的天波雷达多路径数据关联方法，其特征在于：所述步骤一中建立天波雷达运动目标的运动模型和多路径量测模型，将天波雷达运动目标地面坐标系下的运动状态标记为L维的运动状态向量x，天波雷达运动目标雷达坐标系下的量测值标记为M维的量测向量y，根据运动模型获得天波雷达状态误差协方差矩阵P和天波雷达过程噪声协方差矩阵Q，根据量测模型获得量测协方差矩阵R；具体过程为：

在采样时刻k，天波雷达运动目标在地面坐标系下的运动状态向量表示为：

其中，ρ(k)表示天波雷达运动目标在地面坐标系下的径向距离，表示天波雷达运动目标在地面坐标系下的径向距离变化率，b(k)表示天波雷达运动目标在地面坐标系下的方位角，表示天波雷达运动目标在地面坐标系下的方位角变化率，L＝4，分别为ρ(k)、b(k)、 为求的转置；

则天波雷达运动目标运动状态方程表示为：

x(k+1)＝F(k)x(k)+v(k)

其中，F(k)是L×L阶的状态转移矩阵，L维过程噪声v(k)是一个零均值的高斯白噪声，v(k)协方差矩阵为Q(k)；

在采样时刻k，天波雷达运动目标一次量测值在雷达坐标系下的量测向量表示为：

y(k)＝[R_g(k),R_r(k),A_z(k)]′

式中，R_g为天波雷达运动目标在在雷达坐标系下的径向距，R_r为天波雷达运动目标在在雷达坐标系下的多普勒，A_z为天波雷达运动目标在在雷达坐标系下的面方位角，M＝3，分别为R_g、R_r、A_z；[R_g(k),R_r(k),A_z(k)]′为求[R_g(k),R_r(k),A_z(k)]的转置；

则考虑多径效应下天波雷达的***量测模型表示为：

y(k)＝H_i(x(k))+w_i(k),i＝1,...,n_k

其中，H_i(x)是非线性量测方程；w_i(k)是零均值的高斯白噪声，量测协方差矩阵为R_i(k)；n_k是时刻k传播模式的个数；

天波雷达过程噪声协方差矩阵Q为：

Q(k)＝E[v(k)v(k)']

式中，E[v(k)v(k)']为求v(k)v(k)'的期望；

量测协方差矩阵R为：

R_i(k)＝E[w_i(k)w_i(k)']。

3.根据权利要求2所述一种基于帝国竞争算法的天波雷达多路径数据关联方法，其特征在于：所述步骤二中根据步骤一获得的第k时刻天波雷达运动目标在地面坐标系下的运动状态向量x(k)、天波雷达状态误差协方差矩阵P(k)，计算天波雷达运动目标第k+1时刻的运动状态预测值和状态预测协方差矩阵P(k+1|k)，并由此得到天波雷达运动目标第k+1时刻传播模式i＝1,...,n_k+1在雷达坐标系下的量测预测值和量测预测协方差S_i(k+1)；具体为：

当k＝0时，采用传感器获得天波雷达运动目标在雷达坐标系下的测量值向量，初始化天波雷达运动目标在地面坐标系下的运动状态向量x(0)和天波雷达状态误差协方差矩阵P(0)；

当k＞0时，用k时刻的滤波估计值和滤波协方差P(k|k)初始化天波雷达运动目标在地面坐标系下的运动状态向量x(k)和天波雷达状态误差协方差矩阵P(k)，P(k)＝P(k|k)；计算第k+1时刻的运动状态预测值和状态预测协方差矩阵P(k+1|k)，并由此得到天波雷达运动目标第k+1时刻传播模式i＝1,...,n_k+1的量测预测值和量测预测协方差S_i(k+1)；计算方法具体为：

P(k+1|k)＝F(k)P(k|k)F(k)

S_i(k+1)＝J_i(k+1)P(k+1|k)J_i(k+1)'+R_i(k),i＝1,...,n_k+1

其中，J_i(k+1)为量测函数H_i(x)在当前估计值下的雅可比矩阵，J_i(k+1)'为J_i(k+1)的转置，n_k+1是时刻k+1传播模式的个数。

4.根据权利要求3所述一种基于帝国竞争算法的天波雷达多路径数据关联方法，其特征在于：所述步骤三中传感器获取当前时刻天波雷达运动目标在雷达坐标系下的量测值，随机生成N个初始量测-模式关联对，即初始国家，并计算每个国家的权力c_t，t＝1，2，3，...，N，初始化当前迭代次数g＝1，预设算法迭代次数G；N取值为正整数；

具体过程为：

随机形成N个初始国家C＝{C₁,C₂,...,C_N}，并计算每个国家的权力c＝{c₁,c₂,...,c_N}，将国家权力c与似然函数d建立相关性c_n＝1/D_n；

其中c_n，D_n分别表示第n个国家的权力和对应的似然函数，n＝1，2，3，...，N；

似然函数的计算公式为

式中，d_zi表示第z个量测与第i个传播模式关联的似然函数，y_z(k)表示第k时刻第z个量测，表示第k时刻第i个传播模式的量测预测值，S_i(k)表示第k时刻第i个传播模式的新息协方差，为的转置。

5.根据权利要求4所述一种基于帝国竞争算法的天波雷达多路径数据关联方法，其特征在于：所述步骤四中依据国家权力大小从高到低排序，选择前N_imp个作为殖民国家，其余N_col个作为殖民地，按殖民国家的权力大小将殖民地分配给初始殖民国家，形成初始帝国；具体为：

根据国家权力大小从高到低排序，选择前N_imp个作为殖民国家，殖民国家的权力定义为c_imp_f，f＝1,2,…,N_imp；

其余N_col个作为殖民地，殖民地的权力定义为c_col_j，j＝1,2,…,N_col；

然后对殖民国家权力进行标准化：

式中，c_imp_h为第h个殖民国家的权力，h＝1,2,…,N_imp；

根据公式N_imp_f＝round{p_imp_f×N_col}将所有殖民地分配给殖民国家，第f个殖民国家分得N_imp_f个殖民地；

其中，round{p_imp_f×N_col}为对p_imp_f×N_col取整；N_imp_f为殖民地个数。

6.根据权利要求5所述一种基于帝国竞争算法的天波雷达多路径数据关联方法，其特征在于：所述步骤五中殖民地同化作用：殖民国家以自身量测-模式关联对为基础，依据一定同化率随机地替换和改变殖民地路径中的α个量测-模式关联对；具体过程为：

步骤五一、定义ρ为同化率，并满足0＜ρ＜1，在殖民地每个传播模式的位置生成一个[0,1]之间的随机数，分别执行步骤五二、步骤五三和步骤五四；

步骤五二、随机数大于等于ρ的殖民地传播模式的位置处发生同化；

同化方式为：用殖民国家中与随机数大于等于ρ的殖民地传播模式的位置处殖民地传播模式关联的量测作为同化后殖民地相同传播模式的关联量测；

步骤五三、随机数小于ρ的传播模式的位置处，若殖民地与随机数小于ρ的传播模式的关联量测在同化时未出现过，则视作被保持；

保持方式为：殖民地在随机数小于ρ的传播模式的位置处的量测-传播模式关联对直接作为同化后殖民地的量测-传播模式关联对；

步骤五四、随机数小于ρ的传播模式的位置处，若殖民地与随机数小于ρ的传播模式的关联量测已在同化时出现过，则将殖民国家在随机数小于ρ的传播模式的位置处的关联量测替换此位置的关联量测；

步骤五五、将同化后的殖民地替换原殖民地。

7.根据权利要求6所述一种基于帝国竞争算法的天波雷达多路径数据关联方法，其特征在于：所述步骤六中殖民地革命：殖民地按革命概率P_m随机地选取殖民地两个位置的关联量测进行交换，然后计算殖民地革命前后的权力，保留权力大的殖民国家，并重新评估殖民国家和殖民地的权力，若后者大于前者，则取代前者成为新的帝国；具体为：

步骤六一、定义p_m为革命概率：

其中，p_imp_f为殖民地所属殖民国家的标准化权力，p_imp_max为所有殖民国家中最大的标准化权力，g为当前运算迭代数，G为预设算法迭代次数，p₀为初始革命概率，0＜p₀＜1；

步骤六二、每个殖民地生成一个[0,1]之间的随机数，若该随机数大于p_m，则该殖民地发生革命；

革命方式为：随机将殖民地的两个传播模式的关联量测交换，直至每个殖民地的量测交换完，即完成革命操作；

步骤六三、每个殖民地完成革命操作后，重新评估帝国内的殖民国家和殖民地权力并进行排序，若殖民地的权力大于殖民国家，则该殖民地和殖民国家交换身份。

8.根据权利要求7所述一种基于帝国竞争算法的天波雷达多路径数据关联方法，其特征在于：所述步骤七中帝国增强，根据每个殖民国家权力，对权力最小的国家进行增强，增强方式为：随机移除该殖民国家中β个关联量测，重新建立量测-模式关联对，重新计算殖民国家增强前后的权力，保留权力大的殖民国家；具体过程为：

根据每个殖民国家权力，对权力中权力最小的殖民国家进行增强操作，增强方式为：随机移除β个关联量测，从剩余量测集中选取β个量测***移除的位置，重新计算该殖民国家的权力，若替换后的殖民国家的权力大于原来的殖民国家的权力，则增强成功，否则，保留原殖民国家。

9.根据权利要求8所述一种基于帝国竞争算法的天波雷达多路径数据关联方法，其特征在于：所述步骤八中殖民竞争，以帝国为单位，计算帝国内殖民国家和所有殖民地的标准化权力TC_n；通过评估所有帝国权力，将最弱帝国的最弱殖民地分配给最强帝国；当某个帝国内没有任何殖民地时，帝国消失，该帝国的殖民国家也被分配给最强帝国作为殖民地；具体为：

步骤八一、计算帝国权力，帝国的权力由组成该帝国的殖民国家的权力与殖民地平均权力加权的和构成：

其中，TC_f为第f个帝国的权力，c_col_j为第f个帝国第j个殖民地的权力，N_imp_f为当前该帝国殖民地个数，δ为权系数；c_imp_f为第f个殖民国家的权力，f＝1,2,…,N_imp；

步骤八二、对帝国权力进行标准化，第h个帝国的标准化权力TP_h表示为：

且满足建立帝国权力向量TP：

式中，NTC_h、NTC_f为中间变量；TC_h为第h个帝国的权力，h＝1,2,…,N_imp；

步骤八三、建立一个维数与帝国数量相同的，由均匀分布随机数r_h组成的向量：

h＝1,2,…,N_imp；

其中，NI_h表示第h个帝国内拥有国家的数量，U为均匀分布，并定义概率向量:

假设d_h为D中所有元素的最大值，h≤N_imp，则把最弱帝国的最弱殖民地分配给第h个帝国，更新最弱和最强帝国的殖民地数量；

第h个帝国即是最强帝国；

步骤八四、判断是否存在无殖民地的帝国，若是，则该帝国的殖民国家作为殖民地分配给最强帝国，更新帝国即殖民国家的数量N_imp＝N_imp-1。

10.根据权利要求9所述一种基于帝国竞争算法的天波雷达多路径数据关联方法，其特征在于：所述步骤十中根据步骤九输出的关联结果计算当前时刻地面坐标系下的目标状态估计值；具体为：

按照步骤九输出的量测-模式关联对更新目标的状态：

式中，表示第z个量测与第i个传播模式相关联；为第k+1时刻天波雷达运动目标在地面坐标系下的滤波估计值，P(k+1|k+1)为第k+1时刻天波雷达运动目标在地面坐标系下的滤波协方差。