地磁观测数据不仅具有学术研究意义,也是不可或缺的国家基础战略资源。但是,我国地磁台网中近1/3地磁台站正日益受到外界人工电磁干扰源的干扰(谢凡,2012),导致产出的地磁观测数据质量下降,影响数据的进一步分析和利用,必须进行数据预处理。
地磁秒数据预处理仍依靠“中国地震地球物理台网数据处理系统地磁学科版软件”完成,其去除逻辑是,使用软件的放大功能,用肉眼进行时间上的遍历,人工判断噪声,进而借助软件进行形态修正,实现噪声消除。该去噪方法自动化程度不高,且存在大量重复操作,特别是“尖峰状”干扰的去除,消耗大量人力。而且,噪声去除效果因人而异,难以保证数据具有稳定质量。因此,设计开发一套地磁数据噪声处理软件,自动化处理地磁观测数据,有针对性地识别干扰噪声,使得地磁数据预处理工作更加精准和高效,分析预报人员将更多精力投入地磁数据的深度挖掘和分析应用。
1 算法设计 1.1 数学形态学算法谢凡等(2011)认为,目前对地磁信号的常用手段仍是数字滤波器,尤以IIR滤波器、FIR滤波器为主,无法避免存在的时滞、相移等缺陷,并且当信号与噪声差异较小且难以区分,或者干扰信号的频带分布较宽时,数字滤波方法的应用受到限制。数学形态学(Mathematical Morphology,MM)是基于积分几何、随机集合论等数学理论建立起来的一种非线性信号处理方法。由数学家Matheron等(2002)共同创立,最早以图像形态特征作为研究对象,现已成功应用于图像处理、图像分析、计算机视觉以及电能扰动等工程实践领域。在噪声抑制领域,数学形态学应用广泛,主要集中在图像复原、心(脑)电信号的干扰抑制,其最大优点是不会增加局部信息的相关性。腐蚀、膨胀、开运算、闭运算是数学形态学中的4大基本运算。地磁秒数据是1 Hz采样的离散信号,用函数f(n)表示,描述实数域上的离散时间序列,其定义域
$ (f \ominus g(- n))(x) = \mathop {\inf }\limits_{y \in D \cap G} \{ f(y) - g(y - x)\} $ | (1) |
膨胀运算定义为
$ (f \oplus g(- n))(x) = \mathop {\sup }\limits_{y \in D \cap G} \{ f(y) + g(y - x)\} $ | (2) |
开运算定义为
$ (f \circ g)(n) = (f \ominus g \oplus g)(n) $ | (3) |
闭运算定义为
$ (f \bullet g)(n) = (f \oplus g \ominus g)(n) $ | (4) |
一般,开运算将信号中的孤立部分分离,抑制信号中的正脉冲噪声;闭运算的作用是补缺和内部连通,用于抑制信号的负脉冲噪声。
1.2 滤波器设计基于数学形态学去除噪声的基本思想是,利用一个结构元素去探测图像,达到对图像中噪声分析和认识,进而消除噪点的目的。
结构元素对目标信号进行变化和匹配,是识别并抑制不相干信号的“探针”。因此,结构元素选取是滤波器设计的关键,选取合适的结构元素能够提高滤波器性能。结构元素一般有脉冲型、正(余)弦型、扁平型、矩形等。为了消除脉冲状噪声,优化选取脉冲型结构元素,并基于数学形态学的4大基本运算,构建交替混合滤波器,则
$ {f_{{\rm{ah}}}}(n) = \left[ {{f_{{\rm{oc}}}}(n) + {f_{{\rm{co}}}}(n)} \right]/2 $ | (5) |
其中,
软件设计是指,对软件系统功能进行模块划分,主要包括类的设计和界面规划,并结合实际需要,提供相应参数接口。软件项目设计开发的第1个阶段是需求分析,即理解用户需求,就软件功能和用户达成一致,估计软件风险,评估项目代价,最终形成开发计划。软件需满足以下需求:①从数据库读取地磁数据。重庆市地球物理台网台网数据库依照中国地震地球物理台网数据库结构规范进行设计(周克昌等,2010),本软件需要直接连接重庆市地球物理台网数据库,并读取相应地磁数据;②数据可视化。能够清晰、简洁地显示地磁数据,通过查看,用户能大致确定噪点所处位置;③自动化去除噪声数据。基于数学形态学算法,自动计算并消除地磁数据噪声,利用简单参数控制消除效果的强度。
根据需求分析,设计软件开发技术路线,见图 1。
针对用户需求,设计软件系统的3大功能模块,包括数据库读写模块、地磁数据可视化模块、噪声消除计算模块,均以地磁数据可视化模块的人机交互界面为主界面,并在主界面基础上,对用户交互操作进行响应。
(1) 数据库读写模块。该模块是地磁数据的入口和出口。数据库读写由DBReadWrite类实现,通过连接重庆市地球物理台网数据库,获取QZDATA.QZ_312_DYS_02(原始数据)表和QZDATA.QZ_312_DYU_02(预处理数据)表,字段为OBSVALUE的地磁秒数据,将Oracle数据库特有的CLOB类型转化为String数组类型,同时借助Qt框架中执行SQL语句的方法,进行数据库读写操作。DBReadWrite类主要实现querySQL()、saveSQL()、updateSQL()等方法对待执行SQL语句的更新。
(2) 地磁数据可视化模块。DisplayView类是3大模块核心,主要负责从数据库读写模块中获取数据(通过调用readDataFromDB()函数),对数据进行可视化,同时以直观形式为地磁数据的噪声消除提供参考。数据可视化主要包括地磁数据可视化(包括原始数据和预处理数据)、检测的噪声点可视化、去除噪声后的地磁数据可视化。通过导入第三方库QCustomPlot,借助其实例化对象plotWidget,进一步实现对地磁数据显示曲线的无级别放大缩小、拖曳,以及自适应的横纵坐标刻度标尺和刻度值的显示(通过调用drawGraph()、drawXais()、drawYaxis()等函数实现)。
(3) 噪声消除计算模块。算法类FilterAlgorithm实现数学形态学滤波器的算法。考虑到噪声消除计算对效率和潜在可移植性的需求,噪声消除计算模块由C语言编写,并预留接口函数。FilterAlgorithm类实现了数学形态学的4大基本运算:腐蚀运算[corrode(Qvector<double>dat, Qvector<double>g, int op)]、膨胀运算[expand(Qvector<double>dat, Qvector<double>g, int op)]、闭运算[closing(Qvector<double>dat, Qvector<double>g, int op)]、开运算[opening(Qvector<double>dat, Qvector<double>g, int op)],借助形参一致的4大基本运算函数,进一步构建数学形态学的滤波器filter(Qvector<double>dat, Qvector<double>g, int op)。
基于UML思想(Unified Modeling Language,统一建模语言)规划控制类图(范晓平,2005),各模块之间关系如下:DBReadWrite类和DisplayView类是单向关联(Association)关系,二者通信实现噪声滤波前后地磁数据的数据库存取。MainWindow类是系统控件主窗口的实现类,属于软件控制层,与属于视图层的DisplayView类组合(Composition),完成用户界面绘制。FilterAlgorithm类聚合(Aggregation)于主窗口类MainWindow,且通过接口函数通信,降低模块间的耦合性。以DisplayView类为核心,4个类相互协作,协同实现软件的主要功能。软件UML控制类见图 2。
本软件实现对重庆地球物理台网数据库中地磁数据的浏览功能,可筛选日期、台站、测项、测点(图 3)。可视化功能的亮点是,实现地磁数据曲线的无级别缩放,并通过优化代码,使缩放过程无任何卡顿。
本软件能够对数据进行一阶差分情况下的超差检测。在数据正常读取显示前提下,点击功能区“超差检测”按钮,对超差数据点进行标红醒目处理,超差处通常为尖峰噪声的存在之处。在超差检测完成基础上,点击功能区“去除尖峰”按钮,将执行数学形态学滤波算法驱动下的尖峰做剔除处理,完成后数据以绿色线条形式绘制,见图 4。点击“保存到数据库”按钮,即可将噪声处理完成后的数据一键保存至地球物理台网系统数据库。
在重庆地区实际地磁观测中,地磁数据受环境、人为等因素的影响(城市轨道交通、超高压直流输电、仪器背景噪声等)。由时间维度和噪声数量可知,绝大部分干扰属于仪器背景噪声产生的“尖峰状”干扰,其他类型的干扰有限。本软件去除噪声的关键核心是,既要保证“尖峰状”噪声的去除效果,又不破坏地磁数据的有用信息(特别是磁暴、磁扰等地磁形态变化剧烈的空间事件,数据应该被完整保留)。因此,软件使用数学形态学算法优化对应的结构元素,并采用一阶差分和形态学滤波器的融合算法,在精准识别“尖峰状”噪声的同时,实现对其他信息的无损伤保留。
以万州天星2018年8月16日Z分量原始数据为例,使用阈值设定为0.5的一阶差分检测方法,将超限数值点用红色标记标出,通过数学形态学滤波器进行噪声消除(消除效果取决于数学形态学滤波器和构建的结构元素),滤波前后数据曲线见图 5,可见消除结果准确去掉“尖峰状”噪点,且未产生相位畸变和时移,完整保留噪声以外的数据信息。
另一个压力测试实例为涪陵江东2017年8月31日Z分量原始数据,噪声处理前后曲线见图 6。由图 6可见,2017年8月31日5时37分起,Z分量记录到K指数最大为7的急始磁暴,期间Z分量数据形态变化剧烈,出现多个台阶和受磁暴影响而产生的复杂扰动形态,经噪声消除,完整保留地磁磁暴形态信息。
利用Qt环境,设计开发地磁数据噪声处理软件。该软件可在Windows系列系统运行,并通过LGPL授权发布。在软件设计过程中,基于UML标准化建模语言,对软件的结构层次进行规划,遵循MVC设计原则,将Model(模型层)、View(视图层)、Control(控制层)分离,保证功能模块的高内聚、低耦合,使得软件足够健壮并拥有较长的生命周期。
地磁数据噪声处理软件能够直接连接重庆市地球物理台网数据库,将地磁秒数据绘制为可视化图形;利用数学形态学算法和优化选取的结构元素,消除地磁数据噪声并存入数据库。该软件界面清晰,操作简便,节省了以往半自动化去除噪声所需人力。软件使用的数学形态学滤波器方法具有积极的理论意义,形成的软件实体,具有较为广泛的应用前景。
范晓平. UML建模实例详解[M]. 北京: 清华大学出版社, 2005. | |
谢凡, 滕云田, 胡星星. 数学形态滤波在地磁数据干扰抑制中的应用[J]. 地球物理学进展, 2011, 26(1): 147-156. DOI:10.3969/j.issn.1004-2903.2011.01.016 | |
谢凡. 地磁观测中干扰抑制方法的发展及展望[J]. 地球物理学进展, 2012, 27(3): 967-976. | |
周克昌, 蒋春花, 纪寿文, 等. 地震前兆数据库系统设计[J]. 地震, 2010, 30(2): 143-151. DOI:10.3969/j.issn.1000-3274.2010.02.016 | |
Maragos P, Schafer R W. Morphological filters. I. Their set-theoretic analysis and relations to linear shift-invariant filters[J]. IEEE Transactions on Acoustics Speech & Signal Processing, 1987, 35(8): 1153-1169. | |
Matheron G, Serra J. The birth of mathematical morphology[C]//Proc. 6th International Symposium on Mathematical Morphology. Sydney, Australia: CSIRO Publishing, 2002. |