地统计学的基本概念及公式详解
本文对插值、平稳假设、变异函数、克里格等常用的地学计算概念加以介绍，并对相关公式进行推导。
1 引言
最近的几篇博客，分别从遥感的实际应用出发，对影像前期处理与相关算法、反演操作等加以详细介绍。而通过遥感手段获取了丰富的各类地表信息数据后，如何对数据加以良好的数学处理与科学分析，同样是我们需要重视的问题。因此，准备由这一篇博客入手，逐篇地对地学计算方面的内容加以初步总结。
那么首先，我们就由地学计算的几个基本概念入手，对相关理论方面的内容加以一定了解。
需要注意的是，以下内容如果单独来看或许有些不好理解，但一旦将其与实际应用结合，便会豁然开朗。其中，具体的实际应用部分我将会在后面的博客中涉及。
2 空间插值
空间数据的获取是进行空间分析的基础与起源。为了提高研究结论的精度，我们亦总是希望能够获取研究区域内更多、更全面的精确空间属性数据信息。然而，在实际研究、工作中，由于人力、成本、资源等外部条件限制，我们不可能对全部未知区域加以采样与测量，而往往只能得到研究区域内有限数量的采样点及其相关属性数据。因此，往往可以考虑选取合适的空间采样点，利用一定数学模型，依据已知采样点各自对应属性数据对研究区域所有位置的未知属性信息加以预测。
空间插值（Spatial Interpolation）即可以实现这一需求。其是一种将离散采样点测量数据转换为连续数据曲面的常用方法，包括内插（Interpolation）和外推（Extrapolation）两种应用形式。一般地，对样本点范围以内（即所有采样点最大外接矩形内部）的空间进行插值，才可称作“内插”（部分文献亦直接用“插值”代替“内插”）；反之则称“外推”或“预测”，往往认为外推的结果误差较大。
空间插值理论及其方法基于著名的“地理学第一定律（Tobler's First Law of Geography）”，即一般地，距离越近的地物具有越高的相关性。这一至今已产生深远影响的地学定律最早由美籍瑞士地理学家Waldo R. Tobler教授于1970年提出。
在各方法所对应的数学计算原理层面，空间插值一般可以分为确定性插值方法（Deterministic Interpolation）与地统计插值方法（Geostatistics，亦称非确定性插值方法）两种。其中，确定性插值方法基于研究区域内各信息点之间相似程度或整个曲面的平滑程度，从而创建连续的拟合曲面；其依据插值计算时纳入考虑的采样点分布范围，又可进一步分为整体插值法与局部插值法。地统计插值方法则是基于研究区域内各信息点的综合统计学规律，以变异函数（Variogram）理论与结构分析为基础，实现其属性的空间自相关性定量化，从而创建得出连续插值曲面。
在所创建连续插值表面通过全部采样点的与否层面，空间插值一般又可以分为精确性插值与非精确性插值两种。其中，前者预测样点的属性数值与其各自实测值相等，即其采样点属性数据全部落入预测结果曲面；后者预测样点的属性数值则往往不与其各自实测值相等，即其采样点属性数据一般不会落入预测结果曲面。因此，使用非精确性插值方法往往可以避免在预测表面中出现明显的波峰或波谷，整体呈现出平缓态势。
3 几个重要假设
地学计算中，几个重要假设可以说是关键中的关键；而其往往也比较难以理解，不怕，我们慢慢往下看。
3.1 平稳假设
平稳假设（Stationary Assumption）是指，一组观测值的均值是始终固定的，其与观测值所在位置无关；将既定的某个点集由某一研究区域内某处移动至另一处时，随机函数的性质保持不变。即随机函数的分布规律不因位置的改变而改变，具有严格的平稳性。平稳性假设的公式表达为：
其中，F_(x_1,⋯,x_n ) (z_1,⋯,z_n )表示位置在(x_1,⋯,x_n )上的点集(z_1,⋯,z_n )对应的随机函数。
3.2 二阶平稳性假设
二阶平稳性假设（Second Stationary Assumption）又称弱平稳假设，其与协方差函数（Covariance Function）有关。这一假设认为，随机函数的均值为一常数，且任意两个随机变量之间的协方差仅仅依赖于其二者之间的距离与方向，而与其具体位置无关。
将上述两个条件用公式分别表达为：
其中，E[Z(x)]为区域化变量Z(x)的数学期望， Cov[Z(x),Z(x+h)]为区域化变量Z(x)与Z(x+h)所对应的协方差函数，C(h)为仅与h有关的协方差取值，m为一常数，h为滞后距。
上述二阶平稳性假设是针对整个研究区域范围而言。若区域化变量仅仅在整个研究区域内的某个有限区域中满足上述条件，即条件仅在局部区域成立，则称此区域化变量满足准二阶平稳性假设（Quasi Second Stationary Assumption）。准二阶平稳性假设可以视作一种折中方案，既考虑到平稳范围的大小，又顾及到有效数据的数量。
3.3 本征假设
本征假设（Intrinsic Hypothesis）又称内蕴假设，其与变异函数有关。这一假设认为，区域化变量的增量满足以下两个条件：在整个研究区域内，区域化变量增量的数学期望为0；且其方差函数存在，并只依赖于滞后距，而与所处位置无关。
将上述两个条件用公式分别表达为：
当E(x)存在时，上述第一个公式可以写作：
其中，Var[Z(x)-Z(x+h)]为区域化变量Z(x)与Z(x+h)所对应的方差函数，γ(h)为区域化变量在滞后距为h时的变异函数，m为一常数，h为滞后距，其它符号同前述意义。
同样的，上述本征假设亦是针对整个研究区域范围而言。若区域化变量仅仅在整个研究区域内的某个有限区域中满足上述条件，则称此区域化变量满足准本征假设（Quasi Intrinsic Hypothesis）。与准二阶平稳性假设类似，准本征假设亦可视作一种折中方案，其同样既考虑到了本征假设对应范围的大小，又顾及到了有效数据的数量。
此外，本征假设是地统计学中对随机函数的基本假设。
3.4 不同假设对比
结合上述二阶平稳性假设与本征假设相关原理，可以看到两种假设的讨论对象具有一定区别：二阶平稳性假设更多是讨论某一特定研究区域内区域化变量自身【即Z(x)】的特征，而本征假设则是研究区域化变量所对应增量【即Z(x)-Z(x+h)】的特征。
一般认为，二阶平稳性假设对区域化变量的要求较之本征假设更为严格，即若某个研究区域内的某一区域化变量满足二阶平稳性假设，则其一定满足本征假设；反之则反，若仅知区域化变量满足本征假设，则其不一定满足二阶平稳性假设。
再结合平稳假设，上述三种假设的严格程度由大至小依次排列为：平稳假设、二阶平稳性假设与本征假设。
此外，结合二阶平稳性假设的两个条件，还可以推出协方差函数与变异函数之间的关系：
其中，γ(h)为区域化变量所对应变异函数，C(0)为距离为0时此区域化变量所对应的协方差取值【即基台值γ(∞)】，C(h)为距离为h时此区域化变量所对应的协方差取值。
由这一关系可知，用以衡量某一区域化变量在相距为h的两空间位置点分别取值的自相关性的指标——协方差函数与变异函数之间具有相互关系。因此，在满足二阶平稳性假设的条件下，若协方差函数平稳，则可知变异函数平稳，即其取值只与滞后距h有关。
4 变异函数
克里格插值法需要借助空间数据的试验变异函数及其散点图特点，因此变异函数的计算在克里格插值过程中发挥着重要作用；变异函数及其模型拟合对克里格插值结果精度具有较大影响。
变异函数（Variogram），又称为半变异函数、半方差函数（Semi-variogram）等，其用以描述区域化变量的空间变化特征与强度，被定义为区域化变量增量平方的数学期望。
在一维条件下，直接将区域化变量Z在位置(x)与(x+h)处的取值Z(x)与Z(x+h)之差的方差定义为变异函数，其因变量为距离h；而在二维或三维条件下，可以将上述一维中具有单一方向的距离h进一步引申为在任意方向α上的距离|h|。具体公式表达为如下。
其中，γ(x,h)即为变异函数。由于公式中在其前具有一个系数2，因此其亦被称作半变异函数。结合（准）二阶平稳性假设、（准）本征假设等地学基本假设，变异函数取值与区域化变量样点所处位置x无关，仅仅与样点之间的距离h有关，则变异函数可以写作：
其中，〖γ(h)〗^#为区域化变量Z(x)的变异函数，N(h)为区域化变量样点集中距离为h的点对数量，x_i为距离为h时所对应的第i个点。在这里，距离h亦被称作滞后距（Lag Distance）。
一般地，区域化变量变异函数图像往往呈现出“先快速上升，再增速减缓，后趋于平稳”的曲线特征。其具有三个十分重要的相关概念，分别为块金常数（Nugget）、基台值（Sill）与变程（Range）。
块金常数代表区域化变量的随机性大小。由理论角度，在间距为0（即滞后距为零）时，区域化变量采样点数值应当相等；而在间距无限趋近于0时，对应变异函数数值应当亦向0趋近。但是，在实际研究中，试验变异函数在滞后距为0时，其取值并不为0，而是一个大于0的数值。这一数值便称为块金常数。一般地，上述块金效应的产生可以归因于测量误差，或小于采样间隔距离处的空间变化。
基台值用以衡量区域化变量变化幅度的大小。当滞后距无限增大并到达某一程度后，试验变异函数若趋于平稳，则这一平稳水平所对应的数值即为基台值。然而，并不是所有的区域化变量均具有基台值——如无基台值模型对应的变异函数。
变程用以衡量区域化变量自相关范围的大小。当滞后距无限增大并到达某一程度后，试验变异函数若趋于平稳，则此时对应的滞后距即为变程。其中，小于变程的距离所对应的样本位置与空间自相关，而大于变程的距离所对应样本位置不存在空间自相关。
此外，变异函数还有其它相关指标，如基台值与块金常数的差值——偏基台值（Partial Sill），用以衡量空间变异性程度的块金常数与基台值的比值——块金系数等。
基于不同区域化变量对应的变异函数特征，可以将其分为不同类别。依据变异函数基台值的有无，可以将模型分为有基台值模型、无基台值模型与孔穴效应模型。
其中，有基台值模型可以依据变异函数的特征进一步分为纯块金效应模型（Pure Nugget Effect Model）、球状模型（Spherical Model）、指数模型（Exponential Model）、高斯模型（Cubic Model或Gaussian Model）与线性有基台模型（Linear with Sill Model）等。上述几种模型中，较为常用的模型包括球状模型、指数模型与高斯模型。
此外，同样依据变异函数的特征，无基台值模型还可进一步分为线性无基台值模型、幂指数模型与对数模型等。同样的，孔穴效应模型可分为基台值模型和无基台值模型。
同时，针对某种区域化变量而言，其在不同方向、不同滞后距情况下可能受到不同因素影响；套合结构可以很好解决这一问题。套合结构可以表示为多个变异函数之和，每一个变异函数均代表着某种方向或某一尺度中的变异性，从而对区域化变量的特征加以更好概括。
5 克里格插值
克里格插值法（Kriging Method）又称为空间局部插值法，是以上述变异函数理论及其结构分析为基础，在有限区域内对区域化变量进行线性无偏最优估计（Best Linear Unbiased Prediction，BLUP）的一种方法，在地统计学中也被称为空间最优无偏估计器（Spatial BLUP）。其中，上述“线性”是指克里格插值法对未知点属性数值的估计采用线性估计，其公式如下：
其中，(z_0 ) ̂是区域化变量在点(x_0,y_0 )位置处的预测值，λ_i为第i个已知点的权重系数，z_i为第i个已知点的实测值。
上述权重λ_i可以使得各点处的预测值与实测值之间的方差最小，而求取这一权重则为克里格插值的主要内容。
上述“无偏”是指区域化变量在各点上估计量的数学期望等于其在同一位置上的真实值，公式如下：
结合本征假设，无偏性可以表示为：
上述“最优”是指区域化变量在各点上估计量与其在同一位置上的真实值的方差最小，公式如下：
其中，上述方差被称作估计方差或估值方差，是对估值准确程度的一种定量表示；而在克里格插值方法中，又可以称为克里金方差。后续将克里金方差记作σ_k^2。
经过统计学相关推导，可以将克里金方差写作：
由此转换为在无偏条件约束下的最小值求解问题。引入拉格朗日乘子φ，构造拉格朗日函数：
对权值与拉格朗日乘子求一阶偏微分：
偏微分求解结果为：
将上述求和部分展开：
可将上式进一步写为矩阵相乘形式，即化为：
其中，A代表在原有变异函数矩阵中额外添加全1行与全1列（交界处1换为0）后的矩阵，λ代表各权重组成的列向量，φ代表前述分析引入的拉格朗日乘子，B为各位置与待求解位置对应距离的变异函数值组成的列向量，且在列尾增加一个1。
由此，即将上述函数转化为(n+1)个未知数、(n+1)个表达式组成的方程组；通过矩阵求逆，求解方程组即可得到待求解位置与其它已知点的权重。对每一个待插点进行同样操作，完成克里格插值。
6 回归克里格
正如上述分析，普通克里格（Ordinary Kriging）方法更多依靠所选采样点数据对空间加以插值，其插值效果较依赖于采样点的个数、密度，以及其数据的准确性；另一方面，许多空间属性往往受到其它环境变量的影响。例如，土壤有机碳含量与气温、海拔、降水量、坡度等多种环境因子在0.05水平显著相关；近地面气温与海拔、海陆距离、NDVI等环境因子在0.01水平显著相关。由此观之，若简单地忽略环境要素对待插值空间属性的影响，可能会降低最终插值结果精度。
基于这种考虑，可以使用回归克里格（Regression Kriging）方法对环境因素加以考虑。应用回归克里格插值，首先需要确定目标变量与辅助环境变量之间的相关性；确认相关关系符合要求后，建立目标变量与环境变量的一元或多元回归模型，从而剔除空间趋势项。随后，依据采样点实测数据与回归模型计算得出的对应位置数值，求得目标变量的确定性趋势项。运用普通克里格方法，将残差进行插值，并最终将回归预测的趋势项与普通克里格的插值结果相加，从而得到目标变量估测值。
最终空间插值结果的公式表达为：
其中，Z(x_0 )为回归克里格估值在x_0位置的结果，m ̂(x_0 )为确定性趋势项在这一位置的取值，e ̂(x_0 )为残差普通克里格插值在这一位置的取值。对研究区域内每一位置执行同样操作，完成回归克里格插值。
一般地，对于所受外界影响较大的空间属性，回归克里格插值效果优于普通克里格插值。
另一方面，如前所述，尽管回归克里格方法对于残差项的估计仍采用普通克里格插值计算，但由于其对最终目标变量插值结果的影响较之由环境因子得到的空间确定性趋势项而言仍较低，因此回归克里格方法所得结果较之普通克里格方法结果将会粗糙、破碎一些。此外，由方法考虑的因素这一角度来看，回归克里格方法与协同克里格（Cokriging）方法具有一定相同之处，即二者均利用相关环境因子对目标变量的空间分布加以估测；而其不同之处在于，回归克里格的辅助数据为面状分布，如栅格图层；而协同克里格的辅助数据为点状分布。但无论采取何种方法，最终所得插值结果均为一个完整的表面。