阅读说明：本技术 一种考虑地表非朗伯特性的气溶胶光学厚度遥感反演方法 (Aerosol optical thickness remote sensing inversion method considering surface non-Lambert characteristics ) 是由 田信鹏 高志强 于 2019-09-05 设计创作，主要内容包括：本发明涉及一种考虑地表非朗伯特性的气溶胶光学厚度遥感反演方法,包括获取遥感数据进行预处理；利用遥感图像地表反射率数据基于先验知识约束完成半经验核驱动模型的每日核系数估算；通过离散余弦变换的惩罚最小二乘估计时间序列平滑滤波算法构建单位时间内核系数地表先验知识数据集；建立气溶胶反演通用查找表；利用半经验核驱动模型结合核系数地表先验知识数据集,确定非朗伯体辐射传输前向模型的地表反射率参数；利用上述结果反演获取气溶胶光学厚度。本发明可获取高精度的空间分辨率为500m的气溶胶光学厚度产品,能提供空间相对更为连续的气溶胶分布,显示了在支持区域污染精细管控、重点城市污染传输通道监测及污染物溯源等领域的应用潜力。(The invention relates to an aerosol optical thickness remote sensing inversion method considering surface non-Lambert characteristics, which comprises the steps of obtaining remote sensing data for preprocessing; utilizing the remote sensing image earth surface reflectivity data to complete the estimation of the daily nuclear coefficient of the semi-empirical nuclear driving model based on the prior knowledge constraint; constructing a unit-time kernel coefficient earth surface prior knowledge data set by a punished least square estimation time sequence smoothing filtering algorithm of discrete cosine transform; establishing an aerosol inversion general lookup table; determining the earth surface reflectivity parameter of a non-Lambertian radiation transmission forward model by utilizing a semi-empirical kernel driving model in combination with a kernel coefficient earth surface prior knowledge data set; and (5) obtaining the optical thickness of the aerosol by utilizing the result inversion. The invention can obtain aerosol optical thickness products with high-precision spatial resolution of 500m, can provide aerosol distribution with relatively more continuous space, and shows application potential in the fields of supporting regional pollution fine control, monitoring of key urban pollution transmission channels, tracing of pollutants and the like.)

一种考虑地表非朗伯特性的气溶胶光学厚度遥感反演方法

技术领域

本发明涉及大气遥感领域，尤指一种考虑地表非朗伯特性的气溶胶光学厚 度遥感反演方法。

背景技术

气溶胶作为地球气候系统的重要组成部分，是大气中三大活跃成分之一， 其数量与性质在时间及空间尺度上的变化均是巨大的，导致难以用固定模式予 以描述。气溶胶光学厚度(Aerosol Optical Depth,AOD)作为表征气溶胶粒子消 光作用的重要光学参数之一，在地气系统辐射收支平衡、全球气候变化及大气 环境和人类健康等方面具有重要指示作用(张莹等，2013；王中挺等，2016)。 由于大气气溶胶具有较大的时空变异性，卫星遥感成为能够提供可靠气溶胶分 布信息的有效手段。卫星传感器接收到的大气层顶辐射信息是由地表及大气共 同作用、相互耦合的最终结果，并且大气贡献相比地表是十分有限的，尤其是 城市等高反射率区域。遥感定量化过程就是从这单一信息中，区分两个未知量， 其本质是一个病态反演问题。在实际的气溶胶遥感中，需要通过先验知识来确 定气溶胶类型及地表反射率参数，这两者的不确定性也正是气溶胶反演误差的 主要来源。气溶胶反演算法的发展与改进，也正是围绕气溶胶类型的假设及地 表反射率估算这两个方面。

近年来，国内外诸多学者在气溶胶反演中对地表方向反射的估算方面做了 许多努力。Thomas等提出了ORAC(Oxford-RAL Aerosol and Cloud)算法实现了 ATSR-2和AATSR数据的气溶胶参数反演，该算法是使用BRDF/Albedo产品构 建反射率先验知识库实现地气分离。Xue等通过中分辨率成像光谱仪(Moderate Resolution ImagingSpectroradiometer,MODIS)BRDF/反照率产品数据来增加方 程约束条件，基于AATSR星下点及前向两个角度观测数据提出了考虑BRDF效 应的气溶胶协同反演算法。在此基础上，Guang等利用MODIS数据提出了AOD 与地表反射率双星协同反演算法，进一步提高了亮地表上空气溶胶反演精度。 Qin等在假设地表BRDF形状稳定的基础上，基于比值法实现了AATSR数据澳 大利亚地区气溶胶及BRDF的反演。上述诸多算法均需要大气条件不变的时间 序列同步观测数据，导致其普适性较差，且使用的仍是朗伯体前向辐射传输模 型计算表观反射率，未从根本上解决地表非朗伯反射的影响。

发明内容

为解决传统气溶胶算法中地表朗伯反射假定对地表反射率估算带来的不确 定性，本算法发明提出一个耦合地表二向反射特性的气溶胶反演算法，最终构 建一个具有较强普适性、可业务化运行的气溶胶反演算法。

为实现上述目的，本发明所采用的技术方案是：一种考虑地表非朗伯特性 的气溶胶光学厚度遥感反演方法，包括：

获取遥感数据并进行预处理；

地表先验知识数据集的构建：利用长时间序列的遥感图像地表反射率数据， 基于先验知识条件约束完成半经验核驱动模型的每日核系数估算；通过离散余 弦变换的惩罚最小二乘估计DCT-PLS时间序列平滑滤波算法构建单位时间内的 核系数地表先验知识数据集；

查找表构建：构建通用大气参数数据集，再根据MODIS光谱响应函数、对 该大气参数数据集中的数据通过积分运算建立MODIS传感器气溶胶反演查找 表文件；

地表反射率参数的确定：利用RTLS半经验核驱动模型，结合构建的核系数 地表先验知识数据集SANIFD，确定非朗伯体辐射传输前向模型中四个地表反射 率参数；

气溶胶光学厚度反演：根据成像角度信息、地表反射率参数及查找表中大 气参数信息，结合预处理获取的表观反射率数据，基于非朗伯体辐射传输前项 模型，反演获取气溶胶光学厚度。

所述基于先验知识条件约束完成半经验核驱动模型核系数的估算，采用的 代价函数为：J(f)＝(Kf-ρ)^TΣ^T(Kf-ρ)+(f-f₀)^TM^T(f-f₀)；

其中，K为三个核矩阵，f为对应的三个核系数向量，ρ为方向反射率向量， Σ为反射率误差的协方差矩阵，f₀及M为核系数先验知识的均值向量及协方差 矩阵；

上述代价函数的最小二乘解表达式可写成：

f＝(K^TΣ^-1K+M^-1)^-1(K^TΣ^-1ρ+M^-1f₀)，通过该式即可完成核系数的估算。

所述半经验核驱动模型的核系数包括：各向同性散射f_iso、、体散射f_vol及几何 光学散射f_geo三个核系数。

所述单位时间内为每若干天内。

通过DCT-PLS时间序列平滑滤波算法构建核系数地表先验知识数据集，具 体实现公式为：

式中，|| ||表示欧几里得范数；表示待平滑的时间序列变量；W表示最佳 平滑滤波权重；s为控制滤波平滑度参数，随着s值的增大，会越平滑；对于 一维数据的计算，惩罚项P通常可用Laplace差分计算。

所述通用大气参数数据集记录如下数据：

不同观测几何数据：18个太阳天顶角、15个观测天顶角、19个相对方位角；

气溶胶光学厚度AOD条件(0.05-3.0范围内)下16个气溶胶光学厚度；

0.3-2.5μm光谱范围的大气参数信息，包括t_dd(θ_s)、t_dd(θ_v)、t_dh(θ_s)、t_hd(θ_v)、 及ρ₀。

所述非朗伯体辐射传输前向模型中四个地表反射率参数包括：r_dd、r_dh、r_hd、 r_hh，计算公式如下：

其中，K_vol及K_geo分别为体散射核及几何光学核，f_iso、f_vol及f_geo分别表示各 向同性散射、体散射及几何光学散射核核系数值，与波段λ有关；R为地表反射 率值，当分别计算时，r_dd、r_dh、r_hd、r_hh表示地表不同的反射特性，分别为地表 方向-方向反射率、方向-半球反射率、半球-方向反射率及半球-半球反射率。

非朗伯体辐射传输前项模型为：

其中，

矩阵T(θ_s)、T(θ_v)及分别定义为：

T(θ_s)＝[t_dd(θ_s)t_dh(θ_s)]

式中，θ_s、θ_v、分别为太阳天顶角、传感器天顶角、太阳方位角、 传感器方位角及相对方位角；ρ_TOA为表观反射率；ρ₀为大气路径散射反射率，表 示大气程辐射项；为大气球形反照率；t和r分别表示透过率和反射率，下标 h和d分别表示半球(漫射)和方向(直射)；因此，t_hd(θ_v)表征在卫星观测方向 上大气对地表漫射辐射的消光，即大气上行漫射透过率；t_dd(θ_v)表征在卫星观测 方向上大气对地表反射直射辐射的消光，即大气上行直射透过率；同样，t_dh(θ_s) 表征在入射方向上大气对天空漫射辐射的消光，即大气下行漫射透过率；t_dd(θ_s) 则为入射方向上大气对太阳直射辐射的消光，即大气下行直射透过率；r_dd、r_dh、 r_hd、r_hh表示地表不同的反射特性，分别为地表方向-方向反射率、方向-半球反射率、半球-方向反射率及半球-半球反射率。

对遥感数据进行预处理包括云检测、冰/雪以及水体识别，采用改进的动态 阈值法，具体表示为：

1)首先使用地表先验知识数据集，计算获取各波段地表反射率值： 然后，利用不同波段云检测动态阈值模型实现云像元的识别，公式表示为：

2)由于采用的半经验核驱动模型中几何光学核及体散射核不适用于冰/雪及 水体地表类型的反射率计算，在气溶胶反演中需要将冰/雪及水体进行掩膜剔除， 具体如下：

使用遥感图像预处理后的第2、5波段表观反射率值按照公式计算，识别冰/雪；使用遥感图像预处理后的第1、2波段 的地表反射率值按照公式计算，识别水体。

本发明的有益效果是：

本发明可以获取高精度的空间分辨率为500m的气溶胶光学厚度产品，且 能提供空间相对更为连续的气溶胶分布信息，显示了在支持区域污染精细管控、 重点城市污染传输通道监测及污染物溯源等领域的应用潜力。为了实现上述以 及相关目的，本发明的一个或多个方面同时包括后面的详细说明及权利要求书 中特别指出的特征。

附图说明

图1为本发明的气溶胶光学厚度反演方法流程示意图。

图2为本发明地表先验知识数据集构建流程示意图。

图3为本发明气溶胶产品各数据层示例图。

图4为本发明获取的气溶胶光学厚度与MODIS气溶胶产品空间分布对比图。

图5为本发明的方法AOD反演结果与MODIS各气溶胶产品(a：DT算法，b：DB算法，c：DTDB算法，d：MAIAC算法)在香河植被AERONET站点的 交叉对比验证图。

图6为本发明的方法AOD反演结果与MODIS各气溶胶产品(a：DT算法， b：DB算法，c：DTDB算法，d：MAIAC算法)在北京城市AERONET站点的 交叉对比验证图。

具体实施方式

下面将结合附图对本发明进行更全面的说明，将阐述本发明实现过程中涉 及到的重要具体细节。然而，本发明可以体现为多种不同形式，提供这些实施 例，是为使本发明更加全面和完整，并将本发明的范围完全地传达给本领域的 普通技术人员。

本发面采用一种考虑地表各向异性反射特性的快速、准确的辐射传输模型， 该模型将辐射过程分为上行和下行分量，并且每个分量中包含直射和漫射两部 分。这种划分可以很好的将入射到地表及地表反射辐射的各个反射率分量清楚 的区分确定，其表观反射率计算公式如下：

其中，

矩阵T(θ_s)、T(θ_v)及分别定义为：

T(θ_s)＝[t_dd(θ_s) t_dh(θ_s)]

式中，θ_s、θ_v、分别为太阳天顶角、传感器天顶角、太阳方位角、 传感器方位角及相对方位角；ρ_TOA为表观反射率；ρ₀为大气路径散射反射率，表 示大气程辐射项；为大气球形反照率；t和r分别表示透过率和反射率，下标 h和d分别表示半球(漫射)和方向(直射)；因此，t_hd(θ_v)表征在卫星观测方向 上大气对地表漫射辐射的消光，即大气上行漫射透过率；t_dd(θ_v)表征在卫星观测 方向上大气对地表反射直射辐射的消光，即大气上行直射透过率；同样，t_dh(θ_s) 表征在入射方向上大气对天空漫射辐射的消光，即大气下行漫射透过率；t_dd(θ_s) 则为入射方向上大气对太阳直射辐射的消光，即大气下行直射透过率；r_dd、r_dh、 r_hd、r_hh表示地表不同的反射特性，分别为地表方向-方向反射率、方向-半球反射率、半球-方向反射率及半球-半球反射率。

由式(1)可知，在计算表观反射率时需要估算与地表有关的反射率参数及 与大气有关的透过率、程辐射等大气参数。由于地表参数和大气参数的估算是 相互独立的，对于地表各反射率参数的估算可通过半经验核驱动模型实现，大 气参数则可以根据观测几何条件从辐射传输代码预先构建的LUT中获取，关于 这两个参数的获取方法在下面的反演流程中会详细阐述。

如图1所示，为本发明的气溶胶光学厚度反演方法流程示意图。流程如下：

(1)遥感数据的获取：本发明算法基于MODIS实现了高分辨率的气溶胶 光学厚度产品的反演生产，在算法构建及产品生产过程中需要的遥感数据主要 包括：1)MODISLevel-1空间分辨率为500m的表观反射率数据(MOD02 HKM)， 2)MODIS Level-1几何位置数据(MOD03)，3)MODIS Level-2空间分辨率为 500m的地表反射率产品数据(MOD09GA)及BRDF/反照率产品数据 (MCD43A1)，4)MODIS Level-3空间分辨率为500m的植被指数产品数据 (MOD13A1)，5)MODIS Level-2气溶胶产品数据(MOD04、MCD19A2)。以 上所有数据均可从LAADS DAAC免费下载获取 (https://ladsweb.modaps.eosdis.nasa.gov/)。

(2)数据处理：对MOD02HKM数据进行几何校正、云/冰雪检测等预处 理，获取可进行气溶胶反演区域的表观反射率数据。

(3)地表先验知识数据集的构建：在本发明算法中，地表先验知识数据集 为半经验核驱动模型中各向同性散射、体散射及几何光学散射三个核系数值， 构建过程包括两步：1)利用长时间序列的MOD09GA数据，基于先验知识条件 约束算法，实现半经验核驱动模型(公式(8))中f_iso、f_vol及f_geo三个核系数估 算(每日产品)；2)利用第一步估算的逐日核系数数据，基于离散余弦变换的 惩罚最小二乘估计(Discrete Cosine Transform-Penalized Least Square regression, DCT-PLS)时间序列平滑滤波算法合成构建最终每8天的f_iso、f_vol及f_geo核系数 地表先验知识数据集(Surface ANIsotropic FactorDatabase,SANIFD)，构建流程 如图2所示。

(4)查找表构建：基于大气辐射传输程序(MODTRAN 5)构建通用大气 参数数据集，并结合MODIS传感器光谱响应函数，构建其气溶胶反演查找表文 件。

(5)地表反射率的确定：基于第(3)步构建的地表先验知识数据集，结 合半经验核驱动模型来确定非朗伯体辐射传输前向模型中各地表反射率参数。

(6)气溶胶光学厚度反演：根据成像角度信息、地表反射率及查找表中大 气参数信息，结合第(2)步获取的表观反射率数据，反演获取气溶胶光学厚度。

具体反演流程：

1)几何校正及表观反射率计算。首先利用MOD03数据中的经纬度数据集 作为控制点，逐像元插值计算MOD02HKM对应的实际经纬度；然后，针对每 一个扫描带数据进行几何校正；最后，对所有扫描带校正结果进行拼接完成几 何校正，并同时实现数据的辐射定标获取真实表观反射率数据图像。

2)云检测。本发明采用改进的动态阈值云检测算法进行MODIS数据的云 检测，其基本过程为：首先使用地表先验知识数据集，计算获取公式(2)-(4) 中各波段地表反射率值然后，利用不同波段云检测动态阈值 模型模拟出云像元对应的表观发射率最后，利用真实表 观反射率数据与模拟表观反射率对比实现云像元的识别， 公式表示为：

3)冰/雪及水体识别。由于采用的半经验核驱动模型中几何光学核及体散射 核不适用于冰/雪及水体地表类型的反射率计算，在气溶胶反演中需要将这类地 表类型掩膜剔除。算法中采用的利用MODIS第2、5波段的表观反射率对冰/雪 识别(公式(6))。类似的，对于水体像元的识别则利用第1、2波段的表观反 射率(公式(7))：

4)查找表构建。为了实现不同传感器气溶胶反演，利用MODTRAN 5大气 辐射传输程序构建了一个通用的大气参数数据集(General Atmospheric Parameters Dataset,GAtmParas)，该数据集中记录有不同观测几何(18个太阳天 顶角、15个观测天顶角、19个相对方位角)及AOD条件下(0.05-3.0范围内 16个气溶胶光学厚度)0.3-2.5μm光谱范围的大气参数信息，包括t_dd(θ_s)、t_dd(θ_v)、 t_dh(θ_s)、t_hd(θ_v)、及ρ₀六个大气参数。根据MODIS光谱响应函数，利用该大气 参数数据集，通过积分运算即可完成MODIS传感器的气溶胶反演查找表的构建。

5)地表反射率的确定。在实际计算中，可利用半经验核驱动模型(公式(8)) 计算公式(1)中的r_dd、r_dh、r_hd、r_hh四个反射率。

其中，R为地表反射率值，K_vol及K_geo分别为体散射核及几何光学核，f_iso、 f_vol及f_geo分别表示各向同性散射、体散射及几何光学散射核核系数值，与波段λ 有关。

并且，实验证明，在利用RossThick体散射核及LiSparseR几何光学核积分 计算黑空及白空反照率时用式(9)能够很好的拟合。因此，本发明也采用该方 式确定各地表反射率，多项式系数如表1所示：

表1公式(9)中多项式h_k(θ)＝g_0k+g_1kθ²+g_2kθ³系数取值

6)气溶胶产品的生产。首先，根据遥感数据的成像几何条件，通过步骤(4) 构建的查找表查算出对应的大气参数；然后，利用步骤(5)确定的地表反射率 信息，结合公式(1)即可计算获取各气溶胶光学厚度条件下的表观反射率向量； 最后，利用步骤(1)计算的真实表观反射率观测值，求算出对应的气溶胶光学 厚度。本发明生产的气溶胶光学厚度产品示例如图3所示。

为了评估反演方法的性能，交叉对比分析了本发明方法获取的气溶胶光学 厚度与MODIS各气溶胶产品在北京地区的空间分布情况，如图4所示。可以看 出，各算法获取的AOD在空间分布上具有较好的一致性，均表现为西南城市地 区具有较高的AOD载荷，植被地区AOD值较低。本文提出的新算法在植被低 反射率及城市等高反射率地区均能得到AOD结果，空间连续性上较好，且空间 分辨率达到500m，可以提供比MODIS各气溶胶产品更多的细节信息。

此外，为了更加细致地对本发明反演算法精度进行评价，利用2013-2016 年香河(XHS)及北京(BJS)两个AERONET站点地面观测资料，将本发明新 算法获取的气溶胶结果分别与MODIS C6.1 DT、DB、DTDB算法产品及C6 MAIAC算法产品进行交叉对比分析，如图5和6所示。新算法获取气溶胶与 MAIAC在城市及植被均具有较高的精度，在BJS/XHS站点的在期望误差内的 比重(Percentage Within the Expected error,PWE)分别为74.9/82.1％、73.1/81.2％， 能够有效地改善MODIS DT、DB及DTDB算法对AOD的高估，提高AOD反 演精度，尤其是植被地区表现更为明显，高于期望误差内的比重(Percentage Above the Expectederror,PAE)降低了15％左右。与MAIAC算法相比，新算法 与实测数据具有更高的相关系数及更小的RMSE值。此外，由于新算法获取的 AOD空间分辨率为500m，可以提供比MODIS各产品更加详细的细节信息，这 对城市地区的环境监测等更加有利。

综合来看，利用本方法可以获取城市上空高精度、高空间分辨率的气溶胶 光学厚度信息，显示了在支持区域污染精细管控、重点城市污染传输通道监测 及污染物溯源等领域的应用潜力。

以上所述是本发明的优选实施方式，应当指出，对于本技术领域的普通技 术人员来说，在不脱离本发明所述原理的前提下，还可以做出若干改进和润饰， 这些改进和润饰应视为本发明的保护范围。