CN105787283A - 一种基于时空相关性的土遗址监测数据修正拟合方法 - Google Patents

Abstract

本发明公开了一种基于时空相关性的土遗址监测数据修正拟合方法，利用历史数据的预测值和当前检测数据进行时空相关性分析，进而判定当前数据的正确性，并对噪声数据进行修正、缺失值予以拟合。本发明用于土遗址保护领域，实现土遗址的智能状态监测，特别是针对大量的环境数据和土遗址本体数据，能够准确实现噪声数据的发现和修正，缺失数据趋势的拟合，为土遗址风险预警和有效管理提供科学依据。

Description

一种基于时空相关性的土遗址监测数据修正拟合方法

技术领域

本发明属于数据挖掘领域技术领域，具体涉及长城土遗址保护领域中基于时空相关性的土遗址环境数据处理方法，可应用于数据去除噪声和缺失数据拟合，以达到准确评估风险因素的目的。

背景技术

在土遗址保护领域中，为监测土遗址的环境数据和本体数据，通常在土遗址表面布置大量传感器节点，如：雨量传感器节点、土壤温湿度传感器节点、土壤盐分传感器节点等等，由于环境恶劣和人为破坏以及节点本身稳定性等原因会导致在数据采集阶段会产生部分噪声数据和数据缺失。另外，节点采集数据之后，网关会对节点采集的原始数据进行传输，通常采用运营商的移动蜂窝网络传输数据，由于野外信号质量以及网络拥塞的原因会导致在数据传输阶段也会有数据丢失。

噪声数据和缺失数据会影响土遗址的风险量化评估。噪声数据会掩盖土遗址真实数据的变化趋势，可能导致错误的风险评估，阻碍风险源的发现，影响土遗址保护。缺失数据，也可能导致遗漏风险量化评估的关键数据，不能达到土遗址状态检测的目的。

目前的数据噪声修正方法和缺失数据拟合方法存在很多不足：针对噪声数据，当前多数方法采用的是根据经验对不同类型的数据利用相应的阈值过滤；这种处理方法一方面会错误的将孤立点过滤掉，然而孤立点往往是土遗址关键的数据；另一方面，有一部分噪声数据在正常范围内，而用阈值过滤方法无法发现。针对数据缺失，通常由人工根据数据曲线的趋势补全数据，人工补全的方法，由于人的主观性会导致补全的数据有偏差。并且在土遗址数据检测中，数据量为数十万，甚至数百万条，人工补全数据没有可行性。

发明内容

针对在土遗址保护中数据噪声修正方法和缺失数据拟合方法准确性差和效率低的不足，本发明的目的在于，提供一种基于时空相关性的土遗址监测数据修正拟合方法，以解决现有方法准确性差、效率低的缺点，提高了土遗址监测数据的质量，为土遗址风险预警和有效管理提供可靠的科学根据。

为了实现上述任务，本发明采用以下技术方案：

一种基于时空相关性的土遗址监测数据修正拟合方法，包括以下步骤：

步骤一，提取监测数据的历史数据，利用预测算法根据历史数据计算当前的预测值；

步骤二，读取当前的监测值，判断监测值是否在传感器的量程范围内，若监测值缺失，则执行步骤五；若监测值不在范围内，则判定当前监测值为噪声数据；若监测值在正常范围内，则进行下一步处理；

步骤三，计算步骤一得到的预测值与步骤二得到的监测值之间的偏差，分析当前监测值的时间相关性；若符合时间相关性，则当前监测数据为正常值，若不符合时间相关性，则进行下一步处理；

步骤四，计算当前监测值与其空间临近节点所采集数据的偏差，并计算支持率，以分析当前监测值的空间相关性；若符合空间相关性，则判定当前的监测值为离群值；若不符合空间相关性，则判定当前的监测值为错误值，即噪声数据；

步骤五，经过步骤一至步骤四的处理，若当前的监测值为噪声数据，则对噪声数据进行修正；若当前的监测值缺失，则为缺失数据赋予拟合值。

进一步地，步骤一中所述的预测算法包括：

步骤S10，计算平滑系数α，计算公式为：

$\mathrm{\α}=1-{\exp }^{\[\frac{\ln (1-w\%)}{n}\]}---\left(1\right)$

上式中，ln表示自然对数；w表示历史数据在计算中所占权重的百分比；n表示所取历史数据中时间序列包含的数据点的个数；

步骤S11，计算预测参数a_i、b_i和c_i，公式如下：

$\left\{\begin{array}{c}{a}_i=3S_i^{\left(1\right)}-3S_i^{\left(2\right)}+S_i^{\left(3\right)}\\ b_i=\frac{\mathrm{\α}}{2{(1-\mathrm{\α})}^2}\[(6-5\mathrm{\α})S_i^{\left(1\right)}-(10-8\mathrm{\α})S_i^{\left(2\right)}+(4-3\mathrm{\α})S_i^{\left(3\right)}\]\\ c_i=\frac{{\mathrm{\α}}^2}{2{(1-\mathrm{\α})}^2}(S_i^{\left(1\right)}-2S_i^{\left(2\right)}+S_i^{\left(3\right)})\end{array}\right.---\left(2\right)$

步骤S12，将历史数据代入下面公式3的模型中：

$\left\{\begin{array}{c}{S}_i^{\left(1\right)}={\mathrm{\αX}}_i+(1-\mathrm{\α})S_i^{\left(1\right)}\\ S_i^{\left(2\right)}={\mathrm{\αS}}_i^{\left(1\right)}+(1-\mathrm{\α})S_{i-1}^{\left(2\right)}\\ S_i^{\left(3\right)}={\mathrm{\αS}}_i^{\left(2\right)}+(1-\mathrm{\α})S_{i-1}^{\left(3\right)}\end{array}\right.---\left(3\right)$

在公式(2)和公式(3)中，X_i为所取历史数据中时间序列{X_i}的第i期实际数据，为第i期数据的一次指数平滑值；为第i期数据的二次指数平滑值；为第i期数据的三次指数平滑值，

步骤S13，计算第i+m期，即当前的预测值公式如下：

${\hat{X}}_{i+m}=a_i+b_{i}m+c_i{m}^2---\left(4\right)$

上式中，m为预测步长，取值为正整数。

进一步地，所述的步骤三中，分析当前监测值的时间相关性的方法包括：

计算当前的监测值与当前的预测值的绝对偏差D，并记时间相关性阈值为T，若存在D＜T,则判定当前的监测值符合时间相关性，否则判定当前的监测值不符合时间相关性。

进一步地，所述的步骤四中，计算支持率的过程包括：

步骤S40,计算当前节点的监测值所在与其空间相邻的每一个节点采集数据的绝对偏差，公式为：

D_j＇＝Y_j-X_i(6)

上式中，D_j＇为编号为j的节点与当前节点监测值X_i的绝对偏差，Y_j为编号为j的节点采集的数据；

步骤S41，计算支持率：若D＇_j＜T＇，则表示编号为j的节点采集的数据支持当前节点的监测值，若D＇_j＞T＇，则表示编号为j的节点采集的数据不支持当前节点的监测值；其中T＇为偏差阈值；则监测值在相邻节点中的支持率P(X_i)为：

P(X_i)＝n/(n+m)(7)

上式中，n为当前节点的所有空间临近节点中支持监测值的节点数量，m为当前节点的所有空间临近节点中不支持监测值的节点数量。

进一步地，所述的步骤四中，分析当前监测值的方法包括：

步骤S42，判断如果P(X_i)≥0.5，则当前监测符合空间相关性；如果P(X_i)＜0.5，则当前监测值不符合空间相关性。

本发明与现有技术相比具有以下技术特点：

第一，本发明利用通过根据历史数据计算预测值，基于土遗址检测数据的时间相关性，根据近期数据比相对较远期的数据对预测值的影响更大的特点，利用非线性变化趋势的时间序列，准确计算预测值，并可靠地去除和修正噪声数据，拟合缺失数据的变化趋势，更加符合土遗址真实监测状态；

第二，本发明数据修正拟合方法对应的程序，对于土遗址状态监测中海量的环境数据和本体数据，在无监督的情况下，能够准确高效的运行，很好的解决了现有人工处理数据中存在的不足；

第三，本发明用于土遗址保护领域，可实现土遗址的智能状态监测，特别是针对大量的环境数据和土遗址本体数据，能够准确实现噪声数据的发现和修正，缺失数据趋势的拟合，为土遗址风险预警和有效管理提供科学依据。

附图说明

图1为本发明方法的流程图；

图2为表2中数据绘制成的折线图；

图3为表4中数据绘制成的趋势图；

图4为表5中数据绘制成的趋势图；

具体实施方式

结合本发明方法的流程图1，本发明方法的具体内容如下：

一种基于时空相关性的土遗址监测数据修正拟合方法，包括以下步骤：

步骤一，提取当前节点监测数据的历史数据，利用预测算法根据历史数据计算当前的预测值；这里的历史数据可根据实际处理要求选择一段时间的历史数据。具体的一种实现过程包括：

步骤S10，读取历史数据，计算反映不同时期历史数据在指数平滑值中所占的比重，即平滑系数α，计算公式为：

$\mathrm{\α}=1-{\exp }^{\[\frac{\ln (1-w\%)}{n}\]}---\left(1\right)$

上式中，ln表示自然对数；w表示历史数据在计算中所占权重的百分比；n表示所取历史数据中时间序列包含的数据点的个数；例如，在长城土遗址的数据预测中，过去两天的历史数据在计算中所占的权重是95％，由于节点数据是3小时/条，那么w＝95,n＝16,计算得α＝0.17。

步骤S11，计算预测参数a_i、b_i和c_i，公式如下：

$\left\{\begin{array}{c}{a}_i=3S_i^{\left(1\right)}-3S_i^{\left(2\right)}+S_i^{\left(3\right)}\\ b_i=\frac{\mathrm{\α}}{2{(1-\mathrm{\α})}^2}\[(6-5\mathrm{\α})S_i^{\left(1\right)}-(10-8\mathrm{\α})S_i^{\left(2\right)}+(4-3\mathrm{\α})S_i^{\left(3\right)}\]\\ c_i=\frac{{\mathrm{\α}}^2}{2{(1-\mathrm{\α})}^2}(S_i^{\left(1\right)}-2S_i^{\left(2\right)}+S_i^{\left(3\right)})\end{array}\right.---\left(2\right)$

其中平滑初值

步骤S12，将历史数据代入下面公式3的模型中：

$\left\{\begin{array}{c}{S}_i^{\left(1\right)}={\mathrm{\αX}}_i+(1-\mathrm{\α})S_i^{\left(1\right)}\\ S_i^{\left(2\right)}={\mathrm{\αS}}_i^{\left(1\right)}+(1-\mathrm{\α})S_{i-1}^{\left(2\right)}\\ S_i^{\left(3\right)}={\mathrm{\αS}}_i^{\left(2\right)}+(1-\mathrm{\α})S_{i-1}^{\left(3\right)}\end{array}\right.---\left(3\right)$

在公式(2)和公式(3)中，X_i为所取历史数据中时间序列{X_i}的第i期实际数据，为第i期数据的一次指数平滑值；为第i期数据的二次指数平滑值；为第i期数据的三次指数平滑值，α为步骤S10计算的平滑系数，且α∈(0,1)；

步骤S13，计算第i+m期，即当前的预测值公式如下：

${\hat{X}}_{i+m}=a_i+b_{i}m+c_i{m}^2---\left(4\right)$

上式中，m为预测步长，即需要预测的期数与当前期数的间隔数，取值为正整数1，2，3，...。如m＝1时表示一步预测；另外记

步骤二，读取当前节点的实时监测值，判断监测值是否在传感器的量程范围内，若监测值缺失，则执行步骤五；若监测值不在范围内，则判定当前监测值为噪声数据；若监测值在正常范围内，则进行下一步处理；

读取当前传感器节点的实时监测值X_i，而不同类型传感器在量程范围不尽相同，设当前监测值对应的传感器量程范围为[min,max]。若监测值缺失，则进入步骤五；若监测值不在范围内，则该数据为噪声数据；若监测值在正常范围[min,max]内，则进入步骤三；

步骤三，计算步骤一得到的预测值与步骤二得到的监测值之间的偏差，分析当前监测值的时间相关性；若符合时间相关性，则当前监测数据为正常值，若不符合时间相关性，则进行下一步处理；

比较当前数据与历史数据之间是否满足时间相关性，即计算预测值与监测值的绝对偏差D：

$D=|{\hat{X}}_i-X_i|---\left(5\right)$

记时间相关性阈值为T，若存在D＜T,则判定当前的监测值符合时间相关性，该数据为正常值，无需处理；否则判定当前的监测值不符合时间相关性。

步骤四，计算当前监测值与其空间临近节点所采集数据的偏差，并计算支持率，以分析当前监测值的空间相关性；若符合空间相关性，则判定当前的监测值为离群值；若不符合空间相关性，则判定当前的监测值为错误值，即噪声数据；具体过程包括：

步骤S40,计算当前节点的监测值所在与其空间相邻的每一个节点采集数据的绝对偏差，公式为：

D_j′＝Y_j-X_i(6)

上式中，D_j′为编号为j的节点与当前节点监测值X_i的绝对偏差，Y_j为编号为j的节点采集的数据；

步骤S41，计算支持率：若D′_j＜T′，则表示编号为j的节点采集的数据支持当前节点的监测值，若D′_j＞T′，则表示编号为j的节点采集的数据不支持当前节点的监测值；其中T′为偏差阈值；则监测值在相邻节点中的支持率P(X_i)为：

P(X_i)＝n/(n+m)(7)

上式中，n为当前节点的所有空间临近节点中支持监测值的节点数量，m为当前节点的所有空间临近节点中不支持监测值的节点数量。

步骤S42，空间相关性判断：如果P(X_i)≥0.5，则当前监测符合空间相关性，该数据为离群值；如果P(X_i)＜0.5，则当前监测值不符合空间相关性，该数据为错误值，即噪声数据。

步骤五，数据修正和拟合：经过步骤一至步骤四的处理，若当前的监测值为噪声数据，则对噪声数据进行修正，即利用预测值对监测值进行修正，将预测值赋予当前的监测值；若当前的监测值缺失，则为缺失数据赋予拟合值，即利用预测值拟合数据走势。

仿真实验：

一、数据去噪对比

在陕西榆林镇北台长城土遗址监测数据库中截取一段数据，如表1所示，时间201511041915表示2015年11月04日19:15，数据类型为19表示30cm温度数据。分别对比现有的阈值法和本发明方法对表1中数据噪声处理的结果。

表1496号节点30cm温度数据

(1)使用现有的阈值法处理表1中的数据。土温传感器节点实际测量范围为[-35,85]，显然范围过大，根据榆林长期历史温度变化情况和经验适当缩小区间为[-35,50],发现表1没有数据在该区间外。处理结果：表1中没有噪声数据。

(2)使用本发明方法处理表1中的数据。将表中数据代入步骤一中，并分析时间相关性和空间相关性，发现时间为201511052240的数据不符合时空相关性。处理结果：时间为201511052240的数据为噪声数据，用预测值修正该值，其余均为正常数据不做改变，如表2所示。

表2本发明方法处理后的数据

为了分析上述两种方法处理的结果，将表2中的数据绘制成折线图，如图2所示。根据一天中温度变化的规律，最低温度一般为5点-6点的日出前后，而最高温度一般为1点-2点。时间为201511052240的原始数据，在22:40温度却达到一天中的最高值，与客观规律不符，应该为噪声数据，与本发明中方法得出的结论一致。并且对该数据修正后，曲线呈持续下降的趋势，符合一般规律。

二、缺失数据处理

同样，在陕西榆林镇北台长城土遗址监测数据库中截取一段数据如表3，由于网关故障缺失4条数据，对于缺失数据，分别对比人工补全方法和本发明方法处理的结果。

表3496号节点30cm温度数据

(1)人工补全方法根据数据的规律和经验，通过对数据趋势的分析，人工对数据进行补全，如表4：

表4人工补全数据的结果

时间	原始数据	处理后的数据	节点号	数据类型
					201511131915	11	11	496	19
201511132217	9.3	9.3	496	19
					201511140120	8.7	8.7	496	19
201511140423	8.1	8.1	496	19
					201511140726	8.8	8.8	496	19
201511141028	缺失	9.2	496	19
					201511141331	缺失	10.5	496	19
201511141634	缺失	9.7	496	19
					201511141937	缺失	9.5	496	19
201511142240	9.3	9.3	496	19
					201511150142	8.9	8.9	496	19
201511150445	7.6	7.6	496	19
					201511150748	8.5	8.5	496	19

(2)利用本发明方法计算出预测值，将预测值赋予缺失点，如表5所示：

表5本发明处理的结果

为分析两种方法的处理结果，分别画出表4和表5的趋势图，如图3和图4所示。虽然二者补全的值不一样，但曲线的趋势表现一致，基本符合实际情况。但是人工补全方法明显的缺点是对参与人的主观性依赖性强，不同人补全的结果不尽相同，最重要的是在大量的数据处理中，人工补全方法的效率极低，完全无法实施。而本发明方法可以高效的处理大量数据。

Claims (5)

1.一种基于时空相关性的土遗址监测数据修正拟合方法，其特征在于，包括以下步骤：

步骤一，提取监测数据的历史数据，利用预测算法根据历史数据计算当前的预测值；

步骤二，读取当前的监测值，判断监测值是否在传感器的量程范围内，若监测值缺失，则执行步骤五；若监测值不在范围内，则判定当前监测值为噪声数据；若监测值在正常范围内，则进行下一步处理；

步骤三，计算步骤一得到的预测值与步骤二得到的监测值之间的偏差，分析当前监测值的时间相关性；若符合时间相关性，则当前监测数据为正常值，若不符合时间相关性，则进行下一步处理；

步骤四，计算当前监测值与其空间临近节点所采集数据的偏差，并计算支持率，以分析当前监测值的空间相关性；若符合空间相关性，则判定当前的监测值为离群值；若不符合空间相关性，则判定当前的监测值为错误值，即噪声数据；

步骤五，经过步骤一至步骤四的处理，若当前的监测值为噪声数据，则对噪声数据进行修正；若当前的监测值缺失，则为缺失数据赋予拟合值。

2.如权利要求1所述的基于时空相关性的土遗址监测数据修正拟合方法，其特征在于，步骤一中所述的预测算法包括：

步骤S10，计算平滑系数α，计算公式为：

上式中，ln表示自然对数；w表示历史数据在计算中所占权重的百分比；n表示所取历史数据中时间序列包含的数据点的个数；

步骤S11，计算预测参数a_i、b_i和c_i，公式如下：

步骤S12，将历史数据代入下面公式3的模型中：

在公式(2)和公式(3)中，X_i为所取历史数据中时间序列{X_i}的第i期实际数据，为第i期数据的一次指数平滑值；为第i期数据的二次指数平滑值；为第i期数据的三次指数平滑值，

步骤S13，计算第i-m期，即当前的预测值公式如下：

上式中，m为预测步长，取值为正整数。

3.如权利要求1所述的所述的基于时空相关性的土遗址监测数据修正拟合方法，其特征在于，所述的步骤三中，分析当前监测值的时间相关性的方法包括：

计算当前的监测值与当前的预测值的绝对偏差D，并记时间相关性阈值为T，若存在D＜T,则判定当前的监测值符合时间相关性，否则判定当前的监测值不符合时间相关性。

4.如权利要求3所述的所述的基于时空相关性的土遗址监测数据修正拟合方法，其特征在于，所述的步骤四中，计算支持率的过程包括：

步骤S40,计算当前节点的监测值所在与其空间相邻的每一个节点采集数据的绝对偏差，公式为：

D_j′＝Y_j-X_i(6)

上式中，D_j′为编号为j的节点与当前节点监测值X_i的绝对偏差，Y_j为编号为j的节点采集的数据；

步骤S41，计算支持率：若D′_j＜T′，则表示编号为j的节点采集的数据支持当前节点的监测值，若D′_j＞T′，则表示编号为j的节点采集的数据不支持当前节点的监测值；其中T′为偏差阈值；则监测值在相邻节点中的支持率P(X_i)为：

P(X_i)＝n/(n-m)(7)上式中，n为当前节点的所有空间临近节点中支持监测值的节点数量，m为当前节点的所有空间临近节点中不支持监测值的节点数量。

5.如权利要求4所述的所述的基于时空相关性的土遗址监测数据修正拟合方法，其特征在于，所述的步骤四中，分析当前监测值的方法包括：

步骤S42，判断如果P(X_i)≥0.5，则当前监测符合空间相关性；

如果P(X_i)＜0.5，则当前监测值不符合空间相关性。