一种基于克里金插值的PM2.5新增移动站点选址方法与流程

本发明属于地理信息空间分析技术领域，涉及一种基于克里金插值的pm2.5新增移动站点选址方法。

背景技术：

在对城市近地面pm2.5浓度进行监测管理时，需要通过地面监测站获取pm2.5浓度数据并进行分析，为pm2.5浓度掌控管理提供可靠依据。但是随着大气环境保护矛盾的日益突出和对pm2.5研究的迫切需要，监测站网其本身的缺点也日益显露。在理想条件下，监测点位越多，分布范围越广，所反映出的信息就越完整，但是pm2.5监测站点的布设、运行以及维护需要付出大量的人力、物力、财力，因此不可能将监测站点分布至城市的各个角落。而现如今的pm2.5监测站点大多数为人为随机选址，监测站点分布过于稀疏导致不能获取区域的关键信息，或在某一块区域站点分布过于密集导致信息冗余，从而使得监测站网的空间布局缺乏科学性。因此，pm2.5监测站点选址优化工作日益重要。目前只能通过对站点布局进行优化来表征出整个区域的大气环境质量，而目前的固定监测站点的位置是固定不变的，要想利用现有的固定站点的数据对监测站点布局进行优化，可以对新增站点的位置进行研究分析，寻到到最优的新增站点位置，从而提高监测站网整体的插值精度。

综上所述，利用当前pm2.5监测站点所获取的信息进行插值分析可能会对分析结果造成较大误差，通过新增站点可以提高整体监测站网络的插值精度，原始的方法不能寻找到新增移动站点的最佳位置，故提出一种新的站点优化方法。

相关参考文献如下：

[1]sophocleousm,paschettoje,oleara.ground‐waternetworkdesignfornorthwestkansas,usingthetheoryofregionalizedvariables[j].groundwater,1982,20(1):48-58.

[2]cesareld,myersde,posad.estimatingandmodelingspace-timecorrelationstructures[j].statistics&probabilityletters,2001,51(1):9-14.

[3]黄中伟.玉林市大气环境监测点位优化设置研究[j].沿海企业与科技,2004,(s1):42-44.

[4]段玉森,魏海萍,黄嵘,等.上海市tsp-pb空气质量监测点位优化研究[j].环境科学与管理,2007,23(9):131-134.

[5]李莎,舒红,徐正全.东北三省月降水量的时空克里金插值研究[j].水文,2011,31(3):31-35.

[6]sheny,wuy.optimizationofmarineenvironmentalmonitoringsitesintheyangtzeriverestuaryanditsadjacentsea,china[j].ocean&coastalmanagement,2013,73(73):92-100.

[7]秦怡雯,钱瑜,荣婷婷.基于大气特征污染物的监测布点选址优化研究[j].中国环境科学,2015,35(4):1056-1064.

[8]刘莎.空间数据与时空数据的分析方法及比较[d].长安大学,2015.

[9]倪敏捷.基于克里金插值法的海湾水环境监测网络优化研究——以福建泉州湾为例[d].厦门大学,2015.

[10]宣腾,李锦辉,李典庆,宋磊.基于普通克里金法的勘探位置优化方法[j].武汉大学学报(工学版),2016,49(05):714-719+739.

技术实现要素：

本发明的目的在于针对城市中固定监测站点的数目有限且分布缺乏科学性，如果只在这些站点返回的数据基础上进行分析，可能会对分析结果造成较大误差这一问题，提出一种基于克里金插值的pm2.5新增移动站点选址方法。本发明技术方案提供一种基于克里金插值的pm2.5新增移动站点选址方法，包括以下步骤：

步骤a，对初始数据进行前期分析，包括对已设置监测站点的pm2.5浓度数据的预处理、站点间距离统计以及研究区域格网的划分，将每个格网作为一个备选站点；

步骤b，计算普通克里金变异函数，包括选取研究区域某一时刻的若干已设置监测站点的数据并选取最佳变异函数模型进行计算；

步骤c，计算各个网格的pm2.5浓度的预测值，并计算各个网格的克里金估计误差方差，通过误差方差的原理计算研究区域内所有备选网格站点的克里金估计误差方差并划分等级进行比较；

步骤d，计算新增不同等级的备选网格站点后的整体区域方差，比较整体方差与其对应新增站点等级的关系，达到验证作用。

而且，所述步骤b的实现方式如下，

假设区域化变量z(x)满足本征假设和二阶平稳假设，其数学期望为m，协方差函数c(h)及变异函数γ(h)存在，表示如下，

e[z(x)]＝m

c(h)＝e[z(x)z(x+h)]-m²

式中，h表示区域化变量间的空间距离；

克里金插值法是通过估计区域变量z(x)在待插站点x0位置的值来近似接近真实值，预测公式如下，

式中，z^*(x0)是插值后的估计值，z(xi)是采样值，样本总数为n，λi是各空间样本点的权重系数，权重系数λi的计算必须满足以下两个条件，

①z^*(xi)是z(xi)的无偏估计量；

②估计方差最小，

估计方差用协方差函数表示如下，

式中，c(x0,x0)表示待插网格站点x0的协方差，c(xi,xj)表示已有监测站点间的协方差，c(xi,x0)表示已有监测站点与待插网格站点x0间的协方差；

要使估计方差最小，根据拉格朗日乘数法，令假设的拉格朗日函数式中，μ为拉格朗日乘子，为估计方差，令f对λi和μ的偏导数等于0，得到克里金方程组如下，

整理后得式(4)如下，

而且，步骤c实现方式如下，

在变异函数存在的条件下，根据协方差函数c(h)及变异函数γ(h)的关系式：c(h)＝c(0)-γ(h)，用变异函数表示的普通克里金方程组和克里金估计方差如下：

普通克里金方程组:

克里金估计误差方差:

式中，γ(xi,xj)为已有监测站点i和j之间的变异函数，γ(xj,x0)为已有监测站点j和待插网格站点x0之间的变异函数，μ为拉格朗日乘数；通过变换克里金方程组得到克里金方程矩阵，矩阵形式如(7)所示，

kλ＝d

解方程组(7)，求出权重系数λi和拉格朗日乘数μ，代入式(5)和(6)，求出估计值和估计方差。

而且，所述步骤d的实现方式如下，

为了验证新增移动站点的正确性，对研究区域内右上角的格网的克里金估计误差方差按照绝对值从小到大分为1、2、3、4四个等级，在每个等级分别选取站点进行新增；

分别从不同方差等级的区域选取格网作为新增站点并计算新增格网后整个研究区域格网的平均克里金估计误差方差通过对比新增不同方差等级的格网后的整个研究区域的普通克里金估计误差方差得出，在新增格网所在的方差等级越高，新增后区域的平均克里金估计误差的方差就越小，由初始数据经过普通克里金插值后区域内格网的克里金估计误差的方差绝对值的大小与研究区域内的新增移动站点的优先级成正比。

本发明利用对监测站点加密的方法，重点分析在现有固定监测站点的基础上，从新增移动站点和新增固定站点两个角度，对南昌市区pm2.5监测站网进行新增站点选址优化，从而利用较少的监测站点估计更准确的pm2.5浓度信息，从而更加有效地开展大气污染防治工作。本发明提出利用普通克里金插值法为研究区域找到最优的新增移动站点，结果具备一定的科学性，可以提高监测站网的预测精度。

附图说明

图1为本发明实施例的流程图；

图2为本发明涉及的克里金插值法的流程图；

图3为本发明实施例的变异函数模型曲线图。

具体实施方式

下面结合附图和实施例，进一步阐明本发明，应理解实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域的技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。

以下是与本发明相关的一些概念：

空间统计学是将具有地理空间信息特征的事物对象之间的相互作用及变化规律作为研究对象，结合统计学与现代图形的计算技术，通过较为直观的方式，来揭露空间数据中所隐含的空间模式、空间分布以及空间相互作用等一系列的数据特征。在地址统计学中，认为地质现象之间的属性关系在一定的范围内并不是相互独立存在的，而是具备一定的相关性，这种相关性不仅与样本间的距离有关，还与样本间的相对方向有关。

区域化变量，地质统计学学术用语。当一个变量呈现空间分布时，就称之为区域化。这种变量往往反映某种空间特质，用区域化变量来描述的现象，称之为区域化现象。亦称区域化随机变量。

区域化变量具有两个显著的特征：随机性和结构性。就随机性而言，区域化变量是一个随机函数，具有局部的、随机的、异常的性质；而对于结构性，区域化变量具有一般的或平均的结构性质，即变量在点x处的数值z(x)与偏离x空间距离为h的点x+h处的数值z(x+h)具有某种程度的自相关性，这种自相关依赖于两点间的距离向量h及变量特征。

为了更好地理解本发明的技术方案，下面实施例以真实基于克里金插值的新增移动站点选址实验来说明具体实施方式。参见图1，实施例提供的方法包括以下步骤：

1数据前期分析

数据的前期分析包括数据的预处理、站点间距离统计以及区域格网的划分三个部分。实施例中包括获取南昌市区9个国控监测站点的pm2.5浓度、pm2.5浓度数据数据预处理、9个国控监测站点的距离计算以及南昌市区区域网格划分。

实施例中的步骤1具体实现包括以下步骤：

1.1数据预处理

在对数据进行普通克里金插值分析之前，必须对数据进行前期处理，包括对已设置监测站点的pm2.5浓度数据的预处理、站点间距离统计以及研究区域格网的划分，将每个格网作为一个备选站点。在确定该时刻监测站点返回的数据无小于0或非常大的情况下，选取南昌市区9个大气环境监测站点2015年01月23日11：00的近地面pm2.5浓度数据，具体数据如下表1所示。

表1监测站点pm2.5浓度值

1.2站点距离分析

由于南昌市区一共包含9个大气监测站点且站点间最大距离为25km，所以在对变异函数进行计算的过程中，将这9个站点全部参与变异函数的计算，并按站点距离进行统计。这9个站点的分布对于整个研究区域而言具有良好的代表性，根据这9个站点的数据进行计算可以获得的拟合效果较好的变异函数。

在对变异函数的计算过程中，将站点距离进行分类，本文取步长为4km进行分类，共7类，即：

h≤4km

4km＜h≤8km

8km＜h≤12km

12km＜h≤16km

16km≤h≤20km

20km＜h≤24km

24km＜h≤28km

由于南昌目前的站点有9个，所以需要测算站点间的距离36个，具体结果如下表2所示。

表2南昌市区大气监测站点距离统计表

根据测算结果，将距离从0-28km划分为7个区间，具体数据如下表3所示。

表3站点间距离概况

从表中可以看出，对于南昌市区而言，9个大气监测站点间的距离与点对个数之间的关系可以反映出较好的代表性。

综上所述，可以采用2015年01月23日11：00南昌市区9个大气监测站点的数据对南昌市区进行插值分析。

1.3区域格网划分

在对南昌市区进行新增站点选址之前，需要对南昌市区整个区域进行格网化，将每个格网作为一个备选站点进行研究分析，这些格网的中心坐标即为备选站点的坐标。

本文划分出的格网区域的经度范围为115.7°-116.0°，纬度范围为28.5°-28.8°，将其均匀划分成一个31×31个小格网，每个小格网的大小大约为1km×1km，之所以将格网大致大小设为1km×1km是根据南昌市区的实际地域面积以及克里金插值算法的适宜距离等因素综合考虑而决定的。

在设计好格网后，对格网进行编号，便于后续寻找站点。从左下角一个的格网开始编号，在第一列依次向上递增，直至遍历每一列后，即可具体格网的编号与格网的经纬度信息对应起来。部分格网编号及坐标信息如下表4所示。

表4部分格网编号及坐标信息统计表

2计算普通克里金变异函数

步骤2中，计算普通克里金变异函数，包括选取南昌市区某一时刻的9个国控站点数据并选取最佳变异函数模型进行计算；

参见图2中的克里金插值法，本发明提出步骤2实现方式：

假设区域化变量z(x)满足本征假设和二阶平稳假设，其数学期望为m，协方差函数c(h)及变异函数γ(h)存在，即：

e[z(x)]＝m

c(h)＝e[z(x)z(x+h)]-m²

式中，h表示区域化变量间的空间距离，e表示数学期望，m表示常数。

克里金插值法是基于随机变量在空间和时间上的相关结构而建立的。其基本思想是：通过估计区域变量z(x)在待插网格站点x0位置的值来近似接近真实值，预测公式为式(1)。

式中，z^*(x0)是插值后的估计值，z(xi)是采样值，样本总数为n，λi是各空间样本点的权重系数，i＝1,2,…,n，权重系数λi的计算必须满足以下两个条件：

①z^*(xi)是z(xi)的无偏估计量，即时，有

②估计方差最小，也就是z^*(xi)和z(xi)之差的平方和最小。即

最小。

估计方差用协方差函数可表示为式(2)。

式中，c(x0,x0)表示待插网格站点的协方差，c(xi,xj)表示已有监测站点中第i个和第j个站点间的协方差，c(xi,x0)表示已有监测站点与待插网格站点间的协方差。

要使估计方差最小，根据拉格朗日乘数法，令式中，f即为假设的拉格朗日函数，μ为拉格朗日乘子，为估计方差。

令f对λi和μ的偏导数等于0，得到克里金方程组(3)

整理后得式(4)

实施例中，选取南昌市区某一时刻的9个国控站点数据并选取最佳变异函数模型进行计算。在对站点距离进行分类后，根据经验变异函数模型的计算原理，必须计算出站点之间的近地面pm2.5浓度数据的水平增量，各站点近地面pm2.5浓度的水平增量如下表5所示。

表5南昌市区站点pm2.5浓度水平增量统计表单位：ug/m3

tab.5statisticalofincreaseofpm2.5concentrationlevelinnanchangurbansiteunit:ug/m3

在计算过程中，根据计算出的变异函数值与拟合的函数模型之间的差异程度来决定最终使用的变异函数模型，通过实验变异函数值的计算公式所计算出的变异函数云图与球状模型最为吻合，所以利用初始数据选取球状模型来对研究区域进行插值分析是较为合适的，根据对数据点的拟合以及不断的调整参数，在综合考虑变异函数性质的情况下，得到近似最优的变异函数。变异函数公式表示为：

γ(h)＝537×[1.5(h/30)-0.5(h/30)³]

变异函数模型曲线图如图3所示。

3克里金方差计算

步骤3中，计算各个网格的pm2.5浓度的预测值，并计算各个网格的克里金估计误差方差。通过误差方差的原理计算研究区域内所有备选网格站点的克里金估计误差方差并划分等级进行比较。

步骤3实现方式如下，

在变异函数存在的条件下，根据协方差函数c(h)及变异函数γ(h)的关系式：c(h)＝c(0)-γ(h)，式中c(0)表示两点间距离趋于0时的协方差，为定值。用变异函数表示的普通克里金方程组和克里金估计方差，即：

普通克里金方程组:

克里金估计误差方差:

式(5)中γ(xi,xj)为国控点i和j之间的变异函数，式6中γ(xj,x0)为点j和待插站点之间的变异函数，μ为拉格朗日乘数。通过变换克里金方程组得到克里金方程矩阵，矩阵形式如(7)所示。

kλ＝d

式中，γij表示已有监测站点i和j之间的变异函数，γ(xi,x)表示已有监测站点i与任意待插网格站点x的变异函数，k为下方γij所构成的矩阵，d为下方γ(xi,x)所构成的矩阵，λi是各空间样本点的权重系数。

解方程组(7)，求出权重系数λi和拉格朗日乘数μ，代入式(5)和(6)，求出估计值和估计方差。

实施例中，在得到变异函数后，根据现有的南昌市区的9个国控监测站点对每个备选站点进行插值，通过普通克里金矩阵计算得到9个国控站点对于样本点的权重系数，计算各个网格的pm2.5浓度的预测值，并计算各个网格的克里金估计误差方差。通过误差方差的原理计算研究区域内所有备选网格站点的克里金估计误差方差并划分等级进行比较。

所述步骤3的实现方式如下：

在监测站点优化的评价中，主要是根据克里金估计误差的方差来对监测站点的分布进行优化。

由于克里金插值法中提供了克里金估计误差的方差，而这种特性正是克里金插值法的优势所在，其主要特点是：插值过程中，估计误差的方差大小只与已知样本点的数量、位置分布和变异结构有关，与已知样本点的数据值无关。克里金方差可以表示为：

式中，λi为各空间样本点的权重，γ(xi,x0)为已有监测站点i和待插网格点之间的变异函数，μ为拉格朗日乘数。

也就是说，当克里金估计误差的方差变小时，点位的分布更加合理，可以从研究区监测网络获取的信息也就越多。

4结果分析及验证

步骤4中，计算新增不同等级的备选网格站点后的整体区域方差，比较整体方差与其对应新增站点等级的关系，达到验证作用。

实施例中，步骤4包括以下步骤，

4.1克里金估计误差方差分析

通过对南昌市区9个大气监测站点2015年01月23日11：00的pm2.5浓度数据进行普通克里金插值并计算估计误差方差后，将计算出每个备选网格站点的克里金误差方差导入arcgis并生成克里金误差方差的等值线图，具体图文步骤在实审参考资料总提供。

4.2结果验证

为了验证新增移动站点的正确性，对研究区域内右上角的格网的克里金估计误差方差按照绝对值从小到大分为1、2、3、4四个等级，在每个等级分别选取站点进行新增，具体新增不同等级的区域分布图在实审参考资料总提供。

具体实施时，以上流程可采用软件技术实现自动运行。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。