基于薄板样条插值算法的巨幅影像分块并行处理

摘 要： 薄板样条函数是空间插值中的一种重要方法。对于巨幅影像数据使用薄板样条函数进行空间插值时，可能会出现运行时间太长，以及计算机内存空间不足或程序运行无响应的问题。针对这些问题，根据薄板样条函数光滑、连续的特点，基于GDAL开源函数库，提出对巨幅影像数据的分块读取，在块内利用并行技术求解线性方程组，确定薄板样条函数，最后进行空间插值的方法。结果表明，该方法可以有效的解决这些问题。

关键词： 薄板样条函数; 空间插值; GDAL; 分块; 并行

中图分类号：TP399 文献标志码：A 文章编号：1006-8228（2015）07-04-03

Huge image block parallel processing based on thin plate spline interpolation algorithm

Ma Jun1，2， Kong Shuaike1， Zhou Bing1，2， Zhang Tong1

（1. Department of Computer Science and Technology， Henan University， Kaifeng， Henan 475004， China;

2. Data and Knowledge Engineering Research Institute， Henan University）

Abstract： Thin plate spline is an important algorithm of spatial interpolation. For huge image data when using thin plate spline interpolation， the problem of running too long and insufficient computer memory space or program runs no response may occur.. To solve the problem， according to the smooth， continuous feature of thin plate spline， and based on GDAL open source library， proposes a method， in which the huge image data is divided into blocks， parallel technology is used to solve the linear equations in the block， and then interpolates with thin plate spline. The results indicate that it is an effective method to solve the problem.

Key words： thin plate spline; spatial interpolation; GDAL; parallel processing

0 引言

薄板样条函数是一种常用插值函数，是自然样条函数在多维空间的推广，它可以用来表示多维空间曲面。薄板样条函数在各个学科均有广泛的应用[1]。空间插值是根据已知点推求一定区域内任意点的方法[2-3]，是将点数据转换成面数据的一种方法[4]，其任务是基于已知点来为新的点计算出最可能的值[5]。常用的空间差值算法有反距离加权插值法、样条函数以及克里金插值[6]。薄板样条函数将插值问题模拟为一个薄金属板在点约束下的弯曲变形，用简练的代数式表示变形的能量，基于点的非线性变换方法，用于离散点数据插值得到曲面的一种工具[7]。薄板样条函数广泛应用于数字高程模型的建立、等值线的自动制图以及遥感卫星影像的几何校正。薄板样条函数与常用的空间插值方法相比，能够很好的反映地表高程异常变化的特性，并具备样条函数的光滑、连续、弹性好的特点[8]。在地面高程没有突变的地区，用薄板样条函数可以很好地描述其形状，当地形比较复杂时，可将测区适当分块[9]。基于数据分块技术，研究薄板样条函数的高效率实现，具有非常重要的理论与现实意义。

1 薄板样条函数原理

用于空间插值的数据通常是复杂空间变化有限的采样点测量数据，这些已知的测量数据称为“硬数据”。在采样点数据比较少的情况下，可以根据已知的导致某种空间变化的自然过程或现象的信息机理，辅助进行空间插值，这种已知的信息机理，称为“软信息”。但通常情况下，由于不清楚这种自然过程机理，往往不得不对该问题的属性在空间的变化作一些假设，例如假设采样点之间的数据变化是平滑的，并假设服从某种分布概率和统计稳定性关系。采样点的空间位置对空间插值的结果影响很大，理想的情况是在研究区内均匀分布。

空间插值常用于将离散点的测量数据转换为连续的数据曲面。简单来说，空间插值是根据已知点推求一定区域内任意点的方法。实现插值的核心就是利用薄板样条函数通过已知样本点坐标来推求出一定区域内任意点的坐标。对于给定的样本点（xi，yi），薄板样条函数可定义为：

⑴

包涵的条件：

⑵

其中：

⑶

f（xi，yi）=zi ⑷

根据薄板样条函数，在对数据插值处理前，首先选取N（N>=7）个非零样本点，每个样本点可以用它的坐标（x，y）以及该坐标的图像灰度值z来表示成（xi，yi，zi），然后，将得到的N个样本点值带入到薄板样条函数f（x，y）中得到N个对应方程。每个方程中含有N+3个未知量：b1，b2，b3…bn，a0，a1，a2。根据公式⑵得到三个等式，这样便构造出一个N+3元方程组。对得到的N+3元方程组求解，得到N+3个解，对应于薄板样条函数f（x，y）的系数：b1，b2，b3…bn，a0，a1，a2。这样便根据所选取已知样本点的数据确定了薄板样条函数f（x，y）。对于需要插值的点（xi，yi）便能根据公式⑷求得插值后该点的值zi。

2 基于薄板样条插值算法的分块并行插值处理

2.1 分析

根据薄板样条函数的原理，要完成空间插值，可以分为两个步骤，第一步是确定薄板样条函数的系数，即求解线性方程组;第二步是根据第一步确定的系数，计算未知点的值。薄板样条函数被广泛用于空间插值，然而，由于计算效率问题，它对于大型数据集的应用却是很有限。薄板样条函数分析计算的时间复杂度为O（n3），其中n是数据点的个数，因此对于超过2000个数据点集的常规计算变得不可行[10]。对于巨幅影像数据，由于计算机的性能限制，在程序运行时可能会出现运行时间太长和内存不足的问题。以高分1号卫星WFV传感器拍摄的16米影像来说，单幅影像的尺寸大小为12000*13400，占用磁盘空间为2G左右，若以薄板样条函数对其进行几何校正，需要将整幅影像读入到计算机内存中，并且构造的系数矩阵会相当的庞大，对于个人计算机而言，严重影响系统的运行，甚至会出现应用程序无响应的情况。针对这一问题，考虑到薄板样条函数光滑、连续、弹性好的特点，提出将影像分块读取，然后对于每个块内的数据再采用并行运算，确定块内数据的函数的系数值，然后对块内的无效值进行插值运算。这样数据分块可以有效的解决内存问题及方程组系数矩阵庞大的问题，同时块内数据并行运算也很大程度的提高了程序的运行效率。分块读取影像文件的流程如图1所示。

2.2 实现

2.2.1 数据分块读取

假设影像的长度为Length，宽度为Width，考虑到分块结果要进行矩阵运算，所以分块的长度和宽度要保持相等，假设分块长度和宽度都为sideLength，则整幅影像可被分为n块。其中n=LNum*WNum，LNum=Length/sideLength+1，WNum=Width/sideLenght+1。第i个分块的左上角的位置对应于原始影像左上角的偏移量为：水平方向xoff=（i-1）%LNum*sideLength，竖直方向yoff=（i-1）/LNum*sideLength。对于影像边缘位置，若分块已超过影像的大小，则对分块的左上角进行修正：水平方向上，若xoff+sideLength>Length，则xoff=Length-sideLength，竖直方向上，若yoff+sideLength>Width，则yoff=Width-sideLength。在进行分块数据读取时，调用开源的GDAL函数库，具体函数原型如下：

图1 分块读取影像文件流程图

2.2.2 块内数据并行运算

要完成空间数据插值，最主要的便是确定薄板样条函数的系数，即求解线性方程组。因此，块内数据并行运算，最终归结为并行求解线性方程组。对于并行求解线性方程组[11]，我们可以找到大量的求解方法，这里我们采用已经成熟的LU并行分解算法及三角方程组的并行解法，考虑到处理机的负载均衡，两种解法都采取卷帘方式[12]存放数据，n*n的矩阵A在处理机上的存放方式如表1所示，其中m=n/t。

表1 矩阵A在处理机上的存放方式

有大量的试验证明，这两种解法在求解线性方程组方面都取得了很好的效果。对于这两种解法的具体实现，本文不再赘述。

2.3 测试结果

对北京某地区6001*6001的DEM影像数据进行空间插值时，表2和表3分别为分块大小为10*10和分块大小为100*100的相同区域数据（试验条件：windows7 x64系统，4G内存）。

表2 分块大小为10*10插值后DEM值

表3 分块大小为100*100插值后DEM值

通过对比两表数据，误差在允许范围内，分块大小为10*10与分块大小为100*100的插值结果几乎一样。

对于100*100的矩阵A运用串行LU分解法求解方程组时，运行时间达数小时之久，采用LU并行分解算法及三角方程组的并行解法，运行效率得到大大提高，程序运行时间大大缩小。

经过多次测试分析，得出分块大小与算法的精度如图2所示。

图2 分块大小与算法精度关系

当每一块比较小时，块内的有效值就相对比较少，参与构造薄板样条函数的系数矩阵的点就会比较少，因此计算的结果误差较大，而且每一块比较小时，无法充分利用LU并行分解算法及三角方程组的并行解法的性能，分块数n也会比较大，所以在运行时间上也会比较长。当每一块比较大时，块内的有效值也相对比较多，参与构造薄板样条函数的系数矩阵的点就会比较多，因此计算的结果精度就比较高，但是如果每一块特别大，那么块内构造的矩阵就会比较大，求解线性方程组的时间就会相对比较长。由于计算机性能及影像数据的尺寸大小不同，要达到最佳性能，分块大小也不相同，建议分块大小在1000*1000至1500*1500范围内。

3 结束语

本文根据薄板样条函数光滑、连续、弹性好等特点，针对大影像数据，提出分块读取，对分块内数据创建系数矩阵，运用较为成熟的LU并行分解算法及三角方程组的并行解法确定薄板样条函数系数，完成每个数据块内无效值的计算，既解决了对巨幅影像数据空间插值时的计算机内存空间不足的问题，同时也提高了程序运行的效率。

参考文献：

[1] 孙海燕，黄胜.薄板样条函数逐次增加节点的算法[J].测绘工程，

2006.15（3）：12-14

[2] 龚健雅，杜道生.当代地理信息技术[M].科学出版社，2004.

[3] 黎夏，刘凯.GIS与空间分析――原理与方法[M].科学出版社，2006.

[4] Kang-tsung Chang著，陈健飞等译.地理信息系统导论（第3版）[M].

北京清华大学出版社，2009.

[5] Tor Bernhardsen著，王浒，李浩川译.地理信息系统导论[M].机械工

业出版社，2006.

[6] 黄舸.多种空间插值法在连续分布环境要素分析中的应用精度比较[J].

环保科技，2012.3：35-37

[7] 杜国明，贾良文.薄板样条函数在空间数据插值中的应用[J].计算机

工程与应用，2009.45（36）：238-240

[8] 岳云娟，张庆敏，白林.薄板样条函数在城市三维地质建模中的应用[J].

四川理工学院学报（自然科学版），2012.25（2）：28-30

[9] 孙海燕，丁咚.薄板样条函数及复杂曲面的数学表示[J].测绘工程，

2006.15（2）：7-8

[10] P.A.Hancock，M.F.Hutchinson. Spatial interpolation of large

climate data sets using bivariate thin plate smoothing splines[J]. Environmental Modelling and Software，2005.21（12）：1684-1694

[11] 迟学斌.在具有局部内存与共享主存的并行机上并行求解线性方

程组[J].计算数学，1995.17（2）：210-217

[12] 张健飞，姜弘道.对称正定矩阵的并行LDLT分解算法实现[J].计算

机工程与设计，2003.24（10）：75-77