☰
ENVI遥感水质反演实战:从模型原理到叶绿素a与悬浮物制图
2026/10/4 9:06:09 网站建设 项目流程

经常有同行问我,说看到别人用ENVI做水质反演,一键就能出叶绿素a或者悬浮物的分布图,自己照着操作却总是不对劲,要么反演结果出现成片负值,要么模型精度低到不敢用。其实遥感水质反演这件事,听起来高大上,剥开来看就是“找关系”——找水体反射率和水质参数实测值之间的统计关系,再把这个关系推广到整个影像。ENVI只是帮你把中间这些波段运算、回归拟合、制图输出的流程串起来,真正决定成败的,是你对反演模型本身的理解深度。

这篇博文我会老老实实把水质反演模型那点事拆开讲清楚,从模型原理讲到ENVI实操,再讲到我用这些模型踩过的坑。不管你手里是Landsat 8/9还是Sentinel-2,不论你想反演叶绿素a、悬浮物还是透明度,核心思路都是相通的。哪怕你刚接触遥感,只要跟着把每个环节的逻辑理顺,也能在ENVI里跑通一条完整的水质反演流程。

1. 遥感水质反演的基本逻辑:从反射率到水质参数

1.1 水质参数为什么能用遥感来反演

要理解反演,先得知道遥感卫星拍到的是什么。卫星传感器记录的是地物反射太阳辐射的能量,经过量化后变成DN值,再转成大气层顶反射率,经过大气校正后得到地表反射率。水面反射率里面,包含了水体自身、水中悬浮物、浮游植物、溶解性有机物的散射和吸收信息。

水中的叶绿素a会吸收蓝光和红光,在近红外波段反射率较低;悬浮物浓度高时,水体在可见光到近红外的反射率整体抬升;透明度高的水体,蓝绿波段透过率高。这些光谱特征就是反演的依据。说白了,每个水质参数都有对应的光谱指纹,反演模型就是去量化“反射率变化了多少,对应水质参数变了多少”。

还有个概念必须搞清楚,我们反演的不是水质参数本身,而是通过遥感反射率建立统计关系来估算。这里面有物理基础,但更多是经验统计。所以模型的适用范围往往受限于建模样点的水体类型、季节、卫星传感器,离开这些前提去用,结果就会跑偏。

1.2 三大类反演模型:经验、半经验半分析、分析模型

遥感水质反演模型大致分三类。

第一类是纯经验模型(Empirical)。直接把实测水质参数和遥感波段反射率做回归,比如叶绿素a = a × Band3 + b。这种模型简单粗暴,不需要了解辐射传输机理,只要统计关系显著就能用。问题是样本代表性决定一切,换个湖或者换个季节,模型可能就失效了。

第二类是半经验半分析模型(Semi-empirical / Semi-analytical)。它结合了生物光学模型,从辐射传输理论出发,推导出水色参数与固有光学量之间的关系,再选择对目标参数敏感的波段比值或指数,最后用实测数据拟合系数。比如叶绿素a常用近红外与红波段的比值,悬浮物常用红波段与近红外的组合。这类模型物理基础强一些,适用性也更好。

第三类是分析模型(Analytical),也就是基于辐射传输方程,通过反演吸收系数和后向散射系数来估算水质参数。这类模型严格但需要大量实测光谱数据支撑,运算复杂,在工程应用中反而不如前两类普及。

在实际做项目中,绝大多数人用的还是第二类。因为它在ENVI里实现方便,模型结构有物理依据,精度又可以被实测数据约束。我后续讲的内容也以半经验半分析模型为主。

1.3 反演模型构建的核心流程概览

不管用什么模型,整个反演流程都是固定的几步。

第一步,数据准备。获取研究区的遥感影像,完成辐射定标和大气校正,得到地表反射率。这一步决定反射率数据是否可靠。

第二步,样本获取。采集同步或准同步的地面实测水质数据(叶绿素a、悬浮物、透明度、浊度、总磷等),同时记录采样点的GPS坐标。这些点是建模的真值。

第三步,光谱匹配。在影像上提取采样点对应的波段反射率,将反射率与水质参数一一对应,形成建模数据集。

第四步,模型构建。通过相关性分析、散点图观察,选择敏感波段或波段组合,用回归分析拟合模型参数。

第五步,模型应用。用ENVI的Band Math或光谱指数工具,把回归方程写成波段运算表达式,应用到整幅影像上,得到水质参数空间分布图。

第六步,精度验证。留出一部分实测数据不参与建模,用来验证模型预测值和实测值的一致性,计算决定系数、均方根误差等。

这套流程看着简单,每一步都有坑。下面我按步骤展开讲,每一步都说明为什么这么做,以及在ENVI里怎么操作最稳。

2. 数据准备与预处理:反演成败的第一道坎

2.1 影像选择:Landsat和Sentinel-2怎么选

遥感水质反演对影像有两个基本要求:一是空间分辨率要能匹配水面采样点,二是波段设置要包含对水质参数敏感的波段。

Landsat 8/9的OLI传感器有30米空间分辨率,波段覆盖蓝、绿、红、近红外,对内陆水体、近岸海域的叶绿素a和悬浮物反演非常成熟,历史数据从2013年到现在连续可选,适合做长期趋势分析。Sentinel-2的MSI传感器有10米和20米分辨率,而且多了三个红边波段,对富营养化水体的叶绿素a反演特别有用,尤其在浑浊水体里,红边波段比传统近红外波段更敏感。

做反演之前,你得先想清楚研究尺度和目标参数。如果是大范围湖泊水质监测,Landsat完全可以胜任;如果是小水体、河流、或者需要精细识别水华分布,Sentinel-2更合适。个人建议能用Sentinel-2尽量用Sentinel-2,红边波段带来的优势太明显,而且10米分辨率对混合像元的抑制也更好。

2.2 辐射定标与大气校正

很多人拿到影像就直接提取反射率,这是大忌。原始影像的DN值不是反射率,里面叠加了大气散射、吸收和传感器响应的影响。必须经过辐射定标和大气校正。

在ENVI里,辐射定标用Radiometric Calibration工具,输入影像后选择校准类型为“Reflectance”,输出数据类型建议Float,这样保留精度。注意Landsat和Sentinel-2的定标参数有区别。Landsat的Mtl文件自带辐射定标系数,ENVI可以直接读取;Sentinel-2 L1C产品需要先做大气校正,在ENVI中可以用“Sen2Cor”插件处理(独立于ENVI运行),输出L2A地表反射率产品,或者用ENVI自带的QUAC和FLAASH做粗略校正。

大气校正方法选择上,我强烈建议优先用FLAASH或6S这类基于辐射传输模型的工具,精度高,有物理依据。QUAC精度稍低,但在缺少大气参数时可以作为备选。做水质反演时,不推荐跳过大气校正直接反演,因为大气路径辐射对蓝绿波段的干扰非常严重,会导致反演结果系统性偏高。

2.3 水面像元提取与掩膜

拿到地表反射率以后,第一步不是急着建模型,而是把陆地、植被、云和云影从影像里去掉。水下有底质或者水体很浅的地方也建议剔除,否则这些像元的光谱是水底反射和水的混合,反演出来是假的。

在ENVI里做掩膜有几个思路。

一是用NDWI归一化差异水体指数。在Band Math里输入表达式(green - nir) / (green + nir),其中green是绿波段,nir是近红外波段,然后用阈值(一般大于0.1-0.2)提取水体。适用于水体与陆地对比明显的内陆湖泊、水库。

二是用NDVI加阈值。对富含悬浮物的浑浊水体,NDWI容易漏提,可以结合NDVI(植被指数)剔除植被,再结合红外波段设置反射率阈值来分离水陆。

三是直接人工矢量化。如果你的研究区只是一个小水库,而且水面边界清楚,直接在影像上画ROI(感兴趣区)最干净。不要觉得手工画费劲,有时候最笨的方法反而最可靠,尤其是在城市小水体里,自动水体提取算法很容易把阴影和建筑物误判为水。

掩膜生成后,建议保存为mask文件,在后续所有波段计算、统计提取中都使用这个掩膜,保证像元范围一致。

2.4 波段反射率提取

在建模之前,需要把每个采样点对应像元的反射率提取出来。操作路径是:ROI Tool → 打开采样点矢量文件(或手工创建ROI)→ 在影像上生成ROI → 用“Spatial Statistics”或“Extract Values to Points”工具,把每个像元的各波段反射率导出成表格。

千万不要直接把采样点落在影像边界上。由于定位误差,采样点可能落在相邻像元里。稳妥做法是,对每个采样点取周围3×3像元的平均值作为光谱值,这样能有效消除空间配准误差的影响。ENVI里可以用“Buffer Zone”给采样点做缓冲区,再提取缓冲区内像元均值。

提取出的数据整理成一张表,结构一般是:站点ID、叶绿素a实测值、悬浮物实测值、B1_B2_B3...各波段反射率。后面建模和回归分析都用这张表。

3. 在ENVI中构建水质反演模型:实操手记

3.1 建模数据准备:实测水质数据与影像反射率匹配

建模数据质量直接决定模型上限。要特别留意匹配的时间窗口。水体水质变化快,最好是卫星过境当天采样。如果做不到当天,最多放宽到前后1-2天,而且要选期间没有明显降雨、风浪的水体。否则你拿到的实测值是今天的,影像反射率是前天的,两者根本对不上。

另外,采样点要覆盖足够宽的浓度范围。比如你做叶绿素a反演,采样点的实测浓度要能从低到高分布开,最好包含贫营养、中营养、富营养等不同状态的水体。如果所有样本都集中在浓度差不多的范围里,回归模型虽然拟合得不错,但外推到高浓度或低浓度区域时就完全失效。

在ENVI里提取光谱值之前,建议把采样点矢量文件通过File → Open打开,然后叠加在影像上检查每个点是否落在有效水体范围内。位于岸边、芦苇丛、船坞附近的点要手动剔除掉。

3.2 波段组合与敏感波段筛选

建模型之前,先做相关性分析,找出哪些波段或波段组合与水质参数的关系最密切。这一步在Excel、SPSS或者Python里做都行,我是用Excel加Python结合来做的。

将实测叶绿素a浓度作为因变量,各波段反射率作为自变量,逐一计算Pearson相关系数。一般叶绿素a在蓝绿波段呈负相关,在近红外波段可能呈正相关。然后计算波段比值、差值和归一化指数,比如B5 / B4、(B5 - B4) / (B5 + B4),再与叶绿素a计算相关性,选出相关系数绝对值最大的组合。

这里有个经验:相比单一波段,波段比值类指数往往能显著提高相关性。原因是比值运算可以在一定程度上消除大气校正残留误差和水面粗糙度的影响,相当于做了一次归一化。这也是半经验模型中为什么大量使用比值指数的原因。

筛选敏感波段时,还要注意波段之间是否存在严重共线性。如果两个波段相关系数已经达到0.95以上,同时放进多元回归模型会导致系数不稳定。建议先看相关矩阵,再决定选哪几个波段。

3.3 回归模型构建:单波段、比值、多元线性

模型构建在统计软件里做,但ENVI也能做。我说两种路径。

路径一,用ENVI自带回归功能。在工具箱中找“Band Math”虽然能算,但统计回归更适合用ENVI的“Earth Analytics”或第三方扩展。比较传统的方式是导出Excel后在Excel里做线性回归。

路径二,也是我推荐的做法:在Excel或Python里拟合回归方程,然后把方程搬到ENVI的Band Math里应用。这样统计检验更直观,自由度也更高。

模型形式一般从简单到复杂逐个试。

单波段线性模型:y = a * x + b,x为单个敏感波段反射率。适合相关性已经很高的情况,但反演精度波动大。

比值线性模型:y = a * (B5 / B4) + b。这个常用在叶绿素a反演,能有效抑制环境干扰。

对数模型:y = a * ln(x) + b或y = a * ln(B5 / B4) + b。悬浮物浓度高时,反射率与浓度呈非线性关系,对数变换后线性化效果更好。

多元线性模型:y = a1 * x1 + a2 * x2 + a3 * x3 + b,多个波段或指数同时进入。适合单波段解释力不足的水体,但要非常注意过拟合。样本量小于30时,不建议使用超过2个自变量。

拟合完看几个关键指标:决定系数R²,越接近1越好,一般要求大于0.6才有实用价值;p值小于0.05说明关系显著;均方根误差RMSE衡量反演误差。还要看残差分布,如果残差随实测值有规律变化,说明模型形式可能错了,需要换结构。

3.4 模型精度的评价指标:R²、RMSE、MAPE

建模时还有一个很常见的错误:把用于建模的样本拿去检验模型精度。这样做算出来的精度虚高,因为模型已经把样本的“脾气”都记住了。正确做法是把样本分成两份,比如70%用来建模,30%用来验证,或者用交叉验证。

精度评价指标要分清建模精度和验证精度。我一般同时报告三个数。

R²也就是决定系数,反映模型对实测值方差的解释比例。验证集R²高的模型才值得用。

RMSE均方根误差,单位与实测参数相同,直接反映误差大小。比如叶绿素a的RMSE是5 μg/L,说明平均误差在5个单位左右。RMSE对异常值敏感,如果数据里有极端高值,RMSE会被拉得很大。

MAPE平均绝对百分比误差,用百分比表示误差占比。适合比较不同水质参数之间的模型精度,但叶绿素a在低浓度时MAPE会非常大,因为低浓度除出来的百分比数值很高,这点要谨慎解读。

最后,一定要把模型预测值和实测值的散点图做出来,看1:1线附近的数据分布。如果点云明显偏离1:1线,说明模型存在系统性偏差,可能需要在回归模型中增加截距项,或者改用非线性模型。

4. 经典反演模型案例:叶绿素a和悬浮物浓度

4.1 叶绿素a反演:为什么多用近红外与红边

叶绿素a是反映水体富营养化的核心指标。它有两个光谱特征:在蓝紫波段(440-470nm)和红波段(670nm附近)有强吸收峰,在近红外波段(700nm附近)有荧光峰。当叶绿素浓度升高时,红光波段的反射率下降,近红外波段的反射率上升,所以比值模型特别有效。

对于Landsat 8 OLI传感器,常用波段组合是Band5 / Band4(近红外/红)。对于Sentinel-2,我更推荐使用红边波段组合,比如Band5 / Band4或Band6 / Band5甚至Band7 / Band5。红边波段位于叶绿素吸收与散射的过渡区,对叶绿素浓度的敏感性比宽波段近红外更高。

一个典型模型形式是:

Chl-a = a * (B5 / B4) + b

其中B5和B4分别对应影像的近红外和红波段。我在项目里拟合过一个Landsat 8的模型,R²能达到0.75左右,RMSE约4.2 μg/L。换到Sentinel-2用红边比值后,R²直接提升到0.83,说明红边波段的优势很明显。

4.2 悬浮物浓度反演:波段比值与对数模型

悬浮物在光谱上的表现和叶绿素相反,它主要增强水体反射率,波长越长散射越明显。因此红波段和近红外波段的反射率与悬浮物浓度呈正相关。在浑浊水体里,近红外波段的反射率对悬浮物非常敏感,但在清水里反射率很低,容易受噪声干扰。

常用反演模型有两种。

线性比值模型:SS = a * (B4 / B3) + b,红波段与绿波段比值。这个模型对轻微浑浊水体较稳,但在高浓度时会饱和。

对数模型:SS = a * ln(B4) + b或SS = a * ln(B4) + b * ln(B5) + c。因为悬浮物浓度与反射率之间是指数型关系,浓度越高反射率增幅越小,取对数后关系更接近线性。

我在珠江口某项目里用过Sentinel-2的模型:SS = 452.6 * ln(B4) - 93.7,验证集R²为0.79,RMSE约12.6 mg/L。注意这个公式只能用于相近水体和同样传感器,换数据必须重新拟合。

4.3 在ENVI中用Band Math实现模型应用

假设你已经通过回归分析得到叶绿素a的模型为:

Chl-a = 36.5 * (B5 / B4) + 8.2

注意这里的B5、B4是对应波段反射率,也就是你在ENVI里做完大气校正后的反射率文件波段。以Landsat 8为例,Band Math表达式中应该写作:

36.5 * (float(b5) / float(b4)) + 8.2

操作步骤:

在ENVI工具箱搜索“Band Math”,打开对话框,在“Enter an expression”里输入上面的公式。

点击“Add to List”,然后选择波段映射。系统会弹出变量定义界面,你要把b5对应到影像的第五波段,b4对应到第四波段。一定不要选错,很多新手在这里把b5和b4对应反了,结果出来的分布图完全相反。

输出文件选择Float Single类型,并定义好输出路径。如果想避免出现边缘突变,可以勾选“Output Log”或使用掩膜文件限定计算范围。

方程组写完后,点击OK执行。结果影像中每个像元的灰度值就是预测叶绿素浓度。别急着出图,先检查统计值,看最小值和最大值是否合理。如果出现大量负值,通常是对应波段反射率接近零或模型拟合不佳造成的,这个我在第5部分会细说。

如果你用的是Sentinel-2影像,波段顺序可能和Landsat不同,建议先查一下影像元数据,或者在ENVI的Layer Manager里确认波段波长。比如Sentinel-2的Band4是红波段,Band5是植被红边波段,Band8是近红外,不同产品的波段编号虽然一样但波长并不同,一定要用波长来核对。

5. 实战中的常见问题与排查技巧

5.1 反演结果出现负值,怎么处理

反演结果出现负值几乎是每个人都会遇到的。产生负值的场景主要有三种。

第一种,影像反射率本身有负值。大气校正过度校正时,暗像元的蓝绿波段反射率会被压成负值,代入模型后自然得到负浓度。处理办法:在Band Math的输出表达式中加一个条件判断,比如把小于0的设为0或设为无效值。表达式可以写成:

(b1 lt 0) * 0 + (b1 ge 0) * b1

第二种,模型外推导致负值。当水体反射率超出了建模样点的取值范围,回归方程计算出来的结果可能为负。这种情况要检查建模样本是否覆盖了全影像的反射率动态范围,如果影像里有特别清澈或特别浑浊的区域超出样本范围,最好对反演结果做截断处理,比如把小于0的值设为0,大于实测最大值150%的值设为无效。

第三种,波段运算时数据类型问题。如果你用的是整数型影像,除法和减法可能导致截断误差。务必在Band Math里用float()函数把波段转成浮点型。

我个人的处理思路是:先检查负值像元的空间分布和比例。如果只是零星出现在水陆交界处,直接掩膜掉就行;如果成片出现在开阔水面,说明大气校正出了问题,回头重新做大气校正比在结果上修补更靠谱。

5.2 模型精度很差,可能错在哪

精度差一般不是模型公式的锅,而是数据链路上出了问题。我排查的顺序是固定的。

先看实测数据和影像时间是否匹配。采样时间差超过3天,精度很难保证。再检查大气校正是否有效。我有个验证方法:找一块均匀的深水体,查看蓝绿波段反射率是否处于合理范围(一般在0.02-0.06之间),如果反射率在0.1以上,说明大气校正偏得离谱。

然后看采样点的空间代表性。在GPS定位误差大的情况下,比如手持GPS定位精度5-10米,对于Landsat 30米像元还能凑合,但对于Sentinel-2的10米像元就很容易偏移。所以我前面反复强调,一定要取3×3像元均值。

最后再看模型本身。样本量是不是太少了?浓度范围够不够宽?有没有异常值在作怪?画一个箱线图检查实测数据,离群点该删就要删,但删除要有依据,比如该点明显受局部排污影响,或者采样记录里备注过异常情况。

5.3 大气校正选择哪种方法更稳

这里专门说一下大气校正的方法选择,因为这是最多人纠结的点。

ENVI中有FLAASH、QUAC、6S等选项。我按稳定性排序:6S和FLAASH属于物理模型,精度最高,但FLAASH需要输入大气模型参数(水汽、气溶胶类型、能见度等),参数给不对反而会出错;QUAC是快速近似方法,不需要输入大量参数,直接用影像自身信息估算,处理速度快,但精度在蓝波段和浑浊大气条件下会明显下降。

做水质反演,我的建议是优先使用ESA发布的Sen2Cor处理Sentinel-2 L2A产品(可以直接下载L2A级别产品),Landsat则可以用USGS官方提供的LEDAPS或LaSRC处理过的表面反射率产品。直接下载这些官方表面反射率产品,比自己在ENVI里用FLAASH校正更稳。

如果你非要自己处理,记得至少做两步:第一步,检查影像中有无明显薄云和雾霾;第二步,对校正结果做一个简单的波谱曲线检查。方法是在影像中找一块清澈深水区,提取反射率光谱曲线,如果水体在红波段和近红外波段的反射率很低(低于0.03),说明校正结果基本合理。

5.4 样本点不多,怎么提高模型稳健性

很多项目受条件限制,现场采样的样本量只有十几个甚至几个。这种情况还想建一个能用的模型,我有几个实操技巧。

第一,合并多期影像的样本。把同一传感器在相近时段内采集的影像和对应实测数据放到一起建模,相当于扩充了样本量。前提是各期影像的大气校正方法一致,水体状态没有发生剧烈变化(比如没有发生水华暴发)。

第二,采用留一交叉验证。样本量少时,传统划分训练集和验证集会浪费样本。留一法每次用全部样本中N-1个建模,剩下1个验证,循环N次,计算平均精度。虽然过程麻烦点,但对样本利用效率高。

第三,尽量用简单模型。样本量少于20时,单个波段或单波段比值模型远比多元线性模型稳健。不要追求R²高而引入太多自变量,拟合得好不叫好,验证得好才叫好。

第四,用已知文献模型先做一个粗略反演,然后用稀疏实测点做线性校正。举个例子,先引用别人的模型计算初始浓度分布,然后建立“实测值 = a × 模型预测值 + b”的校正方程,把实测点代入求a和b。这种方法在实测点少时很实用,能借用前人的先验知识。

6. 反演结果的可视化出图与后续扩展

6.1 密度分割与专题图制作

反演得到的水质参数分布图,灰度图是没法直接用的,还要做可视化。最常用的是密度分割,就是按照浓度阈值把连续灰度值切成若干等级,每级赋予不同颜色。

在ENVI里可以用“Density Slice”工具,设置阈值范围。比如叶绿素a可以分成小于2、2-5、5-10、10-20、大于20 μg/L五级。阈值要根据研究区水体富营养化评价标准来定,不要随便切。

切完以后,在ENVI中通过“Apply Color”给不同区间赋色,再叠加到研究区底图上,加上图例、比例尺、指北针,导出为TIFF或JPG。如果需要做论文插图,建议导出带坐标系的GeoTIFF,然后在ArcMap或QGIS里添加图例和格网。

6.2 时间序列分析与动态监测

反演模型一旦建立,可以应用到多期影像上,分析水质参数的时空变化趋势。我做过一个水库项目,把过去5年LandSat影像全部做了一遍大气校正和反演,计算每期水面平均叶绿素浓度,再画时间变化曲线,成功识别出夏季水华的高发时段。

做时间序列时有几个细节要注意:不同时期影像的大气条件不同,直接比较反射率可能有偏差。因此最好选择同季节、同传感器的影像,或者对每期影像做同样的相对辐射归一化。另外,模型如果是从某一期影像建出来的,直接套到其他影像上时最好用当地实测数据做一次校正,否则可能会存在系统偏移。

6.3 结合机器学习模型进一步优化

近年来,机器学习在水质反演里用得很火。比起传统回归,随机森林、支持向量机、神经网络能自动捕捉波段与非线性的关系,适合复杂水体。

在ENVI中直接跑机器学习不太方便。我的做法是:先用ENVI完成预处理和光谱提取,把样本数据导出成CSV,然后在Python里用Scikit-learn训练模型,最后把最优模型方程再转回ENVI的Band Math,或者用ENVI的“Save Model”和“Predict”功能备份。当然,如果你熟悉Python,直接全程用Python做遥感影像处理和模型训练也是可行的。

需要提醒的是,机器学习模型更容易过拟合,样本量小的项目慎用。如果样本量只有二三十个,推荐用简单随机森林加交叉验证,但结果要用传统回归模型做对照,避免盲目堆模型复杂度。

最后再分享一个我自己的习惯:每次做完反演,一定会把建模样本的实测值与预测值散点图保存下来,连同模型方程、精度指标、影像预处理参数一起归档。很多项目过了几个月回头要补数据、改模型时,这些记录能省下大量重复工作。水质反演这个活儿,头发掉得多的往往不是建模那几步,而是前期数据没整明白、后期文档找不到。你把过程做得规矩一点,后面会顺很多。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询