馬駿 孔帥可 周兵 張桐

摘 ?要： 薄板樣條函數(shù)是空間插值中的一種重要方法。對于巨幅影像數(shù)據(jù)使用薄板樣條函數(shù)進(jìn)行空間插值時(shí)，可能會(huì)出現(xiàn)運(yùn)行時(shí)間太長，以及計(jì)算機(jī)內(nèi)存空間不足或程序運(yùn)行無響應(yīng)的問題。針對這些問題，根據(jù)薄板樣條函數(shù)光滑、連續(xù)的特點(diǎn)，基于GDAL開源函數(shù)庫，提出對巨幅影像數(shù)據(jù)的分塊讀取，在塊內(nèi)利用并行技術(shù)求解線性方程組，確定薄板樣條函數(shù)，最后進(jìn)行空間插值的方法。結(jié)果表明，該方法可以有效的解決這些問題。

關(guān)鍵詞： 薄板樣條函數(shù); 空間插值; GDAL; 分塊; 并行

中圖分類號：TP399 ? ? ? ? ?文獻(xiàn)標(biāo)志碼：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 引言

薄板樣條函數(shù)是一種常用插值函數(shù)，是自然樣條函數(shù)在多維空間的推廣，它可以用來表示多維空間曲面。薄板樣條函數(shù)在各個(gè)學(xué)科均有廣泛的應(yīng)用[1]。空間插值是根據(jù)已知點(diǎn)推求一定區(qū)域內(nèi)任意點(diǎn)的方法[2-3]，是將點(diǎn)數(shù)據(jù)轉(zhuǎn)換成面數(shù)據(jù)的一種方法[4]，其任務(wù)是基于已知點(diǎn)來為新的點(diǎn)計(jì)算出最可能的值[5]。常用的空間差值算法有反距離加權(quán)插值法、樣條函數(shù)以及克里金插值[6]。薄板樣條函數(shù)將插值問題模擬為一個(gè)薄金屬板在點(diǎn)約束下的彎曲變形，用簡練的代數(shù)式表示變形的能量，基于點(diǎn)的非線性變換方法，用于離散點(diǎn)數(shù)據(jù)插值得到曲面的一種工具[7]。薄板樣條函數(shù)廣泛應(yīng)用于數(shù)字高程模型的建立、等值線的自動(dòng)制圖以及遙感衛(wèi)星影像的幾何校正。薄板樣條函數(shù)與常用的空間插值方法相比，能夠很好的反映地表高程異常變化的特性，并具備樣條函數(shù)的光滑、連續(xù)、彈性好的特點(diǎn)[8]。在地面高程沒有突變的地區(qū)，用薄板樣條函數(shù)可以很好地描述其形狀，當(dāng)?shù)匦伪容^復(fù)雜時(shí)，可將測區(qū)適當(dāng)分塊[9]?；跀?shù)據(jù)分塊技術(shù)，研究薄板樣條函數(shù)的高效率實(shí)現(xiàn)，具有非常重要的理論與現(xiàn)實(shí)意義。

1 薄板樣條函數(shù)原理

用于空間插值的數(shù)據(jù)通常是復(fù)雜空間變化有限的采樣點(diǎn)測量數(shù)據(jù)，這些已知的測量數(shù)據(jù)稱為“硬數(shù)據(jù)”。在采樣點(diǎn)數(shù)據(jù)比較少的情況下，可以根據(jù)已知的導(dǎo)致某種空間變化的自然過程或現(xiàn)象的信息機(jī)理，輔助進(jìn)行空間插值，這種已知的信息機(jī)理，稱為“軟信息”。但通常情況下，由于不清楚這種自然過程機(jī)理，往往不得不對該問題的屬性在空間的變化作一些假設(shè)，例如假設(shè)采樣點(diǎn)之間的數(shù)據(jù)變化是平滑的，并假設(shè)服從某種分布概率和統(tǒng)計(jì)穩(wěn)定性關(guān)系。采樣點(diǎn)的空間位置對空間插值的結(jié)果影響很大，理想的情況是在研究區(qū)內(nèi)均勻分布。

空間插值常用于將離散點(diǎn)的測量數(shù)據(jù)轉(zhuǎn)換為連續(xù)的數(shù)據(jù)曲面。簡單來說，空間插值是根據(jù)已知點(diǎn)推求一定區(qū)域內(nèi)任意點(diǎn)的方法。實(shí)現(xiàn)插值的核心就是利用薄板樣條函數(shù)通過已知樣本點(diǎn)坐標(biāo)來推求出一定區(qū)域內(nèi)任意點(diǎn)的坐標(biāo)。對于給定的樣本點(diǎn)（xi，yi），薄板樣條函數(shù)可定義為：

⑴

包涵的條件：

⑵

其中：

⑶

f（xi，yi）=zi ? ⑷

根據(jù)薄板樣條函數(shù)，在對數(shù)據(jù)插值處理前，首先選取N（N>=7）個(gè)非零樣本點(diǎn)，每個(gè)樣本點(diǎn)可以用它的坐標(biāo)（x，y）以及該坐標(biāo)的圖像灰度值z來表示成（xi，yi，zi），然后，將得到的N個(gè)樣本點(diǎn)值帶入到薄板樣條函數(shù)f（x，y）中得到N個(gè)對應(yīng)方程。每個(gè)方程中含有N+3個(gè)未知量：b1，b2，b3…bn，a0，a1，a2。根據(jù)公式⑵得到三個(gè)等式，這樣便構(gòu)造出一個(gè)N+3元方程組。對得到的N+3元方程組求解，得到N+3個(gè)解，對應(yīng)于薄板樣條函數(shù)f（x，y）的系數(shù)：b1，b2，b3…bn，a0，a1，a2。這樣便根據(jù)所選取已知樣本點(diǎn)的數(shù)據(jù)確定了薄板樣條函數(shù)f（x，y）。對于需要插值的點(diǎn)（xi，yi）便能根據(jù)公式⑷求得插值后該點(diǎn)的值zi。

2 基于薄板樣條插值算法的分塊并行插值處理

2.1 分析

根據(jù)薄板樣條函數(shù)的原理，要完成空間插值，可以分為兩個(gè)步驟，第一步是確定薄板樣條函數(shù)的系數(shù)，即求解線性方程組;第二步是根據(jù)第一步確定的系數(shù)，計(jì)算未知點(diǎn)的值。薄板樣條函數(shù)被廣泛用于空間插值，然而，由于計(jì)算效率問題，它對于大型數(shù)據(jù)集的應(yīng)用卻是很有限。薄板樣條函數(shù)分析計(jì)算的時(shí)間復(fù)雜度為O（n3），其中n是數(shù)據(jù)點(diǎn)的個(gè)數(shù)，因此對于超過2000個(gè)數(shù)據(jù)點(diǎn)集的常規(guī)計(jì)算變得不可行[10]。對于巨幅影像數(shù)據(jù)，由于計(jì)算機(jī)的性能限制，在程序運(yùn)行時(shí)可能會(huì)出現(xiàn)運(yùn)行時(shí)間太長和內(nèi)存不足的問題。以高分1號衛(wèi)星WFV傳感器拍攝的16米影像來說，單幅影像的尺寸大小為12000*13400，占用磁盤空間為2G左右，若以薄板樣條函數(shù)對其進(jìn)行幾何校正，需要將整幅影像讀入到計(jì)算機(jī)內(nèi)存中，并且構(gòu)造的系數(shù)矩陣會(huì)相當(dāng)?shù)凝嫶?，對于個(gè)人計(jì)算機(jī)而言，嚴(yán)重影響系統(tǒng)的運(yùn)行，甚至?xí)霈F(xiàn)應(yīng)用程序無響應(yīng)的情況。針對這一問題，考慮到薄板樣條函數(shù)光滑、連續(xù)、彈性好的特點(diǎn)，提出將影像分塊讀取，然后對于每個(gè)塊內(nèi)的數(shù)據(jù)再采用并行運(yùn)算，確定塊內(nèi)數(shù)據(jù)的函數(shù)的系數(shù)值，然后對塊內(nèi)的無效值進(jìn)行插值運(yùn)算。這樣數(shù)據(jù)分塊可以有效的解決內(nèi)存問題及方程組系數(shù)矩陣龐大的問題，同時(shí)塊內(nèi)數(shù)據(jù)并行運(yùn)算也很大程度的提高了程序的運(yùn)行效率。分塊讀取影像文件的流程如圖1所示。

2.2 實(shí)現(xiàn)

2.2.1 數(shù)據(jù)分塊讀取

假設(shè)影像的長度為Length，寬度為Width，考慮到分塊結(jié)果要進(jìn)行矩陣運(yùn)算，所以分塊的長度和寬度要保持相等，假設(shè)分塊長度和寬度都為sideLength，則整幅影像可被分為n塊。其中n=LNum*WNum，LNum=Length/sideLength+1，WNum=Width/sideLenght+1。第i個(gè)分塊的左上角的位置對應(yīng)于原始影像左上角的偏移量為：水平方向xoff=（i-1）%LNum*sideLength，豎直方向yoff=（i-1）/LNum*sideLength。對于影像邊緣位置，若分塊已超過影像的大小，則對分塊的左上角進(jìn)行修正：水平方向上，若xoff+sideLength>Length，則xoff=Length-sideLength，豎直方向上，若yoff+sideLength>Width，則yoff=Width-sideLength。在進(jìn)行分塊數(shù)據(jù)讀取時(shí)，調(diào)用開源的GDAL函數(shù)庫，具體函數(shù)原型如下：

圖1 ?分塊讀取影像文件流程圖

2.2.2 塊內(nèi)數(shù)據(jù)并行運(yùn)算

要完成空間數(shù)據(jù)插值，最主要的便是確定薄板樣條函數(shù)的系數(shù)，即求解線性方程組。因此，塊內(nèi)數(shù)據(jù)并行運(yùn)算，最終歸結(jié)為并行求解線性方程組。對于并行求解線性方程組[11]，我們可以找到大量的求解方法，這里我們采用已經(jīng)成熟的LU并行分解算法及三角方程組的并行解法，考慮到處理機(jī)的負(fù)載均衡，兩種解法都采取卷簾方式[12]存放數(shù)據(jù)，n*n的矩陣A在處理機(jī)上的存放方式如表1所示，其中m=n/t。

表1 ?矩陣A在處理機(jī)上的存放方式

有大量的試驗(yàn)證明，這兩種解法在求解線性方程組方面都取得了很好的效果。對于這兩種解法的具體實(shí)現(xiàn)，本文不再贅述。

2.3 測試結(jié)果

對北京某地區(qū)6001*6001的DEM影像數(shù)據(jù)進(jìn)行空間插值時(shí)，表2和表3分別為分塊大小為10*10和分塊大小為100*100的相同區(qū)域數(shù)據(jù)（試驗(yàn)條件：windows7 x64系統(tǒng)，4G內(nèi)存）。

表2 ?分塊大小為10*10插值后DEM值

表3 ?分塊大小為100*100插值后DEM值

通過對比兩表數(shù)據(jù)，誤差在允許范圍內(nèi)，分塊大小為10*10與分塊大小為100*100的插值結(jié)果幾乎一樣。

對于100*100的矩陣A運(yùn)用串行LU分解法求解方程組時(shí)，運(yùn)行時(shí)間達(dá)數(shù)小時(shí)之久，采用LU并行分解算法及三角方程組的并行解法，運(yùn)行效率得到大大提高，程序運(yùn)行時(shí)間大大縮小。

經(jīng)過多次測試分析，得出分塊大小與算法的精度如圖2所示。

圖2 ?分塊大小與算法精度關(guān)系

當(dāng)每一塊比較小時(shí)，塊內(nèi)的有效值就相對比較少，參與構(gòu)造薄板樣條函數(shù)的系數(shù)矩陣的點(diǎn)就會(huì)比較少，因此計(jì)算的結(jié)果誤差較大，而且每一塊比較小時(shí)，無法充分利用LU并行分解算法及三角方程組的并行解法的性能，分塊數(shù)n也會(huì)比較大，所以在運(yùn)行時(shí)間上也會(huì)比較長。當(dāng)每一塊比較大時(shí)，塊內(nèi)的有效值也相對比較多，參與構(gòu)造薄板樣條函數(shù)的系數(shù)矩陣的點(diǎn)就會(huì)比較多，因此計(jì)算的結(jié)果精度就比較高，但是如果每一塊特別大，那么塊內(nèi)構(gòu)造的矩陣就會(huì)比較大，求解線性方程組的時(shí)間就會(huì)相對比較長。由于計(jì)算機(jī)性能及影像數(shù)據(jù)的尺寸大小不同，要達(dá)到最佳性能，分塊大小也不相同，建議分塊大小在1000*1000至1500*1500范圍內(nèi)。

3 結(jié)束語

本文根據(jù)薄板樣條函數(shù)光滑、連續(xù)、彈性好等特點(diǎn)，針對大影像數(shù)據(jù)，提出分塊讀取，對分塊內(nèi)數(shù)據(jù)創(chuàng)建系數(shù)矩陣，運(yùn)用較為成熟的LU并行分解算法及三角方程組的并行解法確定薄板樣條函數(shù)系數(shù)，完成每個(gè)數(shù)據(jù)塊內(nèi)無效值的計(jì)算，既解決了對巨幅影像數(shù)據(jù)空間插值時(shí)的計(jì)算機(jī)內(nèi)存空間不足的問題，同時(shí)也提高了程序運(yùn)行的效率。

參考文獻(xiàn)：

[1] 孫海燕，黃勝.薄板樣條函數(shù)逐次增加節(jié)點(diǎn)的算法[J].測繪工程，

2006.15（3）：12-14

[2] 龔健雅，杜道生.當(dāng)代地理信息技術(shù)[M].科學(xué)出版社，2004.

[3] 黎夏，劉凱.GIS與空間分析——原理與方法[M].科學(xué)出版社，2006.

[4] Kang-tsung Chang著，陳健飛等譯.地理信息系統(tǒng)導(dǎo)論（第3版）[M].

北京清華大學(xué)出版社，2009.

[5] Tor Bernhardsen著，王滸，李浩川譯.地理信息系統(tǒng)導(dǎo)論[M].機(jī)械工

業(yè)出版社，2006.

[6] 黃舸.多種空間插值法在連續(xù)分布環(huán)境要素分析中的應(yīng)用精度比較[J].

環(huán)?？萍?，2012.3：35-37

[7] 杜國明，賈良文.薄板樣條函數(shù)在空間數(shù)據(jù)插值中的應(yīng)用[J].計(jì)算機(jī)

工程與應(yīng)用，2009.45（36）：238-240

[8] 岳云娟，張慶敏，白林.薄板樣條函數(shù)在城市三維地質(zhì)建模中的應(yīng)用[J].

四川理工學(xué)院學(xué)報(bào)（自然科學(xué)版），2012.25（2）：28-30

[9] 孫海燕，丁咚.薄板樣條函數(shù)及復(fù)雜曲面的數(shù)學(xué)表示[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] 遲學(xué)斌.在具有局部內(nèi)存與共享主存的并行機(jī)上并行求解線性方

程組[J].計(jì)算數(shù)學(xué)，1995.17（2）：210-217

[12] 張健飛，姜弘道.對稱正定矩陣的并行LDLT分解算法實(shí)現(xiàn)[J].計(jì)算

機(jī)工程與設(shè)計(jì)，2003.24（10）：75-77