Modeling of a 3D Temperature Field by Integrating a Physics-Speciﬁc Model and Spatiotemporal Stochastic Processes

: Engineering thermal management (ETM) is one of the critical tasks for quality control and system surveillance in many industries, and acquiring the temperature ﬁeld and its evolution is a prerequisite for e ﬃ cient thermal management. By harnessing the sensing data from sensor networks, an unprecedented opportunity has emerged for an accurate estimation of the temperature ﬁeld. However, limited resources of sensor deployment and computation capacity pose a great challenge while modeling the spatiotemporal dynamics of the temperature ﬁeld. This paper presents a novel temperature ﬁeld estimation approach to describe the dynamics of a temperature ﬁeld by combining a physics-speciﬁc model and a spatiotemporal Gaussian process. To reduce the computational burden while dealing with a large set of spatiotemporal data, we employ a tapering covariance function and develop an associated parameter estimation procedure. We introduce a case study of grain storage to show the e ﬀ ectiveness and e ﬃ ciency of the proposed approach.


Introduction
Engineering thermal management (ETM) is one of the critical tasks for quality control and system surveillance in many industries. As a prerequisite in ETM, the temperature field distribution should be known. Generally, estimation of a temperature field is of great interest, as accurately estimating a temperature field is required in a variety of research domains, such as ecology [1,2], health care [3,4], and meteorology [5,6]. An accurate estimation of a temperature field may provide information for the guidance of the optimum design of engineering structures to reduce production cost and conserve energy [7].
Affected by a number of extrinsic factors (e.g., ambient temperature) and intrinsic factors (e.g., nature currents and moisture migration), a temperature field is often regarded as a complex system with stochastic physical dynamics and spatiotemporally varying uncertainties. For example, the temperature field in a granary (Figure 1), in which grain temperature at each location varies over time, is considered as a spatiotemporally varying thermodynamic system. However, due to the unknown inner mechanism and varying boundary conditions, accurately estimating a temperature field with spatiotemporal stochastic dynamics remains a challenging task.
In conventional field estimation methods, a temperature field is estimated by physics-specific models. A heat transfer model for a cubic granary was adopted by Wang and Zhang [8] and Rusinek and Kobylka [9]. Jia, Sun, and Cao [10] and Alagusundaram et al. [11] developed an unsteady heat transfer model to characterize a temperature field of a cylindrical granary. Khatchatourian and Oliveira [12] introduced a mathematical approach to model the airflow and the thermal state in a granary. For the In conventional field estimation methods, a temperature field is estimated by physics-specific models. A heat transfer model for a cubic granary was adopted by Wang and Zhang [8] and Rusinek and Kobylka [9]. Jia, Sun, and Cao [10] and Alagusundaram et al. [11] developed an unsteady heat transfer model to characterize a temperature field of a cylindrical granary. Khatchatourian and Oliveira [12] introduced a mathematical approach to model the airflow and the thermal state in a granary. For the stored grain ecosystem, a 3D transient heat, mass, and momentum transfer model was developed by Gastón et al. [13] and Jian et al. [14] to predict temperatures and moisture contents of stored grains. Jian et al. [15] and Lawrence [16] extended the Gastón's and Jian's work by considering a species transfer model, in which data of ambient temperature and humidity, wind speed, and solar radiation were used. Conventional methods for field estimation are mostly achieved by the simulation of physical models, in which sensing observations collected from distributed sensor networks are not considered. These methods may lead to a discrepancy between model predictions and observations because of stochastic dynamics and uncertainties in the temperature fields. For instance, excessive moisture of grain may lead to local overheating of grain, which cannot be easily captured by conventional field estimation methods.
With the development of sensor technologies, wireless sensor networks have emerged as a key tool to collect sensing data and enable the sensing data to measure the fields [17][18][19]. Many researchers have employed statistical methods for the estimation of temperature fields using sensing observations, which aim to capture stochastic dynamics and uncertainties of temperature fields. For instance, Uchiyama, Hamatani, and Higashino developed a health core temperature estimation approach using data collected from wearable sensors [20]. Yin et al. analyzed temperature fields of the stored grain in different warehouses using grain temperature sensing data [21]. Uncertainty quantification to temperature sensing data is an essential part for temperature field estimation. Some research has focused on data uncertainty quantification problems [22,23] For example, a data uncertainty quantification approach proposed by Lermusiaux is for estimating temperature fields in the ocean [24]. Sapsis and Lermusiaux employed dynamically orthogonal field equations to characterize continuous dynamical systems [25]. Some other research work has focused on tackling uncertainties using Gaussian processes. For example, Rasmussen presented the application of Gaussian processes in machine learning [26]. Krause, Singh, and Guestrin used a Gaussian process to achieve near-optimal sensor placements [27]. Sang and Huang proposed an approximation approach to covariance functions of Gaussian process for spatial data [28]. Cressie, Noel, and Wikle applied Gaussian process models to study the spatiotemporal correlation of a 2D or 3D field [29]. However, only limited sensor observations are collected, because only sparse sensors are deployed in sensor networks given the high costs of sensors and limited budget. Due to limited sensing With the development of sensor technologies, wireless sensor networks have emerged as a key tool to collect sensing data and enable the sensing data to measure the fields [17][18][19]. Many researchers have employed statistical methods for the estimation of temperature fields using sensing observations, which aim to capture stochastic dynamics and uncertainties of temperature fields. For instance, Uchiyama, Hamatani, and Higashino developed a health core temperature estimation approach using data collected from wearable sensors [20]. Yin et al. analyzed temperature fields of the stored grain in different warehouses using grain temperature sensing data [21]. Uncertainty quantification to temperature sensing data is an essential part for temperature field estimation. Some research has focused on data uncertainty quantification problems [22,23] For example, a data uncertainty quantification approach proposed by Lermusiaux is for estimating temperature fields in the ocean [24]. Sapsis and Lermusiaux employed dynamically orthogonal field equations to characterize continuous dynamical systems [25]. Some other research work has focused on tackling uncertainties using Gaussian processes. For example, Rasmussen presented the application of Gaussian processes in machine learning [26]. Krause, Singh, and Guestrin used a Gaussian process to achieve near-optimal sensor placements [27]. Sang and Huang proposed an approximation approach to covariance functions of Gaussian process for spatial data [28]. Cressie, Noel, and Wikle applied Gaussian process models to study the spatiotemporal correlation of a 2D or 3D field [29]. However, only limited sensor observations are collected, because only sparse sensors are deployed in sensor networks given the high costs of sensors and limited budget. Due to limited sensing observations and the existence of various data uncertainties, temperature field estimation purely using statistical approaches fails to achieve an accurate estimation result. Therefore, an accurate temperature field estimation requires combining appropriate physical laws and statistical approaches to finely characterize temperature distribution and associated evolution. However, few studies have addressed this issue to the best of our knowledge.
In this paper, we develop a field estimation approach to model spatiotemporal dynamics of temperature fields based on grain storage sensor networks by combining a physics-specific model and spatiotemporal stochastic processes. First, a three-dimensional thermodynamic model, that is, the 3D heat transmission model for the grain temperature field, is employed to describe the global temperature profile under ideal conditions by considering extrinsic factors. Then, we capture local temperature changes caused by intrinsic factors by a spatiotemporal Gaussian process based on the deviance of sensing observations from the global temperature profile. Directly estimating the spatiotemporal dynamics with conventional Gaussian process models will generally be computationally infeasible, as O n 3 operations is required for a dataset of size n. Hence, we borrow a tapering covariance function and integrate it into the Gaussian process model to reduce the computational burden. We conduct a real case of temperature fields in grain storage to validate the model performance of our field estimation method.
The remainder of this article is organized in the following. Section 2 introduces the methodology of the temperature field estimation. Section 3 discusses a real case study in grain storage to validate the proposed method. Section 4 provides the discussion, the conclusion, and the future work.

Formulation
We denote y(s, t) as the response variable of a temperature field at location s and time t, where s ∈ S ⊆ R 3 , t ∈ T ⊆ R + . We describe the temperature field by a mixed effect model framework, which contains three parts: the global temperature profile, the local variability, and the model error.
where g(s, t) denotes the mean function describing the global temperature profile, δ(s, t) denotes the local variability, and ε denotes the model error. As shown in Figure 2, we use a physics specific model, that is, the thermodynamic model, to characterize the global temperature profile under ideal conditions by considering extrinsic factors. There exists a discrepancy between the global temperature profile acquired from the conventional thermodynamic model and real temperature values due to some intrinsic factors, such as the existence of nature currents and moisture migration. To deal with this issue, we model this discrepancy as the local variability. The local variability is affected by these intrinsic factors, which cannot be fully captured by the thermodynamic model. In our method, the sensing data provides us a unique chance to capture the uncertainties of the local variability caused by intrinsic factors. We adopt a spatiotemporal Gaussian process model to describe the local variability using the deviance of the sensing data from the global temperature profile. Therefore, both the effects from extrinsic and intrinsic factors, e.g., heat transfer, nature currents and moisture migration, are substantially incorporated into the mixed-effect model framework.
Appl. Sci. 2019, 9, x 3 of 13 and spatiotemporal stochastic processes. First, a three-dimensional thermodynamic model, that is, the 3D heat transmission model for the grain temperature field, is employed to describe the global temperature profile under ideal conditions by considering extrinsic factors. Then, we capture local temperature changes caused by intrinsic factors by a spatiotemporal Gaussian process based on the deviance of sensing observations from the global temperature profile. Directly estimating the spatiotemporal dynamics with conventional Gaussian process models will generally be computationally infeasible, as ( ) operations is required for a dataset of size . Hence, we borrow a tapering covariance function and integrate it into the Gaussian process model to reduce the computational burden. We conduct a real case of temperature fields in grain storage to validate the model performance of our field estimation method. The remainder of this article is organized in the following. Section 2 introduces the methodology of the temperature field estimation. Section 3 discusses a real case study in grain storage to validate the proposed method. Section 4 provides the discussion, the conclusion, and the future work.

Formulation
We denote ( , ) as the response variable of a temperature field at location and time , where ∈ ⊆ , ∈ ⊆ . We describe the temperature field by a mixed effect model framework, which contains three parts: the global temperature profile, the local variability, and the model error.
where ( , ) denotes the mean function describing the global temperature profile, ( , ) denotes the local variability, and ε denotes the model error. As shown in Figure 2, we use a physics specific model, that is, the thermodynamic model, to characterize the global temperature profile under ideal conditions by considering extrinsic factors. There exists a discrepancy between the global temperature profile acquired from the conventional thermodynamic model and real temperature values due to some intrinsic factors, such as the existence of nature currents and moisture migration.
To deal with this issue, we model this discrepancy as the local variability. The local variability is affected by these intrinsic factors, which cannot be fully captured by the thermodynamic model. In our method, the sensing data provides us a unique chance to capture the uncertainties of the local variability caused by intrinsic factors. We adopt a spatiotemporal Gaussian process model to describe the local variability using the deviance of the sensing data from the global temperature profile. Therefore, both the effects from extrinsic and intrinsic factors, e.g., heat transfer, nature currents and moisture migration, are substantially incorporated into the mixed-effect model framework.

Physics-Specific Model for the Global Temperature Profile
We employ a physics-specific model to describe the global temperature profile of a temperature field over time. To illustrate how to choose a physics-specific model for a specific temperature field, we take the grain temperature field as an example.
We consider a temperature field in a cubic granary. To describe the global temperature profile in a cubic granary, location s in Equation (1) is represented by x, y, z in Cartesian coordinates. We employ a 3D heat transmission in a granary as the mean function as follows: where the parameters α x , α y , and α z are determined on the basis of thermal properties. For stored grains, α x , α y , and α z denote the thermal diffusivities of the grains in the x, y, and z directions, respectively. Alternatively, if the thermal diffusivities of the grains are unknown, we can calculate α x , α y , and α z using the proposed parameter estimation approach in [30]. The finite-difference method is employed to solve Equation (2). We use an explicit marching method to calculate the global temperature profile g(s, t), as shown in Equation (3), where a central difference in x, y, and z directions, and a forward difference in time t are considered. To be specific, where i x , i y , i z , and i t are the indexes of locations x, y, z and the time t. ∆x, ∆y, ∆z, and ∆t are the grid spacings in the x, y, and z directions of the space domain and time t.
To satisfy a highly accurate numerical solution, this method is stable under the condition [31]: Given the grain temperature at the initial time point and the temperature at the walls of the granary as the boundary conditions, the global temperature profiles are obtained. We adopt a Gaussian process model [26] for interpolation using the sensing data on the granary walls to obtain the boundary conditions.

Formulation of the Gaussian Process Model
Local variability of a temperature field can be affected by many intrinsic factors that cannot be fully captured by physics-specific models. For example, excessive moisture of grain leads to a quick local temperature increment, which cannot be directly characterized by conventional physical models. As shown in Figure 3, although the dynamics model is with continuous variables that consider the spatiotemporal correlation, the local variability in spatial and temporal domain caused by intrinsic factors substantially exists and exhibits with various correlation patterns. Given that the thermodynamic model is able to capture nonstationary profiles of a temperature field across the spatial and temporal domains, the local variability is assumed to be stationary. Without loss of generality, we adopt a spatiotemporal Gaussian process model to describe the spatiotemporal dependence structures of local variability using the deviance of the sensing observations from the global temperature profile in the following: where Σ δ is the covariance matrix of the Gaussian process.
where Σ is the covariance matrix of the Gaussian process. The covariance matrix of the Gaussian Process is obtained from a covariance function describing the spatial and temporal correlations. Because various thermal transmission effects exist in different directions of the granary, there are various spatial correlation patterns of local variability in different directions of the granary. Temporal evolution also exists among the time domain. Therefore, we use an anisotropic spatiotemporal covariance function to characterize spatiotemporal correlation between two different spatiotemporal locations as follows: where = { , , , } is the distance metric of two locations in the spatial and temporal domains, respectively, and = { , , , , } is the set of parameters in the covariance function.
In the spatial domain, , , and are the parameters corresponding to the distance in , , and directions respectively, which capture anisotropic spatial correlation structures of a granary; in the temporal domain, a decay parameter corresponding to the time dependence effect is used to capture the temporal correlation of the local variability.

Covariance Tapering
Traditional Gaussian process models, which require the computational complexity of ( ) to estimate parameters with a dataset of size n [32]. In order to reduce the computational burden of Gaussian process models and satisfy the necessity in practice, we adopt a covariance tapering method to characterize the covariance function. Considered that distant pairs of temperature values are uncorrelated, the covariance tapering method aims to model a compactly supported covariance function in which the correlation of distant pairs is defined as zero. The covariance matrix of Gaussian process models obtained by the covariance tapering method is positive definite and can retain most information of the temperature field.
We adopt a tapering covariance function in the Gaussian process model by taking the product of an original covariance function ( ; ) and a tapering function ( ; ) as follows: where ( ; ) is defined as an isotropic correlation function in the spatial domain that is zero where is a range parameter. We use the spherical tapering function as the tapering function in our The covariance matrix of the Gaussian Process is obtained from a covariance function describing the spatial and temporal correlations. Because various thermal transmission effects exist in different directions of the granary, there are various spatial correlation patterns of local variability in different directions of the granary. Temporal evolution also exists among the time domain. Therefore, we use an anisotropic spatiotemporal covariance function to characterize spatiotemporal correlation between two different spatiotemporal locations as follows: where h = δx, δy, δz, δt is the distance metric of two locations in the spatial and temporal domains, respectively, and θ = σ 2 , a x , a y , a z , ϕ is the set of parameters in the covariance function. In the spatial domain, a x , a y , and a z are the parameters corresponding to the distance in x, y, and z directions respectively, which capture anisotropic spatial correlation structures of a granary; in the temporal domain, a decay parameter ϕ corresponding to the time dependence effect is used to capture the temporal correlation of the local variability.

Covariance Tapering
Traditional Gaussian process models, which require the computational complexity of O n 3 to estimate parameters with a dataset of size n [32]. In order to reduce the computational burden of Gaussian process models and satisfy the necessity in practice, we adopt a covariance tapering method to characterize the covariance function. Considered that distant pairs of temperature values are uncorrelated, the covariance tapering method aims to model a compactly supported covariance function in which the correlation of distant pairs is defined as zero. The covariance matrix of Gaussian process models obtained by the covariance tapering method is positive definite and can retain most information of the temperature field.
We adopt a tapering covariance function in the Gaussian process model by taking the product of an original covariance function C(h; θ) and a tapering function c taper (δs; γ) as follows: Appl. Sci. 2019, 9, 2108 6 of 13 where c taper (δs; γ) is defined as an isotropic correlation function in the spatial domain that is zero whenever where γ is a range parameter. We use the spherical tapering function as the tapering function in our model, which is defined by otherwise.

Parameter Estimation for θ
We incorporate the tapering covariance function Equation (8) into the Gaussian process Equation (5). Then, the tapering covariance matrix is where Σ(θ) and A(γ) are the covariance matrices obtained by C(h; θ) and c taper (δs; γ), respectively; the notation • denotes the Schur product that is defined as the element wise matrix product. For example, two m × n matrixes, P and Q, have the Schur product P • Q = p ij q ij . We use maximum-likelihood estimation (MLE) approach to estimate θ. Specifically, the log-likelihood function of θ is as follows where Z = Y − G denotes the difference between the sensing observations and the global temperature profiles, and d is a constant term that is uncorrelated with θ. Maximizing l(θ) is to solve an unbiased estimating equation for θ, i.e., E ∂ ∂θ l(θ) = 0. The derivative of l(θ) is equal to Equation (11) can be used in conjunction with a numerical optimization algorithm, such as conjugate gradients, to obtain the estimation of θ [26].

Selection of Range Parameter γ
The choice of the range parameter γ in the tapering covariance function influences the accuracy of the field estimation and the computational complexity of our proposed model. Given that maximizing the likelihood function in Equation (10) is to solve an unbiased estimating equation for θ, it is suggested that an estimator of the sampling variability ofθ should be estimated using the robust information criterion in [33]. To achieve a trade-off between field estimation accuracy and computational complexity, we adopt the sampling variability (SV) and the model computational time (CT) as a criterion to select an appropriate γ.
Let U(s, t; θ) be an unbiased estimating function for θ. The robust information criterion for U is U is the matrix of derivatives of U with respect to θ, The sampling variability is defined by the diagonal elements of E(U) −1 , namely, the trace of E(U) −1 .
The procedure of selecting γ is as follows (Figure 4): First, we obtain a pilot estimatê θ pilot . For example, we calculate the MLE for a small and randomly chosen subset of sensing observations. Second, we calculate the sampling variability SV(θ pilot , γ i ) = tr{E(G(s, t;θ pilot , γ i )) −1 } and model computational time CT(θ pilot , γ i ) for an available value of γ i . The initial value γ 1 is the minimal range parameter value. Third, we repeat the second step and calculate the sampling variability SV(θ pilot , γ i ) and model computational time CT(θ pilot , γ i ) for γ i with i = 2, 3, . . . , n γ , where n γ is the number of γ; Finally, we select γ * to construct a balance between the sampling variability SV(γ) = {SV(θ pilot , γ i ), i = 1, 2, . . . , n γ } and the model computational time CT(γ) = {CT(θ pilot , γ i ), i = 1, 2, . . . , n γ }, that is, smaller sampling variability and shorter model computational time.

Spatiotemporal Field Estimation
When the covariance parameters are estimated and the range parameter is selected, we obtain the spatiotemporal temperature field by given environmental conditions and sensing observations. Temperature * at a new location * and time * is estimated by * = * + Σ * Σ ( − ), (15) and the variance of * is calculated by where * is the global temperature profile, is the vector of sensing observations, is the vector of global temperature profiles for , Σ * * is the variance matrix of * , Σ * is the covariance matrix between and * , and Σ is the covariance matrix of .

Case Study
We conducted a real case study of a grain temperature field to validate our proposed method. The sensing data in our real case were collected from the grain temperature sensor networks in a national granary in central China at least every seven days from 31 January in one year to 4 March in the second year. Figure 5 presents the grain temperature sensor network of the granary. The sensor

Spatiotemporal Field Estimation
When the covariance parameters θ are estimated and the range parameter γ is selected, we obtain the spatiotemporal temperature field by given environmental conditions and sensing observations. Temperature Y * at a new location s * and time t * is estimated bŷ and the variance ofŶ * is calculated by where g * is the global temperature profile, Y is the vector of sensing observations, G is the vector of global temperature profiles for Y, Σ * * δ is the variance matrix of Y * , Σ * δ is the covariance matrix between Y and Y * , and Σ δ is the covariance matrix of Y.

Case Study
We conducted a real case study of a grain temperature field to validate our proposed method. The sensing data in our real case were collected from the grain temperature sensor networks in a national granary in central China at least every seven days from 31 January in one year to 4 March in the second year. Figure 5 presents the grain temperature sensor network of the granary. The sensor network consisted of 240 sensors, which were uniformly distributed in a traditional granary with the volume of 46 m × 26 m × 6 m. The sensors were deployed at the initial positions of 0.5 m, 0.5 m, and 0.3 m in the space domain and located at 5 m intervals in the x and y directions and 1.8 m intervals in the z direction. There were in total 17,280 sensing observations. We randomly selected training data from sensing data in the first 66 time points from 31 January in one year to 1 February in the second year, and set the remaining data as test data. The ratio of training data set to test data set was set as 3:1.
To validate the prediction capability of our proposed model, we also retained the sensing data from 2 February to 4 March in the second year as the test data.

Material Thermal Diffusivity in , , and Directions
Mixed wheat 5.56, 5.55, 1.50 (× 10 m /s) The physical properties of stored grains were given in Table 1, and the attribute values were further incorporated into the thermodynamic model. In Table 1, thermal diffusivity in the vertical directions is greater than the transverse directions due to the faster airflow enhancing the heat transfer in vertical direction.
The validation of the proposed field estimation approach contained three layers: modeling of the global temperature profile, parameter estimation of the proposed Gaussian process model and spatiotemporal field estimation.

Global Temperature Profile
For the global temperature profile, the parameters in the mean function are determined by the properties of grain that are listed in Table 1. Given the grain temperature sensing data in the real case, we chose ∆ = 0.5 m, ∆ = 0.5 m, ∆ = 0.5 m, and ∆ = 1 day to satisfy the stability condition of the finite-difference method in Equation (4). We acquired the sensing data of the granary walls collected by the sensors on the walls and adopted a Gaussian process model [26] for interpolation using the sensing data on the granary walls to obtain the boundary conditions of the thermodynamic model.

Parameter Estimation of the Proposed Gaussian Process Model
For the parameter estimation of the tapering covariance function in the Gaussian process, we have initially addressed the selection of the range parameter . We obtained a pilot estimate of the covariance parameters using a subset training data of size 540, which was randomly chosen from the training data set. Figure 6 shows the curves of both the sampling variability and the computational time given a series of values of . For the increasing values of , the sampling variability decreases with a substantial range when is small and tends to become flat when is The physical properties of stored grains were given in Table 1, and the attribute values were further incorporated into the thermodynamic model. In Table 1, thermal diffusivity in the vertical directions is greater than the transverse directions due to the faster airflow enhancing the heat transfer in vertical direction. The validation of the proposed field estimation approach contained three layers: modeling of the global temperature profile, parameter estimation of the proposed Gaussian process model and spatiotemporal field estimation.

Global Temperature Profile
For the global temperature profile, the parameters in the mean function are determined by the properties of grain that are listed in Table 1. Given the grain temperature sensing data in the real case, we chose ∆x = 0.5 m, ∆y = 0.5 m, ∆z = 0.5 m, and ∆t = 1 day to satisfy the stability condition of the finite-difference method in Equation (4). We acquired the sensing data of the granary walls collected by the sensors on the walls and adopted a Gaussian process model [26] for interpolation using the sensing data on the granary walls to obtain the boundary conditions of the thermodynamic model.

Parameter Estimation of the Proposed Gaussian Process Model
For the parameter estimation of the tapering covariance function in the Gaussian process, we have initially addressed the selection of the range parameter γ. We obtained a pilot estimate of the covariance parameters θ using a subset training data of size 540, which was randomly chosen from the training data set. Figure 6 shows the curves of both the sampling variability and the computational time given a series of values of γ. For the increasing values of γ, the sampling variability decreases with a substantial range when γ is small and tends to become flat when γ is large. In contrast, the computational time of field estimation rises rapidly as with the increase of γ. We chose γ * = 22 as the appropriate range parameter in the spatial domain by weighting the trade-off between sampling variability and computational time. Then, we used the MLE method to acquire estimated covariance parameters θ * , which are shown in Table 2.

Spatiotemporal Field Estimation
For the spatiotemporal field estimation, we used the mean square error (MSE) as in Equation (17) to evaluate the proposed field estimation method.
where ( , ) denotes a sensing observation in the test data set, and ( , ) denotes the estimated temperature value of ( , ) . We tested the proposed field estimation method by estimating the grain temperature at the sensor sites of test data set by Equations (15) and (16). The MSE is equal to 0.04, which indicates the effectiveness of the estimation results using the proposed field estimation method. We show the examples of temperature profiles by the proposed field estimation method in Figure 7. The blue-star dots represent the training data, the red-plus dots represent the test data, and the blue lines denotes the estimated temperature profiles using our proposed method. The results show both interpolation and prediction capability of the proposed model with acceptable confidence intervals.
To show the superiority of our approach, which integrates the physics-specific model and the spatiotemporal Gaussian process model, we compared it with other existing methods. In Table 3, we compared the model performance by changing mean functions proposed in the relevant literature, that is, a constant and a sinusoidal model. The result shows the superiority of the thermodynamic model as the mean function. In Table 4, our proposed method delivers a better result than the thermodynamic model and consumes less than 50% of the computational time than the conventional Gaussian process model.

Spatiotemporal Field Estimation
For the spatiotemporal field estimation, we used the mean square error (MSE) as in Equation (17) to evaluate the proposed field estimation method.
where Y(s i , t m ) denotes a sensing observation in the test data set, andŶ(s i , t m ) denotes the estimated temperature value of Y(s i , t m ). We tested the proposed field estimation method by estimating the grain temperature at the sensor sites of test data set by Equations (15) and (16). The MSE is equal to 0.04, which indicates the effectiveness of the estimation results using the proposed field estimation method. We show the examples of temperature profiles by the proposed field estimation method in Figure 7.
The blue-star dots represent the training data, the red-plus dots represent the test data, and the blue lines denotes the estimated temperature profiles using our proposed method. The results show both interpolation and prediction capability of the proposed model with acceptable confidence intervals.
To show the superiority of our approach, which integrates the physics-specific model and the spatiotemporal Gaussian process model, we compared it with other existing methods. In Table 3, we compared the model performance by changing mean functions proposed in the relevant literature, that is, a constant and a sinusoidal model. The result shows the superiority of the thermodynamic model as the mean function. In Table 4, our proposed method delivers a better result than the thermodynamic model and consumes less than 50% of the computational time than the conventional Gaussian process model.   Figure 8 shows an example of temperature profiles at different locations estimated by the proposed field estimation method. The grain temperature field obtained by our model changes with seasonal effects. Specifically, grain temperature inside the granary varies more slowly over time and has a smaller fluctuating range than grain temperature at the outer locations. This finding is

Model MSE
Constant + Gaussian process with tapering covariance function 7.86 Sinusoidal model + Gaussian process with tapering covariance function 0.60 Thermodynamic model + Gaussian process with tapering covariance function 0.04  Figure 8 shows an example of temperature profiles at different locations estimated by the proposed field estimation method. The grain temperature field obtained by our model changes with seasonal effects. Specifically, grain temperature inside the granary varies more slowly over time and has a smaller fluctuating range than grain temperature at the outer locations. This finding is consistent with the grain temperature fields in most of the granaries in China. Moreover, it provides a new quantitative understanding of temperature field evaluation in granaries [34]. consistent with the grain temperature fields in most of the granaries in China. Moreover, it provides a new quantitative understanding of temperature field evaluation in granaries [34].

Conclusions
This study proposes an estimation method for a 3D spatiotemporal temperature field in grain storage sensor networks. In conventional field estimation methods, a temperature field could not be estimated with high accuracy by either the pure physical or the statistical model because of the sparsity and spatiotemporal uncertainties of temperature observations. To tackle this problem, we obtain the estimation of a temperature field by combining the physics-specific model and the spatiotemporal stochastic process. First, we employ a physics-specific model to characterize the global temperature profile. Then, we adopt a spatiotemporal Gaussian process to describe the local variability using the deviance of the sensing observations from the global temperature profile. To reduce the computational burden of conventional Gaussian process models, we employ a tapering covariance function by integrating a tapering function with the original covariance function of Gaussian process model. In the case study, we validate the performance of our method by using grain temperature sensing observations collected from grain storage sensor networks. The results of the real case study deepen our understanding of the temperature field in grain storage and provide useful information for the guidance of future sensor placement designs and configurations.
In the future work, we aim to adaptively calibrate the parameters in the spatiotemporal field estimation model for precise engineering thermal management.
Author Contributions: X.Z. conceived the study and contributed to the design of the study, modeling, data interpretation, and revised the manuscript. D.W. contributed to the development of algorithms, data evaluation, data analysis, and drafted the manuscript. Both authors reviewed the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.

Conclusions
This study proposes an estimation method for a 3D spatiotemporal temperature field in grain storage sensor networks. In conventional field estimation methods, a temperature field could not be estimated with high accuracy by either the pure physical or the statistical model because of the sparsity and spatiotemporal uncertainties of temperature observations. To tackle this problem, we obtain the estimation of a temperature field by combining the physics-specific model and the spatiotemporal stochastic process. First, we employ a physics-specific model to characterize the global temperature profile. Then, we adopt a spatiotemporal Gaussian process to describe the local variability using the deviance of the sensing observations from the global temperature profile. To reduce the computational burden of conventional Gaussian process models, we employ a tapering covariance function by integrating a tapering function with the original covariance function of Gaussian process model. In the case study, we validate the performance of our method by using grain temperature sensing observations collected from grain storage sensor networks. The results of the real case study deepen our understanding of the temperature field in grain storage and provide useful information for the guidance of future sensor placement designs and configurations.
In the future work, we aim to adaptively calibrate the parameters in the spatiotemporal field estimation model for precise engineering thermal management.
Author Contributions: X.Z. conceived the study and contributed to the design of the study, modeling, data interpretation, and revised the manuscript. D.W. contributed to the development of algorithms, data evaluation, data analysis, and drafted the manuscript. Both authors reviewed the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.