CN103336093A - 一种区域空间质量分析方法 - Google Patents

Abstract

本发明涉及空间质量分析相关技术领域，特别是涉及一种区域空间质量分析方法，包括：对所述区域进行网格划分，并确定每个网格中心点的经纬度坐标；获取单个时间片上的大气环境监测数据；确定以每个网格中心点为插值点的影响点搜索范围；以所述影响点搜索范围内的监测站点作为影响点搜索范围的影响样本点，确定每个影响样本点的权重值；对每个插值点采用对应的影响样本点的污染物浓度监测值，根据对应的权重值进行插值计算，得到与每个插值点对应的网格中心点的浓度插值结果。本发明的区域空间质量分析方法，适用于中小尺度的区域空气质量空间分析，保证计算速度快的同时，插值精度更高，为公众提供更为准确的、直观的区域污染物空间分布信息。

Description

一种区域空间质量分析方法

技术领域

本发明涉及空间质量分析相关技术领域，特别是涉及一种区域空间质量分析方法。

背景技术

环境问题已成为当今人类发展所面临的最严峻挑战之一。环保部门在环境监测设备维护上投入大量人力和物力，但仍疲于奔命。各地大气环境监测站通常离散、不均匀、数量有限，直接采用在各个站点上标注数据的方式，显示速度快，数据准确，但对区域空气质量的显示不够直观，且仅从大量的数据表现方式上难以快速找出重点污染区。采用有限的大气监测站点进行空间分析，直观地反映区域内的污染物分布信息，是公众关心的热点问题。

现有的成熟的空间分析方法主要包括：反距离插值权重（IDW）、克里金（Kriging）、样条函数（Spline）法。IDW法是基于Tobler定理提出的一种简单插值方法。其原理是通过计算以被估计点和其附近各个观测点的距离的p次幂表示的权重，进行加权平均来进行插值。Kriging法不仅考虑观测点和被估计点的相对位置，而且还考虑了各观测点之间的相对位置关系，该法主要分为两步：第一步是对空间进行结构分析，选取变异函数模型并计算模型参数；第二步是在变异函数模型的基础上进行插值计算。样条函数法（Spline）是采用分段多项式逼近已知数据点，同时保证在各段交接处有一定的光滑性，一般要求二阶光滑。

由于固定式大气环境监测站点数量有限，通常分布不均匀，同时污染物浓度分布受污染物物理化性质、地表面特征、气象及人为等因素的影响，并且插值区域的尺度影响污染物分布的空间自相关性，采用IDW法，插值速度快，但只考虑距离因素，可能出现插值精度不高的问题；采用Kriging法，考虑了空间自相关性，但计算量大，插值速度慢，缺乏空间自相关性时插值效果亦不佳；会导致插值结果不够及时、准确和有效地得到。

发明内容

基于此，有必要针对现有的区域空气质量空间插值精度不高、浓度插值结果不可靠的技术问题，提供一种区域空间质量分析方法。

一种区域空间质量分析方法，包括：

对所述区域进行网格划分，并确定每个网格中心点的经纬度坐标；

获取单个时间片上的大气环境监测数据，提取至少一个大气环境固定监测站点的经纬度坐标和污染物浓度监测值；

对每个网格中心点，根据所述监测站点的经纬度坐标和污染物浓度监测值确定空间分布特征，基于所述空间分布特征确定以每个网格中心点为插值点的影响点搜索范围；

以所述影响点搜索范围内的监测站点作为插值点的影响样本点，确定每个影响样本点的权重值；

对每个插值点采用对应的影响样本点的污染物浓度监测值，根据对应的权重值进行插值计算，得到与每个插值点对应的网格中心点的浓度插值结果。

进一步的，所述对每个网格中心点，根据所述监测站点的经纬度坐标和污染物浓度监测值确定空间分布特征，基于所述空间分布特征确定以每个网格中心点为插值点的影响点搜索范围的步骤，具体包括：

每两个监测站点组成监测站点对，确定监测站点对所呈的位置角度，确定拟合方向个数，并给定预设方向角度误差范围，根据所述监测站点对的位置角度信息，判断监测站点对所属的拟合方向；

计算同一拟合方向上的多个监测站点对的试验变异函数值，所述试验变异函数以监测站点对间的距离为自变量；

根据所述试验变异函数值以及对应的监测站点距离进行拟合，得到每个方向上的拟合函数，并获得每个拟合方向上的变程值；

根据所述每个拟合方向上的变程值，以插值点为原点，拟合为圆或椭圆，以所述圆或椭圆所包括的范围作为影响点搜索范围。

更进一步的，所述计算同一拟合方向上的多个监测站点对间的试验变异函数值，所述试验变异函数以监测站点对间的距离为自变量的步骤，具体包括：

确定每个拟合方向上的监测站点，在拟合方向上的角度的预设误差范围内的监测站点为该拟合方向上的监测站点；

所述试验变异函数为

其中，监测站点x_i和监测站点x_j构成所述监测站点对，Z(x_i)，Z(x_j)分别为监测站点x_i和监测站点x_j两点的污染物浓度监测值；E表示半方差；θ为监测站点x_i和监测站点x_j两点位置间的角度；h为监测站点x_i和监测站点x_j两点间的距离，且

其中，（x_i',y_i'）为监测站点x_i的经纬度坐标；（x_j',y_j'）为监测站点x_j的经纬度坐标。

更进一步的，所述根据所述试验变异函数值以及对应的监测站点距离进行拟合，得到每个方向上的拟合函数，计算每个拟合函数的变程值，得到每个拟合方向上的变程值的步骤，具体包括：

对各个方向的试验变异函数进行球面理论变异函数模型的拟合，得出每个方向对应的球面模型参数；

当所述拟合函数趋于水平变化时的空间距离作为该方向上的变程值；

所述根据所述每个拟合方向上的变程值，以其对应的网格中心点为原点，拟合为圆或椭圆，以所述圆或椭圆所包括的范围作为影响点搜索范围的步骤，具体包括：

当拟合为圆时，以各个方向的试验变异函数的平均变程值为半径，插值点为圆心，绘制的圆的内部区域即为影响点搜索范围；

当拟合为椭圆时，以各个方向的试验变异函数的变程值拟合出来的长半轴及短半轴、顺时针偏转角度，以插值点为轴心，绘制的椭圆内部区域即为影响点搜索范围。

进一步的，所述以所述影响点搜索范围内的监测站点作为插值点的影响样本点，确定每个影响样本点的权重值的步骤，具体包括：

确定至少一种搜索策略的搜索参数，所述搜索参数包括搜藏距离、搜索方向个数以及单个搜索方向上的最小和最大搜索点位个数；

对每一种搜索策略，确定每个搜索方向上的影响点位数，根据对应的搜索参数，选择每个方向上满足搜索条件的监测站点作为插值点的影响样本点；

确定至少一种权重值设定方案，对每个影响样本点根据权重值设定方案计算对应的权重值。

更进一步的，所述搜索条件为：

若搜索方向上影响点位数k小于最小搜索点位数c，且搜索方向上存在的监测站点数大于0，则与插值点最近的c个监测站点满足搜索条件；

若搜索方向上影响点位数k小于最小搜索点位数c，且搜索方向上存在的监测站点数为0，则搜索方向上的所有监测站点均不满足搜索条件；

若搜索方向上最小搜索点位数c<影响点位数k<最大搜索点位数d，则与插值点最近的k个监测站点满足搜索条件；

若搜索方向上影响点位数k>最大搜索点位数d，则与插值点最近的d个监测站点满足搜索条件。

更进一步的，所述确定至少一种权重值设定方案，对每个影响样本点根据权重值设定方案计算对应的权重值的步骤，具体包括：

设定至少一个幂值p，每个权重值设定方案采用一个幂值；

计算

其中λ_i为影响点搜索范围的第i个影响样本点的权重，n为影响点搜索范围的影响样本点总数；d_i为插值点与第i个影响点搜索范围的影响样本点之间的距离。

进一步的，所述对每个插值点采用对应的影响样本点的污染物浓度监测值，根据对应的权重值进行插值计算，得到与每个插值点对应的网格中心点的浓度插值结果的步骤，具体包括：

对每个插值点采用对应的影响样本点的污染物浓度监测值，计算在至少一个搜索策略和权重值设定方案下的污染物浓度属性值

其中，λ_i为影响点搜索范围的第i个影响样本点的权重，Z(x_o)为插值点的污染物浓度属性值。

再进一步的，所述对每个插值点采用对应的影响样本点的污染物浓度监测值，根据对应的权重值进行插值计算，得到与每个插值点对应的网格中心点的浓度插值结果的步骤，还包括：

采用交叉验证法，以均方根误差为评价指标，对获得的至少一个搜索策略和权重值设定方案下的污染物浓度属性值进行验证，选择均方根误差最小的搜索策略和权重值设定方案下的污染物浓度属性值作为与插值点对应的网格中心点的浓度插值结果。

进一步的，还包括：

根据每个网格中心点的浓度插值结果确定所对应的颜色特征；

对每个网格赋予其网格中心点的浓度插值结果所对应的颜色特征；

根据每个网格的颜色特征绘制渲染图形。

本发明的区域空间质量分析方法，适用于中小尺度的区域空气质量空间分析，该方法综合考虑了IDW算法和Kriging的优缺点，提出了改进的IDW算法，同时保持了IDW算法插值速度快和Kriging算法考虑空间结构的插值精度高的优点，在不需要求解方程获得样本点权重系数的基础上，从搜索策略上优化了空间插值过程，使得保证一定计算速度的同时，插值精度更高，为公众提供更为准确的、直观的区域污染物空间分布信息。

附图说明

图1为本发明一种区域空间质量分析方法的工作流程图；

图2为本发明一个搜索范围的拟合圆的例子；

图3为本发明一个搜索范围的拟合椭圆的例子；

图4为本发明一个例子的工作流程图。

具体实施方式

下面结合附图和具体实施例对本发明做进一步详细的说明。

如图1所述为本发明一种区域空间质量分析方法的工作流程图，包括：

步骤S101，对所述区域进行网格划分，并确定每个网格中心点的经纬度坐标；

步骤S102，获取单个时间片上的大气环境监测数据，提取至少一个大气环境固定监测站点的经纬度坐标和污染物浓度监测值；

步骤S103，对每个网格中心点，根据所述监测站点的经纬度坐标和污染物浓度监测值确定空间分布特征，基于所述空间分布特征确定以每个网格中心点为插值点的影响点搜索范围；

步骤S104，以所述影响点搜索范围内的监测站点作为插值点的影响样本点，确定每个影响样本点的权重值；

步骤S105，对每个插值点采用对应的影响样本点的污染物浓度监测值，根据对应的权重值进行插值计算，得到与每个插值点对应的网格中心点的浓度插值结果。

其中，在步骤S101中，对所述区域进行网格划分可以采用如下方法：

步骤S1011，对于指定的研究区域，绘制包含整个区域的矩形范围；

步骤S1012，选取适当的网格精度，对该矩形范围进行进一步网格划分，网格中的格子数目为M*N；

步骤S1013，根据区域尺度类型选取平面或墨卡托坐标系，确定各个网格中心点的经纬度坐标（x_k，y_j），k=0,1...,M，j=0,1...,N。

在其中一个实施例中，步骤S102，具体包括：

选取某一时刻，读取该时刻的大气环境监测基础数据；

提取用于做进一步空间分析的有效数据，包括各个大气环境监测站点的经纬度坐标信息，以及相应站点上的污染物浓度信息，污染物浓度信息包括污染物浓度监测值；此处的经纬度坐标系须与步骤S1033一致。

在其中一个实施例中，步骤S103，具体包括：

步骤S1031，每两个监测站点组成监测站点对，确定监测站点对所呈的位置角度，确定拟合方向个数，并给定预设方向角度误差范围，根据所述监测站点对的位置角度信息，判断监测站点对所属的拟合方向；

步骤S1032，计算同一拟合方向上的多个监测站点之间的试验变异函数值，所述试验变异函数以监测站点之间的距离为自变量；

步骤S1033，根据所述试验变异函数值以及对应的监测站点距离进行拟合，得到每个方向上的拟合函数，计算每个拟合函数的变程值，得到每个拟合方向上的变程值；

步骤S1034，根据所述每个拟合方向上的变程值，以插值点为原点，拟合为圆或椭圆，以所述圆或椭圆所包括的范围作为影响点搜索范围。

在其中一个实施例中，步骤S1031具体包括：将所有的监测站点按照空间位置，沿着x轴从左到右排列编号，如1～N；从1开始，计算监测站点1与监测站点2,3,4,5...N所呈的角度。例如：设监测站点A为基准向量（100,0），计算监测站点B与C所形成的向量BC，记为D(x′_D,y′_D)，则监测站点对的位置角度为向量A与D所呈的角度P：

若y′_D>0，则有

$P=\arccos \left(\frac{A*D}{\left|A\right|*\left|D\right|}\right)*\frac{180}{\mathrm{\&pi;}}=\arccos \left(\frac{x_D^{\&prime;}}{\sqrt{{\left(x_D^{\&prime;}\right)}^2+{\left(y_D^{\&prime;}\right)}^2}}\right)*\frac{180}{\mathrm{\&pi;}}$

若y′_D＜0，则有

$P=\arccos \left(\frac{A*D}{\left|A\right|*\left|D\right|}\right)*\frac{180}{\mathrm{\&pi;}}+180=\arccos \left(\frac{x_D^{\&prime;}}{\sqrt{{\left(x_D^{\&prime;}\right)}^2+{\left(y_D^{\&prime;}\right)}^2}}\right)*\frac{180}{\mathrm{\&pi;}}+180;$

遍历所有点（如N个点），得到所有监测站点对的角度信息。

在其中一个实施例中，步骤S1031前，还可以包括：对已知的大气环境监测数据进行探索性数据分析，主要包括探索数据是否满足正态分布、剔除过大过小异常值、初步寻求及研究数据的空间自相关性。

在其中一个实施例中，步骤S1032，具体包括：

确定每个拟合方向上的监测站点，在拟合方向上的角度的预设误差范围内的监测站点为该拟合方向上的监测站点；

所述试验变异函数为

（公式1），其中，Z(x_i)，Z(x_j)分别为监测站点x_i和监测站点x_j两点的污染物浓度监测值；E表示半方差；h为监测站点x_i和监测站点x_j两点间的距离，且

其中，（x_i',y_i'）为已知影响样本点x_i的经纬度坐标；（x_j',y_j'）为插值点x_j的经纬度坐标。

因为监测站点监测到某固定位置的污染物浓度值是随时间不断变化的，可认为污染物浓度监测事件为随机事件，而期望可以衡量和距离有关的污染物浓度变异程度。可以认为某区域内的不同的站点监测事件对应的区域随机变量服从独立同一分布（iid），另外由于总体分布是未知的，可以通过样本均值来计算。如下所示：

${\mathrm{\&gamma;}}^{*}\left(h\right)={\mathrm{\&gamma;}}^{*}\left(\right|x_i-x_j\left|\right)=\frac{1}{2m}{\mathrm{\&Sigma;}}_{k=1}^m{[Z\left(x_i\right)-Z\left(x_j\right)]}^2---\left(2\right)$

m表示实际计算时距离为h的样本点对(x_i,x_j)的个数。

另外污染物浓度由于受当地的风速、气温等具有地理特征因素的影响，污染物浓度变异程度还和已知样点的地理位置有关。表示如下：

$\mathrm{\&gamma;}(h,\mathrm{\&theta;})=\mathrm{\&gamma;}\left(\right|x_i-x_j|,\mathrm{\&theta;})=\frac{1}{2}E{[Z\left(x_i\right)-Z\left(x_j\right)|\mathrm{\&theta;}]}^2---\left(3\right)$

Θ为样本点x_i和x_j的位置角度。

公式（3）的意思是，若是任意两个站点的空间距离和位置角度相同，污染物浓度的变异值相同。对应的样本均值表达式为

${\mathrm{\&gamma;}}^{*}(h,\mathrm{\&theta;})={\mathrm{\&gamma;}}^{*}\left(\right|x_i-x_j|,\mathrm{\&theta;})=\frac{1}{2m}{\mathrm{\&Sigma;}}_{k=1}^m{[Z\left(x_i\right)-Z\left(x_j\right)|\mathrm{\&theta;}]}^2---\left(4\right)$

设定方向特征的个数N（N=4,8），每一个方向对应的方向角度为

其中p=1，…,N；确定容忍角

并计算每个方向组的角度范围

另外N=0表示无方向，将不同的角度范围代入公式（3）或（4）的θ中，得到每个拟合方向上的试验变异函数。

在其中一个实施例中，步骤S1033，具体包括：

所述根据所述试验变异函数值以及对应的监测站点距离进行拟合，得到每个方向上的拟合函数，计算每个拟合函数的变程值，得到每个拟合方向上的变程值的步骤，具体包括：

对各个方向的试验变异函数进行球面理论变异函数模型的拟合，得出每个方向对应的球面模型参数；

当所述拟合函数趋于水平变化时的空间距离作为该方向上的变程值。

具体地，选取常用的拟合效果较好的球状模型对无方向的变异函数散点进行拟合，球面模型的表达式分别为：

$\mathrm{\&gamma;}\left(h\right)=\left\{\begin{array}{cc}0& \\ C_o+C(\frac{3h}{2a}-\frac{h^3}{{2a}^3})& (0<h\&le;a)\\ C_o+C& (h>a)\end{array}\right.---\left(5\right)$

进而对各个方向的变异函数散点图进行球面理论变异函数模型的拟合，得出每个方向对应的球面模型参数。

根据每个方向对应的球面模型不同变程值，绘制玫瑰图，根据玫瑰图判定污染物的空间分布为各向同性表达还是各向异性表达。若是玫瑰图呈现圆的特性，则认为空间分布具有各向同性的特点，如图2所示；若具有椭圆的特性，则为各向异性表达。实例如图3所示。

步骤S1034，具体包括：

当拟合为圆时，以各个方向的试验变异函数的平均变程值为半径，插值点为圆心，绘制的圆的内部区域即为影响点搜索范围；

根据其描述，以插值点x₀=(x′₀,y′₀)为圆心，平均变程值为半径，得到的圆的表达式为：

$\frac{{(x-x_0^{\&prime;})}^2}{{\stackrel{\&OverBar;}{h}}^2}+\frac{{(y-y_0^{\&prime;})}^2}{{\stackrel{\&OverBar;}{h}}^2}=1.$

当拟合为椭圆时，以各个方向的试验变异函数的变程值拟合出来的长半轴a及短半轴b、顺时针偏转角度θ，以插值点为轴心，绘制的椭圆内部区域即为影响点搜索范围。

根据其描述，以插值点x₀=(x′₀,y′₀)为轴心，得到的椭圆表达式为

$\frac{{(x\cos \mathrm{\&theta;}+y\sin \mathrm{\&theta;}-x_0^{\&prime;})}^2}{a^2}+\frac{{(y\cos \mathrm{\&theta;}-x\sin \mathrm{\&theta;}-y_0^{\&prime;})}^2}{b^2}=1.$

在其中一个实施例中，步骤S104，具体包括：

步骤S1041，确定至少一种搜索策略的搜索参数，所述搜索参数包括搜索距离、搜索方向个数以及单个搜索方向上的最小和最大搜索点位个数；

步骤S1042，对每一种搜索策略，确定每个搜索方向上的影响点位数，根据对应的搜索参数，选择每个方向上满足搜索条件的监测站点作为插值点的影响样本点；

步骤S1043，确定至少一种权重值设定方案，对每个影响样本点根据权重值设定方案计算对应的权重值。

在其中一个实施例中，所述搜索条件为：

若搜索方向上影响点位数k小于最小搜索点位数c，且搜索方向上存在的监测站点数大于0，则与插值点最近的c个监测站点满足搜索条件；

若搜索方向上影响点位数k小于最小搜索点位数c，且搜索方向上存在的监测站点数为0，则搜索方向上的所有监测站点均不满足搜索条件；

若搜索方向上最小搜索点位数c<影响点位数k<最大搜索点位数d，则与插值点最近的k个监测站点满足搜索条件；

若搜索方向上影响点位数k>最大搜索点位数d，则与插值点最近的d个监测站点满足搜索条件。

在其中一个实施例中，所述步骤S1043，具体包括：

设定至少一个幂值p，每个权重值设定方案采用一个幂值；

计算

（6），其中λ_i为影响点搜索范围的第i个影响样本点的权重，n为影响点搜索范围的影响样本点总数；d_i为插值点与第i个影响点搜索范围的影响样本点之间的距离。

在其中一个实施例中，所述步骤S105，具体包括：

对每个插值点采用对应的影响样本点的污染物浓度监测值，计算在至少一个搜索策略和权重值设定方案下的污染物浓度属性值（7），其中，λ_i为影响点搜索范围的第i个影响样本点的权重，Z(x_o)为插值点的污染物浓度属性值。

在其中一个实施例中，步骤S105，还包括：

采用交叉验证法，以均方根误差为评价指标，对获得的至少一个搜索策略和权重值设定方案下的污染物浓度属性值进行验证，选择均方根误差最小的搜索策略和权重值设定方案下的污染物浓度属性值作为与插值点对应的网格中心点的浓度插值结果。

交叉验证的基本思路为：假设存在A个影响样本点，假设第i（i=1,2,3...A）个样本点不存在，用剩余的A-1个影响样本点的污染物浓度属性值对该点进行插值，依次计算出这A个站点的污染物浓度属性估计值，通过计算A个影响样本点用于插值获得的均方根误差，作为评价插值精度。均方根误差的具体计算公式的具体实现过程如下：

(8)其中，Z_a,i为第i个影响样本点的污染物浓度属性值，Z_e,i为第i个影响样本点的污染物浓度属性估计值，m为验证影响样本点总数。

选择均方根误差最小的搜索策略和权重值设定方案下，即确定均方根误差最小时的幂值p、搜索方向个数、各方向最大、最小搜索点位数等参数。

在其中一个实施例中，本发明一种区域空间质量分析方法还包括：

根据每个网格中心点的浓度插值结果确定所对应的颜色特征；

对每个网格赋予其网格中心点的浓度插值结果所对应的颜色特征；

根据每个网格的颜色特征绘制渲染图形，优选采用地理信息系统（GIS）技术绘制。

作为一个例子，如图4所示，包含以下步骤：

S1：根据珠三角区域绘制矩形边框，并选择250m*250m的网格精度，对矩形进行网格划分，格子数目为536*536个网格，并获取各个网格中心坐标（x_k，y_j），k=0,1...,536，j=0,1...,536；

S2：以2012年10月1日珠三角区域62个监测站点的SO₂与PM_2.5监测数据为例，获取的主要信息包括62个监测站点的经纬度坐标、污染物浓度，根据纬度坐标由小到大对监测站点编号。

S3：分析珠三角区域62个监测站点SO₂与PM_2.5的空间自相关性特征，根据空间自相关性特征设定空间影响距离（即变程值），拟合出的椭圆搜索范围如图3所示，同时设置IDW算法的参数，包括幂值、搜索方向个数、各个方向上的最大最小点位个数。空间自相关性的计算公式具体如下：

$\mathrm{\&gamma;}=(h,\mathrm{\&theta;})=\frac{1}{N\left(h\right)}\underset{i=1}{\overset{N\left(h\right)}{\mathrm{\&Sigma;}}}{(Z\left({\mathrm{\&mu;}}_i\right)-Z({\mathrm{\&mu;}}_i+h))}^2$

其中，Z(μ_i)，Z(μ_i+h)分别为监测样点μ_i，μ_i+h的污染物浓度监测值；E表示半方差；θ为监测站点μ_i，μ_i+h两点所成的角度；h为监测站点μ_i，μ_i+h两点间的距离：

$h=\frac{1}{2}{({(x_i^{\&prime;}-x_j^{\&prime;})}^2+{(y_i^{\&prime;}-y_j^{\&prime;})}^2)}^{h={({x_i}^{\&prime;}-{x_j}^{\&prime;})}^2+{({y_i}^{\&prime;}-{y_j}^{\&prime;})}^2}$

其中，（x_i',y_i'）为已知样点μ_i的经纬度坐标；（x_j',y_j'）为监测站点μ_i+h的经纬度坐标；

S4：根据S3中上述参数设定插值过程中的各项参数，包括方向个数（如八方向）、变程值（43.2～68km）、最大、最小点位个数，并计算出各个影响点位的权重系数，具体过程如下：

${\mathrm{\&lambda;}}_i=\frac{{\left(d_i\right)}^{-p}}{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\left(d_i\right)}^{-p}}$

其中i为影响插值点的样本点；n为影响插值点的样本点总个数。

S5：根据影响点位的浓度属性及权重系数，计算536*536个网格中心点的污染物浓度属性，并根据污染物浓度范围赋予不同颜色，采用ArcGIS平台制作渲染图。插值点污染物浓度具体计算公式如下：

$Z\left(x_o\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_{i}Z\left(x_i\right)=\frac{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\frac{1}{{\left(d_i\right)}^p}Z\left(x_i\right)}{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\frac{1}{{\left(d_i\right)}^p}}$

其中，λ_i为已知样点x_i的权重系数；Z(x_o)为插值点的污染物浓度属性值；p为幂值；i表示影响插值点的第i个样点，n分别为影响插值点的样点总数；d_i为插值点x_o与已知样点x_i之间的距离。

采用交叉验证法，以均方根误差（RMSIE）为评价指标，对插值精度进行验证，并确定精度最高时的各项参数。交叉验证的基本思路为：假设存在A个影响样本点，假设第i（i=1,2,3...A）个样本点不存在，用剩余的A-1个影响样本点的属性值对该点进行插值，依次计算出这A个站点的属性估计值，通过计算A个站点的均方根误差，用以评价插值精度。误差（ME）与均方根误差（RMSIE）的具体计算公式的具体实现过程如下：

$\mathrm{RMSIE}=\sqrt{\frac{\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}{(Z_{a,i}-Z_{e,i})}^2}{m}}$

其中，Z_a,i为第i个站点的实际监测浓度值，Z_e,i为第i个站点的预测值，m为验证站点总数。

以PM_2.5为例说明，当搜索范围采用上述椭圆模式，分别设置长轴100km、短轴80km及偏转角度45°，幂值为1，搜索方向为四方向、最大点位数为6、最小点位数为2时，ME与RMSIE值分别为-0.19及11.14（改进前取p=2，无方向，ME与RMSIE分别为0.25及11.38）。上述基于插值优化的区域空气质量空间分析方法在保证简单计算量的同时，有效提高了插值精度。

采用适当的配色方案，结合GIS技术，对上述网格数据进行渲染，获得渲染图。

以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对本发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。

Claims (10)

1.一种区域空间质量分析方法，其特征在于，包括：

对所述区域进行网格划分，并确定每个网格中心点的经纬度坐标；

获取单个时间片上的大气环境监测数据，提取至少一个大气环境固定监测站点的经纬度坐标和污染物浓度监测值；

对每个网格中心点，根据所述监测站点的经纬度坐标和污染物浓度监测值确定空间分布特征，基于所述空间分布特征确定以每个网格中心点为插值点的影响点搜索范围；

以所述影响点搜索范围内的监测站点作为插值点的影响样本点，确定每个影响样本点的权重值；

对每个插值点采用对应的影响样本点的污染物浓度监测值，根据对应的权重值进行插值计算，得到与每个插值点的浓度插值结果。

2.根据权利要求1所述的区域空间质量分析方法，其特征在于，所述对每个网格中心点，根据所述监测站点的经纬度坐标和污染物浓度监测值确定空间分布特征，基于所述空间分布特征确定以每个网格中心点为插值点的影响点搜索范围的步骤，具体包括：

每两个监测站点组成监测站点对，确定监测站点对所呈的位置角度，确定拟合方向个数，并给定预设方向角度误差范围，根据所述监测站点对的位置角度信息，判断监测站点对所属的拟合方向；

计算同一拟合方向上的多个监测站点对间的试验变异函数值，所述试验变异函数以监测站点对间的距离为自变量；

根据所述试验变异函数值以及对应的监测站点距离进行拟合，得到每个方向上的拟合函数，并进一步获得每个拟合方向上的变程值；

根据所述每个拟合方向上的变程值，以其对应的网格中心点为原点，拟合为圆或椭圆，以所述圆或椭圆所包括的范围作为影响点搜索范围。

3.根据权利要求2所述的区域空间质量分析方法，其特征在于，所述计算同一拟合方向上的多个监测站点对间的试验变异函数值，所述试验变异函数以监测站点对间的距离为自变量的步骤，具体包括：

所述试验变异函数为

其中，监测站点x_i和监测站点x_j构成所述监测站点对，Z(x_i)，Z(x_j)分别为监测站点x_i和监测站点x_j两点的污染物浓度监测值；E表示期望值；θ为监测站点x_i和监测站点x_j的位置角度，h为监测站点x_i和监测站点x_j两点间的距离，且

其中，（x_i',y_i'）为监测站点x_i的经纬度坐标；（x_j',y_j'）为监测站点x_j的经纬度坐标。

4.根据权利要求2所述的区域空间质量分析方法，其特征在于，所述根据所述试验变异函数值以及对应的监测站点距离进行拟合，得到每个方向上的拟合函数，计算每个拟合函数的变程值，得到每个拟合方向上的变程值的步骤，具体包括：

对各个方向的试验变异函数进行球面理论变异函数模型的拟合，得出每个方向对应的球面模型参数；

当所述拟合函数趋于水平变化时的空间距离作为该方向上的变程值；

所述根据所述每个拟合方向上的变程值，以其对应的网格中心点为原点，拟合为圆或椭圆，以所述圆或椭圆所包括的范围作为影响点搜索范围的步骤，具体包括：

当拟合为圆时，以各个方向的试验变异函数的平均变程值为半径，插值点为圆心，绘制的圆的内部区域即为影响点搜索范围；

当拟合为椭圆时，以各个方向的试验变异函数的变程值拟合出来的长半轴及短半轴、顺时针偏转角度，以插值点为轴心，绘制的椭圆内部区域即为影响点搜索范围。

5.根据权利要求1所述的区域空间质量分析方法，其特征在于，所述以影响点搜索范围内的监测站点作为插值点的影响样本点，确定每个影响样本点的权重值的步骤，具体包括：

确定搜索策略中的搜索参数，所述搜索参数包括搜索距离、搜索方向个数以及单个搜索方向上的最小和最大搜索点位个数；

根据对应的搜索参数，选择每个方向上满足搜索条件的监测站点作为插值点的影响样本点；

确定至少一种权重值设定方案，对每个影响样本点根据权重值设定方案计算对应的权重值。

6.根据权利要求5所述的区域空间质量分析方法，其特征在于，所述搜索条件为：

若搜索方向上影响点位数k小于最小搜索点位数c，且搜索方向上存在的监测站点数大于0，则与插值点最近的c个监测站点满足搜索条件；

若搜索方向上影响点位数k小于最小搜索点位数c，且搜索方向上存在的监测站点数为0，则搜索方向上的所有监测站点均不满足搜索条件；

若搜索方向上最小搜索点位数c<影响点位数k<最大搜索点位数d，则与插值点最近的k个监测站点满足搜索条件；

若搜索方向上影响点位数k>最大搜索点位数d，则与插值点最近的d个监测站点满足搜索条件。

7.根据权利要求5所述的区域空间质量分析方法，其特征在于，所述确定至少一种权重值设定方案，对每个影响样本点根据权重值设定方案计算对应的权重值的步骤，具体包括：

设定至少一个幂值p，每个权重值设定方案采用一个幂值；

计算

其中λ_i为影响点搜索范围的第i个影响样本点的权重，n为影响点搜索范围的影响样本点总数；d_i为插值点与第i个影响点搜索范围的影响样本点之间的距离。

8.根据权利要求5所述的区域空间质量分析方法，其特征在于，所述对每个插值点采用对应的影响样本点的污染物浓度监测值，根据对应的权重值进行插值计算，得到与每个插值点对应的网格中心点的浓度插值结果的步骤，具体包括：

对每个插值点采用对应的影响样本点的污染物浓度监测值，计算在至少一个搜索策略和权重值设定方案下的污染物浓度属性值

其中，λ_i为影响点搜索范围的第i个影响样本点的权重，Z(x_o)为插值点的污染物浓度属性值。

9.根据权利要求8所述的区域空间质量分析方法，其特征在于，所述对每个插值点采用对应的影响样本点的污染物浓度监测值，根据对应的权重值进行插值计算，得到与每个插值点对应的网格中心点的浓度插值结果的步骤，还包括：

采用交叉验证法，以均方根误差为评价指标，对获得的至少一个搜索策略和权重值设定方案下的污染物浓度属性值进行验证，选择均方根误差最小的搜索策略和权重值设定方案下的污染物浓度属性值作为与插值点对应的网格中心点的浓度插值结果。

10.根据权利要求1所述的区域空间质量分析方法，其特征在于，还包括：

根据每个网格中心点的浓度插值结果确定所对应的颜色特征；

对每个网格赋予其网格中心点的浓度插值结果所对应的颜色特征；

根据每个网格的颜色特征绘制渲染图形。

Images