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）是指利用一组源面单元（source areal units）上的已知属性值，预测另一组目标面单元（target areal units）上未知属性值的过程。地理统计面插值（Geostatistical Areal Interpolation）采用空间预测算法，即克里金（Kriging）算法的变体，该类算法在插值过程中显式纳入了源与目标单元之间的空间自相关与尺度差异。当使用所有可用的源观测数据进行插值时（即采用全局搜索邻域时），地理统计面插值的计算量极大。即便借助Kyriakidis等人（2005）与Liu等人（2006）基于快速傅里叶变换（Fast Fourier Transform, FFT）提出的谱算法所带来的性能大幅提升，此类插值仍需占用海量内存与计算资源：其中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的并行地理统计面插值算法，以解决此类问题的计算瓶颈。该算法包含三个并行流程：（1）通过FFT计算源-源与源-目标协方差矩阵；（2）对源-源协方差矩阵进行QR分解；（3）通过克里金算法计算源-目标权重，并随后计算目标支撑域的预测属性值。基于真实世界数据集的实验（即从美国东部时区及美国大陆地区的县域人口密度，预测流域的人口密度）结果表明，该并行算法将计算时间大幅缩减至实际空间分析应用可接受的范围，并实现了较高的加速比与计算效率。实验同时证明，随着处理器数量与问题规模的提升，该算法均能保持良好的可扩展性。