海洋水色卫星遥感器接收的总信号中大约85%以上的信号来自瑞利散射和气溶胶散射, 而包含水体信息的离水辐射信号不足15%, 即使是很小的大气校正误差也能引起很大的水质参数反演误差, 导致大多数光学遥感图像的可利用性低。因此, 消弱这些因素的影响, 即对获得的影像进行大气校正显得尤为重要。如何将水体信息同大气信息分离, 是水色遥感得以应用的关键^[1-4]。

大气校正是要从混浊水域中把气溶胶和海洋贡献分离。在国内外的水色遥感研究, 前人已经提出了多种大气校正方法。其中, 基于辐射传输理论的有MOTRAN和6S等^[5-7]算法。然而此类方法需要实时的大气参数数据。因此, 基于遥感影像特定波段的算法得到很好的发展。针对这一类水体Gordon等^[8]把两个或两个以上的近红外波段的离水辐亮度假定为零, 提出海洋大气校正算法。对清澈水近红外反射率零是有效的^[9]。然而浑浊水体近红外波段的反射率并不为零, 高浊度海洋辐亮度在近红外波段不是线性关系^[10-11], 需要准确的近红外波段的光谱之间关系来建模^[12-13]。因此该算法在二类水体中得到了进一步的优化, 并出现了两种假设: (1)浑浊水体在近红波段的反射率的比值为常数; (2)浑浊水体在短波红外波段的反射率为零。在GOCI领域内Wang等^[14]已开发高混浊水的模型。海洋水色传感器的短波红外波段和现有的短波段有显着的噪音^[15-17]。在SeaDAS水域处理中非零的近红外反射率模型已运用在低到中度混浊水域^[18-19]。

目前Vanhellemont和Ruddick^{[20, 21]}针对Landsat-8 OLI发展了基于短波红外波段的一种大气校正算法, 并用中分辨率成像光谱仪(MODIS)进行了检验, 已经在比利时的沿岸浑浊水体和泰晤士河口水体中得到了验证。青松等^[22]已验证基于Landsat8 OLI短波红外的Vanhellemont和Ruddick的算法在乌梁素海的适用性。

黄河口位于渤海的西部, 受黄河大量泥沙注入影响, 其光学性质非常复杂, 水体也极度浑浊。因此本文基于Vanhellemont和Ruddick的算法, 进行黄河口Landsat-8 OLI数据的大气校正, 进而检验该算法在浑浊的海洋沿岸水体中的适用性。在此基础上, 利用时间序列OLI数据反演黄河口水体的悬浮物浓度, 并结合输沙量数据和海面风场数据分析其时空分布特征和影响因素。

1 研究数据 1.1 研究区

黄河是我国最有代表性的入海河流之一, 并且黄河是世界上泥沙含量最高的河^[23]。本文共获取黄河口海域16个站位的光谱数据, 采样时间为2015年6月3—11日, 并将其作为研究区大气校正结果的检验数据, 站位分布见图 1。

1.2 实测数据

利用ASD(Analytical Spectral Devices)公司生产的便携式双通道可见光—近红外光谱辐射计Field SpecPro Dual VNIR获取水面遥感反射率。ASD光谱仪波长覆盖范围为350~1 050 nm, 采样间隔为1.4 nm, 记录数据间隔为1 nm, 光谱分辨率为3 nm, 波长精度为±1 nm。它配有两个512通道的分光计, 可同时测量目标物体和参考物体的光谱。

遥感反射率光谱如图 2所示。可以看出光谱曲线形状很相似, 当悬浮物浓度处于不同区间时, 水体光谱具有明显的差异性。现场观测光谱中发现这16个站点光谱曲线已经分成明显的两部分, 上半段五个曲线为黄河口处的9、10、11、12和13号站点, 这五个曲线反射率明显高于其他点, 其在黄河口输沙口, 还有750~900 nm期间这两组曲线区别明显, 可能是因为下半部分对应的区域水域较混浊。

从网上下载得到2013—2015年的黄河利津水文站的月输沙量(黄河泥沙公报网http://www.yellowriver.gov.cn/nishagonggao/)和2013~2015年CCMP卫星遥感海面风场数据(https://www.remss.com/measurements/ccmp), 分析黄河口海域悬浮物浓度分布特征。该海面风场数据空间分辨率为0.25°, 时间间隔为6 h, 给出的是距海面10 m处风场沿经向和纬向的速度分量。

1.3 遥感数据

本文所用的遥感数据为黄河口海域2013—2015年7景(121/34)Landsat-8 OLI影像, (https://glovis.usgs.gov/), OLI数据波段设置见表 1。本文还下载了2013—2015年7景对应的中分辨率成像光谱仪(https://ladsweb.modaps.eosdis.nasa.gov/search/)MODIS数据, 用于大气校正结果的检验。

2 方法

本文利用Vanhellemont和Ruddick^[20-21]的方法, 首先利用Landsat-8 OLI数据对其进行辐射定标、瑞利校正和气溶胶校正即完成大气校正, 获得Landsat-8 OLI数据遥感反射率, 并对MODIS数据采用同样的方法获得MODIS数据的遥感反射率, 地面实地调查获得实测水体光谱数据, 采用地面实测数据和MODIS数据来验证大气校正结果, 最后反演与分析黄河口悬浮物浓度空间分布特征, 并采用输沙量数据和海面风场数据分析悬浮物浓度的影响。

2.1 大气校正

本文参考Vanhellemont和Ruddick的方法基于黄河口区域Landsat-8 OLI、MODIS影像进行大气校正。主要方法如下所述。

2.1.1 大气层顶反射率

首先, 利用如下公式计算遥感影像的辐射亮度值:

$ {L_{{\rm{TOA}}}} = a{D_{\rm{N}}} + b $

(1)

式中, L_TOA为大气顶层辐射亮度, D_N为图像灰度值, a和b分别为定标参数, 由OLI数据的头文件中获取。

其次, 基于公式(2)将辐射亮度转换为大气层顶的反射率:

$ {\rho _{{\rm{TOA}}}} = \frac{{\pi {L_{{\rm{TOA}}}}{d^2}}}{{{F_0}\cos {\theta _0}}} $

(2)

式中, d是日地距离天文单位, F₀为波段平均太阳辐射, θ₀为太阳高度角。

2.1.2 瑞利校正

水体的大气顶层反射率可以表示为瑞利反射率、气溶胶反射率和离水反射率之和。

$ {\rho _{{\rm{TOA}}}} = {\rho _\gamma } + {\rho _{\rm{a}}} + t{\rho _{\rm{w}}} $

(3)

式中, ρ_γ是瑞利反射率, ρ_a是气溶胶反射率, ρ_w是离水反射率, t为大气透射率。

利用如公式(4)计算瑞利反射率:

$ {\rho _\gamma } = \frac{{{T_\gamma }{P_\gamma }({\theta _0}, {\theta _{\rm{v}}}, \Delta \varphi )}}{{4\cos {\theta _0}\cos {\theta _{\rm{v}}}}} $

(4)

式中, T_γ是瑞利光学厚度, P_γ是瑞利散射相位函数, θ₀和θ_v分别是太阳和卫星的天顶角, Δφ是太阳和卫星的相对方位角。

其次, 基于公式(5)计算瑞利校正后的反射率:

$ {\rho _{\rm{c}}} = {\rho _{{\rm{TOA}}}} - {\rho _\gamma } $

(5)

2.1.3 气溶胶校正

Vanhellemont和Ruddick的算法的核心是气溶胶校正。此算法有两个假设: (1)近红外波段的ρ_w值为零, 因此在近红外和短波红外波段处ρ_a=ρ_c; (2)近红外和短波红外波段的气溶胶反射比ε为常数, 可表示为:

$ {\varepsilon ^{(S, L)}} = \frac{{\rho _{\rm{a}}^{(S)}}}{{\rho _{\rm{a}}^{(L)}}} = \frac{{\rho _c^{(S)}}}{{\rho _{\rm{c}}^{(L)}}} $

(6)

式中, S和L代表OLI数据的近红外和短波红外波段, 对应于5、6和7。

由ε^{(S, L)}可以得到第i波段的ε:

$ {\varepsilon ^{(S, L)}} = {({\varepsilon ^{(S, L)}})^{{\sigma _i}}} $

(7)

式中, $ {\sigma _i} = \frac{{L - {\lambda _i}}}{{L - S}} $。

然后第i波段的气溶胶反射率可由式(8)计算得到:

$ \rho _{\rm{a}}^i = {\varepsilon ^{(i, L)}}\rho _{\rm{a}}^L $

(8)

将公式(7)和公式(8)代入公式(5)中, 可以得到OLI数据各波段的离水反射率:

$ \rho _{\rm{w}}^i = \frac{1}{{{t^i}}}[\rho _{\rm{c}}^i - {({\varepsilon ^{(S, L)}})^{{\sigma _i}}}\rho _{\rm{a}}^L] $

(9)

2.2 悬浮物浓度反演

本文利用单波段算法来反演黄河口水体的悬浮物浓度, Nechad等^[24]提出一种基于多传感器数据的浑浊水体悬浮物浓度反演算法, 结合光谱反射率模型, 以实测的欧洲南北海海区水体遥感反射率与悬浮物浓度数据为基础, 拟合得到了适用于SeaWIFS、MODIS、MERIS等多水色卫星传感器的单波段悬浮颗粒物浓度半分析反演算法, 并在欧洲北海海区得到良好的真实性检验效果。基于此算法Vanhellemont和Ruddick等^[20-21]已发展出基于Landsat OLI、MODIS数据的单波段反演算法。

$ {C_{{\rm{SPM}}}} = \frac{{A{\rho _{\rm{w}}}}}{{1 - {\rho _{\rm{w}}}/C}} $

(10)

式中, C_SPM为悬浮物(Suspended particulate matter, SPM)浓度, ρ_w为离水反射率。针对Landsat-8 OLI的655 nm波长使用, A=289.29 g/m³, C=0.168 6, 对于MODIS/Aqua的645 nm使用, A=258.85 g/m³且C=0.164 1。

2.3 精度检验

基于MODIS得到的SPM结果检验了Landsat OLI反演结果, 精度检验使用平均相对误差公式:

$ e = \frac{1}{N}\sum\limits_{i = 1}^N {\frac{{|{x_i} - {y_i}|}}{{{x_i}}}} \times 100\% $

(11)

式中, x_i为Landsat OLI数据SPM, y_i为MODIS数据SPM, N为样本数。

3 结果与讨论 3.1 大气校正算法验证

为了检验大气校正算法的适用性, 对由Vanhelle mont和Ruddick算法获取的OLI反射率(5日成像)与实测离水反射率(4、5、8、9日获取)进行对比(见图 3), 发现大气校正后OLI数据中的反射率与实测值之间较一致, R为0.95, 在483 nm(蓝)、561 nm(绿)、655 nm(红)波段处平均相对误差分别为16.4%、17.3%和25.7%。误差存在的主要原因是实测时间与卫星过境时间存在差异, 水体性质也不同, 将会直接影响表现光学性质, 导致水体反射率的数值发生变化, 但水体反射率的光谱形状不会发生很大变化, 说明本文所用的Vanhellemont和Ruddick的方法适用于对黄河口OLI数据大气校正。

此外, 对由Vanhellemont和Ruddick算法获取的OLI数据蓝、绿、红波段反射率与MODIS数据蓝、绿、红波段反射率进行对比(见图 4), 发现大气校正后OLI数据中的反射率与MODIS数据反射率之间的一致性较好, R分别为0.89、0.91、0.90, 在蓝、绿、红波段处平均相对误差为26.8%、17.8%和18.7%。

对基于OLI数据和MODIS数据计算得到的SPM结果进行对比(见图 5), 发现OLI数据SPM与MODIS数据SPM之间的一致性较好, R为0.85, 平均相对误差为20.5%。

MODIS反射率和SPM存在误差的主要原因是研究区的时间动态观察条件和传感器设计的差异, 对MODIS来说由于观测天顶角较大, 导致图像中的几何变形也很大, 以至于造成与Landsat-8 OLI数据不匹配。这空间错配引入了相当大的错误, 尤其是分辨率不同的传感器(OLI: 30 m, MODIS: 250 m)。此外, 由于MODIS的扫描系统在扫描线之间有显著的重叠和他们之间也存在时间差, 可以造成表面SPM浓度的显着差异。Landsat-8 OLI和其他高分辨率的传感器可以被用来监测海上风力发电机组引起的混浊。悬浮物浓度的定量监测, 红色和近红外波段的气溶胶校正需要信噪比足够高的传感器。说明本文所用的Vanhellemont和Ruddick的方法适用于对黄河口水体OLI悬浮物浓度反演。

3.2 悬浮物浓度反演

将SPM模型(公式(10))应用于2013—2015年黄河口大气校正结果图, 获取了悬浮物浓度时空分布图, 结果如图 6所示。大致反映了春夏秋冬黄河口附近SPM分布变化的规律, 整个黄河口附近的海面, SPM小于60.0 g/m³:冬季 > 夏季, 夏季SPM最低, 这与许多学者^[25-27]研究东海陆架泥质区悬浮体含量有季节性变化(冬季悬浮体含量远高于其他季节, 夏季悬浮体含量低)的结论基本一致。

黄河口附近SPM年际变化发现不同年份中同一月份SPM分布情况总体趋同, 分布情况: 5月、7月SPM值较低; 11月SPM值较大。黄河口表层悬浮物浓度的季节性差异显著, 风场等气候因素应是该季节性差异的主要原因^[28-29]。在黄河口附近选定A、B、C三个区域, 分别统计其SPM均值, A、B区SPM值相对于C区大些, 在3个区域中位于B区悬浮物浓度为相对较高, 可以看出黄河每年携带大量泥沙入海对渤海悬浮物浓度B区的贡献很大, 可能是A区域与B区域离黄河口比较相近, C区域离黄河口较远, 也有可能离悬沙物源较远、以及水深有关, 这与吴昊^[30]等的黄、渤海悬浮物浓度近岸浓度高, 远岸浓度低结果一致。文本用2013—2015年输沙量数据和海面风场数据分析对B区域悬浮物浓度的影响(如图 7所示)。

黄河口悬浮物浓度变化差异主要由风速和输沙量引起的。本文利用6 h时间分辨率CCMP海面风场数据和黄河利津水文站的月输沙量数据进行分析。图 8a为2013—2015年悬浮物浓度与当月至前3月输沙量间的相关系数图, 图 8b为2013—2015年悬浮物浓度与当日至前30 d风速间的相关系数图。通过分析数据所对应时间上的输沙量和风速变化与悬浮物浓度的关系, 能更好地揭示水域不同时期与哪个要素的关系更为密切。从输沙量与悬浮物浓度关系来看(图 8a), 与前3月输沙量相关系数更显著, 其他月份悬浮物浓度与输沙量变化的相关性较小, 从图 8c中可知, 与前1月的输沙量间相关系数为–0.42。从风速与悬浮物浓度关系来看(图 8b), 与前30 d和26 d风速相关系数更大, 从图 8d中可知, 与前30 d相关系数为0.59。对照当月、前1月、前2月和前3月输沙量变化与悬浮物浓度变化的关系, 来揭示输沙量对悬浮物韵律的控制。从变化关系来看(图 9a), 悬浮物浓度与前1月输沙量变化基本一致, 说明悬浮物浓度对前1月输沙量存在滞后性。从与风速变化关系来看(图 8b), 与前30 d风速的相关性明显高于与当日风速变化的相关性, 悬浮物浓度与前30 d风速变化基本一致, 说明悬浮物浓度对风速变化具有明显的滞后性(如图 9b), 这与已有学者研究结果一致^[29]。悬浮物浓度与风速和输沙量的滞后性基本一致。这跟崔廷伟等^[28]的研究, 风速具有滞后性, 大风过程可以显著改变全渤海的悬浮泥沙空间分布格局吻合。

4 结束语

本文基于Vanhellemont和Ruddick的算法, 对2013—2015年黄河口Landsat-8 OLI数据进行大气校正, 利用2015年6月黄河口高光谱实测反射率数据检验了OLI数据的大气校正结果, 经检验发现, 两种反射率之间较一致, 相关系数为0.95, 在蓝、绿、红波段处平均相对误差分别为16.4%、17.3%和25.7%。误差的主要原因是实测时间与卫星过境时间不匹配。利用MODIS数据检验了OLI数据大气校正结果。发现Landsat-8 OLI数据蓝、绿、红波段与MODIS数据蓝、绿、红波段的相关系数分别为0.89、0.91、0.90, 在蓝、绿、红波段处平均相对误差为26.8%、17.8%和18.7%。OLI数据和MODIS数据计算得到的SPM结果进行对比发现OLI数据SPM与MODIS数据SPM之间的一致性较好, 相关系数为0.85, 平均相对误差为20.5%。误差的主要原因是OLI与MODIS的传感器的不同。该算法在黄河口水域具有一定的适用性。本文利用时间序列OLI数据分析黄河口悬浮物浓度的时空分布特征, 黄河口悬浮物浓度具有明显的季节变化特征, 主要受风速和输沙量的影响黄河口悬浮物浓度变化存在差异, 风速和输沙量对悬浮物浓度的影响存在一定的时间滞后现象。