在 Big Bend 国家公园这类半干旱生态保护区里推土机清表造成的地表扰动往往不是一次性事件而是一个从无到有、从边缘到中心的过程。人工巡护很难在第一时间发现偏远地段的植被移除和土壤裸露等肉眼可见时破坏可能已经持续了数天甚至数周。卫星遥感提供了一条可回溯、可量化、可复核的检查路径通过对比扰动前后的多光谱影像把“突然消失的植被”和“突然暴露的土壤”变成一张带坐标的扰动图。下面以 Big Bend 国家公园的典型场景为例介绍如何用 Sentinel-2 影像和 NDVI 差值完成一次最小可复现的保护区地表扰动检测并给出从实验脚本走向生产监测系统的方案。1. 保护区地表扰动为什么适合用遥感发现1.1 人工巡护的盲区保护区面积大地形复杂很多区域车辆无法到达。靠巡护员开车或徒步巡查通常只能覆盖道路附近和定期巡护线路。推土机作业通常不会发生在入口附近而是选在偏远、隐蔽、接近边界的位置。等巡护员发现异常时现场可能已经被简单恢复或者已经错过最佳取证时间。遥感解决的不是“人能不能看见”的问题而是“机器能不能持续盯住每一块土地”的问题。卫星按固定时间重复经过同一区域历史影像可以随时回放。只要有扰动前和扰动后的数据就能在办公室先把疑似区域缩小到几个多边形再安排巡护员带着坐标去现场核查。1.2 NDVI 如何把“扰动”变成可计算的数值NDVI 是归一化植被指数公式是NDVI (NIR - Red) / (NIR Red)对于 Sentinel-2 影像NIR 用 B8 近红外波段Red 用 B4 红光波段。活体植被会大量吸收红光、强烈反射近红外所以 NDVI 高。推土机清表后植被被移除土壤或碎石暴露红光反射上升近红外反射下降NDVI 会明显降低。在 Big Bend 这类半干旱地区天然植被本身就不算茂密NDVI 的绝对值不会像森林那么高。但这不代表方法失效因为扰动检测关注的是“同一个位置从上一期到下一期突然变低”而不是和全球平均值比较。一个像元如果从 0.35 掉到 0.08即使绝对值不高也是一个强烈的变化信号。1.3 扰动检测的关键是找突变而不是找低值很多初次接触遥感的人会直接提取“NDVI 小于某个值”的区域作为扰动图。这在保护区场景里会出问题干涸河床、裸岩、沙地天然就是低 NDVI但它们不是推土机造成的。正确的做法是找“突变”前一期影像中 NDVI 正常。后一期影像中 NDVI 明显下降。同一区域没有降水、火灾、季节更替等自然解释。比较两个时相的 NDVI本质上是在做时间维度的变化检测。突变的幅度越大、连续像元越多越符合机械扰动的特征。注意遥感变化检测的结果只能说明“地表出现了像扰动一样的变化”不能单独作为执法或行政依据。最终结论需要结合历史影像、管理记录和现场核查。2. 数据选型Sentinel-2、波段组合与工具链2.1 Sentinel-2 的优势分辨率、重访周期、免费开放Sentinel-2 是目前做地表扰动检测最常用的开放中分辨率数据源。它有 13 个波段其中红、绿、蓝、近红外 4 个波段空间分辨率达到 10 米适合识别推土机清表产生的大范围地表变化。数据特性数值对保护区监测的意义空间分辨率10 米 / 20 米 / 60 米能识别几十米宽的机械作业带重访周期单星 10 天双星约 5 天可以捕捉周尺度突变辐射定标L2A 地表反射率产品不同日期影像可比性较好波段数量13 个支持 NDVI、BSI、NDMI 等多个指标获取成本免费开放适合长期持续监测实际项目中需要确认拿到的是 L2A 地表反射率产品。如果只有 L1C 产品需要先做大气校正否则不同日期影像之间的 NDVI 差异会受到大气条件干扰。2.2 波段选择和指数公式针对推土机扰动重点使用以下波段波段通用名称近似中心波长分辨率在扰动检测中的作用B2Blue490 nm10 米计算裸土指数 BSIB4Red665 nm10 米计算 NDVI、BSIB8NIR842 nm10 米计算 NDVI、NDMIB11SWIR11610 nm20 米计算 BSI、NDMI需要重采样到 10 米除 NDVI 外常用的辅助指标有两个裸土指数 BSI扰动后土壤暴露越明显BSI 越高。归一化差异水分指数 NDMI植被被移除后地表含水量和植被含水量下降NDMI 通常降低。单看 NDVI 容易漏掉部分变化加入 BSI 和 NDMI 可以让误报判断更可靠。2.3 本案例的坐标系与预处理顺序Big Bend 国家公园位于德克萨斯州西部使用 UTM 投影时通常对应 UTM 13N 分带。实际处理时建议把所有栅格统一到一个投影坐标系例如 EPSG:32613。预处理顺序建议下载同一区域、同一轨道号的 Sentinel-2 L2A 产品。检查前后两个时相是否包含云覆盖。将 20 米波段重采样到 10 米网格。裁剪到保护区边界外扩缓冲区。根据 SCL 波段生成云、云影、水体掩膜。计算 NDVI、BSI、NDMI。不要在还没做云掩膜时就计算差值否则云影会造成大范围 NDVI 突降误报。3. 最小实现两期 NDVI 差值定位推土机扰动3.1 准备目录和依赖下面的示例基于本地 Python 环境。先创建虚拟环境并安装依赖mkdir ndvi-disturbance cd ndvi-disturbance python -m venv venv source venv/bin/activate pip install rasterio numpy pandas matplotlib如果已经有真实影像目录结构可以这样组织ndvi-disturbance/ ├── before/ │ ├── B04.tif │ ├── B08.tif │ └── SCL.tif ├── after/ │ ├── B04.tif │ ├── B08.tif │ └── SCL.tif └── scripts/ └── detect_disturbance.py如果暂时没有真实影像可以先准备两张模拟数组验证流程再替换成真实影像。下面的代码以单波段 tif 文件为输入。3.2 计算 NDVI 并生成云掩膜先写一个读取波段的函数再计算 NDVIimport rasterio import numpy as np def read_band(path): with rasterio.open(path) as src: return src.read(1).astype(float32), src.profile def calc_ndvi(b4_path, b8_path): b4, profile read_band(b4_path) b8, _ read_band(b8_path) denom b8 b4 with np.errstate(divideignore, invalidignore): ndvi (b8 - b4) / np.where(denom 0, np.nan, denom) return np.clip(ndvi, -1, 1), profile执行ndvi_before, profile calc_ndvi(before/B04.tif, before/B08.tif) ndvi_after, _ calc_ndvi(after/B04.tif, after/B08.tif)接着读取 SCL 场景分类波段排除云、云影等干扰像元def make_land_mask(scl_path): with rasterio.open(scl_path) as src: scl src.read(1) # 4 植被5 裸土6 水体 # 3 云影7/8/9/10 云和卷云 return np.isin(scl, [4, 5, 6]) mask_before make_land_mask(before/SCL.tif) mask_after make_land_mask(after/SCL.tif) valid mask_before mask_afterSCL 波段是 L2A 产品自带的场景分类结果。实际项目中如果 SCL 没有覆盖到研究区也可以使用云概率波段或 Fmask 算法生成掩膜。3.3 差值制图与面积统计两期 NDVI 相减再设置变化阈值delta_ndvi np.where(valid, ndvi_after - ndvi_before, np.nan) # 阈值需要根据区域标定-0.25 只是示例 threshold -0.25 candidate np.full(delta_ndvi.shape, 0, dtypeuint8) candidate[(delta_ndvi threshold) (delta_ndvi -1.0)] 255 pixel_area_m2 10.0 * 10.0 pixels int((candidate 255).sum()) area_ha pixels * pixel_area_m2 / 10000.0 print(疑似扰动像元数:, pixels) print(疑似扰动面积(ha):, area_ha)执行后正常输出类似疑似扰动像元数: 3860 疑似扰动面积(ha): 38.60这里的面积是基于投影坐标网格计算的。只要影像统一在 UTM 投影下10 米分辨率的单像元面积就是 100 平方米乘像元数再除以 10000 就是公顷数。3.4 输出扰动候选点坐标将结果写成 GeoTIFFprofile.update(dtypeuint8, count1, nodata0) with rasterio.open(disturbed.tif, w, **profile) as dst: dst.write(candidate, 1)再导出候选点坐标便于导入地图或巡护终端with rasterio.open(disturbed.tif) as src: data src.read(1) transform src.transform rows, cols np.where(data 255) # 抽样输出避免点密度过高 xs, ys rasterio.transform.xy(transform, rows[::50], cols[::50]) for x, y in zip(xs, ys): print(f{x:.2f}, {y:.2f})如果业务系统需要 WGS84 经纬度可以再用 pyproj 把 UTM 坐标转成经纬度。坐标转换应该在输出给巡护人员之前完成。4. 阈值不能乱拍在 Big Bend 场景下调参的思路4.1 半干旱区 NDVI 基线为什么偏低Big Bend 地处奇瓦瓦沙漠生态区植被以灌丛、牧草和稀疏林地为主。很多像元的 NDVI 天然只有 0.2 到 0.4河床和裸岩可能接近 0。如果照搬湿润地区常用的“NDVI 降到 0.3 以下就是扰动”这类规则会得到大量假阳性。正确思路是看变化量而不是绝对值。扰动前 NDVI 0.35扰动后 0.35不一定没问题扰动前 0.35扰动后 0.10即便绝对值不高也应列为疑似扰动。4.2 用统计量替代绝对阈值最简单的方法是把全图 NDVI 变化量的分布画出来看看“异常尾部”在哪里。对低植被区域可以用中位数和绝对中位偏差 MAD 构造动态阈值median np.nanmedian(delta_ndvi) mad np.nanmedian(np.abs(delta_ndvi - median)) # 1.4826 是为了让 MAD 近似正态分布标准差 sigma 1.4826 * mad k 5.5 threshold median - k * sigma print(动态阈值:, threshold)这种方法适合快速探索但不能直接写成生产规则。生产环境需要按土地覆盖类型、季节、历史基线分别标定阈值否则一场正常降水后的植被变绿也可能被当成“变化”。动态阈值的优势在于它不依赖某个固定数值而是依赖同一期影像内部的空间统计关系。4.3 排除火灾、物候和水体干扰NDVI 大幅下降并不一定都是推土机作业。常见干扰包括火烧迹地火灾后植被烧毁NDVI 下降但表面会残留炭化物质短波红外特征和机械翻土不同。物候变化雨季前草场枯黄时NDVI 会自然下降。水体淹没或干涸水体会导致近红外反射率明显变化。因此建议同时计算 BSI 和 NDMI形成多指标判断指标公式机械扰动后的典型变化NDVI(B8 - B4) / (B8 B4)下降BSI((B11 B4) - (B8 B2)) / ((B11 B4) (B8 B2))上升NDMI(B8 - B11) / (B8 B11)下降一次疑似机械扰动应该同时满足NDVI 突降、BSI 上升、NDMI 下降并且变化区域有一定空间连续性。如果只是 NDVI 下降而 BSI 没有明显变化要先考虑物候或云影干扰。在低植被区不要使用其他地区常见的 NDVI 阈值。先画出扰动前影像覆盖区的 NDVI 直方图再决定变化幅度门槛。5. 从单张异常图到疑似扰动工单完整判断链路5.1 一条示例记录应该包含哪些字段从影像中检测出候选像元后不能直接把栅格丢给业务人员。应该按“疑似扰动斑块”输出结构化记录字段示例值说明obs_idBD-20240612-01扰动斑块唯一编号longitude-103.25斑块质心经度latitude29.20斑块质心纬度area_ha2.36斑块面积delta_ndvi-0.31斑块平均 NDVI 变化bsi_change0.12斑块平均 BSI 变化valid_ratio0.97有效像元占比satelliteS2B数据源卫星date_before2024-05-20扰动前影像日期date_after2024-06-14扰动后影像日期confidencehigh置信度review_statuspending复核状态有了这些字段巡护人员拿到的不再是一张灰度图而是一条可以导航、可以排序、可以回传到管理系统的任务单。5.2 现场复核前先做空间过滤栅格检测结果通常包含大量孤立像元。推土机作业不会是单个像元大小通常会有带状或块状的连续区域。建议在生成多边形前做三步过滤去掉小于最小成图单元的斑块。如果使用 10 米分辨率0.5 公顷大约对应 50 个像元可以作为初筛尺度。排除道路缓冲区和已知防火带。这类线性地物也会造成 NDVI 变化但属于正常管理活动。对二值图做形态学开运算消除孤立噪点。空间过滤可以减少现场核查工作量也能避免把一个自然裸地碎斑当成异常报警。5.3 时间序列确认连续影像比单期更可靠两期影像对比只能证明“变化发生了”还不能确认“变化是否突然发生”。更好的做法是把三期、四期影像放在一起看扰动前第 1 期NDVI 0.35。扰动前第 2 期NDVI 0.32。扰动后第 1 期NDVI 0.08。扰动后第 2 期NDVI 0.05。如果前两期稳定第三期突然掉到 0.1 以下这就是非常典型的机械扰动信号。如果四期 NDVI 一直在缓慢下降则更可能是干旱或病虫害。判断时可以用时间序列的“跳变”特征做规则previous_stable abs(ndvi_series[1] - ndvi_series[0]) 0.05 sudden_drop ndvi_series[2] - ndvi_series[1] -0.20 if previous_stable and sudden_drop: print(疑似突然地表扰动)这个规则不复杂但能过滤掉相当多物候和降雨造成的假阳性。生产系统可以在哨兵影像新数据到达后自动追加一期指数形成滚动时间序列。6. 常见问题与排查清单6.1 五个高频问题问题现象常见原因检查方式处理建议云影区域出现大量 NDVI 下降云掩膜未生效或 SCL 不可用查看 SCL 影像和 RGB 真彩色合成补充云概率掩膜或使用时序合成影像整个研究区 NDVI 都下降前后影像季节不同比较两期影像获取日期和物候状态选择同季节、同生育期影像重新对比阈值调高后有大量误报固定阈值不适用于低植被区绘制 NDVI 变化直方图改用中位数 MAD 动态阈值道路和防火带被识别为扰动线性地物也存在 NDVI 下降叠加道路、步道、管理边界矢量增加道路缓冲区排除设置最小成图单元明显的推土区没有被识别变化区域小于混合像元或 NDVI 变化不明显查看高分辨率影像和原始波段组合使用更高分辨率数据或引入 BSI / NDMI 多指标投票6.2 排查顺序从数据到结论遇到异常结果按下面的顺序排查检查前后两期影像的日期、轨道号、覆盖范围是否一致。检查 SCL 云掩膜是否把真实地表像元错误剔除。检查 NDVI 计算后的数值范围是否正常。检查差值分布中疑似扰动的像元是否集中在同一区域。叠加保护区边界、道路、水系等矢量数据。查看历史高分辨率影像确定是否存在真实作业痕迹。最后再到现场核查。不要在第一步就怀疑算法。遥感检测的大部分“神秘误报”最后都能追溯到数据问题或预处理问题。7. 从实验脚本到生产监测系统的落地建议7.1 一个可工作的定时流程实验脚本跑通后生产监测系统可以按如下频率运行每天或每周检查 Sentinel-2 是否有覆盖研究区的新 L2A 产品。下载新影像并裁剪到保护区边界。执行云掩膜、重采样、指数计算。将最新一期 NDVI 与最近 3 到 5 期历史影像比较。如果检测到疑似扰动斑块生成 GeoJSON 或 Shapefile。将结果写入 PostGIS 数据库并推送消息通知。业务人员在 WebGIS 上查看斑块标记“待核查”“已核查”“误报”。小团队不要一开始就上大数据平台。先用定时任务、本地数据库和消息通知跑几个季度确定算法稳定后再考虑迁移到云环境。7.2 学习环境与生产环境的区别维度学习验证生产监测影像获取手工下载几期自动订阅或定时拉取预处理Python 脚本手动执行工作流调度失败重试阈值固定阈值探索按土地类型、季节动态标定存储本地 tif对象存储 数据库告警终端打印面积消息通知 工单系统核查没有支持人工标记与反馈回滚不需要需要版本化模型和配置生产环境的另一个重点是配置外置化。影像目录、阈值、缓冲区大小、通知地址这些参数不要写死在代码里应该放到配置文件或环境变量中方便调整和审计。7.3 结果如何交付给业务人员遥感专业人员眼里的“NDVI 突降区”和巡护人员眼里的“需要去现场看的位置”是两回事。交付结果时要给出可直接使用的信息斑块边界坐标。前后两期影像的缩略图。变化幅度和面积。置信度说明。现场核查需要拍摄的信息点。这样巡护人员才能带着任务出去而不是带着一堆栅格文件出去。8. 还可以继续深挖的方向8.1 时间序列算法两期差值只是入门。保护区地表扰动检测更适合使用时间序列方法例如对每个像元的 NDVI 历史数据计算标准差标记偏离历史基线的时段。使用 BFAST 检测时间序列中的突变断点。使用 LandTrendr 检测植被覆盖的快速损失和长期趋势。这些方法更适合区分“长期退化”和“突然破坏”也能减少季节干扰。8.2 更高分辨率数据Sentinel-2 的 10 米分辨率可以发现较大的机械作业区但对小型违法开挖、车辆路径或局部清除仍然不够。如果有预算可以在疑似区域启用更高分辨率影像或无人机航拍进行二次确认。用“卫星遥感初筛 无人机精查 巡护人员现场核实”的组合比单靠任何一种方式都可靠。8.3 与巡护、无人机、工单系统联动更成熟的保护区监测系统会把遥感告警和巡护任务打通遥感发现疑似扰动后自动生成巡护工单巡护人员在移动端接收坐标到现场拍照并回传系统根据回传结果优化检测阈值。这样遥感模型会越用越准而不是永远停留在实验室阶段。8.4 给新手的练习路径如果刚开始接触遥感变化检测建议按下面的顺序练习找一景 Sentinel-2 L2A 影像计算单景 NDVI。在影像上手动模拟一个“扰动区域”再对比相邻时相。用两期差值定位变化区域。加入 SCL 云掩膜观察误报减少的效果。收集连续 5 期影像做时间序列跳变检测。把检测结果导出成矢量导入 QGIS 和业务系统测试。回看整个流程保护区地表扰动检测真正的门槛不是算法而是数据、坐标、云掩膜和阈值标定这些反复出现的细节。任何一张扰动图都应当被当成“待核查线索”而不是最终结论。下一步可以先找一个熟悉的小区域把 10 景历史影像下载下来跑通从 NDVI 到差分再到候选点导出的闭环再逐步把时间序列、空间过滤和生产调度加进去。