王 林,王 祥,周 超,王新新,孟庆辉,陈艳拢

国家海洋环境监测中心,辽宁 大连 116023

引 言

我国近岸和近海海域的主要污染物80%以上来自陆源排污,每年上百亿吨的工业和生活废水将大量的污染物质携带排放入海,造成近岸海域水质普遍较差;而入海河流在陆源污染物向海洋输送的过程中发挥着重要作用。近几年,随着习近平生态文明思想的深入贯彻,“河长制”、“湾长制”等管理制度的广泛推行,以及污染防治攻坚战的全面打响,入海河流与近岸海域的水质状况得到稳步改善,但尚存在波动情况,污染防治形势依然严峻。连云港地处江苏东北部,海州湾西岸,东濒黄海,入海河流主要包括灌河、古泊善后河及临洪河等,2020年以来水质明显提升,15条入海河流平均水质全部消除劣Ⅴ类,但氮磷等营养盐仍为主要污染物质,部分入海河流叶绿素a(chlorophyll a,Chl a)浓度非常高,呈现明显的富营养化状态。

随着各类卫星影像的不断涌现,遥感技术在海洋、湖泊、水库及河流等水生态环境监测与评价中得到广泛应用,尤其针对海洋和大型湖泊生态环境要素的定量遥感研究已持续多年[1-6],技术相对较为成熟。近几年,高分卫星与无人机遥感技术得到快速发展,针对狭长河流水体生态环境要素的定量遥感研究逐渐出现[7],反演参数包括浊度(turbidity,Tur)[8]、悬浮颗粒物(suspended particulate matter,SPM)[9]、总氮(total nitrogen,TN)[10-11]、总磷(total phosphorus,TP)[10-11]、叶绿素a[12-15]及水体营养状态[16]等。富营养化是入海河流和近岸海域典型的水质污染问题,是水体中营养物质富集导致浮游藻类过度生长的生态异常现象,对区域生态环境及人类生活安全构成严重威胁[17]。因此,对入海河流中叶绿素a和综合营养状态指数(trophic level index,TLI(Σ))的定量遥感研究具有重要的实际应用价值。Shi等[14]采用Gaofen-6和Sentinel-2卫星数据,通过机器学习模型反演了北京温榆河、北运河等水体中叶绿素a浓度,较好反映了内陆河流叶绿素a浓度的空间分布特征。马驰[15]以东北松嫩平原河流、湖泊水体为研究对象,建立了基于Sentinel-2A窄近红外与深蓝波段的叶绿素a浓度一元二次模型,反演结果显示研究区内叶绿素a浓度较低,水质较好。张勇等[16]基于Landsat-8卫星影像,建立了合肥市环城河总氮、总磷及氨氮浓度遥感反演模型,并运用评分法进行了富营养化评价,最终确定了各河段水体的富营养化状态。入海河流下游受海洋潮汐作用影响,存在一定的咸淡水交界区域,与一般的内陆河流水体有一定差别,目前针对入海河流水体叶绿素a浓度及综合营养状态指数的定量遥感研究还相对较少。

研究以连云港市境内灌河、古泊善后河、蔷薇河及临洪河等入海河流为研究对象,采用Sentinel-2 MSI影像及叶绿素a、总氮、总磷等现场实测数据,建立了该区域主要入海河流叶绿素a及综合营养状态指数的遥感定量反演模型。

1 实验部分

1.1 研究区域与调查站点分布

研究区域选择连云港市四条主要的入海河流,包括灌河、古泊善后河、蔷薇河及临洪河。其中,灌河是苏北地区唯一在干流上没有建闸的黄金入海通道,受两岸工业排污和船舶航运等影响,水质较差且水体常年浑浊;古泊善后河作为苏北地区行洪排涝的骨干河道,具备防洪、排涝、供水、灌溉、航运等综合功能,受沿线污染汇入影响,河流水质超标事件时有发生[18];蔷薇河沿线入河支流排口多,主要受农村农业面源污染、酸洗石英砂点源污染等因素影响,河水水质很不稳定[19];临洪河是蔷薇河的下游河段,经临洪闸分界,于临洪河口入海,故将蔷薇河与临洪河作为一条入海河流进行后续研究。

现场调查于2022年6月24、25、27日进行,调查参数包括总氮、总磷及叶绿素a等水质要素。共计23个站点,其中蔷薇河/临洪河9个,古泊善后河7个,灌河7个,研究区域范围及具体站点布设情况如图1所示。

图1 研究区域外业调查站点Fig.1 Field survey stations in the study area

1.2 数据获取

1.2.1 水质参数

叶绿素a浓度测量采用荧光法;总氮浓度测量采用碱性过硫酸钾消解紫外分光光度法;总磷浓度测量采用钼酸铵分光光度法。

1.2.2 卫星影像数据

卫星影像数据从欧洲航天局哥白尼开放数据访问中心网站(https://scihub.copernicus.eu/)下载。采用2022年6月25日Sentinel-2A MSI L2A影像数据与现场实测水质数据进行匹配,提取与实测水质站点对应的光谱反射率数据;建模后应用于2022年6月25日Sentinel-2A MSI L2A影像,得到连云港主要入海河流叶绿素a和综合营养状态指数遥感定量反演结果,并进行空间分布特征分析。

1.3 综合营养状态指数分析方法

采用中国环境监测总站关于“湖泊(水库)富营养化评价方法及分级技术规定”方法进行综合营养状态指数计算,见式(1)和式(2)

(1)

(2)

式(1)和式(2)中:TLI(Σ)(无量纲)为综合营养状态指数;Wj(无量纲)为第j种评价指标营养状态指数的相关权重;TLI(j)(无量纲)为第j种评价指标的营养状态指数。

综合营养状态指数的分级为:TLI(Σ)<30时为贫营养;30≤TLI(Σ)≤50时为中营养;5070时为重度富营养。

河流水体富营养化受多种环境因子影响,其中最为重要的两个参数即是氮、磷,而叶绿素a则是富营养化最重要的表征指标[20],对水体综合营养状态的评估均具有决定作用。因此,将TN、TP及Chl a作为综合营养状态指数的评价指标,分别计算得到各评价指标的营养状态指数TLI(j),确定其相关权重指数Wj分别为0.282 76、0.296 72、0.420 52,进而计算得到综合营养状态指数TLI(Σ)。

1.4 模型评价

为了筛选叶绿素a和综合营养状态指数遥感反演的最佳波段,采用相关系数(r)进行评价,见式(3)

(3)

叶绿素a和综合营养状态指数反演模型的优劣,采用决定系数(R2)、平均绝对百分比误差(mean absolute percentage error,MAPE)和均方根误差(root mean square error,RMSE)进行评价,见式(4)、式(5)和式(6)

(4)

(5)

(6)

2 结果与讨论

2.1 入海河流营养状态背景

现场实测数据分析可知,蔷薇河/临洪河、古泊善后河、灌河水质均较差,总氮超标现象最为突出,大多站点未达到水质Ⅲ类考核要求,劣五类水质站点比例达65%,可能与现场调查时间为丰水期,受沿线污染大量汇入影响较大有关。进一步对比发现,Chl a、TP及TLI(Σ)均值以蔷薇河/临洪河最高,分别达42.83 μg·L-1、0.30 mg·L-1、67.75,古泊善后河次之,灌河最低;TN均值以灌河最高,达2.96 mg·L-1,蔷薇河/临洪河次之,古泊善后河最低。具体调查结果如表1所示。

表1 连云港主要入海河流的营养状态参数统计Table 1 Statistics on trophic status parameters for Lianyungang major seagoing rivers

2.2 遥感反演敏感波段筛选

提取与水质现场实测站点相匹配的水体光谱反射率R(λ),并分别与Chl a、TLI(Σ)进行线性拟合,各波段的相关系数如图2所示。整体看来,400～2 200 nm波段范围内R(λ)与Chl a、TLI(Σ)均存在负相关关系,即R越高水体越浑浊,Chl a、TLI(Σ)则越低,表明浑浊水体光线透射较差,不利于浮游藻类生存生长;且TLI(Σ)与R(λ)的相关性略高于Chl a,与TLI(Σ)计算时除Chl a之外还引入了TP、TN有关;400～750 nm可见光波段R(λ)与Chl a、TLI(Σ)的相关性明显高于其他波段,可见光波段的相关性先增加后减小,490、560、665 nm三个波段的相关性最佳,其中R(λ)与Chl a的相关系数分别为-0.697、-0.681、-0.693,R(λ)与TLI(Σ)的相关系数分别为-0.728、-0.744、-0.706。故将490、560、665 nm三个波段作为连云港主要入海河流叶绿素a和综合营养状态指数遥感定量反演的敏感波段。

图2 Sentinel-2 MSI主要波段反射率与叶绿素a浓度、综合营养状态指数的相关系数Fig.2 Correlation coefficients of Sentinel-2 MSI major band reflectance with Chl a concentration and TLI(Σ)

2.3 遥感反演模型建立

以490、560、665 nm为敏感波段,并选择单波段、双波段、三波段组合形式,在Chl a对数坐标下采用线性、乘幂及指数模型进行拟合分析,R2最高而MAPE、RMSE最低时,则模型最佳,结果如表2所示。可发现,单波段模型中,以R(665)为自变量,Chl a对数坐标下的乘幂函数表现最佳,拟合R2为0.67,MAPE为47.34%,RMSE为12.89 μg·L-1;双波段模型中,常用的波段比值模型表现较差,MAPE均超过60%,不适用于连云港主要入海河流水体,而二元一次函数模型略优于波段比值模型,但整体表现均较差;采用三元一次函数拟合的三波段模型表现一般,也不适用于该区域入海河流水体。因此,将lg(Chl a)=0.267×R(665)-0.982作为连云港主要入海河流叶绿素a浓度遥感定量反演的最优模型。

表2 连云港主要入海河流叶绿素a浓度遥感定量反演模型的精度评价Table 2 Accuracy evaluation of Chl a concentration in Lianyungang major seagoing rivers by using remote sensing inversion models

同样以490、560、665 nm为敏感波段,选择单波段、双波段、三波段组合形式,并采用线性、乘幂函数进行综合营养状态指数模型拟合分析,R2最高而MAPE、RMSE最低时,则模型最佳,结果如表3所示。可发现,单波段模型中,以R(560)为自变量,并采用乘幂函数进行拟合的模型表现最佳,拟合R2为0.61,MAPE为4.36%,RMSE为3.45;双波段模型中,波段比值模型表现较差,MAPE均超过6%,但二元一次函数模型明显优于波段比值模型;采用三元一次函数拟合的三波段模型表现较佳,但逊于单波段模型。因此,将TLI(Σ)=44.898×R(560)-0.241作为连云港主要入海河流综合营养状态指数遥感定量反演的最优模型。

表3 连云港主要入海河流综合营养状态指数遥感定量反演模型的精度评价Table 3 Accuracy evaluation of TLI(Σ) in Lianyungang major seagoing rivers by using remote sensing inversion models

2.4 空间分布特征

将上节中最优模型应用于2022年6月25日Sentinel-2A MSI L2A影像,反演得到连云港主要入海河流Chl a和TLI(Σ),如图3和图4所示。整体看来,蔷薇河/临洪河Chl a最高,古泊善后河次之,灌河最低,且各河流上游Chl a普遍高于下游,可能与涨潮时海水涌入河流下游有关;类似地,蔷薇河/临洪河TLI(Σ)最高,古泊善后河次之,灌河最低,且各河流上游富营养化程度普遍高于下游。进一步分析可知,蔷薇河/临洪河80%以上水域Chl a大于10 μg·L-1,且大于20 μg·L-1水域占比48.68%;而灌河90%以上水域Chl a则小于10 μg·L-1,小于5 μg·L-1水域占比53.89%;古泊善后河85%以上水域Chl a在5～20 μg·L-1范围,5～10 μg·L-1水域占比52.52%。蔷薇河/临洪河、古泊善后河及灌河基本均为富营养水体,其中蔷薇河/临洪河以中度富营养水体为主占比91.15%,重度富营养水体占比8.62%;古泊善后河也以中度富营养水体为主占比85.70%,而轻度富营养水体占比11.85%;灌河则以轻度富营养水体为主占比72.52%,中度富营养水体占比25.04%。具体统计结果如表4所示。

表4 不同区间叶绿素a和综合营养状态指数所占水域面积的比例统计Table 4 Interval statistics of watershed areas with different levels of Chl a and TLI(Σ)

图3 连云港主要入海河流叶绿素a浓度(2022.06.25)Fig.3 Chl a concentration in major seagoing rivers in Lianyungang (2022.06.25)

图4 连云港主要入海河流综合营养状态指数(2022.06.25)Fig.4 TLI(Σ) in major seagoing rivers in Lianyungang (2022.06.25)

以上模型研究发现,分别以665和560 nm为敏感波段的单波段模型更适用于连云港主要入海河流叶绿素a浓度和综合营养状态指数的遥感定量反演,这与湖泊、海洋水体存在一定差别。在湖泊、海洋相关研究中,560和665 nm波段多用于悬浮物浓度的遥感反演,且560 nm波段还作为参比波段用于叶绿素a浓度的遥感反演。分析原因在于,与湖泊或海洋相比,河流水体的浑浊度更高,水体光学特性更为复杂,不同河流或河段光线透射进入水体的深度存在一定差异,导致浮游藻类的生存生长状况呈现较大的不同,进而使得叶绿素a浓度、综合营养状态指数与水体浑浊度具有负相关关系,且常用于浑浊度遥感反演的特征波段成为叶绿素a浓度、综合营养状态指数反演的最优波段。

将反演模型应用于卫星影像,得到灌河、古泊善后河、蔷薇河/临洪河Chl a和TLI(Σ)空间分布结果,发现灌河Chl a、TLI(Σ)均最低,古泊善后河次之,蔷薇河/临洪河最高,这与河流水体的浑浊程度有直接关联。进一步对比分析可知,灌河各调查站点SPM均值为443.81 mg·L-1,古泊善后河为65.41 mg·L-1,蔷薇河/临洪河为33.34 mg·L-1,Chl a、TLI(Σ)与SPM呈负相关关系,表明SPM越高,水体越浑浊,光线透射越差,不利于浮游藻类生长,Chl a越低;TLI(Σ)由Chl a、TN、TP计算而得,其中Chl a是TLI(Σ)最重要的表征参数,权重指数最高,对TLI(Σ)影响最大,就灌河而言,即使TN、TP相对较高,但其Chl a尤其低,综合作用下TLI(Σ)仍最低。

3 结 论

(1)受河流沿线污染汇入影响,蔷薇河/临洪河、古泊善后河、灌河水质均较差。其中,蔷薇河/临洪河Chl a、TP及TLI(Σ)均值最高,古泊善后河次之,灌河最低;灌河TN均值最高,蔷薇河/临洪河次之,古泊善后河最低。

(2)叶绿素a浓度、综合营养状态指数与可见光波段反射率的相关性明显高于其他波段,最佳波段为490、560和665 nm。其中,R(λ)与Chl a的相关系数分别为-0.697、-0.681、-0.693,R(λ)与TLI(Σ)的相关系数分别为-0.728、-0.744、-0.706。

(3)以R(665)为自变量,在叶绿素a浓度对数坐标下的乘幂模型为其遥感定量反演的最优模型:lg(Chl a)=0.267×R(665)-0.982(R2=0.67,MAPE=47.34%,RMSE=12.89 μg·L-1);而以R(560)为自变量的乘幂模型是综合营养状态指数遥感定量反演的最优模型:TLI(Σ)=44.898×R(560)-0.241(R2=0.61,MAPE=4.36%,RMSE=3.45)。

(4)将最优模型应用于2022年6月25日Sentinel-2A MSI L2A影像,得到连云港主要入海河流叶绿素a浓度和综合营养状态指数的遥感反演结果,发现蔷薇河/临洪河中度富营养水体的占比为91.15%,重度富营养水体占比为8.62%;古泊善后河中度富营养水体占比为85.70%,轻度富营养水体占比11.85%;而灌河则以轻度富营养水体为主占比72.52%,中度富营养水体占比25.04%。