基于格拉斯曼流形的状态跟踪方法研究*

赵 成，苏岭东，马祥林

（1.杭州邦友安派智能科技有限公司，浙江 杭州 310012；2.国网徐州供电公司，江苏 徐州 221000；3.常州致科自动化科技有限公司，江苏 常州 213001）

0 引言

近年来，格拉斯曼流形在医疗科学、天文学和机器人控制等领域得到了广泛运用，特别是在雷达信号处理过程中，m个基点观察k个目标，格拉斯曼流形Gk，m-k可以作为一个新的目标群采样空间[1]。而对于观测过程中难以避免的噪声干扰问题，国外学者Y.Chikuse 完成了很多基础性工作并进行了有益探索。文献[2]中提出一种将状态空间模型置于格拉斯曼流行上，在后验跟踪模型的基础上递推估计的方法。但是，该方法将观测模型和噪声都限制在格拉斯曼流形上，使得估计值与实际值存在较大的偏差。本文中尝试利用粒子滤波良好的滤波性能，提出了一种基于格拉斯曼流形的粒子滤波（Grass-Mann Manifolds-Paticle Filter，GM-PF）算法来估计流形上的隐马尔科夫过程。

安排如下：首先，简要介绍必要的背景知识，包括流形上的粒子滤波的概念；其次，提出详细的动态模型和估计方法；再次，进行仿真并分析和讨论；最后，简要总结和讨论今后的研究方向。

1 背景知识

1.1 格拉斯曼流形及其上的概率分布

格拉斯曼流形Gk，m-k可以看作是所有点为k维子空间v的空间集合，即在Rm上包含原点的k维超平面。每个在Gk，m-k上的k维子空间v都对应一个唯一的m×m正交投影矩阵P，幂为k。假如存在 在Vk，m上k列 的m×k矩 阵Y，则YY´=P；Pk，m-k表示全部幂为k的m×m正交映射矩阵的集合。因此，在流形Pk，m-k上进行的统计分析等于在格拉斯曼流形Gk，m-k上进行分析。Pk，m-k的随机矩阵P分布的密度方程为：

pFq为超几何函数矩阵。在文献[3]中的分布上略加修改即可得到Pk，m-k的矩阵Langevin 分布，记为L(p)(m，k;B)。为了简化分析，这里假设B的秩为k，则B的谱分解为：

其中：

其中在B上限制如下条件：trB=b被固定或者rankB=b＜m。

对于分布模式ΓΓ´，当B为半正定时，λi可能影响分布的形式。Chikuse and Watson 在文献[4]中定义并讨论了式（1）的分布在矩阵参数为多项式或超几何函数的情况下的一般性。

为了更好地理解格拉斯曼流形上的概率分布，这里引入不变测度的概念。定义矩阵X∈Rmm，其点集X∈Rmm，关于矩阵变量的实值函数f(x)可视为矩阵元上面的多重积分。M为一黎曼流形，而Ξ是作用在M上的可测子集。那么，如果对于每一个可测量S⊂M，有μ(S)=μ(TS)=μ(ST)，T∈Ξ，μ是关于可测子集的不变测度。若令μ(M)=1，则μ是单位不变测度。

与通常定义在Rn上的Lebesgue 测度dX不同，Gk，m-k上的概率分布被定义表示为[dX]。格拉斯曼流形上的概率分布的数学表达形式为，其中为关于不变测度微分形式[dX]的矩阵变量的密度函数[5]。

1.2 基于格拉斯曼流形的粒子滤波

在经典的欧式空间信号处理过程中，一个典型的滤波框架由离散时间的隐马尔科夫过程xt和带噪声观测的yt=g(xt)+nt构成，其中g为映射，nt通常设为高斯白噪声。概括来说，就是通过在t时刻关于未知状态xt的一些条件函数和观测变量y1，…，yt来估计E(f(xt)|y1，…，yt)。

如果选择噪声方差为高斯密度函数，就是人们熟知的卡尔曼滤波器（Kalman Filter），映射g为线性映射，p(xt|xt-1)为马尔科夫过程[6]。在这种情况下，大量统计得到的后验分布可以接近真实分布。但是，出现非线性非高斯噪声时，常规的递归滤波方法难以处理，使得需寻求一种近似计算的方法。一种常用的近似方法将局部线性化，即扩展卡尔曼滤波（Extended Kalman Filter，EKF）。另一种基于序贯重要性采样的近似方法，即粒子滤波（Particle Filter）。在重要性分布上抽取N个样本（即“粒子”），根据密度p(xt|xt-1)和p(yt|yt)递归更新归一化重要性权值，则后验期望E(f(xt)|y1，…，yt)的估计可由得到。显然，重要性函数的选择直接决定变量的估计和滤波器的性能[7-9]。

然而，在格拉斯曼流形框架下粒子滤波方法同样通过计算加权均值样本来接近期望的函数E(f)。设点位于格拉斯曼流形Gk，m-k上，目前基本上有两种方法计算他们的均值。外蕴均值法将流形中嵌入一个欧式空间，在欧式空间上计算算术平均值，然后将结果点映射回流形上。但是，假如最接近的点不唯一，外蕴均值法会产生多个均值或者均值为空。内蕴均值法充分利用流形的几何结构[10]，可以有效解决由于指数映射导致的低效问题，但同样存在缺陷，需要解决相同尺度下最近点的最小化问题（如Karcher 均值法和质心法[11]）。本文采用外蕴均值法，内蕴均值计算法是今后工作的重点。

2 模型选择

2.1 矩阵的概率分布

一个m×k随机矩阵X的正态分布如下：

M为均值矩阵，Ω为协方差矩阵，上的概率密度。dX为Lebesgue 测度，tr(·)表示迹。

在格拉斯曼流形上一个m×k随机矩阵X和参数矩阵B的von Mises-Fisher 分布（也被称为Langebin 分布）为[12]：

其中C为标准化常数，[dX]为在格拉斯曼流形上的单位不变测度，p(X)[dX]=1[dX]是相应的均匀分布。

矩阵的von Mises-Fisher 分布与欧几里得空间的正态分布密切相关。假如一个始于矩阵正态密度的随机矩阵X，参数，限制条件X´X=I，易得到由此产生的密度Gmpd(X;Ω-1G)。所以，将选择不同的噪声能量对比GM-PF 算法和GM-PMTA 算法比较优劣。

2.2 GM-PF 动态模型

为了解决在统计信号处理过程中的矩阵变量的估计和跟踪问题，这里介绍一种基于格拉斯曼流形的动态空间模型。假定离散时间的隐马尔科夫过程Xt位于格拉斯曼流形Gk，m-k上，则矩阵的von Mises-Fisher 分布为

设观测Yt由Yt=Xt+Nt得到，Nt为均值为零的高斯矩阵，参数，所以观测模型为：

2.3 后验跟踪算法的状态空间模型

与上文不同，文献[2]中的状态空间模型的观测被限制于格拉斯曼流形上，包含如下两个von Mises-Fisher 分布。Xt为在格拉斯曼流形上的隐马尔科夫过程：

观测模型为：

注意到观测误差很有可能是量测噪声干扰的结果，所以这里优先考虑2.2 节提出的动态模型。但是，文献[2]中提出的后验模型跟踪算法可以作为比较的基础，给定观测噪声来仿真两种模型，设β=Ω-1。

3 仿真实验

3.1 算 法

结合2.2 节和2.3 节提出的不同状态空间模型得到详细的滤波算法。这种算法需要在mkℝ 空间上的点Y映射到Gk，m-k上，根据文献[2]中的定义，通过对Y进行极分解得到X，即。

根据2.2 节的模型，假设有N个粒子，初始化权值，i-1，…，N。当t≥1 时粒子滤波过程如下所示[5，13-14]。

算法1：格拉斯曼流形上的粒子滤波

1.得到观测Yt～N(Yt;Xt，Ω，∑)

2.fori=1，…，N

4.是否需要重采样

根据2.3 节中的模型，用χt来表示Xt的估计。文献[2]将χ0=X0，但在这里将χ0直接随机均匀抽取。第一个观测量Y1～Gmpd(Y1;X1β)，则。当t≥2 时，算法过程如下。

算法2：后验模型跟踪算法

3.fori=1，…，收敛

3.2 仿真结果分析

将算法1 用于2.2 节中的模型，其中算法仿真时间步为125，过程噪声方差0.01，观测噪声方差0.25，每个参数都固定实验30 次取平均值。图1 为在不同粒子数的情况下算法的最小均方误差。如图1 所示，随着粒子数的加大，基于格拉斯曼流形的粒子滤波方法的误差逐渐减小，但当粒子数超过100 时，GM-PF 方法付出了巨大的计算代价，而性能提高并不明显。因此，为了兼顾滤波精度和效率，下文仿真中都将设置N=100。

图1 在G1，3上，粒子数N=25 到N=1 000 时，GM-PF 算法的MMSE

在格拉斯曼流形上的马尔科夫过程过程Xt，观测Yt通过式（7）产生。最小均方差（Minimum Mean Square Error，MMSE）状态估计根据算法1 得到。图2 和图3 显示了一个在格拉斯曼流形G1，3上典型的HMM 过程的跟踪效果，，两种算法的跟踪性能比较。图2、图3 比较可以发现，两种算法都能成功跟踪隐马尔科夫过程Xt的真值，但是后验模型跟踪算法的状态估计更为平滑，粒子滤波算法的MMSE 估计更接近真实Xt曲线。此外发现，算法1 计算所需时长是算法2 的10 倍，但是算法1可得到更低的均方误差。

为了将算法2 应用于最大后验（MAP）估计，将观测Yt映射ε(Yt)到格拉斯曼流形上产生“伪观测”。与映射前对比发现，这种简单的映射χt=ε(Yt)也能够产生一定程度的“降噪”效果。

图2 在G1，3上相同隐马尔科夫过程和噪声条件下，GM-PF 滤波算法实验结果

图3 在G1，3上相同隐马尔科夫过程和噪声条件下，GM-PMTA 方法实验结果

图4 根据2.2 节的模型在G1，3上对比几种不同方法的实验结果。首先固定过程参数α，对比各种方法在不同噪声情况下的估计误差，其中Yt和Xt之间的距离用Frobenius 范数处理，过程参数α=10I。仿真时间设为100 时间步，粒子数N=100，ε(Yt)在格拉斯曼流形上为GM-PMTA 算法产生合适的观测量。图4 中的实验结果都是同一参数30 次实验取均值的结果，Ω=σ2I，β=σ-2I。其中，GM-PF为格拉斯曼流形粒子滤波后得到的观测误差，GMPGMA 就是后验模型跟踪算法（PMTA）滤波后得到的观测误差，Projected Observed 为映射后观测到的误差，GM Observed 为在流形G1，3上未进行滤波观测到的误差。

在图4 中可以清楚看出，GM-PF 方法胜过GM-PMTA 方法，特别是在噪声增大的时候（相同的情况同样出现在2.3 节中的模型）。当噪声增大到一定程度后，GM-PMTA 方法的性能迅速下降，甚至会比仅仅将状态估计映射ε(Yt)所产生的降噪效果还要差。

图4 不同噪声能量σ2 情况下，在G1，3上不同方法的MMSE变化曲线

4 结语

本文介绍了一种矩阵变量的动态模型和与之匹配的粒子滤波算法。在低维空间上的仿真表明，相比较确定性模型跟踪算法，该方法显著提高了滤波性能，但付出了较高的计算代价。为了兼顾粒子滤波方法的精度和效率，未来研究的重点放在利用流形内部的几何结构进行内蕴均值计算，或者通过传递矩阵并行传递来替代指数映射，避免反复指数映射带来的大量计算，在保证较高精度的情况下提高效率。另外，着重研究一种在矩阵von Mises-Fisher 分布上高效的采样方法，使得这种基于格拉斯曼流形的粒子滤波方法能够解决高维空间的状态估计问题。