NDRI裸岩指数详解：从波段原理到阈值分割与裸露率计算

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的实现难点不在算法，而在波段匹配。不同的卫星传感器，波段编号和波长范围不一样，需要先确认你的影像里哪个波段是近红外、哪个是短波红外。

我在常用数据源里的对应关系是：

还有一个容易踩的坑：文献中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往极端负值拉，部分像元还可能越过高阈值变成假岩石。

所以我的标准预处理流程是：

1. 选取云量低于10%的影像，优先选择植被生长季和非雨季时相；
2. 用Scene Classification Layer波段（Sentinel-2自带）或CFMask算法结果生成云和云影掩膜；
3. 将掩膜区域在后续分析中设为无效值，不参与统计；
4. 如果研究区地形起伏超过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_path=None): 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(divide="ignore", invalid="ignore"): ndri = (swir - nir) / (swir + nir) ndri = np.where(np.isinf(ndri), np.nan, ndri) profile.update(dtype="float32", count=1) 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.08，20个非岩石样点（草地、裸土、建筑）的最大值是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最大值合成一张“年度最强裸露状态”图。云影和季节性植被遮蔽都会被压低，剩下来的高值区基本就是该年真正见底的位置。这个方法不复杂，却经常被忽略，实测做出来比单景影像稳定得多。