CN102103204B - 基于环境一号卫星的陆地气溶胶光学厚度反演方法 - Google Patents

Abstract

本发明公开了一种基于环境一号卫星的陆地气溶胶光学厚度反演方法，该方法包括步骤：S1.模拟观测几何、大气状况，进行辐射传输计算，获得大气参数，形成查找表；S2.根据辅助数据，提取各像元的观测天顶角以及方位角，并计算太阳天顶角以及方位角，获得观测几何；S3.对HJ-1A/B星的CCD数据进行预处理；S4.利用暗目标算法获得各波段地表反射率信息，将查找表插值到各像元的观测几何，根据经过预处理得到的表观反射率计算不同陆地气溶胶光学厚度下的地表反射率，然后拟合得到陆地气溶胶光学厚度。本发明的方法能够快速、方便的实现大范围的陆地气溶胶监测，为颗粒物污染监测提供了数据源。

Description

基于环境一号卫星的陆地气溶胶光学厚度反演方法

技术领域

本发明涉及大气遥感技术领域，尤其涉及一种基于环境一号卫星的陆地气溶胶光学厚度反演方法。

背景技术

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

目前的卫星遥感陆地气溶胶反演算法多为暗目标算法，该方法多利用短波红外波段识别暗目标并获得地表反射率进而反演得到气溶胶，但是短波红外波段数据的获得复杂，需对观测卫星的遥感传感器的要求很高，需进行特殊的波段设置，另一方便也导致其成本的增加。

发明内容

(一)要解决的技术问题

本发明要解决的技术问题是：提供一种基于环境一号卫星的可快速、方便的实现大范围的陆地气溶胶监测的陆地气溶胶光学厚度反演方法。

(二)技术方案

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

S1.根据环境一号HJ-1A/B星的CCD相机的特征，模拟各像元的观测几何、大气状况，进行辐射传输计算，获得相应的大气参数，存储于数据文件中形成查找表；

S2.根据环境一号HJ-1A/B星的CCD数据的辅助数据，提取各像元的观测天顶角以及观测方位角，并计算太阳天顶角以及太阳方位角，获得各像元的观测几何；

S3.对环境一号HJ-1A/B星的CCD数据进行如下预处理：辐射定标获得相机各波段的表观反射率、去除云覆盖像元、以及海陆分离提取陆地像元；

S4.利用暗目标算法获得各波段地表反射率信息，将所述查找表插值到各像元的观测几何，，根据经过预处理得到的表观反射率计算不同陆地气溶胶光学厚度下的地表反射率，拟合得到陆地气溶胶光学厚度。

其中，步骤S1进一步包括：

S1.1设定不同的观测情况：9个太阳天顶角，12个观测天顶角，16个相对方位角，气溶胶模式为大陆型气溶胶，相对于0.55μm波长处的气溶胶光学厚度设为6个等级，查找表计算的波段为环境一号HJ-1A/B星的CCD的蓝波段以及红波段，海拔设置为0m，所述相对方位角为太阳方位角与观测方位角之差；

S1.2调用辐射传输模式进行模拟运算，得到相应的大气程辐射、大气整层双程透过率、大气向下的半球反射率的三个参数，并存储在数据文件中形成查找表。

其中，步骤S2进一步包括：

S2.1根据环境一号HJ-1A/B星的CCD数据中的几何定位信息，计算每个像元的经纬度，得到太阳天顶角和太阳方位角，太阳天顶角θ_s为：

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

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

太阳方位角φ为：

$\cos \left(\mathrm{\φ}\right)=\tan (\mathrm{\π}/2-\mathrm{\θ})\tan \left(\mathrm{lat}\right)-\frac{\sin \mathrm{\δ}}{\cos {(\mathrm{\π}/2-\mathrm{\θ})}^{*}\cos \left(\mathrm{lat}\right)};$

S2.2根据环境一号HJ-1A/B星的CCD数据的辅助数据中的观测方位角以及所述太阳方位角，计算所述相对方位角。

其中，在步骤S3中，所述辐射定标获得表观反射率的方法包括步骤：

S3.11从所述辅助数据中读取辐射定标系数g和L₀，将遥感影像像元亮度DN值转换为表观辐亮度L：

L＝DN/g+L₀

S3.12根据所述太阳天顶角θ_s，将所述表观辐亮度归一化为各波段的表观反射率：

$\mathrm{\ρ}=\frac{\mathrm{\πL}}{E_{\mathrm{\λ}}\cos {\mathrm{\θ}}_s}$

其中，E_λ为大气层顶太阳辐照度；

其中，在步骤S3中，所述去除云覆盖像元时，选取0.2作为云去除阈值。

其中，在步骤S3中，所述海陆分离提取陆地像元是通过借助归一化植被指数NDVI阈值实现的：

$\mathrm{NVID}=\frac{{\mathrm{\ρ}}_{\mathrm{nir}}-{\mathrm{\ρ}}_{\mathrm{red}}}{{\mathrm{\ρ}}_{\mathrm{nir}}+{\mathrm{\ρ}}_{\mathrm{red}}}$

其中，ρ_nir和ρ_red分别表示环境一号HJ-1A/B星的CCD近红外波段和红波段的地表反射率，判断标准为：水体的NDVI小于0，土壤的NDVI为0.1，植被的NDVI大于0.3。

其中，步骤S4进一步包括：

S4.1使用NDVI识别出暗目标浓密植被；

S4.2将所述查找表插值到各像元的观测几何，得到不同陆地气溶胶光学厚度下插值后的大气参数；

S4.3根据步骤S3.12中得到的蓝波段和红波段的表观反射率以及下式：

μ_S＝cosθ_s，μ_V＝cosθ_v，

与μ_v分别为太阳天顶角以及观测天顶角，r为朗伯体地表反射率，S为大气向下的半球反射率，T为大气整层双程透过率，ρ₀为大气的路径辐射项等效反射率，

为相对方位角；

拟合得到的蓝波段和红波段地表反射率；

S4.4使步骤S4.3得到的蓝波段和红波段地表反射率符合蓝波段和红波段地表反射率的线性关系，从而得到的气溶胶光学厚度；所述暗目标在蓝波段和红波段地表反射率的线性关系表示为：

$r_{\mathrm{red}}^S={\mathrm{kr}}_{\mathrm{blue}}^S$

其中，

分别表示暗目标在红波段和蓝波段的地表反射率，k为红波段和蓝波段地表反射率比率。

其中，所述环境一号HJ-1A/B星的CCD数据为经过几何校正的蓝、绿、红、近红外四个波段，所述环境一号HJ-1A/B星的CCD数据的辅助数据包括：包含标定系数、过境时间辅助信息的xml文件，以及包含观测几何的文本文件。

(三)有益效果

本发明的方法从环境一号HJ-1A/B星的CCD数据中反演得到陆地气溶胶光学厚度，能够快速、方便的实现大范围的陆地气溶胶监测，为颗粒物污染监测提供了数据源。

附图说明

图1为依照本发明一种实施方式的基于环境一号卫星的陆地气溶胶光学厚度反演方法流程图。

具体实施方式

本发明提出的基于环境一号卫星的陆地气溶胶光学厚度反演方法，结合附图及实施例详细说明如下。

HJ-1A/B是我国的环境与灾害监测预报小卫星星座系统，具有较高空间、时间分辨率和宽覆盖特点，其数据不仅能为环境与减灾业务运行系统提供重要保障，还将成为很多部门日常业务的重要数据源。基于环境卫星数据建立的环境与减灾应用系统，对推动遥感卫星业务服务具有重要的示范作用，为我国环境空气质量遥感监测提供新的遥感平台。本发明的方法主要基于中国环境一号A/B星(本申请中统称HJ-1A/B星)CCD数据，利用HJ-1A/B星CCD数据从HJ-1A/B星CCD数据的辅助数据计算得到的太阳天顶角、相对方位角等信息，经过预处理得到红光和蓝光波段表观反射率，最后反演得到的气溶胶光学厚度分布图。从而，本发明能够快速、方便的实现大范围的陆地气溶胶光学厚度的监测，为颗粒物污染监测提供了数据源。

如图1所示，依照本发明一种实施方式的基于HJ-1A/B的陆地气溶胶光学厚度反演方法包括步骤：

S1.利用辐射传输模式(本实施方式中辐射传输模型为SecondSimulation of the Satellite Signal in the Solar Spectrum，6S模型)，针对HJ-1A/B卫星的CCD相机的特征，模拟各像元的观测几何、大气状况等情况，进行辐射传输计算，获得相应的大气参数，存储在数据文件中形成查找表；

S2.根据HJ-1A/B星的CCD相机数据的辅助数据，提取各像元的观测天顶角以及观测方位角，并计算太阳天顶角以及太阳方位角，获得各像元的观测几何；

S3.对HJ-1A/B星的CCD相机数据进行如下预处理：辐射定标获得表观反射率、去除云覆盖像元、以及海陆分离提取陆地像元；

S4.利用暗目标算法获得各波段地表反射率信息，将所述查找表插值到各像元的观测几何，根据经过预处理得到的表观反射率计算不同陆地气溶胶光学厚度下的地表反射率，拟合得到陆地气溶胶光学厚度。

其中，步骤S1进一步包括：

S1.1设定不同的观测情况：9个太阳天顶角，12个观测天顶角，16个相对方位角，气溶胶模式为大陆型气溶胶，相对于0.55μm波长处的气溶胶光学厚度设为6个等级(0、0.25、0.5、1、1.5和1.95)，查找表计算的波段为HJ-1A/B星的CCD相机的第一波段(蓝波段)以及第三波段(红波段)，海拔设置为0m，相对方位角为太阳方位角与观测方位角之差；

S1.2调用辐射传输模式进行模拟运算，得到相应的大气程辐射(大气的路径辐射项等效反射率)、大气整层双程透过率、大气下界的半球反射率的三个参数，并存储在数据文件中形成查找表。

步骤S2进一步包括：

S2.1根据HJ-1A/B星的CCD相机数据中的几何定位信息，计算每个像元的经纬度，得到太阳天顶角和太阳方位角，太阳天顶角θ_s为：

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

太阳方位角φ为：

$\cos \left(\mathrm{\φ}\right)=\tan (\mathrm{\π}/2-\mathrm{\θ})\tan \left(\mathrm{lat}\right)-\frac{\sin \mathrm{\δ}}{\cos {(\mathrm{\π}/2-\mathrm{\θ})}^{*}\cos \left(\mathrm{lat}\right)};---\left(2\right)$

S2.2根据HJ-1A/B星的CCD相机数据的辅助数据中的观测方位角计算其相对方位角。

在步骤S3中，辐射定标获得表观反射率的方法包括步骤：

S3.11从HJ-1A/B星的CCD相机数据的辅助数据中读取辐射定标系数g和L₀，将遥感影像像元亮度值(DN值)转换为表观辐亮度L：

L＝DN/g+L₀(3)

S3.12根据步骤S2得到的太阳天顶角θs，将表观辐亮度归一化为表观反射率：

$\mathrm{\ρ}=\frac{\mathrm{\πL}}{E_{\mathrm{\λ}}\cos {\mathrm{\θ}}_s}---\left(4\right)$

其中，E_λ为大气层顶太阳辐照度。

去除云覆盖像元时，在红光波段卫星观测到的云的反射率一般大于0.3，而一般地物的反射率要小于0.25，选取0.2作为云去除阈值。

海陆分离提取陆地像元是通过借助归一化植被指数(NDVI)阈值实现的：

$\mathrm{NVID}=\frac{{\mathrm{\ρ}}_{\mathrm{nir}}-{\mathrm{\ρ}}_{\mathrm{red}}}{{\mathrm{\ρ}}_{\mathrm{nir}}+{\mathrm{\ρ}}_{\mathrm{red}}}---\left(5\right)$

其中，ρ_nir和ρ_red分别表示HJ-1A/B星的CCD相机第四波段(近红外波段)和第三波段(红波段)的反射率。水体的NDVI一般小于0，土壤的在0.1左右，植被大于0.3。

对于植被密集的(即暗目标)像元红蓝通道反射率较低且与近红外通道反射率与具有很好的线性相关；气溶胶的影响在短波红外波段比在可见光波段小15-30倍。具体如下，在地表朗伯体、大气水平均一的假设条件下，表观反射率ρ为：

其中，μ_S＝cosθ_s，μ_V＝cosθ_v，

与μ_v分别为太阳天顶角以及观测天顶角，r为朗伯体地表反射率，S为大气下界的半球反射率，T为大气整层双程透过率，ρ₀为大气的路径辐射项等效反射率，

为相对方位角。S、ρ₀以及T(μ_S)T(μ_V)是3个代表大气状况的参数，遥感反演中可以从中获取所需的大气参数，即气溶胶光学厚度。

HJ-1A/B星CCD相机数据应用本发明的方法的关键问题是从红、蓝通道地表反射率的线性关系实现地气解耦。步骤S4进一步包括：

S4.1使用NDVI识别出暗目标浓密植被；

S4.2将查找表插值到各像元的观测几何，得到不同陆地气溶胶光学厚度下插值后连续的大气参数；

S4.3根据公式(4)得到的蓝波段和红波段的表观反射率以及公式(6)拟合得到的蓝波段和红波段地表反射率；

S4.4使步骤S4.3得到的蓝波段和红波段地表反射率符合公式(7)的线性关系，进而得到气溶胶光学厚度。

暗目标在蓝波段和红波段地表反射率的线性关系表示为：

$r_{\mathrm{red}}^S={\mathrm{kr}}_{\mathrm{blue}}^S---\left(7\right)$

其中，

分别表示红波段和蓝波段浓密植被的地表反射率，k为红蓝地表反射率比率，要根据HJ-1A/B星CCD相机的特征，结合地面观测数据设定。

上述HJ-1A/B星CCD相机数据为经过几何校正的二级产品，包括蓝、绿、红、近红外四个波段，HJ-1A/B星CCD相机数据的辅助数据包括：包含标定系数、过境时间等辅助信息的xml文件，以及包含观测天顶角、观测方位角等观测几何的文本文件。

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

Claims (6)

1.一种基于环境一号卫星的陆地气溶胶光学厚度反演方法，其特征在于，该方法包括步骤：

S1.根据环境一号HJ-1A/B星的CCD相机的特征，模拟各像元的观测几何、大气状况，进行辐射传输计算，获得相应的大气参数，存储于数据文件中形成查找表；

S2.根据环境一号HJ-1A/B星的CCD数据的辅助数据，提取各像元的观测天顶角以及观测方位角，并计算太阳天顶角以及太阳方位角，获得各像元的观测几何；

S3.对环境一号HJ-1A/B星的CCD数据进行如下预处理：辐射定标获得相机各波段的表观反射率、去除云覆盖像元、以及海陆分离提取陆地像元，其中，所述海陆分离提取陆地像元是通过借助归一化植被指数NDVI阈值实现的：

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

ρ_nir和ρ_red分别表示环境一号HJ-1A/B星的CCD近红外波段和红波段的地表反射率；

S4.利用暗目标算法获得各波段地表反射率信息，将所述查找表插值到各像元的观测几何，根据经过预处理得到的表观反射率计算不同陆地气溶胶光学厚度下的地表反射率，拟合得到陆地气溶胶光学厚度；

在步骤S3中，所述辐射定标获得表观反射率的方法包括步骤：

S3.1从所述辅助数据中读取辐射定标系数g和L₀，将遥感影像像元亮度DN值转换为表观辐亮度L：

L=DN/g+L₀

S3.2根据所述太阳天顶角，将所述表观辐亮度归一化为各波段的表观反射率：

$\mathrm{\ρ}=\frac{\mathrm{\πL}}{E_{\mathrm{\λ}}\cos {\mathrm{\θ}}_s}$

其中，θ_s为所述太阳天顶角，E_λ为大气层顶太阳辐照度；

步骤S4进一步包括：

S4.1使用NDVI识别出暗目标浓密植被；

S4.2将所述查找表插值到各像元的观测几何，得到不同陆地气溶胶光学厚度下插值后的大气参数；

S4.3根据步骤S3.2中得到的蓝波段和红波段的表观反射率以及下式：

μ_S=cosθ_s，μ_v=cosθ_v，μ_s与μ_v分别为太阳天顶角的余弦值以及观测天顶角的余弦值，r为朗伯体地表反射率，S为大气向下的半球反射率，T为大气整层双程透过率，ρ₀为大气的路径辐射项等效反射率，

为相对方位角；

拟合得到的蓝波段和红波段地表反射率；

S4.4使步骤S4.3得到的蓝波段和红波段地表反射率符合蓝波段和红波段地表反射率的线性关系，从而得到气溶胶光学厚度；暗目标在蓝波段和红波段地表反射率的线性关系表示为：

$r_{\mathrm{red}}^S={\mathrm{kr}}_{\mathrm{biue}}^S$

其中，

分别表示暗目标在红波段和蓝波段的地表反射率，k为红波段和蓝波段地表反射率比率。

2.如权利要求1所述的基于环境一号卫星的陆地气溶胶光学厚度反演方法，其特征在于，步骤S1进一步包括：

S1.1设定不同的观测情况：9个太阳天顶角，12个观测天顶角，16个相对方位角，气溶胶模式为大陆型气溶胶，相对于0.55μm波长处的气溶胶光学厚度设为6个等级，查找表计算的波段为环境一号HJ-1A/B星的CCD的蓝波段以及红波段，海拔设置为0m，所述相对方位角为太阳方位角与观测方位角之差；

S1.2调用辐射传输模式进行模拟运算，得到相应的大气程辐射、大气整层双程透过率、大气向下的半球反射率的三个参数，并存储在数据文件中形成查找表。

3.如权利要求2所述的基于环境一号卫星的陆地气溶胶光学厚度反演方法，其特征在于，步骤S2进一步包括：

S2.1根据环境一号HJ-1A/B星的CCD数据中的几何定位信息，计算每个像元的经纬度，得到太阳天顶角和太阳方位角，太阳天顶角θ_s为：

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

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

太阳方位角φ为：

$\cos \mathrm{\φ}=\tan (\mathrm{\π}/2-{\mathrm{\θ}}_s)\tan \left(\mathrm{lat}\right)-\frac{\sin \mathrm{\δ}}{\cos (\mathrm{\π}/2-{\mathrm{\θ}}_s)\cos \left(\mathrm{lat}\right)};$

S2.2根据环境一号HJ-1A/B星的CCD数据的辅助数据中的观测方位角以及所述太阳方位角，计算所述相对方位角。

4.如权利要求1所述的基于环境一号卫星的陆地气溶胶光学厚度反演方法，其特征在于，在步骤S3中，所述去除云覆盖像元时，选取0.2作为云去除阈值。

5.如权利要求4所述的基于环境一号卫星的陆地气溶胶光学厚度反演方法，其特征在于，在步骤S3中，所述海陆分离提取陆地像元是通过借助归一化植被指数NDVI阈值实现的判断标准为：水体的NDVI小于0，土壤的NDVI为0.1，植被的NDVI大于0.3。

6.如权利要求1-5任一项所述的基于环境一号卫星的陆地气溶胶光学厚度反演方法，其特征在于，所述环境一号HJ-1A/B星的CCD数据为经过几何校正的蓝、绿、红、近红外四个波段，所述环境一号HJ-1A/B星的CCD数据的辅助数据包括：包含标定系数、过境时间辅助信息的xml文件，以及包含观测几何的文本文件。

Images