林用智，胡铁泷，李映潭，贺海明，王 杰
(西华师范大学地理科学学院，四川南充 637009)
植被指数是表征地表植被覆盖的简单有效的度量参数，被广泛用于植被类型监测[1]、叶面积指数监测[2-3]、覆盖度监测[4]和生物量估算[5]等研究.常见植被指数有差值植被指数(DVI)[6]、归一化植被指数(NDVI)[7]、增强型植被指数(EVI)[8-9]和土壤调节植被指数(SAVI)[10]等.植被指数主要是通过光学影像的近红外波段和红波段反射率而获取[11]，Sentinel-2影像被广泛应用于植被指数的提取与作物信息监测.Sentinel-2影像为多光谱影像，有13个波段，空间分辨率分为10、20、60 m，时间分辨率为5 d，Sentinel-2影像在红边范围有3个波段(Band 5、Band 6、Band 7)，3个红边波段能有效地监测植被信息[12]，然而，Sentinel-2影像易受云雾污染，云层的存在导致传感器无法获取地物反射信息，使植被指数获取较为困难，因此获取植被指数时，要考虑去除云雾的影响.传统的去云方法通常分为两大类：多时相的光学影像去云和SAR与光学融合去云.国内外学者在去除光学影像中的云雾、获取云层下的植被指数方面做了大量研究，Cerra等人首次将基于稀疏解混的去噪应用于Sentinel-2影像的去云，该算法的速度较快、效率较高[13]；
合成孔径雷达具有全天候、全天时监测的特点，其波长较长，受云雾影响较小，且微波能穿透植被，获取植被的种植结构信息，使其在植被监测中有一定优势.光学和SAR数据各自具有优势，有效地结合它们的优势能提高数据可用性和协作性，因此有学者利用SAR与光学融合去云[14-15]，以机器学习去云为主，机器学习去云包括回归方法和深度学习去云算法，Meraner等人设计了一种融合SAR和光学数据的深度残差神经网络构架去除Sentinel-2影像中的云雾，并最大限度地保留了原始信息[16]；
Filgueiras等人利用8种回归算法，结合Sentinel-1和Sentinel-2影像，验证了4种建模方法，得出8种回归算法中，随机森林方法获取的精度最好[17].
以上研究中，基于多时相的光学影像去云算法需要多景影像数据，影像之间存在的时间间隔，加上降雨、光照等因素的影响，植被的光谱曲线会发生变化，导致多时相去云算法的误差；
机器学习算法虽然能最大程度减小光学与SAR的成像时间间隔，但是其结果通常会丢失部分光谱特征，或者获取的光谱曲线发生巨大的变化，从而导致植被指数获取呈现巨大的误差.如何有效地去除光学影像的云像元并保留云下像元的光谱特征是计算有云影像植被指数的关键.针对上述问题，本文利用Sentinel-1影像构建出6种匹配模型，分别对6种匹配模型采用全局光谱最小距离法去除Sentinel-2影像上的云与阴影像元，获取6景去云后的Sentinel-2影像，恢复云下地物的光谱曲线，然后获取植被指数，最后验证去云与阴影区域的植被指数提取的精度.
1.1 试验区概况
研究区位于新疆维吾尔自治区乌苏市(44°27′～44°39′N，84°29′～84°45′E)，属温带大陆性气候，四季分明，地形平坦，主要农作物为棉花、小麦和玉米等，且种植面积广，植被覆盖面积较大.因而本研究选择该区域作为研究区.
1.2 数据获取及处理
本文的研究数据来自ESA数据共享网站(https://scihub.copernicus.eu/dhus/#/home)，包括SAR数据和光学数据.SAR数据为2020年9月11日的Sentinel-1影像GRD数据；
光学影像为2020年9月3日无云的Sentinel-2B和2020年9月8日有云的Sentinel-2A影像.
本文选择的Sentinel-2影像为L1C数据.使用SNAP 8.0和ENVI 5.3对数据进行预处理，步骤如下：(1)用Sen2cor模块对影像进行大气校正；
(2)再用Sentinel-2 super-resolution模块进行空间分辨率增强；
(3)在ENVI中合成多光谱影像.由于Sentinel-1A、Sentinel-2A两颗卫星过境时间不一致，且研究区的夏季没有时间点完全一致且能使用的影像，而有云的Sentinel-2影像时间为2020年9月8日，SAR影像和光学影像相差3天的时间对影像影响较小，所以本文采用2020年9月11日的Sentinel-1影像.用SNAP软件对SAR影像进行预处理：(1)应用轨道文件；
(2)辐射校正，输出sigma0波段；
(3)斑点滤波；
(4)多普勒地理编码，输出本地入射角(LIA，Local Incidence Angle)；
(5)同时生产双极化的dB波段；
(6)最后将Sentinel-1A、Sentinel-2A影像进行配准.
图1 Sentinel影像Fig.1 Sentinel images
配准完成后，裁剪出研究区域的Sentinel影像.如图1所示，(a)为无云的Sentinel-2B影像(RGB：8A，4，3)，(b)为有云的Sentinel-2A影像(RGB：8A，4，3)，(c)为Sentinel-1影像(VV，VH，VV/VH合成).Sentinel-2A影像经大气校正后，产生云覆盖的概率图，将其转换为二值影像，再手动添加未被探测的被云污染的区域，如图2.
1.3 研究方法
1.3.1 Sentinel-2A影像去云处理 胡铁泷[18]等人通过Sentinel-2A影像去云方法有效地结合了SAR数据和光学数据的优势，达到Sentinel-2A影像去云与阴影的目的，如图3所示，其基本流程如下：通过预处理后的Sentinel-2A影像获取云掩膜图，利用云掩膜图将研究区的Sentinel-1A影像分为云和非云两类像元，然后采用欧式距离获取Sentinel-1A影像中离每个云像元(X′)距离最近的非云像元(Y′)的位置，最后根据获取的非云像元位置，在Sentinel-2A上找到相同位置的非云像元(Y)，再用非云像元的光谱曲线代替距离最近的云像元(X)的光谱.通过以上的算法，寻找Sentinel-1A影像上的最相近的位置来替代Sentinel-2A影像的云光谱曲线，能够有效地恢复地表物质的光谱曲线，相比机器学习算法，其更具有物理意义，因为其影像上大量存在相似光谱曲线.
图2 云与阴影的感兴趣区域Fig.2 Clouds and shadows region of interest图3 去云与阴影流程图Fig.3 Flow figure of removal cloud and shadow
以上算法的核心是实现Sentinel-1A影像上像元的正确匹配，Sentinel-1A影像只有两个极化波段，分别为VV与VH，由于类内光谱曲线呈现较大的变化[19-21]，这样会增大匹配的误差，从而导致后期的植被指数的计算精度较低.为了增加Sentinel-1A影像像元的匹配精度，需要考虑增加波段信息，为此，本文通过不同的特征运算，增加了数个特征，同时对这些特征进行组合，产生六种不同的匹配模型，验证不同波段模型组合的匹配精度.其中，增加的特征运算包括：[1]VV与VH极化波段的dB值计算，如式(1)所示；
[2]Sentinel-1A影像的比值植被指数(RVI)[22-23]，如式(2)所示；
[3]Sentinel-1A影像归一化比值指数(NRPB)[17]，如式(3)所示；
[4]Sentinel-1A影像本地入射角(LIA)，这个波段是对Sentinel-1A影像地理编码后生成的，将其导入即可[17].
(1)
(2)
(3)
1.3.2 植被指数 获取去云与阴影影像之后，计算云与阴影覆盖区域的植被指数，包括归一化植被指数(NDVI)、土壤调节植被指数(SAVI)和增强型植被指数(EVI).NDVI是根据植被在近红外波段和红波段的梯度差得出的公式，NDVI的优势在于能消除大部分地形和大气产生的影响，但易受土壤等因素的干扰，因此，Huete[10]等人在此基础上考虑了土壤的影响，加入了土壤调节参数，提出了SAVI.同时，LIU[9]等人为了解决NDVI存在的问题，提出了EVI，该指数引入的土壤调节参数和大气修正参数减少了土壤和大气的影响.上述植被指数计算公式为
(4)
(5)
(6)
其中：ρNIR和ρred分别为近红外波段反射率和红波段反射率；
ρblue为蓝波段反射率；
C1为大气修正红光校正参数，取6.0；
C2为大气修正蓝光校正参数，取7.5；
L为土壤调节参数，取1[22].
1.3.3 精度评价 本文采用相关系数(R)、平均绝对误差(MAE)和均方根误差(RMSE)3个指标作为模型精度的验证[24-26].
(7)
(8)
(9)
图4 无云的Sentinel-2B影像计算的3种植被指数Fig.4 Cloud-free Sentinel-2B image calculated three kinds of vegetation indices
为了验证6个匹配模型所计算的植被指数的精度，使用2020年9月3日无云的Sentinel-2B影像计算出的EVI、NDVI和SAVI作为真值，其结果如图4.结合去云算法和6个匹配模型的光谱特征去除Sentinel-2影像上的云与阴影像元，获取6景去云与阴影后的Sentinel-2影像，如图5所示，(a)～(f)分别与匹配模型的01～06编号相对应.采用相关系数、均方根误差和平均绝对误差对植被指数的整体精度和6个匹配模型的精度进行评价.首先评价去云算法，对比图1(a)和图5(a)～(f)可以看出，去云算法能较好地修复受云与阴影污染的像元，恢复地物的光谱曲线，清晰地呈现出地块的基本形状.观察发现，图5绿色方框中的水体区域恢复出了部分土壤的像元，这与实际情况不符合，分析图1的Sentinel-1A影像之后，发现其有部分土壤的信号，主要原因是水体较浅且浑浊(最近有一次强的降雨)，导致后向散射系数增大所致，后期可结合前期的Sentinel-2A影像，探测出相关的水体，然后掩膜掉，以提高恢复影像的质量.
图5 去云与阴影后的Sentinel-2A影像Fig.5 The Sentinel-2A image of removal cloud and shadow图6 6种匹配模型计算的EVIFig.6 EVI calculated by six kinds of matching models
对比图4(a)、图4(c)、图6、图8和表1分析发现，在3种植被指数的整体精度中，云与阴影覆盖区域所计算的EVI和SAVI与真值较为接近，且图像色差较小.与真值相比，SAVI整体的均方根误差和平均绝对误差分别小于0.092 71和0.071 40，EVI整体的均方根误差和平均绝对误差分别小于0.124 31和0.093 24，SAVI和EVI的相关系数均大于0.8，说明SAVI的整体精度高于EVI.结合图4(b)和图7可以看出，云与阴影覆盖区域计算的NDVI与真值的图像相比，颜色略有差异，同时分析表1可以得出，去云与阴影后计算的NDVI的RMSE和MAE分别大于0.14和0.11，研究气象资料发现，9月3日至9月8日有降雨，对植被的影响较小，但降雨对土壤影响较大，因此，没有引入土壤调节参数的NDVI的精度要略低于SAVI和EVI.
图7 6种匹配模型计算的NDVIFig.7 NDVI calculated by six kinds of matching models图8 6种匹配模型计算的SAVIFig.8 SAVI calculated by six kinds of matching models
通过表1，综合对比六个匹配模型的评价指标，可以得出04_VVdB_VHdB_RVI匹配模型在三种植被指数中的精度都要高于其余五个匹配模型，在SAVI中的精度最高，其RMSE和MAE分别为：0.09236和0.07117.在六个匹配模型中，01_VV_VH_VVdB匹配模型的精度较低，但与其余五个匹配模型相比，精度相差较小，且六个匹配模型在三个植被指数中的相关系数都大于0.8，说明六个匹配模型均适用于植被指数的计算，且能获取较好的结果.
表1 6种匹配模型计算的植被指数精度对比Tab.1 Comparison of precision of vegetation index calculated by 6 matching models
3.1 结论
为了有效地恢复受云与阴影污染的光学影像，获取云与阴影覆盖区域的植被指数，本文利用Sentinel-1影像构建6种匹配模型，结合去云算法，获取6景去云与阴影后的Sentinel-2影像，再分别利用6景去云与阴影的影像计算3种植被指数，并进行精度评价.通过上述分析，得出结论如下：
(1)去云算法能有效地去除云与阴影对光学影像产生的影响，恢复地表物质的光谱曲线，为计算受云与阴影污染的光学影像的植被指数提供算法支持.
(2)六个匹配模型均适用于植被指数的计算，且能获取较好的结果，其中04_VVdB_VHdB_RVI模型在三种植被指数中的精度均为最高.
3.2 讨论
六个匹配模型中，01、02模型以VV和VH极化为主，结合dB值构建模型，而03～06模型以VVdB和VHdB为主，结合雷达指数、本地入射角等特征构建模型，用三个特征构建模型时，其精度已趋于稳定，模型包含四个特征时，精度无明显提升.随着特征的增加，算法运行的时间也会增加.该算法适用于SAR与光学影像行列数与分辨率相同的遥感数据，基于SAR和光学影像融合的去云算法计算植被指数能获取较好的效果，但仍有改进空间，以后可以考虑对影像进行归一化处理[19-21]，降低类内的光谱变化.去云的过程中，本文采用全局搜索进行匹配，虽然精度较高，但是效率极低，后续考虑地表精细化分类(相邻时间无云影像数据)，以减小搜索的范围.
猜你喜欢植被指数阴影波段你来了，草就没有了阴影文苑(2020年11期)2020-11-19影中国诗歌(2019年6期)2019-11-15基于植被指数选择算法和决策树的生态系统识别农业机械学报(2019年6期)2019-06-27AMSR_2微波植被指数在黄河流域的适用性对比与分析水土保持研究(2018年5期)2018-10-12河南省冬小麦产量遥感监测精度比较研究中国农业信息(2018年2期)2018-07-28基于PLL的Ku波段频率源设计与测试电子制作(2018年2期)2018-04-18让光“驱走”阴影中学生数理化·八年级物理人教版(2017年11期)2017-04-18M87的多波段辐射过程及其能谱拟合高师理科学刊(2016年8期)2016-06-15阴影魔怪数学大王·中高年级(2016年4期)2016-05-14日常维护对L 波段雷达的重要性西藏科技(2015年4期)2015-09-26