Groundwater observation interference recognition model based on decision tree algorithm
-
摘要: 为提高地下水位观测数据中干扰事件的识别效率，利用决策树算法对宝坻等5个台站近5年的水位观测数据进行样本训练和数据验证。结果表明，决策树算法对观测系统干扰和场地环境干扰事件的分类准确率在80%以上。在大量准确的训练样本基础上，决策树算法对于各种水位干扰事件具有良好的识别效果。Abstract: To improve the identification efficiency of disturbance events in groundwater observation data, decision tree algorithm is used to perform sample training and data verification for groundwater data of Baodi and other four stations in recent five years. The results show that the classification accuracy of the decision tree algorithm for observing system interference and environmental interference events is above 80%. Based on a large number of accurate training samples, the decision tree algorithm can identify various water level interference events efficiently.
-
Keywords:
- groundwater observation /
- decision tree /
- interference identification
-
引言
地震前兆观测是分析预报地震的基础，而地下水位观测则是地震前兆观测中重要的一部分内容。经过多年的发展，目前我国地震前兆台网已实现数字化和网络化，观测数据精度和观测设备采样率都有了大幅度提高[1]。但是由于在高频度、高精度的观测数据中存在大量的干扰信息，因此，去除这些干扰因素引起的数据异常变化，是进行前兆数据分析的前提工作[2]。目前我国前兆台网日常工作模式中，需要依靠工作人员结合观测场地干扰情况，进行人工判别和处理。随着前兆观测数据日益增长，依靠人工处理方式从大量的前兆观测数据中剔除异常干扰数据，已经无法满足日常工作的需要，因此，需要研究一种高效实用的前兆观测数据干扰识别的方法。
自20世纪80年代开始，国内外就已经对前兆观测数据的干扰识别展开了一系列的研究，其中应用较为广泛的方法有方差分析方法、滤波分析方法和差分分析方法。方差分析主要是根据前兆观测数据大多具有一定的变化规律的特性，当数据变化量超过N倍均方差时，就可以被认为是干扰数据；滤波分析主要是通过数字信号处理技术提取前兆观测数据中不同周期的数据，当该周期成分的数据与理论值出现较大差距时，就可以被认为是干扰数据[3]；差分分析由于可以突出信号中的短周期成分，因此，对瞬时变化的干扰信号具有较好的检测作用[4]。
近年来不少学者将人工智能技术应用于地震学领域，利用卷积神经网络进行地震震相的自动识别[5]，Marzocchi等[6]利用贝叶斯预估值和集成的方法来分析地震数据。周雯等[7]对海南及邻近地区地震数据利用聚类分析发现了较为规律的地震分布格局。但是，大数据分析技术在地震前兆数据分析中的应用还较为较少。本文利用经典的决策树算法对水位观测数据中存在的主要干扰进行识别，以此提高前兆台网日常工作的处理效率。
1. 数据预处理
由于决策树算法的训练样本需要有一定综合性和历史跨度，因此，选取5个地震前兆台站2013—2017年间的水位观测数据，利用前兆数据跟踪分析产品，提取出工人标识出的干扰数据时段。将水位观测数据分为正常动态数据、自然环境干扰数据、场地环境干扰数据和观测系统故障数据。经过整理后将连续的水位观测数据抽分成2 128条事件数据，表1为具体事件情况。
表 1 水位观测事件总表台站 正常事件 自然环境 观测系统 场地环境 总计 宝坻台 300 50 12 35 397 静乐台 238 69 33 0 340 溧阳苏22井 119 39 9 1 168 温泉台 140 147 76 688 1 051 岫岩台 148 16 8 0 172 总计 945 321 138 724 2 128 对于各类事件数据，需首先进行预处理，之后才能进行数据特征的提取。我们为每个事件数据设置的数据时间窗口为5小时，每个事件数据含有300个观测数据样本。数据预处理的主要步骤包括：
（1）缺失值处理：检测观测数据中是否有单点的缺失数据，如果有，则将该缺失值改为它前后两个数据的均值。这个步骤主要是去除由于观测系统自身的不稳定造成的数据缺失。
（2）突跳值处理：检测观测数据中是否有超出均值很多的单点数据（阈值设为5倍均值），如果有，则将该异常值改为它前后两个数据的均值。因为单点的数据突跳并不能反映真实的观测背景异常。
（3）平稳化处理：该步骤首先将事件窗口内的所有数据点减去该窗口内的数据均值，通过去均值操作使得信号的均值为0，保证了信号在0基线上下浮动。然后通过一阶差分方法对信号进行平稳化，去掉变化缓慢的低频成分，留下变化激烈的高频成分。
2. 决策树算法
决策树算法是根据给定的训练数据集构建一个决策树模型，使它能够对实例进行正确的分类[8]。在构造决策树时，我们首先需要对每个特征进行评估，找到最具有决定性的特征。完成测试后，数据集就被划分为若干个数据子集。如果某个数据子集下的数据属于同一类型，则该数据子集已经被正确分类，如果数据子集中仍然存在不同类型的数据，那么需要重复进行数据子集的划分过程，直到所有均有相同类型的数据均在一个数据子集内。创建决策树的伪代码如下：
检测数据集中的每个子项是否属于同一分类：
如果是，则返回该类标签
如果否
寻找划分数据集最好的特征
划分数据集，创建数据子集
对于划分的每个数据子集
递归调用决策创建代码
返回数据子集标签
以上伪代码是一个递归函数，在倒数第二行直接调用了自己。
从上面的伪代码可以看出，决策树递归对原始数据集进行划分，直至所有具有相同类型的数据均在一个数据子集内。而每次对数据集的划分方法是算法的关键。我们使用参数H=–
$\displaystyle\mathop \sum \limits_{i=1}^n p\left( {{x_i}} \right){\log _2}p\left( {{x_i}} \right)$ 对数据进行划分，参数H表征了数据集的无序程度，参数值越高表示数据集中混合的数据越杂乱，其中n是分类的数目，p(xi)是每个类别数据出现的频率。划分数据集的伪代码如下：循环所有的特征：
根据选择的特征对数据集进行划分；
计算新的数据集参数H；
当参数H最小时，返回特征值；
该方法返回的特征值就是当前数据集划分的依据。
3. 特征选择
对每个事件窗口内的数据设计12个特征量。表2为数据特征表。xi代表每个窗口中第i个数据的数值。其中能量以数据的振幅绝对值来表征。
表 2 数据特征表物理量 计算公式 能量峰值 EMax=max(Exi)，Exi=|xi| 总能量 Eall=$\displaystyle\sum\limits_{i=1}^{300} {{{E_{xi}}}} $，Exi=|xi| 能量波谷值 EMin=min(Exi)，Exi=|xi| 能量中位数 Emid=mid(Exi)，Exi=|xi| 分段能量均值 Emean=$\dfrac{1}{ {300} }\displaystyle\sum\limits_{i=1}^{300} { { {E_{xi} } } } $，Exi=|xi| 能量的方差 ${\sigma ^2}=\dfrac{1}{ {300} }\displaystyle\sum\limits_1^{300} { { {({ {x_i } } - u)}^2} } $，其中u为均值 峰值变化率 LAR1=tmax/Smax，tmax 代表1 s数据段中最大值出现的位置，Smax 代表最大值 正段能量 E＋(i)=$\dfrac{1}{{30}}\displaystyle\sum\limits_{j=1}^{30} X (i*30 + j,1)$, X(i*30＋j, 1)＞0 (i=0，1，2，…，9) 负段能量 E＋(i)=$\dfrac{1}{{30}}\displaystyle\sum\limits_{j=1}^{30} X (i*30 + j,1)$, X(i*30＋j, 1)＜0 (i=0，1，2，…，9) 小波系数 ${C_{j, k} }(f, {\textit{φ} _{j, k} })={2^{ - j/2} }\displaystyle\sum\limits_{n= - \infty }^\infty {f(t){\textit{φ} } ({2^{ - j} }t - k)} $ 短时能量 $X'(i, 1)=\dfrac{1}{ {60} }\displaystyle\sum\limits_{j=1}^{60} X {(i + j, 1)^2},i=1,2,\cdots, 300$ 短时能量峰值 ${\rm {STEM} }(i,1)=\max [X'(i, 1)], i=1, \cdots, 300$ 注：能量峰值是事件窗口内变化能量最大瞬间值；总能量是在事件窗口内能量波的下方面积；能量波谷值是事件窗口内变化能量最小瞬间值；能量中位数指的是在所考虑的时间间隔内，分布中的一个数值，其可将所有不同时间的能量划分为相等的上下两部分；分段能量均值是指将每个时间窗口等分10段后，各段的能量均值；能量的方差是指每个样本值与总样本平均数之差的平方值的平均数；峰值变化率是事件窗口内最大值出现的位置除以最大值；正段能量是指所有xi＞0部分的能力总和；负段能量是指所有xi＜0部分的能力总和；小波系数是采用db5小波将事件信号展开后，得到各分量的系数的平均值；短时能量是时间窗口为1 h的数据能量值；短时能量峰值是每个短时能量窗口中的最大值 4. 多元分类模型
多元模型对所有事件进行分类标注，经过多元模型分类后，所有事件将被分成正常事件、自然环境干扰事件、场地环境干扰事件和观测系统故障事件。我们使用全部2 128条事件样本来作为分析样本，从中随机选取20%的数据作为独立验证样本，另外80%的数据作为训练数据集，使用ID3算法进行决策树构建。表3—表5分别反映了决策树算法在训练数据集、验证数据集和总数据集上的表现。
表 3 训练样本分类情况事件类型 模型分类类型 正常 自然环境 观测系统 场地环境 合计 查全率 实际类型 正常 569 133 14 36 752 79.25% 自然环境 14 221 11 21 267 82.77% 观测系统 0 1 115 0 116 99.14% 场地环境 13 93 17 453 576 78.65% 合计 596 448 157 510 1 711 查准率 95.47% 49.33% 73.25% 88.82% 表 4 检验样本分类情况事件类型 模型分类类型 正常 自然环境 观测系统 场地环境 合计 查全率 实际类型 正常 130 42 5 16 193 67.36% 自然环境 6 42 3 3 54 77.78% 观测系统 0 0 22 0 22 100% 场地环境 4 32 3 109 148 73.65% 合计 140 116 33 128 417 查准率 92.86% 36.21% 66.67% 85.16% 表 5 总样本分类情况事件类型 模型分类类型 正常 自然环境 观测系统 场地环境 合计 查全率 实际类型 正常 699 175 19 52 945 73.97% 自然环境 20 263 14 24 321 81.93% 观测系统 0 1 137 0 138 99.28% 场地环境 17 125 20 562 724 77.62% 合计 736 564 190 638 2 128 查准率 94.97% 46.63% 72.11% 88.09% 从多元模型在训练样本和检验样本的表现上来看，多元模型对于正常事件和场地环境干扰事件的分类准确率较高，而对自然环境事件和观测系统干扰事件的分类准确性较低；而从事件分类结果的完整性来看，虽然观测系统故障事件的查全率较高，但是由于样本数量过少，因此，并不具有普遍的代表性，因此，总体来看各类事件的查全率均不算太高，模型对于事件分类的性能还有提高的空间。
因此，针对每一类干扰事件，构建一个二元分类模型来提升模型效果，每个决策树模型只需要针对特定种类的事件给出是非判断。对每条事件数据均需要通过这3个分类模型给出是否属于该类事件的判断结果，最终将该事件归入概率值最高的分类，以检验模型的效果。
5. 自然环境干扰模型
我们使用全部945个正常事件样本和321个自然环境干扰事件样本作为分析样本，从中随机抽取20%的事件用来检验模型效能，剩余80%的事件作为模型的训练数据。得到检验数据验证的分类情况（表6）。
表 6 自然环境事件模型分类情况实际类别 类别
总数正确率 预测类别 正常事件 N=175 自然环境干
扰N=76正常事件 191 82.72% 158 33 自然环境干扰 60 71.67% 17 43 整体正确率 80.08% 模型对于自然环境干扰和正常事件的区分效果仍然存在一些误差，其中有33条正常事件被模型预测为自然环境干扰事件的情形。通过进一步调取台站的降雨资料进行分析，发现这33条正常事件中有20条存在降雨背景，因此，这20条正常事件可以认为是存在降雨背景的干扰。通过对台站降雨数据的对比分析后，得到检验数据验证的分类情况（表7）。
表 7 自然环境事件模型分类情况实际类别 类别
总数正确率 预测类别 正常事件 N=175 自然环境干扰 N=76 正常事件 171 82.72% 158 13 自然环境干扰 80 71.67% 17 63 整体正确率 80.08% 利用模型可以给出每个特征对于决策树分类结果的重要性，原始能量中位数与自然环境事件的判断关联度最高。从图1可以看出，当原始能量中位数在9.2—14.3之间时，事件为自然环境事件的概率较高。
6. 观测系统干扰模型
我们使用全部945个正常观测样本和138个观测系统事件样本作为分析样本，从中随机抽取20%检验模型效能，剩余80%的事件作为模型的训练数据。得到检验数据验证的分类情况（表8）。
表 8 观测系统事件模型分类情况实际类别 类别
总数正确率 预测类别 正常事件 N=196 观测系统干扰 N=23 正常事件 197 99.49% 196 1 观测系统干扰 22 100.00% 0 22 整体正确率 99.54% 模型对于观测系统故障时间具有非常高的分类准确性，通过对原始数据的分析可以看出，观测系统故障造成的数据变化特征较为明显。首先对于原始数据分析时，多数情况下原始数据的变化频率会明显增大，变化幅度显著区别于正常背景值；而对于预处理数据分析时，由于该类异常数据会被删除，因此，预处理数据会出现连续的数据缺失，导致能量偏低。因此，对于整个模型而言，观测系统故障事件和正常数据是较为容易区分的。
通过模型计算过程产出的各个特征量对于判断结果的重要程度，原始数据的总能量与观测系统干扰事件的判断关联度最高。从图2可以看出，当原始总能量超过678时，事件为观测系统故障的概率可能性很高。
7. 场地环境干扰模型
我们使用全部945个正常观测样本和724个场地环境事件样本作为分析样本，从中随机抽取20%检验模型效能。剩余80%的事件作为模型的训练数据。得到检验数据验证的分类情况（表9），显示模型具有较好的分类识别能力。
表 9 场地环境事件模型分类情况实际类别 类别总数 正确率 预测类别 正常事件 N=187 场地环境干扰 N=149 正常事件 187 94.12% 176 11 场地环境干扰 149 92.62% 11 138 整体正确率 93.45% 模型对于场地环境干扰和正常事件的区分效果较好，这主要是由于场地环境干扰事件相对于其他干扰事件而言，训练的样本数据更加丰富，而且经过对原始数据的分析发现，大部分标注的场地环境干扰事件都是周边的抽水对水位造成的趋势性变化，变化过程和形态较为统一。
通过模型计算过程产出的各个特征量对于判断结果的重要程度，原始能量峰值与场地环境事件的判断关联度最高。从图3可以看出，原始能量峰值小于4时，事件为场地环境事件的概率较高。
8. 组合模型评估
我们使用以上3个二元分类模型对总样本进行打分评估，将事件归入概率值最高的分类。
从结果来看，组合模型（表10）相较于多元分类模型而言，对于总体事件分类的查全率和准确性相差不多，主要区别在于组合分类模型对于自然环境干扰事件的分类准确性较高，但是容易将部分自然环境干扰事件错误地划分至场地环境干扰中；而多元分类模型对于自然环境干扰事件的分类查全率较高，但是容易将部分正常事件错误地划分至自然环境干扰事件中。
表 10 组合模型分类情况事件类型 模型分类 正常 自然环境 观测系统 场地环境 总计 查全率 实际分类 正常 908 4 2 31 945 96.08% 自然环境 85 125 10 101 321 38.94% 观测系统 4 2 119 13 138 86.23% 场地环境 29 24 29 642 724 88.67% 查准率 88.50% 80.65% 74.38% 81.58% 9. 总结
本文借助中国地震前兆台网产出的标签数据，利用决策树算法建立分类模型，对水位观测中常见的异常干扰事件进行自动分类。通过单个多元分类模型和多个二元分类模型组成的组合模型进行数据验证，发现决策树算法在水位观测数据异常事件的自动识别中具有一定的可靠性。
但是从结果来看，模型对于异常事件的判断仍出现了不少错误，这一方面是由于可以使用的数据仍较为有限，本次使用的模型训练数据量非常有限，这在很大程度上也影响了模型的最终性能，另一方面是由于在特征选取方面，没有进一步深入的设计更能表征异常事件变化特点的特征值，这也是影响决策树算法最后效果的关键因素之一。
-
表 1 水位观测事件总表
台站 正常事件 自然环境 观测系统 场地环境 总计 宝坻台 300 50 12 35 397 静乐台 238 69 33 0 340 溧阳苏22井 119 39 9 1 168 温泉台 140 147 76 688 1 051 岫岩台 148 16 8 0 172 总计 945 321 138 724 2 128 
下载: 导出CSV
表 2 数据特征表
物理量 计算公式 能量峰值 EMax=max(Exi)，Exi=|xi| 总能量 Eall=$\displaystyle\sum\limits_{i=1}^{300} {{{E_{xi}}}} $，Exi=|xi| 能量波谷值 EMin=min(Exi)，Exi=|xi| 能量中位数 Emid=mid(Exi)，Exi=|xi| 分段能量均值 Emean=$\dfrac{1}{ {300} }\displaystyle\sum\limits_{i=1}^{300} { { {E_{xi} } } } $，Exi=|xi| 能量的方差 ${\sigma ^2}=\dfrac{1}{ {300} }\displaystyle\sum\limits_1^{300} { { {({ {x_i } } - u)}^2} } $，其中u为均值 峰值变化率 LAR1=tmax/Smax，tmax 代表1 s数据段中最大值出现的位置，Smax 代表最大值 正段能量 E＋(i)=$\dfrac{1}{{30}}\displaystyle\sum\limits_{j=1}^{30} X (i*30 + j,1)$, X(i*30＋j, 1)＞0 (i=0，1，2，…，9) 负段能量 E＋(i)=$\dfrac{1}{{30}}\displaystyle\sum\limits_{j=1}^{30} X (i*30 + j,1)$, X(i*30＋j, 1)＜0 (i=0，1，2，…，9) 小波系数 ${C_{j, k} }(f, {\textit{φ} _{j, k} })={2^{ - j/2} }\displaystyle\sum\limits_{n= - \infty }^\infty {f(t){\textit{φ} } ({2^{ - j} }t - k)} $ 短时能量 $X'(i, 1)=\dfrac{1}{ {60} }\displaystyle\sum\limits_{j=1}^{60} X {(i + j, 1)^2},i=1,2,\cdots, 300$ 短时能量峰值 ${\rm {STEM} }(i,1)=\max [X'(i, 1)], i=1, \cdots, 300$ 注：能量峰值是事件窗口内变化能量最大瞬间值；总能量是在事件窗口内能量波的下方面积；能量波谷值是事件窗口内变化能量最小瞬间值；能量中位数指的是在所考虑的时间间隔内，分布中的一个数值，其可将所有不同时间的能量划分为相等的上下两部分；分段能量均值是指将每个时间窗口等分10段后，各段的能量均值；能量的方差是指每个样本值与总样本平均数之差的平方值的平均数；峰值变化率是事件窗口内最大值出现的位置除以最大值；正段能量是指所有xi＞0部分的能力总和；负段能量是指所有xi＜0部分的能力总和；小波系数是采用db5小波将事件信号展开后，得到各分量的系数的平均值；短时能量是时间窗口为1 h的数据能量值；短时能量峰值是每个短时能量窗口中的最大值
下载: 导出CSV
表 3 训练样本分类情况
事件类型 模型分类类型 正常 自然环境 观测系统 场地环境 合计 查全率 实际类型 正常 569 133 14 36 752 79.25% 自然环境 14 221 11 21 267 82.77% 观测系统 0 1 115 0 116 99.14% 场地环境 13 93 17 453 576 78.65% 合计 596 448 157 510 1 711 查准率 95.47% 49.33% 73.25% 88.82%
下载: 导出CSV
表 4 检验样本分类情况
事件类型 模型分类类型 正常 自然环境 观测系统 场地环境 合计 查全率 实际类型 正常 130 42 5 16 193 67.36% 自然环境 6 42 3 3 54 77.78% 观测系统 0 0 22 0 22 100% 场地环境 4 32 3 109 148 73.65% 合计 140 116 33 128 417 查准率 92.86% 36.21% 66.67% 85.16%
下载: 导出CSV
表 5 总样本分类情况
事件类型 模型分类类型 正常 自然环境 观测系统 场地环境 合计 查全率 实际类型 正常 699 175 19 52 945 73.97% 自然环境 20 263 14 24 321 81.93% 观测系统 0 1 137 0 138 99.28% 场地环境 17 125 20 562 724 77.62% 合计 736 564 190 638 2 128 查准率 94.97% 46.63% 72.11% 88.09%
下载: 导出CSV
表 6 自然环境事件模型分类情况
实际类别 类别
总数正确率 预测类别 正常事件 N=175 自然环境干
扰N=76正常事件 191 82.72% 158 33 自然环境干扰 60 71.67% 17 43 整体正确率 80.08%
下载: 导出CSV
表 7 自然环境事件模型分类情况
实际类别 类别
总数正确率 预测类别 正常事件 N=175 自然环境干扰 N=76 正常事件 171 82.72% 158 13 自然环境干扰 80 71.67% 17 63 整体正确率 80.08%
下载: 导出CSV
表 8 观测系统事件模型分类情况
实际类别 类别
总数正确率 预测类别 正常事件 N=196 观测系统干扰 N=23 正常事件 197 99.49% 196 1 观测系统干扰 22 100.00% 0 22 整体正确率 99.54%
下载: 导出CSV
表 9 场地环境事件模型分类情况
实际类别 类别总数 正确率 预测类别 正常事件 N=187 场地环境干扰 N=149 正常事件 187 94.12% 176 11 场地环境干扰 149 92.62% 11 138 整体正确率 93.45%
下载: 导出CSV
表 10 组合模型分类情况
事件类型 模型分类 正常 自然环境 观测系统 场地环境 总计 查全率 实际分类 正常 908 4 2 31 945 96.08% 自然环境 85 125 10 101 321 38.94% 观测系统 4 2 119 13 138 86.23% 场地环境 29 24 29 642 724 88.67% 查准率 88.50% 80.65% 74.38% 81.58%
下载: 导出CSV
-
[1] 李正媛,熊道慧,刘高川. 基于大数据挖掘的地震前兆台网观测数据跟踪分析[J]. 地震地磁观测与研究,2016,37(3):1-6 [2] 陈敏,郭欣,郑许东. 重庆台倾斜潮汐形变干扰异常与震兆异常的识别[J]. 大地测量与地球动力学,2007,27(1):79-83 [3] 林纪曾. 观测数据的数学处理[M]. 地震出版社, 1979: 39-69 [4] 胡哲,朱涛,王兰炜. 地震地电阻率数据处理方法[J]. 地球物理学进展,2011,26(1):79-88 doi: 10.3969/j.issn.1004-2903.2011.01.008 [5] Thibaut P,Michaël G,Marine D. Convolutional neural network for earthquake detection and location[J]. Sci. Adv.,2018,4:e1700578
[6] Marzocchi W,Zechar J,Jordan T H. Bayesian forecast evaluation and ensemble earthquake forecasting[J]. Bull. Seismol. Soc. Amer,2012,102(6):2574-2584 doi: 10.1785/0120110327
[7] 周雯,焦晨,苏荣托雅. 海南及邻近区域地震活动性聚类分析[J]. 地震地磁观测与研究,2013,4(5/6):75-79 [8] 李孝伟,陈福才,李邵梅. 基于分类规则的C4.5决策树改进算法[J]. 计算机工程与设计,2013,34(12):4321-4325 doi: 10.3969/j.issn.1000-7024.2013.12.046 -
期刊类型引用(0)
其他类型引用(1)
计量
- 文章访问数: 542
- HTML全文浏览量: 166
- PDF下载量: 8
- 被引次数: 1
