1. 裸露岩体监测实地量测的痛点和遥感替代方案1.1 常规实地调查为什么撑不住大面积评估我前几年参与过某矿区的生态修复效果评估项目方最关心一个问题山体上到底还有多少裸岩按现在的治理速度几年能复绿。当时团队先按老办法出方案——布设样线、打样方、带手持GPS沿山脊走每50米记录一次地表类型。真走起来才发现裸露岩区基本都在陡坡和冲沟两侧人爬上去需要绕路一天走完一两条样线就到极限。一个面积一百多平方公里的评估区样线覆盖率连百分之一都不到最后统计出来的裸露率误差大到没法用。无人机航拍能解决一部分问题但大范围航测需要申请空域、拼图、人工勾绘周期按周算。如果只需要一个宏观的“裸露面积占比”还有一种更直接的思路用卫星多光谱影像计算归一化岩石指数NDRI把岩石提取问题变成一次波段运算。NDRI全称Normalized Difference Rock Index很多资料里也直接叫裸岩指数它的输出是每个像元的连续数值阈值化之后就能统计岩石像元占比也就是岩石裸露率。这个指标的优点在于不需要额外采集训练样本就能出图计算成本低而且用公开的遥感数据就能实现适合区域尺度的快速评估。1.2 多光谱指数为什么能代替人眼判断有人会问既然是“裸岩”我在影像上直接目视解译不就行了大范围影像上目视解译确实可行但效率低、标准不统一而且影像受云影、光照和大气影响不同时间片的色调差异很大靠人眼很难稳定地分辨岩石和裸土。多光谱指数的思路是把“颜色差异”变成数值差异不同地物在不同波段上的反射率曲线不同波段比值相当于把曲线形态的差异放大再压制掉光照强度的影响。NDRI能代替人眼本质上是因为岩石、植被、水体在近红外和短波红外两个波段上的反射特征正好相反。基岩和风化面在短波红外区反射率较高而在近红外区相对偏弱植被却是在近红外区反射很强、在短波红外区因为叶片水分吸收而大幅下降。两者一对比指数值天然就分开了。理解了这一点后面所有步骤——选波段、定阈值、排除误判——都有了判断依据。2. NDRI公式拆解近红外与短波红外的光谱博弈2.1 三种地物在两个波段上的反射特征NDRI的标准表达形式是NDRI (SWIR - NIR) / (SWIR NIR)其中NIR是近红外波段SWIR是短波红外波段。这个式子长得和归一化植被指数NDVI一模一样只是把NDVI里的红光波段换成了短波红外因此它本质上继承了比值指数的所有优点和所有坑。先看岩石。干燥岩石表面在短波红外大约1.5到1.8微米区间反射率普遍比较高这跟矿物成分有关硅酸盐矿物、碳酸盐类矿物在这个区间内没有特别强烈的吸收带整体反射能量比近红外区要高一些。植被则相反近红外是叶片细胞结构散射最强的区间反射率往往能到0.4以上而到了短波红外叶片水分含量高水分吸收让能量掉到0.1以下形成明显的“低槽”。至于水体两个波段都在0.05以下的低位更像一个暗目标。把三个特征放到同一对波段里来看岩石在SWIR上高、NIR上中等NDRI产出正值植被在NIR上很高、SWIR偏低NDRI产出负值水体两个波段都很低比例不稳定偶尔会出现噪声导致的正值。这就是后续误判的主要来源后面会专门讲怎么清理。2.2 比值公式的物理意义和它压制光照的特性为什么要用比值而不是直接用某个单波段的亮度值因为卫星传感器的辐射响应受太阳高度角、地形坡向、大气散射的共同影响。同一个岩性、同一条山坡阳面和阴面在单波段上的亮度可以差出一大截直接拿亮度阈值去分割阳面能识别阴面就被漏掉了。比值运算把两个波段的亮度放在同一个像元内比较分子分母同时受光照影响的部分会发生抵消。地貌学里管这叫“归一化”——它不考虑绝对亮度只看两个波段反射率曲线的形态差异。NDRI值域在-1到1之间岩石类地物通常在0.1到0.5区间植被通常在-0.4到-0.05区间。这个区间差非常明显是NDRI能快速分类的核心物理基础。需要提醒的是比值指数对暗像元特别敏感。当SWIR和NIR两个波段的值都很低时比如云影、山体阴影覆盖的像元哪怕反射率只差0.01算出来的NDRI也可能出现很大的正值或负值属于典型的“低信噪比区域放大”。这也是为什么很多教程里直接拿原始影像计算NDRI会出现大量椒盐噪声看起来像模像样一验证全是虚警。2.3 不同传感器的波段对应与公式版本差异NDRI的实现难点不在算法而在波段匹配。不同的卫星传感器波段编号和波长范围不一样需要先确认你的影像里哪个波段是近红外、哪个是短波红外。我在常用数据源里的对应关系是数据源NIR波段SWIR波段建议分辨率注意事项Landsat 8/9 OLIB50.85-0.88微米B61.57-1.65微米30米尽量使用Collection 2 Level-2表面反射率产品Sentinel-2 MSIB80.78-0.90微米B111.56-1.65微米20米也可用B8A替代B8能减少大气吸收影响国产资源系列依据传感器载荷说明确定依据载荷说明确定视分辨率需额外注意辐射定标参数是否完备还有一个容易踩的坑文献中NDRI公式存在两种版本。有的文章采用(SWIR1 - NIR)/(SWIR1 NIR)跟上面一致但也有文献把可见光红波段替代NIR参与计算用来增强裸土和岩石的区分度。不同公式算出来的数值区间不一致直接拿别人文献里的阈值套用往往会错。我的建议是在同一项目内老老实实固定一个公式版本从数据准备到精度验证全程保持一致不要混用。3. 数据准备到波段运算实操前的各项准备3.1 数据源与波段对应关系计算NDRI的第一步是拿到正确的波段栅格。我的经验是优先选择官方发布的Level-2表面反射率产品也就是已经做完大气校正的数据。Landsat 8/9的Level-2产品直接在公开数据平台上按“Surface Reflectance”筛选即可Sentinel-2的L2A产品是官方默认提供的下载时已经自带大气校正结果不用额外再跑一遍辐射传输模型。如果只能拿到Level-1原始数据那大气校正这一步就跑不掉。商业软件里用FLASSH、QUAC都可以开源方案可以选Py6S或者用]如果真的只做NDRI而不做其他定量反演相对简单的简化校正法也能接受但要清楚波段比值虽然能抵消一部分大气程辐射的影响却不能被完全抵消尤其是短波红外通道对水汽吸收比较敏感大气校正不到位的影像NDRI直方图整体会偏移阈值也会跟着飘。3.2 大气校正与云掩膜必须做的两步在我做过的几次NDRI计算里最影响结果稳定性的不是公式本身而是云和云影。云体在短波红外附近反射率很高会把NDRI顶到接近1的异常值云影则正好相反把NDRI往极端负值拉部分像元还可能越过高阈值变成假岩石。所以我的标准预处理流程是选取云量低于10%的影像优先选择植被生长季和非雨季时相用Scene Classification Layer波段Sentinel-2自带或CFMask算法结果生成云和云影掩膜将掩膜区域在后续分析中设为无效值不参与统计如果研究区地形起伏超过30度额外用DEM生成山体阴影图对山影区域做单独标注方便后面对比验证。很多教程跳过了云掩膜这一步直接拿全图算NDRI再统计裸露率。这样统计出来的裸露率在干旱区看起来差别不大到了多雨山区云影区域的假岩石像元能把整体裸岩占比抬高好几个百分点误差直接超出可接受范围。3.3 用QGIS和Python分别算出NDRI计算这一步没有什么玄学关键是控制波段的像元值范围和空间对齐。如果只做单景影像打开QGIS的栅格计算器比如Landsat 8数据直接写(B6 - B5) / (B6 B5)输出为Float32类型就会生成NDRI栅格。注意一点如果下载的是整数型DN值产品比如某些Level-1数据直接做缩放需要先乘以缩放因子转成浮点型反射率再算指数否则分母和分子的数值接近指数会退化成一堆离散值。如果处理多期影像或者要批量做时间序列我建议直接用Python配合rasterio或者GDAL写脚本。下面这段是我常用的计算流程import rasterio import numpy as np def compute_ndri(nir_path, swir_path, out_path, mask_pathNone): with rasterio.open(nir_path) as nir_src: nir nir_src.read(1).astype(float32) profile nir_src.profile with rasterio.open(swir_path) as swir_src: swir swir_src.read(1).astype(float32) with np.errstate(divideignore, invalidignore): ndri (swir - nir) / (swir nir) ndri np.where(np.isinf(ndri), np.nan, ndri) profile.update(dtypefloat32, count1) with rasterio.open(out_path, w, **profile) as dst: dst.write(ndri, 1)脚本里有个小细节分母出现零或接近零的像元会产生inf或者nan不能让这些异常值继续参与后面的阈值和统计所以要在计算后立即用np.where处理掉。4. 阈值分割与裸露率统计从连续值到面积比例4.1 阈值怎么定直方图法和样地统计法NDRI计算出来后是连续的浮点栅格要得到岩石分布必须确定一个阈值把“是岩石”和“不是岩石”分开。阈值怎么定是最常见的问题。最省事的做法是看直方图。把NDRI栅格在QGIS里生成直方图正常情况下会出现两个峰一个在-0.2到0的区间对应植被和土壤一个在0.1到0.4区间对应岩石。取两个峰之间的谷底值作为初始阈值通常落在0.05到0.15之间。直方图法的问题是遇到裸地与岩石光谱重叠区域双峰不明显谷底很宽。这种时候用样地统计的思路更可靠在高分辨率影像或者无人机正射影像上随机挑选几十个岩石/非岩石样点提取对应位置的NDRI值分别统计两组的分布区间用两者的边界值做阈值。比如我在某矿区这样标定过20个岩石样点的NDRI最小值是0.0820个非岩石样点草地、裸土、建筑的最大值是0.16中间重叠了0.08最后取0.12作为阈值。这个值不是固定的换研究区、换影像时相之后要重新标定。4.2 后处理与分区域统计确定阈值后把NDRI栅格转成二值栅格高于阈值为1低于阈值为0。但直接这样用会有一个问题单独的像元噪声很多比如路边一两块碎石的像元、影像噪声点都会形成零星分布的小图斑。我做后处理的常用组合是3×3窗口中值滤波去除椒盐噪声设定最小图斑面积比如30米分辨率的影像小于3个像元的连片区域忽略不计如果研究区有云掩膜和水体掩膜先更新有效区域掩膜保证统计口径统一。分区域的裸露率统计也在这个阶段做。按行政区、子流域或者矿区边界建立网格矢量用分区统计工具求每个区域内岩石像元数占有效像元数的比例。输出的就是各分区岩石裸露率可以直接出专题图。4.3 裸露率计算的两种口径统计岩石裸露率时有两种口径建议根据应用场景选一种是像素口径裸露率 岩石像元数 / 全部有效像元数。这种口径适合拿来评估区域整体的裸岩覆盖情况比如生态修复前后的对比监测。另一种是有效面积口径裸露率 岩石面积 / 研究区有效面积。有效面积需要扣除云掩膜、水体面积避免把湖泊面积拉进分母。两者数值差异在山区往往很小但如果研究区内有大片稳定水体或者大面积云覆盖差异就明显了。严谨的做法是报告里同时写明扣除无效像元数和对应的有效面积。在代码层面用numpy统计是最直接的ndri_flat ndri[valid_mask] # 只保留有效像元的NDRI值 rock_pixels (ndri_flat threshold).sum() valid_pixels valid_mask.sum() bare_rock_ratio rock_pixels / valid_pixels print(f岩石裸露率: {bare_rock_ratio:.2%})5. 精度验证与误判清理水体、阴影和人工地表5.1 验证点位布设与混淆矩阵用NDRI算出来的裸露率如果只看数据不验证等于把结果建立在假设上。每一种岩性、每一种地表覆盖类型的光谱响应都有差异同一个研究区内花岗岩出露区和碎屑岩坡积物就可能对应完全不同的NDRI分布。我在一个项目中采用过的验证方案是在影像覆盖范围内按地貌类型分层随机布设验证点岩石区和非岩石区各布不少于30个点用无人机正射影像或者野外GPS对照判定每个点位的真实属性。最终统计一个简单的混淆矩阵分类结果实际为岩石实际为非岩石判别为岩石正确识别数误判数判别为非岩石漏判数正确排除数从矩阵中可以算出用户精度识别为岩石的像元里有多大比例确实是岩石和生产者精度真实岩石像元有多大比例被识别出来。通常用户精度更重要因为它反映的是你画到地图上的每一片“岩石”是否可靠。5.2 四类常见误判的排查思路在实际操作中NDRI提取结果最常见的误判来源是四个水体、阴影、裸土和人工地表。每类问题的处理逻辑都不一样。水体的问题在于部分浅水或浑浊水体在短波红外和近红外之间也可能出现微弱正差NDRI偶尔冒出正值。我的处理方式是用水体指数或者简单的绿色波段亮度过高的区域生成水掩膜在统计分析前置为无效值。阴影是最容易忽视的。地形阴影覆盖的岩石像元两个波段都被压低但比值不稳定部分出现正值。处理思路有两条一是利用DEM计算山体阴影将坡度较大且光照条件差的区域标记出来对比NDRI结果单独观察二是对两波段反射率绝对值设置下限只有两波段反射率都高于某一阈值时才认定岩石低于阈值的全部归为阴影或水体。裸土和干燥耕地跟岩石在光谱上确实有重叠而且这种重叠难以用单一个指数完全解决。目前比较实用的做法是引入地形辅助天然裸岩多见于坡度较陡、地表粗糙度大的区域裸土多分布在平缓地带用DEM坡度数据做一个约束能减少大部分裸土误判。人工地表比如水泥硬化地面、沥青路面反射形态跟岩石接近这一类没法靠NDRI单独排除。如果研究区范围内有建筑区或道路网络矢量最省事的方式是直接用矢量擦除如果研究区本身是自然山区这类误判占比通常很小人工审查即可。5.3 变化监测的一点扩展NDRI除了算静态裸露率做多时相对比价值更大。同一研究区、同一传感器、同一季节的两期NDRI相减能直接给出裸岩扩张或植被恢复的空间分布。需要注意两点一是两期影像的时间尽量对齐避免物候差异干扰二是阈值的选取如果跨期漂移不能直接比较像元级差值最好先分别按各自阈值提取二值图再比较变化像元比例。我做过多期对比后有一个小体会单期的岩石裸露率只反映“当前状态”真正对决策有用的是连续几年的裸露率曲线。第一年NDRI提取的裸岩区域第二年算出来下降了多少第三年是否反弹这才是生态修复评估需要看的东西。最后再分享一个实用技巧做NDRI时间序列时可以用每年多景影像的NDRI最大值合成一张“年度最强裸露状态”图。云影和季节性植被遮蔽都会被压低剩下来的高值区基本就是该年真正见底的位置。这个方法不复杂却经常被忽略实测做出来比单景影像稳定得多。