{'fraction_non_alphanumeric': 0.06252942099482191, 'fraction_numerical': 0.016671897065746115, 'mean_word_length': 3.89684978870534, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Scalable Track-Before-Detect Method With Poisson/Multi-Bernoulli Model

We propose a scalable track-before-detect (TBD) tracking method based on a Poisson/multi-Bernoulli model. To limit computational complexity, we approximate the exact multi-Bernoulli mixture posterior probability density function (pdf) by a multi-Bernoulli pdf. Data association based on the sum-product algorithm and recycling of Bernoulli components enable the detection and tracking of low-observable objects with limited computational resources. Our simulation results demonstrate a significantly improved tracking performance compared to a state-of-the-art TBD method.


I. INTRODUCTION
Multiobject tracking aims to estimate the time-dependent number and states of multiple objects based on data provided by one or more sensors. The measurements of conventional multiobject tracking methods are produced by a detector that performs thresholding and, possibly, some further preparatory processing of the raw sensor data [1]- [9]. This preparatory processing step reduces data flow and computational complexity, but results in a loss of tracking information. Thus, the conventional detect-than-track approach can lead to a significantly reduced tracking performance, especially in scenarios with low signal-to-noise (SNR) ratio [10]- [17].
In such scenarios, tracking methods that use raw sensor data as measurements, can potentially achieve an improved tracking performance. Many existing track-before-detect (TBD) methods are based on batch processing using maximum likelihood estimation [10], the Hough transform [11], and dynamic programming techniques [12]. However, due to their high computational complexity, they are often unsuitable for realtime operation. A well-established real-time TBD method is the Histogram Probabilistic Multi-Hypothesis Tracker (H-PMHT) [13], which is based on the expectation-maximization algorithm. However, tuning the parameters of the H-PMHT is known to be difficult [17]. Another class of real-time TBD tracking methods is those derived using random finite sets (RFSs). These methods include the (single) Bernoulli filter [14], [15] for single-object tracking and multi-Bernoulli (MB) filters [16], [17] for multiobject tracking. Existing MB filters [16], [17] for multiobject TBD rely on certain heuristics to introduce Bernoulli components for objects that appear in the scene for the first time.
A state-of-the-art conventional tracking method is the trackoriented marginal multi-Bernoulli/Poisson (TOMB/P) filter [5], [6], [8], which models an unknown number of objects and their states by the union of a Poisson RFS and a MB RFS. The TOMB/P filter relies on a high-dimensional marginalization operation performed by means of message passing [6], [18]. A modified version of the TOMB/P filter [19] "recycles" Bernoulli components with low existence probability by transferring them to the Poisson RFS instead of pruning them. This recycling step can significantly improve tracking performance in scenarios with reduced sensor performance [19]. The TOMB/P filter has also been derived and extended to multiple sensors using the framework of factor graphs and the sumproduct algorithm (SPA) [6], [7]. It has been demonstrated that it can outperform existing detect-than-track methods and that it is highly scalable in relevant system parameters [5]- [8].
In this paper, we propose a scalable TBD method that is based on a Poisson/MB model similar to the one of the TOMB/P filter. We first review the general TBD measurement model and discuss a simplification that is suitable for the case where at any discrete time step, the influence of an object is limited to a single "bin" or "pixel" of the raw sensor data. Based on this simplified TBD measurement model, we derive an exact update step. Since the resulting exact posterior probability density function (pdf) is a MB mixture pdf, we approximate it by an MB pdf to limit computational complexity. To guarantee reliable multiobject tracking performance, we recycle Bernoulli components with a low existence probability instead of pruning them. Finally, we verify the performance of the proposed method numerically and compare it to a stateof-the-art TBD tracking method [16].
The main contributions of this paper can be summarized as follows: izations in serif, italic fonts. Vectors and matrices are denoted by bold lowercase and uppercase letters, respectively. For example, x is a random variable, and x is its realization and x is a random vector and x is its realization. Random sets and their realizations are denoted by upright sans serif and calligraphic font, respectively. For example, X is a random set and X is its realization. We denote pdfs by f (·) and probability mass functions (pmfs) by p(·). N (x; µ, Σ) denotes the Gaussian pdf (of random vector x) with mean µ and covariance matrix Σ and R(x; σ) denotes the Rayleigh pdf (of scalar random variable x) with scale parameter σ. The probability of an event is denoted Pr{·}. The symbol ∝ indicates equality up to a normalization factor.

II. RFS FUNDAMENTALS
A RFS X is a random variable whose realizations X are finite sets x (1) , . . . , x (n) of vectors x (i) ∈ R nx . Both the vectors x (i) and their number n = |X| (the cardinality of X) are random. Thus, X consists of a random number n of random vectors x (1) , . . . , x (n) . While the conventional Riemann integral is not defined for sets, one can define the set integral of a realvalued set function g(X ) as [4] g Note that each term of the sum corresponds to one value of the cardinality n = |X |. The statistics of an RFS X can be described by its multiobject pdf f (X ). For any realization X = x (1) , . . . , x (n) , the multiobject pdf f (X ) is given by Here, ρ(n) Pr{|X| = n}, n ∈ N 0 is the pmf of n = |X|, and f n (x (1) , . . . , x (n) ) is a joint pdf of the random vectors x (1) , . . . , x (n) that is invariant to a permutation of its arguments. Note that based on the set integral (1), f (X ) integrates to one.
Next, we will review five classes of RFSs [4]- [6] that are relevant for the derivation of the proposed method. A Poisson RFS X is characterized by a cardinality that is Poisson distributed with mean µ, i.e., ρ(n) = e −µ µ n /n!, n ∈ N 0 and by elements x (1) , . . . , x (n) that are independent and identically distributed (iid) according to the spatial pdf where λ(x) = µf (x) is called the probability hypothesis density (PHD) or intensity function.
A Bernoulli RFS X is represented by an existence probability r and a spatial pdf f (x). It consists of either none or one element with probability 1 − r and r, respectively. According to (2), the multiobject pdf thus reads A MB RFS X is the union of a fixed number J of statistically independent Bernoulli RFSs X (j) , j ∈ {1, . . . , J} with multiobject pdfs f (j) (X ) (cf. (3)) described by the existence probabilities r (j) and the spatial pdfs f (j) (x). The multiobject pdf f MB (X ) of the MB RFS can be obtained by applying the set convolution [4] to the individual Bernoulli pdfs f (j) (X ). For any realization X = {x (1) , . . . , x (n) } with n J, the multiobject pdf f MB (X ) can be evaluated as Here, X (1) ... X (J) =X denotes the sum over all disjoint decompositions of X into sets X (j) , j ∈ {1 . . . , J} such that X (1) ∪ . . . ∪ X (J) = X . For example, for X = x (1) , x (2) and J = 2, the multiobject pdf f MB (X ) can be evaluated as Similarly, for X = x (1) and J = 3, the multiobject pdf f MB (X ) reads Note that for n > J, we have f MB (X ) = 0.
A Poisson/MB RFS X is the union of a Poisson RFS and a MB RFS. The pdf of a Poisson/MB RFS can be obtained by applying the set convolution to the pdf of the Poisson RFS and the pdf of the MB RFS. Let J be the number of components of the MB RFS. For any realization X = {x (1) , . . . , x (n) } with n ∈ N 0 , the multiobject pdf f PMB (X ) can be evaluated as Here, X (0) X (1) = X denotes the sum over all disjoint decompositions of X into two sets X (0) and X (1) such that Finally, a MB mixture RFS is a weighted sum of MB RFSs where without loss of generality, we assume that all MB RFSs have the same number of Bernoulli components J. For any realization X = {x (1) , . . . , x (n) } with n J, the multiobject pdf f MBM (X ) can be evaluated as where I is the number of different MB pdfs and I i = 1 w i = 1. Note that in the implementation of our proposed tracking method, the sums in (4), (5) and (6) are never explicitly evaluated.

III. SYSTEM MODEL
In this section, we describe the system model underlying the proposed algorithm. The multiobject state at time k is represented by an RFS X k = {x . . , M }, that represent the measured intensity of the m'th "bin" or "pixel".

A. State-Transition Model
We use the well-established conventional RFS statetransition model [4]- [6]. At time k − 1, an object with state x k−1 ∈ X k−1 either survives or dies with probabilities p S and 1−p S , respectively. If it survives, its new state x k is distributed according to the single-object state transition pdf f (x k |x k−1 ).
We assume that the states of different objects survive/die and evolve in time independently, i.e., given x k−1 , the single-object state x k is conditionally independent of all the other single-object states x k . Thus, conditioned on the multiobject state X k−1 , the multiobject state of the survived objects X S k , can be modeled by an MB RFS, i.e., k are assumed independent of the newborn objects X B k . Thus, for X k−1 fixed, the overall multiobject state at time k, is given by This model defines the state-transition pdf f (X k |X k−1 ), which can be calculated explicitly via set convolution [4].

B. General Measurement Model
We consider a general superpositional 1 intensity model [17] where the influence of object x This total influence D (m) (X k ) is a sufficient statistic with respect to measurement z n is the variance of the background noise. Conditioned on X k , all measurements are assumed statistically independent. The general joint likelihood function can thus be obtained as Note that according to this general intensity model, every object can potentially contribute to every measurement.

C. Considered Simplified Measurement Model
As an approximation for simplified inference, we consider a model where (i) every object contributes to exactly one distinct measurement z This simplified model is suitable for scenarios where the PSF of each object is highly concentrated around a single measurement, and objects are unlikely to fall in the same pixel.
To describe the associations of objects and measurements at time k, we introduce the vector a k [a k . An object-measurement association vector a k is admissible if each measurement involves contributions of at most one object. All admissible associations form the association alphabet A n,M . For a fixed association vector a k ∈ A n,M , the total influence of the multiobject state X k on measurement z (m) k , m ∈ {1, . . . , M } is thus given by (cf. (7)) . If no object is associated to measurement m, then there is no influence of the multiobject state X k on measurement m.
As we show in Appendix A, the joint likelihood function of z k given X k = {x for n M and according to f (z k |X k ) = 0 for n > M , respectively. Here, M a k = 1, . . . , M \ a consists of the indexes of all measurements that are not associated to any object state x . . , M } for the cases where it is associated to none or one object i ∈ {1, . . . , n}, respectively. For example, following the Swerling 1 model, the pdfs

IV. EXACT ESTIMATION
In the Bayesian sequential estimation framework, the statistics of the state X k at time k, conditioned on all received measurements z 1:k [z 1 . . . z k ] T up to time k, can be described by the posterior pdf f (X k |z 1:k ). This pdf is calculated from the previous posterior pdf f (X k−1 |z 1:k−1 ) via a prediction and an update step. The prediction step calculates the predicted posterior pdf f (X k |z 1:k−1 ) from the previous posterior pdf f (X k−1 |z 1:k−1 ) based on the state-transition pdf f (X k |X k−1 ) discussed in Section III-A. It is identical to the one performed by the TOMB/P filter [5]- [8] and will be skipped. Note that f (X k |z 1:k−1 ) is a Poisson/MB pdf with J k−1 Bernoulli components and intensity function λ k|k−1 (x k ).
The update step calculates the current posterior pdf f (X k |z 1:k ) from the predicted posterior pdf f (X k |z 1:k−1 ) based on the likelihood function f (z k |X k ) in (8). As we show in Appendix B, the updated posterior pdf f (X k |z 1:k ) is no longer a Poisson/MB pdf but a MB mixture pdf (cf. (6)). Contrary to the conventional TOMB/P filter, the updated posterior pdf contains no Poisson part anymore. This is a direct consequence of the fact that the TBD measurement model (cf. Section III-C) does not consider missed detections. However, this is not restrictive for practical scenarios, since objects with very low intensities are also admitted by our model (and thus difficult to distinguish from pure noise). More precisely, as derived in Appendix B, the exact updated posterior pdf can be expressed as where each MB mixture component (i) corresponds to one admissible object-measurement association a k [a ] T , (ii) has J k = J k−1 + M Bernoulli components, and (iii) is weighted by the probability p(a k ). Note that each Bernoulli component represents an object that po-tentially exists and that a k associates each potential object j ∈ {1, . . . , J k } to a measurement m ∈ {1, . . . , M }.
In contrast to the association vector a k used for the derivation of the likelihood function in (8), here a k accounts also for the possible non-existence of objects. More precisely, a k has entries a for a k ∈ A J k ,M and by p(a k ) = 0 for a k / ∈ A J k ,M . In the following, we provide expressions for the existence probabilities r (j,m) k , the spatial pdfs f (j,m) (x k ), and the association weights β (j,m) k used in (11) and (12). For each potential object j ∈ {1, . . . , J k−1 } and measurement m ∈ {1, . . . , M }, we have where we introduced the normalization constant c Here, (13) indicates that the object modeled by Bernoulli component j exists and contributes to measurement m. The state of this object is distributed according to (14) and the likelihood of this event is characterized by (15). Furthermore, for m = 0 we have r Here, r (j,0) k = 0 indicates that the object modeled by Bernoulli component j does not exist. The likelihood of this event is given by (16).
For each new potential object j = J k−1 + m, m ∈ {1, . . . , M }, we get with normalization constant c (17) is the probability that there is a new object that contributes to measurement m (assuming that no other existing object contributed to measurement m). The state of this object is distributed according to (18) and the likelihood of this event is given by (19). Finally, we have β (j,0) k = 1, f (j,0) (x k ) not defined, and r (j,0) k = 0.

V. APPROXIMATE UPDATE STEP
In the following, we will describe the approximate update step used by our TBD method. This update step approximates the exact MB mixture posterior pdf in (11) by a Poisson/MB pdf. This limits computational complexity and, thus, enables real-time tracking of multiple low-observable objects.

A. MB Approximation
First, we approximate the exact MB mixture posterior pdf in (11) by an MB pdf. This approximation is based on approximating the association pmf p(a k ) in (12) by the product of its marginals. We first extend the association alphabet Note thatĀ J k ,M now also contains inadmissible associations. This does not affect p(a k ) since by definition p(a k ) = 0 for A J k ,M \ A J k ,M . Next, we approximate p(a k ) according to with p a Here, ∼ a with j ∈ {1, . . . , J k }\j. Note that a fast and scalable calculation of p a (j) k is enabled by the SPA [5], [6], [18].
Next, we insert (20) into (11), which yields the approximate posterior pdff (X k |z 1:k ) given bỹ By further using the identity By comparing (22) with (4), it can be seen that this approximate posterior pdff (X k |z 1:k ) is now a MB pdf, i.e., it can be rewritten as where the existence probabilities and spatial pdfs of the Bernoulli pdfs f (j) (X (j) k ) for j ∈ {1, . . . , J k−1 } are given by and for j ∈ {J k−1 + 1, . . . , J k } by Note that the idea of approximating a MB mixture pdf by an MB pdf was also used in the derivation of the TOMB/P filter [5] and the SPA-based labeled MB filter [20].

B. Recycling of MB Components
The approximate posterior pdf in (23) consists of J k = J k−1 +M Bernoulli components, i.e., the number of Bernoulli components increases by M at each time step k. This is because newborn objects, modeled by the intensity function λ B (x k ), may appear in the scene and potentially contribute to all of the M measurements. However, many Bernoulli components typically have a very low existence probability and are thus unlikely to represent an existing object.
In order for real-time tracking to remain feasible, the number of Bernoulli components has to be limited. Contrary to most multiobject tracking methods (see e.g. [1]- [4], [6]- [8] ), rather than pruning components with an existence probability below a fixed threshold and potentially discarding valuable tracking information, we employ the concept of recycling [19]. Here, Bernoulli components j ∈ J R k ⊆ {1, . . . , J k } with an existence probability r (j) k below a predefined threshold η R are "transferred" to the Poisson part of the posterior pdf by means of moment matching [19]. This yields the approximated posterior PHD where r (j) k is given by (24) or (26) and f (j) (x k ) by (25) or (27), respectively. After applying this recycling step, the approximate posterior pdf is again a Poisson/MB pdf. In particular, the Poisson part is represented by the approximate posterior PHD λ(x k ) in (28) and the MB part by the existence probabilities r

VI. NUMERICAL STUDY
We consider a two-dimensional (2D) simulation scenario with a region of interest (ROI) of [0m, 64m] × [0m, 64m]. We simulated 10 objects during 200 time steps. The object states consist of 2D position, 2D velocity, and the object's intensity, i.e., x k = [p k γ k ] T with p k = [p k,1 p k,2ṗk,1ṗk,2 ] T . The kinematic part of the object's state p k evolves according to the nearly constant velocity motion model [21] with iid driving noise distributed according to N ( k,p ; 0 4 , 10 −3 I 4 ), and the object's intensity according to a random walk model with iid driving noise distributed according to N ( k,I ; 0, 10 −4 ), respectively. The objects appear at various times before time step 30 and at randomly chosen positions in the area [17m, 48m] × [17m, 48m] at each simulation run, and they disappear at various times after time step 170 or when they leave the ROI. The object's initial velocity is drawn from N (v k ; 0, σ 2 v I 2 ) with variance σ 2 v = 10 −2 . We consider two scenarios. All objects appear with an initial intensity of γ I = 10 and γ I = 4 in scenarios 1 and 2, respectively. While the scenario is not deliberately constructed to cause objects to come into close proximity, this occurs randomly, and the behavior of the compared methods is observed.
The measurement is an image consisting of 64 × 64 cells covering the ROI. Each cell is a square of 1m side length and has a scalar intensity. We use the measurement model described in Section III-B and set the PSF to k is in cell m and according to R z We employ a particle implementation of our proposed TBD algorithm, briefly termed T-TOMB/P due to the close relation to the conventional TOMB/P filter [5] for point measurements. We compare T-TOMB/P to the TBD-based MB filter proposed in [16], briefly referred to as T-MB. Both T-TOMB/P and T-MB represent the spatial pdf of each Bernoulli component by 3,000 particles. T-TOMB/P further represents the posterior PHD by 50,000 particles, the birth PHD by another 50,000, where the resulting 100,000 particles are again reduced to 50,000 after the update step. More precisely, the birth PHD λ B (x k ) = µ B f B (x k ) is represented by drawing particles from the pdf f B (x k ) = f (p k,1 p k,2 )f v ṗ k,1 ,ṗ k,2 f I (γ k ) and by setting µ B = 4/64 2 . Here, f (p k,1 p k,2 ) is uniform over the ROI, f v ṗ k,1 ,ṗ k,2 is N (ṗ k,1 ,ṗ k,2 ; 0 2 , σ 2 v I 2 ), and f I (γ k ) is uniform from 0 to η I = 30, respectively. For the generation of new Bernoulli components of T-MB, we adapt the scheme of [17]. In fact, T-MB generates a new Bernoulli component for each measurement z (m) k−1 whose intensity value is above the threshold η new = 1.5 γ I + σ 2 n . More precisely, the existence probability of each new Bernoulli component is set to 10 −4 and the spatial pdf is represented by drawing particles from Here, the function f z The example shown in Fig. 1 suggests excellent detection and estimation performance of the proposed method. For a quantitative assessment and comparison of the performance of both filters, we computed the mean Euclidean distance based optimal subpattern assignment (MOSPA) metric [22] with cutoff parameter c = 20, order p = 2, and averaged over 1000 simulation Monte Carlo runs. Fig. 2 shows the obtained results of T-TOMB/P and T-MB for γ I = 10 and γ I = 4. It is seen that T-TOMB/P consistently outperforms T-MB. This can be attributed to the excellent behavior of the Bernoulli component generation/recycling scheme of T-TOMB/P and the fact that T-TOMB/P considers data association, i.e., a measurement is allowed to be associated to at most one object and an object to exactly one measurement, respectively. The performance of T-MB depends greatly on the choice of η new . A smaller value of η new results in a faster detection of newly appearing objects, but in a higher number of false tracks, i.e., Bernoulli components modeling not existing objects, on the other hand. Another cause for the performance difference can be understood by considering the case where two objects come close to each other. Since T-MB allows a measurement to be associated to more than one object, one of the two objects can be tracked by both tracks and the other object by none after object separation. This is a direct consequence of the fact that T-MB does not consider data association.

VII. CONCLUSION
We proposed a scalable track-before-detect (TBD) method for the tracking of low-observable objects that relies on a Poisson/multi-Bernoulli model. To limit the computational complexity, we approximated the exact posterior pdf by a multi-Bernoulli pdf. For a reliable tracking in real time, a recycling of Bernoulli components is performed. We demonstrated that with the proposed method, a significant improvement in tracking performance can be achieved compared to a state-ofthe-art TBD tracking method. A possible direction for future research is an extension of the proposed method to a more general measurement model by using sum-product algorithms for data association with extended objects [23], [24].

APPENDIX A
In the following, we derive (8). Based on assumptions (i) and (ii) in Section III-C, the joint likelihood function for X k = {x (1) k , . . . , x (n) k } can be found as [25] f (z k |X k ) ∝ By further inserting (28), we get Expression (29) can be reformulated by using the objectmeasurement association vector a k introduced in Section III-C, which yields the final expression for the joint likelihood function in (8).
APPENDIX B In this appendix, we derive the exact posterior pdf f (X k |z 1:k ) in (11). First, we obtain the posterior pdf by applying Bayes theorem [4] according to Next, we reformulate the likelihood function in (8) as follows. We define the measurement-object association vector . . , n}, which expresses the same information as a k but in a different form [6], [18]. In particular, b (m) k = i ∈ {1, . . . , n} indicates that object i contributes to measurement m and b k = 0 indicates that no object contributes to measurement m. An association b k is admissible if any object contributes to exactly one measurement m, and at most one object contributes to each measurement. All admissible measurement-object associations form the set B M,n .
The combination of both a k and b k is now used to reformulate (8) for each X k = {x with n M . Here, X k is decomposed into arbitrary subsets k,1 } with X k = X k,0 ∪ X k,1 and n = n 0 + n 1 . Further, g z where the Bernoulli pdfs f (j) (X k ) are parametrized by the existence probabilities r (j) k|k−1 and spatial pdfs f (j) k|k−1 (x k ) and the Poisson pdf f P (X k,0 ) is represented by the posterior PHD λ k|k−1 (x k ). By plugging (31) and (32) into (30) and performing certain reformulations, we obtain [25] f (X k |z 1:k ) ∝ X k,0 X k,1 =X k ā k ∈Ā J k−1 ,M h 1 (X Here the association vectorā k [ā In addition, the function h 1 (X k,1 ,ā k ) is given by (34) Here, f (j) (X (j) k,1 ) is the Bernoulli pdf in (32) and k,1 is obtained as Expression (34) is a product of weighted Bernoulli pdfs with parameters as in (13)- (16). Furthermore, the function h 0 (X k,0 , b k ,ā k ) reads h 0 (X (1) k,0 , . . . , X (M ) is the Poisson pdf in (32). Expression (35) is a product of weighted Bernoulli pdfs with the parameters as in (18)- (19). By using the association vector a k as defined in Section IV, we can now reformulate (33) according to Finally, we introduce the association pmf p(a k ) (cf. (12)) and obtain the final expression for the posterior pdf in (11).