CN110909309A - 一种逐小时高分辨率pm2.5数据的获取方法 - Google Patents

一种逐小时高分辨率pm2.5数据的获取方法 Download PDF

Info

Publication number
CN110909309A
CN110909309A CN201911150572.7A CN201911150572A CN110909309A CN 110909309 A CN110909309 A CN 110909309A CN 201911150572 A CN201911150572 A CN 201911150572A CN 110909309 A CN110909309 A CN 110909309A
Authority
CN
China
Prior art keywords
data
concentration
station
aerosol
optical thickness
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.)
Granted
Application number
CN201911150572.7A
Other languages
English (en)
Other versions
CN110909309B (zh
Inventor
阎福礼
韩春明
范湘涛
白林燕
丁翼星
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Aerospace Information Research Institute of CAS
Original Assignee
Institute of Remote Sensing and Digital Earth of CAS
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Institute of Remote Sensing and Digital Earth of CAS filed Critical Institute of Remote Sensing and Digital Earth of CAS
Priority to CN201911150572.7A priority Critical patent/CN110909309B/zh
Publication of CN110909309A publication Critical patent/CN110909309A/zh
Application granted granted Critical
Publication of CN110909309B publication Critical patent/CN110909309B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • 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
    • 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

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Chemical & Material Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Probability & Statistics with Applications (AREA)
  • Evolutionary Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Dispersion Chemistry (AREA)
  • General Health & Medical Sciences (AREA)
  • Pathology (AREA)
  • Immunology (AREA)
  • Algebra (AREA)
  • Health & Medical Sciences (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Operations Research (AREA)
  • Analytical Chemistry (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Biochemistry (AREA)
  • Investigating Or Analysing Materials By Optical Means (AREA)

Abstract

本发明公开一种逐小时高分辨率PM2.5数据的获取方法，包括如下步骤：基于相关计算的气溶胶光学厚度数据填补；基于PM2.5浓度季节变化的遥感反演算法模型；基于空间权重的PM2.5遥感监测数据的多时相插值；基于台站数据与遥感数据的信息融合。高效率的填补缺失的气溶胶光学厚度数据，能够提升大气污染物的连续观测能力；建立能反映季节变化和大气边界层变化规律的地面PM2.5浓度反演模型，获取具有逐小时特征的高空间分辨率PM2.5浓度数据，使得高空间分辨率的遥感监测数据达到了逐小时特征的高时间分辨率，台站、遥感数据的信息融合，能够建立具有逐小时、高空间分辨率特征的精准PM2.5浓度数据集。

Description

一种逐小时高分辨率PM2.5数据的获取方法
技术领域
本发明涉及环境遥感技术领域。具体地说是一种逐小时高分辨率PM2.5数据的获取方法。
背景技术
PM2.5又称细粒、细颗粒、细颗粒物，是指大气环境中空气动力学当量直径小于等于2.5微米的颗粒物，是我国大部分城市环境空气质量的首要污染物(郦嘉诚，2018)。目前的监测手段，主要是通过大量布设地面台站，实测地面PM2.5的浓度，能够获取每小时的大气污染数据，具有监测连续、数据精准的技术优势，但是地面台站仅仅是离散点的观测方式，难以获取高空间分辨率的PM2.5浓度分布；利用遥感手段获取的气溶胶光学厚度(AOD)数据反演大气PM2.5浓度的二维空间分布，具有高空间分辨率的技术优势，但是TERRA、AQUA卫星联合也只能每天最多获取2期数据(TERRA卫星上午10:30前后过境；AQUA卫星下午2:00前后过境)，无法实现逐小时连续的PM2.5监测。要研究并连续观测大气PM2.5浓度的空间变化规律，有必要建立逐小时的高空间分辨率的大气PM2.5数据集。此外，由于遥感获取气溶胶光学厚度(AOD)受云团零星覆盖或完全覆盖等影响严重，造成气溶胶光学厚度数据的局部孔洞或整体缺失，如何对气溶胶光学厚度数据的缺失部分进行有效填补，这是构建连续PM2.5数据集的重要需求和基础。
近年来，利用中分辨率成像光谱仪MODIS和多角度成像光谱仪MISR和Geochem模型反演的年度PM2.5效果良好，但是每天的PM2.5不确定性极大(Donkelaar，2010)。STRAWA(2013)发展了采用OMI和MODIS结合的方式。目前，利用遥感技术反演PM2.5浓度计算方法大致分为四种：单线性模型、多元线性回归模型、人工智能模型和广义加法模型(马宗伟，2015)。这些方法监测范围大、成本低、可长时间动态监测，有效弥补了地面监测的不足，是目前PM2.5浓度监测的主要手段。但是大气污染中细颗粒物的光学特性与颗粒物的物质组成、粒径、温度、湿度等密切相关，相应的与温度、降水、地面秸秆焚烧、工业生产、供暖等都具有显著的季节变化现象，缺乏能够考虑季节差异的PM2.5浓度遥感反演方法。
多年来观测数据与模式数据的同化研究方面发展迅速，在国家气象业务中得到充分论证和应用，成效显著(梁顺林，2013；姜立鹏等，2013)。Terra和AQUA卫星的同类观测数据之间的信息融合工作也进展迅速(陈韵竹，2013)。遥感监测的PM2.5数据受颗粒物物质组成、性质变化等有季节性规律的波动，会有反演精度的较大波动，需要利用地面台站数据进行反演结果的信息融合校准。
站点数据的空间化，最为常用的是台站数据内插。如吴亮(2016)利用ARCGIS采用克里金插值方法对站点PM2.5数据的空间化；大气数据的空间化插值，往往需要考虑地理因素的影响。如廖顺宝(2003)考虑了经纬度、海拔的规律性和其他非规律性影响因素，建立了多元线性回归和反距离权重方法，对气温数据进行的空间化。大气PM2.5的变化也需要充分考虑大气物理模式作用下的浓度空间变化规律，这是台站数据空间化的重要控制性因素。
发明内容
为此，本发明所要解决的技术问题台站观测数据缺乏高分辨率空间属性，遥感观测数据缺乏逐小时的时间特征，提供一种能够获取逐小时高分辨率PM2.5数据的获取方法。
为解决上述技术问题，本发明提供如下技术方案：
一种逐小时高分辨率PM2.5数据的获取方法，包括如下步骤：
(1)基于相关计算的气溶胶光学厚度数据填补；
(2)基于PM2.5浓度季节变化的遥感反演算法模型；
(3)基于空间权重的PM2.5遥感监测数据的多时相插值；
(4)基于台站数据与遥感数据的信息融合。
上述逐小时高分辨率PM2.5数据的获取方法，在步骤(1)中，
(1-1)对零星云团覆盖，气溶胶光学厚度有孔洞的地区，填补缺失部分的数据：对待填补气溶胶数据与存档数据的二维空间相关性计算，通过相关性阈值遴选合适的存档数据进行气溶胶光学厚度孔洞的填补；
(1-2)对于云团完全覆盖、无气溶胶光学厚度数据区域，填补缺失数据：对大面积无值区域，采用待填补时相的台站数据与台站存档数据进行相关性计算，遴选高相关时相的气溶胶光学厚度数据，填补无数据区域；
(1-3)进行加权平均，赋值无数据像元。
上述逐小时高分辨率PM2.5数据的获取方法，在步骤(1-1)中：实行二维空间相关的自适应滑动窗口，在待填补像元位置(Xi，Yi)周边搜索，至少保证10×10有效像元数量，遍历历史存档的气溶胶光学厚度数据，遴选，并确定平均相关系数在0.4-1.0之间的历史存档数据：
二维图像的孔洞像元填补：对于小面积缺失的气溶胶光学厚度孔洞，缺失像元的搜索半径在5-20个像元以内，选择待填补区域的单个像元P(a,b)，a，b为像元行列号，遍历其周围最邻近非零像元，确定大于10×10的像元的行列坐标，读取像元数值X(i，j)；搜索长时间序列气溶胶光学厚度的历史存档数据，找到对应坐标的像元数值Y(i，j)，进行计算：
Figure BDA0002283416950000031
式(Ⅰ)中，R为相关系数；X(i，j)为待填补窗口非零像元的数值(i＝1….10,j＝1…10)；Y(i,j)为存档气溶胶数据中相对应坐标的非零像元数值(i＝1….10,j＝1…10)，
Figure BDA0002283416950000032
分别为待填补窗口中非零像元均值、对应非零像元位置的存档数据非零像元均值；
计算得到的R值位于0.4-1.0之间，平均误差位于0-50％的存档气溶胶光学厚度数据；选出符合标准的多期气溶胶光学厚度数据。
上述逐小时高分辨率PM2.5数据的获取方法，在步骤(1-2)中：台站数据与台站数据的相关性计算，确定卫星过境相同时间的站点数据，获得相关系数在0.4-1.0之间的时相，选定存档的气溶胶光学厚度历史数据；完全无数据区域填补：对于大面积缺失数据的区域，搜索半径距离大于20个像元，利用卫星过境前后30-90分钟以内的n个监测台站的3n个监测数据，对于待填补的卫星过境时间，每个台站有三个监测时刻，三个监测数据；计算3n个监测数据与历史存档数据的相关性和误差大小进行填补；
台站数据相关计算：
Figure BDA0002283416950000041
式(Ⅱ)中，R为相关系数；xi、yi分别为待填补时间的台站监测数据、存档台站数据中监测数值，i＝1….3n,n为监测台站数量，
Figure BDA0002283416950000042
分别为待填补时刻的所有监测数据均值、对应时刻的存档数据均值；
计算得到的R值位于0.4-1.0之间，平均误差位于0-60％的存档数据，记录相应的成像日期，获取该时刻的对应气溶胶光学厚度数据。
上述逐小时高分辨率PM2.5数据的获取方法，在步骤(1-3)中，对步骤(1-1)和步骤(1-2)中的多期次的气溶胶光学厚度数据，进行加权平均计算：
Figure BDA0002283416950000043
式(Ⅲ)中，像元P(a,b)为待填补像元数值；a，b为像元行列号；
Figure BDA0002283416950000044
为步骤(1-1)或步骤(1-2)得到的n期气溶胶光学厚度存档数据的和；n为非零像元的期数。
上述逐小时高分辨率PM2.5数据的获取方法，在步骤(2)中，利用大气边界层变化函数，修正地面气溶胶光学厚度，构建的逐月的PM2.5多时相反演模型的方法和精度变化规律，确定的高敏感月份和算法形式；
具体高空间分辨率PM2.5浓度数据的计算步骤为：基于步骤(1)获取的全覆盖气溶胶光学厚度，利用细颗粒物光学特性季节规律建立PM2.5大气柱浓度模型，结合大气边界层变化规律，计算大气层底部的PM2.5浓度数据。
上述逐小时高分辨率PM2.5数据的获取方法，基于全覆盖气溶胶光学厚度AOD数据和实测PM2.5数据，建立的季节性PM2.5浓度反演模型为：
1月份模型μg/m3:PM2.5＝37.15+78.21×AOD；
2月份模型μg/m3：PM2.5＝23.47+63.18×AOD；
3月份模型μg/m3：PM2.5＝20.90+46.63×AOD；
4月份模型μg/m3：PM2.5＝2.90+65.93×AOD；
5月份模型μg/m3：PM2.5＝21.57+42.00×AOD；
6月份模型μg/m3：PM2.5＝14.14+34.96×AOD；
7月份模型μg/m3：PM2.5＝19.75+33.90×AOD；
8月份模型μg/m3：PM2.5＝17.10+30.66×AOD；
9月份模型μg/m3：PM2.5＝13.34+38.34×AOD；
10月份模型μg/m3：PM2.5＝22.25+46.61×AOD；
11月份模型μg/m3：PM2.5＝21.24+65.97×AOD；
12月份模型μg/m3：PM2.5＝22.12+107.88×AOD；
PM2.5与全覆盖气溶胶光学厚度AOD数据在12个月份的相关系数平均为0.66，决定系数R2为0.43，N＝1009，超过统计学上95％的置信区间，相关显著。
上述逐小时高分辨率PM2.5数据的获取方法，在步骤(3)中，高时间分辨率数据特征的获取，包括：
(3-1)基于步骤(2)获取的PM2.5浓度数据，利用大气模式模型，结合台站监测的风力、风向实测数据、高程数据，建立基于大气模式模拟的PM2.5浓度变化权重方法；
(3-2)通过当天的遥感监测TERRA卫星的PM2.5数据、AQUA卫星的PM2.5数据，联合第二天的TERRA卫星、AQUA卫星监测的PM2.5浓度监测数据，利用步骤(3-1)的大气模式模拟的PM2.5浓度权重变化，构建长时间序列的时相插值方法，计算得到空间分布函数，实现时相插值，建立高空间分辨率数据的PM2.5浓度数据，获得24幅逐小时的PM2.5浓度数据。
上述逐小时高分辨率PM2.5数据的获取方法，所述时相插值方法的模型为：
Figure BDA0002283416950000061
式(Ⅳ)中，PM2.5T时刻插值为T小时的PM2.5浓度分布；PM2.50时刻、PM2.524小时分别为起始0时刻、第24小时的终点时刻由全覆盖气溶胶光学厚度AOD数据反演的大气PM2.5浓度；PM2.5T时刻模拟值、PM2.524小时大气模拟值分别为T小时、第24小时大气模式模拟的大气污染分布数值；其中T＝1，2，3…24小时，PM2.50时刻、PM2.524小时、PM2.5T时刻模拟值、PM2.524小时大气模拟值、PM2.5T时刻插值的单位均为微克/立方米。
上述逐小时高分辨率PM2.5数据的获取方法，在步骤(4)中，基于步骤(3)获取的24幅PM2.5浓度数据，利用24个小时的台站点数据和24幅逐小时PM2.5空间分布数据，实现台站记录的24个小时的点观测数据和24幅PM2.5浓度的二维面数据，进行点-面距离方法的信息融合，实现逐小时数据的空间化过程，提高高空间分辨率数据的精度，形成每小时的PM2.5浓度分布数据；
采用遥感PM2.5浓度图像的空间分布的距离，遥感数据浓度分布隐含的方向角度，作为权重函数为：
PM2.5(融合)(i+Δi，J+Δj)＝PM2.5(卫星)(i，j)+[PM2.5(台站)-PM2.5(卫星)(i+Δi，j+Δj)]W(Δi,Δj)(Ⅴ)；
式(Ⅴ)中：Δi,Δj分别为距离台站位置的横向像元数、纵向像元数；其中w(Δi,Δj)为：
Figure BDA0002283416950000062
i，j分别为台站所在的像元的行、列数；PM2.5(融合)(i+Δi，J+Δj)为距离台站所在像元(i，j)的距离为(Δi，Δj)的PM2.5浓度融合结果；PM2.5(卫星)(i，j)为台站所在位置的卫星观测的PM2.5浓度数据；PM2.5(台站)为台站监测的PM2.5浓度数据；PM2.5(卫星)(i+Δi，j+Δj)为距台站所在像元(i，j)的距离为(Δi，Δj)的卫星观测PM2.5浓度数据；
能够计算距离台站位置5km半径内的点源与面源信息的融合。
本发明的技术方案取得了如下有益的技术效果：
针对上述现有工作进展和研究需求，本申请利用TERRA\AQUA卫星的气溶胶数据，开展了基于空间相关计算的长时间序列AOD数据的填补、基于季节变更的PM2.5浓度反演、最终结合WRF大气模式模拟权重，实现高空间分辨率数据的时相插值、实现点面信息融合，形成具有逐小时特征的高空间分辨率PM2.5监测数据集。
本发明的逐小时高分辨率PM2.5数据的获取方法。在气溶胶光学厚度存在孔洞的区域，通过建立基于二维图像相关性计算，遍历积累的长时间序列气溶胶光学厚度数据，找到气溶胶分布的高相关区域，填补缺失的气溶胶光学厚度数据；在气溶胶光学厚度完全缺失的时相，利用台站数据，建立相关性计算方法，遍历积累的站点数据，确定高相关时相，找到相应多个时相的气溶胶光学厚度数据，计算并代替缺失数据；高效率的填补缺失的气溶胶光学厚度数据，能够提升大气污染物的连续观测能力；根据台站数据、遥感监测的气溶胶光学厚度数据的季节变化，建立能反映季节变化和大气边界层变化规律的地面PM2.5浓度反演模型，并对TERRA、AQUA卫星监测的多时相PM2.5浓度数据，通过大气模式的权重约束下，进行时相插值，获取具有逐小时特征的高空间分辨率PM2.5浓度数据，使得高空间分辨率的遥感监测数据达到了逐小时特征的高时间分辨率；台站、遥感数据的信息融合，能够建立具有逐小时、高空间分辨率特征的精准PM2.5浓度数据集。为源解析、排放清单以及精细大气联防联控政策的制定和成效评价提供了重要数据支撑。
附图说明
图1本发明逐小时高分辨率PM2.5数据的获取方法的流程框示意图；
图2本发明逐小时高分辨率PM2.5数据的获取方法的数据计算流程示意图。
具体实施方式
获取逐小时高分辨率PM2.5数据的方法，利用空间相关计算填补气溶胶光学厚度缺失数据，实现季节性PM2.5浓度反演，实现遥感数据时相内插，完成台站数据和遥感数据信息融合，获取逐小时高分辨率的PM2.5数据。
如图1所示：能够利用存档的气溶胶光学厚度数据填补待填补的气溶胶光学厚度孔洞或大面积缺失数据，通过一系列计算模型和方法，能够获得逐小时高空间分辨率的PM2.5数据。
包括如下步骤：
第一步骤，基于相关计算的气溶胶光学厚度数据填补；
第二步骤，基于PM2.5浓度季节变化的遥感反演算法模型；
第三步骤，基于空间权重的PM2.5遥感监测数据的多时相插值；
第四步骤，基于台站数据与遥感数据的信息融合。
实施例
第一步骤，理解为对待填补气溶胶数据与存档数据的二维空间相关性计算，通过相关性阈值遴选合适的存档数据进行气溶胶光学厚度孔洞的填补；对大面积无值区域，采用待填补时相的台站数据与台站存档数据进行相关性计算，遴选高相关时相的气溶胶光学厚度数据，填补无数据区域。
①二维图像的孔洞像元填补
实行二维空间相关的自适应滑动窗口，在待填补像元位置(Xi，Yi)周边搜索，至少保证10×10有效像元数量，遍历历史存档的气溶胶光学厚度数据，遴选，并确定平均相关系数在0.4-1.0之间的历史存档数据，有选择地进行加权平均，赋值无数据像元；
对于小面积缺失的气溶胶光学厚度孔洞(缺失像元的搜索半径在5-20个像元以内)，选择待填补区域的单个像元P(a,b)(a，b为像元行列号)，遍历其周围最邻近非零像元，确定大于10×10的像元的行列坐标，读取像元数值X(i，j)；搜索长时间序列气溶胶光学厚度的历史存档数据，找到对应坐标的像元数值Y(i，j)，进行计算：
Figure BDA0002283416950000091
式(Ⅰ)中，R为相关系数；X(i，j)、Y(i,j)分别为待填补窗口非零像元的数值、存档气溶胶数据中相对应坐标的非零像元数值(i＝1….10,j＝1…10)，
Figure BDA0002283416950000093
分别为待填补窗口中非零像元均值、对应非零像元位置的存档数据非零像元均值。
计算得到的R值位于0.4-1.0之间，平均误差位于0-50％的存档气溶胶光学厚度数据。选出符合标准的多期气溶胶光学厚度数据。
②全无数据区域填补
台站数据与台站数据的相关性计算，确定卫星过境相同时间的站点数据，获得相关系数在0.4-1.0之间的时相，选定存档历的气溶胶光学厚度史数据，进行加权平均。
对于大面积缺失数据的区域(搜索半径距离大于20个像元)，利用卫星过境前后30-90分钟以内的n个监测台站(对于待填补的卫星过境时间，每个台站有三个监测时刻，三个监测数据)的3n个监测数据，计算其与历史存档数据的相关性和误差大小进行填补。
台站数据相关计算：
Figure BDA0002283416950000092
式(Ⅱ)中，R为相关系数；xi、yi分别为待填补时间的台站监测数据、存档台站数据中监测数值(i＝1….3n,n为监测台站数量)，
Figure BDA0002283416950000104
分别为待填补时刻的所有监测数据均值、对应时刻的存档数据均值。
计算得到的R值位于0.4-1.0之间，平均误差位于0-60％的存档数据，记录相应的成像日期，获取该时刻的对应气溶胶光学厚度数据；
③赋值计算
对多期次的气溶胶光学厚度数据，进行加权平均计算：
Figure BDA0002283416950000101
式(Ⅲ)中，像元P(a,b)为待填补像元数值；a，b为像元行列号；
Figure BDA0002283416950000102
为相关性高、平均误差小的n期气溶胶光学厚度存档数据的和；n为非零像元的期数。
该方法能准确、高效率的填补由于气象原因导致的气溶胶光学厚度数据的缺失，为大气污染的连续观测提供了重要技术方法和数据支撑。
第二步骤，理解为高空间分辨率PM2.5浓度数据的计算步骤，基于第一步骤获取的全覆盖气溶胶光学厚度，利用细颗粒物光学特性季节规律建立PM2.5大气柱浓度模型，结合大气边界层变化规律，计算大气层底部的PM2.5浓度数据。
基于AOD数据和实测PM2.5数据，建立的季节性PM2.5浓度反演模型(表1)。
表1逐月PM2.5反演模型
Figure BDA0002283416950000103
Figure BDA0002283416950000111
PM2.5与AOD在12个月份的相关系数平均为0.66，决定系数R2为0.43，N＝1009，超过统计学上95％的置信区间，相关显著。
第三步骤，理解为高时间分辨率数据特征的获取。基于第二步骤获取的PM2.5浓度数据，利用大气模式和台站风力、风向实测数据、高程数据，计算得到空间分布函数，实现时相插值，获得逐小时的PM2.5浓度数据。
通过当天的遥感监测TERRA卫星的PM2.5数据、AQUA卫星的PM2.5数据，，联合第二天的TERRA卫星、AQUA卫星监测的PM2.5浓度监测数据，即以2天之间的4幅(PM2.5-TerraAOD\PM2.5-Aqua AOD)PM2.5浓度数据为基准，利用大气模式预测的空间分布规律，构建长时间序列的时相插值方法，计算得到空间分布函数，实现时相插值，获取24幅逐小时PM2.5浓度分布数据。
时相插值模型：
Figure BDA0002283416950000112
式(Ⅳ)中，，PM2.5T时刻插值为T小时的PM2.5浓度分布；PM2.50时刻、PM2.524小时分别为起始0时刻、第24小时的终点时刻由全覆盖气溶胶光学厚度AOD数据反演的大气PM2.5浓度；PM2.5T时刻模拟值、PM2.524小时大气模拟值分别为T小时、第24小时大气模式模拟的大气污染分布数值；其中T＝1，2，3…24小时，PM2.50时刻、PM2.524小时、PM2.5T时刻模拟值、PM2.524小时大气模拟值、PM2.5T时刻插值的单位均为微克/立方米。
第四步骤，理解台站数据与遥感数据的信息融合计算。基于第三步骤获取的24幅PM2.5浓度数据，实现台站记录的24个小时的点观测数据和24幅PM2.5浓度的二维面数据，进行点面数据的信息融合，获得具有台站数据精确性和遥感高空间分辨率的PM2.5浓度数据。
利用24个小时的台站点数据和24幅逐小时PM2.5空间分布数据，实现基于点-面距离方法的信息融合，提高高空间分辨率数据的精度。台站数据与遥感数据的相关性，不仅与距离关，也与方向有关，是距离和角度的函数。因此采用遥感PM2.5浓度图像的空间分布的距离，以遥感数据浓度分布隐含的方向角度，作为权重函数：
PM2.5(融合)(i+Δi，J+Δj)＝PM2.5(卫星)(i，j)+[PM2.5(台站)-PM2.5(卫星)(i+Δi，j+Δj)]W(Δi,Δj)(Ⅴ)；
Figure BDA0002283416950000121
式(Ⅴ)中；其中Δi,Δj为距离台站位置的横向、纵向像元数；能够计算距离台站位置5km半径内的点源与面源信息的融合。
i，j分别为台站所在的像元的行、列数；PM2.5(融合)(i+Δi，J+Δj)为距离台站所在像元(i，j)的距离为(Δi，Δj)的PM2.5浓度融合结果；PM2.5(卫星)(i，j)为台站所在位置的卫星观测的PM2.5浓度数据；PM2.5(台站)为台站监测的PM2.5浓度数据；PM2.5(卫星)(i+Δi，j+Δj)为距台站所在像元(i，j)的距离为(Δi，Δj)的卫星观测PM2.5浓度数据。
显然，上述实施例仅仅是为清楚地说明所作的举例，而并非对实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。而由此所引伸出的显而易见的变化或变动仍处于本专利申请权利要求的保护范围之中。

Claims (10)

1.一种逐小时高分辨率PM2.5数据的获取方法，其特征在于，包括如下步骤：
(1)基于相关计算的气溶胶光学厚度数据填补；
(2)基于PM2.5浓度季节变化的遥感反演算法模型；
(3)基于空间权重的PM2.5遥感监测数据的多时相插值；
(4)基于台站数据与遥感数据的信息融合。
2.根据权利要求1所述的逐小时高分辨率PM2.5数据的获取方法，其特征在于，在步骤(1)中，
(1-1)对零星云团覆盖，气溶胶光学厚度有孔洞的地区，填补缺失部分的数据：对待填补气溶胶数据与存档数据的二维空间相关性计算，通过相关性阈值遴选合适的存档数据进行气溶胶光学厚度孔洞的填补；
(1-2)对于云团完全覆盖、无气溶胶光学厚度数据区域，填补缺失数据：对大面积无值区域，采用待填补时相的台站数据与台站存档数据进行相关性计算，遴选高相关时相的气溶胶光学厚度数据，填补无数据区域；
(1-3)进行加权平均，赋值无数据像元。
3.根据权利要求2所述的逐小时高分辨率PM2.5数据的获取方法，其特征在于，在步骤(1-1)中：实行二维空间相关的自适应滑动窗口，在待填补像元位置(Xi，Yi)周边搜索，至少保证10×10有效像元数量，遍历历史存档的气溶胶光学厚度数据，遴选，并确定平均相关系数在0.4-1.0之间的历史存档数据：
二维图像的孔洞像元填补：对于小面积缺失的气溶胶光学厚度孔洞，缺失像元的搜索半径在5-20个像元以内，选择待填补区域的单个像元P(a,b)，a，b为像元行列号，遍历其周围最邻近非零像元，确定大于10×10的像元的行列坐标，读取像元数值X(i，j)；搜索长时间序列气溶胶光学厚度的历史存档数据，找到对应坐标的像元数值Y(i，j)，进行计算：
Figure FDA0002283416940000021
式(Ⅰ)中，R为相关系数；X(i，j)为待填补窗口非零像元的数值(i＝1….10,j＝1…10)；Y(i,j)为存档气溶胶数据中相对应坐标的非零像元数值(i＝1….10,j＝1…10)，
Figure FDA0002283416940000023
分别为待填补窗口中非零像元均值、对应非零像元位置的存档数据非零像元均值；
计算得到的R值位于0.4-1.0之间，平均误差位于0-50％的存档气溶胶光学厚度数据；选出符合标准的多期气溶胶光学厚度数据。
4.根据权利要求2所述的逐小时高分辨率PM2.5数据的获取方法，其特征在于，在步骤(1-2)中：台站数据与台站数据的相关性计算，确定卫星过境相同时间的站点数据，获得相关系数在0.4-1.0之间的时相，选定存档的气溶胶光学厚度历史数据；完全无数据区域填补：对于大面积缺失数据的区域，搜索半径距离大于20个像元，利用卫星过境前后30-90分钟以内的n个监测台站的3n个监测数据，对于待填补的卫星过境时间，每个台站有三个监测时刻，三个监测数据；计算3n个监测数据与历史存档数据的相关性和误差大小进行填补；
台站数据相关计算：
Figure FDA0002283416940000022
式(Ⅱ)中，R为相关系数；xi、yi分别为待填补时间的台站监测数据、存档台站数据中监测数值，i＝1….3n,n为监测台站数量，
Figure FDA0002283416940000024
分别为待填补时刻的所有监测数据均值、对应时刻的存档数据均值；
计算得到的R值位于0.4-1.0之间，平均误差位于0-60％的存档数据，记录相应的成像日期，获取该时刻的对应气溶胶光学厚度数据。
5.根据权利要求2所述的逐小时高分辨率PM2.5数据的获取方法，其特征在于，在步骤(1-3)中，对步骤(1-1)和步骤(1-2)中的多期次的气溶胶光学厚度数据，进行加权平均计算：
Figure FDA0002283416940000031
式(Ⅲ)中，像元P(a,b)为待填补像元数值；a，b为像元行列号；
Figure FDA0002283416940000032
为步骤(1-1)或步骤(1-2)得到的n期气溶胶光学厚度存档数据的和；n为非零像元的期数。
6.根据权利要求1所述的逐小时高分辨率PM2.5数据的获取方法，其特征在于，在步骤(2)中，利用大气边界层变化函数，修正地面气溶胶光学厚度，构建的逐月的PM2.5多时相反演模型的方法和精度变化规律，确定的高敏感月份和算法形式；
具体高空间分辨率PM2.5浓度数据的计算步骤为：基于步骤(1)获取的全覆盖气溶胶光学厚度，利用细颗粒物光学特性季节规律建立PM2.5大气柱浓度模型，结合大气边界层变化规律，计算大气层底部的PM2.5浓度数据。
7.根据权利要求6所述的逐小时高分辨率PM2.5数据的获取方法，其特征在于，基于全覆盖气溶胶光学厚度AOD数据和实测PM2.5数据，建立的季节性PM2.5浓度反演模型为：
1月份模型μg/m3:PM2.5＝37.15+78.21×AOD；
2月份模型μg/m3：PM2.5＝23.47+63.18×AOD；
3月份模型μg/m3：PM2.5＝20.90+46.63×AOD；
4月份模型μg/m3：PM2.5＝2.90+65.93×AOD；
5月份模型μg/m3：PM2.5＝21.57+42.00×AOD；
6月份模型μg/m3：PM2.5＝14.14+34.96×AOD；
7月份模型μg/m3：PM2.5＝19.75+33.90×AOD；
8月份模型μg/m3：PM2.5＝17.10+30.66×AOD；
9月份模型μg/m3：PM2.5＝13.34+38.34×AOD；
10月份模型μg/m3：PM2.5＝22.25+46.61×AOD；
11月份模型μg/m3：PM2.5＝21.24+65.97×AOD；
12月份模型μg/m3：PM2.5＝22.12+107.88×AOD；
PM2.5与全覆盖气溶胶光学厚度AOD数据在12个月份的相关系数平均为0.66，决定系数R2为0.43，N＝1009，超过统计学上95％的置信区间，相关显著。
8.根据权利要求1所述的逐小时高分辨率PM2.5数据的获取方法，其特征在于，在步骤(3)中，高时间分辨率数据特征的获取，包括：
(3-1)基于步骤(2)获取的PM2.5浓度数据，利用大气模式模型，结合台站监测的风力、风向实测数据、高程数据，建立基于大气模式模拟的PM2.5浓度变化权重方法；
(3-2)通过当天的遥感监测TERRA卫星的PM2.5数据、AQUA卫星的PM2.5数据，联合第二天的TERRA卫星、AQUA卫星监测的PM2.5浓度监测数据，利用步骤(3-1)的大气模式模拟的PM2.5浓度权重变化，构建长时间序列的时相插值方法，计算得到空间分布函数，实现时相插值，建立高空间分辨率数据的PM2.5浓度数据，获得24幅逐小时的PM2.5浓度数据。
9.根据权利要求8所述的逐小时高分辨率PM2.5数据的获取方法，其特征在于，所述时相插值方法的模型为：
Figure FDA0002283416940000041
式(Ⅳ)中，PM2.5T时刻插值为T小时的PM2.5浓度分布；PM2.50时刻、PM2.524小时分别为起始0时刻、第24小时的终点时刻由全覆盖气溶胶光学厚度AOD数据反演的大气PM2.5浓度；PM2.5T时刻模拟值、PM2.524小时大气模拟值分别为T小时、第24小时大气模式模拟的大气污染分布数值；其中T＝1，2，3…24小时，PM2.50时刻、PM2.524小时、PM2.5T时刻模拟值、PM2.524小时大气模拟值、PM2.5T时刻插值的单位均为微克/立方米。
10.根据权利要求1所述的逐小时高分辨率PM2.5数据的获取方法，其特征在于，在步骤(4)中，基于步骤(3)获取的24幅PM2.5浓度数据，利用24个小时的台站点数据和24幅逐小时PM2.5空间分布数据，实现台站记录的24个小时的点观测数据和24幅PM2.5浓度的二维面数据，进行点-面距离方法的信息融合，实现逐小时数据的空间化过程，提高高空间分辨率数据的精度，形成每小时的PM2.5浓度分布数据；
采用遥感PM2.5浓度图像的空间分布的距离，遥感数据浓度分布隐含的方向角度，作为权重函数为：
PM2.5(融合)(i+Δi，J+Δj)＝PM2.5(卫星)(i，j)+[PM2.5(台站)-PM2.5(卫星)(i+Δi，j+Δj)]W(Δi,Δj)(Ⅴ)；
式(Ⅴ)中：Δi,Δj分别为距离台站位置的横向像元数、纵向像元数；其中w(Δi,Δj)为：
Figure FDA0002283416940000051
i，j分别为台站所在的像元的行、列数；PM2.5(融合)(i+Δi，J+Δj)为距离台站所在像元(i，j)的距离为(Δi，Δj)的PM2.5浓度融合结果；PM2.5(卫星)(i，j)为台站所在位置的卫星观测的PM2.5浓度数据；PM2.5(台站)为台站监测的PM2.5浓度数据；PM2.5(卫星)(i+Δi，j+Δj)为距台站所在像元(i，j)的距离为(Δi，Δj)的卫星观测PM2.5浓度数据；
能够计算距离台站位置5km半径内的点源与面源信息的融合。
CN201911150572.7A 2019-11-21 2019-11-21 一种逐小时高分辨率pm2.5数据的获取方法 Active CN110909309B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911150572.7A CN110909309B (zh) 2019-11-21 2019-11-21 一种逐小时高分辨率pm2.5数据的获取方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911150572.7A CN110909309B (zh) 2019-11-21 2019-11-21 一种逐小时高分辨率pm2.5数据的获取方法

Publications (2)

Publication Number Publication Date
CN110909309A true CN110909309A (zh) 2020-03-24
CN110909309B CN110909309B (zh) 2021-04-30

Family

ID=69818451

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911150572.7A Active CN110909309B (zh) 2019-11-21 2019-11-21 一种逐小时高分辨率pm2.5数据的获取方法

Country Status (1)

Country Link
CN (1) CN110909309B (zh)

Cited By (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111596093A (zh) * 2020-04-21 2020-08-28 天津大学 一种基于adcp的海水流速数据处理方法
CN111737850A (zh) * 2020-05-15 2020-10-02 中国科学院空天信息创新研究院 一种基于像元尺度上不确定性的多源卫星aod融合方法
CN112800603A (zh) * 2021-01-26 2021-05-14 北京航空航天大学 一种基于集合最优插值算法的大气环境数据同化方法
CN113297527A (zh) * 2021-06-09 2021-08-24 四川大学 基于多源城市大数据的pm2.5全面域时空计算推断方法
CN113344149A (zh) * 2021-08-06 2021-09-03 浙江大学 一种基于神经网络的pm2.5逐小时预测方法
CN113837566A (zh) * 2021-09-08 2021-12-24 广州大学 一种pm2.5人群暴露风险评估方法及装置
CN114018773A (zh) * 2021-11-03 2022-02-08 中科三清科技有限公司 Pm2.5浓度空间分布数据的获取方法、装置、设备及存储介质
CN115859026A (zh) * 2022-11-18 2023-03-28 二十一世纪空间技术应用股份有限公司 一种高分辨率近地面pm2.5浓度遥感反演方法及装置
CN112800603B (zh) * 2021-01-26 2024-05-24 北京航空航天大学 一种基于集合最优插值算法的大气环境数据同化方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20160091474A1 (en) * 2014-09-29 2016-03-31 Tanguy Griffon Method and a System for Determining at Least One Forecasted Air Quality Health Effect Caused in a Determined Geographical Area by at Least One Air Pollutant
CN106979911A (zh) * 2017-03-07 2017-07-25 南京航空航天大学 利用卫星多光谱影像数据进行pm 2.5和pm 10估算的方法
CN107202750A (zh) * 2017-05-17 2017-09-26 河北中科遥感信息技术有限公司 一种大气颗粒物星地综合监测定量遥感融合处理方法
CN110030934A (zh) * 2019-04-30 2019-07-19 中南大学 基于modis卫星传感器的气溶胶光学厚度的获取方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20160091474A1 (en) * 2014-09-29 2016-03-31 Tanguy Griffon Method and a System for Determining at Least One Forecasted Air Quality Health Effect Caused in a Determined Geographical Area by at Least One Air Pollutant
CN106979911A (zh) * 2017-03-07 2017-07-25 南京航空航天大学 利用卫星多光谱影像数据进行pm 2.5和pm 10估算的方法
CN107202750A (zh) * 2017-05-17 2017-09-26 河北中科遥感信息技术有限公司 一种大气颗粒物星地综合监测定量遥感融合处理方法
CN110030934A (zh) * 2019-04-30 2019-07-19 中南大学 基于modis卫星传感器的气溶胶光学厚度的获取方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
TAOXUE ET.AL: "Spatiotemporal continuous estimates of PM2.5 concentrations in China, 2000–2016: A machine learning method with inputs from satellites, chemical transport model, and ground observations", 《HTTPS://DOI.ORG/10.1016/J.ENVINT.2018.11.075》 *
胡泓达: ""利用气溶胶光学厚度遥感估算PM2.5浓度的时空回归克里金方法", 《中国博士学位论文全文数据库 工程科技Ⅰ辑》 *

Cited By (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111596093A (zh) * 2020-04-21 2020-08-28 天津大学 一种基于adcp的海水流速数据处理方法
CN111737850B (zh) * 2020-05-15 2024-04-19 中国科学院空天信息创新研究院 一种基于像元尺度上不确定性的多源卫星aod融合方法
CN111737850A (zh) * 2020-05-15 2020-10-02 中国科学院空天信息创新研究院 一种基于像元尺度上不确定性的多源卫星aod融合方法
CN112800603A (zh) * 2021-01-26 2021-05-14 北京航空航天大学 一种基于集合最优插值算法的大气环境数据同化方法
CN112800603B (zh) * 2021-01-26 2024-05-24 北京航空航天大学 一种基于集合最优插值算法的大气环境数据同化方法
CN113297527A (zh) * 2021-06-09 2021-08-24 四川大学 基于多源城市大数据的pm2.5全面域时空计算推断方法
CN113297527B (zh) * 2021-06-09 2022-07-26 四川大学 基于多源城市大数据的pm2.5全面域时空计算推断方法
CN113344149A (zh) * 2021-08-06 2021-09-03 浙江大学 一种基于神经网络的pm2.5逐小时预测方法
CN113837566A (zh) * 2021-09-08 2021-12-24 广州大学 一种pm2.5人群暴露风险评估方法及装置
CN114018773A (zh) * 2021-11-03 2022-02-08 中科三清科技有限公司 Pm2.5浓度空间分布数据的获取方法、装置、设备及存储介质
CN114018773B (zh) * 2021-11-03 2022-10-04 中科三清科技有限公司 Pm2.5浓度空间分布数据的获取方法、装置、设备及存储介质
CN115859026B (zh) * 2022-11-18 2023-12-05 二十一世纪空间技术应用股份有限公司 一种高分辨率近地面pm2.5浓度遥感反演方法及装置
CN115859026A (zh) * 2022-11-18 2023-03-28 二十一世纪空间技术应用股份有限公司 一种高分辨率近地面pm2.5浓度遥感反演方法及装置

Also Published As

Publication number Publication date
CN110909309B (zh) 2021-04-30

Similar Documents

Publication Publication Date Title
CN110909309B (zh) 一种逐小时高分辨率pm2.5数据的获取方法
Zhai et al. An improved geographically weighted regression model for PM2. 5 concentration estimation in large areas
WO2020043030A1 (zh) 大气污染监测设备数据可信度评价及校准方法
Zeng et al. Daily global solar radiation in China estimated from high‐density meteorological observations: a random forest model framework
Pielke Sr et al. Documentation of uncertainties and biases associated with surface temperature measurement sites for climate change assessment
Greybush et al. The ensemble Mars atmosphere reanalysis system (EMARS) version 1.0
Baker A cluster analysis of long range air transport pathways and associated pollutant concentrations within the UK
CN110428104A (zh) 一种污染贡献率确定方法、装置、电子设备及存储介质
Deeter et al. Evaluation of MOPITT retrievals of lower‐tropospheric carbon monoxide over the United States
Xu et al. Estimation of ground-level PM2. 5 concentration using MODIS AOD and corrected regression model over Beijing, China
Zhou et al. Trends in downward surface shortwave radiation from multi‐source data over China during 1984–2015
CN110411919B (zh) 一种基于卫星多光谱技术的pm2.5浓度遥感估算方法
CN114254802B (zh) 气候变化驱动下植被覆盖时空变化的预测方法
Belehaki et al. Comparison of the topside ionosphere scale height determined by topside sounders model and bottomside digisonde profiles
Lin et al. Spatio-temporal variability of aerosols over East China inferred by merged visibility-GEOS-Chem aerosol optical depth
Wei et al. Spatial interpolation of PM2. 5 concentrations during holidays in south-central China considering multiple factors
Lei et al. Full coverage estimation of the PM concentration across china based on an adaptive spatiotemporal approach
Nalini et al. High‐Resolution Lagrangian Inverse Modeling of CO2 Emissions Over the Paris Region During the First 2020 Lockdown Period
CN110310370B (zh) 一种gps与srtm点面融合的方法
Rhoads et al. Validation of land surface models using satellite‐derived surface temperature
Zeng et al. Distributed modeling of monthly air temperatures over the rugged terrain of the Yellow River Basin
CN114722350A (zh) 一种fy-3d被动微波数据云下地表温度反演与验证方法
Li et al. A deep learning approach to increase the value of satellite data for PM 2.5 monitoring in China
Wang et al. Bias correction and variability attribution analysis of surface solar radiation from MERRA-2 reanalysis
Song et al. Analysis of spatiotemporal PM2. 5 concentration patterns in Changwon, Korea, using low-cost PM2. 5 sensors

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
TA01 Transfer of patent application right

Effective date of registration: 20210409

Address after: No. 9 Dengzhuang South Road, Haidian District, Beijing 100094

Applicant after: Research Institute of aerospace information innovation, Chinese Academy of Sciences

Address before: 100101 Beijing Chaoyang District Andingmen Datun Road No. 20 North

Applicant before: Institute of Remote Sensing and Digital Earth, Chinese Academy of Sciences

TA01 Transfer of patent application right
GR01 Patent grant
GR01 Patent grant