R语言在地统计学中的应用-原作强

原作强中国科学院沈阳应用生态研究所[e-mail]语言在地统计学中的应用UsingRtoimplementgeostaticalanalysisinopensourceenvironment地统计方法的基本原理前提假设区域化变量变异分析空间估值应用实例以geoR包为例概念和应用范围地统计(Geostatistics)又称地质统计，是在法国著名统计学家G.Matheron大量理论研究的基础上逐渐形成的一门新的统计学分支。它是以区域化变量为基础,借助变异函数,研究既具有随机性又具有结构性，或空间相关性和依赖性的自然现象的一门科学凡是与空间数据的结构性和随机性,或空间相关性和依赖性，或空间格局与变异有关的研究,并对这些数据进行最优无偏内插估计,或模拟这些数据的离散性、波动性时,皆可应用地统计学的理论与方法生态学研究中的生物或非生物因子的时空依赖性是普遍存在的，如种群梯度分析、种群格局、环境因子的时空分布的异质性等等地统计学与经典统计学相同点：它们都是在大量采样的基础上，通过对样本属性值的频率分布或均值、方差关系及其相应规则的分析, 确定其空间分布格局与相关关系不同点：地统计学既考虑到样本值的大小，又重视样本空间位置及样本间的距离，弥补了经典统计学忽略空间方位的缺陷，能够对研究对象的空间格局进行检验、模拟和估计地统计分析理论基础—前提假设(1)随机过程所有样本值都不是相互独立(2)正态分布若不符合正态分布的假设，应对数据进行变换，并尽量选取可逆的变换形式(3)平稳性(a)均值平稳(b)与协方差函数有关的二阶平稳和与半变异函数有关的内蕴平稳当一个变量呈现为空间分布时，就称之为区域化变量(regionalizedvariable).这种变量常常反映某种空间现象的特征，用区域化变量来描述的现象称之为区域化现象.区域化变量，亦称区域化随机变量，G.Mathero(1963)将它定义为以空间点x的三个直角坐标为自变量的随机场.区域化变量具有两个最显著，而且也是最重要的特征，即随机性和结构性.地统计分析理论基础—区域化变量协方差又称半方差，表示两随机变量之间的差异。在概率论中，随机变量X与Y的协方差定义为：协方差函数式中：h为两样本点空间分隔距离；Z(xi)为xi处的实测值；Z(xi+h)为Z(xi)在xi处距离偏离h的实测值;N(h)为分割距离是h时样本点的对数地统计分析理论基础—变异分析半变异函数又称半变差函数、半变异矩，是地统计分析的特有函数。区域化变量Z(x)在点x和x+h处的值Z(x)与Z(x+h)差的方差的一半称为区域化变量Z(x)的半变异函数，记为r(h)，2r(h)称为变异函数。协方差函数地统计分析理论基础—变异分析半变异函数和协方差函数把统计相关系数的大小作为一个距离的函数，是地理学相近相似定理定量量化地统计分析理论基础—变异分析Title•. 块金值(Nugget):由于测量误差和空间变异,得两采样点非常接近时，它们的半变异函数值不为0基台值（Sill）：即函数值不随采样点间隔距离而改变时，空间相关性不存在偏基台值(partialSill)：基台值与块金值的差值变程(range):当半变异函数的取值由初始的块金值达到基台值时，采样点的间隔距离空间相关性的强弱可由PartialSill/Sill来反映,越大,间相关性越强Nugget/Sill称为基底效应，表示样本间的变异特征，该值越大，更多得是由随机因素引起的变异函数的理论模型地统计学将变异函数理论模型分为3大类：第1类是有基台值模型，包括球状模型、指数模型、高斯模型、线性有基台值模型和纯块金效应模型；第2类是无基台值模型，包括幂函数模型、线性无基台值模型、抛物线模型；第3类是孔穴效应模型。下面有代表性地介绍几种常见的变异函数理论模型。克立格（Kriging）插值法，又称空间局部估计或空间局部插值法，是地统计学的主要内容之一。克立格法是建立在变异函数理论及结构分析基础之上的，它是在有限区域内对区域化变量的取值进行无偏最优估计的一种方法。克立格法适用的条件是，如果变异函数和相关分析的结果表明区域化变量存在空间相关性。其实质是利用区域化变量的原始数据和变异函数的结构特点，对未采样点的区域化变量的取值进行线性无偏、最优估计。地统计分析理论基础—空间估值克立格插值Ｋ(riginginterpolation)是根据变异函数模型而发展起来的一系列地统计的空间插值方法，包括：普通克立格法(OrdinaryKriging);泛克立格法(UniversalKriging);指示克立格法(IndicatorKriging);析取克立格法(DisjunctiveKriging);协同克立格法(Cokriging)等Kriging&stochasticsimulationKrigingandstochasticsimulationaregeostatisticalmethodswhicharecommonlyusedtocreatecontinuousmaps.Thetwomethodshaveopposinggoals:1.krigingaimsatlocalaccuracythroughminimizingacovariance-basederrorvariance,krigingproducesoneoutputmap.krigingispreferredforlocalestimation2.simulationaimsatreproducingspatialstructurethroughacovariancemodel.stochasticsimulationgeneratesmultiplerealizationsofthespatialdistributionof(soil)attributevaluesanditusesdifferencesamongsimulatedmapsasameasureofuncertainty.simulationisincreasinglypreferredforassessmentofspatialuncertaintyandreproductionofglobalstatistics,riskassessment,flowmodeling,andwaterqualitysimulationmodeling实例•Packages:geoR& sp •Author：RibeiroIJrDiggle,P.J•URL：•Book:Model-basedGeostatistics•URL:数据建立数据导入：s100=as.geodata(soil2009,coords.col=2:3,data.col=4)s100=read.geodata(“soil.txt”,coords.col=1:2,data.col=3,header=T)2.探索性分析summary(soil)plot(s100)3.数据检验1.坐标点的检查dup.coords(s100$coords)2.Shapiro-Wilk检验hist(s100$data)shapiro.test(s100$data)p-value=0.02474长度：3-5000p0.0001拒绝零假设,即不符合正态分布3.数据的转换(Boxcox)boxcox.fit(abc$data)shapiro.test(abc$data)p-value=8.601e-14Fittedparameters:lambdabetasigmasq0.1310.6861.682shapiro.test(((abc$data)^0.131-1)/0.131)P=0.00021834.计算经验方差图1.各项异性计算四个方向(00、450、900and1350)vario.4-variog4(s100,max.dist=1.0)plot(varog.4)2.各项同性计算bin1=variog(s100,trend=cte,lambda=1,uvec=seq(0,1,l=13),max.dis=1)plot(bin1)5.拟和理论方差图画出散点：plot(bin1)参数估计：1.byeye2.byleastsquaresfitofempiricalvariogramsOrdinary(OLS)weighted(WLS)leastsquares3.bylikelihoodbasedmethodsMaximumlikelihood(ML)andrestrictedmaximumlikelihood(REML)4.Bayesianmethods5.拟和理论方差图1.Fittingmodelswithnuggetfixedtozerowls-variofit(bin1,ini=c(1,0.5),fix.nugget=T)ols-variofit(bin1,ini=c(1,0.5),fix.nugget=T,weights=equal)ml-likfit(s100,ini=c(1,0.5),fix.nugget=T)reml-likfit(s100,ini=c(1,0.5),fix.nugget=T,method=RML)2.Fittingmodelswithafixedvalueforthenuggetwls.fn-variofit(bin1,ini=c(1,0.5),fix.nugget=T,nug=0.15)ols.fn-variofit(bin1,ini=c(1,0.5),fix.nugget=T,nug=0.15,weights=equal)ml.fn-likfit(s100,ini=c(1,0.5),fix.nugget=T,nug=0.15)reml.fn-likfit(s100,ini=c(1,0.5),fix.nugget=T,nug=0.15,method=RML)3.Fittingmodelsestimatednuggetwls.n-variofit(bin1,ini=c(1,0.5),nug=0.5)ols.n-variofit(bin1,ini=c(1,0.5),nug=0.5,weights=equal)ml.n-likfit(s100,ini=c(1,0.5),nug=0.5)reml.n-likfit(s100,ini=c(1,0.5),nug=0.5,method=RML)5.拟和理论方差图wls-variofit(bin1,ini=c(0.878,0.333),cov.model=exponential)eyefit(bin1)matern,exponential,gaussian,spherical,circular,cubic,wave,power,powered.exponential,cauchy,gencauchy,gneiting,gneiting.matern,pure.nugget.plot(bin1,main=“长白山)lines(wls)5.包基线零假设：数据没有空间相关关系通过置换数据空间位置,进行99次模拟,取最大和最小值env.mc-variog.mc.env(s100,obj.var=bin1)plot(bin1,env=env.mc,main=)lines(ols.n,col=red,lwd=2)6.数据验证1.交互验证xv.ols=xvalid(s100,model=ols.n)2.部分验证s100.xv.ml=xvalid(s100,model=wls,locations.xvalid=loci,data.xvalid=c(1.0,1.14,0.12,-0.55))plot(s100.xv.ml)names(xv.ols)[1]datapred