一种浑浊水域气溶胶光学厚度与浑浊度的反演算法

技术领域

本发明属于卫星遥感技术领域，具体涉及浑浊水域气溶胶光学厚度与浑浊度的反演算法。

背景技术

气溶胶光学厚度(AOD)是气溶胶最重要的参数之一，代表着介质的消光系数在垂直方向上的积分，它可以表征大气的环境质量和混浊程度。当前AOD的监测主要有地基观测和卫星遥感两种手段。地基观测方法主要使用太阳光度计进行测量。然而地面监测往往只能得到站点周围的信息，没有大范围的整体观测数据，并且地面站点的数量及分布也受到限制。但是气溶胶卫星遥感反演就不存在这些问题，并且空间分辨率高、探测范围广、精度更高、获取数据快(Wei et al.,2019)。近年来，随着遥感卫星技术的发展，对于大气气溶胶的反演研究不断深入；定量反演气溶胶光学厚度也取得了很大的进展，同时气溶胶的精确反演也对于研究气候变化与空气污染防治做出了重大的贡献。

现有气溶胶卫星反演算法都是针对不同地表类型和气溶胶类型组成的不同，以不同的原理来对气溶胶进行反演。目前，AOD反演算法主要有：暗像元法、深蓝算法、SARA算法、结构函数法、多角度遥感法、双星协同法、偏振特性遥感等。其中暗目标算法又称为浓密植被算法，其核心原理是植被覆盖地区的红光波段的地表反射率和蓝光波段的地表反射率与短波近红外波段的地表反射率存在一定的线性关系(Remer et al.,2008)。深蓝算法主要选取蓝光波段作为工作波段(Hsu et al.,2006)。深蓝算法是利用在蓝波段大气反射较强，地表反射较弱的特点，并假设同期的地表反射率不变，将晴天的地表反射率代入公式反演气溶胶。SARA算法在设定地表为朗伯体、近似单次散射和气溶胶单次散射率(SingleScattering Albedo，SSA)及不对称因子(Asymmetry Factor，AF)不存在空间变化前提条件下，引入大气后向散射率、气溶胶散射相函数的近似表达函数(Bilal,Nichol,&Chan,2014)。结构函数方法反演AOD的实质就是先找到该地区的地表反射率的分布规律。假设分布规律保持不变，则表观反射率的空间变化认为是大气气溶胶光学厚度产生的贡献，从而反演出气溶胶光学厚度(Tanre et al.,1988)。

为了方便对水色遥感的研究，根据水体中悬浮和溶解物对水体光学属性影响不同，Morel和Prieur(1977年)把海洋水体分为一类和二类水体，一类水体主要指光学属性变化基本来自于浮游植物影响的水体，通常呈深蓝色，例如大洋。而二类水体不仅受浮游植物和相关的颗粒影响，还受到悬浮颗粒和沉积物的影响，对于一些近海浑浊水域，目前利用卫星遥感技术反演AOD卫星是一个很大的挑战。浑浊水域主要指近岸、河口等受陆源物质排放影响较为严重的地方。作为最主要的气溶胶卫星遥感数据源之一的MODIS产品反演陆地与海洋Ⅰ类水体AOD算法较多而且比较成熟，但是现有的算法都无法反演Ⅱ类水体的气溶胶参数。

在一些海浑浊水域反演AOD十分困难，是因为气溶胶浓度关系，且气溶胶的组成成分也十分复杂且多样，包括沙漠尘埃、海盐物质等，这为AOD的精准遥感反演增加了难度，并且一些沿海由于大量泥沙的堆积导致水体比较复杂。MODIS反演海洋上的AOD目前采用的是MODIS DT-ocean算法，DT-ocean算法(Levy et al.,2013)使用了不同精细气溶胶粒子和粗气溶胶粒子的加权组合在6个波段(0.55、0.65、0.86、1.24、1.63和2.11μm)的光谱反射率，通过建立查找表反演出AOD与FMF精细模式分数。此外，查找表假设除0.55微米(使用0.005的固定离水反射率)之外，所有波段的water-leaving radiance为零。而传统大气校正中采用的近红外波段(0.86μm)water-leaving radiance近似为零的假设却只适用于公海，对于混浊的沿海水域(二类水体)就不再适用。因为在浑浊水域，浅水海底的反射(尤其是在0.66μm波段)和水中悬浮或溶解的颗粒物质的反射(尤其是在0.55、0.86μm波段)会贡献water-leaving radiance。因此，在DT-ocean算法中浑浊水域像素将会被标记，并且不被反演。

目前，一些沿海的AOD卫星产品的精度和适用性还很低。NASA的MODIS气溶胶产品在Ⅱ类水体上空也通常为空白值,或者与实际测量值存在较大的误差而无法使用。Wang(2017年)分析了在全球和区域范围内无云条件下的MODIS AOD反演不可用性，以揭示AOD不能被反演仅仅是因为水的混浊度(而不是其他因素，如云层)，通过分析发现在所有公海上，这种数据可用性几乎是完全的，并且朝着海岸线急剧下降(下降90-100％)，在全球平均水平上，AOD在沿海水域的不可用性约为20％。气溶胶是研究这些区域大气环境质量的重要因子之一和遥感定标的重要参数，并且约60％的人口生活在沿海地区，扩大这些地区上卫星遥感数据集至关重要，因此对Ⅱ类水体开展气溶胶性质反演及分布的研究具有重要意义。

沿海水域的特征通常是由于海床再次悬浮或河流排放的颗粒而产生的高浓度悬浮有机物和无机物，悬浮物浓度直接影响水体的浑浊度和颜色。长江每年平均向河口排放约9×1011m3的水和4×108吨的泥沙,传送的泥沙量主要受枯水期和汛期的影响，汛期(6月至9月)时输送的泥沙量约占全年输送泥沙总量的87％(chenet al.2003年)。夏季，长江汇入东海扩散之后形成了浑浊海域，向东延伸了104-105km2的区域(Zhang et al.2007)。长江约50％的泥沙堆积在河口处(Shen et al.2001)，SSC值跨越两个数量级，范围从20mg/l到2500mg/l甚至更多。

对于沿海悬浮物浓度的卫星遥感分析目前也已经有了一些进展。如Miller和McKee(2004)分析了大气校正后的MODIS第1波段反射率与总悬浮物(TSM)浓度之间的回归关系，发现MODIS Terra 250m波段1(620～670nm)数据与实地测量的TSM数据之间可以建立起良好的线性关系。Doxaranet al.(2002,2009)一种从卫星数据的可见光和近红外(NIR)波长来确定悬浮物浓度的方法，通过这些大量的现场测量，在可见光到近红外波段将遥感反射率和悬浮物浓度之间建立了经验关系，发现近红外波段XS3(790–890nm)的遥感反射率Rrs和可见光波段XS1(500–590nm)与XS2(790–890nm)遥感反射率的比(Rrs(XS3)/Rrs(XS1)和Rrs(XS3)/Rrs(XS2))与SPM之间有很好的相关性。这些算法可以适用于高悬浮物浓度的区域，但是用于低悬浮物浓度区域时会产生较大的误差。

针对以上问题，本发明提出了一个改进的浑浊水域气溶胶光学厚度反演算法，主要利用在浑浊水域近红外波段(0.86μm)的water-leaving radiance不再近似为0，而短波近红外波段(2.1μm)处的water-leaving radiance可近似为0这一特点，选取了2.1μm来反演气溶胶光学厚度，并且不采用传统的近红外波段大气校正的方法来研究泥沙量，而是通过海洋BRDF模式来定义浑浊度，并且反演出0.86μm的浑浊度，从而建立浑浊度与泥沙量之间的关系。

发明内容

本发明的目的在于提供一种能够正确反映大气的环境质量的浑浊水域气溶胶光学厚度与浑浊度的反演算法。

本发明提供的浑浊水域气溶胶光学厚度与浑浊度的反演算法，内容主要包括四部分：

一、通过四种气溶胶粒子(沙尘型，水溶性，海洋性，煤烟型)的不断迭代来确定研究区域的气溶胶类型，利用MODIS数据的2.1μm波段，通过6S辐射传输模型

二、使用海洋BRDF模型

三、再次利用6S辐射传输模式来反演出0.86μm处的浑浊度；

四、结合站点的实测数据验证该算法反演AOD的适用性，并进一步结合泥沙量数据，与反演得到的地表浑浊度进行回归分析，研究他们之间的响应关系。

具体步骤为：

(1)确定研究区域的气溶胶类型；将气溶胶类型看做沙尘型、水溶性、海洋性和煤烟型四种气溶胶粒子的组合；假设四种气溶胶粒子的体积浓度百分比范围后，让四种粒子在各自区间进行迭代，利用6S辐射传输模型计算出理论表观反射率，将研究区域每个格点的理论表观反射率与观测表观反射率最匹配时所对应的四种粒子百分比的组合记录下来，所有格点记录的四种粒子百分比的平均值就代表该地区的气溶胶类型；

(2)选择2.1μm通道进行反演；因为2.1μm处的离水辐射率可近似为0；利用6S辐射传输模型，输入相关参数计算出表观反射率理论值，得到AOD、地表反射率、太阳天顶角、卫星天顶角、相对方位角变化的查找表；根据每个格点上的地表反射率(由MOD09GA得到，这里，MOD09GA提供500m分辨率的MODIS波段1-7的每日表面反射率，并提供1km的观测和地理位置统计数据)、太阳天顶角、卫星天顶角、相对方位角(由MOD02得到，这里MOD02提供1km分辨率的MODIS各波段的表观反射率、地理坐标信息以及几何参数)，对查找表进行插值，找到表观反射率观测值(由MOD02得到)与表观反射率理论值最匹配时所对应的AOD；

(3)使用海洋BRDF模式，定义浑浊度；双向反射率表示遥感器观测到的地球表面反射率不仅取决于目标本身，还取决于太阳和传感器相对于目标的位置；双向反射分布函数(BRDF)模型是双向反射率的数学表达式，用于计算在不同的太阳、传感器几何参数下观察到的反射率；将海洋BRDF模型做出改进，定义一个新的参数：浑浊度；

(4)使用0.86μm反演出浑浊度；依旧利用6S辐射传输模型生成浑浊度、AOD、太阳天顶角、卫星天顶角、相对方位角变化的查找表；根据每个格点上的AOD(由前一步的反演得到)、太阳天顶角、卫星天顶角、相对方位角(由MOD02得到)，对查找表进行插值，找到表观反射率观测值(由MOD02得到)与表观反射率理论值最匹配时所对应的浑浊度；

(5)AOD反演结果的精度检验：对地面站点AOD资料进行时空匹配；将气溶胶反演结果与地基数据观测的AOD进行线性回归分析，分析反演结果与实测值的相关性；

(6)建立浑浊度与泥沙量之间的关系：对泥沙量资料进行时空匹配；通过经验模型分析法、理论模型分析法、半分析模型分析法等多种方法，找到最相匹配的浑浊度与泥沙量之间的经验公式；将经验公式用于其他个例以检验反演浑浊度的精度。

本发明的反演算法中，涉及有关算法原理有：

1、AOD反演原理

当通过卫星传感器进行测量时，大气和地表的贡献取决于大气和地面光学参数以及相对方位角、太阳天顶角以及卫星天顶角；卫星传感器所接收到的地面目标的总辐射亮度，是经大气衰减后的地表面辐射亮度和大气本身的路径辐射之和(Tanré et al.,1996)；大气程辐射是指太阳辐射在大气传输过程中各组分及气溶胶微粒散射后直接到达传感器的辐射。

对于无云且大气水平均一的地-气系统，朗伯地表(各向同性反射)上空卫星观测的TOA方程可以表达为(Kaufman et al.,1988)：

式中，θ

其中，ρ

根据方程(3)可知，表观反射率不仅是气溶胶光学厚度的函数，也是地表反射率的函数。其中，角度数据与表观反射率可由卫星数据得到，而单次散射反照率、散射相函数等参数可通过符合实际情况的大气模式及气溶胶模型确定。实际操作中，我们通过辐射传输模型输入相应参数来获得这些量。因此只要求得地表反射率和气溶胶模型，理论上便可以通过公式(3)反演气溶胶光学厚度。

2、海洋BRDF模型

通过考虑白帽，太阳耀斑和浑浊度的影响来计算海洋表面的BRDF；假定海洋表面的反射率ρ

ρ

其中，ρ

而R

辐照度反射率R

Morel将总后向散射系数分为2个部分：

其中b

其中b是色素的散射系数，计算公式为：

b＝0.3C

总吸收系数写为：

a(λ)＝u(λ)·K

其中u(λ)是与波长有关的函数，计算如下：

K

K

根据Morel的模型，反射率R

3、浑浊度的反演以及与泥沙量建立关系

使用0.86μm反演出浑浊度，因为在干净的水域，0.86μm的water-leavingradiance可近似为0，然而在浑浊水域2.1μm波段的water-leaving radiance可近似为0，但是0.86μm的water-leaving radiance却不能忽略，这其中就包含了水域浑浊度的信息。依旧利用6s辐射传输模型生成浑浊度、AOD、太阳天顶角、卫星天顶角、相对方位角变化的查找表；根据每个格点上的AOD(由前一步的反演得到)、太阳天顶角、卫星天顶角、相对方位角(由MOD02得到)，对查找表进行插值，找到表观反射率观测值(由MOD02得到)与表观反射率理论值最匹配时所对应的浑浊度。接下来对泥沙量资料进行时空匹配，通过经验模型分析法、理论模型分析法、半分析模型分析法等多种方法找到最相匹配的浑浊度与泥沙量之间的经验公式，最后将经验公式用于其他个例以检验反演浑浊度的精度。

附图说明

图1为2020年8月15日MYD04_L2 AOD分布图。

图2为算法流程图。

图3为绿色区域为本文的研究区域，红点为地面观测站点。

图4为SONET观测站2020年8月16日上海站点不同气溶胶组分百分比图。其中，BC：黑碳气溶胶；BrC：有机碳气溶胶；CM：粗粒子气溶胶(矿物沙尘或海盐)；FS：精细模式散射成分(如硫酸盐、硝酸盐和轻质非吸收有机物)；AW：水溶型气溶胶。

图5自定义气溶胶类型下表观反射率理论值与表观反射率观测值的线性回归。

图6为2020年8月15日上海沿海区域的2.1μm通道反射率分布图。其中，(a)为表观反射率分布图，(b)为地表反射率分布图。

图7为2.1μm表观反射率(a)与地表反射率(b)的统计图。

图8为表观反射率对于AOD的敏感性实验。其中，(a)是0.65μm通道，(b)是0.86μm通道。

图9为表观反射率对于AOD的敏感性实验(2.1μm通道)。

图10为表观反射率对于地表反射率的敏感性实验。其中，(a)是2.1μm通道，(b)是0.65μm通道。

图11为MYD04_L2 AOD产品以及反演得到的AOD(红色边框区域为反演区域)。

具体实施方式

下面结合附图和实施例对本算法进一步说明。

实施例1

国内外研究者通过不同的卫星数据来实现对AOD的反演已经取得了一定的成果，但是目前的算法还存在一定的适用性的问题，例如现在常见的海岸气溶胶反演算法只能适用于一类水体，而对于一些二类水体而言，这些算法还需要进行改进。在图1中可以看到MODIS的AOD产品在上海沿海区域就是缺测值。具体技术路线如图2所示。

研究区域和时间的选择：

地点选择上海沿海区域(经度120.9-122.5，纬度30.4-32.3)，也就是长江口到杭州湾这一区域时间选择2020年7,8,9月的晴天。

辐射传输模式的选择：

选择了6S大气辐射传输模型，6S大气辐射传输模型是模拟太阳辐射在地气系统中的传输过程的模型。主要作用是模拟传感器在无云条件下，考虑了非朗伯体情况下，真实大气各组分对于太阳辐射的吸收和散射作用后，卫星传感器理论上应该接受到的辐亮度。改进后的6S模型采用了最新的散射计算方法，考虑了多种散射，包括：分子散射、多次散射、气溶胶散射，使得计算得到的太阳光谱散射值精度更高(Lee et al.,2004)。

卫星数据的选择与处理：

卫星数据选择了MOD021km与MOD09GA；辐射定标：将原始的数字信号转换为反射率和亮温等实际物理参量；几何校正：将数据准确定位到特定的地理坐标系；数据重采样：将不同分辨率下的MODIS波段重采样到相同分辨率；裁剪：将数据裁剪成我们需要的范围。

气溶胶模式的选择：

气溶胶类型可以看做沙尘型，水溶性，海洋性，煤烟型气溶胶粒子的组合。根据Fan等人在2015年对选取的几个中国地区aeronet站点进行了对气溶胶各种组分的百分比的统计。发现每个地区水溶性所占比分最高，均在40％以上，而沙尘性均在40％以下，海洋性和煤烟性最少，均低于10％。同时根据SONET上海站点对不同气溶胶粒子的观测可以发现，水吸收性气溶胶粒子的占比是最高的，达到了80％，而黑碳型气溶胶粒子时占比最少的。

因此可以随机假设四种气溶胶粒子的体积浓度百分比为：沙尘型x1，水溶性x2，海洋性x3，煤烟型x4。对于各组分的取值可以限定为0≤x1≤20，40≤x2≤80，0≤x4≤5，x3＝100-(x1+x2+x4)。将x1,x2,x3,x4在各自区间进行迭代，迭代步长均为1，计算MOD02得到的四个角度的平均值作为输入的角度参数。使用站点资料的AOD作为输入。每迭代一次就运行6S得到F(θs)·T(θv)总透射率、S大气半球反照率、ρa路径辐射。现在共有151×203个格点，对于每个格点，结合地表反射率、F(θs)·T(θv)总透射率、S大气半球反照率、ρa路径辐射计算出理论表观反射率ρ*，然后与观测表观反射率值ρ进行对比。计算出ε＝|ρ*-ρ|，找到ε最小值的情况，并将对应的x1,x2,x3,x4记录下来。所有格点记录的x1,x2,x3,x4的平均值就代表该地区的气溶胶类型。经过计算得到了表1中的自定义气溶胶模型：

表1自定义的气溶胶类型

为了验证自定义的气溶胶模型的精确性，于是将自定义气溶胶模型作为6S模式的输出参数，计算得到这个区域的表观反射率理论值，再和MODIS观测的表观反射率进行线性回归分析，如图5所示，发现相关系数达到了0.8698，表观反射率理论值和观测值拟合得较好，说明自定义的气溶胶模型很好地反映了这个地区的气溶胶类型。

敏感性实验：

在反演之前需要做一个敏感性实验，验证表观反射率是否对气溶胶光学厚度以及地表反射率敏感。

图6显示的是2.1μm表观反射率以及地表反射率分布图，可以发现在长江口以及上海几个港口的下游都呈现表观反射率以及地表反射率数值的大值区，但是总体来说，表观反射率以及地表反射率的数值都在0.15以下。

图7是2.1μm表观反射率和地表反射率的统计分布图，可以看到表观反射率主要集中在0.06到0.08之间，而地表反射率主要集中在0.065到0.08之间，比常用的0.65μm的地表反射率值小得多。于是利用6S辐射传输模型，模拟在0.65μm通道和2.1μm通道下，表观反射率对气溶胶光学厚度以及地表反射率变化的变化。

如图8所示，对于0.65μm波段，在比较清澈的水域时(地表反射率为0.05)，表观反射率随AOD是明显的上升趋势，敏感性较好，但是当水域比较浑浊时(地表反射率为0.2)，表观反射率先是随AOD下降后趋于平滑，敏感性不再显著。对于0.86μm波段，无论是在清澈的水域还是浑浊水域，敏感性都不够显著。如图9所示，对于2.1μm波段，无论是在清澈的水域还是浑浊的水域，表观反射率都随AOD有明显的下降的趋势，呈现出较好的敏感性。因此在浑浊水域，使用2.1μm波段来反演比使用0.65μm波段更为准确。从图10可以看出无论是在0.65μm波段还是在2.1μm波段表观反射率都对地表反射率表现出较好的敏感性，且都呈现出表观反射率随地表反射率增大而增大的趋势。

建立AOD查找表：

本算法选择2.1微米通道进行反演，因为浑浊水域2.1μm处的water-leavingradiance可近似为0。确定好气溶胶模型之后，利用6S辐射传输模型输入相关参数计算出表观反射率理论值，得到AOD、地表反射率、太阳天顶角、卫星天顶角、相对方位角变化的查找表。其中建立的查找表中各变量的范围为：AOD:0.1-0.5(步长为0.01)，地表反射率：0-0.2(步长为0.02)，太阳天顶角：22-25(步长为1)，卫星天顶角：2-20(步长为3)，相对方位角：211-226(步长为3)。

对MODIS表观反射率数据进行气体校正、云掩膜后，根据每个格点上的地表反射率(由MOD09GA得到)、太阳天顶角、卫星天顶角、相对方位角(由MOD02得到)，对查找表进行插值，找到表观反射率观测值(由MOD02得到)与表观反射率理论值最匹配时所对应的AOD。

图11是反演得到的AOD(红色边框区域)以及MODIS AOD产品分布图，趋势大致吻合，根据反演的AOD分布，可以发现在长江口以及上海几个港口的下游地区的AOD都呈现大值，达到了0.45以上，并且这两个地方的地表反射率也是大值区(水域较浑浊)，说明这两个地方受到了沿海工业以及船舶的污染影响。

AOD的对比验证：

对地面站点AOD资料进行时空匹配。将气溶胶反演结果与地基数据观测的AOD进行线性回归分析，分析反演结果与实测值的相关性。

使用BRDF模式定义浑浊度：

通过考虑白帽，太阳耀斑和浑浊度的影响来计算海洋表面的BRDF。假定海洋表面的反射率ρ

ρ

如果我们假设海洋是朗伯反射，则ρ

如公式(4)-(12)所示，R

建立浑浊度查找表：

使用0.86μm反演出浑浊度，因为在干净的水域，0.86μm的water-leavingradiance可近似为0，然而在浑浊水域2.1μm波段的water-leaving radiance可近似为0，但是0.86μm的water-leaving radiance却不能忽略，这其中就包含了水域浑浊度的信息。依旧利用6s辐射传输模型生成浑浊度、AOD、太阳天顶角、卫星天顶角、相对方位角变化的查找表；根据每个格点上的AOD(由前一步的反演得到)、太阳天顶角、卫星天顶角、相对方位角(由MOD02得到)，对查找表进行插值，找到表观反射率观测值(由MOD02得到)与表观反射率理论值最匹配时所对应的浑浊度。

同泥沙量数据建立经验关系：

对泥沙量资料进行时空匹配，通过经验模型分析法、理论模型分析法、半分析模型分析法等多种方法找到最相匹配的浑浊度与泥沙量之间的经验公式，最后将经验公式用于其他个例以检验反演浑浊度的精度。