Quantifying spatial heterogeneity at the landscape scale using variogram models 文献翻译与解读

本人的解读用底色标示。

1-3 引言、文献综述、数据说明

略

4 方法

  量化景观的空间异质性需要几个假设：

• H1：图像的范围相比于要研究的空间特征足够大 （H1的验证标准见4.3节）
• H2：电磁测量误差相对于地表变化小到可以忽略
• H3：尺度小于采样步长的空间变异性可以忽略
• H4：高空间分辨率的影像数据是点观测值（punctual，对应H3，认为像元是均质的，相关名词解释见第8节补充材料 ）

4.1 实验变异函数

  NDVI数据可以认为是区域化变量 $z(x)$ （例如地理位置的函数）的点观测值，它描述了景观植被在影像 $I$ 上的空间分布。实验变异函数（experimental viogram）通过观测数据计算得到。它度量了不同距离下像元值之差平方的平均值：
$${\gamma }_e(\mathbf{h})=\frac{1}{2N(\mathbf{h})}\sum _{||x_{\alpha }-x_{\beta }||\approx \mathbf{h}}(z(x_{\alpha })-z(x_{\beta }){)}^2$$其中$x_{\alpha }、x_{\beta }$是像元值，$\mathbf{h}$是像元之间的距离向量（本文后续内容基于各向同性，不考虑方向，因此$\mathbf{h}$退化为一个表示距离的标量。考虑各项异性时也可据此推导）。
  在距离很大时变异函数在统计意义上是不可靠的(Chilès & Delfiner, 1999)，因此本文设置了计算变异函数的距离阈值 $d_{max}=1500m$，也就是影像边长的一半。实验变异函数通常有几个重要特征(Chilès & Delfiner, 1999)：

• 通常是距离 $||h||$的递增函数(Fig. 1)，因为距离近的像元通常比距离远的像元更相似。
• 距离大到一定程度后它通常会达到基台值（sill）或无限增大。VALERI数据库中大多数实验变异函数在距离达到 $d_{max}$ 之前就已经达到基台值。
• 基台值（sill）也是空间变异性的表现。变程（range）是变异函数达到基台值时的距离。如果距离超过了$d_{max}$ 还没有达到基台值，则认为影像的范围不能完全包含该空间结构 （违反了H1）。
• 变异函数在原点附近的特征表现了研究变量的连续性。距离为0时的变异函数值被称为块金值（nugget），理论上应该为0。如果不为0可能是由于：测量误差或空间结构的长度规模比像元尺寸还要小 （像元内部的异质性）。 VALERI数据中所有的实验变异函数在原点附近都是线性的，没有任何块金效应。这一发现有力地支持了假设H2和H3。

4.2 变异函数建模

  实验变异函数为NDVI空间分布提供了一种二阶统计的经验性描述。使用一种基于概率的模型对区域变量 $z(x)$ 的空间异质性成分进行参数化刻画。$z(x)$ 是所有可能的随机函数 $Z(x)$ 的一种实现。 $Z(x)$ 的二阶平稳性假设它的一阶和二阶矩存在且是平稳的。对于所有$x$和$\mathbf{h}$有：
$$E[Z(x)]=mCov(Z(x),Z(x+\mathbf{h}))=C(\mathbf{h})$$函数 $C(\mathbf{h})$ 是 $Z(x)$ 的协方差函数，它描述了 $Z(x)$ 的空间分布。在二阶平稳性假设下，理论变异函数
$$\gamma (\mathbf{h})=0.5Var[Z(x+\mathbf{h})-Z(x)]$$与协方差函数存在关系
$$\gamma (\mathbf{h})={\sigma }^2-C(\mathbf{h})$$其中 ${\sigma }^2$ 是 $Z(x)$ 的方差。理论变异函数在 $||\mathbf{h}||=0$ 时值为0，最终当 $||\mathbf{h}||$ 趋近于无穷时收敛到基台值（sill）。理论变异函数的变程（range）是它达到基台值时的距离大小。距离超过变程的数据是不相关的。
  理论变异函数通过用一些函数拟合经验变异函数得到。这些函数（也被称为authorized model）必须是conditionally negative function。本研究使用指数函数和球函数，因为它们满足了经验变异函数的主要特征：原点附近接近线性且连续，并且能收敛到一个基台值。Fig 2a和Table 3描绘了这两种基本模型和它们的主要性质。需要注意的是指数模型的变程参数也被称为实际变程，当距离达到实际变程时变异函数的值为基台值的95% （令h=r，exp(-3)约为5%）。


  为了解释数据的多重长度尺度，我们使用以上两种函数的线性组合来对变异函数进行建模。这一变异函数的扩展模型被成为区域化线性模型（linear model of regionalization），是 $l$ 个基本变异函数模型 $g_k(r_h,\mathbf{h}),k=1,...,l$ 的加权之和：
$$\gamma (\mathbf{h})={\sigma }^2\sum _{k=1}^{k=l}{b}_k{g}_k(r_k,\mathbf{h})$$   这一模型尤其适合用来描述一组独立的、有着不同长度尺度和空间变异性的、叠加在同一区域的空间结构。它描述了引言部分定义的空间异质性的两个组成部分（Fig 2b）：

• 影像的空间变异性程度通过基台值 ${\sigma }^2$ 给出
• 影像的空间结构通过每个基本模型 $g_k$ 的相应参数刻画：变程 $r_k$ 和每个变程占总方差的比例 $b_k$。

   $\gamma (\mathbf{h})$ 的参数通过半自动拟合法（Istatis software）估计。必要的基本变异函数的数量 $l$ 和其变程根据视觉调整。我们发现令 $l$ 等于1或2总是合适的。随后，基台值和方差权重使用加权平均平方优化法（Cressie，1985）得到。由于变程超过 $d_{max}$ 的实验变异函数是不可靠的，任何估计出的变程超过这一阈值的变异函数也被认为是不可靠的，因此也不被采用。这种情况下二阶平稳性假设不成立。

4.3 图像空间结构的刻画

  我们用一个参数变程积分 $A$（integral range）来概括变程 $r_k$ 和方差占比 $b_k$ 提供的变异函数的结构信息。一个二阶平稳随机函数 $Z(x)$ 的变程积分定义为（Chilès and Delfiner, 1999, Lantuéjoul, 2002）：
$$A=\frac{1}{{\sigma }^2}{\int }_{\mathbf{h}\in {\Re }^2}({\sigma }^2-\gamma (\mathbf{h}))d\mathbf{h}$$   Table 3给出了球和指数模型的变程积分值（Lantuéjoul, 2002）。对于区域化的线性模型，$A$ 是每个基本模型变程积分 $A_k$ 的加权线性组合：
$$A=\sum _{k=1}^l{b}_k{A}_k,A_k={\int }_{\mathbf{h}\in {\Re }^2}(1-g_k(r_k,\mathbf{h}))d\mathbf{h}$$   变程积分是一种面积矩（Serra, 1982）。遍历性协方差理论（ergodic covariance theorem）（Chilès and Delfiner, 1999, Lantuéjoul, 2002）指出如果图像区域 $I$ 相对于变程积分足够大，那么 $Z(x)$ 在 $I$ 上的空间平均值的方差（记为 $\overline{Z_I}$ ）为：
$$Var(\overline{Z_I})\approx \frac{{\sigma }^{2}A}{|I|}$$    因此，如果与传统的 $M$ 个独立点的算数平均数的方差公式相对照，比率 $M=|I|/A$ 可以被看作 $I$ 中独立值的数量。变程积分 $A$ 可以被视为它们的等价影响面积。因此它的平方根 $D_c$ 是 $M$ 个“独立数据”摆放在规则方形格网上时相互之间的距离。 $D_c$ 将拟合模型的所有的结构参数概括为一个距离特征，它是不同变程的加权平均。我们将其称为NDVI影像的“平均长度尺度”（mean length scale） 。它与影像空间结构的平均大小有关。$A$ 和 $D_c$ 可以被用于评判影像尺度是否足够大以至于能在一张影像中探测出景观的长度尺度（假设H1）。我们规定当 $\overline{Z_I}$ 的方差小到可以忽略时假设H1成立，例如变程积分 $A$ 小于影像大小 $I$ 的5%。因此，对3×3km的影像，如果 $A$ 小于阈值 $A_{T,3km}=4.5\times 10^5{m}^2$时，也就是 $D_c$ 小于 $D_{c,T,3km}=671m$，影像可以认为是足够大，其空间结构可以通过变异函数刻画。在本研究使用的18张影像中，12张满足假设H1。（$D_c$ 是本文核心指标）

4.4 空间异质性的量化

  前文将变异函数用于在影像尺度上描述NDVI的空间分布。变异函数同样可以量化影像子区域的空间异质性。下文中，通过移动通用区域块（block domain）将高空间分辨率的影像分割成 $P$ 个一致的块（block） $v_i$ (Atkinson, 2001; Chilès & Delfiner, 1999; Myers, 1997)（块对应了某种分辨率的像元）。图像的空间变异性被分解到三个尺度：

• 影像尺度： 影像中 $N$ 个高空间分辨率像元值的实验离散方差（dispersion variance）（i.e. 影像方差）
  $$s^2(x|I)=\frac{1}{N}\sum _{\alpha =1}^N(z(x_{\alpha })-\overline{z_I}{)}^2$$其中 $\overline{z_I}$ 是整幅影像上 $z(x)$ 的均值。
• 块间尺度： 影像中 $P$ 个块值的实验离散方差（i.e. 块之间的变异性）（以块为单位，形式和上面以像元为单位一样）
  $$s^2(v|I)=\frac{1}{P}\sum _{i=1}^P(\overline{z_{v_i}}-\overline{z_I}{)}^2$$其中 $\overline{z_{v_i}}$ 是区域块 $v_i$ 中 $z(x)$ 的均值。
• 块内尺度： 区域块内 $n$ 个高分辨率像元值实验离散方差的平均值（i.e. 区域块内部的变异性） （先计算每个块内的平均值，再对所有区域块的平均值求平均）
  $$s^2(x|v)=\frac{1}{P}\sum _{i=1}^P\frac{1}{n}\sum _{\alpha =1}^n(z(x_{\alpha })-\overline{z_{v_i}}{)}^2$$

  它们三者的关系为
$$s^2(x|I)=s^2(v|I)+s^2(x|v)$$

  从Fig. 3可以看出，当块尺寸 $v$ 增大时，块内变异性 $s^2(x|v)$ 增大、块间变异性 $s^2(v|I)$ 减小。

  $s^2(x|v)$ 经验性地量化了当几个像元被合并为一个更粗的分辨率的像元时图像变异性的损失。它可以直接根据高分辨率影像的理论变异函数计算得到。在随机函数二阶平稳性框架下，$s^2(x|v)$ 是随机变量 $S^2(x|v)$ 的一种实现。区域 $v$ 中 $Z(x)$ 的理论离散方差 ${\sigma }^2(x|v)$ 定义为$S^2(x|v)$ 的数学期望(Wack- ernagel, 2003)： （注：论文原文等号右侧是向下取整符号，但查看了引用的文献发现是普通中括号）
$${\sigma }^2(x|v)=E[S^2(x|v)]$$

  理论离散方差与理论变异函数在区域 $v$ 上的二重积分相等(Chilès & Delfiner, 1999)：
$${\sigma }^2(x|v)=\gamma (v,v)=\frac{1}{{|v|}^2}{\int }_{x\in v}{\int }_{y\in v}\gamma ||x-y||dxdy$$ 其中 $||x-y||$ 代表区域 $v$ 中的点 $x$ 和 $y$ 的距离，$|v|$ 是区域 $v$ 的面积。
  $\gamma (v,v)$ 量化了影像中一个一般的子区域 $v$ 内部的空间异质性的程度。给 $v$ 设置了几种尺寸（60，100，200，300，500和1000m）。当 $v$ 增大时，$\gamma (v,v)$ 渐渐趋近于基台值 ${\sigma }^2$。当达到基台值时，数据被完全正则化（completely regularized），也就是说再将像元合并也不会引起影像变异性的损失。（本研究的目标之一是找到使像元内空间变异性 $\gamma (v,v)$ 足够小的最大像元尺寸 $v$。因为在进行地表参数估计时希望像元内部是均质的。） 我们定义在分辨率 $v$ 下数据的正则化比率 ${TH}_v$，即影像空间变异性的损失率。它在几种空间分辨率下的值刻画了 $\gamma (v,v)$ 收敛到基台值的速度。（$T{H}_v$相当于对 $\gamma (v,v)$ 进行了标准化，以便比较）
$${TH}_v=\frac{100\gamma (v,v)}{{\sigma }^2}$$

5 结果

5.1 四种差异明显的景观的变异函数分析

  为了体现变异函数刻画景观空间异质性的能力，首先分析四种差异明显的景观的变异函数：农田(Alpilles01)，封闭灌木丛(Puechabon01)，松林(Nezer01)和热带雨林(Counami01)，如Fig1和Fig4所示。

  （整张影像 ${\sigma }^2$ 的分析）
  这四种景观的影像的空间变异性 ${\sigma }^2$ 大小排序为：农田>灌木丛>松林>热带雨林。

• 松林(Nezer01)和热带雨林(Counami01)基台值小的原因是这些区域树木茂盛且树下还有植被，因此植被覆盖的空间差异小。然而，NDVI对于茂密的植被有饱和效应，它的空间分布可能不能反应这些区域完整的植被覆盖的情况。
• 农田(Alpilles01)的高度变异性可以用高NDVI的田地（种植的冬小麦正是最绿的时候）和低NDVI的裸土（夏季作物还没有种上）的镶嵌来解释。
• 灌木丛(Puechabon01)中等的基台值是因为地中海植被结构固有的变异性（有空地，植被高度和密度的差异等），以及存在某些低NDVI的物体（道路，采石场，岩石等）与植被区域形成了对比。

  （拟合变异函数的分析）
  变异函数的形状为NDVI在影像区域内的空间结构提供了一种理解。Fig1的变异函数通过包含两种空间结构的线性模型建模得到。

• 第一种空间模式

  • 热带雨林(Counami01)的变异函数迅速增加，并且在非常短的变程内（$r_1=57m$）几乎达到了整个影像的方差（$b_1=85\%$）。因此这种景观的空间结构在这种观测尺度下无法得到较好的体现。（对于热带雨林来说所用的分辨率还是太低了？） 其他三种景观的变异函数体现出了更大的空间结构（$r_1$ 在200-300m之间）：

  • 农田(Alpilles01)的第一种空间结构（$r_1=268m$）与田地的镶嵌结构有关，这一数值与田地的平均大小（250-300 m）非常接近。第一种空间结构的方差（$b_1=60.5\%$）刻画了不同田地之间植被覆盖的差异，这一差异是这种景观空间异质性的主要原因。

  • 灌木丛(Puechabon01)的第一种空间结构（$r_1=230m$）的形状比农田更加模糊，这来源于地中海植被固有的长度尺度和景观中不一致物体的存在（采石场，道路等）。

  • 松林(Nezer01)在影像上呈现出矩形林地的空间模式，但其第一种空间模式的变程（$r_1=205m$）小于林地的尺寸（约为500m）。变异函数没有发现矩形林地结构，因为林地之间NDVI的变异性较小。第一种空间模式描述了林地内部林下叶层或种植密度的空间异质性。

• 第二种空间模式
    农田(Alpilles01)、松林(Nezer01)、灌木丛(Puechabon01)的第二种空间模式只解释了景观空间差异的一小部分（小于40%）。

  • 农田(Alpilles01)的第二种空间结构可能是土壤的差异导致的。
  • 灌木丛(Puechabon01)和松林(Nezer01)的变异函数在1500m （设定的最大变程） 之内都没有达到基台值，因此它们的第二种空间结构的变程非常不确定，因为这种结构不能被一张影像完整的包括 （回应H1和4.3节，从数值上看是 $D_c>671m$）

  这四种景观的NDVI变异函数有效地描述了它们植被覆盖的空间异质性。现在我们把分析扩展到数据集中的其他14种景观。

5.2 景观层面空间异质性的描述

  为了理解在景观层面上影响空间异质性的主要因素，我们比较这18种景观的空间异质性特征 ${\sigma }^2$ 和 $D_c$，见Fig 5和Table 4。

5.2.1 空间变异性

• 不同景观NDVI空间变异性之间的差异（${\sigma }^2$ 在0.0001到0.052之间）主要可以用景观类型解释 （主要参见Fig5）
  • 农田异质性最强（${\sigma }^2>0.02$），它的空间变异性可以通过田地之间NDVI值的差异解释，NDVI的差异是由作物的特征和状态造成的。NDVI低的裸土田地和NDVI高的种植了作物的田地镶嵌分布，尤其是当农田里同时由冬季作物和夏季作物时NDVI的空间变异性会急速增大。
  • 自然植被和森林在景观层面比农田更加同质化（大多数地点${\sigma }^2<0.008$）。森林（NDVI约为0.7）中重要的植被包括绿色林下叶层、高密度的数或阔叶，导致NDVI的分布同质化。
  • 中等植被覆盖（NDVI约为0.49，如草地）和低植被覆盖的地点（NDVI约为0.17，如热带稀疏草原和灌木带）在景观层面上是同质的。
• 然而，景观类型并不总能充分地解释其空间变异性。事实上，某些地区具有相对其景观类别反常的基台值（$0.008<{\sigma }^2<0.02$），需要结合其他因素解释。（主要参见Table2）
  • 土地利用可能会增大或减小观察区域的空间变异性。例如Concepcion03 （Table2中的分类是针叶林，与谷歌影像比对发现低NDVI区域是裸土和公路） 和Hirsikangas03 （针叶林） 中，低NDVI（小苗、裸土和水域）和高NDVI区域（森林）的对比增大了这些景观的空间变异性。相反地，在Gilching02 （农田和混合森林，Fig5中的六角星） 中一片同质化的森林区域可以解释它的基台值比其他农田低的原因。
    
  • 其他环境因素诸如水域（Laprida01， Larose03）或土壤性质的变化（Laprida01的土壤盐分）也可能是NDVI空间变异性的原因。

5.2.2 平均长度尺度

  影像中空间结构的平均大小由平均长度尺度（mean length scale） $D_c$ 来量化。它与景观类型没有很强的依赖性。在所研究的18个地区中，$D_c$ 的取值范围是215-1059m。为了阐述清晰，接下来对每种景观类型分别讨论。

• 农田
    $D_c$ 的取值范围是289-664m。它们的长度尺度完全包含在影像的尺度内（$D_c<D_{c,T,3km}$ ）。田地之间NDVI的不连续性造成了一种镶嵌状的空间结构，这一现象主要是源于人为的过程（例如田地大小、播种和收获的时间、作物轮作等）。农田空间结构的大小主要受到田地大小（field size）的影响。大块田地（Fundulea01, $D_c=619m$； Barrax02, $D_c=514m$）和小块田地（Alpilles02，$D_c=289m$）各自不同。然而通过变异函数的变程刻画田地大小不是一个简单的问题。事实上，这会受到具有相似NDVI值的多块田地的干扰，从而制造出更大的空间结构（例如SudOuest02的第二变程）。因此，NDVI的空间结构不仅取决于景观本身的属性，也依赖于季节。例如，Barrax03的影像拍摄于七月，灌溉盘之外都是裸土。由于NDVI描述的是植被覆盖，它的变异函数探测到的是灌溉盘内的空间结构而不是裸土的空间结构。如果拍摄日期不同获取的结果可能会不同。除此之外，辐射度量值可能会比NDVI对土壤属性更敏感（例如红色反射比可能会更好地捕获裸土区域的空间结构）。
  
• 自然植被和森林
    自然植被和森林的 $D_c$ 值在215-1059m。它们的空间结构模式更模糊，NDVI的不连续性比农田更平滑。主要原因包括生态过程（光截获、物种竞争）、地貌因素（微地形、岩区）、季节（林下叶层的变化）和人类因素（森林采伐、草场放牧强度）。
  Gilching02由森林和农田组成，它的变异函数尤为有趣，因为它同时表现了森林和农田的空间结构。第一变程（$r_1=525m$）解释了影像92%的方差，可能体现了田地的结构。更大的变程（$r_2=1150m$）可能与森林内部更模糊的空间结构有关，在影像尺度上的变化更小。这个例子说明田地的空间结构是农田景观层面NDVI变异性的主要原因。
  与农田相对的是，自然植被和森林的特征并不总能在影像尺度表达出来。事实上，NDVI的 $D_c>D_{c,T,3km}$ （Fig5中虚线右侧的6个点，这里主要指热带草原、稀疏植被和针叶林） 的影像不够大，不能完全包含景观的变化。这是由于存在相对于影像范围的大型空间结构，这一结构是影像变异性的重要组成部分。此类的部分地区第二变程超过了 $d_{max}$。如前文所述，二阶平稳性假设 （H1） 在这些地区不成立。Puechabon01影像边缘对比强烈的空间结构可能是其在现有影像尺度不满足平稳性的原因。
• 总结
    这一分析为景观层面空间异质性的影响因素提供了一些理解。变异函数基台值（${\sigma }^2$）是一种景观变异性指标。景观层面空间异质性的主要原因是人类活动导致的土地利用的变异性。农田是异质性最强的地区。田地的空间结构是景观空间变异性的重要来源。自然植被和森林在景观层面同质性更强。本研究中用 $D_c$ 度量典型的长度尺度，取值从215到1059m不等。这一信息将在讨论部分被用于定义可以识别景观主要空间变异性的传感器空间分辨率。

5.3 将空间异质性表达为空间分辨率的函数

  Fig 6展示了同一地点（Puechabon01）的几种NDVI经验变异函数图，分别由拍摄日期接近的Ikonos（4m）、SPOT-HRV（20m）和SPOT-HRV聚合到低分辨率（60，300，500，1000m）得到。变异函数随空间分辨率的变化描述了数据正则化（data regularization）对空间异质性的影响。基台值的降低标志着空间分辨率降低时空间变异性的损失。Ikonos和SPOT-HRV影像NDVI变异性的差异强调了这两种传感器发现了不同的景观特征。Ikonos影像上有更多纯植被或岩石像元，因此增大了NDVI变异性。其他如仪器的PSF、拍摄时间之间大气条件的差异等因素可以解释这两种影像基台值之间的差异。
  此外，当空间分辨率降低时变异函数变得更加规则。4、20、60m分辨率的变异函数有相似的形状，且探测到了同样的长度尺度。但对超过300m的分辨率，230m的第一长度尺度没有被变异函数探测出来。因此，当传感器像元比长度尺度大时空间结构不能被提取。

  本研究的一个重要问题是如何量化低分辨率下空间异质性的损失。
  如4.4节所述，离散方差 $\gamma (v,v)$ 代表了低分辨率像元内部空间异质性的程度。 Fig 7中，离散方差随着像元尺度增大的速率可以用Fig 5中展示的景观异质性特征解释。例如空间结构弱的Counami01（热带雨林，$D_c=215m$），NDVI数据在300m分辨率时几乎已经完全正则化了（见Table 4，$T{H}_{300}=88\%，T{H}_{1000}=97\%$）（TH值可以理解为Fig7中的 $\gamma (v,v)$ 除以 ${\sigma }^2$，起到标准化的效果，以便对${\sigma }^2$不同的图进行比较，代表像元内异质性丢失的比例） 。但Alpilles01（农田）正则化不那么明显（见Table 4，$T{H}_{300}=49\%，T{H}_{1000}=82\%$），因为存在一种长度尺度更大的空间结构（$D_c=664m$）。
  在给定空间分辨率下，离散方差 $\gamma (v,v)$ 取决于 ${\sigma }^2$ 和影像长度尺度与像元大小的比例。 例如在500m分辨率下，尽管Alpilles01（农田，$r_1=268,r_2=1290$）的 ${\sigma }^2$ 比Fundulea01（也是农田, $r_1=781$，没有$r_2$ ）小，但离散方差更大 。同样在这一空间分辨率下，由于Fundulea01的长度尺度 （781 m）仍然能被传感器探测出来（$T{H}_{500}=47\%$）, 像元内部的异质性低。然而对Alpilles01，由于与第一长度尺度（$r_1=268m$）有关的重要空间变异性（$b_1=61\%$）在500m分辨率下丢失了（$T{H}_{500}=64\%$），其500m像元的离散方差大 （Alpilles01的第一种空间结构被包含在一个像元以内了，因此像元内空间异质性更大） 。在1000m分辨率下，Fundulea01的长度尺度不再能被传感器探测出来。因此由于Fundulea01的 ${\sigma }^2$ 比Alpilles01大，Fundulea01的 $\gamma (v,v)$ 更大。
   （总结） 变异函数是量化中分辨率下数据正则化的有效工具。像元内部的空间异质性由离散方差 $\gamma (v,v)$ 量化，并且随着像元的增大迅速增大，直到像元大小超过数据的平均长度尺度，随后渐渐趋近于变异函数的基台值 ${\sigma }^2$ 。 通常而言，本研究使用的数据在1000m分辨率下离散方差的差异可以完全用影像变异性 ${\sigma }^2$ 来解释。

6 讨论

  本节将讨论空间异质性的刻画引发的几个问题。

6.1 影像大小对刻画空间异质性的影响

  我们提出了一种指标来评估3000m影像是否足够刻画和量化空间异质性：平均长度尺度 $D_c$ 必须小于阈值 $D_{c,T,3km}=671m$ 。实验变异函数出现平稳的基台值是一个必要条件，但不是充分条件。前文分析展示了对某些景观，由于存在一种相对于影像较大的空间结构，NDVI空间变异性不能被一张影像完全包含进来。
  Fig 8展示了增大影像大小对实验变异函数的影响。Alpilles02的 $D_c=289m$，比 $D_{c,T,3km}=671m$ 小，变异函数的形状没有随影像大小增大而改变。对Puechabon01而言，在3000×3000m影像上它的实验变异函数在 $d_{max}=1500m$ 之前没有到达基台值：它的 $r_2=1750m$，$D_c=841m$，超过了 $D_{c,T,3km}$。但在7000×7000m的影像上（与之前的影像中心相同）它的变异函数达到了基台值。在这一地区和某些其他地区，增大影像尺寸能够更好地刻画景观的长度尺度。然而，这一规则可能不能适用于所有地区，因为增大影像范围可能会包含进新的空间结构。如果新空间结构的尺寸仍相对于新影像范围较大，它们还是不能完全被变异函数所描述。
  因此，空间异质性的刻画很大程度上依赖于观测尺度。$D_c$ 标准可以作为一种经验性的指标来检查所使用的影像是否足够大，使变异函数能够恰当地量化景观长度尺度。这一标准适用于本工作研究的大多数地区，尤其是异质性的农田地区。

6.2 用变异函数刻画局部空间异质性

  空间变异性（${\sigma }^2$）和平均长度尺度（$D_c$）代表了影像尺度上的平均空间异质性特征。为了测试它们在局部尺度上是否有代表性，将原始影像划分成9个1000×1000m的子区域分别进行建模。在Alpilles02（小块田地）和Fundulea01（大块田地）上研究局部空间异质性特征。Fig 9展示了每个子区域 $v_i,i=1,...,9$ 的局部平均长度尺度 $D_{c,v_i}$ 相对于影像的平均长度尺度 $D_c$ 的比例：
$$D_{c,v_i,rel}=100\frac{D_{c,v_i}-D_c}{D_c}$$ 和局部基台值 ${\sigma }_{v_i}^2$ 相对于影像基台值 ${\sigma }^2$ 的比例：
$${\sigma }_{v_i,rel}^2=100\frac{{\sigma }_{v_i}^2-{\sigma }^2}{{\sigma }^2}$$

• 对于Fundulea01，田地空间结构的大小与1000m子区域接近。因此，九个子区域计算出的特征之间有很强的变异性（variability），并且子区域的平均值和整幅影像的特征之间也有很大差异（如Fig 9 左图）。
• 相对地，对于Alpilles02，影像的变异函数能够很好地刻画局部空间异质性的特点。第一变程（184m）像子区域 $v_7(D_{c,v_7}=170m)$一样描述了小的田地，第二变程（410m）解释了如子区域 $v_9(D_{c,v_9}=370m)$ 的较大的田地。除此之外，子区域的平均空间异质性特征与整幅影像的特征很接近（Fig 9 右图）。

  （总结） 由于变异函数提供的是影像尺度的平均空间异质性特征，这一信息足够从全局角度量化影像子区域的平均空间异质性。（Fig9右图） 然而对局部空间变异性的描述是有限的 （Fig9左图） 。应该使用小波分析等其他方法来探测和量化局部变异性和影像内部的局部长度尺度。

6.3 确定足够的像元大小以表达景观长度尺度

  本研究使用参数 $D_c$ 量化景观的典型长度尺度，这一变量可以用于计算各幅影像能保留主要的空间变异性的最大像元尺度。香农的信息论认为信号合适的采样频率必须超过其最大频率的两倍（Shannon，1984）。我们已经看到，如果假设数据按规则格网分布，那么 $D_c$ 是能分开“等效独立数据”的距离。根据香农的理论，空间采样的频率必须超过 $2/D_c$ 。因此，为了保存大部分的NDVI空间变异性，像元尺寸必须小于 $D_c/2$ 。 将这一标准应用于本研究的数据，$D_c/2$ 的范围在108-530m以内，平均为324m。我们的结果受限于研究数据的数量和特征，但计算出的 18张VALERI数据库影像的 $D_c/2$ 最小值 （108m） 可以作为表达景观尺度植被覆盖的大部分空间变异性的像元尺寸上界的一个指标。
  （总结） 为了限制空间异质性对使用遥感影像进行地表参数非线性估计过程的影响，像元的尺寸必须足够小以表达数据的空间变异性和最小化像元内部的变异性。参数 $T{H}_v$ 度量了在给定分辨率 $v$ 下NDVI空间变异性相对于20m分辨率的损失率。对于本工作研究的大部分地区，100m分辨率下NDVI空间变异性的损失 ${TH}_{100}$ 少于30%（见Table 4）。此外， ${TH}_{100}$ 较大的是在景观尺度具有同质性的地区，因此它们像元内部的异质性 $\gamma (v,v)$ 在100m分辨率下很小。（空间变异性绝对数值本来就小，相对比例大点也没事的意思？） 因此，对于本研究使用的数据，取100m作为像元大小可以减少像元内部空间异质性对非线性估计过程造成的偏差。

7 结论

  本文说明了高分辨率NDVI的变异函数建模是刻画和量化景观级别植被覆盖空间异质性的有力工具。变异函数基台值 ${\sigma }^2$ 度量了景观的变异性。影像的空间结构由变程 和各变程在总方差中的占比共同刻画 （$r、b$） 。另外，我们介绍了变程积分 （$A$） 的概念，它将变异函数模型全部的结构参数总结为一个特征区域。它的平方根 （$D_c$） 是各变程的加权和，量化了数据的平均长度尺度，例如影像空间结构的平均大小。变程积分可以用来判断影像大小是否足以用变异函数合理地度量长度尺度。我们提出它必须小于影像面积的5%。因此对于3000×3000m的影像，其变程积分的平方根必须小于671m。
  对十八种景观NDVI变异函数的建模凸显了土地利用在景观层面上对植被覆盖空间异质性的影响。异质性最大的地区（${\sigma }^2$在0.02到0.05之间）是农田，其田地的空间结构解释了大部分的（60%~100%）NDVI空间变异性。自然植被和森林在景观尺度更加同质化（${\sigma }^2$在0.0001到0.02之间）。然而，特殊物体的存在可能会增加其空间变异性。它们的平均长度尺度在216-1060m。它是人类活动、生态功能和气候等多种过程共同造成的。需要注意的是，3000m影像对于某些景观长度尺度的度量太小了。7000×7000m范围可能会更合适。
  变异函数建模是刻画传感器捕获的空间变异性随空间分辨率降低而损失的有效方式。当空间结构的长度尺度比像元还小时它的空间结构不能被传感器表达出来。通过变异函数模型计算出的理论离散方差 （$\gamma (v,v)$） 量化了像元内的空间异质性。它随着像元的增大迅速增大，直到像元达到平均长度尺度，随后收敛于变异函数的基台值 ${\sigma }^2$。离散方差可以被用作修正遥感数据地表参数非线性估计的额外知识。然而，需要找到无需同时相高分辨率影像就能获得像元内空间异质性评价指标的方法，因为如果已经有了高分辨率影像中分辨率影像就没有用了。一种可能的方法是时间采样或对每种景观的高分辨率影像做空间采样。为了测试这一方法，需要在更大的高分辨率时空数据上做实验。
  最后，我们根据变异函数得出的平均长度尺度提出了刻画景观植被覆盖主要空间变异性的像元尺寸上界 （$D_c/2$） 。根据对VALERI数据库中18种景观的分析，这一数值约为100m。由于对同质性的景观而言，这一空间分辨率下NDVI空间变异性的损失很小，像元内异质性对非线性估计过程造成的偏移能够减小。然而，这一结论受限于本文分析的景观的特征和数量。需要对各景观类型进行更有代表性的采样使这一数值更精细化。由于变异函数描绘了影像空间异质性的平均特征，在某些场景下它对长度尺度的量化可能是粗糙的。需要引入小波分析等其他工具来量化影像更精细的局部空间尺度。
  最佳像元尺度的定义不是一个小问题。它主要取决于研究目标、观测对象和观测技术。首先，像元大小必须足够大以于目标物体和观测技术相匹配。如果空间分辨率太高，小尺度空间结构可能会妨碍表面特征的获取。需要进一步研究来确定能够刻画景观尺度植被覆盖像元大小的合适下界。其次，像元大小必须足够小以捕获数据的空间变异性并最小化像元内的空间变异性。本文提出的像元尺度是在景观层面描述植被覆盖的合适的像元尺寸的上界。最佳像元尺度应该在两个界限之间，但是在未来的对地观测任务中也需要考虑技术和经济的限制等其他因素。

8 补充材料

8.1 二阶平稳性假设、内蕴假设

参考：插值、平稳假设、本征假设、变异函数、基台、块金、克里格、线性无偏最优…地学计算概念及公式推导

8.2 Punctual与Regularized

• 点观测值的变异函数称为punctual variogram。
• 区域平均观测值的变异函数称为regularized variogram。
    regularization是指传感器聚合了一个区域内的信息产生平均观测值的过程。二者的关系：Regularized variogram可以看作是由其背后基于点观测值的punctual variogram在不同聚合尺度下衍生出的。
    遥感影像的regularization area是传感器瞬时视场角（IFOV）和点扩散函数（PSF）决定的。像元可以看作regularization unit。对同一punctual variogram，随着聚合分辨率的降低，regularized variogram的基台值逐渐降低，当分辨率不是特别低时形状仍与punctual variogram相似，但当分辨率变得很低时regularized variogram形状变得很简单（如本文Fig 6和下图所示）。
    本文把20m分辨率的SPOT影像的变异函数作为punctual，不同聚合尺度的VALERI影像的变异函数作为regularized，以此度量不同像元尺度下空间异质性的损失程度。
  

本节参考文献：

• Woodcock, C. E., Strahler, A. H., & Jupp, D. L. B. (1988). The use of variograms in remote sensing: A/ scene models and simulated images. Remote Sensing ofEnvironment, 25, 323–348.
• Jupp, D. L. B., Strahler, A. H., &Woodcock, C. E. (1988). Autocorrelation and regularization in digital images: I. Basic theory. IEEE Transactions on Geoscience and Remote Sensing, 26, 463–473