Sentinel 2A影像去云下的丘陵地区植被覆盖度反演

胡铁泷,蒋良群,王 杰

(西华师范大学 国土资源学院,四川 南充 637009)

植被是连接土壤、大气和水分的纽带,具有减少雨滴击溅、减缓地表径流、增加土壤保固等功能[1]。植被覆盖度是指植被的叶、茎、枝等在地面的垂直投影面积占统计区总面积的百分比[2]，通常用植被覆盖度来表示地表植被的覆盖情况。目前,获取植被覆盖度的方式有地面监测和遥感监测,其中地面人工监测耗时耗力,因此通过遥感数据反演植被覆盖信息的方法被广泛应用[3]。南方丘陵地区常年多云雾,光学影像受到云雾的污染,很难准确反演植被覆盖度,因此利用光学影像进行植被覆盖度的反演时,通常要考虑云雾的影响。合成孔径雷达能适应各种恶劣的天气全天候成像,同时由于波长较长能穿透云雾,可较清晰地呈现地表的景观格局,被广泛地应用于生产实践之中[4]。虽然合成孔径雷达具有较大的优势,但相对于目前在轨的多光谱资源卫星如Landsat 8[5]、Sentinel 2A[6,7]等,极化组合的波段较少、斑点噪声多、人工目视解译难,涉及较复杂的辐射传输模型[8],大大限制了合成孔径雷达的广泛利用。光学与SAR影像各自具有优势,如何有效地集合它们的优势是目前研究的热点[9,10]。通常情况下,这类算法需要多时相的合成孔径雷达数据,但会涉及到变化检测的问题,增加植被覆盖度提取的难度和不确定性。单时相的SAR与光学数据的融合研究非常具有科研价值,它简化了融合算法的复杂程度,可兼顾时效性。植物的光谱曲线随着时间不断变化,同时由于地形、光照、水分等因素的影响,植被光谱的变化较大,增加了植被覆盖度提取的难度。如何有效地减弱地物光谱的变化,近几年大量的光谱变化混合分析模型被提出[11-17]，大致分为两类,即确定性模型与统计性模型。虽然光谱变化混合分析模型能解决端元间的变化,但是对相似度较高的端元内的光谱变化大都很难解决[18,19]。

针对上述研究中的问题,本文基于一景Sentinel 1A和Sentinel 2A影像,避免多时相数据的复杂性与时效性,利用全局光谱最小距离法去除Sentinel 2A影像的云,同时恢复Sentinel 2A云覆盖下的地物光谱曲线,然后采用光谱归一化技术减弱地物内的光谱变化,在此基础上采用扩展线性混合模型反演植被覆盖度。

1 研究区概况

本文选择夏季植被覆盖茂盛、地形复杂的四川丘陵区作为研究区。南充市位于四川省东北部,地处嘉陵江中游,属于亚热带湿润季风气候区,四季分明、雨热同期。研究区位于四川省南充市,属丘陵地貌,介于30°74′—31°9′N、106°08′—106°35′E。研究区夏季植被覆盖茂盛、类别复杂,光学遥感影像常年多云覆盖,适合研究光学遥感影像的去云和植被覆盖度的反演。

2 数据源及预处理

数据源采用Sentinel 1A和Sentinel 2A影像,其中Sentinel 1A影像为干涉宽模式(IW)下的地距产品(GRDH),而Sentinel 2A影像为L1C格式,需要进行大气校正。从欧洲航天局(ESA，https://scihub.copernicus.eu/dhus/)获取南充市附近2018年7月15日的Sentinel 1A、2018年7月21日的Sentinel 2A遥感影像。由于Sentinel-1A和Sentinel-2A两颗卫星过境时间不一致,夏季没有时间点完全一致且能利用的影像,加之在7月份相差6天的时间对影像的影响较小,最终选定了这两景影像。本文的数据处理流程见图1。图1的右部分为Sentinel 2A影像去云处理(结合Sentinel 1A影像),左部分为Sentinel 2A影像的光谱归一化与光谱混合分析,从而获取植被覆盖度,并进行精度评价。

图1 数据处理流程图

Sentinel 2A影像的预处理主要通过ENVI 5.3和ESA提供的SNAP 6.0(Sentinel Application Platform)软件完成。首先，运用SNAP软件的Sen2cor模块对Sentinel 2A L1C数据进行大气校正处理,输出10m、20m、60m分辨率的11个波段(Band 1、Band 2、Band 3、Band 4、Band 5、Band 6、Band 7、Band 8、Band 8a、Band 11和Band 12)的L2A地表反射率数据，利用空间分辨率增强工具(Sentinel-2 super-resolution)将L2A所有波段增强到10m,使影像的波段保持相同的分辨率;然后将增强的影像重采样至20m,以保持与Sentinel 1A预处理后的空间分辨率一致;最后利用ENVI5.3将11个波段数据合成多光谱影像。

运用SNAP软件对Sentinel 1A影像进行预处理。由于地形导致辐射的偏差,本次试验对Sentinel 1A数据采用地形辐射校正。预处理步骤:①应用轨道文件(Apply Orbit File)自动下载精确的轨道信息文件。②对影像进行辐射校正(Calibrate)。输出beta0波段,因后续地形辐射校正算法需要输入beta0波段;对影像进行斑点滤波(Frost Filter),消除影像上的噪声。③进行多视(Multilooking)处理。将影像重采样到20m分辨率,在保持数据有效的情况下,同时提高数据处理的效率,为Sentinel 2A去云处理做准备。④进行辐射地形校正(Radiometric Terrain Flattening)。自动下载30m的数字高程模型。⑤进行距离多普勒地理编码(Range-Doppler Terrain Correction),获取准确的地理位置。为提高后续的匹配运算精度,分别生产VH与VV极化的辐射强度(db)波段。

为减弱两类影像的配准误差,完成Sentinel 2A和Sentinel 1A的预处理后,利用SNAP软件的地理位置配准工具,以Sentinel 1A作为主影像,Sentinel 2A为辅影像,将两类影像进行精确地配准,为后续的去云算法做好数据准备。在Sentinel 2A影像上裁剪出研究区,选择12、8、2这三个少云的波段合成影像，见图2；其他波段大都呈现大量的云,见图3显示的波段1。Sentinel 2A影像经过大气校正后,产生云覆盖的概率图,将其转换为二值影像,手动添加未被探测到的云污染区域,见图4。同时，裁剪出相应区域的Sentinel 1A影像,见图5(VH,VV_db,VV合成)。从图5可见,虽然Sentinel 1A影像经过光谱滤波与多视处理后纹理非常清晰,但是影像斑噪仍然非常明显。

图2 Sentinel 2A影像彩色图

图3 Sentinel 2A波段1

图4 Sentinel 2A影像云掩模图

图5 Sentinel 1A彩色合成图

3 研究方法

3.1 Sentinel 2A影像去云处理

由于光学影像Sentinel 2A上存在大量的云,而微波影像Sentinel 1A无云,因此利用Sentinel 1A影像的光谱特征,将Sentinel 2A的云雾去除。具体算法流程为:首先，通过Sen2cor大气辐射校正模型获取Sentinel 2A影像的地表反射率和云掩模产品,在此基础上,用云掩模产品(二值影像)将相应的Sentinel 1A影像标记成云与非云两类像元；然后在整个Sentinel 1A影像空间采用欧式距离计算离每个云像元距离最短的非云像元的位置；最后根据在Sentinel 1A获取的非云像元位置,在Sentinel 2A上找到相同位置的非云像元,再用找到的Sentinel 2A非云像元位置的光谱曲线代替距离最短的云像元的光谱。以上处理可有效结合两类数据优势,达到Sentinel 2A影像去云的目的,为后续的端元变化解混准备数据源。去云处理主要通过IDL代码来实现,由于篇幅的限制,文中没有罗列去云的相关代码。

去云的核心部分见图6。首先在Sentinel 2A找到云像元的位置,在Sentinel 1A找到相同位置的像元,然后在整个Sentinel 1A影像空间采用欧氏距离计算出距离像元最小距离的非云像元的位置,再在Sentinel 2A上找到与位置相同的非云像元,最后用Sentinel 2A非云像元的光谱曲线替代有云像元的光谱。对Sentinel 2A有云的像元,运用程序依次循环以上的过程,最终达到去云的效果。

图6 去云流程图

3.2 光谱归一化

常见的光谱解混模型是首先固定端元组,然后应用解混算法获取各个地类的光谱曲线。如果端元组内光谱的相关性过高、变化较大,那么选择合理的端元难度增大,从而解混精度会降低[18,19],因此采用光谱归一化技术减弱端元组内的变化,具有十分重要的意义。光谱归一化就是各个端元光谱除以它们各自的范数,在不改变端元之间的相关性和它们在高维特征空间的相对位置的情况下,压缩数据空间,减弱端元变化的绝对量,减弱端元组内的差异。

3.3 光谱混合分析

由于空间分辨率的限制,常常导致植被与其他地物以混合像元形式呈现在影像上。通过前期的光谱归一化技术,虽然能降低植被内的光谱变化,但是不能从根本上解决植被与其他地物的混合状态,因此混合像元分解技术被采用。在此次试验中,扩展线性混合模型[15,20]被用于反演植被的覆盖度。假设包含N个像元的L波段遥感图像包含P个端元,用x∈RL×N、A∈RP×N、S∈RL×N、E∈RL×N分别表示像元矩阵、丰度矩阵、端元矩阵和误差项。那么对单个像元,线性混合模型可表示为:

(1)

式中,xk(k=1…N)为一个像元;ek(k=1…N)为误差向量。通过对丰度进行限制,使其具有较强的物理意义,再增加两个限制条件:

(2)

本文将式(1)与式(2)称为全约束线性最小二乘(FCLS)算法[21]。当考虑每个像元的端元变化时,式(1)可改写为:

(3)

当考虑到光照因子为主要影响因素时,有fpk(s0p)=Ψpksop。式中,sop为参考端元,式(3)可改写为:

(4)

式中,φk=diag(φk),φk∈RP×P为值的对角矩阵,相当于部分约束最小二乘法(SCLSU)[17]。式(4)只是将光照因子作为影响端元变化的主要因素,当同时考虑像元间的空间邻域信息和光谱变异问题时,在保留ELMM框架下可写为:

(5)

=λA(PHh(A)P2,1+PHv(A)P2,1)+LRP×N+(A)+μT(ATlP-lN)

(6)

(7)

4 结果及分析

利用Sentinel 1A和Sentinel 2A结合去云,得到的结果见图7。从图7可见,在上段河流中出现了大量的土壤像元,再结合图5,发现Sentinel 1A预处理后的河流中也有很多强反射率的土壤像元,说明去云算法有效利用了Sentinel 1A影像的特征。虽然去云算法可能会导致恢复的图像出现部分错误,但是未影响植被覆盖度的提取。从图7可见,去云后的Sentinel 2A影像植被出现斑点状且不连续,主要是由Sentinel 1A影像的斑点状引起。为了减弱Sentinel 1A影像的斑点对植被覆盖度提取的影响,再对去云后的Sentinel 2A影像进行中值滤波平滑,以减小Sentinel 2A影像噪声的影响(图8)。从图8可见,中值滤波的影像更趋于平滑,压制了部分噪声。

图7 去云后的Sentinel 2A影像

图8 中值滤波后的Sentinel 2A影像

在进行去云与中值滤波处理后,对Sentinel 2A影像进行光谱归一化,以减弱地物内的光谱差异,结果见图9。对比图8和图9可见,建筑物、土壤、水体、植被呈现大面积的一致性,说明光谱归一化技术能减弱地类内的光谱变化,但其中建筑物与土壤的光谱还存在一定的光谱变化。完成影像光谱归一化后,进行端元提取与光谱混合分析。经过对研究区的实地调研,发现多个水塘有水生植物,可能会影响植被覆盖度反演的精度,为此结合ENVI5.3软件与以往的土地利用分类图提取研究区的水体,掩膜掉Sentinel 2A影像上的水体,只评价陆地上的植被覆盖度。端元提取算法采用了非降维的体积迭代[22],共提取了33个端元,并获取了端元的位置,最终确定了5类端元的光谱曲线,见图10。由于篇幅的限制,本文未对大量的植被光谱归一化前后进行对比,可参考相关文献[18,19]。

图9 光谱归一化后的Sentinel 2A影像

图10 归一化前后所提端元的光谱曲线

采用ELMM模型定量植被的覆盖度,其中ELMM算法需要数据的紧凑度参数、丰度系数参数和尺度比例参数。为了减弱随机输入参数的误差,我们设置了20组参数的组合,取RMSE最小值对应的参数作为最终参数,将其解混结果作为ELMM算法的最终结果。试验发现,其组合(λs=0.1,λA=0.001,λψ=0.0)最优,反演的结果见图11(右)。为了验证植被覆盖度解混的精度,以近几年的土地利用分类图和时间点相近的无云Sentinel 2A影像作为参考,使用Sentinel 2A数据(10m分辨率影像,SNAP软件Sentinel 2A空间分辨率增强运算)进行神经网络分类,其中云与阴影覆盖的区域采用手工矢量化改正,将分类的结果作为参考影像对解混的结果进行精度评价(同时去掉建筑物与裸地等非植被类型),结果见图11(左)。评价指标采用均方根误差(RMSE)、相关系数(R)。精度评价后其RMSE约为0.14,R的值约为0.95,由此可知去云后影像植被覆盖度的反演结果精度较高,接近地表植被覆盖的真实情况。

对框架各部分的程序运行效率进行分析,在电脑处理器为Inter(R) Xeon(R) CPU E5-2690 v4@2.6GHZ,运行内存128G和64位Windows10的操作系统,以及MATLAB版本为2016a 64bit(主要ELMM模型的运行)、ENVI 5.3 64bit(主要进行去云模块的实现)下,对781×1305大小的Sentinel 2A影像进行处理,其中ELMM模型用时大约470秒,而去云模块用时大约为7.5h,说明逐像元循环去云模块的效率极低。今后应考虑多台电脑进行集群与分块计算,以提高处理效率。

图11 植被覆盖度的参考影像与反演结果

5 结论与展望

本文针对南方光学影像多云的特征,提出了结合SAR去除光学影像上的云;针对丘陵地区复杂的地形、水分、光照等因素对植被的影响,对影像和端元进行归一化,以减弱因端元变化导致的解混误差,并采用扩展线性模型进行植被覆盖度反演和精度评价。通过以上分析,得出如下结论:基于去云的Sentinel 2A影像反演植被覆盖度的效果较好,接近地表植被覆盖的真实情况。由于Sentinel 1A影像的特点,部分地物反演的面积不准确,特别是较破碎化的裸地,但对植被的覆盖度影响不大。本次试验仅采用了一景Sentinel 1A影像,以后可考虑多景Sentinel 1A数据进行变化检测,同时进行多波段合并与配准,增加极化波段,提高Sentinel 1A的光谱匹配能力。为了提高程序的运行效率,以后将数百台计算机进行集群与分块计算。