CN112434935A - 一种可选的pm2.5浓度估算方法 - Google Patents

一种可选的pm2.5浓度估算方法 Download PDF

Info

Publication number
CN112434935A
CN112434935A CN202011315404.1A CN202011315404A CN112434935A CN 112434935 A CN112434935 A CN 112434935A CN 202011315404 A CN202011315404 A CN 202011315404A CN 112434935 A CN112434935 A CN 112434935A
Authority
CN
China
Prior art keywords
concentration
variation function
function
kriging
spatial
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
CN202011315404.1A
Other languages
English (en)
Other versions
CN112434935B (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.)
China Institute of Radio Wave Propagation CETC 22 Research Institute
Original Assignee
China Institute of Radio Wave Propagation CETC 22 Research Institute
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 China Institute of Radio Wave Propagation CETC 22 Research Institute filed Critical China Institute of Radio Wave Propagation CETC 22 Research Institute
Priority to CN202011315404.1A priority Critical patent/CN112434935B/zh
Publication of CN112434935A publication Critical patent/CN112434935A/zh
Application granted granted Critical
Publication of CN112434935B publication Critical patent/CN112434935B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q10/00—Administration; Management
    • G06Q10/06—Resources, workflows, human or project management; Enterprise or organisation planning; Enterprise or organisation modelling
    • G06Q10/067—Enterprise or organisation modelling
    • 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
    • G06F17/12—Simultaneous equations, e.g. systems of linear equations
    • 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/15—Correlation function computation including computation of convolution operations
    • 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/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q50/00—Information and communication technology [ICT] specially adapted for implementation of business processes of specific business sectors, e.g. utilities or tourism
    • G06Q50/06—Energy or water supply

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Business, Economics & Management (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Strategic Management (AREA)
  • Human Resources & Organizations (AREA)
  • Economics (AREA)
  • Algebra (AREA)
  • General Engineering & Computer Science (AREA)
  • Software Systems (AREA)
  • Databases & Information Systems (AREA)
  • Marketing (AREA)
  • General Business, Economics & Management (AREA)
  • Health & Medical Sciences (AREA)
  • Computing Systems (AREA)
  • Tourism & Hospitality (AREA)
  • Entrepreneurship & Innovation (AREA)
  • Operations Research (AREA)
  • Quality & Reliability (AREA)
  • Game Theory and Decision Science (AREA)
  • Educational Administration (AREA)
  • Development Economics (AREA)
  • Primary Health Care (AREA)
  • General Health & Medical Sciences (AREA)
  • Water Supply & Treatment (AREA)
  • Public Health (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种可选的PM2.5浓度估算方法，包括如下步骤：步骤1，通过离散变异函数公式计算所有站点PM2.5浓度样本点对的实验变异函数值；步骤2，采用最小二乘支持向量机拟合实验变异函数值，得到理论变异函数模型；步骤3，建立克里金空间插值方程组，根据理论变异函数模型求解克里金权重系数；步骤4，根据克里金权重系数计算待估位置点
Figure DEST_PATH_IMAGE002
处的PM2.5浓度估计值。本发明所公开一种可选的PM2.5浓度估算方法，克服了普通克里金空间插值PM2.5浓度估算方法中传统理论变异函数模型形状固定且未考虑空间变化趋势的问题，以最小二乘支持向量机拟合实验变异函数值，使得PM2.5浓度估算的结果符合本身的空间变化趋势。

Description

一种可选的PM2.5浓度估算方法
技术领域
本发明属于大气环境质量监测领域，特别涉及该领域中的一种可选的PM2.5浓度估算方法。
背景技术
细颗粒物PM2.5是指在空气动力学中直径小于或等于2.5μm的颗粒物，这种颗粒物能够长时间悬浮在空气中，被人体吸附进入支气管和肺泡中，直接影响肺的换气通气等功能，甚至也可以通过毛细血管进入人体血液循环系统，对心脏及心血管造成严重危害，近年来越发引发民众对大气环境质量问题的担忧。因此，如何准确地估算特定范围内的PM2.5浓度空间分布特征，实现PM2.5浓度的监测和预警，一直是相关领域研究的热点和难点。
当前，对于PM2.5浓度的估算主要通过站点观测的手段，即在特定区域内基于有限个PM2.5浓度观测站点数据进行空间插值反演，估算一定范围内的PM2.5浓度，这种方法得到的数据信息可靠性高，且能实现连续的动态观测。这里提及的空间插值是指根据已知观测样本点的PM2.5浓度真实值来估计待估位置点的PM2.5浓度值，其原理是通过已知观测点样本数据构建函数关系，综合空间位置关系以及空间相关性，从而估算其他任意点或任意分区的PM2.5浓度。
常用的空间插值法可以归结为二类，一类是确定性方法，另一类是空间统计方法。确定性方法中最具代表性的是反距离权重法和泰森多边形法。这种传统的概念性模型对PM2.5浓度的空间变化一般采用均化处理方式，用以方便产流计算，但随之带来的问题就是难以客观反映PM2.5浓度的空间分布，一定程度上忽略了PM2.5浓度在空间上的随机性，也未能有效考虑不同观测站点PM2.5浓度在空间上的相关性。空间统计方法中最具代表性的是克里金插值方法，该方法以区域化变量独有的随机性和结构性性质为基础，使得探索PM2.5浓度空间结构和空间变异规律变得有迹可循，同时以变异函数为基本工具对区域化变量的连续性、相关性、尺度性等要素进行空间描述。PM2.5浓度正是具有这种随机性(不确定性)与结构性(相关性)双重特征的区域化变量，应用克里金空间插值实现对PM2.5浓度估算其实质在于通过对已知位置点的PM2.5浓度内插或外推的方式，对待估位置点PM2.5浓度的取值进行无偏、最优估计。
克里金空间插值PM2.5浓度估算过程中，插值模型的精度取决于模型对空间变异性和空间相关性的反映程度。传统方法中插值模型需要用理论变异函数拟合实验变异函数，而理论变异函数模型往往是根据人的经验来选定，因此会导致不同的理论变异函数模型对PM2.5浓度估算结果的优劣程度产生较大的影响。此外，理论变异函数形状固定，有限的已知PM2.5浓度观测站点数据进行拟合时，无法反映实际PM2.5浓度在空间上的相关性和差异性，空间变化趋势被淹没。故如何根据现有PM2.5浓度观测站的数据生成空间计算单元内的PM2.5浓度，且尽可能地客观反映PM2.5浓度区域空间变化特征，是相关领域亟需解决的关键科学技术问题。
发明内容
本发明所要解决的技术问题就是提供一种可选的PM2.5浓度估算方法。
本发明采用如下技术方案：
一种可选的PM2.5浓度估算方法，其改进之处在于，包括如下步骤：
步骤1，通过离散变异函数公式计算所有站点PM2.5浓度样本点对的实验变异函数值；
将PM2.5浓度z(x)在站点x处和站点x+h处方差的一半定义为z(x)在x轴方向上的变异函数，记为γ(x，h)；
Figure BDA0002791198060000021
其中，Var[z(x)-z(x+h)]代表z(x)-z(x+h)的方差；
在二阶平稳假设的情况下，对任意h有：
E[z(x)]＝E[z(x+h)] (2)
其中，E[z(x)]和E[z(x+h)]分别代表z(x)和z(x+h)的PM2.5浓度数学期望，整理得：
Figure BDA0002791198060000022
这时变异函数γ(x，h)依赖于两个变量：站点位置x和PM2.5浓度样本点对的距离h，若变异函数仅仅依赖于距离h而与站点位置x无关，则γ(x，h)可以写为γ(h)：
Figure BDA0002791198060000023
离散样本的实验变异函数计算公式为：
Figure BDA0002791198060000024
其中，i＝1，...，N(h)，h为各站点PM2.5浓度样本点对的距离，N(h)代表站点PM2.5浓度样本点对距离为h时所有样本点对的个数，z(xi)和z(xi+h)分别代表在空间位置点xi和xi+h处的PM2.5浓度真实值；
步骤2，采用最小二乘支持向量机拟合实验变异函数值，得到理论变异函数模型；
步骤21，最小二乘支持向量机LS-SVM输入变量和输出变量分别为距离和实验变异函数值，假定需要拟合的数据集为
Figure BDA0002791198060000031
其中，hi∈Rd，取d＝1，hi表示第i个样本点对的距离，作为自变量，γ(hi)∈R，表示在距离hi下的实验变异函数值，作为因变量，n代表分组后需要拟合的实验变异函数值总个数；回归函数f(h)的基本形式可以用下式(6)表示：
Figure BDA0002791198060000032
其中，h为PM2.5浓度样本点对的距离，ω为权系数向量，即列向量，ωT代表其转置向量；
Figure BDA0002791198060000033
为输入空间到特征空间的映射函数，b为常数项；
步骤22，依据统计学习理论，支持向量机模型的目的是使结构风险和经验风险同时达到最小，将支持向量机模型转换为优化函数
Figure BDA0002791198060000034
Figure BDA0002791198060000035
其中，i＝1，...，N，N为分组后需要拟合的实验变异函数值总个数，
Figure BDA0002791198060000036
表示结构风险，结构风险描述支持向量机模型的复杂度，||ω||2为ω的2-范数平方，
Figure BDA0002791198060000037
Figure BDA0002791198060000038
表示经验风险，经验风险描述支持向量机模型与真实数据的拟合程度，在最小二乘支持向量机方法中，经验风险用误差平方和表示，ei表示误差项，γ表示正则化参数；
在LS-SVM中，误差项ei等于站点PM2.5浓度真实值yi与回归模型计算的PM2.5浓度
Figure BDA0002791198060000039
之差，因此优化函数须满足约束条件：
Figure BDA00027911980600000310
步骤23，利用拉格朗日乘数法将步骤22中式(8)含约束条件的优化函数转化为无约束条件的拉格朗日函数，所述拉格朗日函数
Figure BDA00027911980600000311
Figure BDA00027911980600000312
整理为：
Figure BDA00027911980600000313
其中，αi是拉格朗日乘子，根据KKT条件，所述拉格朗日函数最优解条件为：
Figure BDA00027911980600000314
用线性方程组表示上式，得到：
Figure BDA0002791198060000041
其中，
Figure BDA0002791198060000042
e＝[e1，e2，…，eN]T，y＝[y1，y2，…，yN]T，
Figure BDA0002791198060000043
α＝[α1，α2，…，αN]T，I为单位阵；
消去变量ω和e，并结合Mercer条件定义核函数：
Figure BDA0002791198060000044
其中，K(hi，hj)是核函数；
化简，得到线性方程组为：
Figure BDA0002791198060000045
其中，定义Ω为：
Figure BDA0002791198060000046
线性方程组写成矩阵形式：
Figure BDA0002791198060000047
记
Figure BDA0002791198060000048
由于A1是对称半正定矩阵，
Figure BDA0002791198060000049
存在，则上述方程组的解如下：
Figure BDA00027911980600000410
Figure BDA00027911980600000411
最后得到最小二乘支持向量回归模型：
Figure BDA00027911980600000412
其中，αi表示拉格朗日乘子，K(hi，h)表示核函数，其中，i＝1，...，N，N为分组后需要拟合的实验变异函数值总个数，b为常数项；
步骤3，建立克里金空间插值方程组，根据理论变异函数模型求解克里金权重系数λi；
克里金插值方法需要满足无偏性和估计方差最小；
E[z(x0)-z*(x0)]＝0 (20)
Var[z(x0)-z*(x0)]＝min (21)
其中，z(x0)为待估位置点x0处的PM2.5浓度真实值，z*(x0)为待估位置点x0处的PM2.5浓度估计值，E[z(x0)-z*(x0)]代表z(x0)-z*(x0)的数学期望，Var[z(x0)-z*(x0)]代表z(x0)-z*(x0)的方差；
无偏性表现为：
Figure BDA0002791198060000051
即：
Figure BDA0002791198060000052
估计方差最小表现为：
Figure BDA0002791198060000053
根据拉格朗日乘数原理求条件极值：
Figure BDA0002791198060000054
其中，μ为克里金拉格朗日乘子，求F对λi和μ的偏导数，经整理，用变异函数来表示，即：
Figure BDA0002791198060000055
其中，k为PM2.5浓度样本观测点的总个数，λi为克里金权重系数，表示各空间位置点xi处的PM2.5浓度z(xi)对待估位置点x0的贡献程度，μ为克里金拉格朗日乘子，γ(xi，xj)为空间位置点xi与xj距离下的实验变异函数值，γ(x0，xj)为待估位置点x0与空间位置点xj距离下的实验变异函数值；
上述公式用矩阵表示如下：
Figure BDA0002791198060000056
其中，λ为列向量，代表λi，i＝1，...，k，普通克里格方程组为：
K*λ＝D (28)
解得：
λ＝K-1D (29)
步骤4，根据克里金权重系数λi计算待估位置点x0处的PM2.5浓度估计值z*(x0)：
Figure BDA0002791198060000061
本发明的有益效果是：
本发明所公开一种可选的PM2.5浓度估算方法，克服了普通克里金空间插值PM2.5浓度估算方法中传统理论变异函数模型形状固定且未考虑空间变化趋势的问题，以最小二乘支持向量机拟合实验变异函数值，使得PM2.5浓度估算的结果符合本身的空间变化趋势。本发明直观的效果体现为高精度的插值估算结果，即准确地估计了PM2.5浓度的时空分布特征，提高了PM2.5浓度估算的精度和可靠性，解决了PM2.5浓度估算的问题。对于城市大气环境监测及防护具有重要意义，更深层次的技术效果体现为提高地理数据解释地理现象的能力，这是极具规律性的。
附图说明
图1是本发明方法的流程示意图；
图2是青岛市19个国控空气质量观测站点的分布示意图；
图3是青岛市19个观测站点2020年8月份的PM2.5浓度平均值；
图4是通过普通克里金方法得到的PM2.5浓度估算效果图；
图5是通过本发明方法得到的PM2.5浓度估算效果图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图和实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
实施例1，如图1所示，本实施例公开了一种可选的PM2.5浓度估算方法，包括如下步骤：
步骤1，通过离散变异函数公式计算所有站点PM2.5浓度样本点对的实验变异函数值；
将PM2.5浓度z(x)在站点x处和站点x+h处方差的一半定义为z(x)在x轴方向上的变异函数，记为γ(x，h)；
Figure BDA0002791198060000062
其中，Var[z(x)-z(x+h)]代表z(x)-z(x+h)的方差；
在二阶平稳假设的情况下，对任意h有：
E[z(x)]＝E[z(x+h)] (2)
其中，E[z(x)]和E[z(x+h)]分别代表z(x)和z(x+h)的PM2.5浓度数学期望，整理得：
Figure BDA0002791198060000063
这时变异函数γ(x，h)依赖于两个变量：站点位置x和PM2.5浓度样本点对的距离h，若变异函数仅仅依赖于距离h而与站点位置x无关，则γ(x，h)可以写为γ(h)：
Figure BDA0002791198060000071
离散变异函数γ(h)是克里金空间插值独有的参量，计算公式如下式(5)所述，通过对所有站点PM2.5浓度样本点对的计算，得到该参量值。
离散样本的实验变异函数计算公式为：
Figure BDA0002791198060000072
其中，i＝1，...，N(h)，h为各站点PM2.5浓度样本点对的距离，N(h)代表站点PM2.5浓度样本点对距离为h时所有样本点对的个数，z(xi)和z(xi+h)分别代表在空间位置点xi和xi+h处的PM2.5浓度真实值；
实验变异函数值拟合之前，若需要拟合的实验变异函数值较多，则先进行分组操作，以便降低计算过程中的复杂度。由于本实施例样本点数目较多，拟合实验变异函数之前，进行分组操作，分组后需要拟合的实验变异函数值总个数为20，即n＝20。
步骤2，采用最小二乘支持向量机拟合实验变异函数值，得到理论变异函数模型，通过支持向量机与克里金空间插值相结合，创新性思维解决了PM2.5浓度估算过程中对所有待估位置点PM2.5浓度无偏、最优的估计问题，这也是本发明区别于普通的PM2.5浓度估算方法最明显的特征，具体过程如下：
步骤21，最小二乘支持向量机(简称LS-SVM)输入变量和输出变量分别为距离和实验变异函数值，假定需要拟合的数据集为
Figure BDA0002791198060000073
其中，hi∈Rd，取d＝1，hi表示第i个样本点对的距离，作为自变量，γ(hi)∈R，表示在距离hi下的实验变异函数值，作为因变量，n代表分组后需要拟合的实验变异函数值总个数；根据最小二乘支持向量机模型，该模型通过求解线性方程组解决回归问题，回归函数f(h)的基本形式可以用下式(6)表示：
Figure BDA0002791198060000074
其中，h为PM2.5浓度样本点对的距离，ω为权系数向量(列向量)，ωT代表其转置向量；
Figure BDA0002791198060000075
为输入空间到特征空间的映射函数，b为常数项；
步骤22，依据统计学习理论，支持向量机模型的目的是使结构风险和经验风险同时达到最小，将支持向量机模型转换为优化函数
Figure BDA0002791198060000076
Figure BDA0002791198060000077
其中，i＝1，...，N，N为分组后需要拟合的实验变异函数值总个数，
Figure BDA0002791198060000081
表示结构风险，结构风险描述支持向量机模型的复杂度，||ω||2为ω的2-范数平方，
Figure BDA0002791198060000082
Figure BDA0002791198060000083
表示经验风险，经验风险描述支持向量机模型与真实数据的拟合程度，在最小二乘支持向量机方法中，经验风险用误差平方和表示，ei表示误差项，γ表示正则化参数；
在LS-SVM中，误差项ei等于站点PM2.5浓度真实值yi与回归模型计算的PM2.5浓度
Figure BDA0002791198060000084
之差，因此优化函数须满足约束条件：
Figure BDA0002791198060000085
步骤23，利用拉格朗日乘数法将步骤22中式(8)含约束条件的优化函数转化为无约束条件的拉格朗日函数，所述拉格朗日函数
Figure BDA0002791198060000086
Figure BDA0002791198060000087
整理为：
Figure BDA0002791198060000088
其中，αi是拉格朗日乘子，根据KKT条件，所述拉格朗日函数最优解条件为：
Figure BDA0002791198060000089
用线性方程组表示上式，得到：
Figure BDA00027911980600000810
其中，
Figure BDA00027911980600000811
e＝[e1，e2，…，eN]T，y＝[y1，y2，…，yN]T，
Figure BDA00027911980600000812
α＝[α1，α2，…，αN]T，I为单位阵；
消去变量ω和e，并结合Mercer条件定义核函数：
Figure BDA00027911980600000813
其中，K(hi，hj)是核函数，可取线性、多项式、高斯等核函数；
化简，得到线性方程组为：
Figure BDA0002791198060000091
其中，定义Ω为：
Figure BDA0002791198060000092
线性方程组写成矩阵形式：
Figure BDA0002791198060000093
记
Figure BDA0002791198060000094
由于A1是对称半正定矩阵，
Figure BDA0002791198060000095
存在，则上述方程组的解如下：
Figure BDA0002791198060000096
Figure BDA0002791198060000097
根据拉格朗日乘数法求解，最后得到最小二乘支持向量回归模型：
Figure BDA0002791198060000098
其中，αi表示拉格朗日乘子，K(hi，h)表示核函数，其中，i＝1，...，N，N为分组后需要拟合的实验变异函数值总个数，b为常数项；
本实施例中，最小二乘支持向量机核函数采用了RBF核函数，惩罚参数和核参数优化都是通过交叉验证的方式进行优化，最终得到惩罚参数γ＝3.18×103，核参数σ2＝59，常数b＝13.98。
步骤3，建立克里金空间插值方程组，根据理论变异函数模型求解克里金权重系数λi；
以普通克里金插值方法的原理进行说明，普通克里金插值方法估计需要满足两个条件：无偏性和估计方差最小；
E[z(x0)-z*(x0)]＝0 (20)
Var[z(x0)-z*(x0)]＝min (21)
其中，z(x0)为待估位置点x0处的PM2.5浓度真实值，z*(x0)为待估位置点x0处的PM2.5浓度估计值，E[z(x0)-z*(x0)]代表z(x0)-z*(x0)的数学期望，Var[z(x0)-z*(x0)]代表z(x0)-z*(x0)的方差；
无偏性表现为：
Figure BDA0002791198060000101
即：
Figure BDA0002791198060000102
估计方差最小表现为：
Figure BDA0002791198060000103
根据拉格朗日乘数原理求条件极值：
Figure BDA0002791198060000104
其中，μ为克里金拉格朗日乘子，求F对λi和μ的偏导数，经整理，用变异函数来表示，即普通克里金插值方法在满足无偏性和估计方差最小这两个原则的条件下，建立了含有约束条件的拉格朗日函数。其中，无偏性体现在约束条件上，而方差最小体现在求极值问题上。通过拉格朗日乘数法求解，得到克里金插值方程组可表示为：
Figure BDA0002791198060000105
其中，k为PM2.5浓度样本观测点的总个数，λi为克里金权重系数，表示各空间位置点xi处的PM2.5浓度z(xi)对待估位置点x0的贡献程度，μ为克里金拉格朗日乘子，γ(xi，xj)为空间位置点xi与xj距离下的实验变异函数值，γ(x0，xj)为待估位置点x0与空间位置点xj距离下的实验变异函数值；
上述公式用矩阵表示如下：
Figure BDA0002791198060000106
其中，λ为列向量，代表λi，i＝1，...，k，通过拟合得到的理论变异函数模型，解方程求解克里金权重系数λi。普通克里格方程组为：
K*λ＝D (28)
解得：
λ＝K-1D (29)
步骤4，根据克里金权重系数λi计算待估位置点x0处的PM2.5浓度估计值z*(x0)，即实现待估位置点PM2.5浓度无偏、最优的估计，也是本发明所要达到的目的及技术效果，可用下式(30)进行计算：
Figure BDA0002791198060000111
青悦开放环境数据中心提供的山东省青岛市19个国控空气质量观测站点2020年8月份PM2.5浓度月平均值(单位：μg/m3)，观测站点包括青岛市区9个、西海岸新区2个、即墨区2个、平度市2个、胶州市2个、莱西市2个，其分布如图2所示。青岛市19个观测站点2020年8月份PM2.5浓度平均值如图3所示。
在已知某一区域范围内有限个观测站点PM2.5浓度值的情况下，应用克里金空间插值实现对PM2.5浓度估算的实质在于通过对已知位置点的PM2.5浓度内插或外推的方式，对待估位置点PM2.5浓度的取值进行无偏、最优估计。本发明中提到的一种可选的PM2.5浓度估算方法，从已知PM2.5浓度样本数据点的空间变化趋势出发，通过LS-SVM拟合实验变异函数，这种方式拟合的结果符合PM2.5浓度本身的空间变化趋势，解决了传统方法无法反映实际空间变化趋势的问题，实现了通过现有PM2.5浓度观测站数据生成区域空间计算单元内的PM2.5浓度，并尽可能地客观反映PM2.5浓度区域空间变化特征。
为验证本发明的准确性，采用了普通克里金方法和本发明方法两种方法进行对比验证。图4表示通过普通克里金方法得到的PM2.5浓度估算效果图。图5表示通过本发明方法得到的PM2.5浓度估算效果图。从两图的对比中可以看出，采用本发明得到的结果与采用普通克里金方法得到的PM2.5浓度整体趋势一致，且得到的PM2.5浓度都在13μg/m3～30μg/m3的区间范围内。但更为细致地是，本发明得到的PM2.5浓度呈现东南和东北方向偏低、西北和西南方向偏高的结果，且整体上的空间分布特征从东向西PM2.5浓度由低到高，表现为逐渐上升的趋势；同时，本发明得到的PM2.5浓度空间层次感变化明显，且存在特征变化急剧的位置，体现了LS-SVM模拟拟合变异函数中的空间变化效应。因此，认为本发明具有一定的科学性和准确性，根据其时空分布特征来重点防护与治理是进行研究的意义所在。
综上所述，为详细阐述本发明方法的具体过程，以山东省青岛市19个国控空气质量监测站点2020年8月份PM2.5浓度月平均值数据为例进行说明，通过克里金空间插值实现PM2.5浓度的估算是本发明要解决的典型问题，问题本身契合具体的技术发明领域，具备了良好的代表性，因此选择此数据作为本技术方案的实施例。同时，本发明的技术方案不能脱离技术流程和数学模型而独立存在，各个步骤出现的计算和公式是实现本发明方案的必不可少的技术手段，其技术特征已经体现在数学函数参数所代表的物理意义及具体的应用领域。

Claims (1)

1.一种可选的PM2.5浓度估算方法，其特征在于，包括如下步骤：
步骤1，通过离散变异函数公式计算所有站点PM2.5浓度样本点对的实验变异函数值；
将PM2.5浓度z(x)在站点x处和站点x+h处方差的一半定义为z(x)在x轴方向上的变异函数，记为γ(x，h)；
Figure FDA0002791198050000011
其中，Var[z(x)-z(x+h)]代表z(x)-z(x+h)的方差；
在二阶平稳假设的情况下，对任意h有：
E[z(x)]＝E[z(x+h)] (2)
其中，E[z(x)]和E[z(x+h)]分别代表z(x)和z(x+h)的PM2.5浓度数学期望，整理得：
Figure FDA0002791198050000012
这时变异函数γ(x，h)依赖于两个变量：站点位置x和PM2.5浓度样本点对的距离h，若变异函数仅仅依赖于距离h而与站点位置x无关，则γ(x，h)可以写为γ(h)：
Figure FDA0002791198050000013
离散样本的实验变异函数计算公式为：
Figure FDA0002791198050000014
其中，i＝1，...，N(h)，h为各站点PM2.5浓度样本点对的距离，N(h)代表站点PM2.5浓度样本点对距离为h时所有样本点对的个数，z(xi)和z(xi+h)分别代表在空间位置点xi和xi+h处的PM2.5浓度真实值；
步骤2，采用最小二乘支持向量机拟合实验变异函数值，得到理论变异函数模型；
步骤21，最小二乘支持向量机LS-SVM输入变量和输出变量分别为距离和实验变异函数值，假定需要拟合的数据集为
Figure FDA0002791198050000015
其中，hi∈Rd，取d＝1，hi表示第i个样本点对的距离，作为自变量，γ(hi)∈R，表示在距离hi下的实验变异函数值，作为因变量，n代表分组后需要拟合的实验变异函数值总个数；回归函数f(h)的基本形式可以用下式(6)表示：
Figure FDA0002791198050000016
其中，h为PM2.5浓度样本点对的距离，ω为权系数向量，即列向量，ωT代表其转置向量；
Figure FDA0002791198050000017
为输入空间到特征空间的映射函数，b为常数项；
步骤22，依据统计学习理论，支持向量机模型的目的是使结构风险和经验风险同时达到最小，将支持向量机模型转换为优化函数
Figure FDA0002791198050000018
Figure FDA0002791198050000019
其中，i＝1，...，N，N为分组后需要拟合的实验变异函数值总个数，
Figure FDA0002791198050000021
表示结构风险，结构风险描述支持向量机模型的复杂度，||ω||2为ω的2-范数平方，
Figure FDA0002791198050000022
Figure FDA0002791198050000023
表示经验风险，经验风险描述支持向量机模型与真实数据的拟合程度，在最小二乘支持向量机方法中，经验风险用误差平方和表示，ei表示误差项，γ表示正则化参数；
在LS-SVM中，误差项ei等于站点PM2.5浓度真实值yi与回归模型计算的PM2.5浓度
Figure FDA0002791198050000024
之差，因此优化函数须满足约束条件：
Figure FDA0002791198050000025
步骤23，利用拉格朗日乘数法将步骤22中式(8)含约束条件的优化函数转化为无约束条件的拉格朗日函数，所述拉格朗日函数
Figure FDA0002791198050000026
Figure FDA0002791198050000027
整理为：
Figure FDA0002791198050000028
其中，αi是拉格朗日乘子，根据KKT条件，所述拉格朗日函数最优解条件为：
Figure FDA0002791198050000029
用线性方程组表示上式，得到：
Figure FDA00027911980500000210
其中，
Figure FDA00027911980500000211
e＝[e1，e2，…，eN]T，y＝[y1，y2，…，yN]T，
Figure FDA00027911980500000212
α＝[α1，α2，…，αN]T，I为单位阵；
消去变量ω和e，并结合Mercer条件定义核函数：
Figure FDA00027911980500000213
其中，K(hi，hj)是核函数；
化简，得到线性方程组为：
Figure FDA0002791198050000031
其中，定义Ω为：
Figure FDA0002791198050000032
线性方程组写成矩阵形式：
Figure FDA0002791198050000033
记
Figure FDA0002791198050000034
由于A1是对称半正定矩阵，
Figure FDA0002791198050000035
存在，则上述方程组的解如下：
Figure FDA0002791198050000036
Figure FDA0002791198050000037
最后得到最小二乘支持向量回归模型：
Figure FDA0002791198050000038
其中，αi表示拉格朗日乘子，K(hi，h)表示核函数，其中，i＝1，...，N，N为分组后需要拟合的实验变异函数值总个数，b为常数项；
步骤3，建立克里金空间插值方程组，根据理论变异函数模型求解克里金权重系数λi；
克里金插值方法需要满足无偏性和估计方差最小；
E[z(x0)-z*(x0)]＝0 (20)
Var[z(x0)-z*(x0)]＝min (21)
其中，z(x0)为待估位置点x0处的PM2.5浓度真实值，z*(x0)为待估位置点x0处的PM2.5浓度估计值，E[z(x0)-z*(x0)]代表z(x0)-z*(x0)的数学期望，Var[z(x0)-z*(x0)]代表z(x0)-z*(x0)的方差；
无偏性表现为：
Figure FDA0002791198050000039
即：
Figure FDA00027911980500000310
估计方差最小表现为：
Figure FDA0002791198050000041
根据拉格朗日乘数原理求条件极值：
Figure FDA0002791198050000042
其中，μ为克里金拉格朗日乘子，求F对λi和μ的偏导数，经整理，用变异函数来表示，即：
Figure FDA0002791198050000043
其中，k为PM2.5浓度样本观测点的总个数，λi为克里金权重系数，表示各空间位置点xi处的PM2.5浓度z(xi)对待估位置点x0的贡献程度，μ为克里金拉格朗日乘子，γ(xi，xj)为空间位置点xi与xj距离下的实验变异函数值，γ(x0，xj)为待估位置点x0与空间位置点xj距离下的实验变异函数值；
上述公式用矩阵表示如下：
Figure FDA0002791198050000044
其中，λ为列向量，代表λi，i＝1，...，k，普通克里格方程组为：
K*λ＝D (28)
解得：
λ＝K-1D (29)
步骤4，根据克里金权重系数λi计算待估位置点x0处的PM2.5浓度估计值z*(x0)：
Figure FDA0002791198050000045
CN202011315404.1A 2020-11-20 2020-11-20 一种可选的pm2.5浓度估算方法 Active CN112434935B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011315404.1A CN112434935B (zh) 2020-11-20 2020-11-20 一种可选的pm2.5浓度估算方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011315404.1A CN112434935B (zh) 2020-11-20 2020-11-20 一种可选的pm2.5浓度估算方法

Publications (2)

Publication Number Publication Date
CN112434935A true CN112434935A (zh) 2021-03-02
CN112434935B CN112434935B (zh) 2022-04-22

Family

ID=74693417

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011315404.1A Active CN112434935B (zh) 2020-11-20 2020-11-20 一种可选的pm2.5浓度估算方法

Country Status (1)

Country Link
CN (1) CN112434935B (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116580123A (zh) * 2023-03-29 2023-08-11 云南省生态环境科学研究院 一种区域大气环境光化学网格模型图像生成系统
CN116913410A (zh) * 2023-07-12 2023-10-20 中国科学院地理科学与资源研究所 一种元素浓度值确定方法和装置

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106404620A (zh) * 2015-07-30 2017-02-15 中国科学院遥感与数字地球研究所 地统计插值与卫星遥感联合反演地面pm2.5的方法及系统
CN106407633A (zh) * 2015-07-30 2017-02-15 中国科学院遥感与数字地球研究所 基于时空回归克里金模型估算地面pm2.5的方法及系统
CN109523066A (zh) * 2018-10-29 2019-03-26 东华理工大学 一种基于克里金插值的pm2.5新增移动站点选址方法
CN109856021A (zh) * 2018-12-24 2019-06-07 天津珞雍空间信息研究院有限公司 一种pm2.5反演方法及监测区域分割方法
CN110766191A (zh) * 2019-08-27 2020-02-07 东华理工大学 一种基于时空克里金插值的新增pm2.5固定监测站点选址方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106404620A (zh) * 2015-07-30 2017-02-15 中国科学院遥感与数字地球研究所 地统计插值与卫星遥感联合反演地面pm2.5的方法及系统
CN106407633A (zh) * 2015-07-30 2017-02-15 中国科学院遥感与数字地球研究所 基于时空回归克里金模型估算地面pm2.5的方法及系统
CN109523066A (zh) * 2018-10-29 2019-03-26 东华理工大学 一种基于克里金插值的pm2.5新增移动站点选址方法
CN109856021A (zh) * 2018-12-24 2019-06-07 天津珞雍空间信息研究院有限公司 一种pm2.5反演方法及监测区域分割方法
CN110766191A (zh) * 2019-08-27 2020-02-07 东华理工大学 一种基于时空克里金插值的新增pm2.5固定监测站点选址方法

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116580123A (zh) * 2023-03-29 2023-08-11 云南省生态环境科学研究院 一种区域大气环境光化学网格模型图像生成系统
CN116913410A (zh) * 2023-07-12 2023-10-20 中国科学院地理科学与资源研究所 一种元素浓度值确定方法和装置
CN116913410B (zh) * 2023-07-12 2024-05-17 中国科学院地理科学与资源研究所 一种元素浓度值确定方法和装置

Also Published As

Publication number Publication date
CN112434935B (zh) 2022-04-22

Similar Documents

Publication Publication Date Title
CN112434935B (zh) 一种可选的pm2.5浓度估算方法
CN107908913B (zh) 基于并行计算机的地球动力数模方法
CN109726433B (zh) 基于曲面边界条件的三维无粘低速绕流的数值模拟方法
CN111833450A (zh) 一种融合有限元分析法的超声波三维快速重建及分析方法
CN107204040A (zh) 多点地质统计学建模方法及装置、计算机存储介质
CN112083453A (zh) 一种涉及水汽时空参数的对流层层析方法
CN112257021A (zh) 一种可选的克里金空间插值降雨量估算方法
Cecil et al. Simplex free adaptive tree fast sweeping and evolution methods for solving level set equations in arbitrary dimension
Xue et al. Bias estimation and correction for triangle-based surface area calculations
CN107621637B (zh) 基于单多普勒雷达的切变区风场反演方法
Aiton A radial basis function partition of unity method for transport on the sphere
US9633165B1 (en) Composite design direction
CN113496097B (zh) 一种基于sph的飞机油箱燃油晃动仿真分析方法
CN112507282B (zh) 一种基于速度梯度张量特性的流动显示方法
CN115754007A (zh) 一种基于声发射技术和层析成像技术的损伤检测方法
Shure et al. Physical constraints for the analysis of the geomagnetic secular variation
CN103971411A (zh) 利用三维物体的空间曲面采样点对空间曲面建模的方法
Kosnikov Increasing the information capacity of the interface of the control system for multyparameter objects
Lima et al. A consistent mass-conserving C-staggered method for shallow water equations on global reduced grids
RU2812143C1 (ru) Способ и устройство для измерения характеристик колонки породы для создания модели поровой системы
CN116186859B (zh) 一种大跨索屋盖面组分区的风振动力计算方法及其系统
Hollister et al. Bivariate quantile interpolation for ensemble derived probability density estimates
Pedder et al. The semi‐geostrophic diagnosis of vertical motion. I: Formulation and coordinate transformations
CN118013870B (zh) 一种基于分布函数思想求解ns方程的深度神经网络方法
CN113607920B (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