Distributed Multi-sensor Multi-view Fusion based on Generalized Covariance Intersection

TL;DR

This paper introduces a distributed multi-sensor fusion method for multi-target tracking that combines Generalized Covariance Intersection with clustering to improve robustness and efficiency in sensors with different fields of view.

Contribution

It proposes a novel fusion algorithm integrating GCI with clustering to handle multi-view sensor data and reduce false targets, with detailed GM implementation and analysis.

Findings

Effective fusion of multi-view sensor data demonstrated

Reduces false target detections compared to existing methods

Robustness confirmed through numerical experiments

Abstract

Distributed multi-target tracking (DMTT) is addressed for sensors having different fields of view (FoVs). The proposed approach is based on the idea of fusing the posterior Probability Hypotheses Densities (PHDs) generated by the sensors on the basis of the local measurements. An efficient and robust distributed fusion algorithm combining the Generalized Covariance Intersection (GCI) rule with a suitable Clustering Algorithm (CA) is proposed. The CA is used to decompose each posterior PHD into well-separated components (clusters). For the commonly detected targets, an efficient parallelized GCI fusion strategy is proposed and analyzed in terms of $L_1$ error. For the remaining targets, a suitable compensation strategy is adopted so as to counteract the GCI sensitivity to independent detections while reducing the occurrence of false targets. Detailed implementation steps using a Gaussian Mixture (GM) representation of the PHDs are provided. Numerical experiments clearly confirms the effectiveness of the proposed CA-GCI fusion algorithm.

Tables

Table 1. Table 1: Notations

${𝐗}_k$

$=$

$\{{𝐱}_k^1,{𝐱}_k^2,\dots ,{𝐱}_k^{n_k}\}$

Table 2. Table 2: Target birth and death

Target	Init. Loc.($m$)	Init. Velo.($m/s$)	Birth/Death($s$)
$T_1$	$[1000,400]$	$[-14,0]$	$1/80$
$T_2$	$[1250,400]$	$[-4,-2.5]$	$1/80$
$T_3$	$[500,100]$	$[-8,10]$	$10/60$
$T_4$	$[0,600]$	$[0,-4]$	$10/80$
$T_5$	$[1000,200]$	$[-9,9]$	$10/70$
$T_6$	$[1250,505]$	$[-14,-7]$	$20/60$
$T_7$	$[1000,600]$	$[-12,-7]$	$20/60$
$T_8$	$[250,200]$	$[8,10]$	$20/70$
$T_9$	$[1250,300]$	$[-16,0]$	$30/70$
$T_{10}$	$[-150,500]$	$[32,0]$	$30/70$
$T_{11}$	$[400,600]$	$[12,3]$	$40/80$

Table 3. Table 3: OSPA errors with different detection probabilities and same clutter rate ( λ c = 20 subscript 𝜆 𝑐 20 \lambda_{c}=20 ).

$P_D$	Sensor 1	Sensor 2	C-GM-PHD	Proposed
0.75	23.8612	24.0699	21.6949	21.1739
0.85	21.9649	22.0665	19.4878	19.0956
0.90	20.6366	20.5620	17.2009	16.5622
0.95	19.4173	19.1632	15.5130	14.4411
0.98	18.9021	18.3692	14.3381	13.7560

Table 4. Table 4: OSPA errors with different clutter rates and same detection probability ( p D = 0.95 subscript 𝑝 𝐷 0.95 p_{D}=0.95 ).

${\lambda }_c$	Sensor 1	Sensor 2	C-GM-PHD	Proposed
10	19.1457	18.6978	14.8450	13.3432
20	19.4173	19.1632	15.5130	14.4411
30	19.7298	19.5954	15.8714	15.1374
40	20.2337	19.9615	16.4390	16.0174
50	20.5805	20.8820	17.5043	17.1127

Equations

$${\mathbf{Z}}_k^l$$

$$\begin{array}{rl}& {\pi }_{k|k-1}\left({\mathbf{X}}_k\right)=\int f_{k|k-1}({\mathbf{X}}_k|{\mathbf{X}}_{k-1}){\pi }_{k-1}({\mathbf{X}}_{k-1})\delta {\mathbf{X}}_{k-1}\end{array}$$

$${\pi }_k\left({\mathbf{X}}_k\right)=\frac{g_k\left({\mathbf{Z}}_k|{\mathbf{X}}_k\right){\pi }_{k|k-1}\left({\mathbf{X}}_k\right)}{\int g_k\left({\mathbf{Z}}_k|{\mathbf{X}}_k\right){\pi }_{k|k-1}\left({\mathbf{X}}_k\right)\delta {\mathbf{X}}_k}$$

$$\pi (\mathbf{X})=\exp \left(-{\int }_{\mathbb{X}}v(\mathbf{x})d\mathbf{x}\right)\prod _{\mathbf{x}\in \mathbf{X}}v(\mathbf{x})$$

$$\begin{array}{rl}{v}_{k|k-1}(\mathbf{x})=& \int p_S(\zeta )f_{k|k-1}(\mathbf{x}|\zeta )v_{k-1}(\zeta )d\zeta +{\gamma }_k(\mathbf{x})\end{array}$$

$$\begin{array}{rl}{v}_k(\mathbf{x})=& [1-p_D(\mathbf{x})]v_{k|k-1}(\mathbf{x})\\ & +\sum _{z\in {\mathbf{Z}}_k}\frac{p_D(\mathbf{x})g_k(z|\mathbf{x})v_{k|k-1}(\mathbf{x})}{{\kappa }_k(z)+\int p_D(\zeta )g_k(z|\zeta )v_{k|k-1}(\zeta )d\zeta }\end{array}$$

$$v_k(\mathbf{x})=\sum _{i=1}^N{\alpha }_{k,i}\mathcal{N}(\mathbf{x};{\mathbf{x}}_{k,i},{\mathbf{P}}_{k,i})$$

$${\pi }_k^{1,2}\left(\mathbf{X}\right)=\frac{{\pi }_k^1{\left(\mathbf{X}\right)}^{{\omega }_1}{\pi }_k^2{\left(\mathbf{X}\right)}^{{\omega }_2}}{\int {\pi }_k^1{\left(\mathbf{X}\right)}^{{\omega }_1}{\pi }_k^2{\left(\mathbf{X}\right)}^{{\omega }_2}\delta \mathbf{X}}$$

$${\pi }_k^{1,2}=\arg \underset{\pi }{\inf }\left({\omega }_1{D}_{KL}(\pi \Vert {\pi }_k^1)+{\omega }_2{D}_{KL}(\pi \Vert {\pi }_k^2)\right)$$

$$\frac{D_{KL}(\pi \Vert {\pi }^i)\text{buildrel}\Delta }{=\int \pi (\mathbf{X})\log \frac{\pi (\mathbf{X})}{{\pi }^i(\mathbf{X})}\delta \mathbf{X}.}$$

$$v_k^{1,2}(\mathbf{x})={\left[v_k^1(\mathbf{x})\right]}^{{\omega }_1}{\left[v_k^2(\mathbf{x})\right]}^{{\omega }_2}.$$

$${\left[v_k^l(\mathbf{x})\right]}^{{\omega }_l}\approx \sum _{i=1}^{N^l}{\tilde{\alpha }}_{k,i}^l\mathcal{N}(\mathbf{x};{\tilde{\mathbf{x}}}_{k,i}^l,{\tilde{\mathbf{P}}}_{k,i}^l)$$

$$v_k^{1,2}(\mathbf{x})=\sum _{i=1}^{N^1}\sum _{j=1}^{N^2}{\alpha }_{k,i,j}^{1,2}\mathcal{N}\left(\mathbf{x};{\mathbf{x}}_{k,i,j}^{1,2},{\mathbf{P}}_{k,i,j}^{1,2}\right)$$

$${\mathbf{P}}_{k,i,j}^{1,2}$$

$${\mathbf{m}}_{k,i,j}^{1,2}$$

$${\alpha }_{k,i,j}^{12}$$

$${\tilde{\mathbf{x}}}_{k,i}^l$$

$${\tilde{\mathbf{P}}}_{k,i}^l$$

$${\tilde{\alpha }}_{k,i}^l$$

$$\kappa \left(\omega ,\mathbf{P}\right)\approx \sqrt{\frac{\left|2\pi \frac{\mathbf{P}}{\omega }\right|}{{\left|2\pi \mathbf{P}\right|}^{\omega }}}.$$

$$p_1={\left[400,0\right]}^{\top },p_2={\left[800,0\right]}^{\top }$$

$$v^1(\mathbf{x})$$

$$v^2(\mathbf{x})$$

$${\hat{v}}_p^l(\mathbf{x}){\hat{v}}_q^l(\mathbf{x})\approx 0$$

$${\hat{v}}_p^1(\mathbf{x}){\hat{v}}_q^2(\mathbf{x})\approx 0$$

$${\left[v^l(\mathbf{x})\right]}^{{\omega }_l}={\left[\sum _{p=1}^{M_l}{\hat{v}}_p^l(\mathbf{x})\right]}^{{\omega }_i}\approx \sum _{p=1}^{M_l}{\left[{\hat{v}}_p^l(\mathbf{x})\right]}^{{\omega }_l}$$

$${\left[v^1(\mathbf{x})\right]}^{{\omega }_1}{\left[v^2(\mathbf{x})\right]}^{{\omega }_2}\approx \sum _{p=1}^Q{\left[{\hat{v}}_p^1(\mathbf{x})\right]}^{{\omega }_1}{\left[{\hat{v}}_p^2(\mathbf{x})\right]}^{{\omega }_2}.$$

$${\hat{v}}^{1,2}(\mathbf{x})=\sum _{p=1}^Q{\left[{\hat{v}}_p^1(\mathbf{x})\right]}^{{\omega }_1}{\left[{\hat{v}}_p^2(\mathbf{x})\right]}^{{\omega }_2}.$$

$$\mathbb{X}={\mathbb{X}}_1\cup {\mathbb{X}}_2\cup \cdots \cup {\mathbb{X}}_M$$

$${\mathbb{X}}_p\cap {\mathbb{X}}_q=\varnothing$$

$${\hat{v}}_p^l(\mathbf{x})\approx 0$$

Full text

††journal: Elsevier

1 Introduction

Multitarget tracking (MTT) is a problem of great practical relevance in many different contexts, from traffic monitoring and defense to robotics. In many cases, the information collected from multiple sensors has to combined so as to improve performance. In this context, distributed MTT algorithms have recently gained a considerable attention because of their fault tolerance, flexibility, and reduced computational burden as compared to a centralized fusion framework [1, 2, 3, 4, 5]. In the single-target case, many efficient distributed tracking solutions have been developed over the years both in linear and nonlinear settings. In the former case distributed Kalman filters (KF) are usually adopted [6, 7], while in the latter case solutions are usually based either on nonlinear variant of the KF [8] or on the particle filter (PF) [9]. In the multi-target case, the problem is substantially more challenging due to the many reasons, above all the unknown association between targets and measurements. Traditional solutions include Joint Probabilistic Data Association Filter (JPDAF) [10] and Multiple Hypothesis Tracker (MHT) [11, 12]. However, the JPDAF requires the number of targets to be known a priori while the MHT involves explicit association between targets and measurements, which can involve a large computational load. An alternative solution is the so-called Probability Hypothesis Density (PHD) filter [13, 15], which is based on the theory of Random Finite Sets (RFSs) [13, 14]. The PHD filter can avoid explicit data association and can deal with an unknown and time-varying number of targets. The PHD filter is usually implemented by resorting to either Gaussian Mixtures (GMs) [27] or Sequential Monte Carlo (SMC) approximations [28, 29].

The state of the art approach to distributed multi-sensor MTT is Generalized Covariance Intersection (GCI) [16, 17], also called Exponential Mixture Density (EMD) [18]. The GCI fusion amounts to computing the density that minimizes the sum of the information gains (Kullback-Leibler divergences [19, 7], KLD) from local posteriors, thus avoiding the problem of double-counting of common information [25]. In the past years, any distributed RFS filters based on the GCI fusion rule have been proposed [20, 21, 24, 25, 22, 23].

While the GCI fusion rule performs well in many settings, it has been observed that this fusion rule can be sensitive to a high miss-detection rate [26]. In fact, the GCI fusion rule tends to preserve only tracks that are present in all the local posteriors. This pitfall is exacerbated when the sensors have different fields of views (FoVs). For this reason, recently, some techniques have been proposed for dealing with the problem of different FoVs within the GCI fusion framework. For instance, two possible solutions are proposed in [30] for the SMC-PHD filter. In the first method, particles of different sensors are combined only when they are believed to represent the same target. In the second method, particles corresponding to the same target are hierarchically clustered and used to perform state extraction. However, both methods are prone to estimation errors and can underestimate the target number when nearby targets appear. A distributed fusion algorithm based on the SMC-PHD filter, abandoning the limitation of fully overlapping FoV, is proposed in [31], which classifies the received particles into common particle set and external particle set. As for the GM-PHD filter, [32] proposes a solution to handle different FoVs in the context of simultaneous localization and mapping (SLAM). Specifically, the approach of [32] is based on the idea of using a uniform intensity over the whole region of interest to initialize all the local PHDs, in order to model the total uncertainty on the target (landmark) positions in the unexplored region of the map. A different solution was proposed in [33], which amends the traditional GCI fusion algorithm by considering the distance between Gaussian components (GCs). However, this method cannot solve the false alarm problem and meanwhile the cardinality is overestimated. Similar ideas have recently been exploited to extended the labeled-RFS filter to multi-view sensors [34, 35, 36].

In this paper, we focus on the problem of distributed multi sensor MTT using the GM-PHD filter for a sensor network with different FoVs. An efficient fusion strategy is developed that makes use of a clustering algorithm (CA) to improve robustness and performance of GCI fusion rule.

The main contributions of this paper are summarized as follows.

We analyse the GCI fusion mismatch phenomenon caused by the limited sensor FoVs both from the theoretical point of view and via a simulation experiment. 2. 2.

Motivated by the aforementioned analysis, we propose an effective and robust solution for GCI fusion of PHDs generated from sensors with different FoVs that is immune to the effect of fusion mismatch. The solution includes two parts: a parallelized GCI fusion in the commong FoV, and a compensation strategy outside the common FoV. 3. 3.

We provide an error analysis in terms of $L_{1}$-norm between the parallelized GCI fusion and the traditional GCI fusion algorithm. Specifically, we show that the error tends to zero when the partial intensities generated by the clustering algorithm are well separated. 4. 4.

We implement the proposed fusion algorithm based on the GM-PHD filter with adaptive birth model. In the implementation, the union find set (UFS) [37, 38] algorithm is adopted to ensure that the GC subsets generated during clastering are disjoint, and Optimal SubPattern Assignment (OSPA) metric [39] is used to compute the similarity between different subsets. In addition, Murty’s algorithm [41] is utilized to obtain the association pairs effectively.

The remainder of the paper is organized as follows. Section II introduces the main background notions including the multi-target Bayes filter, the PHD filter, and the GCI fusion rule. The mismatch of the GCI fusion for multi-view sensors is analyzed in Section III. The proposed fusion algorithm is described in Section IV, including a detailed description of the proposed GM implementation. Section V is dedicated to the simulations, while conclusions and perspectives for future work are given in Section VI.

2 Background

This section provides an overview of the background material relevant to our work. Some of the related notations are summarized in Table I.

2.1 Multi-target Bayes filter

In this paper, the multi-target state is modelled as a RFS, i.e. a variable which is random in both the number of elements (i.e. the number of targets) and the values of the elements (i.e. the target states). From a probabilistic point of view, the RFS ${\mathbf{X}}_{k}$ is completely characterized by the multi-object posterior ${\pi_{k}}\left(\mathbf{X}_{k}\right)$, corresponding to the multi-object density of ${\mathbf{X}}_{k}$ conditional to all the observations collected up to and including time $k$.

Given the multi-target posterior at time $k-1$ and the observation set $\mathbf{Z}_{k}$ at time $k$, the multi-target Bayes filter propagates the multi-target posterior in time via the Chapman-Kolmogorov prediction

[TABLE]

and then updates the density using the Bayes formula

[TABLE]

where $f_{k|k-1}(\mathbf{X}_{k}|\mathbf{X}_{k-1})$ and ${g_{k}}({{\mathbf{Z}_{k}}|{\mathbf{X}_{k}}})$ are the multi-target transition density and the multi-target likelihood, respectively. The integrals in (1)-(2) are to be intended in the set-integral sense [14].

2.2 The PHD filter

The multi-target Bayes filter described in the previous section is in general intractable from the computational point of view. A possible approximation is the so-called PHD filter [13, 15] that is based on the idea of representing the multi-target set as a Poisson RFS. The multi-object density $\pi({\mathbf{X}})$ of a Poisson RFS ${\mathbf{X}}$ takes the form

[TABLE]

where ${v}(\mathbf{x})$ is the so-called PHD or intensity function. The PHD ${v}(\mathbf{x})$ can be interpreted as the local density of the number of targets in the single-target state space. In fact, given a region $\mathcal{X}\subseteq\mathbb{X}$, the expected number of targets inside such a region can be computed as $\int_{\mathcal{X}}v(\mathbf{x})\,d\mathbf{x}\,,$ and the total number of expected targets in the whole state space is given immediately by $\int_{\mathbb{X}}v(\mathbf{x})\,d\mathbf{x}\,.$ The PHD filter propagates in time the PHD ${v_{k}}(\mathbf{x})$ of the multi-target state conditioned to all the collected observations. The detailed PHD filter recursions are given as follows.

[TABLE]

[TABLE]

where $\gamma_{k}(\mathbf{x})$ is the prior PHD of target births at time $k$, and $\kappa_{k}(z)$ is the clutter density.

From (4) and (5), we can see that the PHD filter operates in the single-target state space and avoids data associations. In practice, instead of propagating the multi-target posterior density, the PHD filter chooses to propagate the first-order moment of the true posterior by approximating the multi-target set as a Poisson RFS. The most typical implementations of the PHD filter are based either on GMs or SMC. In this paper, we mainly focus on the GM implementation. Accordingly, the posterior PHD ${v_{k}}(\mathbf{x})$ is represented as

[TABLE]

where ${{\cal N}}(\mathbf{x};\mathbf{x}_{k,i},\mathbf{P}_{k,i})$ denotes the Gaussian PDF with mean $\mathbf{m}_{k,i}$ and covariance $\mathbf{P}_{k,i}$.

Concerning the birth model it can be either predefined [27, 28] using prior information or adaptive [42, 43, 44, 45, 46]. The adaptive birth model (ABM) is based on the idea of using the measurements to initialize new components in the PHD, and is preferable when new targets can born in the whole state space. As discussed in [42, 43, 44, 45, 46], the ABM makes use of a uniform (or partially uniform) birth intensity so as to reflect the uncertainty on the location of newborn targets.

2.3 GCI fusion rule

Consider two multi-target posteriors $\pi^{1}_{k}\left({\mathbf{X}}\right)$ and $\pi^{2}_{k}\left({\mathbf{X}}\right)$ conditioned on measurement set sequences from two different sensors. When the correlation between the two measurement set sequences is unknown, the two multi-target posteriors can be fused by resorting to the so-called GCI fusion rule [16]. Under the GCI fusion rule, the two multi-target posteriors are fused into the multi-target posterior

[TABLE]

where the positive scalars $\omega_{1},\omega_{2}$, with $\omega_{1}+\omega_{2}=1$, are weights determining the relative importance of each multi-target posterior. The fusion rule (7) was originally proposed for the fusion of PDFs in the statistics literature [47], and then generalizes to the case of RFS densities by Mahler [16]. As discussed in [48], the fused density given in (7) is the one minimizing the weighted sum of the KLD with respect to the densities to be fused

[TABLE]

where

[TABLE]

See [20] for a proof in the case of RFS densities. While, for the sake of simplicity, only the case of two densities is considered, all the above considerations apply also in the case of more than two sensors.

Let now the multi-target posteriors to be fused be Poisson RFSs with PHDs $v^{1}_{k}(\mathbf{x})$ and $v^{2}_{k}(\mathbf{x})$, respectively, as it happens when each sensor runs locally a PHD filter. Then, application of the GCI fusion rule yields again a Poisson RFS with fused PHD given by [17]

[TABLE]

While the result of the above cannot be computed in closed-form when the PHD is represented by a GM, because the exponentiation of a GM is in general not a GM, various approximation strategies exist. For instance, one possible way of overcoming this drawback amounts to approximating each power $\left[{v^{l}_{k}}(\mathbf{x})\right]^{\omega_{l}}$ as a Gaussian mixture

[TABLE]

and, then, computing the fused PHD as

[TABLE]

where

[TABLE]

As discussed in [20], when the Gaussian Components (GCs) of each local PHD are well separated, the approximation in (11) can be simply obtained as

[TABLE]

where

[TABLE]

More accurate approximations, like the one of [49] based on sigma-points, must instead be used in the case of closely-spaced Gaussian components.

It can be seen from (12) that, there will be $N_{1}N_{2}$ GCs after fusion. In the general case of $S$ sensors, centralized GCI fusion [50] would yield $\prod\nolimits_{s}{{N_{s}}}$ GCs, thus requiring a high computational load. Hence, when $S>2$ a sequential fusion strategy [51] may be preferable in practice.

3 Analysis of the mismatch of GCI fusion for sensors with different FoVs

In this section, we provide a brief overview of the behavior of GCI fusion when applied to sensors with different FoVs, and we discuss the reason why GCI fusion may fail in this case. To this end, we make use of the following example.

Example 1: Consider a distributed sensor network with two sensors, each one running a GM-PHD filter with adaptive birth model (AB-GMPHD). The surveillance region is $[0,1500]m\times[0,1000]m$. The dynamics and observation models are the same as in [52]. The survival probability is $p_{S,k}=0.99$. The positions of the two sensors are

[TABLE]

Each sensor has a limited FoV and, specifically, can detect only targets with relative angle in the interval $[-60^{\circ},60^{\circ}]$. The detection probability is $p_{D}=0.98$ within the FoV, and zero otherwise. The true tracks are shown in Fig. 1 (a) and the simulation results including OSPA errors [39] and cardinality estimates are shown in Fig. 1 (b). The OSPA parameters are set equal to $c=30m$, $p=2$.

It can be seen from Fig. 1 (a) that target 1 is detected only by sensor 1, while target 2 is detected only by sensor 2. Accordingly, the PHD of sensor 1 will not contain GCs representing target 2, while the PHD of sensor 2 will not contain GCs representing target 1. Hence, the GCs of sensor 1 will be far from the GCs of sensor 2. The result is that tracking performance after fusion is substantially worse than that of a single sensor.

The reason for the behavior highlighted in Example 1 is that the essence of GCI fusion rule is the weighted multiplication between target densities, and the fusion process will work well only when both sensors detect the same targets. In fact, when a sensor does not detect a target, its PHD will be close to zero in the corresponding region of the state space. Hence, even if the other sensor is able to detect the target and, hence, has GCs with no-negligible weights in that region, application of the GCI fusion rule will yield a substantial decrease of such weights. Consequently, the target may be lost in the fused multi-target distribution, as it happens in Example 1. This state of affairs can be also understood by looking at the GM-based implementation (12) of GCI fusion. In fact, we can see from (15) that a large distance between GCs results in a small fused weight, since ${\cal N}({\widetilde{\mathbf{x}}}_{k,i}^{1}-{\widetilde{\mathbf{x}}}_{k,j}^{2};0,{{{\widetilde{\mathbf{P}}}_{k,i}^{1}}}+{{{\widetilde{\mathbf{P}}}_{k,j}^{2}}})$ tends to [math] when the distance ${\widetilde{\mathbf{x}}}_{k,i}^{1}-{\widetilde{\mathbf{x}}}_{k,j}^{2}$ increases. .

4 The CA-GCI fusion algorithm

In order to overcome the drawbacks of the traditional GCI fusion algorithm applied in multi-sensor fusion with partial overlapping FoVs, we propose a novel clustering-based GCI fusion algorithm, named CA-GCI fusion algorithm. First, we describe the basic ideas of the proposed fusion algorithm and analyse the approximation error as compared to traditional GCI fusion. Then, we discuss how the fusion rule can be modified so as to deal with the different FoVs. Finally, detailed guidelines on the implementation are provided.

For ease of presentation, the problem of fusing the PHDs coming from two sensors is considered. All the presented techniques can be readily generalized to sensor networks with multiple sensors ($S>2$) by resorting to a sequential fusion strategy.

4.1 Fusion for overlapping FoVs

Let the two PHDs $v^{1}(\mathbf{x})$ and $v^{2}(\mathbf{x})$ to be fused be represented in GM form as in (6). Suppose now that the GCs of each of the two PHDs are grouped together in clusters so that the two PHDs can be decomposed as

[TABLE]

where each divided intensity function (DIF) ${\hat{v}_{p}^{i}(\mathbf{x})}$ consists of the linear combination of GCs. The clusters are constructed so that:

For each sensor $l$, clusters are well separated in the sense that for any $\mathbf{x}\in\mathbb{X}$ at most one of the DIFs ${\hat{v}_{p}^{l}(\mathbf{x})}$ is substantially different from zero, i.e.,

[TABLE]

for $l=1,2\,$ and for any $p\neq q$; 2. 2.

The first $Q$ DIFs of sensors 1 are matched with the first $Q$ DIFs of sensor 2. These DIFs correspond to the regions of the state space in which both sensors have detected targets. Without loss of generality, we suppose that ${\hat{v}_{p}^{1}(\mathbf{x})}$ is matched with ${\hat{v}_{p}^{2}(\mathbf{x})}$, for $p=1,\ldots,Q$, so that

[TABLE]

for any $p\neq q$; 3. 3.

The last $M_{1}-Q$ DIFs of sensor 1 are not matched with any DIF of sensor 2, i.e.,

[TABLE]

for any $p=Q+1,\ldots,M_{1}$ and for any $q$. These DIFs correspond to targets detected only by sensor 1; 4. 4.

The last $M_{2}-Q$ DIFs of sensor 2 are not matched with any DIF of sensor 1, i.e.,

[TABLE]

for any $p$ and for any $q=Q+1,\ldots,M_{2}$. These DIFs correspond to targets detected only by sensor 2.

An efficient algorithm for the construction of the matched clusters based on the similarity between GCs will be provided later on in Section IV.D.

Given the matched clusters and the corresponding DIFs for the two PHDs, an approximate expression for the GCI fusion rule can be easily derived. More specifically, in view of property (a) the DIFs are well separated and, hence, the exponentiations of the two PHDs can be approximated as

[TABLE]

for $l=1,2\,$. Further, by exploiting properties (b)-(d), we also have that

[TABLE]

This means that the matched DIFs can be fused in parallel and the fused PHD $v^{1,2}(\mathbf{x})$ can be approximated by the intensity

[TABLE]

From the above equation, it is clear that only the DIFs corresponding to targets detected by both sensors are preserved after fusion, while the DIFs corresponding to targets detected by either sensor 1 or 2 (see points (c) and (d) above) disappear. Hence, at a basic level, the GCI fusion rule can be interpreted as performing an intersection of the DIFs (and hence of the tracks) of the two sensors. This is consistent with the behavior observed empirically in Section III. Hereafter, the approximated fusion rule (28) will be referred to as parallelized GCI.

Remark 1**.**

We point out that, in the worst case where all targets are so close that it is impossible to partition the PHDs into separated DIFs, there will be only one DIF for each sensor and the parallelized GCI fusion (28) will reduce to the standard GCI fusion algorithm.

4.2 Analysis of the approximation error

In this section, we analyze the error committed when the approximated fusion rule (28) is used. To this end, it is convenient to partition the state space $\mathbb{X}$ into $M=M_{1}+M_{2}-Q$ subsets $\mathbb{X}_{1}$, $\ldots$, $\mathbb{X}_{M}$ so that

[TABLE]

and

[TABLE]

for any $p\neq q$. The partition is built so that:

for $p=1,\ldots,Q$, the set $\mathbb{X}_{p}$ contains the tracks represented by the matched DIFs ${\hat{v}_{p}^{1}(\mathbf{x})}$ and ${\hat{v}_{p}^{2}(\mathbf{x})}$ in the sense that

[TABLE]

for $l=1,2$ and for any $\mathbf{x}\in\mathbb{X}\setminus\mathbb{X}_{p}$; 2. 2.

for $p=Q+1,\ldots,M_{1}$, the set $\mathbb{X}_{p}$ contains the tracks represented by the unmatched DIF ${\hat{v}_{p}^{1}(\mathbf{x})}$ of sensor 1 in the sense that

[TABLE]

for $\mathbf{x}\in\mathbb{X}\setminus\mathbb{X}_{p}$; 3. 3.

for $p=M_{1}+1,\ldots,M$, the set $\mathbb{X}_{p}$ contains the tracks represented by the unmatched DIF ${\hat{v}_{p-M_{1}+Q}^{2}(\mathbf{x})}$ of sensor 2 in the sense that

[TABLE]

for $\mathbf{x}\in\mathbb{X}\setminus\mathbb{X}_{p}$.

Notice that the above construction is possible in view of fact that the DIFs are supposed to be well separated and matched as discussed in points (a)-(d) above. Notice also that this partition need not be actually determined, since it is just instrumental for the analysis, but it is not used in the actual fusion rule.

Let us now decompose the two PHDs as

[TABLE]

for $l=1,2$ , where each partial intensity $v_{p}^{l}(\mathbf{x})$ is defined as

[TABLE]

with ${\bf{1}}_{\mathbb{X}_{p}}(\cdot)$ the indicator function

[TABLE]

Notice that the partial intensities $v_{p}^{l}(\mathbf{x})$ satisfy properties (a)-(b) exactly in that

[TABLE]

for any $l,r$ and any $p\neq q$. As a consequence, the fused PHD $v^{1,2}(\mathbf{x})$ can be equivalently written as

[TABLE]

Notice also that, since the DIFs are determined so as to approximately satisfy properties (a)-(d), we have that

[TABLE]

The above relationships can be quantified by bounding the discrepancies in terms of a scalar $\delta$

[TABLE]

where, given a function $f:\mathbb{X}\rightarrow\mathbb{R}$, $\|f\|_{\infty}$ denotes its infinite norm $\|f\|_{\infty}=\sup_{\mathbf{x}\in\mathbb{X}}|f(\mathbf{x})|$. Clearly, the more well separated are the matched DIFs ${\hat{v}_{p}^{l}(\mathbf{x})}$ the smaller is the scalar $\delta$. In other words, $\delta$ can be interpreted as a quantitative measure of how much each DIF ${\hat{v}_{p}^{l}(\mathbf{x})}$ is concentrated in the corresponding domain of the partition. For the readers’ convenience, Fig. 2 shows the difference between the DIFs (in blue) and the partial intensities (in red) in a simple example.

The following result provides a quantification of the error committed when the parallelized GCI fusion (28) is used in place of the exact GCI fusion (10).

Proposition 1**.**

Let the two PHDs ${\left[{{v^{1}}(\mathbf{x})}\right]^{{\omega_{1}}}}$ and ${\left[{{v^{2}}(\mathbf{x})}\right]^{{\omega_{2}}}}$ to be fused be bounded, and let

[TABLE]

*for any $p$, $q$, where $\|\cdot\|_{1}$ denotes the $L_{1}$-norm, ${\left\|f\right\|_{1}}=\int{\left|{f(\mathbf{x})}\right|}d\mathbf{x}$.

Then, the following bound holds*

[TABLE]

The proof of Proposition 1 is given in the Appendix A. It follows from Proposition 1 that the discrepancy ${\left\|{\hat{v}^{1,2}-v^{1,2}}\right\|_{{1}}}$ between the parallelized GCI fusion (28) and the exact GCI fusion (10) is small when the scalar $\delta$ is small (i.e., when the matched DIFs are well separated).

4.3 Fusion for different FoVs

When the sensors share the same FoV, the local PHDs contain information on the same targets and the GCI fusion performs well. However, as discussed in Section III, when different sensors have different FoVs, the GCI fusion of local posteriors can lead to unsatisfactory results. This is because, at a basic level, GCI fusion can be interpreted as performing an intersection among the tracks of different sensors. In fact, recalling that GCI fusion of two PHDs can be approximated as (28)), there will be $M_{1}-Q$ DIFs for sensor 1, namely $\hat{v}^{1}_{Q+1}(\mathbf{x}),\ldots,\hat{v}^{1}_{M_{1}}(\mathbf{x})$, and $M_{2}-Q$ DIFs for sensor 2, namely $\hat{v}^{2}_{Q+1}(\mathbf{x}),\ldots,\hat{v}^{2}_{M_{2}}(\mathbf{x})$, that do not participate in the fusion process so that the corresponding tracks disappear after fusion. Hereafter, we will refer to these DIFs as unconfirmed DIFs. Then, the idea is to modify the fusion rule so as to preserve, at least partially, the information contained in these unconfirmed DIFs, in order account for the fact that each sensor has only a partial view of the overall surveillance area. Of course, special care has to be taken in this operation since the unconfirmed DIFs, corresponding to tracks present in only one of the sensors, may actually be false targets. With this respect, it is important to distinguish the unconfirmed DIFs contained in the intersection of the FoVs from those contained only in the FoV of one sensor. In fact, while it is reasonable to think that the former ones correspond to false alarms, no similar conclusion can be drawn on the latter ones. To this end, it is convenient to introduce the following concept.

Definition 1**.**

Given an unconfirmed DIF $\hat{v}^{l}_{p}(\mathbf{x})$, we say that $\hat{v}^{l}_{p}(\mathbf{x})$ is observed by sensor $r$ if

[TABLE]

where $\gamma\in(0,1)$ is a threshold.

Let us now introduce, for each DIF of sensor $l$, a binary variable $\beta^{lr}_{p}$ taking value $1$ if $\hat{v}^{l}_{p}(\mathbf{x})$ is not observed by sensor $r$, and value [math] otherwise. Then, the idea is to preserve all the DIFs of one sensors that are not observed by the other. Accordingly, we define the preserved PHDs of sensors $1$ and $2$ as

[TABLE]

We have now to distinguish two cases: namely complete trust, and partial trust. The complete trust case refers to the situation in which the sensor FoVs are time-invariant, and the sensors do not receive any feedback after fusion so that the local PHDs only contain information on the targets detected within the local constant FoV. In this case, it is clear that it makes sense to perform fusion only in the intersection $Fo{V_{1}}\cap Fo{V_{2}}$, whereas the fusion center has to completely trust sensor $1$ in its exclusive field of view $Fo{V_{1}}\setminus Fo{V_{2}}$ (and, similarly, a complete trust should be given to sensor $2$ in $Fo{V_{2}}\setminus Fo{V_{1}}$). Hence, in this case, the total fused intensity can be taken equal to

[TABLE]

On the other hand, in many situations, local PHDs also contain information outside the local FoV. For instance, this can happen when the FoV is time-varying, like in the case of mobile sensors, or when there is a feedback from the fusion center to the local nodes111Even in a centralized setting, feedback is important. In fact, thanks to the information fed back from the fusion center, each sensor node is able to keep track of all the targets currently present in the total surveillance area. For instance, this means that when a target already detected by one sensor first enters the FoV of another sensor, no initialization is required because the predicted density already contains a GC corresponding to that target. , or when the fusion is distributed/decentralized. Clearly, in these cases, it would be inappropriate to completely trust a DIF that is present in one sensor but is not present in the other. In this case, we follow a different approach by introducing a compensation intensity in the PHDs to be fused (a similar idea was proposed in [36] in the context of GCI fusion of labelled RFS densities). This amounts to adding a new term to each local PHD that models the birth of new targets outside the current FoV. The new term, which corresponds to a uniform (uninformative) intensity outside the current FoV, makes it possible for each sensor to initialize new tracks from those of the other sensor (similarly to what happens in adaptive birth model wherein a uniform birth density within the FoV is used to initialize new tracks from the measurements). The value $\Delta>0$ of this additional intensity plays the role of a confidence factor and can be also interpreted as the expected target intensity. Then, the fused intensity can be written as

[TABLE]

where the two additional terms arise from the fusion of the preserved PHD ${{\bar{v}}^{i}}(\mathbf{x})$ of one sensor with the uniform (uninformative) intensity $\Delta$ of the other, and $\bar{\omega}_{1}\in\left[0,1\right]$ and $\bar{\omega}_{2}\in\left[0,1\right]$ are fusion weights for non-overlapping areas, which can be different from $\omega_{1}$ and $\omega_{2}$. Notice that $\Delta$, $\bar{\omega}_{1}$, $\bar{\omega}_{2}$ can be tuned according to the actual sensor network under consideration. For instance, by setting $\Delta=1$, $\bar{\omega}_{1}=1$ and $\bar{\omega}_{2}=1$ we retrieve the case of complete trust (46), whereas by choosing $\Delta=0$ we consider only the information pertaining to the intersection of the FoVs.

4.4 The implementation of the proposed algorithm

In this section, the detailed GM implementation of the proposed CA-GMPHD-GCI algorithm is provided, including the construction of DIFs. The GCI fusion between matched DIFs is efficiently implemented in a parallelized way, while the preserved DIFs are fused using the compensation intensity introduced in the previous section.

Given the PHDs of the two sensors in GM form (6), the following eight implementation steps are carried out.

1. Grouping by pre-clustering: Considering that the number of clusters is unknown, we choose the GCs whose weight $\alpha$ satisfies

[TABLE]

as candidate cluster centers, where $T_{\alpha}$ is a pre-specified threshold. Specifically, $T_{\alpha}$ is suitably tuned so as to ensure that newborn targets can be quickly detected while avoiding that the number of groups becomes too large. The corrected Mahalanobis distance

[TABLE]

is used to measure the distance between two GCs having mean/covariance $({\mathbf{x}}_{1},{\mathbf{P}}_{1})$ and $({\mathbf{x}}_{2},{\mathbf{P}}_{2})$, respectively. For each sensor, clusters are formed by grouping together all GCs whose distance from a cluster center is less than a given threshold $T_{d}$. The same distance criterion is used for the two sensors. After the clustering process is finished, the GCs of the two sensors are decomposed into the following subsets

[TABLE]

where $a$ and $b$ are the numbers of subsets for sensor $1$ and $2$, respectively. Due to the different FoVs, $a$ and $b$ are in general different.

2. Making the clusters disjoint: The subsets resulting from the above-described procedure need not be disjoint, i.e. each GC is in general can be assigned to more than one subset $L_{l}^{i}$. This usually happens when there are targets at close distance. In order to ensure that each GC is associated with a single subset, we resort to the union find set (UFS) algorithm [37, 38]. After application of UFS, a new decomposition is obtained for the GCs of the two sensors

[TABLE]

where $M_{1}\leq a$ and $M_{2}\leq b$ are the new numbers of subsets for sensor $1$ and $2$, respectively, and where, by construction, $L^{l}_{p}\cap L^{l}_{q}=\emptyset$ for any $p\neq q$. Hence, denoting by $N^{l}_{p}$ the number of GCs in the subset $L^{l}_{p}$ and by $N^{l}$ the total number of GCs of sensor $l$, we have ${N^{l}}=\sum\nolimits_{p=1}^{{M_{l}}}{N_{p}^{l}}$. Notice that each of the subset $L^{l}_{p}$ defines a DIF $\hat{v}^{l}_{p}(\mathbf{x})$ obtained as combination of the GCs in $L^{l}_{p}$.

3. Matching subsets: In order to find a matching between the subsets of sensor $1$ and those of sensor $2$, i.e. between the DIFs of sensor $1$ and those of sensor $2$, we use the OSPA distance [39] in order to measure the distance (dissimilarity) between subsets. Then, the following ${M_{1}\times M_{2}}$ distance matrix is constructed

[TABLE]

where $d_{p,q}$ denotes the OSPA distance between $L_{p}^{1}$ and $L_{q}^{2}$.

Then an optimal matching between the subsets $L^{1}$ and $L^{2}$ is found by solving a linear assignment problem [40]. Specifically, Let us denote by $S$ the binary assignment matrix defined so that $S_{p,q}=1$ if $L_{p}^{1}$ is associated to $L_{q}^{2}$, and $S_{p,q}=0$ otherwise. Further let us suppose, without loss of generality, that $M_{1}\geq M_{2}$. Then, we use Murty’s algorithm [41] to find the matrix $S^{*}$ that minimizes

[TABLE]

under the constraint $\sum_{p=1}^{M_{1}}S_{p,q}=1$.

Afterwards, for every associated pair $(p,q)$ such that $S_{p,q}=1$, we check whether the distance $d_{p,q}$ is below or above a given threshold $T_{r}$. In the former case, the association is considered valid. Conversely, if the distance is too large, the association is discarded. As a result, in general only $Q\leq\min\{M_{1},M_{2}\}$ pairs $(p,q)$ are considered valid associations (these are the DIFs corresponding to the targets detected by both sensors), while $M_{1}-Q$ subsets of sensor $1$ and $M_{2}-Q$ subsets of sensor 2 are not associated (these are the unconfirmed DIFs corresponding to the targets detected only by one of the sensors). Finally, the subsets of the two sensors are reordered so that the first $Q$ pairs $(L_{p}^{1},L_{p}^{2})$, $p=1,\ldots,Q$, correspond to matched subsets/DIFs.

5. Selecting the preserved DIFs: After the match procedure is finished, we resort to (43) to determine which of the unconfirmed DIFs should be preserved and which, instead, should be deleted. Specifically, each unmatched subset $L^{1}_{p}$, $p=Q+1,\ldots,M_{1}$ of sensor is preserved (i.e. $\beta^{12}_{p}=1$) if and only if

[TABLE]

If, instead, the above condition is not satisfied, the GCs belonging to the subset $L_{p}^{1}$ are deleted. An analogous procedure is followed for the unmatched subsets $L^{2}_{p}$, $p=Q+1,\ldots,M_{2}$ of sensor 2.

6. Fusion of the matched DIFs: For each association pair ${\left({L_{p}^{1},L_{p}^{2}}\right)}$, $p=1,\ldots,Q$, GM approximations of the powers $[\hat{v}^{1}_{p}(\mathbf{x})]^{\omega_{1}}$ and $[\hat{v}^{2}_{p}(\mathbf{x})]^{\omega_{1}}$ of the corresponding DIFs as described in Section II-C. Then, given the GMs

[TABLE]

for $p=1,\ldots,Q$ and $l=1,2$, the parallelized GCI fusion structure is adopted. Accordingly, the fused PHD for the matched DIFs is obtained as

[TABLE]

where

[TABLE]

Remark 2**.**

In terms of computational burden, the traditional GCI fusion would require to compute and store $\left(\sum\nolimits_{p=1}^{Q}{{N_{p}^{1}}}\right)\cdot\left(\sum\nolimits_{p=1}^{Q}{{N_{p}^{2}}}\right)$ GCs. On the other hand, the parallelized GCI fusion only requires the computation of $\sum\nolimits_{p=1}^{Q}\left({{N_{p}^{1}}\cdot{N_{p}^{2}}}\right)$ GCs, thus significantly reducing the computational burden when $Q>1$.

7. Fusion of the preserved DIFs: For each unmatched subset $L^{l}_{p}$ corresponding to a preserved DIF (as determined in step 5), a GC approximation of the power $[\hat{v}^{l}_{p}(\mathbf{x})]^{\bar{\omega}_{l}}$ is computed as in (57). Then, the total fused density takes the form

[TABLE]

where clearly the terms corresponding to the subsets deleted at step 5 need not be computed because $\beta_{p}^{lr}=0$.

8. Target estimate extraction: Following [27], after the fused density has been obtained, the extraction of multi-target state estimates can be executed in a straightforward way by choosing the GCs whose weight $\alpha$ is greater than a given threshold, e.g., 0.5.

5 Simulation results

In this section, we test the tracking performance of the proposed fusion algorithm (CA-GCI) using the GM-PHD filter with adaptive birth model. The proposed algorithm is compared to the C-GM-PHD filter of [33] using the OSPA error (i.e. the OSPA distance [39] between the true and the estimated target sets) as performance index. The OSPA distance parameters are set to be $c=30m$ and $p=2$. The errors are averaged over 200 independent Monte Carlo runs.

Tracking model and scenario: Consider the problem of tracking an unknown and time-varying number of targets observed in clutter. The single-target state is

[TABLE]

where $(p_{x,k},p_{y,k})$ and $({\dot{p}_{x,k}},{\dot{p}_{y,k}})$ represent the target position and, respectively, velocity in Cartesian coordinates at time $k$. Each target moves according to the following dynamic model

[TABLE]

[TABLE]

where $I_{n}$ and $0_{n}$ denote, respectively, the $n\times n$ identity and zero matrices. $T_{s}=1s$ is the sampling period, $\mathbf{Q}$ is the process noise covariance and $\sigma_{w}=2m/s^{2}$ is the standard deviation of the process noise.

The considered scenario is depicted in Fig. 3. Overall, $11$ targets enter the scenario at different time instants as detailed in Table II. Two sensors located at $(400m,0m)$ and $(800m,0m)$, respectively, provide measurements of the unknown targets. Each sensor has a limited FoV as depicted in Fig. 3. Within the FoV the probability of detection is constant, i.e.

[TABLE]

for $l=1,2$. When a target is detected, a measurement is generated according to the linear noisy model

[TABLE]

where $\mathbf{R}$ is the measurement noise covariance and $\sigma_{\varepsilon}=10m$ is the standard deviation of the measurement noise. Clutter follows a uniformly distributed Poisson RFS with an average of 20 clutter points per scan ($\lambda_{c}=20$).

For the PHD filters, the survival probability $p_{s}$ is constant and equal to $0.99$. Following [43], the covariance of newborn targets is set to a large value so as to ensure they can cover the surveillance area as much as possible. For the CA-GCI fusion algorithm, the pre-clustering threshold and the weight threshold are set equal to $T_{d}=15$ and $T_{\alpha}=0.02$, respectively. The matching threshold is $T_{r}=15$. The fusion weights for the preserved DIFs are $\bar{\omega}_{1}=\bar{\omega}_{2}=0.8$ and the confidence factor is $\Delta=0.9$. The weights for the fusion of the matched DIFs are $\omega_{1}=\omega_{2}=0.5$, which on average achieves nearly optimal results [53].

Simulation results: To assess the effectiveness of the proposed algorithm, we examine the fusion performance in a situation. Comparisons in terms of OSPA errors and cardinality estimates for the considered filters are shown in Figs. 4 and 5, respectively. From Fig. 4, we can see that both the proposed method and the C-GM-PHD filter are able to improve the local sensor performance. Furthermore, the proposed algorithm significantly outperforms the C-GM-PHD filter, especially in terms of cardinality estimate (Fig. 5). In fact, the C-GM-PHD filter substantially overestimates the true cardinality during the time interval $[30,70]s$, because the weights of the fused GCs are computed using the weights of multiple GCs. In contrast, the proposed algorithm is able to closely follow the true cardinality, a part from a small delay in the detection of newborn targets.

In order to verify the robustness of the proposed algorithm, Monte Carlo simulations have been performed with different detection probabilities and clutter rates. A comparison for different $p_{D}$ and same $\lambda_{c}$ is given in Table III. As expected for all the considered filters, the OSPA error increases as the $p_{D}$ increases, but the proposed algorithm provides the best performance. Further, in Table IV a comparison for different $\lambda_{c}$ and same $p_{D}$ is given. Also in this case, the proposed algorithm turns out to be more robust with respect to an increase in the clutter rate.

6 Conclusions and future work

In this paper, an effective and robust fusion algorithm has been proposed in a distributed setting by combining the GCI fusion rule with a suitable CA. Starting from an analysis of the pitfalls of the GCI fusion rule in the case of sensors with different FoVs, a solution to this problem based on the decomposition of the local PHDs into well-separated components have been presented. Then, two different fusion strategies have been presented: a parallelized GCI fusion strategy for combining the components in the common FoV and detected by all sensors; and a suitable compensation strategy to preserve the components outside the common FoV and, hence, not present in all the local PHDs. The effectiveness of the proposed approach has been analyzed by means of simulations in a challenging tracking scenario. An extension of the method to the GM-CPHD filter [20] considering both intensity function and cardinality distribution, is a relevant topic for future research.

Acknowledgment

This work was supported in part by the Chang Jiang Scholars Program, in part by the 111 Project No. B17008, in part by the National Natural Science Foundation of China under Grant 61771110, in part by the Fundamental Research Funds of Central Universities under Grant ZYGX2016J031, and in part by the Chinese Postdoctoral Science Foundation under Grant 2014M550465 and Special Grant 2016T90845.

7 Appendix

Proof.

Combining (28) with (38), we can derive the bound on ${\left\|{\hat{v}^{1,2}-v^{1,2}}\right\|_{{1}}}$. For convenience, the following shorthand notation will be used

[TABLE]

By applying the triangular inequality, the $L_{1}$-norm ${\left\|{{{\hat{v}}^{1,2}}-{v^{1,2}}}\right\|_{1}}$ can be decomposed into two parts as shown in (62). Next, two parts are discussed separately.

Part 1 in (62):

The following upper bound can be obtained

[TABLE]

where

(67) holds because

[TABLE] 2. 2.

(67) is obtained by applying the triangular inequality; 3. 3.

(67) holds because $\|f\,g\|_{1}\leq\|f\|_{\infty}\,\|g\|_{1}$; 4. 4.

(67) can be ontained in a similar way as (67); 5. 5.

As to (68), in view of (40) and (41), we have

[TABLE]

Then, since $\omega_{1}\leq 1$ and $\omega_{2}\leq 1$, we have

[TABLE]

As a result, (68) is obtained.

Part 2 in (62): Similarly to part 1, the following upper bound can be obtained

[TABLE]

Combining (68) and (78) and recalling that $M=M_{1}+M_{2}-Q$, the upper bound

[TABLE]

can be readily obtained. Hence, Proposition 1 holds. ∎