基于色度角的岱海水体营养状态遥感评估*

王 芳，青 松，刘 楠，郝艳玲，包玉海

(1:内蒙古师范大学地理科学学院，呼和浩特 010022)(2:内蒙古大学生态与环境学院，呼和浩特 010021)

湖泊富营养化已成为内陆湖泊面临的最紧迫的水环境问题之一[1-2].湖泊富营养化是在自然因素和人类活动双重影响下使得氮磷等营养盐浓度增加，导致湖泊生产力和生物量过度提高，溶解氧含量和湖泊透明度降低，严重时产生毒素和难忍气味，从而威胁到水生生态系统的可持续性以及人类健康[3].因此，为水体富营养化的防治提供科学依据，加强湖泊富营养化方面的方法研究尤为重要.

“营养状态”是一个与生态系统功能和人类活动对水质影响密切相关的属性[4].自1960s以来，人们尝试用单变量营养指数或多参数方法定量评价湖泊的营养状态[5].Carlson提出了温带湖泊营养状态指数，该指数可以在叶绿素、透明度和总氮的对数尺度上表征营养状态，用以评价湖泊的养分富集状况和初级生产力.其中，水体透明度与悬浮藻类的生物量和生产能力密切相关[6]，是直接反映水上植物初级生产力和水中微生物活动的一个重要参数[7].透明度是湖泊管理机构用于水质评价的一种简单而广泛使用的工具[8]，与湖泊水体光谱辐射观测密切相关[9].

湖泊富营养化遥感评价是通过分析水体反射、吸收和散射光谱特征与水质参数浓度之间的关系，建立富营养化水质参数的定量遥感反演模型[10].随着内陆水色遥感研究的发展，水质参数的监测从定性发展到定量，其中伴随着多光谱遥感数据应用于水质监测而发展起来一系列经验算法.这些方法是基于遥感数据和地面实测数据，选择最优波段或波段组合，与水质参数建立经验关系，进而评价水体营养状态[11-14].

水体颜色作为最古老的时间序列水质数据之一，是水质调查中的一项重要内容.关于水体颜色的早期研究主要是试图建立颜色与水质参数之间的关系.Alföldi等[15]采用Landsat MSS假彩色数据作为CIE标准坐标系统中的X、Y、Z三刺激值，分析了水体颜色与水体状况的关系.自此，随着后续学者的研究和逐步完善，构建了完整的FUI水色指数[16].FUI水色指数是基于光学特性的水色分类，取值为1～21的离散值.Wang等[17]结合实测数据和MODIS产品，建立了FUI与营养状态指数之间的关系，提出了一种评估全球内陆水体营养状态的方法.然而，根据Pitarch等[18]的研究指出，与色度角alpha相比，FUI作为大多数营养海洋水域的监测工具是粗糙的.色度角alpha是CIE色度系统中计算FUI的过程参量，取0°～360°之间的连续值[19]，是自然水体的一个简单客观参数，被称为自然水域的真实颜色[20].van der Woerd等[21]证明了由色度角alpha表示的颜色能够从海洋水色卫星仪器中准确获取，并指出色度角在Sentinel_2 MSI、Landsat_8 OLI和Landsat_7 ETM+等传感器上有很大的推广应用潜力.Xu等[22]基于色度角，揭示了洞里萨湖19年间水体颜色变化趋势及规律，并分析了色度角变化率的空间分布和富营养化状态.然而，到目前为止还没有相应的文献是由连续变量alpha与TSI指数之间的相互关系来评估湖泊营养状态.此外，营养状态指数已广泛应用于温带、热带和亚热带地区湖泊的营养状态分析[23]，但在我国北方半干旱地区湖泊水体中的研究较少.

因此，本文试图探究色度角alpha和TSI之间的相关性，建立适合于北方典型湖泊——岱海的营养状态评估模型.本文以内蒙古岱海为研究区，进行如下研究：(1)基于实测数据建立岱海水体营养状态评估模型；(2)将评估模型应用于Landsat TMETM+OLI和Sentinel_2A MSI遥感数据；(3)分析岱海水体营养状态的时空分布特征及影响因素.

1 数据

1.1 研究区

岱海盆地在新构造断块活动时开始出现，是典型的地堑式断陷盆地，早更新世初开始形成为湖泊[24]，现称其为岱海(40°29′07″～40°37′06″N，112°33′31″～113°46′40″E)，位于内蒙古乌兰察布市凉城县东南部，四周环山，是一个典型的封闭性内陆咸水湖(图1).湖泊水源补给方式主要为降水补给、地表径流和地下补给等[25].岱海地处我国半干旱与半湿润的过渡带，即北方环境敏感地带，属于中温带半干旱大陆性季风气候，冬季长而干冷，夏季短而炎热[26].岱海，是内蒙古的第三大内陆湖，与呼伦湖、乌梁素海统称为 “一湖两海”，在调节气候、生态修复、涵养水源方面起着重要作用，是我国北方生态安全屏障的重要组成部分.岱海是被列入《中国湿地保护行动计划》的第179块国家重要湿地之一，同时也是自治区级湖泊湿地自然保护区[27].多年来，由于降雨量减少、地下水过度开采和环境污染多种因素的叠加，岱海水体水面急速缩减，导致水体富营养化问题日益突出[28-29].

1.2 实测数据

1.2.1 遥感反射率 在2019年10月9日、2020年6月14日、7月19日、8月24日、9月17日和10月12日，在岱海水域共开展了6次野外调查，进行水面光谱测量，共获取了216组遥感反射率(Rrs)光谱数据.Rrs测量采用“水上法”[30]，每个采样点需要测量标准灰板、水体和天空的辐射亮度.本研究中采用美国ASD公司生产的ASD Field Spec 4 地物光谱仪，测量水体光谱.遥感反射率计算公式如下，文中采用400～900 nm波段处的遥感反射率.

(1)

式中，LW、LSky、Lp分别为测量的水体、天空和标准灰板的辐射亮度；r为气-水界面的反射率；ρp为标准灰板的反射率.根据上述处理过程得到216个采样点的遥感反射率光谱曲线，如图2.由于水体中各种物质的影响，水体光谱曲线表现出明显的峰值.水体实测反射率光谱曲线在蓝波段处的反射率较小，无明显峰谷，在波段550 nm处由于悬浮物浓度的影响出现高吸收峰值，在700 nm左右的波段出现较小峰值.

图2 岱海水体实测遥感反射率光谱曲线

1.2.2 透明度 在岱海水域开展的6次野外调查中，与遥感反射率同步采集了216组透明度(SDD)数据.透明度是以塞氏圆盘现场测定，将塞氏圆盘放入水中，以消失在视野中时的下沉深度作为观测点的透明度值[31].本文将透明度作为基本参数计算了水体营养状态指数TSI(详见2.1节).

1.3 遥感数据

本文利用Sentinel_2 MSI和Landsat_5 TM、Landsat_7 ETM+、Landsat_8 OLI遥感数据进行营养状态反演.MSI Level-2A数据是由Sen2cor软件得到.OLI-SR 数据是基于 6S 大气辐射传输模型处理获得[32].TM/ETM+的SR数据由LEDAPS算法生成[33].使用ENVI IDL进行批量裁剪及重采样处理，随后利用MATLAB软件进行NDVI水体提取和遥感反演.

本文获取到1986-2020年间(缺少1988年遥感影像)岱海区域遥感数据，包括1986-1999年间Landsat_5 TM影像共69景，2000-2012年间Landsat_7 ETM+影像共66景，2013-2020年间Landsat_8 OLI影像共47景，2019-2020年间Sentinel_2 MSI 影像共29景.采用gapfill插件进行了ETM+影像条带修补.Landsat系列数据均来自美国地质调查局官网(https://earthexplorer.usgs.gov/)的地表反射率产品，Sentinel_2数据来自欧空局官网(https://scihub.copernicus.eu/)的MSI Level-2A 级地表反射率产品.

此外，从中国气象信息中心(http://data.cma.cn)下载了1986-2020 年间凉城气象站的气温、降水和风速等数据，用于研究湖泊营养状态与气象因素之间的关系.由于岱海湖泊存在结冰期，对气象数据统计采用4-10月份的月平均值以及年平均值.

2 方法

2.1 营养状态指数

卡尔森(1977年)提出了TSI与透明度、叶绿素和总磷等参数相关联的理论解释，开发了基于透明度为核心的营养状态指数.透明度易于获得，是水生生态系统营养状态的综合代理，由水体中的悬浮物、藻类和溶解有机物控制，通过控制水下光场而影响水生植物和生物的初级生产能力[34].本文采用透明度单个水质指标来计算营养状态指数[35-36]，计算公式如下：

(2)

式中，SDD为透明度，ln2为2的自然对数.根据我国环境监测总站颁布的《湖泊(水库)富营养化评价方法及分级技术规定》，可将湖泊营养状态分为3类，TSI<30为贫营养状态、30≤TSI≤50为中营养状态、TSI≥50为富营养状态.富营养状态又可分为，轻度富营养化(5070).

2.2 色度角alpha

与FUI相比，色度角alpha是一个连续的数值变量，用一个数字表示光谱的颜色[22].色度角是(x，y)坐标的极角，通过将3个可见光波段的Rrs转换到国际照明委员会(CIE)色彩空间中的坐标(X，Y，Z)，再将其转化为CIE色度系统内的色度坐标(x，y，z).为了清晰量化CIE系统中的颜色，将坐标原点旋转到(1/3，1/3)处，则每个(x，y)坐标转换为新的坐标(x′，y′)[37-39].CIE-RGB分别选取为665、560、490 nm波段.

CIE-RGB与CIE-XYZ之间的转换公式如下：

X=2.7689R+1.7517G+1.1302B

Y=1.0000R+4.5907G+0.0601B

Z=0.0000R+0.0565G+5.5934B

(3)

CIE色度坐标(x，y，z)是从X、Y、Z通过归一化到0～1之间计算获得.又基于色度坐标(x，y)，建立新的坐标系(x′，y′).公式如下：

x′=y-1/3

y′=x-1/3

(4)

基于新坐标，色度角计算公式如下：

alpha=arctan2(x′,y′)180/π

(5)

式中，arctan2函数表示双变量反正切函数，值域为(0°，360°)，x′和y′表示新坐标系中的坐标.

2.3 经验算法建立

色度角是一个定义明确的物理量，它直接与光学深水的组成和固有光学性质有关[22]，可以表示水体颜色，且不同混浊和营养状态的水在色阶上的呈现不同[37].因此，本文通过最小二乘回归试图建立TSI与alpha之间的经验关系.

TSI=a·alpha+b

(6)

式中，a和b为模型系数.

2.4 精度检验

利用均方根误差(RMSE)、平均百分比误差(MAPE)评价反演模型精度，公式分别为：

(7)

(8)

式中，ximea和xipre分别为模型估计TSI和实测估计TSISDD，N为采样点数.

2.5 数据分析

为了探究岱海水体营养化的主要气候变量，首先利用SPSS软件对数据进行了Z-score标准化，随后使用Amos 24.0软件创建了一个结构方程模型(SEM)，量化了气温、降雨量和风速对岱海水体营养化的影响，并利用绝对拟合值数：拟合优度指数(GFI)、近似均方根误差(RMSEA)和相对拟合值数：比较拟合指标(CFI)进行模型拟合度检验.则GFI和CFI大于0.9，越接近1越好，RMSEA小于0.05时模型拟合良好.

3 结果与分析

3.1 营养状态遥感反演模型

本文试图寻找基于遥感反射率R、G、B波段运算的alpha与TSISDD之间的关系.将alpha作为自变量，TSISDD视为因变量(图3a)，利用随机抽样选取的162组实测透明度和遥感反射率数据，采用最小二乘法确定模型参数，建立岱海湖营养状态遥感评估模型，公式如下：

图3 基于实测数据的TSISDD与色度角alpha散点图(a); 模型检验(b)

TSI=0.65alpha-45.78

(9)

利用剩余的54组数据，检验模型精度(图3b)，可以看出反演TSI与实测TSISDD较均匀分布在1∶1线两侧.根据检验结果显示，R2=0.74，MAPE=4.84%，RMSE=3.66，表明本文建立的经验模型精度较高，可以很好地反映岱海的营养状态.将模型应用于216组实测数据得到岱海营养状态，如表1所示.岱海共216个采样点，其中3个采样点为中营养化，213个采样点为富营养化(其中，属于轻度富营养化的站点最多，为96)，无贫营养化采样点.

表1 基于实测数据的营养状态评估模型评价结果

3.2 营养状态评估模型的遥感应用

3.2.1 Sentinel_2 MSI反射率验证 利用时空匹配的实测遥感反射率数据，验证了 MSI 和 OLI 数据的大气校正结果，见图4.从图4a、b、c和e可以看出，2019年10月9日(R2=0.82，RMSE=0.01 sr-1，MAPE=31.70%)、2020年6月14日(R2=0.92，RMSE=0.01 sr-1，MAPE=49.84%)、2020年7月19日(R2=0.77，RMSE=0.01 sr-1，MAPE=27.37%)和2020年10月12日(R2=0.72，RMSE=0.01 sr-1，MAPE=42.42%)的实测与Sentinel_2 MSI光谱反射率具有较好的一致性.图4d为2020年8月24日实测与Landsat_8 OLI光谱反射率之间的R2=0.57，RMSE=0.01 sr-1，MAPE=31.0%.这些结果表明了MSI和OLI数据大气校正结果的可靠性.

3.2.2 Sentinel_2 MSI反演结果验证 根据同步野外实测日期的时空匹配遥感影像，利用2019年10月8日、2020年6月14日和2020年10月12日的MSI数据，验证了alpha和TSI反演结果(图5).从图5a中可以看出，实测与同一天MSI数据反演得到的色度角alpha具有较好的一致性，R2=0.69，MAPE=2.51%，RMSE=5.43°.图5b中，基于本文经验模型反演得到的TSI与实测TSISDD的验证结果为，R2=0.42，MAPE=4.87%，RMSE=6.31.这些结果表明了基于MSI数据的alpha和TSI反演结果的可靠性.

3.2.3 Landsat _8 OLI与Sentinel_2 MSI反演结果对比 利用 2020 年 8 月 24 日数据检验了OLI 数据反演得到的alpha和TSI结果.从图6可以看出精度检验看出，alpha和TSI都得到了较好的反演结果，检验精度分别为R2=0.84，RMSE=9.28°，MAPE=4.65%和R2=0.85，RMSE=5.76，MAPE=7.64%，说明反演效果较好.

图6 实测alpha(a)、TSI(b)与Landsat_8 OLI遥感反演对比

利用 2019 年 4 月 16 日的 MSI 与 OLI数据进行交叉验证，alpha验证结果如图7a，为R2=0.86，RMSE=1.60°，MAPE=0.70%.TSI验证结果如图7b，为R2=0.86，RMSE=1.04，MAPE=1.18%.从精度检验结果来看 MSI 和 OLI 数据在反演水体营养状态方面有很好的一致性.

图7 Sentinel_2 MSI 与 Landsat_8 OLI 的交叉验证

根据本文建立的经验模型，利用2019年10月8日、2020年6月14日、7月21日、8月24日和10月12日的遥感影像进行岱海营养状态遥感反演(图8).从图中可以看出，在6月14日水体南部富营养现象较严重，7月21日水体东部富营养化现象增加南部稍缓，且8月24日岱海西南部富营养化现象增加，水质整体下降，到10月12日水体富营养化逐渐缓和.2019年10月相较于2020年10月水质状况较好(图8a).如图8b岱海富营养化分级图中可以看出，2019年10月8日、6月14日以及7月21日岱海主要处于中营养状态，8月24日和10月12日处于轻度富养状态，但 8月24日边缘部分富营养化较严重，7月21日相较于其他日期水质状况较好.总体上，岱海水体水陆边界附近营养化存在较严重现象，且逐渐往湖中心趋于缓和.

图8 岱海营养状态遥感反演(a)和富营养化分级(b)

3.2.4 岱海水体营养状态年际变化 本文利用时间序列Landsat_5 TM、Landsat_7 ETM+、Landsat_8 OLI和Sentinel_2A遥感影像，用上述经验模型进行营养状态反演，并除以给定年份的图像总数得到1986-2020年岱海营养状态年平均空间分布图(图9)和岱海富营养化分级空间分布图(图10).

图9 岱海营养状态年平均空间分布

从1986-2020年岱海营养状态年度动态中(图9)可以看出，1986-1993年岱海湖富营养化现象逐渐严重，主要在南部以及西南部富营养化较严重，1992年和1993年湖泊污染最严重，湖泊整体呈现出中度富营养化.从1994年开始富营养化现象逐渐好转，2002年中营养化现象最为明显.到2006年轻度富营养化现象又逐渐加重，2007年南部离岸附近富营养化现象进一步加重.2008年水体中营养化所占区域增大，到2009年又进入轻度富营养化状态，中营养化区域减少，且南部边缘区营养化稍加严重.2010年和2011年水体又转变为中营养化，尤其在2010年出现贫营养化状态，为水质最好时期.2012年又恢复到富营养化状态，且在2013年以及2014年营养化现象好转，中营养化区域逐渐增加但从2015年以来水体富营养化又加重，在2016年开始水体轻度富营养化开始出现在湖岸附近，水体整体以重度富营养化为主，到2018年逐渐转变中度富营养化.2019年水体营养状况好转，但在2020年又进入轻度富营养化.总体来讲，从1986-2020年之间，岱海水体营养化现象大体上从边缘逐渐向湖中心趋于缓和，重度富营养化现象主要出现水陆边界线附近，离岸边越近富营养化现象越严重，趋向湖中心则以中营养为主，可得出岱海水体长期处于轻度富营养化状态，并在中度富营养化或中营养化之间转换(图10).

图10 岱海富营养化分级

3.2.5 岱海水体营养状态月变化 从图11明显可以看出，4月和5月水体富营养化现象最为严重，其次为10月，主要出现在湖泊边缘区域，湖中心则以中度富营养化为主.此外，6月虽然是以中度富营养化为主，但湖中心区域出现中营养化，7月和9月水质情况相同，水体以中营养为主，边缘区域则为轻度富营养化.8月为水质状况最好时期，主要以中营养化为主，湖中心出现贫营养化现象，且湖岸附近明显为轻度富营养化.大致与岱海营养状态年度动态变化相同，湖泊富营养化严重部分主要在水路边界附近，趋向湖中心营养化等级一般则降一个度，水质营养状况最好区域一般出现在湖中心部分.

图11 岱海营养状态月平均空间分布(a)和富营养化分级(b)

3.3 岱海营养状态影响因素分析

为说明岱海湖泊营养状态变化驱动因素，收集了气候因素数据，分别为气温、风速、降水量等.在湖泊的演变过程中，气候起着十分重要的作用，尤其对我国干旱半干旱区域的内陆湖泊，气候变化对环境的影响更为强烈和直接[40].从图12可以看出，1986-2020年凉城县气象站年平均气温呈现上升趋势，降雨量和风速整体无明显趋势，但降水量在相邻年间的波动较大.总体上，气温、风速、降水量等与年平均营养状态无显著相关关系，但在部分时间段上呈现出正相关或负相关关系.利用1986-2020年凉城县气象站获取月平均降水量、月平均气温和月平均风速(图13)用以分析岱海营养状态.可以看出，TSI与风速趋势大致相同，都是6-8月逐渐下降，再从8-10月上升，8月为最低点.TSI与降水量趋势大致相反，与气温之间的相关性不太显著.

图12 营养状态与年平均气象因素趋势

图13 营养状态与月平均气象因素趋势

采用结构方程模型(SEM)方法分析了气象因素对TSI的影响.模型拟合度为χ2=2.252、df=2、P=0.324、GFI=0.965、CFI=0.914、近似均方根误差(RMSEA)=0.062和χ2=9.614、df=2、P=0.008、GFI=0.720、CFI=0.681 和RSEMA=0.797(图14).结构方程模型表明，潜在变量如气候因子对岱海水体年际营养状态变化的解释率为0.13%，气温和风速对其间接影响属于负相关，相关系数分别为-0.26和-0.14，降水量对其间接影响属于正相关，相关系数为0.16.对于月平均营养状态变化的解释率为93%，其中与气温和降雨量为负相关，风速为正相关，间接影响可达93%，与图13展现的相关性一致.气候因子对岱海水体年际营养化影响并不显著，但与月平均营养化的相关性较高.气温对年际变化或月平均TSI的影响都表现为负相关，岱海流域气温呈现出逐年上升趋势，水体温度的升高，影响了水体浮游植物的光合作用[41].风速无明显年变化性，但月变化明显.风速过大，引起湖底泥沙搅动，沉积悬浮物升起使得水体透明度下降，从而影响水体营养状态.此外，另有研究表明，风浪扰动引起悬浮的营养盐形态主要以颗粒态为主，而能被藻类直接利用的溶解浓度变化不大，并且水柱中营养盐浓度随着风速减弱、颗粒物沉降而逐渐降低[42].浅层富营养化湖泊对风速下降极其敏感，随着风速降低，尤其是低风速持续时间延长，湖底间歇性缺氧/厌氧的概率增加，会导致浅水湖泊中底泥溶解性养分的高释放，从而加重水体富营养化[43].降雨量对年平均TSI的影响并不明显，是由于岱海地区降水量少而集中，主要体现在年平均降水量较少，从而对TSI的影响微小呈现负相关；相反，对于月平均降水量主要表现为集中的特点，从而对TSI的影响较大呈现显著正相关关系.

图14 营养状态影响因素结构方程模型(年平均数据和月平均数据)(矩形表示测量变量，而圆圈表示潜在变量.与路径箭头相关联的值表示标准化的路径系数)

4 结论与讨论

本文得到的营养状态反演结果存在一定的不确定性，主要体现在实测数据与遥感数据时间和空间尺度上的不一致，并且有些年份中月影像数据的缺失，以及由于2003年ETM+机载扫描行校正器发生故障，导致之后收集的数据缺失出现条带状，等现象均可能会对结果产生一定的影响.本文研究中采用透明度作为参数，与张运林等的研究相比指标参数过于单一.在进一步的研究中，需要结合其他遥感数据以及参数，需对alpha展开更为深入的研究.

本文利用实测光谱和透明度数据，基于营养状态指数法和色度角之间建立了适于岱海营养状态的经验模型，并对模型进行了精度检验.随后，将算法应用到Sentinel_2 MSI、Landsat_5 TM、Landsat_7 ETM+和Landsat_8 OLI遥感影像，分析了1986-2020年间岱海水体营养状态时空分布特征，并使用结构方程模型分析了水体富营养化的气候影响因素，具体结论如下：

1)利用Sentinel_2 MSI和Landsat_8 OLI数据，基于实测TSI数据与色度角，建立了营养状态评估模型，用于预测北方内陆湖泊-岱海水体的营养状态.采用交叉验证方法进行精度验证，检验精度结果较好(R2=0.70，MAPE=7.4%，RMSE=5.20).

2)将该经验模型应用到Sentinel_2 MSI和Landsat_8 OLI卫星数据，反演水体营养状态.利用Sentinel_2 MSI反演的TSI进行验证，结果表明遥感反演结果较好，TSI的检验精度为R2=0.42，MAPE=4.87%，RMSE=6.31；alpha的检验精度为R2=0.69，MAPE=2.51%，RMSE=5.43°.

3)通过时间序列MSI、TM、ETM+和OLI数据得到了岱海1986-2020年的时空分布图.从1986-2020年之间，岱海水体面积逐年减少，且多数时间处在轻度富营养化状态，并与中营养化或中度富营养化之间切换，水体富营养化现象大体上从边缘逐渐向湖中心趋于缓和，离岸边越近富营养化现象越严重，趋向湖中心则以中营养为主，贫营养化现象极少.

4)岱海的营养状态时空变化与气温、风速和降水量等相关性并不显著，气候因子对其解释率为13%.营养状态月变化与气候因子影响显著，对其解释率为93%.