1 引言(Introduction)

海洋水色遥感技术在获取海洋生态环境数据方面具有监测范围广和时效性强等特点，很好地弥补了传统监测方法获得的数据在时空上分散的缺点，在海洋生态环境研究中越来越受到重视.

卫星传感器所接收到的总辐射亮度中，不仅有包含海水物质光学性质的离水辐亮度，还有许多其他成分，而这些成分所占比例高达80%~90%(刘良明，2005).大气校正的目的就是要从复杂的信号中消除太阳和卫星几何条件、大气衰减作用等影响，从而提取出离水辐亮度.因此，大气校正是利用遥感数据反演物质浓度极其关键的一个步骤.海水通常被划分为一类水体和二类水体两个类型(Morel et al.，1999).对于一类水体(多为开阔大洋)，浮游植物对水体的光学特性起着决定性作用，目前已有了较为成熟的清洁水体大气校正算法，也称近红外波段(NIR)大气校正算法.对于二类水体(多为近岸海域或内陆湖泊)，其水体组分复杂，除浮游植物外的其他光学有关成分(如黄色物质和悬浮无机物等)也会影响水体的光学特性，尚无很好的通用算法.

在诸多的针对二类水体的大气校正模型中，精度较高的为有准确大气剖面辅助数据参与的辐射传输模型法.但这类大气校正模型的精度取决于模型所需的大气剖面辅助数据与实际大气状况的吻合程度，因此，在目前的观测条件下该模型难以适用于所有海域(陈军，2014).除此之外，近红外波段迭代法、邻近象元法和短波红外波段(SWIR)大气校正法也是目前常用的几种算法.近红外波段迭代法有两个假设：①红光和近红外波段的水体总吸收系数主要由纯水的吸收系数决定；②颗粒物后向散射在红光和近红外波段内随光谱是线性变化的(Arnone et al.，1998).对近岸浑浊度较高的水体(如渤海)，以上两个假设可能不成立.邻近象元法假设气溶胶类型在较小的空间范围内(约50~100 km)保持不变，因此，可参照或直接利用邻近清洁水体的气溶胶类型来确定浑浊区域的气溶胶类型(Hu et al.，2000).对于渤海这样的较浑浊水体，其毗邻的黄海也大多为二类水体，很难从其附近找到清洁水体，此方法也不适用.SWIR方法可视为适用于一类水体的NIR方法在二类水体中的版本.在二类水体中，近红外波段的清洁水体的假设不再成立，而短波红外波段为清洁水体假设依旧成立，基于此发展出了以短波红外波段为基准波段的SWIR大气校正算法(Wang，2007).目前，已经有许多学者将此方法应用于中国沿海区域(Wang et al.，2007；Wang et al.，2009；陈军，2014).

现有的SWIR方法假设在对数坐标下，大气校正因子在外推波段与基准波段间按线性关系变化.当波段跨度不大时，这一线性近似与理论结果相符，而当波段跨度较大时，则存在一定差异.针对SWIR方法存在的不足，本文提出一种改进方法，并选取有代表性的渤海海域和应用广泛的MODIS遥感数据对改进方法的效果进行定量评估.

2 SWIR大气校正方法的基本原理(Basic principle of SWIR atmospheric correction algorithm)

卫星辐射计接收到的总辐射亮度包括在大气传输过程中产生的辐亮度、离水辐亮度、海浪白冠的反射、海底反射的辐亮度及太阳耀斑(Gordon et al.，1994；潘德炉等，1997).其中，离水辐亮度是经过大气校正之后希望得到的最终结果.太阳耀斑一般可以通过调整辐射计的轨道设计和探测方向来避免，在较深的一类水体或浑浊的二类水体中海底反射可以忽略.因此，总辐射亮度为大气传输过程中产生的辐亮度、离水辐亮度和海浪白冠的反射这3项之和，即：

(1)

式中，i表示传感器的第i个通道;λ_i表示第i个通道对应的波长(nm)；L_path(λ_i)为海表面辐射在大气传输过程中产生的辐亮度(W·m^-2·μm^-1·sr^-1)；L_w(λ_i)为离水辐亮度(W·m^-2·μm^-1·sr^-1)；L_wc(λ_i)为海浪白冠的反射(W·m^-2·μm^-1·sr^-1)；t(λ_i)为水体到传感器的大气漫透射率.

海浪白冠的反射主要受风所引起的海面波的斜率影响，为风速的函数，有较精确的经验公式可用.大气传输过程中的辐亮度包括大气瑞利散射、大气气溶胶散射(包括气溶胶和大气分子的交互项).瑞利散射有精确方法进行计算，而气溶胶散射由于其复杂性，其计算成为大气校正中的难点.

清洁水体中，近红外波段满足黑体假设，即离水辐亮度近似为0(Gordon et al.，1992)，因此，一类水体处理气溶胶散射的方法是选取近红外波段为基准波段，由此外推出其他波段的气溶胶散射，即清洁水体大气校正法.该方法建立的外推关系式如下：

(2)

(3)

式中，λ_j为基准波段，λ_i为外推波段；ε(λ_i，λ_j)为i波段的大气校正因子；c为常量，用于外推到其他波段大气校正因子；F′₀(λ_i)、F′₀(λ_j)为从大气圈外入射的太阳光中波长为λ_i、λ_j的光通量F₀(λ_i)、F₀(λ_j)两次穿过大气层被衰减后的值；L_a(λ_i)、L_a(λ_j)分别为波长为λ_i、λ_j的气溶胶散射产生的辐亮度(W·m^-2·μm^-1·sr^-1).

针对MODIS，选取第15、16波段(即748和869 nm)为基准波段，即令i=15，j=16，可计算气溶胶散射L_a(748)和L_a(869)；然后利用式(2)和式(3)计算出常数c；进一步固定j=16，利用式(3)可计算出不同波段的大气校正因子ε(λ_i，λ_j)；再代入式(2)可得相应波段的气溶胶散射项；最后利用式(1)即可求得该波段的离水辐亮度.

在浑浊的二类水体中，近红外波段的离水辐亮度不再为0，不宜选取近红外波段为基准波段，但短波红外波段(MODIS的第5~7波段，即1240、1640、2130 nm)的离水辐亮度仍然为0.MODIS传感器的1640 nm波段有4个阵元无法正常工作，使得1640和2130 nm波段组合不能用于SWIR大气校正.鉴于此，对于原始的MODIS短波红外数据而言，仅有1240和2130 nm波段组合可以用于SWIR大气校正.Shi等(2009)的研究表明，1240 nm波段的黑体假设在中等及以下浑浊水体中是有效的，但在非常浑浊的水体中，1240 nm波段水体散射的贡献不能被忽略.然而，没有一种大气校正方法被证实对任何水质类型和区域都有效(Chen et al.，2013).因此，基于1240和2130 nm波段组合的SWIR大气校正方法仍可能是一种合理的选择.

鉴于选取1240或是1640 nm波段对于本文关注的问题并无太大影响，因此，采用第5、7波段(即1240和2130 nm)为基准波段进行分析.Wang(2007)借鉴清洁水体大气校正法，针对二类水体，选取短波红外波段为基准波段建立了SWIR大气校正算法.即在上述计算步骤中，先令i=5，j=7，计算气溶胶散射L_a(1240)和L_a(2130)；然后利用式(2)和式(3)计算出常数c；进一步固定j=7，利用式(3)外推出不同波段的大气校正因子ε(λ_i，λ_j)；再代入式(2)中计算相应波段的气溶胶散射项；最后代回式(1)即可求得该波段离水辐亮度.此方法即为SWIR大气校正法.

3 SWIR大气校正方法的改进(Improved SWIR atmospheric correction algorithm)

在Wang(2007)提出的SWIR大气校正算法中有一个重要的假设，即：在对数坐标下，大气校正因子在外推波段与基准波段间按线性关系变化.应该说，这一假设在整体上与11个气溶胶模型通过辐射传输模型所计算的模拟值是相符的(图 1).但仔细分析图 1、图 2可以看出，随着波长的减小，曲线呈现上扬趋势，并非严格的(对数)线性函数.

对MODIS数据的分析表明，渤海海域的大气校正因子大于1.因此，针对渤海海域，选取图 1中T50、T90、T99、C50、C70、M50、M70、C90、M90种气溶胶类型为观察对象，分析图 1中不同气溶胶类型的曲线.分段计算两个相邻波段之间的曲线斜率，取所有气溶胶类型曲线斜率的平均值和斜率相对变化率的绝对值，结果如表 1所示.其中，斜率S满足式(4)，斜率相对变化率△S的绝对值满足式(5).

(4)

(5)

式中，λ_i、λ_j为两相邻波段，且λ_i>λ_j.

从表 1可以看出，当λ_i＜869 nm时，斜率相对变化率较λ_i>869 nm时明显增大(即从4.382%增大到19.627%)，即从λ_i=869 nm开始，曲线有明显的转折.进一步分析，对λ_i>869 nm部分的曲线进行线性拟合，并与图 1曲线进行对照，从图 2可明显看出，曲线在近红外波短波长(869 nm)附近产生了转折.图 1和图 2中λ₀是基准波段，θ₀为太阳天顶角，θ为卫星天顶角，Δ$\phi $为相对方位角.

对这一差异给离水辐亮度计算结果造成的影响进行分析，计算离水辐亮度的目的是为了得到遥感反射率，从而计算出叶绿素浓度.通常采用不同波段的遥感反射率数据组合进行叶绿素浓度反演(Gurlin et al.，2011；Chen et al.，2012；许士国等，2015)，因此，不同波长的遥感反射率数据之间的差异会影响到后续工作的精度.针对这一问题，一个自然的想法是，做分段处理可以更好地逼近理论曲线，改进计算结果，而又不至于带来过大的额外的计算量.因此，作为初步的尝试，本文采用了两段线性近似，将869 nm选为分界点，即假设：

(6)

(7)

式中，c₁为常数；c₂(λ_i)为随波长变化的变量.其余的做法均与Wang(2007)提出的SWIR校正算法一致.

在实际计算时，若外推的目标波段的波长λ_i>869 nm，则与SWIR方法一致(仅出现一个常数c₁)；若λ_i＜869 nm，则首先外推至λ_i=869 nm(出现常数c₁)，然后再确定ε(λ₈₆₉，λ₂₁₃₀)在图中哪两条曲线之间，以及到其距离的比例关系，将此比例关系应用到外推的下一波段(图 1中是从右向左)，从而计算出c₂(λ_i)(λ_i＜869 nm).

4 改进算法的效果分析及讨论(Analysis and discussion of improved algorithm results)

选取2010年10月4日和2012年4月4日过境渤海(37°~41°N，117.5°~122.5°E)空间分辨率为1 km的MODIS图像作为数据源.以1240和2130 nm波段为基准波段，分别用SWIR算法和改进的SWIR算法计算出了443、488、531、551、667、678 nm波段的离水辐亮度.虽然渤海大多海域是二类水体，但为了保证精度有必要对其水体类型进一步分析，在此引入浑浊度T_ind来判断渤海的水体类型(Shi et al.，2007)：

(8)

通过计算，渤海海峡以东海域浑浊度均小于1.3(大体在1.1~1.3之间)，而渤海海峡以西海域浑浊度均大于1.3.分别用NIR法和SWIR法进行大气校正，可以发现浑浊度在1.1~1.3之间的海域，NIR和SWIR法计算的离水辐亮度只相差不到5%，而在浑浊度大于1.3时二者差异平均值高达40%~50%，这与Wang等(2007)的结论基本一致.由于SWIR波段的信噪比低，在结果相差很小的情况下优先选取NIR法.因此，本文选择T_ind=1.3为阈值，当T_ind>1.3时，采用SWIR方法；当T_ind＜1.3时，采用NIR方法.在此主要针对渤海海峡及其西部海域(即T_ind>1.3)进行分析，来讨论算法的效果.

图 3a和图 3b分别给出了两景图片443、531、667 nm波段的离水辐亮度的分布，为表述方便，图中及后文用SWIR2表示改进的分两段线性拟合的SWIR方法.从图中可以看出，用SWIR方法和SWIR2方法得出的离水辐亮度的分布趋势大致相同：离水辐亮度在近岸区域相对较大、渤海中部较小、渤海海峡处最小的总体特性与已有的研究结果也是一致的(王峰等，2009；Qiao et al.，2010；牟冰，2014).488、551、678 nm波段的结果与此类似，不再赘述.

图 4a和4b比较了两景图片分别用SWIR方法和SWIR2大气校正法计算得出的离水辐亮度.总体而言，在渤海海域SWIR2方法所得结果比SWIR方法偏小(大多数散点位于1∶1线下方)，两者的差异随着外推目标波段与基准波段波长差别的增大而增大(图 1中，基准波段(2130 nm波段)在图的右侧，而外推目标波段按678、667、551、531、488、443 nm的顺序逐渐向左侧移动).从图 1可看出，大部分曲线在近红外波短波长附近(869 nm)明显向上转折，这正是导致出现上述偏差规律的原因.

图 5为两种方法得出的6个波段的离水辐亮度的相对差异在渤海的分布(图片中空白区域是MODIS遥感影像缺失数据的区域，主要是667和678 nm两个波段有数据缺失).可以看出，随着波长的减小，改进方法与原方法的差异是逐渐增大的，这对后续利用波段遥感数据组合反演叶绿素浓度的精度可能有显著的影响.从图 5可以看出，二者的相对差异总体在15%以下.两景图片计算结果的相对差异分布不同，这很可能与不同时期由于各种因素导致的渤海海水组分不同有关.

5 结论(Conclusions)

本文在Wang提出的SWIR大气校正法的基础上，通过分段线性拟合的方法对大气校正因子ε(λ_i，λ_j)表达式中的c进行了改进.基于渤海海域的MODIS图像对比分析了改进后的方法与原SWIR方法所得离水辐亮度的差异，结果表明，改进后方法的结果比原方法整体偏小，在整个渤海的相对差异大多处于1%~15%之间.对于不同时期不同海域的遥感数据，改进方法获得结果的改善程度不同；不同波段的数据改进后的差异也不同(整体趋势是差异随着波长的减小而增大)，而后期计算叶绿素浓度是需要不同波段的遥感反射率数据进行组合，不同波段数据的差异将间接影响到反演出的叶绿素浓度的精度.当然，这只是一个初步的工作，改进方法的最终效果还需根据具体水色参数的反演结果(如叶绿素浓度)和实测数据的对比进行评价，但在卫星遥感数据的反演中关注此问题仍是必要的.

致谢: 本文所用的MODIS卫星遥感数据由美国航空航天局(NASA)水色处理中心提供(http://oceancolor.gsfc.nasa.gov/)，在此表示感谢！