CN112905560B - 一种多源时空大数据深度融合的空气污染预测方法 - Google Patents
一种多源时空大数据深度融合的空气污染预测方法 Download PDFInfo
- Publication number
- CN112905560B CN112905560B CN202110144010.2A CN202110144010A CN112905560B CN 112905560 B CN112905560 B CN 112905560B CN 202110144010 A CN202110144010 A CN 202110144010A CN 112905560 B CN112905560 B CN 112905560B
- Authority
- CN
- China
- Prior art keywords
- data
- variables
- 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.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 85
- 238000003915 air pollution Methods 0.000 title claims abstract description 52
- 230000004927 fusion Effects 0.000 title claims abstract description 29
- 238000012549 training Methods 0.000 claims abstract description 28
- 239000000443 aerosol Substances 0.000 claims abstract description 26
- 238000005457 optimization Methods 0.000 claims abstract description 17
- 238000012795 verification Methods 0.000 claims abstract description 15
- 125000004122 cyclic group Chemical group 0.000 claims abstract description 5
- 239000000872 buffer Substances 0.000 claims description 44
- 238000000605 extraction Methods 0.000 claims description 26
- 239000000809 air pollutant Substances 0.000 claims description 20
- 231100001243 air pollutant Toxicity 0.000 claims description 20
- 238000012544 monitoring process Methods 0.000 claims description 12
- 239000003344 environmental pollutant Substances 0.000 claims description 10
- 231100000719 pollutant Toxicity 0.000 claims description 10
- 238000004458 analytical method Methods 0.000 claims description 9
- 238000006243 chemical reaction Methods 0.000 claims description 9
- 238000011156 evaluation Methods 0.000 claims description 9
- 238000010206 sensitivity analysis Methods 0.000 claims description 8
- 238000007781 pre-processing Methods 0.000 claims description 7
- 238000003908 quality control method Methods 0.000 claims description 7
- 230000002776 aggregation Effects 0.000 claims description 6
- 238000004220 aggregation Methods 0.000 claims description 6
- 238000012952 Resampling Methods 0.000 claims description 5
- 230000008569 process Effects 0.000 claims description 5
- 239000011159 matrix material Substances 0.000 claims description 3
- 230000003139 buffering effect Effects 0.000 claims 4
- 239000004480 active ingredient Substances 0.000 claims 1
- 239000000203 mixture Substances 0.000 claims 1
- 238000012360 testing method Methods 0.000 abstract description 17
- 230000007246 mechanism Effects 0.000 abstract description 10
- 238000005516 engineering process Methods 0.000 abstract description 9
- 230000006870 function Effects 0.000 description 23
- 230000002123 temporal effect Effects 0.000 description 22
- 238000009826 distribution Methods 0.000 description 9
- 238000005259 measurement Methods 0.000 description 8
- CBENFWSGALASAD-UHFFFAOYSA-N Ozone Chemical compound [O-][O+]=O CBENFWSGALASAD-UHFFFAOYSA-N 0.000 description 5
- 238000012937 correction Methods 0.000 description 5
- 238000013135 deep learning Methods 0.000 description 5
- 230000000694 effects Effects 0.000 description 5
- JCXJVPUVTGWSNB-UHFFFAOYSA-N nitrogen dioxide Inorganic materials O=[N]=O JCXJVPUVTGWSNB-UHFFFAOYSA-N 0.000 description 5
- 238000013527 convolutional neural network Methods 0.000 description 4
- 238000010586 diagram Methods 0.000 description 4
- 230000008034 disappearance Effects 0.000 description 3
- 230000007613 environmental effect Effects 0.000 description 3
- 239000000284 extract Substances 0.000 description 3
- 238000011160 research Methods 0.000 description 3
- KRQUFUKTQHISJB-YYADALCUSA-N 2-[(E)-N-[2-(4-chlorophenoxy)propoxy]-C-propylcarbonimidoyl]-3-hydroxy-5-(thian-3-yl)cyclohex-2-en-1-one Chemical compound CCC\C(=N/OCC(C)OC1=CC=C(Cl)C=C1)C1=C(O)CC(CC1=O)C1CCCSC1 KRQUFUKTQHISJB-YYADALCUSA-N 0.000 description 2
- ORILYTVJVMAKLC-UHFFFAOYSA-N Adamantane Natural products C1C(C2)CC3CC1CC2C3 ORILYTVJVMAKLC-UHFFFAOYSA-N 0.000 description 2
- RAHZWNYVWXNFOC-UHFFFAOYSA-N Sulphur dioxide Chemical compound O=S=O RAHZWNYVWXNFOC-UHFFFAOYSA-N 0.000 description 2
- 229910002091 carbon monoxide Inorganic materials 0.000 description 2
- 238000013461 design Methods 0.000 description 2
- 238000011478 gradient descent method Methods 0.000 description 2
- 230000009931 harmful effect Effects 0.000 description 2
- 230000003287 optical effect Effects 0.000 description 2
- 239000013618 particulate matter Substances 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 238000012950 reanalysis Methods 0.000 description 2
- 238000005070 sampling Methods 0.000 description 2
- 230000001932 seasonal effect Effects 0.000 description 2
- 238000012706 support-vector machine Methods 0.000 description 2
- 230000009466 transformation Effects 0.000 description 2
- 239000002699 waste material Substances 0.000 description 2
- MGWGWNFMUOTEHG-UHFFFAOYSA-N 4-(3,5-dimethylphenyl)-1,3-thiazol-2-amine Chemical compound CC1=CC(C)=CC(C=2N=C(N)SC=2)=C1 MGWGWNFMUOTEHG-UHFFFAOYSA-N 0.000 description 1
- 238000012935 Averaging Methods 0.000 description 1
- UGFAIRIUMAVXCW-UHFFFAOYSA-N Carbon monoxide Chemical compound [O+]#[C-] UGFAIRIUMAVXCW-UHFFFAOYSA-N 0.000 description 1
- 206010035664 Pneumonia Diseases 0.000 description 1
- 230000002159 abnormal effect Effects 0.000 description 1
- 230000003044 adaptive effect Effects 0.000 description 1
- 238000007792 addition Methods 0.000 description 1
- 239000000654 additive Substances 0.000 description 1
- 230000000996 additive effect Effects 0.000 description 1
- 238000013528 artificial neural network Methods 0.000 description 1
- 208000006673 asthma Diseases 0.000 description 1
- 239000008277 atmospheric particulate matter Substances 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 235000000332 black box Nutrition 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000013480 data collection Methods 0.000 description 1
- 230000007423 decrease Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 239000011521 glass Substances 0.000 description 1
- 230000036541 health Effects 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 238000012417 linear regression Methods 0.000 description 1
- 230000007774 longterm Effects 0.000 description 1
- 238000010801 machine learning Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 239000010813 municipal solid waste Substances 0.000 description 1
- 238000010606 normalization Methods 0.000 description 1
- 239000010893 paper waste Substances 0.000 description 1
- 238000002360 preparation method Methods 0.000 description 1
- 238000002310 reflectometry Methods 0.000 description 1
- 210000002345 respiratory system Anatomy 0.000 description 1
- 239000010865 sewage Substances 0.000 description 1
- 230000006403 short-term memory Effects 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 238000012732 spatial analysis Methods 0.000 description 1
- 229940044609 sulfur dioxide Drugs 0.000 description 1
- 235000010269 sulphur dioxide Nutrition 0.000 description 1
- 238000012876 topography Methods 0.000 description 1
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-resolutionImaging Spectroradiometer)卫星传感器提取气溶胶光学深度系数(aerosol opticaldepth，简称AOD)，可以从AOD反演PM10(fine particulate matter with aerodynamicdiameter smaller than 10μm)及PM2.5(fine particulate matter with aerodynamicdiameter smaller than 2.5μm)浓度。而从Ozone Monitoring Instrument(OMI)传感器Aura提取的O3(ozone)及NO2(nitrogen dioxide)，可用于反演地面空气污染物O3及NO2的时空变化情况。但卫星数据最大的问题是由于云覆盖、地面的高反射率及传感器本身的缺陷，会导致有大量的缺失值产生，从而严重影响卫星产品的应用。
为解决卫星数据的缺失问题，现有技术设计了基于概化累加模型(generalizedadditive model，简称GAM)的缺值插补的非线性方法，采用基于GEOS 5-FP的气象数据、土地利用、CMAQ(Community Multiscale Air Quality)模拟数据等插补完整中国长三角地区的卫星气溶胶数据MAIAC(Multi-Angle Implementation of Atmospheric Correction)AOD，然后用于估算地面PM2.5的浓度；也有采用了类似的变量，通过普通的前馈神经网络进行缺值插补。其他的一些方法还包括了采用均值、最近邻或其他粗分辨率如CMAQ的模拟AOD代替缺失值。虽然现有技术已采用了非线性方法插补缺失值，比简单的替代方法要好，但是他们的缺值插补方法是基于粗分辨率的气象数据(分辨率25km x 25km)，同反演的地面空气污染浓度的分辨率(1x1km2)差距过大，这时导致估计结果偏差原因之一。
由于影响空气污染的高分辨率的重要气象参数(空气温度、相对湿度、风速及气压等)缺乏，将会导致在缺乏高分辨率气象参数情况下采用卫星与/或其他模拟的气象变量反演地面空气污染物浓度的估计偏差。现有技术提出了融合多源遥感信息与气候环境的综合性气象环境评估方法，其方法基于较粗的卫星分辨率，且只是反演了季节性的空气污染浓度分布并进行评估，缺乏高时间分辨率的结果。现有技术还提出了一种大气颗粒物星地综合定量遥感融合的反演方法，该方法没有对卫星气溶胶的缺失进行有效的处理，会导致相当多的时空点由于卫星数据的缺乏无法进行反演；而反演采用的气象等参数分辨率也较粗。对于采用多源遥感数据融合，基于GEOS-FP(Goddard Earth Observing System-Forward Processing)系列数反演PM2.5污染物，同样其用于反演的地面分辨率较粗，难以反映地面的在细尺度下的变化情况。以及现有的深度森林算法估计城市模型估算，这种方法会受到基于树学习模型的输入离散数据的限制，在样本数较少时会导致表面建模不连续性。虽然目前已有采用深度学习的CNN等进行空气污染的地表参数重构，但由于空气污染同影响因素的复杂非线性关系，难以采用卷积网取得理想的结果，且过深的网络会导致梯度消失问题，影响最后的估计结果。现有技术还提出基于空气质量数据与图像实现多源异构的融合的PM2.5预测模型，该方法需要采集照片估算PM2.5的浓度,评估结果受到环境散射光线的影响，估计精度有限。
发明内容
为了解决上述技术所存在的不足之处，本发明提供了一种多源时空大数据深度融合的空气污染预测方法。
为了解决以上技术问题，本发明采用的技术方案是：一种多源时空大数据深度融合的空气污染预测方法，包括以下步骤：
步骤一、采集多源大数据；
步骤二、对步骤一采集到的数据进行预处理；
步骤三、高时空分辨率气象数据的插补，反演地面统一坐标的地面气象参数；
步骤四、气溶胶参数、NO2遥感参数缺失反演及升尺度；
步骤五、交通变量、土地利用变量、社会经济及POI变量、时空变异变量的提取；
步骤六、将各类时空大数据的协变量数据进行时空融合，形成统一尺度及空间坐标的数据集；
步骤七、反演空气污染浓度地表参数；
步骤八、精度验证及评估；
步骤九、验证达标与否；如若达标，进入步骤十一；如若不达标，进入步骤十；
步骤十、对于不合理或不达标预测，调整预测协变量、超参数及限制性条件循环训练，直至得到合理的模型及预测；
步骤十一、结果输出。
进一步地，步骤二中，数据预处理包括对数据进行质量控制、时空融合及格式转换；质量控制是按照有效数据阈值，或者质量标记对数据进行清理，删除无效数据；时空融合是对多源大数据在时空上进行融合，取得统一的目标时空分辨率；格式转换是将格式不一致的数据，转换成统一的高分空间栅格数据。
进一步地，步骤三中，采用一体化的插值方法进行高时空分辨率气象数据的插补，汇总气象数据中的特征协变量，建立目标变量统一且引入注意力层的全残差深度网络模型；对各特征协变量加权，设计多变量输出层，采用多个目标变量同时输出，便于模型参数间共享；在建模过程中，令特征的注意力权重输入为：C＝{ci}，并由此对输入的特征进行加权：
其中，x＝{xi}为特征值矢量输入，xi则为第i个特征，c＝{ci}为特征值矢量的注意力权重，ci为第i个特征的注意力权重系数；
为被注意力权重加权之后的特征输出，°代表矩阵元素相乘；F为特征的个数；
采用softmax层来实现注意力层，以满足式2的限制条件：
其中，ei为注意力层的第i单元的输入，ci为第i个单元的权重输出，exp(…)为取指数函数；j为特征的索引，ej为注意力层的j索引的输入，F为特征的个数；注意力的权重系数通过模型优化自适应求解；
多变量输出的损失函数为：
其中，N为样本数，
代表空气温度的损失，
代表相对湿度的损失，
代表空气压强的损失，
代表风速的损失；yt、yr、yp及yw分别代表空气温度、相对湿度、压强及风速的正则化后的观察值，而
及
分别代表空气温度、相对湿度、压强及风速的估计值；θW,b为网络模型的需要优化的参数集合，Ω(θW,b)则为参数集θW,b的正则化项，采用弹性网络进行正则化；
最后，对目标变量结果进行反正则化，获得同原数据尺度一致的估计值。
进一步地，引入全残差深度网络模型为基模型的聚集引导聚集建模，训练多个基模型，计算各参数的均值及变差，获得更优的地面气象参数反演结果。
进一步地，步骤四中，气溶胶参数指的是MAIAC AOD，NO2遥感参数指的是OMI-NO2，在步骤三生成的高分气象参数基础上，分别建立全残差深度网络的AOD及OMI-NO2的缺值反演模型；在此缺值反演模型中引入注意力层，提升重要协变量权重，以输入协变量+缺值需要插补的变量共同作为输出的目标变量；
再分别对MAIACAOD及OMI-NO2建立缺值反演及升尺度模型，采用MSE作为目标损失函数：
其中，N为样本数，
表示采用均方误差的损失函数，y代表了地面空气污染监测站点的观察的NO2浓度，
为地面高分辨率的尺度转换后得到的地面的OMI-NO2的代理变量，引入了坡度因子s及截距a，通过其转换成地面NO2，使其同地面之间的观察值误差最低，从而优化参数得到高分辨率的尺度化结果，对于缺失值直接采用尺度转化模型输入协变量即可得到缺失值的估计。
进一步地，步骤五中，交通变量的提取采用最近邻提取或缓冲区分析法；缓冲区分析法计算缓冲区内主干道的长度，通过敏感性分析，选取最优缓冲距离作为最后距离：
其中，d*代表最优缓冲距离，r为相关系数函数，rlen(di)代表以di为缓冲距离时统计缓冲区内的道路长度，x则为目标点空气污染物的浓度，i为测试的缓冲距离索引；
土地利用变量的提取通过计算一定缓冲区范围内土地利用所占的面积比例，将此比例作为污染来源的代理变量；作敏感性分析，确定最优的缓冲距离：
其中，d*代表最优缓冲距离，r为相关系数函数，lprop(di)代表以di为缓冲距离时统计缓冲区内的土地利用面积比例，x则为目标点空气污染物的浓度，i为测试的缓冲距离索引；
社会经济指人口密度及GDP，通过空间叠加分析获得样本相应的人口密度及GDP；POI数据则先选择污染源相关的兴趣点，计算到最近的POI的距离；
时空变异变量的提取包括空间变量、时间变量的提取两方面，其中，空间变量包括高程、坐标及其派生变量，用于捕捉空间变异信息；时间变量包括多尺度时间信息，包括年-天、月份、星期、及年份。
进一步地，步骤七中，基模型采用引入注意力层的全残差深度网络模型，输入变量包括步骤一至步骤五获得的参数数据，设输入变量个数为D个，注意力节点也为D个，对关键变量进行加权，输出变量为M个，即需要估计的变量个数；同时，对输出设定了限制性条件，通过限制性优化方法，取得符合限制性条件的解；输出的损失函数采用下式：
其中，L(θW,b)代表了网络参数集θW,b的损失函数，N为样本数，
及
分别代表NO2及PM2.5正则化后的观察值，而
及
分别代表了针对输入协变量集x的NO2及PM2.5的网络预测值，
及
分别代表NO2及PM2.5的MSE的损失函数，Ω(θW,b)则为参数集θW,b的正则化项，同理，采用弹性网络实现正则化；式9及式10定义了NO2及PM2.5的网络预测值需要满足的条件，其中，
及
分别表示NO2的取值的最小值与最大值，
及
分别表示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种各种污染来源的MERRA2GMI(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同化数据拟合的污染物表面栅格相关变量列表
步骤二、对步骤一采集到的数据进行预处理；
数据预处理主要包括对数据进行质量控制、时空融合及格式转换等；质量控制是按照有效数据阈值，或者质量标记对数据进行清理，删除无效数据；采用统计标注，删除一些对建模无意义的孤立点，采用阈值定义：
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}，并由此对输入的特征进行加权：
其中，x＝{xi}为特征值矢量输入，xi则为第i个特征，c＝{ci}为特征值矢量的注意力权重，ci为第i个特征的注意力权重系数；
为被注意力权重加权之后的特征输出，°代表矩阵元素相乘；F为特征的个数；
采用softmax层来实现注意力层，以满足式2的限制条件：
其中，ei为注意力层的第i单元的输入，ci为第i个单元的权重输出，exp(…)为取指数函数；j为特征的索引，ej为注意力层的j索引的输入，ej与ei一样，均为注意力层的输入，不过j索引表示对所有的ej；F为特征的个数；注意力的权重系数通过模型优化自适应求解；
多变量输出的损失函数为：
其中，N为样本数，
代表空气温度的损失，
代表相对湿度的损失，
代表空气压强的损失，
代表风速的损失；yt、yr、yp及yw分别代表空气温度、相对湿度、压强及风速的正则化后的观察值，而
及
分别代表空气温度、相对湿度、压强及风速的估计值；θ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作为目标损失函数：
其中，N为样本数，
表示采用均方误差(mean square error，即MSE)的损失函数，y代表了地面空气污染监测站点的观察的NO2浓度，
为地面高分辨率的尺度转换后得到的地面的OMI-NO2的代理变量，引入了坡度因子s及截距a，通过其转换成地面NO2，使其同地面之间的观察值误差最低，从而优化参数得到高分辨率的尺度化结果，对于缺失值直接采用尺度转化模型输入协变量即可得到缺失值的估计。
步骤五、交通变量、土地利用变量、社会经济及POI变量、时空变异变量的提取；
交通变量提取：地面空气污染物如NO2同交通紧密相关，需要提取交通相关的变量。由于具体交通流量信息很难获取，本发明基于主要交通干道，包括各级高速公路、主要道路提取2个关键的交通要素：到最近交通主干道的距离、一定缓冲距离范围内道路的长度。对前者，采用最近邻提取；而对于后者，采用缓冲区分析方法，计算缓冲区内交通主干道的长度。通过敏感性分析，即计算一系列的距离缓冲区内道路长度同现有污染物NO2浓度间的相关性，选取最优缓冲距离作为最后距离：
其中，d*代表最优缓冲距离，r为相关系数函数，rlen(di)代表以di为缓冲距离时统计缓冲区内的道路长度，x则为目标点空气污染物的浓度，i为测试的缓冲距离索引。
土地利用变量提取：土地利用也提供了主要的排放源信息。本发明将污染来源相关的数据划分为一类，通过计算一定缓冲区范围内土地利用所占的面积比例，将此比例作为污染来源的代理变量；同交通变量提取类似，作敏感性分析，确定最优的缓冲距离：
其中，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。此处同时对输出设定了限制性条件，通过限制性优化方法，尽量取得符合限制性条件的解。输出的损失函数可以采用下式：
其中，L(θW,b)代表了网络参数集θW,b的损失函数，N为样本数，
及
分别代表NO2及PM2.5正则化后的观察值，而
及
分别代表了针对输入协变量集x的NO2及PM2.5的网络预测值，
及
分别代表NO2及PM2.5的MSE的损失函数，Ω(θW,b)则为参数集θW,b的正则化项，同理，采用弹性网络实现正则化；式9及式10定义了NO2及PM2.5的网络预测值需要满足的条件，即需要分别落在在区间
及
之内，其中
及
分别表示NO2的取值的最小值与最大值，
及
分别表示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°lon x0.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高分辨率的气象参数反演的测试精度
步骤四、气溶胶参数、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个来自MERRA2 GMI变量(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 (6)
1.一种多源时空大数据深度融合的空气污染预测方法，其特征在于：包括以下步骤：
步骤一、采集多源大数据；
步骤二、对步骤一采集到的数据进行预处理；
步骤三、高时空分辨率气象数据的插补，反演地面统一坐标的地面气象参数；
步骤四、气溶胶参数、NO2遥感参数缺失反演及升尺度；
步骤五、交通变量、土地利用变量、社会经济及POI变量、时空变异变量的提取；
步骤六、将各类时空大数据的协变量数据进行时空融合，形成统一尺度及空间坐标的数据集；
步骤七、反演空气污染浓度地表参数；
步骤八、精度验证及评估；
步骤九、验证达标与否；如若达标，进入步骤十一；如若不达标，进入步骤十；
步骤十、对于不合理或不达标预测，调整预测协变量、超参数及限制性条件循环训练，直至得到合理的模型及预测；
步骤十一、结果输出；
在步骤三中，采用一体化的插值方法进行高时空分辨率气象数据的插补，汇总气象数据中的特征协变量，建立目标变量统一且引入注意力层的全残差深度网络模型；对各特征协变量加权，设计多变量输出层，采用多个目标变量同时输出，便于模型参数间共享；在建模过程中，令特征的注意力权重输入为：c＝{ci}，并由此对输入的特征进行加权：
其中，x＝{xi}为特征值矢量输入，xi则为第i个特征，c＝{ci}为特征值矢量的注意力权重，ci为第i个特征的注意力权重系数；
为被注意力权重加权之后的特征输出，
代表矩阵元素相乘；F为特征的个数；
采用softmax层来实现注意力层，以满足式2的限制条件：
其中，ei为注意力层的第i单元的输入，ci为第i个单元的权重输出，exp(…)为取指数函数；j为特征的索引，ej为注意力层的j索引的输入，F为特征的个数；注意力的权重系数通过模型优化自适应求解；
多变量输出的损失函数为：
其中，N为样本数，
代表空气温度的损失，
代表相对湿度的损失，
代表空气压强的损失，
代表风速的损失；yt、yr、yp及yw分别代表空气温度、相对湿度、压强及风速的正则化后的观察值，而
及
分别代表空气温度、相对湿度、压强及风速的估计值；θW,b为网络模型的需要优化的参数集合，Ω(θW,b)则为参数集θW,b的正则化项，采用弹性网络进行正则化；
最后，对目标变量结果进行反正则化，获得同原数据尺度一致的估计值；
在步骤四中，气溶胶参数指的是MAIAC AOD，NO2遥感参数指的是OMI-NO2，在步骤三生成的高分气象参数基础上，分别建立全残差深度网络的AOD及OMI-NO2的缺值反演模型；在缺值反演模型中引入注意力层，提升重要协变量权重，以输入协变量+缺值需要插补的变量共同作为输出的目标变量；
再分别对MAIAC AOD及OMI-NO2建立缺值反演及升尺度模型，采用MSE作为目标损失函数：
其中，N为样本数，lMSE表示采用均方误差的损失函数，y代表了地面空气污染监测站点的观察的NO2浓度，
为地面高分辨率的尺度转换后得到的地面的OMI-NO2的代理变量，引入了坡度因子s及截距a，通过其转换成地面NO2，使其同地面之间的观察值误差最低，从而优化参数得到高分辨率的尺度化结果，对于缺失值直接采用尺度转化模型输入协变量即可得到缺失值的估计；
在步骤五中，交通变量的提取采用最近邻提取或缓冲区分析法；缓冲区分析法计算缓冲区内主干道的长度，通过敏感性分析，选取最优缓冲距离作为最后距离：
其中，d*代表最优缓冲距离，r为相关系数函数，rlen(di)代表以di为缓冲距离时统计缓冲区内的道路长度，x则为目标点空气污染物的浓度，i为测试的缓冲距离索引；
土地利用变量的提取通过计算一定缓冲区范围内土地利用所占的面积比例，将此比例作为污染来源的代理变量；作敏感性分析，确定最优的缓冲距离：
其中，d*代表最优缓冲距离，r为相关系数函数，lprop(di)代表以di为缓冲距离时统计缓冲区内的土地利用面积比例，x则为目标点空气污染物的浓度，i为测试的缓冲距离索引；
社会经济指人口密度及GDP，通过空间叠加分析获得样本相应的人口密度及GDP；POI数据则先选择污染源相关的兴趣点，计算到最近的POI的距离；
时空变异变量的提取包括空间变量、时间变量的提取两方面，其中，空间变量包括高程、坐标及其派生变量，用于捕捉空间变异信息；时间变量包括多尺度时间信息，包括年-天、月份、星期、及年份；
在步骤七中，基模型采用引入注意力层的全残差深度网络模型，输入变量包括步骤一至步骤五获得的参数数据，设输入变量个数为D个，注意力节点也为D个，对关键变量进行加权，输出变量为M个，即需要估计的变量个数；同时，对输出设定了限制性条件，通过限制性优化方法，取得符合限制性条件的解；输出的损失函数采用下式：
其中，L(θW,b)代表了网络参数集θW,b的损失函数，N为样本数，
及
分别代表NO2及PM2.5正则化后的观察值，而
及
分别代表了针对输入协变量集x的NO2及PM2.5的网络预测值，
及
分别代表NO2及PM2.5的MSE的损失函数，Ω(θW,b)则为参数集θW,b的正则化项，采用弹性网络实现正则化；式9及式10定义了NO2及PM2.5的网络预测值需要满足的条件，其中，
及
分别表示NO2的取值的最小值与最大值，
及
分别表示PM2.5的取值的最小值与最大值；
将二污染物合在一起在一个网络中输出，采用bootstrap重抽样方法，对样本及特征进行多次重复抽样，得到不同的样本及特征变量集，对这些不同的样本集合对全残差深度模型结构进行样本数随机化，以获得差异较大的基模型，将这些样本分别训练这些基模型，最后将训练的模型分别预测，预测结果反正则化得到NO2与PM2.5的估计浓度，统计多个模型预测浓度得到估计平均值及其变差。
2.根据权利要求1所述的多源时空大数据深度融合的空气污染预测方法，其特征在于：步骤二中，数据预处理包括对数据进行质量控制、时空融合及格式转换；质量控制是按照有效数据阈值，或者质量标记对数据进行清理，删除无效数据；时空融合是对多源大数据在时空上进行融合，取得统一的目标时空分辨率；格式转换是将格式不一致的数据，转换成统一的高分空间栅格数据。
3.根据权利要求1所述的多源时空大数据深度融合的空气污染预测方法，其特征在于：引入全残差深度网络模型为基模型的聚集引导聚集建模，训练多个基模型，计算各参数的均值及变差，获得更优的地面气象参数反演结果。
4.根据权利要求1所述的多源时空大数据深度融合的空气污染预测方法，其特征在于：步骤八中，采用独立性验证方法对聚集引导得到的结果进行精度验证及结果评估。
5.根据权利要求1所述的多源时空大数据深度融合的空气污染预测方法，其特征在于：步骤九中，验证总体的预测目标是否符合要求，验证各个协变量对总的预测解释是否合理。
6.根据权利要求1所述的多源时空大数据深度融合的空气污染预测方法，其特征在于：步骤十一中，对于得到的合理的优化的训练模型及超参数，保存相应的模型及参数，并将合理的预测结果输出，供下一步的应用使用。
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 CN112905560A (zh) | 2021-06-04 |
| CN112905560B true 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) |
Families Citing this family (22)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113297528B (zh) * | 2021-06-10 | 2022-07-01 | 四川大学 | 一种基于多源大数据的no2高分辨率时空分布计算方法 |
| CN113535510B (zh) * | 2021-06-24 | 2024-01-26 | 北京理工大学 | 一种大规模数据中心数据采集的自适应抽样模型优化方法 |
| CN113434958B (zh) * | 2021-06-25 | 2023-09-01 | 南京工程学院 | 基于异构集成代理模型的磁悬浮飞轮电机多目标优化方法 |
| CN113408918B (zh) * | 2021-06-28 | 2022-04-19 | 哈尔滨工业大学 | 基于多时相遥感分析的村镇垃圾降尺度时空分布反演方法 |
| CN114139719B (zh) * | 2021-11-16 | 2025-01-28 | 中国科学院空天信息创新研究院 | 一种基于机器学习的多源人为热时空量化方法 |
| CN114676121B (zh) * | 2022-03-30 | 2024-06-28 | 重庆大学 | 基于多尺度时空记忆共享网络的缺失监测数据填补方法 |
| CN114609008B (zh) * | 2022-03-30 | 2023-09-01 | 中国科学院地理科学与资源研究所 | 一种气体浓度预测模型的训练方法和气体浓度预测方法 |
| CN114782871B (zh) * | 2022-04-29 | 2022-11-25 | 广东技术师范大学 | 一种基于物联网的海洋异常信息监测方法和装置 |
| CN114896826B (zh) * | 2022-07-13 | 2022-09-09 | 成都信息工程大学 | 基于物理学结合残差注意力网络的行星边界层参数化方法 |
| CN114926749B (zh) * | 2022-07-22 | 2022-11-04 | 山东大学 | 基于遥感图像的近地面大气污染物反演方法及系统 |
| CN115345075B (zh) * | 2022-08-17 | 2023-04-18 | 北京城市气象研究院 | 一体化气溶胶污染气象指数-气溶胶浓度估算方法及系统 |
| CN116596158A (zh) * | 2023-06-14 | 2023-08-15 | 深圳市汉宇环境科技有限公司 | 区域污染源排放总量预测方法 |
| CN117113012A (zh) * | 2023-06-15 | 2023-11-24 | 中国水利水电第十二工程局有限公司 | 一种基于等级控制的平差数据预处理系统 |
| CN116910457B (zh) * | 2023-09-13 | 2023-12-08 | 北京师范大学 | 一种基于区域的污染物反演方法和装置 |
| CN117216480B (zh) * | 2023-09-18 | 2024-06-28 | 宁波大学 | 一种深度耦合地理时空信息的近地表臭氧遥感估算方法 |
| CN117129638B (zh) * | 2023-10-26 | 2024-01-12 | 江西怡杉环保股份有限公司 | 一种区域空气环境质量监测方法及系统 |
| CN117390008B (zh) * | 2023-12-11 | 2024-04-12 | 北京星球空天信息技术有限公司 | 多类型观测仪器的测量数据处理方法和装置 |
| CN117592769B (zh) * | 2024-01-19 | 2024-04-05 | 四川绿豆芽信息技术有限公司 | 一种碳小屋站点管理方法及系统 |
| CN117592005B (zh) * | 2024-01-19 | 2024-04-26 | 中国科学院空天信息创新研究院 | Pm2.5浓度卫星遥感估算方法、装置、设备及介质 |
| CN118196151B (zh) * | 2024-03-13 | 2024-09-24 | 中国人民解放军国防科技大学 | 一种基于气象敏感性分析的光谱反演参数预测方法及装置 |
| CN118094112B (zh) * | 2024-04-18 | 2024-07-30 | 北京市农林科学院智能装备技术研究中心 | 灌溉数据同化方法、装置、电子设备及存储介质 |
| CN118380078A (zh) * | 2024-06-26 | 2024-07-23 | 山东建筑大学 | 一种基于多任务学习的pm2.5和o3污染贡献定量估算方法 |
Family Cites Families (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106909722B (zh) * | 2017-02-10 | 2019-07-26 | 广西壮族自治区气象减灾研究所 | 一种近地面气温的大面积精准反演方法 |
| 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深度学习反演方法 |
-
2021
- 2021-02-02 CN CN202110144010.2A patent/CN112905560B/zh active Active
Also Published As
| Publication number | Publication date |
|---|---|
| CN112905560A (zh) | 2021-06-04 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN112905560B (zh) | 一种多源时空大数据深度融合的空气污染预测方法 | |
| Liao et al. | Deep learning for air quality forecasts: a review | |
| Seng et al. | Spatiotemporal prediction of air quality based on LSTM neural network | |
| CN110738252B (zh) | 空间自相关的机器学习卫星降水数据降尺度方法、系统 | |
| Yi et al. | Inferencing hourly traffic volume using data-driven machine learning and graph theory | |
| CN111260111B (zh) | 基于气象大数据的径流预报改进方法 | |
| Yang et al. | Analyzing historical land use changes using a Historical Land Use Reconstruction Model: A case study in Zhenlai County, northeastern China | |
| CN114936691B (zh) | 一种融合关联度加权与时空注意力的温度预报方法 | |
| CN112699959B (zh) | 基于能量泛函模型的多源多尺度降水数据融合方法和装置 | |
| Ma et al. | Statistical spatial-temporal modeling of ambient ozone exposure for environmental epidemiology studies: A review | |
| Yu et al. | Deep learning-based downscaling of tropospheric nitrogen dioxide using ground-level and satellite observations | |
| CN114254802B (zh) | 气候变化驱动下植被覆盖时空变化的预测方法 | |
| Wu et al. | Seasonal prediction of daily PM2. 5 concentrations with interpretable machine learning: a case study of Beijing, China | |
| Shi et al. | A balanced social LSTM for PM2. 5 concentration prediction based on local spatiotemporal correlation | |
| He et al. | Comparative performance of the LUR, ANN, and BME techniques in the multiscale spatiotemporal mapping of PM 2.5 concentrations in north China | |
| Tiancheng et al. | Application of Improved Naive Bayesian‐CNN Classification Algorithm in Sandstorm Prediction in Inner Mongolia | |
| CN108764527B (zh) | 一种土壤有机碳库时空动态预测最优环境变量筛选方法 | |
| CN117370813A (zh) | 一种基于k线模式匹配算法的大气污染深度学习预测方法 | |
| Li et al. | Deepphysinet: Bridging deep learning and atmospheric physics for accurate and continuous weather modeling | |
| Li et al. | Improving air quality assessment using physics-inspired deep graph learning | |
| Alsaedi et al. | Spatial and temporal data analysis with deep learning for air quality prediction | |
| Li et al. | DeepSAT4D: Deep learning empowers four-dimensional atmospheric chemical concentration and emission retrieval from satellite | |
| CN116223395A (zh) | 近地表痕量气体浓度反演模型及反演方法 | |
| Yang et al. | A novel PM2. 5 concentrations probability density prediction model combines the least absolute shrinkage and selection operator with quantile regression | |
| CN113901348A (zh) | 一种基于数学模型的钉螺分布影响因素识别与预测方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |