基于数据的前兆仪器分辨力检测方法研究

成万里，王秀英

(1.河南省地震局，河南 郑州 450016；2.中国地震局地壳应力研究所，北京 100085)

0 引言

在地震科学领域，特别是“十五”数字地震观测网络项目以来，地震前兆观测数据在时间积累、采样频率、观测布点、学科种类四个维度上激增，各学科产出的海量数据中蕴含着大量价值信息,同时也包含大量的噪声信息，要从噪声数据中分辨出价值信息，必须正确认识仪器分辨力。刘洋君等认为地震仪器的信噪比大小决定测量信号的真实度,因而各种噪声之和的大小决定测量的精度、灵敏度及检测下限。要降低检测下限,提高检测精度,首先应该设法降低各种噪声的水平[1]。目前关于信号降噪的研究很多，包括小波变换、经验模态分解、滤波技术等[2-3]，利用小波方法,将信号分解为不同频段，实现信噪分离。但是,如果对信号的特征频段不够了解，则难以获得满意的结果。经验模态分解是基于信号本身的时间尺度特征，无需选择基函数就可把复杂信号由精细尺度到粗大尺度分解为若干本征模态分量，适合对非平稳、非线性信号进行平稳化处理。对于由偶然因素引起的脉冲干扰噪声很容易通过限幅滤波法、中值滤波法和算术平均法等方法去除。可以看出，上述降噪方法都是从原始信号中去除这些加噪声信息的盲源分离方法。通常情况下，特定周期的环境动态变化与噪声信息的周期或尺度无显著差异时，盲源分离方法很难真正实现信噪分离[4]。因此，噪声对数据应用的影响大小与环境动态变化特征有关，如分钟值尺度的噪声信号可能与小时周期动态变化尺度相当，两者难以分辨，但其可能对日周期环境动态变化无明显影响。据此特点，文章以观测数据短期变化能否有效区分长时间环境变化，来描述噪声信号与环境动态变化速率的差异性，并以此来度量仪器的相对分辨力。这与通常在实验室中测试完成的仪器分辨力有一定的差别。事实上，实验室的测试环境与实际观测环境就有很大差别。

收集2006年至 2017年全国595个测点近1TB的地下水温数据，以观测数据短期变化能否有效区分长时间环境变化为度量依据，尝试用大数据统计分析方法从这巨量数据中估算3类地热观测仪器的分辨力。一方面探索从数据中解决实际监测应用问题的新方法；同时，用全部实测数据参与计算，避免偶然因素对观测数据和计算结果的影响。

1 方法与数据

1.1 方法原理分析

仪器观测数据携带的信息IT，可分解为仪器自身影响的信息ID和环境动态变化影响信息IE。其中，环境影响包含地震、地球物理场及各种干扰等，即：

IT=ID+IE

。

(1)

大多数情况，环境信息的变化需要一定的持续时间，认为在较短时间内观测环境IE不变或变化极小可以忽略不计，则较短时间的变化主要由仪器自身因素引起，式(1)表示为IT≈ID，表明较短时间内数据的变化仅反映仪器自身性能的变化。

较长时间的变化反映出各种因素叠加的环境变化，如果环境变化远大于仪器性能(分辨力)变化，则式(1)表示为IT≈IE，说明仪器可以满足观测环境变化的要求。如果较长时间的环境变化与较短时间的仪器性能变化接近，则仪器无法满足观测环境的要求。

按照上述方法原理，将短时间内观测数据的波动看作主要是由仪器自身分辨力造成。如果在此较短时间观测环境的波动与仪器分辨力接近，得到的仪器分辨力仍可反映仪器自身的特性。

利用该方法计算分辨力时需解决两个问题。第一是确定观测数据满足IT≈ID的时间窗长度；第二是如何避免计算时偶然因素的影响。针对问题一，采用不同的时间窗计算数据波动情况，综合大量数据选择一个合理的计算时间窗；针对问题二，全部数据参与计算，抵消短时偶然因素对数据的影响。如，利用同一测点全部数据参与计算，降低局部可能存在的干扰或其他影响；利用同类观测仪器、不同观测点的全部数据，避免单一台站或测点可能出现的偏差或偶发情况。因此，对于两个问题，最终转化为大量数据计算问题，通过全部数据参与计算，挖掘出数据的动态变化特征，对仪器技术指标进行合理评价。

计算仪器分辨力、环境变化等反映数据动态变化的指标，称为特征量，可以选取基于统计学的均值、标准差、中位数、极差等量，也可以根据具体的研究目标特别设计，视具体任务而定[5]。

1.2 研究数据

深井温度短期变化非常微小，波动通常不大于0.000 1 ℃，年漂移变化幅度远小于0.1 ℃～ 0.2 ℃[6-7]。地下水温的这种变化特点，使式(1)中短期瞬时变化与长期环境变化的关系更清楚。选取地下水温数据作为实验数据，取“十五”前兆数据库中2006年至2017年所有水温观测数据，共计595套仪器约128万条日数据记录，每条日数据记录长度为1440的观测时间序列。主要包含Szw-1a、Szw-1和ZKGD3000三种型号的观测仪器，按表层(4311)、中层(4312)、深层(4313)三段水温划分测项代码。因部分测点停测或数据不可用，整理出有效参与计算的数据信息如表1所示。

由于选取观测时间跨度较长，仪器观测起止时间不统一，缺测情况较多，因此研究采取小时段分段的方法选取数据，不考虑数据的完整率和连续率，按照大数据分析的思想，只要有数据片段均可参与计算。

2 计算时间窗确定

据式(1)，确定一个合理的时间窗，使得在此计算时段内，环境变化尽量足够小。原理上这个时间越短越好，但时间太短，分钟采样数据参与运算的太少，会由于偶然数据问题对计算结果产生较大影响。

表1 全国地下水温观测数据基本信息Table 1 Basic information of underground water temperature observation data in China

因此，选择的时间窗虽尽量短，但确保有足够的数据参与计算。

为确定实际仪器分辨力评估计算时的时间窗，以沂南台的水温观测数据进行计算选取。沂南台的观测环境较好，日变幅、月变幅、年变幅分别不超过0.004 ℃、0.01 ℃和0.2 ℃，对于水温观测环境的要求具有代表性。

确定时间窗需要描述数据动态变化的指标。用一小段数据X的标准差SD(X)来描述仪器自身影响的变化，用滑动平均法(n=3)平滑后X的最大变幅RANG(X)表示水温受环境影响的变化。在数据上下均匀波动无趋势变化的理想状态下， RANG(X)/2略小于SD(X)。根据这一数学特性，当RANG(X)/2明显大于SD(X)时，表明此时环境影响与仪器自身影响的信息可以区分开。通过分别计算不同时间窗(即不同长度的X)下的SD(X)与RANG(X)/2，找到刚好RANG(X)/2大于SD(X)的时间窗。如图1所示，设计沂南台水温数据从10 min～120 min12个Δt的时间窗进行对比。

图1 沂南台SD(X)与RANG(X)/2的变化关系Fig.1 The relationship between SD (X) and RANG (X)/2 in Yinan Station

由图1看出，20 min之前，RANG(X)/2足够小，很难与SD(X)区分开。由于环境变化与仪器自身影响变化速率不同，20 min后，RANG(X)/2明显大于SD(X)。说明沂南台水温数据Δt 取值20 min附近能够较好地区分出仪器自噪声和环境变化。需要说明的是，在计算水温动态变化特征时，时间窗选取采用5年以上的数据全部参与计算，通过聚类方法排除偶然因素后的计算结果。如计算长度为15 min的SD(X)与RANG(X)/2，将全部数据中满足15分钟长度的数据段全部计算出结果，通过聚类方法，将密度最大、个数最多的一类视为正常，将该类结果的均值作为最后结果。图1表明，环境变化和仪器自身影响变化都是一种稳定规律的变化，说明计算结果可信。

3 计算结果

按照如上方法原理，时间窗确定后，还需解决计算结果如何避免偶然因素影响的问题。空间维度上，有些测点的数据变化不正常；时间维度上，一个测点有些时段数据变化不正常。研究通过两步聚类分析来分别筛选掉这些异常数据。聚类分析是用机器学习的方法将数据集中在某些方面相似的数据成员进行自动分类组织的过程。因此无论在空间维度还是时间维度，都可以将绝大部分正常时段的数据聚为一类，将少量异常时段的数据聚为一类。具体过程与结果如下。

以Szw-1a型仪器、测项代码为4313的中层水温209个测点573 250条记录为例(每条记录包含24 h的数据)。分别计算3个相隔长度为6 h的均值差来描述环境变化情况。在每个测点内按照均值差进行均值聚类，记录个数最多的一类作为代表该测点正常环境变化的记录，并求均值作为该测点的环境变化特征值。表2是其中台站代码为16002测点1的3 503条记录的聚类中心。显然，聚类结果最多的第1类，其聚类中心的3个值差异最小，样本个数最多，数据比较符合逻辑，视为正常记录，将该类的聚类中心3个值作为该测点的环境变化特征。

表2 台站代码为16002测点1的均值差聚类中心Table 2 The mean difference clustering center of the measurement points 1 in the station with code 16002

分别筛选出209个测点的环境变化正常时段样本，计算各测点环境变化特征。从空间维度根据测点环境变化特征对209个测点进行二次聚类，其聚类结果如图2所示(图中Y轴仪器相对分辨力反映的是209个测点仪器自身因素影响的变化整体情况，X轴反映的是209个测点6小时内受环境影响变化的整体情况)。二次聚类中1类的113个测点均值差聚类中心为0.000 55，说明该类测点受观测环境影响较小，认为这类测点短时间内环境温度相对稳定。将描述其短时间变化的标准差作为该类测点的相对分辨力。经计算该类113个测点时间窗内的标准差均值为0.000 38。

图2 szw-1a型水温仪观测数据环境变化特征二次聚类结果Fig.2 Quadratic clustering results of the environmental change characteristics of szw-1a type water temperature meter observation data

按照该方法，计算其他类型仪器和测项，结果如图3所示。不同型号的仪器及不同的入水深度在相对分辨力上均有一定的差异性。图中结果表明该计算结果与仪器实际性能较吻合。

图3 三套仪器在不同测项的分辨力Fig.3 Resolution of three sets of instruments in different terms

4 讨论分析

在实验方法设计时，时间窗的选取是重要步骤，其对实验结果有影响，是一个求解分离观测数据的宏观变化与微观变化最优解的过程。研究采取大量实验对比方法选取的最优值，有待于通过设计算法完成的最优解来验证其精确度。

分析表明，受观测环境影响较小的测点和观测时段的样本更适合用于计算仪器的分辨力，这类测点对仪器分辨力要求更高。对于相当一部分环境日变化较小的测点，仪器的分辨力不能分辨出日变化以内的环境变化。从图2看出，测点的分辨力和等距环境变化沿着聚类中心呈一定斜率展布，其中部分至原点斜率大于1的测点，无论分辨力高低，因仪器自噪声水平已经大于长期环境变化，该类测点至少无法满足记录环境日变化的要求。

根据实验结果，发现测点的数量对实验结果有影响。浅层水温ZKGD3000型测点有5个，Szw-1型16个；深层水温ZKGD3000型测点有11个，Szw-1型9个。因这几类测点数量太少，在自动聚类时稳定性不是很好，其结果与实际情况有一定的偏差。如图3所示，尽管结果能将3类仪器差异明显地展示出来，但对传感器入水深度差异的区分效果不明显。

考虑到仪器类型、传感器入水深度、季节变化、地域影响等因素，文章仅对不同仪器和不同测项做差异性假设检验。实验表明，同一类仪器在不同测项下计算的分辨力有差异，同一测项的不同类仪器的测点分辨力也有差异。因此在设计实验时，按照测项进行分类。从计算结果看，每类仪器在不同测项下的分辨力有小幅度的差异，但不影响整体评价。季节变化、地域影响等受测点数量的限制，未做一一分类计算，但这些影响因素的差异性分析也是一个非常有意义的研究方向。

5 结论

从某一类仪器的全部观测数据中挖掘信息，用来分析估计仪器分辨力，是一种新方法的应用尝试。该方法对研究数据质量要求低，抗干扰性强，避免缺数插值、台阶处理等复杂的数据预处理过程，数据的统计模型简单，物理意义明确。相对于实验室对比观测的分辨力研究方法，更容易操作。实验结果表明，得出各类仪器的相对分辨力与其实际性能都较吻合，说明该研究方法有效、可行。