Millimeter-wave Mobile Sensing and Environment Mapping: Models, Algorithms and Validation
Carlos Baquero Barneto, Elizaveta Rastorgueva-Foi, Musa Furkan Keskin, Taneli Riihonen, Matias Turunen, Jukka Talvitie, Henk Wymeersch, Mikko Valkama
11 Radio-based Sensing and Environment Mappingin Millimeter-Wave 5G and Beyond Networks
Carlos Baquero Barneto ,
Student Member, IEEE , Elizaveta Rastorgueva-Foi , Musa Furkan Keskin ,Taneli Riihonen ,
Member, IEEE , Matias Turunen, Jukka Talvitie ,
Member, IEEE ,Henk Wymeersch ,
Senior Member, IEEE , and Mikko Valkama ,
Senior Member, IEEE
Abstract —Integrating efficient connectivity, positioning andsensing functionalities into 5G New Radio (NR) and beyondmobile cellular systems is one timely research paradigm, espe-cially at mm-wave and sub-THz bands. In this article, we addressthe radio-based sensing and environment mapping prospect withspecific emphasis on the user equipment (UE) side. We firstdescribe an efficient (cid:96) -regularized least-squares (LS) approachto obtain sparse range–angle charts at individual measurementor sensing locations. For the subsequent environment mapping,we then describe both grid-based static solution as well as moreadvanced tracking-based dynamic approaches, where interactionmultiple-model extended Kalman filtering and smoothing areutilized. We provide numerical indoor mapping results at 28 GHzband deploying OFDM-based 5G NR uplink waveform with400 MHz channel bandwidth, covering both accurate ray-tracingbased as well as actual RF measurement results. The resultsillustrate the superiority of the dynamic tracking-based solutions,while overall demonstrate the excellent prospects of radio-basedenvironment sensing and mapping in future mm-wave networks. Index Terms —5G New Radio (NR), 6G, joint communicationsand sensing, RF convergence, indoor mapping, mm-waves, sub-THz, OFDM radar.
I. I
NTRODUCTION F IFTH generation (5G) New Radio (NR) mobile cellularsystems provide large improvements in terms of, e.g, peakdata rates, network capacity, number of connected devices, andradio access latency, compared to earlier Long-Term Evolution(LTE) based networks [1]. While the primary purpose ofmobile networks is the efficient data connectivity, they canalso facilitate terrestrial positioning service – with ambitiousaccuracy requirements of around one meter in 5G NR [2]–[5].Additionally, to further leverage the capabilities of millimeter-wave (mm-wave) transceiver systems, the radio-frequency (RF)convergence paradigm refers to integrating also radio-basedsensing capabilities to the networks [6]–[10]. To this end,the development of joint communication and sensing (JCAS)
This work was partially supported by the Academy of Finland (grants systems that can perform both functionalities while sharingthe same transmit waveforms and hardware platforms is avery timely research area [11]–[14], with applications, e.g., invehicular systems and indoor mapping [11], [15]–[18]. Thelarge channel bandwidths available at mm-wave frequencies,together with highly directional antenna arrays, form thetechnical basis for high-accuracy radio positioning and sensing[4], [19], [20]. While 3GPP 5G NR supports currently channelbandwidths up to 400 MHz [21], further increased channelbandwidths are expected through the NR evolution towards thesub-THz regime [22], paving the way to future 6G networks[16].In general, in typical outdoor mm-wave scenarios withdensely deployed network nodes or radio heads, the userequipment (UE) is likely to have line-of-sight (LoS) connectionto the network while also being within the coverage rangeof multiple network nodes — both being issues that furtherfacilitate high-accuracy positioning [2]–[4]. However, in indoor deployment scenarios, the propagation characteristics and thusthe positioning become much more challenging due to the highpath loss as well as the reduced penetration and diffraction [23],[24]. To address such challenges, the concept of simultaneouslocalization and mapping (SLAM) is being studied, e.g., in[25], [26] where the environment sensing and mapping isone key ingredient. To this end, in [19], a personal mobileradar concept for environment sensing or mapping is studied,building on large antenna arrays at mm-waves. The work in[25], in turn, described a SLAM architecture for mm-wavelocalization through obstacle detection and dimensioning forindoor environments. Furthermore, a message passing-basedestimator which jointly estimates the position and orientation ofthe UE while sensing the environment’s reflectors is presentedin [27]. In [28], a JCAS system is pursued for smart homescenarios, investigating different integration approaches fortraditional sensing and WiFi devices. In [26], a map-free indoorlocalization method using ultra-wideband mm-wave signals isproposed, utilizing the estimated scatterer positions of theenvironment as virtual anchors to obtain the device’s location.Due to the notable challenges of performing RF measurementsat mm-wave bands, only very few works have supported theirpositioning/sensing/mapping solutions with empirical results toassess and demonstrate the performance with real transceiverand antenna hardware [24], [26], [29], [30].In this article, we describe a user-centric mm-wave indoormapping system and associated signal processing methods in5G NR UE context. In the considered approach, illustrated a r X i v : . [ ee ss . SP ] F e b TX RX Radar processing I ndoo r m app i ng p r o c e ss i ng UE (a) (b)Fig. 1. (a) Considered joint communication and sensing/radar scenario in mm-wave 5G NR network, for user terminal based mapping applications. (b) Examplemeasurement based illustration of selected indoor mapping challenges related to angular degradation as well as to multi-bounce phenomenon at 28 GHz. conceptually in Fig. 1, the UE senses the surrounding envi-ronment through beamformed uplink transmissions while thenobserving and collecting the target reflections. This is followedby range–angle processing, while the corresponding range–angle charts are then further post-processed in the subsequentmapping stage. We specifically address the moving nature ofthe sensing and measurement device, in the form of advancedtracking-based dynamic mapping solutions, incorporating bothspecular reflections and diffuse scattering. With the proposedtracking approach, besides mapping of fixed structures, it ispossible to track moving targets, adapt to the changes in theenvironment, as well as predict the mobility of a tracked target.Compared to [19], [31] and to the overall state-of-the-artliterature in the field, the contributions and novelty of the articlecan be summarized as follows: • An efficient (cid:96) -regularized least-squares (LS) approach toobtain sparse range–angle charts at individual measure-ment or sensing locations is described; • Advanced tracking-based dynamic mapping approachesare devised where interacting multiple model (IMM)extended Kalman filtering and smoothing are deployed,allowing efficient tracking of individual scatterers overtime; • Multi-bounce or double-reflection challenge of practicalcomplex environments is addressed and a scatterer cross-identification method is devised; • Numerical indoor mapping results at 28 GHz band deploy-ing orthogonal frequency-division multiplexing (OFDM)based 5G NR uplink waveform with 400 MHz channelbandwidth are provided, covering both accurate ray-tracingbased as well as actual RF measurement results.The rest of this article is organized as follows. The fun-damental system model for collecting beamformed sensingmeasurements and the associated measurement system geom-etry are described in Section II. The LS-based range–angleprocessing methods are addressed in Section III, while SectionsIV and V present the considered static and dynamic mappingsolutions, respectively. In Section VI, the evaluation scenario,ray-tracing environment and the experimental RF measurementenvironment and equipment are described. Obtained numericalresults and provided and analyzed in Section VII, while conclusions are provided in Section VIII. Finally, further detailsof the LS-based range–angle processing solution are describedin supplementary material.II. S
YSTEM M ODEL
A. Basics
The sensing and mapping functionality of the NR UE buildson the OFDM-based uplink transmit waveform whose complexI/Q samples are fully known in the device. Furthermore,the mm-wave UE is assumed to be equipped with phasedarray(s) for both transmit (TX) and receive (RX) functions,thus allowing for beamformed or directional transmission andcorresponding measurements of the target reflections in orderto map the environment. This is conceptually illustrated inFig. 1(a). In this work, we specifically focus on sensing theazimuth domain to create a 2D map of the environment,however, the same principles can be applied for sensing in theelevation domain and subsequently for extracting 3D maps.For given sensing direction, we assume that the UE transmitsa sequence of M beamformed OFDM symbols on N activesubcarriers, with x n,m denoting the transmitted frequency-domain symbol at n th subcarrier in m th OFDM symbol. Tothis end, the frequency-domain OFDM-based radar model with K targets or reflections reads [13], [32], [33] y n,m = w H RX A RX ,n Γ n A H TX ,n w TX x n,m + v n,m , (1)where v n,m , w TX ∈ C N TX × and w RX ∈ C N RX × are thecomplex Gaussian noise and the TX and RX array beamformingweights, respectively. Moreover, N TX and N RX denote thenumber of antenna elements in the TX and RX array, and A TX ,n ∈ C N TX × K and A RX ,n ∈ C N RX × K are the steeringmatrices for the TX and RX array, respectively. The multipathcoefficients of a sensed environment with K target reflectionsare modelled through Γ n = diag ( γ ,n , . . . , γ K − ,n ) ∈ C K × K that is a diagonal matrix with its k th element being of the form γ k,n = b k,n e − j πn ∆ fτ k with | b k,n | = λ n σ k (4 π ) d k . (2)In above, τ k , d k , b k,n and σ k model the k th target’s propa-gation delay, propagation distance, effective attenuation andradar cross section (RCS), respectively, stemming from the well-known radar range equation [32]. Here, λ n refers tothe wavelength of the n th subcarrier while ∆ f denotes thesubcarrier spacing of the OFDM waveform.In this work, for notational simplicity, we consider theadoption of uniform linear array (ULA) for TX and RXbeamforming. To this end, the TX array steering vector for a k th target with azimuth angle ϕ TX k can be expressed at subcarrier n as a TX ,n ( ϕ TX k ) = (cid:104) , e j Υ n ( ϕ TX k ) , . . . , e j ( N TX − n ( ϕ TX k ) (cid:105) T , (3)where Υ n ( ϕ TX k ) = 2 π υλ n sin ( ϕ TX k ) is the correspondingelectrical angle of departure (AoD) for the k th target, wherein υ is the antenna element separation. The RX steering vec-tor a RX ,n ( ϕ RX k ) and the corresponding electrical angle ofarrival (AoA) Υ n ( ϕ RX k ) are obtained similarly. As a result,the set of steering vectors for the considered K targetreflections can be expressed with compact matrix notationas A TX ,n = [ a TX ,n ( ϕ RX ) , . . . , a TX ,n ( ϕ RX K − )] and A RX ,n =[ a RX ,n ( ϕ RX ) , . . . , a RX ,n ( ϕ RX K − )] which are deployed in (1). B. Geometry
Next, we consider the specific subtleties or character-istics related to the fundamental system geometry shownin Fig. 2, where the sensing device is located at position p UE l = ( x UE l , y UE l ) with l denoting the sensing instance orlocation index. Like illustrated in the figure, we allow forseparate TX and RX arrays in our system model, to relax theself-interference challenge in OFDM radars [13], [34]. There-fore, with separate TX array position p TX l = ( x TX l , y TX l ) andRX array position p RX l = ( x RX l , y RX l ) , the physical propagationdelay of the k th target reflection τ l,k = ( d TX l,k + d RX l,k ) /c isdetermined cumulatively according to the distance between theTX array and target position d TX l,k , and the distance between theRX array and target position d RX l,k . Consequently, the distance(i.e., range) between the UE position p UE l and the k th targetposition p SC l,k with azimuth angle ϕ UE l,k can be expressed as d UE l,k = ρ ( τ l,k , ϕ UE l,k ) = (cid:112) ( c τ l,k ) − d ant (cid:114) − (cid:16) d ant c τ l,k cos( ϕ UE l,k ) (cid:17) , (4)where c is the speed of light and d ant corresponds to theseparation between TX and RX arrays. For the special caseof co-located or joint TX-RX antenna system, d ant = 0 whilethe above expression is written in a general form. Furthermore,for a given physical propagation delay τ l,k and the systemgeometry shown in Fig. 2, the transmit angle ϕ TX l,k — andsimilarly the receive angle — can be described as ϕ TX l,k = atan2( y UE l − y TX l + ρ ( τ l,k , ϕ UE l,k ) sin( ϕ UE l,k ) ,x UE l − x TX l + ρ ( τ l,k , ϕ UE l,k ) cos( ϕ UE l,k )) , (5)where ρ ( τ l,k , ϕ UE l,k ) is defined in (4) and atan2( y, x ) is a four-quadrant inverse tangent function.Assuming then that the distances of the considered scattererssatisfy the far-field condition d UE l,k (cid:29) ( N TX υ ) λ n and d UE l,k (cid:29) ( N RX υ ) λ n , where N TX υ and N RX υ denote the physical length xy Scatterer
Fig. 2. Geometry of the considered sensing setup, where the illuminationand measurement system located at p UE l = ( x UE l , y UE l ) with orientation ϕ OR l senses the environment, containing a target at distance d UE l,k with respect tothe UE location. of the arrays, we can express the TX and RX gain patternsassociated to a target or scatterer k as g TX ,n ( ϕ TX l,k ) = a H TX ,n ( ϕ TX l,k ) w TX ,g RX ,n ( ϕ RX l,k ) = w H RX a RX ,n ( ϕ RX l,k ) . (6)For simplicity, in the rest of the paper, we assume frequency-flatgain patterns with g TX ,n ( ϕ TX l,k ) = g TX ( ϕ TX l,k ) and g RX ,n ( ϕ RX l,k ) = g RX ( ϕ RX l,k ) . Finally, the combined TX-RX antenna pattern isdefined as g ( ϕ UE l,k ) = g RX ( ϕ RX l,k ) g TX ( ϕ TX l,k ) = w H RX a RX ( ϕ RX l,k ) a H TX ( ϕ TX l,k ) w TX . (7)For generality, we also note that in the JCAS literature,different approaches have been considered to optimize the TXand RX beamforming weights, w TX and w RX , see, e.g., [35]–[38]. Specifically, different advanced techniques facilitating,e.g., multiple simultaneous beams have been proposed in theserecent works. In this paper, however, we consider more ordinarysingle-beam approach due to its implementation simplicity,while defer to potential multi-beam considerations to our follow-up future work.III. R ANGE –A NGLE P ROCESSING
In this section, we start by extending the basic signal modelin (1) to the case of multiple sensing directions. Then, wedefine the range–angle processing problem for the consideredenvironment mapping setup and state the corresponding LSsolution. Finally, we propose an (cid:96) -regularized version of theLS estimation problem to obtain sparse range–angle charts toefficiently facilitate subsequent mapping phases. A. Problem Statement
First, we extend the signal model in (1) to include multiplesensing directions. For the i th sensing direction with azimuthangle ϕ S i , with i = 0 , . . . , I − , and the m th TX symbol, thefrequency-domain vector x i,m = [ x i,m, , . . . , x i,m,N − ] T ∈ C N × is transmitted. The OFDM frequency-domain measure-ment vector observed and collected by the RX antenna array,containing all the target reflections, can then be expressed as y i,m = C R − (cid:88) p =0 C ϕ − (cid:88) q =0 b p,q g ( ˜ ϕ q − ϕ S i ) c (˜ τ p ) (cid:12) x i,m + v i,m , (8)where (cid:12) denotes the Hadamard (element-wise) product, C R and C ϕ are the numbers of range and azimuth cells in the discretizedrange–angle chart, respectively, and b p,q is the reflection coef-ficient at the ( p, q ) th range–angle cell located at (˜ τ p , ˜ ϕ q ) . Thevariables g ( ϕ ) and v i,m = [ v i,m, , . . . , v i,m,N − ] T ∈ C N × refer to the combined TX-RX antenna pattern, defined in (7),and the additive noise vector, respectively. Furthermore, in (8),the frequency-domain steering vector for a specific cell’s delay ˜ τ p reads c (˜ τ p ) (cid:44) (cid:104) , e − j π ∆ f ˜ τ p , . . . , e − j π ( N − f ˜ τ p (cid:105) T . (9)The purpose of the range–angle processing is to estimatethe reflection coefficients b p,q = ( B ) p,q at predefined range–angle grid locations (˜ τ p , ˜ ϕ q ) with p = 0 , . . . , C R − and q =0 , . . . , C ϕ − , while exploiting the knowledge of the combinedantenna pattern. Note that b p,q = 0 if there is no scatterer atthe ( p, q ) th range-azimuth cell. For notational convenience, wecan next rewrite the frequency-domain measurement vector (8)as y i,m = C R − (cid:88) p =0 b Tp g ( ϕ S i ) c (˜ τ p ) (cid:12) x i,m + v i,m = x i,m (cid:12) (cid:32) C R − (cid:88) p =0 c (˜ τ p ) b Tp (cid:33) g ( ϕ S i ) + v i,m = x i,m (cid:12) CBg ( ϕ S i ) + v i,m , (10)where b p (cid:44) (cid:2) b p, , . . . , b p,C ϕ − (cid:3) T , (11a) g ( ϕ S i ) (cid:44) (cid:2) g ( ˜ ϕ − ϕ S i ) , . . . , g ( ˜ ϕ C ϕ − − ϕ S i ) (cid:3) T , (11b) C (cid:44) [ c (˜ τ ) , . . . , c (˜ τ C R − )] ∈ C N × C R , (11c) B (cid:44) [ b , . . . , b C R − ] T ∈ C C R × C ϕ . (11d)Specifically, the matrix B defined in (11d) represents the range–angle chart to be estimated. Aggregating (10) over all I sensingdirections, we have the following observation matrix for the m th OFDM symbol at a given UE location: Y m = X m (cid:12) CBG + V m , (12)for m = 0 , . . . , M − , where Y m (cid:44) [ y ,m , . . . , y I − ,m ] ∈ C N × I , (13a) X m (cid:44) [ x ,m , . . . , x I − ,m ] ∈ C N × I , (13b) G (cid:44) (cid:2) g ( ϕ S ) , . . . , g ( ϕ S I − ) (cid:3) ∈ C C ϕ × I , (13c) V m (cid:44) [ v ,m , . . . , v I − ,m ] ∈ C N × I . (13d)At any given l th UE location, the range–angle estimationproblem can then be formally defined as follows. Given thetransmit symbols { X m } M − m =0 , the steering matrix C (known Algorithm 1
Iterative Shrinkage/Thresholding Algorithm(ISTA) for Sparse Range–Angle Processing in (25)
Input:
Observation y in (20), step size η , regularizationparameter λ . Output:
Range–angle chart (cid:98) B ISTA . Initialization:
Set k = 0 , B (0) = (cid:98) B LS and η = βλ max ( C H C ) λ max ( GG H ) for some β ∈ (0 , . Repeat B ( k +1) = T λη (cid:18) B ( k ) − η C H (cid:20) CB ( k ) G (14) − M M − (cid:88) m =0 ( X ∗ m (cid:12) Y m ) (cid:21) G H (cid:19) ,k = k + 1 , where T α ( B ) r,q = ( | B r,q | − α ) + sgn( B r,q ) . (15) until convergence through the delay grid locations { ˜ τ p } C R − p =0 ), the antenna patternmatrix G (known through the angular grid locations { ˜ ϕ q } C ϕ − q =0 and the sensing directions { ϕ S i } I − i =0 ), estimate the range–anglechart B from the observation { Y m } M − m =0 in (12). B. LS Solution
The LS approach to estimate the range–angle chart B in(12) is formulated as min B M − (cid:88) m =0 (cid:13)(cid:13) Y m − X m (cid:12) CBG (cid:13)(cid:13) F . (16)After vectorization, we have min b M − (cid:88) m =0 (cid:13)(cid:13) y m − x m (cid:12) ( G T ⊗ C ) b (cid:13)(cid:13) , (17)where x m (cid:44) vec ( X m ) , y m (cid:44) vec ( Y m ) and b (cid:44) vec ( B ) .Further simplifying (17) yields min b M − (cid:88) m =0 (cid:13)(cid:13) y m − Φ m b (cid:13)(cid:13) , (18)where Φ m (cid:44) diag ( x m ) (cid:0) G T ⊗ C (cid:1) ∈ C NI × C R C ϕ , (19)is a known matrix. Defining then y = (cid:2) y T , . . . , y TM − (cid:3) T ∈ C NIM × , (20) Φ = (cid:2) Φ T , . . . , Φ TM − (cid:3) T ∈ C NIM × C R C ϕ , (21)(18) becomes min b (cid:13)(cid:13) y − Φb (cid:13)(cid:13) . (22)The LS estimate of the range–angle chart in (22) can thus beobtained as (cid:98) b LS = Φ † y , (23) where ( · ) † denotes matrix pseudo-inverse. The LS solution canbe approximated as (see the Supplementary Material for thederivations): (cid:98) B LS = (cid:0) C H C (cid:1) − C H (cid:124) (cid:123)(cid:122) (cid:125) Range compression M M − (cid:88) m =0 ( X ∗ m (cid:12) Y m ) (cid:124) (cid:123)(cid:122) (cid:125) Coherent integration G H (cid:0) GG H (cid:1) − (cid:124) (cid:123)(cid:122) (cid:125) Pattern correlation , (24)which consists of coherent integration over the M receivedsymbols, range compression/matched filtering and antennapattern correlation operations. C. Regularized LS for Sparse Range–Angle Processing
To harness the sparsity of the mm-wave propagation channels,and thus to obtain sparse range–angle charts, we consider the (cid:96) -regularized version of the LS problem in (22), expressed as min b (cid:13)(cid:13) y − Φb (cid:13)(cid:13) + λ (cid:13)(cid:13) b (cid:13)(cid:13) . (25)To solve the sparse map recovery problem in (25), we resortto the iterative shrinkage/thresholding algorithm (ISTA) [39]–[41], as outlined in Algorithm 1, where λ max ( · ) denotes themaximum eigenvalue of a matrix and ( · ) + yields the positivepart of a real number while sgn( z ) = z/ | z | yields the sign ofa complex number. Finally, the targets’ Cartesian coordinatescan be calculated from the obtained sparse range–angle chartsimilarly as described below in (26).Next, we describe the actual environment mapping solutions,starting with grid-based static mapping in Section IV followedby dynamic tracking-based mapping in Section V.IV. G RID - BASED S TATIC M APPING
As a baseline or reference solution for static environmentswith fixed structures, building on our early work in [19], weconsider a grid-based static mapping approach to create adiscretized representation of the environment by exploitingthe UE sensing capabilities. As described in Sections II andIII, the UE estimates the range and angle information of itssurrounding environment by steering its beam pattern throughdifferent I azimuth angles. For the l th measurement location,the UE provides a range–angle chart, B l , of size C R × C ϕ . Forsimplicity, we assume that the number of azimuth cells is thesame as the number of steering directions C ϕ = I . The actualgrid-based static map is then obtained by collecting the UEmeasurements from different L locations, providing a total of A = L × C R × C ϕ measured points.Considering now a UE location p UE l = (cid:0) x UE l , y UE l (cid:1) withlocation index l ∈ [0 , L − , the targets’ Cartesian coordinatesof the ( p, q ) th range–angle cell located at (˜ τ p , ˜ ϕ q ) can becalculated as x a = x UE l + ρ (˜ τ p , ˜ ϕ q + ϕ OR l ) cos( ˜ ϕ q + ϕ OR l ) ,y a = y UE l + ρ (˜ τ p , ˜ ϕ q + ϕ OR l ) sin( ˜ ϕ q + ϕ OR l ) , (26)where a = { p + qC R + lC R C ϕ , ∀ p, ∀ q, ∀ l } , and ρ ( τ, ϕ ) isdefined in (4). The parameter ϕ OR l denotes the UE orienta-tion angle as illustrated in Fig. 2. Similarly, the reflectioncoefficients for all the UE locations are collected from the estimated range–angle charts as ˜ b a = ( B l ) p,q for ∀ p, ∀ q, ∀ l with a = 0 , . . . , A − .After estimating and collecting all the range–angle chartsfrom all the UE locations, the environment is discretizedcreating a grid of size C y × C x denoted as M α,β where α ∈ { , . . . , C y − } and β ∈ { , . . . , C x − } denote thevertical and horizontal map’s cell indices, respectively. Then,the reflection coefficient of the ( α, β ) th cell is calculated asthe mean of the reflection coefficients of all detected targetswithin that specific cell as M α,β = (cid:80) a ∈M α,β ˜ b a C {M α,β } , (27)where M α,β = (cid:8) a | (cid:0) x a ∈ χ β (cid:1) ∩ ( y a ∈ ψ α ) , ∀ a (cid:9) . The pa-rameters χ β = [ χ β, min , χ β, max ] and ψ α = [ ψ α, min , ψ α, max ] referto the horizontal and vertical ranges of the ( α, β ) th map cell,respectively. The operator C {·} denotes the cardinality of theset and consequently the total number of targets within aspecific cell. Then, the grid-based static map is smoothed byconvolution with a kernel matrix. The smoothed cells can beexpressed as ˜ M α,β = U (cid:88) u = − U V (cid:88) v = − V Ω u,v M α − u,β − v , (28)where Ω refers to the kernel matrix with − U ≤ u ≤ U and − V ≤ v ≤ V . In this work, we assume a square Gaussiankernel matrix with U = V and standard deviation σ for bothmap’s dimensions described as Ω u,v = 12 πσ exp (cid:20) − u + v σ (cid:21) . (29)Finally, each cell is subject to a thresholding stage toemphasize the most relevant targets of the environment,expressed as ˜ M α,β H ≶ H M th , (30)where M th denotes the detection threshold. Specifically, if ˜ M α,β > M th , the detector declares that alternative hypothesis H is fulfilled and subsequently a target is present in the ( α, β ) th cell. In contrast, if ˜ M α,β < M th , the detector declaresthat no target is detected according to the null hypothesis H .This serves as the baseline reference approach, while a moreadvanced tracking-based dynamic mapping solution is pursuednext. V. T RACKING - BASED D YNAMIC M APPING
A. Background
The grid-based mapping approach, discussed in Section IV,provides a straight-forward solution for mapping static environ-ments with fixed structures. However, it lacks e.g. the capabilityof performing real-time tracking with possible moving targets,which inevitably limits the utilization of full potential of radio-based sensing and mapping. In this section, we propose andformulate an IMM filtering framework, which is able to trackand also predict changes in the radio environment, and thusprovides means for both real-time radio-based sensing and environment mapping. Moreover, the considered IMM-basedapproach can continuously adapt the movement model of atracked target and predict the target behavior for upcomingtime steps in the tracking process.In general, the transmitted signal is reflected from thesurrounding scattering objects that range from stationarystructures, such as walls and furniture, to moving targets,like opening-closing doors or people that are walking. Ateach time instant, the environment can be represented as acollection of the discrete scatterers, which correspond to thepoints of interaction of the radio waves with the environment.These interactions include specular reflections of the signaloff the walls and other flat surfaces, and diffuse scattering offthe multifaceted objects or objects with rough surfaces [42].Tracking time evolution of the interaction points’ positions canhelp to reconstruct the geometry of the considered environment[43].
B. Fundamentals and Rationale
Stemming from the basic system model in the previoussections, we consider the sequence of UE or sensing positionsdenoted as p UE l = (cid:0) x UE l , y UE l (cid:1) , l = 0 , . . . , L − , while nowalso associate the time instant t l corresponding to the UEposition p UE l . In the continuation, for notational simplicity,we use the index l to refer to a time-dependent value at timeinstant t l . Following the previous system model, the signalemitted by the transmitter at the position p TX l is reflectedoff the K interaction points with absolute 2D coordinates p SC l,k = ( x SC l,k , y SC l,k ) , k = 0 , . . . , K − . The AoD of the signalfrom the TX to the k th interaction point is ϕ TX l,k . The reflectedsignal from each interaction point k reaches the RX located atthe position p RX l = ( x RX l , y RX l ) from the AoA of ϕ RX l,k .In general, one can expect to observe a mixture of stationaryand moving scatterers. Moving scatterers represent movingobject as well as interaction points corresponding to thespecular reflections. For the specular reflections, there is anunambiguous relationship between the incidence and reflectionangles, as ϕ ref = ϕ inc . Hence, the location of specular reflectionpoint at each time instant t l is a function of the positions ofTX and RX, while the motion of the radar causes the shiftof the specular interaction point. In this work, we considerthe continuous white noise acceleration (CWNA) motionmodel [44], implying that the object moves with constantvelocity perturbed with a white noise process, to represent thekinematics of the moving scatterers.In case of diffuse scattering, on the contrary, a relativelywide distribution of the reflection angles may correspond toeach incident angle. Due to this dispersion of the reflectedwave, the apparent location of the interaction point remainsconstant for a range of transmitter and receiver positions { p TX l k , p RX l k } , l k = l , ..., l ν , where ν is the range of radar motion steps,within which the distributions of the corresponding reflectionangles intersect, i.e., ϕ ref l ∩ ... ∩ ϕ ref l ν (cid:54) = ∅ . To describe thekinematics of the stationary scatterers, we consider scattererposition to be constant, perturbed with a white noise process. Inaccordance with naming convention of [44, Ch. 6], we call thisapproach a continuous white noise velocity (CWNV) modelin the continuation. Since the impinging radio waves of the moving radar mayilluminate different portions of the environment in consecu-tive time instants, depending in general on the environmentgeometry, it is natural to assume that stationary and movingscatterers may also mutually interchange their roles. To accountfor such potential time-evolution of the motion model, weformulate in the following Sub-section V.C an IMM algorithm[44] based approach to track the positions of the scatterers. TheIMM approach is generally known to yield good performancein tracking maneuvering targets by running several differentmodels in parallel, while using the same observation data foran update stage, and fusing the results at each time step. Assub-models in our work, the considered IMM utilizes extendedKalman filters (EKFs) [45], [46] incorporating the noted CWNAand CWNV motion models for different types of trackedscatterers. These are described in Sub-section V.D. Furthermore,in order to cross-identify the scatterers between consecutivetime instants, we describe a proximity-based data associationalgorithm in Sub-section V.E that utilizes the tracked covariancematrix of the scatterer location. Finally, an IMM smoothingstage is described in Sub-section V.F. C. IMM Framework
In the IMM-based approach, each of the K scatterers is ingeneral tracked separately, but we omit the index in some of themath presentation below for notational simplicity. Now, let usassume that there is a set of prior probabilities associated withmodels { X , ..., X W } , µ j = P { X j } , j = 1 , ..., W where W denotes the number of considered (sub)models. The probability µ j tells how likely it is for a motion of a random scatterer k to be described by the model j at the initial time instant.Also we assume that probabilities to switch from model i to model j at the next time instant are known and denotedby p ij = P { X jl | X il − } , i, j = 1 , .., W . Both µ j and p ij areIMM filter design parameters and are chosen to reflect theproperties of environment. The IMM filter can be consideredas a Markovian switching system with probability transitionmatrix [ p ij ] ∈ R W × W . At each time step, one obtains theinitial conditions for each model in { X , ..., X W } using allprevious state estimates (the so-called mixing stage) while thenpropagates the states using the respective EKF filter equationsdescribed in the next sub-section. Also the probabilities of eachmodel are updated at every time step. The IMM combined stateestimate and the covariance are then calculated as a weightedmean of a-posteriori EKF state estimates and the covariancematrices of all the models, where the weights are defined bythe probabilities [46].Specifically, during the so-called interaction step, the mixingprobabilities µ i | j for all model combinations are calculated,including also the normalization coefficient for each model ¯ c i ,expressed as ¯ c j = W (cid:88) i =1 p ij µ il − , and µ i | jl = 1¯ c j p ij µ il − . (31)The mixing probabilities are then used to calculate the mixedinput state s jl − and covariance matrices C jl − for each EKF sub-filter, as s jl − = W (cid:88) i =1 µ i | j s il − , C jl − = W (cid:88) i =1 µ i | j (cid:18) C il − + (cid:104) s il − − s jl − (cid:105) (cid:104) s il − − s jl − (cid:105) T (cid:19) , (32)where s il − and C il − are the a-posteriori state estimate andcovariance, respectively, of the i th model at the time moment l − .In the next filtering stage, each model i runs independently,propagating state s il − and covariance C il − through EKFprediction and update (c.f. Sec. V-D). The a-posteriori estimatesresulting from each model at time instant l are s il and C il ,respectively. The likelihood of each of the measurements isalso computed at this stage, expressed as Λ il = N ( r il ; 0 , S il ) , (33)where r il is the measurement residual calculated at the EKFupdate step, shown in (42), while S il is measurement covarianceat time instant l . This allows to calculate the probability ofeach model µ il as µ il = 1 (cid:80) Wi =1 Λ il ¯ c i Λ il ¯ c i . (34)The combined IMM state s l and covariance C l estimates arethen computed at the next stage, the so-called combination stage, as s l = W (cid:88) i =1 µ il s il , C l = W (cid:88) i =1 µ il (cid:16) C il + (cid:2) s il − s l (cid:3) (cid:2) s il − s l (cid:3) T (cid:17) . (35)While the above presentation is written in general form, wenote that in this paper W = 2 , and as the corresponding IMMsub-models { X , X } we consider two EKFs utilizing CWNVor CWNA motion models for the state dynamics. These aredescribed next. D. EKFs and Motion Models
Let us next denote the radar-based measurement model vectorfor the k -th scatterer obtained at the l -th UE location p UE l as m l,k = [ ϕ UE l,k , d UE l,k ] T , (36)where ϕ UE l,k and d UE l,k are the angle and distance estimatesobtained as described in Section III, respectively. Measurementcovariance R is a function of the accuracy of the angle andrange measurements while the choice of the measurement func-tion is an EKF design parameter. To this end, the measurements [ ϕ UE l,k , d UE l,k ] T are related to the locations of the scatterer and locations of the TX and RX via the non-linear measurementequation of the form ϕ UE l,k d UE l,k = h ( p SC l,k ) = atan2 (cid:16) y SC l,k − y UE l , x SC l,k − x UE l (cid:17) || p TX l − p SC l,k || + || p RX l − p SC l,k || , (37)where h ( p SC l,k ) is a differentiable observation function.The state vector at time instant t l , corresponding to the UElocation p UE l , is generally denoted as s l . Then, depending onthe considered motion model, the state vector includes 2Dposition of scatterer (the CWNV model), or scatterer’s positionand velocity (the CWNA model), written as s X l,k = [ x X l,k , y X l,k ] T , (38) s X l,k = [ x X l,k , y X l,k , ˙ x X l,k , ˙ y X l,k ] T . (39)We note that in the following EKF equations, we omit thesub-model related superscripts X and X , where feasible, forthe simplicity of notation, while specifically differentiate themodels where it is necessary to distinguish between them. Tothis end, in the considered EKFs, the state-transition model islinear while the measurement model is non-linear which weexpress as s l,k = Fs l − ,k + u l,k , m l,k = h ( s l,k ) + e l,k , (40)where F is a state-transition matrix, u is the state-processnoise, u l,k ∼ N (0 , Q ) , and e is the measurement noise, e l,k ∼N (0 , R ) . Expressions for F and state-process noise covariance Q for both models can be obtained from the discretization ofthe continuous state-transition model [45], [46]. The powerspectral density Q c of the state-process noise is again one ofthe EKF design parameters.The EKF prediction step is then propagating the a-posterioristate estimate for the k -th scatterer at time instant t l , ˆ s + l,k , andcovariance estimate ˆ C + l,k in accordance to its correspondingmotion model to obtain the a-priori estimates at time instant t l , written as ˆ s − l,k = F ˆ s + l − ,k , ˆ C − l,k = F ˆ C + l − ,k F T + Q . (41)Note that the sizes of the state-transition and covariance matri-ces as well as those of their estimates are defined by the usedmotion model, i.e., C M , ˆ C M , F M , Q M , R M ∈ R × ; C M , ˆ C M , F M , Q M , R M ∈ R × .Finally, the EKF update step can be expressed as K l,k = ˆ C − l,k H Tl,k (cid:16) H l,k ˆ C − l,k H Tl,k + R (cid:17) − , C + l,k = ( I − K l,k H l,k ) ˆ C − l,k , r l,k = (cid:16) m l,k − h (cid:16) s − l,k (cid:17)(cid:17) , s + l,k = s − l,k + K l,k r l,k , (42)where K l,k denotes the Kalman gain, and H l,k denotes theJacobian matrix of the measurement model (37) evaluatedat the a-priori estimation of the scatterer location, and themeasurement residual is denoted as r l,k . E. Scatterer Measurement Association
We observe that due to the measurement noise, complicatedsystem geometry and the presence of the double reflections,cross-identification of the reflection points between consecutiveobservations may become a challenge. In this paper, we usea proximity argument to conclude whether two scatterersobserved at two consecutive time instants are in fact the samespecular or diffuse reflection point. We consider relatively slowmotion of the radar or frequent measurement updates so thatan apparent position of a reflection point is not expected tochange dramatically during one time step.Let us assume that there are K l − scatterers available at timeinstant l − , and K l at the next time instant. In general, amongthose K l scatterers, there might be a number of scatterersobserved at the previous time instant as well as newly observedones. Likewise, some of the K l scatterers observed at theprevious time instant may have become obstructed, or thereflected signal level is too weak, so that those scatterers areno longer available for tracking at time instant l .Let us assume that { ˆ p l − ,k } , k = 1 , ..., K l − , are the IMMestimates (cf. Eq. (35)) of all scatterer positions that have beenidentified at time instant l − . Let us also consider coordinatesof k -th scatterer that have been crudely estimated from theobserved angle and distance at a time instant l + 1 before usingthe method described in Section III-C, expressed as ˜ p l,k = (˜ x l,k l , ˜ y l,k l ) = (cid:0) d UE l,k l cos( ϕ UE l,k l ) , d UE l,k l sin( ϕ UE l,k l ) (cid:1) , (43)and that there are K l such estimates available, { ˜ p l,k } , k =1 , ..., K l . In order to associate the newly measured scatterersat step l with those observed and tracked at the previous stepand continue or discontinue tracking, an association distanceis calculated for each pair of old and new scatterers as D ak,q = || ˜ p l,k − ˆ p l − ,q || , k = 1 , ..., K l , q = 1 , ..., K l − . (44)For each newly identified scatterer q , its closest proximitycounterpart is chosen as a candidate ”parent”, expressed as k ∗ ( q ) = arg min k ( D ak,q ) . (45)At the second step in scatterer association process, we checkif ˜ p l +1 ,k is within the P -percent confidence ellipse of IMMweighted position estimate ˆ p l,k with covariance ˆ C l,k , where P is a design parameter. For the associated scatterers, we continuethe tracking process. If several new scatterers are associatedwith the same ”parent”, we continue tracking only the closestone, all others are dropped. For the newly identified scattererswhich are not associated with any of the ”parents”, we initializethe IMM EKF filters. If a ”parent” scatterer has no candidates,we use the prior position estimate for the next time instant,and drop it from the tracking process if there is no secondassociation in a row. Also other data association approaches,such as the Hungarian algorithm, were considered, but withoutobserved performance improvement. F. IMM Smoothing
In order to utilize all the available radar-based measurementsfor building a map of the underlying environment at the end of a measurement campaign, we consider an optional smoothingstep for the IMM results. From a variety of available IMMsmoothing algorithms, we have chosen the approach proposedby [47] as it is relatively straightforward and does not requirestoring of the radar distance and angle measurements. In theused approach, the IMM smoothing mimics the behaviour ofthe forward IMM propagation, with each of the interactingIMM sub-models using the so-called Rauch-Tung-Striebelsmoothing recursion [48] combined with the model interaction.This becomes possible due to approximation of the backwardsmodel transition probabilities by the forward model transitionprobabilities, stemming from the Markov properties of themodel transition, as well as approximation of the smoothedprobability density with the Gaussian mixture of W modeconditioned smoothing densities. Such a mixture can beconsidered a sufficient statistics of the measurement set.As the very last stage, a thresholding operation is imposedon the refined covariances of the smoothed scatterer positionestimates in order to remove the unreliable estimates. Wenote that the covariance of the scatterer’s position is relatedto the SNR and to the accuracy of the angle and distancemeasurements. Hence, this approach allows us to use allavailable measurement information, automatically discardingthe unreliable position estimates. This will be concretelyillustrated in Section VII through the processing of both ray-tracing based and the actual RF measurement data.VI. E VALUATION S CENARIO AND E NVIRONMENTS
In this section, the evaluation scenario, ray-tracing environ-ment as well as the experimental RF measurement environmentand equipment, used for the validation of the proposed sensingand mapping algorithms are described. It is also noted thatthe complete I/Q measurement data is openly available athttps://doi.org/10.5281/zenodo.4475160 [49].
A. Scenario Description
The evaluation scenario is an indoor office environment atthe Hervanta Campus of Tampere University, Finland, as shownin Fig. 3(a). The considered environment consist of a corridorof some 2 m wide and 60 m long with different office roomson both sides as illustrated by the line art overlaid to all thefollowing mapping results.The environment is sensed along three parallel trajectories asshown in Fig. 3(a). In the center trajectory, the sensing relatedmeasurements are conducted with distance step of 0.5 m, whilethe side trajectories have a step of 1 m. In addition, thesepositions are deployed in both directions along the corridor,providing a total of six sets of measurements. In Fig. 3(a),the considered measurement locations as well as the mostsignificant targets from the radar perspective are shown. We canhighlight three walls of adjacent corridors that are perpendicularto the system trajectory— A , B and C —located at the left sideof the figure. Moreover, at the right side, three metal lockers— D , E and F —are expected to be the main reflective targets dueto their notable RCS. B. 5G NR Waveform
In both the RF measurements and the ray tracing simulations,OFDM-based NR uplink waveform is used with the widest available channel bandwidth at the mm-wave frequency rangeaccording to [21]. Therefore, a channel bandwidth of 400 MHzwith subcarrier spacing of 120 kHz is utilized. In particular,for each sensing location and scanning direction, we consideran uplink NR frequency-domain resource grid with N = 3168 active subcarriers and M = 28 OFDM symbols correspondingto an observation window of around 0.25 ms. In this case,the M consecutive OFDM symbols are coherently combinedto improve the signal-to-noise ratio (SNR) of the obtainedrange–angle charts as described in (24). According to [50], theconsidered transmit waveform provides a basic range resolutionof about 40 cm.It is also shortly noted that some recent studies, such as [51]–[54], have raised the idea of joint waveform optimization toimprove the sensing performance in JCAS type of systems. Inthis work, however, we deliberately use 3GPP 5G NR standardcompliant uplink waveform with all physical channels andsignal structures involved, to reflect the true waveform of NRUEs as accurately as possible. C. Ray-Tracing Environment
For validation purposes, the RF measurement campaignresults are corroborated by realistic ray-tracing simulationsusing Wireless Insite® [55] software. In these simulations,reproducing the scenario shown in Fig. 3(a), the measurementdevice follows a similar trajectory as in the RF measurementcampaign with 31 test locations as shown in Fig. 5. The TXand RX array operation is simulated through a directive beampattern with 3 dB beam-width of ◦ , similar to the realdirective antenna systems used in the RF measurements. Inaddition, the same antenna separation and height are used, andthe carrier frequency is GHz.The simulator is configured to consider a maximum of 15rays per simulation, and the number of allowed interactions islimited to six reflections and one diffraction. The walls, floorand ceiling are built using the frequency-specific materials,namely, ITU layered dry wall and floor or ceiling board. Somereasonable simplifications in the environment modeling wereallowed, compared to the true physical environment shownin Fig. 3(a), to reduce the computational complexity andsimulation time. However, like the results will also illustrate,the ray tracing environment models accurately the physicalenvironment.
D. RF Measurement Equipment
In the actual RF measurements, we use state-of-the-artmm-wave equipment to emulate the UE operation at the28 GHz band, with selected equipment shown in Fig. 3(b). Thebaseline hardware platform in the measurement setup is theNI vector signal transceiver (PXIe-5840) which implementsthe RF transmitter and receiver functionalities at intermediatefrequency of 3.5 GHz, as well as controls the rest of the devices.In addition, two signal generators (Keysight N5183B–MXG)are used as local oscillators, which together with externalmixers (Marki Microwave T3-1040) up- and down-convert theIF signal to and from the desired mm-wave carrier frequency.To emulate the UE’s phased array beam-steering operation,two directive horn antennas (PE9851A-20) are utilized for the B D C E A F (a) Evaluation scenario TXRX (b) RF measurement equipmentFig. 3. (a) The evaluation scenario including the main sensing locations andthe most important radar targets. (b) The main equipment used in the actualRF measurements at 28 GHz.
TX and RX sides. These antennas are mounted on mechanicalsteering systems which enable to steer and direct the hornsin the whole azimuth plane very accurately. According to thespecifications, both horns provide a nominal gain of 20 dBiwith a 3 dB beam width of 17 ◦ . The antennas are placed atone meter above the floor level with a separation of 60 cmin order to avoid larger mutual coupling between TX and RXchains. In the TX side, two external power amplifiers (PAs)are also used that together with the antenna system facilitatean EIRP of around +20 dBm.VII. R ESULTS
In this section, we provide the ray-tracing and RF measure-ment based results. We start with a short example of range-angleprocessing, utilizing the RF measurement data, while then putmost of the focus on the actual mapping results, covering boththe ray-tracing data and the RF measurement data.
A. Example Range–Angle Processing Results
We start by demonstrating the proposed range–angle process-ing technique described in Section III using the 28 GHz RFmeasurement data. Specifically, the 10 th location of the centertrajectory shown in Fig. 6(e) is considered as a representativeexample. First, Fig. 4(a) illustrates the resulting range–anglechart when applying the considered LS estimation approach.In this case, distances up to 30 m and azimuth angles between (a) Range–angle chart via LS approach(b) Range–angle chart via ISTA approachFig. 4. Example range–angle charts obtained via (a) the LS approach and (b)the ISTA approach, when processing the 28 GHz RF measurement data atlocation highlighted in green in Fig. 6(e). − ◦ to ◦ are investigated, using a total of C R = 391 and C ϕ = 221 range and angle cells, respectively. In addition, aHamming window is used to improve the side-lobe suppressionin the range domain. As it can be observed, the 400-MHzbandwidth used in the measurements facilitates high-accuracyrange resolution, that enables to distinguish targets with mutualdistances down to around 0.4 m. However, we also identifysignificant side-lobes in the angular-domain due to the TX andRX horn antenna patterns that degrade the target detection inthis domain.Then, a sparse representation of this range–angle chart isobtained by applying the proposed ISTA approach, summarizedin Algorithm 1, to better facilitate the subsequent mappingphases. The corresponding results are shown in Fig. 4(b). Asit can be seen, only the most significant target reflectionsremain in the sparse chart, suppressing the possible side-lobesin both range and angular domains. Hence, the ISTA approachis deployed as the main range–angle chart processing enginein the forthcoming mapping results. B. Ray-Tracing based Mapping Results
The proposed grid-based static and tracking-based dynamicmapping methods are next assessed and validated with theray-tracing based evaluations. As already noted, the ray-tracingenvironment resembles the true physical environment and thecorresponding RF measurements as closely as possible – withan exception that we adopt a wider angular scanning range from − ◦ to ◦ in order to seek for the maximum mappingperformance. The obtained results with the different mappingmethods are presented and illustrated in Fig. 5, including alsothe building floor plan for reference. (a) Grid-based static mapping(b) Tracking-based dynamic mapping(c) Smoothed tracking-based dynamic mappingFig. 5. Indoor mapping results with ray-tracing data at 28 GHz, for (a) grid-based static mapping, (b) tracking-based dynamic mapping without smoothingand (c) tracking-based dynamic mapping with smoothing. First, the grid-based static mapping approach, described inSection IV, is shown in Fig. 5(a). This figure represents thefinal grid-based map after applying the different averaging,filtering and thresholding stages. Specifically, we pursue a 2Dmap of the simulated corridor using square cells with size of . × . m which is consistent with the range resolution ofabout 0.4 m, to create the initial average map in (27). Next, aGaussian kernel matrix with parameters U = V = 5 and σ = 1 is applied to smooth the map and reduce the effects of noiseas described in (28). Finally, a thresholding stage is deployedto emphasize the most relevant targets of the environment. Asit can be observed, the grid-based method is able to fairlyaccurately sense the indoor scenario, providing a mappingreconstruction that clearly reflects the true layout. However,it can also be observed and noted that there are some clearartifacts in the constructed map, calling and motivating for themore advanced processing methods.Next, the same ray-tracing scenario and data is processedwith the tracking-based dynamic mapping approach, describedin Section V. The corresponding results are shown in Fig. 5(b)and (c), without and with IMM smoothing, respectively. Theseresults consider the IMM-based method which deploys two (a) Grid-based static mapping (moving from left to right) (b) Grid-based static mapping (moving from right to left)(c) Tracking-based dynamic mapping (moving from left to right) (d) Tracking-based dynamic mapping (moving from right to left)(e) Smoothed tracking-based dynamic mapping (moving from left to right) (f) Smoothed tracking-based dynamic mapping (moving from right to left)Fig. 6. Indoor mapping results with RF measurement data at 28 GHz. Subfigures (a), (c) and (e) illustrate different mapping methods for the center trajectorywhen the measurement equipment is moving from left to right along the corridor. Similarly, subfigures (b), (d) and (f) present mapping results when themeasurement system senses the corridor while moving from right to left. separate EKF filters using the CWNV and CWNA motionmodels for stationary and moving scatterers, respectively. Forthe initial state covariance in the EKF models, we use a varianceof 0.5 for the target x and y coordinates, as well as for thetarget velocities in x and y directions when considering theCWNA model. The used power spectral densities of the statecovariance matrices Q c of the CWNV and CWNA models are1 and 3, respectively. The initial probabilities of the scatterermotion models ( µ j ) for the CWNA and CWNV are given as0.9 and 0.1, respectively.In the forward IMM filter pass, the scatterer positionestimates are tracked together with their covariance matrices.At each step, they are updated using the newly acquired angleand distance measurements. To this end, Fig. 5(b) shows allthe tracked scatterer positions through the UE trajectory. Thepositions corresponding to the same scatterer cross-identifiedbetween multiple time instants are interconnected formingtrajectories that are shown to follow the walls of the corridorwith notable accuracy. Note that the covariance matrices ofthe scatterer positions become smaller as more measurementare available, hence, the position estimates in the beginning of tracking are not as reliable as those closer to the end. Also,we can still observe some apparent scatterers beyond the walls,corresponding, e.g., to significant double reflections.After the UE has finished the tracking process, it provides thetracked targets for post-processing, i.e., to the IMM smootherthat is deployed to further improve the final map representation.During the backwards IMM smoothing pass, the scatterers’position and covariance estimates are updated taking intoaccount all measurements available at the last tracking timeinstant. In the final map processing step, we use the smoothedcovariance of each position estimate as a measure of reliability,allowing thus to extract only the most reliable ones. Thecorresponding final map shown in Fig. 5(c), demonstratesthat a number of scatterers did not pass to the final stage dueto their large covariance, especially those in the beginning ofthe tracking process, those far away from the transmitter andthose corresponding to the multiple reflections. C. RF Measurement based Mapping Results
Finally, the proposed methods are assessed and tested withRF measurements to validate the considered sensing and (a) Smoothed tracking-based dynamic mapping (from left to right)(b) Smoothed tracking-based dynamic mapping (from right to left)Fig. 7. Combined smoothed tracking-based dynamic mapping for three showntrajectories when the measurement equipment moves (a) from left to right and(b) from right to left along the corridor. mapping functionality with real-world equipment and physicalenvironment. We address how the mapping system performanceis subject to the UE orientation by showing mapping results fortwo main UE trajectories separately. To this end, Fig. 6(a), (c)and (e) compare the grid-based static and the tracking-baseddynamic mapping results when the measurement equipment ismoving along the corridor from left to right, while Fig. 6(b),(d) and (f) show the corresponding mapping results for theopposite moving direction, from right to left.It can be seen how the proposed methods very efficientlyand accurately reconstruct the complex physical environmentdespite the fairly limited scanning range of − ◦ to ◦ ofthe available antenna systems, compared to the ray-tracingscenario, and despite the numerous real-world effects andimpairments that are present in the measurement data. Weespecially highlight the locations of the main corridor walls(with A , B and C ) and the metallic lockers (with D , E and F )locations, both in Fig. 3 and in Fig. 6, for better interpretationof the results. We can observe how the tacking-based dynamicapproach is able to accurately track the main targets of theenvironment, and subsequently, improve the map quality –especially when the IMM smoothing stage is deployed.Finally, we showcase that by combining the data fromthe multiple measurement routes, a more rich and accuraterepresentation of the environment can be obtained. In thisregard, Figure 7 illustrates the combined mapping results withIMM smoothing filter for all the three different trajectories.As it can be seen, the combination of multiple UE trajectoriesprovide a more complete representation of the indoor map incomparison with the maps shown in Fig. 6 for a single trajectory.Overall, the mapping results with real-world RF measurementdata demonstrate the applicability of the proposed methodsalso in true complex physical environments. VIII. C ONCLUSIONS
This paper investigated the sensing and environment mappingprospects of mm-wave 5G NR and beyond networks withspecific emphasis on the UE side for indoor mapping applica-tions. In the considered framework, the UE operates as a jointcommunication and sensing system, scanning its surroundingenvironment while steering its beam pattern towards differentdirections and observing the target reflections. First, we deriveda novel LS-based processing technique to estimate sparse range–angle charts to facilitate the subsequent mapping processing.Then, two different mapping approaches were presented –a static grid-based approach and a more advanced dynamictracking-based approach. While the grid-based approach servesas a computationally simple baseline, the dynamic tracking-based approach building on IMM-EKF processing solutionincorporates both specular reflections and diffuse scatteringwhile allowing to track individual scatterers over time. Thistogether with the IMM smoothing and described scatterermeasurement association approach provide an efficient frame-work to reconstruct and map complex physical environments.The applicability of the proposed methods and algorithmswas then assessed and evaluated, through both ray-tracingsimulations and actual RF measurements at the 28 GHz bandin practical indoor office type of an environment. The obtainedresults demonstrate the good applicability of the proposedmethods, while the IMM-EKF based dynamic tracking solutionwas shown to clearly outperform the more simple grid-basedapproach. Our future work will consider extensions to 3Dsensing and mapping, while extending also to the actual SLAMwhere the coordinates of the sensing device are unknown.R
EFERENCES[1] H. Holma, A. Toskala, and T. Nakamura,
5G Technology: 3GPP NewRadio . Wiley, 2019.[2] H. Wymeersch, G. Seco-Granados, G. Destino, D. Dardari, and F. Tufves-son, “5G mmWave positioning for vehicular networks,”
IEEE WirelessCommunications , vol. 24, no. 6, pp. 80–86, 2017.[3] M. Koivisto et al. , “High-efficiency device positioning and location-aware communications in dense 5G networks,”
IEEE CommunicationsMagazine , vol. 55, no. 8, pp. 188–195, 2017.[4] J. Talvitie, M. Koivisto, T. Levanen, M. Valkama, G. Destino, andH. Wymeersch, “High-accuracy joint position and orientation estimationin sparse 5G mmwave channel,” in
Proc. IEEE ICC 2019 , 2019, pp. 1–7.[5] “3GPP TS 22.261 V18.1.1, ”Service requirements for the 5G system”,Tech. Spec. Group Services and System Aspects, Rel. 18,” Jan. 2021.[6] F. Liu, C. Masouros, A. P. Petropulu, H. Griffiths, and L. Hanzo, “Jointradar and communication design: Applications, state-of-the-art, and theroad ahead,”
IEEE Transactions on Communications , vol. 68, no. 6, pp.3834–3862, 2020.[7] B. Paul, A. R. Chiriyath, and D. W. Bliss, “Survey of RF communicationsand sensing convergence research,”
IEEE Access , vol. 5, pp. 252–270,2017.[8] L. Zheng, M. Lops, Y. C. Eldar, and X. Wang, “Radar and communicationcoexistence: An overview: A review of recent methods,”
IEEE SignalProcessing Magazine , vol. 36, no. 5, pp. 85–99, 2019.[9] A. R. Chiriyath, B. Paul, and D. W. Bliss, “Radar-communicationsconvergence: Coexistence, cooperation, and co-design,”
IEEE Trans.Cognitive Commun. Netw. , vol. 3, no. 1, pp. 1–12, 2017.[10] A. Hassanien, M. G. Amin, E. Aboutanios, and B. Himed, “Dual-functionradar communication systems: A solution to the spectrum congestionproblem,”
IEEE Signal Processing Mag. , vol. 36, no. 5, pp. 115–126,2019.[11] M. L. Rahman, J. A. Zhang, X. Huang, Y. J. Guo, and R. W. Heath,“Framework for a perceptive mobile network using joint communicationand radar sensing,”
IEEE Transactions on Aerospace and ElectronicSystems , vol. 56, no. 3, pp. 1926–1941, 2020. [12] M. L. Rahman, J. A. Zhang, K. Wu, X. Huang, Y. J. Guo, S. Chen,and J. Yuan, “Enabling joint communication and radio sensing inmobile networks – a survey,” CoRR , 2020. [Online]. Available:https://arxiv.org/abs/2006.07559[13] C. Baquero Barneto et al. , “Full-duplex OFDM radar with LTE and5G NR waveforms: Challenges, solutions, and measurements,”
IEEETransactions on Microwave Theory and Techniques , vol. 67, no. 10, pp.4042–4054, 2019.[14] C. Baquero Barneto, S. D. Liyanaarachchi, M. Heino, T. Riihonen,and M. Valkama, “Full-duplex radio/radar technology: The enabler foradvanced joint communication and sensing,”
IEEE Wireless Communica-tions , 2020, in press.[15] J. Choi, V. Va, N. Gonzalez-Prelcic, R. Daniels, C. R. Bhat, and R. W.Heath, “Millimeter-wave vehicular communication to support massiveautomotive sensing,”
IEEE Communications Magazine , vol. 54, no. 12,pp. 160–167, 2016.[16] C. De Lima et al. , “Convergent communication, sensing and localiza-tion in 6G systems: An overview of technologies, opportunities andchallenges,”
IEEE Access , pp. 1–24, 2021.[17] M. Alloulah and H. Huang, “Future millimeter-wave indoor systems:A blueprint for joint communication and sensing,”
Computer , vol. 52,no. 7, pp. 16–24, 2019.[18] P. Kumari, J. Choi, N. Gonz´alez-Prelcic, and R. W. Heath, “IEEE802.11ad-based radar: An approach to joint vehicular communication-radar system,”
IEEE Transactions on Vehicular Technology , vol. 67, no. 4,pp. 3012–3027, 2018.[19] F. Guidi, A. Guerra, and D. Dardari, “Personal mobile radars withmillimeter-wave massive arrays for indoor mapping,”
IEEE Transactionson Mobile Computing , vol. 15, no. 6, pp. 1471–1484, 2016.[20] C. Baquero Barneto, M. Turunen, S. D. Liyanaarachchi, L. Anttila,A. Brihuega, T. Riihonen, and M. Valkama, “High-accuracy radio sensingin 5G new radio networks: Prospects and self-interference challenge,” in ,2019, pp. 1159–1163.[21] “3GPP TS 38.104 v15.4.0, ”NR; Base Station (BS) radio transmissionand reception”, Tech. Spec. Group Radio Access Network, Rel. 15,” Dec.2018.[22] T. Levanen, O. Tervo, K. Pajukoski, M. Renfors, and M. Valkama,“Mobile communications beyond 52.6 GHz: Waveforms, numerology,and phase noise challenge,”
IEEE Wireless Communications , pp. 1–8,2020.[23] T. S. Rappaport, G. R. MacCartney, M. K. Samimi, and S. Sun, “Wide-band millimeter-wave propagation measurements and channel modelsfor future wireless communication system design,”
IEEE Transactionson Communications , vol. 63, no. 9, pp. 3029–3056, 2015.[24] A. Guerra, F. Guidi, D. Dardari, A. Clemente, and R. D’Errico, “Amillimeter-wave indoor backscattering channel model for environmentmapping,”
IEEE Transactions on Antennas and Propagation , vol. 65,no. 9, pp. 4935–4940, 2017.[25] A. Yassin, Y. Nasser, A. Y. Al-Dubai, and M. Awad, “MOSAIC:Simultaneous localization and environment mapping using mmwavewithout a-priori knowledge,”
IEEE Access , vol. 6, pp. 68 932–68 947,2018.[26] Y. Ji, J. Hejselbæk, W. Fan, and G. F. Pedersen, “A map-free indoorlocalization method using ultrawideband large-scale array systems,”
IEEEAntennas and Wireless Propagation Letters , vol. 17, no. 9, pp. 1682–1686,2018.[27] R. Mendrzik, H. Wymeersch, and G. Bauch, “Joint localization andmapping through millimeter wave MIMO in 5G systems,” in , 2018, pp. 1–6.[28] Q. Huang, H. Chen, and Q. Zhang, “Joint design of sensing andcommunication systems for smart homes,”
IEEE Network , vol. 34, no. 6,pp. 191–197, 2020.[29] F. Guidi, A. Mariani, A. Guerra, D. Dardari, A. Clemente, andR. D’Errico, “Indoor environment-adaptive mapping with beamsteeringmassive arrays,”
IEEE Transactions on Vehicular Technology , vol. 67,no. 10, pp. 10 139–10 143, 2018.[30] F. Guidi, A. Guerra, D. Dardari, A. Clemente, and R. D’Errico,“Environment mapping with millimeter-wave massive arrays: Systemdesign and performance,” in , 2016, pp. 1–6.[31] C. Baquero Barneto, T. Riihonen, M. Turunen, M. Koivisto, J. Talvitie,and M. Valkama, “Radio-based sensing and indoor mapping withmillimeter-wave 5G NR signals,” in , 2020, pp. 1–5.[32] M. Braun, “OFDM radar algorithms in mobile communication networks,”
Ph. D. dissertation , 2014. [33] C. Sturm and W. Wiesbeck, “Waveform design and signal processingaspects for fusion of wireless communications and radar sensing,”
Proceedings of the IEEE , vol. 99, no. 7, pp. 1236–1259, 2011.[34] I. T. Cummings, J. P. Doane, T. J. Schulz, and T. C. Havens, “Aperture-level simultaneous transmit and receive with digital phased arrays,”
IEEETransactions on Signal Processing , vol. 68, pp. 1243–1258, 2020.[35] J. A. Zhang, X. Huang, Y. J. Guo, J. Yuan, and R. W. Heath, “Multibeamfor joint communication and radar sensing using steerable analog antennaarrays,”
IEEE Transactions on Vehicular Technology , vol. 68, no. 1, pp.671–685, 2019.[36] F. Liu, C. Masouros, A. Li, H. Sun, and L. Hanzo, “MU-MIMO com-munications with MIMO radar: From co-existence to joint transmission,”
IEEE Transactions on Wireless Communications , vol. 17, no. 4, pp.2755–2770, 2018.[37] X. Liu, T. Huang, N. Shlezinger, Y. Liu, J. Zhou, and Y. C. Eldar, “Jointtransmit beamforming for multiuser MIMO communications and MIMOradar,”
IEEE Trans. Signal Processing , vol. 68, pp. 3929–3944, 2020.[38] Y. Luo, J. A. Zhang, X. Huang, W. Ni, and J. Pan, “Multibeamoptimization for joint communication and radio sensing using analogantenna arrays,”
IEEE Transactions on Vehicular Technology , vol. 69,no. 10, pp. 11 000–11 013, 2020.[39] I. Daubechies, M. Defrise, and C. De Mol, “An iterative thresholdingalgorithm for linear inverse problems with a sparsity constraint,”
Com-munications on Pure and Applied Mathematics: A Journal Issued bythe Courant Institute of Mathematical Sciences , vol. 57, no. 11, pp.1413–1457, 2004.[40] A. Beck and M. Teboulle, “A fast iterative shrinkage-thresholdingalgorithm for linear inverse problems,”
SIAM journal on Imaging Sciences ,vol. 2, no. 1, pp. 183–202, 2009.[41] M. A. T. Figueiredo and R. D. Nowak, “An EM algorithm for wavelet-based image restoration,”
IEEE Transactions on Image Processing ,vol. 12, no. 8, pp. 906–916, 2003.[42] C. B. N. Pinel,
Electromagnetic wave scattering from random roughsurfaces : asymptotic models . Wiley, 2013.[43] M. Oren and S. Nayar, “A theory of specular surface geometry,”
International Journal of Computer Vision (IJCV) , vol. 24, 10 2002.[44] T. K. Y. Bar-Shalom, X. Rong Li,
Estimation with Applications toTracking and Navigation: Theory, Algorithms and Software . WileyInterscience, 2001.[45] S. Kay,
Fundamentals of Statistical Signal Processing: Estimation Theory .Prentice-Hall Signal Processing Series, 1993.[46] J. Hartikainen, A. Solin, and S. S¨arkk¨a,
Optimal Filtering with KalmanFilters and Smoothers – a Manual for the Matlab toolbox EKF/UKF .Department of Biomedical Engineering and Computational Science, AaltoUniversity School of Science, 2011.[47] N. Nadarajah, R. Tharmarasa, M. McDonald, and T. Kirubarajan,“IMM forward filtering and backward smoothing for maneuvering targettracking,”
IEEE Transactions on Aerospace and Electronic Systems ,vol. 48, no. 3, pp. 2673–2678, 2012.[48] C. T. S. H. E. Rauch, F. Tung, “Maximum likelihood estimates of lineardynamic systems,”
AIAA Journal , vol. 3, no. 8, p. 1445–1450, 1965.[49] C. Baquero Barneto, E. Rastorgueva-Foi, M. F. Keskin, T. Riihonen,M. Turunen, J. Talvitie, H. Wymeersch, and M. Valkama, “Datasetfor radio-based sensing and environment mapping in millimeter-wave5G networks,” Jan. 2021. [Online]. Available: https://doi.org/10.5281/zenodo.4475160[50] C. Baquero Barneto, S. D. Liyanaarachchi, T. Riihonen, L. Anttila, andM. Valkama, “Multibeam design for joint communication and sensingin 5G new radio networks,” in
ICC 2020 - 2020 IEEE InternationalConference on Communications (ICC) , 2020, pp. 1–6.[51] M. Bic˘a and V. Koivunen, “Radar waveform optimization for targetparameter estimation in cooperative radar-communications systems,”
IEEE Transactions on Aerospace and Electronic Systems , vol. 55, no. 5,pp. 2314–2326, 2019.[52] P. Kumari, S. A. Vorobyov, and R. W. Heath, “Adaptive virtualwaveform design for millimeter-wave joint communication–radar,”
IEEETransactions on Signal Processing , vol. 68, pp. 715–730, 2020.[53] F. Liu, C. Masouros, T. Ratnarajah, and A. Petropulu, “On range sidelobereduction for dual-functional radar-communication waveforms,”
IEEEWireless Communications Letters , vol. 9, no. 9, pp. 1572–1576, 2020.[54] S. D. Liyanaarachchi, C. Baquero Barneto, T. Riihonen, and M. Valkama,“Joint OFDM waveform design for communications and sensing conver-gence,” in
Proc. IEEE ICC 2020 S UPPLEMENTARY M ATERIAL :A PPROXIMATION OF R ANGE -A NGLE
LS S
OLUTION
In this supplementary material, we provide a fast approximatesolution to evaluate (23) using the Kronecker structure in (19).Suppose Φ has full column rank – an assumption that is well-justified due to randomness of data symbols x m and largenumber of OFDM subcarriers/symbols so that N IM > C R C ϕ .Then, we can rewrite (23) as (cid:98) b LS = ( Φ H Φ ) − Φ H y = (cid:32) M − (cid:88) m =0 Φ Hm Φ m (cid:33) − M − (cid:88) m =0 Φ Hm y m = (cid:34)(cid:0) G ∗ ⊗ C H (cid:1) M − (cid:88) m =0 diag ( x m ) H diag ( x m ) (cid:0) G T ⊗ C (cid:1)(cid:35) − × (cid:0) G ∗ ⊗ C H (cid:1) M − (cid:88) m =0 diag ( x m ) H y m = (cid:34)(cid:0) G ∗ ⊗ C H (cid:1) M − (cid:88) m =0 diag (cid:0) {| x i,m,n | } i,n (cid:1) (cid:0) G T ⊗ C (cid:1)(cid:35) − × (cid:0) G ∗ ⊗ C H (cid:1) M − (cid:88) m =0 ( x ∗ m (cid:12) y m ) (46) ≈ (cid:2)(cid:0) G ∗ ⊗ C H (cid:1) M I (cid:0) G T ⊗ C (cid:1)(cid:3) − × (cid:0) G ∗ ⊗ C H (cid:1) M − (cid:88) m =0 ( x ∗ m (cid:12) y m ) (47) = (cid:2)(cid:0) G ∗ G T (cid:1) ⊗ (cid:0) C H C (cid:1)(cid:3) − (cid:0) G ∗ ⊗ C H (cid:1) M M − (cid:88) m =0 ( x ∗ m (cid:12) y m )= (cid:104)(cid:0) G ∗ G T (cid:1) − ⊗ (cid:0) C H C (cid:1) − (cid:105) (cid:0) G ∗ ⊗ C H (cid:1) × M M − (cid:88) m =0 ( x ∗ m (cid:12) y m )= (cid:104)(cid:0) G ∗ G T (cid:1) − G ∗ ⊗ (cid:0) C H C (cid:1) − C H (cid:105) × M M − (cid:88) m =0 ( x ∗ m (cid:12) y m )= vec (cid:18) (cid:0) C H C (cid:1) − C H M M − (cid:88) m =0 ( X ∗ m (cid:12) Y m ) (48) × G H (cid:0) GG H (cid:1) − (cid:19) , which yields (cid:98) B LS = (cid:0) C H C (cid:1) − C H M M − (cid:88) m =0 ( X ∗ m (cid:12) Y m ) G H (cid:0) GG H (cid:1) − . (49)Going from (46) to (47), we use the law of large numbers forsufficiently large M by assuming that the constellation has anaverage magnitude of .As seen from (49), the LS solution corresponds to amatched filtering operation. In particular, the observationmatrix is first multiplied by the conjugate of the OFDM transmit symbols to cancel out their effect. Then, the resultingmatrices are coherently integrated over M symbols relyingon the assumption of negligible Doppler over the durationof M symbols. Next, left multiplication by (cid:0) C H C (cid:1) − C H represents range compression/matched filtering via the range-dependent frequency-domain steering matrix C . Notice thatthe columns of C coincide with those of an N × N DFTmatrix. This means that matched filtering via (cid:0) C H C (cid:1) − C H is essentially a normalized IDFT operation over the columnsof X ∗ m (cid:12) Y m , which is a standard approach in OFDM radars[32], [33]. Similarly, right multiplication by G H (cid:0) GG H (cid:1) −1