杜國(guó)明，曾志芳

(1.中山大學(xué)地理科學(xué)與規(guī)劃學(xué)院 廣東省城市化與地理環(huán)境空間模擬重點(diǎn)實(shí)驗(yàn)室，廣東 廣州 510275;2.中山大學(xué)圖書(shū)館，廣東 廣州 510275)

草地植被與降水關(guān)系緊密，降水的減少是造成草地退化的重要因素[1]，對(duì)二者關(guān)系的研究一直是一個(gè)活躍的研究領(lǐng)域[2-3]。尤其我國(guó)是一個(gè)干旱、缺水嚴(yán)重的國(guó)家，人均水資源只有2 200 m3，僅為世界平均水平的四分之一[4]，是全球13個(gè)人均水資源最貧乏的國(guó)家之一，且分布極不均衡：東南多，西北少。研究降水量的空間分布對(duì)草地資源分布、區(qū)域水文和水資源分析以及區(qū)域水資源管理、旱澇災(zāi)害管理、生態(tài)環(huán)境治理都具有重要意義[5]。受觀測(cè)臺(tái)站數(shù)量及分布的限制，空間插值是十分必要的。常用插值方法有普通克里金、樣條法、反距離加權(quán)插值法等[6-7]，但在常用插值方法中通常忽略了高程值的影響，而草地植被、降水量與高程等地形因子有著密切的相關(guān)性，這樣造成較大的插值誤差，尤其對(duì)我國(guó)地形復(fù)雜多樣、山區(qū)面積廣大的區(qū)域，在研究降水量等氣候要素的空間分布時(shí)更應(yīng)該考慮高程等地形因子的影響。為了彌補(bǔ)這一不足，本研究利用三維薄板樣條函數(shù)(Three-Dimensional Thin Plate Splines，簡(jiǎn)稱3D TPS)進(jìn)行中國(guó)降水量空間數(shù)據(jù)的插值研究。

薄板樣條函數(shù)(TPS)[8-10]作為一種空間插值方法，尋找一個(gè)通過(guò)所有控制點(diǎn)的彎曲最小(由一個(gè)能量函數(shù)定義)的光滑曲面，就像一個(gè)薄鐵板，通過(guò)給定的幾個(gè)“樣條”，使鐵板表面盡量光滑。TPS具有光滑、連續(xù)、彈性好的特點(diǎn)，且不要求測(cè)量點(diǎn)呈規(guī)格網(wǎng)排列，從而廣泛應(yīng)用于各種曲面的擬合以及復(fù)雜曲面的數(shù)學(xué)表示[11]。3D TPS除了考慮位置坐標(biāo)(x,y)，還充分考慮高程值z(mì)，z不是協(xié)變量而是同位置坐標(biāo)(x,y)一樣，以自變量形式存在，這樣，3D坐標(biāo)(x,y,z)共同參與插值過(guò)程，這對(duì)于地形復(fù)雜區(qū)域的空間插值研究具有重要意義。

1 三維薄板樣條函數(shù)插值原理

1.1數(shù)據(jù) 數(shù)據(jù)包括中國(guó)國(guó)界圖、中國(guó)數(shù)字高程模型(DEM)圖和中國(guó)氣象站分布圖。其中，中國(guó)國(guó)界圖為矢量圖，可分為面狀圖層和線狀圖層兩種。中國(guó)DEM圖是在美國(guó)地質(zhì)調(diào)查局(USGS)的地球資源觀測(cè)衛(wèi)星數(shù)據(jù)中心生產(chǎn)的全球DEM基礎(chǔ)上，結(jié)合中國(guó)國(guó)界圖截取生成的DEM，它是GRID格式數(shù)據(jù)，本研究采用的分辨率為2.5分(即0.041 67度)，大約相當(dāng)于5 km。中國(guó)氣象站分布圖以矢量點(diǎn)狀圖shape格式存儲(chǔ)，在圖形文件中包含每個(gè)站點(diǎn)的位置(x,y)信息，在屬性數(shù)據(jù)庫(kù)中存放年均降水量等屬性信息，內(nèi)含1951―2007年年均降水量值，共有454個(gè)氣象站。在ArcMap中將中國(guó)DEM圖和中國(guó)氣象站點(diǎn)分布圖疊加，結(jié)果如圖1所示。

圖1 中國(guó)數(shù)字高程及氣象站空間分布圖Fig.1 Digital elevation model and spatial distribution of meteorological stations in China

1.23D TPS原理介紹 假定有n個(gè)數(shù)據(jù)點(diǎn)的觀測(cè)值r(xi,yi,zi)在三維坐標(biāo)點(diǎn)(xi,yi,zi)處，則觀測(cè)值可用如下公式表示[8-13]：

r(xi,yi,zi)=f(xi,yi,zi)+ε(xi,yi,zi),i=1,2,…,n.

(1)

式中，f(xi,yi,zi)是觀測(cè)值r(xi,yi,zi)對(duì)應(yīng)的估計(jì)函數(shù)，ε(xi,yi,zi)是估計(jì)誤差，它是零均值、非相關(guān)隨機(jī)誤差；(xi,yi,zi)是三維歐氏空間的坐標(biāo)。

若要估計(jì)函數(shù)f最大限度地接近觀測(cè)值，可以通過(guò)最小化公式ξi得出，公式ξi如下：

λJm(f)

(2)

式中，平滑參數(shù)λ>0，Jm(f)是三維空間關(guān)于估計(jì)函數(shù)f的平滑項(xiàng)，被定義為：

(3)

f(xi,yi,zi)可通過(guò)以下公式得出：

(4)

式中，m是維度，在3D TPS中的m=4。φj是多項(xiàng)式，其中φ1=1,φ2=x,φ3=Y,φ4=z。參數(shù)aj和bi可由下式求出：

(5)

式中，a=(a1,a2,…,am)，b=(b1,b2,…,bn)T，v=(v1,v2,…,vn)。E是n×n矩陣，Eij=U(dij)。T是m×n矩陣，Tij=φ(tj)。U是徑向基函數(shù)，dj是3D坐標(biāo)點(diǎn)(x,y,z)和(xi,yi,zi)間的歐氏距離，公式表示為：

(6)

dij是(xi,yi,zi)和(xj,yj,zj)間的歐氏距離，可通過(guò)下式得出：

(7)

1.3方法與步驟 利用ArcGIS中的ArcMap完成數(shù)據(jù)預(yù)處理及成圖，在Visual Basic.Net平臺(tái)上利用ArcGIS Engine進(jìn)行二次開(kāi)發(fā)完成數(shù)據(jù)提取、處理和運(yùn)算等。

1.3.1數(shù)據(jù)預(yù)處理 根據(jù)全球DEM，利用ArcGIS，通過(guò)中國(guó)國(guó)界圖掩膜，得到中國(guó)DEM。將中國(guó)氣象站點(diǎn)以shape格式點(diǎn)狀圖層形式存儲(chǔ)，(x,y)坐標(biāo)存放在圖形文件中，年均降水量等屬性存放在屬性數(shù)據(jù)庫(kù)中。

1.3.2提取氣象站點(diǎn)的3D坐標(biāo)(x,y,z)及年均降水量 在Visual Basic.Net平臺(tái)上利用ArcGIS Engine的二次開(kāi)發(fā)技術(shù)，獲取2D坐標(biāo)(x,y)和年均降水量P。高程值z(mì)可通過(guò)提取屬性數(shù)據(jù)庫(kù)的屬性信息得到，也可通過(guò)DEM獲取。兩種方法之間有一定誤差，由于在插值時(shí)，用DEM提取的高程值代替氣象站點(diǎn)的實(shí)測(cè)海拔可以提高插值精度[14]，所以本研究采用DEM獲取高程值的方法。

1.3.3格網(wǎng)化數(shù)據(jù) 先均勻生成等間距的格網(wǎng)，再對(duì)每個(gè)格網(wǎng)點(diǎn)進(jìn)行空間插值。另外，還可采用中國(guó)DEM作為模板，即采用2.5′ × 2.5′ 作為格網(wǎng)單元，且與中國(guó)DEM的行列數(shù)、起始點(diǎn)坐標(biāo)相同，每個(gè)格網(wǎng)的大小一樣。這樣，插值時(shí)共有1 134×1 479=1 677 186個(gè)格網(wǎng)點(diǎn)。

1.3.43D TPS插值 根據(jù)已知?dú)庀笳军c(diǎn)的3D坐標(biāo)(x,y,z)和年均降水量值P，推求未知點(diǎn)——無(wú)觀測(cè)值的每個(gè)格網(wǎng)點(diǎn)處的年均降水量值，由于格網(wǎng)點(diǎn)間是均勻分布的，利用ArcGIS Engine的二次開(kāi)發(fā)技術(shù)不難求出每個(gè)格網(wǎng)點(diǎn)的3D坐標(biāo)(x′,y′,z′)，再根據(jù)公式(2)可以求出每個(gè)格網(wǎng)點(diǎn)處的年均降水量P′。

1.3.5生成中國(guó)年降水量空間分布圖 可用中國(guó)DEM作為模板，用年均降水量值P′，替換中國(guó)DEM中每個(gè)格網(wǎng)點(diǎn)處的高程值，就可以得到中國(guó)降水量空間分布圖。

1.3.6與其他方法的對(duì)比分析 與二維薄板樣條法(簡(jiǎn)稱2D TPS)、反距離加權(quán)法(Inverse distance weighted，簡(jiǎn)稱IDW)、樣條法(Spline)、普通克里金、泛克里金進(jìn)行比較分析[15-16]。

本研究采用標(biāo)準(zhǔn)誤差(SD)作為評(píng)價(jià)指標(biāo)，即插值平均誤差平方和的算術(shù)平方根。設(shè)n個(gè)估計(jì)值與真實(shí)值的誤差為ε1、ε2、……、εn，則標(biāo)準(zhǔn)誤差公式為：

(8)

標(biāo)準(zhǔn)誤差的值越小，插值精度越高；反之越低。

2 結(jié)果與分析

2.1試驗(yàn)結(jié)果 圖2是利用3D TPS插值的結(jié)果，由于插值過(guò)程中充分考慮高程值對(duì)降水量的影響，插值效果較好。然而，由于本研究區(qū)域較廣，降水量變化很大，從東南到西北，降水量總體呈遞減趨勢(shì)。東南沿海區(qū)域降水量多且地勢(shì)較低，西北降水少。對(duì)于某些局部，如西藏南部降水量明顯增多區(qū)域，圖2中清晰可鑒；而在圖3和圖4中不很明顯。另外，圖2右下角的南海諸島，由于缺乏氣象站點(diǎn)，且島礁相對(duì)大陸面積較小，結(jié)果不是很清楚，圖中環(huán)形區(qū)域僅為示意圖。

圖3是利用2D TPS插值的結(jié)果，由于2D TPS僅考慮了地理坐標(biāo)(x,y)，沒(méi)有高程值參與插值，在地理位置相近區(qū)域，降水量相差不多。因此，許多插值細(xì)節(jié)不能體現(xiàn)出來(lái)。

圖4是利用普通克里金插值的結(jié)果，該方法僅考慮地理位置、沒(méi)有考慮高程值對(duì)降水量的影響，因此，除了與2D TPS方法有類似的問(wèn)題外，普通克里金插值還出現(xiàn)了一些斑塊狀區(qū)域以及亮點(diǎn)，在局部區(qū)域形成個(gè)別“牛眼”現(xiàn)象，說(shuō)明與周?chē)鷧^(qū)域顏色有較大差異。這可能由于采用的氣象站數(shù)量有限且分布不均勻，同時(shí)沒(méi)有考慮高程值影響造成的。

圖2 基于3D TPS插值的中國(guó)年均降水量空間分布圖Fig.2 Spatial distribution map of annual precipitation in China based on 3D TPS interpolation

圖3 基于2D TPS插值的中國(guó)年均降水量空間分布圖Fig.3 Spatial distribution map of annual precipitation in China based on 2D TPS interpolation

圖4 基于普通克里金插值的中國(guó)年均降水量空間分布圖Fig.4 Spatial distribution map of annual precipitation in China based on ordinary Kriging interpolation

2.2精度評(píng)價(jià) 為了評(píng)價(jià)3D TPS插值方法的精度，進(jìn)一步進(jìn)行各種插值方法的對(duì)比分析。除了TPS外，還有反距離加權(quán)插值法(IDW)、樣條法，克里金插值再細(xì)分為普通克里金和泛克里金兩種插值方法。從454個(gè)氣象站點(diǎn)中，將20個(gè)氣象站點(diǎn)觀測(cè)值作為校驗(yàn)值，為了提高可信度，這20個(gè)站點(diǎn)是隨機(jī)選取的，且盡量覆蓋整個(gè)研究區(qū)域；這20個(gè)站點(diǎn)不參與插值，僅用來(lái)驗(yàn)證插值結(jié)果。

3D TPS的SD值最小(SD=113.95)，插值精度最高，2D TPS次之(SD=156.82)，再次是泛克里金(SD=160.97)，與2D TPS相差不多，其他依次是IDW(SD=181.56)和普通克里金(SD=197.19)，精度最低的是樣條法(SD=241.49)(表1)。由于20個(gè)氣象站點(diǎn)的選取是隨機(jī)的，且不參與插值，僅用于驗(yàn)證，因此，結(jié)果的可信度還是較高的。

進(jìn)一步分析得知，絕對(duì)誤差最大的分別是元江氣象站和五臺(tái)山氣象站，如，元江氣象站的絕對(duì)誤差最小的是由3D TPS插值得到390.81；而絕對(duì)誤差最大的是樣條法，高達(dá)707.94。五臺(tái)山氣象站除3D TPS插值得到的絕對(duì)誤差55.67外，其他方法插值得到絕對(duì)誤差都在300左右。這兩個(gè)氣象站共同特點(diǎn)是分別與鄰域氣象站點(diǎn)的高程值有較大差異，如：元江氣象站高程為401 m，而周?chē)鷼庀笳军c(diǎn)高程都在1 000 m以上，大大低于周?chē)鷼庀笳镜母叱讨?；五臺(tái)山氣象站的高程為2 208 m，而周?chē)鷼庀笳镜母叱讨党贁?shù)略高于1 000 m(但在1 500 m以下)外，大多低于1 000 m，遠(yuǎn)高于周?chē)鷼庀笳靖叱讨?。其他絕對(duì)誤差較大的氣象站，如，安慶氣象站等，情況類似。由于3D TPS插值充分考慮了高程值的影響，在地形起伏較大、高程變化明顯的區(qū)域，插值精度最高，優(yōu)勢(shì)較為明顯。

表1 各種插值方法的比較Table 1 The comparison of different interpolation methods

最后，將本研究與其它文獻(xiàn)研究進(jìn)行比較。文獻(xiàn)[17]應(yīng)用趨勢(shì)面模型對(duì)我國(guó)34個(gè)重要城市模擬，得到SD=155.95，僅大于本研究3D TPS模型的SD，即文獻(xiàn)[17]應(yīng)用的趨勢(shì)面模型插值精度低于3D TPS，而高于本研究中其它5種插值方法。由于文獻(xiàn)[17]也考慮了高程值的影響，所以會(huì)比本研究中其它5種方法(沒(méi)有考慮高程影響)插值精度較高。但是，在同樣考慮了高程值影響的情況下，3D TPS插值精度更高。

3 結(jié)束語(yǔ)

我國(guó)是一個(gè)多山的國(guó)家，地形變化復(fù)雜，草地資源分布廣闊。在降水量空間插值研究中，高程值往往被忽略，以致造成較大誤差，本研究采用3D TPS對(duì)降水量進(jìn)行空間插值，除了考慮地理位置外，還將高程值作為重要的插值因子，充分考慮了高程值對(duì)降水量的影響。通過(guò)試驗(yàn)，與其他插值方法進(jìn)行比較分析，得出3D TPS對(duì)降水量進(jìn)行空間插值的精度最高，尤其在高程值變化較大的區(qū)域優(yōu)勢(shì)更為明顯。這為森林、草地植被、水資源、生物量等空間分布的研究提供了重要的參考依據(jù)。

[1] 崔慶虎,蔣志剛,劉季科,等.青藏高原草地退化原因述評(píng)[J].草業(yè)科學(xué),2007,24(5):20-26.

[2] 郭姜寧.森林及草地植被與降水[J].中國(guó)草業(yè)科學(xué),1987,4(4):58-60.

[3] 劉桂霞,苗玉華,李記開(kāi).放牧干擾和年際間降水量變化對(duì)地木耳生長(zhǎng)速度的影響[J].草業(yè)科學(xué),2011,28(9):1649-1652.

[4] 張曉松.我國(guó)人均水資源為世界平均水平1/4[J].國(guó)土經(jīng)濟(jì),2001(1):46.

[5] 朱會(huì)義,賈紹鳳.降雨信息空間插值的不確定性分析[J].地理科學(xué)進(jìn)展,2004,23(2):34-42.

[6] 儲(chǔ)少林,周兆葉,袁雷,等.降水空間插值方法應(yīng)用研究——以甘肅省為例[J].草業(yè)科學(xué),2008,25(6):19-23.

[7] 馬軒龍,李春娥,陳全功.基于GIS的氣象要素空間插值方法研究[J].草業(yè)科學(xué),2008,25(11):13-19.

[8] Duchon J.Splines minimizing rotation-invariant semi-norms in Sobolev spaces[J].Constructive Theory of Functions of Several Variables,1977,571:85-100.

[9] Hutchinson,M.F.Interpolating mean rainfall using thin plate smoothing splines[J].International Journal of Geographic Information Systems,1995,9(4)：385-403.

[10] Wahba G.Spline Models for Observational Data[M].Philadelphia:Society for Industrial and Applied Mathematics,1990:10-99.

[11] 孫海燕,丁咚.薄板樣條函數(shù)及復(fù)雜曲面的數(shù)學(xué)表示[J].測(cè)繪工程,2006,15(2):7-8.

[12] 閻洪.薄板光順樣條插值與中國(guó)氣候空間模擬[J].地理科學(xué),2004,24(2):163-169.

[13] Sinthanayothin C,Bholsithi W,Tharanon W.Tooth alignment of the dental cast using 3D thin plate spline[A].Proceedings of the 3rd IASTED International Conference on Advances in Computer Science and Technology[C].Anaheim,CA,USA:ACTA Press,2007:255-259.

[14] 游松財(cái),李軍.海拔誤差影響氣溫空間插值誤差的研究[J].自然資源學(xué)報(bào),2005,20(1):140-144.

[15] 徐吉宏,柳小妮,張德罡,等.基于ArcGIS的天祝草地綜合順序分類研究[J].草業(yè)科學(xué),2011,28(5):866-870.

[16] 朱猛蒙,孫玉榮,張蓉,等.基于GIS的苜蓿斑蚜區(qū)域化預(yù)測(cè)預(yù)報(bào)技術(shù)初步研究[J].草業(yè)學(xué)報(bào),2011,20(2):163-169.

[17] 尚宗波,高瓊,楊奠安.利用中國(guó)氣候信息系統(tǒng)研究年降水量空間分布規(guī)律[J].生態(tài)學(xué)報(bào),2001,21(5):689-694.