技术领域
[0001] 本
发明涉及一种PM2.5浓度监测方法，尤其涉及一种基于天地一体化信息的PM2.5浓度监测方法。
背景技术
[0002] 大气细颗粒物(PM2.5)指空
气动力学直径小于等于2.5μm的颗粒，它的化学成分主要有有机
碳、元素碳、
硝酸盐和
硫酸盐等，是大气污染的重要成分之一，也是造成灰霾的主要原因。大气细颗粒物极易附带有毒有害物质，长时间在空间中悬浮会被人体呼吸进入
呼吸道、支气管和
肺部，对人体健康造成很大影响。
[0003] 国家和各级环保部
门近年来采取多种措施改善空气
质量，对监测系统也尤为重视。传统地面监测手段存在观测点离散、分布不均匀和建造成本高等缺点，与传统的地面监测手段相比，极轨卫星遥感由于具有大范围空间尺度、连续对地观测成像和高时间
分辨率过境等优势，为近地面PM2.5浓度信息获取提供了多种可能手段。
[0004] 目前，利用遥感的手段，系统化全自动监测PM2.5颗粒物浓度的方法尚未普及，相关环保部门缺少一种基于天地一体化信息的PM2.5颗粒物浓度自动化流程监测方法来支持雾霾治理工作。
发明内容
[0005] 本发明的目的是提供一种基于天地一体化信息的PM2.5浓度监测方法，能够高效准确的实现
指定监测区域内的PM2.5浓度监测。
[0006] 本发明采用下述技术方案：
[0007] 一种基于天地一体化信息的PM2.5浓度监测方法，包括以下步骤：
[0008] A：从极轨卫星接收站实时接收并获取遥感卫星影像数据；
[0009] B：预处理模
块对获取的遥感卫星影像数据进行预处理，将遥感卫星影像数据转化为各级影像产品，然后通过几何校正和辐
亮度计算最终得到亮度
温度值和地表反射率值；
[0010] C：进行大气
气溶胶光学厚度的反演，即首先进行
云识别并剔除干扰像元的干扰，再确定气溶胶模式，利用6S大气
辐射传输模型构建查找表，然后通过步骤B中获得的各级影像产品结合查找表，反演得到已确定的气溶胶模式下的大气气溶胶光学厚度；
[0011] D：构建卫星遥感监测PM2.5浓度的空间变化系数回归模型；得到待估算点i处的PM2.5浓度yi为： i＝1,2,…,n；
[0012] 其中，(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))；
[0013] E：利用步骤D中得到的待估算点i处的PM2.5浓度yi的计算公式，计算得出待监测区域影像中所有的
像素点i处的PM2.5浓度yi，即待监测区域PM2.5浓度分布的反演结果。
[0014] 所述的步骤A中，从极轨卫星接收站实时接收并获取Terra和Aqua的MODIS
载荷多颗遥感卫星影像数据。
[0015] 所述的步骤A包括如下具体步骤：
[0016] A1：接收处理系统进行卫星轨道两行数据下载，完成两行数据格式化分解后进行两行数据解算，生成轨道预报文件；
[0017] A2：根据轨道预报文件中的卫星轨道数据，接收处理系统进行卫星轨道关联和冲突预设消解工作，最终生成
跟踪接收计划；
[0018] A3：接收处理系统根据生成的跟踪接收计划，进行时间符合，在时间符合成功后接收处理系统对运行调度控制指令进行分析，将运行调度控制指令发送至天线设备，完成卫星的捕获与跟踪；同时，接收处理系统根据预先设置的卫星数据接收设备参数，完成卫星数据
锁定与接收；
[0019] A4：卫星的数传频段为X波段，接收处理系统采用26.25Mbs-1和15.0Mbs-1的码速率分别接收Terra和Aqua卫星的X频段
信号，经过下
变频器处理输出得到中频信号，再由接收机解调后输出数字基带信号；
[0020] A5：根据不同卫星数传数据格式要求，接收处理系统通过对数字基带信号中
帧同步字段的搜索，完成数据帧的字节对齐、容错和填充；
[0021] A6：接收处理系统根据卫星数传数据加扰方式，按照多项式生成PN码与数据帧同步字段后面的图像数据进行异或运算进行解扰；再根据卫星数传数据编码方式，对经过解扰后的图像数据进行RS译码以及CRC校验，最终得到影像的原始码流数据，即遥感卫星影像数据。
[0022] 所述的步骤B包括如下具体步骤：
[0023] B1：预处理模块实时监控原始码流数据存储路径，当监测到有新的原始码流数据文件产生，且满足新的原始码流数据文件不被占用同时大小大于0KB的条件时，开始对新的原始码流数据文件进行预处理；
[0024] B2：预处理模块通过解析新的原始码流数据文件的文件名，根据不同卫星载荷名称自动加载不同卫星载荷相应的xml配置文件；
[0025] B3：预处理模块对新的原始码流数据文件按照不同卫星载荷的数据帧格式，进行高级在轨系统帧格式解析，即自动调用新的原始码流数据文件的标准帧并解码，解码包括虚拟通道数据单元分离以及从不同虚拟通道数据单元提取不同信源包，重新定格式后输出虚拟通道数据单元，生成L0级数据产品；
[0026] 其中，L0级数据产品，指经数据重构，且未进行任何处理的原始数据；L0级数据产品中所有的通信信息被移除，通信信息包括数据帧同步字段、通信头和重复数据；
[0027] B4：预处理模块结合辅助数据信息，从L0级数据产品中的CCSDS包解包还原出扫描数据的集合，即预处理模块通过提取辅助数据信息中的星历、
姿态、跳秒和世界标准时间信息，对图像进行校正，完成L1A级数据产品的生产；
[0028] 其中，辅助数据信息包括星历、姿态、跳秒、世界标准时间、辐射定标参数、反射率参数和质量参数信息；L1A级产品，指经数据重构，且具有时间参考及地理坐标参数信息的未进行任何处理的原始数据；
[0029] B5：在L1A级数据产品的
基础上，预处理模块利用辅助数据信息中的辐射定标参数、反射率参数和质量参数信息，结合L1A级产品中的地理坐标参数信息，根据不同卫星载荷进行
定位和定标运算，生成包括经纬度信息、探测器高度和不同通道的大气表观反射率与辐射率的数据集，即L1B级数据产品；
[0030] 其中，L1B级数据产品，指在L1A级数据产品的基础上，进行定位和定标处理后的数据集；
[0031] B6：预处理模块自动调取包含大气表观反射率、辐射率和经纬度信息的L1B级数据产品，通过几何校正和辐亮度计算最终得到亮度温度值和地表反射率值。
[0032] 所述的步骤C中，气溶胶模式为大陆型气溶胶模式。
[0033] 所述的步骤C中包括如下具体步骤：
[0034] C1：云识别模块从步骤B6中得到的不同通道的大气表观反射率值和亮度温度值，设定相应的
阈值分别对干扰像元进行识别并剔除干扰像元的干扰，干扰像元包括云、雾和高反射率地物；
[0035] 其中，使用的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热红外波段的亮度温度；
[0036] 若通道1可见光红光波段的大气表观反射率ρR大于0.25，则判定该区域为云、雾和高反射率地物；若通道1可见光红光波段的大气表观反射率ρR小于等于0.25时，且通道3可见光蓝光波段的大气表观反射率ρB大于0.28，则判定该区域为云、雾和高反射率地物；
[0037] 若通道35热红外波段的亮度温度BT35不小于230且不大于250时，则判定该区域为云；
[0038] 若通道2近红外波段与通道1可见光红光波段的大气表观反射率比值ρNR2/ρR为[0.9，1.6]时，或通道18近红外波段与通道2近红外波段的大气表观反射率比值ρNR18/ρNR2为[0.28，0.5]时，则判定该区域为云；
[0039] 若通道26大气表观反射率ρNR26大于0.05时，则判定该区域为云；
[0040] 若通道20热红外波段的亮度温度BT20与通道22热红外波段的亮度温度BT22的差值为[2.8，7.6]，或通道20热红外波段的亮度温度BT20与通道31热红外波段的亮度温度BT31的差值为[10，19]，则判定该区域为云和雾；
[0041] C2：确定气溶胶模式为大陆型并构造查找表；
[0042] 将确定的大陆型气溶胶模式、气溶胶含量参数、卫星观测几何参数、波段的中心
波长和探测器高度带入6S大气辐射传输模型中，计算得到不同波段及不同气溶胶光学厚度下的大气参数，从而构建并得到查找表；接收数据大气模式、
光谱参数和地表反射率类型从6S辐射传输模型参数设置选项中选择；
[0043] 其中，气溶胶模式为大陆型气溶胶模式，气溶胶含量参数设定为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数据产品中获得；
[0044] C3：根据大气辐射传输原理，当采用标量法进行陆地气溶胶反演时，在地表朗伯体和大气
水平均一假设条件下，大气表观反射率ρTOA表示为：
[0045]
[0046] 式中， θs与θv分别为太阳天顶角与观测天顶角，φ为相对方位角，T为大气透过率；S为大气下界的半球反射率，ρ0是大气的路径辐射项等效反射率，ρs为地表反射率；下角标s表示太阳，v表示观测；
[0047] 从步骤C2中得到的查找表中选取与卫星观测几何最邻近的大气参数，根据L1B级数据产品进行线性插值，得到不同波段和不同气溶胶光学厚度下的大气参数中大气下界的半球反射率S、大气的路径辐射项等效反射率ρ0和大气透过率T的数值，从MODIS多天历史的地表反射率数据中查找晴天状态下的地表反射率数据，构造晴天地表反射率库，得到所设定的时间段内的晴天状态下的地表反射率ρs，通过式(1)采用最小二乘法拟合模拟真实的大气表观反射率ρTOA，对真实的大气表观进行线性插值，最终获得大气气溶胶光学厚度。
[0048] 所述的步骤D中，采用的颗粒物反演模型基于地理加权回归的PM2.5反演模型，构建原理基于近地面PM2.5浓度和大气气溶胶光学厚度AOD之间存在的先验关系，即首先根据近地面PM2.5浓度和大气气溶胶光学厚度AOD之间的先验关系推导得到PM2.5的多元线性回归模型的对数形式；
[0049] 近地面PM2.5浓度和大气气溶胶光学厚度AOD之间的先验关系为：
[0050]
[0051] 其中，ρ为PM2.5的平均
密度，r为有效半径，Q为消光效率因子，AOD为大气气溶胶光学厚度，HPBL为边界层高度，f(RH)为气溶胶散射吸湿增长因子，一般表示为相对湿度的函数f(RH)＝1/(1.0-RH/100)；对上式两边同时取自然对数，得到PM2.5的多元线性回归模型的对数形式为：
[0052]
[0053] 其中，β1，β2，β3为回归系数，β1＝1，β2＝-1，β3＝-1；β0为回归常数；
[0054] 将数据的地理位置，嵌入到PM2.5的多元线性回归模型的对数形式中的回归系数之中，得到卫星遥感监测PM2.5浓度的空间变化系数回归模型；
[0055] 卫星遥感监测PM2.5浓度的空间变化系数回归模型为：
[0056]
[0057] 待估算点i处的PM2.5浓度yi为：
[0058]
[0059] 其中，(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))。
[0060] 还包括步骤F；
[0061] F：利用步骤E中获得的待监测区域PM2.5浓度分布的反演结果，制作待监测区域PM2.5浓度分布专题图。
[0062] 本发明通过获取遥感卫星影像数据，并将遥感卫星影像数据转化为各级影像产品，通过几何校正和辐亮度计算得到亮度温度值和地表反射率值，然后利用得到的亮度温度值和地表反射率值进行大气气溶胶光学厚度的反演，得到大气气溶胶光学厚度AOD；随后本发明通过构建卫星遥感监测PM2.5浓度的空间变化系数回归模型，计算得到待估算点i处的PM2.5浓度yi，并最终得到待监测区域PM2.5浓度分布的反演结果，实现基于天地一体化信息的PM2.5浓度监测。上述过程中，本发明基于影像光谱波段的不同反射特征，设定不同的阈值通过云识别模块把影像中的干扰像元检测出来，有效地降低了云、雾和高反射率地物对气溶胶光学厚度反演
精度的影响，同时，在卫星遥感监测PM2.5浓度的空间变化系数回归模型中，将数据的地理位置嵌入到PM2.5的多元线性回归模型的对数形式中的回归系数之中，使得回归系数随地理空间位置而变化，避免因回归系数单一而导致反演结果产生较大误差，也避免了现有线性回归模型中无法引入不同区域的气象参数和地表类型的
缺陷，进一步提高PM2.5浓度监测的准确度。
附图说明
[0063] 图1为本发明的流程示意图；
[0064] 图2为本发明中云识别模块的识别过程示意图。
具体实施方式
[0065] 以下结合附图和
实施例对本发明作以详细的描述：
[0066] 如图1所示，本发明所述的基于天地一体化信息的PM2.5浓度监测方法，包括以下步骤：
[0067] A：从极轨卫星接收站实时接收并获取遥感卫星影像数据；
[0068] 本发明中，从极轨卫星接收站实时接收并获取Terra和Aqua的MODIS载荷多颗遥感卫星影像数据。Terra和Aqua是美国地球观测系统计划中两颗重要卫星，Terra为上午星，Aqua为下午星，MODIS是搭载在Terra和Aqua卫星上的中分辨率成像光谱仪。
[0069] 本实施例中，步骤A包括如下具体步骤：
[0070] A1：接收处理系统进行卫星轨道两行数据下载，完成两行数据格式化分解后进行两行数据解算，生成轨道预报文件。
[0071] A2：根据轨道预报文件中的卫星轨道数据，接收处理系统进行卫星轨道关联和冲突预设消解工作，最终生成跟踪接收计划。
[0072] A3：接收处理系统根据生成的跟踪接收计划，进行时间符合，在时间符合成功后接收处理系统对运行调度控制指令进行分析，将运行调度控制指令发送至天线设备，完成卫星的捕获与跟踪；同时，接收处理系统根据预先设置的卫星数据接收设备参数，完成卫星数据锁定与接收。
[0073] A4：卫星的数传频段为X波段，接收处理系统采用26.25Mbs-1和15.0Mbs-1的码速率分别接收Terra和Aqua卫星的X频段信号，经过下变频器处理输出得到中频信号，再由接收机解调后输出数字基带信号。
[0074] A5：根据不同卫星数传数据格式要求，接收处理系统通过对数字基带信号中帧同步字段的搜索，完成数据帧的字节对齐、容错和填充。
[0075] A6：接收处理系统根据卫星数传数据加扰方式，按照多项式生成PN码(Pseudor-Noise伪随机码)与数据帧同步字段后面的图像数据进行异或运算进行解扰；再根据卫星数传数据编码方式，对经过解扰后的图像数据进行RS(Reed-Solomon)译码以及CRC校验，最终得到影像的原始码流数据，即遥感卫星影像数据。
[0076] 本实施例中，操作员可以对上述整个跟踪接收计划的生成和下发的工作状态，以及卫星跟踪和接收的设备控制响应情况进行实时监视。
[0077] B：预处理模块对获取的遥感卫星影像数据进行预处理，将遥感卫星影像数据转化为各级影像产品，然后通过几何校正和辐亮度计算最终得到亮度温度值和地表反射率值；
[0078] 本发明中，步骤B包括如下具体步骤：
[0079] B1：预处理模块实时监控原始码流数据存储路径，当监测到有新的原始码流数据文件产生，且满足新的原始码流数据文件不被占用同时大小大于0KB的条件时，开始对新的原始码流数据文件进行预处理；
[0080] B2：预处理模块通过解析新的原始码流数据文件的文件名，根据不同卫星载荷名称自动加载不同卫星载荷相应的xml配置文件；
[0081] B3：预处理模块对新的原始码流数据文件按照不同卫星载荷的数据帧格式，进行高级在轨系统帧格式解析，即自动调用新的原始码流数据文件的标准帧并解码，解码包括虚拟通道数据单元分离以及从不同虚拟通道数据单元提取不同信源包，重新定格式后输出虚拟通道数据单元，生成L0级数据产品；L0级数据产品，指经数据重构，且未进行任何处理的原始数据；L0级数据产品中所有的通信信息被移除，通信信息包括数据帧同步字段、通信头和重复数据。
[0082] B4：预处理模块结合辅助数据信息，从L0级数据产品中的CCSDS包解包还原出扫描数据的集合，即预处理模块通过提取辅助数据信息中的星历、姿态、跳秒和世界标准时间信息，对图像进行校正，完成L1A级数据产品的生产；其中，辅助数据信息包括星历、姿态、跳秒、世界标准时间、辐射定标参数、反射率参数和质量参数信息；L1A级产品，指经数据重构，且具有时间参考及地理坐标参数信息的未进行任何处理的原始数据。
[0083] B5：在L1A级数据产品的基础上，预处理模块利用辅助数据信息中的辐射定标参数、反射率参数和质量参数信息，结合L1A级产品中的地理坐标参数信息，根据不同卫星载荷进行定位和定标运算，生成包括经纬度信息、探测器高度和不同通道的大气表观反射率与辐射率的数据集，即L1B级数据产品。L1B级数据产品，指在L1A级数据产品的基础上，进行定位和定标处理后的数据集。
[0084] B6：预处理模块自动调取包含大气表观反射率、辐射率和经纬度信息的L1B级数据产品，通过几何校正和辐亮度计算最终得到亮度温度值和地表反射率值。
[0085] 要进行大气气溶胶光学厚度的反演，必须分离地表反射率的贡献，通过对云像元进行识别，剔除云、雾和高反射率地物的干扰，然后建立气溶胶模式，构建查找表，从而通过经步骤B预处理后得到的L1B级数据产品反演某种气溶胶模式下的气溶胶光学厚度，其中，高反射率地物指反射率不小于50％的非选择性反射体。因此，本
申请中，步骤C如下所述：
[0086] C：进行大气气溶胶光学厚度的反演，即首先进行云识别并剔除干扰像元的干扰，再确定气溶胶模式，利用6S大气辐射传输模型构建查找表，然后通过步骤B中获得的各级影像产品结合查找表，反演得到已确定气溶胶模式下的大气气溶胶光学厚度；
[0087] 大气气溶胶光学厚度，英文名称为AOD(Aerosol Optical Depth)或AOT(Aerosol Optical Thickness)，定义为介质的消光系数在垂直方向上的积分，是描述气溶胶对光的削减作用的。大气气溶胶光学厚度是气溶胶最重要的参数之一，表征大气浑浊程度的关键的物理量，也是确定气溶胶
气候效应的重要因素。
[0088] 本申请中，云识别模块是通过基于影像光谱波段的不同反射特征，设定不同的阈值把影像中的云检测出来，剔除云、雾和高反射率地物对气溶胶光学厚度反演的影响。
[0089] 气溶胶模式指一组定义好的气溶胶类型的集合，如标准气溶胶模型中的大陆型、海洋型和城市型，每种气溶胶模式是由四种不同的粒子组分按一定的比率组合而成，详见表1。所谓气溶胶类型不同，也就是气溶胶组分或其组合比率不同。实际反演中是用辐射传输模型根据假设的气溶胶模式构建查找表，由于不同大气状况下粒子组分的不同，对最终的反演结果会带来不同的影响，因此气溶胶模式是影响反演精度的一个重要因素。参考研究区的气溶胶吸湿增长和闭合实验等研究结果，我们设定待监测区域的气溶胶模式为大陆型气溶胶模式。
[0090] 表1城市型、大陆型和海洋型气溶胶模式不同粒子的组成对比
[0091]
[0092] 确定气溶胶模式后，首先需要进行查找表的构造。查找表是通过设定不同卫星观测几何参数、大气模式、气溶胶含量参数、探测器高度、光谱参数、地表反射率类型以及观测数据所在的波段中心波长，采用某种辐射传输方程来进行。本发明使用6S辐射传输模型进行辐射传输计算来构造查找表，6S辐射传输模型采用了最新的散射计算方法，使
太阳光谱波段的散射计算精度比5S有明显提高。
[0093] 本发明中，步骤C中包括如下具体步骤：
[0094] C1：云识别模块从步骤B6中得到的不同通道的大气表观反射率值和亮度温度值，设定相应的阈值分别对干扰像元进行识别并剔除干扰像元的干扰，干扰像元包括云、雾和高反射率地物。根据MODIS载荷的特点及中国华北地区多影像分析，通过可见光、近红外和红外的多个波段的大气表观反射率、亮度温度和空间变化特征对阈值进行调整和优化，最终本申请中所使用的云识别模块的识别过程如图2所示。
[0095] 其中，使用的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热红外波段的亮度温度。
[0096] 为了进一步提高云识别模块识别的准确性，尽可能地降低云、雾和高反射率地物对气溶胶光学厚度反演精度的影响，本申请根据不同通道下不同波段的可见光、近红外和红外的特性，对于云、雾和高反射率地物设定不同的大气表观反射率阈值进行判断。
[0097] 本申请利用MODIS多时相影像分析，若通道1可见光红光波段的大气表观反射率ρR大于0.25，则判定该区域为云、雾和高反射率地物；若通道1可见光红光波段的大气表观反射率ρR小于等于0.25时，且通道3可见光蓝光波段的大气表观反射率ρB大于0.28，则判定该区域为云、雾和高反射率地物；上述阈值能够有效判定干扰像元中的厚云、浓雾和高反射率地物的情况；
[0098] 若通道35热红外波段的亮度温度BT35不小于230且不大于250时，则判定该区域为云；上述阈值设定中，为了排除云识别过程时霾对于云的干扰，经分析得出，在陆地上空云区域具有较低的亮度温度值，而霾区域的亮度温度值相对高于云区域的亮度温度值，因此将通道35热红外波段的亮度温度BT35的阈值设定为[230，250]，以准确区分霾区域和云区域；
[0099] 若通道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]，以有效判定干扰像元中存在的薄云、碎云和层积云；
[0100] 若通道26大气表观反射率ρNR26大于0.05时，则判定该区域为云；经分析得出，通道26对于确定卷云具有明显优势，因此本申请中设定ρNR26的阈值为0.05，以有效判定干扰像元中存在的卷云；
[0101] 若通道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]，以有效判定干扰像元中存在的低层云和雾粒子。
[0102] 上述过程中，厚云、浓雾、薄云、碎云、层积云、卷云、低层云和雾粒子均为气象学领域的专业名词。
[0103] C2：确定气溶胶模式为大陆型并构造查找表；
[0104] 将确定的大陆型气溶胶模式、气溶胶含量参数、卫星观测几何参数、波段的中心波长和探测器高度带入6S大气辐射传输模型中，计算得到不同波段及不同气溶胶光学厚度下的大气参数，从而构建并得到查找表；接收数据大气模式、光谱参数和地表反射率类型从6S辐射传输模型参数设置选项中选择。
[0105] 其中，气溶胶模式为大陆型气溶胶模式，气溶胶含量参数设定为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数据产品中获得。
[0106] C3：根据大气辐射传输原理，当采用标量法进行陆地气溶胶反演时，在地表朗伯体和大气水平均一假设条件下，大气表观反射率ρTOA表示为：
[0107]
[0108] 式中， θs与θv分别为太阳天顶角与观测天顶角，φ为相对方位角，T为大气透过率；S为大气下界的半球反射率，ρ0是大气的路径辐射项等效反射率，ρs为地表反射率；下角标s表示太阳，v表示观测；
[0109] 从步骤C2中得到的查找表中选取与卫星观测几何最邻近的大气参数，根据L1B级数据产品进行线性插值，得到不同波段和不同气溶胶光学厚度下的大气参数中大气下界的半球反射率S、大气的路径辐射项等效反射率ρ0和大气透过率T的数值。从MODIS多天历史的地表反射率数据中查找晴天状态下的地表反射率数据，构造晴天地表反射率库，得到所设定的时间段内的晴天状态下的地表反射率ρs，通过式(1)采用最小二乘法拟合模拟真实的大气表观反射率ρTOA，对真实的大气表观进行线性插值，最终获得大气气溶胶光学厚度AOD。
[0110] D：构建卫星遥感监测PM2.5浓度的空间变化系数回归模型；
[0111] 本发明采用的颗粒物反演模型是基于地理加权回归的PM2.5反演模型；其构建原理基于近地面PM2.5浓度和大气气溶胶光学厚度AOD之间存在一定的先验关系，即首先根据近地面PM2.5浓度和大气气溶胶光学厚度AOD之间的先验关系推导得到PM2.5的多元线性回归模型的对数形式；
[0112] 近地面PM2.5浓度和大气气溶胶光学厚度AOD之间的先验关系为：
[0113]
[0114] 其中，ρ为PM2.5的平均密度，r为有效半径，Q为消光效率因子，AOD为大气气溶胶光学厚度，HPBL为边界层高度，f(RH)为气溶胶散射吸湿增长因子，一般表示为相对湿度的函数f(RH)＝1/(1.0-RH/100)；对上式两边同时取自然对数，得到PM2.5的多元线性回归模型的对数形式为：
[0115]
[0116] 其中，β1，β2，β3为回归系数，β1＝1，β2＝-1，β3＝-1；β0为回归常数；
[0117] PM2.5的多元线性回归模型的对数形式是一种全局空间回归模型，β1，β2，β3作为回归系数会随地理空间位置而变化。但现有的基于的区域AOD浓度遥感反演大多都没有考虑到局部参数的变化特征，根据一般经验β1，β2，β3分别取1，-1和-1。因此现有的回归系数单一的线性回归模型会导致反演结果有较大的误差，因此本发明中：
[0118] 将数据的地理位置，嵌入到PM2.5的多元线性回归模型的对数形式中的回归系数之中，得到卫星遥感监测PM2.5浓度的空间变化系数回归模型；
[0119] 卫星遥感监测PM2.5浓度的空间变化系数回归模型为：
[0120]
[0121] 待估算点i处的PM2.5浓度yi为：
[0122]
[0123] 其中，(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))；
[0124] E：利用步骤D中得到的待估算点i处的PM2.5浓度yi的计算公式，计算得出待监测区域影像中所有的像素点i处的PM2.5浓度yi，即待监测区域PM2.5浓度分布的反演结果；
[0125] 本发明中，还包括步骤F，根据获得的待监测区域PM2.5浓度分布的反演结果，制作待监测区域PM2.5浓度分布专题图；
[0126] F：利用步骤E中获得的待监测区域PM2.5浓度分布的反演结果，制作待监测区域PM2.5浓度分布专题图。
[0127] 本发明中，当系统接收到待监测区域PM2.5浓度分布的反演结果时，可自动进行专题图和专题报的生产。
[0128] 专题图制作所需文件包括：待监测区域PM2.5浓度分布的反演结果、专题图中显示范围的矢量文件、
渲染方式配置文件、比例尺、指北针及图例。
[0129] 专题报制作所需文件包括：待监测区域PM2.5浓度分布的反演结果、专题图及专题报模板。
[0130] 专题图制作的具体步骤如下：
[0131] 1、加载待监测区域PM2.5浓度分布的反演结果，加载待监测区域行政区边界数据，读取渲染配置文件；待监测区域行政区边界数据采用shp格式矢量文件；
[0132] 2、依据待监测区域行政区边界数据，对待监测区域PM2.5浓度分布的反演结果进行裁剪；
[0133] 3、根据渲染方式配置文件，对裁剪后的待监测区域PM2.5浓度分布的反演结果进行渲染；
[0134] 4、添加比例尺、指北针、专题图名称及文本标注；
[0135] 5、提取专题图的元数据，注册到
数据库，并输出专题图到指定路径下。
[0136] 专题报制作步骤如下：
[0137] 1、读取待监测区域PM2.5浓度分布的反演结果，依据PM2.5的像素值分为优(0＜像素值≤35)，良(35＜像素值≤75)，轻度污染(75＜像素值≤115)，中度污染(115＜像素值≤150)，重度污染(150＜像素值≤250)和严重污染(250＜像素值)6个污染等级；
[0138] 2、计算基础数据：待监测区域各行政区的有效监测面积，PM2.5浓度均值，以及待监测区域内优、良、轻度污染、中度污染、重度污染和严重污染6个污染等级所占面积；
[0139] 3、计算统计数据：依据步骤2得到的基础数据进行统计分析，统计出各行政区的最高污染等级面积以及分布情况，依据排名规则，对行政区进行排名。
[0140] 4、读取专题报模板和专题图，插入专题图，依据统计数据，制作统计表、生成文字段落。最后输出专题报到指定路径下，提取专题报的元数据，注册到数据库。
[0141] 本发明还可以通过读取不同时间段范围的PM2.5浓度分布的反演结果，进行像素值加权平均，得到PM2.5浓度分布的日、周、月及年均值数据，按照专题图、报制作步骤生成日、周、月和年等不同时间
频率的PM2.5浓度分布及专题报。
[0142] 此外，本发明最后还可对反演的PM2.5精度进行分析。
[0143] 本发明在制作PM2.5浓度分布监测专题图及专题报的同时，对注册到数据库中的反演结果进行精度分析。利用省/市/区域尺度范围内国控地基
站点PM2.5观测值为基准值，计算得到卫星过境时间先后半小时范围中省内各省/市/区域地基站点算术平均值Sn。同时，根据国控站点经纬度，计算得出周边10公里范围内的本发明所测得的有效PM2.5算术平均值，并作为地基站点所在省/市/区域卫星观测平均值Pm。计算省/市/区域PM2.5浓度产品精度公式为：
[0144]
[0145] 通过该公式计算分析，可对本发明所测得的省/市/区域尺度PM2.5浓度的精度结果进行直接验证，有利于相关部门直观的掌握相关省/市/区域的PM2.5污染情况。
