邓智天 孙永华 邱 琦 孙 薇 倪 萍 邢 瑞
(首都师范大学资源环境与旅游学院,北京 100048)
悬浮物浓度(suspended sediment concentration, SSC)是海水最重要的参数之一[1].不同浓度的悬浮物直接导致光散射、折射和吸收的差异,最终导致水的光学性质的差异.超过一定浓度的悬浮物会使海水浑浊,水的透明度降低,从而导致入射到水中的太阳辐射减弱,海水中浮游生物吸收的能量减少,使得光合作用被削弱,并最终影响海水中浮游生物的生长.悬浮物对水质的影响也较大,原因在于悬浮物包括一定数量的胶体物质和黏土矿物.一方面,悬浮物具有一定的吸收能力,可以吸收污染物,起到净化水的作用;另一方面,它成为污染物的迁移和再循环的重要依附物.同时,悬浮物也是研究泥沙输移和地貌演化的重要依据,可以帮助预测沉积物的承载能力,并在沿海水域泥沙输移规律的研究中发挥重要作用.其研究不仅可以探索海水污染的情况,还可以进一步追踪污染物的流动路径和海水中污染物的方向,然后探讨海水中污染物的运动路径和变化规律,为海水污染控制提供技术支持.
遥感可以有效地克服现场测量的不足,具有范围广、效率高的优点,成为了研究海水中悬浮物的主要技术手段.20世纪90年代初,Ritchie和Cooper[2]和Topliss等[3]使用Landsat 3 多光谱扫描仪(multispectral scanner, MSS)进行了水体悬浮物的浓度反演,指出MSS对悬浮物浓度有相关性;Landsat 5在轨稳定运行之后,很多研究者对专题绘图仪(thematic mapper,TM)数据进行了水体悬浮物的研究[4-12],使用了单波段、波段组合等方法构建模型,对水体悬浮物进行了定量研究;还有学者使用了MODIS数据对更大范围的湖泊、海洋等进行了研究[13-19],由于其波段较多且数据较多,包括卡尔曼滤波、机器学习等方法相继应用于水体悬浮物的定量研究;相关学者分别对海景宽视场传感器[20]、机载高光谱扫描仪[21]、小型机载光谱成像仪[22]、中分辨率成像光谱仪[23]、伊科诺斯[24]等水色遥感、高光谱遥感、高分遥感传感器进行悬浮物反演,沿用了MODIS和Landsat的反演方法,并使用了更为专业的传感器,因此在反演精度上有所提高.
QAA由Lee等[25]在2002年提出,用于推导出深水光学的固有光学特性.固有光学特性(inherent optical properties,IOP)是指仅与水体自身组成成分相关,由介质决定,不随光照条件改变而改变的水体光学参量,主要包含海水的吸收系数、散射系数和衰减系数等[26],其相关研究集中在IOP与悬浮物颗粒的特性[27].国际海洋水色遥感组织在2006年提交了IOP算法报告[28],在该报告中,表明QAA算法对于清澈海水和沿海水域都有非常高适用性.直到2018年,QAA已经历了5次更新,最新正式发布的版本是QAAv6[29].这些更新使QAA更加完美,性能也大大提高.在Lee之外,许多学者对算法进行了改进,如Chen和Zhang[30]基于MODIS数据,提出了一种改进的算法QAA-RGR, 直接选择卫星遥感数据而非实测高光谱,证明卫星遥感数据在一定条件下能替代实测高光谱进行QAA算法的计算,并能取得较高的精度,这也开启了近几年使用不同的遥感数据关于QAA算法的研究.
近年来,基于QAA算法悬浮物的研究有所增加.陈莉琼[31]使用实测高光谱进行研究,指出QAA模型比经验模型具有更高的精度和时相稳定性;王建国等[32]、段化杰等[33]分别使用MODIS和Sentinel 3 OLCI数据对悬浮物浓度进行了相关研究,发现其精度与实测高光谱相当.
通过QAA算法,可以实现由遥感反射率数据推导IOP,既减少了IOP数据获取的成本,具有时效性与经济性,又可以对遥感反射率数据进行深层次的数据挖掘,而非局限于反射率数值.而研究卫星多光谱数据QAA算法适用性,探讨其能否在一定条件下替代现场实测高光谱数据,可以进一步放大QAA算法的时效性与经济性.
本研究将以含沙量较高的辽河口三角洲为例,基于QAA算法,使用天宫二号宽波段成像仪数据,选取682.5 nm为参考波长,求出水体总吸收系数a(λ),利用2018年6月7—9日实测水体悬浮物浓度数据,建立基于a(λ)的悬浮物浓度反演模型.
辽河平均径流量95.27亿m3,河道弯曲,呈不规则河型,水系发育,大小支流70余条,中下游河道宽浅,河道宽1 000~2 000 m;铁岭水文站多年平均含沙量3.60 kg/m3,多年平均输沙量2 098万吨;水流缓慢,泥沙淤积;河床质为沙壤土.治理水土流失和开展水土保持也是辽河流域的一个主要问题,50年代初辽河流域水土流失面积约1 220万hm2,占流域面积的55.7%;1985年全流域水土流失面积为951万hm2,占全流域水土流失面积的45.2%;其中老哈、教来河位于冀北辽西山地和黄土丘陵区,植被覆盖率不到30%,水土流失十分严重.
样本采集于6月7—9日沿着双台子河和大辽河进行,如图1所示,图例中绿色三角形为采样位置,底图为全球海图数据,两条河流属于辽河支流.本研究采样方案着重在河道中间深水处,避开滩涂,即图中绿色部分(水深≤0 m);由于海船体积较大,吃水深,行驶时可能会带动底层沉积物及搅浑海水,本研究采样时采用线性采样,沿着行驶的路线等间距采样,并在驶入采样区时逐渐减速,以减少底层沉积物带来的误差;北边采样重点在双台子河的悬浮物浓度;中间部分采样重点在河口三角洲滩涂附近的悬浮物浓度,并在大潮时(水深>1.5 m)对滩涂正上方进行了数次采样;东边采样重点在大辽河的悬浮物浓度.本次样本采集一共采集104个水样进行测量.
试验标准参照中华人民共和国生态环境部确定悬浮物浓度标准《水质 悬浮物的测定 重量法:GB 11901—89》[34].该标准指出水质中的悬浮物是指水样通过孔径为0.45 μm的滤膜,截留在滤膜上并于103~105 ℃烘干至恒重的固体物质.
(1)用扁咀无齿镊子夹取微孔滤膜放于事先恒重的称量瓶里,移入烘箱中于103~105 ℃下烘烤半小时,烘干半小时后取出置干燥器内冷却至室温,称其重量.反复烘干、冷却、称量,直至两次称量的重量差<0.2 mg;
(2)量取充分混合均匀的试样100 mL抽吸过滤;使水分全部通过滤膜;再以每次10 mL蒸馏水连续洗涤三次,继续吸滤以除去痕量水分;
(3)停止吸滤后,仔细取出载有悬浮物的滤膜放在原恒重的称量瓶里,移入烘箱中于103~105 ℃下烘干1小时后移入干燥器中,使冷却到室温,称其重量;反复烘干、冷却、称量,直至2次称量的重量差<0.4 mg为止.
图1 研究区及采样位置示意图
天宫二号空间实验室,是继天宫一号任务完成其使命后,发射的第二个空间实验室,安排了一批体现科学前沿和战略高技术发展方向的科学与应用任务,开展相应的应用与试验,对相关新技术进行体制验证.天宫二号空间实验室发射后,与神舟十一号载人飞船对接,进行人在太空的中期驻留实验;与天舟一号货运飞船进行对接,开展推进剂补给等相关试验.
天宫二号空间实验室搭载了全新的空间应用载荷设备,载荷数量及规模都超过了以往我国各次载人航天任务.开展10余项应用与实验,涉及对地观测和空间地球科学、空间天文、微重力基础物理、微重力流体物理及空间材料科学、空间生命科学和空间环境与空间物理等多个领域.其中对地观测和空间地球科学有宽波段成像仪、三维成像微波高度计和紫外临边成像光谱仪,其中宽波段成像仪的波段设置如表1.
表1 天宫二号宽波段成像仪波段设置
本研究使用影像获取时间与采样时间间隔最小的一幅TG-2数据,2018年6月1日的数据,由于2018年6月初无强降雨台风等天气,故认为2018年6月1日的影像能反映悬浮物在2018年6月7—9日采样时的实际分布.
TG-2 MWI Level-2数据通过载人航天空间应用数据推广服务平台申请数据产品(http://www.msadc.cn/),该Level-2数据产品已经经过几何校正和辐射纠正.Level-2数据包含表观反射率的校准数据,图像的表观反射率由ENVI中的大气校正模块Quick Atmospheric Correction(QUAC)生成.
为去除水面折射对反射率的影响,通过公式(1)将表观反射率Rrs(λ)换算为水面以下反射率rrs(λ):
(1)
后向散射与吸收系数比值u(λ)和总吸收系数a(λ)及后向散射系数bb(λ)的定量关系如公式(2):
(2)
其中,g0=0.089,g1=0.1245.
表2 QAAv6算法
悬浮粒子后向散射系数bbp(λ)和波长λ满足幂函数规律,如公式(3):
(3)
其中λ0代表参考波长.
在海洋中,纯水后向散射bbw(λ)占bb(λ)的1/10[35],如公式(4):
(4)
根据悬浮物反射特征,700 nm处的吸收系数迅速下降,我们将682.5 nm指定为λ0,即VNI Band 7的中心波长.
QAAv6算法步骤见表2,步骤①将水面的表观反射率转换到水面以下的反射率;步骤②计算后向散射与吸收系数比值,经过查找天宫数据元文件,Rrs(670)>0.001 5sr-1,执行右边的公式;步骤③~④为参考波长的固有光学参数估算,其中纯水吸收系数aw(670)根据Deng等[36]所做的室内纯水吸收系数表格查找得出;步骤⑤~⑥为任意波长的固有光学参数估算,通过计算的λ0与λ的系数,估算出任意波长的后向散射系数bbp(λ);步骤⑦高光谱波长的总吸收系数估算,通过前面计算得出的u(λ)、bbw(λ)、bbp(λ),估算总吸收系数a(λ).本研究使用线性回归模型建立a(λ)和SSC之间的反演模型.
本研究构建悬浮物浓度反演模型主要有2个部分:第一部分是输入遥感反射率为,利用QAA估算水体中总吸收系数a(λ);第二部分是根据QAA算法估算的水体中总吸收系数a(λ)构建波段归一化因子与悬浮物浓度值建立回归模型,进而求解水体悬浮物浓度.
悬浮物浓度测定结果见图2,大部分样品浓度位于18~33 mg/L的区间内,其中双台子河口采样区域43个点均值为23.91 mg/L,大辽河口采样区域35个点均值为24.66 mg/L,河口三角洲采样区域26个点均值为26.08 mg/L;有4个离群点,应为样本保存或实验过程中的误差,本研究将该4个离群值去除.
图2 悬浮物浓度测定结果
除去异常点,在该研究中使用100个点;随机选择70个点进行回归分析,其余30个点用于验证.如图3所示,原始光谱曲线在700 nm附近有一个反射峰,这与悬浮物的光谱特性一致.
图3 样本表观反射率Rrs(λ)
bbp(λ)与Rrs(λ)和a(λ)相比,由于计算方法的特性,仅需要一个参考点λ0便可以计算出所需的400~1 000 nm内所有波段,如图4所示,其光谱特征具有一致性,后向散射系数随着波长的增加呈幂函数递减.
图4 样本后向散射系数bbp(λ)
在QAA算法模型的第⑦步之后,得到了总样本的吸收系数,结果如图5所示.此处由于涉及到原始多光谱的相关波段,所以估算出来每一个样品的a(λ)只有14个点,后文a(565)代表着565 nm处的a(λ).
图5 样本总吸收系数a(λ)
在构建模型时,使用显著性检验来辅助波段选择,如图6所示,选择超过显着性水平0.01的波段,即图中黑线外的点.
图6 a(λ)与SSC的相关系数
根据悬浮物的反射特性,本研究选取Band9和Band10波段的水体吸收系数a(620)和a(565)构建波段归一化因子,并将该归一化因子与悬浮物浓度值建立线性回归方程:
(5)
本研究将验证组的30个点的相关数据,通过QAA模型进行演算,对比其模拟值与实测值,评估总吸收系数的反演效果并分析QAA模型准确性.本研究使用均方根误差(average relative error,简称RMSE)[37],平均相对误差(average relative error,ARE),如公式(6)和(7),以及相关系数(R2)来评估.
(6)
(7)
QAA模型的准确度R2为 0.685 9,RMSE为2.814 8,ARE为9.68%.通常,QAA估算得到的a(λ)与绿色(500~550 nm)和红色(600~650 nm)波段一致,但在其他波段,尤其是蓝色波段(400~500 nm),估算结果较差且明显被低估,而在550~600 nm,则有过高估算[38].对于大多数内陆浑浊水域,浮游植物色素吸收系数aφ(λ)和黄色物质吸收系数ag(λ)是相互独立,而海洋或沿海水样中的水质参数浓度或吸收系数相对低于内陆浑水,aφ(λ)和ag(λ)之间存在相关关系.以往的研究表明,ag(λ)在水中呈现出指数衰减规律,导致短波段遥感反射中出现大量复杂的水色信息,从而增加了短波段算法的不确定性.
本研究利用2018年6月1日的TG-2 MWI数据,基于QAA算法,选择682.5 nm作为参考波长λ0,建立了适用于了河口的悬浮物浓度反演模型,将该模型应用于TG-2 WIS数据,其结果如图7所示.
图7 辽河三角洲悬浮物浓度空间分布
从图7可知,悬浮物浓度呈现自北向南逐渐降低的趋势,浓度范围11~32 mg/L,双台子河、大辽河河口浓度都大于 20 mg/L,与实测值相近;码头静止的水的悬浮物浓度较码头外低了10 mg/L 左右,符合实际情况.
第一个误差是由于大多数采样点靠近岸边,其浓度值高于研究区域的其他部分,导致反演结果相对集中在高值.第二个误差在于QAA算法是根据水的辐射度逐步建立光学和生化特征之间的关系,其计算步骤容易传输并累积误差.第三个误差主要来自水面以上的遥感反射转换过程,计算无法直接获得的水下反射率rrs(λ),目前,对不同的测量角度进行精确建模需要考虑风速、阴影等,而本方法使用经验参数来替代.此外,即使排除河床底部反射并假设其为零误差,QAA算法的计算仍存在许多不确定性.如在QAA算法的步骤1中,g0和g1不是常数,但受到太阳天顶角和水体光学特性的影响,并且随水体散射特性而变化,所以参考波长处的吸收系数也是QAA算法中的主要误差源之一.
本研究以含沙量较高的辽河口三角洲为例,利用2018年6月7—9日实测水体悬浮物浓度数据,基于QAA算法,使用天宫二号宽波段成像仪数据,选取682.5 nm为参考波长,求出水体总吸收系数a(λ),建立基于a(λ)的悬浮物浓度反演模型,得到2018年6月辽河三角洲悬浮物浓度空间分布图.主要结论包括:
(1)根据悬浮物反射特征,700 nm处的吸收系数最低,将682.5 nm指定为λ0是合适的;(2)基于遥感反射率与吸收系数推出的总吸收系数a(λ)-悬浮物模型的拟合精度R2为0.685 9,能较为正确的反映辽河三角洲悬浮物浓度变化趋势,但是精度有待提高;(3)QAA算法强大但复杂,为了减少其误差,需要在采样时留意距岸距离、天气等因素,才能更好发挥其作用.
本研究探讨了天宫二号宽波段成像仪数据通过QAA模型演算IOP的适用性,通过将遥感反射率转化为总吸收系数,研究总吸收系数与悬浮物浓度的相关性,选取最佳波段,总结反演模型.受限于作者水平,本研究IOP演算止于总吸收系数,如果能从总吸收系数中分离出悬浮物单位吸收系数会有助于提高反演模型的精度.