杨传训,李 勇,杨 骥,舒思京
(广东省科学院广州地理研究所,广东 广州 510070)
水体悬浮泥沙浓度大小决定水体透度、浑浊度和水色等光学性质,对河口自然环境、生态环境及港航通道等具有重要影响[1-4],河口海岸作为连接淡水水生系统和海洋水生生态系统,是最重要的地球生态系统之一,对全球生物元素循环和气候具有公认的重要性[5-6]。
当前监测悬浮泥沙浓度监测主要以现场实测和遥感反演两种方式为主[7-8]。站点式监测能够在时间尺度上连续、准确获取悬浮泥沙含量,但站点建设和维护成本高,站点数量少、分布稀疏,无法进行大面积的空间连续监测[8-9];遥感反演手段能够实现大范围的空间连续监测,能够较好地体现悬浮泥沙空间连续变化[10-11]。在珠江口悬浮泥沙遥感反演模式研究上,文献[12]在总结和对比前人线性关系、对数关系、Gordon、负指数关系基础上,提出并推导了悬浮泥沙遥感定量的统一模式,并利用 TM 影像数据反演珠江口悬浮泥沙,证明了统一模式与其他模式相比,精度较高;文献[13]建立了TM 辐射亮度值与悬浮泥沙浓度的多波段对数模式;文献[14]推导出悬浮泥沙定量遥感综合模式,该模式在条件简化下包含前人研究的半经验模型,并在珠江口取得成功应用。
当前,机器学习方法在水质遥感监测中已被广泛使用[15-17],本文以珠江口悬浮泥沙水质采样数据及高光谱数据为基础,构建珠江口人工神经网络、支持向量回归机、随机森林悬浮泥沙反演模型,实现对2006年12月19日TM影像、2013年9月10日Landsat OLI影像及2020年4月29日Sentinel-2影像数的反演与分析。
珠江入海口为研究区域,该区域多年平均含沙量约为0.126~0.334 kg/m3,断面最大含沙量为1.23~4.08 kg/m3。流域内年均降水量约2000 mm,具有典型的汛期和枯水期,汛期(4—9月)降水量占全年的67.6%~83.5%。珠江口水域属于典型的不规则半日混合潮型,大约12 h为一个平均潮汐周期,为弱潮河口,高潮位与低潮位之间的潮位差较小,属于典型的低潮差入海口。
水质及光谱数据:为研究珠江口潮汐效应改正的珠江口悬浮泥沙遥感反演与预测分析,共在珠江口的多处地方采集了229个点的水质参数及地面同步高光谱数据。
影像数据:采用的遥感影像数据为2006年12月19日TM影像数据、2013年9月10日Landsat OLI影像数据及2020年4月29日Sentinel-2影像数据。
从229个采集点位中选择出18条典型点位采集的光谱曲线(如图1所示),悬浮泥沙浓度值在1.6~456.2 mg/L中逐渐增大,结果表明在400~900 nm波段范围,水体的反射率值随着悬浮泥沙浓度增加而增加。从实际采样的珠江口高光谱实测数据发现珠江口悬浮泥沙存在两个反射峰值,第一反射峰位置在绿光、红光波段(570~690 nm),第二反射峰位置在近红外波段(780~830 nm),在含沙量较低时,第一反射峰值集中在绿光波段(560 nm)且显著高于第二反射峰值。随含沙量浓度的增加,第一反射峰值波段集中绿光、红光波段(550~670 nm),同时第二反射峰的反射率逐渐升高,第一反射峰值与第二反射峰值差逐步缩小。
图1 ASD地物光谱仪实测的水体光谱数据(N=18)
技术路线如图2所示。
图2 技术路线
2.1.1 人工神经网络
使用多层人工神经网络,该模型包括1个输入层,1~4个隐藏层和1个输出层。输入层是诊断光谱变量,系水质光学特征参数(OAWQ)与原始光谱波段和波段比值之间的高相关系数确定。第一隐藏层神经元的值等于输入参数乘以其连接权重加上偏差。第二层除了使用第一隐藏层神经元的值作为输入参数外,计算过程类似于第一层。输出层为OAWQ。假设输入与第一隐藏层的转换关系为Ai,两个隐藏层的关系为AAi。Ai和AAi可表示为
(1)
式中,B1和B2为激活函数;b1和b2为拟合过程中两层之间的偏差(Biasi);n为节点数;j为前一层的第j个节点;i为当前层中的第i个节点;k为两层之间的第k次迭代计算。在第k次迭代中,是连接输入层第j个节点和第一隐藏层第i个节点之间的权重;是两个隐藏层之间的连接权重。在OAWQi估算中,第二隐藏层和输出层之间的转换关系可表示为
(2)
2.1.2 支持向量回归机
支持向量机是一种相对简单的监督机器学习算法,用于解决分类或回归问题。支持向量回归机在长江口等其他流域区域在悬浮泥沙研究中也都表现出不错的结果[18]。常用的核函数有高斯径向基核、多项式核、sigmoid核、RBF核函数。本文采用RBF核函数,公式为
(3)
式中,γ为核参数。
2.1.3 随机森林模型
随机森林模型是一种基于分类树的机器学习方法,具有可输入大量的变数,分类规则清晰,结果容易理解,准确度高、计算量较小,实现速度快等优点,能够用于回归、分类、预测等分析。
随机森林模型原理是采用bootsrap重抽样方法从原始样本中抽取N个样本,对N个样本分别建立N个决策树模型,通过组成多棵决策树最终通过投票得出模拟结果[19]。
人工神经网络(artificial neural network, ANN)回归算法在悬浮泥沙浓度(total suspended solids, TSS)遥感反演估算结果如图3(a)所示。结果显示训练集精度和测试集精度R2均大于0.85(训练集R2=0.86,测试集R2=0.88),训练集RMSE=7.6 mg/L,MAE=4.76 mg/L,测试集RMSE=8.1 mg/L,MAE=5.54 mg/L,PRD=2.31。结合散点图发现,散点能较均匀分布1∶1线两侧,拥有接近于1的Slope(0.88 vs 0.76),但是训练、测试采样点存在部分饱和现象。
图3 基于ANN、SVR、RF的TSS遥感反演估算
支持向量回归机(support vector regression,SVR)算法在TSS遥感反演估算结果如图3(b)所示。结果显示测试集R2略优于训练集(训练集R2=0.84,测试集R2=0.86)。训练集RMSE=5.96 mg/L,MAE=5.01 mg/L,测试集RMSE=6.53 mg/L,MAE=5.22 mg/L,PRD=2.3。
随机森林(random forest, RF)算法在TSS遥感反演估算如图3(c)所示。结果显示测试集R2均大于0.90(训练集R2=0.91,测试集R2=0.95)。训练集RMSE=4.68 mg/L,MAE=3.89 mg/L,测试集RMSE=5.41 mg/L,MAE=4.18 mg/L,PRD=2.12。结合散点图发现散点能较均匀分布1∶1线两侧,拥有接近于1的Slope(0.90 vs 0.80)。
相对人工神经网络和支持向量回归机方法,随机森林模型预测方法预测精度高。3种机器学习模型精度对比见表1,与ANN相比,R2由0.88提升至0.95,RMSE降低2.69 mg/L,MAE降低1.36 mg/L。与SVR相比,R2由0.86提升到0.95,RMSE降低1.12 mg/L,MAE降低1.04 mg/L,主要原因为ANN、SVR黑箱模型内部的运行机理复杂,对影响因子权重难以确定,加上这两种模型对原始数据量要求高、计算量大等因素。
表1 各模型精度对比
随机森林模型分类规则清晰,结果容易理解,准确度高、计算量较小,实现速度快,且相关系数和预测精度明显优于ANN、SVR等其他模型。与ANN和SVR模型相比,RF模型在TSS水质参数中的模拟精度均有提升。因此,随机森林模型更适用于水体悬浮物监测遥感估算。
结合地面采样的229个光谱数据和水质数据,找到红光、近红外为悬浮泥沙敏感波段,而常用的传感器TM5、ETM+、OLI、Sentinel-2都存在绿光与近红外通道,可用于悬浮泥沙浓度的反演。但是多光谱传感器波段较宽(几十至几百nm),直接将构建的高光谱波段(1~20 nm)模型应用于多光谱传感器可能存在较大误差,因此需要分析宽波段下悬浮泥沙反演模型的适用性及构建。TM5、ETM+、OLI、Sentinel-2光谱响应函数如图4所示。绿光波段整体范围为0.50~0.65 μm,中心波长为0.55 μm。近红波段整体范围为0.71~0.93 μm。TM与ETM+光谱响应函数差异较小,Sentinel-2与OLI差异明显。
图4 绿光与近红外多传感器光谱响应函数(TM5、ETM+、OLI、Sentinel-2)
由于ASD光谱仪采集的数据光谱分辨率较高、通道窄,而不同遥感影像每个通道的光谱分辨率不一致,为了研究多传感器模型效果,在模型校正与验证前,将野外实测的229条反射率采集到TM5、ETM+、OLI、Sentinel-2波段传感器。具体计算公式为
(4)
式中,λ为波长;λmin为通道的起始波长;λmax为通道的终止波长;r(λ)为对应λ波长的地表反射率;f(λ)为光谱响应函数。
随机筛选152条光谱重采样数据,基于经验方法分析不同传感器通道响应对模型的适用性,经验模型校正结果如图5所示。上述模型拟合效果均较好(R2>0.8),说明将高光谱采集到传感器之后近红与绿光波段的比值能够用于悬浮泥沙的反演,模型都满足二次函数关系,模型形式一致,发现近红与绿波段的比值能削弱外界影响,构建模型精度较高,适用多源传感器较稳定。但是4种传感器模型系数差异显著,因此在应用模型时,不同传感器须获取对应的光谱响应函数,将实测高光谱数据重采样到对应传感器通道,对模型进行重新校正。
图5 多传感器最优悬浮泥沙反演模型
将重采样剩余的77个样本数据应用于校正的经验模型的精度评价,评价结果如图6所示,实测悬浮泥沙与模型计算悬浮泥沙R2均大于0.85,RMSE低于27.46,模型都具有较高的精度,其中Sentinel-2与ETM+传感器优于TM5与OLI传感器。验证散点图斜率为0.80~0.85之间,斜率小于1。由于不同传感的响应函数存在差异,将高光谱重采样到对应传感器波段并对模型校正到相应的传感器通道后4种传感器均具有较高的精度,模型计算结果略低于实测结果。
图6 多光谱传感器的经验模型精度评价
在机器学习中选择模型精度较好的随机森林模型中随机抽选40个样本数据进行验证,评价结果如图7所示,随机森林模型的精度与经验半经验模型相比有一定提高,在TM5、TM7、OLI、Sentinel-2 不同传感器上模型精度R2均大于0.9。TM5传感器RMSE=17.86 mg/L,TM7传感器RMSE=7.00 mg/L,传感器RMSE=21.22 mg/L,Sentinel-2传感器RMSE=11.38 mg/L,整体误差较小,可用于悬浮泥沙浓度反演监测。
图7 随机森林模型精度评价
应用模型评价精度最高的随机森林学习算法,基于2006年12月19日TM影像数据、2013年9月10日Landsat OLI影像数据及2020年4月29日Sentinel-2影像数据实现对珠江口悬浮泥沙反演、精度评价及空间制图,反演结果如图8所示。
图8 珠江口悬浮物浓度反演
由图8可看出,2006年12月20日冬季珠江口悬浮泥沙浓度整体平均浓度较低,处于0~152 mg/L范围,平均值76.34 mg/L,最大值152 mg/L。珠江口悬浮泥沙整体分布西高东低,珠江东侧支流较少且东江来水量小,导致虎门到入海口处悬浮泥沙浓度值较小,处于0~40 mg/L范围。珠江口西部支流众多,横门水道、洪奇门水道因受西江来水量大,加上两岸社会经济发达人口众多,最终携带泥沙含量较大堆积至入海口导致入海口附近悬浮泥沙浓度大,珠江口西部横门水道、洪奇门水道、东部交椅湾悬浮泥沙浓度值较高,处于0~160 mg/L范围,最大值为152.5 mg/L。东部交椅湾海岸由于工程施工,植被覆盖遭到破坏,悬浮泥沙受地表径流携带入海,导致交椅湾悬浮泥沙浓度较大。
2013年9月珠江径流量较冬季增多,珠江干流悬浮泥沙浓度整体平均浓度较高,处于0~242 mg/L范围,平均值125.8 mg/L,最大值242 mg/L。珠江干流较2006年12月中间段悬浮泥沙浓度值较高平均为131.3 mg/L,入海口悬浮泥沙浓度平均为48.1 mg/L。珠江口悬浮泥沙浓度的分布在东莞区域呈现出从近岸到离岸逐渐递减趋势。珠江口西部横门水道、洪奇门水道、东部交椅湾悬浮泥沙浓度值仍为范围内高值区域。秋季,珠江口呈现入海口和西部两片悬浮泥沙值较高区域,主要原因是秋季珠江口沿岸上升流减弱,西北季风还未开始,西南侧远离珠江口顶,致使悬浮泥沙浓度未能出现冬季仅西南侧高现象。而珠江口入海口因此季节收到台风、潮汐涌动影响,在珠江口形成高浓度值区域。
2020年4月处于汛期初,珠江口悬浮泥沙浓度整体平均浓度较高,处于0~248.2 mg/L范围,平均值139.4 mg/L,最大值248.2 mg/L。但此季节,珠江干流悬浮泥沙浓度值较高,浓度值均值为131.6 mg/L,深圳西侧海域悬浮泥沙呈现大面积浓度值较高现象,均值高达168.2 mg/L,主要受深中通道开通工程施工影响,悬浮泥沙较大。东部交椅湾因工程竣工植被恢复悬浮泥沙明显较2006年、2013年降低,项目工程施工对悬浮泥沙浓度影响非常大。珠江口西部横门水道、洪奇门水道受自然径流大影响悬浮泥沙浓度值仍较高,横门水道南支受上游输沙率减少悬浮泥沙浓度值也发生较大减少。
本文对比了机器学习人工神经网络、支持向量回归机及随机森林模型,发现随机森林模型在珠江口悬浮泥沙遥感监测模型精度最高。通过高光谱数据重采样到多光谱波段上,不同传感器构建对应的光谱响应函数,将实测高光谱数据重采样到对应传感器通道,能实现多传感器的悬浮泥沙遥感监测。采用模型精度最高的随机森林模型对珠江口2006年Landsat TM5、2013年Landsat ETM、2020年Sentinel-2进行了悬浮泥沙浓度反演,反演结果发现珠江口岸悬浮泥沙浓度呈现西高东低,从近岸到离岸逐渐递减的趋势。主要原因是珠江口地形为喇叭形,在季风和潮汐共同作用下,珠江口顶部区域受潮汐和风向混合作用强烈;在远离珠江口顶部区域地形较为开阔,各支流径流的悬浮泥沙汇聚堆积在珠江口西南区域,导致珠江口西南区域悬浮泥沙浓度高。基于地面实测站点数据对影像反演结果进行了评价,影像反演结果与实测数值相关系数较高,模型具有较高精度;同时发现港珠澳大桥的修建对珠江口悬浮泥沙的波动起着较大影响。