Data Quality Audit and Bivariate Temperature-Humidity Analysis of a Long-Term LoRaWAN Network: Anomaly Detection, PCA, Clustering and Machine Learning

Data Quality Audit and Bivariate Temperature–Humidity Analysis of a Long-Term LoRaWAN Network: Anomaly Detection, PCA, Clustering and Machine Learning arXiv:2610.32438v1 [cs.LG] 23 Oct 2026

Mamadou SOW LIPN-CNRS7030, Université Sorbonne Paris Nord 99 avenue JB Clément, F-93430 Villetaneuse, France mamadou.sow@univ-paris13.fr

Abstract—Long-term environmental IoT/LPWAN deployments often aggregate heterogeneous data streams from diverse devices without systematic auditing of each variable’s actual coverage prior to analysis—an increased risk when multiple sensor types coexist within the same network. This work proposes a C++17 processing pipeline applied to the Tour Perret LoRaWAN dataset (421,937 frames, June 2021–June 2023, 5 devices). A systematic coverage audit integrated into the extraction phase reveals that CO2 is completely missing (0% coverage) and pressure is reported only by an isolated base module (Wyres, 4.24% of the corpus), invalidating the initial hypothesis of a joint four-variable analysis. The analysis is refocused on the Temperature–Humidity pair, jointly available across 205,005 frames (95.76%) from 4 Elsys EMS sensors. The pipeline sequentially performs anomaly detection (modified Zscore, Isolation Forest), bivariate Principal Component Analysis, clustering (K-means, DBSCAN), and a multilayer perceptron predicting humidity from temperature and time, whose residuals serve as a detector for slow drift. On real-world data: Temperature– Humidity correlation 𝑟 = −0.662, two-component PCA (PC1 = 83.1% of variance), Isolation Forest detects 2.12% multivariate anomalies while univariate Z-score detects none, calibrated DBSCAN isolates 0.17% noise across 13 clusters, and the MLP achieves 𝑅 2 = 0.672 (RMSE = 10.54%) in cross-validation. This work demonstrates that per-variable coverage auditing prior to modeling is a necessary methodological step for heterogeneous multi-sensor IoT datasets. Index Terms—LoRaWAN, IoT, data quality audit, anomaly detection, Principal Component Analysis, clustering, machine learning, C++, urban microclimate, heterogeneous sensors.

processing pipelines expecting all target variables—such as Temperature (𝑇), Relative Humidity (𝐻), Atmospheric Pressure (𝑃), and Carbon Dioxide concentration (CO2 )—to be present synchronously in every transmitted frame. However, in multisensor LPWAN deployments, distinct physical sensors often operate on different transmission schedules, report varying payload schemas across firmware revisions, or belong to completely separate device applications aggregated under a single network server instance. B. Central Finding: Systematic Data Coverage Audit

To prevent silent pipeline failures and catastrophic data distortion (such as invalid outer joins producing null-dense matrices or artificial cross-device frame alignment), an explicit, automated Data Coverage Audit was embedded directly into the primary C++ extraction module (01_extract.cpp). The audit evaluated payload completeness across all registered device identifiers (devEUI) and application boundaries (applicationName). The coverage audit revealed two critical structural properties of the raw dataset: 1) Total Absence of Carbon Dioxide (CO2 ): Despite initial project documentation indicating CO2 monitoring for indoor environmental quality, zero frames across the multiyear corpus contained decoded CO2 payload fields. This stems either from a physical deployment omitting the I. Introduction & Coverage Audit CO2 gas sensor module or an unupdated network decoder A. Context and Operational Challenges instance. Long-term environmental monitoring via Low-Power Wide2) Structural Isolation of Barometric Pressure (𝑃): Area Networks (LPWAN)—specifically using the LoRaWAN Pressure measurements originated exclusively from a protocol—presents unique data engineering and data quality secondary device cohort (Wyres sensors), whereas the challenges. In real-world deployments such as the Perret Tower primary sensor fleet consisted of four Elsys EMS (Environ(Tour Perret) dataset [1], raw telemetry streams suffer from mental Monitoring Sensor) nodes. Because the Wyres and network packet loss, heterogeneous payload formats, multiElsys EMS devices transmitted asynchronously without gateway frame duplication, and hardware configuration shifts clock synchronization, barometric pressure was structurally over time. disconnected from the main stream at the individual frame A common pitfall in multivariate anomaly detection and level. predictive modeling is the a priori assumption of variable co-occurrence. Researchers and practitioners often design

Raw LoRaWAN Telemetry Stream (Aggregated JSON-Lines Log Files)

1. Coverage Audit Evaluate schema completeness by Variable & devEUI

Excluded Variables

Retained Bivariate Core Dataset

CO2 : Totally absent (0%) Pressure: Structurally isolated (Wyres stream)

Elsys EMS frames (95.76%) Co-occurring T & H 𝑛 = 205,005 complete frames

2. Preprocessing & Cleaning Deduplication (devEui | fCount) Physical bounds filtering

Fig. 1. LoRaWAN telemetry ingestion and cleaning pipeline: the coverage audit (Step 1) is executed before any downstream analysis, separating structurally unusable variables (CO2 , pressure) from the retained bivariate (T,H) core dataset.

C. The Bivariate Temperature–Humidity Dataset Attempting an inner join across all four initial variables (𝑇, 𝐻, 𝑃, CO2 ) would yield 0 usable rows, causing a complete execution failure. Conversely, performing a naive cross-device alignment would artificially merge environmental readings from physically distinct locations within the tower structure. Consequently, the scope of the study was refocused on the largest, physically consistent, co-occurring payload schema: the bivariate Temperature–Humidity (𝑇, 𝐻) feature space provided by the four Elsys EMS nodes. Following extraction, ISO 8601 timestamp parsing to Unix epoch, multi-gateway frame deduplication (keyed on devEui|fCount), and physical boundary validation (𝑇 ∈ [−40, 85] °C, 𝐻 ∈ [0, 100]%), the clean corpus comprised 𝑛 = 205,005 complete bivariate frames, representing 95.76% of the retained frames telemetry corpus. This rigorous audit ensures that all subsequent statistical, dimensionality reduction, clustering, and machine learning models operate on true, synchronized physical observations. II. Related Work The processing and analysis of long-term Internet of Things (IoT) data streams lie at the intersection of streaming anomaly detection, unsupervised subspace learning, spatial clustering, and supervised drift monitoring. This section reviews the theoretical and empirical literature foundational to our pipeline architecture. A. Anomaly Detection in Edge Telemetry & Sensor Networks Anomaly detection in streaming sensor networks has been extensively surveyed by [2] and [3]. Traditional point-based methods rely on statistical distributional assumptions, such as the univariate Z-score. However, as noted by [4], standard

Z-scores based on sample mean and standard deviation break down in the presence of extreme spikes due to their low breakdown point; robust estimators like the Median Absolute Deviation (MAD) are required for non-Gaussian IoT telemetry. For multivariate streams, [5] introduced the Isolation Forest algorithm, which isolates anomalies by randomly partitioning feature space using axis-aligned trees. Unlike distance or density estimators, Isolation Forest exhibits near-linear time complexity 𝑂 (𝑛 log 𝑛), making it particularly effective for high-volume edge streams. More recently, [6] emphasized the necessity of real-time edge processing architectures capable of detecting anomalies under memory and bandwidth constraints, reinforcing the importance of low-overhead C++ implementations over heavy Python runtime environments. B. Subspace Methods & Dimensionality Reduction Principal Component Analysis (PCA), established by [7] and generalized by [8], remains the standard technique for modeling global linear dependencies among continuous physical variables. In multivariate environmental sensing, PCA transforms correlated measurements into orthogonal variance components. When deploying linear algebra algorithms directly on constrained IoT gateways or embedded systems, external libraries (e.g., LAPACK, Armadillo) introduce overhead and portability barriers. Golub and Van Loan [9] outline numerical diagonalization algorithms, specifically the Jacobi rotation method for symmetric matrices. For bivariate correlation matrices (2 × 2), closed-form algebraic solutions bypass iterative numerical solvers entirely, providing exact eigenvalue decomposition with zero iterative computational overhead.

C. Partitioning and Density-Based Clustering Cluster analysis serves a dual purpose in environmental telemetry: segmenting distinct climatic regimes and isolating unstructured spatial noise. • Partitioning Methods: The classical K-means algorithm [10] minimizes intra-cluster variance but is notoriously sensitive to initial centroid placement. [11] addressed this with K-means++, introducing probabilistic distance-weighted initialization that guarantees an 𝑂 (log 𝐾) approximation bound. To determine the optimal number of clusters 𝐾, [12] proposed the silhouette index, which quantifies cluster cohesion against separation. • Density-Based Methods: Unlike K-means, DBSCAN [13] discovers clusters of arbitrary shapes and explicitly identifies noise points (labeled −1). However, naive DBSCAN scales at 𝑂 (𝑛2 ) complexity and suffers from memory allocation limits (std::bad_alloc) when expanding dense spatial queues on large datasets (𝑛 > 200,000). Mitigating this requires spatial indexing structures—such as uniform grid bins or 𝑘-d trees—and objective hyperparameter calibration via 𝑘-distance plot analysis. D. Predictive Modeling & Drift Identification Supervised neural networks, specifically Multilayer Perceptrons (MLP) trained via backpropagation [14], are widely used to approximate nonlinear physical couplings between environmental variables. In training deep or recurrent networks, [15] demonstrated that gradient clipping and momentum-based stochastic gradient descent (SGD) stabilize convergence against exploding gradients. While point-based univariate and multivariate detectors (Zscore, Isolation Forest, DBSCAN) evaluate each telemetry frame independently, they cannot identify progressive calibration loss or sensor degradation over multi-month horizons. As highlighted by [3], cross-validating a physical predictive model (e.g., predicting relative humidity from temperature and cyclical time encodings) allows large regression residuals (|𝑦 𝑖 − 𝑦ˆ 𝑖 | > 3×RMSE) to serve as a complementary indicator for slow, systemic sensor drifts that remain within static physical bounds. E. Positioning Relative to Existing Multi-Method Pipelines and Data-Quality Auditing Practice

individual techniques available for outlier detection, missingdata imputation, and virtual sensing in environmental and air-quality monitoring networks, no framework in their survey unifies coverage auditing with cross-validated, complementary anomaly detection at the pipeline level — most deployments instead apply a single technique per data-quality dimension. This gap is precisely what the coverage-audit-first design proposed in Section IV-A, combined with the four-tier detection framework of Section III-B, is intended to address: rather than treating variable-coverage verification as a preliminary, disposable step, this work integrates it directly into the extraction pipeline and demonstrates empirically (Table II) that skipping it would have silently invalidated the entire downstream analysis. To our knowledge, no prior work specifically targets LoRaWAN deployments with an integrated coverage-audit-and-multi-tierdetection pipeline of this kind; positioning this contribution against the broader (but still comparatively sparse) IoT dataquality-auditing literature, rather than solely against singletechnique anomaly detectors, is left as recommended future work for a revised version of this paper. III. Theoretical Framework & Methodology A. The Need for Multi-Tiered Anomaly Detection in IoT Streams Anomalies in environmental time-series data do not manifest as a single homogeneous class. In edge-computing and IoT sensor contexts, anomalies span a spectrum of behaviors: • Extreme Univariate Outliers: Sensor spikes, electrical bursts, or out-of-spec hardware failures. • Multivariate Structural Anomalies: Observations where each individual feature remains within normal seasonal bounds, but their joint combination violates fundamental physical laws or dominant environmental correlations. • Density-Based Clustered Anomalies / Noise: Isolated noise points residing in low-density regions of the multivariate feature space. • Systemic Temporal Drifts: Gradual, continuous calibration loss (e.g., aging capacitive humidity elements) that evades point-based detection because individual readings shift slowly over weeks or months. No single detection algorithm can capture all four anomaly types simultaneously. Therefore, this work establishes a fourtier complementary detection framework, combining univariate statistics, random partitioning, geometric subspace projection, and supervised cross-validation.

Sections II-A–II-D each address anomaly detection, subspace projection, clustering, or drift monitoring as largely separate bodies of literature. Multi-method pipelines that combine sev- B. Mathematical Formulation of the Complementary Frameeral of these detectors into a single, complementary framework work 1) Tier 1(a): Point-Based Univariate Detection — Modified specifically for LPWAN/IoT telemetry remain comparatively Z-Score: Classical 𝑍-scores rely on the sample mean (𝑥) ¯ uncommon: most cited works (e.g., [5]; [13]) evaluate one and standard deviation (𝜎), both of which are notoriously detection paradigm in isolation, and surveys of IoT sensornetwork data quality tend to treat missing-data handling, outlier sensitive to extreme outliers (low breakdown point). To ensure detection, and drift correction as distinct tasks addressed by robust estimation over multi-year climate cycles, we employ the disjoint bodies of work rather than a single orthogonal, four-tier Modified Z-Score based on the Median Absolute Deviation architecture of the kind proposed here. [16], in a recent review (MAD) [4]: of data-quality applications for IoT monitoring sensor networks, 𝑥 𝑖 − 𝑥˜ 𝑀𝑖 = 0.6745 · (1) make a closely related observation: despite the wide range of MAD

Clean Bivariate Dataset (T, H) 𝑛 = 205,005

TIER 1(a) Modified Z-Score (Univariate)

TIER 1(b) Isolation Forest (Multivariate)

TIER 2 Bivariate PCA Reconstruction (Subspace 2D)

TIER 3 Cross-Val MLP Residuals (Temporal/Drift)

Identifies extreme singlevariable bounds ( | 𝑀𝑖 | > 3.5)

Identifies abnormal joint feature combinations (𝑠 > 0.62)

Quantifies departure from dominant T-H anticorrelation (𝑒𝑖 = | 𝑦PC2 | )

Identifies progressive calibration sensor drift ( | res | > 3× RMSE)

Fig. 2. Four-tier complementary detection framework: each tier targets a distinct, orthogonal anomaly signature within the clean bivariate (T,H) dataset.

where 𝑥˜ = median(𝑥) and MAD = median𝑖 (|𝑥 𝑖 − 𝑥|). ˜ The constant 0.6745 equates the MAD to the standard deviation under a normal distribution. A point is flagged as a marginal univariate outlier if |𝑀𝑖 | > 3.5. 2) Tier 1(b): Point-Based Multivariate Isolation —Isolation Forest: Univariate statistics cannot detect anomalous feature couplings. The Isolation Forest algorithm [5] exploits the principle that anomalous observations are few and different, making them susceptible to isolation via random recursive axis-aligned hyperplanes. An ensemble of 𝑡 isolation trees partitions the standardized 2D space (𝑇, 𝐻). The anomaly score 𝑠(𝑥, 𝑛) for an observation 𝑥 over a dataset of size 𝑛 is defined as: 𝐸 [ℎ( 𝑥) ]

𝑠(𝑥, 𝑛) = 2− 𝑐 (𝑛)

(2)

where 𝐸 [ℎ(𝑥)] is the average path length across all isolation trees, and 𝑐(𝑛) is the average path length of unsuccessful searches in a Binary Search Tree (BST): 2(𝑛 − 1) 𝑐(𝑛) = 2𝐻 (𝑛 − 1) − , 𝑛 𝐻 (𝑘) = ln(𝑘) + 𝛾 (𝛾 ≈ 0.5772)

𝑒 𝑖 = 𝑥𝑖′ − 𝑉1𝑉1⊤ 𝑥𝑖′ 2

In a 2D orthonormal basis, this reduces to the absolute magnitude of the rejected second principal component: 𝑒 𝑖 = |𝑦 𝑖,PC2 |. Observations with high 𝑒 𝑖 represent severe departures from the physical correlation structure of the atmosphere (e.g., simultaneous high temperature and high humidity). 4) Tier 3: Temporal & Drift Anomaly Detection — Supervised MLP Cross-Validation: Point-based methods (Modified Z-score, Isolation Forest, PCA) evaluate each frame as an isolated tuple, remaining blind to slow sensor drifts where values drift gradually while staying inside valid static bounds. To capture temporal inconsistencies, a Multilayer Perceptron (MLP) is trained to predict 𝐻 using Temperature 𝑇 alongside cyclically encoded temporal features: h xtime = sin(2𝜋 · hour/24), cos(2𝜋 · hour/24), sin(2𝜋 · day/365.25), cos(2𝜋 · day/365.25)

(3)

As 𝑠(𝑥, 𝑛) → 1, the path length ℎ(𝑥) → 0, indicating an easily isolated multivariate anomaly. We enforce an empirical isolation threshold of 𝑠(𝑥, 𝑛) > 0.62. 3) Tier 2: Geometric Subspace Anomaly Detection — Bivariate PCA Reconstruction Error: Principal Component Analysis (PCA) models the global covariance structure. In a standardized 1 2D space (𝑋 ′ ), the covariance matrix 𝑆 = 𝑛−1 𝑋 ′⊤ 𝑋 ′ is diagonalized to yield eigenvectors 𝑉1 , 𝑉2 and eigenvalues 𝜆1 ≥ 𝜆2 . Retaining only the primary component (PC1, representing the dominant thermo-hygrometric opposition axis), the projection matrix is 𝑉1 . The multivariate reconstruction error 𝑒 𝑖 is given by:

(4)

i (5)

The network architecture (16 hidden units, tanh activation) minimizes the Mean Squared Error: 𝑛

𝐿=

1 ∑︁ (𝑦 𝑖 − 𝑦ˆ 𝑖 ) 2 𝑛 𝑖=1

(6)

An observation 𝑖 that is not flagged by Tier 1 or Tier 2, yet exhibits a large prediction residual |𝑦 𝑖 − 𝑦ˆ 𝑖 | > 3× RMSE, serves as a prime indicator of localized sensor drift or unmodeled environmental micro-disturbances. C. Methodological Orthogonality Summary The theoretical necessity for this multi-tiered architecture is summarized in the matrix below:

TABLE I Methodological orthogonality of the four detection tiers. Detection Tier

Underlying Principle

Primary Target

Modified Z-Score

Marginal distance from median, in units of MAD

Extreme single-variable spikes (𝑇 or 𝐻 )

Insensitive To

Multivariate correlation violations within normal marginal bounds Slow, continuous Isolation Forest Random axis-aligned Unusual joint partitioning path combinations of calibration drifts length (𝑇 , 𝐻 ) over long horizons PCA Reconstruction Distance to the Departure from T–H Anomalies aligned primary variance anti-correlation along the primary subspace (𝑒𝑖 = | 𝑦PC2 | ) variance axis (PC1) DBSCAN Noise Local density Sparse spatial Uniform dense Label threshold outliers, isolated clusters from sensor ( 𝑁 𝜀 ( 𝑝) < minPts) noise points quantization MLP Deviation from Slow sensor drift, Sudden point anomalies uncoupled Cross-Validation learned temporal from physical/temporal micro-climate diurnal/seasonal dependency anomalies trends

IV. Data Processing & Analytical Pipeline A. Data Acquisition, Coverage Audit & Cleaning

choice is what makes the coverage audit below possible in the first place — it lets per-variable and per-devEUI presence be measured across the entire corpus without first discarding rows for missing columns. During extraction, an automated Data Coverage Audit evaluates payload completeness across all registered device identifiers (devEUI) and application boundaries (applicationName) in a single pass. TABLE II Coverage per variable and per device (audit on the 421,937 raw frames). Variable

Raw Frames

Retained

Temperature 393,685 214,076 (𝑇) Humidity (𝐻) 384,614 205,005 Pressure (𝑃) 9,071 0 (isolated) Carbon 0 0 (absent) Dioxide (CO2 )

Cov. Cov. Scope (%Tot.) (%Ret.) 93.30 91.15 2.15 0.00

100.00 Elsys + Wyres 95.76 Elsys only 0.00 Wyres only 0.00 —

Nodes 5 4 1 0

The coverage audit confirms, on the full real-world corpus, what the early four-variable extraction attempt had already suggested empirically: a static extraction assuming joint co-occurrence across all four initially planned variables (𝑇, 𝐻, 𝑃, CO2 ) yields zero rows by construction (complete-case filtering on an empty intersection). This audit — integrated into the extractor and executed in a single pass over the file — is itself the methodological step whose necessity this work argues for: without it, the four-variable hypothesis would have failed silently rather than being diagnosed and corrected. Distribution by devEUI (214,076 retained frames): 4 Elsys EMS sensors (108,678, 35,387, 30,800, and 30,140 frames respectively, totaling 205,005 joint T+H frames) and 1 Wyres module (9,071 T+P frames, no humidity). 2) Multi-Gateway Deduplication & Physical Bounds Filtering: In LoRaWAN deployments, a single wireless transmission frame is frequently picked up by multiple surrounding gateways. Frame deduplication is performed using a composite hash key combining the unique device identifier and the frame counter:

The primary data acquisition and preprocessing layer is implemented in 01_extract.cpp. The raw corpus comprises 421,937 JSON-Lines frames transmitted by five distinct LoRaWAN nodes deployed at the Perret Tower (Tour Perret) in Grenoble, France, between June 2021 and June 2023, decoded by a ChirpStack network server [1] (DOI 10.18709/perscido.2023.06.ds395). 1) JSON Schema Robustness & Coverage Audit: Rather than querying rigid, static path keys, 01_extract.cpp executes a recursive, alias-tolerant search across each frame’s JSON payload hierarchy (object.* for physical measurements, rxInfo[] for reception metadata). This handles schema variations across ChirpStack v3/v4, The Things Network (TTN) v3, and flat-format payloads without throwing runtime keyerrors. Two schema-level pitfalls, identified by diagnostic inspection of the raw file prior to full extraction, required specific handling: (i) the timestamp field is neither time nor timestamp but appears as either _date (an ISO 8601 Deduplication Key = devEui ∥ fCount (7) string) or _timestamp (a Unix epoch in milliseconds, Deduplication eliminated 179,609 duplicate frame receipts not seconds — a common source of silent date-arithmetic (42.6% of total ingested logs), reflecting the dense gateway errors if left unconverted); and (ii) Elsys EMS frames in coverage across the Grenoble urban area. Following deduplicathis deployment do not always carry all four target variables tion, physical boundary validation was applied based on the jointly — the first frames inspected during development carried Elsys EMS manufacturer specification sheet [17]: only temperature and humidity (alongside unrelated fields such as dewpoint, accMotion, pulseAbs, vdd, • Temperature: 𝑇 ∈ [−40.0, 85.0]°C waterleak, x/y/z). An early extraction attempt requiring • Relative Humidity: 𝐻 ∈ [0.0, 100.0]% all four variables in the same frame therefore rejected the entire Only three individual out-of-range sensor readings were corpus, including frames where 𝑇 and 𝐻 were both perfectly discarded. The resulting retained core dataset consists of valid — motivating the per-variable extraction strategy below. 𝑛 = 205,005 synchronized, bivariate (𝑇, 𝐻) frames from four Each of the four target variables (𝑇, 𝐻, 𝑃, CO2 ) is conse- Elsys EMS sensors (108,678, 35,387, 30,800, and 30,140 quently extracted independently: a frame is retained as soon as frames respectively), representing 95.76% of the retained it carries at least one target variable and a valid timestamp, with deduplicated corpus. absent columns left empty in a sparse (“long”) intermediate CSV rather than requiring a strict wide-format join. This design

tourperret.log(.gz) Raw LoRaWAN Logs

01_extract.cpp JSON tree search, dedup, physical bounds filter

clean_data.csv 𝑛 = 205,005 bivariate frames

02_anomaly_detection.cpp Z-Score & Isolation Forest

05_mlp_regression.cpp Cyclical encodings, cross-val

03_pca.cpp Closed-form 2 × 2 PCA, 𝑒𝑖

anomaly_report.csv

mlp_report.csv Residuals = 4th signal

pca_scores.csv PC1, PC2, 𝑒𝑖

04_clustering.cpp K-means++ & Grid-DBSCAN

clusters.csv

Fig. 3. Data processing and analytical pipeline: from raw LoRaWAN logs to the four output reports (anomaly, MLP, PCA/cluster) that feed the Results section.

C. Dimensionality Reduction: Closed-Form Bivariate PCA The module 03_pca.cpp performs the Tier-2 subspace The C++ module 02_anomaly_detection.cpp implements the three point/density-based tiers of the complementary projection introduced in Section III-B (reconstruction error 𝑒 𝑖 , detection framework introduced in Section III-B, refocused on Eq. (4)), and derives its exact closed-form solution for the the bivariate (𝑇, 𝐻) space jointly available across the 205,005 bivariate case. 1) Standardization & Closed-Form 2 × 2 Diagonalization: Elsys EMS frames — down from the four dimensions initially Given raw features with incompatible physical units (°C and planned. 1) (a) Point-Based Univariate Outliers: Robust Modified %), features are standardized to zero mean and unit variance: Z-Score: The Modified Z-score (Eq. (1), Section III-B) is 𝑥𝑖′ 𝑗 = (𝑥𝑖 𝑗 − 𝑥¯ 𝑗 )/𝜎 𝑗 . The 2 × 2 sample correlation matrix 𝑆 is computed independently for 𝑇 and 𝐻; observations with |𝑀𝑖 | > expressed as:   3.5 are flagged as univariate marginal outliers. 1 1 𝑟 𝑇 ,𝐻 2) (b) Multivariate Point Outliers: Isolation Forest: The 𝑆= 𝑋 ′⊤ 𝑋 ′ = (8) 𝑟 𝑇 ,𝐻 1 𝑛−1 Isolation Forest anomaly score (Eqs. (2)–(3), Section III-B) is While general multivariate PCA (𝑑 ≥ 3) requires numerical evaluated on standardized bivariate vectors x𝑖 = [𝑇𝑖′ , 𝐻𝑖′ ] ⊤ . The solvers such as the Jacobi rotation algorithm [9], a 2 × 2 forest is instantiated with 𝑡 = 100 isolation trees, each built symmetric matrix admits an exact closed-form algebraic on a subsample of 𝜓 = 256 observations, with maximum tree solution for its eigenvalues 𝜆 1 , 𝜆2 : depth capped at ℎmax = ⌈log2 (256)⌉ = 8. As in Section III-B, √︁ an empirical decision threshold of 𝑠(x, 𝑛) > 0.62 is enforced. tr(𝑆) 2 − 4 det(𝑆) tr(𝑆) ± 3) (c) Global multivariate structure: DBSCAN: DBSCAN 𝜆1,2 = = 1 ± |𝑟 𝑇 ,𝐻 | (9) 2 (04_clustering.cpp, Section IV-D and Section V-D) The corresponding orthonormal eigenvectors v1 , v2 define contributes a third, orthogonal detector: points located in lowthe principal subspace projection matrix 𝑉 = [v1 | v2 ]: density regions of the PCA space (label −1) reflect the global multivariate structure rather than a per-variable score or random    1 sgn(𝑟 𝑇 ,𝐻 ) 1 −sgn(𝑟 𝑇 ,𝐻 ) tree isolation path. v1 = √ , v2 = √ (10) 1 1 2 2 B. Anomaly Detection Framework

Substituting v1 into the reconstruction error of Eq. (4) (Section III-B) and exploiting the 2D orthonormal basis, 𝑒 𝑖 reduces directly to the absolute magnitude of the rejected PC2 coordinate, 𝑒 𝑖 = |𝑦 𝑖,PC2 | — this closed form is what 03_pca.cpp evaluates for every observation. Observations with high 𝑒 𝑖 represent severe departures from the dominant physical correlation structure of the atmosphere. D. Clustering Algorithms & Real-World Scaling Optimizations The module 04_clustering.cpp implements Kmeans++ and DBSCAN for spatial partitioning and density analysis. 1) K-means++ Inertia Minimization & Silhouette Optimization: K-means minimizes intra-cluster variance (inertia) 𝐽: 𝐽=

𝐾 ∑︁ ∑︁

x − 𝝁𝑘 2

(11)

𝑘=1 x∈𝐶𝑘

Using Lloyd’s algorithm [10] with K-means++ seeding [11], initial cluster centers are sampled with probability proportional to their squared distance to the nearest existing 2 center: 𝑃(x) = Í𝐷𝐷(x) . To select the optimal number of clusters (x) 2 𝐾, the silhouette index 𝑠(𝑖) [12] is computed across candidate evaluations 𝐾 ∈ {2, 3, . . . , 8}: 𝑠(𝑖) =

𝑏(𝑖) − 𝑎(𝑖) max(𝑎(𝑖), 𝑏(𝑖))

(12)

where 𝑎(𝑖) represents mean intra-cluster distance and 𝑏(𝑖) denotes mean nearest-neighbor cluster distance. 2) DBSCAN Spatial Indexing & Algorithmic Scaling (𝑛 = 205,005): DBSCAN [13] defines clusters as continuous highdensity regions separated by low-density space. A point p is a core point if its 𝜀-neighborhood 𝑁 𝜀 (p) contains at least minPts points: 𝑁 𝜀 (p) = { q ∈ D | ∥p − q∥ ≤ 𝜀 }, subject to |𝑁 𝜀 (p)| ≥ minPts (13) Points not density-reachable from any core point are flagged as noise (label −1).

Fig. 4. Uniform grid spatial index (2D), cell size = 𝜀: DBSCAN’s 𝜀neighborhood query is restricted to a 3 × 3 neighborhood of grid cells around the query point’s cell, avoiding an 𝑂 (𝑛2 ) full scan.

3) Real-World Scaling Challenges: Executing naive 𝑂 (𝑛2 ) spatial range queries on 𝑛 = 205,005 bivariate points leads to two computational failures: excessive execution times and memory allocation faults (std::bad_alloc) caused by unchecked expansion queues in high-density regions. To resolve these limitations, 04_clustering.cpp incorporates two structural optimizations: 1. Uniform Grid Spatial Indexing: Standardized PCA feature space is partitioned into a uniform 2D grid with square bin width equal to 𝜀. Range queries for point p = (𝑥, 𝑦) are restricted to candidate points residing within its home grid cell (𝑖, 𝑗) and the 8 immediate neighboring cells. Under quasi-uniform spatial distributions, search complexity scales from 𝑂 (𝑛2 ) down to 𝑂 (𝑛). 2. Expansion Queue State Tracking: A boolean visited array (std::vector<bool> in_queue) tracks point assignments during region expansion, preventing identical points from being repeatedly queued and eliminating runaway memory allocation. 3. Hyperparameter Calibration via 𝑘-Distance Plot Rather than using arbitrary default parameters (e.g., 𝜀 = 0.35), parameter selection is driven by 𝑘-distance analysis with 𝑘 = minPts = 10. For every point in the 𝑛 = 205,005 dataset, the Euclidean distance to its 10th-nearest neighbor was calculated and sorted in ascending order. The maximum curvature point (“knee”) of the resulting curve occurs at 𝜀 ≈ 0.071. Setting 𝜀 = 0.35 proved too coarse, merging 100% of points into a single cluster and returning zero noise points.

serve as an independent fourth anomaly indicator targeting progressive sensor drift and micro-climatic disruptions. F. Synthetic Validation Set: Generation Procedure and Role in the Study Design Before being applied to the real-world Tour Perret corpus, every module of the pipeline (Sections IV-B–IV-E) was first exercised against a synthetic bivariate Temperature– Humidity dataset (𝑛 = 20,000 samples) generated with a known, near-deterministic thermo-hygrometric anti-correlation (target 𝑟 ≈ −0.93). This synthetic set was not designed as a formal experimental baseline to be reported as a standalone result; rather, it served as a pre-deployment calibration and software-verification harness, distinct from the numerical cross-check already described for 03_pca.cpp — in which Fig. 5. 𝑘-distance plot (𝑘 = minPts = 10) used to objectively select the the C++ eigen-decomposition is independently recomputed in DBSCAN neighborhood radius 𝜀 on real-world data. Candidate thresholds are overlaid; the knee of the curve, identified via maximum perpendicular distance Python (numpy.linalg.eigh) on the real data to confirm to the chord, is retained at 𝜀 ≈ 0.071. bit-for-bit agreement on PC1 = 83.10%, PC2 = 16.90% (Section IV-C). The synthetic set, by contrast, was used earlier in the pipeline’s development, before the idiosyncrasies of the E. Cross-Validation via Supervised Learning (MLP) The module 05_mlp_regression.cpp implements a real corpus (variable sparsity, actual point density) were fully Multilayer Perceptron (MLP) to model physical dependencies characterized. It played three concrete roles: 1. Default-parameter calibration and its failure mode between variables and detect progressive sensor drift. (Section IV-D). DBSCAN’s initial density threshold was 1) Target Revision & Physical Rationale: In the initial calibrated on the synthetic set at 𝜀 = 0.35. Applied project design, the MLP was specified to predict CO2 conunchanged to the real corpus (𝑛 = 205,005, substantially centration from thermal and pressure covariates. However, the denser in standardized PCA space than the synthetic data coverage audit (Section IV-A) proved that CO2 is entirely sample), this value proved far too permissive: it collapsed absent from the multi-year corpus (0% presence). the entire dataset into a single cluster with 0% noise — To preserve the predictive cross-validation architecture, the an uninformative result that made clustering structure target variable was revised to predict Relative Humidity impossible to assess. This failure is itself an informative (𝐻) using Temperature (𝑇) and temporal features. Relative finding: it demonstrates that DBSCAN hyperparameters Humidity is especially prone to hardware degradation and calibrated on a synthetic proxy do not transfer to the real capacitive element saturation under prolonged humid conditions, deployment, and it directly motivated the data-driven 𝑘making it an ideal target for cross-validation drift analysis. distance recalibration described in Section IV-D, which 2) Cyclical Feature Encoding & Architecture: Time comyielded the retained value 𝜀 ≈ 0.071. ponents are cyclically encoded exactly as defined in Eq. (5) 2. MLP negative control (Section V-E). For the regression (Section III-B), preventing artificial boundary discontinuities at cross-validation, the synthetic humidity target was addidaily (23h → 00h) and annual (Dec 31 → Jan 1) transitions. ⊤ tionally corrupted with deliberately strong multiplicative The input vector x𝑖 = 𝑇𝑖′ , xtime ∈ R5 feeds a feedforward noise, calibrated so that no learnable (𝑇, time) → 𝐻 network comprising one hidden layer with 16 neurons (tanh relationship would remain. The MLP trained on this coractivation) and a single linear output unit predicting 𝐻ˆ 𝑖 . rupted synthetic target achieved 𝑅 2 = 0.019 — statistically Network parameters are optimized using stochastic gradient indistinguishable from a non-informative predictor. This descent (SGD) with momentum (𝜇 = 0.9, learning rate 𝜂 = negative-control result is the reference against which the 0.002) and gradient norm clipping [15], minimizing the MSE real-world performance (𝑅 2 = 0.672, Section V-E) should loss of Eq. (6) (Section III-B): be read: it shows that the architecture does not spuriously report a high 𝑅 2 in the absence of genuine signal, which v𝑡+1 = 𝜇v𝑡 − 𝜂∇𝐿, 𝑊𝑡+1 = 𝑊𝑡 + v𝑡+1 (14) in turn supports treating the real-world 𝑅 2 = 0.672 as 3) High Residuals as Long-Term Drift Signals: Point-based evidence of a real, learnable thermo-hygrometric and detectors evaluate individual frames as independent events, diurnal/seasonal dependency rather than an artifact of making them blind to slow calibration drifts where readings model capacity or overfitting. remain within static physical bounds. By cross-validating 3. Visualization scale reference (Figure 11). The synthetic physical dependencies through the neural network, observations sample size (𝑛 = 20,000) is markedly smaller and less exhibiting large regression residuals: dense than the real corpus (𝑛 = 205,005); marker size and transparency in the PCA projection scatter plots were |𝑦 𝑖 − 𝑦ˆ 𝑖 | > 3 × RMSE (15) adjusted accordingly (smaller, more transparent markers)

once the pipeline moved from the synthetic development set to the real, denser point cloud. Interpretation. Read together, the correlation and 𝑅 2 comparisons in Section V-B and Section V-E are not a same-footing benchmark between two datasets of equal experimental status: the synthetic set encodes a clean, controllable ground truth used to verify that each detector behaves as expected (or, in the case of DBSCAN’s initial 𝜀, to expose a calibration that does not transfer) before the same code is trusted on field data. The lower real-world correlation magnitude (|𝑟 𝑇 ,𝐻 | = 0.662 vs. 0.93 on the synthetic set) reflects genuine microclimatic variability, weather events, and sensor noise absent from the synthetic generator (Section V-B); the higher real-world 𝑅 2 (0.672 vs. 0.019 on the noise-corrupted synthetic set) confirms the MLP is picking up a genuine physical signal rather than an artifact of the fitting procedure (Section V-E). V. Results A. Descriptive Statistics (𝑛 = 205,005, Complete Elsys EMS T+H Frames) TABLE III Statistical summary of the 𝑛 = 205,005 synchronized Elsys EMS bivariate telemetry corpus. Variable

Mean

Std. Dev.

Min

Q1 (25%)

Median

Q3 (75%)

Max

Temperature 𝑇 (°C) Humidity 𝐻 (%)

14.36 72.14

8.72 18.41

−4.30 19.00

7.00 58.00

14.80 74.00

20.70 88.00

46.10 100.00

1) Physical Interpretation & Measurement Artifacts: • Thermal Extremes (𝑇max = 46.10°C): The maximum recorded temperature exceeds ambient climate norms for Grenoble. This reflects localized solar irradiance heating on the unshaded exterior housing of the rooftop node. • Humidity Saturation (𝐻max = 100.00%): Relative humidity exhibits strong left-skewness, with a third quartile of Q3 = 88.00% and frequent saturation at 100.00%. This behavior corresponds to capacitive sensor saturation during sustained rainfall or high atmospheric humidity. B. Bivariate PCA Analysis The Pearson correlation coefficient between standardized Temperature and Humidity is 𝑟 𝑇 ,𝐻 = −0.662. Through the characteristics of Figure 7, this confirms a strong inverse thermo-hygrometric relationship, consistent with the physical property that warmer air exhibits higher saturation vapor capacity.

Fig. 6. Bivariate PCA eigenvector axes in standardized space. Variance PC1 = 83.10%, PC2 = 16.90%. Points are a schematic reconstruction of the standardized (𝑇 ′ , 𝐻 ′ ) cloud preserving the reported correlation 𝑟 = −0.662; √ the plotted axes are exact (±1/ 2 loadings, Eq. (10)).

TABLE IV Eigenvalue decomposition of the 2 × 2 correlation matrix (cf. Figure 8). Component Eigenvalue Explained Cumulative 𝑇 Loading 𝐻 Loading Var. (%) 𝜆𝑗 Var. (%) PC1 PC2

1.6620 0.3380

83.10% 16.90%

83.10% 100.00%

−0.7071 −0.7071

+0.7071 −0.7071

Loadings: T (PC1 = −0.707, PC2 = −0.707), H (PC1 = +0.707, PC2 = −0.707)—PC1 represents, as expected from a strongly correlated bivariate PCA, the classic thermohygrometric opposition axis (bisector), while PC2 acts as its orthogonal complement (sum T+H, which can be interpreted as a global “thermo-hygrometric energy” indicator for the sensor system). 1) Principal Axis Physical Meanings: • PC1 (Thermo-Hygrometric Opposition Axis): Captures 83.10% of total variance, aligning with the line of inverse correlation (𝐻 ′ ≈ −𝑇 ′ ). • PC2 (Thermo-Hygrometric Energy Axis): Captures the remaining 16.90% of variance along the orthogonal bisector (𝐻 ′ ≈ 𝑇 ′ ). Elevated values along PC2 indicate high simultaneous temperature and humidity, signaling severe departures from normal atmospheric behavior. Figure 7 (correlation): 𝑟 𝑇 ,𝐻 = −0.662—a moderate to strong anti-correlation, consistent with expected physics (air capacity to hold water vapor increases with temperature), yet noticeably noisier than on the synthetic validation set (𝑟 = −0.93, see Section IV-F for its generation procedure and role in the study design). This is expected: real-world data incorporates microclimatic variability, weather events, and measurement noise that a synthetic dataset cannot reproduce.

2. Multivariate Detection: The Isolation Forest successfully flags 4,346 frames where temperature and humidity violate their expected inverse coupling (e.g., high humidity occurring at elevated ambient temperatures), despite both variables remaining within valid single-variable ranges.

Fig. 7. Actual Pearson correlation matrix, Temperature–Humidity (𝑛 = 205,005).

Fig. 9. Schematic reconstruction of the temperature/humidity time series for the Perret Tower (2021–2023) with Isolation Forest-flagged anomalies overlaid (Modified Z-score set to zero anomalies on this dataset). Curve shape preserves the reported seasonal range and anomaly rate (2.12%); exact daily values are not reproduced from the original data.

Fig. 8. Scree plot of eigenvalues (bivariate PCA: T, H, real data). PC1 alone captures 83.10% of total variance.

C. Anomaly Detection: Demonstration of Complementarity on Real-World Data TABLE V Comparative performance across the primary point-based anomaly detectors (𝑛 = 205,005).

Anomaly Detector

Primary Target Anomaly Type

Flagged (𝑛) Rate (%)

Modified Z-Score (|𝑀𝑖 | > 3.5)

Extreme single-variable spikes Atypical multivariate combinations Simultaneous univariate & multivariate

0

0.00

4,346

2.12

0

0.00

Isolation Forest (𝑠 > 0.62) Joint Consensus (both methods)

1) Methodological Insights: The contrast between the univariate Modified Z-score (0%) and the multivariate Isolation Forest (2.12%) highlights the necessity of multi-tiered detection architectures: 1. Univariate Robustness: Even extreme single-variable observations (𝑇max = 46.10°C) remain below the statistical outlier threshold (|𝑀 | > 3.5) due to the wide, natural multi-year climate spread.

This result—zero in univariate, non-zero in multivariate—is the clearest empirical demonstration of the methodological complementarity discussed in Section II [3]: the marginal extremes for T (max = 46.1 °C, |𝑀 | = 3.10) and H (min = 19%, |𝑀 | = 2.47) remain, when considered in isolation, below the significance threshold (3.5) given the natural spread across two full seasonal cycles—yet 2.12% of the (T, H) combinations are statistically atypical under multivariate isolation. This reveals unusual thermo-hygrometric couplings (such as high temperatures combined with high humidity, contradicting the dominant anti-correlation) that remain invisible to singlevariable analysis. D. Cluster Structure Analysis Silhouette-based selection of 𝐾 (Table VI) clearly favors 𝐾 = 2 (silhouette = 0.491), with the silhouette score dropping sharply starting at 𝐾 = 3 (0.404) and then stabilizing around 0.34–0.39 for 𝐾 ≥ 4—a sign that the dominant structure of the validation data is a binary opposition between normal regime / extreme events rather than a fine-grained segmentation into multiple climate regimes (the latter being expected on realworld data, where seasonal and diurnal cyclicity should yield a 3–4 cluster structure, see Section VI).

1% for humidity). These coarse measurement steps produce identical repeated values, creating artificially dense point stacks in low-density boundary regions. DBSCAN (𝜀 = 0.071 determined via k-distance plot, Figure 5, minPts=10): 13 clusters, 339 noise points (0.17%). Figures 10–11b reveal a massive main cluster and 12 microclusters localized in the region PC1 ≈ −1, PC2 ≈ −1.5 to −2.5—corresponding to exact duplicate values resulting from sensor quantization (temperature to the nearest 0.1 °C, humidity to the nearest integer), which create strictly identical points in sufficient numbers to locally form their own dense clusters in the DBSCAN sense. This phenomenon, distinct from a genuine climate structure, is discussed in Section VI. Fig. 10. DBSCAN density clustering in standardized PCA space (𝜀 = 0.071, minPts= 10), accelerated by the uniform grid spatial index of Figure 4. Point cloud is a schematic reconstruction preserving the reported cluster/noise proportions. TABLE VI Evaluating K-means++ clustering across 𝐾 ∈ {2, . . . , 8} using the silhouette index.

Cluster Count 𝐾 Inertia 𝐽 Silhouette 𝑠 Structural Interpretation 𝐾=2

165,530

𝐾=3

114,930

𝐾=4

88,253

𝐾=5

68,978

𝐾=6

58,310

𝐾=7

49,814

𝐾=8

44,585

0.4910 Optimal partition (Warm/Dry vs. Cool/Humid) 0.4040 Fragmented intermediate transition zone 0.3760 Arbitrary seasonal sub-division 0.3910 Sub-cluster fragmentation 0.3740 Over-segmented clusters 0.3780 Over-segmented clusters 0.3480 Over-segmented clusters

E. Performance of the Cross-Validation MLP Evaluated on an 80/20 train/test split of the 𝑛 = 205,005 dataset, the MLP achieved: • Root Mean Squared Error (RMSE): 10.54% Relative Humidity • Coefficient of Determination (𝑅 2 ): 0.672 This demonstrates that 67.2% of relative humidity variance is predictable directly from ambient temperature and cyclical temporal encodings. RMSE = 10.54% (humidity), R2 = 0.672. Significantly higher than the R2 = 0.019 obtained on the synthetic validation set (a deliberately noise-corrupted negative control, see Section IV-F), this result confirms that a substantial portion (67%) of humidity variance is indeed predictable solely from temperature and the temporal cycle, consistent with the thermohygrometric coupling highlighted by the PCA (Section V-B). The residual distribution (Figure 12) is slightly asymmetrical (median −1.68%, standard deviation 10.48%, max 54.07%); 0.62% of the observations exhibit a residual greater than 3×RMSE, forming a candidate subset for investigating slow sensor drifts not captured by the point-based detectors in Section V-C. 1) Residual Analysis & Sensor Drift Candidates: The residual distribution displays a median of −1.68% and a standard deviation of 10.48%. A subset of 1,271 observations (0.62% of the corpus) produced prediction residuals exceeding the drift threshold:

𝐾 = 2 clearly dominates (silhouette 0.491), as the main structure of the real-world data is—as shown in Figure 11a—a simple bipartition along PC1 (high T/low H vs. low T/high H), directly reflecting the strong anti-correlation in Section V-C, rather than a fine-grained segmentation into multiple seasonal |𝑦 𝑖 − 𝑦ˆ 𝑖 | > 3 × RMSE =⇒ |𝑦 𝑖 − 𝑦ˆ 𝑖 | > 31.62% (16) regimes. Optimal partitioning occurs at 𝐾 = 2 (𝑠 = 0.4910), reflecting Because these frames pass univariate and multivariate pointa clear separation along PC1 between high-temperature/lowbased filters, they represent primary candidates for long-term humidity regimes and low-temperature/high-humidity regimes. calibration drift and localized hardware degradation. 1) DBSCAN Density Partitioning & Quantization Artifacts: Applying DBSCAN with 𝜀 = 0.071 and minPts = 10 yielded: VI. Discussion • Primary Cluster: 1 main cluster containing 203,791 A. Methodological Value of Systematic Coverage Audits points (99.41% of the dataset). The empirical findings underscore the necessity of perform• Noise Observations: 339 isolated spatial points (0.17%). ing explicit coverage audits in multi-sensor LPWAN deploy• Quantization Micro-Clusters: 12 small, isolated subments. Had the pipeline assumed full variable co-occurrence clusters (875 points total). (𝑇, 𝐻, 𝑃, CO2 ), execution would have failed completely (0 Further analysis revealed that the 12 micro-clusters are retained frames). caused by coarse sensor resolution (0.1°C for temperature,

Fig. 11. Actual PCA projection (PC1–PC2, 𝑛 = 205,005) and partitioning: (a) K-means bipartition along PC1, colored by cluster label; (b) DBSCAN partition, dominated by a single continuous cluster (label −1 = noise), with the 12 quantization micro-clusters and noise points overlaid. Schematic reconstruction preserving reported silhouette and cluster statistics.

3) Density-Based Micro-Clusters: The discovery of 12 DBSCAN micro-clusters demonstrates how coarse measurement quantization (0.1 °C, 1%) can distort density-based clustering. Introducing low-amplitude random dithering prior to spatial density estimation can prevent artificial micro-cluster formation. C. Complementarity of Point-Based vs. Predictive Anomaly Detectors

Fig. 12. MLP regression residual distribution (| 𝑦𝑖 − 𝑦ˆ 𝑖 |), predictive crossvalidation (humidity estimation). A subset of 𝑛 = 1,271 frames (0.62%) exceeds the 3×RMSE drift threshold. Density curve is a schematic reconstruction matching the reported median (−1.68%) and standard deviation (10.48%).

Furthermore, naive cross-device joins without device-level auditing risk merging asynchronous streams from physically separated sensors, introducing artificial correlations into downstream models. B. Sensor Physics, Saturation Limits & Quantization Artifacts 1) Solar Irradiance Exposure: Peak recorded temperatures (46.10°C) highlight solar radiation effects on unshaded sensor enclosures, demonstrating that outdoor node deployments require radiation shields to collect accurate ambient air temperatures. 2) Humidity Saturation Truncation: Persistent saturation at 100% reflects physical limits in capacitive humidity elements under wet conditions. Future pipeline revisions should treat upper-bound saturation as a known sensor limit rather than an extreme anomaly.

No single statistical model captures all anomaly types in longterm IoT telemetry streams. Point-based methods (Modified Z-score, Isolation Forest) excel at detecting sudden spikes and illegal joint states, but remain blind to gradual calibration shifts that stay within valid static bounds. In contrast, cross-validating physical relationships using a supervised model (MLP) isolates slow calibration drifts by flagging sustained regression residuals (|𝑦 − 𝑦ˆ | > 3 × RMSE). D. Operational Recommendations for IoT Gateways 1. In-Pipeline Auditing: Embed automated coverage audits per variable and per device within data ingestion software prior to multivariate modeling. 2. Closed-Form Subspace Processing: Use closed-form 2 × 2 PCA matrix formulations for edge devices to reduce computational overhead. 3. Spatial Acceleration & Memory Caps: Enforce spatial grid indexing and explicit queue caps when running density-based algorithms on edge hardware to prevent memory allocation faults (std::bad_alloc). 4. Multi-Tiered Anomaly Detection: Combine point-based statistical filters with predictive ML models to capture both immediate spikes and gradual hardware degradation over time.

Fig. 13. Actual performance of the cross-validated MLP (humidity prediction). Panel (a): observed vs. predicted humidity. Panel (b): residual histogram. Schematic reconstruction preserving reported RMSE and 𝑅 2 .

VII. Conclusion & Outlook This work presents a C++ pipeline for long-term LoRaWAN telemetry processing, evaluated on 421,937 raw records from the Tour Perret deployment. An automated coverage audit proved that CO2 was completely absent (0%) and barometric pressure was structurally isolated, refocusing analysis on a clean bivariate core dataset of 205,005 Temperature-Humidity frames (95.76% retention). Using closed-form 2 × 2 PCA, grid-accelerated DBSCAN, and an MLP regression model, the pipeline demonstrated the complementary nature of multi-tiered anomaly detection: • Univariate Modified Z-scores flagged 0% extreme spikes. • Multivariate Isolation Forest identified 2.12% atypical joint states. • MLP cross-validation isolated 0.62% candidate frames displaying long-term calibration drift. A. Future Research Directions 1) Hardware verification: coordinate with deployment maintainers to confirm whether CO2 sensors were omitted or require updated firmware decoders. 2) Temporal stream alignment: realign the isolated Wyres pressure stream with Elsys EMS readings using hourly window aggregation to enable trivariate (𝑇, 𝐻, 𝑃) modeling. 3) Sensor dithering and saturation handling: add controlled dither noise to prevent quantization micro-clusters in DBSCAN, and implement explicit truncation handling for humidity saturation at 100%. 4) Node-specific drift tracking: evaluate regression residuals individually per devEUI to isolate hardware degradation on specific sensor nodes.

VIII. Data & Code Availability A. Open-Access Dataset • Citation: Donsez, D., Alphand, O., & Albarel, N. (2023).

The Tour Perret LoRaWAN frames dataset [Data set]. PerSCiDO, Université Grenoble Alpes. • DOI: 10.18709/perscido.2023.06.ds395 B. Source Code & Dependencies • Repository Structure:

– 01_extract.cpp — Data ingestion, coverage audit, deduplication, and bounds filtering. – 02_anomaly_detection.cpp — Univariate Modified Z-score and bivariate Isolation Forest. – 03_pca.cpp — Closed-form 2 × 2 PCA and reconstruction error computation. – 04_clustering.cpp — K-means++ and grid-accelerated DBSCAN clustering. – 05_mlp_regression.cpp — Cyclical feature encoding and MLP regression for drift detection. • Dependencies: Standalone C++17 standard library implementation. Header-only nlohmann/json library used exclusively for raw JSON parsing. C. Processing Pipeline and Source Code Architecture Compilation command for the standalone C++ pipeline: g++ -O2 -std=c++17 0X_file.cpp -o executable 1) Module Descriptions: 1. 01_extract.cpp: Parses JSON-Lines payloads, executes coverage audit, deduplicates frames via devEui|fCount, filters physical bounds (𝑇 ∈

[−40, 85]°C, 𝐻 ∈ [0, 100]%), and exports cleaned solver such as Jacobi rotation [9] (cited in Section II-B as bivariate frames. the general-𝑑 alternative)—a self-contained implementation 2. 02_anomaly_detection.cpp: Evaluates Modified without external linear algebra library dependencies (suitable Z-scores (|𝑀𝑖 | > 3.5) and builds an Isolation Forest for embedded environments / IoT gateways). Projection onto ensemble (𝑡 = 100, 𝜓 = 256, 𝑠 > 0.62). the primary component: 3. 03_pca.cpp: Standardizes features, solves closed-form 2 × 2 covariance matrix eigenvalues/eigenvectors, and 𝑌PC1 = 𝑋 ′ v1 , calculates PC2 reconstruction errors (𝑒 𝑖 = |𝑦 𝑖,PC2 |). ′ reconstruction error 𝑒 𝑖 = ∥x𝑖′ − v1 v⊤ 4. 04_clustering.cpp: Executes K-means++ (𝐾 ∈ 1 x𝑖 ∥ 2 = 𝑦 𝑖,PC2 {2..8}) with silhouette optimization, and grid-accelerated The reconstruction error constitutes a 3rd anomaly score, DBSCAN (𝜀 = 0.071, minPts = 10). multivariate and based on the global covariance structure 5. 05_mlp_regression.cpp: Encodes cyclical tempo- (complementary to the point-based scores from 02). ral features (sin/cos), trains a 5-16-1 MLP using back5) 04_clustering.cpp — Clustering: propagation with momentum and gradient clipping, and • K-means++ [11] + Lloyd [10], with selection of 𝐾 via flags drift candidates where |𝑦 𝑖 − 𝑦ˆ 𝑖 | > 3 × RMSE. silhouette index [12] over 𝐾 = 2..8: All programs are self-contained (C++17 standard library + nlohmann::json for JSON parsing only), operate sequen𝐾 ∑︁ ∑︁ 𝑏(𝑖) − 𝑎(𝑖) tially as follows: 𝐽= ∥𝑥 − 𝜇 𝑘 ∥ 2 , 𝑠(𝑖) = max(𝑎(𝑖), 𝑏(𝑖)) The processing sequence, shown in Figure 14, mirrors the 𝑘=1 𝑥 ∈𝐶 𝑘

pipeline diagram of Figure 3.

• DBSCAN [13]: natively identifies noise (label −1), pro-

2) 01_extract.cpp — Extraction and Cleaning: Parses viding a 4th anomaly detector based on local density the JSON-Lines file via a recursive, alias-tolerant search and without assuming spherical shapes—orthogonal to across each frame’s payload hierarchy, executes the perIsolation Forest. variable/per-device coverage audit (Section IV-A, Table II), re6) 05_mlp_regression.cpp — Cross-Validation via tains Elsys EMS frames with co-occurring temperature and humidity fields (the only pair jointly available at sufficient Supervised ML: Multilayer Perceptron [14], 1 hidden layer (16 scale — CO2 is 0% present and pressure is structurally isolated neurons, tanh), trained via backpropagation + SGD-momentum on the separate Wyres stream, see Section I-B), converts ISO (with gradient clipping, [15]) to predict Relative Humidity 8601 timestamps to Unix epoch, eliminates multi-gateway (𝐻) from Temperature (𝑇) and cyclically encoded hour and reception duplicates (using key devEui|fCount), and dis- day of the year (sin, cos) — the target was revised from the cards values outside the physical boundaries of the sensor originally planned CO2 prediction because the coverage audit (Section IV-A) found CO2 entirely absent from the corpus (see (manufacturer datasheet: 𝑇 ∈ [−40, 85] °C, 𝐻 ∈ [0, 100]%). 3) 02_anomaly_detection.cpp — Anomaly Detec- Section IV-E): tion: Two complementary detectors: 1 ∑︁ (a) Modified Z-score (MAD) — [4]. Robust to extreme 𝐿= (𝑦 𝑖 − 𝑦ˆ 𝑖 ) 2 , 𝑛 𝑖 values, unlike the classic Z-score based on mean and standard deviation: 𝑣 𝑡+1 = 𝜇𝑣 𝑡 − 𝜂∇𝐿, 𝑊𝑡+1 = 𝑊𝑡 + 𝑣 𝑡+1 𝑥𝑖 − 𝑥˜ , MAD = median𝑖 (|𝑥 𝑖 − 𝑥|) ˜ MAD Recommended threshold: |𝑀𝑖 | > 3.5. (b) Isolation Forest — [5]. Constructs 𝑡 isolation trees via random splits; an anomalous point isolates with a shorter average path length: 𝑀𝑖 = 0.6745 ·

A high residual |𝑦 𝑖 − 𝑦ˆ 𝑖 | on an observation unflagged by 02/04 serves as an indicator of slow sensor drift (gradual bias), undetectable by point-based methods—a strong novelty argument for the Discussion. IX. Complete glossary of all technical concepts, algorithms, metrics and acronyms A. Acronyms & IoT Infrastructure

2(𝑛 − 1) 𝑠(𝑥, 𝑛) = 2 , 𝑐(𝑛) = 2𝐻 (𝑛 − 1) − , 𝑛 𝐻 (𝑘) = ln(𝑘) + 𝛾 (𝛾 ≈ 0.5772) ( 𝑥) ] − 𝐸 𝑐[ℎ(𝑛)

𝑠 → 1: strong anomaly; 𝑠 ≈ 0.5: normal point. Selected empirical threshold: 𝑠 > 0.62. 4) 03_pca.cpp — Dimensionality Reduction (PCA): Standardization 𝑥𝑖′ 𝑗 = (𝑥𝑖 𝑗 − 𝑥¯ 𝑗 )/𝜎 𝑗 , covariance matrix 𝑆 = 1 ′⊤ ′ 𝑛−1 𝑋 𝑋 , diagonalized using the closed-form 2 × 2 algebraic solution (Eqs. (8)–(10), Section IV-C) rather than an iterative

• ChirpStack / The Things Network (TTN): LoRaWAN

Network Servers responsible for data ingestion, raw frame decoding, and managing connectivity for IoT devices. • devEUI (Device Extended Unique Identifier): A globally unique 64-bit hardware identifier assigned to each LoRaWAN node. • fCount (Frame Counter): A counter incremented by the sensor node for each transmission. Paired with devEUI, it forms a unique deduplication key to eliminate multigateway duplicate receipts.

tourperret.log(.gz)

01_extract.cpp (JSON parsing, dedup, physical bounds filter)

clean_data.csv

02_anomaly_detection.cpp

05_mlp_regression.cpp 03_pca.cpp mlp_report.csv (residuals = 4th anomaly signal)

anomaly_report.csv pca_scores.csv (PC1, PC2, reconstruction error)

04_clustering.cpp

clusters.csv (K-means + DBSCAN labels)

Fig. 14. Sequential execution flow of the five C++ modules and their intermediate CSV reports.

• LoRaWAN (Long Range Wide Area Network): A low-

power, long-range Media Access Control (MAC) protocol designed for wireless, low-bitrate communication among IoT devices. • LPWAN (Low-Power Wide-Area Network): A category of wireless telecommunication networks optimized for long-range communication and multi-year battery operation. B. Data Engineering & Stream Quality • Coverage

Audit: An automated inspection procedure that evaluates payload completeness and variable cooccurrence per device (devEUI) and application before running multivariate analyses. • Multi-Gateway Deduplication: Filtering out duplicate receipts generated when a single LoRaWAN transmission frame is captured by multiple surrounding gateways, identified via the (devEui, fCount) tuple. • Sensor Drift: A progressive degradation in sensor calibration (e.g., saturation or aging of a capacitive humidity element) causing a systematic bias over time without breaching static physical bounds. • Quantization Micro-Clusters: Artificial high-density clusters formed during DBSCAN execution due to coarse

sensor measurement steps (e.g., 0.1°C or 1%), resulting in large stacks of identical point values. • Sensor Saturation: A physical state where the measured environmental parameter reaches or exceeds the sensor’s maximum readable limit (e.g., 100% relative humidity during heavy rainfall). C. Algorithms & Anomaly Detection • DBSCAN (Density-Based Spatial Clustering of Appli-

cations with Noise): A spatial clustering algorithm that groups points located in dense neighborhoods (𝜀, minPts) and flags sparse boundary points as noise (label −1). • Isolation Forest: An unsupervised anomaly detection algorithm that isolates multivariate outliers using random decision trees. Anomalies require significantly shorter tree path lengths (𝐸 [ℎ(𝑥)] ≪ 𝑐(𝑛)). • K-means++: An optimized variant of K-means clustering that uses probabilistic centroid initialization based on squared distances, avoiding convergence to poor local minima. • MLP (Multilayer Perceptron): A feedforward artificial neural network trained via backpropagation. In this pipeline, it performs predictive cross-validation to estimate humidity from temperature and cyclical time encodings.

• PCA (Principal Component Analysis): A linear dimen-

sionality reduction technique that identifies orthogonal directions of maximum variance. Solved via an exact closed-form matrix formulation in the 2 × 2 bivariate case. • Modified Z-Score: A robust univariate outlier metric computed using the sample median and Median Absolute Deviation (MAD), making it resilient to extreme outliers compared to standard mean-based Z-scores. D. Mathematical & Algorithmic Concepts • Thermo-Hygrometric

Opposition Axis (PC1): The primary principal component capturing 83.10% of total variance along the inverse correlation axis between temperature and relative humidity (𝐻 ′ ≈ −𝑇 ′ ). • Thermo-Hygrometric Energy Axis (PC2): The secondary principal component orthogonal to PC1 (16.90% variance). Magnitude along this axis measures departures from standard atmospheric thermodynamic coupling. • Cyclical Temporal Encoding: A trigonometric transformation (sin, cos) mapping periodic time dimensions (hours of the day, days of the year) onto continuous 2D coordinates without boundary discontinuities (23h → 00h). • Reconstruction Error (𝑒 𝑖 ): The Euclidean distance between an original standardized data vector and its projection onto the primary subspace. In 2D, it equals the absolute coordinate value along the second principal component: 𝑒 𝑖 = |𝑦 𝑖,PC2 |. • 𝑘-Distance Plot: A calibration heuristic that plots sorted distances to each point’s 𝑘th-nearest neighbor (minPts) to identify the maximum curvature point (“knee”) for setting the optimal DBSCAN search radius 𝜀. • Uniform Grid Indexing: A spatial partitioning optimization that divides feature space into square cells of width 𝜀, restricting neighbor queries to the point’s cell and its 8 adjacent neighbors to scale DBSCAN execution time from 𝑂 (𝑛2 ) down to 𝑂 (𝑛). • MAD (Median Absolute Deviation): A robust measure of statistical dispersion defined as MAD = median𝑖 (|𝑥 𝑖 − 𝑥|). ˜ • Silhouette Index: A metric assessing clustering quality by comparing intra-cluster distance against nearest-cluster distance (𝑠 ∈ [−1, 1]). E. Regression Evaluation Metrics 𝑅 2 (Coefficient of Determination): The proportion of variance in the target variable (Relative Humidity) explained by the MLP model (𝑅 2 = 0.672). • RMSE (Root Mean Squared Error): The square root of the average squared residuals, measuring prediction error (RMSE = 10.54%). Observations yielding residuals |𝑦 𝑖 − 𝑦ˆ 𝑖 | > 3 × RMSE are flagged as primary sensor drift candidates. •

References [1] D. Donsez, O. Alphand, and N. Albarel, “The tour perret lorawan frames dataset,” [Data set], 2023. [2] V. Chandola, A. Banerjee, and V. Kumar, “Anomaly detection: A survey,” ACM Computing Surveys, vol. 41, no. 3, pp. 1–58, 2009. [3] C. C. Aggarwal, Outlier Analysis, 2nd ed. Springer, 2017. [4] B. Iglewicz and D. C. Hoaglin, How to Detect and Handle Outliers. ASQC Quality Press, 1993. [5] F. T. Liu, K. M. Ting, and Z.-H. Zhou, “Isolation forest,” in Proceedings of the IEEE International Conference on Data Mining (ICDM), 2008, pp. 413–422. [6] S. Bhatia et al., “Real-time anomaly detection in edge streams,” ACM Transactions on Knowledge Discovery from Data (TKDD), vol. 16, no. 4, pp. 1–25, 2022. [7] K. Pearson, “On lines and planes of closest fit to systems of points in space,” Philosophical Magazine, vol. 2, no. 11, pp. 559–572, 1901. [8] H. Hotelling, “Analysis of a complex of statistical variables into principal components,” Journal of Educational Psychology, vol. 24, no. 6, pp. 417–441, 1933. [9] G. H. Golub and C. F. Van Loan, Matrix Computations, 4th ed. Johns Hopkins University Press, 2013. [10] S. Lloyd, “Least squares quantization in pcm,” IEEE Transactions on Information Theory, vol. 28, no. 2, pp. 129–137, 1982. [11] D. Arthur and S. Vassilvitskii, “k-means++: The advantages of careful seeding,” in Proceedings of the Eighteenth Annual ACM-SIAM Symposium on Discrete Algorithms (SODA), 2007, pp. 1027–1035. [12] P. J. Rousseeuw, “Silhouettes: A graphical aid to the interpretation and validation of cluster analysis,” Journal of Computational and Applied Mathematics, vol. 20, pp. 53–65, 1987. [13] M. Ester, H.-P. Kriegel, J. Sander, and X. Xu, “A density-based algorithm for discovering clusters in large spatial databases with noise,” in Proceedings of the Second International Conference on Knowledge Discovery and Data Mining (KDD), 1996, pp. 226–231. [14] D. E. Rumelhart, G. E. Hinton, and R. J. Williams, “Learning representations by back-propagating errors,” Nature, vol. 323, no. 6088, pp. 533–536, 1986. [15] R. Pascanu, T. Mikolov, and Y. Bengio, “On the difficulty of training recurrent neural networks,” in International Conference on Machine Learning (ICML), 2013, pp. 1310–1318. [16] P. Ferrer-Cid, J. M. Barcelo-Ordinas, and J. Garcia-Vidal, “A review of graph-powered data quality applications for iot monitoring sensor networks,” Journal of Network and Computer Applications, 2025, preprint: arXiv:2410.21006, 2024. [17] Elsys AB, Elsys EMS — Environmental Monitoring Sensor, Technical Specification, Elsys AB.

Related documents

• Data Quality Audit and Bivariate Temperature-Humidity Analysis of a Long-Term LoRaWAN Network: Anomaly Detection, PCA, Clustering and Machine Learning#1447805
• Efficient detection of intrusions in TON-IoT dataset using hybrid feature selection approach.#1232
• Estimação de dados faltantes em séries temporais de irradiação solar: uma abordagem baseada em aprendizado de máquina#1148548
• Adaptive edge-cloud isolation forest-autoencoder pipeline in real-time with delta policy feedback for zero-day IoT smartwatch intrusion detection.#1374151
• Regularities in creep and hot-fatigue data part 1#491731