基于广义时变平滑层策略的多普勒雷达滤波方法及系统

1.本技术实施例涉及雷达数据处理和目标跟踪技术领域，具体而言，涉及一种基于广义时变平滑层策略的多普勒雷达滤波方法及系统。


背景技术：

2.雷达目标跟踪系统，是利用雷达回波中包含的目标距离、方位角、多普勒频偏等信息，对目标数量和位置、速度、加速度等运动状态进行持续有效的估计的系统，广泛应用于自动驾驶、空中交通管制、气象监测、安防、国防军事等领域；多普勒雷达可以观测目标运动的径向速度信息，进而在雷达目标状态估计环节增加有效信息维度，大幅提高目标跟踪性能。
3.多普勒雷达的目标跟踪算法或目标跟踪系统需要预先设定雷达观测模型和目标运动模型；但是由于多普勒雷达的雷达观测模型是高度非线性的，会使得常用的滤波器产生较大的非线性估计误差，甚至导致滤波发散，如卡尔曼滤波器；而目标运动模型具有不确定性，传统的传统贝叶斯滤波方法，包括卡尔曼滤波器等，均依赖于对目标运动模型的精确描述，若滤波器中预设的目标运动模型和真实模型偏差较大，例如目标发生强机动运动，会导致跟踪误差急剧增大，甚至滤波发散丢失目标。
4.平滑变结构滤波器(smooth variable structure filter，svsf)方法是近年来备受关注的解决目标状态转移模型的不确定性问题的方法，能够保证状态估计误差的有界性，同时具有较低的计算复杂度，但是标准的平滑变结构滤波方法需要满足：线性的观测模型和满秩的观测矩阵，不适用于多普勒雷达的非线性并且欠定的观测模型。
5.现有的多普雷雷达序贯平滑变结构滤波器(sequential smooth variable structure filter，ssvsf)方法虽然克服了非线性且欠定观测模型的问题，但是ssvsf方法中存在参数敏感性问题，即ssvsf方法中使用了固定的平滑层参数，与实际应用中存在时变的目标运动模型的不确定度不匹配，导致滤波器的性能对平滑层参数的选择十分敏感，从而使得目标跟踪精度不理想。


技术实现要素：

6.本技术实施例提供一种基于广义时变平滑层策略的多普勒雷达滤波方法及系统，旨在保持目标运动模型不确定条件下的鲁棒估计性能的前提下，解决参数敏感性问题以提高多普勒雷达的目标跟踪精度。
7.第一方面，本技术实施例提供一种基于广义时变平滑层策略的多普勒雷达序贯平滑变结构滤波方法，应用于多普勒雷达，所述方法包括以下步骤：
8.将所述多普勒雷达对目标的位置观测变量，转换为笛卡尔坐标系下的线性位置观测变量；
9.基于所述线性位置观测变量，执行面向多普勒雷达观测模型的广义时变平滑层策略，其中，所述广义时变平滑层策略用于实时地判决目标运动模型的模型不确定度，所述模
型不确定度表征预设模型相比于真实的所述目标运动模型的失配程度，并为不同的所述模型不确定度分配不同的切换增益项；
10.根据当前帧的先验估计，以及所述模型不确定度对应的切换增益项，计算当前帧的第一后验状态估计；
11.利用所述多普勒雷达在当前帧的径向速度观测变量构造得到的距离-径向速度积伪量测变量，对所述当前帧的第一后验状态估计进行更新，得到当前帧的第二后验状态估计，将所述当前帧的第二后验状态估计作为当前帧最终的目标状态估计。
12.可选地，所述广义时变平滑层策略包括：
13.根据当前帧的先验估计以及所述线性位置观测变量，计算所述多普勒雷达观测模型的广义时变平滑层参数；
14.若所述广义时变平滑层参数大于预设的平滑层参数，判决所述目标运动模型在当前帧为高模型不确定度；
15.若所述广义时变平滑层参数小于等于预设的平滑层参数，判决所述目标运动模型在当前帧为低模型不确定度。
16.可选地，根据当前帧的先验估计以及所述线性位置观测变量，计算所述多普勒雷达观测模型的广义时变平滑层参数，包括：
17.计算所述当前帧的先验估计，所述当前帧的先验估计包括目标状态先验估计、目标位置量测的先验估计和先验状态估计协方差阵；
18.根据所述线性位置观测变量、第一转换偏差以及所述目标位置量测的先验估计，计算先验位置观测误差；
19.根据所述先验位置观测误差，计算混合误差项；
20.根据所述先验状态估计协方差阵、所述混合误差项以及所述第一转换误差协方差，计算得到所述当前帧的广义时变平滑层参数。
21.可选地，所述广义时变平滑层策略还包括：
22.所述目标运动模型在当前帧为高模型不确定度时，根据所述预设的平滑层参数，计算标准切换增益项；
23.所述目标运动模型在当前帧为低模型不确定度时，根据所述广义时变平滑层参数，计算最优切换增益项。
24.可选地，根据当前帧的先验估计，以及所述模型不确定度对应的切换增益项，计算当前帧的第一后验状态估计，包括：
25.所述目标运动模型在当前帧为高模型不确定度时，根据所述当前帧的先验估计以及所述标准切换增益项，计算所述当前帧的第一后验状态估计；
26.所述目标运动模型在当前帧为低模型不确定度时，根据所述当前帧的先验估计以及所述最优切换增益项，计算所述当前帧的第一后验状态估计。
27.可选地，计算当前帧的第一后验状态估计，包括：
28.根据所述目标状态先验估计、所述先验状态估计协方差阵、所述先验位置观测误差，以及所述目标运动模型的模型不确定度对应的切换增益项，计算目标状态后验估计与后验状态估计协方差阵；
29.将所述目标状态后验估计与所述后验状态估计协方差阵作为所述当前帧的第一
后验状态估计。
30.可选地，将所述多普勒雷达对目标的位置观测变量，转换为笛卡尔坐标系下的线性位置观测变量，包括：
31.采用修正的无偏量测转换方法，将所述多普勒雷达在三维球坐标或者二维极坐标下的所述位置观测变量，转换为笛卡尔坐标系下的所述线性位置观测变量，并计算第一转换偏差和第一转换误差协方差。
32.可选地，利用所述多普勒雷达在当前帧的径向速度观测变量构造得到的距离-径向速度积伪量测变量，对所述当前帧的第一后验状态估计进行更新，得到当前帧的第二后验状态估计，将所述当前帧的第二后验状态估计作为当前帧最终的目标状态估计，包括：
33.根据所述多普勒雷达在当前帧的距离观测变量以及径向速度观测变量，构造距离-径向速度积伪量测变量，并采用修正的无偏转换方法计算所述距离-径向速度积伪量测变量的第二转换偏差、转换误差方差以及第二转换误差协方差；
34.对所述距离-径向速度积伪量测变量进行去相关处理；
35.将去相关处理后的所述距离-径向速度积伪量测变量作为输入，构造局部近似的线性最小均方误差估计器，对所述当前帧的第一后验状态估计进行更新，将更新后得到的所述当前帧的第二后验状态估计作为当前帧最终的目标状态估计。
36.可选地，所述广义时变平滑层参数的计算公式为：
[0037][0038][0039][0040]
式中，h
p,1
＝i3×3是观测矩阵h
p
(k)的满秩分块子矩阵；是主对角元素为矢量a的对角阵；对预设的状态转移矩阵引入变换矩阵t后的变换状态转移矩阵，s
p
(k)是所述位置观测变量的新息协方差阵；ez和ey为所述混合误差项；和为所述先验状态估计协方差阵的分块子矩阵。
[0041]
第二方面，本技术实施例提供一种基于广义时变平滑层策略的多普勒雷达序贯平滑变结构滤波系统，所述系统包括：
[0042]
线性转换模块，用于将所述多普勒雷达对目标的位置观测变量，转换为笛卡尔坐标系下的线性位置观测变量；
[0043]
策略执行模块，用于基于所述线性位置观测变量，执行面向多普勒雷达观测模型的广义时变平滑层策略，其中，所述广义时变平滑层策略用于实时地判决目标运动模型的模型不确定度，所述模型不确定度表征预设模型相比于真实的所述目标运动模型的失配程度，并为不同的所述模型不确定度分配不同的切换增益项；
[0044]
第一状态估计模块，用于根据当前帧的先验估计，以及所述模型不确定度对应的切换增益项，计算当前帧的第一后验状态估计；
[0045]
第二状态估计模块，用于利用所述多普勒雷达在当前帧的径向速度观测变量构造得到的距离-径向速度积伪量测变量，对所述当前帧的第一后验状态估计进行更新，得到当前帧的第二后验状态估计，将所述当前帧的第二后验状态估计作为当前帧最终的目标状态估计。
[0046]
可选地，所述策略执行模块包括：
[0047]
参数计算单元，用于根据当前帧的先验估计以及所述线性位置观测变量，计算所述多普勒雷达观测模型的广义时变平滑层参数；
[0048]
模型不确定度判决单元，用于在所述广义时变平滑层参数大于预设的平滑层参数时，判决所述目标运动模型在当前帧为高模型不确定度；在所述广义时变平滑层参数小于等于预设的平滑层参数时，判决所述目标运动模型在当前帧为低模型不确定度。
[0049]
可选地，所述参数计算单元包括：
[0050]
第一计算子单元，用于计算所述当前帧的先验估计，所述当前帧的先验估计包括目标状态先验估计、目标位置量测的先验估计和先验状态估计协方差阵；
[0051]
第二计算子单元，用于根据所述线性位置观测变量、第一转换偏差以及所述目标位置量测的先验估计，计算先验位置观测误差；
[0052]
第三计算子单元，用于根据所述先验位置观测误差，计算混合误差项；
[0053]
第四计算子单元，用于根据所述先验状态估计协方差阵、所述混合误差项、所述第一转换误差协方差，计算所述当前帧的广义时变平滑层参数。
[0054]
可选地，所述策略执行模块还包括：
[0055]
标准切换增益项计算单元，用于在所述目标运动模型在当前帧为高模型不确定度时，根据所述预设的平滑层参数，计算标准切换增益项；
[0056]
最优切换增益项计算单元，用于在所述目标运动模型在当前帧为低模型不确定度时，根据所述广义时变平滑层参数，计算最优切换增益项。
[0057]
可选地，所述第一状态估计模块包括：
[0058]
第一估计单元，用于在所述目标运动模型在当前帧为高模型不确定度时，根据所述当前帧的先验估计以及所述标准切换增益项，计算所述当前帧的第一后验状态估计；
[0059]
第二估计单元，用于在所述目标运动模型在当前帧为低模型不确定度时，根据所述当前帧的先验估计以及所述最优切换增益项，计算所述当前帧的第一后验状态估计。
[0060]
可选地，所述第一状态估计模块包括：
[0061]
状态估计计算单元，用于根据所述目标状态先验估计、所述先验状态估计协方差阵、所述先验位置观测误差，以及所述目标运动模型的模型不确定度对应的切换增益项，计算目标状态后验估计与后验状态估计协方差阵；将所述目标状态后验估计与所述后验状态估计协方差阵作为所述当前帧的第一后验状态估计。
[0062]
可选地，所述线性转换模块包括：
[0063]
第一转换单元，用于采用修正的无偏量测转换方法，将所述多普勒雷达在三维球坐标或者二维极坐标下的所述位置观测变量，转换为笛卡尔坐标系下的所述线性位置观测变量，并计算第一转换偏差和第一转换误差协方差。
[0064]
可选地，所述第二状态估计模块包括：
[0065]
第二转换单元，用于根据所述多普勒雷达在当前帧的距离观测变量以及径向速度观测变量，构造距离-径向速度积伪量测变量，并采用修正的无偏转换方法计算所述距离-径向速度积伪量测变量的第二转换偏差、转换误差方差以及第二转换误差协方差；
[0066]
去相关处理单元，用于对所述距离-径向速度积伪量测变量进行去相关处理；
[0067]
最终估计单元，用于将去相关处理后的所述距离-径向速度积伪量测变量作为输入，构造局部近似的线性最小均方误差估计器，对所述当前帧的第一后验状态估计进行更新，将更新后得到的所述当前帧的第二后验状态估计作为当前帧最终的目标状态估计。
[0068]
有益效果：
[0069]
在进行目标状态估计时，执行面向多普勒雷达观测模型的广义时变平滑层策略，广义时变平滑层策略可以实时地判决目标运动模型的模型不确定度，模型不确定度表征预设模型相比于真实的所述目标运动模型的失配程度，并且对于不同的模型不确定度对应分配有不同的切换增益项，根据当前帧的先验估计以及所述模型不确定度对应的切换增益项，计算当前帧的第一后验状态估计后，再利用距离-径向速度积伪量测变量对第一后验状态估计进行更新，可以得到更加准确的第二后验状态估计，并将第二后验状态估计作为当前帧最终的目标状态估计。
[0070]
本方法考虑到真实的目标运动模型与预设模型之间的失配误差具有时变性特点，进而会导致使用固定平滑层参数的切换增益项不准确，从而目标跟踪精度差；因此，为了提高目标跟踪的精度，执行广义时变平滑层策略实时判决目标运动模型的不确定度，在高模型不确定度条件下采用固定平滑层参数的标准切换增益项，在低模型不确定度条件下采用基于广义时变平滑层参数的最优切换增益项；相比于传统的多普勒雷达svsf方法，本方法在保持目标运动模型不确定条件下的鲁棒估计性能的前提下，解决了参数敏感性问题，从而提高了多普勒雷达的目标跟踪精度。
附图说明
[0071]
为了更清楚地说明本技术实施例的技术方案，下面将对本技术实施例的描述中所需要使用的附图作简单地介绍。显而易见地，下面描述中的附图仅仅是本技术的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
[0072]
图1是本技术一个实施例提出的基于广义时变平滑层策略的多普勒雷达序贯平滑变结构滤波方法的步骤流程图；
[0073]
图2是本技术一个实施例提出的广义时变平滑层策略的步骤流程图；
[0074]
图3是本技术一个实施例提出的广义时变平滑层参数的计算方法的步骤流程图；
[0075]
图4是本技术一个实施例提出的基于广义时变平滑层策略的多普勒雷达序贯平滑变结构滤波系统的功能模块图。
具体实施方式
[0076]
下面将结合本技术实施例中的附图，对本技术实施例中的技术方案进行清楚、完整地描述。显然，所描述的实施例是本技术一部分实施例，而不是全部的实施例。基于本申
请中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本技术保护的范围。
[0077]
平滑变结构滤波器(smooth variable structure filter，svsf)可以解决目标运动模型的不确定性问题，可以保证状态估计误差的有界性，但是标准的平滑变结构滤波方法要求线性的观测模型和满秩的观测矩阵；因此为了利用平滑变结构滤波理论的优势，在现有的序贯平滑变结构滤波(sequential smooth variable structure filter，ssvsf)方法中，将多普勒雷达的高度非线性的雷达观测模型生成的位置观测变量进行线性转换，从而解决了非线性且欠定的问题。
[0078]
但是，现有的序贯平滑变结构滤波(sequential smooth variable structure filter，ssvsf)方法中使用的平滑层参数是一个固定且预设的经验值，当参数设定不准确或者面临时变的模型不确定度场景时，会造成较大的估计误差，因此采用预设的平滑层参数的ssvsf方法存在参数敏感性问题，会造成状态估计精度和目标跟踪性能不理想。
[0079]
示例地，如果预设的平滑层参数过小会加剧抖振问题；如果预设的平滑层参数太大，则会损失滤波器对建模误差的鲁棒性甚至导致滤波发散；特别在模型不确定度较小时，采用预设的平滑层参数的ssvsf不具有最优性能，进而导致滤波精度和目标跟踪性能差于传统的贝叶斯滤波器，例如扩展卡尔曼滤波器(expanded kalmanfilter，ekf)。
[0080]
因此，本技术实施例提出了一种基于广义时变平滑层策略的多普勒雷达序贯平滑变结构滤波方法，可以在保持目标运动模型不确定条件下的鲁棒估计性能的同时，有效解决参数敏感性问题，提高多普勒雷达的目标跟踪精度。
[0081]
参照图1，示出了本发明实施例中的一种基于广义时变平滑层策略的多普勒雷达序贯平滑变结构滤波方法的步骤流程图，所述方法应用于多普勒雷达，所述方法具体可以包括以下步骤：
[0082]
定义应用于多普勒雷达目标跟踪的目标运动模型和观测模型分别为：
[0083]
x(k)＝f(k-1)x(k-1)+g(k-1)u(k-1)+γ(k-1)w(k-1)
[0084]
zm(k)＝h(x(k))+v(k)
[0085]
其中，为描述目标的位置、速度等信息的状态变量；是描述目标运动的真实状态转移矩阵，u(k)和g(k)是已知输入项和其输入系数矩阵；w(k)是状态噪声，其系数矩阵为γ(k)，噪声协方差为q(k)。
[0086]
是多普勒雷达的观测变量，其中rm为距离观测变量，θm为方位角观测变量，为俯仰角观测变量，为径向速度观测变量；v(k)是对应各个观测维度的观测噪声，其噪声协方差阵为r(k)。
[0087]
h(
·
)是多普勒雷达的非线性欠定的观测函数，定义为：
[0088][0089]
注意，考虑目标运动模型的不确定性问题，在对目标进行跟踪时使用预设的状态
转移矩阵和输入系数矩阵相比于理想的目标运动模型存在模型误差和
[0090]
s101：将所述多普勒雷达对目标的位置观测变量，转换为笛卡尔坐标系下的线性位置观测变量。
[0091]
根据多普勒雷达的观测变量的结构，将多普勒雷达对目标在当前帧的观测变量分割为位置观测变量和径向速度观测变量，表现形式如下式：
[0092][0093]
式中，是包含有距离观测变量、方位角观测变量和俯仰角观测变量的位置观测变量，是径向速度观测变量，k表示多普勒雷达的帧数，k＝1，2，3
……
。
[0094]
在一种可行的实施方式中，在对非线性的位置观测变量进行线性转换时，采用修正的无偏量测转换方法(modified unbiased converted measurement，mucm)，将所述多普勒雷达在球坐标系或极坐标系下的位置观测变量转换为笛卡尔坐标系下的线性位置观测变量并计算第一转换偏差μ
p
以及第一转换误差协方差r
p
(k)。
[0095]
需要说明的是，本实施例以三维多普勒雷达为例，在其他实施例中，例如对于二维场景可以直接截取有效维度。
[0096]
s102：基于所述线性位置观测变量，执行面向多普勒雷达观测模型的广义时变平滑层策略。
[0097]
所述广义时变平滑层策略用于实时地判决目标运动模型的模型不确定度，所述模型不确定度表征预设模型相比于真实的所述目标运动模型的失配程度，并为不同的所述模型不确定度分配不同的切换增益项。
[0098]
参照图2，示出了本技术实施例提供的广义时变平滑层策略的步骤流程图，在一种可行的实施方式中，所述广义时变平滑层策略包括：
[0099]
a1：根据当前帧的先验估计以及所述线性位置观测变量，计算所述多普勒雷达观测模型的广义时变平滑层参数。
[0100]
参照图3，示出了本技术实施例提供的广义时变平滑层参数的计算方法的步骤流程图，具体地，可以通过以下步骤确定当前帧的广义时变平滑层参数：
[0101]
b1：计算所述当前帧的先验估计，所述当前帧的先验估计包括目标状态先验估计、目标位置量测的先验估计和先验状态估计协方差阵。
[0102]
具体地，目标状态先验估计目标位置量测的先验估计和先验状态估计协方差阵p(k|k-1)，可以根据以下公式确定：
[0103][0104]
[0105][0106]
式中，为状态估计变量，是观测矩阵，n和m分别表示状态变量和观测变量的维度。
[0107]
特别地，和是不同于传统贝叶斯滤波器要求准确已知的状态转移模型的f和g，和相比于真实值允许存在范数有界的模型误差，因此在模型不确定条件下具有鲁棒性优势；同时，由于本方法中应用修正的无偏量测转换方法(mucm)将非线性的位置观测变量转换为线性位置观测变量，进而观测矩阵h
p
(k)＝[i3×3,03×
(n-3)
]是线性时不变的，优于直接使用非线性的位置观测变量的平滑变结构滤波器。
[0108]
b2：根据所述线性位置观测变量、第一转换偏差以及所述目标位置量测的先验估计，计算先验位置观测误差。
[0109]
具体地，先验位置观测误差e
p
(k|k-1)可以根据以下公式确定：
[0110][0111]
其中，为第k帧的线性位置观测变量，μ
p
(k)为第k帧的第一转换偏差，为第k帧的目标位置量测的先验估计。
[0112]
b3：根据所述先验位置观测误差，计算混合误差项。
[0113]
具体地，根据先验位置观测误差e
p
(k|k-1)，混合误差项ez和ey的计算公式如下：
[0114]ez
(k)＝|e
p
(k|k-1)
abs
+γz|e
p
(k-1|k-1)
abs
[0115]ey
(k)＝|φ
22
φ
12-1ep
(k|k-1)|
abs
+γy|φ
12-1ep
(k-1|k-1)|
abs
.
[0116]
其中，γz和γy是取值在(0,1)区间的衰减因子；e
p
(k-1|k-1)为k-1帧的后验位置观测变量误差，|
·
|
abs
表示对矢量逐元素取绝对值；是对预设的状态转移矩阵引入变换矩阵t后的变换状态转移矩阵，
[0117]
b4：根据所述先验状态估计协方差阵、所述混合误差项以及所述第一转换误差协方差，计算得到所述当前帧的广义时变平滑层参数。
[0118]
具体地，根据先验状态估计协方差阵p(k|k-1)、混合误差项ez和ey以及第一转换误差协方差r
p
(k)，计算当前帧的广义时变平滑层参数ψ
gvbl
，可以按照如下公式计算：
[0119][0120][0121][0122]
式中，h
p,1
＝i3×3是观测矩阵h
p
(k)的满秩分块子矩阵；是主对角元素为矢量a的对
角阵；s
p
(k)是位置观测变量的新息协方差阵，计算方式如下：
[0123]sp
(k)＝h
p
(k)p(k|k-1)h
p
(k)
t
+r
p
(k)。
[0124]
参见前述，先验状态估计协方差阵p(k|k-1)可以表示为：
[0125][0126]
为先验状态估计协方差阵p(k|k-1)左上角的分块子矩阵，为先验状态估计协方差阵p(k|k-1)左下角的分块子矩阵。
[0127]
本实施方式中，基于最小均方误差准则推导广义时变平滑层参数，并且广义时变平滑层参数是根据当前帧的先验估计以及当前帧的线性位置观测变量确定的，因此广义时变平滑层参数包含了当前时刻的模型不确定性信息。
[0128]
a2：比较所述广义时变平滑层参数与预设的平滑层参数，确定所述目标运动模型在当前帧的模型不确定度。
[0129]
预设的平滑层参数即为ssvsf方法中用于计算增益项的固定参数ψ
max
，为了解决ssvsf方法对预设的平滑层参数ψ
max
的敏感性问题，对目标运动模型的不确定度进行实时判决，对不同的不确定度采用不同的平滑层参数计算切换增益项。
[0130]
具体地，将实时计算得到的广义时变平滑层参数ψ
gvbl
与预设的平滑层参数ψ
max
按各个状态维度进行比较，进而确定目标运动模型的不确定度。
[0131]
在所述广义时变平滑层参数ψ
gvbl
大于预设的平滑层参数ψ
max
时，判决所述目标运动模型在当前帧为高模型不确定度。
[0132]
在所述广义时变平滑层参数ψ
gvbl
小于等于预设的平滑层参数ψ
max
时，判决所述目标运动模型在当前帧为低模型不确定度。
[0133]
a3：确定所述目标运动模型在当前帧的模型不确定度对应的切换增益项。
[0134]
具体地，所述目标运动模型在当前帧为高模型不确定度时，根据所述预设的平滑层参数ψ
max
，计算标准切换增益项；所述目标运动模型在当前帧为低模型不确定度时，根据所述广义时变平滑层参数ψ
gvbl
，计算最优切换增益项。
[0135]
示例地，在所述目标运动模型为高模型不确定度时，根据所述预设的平滑层参数ψ
max
，计算所述标准切换增益项，计算过程如下式：
[0136][0137][0138][0139]
式中，sat(.)表示饱和函数；此时的和是预设的平滑层参数矢量。
[0140]
在所述目标运动模型为低模型不确定度时，根据广义时变平滑层参数ψ
gvbl
，计算所述最优切换增益项即将ψ
gvbl
代替ψ
max
带入标准切换增益项的计算过程，也就
是说，将预设的平滑层参数矢量ψz与ψy用步骤b4中计算得到的与代替，经过计算化简后，所得到的化简后的可以表示为下式：
[0141][0142]
也就是说，在高模型不确定度下，采用基于固定平滑层参数的标准切换增益项，这是一种基于lyapunov稳定性理论和luenberger降阶观测器原理设计的增益项，能够保证模型不确定条件下的鲁棒状态估计性能，即在高模型不确定度的情况下，优先选择保证鲁棒估计性能。
[0143]
在低模型不确定度时，则采用基于广义时变平滑层参数的最优切换增益项进行状态估计，可以在模型不确定的程度低的情况下实现对目标状态的近似最优估计，避免参数敏感性问题。
[0144]
广义时变平滑层参数ψ
gvbl
和基于广义平滑层参数的最优切换增益是在序贯平滑变结构滤波器(ssvsf)的基础上，采用最小均方误差准则推导的最优解，推导过程如下：
[0145]
给定最小均方误差准则下的优化目标函数为：
[0146][0147]
其中p(k|k)是序贯平滑变结构滤波器的最终后验状态估计误差协方差，ψ是平滑层参数，k
ε
(k)是ssvsf方法的第二级状态估计器的增益项，tr表示矩阵的迹。
[0148]
求解该优化问题需分别考虑如下两个矩阵偏导方程：
[0149][0150][0151]
其中，第一个方程可以由线性最小均方误差估计理论直接给出最优解；第二方程方程需考虑两级序贯估计器的协方差传递过程，具体的求解过程如下：
[0152]
定义：
[0153][0154]
其中，h
ε
(k)是序贯平滑变结构滤波器的第二级状态估计器中的线性化观测矩阵，r
ε
(k)为线性化观测误差方差。
[0155]
由此，可以将优化目标矩阵p(k|k)写为如下的分块矩阵形式：
[0156][0157]
则，优化目标矩阵p(k|k)的迹可以进一步获得为：
[0158]
tr(p(k|k))＝tr(θ
11
(k))+tr(θ
22
(k))+tr(k
ε
(k)r
ε
(k)k
ε
(k)
t
),
[0159]
其中，θ(k)的分块子矩阵分别为：
[0160]
θ
11
(k)＝a
11
(k)p
p,11
(k|k)a
11
(k)
t
+a
12
(k)p
p,21
(k|k)a
11
(k)
t
+a
11
(k)p
p,12
(k|k)a
12
(k)
t
+a
12
(k)p
p,22
(k|k)a
12
(k)
t
[0161]
θ
22
(k)＝a
21
(k)p
p,11
(k|k)a
21
(k)
t
+a
22
(k)p
p,21
(k|k)a
21
(k)
t
+a
21
(k)p
p,12
(k|k)a
22
(k)
t
+a
22
(k)p
p,22
(k|k)a
22
(k)
t
.
[0162]
p
p
(k|k)的分块子矩阵分别为：
[0163]
p
p,11
(k|k)＝p
11
(k|k-1)-ku(k)h
p,1
p
11
(k|k-1)-p
11
(k|k-1)h
p,1tku
(k)
t
+ku(k)s
p
(k)ku(k)
t
[0164]
p
p,12
(k|k)＝p
12
(k|k-1)-ku(k)h
p,1
p
12
(k|k-1)-p
11
(k|k-1)h
p,1tkl
(k)
t
+ku(k)s
p
(k)k
l
(k)
t
[0165]
p
p,21
(k|k)＝p
21
(k|k-1)-k
l
(k)h
p,1
p
11
(k|k-1)-p
21
(k|k-1)h
p,1tku
(k)
t
+k
l
(k)s
p
(k)ku(k)
t
[0166]
p
p,22
(k|k)＝p
22
(k|k-1)-k
l
(k)h
p,1
p
12
(k|k-1)-p
21
(k|k-1)h
p,1tkl
(k)
t
+k
l
(k)s
p
(k)k
l
(k)
t
[0167]
在低模型不确定度的判决假设下，ku(k)和k
l
(k)可以分别简化为：
[0168][0169][0170]
注意到tr(p(k|k))的表达式中第三项与优化变量ψ无关，因此，优化目标函数可以进一步写为：
[0171][0172]
更具体地，可以拆分写为：
[0173][0174][0175]
注意到两个偏导函数均为偏导变量的二次型凸函数，将ku(k)和k
l
(k)、p
p
(k|k)的分块子矩阵代入θ
11
(k)和θ
22
(k)，利用矩阵偏导计算方法可得最优广义平滑层参数：
[0176][0177]
同时，可解得最优切换增益项的子块矩阵：
[0178][0179][0180]
注意到k
p
(k)的原始定义结构，最优切换增益项可写为：
[0181][0182]
至此，完成了低模型不确定度的判决假设下，广义时变平滑层参数及最优切换增益项的推导。
[0183]
根据上述推导过程可知，本方法基于最小均方误差准则，通过求解序贯平滑变结构滤波器的后验状态估计误差协方差的偏导方程，推导得到广义时变平滑层参数。
[0184]
在低模型不确定度的情况下，广义时变平滑层参数能够保证在多普勒雷达的一阶近似的非线性伪量测模型上，取得状态估计误差的最优性，因此，相比于传统的采用固定平滑层参数的平滑变结构滤波器，本方法能够避免不合理的预设参数值引起的性能恶化，特别是在低模型不确定度的情况下实现对目标状态的近似最优估计，从而提高多普勒雷达的目标跟踪精度。
[0185]
s103：根据当前帧的先验估计，以及所述模型不确定度对应的切换增益项，计算当前帧的第一后验状态估计。
[0186]
具体地，根据目标状态先验估计先验状态估计协方差阵p(k|k-1)、先验位置观测误差e
p
(k|k-1)以及所述目标的模型在所述当前帧的模型不确定度对应的切换增益项，计算目标状态后验估计与后验状态估计协方差阵p
p
(k|k)，计算公式如下：
[0187][0188]
p
p
(k|k)＝(i-k
p
(k)h
p
(k))p(k|k-1)(i-k
p
(k)h
p
(k))
t
+k
p
(k)r
p
(k)k
p
(k)
t
[0189]
其中，在低模型不确定度下，
[0190]
同时，根据当前帧的线性位置观测变量当前帧的第一转换偏差μ
p
(k)，以及目标状态后验估计计算后验位置观测变量误差e
p
(k|k)，计算公式如下：
[0191][0192]
式中，h
p
(k)为观测矩阵。
[0193]
将计算得到目标状态后验估计与后验状态估计协方差阵p
p
(k|k)作为当前帧的第一后验状态估计；后验位置观测变量误差e
p
(k|k)用于处理下一帧的观测变量。
[0194]
s104：利用所述多普勒雷达在当前帧的径向速度观测变量构造得到的距离-径向速度积伪量测变量，对所述当前帧的第一后验状态估计进行更新，得到当前帧的第二后验状态估计，将所述当前帧的第二后验状态估计作为当前帧最终的目标状态估计。
[0195]
在本步骤中，采用一个一阶近似的线性最小均方误差准则估计器和距离-径向速度积伪量测变量，对当前帧的第一后验状态估计进行更新。
[0196]
首先，根据多普勒雷达在当前帧的距离观测变量rm与径向速度观测变量构造得到的距离-径向速度积伪量测变量，距离-径向速度积伪量测变量可以定义为：
[0197]
[0198]
采用带有径向速度观测变量的修正的无偏量测转换方法(mucm)，计距离-径向速度积伪量测变量的第二转换偏差μ
η
、转换误差方差r
η
以及第二转换误差协方差r
ηp
，计算公式如下：
[0199][0200][0201][0202]
其中，ρ表示距离观测和径向速度观测误差的相关系数，σr表示距离观测变量的观测误差标准差，表示径向速度观测变量的观测误差标准差，表示俯仰角观测变量的观测误差标准差，σ
θ
表示方位角观测变量的观测误差标准差；在实际实施的过程中，根据多普勒雷达的型号种类，可以查询到每个观测变量的误差标准差以及相关系数。
[0203]
接着，对所述距离-径向速度积伪量测变量进行去相关处理，去相关系数矩阵为：
[0204]
l(k)＝-r
ηp
(k)(r
p
(k))-1
＝[l1(k),l2(k),l3(k)]
[0205]
去相关的伪量测变量的计算公式为：
[0206][0207]
去相关的伪量测的第二转换偏差的计算公式为：
[0208]
μ
ε
(k)＝l(k)μ
p
(k)+μ
η
(k)
[0209]
去相关的伪量测的方差的计算公式为：
[0210]rε
(k)＝r
η
(k)-r
ηp
(k)(r
p
(k))-1
(r
ηp
(k))
t
[0211]
去相关伪量测变量的观测模型可以表示为：
[0212][0213]
式中，h
ε
(
·
)与分别为去相关伪量测变量的非线性观测模型和观测误差。
[0214]
将第一后验状态估计和去相关的伪量测变量作为输入，构造如下局部近似的线性最小均方误差估计器，对当前帧的第一后验状态估计进行更新，得到当前帧的第二后验状态估计：
[0215]kε
(k)＝p
p
(k|k)h
ε
(k)
t
(h
ε
(k)p
p
(k|k)h
ε
(k)
t
+r
ε
(k))-1
[0216][0217][0218]
p(k|k)＝(i-k
ε
(k)h
ε
(k))p
p
(k|k)(i-k
ε
(k)h
ε
(k))
t
+k
ε
(k)r
ε
(k)k
ε
(k)
t
[0219]
其中，非线性函数的定义如下：
[0220]
[0221]
矩阵h
ε
(k)定义为非线性函数(k)定义为非线性函数(k)定义为非线性函数
[0222]
至此，对于多普勒雷达的当前帧的数据处理完毕，得到的当前帧的第二后验状态估计包括最终的目标状态后验估计与最终的后验状态估计协方差阵p(k|k)，当前帧的第二后验状态估计即为当前帧的最终的后验状态估计的结果。
[0223]
最终的目标状态后验估计是一个矢量，其包含各个坐标系维度的位置、速度等状态估计的标量。
[0224]
对k＝1，2，3
……
帧的多普勒雷达的回波数据重复步骤s101-s104的序贯迭代过程，可以获得在模型不确定条件下的非线性鲁棒估计，并且根据不同的模型不确定度自适应地切换增益项，可以进一步提高多普勒雷达的目标跟踪精度。
[0225]
本方法实质涉及了一种基于广义时变平滑层策略的序贯平滑变结构滤波器，整体包括两级序贯处理过程，在第一级序贯处理过程中，将非线性的位置观测变量转换为线性位置观测变量，并实时计算一个广义时变平滑层参数，利用广义时变平滑层参数与预设的平滑层参数比较来判决当前目标运动模型的模型不确定度，根据判决结果自适应地选择切换增益项的计算方法，若为高模型不确定度，则采用标准的广义平滑变结构滤波器的标准切换增益项来保证鲁棒估计性能；反之，若为低模型不确定度，则采用基于广义时变平滑层参数计算的最优切换增益项来获取近似最优性能，第一级序贯处理过程输出第一后验状态估计。
[0226]
在第二级处理过程中，根据去相关的距离-径向速度积伪量测变量与一个一阶近似的线性最小均方误差准则估计器更新第一级序贯处理过程输出的第一后验状态估计，从而获得当前帧的最终的后验状态估计。
[0227]
本技术至少具有以下效果：
[0228]
1、通过计算一个广义时变平滑层参数，实时地判决当前帧的目标运动模型的不确定度，进而自适应地选择切换增益项；
[0229]
2、在高模型不确定度时，采用标准切换增益项进行状态估计，能够保持多普勒雷达序贯平滑变结构滤波方法的鲁棒估计优势；
[0230]
3、在低模型不确定度时，采用基于广义时变平滑层参数的最优切换增益项进行状态估计，能够获得近似最优估计性能，避免参数敏感性问题，因此能够有效提高目标跟踪精度；
[0231]
4、本方法能够在保持模型不确定条件下的鲁棒估计性能的前提下，通过参数优化方法解决现有的多普勒雷达平滑变结构滤波方法的参数敏感性问题，从而提高多普勒雷达系统的目标跟踪性能与精度；
[0232]
5、本方法融合了传统svsf方法在参数最优性和鲁棒性问题上的收益，通过执行广义时变平滑层策略实时地评估目标运动模型的模型不确定度，互补地利用了基于广义时变平滑层参数的最优切换增益和基于固定参数的标准切换增益两者分别在低不确定度和高不确定度情况下的优势，因此取得了比传统方法更好的目标跟踪精度。
[0233]
本实施例还提供一个利用公开的车载雷达数据集radarscenes中的两个场景，对本实施例提供的基于广义时变平滑层的多普勒雷达序贯平滑变结构滤波方法、现有的扩展卡尔曼滤波(ekf)方法、仅利用位置量测的平滑变结构滤波(position-only smooth variable structure filter，po-svsf)方法、采用雅克比矩阵局部线性化的平滑变结构滤波(smooth variable structure filter with local linearization，l-svsf)方法，以及序贯平滑变结构滤波(ssvsf)方法进行对比。
[0234]
场景1包含一个机动运动目标，是一种高模型不确定度的环境；场景2包含郊区公路上大量近似匀速直线运动的目标，是一种低模型不确定度的环境。测试中的评估指标包括：用方均根误差(rmse)对比目标跟踪精度，用航迹中断次数对比航迹连续度。
[0235]
目标的状态变量定义为x-y坐标轴内的位置、速度，即x＝[x y v
x vy]
t
；雷达观测变量定义采用匀速运动模型(cv)对目标状态进行估计，即：
[0236][0237]
第一级状态估计器的量测矩阵为：
[0238][0239]
第二级状态估计器的线性化量测矩阵由说明书给出方法实时计算。
[0240]
不同的滤波方法在两个场景中的跟踪性能评估如表1所示。
[0241]
表1基于车载雷达实测数据的目标跟踪性能评估
[0242][0243]
由表1可知，本实施例提出的基于广义时变平滑层策略的多普勒雷达序贯平滑变
结构滤波方法在两种场景中，位置误差与速度误差均是最小的，即具有最好的目标跟踪精度；航迹中断次数仅是略多于序贯平滑变结构滤波方法(ssvsf)，表现仍优异。
[0244]
本实施例提供的滤波方法在近似匀速直线运动的目标跟踪场景(即在低模型不确定度)中，表现是几种方法中最优异的，这是因为本方法在低模型不确定度时，采用基于广义时变平滑层参数的最优切换增益项进行状态估计，能够获得近似最优估计性能，避免参数敏感性问题，从而提高目标跟踪精度。
[0245]
本实施例提供的滤波方法在目标机动运动时(即在高模型不确定度下)仍大幅提升了跟踪性能，这是由于本方法保持了高模型不确定条件下的鲁棒性优势，相比于只利用位置量测的po-svsf方法，本方法可以有效利用多普勒速度信息，大幅提高了目标状态估计精度。
[0246]
ssvsf方法在低模型不确定度下的目标跟踪精度表现差于传统的ekf，这是因为ssvsf方法采用了固定平滑层参数，不能随实际环境中模型误差的强度而自适应调整，因而性能较为保守。
[0247]
相比之下，本实施例提供的滤波方法由于采用自适应切换策略计算切换增益项，兼顾了鲁棒性和参数最优性，尤其是在低模型不确定度下具有显著的优势。
[0248]
综上所述，本实施例提供的滤波方法适用于多普勒雷达的非线性观测场景，能够获得低模型不确定度下的目标状态的近似最优估计，同时保证高模型不确定度下的鲁棒性，实现更加精确和稳定的雷达目标跟踪性能，具有很好的应用价值。
[0249]
参照图4，示出了本发明实施例中的一种基于广义时变平滑层策略的多普勒雷达序贯平滑变结构滤波系统的功能模块图，所述系统包括：
[0250]
线性转换模块100，用于将所述多普勒雷达对目标的位置观测变量，转换为笛卡尔坐标系下的线性位置观测变量；
[0251]
策略执行模块200，用于基于所述线性位置观测变量，执行面向多普勒雷达观测模型的广义时变平滑层策略，其中，所述广义时变平滑层策略用于实时地判决目标运动模型的模型不确定度，所述模型不确定度表征预设模型相比于真实的所述目标运动模型的失配程度，并为不同的所述模型不确定度分配不同的切换增益项；
[0252]
第一状态估计模块300，用于根据当前帧的先验估计，以及所述模型不确定度对应的切换增益项，计算当前帧的第一后验状态估计；
[0253]
第二状态估计模块400，用于利用所述多普勒雷达在当前帧的径向速度观测变量构造得到的距离-径向速度积伪量测变量，对所述当前帧的第一后验状态估计进行更新，得到当前帧的第二后验状态估计，将所述当前帧的第二后验状态估计作为当前帧最终的目标状态估计。
[0254]
可选地，所述策略执行模块包括：
[0255]
参数计算单元，用于根据当前帧的先验估计以及所述线性位置观测变量，计算所述多普勒雷达观测模型的广义时变平滑层参数；
[0256]
模型不确定度判决单元，用于在所述广义时变平滑层参数大于预设的平滑层参数时，判决所述目标运动模型在当前帧为高模型不确定度；在所述广义时变平滑层参数小于等于预设的平滑层参数时，判决所述目标运动模型在当前帧为低模型不确定度。
[0257]
可选地，所述参数计算单元包括：
[0258]
第一计算子单元，用于计算所述当前帧的先验估计，所述当前帧的先验估计包括目标状态先验估计、目标位置量测的先验估计和先验状态估计协方差阵；
[0259]
第二计算子单元，用于根据所述线性位置观测变量、第一转换偏差以及所述目标位置量测的先验估计，计算先验位置观测误差；
[0260]
第三计算子单元，用于根据所述先验位置观测误差，计算混合误差项；
[0261]
第四计算子单元，用于根据所述先验状态估计协方差阵、所述混合误差项、所述第一转换误差协方差，计算所述当前帧的广义时变平滑层参数。
[0262]
可选地，所述策略执行模块还包括：
[0263]
标准切换增益项计算单元，用于在所述目标运动模型在当前帧为高模型不确定度时，根据所述预设的平滑层参数，计算标准切换增益项；
[0264]
最优切换增益项计算单元，用于在所述目标运动模型在当前帧为低模型不确定度时，根据所述广义时变平滑层参数，计算最优切换增益项。
[0265]
可选地，所述第一状态估计模块包括：
[0266]
第一估计单元，用于在所述目标运动模型在当前帧为高模型不确定度时，根据所述当前帧的先验估计以及所述标准切换增益项，计算所述当前帧的第一后验状态估计；
[0267]
第二估计单元，用于在所述目标运动模型在当前帧为低模型不确定度时，根据所述当前帧的先验估计以及所述最优切换增益项，计算所述当前帧的第一后验状态估计。
[0268]
可选地，所述第一状态估计模块包括：
[0269]
状态估计计算单元，用于根据所述目标状态先验估计、所述先验状态估计协方差阵、所述先验位置观测误差，以及所述目标运动模型的模型不确定度对应的切换增益项，计算目标状态后验估计与后验状态估计协方差阵；将所述目标状态后验估计与所述后验状态估计协方差阵作为所述当前帧的第一后验状态估计。
[0270]
可选地，所述线性转换模块包括：
[0271]
第一转换单元，用于采用修正的无偏量测转换方法，将所述多普勒雷达在三维球坐标或者二维极坐标下的所述位置观测变量，转换为笛卡尔坐标系下的所述线性位置观测变量，并计算第一转换偏差和第一转换误差协方差。
[0272]
可选地，所述第二状态估计模块包括：
[0273]
第二转换单元，用于根据所述多普勒雷达在当前帧的距离观测变量以及径向速度观测变量，构造距离-径向速度积伪量测变量，并采用修正的无偏转换方法计算所述距离-径向速度积伪量测变量的第二转换偏差、转换误差方差以及第二转换误差协方差；
[0274]
去相关处理单元，用于对所述距离-径向速度积伪量测变量进行去相关处理；
[0275]
最终估计单元，用于将去相关处理后的所述距离-径向速度积伪量测变量作为输入，构造局部近似的线性最小均方误差估计器，对所述当前帧的第一后验状态估计进行更新，将更新后得到的所述当前帧的第二后验状态估计作为当前帧最终的目标状态估计。
[0276]
本说明书中的各个实施例均采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似的部分互相参见即可。
[0277]
本领域内的技术人员应明白，本技术实施例的实施例可提供为方法、系统、或计算机程序产品。因此，本技术实施例可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本技术实施例可采用在一个或多个其中包含有计算机可
用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、cd-rom、光学存储器等)上实施的计算机程序产品的形式。
[0278]
本技术实施例是参照根据本技术实施例的方法、终端设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理终端设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理终端设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的系统。
[0279]
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理终端设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令系统的制造品，该指令系统实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
[0280]
这些计算机程序指令也可装载到计算机或其他可编程数据处理终端设备上，使得在计算机或其他可编程终端设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程终端设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
[0281]
尽管已描述了本技术实施例的优选实施例，但本领域内的技术人员一旦得知了基本创造性概念，则可对这些实施例做出另外的变更和修改。所以，所附权利要求意欲解释为包括优选实施例以及落入本技术实施例范围的所有变更和修改。
[0282]
最后，还需要说明的是，在本文中，诸如第一和第二等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者终端设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者终端设备所固有的要素。在没有更多限制的情况下，由语句“包括一个
……”
限定的要素，并不排除在包括所述要素的过程、方法、物品或者终端设备中还存在另外的相同要素。
[0283]
本文中应用了具体个例对本技术的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本技术的方法及其核心思想；同时，对于本领域的一般技术人员，依据本技术的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本技术的限制。