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

本发明专利技术涉及一种考虑地表非朗伯特性的气溶胶光学厚度遥感反演方法，包括获取遥感数据进行预处理；利用遥感图像地表反射率数据基于先验知识约束完成半经验核驱动模型的每日核系数估算；通过离散余弦变换的惩罚最小二乘估计时间序列平滑滤波算法构建单位时间内核系数地表先验知识数据集；建立气溶胶反演通用查找表；利用半经验核驱动模型结合核系数地表先验知识数据集，确定非朗伯体辐射传输前向模型的地表反射率参数；利用上述结果反演获取气溶胶光学厚度。本发明专利技术可获取高精度的空间分辨率为500m的气溶胶光学厚度产品，能提供空间相对更为连续的气溶胶分布，显示了在支持区域污染精细管控、重点城市污染传输通道监测及污染物溯源等领域的应用潜力。

【技术实现步骤摘要】
一种考虑地表非朗伯特性的气溶胶光学厚度遥感反演方法
本专利技术涉及大气遥感领域，尤指一种考虑地表非朗伯特性的气溶胶光学厚度遥感反演方法。
技术介绍
气溶胶作为地球气候系统的重要组成部分，是大气中三大活跃成分之一，其数量与性质在时间及空间尺度上的变化均是巨大的，导致难以用固定模式予以描述。气溶胶光学厚度(AerosolOpticalDepth,AOD)作为表征气溶胶粒子消光作用的重要光学参数之一，在地气系统辐射收支平衡、全球气候变化及大气环境和人类健康等方面具有重要指示作用(张莹等，2013；王中挺等，2016)。由于大气气溶胶具有较大的时空变异性，卫星遥感成为能够提供可靠气溶胶分布信息的有效手段。卫星传感器接收到的大气层顶辐射信息是由地表及大气共同作用、相互耦合的最终结果，并且大气贡献相比地表是十分有限的，尤其是城市等高反射率区域。遥感定量化过程就是从这单一信息中，区分两个未知量，其本质是一个病态反演问题。在实际的气溶胶遥感中，需要通过先验知识来确定气溶胶类型及地表反射率参数，这两者的不确定性也正是气溶胶反演误差的主要来源。气溶胶反演算法的发展与改进，也正是围绕气溶胶类型的假设及地表反射率估算这两个方面。近年来，国内外诸多学者在气溶胶反演中对地表方向反射的估算方面做了许多努力。Thomas等提出了ORAC(Oxford-RALAerosolandCloud)算法实现了ATSR-2和AATSR数据的气溶胶参数反演，该算法是使用BRDF/Albedo产品构建反射率先验知识库实现地气分离。Xue等通过中分辨率成像光谱仪(ModerateResolutionImagingSpectroradiometer,MODIS)BRDF/反照率产品数据来增加方程约束条件，基于AATSR星下点及前向两个角度观测数据提出了考虑BRDF效应的气溶胶协同反演算法。在此基础上，Guang等利用MODIS数据提出了AOD与地表反射率双星协同反演算法，进一步提高了亮地表上空气溶胶反演精度。Qin等在假设地表BRDF形状稳定的基础上，基于比值法实现了AATSR数据澳大利亚地区气溶胶及BRDF的反演。上述诸多算法均需要大气条件不变的时间序列同步观测数据，导致其普适性较差，且使用的仍是朗伯体前向辐射传输模型计算表观反射率，未从根本上解决地表非朗伯反射的影响。
技术实现思路
为解决传统气溶胶算法中地表朗伯反射假定对地表反射率估算带来的不确定性，本算法专利技术提出一个耦合地表二向反射特性的气溶胶反演算法，最终构建一个具有较强普适性、可业务化运行的气溶胶反演算法。为实现上述目的，本专利技术所采用的技术方案是：一种考虑地表非朗伯特性的气溶胶光学厚度遥感反演方法，包括：获取遥感数据并进行预处理；地表先验知识数据集的构建：利用长时间序列的遥感图像地表反射率数据，基于先验知识条件约束完成半经验核驱动模型的每日核系数估算；通过离散余弦变换的惩罚最小二乘估计DCT-PLS时间序列平滑滤波算法构建单位时间内的核系数地表先验知识数据集；查找表构建：构建通用大气参数数据集，再根据MODIS光谱响应函数、对该大气参数数据集中的数据通过积分运算建立MODIS传感器气溶胶反演查找表文件；地表反射率参数的确定：利用RTLS半经验核驱动模型，结合构建的核系数地表先验知识数据集SANIFD，确定非朗伯体辐射传输前向模型中四个地表反射率参数；气溶胶光学厚度反演：根据成像角度信息、地表反射率参数及查找表中大气参数信息，结合预处理获取的表观反射率数据，基于非朗伯体辐射传输前项模型，反演获取气溶胶光学厚度。所述基于先验知识条件约束完成半经验核驱动模型核系数的估算，采用的代价函数为：J(f)＝(Kf-ρ)TΣT(Kf-ρ)+(f-f0)TMT(f-f0)；其中，K为三个核矩阵，f为对应的三个核系数向量，ρ为方向反射率向量，Σ为反射率误差的协方差矩阵，f0及M为核系数先验知识的均值向量及协方差矩阵；上述代价函数的最小二乘解表达式可写成：f＝(KTΣ-1K+M-1)-1(KTΣ-1ρ+M-1f0)，通过该式即可完成核系数的估算。所述半经验核驱动模型的核系数包括：各向同性散射fiso、、体散射fvol及几何光学散射fgeo三个核系数。所述单位时间内为每若干天内。通过DCT-PLS时间序列平滑滤波算法构建核系数地表先验知识数据集，具体实现公式为：式中，||||表示欧几里得范数；表示待平滑的时间序列变量；W表示最佳平滑滤波权重；s为控制滤波平滑度参数，随着s值的增大，会越平滑；对于一维数据的计算，惩罚项P通常可用Laplace差分计算。所述通用大气参数数据集记录如下数据：不同观测几何数据：18个太阳天顶角、15个观测天顶角、19个相对方位角；气溶胶光学厚度AOD条件(0.05-3.0范围内)下16个气溶胶光学厚度；0.3-2.5μm光谱范围的大气参数信息，包括tdd(θs)、tdd(θv)、tdh(θs)、thd(θv)、及ρ0。所述非朗伯体辐射传输前向模型中四个地表反射率参数包括：rdd、rdh、rhd、rhh，计算公式如下：其中，Kvol及Kgeo分别为体散射核及几何光学核，fiso、fvol及fgeo分别表示各向同性散射、体散射及几何光学散射核核系数值，与波段λ有关；R为地表反射率值，当分别计算时，rdd、rdh、rhd、rhh表示地表不同的反射特性，分别为地表方向-方向反射率、方向-半球反射率、半球-方向反射率及半球-半球反射率。非朗伯体辐射传输前项模型为：其中，矩阵T(θs)、T(θv)及分别定义为：T(θs)＝[tdd(θs)tdh(θs)]式中，θs、θv、分别为太阳天顶角、传感器天顶角、太阳方位角、传感器方位角及相对方位角；ρTOA为表观反射率；ρ0为大气路径散射反射率，表示大气程辐射项；为大气球形反照率；t和r分别表示透过率和反射率，下标h和d分别表示半球(漫射)和方向(直射)；因此，thd(θv)表征在卫星观测方向上大气对地表漫射辐射的消光，即大气上行漫射透过率；tdd(θv)表征在卫星观测方向上大气对地表反射直射辐射的消光，即大气上行直射透过率；同样，tdh(θs)表征在入射方向上大气对天空漫射辐射的消光，即大气下行漫射透过率；tdd(θs)则为入射方向上大气对太阳直射辐射的消光，即大气下行直射透过率；rdd、rdh、rhd、rhh表示地表不同的反射特性，分别为地表方向-方向反射率、方向-半球反射率、半球-方向反射率及半球-半球反射率。对遥感数据进行预处理包括云检测、冰/雪以及水体识别，采用改进的动态阈值法，具体表示为：1)首先使用地表先验知识数据集，计算获取各波段地表反射率值：然后，利用不同波段云检测动态阈值模型实现云像元的识别，公式表示为：2)由于采用的半经验核本文档来自技高网...

【技术保护点】
1.一种考虑地表非朗伯特性的气溶胶光学厚度遥感反演方法，其特征在于包括：/n获取遥感数据并进行预处理；/n地表先验知识数据集的构建：利用长时间序列的遥感图像地表反射率数据，基于先验知识条件约束完成半经验核驱动模型的每日核系数估算；通过离散余弦变换的惩罚最小二乘估计DCT-PLS时间序列平滑滤波算法构建单位时间内的核系数地表先验知识数据集；/n查找表构建：构建通用大气参数数据集，再根据MODIS光谱响应函数、对该大气参数数据集中的数据通过积分运算建立MODIS传感器气溶胶反演查找表文件；/n地表反射率参数的确定：利用RTLS半经验核驱动模型，结合构建的核系数地表先验知识数据集SANIFD，确定非朗伯体辐射传输前向模型中四个地表反射率参数；/n气溶胶光学厚度反演：根据成像角度信息、地表反射率参数及查找表中大气参数信息，结合预处理获取的表观反射率数据，基于非朗伯体辐射传输前项模型，反演获取气溶胶光学厚度。/n

【技术特征摘要】
1.一种考虑地表非朗伯特性的气溶胶光学厚度遥感反演方法，其特征在于包括：
获取遥感数据并进行预处理；
地表先验知识数据集的构建：利用长时间序列的遥感图像地表反射率数据，基于先验知识条件约束完成半经验核驱动模型的每日核系数估算；通过离散余弦变换的惩罚最小二乘估计DCT-PLS时间序列平滑滤波算法构建单位时间内的核系数地表先验知识数据集；
查找表构建：构建通用大气参数数据集，再根据MODIS光谱响应函数、对该大气参数数据集中的数据通过积分运算建立MODIS传感器气溶胶反演查找表文件；
地表反射率参数的确定：利用RTLS半经验核驱动模型，结合构建的核系数地表先验知识数据集SANIFD，确定非朗伯体辐射传输前向模型中四个地表反射率参数；
气溶胶光学厚度反演：根据成像角度信息、地表反射率参数及查找表中大气参数信息，结合预处理获取的表观反射率数据，基于非朗伯体辐射传输前项模型，反演获取气溶胶光学厚度。


2.如权利1要求所述的一种考虑地表非朗伯特性的气溶胶光学厚度遥感反演方法，其特征在于，所述基于先验知识条件约束完成半经验核驱动模型核系数的估算，采用的代价函数为：J(f)＝(Kf-ρ)TΣT(Kf-ρ)+(f-f0)TMT(f-f0)；
其中，K为三个核矩阵，f为对应的三个核系数向量，ρ为方向反射率向量，Σ为反射率误差的协方差矩阵，f0及M为核系数先验知识的均值向量及协方差矩阵；
上述代价函数的最小二乘解表达式可写成：f＝(KTΣ-1K+M-1)-1(KTΣ-1ρ+M-1f0)，通过该式即可完成核系数的估算。


3.如权利1要求所述的一种考虑地表非朗伯特性的气溶胶光学厚度遥感反演方法，其特征在于所述半经验核驱动模型的核系数包括：各向同性散射fiso、、体散射fvol及几何光学散射fgeo三个核系数。


4.如权利1要求所述的一种考虑地表非朗伯特性的气溶胶光学厚度遥感反演方法，其特征在于所述单位时间内为每若干天内。


5.如权利1要求所述的一种考虑地表非朗伯特性的气溶胶光学厚度遥感反演方法，其特征在于，通过DCT-PLS时间序列平滑滤波算法构建核系数地表先验知识数据集，具体实现公式为：



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


6.如权利1要求所述的一种考虑地表非朗伯特性的气溶胶光学厚度遥感反演方法，其特征在于所述通用大气参数数据集记录如下数据：
不同观测几何数据：18个太阳天顶角、15个观测天顶角、19个相对方位角；
气溶胶光学厚度AOD条件(0.05-3.0范围内)下16个气溶胶光学厚度；
0.3-2.5μm光谱...

【专利技术属性】
技术研发人员：田信鹏，高志强，
申请(专利权)人：中国科学院烟台海岸带研究所，
类型：发明
国别省市：山东;37