在实际的目标跟踪系统中常常面临非线性甚至是强非线性的估计问题^[1], 在高斯系统的假设条件下, 文献[2]基于贝叶斯滤波理论推导出非线性高斯系统滤波框架, 但是非线性系统状态后验分布的积分通常难以获得解析解. 因此, 将后验分布积分难以解析的问题转化为非线性目标函数优化问题或者近似分布的优化问题成为必然. 因此, 基于优化方法的非线性系统状态一直是一个热门研究课题, 尤其是在目标跟踪^[3-5]、智能运输^[6-7]、环境监测^[8]、 航天器导航^[9-10]等领域.

近年来, 随着非线性估计精度要求不断提高, 基于拟牛顿、高斯−牛顿和梯度上升等确定性优化方法和基于采样的随机优化方法的非线性滤波器相继被提出. 文献[11]通过不断迭代使扩展卡尔曼滤波器(Extended Kalman filter, EKF)中与非线性量测相关雅克比矩阵得到逐步修正以提高量测预测的精度. 文献[12]以最小二乘为目标函数, 采用牛顿法推导给出一种迭代的卡尔曼滤波算法. 文献[13]采用梯度法和牛顿法, 提出一种无噪声条件下的滚动时域估计器, 并证明其渐近收敛性和稳定性. 文献[14]针对非线性状态空间的状态估计及其估计误差协方差的迭代更新, 利用正则化的非线性最小二乘, 提出一种随机增量近端梯度方法. 文献[15]在多种散度最小化条件下迭代获得与后验分布更接近的近似分布, 并给出$ \alpha $散度和KL (Kullback-Leibler)散度的实现形式. 非线性状态估计优化的实质是多维后验概率分布的不断近似, 并且其近似程度不能简单地采用欧氏距离进行度量. 因此, 选取合理的度量准则将有利于提高后验概率分布的近似程度, KL散度作为分布之间“距离”的度量在后验分布近似性衡量中具有天然的优势. 同时, 从信息几何角度出发, 概率分布是统计流形上的点, 在一定条件下两概率分布之间的KL散度与作为统计流形度规的Fisher信息矩阵满足一定的数学关系.

与确定性优化不同, 变分贝叶斯(Variational Bayes, VB)框架下的随机优化一般将难以线性化的部分(例如具有非共轭特性的似然函数期望)看作“黑盒”, 并通过随机采样的方法近似. 文献[16]提出利用随机优化方法进行变分推理, 允许使用更复杂的贝叶斯模型对数据进行快速学习. 文献[17]通过引入辅助变量搭建了马尔科夫链蒙特卡罗(Markov chain Monte Carlo, MCMC)与VB方法相结合的媒介, 从而有效综合二者优势, 较好地兼顾算法实时性和估计精度指标, 并且算法对非线性非共轭系统具有较好的适应性. 在随机优化的基础上, 研究以少量采样粒子权重构建和搜索梯度方向, 从而在保证估计精度的同时提高实时性. 文献[18]提出变分贝叶斯蒙特卡罗方法(Variational Bayesian Monte Carlo, VBMC)通过基于高斯过程的主动采样并求其期望, 以实现难以积分的似然函数的有效逼近, 从而获得后验分布的非参数近似和相应的置信下界. 文献[19-20]分别针对含有噪声的似然函数模型和隐马尔科夫模型提出相应的改进算法, 为VBMC算法在非线性状态估计中的应用提供了良好的借鉴.

上述方法均在欧氏空间中考虑非线性状态估计的优化, 然而非线性状态估计优化的实质是后验分布的近似, 由于欧氏空间的局限性, 难以表征后验分布的近似程度. 在统计流形中, 概率分布被视为流形中的点. 因此, 在统计流形空间中结合信息几何理论从概率分布属性意义上考虑非线性状态估计问题具有可行性. 文献[21]提出自然梯度的优化方法以实现统计流形空间中目标函数的最速梯度下降(或上升), 并证明采用自然梯度的在线学习是Fisher有效的. 文献[21-23]给出状态空间模型下置信下界的自然梯度, 通过Fisher信息矩阵在统计流形上计算的最速下降方向, 是一种在黎曼空间中考虑后验分布近似的方法. 文献[24]基于自然梯度提出一种非线性卡尔曼滤波优化方法, 通过批处理的方式对状态与参数进行联合估计. 在文献[25]中, 自然梯度被应用于传感器网络目标跟踪系统. 文献[26]证明了针对非线性估计问题的自然梯度在克拉美罗下界意义下是渐近最优的. 结合统计流形空间和信息几何理论, 自然梯度被进一步应用于机器学习和参数辨识等领域^[27-28].

上述优化方法一般需要结合特定的优化框架实现对目标状态进行迭代更新, 因此具有迭代渐近收敛特性的优化框架设计与选择受到国内外相关领域学者的广泛关注. 文献[29]在期望最大化(Expectation maximization, EM)框架中引入一种近端点优化算法, 并对其收敛性做了研究. 在此基础上, 文献[30]采用KL散度作为近端点优化算法的约束项, 从几何近似角度接近真实后验分布, 并将其应用于基于数据驱动的机器学习中的数据分类. EM优化框架通过E步和M步交替迭代实现目标状态和未知参数的联合优化, 是变分贝叶斯推断的一种特殊形式. 在变分贝叶斯推断中, 目标状态或未知参数均被视为满足一定超参数分布的随机隐变量, 通过坐标上升方法迭代优化超参数以实现目标状态或未知参数的联合高精度估计. 文献[31-32]给出VB框架下非线性状态估计的实现. 文献[33]借助变分贝叶斯方法能够将难以获得解析解的非线性后验分布转化为迭代优化问题的优势, 结合正则化方法以KL散度为惩罚项, 构建以目标状态估计及其估计误差协方差为超参数的优化函数, 继而通过线性化的方式获得目标状态估计及其估计误差协方差的迭代优化形式.

非线性状态后验分布的近似过程可以看作近似分布(如变分分布)在统计流形空间中, 沿着选定的迭代搜索方向, 向后验分布移动的过程, 如图1所示, 变分分布$ q({{{\boldsymbol{x}}}}_{k}|{\boldsymbol \psi}_{k}) $沿橙色曲线搜索方向, 向状态后验分布 $ p({{{\boldsymbol{x}}}}_{k}|{{{\boldsymbol{z}}}}_{k}) $ 移动, 其中$ {\cal{D}}_{\rm{KL}}[q({{{\boldsymbol{x}}}}_{k}|{\boldsymbol \psi}_{k})||p({{{\boldsymbol{x}}}}_{k}|{{{\boldsymbol{z}}}}_{k})] $表示变分分布与后验分布之间的KL散度. 变分贝叶斯推断利用参数化变分分布迭代逼近未知分布的思想, 使其在处理后验分布难以获得解析解的问题上有显著的优势^[33-34], 特别是针对非线性后验概率密度近似问题具有较好的近似精度^{[3, 11, 35]}. 因此, 在变分贝叶斯框架下, 结合信息几何中Fisher信息矩阵在变量邻域内与KL散度的二阶微分近似关系, 推导置信下界(Evidence lower bound, ELBO)最大化条件下的自然梯度; 进而, 提出一种基于自然梯度的非线性迭代变分贝叶斯估计算法(VBKF-NG). 该方法将非线性状态后验分布难以积分的问题转化为置信下界的最大化问题, 具有以下特点和优势:

图 1  变量分布近似过程中的KL散度示意图

Fig. 1  The KL divergence of the distribution approximation of a variable

下载: 全尺寸图片 幻灯片

1)在统计流形空间中, 结合信息几何理论从概率分布属性意义上考虑非线性状态估计问题;

2)借助于Fisher信息矩阵推导置信下界的自然梯度, 实现在统计流形空间中参数化变分分布对后验分布的“紧密”近似;

3)结合最优估计理论, 推导状态估计及其误差协方差的自然梯度, 并以此为切入点获得状态估计及其误差协方差的迭代更新解析解.

为便于查找文中符号含义, 表1 给出了相应的符号说明.

表 1  文中变量和符号含义

Table 1  The meaning of variables and symbols

变量	符号含义
$ {\boldsymbol x}_k $	$ k $时刻目标状态真实值
$ {\boldsymbol x}_{k|k} $	$ k $时刻目标状态估计值
$ {\boldsymbol P}_{k|k} $	$ k $时刻目标状态估计误差协方差
$ {\boldsymbol z}_k $	传感器在$ k $时刻的量测值
$ {\boldsymbol\omega}_k $	$ k $时刻的系统噪声
$ {\boldsymbol\upsilon}_k $	$ k $时刻的量测噪声
$ {\boldsymbol Q}_k $	$ k $时刻系统噪声方差
$ {\boldsymbol R}_k $	$ k $时刻量测噪声方差
$ d_x $	目标状态向量的维数
$ d_z $	量测向量的维数
$ {\boldsymbol F}_{k|k-1} $	$ k-1 $时刻到$ k $时刻的状态转移矩阵
$ {\boldsymbol H}_{k} $	$ k $时刻量测矩阵
$ {\boldsymbol\psi}_{k} $	变分分布参数
$ p\left({\boldsymbol x}_k|{\boldsymbol z}_{k}\right) $	$ k $时刻目标状态后验分布
$ q\left({\boldsymbol x}_k|{\boldsymbol\psi}_{k}\right) $	以$ {\boldsymbol\psi}_{k} $为参数的变分分布
$ \mathcal{L}\left({\boldsymbol\psi}_{k}\right) $	以$ {\boldsymbol\psi}_{k} $为变分分布参数的置信下界
${\cal{D}}\left(q\left({\boldsymbol x}_k|{\boldsymbol\psi}_{k}\right)|| p\left({\boldsymbol x}_k|{\boldsymbol z}_{k}\right)\right)$	变分分布$ q\left({\boldsymbol x}_k|{\boldsymbol\psi}_{k}\right) $与状态后验分布$ p\left({\boldsymbol x}_k|{\boldsymbol z}_{k}\right) $的KL散度
$ {\boldsymbol J}_{{\boldsymbol\psi}_k} $	以$ {{\boldsymbol\psi}_k} $为参数的 Fisher 信息矩阵
$ \mathcal{M} $	流形空间
$ \mathcal{S} $	流形空间中的概率分布集合
$ \mathcal{F} $	流形空间中的平滑映射函数
${{\boldsymbol v} }_{OP}$	流形空间中的$ O $点处指向$ P $点的切向量
$|{{\boldsymbol v} }_{OP}|$	切向量${{\boldsymbol v} }_{OP}$的模

下载: 导出CSV 

| 显示表格

1.   预备知识

1.1   非线性系统状态估计

考虑一般非线性动态系统方程及其量测方程

其中, ${{{\boldsymbol{x}}}}_{k} \in{{\bf{R}}^{d_x}}$表示$ k $时刻系统的真实状态(比如位置、速度等其他相关量), ${{{\boldsymbol{z}}}}_{k} \in{{\bf{R}}^{d_z}}$表示$ k $时刻系统状态的量测(例如径向距、角度和多普勒频移等), $ d_x $和$ d_z $分别表示系统状态$ {{{\boldsymbol{x}}}}_{k} $和量测$ {{{\boldsymbol{z}}}}_k $的维数. $ {{{\boldsymbol{f}}}}_{k|k-1}\left(\cdot\right) $和$ {{{\boldsymbol{h}}}}_{k}\left(\cdot\right) $分别表示状态转移函数和量测函数, 系统噪声$ {\boldsymbol\omega}_{k} $和量测噪声$ {\boldsymbol\upsilon}_{k} $为满足相互独立的零均值高斯白噪声, 其方差分别表示为$ {{{\boldsymbol{Q}}}}_{k} $和$ {{{\boldsymbol{R}}}}_{k} $. 图2给出非线性动态系统时序状态转移和量测的示意图.

图 2  非线性动态系统状态转移和量测的示意图

Fig. 2  The state transition and measurement in a nonlinear dynamic system

下载: 全尺寸图片 幻灯片

在经典贝叶斯估计框架下, 状态估计过程是在最小均方误差(Minimum mean square error, MMSE)意义下从给定的含有随机误差的量测信息中推断出当前时刻被估计系统的状态信息. 系统状态估计$ {{{\boldsymbol{x}}}}_{k|k} $及其误差协方差$ {{{\boldsymbol{P}}}}_{k|k} $分别表示为

其中, $ {\rm{E}}_p[\cdot] := {\rm{E}}_{p({{{\boldsymbol{x}}}}_{k} \mid {{{\boldsymbol{z}}}}_k)}[\cdot] $, $({{\boldsymbol{x}}}_k- {{{\boldsymbol{x}}}}_{k|k})(\cdot)^{\rm T} = ({{{\boldsymbol{x}}}}_k\;- {{{\boldsymbol{x}}}}_{k|k})({{{\boldsymbol{x}}}}_k- {{{\boldsymbol{x}}}}_{k|k})^{\rm T}$, 后文类似缩写不再赘述. $ p({{{\boldsymbol{x}}}}_{k}|{\boldsymbol{z}}_k) $为状态的后验分布, 并表示为

在 $ {{{\boldsymbol{x}}}}_k $ 服从高斯分布的条件下, 通过利用查普曼−科莫高洛夫(Chapman-Kolmogorov)方程, 式(5)中的状态预测的概率密度函数$ p\left({{{\boldsymbol{x}}}}_k|{\boldsymbol{z}}_{k-1}\right) $表示为

其中, $ {{{\boldsymbol{x}}}}_{k|k-1} $和$ {{{\boldsymbol{P}}}}_{k|k-1} $分别表示$ k $时刻的状态预测和预测协方差. 量测似然函数为

若式(1)和式(2)中的被估计系统满足高斯线性假设, 则可直接采用卡尔曼滤波器进行状态估计, 包括状态预测和量测更新步骤.

状态预测:

量测更新:

其中, $ {{{\boldsymbol{F}}}}_{k|k-1} $和$ {{{\boldsymbol{H}}}}_{k} $分别表示线性系统中的状态转移函数和量测函数.

在非线性系统中, 式(5)分母中的联合概率密度函数的边缘分布$ \int p\left({{{\boldsymbol{x}}}}_{k},{{{\boldsymbol{z}}}}_k\right){\rm{d}}{{{\boldsymbol{x}}}}_{k} $一般难以求解, 导致无法获得后验分布的解析解, 而只能通过优化方法求其近似解. 因此, 采用何种方法, 如何度量后验分布的近似程度有待进行深入研究, 以提高非线性状态后验分布的近似程度.

欧氏距离不能完全表征随机变量的分布近似性, 特别是多维变量分布近似的情况. 如图3所示, 在欧氏空间中, 单变量高斯分布$ q_1\sim {\rm{N}}(-2,\;0.75^2) $和$ q_2\sim {\rm{N}}(2,\;0.75^2) $ 之间的欧氏距离与$ q_3\sim {\rm{N}}(-2, \;2^2) $和$ q_4\sim {\rm{N}}(2,\;2^2) $之间的欧氏距离相同, 即 $ d_{12}= d_{34} $. 但是由于具有不同的方差, 计算其KL散度分别为${\cal{D}}_{\rm KL}[q_1\|q_2] = 10.667$ 和 $ {\cal{D}}_{\rm KL}[q_3\|q_4] = 2 $. 图4表示多变量高斯分布的KL散度差异, 其中, 分布分别表示为$q_5 \sim{\rm{N}}\left(\left[\begin{aligned} 2 \\ 2 \end{aligned}\right], \left[\begin{aligned} 0.2^2 \;\; 0\;\\ 0 \;\;\;\; 1^2 \end{aligned}\right]\right)$, $q_6 \sim{\rm{N}}\left(\left[\begin{aligned} 0 \\ 2 \end{aligned}\right],\right. $ $\left.\left[\begin{aligned} 0.5^2 \;\; 0\;\;\;\\ 0 \;\;\;\; 0.2^2 \end{aligned}\right]\right)$ 和 $q_7 \sim{\rm{N}}\left(\left[\begin{aligned} 0 \\ 2 \end{aligned}\right], \left[\begin{aligned} 0.5^2 \;\; 0\;\\ 0 \;\;\;\; 1^2 \end{aligned}\right]\right)$, 分布$ q_5 $ 和$ q_6 $之间的欧氏距离与$ q_5 $和$ q_7 $之间的欧氏距离均为$ 2 $, 但其KL散度分别为$ {\cal{D}}_{\rm KL}[q_5\|q_6] = 8.496 $, ${\cal{D}}_{\rm KL} [q_5\|q_7] = 9.625\;7.$ 由此可见, 欧氏距离未考虑到总体差异(如方差)对概率分布之间近似性度量的影响, 而KL散度作为流形空间中概率分布相似性的度量更具有一般性.

图 3  单变量高斯分布的欧氏距离示意图

Fig. 3  The Euclidean distance for univariate Gaussian distributions

下载: 全尺寸图片 幻灯片

图 4  多变量高斯分布的欧氏距离示意图

Fig. 4  The Euclidean distance for multivariate Gaussian distributions

下载: 全尺寸图片 幻灯片

从概率论角度, 非线性状态估计一般认为是后验分布近似的过程, 即在给定量测的条件下通过递推使变量估计值逐渐接近其真实值的过程. 在信息几何理论中概率分布被视为统计流形空间中的点, 自然地, 可将变量空间中后验分布的近似转化为统计流形空间中点的近似. 在已知量测序列条件下, 图5表示在不同空间中通过各类优化方法对非线性状态进行迭代估计的含义. 图5(a)表示变量在欧氏空间中通过迭代优化促使估计值$ {\boldsymbol{x}}_{k|k} $逐渐逼近未知变量真实值$ {\boldsymbol{x}}_k $; 图5(b)表示在分布空间中以$ {\boldsymbol\psi}_{k|k} $为参数的变分分布$ q({\boldsymbol{x}}_{k}|{\boldsymbol\psi}_k) $通过迭代优化不断逼近未知变量的后验分布$ p({\boldsymbol{x}}_{k}|{\boldsymbol{z}}_{k}) $; 图5(c)表示在统计流形空间中变分分布$ q({\boldsymbol{x}}_{k}|{\boldsymbol\psi}_k) $作为该空间中的一点向表示真实后验分布$ p({\boldsymbol{x}}_{k}|{\boldsymbol{z}}_{k}) $的另一点移动.

图 5  非线性状态估计在不同空间中的含义示意图

Fig. 5  The meaning of nonlinear state estimation in different spaces

下载: 全尺寸图片 幻灯片

1.2   自然梯度

设以$ {\boldsymbol{w}} $为变量的目标函数表示为$ {{\cal{L}}({\boldsymbol{w}})} $, 梯度下降方法是选取使$ {{\cal{L}}({\boldsymbol{w}})} $减小的方向作为搜索方向对变量$ {\boldsymbol{w}} $进行迭代寻优, 进而使$ {{\cal{L}}({\boldsymbol{w}})} $达到最小值. 微分向量的平方长度表示为

其中, $ {\rm {d}}w_i $表示$ {\rm {d}}{\boldsymbol{w}} $的元素, $ n $表示变量$ {\boldsymbol{w}} $的维数. 上述过程的合理性须满足一个基本的事实: 即目标函数和变量及其变化量均在同一欧氏空间的正交坐标系中被衡量.

在流形空间中, 作为分布的点向另一点移动的距离 $ |{\rm {d}}{\boldsymbol{w}}|^2 $表示为

其中, $g_{i,j}\left({\boldsymbol{w}}\right)=\delta_{ij}= \left\{ \begin{aligned} 1, \; i=j\\ 0,\; i\neq j \end{aligned} \right.$, 令 $ {{{\boldsymbol{G}}}} = \left(g_{i,j}\left({\boldsymbol{w}}\right)\right) $, 即黎曼张量. 在欧氏空间中, $ {\boldsymbol{G}} $退化为相应维数的单位阵.

令$ {\rm {d}}{\boldsymbol{w}} = \varepsilon{\boldsymbol{a}} $, $ \varepsilon $表示一个较小常数, 则目标函数$ {{\cal{L}}({\boldsymbol{w}})} $在$ {\boldsymbol{w}} $处的泰勒展开表示为

假设向量长度约束 $ |{\boldsymbol{a}}|^2 = \sum g_{i,j}{\boldsymbol{a}}_i{\boldsymbol{a}}_j=1 $, 采用拉格朗日法对线性化后的目标函数求极值, 可得

其中, $ \lambda $为拉格朗日乘数. 进而可得目标函数在 $ {\boldsymbol{w}} $的偏导为

由于$ \lambda $为常数, 自然梯度^[21]即表示为

2.   基于自然梯度的迭代变分贝叶斯滤波

变分贝叶斯方法是一种采用简单分布近似未知变量真实后验分布的推断方法^[36], 通常假设未知变量之间相互独立. 变分贝叶斯推断主要有两个目的^[1]: 一是逼近边缘似然函数$ p\left({{\boldsymbol{z}}_k}|{{\boldsymbol{x}}_k}\right) $以实现模型的选择, 即变分超参数的优化; 二是逼近包含所有未知变量的后验分布$ p\left({{\boldsymbol{x}}_k}|{{\boldsymbol{z}}_k}\right) $以实现未知变量的估计. 假设未知变量的全贝叶斯模型中所有参数的先验分布已知, 变分贝叶斯推断的目的是寻找后验分布的一个变分分布使其与真实后验分布之间的KL散度最小. 因此, 变分贝叶斯作为一种优化方法能够通过构建变分分布和闭环迭代实现非线性系统状态后验分布的“紧密”近似, 进而获得相应的高精度状态估计值.

给定系统状态$ {\boldsymbol{x}}_k $的量测$ {\boldsymbol{z}}_k $, 其联合概率密度函数和边缘概率密度函数分别为$ p({\boldsymbol{x}}_k,{\boldsymbol{z}}_k) $和$ p({\boldsymbol{z}}_k) $, 根据概率论知识可知

其中, $ q\left({\boldsymbol{x}}_k|{\boldsymbol{\psi}}_k\right) $表示以$ {\boldsymbol{\psi}}_k $为参数的变分分布. 在变分贝叶斯框架下, 以变分分布$ q\left({\boldsymbol{x}}_k|{\boldsymbol{\psi}}_k\right) $作为桥梁, 上式可转化为

${{{\cal{L}}}}\left({\boldsymbol{\psi}}_k\right)$表示具有单调递增特性的变分置信下界, ${\cal{D}}_{\rm {KL}}\left[q({\boldsymbol{x}}_k|{\boldsymbol{\psi}}_k)\|p({\boldsymbol{x}}_k|{\boldsymbol{z}}_k)\right]$ 表示变分分布 $ q\left({\boldsymbol{x}}_k|{\boldsymbol{\psi}}_k\right) $ 与后验分布$ p({\boldsymbol{x}}_k|{\boldsymbol{z}}_k) $之间的KL散度, 其表达式分别为

在变分优化的过程中, 以KL散度作为变分分布与后验分布的近似程度的度量, 当变分分布越接近真实后验分布, KL散度的值越小; 反之, KL散度的值越大. 当KL散度为$ 0 $时, 此时的变分分布$ q({{\boldsymbol{x}}}_k|{\boldsymbol{\psi}}_k) $等于真实后验分布$ p({\boldsymbol{x}}_k|{\boldsymbol{z}}_k) $. 图1为变分分布$ q({{\boldsymbol{x}}}_k|{\boldsymbol{\psi}}_k) $在向后验分布$ p({\boldsymbol{x}}_k|{\boldsymbol{z}}_k) $移动(即近似)的过程中$ q({{\boldsymbol{x}}}_k|{\boldsymbol{\psi}}_k) $ (即估计的变分分布$ q({{\boldsymbol{x}}}_k|{{\boldsymbol{x}}}_{k|k},{{\boldsymbol{P}}}_{k|k})) $和$ p({\boldsymbol{x}}_k|{\boldsymbol{z}}_k) $ 之间的KL散度的示意图.

在理论上, KL散度最小化是实现高斯系统非线性状态后验分布$ p({\boldsymbol{x}}_k|{\boldsymbol{z}}_k) $近似的最佳选择. 然而, 在实际应用过程中由于上述KL散度需计算对数边缘分布$ \log p({\boldsymbol{z}}_k) $, 结合式(23), KL散度进一步转化为

其中, 上式中所有期望均是在变分分布$ q({\boldsymbol{x}}_k|{\boldsymbol{\psi}}_k) $下的期望. 从上式可以看出KL散度依赖于$ \log p({{{\boldsymbol{z}}}}_k) $. 由式(5)可知, 非线性系统中的$ \log p({\boldsymbol{z}}_k) $难以计算. 因此, 以上述KL散度为目标函数通常难以实现优化目的. 由于在给定量测条件下$ \log p\left({\boldsymbol{z}}_k\right) $为常数, 结合式(21), KL散度的最小化等价于变分ELBO的最大化问题, 即

对于任意变分分布$ q\left({\boldsymbol{x}}_k|{\boldsymbol{\psi}}_k\right) $, 上述不等式均成立, 并且可以看出ELBO是边缘分布的保守估计下界, 它包含未知变量所有的量测信息. 相应的未知变量的优化可表示为

通过最大化ELBO能够将非线性状态后验难以积分问题转化为优化问题, 同时ELBO单调递增的特点有利于通过迭代的方式逼近最优解, 以获取非线性状态真实后验分布的最优近似分布(即变分分布), 从而实现状态的高精度估计.

2.1   置信下界的自然梯度

式(25)建立了以ELBO最大为条件的目标函数, 结合适当的优化方法能够获得ELBO取最大值时变分分布参数$ {\boldsymbol{\psi}}_k $的最优解. 自然梯度通过信息几何的方法寻找目标函数在统计流形空间中的最速上升/下降方向^{[21, 37]}. 结合式(22)中ELBO的定义, 下面从理论上推导给出ELBO关于变分参数$ {\boldsymbol \psi}_{k} $的自然梯度.

根据式(22), 变分ELBO可以进一步分解为

其中, $ {\rm{E}}_{q}\left[\log p\left({\boldsymbol{z}}_{k}|{\boldsymbol{x}}_{k}\right)\right] $ 和 $ {\cal{D}}_{\text{KL}} \left[q\left({\boldsymbol{x}}_{k}|{\boldsymbol{\psi}}_{k}\right)\| p\left({\boldsymbol{x}}_{k}\right)\right] $ 分别表示对数似然期望和变分分布与先验信息之间的KL散度. 前者通过量测信息作用于变分分布参数$ {\boldsymbol{\psi}}_{k} $, 其原理类似于最大似然估计, 后者考虑以KL散度为度量先验分布$ p({\boldsymbol{x}}_{k}) $对变分参数$ {\boldsymbol{\psi}}_{k} $的影响, ELBO从贝叶斯角度平衡量测似然和先验信息, 以实现变分分布对后验分布的近似.

假设第$ i $次迭代的估计结果是第$ i+1 $次迭代的先验, 即上式中的先验信息$ p\left({\boldsymbol{x}}_{k}\right) $可以近似为$ p\left({\boldsymbol{x}}_{k}\right)\approx q({\boldsymbol{x}}_{k}|{\boldsymbol{\psi}}_{k}^i) $, 相应的KL散度为${\cal{D}}_{\rm {KL}}[q({\boldsymbol{x}}_k|{\boldsymbol{\psi}}_k) \|q({\boldsymbol{x}}_k|{\boldsymbol{\psi}}_k^{i})].$ 那么变分参数$ {\boldsymbol{\psi}}_{k} $的最优估计值可以表示为

对式(28)等号右边进行线性化并计算其关于$ {\boldsymbol \psi}_{k} $的偏导数, 通过令其偏导数为$ 0 $, 可获取ELBO的自然梯度和相应的最优$ {\boldsymbol \psi}_k^* $.

引理 1. 假设系统状态演化满足高斯分布特性, 定义变分参数 $ {\boldsymbol \psi}_k $ 的差分$ \Delta {\boldsymbol \psi}_k = {\boldsymbol \psi}_k- {\boldsymbol \psi}_k^i\to 0 $, 若$ \rm KL $散度$ {\cal{D}}_{\text{KL}}\left[q\left({{{\boldsymbol{x}}}}_{k}|{\boldsymbol \psi}_{k}\right)\| q\left({{{\boldsymbol{x}}}}_{k}|{\boldsymbol \psi}_{k}^{i}\right)\right] $在$ {\boldsymbol \psi}_k $点具有二阶偏导, 则有以下近似关系成立,

其中, $ {{{\boldsymbol{J}}}}_{{\boldsymbol \psi}_{k}} $为置信下界关于参数$ {\boldsymbol \psi}_{k} $的Fisher信息矩阵.

证明. 根据系统中噪声服从高斯分布条件下$ \rm KL $散度的定义, 在变分迭代过程中$ \rm KL $散度表示为

首先, 根据泰勒展开将上式中的$ \log {q\left({\boldsymbol{x}}_k|{\boldsymbol{\psi}}_k\right)} $在点$ {\boldsymbol{\psi}}_k={\boldsymbol{\psi}}_k^{i} $处线性化, 其表达式为

其中, $\nabla_{{\boldsymbol{\psi}}_k}\left[\log q\left({{{\boldsymbol{x}}}}_{k}|{\boldsymbol{\psi}}_k\right)\right]$ 表示$\log q\left({{{\boldsymbol{x}}}}_{k}|{\boldsymbol{\psi}}_k\right)$ 在 $ {\boldsymbol{\psi}}_k= {\boldsymbol \psi}_{k}^{i} $处的梯度, ${{{\rm{O}}}}({\Delta{{\boldsymbol \psi}_{k}}}^{3})$表示$ {\Delta{{\boldsymbol \psi}_{k}}}^{3} $的高阶无穷小. 进而, 式(30)中的$ \rm KL $散度表示为

由于变分分布$ q\left({\boldsymbol{x}}_k|{\boldsymbol{\psi}}_k\right) $ 关于$ {\boldsymbol{x}}_k $ 的积分与$ {\boldsymbol{\psi}}_k $ 无关, 式(32)中的第一项可进一步表示为

式(32)中的第二项转化为

由于

并且根据信息矩阵Fisher的定义式

式(34)可进一步表示为

因此, KL散度$ {\cal{D}}_{\text{KL}} \bigl[q\left({\boldsymbol{x}}_k|{\boldsymbol{\psi}}_k\right) \|q\left({\boldsymbol{x}}_k|{\boldsymbol{\psi}}_k^{i}\right)\bigr] $的线性化表达式为

计算 $ \rm KL $ 散度 $ {\cal{D}}_{\text{KL}} \bigl[q\left({\boldsymbol{x}}_k|{\boldsymbol{\psi}}_k\right) \|q\left({\boldsymbol{x}}_k|{\boldsymbol{\psi}}_k^{i}\right)\bigr] $ 关于$ \Delta {\boldsymbol{\psi}}_k $的二阶偏导, 可得到关于$ {\boldsymbol{\psi}}_k $的Fisher信息矩阵的近似表达式

□

注 1. Fisher信息矩阵具有半正定性, 在统计流形空间中也被称为黎曼度量. 引理1建立了变分置信下界与统计流形空间中黎曼度量之间的数学关系. 同时, 在物理意义上Fisher信息矩阵可视为对数似然的曲率.

1)当Fisher信息矩阵较小时, 作为目标函数的ELBO变化较缓, 在最大值附近存在较多与其相接近的次优值, 导致ELBO达到最大时的参数不容易被“发现”; 反之, 当Fisher信息矩阵较大时, 其目标函数曲率较大, 即具有相对窄而高的特点, 最优参数容易被“发现”.

2)在引理1的推导过程中, 在无限小的邻域内(即$ \Delta{\boldsymbol \psi}_k \to 0 $)给出置信下界的Fisher信息矩阵与KL散度之间的近似关系, 并且$ {{{\boldsymbol{J}}}}_{{{\boldsymbol{\psi}}}_k^i} $ 和$ \nabla_{{\boldsymbol \psi}_k}{{\rm{E}}_{q}[\log p({{{\boldsymbol{z}}}}_{k}|{{{\boldsymbol{x}}}}_{k|k}^{i})]} $的计算依赖于参数化的变分分布$ q({{{\boldsymbol{x}}}}_k|{\boldsymbol \psi}_k) $, 因此, 随着变分迭代的进行, Fisher信息矩阵$ {\boldsymbol{J}}_{{\boldsymbol{\psi}}_k} $不断变化.

3)被估计系统模型满足线性高斯的假设条件下, 变分分布和后验分布之间的KL散度具有解析解, 结合引理1能够解析计算变分贝叶斯迭代过程中的Fisher信息矩阵$ {\boldsymbol{J}}_{{\boldsymbol{\psi}}_k} $.

定理 1. 在满足引理1的前提下, 对数似然期望$ {\rm{E}}_{q}\left[\log p\left({\boldsymbol{z}}_{k}|{\boldsymbol{x}}_{k}\right)\right] $在$ {\boldsymbol \psi}_{k} = {\boldsymbol \psi}_{k}^i $的邻域内一阶可导, 则第$ i $次迭代过程中式(28)变分$ {{{\cal{L}}}}\left({\boldsymbol{\psi}}_{k}\right) $关于$ {\boldsymbol \psi}_{k} $的自然梯度为

其中, $ {{{\boldsymbol{J}}}}_{{\boldsymbol \psi}_k^i} $表示在ELBO最大化条件下第$ i $次迭代后$ {{{\cal{L}}}}\left({\boldsymbol{\psi}}_{k}\right) $关于参数$ {\boldsymbol \psi}_k $的Fisher信息矩阵, $\nabla_{{\boldsymbol \psi}_{k}} {\rm{E}}_{q} \left[\log p\left({{{\boldsymbol{z}}}}_{k}|{{{\boldsymbol{x}}}}_{k}\right)\right]_{{\boldsymbol \psi}_{k} = {\boldsymbol \psi}_{k}^{i}}$表示在第$ i $次迭代过程中对数似然期望$ {\rm{E}}_{q}\left[\log p\left({{{\boldsymbol{z}}}}_{k}|{{{\boldsymbol{x}}}}_{k}\right)\right] $在$ {\boldsymbol \psi}_{k}^{i} $处的梯度.

证明. 根据变分参数$ {\boldsymbol \psi}_{k} $的最优估计表达式(28), 通过对其求关于$ {\boldsymbol \psi}_{k} $的偏导数可获得ELBO最大化条件下的极值点. 令$\Delta {\boldsymbol \psi}_{k} := {\boldsymbol \psi}_{k}- {\boldsymbol \psi}_{k}^{i}\to 0$, 根据泰勒展开式(28)中第一项线性化后可表示为

其中, ${{{\rm{O}}}}({\Delta{{\boldsymbol \psi}_{k}}}^{2})$表示 $ {\Delta{{\boldsymbol \psi}_{k}}}^{2} $的高阶无穷小. 结合引理1, 式(28)进一步转化为

计算上式关于$ {\boldsymbol \psi}_{k} $的偏导, 并令其为0, 可获得在第$ i $次变分迭代中ELBO关于$ {\boldsymbol \psi}_{k} $的自然梯度$ \widetilde{\nabla}_{{\boldsymbol \psi}_{k}}^{i} $

□

注2. ELBO的最大化即是在统计流形空间沿基于Fisher信息矩阵的自然梯度向ELBO增大的方向移动, 从而获得状态后验分布的最优近似. 基于ELBO最大化的自然梯度具有以下性质:

1)自然梯度的方向是目标函数在统计流形空间中的最速下降/上升方向, 在本文自然梯度的方向具体表示为使ELBO上升的方向.

2)当$ \Delta{\boldsymbol \psi}_k \to 0 $时, 式(28)中的第二项$ \rm KL $散度近似表示为式(42)中的第二项, 因此, 从局部意义上讲, 变分迭代过程中变分分布与后验分布的近似误差是正定二次型, 即是$ {\Delta{{\boldsymbol \psi}_{k}}}^{2} $的高阶无穷小.

3)在自然梯度的推导过程中, 在无限小的邻域内(即$ \Delta{\boldsymbol \psi}_k \to 0 $)给出Fisher信息矩阵与KL散度之间的近似关系, 并且$ {{{\boldsymbol{F}}}}_{{{\boldsymbol{\psi}}}_k^i} $和$ \nabla_{{\boldsymbol \psi}_k}{{\rm{E}}_{q}[\log p({{{\boldsymbol{z}}}}_{k}|{{{\boldsymbol{x}}}}_{k|k}^{i})]} $的计算依赖于参数化的变分分布$ q({{{\boldsymbol{x}}}}_k|{\boldsymbol \psi}_k) $, 因此, 随着变分迭代的进行, Fisher信息不断变化.

结合式(43)中的自然梯度表达式, 变分分布参数$ {\boldsymbol \psi}_{k} $的更新表示为

由于自然梯度是Fisher有效的^[37-38], 在ELBO最大化条件下基于自然梯度的迭代估计是Fisher有效的估计器. 特别地, 式(44)中变分参数$ {\boldsymbol{\psi}}_k $的迭代估计同时是无偏估计器.

2.2   置信下界自然梯度的信息几何意义

在统计流形空间$ {\cal{M}} $上, 假设概率分布的集合表示为$ {\cal{S}} = \left\{p_{{\boldsymbol \psi}}=p({\boldsymbol{z}}_k;{\boldsymbol{\psi}}_k)\right\} $, 状态向量$ {\boldsymbol{\psi}}_k \in {\bf{R}}^m $. 定义平滑映射${\cal{F}}: {\cal{S}}\longmapsto {\bf{R}}$, 满足参数空间中的$ {\boldsymbol{\psi}}_k $与统计流形空间的$ p_{{\boldsymbol{\boldsymbol \psi}}} $一一对应, 即$ {\boldsymbol{\psi}}_k\equiv {\cal{F}}(p_{{\boldsymbol{\psi}}}) $. 设统计流形空间$ {\cal{M}} $中一点O, 其单位切向量为 $ {{{\boldsymbol{v}}}}_{O} $, 切空间表示为$ {\cal{T}}_O $, 且 $ {{{\boldsymbol{v}}}}_{O}\in {\cal{T}}_O $表示为

其中, $ \{\frac{\partial {{\cal{F}}}({\boldsymbol \psi}_k)}{\partial {{\boldsymbol \psi}_{k,i}}}\}_{i=1,2,\cdots, m} $ 表示统计流形空间 $ {\cal{M}} $中$ m $维坐标系统的一组向量基, $\nabla{{\cal{F}}}({\boldsymbol \psi}_k^O) = $ $ \left(\frac{\partial {{\cal{F}}}({\boldsymbol \psi}_k)}{\partial {{\boldsymbol \psi}_{k,1}}},\; \cdots,\; \frac{\partial {{\cal{F}}}({\boldsymbol \psi}_k)}{\partial {{\boldsymbol \psi}_{k,m}}}\right)^{\rm T}\bigg|_{{\boldsymbol \psi}_k = {{\boldsymbol \psi}_k^O}} $ 表示统计流形空间的梯度, 向量$ {{{\boldsymbol{a}}}} = \{a_1,\; a_2,\; \cdots,\; a_m\}^{\rm T} $, 且满足式(46). 当平滑映射$ {\cal{F}} $表示为对数似然函数期望$ {\rm{E}}[\log p({{{\boldsymbol{z}}}}_k|{\boldsymbol \psi}_k)] $时, 式(46)中相应的期望则表示为式(47)中关于$ {\boldsymbol \psi}_{k} $的Fisher信息矩阵.

统计流形空间中O点处指向其邻域一点P的切向量为$ {{{{\boldsymbol{v}}}}}_{OP} \in {\cal{T}}_O $, 采用单位切向量将其表示为

其中, $ \varepsilon $表示一个无限小量, 统计流形空间中点O沿切向量方向向其邻域中点P移动的示意图如图6所示. 统计流形空间点O和点P之间的相邻函数关系通过一阶泰勒展开表示为

图 6  点O沿切向量方向向点P处移动的示意图

Fig. 6  Point O moves in the tangential direction towards P

下载: 全尺寸图片 幻灯片

同时, 其对应的参数空间变量$ {\boldsymbol \psi}_k^O $ 和 $ {\boldsymbol \psi}_k^P $之间满足

结合式(45) ~ (50), 参数空间中$ {\boldsymbol \psi}_k $变化量的表达式为

文献[22]在式(46)的约束下结合拉格朗日方法给出向量$ {{{\boldsymbol{a}}}} $的表达式,

因此, 参数空间中状态向量从点$ {\boldsymbol \psi}_k^O $到$ {\boldsymbol \psi}_k^P $的移动表示为

其中, 当平滑映射为$ {\rm{E}}[\log p({{{\boldsymbol{z}}}}_k|{\boldsymbol \psi}_k)] $时, $ \nabla{{\cal{F}}}({\boldsymbol \psi}_k^O) = \nabla{\rm{E}}[\log p({{{\boldsymbol{z}}}}_k|{\boldsymbol \psi}_k)] $, 令$\varepsilon = \sqrt{(\nabla {{\cal{F}}}({\boldsymbol \psi}_k))^{\rm T}{{{\boldsymbol{J}}}}_{{\boldsymbol \psi}_k}^{-1}\nabla {{\cal{F}}}({\boldsymbol \psi}_k)}$, 式(53)可转化为

在平滑映射$ {\cal F} $的条件下, 由于统计流形空间中的点(概率分布) $ {{\cal{F}}}({\boldsymbol \psi}_k) $与参数空间中的变量$ {\boldsymbol \psi}_k $一一对应^[22]. 在变分ELBO最大化条件下, $ {{\cal{F}}}({\boldsymbol \psi}_k) $通过沿自然梯度方向移动获得$ {\boldsymbol \psi}_k $的后验分布的近似分布, 从而实现$ {\boldsymbol \psi}_k $的有效估计. 统计流形空间中点$ {{\cal{F}}}({\boldsymbol \psi}_k) $和参数空间中点$ {\boldsymbol \psi}_k $的移动示意图如图7所示.

图 7  变分迭代过程中置信下界自然梯度的示意图

Fig. 7  The natural gradient of ELBO in variational Bayesian iteration

下载: 全尺寸图片 幻灯片

2.3   基于自然梯度的迭代变分贝叶斯滤波实现

针对预备知识部分阐述的非线性状态估计的优化问题, 结合上述ELBO最大化条件下自然梯度的推导, 提出一种非线性状态迭代估计方法, 以寻找一个最优变分分布实现对后验概率密度函数的“紧密”近似. 在高斯假设条件下由于后验分布$p({\boldsymbol{x}}_{k}|{\boldsymbol{z}}_{k})= {\rm{N}}({\boldsymbol{x}}_{k}|{\boldsymbol{x}}_{k|k},{\boldsymbol{P}}_{k|k})$由状态估计值 $ {\boldsymbol{x}}_{k|k} $及其估计误差协方差$ {\boldsymbol{P}}_{k|k} $决定, 因此将变分参数定义为 $ {\boldsymbol \psi}_{k}^{i}:= ({\boldsymbol{x}}_{k|k}^{i},{\boldsymbol{P}}_{k|k}^{i}) .$ 在引理1和定理1的前提下, 结合式(44)中变分分布参数的更新表达式, 第$ i $次变分迭代估计的非线性系统状态估计值$ {\boldsymbol{x}}_{k|k}^{i} $及其误差协方差$ {\boldsymbol{P}}_{k|k}^{i} $可表示为

其中, $ {\boldsymbol{J}}_{{\boldsymbol{x}}_{k|k}^{i}} $和$ {\boldsymbol{J}}_{{\boldsymbol{P}}_{k|k}^{i}} $分别表示在第 $ i $次迭代更新过程中式(44)中$ {{{\boldsymbol{J}}}}_{{\boldsymbol \psi}_k} $关于变分超参数$ {\boldsymbol{x}}_{k|k} $和$ {\boldsymbol{P}}_{k|k} $的Fisher信息矩阵, $ \nabla_{{\boldsymbol{x}}_{k|k}} {\rm{E}}_{q}[\log p({\boldsymbol{z}}_{k}|{\boldsymbol{x}}_{k|k}^i)] $和$\nabla_{{\boldsymbol{P}}_{k|k}} {\rm{E}}_{q} [\log p({\boldsymbol{z}}_k|{\boldsymbol{x}}_{k|k}^i)] $分别表示在第$ i $次迭代更新过程中式(44)中对数似然期望$ {\rm{E}}_q[\log p({\boldsymbol{z}}_k|{\boldsymbol{x}}_{k|k}^i)] $关于$ {\boldsymbol{x}}_{k|k} $和$ {\boldsymbol{P}}_{k|k} $的一阶偏导. 为实现变分迭代过程中状态估计及其误差协方差的更新, 下面给出$ {\boldsymbol{J}}_{{{\boldsymbol{\psi}}}_{k}}^{-1} $和$ \nabla_{{\boldsymbol{\psi}}_k}{{\rm{E}}_{q}\left[\log p\left({{{\boldsymbol{z}}}}_{k}|{\boldsymbol{x}}_k\right)\right]} $关于$ {\boldsymbol{x}}_{k|k} $和$ {\boldsymbol{P}}_{k|k} $的相应具体表达式以及证明过程.

1) Fisher信息矩阵$ {\boldsymbol{J}}_{{\boldsymbol{x}}_{k|k}^{i}} $和$ {\boldsymbol{J}}_{{\boldsymbol{P}}_{k|k}^{i}} $的理论推导

若高斯分布 $q_{1}\sim{\rm{N}}\left({\boldsymbol{\xi}}_{1}; \,{\boldsymbol{\mu}}_{1},{\boldsymbol{C}}_{1}\right)$ 和 $q_{2}\sim {\rm{N}}({\boldsymbol{\xi}}_{2}; \,{\boldsymbol{\mu}}_{2},{\boldsymbol{C}}_{2})$具有相同维数$ d $, 则两分布之间的KL散度可表示为${\cal{D}}_{\text{KL}}[q_{1} \| q_{2}]= \frac{1}{2}\{\ln (|{\boldsymbol{C}}_{2}||{\boldsymbol{C}}_{1}|^{-1}) +\operatorname{tr}({\boldsymbol{C}}_{2}^{-1} \cdot$ ${\boldsymbol{C}}_{1})+({\boldsymbol{\mu}}_{2}-{\boldsymbol{\mu}}_{1})^{\text{T}} {\boldsymbol{C}}_{2}^{-1}({\boldsymbol{\mu}}_{2}-{\boldsymbol{\mu}}_{1}) -d\}$, 其中, 数学符号$ {\rm{tr(\cdot)}} $和 $ \left\vert\cdot\right\vert $分别表示矩阵的迹和矩阵的行列式. 结合引理1, $ \rm ELBO $关于$ {\boldsymbol{x}}_{k|k} $的Fisher信息矩阵可表示为KL散度关于$ {\boldsymbol{x}}_{k|k} $的二阶偏导

相应地, 关于$ {\boldsymbol{P}}_{k|k} $的Fisher信息矩阵表示为KL散度关于$ {\boldsymbol{P}}_{k|k} $的二阶偏导

其中, 哈达玛积$ \circ $表示取矩阵的对角线元素构成单位阵. 由于在估计过程中估计误差协方差的非对角线元素可以忽略不计, 因此, 上式可以近似简化为

根据矩阵的偏导计算知识$ {\boldsymbol{X}} $的逆函数 $ {\boldsymbol{X}}\to {\boldsymbol{X}}^{-1} $的导数表示为

其中, $ m $和$ p $分别表示矩阵$ {\boldsymbol{X}} $的不同行, $ n $和 $ l $分别表示其不同列, 则式(59)中的偏导可表示为

其中, 数学符号$ \otimes $表示克罗内克积, 并且$ ({{{\boldsymbol{B}}}}^{\rm T}\otimes {{{\boldsymbol{A}}}}){{{\boldsymbol{X}}}} = {{{\boldsymbol{A}}}}{{{\boldsymbol{X}}}}{{{\boldsymbol{B}}}} $. 因此, 在迭代优化过程中式(59)关于估计误差协方差的Fisher信息矩阵转化为

□

2)式(55)中偏导$ \nabla_{{\boldsymbol{x}}_{k|k}} {\rm{E}}_{q}[\log p({\boldsymbol{z}}_{k}|{\boldsymbol{x}}_{k|k}^i)] $和式(56)中$ \nabla_{{\boldsymbol{P}}_{k|k}} {\rm{E}}_{q}[\log p({\boldsymbol{z}}_k|{\boldsymbol{x}}_{k|k}^i)] $ 的理论推导

证明. 根据博纳定理^[39]和高斯假设条件下的量测似然表达式, 对数似然期望$ {\rm{E}}_{q}\left[\log p\left({{{\boldsymbol{z}}}}_{k}|{{{\boldsymbol{x}}}}_{k}\right)\right] $关于$ {{{\boldsymbol{x}}}}_{k|k} $的一阶偏导表示为

结合高斯假设条件下似然函数表达式和局部线性化方法, 式(63)可转化为

其中, $ {{{\boldsymbol{H}}}}_k^{i} $表示在第$ i $次迭代过程中非线性量测矩阵的雅克比矩阵, 并且$ {{{\boldsymbol{H}}}}_k^{i} = \frac{\partial {{{\boldsymbol{h}}}}({{{\boldsymbol{x}}}}_{k})}{\partial {{{\boldsymbol{x}}}}_{k}}\big|_{{{{\boldsymbol{x}}}}_{k} = {{{\boldsymbol{x}}}}_{k|k}^{i}} $.

结合普莱斯定理^[39], 一阶偏导 $\nabla_{{{{\boldsymbol{P}}}}_{k|k}}{\rm{E}}_{q} [\log p\cdot \left({{{\boldsymbol{z}}}}_{k}|{{{\boldsymbol{x}}}}_{k}\right)]$表示为

在高斯假设条件下$ \nabla_{{{{\boldsymbol{x}}}}_{k},{{{\boldsymbol{x}}}}_{k}}^{2}\log p\left({{{\boldsymbol{z}}}}_{k}|{{{\boldsymbol{x}}}}_{k}\right) $表示为

由于 $ {\rm{E}}_{q}\left[{{{\boldsymbol{z}}}}_{k}-{{{\boldsymbol{h}}}}_{k}\left({{{\boldsymbol{x}}}}_{k}\right)\right]=0 $, 局部线性化后, 式(65)可表示为

综合上述过程, 在第$ i $次迭代中变分ELBO关于$ {{{\boldsymbol{x}}}}_{k|k} $和$ {{{\boldsymbol{P}}}}_{k|k} $的自然梯度分别表示为

□

因此, 结合变分迭代更新式(55)和(56)关于$ {{{\boldsymbol{x}}}}_{k|k} $、$ {{{\boldsymbol{P}}}}_{k|k} $的具体Fisher信息矩阵式(57)和(61)以及具体一阶偏导式(64)和(67), 变分优化过程中第$ i+1 $次迭代的状态估计$ {{{\boldsymbol{x}}}}_{k|k}^{i+1} $及其估计误差协方差$ {{{\boldsymbol{P}}}}_{k|k}^{i+1} $表示为

综合式(70)和(71)的构建过程, 提出一种基于自然梯度的变分贝叶斯滤波器VBKF-NG, 并以伪代码形式在算法1中给出具体实现步骤. 为了实现算法迭代过程的自适应, 采用以下相对估计误差${{e}}_r$作为判断迭代是否终止的条件.

其中, $ \epsilon $表示一个较小的正实数, 其数值的大小可以根据实际应用场景的精度需求合理设定. 由式(70)和(71)中的状态估计以及估计误差协方差的更新过程可知, 所提算法具有以下特点和优势:

1)状态估计$ {{{\boldsymbol{x}}}}_{k|k}^{i} $的更新以$ {{{\boldsymbol{P}}}}_{k|k}^{i} $为先决条件, 使得$ {{{\boldsymbol{x}}}}_{k|k}^{i} $沿自然梯度方向自适应向真实状态移动. 在滤波初始阶段, $ {{{\boldsymbol{P}}}}_{k|k}^{i} $较大, 增大移动步长以提高新息在迭代更新中的比重, $ {{{\boldsymbol{x}}}}_{k|k}^{i} $向真实位置移动的速度较快, 从而能够较快地到达真实值附近, 即估计器能够较快地收敛; 反之, 当滤波迭代达到收敛时, $ {{{\boldsymbol{P}}}}_{k|k}^{i} $较小, 减小移动步长使 $ {{{\boldsymbol{x}}}}_{k|k}^{i} $的移动速度减缓, 以更精细地寻找后验分布在统计流形中的位置, 从而使估计器能够更好地接近真实值.

2)根据卡尔曼滤波原理, 式(71)中的$ \left({{{\boldsymbol{H}}}}_k^{i}\right)^{\rm{T}} {{{\boldsymbol{R}}}}_k^{-1}{\boldsymbol{H}}_k^i $和$ {{{\boldsymbol{P}}}}_{k|k}^{i} $均为正定矩阵对称矩阵. 因此, 可以判断$ {{{\boldsymbol{P}}}}_{k|k}^{i+1} \leq {{{\boldsymbol{P}}}}_{k|k}^i $, 即在每次迭代过程中估计误差协方差呈减小趋势, 并且由于$ ({{{\boldsymbol{H}}}}_{k}^i)^{\rm{T}} {{{\boldsymbol{R}}}}_k^{-1}{{{\boldsymbol{H}}}}_{k}^i $是海塞阵的期望形式, 表明估计误差协方差$ {{{\boldsymbol{P}}}}_{k|k}^i $具有二阶收敛性.

算法 1. 基于自然梯度的非线性变分贝叶斯滤波器(VBKF-NG)

1) 初始化状态估计$ {\boldsymbol{x}}_{1|1} $、估计误差协方差$ {\boldsymbol{P}}_{1|1} $、最大迭代次数${I}_{\rm ter}$和初始航向角$ \theta_1 $;

2) 计算状态预测和相应的预测误差协方差

其中, $ {\boldsymbol{F}}_{k} = \frac{\partial {\boldsymbol{f}}_{k}\left({\boldsymbol{x}}\right)}{\partial {\boldsymbol{x}}_k}\big|_{{\boldsymbol{x}}={\boldsymbol{x}}_{k-1|k-1}} $;

3) 初始化迭代次数${I}_{\rm ter}$, 状态估计$ {{{\boldsymbol{x}}}}_{k|k}^{1}= {{{\boldsymbol{x}}}}_{k|k-1} $及其估计误差协方差$ {{{\boldsymbol{P}}}}_{k|k}^{1}= {{{\boldsymbol{P}}}}_{k|k-1} $;

4) FOR $i = 1:{I_{\rm{ter}}}$且$ e_r\leq\epsilon $

5) 　根据式(57)和(62)分别计算关于$ {{{\boldsymbol{x}}}}_{k|k}^{i} $和$ {{{\boldsymbol{P}}}}_{k|k}^{i} $的Fisher信息矩阵;

6) 　根据式(63)和(67)分别计算对数似然期望关于$ {{{\boldsymbol{x}}}}_{k|k}^{i+1} $和$ {{{\boldsymbol{P}}}}_{k|k}^{i+1} $的一阶偏导;

7) 　根据式(70)和(71)分别迭代更新状态估计$ {{{\boldsymbol{x}}}}_{k|k}^{i+1} $及其误差协方差$ {{{\boldsymbol{P}}}}_{k|k}^{i+1} $;

8) 　根据式(72)计算相对误差$ e_r $;

9) ENDFOR

10) 输出状态估计值$ {{{\boldsymbol{x}}}}_{k|k}= {{{\boldsymbol{x}}}}_{k|k}^{i+1} $及其估计误差协方差矩阵$ {{{\boldsymbol{P}}}}_{k|k}={{{\boldsymbol{P}}}}_{k|k}^{i+1} $.

3.   仿真算例

3.1   天基光学传感器量测条件下近地轨道卫星定轨跟踪

为验证所提算法的性能, 采用天基光学传感器量测条件下近地轨道卫星定轨跟踪仿真场景对算法进行性能分析. 在J$ 2000.0 $地心惯性坐标系下, 采用三个天基光学传感器对目标进行被动测角, 进而采用算法实现定轨与跟踪.

3.1.1   目标运动模型

近地轨道(Low earth orbit, LEO)卫星轨道高度范围一般在$ 500\sim 2\;000 $ km之间, 不考虑变轨时推力的作用, 其在轨运动主要受地球非球摄动$ {{{\boldsymbol{a}}}}^G(t) $和大气阻力摄动 $ {{{\boldsymbol{a}}}}^D(t) $ 两方面因素影响, 具有典型的强非线性特点. 在$ t $时刻LEO卫星的连续时间动态模型^[40]表示为

其中, $ {{{\boldsymbol{x}}}}(t) = [x(t)\;\; y(t)\;\; z(t)\;\; \dot{x}(t)\;\; \dot{y}(t)\;\; \dot{z}(t)]^{\rm T} $表示在$ t $时刻轨道目标在坐标系中的状态向量, ${{{\boldsymbol{x}}}}^p(t) = [x(t)\;\; y(t)\;\; z(t)]^{\rm T}$ 和 $ \dot {{{\boldsymbol{x}}}}^v(t) = [\dot{x}(t)\;\; \dot{y}(t)\;\; \dot{z}(t)]^{\rm T} $分别表示$ t $时刻轨道目标状态的位置分量和速度分量, 状态转移矩阵${{{\boldsymbol{A}}}} = \left[ \begin{aligned} {{{\boldsymbol{0}}}}_{3\times3} \;\; {{{\boldsymbol{I}}}}_{3\times3} \\ {{{\boldsymbol{0}}}}_{3\times3} \;\; {{{\boldsymbol{0}}}}_{3\times3} \end{aligned}\right]$, 噪声驱动矩阵${{{\boldsymbol{B}}}} = \left[ \begin{aligned} {{{\boldsymbol{0}}}}_{3\times3} \\ {{{\boldsymbol{I}}}}_{3\times3} \end{aligned}\right]$, 其中, $ {{{\boldsymbol{I}}}}_{3\times3} $表示$ 3\times3 $的单位阵, $ {{{\boldsymbol{0}}}}_{3\times3} $表示$ 3\times3 $的零矩阵, $ {\boldsymbol\omega}(t) $表示服从零均值高斯分布系统噪声, 地球非球摄动$ {{{\boldsymbol{a}}}}^G(t) $表示为

其中, $r(t) =\sqrt {x^2(t) + y^2(t) + z^2(t)}$表示轨道目标到地心的距离, $\mu = 3.986\;005\times 10^{14}$ km³/s²表示地球引力系数, $ {C_e} = (3/2){J_2}R_e^2 $, $ {R_e} = 6\;378.137 $ km表示地球赤道半径, $ {J_2} = 1.082 \;63\times10^{-1} $表示修正模型的地球二阶扁率系数. 大气阻力摄动$ {{{\boldsymbol{a}}}}^D(t) $的表达式为

其中, $ ||\dot {{{\boldsymbol{x}}}}^v(t)|| = \sqrt{(\dot{x}(t))^2 + (\dot{y}(t))^2 + (\dot{z}(t))^2} $表示目标速度的大小, $ \rho_0 = 1.29\times 10^{-10} $km/m³表示轨道目标所在轨道高度的大气密度, $ {S}/{m} = 0.02 $ m²/kg表示有效面质比, $ C_D = 2.2 $表示大气阻力系数. 由式(73)、(75)和(76)可以看出轨道目标的连续动态模型与时间相关并具有强非线性特点, 设采样周期为$ {T} $, 相应$ k $时刻离散化的动态模型为

其中, $ {{{\boldsymbol{f}}}}_k({{{\boldsymbol{x}}}}_k) = {{{\boldsymbol{x}}}}_k + {{{\boldsymbol{g}}}}({{{\boldsymbol{x}}}}_k){T} + {{{\boldsymbol{l}}}}\left( {{{\boldsymbol{x}}}}_k\right) {{{\boldsymbol{g}}}}({{{\boldsymbol{x}}}}_k){{T}^2}/{2} $表示目标离散状态转移阵, 并且

其中, 矩阵$ {\boldsymbol{C}}_{3\times3}({{{\boldsymbol{x}}}}_k) $中的元素分别在式(80)中给出, 其中, $ r_k = \sqrt{x_k^2 + y_k^2 + z_k^2} $表示在离散时刻$ k $轨道目标距地心的距离. 矩阵$ {\boldsymbol{D}}_{3\times3}({{{\boldsymbol{x}}}}_k) $表示为式(81), 其中, $ ||\dot {{{\boldsymbol{x}}}}_k^v|| = \sqrt{\dot{x}_k^2 + \dot{y}_k^2 + \dot{z}_k^2} $表示在离散时刻$ k $目标速度的大小. 噪声驱动矩阵表示为

目标的初始状态估计及其误差协方差分别设为$ {{{\boldsymbol{x}}}}_{1|1} = [6.744\,8 \;{\text{m}}\;\; 1.047\,3\;{\text{m}}\; -0.636\,4\;{\text{m}}\; -0.001\,4\;{\text{m/s}} $ $0.007\,6\,{\text{m/s}}\,\, 0.002\,0\,{\text{m/s}}]^{\rm T} \times 10^6\;和$ ${{{\boldsymbol{P}}}}_{1|1} = {\rm diag}\{7\,800\,{\text{m}}^2 \;\; 12\,000\;{\text{m}}^2\; 9\,500\;{\text{m}}^2\; 50\;{\text{m}}^2/{\text{s}}^2\; 80\;{\text{m}}^2/{\text{s}}^2$ $ 80\;{\text{m}}^2/{\text{s}}^2\} $, 过程噪声方差 $ {{{\boldsymbol{Q}}}}_k = {\rm diag}\{2\;{\text{m}}^2\;\;2\;{\text{m}}^2\; \;2\;{\text{m}}^2\;\; 0.3\;{\text{m}}^2/{\text{s}}^2\;$ $ 0.3\;{\text{m}}^2/{\text{s}}^2\; 0.3\;{\text{m}}^2/{\text{s}}^2\} $.

3.1.2   量测模型

假设位于900 km高度的$ 3 $个低轨卫星(星载传感器), 在跟踪时间段内其各传感器均能观测到目标. 三组天基光学传感器对目标进行量测获得其相对方位角和俯仰角信息, 在$ k $时刻的量测为${{{\boldsymbol{z}}}}_{k}= [({{{\boldsymbol{z}}}}_k^1)^{\rm T}\;({{{\boldsymbol{z}}}}_k^2)^{\rm T}\; ({{{\boldsymbol{z}}}}_{k}^3)^{\rm T}]^{\rm T}$, 式(2)中的量测矩阵表示为

其中, $ i = 1,\; 2,\; 3 $, $ {\boldsymbol{x}}_k^{opi} = [x_k^{oi}\;y_k^{oi}\;z_k^{oi}]^{\rm T} $表示第 $ i $个天基光学传感器在坐标系中的位置分量. 三个传感器的初始位置分别为$S^{1}=(6.929 \,3\times10^6 ,\, -2.228\;0\,\times 10^6 , 0)\; \text{m},$ $S^{2}=(6.869\;3\times10^6 ,\; 2.406\;5\times10^6 ,\; 0) \; \text{m}$和$S^{3}=(4.224\;2\times10^6 ,\; 5.927\;5\times10^6 ,\; 0) \; \text{m}$, 采样周期为$ {T} = $1 s. 三个传感器的量测噪声方差相同为 ${{{\boldsymbol{R}}}}_k^1 $ = $ {{{\boldsymbol{R}}}}_k^2 \,=\, {{{\boldsymbol{R}}}}_k^3\, =\, {\rm diag} \{0.64\times 10^{-8}\;\text{rad}^2/\text{s}^2\; 0.64\times\; 10^{-8}\;\text{rad}^2/\text{s}^2\}$.

3.1.3   仿真结果与分析

图8给出天基量测条件下低轨卫星跟踪定轨的仿真场景, 其中黑色曲线表示轨道目标轨迹, 黑色箭头表示其运动方向, 蓝、绿、红三条曲线分别表示天基传感器的运动轨迹, 其量测起始点分别为蓝、绿、红三个三角符号所在位置. 图中地心位置为$ O(0,\; 0,\; 0) $, $ {\boldsymbol{M}}_1 $和${\boldsymbol{T}} $ 分别表示传感器1和目标的示意位置, 向量$ {\boldsymbol{M}}_1{\boldsymbol{T}} $表示传感器1的量测向量, 假设传感器1自身位置是已知量, 目标在地心惯性坐标系下的位置向量表示为$ {\boldsymbol{OT}}_1={\boldsymbol{OM}}_1+ {\boldsymbol{M}}_1{\boldsymbol{T}} $. 同理, 在传感器2和传感器3量测条件下目标的位置向量分别表示为$ {\boldsymbol{OT}}_2={\boldsymbol{OM}}_2+ {\boldsymbol{M}}_2\boldsymbol{T} $和$ {\boldsymbol{OT}}_3= {\boldsymbol{OM}}_3 \,+ {\boldsymbol{M}}_3{\boldsymbol{T}}. $ 在跟踪时间段内, 各传感器均能观测到目标. 仿真中采用EKF、UKF、IEKF作为对比算法, 蒙特卡罗实验次数为$ 1\;000 $. 表2给出目标的轨道根数.

图 8  天基量测条件下LEO跟踪定轨仿真场景

Fig. 8  Scenario of LEO orbit determination and tracking with space-based measurement

下载: 全尺寸图片 幻灯片

表 2  目标的轨道根数

Table 2  The orbital elements of target

半长轴 (km)	离心率	倾角 (deg)	近地点角 (deg)	升交点赤经 (deg)
7500	0.1	15	30	12

下载: 导出CSV 

| 显示表格

图9 ~ 11分别给出所提VBKF-NG算法与EKF、UKF和IEKF算法在坐标系中在x轴、y轴和z轴位置估计RMSE的对比. 根据动态系统连续模型(式(73))和其离散化模型(式(77)), 轨道目标的状态转移具有强非线性特点, IEKF和VBKF-NG因其迭代功能不但能够获得更高精度的量测函数雅克比矩阵, 而且其转移函数的雅克比矩阵也得到有效提升. 从图中可以看出IEKF和VBKF-NG与EKF和UKF相比, 在x轴、y轴和z轴均具有较小的位置估计RMSE. 同时, 由于利用基于自然梯度优化寻找ELBO最大的优势, VBKF-NG能够进一步提升非线性估计的精度, 从图中也可以验证VBKF-NG的位置估计RMSE曲线最低, 具有明显的精度优势.

图 9  x轴位置估计RMSE的对比

Fig. 9  RMSE of position estimation in x axis

下载: 全尺寸图片 幻灯片

图 10  y轴位置估计RMSE的对比

Fig. 10  RMSE of position estimation in y axis

下载: 全尺寸图片 幻灯片

图 11  z轴位置估计RMSE的对比

Fig. 11  RMSE of position estimation in z axis

下载: 全尺寸图片 幻灯片

图12 ~ 14分别给出所提算法与EKF、UKF和IEKF算法在坐标系中在x轴、y轴和z轴速度估计RMSE的对比. 与位置估计RMSE的结果相似, 即在目标速度估计精度方面, VBKF-NG算法与现有典型算法相比能够有效提升状态估计精度.

图 12  x轴速度估计RMSE的对比

Fig. 12  RMSE of velocity estimation in x axis

下载: 全尺寸图片 幻灯片

图 13  y轴速度估计RMSE的对比

Fig. 13  RMSE of velocity estimation in y axis

下载: 全尺寸图片 幻灯片

图 14  z轴速度估计RMSE的对比

Fig. 14  RMSE of velocity estimation in z axis

下载: 全尺寸图片 幻灯片

图15 ~ 20分别给出所提算法与对比算法在坐标轴上的位置和速度估计误差. 为便于查看, 分别对结果曲线进行局部放大. 从图中可以看出, 所提VBKF-NG算法的估计误差最小. 为进一步定量对比分析估计误差, 表3给出1000次蒙特卡罗次数的估计误差均值. 从表中也可以看出所提VBKF-NG算法的估计误差均值最小.

图 15  x轴位置估计误差

Fig. 15  The position estimation errors in x axis

下载: 全尺寸图片 幻灯片

图 16  x轴速度估计误差

Fig. 16  The velocity estimation errors in x axis

下载: 全尺寸图片 幻灯片

图 17  y轴位置估计误差

Fig. 17  The position estimation errors in y axis

下载: 全尺寸图片 幻灯片

图 18  y轴速度估计误差

Fig. 18  The velocity estimation errors in y axis

下载: 全尺寸图片 幻灯片

图 19  z轴位置估计误差

Fig. 19  The position estimation errors in z axis

下载: 全尺寸图片 幻灯片

图 20  z轴速度估计误差

Fig. 20  The velocity estimation errors in z axis

下载: 全尺寸图片 幻灯片

表 3  算法平均估计误差均值的对比

Table 3  Comparison of the mean estimation errors of the algorithm

算法	EKF	UKF	IEKF	VBKF-NG
x 轴位置	6.8731	6.8493	6.8290	6.6025
x 轴速度	−0.0382	0.0418	−0.0368	−0.0241
y 轴位置	2.8793	2.8770	2.8688	2.7563
y 轴速度	−0.0257	−0.0243	−0.0222	−0.0103
z 轴位置	4.7665	4.6759	4.5546	4.3286
z 轴速度	0.1272	0.1097	0.1076	0.1062

下载: 导出CSV 

| 显示表格

为进一步分析算法性能, 以平均单次运行时间$ T_{\rm mean} $作为指标对比算法实时性, 表4给出所提算法与对比算法的平均运行时间. 由于在每次变分贝叶斯迭代过程中需计算量测矩阵的雅克比矩阵和新息方差的逆矩阵, 导致在提升估计精度的同时增加了计算量, 表4中数据表明所提VBKF-NG算法与IEKF相比计算耗时略大.

表 4  算法平均运行时间$( \times 10^{-4}\;{\rm{s}}) $的对比

Table 4  The comparison of the average run time $( \times 10^{-4}\;{\rm{s}}) $ of algorithms

算法	EKF	UKF	IEKF	VBKF-NG
时间	0.2580	0.9542	1.0989	1.1205

下载: 导出CSV 

| 显示表格

3.2   机载被动传感器纯角度目标跟踪场景

为进一步验证所提算法性能, 采用机载红外被动传感器纯角度跟踪仿真场景对算法进行对比和分析. 同时, 采用经典非线性滤波算法 EKF、UKF 和文献[33]中的 VBKF-COD 算法作为对比算法.

3.2.1   目标运动模型

目标相对于传感器的初始位置为$ 5.1 $ km, 初始航向角为顺时针方向$ 220^{\circ} $ (以$ y $轴正方向为$ 0^{\circ} $), 并且以$ 0.123\;5 $ km/min 的速度做匀速直线运动, 状态转移矩阵为

在距传感器5.1 km附近处通过采样方式获得目标初始估计值, 采样方差为初始目标状态估计误差协方差$ {\boldsymbol{P}}_{0|0} $表示为

其中,

其运动轨迹如图21所示.

图 21  纯角度被动跟踪仿真场景

Fig. 21  Scenario of bearing-only passive tracking

下载: 全尺寸图片 幻灯片

3.2.2   传感器运动模型

传感器平台的初始航向角为顺时针$ 140^{\circ} $, 以速度大小为$ 0.154\;3 $ km/min 的固定速度运动, 运动模型分别为 CV 模型 ($ 0 < k < 15, \; 19 \leq k < 34,\; 38 \leq k\leq 50 $) 和匀速转弯 (Constant turn, CT) 模型 $( 15 \leq k < 19,\; 34 \leq k < 38 $).

平台在15 ~ 18时刻发生机动, 然后在19 ~ 33时刻继续进行匀速直线运动, 随后在34 ~ 37时刻再次发生机动, 继而, 在38 ~ 50时刻进行匀速直线运动. 其中, $ \omega = 30\;(^{\circ}/\text{min}) $表示转弯速率, 时间周期为 $ T =1 $ min. 传感器量测模型表示为

其中, $ x_k^t $和$ y_k^t $分别表示目标在坐标轴的位置, $ x_k^o $和$ y_k^o $分别表示传感器在坐标轴的位置. 量测初始方位角$ \phi_0 = 81^{\circ} $, 量测噪声标准差$ \sigma_{0} = \sqrt{R_{0}} = 1^{\circ} $, 径向距标准差$ \Delta r_0 = 2 $ km, 速度标准差$ \Delta v_0 = 61.7\; {\text{m/min}}. $ 仿真实验的蒙特卡罗次数为100.

3.2.3   仿真结果与分析

图22给出在不同滤波算法下的RMSE对比和 PCRLB 曲线. 从图中可以看出: 在传感器发生机动前 ($ k\in [1, 14] $时刻), 系统由于受不可观性的限制使量测累积量不能为滤波器提供足够的信息, 从而导致滤波器EKF、UKF、VBKF-COD和VBKF-NG均不能得到收敛的估计结果. 当传感器发生机动时($ k\in[15, 19] $), 量测在不同角度上累积, 从而系统可观, 并且使滤波具备收敛的前提条件. 从图22(a)和图22(b)上可以看出EKF、UKF、VBKF-COD和VBKF-NG径向距估计和径向速度估计的RMSE均在$ k\in[15,19] $时刻快速下降. 随着更多测量值的累积滤波器均能通过更新过程提取更多的信息, 从而达到滤波的收敛. 从图22(a)和图22(b)中可以看出算法VBKF-NG的径向距估计和径向速度估计的RMSE曲线均低于EKF、UKF和VBKF-COD算法.

图 22  算法的RMSE对比

Fig. 22  Comparison of RMSE for algorithms

下载: 全尺寸图片 幻灯片

在收敛速度方面, 从图22可以看出VBKF-NG算法的RMSE具有较快的收敛速度. 进一步地, 在滤波达到收敛后, VBKF-NG 具有最低的RMSE曲线, 并且最接近于该仿真场景下非线性估计的PCRLB曲线.

为了定量对比算法的估计精度, 表5给出算法达到收敛后(选取$ k\in [25,50] $时刻) 的径向距估计和径向速度估计RMSE平均值的比较. 从表中可以看出: VBKF-NG算法的平均RMSE值最小.

表 5  算法平均RMSE的对比

Table 5  Comparison of average RMSE between algorithms

算法	径向距(km)	径向速度(km/min)
EKF	1.5052	0.0813
UKF	1.1833	0.0670
VBKF-COD	0.9181	0.0432
VBKF-NG	0.2990	0.0161

下载: 导出CSV 

| 显示表格

图23给出算法 VBKF-NG 的归一化 KL 散度在所有时刻随迭代次数的变化, 图中每条曲线代表一个时刻内变分分布与状态后验分布之间的 KL 散度随迭代次数的变化. 从图中可以看出所有 KL 散度曲线均随着迭代次数的增加呈下降趋势, 即在变分贝叶斯框架下通过变分迭代使变分分布逐渐逼近状态后验分布, 从而在滤波过程中达到提升目标状态估计精度的目的. 此外, 随着迭代次数的增加, KL 散度曲线下降的趋势变缓.

图 23  所有时刻的归一化KL散度的变化

Fig. 23  The change in normalized KL divergence at all times

下载: 全尺寸图片 幻灯片

图24给出算法 VBKF-NG 的归一化 ELBO 散度在所有时刻随迭代次数的变化, 图中每条曲线代表一个时刻内 ELBO 随迭代次数的变化. 从图中可以看出所有 ELBO 曲线均随着迭代次数的增加呈递增趋势, 即在每次迭代过程中 ELBO 不断地被最大化, 并且随着迭代次数的增加 ELBO 增加的趋势变缓.

图 24  所有时刻的归一化ELBO的变化

Fig. 24  The change in normalized ELBO at all times

下载: 全尺寸图片 幻灯片

4.   结束语

针对非线性动态系统状态估计问题, 引入在非线性后验分布近似方面具有天然优势的变分贝叶斯方法, 通过构建参数化变分分布将难以获得解析解的非线性状态后验分布的积分问题转化为 ELBO 最大化的优化问题; 继而, 在 ELBO 最大化条件下推导出其自然梯度的一般形式, 并在流形空间中分析其信息几何意义; 在此基础上, 分别推导出 ELBO 关于状态估计及其误差协方差的自然梯度具体表达式, 提出 VBKF-NG 算法以实现非线性系统状态的高精度估计; 最后, 以天基光学传感器量测近地轨道卫星跟踪定轨场景和红外被动传感器量测条件下运动目标跟踪场景进行仿真实验, 实验结果和分析表明: 所提 VBKF-NG 算法与现有典型非线性状态估计算法相比, 具有更高的状态估计精度.

自然梯度与变分贝叶斯理论结合受到了越来越多的国内外学者关注, 随着计算技术和信息处理技术的飞速发展, 其将在目标跟踪、导航与制导和自然语言处理等具有典型复杂系统特征的领域得到充分发展. 以上系统中一般呈现出非线性、参数时变、多模态、高维数及深耦合等特征, 后续将考虑引入变分贝叶斯方法实现系统参数的自学习和状态的高精度估计.