Implementation of a Parallel GPU-Based Space-Time Kriging Framework

: In the study of spatiotemporal geographical phenomena, the space–time interpolation method is widely applied, and the demands for computing speed and accuracy are increasing. For nonprofessional modelers, utilizing the space–time interpolation method quickly is a challenge. To solve this problem, the classical ordinary kriging algorithm was selected and expanded to a spatiotemporal kriging algorithm. Using the OpenCL framework to integrate central processing unit (CPU) and graphic processing unit (GPU) computing resources, a parallel spatiotemporal kriging algorithm was implemented, and three experiments were conducted in this work to verify the results. The results indicated the following: (1) when the size of the prediction point dataset is consistent, the performance of the method is robust with the increasing size of the observation point dataset; (2) the acceleration effect of the parallel method increases with an increased number of predicted points. Compared with the original sequential program, the implementation of the improved parallel framework showed a 3.23 speedup, which obviously shortens the interpolation time; (3) when cross-validating the temperature data in the Beijing Tianjin Hebei region, the space–time acceleration model provides a better ﬁt than traditional pure space interpolation.


Introduction
A variety of spatial interpolation methods have been proposed to allow better understanding of the spatial distribution of an area covered by existing observed points.Basically, spatial interpolation becomes a standard tool for geographic information system (GIS) software and plays an important role in Digital Earth (DE) research [1].Due to the continuous accumulation of historical data and the societal need for the analysis of variables that vary in space and time, such as weather and air quality variables, the collection and processing of spatiotemporal data is rapidly increasing [2].Researchers have begun to pay attention to the space-time interpolation problem to estimate the spatial and temporal variation characteristics of various spatial elements after time dimension integration and to provide support for decision-making.In recent years, the theoretical aspects of spatiotemporal geostatistics have made great progress [3], but there are still some application problems.Spatiotemporal models are more complex than pure spatial models, and the operations involved in spatiotemporal models are much more expensive.
For these reasons, many scholars have studied the space-time interpolation problem.With the help of existing methods, some scholars have tried to improve time and space interpolations.For example, the kriging method is used for the best, unbiased, optimal estimation of regionalized sampling point variables using the structural features of observations with known spatial distributions.Although originally developed for geostatistics, the kriging method is now widely used in geology [4], hydrology [5], meteorology [6], environmental science [7] and other disciplines.However, in the real world, many geographical phenomena have temporal and spatial evolution processes.As a spatial interpolation method, kriging often ignores the important information about the time dimension in the problem, which is not conducive to further improving the accuracy of the interpolation.
Through the evolution of regionalized variables into spatiotemporal regionalized variables, scholars established spatiotemporal variogram models [8], which can achieve the spatial and temporal extension of kriging interpolation.To date, spatiotemporal kriging interpolation has been applied in many fields.Raja et al. [9] used spatiotemporal kriging interpolation to research the trends of spatiotemporal variation in regional precipitation, which can be used to plan and manage water resources in areas that are heavily dependent on precipitation.Yang et al. [10,11] assessed the sources and spatiotemporal trends of heavy metal accumulation in soil using the spatiotemporal kriging method and provided suggestions for the prevention and control of heavy metal pollution in soil.Jost et al. [12,13] used an interpolation method to determine the distribution of a forest ecosystem and the water storage layer near the desert; their method can assist groundwater managers to make correct decisions.Park et al. [14,15] used the spatiotemporal kriging method to more accurately estimate the spread of air pollutants and the distribution of infectious diseases, adding to omissions in the data collection process.
With the expansion of interpolation algorithms in the time dimension, it is necessary to acquire sample points at different times and to obtain the time trend of the study area.As a result, the number of sampling points increases linearly and interpolation algorithms require large computation and storage resources; the algorithms are slow and inefficient on traditional personal computers (PCs), which is not conducive to the use of applications.
With the development of computing technology, high-performance computing (HPC) technology has been gradually integrated into the field of geographic information, which has improved the application, promotion and development of the field.Multiple processors, parallel clustering, grid computing, graphic processing unit (GPU) incorporation and other high-performance computing technologies have attracted much research attention in the geosciences field [16][17][18].In recent years, parallel computing has mainly appeared in the form of multi-core processors [19].A GPU has a large number of cores, which makes up for the shortcomings of the traditional central processing unit (CPU) architecture.
High-performance GIS computing algorithms are a research hotspot that is mainly divided into two aspects.One aspect is the optimization of processing from the technical level, parallel to existing high-density computing, to improve operational efficiency.The second aspect is the exploration of a new spatiotemporal analysis model using modeling and algorithms to develop a more efficient and convenient spatiotemporal analysis model.
To represent the spatial interpolation model by the kriging interpolation method in the technical level of optimization and new space interpolation algorithm improvement, scholars have begun to use GPUs to reduce computing time and introduce mature parallelization schemes.Cheng [20] accelerated the universal kriging algorithm on the NVIDIA Compute Unified Device Architecture (CUDA) platform and achieved a nearly 18-fold speed increase with respect to the sequential program.Liu [21] used the computation power of modern programmable graphics hardware (GPU) for 3D visualization in a reservoir modeling system.In terms of algorithms and model improvements, Liu [22] proposed an algorithm based on the k-d tree method to address the unevenly distributed spatial data.Hu et al. [23] proposed an fast Fourier transform (FFT)-based parallel algorithm to accelerate regression kriging interpolation, which was computed on a GPU device.As far as we know, spatiotemporal kriging is an improved method for spatiotemporal interpolation; however, the model and the algorithm have changed, and the original parallel method cannot adapt.
Although many people pay attention to the general spatial analysis of GIS, the support of spatiotemporal interpolation and efficiency improvement calculations lack in-depth research.To solve the above problems, the spatiotemporal kriging algorithm was analyzed, data associations were decoupled and parallelized, and a parallel algorithm of spatiotemporal kriging interpolation was proposed.The parallel model of the space-time interpolation is tested with meteorological data in this paper.The purpose of this paper is to propose a new space-time interpolation algorithm and technology supported by HPC, which will provide new research ideas for GIS spatiotemporal analysis modeling and the enrichment of GIS research contents.
This article is arranged as follows: in the first and second parts of Section 2, we introduce the spatiotemporal kriging algorithm and the space-time extension in detail.The third part of Section 2 focuses on the design of the parallel spatiotemporal kriging algorithm and implementation by the OpenCL framework.In Section 3, the experimental results and analysis are given.Section 4 presents the conclusion.

Methodology and Implementation
Similar to spatial interpolation, spatiotemporal geostatistics need are dependent on a suitable spatiotemporal covariance model [24].Many kinds of models have been designed, including linear, separable and nonseparable models [25][26][27].We aimed to implement space-time interpolation using a product-sum model with the support of HPC.Other types of interpolation methods can be properly extended or adjusted on the basis of this method.The data selected in this article are meteorological data.Ordinary kriging interpolation can meet the needs of homogeneous regions.Therefore, ordinary kriging interpolation was used in the experiment.

Ordinary Kriging
The kriging method makes an unbiased, optimal estimation of the regionalized variables of the sampling points and uses a semi-variogram to represent the structural features of the region.Kriging is a linear interpolation method.Using the obtained spatial structure and the known sampling points to predict results, we can use Equation (1): where Ẑs 0 is the estimated value, Z(s i ) refers to each known point, N is the measured value number and λ i is the weight value at the position.The ordinary kriging assumption is that the estimator expected value is equal to the true value, and the prediction error variance is minimized.Therefore, to ensure the best linear unbiased prediction of the prediction point, the sum of the weight coefficients of the surrounding sampling points must be equal to 1.This process can be represented by kriging Equation (2): where γ 11 stands for the variogram value between point i and point j, and µ is a Lagrange multiplier.

Spatiotemporal Kriging Algorithm
The spatiotemporal kriging interpolation method is used to extend the time dimension to the kriging interpolation.The regionalized variables are evolved into spatiotemporal regionalized variables, and the spatiotemporal variogram model is used to fit the sampling points.The relationship between the variogram and covariance is: where γ s (h s ), γ t (h t ), and γ st (h s , h t ) are the spatial, temporal, and spatiotemporal variograms, C s (h s ), C t (h t ), and C st (h s , h t ) are the spatial, temporal, and spatiotemporal covariance functions, and h s and h t are space and time distances, respectively.Due to the different dimensions of the spatial domain and the time domain, the spatial and temporal variograms are difficult to calculate directly.According to a type of product-sum model [28], the sum-product model of space-time covariance function is defined as Equation ( 6): The introduction of k 1 , k 2 and k 3 ensures the positivity of the space-time covariance function, C st (h s , h t ).Where h s = 0, h t = 0 can be obtained from Equation ( 7): C st (0, 0), C s (0), and C t (0) are the nugget values of space, time and space-time, respectively.Using the above equations, the spatiotemporal variogram is obtained (Equation ( 8)): According to ordinary kriging prediction formula, the ordinary spatiotemporal kriging interpolation is obtained with Equation ( 9): The number of spatial sampling points is i = 1, 2... n, the time sampling points are j = 1, 2, ... m, and the total number of samples is n × m.Z*(s 0 , t 0 ) is the linear estimation of the interpolation point (s 0 , t 0 ).Z s i , t j is the value of the neighboring sampling point s i , t j .λ (i,j) is the weighting coefficient of the point and the kriging model in Equation (10): M = γ (01,01) , γ (02,01) , . . ., γ (0n,01) , γ (01,02) , . . ., γ (0n,0m) , 1 Specifically, the serial algorithm of the ordinary spatiotemporal kriging interpolation consists of the following four steps:

•
Step 1: Calculate the variograms of the pure time domain and the pure space domain and fit them.At different time points, the space point has a different nugget, sill and range.Calculate the average C st (0, 0), C s (0), and C t (0) here.

•
Step 4: Bring the known points and the spatiotemporal variogram γ st (h s , h t ) into the kriging equations to obtain the weight vector, and then calculate the estimation of an unobserved point.

Design and Framework of the Parallel Algorithm
The spatiotemporal kriging interpolation predicts the interpolation points using the structural features of spatiotemporal sample points.According to the principle of the ordinary spatiotemporal kriging algorithm, the algorithm can be divided into the following four steps:

•
Step 1: Organize the original data and establish the index to facilitate the query retrieval (Spatiotemporal Data Organization Module, STDOM).

•
Step 2: Fit the spatiotemporal variation model of the sample data to obtain the variogram, γ st , of the spatiotemporal structure of the representative sample (Spatiotemporal Variation Model Fitting Module, STVMFM).

•
Step 3: Proceed with the set of points to be predicted, search for the adjacent spatiotemporal points and set up the data sets (Adjacent Spatiotemporal Points Searching Module).

•
Step 4: Use the spatiotemporal variogram, γ st , and the data sets of the neighboring points for the ordinary spatiotemporal kriging interpolation to calculate the estimated values of the unknown points (Ordinary Spatiotemporal Kriging Interpolation Module, OSTKIM), as shown in Figure 1.

Spatiotemporal Data Organization Module (STDOM)
Spatiotemporal data have space coordinates and time information.In the spatiotemporal kriging interpolation algorithm, it is necessary to search for the points neighboring an unknown point in the spatiotemporal sample data.To facilitate this process, an STDOM was designed to index the data.The spatiotemporal index numbered all space points, and then the

Spatiotemporal Data Organization Module (STDOM)
Spatiotemporal data have space coordinates and time information.In the spatiotemporal kriging interpolation algorithm, it is necessary to search for the points neighboring an unknown point in the spatiotemporal sample data.To facilitate this process, an STDOM was designed to index the data.The spatiotemporal index numbered all space points, and then the station number and date were cross-referenced to obtain the index values shown in Table 1.STVMFM is designed to fit different time and space variogram models to obtain the average C s , C t as the γ s , γ t parameters and to calculate k 1 , k 2 , and k 3 to construct the spatiotemporal variogram γ st .To fit the model, according to the spatiotemporal index, data must be divided into the vector matrix of the same time and the vector matrix of the same space (Figure 2).Each vector in the vector matrix is a set of sample points in time or space.Since the solution to the spatiotemporal parameters of each vector is independent, this step can be parallelized.

Adjacent Spatiotemporal Points Searching Module (ASTPSM)
The amount of data for space-time interpolation is much greater than that for the pure spatial interpolation.Therefore, instead of using all the sample points, the points in a certain spatiotemporal range are calculated according to the arrangement obtained by fitting of the variation model.In ASTPSM, using the spatiotemporal index established before, a neighbor points search method based on the KD-tree is used.The KD-tree algorithm is used to find the adjacent station number in the neighborhood and this is crossed over with the neighboring date to obtain the target index.

Ordinary Spatiotemporal Kriging Interpolation Module (OSTKIM)
The OSTKIM module was used to calculate Equation (9).Parallel calculation of the matrix can greatly speed up the interpolation process.In this step, the coefficient matrix [] needs to be inversed, and the weight vector [] is calculated by matrix-vector multiplication.Much progress has been made on the parallelization of matrix operations.Here, an efficient solution that decomposes the coefficient matrix into a lower triangular matrix and an upper triangular matrix by LU decomposition has been chosen before solving a linear system.Then, the system of linear equation decomposition is solved directly by forward and backward substitution without directly calculating the inverse matrix.

OpenCL-Based Implementation
Open Computing Language (OpenCL) is a framework for writing programs on

Adjacent Spatiotemporal Points Searching Module (ASTPSM)
The amount of data for space-time interpolation is much greater than that for the pure spatial interpolation.Therefore, instead of using all the sample points, the points in a certain spatiotemporal range are calculated according to the arrangement obtained by fitting of the variation model.In ASTPSM, using the spatiotemporal index established before, a neighbor points search method based on the KD-tree is used.The KD-tree algorithm is used to find the adjacent station number in the neighborhood and this is crossed over with the neighboring date to obtain the target index.

Ordinary Spatiotemporal Kriging Interpolation Module (OSTKIM)
The OSTKIM module was used to calculate Equation (9).Parallel calculation of the matrix can greatly speed up the interpolation process.In this step, the coefficient matrix [K] needs to be inversed, and the weight vector [λ] is calculated by matrix-vector multiplication.Much progress has been made on the parallelization of matrix operations.Here, an efficient solution that decomposes the coefficient matrix into a lower triangular matrix and an upper triangular matrix by LU decomposition has been chosen before solving a linear system.Then, the system of linear equation decomposition is solved directly by forward and backward substitution without directly calculating the inverse matrix.

OpenCL-Based Implementation
Open Computing Language (OpenCL) is a framework for writing programs on heterogeneous platforms [29].OpenCL assumes that computing systems consist of multiple computing devices, which can be central processing units (CPUs) or hardware accelerators, such as graphics processing units (GPUs) connected to a host processor.OpenCL provides a standard interface for parallel computing based on task and data parallelism, as shown in Figure 3:  The OpenCL execution model consists of the following two components: the host program and the kernel.The host program is executed on the host to define the device context, queue the kernel, and execute kernel instances using the command queue.The kernel is an executable code based on the OpenCL standard that runs on the computing device, similar to data or taskparallel functions.In the spatiotemporal kriging interpolation, the calculation of each vector in the STVMFM and the solution to each unknown point in the OSTKIM are independent.So, parallel kernel functions are written, the calculation is distributed to computing devices, and the computing speed is improved.
As the hardware continues to evolve, CPUs and GPUs vary widely in architecture due to their different positioning.Taking the Intel Core i7-7700 and NVIDIA GTX1060 6 G as examples, the former has a base frequency speed of 3.60 GHz and four physical cores with 8 threads per core, and the latter has a base frequency speed of 1.50 GHz with 1280 cores.This shows that the main difference between the two is that the CPU has a faster clock speed and processes a single calculation faster, and the GPU has more cores and is better at handling multiple calculations.Based on these characteristics, this paper implemented the spatiotemporal kriging interpolation algorithm as shown in Figure 4.This algorithm uses the CPU as a host to achieve an STDOM module, an ASTPSM module as well as to create a context, task queue and data division distributed to the computing device.At the same time, the GPU is used as a computing device to achieve an STVMFM module, an OSTKIM module and to speed up the interpolation process.
When using the OSTKIM, the coefficient matrix [] is the same for the prediction points with the same neighboring points.Memory resources are invaluable in GPU operations.Therefore, for the unknown points that have the same neighboring points, the coefficient matrix [] is stored once, and the vectors [] of the unknown points with the same coefficient matrix are used to construct the matrix for matrix multiplication.This step reduces memory redundancy and avoids high data communication costs, as shown in Figure 5.The OpenCL execution model consists of the following two components: the host program and the kernel.The host program is executed on the host to define the device context, queue the kernel, and execute kernel instances using the command queue.The kernel is an executable code based on the OpenCL standard that runs on the computing device, similar to data or task-parallel functions.In the spatiotemporal kriging interpolation, the calculation of each vector in the STVMFM and the solution to each unknown point in the OSTKIM are independent.So, parallel kernel functions are written, the calculation is distributed to computing devices, and the computing speed is improved.
As the hardware continues to evolve, CPUs and GPUs vary widely in architecture due to their different positioning.Taking the Intel Core i7-7700 and NVIDIA GTX1060 6 G as examples, the former has a base frequency speed of 3.60 GHz and four physical cores with 8 threads per core, and the latter has a base frequency speed of 1.50 GHz with 1280 cores.This shows that the main difference between the two is that the CPU has a faster clock speed and processes a single calculation faster, and the GPU has more cores and is better at handling multiple calculations.Based on these characteristics, this paper implemented the spatiotemporal kriging interpolation algorithm as shown in Figure 4.This algorithm uses the CPU as a host to achieve an STDOM module, an ASTPSM module as well as to create a context, task queue and data division distributed to the computing device.At the same time, the GPU is used as a computing device to achieve an STVMFM module, an OSTKIM module and to speed up the interpolation process.
When using the OSTKIM, the coefficient matrix [K] is the same for the prediction points with the same neighboring points.Memory resources are invaluable in GPU operations.Therefore, for the unknown points that have the same neighboring points, the coefficient matrix [K] is stored once, and the vectors [M] of the unknown points with the same coefficient matrix are used to construct the matrix for matrix multiplication.This step reduces memory redundancy and avoids high data communication costs, as shown in Figure 5.

Computational Environment
To facilitate further research, considering the mainstream configuration of computers on the market, we chose the computing platform with the following parameters to carry out the experiment, as shown in Table 2.The experimental data were the monthly average temperatures of 57 sites from 2007 to 2016 in Beijing, Tianjin, Hebei and its surrounding areas (Figure 6), collected by the meteorological data center of the China Meteorological Administration (http://data.cma.cn/data/cdcdetail/dataCode/B.0021.0002.html).

Computational Environment
To facilitate further research, considering the mainstream configuration of computers on the market, we chose the computing platform with the following parameters to carry out the experiment, as shown in Table 2.The experimental data were the monthly average temperatures of 57 sites from 2007 to 2016 in Beijing, Tianjin, Hebei and its surrounding areas (Figure 6

Experimental Methodology Design
In the above hardware environment, we compared the efficiencies of a CPU sequential program and a GPU-accelerated parallel algorithm.To this end, the following three experiments were designed.
1. Experiment 1: four observation point dataset sizes of 2000, 3000, 4000 and 5000 points were used to interpolate a regular grid (1000 × 1000) at a specified time, and the influence of the number of sampling points was tested to determine the algorithm's efficiency.

Experimental Methodology Design
In the above hardware environment, we compared the efficiencies of a CPU sequential program and a GPU-accelerated parallel algorithm.To this end, the following three experiments were designed.

1.
Experiment 1: four observation point dataset sizes of 2000, 3000, 4000 and 5000 points were used to interpolate a regular grid (1000 × 1000) at a specified time, and the influence of the number of sampling points was tested to determine the algorithm's efficiency.

2.
Experiment 2: the dataset of observations with 5000 points was used to determine the effect of the number of interpolation points on the efficiency of the algorithm for the interpolation of five types of prediction point datasets including space images of one moment with resolutions of 500 × 500, 800 × 800, 1000 × 1000 and 2000 × 2000.

3.
Experiment 3: the results of the space-time interpolation were verified by cross-validation.The detailed cross-validation process was as follows: • Step 1: Delete the first observation value, Z 1 , from the data set.

•
Step 2: Use other observations, the variogram model and the kriging method to predict the point value, Z * 1 .

•
Step Cross-validation is a gradual process, and the best parameters and variation models can be obtained through many comparisons and practices.
The speed increase was calculated according to Equation ( 14), and the interpolation performances were compared.S p is the speed up, T 1 represents the run time serial executed by CPU, and T 2 is the execution time of the parallel accelerated model in the parallel system.The performance of the proposed parallel algorithm was evaluated with S p :

Interpolation Acceleration Effect and Analysis
The results of experiment 1 are shown in Table 3.As the number of observation points increased, the spatiotemporal kriging parallel algorithm did not significantly increase the time consumption, and the acceleration ratio was approximately 3:1.The parallel interpolation method was able to effectively adapt to changes in the number of sampling points and had good stability.With CPU serial operation, because the number of interpolation points increased, the computation time of the 2000 × 2000 prediction point set at a single time was 11 times that of the 500 × 500 prediction point set.When the parallel framework was used, the time consumption for calculating the five kinds of prediction point datasets significantly decreased.When the prediction point dataset was 2000 × 2000, the speed up ratio was up to 3.23.Compared with an existing solution in R gstat [2], in the same experimental environment, when using the dataset of 3600 observation points to interpolate a 1000 × 1000 grid at one moment, it required 89.6 s (Table 4).Obviously, the parallel framework provided faster computing speeds.
Upon further analysis, the interpolation process included the serial part and the parallel part as well as the communication costs.With the increase in the number of parallel parts, the acceleration effect gradually increased.
Based on the two experiments, the parallel algorithm and small mesh interpolation parallel model were verified through data decoupling, increasing the concurrent number and accelerating operation.
This configuration was able to cope with the changes in sampling points and interpolation grids, provide fast operation and solve practical application problems.Experiment 3 shows the cross-validation results of the spatiotemporal kriging interpolation.Only some of the cross-validation results are displayed because of the large number of time and space points.A comparison of the spatiotemporal kriging interpolation model and the pure spatial kriging interpolation model is shown in Figure 7.The actual observed values, compared to the root mean square error (RMSE) of the spatiotemporal kriging interpolation, was 5.38-less than the kriging interpolation space (6.17).As a result, the spatiotemporal kriging interpolation results were closer to the actual measured values.The use of the proposed parallel spatiotemporal kriging acceleration method for interpolation of the Beijing Tianjin Hebei region in September 2009 was examined.
Based on the two experiments, the parallel algorithm and small mesh interpolation parallel model were verified through data decoupling, increasing the concurrent number and accelerating operation.This configuration was able to cope with the changes in sampling points and interpolation grids, provide fast operation and solve practical application problems.

Conclusions
In the study of large-scale spatiotemporal phenomena, there are often missing data.The application of spatiotemporal kriging interpolation to spatiotemporal data interpolation can complement the missing historical data.It is beneficial to obtain the evolution process of temporal and spatial phenomena, which can improve the cognitive level of scholars and establish a spatiotemporal database.
In this paper, to address the shortage of space-time interpolation computing resource consumption and slow speed, a parallel accelerated algorithm based on OpenCL for spatiotemporal kriging interpolation is proposed, which solves the problem.The time dimension of the data is added to rationally and realistically improve the original spatial interpolation.

Conclusions
In the study of large-scale spatiotemporal phenomena, there are often missing data.The application of spatiotemporal kriging interpolation to spatiotemporal data interpolation can complement the missing historical data.It is beneficial to obtain the evolution process of temporal and spatial phenomena, which can improve the cognitive level of scholars and establish a spatiotemporal database.
In this paper, to address the shortage of space-time interpolation computing resource consumption and slow speed, a parallel accelerated algorithm based on OpenCL for spatiotemporal kriging interpolation is proposed, which solves the problem.The time dimension of the data is added to rationally and realistically improve the original spatial interpolation.
The experiments described showed that the method proposed in this paper can greatly improve the interpolation speed, and it can achieve fast computation of spatiotemporal kriging interpolation on a standard PC when compared with some of the existing solutions.The proposed method has great improvement in usability, experience and universality.

Figure 2 .
Figure 2. Construction of the vector matrix of different spaces and dates.

Figure 2 .
Figure 2. Construction of the vector matrix of different spaces and dates.

Figure 4 .
Figure 4. Overall framework for parallel ordinary spatiotemporal kriging implementation.

Figure 4 .
Figure 4. Overall framework for parallel ordinary spatiotemporal kriging implementation.

Figure 4 .
Figure 4. Overall framework for parallel ordinary spatiotemporal kriging implementation.

Figure 6 .
Figure 6.57 Meteorological stations in Beijing, Tianjin, Hebei and surrounding areas.

Figure 6 .
Figure 6.57 Meteorological stations in Beijing, Tianjin, Hebei and surrounding areas.

Experiment 3
shows the cross-validation results of the spatiotemporal kriging interpolation.Only some of the cross-validation results are displayed because of the large number of time and space points.A comparison of the spatiotemporal kriging interpolation model and the pure spatial kriging interpolation model is shown in Figure7.The actual observed values, compared to the root mean square error (RMSE) of the spatiotemporal kriging interpolation, was 5.38-less than the kriging interpolation space (6.17).As a result, the spatiotemporal kriging interpolation results were closer to the actual measured values.The use of the proposed parallel spatiotemporal kriging acceleration method for interpolation of the Beijing Tianjin Hebei region in September 2009 was examined.

Figure 7 .
Figure 7.Comparison of the cross-validation results of the two interpolation models.

Figure 7 .
Figure 7.Comparison of the cross-validation results of the two interpolation models.

Table 1 .
Definition of the spatiotemporal index.

Table 2 .
Detailed configurations of the experimental platforms.

Table 2 .
Detailed configurations of the experimental platforms.
3: Put Z 1 into the dataset and repeat Step 1 to Step 4 until the other predicted values, Through comparison of the original data, Z 1 , Z 2 , . . ., Z n , and predicted data, Z * 1 , Z * 2 , . . ., Z * n , the error is calculated, and the effect of interpolation can be judged.

Table 3 .
Comparison of acceleration of different data sizes, the serial spatiotemporal kriging (SSTK) and the parallel spatiotemporal kriging (PSTK).

Table 4 .
Acceleration at different resolutions.

Table 4 .
Acceleration at different resolutions.