一种用于异步多基地雷达系统的分布式批估计融合方法与流程

本发明涉及雷达数据融合技术领域，特别涉及异步多基地雷达数据融合技术。

背景技术：

近年来，由于电子对抗环境复杂程度的提高，多基地雷达得到了越来越多的关注，尤其是现代通信技术和多雷达数据融合技术的发展，为多基地雷达的发展提供了良好的契机。与普通单基地雷达相比，多基地雷达具有作用距离远、抗低空突防、隐蔽性和抗干扰性强等众多优点，因此广泛地应用到了目标检测、目标定位和目标跟踪等众多领域。其中，利用多基地雷达对多目标进行跟踪作为一个极为重要的研究课题得到了十分深入的研究。

对于多目标跟踪应用，采用发射天线宽波束泛光照射，接收天线窄波束扫描(宽发窄收)的工作模式可以较好地克服雷达覆盖率与分辨率之间存在矛盾的问题。但是，当多基地雷达工作在该模式下，由于各天线初始指向、扫描周期不同，目标运动等原因，多基地雷达间的数据通常是异步的。因此，利用多雷达数据融合改善目标跟踪精度，降低虚警概率，提高航迹完整度面临异步数据融合问题。文献“targettrackingwithasynchronousmeasurementsbyanetworkofdistributedmobileagents,”inproc.ieeeicassp,3857-3860，mar2012”提出了一种先假设多传感器数据同步，再利用同步数据融合方法融合多传感器数据的错误同步估计融合方法，但该方法没有考虑多雷达数据的异步特性，融合精度较低，仅适用于多基地雷达数据异步程度较弱的情况，特别当异步程度较强时，目标甚至无法被成功跟踪，航迹丢失现象十分严重；文献“distributedsequentialestimationinasynchronouswirelesssensornetworks，ieeesignalprocessletters.spl-22，1965-1969，nov2015”提出了一种首先外推对齐异步数据，然后融合对齐后数据的分布式序贯估计融合方法，它有效地解决了异步数据融合的问题，但该方法需要在每个数据到达的时刻都进行递推和更新操作，当多基地雷达数据到来十分频繁时，系统消耗的计算资源通常难以承受，因此该方法并不适用于数据更新速度快且计算资源有限的多基地雷达系统。

技术实现要素：

本发明的发明目的在于：针对上述存在的问题，提供一种用于异步多基地雷达系统的分布式批估计数据融合方法，以有效地解决异步多基地雷达系统中由于采样速率和初始偏差不同导致的异步数据难以融合的技术问题。

本发明的用于异步多基地雷达系统的分布式批估计数据融合方法，包括下列步骤：

步骤1：设置批估计更新周期，并初始化当前迭代时间t＝0；可以根据本地雷达的采样速率或数据更新的实际需求对批估计更新周期进行设置；

步骤2：判断t值是否等于本地雷达接收到当前量测的时间；若是则执行步骤4；否则执行步骤3；

步骤3：在一个批估计更新周期内，依次接收和存储当前迭代对应的其他雷达(非本地雷达)的高斯参数，并更新迭代时间t＝t+δt，δt为预设设定的时间间隔，然后返回步骤2；

优选的，将时间间隔δt设置为：δt＝批估计更新周期/p，其中p为整数。

步骤4：从本地雷达读取当前量测，并利用粒子滤波算法进行本地滤波，将滤波输出作为本地雷达当前组粒子样本的后验估计结果和本地预测密度函数，然后执行步骤5；

步骤5：采用高斯近似模型将步骤4得到的本地雷达当前组粒子样本的后验估计结果和本地预测密度函数进行参数化表征，得到高斯参数；然后将本地雷达的高斯参数向外发送，使得非本地雷达可以接收到该高斯参数；为了达到减少通信量的目的，用高斯参数近似的表征本地雷达当前组粒子样本的后验概率密度函数和预测的概率密度；

步骤6：将步骤3存储的非本地雷达的高斯参数和步骤5得到的本地雷达的高斯参数进行批估计数据融合，得到本地雷达当前组粒子样本的权值；

步骤7：根据本地雷达当前组粒子样本的权值输出目标状态，并执行步骤8；

步骤8：执行k＝k+1，然后判断k是否小于或等于k，若是则结束；否则执行步骤2；其中k表示本地雷达在系统观测总时间内接收到的总的量测个数，量测数k的初始值为1。

进一步的，步骤6中，在批估计数据融合时，分别对所有本地雷达的高斯参数进行利用目标状态转移特性进行递推的对齐处理，得到对齐结果；再将当前组粒子样本的各对齐结果的连乘作为本地雷达当前组粒子样本的更新参量其中上标(·)表示粒子样本标识符，n表示本地雷达的雷达标识符；基于更新参数得到本地雷达当前组粒子样本的权值其中ω^(q)(tn,k-1)表示对应上一个量测的本地雷达当前组粒子样本的权值；

其中，利用目标状态转移特性进行递推的对齐处理为：

对本地雷达的高斯参数的对齐处理为：由得到本地雷达的对齐结果，其中表示对应粒子样本的后验估计结果的高斯参数，μn,k、pn,k分别表示对应粒子样本的后验估计结果的均值、方差，表示对应粒子样本的预测密度函数的高斯参数，vn,k、rn,k分别表示对应粒子样本的预测密度函数的均值、方差，下标n为本地雷达标识符，k为雷达n的当前量测标识符；

对非本地雷达的高斯参数的对齐处理为：

由得到非本地雷达的对齐结果p(z(ξj,k)|x^(q)(tn,k-1),x^(q)(tn,k))；

其中ξj,k表示非本地雷达的高斯参数的到达时间，下标j表示到达时间ξj,k按照到达顺序的编号，ξj,k的第二个下标为量测标识符；和表示对应到达时间ξj,k的高斯参数，表示对应粒子样本的后验估计结果高斯参数，μi,l、pi,l分别表示对应粒子样本的后验估计结果的均值、方差，表示对应粒子样本的预测密度函数的高斯参数，vn,k、ri,l分别表示对应粒子样本的预测密度函数的均值、方差，下标i用于区分不同的非本地雷达，l为非本地雷达i的当前量测标识符；

其中，p(x(ξj,k)|x^(q)(tn,k-1),x^(q)(tn,k))表示服从均值为e(x(ξj,k)|x(tn,k-1),x(tn,k))、方差为cov(x(ξj,k)|x(tn,k-1),x(tn,k))的高斯分布，均值e(x(ξj,k)|x(tn,k-1),x(tn,k))和方差cov(x(ξj,k)|x(tn,k-1)，x(tn,k))具体为：

e(x(ξj,k)|x(tn,k-1)，x(tn,k))＝e(x(ξj,k)|x(tn,k-1))+ψs^-1(x(tn,k)-e(x(tn,k)|x(tn,k-1)))；

cov(x(ξj,k)|x(tn,k-1),x(tn,k))＝cov(x(ξj,k)|x(tn,k-1))-ψs^-1ψ′；

s＝cov(x(tn,k)|x(tn,k-1))，ψ＝cov(x(ξj,k),x(tn,k)|x(tn,k-1))，

参量e(x(ξj,k)|x(tn,k-1))、e(x(tn,k)|x(tn,k-1))、cov(x(ξj,k)|x(tn,k-1))、s、ψ分别通过目标状态转移方程计算得到；

x(ξj,k)、x(tn,k-1)、x(tn,k)分别表示在时刻ξj,k、tn,k-1、tn,k的目标状态。

发明采用分布式批估计融合方法，首先根据本地雷达的采样速率或数据更新的实际需求等预设批估计更新周期，并利用粒子滤波算法在多基地雷达进行本地滤波，获得多基地雷达在更新周期内的局部后验估计结果与本地预测密度函数，然后采用高斯近似模型将局部后验估计结果和本地预测密度函数近似为高斯分布，得到高斯分布参数(简称高斯参数)，并在多雷达间通信交互对应的高斯参数，最后通过引入一种利用目标状态转移特性进行递推的对齐策略将多雷达的异步高斯参数进行对齐，并采用分布式批估计融合方法融合对齐后的高斯参数，从而有效地解决了多基地雷达系统在宽发窄收模式下由于各天线初始指向不同，各天线扫描周期不同，目标运动导致的异步数据难以融合的技术问题，并带来如下技术效果：

(1)本发明的方法操作简单，只需根据一个简单的对齐策略将异步数据同步即可；

(2)与现有的一些序贯估计方法相比，本发明的方法不需要频繁地进行递推和更新操作，因此计算量较低；

(3)本发明的方法是对一段时间内的数据进行处理的(批处理)，在这一段时间内的异步数据是怎样的时序对于融合精度的影响较小，因此本发明的方法可以对任意时序的量测都有较高的融合精度，鲁棒性较强。

(4)本发明将多基地雷达的本地后验估计结果和本地预测密度函数被近似为高斯分布，之后通过在多雷达间传输高斯参数代替传输高维的粒子样本，因此具有通信量小的优点。

附图说明

图1为本发明的流程图。

图2为本发明实施例提供的多基地雷达场景图。

图3为本发明与错误同步估计融合方法、序贯估计融合方法、真实同步估计融合方法以及单接收机雷达的跟踪精度对比图。

图4为本发明与错误同步估计融合方法以及真实同步融合方法在多个场景下的跟踪精度对比图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面结合实施方式和附图，对本发明作进一步地详细描述。

参见图1，本发明的用于异步多基地雷达系统的分布式批估计数据融合方法具体实现过程如下：

首先初始化系统参数，包括：

监测平面大小；发射机雷达个数m；接收机雷达个数n；发射机雷达m,m＝1,2,...,m的安放位置接收机雷达n,n＝1,2,...,n的安放位置观测总时间ttotal；量测标识符k＝1；t＝0s；目标在t＝0s的初始状态其中(x(0),y(0))目标的初始位置，表示目标的初始速度；目标初始状态偏差服从高斯分布

如图2所示实施例，其对应的监测平面大小为200km×200km；发射机雷达个数m＝2；接收机雷达个数n＝3；发射机雷达的安放位置分别为(120km,170km)和(170km,120km)，接收机雷达的安放位置分别为(180km,80km)，(150km,150km)，和(80km,180km)；观测总时间ttotal＝400s；t＝0s；目标在t＝0s的初始状态x(0)＝(30,0.3,25,0.3)′；目标初始状态偏差服从高斯分布并且c0＝diag(1,0.05,1,0.05)，符号(·)′表示括号中矩阵的转置，下同。

基于实施例2所给出的对应参数，实现本发明的分布式批估计融合处理的步骤为：

步骤1：根据本地雷达的采样速率或数据更新的实际需求等设定批估计更新周期(简称更新周期)。

本实施例中根据本地雷达的采样速率设定更新周期，具体以本地雷达相邻两个量测的间隔为更新周期，即(tn,k-1,tn,k]。

步骤2：判断t值是否等于本地雷达接收到当前量测的时间tn,k；若是则执行步骤4；否则执行步骤3；

其中，tn,k表示本地雷达n接收到第k个量测的时间，第k个量测即为当前量测。1≤n≤n，1≤k≤k，k表示接收机雷达n在ttotal时间内接受到的总的量测个数。初始时t＝0s，因此本申请首先执行步骤3，接收存储其他雷达(非本地雷达)的高斯参数；然后再读取本地量测。

步骤3：在一个批估计更新周期内，依次接收和存储当前迭代时间t对应的其他雷达的高斯参数、高斯参数的到达时间，并更新t值，然后返回步骤2；

在一个批估计更新周期内，即(tn,k-1,tn,k]时间内，tn,k-1表示本地雷达n接收到第k-1个量测的时间，本地雷达n共接收到nk-1组其他雷达发送来的用于近似的表征后验估计结果和预测密度函数的高斯参数，由于后验估计结果和预测密度函数产生的时间与他们对应的量测产生的时间相同；由于量测异步，因此所接收到的高斯参数也是异步的，即存在时间先后顺序；按照其到达本地雷达n的时间先后顺序有：j＝1,2,...,nk-1；其中ξj,k表示接收机雷达n的第j组高斯参数的到达时间，另外，由于时间ξj,k事实上也是接收机雷达i(生成第j组高斯参数的雷达)接收到第l个量测的时间ti,l(即第l个本地量测到达接收机雷达i的到达时间)，因此ξj,k＝ti,l，其中n,i∈{1,2,...,n}且n≠i，1≤l≤k。

所述更新t值具体为：t＝t+δt，δt为设定的时间间隔；δt＝批估计更新周期/p，p为整数。

步骤4：从本地雷达读取当前量测z(tn,k)，并利用粒子滤波算法进行本地滤波，将输出结果作为本地雷达当前组粒子样本的后验估计结果和预测密度函数，然后执行步骤5；

其中，本地雷达的量测为z(tn,k)＝[r(tn,k),θ(tn,k)]′，r(tn,k)和θ(tn,k)分别表示在tn,k时刻的距离和角度量测。x(ti,l)表示在ti,l时刻的目标状态；x^(q)(tn,k)表示目标状态x(tn,k)的第q个粒子样本，ω^(q)(tn,k)表示粒子样本对应的权值，q＝1,2,...,q，q表示总的粒子样本数量；

若tn,k＝0s，对于粒子样本有：

ω^(q)(tn,k)＝1/q；

其中，～表示采样，即从高斯分布中采样出x^(q)(ti,l)；

若tn,k>0s，对于粒子样本有：

x^(q)(tn，k)～p(x(tn，k)|x(tn，k-1))

ω^(q)(tn，k)＝ω^(q)(tn，k-1)·p(z(tn，k)|x^(q)(tn，k))

其中，p(x(tn,k)|x(tn,k-1))表示状态转移概率密度函数，即时刻tn,k-1的目标状态x(tn,k-1)转移到时刻tn,k的目标状态x(tn,k)的状态转移概率密度函数，p(x(tn,k)|x(tn,k-1))可以通过目标的状态转移方程可得。p(z(tn,k)|x^(q)(tn,k))表示似然函数，即目标状态为x^(q)(tn,k)时，量测z(tn,k)的条件概率密度函数，p(z(tn,k)|x^(q)(tn,k))可以通过观测方程得到。

最后获得本地雷达当前组粒子样本近似的后验估计结果、预测密度函数：

其中δ(·)表示标准的狄利赫利函数。

在实际处理时，直接将粒子滤波的输出结果作为本地雷达当前组粒子样本的后验估计结果、预测密度函数并执行步骤5。

步骤5：采用高斯近似模型将步骤4得到的本地雷达当前组粒子样本的后验估计结果p(x^(q)(tn,k)|z(tn,1:k))、预测密度函数p(x^(q)(tn,k)|z(tn,1:k-1))进行参数化表征，得到高斯参数；然后将本地雷达的高斯参数发送至其他雷达；

利用高斯近似模型：

基于上述高斯近似模型，将本地雷达当前组粒子样本的后验估计结果、预测密度函数近似为高斯分布和然后将本地雷达的高斯参数和向外发送，使得其他雷达可以接收到该高斯参数。

步骤6：利用步骤3存储的其他雷达的高斯参数及到达时间和步骤5得到的本地雷达的高斯参数进行批估计数据融合，得到本地雷达当前组粒子样本的权值：

引入一种利用目标状态转移特性进行递推的对齐策略将多雷达的异步高斯参数进行对齐。

多基地雷达异步数据估计后验概率密度函数p(x(tn,k)|z(tn,1:k))包含预测与更新两个步骤。其中，

预测步骤可以表示为：p(x(tn,k)|z(tn,1:k))＝∫p(x(tn,k)|z(tn,1:k))p(x(tn,k)|x(tn,k-1))dx(tn,k-1)；

更新步骤可以表示为：p(x(tn,k)|z(tn,1:k))∝p(zn,k|x(tn,k),zn,1:k-1)p(x(tn,k)|z(tn,1:k-1))；其中，∝表示正比于。考虑到p(zn,k|x(tn,k),zn,1:k-1)与同步估计方法中的似然函数相似，本发明中，将其称为异步似然函数。

由于量测zn,k不仅仅取决于目标状态x(tn,k)，还取决于中间目标状态x(ξj,k)，其中ξj,k＝ti,l。因此，在一阶马尔可夫的假设下，异步似然函数可以表达为：p(zn,k|x(tn,k),zn,1:k-1)＝p(zn,k|x(tn,k-1),x(tn,k))。同时，在多雷达量测相互独立的条件下，异步似然函数可以进一步分解为：

其中因子由本地雷达产生，p(z(ξj,k)|x(tn,k-1),x(tn,k))不能通过滤波方法直接得到，因此根据条件概率的基本性质，该因子进而可以写成：

p(z(ξj,k)|x(tn,k-1),x(tn,k))＝∫p(z(ξj,k)|x(ξj,k))p(x(ξj,k)|x(tn,k-1),x(tn,k))dx(ξj,k)

其中，p(z(ξj,k)|x(ξj,k))是接收雷达i的本地似然函数，而p(x(ξj,k)|x(tn,k-1),x(tn,k))是服从均值为e(x(ξj,k)|x(tn,k-1),x(tn,k))、方差为cov(x(ξj,k)|x(tn,k-1),x(tn,k))的高斯分布。

因此，本发明的批估计融合准则可以表示为：

此外，本地似然还可以进一步写成：其中，p(x(ξj,k)|z(ξj,k),zi,1:(l-1))是本地后验概率密度函数，p(x(ξj,k)|zi,1:(l-1))是根据先验结果预测得到的概率密度函数。则基于上述对本地似然函数的变形，本发明的批估计融合准则可以表示为：

可以看到，其引入了一种根据目标状态转移特性将本地滤波结果递推对齐的策略，即p(x(ξj,k)|x(tn,k-1),x(tn,k))，并联合nk个雷达的本地滤波结果p(x(ξj,k)|z(ξj,k),zi,1:(l-1))/p(x(ξj,k)|zi,1:(l-1))对目标状态x(tn,k)进行估计。

而在本发明中，雷达间相互间是通过交互高斯参数和)来实现对本地的后验估计结果、预测密度函数的交互，因此对于任意接收雷达n，当其联合nk个雷达的本地滤波结果对目标状态x(tn,k)进行估计时，首先通过利用目标状态转移特性进行递推的对齐策略分别将非本地雷达、本地雷达的高斯参数进行对齐，获取对齐结果：

对于本地雷达，直接根据得到本地雷达的对齐结果；

对于非本地雷达，根据下式计算对齐后结果(p(z(ξj,k)|x^(q)(tn,k-1),x^(q)(tn,k)))：

其中，p(x(ξj，k)|x^(q)(tn，k-1)，x^(q)(tn，k))表示服从均值为e(x(ξj，k)|x(tn，k-1)，x(tn，k))、方差为cov(x(ξj，k)|x(tn，k-1),x(tn，k))的高斯分布。均值e(x(ξj，k)|x(tn，k-1)，x(tn，k))和方差cov(x(ξj,k)|x(tn,k-1)，x(tn,k))具体为：

e(x(ξj,k)|x(tn,k-1)，x(tn,k))＝e(x(ξj,k)|x(tn,k-1))+ψs^-1(x(tn,k)-e(x(tn,k)|x(tn,k-1)))；

cov(x(ξj,k)|x(tn,k-1),x(tn,k))＝cov(x(ξj,k)|x(tn,k-1))-ψs^-1ψ′；

s＝cov(x(tn,k)|x(tn,k-1))，ψ＝cov(x(ξj,k),x(tn,k)|x(tn,k-1))，

参量e(x(ξj,k)|x(tn,k-1))、e(x(tn,k)|x(tn,k-1))、cov(x(ξj,k)|x(tn,k-1))、s、ψ分别通过目标状态转移方程计算得到。

然后，基于nk-1个异步似然函数p(z(ξj,k)|x^(q)(tn,k-1),x^(q)(tn,k))、以及本地似然函数(即p(x^(q)(tn,k)|z(tn,1:k)))的连乘对本地雷达当前组粒子样本的权值进行更新，得到更新后的当前组粒子样本的权值ω^(q)(tn,k)：

其中，当j＝nk时，

步骤7：根据步骤6得到的本地雷达当前组粒子样本的权值输出目标状态；然后执行步骤8；其中tn,k时刻目标的状态为

步骤8：执行k＝k+1，然后判断k是否小于或等于k，若是则结束；否则执行步骤2；

图3给出了不同估计方法融合精度的对比结果。从图3可以看到，由于分布式批估计方法有效地融合了多雷达的异步数据，因此获得了相比单接收机雷达更高的跟踪精度。同时，分布式批估计方法的均方根误差曲线距离理想的真实同步估计方法曲线非常近，这也反映了该方法跟踪精度较高的优点。除此之外，从图3还可以看到，由于错误同步方法简单地将异步数据进行同步融合，导致了其结果并不收敛。事实上，它的效果甚至比单部雷达的跟踪精度更低。而与之相反，由于分布式批估计方法采用了本发明提出的“根据状态转移特性设计的对齐策略”将多雷达数据进行了对齐，考虑了多基地雷达数据的异步特性，因此获得了很好的跟踪效果，其融合精度相比错误同步估计方法要高得多。

图4给出了不同异步程度下不同方法的融合精度。其中，对于三个场景的异步程度有场景1>场景2>场景3。从图4可以看到，无论异步程度如何变化，分布式批估计融合方法的跟踪精度始终比错误同步融合方法更高。并且还可以看出，多雷达的异步程度越低，跟踪精度也就越高。最为重要的是，错误同步估计融合方法的跟踪精度受到了异步程度很大的影响。相反，由于分布式批估计融合方法的批处理特性，它的均方根误差曲线在不同异步程度的场景下基本没有波动，受到的影响极小。因此，我们可以得出这样的结论，即相比错误同步估计方法，分布式批估计方法还具有鲁棒性强的优点。

不同估计方法的计算消耗时间的对比结果如表1所示：

表1

当将仿真运行时间作为评估标准。从表1可以看出，不论接收机雷达的个数为多少，真实同步估计方法的计算消耗始终最小。但是在n＝3时，异步序贯估计方法的计算消耗比分布式批估计方法更低。而当n＝8,15,25时，后者的计算消耗则更低。这是因为，对于异步序贯估计方法，当接收机雷达数量较少时，在更新周期内需要进行的递推和更新操作并不频繁，它们对于计算消耗的影响还不明显。相反，由于分布式批估计方法需要在固定时刻进行积分近似、连乘等操作，它在单个时刻需要消耗的计算资源甚至比异步序贯估计方法进行几次递推和更新操作消耗的计算资源更多，因此在雷达数量较少时，可能会出现异步序贯估计方法的计算消耗比分布式批估计方法更低的结果。但是，当接收机雷达的数量一旦增加，如n＝8,15,25时，更新周期内的量测数量也会大大增加，此时过于频繁的递推和更新操作导致异步序贯估计方法计算消耗远远大于分布式批估计方法单个时刻需要消耗的计算资源。因此我们可以得出这样的结论，即在雷达数量较多，数据更新速度快且计算资源有限的异步多基地雷达场景下，相比异步序贯估计方法，分布式批估计方法存在计算量小的优势。

以上所述，仅为本发明的具体实施方式，本说明书中所公开的任一特征，除非特别叙述，均可被其他等效或具有类似目的的替代特征加以替换；所公开的所有特征、或所有方法或过程中的步骤，除了互相排斥的特征和/或步骤以外，均可以任何方式组合。