A parallel computing approach to fast geostatistical areal interpolation

Areal interpolation is the procedure of using known attribute values at a set of (source) areal units to predict unknown attribute values at another set of (target) units. Geostatistical areal interpolation employs spatial prediction algorithms, that is, variants of Kriging, which explicitly incorporate spatial autocorrelation and scale differences between source and target units in the interpolation endeavor. When all the available source measurements are used for interpolation, that is, when a global search neighborhood is adopted, geostatistical areal interpolation is extremely computationally intensive. Interpolation in this case requires huge memory space and massive computing power, even with the dramatic improvement introduced by the spectral algorithms developed by Kyriakidis et al. (2005. Improving spatial data interoperability using geostatistical support-to-support interpolation. In: Proceedings of geoComputation. Ann Arbor, MI: University of Michigan) and Liu et al. (2006. Calculation of average covariance using fast Fourier transform (FFT). Menlo Park, CA: Stanford Center for Reservoir Forecasting, Petroleum Engineering Department, Stanford University) based on the fast Fourier transform (FFT). In this study, a parallel FFT-based geostatistical areal interpolation algorithm was developed to tackle the computational challenge of such problems. The algorithm includes three parallel processes: (1) the computation of source-to-source and source-to-target covariance matrices by means of FFT; (2) the QR factorization of the source-to-source covariance matrix; and (3) the computation of source-to-target weights via Kriging, and the subsequent computation of predicted attribute values for the target supports. Experiments with real-world datasets (i.e., predicting population densities of watersheds from population densities of counties in the Eastern Time Zone and in the continental United States) showed that the parallel algorithm drastically reduced the computing time to a practical length that is feasible for actual spatial analysis applications, and achieved fairly high speed-ups and efficiencies. Experiments also showed the algorithm scaled reasonably well as the number of processors increased and as the problem size increased.面插值（Areal Interpolation）指利用一组（源）面状单元的已知属性值，预测另一组（目标）面状单元的未知属性值的过程。地质统计学面插值采用空间预测算法，即克里金（Kriging）的变体，该算法在插值过程中显式纳入空间自相关与源区、目标单元间的尺度差异。当使用全部可用源数据进行插值（即采用全局搜索邻域）时，地质统计学面插值的计算量极为庞大。即便采用Kyriakidis等人（2005年，《Improving spatial data interoperability using geostatistical support-to-support interpolation》，载于geoComputation会议论文集，密歇根州安娜堡：密歇根大学）与Liu等人（2006年，《Calculation of average covariance using fast Fourier transform (FFT)》，加利福尼亚州门洛帕克：斯坦福大学石油工程系储层预测中心）提出的基于FFT的谱算法及相关改进，此类插值仍需海量内存与超强算力。本研究开发了一种基于并行FFT的地质统计学面插值算法，以解决此类问题的计算难题。该算法包含三个并行流程：（1）通过FFT计算源区-源区及源区-目标区协方差矩阵；（2）对源区-源区协方差矩阵进行QR分解；（3）通过克里金法计算源区-目标区权重，并随后计算目标支撑单元的预测属性值。基于真实世界数据集的试验（即从美国东部时区及美国本土的县人口密度预测流域人口密度）表明，该并行算法将计算时间大幅缩减至实际空间分析应用可接受的合理时长，并实现了较高的加速比与运算效率。试验同时证实，随着处理器数量与问题规模的增长，该算法的扩展性表现良好。