基于多源传感器的FUI水色指数提取与应用
2021-09-24王一飞宫兆宁
王一飞，宫兆宁，张 园，苏 朔
(1.首都师范大学资源环境与旅游学院，北京 100048；2.三维信息获取与应用教育部重点实验室，北京 100048；3.资源环境与地理信息系统北京市重点实验室，北京 100048；4.北京市城市环境过程与数字模拟国家重点实验室培育基地，北京 100048;5.生态环境部卫星环境应用中心，北京 100094)
0 引言
湖库是地表系统最重要的组成部分之一，不仅可为人类提供饮用水资源，还具有关键的生态功能[1-2]。在过去近百年中，密集的人类活动给湖库水体带来一系列的环境危机，包括水质恶化、渔业减产和生物多样性丧失[3]。因此，研究湖库水质现状和变化趋势具有重要的科学和现实意义。而与传统的野外站点水质监测相比，遥感技术具有快速、连续、动态、大范围对地观测的特点，是一种有效的监测手段[4-5]。
影响水体富营养状态的主要是叶绿素a浓度、水体透明度、总氮、总磷和化学需氧量等。其中，叶绿素a浓度和水体透明度可以通过遥感技术进行监测。水体营养状态通常分为贫营养、中营养和富营养3个级别。通过与水体营养状态有关的指标以及指标间的相互关系，对水体的营养状态做出准确的判断就是水体营养状态评价[6]。Carlson[7]于1977年提出了基于浮游植物生物量的营养状态指数(trophic status index，TSI)，对水体富营养化程度进行了 0～100 的量化，通过水中叶绿素 a 浓度、透明度或总磷含量进行计算。后续研究发现，基于叶绿素 a浓度的 TSI指数模型较为可靠[8-9]。TSI指数因为其计算简便结果可靠而被广泛应用[10]。
水体中的3种主要光学物质——叶绿素、非藻类悬浮物和有色可溶性有机物(colored dissolved organic matter,CDOM)以及水分子本身的吸收和散射作用共同决定了水体呈现出的颜色[11]。水体颜色是传统水质监测的一项重要指标，已被世界气象组织(World Meteorological Organization，WMO)列为湖泊的基本气候变量之一。天然水的颜色分类始于19世纪底Forel-Ule比色表的引入。1890年，Forel 提出了一种天蓝计和黄原计的组合的颜色标准，用于从蓝色水体中细分出绿色水。1892年，Ule 补充了蓝绿色到棕色，将天然水体由深蓝到红棕色划分为21个颜色级别，这种颜色区分表被称为 Forel-Ule 量表[12]。
基于遥感影像的 FUI 水色指数(Forel-Ule index)的提取研究始于 2012 年欧洲Citclops 项目。该项目中，Wernand等使用有色化学溶液配制了 Forel-Ule 比色表的不同颜色，完成了对不同颜色的光谱分析，构建出完整的 FUI 指数对应的 CIE-XYZ 色度坐标点[13]；基于历史存留的FUI水色现场测量数据，建立了大洋水体 FUI 水色指数与叶绿素浓度的关系，分析了大洋水体1889—2000年间的叶绿素浓度变化[14]；提出了基于 MERIS 遥感影像的水体 FUI提取算法[15]。在国内的研究中，Li等[16]开发了基于 MODIS 地表反射率产品的 FUI 水色指数反演算法，并应用于中国十大湖泊的水色和营养状态的监测；Chen等[17]使用基于Landsat8的FUI反演算法提取了长江流域2013—2018年的FUI 水色并分析了长江流域富营养状态的变化。
这些研究表明FUI水色指数在湖库环境监测上具有很好的研究潜力和优势。本文以官厅水库为研究对象，收集了2016—2020年的Landsat8 OLI和Sentinel-2数据，并使用不同的反演算法提取其对应的FUI水色参数，分析FUI水色指数在不同遥感影像反演结果的可比性。分别使用实测叶绿素a数据和国际海洋水色协调工作组(The International Ocean-Colour Coordinating Group，IOCCG)模拟数据计算所得的TSI 营养状态指数与FUI水色指数建模，分析FUI水色指数表征水体富营养状态的可行性。
1 研究区概况及数据源
1.1 研究区概况
本文以官厅水库作为重点研究区，官厅水库位于河北省张家口市怀来县和北京市延庆县界内 (图1)，地处N40°18′～40°26′，E115°37′～115°51′，是新中国成立修建的第一座水库，设计总库容41.6亿m3，水库流域总面积为4.34万km2，平均水深约7.6 m[18]。20世纪60年代，国家大力发展工业，在官厅水库上游的桑干河、洋河以及永定河流域建设了大批化工、冶金、造纸工厂，大量污水排入河道，官厅水库初次被污染。20世纪70年代官厅水库经过治理之后，上游流域严重的水污染基本得到了遏制，水库水质恢复到了饮用水源地水质标准[19]。但是到了20世纪80—90年代，乡镇企业扩大规模，地方工业发展迅速，生产的污水直接入河，官厅水库再一次被污染。2015年4月16日，北京市环保局发布《2014年北京市环境状况公报》，公报显示官厅水库水质不符合饮用水要求，水库中的化学需氧量、氟化物和高锰酸盐含量超标[20]。
图1 研究区概况(Landsat8 B5(R)，B6(G)，B4(B)合成影像)Fig.1 Map of the study area
1.2 数据源
1.2.1 实测数据
实测数据包括水质数据和水面光谱数据。各项水质参数使用Hydrolab MS5水质多功能探头测得，MS5多功能探头能同时测得水温、pH值、总溶解固体(total dissolved solids ,TDS)、叶绿素a等参数。在每个采样点间隔10 s测一次数据，同一个点测3组数据求平均值得到最终数据。2019年9月5日官厅水库实地采样点分布如图2所示。
图2 官厅水库采样点(Sentinel-2 B8(R)，B4(G)，B3(B)合成影像)Fig.2 Guanting Reservoir in-situ measured points
光谱测量仪器使用的是FieldSpec 3背挂式野外高光谱辐射仪，光谱范围为350～2 500 nm，350～1 000 nm间光谱的采样间隔为1.4 nm，1 000～2 500 nm光谱采样间隔为2 nm。数据采样时间选择在2019年9月初，该段时间光照充足，采样点晴朗无云，非常适合光谱采集。测量时间选择光照最为充足的10:00—14:00。在每个采样点均对标准白板重新定标优化，依次采集标准板、天空光和水体光谱，其中水体光谱的测量使用水面以上测量法[21]，公式为：
(1)
式中：Rrs为水体光谱反射率值；Ssw,Ssky和Sp分别为光谱仪面向水体、天空和标准板时的测量信号值；ρp为标准板的反射率；r为气水界面反射率，查阅文献可知，平静水面可取r=2.2 %[22]。
1.2.2 遥感数据及其预处理
本文主要选用Sentinel-2影像为研究数据，选用Landsat8 OLI为对比数据，用来验证FUI水色在不同传感器间的可比性。
选取北京市2016—2020年不同季度的20景Sentinel-2和1景Landsat8 OLI影像作为卫星研究数据，均挑选研究区无云遮盖的影像。Sentinel-2影像均从欧州航天局网站(https://scihub.copernicus.eu/dhus/#/home)下载，Landsat8 OLI 影像从美国地质调查局网站(https://earthexplorer.usgs.gov/?tdsourcetag=s_pcqq_aiomsg)下载。
在欧州航天局网站下载到的Sentinel-2影像是level-1C数据，已经过辐射校正和几何纠正，只需对其进行大气校正。实验证明，在针对Sentinel-2的几种大气校正方法中，与SMAC模型、6s校正模型相比，Sen2Cor模型校正效果最好，校正后的光谱曲线与地面实测的光谱曲线具有更高的拟合度[23]，因此本文所使用影像均使用Sen2cor 方法进行大气校正，为了后续的波段计算，将各波段重采样至同一空间分辨率，10 m和60 m空间分辨率的波段均重采样至20 m。Sentinel-2 影像的预处理和后续计算使用SNAP 7.0软件。
Landsat8 OLI影像的大气校正采用Envi 5.3软件里的FLAASH大气校正模块，FLAAAH模型结合了MODTRAN4的大气辐射传输编码，可以方便地选用不同的大气模型和气溶胶类型，对地表反射率进行计算[24]。
为了验证影像大气校正的准确度，使用ASD实测光谱，根据影像的光谱响应函数，即基于影像各个波段的最大值一半处的波宽 (FWHM)的波长中心和高斯滤波函数,将1 nm间隔的野外ASD波谱重采样到影像对应的波段。重采样公式如下[25]：
(2)
2 研究方法
2.1 水体提取
归一化水体指数(normalized difference water index，NDWI)是基于归一化植被指数NDVI改进的[26]，根据水体在可见光波段和近红外波段的反射强度，尽量抑制图像中的植被信息从而增强水体信息。NDWI的表达式为：
(3)
式中：Green和NIR分别为绿光和近红外波段的反射率。
基于NDWI指数和决策树方法，通过找到合适的NDWI阈值来分离水体和非水体，然后转化为水体矢量对影像进行裁剪。提取流程如图3所示。
图3 水体提取方法Fig.3 Water extraction method
2.2 FUI水色指数提取
2.2.1 CIE色度空间与FUI查找表
色彩空间是指用一种客观的方式描述人眼对颜色的感觉，通常需要首先定义3种主要颜色，然后利用颜色叠加模型，即可叙述各种颜色。为了对颜色进行量化表示，国际照明委员会(Commission Internationale De L’Eclairage，CIE)发展了一套CIE-XYZ颜色系统[27]。CIE-XYZ系统以X，Y，Z代替了R，G，B，使色度系统中的光谱三刺激值XYZ全部为正。XYZ-RGB之间的转换关系为：
(4)
CIE-XYZ系统中三刺激值计算公式为：
(5)
CIE-XYZ系统中三刺激值对定义颜色很有用，但是无法直观地对应自然界中的不同颜色。因此，CIE在1931年规定了二维色度空间图，二维色度图上的颜色与亮度无关。色度图上的坐标x,y,z从三刺激值X，Y，Z计算得来，公式为：
(6)
x+y+z=1，所以用x,y两个值就可以确定一种颜色，因此CIE-xy色度空间(图4)可以表示可见光范围内的任一颜色，每种颜色都能在色度空间中找到其对应的色度坐标(xM,yM)。图中白点W
随着社会的发展，老婆饼传到了大都市上海，成为当地一味名食，著名京剧表演艺术大师梅兰芳生前喜食老婆饼，曾经写下：“茶食泰斗”四个大字，以赞扬老婆饼。
图4 CIE-xy色度图与FUI水色指数21个色度坐标划分示意图(色度坐标数据来自于文献[13])Fig.4 Schematic diagram of 21 chromaticity coordinates of CIE-xy chromaticity diagram and FUI water color index
的色度坐标为(0.333 3,0.333 3)，表示三原色的等量混合，称为等能白光点。将等能白光点视为坐标原点，然后计算等能白光点与P(xM,yM)间的向量和x轴(y-yW=0处)正方向之间的角度αM,计算公式为：
αM=arctan(yM-yW,xM-xW)modulus2π 。
(7)
式中，αM从x正轴方向逆时针转动逐渐增大。将计算后的αM与表1中的αi对比，可得出P点的FUI数值。将计算所得的αM与表1中的αi对比，当αM大于αi时，色度角αM重分类得到的FUI数值为i。
表1 FUI指数查找表Tab.1 FUI Index lookup table
由于遥感图像本身波段离散的特点，与高光谱波谱的积分结果相比必然会给色度角αM的计算带来偏差。该偏差为高光谱积分的色度角αhyper与多光谱波谱积分得到的色度角αmulti之差Δα，即
Δα=αhyper-αmulti。
(8)
可建立Δα与αmulti的函数关系得多一个多项式公式，通过多光谱图像计算得到的色度角αmulti后加上系统偏差Δα可以达到消除偏差的效果。
2.2.2 基于离水反射率的FUI水色指数提取
水体的离水反射率Rrs包含了水体表面向上辐射的信息，可以反映水体种的物质组成成分，是水色遥感中最常用的表观光学参量。基于Rrs的水色提取中，将太阳光照视为一个常量，因此可以忽略太阳光照变化对水体颜色的影响。计算水体颜色的关键在于获取水体离水反射率Rrs，但是目前大多数卫星传感器在可见光范围内只有少数几个离散波段。因此，针对不同传感器需要进行不同的波段设置来计算色度角αM和FUI水色指数。
若要通过离散波段的遥感数据计算水色指数，首先要对缺失的波段进行插值。Woerd等[28]通过提取线性加权插值方法来解决这一问题。针对Sentinel-2和Landsat8 OLI的提取公式分别为：
(9)
(10)
分别以式 (9)和(10)计算Sentinel-2 与Landsat8 OLI的三刺激值XYZ，进而求得色度角αM，如Wored等所述，Landsat8 OLI与Sentinel-2所计算的色度角由于光谱带偏移，会与高光谱计算的色度角有介于-5～20的偏移量[29]。可以使用五阶多项式来校正αM，a为αM除以100，五阶多项式的系数如表2所示。
表2 基于波段线性插值的Sentinel-2,Landsat8 OLI 传感器色度角αM的Δα偏差校正多项式系数(a=αM/100)Tab.2 Sentinel-2,Landsat8 OLI sensor chromaticity angle deviation correction polynomial coefficients based on linear interpolation of the band (a=αM/100)
2.3 水体叶绿素a浓度与TSI指数
本文使用Hydrolab MS5(USA,HACH)多功能水质探头分别测得2016年、2017年与2019年官厅水库的叶绿素a浓度，剔除部分异常数据后，共有136组叶绿素a浓度数据。
基于叶绿素 a浓度的 TSI指数模型最为可靠，计算公式如下:
(11)
式中Chla为叶绿素a浓度。
通过计算得到的TSI指数，将水体营养状态分为3级，当TSI<30时，水体为贫营养；30
3 结果与分析
3.1 基于FUI的水体营养状态评价
本文将TSI指数分别与色度角α与FUI水色指数建模，并分别计算数据集间的Pearson相关系数。Pearson相关系数表示两组数据间的相关关系强弱，相关系数绝对值越大，则相关性越强[30]。计算公式如下：
(12)
从图5可以看到，官厅水库TSI数值介于18～41之间，大部分水体处于贫营养与中营养状态。随着TSI的增长，色度角α随之减小，FUI水色指数随之增大。即水体营养状态指数越高，水体颜色越偏黄。二者之间强相关。TSI与α间的拟合模型相关性较强，决定系数R2=0.72，Pearson相关系数为 -0.85。TSI与FUI间的拟合模型决定系数R2=0.693，Pearson相关系数为0.799。
(a)TSI指数与色度角α的拟合 (b)TSI指数与FUI指数的拟合
图6为基于 IOCCG模拟数据集[31]的 FUI 与 TSI 散点图。IOCCG数据集包括 500 条模拟数据，数据集中同时包含水体固有光学量(inherent optical properties,IOPs)和表观光学量(apparent optical properties，AOPs)数据。500条数据中叶绿素浓度范围为 0.05到 30 mg/m3，平均浓度为 6.08 mg/m3。
图6 IOCCG模拟数据集(N=500)中的FUI与TSI之间的关系[32]Fig.6 Relationships between FUI and TSI from the IOCCG simulated dataset (N=500)
从图5和图6中可以看到，当FUI水色从蓝色变为黄色时，水体富营养状态逐渐由贫营养向富营养过渡。尽管实地测量的数据中缺乏富营养水域，但是实地测量的数据集和IOCCG模拟数据集呈现了大致相似的总体趋势。总结两个数据集中FUI与TSI关系规律，可以初步通过FUI分段对水体营养状态进行划分，即当FUI≤6时，TSI<30，水体贫营养；当7≤FUI≤9时，30
由于FUI水色指数与水体光学组分(如叶绿素a和透明度)有着密切的联系，可以通过多源遥感卫星影像获得，因此FUI水色指数具有进行水体营养状态监测的潜力。然而，影响水体营养状态的不仅仅是水体光学成分，部分无色物质如总磷和总氮均对水体综合营养状态有一定的影响。因此，对于氮磷污染较为严重的水域中，仅以FUI水色指数对水体进行营养状态评估是不够的，还需要结合实地采样测量的总磷和总氮等物质浓度综合计算。大多数情况下，传统野外站点监测耗时费力，得到连续常规水质监测的湖库只占很少一部分。针对这种情况，FUI水色指数可以提供大面积区域水体的长时序监测，与传统监测站点优势互补，获取更丰富的水质信息，为决策部门提供数据支撑。
3.2 基于Sentinel-2的FUI水色指数提取结果
使用20景Sentinel-2影像对官厅水库2016—2020年不同季节水色进行反演，结果如图7和图8所示。可以看到，官厅水库水色在空间上有不均匀性，不同区域水色有着明显不同。总的来说，水库中心处FUI数值较低，水体呈现蓝绿色。而在水库边缘处由于人类活动较为频繁，有大量船只往来，水体较为浑浊，水色上就表现得比较偏向于黄棕色，这一点也和作者现场观测的结论一致。
(a)2016-01-04 (b)2016-04-13 (c)2016-07-02 (d)2016-10-10
(i)2018-01-13 (j)2018-04-18 (k)2018-08-01 (l)2018-10-30
图8 官厅水库各季度FUI平均数值Fig.8 The average FUI value of Guanting Reservoir in each quarter
从图7和图8来看，官厅水库FUI水色明显受到季节变化的影响。在每年的1月和10月，由于水位下落，水库水体面积总体减小，浮游植物死亡，水体颜色偏向于黄棕色，FUI数值较高。2020年4月官厅水库FUI数值为6.61，是近几年来FUI数值最低的。近3 a的柱状图呈现出1月FUI数值较高，4月大幅度下降，7月和10月小幅度上升。总体来说，近3 a的FUI数值比前2 a低，这一点得益于北京市政府对官厅水库的重视和有效治理。
3.3 Sentinel-2与Landsat8 反演结果对比
为了验证FUI水色指数在不同传感器之间的可比性，在不同区域的广泛适用性。挑选Landsat8与Sentinel-2同一天过境且云量较小的数据，通过筛选，选择2019年10月20日的数据作为Sentinel-2与Landsat8的对比数据。对两景影像分别进行FUI水色反演，反演结果如图9所示。基于Sentinel-2的反演结果加权平均得到数值为13.04，基于Landsat8 OLI的结果为13.16，二者结果很接近。总的来说，FUI水色指数在不同传感器影像之间的反演所得数据偏差不大，说明FUI在不同传感器之间具有可比性。基于此特性，FUI水色指数的提取可以使用多源传感器来弥补单一传感器数据覆盖范围偏小、重返周期较长的不足。
(a)基于Sentinel-2的FUI反演 (b)基于Landsat8的FUI反演图9 Sentinel-2与Landsat8 OLI反演FUI指数结果对比Fig.9 Sentinel-2 and Landsat8 OLI inversion FUI index results comparison
4 结论与展望
本文介绍了基于CIE-XYZ色彩空间的水体颜色参量色度角α和FUI水色指数的遥感提取算法。分别使用Landsat8 OLI与Sentinel-2影像提取官厅水库FUI水色，并利用色度角α对提取结果进行校正。提取了官厅水库2016—2020年各季相的水色数据。从时间尺度来看，官厅水库FUI水色受季节影响较大，总体呈现出1月最高、4月最低的趋势，FUI平均数值逐年降低。从空间尺度来看，水库中心区域FUI数值较低，而边缘区域数值较高。挑选两景同一天过境的Landsat8 OLI和Sentinel-2影像进行FUI水色提取，反演结果十分接近。证明FUI水色指数在不同传感器间具有可比性。FUI水色指数提取可以使用多源传感器来弥补单一传感器在时空分辨率上的不足。
基于原位测量数据和IOCCG模拟数据建立了FUI-TSI相关关系，构建了基于FUI的水体富营养状态评估模型。由于FUI水色指数与水体光学参量间具有强相关性，因此FUI水色指数具有评价水体营养状态的潜力，可与传统野外站点监测数据结合，综合评价水体富营养状态。
FUI水色数据还有更大的潜力有待挖掘。后续可以深入研究FUI与各项水质参数的关系，从水体光谱的机理着手，展开FUI水色与水体透明度、悬浮物浓度等参数的研究。同时，随着高分系列卫星的发射升空，可以将FUI水色指数的提取算法移植到新的传感器上，以满足更多分辨率层次的水体研究。
