CN113740263B - 一种气溶胶光学厚度反演方法及大气颗粒物遥感反演方法 - Google Patents
一种气溶胶光学厚度反演方法及大气颗粒物遥感反演方法 Download PDFInfo
- Publication number
- CN113740263B CN113740263B CN202110266804.6A CN202110266804A CN113740263B CN 113740263 B CN113740263 B CN 113740263B CN 202110266804 A CN202110266804 A CN 202110266804A CN 113740263 B CN113740263 B CN 113740263B
- Authority
- CN
- China
- Prior art keywords
- aerosol
- optical thickness
- remote sensing
- pixel
- atmospheric
- 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 125
- 230000003287 optical effect Effects 0.000 title claims abstract description 89
- 238000000034 method Methods 0.000 title claims abstract description 72
- 238000012544 monitoring process Methods 0.000 claims abstract description 31
- 238000011160 research Methods 0.000 claims abstract description 22
- 238000002310 reflectometry Methods 0.000 claims abstract description 19
- 238000012937 correction Methods 0.000 claims description 25
- 230000008033 biological extinction Effects 0.000 claims description 15
- CURLTUGMZLYLDI-UHFFFAOYSA-N Carbon dioxide Chemical compound O=C=O CURLTUGMZLYLDI-UHFFFAOYSA-N 0.000 claims description 14
- 230000005540 biological transmission Effects 0.000 claims description 13
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 13
- 230000005855 radiation Effects 0.000 claims description 11
- CBENFWSGALASAD-UHFFFAOYSA-N Ozone Chemical compound [O-][O+]=O CBENFWSGALASAD-UHFFFAOYSA-N 0.000 claims description 9
- 238000010521 absorption reaction Methods 0.000 claims description 9
- 238000005259 measurement Methods 0.000 claims description 9
- 229910002092 carbon dioxide Inorganic materials 0.000 claims description 7
- 239000002245 particle Substances 0.000 claims description 7
- 239000001569 carbon dioxide Substances 0.000 claims description 6
- 238000005070 sampling Methods 0.000 claims description 6
- 230000001360 synchronised effect Effects 0.000 claims description 6
- 238000002834 transmittance Methods 0.000 claims description 6
- 239000005427 atmospheric aerosol Substances 0.000 claims description 5
- 239000004576 sand Substances 0.000 claims description 4
- 239000000779 smoke Substances 0.000 claims description 4
- 238000010998 test method Methods 0.000 claims description 3
- 230000000873 masking effect Effects 0.000 claims description 2
- 238000012545 processing Methods 0.000 claims description 2
- 239000008277 atmospheric particulate matter Substances 0.000 abstract description 8
- 239000011324 bead Substances 0.000 description 8
- 238000004364 calculation method Methods 0.000 description 4
- 239000000428 dust Substances 0.000 description 4
- 238000004088 simulation Methods 0.000 description 3
- 238000004458 analytical method Methods 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 230000007613 environmental effect Effects 0.000 description 2
- 238000011835 investigation Methods 0.000 description 2
- 239000000126 substance Substances 0.000 description 2
- 238000003915 air pollution Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000004422 calculation algorithm Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000010354 integration Effects 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 239000013618 particulate matter Substances 0.000 description 1
- 230000001932 seasonal effect Effects 0.000 description 1
- 239000002689 soil Substances 0.000 description 1
- 230000003595 spectral effect Effects 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
Images
Classifications
-
- 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
-
- 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
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
- G01S19/38—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system
- G01S19/39—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system the satellite radio beacon positioning system transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
-
- 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
-
- 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)
- Chemical & Material Sciences (AREA)
- Remote Sensing (AREA)
- Radar, Positioning & Navigation (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- General Health & Medical Sciences (AREA)
- Biochemistry (AREA)
- Analytical Chemistry (AREA)
- Immunology (AREA)
- Pathology (AREA)
- Computer Networks & Wireless Communication (AREA)
- Dispersion Chemistry (AREA)
- Investigating Or Analysing Materials By Optical Means (AREA)
- Image Processing (AREA)
- Length Measuring Devices By Optical Means (AREA)
Abstract
本申请属于空气质量监测技术领域，具体涉及一种气溶胶光学厚度反演方法及大气颗粒物遥感反演方法；气溶胶光学厚度反演方法包括如下步骤：获取研究区MODIS L1B遥感数据并进行订正；建立超像元窗口；获得2.13μm像元相关暗像元集合；获得0.66μm像元相关像元集合；求交得到有效暗像元集合；计算出波段为0.47μm和0.66μm的地表反射率；建立四维查找表；反演出各像元在对应波段的气溶胶光学厚度值；反演0.55μm波段的气溶胶光学厚度。本申请的气溶胶光学厚度反演方法，通过改进的暗像元法进行反演，能够有效减少反演空白值，提高反演精度和效率；另外，大气颗粒物遥感反演方法中能够准确获取500米分辨率的大气颗粒物地面浓度空间分布。
Description
技术领域
本申请属于空气质量监测技术领域，具体涉及一种气溶胶光学厚度反演方法及大气颗粒物遥感反演方法。
背景技术
卫星遥感技术具有广泛的空间覆盖及对地球和大气的重复观测的优势，用卫星资料反演气溶胶产品为空气质量监测提供了一种新思路。气溶胶光学厚度(Aerosol opticalthickness，AOT)是指沿辐射传输路径单位界面上所有吸收和散射物质在各种环境条件下(湿度、风速、垂直廓线等)产生的总消光系数在垂直方向上的积分，是评估气溶胶引起的大气污染和气候效应的一个重要因子。已有研究表明，AOT与空气污染指数、颗粒物实测浓度之间存在显著的相关性。颗粒物的遥感反演是利用这种相关性建立卫星反演的AOT与地面观测的颗粒物浓度之间的模型，并将模型应用到一定的时空范围内，反演出近地面颗粒物的空间分布。
目前使用遥感技术反演大气颗粒物(PM)空间分布存在以下几个问题：①气溶胶光学厚度(AOT)产品数据精度低，局部地区反演结果精度较低，且空间分辨率在10km左右甚至更低，不适合小区域尺度上的空气质量监测，应用性较弱；②订正AOT的气象数据多为站点监测值或再分析资料，空间分辨率低；③监测站点和卫星遥感数据的空间尺度存在差异；④大多AOT-PM关系模型仅针对特定地点和有限时间范围，时空适用性较弱。
发明内容
为了解决现有技术中存在的至少一个技术问题，本发明提供了一种气溶胶光学厚度反演方法及大气颗粒物遥感反演方法，其综合土地覆盖和季节因素，能较为准确地获取500米分辨率的大气颗粒物空间分布。
第一方面，本申请公开了一种气溶胶光学厚度反演方法，包括如下步骤：
步骤一、获取研究区的无云或相对少云覆盖、分辨率为500米的MODIS L1B遥感数据，并对所述MODIS L1B遥感数据进行订正处理；
步骤二、在所述MODIS L1B遥感数据中，以左上角第一个像元为被计算像元，建立分辨率为10*10km且包含N*N个像元的超像元窗口；
步骤三、在所述超像元窗口中，针对2.13μm波段像元，以表观反射率大于等于0.01且小于等于0.25为条件从中提取像元，以得到第一个暗像元集合；
步骤四、在所述超像元窗口中，针对0.66μm波段像元，去掉该波段亮度的前50％和后20％，以形成第二个像元集合；
步骤五、将所述第一个暗像元集合和所述第二个像元集合求交，得到有效暗像元集合；
步骤六、判断所述有效暗像元集合中的有效暗像元数量是否超过M个；如果超过，则进入步骤七；否则，返回步骤二，且将被计算像元的气溶胶光学厚度值设为-1；
步骤七、在所述有效暗像元集合中，计算波段为0.47μm、0.66μm、1.24μm和2.13μm的表观反射率以及NDVI值，并计算出波段为0.47μm和0.66μm的地表反射率；
步骤八、根据6S辐射传输模型、MODIS L1B遥感数据的几何参数以及研究区的气溶胶类型参数，建立四维查找表；
步骤九、将波段为0.47μm和0.66μm的地表反射率代入到所述四维查找表中，逐像元移动超像元窗口，反演出各像元在对应波段的气溶胶光学厚度值；
步骤十、在气溶胶光学厚度与波长的如下公式(4)中，先引入波段0.47μm和0.66μm的气溶胶光学厚度，以得到α和β值，再通过该公式反演出分辨率为500米、波段为0.55μm的气溶胶光学厚度：
τ＝βλ-α (4)；
其中τ为气溶胶光学厚度，λ为波长，β为Angstrom混浊系数，代表大气中气溶胶的浓度，α为Angstrom波长指数。
根据本申请的至少一个实施方式，在所述步骤一中，对所述MODIS L1B遥感数据进行的所述订正处理包括对云和水域进行掩膜处理，以及通过如下公式(1)-(3)进行水汽、臭氧、CO2订正处理：
其中，
表示表示水汽透过率；/>
表示臭氧透过率；/>
表示二氧化碳透过率；G表示大气质量因子；/>
表示水汽吸收系数；/>
表示臭氧吸收系数；/>
表示二氧化碳吸收系数。
根据本申请的至少一个实施方式，在所述步骤二中，是建立包含20*20个像元的超像元窗口，且在所述步骤六中，是判断有效暗像元数量是否超过12个。
根据本申请的至少一个实施方式，在所述步骤八中，所述6S辐射传输模型所需参数变量包括通道、高程、大气顶层反射率以及光学厚度；以及
所述MODIS L1B遥感数据的几何参数包括太阳高度角和太阳天顶角；以及
所述研究区的气溶胶类型参数包括沙尘、水溶性、海洋、烟尘。
第二方面，本申请还公开了一种大气颗粒物遥感反演方法，包括如下步骤：
步骤一、采用如第一方面中任一项所述的气溶胶光学厚度反演方法获得研究区的气溶胶光学厚度；
步骤二、对气溶胶光学厚度进行垂直和湿度订正，并根据以下公式(6)计算干空气状态下的地表气溶胶消光系数：
kAOT，DRY(λ)＝τa(λ)/[HA·f(RH)] (6)；
其中，f(RH)代表相对湿度的校正系数，f(RH)＝(1-RH/100)-1；kAOT，DRY(λ)代表干空气状态下的地表气溶胶消光系数，τa(λ)代表通过遥感影像反演的气溶胶光学厚度，RH为相对湿度的百分数，HA为大气气溶胶标高，即假定气溶胶浓度随高度分布保持不变时的气溶胶层等效厚度；
步骤三、将所述干空气状态下的地表气溶胶消光系数与大气颗粒物站点监测数据结合，建立不同类型土地覆盖和不同季节的大气颗粒物遥感反演模型(7)：
Y＝aX+b (7)；
其中，X为对应站点所在像元的订正后气溶胶光学厚度值，即kAOT，DRY(λ)；Y为PM10站点监测值；a、b为系数，是根据具体类型土地覆盖以及具体季节时的X、Y带入值求得；另外，所述大气颗粒物遥感反演模型的分辨率为500米；
步骤四、在不同类型土地覆盖地区进行气溶胶光学厚度、大气颗粒物浓度和气象数据的多点同步地面量测，并通过量测结果对大气颗粒物遥感反演模型进行校正和验证；
步骤五、使用校正和验证后的大气颗粒物遥感反演模型反演研究区的大气颗粒物浓度空间分布。
根据本申请的至少一个实施方式，在所述步骤三至步骤五中的所述大气颗粒物遥感反演模型是PM10的大气颗粒物遥感反演模型。
根据本申请的至少一个实施方式，在所述步骤二中，大气气溶胶标高HA来自于Calipso激光雷达数据的大气边界层高度数据，相对湿度的百分数RH来自于中尺度气象模式MM5模拟的气象数据。
根据本申请的至少一个实施方式，在所述步骤三中，所述大气颗粒物站点监测数据是通过Block-Kriging插值处理后的数据。
根据本申请的至少一个实施方式，在所述步骤三中，在建立不同类型土地覆盖和不同季节的大气颗粒物遥感反演模型步骤中，还包括：
将研究区内预定数量的大气颗粒物监测站点建立1km*1km的缓冲区，用缓冲区图层裁剪MODIS土地覆盖类型图，将缓冲区内所占面积最大的土地覆盖类型指定给对应的大气颗粒物监测站点，并计算出多年内各类型土地覆盖的平均大气颗粒物浓度值；以及
采取Tukey’s多重检验法比较各类型土地覆盖的平均大气颗粒物浓度值的两两差异，将无显著差异的土地覆盖类型合并，提取出n类土地覆盖，并选取有代表性的小区域，使用高分辨率卫星遥感数据进行图像解译，从而精确识别土地覆盖类型。
根据本申请的至少一个实施方式，在所述步骤四中，是根据Angstrom波长指数反映的研究区粒径分布，在不同类型土地覆盖地区选取多个采样点和多条采样线，进行气溶胶光学厚度、大气颗粒物浓度和气象数据的多点同步地面量测。
本申请至少存在以下有益技术效果：
本申请的气溶胶光学厚度反演方法及大气颗粒物遥感反演方法中，通过改进的暗像元法进行反演，在反演时，对暗像元进行了改进，使用包含预定数量像元的超像元窗口，通过改变窗口的移动方式、使用适合研究区的气溶胶类型以及使用四维查找表的存储方式等，增加了窗口内参与计算的暗像元个数，有效减少反演空白值，提高500米分辨率气溶胶光学厚度的反演精度和效率；
进一步，在大气颗粒物遥感反演方法中，由于采用上述反演的高精度的500米分辨率气溶胶光学厚度，从而更加适合小区域尺度上的空气质量监测，应用性强；
进一步，在大气颗粒物遥感反演方法中，通过对气溶胶光学厚度进行垂直和湿度订正，从而得到干空气状态下的地表气溶胶消光系数，再与大气颗粒物站点监测数据结合，以建立大气颗粒物遥感反演模型，从而能够消除监测站点与卫星遥感数据的空间尺度差异，提高500米分辨率的大气颗粒物遥感反演精度；并且，对气溶胶光学厚度进行垂直和湿度订正的分别来自Calipso激光雷达数据以及中尺度气象模式MM5模拟的气象数据，空间分辨率更高，同样能够提高500米分辨率的大气颗粒物遥感反演精度；
进一步，在大气颗粒物遥感反演方法中，能够准确建立不同土地利用类型和不同季节里的AOT-PM关系模型，从而使得大气颗粒物遥感反演的时空适用性更强，反演精度更高。
附图说明
图1是本申请气溶胶光学厚度反演方法的技术流程图；
图2是本申请气溶胶光学厚度反演方法中超像元示意图；
图3a-3f是通过本申请气溶胶光学厚度反演方法得到的，珠三角地区2010年11-12月期间6个不同时段的气溶胶光学厚度结果图；
图4是本申请气溶胶光学厚度反演方法得到的AOD(其中，AOD等于AOT)与太阳光度计实地监测数据对比结果图；
图5a和图5b分别是本申请气溶胶光学厚度反演方法得到的AOD与不同MODIS AOD产品对比结果图；
图6是本申请大气颗粒物遥感反演方法的技术流程图；
图7是本申请大气颗粒物遥感反演方法中AOT与地面监测PM10的相关性示意图；
图8是通过本申请大气颗粒物遥感反演方法得到的500米分辨率的珠三角地区年均PM10浓度空间分布图；
图9是通过本申请大气颗粒物遥感反演方法得到的500米分辨率的珠三角15个站点的PM10实测年均值和模型回归年均值的对比图。
具体实施方式
为使本申请实施的目的、技术方案和优点更加清楚，下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行更加详细的描述。在附图中，自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。所描述的实施例是本申请一部分实施例，而不是全部的实施例。下面通过参考附图描述的实施例是示例性的，旨在用于解释本申请，而不能理解为对本申请的限制。基于本申请中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本申请保护的范围。下面结合附图对本申请的实施例进行详细说明。
下面结合附图1-图9对本申请的气溶胶光学厚度反演方法及大气颗粒物遥感反演方法做进一步详细说明。
第一方面，本申请公开了一种气溶胶光学厚度反演方法，为了能较为准确地获取(较高分辨率)500米分辨率的大气颗粒物空间分布，对其中的MOD04-C0005暗像元算法做了改进：对每一个被计算像元取超像元，超像元被计算出来的光学厚度值作为目标像元的光学厚度值。
具体的，本申请的气溶胶光学厚度反演方法可以包括如下步骤：
步骤一、获取研究区(具体地区可以根据需要进行适合的选择)的无云或相对少云覆盖、分辨率为500米的MODIS L1B遥感数据(Terra和Aqua)，并对MODIS L1B遥感数据进行订正处理；
其中，对MODIS L1B遥感数据进行的具体订正处理可以根据需要进行适合的选择，本实施例中，对MODIS L1B遥感数据进行订正处理包括：对云和水域进行掩膜处理，以及通过如下公式(1)-(3)进行水汽、臭氧、CO2订正处理：
其中，
表示表示水汽透过率；/>
表示臭氧透过率；/>
表示二氧化碳透过率；G表示大气质量因子；/>
表示水汽吸收系数；/>
表示臭氧吸收系数；/>
表示二氧化碳吸收系数。
步骤二、在MODIS L1B遥感数据中，以左上角第一个像元为被计算像元，建立分辨率为10*10km且包含N*N个像元的超像元窗口；
其中，本申请中出现的“*”号，表示相乘，后续不再赘述；另外，N的数值可以根据需要进行选择，本实施例中，如图2所示，N为20，即建立包含20*20个像元的超像元窗口。
步骤三、在超像元窗口中，针对2.13μm波段像元，以表观反射率大于等于0.01且小于等于0.25为条件从中提取像元，以得到第一个暗像元集合。
步骤四、在超像元窗口中，针对0.66μm波段像元，去掉该波段亮度的前50％和后20％，以形成第二个像元集合。
步骤五、将步骤三中的第一个暗像元集合和步骤四中的第二个像元集合求交，得到有效暗像元集合；以20*20个像元的超像元窗口为例，这样求交后超像元中剩下来最多120个有效像元。
步骤六、判断所述有效暗像元集合中的有效暗像元数量是否超过M个；如果超过，则进入步骤七；否则，返回步骤二，且将被计算像元的气溶胶光学厚度值设为-1；其中，M的数值同样可以根据需要进行选择，本实施例中优选为12。
步骤七、在有效暗像元集合中，计算波段为0.47μm、0.66μm、1.24μm和2.13μm的各暗像元的表观反射率以及NDVI值，并通过关系模型(参见图1左下角部分记载内容)计算出波段为0.47μm和0.66μm的地表反射率。
步骤八、根据6S辐射传输模型、MODIS L1B遥感数据的几何参数以及研究区的气溶胶类型参数，建立四维查找表；
其中，6S(Second Simulation of the Satellite Signal in theSolarSpectrum)辐射传输模型是一种比较成熟的大气订正模型，采用SOS(successiveorder of scattering)的散射计算方法，准确模拟大气如何影响辐射在太阳一目标物一遥感器路径之间的传输，应用范围较广，主要包括以下几个部分参数：太阳、地物与传感器之间的几何参数；大气模式、气溶胶模式、传感器的光谱特性、地表反射率。应用该模型进行大气校正的工作流程是：按顺序将几何参数、大气模式、气溶胶模式、可见度、海拔高度、观测波段等参数输入6S模型，运行6S模型得出大气校正系数Xa、Xb、Xc，然后根据大气校正系数Xa、Xb、Xc计算得到地表反射率。将设置的不同参数输入到该模型中可以得到不同条件下的地表反射率查找表，本发明通过将其它途径计算得到的地表反射率带入该模型中反演出模型中对应的其它参数。
在实际应用中，6S辐射传输模型的参数变量、MODIS L1B遥感数据的几何参数以及研究区的气溶胶类型参数均可以根据需要进行适合的选择；本发明中，根据研究区情况，6S辐射传输模型所需参数变量参见如下表1所示，包括通道、高程、大气顶层反射率以及气溶胶光学厚度：
表1利用6S辐射传输模式建立查找表所用参数变量
MODIS L1B遥感数据的几何参数包括太阳高度角和太阳天顶角。本发明未包括卫星高度角，因为本申请中所用的数据是MODIS的Level1B数据，而Level1B已经进行了辐射定标和几何校正，因此只需要太阳高度角和方位角，这样一来大大削减了查找表的容量，增加了运算的执行效率。
研究区的气溶胶类型如下表2所示，包括沙尘、水溶性、海洋、烟尘：
表2暗像元法中使用的华南地区气溶胶类型参数
| 沙尘 | 水溶性 | 海洋 | 烟尘 | |
| 组分比例 | 0.13 | 0.7 | 0.15 | 0.005 |
| 标准差 | 0.1 | 0.12 | 0.11 | 0.004 |
将通过上述步骤得到的参数及依据这些参数运行6S模型所得的地表反射率数据用四维数组存储得到四维查找表(读取数组中的数据，明显比从数据库中检索并读取表中数据要快很多)。
步骤九、将波段为0.47μm和0.66μm的地表反射率代入到四维查找表中，逐像元移动超像元窗口，反演出各像元在对应波段的气溶胶光学厚度值。
步骤十、在如下气溶胶光学厚度与波长的公式(4)，分别代入波长0.47μm和0.66μm的气溶胶光学厚度，计算出α和β值。再通过公式(4)反演出分辨率为500米、波段为0.55μm的气溶胶光学厚度：
τ＝βλ-α (4)；
其中τ为气溶胶光学厚度，λ为波长，β为Angstrom混浊系数，代表大气中气溶胶的浓度，α为Angstrom波长指数。
综上，本申请的气溶胶光学厚度反演方法在反演时，对暗像元进行了改进，使用了一个20*20个像元的超像元窗口，通过改变窗口的移动方式、使用适合研究区的气溶胶类型以及使用四维查找表的存储方式等，增加窗口内参与计算的暗像元个数，有效减少反演空白值，提高反演精度和效率。
参见图3a-图3f所示，是通过本申请气溶胶光学厚度反演方法得到的珠三角地区2010年11-12月期间6个不同时段的气溶胶光学厚度结果图；进一步参见图4以及图5a、图5b，通过本申请气溶胶光学厚度反演方法得到的珠三角地区气溶胶光学厚度(AOD)与MICROTOPSII太阳光度计的地面监测数据及MODIS AOT产品进行对比验证，结果显示本申请气溶胶光学厚度反演方法的反演结果精度较高。
第二方面，如图6所示，本申请还公开了一种大气颗粒物遥感反演方法，可以包括如下步骤：
步骤一、采用第一方面中任一项所述的气溶胶光学厚度反演方法获得研究区(具体地区可根据反演需求而定)的气溶胶光学厚度。
步骤二、对气溶胶光学厚度进行垂直和湿度订正，并通过如下公式(6)计算干空气状态下的地表气溶胶消光系数。
其中，AOT反映的是整个大气层的气溶胶光学厚度，而大气颗粒物(PM)浓度常为地表水平，因此两者之间的关系与气溶胶垂直分布有关，也与影响气溶胶散射系数的地面相对湿度有关。具体的，AOT与PM成正相关关系，并且这种关系通过垂直和相对湿度的方法订正后会有更大提高。
另外，AOT和PM这两个属性与大气剖面、环境条件、以及气溶胶的粒径分布和化学组分有关，且在空间和时间上存在异质性，因此直接通过气溶胶光学厚度来估算大气颗粒物PM10浓度有很大的不确定性。为了减小这种不确定性，本发明中大气边界层高度和相对湿度可按照公式(5)对气溶胶光学厚度AOT进行校正：
f(RH)＝(1-RH/100)-1 (5)；
kAOT,DRY(λ)＝τa(λ)/[HA·f(RH)] (6)；
其中，f(RH)代表相对湿度的校正系数，kAOT,DRY(λ)代表干空气状态下的地表气溶胶消光系数，τa(λ)代表通过遥感影像反演的气溶胶光学厚度，RH为相对湿度的百分数，HA为大气气溶胶标高，即假定气溶胶浓度随高度分布保持不变时的气溶胶层等效厚度。
进一步的，由于目前订正AOT的气象数据多为站点监测值或再分析资料，空间分辨率低；为此，本申请是使用Calipso激光雷达数据的大气边界层高度数据和中尺度气象模式MM5模拟的气象数据(主要涉及相对湿度数据)对气溶胶光学厚度进行垂直和湿度订正，以提高订正精度。
步骤三、将步骤二计算得到的干空气状态下的地表气溶胶消光系数与大气颗粒物站点监测数据结合，根据如下一元线性方程(7)，建立不同类型土地覆盖(例如郊区、城区)和不同季节(干湿季节)的大气颗粒物遥感反演模型(即AOT-PM模型，参见图7所示例)，其中，大气颗粒物遥感反演模型的分辨率为500米：
Y＝aX+b (7)；
其中，X为对应站点所在像元的订正后气溶胶光学厚度值(也即步骤二计算得到的干空气状态下的地表气溶胶消光系数kAOT,DRY(λ))；Y为PM10站点监测值；a、b为系数，是根据具体土地覆盖类型以及具体季节时的X、Y带入值求得。
另外，本实施例中大气颗粒物遥感反演模型采用PM10的大气颗粒物遥感反演模型，并且大气颗粒物站点监测数据是通过Block-Kriging插值处理后的数据，以减少数据的空间尺度差异给模型带来的误差。
进一步，为提高土地覆盖类型识别精确，在本步骤在建立不同类型土地覆盖和不同季节的大气颗粒物遥感反演模型过程中，还可以包括如下步骤：
将研究区内预定数量的大气颗粒物监测站点建立1km*1km的缓冲区，用缓冲区图层裁剪MODIS土地覆盖类型图，将缓冲区内所占面积最大的土地覆盖类型指定给对应的大气颗粒物监测站点，并计算出多年(例如5年)内各类型土地覆盖的平均大气颗粒物浓度值；
以及采取Tukey’s多重检验法比较各类型土地覆盖的平均大气颗粒物浓度值的两两差异，将无显著差异的土地覆盖类型合并，提取出n类土地覆盖，并选取有代表性的小区域，使用高分辨率卫星遥感数据进行图像解译，从而精确识别土地覆盖类型。
步骤四、在不同类型土地覆盖地区进行气溶胶光学厚度、大气颗粒物浓度和气象数据的多点同步地面量测，并通过量测结果对大气颗粒物遥感反演模型进行校正和验证。
具体的，是根据Angstrom波长指数反映的研究区粒径分布，在不同类型土地覆盖地区选取多个采样点和多条采样线，进行气溶胶光学厚度、大气颗粒物浓度和气象数据的多点同步地面量测。
步骤五、使用校正和验证后的大气颗粒物遥感反演模型，反演得到研究区的大气颗粒物浓度空间分布。
参见图8所示，是通过本申请大气颗粒物遥感反演方法得到分辨率为500米的珠三角地区年均PM10浓度空间分布图；进一步，参见通过本申请反演方法得到，分辨率为500米的珠三角15个站点的PM10实测年均值和模型回归年均值对比表3以及曲线对比图图9：
表3珠三角15个站点的PM10实测年均值和模型回归年均值对比
结果显示，本申请大气颗粒物遥感反演方法得到的PM10的回归值(即反演得到的值)与实测值偏差较小，说明通过本申请大气颗粒物遥感反演方法得到的反演结果较为准确。
综上，本申请的大气颗粒物遥感反演方法中，由于采用上述反演的高精度的500米分辨率气溶胶光学厚度，从而更加适合小区域尺度上的空气质量监测，应用性强；
进一步，在大气颗粒物遥感反演方法中，通过对气溶胶光学厚度进行垂直和湿度订正，从而得到干空气状态下的地表气溶胶消光系数，再与大气颗粒物站点监测数据结合，以建立大气颗粒物遥感反演模型，从而能够消除监测站点与卫星遥感数据的空间尺度差异，提高500米分辨率的大气颗粒物遥感反演精度；并且，对气溶胶光学厚度进行垂直和湿度订正的分别来自Calipso激光雷达数据以及中尺度气象模式MM5模拟的气象数据，空间分辨率更高，能够进一步提高500米分辨率的大气颗粒物遥感反演精度；
进一步，在大气颗粒物遥感反演方法中，能够准确建立不同土地覆盖类型和不同季节里的AOT-PM关系模型(参见图7所示，模型相关性良好，R2＝0.67)，从而使得大气颗粒物遥感反演的时空适用性更强，反演精度更高，以准确获取500米分辨率的大气颗粒物地面浓度空间分布。
以上所述，仅为本申请的具体实施方式，但本申请的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本申请揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本申请的保护范围之内。因此，本申请的保护范围应以所述权利要求的保护范围为准。
Claims (9)
1.一种气溶胶光学厚度反演方法，其特征在于，包括如下步骤：
步骤一、获取研究区的无云或相对少云覆盖、分辨率为500米的MODIS L1B遥感数据，并对所述MODIS L1B遥感数据进行订正处理；
步骤二、在所述MODIS L1B遥感数据中，以左上角第一个像元为被计算像元，建立分辨率为10*10km且包含N*N个像元的超像元窗口；
步骤三、在所述超像元窗口中，针对2.13μm波段像元，以表观反射率大于等于0.01且小于等于0.25为条件从中提取像元，以得到第一个暗像元集合；
步骤四、在所述超像元窗口中，针对0.66μm波段像元，去掉该波段亮度的前50％和后20％，以形成第二个像元集合；
步骤五、将所述第一个暗像元集合和所述第二个像元集合求交，得到有效暗像元集合；
步骤六、判断所述有效暗像元集合中的有效暗像元数量是否超过M个；如果超过，则进入步骤七；否则，返回步骤二，且将被计算像元的气溶胶光学厚度值设为-1；
步骤七、在所述有效暗像元集合中，计算波段为0.47μm、0.66μm、1.24μm和2.13μm的表观反射率以及NDVI值，并计算出波段为0.47μm和0.66μm的地表反射率；
步骤八、根据6S辐射传输模型、MODIS L1B遥感数据的几何参数以及研究区的气溶胶类型参数，建立四维查找表；
步骤九、将波段为0.47μm和0.66μm的地表反射率代入到所述四维查找表中，逐像元移动超像元窗口，反演出各像元在对应波段的气溶胶光学厚度值；
步骤十、在气溶胶光学厚度与波长的如下公式(4)中，先引入波段0.47μm和0.66μm的气溶胶光学厚度，以得到α和β值，再通过该公式反演出分辨率为500米、波段为0.55μm的气溶胶光学厚度：
τ＝βλ-α (4)；
其中τ为气溶胶光学厚度，λ为波长，β为Angstrom混浊系数，代表大气中气溶胶的浓度，α为Angstrom波长指数；
其中在所述步骤一中，对所述MODIS L1B遥感数据进行的所述订正处理包括对云和水域进行掩膜处理，以及通过如下公式(1)-(3)进行水汽、臭氧、CO2订正处理：
其中，
表示水汽透过率；/>
表示臭氧透过率；/>
表示二氧化碳透过率；G表示大气质量因子；/>
表示水汽吸收系数；/>
表示臭氧吸收系数；/>
表示二氧化碳吸收系数。
2.根据权利要求1所述的气溶胶光学厚度反演方法，其特征在于，在所述步骤二中，是建立包含20*20个像元的超像元窗口，且在所述步骤六中，是判断有效暗像元数量是否超过12个。
3.根据权利要求1所述的气溶胶光学厚度反演方法，其特征在于，在所述步骤八中，所述6S辐射传输模型所需参数变量包括通道、高程、大气顶层反射率以及光学厚度；以及
所述MODIS L1B遥感数据的几何参数包括太阳高度角和太阳天顶角；以及
所述研究区的气溶胶类型参数包括沙尘、水溶性、海洋和烟尘。
4.一种大气颗粒物遥感反演方法，其特征在于，包括如下步骤：
步骤一、采用如权利要求1-3任一项所述的气溶胶光学厚度反演方法获得研究区的气溶胶光学厚度；
步骤二、对气溶胶光学厚度进行垂直和湿度订正，并根据以下公式(6)计算干空气状态下的地表气溶胶消光系数：
kAOT,DRY(λ)＝τa(λ)/[HA·f(RH)] (6)；
其中，f(RH)代表相对湿度的校正系数，f(RH)＝(1-RH/100)-1；kAOT,DRY(λ)代表干空气状态下的地表气溶胶消光系数，τa(λ)代表通过遥感影像反演的气溶胶光学厚度，RH为相对湿度的百分数，HA为大气气溶胶标高，即假定气溶胶浓度随高度分布保持不变时的气溶胶层等效厚度；
步骤三、将所述干空气状态下的地表气溶胶消光系数与大气颗粒物站点监测数据结合，建立不同类型土地覆盖和不同季节的大气颗粒物遥感反演模型(7)：
Y＝aX+b (7)；
其中，X为对应站点所在像元的订正后气溶胶光学厚度值，即kAOT,DRY(λ)；Y为PM10站点监测值；a、b为系数，是根据具体类型土地覆盖以及具体季节时的X、Y带入值求得；另外，所述大气颗粒物遥感反演模型的分辨率为500米；
步骤四、在不同类型土地覆盖地区进行气溶胶光学厚度、大气颗粒物浓度和气象数据的多点同步地面量测，并通过量测结果对大气颗粒物遥感反演模型进行校正和验证；
步骤五、使用校正和验证后的大气颗粒物遥感反演模型反演研究区的大气颗粒物浓度空间分布。
5.根据权利要求4所述的大气颗粒物遥感反演方法，其特征在于，在所述步骤三至步骤五中的所述大气颗粒物遥感反演模型是PM10的大气颗粒物遥感反演模型。
6.根据权利要求5所述的大气颗粒物遥感反演方法，其特征在于，在所述步骤二中，大气气溶胶标高HA来自于Calipso激光雷达数据的大气边界层高度数据，相对湿度的百分数RH来自于中尺度气象模式MM5模拟的气象数据。
7.根据权利要求4至6任一项所述的大气颗粒物遥感反演方法，其特征在于，在所述步骤三中，所述大气颗粒物站点监测数据是通过Block-Kriging插值处理后的数据。
8.根据权利要求6所述的大气颗粒物遥感反演方法，其特征在于，在所述步骤三中，在建立不同类型土地覆盖和不同季节的大气颗粒物遥感反演模型步骤中，还包括：
将研究区内预定数量的大气颗粒物监测站点建立1km*1km的缓冲区，用缓冲区图层裁剪MODIS土地覆盖类型图，将缓冲区内所占面积最大的土地覆盖类型指定给对应的大气颗粒物监测站点，并计算出多年内各类型土地覆盖的平均大气颗粒物浓度值；以及
采取Tukey’s多重检验法比较各类型土地覆盖的平均大气颗粒物浓度值的两两差异，将无显著差异的土地覆盖类型合并，提取出n类土地覆盖，并选取有代表性的小区域，使用高分辨率卫星遥感数据进行图像解译，从而精确识别土地覆盖类型。
9.根据权利要求8所述的大气颗粒物遥感反演方法，其特征在于，在所述步骤四中，是根据Angstrom波长指数反映的研究区粒径分布，在不同类型土地覆盖地区选取多个采样点和多条采样线，进行气溶胶光学厚度、大气颗粒物浓度和气象数据的多点同步地面量测。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110266804.6A CN113740263B (zh) | 2021-03-10 | 2021-03-10 | 一种气溶胶光学厚度反演方法及大气颗粒物遥感反演方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110266804.6A CN113740263B (zh) | 2021-03-10 | 2021-03-10 | 一种气溶胶光学厚度反演方法及大气颗粒物遥感反演方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113740263A CN113740263A (zh) | 2021-12-03 |
| CN113740263B true CN113740263B (zh) | 2023-06-16 |
Family
ID=78728243
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110266804.6A Active CN113740263B (zh) | 2021-03-10 | 2021-03-10 | 一种气溶胶光学厚度反演方法及大气颗粒物遥感反演方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113740263B (zh) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117347282B (zh) * | 2023-08-22 | 2024-05-28 | 中南大学 | 星基气溶胶光学厚度反演方法、装置及系统和存储介质 |
Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101598543A (zh) * | 2009-07-29 | 2009-12-09 | 中国科学院对地观测与数字地球科学中心 | 一种实用的遥感影像大气校正方法 |
| CN106407656A (zh) * | 2016-08-29 | 2017-02-15 | 中国科学院遥感与数字地球研究所 | 一种基于高分辨率卫星影像数据的气溶胶光学厚度反演方法 |
| CN109030301A (zh) * | 2018-06-05 | 2018-12-18 | 中南林业科技大学 | 基于遥感数据的气溶胶光学厚度反演方法 |
| CN110186823A (zh) * | 2019-06-26 | 2019-08-30 | 中国科学院遥感与数字地球研究所 | 一种气溶胶光学厚度反演方法 |
| CN110261873A (zh) * | 2019-05-29 | 2019-09-20 | 电子科技大学 | 一种基于分段统计的大气气溶胶反演方法 |
| CN110673229A (zh) * | 2019-10-23 | 2020-01-10 | 新亚优华科技有限公司 | 一种基于热点网格技术的大气污染物扩散轨迹追踪方法 |
| CN111753439A (zh) * | 2020-07-09 | 2020-10-09 | 中国科学院空天信息创新研究院 | 一种国产多角度偏振卫星传感器的气溶胶光学厚度反演方法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105674901B (zh) * | 2016-04-08 | 2017-12-26 | 电子科技大学 | 基于统计学分段的大气气溶胶反演方法 |
-
2021
- 2021-03-10 CN CN202110266804.6A patent/CN113740263B/zh active Active
Patent Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101598543A (zh) * | 2009-07-29 | 2009-12-09 | 中国科学院对地观测与数字地球科学中心 | 一种实用的遥感影像大气校正方法 |
| CN106407656A (zh) * | 2016-08-29 | 2017-02-15 | 中国科学院遥感与数字地球研究所 | 一种基于高分辨率卫星影像数据的气溶胶光学厚度反演方法 |
| CN109030301A (zh) * | 2018-06-05 | 2018-12-18 | 中南林业科技大学 | 基于遥感数据的气溶胶光学厚度反演方法 |
| CN110261873A (zh) * | 2019-05-29 | 2019-09-20 | 电子科技大学 | 一种基于分段统计的大气气溶胶反演方法 |
| CN110186823A (zh) * | 2019-06-26 | 2019-08-30 | 中国科学院遥感与数字地球研究所 | 一种气溶胶光学厚度反演方法 |
| CN110673229A (zh) * | 2019-10-23 | 2020-01-10 | 新亚优华科技有限公司 | 一种基于热点网格技术的大气污染物扩散轨迹追踪方法 |
| CN111753439A (zh) * | 2020-07-09 | 2020-10-09 | 中国科学院空天信息创新研究院 | 一种国产多角度偏振卫星传感器的气溶胶光学厚度反演方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113740263A (zh) | 2021-12-03 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| He et al. | Validation of MODIS derived aerosol optical depth over the Yangtze River Delta in China | |
| CN111795936B (zh) | 一种基于查找表的多光谱遥感影像大气校正系统、方法及存储介质 | |
| CN111579504B (zh) | 基于光学遥感的大气污染成分垂直分布反演方法 | |
| Leon et al. | Large‐scale advection of continental aerosols during INDOEX | |
| CN110411918B (zh) | 一种基于卫星偏振技术的pm2.5浓度遥感估算方法 | |
| CN112434617A (zh) | 一种基于多源遥感数据的森林生物量变化监测方法及系统 | |
| CN113744249B (zh) | 一种海洋生态环境损害调查方法 | |
| CN111415309A (zh) | 一种基于最小反射率法的高分辨率遥感影像大气校正方法 | |
| CN111079835B (zh) | 一种基于深度全连接网络的Himawari-8大气气溶胶反演方法 | |
| CN111191380B (zh) | 一种基于地基光谱仪测量数据的大气气溶胶光学厚度估算方法和装置 | |
| Lee et al. | Analytical approach to estimating aerosol extinction and visibility from satellite observations | |
| CN113740263B (zh) | 一种气溶胶光学厚度反演方法及大气颗粒物遥感反演方法 | |
| CN112131789A (zh) | 一种基于随机森林算法的多光谱降水检测系统及方法 | |
| CN116664947A (zh) | 一种基于卫星观测数据的蓝藻水华监测方法及系统 | |
| CN111650128B (zh) | 一种基于地表反射率库的高分辨率大气气溶胶反演方法 | |
| Burley et al. | A fast two-stream-like multiple-scattering method for atmospheric characterization and radiative transfer | |
| Lyapustin et al. | A comparison of multi-angle implementation of atmospheric correction and MOD09 daily surface reflectance products from MODIS | |
| CN111723525A (zh) | 一种基于多源数据和神经网络模型的pm2.5反演方法 | |
| Choi et al. | GOCI Yonsei aerosol retrieval version 2 aerosol products: improved algorithm description and error analysis with uncertainty estimation from 5-year validation over East Asia | |
| CN116558652A (zh) | 一种气象卫星热红外通道在轨辐射定标方法 | |
| Richardson et al. | Boundary layer water vapour statistics from high-spatial-resolution spaceborne imaging spectroscopy | |
| CN108132096B (zh) | 一种基于激光雷达的林窗太阳辐射监测方法 | |
| CN108519058B (zh) | 一种光学遥感图像估算一年生海冰厚度的方法 | |
| He et al. | Modis aerosol optical thickness product algorithm verification and analysis | |
| CN113034519B (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 |