CN102636143B - 一种气溶胶光学厚度遥感反演方法 - Google Patents

Abstract

本发明公开了一种气溶胶光学厚度遥感反演方法，具体为：1)根据辐射传输模型建立查找表；2)表观反射率计算；3)依据纯像元指数提取遥感图像中的纯像元，生成纯像元掩膜；4)确定纯像元地表反射率；5)根据查找表，将遥感观测的辐射亮度反演为气溶胶光学厚度。本发明的气溶胶光学厚度遥感反演方法可以有效的用于城区气溶胶监测，除了受季节影响较小外，还有许多优势，如不受城市中高大建筑物的遮挡、不受混合像元影响、更容易获得纯像元的双向反射特征等等。

Description

一种气溶胶光学厚度遥感反演方法

技术领域

本发明涉及遥感领域，尤指一种气溶胶光学厚度遥感反演方法。

背景技术

利用卫星数据进行气溶胶的反演工作已经开展了三十多年。由于海洋表面相对均一，反射率很小并且近似为常数，因此，早期气溶胶直接反演研究和应用主要集中海洋和大的水体表面上空，并投入NOAA/AVHRR产品的业务应用，目前已经发展了两代算法，正在发展第三代算法(Husaret al.，1997)。而在陆地上空，由于气溶胶和地表反射率在时间和空间上的高度不均一性，且陆地地表反射率相对来说较大。另外，大气层顶辐射对气溶胶和地表反射率都有较强的敏感性，因此，很难从卫星对地观测信号中把气溶胶和地表的贡献定量区别开来，以提取气溶胶的光学厚度和地表反射率。

应用于MODIS传感器多波段数据的浓密植被气溶胶反演方法自TERRA卫星发射后已经广泛应用于全球气溶胶光学厚度的观测。其算法原理是：MODIS的通道1(红620～670nm)和通道3(蓝459～479nm)的地表反射率与通道7(近红外2105～2155nm)观测到的表观反射率在密集植被地区呈现良好的线性相关，而且近红外通道7的观测基本不受气溶胶的影响，因此，利用这一通道的数据区分浓密植被，并得到通道1和通道3在浓密植被地表的反射率。2007年MODIS算法(C5)对可见光与近红外波段的地表反射率关系进行了修正，以及在全球范围内使用了更多样化更精确的气溶胶模式后(基于全球气溶胶地基观测网AERONET得到)，反演精度较之前的C4版本有了大幅度的提升(见下图)，但是，升级后的算法在部分城市地区的反演结果并不理想。图1a和图1b、图2a和图2b是Li etal.，(2007)针对MODIS的C4算法和新C5算法在中国不同地区做的验证分析。左图是香河的验证结果，右图是北京城区的验证结果。从图中可以看出，新算法在非城区的精度相当高，拟合斜率达到了1.0，但是在城市地区却较差，拟合斜率只有0.38，相关性也较差。

该算法在城市地区失效是由于在城市地区难以找到足够的浓密植被像元，尤其是在冬季，算法中可见光与中红外地表反射率之间的关系已经不存在。而且，即使浓密植被算法在城市地区是可行的，但是，由于环境一号卫星CCD相机没有中红外波段，无法精确得到浓密植被的地表反射率，因此，该算法不可用于环境一号卫星。

发明内容

针对现有技术存在的问题，本发明的目的在于提供一种气溶胶光学厚度遥感反演方法。

为实现上述目的，本发明的气溶胶光学厚度遥感反演方法，具体为：1)根据辐射传输模型建立查找表；2)表观反射率计算；3)依据像元纯度指数提取遥感图像中的纯像元，生成纯像元掩膜；4)确定纯像元地表反射率；5)根据查找表，将遥感观测的辐射亮度反演为气溶胶光学厚度。

进一步，所述遥感图像由HJ-1 CCD拍摄。

进一步，该方法还包括步骤5)对整幅影像进行内插，得到光学厚度图。

进一步，步骤1)采用6S辐射传输程序进行辐射传输计算得出多组不同的大气参数(R_r+a、T_g(θ_s，θ_v)、T(θ_s)·T(θ_v)和S)组合而成的查找表。

进一步，步骤3)具体为：选取一个时间段内的晴朗日影像和污染日影像，假设在这个时间段内纯像元地表反射特性变化较小(Liang et al.，2006)，对晴朗天气的影像进行简单的大气校正，获取纯像元的地表反射率。

进一步，步骤3)选取的时间段小于3个月。

进一步，步骤4)具体为：根据读取的太阳天顶角、太阳方位角、观测天顶角、观测方位角，在查找表选取相应的数据，进行线性插值，得到不同波段、不同气溶胶光学厚度的R_r+a、T_g(θ_s，θ_v)、T(θ_s)·T(θ_v)和S大气参数；结合步骤3)得到的纯像元地表反射率，按照公式

$R^{*}({\mathrm{\θ}}_s,{\mathrm{\θ}}_v,{\mathrm{\φ}}_s-{\mathrm{\φ}}_v)=T_g({\mathrm{\θ}}_s,{\mathrm{\θ}}_v)[R_{r+a}({\mathrm{\θ}}_s,{\mathrm{\θ}}_v,{\mathrm{\φ}}_s-{\mathrm{\φ}}_v)+T\left({\mathrm{\θ}}_s\right)T\left({\mathrm{\θ}}_v\right)\frac{R}{1-\mathrm{RS}}]$ 计算每个预设气溶胶光学厚度对应的表观反射率，与卫星真实测量的表观反射率进行比较，相差最小的气溶胶光学厚度值即为所求，θ_s是太阳天顶角，φ_s传感器天顶角，θ_v是太阳方位角，φ_v是传感器方位角，T_g是气体吸收透过率，R_r+a是由分子散射加气溶胶散射所构成的路径程辐射反射率，T(θ_s)、T(θ_v)分别为太阳——目标、目标——传感器大气路径透射率，S为大气下届半球反照率，R是地物目标反射率。

本发明的气溶胶光学厚度遥感反演方法可以有效的用于城区气溶胶监测，除了受季节影响较小外，还有许多优势，如不受城市中高大建筑物的遮挡、不受混合像元影响、更容易获得纯像元的双向反射特征等等。

附图说明

图1a为MODIS C4算法验证图；

图1b为MODIS C5算法验证图；

图2a为MODIS C5算法在非城区(香河)验证图；

图2b为MODIS C5算法在城区(北京)验证图；

图3为对整幅影像进行内插，得到光学厚度图(2009年10月10日气溶胶光学厚度图(550nm))；

图4为本发明的方法AOD反演结果与AERONET对比验证图；

图5为MODIS气溶胶产品与AERONET对比验证图。

具体实施方式

与现有的气溶胶反演的传感器相比，环境一号卫星CCD相机具有极高的空间分辨率(30米)，这样，相比于浓密植被像元，在影像上能找到更多的纯像元，这些纯像元可能包含植被，但同时也包括城市中的大型广场(如北京的天安门广场)、楼顶、道路以及湖水，而后面这些纯像元在冬季也都是存在的。环境一号卫星CCD相机的另一个优势是它的高时间分辨率，尤其是将A星和B星结合起来的时候这一优势更加明显。因此，针对同一目标在一定时间段内可获取一系列观测值，其中包含了晴朗天气的观测以及污染天气的观测。已有多种气溶胶光学厚度反演算法假设地表反射率在一定时间段内(如三个月内，如果有足够的观测可考虑一个月内)变化较小，通过该时间段内的晴朗清晰影像获取目标的地表反射率，如结构函数法(Tanréet al.，1988)、深蓝算法(Hsu et al.，2004)、及Liang et al.，(2006)提出的亮地表地区MODIS改进算法。目标的地表反射率可通过晴朗天气的观测值通过大气校正得到，因此污染天气的气溶胶可基于该地表反射率反演得到。使用纯像元除了受季节影响较小外，还有许多优势，如不受城市中高大建筑物的遮挡、不受混合像元影响、更容易获得纯像元的双向反射特征等等。

因此，本方法充分发挥环境一号卫星CCD相机数据高空间分辨率(30m)的优势，利用图像中的纯像元，从观测信号中分离出地表的贡献，得到气溶胶光学厚度。

假定地表为朗伯体的情况下，传感器接收到的表观反射率R^*定义为：

$R^{*}=\frac{\mathrm{\πL}}{F_0\cos {\mathrm{\θ}}_s}---\left(1\right)$

式中，L为表观辐亮度，F₀外大气层外太阳辐照度。

传感器接收到的反射率为(Vermote E.et al，1997)：

$R^{*}({\mathrm{\θ}}_s,{\mathrm{\θ}}_v,{\mathrm{\φ}}_s-{\mathrm{\φ}}_v)=T_g({\mathrm{\θ}}_s,{\mathrm{\θ}}_v)[R_{r+a}({\mathrm{\θ}}_s,{\mathrm{\θ}}_v,{\mathrm{\φ}}_s-{\mathrm{\φ}}_v)+T\left({\mathrm{\θ}}_s\right)T\left({\mathrm{\θ}}_v\right)\frac{R}{1-\mathrm{RS}}]---\left(2\right)$

式中：θ_s是太阳天顶角，φ_s传感器天顶角，θ_v是太阳方位角，φ_v是传感器方位角，T_g是气体吸收透过率，R_r+a是由分子散射加气溶胶散射所构成的路径程辐射反射率，T(θ_s)、T(θ_v)分别为太阳——目标、目标——传感器大气路径透射率，S为大气下届半球反照率，R是地物目标反射率。

反演流程：

1)根据辐射传输模型建立查找表；查找表是通过设定不同卫星观测几何参数(9个太阳天顶角、6个观测天顶角、16个太阳与卫星之间的相对方位角)，不同的大气气溶胶参数(大陆型气溶胶模式、6个0.55μm处大气气溶胶光学厚度值：即0、0.25、0.5、1.0、1.5和1.95)，考虑要观测数据所在的波段，并考虑不同地表类型等参数，使用6S辐射传输程序进行辐射传输计算得出多组不同的大气参数(R_r+a、T_g(θ_s，θ_v)、T(θ_s)·T(θ_v)和s)组合而成的查找表。

2)表观反射率计算；根据定标系数将影像DN值转换为表观幅亮度，进而根据公式(1)得到表观反射率。

3)依据纯像元指数提取纯像元，生成纯像元掩膜。使用ENVI软件中的像元纯度指数(PPI，Pixel Purity Index)来提取影像上的纯像元。PPI是在高光谱数据中寻找“光谱最纯”像元的方法。传感器测量的光谱辐射值是非负的，由这些离散的辐射光谱构成的向量位于多维空间中，在该空间中形成一个凸的多维区域。凸多维区域的边界对应于纯像元(混合模型的最终光谱单元)，其它位于内部的向量可以认为是由这些边界单元线性混合而成。具体做法是，将N维散点图重复投影到单位随机向量上去。找出每次投影的极值像元，并对这些极值像元累计计数，产生一个像元纯度指数图，图中像元的值代表该像元在投影过程中为极值的次数。纯像元就是像元纯度指数图中数值较大的那些像元(Boardman et al.，1993)。该方法也被成功用于提取多光谱影像上的纯像元(胡妹婧等，2010)。

4)确定纯像元地表反射率：选取一个时间段内(小于3个月)的晴朗日影像和污染日影像，假设在这个时间段内纯像元地表反射特性变化较小，对晴朗天气的影像进行大气校正，获取纯像元的地表反射率，根据地表反射率反演污染日影像纯像元的气溶胶光学厚度；

5)依据辐射传输查找表，将遥感观测的辐射亮度反演为气溶胶光学厚度。根据读取的太阳天顶角、太阳方位角、观测天顶角、观测方位角，在查找表选取相应的数据，进行线性插值，得到不同波段、不同气溶胶光学厚度的R_r+a、T_g(θ_s，θ_v)、T(θ_s)·T(θ_v)和s等大气参数。结合上一步得到的纯像元地表反射率，按照公式(2)计算每个预设气溶胶光学厚度对应的表观反射率，与卫星真实测量的表观反射率进行比较，相差最小的气溶胶光学厚度值即为所求。

6)对整幅影像进行内插，得到光学厚度图。

如图3的2009年10月10日气溶胶光学厚度图显示，北京南部地区的气溶胶光学厚度比北部大，这一点与其他学者给出的结论是一致的。图中还可以看出，城中心的气溶胶污染比朝阳、海淀、通州等地区更严重，这可能与城中心人口更密集、汽车尾气排放更多有关系。AERONET北京站点这一天的地基观测AOD是0.8838(550nm)，利用本方法反演得到的AOD为0.9318(过境时间与地基观测时间相差20分钟以内)，绝对误差小于0.05。结果图有效的刻画了城市地区气溶胶污染源的空间变化情况，例如，图中右上角有一处气溶胶光学厚度较大，该污染源对应的是首都国际机场，显示机场上空严重的大气污染。

为了评估反演方法的性能，使用AERONET的地基观测数据进行对比验证。同时为了比对，也开展了过境时间接近的MODIS气溶胶产品与AERONET数据的对比验证。所有地基数据与卫星数据的时间相差都在30分钟以内，以保证地基观测与卫星过境之间大气状况的稳定。采用2009年8月至10月一共约3个月的数据进行验证，其中有32天无AERONET测量，另有32天没用无云的HJ星观测，剩下的16天中有4天该验证点无MODIS气溶胶产品，最后得到12天的数据用于验证。利用本方法反演的气溶胶AOD和MODIS气溶胶产品与AERONET的对比验证拟合结果分别见下图4和图5。验证图中给出了线性拟合的斜率和截距，可以看出本次验证数据中，本方法的拟合斜率与方差均优于MODIS。但是相关系数(R)略低，分别为0.84和0.88。本方法反演的AOD的均方根误差为0.15，低于MODIS的均方根误差(0.27)。验证数据显示，针对本次验证数据，本方法反演得到的AOD的误差均小于MODIS。综合来看，利用本方法，环境一号卫星CCD数据可以有效的用于城区气溶胶监测，为大气校正提供输入。

Claims (4)

1.一种气溶胶光学厚度遥感反演方法，具体为：1)根据辐射传输模型建立查找表；2)表观反射率计算；3)依据纯像元指数提取遥感图像中的纯像元，生成纯像元掩膜；4)确定纯像元地表反射率；5)根据查找表，将遥感观测的辐射亮度反演为气溶胶光学厚度；

其中，步骤3)具体为：选取一个时间段内的晴朗日影像和污染日影像，假设在这个时间段内纯像元地表反射特性变化较小，对晴朗天气的影像进行简单的大气校正，获取纯像元的地表反射率；

步骤4)具体为：根据读取的太阳天顶角、太阳方位角、观测天顶角、观测方位角，在查找表选取相应的数据，进行线性插值，得到不同波段、不同气溶胶光学厚度的R_r+a、T_g(θ_s,θ_v)、T(θ_s)·T(θ_v)和S大气参数；结合步骤3)得到的纯像元地表反射率，按照公式

$R^{*}({\mathrm{\θ}}_s,{\mathrm{\θ}}_v,{\mathrm{\φ}}_s-{\mathrm{\φ}}_v)=T_g({\mathrm{\θ}}_s,{\mathrm{\θ}}_v)[R_{r+a}({\mathrm{\θ}}_s,{\mathrm{\θ}}_v,{\mathrm{\φ}}_s-{\mathrm{\φ}}_v)+T\left({\mathrm{\θ}}_s\right)T\left({\mathrm{\θ}}_v\right)\frac{R}{1-\mathrm{RS}}]$

计算每个预设气溶胶光学厚度对应的表观反射率，与卫星真实测量的表观反射率进行比较，相差最小的气溶胶光学厚度值即为所求，θ_s是太阳天顶角，φ_s传感器天顶角，θ_v是太阳方位角，φ_v是传感器方位角，T_g是气体吸收透过率，R_r+a是由分子散射加气溶胶散射所构成的路径程辐射反射率，T(θ_s)、T(θ_v)分别为太阳——目标、目标——传感器大气路径透射率，S为大气下届半球反照率，R是地物目标反射率。

2.如权利要求1所述的气溶胶光学厚度遥感反演方法，其特征在于，该方法还包括步骤5)对整幅影像进行内插，得到光学厚度图。

3.如权利要求1所述的气溶胶光学厚度遥感反演方法，其特征在于，步骤1)采用6S辐射传输程序进行辐射传输计算得出多组不同的大气参数(R_r+a、T_g(θ_s,θ_v)、T(θ_s)·T(θ_v)和S)组合而成的查找表。

4.如权利要求1所述的气溶胶光学厚度遥感反演方法，其特征在于，步骤3)选取的时间段小于3个月。