基于高分一号卫星的陆地气溶胶光学性质反演方法

技术领域

本发明涉及大气遥感技术领域，尤其涉及一种基于高分一号卫星宽覆盖CCD 数据的溶胶光学性质反演方法。

背景技术

卫星遥感气溶胶的研究始于上世纪七十年代中期，在上个世纪最后的30多年 里，用于研究和形成业务化使用的卫星遥感气溶胶包括在火山喷发气溶胶监测、 NOAA(美国国家海洋和大气管理局)系列卫星对深海上空来自沙尘暴和森林火 灾造成的气溶胶层光学厚度遥感、红外遥感沙尘性气溶胶遥感、臭氧总量成像光 谱仪(TOMS)紫外波段对吸收性气溶胶的遥感等几个方面。美国地球观测系统 计划(EOS)Terra和Aqua携带的中分辨率成像光谱仪(MODIS)对全球发布了 10千米气溶胶光学厚度的全球分布产品；法国自1996年起推动了偏振多角度相机 (POLDER)探测海洋和陆地气溶胶的研究。到目前为止，利用卫星遥感气溶胶 已经形成了一定的全球气溶胶探测体系。

目前的较为成熟和应用广泛的陆地气溶胶卫星遥感反演算法为暗目标方法。 该方法多利用短波红外波段(如2.1微米附近)识别暗目标并获得地表反射率， 进而依据该波段与红、蓝波段的地表反射率的经验关系，反演气溶胶光学厚度。 然而GF1等国产卫星的CCD数据缺少短波红外波段，难以适用上述暗目标像元方 法。本发明提出一种适用于高分一号(下简称GF1)卫星宽覆盖CCD数据的气溶 胶光学厚度反演方法，可突破短波红外数据的限制，为基于国产高分辨率卫星进 行空气质量监测提供技术支持。

发明内容

本发明要解决的技术问题是：提供一种基于GF1卫星宽覆盖CCD数据的运行 高效、精度较好的区域陆地气溶胶遥感反演的方法。

为解决上述问题，本发明提供了一种基于GF1卫星的陆地气溶胶光学性质反 演方法，该方法包括步骤：

S1、生成大气光学参数查找表：根据我国不同区域及季节下适用的多种气溶 胶模式和大气模型、及GF1卫星宽覆盖CCD数据红、蓝波段的光谱响应函数， 用大气辐射传输模拟方法计算不同气溶胶光学厚度、不同观测几何条件以及GF1 卫星宽覆盖CCD数据红、蓝波段对应的大气光学参数，形成多维查找表；

S2、生成地表反射率先验数据库：基于其他卫星获取的覆盖全国、随时间更 新的大区域地表反射率产品，计算得出红、蓝波段先验地表反射率比值，形成地 表反射率先验数据库；

S3、选取适宜反演的暗目标像元：基于GF1卫星宽覆盖CCD数据的辐射定 标系数计算红、蓝波段表观反射率；基于云像元判据去除云像元、基于水体像元 判据去除水体像元、基于暗目标像元判据提取出适宜气溶胶光学厚度反演的暗目 标像元；

S4、计算暗目标像元观测几何条件：基于GF1卫星宽覆盖CCD数据的辅助 定位数据，计算每个暗目标像元的观测几何条件；

S5、插值获取暗目标像元的大气光学参数：根据待反演的GF1卫星宽覆盖 CCD数据影像的过境时间和区域，确定匹配的气溶胶模式及大气模型；然后针 对每一个暗目标像元，结合其波段、观测几何等索引条件从S1所述查找表中插 值得出不同气溶胶光学厚度对应的大气光学参数；

S6、基于辐射传输反演气溶胶光学厚度：针对S3所述每个暗目标像元，将 S3所述红、蓝波段表观反射率、S4所述观测几何条件、以及S5所述大气光学 参数带入辐射传输方程，计算不同气溶胶光学厚度下红、蓝波段地表反射率比值； 拟合建立该像元气溶胶光学厚度与辐射传输计算的红、蓝波段地表反射率比值的 经验关系；从S2所述数据库中提取所述暗目标像元的所述红、蓝波段先验地表 反射率比值，依所述经验关系插值得到所述暗目标像元此时的气溶胶光学厚度。

在上述方法中，优选地，所述GF1卫星宽覆盖CCD数据为经过几何粗校正 的蓝、绿、红、近红外等四个波段，其辅助数据包括：辐射标定系数、过境时间、 空间定位及观测几何等信息。

在上述方法中，优选地，所述多维查找表包含气溶胶模式×大气模型×观测 波段×观测几何×光学厚度的多维度索引条件，及每一组索引条件对应的大气光 学参数，并存储在数据文件中。

优选地，所述随时间更新的大区域地表反射率产品每次更新间隔的时间不超 过30天；所述大区域地表反射率产品所涉及的波段包含GF1卫星红、蓝波段， 或与所述GF1卫星红、蓝波段相近。优选地，所述地表反射率先验数据库包含 具有常见投影坐标的格网化数据，所述数据包含每个格网的红、蓝波段地表反射 率比值。

为提高数据处理质量，本方法的另一方案进一步包含卫星数据质量检查和预 处理的步骤，包括，针对GF1卫星宽覆盖CCD数据，检查卫星影像缺失或条带 影响情况，影像中随机噪声分布及程度，以及是否有明显的几何畸变；基于数据 质量标识去除数据缺失、条带和其他噪声像元数据。本发明的方法适用于GF1 卫星宽覆盖CCD的数据特点，运行高效且反演效果较好，有利于实现区域高分 辨率大气气溶胶光学性质的反演，为颗粒物污染监测提供了有益的数据源。

附图说明

图1是基于GF1卫星宽覆盖CCD数据的陆地气溶胶光学性质反演方法的数据 处理流程图

具体实施方式

本发明提出的基于GF1卫星宽覆盖CCD数据的气溶胶光学性质反演方法， 结合附图及实施例详细说明如下。

步骤S1、生成大气光学参数查找表：根据不同季节、区域，基于气溶胶物理 化学特性等先验知识，确定适用于我国的多种气溶胶模式和大气模型；通过大气 辐射传输软件模拟计算不同气溶胶光学厚度、不同观测几何条件以及GF-1卫星 宽覆盖CCD数据红、蓝波段对应的大气光学参数，存储于数据文件中形成多维 查找表。步骤S1具体包括：

S11、根据我国不同区域气溶胶来源、大气状况及其季节变化，结合气溶胶 物理、光学性质等先验知识，划定我国不同区域及季节下适宜的气溶胶模式与大 气模型；其中，气溶胶模式包括：大陆型、海洋型、城市型、煤烟型、沙尘型、 混合型等；大气模型包括热带、中纬度夏季、中纬度冬季、高纬度夏季、高纬度 冬季等；根据GF1卫星宽覆盖CCD数据红、蓝波段的光谱响应函数，在合理的 值域范围内设定不同观测几何条件(本实施例中为观测天顶角、太阳天顶角以及 相对方位角)以及不同的气溶胶光学厚度，如：1)设置查找表计算的波段为GF1 卫星宽覆盖CCD的红、蓝波段；2)设定不同的观测条件：选定10个太阳天顶 角(0°～90°，间隔10°)，10个观测天顶角(0°～90°，间隔10°)，19个相对方位 角(0°～180°，间隔10°)；3)设定不同的气溶胶光学厚度：相对于0.55μm波长 处的气溶胶光学厚度设为10个等级(0、0.25、0.5、0.75、1、1.25、1.5、1.75、 2.5和3)；

S12、选用适宜的辐射传输模式(如Second Simulation of the Satellite Signal in  the Solar Spectrum，6S)，对S11和S12中各类参数形成的索引条件进行模拟运 算，得到每一组索引条件对应的大气光学参数(本实施例中为大气路径辐射ρ₀、 大气整层双程透过率T、大气下行半球反射率S的三个参数)，并存储在数据文 件中形成查找表。

步骤S2、生成地表反射率先验数据库：基于其他卫星获取的大区域地表反射 率产品计算得出可覆盖全国、时间更新频次较高的红、蓝波段地表反射率比值先 验数据库。首先，选用适宜的卫星地表反射率产品以满足：1)可完整覆盖我国， 2)更新频次较高(如不超过30天)，3)包含GF-1卫星红、蓝波段或与其相近 (如选用MODIS的地表反射率产品MOD09)；然后，根据所选卫星地表反射率 产品的空间定位信息，将其重新处理具有常见投影坐标(如等经纬投影)的格网 化数据(格网大小为0.01度)；计算每个格网的红、蓝波段先验地表反射率比值， 形成数据库。

步骤S3、选取适宜反演的暗目标像元：对GF-1卫星宽覆盖CCD数据进行 如下质量检查与预处理；基于经验阈值提取出适宜气溶胶光学厚度反演的暗目标 像元。步骤S3具体包括：

S31、首先，对卫星影像进行质量检查，若满足：1)卫星数据缺失或条带影 响的像元数不超过整景图像的20％，2)数据中无显著的随机噪声，3)数据无明 显的几何畸变，则保留该景数据进入后续反演环节；然后，基于GF1卫星宽覆 盖CCD数据的质量标识等辅助数据，去除数据缺失或受噪声影响严重的像元； 最后，基于辅助数据中的辐射定标系数计算各波段表观反射率：从辅助数据中读 取辐射定标系数g和L₀，将遥感影像像元亮度DN值转换为表观辐亮度L：

L＝DN/g+L₀

然后将所述表观辐亮度L转化为该波段的表观反射率：

$\rho =\frac{\mathrm{\pi L}}{E_{\lambda }\cos {\theta }_s}$

其中，E_λ为大气层顶的太阳辐照度，可通过查询太阳常数获取；θ_s为太阳 天顶角，可从GF1卫星宽覆盖CCD的辅助数据中获得。

S32、依据特定的经验阈值，基于GF1卫星宽覆盖CCD的多波段表观反射 率数据去除云覆盖像元，所用云像元判据为：

ρ_red＜0.25且ρ_nir＜0.2

其中，ρ_red和ρ_nir分别表示GF1卫星宽覆盖CCD红、近红外两波段表观反 射率；

去除水体像元，所用水体像元判据为：

ρ_nir<ρ_red

保留的像元被认为是有效陆地像元，将进入后续筛选环节。

S33、基于暗目标像元判据筛选出适用于气溶胶光学厚度反演的陆地暗目标 像元。作为实施例，使用对浓密植被等暗目标像元敏感的归一化植被指数 (normalized difference vegetation index,NDVI)来设定阈值：

$\mathrm{NDVI}=\frac{{\rho }_{\mathrm{nir}}-{\rho }_{\mathrm{red}}}{{\rho }_{\mathrm{nir}}+{\rho }_{\mathrm{red}}}$

筛选出适用于气溶胶光学厚度反演的陆地暗目标像元的判据：NDVI大于 0.3。

步骤S4、计算暗目标像元观测几何条件：基于GF1卫星宽覆盖CCD的辅助 定位数据，计算每个暗目标像元的观测几何条件；步骤S4具体包括：

从辅助数据中获取每个暗目标像元的经纬度，在此基础上计算太阳天顶角和 太阳方位角。如太阳天顶角θ_s可依下式计算：

cosθ_s＝sin(lat)sinδ+cos(lat)cosδcost

太阳方位角可依下式计算：

其中，lat为该暗目标像元的纬度，δ为太阳光与地球赤道平面的夹角，t为 太阳的时角，定义纬地方时12点的时间为0，6点时为-π/2，18点时为π/2；

GF1卫星宽覆盖CCD的辅助数据中的观测方位角以及上面计算的太阳方 位角计算二者相对方位角：

步骤S5、插值获取暗目标像元的大气光学参数：根据待反演的GF1卫星宽 覆盖CCD影像的过境时间和区域，确定匹配的气溶胶模式及大气模型；然后针 对每一个暗目标像元，结合其波段、观测几何等索引条件从S3所述查找表中插 值得出不同气溶胶光学厚度对应的大气光学参数。步骤S5具体包括：

根据所在GF1卫星宽覆盖CCD数据覆盖区域和过境时间，结合后向轨迹等 气象信息，确定适宜反演的气溶胶模式与大气模型；然后针对每一个暗目标像元， 依据波段、太阳天的天顶角、观测天顶角以及相对方位角，从S13所述查找表中 插值提取出各个不同气溶胶光学厚度对应的大气光学参数(如本实施例包括大气 路径辐射、大气整层双程透过率和大气下行半球反射率)。

步骤S6、基于辐射传输反演气溶胶光学厚度：针对每个暗目标像元，将S3 所述红、蓝波段表观反射率以及S5所述大气光学参数带入辐射传输方程，计算 不同气溶胶光学厚度下红、蓝波段地表反射率比值；从S2所述数据库中提取该 暗目标像元红、蓝波段先验地表反射率比值，拟合得到该暗目标像元此时的气溶 胶光学厚度。步骤S6具体包括：

S61、针对每一个暗目标像元，依据S32中蓝、红波段表观反射率以及S5中 获取的大气光学参数，在平面平行大气假设下，根据下式计算该像元在不同气溶 胶光学厚度下蓝、红波段的地表反射率：

$\rho ={\rho }_0({\mu }_s,{\mu }_v,\mathrm{\Delta \phi })+\frac{T\left({\mu }_s\right)T\left({\mu }_v\right)r}{[1+\mathrm{rS}]}$

其中，μ_S＝cosθ_s，μ_V＝cosθ_v，θ_s与θ_v分别为太阳天顶角以及观测天顶角， Δφ为相对方位角；ρ为表观反射率，r为该波段地表反射率，S为大气向下的半 球反射率，T为大气整层双程透过率，ρ₀为大气的路径辐射项等效反射率；计算 得出该像元在不同气溶胶光学厚度下红、蓝波段的地表反射率，进而求得二者比 值：

k_i＝r_red,i/r_blue,i

其中，r_red,i、r_blue,i分别为第i个气溶胶光学厚度τ_i对应的红、蓝波段地表反 射率，k_i为二者比值。

S62、拟合建立该像元τ_i与k_i经验关系τ＝f(k)，从S22所述数据库中选出 该像元所在格网的红、蓝波段先验地表反射率比值k’，结合该拟合经验关系插值 得出该像元此刻的光学厚度τ′。

其中，所述GF1卫星宽覆盖CCD数据为经过几何粗校正的蓝、绿、红、近 红外等四个波段，存储为Tiff格式的影像文件；其辅助数据包括：辐射标定系数、 过境时间、空间定位及观测几何等信息，存储于与各波段影像文件对应的XML 文件中。

以上实施方式仅用于说明本发明，而并非对本发明的限制，有关技术领域的 普通技术人员，在不脱离本发明的精神和范围的情况下，还可以做出各种变化和 变型，因此所有等同的技术方案也属于本发明的范畴，本发明的专利保护范围应 由权利要求限定。