珠江口河网水体非光学活性水质参数高光谱反演

蔡建楠 , 刘海龙 , 姜波 , 陈吟晖 , 李杰鸿 , 吴思晓 , 梁建霞 , 黄华, 邢前国

1. 中山市环境监测站, 广东 中山 528403;

2. 中山市生态环境局, 广东 中山 528403;

3. 中国科学院烟台海岸带研究所, 中国科学院海岸带环境过程与生态修复重点实验室, 山东 烟台 264003;

4. 广东环境保护工程职业学院, 广东省无人机环保应用工程技术研究中心, 广东 佛山 528216

珠江口地区河网密布, 区域内河流众多, 随着城市化过程的不断推进, 水质问题越来越突出, 直接影响到居民生活质量, 与粤港澳大湾区生态环境建设的目标不符。为此, 珠江口各地投入了巨大的社会资源, 广泛开展了河流整治工作, 其中长期、密集的河流水质监测工作对于城市未达标水体的精准识别和对河流整治效果的综合评估具有重要意义(王晓 等, 2017; 张雪 等, 2019)。但传统的水质监测分析方法存在成本高昂, 采样及分析过程耗时长, 且仅代表特定区域和时段等缺点, 监测结果往往难以准确反映监测对象整体的水质状况, 难以满足当前形势下水环境质量评价的需要(Wang et al, 2004)。

高光谱水质遥感监测技术主要基于水体实测光谱, 通过建立光谱特征和水质指标的关系模型, 实现对水质指标的反演(王桂芬 等, 2010; 彭建 等, 2018)。高光谱遥感数据的光谱谱段丰富, 在水质参数定量化监测中具有良好的应用前景, 可作为传统水质监测方法的有效补充, 近年来得到迅速的发展(林剑远 等, 2019), 其中浊度、悬浮物、叶绿素a 等光学活性水质参数由于具有机理明确的光学响应特征而得到广泛的研究(Wang et al, 2019)。但水环境管理与整治工作中主要监控的指标为碳组分、营养盐等非光学活性水质参数, 该类参数具有光学响应波段不显著和复杂多样的特点, 因而开展该类参数的光谱特征研究及水质反演具有较高的科学意义和现实意义(邢前国, 2007)。近年来, 国内外学者采取了间接反演、半经验模型等方法开展了水体中总磷(TP)、总氮(TN)、溶解性有机碳(DOC)、氨氮(NH3-N)等非光学活性水质的高光谱反演研究(Song et al, 2012; Sudduth et al, 2015; Fichot et al, 2016; 张海威 等, 2017), 然而相关研究主要集中于湖泊水库, 应用于河流特别是水质复杂河网地区的研究报道较少。本文通过珠江口中山市各类型河流的高光谱监测和同步水质分析, 研究受测水体的高光谱特征及其与河流化学需氧量(CODCr)和总磷(TP)这两项水质指标浓度的关系, 并尝试建立高光谱数据与水质指标浓度的反演模型, 进而为河流非光学活性水质参数高光谱监测的实施及相关高光谱传感器的设计提供参考。

1 研究区概况

中山市位于珠江口西岸, 地处粤港澳大湾区的几何中心, 北连广州, 毗邻港澳, 陆域总面积约1783km2, 多年平均降雨量为1748.7mm, 多年平均水资源总量17.38 亿m3。中山市辖区内河流众多, 除磨刀门水道、鸡鸦水道、小榄水道、横门水道和石岐河五大主要河流外, 还有各种规模的河流超过1000 条(图1)。中山市经济发达, 产业集群密布, 近年来随着城市的快速发展和人口的不断增加, 部分河流水质受到了一定的污染。

图1 研究区河流水系示意图 Fig.1 Schematic diagram of river system in the study area

2 数据与方法

2.1 光谱测量和水样分析

2019 年7 月至9 月在中山市实施了多次河流现场光谱测量。测定点位的选取以图1 所示的中山市主要河流及其支流为主要对象, 兼顾考虑入河污染特征、流经区域类型、周边环境特点的多样性和代表性, 以尽可能覆盖中山市各种水质类型的河流, 最终获得151 个现场实测光谱数据和水质样品(代表性河流照片见图2)。使用的测量设备为美国Ocean Optics 公司的USB4000 光谱仪, 测量光谱范围为360～1000nm, 光谱分辨率<0.2nm。测量时选择晴天无云或少云天气状况下距正午3h 的时段内进行。在每个点位离河流水面高度0.5m 处, 测量河流水体的上行辐射亮度(Lu)、天空辐射亮度(Ls)和太阳下行辐照度(Ed)。每个点位测量前对仪器进行暗电流校正, 并使用标准灰板进行调试, 使其信号基本保持在80%量程范围内。进行3 次测量后取平均值作为光 谱测定结果, 利用所测数据代入公式(1), 即可计算测点的光谱反射率R(无量纲)。

式中: ρ为天空光在水体表面的菲涅尔反射率, 通常取经验值0.028(Mobley, 1999)。

监测点的水样采集与光谱测量同步进行, 水样分析项目为CODCr、TP。水样取河流表层0～20cm处的混合样, 按照相应的国家技术规范要求保存并送至实验室进行分析。实验室分析过程采用平行样、标准样品、加标回收等质量控制方法, 以确保分析的准确度和精密度。全部样品测得的水质指标特征描述见表1。

图2 代表性监测河流照片 黄色箭头标识处为监测点位置 Fig.2 Photos of typical rivers (The position of the arrow indicates the monitoring points)

表1 样品水质特征描述 Tab. 1 Description of characteristics of water quality

2.2 高光谱数据预处理

由于本研究获得的原始高光谱数据的波段数较多, 为提高计算效率和反演精度, 且有利于与其他研究成果相比较, 数据处理过程中仅保留400～900nm 波段范围内的原始数据, 并将其光谱分辨率重采样为1nm, 以降低信息冗余, 重采样后的数据将作为本次研究的数据。光谱分辨率重采样采用分段三次 Hermite 插值多项式(Piecewise Cubic Hermite Interpolating Polynomial, PCHIP)法, 该方法可保持原始数据的形状和相应的单调性(Iqbal et al, 2013), 符合本研究后续数据分析需要。

2.3 相关分析与模型评价

2.3.1 水质指标与光谱反射率相关系数的计算

各水质指标浓度值与水体不同波段高光谱反射率的相关性分析有助于识别对水质指标变化敏感的特征波段, 从而进一步筛选可用于水质反演的最优特征波段及其组合。本研究中各水质指标浓度值与水体高光谱原始反射率相关系数的计算方法为:

2.3.2 反演模型的建立与评价

以相关分析所确定的特征波段反射率及其组合作为自变量, 以水质指标实测浓度值作为因变量, 建立相关水质反演模型。本研究从151 个光谱及水质样品中随机选择130 个样本作为建模数据集, 并计算回归模型的决定系数 R2和均方根误差 RMSE(Root Mean Squared Error), 以评价拟合精度并确定最优模型。在确定最优模型后, 利用未参与建模的21 个样本作为验证数据集, 进行估算结果与实测值的对比分析, 以绝对误差 AE(Absolute Error)和平均相对误差MRE(Mean Relative Error)评价模型的泛化能力。R2、RMSE、AE 和MRE 的计算公式分别为:

3 结果与讨论

3.1 河流水体高光谱反射率特征

本研究中的实测城市河流水体在400～900nm 波段范围内各波段的光谱反射率均值及第10、第90百分位数如图3 所示。所测河流水体高光谱总体特征表现为: 400～500nm 范围内的水体光谱反射率总体较低; 随后水体光谱反射率随波长增加而上升, 并呈现“三峰二谷”的光谱形态特征, 即在585nm、710nm、810nm 附近有明显的反射峰, 在670nm、760nm 附近有明显的反射低谷。根据相关研究, 地表水在585nm、710nm 附近的反射峰主要由后向散射和在该波段处较弱的吸收共同作用所致; 810nm附近的反射峰则是由于该处水分子吸收系数较小, 水体悬浮颗粒物散射程度加强而进一步提高了水体的反射能力(陈雯扬 等, 2010); 670nm、760nm 附近的反射低谷分别由藻类和水分子吸收所致(陈楚群 等, 2001; 刘志国 等, 2007; 杨锦坤 等, 2007)。上述分析表明各河流水体光谱反射率主要受悬浮颗粒物和叶绿素a 影响(刘志国 等, 2007)。

图3 实测水体高光谱反射率 Fig.3 Hyperspectral reflectance of sampled water

3.2 水质指标与水体光谱相关性分析

3.2.1 水质指标与高光谱原始反射率的相关性

根据图4 显示, 各水体样品中CODCr、TP 浓度的自然对数值与水体光谱反射率的相关性特征基本趋于一致, 即在400～900nm 波段上总体呈负相关, 该负相关性先增大后减少; 在550～680nm 波段范围内具有较高的相关性, CODCr、TP 浓度对数值与光谱反射率相关性最高, 分别可达–0.55 和–0.46; 相关程度在680nm 之后的波段中迅速下降, 在720nm以上波长范围内相关系数均不到–0.15。

图4 不同波段光谱反射率与CODCr 和TP 浓度对数值的相关性 Fig.4 Relationships between the logarithm of CODCr, TP concentration and hyperspectral reflectance

3.2.2 水质指标与光谱反射率比值的相关性

对任意两个波段的原始光谱反射率数据进行比值计算, 形成新的光谱变量, 并与水质指标进行回归分析, 进而突出光谱的吸收与反射特征, 这一方法已被广泛应用于水质高光谱反演经验模型的构建中(Song et al, 2012)。本研究将每个样本光谱数据中501 个波段的反射率进行两两比值计算, 形成151个501×501 的比值数据矩阵后, 再按照公式(2)分别计算矩阵中对应光谱反射率比值与CODCr和TP 浓度对数值的相关系数, 以识别出与水质指标相关性最密切的波段组合, 相关性矩阵计算结果见图5。

CODCr、TP 浓度对数值与反射率比值的相关系数最高分别为0.795 和0.731, 对应的波段组合分别为R675/R794、R690/R815。结合图3 的光谱形态分析, 该最优波段组合基本上与光谱反射率曲线吸收峰和近红外反射峰分别对应。这表明与直接利用单波段原始数据相比, 反射率比值处理更加强化了光谱的差异特征, 有助于进一步筛选与CODCr、TP 浓度值相关性高的波段组合。

根据图5 所示, 与CODCr、TP 浓度对数值相关性较高的主要是红波段与近红外反射率间的比值, 如R675/R794、R690/R815, 其波长与叶绿素的特征吸收、泥沙的反射峰位置一致, 说明非光学活性参数CODCr、TP 与光学敏感性参数中的叶绿素及悬浮物有一定的相关性: 由叶绿素表征的浮游植物本身就是CODCr的重要来源, TP 则可能主要以悬浮颗粒形式出现, 且浮游植物本身也对悬浮物浓度有贡献。

图5 CODCr(a)和TP(b)浓度对数值与光谱反射率比值的相关性矩阵 Fig.5 Correlation matrix of relationships between the logarithm of CODCr (a) and TP (b) concentrations and spectral reflection ratio indexes

3.3 基于反射率比值的水质反演模型

3.3.1 模型的建立

基于上述反射率比值的筛选结果, 选择R675/R794和R690/R815作为CODCr、TP 反演模型的自变量, 分别采用线性、指数、幂函数三种拟合方式建立回归模型。建立的水质指标光谱反演回归模型及结果评价见表2、图6。各拟合方程参数均通过P<0.05 的显著性检验。其中, CODCr、TP 最优拟合模型为幂函数模型, 模型决定系数 R2分别达到0.661 和0.569; 总体上CODCr拟合效果略优于TP。

表2 CODCr 和TP 浓度反演模型拟合结果比较 Tab. 2 Fitting results of CODCr and TP concentration retrieval models

图6 CODCr(a)、TP(b)浓度最优反演模型拟合结果 Fig.6 Fitting results of the best CODCr (a) and TP (b) concentration retrieval models

3.3.2 模型的精度验证及评价

基于上述CODCr、TP 浓度高光谱反演最优拟合模型, 利用验证数据集进行模型估算结果与实测值的对比分析。结果(表 3、图 7)显示, CODCr、TP 浓度反演模型的估算精度良好, 两者的MRE 分别为27.2%和32.1%, 表明采用最优反射率比值作为自变量构建的幂函数模型具有良好的适用性, 可较好地反演水体中的CODCr、 TP 浓度值。

表3 CODCr 和TP 浓度估算结果评价 Tab. 3 Evaluation of estimation results for CODCr and TP concentrations

图7 CODCr(a)、TP(b)浓度实测值与模型估算值比较 Fig.7 Comparison of measured and estimated values of CODCr (a) and TP (b) concentrations

4 结论

本文通过对广东省中山市151 个河网水体监测点的高光谱监测和水质分析, 获得了水体的高光谱反射率特征, 并构建了CODCr和TP 浓度反演模型, 得到以下结论:

1) 各河流水体光谱反射率主要受悬浮颗粒物和叶绿素a 影响。在500～680nm 波段范围内, CODCr和TP 浓度对数值与单波段光谱反射率存在较强的 负相关关系。与单波段相比, 特定波段的反射率比值与CODCr、TP 浓度值的相关性较高, 与CODCr和 TP 浓度值相关性最大的波段组合分别为R675/R794、R690/R815。

2) 基于上述反射率波段比值组合建立的CODCr和TP 浓度幂函数反演模型的决定系数R2分别为 0.661 和 0.569, 估算平均相对误差分别为27.2%和32.1%, 表明高光谱技术在珠江口河网水体CODCr、TP 浓度的水质反演中具有一定的应用潜景。