“地统计学方法”课程实验设计

杨 勇, 杨 雪, 刘诗诗

(华中农业大学 资源与环境学院, 湖北 武汉 430070)

地统计学(Geostatistics)是以区域化变量理论为基础,以变异函数为主要工具,研究那些在空间分布上既有随机性，又有结构性,或空间相关性和依赖性的自然现象的科学[1]。它起源于20世纪50年代的南非采矿业中,目前已被广泛应用于地质[2-3]、土壤[4-5]、生态环境[6-7]、气象[8-9]等涉及地理属性空间/时空变异及插值的多个领域,是集空间数据采集、统计分析、建模与计算、地图制图等多个技术环节理论与方法的集合。其核心技术包括基于空间样点的地理变量试验变异函数技术、理论变异函数拟合和地统计插值等,最常用的方法为普通克里格方法。此外,随着各个领域应用的延伸,各种衍生方法不断推出,如协同克里格[10]、回归克里格[11]、泛克里格[12]、指示克里格[13]、析取克里格[14]以及各种随机模拟方法[15-16]等,极大地丰富了地统计学的方法体系和内涵。

目前,有多个软件包含地统计功能,如ArcGIS的Geostatistics模块、GS+、GSLIB、Surfer等,各个软件均能执行部分地统计学方法。这些软件降低了使用门槛的同时,也促进了地统计学方法的普及,使得用户在不了解地统计理论和技术的情况下,也能利用地统计学方法对所分析的地理属性进行空间插值或不确定性评估。但从已发表的文献来看,由于不掌握地统计学的基本理论知识,造成一些失败的建模和计算结果,具体表现在以下几个方面:

(1) 变异模型精度不高,造成最终克里格插值结果精度较低,这可能是在计算试验变异函数时,选择了不合适的滞后间距,或是忽略了变异函数不同方向的各向异性等原因造成的；

(2) 只依赖样点数据,而未能有效地利用其他辅助数据来提供插值精度；

(3) 未能根据数据特点,选择合适的地统计学方法。

鉴于地统计学方法的广泛应用及其在使用过程中容易出现的失误,为普及地统计学基本理论,介绍地统计学方法体系,我校为地理信息科学和生态学专业的本科生和研究生开设了“地统计学方法”课程。由于该课程内容多、算法复杂、理论性较强,单独依靠理论课教学手段难以使学生理解抽象的概念和复杂的算法流程,因此,开设了相关实验教学环节来加深学生对理论和算法的理解并提高学生学习兴趣。

1 地统计学方法实验设计

地统计学的特点之一就是方法多样,但目前尚没有一个软件囊括了所有的方法,且有的方法实施过程(如回归克里格)需要借助多个软件完成,因此课程的实验过程也是多个软件(还包括非GIS或地统计软件)综合利用的过程。另外,本课程的主要授课对象是地理信息科学本科生,更注重学生对基本理论和方法的理解,若实验过程只依赖已发布的地统计软件或模块,这些理论和算法对学生来讲就是暗箱操作,因此,在部分实验环节,设计了程序设计环节,由学生利用基本的程序设计语言完成关键算法。

为了便于方法比较,本实验课程使用一套数据作为实验对象,数据为湖北省某县土壤样点数据,数据项包括:[x,y, 有机质(z1),土壤全氮(z2),pH值(z3)],其中前2项为样点的坐标数据,后3项为待分析的地理属性项。本数据共计658个样点,另外还包括该县的行政边界地形数据等。根据地统计学授课内容、实验数据和上述分析,设计了6个实验内容。

1.1 数据的基本统计分析和预处理

通过数据的基本统计分析能够了解待插值地理属性的基本统计特征,了解其概率分布特点,为正式计算前是否要对数据做某些预处理提供参考依据。基本统计分析内容包括最大值、最小值、均值、中值、标准差、变异系数、峰度系数、偏度系数、K-S检验等。其中,联合峰度系数、偏度系数、K-S检验的结果判断数据是否接近正态分布,从而是否要对属性数据进行转换,以使转换后的数据接近正态分布。数据预处理主要包括2个部分:一是待分析属性数据的异常值检验及处理,鉴于实验数据量大,推荐使用2倍标准差或3倍标准差法；二是数据转换,即若上一步判断原始数据不接近正态分布,则要对数据进行转换,转换方法包括对数、平方、反正弦转换等。本实验利用Excel和SPSS软件联合完成,要求学生以图表的形式提供结果,并为下一个实验提供接近正态分布的属性数据。

1.2 基于样点的试验变异函数计算

定义z(s)为在位置s=(x,y)处的地理属性值,即地理属性Z在位置s处的观测值,则试验变异函数的计算公式为:

(1)

其中h为给定的空间距离,N(h)为样点中,符合空间距离h的点对数量,z(si)和z(si+h)为空间距离为h的点对。取若干个等间隔的h,则能计算得到若干个γ(h)*,形成[h,γ(h)*]数据对,这就是基于样点的试验变异函数计算结果。当样点为规则取样时,只要以取样间距为距离间距,则能很容易地获取各距离h下的若干点对。但大部分的样点并非规则布设,实际计算时若严格按给定的距离寻找样点对,则出现样点对极少甚至可能找不到的可能。因此指导学生按照以下步骤进行计算:

(1) 假设步长为lag,当前滞后级别为n(n为正整数),则h=n*lag；

(2) 对所有样点,找到点对(si,sj),其符合条件:(n-0.5)*lag

(3) 计算[z(si)-z(sj)]2,记为Si；

(6) 将各个滞后距级别上的(havg,γ(havg)*)绘制成图,形成试验变异(半方差)图。

此实验要求学生根据上述算法,基于程序开发语言(如C#或Matlab)完成。

1.3 理论变异函数拟合

基于试验变异函数散点图,由学生根据图形特点选择一个合适的理论变异函数模型(如球状、高斯、指数、线性模型等)后,利用最小二乘法估计模型中的参数。但由于大部分理论模型均不是线性模型,因此,在参数估计前,需进行模型变换,即将非线性模型变换成线性模型,待拟合好线性模型中的参数后,再反算得出理论模型中的参数解[1]。另外,若选择套合模型作为理论模型,由于待估计的参数数量较多,可引导学生使用智能算法(如遗传算法)进行模型参数的求解[17]。模型求解后,要求学生利用试验变异函数散点的实际值及模型拟合值之间的相关系数和均方根误差进行精度验证,若精度较低,则重新进行参数估计,直到满足给定的精度要求为止。此实验要求学生根据算法,利用程序开发语言(如C#或Matlab)完成。

1.4 普通克里格算法实现及地图制图

普通克里格方法是目前运用最广泛的克里格方法,也是目前各领域最常用的地理属性空间插值算法,其核心思想是利用待估位置周围的若干个已知样点对地理属性值进行线性无偏最优估计,并给出估计方差。利用这个思想,推导出克里格方程组:

(2)

(1) 设定待插值栅格矩阵,主要是确定栅格单元大小,为制图美观,栅格单元大小不应大于50 m。设定好栅格矩阵后,每个栅格单元的中心点坐标作为该栅格的待插值位置,其插值结果代表了整个栅格单元的克里格估值结果；

(2) 按一定顺序遍历每个栅格,寻找离栅格单元的中心点最近的n个样点(一般为8个),构建克里格方程组,解该方程组,得到各样点的权重系数,进而得到待估值位置的克里格插值结果；

(3) 对栅格矩阵中每个栅格估值完后,将结果存储成ASC格式文件,导入到ArcGIS软件中,成为栅格图层。利用区域边界对栅格图层进行裁剪,只保留区域边界之内的部分,再添加指北针、比例尺、图例等地图整饰要素,最终形成区域土壤属性空间分布专题地图。

1.5 指示克里格实验

与其他克里格方法所得结果为地理属性值不同,指示克里格方法所得结果为地理属性超过或不超过某阈值概率的空间分布,常用于环境风险评价等领域。指示克里格方法步骤与普通克里格法大体相同,区别在于计算前,要将样点数据按照某个阈值做指示变换,具体实验步骤如下:

(1) 设定阈值c,将原始样点的属性数据做如下变换:

(3)

即将原始地理属性含量值按照阈值c变换成指示值；

(2) 基于上一步获得的指示值,计算试验变异函数,并拟合理论变异函数模型；

(3) 构建区域栅格矩阵,基于指示值和理论变异函数模型,对各栅格中心位置进行普通克里格估值,得到每个栅格超过阈值c的概率；

(4) 利用区域边界数据进行边界裁剪,并添加指北针、比例尺、图例等地图整饰。

该实验可在ArcGIS中的Geostatistics模块上完成。

1.6 回归克里格实验

(1) 在ArcGIS中,提取样点上各环境变量的值；

(2) 在SPSS软件中,计算土壤有机质与各环境变量的相关系数；

(3) 在SPSS中,用逐步回归法建立土壤有机质与各环境变量的线性回归方程,即得到趋势模型m(s)；

(4) 使用ArcGIS的地图代数功能,利用上一步得到的趋势模型,计算各位置上的趋势值；

(5) 提取趋势值到样点图层中,可在Excel中得到各样点的残差值；

(6) 利用ArcGIS中的Geostatistics模块完成样点残差的克里格插值；

(7) 使用ArcGIS的地图代数功能,将各位置上的趋势值与残差的克里格插值结果相加,得到回归克里格预测结果；

(8) 进行边界裁剪,添加地图整饰,完成地图制图并提交。

2 总结与讨论

综上所述,地统计学实验课程分为6个内容,其中第2、3、4个内容需要学生利用程序设计语言自行开发,通过这些实验过程,能加深学生对地统计学最基本理论和最核心算法的理解和掌握,但也对学生提出了更高的要求,即需要有良好的程序设计基础。第1、5、6个实验需要学生综合利用多个软件完成,其中第5、6个实验可以看作是对基本的统计学方法的扩展,拓宽学生视野的同时,锻炼了学生多种工具软件的综合使用能力。而在第6个实验的教学过程中,只介绍回归克里格的思想,但不明确给出具体的实验步骤和所使用的软件,由学生自行探索,实验完成后,提升了学生解决问题的能力。另外,第4、5、6个实验中地图制图的内容能够巩固地理信息科学专业的学生地图学的知识和技能。

此实验方案已实施3年,并不断地修改和完善,根据实验后学生的反馈,通过实验获得的收获包括:

(1) 通过对变异函数和普通克里格插值的算法实现,深刻地理解了地统计学的基本理论和算法,同时在空间数据结构、程序设计语言、数据和文件间的I/O等方面均有深刻体会和提升；

(2) 同归指示克里格和回归克里格实验,印证了课堂的理论教学内容,拓展了视野；

(3) 通过地图制图,巩固了地图学中关于地图整饰和地图符号设计知识；

(4) 通过观察所得土壤属性空间分布结果,理解了地理属性空间分异规律和特征。

但该实验方案还存在一些问题,主要是第2、3、4内容难度太大,对一些程序设计基础不好的学生来讲,难以开展。这需要在实验前对学生的技术能力进行摸底调查,如果大部分学生技术力量薄弱,那么可以通过实验分组来解决。即将多个学生分成一个组,每个学生只负责其中的一部分,共同协调合作完成这部分实验内容。