CN112905560A - 一种多源时空大数据深度融合的空气污染预测方法 - Google Patents

一种多源时空大数据深度融合的空气污染预测方法 Download PDF

Info

Publication number
CN112905560A
CN112905560A CN202110144010.2A CN202110144010A CN112905560A CN 112905560 A CN112905560 A CN 112905560A CN 202110144010 A CN202110144010 A CN 202110144010A CN 112905560 A CN112905560 A CN 112905560A
Authority
CN
China
Prior art keywords
variables
data
space
model
time
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
CN202110144010.2A
Other languages
English (en)
Other versions
CN112905560B (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.)
Institute of Geographic Sciences and Natural Resources of CAS
Original Assignee
Institute of Geographic Sciences and Natural Resources 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 Geographic Sciences and Natural Resources of CAS filed Critical Institute of Geographic Sciences and Natural Resources of CAS
Priority to CN202110144010.2A priority Critical patent/CN112905560B/zh
Publication of CN112905560A publication Critical patent/CN112905560A/zh
Application granted granted Critical
Publication of CN112905560B publication Critical patent/CN112905560B/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
    • G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
    • G06F16/20—Information retrieval; Database structures therefor; File system structures therefor of structured data, e.g. relational data
    • G06F16/21—Design, administration or maintenance of databases
    • G06F16/211—Schema design and management
    • G06F16/212—Schema design and management with details for data modelling support
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N33/00—Investigating or analysing materials by specific methods not covered by groups G01N1/00 - G01N31/00
    • G01N33/0004—Gaseous mixtures, e.g. polluted air
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
    • G06F16/20—Information retrieval; Database structures therefor; File system structures therefor of structured data, e.g. relational data
    • G06F16/24—Querying
    • G06F16/245—Query processing
    • G06F16/2458—Special types of queries, e.g. statistical queries, fuzzy queries or distributed queries
    • G06F16/2474—Sequence data queries, e.g. querying versioned data
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/04—Architecture, e.g. interconnection topology
    • G06N3/045—Combinations of networks
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Chemical & Material Sciences (AREA)
  • Health & Medical Sciences (AREA)
  • Data Mining & Analysis (AREA)
  • Databases & Information Systems (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • General Engineering & Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Computational Linguistics (AREA)
  • General Health & Medical Sciences (AREA)
  • Software Systems (AREA)
  • Molecular Biology (AREA)
  • Biophysics (AREA)
  • Computing Systems (AREA)
  • Biomedical Technology (AREA)
  • Fuzzy Systems (AREA)
  • Probability & Statistics with Applications (AREA)
  • Evolutionary Computation (AREA)
  • Artificial Intelligence (AREA)
  • Combustion & Propulsion (AREA)
  • Food Science & Technology (AREA)
  • Medicinal Chemistry (AREA)
  • Analytical Chemistry (AREA)
  • Biochemistry (AREA)
  • Immunology (AREA)
  • Pathology (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种多源时空大数据深度融合的空气污染预测方法，采集多源大数据并预处理；对气象数据反演获得高分辨率的地面气象参数；气溶胶参数、NO2遥感参数缺失反演及升尺度；提取交通变量、土地利用变量、社会经济及POI变量、时空变异变量；各类时空大数据的协变量数据时空融合，形成统一尺度及空间坐标的数据集；反演空气污染浓度高分辨率的地表参数；精度验证及评估；如若达标，输出结果；如若不达标，调整并循环训练，直至得到合理的模型及预测。本发明时空覆盖度大，通过高级的优化技术提高气象数据的栅格建模及卫星参数的插补，取得了较高的测试精度、较高的泛化性，且结果验证及循环的建模机制减少估计偏差，提高实际应用的效率。

Description

一种多源时空大数据深度融合的空气污染预测方法
技术领域
本发明涉及一种空气污染预测方法，尤其涉及一种多源时空大数据深度融 合的空气污染预测方法。
背景技术
研究表明空气污染对身体健康产生有害的影响，短期会导致哮喘、肺炎 等，长期对呼吸系统及循环系统产生有害影响，同肺癌、心血管疾病紧密相 关，也会对孕妇及儿童发育产生不利影响。在经济及科技不断发展的今天，如 何监测并有效地减少空气污染是当前重要的环境问题。当前，虽然许多地区都 设置有空气污染监测站点，用于监测空气污染水平；但空气污染监测站点整体 设置数量有限，有限的监测数据对于区域及人口时空分异大的实际情况是远远 不够的，如何根据有限的监测数据准确地反演空气污染物地面时空分布情况， 对于空气污染的监控是必要且重要的。
根据测量数据反演未抽样地点的空气污染浓度，最简单的就是线性回归模 型及空间插值如kriging即克里格方法，但是这不能很好地建模空气污染物浓度 同协变量之间的复杂关系，所以非线性时空建模是反演空气污染浓度时空变化 的主要方法。交通作为空气污染主要的排放来源之一，也逐渐被提取用于空气 污染浓度的时空变化的建模之中。近年以来，遥感技术也逐渐用于提取重要的 空气污染协变量，如从MODIS(Moderate-resolution Imaging Spectroradiometer) 卫星传感器提取气溶胶光学深度系数(aerosol optical depth，简称AOD)，可以从 AOD反演PM10(fine particulate matterwith aerodynamic diameter smaller than 10 μm)及PM2.5(fine particulate matterwith aerodynamic diameter smaller than 2.5μm) 浓度。而从Ozone MonitoringInstrument(OMI)传感器Aura提取的O3(ozone)及NO2(nitrogen dioxide)，可用于反演地面空气污染物O3及NO2的时空变化情况。 但卫星数据最大的问题是由于云覆盖、地面的高反射率及传感器本身的缺陷， 会导致有大量的缺失值产生，从而严重影响卫星产品的应用。
为解决卫星数据的缺失问题，现有技术设计了基于概化累加模型 (generalizedadditive model，简称GAM)的缺值插补的非线性方法，采用基于 GEOS 5-FP的气象数据、土地利用、CMAQ(Community Multiscale Air Quality)模拟数据等插补完整中国长三角地区的卫星气溶胶数据MAIAC (Multi-Angle Implementation ofAtmospheric Correction)AOD，然后用于估算 地面PM2.5的浓度；也有采用了类似的变量，通过普通的前馈神经网络进行缺值 插补。其他的一些方法还包括了采用均值、最近邻或其他粗分辨率如CMAQ的 模拟AOD代替缺失值。虽然现有技术已采用了非线性方法插补缺失值，比简单 的替代方法要好，但是他们的缺值插补方法是基于粗分辨率的气象数据(分辨 率25km x 25km)，同反演的地面空气污染浓度的分辨率(1x1km2)差距过大， 这时导致估计结果偏差原因之一。
由于影响空气污染的高分辨率的重要气象参数(空气温度、相对湿度、风 速及气压等)缺乏，将会导致在缺乏高分辨率气象参数情况下采用卫星与/或 其他模拟的气象变量反演地面空气污染物浓度的估计偏差。现有技术提出了融 合多源遥感信息与气候环境的综合性气象环境评估方法，其方法基于较粗的卫 星分辨率，且只是反演了季节性的空气污染浓度分布并进行评估，缺乏高时间 分辨率的结果。现有技术还提出了一种大气颗粒物星地综合定量遥感融合的反 演方法，该方法没有对卫星气溶胶的缺失进行有效的处理，会导致相当多的时 空点由于卫星数据的缺乏无法进行反演；而反演采用的气象等参数分辨率也较 粗。对于采用多源遥感数据融合，基于GEOS-FP(Goddard Earth ObservingSystem-Forward Processing)系列数反演PM2.5污染物，同样其用于反演的地面分 辨率较粗，难以反映地面的在细尺度下的变化情况。以及现有的深度森林算法 估计城市模型估算，这种方法会受到基于树学习模型的输入离散数据的限制， 在样本数较少时会导致表面建模不连续性。虽然目前已有采用深度学习的CNN 等进行空气污染的地表参数重构，但由于空气污染同影响因素的复杂非线性关 系，难以采用卷积网取得理想的结果，且过深的网络会导致梯度消失问题，影 响最后的估计结果。现有技术还提出基于空气质量数据与图像实现多源异构的 融合的PM2.5预测模型，该方法需要采集照片估算PM2.5的浓度,评估结果受到环 境散射光线的影响，估计精度有限。
发明内容
为了解决上述技术所存在的不足之处，本发明提供了一种多源时空大数据 深度融合的空气污染预测方法。
为了解决以上技术问题，本发明采用的技术方案是：一种多源时空大数据 深度融合的空气污染预测方法，包括以下步骤：
步骤一、采集多源大数据；
步骤二、对步骤一采集到的数据进行预处理；
步骤三、高时空分辨率气象数据的插补，反演地面统一坐标的地面气象参 数；
步骤四、气溶胶参数、NO2遥感参数缺失反演及升尺度；
步骤五、交通变量、土地利用变量、社会经济及POI变量、时空变异变量 的提取；
步骤六、将各类时空大数据的协变量数据进行时空融合，形成统一尺度及 空间坐标的数据集；
步骤七、反演空气污染浓度地表参数；
步骤八、精度验证及评估；
步骤九、验证达标与否；如若达标，进入步骤十一；如若不达标，进入步 骤十；
步骤十、对于不合理或不达标预测，调整预测协变量、超参数及限制性条 件循环训练，直至得到合理的模型及预测；
步骤十一、结果输出。
进一步地，步骤二中，数据预处理包括对数据进行质量控制、时空融合及 格式转换；质量控制是按照有效数据阈值，或者质量标记对数据进行清理，删 除无效数据；时空融合是对多源大数据在时空上进行融合，取得统一的目标时 空分辨率；格式转换是将格式不一致的数据，转换成统一的高分空间栅格数 据。
进一步地，步骤三中，采用一体化的插值方法进行高时空分辨率气象数据 的插补，汇总气象数据中的特征协变量，建立目标变量统一且引入注意力层的 全残差深度网络模型；对各特征协变量加权，设计多变量输出层，采用多个目 标变量同时输出，便于模型参数间共享；在建模过程中，令特征的注意力权重 输入为：C＝{ci}，并由此对输入的特征进行加权：
Figure BDA0002929512610000041
满足足
Figure BDA0002929512610000042
其中，x＝{xi}为特征值矢量输入，xi则为第i个特征，c＝{ci}为特征值矢 量的注意力权重，ci为第i个特征的注意力权重系数；
Figure BDA0002929512610000043
为被注意力权重 加权之后的特征输出，°代表矩阵元素相乘；F为特征的个数；
采用softmax层来实现注意力层，以满足式2的限制条件：
Figure BDA0002929512610000044
其中，ei为注意力层的第i单元的输入，ci为第i个单元的权重输出，exp(…) 为取指数函数；j为特征的索引，ej为注意力层的j索引的输入，F为特征的个 数；注意力的权重系数通过模型优化自适应求解；
多变量输出的损失函数为：
Figure BDA0002929512610000051
其中，N为样本数，
Figure BDA0002929512610000052
代表空气温度的损失，
Figure BDA0002929512610000053
代表相对湿度的损 失，
Figure BDA0002929512610000054
代表空气压强的损失，
Figure BDA0002929512610000055
代表风速的损失；yt、yr、yp及yw分 别代表空气温度、相对湿度、压强及风速的正则化后的观察值，而
Figure BDA0002929512610000056
Figure BDA0002929512610000057
及
Figure BDA0002929512610000058
分别代表空气温度、相对湿度、压强及风速的估 计值；θW,b为网络模型的需要优化的参数集合，Ω(θW,b)则为参数集θW,b的正 则化项，采用弹性网络进行正则化；
最后，对目标变量结果进行反正则化，获得同原数据尺度一致的估计值。
进一步地，引入全残差深度网络模型为基模型的聚集引导聚集建模，训练 多个基模型，计算各参数的均值及变差，获得更优的地面气象参数反演结果。
进一步地，步骤四中，气溶胶参数指的是MAIACAOD，NO2遥感参数指 的是OMI-NO2，在步骤三生成的高分气象参数基础上，分别建立全残差深度网 络的AOD及OMI-NO2的缺值反演模型；在此缺值反演模型中引入注意力层，提 升重要协变量权重，以输入协变量+缺值需要插补的变量共同作为输出的目标 变量；
再分别对MAIACAOD及OMI-NO2建立缺值反演及升尺度模型，采用MSE 作为目标损失函数：
Figure BDA0002929512610000059
其中，N为样本数，lMSE表示采用均方误差的损失函数，y代表了地面空气 污染监测站点的观察的NO2浓度，
Figure BDA00029295126100000510
为地面高分辨率的尺度转换后得到 的地面的OMI-NO2的代理变量，引入了坡度因子s及截距a，通过其转换成地面NO2，使其同地面之间的观察值误差最低，从而优化参数得到高分辨率的尺度 化结果，对于缺失值直接采用尺度转化模型输入协变量即可得到缺失值的估 计。
进一步地，步骤五中，交通变量的提取采用最近邻提取或缓冲区分析法； 缓冲区分析法计算缓冲区内主干道的长度，通过敏感性分析，选取最优缓冲距 离作为最后距离：
Figure BDA0002929512610000061
其中，d*代表最优缓冲距离，r为相关系数函数，rlen(di)代表以di为缓冲 距离时统计缓冲区内的道路长度，x则为目标点空气污染物的浓度，i为测试的 缓冲距离索引；
土地利用变量的提取通过计算一定缓冲区范围内土地利用所占的面积比 例，将此比例作为污染来源的代理变量；作敏感性分析，确定最优的缓冲距 离：
Figure BDA0002929512610000062
其中，d*代表最优缓冲距离，r为相关系数函数，lprop(di)代表以di为缓冲 距离时统计缓冲区内的土地利用面积比例，x则为目标点空气污染物的浓度，i 为测试的缓冲距离索引；
社会经济指人口密度及GDP，通过空间叠加分析获得样本相应的人口密度 及GDP；POI数据则先选择污染源相关的兴趣点，计算到最近的POI的距离；
时空变异变量的提取包括空间变量、时间变量的提取两方面，其中，空间 变量包括高程、坐标及其派生变量，用于捕捉空间变异信息；时间变量包括多 尺度时间信息，包括年-天、月份、星期、及年份。
进一步地，步骤七中，基模型采用引入注意力层的全残差深度网络模型， 输入变量包括步骤一至步骤五获得的参数数据，设输入变量个数为D个，注意 力节点也为D个，对关键变量进行加权，输出变量为M个，即需要估计的变量 个数；同时，对输出设定了限制性条件，通过限制性优化方法，取得符合限制 性条件的解；输出的损失函数采用下式：
Figure BDA0002929512610000071
Figure BDA0002929512610000072
Figure BDA0002929512610000073
其中，L(θW,b)代表了网络参数集θW,b的损失函数，N为样本数，
Figure BDA0002929512610000074
及
Figure BDA0002929512610000075
分别代表NO2及PM2.5正则化后的观察值，而
Figure BDA0002929512610000076
及
Figure BDA0002929512610000077
分别代 表了针对输入协变量集x的NO2及PM2.5的网络预测值，
Figure BDA0002929512610000078
及
Figure BDA0002929512610000079
分别代表NO2及PM2.5的MSE的损失函数，Ω(θW,b)则为参数集θW,b的正则化项，同理，采用 弹性网络实现正则化；式9及式10定义了NO2及PM2.5的网络预测值需要满足的 条件，其中，
Figure BDA00029295126100000710
及
Figure BDA00029295126100000711
分别表示NO2的取值的最小值与最大值，
Figure BDA00029295126100000712
及
Figure BDA00029295126100000713
分别表示PM2.5的取值的最小值与最大值；
将二污染物合在一起在一个网络中输出，采用bootstrap重抽样方法，对样 本及特征进行多次重复抽样，得到不同的样本及特征变量集，对这些不同的样 本集合对全残差深度模型结构进行样本数随机化，以获得差异较大的基模型， 将这些样本分别训练这些基模型，最后将训练的模型分别预测，预测结果反正 则化得到NO2与PM2.5的估计浓度，统计多个模型预测浓度得到估计平均值及其 变差。
进一步地，步骤八中，采用独立性验证方法对聚集引导得到的结果进行精 度验证及结果评估。
进一步地，步骤九中，验证总体的预测目标是否符合要求，验证各个协变 量对总的预测解释是否合理。
进一步地，步骤十一中，对于得到的合理的优化的训练模型及超参数，保 存相应的模型及参数，并将合理的预测结果输出，供下一步的应用使用。
本发明公开了一种多源时空大数据深度融合的空气污染预测方法，具有以 下有益效果：1)考虑的时空大数据更全面，影响因子从基本因素扩展到污染 源因素，覆盖的时空范围更广泛，时空分辨率较高，这些大样本及其协变量， 能充分代表研究对象的总体情况，充分考虑了影响空气污染物浓度分布的各要 素；2)不受粗分辨率气象变量对预测的影响，通过采用地面实测气象数据建 立高精度高分辨率的气象地表参数反演模型，更准确拟合高分辨率下地表气象 因子的时空分布，为后续的遥感参数及空气污染物浓度建模奠定基础；3)基 于气象参数反演结果建立高效的遥感气溶胶缺值插补及OMI-NO2的升尺度模型，获得全时空覆盖的高分辨率的遥感参数；4)采用最新的融合注意力的全 残差深度学习时空模型，充分输入变量的值范围，采用共享参数输出，高效的 限制性优化学习可使得本专利在高分辨率的气象地表参数估计、卫星参数缺值 插补及空气污染反演建模取得高的精度。
综合以上优点，本专利同现有方法相比时空覆盖度大，考虑的影响要素更 全面，使用了新的深度学习建模技术，通过高级的优化技术提高气象数据的栅 格建模及卫星参数的插补，从而使得本发明取得了较高的测试精度、较高的泛 化性，而本发明也通过结果验证及循环的建模机制减少估计偏差，提高实际应 用的效率。
附图说明
图1为本发明的流程示意图。
图2为本发明引入注意力机制的全残差深度网络的天模型结构图。
图3为本发明实施例预测的NO2栅格图。
图4为本发明实施例预测的PM2.5栅格图。
具体实施方式
下面结合附图和具体实施方式对本发明作进一步详细的说明。
现有的空气污染物浓度时空估算方法，基于协变量采用回归模型估计，但 采用的气象协变量空间分辨率有限，遥感数据存在大量缺失值，导致了时空不 能全覆盖；而捕捉污染物的来源及其时空分布的协变量有限；许多方法对训练 模型采用黑箱模型，缺乏有效性验证及不确定性度量，无对预测结果的纠偏机 制。在此背景下，本发明提出了一种多源时空大数据深度融合的空气污染预测 方法，采集了气象数据、卫星遥感气溶胶数据、同化数据、交通、土地利用等 多源大数据，设计了一体化的高精度高分辨率气象及空气污染的反演全残差深 度学习方法，通过学习-解译-提高循环式的反演方法，获得更为可靠的空气污 染的地表参数的反演结果。
如图1所示为本发明的流程示意图，主要包括如下步骤：
步骤一、采集多源大数据；
采集的多源大数据包括几个方面，如每小时空气污染物(PM2.5,PM10,SO2 (sulfurdioxide),NO2,O3,CO(carbon monoxide),AQI(air quality index)数 据；气象测量数据包括：温度、相对湿度、气压、风速；气象再分析数据：行 星边界层高度(planetary boundarylayer height,PBLH)、云比例、大气臭氧总 量，地表气温、气压、相对湿度和风速，以及26种各种污染来源的MERRA2 GMI(The Modern-Era Retrospective analysis for Researchand Applications, Version 2，the Global Modeling Initiative)同化数据(表1，来源：https://acd- ext.gsfc.nasa.gov/Projects/GEOSCCM/MERRA2GMI)；交通数据：交通路网；土地 利用：开源的数据；社会经济：POI、GDP(Gross domestic product)、人口密 度；其他数据：中国地区1km数字高程模型的高程、x及y坐标、AERONET (Aerosol Robotic Network)AOD地面验证数据。
表1来自MERRA2 GMI同化数据拟合的污染物表面栅格相关变量列表
Figure RE-GDA0002989988400000101
步骤二、对步骤一采集到的数据进行预处理；
数据预处理主要包括对数据进行质量控制、时空融合及格式转换等；质量 控制是按照有效数据阈值，或者质量标记对数据进行清理，删除无效数据；采 用统计标注，删除一些对建模无意义的孤立点，采用阈值定义：
x∈[Q1(x)-5*IQR(x),Q3(x)+5*IQR(x)] 式1
其中，x∈x＝{xi(i＝1,...,N)}，N为数据总数，Q1(x)与Q3(x)分别表示对数 据集取第一及第三分位数，而IQR(x)表示对数据集x取四分位间距(interquartile range)。
时空融合是对多源大数据在时空上进行融合，以便取得统一的目标时空分 辨率；具体而言，对空间覆盖率范围及分辨率不一致的，将采取空间分析的系 列操作，包括合并、裁切、重投影、重采样及平均等统一空间坐标及空间分辨 率；对时间分辨率不一样的，采样双线性插值或平均化等统一空间分辨率。
格式转换是将格式不一致的(如文本格式)，转换成统一的高分空间栅格 数据。
步骤三、高时空分辨率的气象数据的插补，反演地面统一坐标的地面气象 参数；
采用地面气象监测点，结合地面坐标、高程、气象同化数据(PBLH、各 气象要素)等反演地面统一坐标的地面气象参数。采样一体化的插值方法进行 高分气象数据插补，考虑到气象要素间相互影响，设计了一体化全残差气象要 素集成化反演方法。
模型的输入包括了气象同化数据的气温、风速(包括地面、2米及10米的 高东西及南北向风速)、相对湿度、气压、臭氧、PBLH，以及地面坐标系列 变量，包括x、y、x2、y2及xy，高程，总计输入16个变量，输出为地面需要反演 的气温、气压、相对湿度及风速4个目标变量，同地面的监测站点的四个变量 相对应。基于全残差深度网络模型，引入了注意力层，对各特征协变量加权， 提高重要协变量的预测力，设计多变量输出层(4个目标变量同时输出)。由 于气温、气压、相对湿度及风速间相互影响，故采用多目标变量输出便于模型 参数间共享，可减少过拟合。模型结构框架参见图2(D＝16，M＝4)。注意力 权重的引入主要是为了加权输入特征，提高一些关键特征在预测中的重要性， 是一种自适应优化的建模方法。
在建模过程中，令特征的注意力权重输入为：C＝{ci}，并由此对输入的特 征进行加权：
Figure BDA0002929512610000111
满足足
Figure BDA0002929512610000112
其中，x＝{xi}为特征值矢量输入，xi则为第i个特征，c＝{ci}为特征值矢 量的注意力权重，ci为第i个特征的注意力权重系数；
Figure BDA0002929512610000113
为被注意力权重 加权之后的特征输出，°代表矩阵元素相乘；F为特征的个数；
采用softmax层来实现注意力层，以满足式2的限制条件：
Figure BDA0002929512610000121
其中，ei为注意力层的第i单元的输入，ci为第i个单元的权重输出，exp(…) 为取指数函数；j为特征的索引，ej为注意力层的j索引的输入，ej与ei一样，均 为注意力层的输入，不过j索引表示对所有的ej；F为特征的个数；注意力的权 重系数通过模型优化自适应求解；
多变量输出的损失函数为：
Figure BDA0002929512610000122
其中，N为样本数，
Figure BDA0002929512610000123
代表空气温度的损失，
Figure BDA0002929512610000124
代表相对湿度的损 失，
Figure BDA0002929512610000125
代表空气压强的损失，
Figure BDA0002929512610000126
代表风速的损失；yt、yr、yp及yw分 别代表空气温度、相对湿度、压强及风速的正则化后的观察值，而
Figure BDA0002929512610000127
Figure BDA0002929512610000128
及
Figure BDA0002929512610000129
分别代表空气温度、相对湿度、压强及风速的估 计值；θW,b为网络模型的需要优化的参数集合，Ω(θW,b)则为参数集θW,b的正 则化项，采用弹性网络进行正则化；
最后，对目标变量结果进行反正则化，获得同原数据尺度一致的估计值。
同时，本模型通过随机梯度下降法获得最优参数，并通过一个模型估计得 到了4个气象参数的估计值，提高了网络参数的共享，满足测试精度的同时提 高了计算效率。
为进一步提高参数估计的精准性，引入全残差深度网络模型为基模型的聚 集引导聚集建模，训练多个(100个)基模型，计算各参数的均值及变差，提 高模型的泛化性及气象参数预测的稳定性，获得更优的地面气象参数反演结 果。
步骤四、气溶胶参数、NO2遥感参数缺失反演及升尺度；
气溶胶参数指的是MAIAC AOD，采用了高分的MAIAC AOD来作为PM2.5预测主因子之一；NO2遥感参数指的是OMI-NO2，采用了OMI-NO2来作为地面 的NO2的主因子变量之一。
对于MAIAC AOD及OMI-NO2存在大量的缺失值，缺失值达到50％以上。在 步骤三生成的高分气象参数基础上，加上了气象同化数据的较粗分辨率的协变 量(包括背景AOD或NO2、PBLH、云比例)、坐标(x与y)及其派生变量 (x2，y2及xy)及高程总共9个变量，分别建立全残差深度网络的AOD及OMI- NO2的缺值反演模型，采用了类似图2的网络结构；在此缺值反演模型基础上， 同样按照注意力机制引入注意力层，提升重要协变量权重，而在输出变量方 面，以输入协变量+缺值需要插补的变量(即10个输出变量)共同作为输出的 目标变量；之所以加入输入变量也作为输出变量，主要是加强变量间的数共 享，对于采用点缺值建模样本而言样本数足够大，将输入也作为输出一部分同 目标变量一起输出，起到了参数正则化效果，可取得较好的拟合精度。
同时，本发明分别对MAIAC AOD及OMI-NO2建立缺值反演及升尺度模 型；在升尺度方面，直接考虑了升尺度后的OMI-NO2同地面NO2之间的相关 性，采集连续3天的样本合成训练及测试样本，训练全残差的深度网络，以提 高同目标变量代理变量即地面NO2间的相关性为损失函数，但相关系数作为损 失函数训练可能不能优化，采用了他们之间MSE作为目标损失函数：
Figure BDA0002929512610000131
其中，N为样本数，lMSE表示采用均方误差(mean square error，即MSE) 的损失函数，y代表了地面空气污染监测站点的观察的NO2浓度，
Figure BDA0002929512610000132
为地 面高分辨率的尺度转换后得到的地面的OMI-NO2的代理变量，引入了坡度因子 s及截距a，通过其转换成地面NO2，使其同地面之间的观察值误差最低，从而 优化参数得到高分辨率的尺度化结果，对于缺失值直接采用尺度转化模型输入 协变量即可得到缺失值的估计。
步骤五、交通变量、土地利用变量、社会经济及POI变量、时空变异变量 的提取；
交通变量提取：地面空气污染物如NO2同交通紧密相关，需要提取交通相 关的变量。由于具体交通流量信息很难获取，本发明基于主要交通干道，包括 各级高速公路、主要道路提取2个关键的交通要素：到最近交通主干道的距 离、一定缓冲距离范围内道路的长度。对前者，采用最近邻提取；而对于后 者，采用缓冲区分析方法，计算缓冲区内交通主干道的长度。通过敏感性分 析，即计算一系列的距离缓冲区内道路长度同现有污染物NO2浓度间的相关 性，选取最优缓冲距离作为最后距离：
Figure BDA0002929512610000141
其中，d*代表最优缓冲距离，r为相关系数函数，rlen(di)代表以di为缓冲 距离时统计缓冲区内的道路长度，x则为目标点空气污染物的浓度，i为测试的 缓冲距离索引。
土地利用变量提取：土地利用也提供了主要的排放源信息。本发明将污染 来源相关的数据划分为一类，通过计算一定缓冲区范围内土地利用所占的面积 比例，将此比例作为污染来源的代理变量；同交通变量提取类似，作敏感性分 析，确定最优的缓冲距离：
Figure BDA0002929512610000142
其中，d*代表最优缓冲距离，r为相关系数函数，lprop(di)代表以di为缓冲 距离时统计缓冲区内的土地利用面积比例，x则为目标点空气污染物的浓度，i 为测试的缓冲距离索引。
社会经济及POI提取：此处的社会经济指人口密度及GDP，通过空间叠加 分析获得样本相应的人口密度及GDP，通常更高的GDP意味着更高的排放源； POI数据则先选择污染源相关的兴趣点，计算到最近的POI的距离。
时空变异变量提取：包括空间变量、时间变量的提取两大方面，其中，空 间变量包括高程、坐标及其派生变量，即x、y、x2、y2及xy，用于捕捉空间变异 信息；时间变量包括多尺度时间信息，包括年-天、月份、星期、及年份， “年-天”捕捉天尺度的时间变化，“月份”捕捉粗尺度季节性变化，“星 期”用于捕捉是否周末信息，“年份”用于捕捉年纪变化。
步骤六、将各类时空大数据的协变量数据进行时空融合，形成统一尺度及 空间坐标的数据集；
如数据在粗细尺度不一致需采用线性插值或非线性的升尺度方法统一尺 度，对于空间变量(如坐标)，需要衍生到各个时间点；对于多尺度时间变量 (如年天)，也需要衍生到各个空间点；最后形成统一的数据集。
步骤七、反演空气污染浓度地表参数；
输入变量包括步骤一至步骤五获得的参数数据，包括气象反演参数4个 (温度、气压、相对湿度及风速)，卫星参数2个(AOD及OMI-NO2)，气象 同化参数30个(PBLH、背景气溶胶、背景O3、背景NO2、及其他如表1所列的 26个污染源数据)，交通变量2个，土地利用变量1个，社会经济变量2个，POI 变量1个，空间变量6个，多尺度时间变量4个，总计达到52个变量。变异模型 的结构如图2所示(D＝52,M＝2)，基模型采用包括注意力层的全残差深度网络模型，输入变量个数D＝52个，注意力节点也为52个，对关键变量进行加权， 输出变量为M个，即需要估计的变量个数；例如：估计PM2.5及NO2，则M＝2。 此处同时对输出设定了限制性条件，通过限制性优化方法，尽量取得符合限制 性条件的解。输出的损失函数可以采用下式：
Figure BDA0002929512610000161
Figure BDA0002929512610000162
Figure BDA0002929512610000163
其中，L(θW,b)代表了网络参数集θW,b的损失函数，N为样本数，
Figure BDA0002929512610000164
及
Figure BDA0002929512610000165
分别代表NO2及PM2.5正则化后的观察值，而
Figure BDA0002929512610000166
及
Figure BDA0002929512610000167
分别代 表了针对输入协变量集x的NO2及PM2.5的网络预测值，
Figure BDA0002929512610000168
及
Figure BDA0002929512610000169
分别代表NO2及PM2.5的MSE的损失函数，Ω(θW,b)则为参数集θW,b的正则化项，同理，采用 弹性网络实现正则化；式9及式10定义了NO2及PM2.5的网络预测值需要满足的 条件，即需要分别落在在区间
Figure BDA00029295126100001610
及
Figure BDA00029295126100001611
之内，其中
Figure BDA00029295126100001612
及
Figure BDA00029295126100001613
分别表示NO2的取值的最小值与最大值，
Figure BDA00029295126100001614
及
Figure BDA00029295126100001615
分别表示PM2.5的 取值的最小值与最大值；
将二污染物合在一起在一个网络中输出，主要考虑PM2.5与NO2关系密切， 派生出NO2的交通污染也是PM2.5的来源之一，将二者同时作为一个模型输出， 有助于网络参数共享，减少过拟合，提高泛化性。
对于加入的限制性条件，采用投影梯度下降进行求解，虽不能保证全局最 优解(深度模型不是一个凸模型)，但局部最优解也基本满足要求。
为进一步提高反演稳定性，采用bootstrap重抽样方法，对样本及52个特征 进行多次(100次)重复抽样，得到100个不同的样本及特征变量集，对这些不 同的样本集合对全残差深度模型结构进行样本数随机化，以获得差异较大的基 模型，将这些样本分别训练这些基模型，最后将训练的模型分别预测，将多个 模型预测结果计算平均及变差，将结果反正则化得到NO2与PM2.5的估计浓度。
步骤八、精度验证及评估；
采用独立性验证方法对聚集引导得到的结果进行精度验证及结果评估。主 是对训练完成的模型对独立性样本计算R2及RMSE，获得精度评价标准。采用 SHAP(ShapleyAdditive exPlanations)方法对模型解译，计算各变量对预测贡 献。
步骤九、验证达标与否；
验证总体的预测目标是否符合要求，验证各个协变量对总的预测解释是否 合理。如模型解译表明交通路线促进污染物的下降说明是不合常理的，需要调 整模型进行重新建模训练，进入步骤十，否则保存模型，进入到步骤十一。
步骤十、对于不合理或不达标预测，调整预测协变量、超参数及限制性条 件循环训练，直至得到合理的模型及预测；
检查预测协变量提取是否有问题，是否真正反应两变量间关系，对模型训 练的超参数进行调整，包括学习率、训练的小批次样本大小(mini batch size)、网络的深度及节点数，采用启发式优化搜索方法，获取最优超参数进 入第七步进行下一步，循环直至得到合理的模型及预测。
步骤十一、结果输出。
对于得到的合理的优化的训练模型及超参数，保存相应的模型及参数，便 于以后的空气污染物地表浓度的反演及应用；并将合理的预测结果输出，供下 一步的应用使用。
由此，对于本发明所公开的一种多源时空大数据深度融合的空气污染预测 方法，具有以下技术优势：
1)涵盖了大的研究区域及多年的时段，采集了多种时空大数据。同其他 同类方法比较，本专利空间及时间覆盖范围宽，样本采集充足，能充分体现建 模总体的时空变化；采集的要素全面，包括气象测量、气象及空气污染同化数 据、卫星遥感气溶胶、高程、交通、土地利用、POI、社会经济及坐标等多达 52个变量，尤其是采集了NASA的多达26种污染源栅格背景同化数据。宽泛的 时空覆盖度及充足的预测因子，可以充分捕捉空气污染物的来源及时空分布， 大幅减少预测时的偏差。
2)高精度高分辨率气象地面参数反演，作为影响空气污染时空分布的关 键参数，现有方法通常使用了较粗分辨率的背景气象参数，没有经过地面测量 数据的矫正，空间分辨率过于粗糙，难以反映高空间分辨率如1km下的气象参 数的尺度变化，本发明采集了地面实测气象资料，用气象再分析背景数据、坐 标、高程等数据进行反演，获得了高精度高时空分辨率的气象数据，极大地提 高了空气污染时空建模的效果。
3)卫星气溶胶等数据的缺值插补，针对卫星气溶胶等参数大量的缺失 值，许多现有方法无法直接处理卫星缺失值，采用最近的值或均值替代，有的 方法采用了插补的非线性回归，但精度有限。本发明依托高精度高分辨率的气 象数据，加上背景同化数据、坐标及高程等，实现了高精度的全残差深度网络 缺值反演模型，取得了高精度的反演结果。对于粗尺度缺失的关键变量，采用 地面代理变量优化，取得了升尺度的较好的结果。同现有方法相比，本发明插 补完成了可靠的卫星气溶胶等数据，在空气污染地表参数反演时无需担心缺失 值，可实现时空全覆盖的估计。
4)在空气污染物浓度地表参数反演的建模方法进行了改进，同现有非线 性方法相比，在全残差深度模型基础上，加入了注意力层提高重要影响因子的 权重，采用了多变量输出方式使得变量之间实现参数共享，减少模型训练过程 的过拟合，加入限制性条件，通过限制优化方法求得更符合预期的解。因此本 发明采用的的全残差方法有效解决了深度网络导致的梯度消失问题，加入注意 力层加强重要因子的权重，多变量输出及限制性优化提高计算效率同时取得更 为合理的解。此方法也高效地用到气象参数重构及遥感参数缺值反演之中去。
5)注重对预测结果的解译及有效性验证，采用了可解释的机器学习技术 提取各影响因子对预测的贡献，对于不合理的预测进行干预，通过模型的回调 及再训练，消除预测结果中的偏差，更全面地获得预测的预测结果。同现有 的方法相比，本发明可以分析结果偏差的原由，通过循环纠偏过程，提高模型 应用的有效性。
【实施例】
下面结合具体的实施例，对本发明所公开的多源时空大数据深度融合的空 气污染预测方法做进一步详细的介绍。
本实施例以覆盖中国大陆地区的空间范围，时间覆盖2015-2018年4年，目 标空间分辨率为1x1km2，为时间分辨率为天，以此时空大数据地表参数反演中 国大陆局部地区的时空污染地表浓度。
步骤一、数据采集：采集覆盖中国大陆地区的高时空分辨率时空大数据集 (区域：中国大陆；时间：2015-2018年；空间分辨率1x1km；时间分辨率： 天)。具体包括：
地面测量数据：环境监测站点的空气污染数据从数据共享网站 https://quotsoft.net/air/下载得到；从中国气象数据网下载覆盖中国大陆地区的气 象参数(包括温度、相对湿度、气压、风速)测量值；
气溶胶及OMI-NO2数据：从MODIS Land Team网站获取了MAIAC AOD的 气溶胶光学深度数据(空间分辨率：1x1km，时间分辨率：天)；从 https://aeronet.gsfc.nasa.gov网站获得AERONET AOD数据集，用于验证MAIAC AOD缺值反演精度；从Earth Observing System,Aura网站获得了OMI-NO2数据 (空间分辨率：0.25°lonx0.25°lat，时间分辨率：天)；
气象参数同化背景数据：从全球土地数据同化系统(Global Land DataAssimilation System，GLDAS)获得了粗分辨率的气象同化的数据(温度、相 对湿度、风速；空间分辨率：0.25°lon x0.25°lat；时间分辨率：3小时)；从戈 达德地球观测系统-前向处理(Goddard Earth Observing System-Forward Processing，GEOS-FP)获得臭氧及PBLH数据(空间分辨率：0.25°lon x0.25°lat；时间分辨率：3小时)；从MERRA2 GMI同化数据源获得了26种污 染来源的栅格数据(表1；空间分辨率：0.625°(经度)x 0.5°(维度)；时间 分辨率：3小时)；
高程数据来自资源环境数据云平台的500m空间分辨率的DEM数据，其数 据源自对航天飞机雷达地形测绘任务(Shuttle Radar Topography Mission， SRTM)重采样获得；
交通、土地利用及POI，来自于OpenStreet网站；
GDP(Gross domestic product)及人口密度：GDP来自于资源环境科学与 数据中心；人口密度来自于国家地球系统科学数据中心。
坐标数据x,y,x2,y2,xy根据生成的目标坐标网格中心点直接提取，目标 网格采用1954北京坐标系，空间分辨率为1x1km2。
步骤二、数据预处理：先对采集的数据进行预处理，包括删除无效的测量 数据，如根据MAIAC AOD提供的质量控制标记确定无效的像素值，根据阈值 定义删除极值点；之后，进行时空融合，通过合并、裁切、重投影、重采样及 平均等多种操作统一各类来源数据的空间坐标及空间分辨率。
步骤三、高精度高分辨率的气象栅格数据的反演。根据采集得到的主要气 象测量数据及气象同化数据，总共汇总得到16个协变量，建立气温、气压、相 对湿度及风速多变量输出的统一的引入了注意力层的深度残差回归网络，训练 100个基模型，最后得到预测输出。
表2报告各个气象要素的测试精度(R2及RMSE)，高精度的测试结果 (R2≥0.81)表明一体化的全残差深度网络取得了较好的反演效果。全残差深度 模型均采用了基于Tensorflow的Keras建立模型，基于Python语言实现。而训练 模型采用的节点数依次是[16,96,64,32,16,8,16,32,64,96,16,4]，其中包括16 个输入，4个目标变量的输出(包括温度、相对湿度、风速、气压)，选取的 训练超参数：批学习样本大小为1024，初始的学习率为0.1，采用了Adam的梯 度下降法优化方法。
表2高分辨率的气象参数反演的测试精度
Figure BDA0002929512610000211
步骤四、气溶胶参数、NO2遥感参数缺失反演及升尺度；气溶胶MAIAC AOD虽然提供了同目标坐标一致的空间分辨率，即1x1km2，但包含大量的缺失 值。采用了气象同化数据(背景AOD、PBLH、云比例)、坐标及其派生变 量、高程共9个输入变量，建立引入注意力机制的全残差深度网络的天模型 (图2)，反演MAIAC AOD同这些输入参数间的关系，最后得到缺值插补的天 模型，采用训练得到的天模型，完成缺值插补的功能，总计每天一个模型，总结从2015-2018年建立了1461个模型，平均的R2为0.90，最小的为测试R2为 0.76。对于OMI-NO2，有大量的缺失值且分辨率较粗，如前所述，采用地面实 测的空气污染物NO2作为矫正变量，以提高二者间相关性为目标，建立引入注 意力机制的全残差深度网络升尺度模型，采用可靠的空间点得到训练及测试样 本，训练模型，将训练后的模型进行升尺度，同时进行缺值插补，得到与目标 尺度一致的插补完成的OMI-NO2地面代理变量值，同理建立了1461个升尺度模 型，获得平均的R2为0.88，最小的R2为0.74，表明本升尺度方法取得了理想的效果，同时提高了OMI-NO2地面代理变量同预测目标即地面NO2之间的相关关 系。训练的网络模型采用的节点数依次是[9,128,64,32,16,8,16,32,64,128, 15,10]，其中包括9个输入，10个输出中包括原9个输入加1个目标变量 (MAIAC AOD或OMI-NO2)，选取的训练超参数：批学习样本大小为1024， 初始的学习率为0.1，采用了Adam的梯度下降法优化方法。
步骤五、交通变量的提取，用于提取的道路包括国道及主高速公路，其他 次要的道路未包括在内。提取最近的距离采用基于KD树实现最近邻的快速检 索，采用Python的scipy包的cKDTree实现快速的交通变量的最近邻提取并计算 最近邻的距离。统计一定缓冲距离内的国道或高速公路的道路长度，采用了R 统计软件的包rgeos的gbuffer函数及raster包的Intersect联合实现一定距离缓冲区 国道及主干道路的长度。采用1公里到10公里(间隔1公里)的敏感性分析，得 到当缓冲区距离取10公里时取得道路长度同地面NO2最大的相关系数0.51，因 此取10公里作为最佳的缓冲区距离。
步骤六、土地利用变量的提取，选取了居民地、工业用地、商业用地及垃 圾处理场4类用地，计算一定缓冲范围内这些用地面积占总面积的比例。采用 了类似步骤五中提取交通距离类似的函数，敏感性分析表明10km的缓冲距离可 取得这些土地利用同地面NO2最大的相关性0.47。
步骤七、社会经济及POI提取，社会经济即GDP及人口密度直接是1公里栅 格数据，只需要作简单坐标统一后叠加到协变量数据集中。而对POI，提取了 污水处理厂、垃圾处理厂、废旧玻璃处理厂及废纸处理场等的POI，采用 cKDTree计算目标点到POI的最短距离作为POI的协变量。
步骤八、时空变异变量提取，提取坐标数据，取每个栅格点的中心点的坐 标即x与y，以及其派生变量(x2,y2及xy)，根据坐标提取高程数据；对时间变 量提取年-天、月份、星期及年份4个多尺度时间变量信息。
步骤九、空气污染物浓度地表参数反演，由以上九个步骤获得时空一致的 数据集，建立注意力机制的全残差深度时空网络，如图2所示，其中输入变量 包括反演的高分气象参数4个、卫星变量2个、气象同化参数30个、交通变量2 个、土地利用变量1个、社会经济2个、POI变量1个、空间变量6个、多尺度时 间变量4个，总计52个输入变量，输出包括取log对数的PM2.5及NO2，所有数据 均经过standard scalar正则化处理，预测结果最后反正则化及取指数还原原数据 尺度得到预测结果。训练的网络模型采用14层，包括1个输入层、5个编码层、 1个特征表征层(中间层)、6个解码层及1个输出层，节点数从输入到输出依 次是[52,256,128,64,32,16,8,16,32,64,128,256,52,2]，其中包括52个输入，2个输出(地面的NO2及PM2.5)，选取的训练超参数：批学习样本大小为 1024，初始的学习率为0.1。对预测目标NO2及PM2.5进行了限制，即：0≤NO2≤300μg/m3及0≤PM2.≤2000μg/m3，带限制性条件的优化采用了投影梯度 下降法，采用了基于Tensorflow的软件包TensorFlow ConstrainedOptimization (TFCO)来实现。
步骤十一、精度验证及评估。初步的训练得到测试精度R2：0.84(NO2)及 0.89(PM2.5)；测试RMSE：8.3μg/m3(NO2)及22.34μg/m3(PM2.5)。同时采用 SHAP计算模型各个协变量贡献，计算得到贡献最大的前10个协变量包括 MAIAC AOD及OMI-NO2、3个来自MERRA2GMI变量(CO、NO2及PM2.5)、 交通道路长度、到POI最短距离、风速、坐标，这表明了包括污染源的 MERRE2 GMI的变量的主要贡献。
步骤十二、验证达标与否。总体上如步骤十一报告的NO2及PM2.5的R2及 RMSE基本达到要求。该方法用于预测京津塘地区2015年地面的NO2及PM2.5， 地面预测结果查看发现北部地区有6个点的NO2浓度的预测结果偏高，与其周边 预测结果差距较大，表明这些点预测存在偏差，转入步骤十三。
步骤十三、对于不合理或不达标的预测，调整预测协变量，超参数及限制 性条件循环训练。问题检查表明北部地区这6个点交通的协变量提取错误，对 此进行修正。返回步骤十重新训练。
步骤十四、训练得到模型精度几乎一样，异常点较少，对总的训练精度影 响不大。但修正后的结果更符合要求，最终将预测结果的栅格估计表面输出。 图3展示了2015年12月30日预测的NO2栅格图；图4展示了2015年12月30日预测 的PM2.5栅格图。
通过本实施例可知，本发明针对当前空气污染时空反演方法的主要缺点， 提出了相应的解决办法，具体有：
1)用于建模的数据大部分基于局部区域或时段，样本时空覆盖度不足。 过小的采样域会对总体估计值造成偏差。对此本专利数据的准备覆盖了中国大 陆地区，取多年的天数据，所以本专利具有更为宽泛的时空覆盖度，用于建模 的数据样本更为全面，更能体现总体的分布情况。
2)缺乏高分辨率的气象协变量数据。气象协变量数据是影响空气污染的 关键变量，但全国范围内缺乏可靠的高分辨率栅格气象数据，现有的研究大部 分采用了分辨率较粗的气象同化数据，难以反映在精细尺度上要素的变异情 况，这严重影响空气污染反演效果。本专利结合国家气象局地面实测数据，对 气象变量进行高分辨率及高精度的反演，获得了可靠的结果。
3)现有的方法采用了CTM等及卫星气溶胶数据等提取空气污染的分布， 但是空气污染源于多种不同的来源，现有的大部分方法难以全面的包括各种污 染来源数据，是导致估计有偏的原因之一。本专利考虑得更全面，除了常规的 气溶胶AOD数据，也包括了最新的MERRA2 GMI的多达26种污染源的同化数 据，提取了地面交通、土地利用及POI(Point ofInterest)的多个潜在污染来源 的协变量数据，从而使得本专利的影响要素考虑更为全面。
4)现有的方法许多没有处理卫星遥感气溶胶及OMI NO2数据大量缺失 值，或者由于气象粗分辨率协变量及建模方法的原因反演缺失值的精度有限。 为此本专利提出了基于全残差深度学习的方法将缺失的气溶胶及OMI NO2缺失 数据插补完整，使得研究区域内无缺失值，从而使得本专利的评估范围更为宽 泛，不受到卫星数据缺失值的影响。
5)现有的方法采用了GAM、传统的多层感知机及克里格等，这些传统方 法学习能力有限，虽然有的也采用了深层CNN(convolutional neural network)、LSTM(long short-term memory)及支持向量机等，但支持向量机 需要复杂的特征提取，效率较低，而CNN本身受到深层网络导致梯度消失的影 响等。本专利在建模方法上采用了最新的融入了注意力机制的全残差深度网络 模型，大幅提高学习测试精度，而连续变量的输入保全输入信息，在实际中提 高了其泛化性，是本专利成果实施的基础模型。
6)现有方法大部分对结果大部分采用测试的精度，缺乏对实际结果的有 效性验证及解译，对此本专利采用有效解译方法，解译每个预测变量对结果的 贡献，并建立了回调循环机制，通过多次的调整训练获得更可信结果。
上述实施方式并非是对本发明的限制，本发明也并不仅限于上述举例，本 技术领域的技术人员在本发明的技术方案范围内所做出的变化、改型、添加或 替换，也均属于本发明的保护范围。

Claims (10)

1.一种多源时空大数据深度融合的空气污染预测方法，其特征在于：包括以下步骤：
步骤一、采集多源大数据；
步骤二、对步骤一采集到的数据进行预处理；
步骤三、高时空分辨率气象数据的插补，反演地面统一坐标的地面气象参数；
步骤四、气溶胶参数、NO2遥感参数缺失反演及升尺度；
步骤五、交通变量、土地利用变量、社会经济及POI变量、时空变异变量的提取；
步骤六、将各类时空大数据的协变量数据进行时空融合，形成统一尺度及空间坐标的数据集；
步骤七、反演空气污染浓度地表参数；
步骤八、精度验证及评估；
步骤九、验证达标与否；如若达标，进入步骤十一；如若不达标，进入步骤十；
步骤十、对于不合理或不达标预测，调整预测协变量、超参数及限制性条件循环训练，直至得到合理的模型及预测；
步骤十一、结果输出。
2.根据权利要求1所述的多源时空大数据深度融合的空气污染预测方法，其特征在于：步骤二中，数据预处理包括对数据进行质量控制、时空融合及格式转换；质量控制是按照有效数据阈值，或者质量标记对数据进行清理，删除无效数据；时空融合是对多源大数据在时空上进行融合，取得统一的目标时空分辨率；格式转换是将格式不一致的数据，转换成统一的高分空间栅格数据。
3.根据权利要求1所述的多源时空大数据深度融合的空气污染预测方法，其特征在于：步骤三中，采用一体化的插值方法进行高时空分辨率气象数据的插补，汇总气象数据中的特征协变量，建立目标变量统一且引入注意力层的全残差深度网络模型；对各特征协变量加权，设计多变量输出层，采用多个目标变量同时输出，便于模型参数间共享；在建模过程中，令特征的注意力权重输入为：C＝{ci}，并由此对输入的特征进行加权：
Figure FDA0002929512600000021
其中，x＝{xi}为特征值矢量输入，xi则为第i个特征，c＝{ci}为特征值矢量的注意力权重，ci为第i个特征的注意力权重系数；
Figure FDA0002929512600000022
为被注意力权重加权之后的特征输出，°代表矩阵元素相乘；F为特征的个数；
采用softmax层来实现注意力层，以满足式2的限制条件：
Figure FDA0002929512600000023
其中，ei为注意力层的第i单元的输入，ci为第i个单元的权重输出，exp(…)为取指数函数；j为特征的索引，ej为注意力层的j索引的输入，F为特征的个数；注意力的权重系数通过模型优化自适应求解；
多变量输出的损失函数为：
Figure FDA0002929512600000024
其中，N为样本数，
Figure FDA0002929512600000025
代表空气温度的损失，
Figure FDA0002929512600000026
代表相对湿度的损失，
Figure FDA0002929512600000027
代表空气压强的损失，
Figure FDA0002929512600000028
代表风速的损失；yt、yr、yp及yw分别代表空气温度、相对湿度、压强及风速的正则化后的观察值，而
Figure FDA0002929512600000029
Figure FDA00029295126000000210
及
Figure FDA00029295126000000211
分别代表空气温度、相对湿度、压强及风速的估计值；θW,b为网络模型的需要优化的参数集合，Ω(θW,b)则为参数集θW,b的正则化项，采用弹性网络进行正则化；
最后，对目标变量结果进行反正则化，获得同原数据尺度一致的估计值。
4.根据权利要求3所述的多源时空大数据深度融合的空气污染预测方法，其特征在于：引入全残差深度网络模型为基模型的聚集引导聚集建模，训练多个基模型，计算各参数的均值及变差，获得更优的地面气象参数反演结果。
5.根据权利要求1所述的多源时空大数据深度融合的空气污染预测方法，其特征在于：步骤四中，气溶胶参数指的是MAIACAOD，NO2遥感参数指的是OMI-NO2，在步骤三生成的高分气象参数基础上，分别建立全残差深度网络的AOD及OMI-NO2的缺值反演模型；在缺值反演模型中引入注意力层，提升重要协变量权重，以输入协变量+缺值需要插补的变量共同作为输出的目标变量；
再分别对MAIACAOD及OMI-NO2建立缺值反演及升尺度模型，采用MSE作为目标损失函数：
Figure FDA0002929512600000031
其中，N为样本数，lMSE表示采用均方误差的损失函数，y代表了地面空气污染监测站点的观察的NO2浓度，
Figure FDA0002929512600000032
为地面高分辨率的尺度转换后得到的地面的OMI-NO2的代理变量，引入了坡度因子s及截距a，通过其转换成地面NO2，使其同地面之间的观察值误差最低，从而优化参数得到高分辨率的尺度化结果，对于缺失值直接采用尺度转化模型输入协变量即可得到缺失值的估计。
6.根据权利要求1所述的多源时空大数据深度融合的空气污染预测方法，其特征在于：步骤五中，交通变量的提取采用最近邻提取或缓冲区分析法；缓冲区分析法计算缓冲区内主干道的长度，通过敏感性分析，选取最优缓冲距离作为最后距离：
Figure FDA0002929512600000041
其中，d*代表最优缓冲距离，r为相关系数函数，rlen(di)代表以di为缓冲距离时统计缓冲区内的道路长度，x则为目标点空气污染物的浓度，i为测试的缓冲距离索引；
土地利用变量的提取通过计算一定缓冲区范围内土地利用所占的面积比例，将此比例作为污染来源的代理变量；作敏感性分析，确定最优的缓冲距离：
Figure FDA0002929512600000042
其中，d*代表最优缓冲距离，r为相关系数函数，lprop(di)代表以di为缓冲距离时统计缓冲区内的土地利用面积比例，x则为目标点空气污染物的浓度，i为测试的缓冲距离索引；
社会经济指人口密度及GDP，通过空间叠加分析获得样本相应的人口密度及GDP；POI数据则先选择污染源相关的兴趣点，计算到最近的POI的距离；
时空变异变量的提取包括空间变量、时间变量的提取两方面，其中，空间变量包括高程、坐标及其派生变量，用于捕捉空间变异信息；时间变量包括多尺度时间信息，包括年-天、月份、星期、及年份。
7.根据权利要求1所述的多源时空大数据深度融合的空气污染预测方法，其特征在于：步骤八中，基模型采用引入注意力层的全残差深度网络模型，输入变量包括步骤一至步骤五获得的参数数据，设输入变量个数为D个，注意力节点也为D个，对关键变量进行加权，输出变量为M个，即需要估计的变量个数；同时，对输出设定了限制性条件，通过限制性优化方法，取得符合限制性条件的解；输出的损失函数采用下式：
Figure FDA0002929512600000051
Figure FDA0002929512600000052
Figure FDA0002929512600000053
其中，L(θW,b)代表了网络参数集θW,b的损失函数，N为样本数，
Figure FDA0002929512600000054
及
Figure FDA0002929512600000055
分别代表NO2及PM2.5正则化后的观察值，而
Figure FDA0002929512600000056
及
Figure FDA0002929512600000057
分别代表了针对输入协变量集x的NO2及PM2.5的网络预测值，
Figure FDA0002929512600000058
及
Figure FDA0002929512600000059
分别代表NO2及PM2.5的MSE的损失函数，Ω(θW,b)则为参数集θW,b的正则化项，同理，采用弹性网络实现正则化；式9及式10定义了NO2及PM2.5的网络预测值需要满足的条件，其中，
Figure FDA00029295126000000510
及
Figure FDA00029295126000000511
分别表示NO2的取值的最小值与最大值，
Figure FDA00029295126000000512
及
Figure FDA00029295126000000513
分别表示PM2.5的取值的最小值与最大值；
将二污染物合在一起在一个网络中输出，采用bootstrap重抽样方法，对样本及特征进行多次重复抽样，得到不同的样本及特征变量集，对这些不同的样本集合对全残差深度模型结构进行样本数随机化，以获得差异较大的基模型，将这些样本分别训练这些基模型，最后将训练的模型分别预测，预测结果反正则化得到NO2与PM2.5的估计浓度，统计多个模型预测浓度得到估计平均值及其变差。
8.根据权利要求1所述的多源时空大数据深度融合的空气污染预测方法，其特征在于：步骤八中，采用独立性验证方法对聚集引导得到的结果进行精度验证及结果评估。
9.根据权利要求1所述的多源时空大数据深度融合的空气污染预测方法，其特征在于：步骤九中，验证总体的预测目标是否符合要求，验证各个协变量对总的预测解释是否合理。
10.根据权利要求1所述的多源时空大数据深度融合的空气污染预测方法，其特征在于：步骤十一中，对于得到的合理的优化的训练模型及超参数，保存相应的模型及参数，并将合理的预测结果输出，供下一步的应用使用。
CN202110144010.2A 2021-02-02 2021-02-02 一种多源时空大数据深度融合的空气污染预测方法 Active CN112905560B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110144010.2A CN112905560B (zh) 2021-02-02 2021-02-02 一种多源时空大数据深度融合的空气污染预测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110144010.2A CN112905560B (zh) 2021-02-02 2021-02-02 一种多源时空大数据深度融合的空气污染预测方法

Publications (2)

Publication Number Publication Date
CN112905560A true CN112905560A (zh) 2021-06-04
CN112905560B CN112905560B (zh) 2022-10-11

Family

ID=76121561

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110144010.2A Active CN112905560B (zh) 2021-02-02 2021-02-02 一种多源时空大数据深度融合的空气污染预测方法

Country Status (1)

Country Link
CN (1) CN112905560B (zh)

Cited By (22)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113297528A (zh) * 2021-06-10 2021-08-24 四川大学 一种基于多源大数据的no2高分辨率时空分布计算方法
CN113408918A (zh) * 2021-06-28 2021-09-17 哈尔滨工业大学 基于多时相遥感分析的村镇垃圾降尺度时空分布反演方法
CN113434958A (zh) * 2021-06-25 2021-09-24 南京工程学院 基于异构集成代理模型的磁悬浮飞轮电机多目标优化方法
CN113535510A (zh) * 2021-06-24 2021-10-22 北京理工大学 一种大规模数据中心数据采集的自适应抽样模型优化方法
CN114609008A (zh) * 2022-03-30 2022-06-10 中国科学院地理科学与资源研究所 一种气体浓度预测模型的训练方法和气体浓度预测方法
CN114676121A (zh) * 2022-03-30 2022-06-28 重庆大学 基于多尺度时空记忆共享网络的缺失监测数据填补方法
CN114782871A (zh) * 2022-04-29 2022-07-22 广东技术师范大学 一种基于物联网的海洋异常信息监测方法和装置
CN114896826A (zh) * 2022-07-13 2022-08-12 成都信息工程大学 基于物理学结合残差注意力网络的行星边界层参数化方法
CN114926749A (zh) * 2022-07-22 2022-08-19 山东大学 基于遥感图像的近地面大气污染物反演方法及系统
CN115345075A (zh) * 2022-08-17 2022-11-15 北京城市气象研究院 一体化气溶胶污染气象指数-气溶胶浓度估算方法及系统
CN116596158A (zh) * 2023-06-14 2023-08-15 深圳市汉宇环境科技有限公司 区域污染源排放总量预测方法
CN116910457A (zh) * 2023-09-13 2023-10-20 北京师范大学 一种基于区域的污染物反演方法和装置
CN117113012A (zh) * 2023-06-15 2023-11-24 中国水利水电第十二工程局有限公司 一种基于等级控制的平差数据预处理系统
CN117129638A (zh) * 2023-10-26 2023-11-28 江西怡杉环保股份有限公司 一种区域空气环境质量监测方法及系统
CN117216480A (zh) * 2023-09-18 2023-12-12 宁波大学 一种深度耦合地理时空信息的近地表臭氧遥感估算方法
CN117390008A (zh) * 2023-12-11 2024-01-12 北京星球空天信息技术有限公司 多类型观测仪器的测量数据处理方法和装置
CN117592769A (zh) * 2024-01-19 2024-02-23 四川绿豆芽信息技术有限公司 一种碳小屋站点管理方法及系统
CN117592005A (zh) * 2024-01-19 2024-02-23 中国科学院空天信息创新研究院 Pm2.5浓度卫星遥感估算方法、装置、设备及介质
CN118094112A (zh) * 2024-04-18 2024-05-28 北京市农林科学院智能装备技术研究中心 灌溉数据同化方法、装置、电子设备及存储介质
CN118196151A (zh) * 2024-03-13 2024-06-14 中国人民解放军国防科技大学 一种基于气象敏感性分析的光谱反演参数预测方法及装置
CN118380078A (zh) * 2024-06-26 2024-07-23 山东建筑大学 一种基于多任务学习的pm2.5和o3污染贡献定量估算方法
CN118196151B (zh) * 2024-03-13 2024-09-24 中国人民解放军国防科技大学 一种基于气象敏感性分析的光谱反演参数预测方法及装置

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107798425A (zh) * 2017-10-16 2018-03-13 中国科学院地理科学与资源研究所 一种基于大数据的时空混淆暴露度评估系统及方法
US20190057171A1 (en) * 2017-02-10 2019-02-21 Guangxi Institute of Meteorlogical Disaster- Reducing Research Novel nonlinear method for area-wide near surface air temperature precision retrieval
CN110186820A (zh) * 2018-12-19 2019-08-30 河北中科遥感信息技术有限公司 多源数据融合与环境污染源及污染物分布分析方法
CN110287455A (zh) * 2019-05-28 2019-09-27 武汉大学 一种结合遥感数据与社会感知数据的pm2.5深度学习反演方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20190057171A1 (en) * 2017-02-10 2019-02-21 Guangxi Institute of Meteorlogical Disaster- Reducing Research Novel nonlinear method for area-wide near surface air temperature precision retrieval
CN107798425A (zh) * 2017-10-16 2018-03-13 中国科学院地理科学与资源研究所 一种基于大数据的时空混淆暴露度评估系统及方法
CN110186820A (zh) * 2018-12-19 2019-08-30 河北中科遥感信息技术有限公司 多源数据融合与环境污染源及污染物分布分析方法
CN110287455A (zh) * 2019-05-28 2019-09-27 武汉大学 一种结合遥感数据与社会感知数据的pm2.5深度学习反演方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
吴健生等: "基于AOD数据反演地面PM_(2.5)浓度研究进展", 《环境科学与技术》 *
杨立娟等: "MODIS卫星遥感估计福州地区近地面PM_(2.5)浓度", 《遥感学报》 *

Cited By (35)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113297528B (zh) * 2021-06-10 2022-07-01 四川大学 一种基于多源大数据的no2高分辨率时空分布计算方法
CN113297528A (zh) * 2021-06-10 2021-08-24 四川大学 一种基于多源大数据的no2高分辨率时空分布计算方法
CN113535510A (zh) * 2021-06-24 2021-10-22 北京理工大学 一种大规模数据中心数据采集的自适应抽样模型优化方法
CN113535510B (zh) * 2021-06-24 2024-01-26 北京理工大学 一种大规模数据中心数据采集的自适应抽样模型优化方法
CN113434958B (zh) * 2021-06-25 2023-09-01 南京工程学院 基于异构集成代理模型的磁悬浮飞轮电机多目标优化方法
CN113434958A (zh) * 2021-06-25 2021-09-24 南京工程学院 基于异构集成代理模型的磁悬浮飞轮电机多目标优化方法
CN113408918A (zh) * 2021-06-28 2021-09-17 哈尔滨工业大学 基于多时相遥感分析的村镇垃圾降尺度时空分布反演方法
CN113408918B (zh) * 2021-06-28 2022-04-19 哈尔滨工业大学 基于多时相遥感分析的村镇垃圾降尺度时空分布反演方法
CN114609008A (zh) * 2022-03-30 2022-06-10 中国科学院地理科学与资源研究所 一种气体浓度预测模型的训练方法和气体浓度预测方法
CN114676121A (zh) * 2022-03-30 2022-06-28 重庆大学 基于多尺度时空记忆共享网络的缺失监测数据填补方法
CN114609008B (zh) * 2022-03-30 2023-09-01 中国科学院地理科学与资源研究所 一种气体浓度预测模型的训练方法和气体浓度预测方法
CN114782871A (zh) * 2022-04-29 2022-07-22 广东技术师范大学 一种基于物联网的海洋异常信息监测方法和装置
CN114896826B (zh) * 2022-07-13 2022-09-09 成都信息工程大学 基于物理学结合残差注意力网络的行星边界层参数化方法
CN114896826A (zh) * 2022-07-13 2022-08-12 成都信息工程大学 基于物理学结合残差注意力网络的行星边界层参数化方法
CN114926749B (zh) * 2022-07-22 2022-11-04 山东大学 基于遥感图像的近地面大气污染物反演方法及系统
CN114926749A (zh) * 2022-07-22 2022-08-19 山东大学 基于遥感图像的近地面大气污染物反演方法及系统
CN115345075A (zh) * 2022-08-17 2022-11-15 北京城市气象研究院 一体化气溶胶污染气象指数-气溶胶浓度估算方法及系统
CN115345075B (zh) * 2022-08-17 2023-04-18 北京城市气象研究院 一体化气溶胶污染气象指数-气溶胶浓度估算方法及系统
CN116596158A (zh) * 2023-06-14 2023-08-15 深圳市汉宇环境科技有限公司 区域污染源排放总量预测方法
CN117113012A (zh) * 2023-06-15 2023-11-24 中国水利水电第十二工程局有限公司 一种基于等级控制的平差数据预处理系统
CN116910457A (zh) * 2023-09-13 2023-10-20 北京师范大学 一种基于区域的污染物反演方法和装置
CN116910457B (zh) * 2023-09-13 2023-12-08 北京师范大学 一种基于区域的污染物反演方法和装置
CN117216480A (zh) * 2023-09-18 2023-12-12 宁波大学 一种深度耦合地理时空信息的近地表臭氧遥感估算方法
CN117129638B (zh) * 2023-10-26 2024-01-12 江西怡杉环保股份有限公司 一种区域空气环境质量监测方法及系统
CN117129638A (zh) * 2023-10-26 2023-11-28 江西怡杉环保股份有限公司 一种区域空气环境质量监测方法及系统
CN117390008A (zh) * 2023-12-11 2024-01-12 北京星球空天信息技术有限公司 多类型观测仪器的测量数据处理方法和装置
CN117390008B (zh) * 2023-12-11 2024-04-12 北京星球空天信息技术有限公司 多类型观测仪器的测量数据处理方法和装置
CN117592769A (zh) * 2024-01-19 2024-02-23 四川绿豆芽信息技术有限公司 一种碳小屋站点管理方法及系统
CN117592005A (zh) * 2024-01-19 2024-02-23 中国科学院空天信息创新研究院 Pm2.5浓度卫星遥感估算方法、装置、设备及介质
CN117592769B (zh) * 2024-01-19 2024-04-05 四川绿豆芽信息技术有限公司 一种碳小屋站点管理方法及系统
CN117592005B (zh) * 2024-01-19 2024-04-26 中国科学院空天信息创新研究院 Pm2.5浓度卫星遥感估算方法、装置、设备及介质
CN118196151A (zh) * 2024-03-13 2024-06-14 中国人民解放军国防科技大学 一种基于气象敏感性分析的光谱反演参数预测方法及装置
CN118196151B (zh) * 2024-03-13 2024-09-24 中国人民解放军国防科技大学 一种基于气象敏感性分析的光谱反演参数预测方法及装置
CN118094112A (zh) * 2024-04-18 2024-05-28 北京市农林科学院智能装备技术研究中心 灌溉数据同化方法、装置、电子设备及存储介质
CN118380078A (zh) * 2024-06-26 2024-07-23 山东建筑大学 一种基于多任务学习的pm2.5和o3污染贡献定量估算方法

Also Published As

Publication number Publication date
CN112905560B (zh) 2022-10-11

Similar Documents

Publication Publication Date Title
CN112905560B (zh) 一种多源时空大数据深度融合的空气污染预测方法
CN113297528B (zh) 一种基于多源大数据的no2高分辨率时空分布计算方法
Zhao et al. High-resolution daily AOD estimated to full coverage using the random forest model approach in the Beijing-Tianjin-Hebei region
CN112884079A (zh) 一种基于Stacking集成模型的近地面二氧化氮浓度估算方法
CN112699959B (zh) 基于能量泛函模型的多源多尺度降水数据融合方法和装置
CN107273995A (zh) 空气质量预报方法
Heo et al. Digital elevation model-based convolutional neural network modeling for searching of high solar energy regions
Kong et al. High-resolution (0.05× 0.05) NO x emissions in the Yangtze River Delta inferred from OMI
Ma et al. Statistical spatial-temporal modeling of ambient ozone exposure for environmental epidemiology studies: A review
CN114254802B (zh) 气候变化驱动下植被覆盖时空变化的预测方法
Bagheri A machine learning-based framework for high resolution mapping of PM2. 5 in Tehran, Iran, using MAIAC AOD data
CN114004163B (zh) 一种基于modis和长短时记忆网络模型的pm2.5反演方法
Yu et al. Deep learning-based downscaling of tropospheric nitrogen dioxide using ground-level and satellite observations
Jiang et al. A Municipal PM2. 5 Forecasting Method Based on Random Forest and WRF Model.
CN113011455B (zh) 一种空气质量预测svm模型构建方法
CN112329265A (zh) 卫星遥感降水精细化空间估算方法及系统
CN115544889A (zh) 基于深度学习的数值模式降水偏差校正方法
CN116223395A (zh) 近地表痕量气体浓度反演模型及反演方法
CN116340863A (zh) 空气污染物预测方法、装置、电子设备及可读存储介质
Li et al. Combing GOME-2B and OMI Satellite Data to Estimate Near-Surface NO 2 of Mainland China
Fan et al. Simulation of average monthly ozone exposure concentrations in China: A temporal and spatial estimation method
Sun et al. Construction of a spatially gridded heat flux map based on airborne flux Measurements using remote sensing and machine learning methods
CN115267066A (zh) 基于卫星观测污染气体浓度的化石燃料二氧化碳排放计算方法
CN112905561B (zh) 一种基于深度双模态的气象参数精细尺度转化方法
Hu et al. An Interpolation and Prediction Algorithm for XCO2 Based on Multi-Source Time Series Data

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