陈嘉星,张霄宇,2,3*,黄国容, 孙 杰
(1.浙江大学 地球科学学院,浙江 杭州310027;2.浙江省海洋观测—成像试验区重点实验室,浙江 舟山 316000;3.浙江大学 海洋研究院,浙江 舟山 31600;4.上海卫星工程研究所,上海 201109)
长江是我国的第一大河流,每年输入东海的悬浮泥沙含量可达1.81 亿t(据大通站2000-2011 年输沙量资料)[1],巨量的悬浮泥沙在口门堆积,造成长江口水下三角洲地形地貌的不停发展和演变,对于河口生态环境、航道安全、近岸工程等都具有重要意义。近年来,机载激光雷达技术的兴起为在河口近岸、岛礁以及船只无法到达的水域开展测深提供了新的解决方案。激光雷达同时具有受水深影响小、测量精度高、覆盖面广、测量周期短、测点密度大等优势,但同时存在着航空管制、费用昂贵等问题。因此,开展激光雷达测深的可行性预评估对于有效开展激光雷达测深具有重要意义。
然而,长江口水动力环境复杂,浮游植物旺发,造成了长江口复杂的水下光学场,给激光雷达的水下地形地貌探测带来很大的不确定因素。因此,利用卫星遥感的大面积快速观测能力,采用水体漫衰减系数Kd开展水下光场预评估,可为航空激光雷达测深提供先验知识。丁凯等[2]采用MODIS 的Kd(490)产品开展了对南海北部区域机载激光雷达测深的可行性评估,并指出了开展机载激光雷达测深的合理化时间节点。
水体漫衰减系数Kd是表征透射进入海水中的光能在海水及其悬浮物质的吸收作用和散射作用下,其下行辐照度随海水深度的增加而呈指数衰减的性质,水体漫衰减系数常用于估算水体深度、水体透明度以及真光层深度[3-11]。但是不同时间、不同区域水体的漫衰减系数呈现出一定的时空差异性,因此,选择合适的水体漫衰减系数反演模型是获得准确漫衰减系数的必要前提,也有助于深入研究水体的光学特征、光学类型及其对水生态系统的影响。
目前,已有不少学者针对长江口水体组分特征及其与吸收、散射之间的关系进行了大量研究,对长江口水体漫衰减系数的特性及其影响因素有了比较深入的理解,并以此为基础,构建了不同类型的漫衰减系数反演模型。如,王晓梅等[8]针对黄东海二类水体,利用490 nm、550 nm、670 nm 处的遥感反射比,建立了遥感反射比和漫衰减系数之间的经验统计关系,并利用同年春季实测数据进行检验,并指出该模型对季节变化具有较好的适应性。吴婷婷等[9]根据黄东海及珠江口附近海域的实测水下光谱剖面数据,研究发现漫衰减系数与光谱比值Rrs(555)/Rrs(490)存在分段线性关系,经相关性分析后加入665 nm 处遥感反射比,建立了分波段的Kd(490)反演算法。陈雨等[10]利用2012 年和2013 年冬季长江口邻近海域的实测数据,对漫衰减系数光谱特征、主要影响因子及空间分布进行分析,选取510 nm、590 nm、670 nm 处的遥感反射率建立了长江口冬季Kd(490)反演模型,并由此构建了长江口Kd(490)和Kd(PAR)的关系。Yu Xiaolong 等[11]利用2013 年3 月长江口现场实测数据,对双色模型进行了参数优化改进。因此,目前针对长江口漫衰减系数的反演算法较多,主要方法是基于船点实测数据开展相关敏感波段之间的关联分析。从应用角度来看,主要侧重于对长江口年度、季节变化特征和差异分析,由于构建模型采用数据集不同,使得不同的模型在应用上均有一定的局限性,且受卫星观测时像的影响,目前采用高时相卫星遥感影像进行大面积、动态的长江口漫衰减系数的日变化研究还相对较少。
有鉴于此,本次研究采用世界上第一颗静止水色卫星COMS(Communication Ocean and Meteorological Satellite)/GOCI(Geostationary Ocean Color Imager)为数据源,采用ENVI 自带的QUAC 方法,进行以下4 方面的研究:(1)选择目前研究区域已经采用的多种反演模型进行反演,从水体的性质对算法的适用性进行分析评估,选择适用于长江口的Kd(490)反演模型进行反演;(2)分析研究区域先涨潮再退潮期间Kd(490)的时空变化特征,并对其影响因素进行分析;(3)分析研究区域Kd(490)的季节变化特征,并对其影响因素进行分析;(4)基于Kd(490)的季节变化特征和潮汐变化特征,开展在长江口机载激光雷达测深能力评估。研究将系统地揭示长江口漫衰减系数的时空变化规律,对于精确反演透明度等水体光学参数具有重要意义,同时也为机载激光雷达测深提供背景信息。
长江口呈喇叭状向入海口展开,窄口端宽度远小于宽口江面宽度。距今2000 年伊始,长江河口南岸逐渐向海推进,河口北岸的部分沙岛相继并岸,口门宽度遂逐渐减小,河槽逐年加深,主槽向南偏,逐渐演变为多级分汊的三角洲河口。自1998 年以来,随着长江口北槽深水航道整治工程的建设,长江河口依托实体工程,以长堤的方式向海突出[12]。近年来,受长江口流域生态治理和三峡筑坝工程的影响,长江冲淡水流量、携沙量均呈现出逐年降低的趋势,同时长江口水下三角洲已经出现了明显的补偿冲刷现象,因此,动态观测长江口水深变化,对于研究自然和人类双重影响下长江口水下地形地貌的发展和演化及其生态环境等方面均具有重要意义。研究区域遥感影像来源于http://kosc.kordi.re.kr,如图1 所示。
图1 研究区域遥感影像图
长江口口外海域受苏北沿岸流、江浙沿岸流和台湾暖流的共同影响,存在着明显的季节差异。夏季,长江冲淡水流量大,在台湾暖流挤压下可向东甚至东北方向抵达达济州岛附近。冬季,长江冲淡水流量小,主轴流向为东南方向,只能由长江口到达杭州湾北部[13-14]。
长江口水富沙丰,大量泥沙自长江口入海,导致海岸带水体悬浮颗粒物浓度较高,陆域输入溶解有机物质丰富,加上浮游植物生长繁盛,导致水体光学特征复杂,是典型的二类水体[15]。长江口自徐六泾以下至杭州湾,潮流作用逐渐增强,径流作用逐渐减弱,悬浮泥沙浓度呈现不断增加的趋势。风浪、潮流冲刷引起的泥沙再悬浮,是长江口外悬浮泥沙浓度急剧增高的主要原因。悬浮泥沙浓度存在明显的季节差异,口门及口外夏低冬高,口内夏高冬低[16]。
2010 年,韩国发射了世界上第一颗地球静止气象卫星COMS,GOCI 为其上搭载的水色遥感器,波段设置为可见光到近红外共8 个波段,如表1 所示。时间分辨率可达1 h,每天可提供目标区域共8 个时间段(8:00-15:00)的观测数据,实现对以韩国为中心包括韩国、俄罗斯、日本、朝鲜及中国约6.35×106km2的目标区域的实时监测[17]。由于其较高的时间分辨率,减少了云层遮挡造成的误差与数据缺失的概率,对日间动态观测海洋和监测海洋突发事件具有突破意义。经过筛选,用于潮汐变化研究的数据共选取8 景,为2015 年8 月2 日的全天数据;用于季节变化研究的数据共选取4 景,拍摄时间如下:2015 年1 月19 日13:28;2015 年4 月21 日12:28;2015 年8 月2 日15:28;2015 年10 月11 日13:28。
表1 GOCI 波段及参数信息
国内外已有众多学者就漫衰减系数反演算法进行了大量研究,对比不同学者提出的Kd(490)反演模型,主要包括半分析算法和经验算法,如表2 所示。
半分析算法考虑辐射传输过程,加入吸收系数、后向散射系数等固有光学参量。如Lee 等[18-20]基于生物光学模型提出的QAA 半分析算法,此算法具有一定的普适性,但模型引入的参数较多且难以获得;Wang 等[21]针对切萨皮克湾提出的半分析算法,此算法可用于二类水体,但模型引入的固有光学量较难测量。
经验算法直接建立漫衰减系数与遥感反射比等表观光学量之间的统计学关系,如针对大洋一类水体,Austin 等[22]针对海岸带水色扫描仪CZCS 数据提出的经验算法,Mueller 等[23]针对第二代水色传感器SeaWiFS 提出的经验算法,这些算法主要是根据实测数据建立Kd(490)和蓝绿波段(443 nm、550 nm)的离水辐亮度比值之间的经验关系;针对近岸二类水体以及内陆富营养化水体的算法,金鑫等[24]采用2009 年6 月的巢湖水体的光谱数据、水质参数、上行辐亮度和下行辐亮度等实测数据,分析水体漫衰减系数光谱特征,建立了Kd(490)的多元线性回归模型。Moon J E 等[25]通过对比原位船载数据和GOCI卫星数据,认为在近岸高浊度海域,漫衰减系数与421 nm 和443 nm 处的遥感反射率Rrs的相关性较低,与490~865 nm 之间的遥感反射率Rrs的相关性较高。
另外还有,王晓梅等[8]、崔廷伟等[26]、吴婷婷等[9]、陈雨等[10]提出的经验算法,这些算法除了考虑到蓝绿波段的遥感发射率比值与Kd(490)之间的相关关系,同时引入红或黄波段的遥感反射率,这更加符合二类水体的特点。
表2 漫衰减系数的反演算法
考虑到长江口水体的浑浊特性:口内大部分水域的悬浮泥沙浓度在400~600 mg/L,高值区在镇海海域岸边可达1 000 mg/L,最低值10 mg/L 以下[27-28],且漫衰减系数与GOCI 波段中490~865 nm 之间的遥感反射率具有更高相关性的特点,因此本文选取吴婷婷等[9]提出的分段式漫衰减系数反演算法进行了长江口海域漫衰减系数的反演。
实验操作主要通过软件ENVI (ITT Visual Information Solutions 公司,5.3)和GDPS(GOCI Data Processing System,韩国卫星中心研发的GOCI 卫星配套的图像处理软件)进行,由于ENVI 目前还不能处理GOCI 提供的he5 格式的数据,所以采用GDPS将he5 格式的数据转换为ENVI 可处理的img 格式的数据。
GDPS(GOCI Data Processing System,韩国卫星中心研发的GOCI 卫星配套的图像处理软件)采用的大气校正算法是基于SSMM[29]改进的,在二类水体海域会出现过矫正,因此不适合应用于长江口邻近海域。ENVI 自带的快速大气校正QUAC(Quick Atmosphere Correction)则是一种不需要配套信息的大气校正方法[30-31],不需要考虑太阳天顶角、卫星的天顶角等光谱参数及大气气溶胶的厚度、当天的风速等气候学信息,而是基于图像自身的数据对其进行校正,即根据波段的中心波长信息自动从影像上收集不同物质的波谱信息,获得经验值从而完成快速大气校正。江彬彬等[27]将经QUAC 校正后的GOCI不同波段上的归一化离水辐射率的变化趋势与He等[32]采用紫外光谱校正方法以及实测值进行比较分析,GOCI 不同波段上的归一化离水辐射率的变化趋势与实测值一致,发现在490 nm、555 nm 处略小于实测值,但比值变化不大,在680 nm、745 nm 处略大于实测值,因此采用以490 nm、555 nm 处遥感反射比比值为界的分段式Kd(490)反演算法可能会造成低比值区Kd(490) 反演值偏低,高比值区Kd(490)反演值偏高,造成激光测深系统可探测深度近岸偏浅、远岸偏深。综合考虑运行效率以及作为激光测深的预评估的精度要求,本次研究最终采用了QUAC 法进行GOCI 遥感影像的大气快速校正处理。
除转换数据格式外,反演前处理还包括图像裁剪、快速大气校正、求地表反射率、NDWI 法提取掩膜,以便进行随后的漫衰减系数反演和激光雷达测深能力预评估。
我国近岸水体以二类水体为主,考虑到加拿大Optech 公司研制的CZMIL(Coastal Zone Mapping and Imaging Lidar)系统[33]专门针对浑浊浊水体研发了浊水处理算法模块(HydroFusion Turbid Water Module)。故本次研究选取CZMIL 系统的Kd(532)与最大测深之间的关系为例进行测深评估。
据文献[33-34],CZMIL 系统在海道测量模式下(白天),当底部反射率大于15%时,单脉冲最大测深深度Dmax(单位为m)与532 nm 处漫衰减系数Kd(532)(单位为m-1)之间的关系可近似为:
从式(1)可知,评估CZMIL 系统测深能力的漫衰减系数为Kd(532),而我们反演得到的漫衰减系数为Kd(490),故需将Kd(490)转换为Kd(532)。
CZMIL 系统采用Nd:YAG 激光器,其中用于探测水底信息的激光为蓝绿波段,波长为532 nm,该波段激光具有较强的穿透海水能力。因此分析CZMIL 的测深能力,首先要获得该波段处的漫衰减系数Kd(532)。当前国际上水色遥感的标准数据产品是Kd(490),故需要将Kd(490)转化为Kd(532)。根据陈雨等[35]关于长江口漫衰减系数的研究结果,400~700 nm波段范围内其他波段的漫衰减系数与Kd(490)的相关性在0.97 以上,故此波段内的漫衰减系数可以通过Kd(490)转化得到,两者之间的关系可以表示为:
式中:k(λ)=1.158×10-5λ2-0.01561×λ+5.8573;b(λ)=1.408×10-5λ2-0.01363×λ+3.3054。
由此,532 nm 处的漫衰减系数与490 nm 处的漫衰减系数之间的关系为:
3.1.1 长江口Kd(490)的日变化分布特征 对2015年8 月2 日的8 景GOCI 遥感影像进行长江口及其邻近海域的Kd(490)日变化特征分析,结果显示如图2。总体来看,每个潮位的Kd(490)最高值均在舟山本
岛附近,最高值在(2.8±0.2)m-1,最低值在舟山本岛的南东侧东海海域,最低值基本稳定在(0.10±0.02)m-1。从河口到海,随着离岸距离的增大,基本呈现出先增大再减小的趋势。
图2 2015 年8 月2 日全天长江口及邻近海域Kd(490)反演
本文选取的2015 年8 月2 日的8 景影像在拍摄时间内处于先涨潮再退潮期间。涨潮期间,长江口和杭州湾的水流方向近于SE-NW;退潮期间,长江口和杭州湾的水流方向近于NW-SE。涨潮期间,长江口及其邻近海域的Kd(490)值总体呈现逐渐增大的变化趋势。退潮期间,长江口及其邻近海域的Kd(490)值总体呈现逐渐减小的变化趋势。
结合表3 及长江口悬浮泥沙浓度日变化的相关文献进行分析[27],长江口的漫衰减系数与其他浑浊的河口二类水体的漫衰减系数相当,并且其主要影响因素均为悬浮物。3.1.2 长江口Kd(490)日变化分布特征的影响因素分析 如图2 所示,Kd(490)反演值分布图表明:Kd(490)值随着离岸距离的增大呈现出先增大再减小的趋势,这与研究区域内悬浮泥沙浓度分布规律[27]基本吻合,舟山岛附近出现了漫衰减系数Kd(490)的最大值(2.8±0.2)m-1,此处悬浮泥沙浓度最大,这从侧面说明了悬浮泥沙是Kd(490)的较大影响因素之一。
表3 其他水体的Kd(490)分布
涨潮期间,长江口和杭州湾的水流方向近于SE-NW;退潮期间,长江口和杭州湾的水流方向近于NW-SE,初期杭州湾内漫衰减系数Kd(490)值空间变化相对较小,但随着潮水逐渐退潮,漫衰减系数Kd(490)高值区沿着杭州湾顶逐渐由湾内转移到舟山本岛上,这与研究区域内悬浮泥沙浓度在落潮期间的变化规律[27]很相似,因此有理由认为潮水的潮位通过影响长江口及其邻近海域的悬浮泥沙浓度日变化,进而影响漫衰减系数Kd(490)的日变化。
以往有关研究在采用遥感技术分析我国近岸漫衰减系数的季节变化时,常常采用月平均值对漫衰减系数的季节变化特征进行研究,然而,一个潮周期内可造成漫衰减系数的显著变化,因此,采用月平均值无法对研究区域内水动力环境的动荡变化进行精细描述。本次研究选取2015 年1 月19日、4 月21 日、8 月2 日、10 月11 日各一景共4 景GOCI 影像,这4 景影像的潮位相近,故而可以在相同潮位下研究长江口及其邻近海域的漫衰减系数分布和变化的特征及其影响因素,反演结果如图3所示。
图3 2015 年四季典型潮位长江口及其邻近海域Kd(490)反演
如图3 所示,结果表明,长江口及其邻近海域Kd(490)值夏季(8 月)低于冬季(1 月),春季(4 月)和秋季(10 月)居中,这种变化特征与研究区域内的季节性流量变化密切相关。8 月正值夏季,是长江的丰水期,因此长江冲淡水流量很大,长江冲淡水携带着大量悬浮泥沙冲出长江口外,这个流动过程使得夏季水体中悬浮泥沙被稀释导致其含量减小,而冬季长江冲淡水流量减小,同时受季风影响,底部沉积物有可能受到风浪扰动[13-14],使得冬季水体中悬浮泥沙含量增大,长江口悬浮泥沙含量与冬季相比较低,且相同区域的漫衰减系数Kd(490)值夏季较冬季[10]要小,因此可以认为长江冲淡水流量通过影响长江口及其邻近海域悬浮泥沙含量的季节性变化,进而影响漫衰减系数Kd(490)的季节性变化。
首先,长江口及其邻近海域的Kd(490)值夏低冬高,春秋居中。因此,只需要选择夏季和冬季的两幅影像对CZMIL 测深能力进行评估,评估结果如图4所示。
图4 长江口冬夏两季测深能力空间分布
由图4 可知,夏季较冬季更适合CZMIL 系统进行测深,长江口及杭州湾湾内水域可探测深度在约在5 m 以下,口外清洁水体可探测深度在22 m 左右,探测季节以夏季为宜。
长江口及其邻近海域的Kd(490)值总体呈现出随潮位的增高而增大,随潮位的减小而逐渐减小的变化趋势。因此,只需选择夏季潮位最高和退潮低潮位的两幅影像对CZMIL 测深能力进行评估,评估结果如图5 所示。
图5 长江口夏季最高潮位和退潮低潮位测深能力空间分布
将本次研究结果和李凯等[38]关于在黄海、东海区域采用CZMIL 系统开展的测深能力评估等深线图比对,基本一致:长江口及其邻近海域可探测深度在0~22 m 左右,长江口及杭州湾内可探测深度约在5 m 以下,口外清洁水体可探测深度在22 m左右,探测时间选择退潮低潮位的时刻为宜。
(1)先涨潮后退潮期间,长江口及其邻近海域的Kd(490)值总体表现为先增大后减小,潮位是影响长江口及其邻近海域水体的Kd(490)值日变化的主要因素。
(2)结合长江口的悬浮泥沙浓度日变化的研究,Kd(490)的分布规律和变化特征与悬浮泥沙浓度的变化规律基本一致,说明悬浮泥沙是Kd(490)的较大影响因素之一。
(3)长江口及其邻近海域的Kd(490)值的季节变化特征表现为冬高夏低,春秋居中,且影响其季节变化的主要因素是长江冲淡水流量和季风。
(4)根据反演结果,夏季是最适合开展CZMIL测深系统作业的季节,退潮低潮位比其它潮位更适合开展作业。CZMIL 测深系统在长江口及其邻近海域可探测深度在0~22 m 左右,长江口及杭州湾内可探测深度约在5 m 以下,口外清洁水体可探测深度在22 m 左右。
研究表明,GOCI 8 景/d,1 景/h 的分辨率可以实现对不正规半日潮潮周期内Kd(490)的动态变化监测,而且可以实现在相同潮位下更为合理地描述Kd(490)值的季节变化,为进一步开展机载激光雷达探测提供水体环境观测技术保障和支持。
今后将进一步开展长江口及邻近海域内全潮周期、月度以及年度的漫衰减系数变化规律研究,并对期间的机载激光测深可行性进行较系统的评测。