CN113901384A - 顾及全局空间自相关性和局部异质性的地面pm2.5浓度建模方法 - Google Patents

顾及全局空间自相关性和局部异质性的地面pm2.5浓度建模方法 Download PDF

Info

Publication number
CN113901384A
CN113901384A CN202111121091.0A CN202111121091A CN113901384A CN 113901384 A CN113901384 A CN 113901384A CN 202111121091 A CN202111121091 A CN 202111121091A CN 113901384 A CN113901384 A CN 113901384A
Authority
CN
China
Prior art keywords
spatial
concentration
local
global
space
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.)
Pending
Application number
CN202111121091.0A
Other languages
English (en)
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.)
Wuhan University WHU
Original Assignee
Wuhan University WHU
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 Wuhan University WHU filed Critical Wuhan University WHU
Priority to CN202111121091.0A priority Critical patent/CN113901384A/zh
Publication of CN113901384A publication Critical patent/CN113901384A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N15/00—Investigating characteristics of particles; Investigating permeability, pore-volume, or surface-area of porous materials
    • G01N15/06—Investigating concentration of particle suspensions

Abstract

本发明提出了一种顾及全局空间自相关性和局部异质性的地面PM2.5浓度建模方法。本发明针对PM2.5监测站点数量有限且分布不均，且传统PM2.5回归建模方法缺少对局部空间影响和全局空间影响同时考虑的问题，提出了一种顾及全局空间自相关性和局部异质性的地面PM2.5浓度建模方法，同时考虑全局和局部空间效应，消除空间自相关因素和空间异质因素的影响。分别通过构建全局空间权重矩阵和局部空间权重矩阵，提取出全局空间影响因子和局部空间影响因子。然后再根据两者构建回归模型，得到地面PM2.5浓度估计模型，以供后续研究和分析。

Description

顾及全局空间自相关性和局部异质性的地面PM2.5浓度建模 方法
技术领域
本发明属于空间统计分析服务应用技术领域，特别涉及一种顾及全局空间自相关性和局部空间异质性的地面PM2.5浓度建模方法。
背景技术
PM2.5是空气污染的主要污染物，给人们生产和生活造成巨大影响的同时还对人体健康有着极大的危害。近年来空气污染问题越来越被重视，地面PM2.5浓度监测站在全国范围内建立了起来，但监测站点数量有限且分布不均。随着遥感技术的发展，有学者开始采用遥感数据来反演PM2.5浓度。常用的PM2.5反演数据有气溶胶光学厚度(AOD)、归一化植被指数(NDVI)、数字高程模型(DEM)，以及各种气象因子和污染源数据等。利用遥感数据进行PM2.5反演常用模型有普通最小二乘回归(OLS)、土地利用回归(LUR)、地理加权回归(GWR)、空间滤值(ESF)和空间误差(SEM)等模型。
PM2.5的分布受到空间效应的影响，具有空间自相关性和空间异质性。PM2.5在空间上的分布是连续的，存在空间自相关性；而PM2.5又受到地形、气象因子、污染源等因素的影响，这些因素又随着位置的改变而改变，具有空间异质性。传统PM2.5回归建模方法只单独考虑全局空间自相关性或局部空间异质性，没有同时考虑全局和局部的空间影响，即回归模型只在一定程度上削弱了空间效应的影响。
为了同时消除全局和局部空间效应的影响，本发明提出了一种顾及全局空间自相关性和局部空间异质性的地面PM2.5浓度建模方法。本发明认为局部空间影响可以分为两部分，一部分是局部环境条件本身造成的空间影响；另一部分则是全局空间影响以及该全局影响在局部的表现。这就是全局和局部空间回归建模的结合点。本发明分别通过全局空间权重矩阵和局部空间权重矩阵提取出全局空间影响因子和局部空间影响因子。在提取局部空间影响因子时，把全局空间影响因子作为因子数据，同气溶胶厚度、气象因子、污染源等因子数据一起提取出局部空间影响因子，然后再与全局空间影响因子一起进行回归建模，得到地面PM2.5浓度估算模型。
发明内容
基于遥感数据和污染源数据，对地面站点PM2.5浓度进行建模，并估计无站点位置的PM2.5浓度。发明旨在同时考虑全局和局部的空间效应，结合全局和局部回归思想，消除全局空间自相关因素和局部异质因素的影响，进一步提高PM2.5建模精度和估算效果。
参见图1，本发明所采用的技术方案包括如下步骤：
步骤1：数据获取与处理。获取气象站点的PM2.5地面监测值，并从各数据源获取与地面PM2.5相关的因子数据。气溶胶厚度产品与PM2.5浓度显著相关，为必需自变量。气象数据、污染源数据以及其他因子数据根据研究区气候特征和数据质量选取。对所有数据进行预处理，包括删除空值和异常值、投影转换、时间和空间分辨率的统一。
步骤2：构建全局空间权重矩阵。根据气象站点数据的投影坐标，构建研究区域内所有站点的全局空间权重矩阵，并中心化处理为对称矩阵，将空间关系转化为权重值。构建方法采用：基于邻接关系构建；或者基于反距离加权构建；或者复合型。
步骤3：全局空间因子提取。基于全局空间权重矩阵，根据空间自相关程度大小提取出全局空间影响因子，记为Eα。
步骤4：构建局部空间权重矩阵。选取空间权函数，构建每一个站点与该站点局部空间影响范围内其余站点的局部空间权重矩阵，局部空间影响范围在步骤5中优化确定。
步骤5：确定局部空间影响范围。局部空间影响范围也称为带宽，权函数对带宽的选择十分敏感，需要选择出最优的带宽。常用的最优带宽选择的方法有交叉验证法、AIC准则、贝叶斯信息准则以及平稳指数等。
步骤6：局部空间因子提取。根据局部空间权重矩阵和最优带宽，进行权重确定的局部最小二乘回归，提取出局部空间影响因子，记为βx。
步骤7：构建全局-局部空间回归模型。将步骤3中提取出的全局空间影响因子Eα和步骤6提取的局部空间影响因子βx进行回归建模,具体是将代表全局空间自相关的全局空间影响因子Eα和代表局部空间异质性的局部空间影响因子βx，利用全局空间因子和局部空间因子构建回归模型，构成如下模型：
Figure BDA0003277183590000031
其中，E是全局空间影响因子空间特征向量的集合，α为全局空间影响因子空间特征向量集的系数；(ui，vi)为第i个采样点的坐标系的横纵坐标，β0(ui，vi)为第i个采样点的回归方程的截距项；βk(ui，vi)为第i个采样点上第k个回归参数，xik是第i个采样点上第k个自变；αi是第i个采样点的空间特征向量集合的系数；εi为第i个区域的随机误差，满足零均值、同方差、相互独立等基本假定。
步骤8：评价与分析。计算模型精度评价指标。评价指标包括模型的拟合优度(R2)、调整后拟合优度(Adj.R2)、均方根误差(RMSE)、平均绝对百分比误差(MAPE)、平均绝对误差(MAE)以及残差的莫兰指数(Moran’s I)。并对模型进行鲁棒性检验，常采用交叉验证的方法计算上述指标。
在上述的顾及全局空间自相关性和局部异质性的地面PM2.5浓度建模方法，步骤1中，统一因子的时间和空间分辨率，并进行投影转换，与AOD数据的时空分辨率和空间参考系一致。其中的污染源数据，如工厂位置、道路网密度等点、线数据，计算它们的核密度从而得到覆盖整个研究区域的栅格影像图。
在上述的顾及全局空间自相关性和局部异质性的地面PM2.5浓度建模方法，步骤2具体包括
步骤2.1：
(1)基于邻接关系构建空间权重矩阵W1。构建站点的泰森多边形，可以根据Bishop邻接、Rook邻接和Queen邻接方式来构建全局空间权重矩阵。
(2)基于反距离加权构建空间权重矩阵W1。根据站点间的距离构建邻接关系，通过核函数将距离转为权重值，得到全局空间权重矩阵。核函数可分为指数、高斯、球状模型具体公式如下：
Figure BDA0003277183590000032
Figure BDA0003277183590000033
Figure BDA0003277183590000034
其中，公式(1)、(2)、(3)分别为指数、高斯、球状模型。i，j分别表示位置点i和位置点j。Wi，j表示位置点i和j之间的邻接性(权重)。r表示所有站点的最小生成树中边的最大距离。
(3)基于复合型构建空间权重矩阵W1。如K-最近邻矩阵和基于邻接关系和距离构建权重矩阵。K-最近邻矩阵，首先确定一个阈值K，搜寻每一个站点最近邻的K个站点，构建邻接矩阵，Wi，j表示位置点i和j之间的邻接性(权重)；基于邻接关系和距离构建权重矩阵，首先确定一个阈值d，搜寻每一个站点距离不超过阈值d的其他站点，认为两者邻接，Wi，j＝1，否则就不邻接，Wi，j＝0，构建权重矩阵。
步骤2.2：矩阵中心化。对上述步骤得到的矩阵W1进行中心化处理为对称矩阵W。中心化公式如下：
Figure BDA0003277183590000041
其中，I表示n维单位矩阵，11T是一个所有元素均为1的n×n矩阵，n表示研究区域内监测站点的数量。
在上述的顾及全局空间自相关性和局部异质性的地面PM2.5浓度建模方法，步骤3具体包括
步骤3.1：计算空间权重矩阵W的特征值和特征向量，并进行初筛。求解特征值，得到特征值λi和特征向量Ei。初筛的条件有两个：一是其特征值λi＞0；二是根据经验模型要求满足如下条件：
Figure BDA0003277183590000042
其中λi表示待筛选的特征值，λmax表示最大的特征值。初步筛选出符合条件的空间特征向量。
步骤3.2：空间特征向量筛选。根据逐步回归法筛选出使评价结果最优的空间特征向量。具体步骤如下：
步骤3.2.1：空间特征向量排序。把初筛后的空间特征向量按照特征值大小从小到大排列，得到空间特征向量集E。
步骤3.2.2：把空间特征向量集E的向量依次加入PM2.5与各个因子的普通最小二乘回归方程作为自变量，计算回归方程的AIC值。
步骤3.2.3：从步骤3.2.3中选择AIC值最小的回归方程Yk，对应的空间特征向量为Ek，然后从空间特征向量集E中将除Ek外的特征向量依次加入回归方程Yk，并计算回归方程的AIC值。
步骤3.2.4：判断AIC值是否明显减小，若AIC值明显减小，则重复步骤3.2.3；若AIC值几乎没减小则停止将特征向量加入回归方程，并记录已经加入到回归方程的特征向量。这些特征向量即为筛选出的使评价标准最优的特征向量集。
在上述的顾及全局空间自相关性和局部异质性的地面PM2.5浓度建模方法，步骤4中，
空间权函数有距离阈值函数、指数函数、高斯函数、双重平方函数、三次立方函数，分别如公式6～11所示。
Figure BDA0003277183590000051
Figure BDA0003277183590000052
Figure BDA0003277183590000053
Figure BDA0003277183590000054
Figure BDA0003277183590000055
其中，Wi，j表示位置点i和j之间的邻接性(权重)，dij表示位置点i和j之间的距离，其中b是用于描述权重与距离之间函数关系的非负衰减参数，称为带宽(bandwidth)。带宽越大，权重随距离增加衰减得越慢；带宽越小，权重随距离增加衰减得越快。
在上述的顾及全局空间自相关性和局部异质性的地面PM2.5浓度建模方法，步骤5中，采用交叉验证(cross validation，CV)法计算出最优带宽。
Figure BDA0003277183590000056
Figure BDA0003277183590000061
表示回归参数估计不包括回归点本身，即只根据回归点周围点进行回归计算。CV值最小的b则为最优带宽。
在上述的顾及全局空间自相关性和局部异质性的地面PM2.5浓度建模方法，步骤8中交叉验证有两种方法。(1)K-折交叉验证。K值根据数据量来确定，K通常选择5、8或10。(2)子抽样验证。将数据分为两份，一份作为建模集，一份作为测试集。
在上述的顾及全局空间自相关性和局部异质性的地面PM2.5浓度建模方法，步骤8具体包括：
步骤8.1：模型精度评价。计算所得模型的R2、调整后R2即Adj.R2、均方根误差(RMSE)、平均绝对误差(MAE)以及残差的莫兰指数(Moran’s I)等作为评价指标，以验证所提出的顾及全局空间自相关性和局部空间异质性的地面PM2.5浓度估算模型的精度。
Figure BDA0003277183590000062
其中yi是站点i的PM2.5浓度观测值，
Figure BDA0003277183590000063
是观测数据的平均值，
Figure BDA0003277183590000064
是模型预测的站点i的PM2.5浓度，n是监控站点的个数。
Figure BDA0003277183590000065
其中p是自变量的个数；R2和Adj.R2的取值范围是0～1，值越大说明模型精度越高。
Figure BDA0003277183590000066
Figure BDA0003277183590000067
式中参数含义同上，RMSE和MAE越小说明模型精度越高
Figure BDA0003277183590000068
其中ei是由模型得到的站点i的PM2.5浓度残差，
Figure BDA0003277183590000069
是平均值，cij是站点i和j之间的反距离空间权重。I的取值范围是-1～1，值越接近于0，残差空间自相关性越弱，模型越可靠。
步骤8.2：交叉验证。交叉验证采用的方法为8折交叉验证法，评估模型对非站点区域的PM2.5浓度的预测精度。将步骤1获得的PM2.5站点数据和自变量数据以及步骤3获得的筛选出的空间特征向量随机均分为8份，每次选7份作为训练集，按前述步骤进行建模，剩余的1份作为测试集，将测试集的自变量带入训练集所得模型中，计算测试集的均方根误差；每份数据都做过一次测试集后，计算8个均方根误差RMSE的均值，即为交叉验证的结果。RMSE越小，模型的预测精度越高，其鲁棒性越强，适用性越强。
步骤8.3：栅格计算。步骤7得到了地面PM2.5浓度的数学统计模型，也得到了整个研究区域各个位置点的各个因子的相关系数。将研究区域的各位置点转换为栅格影像，栅格值分别设为各因子的系数值。将提取的空间特征向量插值为覆盖整个研究区域的具有与AOD相同分辨率的栅格图像，然后进行栅格计算，用因子的相关系数栅格影像乘以因子栅格影像然后相加，最后加上截距项，得到地面PM2.5浓度估算值。
步骤8.4：异常值处理。在步骤10中得到的PM2.5浓度分布图，由于误差的存在，使得极少部分区域的PM2.5浓度值可能为负值，同时也可能存在少量的异常的高PM2.5浓度值，需要对这些异常值进行处理。负的PM2.5浓度值表征该点的PM2.5浓度低，故PM2.5浓度值设为0，同时给PM2.5浓度值设定一个阈值，超过该阈值的PM2.5浓度值均修改为该阈值，进行可视化展示，最终得到连续的地面PM2.5浓度分布图。
因此，本发明具有如下优点：
1、有效消除全局空间自相关因素和局部异质因素的影响。同时考虑了空间自相关性和空间异质性，提取出了全局空间影响因子和局部空间影响因子，提高模型精度。
2、PM建模方法和全局空间影响因子选择十分灵活。全局空间因子的选择有多种，如空间权重矩阵的特征向量，空间滞后模型的滞后项等，局部回归模型的选择也有多种，除地理加权回归模型外还可以选择其他空间变系数回归模型，可以根据研究区域和研究数据质量来选取合适的全局和局部的回归建模方法。
3、PM2.5浓度分布图可以覆盖整个研究区域。根据已有PM2.5监测站点，得到无监测站点处的PM2.5浓度值，从而得到连续的PM2.5空间分布图，便于后续相关分析。
附图说明
图1为本发明技术方案的流程图。
图2为本发明实施例的流程图。
图3为本发明实施例步骤2中构建全局空间权重矩阵的子流程图。
具体实施方法
为了便于本领域普通技术人员理解以及实施本发明，下面结合附图及实施例对本发明作进一步阐述，应当理解，此处所述的实施例仅用于说明和解释本发明，并不用于限定本发明。
本发明要解决的问题是：PM2.5分布受到空间效应影响，传统PM2.5回归建模方法只考虑全局空间自相关性或者局部空间异质性，缺少对局部空间影响和全局空间影响同时的考虑。针对这些问题，本发明提出顾及全局空间自相关性和局部异质性的地面PM2.5浓度建模方法，构建地面PM2.5浓度估算模型，进而制作PM2.5浓度空间分布图。
参见图2，本发明提供的顾及全局空间自相关性和局部异质性的地面PM2.5浓度建模方法，包括如下步骤：
步骤1：数据获取与处理。PM2.5数据来源于全国城市空气质量实时发布平台，平台每小时发布各监测站点的实时PM2.5浓度数据。对PM2.5站点数据进行质量检查，剔除明显异常值。遥感气溶胶光学厚度(AOD)数据可以从现有气溶胶产品中获取。气象因子如气温、气压、相对湿度、行星边界层高度、降水量、风速等，会对PM2.5的生成与扩散过程产生影响，进而影响地面PM2.5浓度空间分布；其他因子如植被覆盖情况、土地利用状况、人口密度、高程、工厂和道路分布密度等，也在一定程度上影响地面PM2.5浓度的空间分布。可根据研究区域的气候特征和数据质量选取因子数据，再对所有数据进行预处理。原始PM2.5数据为每小时平均浓度，可依次计算出日均、月均、季均和年均浓度值。其余影响因子也将时间分辨率统一为与PM2.5数据相同的时间分辨率。在空间分辨率上，其余影响因子需要转化为与AOD数据一致的空间分辨率，可采用重采样或空间插值等方法进行空间分辨率的调整。工厂数据和道路网数据分别为点状数据和线状数据，可以通过计算其核密度转换为覆盖整个研究区域的栅格数据。与气溶胶光学厚度空间参考系不同的数据需要进行投影变换，空间参考系要与气溶胶光学厚度相同。预处理完成后，提取PM2.5监测站点处的各因子数据。
步骤2：构建全局空间权重矩阵。根据气象站点数据的投影坐标可采用三种方法构建空间权重矩阵。(1)基于邻接关系构建；(2)基于反距离加权构建；(3)复合型。参见图2具体步骤如下：
步骤2.1：
(1)基于邻接关系构建空间权重矩阵W1。构建站点的泰森多边形，可以根据Bishop邻接、Rook邻接和Queen邻接方式来构建全局空间权重矩阵。
(2)基于反距离加权构建空间权重矩阵W1。根据站点间的距离构建邻接关系，通过核函数将距离转为权重值，得到全局空间权重矩阵。核函数可分为指数、高斯、球状模型具体公式如下：
Figure BDA0003277183590000091
Figure BDA0003277183590000092
Figure BDA0003277183590000093
其中，公式(1)、(2)、(3)分别为指数、高斯、球状模型。i，j分别表示位置点i和位置点j。Wi，j表示位置点i和j之间的邻接性(权重)。r表示所有站点的最小生成树中边的最大距离。
(3)基于复合型构建空间权重矩阵W1。如K-最近邻矩阵和基于邻接关系和距离构建权重矩阵。K-最近邻矩阵，首先确定一个阈值K，搜寻每一个站点最近邻的K个站点，构建邻接矩阵，Wi，j表示位置点i和j之间的邻接性(权重)；基于邻接关系和距离构建权重矩阵，首先确定一个阈值d，搜寻每一个站点距离不超过阈值d的其他站点，认为两者邻接，Wi，j＝1，否则就不邻接，Wi，j＝0，构建权重矩阵。
步骤2.2：矩阵中心化。对上述步骤得到的矩阵W1进行中心化处理为对称矩阵W。中心化公式如下：
Figure BDA0003277183590000101
其中，I表示n维单位矩阵，11T是一个所有元素均为1的n×n矩阵，n表示研究区域内监测站点的数量。
步骤3：全局空间因子提取。以空间权重矩阵的特征向量来表示所有的全局空间因子，通过特征向量对应的特征值来筛选，并利用逐步回归的方式将显著的特征向量作为全局空间因子提取出来。参见图3，具体步骤如下：
步骤3.1：计算空间权重矩阵W的特征值和特征向量，并进行初筛。求解特征值，得到特征值λi和特征向量Ei。初筛的条件有两个：一是其特征值λi＞0；二是根据经验模型要求满足如下条件：
Figure BDA0003277183590000102
其中λi表示待筛选的特征值，λmax表示最大的特征值。初步筛选出符合条件的空间特征向量。
步骤3.2：空间特征向量筛选。根据逐步回归法筛选出使评价结果最优的空间特征向量。具体步骤如下：
步骤3.2.1：空间特征向量排序。把初筛后的空间特征向量按照特征值大小从小到大排列，得到空间特征向量集E。
步骤3.2.2：把空间特征向量集E的向量依次加入PM2.5与各个因子的普通最小二乘回归方程作为自变量，计算回归方程的AIC值。
步骤3.2.3：从步骤3.2.3中选择AIC值最小的回归方程Yk，对应的空间特征向量为Ek，然后从空间特征向量集E中将除Ek外的特征向量依次加入回归方程Yk，并计算回归方程的AIC值。
步骤3.2.4：判断AIC值是否明显减小，若AIC值明显减小，则重复步骤3.2.3；若AIC值几乎没减小则停止将特征向量加入回归方程，并记录已经加入到回归方程的特征向量。这些特征向量即为筛选出的使评价标准最优的特征向量集。
步骤4：构建局部空间权重矩阵。选取合适的空间权函数，构建每一个站点与该站点局部空间影响范围内其余站点的局部空间权重矩阵。空间权函数有距离阈值函数、指数函数、高斯函数、双重平方函数、三次立方函数，分别如公式6～11所示。
Figure BDA0003277183590000111
Figure BDA0003277183590000112
Figure BDA0003277183590000113
Figure BDA0003277183590000114
Figure BDA0003277183590000115
其中，Wi，j表示位置点i和j之间的邻接性(权重)，dij表示位置点i和j之间的距离，其中b是用于描述权重与距离之间函数关系的非负衰减参数，称为带宽(bandwidth)。带宽越大，权重随距离增加衰减得越慢；带宽越小，权重随距离增加衰减得越快。
步骤5：最优带宽计算。最优带宽计算。带宽的选择十分敏感。带宽过大会导致回归参数估计的偏差过大，而带宽过小又会导致回归参数估计的方差过大。采用交叉验证(cross validation，CV)法计算出最优带宽。
Figure BDA0003277183590000116
Figure BDA0003277183590000117
表示回归参数估计不包括回归点本身，即只根据回归点周围点进行回归计算。CV值最小的b则为最优带宽。
步骤6：回归建模，求取因子系数。构建顾及全局空间自相关性和局部异质性的地面PM2.5浓度估计模型。具体公式如下：
Figure BDA0003277183590000118
其中，(ui，vi)为第i个采样点的坐标；βk(ui，vi)为第i个采样点上第k个回归参数，E是全局空间影响因子空间特征向量的集合，α为全局空间影响因子空间特征向量集的系数，αi是第i个采样点的空间特征向量集合的系数。εi为第i个区域的随机误差，满足零均值、同方差、相互独立等基本假定。
步骤7：PM2.5浓度预测。根据步骤6中建立的模型，估算出整个研究区域无PM2.5监测站点的PM2.5浓度。
步骤8：模型精度评价。计算所得模型的R2、调整后R2(Adi.R2)、均方根误差(RMSE)、平均绝对误差(MAE)以及残差的莫兰指数(Moran’s I)等作为评价指标，以验证所提出的顾及全局空间自相关性和局部空间异质性的地面PM2.5浓度估算模型的精度。
Figure BDA0003277183590000121
其中yi是站点i的PM2.5浓度观测值，
Figure BDA0003277183590000122
是观测数据的平均值，
Figure BDA0003277183590000123
是模型预测的站点i的PM2.5浓度，n是监控站点的个数。
Figure BDA0003277183590000124
其中p是自变量的个数；R2和Adj.R2的取值范围是0～1，值越大说明模型精度越高。
Figure BDA0003277183590000125
Figure BDA0003277183590000126
式中参数含义同上，RMSE和MAE越小说明模型精度越高
Figure BDA0003277183590000127
其中ei是由模型得到的站点i的PM2.5浓度残差，
Figure BDA0003277183590000128
是平均值，cij是站点i和j之间的反距离空间权重。I的取值范围是-1～1，值越接近于0，残差空间自相关性越弱，模型越可靠。
步骤9：交叉验证。交叉验证采用的方法为8折交叉验证法，评估模型对非站点区域的PM2.5浓度的预测精度。将步骤1获得的PM2.5站点数据和自变量数据以及步骤3获得的筛选出的空间特征向量随机均分为8份，每次选7份作为训练集，按前述步骤进行建模，剩余的1份作为测试集，将测试集的自变量带入训练集所得模型中，计算测试集的均方根误差；每份数据都做过一次测试集后，计算8个均方根误差RMSE的均值，即为交叉验证的结果。RMSE越小，模型的预测精度越高，其鲁棒性越强，适用性越强。
步骤10：栅格计算。步骤6得到了地面PM2.5浓度的数学统计模型，也得到了整个研究区域各个位置点的各个因子的相关系数。将研究区域的各位置点转换为栅格影像，栅格值分别设为各因子的系数值。将提取的空间特征向量插值为覆盖整个研究区域的具有与AOD相同分辨率的栅格图像，然后进行栅格计算，用因子的相关系数栅格影像乘以因子栅格影像然后相加，最后加上截距项，得到地面PM2.5浓度估算值。
步骤11：异常值处理。在步骤10中得到的PM2.5浓度分布图，由于误差的存在，使得极少部分区域的PM2.5浓度值可能为负值，同时也可能存在少量的异常的高PM2.5浓度值，需要对这些异常值进行处理。负的PM2.5浓度值表征该点的PM2.5浓度低，故PM2.5浓度值设为0，同时给PM2.5浓度值设定一个阈值，超过该阈值的PM2.5浓度值均修改为该阈值，进行可视化展示，最终得到连续的地面PM2.5浓度分布图。
应当理解的是，上述针对本发明中较佳实施例的表述较为详细，但不能因此认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在本发明权利要求保护范围内，可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围以所附权利要求书为准。

Claims (8)

1.顾及全局空间自相关性和局部异质性的地面PM2.5浓度建模方法，其特征在于，包括以下步骤：
步骤1：数据获取与处理；获取气象站点的PM2.5地面监测值，并从各数据源获取与地面PM2.5相关的因子数据；气溶胶厚度产品与PM2.5浓度显著相关，为必需自变量；气象数据、污染源数据以及其他因子数据根据研究区气候特征和数据质量选取；对所有数据进行预处理，包括删除空值和异常值、投影转换、时间和空间分辨率的统一；
步骤2：构建全局空间权重矩阵；根据气象站点数据的投影坐标，构建研究区域内所有站点的全局空间权重矩阵，并中心化处理为对称矩阵，将空间关系转化为权重值；构建方法采用：基于邻接关系构建；或者基于反距离加权构建；或者复合型；
步骤3：全局空间因子提取；基于全局空间权重矩阵，根据空间自相关程度大小提取出全局空间影响因子，记为Eα；
步骤4：构建局部空间权重矩阵；选取空间权函数，构建每一个站点与该站点局部空间影响范围内其余站点的局部空间权重矩阵，局部空间影响范围在步骤5中优化确定；
步骤5：确定局部空间影响范围；局部空间影响范围也称为带宽，权函数对带宽的选择十分敏感，需要选择出最优的带宽；常用的最优带宽选择的方法有交叉验证法、AIC准则、贝叶斯信息准则以及平稳指数等；
步骤6：局部空间因子提取；根据局部空间权重矩阵和最优带宽，进行权重确定的局部最小二乘回归，提取出局部空间影响因子，记为βx；
步骤7：构建全局-局部空间回归模型；将步骤3中提取出的全局空间影响因子Eα和步骤6提取的局部空间影响因子βx进行回归建模,具体是将代表全局空间自相关的全局空间影响因子Eα和代表局部空间异质性的局部空间影响因子βx，利用全局空间因子和局部空间因子构建回归模型，构成如下模型：
Figure FDA0003277183580000011
其中，E是全局空间影响因子空间特征向量的集合，α为全局空间影响因子空间特征向量集的系数；(ui,vi)为第i个采样点的坐标系的横纵坐标，β0(ui,vi)为第i个采样点的回归方程的截距项；βk(ui,vi)为第i个采样点上第k个回归参数，xik是第i个采样点上第k个自变；αi是第i个采样点的空间特征向量集合的系数；εi为第i个区域的随机误差，满足零均值、同方差、相互独立等基本假定；
步骤8：评价与分析；计算模型精度评价指标；评价指标包括模型的拟合优度(R2)、调整后拟合优度(Adj.R2)、均方根误差(RMSE)、平均绝对百分比误差(MAPE)、平均绝对误差(MAE)以及残差的莫兰指数(Moran’s I)；并对模型进行鲁棒性检验，常采用交叉验证的方法计算上述指标。
2.如权利1要求所述的顾及全局空间自相关性和局部异质性的地面PM2.5浓度建模方法，其特征在于：步骤1中，统一因子的时间和空间分辨率，并进行投影转换，与AOD数据的时空分辨率和空间参考系一致；其中的污染源数据，如工厂位置、道路网密度等点、线数据，计算它们的核密度从而得到覆盖整个研究区域的栅格影像图。
3.如权利1要求所述的顾及全局空间自相关性和局部异质性的地面PM2.5浓度建模方法，其特征在于：步骤2具体包括
步骤2.1：
(1)基于邻接关系构建空间权重矩阵W1；构建站点的泰森多边形，可以根据Bishop邻接、Rook邻接和Queen邻接方式来构建全局空间权重矩阵；
(2)基于反距离加权构建空间权重矩阵W1；根据站点间的距离构建邻接关系，通过核函数将距离转为权重值，得到全局空间权重矩阵；核函数可分为指数、高斯、球状模型具体公式如下：
Figure FDA0003277183580000021
Figure FDA0003277183580000022
Figure FDA0003277183580000023
其中，公式(1)、(2)、(3)分别为指数、高斯、球状模型；i，j分别表示位置点i和位置点j；Wi,j表示位置点i和j之间的邻接性(权重)；r表示所有站点的最小生成树中边的最大距离；
(3)基于复合型构建空间权重矩阵W1；如K-最近邻矩阵和基于邻接关系和距离构建权重矩阵；K-最近邻矩阵，首先确定一个阈值K，搜寻每一个站点最近邻的K个站点，构建邻接矩阵，Wi,j表示位置点i和j之间的邻接性(权重)；基于邻接关系和距离构建权重矩阵，首先确定一个阈值d，搜寻每一个站点距离不超过阈值d的其他站点，认为两者邻接，Wi,j＝1，否则就不邻接，Wi,j＝0，构建权重矩阵；
步骤2.2：矩阵中心化；对上述步骤得到的矩阵W1进行中心化处理为对称矩阵W；中心化公式如下：
Figure FDA0003277183580000031
其中，I表示n维单位矩阵，11T是一个所有元素均为1的n×n矩阵，n表示研究区域内监测站点的数量。
4.如权利1要求所述的顾及全局空间自相关性和局部异质性的地面PM2.5浓度建模方法，其特征在于：步骤3具体包括
步骤3.1：计算空间权重矩阵W的特征值和特征向量，并进行初筛；求解特征值，得到特征值λi和特征向量Ei；初筛的条件有两个：一是其特征值λi>0；二是根据经验模型要求满足如下条件:
Figure FDA0003277183580000032
其中λi表示待筛选的特征值，λmax表示最大的特征值；初步筛选出符合条件的空间特征向量；
步骤3.2：空间特征向量筛选；根据逐步回归法筛选出使评价结果最优的空间特征向量；具体步骤如下：
步骤3.2.1：空间特征向量排序；把初筛后的空间特征向量按照特征值大小从小到大排列，得到空间特征向量集E；
步骤3.2.2：把空间特征向量集E的向量依次加入PM2.5与各个因子的普通最小二乘回归方程作为自变量，计算回归方程的AIC值；
步骤3.2.3：从步骤3.2.3中选择AIC值最小的回归方程Yk，对应的空间特征向量为Ek，然后从空间特征向量集E中将除Ek外的特征向量依次加入回归方程Yk，并计算回归方程的AIC值；
步骤3.2.4：判断AIC值是否明显减小，若AIC值明显减小，则重复步骤3.2.3；若AIC值几乎没减小则停止将特征向量加入回归方程，并记录已经加入到回归方程的特征向量；这些特征向量即为筛选出的使评价标准最优的特征向量集。
5.如权利1要求所述的顾及全局空间自相关性和局部异质性的地面PM2.5浓度建模方法，其特征在于：步骤4中，
空间权函数有距离阈值函数、指数函数、高斯函数、双重平方函数、三次立方函数，分别如公式6～11所示；
Figure FDA0003277183580000041
Figure FDA0003277183580000042
Figure FDA0003277183580000043
Figure FDA0003277183580000044
Figure FDA0003277183580000045
其中，Wi,j表示位置点i和j之间的邻接性(权重)，dij表示位置点i和j之间的距离，其中b是用于描述权重与距离之间函数关系的非负衰减参数，称为带宽(bandwidth)；带宽越大，权重随距离增加衰减得越慢；带宽越小，权重随距离增加衰减得越快。
6.如权利1要求所述的顾及全局空间自相关性和局部异质性的地面PM2.5浓度建模方法，其特征在于：步骤5中，采用交叉验证(cross validation,CV)法计算出最优带宽；
Figure FDA0003277183580000046
Figure FDA0003277183580000047
表示回归参数估计不包括回归点本身，即只根据回归点周围点进行回归计算；CV值最小的b则为最优带宽。
7.如权利1要求所述的顾及全局空间自相关性和局部异质性的地面PM2.5浓度建模方法，其特征在于：步骤8中交叉验证有两种方法；(1)K-折交叉验证；K值根据数据量来确定，K通常选择5、8或10；(2)子抽样验证；将数据分为两份，一份作为建模集，一份作为测试集。
8.如权利1要求所述的顾及全局空间自相关性和局部异质性的地面PM2.5浓度建模方法，其特征在于：步骤8具体包括:
步骤8.1：模型精度评价；计算所得模型的R2、调整后R2即Adj.R2、均方根误差(RMSE)、平均绝对误差(MAE)以及残差的莫兰指数(Moran’s I)等作为评价指标，以验证所提出的顾及全局空间自相关性和局部空间异质性的地面PM2.5浓度估算模型的精度；
Figure FDA0003277183580000051
其中yi是站点i的PM2.5浓度观测值，
Figure FDA0003277183580000052
是观测数据的平均值，
Figure FDA0003277183580000053
是模型预测的站点i的PM2.5浓度，n是监控站点的个数；
Figure FDA0003277183580000054
其中p是自变量的个数；R2和Adj.R2的取值范围是0～1,值越大说明模型精度越高；
Figure FDA0003277183580000055
Figure FDA0003277183580000056
式中参数含义同上，RMSE和MAE越小说明模型精度越高
Figure FDA0003277183580000057
其中ei是由模型得到的站点i的PM2.5浓度残差，
Figure FDA0003277183580000058
是平均值，cij是站点i和j之间的反距离空间权重；I的取值范围是-1～1，值越接近于0，残差空间自相关性越弱，模型越可靠；
步骤8.2：交叉验证；交叉验证采用的方法为8折交叉验证法，评估模型对非站点区域的PM2.5浓度的预测精度；将步骤1获得的PM2.5站点数据和自变量数据以及步骤3获得的筛选出的空间特征向量随机均分为8份，每次选7份作为训练集，按前述步骤进行建模，剩余的1份作为测试集，将测试集的自变量带入训练集所得模型中，计算测试集的均方根误差；每份数据都做过一次测试集后，计算8个均方根误差RMSE的均值，即为交叉验证的结果；RMSE越小，模型的预测精度越高，其鲁棒性越强，适用性越强；
步骤8.3：栅格计算；步骤7得到了地面PM2.5浓度的数学统计模型，也得到了整个研究区域各个位置点的各个因子的相关系数；将研究区域的各位置点转换为栅格影像，栅格值分别设为各因子的系数值；将提取的空间特征向量插值为覆盖整个研究区域的具有与AOD相同分辨率的栅格图像，然后进行栅格计算，用因子的相关系数栅格影像乘以因子栅格影像然后相加，最后加上截距项，得到地面PM2.5浓度估算值；
步骤8.4：异常值处理；在步骤10中得到的PM2.5浓度分布图，由于误差的存在，使得极少部分区域的PM2.5浓度值可能为负值，同时也可能存在少量的异常的高PM2.5浓度值，需要对这些异常值进行处理；负的PM2.5浓度值表征该点的PM2.5浓度低，故PM2.5浓度值设为0，同时给PM2.5浓度值设定一个阈值，超过该阈值的PM2.5浓度值均修改为该阈值，进行可视化展示，最终得到连续的地面PM2.5浓度分布图。
CN202111121091.0A 2021-09-24 2021-09-24 顾及全局空间自相关性和局部异质性的地面pm2.5浓度建模方法 Pending CN113901384A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111121091.0A CN113901384A (zh) 2021-09-24 2021-09-24 顾及全局空间自相关性和局部异质性的地面pm2.5浓度建模方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111121091.0A CN113901384A (zh) 2021-09-24 2021-09-24 顾及全局空间自相关性和局部异质性的地面pm2.5浓度建模方法

Publications (1)

Publication Number Publication Date
CN113901384A true CN113901384A (zh) 2022-01-07

Family

ID=79029333

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111121091.0A Pending CN113901384A (zh) 2021-09-24 2021-09-24 顾及全局空间自相关性和局部异质性的地面pm2.5浓度建模方法

Country Status (1)

Country Link
CN (1) CN113901384A (zh)

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114936957A (zh) * 2022-05-23 2022-08-23 福州大学 基于移动监测数据的城市pm25浓度分布模拟及场景解析模型
CN114974459A (zh) * 2022-05-25 2022-08-30 武汉大学 Pm2.5浓度估算模型的构建方法
CN115345075A (zh) * 2022-08-17 2022-11-15 北京城市气象研究院 一体化气溶胶污染气象指数-气溶胶浓度估算方法及系统
CN115599774A (zh) * 2022-12-15 2023-01-13 深圳市规划和自然资源数据管理中心（深圳市空间地理信息中心）（Cn） 基于局部时空树回归模型的时空非平稳性分析方法及系统
CN115630870A (zh) * 2022-11-01 2023-01-20 中国矿业大学 地质碳封存区域大气co2时空分异特征及影响因子分析方法

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114936957A (zh) * 2022-05-23 2022-08-23 福州大学 基于移动监测数据的城市pm25浓度分布模拟及场景解析模型
CN114974459A (zh) * 2022-05-25 2022-08-30 武汉大学 Pm2.5浓度估算模型的构建方法
CN114974459B (zh) * 2022-05-25 2024-04-16 武汉大学 Pm2.5浓度估算模型的构建方法
CN115345075A (zh) * 2022-08-17 2022-11-15 北京城市气象研究院 一体化气溶胶污染气象指数-气溶胶浓度估算方法及系统
CN115345075B (zh) * 2022-08-17 2023-04-18 北京城市气象研究院 一体化气溶胶污染气象指数-气溶胶浓度估算方法及系统
CN115630870A (zh) * 2022-11-01 2023-01-20 中国矿业大学 地质碳封存区域大气co2时空分异特征及影响因子分析方法
CN115630870B (zh) * 2022-11-01 2024-03-22 中国矿业大学 地质碳封存区域大气co2时空分异特征及影响因子分析方法
CN115599774A (zh) * 2022-12-15 2023-01-13 深圳市规划和自然资源数据管理中心（深圳市空间地理信息中心）（Cn） 基于局部时空树回归模型的时空非平稳性分析方法及系统

Similar Documents

Publication Publication Date Title
CN113901384A (zh) 顾及全局空间自相关性和局部异质性的地面pm2.5浓度建模方法
CN108761574B (zh) 基于多源信息融合的降雨量估算方法
US20220043182A1 (en) Spatial autocorrelation machine learning-based downscaling method and system of satellite precipitation data
CN111665575B (zh) 一种基于统计动力的中长期降雨分级耦合预报方法及系统
Mohammadi et al. Possibility investigation of tree diversity mapping using Landsat ETM+ data in the Hyrcanian forests of Iran
CN113496104B (zh) 基于深度学习的降水预报订正方法及系统
CN105974495B (zh) 利用分类拟合法预判目标区域未来平均云量的方法
CN108764643B (zh) 大范围作物病害风险评估方法
CN110927120B (zh) 一种植被覆盖度预警方法
CN111738175A (zh) 一种基于遥感影像和卷积神经网络的农业干旱监测系统
CN113704693B (zh) 一种高精度的有效波高数据估计方法
CN115357847B (zh) 一种基于误差分解的日尺度星地降水融合方法
CN108764527B (zh) 一种土壤有机碳库时空动态预测最优环境变量筛选方法
CN114463947B (zh) 一种基于时空网络卷积模型的对流性致灾强风预警预报方法
CN114781501A (zh) 一种基于主成分回归的多源降水融合方法
CN113901348A (zh) 一种基于数学模型的钉螺分布影响因素识别与预测方法
CN115544706A (zh) 一种小波和XGBoost模型集成的大气细颗粒物浓度估算方法
CN114254802B (zh) 气候变化驱动下植被覆盖时空变化的预测方法
CN114819737B (zh) 公路路域植被的碳储量估算方法、系统及存储介质
CN115598027A (zh) 一种基于遥感和机器学习技术的pm2.5反演方法
CN113191536A (zh) 基于机器学习的近地面环境要素预测模型训练和预测方法
CN114047563A (zh) 一种红外高光谱的全天候同化方法
CN113343783A (zh) 一种农作物智能识别与长势预测方法及系统
Neykov et al. Linking atmospheric circulation to daily precipitation patterns over the territory of Bulgaria
CN113049606A (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