CN116519557B - 一种气溶胶光学厚度反演方法 - Google Patents
一种气溶胶光学厚度反演方法 Download PDFInfo
- Publication number
- CN116519557B CN116519557B CN202310818578.7A CN202310818578A CN116519557B CN 116519557 B CN116519557 B CN 116519557B CN 202310818578 A CN202310818578 A CN 202310818578A CN 116519557 B CN116519557 B CN 116519557B
- Authority
- CN
- China
- Prior art keywords
- image
- red
- blue
- wave band
- earth surface
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 239000000443 aerosol Substances 0.000 title claims abstract description 31
- 238000000034 method Methods 0.000 title claims abstract description 29
- 230000003287 optical effect Effects 0.000 title claims abstract description 17
- 238000002310 reflectometry Methods 0.000 claims abstract description 111
- 230000005855 radiation Effects 0.000 claims abstract description 13
- 238000007781 pre-processing Methods 0.000 claims abstract description 9
- 238000004364 calculation method Methods 0.000 claims description 12
- 238000001514 detection method Methods 0.000 claims description 12
- 238000001308 synthesis method Methods 0.000 claims description 10
- 230000015572 biosynthetic process Effects 0.000 claims description 8
- 238000003786 synthesis reaction Methods 0.000 claims description 8
- 230000005540 biological transmission Effects 0.000 claims description 6
- RGNPBRKPHBKNKX-UHFFFAOYSA-N hexaflumuron Chemical compound C1=C(Cl)C(OC(F)(F)C(F)F)=C(Cl)C=C1NC(=O)NC(=O)C1=C(F)C=CC=C1F RGNPBRKPHBKNKX-UHFFFAOYSA-N 0.000 claims description 6
- 238000012795 verification Methods 0.000 claims description 5
- 238000012937 correction Methods 0.000 claims description 4
- 238000011156 evaluation Methods 0.000 claims description 3
- 230000002194 synthesizing effect Effects 0.000 claims description 3
- 239000002131 composite material Substances 0.000 claims 1
- 239000004973 liquid crystal related substance Substances 0.000 description 14
- 238000012544 monitoring process Methods 0.000 description 5
- 238000002834 transmittance Methods 0.000 description 5
- 238000005516 engineering process Methods 0.000 description 4
- 230000003595 spectral effect Effects 0.000 description 3
- 239000005427 atmospheric aerosol Substances 0.000 description 2
- 238000010276 construction Methods 0.000 description 2
- 239000002245 particle Substances 0.000 description 2
- 230000010287 polarization Effects 0.000 description 2
- 238000012892 rational function Methods 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- PEDCQBHIVMGVHV-UHFFFAOYSA-N Glycerine Chemical compound OCC(O)CO PEDCQBHIVMGVHV-UHFFFAOYSA-N 0.000 description 1
- 102000006463 Talin Human genes 0.000 description 1
- 108010083809 Talin Proteins 0.000 description 1
- 230000004913 activation Effects 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 239000000356 contaminant Substances 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 230000008030 elimination Effects 0.000 description 1
- 238000003379 elimination reaction Methods 0.000 description 1
- 239000010419 fine particle Substances 0.000 description 1
- 238000003908 quality control method Methods 0.000 description 1
- 210000001525 retina Anatomy 0.000 description 1
- 230000001932 seasonal effect Effects 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 230000001131 transforming effect Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N15/00—Investigating characteristics of particles; Investigating permeability, pore-volume or surface-area of porous materials
- G01N15/06—Investigating concentration of particle suspensions
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
- G01N21/17—Systems in which incident light is modified in accordance with the properties of the material investigated
- G01N21/25—Colour; Spectral properties, i.e. comparison of effect of material on the light at two or more different wavelengths or wavelength bands
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V10/00—Arrangements for image or video recognition or understanding
- G06V10/10—Image acquisition
- G06V10/12—Details of acquisition arrangements; Constructional details thereof
- G06V10/14—Optical characteristics of the device performing the acquisition or on the illumination arrangements
- G06V10/143—Sensing or illuminating at different wavelengths
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V10/00—Arrangements for image or video recognition or understanding
- G06V10/70—Arrangements for image or video recognition or understanding using pattern recognition or machine learning
- G06V10/764—Arrangements for image or video recognition or understanding using pattern recognition or machine learning using classification, e.g. of video objects
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V20/00—Scenes; Scene-specific elements
- G06V20/10—Terrestrial scenes
- G06V20/13—Satellite images
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V20/00—Scenes; Scene-specific elements
- G06V20/10—Terrestrial scenes
- G06V20/188—Vegetation
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N15/00—Investigating characteristics of particles; Investigating permeability, pore-volume or surface-area of porous materials
- G01N15/06—Investigating concentration of particle suspensions
- G01N15/075—Investigating concentration of particle suspensions by optical means
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
- G01N21/17—Systems in which incident light is modified in accordance with the properties of the material investigated
- G01N2021/1793—Remote sensing
- G01N2021/1795—Atmospheric mapping of gases
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- General Health & Medical Sciences (AREA)
- Health & Medical Sciences (AREA)
- Multimedia (AREA)
- Chemical & Material Sciences (AREA)
- Analytical Chemistry (AREA)
- Biochemistry (AREA)
- Life Sciences & Earth Sciences (AREA)
- Immunology (AREA)
- Pathology (AREA)
- Evolutionary Computation (AREA)
- Astronomy & Astrophysics (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Remote Sensing (AREA)
- Artificial Intelligence (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Computing Systems (AREA)
- Databases & Information Systems (AREA)
- Dispersion Chemistry (AREA)
- Medical Informatics (AREA)
- Software Systems (AREA)
- Computer Hardware Design (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- Image Processing (AREA)
Abstract
本发明提供一种气溶胶光学厚度反演方法，包括：S1获取多光谱影像和对应的MODIS影像，并进行预处理；S2对多光谱影像进行辐射定标，得到各波段表观反射率；S3剔除云像元及云阴影像元，得到待反演影像；S4计算待反演影像中各像元的归一化植被指数，按暗像元区和亮像元区分类；S5在暗像元区利用暗目标法估算红蓝波段地表反射率比值；在亮像元区，根据MODIS影像，构建逐月红蓝波段地表反射率比值库；S6建立AOD反演查找表；S7利用AOD查找表，在不同AOD值下迭代计算红蓝波段地表反射率理论比值，与逐月红蓝波段地表反射率比值库中各数值作差，将差值最小时的AOD值作为反演结果；本发明能同时获取暗、亮像元的高分辨率AOD反演结果，准确描述城市地区气溶胶变化。
Description
技术领域
本发明涉及大气遥感技术领域，尤其涉及一种气溶胶光学厚度的反演方法。
背景技术
大气气溶胶是悬浮在空气中的分散颗粒，包括初级气溶胶（以颗粒的形式直接排放到大气中）和次级气溶胶（由大气中的主要污染物转化而来），直径为0.001~100μm。大气气溶胶的浓度变化直接影响人体健康与空气质量，并通过辐射强迫的直接或间接效应影响地气收支平衡和气候变化，因此，对气溶胶的空间分布和变化进行监测非常重要。
气溶胶光学厚度（Aerosol Optical Depth，AOD）是描述气溶胶特性的重要参数之一，可表示大气的浑浊度。传统气溶胶监测多为地基监测，然而由于站点布设有限，难以获取大范围的气溶胶空间分布信息。近年来，利用具有高时效性和大尺度观测能力的卫星遥感技术可获取空间分布连续的气溶胶数据，在宏观环境和污染分布监测上具有较大潜力。AOD反演的关键是确定地表反射率信息，以便于实现地气解耦，反演算法包括暗目标算法、深蓝算法、结构函数法、偏振算法、多角度算法等。其中，结构函数法需找到“清洁日”影像作为基准且对几何校正要求较高；偏振算法仅能应用于反演细粒子气溶胶；多角度算法需要特定的传感器支持。上述3种算法的业务化应用均较为困难。近年来暗目标算法不断得到改进与应用，但其仍具有明显局限性，不仅需要短波红外波段的支持，且不适用于亮像元地区（城市、沙漠等区域）。深蓝算法可以反演出亮像元地区的AOD但其精度低于暗目标算法，而且深蓝算法反演结果受外部地表反射率数据空间分辨率限制。
传统的气溶胶反演方法中暗、亮区域空间分辨率无法兼顾。因此，如何能在不损失原有数据空间分辨率的前提下，同时反演出暗、亮像元地区气溶胶光学厚度，是当前气溶胶光学厚度反演研究的难点。
发明内容
鉴于上述问题，本发明提供一种气溶胶光学厚度的反演方法，通过构建逐月红蓝波段地表反射率比值库采用红蓝波段比值法，同时获取包括暗、亮像元的高空间分辨率AOD反演结果，空间连续性较好，能够更准确描述城市地区气溶胶空间变化情况。
本发明的技术方案如下：
S1 获取多光谱影像和对应的MODIS影像，对多光谱影像进行预处理；
S2 对预处理后的多光谱影像进行辐射定标，得到各波段的表观反射率数据，所述波段包括蓝光波段、绿光波段、红光波段、近红外波段；
S3 对预处理后的多光谱影像进行剔除云像元及云阴影像元处理，得到待反演影像；
S4 计算待反演影像中各像元的归一化植被指数，根据预设的阈值将待反演影像中各像元按照归一化植被指数进行分类，分为暗像元区和亮像元区；
S5在暗像元区利用暗目标法估算第一红蓝波段地表反射率比值；在亮像元区，通过谷歌地球引擎平台，根据MODIS影像中各波段的地表反射率数据，利用次最小值合成法构建逐月红蓝波段地表反射率比值库；
S6 基于大气辐射传输模型对不同卫星传感器在不同条件下探测到的表观反射率进行模拟，建立AOD查找表；
S7 利用AOD查找表，在不同AOD值下迭代计算红蓝波段地表反射率理论比值，将所述理论比值分别与第一红蓝波段地表反射率比值和逐月红蓝波段地表反射率比值库中各数值作差，取差值最小时的红蓝波段地表反射率理论比值所对应的AOD查找表的AOD值作为多光谱影像的反演结果；
具体地，步骤S5中在暗像元区，利用暗目标算法估算第一红蓝波段地表反射率比值，包括以下步骤：
在待反演影像的暗像元区，通过线性关系使用短波红外波段表观反射率估算红蓝波段地表反射率，计算公式为：
其中，为蓝光波段地表反射率，/>为红光波段地表反射率，/>为2.1μm短红外波段TOA表观反射率；
将蓝光波段地表反射率除以红光波段地表反射率，得到第一红蓝波段地表反射率比值 为固定值2。
具体地，步骤S5中在亮像元区，利用次最小值合成法构建逐月红蓝波段地表反射率比值库，包括以下步骤：
通过谷歌地球引擎平台获取MODIS影像；
对MODIS影像进行预处理，按最小值合成的方式将长时间序列内同日的多景影像合成为一景影像，并裁剪出与多光谱影像时空相匹配的区域；
获取每个月份的地表反射率数据，利用次最小值合成法构建逐月红蓝波段地表反射率比值库：
其中，代表构建逐月红蓝波段地表反射率比值图像库；/>、/>代表一个月MODIS影像的红蓝波段地表反射率比值图像。
具体地，步骤S4包括：
S41 根据多光谱影像各像元的红光波段、近红外波段的表观反射率，计算各像元的归一化差值植被指数NDVI；
S42将NDVI＞0.55的像元，标记待反演影像上的对应像元为暗像元区；
S43 将NDVI≦0.55的像元，标记待反演影像上的对应像元为亮像元区。
1. 具体地，在步骤S7中在不同AOD值下迭代计算，不同AOD值的变化在0～2.0范围内。
具体地，步骤S1中所述预处理包括正射、几何精校正。
具体地，步骤S3包括：
S31 通过计算得到每幅多光谱影像各像元的波段值；
S32 利用长时间序列的遥感样本，构建含有蓝光波段的晴空背景场，根据含有蓝光波段的晴空背景场以及多光谱影像的各像元的蓝光波段反射率，计算得到云检测阈值；
S33 利用云检测阈值对待反演影像进行云检测，剔除云及云阴影像元，得到待反演影像。
具体地，步骤S8包括：使用AERONET和SONET站点数据对反演所得AOD值进行精度验证，评价指标包括均方根误差、期望误差，计算公式为：
其中，为均方根误差，/>为期望误差，/>为反演AOD值，/>为站点AOD值，n为验证样本数。
本发明的有益效果为：
（1）本发明将MODIS影像分为暗像元区域和亮像元区域，针对不同区域采取不同的计算方法，能够同时获取包括暗、亮像元的高空间分辨率AOD反演结果，空间连续性较好，能够更准确描述城市地区气溶胶空间变化情况；
（2）在亮像元区域，本发明利用次最小值合成法构建比值库，该方法假设大多数地物的地表反射率在一个月内保持不变，剔除了影像中的负值和无效值，受云、山体、建筑物及其阴影的影响较小，提高了反演结果的精度；
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明实施例的方法示意图；
图2为本发明实施例的技术流程图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员所获得的所有其他实施例，都属于本发明保护的范围。
为了实现上述目的，本发明提供一种气溶胶光学厚度反演方法，请参照图1和图2，该方法包括以下步骤：
S1 获取多光谱影像和对应的MODIS影像，对多光谱影像进行预处理；
本发明的实施例中，多光谱影像来自于Landsat-8卫星搭载的OLI（OperationalLand Imager）传感器的2022全年度粤港澳大湾区的的影像数据，利用卫星图像自带的有理多项式系数(Rational PolynomialCoefficients，RPC)，通过有理函数模型(RationalFunction Model，RFM)对影像进行正射、几何校正。MODIS影像为MOD09A1数据，数据是经大气校正后的地表反射率数据，包含红、蓝波段的信息。
S2 对预处理后的多光谱影像进行辐射定标，得到各波段的表观反射率数据，所述波段包括蓝光波段、绿光波段、红光波段、近红外波段；
本发明的实施例中，首先获取定标信息，包括最大最小辐射亮度值、像素的最大最小DN值、传感器的状态参数，以及太阳几何信息，包括太阳天顶角、太阳方位角。对卫星传感器探测信号的辐射定标，将图像的DN值转换为各波段的表观反射率，具体为根据定标公式和太阳几何信息计算表观反射率，定标公式为：
其中：
：传感器获取的光谱辐射亮度，单位是W/(m2sr.μm)；
：像素对应的DN值；
和/>：分别为像素最大和最小DN值；
和/>：分别为对应像素最大和最小DN值的光谱辐射亮度，单位是W/(m2sr.μm)。
光谱辐射亮度到表观反射率的转换公式为：
其中：
：表观反射率，无单位；
：圆周率，大约等于3.14159，无单位；
：日地距离，天文单位；
：平均天顶太阳辐照度，单位W/(m2.μm)；
：太阳天顶角，单位度。
S3 对预处理后的多光谱影像进行剔除云像元及云阴影像元处理，得到待反演影像；
本发明的实施例中，步骤S3包括：
S31 通过计算得到每幅多光谱影像各像元的波段值；
S32 利用长时间序列的遥感样本，构建含有蓝光波段的晴空背景场，根据含有蓝光波段的晴空背景场以及多光谱影像的各像元的蓝光波段反射率，计算得到云检测阈值；
S33 利用云检测阈值对多光谱影像进行云检测，剔除云及云阴影像元，得到待反演影像。
本发明实施例中，首先获取一组历史遥感样本，该遥感样本为长时间序列，对遥感样本中的影像进行云像元识别，并基于泛洪算法对云阴影像元进行识别，对云及云阴影像元进行剔除，筛选出晴空像元，得到T幅晴空影像，对这T幅晴空影像的相同位置的像元，获取蓝光波段反射率的最小值，将所有的最小值按照对应的像元排列顺序进行排列，得到含有蓝光波段的晴空背景场。
基于元数据计算得到每个遥感数据的各个像元(x,y)的波段值，包括：蓝光波段反射率Band1(x,y)、绿光波段反射率Band2(x,y)、红光波段反射率Band3(x,y)、近红外波段反射率Band4(x,y)和第十五波段的亮温Band15(x,y)。
获取晴空背景场各个像元的蓝光波段的地表信息Band1clear-sky(x,y)，将每个遥感数据的像元(x,y)的蓝光波段反射率与Band1clear-sky(x,y)作差，得到差值，对差值进行带色彩活肤的多尺度视网膜增强，得到反射率差值，通过大津算法从反射率差值图像中提取云检测阈值。
根据云检测阈值判断遥感数据中每个像元是否为云层像元，对云进行识别，同时用泛洪算法识别云阴影，对遥感数据的云及云阴影进行识别并剔除。
S4 计算待反演影像中各像元的归一化差值植被指数NDVI，根据预设的阈值将待反演影像中各像元按照归一化差值植被指数进行分类，分为暗像元区和亮像元区；
本发明实施例中，步骤S4包括：
S41 根据多光谱影像各像元的红光波段、近红外波段的表观反射率，计算各像元的归一化植被指数：
其中，为红光波段的表观反射率，/>为近红外波段的表观反射率，NDVI为所述归一化植被指数；
S42将NDVI＞0.55的像元，标记待反演影像上的对应像元为暗像元区；
S43 将NDVI≦0.55的像元，标记待反演影像上的对应像元为亮像元区。
本发明按照预设NDVI阈值，对待反演影像做划分，其中第一分类集中的数据的NDVI值均＞0.55，该分类集中的数据包含大部分的植被区域，且植被浓密，即这些数据中的场景可看成暗像元区域。第二分类集中的数据的NDVI≦0.55，该分类集中的数据包含大部分的复杂地表，即这些数据中的场景可看成亮像元区域。
S5在暗像元区利用暗目标法估算第一红蓝波段地表反射率比值；在亮像元区，通过谷歌地球引擎平台，根据MODIS影像中各波段的地表反射率数据，利用次最小值合成法构建逐月红蓝波段地表反射率比值库；
本发明的实施例中，步骤S5中在暗像元区，利用暗目标算法估算红蓝波段地表反射率比值，包括以下步骤：
在待反演影像的暗像元区，通过线性关系使用短波红外波段表观反射率估算红蓝波段地表反射率，计算公式为：
其中，为蓝光波段地表反射率，/>为红光波段地表反射率，/>为2.1μm短红外波段TOA表观反射率；
将蓝光波段地表反射率除以红光波段地表反射率，得到第一红蓝波段地表反射率比值 为固定值2。
在亮像元区域诸如城市、裸土等，暗目标算法所使用的红、蓝波段和短波红外波段之间通常不满足线性关系，所以不可以使用暗目标算法对非暗像元地表进行气溶胶光学厚度反演。亮像元区域，地表反射率随时间的变化并不如暗像元区具有明显的季节和年际变化。因此可以基于亮像元区域在一定时期内地表反射率变化较小的假设，采用最小值合成技术构建红蓝波段的月尺度地表反射率比值库。由于建成区高大建筑等的影响，会造成地表反射率最小值不正确的情况出现，因此次最小值合成（第二最小值合成）法被用于地表反射率比值库的构建。
谷歌地球引擎（Google Earth Engine, GEE）是一个基于云的大尺度地理空间分析和自然环境监测平台，它可利用 Google 的海量计算能力来处理PB级别的遥感影像数据。GEE对众多遥感影像数据进行了存档，包括 MODIS、Landsat、Sentinel 等。
本发明的实施例中，步骤S5中在亮像元区，利用次最小值合成法构建逐月红蓝波段地表反射率比值库，包括以下步骤：
通过谷歌地球引擎平台获取MODIS影像；
对MODIS影像进行预处理，按最小值合成的方式将长时间序列内同日的多景影像合成为一景影像，并裁剪出与多光谱影像时空相匹配的区域；
获取每个月份的地表反射率数据，利用次最小值合成法构建逐月红蓝波段地表反射率比值库：
其中，代表构建逐月红蓝波段地表反射率比值图像库；/>、/>代表一个月MODIS影像的红蓝波段地表反射率比值图像。
本发明使用的是GEE存档的MOD09A1经过大气校正的地表反射率数据。首先通过GEE ImageCollection获取粤港澳大湾区每个月的地表反射率数据集，然后使用map函数通过影像QA质量控制文件对数据集进行云及云阴影像元剔除除操作，接着采用次最小合成技术构建逐月红蓝波段地表反射率比值库，最后导出到谷歌云盘中进行下载，用于亮像元区域的 AOD 反演。
由于MOD09A1重访周期较长，又经常受到云、云阴影等的影响，导致采用最小值合成技术会导致地表反射率库的较多缺失值出现，因此我们采用了基于GEE 的次最小合成技术用于红蓝波段地表反射率库比值构建，并将分辨率500m降尺度到4km，实现原理如下：
其中，代表构建逐月红蓝波段地表反射率比值图像库；/>、/>代表每个月粤港澳大湾区MOD09A1影像红蓝波段地表反射率比值图像。
S6 基于大气辐射传输模型对不同卫星传感器在不同条件下探测到的表观反射率进行模拟，建立AOD查找表；
本发明点实施例中，利用6S大气辐射传输模型，建立AOD反演查找表，其参数设置如下：
1个卫星天顶角：多光谱影像元数据获取，单位为度；
1个太阳天顶角：多光谱影像元数据获取，单位为度；
1个相对方位角：多光谱影像元数据获取，单位为度；
17个AOD：0.001, 0.01, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8,0.9, 1.0, 1.15, 1.4, 1.65, 2.0；
2个大气类型：中纬度夏季和中纬度冬季；
1个气溶胶类型：大陆型气溶胶；
输出参数：大气程辐射反射率、入射光线从大气顶层到达地面的总透过率、向上进入卫星传感器视场方向的总透过率、大气底层向下的半球反射率；
针对每一景多光谱影像生成一个查找表，这景影像的角度是确定的，因此可以多光谱元数据中获取。
S7 利用AOD查找表，在不同AOD值下迭代计算红蓝波段地表反射率理论比值，将所述理论比值分别与第一红蓝波段地表反射率比值和逐月红蓝波段地表反射率比值库中各数值作差，取差值最小时的红蓝波段地表反射率理论比值所对应的AOD查找表的AOD值作为多光谱影像的反演结果；
当假设陆地表面为均匀的朗伯体表面时，卫星传感器所测量的大气顶层的表观反射率可以表示为：
其中：
θs：太阳天顶角；
θv：卫星天顶角；
φ：相对方位角；
：表观反射率；
：大气程辐射反射率；
：总气体分子透过率；
：入射光线从大气顶层到达地面的总透过率；
：向上进入卫星传感器视场方向的总透过率；
S：大气底层向下的半球反射率，可由6S辐射传输模型计算得到；
：地表反射率。
在本发明实施例中，步骤S7中，在不同AOD值下迭代计算，不同AOD值的变化在0～2.0范围内，步长0.01。
根据上述计算的表观反射率、以及6S大气辐射传输模型输出的三个参数，以0.001作为初始AOD值进行计算，步长设置为0.01，通过下式迭代计算不同AOD条件下的红、蓝波段地表反射率，将上述公式变换可得红、蓝波段地表反射率的求解公式：
根据上述得到的红、蓝波段理论地表反射率，计算其比值，分别与红蓝波段地表反射率比值和逐月红蓝波段地表反射率比值库中各数值作差，取差值绝对值最小时当前迭代的AOD值作为多光谱影像的反演结果；
验证数据来自于AERONET和SONET，SONET是由中国科学院联合中国多所高等院校和科研所等单位在中国典型区域建立的观测网，SONET站点在中国境内分布更均匀，覆盖全国多个省域，目前共有23个站点，本发明的实施例中，选用粤港澳大湾区的Level 1.5级观测数据，用于验证AOD反演结果的精度。
S8利用站点地面观测数据验证反演结果的精度。
本发明的实施例中，使用AERONET和SONET站点数据对反演所得AOD值进行精度验证，评价指标包括均方根误差（RMSE）、期望误差（EE），计算公式如下：
以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以权利要求的保护范围为准。
Claims (6)
1.一种气溶胶光学厚度反演方法，其特征在于，包括以下步骤：
S1获取多光谱影像和对应的MODIS影像，对多光谱影像进行预处理；
S2对预处理后的多光谱影像进行辐射定标，得到各波段的表观反射率数据，所述波段包括蓝光波段、绿光波段、红光波段、近红外波段；
S3对预处理后的多光谱影像进行剔除云像元及云阴影像元处理，得到待反演影像；
S4计算待反演影像中各像元的归一化差值植被指数，根据预设的阈值将待反演影像中各像元按照归一化差值植被指数进行分类，分为暗像元区和亮像元区；
S5在暗像元区利用暗目标法估算第一红蓝波段地表反射率比值；在亮像元区，通过谷歌地球引擎平台，根据MODIS影像中各波段的地表反射率数据，利用次最小值合成法构建逐月红蓝波段地表反射率比值库；
S6基于大气辐射传输模型对不同卫星传感器在不同条件下探测到的表观反射率进行模拟，建立AOD查找表；
S7利用AOD查找表，在不同AOD值下迭代计算红蓝波段地表反射率理论比值，将所述理论比值分别与第一红蓝波段地表反射率比值和逐月红蓝波段地表反射率比值库中各数值作差，取差值最小时的红蓝波段地表反射率理论比值所对应的AOD查找表的AOD值作为多光谱影像的反演结果；
S8利用站点地面观测数据验证反演结果的精度；
其中，步骤S5中在暗像元区，利用暗目标算法估算第一红蓝波段地表反射率比值，包括以下步骤：
在待反演影像的暗像元区，通过线性关系使用短波红外波段表观反射率估算红蓝波段地表反射率，计算公式为：
其中，ρBLUE为蓝光波段地表反射率，ρRED为红光波段地表反射率，ρSWIR2为2.1μm短红外波段TOA表观反射率；
将红光波段地表反射率除以蓝光波段地表反射率，得到第一红蓝波段地表反射率比值为固定值2；
步骤S5中在亮像元区，利用次最小值合成法构建逐月红蓝波段地表反射率比值库，包括以下步骤：
通过谷歌地球引擎平台获取MODIS影像；
对MODIS影像进行预处理，按最小值合成的方式将长时间序列内同日的多景影像合成为一景影像，并裁剪出与多光谱影像时空相匹配的区域；
获取每个月份的地表反射率数据，利用次最小值合成法构建逐月红蓝波段地表反射率比值图像库ρ(mon,wrs)：
ρ(mon,wrs)＝2ndMin_composite(ρ1(mon,wrs),...ρn(mon,wrs))
其中，ρ(mon,wrs)代表构建逐月红蓝波段地表反射率比值图像库；ρ1,...,ρn代表一个月MODIS影像的红蓝波段地表反射率比值图像。
2.根据权利要求1所述的一种气溶胶光学厚度反演方法，其特征在于，步骤S4包括：
S41根据多光谱影像各像元的红光波段、近红外波段的表观反射率，计算各像元的归一化差值植被指数NDVI；
S42将NDVI＞0.55的像元，标记待反演影像上的对应像元为暗像元区；
S43将NDVI≦0.55的像元，标记待反演影像上的对应像元为亮像元区。
3.根据权利要求1所述的一种气溶胶光学厚度反演方法，其特征在于，步骤S1中所述预处理包括正射、几何精校正。
4.根据权利要求1所述的一种气溶胶光学厚度反演方法，其特征在于，步骤S3包括：
S31获取每幅多光谱影像各像元的波段值；
S32利用长时间序列的遥感样本，构建含有蓝光波段的晴空背景场，根据含有蓝光波段的晴空背景场以及多光谱影像的各像元的蓝光波段反射率，计算得到云检测阈值；
S33利用云检测阈值对待反演影像进行云检测，剔除云及云阴影像元，得到待反演影像。
5.根据权利要求1所述的一种气溶胶光学厚度反演方法，其特征在于，在步骤S7中在不同AOD值下迭代计算，不同AOD值的变化在0～2.0范围内。
6.根据权利要求1所述的一种气溶胶光学厚度反演方法，其特征在于，步骤S8包括：使用AERONET和SONET站点数据对反演所得AOD值进行精度验证，评价指标包括均方根误差、期望误差，计算公式为：
EE＝±(0.05+0.2×τ0)
其中，RMSE为均方根误差，EE为期望误差，τ为反演AOD值，τ0为站点AOD值，n为验证样本数。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310818578.7A CN116519557B (zh) | 2023-07-05 | 2023-07-05 | 一种气溶胶光学厚度反演方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310818578.7A CN116519557B (zh) | 2023-07-05 | 2023-07-05 | 一种气溶胶光学厚度反演方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116519557A CN116519557A (zh) | 2023-08-01 |
| CN116519557B true CN116519557B (zh) | 2023-09-19 |
Family
ID=87398025
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310818578.7A Active CN116519557B (zh) | 2023-07-05 | 2023-07-05 | 一种气溶胶光学厚度反演方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116519557B (zh) |
Families Citing this family (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117347282B (zh) * | 2023-08-22 | 2024-05-28 | 中南大学 | 星基气溶胶光学厚度反演方法、装置及系统和存储介质 |
| CN117607919B (zh) * | 2023-11-17 | 2024-06-21 | 中国科学院大气物理研究所 | 一种基于城市建筑物阴影的气溶胶卫星遥感反演方法 |
| CN117555047A (zh) * | 2023-12-07 | 2024-02-13 | 中国人民解放军国防科技大学 | 基于人工智能的电力设备气象监测预警方法 |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106407656A (zh) * | 2016-08-29 | 2017-02-15 | 中国科学院遥感与数字地球研究所 | 一种基于高分辨率卫星影像数据的气溶胶光学厚度反演方法 |
| CN109030301A (zh) * | 2018-06-05 | 2018-12-18 | 中南林业科技大学 | 基于遥感数据的气溶胶光学厚度反演方法 |
| CN110186823A (zh) * | 2019-06-26 | 2019-08-30 | 中国科学院遥感与数字地球研究所 | 一种气溶胶光学厚度反演方法 |
| CN114113001A (zh) * | 2022-01-27 | 2022-03-01 | 航天宏图信息技术股份有限公司 | 一种气溶胶光学厚度反演方法 |
| CN116008140A (zh) * | 2022-09-28 | 2023-04-25 | 武汉大学 | 基于多波段卫星数据的气溶胶光学厚度反演方法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110261873A (zh) * | 2019-05-29 | 2019-09-20 | 电子科技大学 | 一种基于分段统计的大气气溶胶反演方法 |
-
2023
- 2023-07-05 CN CN202310818578.7A patent/CN116519557B/zh active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106407656A (zh) * | 2016-08-29 | 2017-02-15 | 中国科学院遥感与数字地球研究所 | 一种基于高分辨率卫星影像数据的气溶胶光学厚度反演方法 |
| CN109030301A (zh) * | 2018-06-05 | 2018-12-18 | 中南林业科技大学 | 基于遥感数据的气溶胶光学厚度反演方法 |
| CN110186823A (zh) * | 2019-06-26 | 2019-08-30 | 中国科学院遥感与数字地球研究所 | 一种气溶胶光学厚度反演方法 |
| CN114113001A (zh) * | 2022-01-27 | 2022-03-01 | 航天宏图信息技术股份有限公司 | 一种气溶胶光学厚度反演方法 |
| CN116008140A (zh) * | 2022-09-28 | 2023-04-25 | 武汉大学 | 基于多波段卫星数据的气溶胶光学厚度反演方法 |
Non-Patent Citations (3)
| Title |
|---|
| GF-4增强型地表反射率库支持法的气溶胶光学厚度反演;马小雨;陈正华;宿鑫;于会泳;贾丹丹;姚焕玫;;遥感学报(第05期);全文 * |
| 基于Himawari-8的气溶胶反演研究;韦海宁;王维真;黄广辉;徐菲楠;冯姣姣;董磊磊;;遥感技术与应用(第04期);全文 * |
| 资源一号04星WFI数据气溶胶反演与验证――以北京市及北京周边地区为例;丁宇;汪小钦;王峰;;遥感信息(第03期);全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116519557A (zh) | 2023-08-01 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109581372B (zh) | 一种生态环境遥感监测方法 | |
| CN116519557B (zh) | 一种气溶胶光学厚度反演方法 | |
| De Keukelaere et al. | Atmospheric correction of Landsat-8/OLI and Sentinel-2/MSI data using iCOR algorithm: validation for coastal and inland waters | |
| Webster et al. | Three-dimensional thermal characterization of forest canopies using UAV photogrammetry | |
| Li et al. | Evaluation of Sentinel-2A surface reflectance derived using Sen2Cor in North America | |
| Eaton et al. | Reflected irradiance indicatrices of natural surfaces and their effect on albedo | |
| Goslee | Analyzing remote sensing data in R: the landsat package | |
| CN109974665B (zh) | 一种针对缺少短波红外数据的气溶胶遥感反演方法及系统 | |
| Mohamed et al. | Land surface temperature and emissivity estimation for Urban Heat Island assessment using medium-and low-resolution space-borne sensors: A review | |
| CN102854513B (zh) | 环境一号hj-1a/b星ccd数据的云检测方法 | |
| Mitchell et al. | Use of hemispherical piiotographs in forest ecology | |
| Rautiainen et al. | Coupling forest canopy and understory reflectance in the Arctic latitudes of Finland | |
| CN114113001B (zh) | 一种气溶胶光学厚度反演方法 | |
| CN111832518A (zh) | 基于时空融合的tsa遥感影像土地利用方法 | |
| CN109406361B (zh) | 一种基于遥感技术的干旱区灰霾污染预警方法 | |
| CN114970214A (zh) | 一种气溶胶光学厚度反演方法及装置 | |
| Mamaghani et al. | An initial exploration of vicarious and in-scene calibration techniques for small unmanned aircraft systems | |
| CN116822141A (zh) | 利用卫星微光遥感反演夜间大气气溶胶光学厚度的方法 | |
| Zhang et al. | Development of the direct-estimation albedo algorithm for snow-free Landsat TM albedo retrievals using field flux measurements | |
| He et al. | Direct estimation of land surface albedo from simultaneous MISR data | |
| Zheng et al. | Ice/snow surface temperature retrieval from chinese FY-3d MERSI-II data: algorithm and preliminary validation | |
| CN116185616A (zh) | 一种fy-3d mersi l1b数据自动化再处理方法 | |
| CN115452167A (zh) | 基于不变像元的卫星遥感器交叉定标方法和装置 | |
| Sartika et al. | Determining the Precision of Spectral Patterns Arising from Atmospheric Correction Utilizing MODTRAN-FLAASH and 6S Approaches on High-Resolution SPOT-6 Imagery | |
| CN118090636A (zh) | 一种地基可见光高光谱成像仪对月观测数据处理方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |