Wireless localization for mmWave networks in urban environments

Millimeter wave (mmWave) technology is expected to be a major component of 5G wireless networks. Ultra-wide bandwidths of mmWave signals and the possibility of utilizing large number of antennas at the transmitter and the receiver allow accurate identification of multipath components in temporal and angular domains, making mmWave systems advantageous for localization applications. In this paper, we analyze the performance of a two-step mmWave localization approach that can utilize time-of-arrival, angle-of-arrival, and angle-of-departure from multiple nodes in an urban environment with both line-of-sight (LOS) and non-LOS (NLOS) links. Networks with/without radio-environmental mapping (REM) are considered, where a network with REM is able to localize nearby scatterers. Estimation of a UE location is challenging due to large numbers of local optima in the likelihood function. To address this problem, a gradient-assisted particle filter (GAPF) estimator is proposed to accurately estimate a user equipment (UE) location as well as the locations of nearby scatterers. Monte-Carlo simulations show that the GAPF estimator performance matches the Cramer-Rao bound (CRB). The estimator is also used to create a REM. It is seen that significant localization gains can be achieved by increasing beam directionality or by utilizing REM.


Introduction
The demand for wireless broadband communication has been growing rapidly, which has been the driving force for the emergence of 5G cellular networks. It has recently been shown in the literature that millimeter wave (mmWave) technology is not only feasible for dynamic outdoor cellular networks, but can facilitate a thousand fold increase in data capacity [1][2][3][4]. The mmWave cellular networks are expected to first be deployed in dense urban environments where the Global Positioning System (GPS) signal may typically be unavailable, and the demand for large data rates is high. With coverage ranges that can extend to hundreds of meters, the network must be solely responsible for localization, while also simultaneously achieving high data rates in such scenarios [1]. Communication performance in 5G networks will be limited by the amount of time required to align the highly directional beams of the communicating nodes, particularly for exhaustive beam searches, which are costly to capacity [5]. A network that is able to localize node positions can significantly reduce the time spent on beam alignment and increase capacity [6]. Thus, it is important to characterize the performance of mmWave network localization for urban scenarios.
Wireless localization with strictly non-line-of-sight (NLOS) paths is achievable for omni-directional antennas by exploiting the time-of-arrival (TOA), angle-of-arrival (AOA), and angle-of-departure (AOD) measurements [7][8][9]. While accurately measuring AOA and AOD is relatively difficult at lower frequencies due to the rich scattering and poor path separability, this is much easier to achieve in mmWave channels leading to improved localization performance [10]. The mmWave frequencies also allow the use of ultra-wide bandwidths larger than 1 GHz, which helps in providing precise TOA estimates. Figure 1 highlights channel characteristics of mmWave frequencies that offer advantages over traditional microwave frequencies for localization.

Literature review on mmWave localization
There have been recent studies in the literature that evaluate mmWave localization performance in various scenarios. Localization with received signal strength (RSS), TOA, and AOA are analyzed in [11] (separately and jointly for different measurement parameters) for LOS and NLOS scenarios, showing promising results with TOA and AOA, and less reliable results with RSS. A log-normal path loss model is used to evaluate RSS, time-difference-of-arrival (TDOA), and AOA localization methods for LOS paths in [12]. A mobile's location and orientation are estimated jointly in [13,14] for mmWave systems. It is shown in these papers that a single fixed equipment (FE) is sufficient to localize a user equipment (UE), but only LOS paths are considered. A direct localization approach for a user connected with multiple FEs is introduced in [15], but NLOS paths are treated as interference and LOS paths are still required.
At lower frequencies NLOS paths are treated as interference. A major advantage of mmWave frequencies is that very few paths have significant received signal strength, which results in channel sparseness [16]. This property enables NLOS paths to no longer be considered as interference, but to instead be exploited as paths with additional information that can improve localization performance [17]. The work in [10] shows that it is possible to use LOS and NLOS paths to determine the orientation and position of a node communicating with a single transmitter at mmWave frequencies under certain conditions. Specifically, it is shown that sufficient conditions for position and orientation estimation require at least one LOS path or three NLOS paths. The Cramer-Rao bound for localization and orientation is derived, and a localization algorithm is proposed that exploits channel sparseness to estimate AOD/AOA/TOA, which is then used to estimate user position and orientation. The work in [18] extends the Cramer-Rao bound to three dimensions.

Summary of the proposed mmWave localization technique
The algorithm proposed in [10] is only suited for single transmitter scenarios. However, there may be scenarios where a single transmitter is unable to establish enough paths to meet the sufficient conditions for localization. In this work, we consider scenarios that are not necessarily limited to a single transmitter. Instead, an initial access or beam alignment stage is used to obtain rough AOD/AOA/TOA estimates for LOS and NLOS paths from one or multiple transmitters. Then, the user position is estimated using the AOD/AOA/TOA estimates. This step requires optimizing a non-convex cost function with many local maxima. To accomplish this, we propose a gradient-assisted particle filter (GAPF) estimator. While separately estimating TOA, AOD, and AOA is suboptimal, it significantly reduces computational effort and allows links from multiple transmitting FEs to be used for localization.
It should be noted that the work in [10] studies methods to jointly estimate and refine AOD/AOA/TOA for LOS and NLOS paths from a single FE. This process can be used to obtain improved AOD/AOA/TOA estimates for one FE at a time in place of the proposed reduced complexity first step. Then, these estimates can be combined and used in the second step for improved localization accuracy over the proposed reduced complexity approach. The reduced complexity approach is only considered in this work because the minimal processing approach to estimating AOD/AOA/TOA greatly reduces computation time and still achieves adequate localization accuracy. The proposed method is best suited for environments where many nodes are connected, such as a 5G network with many FEs in an urban environment where a single FE may not be able to establish enough paths to meet the sufficient conditions for localization.
Urban environments with many buildings and large structures contain scatterers that remain fixed in space. The three dimensional spatial characteristics of the urban environment can be captured by radio-environmental mapping (REM) and later exploited to estimate scatterer locations for NLOS paths to improve localization [19]. The knowledge provided by localization and REM can be used to relax initial access requirements and improve capacity for 5G communication systems [20]. In this work, we study REM-assisted localization, which assumes scatterer locations are estimated a priori and can be used to improve localization performance.

Methods/experiments and contributions of this work
To our best knowledge, mmWave network localization that exploits LOS and NLOS paths from multiple FEs in an urban environment utilizing TOA/AOA/AOD with/without REM has not been studied in the literature. This paper analyzes the localization performance of 5G mmWave communication networks specifically considering urban canyon and urban corner environments where directional beams are captured at the UE in the downlink from one or multiple FEs. The contributions of this paper are summarized as follows: • A reduced complexity two-step localization scheme is introduced where in the first step the receiver obtains rough estimates of AOD/AOA/TOA for LOS and NLOS paths from initial access or beam alignment stages from multiple FEs [10,21]. The AOD estimate is obtained from a transmitted reference signal where the transmitter embeds the center beam direction and the AOA/TOA are estimated non-jointly by the receiver. The second step uses the AOD/AOA/TOA estimates to estimate the receiver and NLOS scatterer location coordinates. • A gradient-assisted particle filter (GAPF) estimator is proposed as a maximum likelihood (ML) estimator to estimate the UE position and scatterer coordinates over a non-convex space. It is shown to have performance that matches the Cramer-Rao bound (CRB) through Monte-Carlo simulations. • Localization performance is analyzed in urban canyon and urban corner scenarios where a UE is connected with one FE or two FEs. • The localization accuracy of REM-assisted and non-REM-assisted network performance is analyzed. The performance of a perfect REM system that has perfect knowledge of a scatterer locations (with no localization error) for each observed path is used to bound realistic REM systems.
• The scatterer locations that are extracted from the proposed localization approach are used to create an REM.
The rest of this paper is organized as follows. Section 2 introduces the mmWave localization model and studies how LOS and NLOS paths can be separated for localization purposes. Section 3 discusses estimation methods and introduces the GAPF estimator for localization. Section 4 derives the CRB for mmWave localization performance. Section 5 compares estimator performance to the CRB, analyzes mmWave localization performance in urban canyon/corner environments with/without REM, and examines the trade-off between system complexity and localization accuracy. Subsequently, the proposed estimator is used to create an REM. Finally, Section 6 provides concluding remarks.

mmWave localization system model
This section introduces the downlink system model in a mmWave network from multiple FE nodes to a single UE node, where the UE estimates its own position. Without loss of generality, our results can also be extended to network-based localization relying on uplink mmWave signals from the UE. First, the 2D system model for mmWave urban localization is introduced. Then, differentiating LOS from NLOS paths is studied. Additionally, raytracing techniques are used to further study the effects of directional beams at mmWave frequencies and justify that the results in the included 2D model are representative of a full 3D model.

Localization model
While the position of a UE in a wireless network can be estimated directly, it is often more practical to implement a two-step positioning approach that first determines a set of parameters such as TOA, AOD, or AOA, which are then used to estimate the position [22].
The considered scenario utilizes the periodic beam training stage or initial access for mmWave network communications to collect AOA/AOD/TOA for a variety of paths [6,23]. The beam training stage searches over possible beams at the FE and UE where each directional beam points in a different direction. The beams with the strongest signals are then used to estimate AOA/AOD/TOA for the paths associated with those beams. A propagation channel is known to be sparse at mmWave frequencies [16,23], for which one major reason is the use of highly directional beams enabled at mmWave frequencies with many antenna elements. Hence, each beam will typically contain a LOS path or a strong single bounce NLOS path [16,24] along with multiple bounce NLOS paths that will have large attenuation and much lower signal strength [23]. Since only the strongest beams for each UE/FE pair are selected, multiple bounce paths are relatively unlikely, which allows the model to be simplified to LOS paths and single bounce NLOS paths [23].
Assuming a synchronized network, ultra-wide bandwidths of mmWave frequencies provide precise TOA estimates [25]. Measuring AOA and AOD is not feasible at lower frequencies because of large numbers of scattered paths. However, the arrays with large numbers of antenna elements at mmWave frequencies easily fit on chip enabling highly directional beams. Arrays at the UE provide precise AOA measurements, where it is assumed, the orientation of the UE is known so that the AOA relative to the overall coordinate system is determined from the AOA relative to the UE array. An AOD measurement is obtained from the FE, which transmits the beam's quantized AOD relative to the overall coordinate system. This is easily calculated from the AOD of the FE relative to the antenna array under the assumption that the FE orientation is known. Additionally, the FE is assumed to broadcast its position coordinates. These low-complexity methods of obtaining AOA/AOD/TOA make mmWave an ideally suited technology for localization.
From the beam training stage, we assume that the AOA/AOD/TOA are measured for paths corresponding with beams from N FE FEs. Then, each path is identified as LOS or NLOS, which is discussed in 2.3. It is assumed there are N L ≤ N FE LOS paths and N N NLOS paths, which are used to estimate the location of a UE at position ( 1 ) Figure 2 shows NLOS and LOS paths along with each path's associated parameters, which are explained in the following paragraphs. Let the vector q ∈ R 2×N FE be the two-dimensional coordinate vector for N FE transmitting FEs and q(k j ) = q x k j q y k j T (2) be the location of FE k j that transmits path j. The vector s ∈ R 2×N N is the two-dimensional coordinate vector for the locations of all the scatterers of the NLOS paths, and is the location of scatterer i j that reflects path j.
Let α(θ) ∈ R (N L +N N )×1 contain the AOA for all LOS and NLOS paths, β(θ ) ∈ R (N L +N N )×1 contain the AOD for all LOS and NLOS paths, and d(θ ) ∈ R (N L +N N )×1 contain the total traveled path distance (which TOA will measure) for all LOS and NLOS paths. A NLOS path j transmitted from FE k j and reflected from scatterer i j will have AOA α j,k j ,i j , AOD β j,k j ,i j , and distance d j,k j ,i j . A LOS path j transmitted from FE k j will have AOA α j,k j ,0 , AOD β j,k j ,0 , and distance d j,k j ,0 where the index 0 signifies a LOS path.
For consistency, let the first N L elements of α(θ), β(θ), and d(θ ) correspond to LOS paths so that these elements have indices j = 1, ..., N L ; k j = 1, ..., N L ; and i j = 0. The remaining N N elements are the remaining NLOS paths with indices j = N L + 1, ..., N L + N N ; k j corresponds to the FE locations from which NLOS path j originated; and i j = 1, ..., N N corresponds to the scatterer location from which path j reflects. Then, similar to [8], the measured/observed parameters for LOS paths to be used for localization can be written as On the other hand, NLOS paths have the parameters where atan2(·) is the four quadrant inverse tangent: atan2 and atan(·) is the ordinary arctangent function.
Our goal is to estimate p = p x p y T from noisy measurements of the parameters in (4)-(9) from multiple links with one or more FEs. The noise of the AOA/AOD/TOA path parameters is assumed to be zero-mean Gaussian throughout this paper, and the estimation algorithm is designed under the Gaussian noise assumption. Though, the Gaussian distribution is just one possibility for parameter noise distributions. The proposed algorithms can still be used if the parameters have non-Gaussian noise distributions, but with degraded performance. The localization estimator is optimal under the Gaussian noise scenario, but scenarios with non-Gaussian distributions are still expected to have reasonable performance [26]. Nonetheless, the Gaussian assumption is reasonable for AOA/AOD/TOA at mmWave frequencies as supported by [27][28][29], where [29] uses a measurement campaign to show that Gaussian noise is a good fit for AOD and AOA at mmWave frequencies in urban environments. TOA estimation can be used to determine the total distance in (6) and (9) traveled by a path for LOS and NLOS links, respectively. The measured distance for a path is with n where l = L for LOS and l = N for NLOS paths. The center direction of the transmitted beam can be transmitted by the FE and used as an AOD estimate. The measured AOD is with n (β) j,k j ,i j ∼ N 0, σ 2 β l where l = L for LOS and l = N for NLOS paths. The AOA can be measured by a receiver array so that the measured AOA is with n (α) j,k j ,i j ∼ N 0, σ 2 α l where l = L for LOS and l = N for NLOS paths.
Without loss of generality, the variance of the parameters have been assumed to only depend on whether the path is LOS or NLOS. The parameters for the paths in (4)-(9) depend on the unknown position of the UE p and the unknown scattering locations s, which are nuisance parameters. The vector of unknown parameters that must be estimated is where ∈ R 2+2N N is the unknown parameter space. We also consider networks that have REM capabilities. REM provides information about the environment such as building locations and path link information, which can be used to determine scatterer locations and improve localization performance. The best case scenario is when the REM can provide sufficient information so that the scatterer locations are known perfectly. Then, the estimator no longer has to estimate the scatterer locations for each NLOS path and only the UE position p needs to be estimated. Therefore, the performance of REM aided localization is bounded by the case where scatterer locations are known and constant so that the UE position coordinates are the only parameters that need to be estimated and the unknown parameter vector becomes: where ∈ R 2 . With or without REM, the mmWave parameter measurements in (4)-(9) can be characterized using a nonlinear Gaussian model, which fits the general form: where the nonlinear function h : R N → R M and observations z are given by: Localization then requires estimating the UE position p from the measurements z.

Statistics of AOA, AOD, and TOA
The authors in [29] use measurements and simulations in urban environments to characterize the statistics of AOA, AOD, and TOA (needed for (17)) for a system with 2.5 ns multipath resolution and 800 MHz bandwidth, at center frequencies 28 and 73 GHz. Horn antennas are used with 10.9°/8.9°azimuthal/elevation 3 dB beamwidths at 28 GHz and 7°/7°azimuthal/elevation 3 dB beamwidths at 73 GHz. It has been shown that the parameters given closely capture the statistics of measurements in urban environments at these frequencies. We use these parameters for each receiver in simulations, which are summarized in Table 1. Using the model from [29] provides an accurate representation of an urban environment so that ray tracing is not necessary, which greatly simplifies simulations. It should be noted that the source of noise from the experimental campaign in [29], which leads to the variances in Table 1 is from multipath in the same timing bin and the inability to separate them. Antennas with more directionality will reduce the number of multipaths and reduce the noise variances from Table 1 and improve localization accuracy. The included AOD and AOA noise variances in Table 1 that are extracted from mmWave measurements in [29] are assumed to capture the effects such as NLOS reflection power loss, beamwidth, and interfering multipath components. It should also be noted that these parameters assume the same receive SINR for both central frequencies of 28 and 73 GHz as there are no interfering transmitters using the same bands; in other words, impact of central frequency is limited to the statistics reported on Table 1, and interfering path behavior at these two different mmWave bands are not explicitly taken into account. On the other hand, mmWave frequencies have little interference due to highly directional communications and highly attenuated reflections [1].

Differentiating LOS and NLOS paths
A challenging aspect for mmWave localization will be LOS/NLOS path identification and path separation, which is needed for generating the vector of unknown parameters in (14) (includes one unknown scatterer location for each NLOS path, and no scatterer for a LOS path), the covariance matrix in (17) and hence the likelihood function to be defined in (20). To this end, beam directionality can play an important role in LOS/NLOS differentiation and path separation since it significantly impacts the multipath characteristics in mmWave systems. To gain further insights, an example with Remcom's Wireless Insite mmWave ray-tracing simulator is used to simulate path separation with different directionalities of beams. A center frequency of 73 GHz is used, and a geometry is chosen that has multipath from ground bounce and building reflection as seen in Fig. 3a. The FE is at a height of 10 m, and the UE is at a height of 2 m. The geometry is held fixed, and two different beamwidths are used. The first beam is semi-directional with a 80°3 dB beamwidth, and the second is a directional beam with a 28°3 dB Table 1 Standard deviation of the observation parameters that are used in the simulations beamwidth. Figure 3b shows the TOA for the 80°3 dB beamwidth, where two clusters are observed, one of which has the LOS path along with multiple NLOS paths including ground bounce and reflected building paths. In this case, the LOS path is inseparable from others as mmWave systems will have a typical multipath resolution of 2.5 ns. This leads to a noisy AOA measurement and reduced accuracy. On the other hand, Fig. 3c shows the TOA for the beam with the 28°3 dB beamwidth. In this case, the directional beam gain causes all multipath including the ground bounce path to be negligible and only the LOS path is detected. Since the 3D analysis of directional beams are able to separate out ground bounce and multipath from the main path, for simplicity, it is reasonable to use a 2D model as path separation and localization results will be similar to a 3D model. Directional beams and unique properties of mmWave such as signal attenuation, small RMS delay spread, and sparsity make differentiating individual paths and separating LOS path from NLOS paths feasible [10,24]. The work in [10] uses sparse estimation to separate and identify paths, but with added computational complexity. Simpler methods also exist that can separate LOS from NLOS paths. For example, the AOD and AOA can be compared using the overall coordinate system, which can be computed from the AOD and AOA relative to the antenna arrays if the FE and UE orientation are known. Then, a LOS path will have the parameters (4)-(9) with an AOD β and AOA α such that |α − β| = π. However, noisy versions of the parameters are observed with measured AOD β and measured AOA α . Therefore, a threshold ξ can be introduced that identifies the paths as LOS or NLOS for a desired probability of error: where I LOS (α , β ) is an indicator function, which is 1 if LOS, and 0 otherwise. Another approach to separate LOS and NLOS paths is to exploit TOA and RSS. Urban scenario path loss for 5G systems is studied in [30] and [2] where it is shown that LOS and NLOS paths have different RSS statistics and a reflection from a NLOS path will result in 4-6 dB of power loss. Therefore, the path distance from a TOA estimate can be compared to the theoretical RSS path loss in free space to differentiate between LOS and NLOS paths.

mmWave location estimation
In this section, we will first describe the likelihood function for a UE's location based on AOA, AOD, and TOA observations. Subsequently, we will study gradient methods, particle filter methods, and a gradient-assisted particle filter technique to solve the two-step mmWave localization problem.

Maximizing the likelihood function
It is assumed that the distributions of the parameter observations (11)-(13) are known and used to calculate the covariance R as in (17). Then, an estimator for an unknown UE location can be obtained by maximizing the likelihood of a measurement (16) as follows: where the likelihood for a measurement is and enforces the difference between angular parameters (AOD and AOA) to be in the range [ −π, π ]. This is achieved with the following modulus function that forces its argument into the interval [ −π, π]: This then prevents a linear treatment of angles, which have a circular coordinate system and two angles cannot have a magnitude difference greater than π . Without the modulus, the likelihood function will have discontinuities that result from the atan2 function in (4), (5), (7), and (8), which results in estimators with inherent bias [31].

Non-REM-assisted localization
Networks without REM must estimate the UE coordinates p jointly with the scatterer locations (nuisance parameter) s. Thus, maximization of the likelihood function as in (19) is required over the 2 + 2N N dimensions of θ = p T s T T to localize a UE. The nonlinearity of (16) requires a nonlinear estimation algorithm to maximize (20). Furthermore, localization without REM is particularly challenging as it leads to a non-convex likelihood function that must be maximized over many dimensions. The global maximum of the likelihood function is the optimal estimate.

REM-assisted localization
The REM provides a map of the estimated building and scatterer locations, which can be used to estimate a scatterer location for an NLOS path. It should be noted that a realistic implementation of the REM requires addressing many challenging aspects. A major challenge in a realistic REM implementation is linking a measured NLOS path's AOA/AOD/TOA parameters to a particular scatterer location. One method of achieving REM is to use AOA/AOD/TOA to compute and store the scatterer location for each NLOS path received by a node in a communicating network. This allows a database to be created that serves as an REM map, which stores scatterer locations as well as AOA/AOD/TOA and other path information for each observed path. Then, the scatterer location from an unknown NLOS path can be linked with an observed AOA/AOD/TOA by searching the REM map database for scatterers that are associated with similar path parameters. The authors in [19] and [20] provide algorithms that use REM to reduce the search space for UE path tracking. Further details on REM implementation is left as future work.
Instead of focusing on REM implementation, we use a perfect REM system with zero scatterer location estimation error to provide bounds on the limits of REM performance. This serves as a performance bound to any realistic REM system that can be implemented. Thus, the performance of a real REM system will be in between the perfect REM bound and the bound of a system without REM. A perfect REM system has knowledge of scatterer locations s so that s can be treated as a constant and the dimensionality of the likelihood function is reduced to two dimensions so that a UE is localized by maximizing the likelihood function over θ =[ p].

Gradient methods
Gradient optimization methods such as gradient descent, Newton's method, Gauss-Newton method, and Levenberg-Marquardt method are often used for maximizing likelihood functions. These algorithms all rely on the first and second derivatives of the nonlinear function to iteratively find the maximum of (20). However, these optimization methods alone are not guaranteed to find the global maximum or even converge. If (20) has many local maxima, these methods will converge to the local maximum that is closest to the initial estimate. An advantage of these methods is that they are computationally efficient. The function lsqnonlin in MATLAB implements the "trust-region-reflective" algorithm, which is a variant of Newton's method. It will be referred to as "gradient method" for the rest of this paper as any of the other mentioned gradient-based algorithms will give similar performance.

Particle filters
As an alternative to gradient methods, an exhaustive search estimator is guaranteed to find the global maximum, but is not computationally feasible for maximization over high dimensional functions. Particle filters can be used for multivariate nonlinear estimation and attempt to approach the performance of exhaustive searches with less computational effort [32]. Each particle corresponds to a point where the likelihood function is evaluated. Instead of computing the cost function over the entire space, particle filters use random particle searches to focus on likelihood maxima. However, the number of particles required to ensure that the global maximum is reached increases with search space dimensionality. The sequential importance resampling (SIR) particle filter [33] is a simple technique compared to other types of particle filters and is modified in our implementation.
SIR particle filters are typically used to track a target state vector θ k ∈ R N in a nonlinear time-dependent system where t = t k is the time instant at sample index k ∈ N. The target state evolves in time by where f k−1 (·) is a nonlinear function from the target model and v k−1 is process noise, which is included for errors in the model. A measurement of the target state at time t = t k is given by: where h k is a nonlinear measurement function and w k is measurement noise. Typically, a particle filter is used to estimate the target state θ k as it evolves in time (as opposed to h(θ ) in (16)) from the mesurements z k . However, it can also be used as a maximum-likelihood estimator. Instead of letting k represent a state in time, let it represent an iteration and the target not evolve; in other words, let the observations be given by: Then, instead of estimating an evolving target state at each instant of time with an additional measurement z k , the target state estimate is improved with each iteration k from a single measurement from (16). The target state θ k then represents a collection of N particles, which are evaluations of the likelihood function where θ i k is used to represent an individual particle. Under these modifications, the measurement and process noise from (24) and (25) Since each iteration uses the same measurement (z k = z), it is seen that (26) is identical to (20) if the covariance R k = R and state θ k = θ. The proposed GAPF estimation algorithm is a joint particle filter gradient method algorithm that is able to overcome the weaknesses of the individual gradient-based and particle filter-based algorithms. The gradient method is used to find likelihood function maxima nearest to particles, which reduces the number of particles required to search over the parameter space. The particle filter enables random searches of the parameter space in search of other maxima, which aims to eliminate convergence to local maxima.

Gradient-assisted particle filter
Each iteration of the GAPF estimator in algorithmic form [33] is given in Algorithm 1. First, the particles are resampled using the particles θ i k−1 and weights w i k−1 from the last iteration. The resampling function is taken from [33], which draws particles in proportion to the weights so that particles that correspond to small likelihood evaluations are replaced by duplicate particles with larger likelihood function evaluations. Then, the particles are drawn according to θ i k ∼ p θ k |θ i k−1 from (27), which spreads out the particles randomly with the process noise covariance Q k−1 defined by the user. Subsequently, the gradient method is applied to these drawn particles to find the nearest local maxima. The likelihood of each particle is evaluated according to (26) to determine the weight for that iteration. The particle with the greatest weight is the target parameter estimate. An iteration of the particle filter estimator is visualized in Fig. 4. Each iteration of the particle filter begins on the peaks of a set of local maxima. The particles are randomly spread out to surrounding maxima, and the gradient method finds the peaks of the surrounding maxima. Typically, mmWave localization leads to a likelihood function with a cluster of peaks. The estimator must iterate until all of the peaks have been located to find the global maximum.

Initialization
An initial grid search is required for estimation where the grid points/particles for the initial search θ i 0 are fed into the first iteration of the estimator from Algorithm 2. The grid must be extensive enough to guarantee that at least one of the initial particles in the grid is close to the likelihood function maximum, but not so extensive that it adds unnecessary computation. Sufficient conditions for localization require a minimum of one LOS path or two NLOS paths. It is noted that [10] requires three NLOS paths. The discretion occurs because we assume the UE orientation is known and [10] additionally estimates orientation.
Each parameter being estimated adds another dimension that the initial grid must cover. Scenarios with only LOS paths only require a 2D grid to estimate the UE coordinates. On the other hand, NLOS scenarios require the 2D UE coordinates to be estimated in addition to scatterer coordinates for each NLOS path. This adds an extra two dimensions that must be searched over for each NLOS path, making it difficult to cover the entire search space. One option for scenarios with many NLOS paths is to jointly use all paths, which becomes computationally expensive as it requires many particles to cover the high dimensionality that results from all of the unknown scatterer locations. A second option takes advantage of the fact that only two NLOS paths are required to localize the UE and estimate both scatterers. Then, two or three NLOS paths can be separated from the other paths and used to estimate the UE location and scatterer locations for those paths. This can be done for different sets of paths until each path has been used, which results in scatterer location estimates for all paths with significantly less numbers of required overall particles. Then, the UE location estimates from every set of paths can be averaged to give a better estimate for the UE location.
Instead of a full exhaustive grid search in each dimension for the initial particles θ i 0 , an alternative initialization is proposed as seen in Algorithm 2 that performs a separate grid search that significantly reduces the number of grid points required for scenarios with many NLOS paths. Algorithm 2 reduces the dimensionality of the search space by selecting paths such that a grid search is performed with a series of 2D grids rather than the entire grid space at once. For example, a scenario with two LOS paths and two NLOS paths is initialized by first only using the two LOS paths to search over the 2D UE coordinates for an initial estimate. Then, the estimated UE coordinates are held fixed and each of the NLOS paths are analyzed individually with a 2D grid search for the scatterer location. It should be noted that NLOS scenarios require the use of at least two paths in order to have enough measurements to calculate the UE coordinates, which also require estimating the scatterer locations. A single grid point is generated from initialization, which is fed into the particle filter as θ 0 . The initial grid point gives a rough estimate of the peak, and the GAPF estimator searches around this point to find a better estimate.

Algorithm 2 Grid initialization: N L LOS paths and N N NLOS paths
if N L > 0 then • Search over UE coordinates (p x , p y ) using N L LOS paths as measurements.
• Let the coordinates that maximize the likelihood be (p x ,p y ). if N N > 0 then for j=1:N N do • Hold UE coordinates fixed at (p x ,p y ) and search over scatterer j coordinates (s x (j), s y (j)) using all LOS paths and NLOS path j as measurements.
• Let the coordinates for scatterer j that maximize the likelihood be (ŝ x (j),ŝ y (j)). end for end if else • Search over UE coordinates (p x , p y ) and scatterer coordinates (s x (m), s y (m), s x (n), s y (n)) using NLOS paths m and n as measurements.
• Let the grid point that maximizes the likelihood have UE coordinates (p x ,p y ) and scatterer coordinatesŝ = (ŝ x (m),ŝ y (m),ŝ x (n),ŝ y (n)). for j=1:N N (j = m, n) do • Hold UE coordinates fixed at (p x ,p y ) with other scatters fixed atŝ and search over scatterer j coordinates (s x (j), s y (j)) using NLOS path j as measurements.
• Let the coordinates for scatterer j that maximize the likelihood be (ŝ x (j),ŝ y (j)) and add toŝ.

Fundamental lower bounds for mmWave localization
Lower bounds on mmWave localization performance can provide insight into the limits of 5G networks and are helpful in identifying key factors when designing a network to meet certain specifications. The CRB is often used to bound localization performance as it provides a lower bound on the covariance Cθ of any unbiased estimator θ that satisfies the regularity conditions in [26]: where I(θ) is the Fisher information matrix and ≥ 0 represents a positive semidefinite matrix. The covariance on estimatorθ is defined as, Additionally, using (28), it can be shown that the variance for element m ofθ is bounded by, for all m. Since the observations are assumed Gaussian z ∼ N (h(θ), R) and the covariance does not depend on θ , each element of Fisher information matrix is given by [26]: where 1 ≤ m, n ≤ M and ∂h(θ) ∂θ m is the mth column of the Jacobian: It should be noted that the estimator only approaches the performance bound if the observation noise is actually Gaussian. However, the CRB still bounds non-Gaussian scenarios as other noise distributions will lead to degraded performance [26].
For completeness, the individual elements of the Jacobian matrix are calculated in "Appendix: Elements of the Jacobian matrix. " For an unbiased estimator, we need to have E[θ ] = θ and the mean square error (MSE) for the mth element of θ is equivalent to the mth diagonal element of the covariance, given by Using (30), the MSE for the mth element of the unbiased estimator is bounded by the mth element of the inverse of Fisher information matrix, which is substituted from (31): An estimator with MSE that achieves the equality in (33) is referred as an efficient estimator.
In the rest of this paper, the root-MSE (RMSE) will be used instead of the MSE where RMSE θ m = MSE θ m so that the RMSE of a UE position estimator is in meters (m). From (33) the RMSE of the twodimensional estimate of the UE coordinates (RMSE est (θ )) are lower-bounded by RMSE CRB (θ), where is the RMSE for the unbiased estimatorθ and is the CRB. The UE coordinate estimator RMSE in (34) cannot be evaluated exactly, but it can be approximated with Monte-Carlo (MC) simulation, where N sim is the number of Monte-Carlo simulations and p x i andp y i are the ith Monte-Carlo simulation estimates of the UE position coordinates p x and p y . Another bound that was considered was the periodic CRB (PCRB) [34]. The likelihood function (20) prevents noisy angular measurements differing from the actual angle by more than π. This is a cyclic Gaussian process where very noisy measurements can wrap around and be closer to the actual estimate, which needs to be taken into account in bounding the parameters. The PCRB bounds such processes and was implemented. However, it was found that the angular noise required for the PCRB to take effect was larger than what we consider.

Numerical results and discussion
In this section, we evaluate the localization performance of mmWave networks in urban environments. First, Monte-Carlo simulations of the GAPF estimator and the CRB are used to analyze localization performance as a function of beamwidth. Then, mmWave localization performance is studied in urban canyon and urban corner scenarios with one and two FEs using both LOS and NLOS observations, where the urban canyon is represented by two parallel walls and an urban corner is represented by two intersecting orthogonal walls.

Localization performance as a function of beamwidth
The trade-off between beamwidth and localization performance is examined in an urban corner scenario with one LOS and two NLOS paths available as seen in Fig. 5a. The urban corner is simulated with a vertical building wall at x = 0 (m), a horizontal building wall at y = 0 (m), an FE at (18, 10) (m), and a UE at (8, 35) (m). A Monte-Carlo simulation is run where the TOA distance noise standard deviation is fixed at σ d L = σ d N = 0.75 (m). The AOD and AOA noise standard deviations are assumed equal to σ angle = σ α L = σ β L = σ α N = σ β N and increased. Equating all of the angular noise standard deviations and varying σ angle is similar to what occurs if a mmWave antenna varies its beamwidth [35]. A larger beamwidth results in a larger AOD variance since the it is more likely the path did not travel from the center of the beam, which is what is transmitted as the AOD estimate. Larger beamwidth also results in more multipath energy, which corrupts the AOA estimate at the receiver and increases the AOA variance. Thus, varying σ angle has similar effects as varying beamwidth and is useful in providing insight into beamwidth effects on localization performance. with the GAPF estimator. The RMSE is shown for both REM and non-REM systems utilizing different subsets of the available paths from Fig. 5a. Monte-Carlo simulations show the estimator to be closely aligned with the CRB, providing evidence that the GAPF estimator is an efficient estimator. It should be noted that implementation of the particle filter or the gradient method individually leads to RMSE results that are not near the CRB. This results from the convergence of the estimator to local maxima and leads to RMSE values that are larger than those given by the CRB. The individual methods only approach the CRB under certain conditions if initialization is able to provide initial estimates in the immediate vicinity of the global maximum, which is not normally expected. For systems without REM, it is observed that increased beam directionality leads to linear improvement in RMSE values and localization performance. On the other hand, more directionality comes at a cost because it is achieved by adding antenna elements at the transmitter. It may not be feasible to fit the required number of antenna elements at the transmitter to reach a desired directionality. Furthermore, the cost of the transmitter grows with each additional antenna element. Thus, transmitter size constraints and affordability limit the achievable directionality. Additionally, it is seen that utilizing only the LOS path has much better localization performance than using only the NLOS paths while the addition of NLOS paths to the LOS path only provides modest improvements. The REM system RMSE curves reach a threshold where the RMSE does not increase with larger beamwidths. This occurs because REM provides knowledge of the scatterer locations. At smaller beamwidths, angles provide precise information and reduce the area of likely UE positions. Eventually, as the beamwidths increase, a threshold is reached where AOD/AOA measurements become less relevant and the angular uncertainty allows the spread of likely UE positions that are too disparate from the TOA measured path distances. This threshold has important implications in designing mmWave systems for localization. High accuracy localization can be achieved by two methods. Highly directional beams from antennas with many elements can be implemented, but this can be costly and increases the difficulty of beam alignment. On the other hand, systems can use REM, which has a high computational cost. From this example, it is apparent that environments that are mainly LOS, such as rural areas must rely on increasing beam directionality to improve localization performance. Scenarios with many scatterers and NLOS paths, such as dense urban environments can most easily improve localization performance with the use of REM if the network can handle the computational load. It should be noted that these calculations are under the assumption of perfect REM where the scatterer for each NLOS path is known without any estimation error. The plotted results are a bound on the performance of any REM system. Realistic localization with REM will have nonzero errors in scatterer location errors and has many challenges that must be addressed before obtaining small error values. Thus, an actual REM system will have RMSE values between the non-REM CRB and REM CRB with values depending on the scatterer location estimation error that the REM localization system is able to achieve.

Localization in urban environments
In this section, we consider mmWave localization scenarios in urban canyon and urban corner environments where 5G mmWave cellular networks are expected to first get deployed and there may be one or two FEs available to be used for localization purposes. The performance with/without the use of REM for localization is analyzed, where with REM, we assume the scatterer locations are available as in (15

Urban canyon
An urban canyon is simulated with parallel walls at x = 0 (m) and x = 20 (m) as shown in Fig. 6a with an exam-ple set of paths for an FE on the left wall at (−1, 2) (m) and a UE located at (10, 40) (m). In this case, a single LOS and a single NLOS path reflected by the right wall are received by the UE from the FE. From Fig. 6b, it is seen that localization performance is improved with the addition of REM as expected. As with all scenarios, 73 GHz performs slightly better than 28 GHz. Referring to Table 1, 73 GHz is equivalent to 28 GHz in timing estimation, but has better estimation of AOA and AOD. As it is at a higher frequency, more antennas can be fit on a chip, allowing a more directional beam. The more directional beam captures less multipath and reduces the variance on AOA and AOD estimation, enabling improved performance. Figure 6c- Another scenario adds a second FE to the urban canyon at (21, 48) (m), which has an additional LOS link and a NLOS link reflected from the left wall as seen in Fig. 7a. The RMSE curve seen in Fig. 7b shows that REM again greatly improves localization performance. We also realize from Fig. 7b that adding a second mmWave BS significantly improves the localization accuracy with and without REM. Figures 7c-e show that REM offers more performance gain as the UE travels further from either FE.

Urban corner
An urban corner is simulated with a vertical building wall at x = 0 (m) and a horizontal building wall at y = 0 (m) as seen in Fig. 8a with an example set of paths for a FE at (18, 10) (m) and a UE at (15,15) (m). In this case, two LOS paths and two NLOS paths, one reflected from each wall are received by the UE from the FE. Figure 8b shows that the geometry of the corner allows better localization accuracy relative to the canyon as it allows more NLOS paths. From Figs. 8c-e, it is seen that REM provides the most gain far from the FE, which is most significant close to the left reflected wall.
A second FE can be added to the urban corner scenario at (18, 48) (m) as in Fig. 9a. The geometry of the corner allows four NLOS paths in addition to two LOS paths, providing large amounts of information for localization. It is seen from Fig. 9b that the use of many paths results in high localization accuracy. Figure 9c-e shows that localization performance does not only depend on the distance from the FEs. The midpoint between the two FEs has worse a b c d e Fig. 9 a Urban corner with two FEs, two LOS paths, and two NLOS paths. b CDF of RMSE. c RMSE at 73 GHz without REM. d RMSE at 73 GHz with REM. e Improvement in RMSE with REM at 73 GHz performance than points further from both FEs that are closer to the left wall. This scenario is similar to what is expected for the first deployed 5G networks where there may be multiple FE in range and many LOS/NLOS links may be available and exploited for significant gains in localization performance. These results provide evidence that 5G networks using two-step localization in environments rich in available links can expect sub-meter localization accuracy even without REM.

Building REM from localization outcomes
The scatterer locations are nuisance parameters in the estimation of the UE position. However, the scatterer locations are valuable information and can be used to create an REM. Figure 10 plots the estimated scatterer locations (red) obtained from NLOS paths for a given UE trajectory (green dots). Figure 10a shows the estimated scatterer locations from the urban canyon in the scenario shown in Fig. 7. It is seen that scatterer estimates are along the entire length of both the left and right building as both of the FEs are able to establish paths that reflect from both buildings, providing good coverage of the environment. Figure 10b shows the estimated scatterer locations from the urban corner in the scenario shown in Fig. 9. In this scenario, each FE is able to establish two NLOS paths for each UE location. However, it is seen that the very interior of the corner has no scatterer coverage. Coverage is possible if a UE travels further into the corner, but it should be noted that REM may not always provide complete coverage of the environment. The scatterer locations can be stored with the AOD/AOA/TOA parameters associated with them to create an REM. An algorithm that uses the information from REM to assist in localization is left for future work.

Conclusions
In this paper, we analyze the localization performance of a reduced complexity method for 5G mmWave urban networks with multiple available LOS/NLOS links with one or more FEs. Specifically, we study the localization performance in urban canyon and urban corner settings utilizing AOA, AOD, and TOA measurements at a UE from one or more FEs. We consider scenarios with and without REM, where in the latter case, all the scatterer locations for NLOS links are also simultaneously estimated along with the UE location. This results in a high dimensional unknown parameter space. As a consequence, mmWave localization requires processing of a likelihood function for the unknown parameter vector with many local maxima, making it difficult to find the global optimum solution. To address Fig. 10 Mapped scatterer locations (red) generated from GAPF estimator at UE locations (green) on a path through a the urban canyon and b the urban corner this problem, we propose a GAPF estimator which combines particle filter and gradient method, and Monte-Carlo simulations show that its performance matches that of the CRB in terms of localization RMSE. The scatterer locations that are estimated using the proposed algorithm can be used to create an REM for an urban environment.
Our results show that when REM is used with the proposed two-step approach, sub-meter localization accuracy is feasible in mmWave urban networks using even a single FE. However, further research is required for evaluating the feasibility of realizing REM effectively. Without REM (when scatterer locations need to be estimated simultaneously), median RMSE lower than two meters is possible with a single FE, and lower than a meter is possible with two FEs. Results show that the urban corner provides better localization performance due to availability of larger number of NLOS paths. Thus, dense urban environments with non-trivial building structures and many scatterers are best suited for mmWave localization. It is shown that localization is improved by increasing the number of antenna elements to increase beam directionality or by implementing REM to assist in estimating NLOS scatterer locations, both of which come with the trade-offs of higher cost and complexity.

Appendix: Elements of the Jacobian matrix
The Jacobian in (32) is calculated as: