CN105652250B - 一种基于双层期望最大化的机动目标跟踪技术 - Google Patents
一种基于双层期望最大化的机动目标跟踪技术 Download PDFInfo
- Publication number
- CN105652250B CN105652250B CN201610027471.0A CN201610027471A CN105652250B CN 105652250 B CN105652250 B CN 105652250B CN 201610027471 A CN201610027471 A CN 201610027471A CN 105652250 B CN105652250 B CN 105652250B
- Authority
- CN
- China
- Prior art keywords
- mrow
- msub
- mover
- maneuvering target
- additivity
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
- G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
- G01S7/36—Means for anti-jamming, e.g. ECCM, i.e. electronic counter-counter measures
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
- G01S13/66—Radar-tracking systems; Analogous systems
- G01S13/72—Radar-tracking systems; Analogous systems for two-dimensional tracking, e.g. combination of angle and range tracking, track-while-scan radar
- G01S13/723—Radar-tracking systems; Analogous systems for two-dimensional tracking, e.g. combination of angle and range tracking, track-while-scan radar by using numerical data
Abstract
本发明属于目标跟踪技术领域，公开了一种基于双层期望最大化的机动目标跟踪技术。包括：一，通过N个雷达分别实时测量对应得到机动目标的N个雷达量测向量雷达量测向量y包括机动目标与雷达之间的距离、方位角、机动目标与雷达之间的距离变化率；二，通过并行执行第一层期望最大化算法得到机动目标状态向量x估计集合和加性未知干扰伪量测θ集合并将传输给第二层期望最大化算法；三，第二层期望最大化算法接收到后，用混合高斯分布去拟合的前一二阶距，得到加性未知干扰伪量测集合的均值和协方差四，利用的均值和协方差信息，并通过kalman滤波得到状态估计值该技术可保证参数辨识的解析性和收敛性，且该技术可提高目标状态估计精度。
Description
技术领域
本发明属于目标跟踪技术领域，特别涉及一种基于双层期望最大化的机动目标跟踪技术。
背景技术
目标跟踪在军事领域中占据着重要的位置，只有可靠且精确的跟踪才能有效地对目标实施打击。随着国防军事领域的发展，战斗机的机动性能得到了很大的提高。目标的机动性加大了目标跟踪的难度，主要表现在目标机动状态的不确定性引起数学模型的建模误差，该建模误差可表示为加性未知干扰且通常与状态相关。在此背景下，如何提高机动目标的状态跟踪精度成为一个越来越重要的问题，因此迫切需要研究性能更为优越的跟踪滤波方法。
现有针对加性未知干扰存在下的机动目标跟踪问题的解决方法主要有扩维法，鲁棒法，交互多模型(Interacting Multiple Model,IMM)以及传统的期望最大化(Expectation-maximuzation，EM)算法等技术方案。扩维法将加性未知干扰当作状态并扩充到状态向量中从而将加性未知干扰与状态一起进行估计；鲁棒法一定程度上消除了模型误差对估计结果带来的影响，但其结果具有保守性；IMM方法是一种次优的具有高费效比的估计技术，它使用多个不同的运动模型分别匹配目标的不同运动状态,但计算复杂度较高；传统EM方法将状态估计与加性未知干扰辨识看作联合优化问题，然后通过E-步更新状态估计与M-步优化辨识干扰的迭代计算来实现。但上述技术方案都将加性未知干扰看作常值，并未考虑加性未知干扰与状态相关的情况，因此它们都仅用加性未知干扰的均值去校正状态估计。若干扰与状态相关，此时干扰具有均值和协方差特性，则上述技术方案的估计精度并不能达到要求。因此，本发明提出一种基于双层EM实现联合状态估计与加性未知干扰一二阶矩(均值和协方差)拟合辨识的机动目标跟踪技术，并使用加性未知干扰的均值和协方差辨识结果同时去校正状态从而提高目标跟踪估计精度。
发明内容
本发明的目的是提供一种基于双层期望最大化的机动目标跟踪技术，该技术可提高机动目标跟踪精度，解决机动目标跟踪中加性未知干扰与状态耦合的问题。
为达到以上目的，本发明采用以下技术方案予以实现。
一种基于双层期望最大化算法的机动目标跟踪技术，其特征在于：包括以下步骤：
步骤一，通过N个雷达分别实时测量对应得到机动目标的N个雷达量测向量k表示时刻值，k＝1、2、3……；雷达量测向量y包括机动目标与雷达之间的距离、方位角、机动目标与雷达之间的距离变化率；
步骤二，所述N个雷达量测向量通过并行执行第一层期望最大化算法得到机动目标状态向量x估计集合和加性未知干扰伪量测θ集合并将加性未知干扰伪量测集合传输给第二层期望最大化算法；
步骤三，第二层期望最大化算法接收到第一层期望最大化算法传输的加性未知干扰伪量测集合后，用混合高斯分布去拟合加性未知干扰伪量测集合的前一二阶距，得到加性未知干扰伪量测集合的均值和协方差
步骤四，利用加性未知干扰伪量测集合的均值和协方差信息，并通过kalman滤波得到状态估计值
上述技术方案的特点和进一步改进：
(1)进一步的，在步骤二中：
2a、将机动目标运动模型及其中一个雷达量测模型建模为下述带有加性未知干扰a、b的线性模型，表示为：
其中，xk为k时刻的机动目标状态向量，yk为k时刻的雷达量测向量，Fk表示k时刻的状态转移矩阵，Hk为k时刻的观测矩阵，wk表示k时刻的系统噪声，vk为k时刻的量测噪声，且即系统噪声wk的方差为Qk，量测噪声vk的方差为Rk，Mk、Nk为k时刻的适当维数的已知矩阵，加性未知干扰a、b是关于机动目标状态向量x的函数，并令θ＝{a,b}；
第一层期望最大化算法将机动目标状态向量x看作隐含或缺失数据，并计算机动目标状态向量x的期望值其中，r表示迭代次数，表示第r次迭代所辨识出的加性未知干扰伪量测；
期望值的具体计算过程如下：根据贝叶斯规则和模型一阶Markov性质对完全数据的对数似然函数进行等价分解，
其中，
l表示窗长，取值为1至k之间的任意整数；
然后对计算关于概率密度函数的数学期望，该数学期望即为是待辨识加性未知干扰a,b的函数；
2b、通过前向-后向通道的固定区间平滑器计算在当前未知干扰伪量测辨识结果下机动目标状态向量x在区间[k-l,k]内的平滑值和对应的协方差Pi,i|k-l:k，其中i表示区间[k-l,k]内的时刻值，其取值：i＝{k-l,…,k}，前向-后向通道的固定区间平滑器具体实现如下，
该平滑器中的前向估计值对应的协方差Pi,i|k-l:i和后向估计值对应的协方差Pi,i|i+1:k均可通过截断区间Kalman滤波器估计得到；
2c、然后直接计算数学期望关于θ的导数并令它为零，即和可解析得到第r+1次迭代所辨识出的加性未知干扰伪量测是第一层期望最大化算法基于真实量测信息对加性未知干扰的伪量测，辨识结果如下：
式中，Mi-1表示i-1时刻矩阵M，表示i-1时刻矩阵M的转置；
Ni表示i时刻矩阵N，表示i时刻矩阵N的转置；
表示i-1时刻矩阵Q的逆，
Fi-1表示i-1时刻的状态转移矩阵，
表示i时刻矩阵R的逆，
Hi表示i时刻的观测矩阵；
2d、N个雷达量测向量通过并行迭代执行上述2a-2c直至满足迭代终止条件可计算得到机动目标状态向量xk估计集合和加性未知干扰伪量测θ集合机动目标状态向量估计值等于在k时刻最后一次迭代所得的机动目标状态向量平滑值即且其中r max表示最大迭代次数；将机动目标状态向量xk估计集合采用联邦滤波器得到基于标准期望最大化算法的机动目标状态向量xk估计融合结果
(2)进一步的，在步骤2d中所述迭代终止条件为连续两次迭代的完全数据对数似然函数值相差不超过设定的门限值或者达到设定的最大迭代次数。
(3)进一步的，在步骤三中：
3a、推导数学期望的表达式其中表示M个高斯成分中每个高斯分项的权值αj、均值μj和协方差Σj的集合，式中 表示第m次迭代辨识出的参数集合，式中：j表示高斯成分，j＝1、2、3…M；
取作为观测数据，并取表示隐含或缺失数据，其中zi∈{1,2,…,M},i＝1,2,…,N且zi＝j表示样本由第j个高斯成分产生，则完全数据的似然函数表示为然后计算似然函数关于概率密度函数的)数学期望，则可得数学期望 表示对计算关于概率密度函数的期望；
3b、计算概率密度函数 表示在第m次迭代所辨识出参数的前提下由第j个混合高斯分项产生的概率，表达式中的概率密度函数服从多维正态分布，即表示在第m次迭代所辨识出的第j个高斯分项的权值；表示在第m次迭代所辨识出的第j个高斯分项的均值和协方差；表示在条件下由第j个混合高斯分项产生的概率，该概率服从多维正态分布；
3c、辨识出使期望取最大值的参数集对于每个高斯分项的均值μj和协方差Σj的辨识，可直接用对均值μj和协方差Σj求导并取极值，即 得到，
对于权值αj的辨识，用拉格朗日乘子法构造新的目标函数β表示拉格朗日乘子，该目标函数对αj求导并结合约束条件可辨识出αj，每个高斯分项的权值αj、均值μj和协方差Σj的辨识结果分别为：
权值
均值
协方差
其中，为的转置；
3d、通过并行迭代执行上述步骤3b和3c直至满足迭代终止条件可辨识出混合高斯中每个高斯分项的权值αj、均值μj和协方差Σj，的概率密度函数表示为根据均值和协方差的定义可推导出组数据混合高斯分布拟合的均值和协方差
(4)进一步的，在步骤3d中所述迭代终止条件为连续两次迭代的完全数据对数似然函数值相差不超过设定的门限值或达到设定的最大迭代次数。
本发明的基于双层期望最大化算法的机动目标跟踪技术，其中第一层EM为联合状态估计与加性未知干扰伪量测辨识，第二层EM为混合高斯分布拟合加性未知干扰的前一二阶矩(均值和协方差)；本发明的上述两个主要部分均是基于EM框架，且两层EM在辨识加性未知干扰和拟合加性未知干扰一二阶矩时的期望最大化均为凸优化过程，这就保证了双层EM的参数辨识解析性和收敛性。由于双层EM技术利用干扰的均值和协方差信息同时去校正状态估计，针对加性未知干扰明显具有协方差特性的目标模型，该技术可提高目标状态估计精度。
附图说明
图1为本发明的一种基于双层期望最大化的机动目标跟踪技术的流程图；
图2a为机动目标水平加速度辨识结果，图2b为机动目标垂直加速度辨识结果；
图3为混合高斯中各高斯分项的拟合结果图；
图4a为标准EM与双层EM对机动目标状态在水平位置方向估计精度比较；图4b为标准EM与双层EM对机动目标状态在水平速度方向估计精度比较；图4c为标准EM与双层EM对机动目标状态在垂直位置方向估计精度比较；图4d为标准EM与双层EM对机动目标状态在垂直速度方向估计精度比较。
具体实施方式
下面结合附图和具体实施方式对本发明进行详细说明。
参照图1，为发明的一种基于双层期望最大化的机动目标跟踪技术的流程图。本发明的基于双层期望最大化算法的机动目标跟踪技术，包括以下步骤：
步骤一，通过N个雷达分别实时测量对应得到机动目标的N个雷达量测向量k表示时刻值，k＝1、2、3……；雷达量测向量y包括机动目标与雷达之间的距离、方位角、机动目标与雷达之间的距离变化率。
步骤二，所述N个雷达量测向量通过并行执行第一层期望最大化算法得到机动目标状态向量x估计集合和加性未知干扰伪量测θ集合并将加性未知干扰伪量测集合传输给第二层期望最大化算法。
2a、联合状态估计与未知干扰伪量测EM框架中期望表达式的推导：将机动目标运动模型及其中一个雷达量测模型建模为下述带有加性未知干扰a、b的线性模型，表示为：
其中，xk为k时刻的机动目标状态向量，yk为k时刻的雷达量测向量，Fk表示k时刻的状态转移矩阵，Hk为k时刻的观测矩阵，wk表示k时刻的系统噪声，vk为k时刻的量测噪声，且即系统噪声wk的方差为Qk，量测噪声vk的方差为Rk，Mk、Nk为k时刻的适当维数的已知矩阵，加性未知干扰a、b是关于机动目标状态向量x的函数，并令θ＝{a,b}。
EM可用于解决含有不完全数据的最大似然估计问题，则EM框架适用于联合估计与辨识问题的求解。为实现EM求解该联合优化问题，第一层期望最大化算法将机动目标状态向量x看作隐含或缺失数据，并计算机动目标状态向量x的期望值其中，r表示迭代次数，表示第r次迭代所辨识出的加性未知干扰伪量测。
期望值的具体计算过程如下：首先，根据贝叶斯规则和模型一阶Markov性质对完全数据的对数似然函数进行等价分解，
其中，
l表示窗长，取值为1至k之间的任意整数；
然后，对计算关于概率密度函数的数学期望，该数学期望即为其是待辨识加性未知干扰a,b的函数；这是实现EM优化的关键和基础。
2b、状态估计：为计算步骤2a推导出的期望值需首先计算在当前加性未知干扰伪量测辨识结果下状态在区间[k-l,k]内的估计值和对应的协方差Pi,i|k-l:k，该问题是一个状态平滑问题，可通过固定区间平滑器实现。通过前向-后向通道的固定区间平滑器计算在当前未知干扰伪量测辨识结果下机动目标状态向量x在区间[k-l,k]内的平滑值和对应的协方差Pi,i|k-l:k，其中i表示区间[k-l,k]内的时刻值，其取值：i＝{k-l,…,k}，前向-后向通道的固定区间平滑器具体实现如下：
式中，机动目标状态向量x在区间[k-l,i]内的估计值为对应的协方差为Pi,i|k-l:i；为机动目标状态向量x在区间[i+1,k]内的估计值，对应的协方差为Pi,i|i+1:k。
该平滑器中的前向估计值对应的协方差Pi,i|k-l:i和后向估计值对应的协方差Pi,i|i+1:k均可通过截断区间Kalman滤波器估计得到。
2c、加性未知干扰伪量测辨识：根据步骤2b所得的状态估计值，可以计算出步骤2a推导的期望然后直接计算数学期望关于θ的导数并令它为零，即和可解析得到第r+1次迭代所辨识出的加性未知干扰伪量测是第一层期望最大化算法基于真实量测信息对加性未知干扰的伪量测，辨识结果如下：
式中，Mi-1表示i-1时刻矩阵M，表示i-1时刻矩阵M的转置；
Ni表示i时刻矩阵N，表示i时刻矩阵N的转置；
表示i-1时刻矩阵Q的逆，
Fi-1表示i-1时刻的状态转移矩阵，
表示i时刻矩阵R的逆，
Hi表示i时刻的观测矩阵；
文中，带下标的Mk、Nk表示矩阵，而M、N为标量。
2d、基于标准EM的状态估计融合：为了在步骤三中能够拟合出的均值和协方差需要从步骤2c中获得多组加性未知干扰伪量测。所以我们提出利用多雷达量测，N个雷达量测向量通过并行迭代执行上述2a-2c直至满足迭代终止条件可计算得到机动目标状态向量xk估计集合和加性未知干扰伪量测θ集合机动目标状态向量估计值等于在k时刻最后一次迭代所得的机动目标状态向量平滑值即且其中r max表示最大迭代次数；将机动目标状态向量xk估计集合采用联邦滤波器可得到基于标准期望最大化算法的机动目标状态向量xk估计融合结果显然，这种方式的融合结果是仅用均值去校正状态的估计结果。
在步骤2d中所述迭代终止条件为连续两次迭代的完全数据对数似然函数值相差不超过设定的门限值或者达到设定的最大迭代次数；所设定的门限值和最大迭代次数可用于平衡估计精度和计算量，可根据具体指标而设定。
步骤三，第二层期望最大化算法接收到第一层期望最大化算法传输的加性未知干扰伪量测集合后，用混合高斯分布去拟合加性未知干扰伪量测集合的前一二阶距，得到加性未知干扰伪量测集合的均值和协方差混合高斯是用来逼近任意分布的一项重要技术。
3a、为了在EM框架下实现混合高斯拟合，首先推导数学期望的表达式其中表示M个高斯成分中每个高斯分项的权值αj、均值μj和协方差Σj的集合，式中表示第m次迭代辨识出的参数集合，式中：j表示高斯成分，j＝1、2、3…M。
取作为观测数据，并取表示隐含或缺失数据，其中zi∈{1,2,…,M},i＝1,2,…，N且zi＝j表示样本由第j个高斯成分产生，则完全数据的似然函数表示为然后计算似然函数关于概率密度函数的)数学期望，则可得数学期望 表示对计算关于概率密度函数的期望。
3b、为了计算步骤3a中推导出的期望值需要计算概率密度函数 表示在第m次迭代所辨识出参数的前提下由第j个混合高斯分项产生的概率，表达式中的概率密度函数服从多维正态分布，即 表示在第m次迭代所辨识出的第j个高斯分项的权值；表示在第m次迭代所辨识出的第j个高斯分项的均值和协方差；表示在条件下由第j个混合高斯分项产生的概率，该概率服从多维正态分布。可根据多维正态分布的定义直接写出其表达式，进而计算出的值。
3c、辨识出使期望取最大值的参数集对于每个高斯分项的均值μj和协方差Σj的辨识，可直接用对均值μj和协方差Σj求导并取极值，即 得到，
对于权值αj的辨识，由于等值约束的存在，直接求导的方法不可行，用拉格朗日乘子法构造新的目标函数β表示拉格朗日乘子，该目标函数对αj求导并结合约束条件可辨识出αj，第m+1次迭代所辨识出的每个高斯分项的权值αj,m+1、均值μj,m+1和协方差Σj,m+1的辨识结果分别为：
权值
均值
协方差
其中，为的转置；
3d、通过并行迭代执行上述步骤3b和3c直至满足迭代终止条件可辨识出混合高斯中每个高斯分项的权值αj、均值μj和协方差的概率密度函数表示为根据均值和协方差的定义可推导出组数据混合高斯分布拟合的均值和协方差
在步骤3d中所述迭代终止条件为连续两次迭代的完全数据对数似然函数值相差不超过设定的门限值或达到设定的最大迭代次数；所设定的门限值和最大迭代次数可平衡估计精度和计算量，可根据具体指标而设定。
步骤四，同时利用加性未知干扰伪量测集合的均值和协方差信息，并通过kalman滤波可得到基于双层期望最大化的目标状态估计值
步骤二和步骤三中两个EM框架的顺序组合，就是所提出的双层EM技术，它不仅可辨识出加性未知干扰伪量测的值还可拟合出其均值和协方差，相比于传统EM技术仅用加性未知干扰的均值去校正状态估计，我们发明的双层EM技术在加性未知干扰与状态相关时，同时利用被拟合出的加性未知干扰均值和协方差去联合校正状态的均值和协方差，从而提高机动目标跟踪精度。具体实现是同时利用加性未知干扰伪量测的均值和协方差进行kalman滤波，其滤波的结果即基于双层EM技术的状态估计结果。
实施例：机动目标跟踪
步骤1：机动目标跟踪的模型可表示为：
动态模型：
其中，分别表示水平和垂直方向的位置和速度；表示水平和垂直方向上的加速度，对于机动目标尤其是非合作机动目标a未知，可建模为未知干扰；Ts为采样时间。
雷达量测模型：
其中，雷达量测向量yk包括机动目标与雷达之间的距离、方位角、机动目标与雷达之间的距离变化率；
步骤2：第一层EM：仅用均值校正机动目标状态并辨识出机动目标加速度。
图2a反映了第一层EM辨识出的机动目标水平加速度对真实值的跟踪效果。
图2b反映了第一层EM辨识出的机动目标垂直加速度对真实值的跟踪效果。
步骤3：第二层EM：用混合高斯拟合机动目标加速度的均值和协方差。
图3表示了双层EM技术中高斯分项的拟合结果。
步骤4：用机动目标加速度的均值和协方差联合校正目标状态。
从图4a、4b、4c、4d可看出，基于所发明的双层EM技术估计的机动目标状态在速度分量上的精度与基于标准EM技术估计的结果基本一致，但双层EM技术对位置分量的估计结果明显优于标准EM技术。原因在于发明的双层EM技术利用未知干扰的均值和协方差特性联合校正机动目标状态，从而提高了目标跟踪精度。
尽管以上结合附图对本发明的实施方案进行了描述，但是本发明并不局限于上述的具体实施方案和应用领域，上述的具体实施方案仅仅是示意性的、指导性的，而不是限制性的。本领域的普通技术人员在说明书的启示下，在不脱离本发明权利要求所保护的范围的情况下，还可以做出很多种的形式，这些均属于本发明保护之列。
Claims (5)
1.一种基于双层期望最大化算法的机动目标跟踪方法，其特征在于：包括以下步骤：
步骤一，通过N个雷达分别实时测量对应得到机动目标的N个雷达量测向量k表示时刻值，k＝1、2、3……；雷达量测向量y包括机动目标与雷达之间的距离、方位角、机动目标与雷达之间的距离变化率；
步骤二，所述N个雷达量测向量通过并行执行第一层期望最大化算法得到机动目标状态向量x估计集合和加性未知干扰伪量测θ集合并将加性未知干扰伪量测集合传输给第二层期望最大化算法；
步骤三，第二层期望最大化算法接收到第一层期望最大化算法传输的加性未知干扰伪量测集合后，用混合高斯分布去拟合加性未知干扰伪量测集合的前一二阶矩，得到加性未知干扰伪量测集合的均值和协方差
步骤四，利用加性未知干扰伪量测集合的均值和协方差信息，并通过kalman滤波得到状态估计值
2.如权利要求1所述的一种基于双层期望最大化算法的机动目标跟踪方法，其特征在于：
在步骤二中：
2a、将机动目标运动模型及其中一个雷达量测模型建模为下述带有加性未知干扰a、b的线性模型，表示为：
<mrow>
<mfenced open = "{" close = "">
<mtable>
<mtr>
<mtd>
<mrow>
<msub>
<mi>x</mi>
<mi>k</mi>
</msub>
<mo>=</mo>
<msub>
<mi>F</mi>
<mi>k</mi>
</msub>
<msub>
<mi>x</mi>
<mrow>
<mi>k</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
<mo>+</mo>
<msub>
<mi>w</mi>
<mrow>
<mi>k</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
<mo>+</mo>
<msub>
<mi>M</mi>
<mrow>
<mi>k</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
<msub>
<mi>a</mi>
<mrow>
<mi>k</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<msub>
<mi>x</mi>
<mrow>
<mi>k</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msub>
<mi>y</mi>
<mi>k</mi>
</msub>
<mo>=</mo>
<msub>
<mi>H</mi>
<mi>k</mi>
</msub>
<msub>
<mi>x</mi>
<mi>k</mi>
</msub>
<mo>+</mo>
<msub>
<mi>v</mi>
<mi>k</mi>
</msub>
<mo>+</mo>
<msub>
<mi>N</mi>
<mi>k</mi>
</msub>
<msub>
<mi>b</mi>
<mi>k</mi>
</msub>
<mrow>
<mo>(</mo>
<msub>
<mi>x</mi>
<mi>k</mi>
</msub>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>,</mo>
</mrow>
其中，xk为k时刻的机动目标状态向量，yk为k时刻的雷达量测向量，Fk表示k时刻的状态转移矩阵，Hk为k时刻的观测矩阵，wk表示k时刻的系统噪声，vk为k时刻的量测噪声，且即系统噪声wk的方差为Qk，量测噪声vk的方差为Rk，Mk、Nk为k时刻的适当维数的已知矩阵，加性未知干扰a、b是关于机动目标状态向量x的函数，并令θ＝{a,b}；
第一层期望最大化算法将机动目标状态向量x看作隐含或缺失数据，并计算机动目标状态向量x的期望值其中，r表示迭代次数，表示第r次迭代所辨识出的加性未知干扰伪量测；
期望值的具体计算过程如下：根据贝叶斯规则和模型一阶Markov性质对完全数据的对数似然函数进行等价分解，
其中，
l表示窗长，取值为1至k之间的任意整数；
然后对计算关于概率密度函数的数学期望，该数学期望即为是待辨识加性未知干扰a,b的函数；
2b、通过前向-后向通道的固定区间平滑器计算在当前未知干扰伪量测辨识结果下机动目标状态向量x在区间[k-l,k]内的平滑值和对应的协方差Pi,i|k-l:k，其中i表示区间[k-l,k]内的时刻值，其取值：i＝{k-l,…,k}，前向-后向通道的固定区间平滑器具体实现如下：
<mrow>
<mfenced open = "{" close = "">
<mtable>
<mtr>
<mtd>
<msub>
<mover>
<mi>x</mi>
<mo>^</mo>
</mover>
<mrow>
<mi>i</mi>
<mo>|</mo>
<mi>k</mi>
<mo>-</mo>
<mi>l</mi>
<mo>:</mo>
<mi>k</mi>
</mrow>
</msub>
<mo>=</mo>
<msub>
<mi>P</mi>
<mrow>
<mi>i</mi>
<mo>,</mo>
<mi>i</mi>
<mo>|</mo>
<mi>k</mi>
<mo>-</mo>
<mi>l</mi>
<mo>:</mo>
<mi>k</mi>
</mrow>
</msub>
<mo>&lsqb;</mo>
<msup>
<mrow>
<mo>(</mo>
<msub>
<mi>P</mi>
<mrow>
<mi>i</mi>
<mo>,</mo>
<mi>i</mi>
<mo>|</mo>
<mi>k</mi>
<mo>-</mo>
<mi>l</mi>
<mo>:</mo>
<mi>i</mi>
</mrow>
</msub>
<mo>)</mo>
</mrow>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msup>
<msub>
<mover>
<mi>x</mi>
<mo>^</mo>
</mover>
<mrow>
<mi>i</mi>
<mo>|</mo>
<mi>k</mi>
<mo>-</mo>
<mi>l</mi>
<mo>:</mo>
<mi>i</mi>
</mrow>
</msub>
<mo>+</mo>
<msup>
<mrow>
<mo>(</mo>
<msub>
<mi>P</mi>
<mrow>
<mi>i</mi>
<mo>,</mo>
<mi>i</mi>
<mo>|</mo>
<mi>i</mi>
<mo>+</mo>
<mn>1</mn>
<mo>:</mo>
<mi>k</mi>
</mrow>
</msub>
<mo>)</mo>
</mrow>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msup>
<msub>
<mover>
<mi>x</mi>
<mo>^</mo>
</mover>
<mrow>
<mi>i</mi>
<mo>|</mo>
<mi>i</mi>
<mo>+</mo>
<mn>1</mn>
<mo>:</mo>
<mi>k</mi>
</mrow>
</msub>
<mo>&rsqb;</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>P</mi>
<mrow>
<mi>i</mi>
<mo>,</mo>
<mi>i</mi>
<mo>|</mo>
<mi>k</mi>
<mo>-</mo>
<mi>l</mi>
<mo>:</mo>
<mi>k</mi>
</mrow>
</msub>
<mo>=</mo>
<msup>
<mrow>
<mo>&lsqb;</mo>
<mrow>
<msup>
<mrow>
<mo>(</mo>
<msub>
<mi>P</mi>
<mrow>
<mi>i</mi>
<mo>,</mo>
<mi>i</mi>
<mo>|</mo>
<mi>k</mi>
<mo>-</mo>
<mi>l</mi>
<mo>:</mo>
<mi>i</mi>
</mrow>
</msub>
<mo>)</mo>
</mrow>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msup>
<mo>+</mo>
<msup>
<mrow>
<mo>(</mo>
<msub>
<mi>P</mi>
<mrow>
<mi>i</mi>
<mo>,</mo>
<mi>i</mi>
<mo>|</mo>
<mi>i</mi>
<mo>+</mo>
<mn>1</mn>
<mo>:</mo>
<mi>k</mi>
</mrow>
</msub>
<mo>)</mo>
</mrow>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msup>
</mrow>
<mo>&rsqb;</mo>
</mrow>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msup>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>,</mo>
</mrow>
该平滑器中的前向估计值对应的协方差Pi,i|k-l:i和后向估计值对应的协方差Pi,i|i+1:k均可通过截断区间Kalman滤波器估计得到；
2c、直接计算数学期望关于θ的导数并令它为零，即和可解析得到第r+1次迭代所辨识出的加性未知干扰伪量测 是第一层期望最大化算法基于真实量测信息对加性未知干扰的伪量测，辨识结果如下：
<mrow>
<msub>
<mover>
<mi>a</mi>
<mo>^</mo>
</mover>
<mrow>
<mi>r</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</msub>
<mo>=</mo>
<msup>
<mrow>
<mo>(</mo>
<mrow>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mi>k</mi>
<mo>-</mo>
<mi>l</mi>
</mrow>
<mi>k</mi>
</munderover>
<msubsup>
<mi>M</mi>
<mrow>
<mi>i</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
<mi>T</mi>
</msubsup>
<msubsup>
<mi>Q</mi>
<mrow>
<mi>i</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msubsup>
<msub>
<mi>M</mi>
<mrow>
<mi>i</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
</mrow>
<mo>)</mo>
</mrow>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msup>
<mo>&CenterDot;</mo>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mi>k</mi>
<mo>-</mo>
<mi>l</mi>
</mrow>
<mi>k</mi>
</munderover>
<msubsup>
<mi>M</mi>
<mrow>
<mi>i</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
<mi>T</mi>
</msubsup>
<msubsup>
<mi>Q</mi>
<mrow>
<mi>i</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msubsup>
<mrow>
<mo>(</mo>
<msub>
<mover>
<mi>x</mi>
<mo>^</mo>
</mover>
<mrow>
<mi>i</mi>
<mo>|</mo>
<mi>k</mi>
<mo>-</mo>
<mi>l</mi>
<mo>:</mo>
<mi>k</mi>
</mrow>
</msub>
<mo>-</mo>
<msub>
<mi>F</mi>
<mrow>
<mi>i</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
<msub>
<mover>
<mi>x</mi>
<mo>^</mo>
</mover>
<mrow>
<mi>i</mi>
<mo>-</mo>
<mn>1</mn>
<mo>|</mo>
<mi>k</mi>
<mo>-</mo>
<mi>l</mi>
<mo>:</mo>
<mi>k</mi>
</mrow>
</msub>
<mo>)</mo>
</mrow>
<mo>,</mo>
</mrow>
<mrow>
<msub>
<mover>
<mi>b</mi>
<mo>^</mo>
</mover>
<mrow>
<mi>r</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</msub>
<mo>=</mo>
<msup>
<mrow>
<mo>(</mo>
<mrow>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mi>k</mi>
<mo>-</mo>
<mi>l</mi>
</mrow>
<mi>k</mi>
</munderover>
<msubsup>
<mi>N</mi>
<mi>i</mi>
<mi>T</mi>
</msubsup>
<msubsup>
<mi>R</mi>
<mi>i</mi>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msubsup>
<msub>
<mi>N</mi>
<mi>i</mi>
</msub>
</mrow>
<mo>)</mo>
</mrow>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msup>
<mo>&CenterDot;</mo>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mi>k</mi>
<mo>-</mo>
<mi>l</mi>
</mrow>
<mi>k</mi>
</munderover>
<msubsup>
<mi>N</mi>
<mi>i</mi>
<mi>T</mi>
</msubsup>
<msubsup>
<mi>R</mi>
<mi>i</mi>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msubsup>
<mrow>
<mo>(</mo>
<msub>
<mi>y</mi>
<mi>i</mi>
</msub>
<mo>-</mo>
<msub>
<mi>H</mi>
<mi>i</mi>
</msub>
<msub>
<mover>
<mi>x</mi>
<mo>^</mo>
</mover>
<mrow>
<mi>i</mi>
<mo>|</mo>
<mi>k</mi>
<mo>-</mo>
<mi>l</mi>
<mo>:</mo>
<mi>k</mi>
</mrow>
</msub>
<mo>)</mo>
</mrow>
<mo>,</mo>
</mrow>
式中，Mi-1表示i-1时刻矩阵M，表示i-1时刻矩阵M的转置；
Ni表示i时刻矩阵N，表示i时刻矩阵N的转置；
表示i-1时刻矩阵Q的逆，
Fi-1表示i-1时刻的状态转移矩阵，
表示i时刻矩阵R的逆，
Hi表示i时刻的观测矩阵；
2d、N个雷达量测向量通过并行迭代执行上述2a-2c直至满足迭代终止条件可计算得到机动目标状态向量xk估计集合和加性未知干扰伪量测θ集合机动目标状态向量估计值等于在k时刻最后一次迭代所得的机动目标状态向量平滑值即且其中rmax表示最大迭代次数；将机动目标状态向量xk估计集合采用联邦滤波器得到基于标准期望最大化算法的机动目标状态向量xk估计融合结果
3.如权利要求2所述的一种基于双层期望最大化算法的机动目标跟踪方法，其特征在于：在步骤2d中所述迭代终止条件为连续两次迭代的完全数据对数似然函数值相差不超过设定的门限值或者达到设定的最大迭代次数。
4.如权利要求1-3任一项所述的一种基于双层期望最大化算法的机动目标跟踪方法，其特征在于：
在步骤三中：
3a、推导数学期望的表达式其中表示M个高斯成分中每个高斯分项的权值αj、均值μj和协方差Σj的集合，式中 表示第m次迭代辨识出的参数集合，式中：j表示高斯成分，j＝1、2、3…M；
取作为观测数据，并取表示隐含或缺失数据，其中zi∈{1,2,…,M},i＝1,2,…,N且zi＝j表示样本由第j个高斯成分产生，则完全数据的似然函数表示为然后计算似然函数关于概率密度函数的数学期望，则可得数学期望 表示对计算关于概率密度函数的期望；
3b、计算概率密度函数 表示在第m次迭代所辨识出参数的前提下由第j个混合高斯分项产生的概率，表达式中的概率密度函数服从多维正态分布，即 表示在第m次迭代所辨识出的第j个高斯分项的权值；表示在第m次迭代所辨识出的第j个高斯分项的均值和协方差；表示在条件下由第j个混合高斯分项产生的概率，该概率服从多维正态分布；
3c、辨识出使期望取最大值的参数集对于每个高斯分项的均值μj和协方差Σj的辨识，可直接用对均值μj和协方差Σj求导并取极值，即得到，
对于权值αj的辨识，用拉格朗日乘子法构造新的目标函数β表示拉格朗日乘子，该目标函数对αj求导并结合约束条件可辨识出αj，每个高斯分项的权值αj、均值μj和协方差Σj的辨识结果分别为：
权值
均值
协方差
其中，为的转置；
3d、通过并行迭代执行上述步骤3b和3c直至满足迭代终止条件可辨识出混合高斯中每个高斯分项的权值αj、均值μj和协方差Σj，的概率密度函数表示为根据均值和协方差的定义可推导出混合高斯分布拟合的均值和协方差
<mrow>
<msubsup>
<mi>&mu;</mi>
<mi>k</mi>
<mi>&theta;</mi>
</msubsup>
<mo>=</mo>
<mo>&Integral;</mo>
<msub>
<mover>
<mi>&theta;</mi>
<mo>^</mo>
</mover>
<mi>k</mi>
</msub>
<mi>p</mi>
<mrow>
<mo>(</mo>
<msub>
<mover>
<mi>&theta;</mi>
<mo>^</mo>
</mover>
<mi>k</mi>
</msub>
<mo>)</mo>
</mrow>
<mi>d</mi>
<msub>
<mover>
<mi>&theta;</mi>
<mo>^</mo>
</mover>
<mi>k</mi>
</msub>
</mrow>
<mrow>
<msubsup>
<mi>P</mi>
<mi>k</mi>
<mi>&theta;</mi>
</msubsup>
<mo>=</mo>
<mo>&Integral;</mo>
<mrow>
<mo>(</mo>
<msub>
<mover>
<mi>&theta;</mi>
<mo>^</mo>
</mover>
<mi>k</mi>
</msub>
<mo>-</mo>
<mi>&mu;</mi>
<mo>)</mo>
</mrow>
<msup>
<mrow>
<mo>(</mo>
<msub>
<mover>
<mi>&theta;</mi>
<mo>^</mo>
</mover>
<mi>k</mi>
</msub>
<mo>-</mo>
<mi>&mu;</mi>
<mo>)</mo>
</mrow>
<mi>T</mi>
</msup>
<mi>p</mi>
<mrow>
<mo>(</mo>
<msub>
<mover>
<mi>&theta;</mi>
<mo>^</mo>
</mover>
<mi>k</mi>
</msub>
<mo>)</mo>
</mrow>
<mi>d</mi>
<msub>
<mover>
<mi>&theta;</mi>
<mo>^</mo>
</mover>
<mi>k</mi>
</msub>
<mo>.</mo>
</mrow>
5.如权利要求4所述的一种基于双层期望最大化算法的机动目标跟踪方法，其特征在于：在步骤3d中所述迭代终止条件为连续两次迭代的完全数据对数似然函数值相差不超过设定的门限值或达到设定的最大迭代次数。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610027471.0A CN105652250B (zh) | 2016-01-15 | 2016-01-15 | 一种基于双层期望最大化的机动目标跟踪技术 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610027471.0A CN105652250B (zh) | 2016-01-15 | 2016-01-15 | 一种基于双层期望最大化的机动目标跟踪技术 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN105652250A CN105652250A (zh) | 2016-06-08 |
| CN105652250B true CN105652250B (zh) | 2018-01-05 |
Family
ID=56486737
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201610027471.0A Active CN105652250B (zh) | 2016-01-15 | 2016-01-15 | 一种基于双层期望最大化的机动目标跟踪技术 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN105652250B (zh) |
Families Citing this family (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106931966B (zh) * | 2017-02-24 | 2019-07-26 | 西北工业大学 | 一种基于泰勒高阶余项拟合的组合导航方法 |
| CN107390199B (zh) * | 2017-09-20 | 2019-06-18 | 哈尔滨工业大学（威海） | 一种雷达机动目标跟踪波形设计方法 |
| CN107730537B (zh) * | 2017-09-29 | 2020-07-07 | 桂林电子科技大学 | 基于箱粒子概率假设密度滤波的弱目标检测与跟踪方法 |
| CN108226887B (zh) * | 2018-01-23 | 2021-06-01 | 哈尔滨工程大学 | 一种观测量短暂丢失情况下的水面目标救援状态估计方法 |
| EP3575827A1 (en) * | 2018-06-01 | 2019-12-04 | Aptiv Technologies Limited | Method for robust estimation of the velocity of a target using a host vehicle |
| WO2020113353A1 (zh) * | 2018-12-03 | 2020-06-11 | 深圳大学 | 一种机动目标的跟踪方法及系统 |
| CN110208792B (zh) * | 2019-06-26 | 2020-02-11 | 哈尔滨工业大学 | 同时估计目标状态和轨迹参数的任意直线约束跟踪方法 |
| CN110646790B (zh) * | 2019-09-30 | 2021-08-03 | 电子科技大学 | 一种用于雷达组网失序量测集中式融合的目标跟踪方法 |
| CN113534133B (zh) * | 2021-07-21 | 2022-10-21 | 西安电子科技大学 | 基于期望最大迭代算法的雷达多目标联合检测及跟踪方法 |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104794735A (zh) * | 2015-04-02 | 2015-07-22 | 西安电子科技大学 | 基于变分贝叶斯期望最大化的扩展目标跟踪方法 |
-
2016
- 2016-01-15 CN CN201610027471.0A patent/CN105652250B/zh active Active
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104794735A (zh) * | 2015-04-02 | 2015-07-22 | 西安电子科技大学 | 基于变分贝叶斯期望最大化的扩展目标跟踪方法 |
Non-Patent Citations (4)
| Title |
|---|
| EM框架下实现被动测距的状态估计和参数学习;王万平 等;《红外与激光工程》;20120731;第41卷(第7期);1708-1713 * |
| Solutions to the Linear Smoothing Problem;H. E. Rauch;《IEEE TRANSACTIONS ON AUTOMATIC CONTROL》;19631231;371-372 * |
| Using EM To Estimate A Probablity Density With A Mixture Of Gaussians;Aaron A. D Souza;《Technical note》;19991231;1-11 * |
| 未知杂波环境下的多目标跟踪算法;连峰 等;《自动化学报》;20090731;第35卷(第7期);851-858 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN105652250A (zh) | 2016-06-08 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN105652250B (zh) | 一种基于双层期望最大化的机动目标跟踪技术 | |
| CN110223348A (zh) | 基于rgb-d相机的机器人场景自适应位姿估计方法 | |
| CN106153048A (zh) | 一种基于多传感器的机器人室内定位及制图系统 | |
| CN109459019A (zh) | 一种基于级联自适应鲁棒联邦滤波的车载导航计算方法 | |
| CN108490433B (zh) | 基于序贯滤波的空时偏差联合估计与补偿方法及系统 | |
| CN104713560B (zh) | 基于期望最大化的多源测距传感器空间配准方法 | |
| CN107396322A (zh) | 基于路径匹配与编码译码循环神经网络的室内定位方法 | |
| CN104793201B (zh) | 一种跟踪临近空间高超声速目标的修正变结构网格交互多模型滤波方法 | |
| CN106599368A (zh) | 基于改进粒子提议分布和自适应粒子重采样的FastSLAM方法 | |
| CN104331623B (zh) | 一种机动策略自适应的目标跟踪信息滤波方法 | |
| CN103940433B (zh) | 一种基于改进的自适应平方根ukf算法的卫星姿态确定方法 | |
| CN104794735B (zh) | 基于变分贝叶斯期望最大化的扩展目标跟踪方法 | |
| CN110503071A (zh) | 基于变分贝叶斯标签多伯努利叠加模型的多目标跟踪方法 | |
| CN108645415A (zh) | 一种船舶航迹预测方法 | |
| CN108896986A (zh) | 一种基于预测值的量测转换序贯滤波机动目标跟踪方法 | |
| CN105719314A (zh) | 基于单应性估计和扩展卡尔曼滤波的无人机位置估计方法 | |
| CN110514567B (zh) | 基于信息熵的气体源搜索方法 | |
| CN111595334A (zh) | 基于视觉点线特征与imu紧耦合的室内自主定位方法 | |
| CN107727097A (zh) | 基于机载分布式位置姿态测量系统的信息融合方法和装置 | |
| CN103792515B (zh) | 一种异平台2维雷达与红外传感器量测数据合成方法 | |
| CN103296995B (zh) | 任意维高阶(≥4阶)无味变换与无味卡尔曼滤波方法 | |
| CN103839280B (zh) | 一种基于视觉信息的人体姿态跟踪方法 | |
| CN111812978A (zh) | 一种多无人机协作slam方法与系统 | |
| CN111711432B (zh) | 一种基于ukf和pf混合滤波的目标跟踪算法 | |
| Knuth et al. | Distributed collaborative 3D pose estimation of robots from heterogeneous relative measurements: an optimization on manifold approach |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |