CN109239704A - 一种基于序贯滤波交互式多模型的自适应采样方法 - Google Patents

Abstract

本发明公开了一种基于序贯滤波交互式多模型的自适应采样方法，涉及相控阵雷达的时间资源管理领域，主要针对具有多普勒量测的相控阵雷达机动目标跟踪时采样周期的自适应问题。该方法将预测协方差门限法嵌入到序贯滤波交互式多模型方法中，模型的更新概率由各模型的预测概率采用基于位置量测的滤波器组和基于伪量测的滤波器组输出的预测模型概率的均值获得。该方法以交互式多模型算法为基本框架，各模型滤波器采用序贯相接的结构完成对位置量测与多普勒量测的处理。

Description

一种基于序贯滤波交互式多模型的自适应采样方法

技术领域

本发明涉及相控阵雷达的时间资源管理领域，主要针对具有多普勒量测的相控阵雷达机动目标跟踪时采样周期的自适应问题。

背景技术

在现代军事战场及民用领域内，相控阵雷达作为一种先进的有源电扫阵列多功能设备得到广泛的应用，区别于以前的机械雷达靠天线的机械转动来实现波束扫描的传统方式，其采用电子方式改变天线的工作频率或相位从而改变雷达波束的指向来完成扫描(见文献：张光义，赵玉洁.相控阵雷达技术[M].北京：电子工业出版社,2006)。由于相控阵雷达具有灵活迅捷的特点，可对其系统资源进行灵活的配置。其中对相控阵雷达的时间资源管理，关键就是控制目标跟踪时的采样周期。目前对时间资源管理也诸多相关的研究成果：Cohen提出位置残差控制采样周期的方法，利用递推公式来反映采样周期的前后变换情况(见文献：Cohen S A.Adaptive variable update rate algorithm for trackingtargets with a phased array radar[J].IEE Proceedings F-Communications,Radarand Signal Processing,2008,133(3):277-280.)； Van Keuk提出了公式法给出了采样间隔与机动参数的函数关系，并利用期望精度来控制采样周期。(见文献：Van Keuk G,Blackman S S.On phased-array radar tracking and parameter control[J].IEEETransactions on Aerospace&Electronic Systems Aes,1993, 29(1):186-194.)；基于协方差控制的预测误差协方差门限法通过将目标的预测误差标准差与所设门限进行比较，筛选出满足条件的采样周期(见文献：Watson G A,Blair W D.Tracking performance of aphased array radar with revisit time controlled using the IMM algorithm[C]Radar Conference,1994.Record of the 1994 IEEE National.IEEE, 1994:160-165.)。Benoudnine在协方差门限的基础上提出了一种快速自适应IMM算法(Fast AdaptationIMM,FAIMM)，其定义了采样间隔的最大和最小值，当目标机动时选择较小的采样周期，无机动时选择最大的采样周期，目标的机动与否由IMM输出的模型概率确定(见文献：Benoudnine H,Keche M,Ouamri A,et al.Fast Adaptive Update Rate for PhasedArray Radar Using IMM Target Tracking Algorithm[C]IEEE InternationalSymposium on Signal Processing and Information Technology.IEEE,2007:277-282.)。

由上述自适应算法方法可知，采样周期的变化是基于目标跟踪滤波器输出的各参数来确定的。因此，高性能的目标跟踪方法是自适应采样周期的基础，将自适应采样周期的思想融入选取的目标跟踪算法，才能实现有效的自适应采样周期目标跟踪，对于相控阵雷达目标跟踪，其目标跟踪算法要综合考虑两个方面的问题，即量测的非线性特性和目标的机动特性。首先，相控阵雷达的量测参数包括了距离、角度以及径向速度等，但这些量测参数与目标的状态之间往往存在着非线性关系。因此，在进行目标跟踪时需要解决量测与状态的非线性问题。最常见的一种方法是采用量测转换卡尔曼滤波算法(ConvertedMeasurement Kalman Filter,CMFK)。其中，传统转换量测卡尔曼滤波算法(ConventionalCMKF)最早被提出。针对传统方法存在偏差较大的问题，去偏量测转换(Debiased CMKF)算法和无偏量测转换(Unbiased CMKF)相继被提出(见文献Lerro D,Bar-Shalom Y.Trackingwith debiased consistent converted measurements versus EKF[J].Aerospace&Electronic Systems IEEE Transactions on,1993, 29(3):1015-1022；Song X,Zhou Y,Bar-Shalom Y.Unbiased converted measurements for tracking[J].Aerospace&Electronic Systems IEEE Transactions on,1998, 34(3):1023-1027.)二者保证了量测转换的无偏性，但UCMKF在对转换误差推导时存在前后不兼容问题。改进的算法MUCMKF(Modified Unbiased CMKF)对此进行了修正(见文献：Duan Z, Han C,Li X R.Commentson"Unbiased converted measurements for tracking"[J].Aerospace &ElectronicSystems IEEE Transactions on,2004,40(4):1374.)。但以上算法均采用量测值为先验条件进行推导，间接地引入了量测噪声使得滤波估计与量测误差之间存在相关性。

发明内容

本发明主要解决相控阵雷达自适应时间资源管理问题，提出了一种基于序贯滤波交互式多模型的自适应采样方法。该方法将预测协方差门限法嵌入到序贯滤波交互式多模型方法中，模型的更新概率由各模型的预测概率采用基于位置量测的滤波器组和基于伪量测的滤波器组输出的预测模型概率的均值获得。该方法以交互式多模型算法为基本框架，各模型滤波器采用序贯相接的结构完成对位置量测与多普勒量测的处理。

本发明的技术方案为一种基于序贯滤波交互式多模型的自适应采样方法,该方法包括：

步骤1：输入交互估计

首先假设IMM交互式多模型系统中存在N个模型，运动模型集合表示为

式中，M_k表示t_k时刻的运动模型；

按下式计算各模型滤波器输入的状态估计和协方差P^0j t_k-1：

其中，表示在t_k-1时刻模型i的状态估计，表示第l个模型的误差协方差，表示模型混合概率，其计算方式如下：

式中，π_ij表示由模型i转移到模型j的概率；μⁱ(t_k-1)表示t_k-1时刻运动模型i的更新概率， C^j表示模型j的归一化常数；

步骤2：计算各模型位置状态预测和预测协方差

设置一组采样间隔集合其中包含了具有n个不同取值的元素，从大至小选取采样间隔，选择第l个采样周期T_l，代入下式

其中，和P^0j(t_k-1)分别为步骤1中获得的位置估计输入估计与协方差交互输入估计的初始值；F^j,G^j表示模型j的状态转移矩阵与噪声驱动矩阵，Q^j(t_k-1)表示上一时刻的过程噪声自相关矩阵；

步骤3：计算预测位置协方差；

以二维场景为例，预测位置协方差在极坐标下表示为

其中，量测矩阵C＝diag{[1 0 0],2}_2×6，极坐标系下的变换矩阵J_p中的各元素由对应位置偏导数构成

表示采样周期T_l为时对应的x方向上的位置状态，表示表示采样周期为 T_l时对应的y方向上的位置状态；

组合预测协方差计算方法如下：

式中，将t_k时刻的模型j预测概率表示为可由各位置滤波器组的预测概率与伪量测预测概率均值获得，其计算方式如下：

式中，π_ij为概率转移矩阵，分别表示模型j的各位置量测滤波器与伪量测滤波器的模型预测概率，μ^p,i t_k-1，μ^ε,it_k-1分别表示上一时刻模型i的各位置量测滤波器与伪量测滤波器的模型概率；

进而可以获得采样周期为T_l时所对应的距离的预测误差标准差σ_r(T_l)和方位角的预测误差标准差σ_b(T_l)为：

和P_pol(T_l)_2,2的下标(n,n),n＝1,2,3表示矩阵中行和列对应位置；

步骤4：计算协方差门限值P_th：

距离对应的门限预测误差标准差σ_r,th与方位角对应的门限预测误差标准差σ_b,th计算方式如下：

其中，u_0.5a是置信度为P_CL时对应的双侧分位数；为预测方位角对应的波束宽度， P_CL表示概率置信度，L_g表示距离波门宽度，计算方式如下：

式中，M为阵元总数；λ为发射信号的波长；d为各天线单元间距离；为预测方位角；

步骤5：当且仅当T_l满足如下关系

σ_r(T_l)≤σ_r,thσ_b(T_l)≤σ_b,th (13)

选取第l个采样间隔T_l，作为t_k时刻的采样间隔即T(t_k)＝T_l；

否则继续取出中的下一个元素进行测试，直至选取到最小的采样周期T_min；定义上一时刻为t_k-1，则t_k时刻点为t_k＝t_k-1+T(t_k)；

步骤6：利用步骤5中获取的采样周期T(t_k)，产生t_k时刻的量测位置量测值和伪量测

其中，Z_c(t_k)表示t_k时刻的量测转换值；

非线性变换函数x(t_k)表示x方向上的位置状态，表示x方向上的速度，y(t_k)表示y方向上的位置状态，表示y方向上的速度，量测转换误差目标在直角坐标系下的状态为X(t_k)；分别为距离量测、方位角量测和多普勒量测；对应的量测噪声为方差分别为且和不相关，和的相关系数为ρ；λ为偏差补偿因子；

步骤7：基于预测值计算量测转换误差协方差矩阵R^j(t_k)；为简洁表述矩阵中的元素，约定t_k时刻所对应的参量，以下标k的形式表示

协方差矩阵R^j中的各元素的计算方法如下：

其中，极坐标系下的预测状态r_t,θ_t,通过步骤2中预测状态获得，参见式(22)；预测误差方差通过预测误差协方差矩阵的元素获得，参见式(23)。

式中，为中对应的目标预测位置与速度元素，矩阵J的形式如下

步骤8：各模型序贯滤波(SQ Filtering)

将步骤6，7中获得t_k时刻的量测转换值Z_c(t_k)和量测误差协方差R^j(t_k)及步骤2中获得的t_k-1时刻的模型j的输入估计X^0j t_k-1,P^0j t_k-1，代入当前模型的序贯滤波器中。具体过程可参见参考文献(Duan,Zhansheng,C.Han,and X.R.Li."Sequential NonlinearTracking Filter with Range-rate Measurements in Spherical Coordinates."InProceedings of the 7 th International Conference on Information Fusion2004:599--605.)

步骤9：模型概率更新

模型j在t_k时刻的更新概率μ^j(t_k)表示为位置量测模型概率μ^p,j(t_k)和伪量测模型概率μ^ε,j(t_k)的均值

其中，上标j表示运动模型，p表示与位置量测相关的信息，ε表示与伪量测相关的信息；各滤波器的似然函数Λ^p,j,Λ^ε,j计算公式如下：

式中e^p,j为量测残差，S^p,j为残差的自相关矩阵；C^p,j和C^ε,j分别表示模型j的位置量测和伪量测归一化常数；

步骤10：状态估计融合

将步骤7中获得的t_k时刻的序贯滤波估计结果与模型更新概率μ^j(t_k)进行融合

步骤11：将各模型的滤波结果代回到步骤1的(2)式中，计算下一时刻 t_k+1的各模型滤波器输入的状态估计和协方差。

从上述步骤可见，本发明保证雷达能处理目标位置和多普勒信息的同时，能够达到时间资源合理分配的目的。算法采用满足协方差门限的最大采样周期确定下一时刻采样点，合理控制量测点迹的录取，减少了雷达在时间资源上的消耗。IMM各嵌套子滤波器采用序贯相接结构来处理目标的速度信息，其中基于预测值作量测转换方法消除了量测噪声所带来的影响；各模型的更新概率由位置量测和伪量测二者的模型概率的均值共同确定，充分考虑了位置量测与伪量测因素的影响，保证了模型切换时的准确性。可见，ASP-SQ-IMM算法能够实现对带非线性量测信息进行机动目标跟踪的同时也对时间资源进行了合理配置。值得一提，上述阐述过程是以二维平面内的跟踪问题为背景进行的，但该思想同样可用于三维坐标系下的机动目标跟踪。

附图说明

图1为本发明ASP-SQ-IMM算法结构框图；

图2为本发明自适应采样周期算法流程图；

图3为本发明交互式多模型算法流程图；

图4为本发明具体实施方式目标真实航迹；

图5为本发明具体实施方式中场景1采样周期变化、模型概率转移曲线；

图6为本发明具体实施方式中场景2采样周期变化、模型概率转移曲线；

图7为本发明具体实施方式中场景1位置RMSE曲线；

图8为本发明具体实施方式中场景2位置RMSE曲线；

具体实施方式

在此次实施方案中，考虑一部多普勒雷达对二维平面内的某一机动运动目标进行跟踪。假设目标初始值位置为x＝6000m,y＝5000m，设v_x＝350m/s,v_y＝0m/s为初始速度，分别在1～60s、 75～100s和115～200s间进行匀速运动，在60～75s和100～115s时进行匀加速运动，60～75s时，两个方向的加速度分别为：a_x＝-23m/s²,a_y＝-25m/s²；在100～115s时两个方向的加速度分别为：24m/s²,a_y＝25m/s²。量测噪声的统计特性由极坐标系下的方位角、距离及径向速度标准差描述，采样周期集合为T_se＝{0.1,0.2,…,2.5},采用ASP-SQ-IMM算法实现对该目标的跟踪，模型集合由匀速运动(CV)和匀加速运动(CA)模型构成，各模型初始概率均为0.5，概率转移矩阵为：其他参数包括：阵元总数M＝2048，跟踪总时长T_tot＝200s，距离波门宽度Lg＝1575m，最小采样间隔T_min＝0.1s，双侧分位数(置信度为0.99时)

u_0.5a＝2.5758(P_CL＝0.99)，相关系数ρ＝0.9。

两个场景下的量测误差如下：

场景1：传感器的距离量测误差为方位角的量测误差为度，多普勒径向速度的量测误差为0.2m/s；

场景2：传感器的距离量测误差为方位角的量测误差为度，多普勒径向速度的量测误差为0.5m/s；

在上述两个场景下，平均采样间隔及平均位置均方根误差的计算公式为：

其中，N为蒙特卡洛仿真总次数，a,b表示划分区间的起始点和终点，为第j次仿真中的t_k时刻的采样间隔。为第j次仿真中的t_k时刻的估计误差，计数区间总数：

上述平滑处理选出在平滑间隔中有采样点的蒙特卡罗仿真，先对每次仿真计算出该间隔的平均采样周期，再把这些值在满足式(28)的仿真之间进行平均，仿真中平滑间隔取为4s，仿真次数为500。

ASP-SQ-IMM算法在自适应时间资源分配时的采样周期、模型概率及加速度变化结果如图3、图4所示。从图中可以看出，目标在整个运动期间发生了两次大的机动，CA、CV模型概率的变化反映了ASP-SQ-IMM跟踪算法能够有效地检测出这些机动，并能够与加速度的变化相对应。采样间隔曲线的变化表明自适应时间资源分配策略会在目标机动期间自动使用较小的采样间隔对目标进行跟踪，按照目标发生机动次数，采样间隔(周期)曲线出现2个凹口，正如图中采样周期变化所示，对比图4知，在量测误差较大的场景中也能够实现自适应时间资源分配。

针对自适应时间资源分配对跟踪所带来的性能提升问题，考虑将ASP-SQ-IMM算法与固定采样周期的SQ-IMM算法进行对比。首先通过对自适应变化的采样周期取平均值，得到以下两个场景下的平均采样间隔，并将其设置为固定采样间隔再次对目标进行跟踪，两个量测场景中的位置均方根误差(RMSE)曲线如图5，图6所示。

实施结果分析：ASP-SQ-IMM算法比固定采样周期的SQ-IMM算法有较高的跟踪性能，并且通过均分根误差统计表1中可知，ASP-SQ-IMM算法的平均RMSE较SQ-IMM算法小；可见 ASP-SQ-IMM算法具有较高的跟踪精度(在高噪声环境下也能够保持较好的特性)，且能够合理地分配系统的时间资源，相比于固定采样周期的SQ-IMM算法有明显优势。

表1位置平均均方根误差

Claims (1)

1.一种基于序贯滤波交互式多模型的自适应采样方法,该方法包括：

步骤1：输入交互估计

首先假设IMM交互式多模型系统中存在N个模型，运动模型集合表示为

式中，M_k表示t_k时刻的运动模型；

按下式计算各模型滤波器输入的状态估计和协方差P^0j(t_k-1)：

其中，表示在t_k-1时刻模型i的状态估计，表示第l个模型的误差协方差，表示模型混合概率，其计算方式如下：

式中，π_ij表示由模型i转移到模型j的概率；μⁱ(t_k-1)表示t_k-1时刻运动模型i的更新概率，C^j表示模型j的归一化常数；

步骤2：计算各模型位置状态预测和预测协方差

设置一组采样间隔集合其中包含了具有n个不同取值的元素，从大至小选取采样间隔，选择第l个采样周期T_l，代入下式

其中，和P^0j(t_k-1)分别为步骤1中获得的位置估计输入估计与协方差交互输入估计的初始值；F^j,G^j表示模型j的状态转移矩阵与噪声驱动矩阵，Q^j(t_k-1)表示上一时刻的过程噪声自相关矩阵；

步骤3：计算预测位置协方差；

以二维场景为例，预测位置协方差在极坐标下表示为

其中，量测矩阵C＝diag{[1 0 0],2}_2×6，极坐标系下的变换矩阵J_p中的各元素由对应位置偏导数构成

表示采样周期T_l为时对应的x方向上的位置状态，表示表示采样周期为T_l时对应的y方向上的位置状态；

组合预测协方差计算方法如下：

式中，将t_k时刻的模型j预测概率表示为可由各位置滤波器组的预测概率与伪量测预测概率均值获得，其计算方式如下：

式中，π_ij为概率转移矩阵，分别表示模型j的各位置量测滤波器与伪量测滤波器的模型预测概率，μ^p，i(t_k-1)，μ^ε，i(t_k-1)分别表示上一时刻模型i的各位置量测滤波器与伪量测滤波器的模型概率；

进而可以获得采样周期为T_l时所对应的距离的预测误差标准差σ_r(T_l)和方位角的预测误差标准差σ_b(T_l)为：

和P_pol(T_l)_2,2的下标(n,n),n＝1,2,3表示矩阵中行和列对应位置；

步骤4：计算协方差门限值P_th：

距离对应的门限预测误差标准差σ_r,th与方位角对应的门限预测误差标准差σ_b,th计算方式如下：

其中，u_0.5a是置信度为P_CL时对应的双侧分位数；为预测方位角对应的波束宽度，P_CL表示概率置信度，L_g表示距离波门宽度，计算方式如下：

式中，M为阵元总数；λ为发射信号的波长；d为各天线单元间距离；为预测方位角；

步骤5：当且仅当T_l满足如下关系

σ_r(T_l)≤σ_r,thσ_b(T_l)≤σ_b,th (13)

选取第l个采样间隔T_l，作为t_k时刻的采样间隔即T(t_k)＝T_l；

否则继续取出中的下一个元素进行测试，直至选取到最小的采样周期T_min；定义上一时刻为t_k-1，则t_k时刻点为t_k＝t_k-1+T(t_k)；

步骤6：利用步骤5中获取的采样周期T(t_k)，产生t_k时刻的量测位置量测值和伪量测

其中，Z_c(t_k)表示t_k时刻的量测转换值；

非线性变换函数x(t_k)表示x方向上的位置状态，表示x方向上的速度，y(t_k)表示y方向上的位置状态，表示y方向上的速度，量测转换误差目标在直角坐标系下的状态为X(t_k)；r_m(t_k),θ_m(t_k),分别为距离量测、方位角量测和多普勒量测；对应的量测噪声为方差分别为且和不相关，和的相关系数为ρ；λ为偏差补偿因子；

步骤7：基于预测值计算量测转换误差协方差矩阵R^j(t_k)；为简洁表述矩阵中的元素，约定t_k时刻所对应的参量，以下标k的形式表示

协方差矩阵R^j中的各元素的计算方法如下：

其中，极坐标系下的预测状态r_t,θ_t,通过步骤2中预测状态获得，参见式(22)；预测误差方差通过预测误差协方差矩阵的元素获得，参见式(23)。

式中，x_t,y_t,为中对应的目标预测位置与速度元素，矩阵J的形式如下

步骤8：各模型序贯滤波

将步骤6，7中获得t_k时刻的量测转换值Z_c(t_k)和量测误差协方差R^j(t_k)及步骤2中获得的t_k-1时刻的模型j的输入估计X^0j(t_k-1)，P^0j(t_k-1)，代入当前模型的序贯滤波器中；

步骤9：模型概率更新

模型j在t_k时刻的更新概率μ^j(t_k)表示为位置量测模型概率μ^p,j(t_k)和伪量测模型概率μ^ε,j(t_k)的均值

其中，上标j表示运动模型，p表示与位置量测相关的信息，ε表示与伪量测相关的信息；各滤波器的似然函数Λ^p,j,Λ^ε,j计算公式如下：

式中e^p,j为量测残差，S^p,j为残差的自相关矩阵；C^p,j和C^ε,j分别表示模型j的位置量测和伪量测归一化常数；

步骤10：状态估计融合

将步骤7中获得的t_k时刻的序贯滤波估计结果P^ε,j(t_k)与模型更新概率μ^j(t_k)进行融合

步骤11：将各模型的滤波结果P^p,j(t_k)代回到步骤1的(2)式中，计算下一时刻t_k+1的各模型滤波器输入的状态估计和协方差。