Kriging插值组件的混合编程实现

郑兴文,胡泓达,许家琦,舒红

(武汉大学测绘遥感信息工程国家重点实验室,湖北武汉 430079)

Kriging插值组件的混合编程实现

郑兴文∗,胡泓达,许家琦,舒红

(武汉大学测绘遥感信息工程国家重点实验室,湖北武汉 430079)

在.NET技术框架下，提出使用C#语言和R语言混合编程模式实现普通Kriging插值功能组件，并应用该组件完成东北三省气温插值分析。实验效果表明，在.Net技术框架下，使用R和C#混合编程模式开发普通Kriging插值功能的周期短成本低且容易扩展到其他空间统计分析功能实现。

R语言；.NET组件技术；普通Kriging；混合编程

1 前 言

Kriging插值是空间统计分析的重要方法之一,被广泛用于地质矿产[1]、气象[2]、地理信息[3]、计算机图形学[4]等方向。常见的Kriging插值模型包括普通Kriging模型、简单Kriging模型、块Kriging模型、指示Kriging模型、协同Kriging模型等[5]。当区域化变量Z(x)的数学期望E[Z(x)]=m为未知常数时,常采用普通Kriging插值模型[6]进行插值。目前仅有少量软件提供普通Kriging插值功能,例如ArcGIS、MapGIS等。这些价格昂贵的商业软件中Kriging插值实现方法一般人无法获知,不利于Kriging扩展的研发。而本文采用C#与开源语言R混合编程模式开发Kriging插值等空间统计分析功能,不仅成本低,而且更利于体现插值的原理并方便进行扩展实现。运用R语言实现Kriging插值功能具有诸多优势[7],R语言开源且拥有多种空间统计分析功能函数包可直接调用,使开发周期大大缩短。但R语言存在操作界面不友好、功能的跨平台跨语言的可重用性不强等缺点,利用.NET技术实现R与C#混合编程可以弥补这些缺点,实现优势互补[8],便于复杂空间统计分析功能的开发。

2 混合开发模式及Kriging插值原理

2.1 R语言与.NET技术

R是一个开放的统计编程环境,是一种语言,是S语言的一种实现。S语言是由AT＆T Bell实验室的Rick Becker,John Chambers和Allan Wilks开发的一种进行数据探索、统计分析、作图的解释型语言[9]。

R是完全开源的,且具有强大的图形处理和统计分析功能,R软件开发小组不断提供大量的功能库函数,这些优势促使R软件受到广大科研工作者的喜爱,本文选择R进行底层普通Kriging插值编程实现正是基于R的这些优点。

.NET技术是Microsoft公司继.COM技术后新推出的组件技术。.NET组件技术较原有的.COM组件技术大大改善[10]。基于.NET组件编写的软件可以大大提高其复用性和开发效率。.NET核心部分NET Framework 4.0提供了CLR(公共语言运行库),在CLR控制下运行的代码为托管代码(managed code)[11],通过Visual Studio.NET开发工具在公共语言规范(CLS)下编译中间语言,生成动态链接库DLL提供组件接口。其他高级语言可以引用此DLL并调用接口函数实现功能调用。R函数被C#调用的实现代码就是在DLL构建的过程中完成的。

2.2 C#调用R函数方法

R.NET是以.NET框架与R统计相结合的一种混合编程技术。它是R到.NET Framework的中间转换类型库,作为一个桥梁连接起R语言和.NET技术。具体的语言类型转换见网站:http://rdotnet.codeplex. com/documentation。R函数被C#调用必须首先初始化R环境,然后把R编写的功能函数的变量进行类型转换并调用提供函数接口。具体操作是首先用Microsoft visual studio 2010创建类CSTlib和接口ISTlib,并添加引用动态库RDotNet.dll,对注册表进行相应的设置并创建密钥。R环境在C#中初始化的代码如下。

‘声明R引擎变量

static private REngine engine=null；

‘选择R文件夹路径，设环境变量

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

‘变量赋R引擎实例

engine=REngine.CreateInstance("RDotNet")；

‘R初始化

engine.Initialize()；

2.3 Kriging插值原理

(1)变异函数

由美国地理学家W.R.Tobler在1970年提出的地理学第一定律[12]认为:任何事物都相关,只是相近的事物关联更相关。变异函数是描述区域化变量的空间结构和随机性变化的基本工具。可以通过变异函数曲线分析变量随距离变化的统计相关性程度,从而反映数据整体空间变化特征。其变异函数的计算公式为:

式中r(h)为变异函数；h为样点空间间隔距离,称为步长；N(h)为间隔距离为h的样点数；Z(xi)和Z(xi+h)分别是区域变量Z(x)在空间位置xi和xi+h的实测值。常用的基本变异函数模型有球形模型、指数模型、高斯模型、纯块金模型等,具体选择何种模型要根据变异函数曲线拟合效果而定,一般采用交叉验证法[11]进行判断变异函数模型选择的有效程度。普通Kriging插值要求数据要符合固有假设:

即数据Z(x)均值存在且并不取决于x；对于任何距离h变量[Z(x+h)-Z(x)]具有一个有限的方差,此方差不取决于x只与距离h有关。首先对原始数据进行预处理,针对符合固有假设的数据行插值会获得较好的插值精度。

(2)普通Kriging插值

一个随机变量在空间上往往观测值是离散的,在未进行观测的空间位置我们往往采用插值处理。普通Kriging[5]就是插值中最为普遍使用的一种方法,Kriging插值法是一种有效的线性无偏估计方法,估计误差方差最小是Kriging插值法的显著特点。普通Kriging的估计公式为:

Z∗(x0)为x0处属性的预测值,λi为权重矩阵, Z(xi)为已知空间点上的属性观测值。根据无偏性和估值方差最小可推导出权重矩阵计算公式如下:

上式中γij是将xi和xj的距离代入理论变异函数计算公式γ(h)得到的,μ为拉格朗日乘数,解算上式即可获得权重矩阵,从而得到插值估值。

3 程序实现过程

3.1 R进行Kriging插值编程

要实现普通Kriging需要用R编写数据读取功能代码,去趋势功能代码,变异函数拟合功能代码,普通Kriging插值功能代码,普通Kriging插值图转换为等高线图功能代码,绘制散点图功能代码,结果保存功能代码等。普通Kriging插值功能代码:

将经验变异函数赋值给v

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

将拟合好的模型(变异函数表达式)赋给m(alt为变异函数拟合模型)

m〈-fit.variogram(v，vgm(as.numeric(偏基台值)，alt，as.numeric(变程值)，as.numeric(块金值)))

克里金插值结果赋值给lzn.kr(krige.data是去趋势后的属性数据；data.grid为数据框形式存储的栅格数据)

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

将初始插值结果写成数据框赋值给krdata

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

经过趋势补回后的插值结果赋值给lzn.kr$var1.pred lzn.kr$var1.pred〈-pred1误差方差写成标准差形式

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

绘制预测值图

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

绘制标准误差图

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

3.2 .NET组件编写

依据R.NET技术规范,编写普通Kriging功能函数接口,设计函数属性并生成秘钥,在.NET技术框架下完成组件编写。代码如下:

public class CSTlib：ISTlib

{

设计普通Kriging功能函数接口

public void okrige(string alt，int n，……参数)

{

运用R.NET将C#变量类型转换为R识别变量类型

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

engine.SetSymbol("alt"，Alt)；

调用底层R脚本功能函数

engine.Evaluate(＠"source('底层R函数脚本路径')；okrige

(alt，n，……参数)；")；

}

}

经过.NET组件功能函数的编写,将普通Kriging插值功能封装在类CSTlib中并提供组件接口ISTlib,组件内部不被用户所知,用户根据组件接口定义形式使用接口函数。

3.3 应用程序调用组件接口函数

运用C#编写Windows窗体应用程序,调用编写好的.NET组件,可视化的展现整个普通克里金插值的过程,核心代码如下:

组件接口函数实例化

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

动态库接口函数调用

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

用户首先添加引用组件动态链接库DLL,对接口函数实例化,调用组件内所需接口函数。根据用户需求导入参数,实现对组件接口函数的调用。

4 Kriging组件用于东北三省气温插值分析

4.1 研究区域概况及数据分析

实验数据从中国气象科学数据共享服务网中下载获取,整理获得东北三省1970年～2011年年平均气温站点数据(共80个站点),站点分布如图1所示。站点坐标数据采用北京1954坐标系。东北三省的背景栅格地图(如图2所示)是下载的modis数据经高斯-克吕格投影、矢量边界裁剪并重采样获得的栅格单元为8 000 m×8 000 m栅格图。软件通过导入已知站点数据和背景栅格地图进行普通Kriging插值。通过插值结果可以了解东北三省年平均气温的分布情况,从而可进一步分析东北三省气象特征。

图1 东北三省气温观测站点分布

图2 东北三省栅格图

4.2 插值过程及成果显示

利用C#设计用户交互界面(如图3所示),界面包括已知站点数据导入、插值的背景栅格地图导入、去趋势统计分析、变异函数拟合、交叉验证、普通Kriging插值、插值结果保存、绘制等高线和散点图等功能。

图3 普通Kriging插值界面设计

首先导入东北三省站点坐标以及气温数据,并导入插值的背景栅格地图(如图2所示),然后对原始站点气象数据进行趋势分析(如图4所示),分析获得经二次项去趋势后的R2=93.95%说明趋势面反映了原始数据中94%的变异性,还有6%的变化未能在趋势面中得到表现而成为剩余；p〈0.05表明回归方程通过了F检验,是显著的。

图4 去趋势统计分析图

根据交叉验证结果(如图5(A)所示)选择最佳变异函数模型为高斯模型,变异函数拟合曲线(如图5 (B)所示)显示了块金值、变程和基台值。

图5 交叉验证结果(A)和变异函数拟合曲线(B)

进行普通Kriging插值得到东北三省气温插值预测图(图6(A)所示)和误差标准差图(图6(B)所示)。从图中可以直观看到东北三省1970年～2011年年均值气温总体分布呈现状况:气温分布南高北低,与纬度越低气温越高(不考虑海洋、高山等地表形态)基本吻合,由于不考虑地表形态等其他因素对温度的影响,因此插值温度分布曲线平滑,离站点越远,插值误差越大,插值精度在站点密集分布合理区域较高,在周围站点稀疏区域精度较低。

图6 东北三省气温普通Kriging插值图(A)和误差标准差图(B)

插值结果图可以被保存为shp格式便于ArcGIS调用分析,也可以转换为等温线图和三维散点图(如图7所示)进行显示。

图7 普通Kriging插值后转换成三维散点图

(图7中x，y表示坐标位置，单位为km，是通过WGS84坐标经高斯-克吕格法以3°带东经111°为中心投影获得。SMTem轴表示温度值，单位为0.1℃。散点图直观表示温度在东北三省的分布情况。)

5 结 论

本文给出了基于.NET技术与R语言混合编程的普通Kriging插值实现过程。通过对东北三省气温的Kriging插值显示分析,表明在结合R统计语言和.NET组件技术进行地理时空统计组件功能软件开发不仅周期短,成本低,而且跨平台跨语言系统开发伸展性强。

普通Kriging是地学数据统计分析的重要工具,我们在阐述软件开发和程序实现过程中也给出了普通Kriging插值模型和算法。此编程开发技术方法可进一步推广对其他空间统计分析功能的开发应用。

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

[2] 李莎，舒红，徐正全.利用时空Kriging进行气温插值研究[J].武汉大学学报·信息科学版，2012，37(2).

[3] 刘志平，何秀凤，张淑辉.多测度加权克里金法在高边坡变形稳定性分析中的应用[J].水利学报，2009，40(6).

[4] 李小斌，钱建生，赵志凯.基于克里金插值的脑电图成像系统[J].计算机应用于软件，2010，27(7).

[5] 张仁铎.空间变异理论及应用[M].北京：科学出版社，2005：27～35.

[6] 刘湘南，黄方，王平.GIS空间分析原理与方法[M].北京：科学出版社，2008：200～201.

[7] 杨中庆.基于R语言的空间统计分析研究与应用[D].广州：暨南大学，2006：17～18.

[8] 赵毅，史权，赵锁奇等.R语言与.NET混合编程在重质油数据管理分析中的应用[J].计算机与应用化学，2012，29 (4).

[9] 薛毅，陈立萍.统计建模与R软件[M].北京：清华大学出版社，2007.

[10] 李志毅，赵政.软件复用与COM及.NET组件技术[J].微处理机，2006(6).

[11] 覃钊.基于.NET的MATLAB与Visual Basic混合编程的研究[J].城市勘测，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—),男,硕士研究生,研究方向:时空统计分析。

国家高技术研究发展计划(863计划)项目(2011AA0105002)；国家自然科学基金项目(41171313)；国家高技术研究发展计划项目(2012AA1214002)