Yifan Du
email: yifan.du@inria.fr

Françoise Sailhan
email: francoise.sailhan@cnam.fr

Valérie

IAM -Interpolation and Aggregation on the Move: Collaborative Crowdsensing for Spatio-temporal Phenomena

Keywords: Crowdsensing, Ubiquitous Sensing, Pervasive Computing, Aggregation, Interpolation, Opportunistic Relay

Crowdsensing allows citizens to contribute to the monitoring of their living environment using the sensors embedded in their mobile devices, e.g., smartphones. However, crowdsensing at scale involves significant communication, computation, and financial costs due to the dependence on cloud infrastructures for the analysis (e.g., interpolation and aggregation) of spatio-temporal data. This limits the adoption of crowdsensing by activists although sorely needed to inform our knowledge of the environment. As an alternative to the centralized analysis of crowdsensed observations, this paper introduces a fully distributed interpolation-mediated aggregation approach running on smartphones. To achieve so efficiently, we model the interpolation as a distributed tensor completion problem, and we introduce a lightweight aggregation strategy that anticipates the likelihood of future encounters according to the quality of the interpolation. Our approach thus shifts the centralized postprocessing of crowdsensed data to distributed pre-processing on the move, based on opportunistic encounters of crowdsensors through state-of-the-art D2D networking. The evaluation using a dataset of quantitative environmental measurements collected from 550 crowdsensors over 1 year shows that our solution significantly reduces -and may even eliminate-the dependence on the cloud infrastructure, while it incurs a limited resource cost on end devices. Meanwhile, the overall data accuracy remains comparable to that of the centralized approach.

INTRODUCTION

Mobile crowdsensing is an attractive sensing paradigm to monitor the urban environment such as noise level, ambient temperature, luminance, air pressure, relative humidity, etc. [START_REF] Issarny | Dos and don'ts in mobile phone sensing middleware: Learning from a large-scale experiment[END_REF][START_REF] Kang | Enhance the quality of crowdsensing for fine-grained urban environment monitoring via data correlation[END_REF][START_REF] Qin | NoiseSense: A crowd sensing system for urban noise mapping service[END_REF][START_REF] Xu | Urban noise mapping with a crowd sensing system[END_REF]. Citizens may contribute valuable spatio-temporal observations using the low-cost, yet powerful, sensors embedded in their smartphones/tablets, with the GPS for positioning and the Internet access for uploading. Furthermore, pollution is an increasing societal concern, and through opportunistic crowdsensing, people are able to be aware of it in a mutually beneficial way. However, the collected observations often cover the urban space and time unevenly and sparsely, due to the dependence on the mobility of citizens throughout the city [START_REF] Wang | Sparse mobile crowdsensing: challenges and opportunities[END_REF]. State-of-the-art crowdsensing systems address this shortcoming through the centralized analysis -incl. aggregation and interpolation-of observations provided to cloud infrastructure servers [START_REF] Girolami | Sensing interpolation strategies for a mobile crowdsensing platform[END_REF][START_REF] Kong | Data loss and reconstruction in sensor networks[END_REF][START_REF] Mendez | Data interpolation for participatory sensing systems[END_REF][START_REF] Ventura | Assimilation of mobile phone measurements for noise mapping of a neighborhood[END_REF]. The implemented centralization then severely limits the adoption of crowdsensing for environmental monitoring due to the resulting resource and financial costs, and also introduces user privacy leak (e.g., mobility inference) [START_REF] Shi | Prisense: privacy-preserving data aggregation in people-centric urban sensing systems[END_REF].

We argue that enabling fully decentralized crowdsensing systems, including the underlying large-scale data analysis, is key to the democratization of environmental monitoring using crowdsensors. This is the focus of our paper that introduces the IAM (Interpolation and Aggregation on the Move) distributed solution, which builds upon the two following trends to support crowdsensingbased environmental monitoring at scale: (i) IAM views crowdsensors as "social sensors" that outfit human, as opposed to mere physical equipment that senses the environment and transfers data. In particular, mobile crowdsensors often encounter each other during sensing campaigns or as part of the daily routine of their owners [START_REF] Lee | Comon: Cooperative ambience monitoring platform with continuity and benefit awareness[END_REF]. We focus on opportunistic crowdsensing, where people follow their daily routine without being directly involved in the sensing task [START_REF] Du | Let opportunistic crowdsensors work together for resource-efficient, quality-aware observations[END_REF][START_REF] Issarny | Dos and don'ts in mobile phone sensing middleware: Learning from a large-scale experiment[END_REF]. At the micro level, behavioral signatures (i.e., routines) as well as recurrent encounter patterns reflect the underlying relational dynamics of organizations or communities to which the user is affiliated. Furthermore, end-users are more likely to share their data, especially in short encounters, as the risk of losing their anonymity is lower [START_REF] Juhlin | Time to meet face-to-face and deviceto-device[END_REF]. (ii) The collaborative and ubiquitous processing of crowdsensed observations improves the overall efficiency of the system in terms of resource consumption [START_REF] Du | Let opportunistic crowdsensors work together for resource-efficient, quality-aware observations[END_REF][START_REF] Sailhan | Opportunistic multiparty calibration for robust participatory sensing[END_REF][START_REF] Xiao | Lowering the barriers to large-scale mobile crowdsensing[END_REF], regarding both the infrastructure and contributing devices. It also allows end users to collaborate and share knowledge effortlessly and with little or no cost. In particular, D2D collaboration enhances mobile edge computing by enabling the sharing of heterogeneous computing and communication resources between powerful mobile end devices [START_REF] Xu Chen | Exploiting massive D2D collaboration for energy-efficient mobile edge computing[END_REF][START_REF] Lee | Comon: Cooperative ambience monitoring platform with continuity and benefit awareness[END_REF].

Concretely, IAM allows a fully distributed, collaborative approach to crowdsensing: crowdsensors interpolate their observation of the phenomenon, and they aggregate the respective data in an opportunistic way. The intent is to overcome the spatio-temporal sparsity and to limit -or even avoid-the use of a centralized infrastructure server. There are many interpolation methods for inferring spatio-temporal phenomena, and the smartphone is becoming increasingly powerful to perform such advanced tasks. At the same time, P2P wireless ad hoc network technology (e.g., Wi-Fi Direct) enables the discovery of nearby mobile devices and the exchange of data between peers, making the collaboration possible [START_REF] Du | Let opportunistic crowdsensors work together for resource-efficient, quality-aware observations[END_REF][START_REF] Lee | Comon: Cooperative ambience monitoring platform with continuity and benefit awareness[END_REF]. Some crowdsensing systems already exploit the P2P collaboration of crowdsensors as they meet [START_REF] Do | Crowdsourced mobile data transfer with delay bound[END_REF][START_REF] Wang | ecoSense: Minimize participants' total 3G data cost in mobile crowdsensing using opportunistic relays[END_REF][START_REF] Yang | EdgeSense: Edge-mediated spatial-temporal crowdsensing[END_REF]. However, the collaboration primarily deals with handling the relay of data, while deployed static edge servers are in charge of the distributed data aggregation.

Our approach leverages the advantage of the former and overcomes the disadvantage of the latter: it implements an opportunistic data relay and analysis on the move, across the crowdsensors.

In summary, the paper makes the following contributions to enable ubiquitous and collaborative crowdsensing:

(1) A fully distributed approach to the aggregation and interpolation of crowdsensing data on the move ( § 3), which exploits the smartphone's capability to perform 3D tensor completion. In particular, we thoroughly analyze and compare eligible stateof-the-art interpolation methods, which leads us to leverage the Gaussian Process Regression ( § 3.2) that produces the most cost-effective inference for spatio-temporal phenomena along with an estimation of the inference uncertainty.

Another advantage compared to alternative interpolation approaches (i.e., ordinary kriging & tensor decomposition) is that all the specified parameters are learned and the most relevant kernel was selected following related assessment. Finally, only a small portion of the interpolation is running on each crowdsensor. (2) A distributed lightweight and quality-aware aggregation strategy based only on linear operations ( § 3.3) that combines the tensors that each crowdsensor establishes autonomously. Such a linear aggregation takes place opportunistically across the end devices, and possibly at the server if the D2D communication does not allow covering all the devices over the given time window. That is, when crowdsensor peers get in a shared D2D communication range, our algorithm selects one of them to aggregate their interpolated tensors and further relay the new tensor. At the end of a predefined time window, the crowdsensor uploads its tensor to the server unless it has previously met a peer that takes care of the relay. A key aspect of the proposed approach is that the aggregation is much less resource consuming (in time and space) than the interpolation, and the server performs only the aggregation needed to fill the gap between the network islets that the contributing crowd forms through D2D communication.

(3) A prototype implementation ( § 4) of the IAM solution and its performance evaluation using a 1-year crowdsensing dataset ( § 5). The evaluation shows that IAM aggregates a global knowledge that is both robust and accurate compared to the centralized approach. Based on our empirical evaluation, we select the best kernel, which is characterized by a fairly good balance between accuracy and resource consumption. Most importantly, compared to the centralized approach and baseline relay-based aggregation mechanisms, IAM significantly lowers the communication, computation and financial costs of crowdsensing-based environmental monitoring.

RELATED WORK

Sensing data fusion deals with the combination of observations so as to enhance the quality of the knowledge we gather from the sensors. In the specific case of fusing data from mobile crowdsensors and assuming that all the crowdsensors are trusted to provide equally accurate measurements, the supporting algorithms serve the two following functions: (i) aggregating together the measurements associated with related observations, and (ii) interpolating the missing measurements to overcome the sparse contribution coverage.

Data Aggregation

Related crowdsensing observations are usually aggregated using an average [START_REF] Koukoutsidis | Estimating spatial averages of environmental parameters based on mobile crowdsensing[END_REF][START_REF] Shi | Prisense: privacy-preserving data aggregation in people-centric urban sensing systems[END_REF] or a weighted average [START_REF] Jin | Incentive mechanism for privacyaware data aggregation in mobile crowd sensing systems[END_REF][START_REF] Sailhan | Opportunistic multiparty calibration for robust participatory sensing[END_REF] value, although more complex functions may be found in the literature depending on the observed phenomenon [START_REF] Jesus | A survey of distributed data aggregation algorithms[END_REF]. The aggregation may be executed either on the cloud, or in a distributed way -at least partiallyfor which most crowdsensors provide the necessary computing resources. We undertake the latter decentralized approach in our work so as to benefit from ubiquitous computing and in particular limit the dependence on -and related resource and financial costs due to the usage of-a cloud/edge infrastructure. The challenge is therefore to perform a distributed aggregation in a way that both (i) delivers an overall accurate knowledge, and (ii) incurs bearable resource consumption for the end devices. To achieve so, we use the Average aggregate function that is duplicate-sensitive and decomposable; in particular, a batch aggregation is equivalent to several pairwise aggregations. The distributed aggregation protocol may rely on either a structured (e.g., hierarchical) or an unstructured (e.g., flooding) routing; it may even implement a combination of the two. Our solution is based on the unstructured routing that matches the mobility behavior of opportunistic crowdsensors. That is, a crowdsensor exchanges its observations with another crowdsensor as they get in the D2D communication range of each other, so that one of them aggregates the two sets of observations, and in turn repeats the process as it meets a new crowdsensor. The protocol stops when reaching a predefined criterion, which is a given time-window in our case (see § 4). Still, the distributed aggregation protocol must be complemented with the interpolation of missing values [START_REF] Wang | Learning-assisted optimization in mobile crowd sensing: A survey[END_REF].

Data Interpolation

There exist various eligible approaches to the interpolation of crowdsensed observations, wherein the challenge is to overcome both the related data sparsity and computing complexity. In statistical geography, multivariate interpolation and spatial interpolation play an important role as they enable modeling a large-scale phenomenon (and producing a digital elevation model) provided a set of observations/points. Many interpolation techniques may be applied, depending on the characteristics of the observed data points.

Ordinary kriging [START_REF] Girolami | Sensing interpolation strategies for a mobile crowdsensing platform[END_REF][START_REF] Mendez | Data interpolation for participatory sensing systems[END_REF] is one of the methods that is widely applied as it supports a fine-grained interpolation at each location over a 2D space. It estimates the value at an unobserved location based on values at observed locations and it performs well as long as observations are uniformly distributed. In particular, performance greatly degrades for large amounts of missing data and the computational cost gets prohibitive [START_REF] Eldrandaly | Spatio-temporal Interpolation: Current Practices and Future Prospects[END_REF].

Compressive sensing is a recent alternative approach to the interpolation of data for the production of observation maps, while dealing with sparse observations [START_REF] Eldrandaly | Spatio-temporal Interpolation: Current Practices and Future Prospects[END_REF]. It has been used to infer/interpolate urban-scale physical phenomena from crowdsensed observations stored in a 2D matrix [START_REF] Girolami | Sensing interpolation strategies for a mobile crowdsensing platform[END_REF][START_REF] Wang | Sparse mobile crowdsensing: challenges and opportunities[END_REF][START_REF] Xu | More with less: Lowering user burden in mobile crowdsourcing through compressive sensing[END_REF]. However, physical phenomena often have more than two modes of variation and are therefore best represented as multi-dimensions observations, which leads to using 3D tensors [START_REF] Kang | Data correlation based crowdsensing enhancement for environment monitoring[END_REF][START_REF] Qin | NoiseSense: A crowd sensing system for urban noise mapping service[END_REF][START_REF] Zheng | Diagnosing New York city's noises with ubiquitous data[END_REF]. The matrix/tensor completion problem is usually solved by matrix/tensor decomposition and one major drawback is that it trades uncertainty for efficiency, as there is no confidence interval associated with the inference.

Multiple regression is another approach to infer missing values [START_REF] Harris | The use of geographically weighted regression for spatial prediction: an evaluation of models using simulated data sets[END_REF][START_REF] Ohashi | Spatial interpolation using multiple regression[END_REF], which transforms the interpolation problem into regressions. This inspired us to use the Gaussian Process Regression (GPR) to solve the tensor completion problem. GPR is a well-known and general approach applied in machine learning [START_REF] Liu | When Gaussian process meets big data: A review of scalable GPs[END_REF], which supports a non-parametric and interpretable Bayesian model. GPR naturally supports multi-dimension and it does not require the data to follow Gaussian distribution. GPR enables estimating the level of uncertainty associated with the produced model. We compare the three interpolation approaches in the evaluation section ( § 5) where we show that GPR is the best suited to support distributed interpolation across mobile crowdsensors.

Centralized vs Distributed Data Fusion in Crowdsensing Platforms

The great majority of crowdsensing platforms applies centralized analysis to the sensing data: the cloud/edge server first aggregates the raw crowdsensor measurements and then interpolates the missing observations. Existing platforms adopt various interpolation and aggregation methods, while considering different types of sensors (static vs mobile). For instance, one of the solutions leverages ordinary kriging on the cloud to generate a map from the measurements contributed by a combination of static sensors and mobile crowdsensors [START_REF] Girolami | Sensing interpolation strategies for a mobile crowdsensing platform[END_REF]. Some platforms [START_REF] Kang | Data correlation based crowdsensing enhancement for environment monitoring[END_REF][START_REF] Kang | Enhance the quality of crowdsensing for fine-grained urban environment monitoring via data correlation[END_REF][START_REF] Qin | NoiseSense: A crowd sensing system for urban noise mapping service[END_REF][START_REF] Xu | Urban noise mapping with a crowd sensing system[END_REF] exploit additional urban data sources (e.g., road networks, check-in data) to infer the phenomenon, which is represented as a 3D tensor completion using centralized tensor decomposition. Rather than integrating more data sources, the analysis may leverage historical data to improve the accuracy of the generated map. A multi-Output Gaussian Process serves as a unified map generation model, which takes multiple instances (a current sparse instance and an appropriate historical dense instance) to generate an improved air quality map in [START_REF] Cheng | MapTransfer: Urban Air Quality Map Generation for Downscaled Sensor Deployments[END_REF].

The distributed interpolation and aggregation of sensing data have deserved less attention than the centralized counterparts, with most solutions targeting Wireless Sensor Networks (WSN). A localized, distributed interpolation & aggregation scheme based on kriging is introduced in [START_REF] Umer | Spatial interpolation in wireless sensor networks: localized algorithms for variogram modeling and Kriging[END_REF] for a tree-structured WSN; it allows inferring a phenomenon over the holes that the static WSN does not cover. In [START_REF] Xu | Hierarchical data aggregation using compressive sensing (HDACS) in WSNs[END_REF], a hierarchical WSN is organized such that sensors send the measurements to their respective cluster heads, which in turn encode the sparse measurements and use compressive sensing (matrix decomposition) to interpolate the overall phenomenon.

The work that is the most related to ours is the edge-mediated spatial-temporal crowdsensing proposed in [START_REF] Yang | EdgeSense: Edge-mediated spatial-temporal crowdsensing[END_REF]. The solution relies on a trusted edge server (e.g., a deployed cloudlet) that coordinates a few crowdsensors. In practice, crowdsensors independently perform a part of matrix decomposition using a stochastic gradient descent and they exchange factor vectors with the crowdsensors that are in the same WLAN. Each edge server is responsible for a batch of crowdsensors that are fully connected to perform iterative optimization, and ultimately recovers an interpolated map for the sub-area. The set of crowdsensors and the edge server must remain connected and communicate over multiple rounds to establish a single interpolation. Furthermore, as stated previously, the tensor decomposition that is applied is less accurate than GPR, which we leverage within IAM.

Different to previous work, the IAM solution: (1) leverages 3D tensors that embed more information than the 2D matrix data model, and efficient Gaussian Process Regression for interpolation, (2) supports both interpolation and aggregation at the end device in a resource-efficient way, so as to limit the dependence on the infrastructure (edge/cloud server) at a bearable additional resource (including energy) cost for the users' crowdsensing devices, and (3) implements opportunistic P2P aggregation, which benefits from the encounter with other crowdsensors and is hence not constrained by any hierarchical/tree network structure.

LIGHTWEIGHT DECENTRALIZED CROWDSENSING DATA FUSION

The IAM solution to decentralized crowdsensing takes benefits of today's smartphones capacity, while limiting the resulting additional resource consumption on devices. That is, IAM implements lightweight collaborative sensor data fusion across the participating crowdsensing devices as they are in D2D communication range.

Problem Statement

Let IAM be deployed over 𝑚 mobile crowdsensors (e.g., smartphones) embedding (built-in or connected) sensors providing measurements related to the physical phenomenon P. We further assume that the 𝑚 crowdsensors are all trustworthy and provide equally accurate measurements (i.e., they underwent appropriate calibration [START_REF] Sailhan | Opportunistic multiparty calibration for robust participatory sensing[END_REF] prior to contribute measurements). IAM supports the periodic monitoring of P over (a possibly large) area A and a given time period D using the contributions of the 𝑚 crowdsensors.

We represent the data that each crowdsensor collects as a concise 3D tensor where the first two dimensions refer to the spatial space (i.e., latitude and longitude) and the third one to the temporal domain (time). In particular, we discretize the target region A into 𝐼 × 𝐽 areas, which are cells of equal spatial size. We also discretize D into 𝐾 time slots of equal durations. We denote Y 𝑠 ∈ R 𝐼 ×𝐽 ×𝐾 the tensor that crowdsensor 𝑠 (1 ≤ 𝑠 ≤ 𝑚) maintains. The entry 𝑦 𝑠 (𝑥) ∈ R at position 𝑥 := (𝑖, 𝑗, 𝑘) ∈ R 3 is the average of the measurements collected by 𝑠 over the area cell indexed by (𝑖, 𝑗) during the time interval indexed by 𝑘. The value 𝑦 𝑠 (𝑥) is null if 𝑠 does not sense at position 𝑥. In other words, any crowdsensor 𝑠 contributes a tensor Y 𝑠 that provides a sparse/incomplete observation of the physical phenomenon according to its behavior.

IAM allows the opportunistic combination of the various tensors Y 𝑠 from the contributing crowdsensors 𝑠, using interpolation and aggregation, so as to compute an overall Y for a spatio-temporal characterization of phenomenon P over area A and duration D. As a base design choice, crowdsensors first apply interpolation over their local tensor Y 𝑠 prior to engage in the collaborative aggregation, which is key to the -both local and global-resource-efficiency of the IAM approach, while supporting the computation of a globally accurate knowledge.

Spatio-temporal Interpolation using Gaussian Process Regression

Given a sparse tensor Y 𝑠 resulting from the averaging of the local measurements collected at crowdsensor 𝑠 over area A and duration D, interpolation allows completing the tensor by estimating missing cells. The resulting (denser) tensor is denoted Ŷ𝑠 . The quality of the estimation can be established based on the approximation error (i.e., residual). The overall residual is then given by E = ∥Y -Ŷ ∥.

Let Ω be the set of observed cells by a given crowdsensor, that is, the cells to which the crowdsensor contributed observations. The Boolean mask tensor M ∈ B 𝐼 ×𝐽 ×𝐾 is defined such that 𝑚(𝑥) = 1 if there is a corresponding value at point 𝑥 ∈ Ω, and 𝑚(𝑥) = 0 otherwise. Thus, M * Y provides values resulting from actual observations (i.e., ground truth as sensed). When estimating Ŷ based on a sparse tensor Y with mask M, we seek to minimize the following loss function, which is associated with Ŷ:

𝐽 (Y, Ŷ) := 1 2 𝑥 ∈Ω 𝑒 (𝑥) 2 = 1 2 ∥M * ( Ŷ -Y)∥ 2
where: 𝑒 (𝑥) is the residual at point 𝑥, ∥ • ∥ denotes the Euclidean norm of a tensor, and * represents the element-wise multiplication.

Recall that the function 𝑦 : R 3 ↦ → R maps an arbitrary point 𝑥 := (𝑖, 𝑗, 𝑘) to its cell value 𝑦 (𝑥). Following the assessment of the eligible interpolation methods summarized in Section 2.2, we leverage the Gaussian Process Regression [START_REF] Liu | When Gaussian process meets big data: A review of scalable GPs[END_REF] to compute Ŷ out of Y. That is, we assume that 𝑦 follows a Gaussian Process (Gaussian distribution over functions), i.e.: 𝑦 (𝑥) ∼ GP (𝜇 (𝑥), 𝑘 (𝑥, 𝑥 ′ ))

where: 𝜇 (𝑥) = E[𝑦 (𝑥)] refers to the mean function, and 𝑘 (𝑥, 𝑥 ′ ) is the covariance matrix, i.e., the kernel of the GPR, which verifies

𝑘 (𝑥, 𝑥 ′ ) = E[(𝑦 (𝑥) -𝜇 (𝑥))(𝑦 (𝑥 ′ ) -𝜇 (𝑥 ′ ))].
The kernel is a crucial ingredient of GPR as it encodes the notion of similarity between two nearby data points 𝑥 and 𝑥 ′ on the basis that observations that are close to each other (Euclidean distance) are likely to have higher correlation. Thereby, the actual measurements that are close to an approximated observation are assumed to be highly informative for the inference at that point. Various families of kernels exist (see [START_REF] Edward | Gaussian processes in machine learning[END_REF] for an overview). In our case, the Matérn [START_REF] Minasny | The Matérn function as a general model for soil variograms[END_REF] induced the lowest error and execution time compared to alternative kernels (see §5).

In summary, once a GPR model is trained, the inferred mean value μ (𝑥) and its variance σ (𝑥) of any input point 𝑥 are generated (as regression process). The complete approximation tensor Ŷ is produced using μ (𝑥) as ŷ (𝑥) cell value, while a variance tensor Σ2 ∈ R 𝐼 ×𝐽 ×𝐾 is also maintained, in which each element σ2 (𝑥) refers to the variance of corresponding ŷ (𝑥).

GPR is computationally demanding as the training scales in 𝑂 (𝑛 3 ) with 𝑛 being the number of observed cells. Thus, applying such a regression over the overall dataset on the cloud incurs significant computation costs. As an alternative, IAM distributes the training and inference load over the crowdsensors, which further results in a limited computational cost on the device due to the relative low number of contributed cells (see § 5).

Opportunistic P2P Aggregation

Upon the P2P meeting (i.e., discovery through D2D communication) of two crowdsensors 𝑠 and 𝑠 ′ , IAM selects one of them to aggregate their respective tensors Ŷ𝑠 and Ŷ𝑠 ′ . Assuming the selected crowdsensor is 𝑠, then Ŷ𝑠 is updated as the aggregation result of Ŷ𝑠 and Ŷ𝑠 ′ , and Ŷ𝑠 ′ is set to null. Then, 𝑠 may relay Ŷ𝑠 when it meets another crowdsensor, or till the current time window 𝑇 +D expires.

The P2P meetings that IAM fosters correspond to a stochastic process since crowdsensors meet each other in an opportunistic way based on their own mobility. Upon such a meeting, a straightforward aggregation approach would consist in selecting randomly one of the two crowdsensors to take in charge the aggregation and relay the resulting tensor, where we assume that all the crowdsensors have equal resource budgets. Instead, we use the respective quality of inference of the crowdsensor tensors for the relay selection. The evaluation results confirm the relevance of the criterion, and suggest it is an indicator of future meeting occurrences (see § 5). Precisely, we leverage the inference quality as defined by the following asymmetric and positive loss function: Y 𝑠 , M 𝑠 , 𝑛 𝑠 , Σ 2 𝑠 ← 𝑛𝑢𝑙𝑙 5: end if Algorithm 1 introduces the aggregation procedure that crowdsensor 𝑠 runs upon meeting with crowdsensor 𝑠 ′ (with 𝑠 ′ running the same algorithm). Only one of the two should perform the actual aggregation (merge the tensors) and act as the relay node, which we call the mainstay. The crowdsensor with the lowest loss function selects itself (Lines 1-2) as the mainstay. The other crowdsensor no longer maintains its local data (Lines 3-4) so that there is no duplicated uploading. We highlight that the selection of the mainstay aims at optimizing the quality of the data delivered by the distributed aggregation process for which we consider the inference quality as criterion. It is area for future work to increase the overall robustness of the opportunistic aggregation protocol by taking into account additional criteria (e.g., mobility behavior, available resource, and fault tolerance), for which we may leverage state of the art algorithms [START_REF] Wang | ecoSense: Minimize participants' total 3G data cost in mobile crowdsensing using opportunistic relays[END_REF], while keeping the process energy-efficient.

𝐷 ( Ŷ𝑠 , Ŷ𝑠 ′ ) := ∥(M 𝑠 ′ * ¬M 𝑠 ) * ( Ŷ𝑠 ′ -Ŷ𝑠 )∥ 2 2∥M 𝑠 ′ *
Algorithm 2 details the data merge function that the mainstay 𝑠 runs, provided the tensor Y 𝑠 ′ from 𝑠 ′ , to compute the new tensor Y 𝑠𝑠 ′ . The algorithm distinguishes whether the values from 𝑦 𝑠 (𝑥)

Y 𝑠𝑠 ′ ← 0 𝐼 ×𝐽 ×𝐾 2: Y 𝑠𝑠 ′ + = 𝑛 𝑠 Y 𝑠 +𝑛 𝑠 ′ Y 𝑠 ′ 𝑛 𝑠 +𝑛 𝑠 ′ * (M 𝑠 * M 𝑠 ′ ) 3: Y 𝑠𝑠 ′ + = Y 𝑠 * (M 𝑠 * ¬M 𝑠 ′ ) 4: Y 𝑠𝑠 ′ + = Y 𝑠 ′ * (M 𝑠 ′ * ¬M 𝑠 ) 5: Y 𝑠𝑠 ′ + = (𝛽 𝑠 Y 𝑠 /Σ 2 𝑠 + 𝛽 𝑠 ′ Y 𝑠 ′ /Σ 2 𝑠 ′ )/(𝛽 𝑠 Σ -2 𝑠 + 𝛽 𝑠 ′ Σ -2 𝑠 ′ ) * (¬M 𝑠 * ¬M 𝑠 ′ ) 6: 𝚺 2 𝑠𝑠 ′ ← (𝛽 𝑠 𝚺 -2 𝑠 + 𝛽 𝑠 ′ 𝚺 -2 𝑠 ′ ) -1 7: 𝑛 𝑠𝑠 ′ ← 𝑛 𝑠 + 𝑛 𝑠 ′ 8: M 𝑠𝑠 ′ ← M 𝑠 ∨ M 𝑠 ′ 9: return Y 𝑠𝑠 ′ , M 𝑠𝑠 ′ , 𝑛 𝑠𝑠 ′ , 𝚺 2 𝑠𝑠 ′
and 𝑦 ′ 𝑠 (𝑥) at 𝑥 result from actual sensor measurements or from interpolation, as known from the masks M 𝑠 and M 𝑠 ′ :

• Line 2 -The two values result from actual sensor measurements: a merged incremental average is applied. • Lines 3 & 4 -One value results from actual sensor measurements and the other is inferred: the actual sensor measurement is considered to be the ground truth and thus it is selected over the inference. • Line 5 -The two values result from two inferences: the aggregated value is then computed using the Generalized Productof-Expert of GPR, as detailed below.

Aggregating Multiple GPR Inferences

The Generalized Product-of-Expert is a method that allows combining estimated that have been inferred by several experts (e.g., interpolations on several crowdsensors). In particular, it enables weighting the respective importance of the experts according to the reliability of their inference. Let 𝑝 𝑠 (𝑦 * |𝑥 * , 𝑋 𝑠 , 𝑌 𝑠 ) denote the distribution of the measurements for point 𝑥 * , which is inferred by the crowdsensor 𝑠, knowing the observed cell values 𝑌 𝑠 at points 𝑋 𝑠 .

Assuming that 𝑚 crowdsensors aggregate their inference results, the Product-of-Expert for a GPR estimates a value 𝑦 * at point 𝑥 * according to the following joint distribution [START_REF] Peter | Distributed Gaussian processes[END_REF]:

𝑝 (𝑦 * |𝑥 * , 𝑋, 𝑌 ) = 𝑚 𝑠=1 𝑝 𝛽 𝑠 (𝑥 * ) 𝑠 (𝑦 * |𝑥 * , 𝑋 𝑠 , 𝑌 𝑠 )
where 𝛽 𝑠 is a weighting parameter that allows tuning the relative importance of a crowdsensor 𝑠 according to its inference. The aggregation of multiple GPR inferences is a generalized Product-of-Expert, which accounts for multiple inference distributions 𝑝 𝑠 of an arbitrary point 𝑥 * . According to [START_REF] Cao | Generalized product of experts for automatic and principled fusion of Gaussian process predictions[END_REF], it combines many Gaussian distributions with mean μ𝑠 (𝑥 * ) and variance σ2 𝑠 (𝑥 * ) from any crowdsensor 𝑠, and the aggregation result is defined as:

μ (𝑥 * ) = σ2 (𝑥 * ) 𝑚 𝑠=1 𝛽 𝑠 (𝑥 * ) σ-2 𝑠 (𝑥 * ) μ𝑠 (𝑥 * ) σ2 (𝑥 * ) = [ 𝑚 𝑠=1 𝛽 𝑠 (𝑥 * ) σ-2 𝑠 (𝑥 * )] -1
A generalized Product-of-Expert of a GPR allows merging several inferences (i.e., 𝑚 = 2 inferences in our case) in a cost-effective way, as it is characterized by a O (𝑛) time complexity with 𝑛 being the number of merged points. The opportunistic aggregation on the move is asymmetric, as captured by the loss function 𝐷. Our evaluation shows that assigning a greater 𝛽 to the crowdsensor acting as the mainstay (i.e., resulting in the lesser loss function) leads to a higher aggregation accuracy (see § 5).

IAM-BASED CROWDSENSING SYSTEM DESIGN

The IAM solution allows mapping quantitative spatio-temporal physical phenomena through opportunistic data interpolation and aggregation across the participating mobile crowdsensors. Any mobile crowdsensing application dealing with urban environmental monitoring (e.g., noise level, air quality, temperature, etc.) may build upon IAM to produce such a knowledge in a fully decentralized way. Still, the actual D2D meeting of the contributing crowdsensors within the area under monitoring depends on the size of the target area and the density of the contributing crowd. In practice, the monitoring of large areas requires running an ultimate aggregation process at a server to connect the islets that the crowd covers. Figure 1 illustrates the resulting operation of a IAM-based crowdsensing system: crowdsensors sense, pre-process and interpolate the data, then they relay and aggregate the sensing data in a P2P way using wireless D2D communication (e.g., Wi-Fi Direct or Bluetooth technologies), so as to favor ubiquitous computation and thereby limit the dependence on the server infrastructure. At the end of the monitoring period D, the remaining mainstays send their tensors to the server, which composes the tensors. Focusing on the illustrative figure: C aggregates data from B and D, and then uploads to the server the results of the three aggregations across A, B, D and itself. Similarly, the distant mainstay E, which aggregates its contributions with the ones of H, G and F, provides the resulting tensor to the server. At the end of the current time window, the spatio-temporal interpolation is applied on the associated local tensor to infer missing values (Line 6 -See detail in § 3.2). We highlight that the interpolation runs locally, only once and prior to the aggregation process running over the next time window. This allows: (i) minimizing the number of interpolation occurrences and thereby the resource cost on the device, and (ii) leveraging the locally completed tensor to assess the quality of the local measurements against the ones of the peers that the node meets, which determines the mainstay.

Finally, still at the end of the current time window 𝑇 +D , the node sends its local tensor to the server, unless it was relayed to another crowdsensor (Lines 7-9).

Prototype Implementation

The IAM solution assumes an ad hoc framework supporting the opportunistic meeting and collaboration of peer nodes using D2D communication, such as our middleware presented in [START_REF] Du | Let opportunistic crowdsensors work together for resource-efficient, quality-aware observations[END_REF] or other solutions in [START_REF] Higuchi | A neighbor collaboration mechanism for mobile crowd sensing in opportunistic networks[END_REF][START_REF] Lee | Comon: Cooperative ambience monitoring platform with continuity and benefit awareness[END_REF][START_REF] Liu | Scents: Collaborative sensing in proximity iot networks[END_REF]. The IAM prototype is thus specifically focused on the implementation of the distributed, collaborative interpolation and aggregation, further targeting Android smartphones/tablets as crowdsensors.

The prototype is available at https://github.com/sensetogether/ IAM; it requires a Python 3 environment as well as the following third-party packages for data analytics: NumPy (https://numpy.org) handling multi-dimensional arrays (tensors) and Scikit-learn (https: //scikit-learn.org) implementing various machine learning algorithms that is used for GPR training/inference. The key components of our prototype implementation are:

• Pre-processing reads the sensed data that is stored in a local file, and creates the corresponding tensors Y and M. The tensor size is an application-specific parameter that is configured so as to map the target physical phenomenon P over the chosen A and D. Precisely, the tensor size depends on the required sensing resolution, the geographical space that needs to be covered, and the time window.

EVALUATION

IAM supports the opportunistic aggregation of the sensing data along with the interpolation of a physical phenomenon, provided relevant measurements from the crowdsensors. We assess the effectiveness of the IAM distributed approach using a dataset collected with the Ambiciti crowdsensing application (available on Google Play and App Store) for environmental monitoring.

Experiment Setup and Dataset

The experiment supporting our evaluation focuses on urban noise monitoring using crowdsensing. The dataset contains the measurements collected by the Ambiciti application in an opportunistic way [START_REF] Issarny | Dos and don'ts in mobile phone sensing middleware: Learning from a large-scale experiment[END_REF]. The partner company Ambiciti provided us the dataset used for evaluation, which specifically relates to the contributions gathered in Paris over year 2017 and includes about 950k entries from 550 crowdsensors. The dataset entries are tuples of the form: <device ID, latitude index, longitude index, timestamp index, observation value> with the observation value being the average sound level expressed in dB(A). Note that a sound level in dB(A) is a logarithmic quantity and hence sound levels cannot be simply averaged. Instead, the sound levels in dB(A) are first converted into their energy equivalents, and then the energy equivalents are averaged algebraically, and finally the resulted energy equivalent is converted back to its dB(A) value.

Preprocess

Interpolate Aggregate 2: Producing a noise map using IAM (zoomed scale)

In the above context, IAM manages tensors that deal with the monitoring of the noise level over the whole area A of the city during D = 24 hours. We decompose the city area into a 100 × 100 grid, and the sensing data that are collected during 1 hour are stored in the dedicated cell. Over 1 day, this results in a 100×100×24 tensor, which contains at most 240𝑘 of data entries. Figure 2 illustrates the data analysis procedure for the computation of the urban noise map, using the IAM prototype. The figure zooms on the 5 × 5 × 1 snapshot for which the crowdsensors provide sensor measurements. This highlights that the dataset is very sparse, which is the case with most crowdsensing applications and our reference application in particular.

We use the Root Mean Square Error (RMSE) to evaluate the accuracy of the interpolated and aggregated tensors. For crossvalidation, we run 100 rounds of training followed by tests. At each round, both the training and evaluation sets are randomly shuffled, that is, for interpolation, 70% of the dataset is used for training (i.e., as actual observations to complete the tensor) and 30% to test (i.e., to assess the estimated values against the ground truth). Regarding aggregation, the entire approximation tensor is used for the evaluation.

In the following evaluation, the experiments are run either on a DELL Precision 7520 workstation, used both as the centralized server ( § 5.2) and for simulation ( § 5.3), or on Android smartphones as end device testbed ( § 5.4). The D2D communication is assumed to be supported by the WiFi-Direct protocol, and we use the distance to detect P2P encounters [START_REF] Du | Let opportunistic crowdsensors work together for resource-efficient, quality-aware observations[END_REF].

Assessment of the Interpolation Methods

Figure 3 compares the robustness of the three inference strategies that are commonly used to interpolate physical phenomena (see § 2.2): Ordinary Kriging with Gaussian variogram model (OK-Gaussian), CP decomposition with Alternating Least Square (CP-ALS), and Gaussian Process Regression with Matern kernel (GPR-Matern). The interpolation is performed at the server (without involving any aggregation), using the dataset from which we selected the day during which the largest amount of crowdsensing data were collected. The same experiments were run using the whole dataset and the same trends were observed. On the figures, the box corresponds to the interquartile range, the orange line is the median and the green triangle is the mean. At first sight, OK-Gaussian seems to be accurate and hence promising, given the low RMSE interquartile range and median. However, some wrong inferences lead to abnormal values as illustrated by the high RMSE mean value of 47. Similarly but to a lower extent, CP-ALS shows some abnormal RMSE. Instead, GPR-Matern provides both an accurate and robust inference: a stable RMSE without outliers -hence characterized by the lowest variance-is observed.

Figure 4 shows the execution time of the three interpolation approaches. We run the experiments over the 365 days of our dataset, where the number of observed cells varies every day. The execution time of CP-ALS is constant in O (𝑅𝐼 𝐽 𝐾) regardless of the number of available observed cells since the computation applies on the entire fixed-size tensor. Both OK-Gaussian and GPR-Matern have a time complexity in O (𝑛 3 ) with 𝑛 being the number of observed cells used to fit the model. The execution time of GPR-Matern is lower than OK-Gaussian, and below CP-ALS when the number of observed cells is less than 2800. Note that in our dataset, the number of observed cells collected by crowdsensors per day remains lower than 1500. We further evaluated the efficiency of the three interpolation methods in terms of memory consumption. GPR-Matern consumes the least memory: around 3.114𝑀𝐵, with a variance of 1.718. While the memory consumption associated with CP-ALS (resp. OK-Gaussian) is of 3.258𝑀𝐵 with a variance of 0.422 (resp. 4.644𝑀𝐵 with a variance of 1.437). Note that the memory consumption is stable and does not depend on the number of observed cells since our approach always uses a fixed-size tensor that is filled with zeros in the absence of observations. Overall, GPR-Matern is the most efficient in terms of accuracy and robustness, while its execution time is also relatively lower.

Focusing on GPR, we assess the accuracy of the following kernels in Figure 5: constant, RBF, rational quadratic, and Matérn. The rational quadratic and Matérn kernels are the most accurate, while the former slightly outperforms the latter. However, the execution time of the quadratic kernel is twice as much as that of the Matérn kernel. We therefore leverage GPR with Matern kernel within IAM.

Distributed vs Centralized Aggregation

We compare the overall performance of the distributed vs centralized interpolation-mediated aggregation using our 1-year dataset. In the centralized aggregation case, the server collects all the sensing data and performs the aggregation and interpolation based on the whole dataset. As for the distributed aggregation (see Figure 2), we consider the following methods:

• Ideal iterative aggregation is a theoretical and sequential case in which the aggregation starts at the first crowdsensor that aggregates its tensor with the next crowdsensor and the aggregation process repeats with the following crowdsensors until the last crowdsensor is reached. This is the ideal case for which we ignore the actual locations of the crowdsensors. • Base stochastic aggregation represents the real-life scenario: an aggregation occurs when at least two crowdsensors meet, as detected using the actual location and time proximity available from the dataset. The aggregation process thus depends on the mobility of the contributing users. Upon a meeting, the mainstay is selected randomly and 𝛽 = 𝛽 ′ = 1 for the generalized Product-of-Expert (see Algorithm 2). Ultimately, all the tensors are uploaded to and merged at the server, either directly or via a relay depending on the crowdsensors' P2P meetings. • IAM opportunistic aggregation is similar to the above stochastic aggregation with the exception of the selection of the mainstay and the chosen 𝛽 values. It follows our Algorithms 1 & 2, and we set 𝛽 𝑠 = 1.5 (resp. 𝛽 𝑠 ′ = 0.5) for crowdsensor 𝑠 with lower 𝐷 (resp. 𝑠 ′ with higher 𝐷).

As illustrated in Figure 6, the execution time associated with the centralized aggregation (and interpolation) at the server significantly increases when the number of crowdsensors gets high. Instead, when the interpolation and aggregation are mainly performed by crowdsensors, the server execution time is almost negligible regardless of the number of crowdsensors. In addition, the storage requirement is minimized on the server because the data tensor size is always unchanged when aggregating new incoming data (via linear operations).

We further assess the benefit of the decentralized IAM approach to crowdsensing from a financial perspective. We specifically compare the IAM solution with the more classical centralized one, which often relies on a cloud platform for data analysis and storage. This is in particular the configuration of the crowdsensing system of the Ambiciti company that provided us the dataset: the system initially used the Google Cloud Platform (GCP, https: //cloud.google.com/products/), which we consider as an illustrative candidate for estimating the budget associated with a cloud-based configuration. In a nutshell, the price depends on the amount of network traffic, the amount of storage needed, and the load associated with the computation (e.g., execution time associated with the analysis). Figure 7 estimates the monthly financial cost associated with running centralized vs distributed crowdsensing while using the GCP platform as the server and assuming that each crowdsensor sends a 2𝑀𝐵 packet every day. The cost associated with the centralized approach increases linearly because the number of uploads and the computation involved to interpolate the phenomena are both high. Instead, the costs of the stochastic and opportunistic approaches remain low because communication toward the cloud is reduced and only lightweight aggregation is performed on the cloud. Our IAM opportunistic aggregation outperforms the stochastic approach because we further reduce the computing load on the cloud. In general, the IAM decentralized solution significantly alleviates the dependence on the infrastructure server, thereby enabling the wider adoption of crowdsensing systems by communities of users concerned with environmental monitoring. Figure 8 provides the RMSE of the three distributed aggregations compared to the centralized approach that serves as reference. Although the ideal iterative aggregation avoids the processing of interpolation and aggregation on the cloud, the RMSE mean (resp. median) equals 2.213 (resp. 2.345). The accuracy of the stochastic aggregation is quite similar, with a RMSE mean of 2.212, and a RMSE median of 2.341. Our opportunistic aggregation performs better than the other distributed approaches with a RMSE mean (resp. median) of 2.102 (resp. 2.225). Still, as expected, the decentralization impacts on the overall aggregation result, which is to be compared to the resulting resource gains. It is part of our future work to investigate further enhancement of the distributed interpolationmediated aggregation by, e.g., accounting for the significance of the measurements gathered at a node when interpolating.

IAM Impact on the Device Resources

We now focus on the resource consumption of IAM on the end device, for which we analyze the execution time (depending on the number of observed cells and of aggregations) and power consumption (depending on the execution time and D2D protocol). We empirically assess the performance associated with the IAM prototype in terms of execution time and energy consumption, using Android smartphones.

Figure 9 compares the amount of traffic uploaded to the server in the centralized vs distributed cases. The traffic is evaluated based on the number of actual P2P aggregations (within relays). As expected, the distributed aggregation reduces the amount of traffic uploaded to the server and hence the cellular network occupancy is kept to a minimum. Furthermore, the IAM opportunistic aggregation drastically reduces the uploading to the server by 54.2% compared to the stochastic aggregation. A portion of the traffic sent to the server is replaced by the D2D forwarding among crowdsensors; Then, we run experiments on a SAMSUNG GALAXY S7 smartphone embedding a 3000 mAh battery capacity. Figure 10 shows the interpolation execution time depending on the number of cells (i.e., area at the target scale) covered by the crowdsensor, and the aggregation execution time depending on the number of aggregations. Note that the figure shows no more than 500 entries, which is in practice a very high number of cells contributed/visited by one crowdsensor per day. As expected, the interpolation is computationally intensive compared to the aggregation, whose execution time is comparatively negligible: the interpolation takes a couple of minutes when the number of observed cells is greater than 500, while the aggregation takes less than 100 seconds for a number of aggregations below 500 and for a number of observed cells per crowdsensor varying from 1 to 500. The aggregation shows a linear time complexity.

Recall that each crowdsensor executes the interpolation only once, which is the most computation-intensive operation. Assuming the crowdsensor has 500 observed cells, an interpolation consumes the most energy with 88mAh, and an aggregation consumes only 2mAh. Figure 11 estimates the energy consumed by a smartphone that implements the interpolation and opportunistic aggregation when the P2P meeting frequency varies: the more frequent is the aggregation, the more energy is consumed. Nevertheless, the related energy consumption remains under control because in practice the crowdsensor usually has already relayed/aggregated its data before encountering around 8 crowdsensors for a single day. In summary, the interpolation of 500 observed cells and 8 aggregations for a day cost only 2.88% of the battery capacity. With respect to communications, based on the power assessment in [START_REF] Kalic | Energy consumption in android phones when using wireless communication technologies[END_REF], uploading via cellular network consumes 8.9 (resp. 4) times of the D2D relay energy via Bluetooth (resp. Wi-Fi).

DISCUSSION

To the best of our knowledge, our work is the first to investigate and design a fully distributed interpolation and aggregation based on the opportunistic encounters of crowdsensors. Our key hypothesis is that the number of encounters -both past and future-is correlated to the number of observations and thus the inference quality.

Other common criteria such as device resource/status may be further considered when selecting the mainstay. Still, our evaluation does not analyze the power consumption due to D2D communication, as we rely on the results of previous studies that show that D2D networking is cheaper than cellular networking [START_REF] Xu Chen | Exploiting massive D2D collaboration for energy-efficient mobile edge computing[END_REF][START_REF] Do | Crowdsourced mobile data transfer with delay bound[END_REF][START_REF] Wang | ecoSense: Minimize participants' total 3G data cost in mobile crowdsensing using opportunistic relays[END_REF]. It is part of our future work to further investigate the overall effectiveness of IAM in terms of resource efficiency vs data accuracy, compared to the centralized approach. In addition to the required energy efficiency and accuracy for any crowdsensing system, ensuring privacy is another key concerns for the end-users. Here, we claim that the fully decentralized approach of IAM outperforms the centralized approach in terms of privacy, further considering the opportunistic, impromptu encounters of the crowdsensors, which subsequently share aggregated & interpolated tensors about environmental phenomena. The opportunistic approach then raises the potential issue of un-trustworthy crowdsensors that may contribute malicious data. While the IAM solution presented in the paper assumes trustworthy crowdsensors providing equally accurate measurements, our aim in the near future is to investigate mechanisms that filter anomalous data (e.g., outliers) to deal with untrustworthy contributors.

CONCLUSION

One of the major benefits of crowdsensing is the possibility to monitor environmental phenomena at the urban scale, simply leveraging the abundance and capacity of people's smartphones. In practice, the people's mobility makes the crowdsensing contribution unevenly distributed over space and time, which requires the analysis of the contributed observations that is in general performed at a central, often cloud-based, server. However, as the number of contributors grows, the increasing number of observations that the crowdsensing systems must process gets challenging: the high network and financial cost associated to a cloud-centric system hinders the widespread deployment of crowdsensing, and the high computational cost due to the large amount of data makes intractable the modeling of the environmental phenomenon.

We tackle the above issues by exploiting the increasing computing capacity of today's smartphones, that is, we distribute the interpolation and aggregation associated with the sensing data at the powerful end devices. To do so, we introduce IAM that runs on the smartphone to capture complex relationships among the collected observations across both space and time by relying on Gaussian Process Regression and 3D tensors. Then, the resulting tensors are opportunistically combined together following a stochastic process based on the physical encounters of people. The benefit of our approach is threefold: (i) each crowdsensor (i.e., expert) independently establishes an interpolation of the region it covered; (ii) the aggregation resulting from the Product-of-Experts is sharper than any of the individual tensor and renders much more tractable the establishment of the overall tensor; and (iii) the computation achieved on the device is limited, and thus not energyexhausting. Indeed, the evaluation using a real-world dataset shows that our approach significantly reduces the transmission to, and the computing resource consumed on, the infrastructure server, compared to the centralized approach.

Figure 1 : 3 3 :

 133 Figure 1: A IAM-based crowdsensing system

•

  Interpolation creates a GPR model, trains the model based on the M * Y observed cells, and uses the trained model to infer and produce the approximation tensor Ŷ along with the variance tensor Σ 2 . The interpolation has a O (𝑛 3 ) time complexity with 𝑛 being the number of observed cells.• Aggregation computes the loss function of two tensors from a pair of crowdsensors, makes the aggregation decision, then merges the two tensors into one, and updates the current data, mask and variance tensors Ŷ, M and Σ 2 , respectively, following Algorithm 1. The aggregation process on each crowdsensor has a O (𝑝) time complexity, where 𝑝 is the number of P2P meetings.

  Figure 3: Interpolation accuracy -RMSE

Figure 6 :

 6 Figure 6: Execution time at the server only

Figure 7 :

 7 Figure 7: Monthly financial cost of a cloud-based deployment

Figure 8 :

 8 Figure 8: Aggregation accuracy -RMSE

Figure 11 :

 11 Figure 9: Directly uploaded messages

Y

  ¬M 𝑠 ∥ 2 where: ¬ corresponds to the NOT Boolean operation; and M 𝑠 and M 𝑠 ′ correspond to the mask tensors of 𝑠 and 𝑠 ′ . Require: local data tensor Y 𝑠 , local mask tensor M 𝑠 , local variance tensor Σ 2 𝑠 , local merge count 𝑛 𝑠 Input: remote data tensor Y 𝑠 ′ , remote mask tensor M 𝑠 ′ , remote variance tensor Σ 2 𝑠 ′ , remote merge count 𝑛 𝑠 ′ 1: if 𝐷 (Y 𝑠 , Y 𝑠 ′ ) < 𝐷 (Y 𝑠 ′ , Y 𝑠 ) then 𝑠 , M 𝑠 , 𝑛 𝑠 , Σ 2 𝑠 ← Crowdsensing data tensors merger (𝑠, 𝑠 ′ )

	Algorithm 1
	Asymmetric P2P aggregation at 𝑠
	2:
	-see Algorithm 2
	3: else
	4:

  Crowdsensing data tensors merger (𝑠, 𝑠 ′ ) at the mainstay 𝑠 Input: data tensor Y 𝑠 , mask tensor M 𝑠 , variance tensor Σ 2 𝑠 , merge count 𝑛 𝑠 Input: data tensor Y 𝑠 ′ , mask tensor M 𝑠 ′ , variance tensor Σ 2 𝑠 ′ , merge count 𝑛 𝑠 ′ Output: aggregated data tensor Y 𝑠𝑠 ′ , mask tensor M 𝑠𝑠 ′ , variance tensor Σ 2 𝑠𝑠 ′ , merge count 𝑛 𝑠𝑠 ′ 1:

	Algorithm 2


Thanh Trung Nguyen

Gérard Meunier

Jean-Michel Guichon

Olivier Chadebec

T S Nguyen

T-T Nguyen

J-M Guichon

T-S Nguyen

An Integral Formulation for the Computation of 3D Eddy Current Using Facet Elements

Keywords: Eddy currents, Facet element, Generalized PEEC method, Integral equation method

d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

HE PEEC method (Partial Element Equivalent Circuit) is a well-known integral equation technique leading to an equivalent circuit representation of a electromagnetic device. It is mainly used for the modeling of complex interconnect problems and can be favorably applied to devices where the air region is dominant [START_REF] Ruehli | Equivalent circuit models for three dimensional multiconductor systems[END_REF]. However, the classical PEEC method does not enable easily the modeling of 3D conductive media. Moreover, it requires the generation of a structured mesh associated with uniform current density on each element. Thus, it use is still limited for general geometries.

Previously, eddy currents integral formulations based on general Whitney elements have been proposed [START_REF] Albanese | Integral formulation for 3D eddy current computation using edge elements[END_REF] [START_REF] Albert | Comparision of a direct and a vector potential integral equation method for the computation of eddy currents[END_REF]. Our approach has the particularity of being based on the use of facet elements [START_REF] Bossavit | Whitney forms: a class of finite elements for three dimensional computations in electromagnetism[END_REF]. The reliability of such elements has already been shown for finite element magnetostatic formulations [START_REF] Le Menach | Determination and utilization of the source field in 3-D magnetostatic problems[END_REF], for integral magnetostatic formulation [START_REF] Vishnevsky | Modeling magnetization using Whitney facet elements[END_REF] and also for steady conduction problems [START_REF] Vishnevsky | Boundary integral computation of electric fields in multielectrode galvanic systems using normally continuous elements[END_REF]. Different quantities have been approximated, such as the magnetization in [START_REF] Vishnevsky | Modeling magnetization using Whitney facet elements[END_REF] and the current distribution in [3][5].

This paper proposes a generalization of the classical inductive PEEC formulation using facet elements enabling the use of general mesh. A similar approach has been recently proposed by Alotto et al. but limited to the modeling of surface regions [START_REF] Alotto | Dual-PEEC modeling of a two-port TEM cell for VHF applications[END_REF]. In our approach, the current density is linearly interpolated with first order facet elements. An equivalent electrical circuit representation, whose branches are the facets and the nodes of these branches are the centroids of elements, is proposed (i.e. an equivalent electrical circuit build on the dual mesh). The integral inductive PEEC formulation is then adapted to this configuration. The solenoidality of the current distribution is ensured thanks an independent loops search technique.

II. FORMULATION

A. PEEC Integral equation

From Maxwell-Faradays equation, a magnetic vector A and a scalar electric potential U exist such as:

U j grad A E     ( 1
)
where E is the electric field and ω the pulsation. We consider the conductive and non-magnetic region ΩC representing a set on conductors. The electric field E in the conductor satisfies Ohm's law

U j grad A J      (2)
with σ is the material conductivity and J is the current density. Without any magnetic material, we can write the magnetic vector potential A as

C d r C     J A   4 0 ( 3
)
where μ0 is the vacuum permeability. By introducing this expression in (2), we get

U d r j C C grad J J          4 0 (4)
The current density has to satisfy the several conditions. In the conductive region ΩC, div J = 0 [START_REF] Vishnevsky | Boundary integral computation of electric fields in multielectrode galvanic systems using normally continuous elements[END_REF] and on the boundary of ΩC,

n.J = 0 ( 6
)
where n is the external normal.

B. Facet Element Interpolation

The current density is interpolates with first order facet elements such as

  j j j I w J (7)
T where wj is shape function and Ij is flux across the jth facet. Expressions for wj can be easily computed for different kind of element [START_REF] Bossavit | Whitney forms: a class of finite elements for three dimensional computations in electromagnetism[END_REF]. The main property of such element family is that the normal component of wj is conserved through each facet, so also ensuring the conservation of current. Let us remember some the properties for these shape functions:

j j s 1 .   n w ( 8
) e j v 1 div   w ( 9
)
where sj is surface of facet j and ve is volume of element e which contains the facet j. The sign (  ) depends on the facet orientation.

C. System Assembly Using Facet Elements

Applying the Galerkin method to (4) by using wi as projection functions, a system of linear equations is obtained

        b f f b U I L R I Z    jω ( 10
)
with

C j i ij d σ C     w w . R ( 11
) C C j i ij d d r π μ C C       w w . 4 L 0 ( 12
) C i bi d U C      grad w . U ( 13
    f f f f Ud s 1 U ( 15
)
with ve is the volume of element e and sf is surface of facet f. Let us apply divergence theorem on equation ( 13), we get

         C C i i bi d U d U ) ( div ) . ( U w n w ( 16
)
where Г is the boundary of region ΩC.

Let us now consider the computation of Ubi for an internal facet i (see Fig. 2.). The first term in ( 16) is null because the value of wi vanishes on all others facets of both adjacent elements sharing facet i. Thus, from ( 9) and ( 14), we can deduce that Ubi is the difference between the averaged voltages of both elements sharing i.

Let us now consider Ubi for a border facet (i.e. belonging to the boundary of region ΩC). From ( 8) and ( 15) the value of the first term in (16) becomes equal to the averaged potential on the facet. The second term becomes equal to the volume averaged voltage on the only elements to whom facet i belongs. Of course, facet orientations have to be considered in order to properly compute the values of Ubi. From previous considerations, an equivalent electrical circuit (10) can be generated. The branches of this circuit are represented by the facets of initial mesh. Each element of the mesh can then be seen as a node of this circuit. This electrical circuit is in fact an equivalent representation of the dual mesh.

D. Resolution of the Electrical Circuit

In order to ensure the solenoidality of J, we can use a similar technique to which used for classical structured PEEC method. In our approach, we have decided to use a circuit solver resolution based on the determination of independent loops. The general principle consist in find small topological loops (minimum number of branches), allowing a very quick research of independent loops. Furthermore, the choice to maximize the number of small loops allows significantly improving the matrix condition number of the final system of equations.

In addition, we can notice that the boundary condition ( 5) can be easily taken into account in the formulation, by simply suppressing the associated degree of freedom in the equivalent circuit. If the meshed conductor is electrically coupled with an external circuit, the electric coupling can naturally be imposed in the formulation without any difficulty.

Fundamental circuit equations are expressed in the presence of external source voltages.

   0   s b U U M (16)
with [M] is the branch-fundamental independent loop transition matrix where the value of each element can be -1, 0 or 1 and Us the vector of external voltage sources (most part of time equal to 0). We can write a new system of linear equations where unknowns Im are currents flowing in independent loops.

      S m t b U M I M Z M   ( 17
)
Once linear system has been solved, we obtain the currents flowing on each branch by applying the following equation.

E. Addition external source field

In presence of an external source like (a coil for instance) in which the current density is imposed, equation (4) has to be adapted.

U j d r j C C grad A J J        0 0 4     (19)
where A0 is the source magnetic vector potential created by the external source. Equation [START_REF] Phung | Automatic cuts for magnetic scalar potential formulation[END_REF] becomes

  0 b f b U U I Z   (20)
where

     C C i i d j 0 0 . U A w  (21)
Independent loops search technique and resolution is then applied to (20) as previously.

III. NUMERICAL EXAMPLE

In order to validate the proposed formulation and to show its performances, we consider the Problem 7 of the TEAM workshop benchmark problems [START_REF] Turner | Workshops and problems benchmarking eddy current codes[END_REF]. The results will be compared to those obtained with FLUX [START_REF] Cedrat | [END_REF], a commercial Finite Element Method (FEM) program.

The geometry of the Problem 7 is presented in Fig. 3. The problem consists in an asymmetrical conductor with a hole and an exciting coil. The coil is placed 30mm above the plate and excited with 2472 Ampere-turns at 50Hz and 200Hz. The driving current reaches the maximum at a phase of 0 degrees. In this problem, we have to add an additional source magnetic vector potential representing the contribution of the coils in which the current distribution is known. We have

0 0 0 0 d r 4 0     J A   (22)
where, J0 is the current density vector of the coil and Ω0 is coil region. A0 is evaluated numerically thank a tetrahedral discretization of the coil and a Gauss point quadrature. We focus on the computed eddy current distribution and losses in the plate (Fig. 4(a), Fig. 4(b) and Tab. 1.). The results show that the relative error of the formulation in comparison with the 3D FEM associated to a very high mesh density is good (<5 %). Figure 5 Another advantage of the new formulation is its capacity to treat any multiply connected problems thanks to the efficient determination of independent loops of the equivalent electrical circuit. This point is really valuable for the modeling of industrial problem, where the techniques of automatic cuts search can be poorly efficient [START_REF] Phung | Automatic cuts for magnetic scalar potential formulation[END_REF]. In this paper, we have presented an original eddy currents volume integral formulation using facet elements dedicated to the modeling of non magnetic conductors. The formulation has a small relative error in comparison with a 3D reference FEM with a high mesh density. Moreover, the convergence of the results is quickly reached and results with very coarse meshes remain acceptable.

This new formulation seems very attractive because it enables an easy treatment of the multiply connected problems without cuts technique and without meshing the air. Moreover, the coupling with the classical PEEC formulation is natural, both formulation having same theoretical bases. It opens interesting perspective for the modeling for multiply connected volume/surface regions with direct electrical connections a complex wiring system.

The obtained matrix is fully dense but this problem can be overcome by the use of a compression technique like Fast Multipole Method [START_REF] Nguyen | Ships magnetic anomaly computation with integral equation and fast multipole method[END_REF] enabling the save of memory and reducing the computation time. It remains to study the robustness of this approach in this context where convergences can sometimes be difficult to reach.

Fig. 1 .

 1 Fig.1. Representation of shape functions of the faces {i, j, k} and {i, j, k, l} for the reference tetrahedral and reference hexahedra.

  ) Matrix [Zb] can be seen as the impedance matrix of the electrical equivalent circuit generated, If is vector of currents through the facets, [R] is matrix of the resistive terms, and [L] is the matrix of the mutual inductances. Let us notice that the resistive matrix [R] is sparse and similar to a matrix obtained during a classical FEM assembly. On the other hand, the inductive part [L] is fully dense and more representative of integral method-based assembly. To avoid integral singularity inaccuracies while the computation of [L], different Gauss points repartition are considered especially for the computation of the self inductance (i.e. Lii).Let us define Ue the average values of the voltage on each finite volume element Ωe and Uf its average value on each face Γf of the mesh. We have:

Fig. 2 .

 2 Fig. 2. Primal and dual meshes, the black points are the centroid of elements, red ones are the centroid of the faces on the boundary. Both Ubi configurations are represented (internal and border facets).

Fig. 3 .

 3 Fig. 3. Geometry of TEAM problem 7: coils and conducting plate with hole (all dimension are in mm).

  (a) and 5(b) present the convergence of the Joules losses value of two methods follows the number of elements in the plate at 50 Hz and 200 Hz. Results provided by our formulation are very encouraging, the convergence being reached with a very few number of elements in comparison with the FEM one.

Fig. 4 .Fig. 5 .

 45 Eddy current distribution in the plate, (a) at 50 Hz, (b) at 200 Convergence of Joules losses value in the plate, (a) at 50 Hz, (b) at 200 Hz IV. CONCLUSION