-
摘要: 本文研究了一类存在量测信息缺失情况的目标跟踪问题,提出了一种高斯渐进框架下的目标跟踪方法以实现移动目标的跟踪.考虑可能存在的传感器故障或失效问题,采用假设检验方式以删选错误的量测信息.针对非线性滤波问题,量测信息的缺失将可能引起线性化误差、数值计算误差的增大,从而破坏目标跟踪估计器的稳定性和收敛性.为此,对渐进无迹卡尔曼滤波(Progressive unscented Kalman filter,PUKF)方法进行改进,使其更好地处理量测信息缺失引起的线性化误差、数值计算误差增大的问题.另外,通过对改进PUKF(Modified PUKF,MPUKF)方法的理论分析,证明其可保证渐进过程中的状态估计误差有界.最后,通过一个目标跟踪仿真实例表明,MPUKF方法比传统的IUKF方法和PUKF方法具有更高的跟踪精度.Abstract: This paper is concerned with target tracking problems in the case of incomplete measurements, and a target tracking method presents in the framework of Gaussian progressive filtering. Error measurements are pruned by hypothesis testing under the consideration of possible sensor faults. In nonlinear fitering, stability and convergence of estimator are not guaranteed due to incomplete measurements that may lead to the increase of linearization and numerical calculation errors. Thus, a modified progressive unscented Kalman filter (MPUKF) is proposed to deal with the problem of increase in linearization and numerical calculation errors. Additionally, by theoretical analysis of the MPUKF, it is proved that the estimation errors remain bounded in the progressive process. Simulation of a target tracking example demonstrates that the MPUKF has higher tracking precision than the standard iterated Kalman filter.
-
目标跟踪在军事国防、环境监测、城市交通、家庭服务等领域发挥着重要的作用.随着微电子技术、通信技术的发展, 无线传感器网络(Wireless sensor networks, WSNs)在移动目标跟踪或定位中的应用得到了学术界和工业界的广泛关注[1-5]. WSNs利用大量分散节点对移动目标进行协同感知, 并提供丰富的环境信息以及准确的定位服务.
在移动目标的跟踪过程中, 通常采用无迹卡尔曼滤波(Unscented Kalman filter, UKF)方法[6-7]来处理系统中存在的非线性滤波问题.然而, 由于UKF方法能拟合的阶数十分有限, 在滤波过程中会引入较大的线性化误差, 从而影响滤波器的性能.文献[8]提出了迭代无迹卡尔曼滤波(Iterated UKF, IUKF)方法, 其通过多次量测迭代可一定程度地减小线性化误差以提高系统的滤波精度.此外, 文献[9]证明了迭代卡尔曼滤波(Iterated Kalman filter, IKF)方法中的量测迭代过程可看作是一种高斯-牛顿迭代过程, 因此只有当初始的状态估计值足够接近真实值时, 才能保证滤波器全局收敛.为减少线性化误差与数值计算误差的影响, 文献[10]提出了一种渐进高斯滤波(Progressive Gaussian filter, PGF)方法, 该方法根据贝叶斯法则构造系统状态的同伦函数, 在量测更新过程中, 渐进地引入当前观测信息, 进而得到系统的后验状态.进一步, 文献[11]将文献[10]提出的方法推广至多维情况, 采用混合Dirac模型对连续的概率密度函数(Probability density function, PDF)进行离散化.该方法无需假设系统状态与量测间的联合概率密度函数服从高斯分布, 可取得比线性回归卡尔曼滤波(Linear regression Kalman filter, LRKF) [12-13]方法精度更高的估计结果.然而, 它的滤波性能取决于粒子数目, 其时间复杂度与空间复杂度都将随粒子数目增加而增加.针对高斯滤波问题, 通常可以采用确定性采样的方式对其进行近似, 从而有效地减少系统的计算复杂度.这类滤波方法称为高斯近似滤波(Gaussian approximate filter, GAF)方法, 典型的有无迹卡尔曼滤波方法、容积卡尔曼滤波[14-15]方法等.在高斯渐进滤波框架下, 文献[16]提出了一种渐进高斯近似滤波算法, 其估计精度高于现有的IKF方法和GAF方法.尽管通过多次"量测迭代"可一定程度地减小线性化误差与数值计算误差的影响, 但现有的方法并没有充分考虑到线性化误差或数值计算误差的补偿问题, 易导致量测迭代的次数过多, 进而产生过自信的估计结果.
另一方面, 传感器的故障、失效等情况将引起量测信息缺失.特别地, WSNs环境下的目标跟踪系统不可避免地存在时延、丢包等通信不确定性问题, 这将加剧量测信息的不确定性.针对非线性的目标跟踪问题, 某一采样时刻的量测信息缺失将使得系统的线性化误差与数值计算误差增大, 导致滤波器性能下降甚至滤波发散.文献[17]对扩展卡尔曼滤波算法的局部收敛性进行分析并给出滤波器渐进收敛的充分条件, 同时指出合理的量测噪声协方差设计可扩大其吸引域进而有助于滤波器的收敛.进一步, 针对一般的离散时间系统, 文献[18]由李雅普诺夫的递减条件导出一个线性矩阵不等式(Linear matrix inequality, LMI), 通过选取满足LMI的过程噪声协方差以确保滤波器的稳定性.对一类仅状态模型非线性的系统, 文献[19]提出了一种改进的UKF方法并证明其状态估计误差在均方意义下有界, 同时指出偏大的过程噪声协方差有利于滤波器的稳定.然而, 在迭代卡尔曼滤波方法中, 其量测迭代次数往往很难控制.特别地, 过多的量测迭代次数反而破坏滤波器的稳定性.因此, 量测迭代过程中滤波器的稳定性或收敛性问题有待进一步地研究.
本文考虑一类WSNs环境中存在量测信息缺失的目标跟踪问题, 提出了一种高斯渐进框架下的目标跟踪方法.本文的主要工作在于: 1)为避免错误的量测信息对系统的不利影响, 通过假设检验[20-21]对量测信息进行有效地筛选; 2)分析了高斯渐进滤波框架中的渐进过程, 采用MPUKF方法以处理由量测信息缺失引起的线性化误差、数值计算误差增大的问题; 3)通过对MPUKF方法的稳定性分析, 证明其渐进过程中的状态估计误差在均方意义下有界.仿真结果表明, 相比于IUKF方法与PUKF方法, MPUKF方法具有更好的跟踪性能.
1. 问题描述
考虑一类WSNs环境下的多传感器协同目标跟踪问题.如图 1所示, 将WSNs应用于目标跟踪系统, 一方面可使得多个传感器协同工作, 弥补了单个传感器感知范围有限的缺点, 另一方面也带来了多个传感器覆盖范围切换的问题.如图 2所示, 在目标跟踪过程中, 可能会存在传感器间的相互干扰以及故障等问题.同时, 由于WSNs的引入还带来了时延、丢包等通信不确定性问题.如图 3所示, 无论是网络还是传感器节点自身所引起的量测不确定性问题, 都可能使系统的状态估计误差增大, 从而加剧线性化误差.特别地, 当线性化误差超过滤波器承受范围时, 将导致滤波器性能下降甚至滤波发散.
假设移动目标的运动学模型可以描述如下:
$ \begin{equation} {\pmb x_k} = g\left( {\pmb x}_{k - 1} \right) + {\pmb w_k} \end{equation} $
(1) 其中, 为移动目标状态. ${x_{p, k}}$和${y_{p, k}}$分别为$k$时刻移动目标在$x$轴和$y$轴上的位置, ${x_{v, k}}$和${y_{v, k}}$为其对应的速度. ${\pmb w_k}$为零均值且协方差为${Q_k}$的高斯噪声.
假设WSNs中有${N}$个无线传感器节点, 其观测模型可描述为
$ \begin{align} %\begin{array}{*{20}{l}} z_k^i =& {\beta ^i}\left\{h\left( { \pmb x_k, {\pmb x^i}} \right) + v_k^i \right\} =\nonumber\\ & {\beta ^i}\Bigg\{ {\sqrt {{{\left( {{x_{p, k}} - x_p^i} \right)}^2} + {{\left( {{y_{p, k}} - y_p^i} \right)}^2}}+ } \nonumber\\ &{ v_k^i + {\alpha ^i}\psi _k^i} \Bigg\}, \;\;\;i = 1, 2, \cdots , N \label {eq:2} %\end{array} \end{align} $
(2) 其中, $z_k^i$表示$k$时刻传感器节点$i$的量测值, 即其与移动目标间的距离值. ${\pmb x^i}$为传感器节点$i$的位置坐标, $x_p^i$, $y_p^i$分别为其在$x$轴与$y$轴上的坐标值.量测噪声的均值为零, 协方差为$R_k^i$, 并且与过程噪声$\pmb w_k$无关. $\psi _k^i$为干扰信号, 用于表示移动目标超出传感器节点感知范围以及传感器故障等情况. $\alpha $表示发生故障的概率, $\beta $表示发生量测信息丢失的概率, 其服从$0-1$伯努利分布, $\Pr \left\{ {\beta = 1} \right\} = \bar \beta $和$\Pr \left\{ {\alpha = 1} \right\} = \bar \alpha $, $0 \le \bar \alpha , \bar \beta \le 1$.
注 1. 在实际应用中, 传感器的感知范围往往有限, 如超声波传感器、激光传感器等.若移动目标超出传感器节点的感知范围, 传感器节点将可能返回无效的量测数据.另一方面, WSNs带来的时延、丢包等通信不确定性问题易导致量测信息的缺失.为此, 在观测模型中引入随机变量$\beta$、$\alpha$来描述可能存在的量测信息错误、失效以及缺失等情况.
2. MPUKF滤波算法
2.1 高斯渐进滤波方法
在目标跟踪过程中, 考虑存在传感器故障或相互干扰等情况, 采用假设检验剔除错误的量测信息以提高滤波器的稳定性.记量测新息序列为
$ \begin{equation} \label {eq:3} \pmb e_k = \pmb z_k - \hat {\pmb {z}}_{k\left| {k - 1} \right.} \end{equation} $
(3) 其相应的协方差为
$ \begin{equation} \label {eq:4} {P_{zz, k\left| k \right.}} = \iint {\pmb e_k\pmb e_k^{\rm T} f\left( {\pmb z_k \left| \pmb x_k \right.} \right)f\left( {\pmb x_k \left| {{Z_{k-1}}} \right.} \right){\rm{d}}\pmb z_k {\rm{d}}\pmb x_k } \end{equation} $
(4) 并定义马氏距离的平方$\gamma _k$为
$ \begin{equation} \label {eq:5} \gamma _k = {( { {{{ \pmb e_k )}^{\rm{T}}}{{( {P_{zz, k\left| k \right.}} )}^{ - 1}}\pmb e_k } } } \end{equation} $
(5) 根据新息序列$\pmb e_k $的高斯性, 马氏距离的平方服从相应维度的卡方分布, 即
$\begin{equation} \label {eq:6} \Pr \left( {\gamma _k < \chi _\alpha ^2} \right) = 1 - \alpha \end{equation} $
(6) 其中, $\Pr \left( \cdot \right)$表示随机事件发生的概率, 为显著性水平, $\chi _\alpha ^2$为$1 - \alpha $的置信界.当零假设被拒绝或新息序列落在$1 - \alpha $置信界外时, 则认为发生传感器故障等情况.
若当前时刻的量测信息缺失或错误信息被拒绝, 滤波器仅能根据先验信息对移动目标的状态进行估计.在下一个估计周期, 二次预测将可能导致系统的状态估计误差增大, 即.进一步, 对量测值进行泰勒展开, 即
$ \begin{align} \label {eq:7} \pmb z_k = &h\left( {{\hat{\pmb x}_{k\left| {k - 1} \right.}}} \right) + \nabla h\left({{\hat{\pmb x}_{k\left| {k - 1} \right.}}} \right){\tilde{\pmb x}_{k\left| {k - 1} \right.}}+\nonumber\\ & \frac{1}{2}{\nabla ^2}h\left({{\hat{\pmb x}_{k\left| {k - 1} \right.}}} \right){\tilde{\pmb x}_{k\left| {k - 1} \right.}}^2 + \cdots + \pmb v_k \end{align} $
(7) 不难发现, 多次迭代预测易引起线性化误差的增大.同样地, 如图 4所示, 这将导致系统的先验概率密度函数更加偏离似然函数, 使其重合部分变小.特别地, 若量测噪声协方差较小, 会进一步地减少其重合部分.这样, 仅有少量的采样粒子对积分过程有作用[10-11], 使系统的数值计算误差增大.为此, 引入高斯渐进滤波结构来有效地增加"重合部分", 从而减小数值积分过程中的计算误差.
根据贝叶斯法则, 系统状态的后验分布可表示为:
$ \begin{align} \label {eq:8} {f^e}\left( {\pmb x_k} \right) = &f\left( {\pmb x_k\left| {{Z_k}} \right.} \right) = \frac{{f\left( {\left. {\pmb x_k, \pmb z_k} \right|{Z_{k - 1}}} \right)}}{{f\left( {\left. {\pmb z_k} \right|{Z_{k - 1}}} \right)}} = \nonumber\\ & \frac{{f\left( {\pmb x_k\left| {{Z_{k - 1}}} \right.} \right)f\left( {\pmb z_k\left| {\pmb x_k} \right.} \right)}}{{\int {f\left( {\pmb x_k\left| {{Z_{k - 1}}} \right.} \right)f\left( {\pmb z_k\left| {\pmb x_k} \right.} \right){\rm{d}}\pmb x_k} }} \end{align} $
(8) 引入渐进参数$\lambda $, 由上式构造如下的同伦函数:
$ \begin{equation} \begin{split} {{f^e}\left( {\pmb x_k, \lambda } \right) = c\left( \lambda \right)f\left( {\pmb x_k\left| {{Z_{k{\rm{ - }}1}}} \right.} \right){f^\lambda }\left( {\pmb z_k\left| {\pmb x_k} \right.} \right)} \end{split} \end{equation} $
(9) 其中, 为渐进过程中的后验概率密度函数, $c\left( \lambda \right)$为归一化常数.这里, 系统渐进过程中的联合概率密度函数可表示为
$ \begin{align} &f\left( {\pmb x_k , \pmb z_k, \lambda + \Delta } \right) = \frac{1}{\sqrt {2\pi {({R_k}/\Delta)}}}\times\nonumber\\ &\qquad \exp \Bigg\{ {\left( \pmb z_k - {h\left( {\pmb x_k } \right)} \right)}^{\rm{T}} {{\left( {\frac{{{R_k}}}{\Delta }} \right)}^{ - 1}}\nonumber\\ &\qquad\left( {\pmb z_k - h\left( {\pmb x_k} \right)} \right) \Bigg\}{f^e}\left( {\pmb x_k, \lambda } \right)=\nonumber\\ &\qquad {f^{ml}}( {\pmb x_k} ){f^e}( {\pmb x_k, \lambda } ) \end{align} $
(10) 当$\lambda $从0到1连续变化时, 同伦函数定义了从先验分布到后验分布变化过程中的概率分布.其渐进过程又可表示成
$ \begin{align} &{f^e}\left( {\pmb x_k , n \times \Delta } \right) \propto f\left( {\pmb x_k\left| {{Z_{k - 1}}} \right.} \right)\prod\limits_{i = 1}^n {{f^\Delta }\left( {\pmb z_k \left| {\pmb x_k} \right.} \right)} = \nonumber\\ &\qquad {f^e}\left( {\pmb x_k, \left( {n - 1} \right) \Delta } \right){f^\Delta }\left( {\pmb z_k\left| {\pmb x_k} \right.} \right) \end{align} $
(11) 其中, $\Delta =1/ N$, 经过$N$次迭代之后, 最终得到系统的后验概率分布.如图 5所示, 系统通过同伦函数递推的方式, 可渐进地引入量测信息对先验概率密度函数进行修正, 有效地利用中间后验概率密度函数, 从而抑制数值计算误差的增大.
2.2 MPUKF方法设计
MPUKF方法主要由假设检验与系统状态迭代更新两个部分组成.在状态迭代更新过程中, 系统渐进地引入量测信息对当前状态进行修正, 即通过多次量测迭代得到对应时刻的后验状态.为减少迭代过程中的计算量, 假设系统的先验概率密度函数与似然函数服从高斯分布, 进而根据贝叶斯法则推导出其后验概率密度函数服从高斯分布.类似地, 在系统量测渐进过程中, 同样假设系统状态与量测间的联合概率密度函数服从高斯分布, 得到$\lambda + \Delta $时刻的联合高斯分布为
$ \begin{equation} f\left( {\pmb x_k , \pmb z_k, \lambda + \Delta } \right) = {\rm N}\left( {\left[ {\begin{array}{*{20}{c}} \pmb x_k \\ \pmb z_k \end{array}} \right];{m_{xz}}, {P_{xz}}} \right) \end{equation} $
(12) 其中, , . ${\hat{\pmb x}_{k\left| k \right.}^\lambda }$与$P_{xx, k\left| {k } \right.}^\lambda$为后验概率密度函数${f^e}\left( {\pmb x_k, \lambda } \right)$的前二阶距, $ {\hat{\pmb z}_{k\left| k \right.}^{\lambda + \Delta }}$与$P_{zz, k\left| k \right.}^{\lambda + \Delta }$为量测预测概率密度函数${f^z}\left( {\pmb z_k, \lambda + \Delta } \right)$的前二阶距, $P_{xz, k\left| k \right.}^{\lambda + \Delta }$为系统状态与量测值间的互协方差.
根据量测预测值与协方差$P_{zz, k\left| k \right.}^{\lambda + \Delta }$, $P_{xz, k\left| k \right.}^{\lambda + \Delta }$的定义并结合UT变换可得
$ \begin{align} \label {eq:12} {\hat{\pmb z}_{k\left| k \right.}^{\lambda + \Delta } } = &\int {\int {\pmb z_k f\left( {\pmb x_k , \pmb z_k , \lambda + \Delta } \right){\rm{d}}\pmb z_k{\rm{d}}\pmb x_k} }=\nonumber\\ &{ \int {\left\{ {\int {\pmb z_k{f^{ml}}\left( {\pmb x_k} \right){\rm{d}}\pmb z_k} } \right\}{f^e}\left( {\pmb x_k, \lambda } \right){\rm{d}}\pmb x_k} }=\nonumber\\ &{ \int {h\left( {\pmb x_k} \right){f^e}\left( {\pmb x_k, \lambda } \right){\rm{d}}\pmb x_k} }\approx { \sum\limits_{j = 0}^{2L} {w_j^{\left( m \right)}h\left( {\chi _{j, k\left| k \right.}^\lambda } \right)} \;} \end{align} $
(13) $ \begin{align} \label {eq:13} P_{zz, k\left| k \right.}^{\lambda + \Delta } = &\int {\int {\left( {\pmb z_k - { \hat{\pmb z}_{k\left| k \right.}^{\lambda + \Delta }}} \right) {{\left( {\pmb z_k -{\hat{\pmb z}_{k\left| k \right.}^{\lambda + \Delta }}} \right)}^{\rm{T}}}} } \times\nonumber\\ & f\left( {\pmb x_k, \pmb z_k, \lambda + \Delta } \right) {\rm{d}}\pmb z_k{\rm{d}}\pmb x_k =\nonumber \\ & \int {\int {\pmb z_k\pmb z_k^{\rm{T}}} } f\left( {\pmb x_k, \pmb z_k, \lambda + \Delta } \right){\rm{d}}\pmb z_k {\rm{d}}\pmb x_k-\nonumber\\ & {\hat{\pmb z}_{k\left| k \right.}^{\lambda + \Delta }} {\left({{\hat{\pmb z}_{k\left| k \right.}^{\lambda + \Delta }}} \right)^{\rm{T}}}=\nonumber\\ & \int {\left[ {h\left( {\pmb x_k} \right){h^{\rm{T}}}\left( {{\pmb x_k}} \right) + \frac{{{R_k}}}{\Delta }} \right]} {f^e}\left( {{\pmb x_k}, \lambda } \right){\rm{d}}{\pmb x_k}-\nonumber\\ & \hat {\pmb z}_{k\left| k \right.}^{\lambda + \Delta }{\left( {\hat {\pmb z}_{k\left| k \right.}^{\lambda + \Delta }} \right)^{\rm{T}}} \approx \sum\limits_{j = 0}^{2L} {w_j^{\left( c \right)}} \left[ {h\left( {\chi _{j, k\left| k \right.}^{\lambda + \Delta }} \right)} \right. -\nonumber\\ &\left. { \hat {\pmb z}_{k\left| k \right.}^{\lambda + \Delta }} \right]{\left[ {h\left( {\chi _{j, k\left| k \right.}^{\lambda + \Delta }} \right) - \hat {\pmb z}_{k\left| k \right.}^{\lambda + \Delta }} \right]^{\rm{T}}} + \frac{{{R_k}}}{\Delta } \end{align} $
(14) $ \begin{align} \label {eq:14} P_{xz, k\left| k \right.}^{\lambda + \Delta } = &\int {\int {\left( {{\pmb x_k} - \hat {\pmb x}_{k\left| k \right.}^\lambda } \right)} } {\left( {{\pmb z_k} - \hat {\pmb z}_{k\left| k \right.}^{\lambda+ \Delta} } \right)^{\rm{T}}} \times\nonumber\\ & f\left( {{\pmb x_k}, } \right.\left. {{\pmb z_k}, \lambda + \Delta } \right){\rm{d}}{\pmb z_k}{\rm{d}}{\pmb x_k} = \nonumber\\ &\int {\int {{\pmb x_k} \pmb z_k^{\rm{T}}} } \left. {f\left( {{\pmb x_k}, } \right.{\pmb z_k}, \lambda + \Delta } \right)\times\nonumber\\ & {\rm{d}}{\pmb z_k}{\rm{d}}{\pmb x_k} -\hat {\pmb x}_{k\left| k \right.}^\lambda {\left( {\hat {\pmb z}_{k\left| k \right.}^{\lambda + \Delta }} \right)^{\rm{T}}}=\nonumber\\ &\int {{\pmb x_k}} {h^{\rm{T}}}\left( {{\pmb x_k}} \right) {f^e}\left( {{\pmb x_k}, \lambda } \right){\rm{d}}{\pmb x_k} - \nonumber\\ & \hat {\pmb x}_{k\left| k \right.}^\lambda {\left( {\hat {\pmb z}_{k\left| k \right.}^{\lambda + \Delta }} \right)^{\rm{T}}} \approx \nonumber\\ & \sum\limits_{j = 0}^{2L} {w_j^{\left( c \right)}} \;\left[ {\chi _{j, k\left| k \right.}^{\lambda + \Delta }} \right. - \left. {\hat {\pmb x}_{k\left| k \right.}^\lambda } \right] \times \nonumber\\ & {\left[ {h\left( {\chi _{j, k\left| k \right.}^{\lambda + \Delta }} \right) - \hat {\pmb z}_{k\left| k \right.}^{\lambda + \Delta }} \right]^{\rm{T}}} \end{align} $
(15) 其中, ${\chi _{j, k\left| {k } \right.}^{\lambda + \Delta }}$与$w_j$分别为对应的sigma点与权值.最后由式$(13)\sim(15)$可得
$ \begin{equation} \label {eq:15} \begin{split} {{f^e}\left( {{\pmb x_k}, \lambda + \Delta} \right)} = {\rm N}\left( {{\pmb x_k};\hat {\pmb x}_{k\left| k \right.}^{\lambda + \Delta }, P_{k\left| k \right.}^{\lambda + \Delta }} \right) \end{split} \end{equation} $
(16) 其中,
$ \begin{equation} \label {eq:16} \hat {\pmb x}_{k\left| k \right.}^{\lambda + \Delta } = \hat {\pmb x}_{k\left| k \right.}^\lambda + K_k^{\lambda + \Delta }\left( {{\pmb z_k} - \hat {\pmb z}_{k\left| {k - 1} \right.}^{\lambda + \Delta }} \right) \end{equation} $
(17) $ \begin{equation} \label {eq:17} P_{k\left| k \right.}^{\lambda + \Delta } = P_{k\left| k \right.}^\lambda - K_k^{\lambda + \Delta }P_{zz, k\left| k \right.}^{\lambda + \Delta }(K_k^{\lambda + \Delta})^{\rm T} \end{equation} $
(18) $ \begin{equation} \label {eq:18} K_k^{\lambda + \Delta } = P_{xz, k\left| k \right.}^{\lambda + \Delta }{\left( {P_{zz, k\left| k \right.}^{\lambda + \Delta }} \right)^{ - 1}} \end{equation} $
(19) 选取迭代步长$\Delta = 1/N$, 用$i+1$与$i$分别表示$\lambda + \Delta $与$\lambda$, 并通过状态迭代更新得到对应时刻的后验状态.
然而, 系统中增大的线性化误差与数值计算误差将引起量测信息的"可信度"降低.特别地, 当量测噪声协方差较小时, 线性化误差与数值计算误差更易破坏系统的稳定性, 从而产生不相容[22-23]的估计结果.因此, 为防止出现过估计的情况, 引入判定条件
$ \begin{equation} \label {eq:51} \begin{split} {\rm E}\left\{ {{{\left( {\pmb e_k^{i + 1}} \right)}^{\rm{T}}}\pmb e_k^{i + 1}} \right\} \le {\rm E}\left\{ {{{\left( {\pmb e_k^i} \right)}^{\rm{T}}}\pmb e_k^i} \right\} \end{split} \end{equation} $
(20) 其中, ${\pmb e_k^i}$表示$k$时刻第$i$次迭代对应的量测新息序列.判定条件的引入可改变渐进过程中的迭代次数, 间接地调整量测噪声协方差, 进而补偿由线性化误差、数值计算误差引起的不确定性.记$k$时刻的状态预测及其协方差为与$P_{k\left| k \right.}^0 $, MPUKF方法的具体算法流程可描述如下:
算法1. MPUKF算法
1) 初始化;
2) while
3) 计算与$P_{k\left| k \right.}^0$;
4) 式$(3)\sim(5)$;
5) if
6) break;
7) end if;
8) for $(i = 1:N)$
9) 式$(13)\sim(15)$;
10) if
11) break;
12) end if;
13) 式$(17)\sim(19)$;
14) end for;
15) $k = k+1$;
16) end while.
注 2. 为了便于非线性卡尔曼滤波器的实现, 通常假设系统的状态预测误差正交于量测噪声$\pmb v_k$, 即.类似地, 在MPUKF方法的渐进过程中, 同样认为系统状态估计误差${\tilde {\pmb x}_{k\left| k \right.}^i}$正交于量测噪声, 即.
注 3. 在非线性滤波过程中, 线性化误差通常是不可避免的.过程噪声大小、采样间隔以及非线性强度都将影响线性化误差的大小.特别地, 当某一时刻的量测信息缺失时(可看作采样间隔增大), 将可能造成线性化误差的增大. MPUKF方法可有效地减小线性化误差, 同时提高系统对不同线性化误差的自适应能力, 进而改善滤波器性能.此外, 该方法同样适用于其他因素引起的线性化误差过大而导致滤波器性能下降的情况.
2.3 系统稳定性分析
本节将通过李雅普诺夫方法对渐进量测更新过程中的稳定性问题进行分析.在这之前,首先介绍一些重要的变量和等式.将进行泰勒展开, 根据UT变换并结合文献[6]中的式$(69)\sim(78)$可得
$ \begin{align} \label {eq:21} {{\hat {\pmb z}}_{k\left| {k - 1} \right.}} = & \frac{1}{{L + \kappa }} \Bigg\{ \kappa h\left( {{{\hat{\pmb x}}_{k\left| {k{\rm{ - }}1} \right.}}} \right) + \nonumber\\ &{ \frac{1}{2}\sum\limits_{i = 1}^L {h\left[ {{{\hat {\pmb x}}_{k\left| {k{\rm{ - }}1} \right.}} + {{\left( {\sqrt {\left( {L + \kappa } \right){P_{xx}}} } \right)}_i}} \right]} } + \nonumber\\ &\frac{1}{2}\sum\limits_{i = L + 1}^{2L} h\Bigg[ {{\hat {\pmb x}}_{k\left| {k{\rm{ - }}1} \right.}} - \nonumber\\ &{{\left( {\sqrt {\left( {L + \kappa } \right){P_{xx}}} } \right)}_i} \Bigg] \Bigg\}=h\left( {{{\hat {\pmb x}}_{k\left| {k{\rm{ - }}1} \right.}}} \right) +\nonumber\\ & \frac{1}{2}{\nabla ^2}h\left( {{{\hat {\pmb x}}_{k\left| {k{\rm{ - }}1} \right.}}} \right){P_{xx}} + \cdots \end{align} $
(21) 其中, 表示对称阵的Cholesky分解, $i$为分解阵的第$i$列; $P_{xx}$为对应的协方差矩阵; $\kappa$描述系统的分布信息, $L$为系统状态的维度.
记第$i$次量测迭代后的状态估计误差为
$ \begin{equation} \label {eq:50} \tilde {\pmb x}_{k\left| k \right.}^i = {\pmb x_k} - \hat {\pmb x}_{k\left| k \right.}^i \end{equation} $
(22) 并由式$(3)$可得
$ \begin{equation} \begin{split} &\pmb e_k^i = {\pmb z_k} - \hat {\pmb z}_{k\left| {k - 1} \right.}^i \approx \nabla h\left( {\hat {\pmb x}_{k\left| k \right.}^{i - 1}} \right)\tilde {\pmb x}_{k\left| k \right.}^{i - 1} + {\pmb v_k} \end{split} \end{equation} $
(23) 令, 为$h\left( \pmb x \right)$在处的一阶偏导.为刻画系统的线性化误差, 引入对角矩阵, 使得
$ \begin{equation} \begin{split} {\pmb e}_k^i = \beta _k^{i-1}H_k^{i-1}\tilde {\pmb x}_{k\left| k \right.}^{i - 1} + {\pmb v_k} \end{split} \end{equation} $
(24) 下面将通过定理1证明量测渐进过程中的估计误差有界.
定理1. 若存在实数, 使未知的对角矩阵${\beta _k^i}$满足
$ \begin{equation} \label {eq:25} \begin{split} {\Xi _{\min }}I \le \left| {{{\left( {\beta _k^iH_k^i} \right)}^{\rm{T}}}\beta _k^iH_k^i} \right| \le {\Xi _{\max }} I \end{split} \end{equation} $
(25) 同时满足条件$(20)$, 则系统的状态估计误差在均方意义下有界.
证明. 选取李雅普诺夫函数为
$ \begin{equation} \label {eq:24} \begin{split} {V}\left( {\tilde {\pmb x}_{k\left| k \right.}^i} \right) = {\left( {\tilde {\pmb x}_{k\left| k \right.}^i} \right)^{\rm{T}}}\Xi _k^i\tilde {\pmb x}_{k\left| k \right.}^i \end{split} \end{equation} $
(26) 其中, , 根据式$(25)$可知, ${\Xi _k}$有界, 即可得
$ \begin{equation} \label {eq:26} \begin{split} {\Xi _{\min }}{\left\| {\tilde {\pmb x}_{k\left| k \right.}^i} \right\|^2} \le \left\| {{V}\left( {\tilde {\pmb x}_{k\left| k \right.}^i} \right)} \right\| \le {\Xi _{\max }}{\left\| {\tilde {\pmb x}_{k\left| k \right.}^i} \right\|^2} \end{split} \end{equation} $
(27) 同时, 由条件$(20)$可有
$ \begin{align} &{\rm E}\left\{ {{{\left( {\beta _k^iH_k^i\tilde {\pmb x}_{k\left| k \right.}^i + {\pmb v_k}} \right)}^{\rm{T}}}\left( {\beta _k^iH_k^i\tilde {\pmb x}_{k\left| k \right.}^i + {\pmb v_k}} \right)} \right\}\le\nonumber\\ &\quad~~~ {\rm E}\Bigg\{ {{{\left( {\beta _k^{i - 1}H_k^{i - 1}\tilde {\pmb x}_{k\left| k \right.}^{i - 1} + {\pmb v_k}} \right)}^{\rm{T}}}} \times\nonumber\\ &\quad~~~{ \left( {\beta _k^{i - 1}H_k^{i - 1}\tilde {\pmb x}_{k\left| k \right.}^{i - 1} + {\pmb v_k}} \right)}\Bigg\} \end{align} $
(28) 整理上式可得
$ \begin{align} &{\rm E}\left\{ {{{\left( {\tilde {\pmb x}_{k\left| k \right.}^i} \right)}^{\rm{T}}}\Xi _k^i\tilde {\pmb x}_{k\left| k \right.}^i} \right\} \le\nonumber \\ &\;\;\;\;\;\;\;\;\;\;\;\; ~~~ {\rm E}\left\{ {{{\left( {\tilde {\pmb x}_{k\left| k \right.}^{i - 1}} \right)}^{\rm{T}}}\Xi _k^{i - 1}\tilde {\pmb x}_{k\left| k \right.}^{i - 1}} \right\} \end{align} $
(29) 进而有
$ \begin{align} \label {eq:29} &{\rm E}\left\{ {{V}\left( {\tilde {\pmb x}_{k\left| k \right.}^i} \right)\left| {\tilde {\pmb x}_{k\left| k \right.}^{i - 1}} \right.} \right\} \le\nonumber\\ &\;\;\;\;\;\;\;\;\;\;\;\;{\rm E}\left\{ {{V}\left( {\tilde {\pmb x}_{k\left| k \right.}^{i - 1}} \right)\left| {\tilde {\pmb x}_{k\left| k \right.}^{i - 1}} \right.} \right\} = {V}\left( {\tilde {\pmb x}_{k\left| k \right.}^{i - 1}} \right) \end{align} $
(30) 由文献[24]中的定理2可知
$ \begin{align} \label {eq:30} {\rm E}\left\{ {{V}\left( {\tilde {\pmb x}_{k\left| k \right.}^i} \right)\left| {\tilde {\pmb x}_{k\left| k \right.}^{i - 2}} \right.} \right\} =\nonumber\\ {\rm E}\left\{ {{\rm E}\left\{ {{V}\left( {\tilde {\pmb x}_{k\left| k \right.}^i} \right)\left| {\tilde {\pmb x}_{k\left| k \right.}^{i - 1}} \right.} \right\}\left| {\tilde {\pmb x}_{k\left| k \right.}^{i - 2}} \right.} \right\} \end{align} $
(31) 并结合式$(30)$有
$ \begin{equation} \begin{split} {{\rm E}\left\{ {\left. {{V}\left( {\tilde {\pmb x}_{k\left| k \right.}^i} \right)} \right|\tilde {\pmb x}_{k\left| k \right.}^{i - 2}} \right\} \le {\rm E}\left\{ {\left. {{V}\left( {\tilde {\pmb x}_{k\left| k \right.}^{i - 1}} \right)} \right|\tilde {\pmb x}_{k\left| k \right.}^{i - 2}} \right\}} \end{split} \end{equation} $
(32) 进一步, 整理上式有
$ \begin{equation} \begin{split} {\rm E}\left\{ {{V}\left( {\tilde {\pmb x}_{k\left| k \right.}^i} \right)\left| {\tilde {\pmb x}_{k\left| k \right.}^{i - 2}} \right.} \right\} \le {V}\left( {\tilde {\pmb x}_{k\left| k \right.}^{i - 2}} \right) \end{split} \end{equation} $
(33) 以此类推, 最终可得
$ \begin{equation} \begin{split} {\rm E}\left\{ {\left. {{V}\left( {\tilde {\pmb x}_{k\left| k \right.}^i} \right)} \right|\tilde {\pmb x}_{k\left| k \right.}^0} \right\} \le {V}\left( {\tilde {\pmb x}_{k\left| k \right.}^0} \right) \end{split} \end{equation} $
(34) 由式$(27)$可知
$ \begin{equation} \begin{split} {\rm E}\left\{ {{{\left\| {\tilde {\pmb x}_{k\left| k \right.}^i} \right\|}^2}} \right\} \le \frac{{{\Xi _{\max }}}}{{{\Xi _{\min }}}}{\left\| {\tilde {\pmb x}_{k\left| k \right.}^0} \right\|^2} \end{split} \end{equation} $
(35) 因此, 在线性化误差有界的情况下, 条件$(20)$可保证李雅普诺夫函数单调递减, 使得状态估计误差在均方意义下有界, 从而证明了渐进量测更新过程的稳定.
3. 仿真结果与分析
为验证MPUKF方法的合理性与有效性, 本文通过一个目标跟踪的仿真实例, 在由3个测距传感器组成的WSNs环境中, 对比IUKF, PUKF以及MPUKF方法的滤波效果.考虑在目标跟踪过程中可能存在传感器故障、失效以及相互干扰等情况, 采用如式$(2)$的观测模型并假设移动目标的运动学模型如下
$ \begin{equation} \begin{split} {\pmb x_k} = \left[ {\begin{array}{*{20}{c}} 1&{\Delta t}&0&0\\ 0&1&0&0\\ 0&0&1&{\Delta t}\\ 0&0&0&1 \end{array}} \right]{\pmb x_{k - 1}} + {\pmb w_k} \end{split} \end{equation} $
(36) 其中, 系统的周期采样时间$\Delta t=1$ s, 过程噪声与量测噪声$\pmb v_k$对应的协方差分别为, ${R_k} = 2 \times {10^{ - 3}}{ I_{3 \times 3}}$.设初始真实状态向量${\pmb x_{0\left| 0 \right.}} =$ 及其初始误差协方差.同时, 状态估计的初始值${\hat {\pmb x}_0}$根据随机选取并取发生故障的概率和发生量测信息丢失的概率.
为便于仿真结果分析与比较, 定义位置与速度的误差指标(Logarithmic mean square erros, LMSEs)为:
$ \begin{align} {\rm LMSE} = {\lg}\left\{ {\frac{1}{M}\sum\limits_{j = 1}^M {\left[ {{{\left( {{x_k} - {{\hat x}_{k\left| k \right.}}} \right)}^2}} \right.} } \right. + \nonumber\\ \left. {\left. {\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; {{\left( {{y_k} - {{\hat y}_{k\left| k \right.}}} \right)}^2}} \right]} \right\} \end{align} $
(37) 其中, $M$为仿真次数, ${x_k}$和${y_k}$分别为移动目标在$x$轴和$y$轴上的真实值, 和${\hat y_{k\left| k \right.}}$为对应的估计值.
图 6、图 7分别给出IUKF, PUKF与MPUKF方法的MC仿真与${\rm LMSE_{vel}}$的误差分布曲线.显而易见, MPUKF方法的跟踪精度高于IUKF, PUKF方法.在目标移动的过程中, 增大的估计误差会引入较大的线性化误差, 从而导致滤波器的性能下降.相比于IUKF方法与PUKF方法, MPUKF方法能够很好地处理线性化误差、数值计算误差增大的问题.因此, 随仿真步数的增加, MPUKF方法的性能明显优于IUKF方法与PUKF方法.
进一步, 表 1给出了各个滤波器对应的LMSEs均值, 其中MPUKF方法的LMSEs均值最小.此外, 将上述算例中的IUKF、PUKF以及MPUKF方法的执行时间进行对比.由表 2可知, 在相同的迭代次数下, MPUKF方法的耗时明显小于IUKF方法与PUKF方法, 具有更高的计算效率.
表 1 各滤波器${\rm LMSE_{pos}}$与的均值Table 1 Mean of LMSE$ _{\rm{pos}}$ and LMSE$_{\rm {vel}}$估计方法 LMSE$_{\rm{pos}}$ LMSE$_{\rm{vel}}$ IUKF $\left( {N = 30, \Delta = 1} \right)$ $-$4.372 $-$7.894 PUKF $\left( {N = 30, \Delta = 1{\rm{/}}30} \right)$ $-$4.981 $-$8.592 MPUKF $\left( {N = 30, \Delta = 1{\rm{/}}30} \right)$ $-$5.191 $-$8.787 表 2 各滤波器平均执行时间Table 2 Average running time估计方法(执行参数) 执行时间/s IUKF $\left( {N = 30, \Delta = 1} \right)$ 0.126 PUKF $\left( {N = 30, \Delta = 1{\rm{/}}30} \right)$ 0.156 MPUKF $\left( {N = 30, \Delta = 1{\rm{/}}30} \right)$ 0.108 4. 结论
本文提出了一种高斯渐进滤波框架下的目标跟踪方法.该方法通过假设检验对量测信息进行筛选, 避免错误的量测信息对系统产生不利影响.采用MPUKF方法可有效地减小线性化误差与数值计算误差, 同时提高系统对不同线性化误差的自适应能力.另外, 通过对滤波器稳定性的分析, 证明了在线性化误差有界的情况下, MPUKF方法可保证渐进过程中的状态估计误差有界.
-
表 1 各滤波器${\rm LMSE_{pos}}$与${\rm LMSE_{vel}}$的均值
Table 1 Mean of LMSE$ _{\rm{pos}}$ and LMSE$_{\rm {vel}}$
估计方法 LMSE$_{\rm{pos}}$ LMSE$_{\rm{vel}}$ IUKF $\left( {N = 30, \Delta = 1} \right)$ $-$4.372 $-$7.894 PUKF $\left( {N = 30, \Delta = 1{\rm{/}}30} \right)$ $-$4.981 $-$8.592 MPUKF $\left( {N = 30, \Delta = 1{\rm{/}}30} \right)$ $-$5.191 $-$8.787 表 2 各滤波器平均执行时间
Table 2 Average running time
估计方法(执行参数) 执行时间/s IUKF $\left( {N = 30, \Delta = 1} \right)$ 0.126 PUKF $\left( {N = 30, \Delta = 1{\rm{/}}30} \right)$ 0.156 MPUKF $\left( {N = 30, \Delta = 1{\rm{/}}30} \right)$ 0.108 -
[1] 杨旭升, 张文安, 俞立.适用于事件触发的分布式随机目标跟踪方法.自动化学报, 2017, 43(8):1393-1401 http://www.aas.net.cn/CN/abstract/abstract19113.shtmlYang Xu-Sheng, Zhang Wen-An, Yu Li. Distributed tracking method for maneuvering targets with event-triggered mechanism. Acta Automatica Sinica, 2017, 43(8):1393-1401 http://www.aas.net.cn/CN/abstract/abstract19113.shtml [2] Yang X S, Zhang W A, Yu L, Xing K X. Multi-Rate distributed fusion estimation for sensor network-based target tracking. IEEE Sensors Journal, 2016, 16(5):1233-1242 doi: 10.1109/JSEN.2015.2497464 [3] 熊志刚, 黄树彩, 赵炜, 苑智玮, 徐晨洋.均方根嵌入式容积粒子PHD多目标跟踪方法.自动化学报, 2017, 43(2):238-247 http://www.aas.net.cn/CN/abstract/abstract19002.shtmlXiong Zhi-Gang, Huang Shu-Cai, Zhao Wei, Yuan Zhi-Wei, Xu Chen-Yang. Square-root imbedded cubature particle PHD multi-target tracking algorithm. Acta Automatica Sinica, 2017, 43(2):238-247 http://www.aas.net.cn/CN/abstract/abstract19002.shtml [4] Zhang W A, Yang X S, Yu L. Sequential fusion estimation for RSS-based mobile robots localization with event-driven WSNs. IEEE Transactions on Industrial Informatics, 2016, 12(4):1519-1528 doi: 10.1109/TII.2016.2585350 [5] 黄玉龙, 张勇刚, 李宁, 赵琳.一种改进的高斯近似滤波方法.自动化学报, 2016, 42(3):385-401 http://www.aas.net.cn/CN/abstract/abstract18828.shtmlHuang Yu-Long, Zhang Yong-Gang, Li Ning, Zhao Lin. An improved Gaussian approximate filtering method. Acta Automatica Sinica, 2016, 42(3):385-401 http://www.aas.net.cn/CN/abstract/abstract18828.shtml [6] Julier S, Uhlmann J K. A general method for approximating nonlinear transformations of probability distributions. Technical report, Robotics Research Group, Department of Engineering Science, University of Oxford, 1996 [7] Julier S J, Uhlman J K. Unscented filtering and nonlinear estimation. Proceedings of the IEEE, 2004, 92(3):401-422 http://d.old.wanfangdata.com.cn/Periodical/xtgcydzjs-e200801002 [8] Zhan R H, Wan J W. Iterated unscented Kalman filter for passive target tracking. IEEE Transactions on Aerospace and Electronic Systems, 2007, 43(3):1155-1163 doi: 10.1109/TAES.2007.4383605 [9] Bell B M, Cathey F W. The iterated Kalman filter update as a Gauss-Newton method. IEEE Transactions on Automatic Control, 1993, 38(2):294-297 doi: 10.1109/9.250476 [10] Hanebeck U D, Briechle K, Rauh A. Progressive Bayes: a new framework for nonlinear state estimation. In: Proceedings Volume 5099, Multisensor, Multisource Information Fusion: Architectures, Algorithms, and Applications 2003. Orlando, Florida USA: SPIE, 2003. 256-267 [11] Hanebeck U D, Steinbring J. Progressive Gaussian filtering based on Dirac Mixture approximations. In: Proceedings of the 15th International Conference on Information Fusion. Singapore: IEEE, 2012. 1697-1704 [12] Steinbring J, Hanebeck U D. LRKF revisited:the smart sampling Kalman filter (S2KF). Journal of Advances in Information Fusion, 2014, 9(2):106-123 [13] Chlebek C, Hanebeck U D. Bayesian approach to direct pole estimation. In: Proceedings of the 2014 European Control Conference. Strasbourg, France: IEEE, 2014. 1061-1068 [14] Arasaratnam I, Haykin S. Cubature Kalman filters. IEEE Transactions on Automatic Control, 2009, 54(6):1254-1269 doi: 10.1109/TAC.2009.2019800 [15] Zhang L, Li S, Zhang E Z, Chen Q W. Robust measure of non-linearity-based cubature Kalman filter. IET Science, Measurement & Technology, 2017, 11(7):929-938 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=ac386da09116b4daa58363316e178f0f [16] Huang Y L, Zhang Y G, Li N, Zhao L. Gaussian approximate filter with progressive measurement update. In: Proceedings of the 2015 IEEE 54th Annual Conference on Decision and Control. Osaka, Japan: IEEE, 2016. 4344-4349 [17] Boutayeb M, Rafaralahy H, Darouach M. Convergence analysis of the extended Kalman filter used as an observer for nonlinear deterministic discrete-time systems. IEEE Transactions on Automatic Control, 1997, 42(4):581-586 doi: 10.1109/9.566674 [18] Boutayeb M, Aubry D. A strong tracking extended Kalman observer for nonlinear discrete-time systems. IEEE Transactions on Automatic Control, 1999, 44(8):1550-1556 doi: 10.1109/9.780419 [19] Xiong K, Zhang H Y, Chan C W. Performance evaluation of UKF-based nonlinear filtering. Automatica, 2006, 42(2):261-270 doi: 10.1016/j.automatica.2005.10.004 [20] Chang G B. Robust Kalman filtering based on Mahalanobis distance as outlier judging criterion. Journal of Geodesy, 2014, 88(4):391-401 doi: 10.1007/s00190-013-0690-8 [21] Gibbs R G. New Kalman filter and smoother consistency tests. Automatica, 2013, 49(10):3141-3144 doi: 10.1016/j.automatica.2013.07.013 [22] Huang S D, Dissanayake G. Convergence and consistency analysis for extended Kalman filter Based SLAM. IEEE Transactions on Robotics, 2007, 23(5):1036-1049 doi: 10.1109/TRO.2007.903811 [23] Uhlmann J K. Covariance consistency methods for fault-tolerant distributed data fusion. Information Fusion, 2003, 4(3):201-215 doi: 10.1016/S1566-2535(03)00036-8 [24] Tarn T J, Rasis Y. Observers for nonlinear stochastic systems. IEEE Transactions on Automatic Control, 1976, 21(4):441-448 doi: 10.1109/TAC.1976.1101300 期刊类型引用(9)
1. 杨旭升,李唯诣,张文安. 面向RTK定位的整数约束型渐进高斯滤波方法. 自动化学报. 2025(02): 366-375 . 本站查看
2. 杨旭升,吴江宇,胡佛,张文安. 基于渐进高斯滤波融合的多视角人体姿态估计. 自动化学报. 2024(03): 607-616 . 本站查看
3. 张文安,沈嘉俊,史秀纺,杨旭升,王军. 基于自适应高斯渐进滤波的工程车GNSS/INS紧组合定位. 传感技术学报. 2024(04): 620-628 . 百度学术
4. 杨旭升,李福祥,胡佛,张文安. 基于肌电-惯性融合的人体运动估计:高斯滤波网络方法. 自动化学报. 2024(05): 991-1000 . 本站查看
5. 张文安,林安迪,杨旭升,俞立,杨小牛. 融合深度学习的贝叶斯滤波综述. 自动化学报. 2024(08): 1502-1516 . 本站查看
6. 杨旭升,王雪儿,汪鹏君,张文安. 基于渐进无迹卡尔曼滤波网络的人体肢体运动估计. 自动化学报. 2023(08): 1723-1731 . 本站查看
7. 徐瑞龙,石琳. 融合凝聚耦合度与航迹成长值的船舶跟踪算法. 计算机与数字工程. 2021(08): 1574-1579+1665 . 百度学术
8. 陈立军,郝保雷,郑剑. AUV平台的水下机动目标跟踪. 网络新媒体技术. 2020(05): 48-54 . 百度学术
9. 徐从裕,徐俊,高雨婷,胡宗久,杨雅茹. 面向在线测量的亚像素提取与实验验证. 电子测量与仪器学报. 2019(08): 23-29 . 百度学术
其他类型引用(7)
-