Uncertainty-aware temperature interpolation for measurement rooms using ordinary Kriging

Accurate measurements in dimensional metrology necessitate strict controls on spatial and temporal variations in the measurement room temperature. Due to limitations in the number of sensors that can be placed in a given room, interpolation methods that leverage information from multiple sensors are necessary to assess conditions at unsampled locations. In this contribution, Kriging is used to spatially interpolate room temperatures from a limited number of sensors with different measurement uncertainties in a temperature controlled room housing two coordinate measurement machines. A novel method to propagate sensor uncertainties to the interpolated values using a Monte Carlo simulation is also demonstrated. The uncertainty propagation is considered for the explicitly heteroskedastic, i.e. a constituent network of sensors with different measurement uncertainties. The influence of a localized disturbance in the form of a movable heating element in the room is also investigated.


Introduction
The rapid growth of automation in manufacturing brought about by a widespread adoption of interconnected cyberphysical systems and large-scale sensor networks, often bolstered using machine learning methods, is a key aspect of the industrial internet of things, or IIoT, paradigm [1]. The unique set of challenges posed by the aforementioned * Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. developments are intimately connected to the goals of interoperability and decentralization which are central to the ongoing fourth industrial revolution or Industry 4.0 [2,3]. The use of sensor networks presents both unique advantages and challenges due to their versatility [4]. At the forefront of these challenges is the generation of a large amount of data which in turn necessitates the use of methods to optimally use the resulting information. Often, a form of sensor fusion [5] must be employed to combine measurement results from disparate sensors. By consolidating multi-sensor information in this way, insights otherwise unavailable to individual sensors become accessible. A key application of multi-sensor information is the interpolation of sensor data to estimate the value of a physical property at unsampled locations given a finite number of sensors placed at discrete locations [6]. In this contribution a method for estimating the room temperature using Kriging is presented along with a Monte Carlo based method to estimate the associated measurement uncertainty. The sensor readings correspond to climate-controlled measurement rooms for large coordinate measuring machines (CMMs) that must conform to strict temperature requirements in order to ensure the accuracy of their measurement results.

Motivation
Dimensional metrology is a subbranch of metrology that deals with the measurement and quantification of physical size, shape and distance and the development of methods to aid in this endeavor. CMMs are devices that measure the geometrical properties of an object by means of probes that sense discrete points on its surface [7]. CMMs are required to provide highly precise coordinate measurements and in turn subject their measuring rooms to very stringent requirements with respect to temperature variations both across the dimensions of the room as well as in time in order to ensure the accuracy of their measurement results (see, for instance, [8]). In previous work [9], the interpolation of temperatures for the walls of cold-storage rooms and transport containers yielded promising results with Kriging in comparison to other interpolation methods. Kriging has also been successfully used for temperature interpolation in meteorological contexts [10,11]. In contrast, readings in measurement rooms have to be taken with high accuracy in all three spatial dimensions and, consequently, require a good control of temperature variations in all directions. For reasons of practicality as well as cost, only a limited number of sensors can be placed in the room at a given time. As a result, it becomes necessary to provide accurate estimates for the temperature at locations without nearby sensors. Typically this is achieved using appropriate interpolation methods and in turn necessitates an estimation of the uncertainty associated with the interpolated value in order to quantify its trustworthiness. In contrast to deterministic interpolation methods such as splines, Kriging provides the best linear unbiased estimator [12] for an interpolated value as a linear combination of the known interpolation points based on a prior covariance. As a result, the interpolated value takes the form of a random variable with an associated mean and variance. The relation of this variance to the measurement uncertainty is, however, unclear. The aim of the present work is therefore to demonstrate the applicability of Kriging to 3d interpolation in measurement rooms and to propose a method to compute an uncertainty for the interpolated values consistent with metrological principles.

Use case
The system under study is that of a network of 24 Pt100 platinum resistance thermometers which used to monitor the temperature in a cuboidal measuring room of dimensions 7.51 m × 4.26 m × 2.9 m. The sensors are essential to minimize errors arising from the expansion or contraction of materials due to temperature changes [13]. As illustrated in figure 1, the positioning of the aforementioned sensors is dictated by the presence of essential measuring instruments with fixed locations. The floor positions of the CONTURA CMM and the RONDCOM form tester are also indicated in the figure along with that of the PC housed in the measuring room. The readings of the main group of 24 sensors were recorded at a rate of roughly once every 180 s over a period of 130 d from 9 June 2022 to 17 October 2022 resulting in 62 335 temperature measurements for each of the 24 measuring room sensors. The individual readings along with the respective hourly and daily moving average temperatures are plotted for the entire measurement period in figure 2. An increase in measured values corresponding to higher ambient temperatures in the hotter months of July and August can be clearly seen. A sample from the temperature dataset corresponding to the 200th time point is shown in table 1.

Heating element
In addition to the 24 sensors which are considered as main constituents of the measurement setup, the room was also provided with a heating element and two additional temperature sensors placed on a movable trolley starting on 29 July 2022. The heating element, which was placed on a movable trolley along with two temperature sensors as shown in figure 3, is used to assess the influence of localized sources of disturbances in the measurement room. The first of these sensors is placed directly under the element, while the second sensor is attached to the trolley at a distance of 10 cm from the heating element.
The trolley along with the heating element was placed in three different positions (1, 2 & 3) as illustrated in figure 1 on the floor of the measurement room and the heating element was switched on continuously for three distinct periods given by The trolley setup is present in the room between and outside the aforementioned intervals, although the heating element is switched off. After the heating element has been switched off, the trolley remains in its current position and is moved to its new location shortly before the element is switched on again. Starting on 29 July 2022, both the sensors also record temperatures every three minutes, such that Sensor 25 provides 38 375 and Sensor 26 provides 38 365 measurements respectively. The difference in the number of readings between the two sensors is due to a small delay in recording the temperature values from Sensor 26 and has a negligible effect on the present study. The temperature readings of the two sensors, denoted by Sensor 25 and Sensor 26 are depicted in figure 4. The sensor placed directly on the heating element (Sensor   25) records a temperature of roughly 76.6 • C when the heating element is active. The main purpose served by Sensor 25 is to clearly demarcate the periods during which the heating element is active (or switched on). Sensor 26, on the other hand, will later play a more central role in assessing the interpolation model.

Methods
The present section presents a brief overview of ordinary Kriging and it is applicability and limitations with respect to the problem of interpolating temperatures in measurement rooms, particularly with respect to the heteroscedastic case.

Ordinary Kriging
Kriging [14] is a spatial interpolation technique based on Gaussian processes with origins in geostatistics that, given a finite number of samples, provides predictions for the value of a function at unsampled locations under appropriate assumptions about the mean and covariance of the underlying process. Different variations of Kriging are possible with respect to the chosen constraints on the mean and variance of the underlying process. For instance, while simple and ordinary Kriging assume a known and unknown mean respectively, universal Kriging assumes a general polynomial trend. As no a priori knowledge about the temperature variation with time can be assumed, the mean is taken to be unknown. Furthermore, we assume no knowledge of the underlying physical processes governing the room temperature. As a result, we use ordinary Kriging to determine the temperature at the unsampled points in the room.
Given N observed values (for e.g. sensor readings) with so as to satisfy the unbiasedness condition. The Kriging estimate is determined by minimizing the expected value of the estimation error ϵ(x 0 ) where T 0 is the true value at x 0 and σ ϵ is referred to as the Kriging variance. The above equation can be rewritten as where is the covariance matrix corresponding to the known samples, D i = Cov(x i , x 0 ) a vector whose values are the covariances between the observed values and the function value to be estimated at x 0 and Var sp (T 0 ) = Cov(x 0 , x 0 ) the spatial covariance of the T 0 . The coefficients w that minimize the Kriging variance σ ϵ under the constraint of unbiasedness (2) are the solution of the ordinary Kriging system and are given by [15] The minimal ordinary Kriging variance is given by The final step is to select an appropriate prior covariance function Cov(x i , x j ) in order to determine C and D. Typically, the underlying process is assumed to be wide-sense stationary, , or the covariance between the values at two points only depends on their distance and that the mean of the underlying process is constant. As a result, the variance Var(T 0 ) is given by the value of the spatial covariance at distance 0, so that equation (7) can be rewritten as The covariance function is chosen such that nearby points have a stronger correlation. The covariance is estimated using the variogram γ, which is related to the spatial covariance by and corresponds to the variance of the difference between values separated by a distance h. The variogram is approximated from the sampled values by fitting a chosen function to the empirical variogram where N(h) denotes the set of all samples located at a distance h away from each other within a given tolerance. The discrete values of γ(h) are binned and fitted to a parametric variogram function. Common choices for such a function are linear, Gaussian, exponential and spherical.

Kriging variance and measurement uncertainty
As the variogram and consequently the covariance is fitted directly to the current observed values or the sensor readings, the effect of the individual sensor measurement uncertainty determined from its calibration on the reliability of the interpolated value is not directly taken into account. While the Kriging variance is a reasonable first measure of the reliability of the interpolation, it is an incomplete measure of uncertainty from a metrological standpoint and cannot, by itself, be interpreted as such. A reliability statement obtained by propagating the known uncertainties of all sensors used to determine the kriged estimates to the interpolated values would, however, complement the Kriging variance. In principle, the sensor measurement uncertainty can be incorporated in the variogram (9) via the nugget effect, i.e. a discontinuity in γ at h = 0. For, example, an exponential variogram model is given by The parameters c 0 c and R are determined by fitting the above function to the empirical variogram (10) such that Cov = c 0 + c is the measurement uncertainty. The value of the variogram at h → 0 + , i.e. γ = c 0 + c is referred to as the nugget and the value at h → ∞ is referred to as the sill. The value of c 0 c and R is determined by fitting the chosen variogram model to the empirical variogram described in the preceding section; cf (10). Heuristic methods to incorporate the individual sensor uncertainties in the nugget effect by imposing γ(0) = u 0 , where all measured points have the same uncertainty u 0 , have been implemented [16,17]. These methods however do not conform with the standards for uncertainty propagation put forth by the guide to the expression of uncertainty in measurement (GUM) [18] and moreover cannot account for the case where the sampled locations have different associated uncertainties, i.e. for the heteroscedastic case. As a result, the measurement uncertainty for the interpolated values has to be estimated using other means.

Uncertainty propagation
The GUM prescribes methods for the expression of uncertainty for the measurement of physical quantities and for the propagation of uncertainty to quantities expressed using mathematical models applied to physically measured values. For a quantity . , x N , the "combined" standard uncertainty u c (y) associated with y is given up to second order in u 2 (x i ) by using the Taylor series expansion of f such that, where u(x i ) is the standard uncertainty associated with the quantity x i and u(x i , x j ) is the covariance of the quantities x i and x j . In the absence of significant nonlinearities in f and small uncertainties in x i , the first order term is sufficient to express the uncertainty in y and in the case of uncorrelated inputs, u(x i , x j ) vanishes.
In the context of interpolation using Kriging, the independent measured quantities correspond to the temperature sensor readings in the measurement room and their corresponding standard uncertainties are obtained from the appropriate technical specifications. The task at hand is, therefore, to propagate the aforementioned uncertainties to the interpolated temperatures at each point in the grid. Although the interpolated values are, at first glance (cf equation (1)), linear combinations of the sensor readings, it must be remembered that the coefficients themselves are determined from the individual measurements (see equations (4)- (6)). In other words, the Kriging model is inherently nonlinear in the measured values. Moreover, due to increasing complexity of the model with the number of measurements, uncertainty propagation by means of computing the derivatives as shown in equation (12) is no longer practical.

Monte Carlo uncertainty estimation.
For cases such as the aforementioned, an estimated of the propagated measurement uncertainties can be obtained by using the Monte Carlo method in accordance with the first supplement of GUM [19]. Instead of computing the propagated uncertainty u c (y) by computing the Taylor series expansion of y = f(x 1 , x 2 , . . . , x N ) (12), random perturbations of the measurands x 1 , x 2 , . . . , x N are generated as samples from a multivariate normal distribution X ∼ N (x, C), where x 2 , . . . , x N ) and (13) are, respectively, the mean and covariance of the aforementioned multivariate normal distribution and u(x i ) is the known measurement uncertainty of the measurand x i . This is repeated for multiple iterations or trials such that for a given trial i, the value of the function f for a samplex i ∼ N (x, C) drawn from the multivariate normal distribution N is given bỹ The (unbiased) mean and variance of y after N trials Monte Carlo trials is then calculated as In the present case of interpolation using Kriging, the function f corresponds to using the known temperature readings T i to determine the interpolated temperatureT 0 at a given point in the measuring room (cf equation (1)). However, as shown in section 3.1 Kriging also provides an estimate for the reliability of the interpolated value in the form of a Kriging variance (8).
The mean and variance of the Kriging variance σ 2 OK can itself be determined in an analogous way to the Monte Carlo estimation of the mean and variance of the interpolated temperature in equations (16) and (17).

Total variance as measurement uncertainty.
In order to obtain a meaningful value for the propagated measurement uncertainty for the kriged estimate of the room temperature, it is necessary to incorporate the Kriging variance into the measurement uncertainty determined from calibration. Following the derivation in section 3.1, we argue that for a given set of sensor readings, the estimates for the mean and variance of the interpolated temperature are, in effect, a conditional mean and variance given a particular set of sensor measurements and appropriate assumptions about the spatial variance represented by the variogram. In other words, the Kriging mean and variance (see equations (1) and (4)) can expressed aŝ The variance ofT 0 or, equivalently, its uncertainty u 2 MC (T 0 ) can be estimated using the Monte Carlo method outlined in equations (13)-(17) by drawing samplesT from a multivariate normal distribution N (T, C). The mean value of σ 2 OK can be estimated in a similar manner. These two values can now be combined using the law of total variance [20, p 222] such that, The two terms on the right-hand side respectively correspond to the unexplained and explained components of the variance. In doing so, we have combined the inherent randomness arising from the measurement uncertainty of each sensor with the spatial variation arising from our lack of knowledge about the temperature distribution in the measuring room. In contrast to inserting the measurement uncertainty as part of the nugget effect, the present method supports the heteroscedastic case as no restrictions are placed on the individual sensor uncertainties.

Implementation details
The coordinate reference frame was chosen such that floor of the bottom left corner in figure 1 corresponds to the origin (0, 0, 0). The coordinates of the 24 fixed sensors in our convention is given in table 2. The positions of sensor 25 (attached to the heating element) and sensor 26 on the three time intervals specified in section 2.1 are provided in table 3. The analysis of the temperature data and the numerical interpolation was carried out using the Python programming language. In particular, the tasks directly related to Kriging such as the generation of the variogram, the estimation of the kriged values as well as the Kriging variances were implemented using the PyKrige [21] library. In particular the OrdinaryKriging3D class is used to fit a Gaussian variogram to the empirical variogram (10) from the 24 sensor outputs at a given point in time.
The kriged values and the resulting variances were estimated on a grid of size 30 × 17 × 12 defined in the volume illustrated in figure 1. The grid-size was chosen to ensure a nearly equal resolution in all three directions.

Monte Carlo implementation.
In each trial, a random sample is drawn from the distribution modeling the uncertainty knowledge about the 24 sensor readings. To do so, a multivariate normally distributed random variable with zero mean is used whose covariance is determined from the standard uncertainty of the individual sensors. In reality, the sensors used in the setup are identical and have a measurement uncertainty of 100 mK. In order to demonstrate the propagation of uncertainty for the heteroscedastic case, we assume that three out of the 24 sensors have a lower uncertainty. In particular Sensors 16 and 20 are assigned a measurement uncertainty of 25 mK, while Sensor 6 is assumed to have an uncertainty of 50 mK (cf figure 1 for sensor  T krig , σ 2 krig ← OK3D(X grid , Y grid , Z grid ); 8 Tmean ← Tmean + T krig ; 9 var cond ← var cond + T 2 krig ; 10 vartot ← vartot + σ 2 krig ; 11 end 12 var cond ← var cond − T 2 mean /N trials /(N trials − 1); 13 vartot ← var cond +vartot/N trials ; 14 Tmean ← Tmean/N trials ;

Results
In our subsequent analysis, the exponential variogram will be used to perform the Kriging interpolation. In general, no significant effect on the interpolated values was observed upon using Gaussian or spherical variograms. However, the exponential variogram was chosen due to its relatively straightforward functional form and as it leads to a better illustration of the uncertainty propagation.

Interpolated temperature
The interpolated temperature across the volume of the measurement room is illustrated for six different cases in figures 5 and 6 for sensor readings taken at midday, i.e. 12 PM. In the former set of three figures, the values were chosen to represent three of the colder days from the dataset under consideration, while the latter shows interpolated temperatures for three of the warmer days. The two sets of examples were chosen to increase the clarity of the illustrations, given the differences in their temperature scales. In both cases a region of higher temperature is observed at the bottom right corner of the grid, consistent with the presence of the PC in the same location (cf figure 1). This effect is most pronounced on the colder day of 5 October, as the effect of the heat emanating from the PC is significant compared to the effect of the ambient temperature. Moreover, the heating element is located at Position 3 on the aforementioned date (see section 2.1) and a slightly higher interpolated temperature can be observed at its position. The trolley supporting the heating element with the attached sensors is not present in the room at the time of the first two cases, i.e. on 10 June and 1 July 2022. In contrast to the colder days, the temperature on the warmer days is more significantly influenced by the ambient summer temperatures than by the heat emitted by the PC. The heating element is switched on and located at positions 1 and 2 for the first two cases on 3 and 20 August, while it is switched off at position 2 for the third case, i.e. on 1 September. In all cases, the variation of temperature in the room is small due to the constantly running air-conditioning system and the interpolated temperature  at a point is most strongly determined by that of the nearest sensor.

Effect of the heating element.
In order to assess the performance the Kriging interpolation model at locations without sensors, the temperature determined at the physical location of Sensor 26 attached to the movable trolley is compared with the temperature reading shown by the sensor itself. In figure 7, the interpolated temperature is plotted against the reading of Sensor 26 for all three positions of the trolley, both when the heating element is switched on and off. For reference, the reading of Sensor 26 is also plotted alongside the interpolated value at its location for the entire duration of data collection. The sensor measurements corresponding to the time period before the heating element was placed in the measurement room have been excluded. The interpolated temperature values show a good agreement with the readings from Sensor 26 when the heating element is switched off, while measured temperatures are marginally higher when the element is active. A possible reason for this is that the abrupt increase in temperature in the vicinity of the active heating element requires a higher sensor-density in order to be accurately represented by the model. Note that our dataset does not contain readings from the sensors when the heating element is located in position 3 while inactive. The RMS deviation of the interpolated value from the reading of sensor 26 is provided in table 4. It can be seen that the deviation of the interpolated value from the measured one is higher when the heating element is on compared to when it is off. In particular, the deviation in Position 1 is higher by a factor of 8, while that in Position 2 is higher by a factor of 6.

Propagated uncertainty
The convergence of the mean interpolated temperature as well as the variance of the interpolated temperature at the locations of the 24 sensors determined using Kriging, as a function of the number of Monte Carlo trials, N trials , is illustrated in figure 8 for up to N trials = 10 000. We see that, the empirical mean interpolated temperature determined using Kriging at the locations of sensors converges to the sensor readings while the empirical variance of the interpolated temperatures   converges to the individual sensor uncertainties. The convergence of the mean values is significantly faster (within 1% of the true values after 1000 trials) than that of the variances (within 1% of the true values after 10 000 trials for sensors with a higher uncertainty). The variation of the propagated uncertainty determined using the Monte Carlo approach presented in section 3.3.1 for the 200th datapoint (see table 1) across the volume of the measurement room is illustrated in figure 9. The components used in determining the total standard deviation as proposed in section 3.3.1 (see equations (13)- (20)) have also been provided. At locations with a relatively low sensor density, for instance in the center of the room between the RONDCOM and CONTURA devices, a region of high uncertainty can be observed both for the total standard deviation as well as for the variance of the Kriging mean. Similarly, a region of high uncertainty is also evident at the right side of the room, i.e. at the wall corresponding to x = 7.51 m-a region with no nearby sensors. In contrast, the locations with a high concentration of sensors show have a much lower associated uncertainty. However, the variation of the Kriging mean shows counterintuitively low values (comparable or lesser than the sensor uncertainties) at certain unsampled locations, particularly at those in the vicinity of RONDCOM. In contrast, the mean value of the Kriging variance is close to the sampled points and increases at locations far from the sensors. The combination of the values determined using the Monte Carlo approach, however, show a more physically consistent behavior. In other words, the variation arising due to the inherent uncertainty of the sensor measurements is combined with the spatial variation resulting from the Kriging interpolation to generate a more complete measure of the propagated uncertainty for the  interpolated value. In figure 10, the propagated uncertainty determined using the Monte Carlo approach is plotted along side the mean temperature determined using Kriging for the 1D section of the room containing Sensors 09, 11, 13, 17 and 19 (see figure 1). Similar plots for the mean Kriging variance and the total standard deviation are also provided. From the plots it is apparent that Monte Carlo approach reproduces the sensor uncertainties at their individual locations. However, the uncertainties determined in this way for the interpolated points is lower than those for the measured sensor values. Combining the Monte Carlo uncertainties with the Kriging standard deviation according to the law of total variance (cf equation (20)) results in physically consistent uncertainty values at points far away from the sensors. In doing so, the Kriging variance represents the unexplained component of the total uncertainty, whereas the Monte Carlo propagated uncertainty represents to the explained component.

Conclusions
A method for interpolating room temperatures and computing the associated uncertainty by combining the effects of the sensor measurement uncertainty determined from calibration along with the estimated spatial variation arising from the interpolation model has been presented. In particular, Ordinary Kriging was used to interpolate the temperature readings from a limited number of sensors in a measurement room housing CMMs (see figure 1). Readings from the 24 sparselydistributed sensors corresponding to the measurement setup were taken over a period of four months from June to October 2022. As a result, periods with both relatively warmer and colder days were part of the experiment. Moreover, a moveable trolley with an attached heating element and two additional sensors was placed in the room in three different positions for three distinct periods in order to assess the effects of localized disturbances on the measurement system. The performance of the Kriging interpolation was assessed by comparing the temperature measured by one of the sensors moving with the heating element (Sensor 26) with the interpolated value at its physical location determined using the fixed sensors (Sensors 01-24). It was found that the presence of the heating element negatively affects the interpolation, i.e. the average deviation of the interpolated value from the readings of Sensor 26 increase at least by a factor of 6 when the heating element is switched on. The measurement uncertainty propagated from the sensors to a given point in the volume of the measurement room was quantified using a combination of two values determined using a Monte Carlo approach. The first value corresponds to the variance of the mean interpolated value determined using Kriging, whereas the latter denotes the mean value of the Kriging variance across the Monte Carlo trials (cf section 3.3.1). The two aforementioned values were, respectively, identified with the unexplained and explained components of the total variance of a random variable (cf equation (20)), such that the combination of the two can be interpreted as the uncertainty of the interpolated value. The interpolated temperatures and the corresponding uncertainties depend most strongly on those of the nearest sensor such that the propagated uncertainty is smaller in the vicinity of the sensor with lower uncertainty. The present study dealt with the problem of interpolation and uncertainty propagation for a given instance of sensor readings and the influence from anisotropies arising from the operation of air-conditioning units in the measurement rooms has not been considered. In future research, the effect of cooling as well as of adding and removing sensors to the setup will be investigated. Furthermore, the viability utilizing the temporal correlations between sensor values to improve and, potentially, speed-up the Monte Carlo uncertainty propagation will be examined.

Data availability statement
The data cannot be made publicly available upon publication because they contain commercially sensitive information. The data that support the findings of this study are available upon reasonable request from the authors.