高邮湖位于江苏省中部,是我国第六大淡水湖,具有渔业生产、供水、蓄水、泄洪排涝等功能,对当地的生态环境和经济发展都有重要影响。自2021年以来,高邮湖菹草大面积消失,水质参数如总磷、总氮、叶绿素a(Chl.a)浓度都有增长趋势,富营养化威胁日益严重。因此,能够快速、周期性地监测高邮湖水质状况,评价富营养化程度显得尤为重要。水体叶绿素a浓度是衡量湖泊富营养化程度的重要参数,总悬浮物(TSM)浓度分布与水华分布有一定的相关性,因此也是湖泊水质监测的重要参数[1]。相比于耗时长、监测代表范围有限的传统水质监测方法,卫星遥感监测技术具有范围广、周期短、成本低的特点。欧洲航天局发射的哨兵-2号卫星有A、B双星,重访周期5 d,其上搭载的多光谱成像仪(MSI)在可见光-近红外波段的空间分辨率为10~20 m,适合高邮湖这类中小型湖泊的水质监测对时间和空间分辨率的需求。
目前,针对水质参数遥感的模型估算方法分为经验法、半经验法、分析法和半分析法等。其中经验法仅对Ⅰ类水体有效,而内陆浑浊的Ⅱ类水体,其主导光学因子往往是由水体叶绿素a、总悬浮物颗粒、黄色物质等共同主导的,单纯用经验法反演普适性不够[2]。而利用半经验法建立的水质参数反演方法在当年表现良好,却难以在不同年份推行[3]。分析法需要大量准同步气象水体固有光学特性等数据,较难获取。相比之下,半分析法是目前应用较为广泛的方法。针对Chl.a常用的半分析法模型包括两波段法[4]、三波段法[5]、归一化叶绿素指数(NDCI)法[6]、基线荧光高度(FLH)法[7]等。针对TSM,前人研究发现,光谱波长范围为700~850 nm时,水体中色素颗粒物和非色素颗粒物吸收度都较弱,此时单波段反射率与TSM相关性较好[2, 8-10]。
本研究对高邮湖进行地面光谱和水样实测,从几种常用的半分析模型中提取出可应用于哨兵-2 MSI数据运算的波段组合,与实测值模拟出可适用于哨兵-2 MSI卫星影像数据的高邮湖Chl.a和TSM反演模型,以期对高邮湖的富营养化程度进行可持续、周期性、低成本的监测。
1 研究区概况和数据处理 1.1 研究区概况高邮湖总面积为760.67 km2(其中水域面积为648 km2,苇滩和堤坝面积为112.67 km2),其西北方有淮河入江水道贯穿而过,是典型的过水性湖泊。每年淮河的入湖水量占高邮湖总入湖水量的95%。高邮湖南接邵伯湖,北部原有大汕子河与宝应湖连接,后上方水闸关闭,形成半封闭水域,有大面积围网养殖现象。受夏季季风和长江淮河流域汛期影响,在丰水期(5—10月),高邮湖最高水位>5.9 m,在平水期(11月至次年4月),高邮湖最低水位<5.3 m[11]。
1.2 数据采集及处理 1.2.1 地面实测数据在2019年3月和9月,2022年4月和10月对高邮湖进行4次采样,其中2019年采样点位24个,2022年采样点位64个,共获得实测数据176组,研究区及采样点位分布见图 1。剔除部分水样Chl.a低于标准检出限的采样点,以及部分湖底水草影响光谱测量的采样点后,可用于模型模拟和验证的采样点共160组。
现场测量水深、风速和风向,同时采集水样放入冷藏箱保存,带回实验室分析。根据《水质叶绿素a的测定分光光度法》(HJ 897—2017)[12]分析得到水体Chl.a质量浓度,根据《水质悬浮物的测定重量法》(GB 11901—89)[13]分析得到水体TSM质量浓度。
采用水面以上光谱测量法,测量采样点的水体反射光谱曲线[14]。使用手持式地物光谱仪(ASDFieldSpec© HandHeldTM 2,HH 2 pro 325-1075,美国ASD公司),在船头位置分别测量采样点25%灰板(即平均反射率为25%的漫反射标准参照板)、水体和天空光的辐射光谱,通过公式(1)计算得到水面光谱反射曲线。
$ R_{\mathrm{rs}}=\left(S_{\mathrm{w}}-r S_{\mathrm{sky}}\right) \rho_{\mathrm{p}} / \pi S_{\mathrm{p}} $ | (1) |
式中:Rrs——采样点位水面离水辐射对应水面入射总辐射的遥感反射率;Sw、Ssky、Sp——仪器观测水体、天空光和灰板接收到的亮度值,即DN值;ρp——经严格定标后的灰板的反射率,采用的灰板为400~1 080 nm有标定的方向-半球反射比数据,可取平均值0.25;r——气-水表面反射率,一般情况下,在平静水面可取值0.022,当风速为5 m/s左右时可取值0.025,当风速为10 m/s左右时可取值0.026~0.028[15]。为减少仪器噪声造成的毛刺误差,选择Savitzky-Golay平滑算法对反射率曲线进行平滑处理。
1.2.2 哨兵-2 MSI数据哨兵-2 MSI数据来自欧洲航天局官网(https://dataspace.copernicus.eu/),下载与实测数据日期最接近的准同步数据(2019年3月13日与9月29日,2022年4月6日与10月3日)。哨兵-2 MSI数据中L2A级数据是对原始影像完成了大气校正的数据,将L2A可见光-近红外波段数据与地面实测光谱进行对比,结果见图 2(以2022年10月3日影像为例)。由图 2可见,波段数据与实测光谱曲线趋势一致,表明其大气校正结果可以接受。利用官网提供的snap软件对L2A数据进行拼接裁剪后,再利用归一化水指数(NDWI) [16]提取高邮湖水体部分。
高邮湖实测ρ(Chl.a)和ρ(TSM)统计结果见图 3(a)(b)。由图 3可见,ρ(Chl.a)均值在丰水期较高,平水期较低,因为丰水期与藻类生长季时间重合。相同时期,2022年ρ(Chl.a)均值要高于2019年。而ρ(TSM)在不同年份差异较大。在2019年,丰水期均值高于平水期;在2022年,丰水期与平水期均值相当,与2019年丰水期均值接近,都高于2019年平水期均值。高邮湖属于淮河流域,根据扬州市水资源公报[17-18],扬州市属淮河流域的地区年平均降水量在2019年较多年平均降水量少23.3%,其中平水期的降水量占全年的40%左右。而在2022年较多年平均降水量多30.9%,其中平水期的降水量占全年的49%左右。由此可见,位于淮河流域的高邮湖在2022年平水期的降水量高于多年平均降水量,更远远高于2019年平水期的降水量。过多的降水会导致更多的浑浊地表径流经淮河流入高邮湖,导致2022年平水期的ρ(TSM)升高。由此可见,在不同年份、不同水文环境下,高邮湖各水质参数浓度分布情况也不同,因此反演模型在高邮湖不同年份的适应性非常重要。
高邮湖实测水体光谱反射曲线见图 4。由图 4可见,由于水体中Chl.a的反射和TSM的散射作用,在绿光波段(568 nm左右)出现1个反射峰,随后在红光波段(600~680 nm),由于纯水的吸收导致反射率逐渐降低而出现吸收谷,到近红外波段(704和810 nm左右),反射率曲线又再次出现峰值,这是Ⅱ类水体中含有叶绿素a和悬浮泥沙的典型特征。而哨兵-2 MSI数据的可见光-近红外波段基本上涵盖了上述水体在蓝光、绿光、红光以及近红外波段表现出的光谱特性,因此可以用来构建模型,反演水体ρ(Chl.a)和ρ(TSM)。
高邮湖水面典型采样点实测光谱反射曲线见图 5(a)(b)。由图 5可见,水面光谱反射曲线的反射峰、吸收谷的反射率大小主要受ρ(Chl.a)大小的影响,而水体的总反射率大小主要与TSM有关。
哨兵-2 MSI数据为多光谱遥感数据,实地测量的光谱反射曲线的光谱分辨率为1 nm,为减少两者差异,便于反演模型在卫星遥感数据中的应用,本研究结合哨兵-2 MSI的光谱响应函数,将实测光谱反射曲线加权相加,拟合成哨兵-2 MSI波段数据。
采用半分析法,提取出Chl.a和TSM的波段组合公式,建立波段组合与实测水质数据的相关关系,选出与实地测量值相关性较好的模型方法,再通过模拟一元回归方程和多元回归方程,对比验证模型精度,选出更优的反演模型。采用拟合度(R2)检验模型模拟的吻合程度,用均方根误差(RMSE)和平均相对误差(MRE)进行模型精度验证和比较,三者公式见式(2)—(4)。
$ R^2=\frac{\sum(\hat{y}-\bar{y})^2}{\sum(y-\bar{y})^2} $ | (2) |
$ \text { RMSE }=\sqrt{\frac{\sum(\hat{y}-y)^2}{n}} $ | (3) |
$ \text { MRE }=\frac{1}{n} \sum \frac{|\hat{y}-y|}{y} \times 100 \% $ | (4) |
式中:
采用半分析法提取哨兵-2 MSI波段运算方法。两波段法[2, 19-20]:一般利用蓝绿波段、红-近红外波段反射率进行比值或差值运算,参与模型模拟。波段运算公式为蓝绿波段比值RB3/RB2,红-近红外比值RB5/RB4和RB5/RB6,以及蓝绿波段差值(RB3-RB2),红-近红外波段差值(RB5-RB4)和(RB5-RB6)。三波段法[21-22]:3个波段需要分别选择在对浮游植物色素吸收最敏感的波段λ1,不受黄色物质、非色素颗粒物吸收和后向散射干扰的波段λ2,纯水吸收处于主导的近红外波段λ3,在哨兵-2 MSI数据中对应波段为B4、B5和B8,波段运算公式为[(RB4-1-RB5-1)×RB8]。
NDCI法[23]须涉及2个波段,分别位于665~675 nm处(因Chl.a吸收而成的吸收峰)和700 nm附近(对Chl.a变化敏感的反射峰)。在哨兵-2 MSI数据中,对应波段为B4和B5,波段运算公式为(RB5-RB4)/(RB5+RB4)。
经验组合:根据之前的研究,发现在高邮湖蓝绿波段与近红外宽波段差值的比值与Chl.a有较好的相关性,波段运算公式为(RB3-RB8)/(RB2-RB8)[3]。
建立波段组合与实测ρ(Chl.a)的相关关系,通过对相关性的分析发现,在高邮湖,不论平水期和丰水期,ρ(TSM)的变化会影响各波段组合与实测值的相关性[24]。不同ρ(TSM)下各波段组合与实测ρ(Chl.a)之间的相关关系见图 6。由图 6可见,波段组合RB5/RB4,RB5-RB4,三波段法和NDCI法在ρ(TSM)较高的环境中与实测ρ(Chl.a)的相关性较高,且受ρ(TSM)大小变化的影响较小。在ρ(TSM)较低的环境中,经验组合与实测ρ(Chl.a)的相关性较高。综合考虑,选择ρ(TSM)=20 mg/L为区分阈值,当ρ(TSM)<20 mg/L时,选择经验组合与实测值建立一元回归Chl.a反演模型,同时利用偏最小二乘法对经验组合,RB5-RB4和RB7建立多元回归Chl.a反演模型;当ρ(TSM)≥20 mg/L时,选择与实测ρ(Chl.a)相关性最高的三波段法建立一元回归Chl.a反演模型,同时利用偏最小二乘法综合RB5/RB4,三波段法和NDCI法建立多元回归Chl.a反演模型。Chl.a反演模型及精度见表 1。
实测数据符合ρ(TSM)<20 mg/L的模型共有28组,用于模型模拟的有19组(2019年8组,2022年11组),用于模型验证的有9组(2019年4组,2022年5组);实测数据符合ρ(TSM)≥20 mg/L的模型共有132组,用于模型模拟的共有89组(2019年17组,2022年72组),用于模型验证的共有43组(2019年8组,2022年35组)。Chl.a实测值与估算值对比见图 7(a)(b)。由图 7可见,设立阈值ρ(TSM)=20 mg/L,使用一元回归模型和多元回归模型在2019年和2022年反演Chl.a都可以得到较理想的结果,表明这2组Chl.a反演模型在不同年份间有一定适应性。ρ(TSM)<20 mg/L,部分多元回归与一元回归模型反演精度相当;ρ(TSM)≥20 mg/L,部分多元回归反演模型精度更高。总体来看,一元回归模型的RMSE=8.67 μg/L,MRE=20.4%;多元回归模型的RMSE=8.02 μg/L,MRE=18.4%。
选择单波段法、比值法与差值法,分平水期和丰水期建立波段组合与实测TSM的相关关系[25-26],结果见图 8。
由图 8可见,平水期时,波段比值RB4/RB3,RB4/RB2和RB3/RB8与实测TSM相关性较好;丰水期时,近红外单波段RB6—RB8A和波段比值RB6/RB5,RB3/RB8与实测TSM相关性更好。因此在平水期,选择相关性最高的RB4/RB2建立一元回归反演模型,同时利用偏最小二乘法综合RB4/RB3,RB4/RB2和RB3/RB8建立多元回归反演模型;在丰水期,选择相关性最高的单波段RB8建立一元回归反演模型,同时利用偏最小二乘法综合RB8,RB6/RB5和RB3/RB8建立多元回归反演模型。TSM反演模型及精度见表 2。
用于平水期TSM反演的实测数据共有81组,其中用于模型模拟的共有55组(2019年14组,2022年41组),用于模型验证的共有26组(2019年6组,2022年20组);用于丰水期TSM反演的实测数据共有79组,用于模型模拟的共有53组(2019年15组,2022年38组),用于模型验证的共有26组(2019年6组,2022年20组)。TSM实测值与估算值对比见图 9(a)(b)。由图 9可见,一元回归模型和多元回归模型在2019年和2022年反演TSM都可以得到较理想的精度,表明2组反演模型在不同年份间有一定适应性。在平水期和丰水期,使用多元回归方程反演的精度较高。总体来看,一元回归模型的RMSE=17.0 mg/L,MRE=25.7%;多元回归模型的RMSE=16.2 mg/L,MRE=23.3%。
首先,选择2.3.3节的多元回归模型,分平水期和丰水期分别模拟高邮湖的TSM,之后再以ρ(TSM)=20 mg/L为阈值,根据2.3.2节的多元回归模型,分别模拟高邮湖的Chl.a。与实测数据的采样时间进行比较,在2019年9月24日和2022年4月6日,共有9组准同步实测数据可以与遥感反演数据匹配。利用这9组数据对卫星反演结果进行模型精度验证,Chl.a的RMSE=14.5 μg/L,MRE=25.7%,TSM的RMSE=23.4 mg/L,MRE=31.9%,该反演精度可以接受。高邮湖Chl.a和TSM估算值与实测值对比见图 10(a)(b)。
根据3.1节的方法,利用哨兵-2 MSI反演高邮湖水质参数,得到高邮湖Chl.a和TSM的空间分布结果,见图 11(a)—(h)。由图 11可见,高邮湖Chl.a和TSM的空间分布与其围网养殖区域和地表入湖径流有一定关系。在2019年3月13日和2022年4月6日(平水期),高邮湖Chl.a和TSM总体值较低。在2022年4月6日,从淮河入江口到邵伯湖出水口的水流通道ρ(TSM)明显升高,表明有径流携带大量泥沙流入高邮湖,又从南部流出至邵伯湖,与2.1节中对实测值的分析一致。同时期ρ(Chl.a)没有明显变化,表明此时从淮河流入的水体ρ(Chl.a)不高。在2019年9月24日和2022年10月3日(丰水期),高邮湖ρ(Chl.a)和ρ(TSM)在有围网和靠近围网部分有所升高。在2019年9月24日,北部高邮湖ρ(Chl.a)和ρ(TSM)同时升高,疑似由同一股径流导致,可能是北方大汕子河长期关闭的闸门打开,导致一股长期封闭且富营养化的水体流入高邮湖。
(1) 在高邮湖用实测数据模拟哨兵-2 MSI数据,利用半分析法提取适宜波段组合,再与实测值拟合建立Chl.a和TSM反演模型。经验证,模型误差较小,且在不同年份有一定适应性,有实际应用价值。
(2) 应用准同步哨兵-2 MSI数据与反演模型对Chl.a和TSM进行反演,通过准同步实测数据验证,反演误差较小,表明本研究建立的模型可以在高邮湖进行Chl.a和TSM参数遥感反演。
(3) 采用本研究的Chl.a和TSM反演模型,通过对结果时空变化的初步分析,发现高邮湖水质参数的分布在一定程度上与围网养殖区域和入湖径流有关,后续须结合详细水文气象、渔业管理等数据进行对比分析,以实现对高邮湖水质变化情况的监测预警,为更科学地管理和保护高邮湖生态环境提供依据。
[1] |
吴传庆, 王桥, 杨志峰, 等. 太湖水华的遥感分析[J]. 中国环境监测, 2007(3): 52-56. |
[2] |
马荣华, 段洪涛, 唐军武, 等. 湖泊水环境遥感[M]. 北京: 科学出版社, 2010.
|
[3] |
陈宇洁, 陈志芳, 马德高, 等. 基于哨兵2号MSI数据的高邮湖水质参数遥感估算研究[J]. 环境科学与管理, 2022, 47(7): 77-81. |
[4] |
DEKKER A G. Detection of optical water quality parameters for eutrophic waters by high resolution remote sensing[D]. Doctorate Thesis Earth and Life Sciences, Amsterdam: Vrijie University, 1993.
|
[5] |
GITELSON A A, DALL'OLMO G, MOSES W, et al. A simple semi-analytical model for remote estimation of chlorophyll-a in turbid waters: Validation[J]. Remote Sensing of Environment, 2008, 112(9): 3582-3593. DOI:10.1016/j.rse.2008.04.015 |
[6] |
MISHRA S, MISHRA D R. Normalized difference chlorophyll index: A novel model for remote estimation of chlorophyll-a concentration in turbid productive waters[J]. Remote Sensing of Environment, 2012, 117: 394-406. DOI:10.1016/j.rse.2011.10.016 |
[7] |
孙乐成, 王娟, 王林. 基于实测数据的北海区海水表层叶绿素a浓度遥感定量反演[J]. 海洋开发与管理, 2019, 36(4): 34-38. |
[8] |
刘忠华, 李云梅, 檀静, 等. 太湖、巢湖水体总悬浮物浓度半分析反演模型构建及其适用性评价[J]. 环境科学, 2012, 33(9): 3000-3008. |
[9] |
姜广甲, 刘殿伟, 宋开山, 等. 基于半分析模型的石头口门水库总悬浮物浓度反演研究[J]. 遥感技术与应用, 2010, 25(1): 107-111. |
[10] |
高晨, 徐健, 高丹, 等. 基于GF-1与实测光谱数据鄱阳湖丰水期总悬浮物浓度反演[J]. 国土资源遥感, 2019, 31(1): 101-109. |
[11] |
陈玥, 管仪庆, 苗建中, 等. 基于长期水文变化的苏北高邮湖生态水位及保障程度[J]. 湖泊科学, 2017, 29(2): 398-408. |
[12] |
环境保护部. 水质叶绿素a的测定分光光度法: HJ 897—2017[S]. 北京: 中国环境科学出版社, 2017.
|
[13] |
国家环境保护局. 水质悬浮物的测定重量法: GB/T 11901—1989)[S]. 北京: 中国标准出版社, 1989.
|
[14] |
唐军武, 田国良, 汪小勇, 等. 水体光谱测量与分析Ⅰ: 水面以上测量法[J]. 遥感学报, 2004(1): 37-44. |
[15] |
MOREL A, GENTILI B. Diffuse reflectance of oceanic waters. Ⅲ. Implication of bidirectionality for the remote-sensing problem[J]. Applied Optics, 1996, 35(24): 4850-4862. DOI:10.1364/AO.35.004850 |
[16] |
MCFEETERS S K. The use of the Normalized Difference Water Index(NDWI) in the delineation of open water features[J]. International Journal of Remote Sensing, 1996, 17(7): 1425-1432. DOI:10.1080/01431169608948714 |
[17] |
扬州市水利局. 扬州市水资源公报(2019年)[R]. 2019.
|
[18] |
扬州市水利局. 扬州市水资源公报(2022年)[R]. 2022.
|
[19] |
杨敏, 商少凌, 汪文琦, 等. 水库水体近红外反射峰与叶绿素含量之间的关系初探[J]. 湖泊科学, 2009, 21(2): 228-233. |
[20] |
张雨萌, 王泉, 段春钰, 等. 基于Sentinel-2卫星的星云湖叶绿素a遥感估算研究[J]. 环境保护前沿, 2020, 10(1): 20-31. |
[21] |
周琳, 马荣华, 段洪涛, 等. 浑浊Ⅱ类水体叶绿素a浓度遥感反演(Ⅰ): 模型的选择[J]. 红外与毫米波学报, 2011, 30(6): 531-536. |
[22] |
徐祎凡, 施勇, 李云梅. 基于环境一号卫星高光谱数据的太湖富营养化遥感评价模型[J]. 长江流域资源与环境, 2014, 23(8): 1111-1118. |
[23] |
安如, 刘影影, 曲春梅, 等. NDCI法Ⅱ类水体叶绿素a浓度高光谱遥感数据估算[J]. 湖泊科学, 2013, 25(3): 437-444. |
[24] |
陈瑜丽, 沈芳. 长江口及邻近海域悬浮颗粒物对叶绿素a遥感反演算法的影响分析[J]. 遥感技术与应用, 2016, 31(1): 126-133. |
[25] |
高晨, 徐健, 高丹, 等. 基于GF-1与实测光谱数据鄱阳湖丰水期总悬浮物浓度反演[J]. 国土资源遥感, 2019, 31(1): 101-109. |
[26] |
侯琳琳, 马安青, 胡娟, 等. 胶州湾水体悬浮物浓度遥感反演模式优化研究[J]. 中国海洋大学学报(自然科学版), 2018, 48(10): 98-108. |