环评项目里最耗时间的往往不是预测模型本身而是生态专题那一摞图件。我印象很深的一次某公路项目生态评价专章要求出土地利用现状图、植被覆盖度图、生态系统类型图、景观格局指数表外加一张几十年的趋势对比图。当时还是老思路靠ArcGIS一亩一亩人工勾绘光两个图斑就折腾了两天图例顺序和颜色体系还得反复调。后来我把整个工作流换成了ENVIRFragstats三件套ENVI做遥感解译Fragstats算景观格局指数R负责数据汇总和批量出图整个制图环节从一周压缩到一天半图表质量反而更稳定。这篇就把这套流程完整捋一遍适合正在编制环评报告书、尤其是生态影响评价专章的同行参考。想直接跑通的人按下面的章节顺序操作就行每步我都写了参数和避坑点。1. 环评图件到底有哪些先拆清单再定工具分工1.1 报告书里躲不开的图件从基础信息图到生态专题图环评报告书里的图件数量新手往往低估了。按工作内容拆开看其实就三类。第一类是基础信息图。包括项目地理位置图、评价范围图、环境监测点位图、敏感目标分布图、总平面布置图。这些图的核心是矢量化底图加注记ArcGIS或QGIS就能干跟遥感解译关系不大图面整饰也就那样不再展开。第二类是现状评价图这是图件工作量的大头。土地利用现状图、植被类型图、植被覆盖度图、生态系统类型图、土壤侵蚀等级图基本都依赖遥感影像解译得到基础数据。这类图恰恰是ENVI的主场从影像预处理到分类出图一条链路走下来数据来源清清楚楚。第三类是预测和分析类图。大气浓度等值线图、声环境预测图、地下水和土壤影响图以及生态影响评价里的景观格局对比图、生态质量变化趋势图。这些图的核心是统计计算和可视化R语言做这类生产工具最顺手尤其当你要对十几年多期数据做批量对比时Excel根本撑不住。三类图放一起看真正消耗人力的就是现状评价图和趋势对比图而这部分正好能交给ENVIRFragstats这套组合。很多人一上来就埋头在ArcGIS里画画完了数据来源说不清、精度验证没做评审会上被问一句就哑火。我的建议是先列一张图件清单按数据生产链路去分工而不是按图件类型去逐个手画。1.2 为什么是ENVIRFragstats三条流水线如何分工先明确一个逻辑这套三件套不是三个软件各干各的而是一条生产流水线。ENVI负责第一段把原始影像变成分类栅格。它的辐射定标、大气校正、监督分类、精度评价都是菜单化操作比命令行写GDAL门槛低很多。环评工程师的时间不应该耗在写底层代码上ENVI在学习成本和处理效率之间平衡得最好。如果只是做常规生态本底调查不需要复杂编程ENVI足够。Fragstats负责第二段把分类栅格变成景观格局指数。这个软件几乎是景观生态学领域的行业标准指数全、算法公开评审专家认可度高。环评导则里要求描述生态系统结构、功能和破碎化程度时用Fragstats算出来的NP、PD、SHDI这些指标比你自己在Excel里拍脑袋算的要有说服力得多。R语言负责第三段把所有数据汇总成报告能用的图表。它的统计能力和可视化能力正好补上前面两个软件的短板把多年土地利用面积堆在一起做对比图、把Fragstats输出的一堆表格清洗成一张主表、批量输出300dpi的PNG这些都是R的强项。下面三章就按这个流水线顺序展开每一步的实操细节我都会写到。2. ENVI出分类图从原始影像到一份能进报告的土地利用现状图2.1 影像数据源选择Landsat、哨兵与高分的取舍环评制图的第一步不是打开ENVI而是选影像。选错数据源后面做得再漂亮都是白费。判定标准主要看评价范围尺度和精度要求。大区域的线性工程比如公路、铁路、管道评价范围动辄几十平方公里甚至上百用30m分辨率的Landsat或者10m分辨率的哨兵2号就够了两者都可以免费下载。Landsat的优势是历史数据长可以追溯到上世纪八十年代做生态趋势对比非常好用哨兵2号分辨率更高更新频率快适合近几年的本底调查。小范围的厂区、矿山、风电场评价范围只有几平方公里这时必须上高分影像。GF-1、GF-2或者商业的高分影像能做到亚米级或米级分辨率能看清小尺度的地类边界。这类影像成本高一些但生态底图的精度是评审专家盯着的点不能省。还有一类特殊情况矿区的生态影响评价。如果项目涉及尾矿库、采空区需要提取植被胁迫或异常信息可以考虑使用GF-5高光谱数据做蚀变信息提取这是ENVI的光谱分析强项但常规环评图件用不上属于进阶技能这里先带过。选影像还要注意两个细节。一是时相尽量选6到9月的生长季植被光谱特征最明显分类精度最高。二是云量控制在10%以下有云的区域直接放弃重选。跨景拼接前先确认各景影像是否已经做过几何校正时相对不齐的影像拼出来容易出现明显接边痕迹。2.2 预处理四连辐射定标、大气校正、几何校正与裁剪拿到原始影像别急着分类先把预处理走完。第一是辐射定标。原始影像存的DN值是个无量纲数字需要变成辐射亮度或反射率才能做后续计算。ENVI里的路径是工具箱Radiometric Calibration选择多光谱文件后直接输出定标结果这个操作没有太多参数要调但一定要做不做后面的大气校正就无从谈起。第二是大气校正。大气中的水汽和气溶胶会吸收和散射电磁波导致地物反射率失真尤其是植被指数计算不做大气校正NDVI值普遍偏低。ENVI里常用FLASH模块需要填的传感器类型、影像中心经纬度、飞行高度、成像时间以及大气模型和气溶胶模型。大气模型按成像季节和纬度选中纬度地区夏季选Mid-Latitude Summer冬季选Mid-Latitude Winter。如果时间紧不想细调参数也可以用QUAC快速大气校正精度会差一些但环评本科底图够用。第三是几何校正。Landsat L2级产品自带精确几何定位不需要额外处理。高分影像如果发现地物边界和矢量底图有偏移就要手动选控制点做一次几何校正控制点尽量选在道路交叉口、水库坝体这类特征明显的点上均匀分布RMSE控制在0.5个像元以内。最后是裁剪。用评价范围的shp边界在ENVI里执行Subset Data from ROIs/Shapefile输出自带坐标信息的GeoTIFF。裁剪范围建议比评价范围外扩1到2公里给后续制图留出修边余地否则评价边界紧贴影像边缘图面上很难看。2.3 监督分类与精度验证ROI怎么选、SVM怎么用预处理完成后进入核心环节土地利用分类。我推荐监督分类里用SVM支持向量机RBF核小样本情况下比最大似然分类稳定对椒盐噪声的抑制也好一些。分类体系要提前定好结合环评导则和项目特点常见的类别是耕地、林地、草地、建设用地、水域、裸地这六类有特殊项目再增加园地、沼泽、盐田等。类别编码必须从一开始就固定耕地1、林地2、草地3、建设用地4、水域5、裸地6后期多期对比全靠编码对齐。ROI选取是整个分类流程里最考验经验的一步。每个类别至少选取30到50个样本样本要覆盖同一地类的不同光谱表现。比如耕地里有水田和旱地光谱差异很大不能只在一个区域点样本就完事。ROI尽量画成闭合多边形不要只点一个像元点选像元容易选到混合像元分类时容易带偏。训练样本和验证样本要分开。ENVI里可以通过样本分割工具把ROI按比例拆成训练集和验证集我一般按3比1拆分。分类完成后用验证样本计算混淆矩阵重点看总体精度和Kappa系数。环评报告里总体精度85%以上、Kappa大于0.7基本就能应对评审。如果达不到不要急着调分类器先回去补样本样本质量上来了精度自然就上来了。分类完成后还要做后处理。ENVI里的Majority/Minority分析3x3窗口把分类图里孤立的小斑块归到周围主要类别里能有效去掉椒盐效应。跑完之后建议再人工快速扫视一遍重点检查大水体、城镇边界有没有明显错分。2.4 多期分类的统一性问题分类体系、时相与分辨率对齐做生态趋势分析的同行一定会遇到多期影像分类的场景这里有几个坑必须提前踩平。最大的坑是分类体系不一致。不同年份的影像由不同的人分类别定义稍有差异后面对比就全是噪声。我的一般做法是在一张Excel里把每个年份的类别编码表固定下来分类时严格按统一编码执行任何人接手都先看这张表。第二个坑是分辨率不一致。早期Landsat影像只有30m后来用了10m的哨兵如果直接对比面积地类边界会因像元大小不同产生系统偏差。处理办法是统一重采样到相同分辨率我通常以较低分辨率30m为基准高分辨率影像重采样下来再参与面积统计。第三个坑是时相不对齐。不同年份影像如果成像月份差很多植被光谱差异会被误判成土地利用变化。做多年对比时尽量选同一季节的影像如果实在找不到同期影像报告里要明确说明时相差异带来的不确定性别等评审专家来问。每期分类完成后记得导出三个东西分类GeoTIFF、分类结果图以及分类前的原始影像路径和成像日期。后面写报告的数据来源说明时这些东西都是要用的。3. Fragstats算指数参数设错整张景观格局表就废了3.1 输入文件准备把ENVI分类结果喂给Fragstats之前必须做的事ENVI分类结果并不能直接丢给Fragstats中间还差一步格式整理。Fragstats 4.2支持GeoTIFF格式输入但要求像元值是连续的整数类别编码。ENVI分类输出的栅格像元值一般已经是从1开始的整数但可能存在一些非分类像元比如背景区域被赋值成0或者255。最关键的一步就是把这些背景值单独处理掉我会建议统一把背景设为0并且保证类别编码从1开始连续编号类别数量不要超过几十个。这一步在ENVI里可以用Band Math或重分类工具做但我个人更习惯在R里用terra包处理因为可以写成脚本批量跑多期数据一次搞定。比如用subst函数把特定的旧值替换成新的类别编码再统一写回GeoTIFF。处理完的栅格在Fragstats里加载时还要顺手看一眼属性表确认没有出现类别值跳号的情况。3.2 三个关键参数背景值、邻域规则和边缘深度Fragstats主界面的参数看着简单但三个参数一旦设错算出来的指数全盘作废这个坑我踩过太多次了。第一个是背景值设置。在Class Properties里明确指定哪个值代表背景。如果不设背景像元会被当成一个类别参与面积和密度计算NP、PD会全部偏大LPI也可能被背景斑块主导。这是我见过同行最常犯的错误很多人跑完结果就直接粘到报告里完全没意识到指数已经被背景污染了。第二个是邻域规则。4邻域只考虑上下左右四个方向8邻域把对角方向也算进来。这个选择直接影响聚合度、连接度、散布与并列指数等一大类指标的结果。环评里我一般用8邻域因为像元间的生态过程往往是八方向连通的但报告的方法部分必须把这个选择写清楚评审专家会看。第三个是边缘深度。默认设置是1个像元够用。但如果要算核心区面积、核心区密度这些核心区指标就必须根据真实生态影响距离来设置。比如某生态敏感鸟类的影响距离是50m30m影像下50m换算成边缘深度就是2个像元那就要把边缘深度设为2。设置不当核心区面积会严重高估。参数设完后输出层面建议同时选class和landscape两个层级格式选CSV这样每个指标按类别算了一份同时也有一份全景观层面的汇总。导出的.land、.class文件直接用Excel或R打开即可。3.3 报告里真正用得到的景观指数一张表讲清生态含义Fragstats能算几十个指数但环评报告里真正用得到的其实就七八个。选指数的逻辑不是越多越好而是每个指数都要能对应到生态学含义并且能说出生态过程。下面这张表是我在项目里经常用到的核心指数清单指数缩写全称生态学含义环评里的用法NP斑块数量景观破碎化程度数值上升说明生境被切割PD斑块密度单位面积斑块数与NP配合消除面积影响LPI最大斑块指数优势斑块占比反映基质类型和优势度LSI景观形状指数斑块形状复杂度形状越复杂边缘效应越强CONTAG蔓延度景观连通程度下降通常说明连通性变差AI集聚度同类斑块聚集程度说明类别空间分布偏好SHDI香农多样性景观异质性上升说明类型变多、结构复杂SHEI香农均匀度各类型分布均匀程度辅助SHDI解释我在实际报告中写生态影响评价时一般这样组织先把项目区各类用地的面积占比列出来再用NP、PD说明破碎化用CONTAG、AI说明连通性用SHDI、SHEI说明景观多样性变化。每写一个指数后面都要跟一句生态过程解释。比如PD上升说明生境斑块被进一步切割这对于移动能力弱的物种潜在影响更大不能光罗列数字。3.4 多期景观格局对比的实操顺序多期对比的实操顺序其实很简单把每期分类图按3.1节的方法处理成标准输入然后在完全相同的参数设置下分别计算指数。注意参数一定要一致背景值、邻域规则、边缘深度都设成相同值否则不同年份的指数差异包含了方法差异对比就没有意义了。算完多期指数后我一般会建一张主表行是年份列是指标再辅助上一张R出的折线图或柱状图把NP、SHDI这些关键指标的变化趋势直观呈现出来。Fragstats本身不提供可视化这正好接上下一章R语言的内容。评审会上被问得最多的就是这些指数变化的生态学意义是什么所以我建议在出图前先把生态过程解释写好再让图表配合文字而不是先出图再想解释。这也是整个Fragstats环节最核心的写作逻辑。4. R把三件套串成流水线数据整合、ggplot2批量出图与内存坑4.1 R环境与包管理从安装到没有程辑包报错R环境的搭建没有捷径但也就花二十分钟。先在CRAN官网下载安装R然后装一个RStudio日常操作都在RStudio里完成。装包是新人第一道坎。install.packages(包名)是最常用的命令如果遇到不存在叫‘xxx’这个名字的程辑包的报错先别怀疑人生按三步排查第一步查包名拼写R里下划线、大小写都区分第二步查镜像设置Tools - Global Options - Packages里把CRAN镜像换到国内镜像不然下载经常断第三步确认这个包是不是在GitHub上发布的如果是先安装devtools再用install_github(作者/包名)安装。RStudio里我用得最勤的快捷键是CtrlEnter运行当前行和CtrlShiftR插入代码段这两个键记住了写长脚本会顺很多。环评生态分析常用的包有terra和sf处理栅格矢量dplyr和tidyr做数据清洗ggplot2和patchwork出图readxl和writexl读写Excel。少数生态学专用包在Bioconductor仓库里发布需要先安装BiocManager再用BiocManager::install(包名)。这一整套装完后面几乎不会再碰包安装问题。4.2 数据整合栅格统计、Fragstats指数合并成一张总表数据整合是R在这套流水线里最核心的活儿。我用一个具体场景演示三个年份的土地利用分类栅格需要统计每个年份各类别面积再把Fragstats算出的景观指数一并合并。先看栅格统计。用terra包批量读取多期分类结果统计每个类别像元数乘以像元面积换算成平方公里代码很简洁library(terra) library(dplyr) # 多期分类栅格文件路径 files - c( class_2000.tif, class_2010.tif, class_2020.tif ) area_list - list() for (i in seq_along(files)) { r - rast(files[i]) tab - freq(r) # 像元频数统计 cell_km2 - res(r)[1] * res(r)[2] / 1e6 # 像元面积换算成km² area_list[[i]] - data.frame( Year c(2000, 2010, 2020)[i], Class tab$value, Area_km2 tab$count * cell_km2 ) } area_data - bind_rows(area_list)再把Fragstats导出的.land或.class文件读进来。R里直接read.csv()就能读因为Fragstats导出的CSV第一列是类别或景观的名称后面是各类指数值。多个年份的指数文件可以用同样的方式循环读取然后横向合并frag_2000 - read.csv(fragstats_2000.land.csv) frag_2010 - read.csv(fragstats_2010.land.csv) frag_2020 - read.csv(fragstats_2020.land.csv) frag_all - bind_rows( list(Year2000, frag_2000), list(Year2010, frag_2010), list(Year2020, frag_2020) )到这里一个年份 x 类别面积 景观指数的总表就建好了后面所有图表都从这张表取数。这个整合过程才是R在这套流程里真正的价值——它让三套软件的成果汇到一个地方也让你在修改数据时只需改一处不用来回切软件。4.3 ggplot2出图模板与300dpi输出数据整理好了出图就是水到渠成的事。ggplot2的出图逻辑一句话概括数据框里的列映射到图形的视觉元素上。环评里最常用的是两个模板。第一个是土地利用面积变化堆叠柱状图看各类别面积随时间怎么变。第二个是景观指数变化折线图看破碎化、多样性怎么变。两个模板够覆盖绝大多数生态专题图件。library(ggplot2) # 土地利用面积变化图 p1 - ggplot(area_data, aes(x factor(Year), y Area_km2, fill Class)) geom_col(position stack) scale_fill_brewer(palette Set3) theme_minimal(base_size 14) labs(x NULL, y 面积km², title 土地利用面积变化) ggsave(图1 土地利用变化图.png, p1, width 8, height 5, dpi 300)出图时有三个细节必须注意。第一是中文字体问题。ggplot2默认字体在中文环境下经常显示成方块我习惯在脚本开头用showtext包加载系统中文字体比如微软雅黑或思源黑体一行代码解决乱码。第二是出图尺寸和分辨率。环评报告印刷一般要求300dpi宽度控制在8英寸左右这样图面既清晰又不会被排版软件压缩变形。第三是批量输出。多个年份、多个指标的图用for循环配合ggsave循环保存文件名用拼接字符串比如图2 NP变化趋势.png、图3 SHDI变化趋势.png一套脚本全出。如果要做多张图拼图推荐patchwork包操作非常直观library(patchwork) combined - p1 p2 plot_layout(ncol 1) ggsave(图4 组合图.png, combined, width 8, height 10, dpi 300)4.4 大数据栅格处理与经典内存报错最后聊几个R在处理大影像时常见的报错这些坑虽然不起眼但撞上一次就能卡住你半天。第一个是protect(): protection stack overflow。这个报错常见于一次性把大量高分辨率栅格转成向量或者data.frame的场景。我处理高分影像转点时撞过解决思路无非三条一是调大R的内存限制Linux下启动前用环境变量R_MAX_VSIZE设置上限二是把栅格分块处理terra里用chunk方式分批计算后再汇总三是直接把栅格重采样到更低分辨率再分析。环评项目的分析精度要求下第三种办法效率最高。第二个是cannot allocate vector of size。这个报错一般意味着数据量超出了物理内存尤其在读入0.5m高分影像时很容易爆。稳妥做法是统计栅格时直接用terra的freq函数不要让R先把整个栅格转成向量再统计等于少走一步内存翻倍的路。第三个是RStduio运行时偶尔出现的服务权限类报错比如Docker里挂载R服务时提示权限不足。这种问题多半是容器目录权限没配好和环评数据本身没关系检查一下目录权限即可不要花时间在R里找原因。我的原则是能用terra直接在栅格层面算的绝不先转向量再算能分块处理的绝不一次性全部读入内存。数据量大不是问题问题是你有没有用对处理方式。最后再分享一个我自己的习惯每出一个版本的图件我都会把原始影像时相、分类精度表、Fragstats参数设置、R脚本一起存成一份图件数据清单附在报告附件里。评审会上被问这版图的数据怎么来的打开清单一分钟就能讲清楚。这套流程跑顺之后你会发现环评制图真正拼的不是手速是整个过程的可追溯性。祝大家早日摆脱熬夜勾图的循环。
环评生态制图三件套：ENVI遥感解译+Fragstats景观指数+R批量出图
环评项目里最耗时间的往往不是预测模型本身而是生态专题那一摞图件。我印象很深的一次某公路项目生态评价专章要求出土地利用现状图、植被覆盖度图、生态系统类型图、景观格局指数表外加一张几十年的趋势对比图。当时还是老思路靠ArcGIS一亩一亩人工勾绘光两个图斑就折腾了两天图例顺序和颜色体系还得反复调。后来我把整个工作流换成了ENVIRFragstats三件套ENVI做遥感解译Fragstats算景观格局指数R负责数据汇总和批量出图整个制图环节从一周压缩到一天半图表质量反而更稳定。这篇就把这套流程完整捋一遍适合正在编制环评报告书、尤其是生态影响评价专章的同行参考。想直接跑通的人按下面的章节顺序操作就行每步我都写了参数和避坑点。1. 环评图件到底有哪些先拆清单再定工具分工1.1 报告书里躲不开的图件从基础信息图到生态专题图环评报告书里的图件数量新手往往低估了。按工作内容拆开看其实就三类。第一类是基础信息图。包括项目地理位置图、评价范围图、环境监测点位图、敏感目标分布图、总平面布置图。这些图的核心是矢量化底图加注记ArcGIS或QGIS就能干跟遥感解译关系不大图面整饰也就那样不再展开。第二类是现状评价图这是图件工作量的大头。土地利用现状图、植被类型图、植被覆盖度图、生态系统类型图、土壤侵蚀等级图基本都依赖遥感影像解译得到基础数据。这类图恰恰是ENVI的主场从影像预处理到分类出图一条链路走下来数据来源清清楚楚。第三类是预测和分析类图。大气浓度等值线图、声环境预测图、地下水和土壤影响图以及生态影响评价里的景观格局对比图、生态质量变化趋势图。这些图的核心是统计计算和可视化R语言做这类生产工具最顺手尤其当你要对十几年多期数据做批量对比时Excel根本撑不住。三类图放一起看真正消耗人力的就是现状评价图和趋势对比图而这部分正好能交给ENVIRFragstats这套组合。很多人一上来就埋头在ArcGIS里画画完了数据来源说不清、精度验证没做评审会上被问一句就哑火。我的建议是先列一张图件清单按数据生产链路去分工而不是按图件类型去逐个手画。1.2 为什么是ENVIRFragstats三条流水线如何分工先明确一个逻辑这套三件套不是三个软件各干各的而是一条生产流水线。ENVI负责第一段把原始影像变成分类栅格。它的辐射定标、大气校正、监督分类、精度评价都是菜单化操作比命令行写GDAL门槛低很多。环评工程师的时间不应该耗在写底层代码上ENVI在学习成本和处理效率之间平衡得最好。如果只是做常规生态本底调查不需要复杂编程ENVI足够。Fragstats负责第二段把分类栅格变成景观格局指数。这个软件几乎是景观生态学领域的行业标准指数全、算法公开评审专家认可度高。环评导则里要求描述生态系统结构、功能和破碎化程度时用Fragstats算出来的NP、PD、SHDI这些指标比你自己在Excel里拍脑袋算的要有说服力得多。R语言负责第三段把所有数据汇总成报告能用的图表。它的统计能力和可视化能力正好补上前面两个软件的短板把多年土地利用面积堆在一起做对比图、把Fragstats输出的一堆表格清洗成一张主表、批量输出300dpi的PNG这些都是R的强项。下面三章就按这个流水线顺序展开每一步的实操细节我都会写到。2. ENVI出分类图从原始影像到一份能进报告的土地利用现状图2.1 影像数据源选择Landsat、哨兵与高分的取舍环评制图的第一步不是打开ENVI而是选影像。选错数据源后面做得再漂亮都是白费。判定标准主要看评价范围尺度和精度要求。大区域的线性工程比如公路、铁路、管道评价范围动辄几十平方公里甚至上百用30m分辨率的Landsat或者10m分辨率的哨兵2号就够了两者都可以免费下载。Landsat的优势是历史数据长可以追溯到上世纪八十年代做生态趋势对比非常好用哨兵2号分辨率更高更新频率快适合近几年的本底调查。小范围的厂区、矿山、风电场评价范围只有几平方公里这时必须上高分影像。GF-1、GF-2或者商业的高分影像能做到亚米级或米级分辨率能看清小尺度的地类边界。这类影像成本高一些但生态底图的精度是评审专家盯着的点不能省。还有一类特殊情况矿区的生态影响评价。如果项目涉及尾矿库、采空区需要提取植被胁迫或异常信息可以考虑使用GF-5高光谱数据做蚀变信息提取这是ENVI的光谱分析强项但常规环评图件用不上属于进阶技能这里先带过。选影像还要注意两个细节。一是时相尽量选6到9月的生长季植被光谱特征最明显分类精度最高。二是云量控制在10%以下有云的区域直接放弃重选。跨景拼接前先确认各景影像是否已经做过几何校正时相对不齐的影像拼出来容易出现明显接边痕迹。2.2 预处理四连辐射定标、大气校正、几何校正与裁剪拿到原始影像别急着分类先把预处理走完。第一是辐射定标。原始影像存的DN值是个无量纲数字需要变成辐射亮度或反射率才能做后续计算。ENVI里的路径是工具箱Radiometric Calibration选择多光谱文件后直接输出定标结果这个操作没有太多参数要调但一定要做不做后面的大气校正就无从谈起。第二是大气校正。大气中的水汽和气溶胶会吸收和散射电磁波导致地物反射率失真尤其是植被指数计算不做大气校正NDVI值普遍偏低。ENVI里常用FLASH模块需要填的传感器类型、影像中心经纬度、飞行高度、成像时间以及大气模型和气溶胶模型。大气模型按成像季节和纬度选中纬度地区夏季选Mid-Latitude Summer冬季选Mid-Latitude Winter。如果时间紧不想细调参数也可以用QUAC快速大气校正精度会差一些但环评本科底图够用。第三是几何校正。Landsat L2级产品自带精确几何定位不需要额外处理。高分影像如果发现地物边界和矢量底图有偏移就要手动选控制点做一次几何校正控制点尽量选在道路交叉口、水库坝体这类特征明显的点上均匀分布RMSE控制在0.5个像元以内。最后是裁剪。用评价范围的shp边界在ENVI里执行Subset Data from ROIs/Shapefile输出自带坐标信息的GeoTIFF。裁剪范围建议比评价范围外扩1到2公里给后续制图留出修边余地否则评价边界紧贴影像边缘图面上很难看。2.3 监督分类与精度验证ROI怎么选、SVM怎么用预处理完成后进入核心环节土地利用分类。我推荐监督分类里用SVM支持向量机RBF核小样本情况下比最大似然分类稳定对椒盐噪声的抑制也好一些。分类体系要提前定好结合环评导则和项目特点常见的类别是耕地、林地、草地、建设用地、水域、裸地这六类有特殊项目再增加园地、沼泽、盐田等。类别编码必须从一开始就固定耕地1、林地2、草地3、建设用地4、水域5、裸地6后期多期对比全靠编码对齐。ROI选取是整个分类流程里最考验经验的一步。每个类别至少选取30到50个样本样本要覆盖同一地类的不同光谱表现。比如耕地里有水田和旱地光谱差异很大不能只在一个区域点样本就完事。ROI尽量画成闭合多边形不要只点一个像元点选像元容易选到混合像元分类时容易带偏。训练样本和验证样本要分开。ENVI里可以通过样本分割工具把ROI按比例拆成训练集和验证集我一般按3比1拆分。分类完成后用验证样本计算混淆矩阵重点看总体精度和Kappa系数。环评报告里总体精度85%以上、Kappa大于0.7基本就能应对评审。如果达不到不要急着调分类器先回去补样本样本质量上来了精度自然就上来了。分类完成后还要做后处理。ENVI里的Majority/Minority分析3x3窗口把分类图里孤立的小斑块归到周围主要类别里能有效去掉椒盐效应。跑完之后建议再人工快速扫视一遍重点检查大水体、城镇边界有没有明显错分。2.4 多期分类的统一性问题分类体系、时相与分辨率对齐做生态趋势分析的同行一定会遇到多期影像分类的场景这里有几个坑必须提前踩平。最大的坑是分类体系不一致。不同年份的影像由不同的人分类别定义稍有差异后面对比就全是噪声。我的一般做法是在一张Excel里把每个年份的类别编码表固定下来分类时严格按统一编码执行任何人接手都先看这张表。第二个坑是分辨率不一致。早期Landsat影像只有30m后来用了10m的哨兵如果直接对比面积地类边界会因像元大小不同产生系统偏差。处理办法是统一重采样到相同分辨率我通常以较低分辨率30m为基准高分辨率影像重采样下来再参与面积统计。第三个坑是时相不对齐。不同年份影像如果成像月份差很多植被光谱差异会被误判成土地利用变化。做多年对比时尽量选同一季节的影像如果实在找不到同期影像报告里要明确说明时相差异带来的不确定性别等评审专家来问。每期分类完成后记得导出三个东西分类GeoTIFF、分类结果图以及分类前的原始影像路径和成像日期。后面写报告的数据来源说明时这些东西都是要用的。3. Fragstats算指数参数设错整张景观格局表就废了3.1 输入文件准备把ENVI分类结果喂给Fragstats之前必须做的事ENVI分类结果并不能直接丢给Fragstats中间还差一步格式整理。Fragstats 4.2支持GeoTIFF格式输入但要求像元值是连续的整数类别编码。ENVI分类输出的栅格像元值一般已经是从1开始的整数但可能存在一些非分类像元比如背景区域被赋值成0或者255。最关键的一步就是把这些背景值单独处理掉我会建议统一把背景设为0并且保证类别编码从1开始连续编号类别数量不要超过几十个。这一步在ENVI里可以用Band Math或重分类工具做但我个人更习惯在R里用terra包处理因为可以写成脚本批量跑多期数据一次搞定。比如用subst函数把特定的旧值替换成新的类别编码再统一写回GeoTIFF。处理完的栅格在Fragstats里加载时还要顺手看一眼属性表确认没有出现类别值跳号的情况。3.2 三个关键参数背景值、邻域规则和边缘深度Fragstats主界面的参数看着简单但三个参数一旦设错算出来的指数全盘作废这个坑我踩过太多次了。第一个是背景值设置。在Class Properties里明确指定哪个值代表背景。如果不设背景像元会被当成一个类别参与面积和密度计算NP、PD会全部偏大LPI也可能被背景斑块主导。这是我见过同行最常犯的错误很多人跑完结果就直接粘到报告里完全没意识到指数已经被背景污染了。第二个是邻域规则。4邻域只考虑上下左右四个方向8邻域把对角方向也算进来。这个选择直接影响聚合度、连接度、散布与并列指数等一大类指标的结果。环评里我一般用8邻域因为像元间的生态过程往往是八方向连通的但报告的方法部分必须把这个选择写清楚评审专家会看。第三个是边缘深度。默认设置是1个像元够用。但如果要算核心区面积、核心区密度这些核心区指标就必须根据真实生态影响距离来设置。比如某生态敏感鸟类的影响距离是50m30m影像下50m换算成边缘深度就是2个像元那就要把边缘深度设为2。设置不当核心区面积会严重高估。参数设完后输出层面建议同时选class和landscape两个层级格式选CSV这样每个指标按类别算了一份同时也有一份全景观层面的汇总。导出的.land、.class文件直接用Excel或R打开即可。3.3 报告里真正用得到的景观指数一张表讲清生态含义Fragstats能算几十个指数但环评报告里真正用得到的其实就七八个。选指数的逻辑不是越多越好而是每个指数都要能对应到生态学含义并且能说出生态过程。下面这张表是我在项目里经常用到的核心指数清单指数缩写全称生态学含义环评里的用法NP斑块数量景观破碎化程度数值上升说明生境被切割PD斑块密度单位面积斑块数与NP配合消除面积影响LPI最大斑块指数优势斑块占比反映基质类型和优势度LSI景观形状指数斑块形状复杂度形状越复杂边缘效应越强CONTAG蔓延度景观连通程度下降通常说明连通性变差AI集聚度同类斑块聚集程度说明类别空间分布偏好SHDI香农多样性景观异质性上升说明类型变多、结构复杂SHEI香农均匀度各类型分布均匀程度辅助SHDI解释我在实际报告中写生态影响评价时一般这样组织先把项目区各类用地的面积占比列出来再用NP、PD说明破碎化用CONTAG、AI说明连通性用SHDI、SHEI说明景观多样性变化。每写一个指数后面都要跟一句生态过程解释。比如PD上升说明生境斑块被进一步切割这对于移动能力弱的物种潜在影响更大不能光罗列数字。3.4 多期景观格局对比的实操顺序多期对比的实操顺序其实很简单把每期分类图按3.1节的方法处理成标准输入然后在完全相同的参数设置下分别计算指数。注意参数一定要一致背景值、邻域规则、边缘深度都设成相同值否则不同年份的指数差异包含了方法差异对比就没有意义了。算完多期指数后我一般会建一张主表行是年份列是指标再辅助上一张R出的折线图或柱状图把NP、SHDI这些关键指标的变化趋势直观呈现出来。Fragstats本身不提供可视化这正好接上下一章R语言的内容。评审会上被问得最多的就是这些指数变化的生态学意义是什么所以我建议在出图前先把生态过程解释写好再让图表配合文字而不是先出图再想解释。这也是整个Fragstats环节最核心的写作逻辑。4. R把三件套串成流水线数据整合、ggplot2批量出图与内存坑4.1 R环境与包管理从安装到没有程辑包报错R环境的搭建没有捷径但也就花二十分钟。先在CRAN官网下载安装R然后装一个RStudio日常操作都在RStudio里完成。装包是新人第一道坎。install.packages(包名)是最常用的命令如果遇到不存在叫‘xxx’这个名字的程辑包的报错先别怀疑人生按三步排查第一步查包名拼写R里下划线、大小写都区分第二步查镜像设置Tools - Global Options - Packages里把CRAN镜像换到国内镜像不然下载经常断第三步确认这个包是不是在GitHub上发布的如果是先安装devtools再用install_github(作者/包名)安装。RStudio里我用得最勤的快捷键是CtrlEnter运行当前行和CtrlShiftR插入代码段这两个键记住了写长脚本会顺很多。环评生态分析常用的包有terra和sf处理栅格矢量dplyr和tidyr做数据清洗ggplot2和patchwork出图readxl和writexl读写Excel。少数生态学专用包在Bioconductor仓库里发布需要先安装BiocManager再用BiocManager::install(包名)。这一整套装完后面几乎不会再碰包安装问题。4.2 数据整合栅格统计、Fragstats指数合并成一张总表数据整合是R在这套流水线里最核心的活儿。我用一个具体场景演示三个年份的土地利用分类栅格需要统计每个年份各类别面积再把Fragstats算出的景观指数一并合并。先看栅格统计。用terra包批量读取多期分类结果统计每个类别像元数乘以像元面积换算成平方公里代码很简洁library(terra) library(dplyr) # 多期分类栅格文件路径 files - c( class_2000.tif, class_2010.tif, class_2020.tif ) area_list - list() for (i in seq_along(files)) { r - rast(files[i]) tab - freq(r) # 像元频数统计 cell_km2 - res(r)[1] * res(r)[2] / 1e6 # 像元面积换算成km² area_list[[i]] - data.frame( Year c(2000, 2010, 2020)[i], Class tab$value, Area_km2 tab$count * cell_km2 ) } area_data - bind_rows(area_list)再把Fragstats导出的.land或.class文件读进来。R里直接read.csv()就能读因为Fragstats导出的CSV第一列是类别或景观的名称后面是各类指数值。多个年份的指数文件可以用同样的方式循环读取然后横向合并frag_2000 - read.csv(fragstats_2000.land.csv) frag_2010 - read.csv(fragstats_2010.land.csv) frag_2020 - read.csv(fragstats_2020.land.csv) frag_all - bind_rows( list(Year2000, frag_2000), list(Year2010, frag_2010), list(Year2020, frag_2020) )到这里一个年份 x 类别面积 景观指数的总表就建好了后面所有图表都从这张表取数。这个整合过程才是R在这套流程里真正的价值——它让三套软件的成果汇到一个地方也让你在修改数据时只需改一处不用来回切软件。4.3 ggplot2出图模板与300dpi输出数据整理好了出图就是水到渠成的事。ggplot2的出图逻辑一句话概括数据框里的列映射到图形的视觉元素上。环评里最常用的是两个模板。第一个是土地利用面积变化堆叠柱状图看各类别面积随时间怎么变。第二个是景观指数变化折线图看破碎化、多样性怎么变。两个模板够覆盖绝大多数生态专题图件。library(ggplot2) # 土地利用面积变化图 p1 - ggplot(area_data, aes(x factor(Year), y Area_km2, fill Class)) geom_col(position stack) scale_fill_brewer(palette Set3) theme_minimal(base_size 14) labs(x NULL, y 面积km², title 土地利用面积变化) ggsave(图1 土地利用变化图.png, p1, width 8, height 5, dpi 300)出图时有三个细节必须注意。第一是中文字体问题。ggplot2默认字体在中文环境下经常显示成方块我习惯在脚本开头用showtext包加载系统中文字体比如微软雅黑或思源黑体一行代码解决乱码。第二是出图尺寸和分辨率。环评报告印刷一般要求300dpi宽度控制在8英寸左右这样图面既清晰又不会被排版软件压缩变形。第三是批量输出。多个年份、多个指标的图用for循环配合ggsave循环保存文件名用拼接字符串比如图2 NP变化趋势.png、图3 SHDI变化趋势.png一套脚本全出。如果要做多张图拼图推荐patchwork包操作非常直观library(patchwork) combined - p1 p2 plot_layout(ncol 1) ggsave(图4 组合图.png, combined, width 8, height 10, dpi 300)4.4 大数据栅格处理与经典内存报错最后聊几个R在处理大影像时常见的报错这些坑虽然不起眼但撞上一次就能卡住你半天。第一个是protect(): protection stack overflow。这个报错常见于一次性把大量高分辨率栅格转成向量或者data.frame的场景。我处理高分影像转点时撞过解决思路无非三条一是调大R的内存限制Linux下启动前用环境变量R_MAX_VSIZE设置上限二是把栅格分块处理terra里用chunk方式分批计算后再汇总三是直接把栅格重采样到更低分辨率再分析。环评项目的分析精度要求下第三种办法效率最高。第二个是cannot allocate vector of size。这个报错一般意味着数据量超出了物理内存尤其在读入0.5m高分影像时很容易爆。稳妥做法是统计栅格时直接用terra的freq函数不要让R先把整个栅格转成向量再统计等于少走一步内存翻倍的路。第三个是RStduio运行时偶尔出现的服务权限类报错比如Docker里挂载R服务时提示权限不足。这种问题多半是容器目录权限没配好和环评数据本身没关系检查一下目录权限即可不要花时间在R里找原因。我的原则是能用terra直接在栅格层面算的绝不先转向量再算能分块处理的绝不一次性全部读入内存。数据量大不是问题问题是你有没有用对处理方式。最后再分享一个我自己的习惯每出一个版本的图件我都会把原始影像时相、分类精度表、Fragstats参数设置、R脚本一起存成一份图件数据清单附在报告附件里。评审会上被问这版图的数据怎么来的打开清单一分钟就能讲清楚。这套流程跑顺之后你会发现环评制图真正拼的不是手速是整个过程的可追溯性。祝大家早日摆脱熬夜勾图的循环。
