引用本文
于文涛, 李静, 柳钦火, 曾也鲁, 尹高飞, 赵静, 徐保东. 中国地表覆盖异质性参数提取与分析.地球科学进展, 2016, 31(10): 1067-1077
Yu Wentao, Li Jing, Liu Qinhuo, Zeng Yelu, Yin Gaofei, Zhao Jing, Xu Baodong. Extraction and Analysis of Land Cover Heterogeneity over China. Advances in Earth Science, 2016, 31(10): 1067-1077  
Doi:10.11867/j.issn.1001-8166.2016.10.1067
Permissions
Copyright©2016, 地球科学进展 编辑部
中国地表覆盖异质性参数提取与分析
于文涛1,3, 李静1,2*,*, 柳钦火1,2, 曾也鲁1,3, 尹高飞1,3, 赵静1,2, 徐保东1,3
1.中国科学院遥感与数字地球研究所 遥感科学国家重点实验室,北京 100101
2.全球变化研究协同创新中心, 北京 100875
3.中国科学院大学 资源与环境学院,北京 100049
*通信作者:李静(1978-),女,黑龙江齐齐哈尔人,副研究员,主要从事植被辐射传输模型,叶面积指数反演等研究.E-mail:lijing01@radi.ac.cn

作者简介:于文涛(1995-),男,安徽五河人,硕士研究生,主要从事复杂地表的叶面积指数反演方法研究.E-mail:1096392329@qq.com

基金:国家自然科学基金项目“非均质混合像元遥感反射波谱模型构建及叶面积指数反演方法研究”(编号:41271366); 国家重点基础研究发展计划项目“复杂地表遥感信息动态分析与建模”(编号:2013CB733401)资助
摘要

地表异质性广泛存在于陆地表面各个尺度,是地表参数遥感反演不确定性的主要来源之一。基于高分辨率地表分类参考图,提取出低分辨率混合像元的端元数量和边界长度指标来描述地表异质性。然后以中国地区为例,使用全国30 m空间分辨率GlobalLand 30地表分类数据集提取出1 km尺度像元的描述混合结构和破碎程度的异质性指标。并基于提取出的异质性指标分析了中国区域在1 km尺度上非均质地表地物类型的组合特征、斑块特征和不同生态群系内部异质性特征。发现山地和生态交错区是主要的高异质性区域,稀树草原生物群系内部异质性最大(平均边界长度为7 426 m),其次依次为森林(4 323 m)、耕地/草地(3 160 m)和灌丛(1 779 m)。

关键词: 空间异质性; 地表参数反演; 地表分类图; 辐射传输
中图分类号:P237 文献标志码:A 文章编号:1001-8166(2016)10-1067-11
Extraction and Analysis of Land Cover Heterogeneity over China
Yu Wentao1,3, Li Jing1,2,*, Liu Qinhuo1,2, Zeng Yelu1,3, Yin Gaofei1,3, Zhao Jing1,2, Xu Baodong1,3
1.State Key Laboratory of Remote Sensing Science, Institute of Remote Sensing and Digital Earth, Chinese Academy of Sciences, Beijing 100101, China
2.Joint Centre for Global Change Studies (JCGCS), Beijing 100875, China
3.College of Resources and Environment, University of Chinese Academy of Sciences, Beijing 100049, China

First author:Yu Wentao(1995-), male, Wuhe County, Anhui Province, Master student. Research areas include LAI inversion methods of complex land surface.E-mail:1096392329@qq.com

*Corresponding author:Li Jing(1978-), female, Qiqihar City, Helongjiang Province, Associate Professor. Research areas include vegetation radiative transfer model and LAI inversion.E-mail:lijing01@radi.ac.cn

Fund:Project supported by the National Natural Science Foundation of China “Construction of remote-sensed reflectance spectrum model and research on LAI inversion method over inhomogeneous mixed pixel”(No.41271366); The State Key Development Program for Basic Research of China “Dynamic analysis and modeling of remote sensing information over complex surface” (No.2013CB733401)
Abstract

Spatial heterogeneity exists in land surface at every scale, and it is one of key factors to bring uncertainty to land parameter retrieval from remote-sensed data. This paper proposed a methodology to use the boundary length among different land cover types to characterize and quantify land surface heterogeneity based on high-resolution land cover images. Then the heterogeneity feature at 1 kilometer scale in China was extracted from “GlobalLand30” land cover datasets with the spatial resolution of 30 m. The mixed structure, degree of fragmentation and intra-heterogeneity of eight main vegetation biomes from MODIS land cover product over heterogeneous surface in china were analyzed. Mountain area and ecotone are more heterogeneous than other regions. Savanna biome (average boundary length is 7 426 meters) is the most heterogeneous zone followed by forest, grass/crop and shrub biome with average boundary length of 4 323, 3 160, 1 779 meters, respectively.

Keyword: Spatial heterogeneity; Parameter inversion; Land cover dataset; Radiative transfer.
1 引言

空间异质性广泛存在于陆地表面的各个尺度, 是描述地表场景复杂程度的一个重要指标。大区域范围地表植被特征参数(如叶面积指数、植被覆盖度等)的获取及应用(如碳储量、作物产量估算等)常依赖于遥感技术。目前地表植被参数遥感反演中使用的遥感影像的空间分辨率从几十米到几千米。在此分辨率下, 一个像元表示的地表范围通常为多种地物的混合, 像元内部具有异质性。而目前全球尺度的地表植被参数产品反演算法都不考虑像元的混合特性, 均假设地表为单一植被类型, 使用特定结构假设的植被模型进行反演, 如MODIS/LAI和CYCLOPES/LAI产品等。这种方法在均匀地表是合理的, 可以保持较高精度, 而在异质性大的区域会带来很大的不确定性。因此, 研究并构建对地表植被辐射传输过程有影响的异质性特征并对其参数化, 进而对其在地表的分布特征进行分析, 对地表参数遥感反演、参数产品精度评价以及农业产量估算、碳储量估算等具有重要意义。

目前, 人们对于异质性仍没有一个普适性的定义, 而是与研究的目标和关注的尺度相关。Hintz等[1]认为地表异质性可以表现为信息总量、内部相关性或地表参数的各向异性。近年来, 相关学者从不同角度对地表异质性进行了研究和量化。Haralick等[2]通过提取经验指数描述图像的纹理特征来表征地表异质性。De Cola[3]从分形的角度描述地表异质性。Garrigues等[4]基于变异函数模型从空间变异性和空间结构2个方面对地表异质性进行了定量化。Hintz等[1]引入了信息熵光谱的概念来确定哪个尺度上的异质性对地— 气交互影响最重要。这些模型均是从地理统计学的角度分析地表的异质性, 从地表不同地物间的物理散射机理角度考虑异质性的模型不多。在遥感辐射建模和叶面积指数反演研究中, 通常认为空间异质性包括不同属性地块拼接造成的地表不连续性(Discontinuity)和地块内部属性具有的密度变化(Density change)2个方面[5, 6], 即从地表的景观格局出发, 将空间异质性划分为2个尺度, 在景观尺度上认为地表由属性相对均一的斑块相互拼接而成, 到斑块尺度上又认为斑块内部属性在空间上具有一定的抖动。这种划分方法主要源于地表异质性对植被辐射传输过程的影响。

地表的异质性特征对混合像元的辐射传输模型的模拟精度有着至关重要的影响。Gilabert等[7]和Liang[8]指出混合像元的双向反射分布函数(Bidirectional Reflectance Distribution Function, BRDF)满足线性加权规律。然而, Qin等[9]发现, 线性混合规律并非在所有复杂场景和尺度下都适用, 景观尺度上各景观斑块的BRDF近似满足线性加权规律, 当尺度降低到单株树时, 线性加权会给场景的BRDF模拟带来较大的误差。线性混合的不确定性主要来自于斑块之间的非线性交互作用, 即交叉辐射[10]。Marshak等[11]发现, 由于交叉辐射的存在, 辐射场的空间分布比介质的空间分布更平滑, 这一现象被称为“ 辐射平滑” (Radiative smoothing)。Kobayashi等[12]在模拟异质性对冠层中辐射传输过程的影响时发现, 忽略交叉辐射会造成光合有效辐射截获量(Fraction of Absorbed Photosynthetically Active Radiation, FAPAR)的低估, 极端情况下甚至会达到52%。交叉辐射的影响强度与斑块大小密切相关, 斑块越破碎, 即斑块间的边界越长, 边界处交叉辐射的影响越大[11, 13, 14], 而当斑块大小达到某特征尺度时, 交叉辐射的影响可近似忽略, 此时线性加权规律近似成立, Widlowski等[10]通过计算机模拟发现, 该特征尺度在百米级左右。因此, 当地表结构复杂, 斑块平均大小在百米或百米级以下时, 必须考虑交叉辐射的影响。Zeng等[15]考虑混合像元的边界长度和边界高差, 针对典型的农林交错带, 构建了考虑边界交叉辐射的农林交错带二向反射模型(RTAF模型)。

地表异质性也是目前中低分辨率定量遥感产品误差的主要来源之一。在叶面积指数反演研究中也有指出, 相比于斑块内的不均一, 斑块混合对叶面积指数反演精度的影响更大[16]。Chen[5]在分析遥感反演叶面积指数的尺度效应时发现, 当像元存在2个或多个差异很大的地块拼接时, 反演值会出现很大的负偏差, 在森林和水体的混合场景中, 偏差最高可以达到真值的45%。Tian等[17]通过聚合1 km分辨率AVHRR数据分析尺度效应时发现, 低分辨率影像叶面积指数(Leaf Area Index, LAI)的反演误差和像元内主要地表覆盖类型的比例呈负相关。Fang等[18]通过对MODIS 地表覆盖和LAI产品的统计发现, 地表覆盖分类不准会导致稀树草原群系的LAI过估计, 而在森林群系出现欠估计。可见, 像元内不同类型斑块间的拼接特征是对异质性地表的辐射传输建模及叶面积指数反演研究有重要影响的典型性特征。

综上所述, 在异质性地表中斑块边界是描述边界交叉辐射的一个重要参数, 也是决定像元的异质性特征是否对辐射传输过程带来较大影响进而使线性混合规律不再满足的一个重要衡量指标。目前针对辐射传输过程的像元异质性没有定量化的描述方法, 需要对其进行参数化并分析其空间分布规律。本文提出基于地物种类数和边界长度的地表异质性特征参数化方法。目的是为非均质地表的建模及反演提供地表异质性特征的定量化描述, 为进一步提升异质地表场景的参数反演、产量估算、碳汇估算等的精度提供辅助信息。本文以中国区域为研究区, 使用30 m分辨率的地表参考分类图提取并分析了1 km空间分辨率的像元异质性特征, 并简要说明了地表异质性信息的应用方法与潜力。

2 地表覆盖异质性特征描述方法

本文中提出的地表异质性特征描述方法是面向地表的辐射传输过程。如上所述, 像元内斑块拼接对异质性地表的辐射传输建模及叶面积指数反演研究有重要影响。像元内斑块特征明显, 不同端元的边界长, 则边界处的交叉辐射明显; 同时, 当边界处相邻2种植被或非植被的高差明显时, 边界处的交叉辐射更加强烈。在这些条件下, 需要考虑边界处的交叉辐射, 简单的线性混合模型不再适用并会带来较大误差, 需要对其进行参数化, 进而对不同的异质性程度进行区分, 以提高建模和反演的精度。本文的异质性特征描述方法主要考虑斑块混合特征, 也是基于边界长度和边界高差以及端元种类数。由于分析显示斑块内部的不均一性对辐射传输过程和叶面积指数反演的精度影响不大, 因此本文假设斑块内部是均一的。

2.1 端元信息

混合像元在遥感影像中普遍存在, 端元信息可以准确地描述混合像元的组合特征。端元信息包括端元种类数和端元丰度。端元种类数是指混合像元内所包含的端元个数。端元丰度指各端元面积所占混合像元总面积的比例。基于端元种类数和端元丰度, 可以准确地分析混合像元的混合特征。

端元个数越多, 混合像元的斑块特征会越明显。当端元个数少时, 即植被类型少, 同一植被类型可能以多个斑块特征存在, 也可能具有明显的斑块特征。因此, 端元种类数是一个辅助的描述参数。端元种类数一定程度上描述了区域地表的复杂程度, 而像元内的斑块特征是否明显还需要边界信息来描述。

2.2 边界信息

边界长度指不同植被类型的斑块间形成的边界的长度, 是描述边界信息的一个最主要参数。边界长度长, 则斑块特征明显; 边界长度为0时, 说明像元内为同种地物类型, 即纯像元。在考虑边界长度时, 边界处植被高差和边界朝向对地表的辐射传输过程也有较大影响。当边界处相邻2种植被或非植被的高差明显时, 则边界处的交叉辐射会更强烈。如果边界延伸方向与太阳入射的主平面方向平行, 则边界处的交叉辐射弱; 反之, 如果边界延伸方向与太阳入射主平面垂直, 则边界处的交叉辐射最为强烈。

端元信息和边界信息描述的特征相对独立, 混合像元中端元种类数少, 端元边界不一定少, 两者是不相关的。在异质性地表的辐射传输建模及叶面积指数反演研究中, 端元信息可以提供是否包含对反演有严重影响的地物类型及其风度, 比如水体; 对非均质地表建模, 端元信息量化了模型中的线性(一阶)部分。边界信息主要反映斑块的特征, 包括斑块是否破碎、斑块拼接处的相互遮挡及阴影效应, 量化了遥感反演建模中的高阶部分。端元信息和边界信息相辅相成共同描述了混合像元对异质性地表建模及反演有重要影响的特征性信息。

3 研究区与数据

本文以中国区域为研究区, 分析1 km分辨率的像元异质性特征。1 km像元内的地表类型分布采用GlobalLand30地表覆盖分类产品。

3.1 GlobalLand30地表覆盖分类产品

GlobalLand30-2010是以2010年为基准年生产的30 m分辨率全球地表覆盖分类产品[19]。该数据集以Landsat ETM/TM影像作为原始数据, 联合国产减灾卫星(HJ-1)和北京一号小卫星(BJ-1)等资料, 将全球陆地地表分为耕地、森林、草地、灌木地、湿地、水体、苔原、人造地表、裸地、冰川和永久积雪10个类别。根据第三方验证结果, 该产品的总体分类精度达到83%[20]。根据产品的条带编码规则, 挑选出中国区域的地表覆盖影像进行下一步处理。

3.2 MODIS Land Cover产品MCD12Q1

为分析1 km尺度上不同植被类型内部的异质性分布特征, 本文采用了MODIS 2010年的地表覆盖图(MCD12Q1)来区分地表的植被类型。MCD12Q1是依据Terra星和Aqua星的MODIS数据生产的全球地表覆盖产品, 该产品的空间分辨率为500 m, 使用过程中使用最近邻法重采样到1 km分辨率。时间分辨率为1年, 总体的分类精度约为75%[21]。MCD12Q1产品共有5套分类系统, 本文分类系统采用针对LAI和植被冠层吸收的光合有效辐射比率(Fraction of Photosynthetic Active Radiation, FPAR)反演的第三套分类系统。该分类系统将自然地表分为11个类别, 其中植被类型有8个群系, 分别为草地/谷物作物、灌木、阔叶作物、稀树草原、常绿阔叶林、落叶阔叶林、常绿针叶林和落叶针叶林(数据来源:http:∥e4ftl01.cr.usgs.gov/MOTA/MCD12Q1.051/2010.01.01)。该数据产品与GlobalLand30的分类系统在理想状况(即分类准确)下的对应关系见表1。MODIS产品分类结果是像元内的最主要地物类型, 而像元内可能存在的其他地物类型(内部异质性)被忽略。下文中将用30 m的分类数据对不同的植被群系内部异质性进行分析。

表1 MODIS分类系统和GlobalLand30分类系统对应关系 Table 1 Correspondence between MODIS classification system and GlobalLand30 classification system
4 结果分析与讨论
4.1 中国地表异质性特征分布图

中国幅员辽阔, 地表类型复杂多变, 不同的生态景观呈现出不同的异质性表现。图1和图2分别为1 km尺度中国区域端元种类数分布和边界长度分布图。在本文分析中, 不针对特定的太阳入射及观测的几何条件, 因此, 边界长度是指所有边界高差和朝向的边界之和。分析比较发现, 两特征参量生成的产品总体上具有很高的一致性。端元种类数为1~7, 因此, 其分布图体现细节有限。相比之下, 图2能够体现更多细节。图2中显示, 异质性较大的区域主要集中在山区, 包括西南部山地区域、大兴安岭、太行山脉、武夷山脉等。在平原、盆地等区域, 异质性特征不明显。

图1和图2中显示出较明显的块状特征, 即明显的不连续的边界。这主要由分类图决定。本文采用30 m的GlobalLand30地表覆盖分类图产品。产品是基于Landsat/TM和ETM产品生产, 由于天气原因及时间分辨率限制, 不同区域通常基于不同时相的TM数据分类, 导致分类结果的边界不连续。

图1 1 km尺度中国混合像元端元种类特征参量图Fig.1 China endmember types of mixed pixel image at 1 km resolution

图2 1 km尺度中国边界长度特征参量图Fig.2 China boundary length image at 1 km resolution

4.2 非均质地表的地物类型混合特征分析

中国的异质性分布特征显示, 山区是比平原、盆地、高原等异质性更大的区域。那么是什么混合方式导致了山地的高异质性?本节主要分析非均质地表的地物类型混合特征。图3展示了中国区域所有陆地区域和仅植被覆盖区域的地物种类直方图。仅植被覆盖区域根据MODIS分类图判断。通过图3可以看出, 中国区域1 km单位像元中, 纯净像元和2种地物混合的像元所占比例最高, 分别为34%和33%; 其次为3种地物混合的像元, 约占22.5%; 前3种占89.5%。对于MODIS分类图判断得到的植被区域, 所占比例最高的为2种地物的混合像元, 约为36%; 其次为纯净像元和3种地物混合像元, 分别占23.6%和27.4%; 前3种占87%。4种地物和5种地物混合的像元是构成异质性较大区域的主要类型。图3b显示目前1 km分类图中标为植被的区域, 根据30 m分辨率分类图的分析, 仅23.6%为单一植被类型, 多数都为混合像元。为了解1 km分辨率的植被像元的混合特征, 下面将分别分析在MODIS分类图中分为植被的区域的不同混合种类数对应的植被类型的组合特征, 即主要是由哪些植被类型混合得到的。

图3 中国区域1 km像元端元种类数直方图分布
(a)植被与非植被区域; (b) 植被区域Fig.3 Frequency histograms of endmember types in china at 1 km pixels
(a) Vegetation and non-vegetated region; (b) Vegetation zone

为便于分析, 同时考虑地表类别的反射特征和植被的高度特征, 将端元的类型由GlobalLand 30数据集的10种地表覆盖类型简化为7种, 具体的对应关系见表2。

表2 端元类型与地表覆盖类型对应表 Table 2 Look-up table of endmember type and land cover type

4.2.1 端元种类数为2

对于2种端元的混合像元, 共有21种混合方式。以不同的组合方式作为横轴, 出现的频率作为纵轴, 绘制频率直方图如图4a所示。横轴代表的组合形式依次为:AB, AC, AD, AE, AF, AG, BC, BD, BE, BF, BG, CD, CE, CF, CG, DE, DF, DG, EF, EG, FG。比例最高的组合类型依次为森林— 草地的29.33%, 耕地— 非植被的16.42%和草地— 裸地的13.70%。其中比例超过8%的组合形式还有森林— 耕地、灌丛— 草地、耕地— 草地。其余的组合形式所占比例均低于3%。此外, 有水体参与混合的比例为4.72%。

4.2.2 端元种类数为3

对于3种端元的混合像元, 共有35种混合方式。以不同的组合方式作为横轴, 出现的频率作为纵轴, 绘制频率直方图如图4b所示。横轴代表的组合形式依次为ABC, ABD, ABE, ABF, ABG, ACD, ACE, ACF, ACG, ADE, ADF, ADG, AEF, AEG, AFG, BCD, BCE, BCF, BCG, BDE, BDF, BDG, BEF, BEG, BFG, CDE, CDF, CDG, CEF, CEG, CFG, DEF, DEG, DFG, EFG。比例最高的组合类型为森林— 耕地— 草地, 占56.88%, 超过3种地物端元混合比例的一半, 比例相对较高的有森林— 灌丛— 草地, 占7.75%, 灌丛— 草地— 裸地, 占5.97%。其他的组合形式占比均小于5%。其中包含水体的混合像元的比例为11.85%。

图4 异质地表地物组合类型频率直方图Fig.4 Frequency histogram of feature combination types over heterogeneous land surface

4.2.3 端元种类数为4

对于4种地物混合的像元, 7种端元类型35种混合方式。横轴为不同的组合形式, 纵轴为出现的频率, 绘制直方图如图4c所示。横轴代表的组合形式依次为:ABCD, ABCE, ABCF, ABCG, ABDE, ABDF, ABDG, ABEF, ABEG, ABFG, ACDE, ACDF, ACDG, ACEF, ACEG, ACFG, ADEF, ADEG, ADFG, AEFG, BCDE, BCDF, BCDG, BCEF, BCEG, BCFG, BDEF, BDEG, BDFG, BEFG, CDEF, CDEG, CDFG, CEFG, DEFG。比例较高的依次是森林— 灌丛— 耕地— 草地的31.5%, 森林— 耕地— 草地— 水体的20.29%, 森林— 耕地— 草地— 其他非植被的20.84%, 共占72.63%。其他混合形式占比均小于5%。其中包含水体的混合像元占30%。

4.2.4 端元种类数为5

对于5种地物混合的像元, 7种端元类型共有21种混合类型。横轴为不同的组合形式, 纵轴为出现的频率, 绘制直方图如图4d所示。横轴代表的组合形式依次为:ABCDE, ABCDF, ABCDG, ABCEF, ABCEG, ABCFG, ABDEF, ABDEG, ABDFG, ABEFG, ACDEF, ACDEG, ACDFG, ACEFG, ADEFG, BCDEF, BCDEG, BCDFG, BCEFG, BDEFG, CDEFG。比例较高的依次为森林— 灌丛— 水体— 裸地— 其他非植被的32.41%, 耕地— 草地— 水体— 裸地— 其他非植被的19.34%, 灌丛— 耕地— 水体— 裸地— 其他非植被的15.12%和灌丛— 草地— 水体— 裸地— 其他非植被的11.90%, 总共占78.77%, 其他混合形式的比例均小于5%。其中含有水体的混合像元占比91%。

基于以上统计结果, 可以看出在异质性较高的地表场景主要为生态交错区, 即不同植被类型的混合场景。其中, 在每种端元混合数条件下, 频率最高的均是含有森林的混合结构, 而且我国的森林主要分布在山地区域, 这与4.1节中山区的异质性最高的结论一致。此外, 含有水体的混合像元比例也相对较高, 占中国区域所有混合植被像元的12.5%, 且随着混合像元内地物种类的增多, 水体参与的混合类型比例也随之上升。在异质性较高的4种和5种地物类型的混合场景中, 包含水体的混合像元比例超过46%, 因此在高异质性场景中, 水体的影响不可忽略。在异质性地表的遥感反演建模中, 组合特征准确地描述了场景内部的线性混合结构, 为描述异质性场景的宏观特征提供了必要的信息。

4.3 基于边界长度的斑块特征分析

像元内斑块间的拼接特征是地表异质性的典型特征。拼接特征由斑块间的属性差异和斑块的破碎程度决定。属性差异用分类图中的不同类别表征, 破碎程度则可以用像元内的边界长度描述。以2种地物混合的场景为例, 选取不同边界长度对应的地表分类图如图5所示。从图5可以看出, 随着像元内边界长度的增长, 斑块间的交叉逐步增多, 对应的实际地表也越来越破碎。

图5 边界长度对应分类图Fig.5 Corresponding land-cover images of different boundary length

图6a~d分别表示混合像元端元数为2, 3, 4, 5时, 像元内边界长度的频率直方图及边界长度约为20 000 m时对应的地表分类图, 并统计了不同情形下边界长度的均值(m)和标准差(std)。随着端元种类数的改变, 混合像元边界长度的取值区间(0~20 000 m)基本相同, 说明当端元的数目不同时, 混合像元的斑块均可能处于相同的破碎程度。在仅有2个端元混合的场景中, 混合像元出现频率随着边界长度的增长呈迅速下降趋势, 边界长度集中在2 000 m之前; 在3种端元混合场景中, 混合像元边界长度依然集中在2 000 m之前, 之后频率随着边界长度的增加缓慢下降。在4种和5种端元混合场景中, 混合像元边界长度分布更加接近正态分布, 峰值分别为5 000 m和5 500 m。

同时, 观察边界长度的统计值(均值和标准差), 随着像元内参与混合端元数的增加, 边界长度的平均值由2种端元混合的2 292 m逐步增长到5种端元混合时的7 571 m, 边界长度的方差除在2种端元混合时较低外(标准差为3 096), 则呈现了稳定的趋势(3, 4, 5种端元混合时, 标准差均在4 200左右)。这也说明在实际地表场景中, 随着地物种类数的增多, 地表的破碎程度整体上会逐渐增大而波动则相对稳定。

图6 不同端元组合时的边界长度频率分布和边界长度约20 000 m时对应地表分类图Fig.6 Frequency distribution of boundary length at different endmember types and corresponding land cover image when boundary length is 20 000 m

4.4 1 km分类图的异质性特征分析

1 km分类图广泛用于各种应用中。在1 km尺度下, 分类图会存在明显的类别内的异质性特征。本节将分析不同植被类型内部的异质性状况。不同种类的生态系统通常具有不同的异质性特征, 在1 km尺度上不同类型的植被的异质性特征也会有所差异。基于MODIS 地表覆盖图, 我们分别统计了草地/谷物作物、灌木、阔叶作物、稀树草原、常绿阔叶林、落叶阔叶林、常绿针叶林和落叶针叶林边界长度的频率直方图并统计了不同植被类型的边界长度的均值(m)和标准差(std)(图7)。均值代表了该植被类型的边界长度的总体水平, 方差表示该植被类型内边界长度的离散程度。

根据统计结果, 在中低高度植被覆盖区, 草地/谷物作物和灌丛的边界长度总体较低, 均值分别为2 565 m和1 779 m。这2种区域内部的地表异质性较低, 植被尤其是草地生态群系多不与其他高度植被共生。阔叶作物的边界长度均值为3 755 m, 处于较高水平。我国阔叶作物主要分布在东北和华北区域, 农林混交的场景普遍导致地表异质性增大。稀树草原区域是我国平均异质性最高的区域(边界长度均值达7 426 m), 主要分布在贵州、广西西部和云南东南部, 该区域为草地生态系统和东南方向森林生态系统的过渡区域, 树木和灌丛零星分布在草地下垫面上, 这种分布模式会显著增加地表的异质性。在森林区域, 常绿阔叶林、落叶阔叶林、常绿针叶林和落叶针叶林的边界长度均值分别为5 110, 4 486, 4 515和3 184 m。总体上4种森林类型边界长度的分布模式相近, 同时可以看出常绿树木组成的森林内部异质性高于落叶树木形成的森林, 阔叶树木形成的森林的内部异质性高于针叶树木形成的森林。这和不同类型的森林分布的空间位置有关, 落叶、针叶林主要分布中国北部黑龙江的西部, 人口密度低, 植被种类较单一。而常绿、阔叶林多分布在南方山地丘陵区域, 人口密度高, 地形环境复杂导致地表异质性偏高。

图7 1 km尺度下不同生态群系边界长度分布图Fig.7 Frequency distribution of different biomes at 1 km resolution

5 结论

本文提出一种基于高空间分辨率地表覆盖图计算的边界长度并辅以混合像元端元数和端元丰度来定量描述地表的异质性特征的方法。本方法提供了一种定量描述中低空间分辨率像元内部由于斑块拼接特征引起的地表异质性信息的方法, 为异质性地表参数反演方法提供了辅助参数以减少由异质性导致的偏差, 从而提高异质性地表像元的反演精度。组合特征分析提供了混合地表场景的混合形式, 为了解异质场景中植被的高度特征、离散连续特性等提供了辅助信息, 为现有非均质模型的应用提供了参数基础。斑块特征分析提供了地表场景的破碎程度, 可以依据场景的破碎程度来判断采用现有线性混合模型还是异质性地表模型, 为准确决定异质性地表的叶面积指数等地表参数的反演方法及模型提供基础, 以提高参数反演精度。

本文通过对中国区域的地表覆盖类型异质性分析, 发现:

(1) 在中国, 山区具有比平原、盆地、高原等更高的异质性, 山区森林是混合比例最高的一种植被类型。

(2) 中国地表最主要的地物混合方式由不同植被类型混合, 除此以外, 包含水体的混合像元也是一种主要类型, 占中国区域植被覆盖的混合像元的12.5%, 且包含水体端元的混合像元的比例随着地表异质性的增长而增加(由2种地物混合时的4.72%增长到5种地物混合时的91%)。

(3) 边界长度可以有效地描述一定区域内斑块的破碎程度。当边界长度越长时, 地表斑块也就更加破碎。分析MODIS分类图中常见的8种生物群系的边界长度分布特征, 发现稀树草原群系的内部异质性最大, 边界长度均值达到7 426 m, 然后依次为森林群系(平均边界长度为4 323 m)和耕地、草地群系(边界长度均值为3 160 m), 内部异质性最低为灌丛群系, 平均边界长度仅为1 779 m。

本文分析的是地表破碎程度及其特征, 因此是计算所有边界的总和。在用于模拟和反演应用中, 由于只有当太阳入射平面与边界平面垂直时, 边界的交叉辐射最明显, 因此, 需要针对特定传感器过境时的太阳角等分析计算更加敏感的有效边界长度。计算边界长度的地表覆盖分类图的空间分辨率越高, 计算的边界长度越准确。本文使用的30 m分辨率的地表覆盖分类图, 是目前能获取到的空间分辨率最高的全国地表覆盖分类图。随着地表覆盖分类图的空间分辨率提高, 地表覆盖异质性的计算精度也会得到提高。

异质性地表的参数化在异质性地表的建模与叶面积指数反演等领域有着重要的应用空间。如在全球叶面积指数遥感反演中, 目前的算法对地表为均匀假设, 这与实际情况明显不符。需要构建考虑地表异质性的叶面积指数反演算法体系。在本文异质性参数化基础上, 可以实现在多大异质性时, 可以使用简单的线性混合; 在异质性多大时, 需要考虑交叉辐射, 以提高叶面积指数反演精度。需要在本文研究基础上开展更深入的研究。

The authors have declared that no competing interests exist.

参考文献
[1] Hintz M, Lennartz Sassinek S, Liu S, et al. Quantification of land -surface heterogeneity via entropy spectrum method[J]. Journal of Geophysical Research, 2014, 119(14): 8 764-8 777. [本文引用:2]
[2] Haralick R M, Shanmugam K S. Combined spectral and spatial processing of ERTS imagery data[J]. Remote Sensing of Environment, 1974, 3(1): 3-13. [本文引用:1]
[3] De Cola L. Fractal analysis of a classified Land sat scene[J]. Photogrammetric Engineering and Remote Sensing, 1989, 55(5): 601-610. [本文引用:1]
[4] Garrigues S, Allard D, Baret F, et al. Quantifying spatial heterogeneity at the land scape scale using variogram models[J]. Remote Sensing of Environment, 2006, 103(1): 81-96. [本文引用:1]
[5] Chen J M. Spatial scaling of a remotely sensed surface parameter by contexture[J]. Remote Sensing of Environment, 1999, 69(1): 30-42. [本文引用:2]
[6] Wu J, Jelinski D E, Luck M, et al. Multiscale analysis of land scape heterogeneity: Scale variance and pattern metrics[J]. Geographic Information Sciences, 2000, 6(1): 6-19. [本文引用:1]
[7] Gilabert M A, Garcia-Haro F J, Meli J. A mixture modeling approach to estimate vegetation parameters for heterogeneous canopies in remote sensing[J]. Remote Sensing of Environment, 2000, 72(3): 328-345. [本文引用:1]
[8] Liang S. Numerical experiments on the spatial scaling of land surface albedo and leaf area index[J]. Remote Sensing Reviews, 2000, 19(1/4): 225-242. [本文引用:1]
[9] Qin W, Gerstl S A. 3-D scene modeling of semidesert vegetation cover and its radiation regime[J]. Remote Sensing of Environment, 2000, 74(1): 145-162. [本文引用:1]
[10] Widlowski J-L, Pinty B, Lavergne T, et al. Horizontal radiation transport in 3-D forest canopies at multiple spatial resolutions: Simulated impact on canopy absorption[J]. Remote Sensing of Environment, 2006, 103(4): 379-397. [本文引用:2]
[11] Marshak A, Davis A. 3D Radiative Transfer in Cloudy Atmospheres[M]. Heidelberg: Springer Berlin Heidelberg, 2005. [本文引用:2]
[12] Kobayashi H, Suzuki R, Nagai S, et al. Spatial scale and land scape heterogeneity effects on FAPAR in an open-canopy black spruce forest in Interior Alaska[J]. IEEE Geosci Remote Sensing Letters, 2014, 11(2): 564-568. [本文引用:1]
[13] Titov G A. Radiative horizontal transport and absorption in stratocumulus clouds[J]. Journal of the Atmospheric Sciences, 1998, 55(15): 2 549-2 560. [本文引用:1]
[14] Tian Y, Woodcock C E, Wang Y, et al. Multiscale analysis and validation of the MODIS LAI product: I. Uncertainty assessment[J]. Remote Sensing of Environment, 2002, 83(3): 414-430. [本文引用:1]
[15] Zeng Y, Li J, Liu Q, et al. A radiative transfer model for heterogeneous agro-forestry scenarios[J]. IEEE Transactions on Geoscience and Remote Sensing, 2016, 54(8): 4 613-4 628. [本文引用:1]
[16] Yao Yanjuan, Liu Qiang, Liu Qinhuo, et al. LAI inversion uncertainties in heterogeneous surface[J]. Journal of Remote Sensing, 2007, 11(6): 763-770.
[姚延娟, 刘强, 柳钦火, 等. 异质性地表的叶面积指数反演的不确定性分析[J]. 遥感学报, 2007, 11(6): 763-770. ] [本文引用:1]
[17] Tian Y, Wang Y, Zhang Y, et al. Radiative transfer based scaling of LAI retrievals from reflectance data of different resolutions[J]. Remote Sensing of Environment, 2003, 84(1): 143-159. [本文引用:1]
[18] Fang H, Li W, Myneni R. The impact of potential land cover misclassification on MODIS Leaf Area Index (LAI) estimation: A statistical perspective[J]. Remote Sensing, 2013, 5(2): 830-844. [本文引用:1]
[19] Chen Jun, Chen Jin, Liao Anping, et al. Concepts and key techniques for 30 m global land cover mapping[J]. Acta Geodaetica et Cartographica Sinica, 2014, 43(6): 551-557.
[陈军, 陈晋, 廖安平, 等. 全球30 m地表覆盖遥感制图的总体技术[J]. 测绘学报, 2014, 43(6): 551-557. ] [本文引用:1]
[20] Chen J, Chen J, Liao A, et al. Global land cover mapping at 30 m resolution: A POK-based operational approach[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2015, 103: 7-27. [本文引用:1]
[21] Friedl M A, Sulla-Menashe D, Tan B, et al. MODIS collection 5 global land cover: Algorithm refinements and characterization of new datasets[J]. Remote Sensing of Environment, 2010, 114(1): 168-182. [本文引用:1]