前阵子一个师弟抱着电脑来找我,说他在做景观格局分析,数据进了Fragstats之后算出来的指数全是0或者Null,连PLAND这种基础的指标都不对,怀疑是软件装的版本有问题。我看了一眼他准备输入的栅格图层属性,差点没笑出声——坐标系还是GCS_WGS_1984,像元大小是-0.0002度,栅格范围外圈带着一整圈NoData,土地利用类型值从1排到47,里面还有几个0值的背景。这种数据直接送进Fragstats,指数能算对才有鬼。
其实这类问题在景观生态学和生态风险评价方向太常见了。很多人拿到一套土地利用数据,预处理随便做做,然后一头扎进Fragstats把指标一顿跑,最后却不知道怎么解释结果,更不知道这些指数怎么支撑后续的生态风险分析。这篇文章我想把这套完整流程从头到尾捋一遍,从ArcGIS的数据预处理,到Fragstats参数设置和指数计算,再到生态风险指数构建与空间化,每一环都讲清楚“为什么要这么做”,也把这些年替学生、替自己做项目踩过的坑一并交代。如果你正在写这方面的论文,或者在做区域生态规划、国土空间格局优化这类项目,这篇东西应该能帮你省下不少时间。
1. 先把整条链路想清楚:从一张土地利用图到一张生态风险图,中间到底发生了什么
很多新手拿到任务就急着开ArcGIS,但我觉得第一步最应该做的反而不是操作,而是把整条分析链路在脑子里立起来。景观格局指数计算和生态风险分析不是两个独立的模块,它是同一套数据在不同抽象层次上的持续加工。你得先清楚每一环节的输入是什么、输出是什么,才知道自己在哪一步出了问题。
1.1 这套流程到底在解决什么问题
如果你只是想算几个景观指数看看某片区域的破碎化程度,那这个流程其实可以很轻:土地利用栅格进去,选几个指标,Fragstats出结果,Excel里整理一下,完事。
但当你把标题加上“生态风险分析”这几个字,事情就不一样了。生态风险分析在这个语境下,不是指污染场地的人体健康风险,而是指_区域生态风险_——在人类活动和自然因素共同干扰下,区域生态系统结构和功能受到威胁的可能性与程度。景观生态学提供了一个非常实用的视角:土地利用格局本身就是生态过程的载体,格局改变会直接影响物质流、能量流、物种流,因此可以从景观格局指数出发,间接刻画生态风险的空间分布。
换句话说,整条链路可以拆成三层:
- 第一层是“数据层”:土地利用/覆被数据经过ArcGIS预处理,变成标准化、可计算的栅格数据。
- 第二层是“指数层”:Fragstats基于栅格数据计算景观格局指数,量化斑块、类型和景观三个尺度的格局特征。
- 第三层是“风险模型层”:把景观指数作为自变量,通过干扰度、脆弱度或者多指标综合模型,构建生态风险指数,并落到空间单元上呈现。
这三层相互依赖,前一层没做干净,后一层全是垃圾。很多人恰恰是在第一层偷了懒,导致后面几百块的模型结果根本没法解释。
1.2 数据基础怎么选:30米分辨率几乎是这条路的默认起点
我见过不少用MCD12Q1或者GlobeLand30做这类分析的文章。坦白讲,对于中小尺度的区域生态风险评价,30米分辨率基本是底线。15米甚至更高分辨率的土地覆被数据当然更好,但受限于数据获取难度和软件算力,30米是性价比最高的选择。
需要特别提醒的是:分辨率选择必须写进方法部分,并且保持一致。你前后两个时期做对比分析时,如果一期用的是30米数据,另一期用的是250米数据,计算出来的几乎所有面积相关、密度相关、形状相关的指数都不可比。这不是数据精度问题,是尺度本质不同。Fragstats不关心你数据真假,它只会按照你给它的栅格如实计算,但结果的生态学意义全掌握在你自己手里。
1.3 先把“研究区范围”和“类别体系”定死,后面能少改很多图
再啰嗦一个细节:在开始所有预处理之前,先把研究区边界确定好,并且从头到尾使用同一个边界矢量文件。有的同学中途觉得研究区边界不妥,又换了一个,结果所有栅格裁剪范围变了,像元数量变了,面积和周长全部跟着变,前面跑过的Fragstats结果等于白做。
类别体系同样要提前定。你的土地利用数据可能来自某个分类产品,自带一级类、二级类,比如6大类或10小类。但景观格局分析未必需要你把所有类别都保留——比如耕地里细分的水田和旱地,从生态风险脆弱度角度看它们差异未必大,从格局指数角度它们却会显著影响斑块密度和多样性指数。做合并会损失信息,不合并会稀释格局信号,这个权衡没有绝对标准,但必须在分析前定下来,并且写清楚合并依据。
2. ArcGIS侧的数据预处理:投影、重分类、栅格尺度,任何一个没弄好都是给Fragstats埋雷
这章是整条流程的重中之重。我自己判过不少学生的中期材料,发现绝大多数Fragstats“算不出”或者“算出来明显不对”的问题,根源根本不在Fragstats,而在ArcGIS预处理阶段。咱们把该做的步骤一个个过一遍。
2.1 第一步永远是投影:地理坐标系直接拿来算面积会闹笑话
我那个师弟的数据麻烦就在这儿。他的土地利用栅格坐标系是GCS_WGS_1984,也就是经纬度地理坐标系。Fragstats在计算斑块面积、周长、密度时,默认空间单位是米——如果你喂给它的是以度为单位的栅格,它计算出的每一块斑块面积都是以“平方度”为单位的,SPSS里看着数字挺大,实际上完全不能解释。
正确做法是投影到适合你研究区位置的等积投影。对国内大部分区域,可以用阿尔伯斯等积圆锥投影(Albers Equal Area Conic),中央经线按研究区中心设置,双标准纬线按区域跨度设置;如果是小范围研究区,用UTM对应分带也很省事。等积投影能最大程度保证面积计算准确性,这对面积类指数是基础性的。
在ArcGIS里操作时,用“Project Raster”工具而不是“Define Projection”。Project Raster是真正的重投影,会重新计算坐标并重采样像元;Define Projection只是给数据打一个坐标标记,如果你把用于纠正数据的坐标系声明错了,后面全乱套。
2.2 重采样与像元大小:让面积指数不再“虚胖”
投影转换过程中,ArcGIS会要求你选择重采样方法。对类别型(土地利用)数据,我强烈建议选最邻近法(Nearest),因为双线性或三次卷积会在类别边界插值出不存在的“新类别”,破坏类型取值的原始语义。
像元大小也要在这个步骤一并定下来。如果你的多个时期数据原始分辨率不同,一定要统一。30米就统一成30米,建议直接把Cell Size设成30,不用在意研究区边界是否完全对齐。Fragstats对栅格尺寸和像元对齐的要求没有你想象的那么严格,但计算面积、周长时要设定换算因子(见下一节),心里得有数。
2.3 类别合并与重分类:47个类别进Fragstats,你什么都解释不了
土地利用分类数据往往有几十个类别。Fragstats处理能力和类别数没有必然上限,但生态学解释层面,你把47个类别全部塞进去,结果表密密麻麻,类型水平上根本没法逐类讨论,景观水平的多样性指数也会被过度稀释。
我通常的做法是合并到5-8类,比如:耕地、林地、草地、水域、建设用地、未利用地。用“Reclassify”工具,把原始类别值映射到新值,比如耕地=1,林地=2,草地=3,水域=4,建设用地=5,未利用地=6。这一步有一个极其容易踩的坑:重分类后的栅格一定要看一下属性表,确认没有残留的0值、NoData或原始类别漏映射。Fragstats中的背景值(Background)会被排除在计算范围之外,但如果你背景值设了0,而栅格里面真实存在大量0像元,那么这些0值像元就会被当成背景忽略掉,周围的斑块就会被不合理地切开,直接影响斑块数目和聚合度指数。
所以在新栅格中,任何未赋值的区域都应该是NoData,而不是0。这个区别太重要了,后面我还会专门讲。
2.4 矢量转栅格与裁剪:用研究区掩膜把范围卡死
如果你手里的土地利用数据是矢量格式,则必须先转栅格。用“Polygon to Raster”时,唯一值得留意的是Cellsize设置,一定手动填30或你统一好的分辨率,千万别默认。另外Value Field要选到类别编码字段,而不是面积或编号字段。
裁剪操作建议统一用研究区矢量作为掩膜,用“Extract by Mask”处理一次即可。如果直接在ArcGIS里Display->Data->Export Data的方式导出,容易把环绕整图的“背景区域”一并导出,这等于把Fragstats的背景值范围做大了,NoData处理层面的问题就会卷土重来。
2.5 Fragstats前的最后一步:把栅格导出成干净的无背景格式
好,现在你已经有了一个投影正确、分辨率统一、类别值合理、范围贴着研究区边界的栅格。最后一步是导出。ArcGIS的GRID格式Fragstats也能读,但我个人更习惯导出GeoTIFF,压缩选项选LZW,NoData值用-9999显式标注。这一步不要直接存成.img或者.gdb里的栅格,虽然软件都能读,但后期你可能会换电脑、换软件版本,跨平台读GeoTIFF是最不易出问题的。
导出的同时,顺手在ArcGIS里给栅格做一个属性表检查:统计每个类别的像元数量,乘以像元面积,看看是否符合研究区各类型的面积量级。这一步相当于“数据体检”,能提前发现重分类串号、掩膜偏移之类的问题。
3. Fragstats指数计算:参数是窗体上那几个下拉框,但背后是景观格局的三种尺度
当你的栅格干干净净进入Fragstats之后,计算本身反而很快。但窗口上那些参数不是随手选的,每改一个下拉框,最终算出来的数学含义就变一个维度。
3.1 创建Image File:背景值设错的后果是灾难性的
Fragstats读的不是TIFF本身,而是要先通过它的“New Image”按向导创建一个.image文件描述。在这个界面里有两处必须手工确认,一是Background Value,二是指定输入的栅格格式和图层路径。
背景值建议设成NoData对应的值,比如-9999。为什么不能说设成0就万事大吉?因为如果你的类别值里设置了0类(比如水体类编码为0),背景值如果也是0,整片水体就被当成背景剔除了。更麻烦的是,有些栅格导出时NoData已经被烧成了0(这是ArcGIS导出设置不当的常见结果),你设背景值为0,恰好把NoData和真实类别全部抹平。所以预处理阶段才反复强调:NoData就是NoData,类别就是类别,别让它们在数值上混淆。
3.2 邻域规则:8邻域是普遍选择,别稀里糊涂用4邻域
Fragstats里有一个Contiguity Rule选项,默认是8邻域(上下左右加四个对角),也支持4邻域。这两个选项在斑块识别时的结果有系统性差异:4邻域会把对角相邻的区域断成两个斑块,8邻域则把它们归为一块。
生态学上,八邻域通常更符合物种移动和生态过程的连通性理解,因为斑块之间的生态流不严格沿加州经纬方向。但如果你研究的是某些特定扩散方式非常受限的物种,4邻域也有其合理性。多数文献用的还是8邻域,我建议没有特殊理由都选8邻域,并在论文方法里明确标注,因为这个细节会影响一切和斑块计数、连接度相关的指数。
3.3 Patch、Class、Landscape三个层级:你关注的核心科学问题决定了算哪一层
Fragstats输出分为三个层级,这是整个软件最核心的概念框架:
- 斑块水平(Patch level):每个独立斑块一个值,对应“最小空间单元的特征”。
- 类型水平(Class level):同一类别所有斑块汇总后的统计值,对应“某一土地覆被类型的结构与破碎度”。
- 景观水平(Landscape level):整个研究区作为一个整体,对应“区域景观整体异质性与多样性”。
生态风险分析真正常用的是后两者。比如你想回答“建设用地扩张在研究区造成的破碎化效应是什么”,看建设用地的类型水平指数比较有意义;你想回答“整个研究区的景观多样性变化,空间上哪里风险高”,那景观水平指数配合空间格网是关键。
3.4 指数组合别贪多:常规组合已经能讲清楚90%的问题
Fragstats能算的指数上百个,但实际分析中,完全不需要把它们全勾上。我给自己定过一套常规组合,适用于大多数土地利用/生态风险场景,你可以直接参考:
| 尺度 | 指数 | 缩写 | 说明 |
|---|---|---|---|
| 斑块 | 斑块面积 | AREA | 基础统计量 |
| 类型 | 斑块类型面积百分比 | PLAND | 类型面积占比 |
| 类型 | 斑块密度 | PD | 破碎化程度 |
| 类型 | 边缘密度 | ED | 边界效应 |
| 类型 | 平均形状指数 | SHAPE_MN | 形状复杂度 |
| 景观 | 蔓延度 | CONTAG | 聚集与蔓延趋势 |
| 景观 | 散布与并列指数 | IJI | 类型间相邻关系 |
| 景观 | 香农多样性指数 | SHDI | 多样性水平 |
| 景观 | 香农均匀度指数 | SHEI | 均匀度水平 |
这些指标基本覆盖了景观生态风险分析最常见的几个维度:组分、密度、形状、聚集和多样性格局。你后面构建生态风险指数时也主要从这些维度里挑变量。真的不需要在一篇论文里堆15个指数,审稿人看到一堆弱相关指数的简单罗列并不会觉得你厉害,反而会让方法部分显得缺乏筛选逻辑。
3.5 批量运行:同一个模板对多个时期/多期栅格批量算
如果你做的是两个或三个时期的对比,不必每个栅格重新走一遍向导。在Fragstats的Input Layer窗口里可以同时添加多个图像文件,同样设置下一次性批量计算。输出时它会为每个图层分别生成结果,文件名自动带后缀,不会相互覆盖。
还能在Output选项里勾选同层水平汇总表(Class Statistic和Landscape Statistic),结果直接输出成标准文本或CSV,导出后在Excel里整理即可。注意Fragstats输出的面积默认单位可以在Preferences设置里调成公顷或平方千米,这一点要和你的格子尺度保持统一,后面做风险小区合并时才会顺手。
4. 生态风险指数构建与空间化:网格化、权重、插值,一步步把格局“翻译”成风险
指数算完,只是拿到了“格局描述”,距离“风险评价”还得再铺一层桥梁。这一章聊聊后面这半条路,也是很多论文里方法差异最大的部分。
4.1 研究区网格化:用Fishnet把连续景观切成可统计的空间单元
景观水平指数本身是全区一个值,既没办法落到空间上看分布,也没法和具体的行政单元或生态单元做叠加分析。所以一个常见做法是把研究区分成若干格网,每个格网作为一个“风险小区”,格网内部的景观指数聚合后,作为该格网的风险表征。
在ArcGIS里用“Create Fishnet”工具操作。这个格网尺寸其实没有统一规定,一般视研究区面积而定。我给中小尺度区域做时,常用2km×2km的格网;如果是省份级别的宏观评价,可以用5km甚至10km。格子太大,风险图太粗糙,看不出内部差异;格子太小,很多格子里面只有零星几块地,统计出来的指数不稳定,噪声非常大。
一个经验法则是:格子边长大约是主要景观斑块平均面积的10到20倍,或者让研究区内格网数量保持在200到500个之间,这个范围内的统计比较稳定,出图效果也比较好。
4.2 把Fragstats结果按格网聚合:连接字段是唯一标准,别手动画框选
每个格网内部的景观指数怎么取得?一种比较费劲但有些人真这么干的错误做法是,把Fragstats的斑块水平输出导到ArcGIS里按位置连接——这没有必要。实际操作可以用“Tabulate Area”或者把斑块水平输出表和格网唯一标识做空间连接(Spatial Join),按格网ID聚合统计出PLAND、PD等值。
如果你更习惯用景观水平指数作为风险变量,更简单的做法是:用“Extract by Mask”把每个格网区域单独裁出来,依次跑Fragstats,再拼回一张表。但这样效率太低了,强烈建议把栅格叠加到格网上,输出斑块级结果后用Pandas、Excel或者ArcGIS表格连接汇总到格网,效率和结果一致性都高得多。
4.3 核心模型一:经典“干扰度+脆弱度”加权法
这是生态风险评价里应用最广的一类方法,逻辑非常直观:某个格网的生态风险 = 景观干扰指数 × 景观脆弱度指数。
景观干扰指数(有些文献叫景观损失度指数)通常是三个分量的加权合成:
- 景观破碎度:表征格局被切割的碎化程度,可以用斑块密度PD或分割度表达;
- 景观分离度:表征同类斑块在空间上被隔开多远;
- 景观优势度:表征某类景观对整体格局的控制程度;
- 也可引入分维数倒数来表征形状复杂度,用于刻画人为活动干扰引起的形状规整化倾向。
具体公式形式不同文献略有差别,常见形式为:
Di = a·Ci + b·Si + c·Fi其中Ci为破碎度,Si为分离度,Fi为分维倒数(或优势度),a、b、c为权重,常取0.5、0.3、0.2,代表研究者的主观判断。这一组权重在很多经典文献里都出现过,但严格讲它仍是先验权重,你需要在自己的数据上下文里给出取舍理由。
脆弱度指数则与土地覆被类型本身的生态学属性有关:水域、滩涂这种对干扰敏感的生境脆弱度高,建设用地和裸地脆弱度低。常用做法是先给各类别赋一个脆弱度评分,比如未利用地1、林地2、草地3、耕地4、建设用地5、水域6(顺序可按研究需要调整),然后归一化到0到1之间。最后生态风险指数ERI就是干扰度和脆弱度的乘积对面积加权汇总。
这种方法的优点是机理清晰、计算简单、结果容易被审稿人接受;缺点是权重的主观性较强,需要结合区域特征做论证。
4.4 核心模型二:基于多指数的综合风险指数
另一种近年很流行的思路是:不预设权重,而是先选一组景观格局指数(如PD、ED、SHDI、CONTAG等),用主成分分析(PCA)提取主成分,再按主成分方差贡献率加权合成生态风险指数。也可以进一步用熵权法确定权重,避免太多主观成分。
熵权法的大体计算套路是这样的:
- 对每个格网的每个指数做标准化(正向指数用最大值归一化,负向指数用最小值归一化或取倒数);
- 计算每个指数的信息熵Ej = -(1/ln n)·Σ pij·ln pij,其中pij是格网i在指数j下的比重;
- 权重Wj = (1-Ej) / Σ(1-Ej)。信息熵越小,说明该指数在不同格网间差异越大,对区分风险空间的贡献越高,权重自然就大。
这样算下来的风险指数同样是一个格网一个值。它的优点是少了很多人为拍脑袋的权重,可复现性更强;缺点是你需要在方法部分交代清楚对每个指数做的是正向还是负向处理。比如CONTAG值越大通常代表景观越连续聚集,但你如果拿它作为“风险正向指标”处理,就必须说明为什么聚集本身在该研究区代表着高风险。
4.5 空间化表达:插值、分级与出图
格网化的风险指数本质上是一张离散的点或面图层。如果格子足够小,把中心点提取出来,用反距离权重插值或普通克里金插值,就能得到连续的风险面。一般论文里做这一步是让风险图更好看、更好解释,不要过度迷信插值的统计学意义。在做插值前,把格网中心点导出,ArcToolbox里IDW或Kriging跑一下,搜索半径设到能够覆盖相邻4-6个格网的尺度即可。
风险分级常用自然间断点法(Jenks),因为它能最大化类间差异,是风险分区制图最稳妥的默认选择。分5个等级,低、较低、中、较高、高,每个等级的面积占比在结果表里统计出来,写进结论部分当核心量化依据。
这个时候你做出来的生态风险图才真正可以拿去跟土地利用图叠加,分析“高风险区主要分布在哪些乡镇”“高风险的驱动因素是什么”——到了这一步,前面所有准备工作的价值就全都体现出来了。
5. 我替你们踩过的坑:Null结果、指数失真、误读联动的排查过程
最后这张就当排错手册来写。下面几个问题,我不敢说每个做景观分析的人都遇到,但概率相当高,而且一旦出现,往往排查要耗费大半天。我把排查思路也一并说透。
5.1 Fragstats结果全是Null或大面积0:先看背景值,再看栅格数值范围
这是最经典的一个坑。症状是计算结果表里,某个类型的PD和PLAND全是0,或者整个景观水平输出是Null。
排查顺序不要乱:
- 第一步,在ArcGIS里打开栅格属性表,看极值。如果最大值是255这类默认背景值,说明导出时NoData被自动写成了255或0,而你在Fragstats里设置的背景值对不上;
- 第二步,在Fragstats的Image File设置里,重新指定背景值,让它和栅格属性表里的实际NoData数值一致;
- 第三步,重跑。八成问题就解决了。
如果属性表里压根没有NoData,但边界一圈确实是空白,那就说明预处理阶段栅格本身带了背景区域(比如扩展到了矩形外包范围),直接裁到研究区再重跑。这一点我在第2章反复强调过,做的时候不以为然,出错的时候是真急人。
5.2 面积类指数大得离谱:十有八九是投影坐标系没换干净
如果你算出来的平均斑块面积(AREA_MN)是几百万的值,先别想着是不是软件单位设置错。去ArcGIS里右键图层属性看来源,确认坐标系是不是仍然是经纬度,或者投影虽然转了,但栅格导出时选择了“Source”而不是“Projected”选项。
很多看似蹊跷的结果,最后查来查去,都是坐标系这类最基础的问题。所以我建议每做完一版预处理,花十秒钟看一次属性,养成习惯。
5.3 边界效应:研究区边缘的斑块形状被切得不自然,怎么办
研究区当然不可能都是规则矩形,实际边界会把一些斑块从中间截断,导致边缘区域的形状指数、面积统计被扭曲,这就是边界效应。Fragstats本身不提供缓冲消除机制,但你可以做一个缓冲区实验:把研究区向内缓冲一个像元(比如30米),重新算一遍,对比结果变化幅度。如果变化很大,说明你的边界切割确实严重,就应该在方法部分明确声明,并将内缓冲结果作为稳健性检验。
这个方法写进论文里反而是加分项,能显示你对尺度效应的理解和严谨性。
5.4 安装与License相关的“启动没反应”问题:一句话带过的解决方案
有同学问ArcGIS Desktop 10.8在Win11上双击没反应,或者License Server启动不起来。这个问题的根源多半是服务没有正常启动或者兼容模式设置问题。做一个简单的处理方向:以管理员身份运行License Server Administrator,服务列表里找到ArcGIS License Service,手动启动,再以管理员身份运行ArcMap;如果还是不行,检查安装目录下是否有防火墙拦截,或者用兼容模式运行。不要急着重装软件,这种问题九成是服务层面的,不是安装包的锅。
但我也要说句实在话:ArcMap本身已经停止更新,ArcGIS Pro对新电脑和新系统兼容性好得多,新项目建议直接上Pro。Fragstats和Pro之间的数据交互也很成熟,没有任何障碍。
5.5 指数解读的禁忌:SHDI高不等于风险高,别把生态学意义搞反
最后一个坑是方法层面的,但比上面任何一个都重要。
景观多样性指数SHDI或SHEI高,只代表景观组成更丰富、分布更均匀,它本身并不同步等价于“生态风险高”。有些生态风险评价论文里直接拿SHDI当风险正向指标用,这在概念上是有问题的——一个天然草甸、农田、水域镶嵌的异质性景观,SHDI可能很高,但生态风险未必高;反过来,极度破碎的建设用地基质的城市景观,SHDI可能中等,局部风险却很高。
所以在构建综合风险指数时,每引入一个指数,先问自己一句:它在我的研究区语境下,生态风险升高到底是这个指数升高还是降低?想不清楚就不要硬塞,去翻前人针对同类区域的做法,或者回到景观生态学的基本概念里找依据。
最后说句实在话
整套流程走下来,真正花时间的从来不是Fragstats跑指数的那几秒,而是ArcGIS里一遍遍的投影检查、重分类核对、背景值确认,以及后面构建风险指数时对每个指标生态学含义的反复推敲。我自己的经验是,前期的数据预处理占到这个工作量的六成以上,但这也是最不容出错、最值得花时间的部分。很多刚接触的人恨不得一步到位算出个风险图,但越急越容易在第一步形态上翻车。把前面的功夫下足,Fragstats只是十几分钟的事,后面的风险指数解释才真正开始体现功底。
如果你正准备用这套流程做自己的项目,我的建议是:第一次走全流程时,用一个自己熟悉的小研究区试跑,把所有参数记录下来,做成一个自己的模板,后面再换大区域、多时期的时候,效率会快非常多。顺手一个小技巧:给每个栅格类别编码用连续正整数,1、2、3这样排下去,背景单独用一个肯定不冲突的负值或NoData表示,这个习惯会让你在Fragstats和ArcGIS之间的来回调试省下大量时间。