摘 要： 薄板样条函数是空间插值中的一种重要方法。对于巨幅影像数据使用薄板样条函数进行空间插值时，可能会出现运行时间太长，以及计算机内存空间不足或程序运行无响应的问题。针对这些问题，根据薄板样条函数光滑、连续的特点，基于GDAL开源函数库，提出对巨幅影像数据的分块读取，在块内利用并行技术求解线性方程组，确定薄板样条函数，最后进行空间插值的方法。结果表明，该方法可以有效的解决这些问题。
关键词： 薄板样条函数; 空间插值; GDAL; 分块; 并行
Huge image block parallel processing based on thin plate spline interpolation algorithm
（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 spanided 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 薄板样条函数原理
用于空间插值的数据通常是复杂空间变化有限的采样点测量数据，这些已知的测量数据称为“硬数据”。在采样点数据比较少的情况下，可以根据已知的导致某种空间变化的自然过程或现象的信息机理，辅助进行空间插值，这种已知的信息机理，称为“软信息”。但通常情况下，由于不清楚这种自然过程机理，往往不得不对该问题的属性在空间的变化作一些假设，例如假设采样点之间的数据变化是平滑的，并假设服从某种分布概率和统计稳定性关系。采样点的空间位置对空间插值的结果影响很大，理想的情况是在研究区内均匀分布。
空间插值常用于将离散点的测量数据转换为连续的数据曲面。简单来说，空间插值是根据已知点推求一定区域内任意点的方法。实现插值的核心就是利用薄板样条函数通过已知样本点坐标来推求出一定区域内任意点的坐标。对于给定的样本点（xi，yi），薄板样条函数可定义为：
⑴
包涵的条件：
⑵
其中：
⑶
f（xi，yi）=zi ⑷
2.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+sideLengthLength，则xoff=Length-sideLength，竖直方向上，若yoff+sideLengthWidth，则yoff=Width-sideLength。在进行分块数据读取时，调用开源的GDAL函数库，具体函数原型如下：
图1 分块读取影像文件流程图
2.2.2 块内数据并行运算
表1 矩阵A在处理机上的存放方式
有大量的试验证明，这两种解法在求解线性方程组方面都取得了很好的效果。对于这两种解法的具体实现，本文不再赘述。
2.3 测试结果
表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].测绘工程，
[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].
[9] 孙海燕，丁咚.薄板样条函数及复杂曲面的数学表示[J].测绘工程，
[10] P.A.Hancock，M.F.Hutchinson. Spatial interpolation of large
[11] 迟学斌.在具有局部内存与共享主存的并行机上并行求解线性方
[12] 张健飞，姜弘道.对称正定矩阵的并行LDLT分解算法实现[J].计算
机工程与设计，2003.24（10）：75-77
