基于HY-1C CZI影像光谱指数重构数据MNF变换的红树林提取

梁超，刘利，刘建强，邹斌，邹亚荣，崔松雪

(1. 自然资源部 国家卫星海洋应用中心，北京 100081；2. 自然资源部 空间海洋遥感与应用研究重点实验室，北京 100081；3. 中国科学院 空天信息创新研究院，北京 100094)

1 引言

红树林是指生长在热带、亚热带近岸潮间带上部滩涂浅滩，以红树植物为主体的常绿灌木或乔木组成的潮滩湿地木本生物群落，是陆地向海洋过渡的特殊生态系统。红树林具有防风消浪、促淤保滩、净化海水、保持生物多样性等重要功能。近几十年来，由于围海造地、围海养殖、砍伐等人为不合理开发活动等因素影响，使得红树林面积逐年锐减。红树林生态环境保护工作对长期监测红树林的分布范围、生物量、健康水平等提出了要求，其中红树林分布范围的监测是基础。利用遥感手段可以大范围、快速、重复监测红树林时空分布信息，这对于研究保护红树林生态系统具有重要意义。

国内外红树林遥感监测方面开展了大量研究，所采用的数据源通常为NOAA-AVHRR、Landsat、SPOT、Quickbird、ALOS/PALSAR等卫星以及诸如HyMap等航空高光谱遥感影像，主要方法包括人工解译、监督分类、植被指数及纹理特征、面向对象分类等[1-14]。HY-1C是我国第一代海洋水色卫星HY-1A/B的后继星，于2018年9月7日在太原卫星发射中心发射成功，是我国首颗业务化运行、全球开机、连续观测的海洋水色卫星，载荷配置和性能指标、探测能力等均有明显改进，可获取重访周期更短、观测范围更大、精度更高的水色遥感产品。HY-1C搭载的海岸带成像仪（Coastal Zone Imager，CZI）用于获取海陆交互区域实时数据，监测近海、海岛、海岸带等信息，其幅宽、分辨率相比上一代载荷分别提升近2倍和5倍。本文拟以广西山口红树林生态自然保护区为研究区，综合利用红树林与一般陆地植被光谱特征及空间分布特征，选用适当的光谱指数重构原始数据，经最小噪声分离变换（Minimum Noise Fraction Rotation，MNF）后，建立决策树提取红树林，据此探讨HY-1C卫星CZI数据在红树林监测工作中的适用性。

2 数据与方法

2.1 研究区概况

本文以广西山口红树林生态自然保护区为研究区（图1），该保护区由位于广西合浦县东南部沙田半岛东、西两侧的海域、陆域及全部滩涂组成，地处亚热带，是我国第二个国家级的红树林自然保护区，分布着发育良好、结构典型、连片较大、保存较完整的天然红树林，有红海榄、木榄、秋茄、桐花树等15种红树林植物。山口红树林国家级自然保护区具典型的大陆红树林海岸生态系统特征，红树林中还栖息着多种海洋生物和鸟类，具有重要的科学价值。

图1 研究区位置与HY-1C卫星CZI原始影像Fig. 1 Study area location and the CZI original image of HY-1C satellite

2.2 数据源

采用2018年10月31日11时25分获取的HY-1C卫星CZI传感器L1B级数据，CZI主要技术指标见表1。

表1 HY-1C卫星CZI传感器主要技术指标Table 1 The main technical indicators of HY-1C satellite CZI sensor

2.3 研究方法

本文采用波段重构与光谱变换相结合的方法，以增强CZI影像上红树林与一般陆地植被的特征差异，提高二者可分性。图2为研究方法流程图，首先对CZI数据预处理，然后基于红树林与一般陆地植被光谱特征构建多种光谱指数，通过相关性分析筛选最佳指数组合重构原始数据，继而采用MNF进一步增强重构数据中红树林与一般陆地植被的光谱及空间分布特征差异，最后建立决策树提取红树林。

图2 本文研究方法流程图Fig. 2 The research flow chart

2.3.1 数据预处理

数据预处理包括几何校正和辐射定标，CZI数据采用HDF5格式，几何校正利用数据自带的RPC参数。L1B数据像元值为辐射亮度，据公式（1）转换为表观反射率[15]：

式中，ρ为大气层顶表观反射率（无量纲），Lλ为大气层顶进入卫星传感器的光谱辐射亮度，单位为W/(m2·sr·μm)，D为日地距离 (天文单位)，ESUNλ为大气层顶的平均太阳光谱辐照度，单位为W/(m2·μm)，θ为太阳天顶角，相关参数可由HDF5数据集属性参量中查找。

2.3.2 光谱指数

红树林与一般陆地植被光谱特征相似，在CZI原始影像难以区分。为突出植被、水体等信息并增强红树林与一般陆地植被的差异，本文拟选用适当的光谱指数重构CZI原始数据。综合考虑大气、红树林分布环境等因素的影响，本文选用了常见的8种植被指数[16-19]和1种水体指数，为进一步突出CZI影像上植被光谱特征，构造了3个波段比指数和1个光谱斜率比指数，其中波段比指数用来表达目标物在不同波段之间的归一化差值比，可见光光谱斜率比指数用来表达不同地物在可见光波段的光谱曲线形态差异。表2列出了本文所采用的13种光谱指数。

2.3.3 最小噪声分离变换

MNF最早由Green等[20]提出，是在主成分变换基础上提出的一种线性变换算法。主成分变换各分量按照方差降序排列，影像主要信息集中在前几个分量，但不能保证影像质量（信噪比）也按照降序排列，MNF用以改进主成分变换在噪声消除及影像增强上的不足。Boardman和Kruse[21]证明了MNF等效于连续两次主成分变换，设多光谱数据由加性噪声和信号两部分组成，首先对影像噪声协方差矩阵进行变换，使变换后噪声协方差矩阵为单位矩阵且波段间不相关，然后对变换后的数据集作标准主成分变换[22]。

第一步，对多光谱影像X的噪声协方差矩阵进行对角化：

式中，Cxn为原始影像噪声协方差矩阵，Λ为其特征值组成的对角阵，E为由其特征向量组成的正交矩阵。

上式前乘以 (Λ-0.5)T，后乘以Λ-0.5，则

式中，I为单位矩阵，令F=EΛ-0.5，则上式转换为

原始影像经Y=FTX变换后的协方差矩阵为

表2 本文使用的植被光谱指数列表Table 2 List of vegetation spectral indices used in this paper

式中，Cx为原始影像协方差矩阵。

第二步，对Y应用标准主成分变换：

式中，ΛY是CY特征值降序排列的对角阵，B是其特征向量组成的正交矩阵。

综上，MNF矩阵为

多光谱数据经过MNF后各分量按照信噪比降序排列，影像有效信息集中在前几个分量，噪声主要存在于后面的分量中，实现了光谱增强的同时分离数据中的噪声，改善多光谱处理结果。

MNF最重要的一步是估计原始影像噪声协方差，其方法一般借助影像的空间特征，据Canty[23]：

式中，x=(i,j)T表示影像X的一个像元坐标，h=(h1,h2)T表示坐标的偏移量，即噪声协方差可以用原始影像与偏移影像差异的协方差来估计。

3 结果分析

3.1 CZI影像上植被光谱特征分析

通过建立感兴趣区，提取红树林与一般陆地植被样本的表观反射率光谱曲线（图3），可见在CZI影像上二者均表现出典型的绿色植被光谱特征，即绿波段强反射、红波段强吸收以及近红外波段的高反射率，蓝波段由于大气散射的影响亦表现出较强的反射，而近红外波段上二者光谱曲线均较为发散，这是由于中低分辨率影像上混合像元效应的影响。红树林与一般陆地植被光谱曲线具有较高的相似性，使得二者很难在4波段CZI数据中被区分。图4为CZI影像上红树林与一般陆地植被表观反射率分布直方图，可见二者在各波段具有大致相近的直方图形态和峰值分布区间，虽然在第4波段，红树林与一般陆地植被直方图峰值位置略有差异，但如前述，红树林与陆地植被在该波段光谱曲线均较为发散，直方图上分布范围均较宽，从而存在较大的重叠区。综合考虑光谱曲线及直方图可知，单纯依靠CZI表观反射率数据很难区分红树林与一般陆地植被。

图3 CZI影像上红树林（a）与一般陆地植被（b）光谱曲线Fig. 3 Spectral curves of mangrove (a) and general terrestrial vegetation (b) on CZI image

图4 红树林（a）与一般陆地植被（b）表观反射率分布直方图Fig. 4 Histograms of apparent reflectance distribution of mangroves (a) and general terrestrial vegetation (b)

3.2 基于光谱指数重构多光谱数据

红树林与一般陆地植被在CZI原始光谱空间中较高的相似性对红树林检测带来了难度，传统的监督/非监督分类等方法提取红树林效果均不甚理想。本文通过构建对植被和水体信息敏感的光谱指数组合代替原始波段数据，以增强影像中典型植被的光谱信息及其差异性。

本文以植被遥感中最常用的归一化差异植被指数NDVI为参考基准，开展了NDVI与其他12种光谱指数的相关性分析，表3为各指数之间的相关系数。首先剔除与NDVI相关系数极低的光谱指数，如CBRI，因为极低的相关系数意味着该类指数对目标植被信息不敏感；其次，剔除与NDVI具有较高的正相关系数的光谱指数，如 SIPI、CNBI、MSAVI、EVI等，该类指数在表征植被信息方面与NDVI作用相当，从避免冗余考虑予以剔除；最后，选择与NDVI具有中、低相关性及最大负相关性的光谱指数参与重构，如 RI、NDWI、NDGI、CBGI、ARVI和 CVSSR 等，该类指数既能从不同程度反映植被信息，又避免了信息冗余，确保重构后的新光谱空间既能最大化地增强植被及水体信息，又能最大程度减少波段之间的相关性。据此原则并经多次试验，本文最终选取NDVI、NDWI、ARVI和CVSSR为最佳指数组合，并以之依次代替CZI 4个原始波段组成新的多光谱数据。

表3 NDVI与其他光谱指数间的相关系数Table 3 Correlation coefficients between NDVI and the other spectral indices

图5为重构后多光谱数据中红树林与陆地植被的直方图，对比图4可见在新光谱空间中，第1波段上红树林与陆地植被分布大致相近，但红树林分布范围更宽一些；第2波段上红树林的分布峰值位置约为-0.50，而一般陆地植被则小于-0.50；第3波段上二者分布特征最为接近；第4 波段上红树林分布峰值位置约为0.70，较之一般陆地植被偏大。这说明在重构数据光谱空间中红树林与一般陆地植被的差异性得到一定程度增强。

3.3 最小噪声分离变换

对重构后数据进行MNF，图6为变换结果的4个分量图像，根据MNF原理，变换后4个分量按照噪声水平升序排列，第1、2、3分量包含了主要的空间信息，而第4分量则包含较多的噪声信号。特别地，在MNF第3分量上，红树林与一般陆地植被表现出较明显的目视差异，前者具有更高的像元值。由图7可知，MNF后，红树林与一般植被直方图分布形态差异性得到提升，峰值分布区间的重合度进一步减小。具体表现为：第1分量上，红树林分布范围更宽，峰值位置约在9.30附近，陆地植被分布范围集中在10以上，峰值位于11.52附近；第2分量上，红树林分布范围亦较宽，峰值位置约在-1.23附近，陆地植被接近均值为-2.13的正态分布，主要范围在-3.77至1.77之间；第3分量上红树林与一般陆地植被分布差异最为明显，红树林分布中心偏右，峰值位置约为5.13，陆地植被分布峰值位置大致为3.62，且绝大部分分布在约4.36以下，而红树林一般高于此值（见图7中黑色虚线）；第4分量包含更多噪声信号，二者区分不大。这说明，基于光谱指数重构数据，MNF在分离了噪声信号的同时，由于利用了影像空间分布信息（公式（8）），进一步增强了目标地物的波谱差异，提高了红树林与一般陆地植被的可区分度。

图5 光谱指数重构数据样本直方图Fig. 5 Histograms of the spectral index reconstruction data

图6 重构后数据MNF的4个分量Fig. 6 The 4 MNF components of the reconstruction data

图7 光谱指数重构数据MNF结果样本直方图Fig. 7 Histograms of MNF components of the spectral index reconstruction data

通过上述分析，基于光谱指数重构数据MNF结果可以较容易地构建分类决策树，实现红树林提取。图8为本文构建的决策树，其中：

规则 1：NDVI ＞ 0.50，区分植被与非植被。

规则2：MNF3 ＞ 4.36，区分红树林与一般陆地植被。

国家卫星海洋应用中心在诸如“908”专项调查、广西海域使用本底库建设等工作中在该区域开展过多次遥感调查及现场核查，图9a是基于历史调查成果和专家经验支持下的山口红树林分布范围人工解译结果，图9b是基于本文决策树的红树林自动提取结果，可见绝大部分红树林均能被正确提取，其中专家解译的面积为9.44 km2，本文方法提取红树林的面积为9.65 km2，面积检测相对准确率可达97%。进一步在研究区内随机选取2 000个验证样本点，其中覆盖红树林样本点50个，本文方法检测出样本点44个，从而红树林总体检测精度约为88%。

图8 红树林提取决策树Fig. 8 Decision tree for mangrove extracting

4 总结

本文基于HY-1C卫星CZI 4波段影像，以广西山口红树林自然保护区为研究区，采用光谱指数重构CZI数据，经MNF建立决策树，实现了研究区红树林自动提取，研究结果表明：利用NDVI、NDWI、ARVI、CVSSR 4种指数重构CZI数据，减少了地形、大气、传感器等因素对植被光谱信息的影响，且重构后数据对植被、水体等更加敏感。MNF则进一步利用了影像上红树林与陆地植被的波谱及空间分布特征，变换结果中二者差异被显著增强，最终通过建立决策树提取了红树林分布信息。本文结合光谱指数重构和MNF方法，较好地解决了在CZI原始光谱空间中红树林与陆地植被难以区分的问题。因此，HY-1C CZI数据可以有效地用于红树林空间分布监测。

HY-1C CZI数据优势在于重访周期短、幅宽大，可快速重复获取大范围中等分辨率光学影像数据，但因波段较少，缺少对研究土壤与水分情况较重要的短波红外、中波红外或热红外等波段，且波段较宽，对典型植被波谱特征细节描述能力稍弱。本文研究方法可弥补CZI数据用于红树林监测的上述不足，后续还可通过形态学方法对提取结果进一步优化，提升自动检测准确率。需要指出的是，本文决策树规则阈值对于不同地区和时相的数据可能会有所调整，目前HY-1C卫星尚处于在轨测试阶段，相信卫星正式交付后，传感器数据的定标精度和质量稳定性会得到进一步保障，有利于后续对本文方法的验证和完善。

图9 研究区红树林提取结果分布Fig. 9 Mangrove extraction results