CN107610021A - 环境变量时空分布的综合分析方法 - Google Patents

Abstract

本发明公开了一种环境变量时空分布的综合分析方法，首先计算各时空滞后距上的试验变异函数值，并拟合理论变异函数模型；然后联合时空样点数据，进行时空克里格插值，估算未测时空位置的地理属性值；并建立环境变量值与区域位置的定量关系，形成环境变量值时空分布的趋势体，得到基于时空趋势克里格预测结果；提出多种时空不确定性评估方法；最后直观地对环境变量的时空分布状况进行评价。本发明充分考虑了环境变量时空结构性和连续性的特点，从建模、预测、趋势分析、不确定性分析、时空分析多个角度对环境变量时空分布特点进行了充分的模拟和分析，为区域环境评估和相关部门进行时空决策及辅助分析提供了方法依据。

Description

环境变量时空分布的综合分析方法

技术领域

本发明涉及资源生态环境建模与评估技术领域，具体涉及一种环境变量时空分布的综合分析方法。

背景技术

资源生态环境属性(如土壤水盐、重金属含量，大气颗粒物浓度等)普遍存在时空相关性，表现为在空间上分布的结构性，以及随时间变化的连续性。目前，随着监测手段和设备的发展，使得对区域资源环境属性的连续监测和测量成为可能，能够获取多时期的环境变量数据。基于这些数据，获取被研究对象的时空分布特征，挖掘其中蕴含的更多信息，是探求其时空演变规律，进而揭示其中蕴含的驱动机制，是资源环境评估和治理的基础，能够为区域范围内预防环境污染、出台治理措施提供科学依据。

基于点状量测数据，地统计学方法是资源环境属性空间建模和预测的常用方法。其中，空间变异模型建立了距离与环境变量变异的定量关系，表征其空间结构性，在此基础上，多种克里格方法可用于对未测位置环境变量值得预测，最终形成区域资源环境属性含量或不确定性的空间分布图。由于较高的预测精度和成熟的软件支持，近几十年，地统计方法被广泛应用于资源生态环境的多个领域。与此同时，为应对实际应用中存在的问题或改善预测效果，地统计方法本身也在不断地发展壮大，产生了一系列的衍生方法和技术。其中，作为地统计方法的重要分支，时空地统计方法越来越得到重视，其应用领域也越来越广，为资源环境属性的时空分析提供了有效工具和新的思路。

然而，现阶段时空地统计方法及其应用存在诸多问题，使其在解决实际问题和推广中存在障碍，主要表现在以下几个方面：

(1)时空地统计方法虽然早已提出，但具体计算步骤仍不明确，如理论模型参数拟合、时空邻近点确定等。这阻碍了方法的应用及推广。

(2)缺乏对时空分布趋势的考虑：目前，大部分研究和应用是基于时空离散量测点进行时空建模及预测。但大部分的地理属性存在时空趋势性，而在建模和预测时剥离这种趋势能够有效地提高时空预测精度。

(3)缺乏时空不确定性分析方法：在空间地统计中，不确定性分析与制图能够获取研究对象超过或不超过某阈值的概率空间分布图，是一些领域(如环境污染物阈值预警)关注的重点。由于只有一个时期数据，空间维度制图只能获取采样时期的不确定性空间分布图。而在多时期数据的支持下，不确定性分析能够明确一段时期内，地理属性均超过或不超过某阈值的概率分布。这对于分析长时间跨度环境污染物的严重程度，制定相应的治理和应对措施有重要参考意义。而目前国内外尚缺乏这方面的研究。

(4)需要更多的时空分析方法和定量化指标：基于时空插值获得的时空立方体数据，尚缺乏足够的分析方法和定量化指标，不能挖掘出更多的有效信息。

发明内容

本发明的目的在于提供一种环境变量时空分布的综合分析方法，该方法充分考虑了环境变量时空结构性和连续性的特点，从建模、预测、趋势分析、不确定性分析、时空分析多个角度对环境变量时空分布特点进行了充分的模拟和分析，为区域环境评估和相关部门进行时空决策及辅助分析提供了方法依据。

为解决上述技术问题，本发明所设计的环境变量时空分布的综合分析方法，包括如下步骤：

步骤S1：基于环境变量时空样点，计算各个时空滞后距上的试验变异函数值，并拟合得到理论时空变异函数模型；

步骤S2：根据理论时空变异函数模型并采用所述环境变量时空样点数据，进行时空克里格插值，并根据时空克里格插值的结果来估算未测时空位置的环境变量值，所述时空克里格插值的结果为研究范围内的环境变量时空立方体数据；

步骤S3：依据环境变量与时空位置的关系，以时空位置为自变量，环境变量为因变量，建立环境变量值与区域位置的定量关系，形成环境变量值时空分布趋势体，再对环境变量时空样点剥离趋势后的残差进行时空克里格插值，插值结果与环境变量值时空分布趋势体的趋势值相加，得到基于时空趋势克里格估值的预测结果；

步骤S4：采用时空序贯指示模拟方法对环境变量进行模拟得到时空序贯指示模拟结果，并根据时空不确定性评估方法，对环境变量时空分布进行不确定性评估，得到地理属性超过或不超过某阈值的概率分布；

步骤S5：基于时空克里格插值获得的研究范围内的环境变量时空立方体数据，采用时空立方体数据的信息挖掘方法对环境变量的时空分布状况进行评价，获取时空分布规律。

所述步骤S1中的理论时空变异函数模型采用理论时空变异函数的模型库中的函数模型。

所述理论时空变异函数的模型库中的函数模型是根据试验变异函数散点图形状特征进行选择。

所述函数模型参数采用最小二乘法或遗传算法估计进行确定。

所述步骤2中，估算未测时空位置的环境变量值的具体方法为：

S21：构建时空插值网格，即确定时空克里格插值的空间和时间范围，以及时空克里格插值的时空粒度，所述空克里格插值的时空力度包括时空网格的空间边长及时间间隔；

S22：遍历每个时空插值网格，以时空插值网格中心为待估时空位置，搜索周围若干已知时空样点作为插值邻近点，构建克里格矩阵，解克里格矩阵，得到未测时空位置的环境变量值的时空预测结果。

所述步骤S3中得到基于时空趋势克里格估值的预测结果的具体方法为：

步骤S31：将时空分布趋势体表示为基于时空位置的方程；

步骤S32：对每个环境变量时空样点数据，减去时空分布趋势体，得到时空样点残差；利用时空克里格插值对时空样点残差进行时空插值，预测未采样时空位置的环境变量值，得到环境变量时空立方体数据；

步骤S33：对每个环境变量时空样点数据的时空分布趋势体与时空样点残差的插值结果进行相加，得到时空趋势克里格值的预测结果。

所述S4中采用时空序贯指示模拟方法对环境变量进行模拟得到时空序贯指示模拟结果的具体步骤为：

步骤S41：将环境变量时空样点数据转换成K个指示值；

步骤S42：对转换后的K个指示值，分别计算试验变异函数值，并拟合理论变异函数模型；

步骤S43：定义时空网格矩阵，采用时空序贯指示算法，对环境变量进行时空模拟；

步骤S44：设置模拟次数进行模拟，得到模拟结果。

所述步骤4中，对环境变量时空分布进行不确定性评估，得到地理属性超过或不超过某阈值的概率分布的具体范围为：

步骤S45：根据步骤S44的模拟结果，定义4种时空不确定性评估方案，所述4种时空不确定性评估方案包括单点单时刻不确定性方案、单点多时刻不确定性方案、多点单时刻不确定性方案和多点多时刻不确定性方案，得到地理属性超过或不超过某阈值的概率分布。

所述步骤S41中，将环境变量时空样点数据转换成K个指示值的具体方法为：

式中z_c为环境变量样点的若干几个截断值，I(s，t；z_c)为指示值，z(s,t)是在时空位置s,t处的环境变量值；

所述步骤S43中定义时空网格矩阵，采用时空序贯指示算法，对环境变量进行时空模拟的方法为；

步骤S43a：在定义的时空网格矩阵内随机生成模拟路线，路线按随机顺序经过所有时空网格；

步骤S43b：对在随机路线上的每一个时空网格节点的K个指示变量进行普通时空克里格插值，将得到该位置上环境变量值小于上述K个截断值的概率值，对K个概率值进行线性内插或外推，得到该位置上环境变量值的连续型累积概率密度分布函数；

步骤S43c：从上述累积概率密度函数中随机取1值作为该时空点上本次模拟的模拟值，将该模拟值按照公式4转换成指示值，并作为已测点，加入到计算随机路线上下一个时空节点的模拟计算中去；

步骤S43d：沿着随机路径，重复步骤S43b和步骤S43c，直到所有的时空网格节点都有模拟值，即完成一次模拟。

所述步骤5的具体方法为基于时空插值获得的环境变量时空立方体数据，根据查询条件和对象，将时空查询进行分类，根据分类模式和地理对象特征设置查询问题，得到时空分布规律。

本发明的有益效果：

本发明充分考虑了环境变量时空结构性和连续性的特点，从建模、预测、趋势分析、不确定性分析、时空分析多个角度对环境变量时空分布特点进行了充分的模拟和分析，为区域环境评估和相关部门进行时空决策及辅助分析提供了方法依据。从数据验证结果来看，时空预测方法精度要比只使用单一时期数据的空间插值方法精度高；精度高的趋势模型能带来更好的时空插值结果；从随机模拟结果的精度比较来看，时空序贯指示模拟方法可利用多时期数据对环境变量进行时空不确定性评估，并可取得空间随机模拟方法更高的模拟精度。而由于时间维的加入，能获得较空间不确定性评估更为丰富的结果和图件，如单点多时期不确定性和多点多时期不确定性，具体化为某时段内，环境变量超过或不超过某阈值概率的空间分布，这为长时期环境变量的时空分布及预警评估提供了方法；时空查询与分析能给用户挖掘更多的时空信息，使得时空数据更易于大众所理解。

本方法将细化最常用的时空普通克里格方法，明确使用步骤和计算公式，对环境变量的时空预测有积极意义；

本发明提出基于样点的时空趋势建模，及相应的时空回归克里格法，对提高时空预测精度有积极意义；

本发明扩展空间领域的不确定性分析方法，提出时空不确定性分析方法，一方面能够弥补这方面的空白，另一方面对区域长时间跨度的污染情况分析有重要意义

本发明提出一套针对时空立方体数据的分析方法和数值化的指标体系，能更直观地对地理属性的时空分布状况进行评价，有利于分析其时空分布规律，进而揭示其中蕴含的驱动机制。

附图说明

图1为本发明的环境变量时空分布的综合分析方法流程图；

图2为本发明的基于样点的武汉市青山区土壤重金属Cd的四种时空趋势模型；

图3为本发明中基于时空克里格和顾及趋势的时空克里格武汉市青山区土壤Cd时空分布图；

图4为本发明的4指示变量的时空试验变异散点(黑点)及理论时空变异模型拟合结果(曲面)图；

图5为本发明中普通时空克里格预测结果(a)与3次时空序贯指示模拟结果(b,c,d)对比图；

图6为山东省2014年第1(a)、100(b)、200(c)和300(d)天PM2.5浓度超过超过75μg·m^-3的概率分布；

图7为山东省2014年各月份内PM2.5浓度超过超过25μg·m^-3的概率分布；

图8为山东2014年的第1、100、200和300天PM2.5浓度超过75μg·m^-3概率大于0.8和0.9的空间分布；

图9为山东省2014年各月份内PM2.5浓度超过超过25μg·m^-3概率大于0.5空间分布。

图10为山东省在2014年内，PM2.5浓度超过给定阈值(a)75μg·m^-3、(b)115μg·m^-3、(c)150μg·m^-3、(d)250μg·m^-3的天数的空间分布；

图11为山东省2014年各季度和全年PM2.5浓度均值空间分布；

图12为依据PM2.5浓度的山东省2014年宜居性类别空间分布。

具体实施方式

以下结合附图和具体实施例对本发明作进一步的详细说明：

时空随机域的定义如下：定义时空随机域Z(p)＝{Z(s，t)|s∈S，t∈T}，其中S∈R²，表示空间域，通常有S＝{s₁，s₂}，s₁，s₂用于表示空间地理坐标；T∈R，表示时间域。举例来说，在克里格方法框架下，已知周围n个时空量测位置的镉含量值z(p_i)，p_i＝(s_i，t_i)，i＝1，...，n，如要预测土壤重金属镉在时空位置p₀＝(s₀，t₀)处的含量值z(p₀)。在这里，假设z(p_i)是时空随机域Z(p)的一个实现。

步骤S1的实现

在平稳假设(本征假设)条件下，时空地理属性残差部分R(p)＝R(s，t)的时空经验变异函数可表示成：

其中，γ表示变异函数值，h_S和h_T分别表示空间和时间滞后距，N(h_S，h_T)表示时空滞后距(h_S，h_T)条件下，配对的时空样点对数量。基于时空样点，z(s_i，t_i)和z(s_i+h_S，t_i+h_T)分别表示在时空位置(s_i，t_i)和(s_i+h_S，t_i+h_T)的环境变量值。获取经验变异函数后，需要拟合理论时空变异模型，部分理论时空变异模型见参考文献(Kolovos A.,Christakos G.,Hristopulos D.T.and Serre M.L.,2004.“Methods for generating non-separablespatiotemporal covariance models with potential environmental applications”.Advances in Water Resources 27:815-830)。

依据经验变异散点图，确定理论模型形式后，需要拟合模型中的参数值。对比于空间理论变异函数模型，时空理论变异函数模型更为复杂，有更多的参数需要估计，常规最小二乘法难以解决。可考虑遗传算法等智能搜索算法进行估计，思路如下：

(a)确定染色体编码及解码规则：

假设需要估计m(m≤2n+1)个参数，每个参数的最小值、最大值和估值精度分别是Umin_i，Umax_i和Q_i，则将m个参数分别以L₁,L₂,...,L_m为长度进行二进制编码，其中：

L_i＝ceil(log₂((Umax_i-Umin_i)/Q_i))， (公式2)

其中L_i表示第i个参数对应的染色体长度，ceil表示对括号内计算值向上取整。则每条染色体长度为染色体中第i个参数编码对应码公的解码公式为：

k表示第i个参数编码的二进制位数序号，b_k表示二进制第k位的值，以这种编码策略随机产生T组染色体。

(b)确定个体适应度评价函数

考虑到实际变异函数曲面上个点的重要性是不同的，往往滞后距小的点的重要性要大于滞后距大的点，因此采用以滞后距的倒数作为权系数参与适应度评价函数的构建。

(c)遗传操作：

确定合适的交叉概率和变异概率，依次进行选择、交叉和变异操作，从而产生新一代种群。设置最大进化次数为500或1000，完成后，选择最优染色体，进行解码，得到每个参数的估计值。

步骤S2的实现

构建时空插值网格，即确定插值的空间和时间范围，以及时空粒度(时空网格的空间边长及时间间隔)。遍历每个时空网格，以网格中心为待估时空位置，搜索周围若干已知时空样点作为插值邻近点，构建克里格矩阵：

矩阵中，γ表示变异函数值，p₀为待估时空位置，p_i(i＝1，2...n)为p₀周围n个已测时空样点的位置，λ_i为第i个样点所分配的权重系数。解此方程，求得λ_i(i＝1，2...n)，则p₀处地理属性的克里格估计值为：

式中，z(p_i)为在p_i处的样点属性值，λ_i为公式4计算得到的第i个样点所分配的权重值，z^*(p₀)为时空克里格法在p₀时空位置处的预测值。值得说明的是，在克里格估值中，需找到待估点附近的若干各邻近点构成克里格矩阵，这就要计算各采样点与待估点的距离。在空间克里格中，此距离等于欧式距离，但在时空克里格中，此距离应为时空距离，既要考虑两点间的空间距离，又要考虑两点间的时间距离。因此，如何确定时空距离的量算方式对临近点的选取十分重要。针对这个问题，提出两个解决方案：(a)如果时空变异模型中，包含了时空几何异向比率之类的参数，则将此参数直接使用，按照时空变异模型中的方式计算时空距离；(b)如果时空变异模型中没有相应参数，则定义时空距离计算公式为：

h_ST＝h_S+αh_T，其中

式中，h_S，h_T和h_ST分别表示时间距离、空间距离和时空距离，S_range为区域内空间跨度，即区域内最远两点的空间距离，T_range为数据的时间跨度，即所使用数据时间上的最长间隔，k为可调节参数，取1时，则时空数据三维形状将接近立方体。当k较小时，时空立方体被压缩，则为待估点搜索到相邻时期的样点作为邻近点的概率越大；反之，而当k较大时，时空立方体膨胀，则这种概率变小。因此，可根据单一时期所采样点密度调节k值，密度越大，则k值越大，反之，k值越小。α称为时空异向比例系数。

步骤S3的实现

定义时空随机域由以下两部分组成，即：

Z(p)＝m(p)+R(p) (公式7)

其中，Z(p)为在p位置处的地理属性值，m(p)为在p位置处的时空趋势值，R(p)为在p位置处去除趋势后的残差。在时空随机域框架下，如无其他相关辅助数据，那么时空趋势可表示为一系列的基于时空位置的方程，形式如下：

其中，μ和ν分别为空间域和时间域的最大次幂数，f_ρζ(s，t)表示μ×ν个已知基础函数项，而b_ρζ表示各已知基础函数项的系数，需要用已量测数据拟合，m(s，t)表达地理属性在时空域中的均值变化趋势。因此基础函数项又可表示为一系列的时空多项式

根据以往研究成果，μ＝ν＝2已足够表达大多数地理属性的实际时空趋势情况。因此，本发明中ρ≤μ＝2，ξ≤ν＝2。则根据ρ和ζ的不同选择和组合，时空趋势方程的阶数可有4中组合，分别是：

{ρ＝ζ＝1；ρ＝1，ξ＝2；ρ＝2，ξ＝1；ρ＝ξ＝2}。

样点剥离趋势后，对样点残差按照公式1计算残差的试验变异函数值，并按照S1中的步骤拟合理论变异函数模型，然后按照S2中的步骤对样点残差进行时空普通克里格插值。

步骤S4的实现

将原始时空样点数据转换成K个指示值，转换方法如下：

式中z_c为环境变量样点的若干几个截断值，一般取全部样点数据的20％、40％、60％和80％分位数对应的值，指示转换后，原始样点数据将转换为K组指示数据，I(s，t；z_c)为指示值，z(s,t)是在时空位置s,t处的环境变量值。对转换后的K组指示值，按照公式1计算残差的试验变异函数值，并按照S1中的步骤拟合理论变异函数模型。

定义时空网格矩阵，提出时空序贯指示算法，对环境变量进行时空模拟：

(a)在定义的时空网格矩阵内随机生成模拟路线，路线按随机顺序经过所有时空网格；

(b)对在随机路线上的每一个时空网格节点的K个指示变量进行普通时空克里格插值，将得到该位置上环境变量值小于上述K个截断值的概率值，对K个概率值进行线性内插或外推，得到该位置上环境变量值的连续型累积概率密度分布函数；

(c)从上述累积概率密度函数中随机取1值作为该时空点上本次模拟的模拟值，将该模拟值按照公式4转换成指示值，并作为已测点，加入到计算随机路线上下一个时空节点的模拟计算中去；

(d)沿着随机路径，重复步骤(b)和(c)，直到所有的时空网格节点都有模拟值，即完成一次模拟。

模拟次数设置为1000次或500次，基于这样的多次模拟结果，可以得到以下5种结果：

(1)每次模拟都可看作是一次时空插值的实现，与克里格方法相比，模拟的结果更注重对环境变量变异性的体现，在统计特征上更接近原始样点数据，而克里格方法结果更趋于平滑，追求的是最小预测方差；

(2)单点单时刻不确定性：设定阈值z_s，则时空不确定性，即在时空位置z(p′)＝z(s′，t′)上大于该阈值的概率，用P_ST表示，计算公式为：

式中，n(p′)为在1000次模拟实现中，时空位置p′处模拟结果大于该阈值的次数；

(3)单点多时刻不确定性：即计算每个空间位置上，多个时期(如一年或某几个月)均超过给定阈值的概率，用P_STT表示，计算公式为：

式中，t₁...t_q为模拟的q个时期，n_t(p′)为在空间位置p′处的1000次模拟中，q个时期模拟值均超过阈值的次数；

(4)多点单时刻不确定性：计算单一时刻，多个位置同时大于阈值的概率，用P_SST表示，计算公式为：

式中，n(p′₁，....，p′_m)为m个位置的1000次模拟实现中模拟值均大于阈值的次数。另外，可以给定一个判定概率p_c，可以划定超过一定污染风险概率的时空区域，公式为：

P_SST[z(p′)＞z_s]≥p_c(公式14)

(5)多点多时刻不确定性：计算多个位置在多个时期同时大于阈值的概率，用P_SSTT表示，计算公式为：

n(p′₁，....，p′_m)为在1000次模拟中，m个位置在q个时期，模拟值均大于阈值的次数，t₁-t_q表示从t₁时刻到t_q时刻，表示t₁时刻到t_q时刻括号内位置的环境变量属性值。同样的，可以通过定义判定概率p_ct来划定多个时期都超过一定污染风险概率的区域，公式为：

步骤S5的实现

(1)时空查询与分析

基于查询条件和对象，从根本上可将时空查询分为三类模式，根据这三类模式和研究领域的具体地理对象特征，派生出各种查询问题：

(a)What+Where→When，即由事件和地点作为条件，查询时间对象，如可查询在2014年内，PM2.5浓度超过给定阈值的天数的空间分布。

(b)When+Where→What，根据时间和地点查询事件或对象状态，如2014年全年和各季节山东省PM2.5的平均浓度和变异系数。

(c)When+What→Where，由事件或对象和时间查询地点，如按照一定规则划分2014年山东省宜居性指数。如我们将宜居性指数定义为PM2.5浓度少于75μg·m^-3的天数与全年PM2.5浓度均值的比值。

为演示S1，S2和S3，示例数据一为土壤重金属Cd，采自武汉市青山区东部的城乡交错地区，面积约为30km²，2010—2014年，在该区域分别采集样点124、45、48、55和48个，随机布点。

首先利用时空建模点构建了4个基于时空位置的土壤重金属Cd时空分布的趋势模型，并利用最小二乘法拟合了各模型中的参数，结果如图2所示。然后对样点残差计算试验变异函数值，并拟合理论变异函数模型，基于这4个时空变异理论模型和样点残差进行时空克里格插值，再加上时空趋势值，得到基于时空趋势克里格预测结果。在利用原始数据，进行时空普通克里格插值，结果如图3所示。

为演示S4和S5，示例数据二为2014年度96个监测站的PM2.5日监测平均值，将原始样点转换成4个指示值，截断值取全部样点数据的20％、40％、60％和80％分位数对应的值，对于本例，分别为34μg·m^-3(20％percentile),53μg·m^-3(40％percentile),74μg·m^-3(60％percentile),和106μg·m^-3(80％percentile)。指示转换后，原始样点数据将转换为4组指示数据。对转换后的4组指示值，分别计算其试验变异函数值，拟合理论时空变异函数模型，如图4所示。

设置网格单元大小为5000m×5000m×1day，模拟次数设置为1000次，得到1000个模拟结果，其中3次和普通时空克里格插值结果的对比如图5所示。基于上述算法获得的1000次模拟结果，得到4种时空不确定性评估结果：

(1)单点单时刻不确定性：如图6为PM2.5浓度在2014年的第1、100、200和300天超过75μg·m^-3的概率分布；

(2)单点多时刻不确定性：如图7为2014年各个月份内山东省PM2.5浓度均超过25μg·m^-3(国际大气质量达标标准)的概率空间分布；

(3)多点单时刻不确定性：如图8为山东省在2014年的第1、100、200和300天PM2.5浓度超过75μg·m^-3的概率大于0.8和0.9的空间分布；

(4)多点多时刻不确定性：如图9为2014年各个月份内山东省PM2.5浓度均超过25μg·m^-3概率大于0.5空间分布。

在时空查询与分析方面，按照S5中的三种模式，进行了实现，如图10为2014年山东省各地超过各种PM2.5污染等级的天数空间分布状况；图11为2014年全年和各季节山东省PM2.5的平均浓度；图12为2014年山东省宜居性指数空间分布；

本说明书未作详细描述的内容属于本领域专业技术人员公知的现有技术。

Claims (10)

1.一种环境变量时空分布的综合分析方法，其特征在于，它包括如下步骤：

步骤S1：基于环境变量时空样点，计算各个时空滞后距上的试验变异函数值，并拟合得到理论时空变异函数模型；

步骤S2：根据理论时空变异函数模型并采用所述环境变量时空样点数据，进行时空克里格插值，并根据时空克里格插值的结果来估算未测时空位置的环境变量值，所述时空克里格插值的结果为研究范围内的环境变量时空立方体数据；

步骤S3：依据环境变量与时空位置的关系，以时空位置为自变量，环境变量为因变量，建立环境变量值与区域位置的定量关系，形成环境变量值时空分布趋势体，再对环境变量时空样点剥离趋势后的残差进行时空克里格插值，插值结果与环境变量值时空分布趋势体的趋势值相加，得到基于时空趋势克里格估值的预测结果；

步骤S4：采用时空序贯指示模拟方法对环境变量进行模拟得到时空序贯指示模拟结果，并根据时空不确定性评估方法，对环境变量时空分布进行不确定性评估，得到地理属性超过或不超过某阈值的概率分布；

步骤S5：基于时空克里格插值获得的研究范围内的环境变量时空立方体数据，采用时空立方体数据的信息挖掘方法对环境变量的时空分布状况进行评价，获取时空分布规律。

2.根据权利要求1所述的环境变量时空分布的综合分析方法，其特征在于：所述步骤S1中的理论时空变异函数模型采用理论时空变异函数的模型库中的函数模型。

3.根据权利要求2所述的环境变量时空分布的综合分析方法，其特征在于：所述理论时空变异函数的模型库中的函数模型是根据试验变异函数散点图形状特征进行选择。

4.根据权利要求3所述的环境变量时空分布的综合分析方法，其特征在于：所述函数模型参数采用最小二乘法或遗传算法估计进行确定。

5.根据权利要求1所述的环境变量时空分布的综合分析方法，其特征在于：所述步骤2中，估算未测时空位置的环境变量值的具体方法为：

S21：构建时空插值网格，即确定时空克里格插值的空间和时间范围，以及时空克里格插值的时空粒度，所述空克里格插值的时空力度包括时空网格的空间边长及时间间隔；

S22：遍历每个时空插值网格，以时空插值网格中心为待估时空位置，搜索周围若干已知时空样点作为插值邻近点，构建克里格矩阵，解克里格矩阵，得到未测时空位置的环境变量值的时空预测结果。

6.根据权利要求1所述的环境变量时空分布的综合分析方法，其特征在于：所述步骤S3中得到基于时空趋势克里格估值的预测结果的具体方法为：

步骤S31：将时空分布趋势体表示为基于时空位置的方程；

步骤S32：对每个环境变量时空样点数据，减去时空分布趋势体，得到时空样点残差；利用时空克里格插值对时空样点残差进行时空插值，预测未采样时空位置的环境变量值，得到环境变量时空立方体数据；

步骤S33：对每个环境变量时空样点数据的时空分布趋势体与时空样点残差的插值结果进行相加，得到时空趋势克里格值的预测结果。

7.根据权利要求1所述的环境变量时空分布的综合分析方法，其特征在于：所述S4中采用时空序贯指示模拟方法对环境变量进行模拟得到时空序贯指示模拟结果的具体步骤为：

步骤S41：将环境变量时空样点数据转换成K个指示值；

步骤S42：对转换后的K个指示值，分别计算试验变异函数值，并拟合理论变异函数模型；

步骤S43：定义时空网格矩阵，采用时空序贯指示算法，对环境变量进行时空模拟；

步骤S44：设置模拟次数进行模拟，得到模拟结果。

8.根据权利要求7所述的环境变量时空分布的综合分析方法，其特征在于：所述步骤4中，对环境变量时空分布进行不确定性评估，得到地理属性超过或不超过某阈值的概率分布的具体范围为：

步骤S45：根据步骤S44的模拟结果，定义4种时空不确定性评估方案，所述4种时空不确定性评估方案包括单点单时刻不确定性方案、单点多时刻不确定性方案、多点单时刻不确定性方案和多点多时刻不确定性方案，得到地理属性超过或不超过某阈值的概率分布。

9.根据权利要求7所述的环境变量时空分布的综合分析方法，其特征在于：所述步骤S41中，将环境变量时空样点数据转换成K个指示值的具体方法为：

<mrow> <mi>I</mi> <mrow> <mo>(</mo> <mi>s</mi> <mo>,</mo> <mi>t</mi> <mo>;</mo> <msub> <mi>z</mi> <mi>c</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <mn>1</mn> <mo>,</mo> </mrow> </mtd> <mtd> <mrow> <mi>z</mi> <mrow> <mo>(</mo> <mi>s</mi> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>&amp;le;</mo> <msub> <mi>z</mi> <mi>c</mi> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mn>0</mn> <mo>,</mo> </mrow> </mtd> <mtd> <mrow> <mi>o</mi> <mi>t</mi> <mi>h</mi> <mi>e</mi> <mi>r</mi> <mi>w</mi> <mi>i</mi> <mi>s</mi> <mi>e</mi> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>,</mo> <mi>c</mi> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mo>...</mo> <mo>,</mo> <mi>K</mi> </mrow>

式中z_c为环境变量样点的若干几个截断值，I(s，t；z_c)为指示值，z(s,t)是在时空位置s,t处的环境变量值；

所述步骤S43中定义时空网格矩阵，采用时空序贯指示算法，对环境变量进行时空模拟的方法为；

步骤S43a：在定义的时空网格矩阵内随机生成模拟路线，路线按随机顺序经过所有时空网格；

步骤S43b：对在随机路线上的每一个时空网格节点的K个指示变量进行普通时空克里格插值，将得到该位置上环境变量值小于上述K个截断值的概率值，对K个概率值进行线性内插或外推，得到该位置上环境变量值的连续型累积概率密度分布函数；

步骤S43c：从上述累积概率密度函数中随机取1值作为该时空点上本次模拟的模拟值，将该模拟值按照公式4转换成指示值，并作为已测点，加入到计算随机路线上下一个时空节点的模拟计算中去；

步骤S43d：沿着随机路径，重复步骤S43b和步骤S43c，直到所有的时空网格节点都有模拟值，即完成一次模拟。

10.根据权利要求7所述的环境变量时空分布的综合分析方法，其特征在于：所述步骤5的具体方法为基于时空插值获得的环境变量时空立方体数据，根据查询条件和对象，将时空查询进行分类，根据分类模式和地理对象特征设置查询问题，得到时空分布规律。