地下流体数据处理与产品加工软件系统设计

刘春国 李正媛 王建国 张彬 樊春燕 陈华静

1)中国地震台网中心，北京市西城区三里河南横街5号 100045

2)天津市地震局，天津 300201

3)中国地震局地壳应力研究所，北京 100085

0 前言

“十五”地震数字网络项目完成后，国家地下流体台网中心(以下简称地下流体中心)部署有前兆台网数据管理系统、前兆数据处理系统、前兆数据评价系统、地下流体台网管理与质量监控系统、地下流体专业网站等5个软件系统，构成了地下流体中心数据汇集、数据处理、质量监控与台网管理、数据服务的基础业务工作平台。

随着我国地下流体台网整合的完成(周克昌等，2013)，地下流体台网已经从建设期转入相对稳定的运行期，台网的产品产出能力日益引起关注，要求地下流体中心能够快速加工产出高质量的地下流体专业产品，以满足震情监测、地震科研以及其他防震减灾工作的迫切需要。现有的软件系统已经无法满足上述需求，在此背景下，配合“十一·五”背景场探测项目“前兆数据处理软件系统”的实施，对地下流体中心数据处理与产品加工软件系统的功能、处理模型和数据库进行了规划和设计。

系统设计的目标是满足背景场探测项目完成后地下流体中心所有地下流体测项常规产出和应急产品的数据处理、产品加工和制作的需求。处理的地下流体测项包括井水位(以下简称水位)、井(泉)水温(以下简称水温)、水(气)氡(以下简称氡)、水(气)汞(以下简称汞)、气体及离子等;根据测项的性质不同，这些测项可分为物理量测项和化学量测项两大类，物理量包括水位和水温等，化学量包括氡、汞、气体及离子等。

1 系统总体设计

1.1 设计基本原则

地下流体数据处理与产品加工软件系统是在背景场探测项目“前兆数据处理软件系统”总体设计框架下设计的。

背景场探测项目“前兆数据处理软件系统”的构成如图1所示，包括数据处理与产品加工和产品存储、交换与管理以及产品可视化展示、服务等3大分系统。其中数据处理与产品加工分系统包括地下流体、电磁和形变的数据处理与产品加工系统，实现3大学科台网数据处理与加工、产品制作功能;产品存储、交换与管理分系统实现产品数据管理和产品数据从5个学科台网中心(国家地下流体、地电、地磁、形变和重力台网中心)到国家前兆台网中心的自动汇集;产品可视化展示、服务分系统实现加工产品的可视化展示与对外服务功能。

地下流体数据处理与产品加工软件系统设计遵循统一规划、标准先进、模块化设计等原则。

图1 背景场前兆数据处理软件系统

1.2 系统构成与加工流程

系统的构成如图2所示，其中包括地下流体数据处理、数据产品加工和图形产品制作等3大软件模块和1个地下流体产品数据库。

系统加工流程为:通过数据处理模块，完成地下流体观测数据的质量检查、专业处理和质量分类工作，产出浅加工产品数据入产品数据库保存;数据产品加工模块从产品数据库中获取浅加工产品数据(经过数据处理、且质量等级满足加工需求)进行进一步加工处理，生成深加工产品数据;图形产品制作模块采用可视化展示技术对产品数据进行图形产品制作。

图2 地下流体数据处理与产品加工系统构成与加工基本流程

2 地下流体产品数据库设计

地下流体产品数据库用于存储地下流体数据与加工系统加工的各类产品数据及其辅助信息，是前兆产品数据库的组成部分。地下流体产品数据库中的数据表可分为3大类:

(1)前兆数据基本表

包括台站、测点、观测仪器、观测井泉等基础信息表，分钟值、整点值、日均值等数据表，以及质量评价信息表等。这些表已包括在现有的前兆数据库表结构(中国地震局，2012)中，不再设计，直接采用。

(2)新增学科通用表

包括月均值及年均值产品数据、时序图产品信息、空间图产品信息、背景图层信息、数据处理日志等数据表。这类各学科具有共性的数据表不需自行设计，采取与前兆产品数据库统一设计的表结构。

(3)新增流体自用表

这类表在现有前兆数据库表结构中没有，且只有地下流体的数据处理与加工才需要，因此需要自行设计。设计的基本原则是严格按照3NF范式进行规范化设计，考虑的主要因素包括数据的读写速度与数据冗余、数据表的可扩展性与动态参数的管理等。新增流体自用表如表1所示。

表1 表1 新增流体自用表

3 功能设计

系统共设计了数据处理、数据产品加工和图形产品制作等3个一级模块、16个子模块和37个功能小模块。

3.1 数据处理模块

数据处理模块主要完成数据产品加工前观测数据检查、问题数据处理，数据质量的审核与数据分类，以及产品加工前的一些必要的专业处理。通过数据处理后，可以生成带质量分类信息的观测数据产品和包括水位高程校正值、标定误差消除值等浅加工产品，并将其保存在产品数据库中。

数据处理模块包括数据检查与改正、质量审核和专业处理等3个子模块。

数据检查与改正子模块包含6个主要功能小模块:数据初检、错误数据处理、突跳数据处理、台阶数据处理、测量单位修正、水位值转换等，以实现对产品加工源数据的快速检查与地下流体常见的错误数据、突跳数据、台阶数据、非标准单位数据和未转换成水位值的水柱高度数据等的处理功能。

质量审核由质量审核参数管理和质量审核分类等2个功能小模块组成，以实现质量审核分类参数的动态管理、质量指标的自动计算和数据质量分级功能。

专业处理由标定误差分配、水位高程校正、缺数插值等3个功能小模块组成。其中标定误差分配功能是为了消除某些测项时序产品中出现的因仪器启用新K值的台阶而设计的;水位高程校正是将基于不同基准面的水位值归一到海平面，为水位的空间产品加工提供可能;缺数插值功能主要用于水位潮汐产品加工数据的处理。

3.2 数据产品加工

数据产品加工模块主要完成对经过数据处理且数据质量满足加工需求(质量为优秀或合格)的观测数据进行产品加工，以生成产品数据。数据产品加工模块由8个功能子模块组成:

(1)水位气压改正:对受气压影响明显的水位观测序列气压改正。

(2)水位逐日潮汐分析:对有固体潮效应的水位观测序列逐日计算日潮汐参数及其中的误差。

(3)均值与均方差:台网地下流体测项日均值、五日均值、旬均值、月均值、年均值及其均方差的自动计算。

(4)物理量变差:台网物理量年变差和月变差值的自动计算。

(5)化学量相对变差:台网化学量的月相对变差、年相对变差的自动计算。

(6)地震动态响应识别:台网中有地震响应变化的水位、水温、气氡、气汞测点的自动筛选识别，地震响应参数的提取。

(7)月潮汐因子变化:台网水位月潮汐因子变化值自动计算。

(8)产品数据归一化:对于需要在空间分布图上展示的产品数据，当数据变化范围较大时，可以采取该模块将数据归一到某一个范围。

3.3 图形产品制作

图形产品制作模块主要利用数据产品加工模块生成的产品数据制作图形产品。包括5个子模块及18个功能小模块:

(1)背景值分布图:采用台网地下流体主要测项的年(月)均值或均值归一化值产品数据，绘制背景值空间分布图。

(2)观测量变化分布图:采用台网物理量年(月)变差值及归一化年(月)变差值、化学量年(月)相对变差值产品数据，绘制各观测量变化空间分布图。

(3)同震响应分布图:采用台网水位、水温、气氡、气汞等同震响应参数(响应类型、响应幅度)，绘制同震相应分布图。

(4)水位潮汐因子变化分布图:采用台网水位潮汐因子变化值产品数据，绘制潮汐因子变化分布图。

(5)时序产品绘制:采用水位日潮汐参数、水位气压改正数据、各测项月变差值及年变差值等产品数据，绘制水位潮汐变化时序图、水位气压改正时序图、各测项月变差及年变差时序图等。

此外，(1)～(4)的分布图制作功能还包含以下功能:

可以在1幅图上，绘制1种测项或多种测项，用符号区分;通过颜色、符号大小显示物理量的大小;可标注地震;可叠加1:20区域水文地质图、断层分布图(绘制全国图时采用块体边界图替代断层分布图)，点击指定台站，能够显示台站物理量及其站点基础信息等。

时序产品绘制功能子模块还包含多条对比曲线、均值线、异常上下限的绘制，地震的标注，以及观测日志、站点基础信息的联机查询等功能。在水位潮汐变化时序图中提供了潮汐参数中误差(误差棒的形式)的标识功能。

4 处理模型设计

处理模型设计包括功能模块的输入、输出、处理流程与算法设计。基本原则是尽量采用在长期地下流体数据分析处理和运行评价中形成的成熟的数学模型，适当引入近年来地下流体监测预报领域使用的分析研究方法。

4.1 引用模型的改进

表2为本系统引用现有软件系统处理模型的基本情况。引用的处理模型只涉及到数据检查与改正(包括5个功能子模块)、质量审核和水位气压改正模块。这些模块并不能完全满足本系统的需求，必须加以改进。

(1)数据检查与改正

为了提高数据检查、改正的效率，在数据检查与改正模块中增加数据初检模块:采取自动扫描的方法，先筛选出疑似有问题的台项，再提交给错误数据处理、突跳数据处理、台阶数据处理、测量单位修正、水位值转换功能小模块进行处理与改正。筛选模型主要包括合理值识别、突变数据识别:

①合理值识别方法:根据各测项的合理值域范围来判定，超出范围的观测值被认为是错误数据。

②突变数据识别:采用一阶差分序列的超差识别方法，当差分值ΔXi满足

时，认为该点数据出现突变，有可能是突跳或台阶。式中σ为一阶差分序列的均方差，K为阈值。

(2)质量审核模块

为了对地下流体数据质量进行科学分类，我们在分析2010～2012年地下流体台网月评比的评价指标及得分标准、质量状况的基础上，筛选出内在质量、动态规律、动态稳定性等3个指标来衡量数据质量，数据质量分级模型为

式中，S动态规律、S稳定性、S内在质量分别为动态规律、稳定性、内在质量指标的月评比得分(均转换为100分制);α、β、γ为权重指标，α+β+γ=1，不同测项的 α、β、γ可以不同。根据S值的大小和分级标准，将观测数据质量分为合格、基本合格和不合格3个等级。

(3)水位气压改正处理模型

在水位气压改正之前，采用别尔采夫潮汐滤波方法(国家地震局科技监测司，1995)对水位和气压整点值数据进行了滤除潮汐引力的处理;采用对气压向前滑动的水位-气压回归分析法，求出回归分析相关系数最大且经过相关显著性检验的线性方程，进而求出滞后时间和气压系数。

4.2 新处理模型的提出

这里的“新”指的是在现有的前兆监测预报专业软件中尚未出现过。新处理模型主要包括标定误差分配、水位高程校正、水位逐日潮汐分析、月潮汐因子变化、物理量变差或化学量相对变差、地震动态响应识别和产品归一化等。

(1)标定误差分配

在目前的地下流体观测中经常出现仪器标定后启用新K值而出现台阶的现象。K值的变化过程很复杂，为了简化问题，认为台阶是误差线性累积造成的。将标定误差平均分配在两个标定周期间的观测数值中，对标定K值变化引起的台阶进行消除。算法如下

上一次标定结束后的观测序列{Xi，i=0，1，2…n}，Xn为此次标定前的最后一个观测数据，因本次标定造成观测值出现台阶，阶变量为ΔR，消除台阶的改正序列{Yi}可以用

进行计算。以某台站数字化气氡观测为例。2011年12月24日10点(TA)至2012年12月24日10点(TB)为上次标定结束后的观测时段(共计8785个数，n=8784)，2012年12月24日11点至2012年12月25日21点为本次标定时段，2012年12月25日22点恢复观测，测值与标定前测值出现台阶ΔR=5Bq/L，采用式(3)消除台阶。观测曲线、改正曲线如图3所示。从图3中可看出误差累积线与观测曲线叠加即为改正曲线，改正曲线显示台阶已经消除。

图3 某台站气氡观测时间序列标定误差分配示意图

(2)水位高程校正

水位观测按照是否有泄流，分为静水位观测和动水位观测。静水位一般是以井口为基准面，值越大表示水位上升;而动水位一般是以泄流口为基准面，值越大表示水位下降。

水位的空间产品要求空间分布的所有水位观测站点采用统一的基准面，一般选取海平面作为统一基准面。水位高程校正公式为

其中，H为水位校正值;G为观测井口或泄流口的海拔高程;h为水位值，当基准面为泄流口时h为正，为井口时h为负。

(3)水位逐日潮汐分析与月潮汐因子变化

水位逐日潮汐分析方法:在选定时段内，以一定的滑动步长(可以为1天、2天或任意给定天数)和窗长(一般为30天)，进行逐日水位固体潮汐分析。采用Venedikov调和分析方法(周克昌等，2011)，计算每次滑动后的水位M2波潮汐因子及其中误差、相位滞后及其中误差，将其作为某一天的潮汐参数。滑动结束后，计算得到水位M2波潮汐因子及其中误差、相位滞后及其中误差时间序列。

月潮汐因子变化自动计算具体算法如下:①首先根据上一年度评比计算结果，自动筛选出水位观测网中观测质量为合格以上，M2潮汐因子相对误差＜20%的观测站作为计算对象;②根据月评比计算的当月M2潮汐因子值以及上一年12个月的M2潮汐因子值，计算每个观测站上一年度水位月M2潮汐因子的平均值和均方差;③剔除超过3倍均方差的潮汐因子，剩余值的平均值作为背景值，第i观测站水位月M2潮汐因子变化值 Δi为

其中δi第i观测站当月M2潮汐因子;Vi为第i观测站的背景值。

(4)物理量变差或化学量相对变差

根据月均值和年均值的计算结果，自动计算质量等级合格以上的所有观测站物理量变差值和化学量相对变差值。

物理量变差包括月变差和年变差，月变差反映的是当月与上月相比的变化量，年变差反映的是当年与上年相比的变化量。计算公式为

其中X为水位高程校正值或水温观测值的月均值或年均值;ΔX为变差值。

化学量月变差反映的是化学量当月与上月相比的相对变化量，年变差反映的是当年与上年相比的相对变化量。相对变差计算公式为

其中，ΔX月=X当月－X上月，A月=(X当月+X上月)/2，ΔX年=X当年－X上年，A年=(X当年+X上年)/2，X为化学量月均值或年均值，ΔX为变差，A为平均值，CΔX为相对变差。

(5)地震动态响应识别

采用差分-均方差法与人机对话相结合的方式，快速识别有地震动态响应的台项、提取地震动态响应参数。具体步骤如下:

1)设原始观测序列 {Xi}，一阶差分序列{ΔXi}，对应的时间序列{ti}(一般不大于3天，应包含发震时刻td)，σ为{ΔXi}的均方差，阈值K(一般为3)，当某个观测站测项满足

初步认定其动态对地震有响应，响应开始时间为满足上述条件的最小时间ti。

2)自动筛选所有有地震动态响应的台项，并记录响应开始时间。

3)逐一浏览筛选出的台项的观测曲线，确认台项的地震响应判断和自动记录的响应开始时间是否准确，无误后再判断响应类型与形态，拉框自动计算响应幅度、提取响应结束时间。

(6)产品归一化

归一化方法在数据挖掘中也叫数据的正规化，将数据属性值从原来的取值区间映射到一个更适当的区间。常见的归一化方法有最小最大值、零均值和小数尺度法等3种。最小最大值法就是将数据归一至(0，1)(李爱国等，2012)。

考虑到产品数据有正、负值，我们对最小最大值法进行了改良，将正数归一到(0，1)，负数归一化到(0，－1)，计算方法如下:

设需要归一化的一组产品数据为{Xi}i=0，1，2…，n。

其中Xmin为{|xi|}的最小值，Xmax为{|xi|}的最大值;Yi为归一化值，sign(Xi)为产品数据的符号(正号或负号)。

4.3 处理模型测试与应用

我们重点对改进的处理模型、新提出的处理模型进行了测试。

采用C#开发语言，编制了基于C/S架构的测试软件。利用测试软件，连接地下流体中心地下流体数据库获取数据进行测试，测试结果显示处理模型算法正确、合理。

测试软件中的地震响应识别、水位气压改正、水位逐日潮汐分析等模块已经在地震应急和分析研究中使用。

采用水位气压改正模块对河南兰考豫11井2011年9月1～30日水位数据及其气压数据进行分析处理，计算结果表明该井水位与气压变化基本同步，气压系数为0.0075m/hPa，水位、气压观测曲线及其水位气压改正后曲线见图4所示。

图4 水位气压改正示例

在2013年4月芦山地震应急工作中，我们采用地震响应识别模块(图5)，对全国水位、水温观测网的同震效应进行了快速扫描，在1个小时内快速产出了水位同震分布图、水温同震分布图和同震时序图集，而在以前至少需要1天时间。

图5 地震响应识别示例

在2013年度地下流体学科会商水位资料处理中，应用水位逐日潮汐分析模块对某些台站的水位潮汐因子、相位滞后的变化情形进行分析研究。图6显示了重庆石柱鱼池水位在2012年～2013年8月潮汐因子的变化图像，图6中同时绘制了潮汐因子中误差棒，直观地给出了潮汐因子计算结果的可靠性。此外绘制的潮汐因子均值、2倍均方差上下限有助于潮汐因子异常的判断。

5 结束语

本系统设计是在对现有产品的梳理、新增产品的实用化研究、产品加工模型的调研等工作基础上完成，软件设计已经通过地下流体学科专家组的评审，并被背景场项目实施组采纳。本系统软件设计的创新点在于:

(1)本系统设计从数据检查、数据处理、数据浅加工、数据深加工到产品产出的数据加工流程，与实际的人工产品加工流程相吻合，从现有的软件中集成了与产品加工相关的部分功能，改变了以往数据处理、质量评价与产品加工软件分离的现状，确保加工生成的产品的正确性和合理性。

(2)建立了观测数据质量分级模型，为产品加工源数据的筛选、观测数据产品分级、对外提供数据共享服务奠定了基础。

(3)在多个模块中设计了自动搜索筛选及批量处理功能，为产品的快速产出提供了可能。引入归一化模型，使得物理量参差不齐、大小差异显著的空间图更加可视化。

图6 水位逐日潮汐分析示例

(4)除了满足地下流体中心常规产出与应急产出以外，还吸纳一些近年来被流体预测研究人员所关注的研究产品如水位日潮汐变化时序图、月M2潮汐因子变化分布图等，为地下流体中心开展分析研究提供便捷工具。

当然，地下流体台网产出产品水平还处在一个相对初级的阶段，物理意义明确又具有实际应用的产品还很少，一些产品的加工处理模型还不成熟，地下流体台网产品研发工作有待加强，产品加工处理模型有待于随着地下流体学科的发展进一步补充和完善。