CN110726653A - 基于天地一体化信息的pm2.5浓度监测方法 - Google Patents
基于天地一体化信息的pm2.5浓度监测方法 Download PDFInfo
- Publication number
- CN110726653A CN110726653A CN201910910189.0A CN201910910189A CN110726653A CN 110726653 A CN110726653 A CN 110726653A CN 201910910189 A CN201910910189 A CN 201910910189A CN 110726653 A CN110726653 A CN 110726653A
- Authority
- CN
- China
- Prior art keywords
- data
- atmospheric
- satellite
- reflectivity
- concentration
- 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
Links
- 238000012544 monitoring process Methods 0.000 title claims abstract description 40
- 238000000034 method Methods 0.000 title claims abstract description 32
- 230000010354 integration Effects 0.000 title claims abstract description 15
- 239000000443 aerosol Substances 0.000 claims abstract description 65
- 230000003287 optical effect Effects 0.000 claims abstract description 42
- 238000007781 pre-processing Methods 0.000 claims abstract description 29
- 239000005427 atmospheric aerosol Substances 0.000 claims abstract description 28
- 238000004364 calculation method Methods 0.000 claims abstract description 19
- 238000010276 construction Methods 0.000 claims abstract description 7
- 238000002310 reflectometry Methods 0.000 claims description 96
- 238000012545 processing Methods 0.000 claims description 30
- 238000009826 distribution Methods 0.000 claims description 29
- 230000005540 biological transmission Effects 0.000 claims description 28
- 230000005855 radiation Effects 0.000 claims description 23
- 238000012417 linear regression Methods 0.000 claims description 13
- 238000004519 manufacturing process Methods 0.000 claims description 11
- 230000008033 biological extinction Effects 0.000 claims description 9
- 238000004891 communication Methods 0.000 claims description 8
- 238000012937 correction Methods 0.000 claims description 7
- 239000013618 particulate matter Substances 0.000 claims description 6
- 239000004576 sand Substances 0.000 claims description 6
- 230000001419 dependent effect Effects 0.000 claims description 5
- 230000003595 spectral effect Effects 0.000 claims description 5
- 238000002834 transmittance Methods 0.000 claims description 5
- 238000010521 absorption reaction Methods 0.000 claims description 4
- 150000001875 compounds Chemical class 0.000 claims description 3
- 238000013500 data storage Methods 0.000 claims description 3
- 238000000354 decomposition reaction Methods 0.000 claims description 3
- 238000000605 extraction Methods 0.000 claims description 3
- 239000003102 growth factor Substances 0.000 claims description 3
- 210000005104 human peripheral blood lymphocyte Anatomy 0.000 claims description 3
- 239000000523 sample Substances 0.000 claims description 3
- 238000000926 separation method Methods 0.000 claims description 3
- 239000000853 adhesive Substances 0.000 claims 1
- 230000001070 adhesive effect Effects 0.000 claims 1
- 239000002245 particle Substances 0.000 description 7
- 230000008569 process Effects 0.000 description 6
- 238000004458 analytical method Methods 0.000 description 5
- 230000008901 benefit Effects 0.000 description 4
- 238000009877 rendering Methods 0.000 description 4
- OKTJSMMVPCPJKN-UHFFFAOYSA-N Carbon Chemical compound [C] OKTJSMMVPCPJKN-UHFFFAOYSA-N 0.000 description 2
- 229910052799 carbon Inorganic materials 0.000 description 2
- 230000008859 change Effects 0.000 description 2
- 238000009795 derivation Methods 0.000 description 2
- 230000007613 environmental effect Effects 0.000 description 2
- 238000010191 image analysis Methods 0.000 description 2
- 238000003384 imaging method Methods 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 238000001228 spectrum Methods 0.000 description 2
- 238000010561 standard procedure Methods 0.000 description 2
- 229910002651 NO3 Inorganic materials 0.000 description 1
- NHNBFGGVMKEFGY-UHFFFAOYSA-N Nitrate Chemical compound [O-][N+]([O-])=O NHNBFGGVMKEFGY-UHFFFAOYSA-N 0.000 description 1
- QAOWNCQODCNURD-UHFFFAOYSA-L Sulfate Chemical compound [O-]S([O-])(=O)=O QAOWNCQODCNURD-UHFFFAOYSA-L 0.000 description 1
- 102100029469 WD repeat and HMG-box DNA-binding protein 1 Human genes 0.000 description 1
- 101710097421 WD repeat and HMG-box DNA-binding protein 1 Proteins 0.000 description 1
- 238000003915 air pollution Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 210000000621 bronchi Anatomy 0.000 description 1
- 238000011109 contamination Methods 0.000 description 1
- 125000004122 cyclic group Chemical group 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 239000000428 dust Substances 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 239000010419 fine particle Substances 0.000 description 1
- 239000003292 glue Substances 0.000 description 1
- 230000036541 health Effects 0.000 description 1
- 210000004072 lung Anatomy 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 210000002345 respiratory system Anatomy 0.000 description 1
- 230000004044 response Effects 0.000 description 1
- 238000005096 rolling process Methods 0.000 description 1
- 238000007619 statistical method Methods 0.000 description 1
- 239000000126 substance Substances 0.000 description 1
- 231100000331 toxic Toxicity 0.000 description 1
- 230000002588 toxic effect Effects 0.000 description 1
Images
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
- 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
- Chemical & Material Sciences (AREA)
- Dispersion Chemistry (AREA)
- Physics & Mathematics (AREA)
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Analytical Chemistry (AREA)
- Biochemistry (AREA)
- General Health & Medical Sciences (AREA)
- General Physics & Mathematics (AREA)
- Immunology (AREA)
- Pathology (AREA)
- Image Processing (AREA)
Abstract
本发明公开了一种基于天地一体化信息的PM2.5浓度监测方法，包括以下步骤：A：从极轨卫星接收站实时接收并获取遥感卫星影像数据；B：预处理模块将遥感卫星影像数据转化为各级影像产品，计算最终得到亮度温度值和地表反射率值；C：进行大气气溶胶光学厚度的反演，得到已确定的气溶胶模式下的大气气溶胶光学厚度AOD；D：构建卫星遥感监测PM2.5浓度的空间变化系数回归模型；得到待估算点处的PM2.5浓度计算公式；E：计算得出待监测区域影像中所有的像素点处的PM2.5浓度yi。本发明能够高效准确的实现指定监测区域内的PM2.5浓度监测。
Description
技术领域
本发明涉及一种PM2.5浓度监测方法，尤其涉及一种基于天地一体化信息的PM2.5浓度监测方法。
背景技术
大气细颗粒物(PM2.5)指空气动力学直径小于等于2.5μm的颗粒，它的化学成分主要有有机碳、元素碳、硝酸盐和硫酸盐等，是大气污染的重要成分之一，也是造成灰霾的主要原因。大气细颗粒物极易附带有毒有害物质，长时间在空间中悬浮会被人体呼吸进入呼吸道、支气管和肺部，对人体健康造成很大影响。
国家和各级环保部门近年来采取多种措施改善空气质量，对监测***也尤为重视。传统地面监测手段存在观测点离散、分布不均匀和建造成本高等缺点，与传统的地面监测手段相比，极轨卫星遥感由于具有大范围空间尺度、连续对地观测成像和高时间分辨率过境等优势，为近地面PM2.5浓度信息获取提供了多种可能手段。
目前，利用遥感的手段，***化全自动监测PM2.5颗粒物浓度的方法尚未普及，相关环保部门缺少一种基于天地一体化信息的PM2.5颗粒物浓度自动化流程监测方法来支持雾霾治理工作。
发明内容
本发明的目的是提供一种基于天地一体化信息的PM2.5浓度监测方法，能够高效准确的实现指定监测区域内的PM2.5浓度监测。
本发明采用下述技术方案：
一种基于天地一体化信息的PM2.5浓度监测方法，包括以下步骤：
A：从极轨卫星接收站实时接收并获取遥感卫星影像数据；
B：预处理模块对获取的遥感卫星影像数据进行预处理，将遥感卫星影像数据转化为各级影像产品，然后通过几何校正和辐亮度计算最终得到亮度温度值和地表反射率值；
C：进行大气气溶胶光学厚度的反演，即首先进行云识别并剔除干扰像元的干扰，再确定气溶胶模式，利用6S大气辐射传输模型构建查找表，然后通过步骤B中获得的各级影像产品结合查找表，反演得到已确定的气溶胶模式下的大气气溶胶光学厚度；
D：构建卫星遥感监测PM2.5浓度的空间变化系数回归模型；得到待估算点i处的PM2.5浓度yi为：
i＝1,2,…,n；
其中，(ui,vi)为待估算点i的经纬度坐标；β0(ui,vi)为回归常数；βd(ui,vi)表示为观测点(ui,vi)处的第d个回归系数，d＝1，2，3…；βd(ui,vi)是地理位置的函数，xid代表待估算点i的第d个与PM2.5相关的因子，xid包括近地面消光系数、边界层高度、相对湿度、近地面温度和风速；随机误差项服从正态分布εi～N(0,σ2)，随机误差相互独立，不存在序列相关Cov(εi,εj)＝0(i≠j，(i，j＝1,2,…,n))；
E：利用步骤D中得到的待估算点i处的PM2.5浓度yi的计算公式，计算得出待监测区域影像中所有的像素点i处的PM2.5浓度yi，即待监测区域PM2.5浓度分布的反演结果。
所述的步骤A中，从极轨卫星接收站实时接收并获取Terra和Aqua的MODIS载荷多颗遥感卫星影像数据。
所述的步骤A包括如下具体步骤：
A1：接收处理***进行卫星轨道两行数据下载，完成两行数据格式化分解后进行两行数据解算，生成轨道预报文件；
A2：根据轨道预报文件中的卫星轨道数据，接收处理***进行卫星轨道关联和冲突预设消解工作，最终生成跟踪接收计划；
A3：接收处理***根据生成的跟踪接收计划，进行时间符合，在时间符合成功后接收处理***对运行调度控制指令进行分析，将运行调度控制指令发送至天线设备，完成卫星的捕获与跟踪；同时，接收处理***根据预先设置的卫星数据接收设备参数，完成卫星数据锁定与接收；
A4：卫星的数传频段为X波段，接收处理***采用26.25Mbs-1和15.0Mbs-1的码速率分别接收Terra和Aqua卫星的X频段信号，经过下变频器处理输出得到中频信号，再由接收机解调后输出数字基带信号；
A5：根据不同卫星数传数据格式要求，接收处理***通过对数字基带信号中帧同步字段的搜索，完成数据帧的字节对齐、容错和填充；
A6：接收处理***根据卫星数传数据加扰方式，按照多项式生成PN码与数据帧同步字段后面的图像数据进行异或运算进行解扰；再根据卫星数传数据编码方式，对经过解扰后的图像数据进行RS译码以及CRC校验，最终得到影像的原始码流数据，即遥感卫星影像数据。
所述的步骤B包括如下具体步骤：
B1：预处理模块实时监控原始码流数据存储路径，当监测到有新的原始码流数据文件产生，且满足新的原始码流数据文件不被占用同时大小大于0KB的条件时，开始对新的原始码流数据文件进行预处理；
B2：预处理模块通过解析新的原始码流数据文件的文件名，根据不同卫星载荷名称自动加载不同卫星载荷相应的xml配置文件；
B3：预处理模块对新的原始码流数据文件按照不同卫星载荷的数据帧格式，进行高级在轨***帧格式解析，即自动调用新的原始码流数据文件的标准帧并解码，解码包括虚拟通道数据单元分离以及从不同虚拟通道数据单元提取不同信源包，重新定格式后输出虚拟通道数据单元，生成L0级数据产品；
其中，L0级数据产品，指经数据重构，且未进行任何处理的原始数据；L0级数据产品中所有的通信信息被移除，通信信息包括数据帧同步字段、通信头和重复数据；
B4：预处理模块结合辅助数据信息，从L0级数据产品中的CCSDS包解包还原出扫描数据的集合，即预处理模块通过提取辅助数据信息中的星历、姿态、跳秒和世界标准时间信息，对图像进行校正，完成L1A级数据产品的生产；
其中，辅助数据信息包括星历、姿态、跳秒、世界标准时间、辐射定标参数、反射率参数和质量参数信息；L1A级产品，指经数据重构，且具有时间参考及地理坐标参数信息的未进行任何处理的原始数据；
B5：在L1A级数据产品的基础上，预处理模块利用辅助数据信息中的辐射定标参数、反射率参数和质量参数信息，结合L1A级产品中的地理坐标参数信息，根据不同卫星载荷进行定位和定标运算，生成包括经纬度信息、探测器高度和不同通道的大气表观反射率与辐射率的数据集，即L1B级数据产品；
其中，L1B级数据产品，指在L1A级数据产品的基础上，进行定位和定标处理后的数据集；
B6：预处理模块自动调取包含大气表观反射率、辐射率和经纬度信息的L1B级数据产品，通过几何校正和辐亮度计算最终得到亮度温度值和地表反射率值。
所述的步骤C中，气溶胶模式为大陆型气溶胶模式。
所述的步骤C中包括如下具体步骤：
C1：云识别模块从步骤B6中得到的不同通道的大气表观反射率值和亮度温度值，设定相应的阈值分别对干扰像元进行识别并剔除干扰像元的干扰，干扰像元包括云、雾和高反射率地物；
其中，使用的MODIS的通道有1、2、3、18、20、22、26、31和35，ρR为通道1可见光红光波段的大气表观反射率，ρB为通道3可见光蓝光波段的大气表观反射率，ρNR2为通道2近红外波段的大气表观反射率，ρNR18为通道18近红外波段的大气表观反射率，ρNR26为通道26近红外波段的大气表观反射率，BT20、BT22、BT31和BT35分别为通道20、22、31和35热红外波段的亮度温度；
若通道1可见光红光波段的大气表观反射率ρR大于0.25，则判定该区域为云、雾和高反射率地物；若通道1可见光红光波段的大气表观反射率ρR小于等于0.25时，且通道3可见光蓝光波段的大气表观反射率ρB大于0.28，则判定该区域为云、雾和高反射率地物；
若通道35热红外波段的亮度温度BT35不小于230且不大于250时，则判定该区域为云；
若通道2近红外波段与通道1可见光红光波段的大气表观反射率比值ρNR2/ρR为[0.9，1.6]时，或通道18近红外波段与通道2近红外波段的大气表观反射率比值ρNR18/ρNR2为[0.28，0.5]时，则判定该区域为云；
若通道26大气表观反射率ρNR26大于0.05时，则判定该区域为云；
若通道20热红外波段的亮度温度BT20与通道22热红外波段的亮度温度BT22的差值为[2.8，7.6]，或通道20热红外波段的亮度温度BT20与通道31热红外波段的亮度温度BT31的差值为[10，19]，则判定该区域为云和雾；
C2：确定气溶胶模式为大陆型并构造查找表；
将确定的大陆型气溶胶模式、气溶胶含量参数、卫星观测几何参数、波段的中心波长和探测器高度带入6S大气辐射传输模型中，计算得到不同波段及不同气溶胶光学厚度下的大气参数，从而构建并得到查找表；接收数据大气模式、光谱参数和地表反射率类型从6S辐射传输模型参数设置选项中选择；
其中，气溶胶模式为大陆型气溶胶模式，气溶胶含量参数设定为6个在波长0.55μm处的大气气溶胶光学厚度值：即0、0.25、0.5、1.0、1.5和1.95；卫星观测几何参数设定选取0°、6°、12°、24°、35.2°、48°、54°、60°和66°共9个太阳天顶角，从[0°，66°]范围内每隔6°的12个观测天顶角，从[0°，180°]范围内每隔12°的16个相对方位角；波段的中心波长取0.47μm和0.66μm；探测器高度从步骤B5中得到的L1B数据产品中获得；
C3：根据大气辐射传输原理，当采用标量法进行陆地气溶胶反演时，在地表朗伯体和大气水平均一假设条件下，大气表观反射率ρTOA表示为：
式中，
θs与θv分别为太阳天顶角与观测天顶角，φ为相对方位角，T为大气透过率；S为大气下界的半球反射率，ρ0是大气的路径辐射项等效反射率，ρs为地表反射率；下角标s表示太阳，v表示观测；
从步骤C2中得到的查找表中选取与卫星观测几何最邻近的大气参数，根据L1B级数据产品进行线性插值，得到不同波段和不同气溶胶光学厚度下的大气参数中大气下界的半球反射率S、大气的路径辐射项等效反射率ρ0和大气透过率T的数值，从MODIS多天历史的地表反射率数据中查找晴天状态下的地表反射率数据，构造晴天地表反射率库，得到所设定的时间段内的晴天状态下的地表反射率ρs，通过式(1)采用最小二乘法拟合模拟真实的大气表观反射率ρTOA，对真实的大气表观进行线性插值，最终获得大气气溶胶光学厚度。
所述的步骤D中，采用的颗粒物反演模型基于地理加权回归的PM2.5反演模型，构建原理基于近地面PM2.5浓度和大气气溶胶光学厚度AOD之间存在的先验关系，即首先根据近地面PM2.5浓度和大气气溶胶光学厚度AOD之间的先验关系推导得到PM2.5的多元线性回归模型的对数形式；
近地面PM2.5浓度和大气气溶胶光学厚度AOD之间的先验关系为：
其中，ρ为PM2.5的平均密度，r为有效半径，Q为消光效率因子，AOD为大气气溶胶光学厚度，HPBL为边界层高度，f(RH)为气溶胶散射吸湿增长因子，一般表示为相对湿度的函数f(RH)＝1/(1.0-RH/100)；对上式两边同时取自然对数，得到PM2.5的多元线性回归模型的对数形式为：
其中，β1，β2，β3为回归系数，β1＝1，β2＝-1，β3＝-1；β0为回归常数；
将数据的地理位置，嵌入到PM2.5的多元线性回归模型的对数形式中的回归系数之中，得到卫星遥感监测PM2.5浓度的空间变化系数回归模型；
卫星遥感监测PM2.5浓度的空间变化系数回归模型为：
待估算点i处的PM2.5浓度yi为：
其中，(ui,vi)为待估算点i的经纬度坐标；β0(ui，vi)为回归常数；βd(ui,vi)表示为观测点(ui,vi)处的第d个回归系数，d＝1，2，3…；βd(ui,vi)是地理位置的函数，xid代表待估算点i的第d个与PM2.5相关的因子，xid包括近地面消光系数、边界层高度、相对湿度、近地面温度和风速；随机误差项服从正态分布εi～N(0,σ2)，随机误差相互独立，不存在序列相关Cov(εi,εj)＝0(i≠j，(i，j＝1，2，…,n))。
还包括步骤F；
F：利用步骤E中获得的待监测区域PM2.5浓度分布的反演结果，制作待监测区域PM2.5浓度分布专题图。
本发明通过获取遥感卫星影像数据，并将遥感卫星影像数据转化为各级影像产品，通过几何校正和辐亮度计算得到亮度温度值和地表反射率值，然后利用得到的亮度温度值和地表反射率值进行大气气溶胶光学厚度的反演，得到大气气溶胶光学厚度AOD；随后本发明通过构建卫星遥感监测PM2.5浓度的空间变化系数回归模型，计算得到待估算点i处的PM2.5浓度yi，并最终得到待监测区域PM2.5浓度分布的反演结果，实现基于天地一体化信息的PM2.5浓度监测。上述过程中，本发明基于影像光谱波段的不同反射特征，设定不同的阈值通过云识别模块把影像中的干扰像元检测出来，有效地降低了云、雾和高反射率地物对气溶胶光学厚度反演精度的影响，同时，在卫星遥感监测PM2.5浓度的空间变化系数回归模型中，将数据的地理位置嵌入到PM2.5的多元线性回归模型的对数形式中的回归系数之中，使得回归系数随地理空间位置而变化，避免因回归系数单一而导致反演结果产生较大误差，也避免了现有线性回归模型中无法引入不同区域的气象参数和地表类型的缺陷，进一步提高PM2.5浓度监测的准确度。
附图说明
图1为本发明的流程示意图；
图2为本发明中云识别模块的识别过程示意图。
具体实施方式
以下结合附图和实施例对本发明作以详细的描述：
如图1所示，本发明所述的基于天地一体化信息的PM2.5浓度监测方法，包括以下步骤：
A：从极轨卫星接收站实时接收并获取遥感卫星影像数据；
本发明中，从极轨卫星接收站实时接收并获取Terra和Aqua的MODIS载荷多颗遥感卫星影像数据。Terra和Aqua是美国地球观测***计划中两颗重要卫星，Terra为上午星，Aqua为下午星，MODIS是搭载在Terra和Aqua卫星上的中分辨率成像光谱仪。
本实施例中，步骤A包括如下具体步骤：
A1：接收处理***进行卫星轨道两行数据下载，完成两行数据格式化分解后进行两行数据解算，生成轨道预报文件。
A2：根据轨道预报文件中的卫星轨道数据，接收处理***进行卫星轨道关联和冲突预设消解工作，最终生成跟踪接收计划。
A3：接收处理***根据生成的跟踪接收计划，进行时间符合，在时间符合成功后接收处理***对运行调度控制指令进行分析，将运行调度控制指令发送至天线设备，完成卫星的捕获与跟踪；同时，接收处理***根据预先设置的卫星数据接收设备参数，完成卫星数据锁定与接收。
A4：卫星的数传频段为X波段，接收处理***采用26.25Mbs-1和15.0Mbs-1的码速率分别接收Terra和Aqua卫星的X频段信号，经过下变频器处理输出得到中频信号，再由接收机解调后输出数字基带信号。
A5：根据不同卫星数传数据格式要求，接收处理***通过对数字基带信号中帧同步字段的搜索，完成数据帧的字节对齐、容错和填充。
A6：接收处理***根据卫星数传数据加扰方式，按照多项式生成PN码(Pseudor-Noise伪随机码)与数据帧同步字段后面的图像数据进行异或运算进行解扰；再根据卫星数传数据编码方式，对经过解扰后的图像数据进行RS(Reed-Solomon)译码以及CRC校验，最终得到影像的原始码流数据，即遥感卫星影像数据。
本实施例中，操作员可以对上述整个跟踪接收计划的生成和下发的工作状态，以及卫星跟踪和接收的设备控制响应情况进行实时监视。
B：预处理模块对获取的遥感卫星影像数据进行预处理，将遥感卫星影像数据转化为各级影像产品，然后通过几何校正和辐亮度计算最终得到亮度温度值和地表反射率值；
本发明中，步骤B包括如下具体步骤：
B1：预处理模块实时监控原始码流数据存储路径，当监测到有新的原始码流数据文件产生，且满足新的原始码流数据文件不被占用同时大小大于0KB的条件时，开始对新的原始码流数据文件进行预处理；
B2：预处理模块通过解析新的原始码流数据文件的文件名，根据不同卫星载荷名称自动加载不同卫星载荷相应的xml配置文件；
B3：预处理模块对新的原始码流数据文件按照不同卫星载荷的数据帧格式，进行高级在轨***帧格式解析，即自动调用新的原始码流数据文件的标准帧并解码，解码包括虚拟通道数据单元分离以及从不同虚拟通道数据单元提取不同信源包，重新定格式后输出虚拟通道数据单元，生成L0级数据产品；L0级数据产品，指经数据重构，且未进行任何处理的原始数据；L0级数据产品中所有的通信信息被移除，通信信息包括数据帧同步字段、通信头和重复数据。
B4：预处理模块结合辅助数据信息，从L0级数据产品中的CCSDS包解包还原出扫描数据的集合，即预处理模块通过提取辅助数据信息中的星历、姿态、跳秒和世界标准时间信息，对图像进行校正，完成L1A级数据产品的生产；其中，辅助数据信息包括星历、姿态、跳秒、世界标准时间、辐射定标参数、反射率参数和质量参数信息；L1A级产品，指经数据重构，且具有时间参考及地理坐标参数信息的未进行任何处理的原始数据。
B5：在L1A级数据产品的基础上，预处理模块利用辅助数据信息中的辐射定标参数、反射率参数和质量参数信息，结合L1A级产品中的地理坐标参数信息，根据不同卫星载荷进行定位和定标运算，生成包括经纬度信息、探测器高度和不同通道的大气表观反射率与辐射率的数据集，即L1B级数据产品。L1B级数据产品，指在L1A级数据产品的基础上，进行定位和定标处理后的数据集。
B6：预处理模块自动调取包含大气表观反射率、辐射率和经纬度信息的L1B级数据产品，通过几何校正和辐亮度计算最终得到亮度温度值和地表反射率值。
要进行大气气溶胶光学厚度的反演，必须分离地表反射率的贡献，通过对云像元进行识别，剔除云、雾和高反射率地物的干扰，然后建立气溶胶模式，构建查找表，从而通过经步骤B预处理后得到的L1B级数据产品反演某种气溶胶模式下的气溶胶光学厚度，其中，高反射率地物指反射率不小于50％的非选择性反射体。因此，本申请中，步骤C如下所述：
C：进行大气气溶胶光学厚度的反演，即首先进行云识别并剔除干扰像元的干扰，再确定气溶胶模式，利用6S大气辐射传输模型构建查找表，然后通过步骤B中获得的各级影像产品结合查找表，反演得到已确定气溶胶模式下的大气气溶胶光学厚度；
大气气溶胶光学厚度，英文名称为AOD(Aerosol Optical Depth)或AOT(AerosolOptical Thickness)，定义为介质的消光系数在垂直方向上的积分，是描述气溶胶对光的削减作用的。大气气溶胶光学厚度是气溶胶最重要的参数之一，表征大气浑浊程度的关键的物理量，也是确定气溶胶气候效应的重要因素。
本申请中，云识别模块是通过基于影像光谱波段的不同反射特征，设定不同的阈值把影像中的云检测出来，剔除云、雾和高反射率地物对气溶胶光学厚度反演的影响。
气溶胶模式指一组定义好的气溶胶类型的集合，如标准气溶胶模型中的大陆型、海洋型和城市型，每种气溶胶模式是由四种不同的粒子组分按一定的比率组合而成，详见表1。所谓气溶胶类型不同，也就是气溶胶组分或其组合比率不同。实际反演中是用辐射传输模型根据假设的气溶胶模式构建查找表，由于不同大气状况下粒子组分的不同，对最终的反演结果会带来不同的影响，因此气溶胶模式是影响反演精度的一个重要因素。参考研究区的气溶胶吸湿增长和闭合实验等研究结果，我们设定待监测区域的气溶胶模式为大陆型气溶胶模式。
表1城市型、大陆型和海洋型气溶胶模式不同粒子的组成对比
确定气溶胶模式后，首先需要进行查找表的构造。查找表是通过设定不同卫星观测几何参数、大气模式、气溶胶含量参数、探测器高度、光谱参数、地表反射率类型以及观测数据所在的波段中心波长，采用某种辐射传输方程来进行。本发明使用6S辐射传输模型进行辐射传输计算来构造查找表，6S辐射传输模型采用了最新的散射计算方法，使太阳光谱波段的散射计算精度比5S有明显提高。
本发明中，步骤C中包括如下具体步骤：
C1：云识别模块从步骤B6中得到的不同通道的大气表观反射率值和亮度温度值，设定相应的阈值分别对干扰像元进行识别并剔除干扰像元的干扰，干扰像元包括云、雾和高反射率地物。根据MODIS载荷的特点及中国华北地区多影像分析，通过可见光、近红外和红外的多个波段的大气表观反射率、亮度温度和空间变化特征对阈值进行调整和优化，最终本申请中所使用的云识别模块的识别过程如图2所示。
其中，使用的MODIS的通道有1、2、3、18、20、22、26、31和35。ρR为通道1可见光红光波段的大气表观反射率，ρB为通道3可见光蓝光波段的大气表观反射率，ρNR2为通道2近红外波段的大气表观反射率，ρNR18为通道18近红外波段的大气表观反射率，ρNR26为通道26近红外波段的大气表观反射率，BT20、BT22、BT31和BT35分别为通道20、22、31和35热红外波段的亮度温度。
为了进一步提高云识别模块识别的准确性，尽可能地降低云、雾和高反射率地物对气溶胶光学厚度反演精度的影响，本申请根据不同通道下不同波段的可见光、近红外和红外的特性，对于云、雾和高反射率地物设定不同的大气表观反射率阈值进行判断。
本申请利用MODIS多时相影像分析，若通道1可见光红光波段的大气表观反射率ρR大于0.25，则判定该区域为云、雾和高反射率地物；若通道1可见光红光波段的大气表观反射率ρR小于等于0.25时，且通道3可见光蓝光波段的大气表观反射率ρB大于0.28，则判定该区域为云、雾和高反射率地物；上述阈值能够有效判定干扰像元中的厚云、浓雾和高反射率地物的情况；
若通道35热红外波段的亮度温度BT35不小于230且不大于250时，则判定该区域为云；上述阈值设定中，为了排除云识别过程时霾对于云的干扰，经分析得出，在陆地上空云区域具有较低的亮度温度值，而霾区域的亮度温度值相对高于云区域的亮度温度值，因此将通道35热红外波段的亮度温度BT35的阈值设定为[230，250]，以准确区分霾区域和云区域；
若通道2近红外波段与通道1可见光红光波段的大气表观反射率比值ρNR2/ρR为[0.9，1.6]时，或通道18近红外波段与通道2近红外波段的大气表观反射率比值ρNR18/ρNR2为[0.28，0.5]时，则判定该区域为云；经分析得出，通道1、通道2和通道18对于确定薄云、碎云和层积云具有明显优势，因此本申请中设定ρNR2/ρR的阈值为[0.9，1.6]，ρNR18/ρNR2的阈值为[0.28，0.5]，以有效判定干扰像元中存在的薄云、碎云和层积云；
若通道26大气表观反射率ρNR26大于0.05时，则判定该区域为云；经分析得出，通道26对于确定卷云具有明显优势，因此本申请中设定ρNR26的阈值为0.05，以有效判定干扰像元中存在的卷云；
若通道20热红外波段的亮度温度BT20与通道22热红外波段的亮度温度BT22的差值为[2.8，7.6]，或通道20热红外波段的亮度温度BT20与通道31热红外波段的亮度温度BT31的差值为[10，19]，则判定该区域为云和雾；经分析得出，利用通道20和通道22以及通道20和通道31的热红外波段的亮度温度的差值进行识别，对于确定低层云和雾粒子具有明显优势，因此本身在设定BT20与BT22的差值阈值为[2.8，7.6]或BT20与BT31的差值阈值为[10，19]，以有效判定干扰像元中存在的低层云和雾粒子。
上述过程中，厚云、浓雾、薄云、碎云、层积云、卷云、低层云和雾粒子均为气象学领域的专业名词。
C2：确定气溶胶模式为大陆型并构造查找表；
将确定的大陆型气溶胶模式、气溶胶含量参数、卫星观测几何参数、波段的中心波长和探测器高度带入6S大气辐射传输模型中，计算得到不同波段及不同气溶胶光学厚度下的大气参数，从而构建并得到查找表；接收数据大气模式、光谱参数和地表反射率类型从6S辐射传输模型参数设置选项中选择。
其中，气溶胶模式为大陆型气溶胶模式，气溶胶含量参数设定为6个在波长0.55μm处的大气气溶胶光学厚度值：即0、0.25、0.5、1.0、1.5和1.95；卫星观测几何参数设定选取0°、6°、12°、24°、35.2°、48°、54°、60°和66°共9个太阳天顶角，从[0°，66°]范围内每隔6°的12个观测天顶角，从[0°，180°]范围内每隔12°的16个相对方位角；波段的中心波长取0.47μm和0.66μm；探测器高度从步骤B5中得到的L1B数据产品中获得。
C3：根据大气辐射传输原理，当采用标量法进行陆地气溶胶反演时，在地表朗伯体和大气水平均一假设条件下，大气表观反射率ρTOA表示为：
式中，
θs与θv分别为太阳天顶角与观测天顶角，φ为相对方位角，T为大气透过率；S为大气下界的半球反射率，ρ0是大气的路径辐射项等效反射率，ρs为地表反射率；下角标s表示太阳，v表示观测；
从步骤C2中得到的查找表中选取与卫星观测几何最邻近的大气参数，根据L1B级数据产品进行线性插值，得到不同波段和不同气溶胶光学厚度下的大气参数中大气下界的半球反射率S、大气的路径辐射项等效反射率ρ0和大气透过率T的数值。从MODIS多天历史的地表反射率数据中查找晴天状态下的地表反射率数据，构造晴天地表反射率库，得到所设定的时间段内的晴天状态下的地表反射率ρs，通过式(1)采用最小二乘法拟合模拟真实的大气表观反射率ρTOA，对真实的大气表观进行线性插值，最终获得大气气溶胶光学厚度AOD。
D：构建卫星遥感监测PM2.5浓度的空间变化系数回归模型；
本发明采用的颗粒物反演模型是基于地理加权回归的PM2.5反演模型；其构建原理基于近地面PM2.5浓度和大气气溶胶光学厚度AOD之间存在一定的先验关系，即首先根据近地面PM2.5浓度和大气气溶胶光学厚度AOD之间的先验关系推导得到PM2.5的多元线性回归模型的对数形式；
近地面PM2.5浓度和大气气溶胶光学厚度AOD之间的先验关系为：
其中，ρ为PM2.5的平均密度，r为有效半径，Q为消光效率因子，AOD为大气气溶胶光学厚度，HPBL为边界层高度，f(RH)为气溶胶散射吸湿增长因子，一般表示为相对湿度的函数f(RH)＝1/(1.0-RH/100)；对上式两边同时取自然对数，得到PM2.5的多元线性回归模型的对数形式为：
其中，β1，β2，β3为回归系数，β1＝1，β2＝-1，β3＝-1；β0为回归常数；
PM2.5的多元线性回归模型的对数形式是一种全局空间回归模型，β1，β2，β3作为回归系数会随地理空间位置而变化。但现有的基于的区域AOD浓度遥感反演大多都没有考虑到局部参数的变化特征，根据一般经验β1，β2，β3分别取1，-1和-1。因此现有的回归系数单一的线性回归模型会导致反演结果有较大的误差，因此本发明中：
将数据的地理位置，嵌入到PM2.5的多元线性回归模型的对数形式中的回归系数之中，得到卫星遥感监测PM2.5浓度的空间变化系数回归模型；
卫星遥感监测PM2.5浓度的空间变化系数回归模型为：
待估算点i处的PM2.5浓度yi为：
其中，(ui，vi)为待估算点i的经纬度坐标；β0(ui，vi)为回归常数；βd(ui，vi)表示为观测点(ui，vi)处的第d个回归系数，d＝1，2，3…；βd(ui,vi)是地理位置的函数，在本发明中该系数会随着变量地理位置的不同而不同，能够有效避免反演结果出现较大的误差；xid代表待估算点i的第d个与PM2.5相关的因子，xid包括近地面消光系数、边界层高度、相对湿度、近地面温度和风速，从地面观测中得到；随机误差项服从正态分布εi～N(0，σ2)，随机误差相互独立，不存在序列相关Cov(εi，εj)＝0(i≠j，(i，j＝1，2,…,n))；
E：利用步骤D中得到的待估算点i处的PM2.5浓度yi的计算公式，计算得出待监测区域影像中所有的像素点i处的PM2.5浓度yi，即待监测区域PM2.5浓度分布的反演结果；
本发明中，还包括步骤F，根据获得的待监测区域PM2.5浓度分布的反演结果，制作待监测区域PM2.5浓度分布专题图；
F：利用步骤E中获得的待监测区域PM2.5浓度分布的反演结果，制作待监测区域PM2.5浓度分布专题图。
本发明中，当***接收到待监测区域PM2.5浓度分布的反演结果时，可自动进行专题图和专题报的生产。
专题图制作所需文件包括：待监测区域PM2.5浓度分布的反演结果、专题图中显示范围的矢量文件、渲染方式配置文件、比例尺、指北针及图例。
专题报制作所需文件包括：待监测区域PM2.5浓度分布的反演结果、专题图及专题报模板。
专题图制作的具体步骤如下：
1、加载待监测区域PM2.5浓度分布的反演结果，加载待监测区域行政区边界数据，读取渲染配置文件；待监测区域行政区边界数据采用shp格式矢量文件；
2、依据待监测区域行政区边界数据，对待监测区域PM2.5浓度分布的反演结果进行裁剪；
3、根据渲染方式配置文件，对裁剪后的待监测区域PM2.5浓度分布的反演结果进行渲染；
4、添加比例尺、指北针、专题图名称及文本标注；
5、提取专题图的元数据，注册到数据库，并输出专题图到指定路径下。
专题报制作步骤如下：
1、读取待监测区域PM2.5浓度分布的反演结果，依据PM2.5的像素值分为优(0＜像素值≤35)，良(35＜像素值≤75)，轻度污染(75＜像素值≤115)，中度污染(115＜像素值≤150)，重度污染(150＜像素值≤250)和严重污染(250＜像素值)6个污染等级；
2、计算基础数据：待监测区域各行政区的有效监测面积，PM2.5浓度均值，以及待监测区域内优、良、轻度污染、中度污染、重度污染和严重污染6个污染等级所占面积；
3、计算统计数据：依据步骤2得到的基础数据进行统计分析，统计出各行政区的最高污染等级面积以及分布情况，依据排名规则，对行政区进行排名。
4、读取专题报模板和专题图，***专题图，依据统计数据，制作统计表、生成文字段落。最后输出专题报到指定路径下，提取专题报的元数据，注册到数据库。
本发明还可以通过读取不同时间段范围的PM2.5浓度分布的反演结果，进行像素值加权平均，得到PM2.5浓度分布的日、周、月及年均值数据，按照专题图、报制作步骤生成日、周、月和年等不同时间频率的PM2.5浓度分布及专题报。
此外，本发明最后还可对反演的PM2.5精度进行分析。
本发明在制作PM2.5浓度分布监测专题图及专题报的同时，对注册到数据库中的反演结果进行精度分析。利用省/市/区域尺度范围内国控地基站点PM2.5观测值为基准值，计算得到卫星过境时间先后半小时范围中省内各省/市/区域地基站点算术平均值Sn。同时，根据国控站点经纬度，计算得出周边10公里范围内的本发明所测得的有效PM2.5算术平均值，并作为地基站点所在省/市/区域卫星观测平均值Pm。计算省/市/区域PM2.5浓度产品精度公式为：
通过该公式计算分析，可对本发明所测得的省/市/区域尺度PM2.5浓度的精度结果进行直接验证，有利于相关部门直观的掌握相关省/市/区域的PM2.5污染情况。
Claims (8)
1.一种基于天地一体化信息的PM2.5浓度监测方法，其特征在于，包括以下步骤：
A：从极轨卫星接收站实时接收并获取遥感卫星影像数据；
B：预处理模块对获取的遥感卫星影像数据进行预处理，将遥感卫星影像数据转化为各级影像产品，然后通过几何校正和辐亮度计算最终得到亮度温度值和地表反射率值；
C：进行大气气溶胶光学厚度的反演，即首先进行云识别并剔除干扰像元的干扰，再确定气溶胶模式，利用6S大气辐射传输模型构建查找表，然后通过步骤B中获得的各级影像产品结合查找表，反演得到已确定的气溶胶模式下的大气气溶胶光学厚度；
D：构建卫星遥感监测PM2.5浓度的空间变化系数回归模型；得到待估算点i处的PM2.5浓度yi为：
其中，(ui，vi)为待估算点i的经纬度坐标；β0(ui,vi)为回归常数；βd(ui，vi)表示为观测点(ui，vi)处的第d个回归系数，d＝1，2，3…；βd(ui,vi)是地理位置的函数，xid代表待估算点i的第d个与PM2.5相关的因子，xid包括近地面消光系数、边界层高度、相对湿度、近地面温度和风速；随机误差项服从正态分布εi～N(0，σ2)，随机误差相互独立，不存在序列相关Cov(εi，εj)＝0(i≠j，(i,j＝1，2，…，n))；
E：利用步骤D中得到的待估算点i处的PM2.5浓度yi的计算公式，计算得出待监测区域影像中所有的像素点i处的PM2.5浓度yi，即待监测区域PM2.5浓度分布的反演结果。
2.根据权利要求1所述的基于天地一体化信息的PM2.5浓度监测方法，其特征在于：所述的步骤A中，从极轨卫星接收站实时接收并获取Terra和Aqua的MODIS载荷多颗遥感卫星影像数据。
3.根据权利要求2所述的基于天地一体化信息的PM2.5浓度监测方法，其特征在于，所述的步骤A包括如下具体步骤：
A1：接收处理***进行卫星轨道两行数据下载，完成两行数据格式化分解后进行两行数据解算，生成轨道预报文件；
A2：根据轨道预报文件中的卫星轨道数据，接收处理***进行卫星轨道关联和冲突预设消解工作，最终生成跟踪接收计划；
A3：接收处理***根据生成的跟踪接收计划，进行时间符合，在时间符合成功后接收处理***对运行调度控制指令进行分析，将运行调度控制指令发送至天线设备，完成卫星的捕获与跟踪；同时，接收处理***根据预先设置的卫星数据接收设备参数，完成卫星数据锁定与接收；
A4：卫星的数传频段为X波段，接收处理***采用26.25Mbs-1和15.0Mbs-1的码速率分别接收Terra和Aqua卫星的X频段信号，经过下变频器处理输出得到中频信号，再由接收机解调后输出数字基带信号；
A5：根据不同卫星数传数据格式要求，接收处理***通过对数字基带信号中帧同步字段的搜索，完成数据帧的字节对齐、容错和填充；
A6：接收处理***根据卫星数传数据加扰方式，按照多项式生成PN码与数据帧同步字段后面的图像数据进行异或运算进行解扰；再根据卫星数传数据编码方式，对经过解扰后的图像数据进行RS译码以及CRC校验，最终得到影像的原始码流数据，即遥感卫星影像数据。
4.根据权利要求3所述的基于天地一体化信息的PM2.5浓度监测方法，其特征在于，所述的步骤B包括如下具体步骤：
B1：预处理模块实时监控原始码流数据存储路径，当监测到有新的原始码流数据文件产生，且满足新的原始码流数据文件不被占用同时大小大于0KB的条件时，开始对新的原始码流数据文件进行预处理；
B2：预处理模块通过解析新的原始码流数据文件的文件名，根据不同卫星载荷名称自动加载不同卫星载荷相应的xml配置文件；
B3：预处理模块对新的原始码流数据文件按照不同卫星载荷的数据帧格式，进行高级在轨***帧格式解析，即自动调用新的原始码流数据文件的标准帧并解码，解码包括虚拟通道数据单元分离以及从不同虚拟通道数据单元提取不同信源包，重新定格式后输出虚拟通道数据单元，生成L0级数据产品；
其中，L0级数据产品，指经数据重构，且未进行任何处理的原始数据；L0级数据产品中所有的通信信息被移除，通信信息包括数据帧同步字段、通信头和重复数据；
B4：预处理模块结合辅助数据信息，从L0级数据产品中的CCSDS包解包还原出扫描数据的集合，即预处理模块通过提取辅助数据信息中的星历、姿态、跳秒和世界标准时间信息，对图像进行校正，完成L1A级数据产品的生产；
其中，辅助数据信息包括星历、姿态、跳秒、世界标准时间、辐射定标参数、反射率参数和质量参数信息；L1A级产品，指经数据重构，且具有时间参考及地理坐标参数信息的未进行任何处理的原始数据；
B5：在L1A级数据产品的基础上，预处理模块利用辅助数据信息中的辐射定标参数、反射率参数和质量参数信息，结合L1A级产品中的地理坐标参数信息，根据不同卫星载荷进行定位和定标运算，生成包括经纬度信息、探测器高度和不同通道的大气表观反射率与辐射率的数据集，即L1B级数据产品；
其中，L1B级数据产品，指在L1A级数据产品的基础上，进行定位和定标处理后的数据集；
B6：预处理模块自动调取包含大气表观反射率、辐射率和经纬度信息的L1B级数据产品，通过几何校正和辐亮度计算最终得到亮度温度值和地表反射率值。
5.根据权利要求1所述的基于天地一体化信息的PM2.5浓度监测方法，其特征在于：所述的步骤C中，气溶胶模式为大陆型气溶胶模式。
6.根据权利要求4所述的基于天地一体化信息的PM2.5浓度监测方法，其特征在于：所述的步骤C中包括如下具体步骤：
C1：云识别模块从步骤B6中得到的不同通道的大气表观反射率值和亮度温度值，设定相应的阈值分别对干扰像元进行识别并剔除干扰像元的干扰，干扰像元包括云、雾和高反射率地物；
其中，使用的MODIS的通道有1、2、3、18、20、22、26、31和35，ρR为通道1可见光红光波段的大气表观反射率，ρB为通道3可见光蓝光波段的大气表观反射率，ρNR2为通道2近红外波段的大气表观反射率，ρNR18为通道18近红外波段的大气表观反射率，ρNR26为通道26近红外波段的大气表观反射率，BT20、BT22、BT31和BT35分别为通道20、22、31和35热红外波段的亮度温度；
若通道1可见光红光波段的大气表观反射率ρR大于0.25，则判定该区域为云、雾和高反射率地物；若通道1可见光红光波段的大气表观反射率ρR小于等于0.25时，且通道3可见光蓝光波段的大气表观反射率ρB大于0.28，则判定该区域为云、雾和高反射率地物；
若通道35热红外波段的亮度温度BT35不小于230且不大于250时，则判定该区域为云；
若通道2近红外波段与通道1可见光红光波段的大气表观反射率比值ρNR2/ρR为[0.9，1.6]时，或通道18近红外波段与通道2近红外波段的大气表观反射率比值ρNR18/ρNR2为[0.28，0.5]时，则判定该区域为云；
若通道26大气表观反射率ρNR26大于0.05时，则判定该区域为云；
若通道20热红外波段的亮度温度BT20与通道22热红外波段的亮度温度BT22的差值为[2.8，7.6]，或通道20热红外波段的亮度温度BT20与通道31热红外波段的亮度温度BT31的差值为[10，19]，则判定该区域为云和雾；
C2：确定气溶胶模式为大陆型并构造查找表；
将确定的大陆型气溶胶模式、气溶胶含量参数、卫星观测几何参数、波段的中心波长和探测器高度带入6S大气辐射传输模型中，计算得到不同波段及不同气溶胶光学厚度下的大气参数，从而构建并得到查找表；接收数据大气模式、光谱参数和地表反射率类型从6S辐射传输模型参数设置选项中选择；
其中，气溶胶模式为大陆型气溶胶模式，气溶胶含量参数设定为6个在波长0.55μm处的大气气溶胶光学厚度值：即0、0.25、0.5、1.0、1.5和1.95；卫星观测几何参数设定选取0°、6°、12°、24°、35.2°、48°、54°、60°和66°共9个太阳天顶角，从[0°，66°]范围内每隔6°的12个观测天顶角，从[0°，180°]范围内每隔12°的16个相对方位角；波段的中心波长取0.47μm和0.66μm；探测器高度从步骤B5中得到的L1B数据产品中获得；
C3：根据大气辐射传输原理，当采用标量法进行陆地气溶胶反演时，在地表朗伯体和大气水平均一假设条件下，大气表观反射率ρTOA表示为：
式中，
θs与θv分别为太阳天顶角与观测天顶角，φ为相对方位角，T为大气透过率；S为大气下界的半球反射率，ρ0是大气的路径辐射项等效反射率，ρs为地表反射率；下角标s表示太阳，v表示观测；
从步骤C2中得到的查找表中选取与卫星观测几何最邻近的大气参数，根据L1B级数据产品进行线性插值，得到不同波段和不同气溶胶光学厚度下的大气参数中大气下界的半球反射率S、大气的路径辐射项等效反射率ρ0和大气透过率T的数值，从MODIS多天历史的地表反射率数据中查找晴天状态下的地表反射率数据，构造晴天地表反射率库，得到所设定的时间段内的晴天状态下的地表反射率ρs，通过式(1)采用最小二乘法拟合模拟真实的大气表观反射率ρTOA，对真实的大气表观进行线性插值，最终获得大气气溶胶光学厚度。
7.根据权利要求6所述的基于天地一体化信息的PM2.5浓度监测方法，其特征在于：所述的步骤D中，采用的颗粒物反演模型基于地理加权回归的PM2.5反演模型，构建原理基于近地面PM2.5浓度和大气气溶胶光学厚度AOD之间存在的先验关系，即首先根据近地面PM2.5浓度和大气气溶胶光学厚度AOD之间的先验关系推导得到PM2.5的多元线性回归模型的对数形式；
近地面PM2.5浓度和大气气溶胶光学厚度AOD之间的先验关系为：
其中，ρ为PM2.5的平均密度，r为有效半径，Q为消光效率因子，AOD为大气气溶胶光学厚度，HPBL为边界层高度，f(RH)为气溶胶散射吸湿增长因子，一般表示为相对湿度的函数f(RH)＝1/(1.0-RH/100)；对上式两边同时取自然对数，得到PM2.5的多元线性回归模型的对数形式为：
其中，β1，β2，β3为回归系数，β1＝1，β2＝-1，β3＝-1；β0为回归常数；
将数据的地理位置，嵌入到PM2.5的多元线性回归模型的对数形式中的回归系数之中，得到卫星遥感监测PM2.5浓度的空间变化系数回归模型；
卫星遥感监测PM2.5浓度的空间变化系数回归模型为：
待估算点i处的PM2.5浓度yi为：
其中，(ui，vi)为待估算点i的经纬度坐标；β0(ui，vi)为回归常数；βd(ui，vi)表示为观测点(ui，vi)处的第d个回归系数，d＝1，2，3…；βd(ui，vi)是地理位置的函数，xid代表待估算点i的第d个与PM2.5相关的因子，xid包括近地面消光系数、边界层高度、相对湿度、近地面温度和风速；随机误差项服从正态分布εi～N(0，σ2)，随机误差相互独立，不存在序列相关Cov(εi，εj)＝0(i≠j，(i,j＝1,2,…,n))。
8.根据权利要求7所述的基于天地一体化信息的PM2.5浓度监测方法，其特征在于：还包括步骤F；
F：利用步骤E中获得的待监测区域PM2.5浓度分布的反演结果，制作待监测区域PM2.5浓度分布专题图。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910910189.0A CN110726653B (zh) | 2019-09-25 | 2019-09-25 | 基于天地一体化信息的pm2.5浓度监测方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910910189.0A CN110726653B (zh) | 2019-09-25 | 2019-09-25 | 基于天地一体化信息的pm2.5浓度监测方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN110726653A true CN110726653A (zh) | 2020-01-24 |
| CN110726653B CN110726653B (zh) | 2022-03-04 |
Family
ID=69219371
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910910189.0A Active CN110726653B (zh) | 2019-09-25 | 2019-09-25 | 基于天地一体化信息的pm2.5浓度监测方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110726653B (zh) |
Cited By (14)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111650128A (zh) * | 2020-06-08 | 2020-09-11 | 电子科技大学 | 一种基于地表反射率库的高分辨率大气气溶胶反演方法 |
| CN112378828A (zh) * | 2020-12-11 | 2021-02-19 | 中科三清科技有限公司 | 基于卫星遥感数据反演大气细颗粒物浓度的方法及装置 |
| CN112504144A (zh) * | 2020-12-04 | 2021-03-16 | 南京大学 | 一种海冰表面积雪厚度的遥感估算方法 |
| CN112798482A (zh) * | 2020-12-30 | 2021-05-14 | 联雨科技(天津)有限公司 | 一种基于卫星遥感的pm2.5与pm10估算方法 |
| CN113065098A (zh) * | 2021-03-23 | 2021-07-02 | 北京航天创智科技有限公司 | 基于卫星遥感和地气***数据的pm2.5浓度检测方法及装置 |
| CN114324074A (zh) * | 2021-12-22 | 2022-04-12 | 陕西智星空间科技有限公司 | 一种基于后向轨迹模型与遥感影像的大气监测走航方法 |
| CN114943142A (zh) * | 2022-04-29 | 2022-08-26 | 中国科学院空天信息创新研究院 | 高光谱地表反射率和大气参数一体化反演方法及装置 |
| CN115859026A (zh) * | 2022-11-18 | 2023-03-28 | 二十一世纪空间技术应用股份有限公司 | 一种高分辨率近地面pm2.5浓度遥感反演方法及装置 |
| CN116297288A (zh) * | 2023-05-24 | 2023-06-23 | 上海航天空间技术有限公司 | 大气甲烷干空气混合比的快速遥感反演方法及*** |
| CN116466368A (zh) * | 2023-06-16 | 2023-07-21 | 成都远望科技有限责任公司 | 基于激光雷达和卫星资料的沙尘消光系数廓线估算方法 |
| CN116698688A (zh) * | 2023-04-20 | 2023-09-05 | 兰州大学 | 一种基于风云4号双星协同估算大气颗粒物浓度的方法 |
| CN117390971A (zh) * | 2023-12-08 | 2024-01-12 | 北京英视睿达科技股份有限公司 | 一种基于供暖季的no2浓度的预估方法及*** |
| CN117555047A (zh) * | 2023-12-07 | 2024-02-13 | 中国人民解放军国防科技大学 | 基于人工智能的电力设备气象监测预警方法 |
| CN117576570A (zh) * | 2024-01-16 | 2024-02-20 | 山东大学 | 一种全天候水汽反演弹性方法、***、设备及介质 |
Citations (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101419282A (zh) * | 2008-12-05 | 2009-04-29 | 航天恒星科技有限公司 | 一体化高速遥感数据接收处理设备 |
| CN102103204A (zh) * | 2011-01-26 | 2011-06-22 | 环境保护部卫星环境应用中心 | 基于环境一号卫星的陆地气溶胶光学厚度反演方法 |
| KR101382507B1 (ko) * | 2012-10-19 | 2014-04-10 | 사단법인대기환경모델링센터 | 대기질 예측 및 관리 시스템 |
| CN104535979A (zh) * | 2014-12-23 | 2015-04-22 | 中国科学院遥感与数字地球研究所 | 一种陆地云光学厚度的遥感反演方法及*** |
| CN105092575A (zh) * | 2014-04-22 | 2015-11-25 | 株式会社日立制作所 | 评估沙尘暴强度的方法和装置 |
| CN105678085A (zh) * | 2016-01-12 | 2016-06-15 | 环境保护部卫星环境应用中心 | 一种pm2.5浓度的估算方法及*** |
| CN109001091A (zh) * | 2018-07-18 | 2018-12-14 | 北京航天宏图信息技术股份有限公司 | 大气污染的卫星遥感监测方法、装置和计算机可读介质 |
| CN109030301A (zh) * | 2018-06-05 | 2018-12-18 | 中南林业科技大学 | 基于遥感数据的气溶胶光学厚度反演方法 |
| CN110163035A (zh) * | 2018-02-11 | 2019-08-23 | 青岛星科瑞升信息科技有限公司 | 一种先验数据支持的云阴影识别方法 |
| CN110186823A (zh) * | 2019-06-26 | 2019-08-30 | 中国科学院遥感与数字地球研究所 | 一种气溶胶光学厚度反演方法 |
-
2019
- 2019-09-25 CN CN201910910189.0A patent/CN110726653B/zh active Active
Patent Citations (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101419282A (zh) * | 2008-12-05 | 2009-04-29 | 航天恒星科技有限公司 | 一体化高速遥感数据接收处理设备 |
| CN102103204A (zh) * | 2011-01-26 | 2011-06-22 | 环境保护部卫星环境应用中心 | 基于环境一号卫星的陆地气溶胶光学厚度反演方法 |
| KR101382507B1 (ko) * | 2012-10-19 | 2014-04-10 | 사단법인대기환경모델링센터 | 대기질 예측 및 관리 시스템 |
| CN105092575A (zh) * | 2014-04-22 | 2015-11-25 | 株式会社日立制作所 | 评估沙尘暴强度的方法和装置 |
| CN104535979A (zh) * | 2014-12-23 | 2015-04-22 | 中国科学院遥感与数字地球研究所 | 一种陆地云光学厚度的遥感反演方法及*** |
| CN105678085A (zh) * | 2016-01-12 | 2016-06-15 | 环境保护部卫星环境应用中心 | 一种pm2.5浓度的估算方法及*** |
| CN110163035A (zh) * | 2018-02-11 | 2019-08-23 | 青岛星科瑞升信息科技有限公司 | 一种先验数据支持的云阴影识别方法 |
| CN109030301A (zh) * | 2018-06-05 | 2018-12-18 | 中南林业科技大学 | 基于遥感数据的气溶胶光学厚度反演方法 |
| CN109001091A (zh) * | 2018-07-18 | 2018-12-14 | 北京航天宏图信息技术股份有限公司 | 大气污染的卫星遥感监测方法、装置和计算机可读介质 |
| CN110186823A (zh) * | 2019-06-26 | 2019-08-30 | 中国科学院遥感与数字地球研究所 | 一种气溶胶光学厚度反演方法 |
Non-Patent Citations (4)
| Title |
|---|
| 吴黎: "《MODIS遥感信息处理方法及应用》", 30 June 2017, 哈尔滨工程大学出版社 * |
| 夏安全: "基于MODIS影像的济南市气溶胶反演研究与***开发", 《中国优秀博硕士学位论文全文数据库 工程科技Ⅰ辑》 * |
| 邵彦川 等: "基于地理加权回归克里金的中国PM2.5浓度空间制图方法", 《遥感技术与应用》 * |
| 陈辉 等: "基于地理加权模型的我国冬季PM2.5遥感估算方法研究", 《环境科学学报》 * |
Cited By (24)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111650128B (zh) * | 2020-06-08 | 2021-07-13 | 电子科技大学 | 一种基于地表反射率库的高分辨率大气气溶胶反演方法 |
| CN111650128A (zh) * | 2020-06-08 | 2020-09-11 | 电子科技大学 | 一种基于地表反射率库的高分辨率大气气溶胶反演方法 |
| CN112504144A (zh) * | 2020-12-04 | 2021-03-16 | 南京大学 | 一种海冰表面积雪厚度的遥感估算方法 |
| CN112504144B (zh) * | 2020-12-04 | 2021-10-29 | 南京大学 | 一种海冰表面积雪厚度的遥感估算方法 |
| CN112378828A (zh) * | 2020-12-11 | 2021-02-19 | 中科三清科技有限公司 | 基于卫星遥感数据反演大气细颗粒物浓度的方法及装置 |
| CN112378828B (zh) * | 2020-12-11 | 2021-09-17 | 中科三清科技有限公司 | 基于卫星遥感数据反演大气细颗粒物浓度的方法及装置 |
| CN112798482A (zh) * | 2020-12-30 | 2021-05-14 | 联雨科技(天津)有限公司 | 一种基于卫星遥感的pm2.5与pm10估算方法 |
| CN112798482B (zh) * | 2020-12-30 | 2022-09-06 | 联雨科技(天津)有限公司 | 一种基于卫星遥感的pm2.5与pm10估算方法 |
| CN113065098A (zh) * | 2021-03-23 | 2021-07-02 | 北京航天创智科技有限公司 | 基于卫星遥感和地气***数据的pm2.5浓度检测方法及装置 |
| CN114324074A (zh) * | 2021-12-22 | 2022-04-12 | 陕西智星空间科技有限公司 | 一种基于后向轨迹模型与遥感影像的大气监测走航方法 |
| CN114943142B (zh) * | 2022-04-29 | 2023-11-28 | 中国科学院空天信息创新研究院 | 高光谱地表反射率和大气参数一体化反演方法及装置 |
| CN114943142A (zh) * | 2022-04-29 | 2022-08-26 | 中国科学院空天信息创新研究院 | 高光谱地表反射率和大气参数一体化反演方法及装置 |
| CN115859026B (zh) * | 2022-11-18 | 2023-12-05 | 二十一世纪空间技术应用股份有限公司 | 一种高分辨率近地面pm2.5浓度遥感反演方法及装置 |
| CN115859026A (zh) * | 2022-11-18 | 2023-03-28 | 二十一世纪空间技术应用股份有限公司 | 一种高分辨率近地面pm2.5浓度遥感反演方法及装置 |
| CN116698688A (zh) * | 2023-04-20 | 2023-09-05 | 兰州大学 | 一种基于风云4号双星协同估算大气颗粒物浓度的方法 |
| CN116297288B (zh) * | 2023-05-24 | 2023-08-04 | 上海航天空间技术有限公司 | 大气甲烷干空气混合比的快速遥感反演方法及*** |
| CN116297288A (zh) * | 2023-05-24 | 2023-06-23 | 上海航天空间技术有限公司 | 大气甲烷干空气混合比的快速遥感反演方法及*** |
| CN116466368A (zh) * | 2023-06-16 | 2023-07-21 | 成都远望科技有限责任公司 | 基于激光雷达和卫星资料的沙尘消光系数廓线估算方法 |
| CN116466368B (zh) * | 2023-06-16 | 2023-08-22 | 成都远望科技有限责任公司 | 基于激光雷达和卫星资料的沙尘消光系数廓线估算方法 |
| CN117555047A (zh) * | 2023-12-07 | 2024-02-13 | 中国人民解放军国防科技大学 | 基于人工智能的电力设备气象监测预警方法 |
| CN117390971A (zh) * | 2023-12-08 | 2024-01-12 | 北京英视睿达科技股份有限公司 | 一种基于供暖季的no2浓度的预估方法及*** |
| CN117390971B (zh) * | 2023-12-08 | 2024-04-12 | 北京英视睿达科技股份有限公司 | 一种基于供暖季的no2浓度的预估方法及*** |
| CN117576570A (zh) * | 2024-01-16 | 2024-02-20 | 山东大学 | 一种全天候水汽反演弹性方法、***、设备及介质 |
| CN117576570B (zh) * | 2024-01-16 | 2024-04-12 | 山东大学 | 一种全天候水汽反演弹性方法、***、设备及介质 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN110726653B (zh) | 2022-03-04 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110726653B (zh) | 基于天地一体化信息的pm2.5浓度监测方法 | |
| Martens et al. | Evaluating the land-surface energy partitioning in ERA5 | |
| Kumar et al. | Wintertime characteristics of aerosols at middle Indo-Gangetic Plain: Impacts of regional meteorology and long range transport | |
| Jackson et al. | Suomi‐NPP VIIRS aerosol algorithms and data products | |
| Nakajima et al. | Overview of the Atmospheric Brown Cloud East Asian Regional Experiment 2005 and a study of the aerosol direct radiative forcing in east Asia | |
| CN102539336B (zh) | 基于环境一号卫星的可吸入颗粒物估算方法及*** | |
| CN102103204B (zh) | 基于环境一号卫星的陆地气溶胶光学厚度反演方法 | |
| Gupta et al. | Applying the Dark Target aerosol algorithm with Advanced Himawari Imager observations during the KORUS-AQ field campaign | |
| Shikwambana et al. | Global distribution of aerosol optical depth in 2015 using CALIPSO level 3 data | |
| Thieuleux et al. | Remote sensing of aerosols over the oceans using MSG/SEVIRI imagery | |
| Dianov-Klokov et al. | A spectroscopic study of the global space-time distribution of atmospheric CO | |
| Popp et al. | Remote sensing of aerosol optical depth over central Europe from MSG‐SEVIRI data and accuracy assessment with ground‐based AERONET measurements | |
| Zavody et al. | The ATSR data processing scheme developed for the EODC | |
| Glantz et al. | Trends in MODIS and AERONET derived aerosol optical thickness over Northern Europe | |
| Barnes et al. | An assessment of diurnal and seasonal cloud cover changes over the Hawaiian Islands using Terra and Aqua MODIS | |
| CN116822141A (zh) | 利用卫星微光遥感反演夜间大气气溶胶光学厚度的方法 | |
| Verstraete et al. | Generating 275-m resolution land surface products from the Multi-Angle Imaging Spectroradiometer data | |
| Kim et al. | On air pollutant variations in the cases of long-range transport of dust particles observed in central Korea in the leeside of China in 2010 | |
| CN111914396B (zh) | 基于高分辨率dem数据的次网格地形三维地表太阳辐射强迫效应快速参数化方法 | |
| Ma Lu et al. | Photosynthetically active radiation separation model for high-latitude regions in agrivoltaic systems modeling | |
| Whitlock et al. | WCRP surface radiation budget shortwave data product description, version 1.1 | |
| Huh | Limitations and capabilities of the NOAA satellite advanced very high resolution radiometer (AVHRR) for remote sensing of the earth's surface | |
| Gao et al. | Vortex streets to the lee of Madeira in a kilometre-resolution regional climate model | |
| Gupta et al. | Retrieval of aerosols over Asia from the Advanced Himawari Imager: Expansion of temporal coverage of the global Dark Target aerosol product | |
| Peng et al. | An evaluation of the NOAA global daily gap-filled VIIRS surface albedo |
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 | ||
| TR01 | Transfer of patent right | ||
| TR01 | Transfer of patent right |
Effective date of registration: 20230828 Address after: 450000 Boxue Road, Zhengdong New District, Zhengzhou City, Henan Province, 36 Patentee after: NO.27 RESEARCH INSTITUTE, CHINA ELECTRONICS TECHNOLOGY Group Corp. Patentee after: Henan Fangda Space Information Technology Co.,Ltd. Address before: 450000 Boxue Road, Zhengdong New District, Zhengzhou City, Henan Province, 36 Patentee before: NO.27 RESEARCH INSTITUTE, CHINA ELECTRONICS TECHNOLOGY Group Corp. |