Member, IEEE Francesco Guidi
email: francesco.guidi@cea.fr

Student, IEEE Anna Guerra
email: anna.guerra3@unibo.it

Senior Member, IEEE ✦ Davide Dardari
email: davide.dardari@unibo.it

Personal Mobile Radars with Millimeter-Wave Massive Arrays for Indoor Mapping

Keywords: Personal radar, indoor mapping, millimeter-wave, massive antenna arrays, Sensors

The adoption of millimeter-wave technology could open the possibility to integrate massive antenna arrays inside future 5G user mobile devices, with the possibility to enable new interesting applications. Within this context, in this paper we put forth the concept of a personal mobile radar operating at millimeter-waves and consisting of a massive array for accurate environmental mapping. Frequency selectivity and phase quantization effects are accounted for to characterize the achievable angle and range resolution necessary to collect environmental information. Successively, we propose an effective grid-based Bayesian mapping approach by introducing a new state-space model, which profits of the beneficial effects of the massive antenna array characteristics. Numerical results show that the idea herein investigated is feasible, and that a significant mapping performance is attainable even employing coarse antenna arrays provided that the number of antenna elements is sufficiently high.

INTRODUCTION

W ITHIN the concept of the Internet of Things, ob- jects and devices are expected to be connected and mapped into the Internet space to enable environmentrelated applications [START_REF] Floerkemeier | The Internet of Things, First Int. Conf[END_REF]. In this context one key step is the digital representation (mapping) of the physical world into the so called smart space (i.e. a space that enables the cooperation and interaction among smart objects). While in the last years this process has rapidly taken place in outdoor environments thanks to the availability of GPS information and the introduction of locationbased applications, other techniques have to be adopted in indoor environments [START_REF] Gjengset | Phaser: Enabling phased array signal processing on commodity WiFi access points[END_REF], [START_REF] Adib | See through walls with wifi![END_REF], whereas manual mapping and updating of all buildings in the world would Fig. 1: Users, equipped with massive arrays on their smartphones, in participatory or not manner scan indoor environments with the opportunity to share their maps. require an unaffordable cost and effort. Traditionally, simultaneous localization and mapping (SLAM) is based on the concept that a robot, moving in an unknown environment, recognizes the surrounding objects being then able to reconstruct a 2D/3D map of the area. Current accurate SLAM, for automatic creation of indoor maps [4]- [START_REF] Durrant-Whyte | Simultaneous localisation and mapping (SLAM): Part I the essential algorithms[END_REF], requires both high-definition distance estimates (ranging) and very narrow steering beams (angle resolution), traditionally accomplished by laser technology. In addition, the exploitation of WiFi and sensors integrated in smartphones for zero-effort automatic indoor positioning and mapping applications relate very poorly to the actual mobile position, and topological environment information can be hardly inferred.

Looking one step ahead, of particular interest for next generation (5G) mobile wireless communications is the millimeter-wave (mmW) technology, for example that working in the 60 GHz band, which allows the exploitation of a large clean unlicensed bandwidth (up to 7 GHz) [START_REF] Rappaport | Millimeter Wave Wireless Communications[END_REF]- [START_REF] Boccardi | Five disruptive technology directions for 5G[END_REF]. The reduced wavelength at mmW paves the way for packing massive arrays with electronic steering capabilities into a small area, as done in [START_REF] Hong | Study and prototyping of practically largescale mmWave antenna systems for 5G cellular devices[END_REF]. Authors in [START_REF] Kaouach | Wideband low-loss linear and circular polarization transmit-arrays in V-Band[END_REF] studied how to realize mmW massive arrays, with dimensions which are becoming closer and closer similar to that of a tablet, even though several technological issues have still to be solved, such as energy consumption.

The idea to adopt mmW radars to overcome the short-comings of laser is not new and it has been considered for outdoor applications in some previous works, such as [START_REF] Dissanayake | A solution to the simultaneous localization and map building (SLAM) problem[END_REF]- [START_REF] Jose | Predicting Millimeter Wave Radar Spectra for Autonomous Navigation[END_REF], where a dedicated very high-directional antenna is mechanically steered. The consequent nearpencil beam returns a precise angle and range information thus making the modeling and characterization of the environment with mmW radars very similar to that based on laser. Obviously, such systems can not be integrated into smartphones or tablets. Stimulated by this possibility, in this paper we put forth the idea of placing massive antenna arrays in smartphones or tablets thus realizing a high-definition and low-cost personal mobile radar [START_REF] Guidi | Millimeter-wave massive arrays for indoor slam[END_REF]- [START_REF] Wu | Smartphones based crowdsourcing for indoor localization[END_REF]. In fact, the availability of a large number of antennas on mobile devices opens the possibility of realizing narrow beams with electronic-driven steering capabilities that are fundamental for high-accuracy mapping as a cheap, energy efficient alternative to laser technology (see Fig. 1) [START_REF] Guerra | Application of transmitarray antennas for indoor mapping at millimeter-waves[END_REF]. In our vision the pervasive availability of personal radar capabilities on smartphones will enable high-definition, zero-effort and automatic indoor building mapping by exploiting the possibility of crowd sourcing (crowd mapping) [START_REF] Alzantot | Crowdinside: automatic construction of indoor floorplans[END_REF]- [START_REF] Faragher | Opportunistic radio SLAM for indoor navigation using smartphone sensors[END_REF]. In addition, other new applications are possible such as blind people aid or the pictorial 3D environmental mapping where the third dimension results by the integration of the angle/range information collected by the personal radar with the 2D picture coming from the conventional embedded camera. From a complementary perspective, in [25], a system of stereo cameras and integrated sensors is used to have the depth information and to create a 3D image of the environment. The main drawback of this solution is that cameras are not able neither to automatically scan the environment as they do not perform a beamforming operation (i.e. the user must be participatory as he/she has to "pilot" the scan process by properly moving the camera) nor to retrieve information about the electromagnetic properties of materials. Our approach outperforms these weakness thanks to the adoption of mmW technology and of a multi-antenna radar system, which can scan the environment even if the smartphone is kept in the user pocket, i.e. without having the user participatory.

In our proposed mmW personal radar, wideband signals are employed to assure high ranging performance [START_REF] Guerra | Position and orientation error bound for wideband massive antenna arrays[END_REF]. Unfortunately, technology constraints prevent the use of wideband time-delay components necessary to realize the required signals alignment from different array elements and only simple digitally-controlled phase shifters can be realistically adopted. Therefore, loss in steering and ranging accuracy due to quantized phase shifters and frequency selectivity might arise [START_REF] Ruze | The effect of aperture errors on the antenna radiation pattern[END_REF], [START_REF] Rondinelli | Effects of random errors on the performance of antenna arrays of many elements[END_REF]. In addition, even in case massive arrays are adopted, laserlike beams are not obtainable and each measurement might be affected by contributions coming from different directions due to the effect of an increased main lobe width and side-lobes level in the antenna radiation pattern. Thus, dedicated mapping approaches are needed to overcome these limitations and to make massive arrays an attractive candidate for environment mapping. To validate our idea, the main contributions of the paper can be summarized as follows:

• We analyse the steering and ranging capabilities of the proposed mmW personal radar concept as a function of the signal bandwidth and the number of antennas composing the array. Specifically, the impact of frequency selectivity and quantization phase errors on the beam pointing error and beam width spread is jointly characterized. • We propose a low-complexity personal radar transceiver employing a massive antenna array and based on the transmission of wideband signals with non-coherent processing of the corresponding backscatter response of the environment. Our approach is different from the current state-of-theart. In fact, a classic two-step approach is typically followed [START_REF]An augmented state SLAM formulation for multiple lineof-sight features with millimetre wave radar[END_REF], [START_REF] Jose | Predicting Millimeter Wave Radar Spectra for Autonomous Navigation[END_REF]: first, for each test direction the obstacle is detected and its distance measured and, subsequently, this "hard" decision becomes the input for the mapping algorithm. Differently, we propose a new observation model that considers all the available raw measurements ("soft" decision) and hence fully exploits all the information coming from the antenna array. • The environment mapping capabilities of the personal radar are investigated by introducing a new Bayesian state-space model which considers a grid representation of the environment based on the electromagnetic characteristics of each space element and incorporating the new observation model able to properly account for the intrinsic characteristic of the massive antenna array and scanning process. The rest of the paper is organized as follows. The impact of errors on pointing and angle resolution in massive antenna arrays due to frequency selectivity and phase quantization noise is characterized in Sec. 2 and 3, where the trade-off between ranging and steering accuracy as a function of the number of elements composing the antenna array is also investigated. The personal radar mapping capabilities are modeled and characterized in Sec. 4, whereas in Sec. 5 numerical results related to a typical indoor scenario are provided. Finally, the conclusions are drawn in Sec. 6.

SOURCES OF ERROR IN MASSIVE ANTENNA ARRAYS

Planar Phased Arrays

The planar array configuration appears well suited for massive antennas array integration in smartphones and tablets. Consider the scheme reported in Fig. 2, where N array = M N antenna elements are placed in a rectangular area of dimension L x × L y , with L x = M d x and L y = N d y . The elements are spaced apart of d x (d y ) in the x (y) dimension. According to [START_REF] Balanis | Antenna Theory[END_REF], and considering a uniform array, the array factor (AF) at frequency f is (1)

with

Ψ x = kd x sin(θ) cos(φ) + β x Ψ y = kd y sin(θ) sin(φ) + β y ( 2
)
with Θ = (θ, φ) being the considered direction identified by the elevation and the azimuthal angles, θ and φ, respectively, and k = 2πf /c the wavenumber, with c indicating the speed of light. To steer the main lobe towards a specific direction Θ 0 , parameters β x and β y have to be set to

β x = -k 0 d x sin(θ 0 ) cos(φ 0 ) β y = -k 0 d y sin(θ 0 ) sin(φ 0 ) (3)
with k 0 = 2πf 0 /c, and f 0 being the reference frequency for which the steering parameters are designed. In the steering direction, for f = f 0 and in the absence of array non-idealities, it is AF(Θ 0 , f 0 ) = N array .

The complex beam-steering weights determining the beam shape and pointing direction can be expressed as

ω mn = e jϕmn = e j[(m-1)βx+(n-1)βy] (4)
where ϕ mn is the weight phase, and β x and β y are defined in (3).

Sources of Error in Antenna Arrays

There are several issues that arise when massive wideband phased arrays are adopted and both accurate beam steering and ranging have to be guaranteed. They are even more pronounced when working at very high frequencies due to technological and cost constraints [START_REF] Krieger | Dense delta-sigma phased arrays[END_REF]. This section deals with the effects of phase quantization errors and large signal bandwidth. The joint impact of such sources of error results in beam pointing errors and beamwidth spread that must be properly characterized when analyzing the performance of mapping schemes, as done in Sec. 4.

Large Signal Bandwidth

A large bandwidth is in general desirable thanks to the corresponding achievable high ranging resolution [START_REF] Dardari | Ranging with ultrawide bandwidth signals in multipath environments[END_REF]. However, in a wideband system, the received signal arrives at each antenna element with a delay that is not negligible compared to the signal duration, and hence it cannot be compensated by adopting only phase shifters, as typically done in narrowband systems. Unfortunately, the adoption of a huge number of time delay circuits represents a high-cost solution, especially at millimeter frequencies, then phase shifters remain at the moment 1. For simplicity and without loss of generality, we consider isotropic array elements. the most viable solution [START_REF] Baltus | Systems and Architectures for Very High Frequency Radio Links[END_REF]. Consequently, in the absence of time delays, the accumulation of the component received by each branch will result in a signal shape distortion. In other words, the actual steering direction and beamwidth would depend on frequency (beam squinting effect). For this reason the effects of signal bandwidth have to be carefully accounted for.

Array Weights Errors

Even the adoption of high resolution analog phase shifters might still not be a feasible solution for mmW massive antenna arrays if we aim at keeping a low system cost. Digitally controlled phase shifters implementing a discrete set of phase shifts represent a cheaper alternative at the expense of quantization errors [START_REF] Guerra | Application of transmitarray antennas for indoor mapping at millimeter-waves[END_REF]. Some new solutions have been investigated in [START_REF] Krieger | Dense delta-sigma phased arrays[END_REF] to find a compromise between the number of elements and the phase shifters accuracy.

Starting from the general model for the quantized version ωmn of weights ω mn proposed in [START_REF] Rondinelli | Effects of random errors on the performance of antenna arrays of many elements[END_REF], we express ωmn = ω mn e jδmn , with δ mn ∼ U(-∆/2, ∆/2) uniformly distributed random variables (RVs) in the quantization step ∆. Note that array weights errors affect the beam steering accuracy (see Fig. 2) and, consequently, the beamwidth characteristics. In addition, the misalignment of pulses at each array element causes temporal spread if it is not well counteracted. These effects can reduce the angle and ranging resolution as it will be detailed in the next section.

IMPACT OF ERRORS ON POINTING AND AN-GLE RESOLUTION

We now investigate the impact of errors on the array characteristics. In particular, denote with AF(Θ, f ) the actual AF in the presence of the aforementioned errors. Expression (1) becomes

AF(Θ, f ) = M m=1 N n=1
e j[(m-1)Ψx+(n-1)Ψy] e jδmn [START_REF] Morris | 3D indoor mapping for micro-uavs using hybrid range finders and multi-volume occupancy grids[END_REF] where δ mn 's are the phase errors previously defined. Now denote with Θ 0 = ( θ0 , φ0 ) the actual beam pointing direction at frequency f and in presence of arrays non idealities. We can relate the actual beam pointing error to the desired beam pointing direction as Θ 0 = Θ 0 +∆ Θ 0 +δΘ 0 , with ∆ Θ 0 +δΘ 0 indicating the shift from the desired pointing direction consisting of both a deterministic ∆ Θ 0 and random quantity δΘ 0 = (δθ 0 , δφ 0 ) due to frequency selectivity and phase quantization, respectively (Fig. 2). Note that Θ 0 = ( θ0 , φ0 ) = Θ 0 + ∆ Θ 0 is the actual beam pointing direction at frequency f and in absence of quantization error, with

θ0 = sin -1 f 0 f sin(θ 0 ) , φ0 = φ 0 . (6)
For what follows it is convenient to derive the squared AF

| AF(Θ, f )| 2 = mnpq e j[(m-p)Ψx+(n-q)Ψy] e j(δmn-δpq) = |AF(Θ, f )| 2 + ζ err (Θ, f ) (7)
with

ζ err (Θ, f ) = mnpq e j[(m-p)Ψx+(n-q)
Ψy] e j(δmn-δpq) -1 ( 8)

and mnpq = M m N n M p N q
. Equation ( 7) puts in evidence the deterministic and random terms composing the squared AF.

Considering δ mn 's as independent, identically distributed (i.i.d.) RVs, the average squared AF is

E | AF(Θ, f )| 2 = N array + |AF(Θ, f )| 2 -N array χ 2 (9)
where χ = E e jδmn = 2 ∆ sin ∆ 2 . In the following we will discuss about the two effects that arise when phase errors and frequency selectivity are present. First, the beam pointing error phenomenon is investigated. Second, we address the beamwidth deviations caused by the joint influence of phase errors and squint effect.

Beam Pointing Error

We now evaluate the impact of phase errors on the beam pointing direction. In particular we derive the mean and the variance of the mean beam pointing error.

The mean square error (MSE) of the beam pointing error is given by

MSE θ (f, θ 0 ) = E δθ 0 +∆ θ 0 2 = ∆ θ 0 2 + E δθ 2 0 ( Θ 0 , f ) MSE φ (f, φ 0 ) = E δφ 0 +∆ φ 0 2 = ∆ φ 0 2 + E δφ 2 0 ( Θ 0 , f ) (10)
with

∆ φ 0 2 =0, ∆ θ 0 2 = sin -1 f 0 f sin(θ 0 ) -θ 0 2 . (11)
In Appendix A it is shown that E [δθ 0 ] = E [δφ 0 ] = 0, and in Appendix B explicit expressions for E δθ 2 0 ( Θ 0 , f ) and E δφ 2 0 ( Θ 0 , f ) are derived and expressed by (58).

Given the transmitted signal power spectral density (PSD) P t (f ), with bandwidth W , we finally define the pointing root mean square error (PRMSE) as

PRMSE(Θ 0 ) = W P t (f ) [MSE θ (f, θ 0 )+MSE φ (f, φ 0 )]df W P t (f )df . ( 12
)
It will be shown in Sec. 5 that small bandwidths W are better suited for high beam pointing accuracy. This is in contrast with the need to have large bandwidth signals for accurate ranging and a trade-off has to be found as will be reported in the numerical results.

Beamwidth Spread

We are now interested in finding the 3 dB beamwidth spread Θ 3dB = θ0 ± δθ, φ0 ± δφ , which varies with the frequency even if this is not explicitly indicated for notation convenience. In particular, the half power beamwidth (HPBW) can be evaluated by solving the following equation

|AF( Θ 3dB , f )| 2 = 1 2 |AF( Θ 0 , f )| 2 . (13)
It follows that mnpq ŵmnpq e j(m-p)Ψ δx e j(n-q)Ψ δy e j(δmn-δpq) = 1 2

|AF( Θ 0 ,f )| 2 (14)
where ŵmnpq = e j (m-p)Ψx| Θ= Θ 0 +(n-q)Ψy| Θ= Θ 0 and

Ψ δx = kd x cos( θ0 ) cos( φ0 )δθ -sin( θ0 ) sin( φ0 )δφ Ψ δy = kd y sin( θ0 ) cos( φ0 )δφ + cos( θ0 ) sin( φ0 )δθ . (15)
In Appendix C a practical expression for the numerical evaluation of ( 13) is derived. From ( 62) and ( 63) in Appendix C we introduce the beam solid area (BSA)

BSA = Ω A (Θ 0 , f ) = 4 E δ θx0 ( Θ 0 , f ) E δ θy0 ( Θ 0 , f ) sec θ0 . (16)
In addition, we define the effective BSA as

Ω A, eff (Θ 0 ) = W P t (f )Ω 2 A (Θ 0 , f )df W P t (f )df . ( 17
)
This parameter is an indicator of the angle resolution of the scanning process at direction Θ 0 . It will be shown in Sec. 5 through simulations that a large number of elements lets to reduce the impact of quantization errors as well as to narrow the beamwidth.s

Ranging and Steering Accuracy

We now characterize the ranging and angle resolution of massive antenna arrays accounting for bandwidth and sources of errors discussed in Sec. 2.

Target distance estimate (ranging) during each step of the scanning process is obtained through the measurement of the signal round-trip time (RTT). The fundamental limit of ranging accuracy of a target at distance d depends on the signal effective bandwidth W eff and the received signal-to-noise ratio (SNR) as given by the Cramér-Rao lower bound (CRLB) [START_REF] Dardari | Ranging with ultrawide bandwidth signals in multipath environments[END_REF]

CRLB = c 2 8π 2 SNR W 2 eff (18)
where W eff is defined as

W eff +∞ -∞ f 2 P t (f -f 0 ) df +∞ -∞ P t (f -f 0 ) df ≈ W . ( 19
)
The mmW radar range equation related to a target at distance d, a test steering direction Θ 0 , and considering a waveform with incident direction Θ i = (θ i , φ i ) with respect to the target, can be expressed as

P r (Θ 0 , Θ i , d) = W P t (f )c 2 G 2 (Θ 0 , f ) M (Θ 0 , Θ i , f, d) f 2 (4π) 3 d 4 df (20)
with G(Θ 0 , f ) ≈ π 2 /Ω A (Θ 0 , f ), being the array gain [START_REF] Balanis | Antenna Theory[END_REF] and M (Θ 0 , Θ i , f, d) the target radar cross section (RCS). Under the narrow beamwidth hypothesis, i.e. high N array , we approximate the illuminated area on the target surface with the beam area of the incident wave generated by the massive array (footprint area). Thus, the RCS of a target at mmW is given by [START_REF] Wagner | Radiometric calibration of small-footprint fullwaveform airborne laser scanner measurements: Basic physical concepts[END_REF], [START_REF] Sarabandi | Modeling and measurements of scattering from road surfaces at millimeter-wave frequencies[END_REF]

M (Θ 0 , Θ i , f, d) = 4πρA fp (Θ 0 , f, d)ζ(Θ i ) ( 21
)
where ρ is a parameter which depends on the type of scattering and the target material, and ζ(Θ i ) is the incident/reflection coefficient of the target in the direction Θ i . The footprint area is

A fp (Θ, f, d) = π 4 d 2 Ω A (Θ, d) (22)
with Ω A (Θ, f ) given by [START_REF] Jose | Predicting Millimeter Wave Radar Spectra for Autonomous Navigation[END_REF]. Thus [START_REF] Guerra | Application of transmitarray antennas for indoor mapping at millimeter-waves[END_REF] results

P r (Θ 0 , Θ i , d) = πc 2 64 W EIRP(f ) ρ ζ(Θ i ) d 2 f 2 df ( 23
)
where EIRP(f ) = P t (f ) G(Θ 0 , f ) is the effective radiated isotropic power (EIRP) typically constrained by regulatory power emission limits. Note the quadratic dependency of P r (Θ 0 , Θ i , d) on the mobile-target distance d.

The actual SNR can be improved by increasing the transmission power or the observation time (e.g., by transmitting several pulses/symbols). However, both parameters have some constraints on spectrum emission regulations and the maximum tolerable scanning time, respectively. For what the latter is regarded, suppose we have to scan half of the full solid angle or half of the full angle, i.e. Ω o = 2π [sr] and

Ω o = π [rad] in 3D and
2D scenarios, respectively. The number N steer of steering test directions to complete the scan process is strictly related to the array beamwidth, and thus to the number of antenna elements N array , by the relation 2

N steer = Ω o Ω mean ( 24
)
where Ω mean denotes the mean BSA HPBW averaged over all the test directions Θ 0 given by

Ω mean = 1 Ω o Ωo Ω A,eff (Θ) dΩ . ( 25
)
Note that Ω mean provides an indication on the scanning average angle resolution. As a consequence, the scanning time is T scan = T ob N steer , with T ob representing the observation time for each test direction. Note the dependency on the BSA: a narrow beam increases the scanning time, while the scanning angle resolution is improved. Under scanning time constraints, the SNR is given by

SNR = P r (Θ 0 , Θ i , d)T ob N 0 = P r (Θ 0 , Θ i , d)T scan N steer N 0 (26)
where

N 0 = N array κ • T 0 • F , being κ the Boltzmann constant, T 0 = 290 [K]
the receiver temperature, and F the receiver noise figure. Eq. ( 26) can be used in [START_REF] Guerra | Millimeter-wave personal radars for 3D environment mapping[END_REF] to derive the fundamental limit of ranging as a function of the average BSA and scanning time constraint as it will be investigated in the numerical results.

PERSONAL RADAR MAPPING CAPABILITY

The previous analysis has highlighted how angle resolution, scanning time, signal bandwidth and ranging accuracy are tightly coupled.

We investigate now the mapping capabilities of the personal radar integrating massive arrays working at mmW. We consider a person with a mobile terminal (e.g., a smartphone or a tablet) equipped with a mmW massive antenna array who walks with constant speed along a pre-determined trajectory and whose position and orientation are constantly determined and updated by sensors integrated on the terminal or by external localization technologies [START_REF] Dardari | Indoor tracking: Theory, methods, and technologies[END_REF]. The main task of personal radar is to incrementally build a consistent 2D/3D map of the surrounding environment. To this purpose, first we propose a mmW radar collecting a rich set of measurements from the massive antenna array. Then a grid-based probabilistic mapping approach is adopted to fully exploit such measurements in a static environment. Hereafter, for simplicity of notation, we will restrict our analysis to a 2D scenario even if the 3D extension is straightforward.

2. ⌊x⌋ denotes the biggest integer smaller than x.

Massive Antenna Array Radar

The purpose of the radar is to provide a rich set of observations for the subsequent mapping phase starting from the analysis of the signals backscattered by the surrounding environment.

Let

θ b = -π/2 Nsteer-1
Nsteer + (b -1)π/N steer , b = 1, 2, . . . N steer , be the set of N steer steering directions considered during the scanning process defined, and θ v the mobile orientation.

Consider a generic interrogation signal, for each steering direction θ b , composed of N p wideband pulses p(t)

g(t) = Np-1 i=0 p(t -iT f ) (27)
with T f being the time frame chosen so that all signals backscattered by the environment are received by the receiver before the transmission of the successive pulse, thus avoiding inter-frame interference. The consequent scanning time is T scan = T ob N steer , where T ob = N p T f . Since each pulse is backscattered by the surrounding targets populating the environment, for the steering direction θ b the received signal can be expressed as

r(t, θ b ) = Np-1 i=0 s(t -iT f , θ b ) + n(t) (28)
where s(t, θ b ) = p(t)⊗h C (t, θ b ) is the channel response to the transmitted pulse p(t) at direction θ b , with h C (t, θ b ) being the channel impulse response (CIR), de-embedded of the propagation time, taking the joint environment and antenna array response into account, 3 and with n(t) being the additive white Gaussian noise (AWGN) with two-sided power spectral density N 0 /2. The received signal is first passed through an ideal bandpass filter with center frequency f 0 to eliminate outof-band noise. 4 The filtered signal is denoted by

y(t, θ b ) = Np-1 i=0 x(t -iT f , θ b ) + z(t) (29)
where

x(t, θ b ) = s(t, θ b )⊗h F (t) ∼ = s(t, θ b ) and z(t) = n(t)⊗ h F (t), with h F (t)
being the impulse response of the filter.

To conjugate the need of having a manageable number of measurements and a low complexity receiver, we consider a non-coherent approach based on energy measurements in a discretized time scale. This approach also accounts for the complete uncertainty on the received waveform shape deriving from the a priori ignorance on environment electromagnetic characteristics. Specifically, energy measurements are taken during the time frame T f after the transmission of each pulse by subdividing the time frame into N bin = ⌊T f /T ED ⌋ time slots (bins) of duration T ED . Energy measurements are accumulated for each time bin over the N p frames of the interrogation 3. ⊗ is the convolutional operator. 4. This operation is necessary since the receiver is energy-based, as it will be described later.

signal. The accumulated energy measurement at the sth time bin and bth steering angle θ b is

e bs = Np-1 k=0 s TED (s-1) TED y 2 (t + kT f , θ b ) dt (30)
with s = 1, 2, . . . , N bin and b = 1, 2, . . . N steer . Note that T ED must be chosen to accommodate most of the energy of the received pulse, i.e. T ED ≈ 1/W .

According to [START_REF] Urkowitz | Energy detection of unknown deterministic signals[END_REF], for each energy bin, the normalized energy measurement output can be well approximated by

Λ bs = 2 N 0 e bs ≃ 1 σ 2 Np-1 k=0 sN d i=(s-1)N d x i (θ b ) + z (k) i 2 (31)
where

N d = 2W T ED , σ 2 = N 0 W is the noise variance, and
z (k) i
are for odd i (even i) the samples of the real (imaginary) part of the equivalent low-pass of z(t + kT f ), k = 1, 2, . . . , N p , taken at Nyquist rate W in each interval T ED . In (31) we used the property x(t+kT f , θ b ) = x(t, θ b ), with k = 1, 2, . . . , N p , so that x i (θ b ) represents for odd i (even i) the samples of the real (imaginary) part of the equivalent low-pass of x(t, θ b ), taken at Nyquist rate W in each interval T ED , which does not depend on k.

It turns out that Λ bs is a non-central Chi-square distributed RV with N = N p N d degrees of freedom. The non-centrality parameter is λ bs = 2γ bs [START_REF] Urkowitz | Energy detection of unknown deterministic signals[END_REF] , where γ bs is given by

γ bs = N p 2σ 2 sN d i=(s-1)N d x 2 i (θ b ) ≃ N p N 0 s TED (s-1) TED x 2 (t, θ b ) dt . ( 32
)
Note that [START_REF] Baltus | Systems and Architectures for Very High Frequency Radio Links[END_REF] represents the accumulated SNR corresponding to the sth time bin and the bth steering direction which increases with the number N p of pulses at the expense of a longer observation time T ob . Furthermore, each energy measurement e bs does not depend only on the signal backscattered in the steering direction θ b , but also on all contributions coming from different directions according to the shape of the array radiation pattern which affects h C (t, θ b ) and hence x(t, θ b ). During the mapping process these contributions could create ambiguities but, at the same time, permit to discriminate between different array configurations.

Probabilistic State-Space Model for Mapping

In the following, we introduce a state-space model for probabilistic mapping able to exploit the measurements from the massive array radar under a Bayesian filtering approach [START_REF] Guerra | Millimeter-wave personal radars for 3D environment mapping[END_REF]. In a probabilistic mapping framework, the position and orientation of the mobile is known or partially known if resulting from an external localization system exploiting, for example, inertial measurement unit (IMU) and RF sensors measurements. The same mmW technology could be in principle exploited to provide extremely high accuracy positioning down to a few cm when interacting with base stations [START_REF] Guerra | Position and orientation error bound for wideband massive antenna arrays[END_REF].

Different models can be adopted to represent the map. Feature-based models are privileged when a set of landmarks characterizes the environment and the sensors on the mobile device are able to well-discriminate them (i.e. using laser and video technologies). On the contrary, grid-based models for map representation are more suitable when sensors provide imprecise measurements from which it is difficult to discriminate between different features, such as when using radar [START_REF] Wurm | Bridging the gap between feature-and grid-based SLAM[END_REF]. In gridbased approaches, the environment is discretized in cells whose position is a-priori known and hence the number of states is fixed to the number of cells. In most of works, the state is a binary RV and represents the presence or not of any physical element in the corresponding cell (occupancy grid) [START_REF] Ozisik | An occupancy grid based SLAM method[END_REF].

In our work we consider RCS-based description of cells in a map-grid representation context, and the adoption of massive array radars with electronically steering beams instead of mechanically steering devices. In particular, a 2D grid of N L = X grid Y grid cells representative of the environment is considered. Differently from previous works, we propose a one-step procedure where a larger observation vector is introduced accounting for the specific massive antenna array radiation pattern, investigated in Secs. 2 and 3, and whose response could depend, in principle, on several cells, as illustrated in Fig. 3. Therefore, the mapping process works with a richer information set and no specific radar detection scheme and associated detection threshold evaluation strategy have to be designed, as required instead in conventional two-step approaches [START_REF]An augmented state SLAM formulation for multiple lineof-sight features with millimetre wave radar[END_REF]. In the following, the adopted scheme is described, and the differences with classic Bayesian SLAM are highlighted.

State Vector

We define the state vector at time k as

x(k) = m(k) = [m 1 (k), . . . , m i (k), . . . , m N L (k)] T (33)
where m i (k) indicates the root radar cross section (RRCS) (with sign) of the ith cell of the grid. Obviously, if the cell is empty (only air), the corresponding RRCS is zero. Note that in the model herein adopted we neglect the dependency of m i (k) from the frequency and the environment is assumed stationary. T [START_REF] Sarabandi | Modeling and measurements of scattering from road surfaces at millimeter-wave frequencies[END_REF] where e bs (k) is given by [START_REF] Krieger | Dense delta-sigma phased arrays[END_REF]. As shown in the previous section, e bs (k) is proportional to a non-central Chisquare distributed RV that makes the formulation of the observation model challenging. However, for large N (typically > 50), i.e. for large N p , the non-central Chi-square distribution can be approximated with a Gaussian distribution having mean λ bs +N and variance 2(N + 2λ bs ) [START_REF] Chun | Normal versus noncentral chi-square asymptotics of misspecified models[END_REF]. We exploit this property to construct a tractable Gaussian observation model of the state-space model

z(k) = [z 11 (k), . . . , z bs (k), . . . , z NsteerN bin (k)] T = h(x(k)) + v(k) (35)
where z bs (k) is the energy observation corresponding to the sth time bin and the bth steering angle, h(•) is the observation function that relates the state vector to the observations, and

v(k) = [v 11 (k), . . . , v bs (k), . . . , v NsteerN bin (k)] T ( 36
)
is the vector of uncorrelated Gaussian observation errors with zero mean and covariance matrix R(k) = diag σ 2 11 (k), . . . , σ 2 bs (k), . . . , σ 2 NsteerN bin (k) with σ 2 bs (k) = Var (e bs (k)). Note that the length N m = N steer N bin of z(k) is related either to the number of the steering directions, depending on the antenna array pattern G(θ, f ), and to the number of time bins depending on the signal bandwidth W through T ED .

The generic element of z(k) can be written as

z bs (k) = h bs (x(k)) + v bs (k) = E [e bs (k)] + v bs (k) ( 37
)
with b = 1, 2, . . . , N steer , s = 1, 2, . . . , N bin . From [START_REF] Dardari | Ranging with ultrawide bandwidth signals in multipath environments[END_REF] and [START_REF] Baltus | Systems and Architectures for Very High Frequency Radio Links[END_REF] it is E [e bs (k)] = T f P bs (k) whereas v bs (k) is the measurement noise, as defined in [START_REF] Urkowitz | Energy detection of unknown deterministic signals[END_REF], with variance σ 2 bs (k) = N 0 (N p W T ED N 0 + 2T f P bs (k)). Using (21), the term P bs (k) depends on x(k) and represents the average received power at the bth steering angle and the sth time bin given by

P bs (k) = E [e bs (k)] T f = P w + i∈R(s) W P t (f ) c 2 m 2 i (k) G 2 (θ i -θ v -θ b , f )N p f 2 (4π) 3 d 4 i df ( 38
)
where

P w = σ 2 T ED Np T f
is the noise power, d i and θ i are, respectively, the distance and angle between the mobile and the ith cell. Function R(s) returns the set of cells indexes contributing to the sth time bin (i.e. all the cells located at a distance close to s • d bin with

d bin = c • T ED /2).
Note that in [START_REF] Ozisik | An occupancy grid based SLAM method[END_REF] the actual footprint area corresponds to the aggregate area of the cells illuminated by the beam. In fact, m i (k) is evaluated using [START_REF] Alzantot | Crowdinside: automatic construction of indoor floorplans[END_REF], by setting the footprint area A fp to the area of a cell.

EKF-Mapping

There are several recursive approaches belonging to the Bayesian filtering theory to evaluate the posterior distribution p(x(k)|z(1 : k)) of x(k) given the set of measurements z(1 : k) = {z(1), z(2), . . . , z(k)} collected by the mobile device until time instant k. Among them, here we employ the extended Kalman-Filter (EKF) which offers a good trade-off between performance and complexity considering the intrinsic large state vector of dimension N L resulting from the grid-based map representation.

The mean and covariance state matrix of the predicted state are, in linear, x(k + 1|k) = x(k|k) and P(k + 1|k) = P(k|k) respectively, due to the stationarity of the environment. 5 The prediction for the observation model is given by ẑ(k + 1|k) = h(x(k + 1|k)). Once the new measurement vector e(k + 1) is available, the innovation in the model is

ν(k + 1) = e(k + 1) -ẑ(k + 1|k) . ( 39
)
Consequently, the correction step performed by the EKF is described by

x(k + 1|k + 1) = x(k + 1|k) + W(k + 1) ν(k + 1)

P(k+1|k+1) = P(k + 1|k) -W(k + 1)S(k + 1) W T (k+1) (40)
where

S(k + 1) = ∇h P(k + 1|k)∇h T + R(k + 1) W(k + 1) = P(k + 1|k)∇h T S -1 (k + 1) (41)
with ∇h being the Jacobian of h(•) evaluated in x(k + 1|k), which gives using ( 37) and ( 38)

∂h bs ∂m i = 2 m i W Pt(f ) c 2 G 2 (θi-θv-θ b ,f )T f Np f 2 (4π) 3 d 4 i df if i ∈ R(s) 0 otherwise. ( 42
)
In the numerical results a case study is provided implementing the state-space model and the EKF mapping approach here introduced.

NUMERICAL RESULTS

We first evaluate the impact of frequency selectivity and quantization errors on antenna array characteristics in order to investigate the ranging and steering accuracy as a function of N array . Successively, mapping capability is analyzed by exploiting the new state-space model proposed in Sec. 4.

System Parameters

According to the wavelength λ 0 and the object size L, there are different scattering regimes. Coherently with [START_REF] Sarabandi | Modeling and measurements of scattering from road surfaces at millimeter-wave frequencies[END_REF], we assume here to have λ 0 ≪ L (optics regime) even if in some cases it could result in a rough approximation. This assumption is reasonable in indoor environments where most of targets are constituted by walls that have an extension much larger than the wavelength. We consider the approximation of a Lambertian scattering [START_REF] Degli-Esposti | Measurement and modelling of scattering from buildings[END_REF], adopting the laser model proposed in [START_REF] Wagner | Radiometric calibration of small-footprint fullwaveform airborne laser scanner measurements: Basic physical concepts[END_REF] for the monostatic laser RCS, where ρ is defined as ρ = ρ d /π, with ρ d being the diffuse reflectivity (albedo), and Θ i normal to the target. In addition, we assume a target made of aerate concrete by approximating the albedo with the power reflection coefficient, evaluated by setting the concrete relative permittivity to ǫ r = 2.26 and its loss tangent to 0.0491 [START_REF] Correia | Estimation of materials characteristics from power measurements at 60 GHz[END_REF]. We also considered a scanning time T scan = 0.1 ms, F = 4 dB and the transmission of a root raised cosine (RRC) signal compliant with the Federal Communications Commission (FCC) mask at 60 GHz [START_REF]Revision of Part 15 of the Commissions Rules Regarding Operation in the 57-64 GHz Band[END_REF], but with an average EIRP limited to 30 dBm to make the system suitable for battery-powered handset devices.

Ranging and Steering Accuracy Trade-off

In this section, an example of trade-off between ranging and angle resolution is given considering a planar array working at 60 GHz with d x = d y = λ 0 /2 and M = N . The RMS CRLB on ranging is reported in Fig. 4, obtained for different mobile-target distances and for different values of array elements N array , where the transmitted signal bandwidth is set to W = 500 MHz. When N array increases, the SNR decreases due to the reduced BSA and the constraint on EIRP according to [START_REF] Guerra | Position and orientation error bound for wideband massive antenna arrays[END_REF]. 6 In fact, N array enters into N 0 and the array gain, which is inversely proportional to Ω A (Θ, f ). Thus, large 6. Note that here performance is strictly related to the adopted model, which is valid only for a high number of antenna elements, which justify the performance degradation when Narray increases. values of N array improve the scanning angle resolution, but worsen the SNR affecting [START_REF] Guerra | Millimeter-wave personal radars for 3D environment mapping[END_REF]. As an example, if we set a target ranging error of 5 cm at 5 m distance, from Fig. 4 it must be N array ≤ 100. On the contrary, the bandwidth improves the ranging accuracy, as shown in Fig. 5, but it has a detrimental effect on the PRMSE, as can be noticed in Fig. 6. From the same figure it can be seen that PRMSE is quite sensitive to the test steering elevation angle θ 0 and to phase quantization errors. Finally, Fig. 7 reports the BSA, evaluated according to (62) and (63), whose expectation is computed through Monte Carlo simulation. Note the weak dependance of the BSA on the frequency, and how large N array results in lower BSA and lower sensitivity to quantization errors. In any case the effect of quantization errors on BSA can be neglected with good approximation. Moreover, if we set BSA ≤ 0.04 [sr], it must be N array ≥ 64.

Thus, we assist to a strong impact of N array on the BSA, which affects both the SNR and the ranging accuracy, and hence it needs to be properly taken into account during the system design.

Furthermore, there is also the trade-off in the choice of the bandwidth: the ranging accuracy improves with a larger bandwidth but the array pointing error increases with W due to frequency selectivity effects. The bandwidth and the phase errors slightly affect the effective BSA and hence the scanning angle resolution.

In our numerical example, N array = 100 and W = 1 GHz represent a reasonable compromise.

Mapping Performance

We implemented the previously described mapping algorithm in MATLAB in order to validate its theoretical performance.

The case study proposed refers to a typical indoor office environment of 20 × 20 m 2 , whose layout is shown in Fig. 8, which has been discretized in a grid of cells having area 0.2 × 0.2 m 2 each, according to the high arrays resolution, and accounting for a wall attenuation, due to its thickness, of 40 dB [START_REF] Correia | Estimation of materials characteristics from power measurements at 60 GHz[END_REF]. The mobile device with personal radar capabilities moves along a pre-determined trajectory and speed of 1m/s, and its reference orientation is always towards the direction of movement. Radar measurements are taken every second, and black circles refer to the true radar positions while the green squares to the positions perceived by the mobile device.

We compare the personal radar performance in three cases, N array = 15×15 , N array = 10×10 and N array = 4×4, respectively. In the previous numerical results, the array with N array = 10 × 10 elements resulted as a value corresponding to a good compromise between ranging and steering accuracy. A time frame of T f = 100 ns is considered with bandwidth W = 1 GHz and, consequently, T ED = 1 ns if not otherwise indicated. The constraint on the scanning time is fixed to T scan = 80µs while the number of steering directions and of pulses per interrogation signal varies with the number of antenna elements, as N steer increases with N array . In fact, for N array = 16, N array = 100 and N array = 225 antenna elements, it is N p = 100, N p = 40 and N p = 27, respectively.

In the following results, the reconstructed map is represented in terms of RCS estimated through the EKF algorithm and it is juxtaposed to the true map of the office environment. Obviously, only the walls resulting in line-of-sight during the mobile movement are expected to be identified.

Performance in Absence of System Non-Idealities

Fig. 8 shows the estimated maps obtained with ideal (i.e. quantization error-free) antenna arrays of 4 × 4 (left), 10 × 10 (right) and 15 × 15 (bottom) elements, and with perfect knowledge of the mobile position and orientation. The higher number of cells that are correctly identified using the 10 × 10 array with respect to 4 × 4 and 15 × 15 array puts in evidence that an enhanced performance in terms of map estimation accuracy can be found as a compromise between ranging accuracy and angular resolution, in agreement with the results of Sec. 5.2. In fact, for N array = 15 × 15, the map is well defined at the prize of a reduced SNR, and thus part of the map is not reconstructed, while for N array = 4 × 4 the higher number of pulses improves the SNR enabling the full reconstruction of the map, but with a lower angle resolution. mapping performance improves due to the increased time resolution at the expense of an increased number N bin of bins, and thus of receiver complexity.

Performance with System Non-Idealities

The effect of phase quantization errors is shown in Fig. 10 where antenna arrays with very rough phasing capability of step ∆ = π/2 are considered. When compared to Fig. 8, it can be noted that, even if quantization errors might affect the array performance, their negative effect on map reconstruction can be significantly mitigated if their impact on radiation pattern is properly characterized in the measurement model ( 37) by adopting the analysis of Sec. 3. In Fig. 11 we include also mobile orientation and position errors caused by an imperfect external positioning system. In particular, we have taken into account orientation estimation errors with standard deviation of 5 • , and position estimation errors with 0.1 m standard deviation [START_REF] Xiong | Arraytrack: A fine-grained indoor location system[END_REF], both along the xand the y-axis. It is evident from the figure how estimated maps become much more degraded even if the performance of 10 × 10 antenna array is still acceptable.

Despite future indoor location systems are expected to guarantee a self-localization with errors below 0.5 m, current technologies provide errors typically greater than 0.5 m. Thus, in Fig. 12 we considered a mobile position error of 0.5 and 1 m, respectively. Performance is evidently degraded, but it could be improved, for instance, through crowd sensing mechanism where users share their measurements in a cooperative way. Fig. 13 reports the mapping results considering an augmented orientation estimation error of 10 • and 20 • , respectively. As we can see, the orientation estimation accuracy has a less detrimental effect on the mapping reconstruction with respect to the position.

Finally, in Fig. 14 we reported the results obtained when high position and orientation errors are present. As shown, the mapping quality is highly degraded and a crowd-sensing mechanism becomes essential to refine the map estimation.

CONCLUSIONS

We put forth the idea of adopting massive antenna arrays at mmW as an effective way to embed personal radar capabilities in future smartphones. The quantization and frequency selectivity effects have been accounted for to derive the pointing error, the beamwidth spread and the ranging accuracy of the personal radar as a function of the number of antenna elements in the array. Our analysis has highlighted how angular resolution, scanning time, signal bandwidth and ranging accuracy are tightly coupled and how a suitable trade-off has to be obtained. Successively, the environment mapping capabilities of the proposed personal radar have been investigated through the introduction of a non-coherent radar architecture and a new state-space model, within a grid-based Bayesian framework, in which each cell of the grid has been characterized by a specific RCS-based occupancy value. The peculiarities of the antenna array have been properly exploited by a direct integration of massive array radar measurements with the radar and mapping processes.

To assess the feasibility of the system, simulations have been performed in a realistic indoor scenario. Results have shown that an increased number of array elements guarantees more robustness to non-ideal effects and allows an improved performance in terms of map reconstruction, but it might suffer from a reduced SNR when a constraint on the scanning time is given. Thus, a tradeoff has to be found. In addition, it has been assessed that map reconstruction accuracy is more sensitive to mobile position and orientation perception errors than antenna quantization errors. The latter leads to important implications for the feasibility of indoor mapping with low-cost massive array antennas embedded in nextgeneration smartphones and toward the diffusion of crowd sourcing indoor mapping applications.

APPENDIX A

To evaluate E [δΘ 0 ], consider that in the presence of nonideal effects, the actual pointing direction be derived by solving

S( Θ 0 , f ) = ∂| AF(Θ, f )| 2 ∂θ Θ= Θ0 = 0 F ( Θ 0 , f ) = ∂| AF(Θ, f )| 2 ∂φ Θ= Θ0 = 0 . ( 43
)
Considering the expansion with the first order of the Taylor series evaluated in Θ 0 , it is S( Θ 0 ,f ) = S( Θ 0 + δΘ 0 , f ) = S( Θ 0 , f )+δθ 0 ∂S(Θ, f ) ∂θ

Θ= Θ0 = 0 F ( Θ 0 ,f )= F ( Θ 0 + δΘ 0 ,f ) = F ( Θ 0 , f )+δφ 0 ∂F (Θ, f ) ∂φ Θ= Θ0 = 0 . (44)
Thus, considering δΘ 0 = (δθ 0 ( Θ 0 , f ), δφ 0 ( Θ 0 , f )), and the approach followed in [START_REF] Rondinelli | Effects of random errors on the performance of antenna arrays of many elements[END_REF], where the derivatives are approximated with their respective mean values, we have from which E S( Θ 0 , f ) = 0. Analogously, we can make the same considerations for φ, obtaining E δθ 0 ( Θ 0 , f ) = E δφ 0 ( Θ 0 , f ) = 0, and consequently E [δΘ 0 ] = 0.

APPENDIX B

We express E δθ 2 0 ( Θ 0 , f ) and E δφ 2 0 ( Θ 0 , f ) as

E δθ 2 0 ( Θ 0 , f ) = E S 2 (Θ, f ) E ∂S(Θ,f ) ∂θ 2 Θ= Θ0 E δφ 2 0 ( Θ 0 , f ) = E F 2 (Θ, f ) E ∂F (Θ,f ) ∂φ 2 Θ= Θ0 . (51)
For what the denominator of (51) is regarded, by taking the derivatives of (47), we have: where

Ω mnpq = Λ mnpq ( Θ 0 , f )+j ψ 2 mnpq ( Θ 0 , f ) . ( 55
)
Regarding the numerator of (51) it is S 2 ( Θ 0 , f ) = mnpq rstu ψ mnpq ( Θ 0 , f )ψ rstu ( Θ 0 , f )e j(δmn-δpq) e j(δrs-δtu)

APPENDIX C

Equation ( 14) can be rearranged by developing the second-order series expansion of the term e j[(m-p)Ψ δx +(n-q)Ψ δy ] as where ∆ mnpq = (m -p)Ψ δx + (n -q)Ψ δy , and δ mnpq = δ mn -δ pq . To derive the BSA we adopt an approach similar to that proposed in [START_REF] Balanis | Antenna Theory[END_REF] and [START_REF] Elliot | Antenna Theory and Design[END_REF] where the beam is considered alternatively lying on the XZand YZplane. Specifically, consider the case in which φ0 = 0 and δφ = 0 so that the beam lies on the XZ-plane, with deviation ±δ θx0 from the beam pointing direction in that plane. We have that ∆ mnpq = (m -p)kd x δ θx0 cos( θ0 ) and (60) can be written as

and analogously

E δ θy0 ( Θ 0 , f ) ≈ E   ± |AF( Θ 0 , f )| 2
mnpq ŵmnpq e jδmnpq (n -q) 2 k 2 d 2 x cos 2 ( θ0 )

1 /2   . (63)
Solving ( 62) and (63) appears very complicated because the expected value regards the ratio of dependent RVs. Nonetheless, they are useful to drastically facilitate numerical calculation of the solution, as it will be done in Sec. 5.2.

  (m-1)Ψx+(n-1)Ψy]

Fig. 2 :

 2 Fig. 2: Array structure (a) and beam pointing error scheme (b).

Fig. 3 :

 3 Fig. 3: Grid-based mapping performed with the personal radar in three different positions at the instants k -1, k and k + 1.

4. 2 . 2

 22 Observation Model Define e(k) the vector containing the accumulated measured energy at the output of the receiver at time k e(k) = [e 11 (k), . . . , e 1N bin (k), . . . , e bs (k), . . . , e Nsteer1 (k), . . . e NsteerN bin (k)]

5 .

 5 The notation (m|n) indicates at time m given the observations until time n.

N

  Fig. 4: RMS CRLB vs. mobile-target distance for different values of Narray, Θ0 = (45 • , 45 • ), and W = 500 MHz.

NFig. 5 :

 5 Fig. 5: RMS CRLB vs. transmitted signal bandwidth W for a mobile-target distance of 10 m, and for Θ0 = (45 • , 45 • ).

2 Fig. 6 :

 26 Fig. 6: PRMSE vs. bandwidth for different values of Narray, different steering angles and in presence/absence of quantization effects.

N array = 36 Fig. 7 :

 367 Fig. 7: BSA vs.frequency for different values of Narray and quantization errors.

Fig. 8 :

 8 Fig. 8: 2D mapping using 4 × 4 (top-left), 10 × 10 (top-right) and 15 × 15 (bottom) antenna arrays with a perfect knowledge of mobile position and no quantization errors.

Fig. 9

 9 reports the maps estimated through the adoption of an ideal antenna array equipped with a 10 × 10 antenna array for different signal bandwidths. In particular, from the left to the right, the results obtained with W = 500 MHz (i.e. T ED = 2 ns), W = 2 GHz (i.e. T ED = 0.5 ns) and W = 3 GHz (i.e. T ED = 0.33 ns) are shown. As expected, when increasing the bandwidth the

Fig. 9 :

 9 Fig. 9: 2D mapping using 10 × 10 antenna arrays with W = 500 MHz (top-left), W = 2 GHz (top-right) and W = 3 GHz (bottom) with no quantization errors and a perfect knowledge of mobile position.

Fig. 10 :

 10 Fig. 10: 2D mapping using 4 × 4 (top-left), 10 × 10 (top-right) and 15 × 15 (bottom) antenna arrays with a perfect knowledge of mobile position and with array quantization errors.

Fig. 11 :

 11 Fig. 11: 2D mapping using 4 × 4 (top-left), 10 × 10 (top-right) and 15×15 (bottom) antenna arrays and with quantization and mobile position (10 cm) and orientation errors (5 • ).

Fig. 12 :

 12 Fig. 12: 2D mapping using 10 × 10 antenna arrays with mobile position errors of 0.5 m (left) and 1 m (right) and W = 1GHz.

Fig. 13 :

 13 Fig. 13: 2D mapping using 10 × 10 antenna arrays with mobile orientation errors of 10 • (left) and 20 • (right) and W = 1GHz.

Fig. 14 :

 14 Fig. 14: 2D mapping using 10×10 antenna arrays with position and mobile orientation errors 0.5 m+20 • (left) and 1 m+20 • and W = 1GHz.

  allow us to write the mean components of δΘ 0 asE δθ 0 ( Θ 0 , f ) = -E [S(Θ, f )] E ∂S(Θ,f ) ∂θ Θ= Θ0 E δφ 0 ( Θ 0 , f ) = -E [F (Θ, f )] E ∂F (Θ,f ) 46), let express S(Θ, f ) as S(Θ, f ) = ∂| AF(Θ, f )| 2 ∂θ = ∂|AF(Θ, f )| 2 ∂θ + ∂ζ err (Θ, f ) ∂θ . =k cos(θ) [(m -p)d x cos(φ)+(n -q)d y sin(φ)] = ψ mnpq (Θ, f )(48)it isE ∂ζ err (Θ, f ) ∂θ = j χ 2 -1 mnpq ψ mnpq (Θ, f ) e jαmnpq (49)where ψ mnpq = 0 when m = p, n = q. Consequently we haveE [S(Θ, f )] = jχ 2 mnpq ψ mnpq (Θ, f )e jαmnpq(50)

∂ 2 α mnpq ∂θ 2 =

 2 -k sin(θ) [(m -p)d x cos(φ)+(n -q)d y sin(φ)] = Λ mnpq (Θ, f ) .

E S 2 (

 2 Θ 0 , f ) = χ 4 mnpq rstu mn =pq =rs =tu ψ mnpq ( Θ 0 , f )ψ rstu ( Θ 0 , f ) + υ 2 -1 mnpq ψ 2 mnpq ( Θ 0 , f ) + 2 υ χ 2 mnpq rs mn =rs ψ mnpq ( Θ 0 , f )ψ rspq ( Θ 0 , f ) + 2χ 2 mnpq rs pq =rs ψ mnpq ( Θ 0 , f )ψ rsmn ( Θ 0 , f )(57)where υ = E e -j2δmn = E e j2δmn 1 ∆ sin(∆). We finally obtainE δθ 2 0 ( Θ 0 , f ) = E S 2 ( Θ 0 , f )

mnpq ŵmnpq e jδmnpq ∆ 2

 2

  mnpq ŵmnpq e jδmnpq (m-p) 2 k 2 d 2 x δ θ2 x0 cos 2 ( θ0 ) = |AF( Θ 0 , f )| 2 . (61)It is possible to writeE δ θx0 ( Θ 0 , f ) ≈ E   ± |AF( Θ 0 , f )| 2mnpq ŵmnpq e jδmnpq (m -p) 2 k 2 d 2 x cos 2 ( θ0 )

and in part by the H2020 IF-EF Marie-Curie project MAPS (Grant 659067). This paper

Davide Dardari (M'95-SM'07) received the Laurea degree in electronic engineering (summa cum laude) and the PhD degree in electronic engineering and computer science from the University of Bologna, Italy, in 1993 and 1998, respectively.

He is an Associate Professor at the University of Bologna at Cesena, Italy, where he participates with WiLAB (Wireless Communications Laboratory).


Yang Liu
email: liuyang@bucea.edu.cn

Mingyi Du

Wenquan Zhu

Examination Method and Implementation for Field Survey Data of Crop Types Based on Multi-resolution Satellite Images

Keywords: Examination method and implementation, Stratified random sampling, Automatic processing, Crop types, Multi-resolution satellite images

In order to examine the accuracy of large amount of the field survey data with less accurate, an examination method based on multi-resolution satellite images was proposed in this paper. As there were so large amount of data, stratified random sampling was used to obtain effective samples. Firstly, vegetation index derived from low-resolution satellite images at different times has been adopted as analysis factor. And wave curve charts were drawn with the vegetation index. From those charts, the statistics law of wave curves for different crop types was recognized using for crop types' classification. Secondly, high-resolution satellite images were used to correct the area of crop types to get the final classification results. Finally, the accuracy of the field survey data can be calculated by comparing the original survey data with the final classification results. Moreover, for convenience using, a software has been developed according to the above examination method.

Introduction

Planting area and yields of crops are the important basis for government's economic policy making. For a long time, the two kinds of traditional methods were used for the statistics of panting area of crops [START_REF] Wu | Crop Acreage Estimation Using Two Individual Sampling Frameworks with Stratification[END_REF]. One way is the comprehensive statistical report coming from statistical and administrative units at various levels step by step [START_REF] Wang | Spatial Sampling Design for Monitoring the Area of Cultivated Land[END_REF].

Another way is to sample field survey of reported data. This method is suitable for a large degree on the discrete variables; frequency distribution was highly skewed socio-economic phenomena for investigation [START_REF] Li | Spatial Sampling Optimized Decision-making of National Land Remote Sensing Investigation[END_REF]. But, no matter whichever method was used, one step can not be omitted, that is examining the accuracy of large amount data of the field survey. Only in this way, the data can become the reliable argument for the decision analysis.

In recent years, with the extensive application of remote sensing technology, the study of estimation and examination of crop area using remote sensing techniques has made significant progress, and is steadily moving towards the direction of business [START_REF] Narayanan | Effect of Spatial Resolution on Information Content Characterization in Remote Sensing Imagery based on Classification Accuracy[END_REF] [START_REF] Chen | Bulletin of Wheat, Corn and Rice Yields Estimate Technique with Remote Sensing[END_REF]. Compared with traditional ground survey and verification, the demand for manpower and costs are significantly reduced, using remote sensing image data [START_REF] Woodcock | The Factor of Scale in Remote Sensing[END_REF]. However, from the level of remote sensing technology for carrying out, large-scale remote sensing examination of crop area is facing three problems, namely, precision, efficiency and cost problems [START_REF] Liu | Sampling Method With Remote Sensing for Monitoring of Cultivated Land Changes on Large Scale[END_REF] [START_REF] He | Effects of Spatial Aggregation Approaches on Classified Satellite Imagery[END_REF].

Research Process

This article aims to use multi-resolution remote sensing imaging technology(two kinds of resolution images: the high-resolution images, such as ASTER images, and the low resolution images , such as MODIS images, Utilizing their own advantage, combined with the traditional random, stratified sampling method [START_REF] Foody | Status of Land Cover Classification Accuracy Assessment[END_REF], the examination scheme for field survey data of crop types was designed, In the scheme, stratified random sampling was used to ensure minimum sample size and stability, and the difference spectrum of low resolution remote sensing images was used to match the types of corps, and high-resolution remote sensing images was used to the visual identification of the types and area of corps. We try to solve the confliction of precision and cost, and try to provide methodological guidance for a large-scale examination for field survey data of the types and area [START_REF] Jiao | Paddy Rice Area Estimation Using a Stratified Sampling Method with Remote Sensing in China[END_REF].

The flow chart of technology route in this work, see Fig. 1, is expressed as follows: In order to make the multi-resolution satellite images matched each other strictly and be in correspondence with the field survey data, the procedure of remote sensing data processing should be done. And stratified random sampling was used to obtain effective samples of field survey data of crop types. Vegetation index derived from low-resolution satellite images at different times has been adopted as analysis factor.

After all the data have processed, the identification and verification should be done through visual observation method. The statistics law of wave curves of vegetation index at different times was recognized using for crop types' classification of the samples. While, high-resolution satellite images were used to correct the classified crop area of samples to get the classification results. Finally, the accuracy of the survey data of crop types can be calculated by comparing the original survey data with the final classification results using statistical analysis methods.

In this paper, we focus on stratified random sampling in this work.

Stratified Random Sampling

The effectiveness of stratified random sampling can be evaluated along two dimensions. First, it might reduce the total quantity of measurements required. Second, it can do a reduction of measurements would imply cost savings [START_REF] Chen | The Theoretic Structure and Computerized Realization of Stratified Audit Sampling[END_REF].

Overview

In the sampling method, stratified sampling is one of the most effective methods.

Compared with simple random sampling, stratified sampling has its several great advantages, such as the fewer number of samples, higher sampling precision and the lower cost. It is the effective way of large-scale statistical sampling surveys [START_REF] Wang | Probability Theory and Mathematical Statistic[END_REF]. So stratified sampling was adopted in the paper. Stratified sampling, also known as type sampling, is one of the most commonly used sampling techniques in practical work [START_REF] Brus | Designing Spatial Coverage Samples Using the K-means Clustering Algorithm[END_REF] [START_REF] Gallego | Stratified Sampling of Satellite Images with a Systematic Grid of Points[END_REF]. Firstly, overall sample is divided into several strata (groups) in accordance with certain rules in stratified sampling. Secondly, sampling was done within each stratum independently. The resulting sample is called stratified sample. Accordingly, the sample of each stratum is independent also [START_REF] Minasny | The Variance Quadtree Algorithm: Use for Spatial Sampling Design[END_REF]. Furthermore, if the sampling method of each stratum is simple random sample, the sample is called stratified random sampling. Thus the resulting samples are called stratified random samples.

Determination of Stratified Number

In practice, the number of strata can not exceed half of the sample size, because of ensuring that each stratum has at least two samples and having the need to calculate the standard deviation of each stratum in the method [START_REF] Kadilar | A New Ratio Estimator in Stratified Random Sampling[END_REF].In this paper, we stratify the population into six strata.

Determination of Stratified Boundaries

It is usually to determinate the boundaries of strata in accordance with characteristics or features of population in the stratification [START_REF] Tetsuji Ohyam | Estimating population Characteristics by Incorporating Prior Values in Stratified Random Sampling/Ranked Set Sampling[END_REF]. However, when the obvious characteristics or features is not easily recognized, the general way to determinate the boundaries of strata is studying general distribution of some properties in population or studying a relative relationship of the variable X and the variable Y which are subfeatures of population by means of some mathematical methods.

In this paper, the Optimal Stratification Method(also known as Accumulative Square Root Method) is used to determinate the boundaries of strata, which was proposed by Dalenius and Hodges [START_REF] Kadilar | Ratio Estimators for the Population Variance in Simple and Stratified Random Sampling[END_REF], Firstly, a characteristic property of population should be determined, according to which population is stratified. Secondly, population is sorted by the values of the property of population from small to large. Then, the sum of the values of the property of population is calculated, and the square root of the sum should be obtained [START_REF] Wu | Empirical Likelihood Inference Under Stratified Random Sampling in the Presence of Measurement Error[END_REF]. Finally, population was divided into several strata by equal division method of the square root.

Calculation of Sample Size

There are two main steps in the calculation of sample size of each stratum as show in the below.

Computation of the total sample size. The total sample size calculations generally use the following:

∑ ∑ = = + = L h h h L h h h h S W N t Y r w S W n 1 2 2 1 2 2 1 ) ( (1)
In Equation [START_REF] Wu | Crop Acreage Estimation Using Two Individual Sampling Frameworks with Stratification[END_REF], n is the total sample size, and h is the strata-specific variable which ranges from 1 to L. L is number of strata. W h is the Weight of stratum h, which can be calculate by means of N h (the population for stratum h) divide by the total stratified sample size N. S h 2 is the population deviation of stratum h, which is calculated using Equation [START_REF] Wang | Spatial Sampling Design for Monitoring the Area of Cultivated Land[END_REF]. And w h is the sampling radio of stratum h, which is calculated using Equation (3). And r is the relative error limit (also as known the confidence interval) for (general set 95%). While, Y is the population mean, and t is the percentile of the standard normal distribution (z=0.05 for 95% confidence).

S h 2 is determined by:

∑ = - = h N i h hi h Y y N S 1 2 2 ) ( 1 (2)
Note that the value of sample i of stratum h is marked with a hat as y hi , h Y is the population mean for stratum h. Other symbols are the same as Equa tion [START_REF] Wu | Crop Acreage Estimation Using Two Individual Sampling Frameworks with Stratification[END_REF]. And w h is determined by:

∑ = = L h h h h h h S W S W w 1 (3)
Where, S h is the population standard deviation of stratum h, namely the square root of S h 2 . Also, other symbols are the same as Equation (1) or Equation [START_REF] Wang | Spatial Sampling Design for Monitoring the Area of Cultivated Land[END_REF].

Sample allocation in each stratum. For stratified sampling, the sample size of each stratum still need to determine after the total sample size is fixed. When doing the population estimation, the population variance which can be estimated is related not only to the variance of each stratum, but also to the sample size of each stratum.

There are a lot of sample allocation methods in practice [START_REF] Wright | Equal-Precision Allocations and Other Constraints in Stratified Random Sampling[END_REF]. You can allocate the samples in accordance with the radio of between the sample size of each stratum and the sample size of population, or according to desired overall confidence interval with minimum total sample size.

The sample allocation of each stratum was done using the Optimal Sample Allocation, which was proposed by Neyman [START_REF] Bosch | Optimum Allocation of Stratified Random Samples Designed for Multiple Mean Estimates and Multiple Observed Variables[END_REF]. It is expressed mathematically as follows:

∑ = × = L h h h h h h S N S N n n 1 (n h shall be as an integer) (4)
where, N h is the sample size of stratum h. Also, other symbols are the same as Equation (1) or Equation (3).

However, the two problems should be noted in the Optimal Sample Allocation. There are:

One is in the case of the sampling ratio f = n / N is very large. This situation leads to S h (the standard deviation of each stratum) is relatively large, and n h (the sample size of each stratum) may be larger than N h (the population for stratum h). At this situation, it required for 100% sampling to the population for stratum h.

The other case is that n h is less than 2 after the calculation using the Optimal Sample Allocation. The value of n h needs to be set to 2 to reduce the impact of random error on the results and the stability when the sample size is 1.

Examination

The process of examination for field survey data of crop types can be mainly divided into the following four steps. Note that, the experiment of examination itself only use a small part of the test data.

4.1

Fig. 2. This shows the experimental data in this work. From left to right is successively vegetation index at different times, high-resolution satellite images, field survey data and the overlay of the above three data which share the same region in the rectangular box.

We have collected the latest report information of field survey data of crop types of county-level, the multi-resolution satellite images which can cover the survey region in this step and prepare for the next step of Sampling.

The field survey data of crop types are Shape-file data, which should include the type and area property fields. High-resolution satellite data now is ASTER images. But if required, it can be higher resolution images. Vegetation index at different times were derived from MODIS time series data, which were formed as a multi-band image (namely one time corresponds to one band). And all test data was matched strictly

Sampling

Fig. 3 shows sampling result comparing with raw data, and the raw data was covered by the sampling data. So the visible parts of the raw data are not sampled. We also do the statistics to the sampling result which is convenient for customers to find out sampling information to decide whether to re-sampling, shown as the Table 1 below.

Crop Types' Identification and Areas' Correction

In this step, multi-resolution satellite images were used to determine crop types and area of the samples by visual observation method. So there are two factors as following:

Crop types' identification. Wave curve charts were drawn with the index of bands as X-axis and the value of vegetation index as Y axis, as shown in Fig. 4. There are 17 bands in this work. Fig. 5 show wave curve charts of four different positions. From those charts, the statistics law of wave curves for different crop types was identified using for crop types' classification. Crop areas' verification. It provides a tool to draw polygons and calculate area of those. By intercomparison between high-resolution satellite images and samples of the same region, you can get the area of samples by measuring. The Area of all samples should be measured, which has been regarded as hypothetical real value.

Statistical Analysis

It needs to get error situation of the field survey data by doing statistics for result of type identification and area verification of samples.

As shown in Table 2, the number of qualified area, the number of qualified type and the number of qualified both of each stratum was obtained by statistics. Also, qualified rate of area, qualified rate of type, and overall qualified rate (area and type are both qualified) of samples were calculated by anti-deduction. In statistics, we can find out the abnormal samples which are the types is right, but the rate of area change exceeding a certain threshold, as shown in Figure 6.In this works, the threshold was set to 10%. These abnormal samples should be report to the related departments of doing field survey to redo site verification.

System Implementation

An Examination System for Field Survey Data (ESFSD) has been developed based on the above examination method. The system (ESFSD) provides integrated and processoriented functions, such as data processing, stratified random sampling, sample identification and verification, statistical analysis of samples, and so on. Using the system, you can easily complete the entire process of examination for field survey data of crop types.

The system (ESFSD) adopts C/S (client/server) software system structure, C# as the programming language and SQL Server 2000 as back database server. Using ArcEngine as development tools, the application program developed has been closely integrated with system by UG/Open MenuScript. Figure 7 presents the interface of ESFSD.

ESFSD also provides data management statistics information of examination result, and supply better decision supports for re-sampling of samples. Fig. 7. This shows the system interface.

Conclusion

In order to examine the accuracy of large amount of the field survey data, an examination method based on multi-resolution satellite images was proposed in this paper. For reducing the number of data and cost in the examination, stratified random sampling was used to obtain effective samples. For determining crop types and area of the samples correctly and quickly, multi-resolution satellite images were used utilizing their own advantage. Low-resolution satellite images are used to identify the corps type of samples, and high-resolution satellite images were used to verify the area of samples. Finally, the accuracy of the field survey data can be calculated by comparing the original survey data with the final examination results by statistics. By the testing of actual data, the examination method is very effective to solve the confliction of precision and cost. And it also provides methodological guidance for a large-scale examination for field survey data of the types and area. Moreover, for convenience using, the system (ESFSD) has been developed according to the above examination method. The results of system testing show ESFSD is comprehensive processing platform which has realized integrated and process-oriented functions for the examination. And it has strong robustness and stability.

Fig. 1 .

 1 Fig. 1. This shows the flow chart of technology route.

Fig. 3 .

 3 Fig. 3. This shows the sampling result comparing with raw data.

Fig. 4 .

 4 Fig. 4. This shows the wave curve charts of vegetation index.

Fig. 5 .

 5 Fig. 5. This shows the schematic diagram of comparison with different wave curves.

Fig. 6 .

 6 Fig. 6. This shows the abnormal samples.

Table 1 .

 1 Statistics of sampling data.

	Index of Strata	Boundaries	Num of Population	Num of Samples	Standard Deviations	Means
	1	100.00 -890.13	109	2	119.42	388.07
	2	890.13 -1780.25	41	2	81.26	1212.19
	3	1780.25 -2670.38	22	2	61.56	2136.36
	4	2670.38 -3560.50	12	2	33.11	2974.99
	5	3560.50 -4450.63	2	2	14.76	4000.00
	6	4450.63 -3217100.00	181	166	273122.73 108427.62
	Total is	100.00 -3217100.00	367

Table 2 .

 2 Examination result.

	Index Of	Boundaries	Num of	Num of	Num of	Num of	Num of
	Strata		Population	Samples	qualified	qualified	qualified
					area	type	Both
	1	100.00 -890.13	109	2	2	1	1
	2	890.13 -1780.25	41	2	2	2	2
	3	1780.25 -2670.38	22	2	1	1	0
	4	2670.38 -3560.50	12	2	2	2	2
	5	3560.50 -4450.63	2	2	1	1	0
	6	4450.63 -3217100.00	181	166	83	81	23
	Total is	100.00 -3217100.00	367	176	91	88	28
	Conclusion is Qualified rate of area is 72.07% Qualified rate of	Overall qualified rate is
			type is 56.63%	36.12%

Acknowledgments. The research performed is supported by Beijing Natural Science Foundation (No: 4102019), and is partially funded by Jurisdiction of Beijing Municipality under Grant No. PHR200906138, PHR200907127 and PHR 20070101. The work was performed under Beijing University of Civil Engineering and Architecture (No:100901405). We wish to thank Guoyin Cai for his careful proofreading of the manuscript.

Data

The entire test data used in this paper was shown in Fig. 2.