基于GOCI 数据黄海气溶胶光学厚度反演方法研究

郑小慎，张亚男，刘雨华

（天津科技大学 海洋与环境学院，天津 300457）

大气气溶胶是大气中的微量成分，在大气中气溶胶的含量虽然不高，但对调控地球辐射收支平衡及全球气候模式有着十分重要的作用，气溶胶成为研究大气污染和气候变化的重要参数。海洋上空气溶胶光学性质的研究对于估计电磁波的辐射、分析全球气候变化和提高定量遥感精度有重要意义（刘玉光，2009；沈永平 等，2013）。Griggs（1975）和Fraser（1976）利用单通道原理建立的海洋气溶胶光学厚度（Aerosol Optical Depth，AOD）反演算法，对于海表反射率较大的区域反演效果并不理想。多通道遥感算法反演气溶胶主要是利用通道反射特征避免下垫面反射的影响来提高反演精度。Durkee 等（1986）用AVHRR 的630 nm 和830 nm波段来反演AOD，Remer 等（2005）验证了660 nm和870 nm MODIS 海上AOD 在全球范围内的反演有效性，结果证明MODIS 数据的反演效果很好。刘毅等（2008）利用MODIS 气溶胶业务产品验证渤海、黄海、东海和日本以南海域的反演结果，分析结果表明基于MODIS 数据的双通道算法和业务算法对黄海的反演结果并不理想。多角度气溶胶遥感的原理是基于大气和地表对大气顶层卫星信号贡献的比率随观测角度的不同而不同的特点，将地表与大气信号区分开来。Veefkind 等（1998；2000）利用双角度算法（ATSR-DV）算法反演了不同区域的气溶胶并得到了较好的反演结果。

针对黄海海域，气溶胶光学特性的反演主要受到以下三点因素的影响，（1）黄海海域上空不但受到海洋水汽蒸发的影响，也受到包括沿岸工业排放、化石燃料及西北内陆黄沙等陆源物质的影响。因此，其上空的气溶胶包括非吸收性气溶胶和吸收性气溶胶。（2）区域下垫面受到内陆河流影响较为明显，近岸海域为典型的浑浊水体（牟冰，2014；王正，2013），下垫面的反射影响反演精度。（3）由于气溶胶的空间分布不均匀，时间变化较大，低时空分辨率的遥感数据不能在实际观测应用中起到作用。因此，减少以上因素的影响有助于提高黄海上空气溶胶光学特性的反演精度。在此背景之下，本文提出了一种针对黄海气溶胶光学厚度的新反演算法。该算法基于高时空分辨率卫星GOCI数据，通过增加吸收性气溶胶模态，减小由于水体反射引起的误差，达到提高反演精度的目的。

1 数据来源

1.1 地基数据

AFRONFT（Aerosol Robotic Network） 是一个由NASA 和PHOTONS 共同建立的气溶胶地基遥感观测网络（郭阳洁，2014），该网络现已经覆盖了全球主要的区域，目前全球共有超过500 个站点。通常，将AFRONFT 所测得的气溶胶光学厚度作为真值，来检验遥感反演的气溶胶光学厚度的准确性，图1 为AFRONFT 全球站位图。

图1 AERONET 全球站位分布图

AFRONFT 气溶胶数据包括1 级、1.5 级和2级数据产品， 1 级产品未经筛选和最终的校正；1.5 级产品经过了云污染筛选，未经过最终校正；2级产品经过了校正、云污染筛选和质量控制。黄海AFRONFT 站位少，且部分站位数据的可用性差，实验选择有数据质量保证且位于GOCI 刈幅内的Anmyon 站位（36.54°N，126.33°F），能够验证反演算法对浑浊水体的适用性。下载网址为https://aeronet.gsfc.nasa.gov。

1.2 GOCI 数据

2010 年韩国成功发射地球静止轨道卫星COMS，这颗卫星上搭载着世界上第一个静止轨道GOCI 海洋水色传感器，每天能够提供从00:16 UTM到7:16 UTM 的分辨率为500 m 的8 幅卫星遥感图像。由于时空分辨率较高，GOCI 数据对于监测黄海的海洋环境变化十分有意义，本文选用的GOCI遥感数据包括GOCI L1B 数据以及GDPS 处理的L2A 级数据，GOCI 数据直接从韩国COMS 卫星官网上免费下载，网站为http://kosc.kiost.ac.kr/。

2 反演方法

2.1 GPRO 算法

GOCI 传感器的优势在于利用静止轨道卫星获取高时空分辨率的遥感影像，GPRO 算法为AOD反演的业务化方法，通过去除离水反射率的影响反演AOD。反演流程如图2 所示。

图2 GPRO 算法流程图

对于浑浊水体，GPRO 仅利用去除离水反射率的490 nm 反射率反演550 nm 处的AOD。与AFRONFT 数据比较，GPRO 反演算法相对于MODIS 业务算法的反演精度更高（Lee et al，2010）。

GPRO 反演算法主要有三点不足：第一，对浑浊水体采用30 d 清洁窗确定离水反射率存在一定的误差；第二，GPRO 算法中并没有具体给出气溶胶Type；第三，GOCI 业务化算法对云去除过程过于严格导致在浑浊水体上空反演会出现偏差。

2.2 GNFW 算法

针对GPRO 算法中对离水反射率的处理以及查找表中气溶胶模态设置中存在的缺陷，本文提出一种基于GOCI 数据的气溶胶光学特性反演算法（GNFW 算法）。该算法的主要流程如图3 所示。

该算法主要由三部分组成，第一部分为表观反射率与辐亮度之间的相互转换反演，为了去除由于太阳天顶角及观测点与太阳之间距离的变化而引起的误差，该算法首先将GOCI 的辐亮度数据转化为表观反射率。该算法利用去除离水反射率的反射数据来减少由于水体辐射信号引起的干扰。第二部分为离水反射率除噪声反演，针对离水反射率产品存在明显条带噪声的问题，该算法提出了插值法去条带噪声以提高离水反射率的准确性和可用性。第三部分为查找表的建立和搜索。该算法的输出结果为气溶胶光学厚度（AOD），细颗粒比例（Fraction of AOD contributed by fine dominated model，FMF） 及气溶胶模态（Type）。

图3 基于GOCI 数据的气溶胶光学特性反演流程图

2.2.1 表观反射率的反演

基于下载的GOCI L1B 数据，根据大气校正方程（张民伟，2009）将卫星接收的辐亮度转化为表观反射率，式（1） 为表观反射率的计算公式（Teillet et al，2001），

其中，ρ（λ）表示λ 波长处的表观反射率；d 表示日地距离（天文单位）；θ 表示太阳天顶角；ESUN（λ）表示λ 波长处平均太阳辐照度；L（λ）表示λ 波长处传感器接收到的辐亮度，单位为W·m-2·sr-1·μm-1。

上式中日地距离是遥感影像成像时间的函数，表达式如下：

其中，JD 表示遥感成像儒略日。

波段平均太阳辐照度是大气层外的太阳光谱辐照度与传感器各波段光谱响应函数的积分，表达式如下：其中，S（λ）表示λ 波长处的光谱响应函数，BMSI（λ）表示λ 波长处大气层外太阳光谱辐照度。

采用Wehrli85 太阳光谱曲线计算GOCI 接收的波段平均太阳辐照度（张楠，2013；张璐 等，2014）。Wehrli85 太阳光谱范围为0.2～10 μm，GOCI 的波谱响应区间主要集中在300～1 000 nm 范围内，因此，对300～1 000 nm 区间内的大气层外太阳光谱辐照度和光谱响应函数进行积分，最终获得GOCI 传感器接收的大气层外平均太阳光谱辐照度（Mean solar exoatmospheric irradiances，FSUN），如表1 所示。

表1 GOCI 接收的平均太阳辐照度（W·m-2·sr-1·μm-1）

2.2.2 离水反射率反演

大洋水体下垫面比较均匀，离水辐射在可见光和近红外波段受到水体的影响较小，针对大洋水体的气溶胶反演已经很成熟并且反演精度较高。但是浑浊水体受到的影响比较复杂，影响其上空气溶胶散射的主要因素是离水反射率，计算离水反射率成为气溶胶光学特性反演的重要过程。

2.2.2.1 NIR-SWIR 算法反演离水反射率

对于浑浊水体区域，NIR 离水反射率算法处理后会出现明显的误差（Lavender et al，2005；Ruddick et al，2000；Siegel et al，2000），利 用SWIR 在近岸浑浊水体反射率基本为零的特性来反演离水反射率更为准确。GOCI 没有SWIR 波段，用MODIS 数据来弥补GOCI 波段设置的缺陷，使用NIR-SWIR 算法提高反演离水反射率精度，间接提高AOD 的反演精度。利用浑浊水体系数对水体进行区分，若为浑浊水体采用SWIR 波段进行处理，若为非浑浊水体采用NIR 波段进行处理。

GNFW 算法中使用的是MODIS 离水反射率每天的数据产品，这相比于GPRO 算法中通过30 天影像处理，而获取离水反射率更高效且精度更高。MODIS 数据的处理需要与GOCI 数据在时间上保持对应，GOCI 中心扫描时间为0:16、1:16、2:16、3:16、4:16、5:16、6:16、7:16 UTC，MODIS Terra 通过黄海上空的时间为1:30～3:30 UTC 之间，因此，实验采用MODIS Terra 1:30～3:30 UTC 之间经过黄海上空的遥感数据进行处理分析。

2.2.2.2 离水反射率反演去噪声处理

MODIS NIR-SWIR 算法处理后的离水反射率表现出明显的条带状噪声。采用插值法对MODIS离水反射率影像上的条带噪声进行处理，利用式（4）确定最终的噪声像元位置。

其中，l 表示数据矩阵的行，*表示数据矩阵l 行的所有列元素，fl 表示最近非噪声行，T 为检测噪声的阈值。

当确定输入图像的噪声像元和非噪声像元位置之后，对原始数据按照5*5 窗口进行搜索，若目标像元为噪声像元，利用窗口内非噪声像元的均值代替噪声像元原有的值。若目标像元为非噪声像元仍然维持原值。

2.2.3 辐亮度反演

将处理后反射率转化为辐亮度，输出结果将为后续利用查找表确定最终反演结果做准备。

2.2.4 气溶胶光学厚度反演

通过以上操作可以得到GOCI 数据8 个波段对应的去掉离水反射率之后转换得到的辐亮度结果。

卫星遥感的辐亮度大气校正方程可表达为

其中，L（λ）表示卫星探测的辐亮度；LR（λ）表示由于大气中分子瑞利散射引起的辐亮度；LA（λ）表示由于气溶胶散射引起的辐亮度；LRA表示由于大气中空气分子与气溶胶之间多次散射引起的辐亮度；Lr（λ）表示由于太阳耀斑引起的辐亮度；T（λ，θ）表示大气的直接透射率；Lwc（λ）表示由于白帽引起的辐亮度；t（λ，θ）表示大气的漫透射率；Lw（λ）表示离水辐亮度；λ 表示传感器波长；θ 表示卫星天顶角。

对于GOCI 数据8 个可见光波段的单次散射来说，如果不考虑多次散射引起的辐亮度，消除了每个波段中瑞利散射，太阳耀斑，白帽和离水辐亮度，气溶胶散射引起的辐亮度简化为：

其中，ωA（λ）为气溶胶单次散射比，ωA（λ）=KAS（λ）/KAE（λ），KAS为气溶胶散射系数，KAE为气溶胶消光系数，取经验值0.8。τA（λ）表示气溶胶光学厚度，FS（λ）为衰减的太阳辐亮度，PA（λ，θ，θs）为气溶胶相函数，θs为太阳天顶角，θ 为观测角。

2.2.5 查找表搜索

2.2.5.1 数据输入

Rstar5b 是由 Nakajima 和Tanaka（Nakajima et al，1986；1988；Stamnes et al，1988）提出的用于模拟大气—陆地—海洋之间辐射传输过程的辐射传输模型。该模型曾用于ADFOS/OCTS 及ADFOS-Ⅱ/GLI 产品反演算法中。Rstar5b 中假定陆地或海洋之上的平行平面为若干均匀层，模拟波段范围为0.2～200 μm。

2.2.5.2 查找表的建立和搜索

Rstar5b 中包括水蒸气、冰粒子、类尘埃、水溶性粒子、海盐粒子、烟尘粒子、火山灰粒子、75 %H2SO4、黄沙在内的9 种基本的气溶胶颗粒。气溶胶基础模态由这9 种基本气溶胶颗粒按照不同比例进行内部混合，Rstar5b 中气溶胶基本模态组成比例如表2 所示。

表2 Rstar5b 气溶胶基本模态组成

由表2 可以看出，气溶胶模态在内部混合比例上为固定值。

实验中针对黄海上空气溶胶模态的分析，选用类尘埃、烟尘、海盐、城市、黄沙这五种气溶胶基本模态进行外部混合。气溶胶的外部混合方式按照表3 所示的方式，其中√表示对该基本模态进行混合，空白表示不进行处理。

如表所示，分别用代码3、4、8、9、11 表示类尘埃、烟尘、海盐、城市、黄沙气溶胶模态，外部混合方式为将代码依次连接起来，依次表示为Type=38，Type=39，Type=48，Type=49，Type=98，Type=811，Type=911。

表3 Rstar5b 气溶胶模态外部混合组成

模型中用FMF 表示混合模态（Type=AB）中A 的混合比例，也就说，FMF 越大，混合模态中A的比例越高。

因此，Rstar5b 模型的输入变量可根据以上说明建立，特别说明的是其中方位信息可根据图像角度文件获得，对于输入变量中的自变量AOD＇及FMF 的组成分别如表4 所示。

表4 AOD'和FMF 的输入数据

如表5 所示，Rstar5b 的输入文件个数取决于太阳天顶角、出射光线角、相对方位角、AOD＇以及FMF 的个数。由于以上参数进行组合后数据量极大，因此实验中对原Rstar5b 模型进行修改，使得模型可以循环运行并逐个输出对应文件。

查找表（Look up table，LUT） 建立时对于给定角度参数，每种LUT 包括11 种不同的FMF 参数，每种FMF 中包括10 种不同的AOD＇。查找表搜索流程如图4 所示。在确定像元在LUT 中的起始位置后，对每种LUT 提取相应的110 行进行逐个计算，确定最佳LUT。实验中设计的气溶胶模态有7 种，将上述过程循环7 次获取7 个相应的FMF、490 nm AOD 及490 nm 辐亮度值。将每种LUT 对应的490 nm 处辐亮度插值结果与观测辐亮度进行对比，将最小差值作为最佳的LUT，同时输出相应LUT 的FMF 值和AOD 值。因此，输出结果包括每个像元的FMF、AOD 及气溶胶模态。

图4 LUT 搜索流程图

3 反演结果验证

基于Anmyon 站位的AFRONFT 数据对GNFW算法、GPRO 算法的反演结果进行了验证，其中GNFW 为本文提出的基于GOCI 数据气溶胶光学厚度反演算法，GPRO 算法为GOCI 数据业务化AOD反演算法。基于GOCI 数据时间与质量，选择Anmyon 站位无云天气状况下，2014 年和2015 年夏季（6-8 月）、2015 年和2014 年冬季（12-2 月）1.5 级AFRONFT 数据进行验证。

在Anmyon 站位基于AFRONFT 数据分别对比GPRO 反演数据和GNFW 反演数据，结果如图5所示，（a）图为夏季，（b）图为冬季。

图5 Anmyon 站位AERONET 数据与GPRO 算法和GNEW 算法反演结果对比

如图5 所示，在夏季，GNFW 算法相较于GPRO 算法更接近AFRONFT 数据，并且GPRO 算法在AOD 大于0.25 时有明显的低估。在冬季，在整个监测区间内GNFW 算法比GPRO 算法更接近AFRONFT 数据，并且主要的气溶胶光学厚度主要集中在0.02～0.3 区间内。冬季当气溶胶光学厚度在0.02～0.15 区间时，GNFW 算法与GPRO 算法都存在高估的情况，但是GNFW 算法仍更接近AFRONFT 数据。当AOD 大于0.15 时，GNFW 算法与GPRO 算法都存在低估的情况，但是GNFW算法仍更接近AFRONFT 数据。总体而言，GNFW算法相比于GPRO 算法更接近AFRONFT 数据。

表5 GNEW 算法和GPRO 算法的相关统计数据

表5 列举了GNFW 算法和GPRO 算法在夏季和冬季反演结果的相关统计数据，由于验证GOCI和AFRONFT 的数据是随机选择的，因此结果可以说明GNFW 算法在冬季和夏季都表现出较GPRO算法更好的反演效果。综上所述，GNFW 算法在浑浊水体上空能较为准确地反演AOD，同时GNFW算法的反演精度相比于GPRO 算法有提高。

GNFW 算法误差产生的可能原因，主要包括以下三点，第一，辐射传输模型中大气及气溶胶颗粒参数设定不准确，这其中包含大气组分浓度、大气压力廓线、大气质量、温度廓线、粒子尺度分布等（王欣睿等，2016）；第二，由于水体环境变化使得离水反射率随时间出现变化。第三，去除离水反射率时对大气漫透射率的估计出现误差。

4 结论

黄海海域气溶胶光学特性反演精度受到陆源物质的影响显著，近岸水体上空的气溶胶模态为吸收性气溶胶和非吸收性气溶胶共存的状态。本文在GOCI 数据业务化气溶胶反演算法GPRO 算法基础上，提出了改进的GNFW 算法。针对黄海海域气溶胶光学特性的反演，该算法基于Rstar5b 气溶胶模态内部比例，选用类尘埃、烟尘、海盐、城市、黄沙这五种气溶胶基本模态进行外部混合形成黄海上空气溶胶模态类型，兼顾了黄海上空非吸收性气溶胶和吸收性气溶胶共存的状态。对于浑浊水体，利用MODIS 波段设置的优势，进行离水反射率反演，进行GOCI 数据中离水反射率的去除，对于气溶胶反演精度的提高有重要作用。GOCI 数据的高时空分辨率在气溶胶观测中起到了重要作用。该算法能较好地反演浑浊水体上空的气溶胶光学厚度，相比于GPRO 算法在反演精度上有所提高。