本发明属于多传感器数据融合技术领域，具体涉及一种非线性系统任意延迟步数的单个失序量测集中式融合方法。
背景技术：
利用多个传感器协同工作可以提高系统可靠性和鲁棒性，扩展时间和空间覆盖率，改善目标探测性能，随着信息处理和网络通信等技术的迅速发展，多传感器数据融合技术受到了广泛关注。其中，在集中式数据融合方法中，各本地传感器将局部量测送至融合中心进行处理，不存在信息损失，可以获得最优的全局估计值，因而在区域监视、目标跟踪、目标检测和定位等众多领域得到了广泛的应用。
集中式融合系统中，由于各传感器的数据处理时间和通信传输时间不同，到达融合中心的各传感器局部量测存在随机时延，导致来自同一目标的量测失序，即产生失序量测问题。一种常用的解决方法是忽略失序量测，但是这种方法通常会降低跟踪性能，尤其当失序量测频繁出现或者失序量测来自精度更高的局部传感器时，直接忽略失序量测将严重影响融合精度。文献“一种最佳多延迟无序量测处理算法，系统工程与电子技术，2010，31(11)：2592-2596”基于最优线性固定点平滑思想，直接使用失序量测更新当前时刻的状态估计，有效地解决了多步延迟的失序量测难以更新融合的问题，但是该方法基于卡尔曼滤波算法，仅适用于线性或非线性强度很弱的动态系统；文献“Out-of-sequence measurement processing for particle filter：exact bayesian solution，IEEE Trans.Aersp.Electron.Syst.，vol.48，no.4，pp.2818-2831，2012”提出了一种基于确定贝叶斯方法的粒子滤波算法(A-PF算法)，该方法适用于一般地非线性系统，并且当粒子数目足够大时该方法近似最优，跟踪精度高，然而由于该方法中使用粒子平滑方法计算平滑密度函数，其计算量过大，存储要求过高，仅适用于低维或非实时场景。
技术实现要素：
本发明的发明目的是：为了解决现有技术中存在的以上问题，本发明提出了一种非线性系统任意延迟步数的单个失序量测集中式融合方法。
本发明的技术方案是：一种非线性系统任意延迟步数的单个失序量测集中式融合方法，其特征在于，包括以下步骤：
A、对融合中心进行初始化，并初始化粒子样本及其相应的权值、初始化更新时间tk＝0，设初始更新帧数k＝0；
B、新的局部量测到达融合中心，判断当前量测的时戳tz是否小于或等于系统观测总时间ttotal；若是，则执行步骤C；若否，则操作结束；
C、判断当前量测的时戳tz是否大于或等于上一个更新时间tk；若是，则执行步骤D；若否，则执行步骤E；
D、利用粒子滤波算法获得k+1时刻的顺序后验权值，得到k+1时刻的目标状态估计值，并将以粒子样本表征的后验状态近似为高斯参数进行存储，再返回步骤B；
E、采用广义高斯平滑方法求解失序量测产生时刻的平滑概率密度函数；
F、以tz时刻的平滑概率密度函数作为重要性采样函数进行二次采样，联合局部似然函数求解异步似然函数；
G、利用异步似然函数更新k时刻的后验权值，再通过重采样得到k时刻更新后以粒子样本表征的融合后验状态，得到k时刻的目标状态估计值，并将以粒子样本表征的后验状态近似为高斯参数进行存储，再返回步骤B。
进一步地，所述步骤A还包括初始化系统参数，所述系统参数包括监测平面大小、传感器总个数、观测总时间以及目标初始状态。
进一步地，所述步骤D中k+1时刻的目标状态估计值表示为：
其中，表示k+1时刻的目标状态估计值，表示k+1时刻第q个粒子样本的目标状态，表示k+1时刻第q个粒子样本对应的顺序后验权值，q＝1,2,…,Q表示粒子样本标号，Q为粒子数目。
进一步地，所述步骤E采用广义高斯平滑方法求解失序量测产生时刻的平滑概率密度函数，具体包括以下分步骤：
E1、根据融合中心第k个更新时刻对应的更新时间tk，通过条件tk-l＜tz＜tk-l+1确定失序量测的延迟步数l，判断延迟步数l是否等于1；若延迟步数l＝1，则执行步骤E2；若否，则执行步骤E3；
E2、利用广义高斯平滑方法求解tz时刻的平滑概率密度函数；
E3、利用广义高斯平滑方法求解中间统计量，通过对中间统计量求解联合积分，获得tz时刻的平滑概率密度函数。
进一步地，所述步骤E2中平滑概率密度函数表示为：
其中，τ表示tz时刻，x(τ)表示tz时刻的目标状态，分别表示平滑概率密度的均值和协方差。
进一步地，所述步骤E3利用广义高斯平滑方法求解中间统计量，通过对中间统计量求解联合积分，获得tz时刻的平滑概率密度函数，具体包括以下分步骤：
E31、进行迭代初始化，以时刻k-1代替时刻τ，得到k-1时刻的高斯近似的平滑密度函数，表示为
E32、当延迟步数l＞2时进行迭代计算，令n＝k-2,…,k-l+1，得到n时刻的平滑概率密度，表示为
迭代完成后，获得k-l+1时刻的平滑概率密度函数，表示为
E33、利用广义高斯平滑方法计算τ时刻的一步平滑概率密度函数，表示为
E34、对步骤E32中k-l+1时刻的平滑概率密度函数和步骤E33中τ时刻的一步平滑概率密度函数求解联合积分，表示为
p(x(τ)|xk,z1:k)＝∫p(x(τ)|xk-l+1,z1:k-l)p(xk-l+1|xk,z1:k)dxk-l+1
得到tz时刻的平滑概率密度函数。
进一步地，所述步骤F中异步似然函数表示为
其中，z(τ)表示τ时刻的量测，j＝1,2,…,S表示二次采样得到的粒子样本标号，S表示二次采样的粒子总数。
进一步地，所述步骤G中利用异步似然函数更新k时刻的后验权值具体为：
进一步地，所述步骤G中k时刻的目标状态估计值表示为
本发明的有益效果：本发明采用广义高斯平滑方法求解任意非线性系统的平滑概率密度函数，并利用粒子滤波算法实现更新，首先判断量测是否是顺序量测，若是则利用粒子滤波算法获得顺序后验状态，否则进行失序量测融合，即先利用广义高斯平滑方法求解失序量测产生时刻的平滑概率密度函数，然后联合平滑概率密度函数和局部似然函数求解异步似然函数，最后利用异步似然函数更新上一时刻的顺序后验状态，获得融合后验状态，并将以粒子样本表征的融合后验状态近似为高斯参数进行存储，有效地解决了异步传感器网络中由于数据处理和通信延迟等原因导致的失序量测数据难以融合的问题；并且具备以下优点：
(1)本发明的方法操作简单，只需要根据失序量测的延迟周期利用广义高斯平滑方法求解平滑概率密度函数，进而得到异步似然函数；
(2)与现有的一些非线性系统下失序量测数据融合方法相比，本发明的方法从确定地贝叶斯方法出发，理论最优，融合精度高；
(3)本发明的方法适用于任意延迟周期的失序量测；
(4)与利用粒子平滑方法计算平滑概率密度的最优粒子滤波算法(A-PF)相比，本发明的方法通过高斯近似，采用广义高斯平滑方法计算平滑密度函数，大大减小了计算复杂度；
(5)本发明的方法存储后验状态的高斯参数，相比于存储表征后验状态的粒子及相应权值，存储要求更低。
附图说明
图1是本发明的非线性系统任意延迟步数的单个失序量测集中式融合方法的流程示意图；
图2是本发明实施例中失序量测延迟一步的时序示意图；
图3是本发明实施例中失序量测延迟两步的时序示意图；
图4是本发明实施例中失序量测延迟三步的时序示意图；
图5是本发明实施例中传感器网络场景与目标轨迹示意图；
图6是本发明实施例中失序量测延迟一步的场景下的跟踪精度对比图；
图7是本发明实施例中失序量测延迟二步的场景下的跟踪精度对比图；
图8是本发明实施例中失序量测延迟三步的场景下的跟踪精度对比图；
图9是本发明实施例中失序量测延迟一步、两步、三步场景下的计算效率提升倍数对比图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。
如图1所示，为本发明的非线性系统任意延迟步数的单个失序量测集中式融合方法的流程示意图。一种非线性系统任意延迟步数的单个失序量测集中式融合方法，包括以下步骤：
A、对融合中心进行初始化，并初始化粒子样本及其相应的权值、初始化更新时间tk＝0，设初始更新帧数k＝0；
B、新的局部量测到达融合中心，判断当前量测的时戳tz是否小于或等于系统观测总时间ttotal；若是，则执行步骤C；若否，则操作结束；
C、判断当前量测的时戳tz是否大于或等于上一个更新时间tk；若是，则执行步骤D；若否，则执行步骤E；
D、利用粒子滤波算法获得k+1时刻的顺序后验权值，得到k+1时刻的目标状态估计值，并将以粒子样本表征的后验状态近似为高斯参数进行存储，再返回步骤B；
E、采用广义高斯平滑方法求解失序量测产生时刻的平滑概率密度函数；
F、以tz时刻的平滑概率密度函数作为重要性采样函数进行二次采样，联合局部似然函数求解异步似然函数；
G、利用异步似然函数更新k时刻的后验权值，再通过重采样得到k时刻更新后以粒子样本表征的融合后验状态，得到k时刻的目标状态估计值，并将以粒子样本表征的后验状态近似为高斯参数进行存储，再返回步骤B。
在本发明的一个可选实施例中，本发明首先将用于处理非线性系统中失序量测的高效粒子滤波算法初始化系统参数，包括：观测平面大小，传感器个数Ni，观测总时间ttotal，目标的初始状态其中(x(0),y(0))表示目标的初始位置，表示目标的初始速度，目标初始状态偏差服从高斯分布
如图5所示，为本发明实施例中传感器网络场景与目标轨迹示意图，其中N＝2部传感器对一个观测平面大小为200km×200km的二维平面区域一运动目标进行监视，目标的初始状态x(0)＝(200,-1,200,-1.2)′，即目标的初始位置是(200,200)，并以(-1,-1.2)的速度运动，图5给出了观测总时间ttotal＝35s的跟踪结果。目标初始状态偏差服从高斯分布并且C0＝diag(1,0.01,1,0.01)。
本发明设定传感器2由于系统数据处理和通信传输原因，所有量测均延迟某一周期数；本实施例设置了三个延迟场景，设在场景1、2、3中，传感器2的所有局部量测分别延迟一步、两步、三步。如图2、3、4所示，分别为本发明实施例中失序量测延迟一步、两步、三步的时序示意图。
在本发明的一个可选实施例中，上述步骤A对融合中心进行初始化，并初始化粒子样本及其相应的权值、初始化更新时间tk＝0，设初始更新帧数k＝0，产生Q个初始粒子样本服从其相应地权值为其中，q＝1,2,…,Q表示粒子样本标号，Q为粒子数目。
在本发明的一个可选实施例中，上述步骤C判断当前量测的时戳tz是否大于或等于上一个更新时间tk；若当前量测的时戳tz大于或等于上一个更新时间tk，则表明该当前量测是顺序量测，那么执行步骤D，利用粒子滤波算法顺序更新，获得当前时刻的后验状态；若当前量测的时戳tz小于上一个更新时间tk，则表明该当前量测是失序量测，那么执行步骤E，对该失序量测进行融合处理，以校正上个更新时刻的后验状态。
在本发明的一个可选实施例中，上述步骤D利用粒子滤波算法获得k+1时刻的顺序后验状态，得到k+1时刻的目标状态估计值，并将以粒子样本表征的后验状态近似为高斯参数进行存储，再返回步骤B；
已知粒子样本对于所求粒子样本有
其中，表示k时刻第q个粒子样本的目标状态，表示k时刻第q个粒子样本对应的顺序后验权值；～表示采样，即从高斯分布中采样得到粒子样本p(xk+1|xk)表示k时刻的目标状态xk转移到k+1时刻的目标状态xk+1的条件概率密度函数，可以由目标状态转移方程得到；zk+1表示k+1时刻到达融合中心的量测，p(zk+1|xk+1)表示似然函数，可以由量测方程得到。
从而得到k+1时刻的目标状态估计值，表示为：
其中，表示k+1时刻第q个粒子样本对应的顺序后验权值。
在本发明的一个可选实施例中，上述步骤E采用广义高斯平滑方法求解失序量测产生时刻的平滑概率密度函数，具体包括以下分步骤：
E1、根据融合中心第k个更新时刻对应的更新时间tk，通过条件tk-l＜tz＜tk-l+1确定失序量测的延迟步数l，判断延迟步数l是否等于1；若延迟步数l＝1，则执行步骤E2；若否，则执行步骤E3；
E2、利用广义高斯平滑方法求解tz时刻的平滑概率密度函数；
E3、利用广义高斯平滑方法求解中间统计量，通过对中间统计量求解联合积分，获得tz时刻的平滑概率密度函数。
在本发明的一个可选实施例中，上述步骤E1中对于一步延迟失序量测(l＝1)，利用广义高斯平滑方法直接得到tz时刻的平滑概率密度函数；对于多步延迟失序量测(l≥2)，利用广义高斯平滑方法计算中间统计量，通过中间统计量联合积分求解tz时刻的平滑概率密度函数。
在本发明的一个可选实施例中，上述步骤E2利用广义高斯平滑方法求解tz时刻的平滑概率密度函数，具体为：
对tz时刻的平滑密度进行高斯近似，表示为
其中τ表示tz时刻，x(τ),xτ均表示tz时刻的目标状态，z1:k表示从第一个时刻开始到第k时刻所有量测的集合，分别表示时刻τ对应的平滑概率密度的均值和协方差。
设目标的状态转移方程为
xk＝g(xk-1)+wk,k-1
其中，xk,xk-1分别表示离散k时刻和k-1时刻的目标状态，g(·)表示一般地非线性状态转移函数，wk,k-1表示从k-1时刻到第k时刻积累的过程噪声，服从高斯分布其中Qk,k-1表示从时刻k-1到时刻k积累的过程噪声的协方差矩阵。
将τ时刻的平滑概率密度函数表示为
将视作量测来更新预测概率密度p(xτ|z1:k-1)，gk|τ(·)表示从时刻τ到时刻k的目标状态转移函数。
从而通过利用广义高斯平滑方法计算
其中，m表示概率密度函数的均值，其上下标表示不同概率密度函数的均值，表示τ时刻预测概率密度函数的均值，表示τ时刻平滑概率密度函数的均值；P表示协方差矩阵，其上下标表示不同的协方差矩阵，表示τ时刻预测概率密度函数的协方差矩阵，表示τ时刻平滑概率密度函数的协方差矩阵，这些协方差矩阵与其带相同上下标的均值作为一组参数可以刻画一个概率密度函数；E[·]表示均值运算，即对括号内的概率密度函数求均值；Qτ|k-l表示从时刻k-l到时刻τ积累的过程噪声的协方差矩阵；表示τ时刻第q个粒子对应的权值；和xτ均表示τ时刻的目标状态，而xτ是一个随机变量，是第q个粒子的状态向量，是一个确定的向量，在Q个粒子样本，随机变量xτ可能取到这些中的任意向量值，而取到哪一个的概率大小就由来刻画；Dk表示交互协方差矩阵；Gτ表示一个增益矩阵；上标“'”表示矩阵的转置。
得到tz时刻平滑概率密度函数的高斯近似，表示为：
在本发明的一个可选实施例中，上述步骤E3利用广义高斯平滑方法求解中间统计量，通过对中间统计量求解联合积分，获得tz时刻的平滑概率密度函数，具体包括以下分步骤：
E31、进行迭代初始化，以时刻k-1代替步骤E2中的时刻τ，得到k-1时刻的高斯近似的平滑密度函数，表示为
E32、当延迟步数l＞2时进行迭代计算，令n＝k-2,…,k-l+1，
其中，mn表示n时刻顺序后验概率密度函数的均值，表示n时刻一步平滑概率密度函数的均值，Pn表示n时刻顺序后验概率密度函数的协方差矩阵，表示n时刻一步平滑概率密度函数的协方差矩阵。
得到n时刻的平滑概率密度，用高斯参数表示为
迭代完成后，获得k-l+1时刻的平滑概率密度函数，表示为
E33、利用广义高斯平滑方法计算τ时刻的一步平滑概率密度函数，具体为
平滑概率密度函数p(x(τ)|xk-l+1，z1：k-l)可以用高斯参数表示为
E34、对步骤E32中k-l+1时刻的平滑概率密度函数和步骤E33中τ时刻的一步平滑概率密度函数求解联合积分，表示为
p(x(τ)|xk,z1:k)＝∫p(x(τ)|xk-l+1,z1:k-l)p(xk-l+1|xk,z1:k)dxk-l+1
对tz时刻的平滑概率密度函数进行高斯近似，得到tz时刻的平滑概率密度函数，表示为
其中，
在本发明的一个可选实施例中，上述步骤F以tz时刻的平滑概率密度函数作为重要性采样函数进行二次采样，联合局部似然函数求解异步似然函数，具体为
以作为重要性采样函数，进行二次采样，表示为
其中，j＝1,2,…,S，表示二次采样得到的粒子样本标号，S表示二次采样的粒子总数。
联合τ时刻的似然函数p(z(τ)|x(j)(τ))，得到异步似然函数，表示为
其中，z(τ)表示τ时刻的量测，即需要处理的失序量测。
在本发明的一个可选实施例中，上述步骤G利用异步似然函数更新k时刻的顺序后验状态，更新表达式为
其中，表示利用顺序量测更新得到的k时刻的顺序后验权值，表示融合k时刻顺序量测和失序量测进行更新得到的k时刻的融合后验权值。
再通过重采样，得到获得融合了该失序量测后的k时刻的融合后验状态估计，表示为
如图6、7、8所示，分别为失序量测延迟1步、2步、3步场景下忽略失序量测方法、贝叶斯最优实现方法和本发明的跟踪精度对比图。其中，跟踪精度通过估计轨迹与真实轨迹的均方根误差进行描述。由图可以看出，由于忽略失序量测方法直接忽略了失序量测，它有最差的跟踪性能；由于贝叶斯最优实现方法是确定贝叶斯方法近似最优的实现方法，所以有着很好的跟踪性能；而贝叶斯高效实现方法，相比于忽略失序量测方法有更好的跟踪结果，这说明了本发明可以有效地处理失序量测，同时，它的均方根误差曲线相近于贝叶斯最优实现方法的均方根误差曲线，进一步说明了本发明的有效性。
如图9所示，为失序量测延迟一步、两步、三步场景下，本发明在计算效率上相较于贝叶斯最优实现方法提升倍数曲线图，其通过贝叶斯最优实现方法计算时间除以本发明方法计算时间的比值来描述。本发明使用广义高斯平滑方法计算平滑概率密度，以此代替了粒子滤波平滑方法，这在计算效率上相比于贝叶斯最优实现方法非常有效。由图9可以看到，本发明对多步延迟失序量测的处理效率尤为有效，当l＝3时，贝叶斯最优实现方法的计算时间几乎是本发明的180倍，也就是说，本发明相较于贝叶斯最优实现方法节省了将近180倍的计算资源，表明了本发明计算复杂度低，适用于高维或实时场景。
本发明通过综合考虑目标跟踪精度、计算复杂度以及存储需求等问题，提出的非线性系统中任意延迟步数的单个失序量测集中式融合方法相比于忽略失序量测方法和贝叶斯最优实现方法存在优势；相比于直接忽略失序量测的融合精度更高；将后验和平滑密度近似为高斯参数，较存储粒子样本及相应的权值节省了很大的存储空间；相比于采用粒子平滑方法计算平滑概率密度，利用广义高斯平滑方法计算平滑密度的计算量更小。
本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。