SUBMITTED TO IEEE INTERNET OF THINGS JOURNAL, SEPTEMBER 2026
1
Radio Map Construction with Post-Hoc Location Calibration under Quasi-Static Positioning Errors: Joint Estimation, Performance Bounds, and GNSS-Based Evaluation
arXiv:2609.11142v1 [eess.SP] 10 Sep 2026
Koki Kanzaki, Graduate Student Member, IEEE, Katsuya Suto, Senior Member, IEEE and Koya Sato, Senior Member, IEEE
Abstract—Radio maps enable environment-aware wireless and Internet-of-Things applications and can be constructed from location-tagged received signal strength (RSS) measurements collected by mobile devices. In urban environments, temporally correlated GNSS errors can shift an entire sensing trajectory, causing systematic spatial misregistration that is not mitigated by collecting more measurements. This paper presents a radiomap construction framework that uses the radio measurements themselves to calibrate erroneous location tags after data collection. The dominant positioning error is modeled as a sensor-specific quasi-static offset, which is jointly estimated with radio-propagation parameters in a Gaussian process regression (GPR) framework by exploiting complementary spatial information from distance-dependent path loss and spatially correlated shadowing. We establish lower and upper bounds on the conditional Bayes risk and show that, under a translationinvariant trajectory model, trajectory information alone cannot identify the quasi-static offset, thereby motivating the use of RSS-derived spatial information for calibration. Numerical evaluations across propagation conditions show that the proposed method reduces the mean squared error (MSE) gap from ideal GPR to approximately 3.26 dB2 , compared with about 10 dB2 for position-error-agnostic and noisy-input GPR baselines. Evaluation using positioning-error models derived from smartphone GNSS measurements shows that the proposed method outperforms a KF–RTS trajectory-smoothing baseline despite unmodeled time-varying positioning errors, remaining within approximately 5 dB2 of ideal GPR at the median MSE. These results demonstrate that RSS measurements can serve not only as observations for radio-map reconstruction but also as spatial cues for post-hoc calibration of imperfectly geotagged sensing data. Index Terms—Gaussian process, radio map construction, quasi-static error, bound analysis
I. I NTRODUCTION
E
NVIRONMENT-aware wireless communication requires site-specific knowledge of radio-propagation conditions. Radio maps provide such spatial knowledge by associating A part of this work was presented at the IEEE GLOBECOM Workshops [1] (DOI: 10.1109/GCWkshps68340.2025.11591047). This work was supported by JST, PRESTO under Grant Number JPMJPR23P3, CRONOS under Grant Number JPMJCS24N1, ASPIRE under Grant JPMJAP2346, BOOST under Grant Number JPMJBS2415, and JSPS KAKENHI under Grant Number 25K00138. Corresponding author: Koki Kanzaki. Koki Kanzaki and Koya Sato are with the Artificial Intelligence eXploration Research Center, The University of Electro-Communications, Tokyo 1828585, Japan (e-mail: [email protected] and k [email protected]). Katsuya Suto is with the Faculty of Information Science and Technology, Hokkaido University, Sapporo, Hokkaido 060-0808, Japan (e-mail: [email protected]).
radio-propagation quantities, such as received signal strength (RSS), path loss, and channel gain, with geographical locations [2]. They have been widely investigated as a basis for wireless resource management, localization, network optimization, and other environment-aware applications [3]–[5]. More broadly, spatial representations of radio-environment knowledge are also important components of network digital twins and channel knowledge maps (CKMs) for next-generation wireless systems [6]–[8]. Measurement-based radio-map construction is particularly attractive because it can directly capture site-specific propagation characteristics. In a typical mobile-sensing architecture, terminals such as smartphones collect radio measurements while moving, associate each measurement with an estimated location, and upload the resulting data to a remote server. The server aggregates measurements collected by multiple terminals and estimates the radio conditions at unobserved locations through spatial interpolation. A variety of modeldriven and data-driven estimators have been developed for this purpose [9]–[11]. Among them, Gaussian process regression (GPR) [12] provides a probabilistic interpolation framework that can explicitly exploit the spatial correlation of shadowing to reconstruct radio maps from sparse measurements. A fundamental difficulty in such measurement-based mapping is that the measurement coordinates themselves can be inaccurate. This is particularly relevant to outdoor mobile sensing, where measurement locations are commonly obtained using a Global Navigation Satellite System (GNSS). In urban environments, surrounding buildings, non-line-of-sight (NLOS) reception, multipath propagation, and satellite geometry can substantially degrade GNSS positioning accuracy [13]. Moreover, GNSS positioning errors are not necessarily independent across time. When the propagation conditions of GNSS signals or the geometry of the received satellites remain similar over a certain period, the resulting positioning displacement can persist across consecutive measurements. Indeed, positioning errors observed along mobile trajectories have been reported to exhibit strong temporal autocorrelation [14], [15]. Consequently, a sequence of measurements collected by the same terminal can be coherently displaced from the true trajectory rather than independently scattered around it. In this paper, we focus on the dominant quasi-static component of such temporally correlated positioning errors, which we refer to as a quasi-static positioning bias. Specifically,
10 True Position 9 Observed Position 8 7 6 5 4 3 2 1 00 1 2 3 4 5 6 7 8 9 10
𝑥 [m]
1 (a)
𝑦 [m]
𝑦 [m]
SUBMITTED TO IEEE INTERNET OF THINGS JOURNAL, SEPTEMBER 2026
10 True Position 9 Observed Position 8 7 6 5 4 3 2 1 00 1 2 3 4 5 6 7 8 9 10
𝑥 [m]
1 (b)
Fig. 1: Illustrative examples of positioning errors: (a) i.i.d. errors and (b) trajectory-wise quasi-static bias.
we model this component as a position offset shared by the measurements along each trajectory. In the mobile-sensing setting considered in this paper, each sensor generates one measurement trajectory and is therefore associated with one shared positioning offset. Fig. 1 highlights an important distinction between pointwise positioning uncertainty and the trajectory-wise bias. In Fig. 1(a), independent positioning errors perturb individual measurements around the true trajectory. In contrast, in Fig. 1(b), a sequence of reported positions is coherently displaced in a common direction. Such a shared offset is not averaged out by collecting additional measurements along the same trajectory. Instead, multiple radio measurements remain systematically registered at incorrect spatial locations, which can directly produce spatial misregistration in the reconstructed radio map. Two classes of existing approaches are particularly relevant to this problem. The first accounts for uncertainty in measurement locations within the regression model. For example, noisy-input Gaussian processes (NIGPs) approximate the effect of Gaussian input noise as additional output uncertainty [16], while more general GP-based methods explicitly incorporate uncertain or latent inputs into the inference procedure [17], [18]. For radio-map construction, positional uncertainty has also been incorporated into GP inference through variational methods [19]. These approaches provide principled mechanisms for handling uncertainty in individual input locations. However, the trajectory-wise quasi-static bias is qualitatively different: a common offset coherently translates a sequence of measurements along the same trajectory. Thus, treating such a shared displacement only as pointwise input uncertainty does not remove the resulting trajectory-level spatial misregistration. The second class of approaches attempts to improve the positioning trajectory before radio-map construction. Trajectoryestimation methods based on the Kalman filter (KF), the Rauch–Tung–Striebel (RTS) smoother, and GNSS/IMU fusion can exploit motion models and temporal error correlation to suppress time-varying positioning errors [20]–[22]. However, trajectory information alone does not necessarily determine a trajectory-wise constant positioning offset. Without sufficiently informative absolute-location constraints, a constant translation of the underlying trajectory can be compensated for by
2
an opposite change in the positioning bias while producing the same reported trajectory. Thus, temporal smoothing can reduce time-varying positioning errors while leaving a residual absolute spatial misregistration that directly affects subsequent radio-map construction. The key observation of this work is that the radio measurements themselves provide spatial information that is unavailable to trajectory-only correction. The mean RSS depends on the transmitter–receiver distance through distance-dependent path loss and therefore provides information about the absolute placement of a measurement trajectory relative to a known base station (BS). At the same time, spatially correlated shadowing provides pairwise information about the relative placement of measurements collected by different sensors. Hence, RSS observations can provide additional spatial constraints for estimating trajectory-wise positioning offsets. From this perspective, radio-map construction and location calibration are coupled inference problems: the measurements used to reconstruct the radio environment can also help determine where those measurements were actually collected. Motivated by this observation, we propose a GPR-based radio-map construction framework with post-hoc location calibration under quasi-static positioning biases. The proposed method models the dominant positioning error as an unknown offset shared by the measurements along each trajectory and jointly estimates these offsets and the radio-propagation parameters from the RSS observations. The estimation exploits both the distance-dependent mean and the spatial covariance structure of RSS. The estimated offsets are then used to calibrate the measurement coordinates, and the calibrated locations are used as the GPR inputs for radio-map construction. In this way, the same RSS measurements serve both as observations for reconstructing the radio map and as spatial information for post-hoc calibration of their associated location tags. The main contributions of this work are summarized as follows. • We develop a GPR-based radio-map construction framework that jointly performs post-hoc location calibration and radio-map estimation under trajectory-wise quasi-static positioning biases. By exploiting distancedependent path loss and spatially correlated shadowing, the proposed method jointly estimates the positioning offsets and radio-propagation parameters and reconstructs the radio map using the calibrated measurement locations. • We characterize the fundamental role of radio observations in resolving trajectory-wise positioning offsets. We show that, under a translation-invariant trajectory law, the reported positioning trajectory alone does not update the prior distribution of a trajectory-wise constant offset, clarifying why additional spatial information is required for post-hoc calibration. We further establish a bounding framework based on a conditional minimummean-square-error (MMSE) benchmark, deriving lower and upper bounds on the conditional Bayes risk and thereby bounding the excess risk of an arbitrary radiomap estimator without explicitly evaluating the Bayesoptimal estimator.
SUBMITTED TO IEEE INTERNET OF THINGS JOURNAL, SEPTEMBER 2026
•
We evaluate the proposed method over a broad range of radio-propagation conditions and further assess its robustness to positioning-error model mismatch. Using positioning-error models derived from measured smartphone GNSS data, we consider temporally correlated error components that are not explicitly represented by the proposed quasi-static offset model. In addition to position-error-agnostic and uncertain-input methods, we compare the proposed method with a KF–RTS-based trajectory-correction baseline that explicitly exploits temporal error correlation, thereby assessing the benefit of using the spatial structure of the radio observations themselves for location calibration.
Sec. II describes the system model. Secs. III and IV present the baseline and proposed radio-map construction methods, respectively. Sec. V presents the bounding framework for evaluating the deviation from Bayes-optimal performance. Sec. VI evaluates the proposed method under various radiopropagation conditions. Sec. VII further evaluates its robustness to temporally correlated positioning-error components that are not explicitly represented by the quasi-static offset model, using error models derived from measured GNSS data. Finally, Sec. VIII concludes this paper. Notations: The transpose and inverse are denoted by (·)⊤ and (·)−1 , respectively, and |A| denotes the determinant of a square matrix A. The Euclidean norm is denoted by ∥·∥, and tr(·) denotes the trace. E[·] and Cov[·] denote expectation and covariance, respectively. In Sec. V, Eθ [·] and Covθ (·) denote expectation and covariance under the generative model with θ fixed at its true value. II. S YSTEM M ODEL Fig. 2 provides an overview of the system model. We consider the construction of a radio map that estimates the received power from a single base station (BS) at each location. In the target area, N mobile devices measure the received power while moving and report the measurements to a remote server together with their corresponding observed coordinates. After collecting sufficient measurements, the server constructs a radio map by spatially interpolating the RSS at unobserved locations. The radio propagation parameters are also estimated during this process. Let the target area be a two-dimensional region A ⊂ R2 , and let xTx ∈ A denote the coordinates of the BS, which are assumed to be known. Assuming that the effects of small-scale fading are largely suppressed through signal averaging during the measurement process, the received power PRx [dBm] at x ∈ A can be modeled as [23], ∥xTx − x∥ + W (x) (1) PRx (x) = PTx − 10η log10 d0 ≜ f (x), (2) where PTx [dBm] denotes the transmit power of the BS, η denotes the path-loss index, ∥·∥ denotes the Euclidean distance, and W (·) denotes the shadowing term. The constant d0 is a reference distance. Following Gudmundson’s shadowing
3
Client-side processing 1. Observe RSS and position 2. Upload to the server Observed trajectory
Observed RSS
Quasi-static error
Upload to the Server Remote Server Processing 1. Collect data 2. Correct positions and perform GPR Corrected trajectory
Constructed radio map
Error correction
Fig. 2: Overview of the system model.
correlation model [24], the shadowing component can be modeled as W (·) ∼ GP (0, kexp (·, ·)) , (3) where kexp (·, ·) is the exponential kernel defined as ∥x − x′ ∥ ln 2 kexp (x, x′ ) = σf2 exp − . dcor
(4)
Here, σf2 is the variance of the shadowing term, and dcor is the correlation distance. Let M (i) denote the number of measurements collected by (i) the i-th sensor, and let xj denote its location at time index (i) j, for j = 1, · · · , M . The corresponding RSS measurement, (i) denoted by pj , is modeled as (i) (i) (i) (i) pj = f xj + ϵp,j , ϵp,j ∼ N 0, σp2 , (5) (i)
where ϵp,j denotes Gaussian measurement noise that captures thermal noise and residual small-scale fading remaining after signal averaging and σp2 is its variance. We next model the errors in the reported coordinates. As discussed in Sec. I, GNSS positioning errors can exhibit strong temporal persistence. As a tractable approximation, we model their dominant quasi-static component as a sensorspecific constant offset over the observation interval. Residual time-varying positioning errors that are not captured by this (i) model are considered in Sec. VII. Let x̃j denote the position reported by sensor i at time index j. Using the true position (i) xj , the observed position is modeled as (i)
(i)
x̃j = xj + e(i) .
(6)
SUBMITTED TO IEEE INTERNET OF THINGS JOURNAL, SEPTEMBER 2026
Under this approximation, e(i) is assumed to remain constant over the observation interval of the i-th sensor and to be independent across sensors. Further, it follows e(i) ∼ N (µs , Σs ),
(7)
where µs ∈ R2 and Σs ∈ R2×2 denote the mean vector and the covariance matrix of the positioning error, respectively1 . We further assume that {e(i) }N i=1 is independent of the shadowing field, of the measurement noise, and of the true measurement locations. Based on the above observation model, we define the measurement data D(i) collected by the i-th sensor as o n (i) (i) (8) D(i) ≜ j = 1, · · · , M (i) . x̃j , pj After all sensors upload their local dataset, the server aggregates the local datasets as oN n . (9) D ≜ D(i)
4
the n-th measurement, let in denote the sensor index and jn denote the time index within the measurement sequence of sensor in . Then, the reported location x̃n and the RSS measurement pn are defined as x̃n ≜ x̃jnn ,
(i )
(14)
(i ) pn ≜ pjnn .
(15)
The reported locations and RSS measurements are vectorized as ⊤
X̃ ≜ [x̃1 , · · · , x̃MN ] ∈ RMN ×2 , ⊤
MN
p ≜ [p1 , · · · , pMN ] ∈ R
,
(16) (17)
respectively. The observation vector conditioned on the reported locations and model parameters follows p | X̃, θ ∼ N mθ,X̃ , Cθ,X̃ , (18)
i=1
The server is assumed to know the transmitter coordinates xTx , the transmit power PTx , and the prior statistics µs and Σs of the positioning errors. Let X∗ ⊆ A denote the set of query points at which the radio map is evaluated, and let x∗ ∈ X∗ denote a generic query point. Given these quantities and the dataset D, the server estimates the RSS at x∗ . III. R ADIO M AP C ONSTRUCTION WITHOUT ACCOUNTING FOR P OSITION E RRORS We first establish a position-error-agnostic GPR baseline that treats the reported measurement locations as the true GP inputs. This approach treats the reported measurement (i) (i) locations as the true locations; i.e., xj ≈ x̃j . Ignoring positioning error, the unknown model parameters in the system model are collected into ⊤ θ ≜ η, σf2 , dcor , σp2 . (10) Because the shadowing component is modeled as a zero-mean GP, the latent field of RSS can be expressed as f (x) ∼ GP (mθ (x), kexp,θ (x, x′ )) , where the mean function is given by ∥xTx − x∥ mθ (x) ≜ PTx − 10η log10 1 + , d0
(11)
N X
M (i) .
⊤
mθ,X̃ ≜ [mθ (x̃1 ), · · · , mθ (x̃MN )] ,
(19)
and the covariance matrix is Cθ,X̃ ≜ Kθ,X̃ + σp2 I.
(20)
Here, the kernel matrix Kθ,X̃ has entries h i Kθ,X̃ = kexp,θ (x̃m , x̃n ) , m, n = 1, · · · , MN . (21) m,n
The resulting log marginal likelihood is MN 1 log(2π) log L p | X̃, θ = − log Cθ,X̃ − 2 2 ⊤ 1 − p − mθ,X̃ C−1 p − m θ,X̃ . θ,X̃ 2
(22)
The model parameters are estimated by maximizing the log marginal likelihood: θ̂ = arg max log L p | X̃, θ , (23) θ∈Θ
(12)
and the exponential kernel is given by Eq. (4) with the corresponding parameters specified by θ 2 . Let us define the total number of measurements collected from all sensors as MN ≜
where the mean vector is
where Θ denotes the feasible parameter set. We can maximize this function based on a gradient-based method; this paper solves this problem using Adam [25]. Using the estimated parameters, the cross-covariance vector between the query point x∗ and the reported measurement locations is given by h i kθ̂,X̃ (x∗ ) = kexp,θ̂ (x∗ , x̃n ) , n = 1, · · · , MN (24) n
(13)
i=1
To express the GP likelihood in a vector form, we enumerate all measurements using a global index n = 1, · · · , MN . For 1 The prior covariance can be obtained from positioning-quality indicators or estimated empirically from positioning-error statistics. 2 The additive constant 1 is introduced for analytical convenience. Its effect is negligible over the communication distance considered in this paper, for which ∥xTx − x∥/d0 ≫ 1.
The RSS at x∗ is then estimated by the posterior mean of the latent field of RSS: h i p̂(x∗ ) ≜ E f (x∗ ) | p, X̃, θ̂ (25) = mθ̂ (x∗ ) + k⊤ (x )C−1 p − mθ̂,X̃ . (26) θ̂,X̃ ∗ θ̂,X̃ Evaluating this posterior mean over the query points yields the radio map. The complete procedure is summarized in Alg. 1.
SUBMITTED TO IEEE INTERNET OF THINGS JOURNAL, SEPTEMBER 2026
Algorithm 1 Radio map construction without accounting for positioning errors Dataset D, transmitter location xTx , transmit power PTx , and query-point set X∗ 2: Output: Radio map {p̂(x∗ )}x∗ ∈X∗ and estimated model parameters θ̂ 3: Form the reported-location matrix X̃ and RSS observation vector p from D 4: Estimate θ̂ by solving Eq. (23) 5: Form mθ̂,X̃ using Eq. (19) 6: Form Kθ̂,X̃ using Eq. (21) and Cθ̂,X̃ using Eq. (20) 7: for all x∗ ∈ X∗ do 8: Form kθ̂,X̃ (x∗ ) using Eq. (24) 9: Compute p̂(x∗ ) using Eq. (26) 10: end for 11: return {p̂(x∗ )}x∗ ∈X∗ and θ̂ 1: Input:
IV. R ADIO M AP C ONSTRUCTION AND J OINT E STIMATION OF P ROPAGATION PARAMETERS AND P OSITION E RRORS Alg. 1 treats the reported measurement locations as exact GP inputs. We now relax this assumption by treating the sensorspecific positioning offsets as unknown variables jointly estimated with the propagation parameters. Because these offsets shift the GP inputs, they affect both the distance-dependent mean and the spatial covariance of the RSS observations. The proposed method exploits these two spatial structures to jointly estimate the model parameters and the sensor-specific positioning errors. A. Joint Estimation of Propagation Parameters and Position Errors From the positioning error model in Eq. (6), the true location of the j-th measurement collected by sensor i is given by (i)
(i)
xj = x̃j − e(i) .
(27)
We collect the sensor-specific positioning errors into the vector ⊤ ⊤ ⊤ evec ≜ e(1) , · · · , e(N ) ∈ R2N . (28) Under the global measurement index introduced in the previous section, let
5
Its n-th element can be written explicitly as ∥xTx − xn (evec )∥ . mθ,X(evec ) n = PTx − 10η log10 1 + d0 (33) The corresponding covariance matrix is Cθ,X(evec ) ≜ Kθ,X(evec ) + σp2 I
(34)
where the kernel matrix Kθ,X(evec ) has entries Kθ,X(evec ) m,n = kexp,θ (xm (evec ), xn (evec )) ∥(x̃m − em ) − (x̃n − en )∥ ln(2) = σf2 exp − , (35) dcor for m, n = 1, · · · , MN . The resulting log marginal likelihood is MN 1 log(2π) log L p | X̃, θ, evec = − log Cθ,X(evec ) − 2 2 ⊤ 1 p − mθ,X(evec ) . (36) p − mθ,X(evec ) C−1 − θ,X(e ) vec 2 To clarify the spatial information provided by the mean and covariance, define rn ≜ ∥xn (evec ) − xTx ∥ ,
(37)
ρmn ≜ ∥xm (evec ) − xn (evec )∥ .
(38)
and Let δa,b denote the Kronecker delta. For rn > 0, the gradient of the n-th mean component with respect to the positioning error of sensor i is 10η xn (evec ) − xTx (39) ∇e(i) mθ,X(evec ) n = δi,in ln 10 rn (d0 + rn ) Thus, the mean provides spatial information along the radial direction from the transmitter to each corrected measurement location. For m ̸= n with ρmn > 0, the corresponding gradient of the kernel entry is ∇e(i) Kθ,X(evec ) m,n ln 2 xm (evec ) − xn (evec ) = Kθ,X(evec ) m,n (δi,im − δi,in ) dcor ρmn (40)
Hence, the covariance provides information along pairwise directions between corrected measurement locations. In parThe corrected location associated with the n-th measurement ticular, this gradient vanishes for pairs of measurements colis then lected by the same sensor, because their common positioning xn (evec ) ≜ x̃n − en . (30) offset does not change their pairwise distance. Therefore, the covariance primarily constrains the relative offsets between Collecting these locations gives different sensors, whereas the mean constrains the absolute ⊤ X(evec ) ≜ [x1 (evec ), · · · , xMN (evec )] (31) placement of the corrected measurement locations through their distances from the known transmitter. Consequently, Conditioned on the reported locations, model parameters, and when the measurement trajectories span only a limited range positioning errors, the RSS observation vector is Gaussian with of directions relative to the transmitter, or when the directions mean of inter-sensor measurement pairs are insufficiently diverse, different positioning error vectors may produce similar mean ⊤ mθ,X(evec ) ≜ [mθ (x1 (evec )) , · · · , mθ (xMN (evec ))] (32) vectors and covariance matrices. en ≜ e(in )
n = 1, · · · , MN .
(29)
SUBMITTED TO IEEE INTERNET OF THINGS JOURNAL, SEPTEMBER 2026
To regularize these weakly identifiable directions, we adopt a MAP-based formulation for the positioning errors by incorporating the prior distribution introduced in Sec. II. Because the positioning-error vectors are independent across sensors, their joint prior density is N Y pe (evec ) = ϕ e(i) ; µs , Σs ,
(41)
i=1
where 1 p 2π |Σs | 1 ⊤ × exp − (x − µs ) Σ−1 (x − µ ) . s s 2
ϕ (x; µs , Σs ) =
(42)
The prior also provides an absolute reference in the positioning-error space. Specifically, its gradient with respect to the positioning error of sensor i is (43) e(i) − µs . ∇e(i) log pe (evec ) = −Σ−1 s Thus, the prior anchors each sensor-specific positioning error around the prior mean µs and regularizes directions that are weakly constrained by the RSS likelihood. The resulting MAPbased objective is J (θ, evec ) ≜ log L p | X̃, θ, evec + log pe (evec ) (44) N log 4π 2 |Σs | = log L p | X̃, θ, evec − 2 N ⊤ 1 X (i) − e − µs Σ−1 e(i) − µs . s 2 i=1 (45) The model parameters and positioning errors are jointly estimated by solving θ̂, êvec = arg max J (θ, evec ) , (46) θ∈Θ, evec ∈E
where E denotes a finite box constraint set for the sensorspecific positioning errors. Equation (46) corresponds to MAP estimation of the positioning errors jointly with maximumlikelihood estimation of θ under the conditional model that treats the reported locations X̃ as fixed. Unlike the full generative posterior considered later in Sec. V, this conditional formulation does not include the probability law of the true measurement trajectories. Because the objective is generally nonconvex with respect to both θ and evec , the numerical optimizer is not guaranteed to attain the global maximum. In our implementation, we therefore obtain a numerical solution of Eq. (46) by minimizing the negative objective −J using the Adam optimizer [25], with the constrained parameters reparameterized so that the optimization is performed in an unconstrained space. The resulting predictor uses the estimated positioning errors and model parameters as plug-in estimates rather than marginalizing over their posterior uncertainty.
6
Algorithm 2 Radio map construction and joint estimation of propagation parameters and position errors Dataset D, transmitter location xTx , transmit power PTx , error-prior mean µs , error-prior covariance Σs , and query-point set X∗ 2: Output: Radio map {p̂(x∗ )}x∗ ∈X∗ , estimated model parameters θ̂, estimated positioning errors êvec , and corrected locations X̂ 3: Form X̃ and p from D 4: Initialize θ and evec 5: Obtain θ̂ and êvec by numerically solving Eq. (46) using Adam to minimize −J (θ, evec ) 6: Form the corrected-location matrix X̂ using Eq. (50) 7: Form mθ̂,X̂ using Eq. (32) 8: Form Kθ̂,X̂ using Eq. (35) and Cθ̂,X̂ using Eq. (51) 9: for all x∗ ∈ X∗ do 10: Form kθ̂,X̂ (x∗ ) using Eq. (52) 11: Compute p̂(x∗ ) using Eq. (53) 12: end for 13: return {p̂(x∗ )}x∗ ∈X∗ , θ̂, êvec , and X̂ 1: Input:
B. Radio Map Construction Let ê(i) denote the estimated positioning error of sensor i. The corrected estimate of its j-th measurement location is (i)
(i)
x̂j = x̃j − ê(i) .
(47)
Using the global measurement index, let ên ≜ ê(in ) ,
(48)
x̂n ≜ x̃n − ên
(49)
The corrected locations are collected into ⊤
X̂ ≜ [x̂1 , · · · , x̂MN ] ∈ RMN ×2 .
(50)
Using the estimated parameters and corrected locations, the covariance matrix is Cθ̂,X̂ = Kθ̂,X̂ + σ̂p2 I.
(51)
The cross-covariance vector between x∗ and the corrected measurement locations has entries h i kθ̂,X̂ (x∗ ) = kexp,θ̂ (x∗ , x̂n ) , n = 1, · · · , MN (52) n
The RSS at x∗ is estimated by the posterior mean −1 p̂(x∗ ) = mθ̂ (x∗ ) + k⊤ (x )C p − m θ̂,X̂ . θ̂,X̂ ∗ θ̂,X̂
(53)
Evaluating Eq. (53) over the query points yields the radio map. The complete procedure for radio map construction is summarized in Alg. 2. V. P ERFORMANCE B OUND A NALYSIS This section quantifies how closely a radio-map estimator approaches the minimum MSE achievable from the observed data. We first define the conditional Bayes risk and decompose the estimator risk into irreducible and excess components. Because the Bayes risk is generally intractable, we bound it using an oracle lower bound and a working-model upper
SUBMITTED TO IEEE INTERNET OF THINGS JOURNAL, SEPTEMBER 2026
bound, which in turn provide bounds on the excess risk. The working-model analysis also clarifies the identifiability of trajectory-wise positioning offsets from trajectory information alone.
7
positioning-error vector evec , the GP posterior mean at the evaluation points is defined as µGP (evec ; D, θ) = mθ,X∗ + K∗X(evec ) C−1 θ,X(evec ) [p − mθ,X(evec ) ].
A. Bayes Risk and Excess Risk Throughout this section, the propagation parameter vector θ is not treated as a random variable, and the estimation performance is evaluated conditional on its true, fixed value. Hereafter, Eθ [·] denotes expectation under the generative model with θ held fixed, where the expectation is taken over the true measurement trajectories, positioning errors, shadowing field, measurement noise, and, when applicable, the internal randomness of an estimator. The evaluation-point set X∗ = {x∗,1 , · · · , x∗,M∗ } is assumed to be deterministic, and the evaluation points themselves are free of positioning errors. The true radio map values at the evaluation points are defined as ⊤
f∗ ≜ [f (x∗,1 ), · · · , f (x∗,M∗ )] ∈ RM∗ .
(54)
Let the matrix containing all true measurement locations be ⊤
X = [x1 , · · · , xMN ] ∈ RMN ×2 ,
(55)
where X is assumed to be a random variable following a known trajectory-generation law that is independent of the propagation environment. Its joint probability density is denoted by pX (X). The joint prior density of the positioningerror vector evec is given by pe (evec ) ≜
N Y ϕ e(i) ; µs , Σs .
(56)
i=1
For a given measurement-location matrix X and θ, let pGP (p | X, θ) denote the multivariate Gaussian probability density function of the RSS observation vector p induced by the GP observation model. An arbitrary estimator constructed from the observed data D and an internal random variable Zδ , which is independent of both D and the true generative process, is written as f̂δ = δ(D, Zδ ).
(57)
The risk of estimator δ is defined as the MSE per evaluation point: 2 1 Eθ f∗ − f̂δ . (58) Rδ (θ) ≜ M∗ 2 For fixed reported locations X̃, the true measurement locations corresponding to a positioning-error vector evec are determined by X(evec ) defined in Sec. IV. Therefore, under the true generative model, the posterior density of the positioning errors is given by
(60)
Under squared-error loss, the Bayes-optimal estimator conditional on the fixed value of θ is f̂B ≜ Eθ [f∗ | D],
(61)
and, by the law of iterated expectations, it can be written as Z f̂B = µGP (evec ; D, θ)πtrue (evec | D, θ)devec . (62) The Bayes risk achieved by this estimator is defined as 2 1 (63) Eθ f∗ − f̂B LB (θ) ≜ M∗ 2 1 = Eθ [tr{Covθ (f∗ | D)}] . (64) M∗ Thus, LB is the minimum MSE achievable from the observed data under the assumed generative model. Because θ is fixed throughout the analysis, LB (θ) is a θ-conditional Bayes risk rather than a Bayes risk obtained by assigning a prior distribution to θ. Since πtrue contains the trajectory density pX and the positioning errors affect both the GP mean and covariance nonlinearly, f̂B and LB are generally difficult to evaluate in closed form. Because f̂B is the conditional mean of f∗ given D, the orthogonality property of conditional expectation gives h i Eθ (f∗ − f̂B )⊤ (f̂B − f̂δ ) = 0. (65) Therefore, for any estimator δ, Rδ (θ) = LB (θ) + Gδ (θ), where Gδ (θ) ≜
1 Eθ M∗
2
f̂δ − f̂B
(66)
≥0
2
(67)
is the excess risk of estimator δ. Thus, Eq. (66) decomposes the total MSE into the Bayes risk LB , which is the irreducible component given the observed data, and the estimatordependent excess risk Gδ , which is the additional, in-principle avoidable component. This excess risk includes additional losses relative to the conditional Bayes-optimal estimator, such as those arising from estimating the unknown propagation parameters from the observations and from model mismatch. B. Oracle Lower Bound
Because LB is generally intractable under the true generative model, we introduce an oracle estimator that can access to the true measurement locations X in addition to the observed πtrue (evec | D, θ) ∝ pGP (p | X(evec ), θ)pe (evec )pX (X(evec )). data D, and use its risk as a lower bound on LB . The oracle (59) estimator is defined as Thus, in addition to the likelihood of the RSS observations f̂or ≜ Eθ [f∗ | D, X] . (68) and the prior distribution of the positioning errors, the true posterior depends on how plausible the corrected measurement Under the assumed independence relations, once X and the trajectory is under the trajectory-generation law. For a fixed RSS observations p are given, additionally conditioning on
SUBMITTED TO IEEE INTERNET OF THINGS JOURNAL, SEPTEMBER 2026
the reported locations X̃ does not change the conditional distribution of f∗ , and hence f̂or = Eθ [f∗ | X, p] .
(69)
For fixed θ and X, f∗ and p are jointly Gaussian. Define the kernel matrices among the evaluation points, between the evaluation and measurement points, and among the measurement points, respectively, as M
∗ K∗∗ ≜ [kexp,θ (x∗,m , x∗,n )]m,n=1 ,
(70)
K∗X ≜ [kexp,θ (x∗,m , xn )]m=1,··· ,M∗ ,
(71)
n=1,··· ,MN MN KXX ≜ [kexp,θ (xm , xn )]m,n=1 , KX∗ ≜ K⊤ ∗X .
(72) (73)
By Gaussian conditioning, the posterior covariance of the oracle estimator is given by Σor (θ; X) = K∗∗ − K∗X (KXX + σp2 I)−1 KX∗ .
(74)
Therefore, the conditional MSE for a given true measurementlocation matrix X is 1 ℓor (θ; X) ≜ tr [Σor (θ; X)] . (75) M∗ Averaging this conditional MSE under the generative model with θ fixed at its true value gives the oracle risk as Lor (θ) ≜ Eθ [ℓor (θ; X)] .
(76)
Because f̂or is a conditional-mean estimator with access to more information than D alone, the projection property of conditional expectation gives Lor (θ) ≤ LB (θ).
(77)
Because f̂B is the conditional-mean estimator based on D, 2 1 . (81) Uwork (θ) = LB (θ) + Eθ f̂work − f̂B M∗ 2 Therefore, LB (θ) ≤ Uwork (θ),
C. Working-Model Upper Bound and Its Tightness The risk of any estimator constructed solely from the observable information D is no smaller than LB . We therefore introduce a working-model estimator as a reference estimator for the analysis. In the working-model, the reported locations X̃ are treated as a fixed design, and the trajectory density pX (X(evec )) appearing in the true posterior is not used. Accordingly, the working-model posterior is defined as (78)
(82)
and Uwork provides an upper bound on the Bayes risk. This upper-bound relation does not require the working model to coincide with the true generative model or f̂work to be close to the true Bayes-optimal estimator. Combining the true posterior defined above with the working-model posterior in Eq. (78), their density ratio can be written as πtrue (evec | D, θ) = c(D, θ)pX (X(evec )) , (83) πwork (evec | D, θ) where c(D, θ) is independent of evec . Hence, the discrepancy between the Bayes-optimal and working-model predictors arises from the trajectory-density factor omitted in the working model. In particular, if pX (X(evec )) is constant for πwork almost every evec , then πtrue = πwork , and consequently f̂B = f̂work and LB (θ) = Uwork (θ). One sufficient condition is a sensor-wise translation-invariant trajectory law with no active boundary or other absolute-location constraints over the relevant support. On a bounded domain, however, boundary effects can break this invariance; therefore, the equality is not assumed in the numerical evaluation. The same translation invariance also clarifies the information available for estimating the sensor-specific constant offsets from the reported trajectories alone. Without using the RSS observations, the posterior of the positioning-error vector conditioned on the reported locations satisfies p(evec | X̃) ∝ pe (evec )pX (X(evec )) .
Thus, Lor provides a lower bound on the Bayes risk.
πwork (evec | D, θ) ∝ pGP (p | X(evec ), θ)pe (evec ).
8
(84)
Thus, when the trajectory law is translation invariant over the relevant support, the reported trajectory does not update the prior distribution of the constant offset. Absolute offset information must then be supplied by absolute-location constraints or by additional observations. In the proposed method, the RSS observations provide such information through the distancedependent mean and the inter-sensor spatial covariance structure described in Sec. IV. D. Bounds on the Bayes Risk and Excess Risk From Eq. (77) and Eq. (82),
The predictor obtained by marginalizing the positioning errors under this posterior is defined as Z f̂work ≜ µGP (evec ; D, θ)πwork (evec | D, θ)devec . (79) Because this integral generally cannot be evaluated in closed form, it is approximated numerically as described in Sec. VI. For a fixed true value of θ, f̂work is an estimator determined solely by D. Its risk evaluated under the true generative model is defined as 2 1 Uwork (θ) ≜ Eθ f∗ − f̂work . (80) M∗ 2
Lor (θ) ≤ LB (θ) ≤ Uwork (θ)
(85)
holds. Therefore, since Gδ (θ) = Rδ (θ) − LB (θ) for any estimator δ, +
(Rδ (θ) − Uwork (θ))
≤ Gδ (θ) ≤ Rδ (θ) − Lor (θ)
(86)
is obtained, where (·)+ denotes the positive part. Thus, without directly evaluating the true Bayes-optimal estimator f̂B or the Bayes risk LB , the deviation of an arbitrary estimator constructed from the observed data D from the Bayes-optimal performance can be bounded.
SUBMITTED TO IEEE INTERNET OF THINGS JOURNAL, SEPTEMBER 2026
9
TABLE I: Reference condition. Parameter Simulation area A Reference distance d0 BS location xTx Transmit power PTx Path-loss index η Shadowing variance σf2 Correlation distance dcor Measurement-noise variance σp2 Number of sensors N Measurement interval Observation duration per sensor Learning rate of θ Learning rate of evec Lévy walk parameter α Lévy walk parameter β
Value 300 × 300 [m] 1.0 [m] [0, 150] 10 [dBm] 3.0 64 [dB2 ] 20 [m] 1.0 [dB2 ] 20 20 [s] 1800 [s] 0.24 0.6 0.5 1.0
To specialize Eq. (86) to the two estimators considered in this study, let Rprop (θ) and Gprop (θ) denote the risk and excess risk, respectively, of the proposed method in Sec. IV, and let Rbase (θ) and Gbase (θ) denote those of the positionerror-agnostic baseline in Sec. III. Applying Eq. (86) to these two estimators yields their corresponding excess-risk intervals. These intervals are evaluated in Sec. VI. VI. N UMERICAL E VALUATION This section evaluates the radio map estimation performance of the proposed method through 1000 independent numerical trials. Using the Oracle lower bound and working-model upper bound derived in Sec. V, we also evaluate the excess risks of the proposed method and the position-error-agnostic baseline. The parameters listed in Table I are used as the reference condition3 . The true measurement trajectories are generated using Lévy walks [26], and the sensor-specific quasi-static positioning errors are generated with µs = 0 and Σs = diag[100, 100]. The MSE is evaluated on a 50 × 50 grid covering the central 150 m × 150 m region. Before applying any of the estimation methods, the observations were sorted by sensor index and time. To avoid numerical instability caused by closely spaced GP inputs, we apply global minimumdistance thinning with a threshold of 7.5 m. Sensors with fewer than four retained measurements after thinning are excluded from radio-map construction4 . Fig. 3 shows an example of the simulations. To evaluate the dependence on propagation conditions, we conduct experiments under various settings in which the propagation parameters are varied from their reference values. The proposed method, the position-error-agnostic baseline, NIGP, which treats input-location uncertainty as effective noise, and Ideal GPR, which uses the true measurement locations, are mainly considered in the evaluation. The working-model estimator introduced in Sec. V is numerically evaluated using 50 trials randomly selected from the 1000 trials for each propagation condition. In the working 3 In practical scenarios, such as constructing public Wi-Fi radio maps, there are many cases where radio maps must be constructed for locations that are sufficiently far from the base station. Therefore, the base station was placed at [0, 150], and the central region was designated as the evaluation region. 4 Since 7.5 m is sufficiently small relative to the [300 × 300] m simulation domain, the effect of thinning is limited in practice.
Evaluation region
(a)
(b)
(c)
Fig. 3: Example simulation: (a) sensor trajectories, (b) observed RSS, and (c) ground-truth radio map in the evaluated region. model, the true propagation parameters θ are assumed to be known, and the posterior πwork in Eq. (78) is constructed by treating the reported locations as a fixed design. Because πwork is a non-Gaussian distribution over R2N , the marginalization integral defining f̂work cannot be evaluated in closed form. We therefore introduce the tempering path from the positioningerror prior pe to πwork as λ
πλ (evec ) ∝ pGP (p | X(evec ), θ) pe (evec ),
λ ∈ [0, 1], (87) and approximate f̂work using an adaptive tempered sequential Monte Carlo (SMC) sampler [27], [28]. Since evec is a fixeddimensional latent variable rather than a time-series state, the particle population is updated along increasing values of λ following the sequential-sampling framework for static models [29]. Each SMC replica uses Np = 256 particles. The temperature sequence {λt } is adaptively selected to target a conditional effective sample size of 0.9Np and systematic resampling is performed whenever the ordinary effective sample size falls below 0.5Np For each trial, eight independent replicas are executed, and their results are averaged to estimate Uwork . The Oracle risk Lor , the working-model risk Uwork , and the risks Rprop and Rbase are all averaged over the same 50 trials, and each endpoint of the excess-risk intervals derived in Sec. V is evaluated as a paired difference on these common trials. A. Radio Map Accuracy Performance The cumulative distribution function (CDF) of the radio map estimation MSE under the reference condition is shown in Fig. 4. Ideal GPR denotes radio map estimation using the true measurement locations, Proposed denotes the proposed method, NIGP denotes noisy-input GP, Without Correction denotes the position-error-agnostic baseline, and Oracle Tuned Without Correction denotes the same position-error-agnostic GP model whose objective is directly chosen to minimize the MSE against the true radio map5 . Compared with Ideal GPR, the proposed method exhibits an increase in MSE of approximately 3.26 dB2 , whereas Without Correction and NIGP exhibit increases of approximately 10.0 dB2 . Even Oracle Tuned Without 5 This method has access to the true radio map and therefore involves data leakage. Although it cannot be implemented in practice, it is included as a reference for assessing the representational limitation of the position-erroragnostic model.
SUBMITTED TO IEEE INTERNET OF THINGS JOURNAL, SEPTEMBER 2026
distance-dependent mean and the pairwise covariance structure used for position calibration. Consequently, the advantage of the proposed method gradually decreases with measurement noise and almost vanishes at σp = 6 dB (i.e., σp2 = 36 dB2 ). Overall, these results indicate that jointly exploiting the mean and spatial covariance structures enables effective position calibration over a broad range of propagation conditions, provided that the RSS measurements retain sufficient spatial information.
1.0 0.8 CDF
10
0.6 0.4 0.2
B. Deviation from the Bayes-Optimal Performance
0.0
20
30
40
MSE [dB2 ] Proposed Without Correction NIGP
Ideal GPR Oracle Tuned Without Correction
Fig. 4: CDF of the radio-map 1 estimation MSE under the reference condition.
Correction, which is allowed to exploit data leakage, exhibits an increase of approximately 9.47 dB2 , leaving a clear gap from the proposed method. These results demonstrate the effectiveness of the proposed method under the reference condition. Next, Fig. 5 shows the median MSE when η, σf , dcor , and σp are varied from the reference condition. The proposed method consistently outperformed Without Correction over almost all evaluated conditions. In Fig. 5(a), the performance of the proposed method is relatively insensitive to η. The sensitivity of the distance-dependent mean to a positioning offset is proportional to η; hence, for a small η, the mean provides weaker absolute-location information and its variation is less distinguishable from stochastic shadowing. The relatively stable performance of the proposed method therefore indicates that the pairwise spatial covariance provides complementary information for location calibration when the mean structure is less informative. In Fig. 5(b), the MSE increases with σf even for Ideal GPR, reflecting the increased intrinsic variation of the shadowing field. A larger shadowing amplitude also increases the effect of spatial misregistration, which explains the substantial degradation of the methods that do not explicitly correct the measurement locations. As shown in Fig. 5(c), increasing dcor improves the performance of all methods because shadowing observations remain informative over a wider spatial range. Moreover, a larger dcor increases the number of inter-sensor measurement pairs with non-negligible covariance, providing additional pairwise constraints on the relative sensor offsets exploited by the proposed method. Finally, in Fig. 5(d), increasing σp increases the diagonal measurement-noise term σp2 I in the observation covariance matrix. When this term dominates the spatial covariance, the RSS observations become less informative about both the
Fig. 6 complements the absolute-MSE results in Fig. 5 by quantifying what fraction of each method’s MSE is attributable to excess risk, i.e., error that could in principle be reduced by a Bayes-optimal use of the same observations. Under the reference condition (Fig. 6(a), η = 3), the bound interval for Gprop /Rprop is approximately 5–16%, whereas that of Gbase /Rbase is approximately 32–40%. Thus, only a relatively small fraction of the proposed method’s MSE is attributable to estimator suboptimality, whereas roughly one third or more of the baseline MSE remains excess error. This result indicates that the proposed method not only achieves a lower absolute MSE but also exploits the information contained in the observations substantially more effectively. To interpret the ordinate, recall that the total risk of an estimator δ is decomposed as Rδ = LB + Gδ . Here, LB is the irreducible error remaining even when the observed data are used optimally under the assumed model, whereas Gδ is the additional, in-principle avoidable error caused by deviation from the Bayes-optimal estimator. Hence, Gδ /Rδ represents the fraction of the total MSE attributable to this avoidable component, and a smaller value indicates performance closer to the Bayes-optimal benchmark. Because LB cannot be evaluated directly, we use Eq. (86): the lower and upper bounds of Gδ /Rδ are (Rδ − Uwork )+ /Rδ and (Rδ − Lor )/Rδ , respectively. Accordingly, the shaded regions in Fig. 6 represent bound intervals rather than statistical confidence intervals. The normalized excess risk of the proposed method remains relatively stable with η in Fig. 6(a), indicating that the spatial covariance provides complementary information when the distance-dependent mean is less informative at small η. Similarly, in Fig. 6(b), the normalized excess-risk intervals change little with σf , even though the absolute MSE increases substantially, showing that larger shadowing variance primarily changes the difficulty of the estimation problem rather than the relative efficiency of the estimators. In Fig. 6(c), the difference between the two methods becomes larger with dcor because the wider correlation range provides more intersensor measurement pairs that carry useful information about relative positioning offsets, which is exploited by the proposed method. Finally, in Fig. 6(d), the two intervals converge as σp increases. This convergence does not indicate improved absolute accuracy of the baseline; rather, strong measurement noise both suppresses the spatial information available for position calibration and increases the irreducible component of the estimation error. Consequently, the benefit of position correction becomes negligible in the noise-dominated regime.
SUBMITTED TO IEEE INTERNET OF THINGS JOURNAL, SEPTEMBER 2026
(a)
(b)
11
(c)
(d)
Fig. 5: Effects of radio-propagation parameters on the median MSE: (a) η, (b) σf , (c) dcor , and (d) σp .
(a)
(b)
(c)
(d)
Fig. 6: Effects of radio propagation parameters on the excess risk ratio: (a) η, (b) σf , (c) dcor , and (d) σp . Overall, Figs. 5 and 6 together show that the proposed method is particularly effective when the RSS observations retain informative spatial structure: it achieves both a lower total MSE and a substantially smaller fraction of avoidable error than the position-error-agnostic baseline. VII. E VALUATION BASED ON M EASURED GNSS E RROR M ODELS This section evaluates the robustness of the proposed method to temporally correlated positioning-error components that are not explicitly represented by the sensor-specific constant-offset model. To construct measurement-informed positioning-error processes, we identify AR(2) models from GNSS positioning errors measured with two smartphones and use the resulting models to generate GNSS-based observation trajectories. Because the numerical experiments in Sec. VI use a 20 s measurement interval, the AR(2) models identified from the 1 Hz measurements are converted to approximate 20 s models before simulation. We also compare the proposed method with a KF–RTS trajectory-smoothing baseline that explicitly exploits the temporal correlation of the positioning errors. Details of the GNSS measurements and model validation, the conversion of the 1 Hz models to model AR(2) models sampled at 20 s intervals, and the complete state-space formulation of the KF–RTS baseline are given in Appendices A and B. Note that the RSS observations and true trajectories remain simulated; the measured GNSS data are used to construct the positioning-error models.
A. GNSS-Derived Error Model and Evaluation Setup Let xref,k ∈ R2 and xobs,k ∈ R2 denote the reference and smartphone-reported GNSS positions, respectively, at the 1 Hz time index k. The measured positioning error is decomposed into a quasi-static component a and a temporally varying residual uk as eobs,k ≜ xobs,k − xref,k = a + uk .
(88)
For each axis q ∈ {x, y}, the residual is modeled by an AR(2) process, ϵq,k ∼ N (0, σq2 ). (89) The quasi-static component a corresponds to the sensorspecific positioning offset assumed by the proposed method, whereas uk represents the time-varying component that is not explicitly modeled by the proposed method. The error models were identified from GNSS measurements collected using an RTK reference receiver and two smartphones, a motorola edge 40 neo and a moto g24. The detailed measurement configuration, observed trajectories, identified AR(2) parameters, and validation results are reported in Appendix A. Although the fitted AR(2) models do not reproduce the complete marginal distributions of the measured errors, they approximately reproduce the strength and decay characteristics of the measured autocorrelation, which are the temporal characteristics of interest in this evaluation. uq,k = ωq,1 uq,k−1 + ωq,2 uq,k−2 + ϵq,k ,
SUBMITTED TO IEEE INTERNET OF THINGS JOURNAL, SEPTEMBER 2026
To make the identified error models consistent with the 20 s measurement interval used in Sec. VI, each 1 Hz AR(2) model is approximated by a 20 s AR(2) model that matches the stationary autocovariances at lags 0, 20, and 40 s. The corresponding Yule–Walker construction is given in Appendix A. As in Sec. VI, the true trajectories of N = 20 sensors are generated independently in each trial. Because the measured values of a are particular realizations associated with the individual devices, they are not assigned identically to all simulated sensors. Instead, the quasi-static positioning bias of sensor i is generated as b(i) ∼ N (0, 100I) .
(90)
12
(a)
(b)
Fig. 7: CDF of radio-map MSE with time-correlated positioning errors: (a) moto g24 and (b) motorola edge 40 neo.
(i)
Let ūn denote an independent residual sequence generated from the corresponding 20 s AR(2) model. The observed position is then generated as (i) (i) x̃(i) + u(i) n . n = xn + b
(91)
The AR(2) states are initialized from their stationary distributions. The prior distribution of the constant positioning bias used by the proposed method is also set to N (0, 100I). Thus, the quasi-static component is statistically matched between the data-generation model and the proposed method, whereas the additional time-varying AR(2) residual constitutes deliberate positioning-error model mismatch. For comparison, we use a KF–RTS baseline whose linear Gaussian state-space model includes the sensor position and velocity, a sensor-specific constant bias, and the AR(2) residual. The baseline is provided with the same 20 s AR(2) coefficients and innovation variances as those used to generate the positioning errors, but the true positions, realized biases, and realized residuals are not provided. Only the GNSS observation trajectories are used in the KF–RTS smoothing stage; the RSS observations are introduced only in the subsequent GPR stage. The complete state-space model and initialization are given in Appendix B. B. Radio Map Estimation Performance under Temporally Correlated Positioning Errors Using the setup in Sec. VII-A, we conduct 1000 independent trials for each device under the same propagation parameters and true trajectory-generation conditions as the reference condition in Sec. VI. Because the proposed method estimates only a sensor-specific constant positioning bias, the time-varying AR(2) residual in Eq. (91) constitutes deliberate model mismatch. The CDFs of the radio map MSE obtained by the different methods are shown in Fig. 7. As shown in Fig. 7, although KF + RTS Smoother is provided with the AR(2) dynamics used to generate the positioning errors, its radio map estimation performance under the present evaluation conditions is close to that of Without Correction. This behavior is consistent with the translation ambiguity described in Sec. V-C and formalized by the state-space model in Appendix B. The KF + RTS smoother can exploit the temporal correlation to estimate the time-varying AR(2) residual, but the GNSS trajectory likelihood cannot separately identify a sensor-specific quasi-static bias from a constant translation
of the true trajectory. Consequently, the separation of these two components is determined only by the absolute-position and bias priors, and a residual sensor-wise translation can remain in the smoothed trajectory. Since the resulting positions are subsequently used as GPR inputs, this residual spatial misregistration limits the improvement over Without Correction. In contrast, although the proposed method does not explicitly model the time-varying residual, it achieves lower MSE than both Without Correction and KF + RTS Smoother for the error models of both devices. The proposed method uses spatial information from the RSS observations that is unavailable to the trajectory-only smoother: the distance-dependent mean provides information about absolute placement relative to the known BS, while the inter-sensor covariance provides information about relative sensor offsets. These RSS-derived constraints provide additional information for estimating the quasi-static positioning offsets even in the presence of the unmodeled time-varying residual. At the median of the CDF, the MSE of the proposed method increases by only approximately 5 dB2 relative to Ideal GPR, whereas the other non-ideal methods exhibit increases of approximately 10 dB2 or more. These results indicate that, under the trajectory-generation conditions and temporalcorrelation models considered here, modeling the temporal structure of the positioning errors alone is insufficient to reliably remove the sensor-specific quasi-static displacement, whereas incorporating the spatial structure of the RSS observations provides effective additional information for post-hoc location calibration. Importantly, the GNSS trajectory likelihood of this baseline (i) (i) is invariant to a sensor-wise transformation pn 7→ pn + c(i) and b(i) 7→ b(i) − c(i) for any constant c(i) ∈ R2 . Hence, except for the absolute-location information provided by the initial-position and bias priors, the GNSS trajectory alone cannot distinguish a constant translation of the trajectory from the corresponding constant positioning bias. This structural ambiguity is consistent with the trajectory-only identifiability discussion in Sec. V-C. VIII. C ONCLUSION This paper addressed radio map construction in the presence of temporally persistent, sensor-specific positioning errors. We
SUBMITTED TO IEEE INTERNET OF THINGS JOURNAL, SEPTEMBER 2026
A PPENDIX A D ETAILS OF THE GNSS E RROR M ODEL AND I TS 20 S A PPROXIMATION This appendix details the GNSS measurements and the construction of the 20-s AR(2) error models used in Sec. VII. We first describe the GNSS measurement procedure and validate the AR(2) models identified from the measured positioning errors. Then, the construction of the approximate 20 s AR(2) models, used in the numerical evaluation, is described. A. GNSS Measurements and AR(2) Model Validation The GNSS measurements were collected on the campus of The University of Electro-Communications on June 6, 2025. The reference trajectory was obtained by RTK positioning using a u-blox ZED-F9R, while the GNSS observation trajectories were collected using two Motorola smartphones, a
moto g24
ZED-F9R and antenna
motorola edge 40 neo
Fig. 8: Measurement setup. −37800
Northing [m]
developed a GPR-based radio map construction method with post-hoc location calibration, which jointly estimates sensorspecific position offsets and radio propagation parameters from the RSS observations and constructs the radio map using the calibrated measurement locations. By exploiting both the distance-dependent path-loss structure and the spatial correlation of shadowing, the proposed method uses the radio measurements themselves as additional information for correcting the locations at which they were collected. Numerical evaluations demonstrated the effectiveness of the proposed method over a wide range of propagation conditions. Under the reference condition, the proposed method exhibited an MSE gap of approximately 3.26 dB2 relative to Ideal GPR, whereas the position-error-agnostic and NIGP methods exhibited gaps of approximately 10 dB2 . The benefit of location calibration was particularly pronounced in low-tomoderate measurement-noise regimes, where the RSS observations retain sufficient spatial information for estimating the positioning offsets. Further, results of the theoretical analysis based on a conditional MMSE benchmark showed that the fraction of MSE attributable to excess risk was approximately 24–27 percentage points lower for the proposed method than for the position-error-agnostic baseline. Finally, we evaluated robustness to positioning-error model mismatch using GNSS-based error sequences derived from measured GNSS data. Although the proposed method assumes a constant positioning offset for each sensor and does not explicitly model the remaining time-varying error component, it consistently outperformed both the position-erroragnostic baseline and the KF + RTS smoother baseline for the GNSS error models considered. These results demonstrate that exploiting the spatial structure of RSS for both radiomap construction and post-hoc location calibration is effective under temporally correlated positioning errors. An interesting direction for future work is to incorporate absolute-location constraints, such as road networks, walkable regions, or known anchors, through pX and the corresponding trajectory-density term in the calibration objective, thereby breaking the translation symmetry and providing additional information about the sensor-specific offsets.
13
−37900
−38000
−38100 −26500 −26400 −26300 −26200
Easting [m]
Buildings u-blox ZED-F9R (RTK)
motorola edge 40 neo moto g24
Fig. 9: Observed trajectory (the1building data is licensed under ©OpenStreetMap contributors).
motorola edge 40 neo and a moto g24. The motorola edge 40 neo is equipped with a dual-frequency GNSS module, whereas the moto g24 uses a single-frequency GNSS module. Positions were recorded at 5 Hz by the ZED-F9R and at 1 Hz by the smartphones, yielding 1249 observed positions from the motorola edge 40 neo and 1251 from the moto g24. The measurement setup and observed trajectories are shown in Figs. 8 and 9, respectively. The parameters of the error model in Eqs. (88) and (89) identified for the two devices are summarized in Table II. Here, a denotes the quasi-static offset obtained from each measured sequence, whereas the AR(2) coefficients and innovation standard deviations characterize the 1 Hz time-varying residual. The histograms of the measured positioning errors and those generated from the identified AR(2) models are shown in Fig. 10, and the corresponding autocorrelation functions (ACFs) are shown in Fig. 11. Although differences between the measured and generated marginal distributions can be observed, the AR(2) models approximately reproduce both the strength and decay characteristics of the measured autocorrelation. Thus, the AR(2) models are not intended to reproduce the complete
SUBMITTED TO IEEE INTERNET OF THINGS JOURNAL, SEPTEMBER 2026
motorola edge 40 neo [−6.37, 1.00]⊤ 1.183 −1.947 × 10−1 4.836 × 10−1 1.246 −2.621 × 10−1 5.258 × 10−1
moto g24 [−1.91, 1.99]⊤ 1.491 −5.084 × 10−1 3.130 × 10−1 1.405 −4.144 × 10−1 4.625 × 10−1
1.0
0.8
0.8
0.6
0.6
ACF
Parameter a [m] ωx,1 ωx,2 σx [m] ωy,1 ωy,2 σy [m]
1.0
ACF
TABLE II: Parameters of 1 Hz AR(2) models identified from measured GNSS positioning errors.
14
0.4
0.4
0.2
0.2
0.00
5
10
Lag
Observed 𝑥 ACF
15
20
0.00
5
1
0.10
0.00 −20
0
Error 𝑥 [m]
Observed 𝑥 Error
10
20 0.00
−10
0
Error 𝑦 [m]
Observed 𝑦 Error
Simulated 𝑥 Error
0.8
0.6
0.6 0.4
0.2
0.05 −10
1.0
0.8
0.4
0.10
0.05
1.0
ACF
0.15
0.2
0.00
10
5
10
Lag
Observed 𝑥 ACF
Simulated 𝑦 Error
15
20
0.00
5
20
Simulated 𝑦 ACF
1
(c)
(b)
15
10
Lag
Observed 𝑦 ACF
Simulated 𝑥 ACF
1
1
1
(a) 0.20
(b)
ACF
0.15
20
Simulated 𝑦 ACF
1
(a) 0.20
15
10
Lag
Observed 𝑦 ACF
Simulated 𝑥 ACF
(d)
Fig. 11: ACF of positioning error (a) motorola edge 40 neo (x-axis), (b) motorola edge 40 neo (y-axis), (c) moto g24 (xaxis), (d) moto g24 (y-axis).
0.15
0.15
0.10 0.10
0.05
0.05 0.00 −20
−10
0
Error 𝑥 [m]
Observed 𝑥 Error
10
20 0.00
Simulated 𝑥 Error
−10
0
10
Error 𝑦 [m]
Observed 𝑦 Error
20
Simulated 𝑦 Error
1
1
(c)
(d)
Fig. 10: Histograms of positioning error: (a) motorola edge 40 neo (x-axis), (b) motorola edge 40 neo (y-axis), (c) moto g24 (x-axis), (d) moto g24 (y-axis). distribution of the GNSS positioning errors; rather, they provide measurement-informed approximations of the temporal correlation considered in Sec. VII. The measured sequences also show that strongly temporally correlated components remain after accounting for the quasi-static offset, motivating the model-mismatch evaluation in that section. B. Construction of the 20 s AR(2) Model The AR(2) models in Table II were identified from GNSS positioning-error sequences sampled at 1 Hz, whereas the numerical evaluations in Sec. VI use a measurement interval of 20 s. Directly using the 1 Hz AR(2) coefficients at the 20 s sampling interval would not preserve the temporal correlation of the original process. We therefore construct an approximate AR(2) model at 20 s intervals by matching selected stationary autocovariances of the 1 Hz model. For each axis q ∈ {x, y}, let the stationary autocovariance of the 1 Hz AR(2) process be γq (ℓ) ≜ E [uq,k uq,k−ℓ ] .
(92)
The approximate AR(2) process at 20 s intervals is written as ūq,n = ω q,1 ūq,n−1 + ω q,2 ūq,n−2 + ϵq,n ,
ϵ̄q,n ∼ N (0, σ̄q2 ). (93)
To preserve the stationary autocovariances at lags 0, 20, and 40 s, the AR coefficients are obtained from the Yule–Walker equations [30], γq (0) γq (20) ω q,1 γ (20) = q . (94) γq (20) γq (0) ω q,2 γq (40) The corresponding innovation variance is σ̄q2 = γq (0) − ω q,1 γq (20) − ω q,2 γq (40).
(95)
This construction is a moment-matching approximation and does not assume that exact 20 s subsampling of the original 1 Hz AR(2) process remains an AR(2) process. The resulting 20 s models are used to generate the residual sequences u(i) n in Eq. (91). A PPENDIX B S TATE -S PACE F ORMULATION OF THE KF–RTS BASELINE This appendix gives the complete state-space formulation of the KF–RTS baseline used in Sec. VII. The baseline explicitly represents the temporally correlated AR(2) residual in a linear Gaussian state-space model, while using only the GNSS observation trajectory during the smoothing stage. The resulting smoothed positions are subsequently used as the training locations for GPR. For sensor i, the state vector at time n is h (i) (i) (i) (i) (i) (i) s(i) n = px,n , vx,n , py,n , vy,n , bx , by , i⊤ (i) (i) (i) u(i) . (96) x,n , ux,n−1 , uy,n , uy,n−1 (i)
(i)
(i)
(i)
Here, pn = [px,n , py,n ]⊤ and vn denote the true position and velocity, respectively, b(i) denotes the sensor-specific (i) constant positioning bias, and un denotes the AR(2) residual.
SUBMITTED TO IEEE INTERNET OF THINGS JOURNAL, SEPTEMBER 2026
For each axis q ∈ {x, y}, the true position and velocity follow the constant-velocity model " # " (i) # (i) pq,n 1 ∆t pq,n−1 (i) = + wq,n . (97) (i) (i) 0 1 vq,n vq,n−1 The process noise is modeled as 3 ∆t /3 (i) wq,n ∼ N 0, qmot ∆t2 /2
∆t2 /2 ∆t
,
(98)
where qmot = 0.30 m2 /s3 is used in this evaluation. The constant positioning bias evolves according to (i)
b(i) n = bn−1 .
(99)
(i)
(i)
(i)
For the AR(2) residual, define rq,n = [uq,n , uq,n−1 ]⊤ . Using the 20 s AR(2) parameters obtained in Appendix A, we write ω̄q,1 ω̄q,2 (i) (i) rq,n−1 + ξ (i) (100) rq,n = q,n , 1 0 {z } | Aq
where ξ (i) q,n ∼ N (0, Qq ) ,
Qq ≜
2 σ̄q 0
0 . 0
(101)
When thinning increases the interval between consecutive retained measurements to ∆t = 20K s, the state-transition matrix is AK q , and the accumulated process covariance is Qq,K =
K−1 X
Ajq Qq Ajq
⊤
.
(102)
j=0
The observation equation is # " (i) (i) (i) px,n + bx + ux,n (i) + ν (i) x̃n = (i) n . (i) (i) py,n + by + uy,n
(103)
Because the simulated observations in Eq. (91) contain no additional white positioning noise independent of the AR(2) (i) residual, we set ν n = 0. For the initial position, the uniform distribution over the 300 m × 300 m simulation area is approximated by a Gaussian distribution through moment matching, E[p0 ] = [150, 150]⊤ ,
Cov(p0 ) = 7500I.
(104)
We further use b(i) ∼ N (0, 100I), (105) and initialize the AR(2) states from their stationary distributions. A forward KF is first applied to the state-space model, followed by fixed-interval RTS smoothing over the entire (i) observation sequence. The smoothed positions p̂n,RTS are then used as the training locations of the subsequent GPR. In this comparison, the KF–RTS smoother is provided with the same 20 s AR(2) coefficients and innovation variances as those used to generate the positioning errors. However, the true positions, true constant biases, and realized AR(2) residuals in each trial are not provided to the baseline. Only the GNSS observation E[v0 ] = 0,
Cov(v0 ) = 2I,
15
trajectories are used for position smoothing, and the RSS observations are introduced only in the subsequent GPR stage. The transition and observation likelihoods of this statespace model are invariant to the sensor-wise transformation (i) (i) pn 7→ pn + c(i) and b(i) 7→ b(i) − c(i) for any constant (i) vector c ∈ R2 . Therefore, apart from the absolute-position information introduced through the initial-position and bias priors, the GNSS trajectory likelihood alone cannot distinguish a constant translation of the trajectory from the corresponding sensor-specific positioning bias. Even when the temporal AR(2) dynamics are correctly specified, temporal smoothing alone therefore provides no additional likelihood information for resolving this quasi-static translation ambiguity. This is consistent with the trajectory-only identifiability discussion in Sec. V-C. R EFERENCES [1] K. Kanzaki and K. Sato, “Joint Ex-Post Location Calibration and Radio Map Construction under Biased Positioning Errors,” in 2025 IEEE Globecom Workshops, Dec. 2025, pp. 2071–2076. [2] D. Romero and S.-J. Kim, “Radio Map Estimation: A data-driven approach to spectrum cartography,” IEEE Signal Process. Mag., vol. 39, no. 6, pp. 53–72, Nov. 2022. [3] S. Bi, J. Lyu, Z. Ding, and R. Zhang, “Engineering Radio Maps for Wireless Resource Management,” IEEE Wireless Commun., vol. 26, no. 2, pp. 133–141, Apr. 2019. [4] X. Zhang, W. Sun, J. Zheng, A. Lin, J. Liu, and S. S. Ge, “Wi-FiBased Indoor Localization With Interval Random Analysis and Improved Particle Swarm Optimization,” IEEE Trans. Mobile Comput., vol. 23, no. 10, pp. 9120–9134, Oct. 2024. [5] A. C. Suarez Rodriguez, N. Haider, Y. He, and E. Dutkiewicz, “Network Optimisation in 5G Networks: A Radio Environment Map Approach,” IEEE Trans. Veh. Technol., vol. 69, no. 10, pp. 12 043–12 057, Oct. 2020. [6] D.-H. Tran, N. Waheed, Y. M. Saputra, X. Lin, C. T. Nguyen, T. S. Abdu, V. N. Vo, V.-Q. Pham, M. Alsenwi, A. B. M. Adam, S. Chatzinotas, E. Lagunas, H. Tran, T. H. Dac, and N. V. Huynh, “Network Digital Twin for 6G and Beyond: An End-to-End View Across Multi-Domain Network Ecosystems,” IEEE Open J. Commun. Soc., vol. 6, pp. 6866– 6911, 2025. [7] H. Wang, J. Zhang, G. Nie, L. Yu, Z. Yuan, T. Li, J. Wang, and G. Liu, “Digital Twin Channel for 6G: Concepts, Architectures and Potential Applications,” IEEE Commun. Mag., vol. 63, no. 3, pp. 24–30, Mar. 2025. [8] Y. Zeng, J. Chen, J. Xu, D. Wu, X. Xu, S. Jin, X. Gao, D. Gesbert, S. Cui, and R. Zhang, “A Tutorial on Environment-Aware Communications via Channel Knowledge Map for 6G,” IEEE Commun. Surveys Tuts., vol. 26, no. 3, pp. 1478–1519, 2024. [9] G. Chen, Y. Liu, J. Zhang, T. Zhang, K. Liu, and J. Yang, “GPRT: A Gaussian Process Regression-Based Radio Map Construction Method for Rugged Terrain,” IEEE Internet Things J., vol. 12, no. 13, pp. 23 905– 23 920, Jul. 2025. [10] H. Sun and J. Chen, “Propagation Map Reconstruction via Interpolation Assisted Matrix Completion,” IEEE Trans. Signal Process., vol. 70, pp. 6154–6169, 2022. [11] Y. Teganya and D. Romero, “Deep Completion Autoencoders for Radio Map Estimation,” IEEE Trans. Wireless Commun., vol. 21, no. 3, pp. 1710–1724, Mar. 2022. [12] C. E. Rasmussen and C. K. I. Williams, Gaussian Processes for Machine Learning. Cambridge, Mass.: MIT Press, 2005. [13] N. Zhu, J. Marais, D. Bétaille, and M. Berbineau, “GNSS Position Integrity in Urban Environments: A Review of Literature,” IEEE Trans. Intell. Transp. Syst., vol. 19, no. 9, pp. 2762–2778, Sep. 2018. [14] H. Xiong, R. Bian, Y. Li, Z. Du, and Z. Mai, “Fault-Tolerant GNSS/SINS/DVL/CNS Integrated Navigation and Positioning Mechanism Based on Adaptive Information Sharing Factors,” IEEE Syst. J., vol. 14, no. 3, pp. 3744–3754, Sep. 2020. [15] P. Ranacher, R. Brunauer, W. Trutschnig, S. Van der Spek, and S. Reich, “Why GPS makes distances bigger than they are,” Int. J. Geogr. Inf. Sci., vol. 30, no. 2, pp. 316–333, Feb. 2016.
SUBMITTED TO IEEE INTERNET OF THINGS JOURNAL, SEPTEMBER 2026
[16] A. Mchutchon and C. Rasmussen, “Gaussian Process Training with Input Noise,” in Adv. Neural Inf. Process. Syst., vol. 24. Curran Associates, Inc., 2011. [17] A. Girard, C. Rasmussen, J. Q. Candela, and R. Murray-Smith, “Gaussian Process Priors with Uncertain Inputs Application to Multiple-Step Ahead Time Series Forecasting,” in Adv. Neural Inf. Process. Syst., vol. 15. MIT Press, 2002. [18] A. C. Damianou, M. K. Titsias, and N. D. Lawrence, “Variational Inference for Latent Variables and Uncertain Inputs in Gaussian Processes,” J. Mach. Learn. Res., vol. 17, no. 42, pp. 1–62, 2016. [19] P. Zhen, B. Zhang, Y.-Q. Xu, Z. Chen, H. Wang, and D. Guo, “Radio Environment Map Construction Based on Gaussian Process With Positional Uncertainty,” IEEE Wireless Commun. Lett., vol. 11, no. 8, pp. 1639–1643, Aug. 2022. [20] R. E. Kalman, “A New Approach to Linear Filtering and Prediction Problems,” J. Basic Eng., vol. 82, no. 1, pp. 35–45, Mar. 1960. [21] S. Särkkä and L. Svensson, Bayesian Filtering and Smoothing. Cambridge Univ. Press, 2023. [22] K. Iyer, A. Dey, B. Xu, N. Sharma, and L.-T. Hsu, “Enhancing Positioning in GNSS Denied Environments Based on an Extended Kalman Filter Using Past GNSS Measurements and IMU,” IEEE Trans. Veh. Technol., vol. 73, no. 6, pp. 7908–7924, Jun. 2024. [23] A. Goldsmith, Wireless Communications. Cambridge Univ. Press, Aug. 2005. [24] M. Gudmundson, “Correlation model for shadow fading in mobile radio systems,” Electron. Lett., vol. 27, no. 23, pp. 2145–2146, Nov. 1991. [25] D. P. Kingma and J. Ba. Adam: A Method for Stochastic Optimization. [Online]. Available: http://arxiv.org/abs/1412.6980 [26] I. Rhee, M. Shin, S. Hong, K. Lee, S. J. Kim, and S. Chong, “On the Levy-Walk Nature of Human Mobility,” IEEE/ACM Trans. Netw., vol. 19, no. 3, pp. 630–643, Jun. 2011. [27] R. M. Neal, “Annealed importance sampling,” Stat. Comput., vol. 11, no. 2, pp. 125–139, Apr. 2001. [28] P. Del Moral, A. Doucet, and A. Jasra, “Sequential Monte Carlo Samplers,” J. R. Stat. Soc. Ser. B Stat. Methodol., vol. 68, no. 3, pp. 411–436, Jun. 2006. [29] N. Chopin, “A Sequential Particle Filter Method for Static Models,” Biometrika, vol. 89, no. 3, pp. 539–551, 2002. [30] B. Friedlander and B. Porat, “The Modified Yule-Walker Method of ARMA Spectral Estimation,” IEEE Trans. Aerosp. Electron. Syst., vol. AES-20, no. 2, pp. 158–173, Mar. 1984.
Koki Kanzaki (Graduate Student Member, IEEE) received the B.E. and M.E. degrees in engineering from The University of Electro-Communications in 2024 and 2026, respectively. He is currently pursuing a Ph.D. degree at the same university. He received first place in the ITU AI/ML in 5G Challenge 2025. His research interests include spatio-temporal statistics.
Katsuya Suto (Senior Member, IEEE) received the B.Sc. degree in computer engineering from Iwate University, Morioka, Japan, in 2011, and the M.Sc. and Ph.D. degrees in information science from Tohoku University, Sendai, Japan, in 2013 and 2016, respectively. He has worked as a Post-Doctoral Fellow for Research Abroad, Japan Society for the Promotion of Science, at the Broadband Communications Research Laboratory, University of Waterloo, Waterloo, ON, Canada, from 2016 to 2018. He is an Associate Professor with the Faculty of Information Science and Technology, Hokkaido University, Sapporo, Hokkaido. His research interests include semantic communication, radio propagation, deep learning, and graph representation. Dr. Suto is a member of IEICE. He received the Best Paper Award from the IEEE VTC2013-spring, IEEE/CIC ICCC2015, IEEE ICC2016, and IEEE Transactions on Computers in 2018.
16
Koya Sato (Senior Member, IEEE) received the B.E. degree in electrical engineering from Yamagata University, in 2013, and the M.E. and Ph.D. degrees from The University of Electro-Communications, in 2015 and 2018, respectively. From 2018 to 2021, he was an Assistant Professor with the Tokyo University of Science. He is currently an Assistant Professor with the Artificial Intelligence eXploration Research Center, The University of ElectroCommunications. His current research interests include wireless communication, distributed machine learning, and spatial statistics. Since 2026, he has served as an Associate Editor of IEEE Open Journal of the Communications Society and IEICE Transactions on Communications.