一种基于MODIS和长短时记忆网络模型的PM2.5反演方法与流程

一种基于modis和长短时记忆网络模型的pm2.5反演方法
技术领域
1.本发明涉及遥感影像处理技术领域，具体涉及一种基于modis和长短时记忆网络模型的pm2.5反演方法。本方法可用于反演整个研究区域的pm2.5浓度，对环境污染治理和pm2.5监测方面具有一定的参考意义。


背景技术：

2.细颗粒物是指空气动力学中直径小于或等于2.5微米的大气悬浮颗粒物，主要来自自然过程以及大量人类活动的排放。颗粒物污染与人们的健康息息相关，人类通过呼吸将细颗粒物带入肺泡，容易造成心血管疾病及呼吸性疾病。近年来，中国经济飞速发展的同时给环境带来了巨大的压力，各种工业尾气和机动车尾气的排放引发了严重的空气污染问题，pm2.5污染一直是主要的环境问题，并且已经成为公众关心的主要问题。由于pm2.5监测站数量有限，且分布不均匀，难以提供长期连续、覆盖面较广的pm2.5数据，于是利用pm2.5和大气气溶胶之间的关系，使用气溶胶卫星遥感反演产品来间接估算pm2.5浓度已经成为研究热点。气溶胶光学厚度(aod)是指整个大气层中介质消光系数在垂直方向上的积分，用于描述大气中颗粒物对太阳辐射的削弱程度，反映了该大气层中颗粒物含量或空气污染程度。利用aod和pm2.5之间的关系可以进行全天候、大范围的连续观测，此外，遥感技术不断向高时间分辨率、高空间分辨率方向发展，气溶胶反演产品的精度也随着反演算法的改进得到提高，使得卫星遥感技术逐渐成为反演pm2.5浓度的有效手段。
3.aod数据目前已经被广泛用于pm2.5浓度的估计，但由于云雨天气、云层覆盖、极高的气溶胶负荷等原因引起aod数据经常存在大范围缺失问题，人们通常在进行pm2.5浓度估计时会直接忽略这些缺失，但这样会导致pm2.5浓度估计的结果不够准确，因此填补aod数据的缺失可以大大提高pm2.5浓度估计的准确率。为此提出了多种方法来处理缺失的aod数据。一种方法是融合多种卫星数据来填补aod数据，例如将terra卫星与aqua卫星的数据结合起来，利用其平均值、融合值、季节平均值、多重插补方法来减少aod数据的缺失，或者可以将aod数据与aeronet数据合并来提高aod覆盖率；另外也可以结合使用插值方法如克里金插值，反距离加权法，然而这些方法不能在减少缺失率的同时提高反演精度，在大范围缺失的情况下通常是不够准确的。近年来，人们利用一些先进的统计模型来表示pm2.5和各种预测因子之间的关系，如地理加权回归模型，两阶段模型。然而，这些模型可能不能够完全描述pm2.5和预测因子之间的关系。
4.深圳先进技术研究院申请的专利“一种基于modis的pm2.5遥感反演方法”(专利申请号201711398781.4)中公开了一种基于modis的pm2.5遥感反演方法，该方法将pm2.5数据插值成pm2.5影像，并进行训练集和测试集的划分，将训练集用于机器学习算法的训练，最后选取一个最优模型并用于整幅modis影像，最后得到整幅modis的pm2.5浓度反演图。该方法不足之处在于没有考虑modis影像的缺失情况，直接建立其与pm2.5之间的关系，会产生一定的误差，造成实验结果不够准确。
5.针对现有pm2.5反演技术的不足，目前需要一种可以充分利用所获取的modis数据
和pm2.5数据的反演方法，本发明针对modis aod数据进行机器学习填补，随后再进行modis与pm2.5浓度的关系建立，以便于提高pm2.5反演精度。


技术实现要素：

6.为了克服上述现有技术的不足，本发明提出一种基于modis和长短时记忆网络模型的pm2.5反演方法。本发明与现有技术相比通过使用机器学习方法对modis aod数据采集中的缺失值进行填补，增加了数据的完整性，同时对辅助数据进行高湿订正，避免一些气象条件的干扰，采用长短时记忆网络模型提高了pm2.5反演精度。
7.本发明解决上述问题的技术方案如下：一种基于modis和长短时记忆网络模型的pm2.5反演方法，包括以下步骤：
8.步骤1：获取需要反演的pm2.5研究时段的modis影像、用于计算aod的气象数据和研究区域pm2.5监测站数据；气象数据包括风速、温度、边界层高度、相对湿度和地表压力；对各时段的modis影像进行融合并裁剪出研究区域对应的部分；对气象数据进行重采样及归一化处理；对监测站数据求取9-15点均值作为当天pm2.5值；
9.步骤2：对气象数据按如下式(1)(2)进行高湿订正，消除高度和湿度的影响；
10.aod＝k
a,0
(λ)blh(1)
11.式中，k
a,0
(λ)表示近地层气溶胶系数，blh表示边界层高度；
12.f(rh)＝1/(1-rh/100)(2)
13.式中，rh代表相对湿度；
14.步骤3：使用双线性插值法和重采样技术将所有的气象数据变量均重采样到0.01
×
0.01分辨率；
15.步骤4：训练随机森林模型；随机森林模型学习到带有aod值的网格与对应位置的气象数据之间的函数关系，利用此函数关系估计缺失位置的aod值，得到覆盖整个区域的连续、大范围的aod值；最后在模型拟合之后计算每个气象数据变量的特征重要性；
16.在随机森林模型中，使用bootstrapsample方法从步骤3得到的气象数据中提取得到样本z，然后为每个样本生成一棵带有多个预测因子的未回归树，在未回归树的每个节点选择最优值进行最佳分割得到回归树，最后通过所有回归树的平均值来进行预测；每棵未回归树三分之一的训练样本不参与回归树f(x)的生成；
[0017][0018][0019]
z1(m,n)＝{x|xj≤n}and z2(m,n)＝{x|xj＞n}(5)
[0020][0021]
[0022]
其中(xi,yi)是z区域(r1,r2,...,rz)中i＝1,2,...,n的样本，cm是对模型的响应，是最优值，m是分离变量，n是分离点，i是示性函数；
[0023]
步骤5：从步骤3得到的气象数据和步骤4得到的覆盖整个区域的aod数据中提取pm2.5监测站点坐标对应的气象数据值以及aod值，按照时间和空间将提取出来的aod值和气象数据值与pm2.5站点值进行匹配，匹配后的数据按7:3的比例分割成训练集和测试集；构建第f个时间节点的特征向量基数tf：
[0024]
tf＝[x1,x2,...,xn,...,xn] (8)
[0025]
其中，x1表示t个时间节点中第f个时间节点的pm2.5的值，xn表示t个时间节点中第f个时间节点的第n个空气参数的值；
[0026]
步骤6：建立长短时记忆网络模型；长短时记忆网络模型包括lstm层和dense层，初始化长短时记忆网络模型的参数并将覆盖整个研究区域的aod数据及完整的气象变量输入，得到输出预测值w，将输出预测值w与监测站点真实值y进行对比，根据误差微调网络参数，实现降低误差；遗忘门、输入门、循环门和输出门为长短时记忆网络模型的主要结构：
[0027]
a.遗忘门：将信息进行选择性遗忘，减少模型的记忆负担；
[0028]ft
＝σ(wf·
[h
t-1
,x
t
] bf) (9)
[0029]
b.输入门：将信息进行选择性记忆：
[0030]it
＝σ(wi·
[h
t-1
,x
t
] bi) (10)
[0031]
c.循环门：更新当前时刻的状态：
[0032][0033]
d.输出门：完成模型输出；
[0034]ot
＝σ(wo·
[h
t-1
,x
t
] bo) (12)
[0035]ht
＝o
t
*tanh(c
t
) (13)
[0036]
式中，σ为sigmoid函数，wo和bo作为待定系数在后续进行训练学习，h
t-1
是上一阶段的输出，x
t
是当前阶段的输入，式中的c
t-1
记录了上一时间节点的信息；
[0037]
步骤7：将得到的输出预测值进行可视化，得到整个研究区域的pm2.5浓度反演图。
[0038]
与现有技术相比，本发明的有益效果：本发明着眼于modis和长短时记忆网络模型进行pm2.5浓度反演。第一，利用多种辅助数据，建立aod和辅助数据之间的关系，利用机器学习方法对aod缺失值部分进行填补，减少了缺失值带来的反演不确定性，增加了反演精度。第二，对数据进行高度订正和湿度订正，减少了高度和湿度带来的误差。第三，使用长短时记忆网络进行pm2.5浓度反演，考虑了时间和空间的依赖关系，大大提高了反演精度。
附图说明
[0039]
图1是基于modis和长短时记忆网络模型的pm2.5浓度反演流程图。
[0040]
图2是研究区域内pm2.5浓度估计值和观察值的每日pm2.5浓度散点图。
[0041]
图3是研究区域内pm2.5浓度反演图。
具体实施方式
[0042]
以下结合技术方案和附图详细叙述本发明的具体实施方式。
[0043]
参照附图1，一种基于modis和长短时记忆网络模型的pm2.5反演方法，包括以下步骤：
[0044]
步骤1：获取需要反演的pm2.5研究时段的modis影像、用于计算aod的气象数据和研究区域pm2.5监测站数据；气象数据包括风速、温度、边界层高度、相对湿度和地表压力；对各时段的modis影像进行融合并裁剪出研究区域对应的部分；对监测站数据求取9-15点均值作为当天pm2.5值；
[0045]
步骤2：对气象数据进行重采样及归一化处理和高湿订正处理；按如下式(1)(2)进行高湿订正，消除高度和湿度的影响；
[0046]
aod＝k
a,0
(λ)blh (1)
[0047]
式中，k
a,0
(λ)表示近地层气溶胶系数，blh表示边界层高度；
[0048]
f(rh)＝1/(1-rh/100) (2)
[0049]
式中，rh代表相对湿度；
[0050]
步骤3：使用双线性插值法和重采样技术将所有的气象数据变量均重采样到0.01
×
0.01分辨率；
[0051]
步骤4：训练随机森林模型；随机森林模型学习到带有aod值的网格与对应位置的气象数据之间的函数关系，利用此函数关系估计缺失位置的aod值，得到覆盖整个区域的连续、大范围的aod值；最后在模型拟合之后计算每个气象数据变量的特征重要性；
[0052]
在随机森林模型中，使用bootstrapsample方法从步骤3得到的气象数据中提取得到样本z，然后为每个样本生成一棵带有多个预测因子的未回归树，在树的每个节点选择最优值进行最佳分割得到回归树，最后通过所有回归树的平均值来进行预测；每棵决策树三分之一的训练样本不参与第k个树和回归树f(x)的生成；
[0053][0054][0055]
z1(m,n)＝{x|xj≤n}and z2(m,n)＝{x|xj＞n} (5)
[0056][0057][0058]
其中(xi,yi)是z区域(r1,r2,...,rz)中i＝1,2,...,n的样本，cm是对模型的响应，是最优值，m是分离变量，n是分离点，i是示性函数；
[0059]
步骤5：从步骤3得到的气象数据和步骤4得到的覆盖整个区域的aod数据中提取pm2.5监测站点坐标对应的气象数据值以及aod值，按照时间和空间将提取出来的aod值和气象数据值与pm2.5站点值进行匹配，匹配后的数据按7:3的比例分割成训练集和测试集；构建第f个时间节点的特征向量基数tf：
[0060]
tf＝[x1,x2,...,xn,...,xn] (8)
[0061]
其中，x1表示t个时间节点中第f个时间节点的pm2.5的值，xn表示t个时间节点中第f个时间节点的第n个空气参数的值；
[0062]
步骤6：建立长短时记忆网络模型；长短时记忆网络模型包括lstm层和dense层，初始化长短时记忆网络模型的参数并将覆盖整个研究区域的aod数据及完整的气象变量输入，得到输出预测值w，将输出预测值w与监测站点真实值y进行对比，根据误差微调网络参数，实现降低误差；遗忘门、输入门、循环门和输出门为长短时记忆网络模型的主要结构：
[0063]
a.遗忘门：将信息进行选择性遗忘，减少模型的记忆负担；
[0064]ft
＝σ(wf·
[h
t-1
,x
t
] bf) (9)
[0065]
b.输入门：将信息进行选择性记忆：
[0066]it
＝σ(wi·
[h
t-1
,x
t
] bi) (10)
[0067]
c.循环门：更新当前时刻的状态：
[0068][0069]
d.输出门：完成模型输出；
[0070]ot
＝σ(wo·
[h
t-1
,x
t
] bo) (12)
[0071]ht
＝o
t
*tanh(c
t
) (13)
[0072]
式中，σ为sigmoid函数，wo和bo作为待定系数在后续进行训练学习，h
t-1
是上一阶段的输出，x
t
是当前阶段的输入，式中的c
t-1
记录了上一时间节点的信息；
[0073]
步骤7：
[0074]
将得到的输出预测值进行可视化，得到整个研究区域的pm2.5浓度反演图。
[0075]
下面结合附图2和附图3对本发明的效果进一步说明。
[0076]
仿真实验：
[0077]
本发明仿真实验的硬件环境为：intel core i7处理器，2.5ghz主频，8g内存；软件环境为：windows10操作系统，matlab r2016b。
[0078]
本发明研究区域为北京市全域，图像获取时间为2018年1月1日到2018年12月31日，分辨率为1km。
[0079]
通过以下指标定量评估本发明和对比实验的反演结果：建模过程中使用的评估指标包括r2和rmse。
[0080]
rmse均方根误差：
[0081][0082]
r2：
[0083][0084]
n为数据个数；为真实值，yi为预测值，为预测数据的均值
[0085]
图2是研究区域内pm2.5浓度估计值和观察值的每日pm2.5浓度散点图，实线是根
据估计值与观察值的拟合回归线，虚线是1：1参考线。
[0086]
图3是研究区域内pm2.5浓度反演图。
[0087]
仿真结果分析：
[0088]
表1 pm2.5反演精度表
[0089][0090]
结果表明，本方法的结果较之多元线性回归模型(mlr)和混合效应模型(lme)效果有一定程度的提升。aod插补对于pm2.5浓度估计结果是有效的。利用该长短时记忆网络模型进行计算，得到北京pm2.5浓度空间分布模拟图如图3所示。北京pm2.5浓度总体呈自北向南增加的趋势。北京西部和北部山区pm2.5浓度明显偏低。这主要是因为市区主要集中在东南部，受工业生产和人类活动的影响。空气污染较为严重，pm2.5浓度偏低主要发生在北方山区和农村地区。
[0091]
以上所述实施例仅表达本发明的实施方式，但并不能因此而理解为对本发明专利的范围的限制。应当指出，对于本领域的技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替换，都应当视为属于本发明的保护范围。