Robust Received Signal Strength Indicator (RSSI)-Based Multitarget Localization via Gaussian Process Regression

We consider the robust localization, via Gaussian process regression (GPR), of multiple transmitters/targets based on received signal strength indicator (RSSI) data collected by fixed sensors distributed in the environment. For such a scenario and approach, we contribute both with a novel noise robust procedure to train the parameters of the GPR model, which is achieved via a mini-batch stochastic gradient descent (SGD) scheme with gradients given in closed form, and with a pair of corresponding robust marginalization procedures for the estimation of target locations. Simulation results validate the contributions by showing that the proposed methods significantly outperform the best related state-of-the-art (SotA) alternative and approach the performance of a genie-aided (GA) scheme.


I. INTRODUCTION
W IRELESS localization technology, both indoors [2] and outdoors [3], has gained great attention in the last few decades as an essential enabling technology for various systems and applications ranging from vehicular localization and tracking [4], to near-field radio frequency identification (RFID) positioning [5], to Internet-of-Things (IoT) networks in smart factories/homes [6] and more.
Traditional localization techniques have typically relied on global navigation satellite systems (GNSS) for the satellitebased geolocation and time information, but such methods exhibit poor power-efficiency, precision, latency, and robustness in dense urban scenarios, especially in comparison to the requirements of beyond fifth-generation (B5G) applications.More recent research in this area has, therefore, focused on the utilization of wireless signals from the immediate environment, relying on wireless sensors, devices and access points to extract information such as time of arrival (ToA), angle of arrival (AoA), Niclas Führling, Hyeon Seok Rou, and Giuseppe Thadeu Freitas de Abreu are with the School of Computer Science and Engineering, Constructor University (previously Jacobs University Bremen), 28759 Bremen, Germany (e-mail: nfuehrling@constructor.university; hrou@constructor.university;gabreu@constructor.university).
Digital Object Identifier 10.1109/JISPIN.2023.3332033 and received signal strength indicator (RSSI) in order to estimate their locations.
Due to the distributed nature of such localization scenarios, resulting algorithms generally consist of multidimensional, multivariate optimization problems based on received signal metrics obtained at the sensors.And in order to alleviate the challenges and complexity of solving such optimization problems, recent state-of-the-art (SotA) methods have proposed solutions leveraging machine learning (ML) techniques.
To cite some examples, in [7], a supervised ML method was used to evaluate Taylor series-linearized time difference of arrival (TDOA) measurements and localize autonomous unmanned aerial vehicles.In [8], a deep residual network was used to process AoA measurements in a singlesensor non-line-of-sight (NLOS) indoor localization scheme.And a method to incorporate heterogeneous sensor information was proposed in [9], where a K-nearest-neighbor (KNN) algorithm was used to adaptively select and combine various signal radio frequency (RF) features for indoor localization.
Besides (and partially thanks to) the recent popularity of ML approaches, an increasing trend is to operate with RSSI only, which has the advantage of not requiring the embedding into wireless signals information tailored to the positioning application itself, leading to added flexibility that enables the design of multifunctional systems [10].Three excellent examples are the work in [11], in which multitarget RSSI localization over a discrete grid is initially estimated via a sparse dictionary updating and a K-means clustering and later refined via dictionary updates; the work in [12], where the ML is used to refine the RSSI data itself; and the scheme of [13], where the location of multiple targets are obtained from RSSI values based on a GPR method.Inspired by the aforementioned works, we propose in this article a noise-robust improvement of the RSSI-based localization method via GPR first proposed in [13], offering in particular the following contributions.
1) A noise-robust mini-batch SGD method to train the GPR model is introduced, which removes the requirement of error-free training data.2) Closed-form expressions of the Hessian of the objective optimized during the training stage of the algorithm are derived, which reduces the complexity of implementation by avoiding numerical evaluation of matrix derivatives.3) A noise-robust variation of the RSSI-based localization method via GPR first proposed in [13] is given which, besides training, also incorporate robustness both in the computation of covariance matrices of the location estimates.Comparisons of the proposed algorithm against the SotA method confirm that independent improvements due to the robust training and the robust location estimate computation are achieved, which combined yield substantial gains over the preceding technique.

II. SIGNAL MODEL AND LOCALIZATION SCHEME
Consider a multi-user (MU) distributed massive multipleinput multiple-output (DM-MIMO) scenario, consisting of K single-antenna transmit devices (target nodes) which are served by M single-antenna sensor nodes centrally connected to a computing unit (CU) via error-free fronthaul links with infinite rate, as depicted in Fig. 1.
The 2-D positions of the mth sensor and kth target nodes, with m ∈ M {1, . . ., M} and k ∈ K {1, . . ., K}, are respectively described by the vectors where (x SN m , y SN m ) and (x TG k , y TG k ) are respectively the xand ycoordinates of the mth sensor and kth target nodes.
Following the above, the Euclidean distance between the mth sensor node and the kth target node is given by A. RSSI Model The received signal vector at the mth sensor node over T consecutive transmission instances is given by where h mk ∈ C is the channel gain from the mth sensor node to the kth transmitter, assumed to remain constant during T transmissions, ρ k ∈ R and ω k ∈ C T ×1 are respectively the transmit power 1 and symbol vector of the kth transmitter, while is the additive white Gaussian noise (AWGN) vector at the m-th sensor node.
The flat-fading uplink channel gain coefficient h mk in (3) is further modeled by where q mk ∼ CN (0, 1) is the small-scale fading factor, and g mk is the large-scale fading factor given by where b 0 ∈ R and η ∈ R are respectively the reference path loss coefficient and the path loss exponent determined by environmental assumptions, and d mk ∈ R and z mk ∼ N 0, σ 2 z mk are respectively the distance defined in (2) and the channel gain due to random shadowing of the path between the mth sensor node and the kth transmitter.
It follows from the above that the RSSI of each kth transmitter at the mth sensor is given by which can be estimated -either semiblindly, by leveraging short orthogonal pilot sequences together with independent payload data, as shown in [13] and [14]; or blindly by exploiting the sparsity resulting from intermittent user activity, such that only a relatively small random subset of transmitters share the channel at each transmission instance, as shown in [15] -with basis on the relation where T respectively denote a matrix collecting the transmit signals and the powers from all K transmitters at the mth sensor.
For future convenience, we shall express the powers p mk in decibel (dB), namely where p dB Finally, with pardon for the slight abuse of notation, the vector containing the RSSI (in dB) from all the transmitters aggregated from all the M sensors at the central processing unit (CPU) can be defined as

B. Location Estimation Via GPR
Following [13], consider now arbitrary functions f x (•) and f y (•) that map the received signal power vector of any kth transmitter to its xand y-coordiantes, respectively, i.e., where it is implied that if the functions f x (•) and f y (•) are known at the CU, the positions of any transmitter can be obtained, given the aggregated RSSI.
Remark: As the expressions for the xand y-coordinates are identical, we hereafter omit the subscripts (•) x and (•) y in corresponding expressions, such that the derivations to follow apply to both xand y-coordinates, unless stated otherwise.Occasionally, we shall use the letter c to indicate either, interchangeably.
Under the GPR approach, the latter coordinate mapping function is drawn from a Gaussian process (GP) prior, i.e., where GP(0, Φ) denotes the GP prior with zero mean and covariance matrix Φ ∈ R K×K , whose element at the kth row and qth column is given by the covariance function φ(p k , p q ), dependent on the RSSI between the kth and qth transmitters.
To accurately model the relationship of the RSSI and the coordinates, a suitable covariance function must be designed to carefully capture the differences in power between different transmitters in the region of interest (ROI), which as argued in [13] is given by (12) where p k ∈ R M ×1 and p q ∈ R M ×1 are respectively the RSSI vectors of the kth and qth transmitters, σ2 er ∈ R is the variance of the measurement error, δ kq is the indicator variable taking on the value 1 if k = q and 0 if k = q, and the parameters α, B, γ are the weight parameters to be learned, with where β m ∈ R ∀m ∈ M are the weights corresponding to each RSSI vector in the exponential term of the covariance function defined in (12).For future convenience, let us define the vector collecting the M + 2 unknown weight parameters of the model, namely Let us also define the sets of coordinates c ∈ R K×1 corresponding to K training locations, and the associated set of training RSSI vectors p k ∈ R M ×1 , with k = 1, • • • K. Accordingly, we define as well the training covariance matrix Ψ ∈ R K×K associated with the K training locations and corresponding RSSI values, and the cross-covariance matrix Σ ∈ R K×K between the K training and the K target locations, constructed from the corresponding sets of RSSI values, respectively.Then, the joint distribution of the coordinate vectors of the training and target locations, according to a conventional GP approach, is given by, where the kth row and kth column of Σ ∈ R K×K is given by φ(p k , p k ), as obtained from (12).
In possession of the joint distribution in (15), the conditional distribution of c is given by where the conditional mean vector μ ∈ R K×1 and covariance matrix C ∈ R K×K is obtained via Finally, the marginal distribution of the individual target node coordinate can be obtained as where c k and ν k are respectively the marginal mean and variance of the estimate of the coordinate of the kth target, given by where [ • ] k and [ • ] k,i respectively denote the kth element of a vector and the matrix element at the kth row and ith column.Since c k is Gaussian distributed, the marginalized mean c k is directly the maximum-a-posteriori estimate of x k , and hence the final predicted coordinate of the kth target.
It can be inferred from the procedure described above that the performance of the GPR-based RSSI localization algorithm is highly dependent on the accuracy of the covariance matrix model given in (12), which in turn is fundamentally determined by the parameter vector θ.In other words, a core step of the method is to optimally determine θ, given a certain amount of training data.
In [13], this is achieved via a maximum likelihood fitting of the distribution of the training coordinates c, given the corresponding set of K RSSI vectors p k , assumed to be free of errors.Such an assumption is not only virtually impossible to meet in practice, 2 but also a cause of poor performance in real-life applications, since the inevitable presence of errors in RSSI values collected is not mitigated by design.
In view of limitation, we introduce in the next subsection the key contribution of the article, which can now be stated clearly as a method to determine the parameter vector θ using noisy training data, lending both feasibility and robustness to the overall approach.

III. PROPOSED METHOD: ROBUST RSSI-BASED LOCALIZATION VIA GPR
In this section, we propose a noise-robust training solution to improve the multitarget localization problem described above, which effectively translates to an optimization problem over the K-pairs of real-valued 2-D coordinates, given the aggregated RSSI at the M sensors of known positions, solved via SGD.

A. Robust Model Training Method: Fundamentals
Let us first revisit the relationship between the information contained in the RSSI vectors and the coordinates of the training points.In particular, unlike SotA methods [12], [13] in which training RSSI values assumed to match perfectly with the model in ( 4)-( 6) are required, we accommodate for the fact that random shadowing noise affects also training data.In other words, similarly to measured target RSSI values utilized in the online (localization) stage of the algorithm, the training RSSI vectors used to obtain the optimal parameter vector θ are also modeled as with z mk ∼ N 0, σ 2 z mk .This is equivalent to treating the coordinates c of the training points not as deterministic, but rather as random variables with a distribution that, under the GPR model can be assumed to be approximated by where the elements of the covariance matrix Ψ ∈ R K×K are determined via (12), that is Given the distribution in (21), the optimal weight parameters θ can then be determined as the solution of the maximum likelihood problem where we implicitly defined the objective g(θ; p 1 , . . ., p K ), which for the purpose of training is a function of the parameter vector θ, as highlighted by the notation, and which in view of the Gaussian approximation in (21), can be modeled as where | • | denotes the determinant of the argument matrix.
It is evident from (12), (22), and (24), that g(θ; p 1 , . . ., p K ) is not convex on the optimization variable θ.In fact, the determinant and inversion operations onto Ψ, and the nonlinear dependence of the latter on the parameters β m gathered in the matrix B, make the problem highly intractable from the perspective of optimization theory [16], [17].In [13], the authors turn their attention to other matters and leave the issue largely unaddressed, limiting the discussion to citing relevant gradient-based methods and classic literature [18].
We therefore address this issue here via a SGD approach [19].To this end, first, consider the partial derivative of g(θ; p 1 , . . ., p K ) with respect to α, which is given by Next, we leverage the matrix derivative identities [20] (26) to simplify the expression in (25) to Observing that the partial derivatives of ( 24) with respect to β m and γ are identical in form to the expression in (27), we omit details of the derivations corresponding to the latter parameters and limit ourselves to providing expressions for each element of the partial derivatives of Ψ with respect to α, β m , and γ, which are respectively given by where The gradient of the objective function in (24) can then finally be put together, yielding (29) such that the optimal parameter vector θ can be optimized via gradient descent (GD), according to the update equation where θ (i−1) and θ (i) are the parameter vectors at the (i−1)th and ith SGD iterations, respectively, while λ is the learning rate, which for convergence guarantees is bounded by the Lipschitz condition with L denoting the Lipschitz constant, given by where max eig(•) denotes the operator returning the maximum eigenvalue of the argument matrix.
It is often the case, that the training procedure described above is performed over a large data set, consisting of multiple snapshots of data collected for a large number of training locations (see e.g., Fig. 2).And since the complexity of the training algorithm is fundamentally determined by the inversion of the covariance matrix Ψ, as evident from (29), a mini-batch SGD-based variation of the training scheme outlined above is described in the sequel, with aim at lowering the complexity of its implementation [21].Next, consider the covariance matrix with structure similar to that of (22), but constructed only with respect to the training location vector cb , i.e.,

B. Robust Model Training Method: Mini-batch SGD Approach
Then, under the SGD method, the update of the parameter vector θ for the mini-batch D b is obtained by the corresponding variations of ( 29) and (31), namely A schematic illustration of the noise-robust parameter θ optimization method via mini-batch SGD is given in Fig. 3, and the corresponding pseudocode is offered in Algorithm 1. Recall that the algorithm is to be run independently for each coordinate Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.x and y, such that its execution yields the optimized (trained) parameters θ * x and θ * y .

C. Robust RSSI-Based Localization GPR
In possession of a parameter vector θ obtained through model training, the SotA RSSI-based localization algorithm via GPR of [13] reduces to evaluating (19a) for each kth target, given the corresponding RSSI measurements p k obtained online, the noise-free training RSSI vectors p k obtained offline, and the associated covariance matrix Ψ.
Besides the parameter vector θ * obtained through the noiserobust approach described above, another distinction between our proposed robust method and the SotA is that instead of noise-free training RSSI vectors we have an entire data-cube D containing multiple snapshots of noisy training RSSI vectors p (s) k , such that different alternatives exist to computing the GPRbased location estimates.One alternative would be, for instance, to treat each noisy RSSI vector p (s) k in the data cube D as a distinct noise-free data point and construct the corresponding KS × KS covariance matrix described by where (36) Then, stacking the training location vector c onto itself S times, which can be concisely represented via the Kronecker product 1 S ⊗ c, the target locations could be obtained from the following variation of (19a) Clearly, the complexity of evaluating (37) is of order O((KS) 3 ), which becomes quickly prohibitive as the number  of training points K and/or snapshots S grows, such that this approach, described here only for completeness, is not recommended and will not be pursued any further hereafter.
A lower-complexity alternative to the latter approach can be obtained by performing location estimation independently using each snapshot of the data cube and averaging the results, which will be hereafter referred to as the "estimate-then-average" (EA) approach and can be concisely expressed as Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.The approach summarized by (38) requires the inversion of S snapshot K × K covariance matrices Ψ(s,s) , such that its complexity is of order O(S • K 3 ), which although substantially lower than that of (37), may still be too high for large S.
In that case, an alternative of lower complexity -of order O(K 3 ) -which will be hereafter referred to as the "averagethen-estimate" (AE) approach, is to utilize the sample covariance matrix and compute the location estimates via Output: Estimated target x-, y-coordinate vector c (AE) or c (EA) .Repeat ∀k ∈ K: The proposed RSSI-based localization scheme via the AE-GPR and EA-GPR approaches are summarized in the form of pseudo-code in Algorithm 2(a) and (b), respectively.

IV. SIMULATION RESULTS
In this section, the effectiveness of the proposed SGD-based robust training mechanism and the corresponding RSSI-based robust localization method via EA/AE-GPR are evaluated via computer simulations and comparisons the SotA [13].
Our first results are given in Fig. 4, which compares the average convergence behavior, as a function of GD epochs, of the proposed mini-batch SGD-based robust training algorithm with different mini-batch sizes K.In order to keep the comparison fair in terms of the total amount of training data used to optimize the parameter vector θ, the number of mini-batches B is varied so that K • B is the same for all curves.
As expected, the results show that smaller mini-batch sizes lead to slower convergence but also to lower points of the objective function, while larger mini-batches yield faster convergence but to higher local minima, which is a typical tradeoff of convergence speed and optimality found in mini-batch-based SGD schemes [21].
Next, we turn our attention to the proposed EA/AE-GPR localization algorithms themselves, comparing their performances to that of the SotA RSSI-based GPR scheme of [13].
The considered scenario is a street intersection containing K = 3 targets, with a ROI of 100 m × 100 m, serviced by M = 56 sensors placed along the roadsides.Training is Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

TABLE I SIMULATION PARAMETERS
conducted over a grid of K = 204 training locations distributed in a regular grid within the area, as illustrated in Fig. 5.
Environment parameters such as the reference path loss coefficient and exponent are set according to the 3GPP urban micro propagation model described in [22]; radio parameters such as transmit/noise powers and radio sensitivity are set according to the LTE standard [23]; and finally SGD parameters are set in light of the results of Fig. 4, which have shown that the optimized objective function reaches lower convergence points with minibatches of smaller sizes.These parameters are summarized in Table I, and will be consistently utilized hereafter, unless when otherwise stated.
Our first results are offered in Fig. 5, which shows the estimates obtained from multiple realizations of the proposed and approaches.To the advantage of the SotA technique [13], and for the sake of visibility, the proposed method is represented only by the lower-complexity robust AE-GPR variation, which shall be later slightly worse performance than EA-GPR alternative.It can be seen from the cloud of points formed by the location estimates that the proposed robust approach is indeed superior to the SotA.
Next, in order to offer a more quantitative assessment of the gains achieved by the proposed schemes over the SotA method, we compare their performances in terms of the location root mean square error (RMSE), defined as For the sake of providing a lower-bounding reference, results corresponding to the location estimates c (GA) obtained via an ideal GA GPR scheme are also included, which are computed via the expression where pk are noise-free RSSI vectors generated via and Ψ is the corresponding noise-free covariance matrix constructed using the optimized parameter vector θ * obtained with the proposed noise-robust training method described in Section III-B, which is given as The results are given in Fig. 6, and clearly show a wide gap between the proposed robust EA/AE-GPR methods and the SotA technique of [13], with the alternative of lowest computational complexity, namely the AE-GPR scheme, slightly outperformed Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.by the also slightly more computationally expensive EA-GPR method.Interestingly, it is also found that both the AE-GPR and EA-GPR robust proposed algorithms achieve RMSE performances that come very close to that of the (unfeasible) genie-aided lower-bounding reference.
Since the GA-GPR scheme employs the parameter vector θ * obtained with the proposed noise-robust training method of Section III-B, the results indicate that most of the gains achieved are in fact due to the proposed SGD-based training approach, with the remaining gap between the EA/AE-GPR methods and the GA-GPR reference due merely to the available amount of data.
The latter finding, namely, that both the AE-GPR and EA-GPR robust proposed algorithms asymptotically approach the same performance of the genie-aided is corroborated by the results of Fig. 6(b), which show that indeed the performance of these schemes are very similar when larger mini-batch sizes K are employed.
To further evaluate the proposed methods, the CDFs of the RMSE of location estimates are given in Fig. 7, for both K = 5 and K = 50, and with a noise variance of σ 2 z = 1.We emphasize that choice of σ 2 z is to the disadvantage of the proposed scheme, since the gains of the latter over the SotA alternative were shown in Fig. 6 to increase with that parameter.Still, the results clearly show that the improvement achieved with the proposed method is very significant not only in average terms, but rather for the entire range of RMSE values.It is particularly remarkable that with larger minibatch sizes the performance of the proposed methods at the low RMSE range overlaps with that of the ideal GA reference.
Finally, we assess the impact of the number of sensors M onto the performance of the GPR-based methods.To that end, in order to isolate the effect of M from other factors such as the geometry of the scenario [24], we test the schemes under a highly symmetric single-target scenario as depicted in Fig. 8, where the target is located at the center of the ROI, and the sensor are positioned symmetrically in the surroundings, as they if they were installed on the roadsides of a crossing.
The results, given in terms of the average RMSE as a function of the number of sensors M , are shown in Fig. 9, and again demonstrate that the proposed methods, not only consistently outperforms the SotA, but also tend to approach the GA performance -in particular the EA-GPR scheme of Algorithm 2(b) -as the number of sensors increases.
All in all, simulations thoroughly corroborate the improvement achieved by the proposed methods.

V. CONCLUSION
We proposed a mini-batch SGD-based noise-robust training procedure to optimize the parameters of a GPR model for the RSSI-based multitarget localization method, along with a pair (AE and EA) of marginalization procedures for the computation of corresponding location estimates.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
Simulation results demonstrated that the proposed robust AE-GPR and EA-GPR schemes significantly outperform the best SotA alternative utilizing a similar GPR-based method, while approaching the performance of an GA scheme that employs the proposed training method combined with an unfeasible marginalization procedures based on the true covariance matrix of the RSSI data.

Manuscript received 14
July 2023; revised 9 October 2023; accepted 6 November 2023.Date of publication 10 November 2023; date of current version 29 November 2023.This work was supported by Continental AG.An earlier version of this paper was presented in part at the 2023 IEEE International Workshop on Computational Advances in Multi-Sensor Adaptive Processing Search form Search (CAMSAP 2023).(Corresponding author: David González G.)

Fig. 1 .
Fig. 1.Distributed massive MIMO scenario considered in the proposed multi-user RSSI-based localization method.

Fig. 2 .
Fig. 2. Illustration of a data cube D corresponding to a case where the number of sensors is M = 5, the number of training locations is K = 7, the total number of snapshots of RSSI vectors per training location is S = 8.An example of a mini-batch D b with K = 4 is also illustrated (gray columns), which consists of a random selection of RSSI vectors corresponding to the training locations c = {x 1 , x 2 , x 4 , x 7 }, such that k = {1, 2, 4, 7}, respectively taken at snapshots s = {2, 8, 5, 8}, yielding D b = p (2)1 , p(8)  2 , p(5)  4 , p(8)   7

Fig. 4 .
Fig. 4. Convergence of proposed robust SGD-based training with different mini-batch sizes K but same total amount of training data ( K • B = 1000, identical for all curves).

Fig. 5 .
Fig. 5. Snapshot the multitarget location estimates obtained with proposed robust and SotA RSSI-based GPR schemes.Measurement noise variance is σ 2 z = 1, and the proposed AE-GPR algorithm employs a parameter vector θ * trained as described in Section III-B with minibatches of size K = 5 and B = 200.

Algorithm 2 :
Robust RSSI-based Localization via GPR.Input: Data cube D ∈ R M ×K×S with S snapshots of each and all RSSI vectors p (s) k , with k ∈ {1, • • •, K} and s ∈ {1, • • •, S}; Measured RSSI values p k from all targets; and Optimized parameter vector θ * .

]
T denotes the true coordinates of the kth target, while and [x Est k , y Est k ] T are the corresponding estimates obtained by the respective algorithm.

Fig. 9 .
Fig. 9. RMSE of the estimated coordinates of the SotA and the proposed algorithm, with varying localization methods, for a varying amount of sensors M , for K = 5 and σ 2 z = 1.