鄭興文,胡泓達(dá),許家琦,舒紅

(武漢大學(xué)測(cè)繪遙感信息工程國(guó)家重點(diǎn)實(shí)驗(yàn)室,湖北武漢 430079)

Kriging插值組件的混合編程實(shí)現(xiàn)

鄭興文?,胡泓達(dá),許家琦,舒紅

(武漢大學(xué)測(cè)繪遙感信息工程國(guó)家重點(diǎn)實(shí)驗(yàn)室,湖北武漢 430079)

在.NET技術(shù)框架下，提出使用C#語言和R語言混合編程模式實(shí)現(xiàn)普通Kriging插值功能組件，并應(yīng)用該組件完成東北三省氣溫插值分析。實(shí)驗(yàn)效果表明，在.Net技術(shù)框架下，使用R和C#混合編程模式開發(fā)普通Kriging插值功能的周期短成本低且容易擴(kuò)展到其他空間統(tǒng)計(jì)分析功能實(shí)現(xiàn)。

R語言；.NET組件技術(shù)；普通Kriging；混合編程

1 前 言

Kriging插值是空間統(tǒng)計(jì)分析的重要方法之一,被廣泛用于地質(zhì)礦產(chǎn)[1]、氣象[2]、地理信息[3]、計(jì)算機(jī)圖形學(xué)[4]等方向。常見的Kriging插值模型包括普通Kriging模型、簡(jiǎn)單Kriging模型、塊Kriging模型、指示Kriging模型、協(xié)同Kriging模型等[5]。當(dāng)區(qū)域化變量Z(x)的數(shù)學(xué)期望E[Z(x)]=m為未知常數(shù)時(shí),常采用普通Kriging插值模型[6]進(jìn)行插值。目前僅有少量軟件提供普通Kriging插值功能,例如ArcGIS、MapGIS等。這些價(jià)格昂貴的商業(yè)軟件中Kriging插值實(shí)現(xiàn)方法一般人無法獲知,不利于Kriging擴(kuò)展的研發(fā)。而本文采用C#與開源語言R混合編程模式開發(fā)Kriging插值等空間統(tǒng)計(jì)分析功能,不僅成本低,而且更利于體現(xiàn)插值的原理并方便進(jìn)行擴(kuò)展實(shí)現(xiàn)。運(yùn)用R語言實(shí)現(xiàn)Kriging插值功能具有諸多優(yōu)勢(shì)[7],R語言開源且擁有多種空間統(tǒng)計(jì)分析功能函數(shù)包可直接調(diào)用,使開發(fā)周期大大縮短。但R語言存在操作界面不友好、功能的跨平臺(tái)跨語言的可重用性不強(qiáng)等缺點(diǎn),利用.NET技術(shù)實(shí)現(xiàn)R與C#混合編程可以彌補(bǔ)這些缺點(diǎn),實(shí)現(xiàn)優(yōu)勢(shì)互補(bǔ)[8],便于復(fù)雜空間統(tǒng)計(jì)分析功能的開發(fā)。

2 混合開發(fā)模式及Kriging插值原理

2.1 R語言與.NET技術(shù)

R是一個(gè)開放的統(tǒng)計(jì)編程環(huán)境,是一種語言,是S語言的一種實(shí)現(xiàn)。S語言是由AT＆T Bell實(shí)驗(yàn)室的Rick Becker,John Chambers和Allan Wilks開發(fā)的一種進(jìn)行數(shù)據(jù)探索、統(tǒng)計(jì)分析、作圖的解釋型語言[9]。

R是完全開源的,且具有強(qiáng)大的圖形處理和統(tǒng)計(jì)分析功能,R軟件開發(fā)小組不斷提供大量的功能庫函數(shù),這些優(yōu)勢(shì)促使R軟件受到廣大科研工作者的喜愛,本文選擇R進(jìn)行底層普通Kriging插值編程實(shí)現(xiàn)正是基于R的這些優(yōu)點(diǎn)。

.NET技術(shù)是Microsoft公司繼.COM技術(shù)后新推出的組件技術(shù)。.NET組件技術(shù)較原有的.COM組件技術(shù)大大改善[10]?；?NET組件編寫的軟件可以大大提高其復(fù)用性和開發(fā)效率。.NET核心部分NET Framework 4.0提供了CLR(公共語言運(yùn)行庫),在CLR控制下運(yùn)行的代碼為托管代碼(managed code)[11],通過Visual Studio.NET開發(fā)工具在公共語言規(guī)范(CLS)下編譯中間語言,生成動(dòng)態(tài)鏈接庫DLL提供組件接口。其他高級(jí)語言可以引用此DLL并調(diào)用接口函數(shù)實(shí)現(xiàn)功能調(diào)用。R函數(shù)被C#調(diào)用的實(shí)現(xiàn)代碼就是在DLL構(gòu)建的過程中完成的。

2.2 C#調(diào)用R函數(shù)方法

R.NET是以.NET框架與R統(tǒng)計(jì)相結(jié)合的一種混合編程技術(shù)。它是R到.NET Framework的中間轉(zhuǎn)換類型庫,作為一個(gè)橋梁連接起R語言和.NET技術(shù)。具體的語言類型轉(zhuǎn)換見網(wǎng)站:http://rdotnet.codeplex. com/documentation。R函數(shù)被C#調(diào)用必須首先初始化R環(huán)境,然后把R編寫的功能函數(shù)的變量進(jìn)行類型轉(zhuǎn)換并調(diào)用提供函數(shù)接口。具體操作是首先用Microsoft visual studio 2010創(chuàng)建類CSTlib和接口ISTlib,并添加引用動(dòng)態(tài)庫RDotNet.dll,對(duì)注冊(cè)表進(jìn)行相應(yīng)的設(shè)置并創(chuàng)建密鑰。R環(huán)境在C#中初始化的代碼如下。

‘聲明R引擎變量

static private REngine engine=null；

‘選擇R文件夾路徑，設(shè)環(huán)境變量

System.Environment.SetEnvironmentVariable("PATH"，＠" C：Program FilesRR-3.0.1ini386")；

‘變量賦R引擎實(shí)例

engine=REngine.CreateInstance("RDotNet")；

‘R初始化

engine.Initialize()；

2.3 Kriging插值原理

(1)變異函數(shù)

由美國(guó)地理學(xué)家W.R.Tobler在1970年提出的地理學(xué)第一定律[12]認(rèn)為:任何事物都相關(guān),只是相近的事物關(guān)聯(lián)更相關(guān)。變異函數(shù)是描述區(qū)域化變量的空間結(jié)構(gòu)和隨機(jī)性變化的基本工具。可以通過變異函數(shù)曲線分析變量隨距離變化的統(tǒng)計(jì)相關(guān)性程度,從而反映數(shù)據(jù)整體空間變化特征。其變異函數(shù)的計(jì)算公式為:

式中r(h)為變異函數(shù)；h為樣點(diǎn)空間間隔距離,稱為步長(zhǎng)；N(h)為間隔距離為h的樣點(diǎn)數(shù)；Z(xi)和Z(xi+h)分別是區(qū)域變量Z(x)在空間位置xi和xi+h的實(shí)測(cè)值。常用的基本變異函數(shù)模型有球形模型、指數(shù)模型、高斯模型、純塊金模型等,具體選擇何種模型要根據(jù)變異函數(shù)曲線擬合效果而定,一般采用交叉驗(yàn)證法[11]進(jìn)行判斷變異函數(shù)模型選擇的有效程度。普通Kriging插值要求數(shù)據(jù)要符合固有假設(shè):

即數(shù)據(jù)Z(x)均值存在且并不取決于x；對(duì)于任何距離h變量[Z(x+h)-Z(x)]具有一個(gè)有限的方差,此方差不取決于x只與距離h有關(guān)。首先對(duì)原始數(shù)據(jù)進(jìn)行預(yù)處理,針對(duì)符合固有假設(shè)的數(shù)據(jù)行插值會(huì)獲得較好的插值精度。

(2)普通Kriging插值

一個(gè)隨機(jī)變量在空間上往往觀測(cè)值是離散的,在未進(jìn)行觀測(cè)的空間位置我們往往采用插值處理。普通Kriging[5]就是插值中最為普遍使用的一種方法,Kriging插值法是一種有效的線性無偏估計(jì)方法,估計(jì)誤差方差最小是Kriging插值法的顯著特點(diǎn)。普通Kriging的估計(jì)公式為:

Z?(x0)為x0處屬性的預(yù)測(cè)值,λi為權(quán)重矩陣, Z(xi)為已知空間點(diǎn)上的屬性觀測(cè)值。根據(jù)無偏性和估值方差最小可推導(dǎo)出權(quán)重矩陣計(jì)算公式如下:

上式中γij是將xi和xj的距離代入理論變異函數(shù)計(jì)算公式γ(h)得到的,μ為拉格朗日乘數(shù),解算上式即可獲得權(quán)重矩陣,從而得到插值估值。

3 程序?qū)崿F(xiàn)過程

3.1 R進(jìn)行Kriging插值編程

要實(shí)現(xiàn)普通Kriging需要用R編寫數(shù)據(jù)讀取功能代碼,去趨勢(shì)功能代碼,變異函數(shù)擬合功能代碼,普通Kriging插值功能代碼,普通Kriging插值圖轉(zhuǎn)換為等高線圖功能代碼,繪制散點(diǎn)圖功能代碼,結(jié)果保存功能代碼等。普通Kriging插值功能代碼:

將經(jīng)驗(yàn)變異函數(shù)賦值給v

v〈-variogram(vgm.data～1，～longt+lat，data)

將擬合好的模型(變異函數(shù)表達(dá)式)賦給m(alt為變異函數(shù)擬合模型)

m〈-fit.variogram(v，vgm(as.numeric(偏基臺(tái)值)，alt，as.numeric(變程值)，as.numeric(塊金值)))

克里金插值結(jié)果賦值給lzn.kr(krige.data是去趨勢(shì)后的屬性數(shù)據(jù)；data.grid為數(shù)據(jù)框形式存儲(chǔ)的柵格數(shù)據(jù))

lzn.kr〈-krige(formula=krige.data～1，data，data.grid，model =m)

將初始插值結(jié)果寫成數(shù)據(jù)框賦值給krdata

krdata〈-as.data.frame(lzn.kr)

經(jīng)過趨勢(shì)補(bǔ)回后的插值結(jié)果賦值給lzn.kr$var1.pred lzn.kr$var1.pred〈-pred1誤差方差寫成標(biāo)準(zhǔn)差形式

lzn.kr$var1.var〈-sqrt(lzn.kr$var1.var)……

繪制預(yù)測(cè)值圖

pl1〈-spplot(lzn.kr["var1.pred"]，main="The result of ordinary Kriging"，col.regions=topo.colors)

繪制標(biāo)準(zhǔn)誤差圖

pl2〈-spplot(lzn.kr["var1.var"]，main="The standard error of ordinary Kriging")

3.2 .NET組件編寫

依據(jù)R.NET技術(shù)規(guī)范,編寫普通Kriging功能函數(shù)接口,設(shè)計(jì)函數(shù)屬性并生成秘鑰,在.NET技術(shù)框架下完成組件編寫。代碼如下:

public class CSTlib：ISTlib

{

設(shè)計(jì)普通Kriging功能函數(shù)接口

public void okrige(string alt，int n，……參數(shù))

{

運(yùn)用R.NET將C#變量類型轉(zhuǎn)換為R識(shí)別變量類型

CharacterVector Alt=engine.CreateCharacterVector(new string []{alt})；

engine.SetSymbol("alt"，Alt)；

調(diào)用底層R腳本功能函數(shù)

engine.Evaluate(＠"source('底層R函數(shù)腳本路徑')；okrige

(alt，n，……參數(shù))；")；

}

}

經(jīng)過.NET組件功能函數(shù)的編寫,將普通Kriging插值功能封裝在類CSTlib中并提供組件接口ISTlib,組件內(nèi)部不被用戶所知,用戶根據(jù)組件接口定義形式使用接口函數(shù)。

3.3 應(yīng)用程序調(diào)用組件接口函數(shù)

運(yùn)用C#編寫Windows窗體應(yīng)用程序,調(diào)用編寫好的.NET組件,可視化的展現(xiàn)整個(gè)普通克里金插值的過程,核心代碼如下:

組件接口函數(shù)實(shí)例化

Rlib.CSTlib sta=new Rlib.CSTlib()；

動(dòng)態(tài)庫接口函數(shù)調(diào)用

sta.okrige(model，order，x，y，column，psill，range，nugget，filename)；

用戶首先添加引用組件動(dòng)態(tài)鏈接庫DLL,對(duì)接口函數(shù)實(shí)例化,調(diào)用組件內(nèi)所需接口函數(shù)。根據(jù)用戶需求導(dǎo)入?yún)?shù),實(shí)現(xiàn)對(duì)組件接口函數(shù)的調(diào)用。

4 Kriging組件用于東北三省氣溫插值分析

4.1 研究區(qū)域概況及數(shù)據(jù)分析

實(shí)驗(yàn)數(shù)據(jù)從中國(guó)氣象科學(xué)數(shù)據(jù)共享服務(wù)網(wǎng)中下載獲取,整理獲得東北三省1970年～2011年年平均氣溫站點(diǎn)數(shù)據(jù)(共80個(gè)站點(diǎn)),站點(diǎn)分布如圖1所示。站點(diǎn)坐標(biāo)數(shù)據(jù)采用北京1954坐標(biāo)系。東北三省的背景柵格地圖(如圖2所示)是下載的modis數(shù)據(jù)經(jīng)高斯-克呂格投影、矢量邊界裁剪并重采樣獲得的柵格單元為8 000 m×8 000 m柵格圖。軟件通過導(dǎo)入已知站點(diǎn)數(shù)據(jù)和背景柵格地圖進(jìn)行普通Kriging插值。通過插值結(jié)果可以了解東北三省年平均氣溫的分布情況,從而可進(jìn)一步分析東北三省氣象特征。

圖1 東北三省氣溫觀測(cè)站點(diǎn)分布

圖2 東北三省柵格圖

4.2 插值過程及成果顯示

利用C#設(shè)計(jì)用戶交互界面(如圖3所示),界面包括已知站點(diǎn)數(shù)據(jù)導(dǎo)入、插值的背景柵格地圖導(dǎo)入、去趨勢(shì)統(tǒng)計(jì)分析、變異函數(shù)擬合、交叉驗(yàn)證、普通Kriging插值、插值結(jié)果保存、繪制等高線和散點(diǎn)圖等功能。

圖3 普通Kriging插值界面設(shè)計(jì)

首先導(dǎo)入東北三省站點(diǎn)坐標(biāo)以及氣溫?cái)?shù)據(jù),并導(dǎo)入插值的背景柵格地圖(如圖2所示),然后對(duì)原始站點(diǎn)氣象數(shù)據(jù)進(jìn)行趨勢(shì)分析(如圖4所示),分析獲得經(jīng)二次項(xiàng)去趨勢(shì)后的R2=93.95%說明趨勢(shì)面反映了原始數(shù)據(jù)中94%的變異性,還有6%的變化未能在趨勢(shì)面中得到表現(xiàn)而成為剩余；p〈0.05表明回歸方程通過了F檢驗(yàn),是顯著的。

圖4 去趨勢(shì)統(tǒng)計(jì)分析圖

根據(jù)交叉驗(yàn)證結(jié)果(如圖5(A)所示)選擇最佳變異函數(shù)模型為高斯模型,變異函數(shù)擬合曲線(如圖5 (B)所示)顯示了塊金值、變程和基臺(tái)值。

圖5 交叉驗(yàn)證結(jié)果(A)和變異函數(shù)擬合曲線(B)

進(jìn)行普通Kriging插值得到東北三省氣溫插值預(yù)測(cè)圖(圖6(A)所示)和誤差標(biāo)準(zhǔn)差圖(圖6(B)所示)。從圖中可以直觀看到東北三省1970年～2011年年均值氣溫總體分布呈現(xiàn)狀況:氣溫分布南高北低,與緯度越低氣溫越高(不考慮海洋、高山等地表形態(tài))基本吻合,由于不考慮地表形態(tài)等其他因素對(duì)溫度的影響,因此插值溫度分布曲線平滑,離站點(diǎn)越遠(yuǎn),插值誤差越大,插值精度在站點(diǎn)密集分布合理區(qū)域較高,在周圍站點(diǎn)稀疏區(qū)域精度較低。

圖6 東北三省氣溫普通Kriging插值圖(A)和誤差標(biāo)準(zhǔn)差圖(B)

插值結(jié)果圖可以被保存為shp格式便于ArcGIS調(diào)用分析,也可以轉(zhuǎn)換為等溫線圖和三維散點(diǎn)圖(如圖7所示)進(jìn)行顯示。

圖7 普通Kriging插值后轉(zhuǎn)換成三維散點(diǎn)圖

(圖7中x，y表示坐標(biāo)位置，單位為km，是通過WGS84坐標(biāo)經(jīng)高斯-克呂格法以3°帶東經(jīng)111°為中心投影獲得。SMTem軸表示溫度值，單位為0.1℃。散點(diǎn)圖直觀表示溫度在東北三省的分布情況。)

5 結(jié) 論

本文給出了基于.NET技術(shù)與R語言混合編程的普通Kriging插值實(shí)現(xiàn)過程。通過對(duì)東北三省氣溫的Kriging插值顯示分析,表明在結(jié)合R統(tǒng)計(jì)語言和.NET組件技術(shù)進(jìn)行地理時(shí)空統(tǒng)計(jì)組件功能軟件開發(fā)不僅周期短,成本低,而且跨平臺(tái)跨語言系統(tǒng)開發(fā)伸展性強(qiáng)。

普通Kriging是地學(xué)數(shù)據(jù)統(tǒng)計(jì)分析的重要工具,我們?cè)陉U述軟件開發(fā)和程序?qū)崿F(xiàn)過程中也給出了普通Kriging插值模型和算法。此編程開發(fā)技術(shù)方法可進(jìn)一步推廣對(duì)其他空間統(tǒng)計(jì)分析功能的開發(fā)應(yīng)用。

[1] DAVISJC.Statistics and Data Analysis in Geology[M].New York John Wiley＆Sons，2002，57～61.

[2] 李莎，舒紅，徐正全.利用時(shí)空Kriging進(jìn)行氣溫插值研究[J].武漢大學(xué)學(xué)報(bào)·信息科學(xué)版，2012，37(2).

[3] 劉志平，何秀鳳，張淑輝.多測(cè)度加權(quán)克里金法在高邊坡變形穩(wěn)定性分析中的應(yīng)用[J].水利學(xué)報(bào)，2009，40(6).

[4] 李小斌，錢建生，趙志凱.基于克里金插值的腦電圖成像系統(tǒng)[J].計(jì)算機(jī)應(yīng)用于軟件，2010，27(7).

[5] 張仁鐸.空間變異理論及應(yīng)用[M].北京：科學(xué)出版社，2005：27～35.

[6] 劉湘南，黃方，王平.GIS空間分析原理與方法[M].北京：科學(xué)出版社，2008：200～201.

[7] 楊中慶.基于R語言的空間統(tǒng)計(jì)分析研究與應(yīng)用[D].廣州：暨南大學(xué)，2006：17～18.

[8] 趙毅，史權(quán)，趙鎖奇等.R語言與.NET混合編程在重質(zhì)油數(shù)據(jù)管理分析中的應(yīng)用[J].計(jì)算機(jī)與應(yīng)用化學(xué)，2012，29 (4).

[9] 薛毅，陳立萍.統(tǒng)計(jì)建模與R軟件[M].北京：清華大學(xué)出版社，2007.

[10] 李志毅，趙政.軟件復(fù)用與COM及.NET組件技術(shù)[J].微處理機(jī)，2006(6).

[11] 覃釗.基于.NET的MATLAB與Visual Basic混合編程的研究[J].城市勘測(cè)，2012(6).

[12] Waldo R.Tobler.A computer movie simulating urban growth in the Detroit region[J].Economic Geography，1970，46(2)：243～240.

The M ixed Programm ing com ponent of Ordinary K riging

Zheng Xingwen,Hu Hongda,Xu Jiaqi,Shu Hong
(State Key Laboratory of Information Engineering in Surveying,Mapping and Remote Sensing,Wuhan University,Wuhan 430079,China)

The paper proposes amix-programmingmodelwith C#and R within the.NET technology framework and analyzes the result of the Ordinary Kriging with the temperature data of the northeastern China.Using thismodel we can obtain the software quickly and inexpensively.What’smore,we can add more new functions into the software.

R；.NET component technology；ordinary kriging；mixed programming

1672-8262(2013)06-139-05

P209

A

2013—06—22

鄭興文(1988—),男,碩士研究生,研究方向:時(shí)空統(tǒng)計(jì)分析。

國(guó)家高技術(shù)研究發(fā)展計(jì)劃(863計(jì)劃)項(xiàng)目(2011AA0105002)；國(guó)家自然科學(xué)基金項(xiàng)目(41171313)；國(guó)家高技術(shù)研究發(fā)展計(jì)劃項(xiàng)目(2012AA1214002)