留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

连续重力观测异常模式的多分辨率识别算法

刘子维 张晓彤 张锐 江颖 韦进 张坤

刘子维, 张晓彤, 张锐, 江颖, 韦进, 张坤. 连续重力观测异常模式的多分辨率识别算法[J]. 武汉大学学报 ● 信息科学版, 2018, 43(6): 840-846, 942. doi: 10.13203/j.whugis20160173
引用本文: 刘子维, 张晓彤, 张锐, 江颖, 韦进, 张坤. 连续重力观测异常模式的多分辨率识别算法[J]. 武汉大学学报 ● 信息科学版, 2018, 43(6): 840-846, 942. doi: 10.13203/j.whugis20160173
LIU Ziwei, ZHANG Xiaotong, ZHANG Rui, JIANG Ying, WEI Jin, ZHANG Kun. Multiresolution Recognition Algorithm for Abnormal Patterns of Continuous Gravity Observation[J]. Geomatics and Information Science of Wuhan University, 2018, 43(6): 840-846, 942. doi: 10.13203/j.whugis20160173
Citation: LIU Ziwei, ZHANG Xiaotong, ZHANG Rui, JIANG Ying, WEI Jin, ZHANG Kun. Multiresolution Recognition Algorithm for Abnormal Patterns of Continuous Gravity Observation[J]. Geomatics and Information Science of Wuhan University, 2018, 43(6): 840-846, 942. doi: 10.13203/j.whugis20160173

连续重力观测异常模式的多分辨率识别算法

doi: 10.13203/j.whugis20160173
基金项目: 

国家重大科学仪器设备开发专项 2012YQ10022506

国家自然科学基金 41374088

国家自然科学基金 41304018

国家自然科学基金 41404065

国家自然科学基金 41404064

地震科技星火计划 XH15028

详细信息
    作者简介:

    刘子维, 博士, 副研究员, 主要研究方向为重力观测技术和数据分析。lzw@eqhb.gov.cn

    通讯作者: 张锐, 高级工程师。rzhang@neis.gov.cn
  • 中图分类号: P223.7

Multiresolution Recognition Algorithm for Abnormal Patterns of Continuous Gravity Observation

Funds: 

The National Key Scientific Instrument and Equipment Development Projects 2012YQ10022506

the National Natural Science Foundation of China 41374088

the National Natural Science Foundation of China 41304018

the National Natural Science Foundation of China 41404065

the National Natural Science Foundation of China 41404064

Science for Earthquake Resilience of China XH15028

More Information
    Author Bio:

    LIU Ziwei, PhD, associate researcher, specializes in the gravity observation technique and gravity data analysis. E-mail: lzw@eqhb.gov.cn

    Corresponding author: ZHANG Rui, senior engineer. E-mail:rzhang@neis.gov.cn
  • 摘要: 在小波多孔算法的基础上,提出了一种综合信号频率信息和幅值信息的连续重力观测数据多分辨率异常模式识别算法,利用小波分解得到高频区域的能量作为频率指标,与幅值相结合,对信号及其多孔小波分解结果进行多分辨率异常模式识别。利用模拟数据和实际超导重力观测数据对算法的有效性进行了验证,结果表明,该算法能够准确地在带有噪声的信号中识别模拟数据的异常模式,可应用于连续重力观测台网数据分析处理,对于提升台网观测数据质量以及地震预测等实际应用都具有重要意义。用此方法分析拉萨和武汉的3台超导重力仪2015-04-25尼泊尔地震前一天的秒采样数据后,得到一段27 min的在频率指标上有超过90%相似性的异常模式,这一结果的更深层次物理解释仍需要进一步研究。
  • 图  1  算法流程图

    Figure  1.  Flow Chart of Algorithm

    图  2  模拟信号

    Figure  2.  Simulated Testing Data

    图  3  模拟信号搜寻结果

    Figure  3.  Search Results of Simulated Data

    图  4  合成信号多孔小波分解结果

    Figure  4.  Wavelet Decomposition Results of Synthetic Signals

    图  5  合成信号多孔小波分解第2层和第3层搜索结果

    Figure  5.  The Second Layer and the Third Layer Search Results of Synthetic Signals

    图  6  拼接数据

    Figure  6.  Splicing Data

    图  7  拼接数据搜索结果

    Figure  7.  Part of Search Results of Splicing Data

    表  1  不同模拟信号搜寻结果/%

    Table  1.   Search Results of Simulated Testing Data with Different SNR/%

    信噪比 周期/s 总比例
    16 32 64 128
    正弦 衰减 正弦 衰减 正弦 衰减 正弦
    1 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
    2 0.00 0.00 0.00 0.00 81.25 29.69 0.00 20.17
    4 0.00 32.81 100.00 60.94 96.09 50.78 43.55 58.66
    8 87.50 56.25 100.00 91.41 100.00 63.28 93.75 87.71
    下载: 导出CSV

    表  2  模拟信号分层搜寻结果/%

    Table  2.   Search Results at Different Decomposition Result of Simulated Testing Data/%

    周期/s 总比例
    16 32 64 128
    正弦 衰减 正弦 衰减 正弦 衰减 正弦
    原始信号 0.00 0.00 0.00 0.00 81.25 29.69 0.00 20.17
    第1层 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
    第2层 60.94 90.63 100.00 100.00 99.22 72.66 58.59 77.63
    第3层 98.44 40.63 100.00 100.00 100.00 100.00 98.83 96.80
    下载: 导出CSV

    表  3  拼接数据相似对

    Table  3.   Search Results of Splicing Data

    搜索结果 搜索结果
    1 262~297 3 096~3 131
    2 387~423 2 391~2 427
    3 749~780 2 879~2 910
    4 812~844 4 278~4 310
    5 950~982 3 385~3 417
    6 950~981 3 386~3 417
    7 951~983 3 385~3 417
    8 1 006~1 044 4 191~4 229
    9 1 137~1 168 3 001~3 032
    10 1 408~1 444 3 381~3 417
    11 1 408~1 443 3 382~3 417
    12 1 409~1 445 3 381~3 417
    13 1 409~1 442 3 384~3 417
    14 1 412~1 446 3 383~3 417
    15 1 464~1 498 2 904~2 938
    16 1 567~1 603 4 181~4 217
    17 1 851~1 885 3 291~3 325
    18 2 387~2 419 3 827~3 859
    下载: 导出CSV

    表  4  拼接数据分层搜索结果

    Table  4.   Search Results at Different Decomposition Result of Splicing Data

    搜索结果 搜索结果
    第1层 640~6741
    442~1 474
    2 608~2 6423
    528~3 560
    368~416 2 664~ 2712
    406~464 3 306~3 364
    第2层 526~575 3 426~3 475
    940~991 3 156~3 207
    1 131~1 178 2 990~3 037
    1 312~1 359 3 260~3 307
    26~75 3 040~3 089
    47~94 1 960~2 007
    232~282 3 293~3 343
    233~280 3 293~3 340
    514~568 3 085~3 139
    第3层 595~648 3 593~3 646
    672~725 2 799~2 852
    762~811 3 432~3 481
    1 023~1 081 2 974~3 032
    1 043~1 093 4 229~4 279
    1 847~1 913 3 593~3 659
    1 877~1 928 3 317~3 368
    17~83 1 963~2 029
    19~83 1 964~2 028
    19~97 2 899~2 977
    21~98 2 902~2 979
    65~131 1 957~2 023
    第4层 1 282~1 349 3 406~3 473
    1 400~1 466 3 136~3 202
    1 957~2 025 4 045~4 113
    1 958~2 034 4 045~4 121
    2 390~2 461 4 122~4 193
    2 616~2 679 4 116~4 179
    530~660 3 735~3 865
    1 340~1 504 3 033~3 197
    第5层 1 340~1 502 3 034~3 196
    1 340~1 485 3 035~3 180
    1 545~1 695 2 990~3 140
    1 545~1 695 2 991~3 141
    下载: 导出CSV
  • [1] 刘子维, 李辉, 徐中华, 等. gPhone重力仪数据采集系统性能的改进[J].大地测量与地球动力学, 2010, 30(B11):102-104 http://mall.cnki.net/magazine/Article/DKXB2010S2024.htm

    Liu Ziwei, Li Hui, Xu Zhonghua, et al. Improvement in Performance of Data Acquisition System of gPhone Gravitymeters[J]. Journal of Geodesy and Geodynamics, 2010, 30(B11):102-104 http://mall.cnki.net/magazine/Article/DKXB2010S2024.htm
    [2] Sun Heping, Chen Xiaodong, Liu Ming, et al. Analysis and Comparison of the Tidal Gravity Observations Obtained with LCR-ET20 Spring Gravimeter[J]. Acta Seismologica Sinica, 2002, 15(5):533-539 doi:  10.1007/s11589-002-0022-1
    [3] Banka D, Crossley D. Noise Levels of Superconducting Gravimeters at Seismic Frequencies[J]. Geophysical Journal International, 1999, 139(1):87-97 doi:  10.1046/j.1365-246X.1999.00913.x
    [4] Tang K W, An D. A New Framework for Pattern Recognition of Time-series Data[D]. New York: Stony Brook University, 2004
    [5] Duda R O, Hart P E, Stork D G. Pattern Classification[M]. New York:Wiley, 1973
    [6] Duda R O, Hart P E, Stork D G. Pattern Classification[M]. 2nd ed. Hoboken, NJ:Wiley, 2001
    [7] Kaufman L, Rousseeuw P J. Finding Groups in Data:An Introduction to Cluster Analysis[M]. Hoboken, NJ:Wiley, 1990
    [8] 杨艳林, 叶枫, 吕鑫, 等.一种基于DTW聚类的水文时间序列相似性挖掘方法[J].计算机科学, 2016, 43(2):245-249 doi:  10.11896/j.issn.1002-137X.2016.02.051

    Yang Yanlin, Ye Feng, Lv Xin, et al. DTW Clustering-based Similarity Mining Method for Hydrological Time Series[J]. Computer Science, 2016, 43(2):245-249 doi:  10.11896/j.issn.1002-137X.2016.02.051
    [9] Carpenter G A, Grossberg S. ART 2:Stable Self-organization of Pattern Recognition Codes for Analog Input Patterns[J]. Applied Optics, 1987, 26(23):4919-4930 doi:  10.1364/AO.26.004919
    [10] Vargas F, Lettnin D, Felippettod C M C, et al. Electrocardiogram Pattern Recognition by Means of MLP Network and PCA: A Case Study on Equal Amount of Input Signal Types[C]. SBRN 2002 Proceedings VⅡ Brazilian Symposium on IEEE, Pernambuco, 2002
    [11] Vedam H, Venkatasubramanian V. A Wavelet Theory-based Adaptive Trend Analysis System for Process Monitoring and Diagnosis[C]. American Control Conference, IEEE, Albuquerque, 1997
    [12] 郑诚, 蔡庆生.一种多尺度的时间序列相似模式匹配算法[J].小型微型计算机系统, 2003, 24(3):546-549 http://www.cqvip.com/QK/95659X/2003003/7473560.html

    Zheng Cheng, Cai Qingsheng. A Muti-scale Similar Pattern Match Approach for Times Series Databases[J]. Mini-Micro System, 2003, 24(3):546-549 http://www.cqvip.com/QK/95659X/2003003/7473560.html
    [13] 张海勤, 蔡庆生.基于小波变换的时间序列相似模式匹配[J].计算机学报, 2003, 26(3):373-377 http://www.cqvip.com/qk/90818x/2003003/7788269.html

    Zhang Haiqin, Cai Qingsheng. Time Series Similar Pattern Matching Based on Wavelet Transform[J]. Chinese Journal of Computers, 2003, 26(3):373-377 http://www.cqvip.com/qk/90818x/2003003/7788269.html
    [14] Carpenter G A, Grossberg S, Rosen D B. Fuzzy ART:Fast Stable Learning and Categorization of Analog Patterns by an Adaptive Resonance System[J]. Neural Networks, 1991, 4(6):759-771 doi:  10.1016/0893-6080(91)90056-B
    [15] Wong S K M. The Relational Structure of Belief Networks[J]. Journal of Intelligent Information Systems, 2001, 16(2):117-148 doi:  10.1023/A:1011237717300
    [16] Looney C G. Pattern Recognition Using Neural Networks:Theory and Algorithms for Engineers and Scientists[M]. New York:Oxford University Press, 1997
    [17] Dutilleux P. An Implementation of the "algorithme à trous" to Compute the Wavelet Transform[M]. Berlin Heidelberg:Springer, 1989
    [18] Percival D B, Walden A T. Wavelet Methods for Time Series Analysis[M]. Cambridge Cambridge University Press, 2002
    [19] 牟力, 陈召曦.重力资料多尺度分析最优小波基的选择[J].物探与化探, 2015, 39(5):1013-1019 http://ch.whu.edu.cn/CN/abstract/abstract1759.shtml

    Mou L, Chen Z X.The Optimal Choice of Wavelet Bases in Gravity Data Multi-scale Analysis[J]. Geophysical and Geochemical Exploration, 2015, 39(5):1013-1019 http://ch.whu.edu.cn/CN/abstract/abstract1759.shtml
    [20] 徐华君, 柳林涛, 许厚泽, 等.重力极潮的小波分析[J].武汉大学学报·信息科学版, 2008, 33(11):1114-1117 http://www.cqvip.com/QK/92848A/200811/28650110.html

    Xu Huajun, Liu Lintao, Xu Houze, et al. Wavelet Approach to Study Gravity Pole Tide[J]. Geomatics and Information Science of Wuhan University, 2008, 33(11):1114-1117 http://www.cqvip.com/QK/92848A/200811/28650110.html
    [21] 许闯, 罗志才, 林旭, 等.重力固体潮观测数据的自动化预处理[J].武汉大学学报·信息科学版, 2013, 38(2):157-161 http://ch.whu.edu.cn/CN/abstract/abstract6095.shtml

    Xu Chuang, Luo Zhicai, Lin Xu, et al. Automatic Preprocessing of Tidal Gravity Observation Data[J]. Geomatics and Information Science of Wuhan University, 2013, 38(2):157-161 http://ch.whu.edu.cn/CN/abstract/abstract6095.shtml
    [22] 徐华君, 柳林涛, 许厚泽, 等.利用小波分析重力的长期变化[J].地球物理学报, 2008, 51(3):735-742 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=dqwlxb200803014

    Xu Huajun, Liu Lintao, Xu Houze, et al. Wavelet Approach to Study the Secular Gravity Variation[J]. Chinese Journal of Geophysics, 2008, 51(3):735-742 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=dqwlxb200803014
  • [1] 姚彦吉, 柳林涛, 王国成, 沈聪, 彭钊, 邵永谦.  地震事件自动识别的标准时频变换方法 . 武汉大学学报 ● 信息科学版, 2022, 47(5): 780-788. doi: 10.13203/j.whugis20190432
    [2] 张立福, 王飒, 刘华亮, 林昱坤, 王晋年, 朱曼, 高了然, 童庆禧.  从光谱到时谱——遥感时间序列变化检测研究进展 . 武汉大学学报 ● 信息科学版, 2021, 46(4): 451-468. doi: 10.13203/j.whugis20200666
    [3] 戴海亮, 孙付平, 姜卫平, 肖凯, 朱新慧, 刘婧.  小波多尺度分解和奇异谱分析在GNSS站坐标时间序列分析中的应用 . 武汉大学学报 ● 信息科学版, 2021, 46(3): 371-380. doi: 10.13203/j.whugis20190107
    [4] 吴浩, 卢楠, 邹进贵, 郭世泰.  GNSS变形监测时间序列的改进型3σ粗差探测方法 . 武汉大学学报 ● 信息科学版, 2019, 44(9): 1282-1288. doi: 10.13203/j.whugis20170338
    [5] 郝国成, 李飞, 白雨晓, 王巍.  基于NDSST的非平稳信号时频分析算法 . 武汉大学学报 ● 信息科学版, 2019, 44(6): 941-948. doi: 10.13203/j.whugis20170271
    [6] 程子桉, 庞小平, 赵羲, 季青.  1978~2014南极海冰边缘线长度时间序列变化 . 武汉大学学报 ● 信息科学版, 2016, 41(11): 1463-1468. doi: 10.13203/j.whugis20150263
    [7] 吴文豪, 王明洲, 李沙, 侯爱羚.  滑动聚束模式SAR影像干涉处理方法 . 武汉大学学报 ● 信息科学版, 2015, 40(12): 1588-1593. doi: 10.13203/j.whugis20140012
    [8] 裴媛媛, 廖明生, 王寒梅.  时间序列SAR影像监测堤坝形变研究 . 武汉大学学报 ● 信息科学版, 2013, 38(3): 266-269.
    [9] 魏二虎, 李智强, 龚光裕, 张帅.  极移时间序列模型的拟合与预测 . 武汉大学学报 ● 信息科学版, 2013, 38(12): 1420-1424.
    [10] 罗三明, 杨国华, 李陶, 袁油新.  PSInSAR方法探测意大利拉奎拉地震形变过程分析 . 武汉大学学报 ● 信息科学版, 2012, 37(5): 602-605.
    [11] 张燕, 吴云, 段维波, 吕品姬.  GPS长趋势变化与大地震的关系 . 武汉大学学报 ● 信息科学版, 2012, 37(6): 675-678.
    [12] 何敏, 何秀凤.  利用时间序列干涉图叠加法监测江苏盐城地区地面沉降 . 武汉大学学报 ● 信息科学版, 2011, 36(12): 1461-1465.
    [13] 归庆明, 李涛, 衡广辉.  时间序列异常值探测的Bayes方法及其在电离层VTEC数据处理中的应用 . 武汉大学学报 ● 信息科学版, 2011, 36(7): 802-806.
    [14] 张鹏, 蒋志浩, 秘金钟, 党亚民.  我国GPS跟踪站数据处理与时间序列特征分析 . 武汉大学学报 ● 信息科学版, 2007, 32(3): 251-254.
    [15] 刘飞, 王新洲, 张鹏林, 余旭.  基于Lebesgue积分理论的时间序列信号中的背景重建 . 武汉大学学报 ● 信息科学版, 2006, 31(5): 462-465.
    [16] 桂任舟, 杨子杰.  基于信号分解的时频分析方法在高频地波雷达目标监测中的应用研究 . 武汉大学学报 ● 信息科学版, 2006, 31(7): 653-656.
    [17] 杜国明, 龚健雅, 朱家松.  时间序列中反向查询算法的研究 . 武汉大学学报 ● 信息科学版, 2004, 29(1): 52-54,62.
    [18] 王智均, 李德仁, 李清泉.  Wallis变换在小波影像融合中的应用 . 武汉大学学报 ● 信息科学版, 2000, 25(4): 338-342.
    [19] 程建权, 黄经南.  一种基于时序数据的动态聚类分析方法 . 武汉大学学报 ● 信息科学版, 1998, 23(3): 194-198.
    [20] 徐培亮.  应用时间序列方法作大坝变形预报 . 武汉大学学报 ● 信息科学版, 1988, 13(3): 23-31.
  • 加载中
图(7) / 表(4)
计量
  • 文章访问数:  1179
  • HTML全文浏览量:  97
  • PDF下载量:  229
  • 被引次数: 0
出版历程
  • 收稿日期:  2016-07-27
  • 刊出日期:  2018-06-05

连续重力观测异常模式的多分辨率识别算法

doi: 10.13203/j.whugis20160173
    基金项目:

    国家重大科学仪器设备开发专项 2012YQ10022506

    国家自然科学基金 41374088

    国家自然科学基金 41304018

    国家自然科学基金 41404065

    国家自然科学基金 41404064

    地震科技星火计划 XH15028

    作者简介:

    刘子维, 博士, 副研究员, 主要研究方向为重力观测技术和数据分析。lzw@eqhb.gov.cn

    通讯作者: 张锐, 高级工程师。rzhang@neis.gov.cn
  • 中图分类号: P223.7

摘要: 在小波多孔算法的基础上,提出了一种综合信号频率信息和幅值信息的连续重力观测数据多分辨率异常模式识别算法,利用小波分解得到高频区域的能量作为频率指标,与幅值相结合,对信号及其多孔小波分解结果进行多分辨率异常模式识别。利用模拟数据和实际超导重力观测数据对算法的有效性进行了验证,结果表明,该算法能够准确地在带有噪声的信号中识别模拟数据的异常模式,可应用于连续重力观测台网数据分析处理,对于提升台网观测数据质量以及地震预测等实际应用都具有重要意义。用此方法分析拉萨和武汉的3台超导重力仪2015-04-25尼泊尔地震前一天的秒采样数据后,得到一段27 min的在频率指标上有超过90%相似性的异常模式,这一结果的更深层次物理解释仍需要进一步研究。

English Abstract

刘子维, 张晓彤, 张锐, 江颖, 韦进, 张坤. 连续重力观测异常模式的多分辨率识别算法[J]. 武汉大学学报 ● 信息科学版, 2018, 43(6): 840-846, 942. doi: 10.13203/j.whugis20160173
引用本文: 刘子维, 张晓彤, 张锐, 江颖, 韦进, 张坤. 连续重力观测异常模式的多分辨率识别算法[J]. 武汉大学学报 ● 信息科学版, 2018, 43(6): 840-846, 942. doi: 10.13203/j.whugis20160173
LIU Ziwei, ZHANG Xiaotong, ZHANG Rui, JIANG Ying, WEI Jin, ZHANG Kun. Multiresolution Recognition Algorithm for Abnormal Patterns of Continuous Gravity Observation[J]. Geomatics and Information Science of Wuhan University, 2018, 43(6): 840-846, 942. doi: 10.13203/j.whugis20160173
Citation: LIU Ziwei, ZHANG Xiaotong, ZHANG Rui, JIANG Ying, WEI Jin, ZHANG Kun. Multiresolution Recognition Algorithm for Abnormal Patterns of Continuous Gravity Observation[J]. Geomatics and Information Science of Wuhan University, 2018, 43(6): 840-846, 942. doi: 10.13203/j.whugis20160173
  • 经过多年的建设,中国已建成一个覆盖中国大陆,集数据自动汇集与数据库式管理、仪器远程监控和数据图形化分析处理为一体的连续重力观测台网[1],观测仪器包括GWR超导、PET/gPhone、DZW、GS15和TRG-1在内共计80余台。产出的观测数据不仅可用于活动断层、地壳运动、地震监测和海潮等地球动力学现象的深入研究[2],并且随着仪器采样率的提高,观测数据中高频信息的丰富,使得区域性的地球动力学现象和地球动力变化的微小过程研究成为可能[3]

    连续重力观测数据是典型的时间序列数据,如何利用有效的时间序列异常模式识别算法在观测数据中自动、精确地定位异常,对于提升观测数据质量以及地震预测等实际应用都具有重要意义。多年来针对时间序列异常模式识别,国内外学者开展了广泛的研究,包括聚类分析[4-8]、自组织映射[9]、主成分分析[10]、小波变换[11-13]、神经网络[14-16]等多种算法被广泛应用。聚类分析方法将数据按一定准则进行划分,使得不同类数据差异性最大化,其最大优点在于直观性强;主成分分析方法用包含原始信号大部分信息的几个变量来代替原始多个变量,且计算较为规范,容易实现;小波分析能够较好地反映信号的瞬态或变异点;神经网络具有并行处理、自适应组织、联想记忆等特点。但这些研究大都只在时间域上考虑了幅值的相似性,而没有在频域上考虑信号频谱特征的相似性,另外大多数算法需要首先提供参考异常模式。本文提出一种综合信号频率和幅值信息的多分辨率连续重力观测数据异常模式识别算法。

    • 考虑两个成分完全一致但在每一时刻对应不同振幅的信号,如果仅通过幅值来比较,它们将被区分为不同的信号,相反对于两个成分不同但幅值几乎相同的信号,反而归为相似的信号。为了克服这两种情形,首先利用多孔算法对原始信号进行快速分解,然后针对原始信号和每一层分解结果,综合观测信号的局部频率信息(通过Haar小波分解计算)和幅值信息对相似信号片段进行匹配,最终实现多分辨率观测数据异常模式的自动识别。

    • 多孔算法是一种利用Mallat算法结构计算小波变换的快速算法,通过在滤波器中补零产生空洞得到的滤波器组对信号进行卷积计算,取消传统小波变换中的二次抽样以保留信号中的大量细节,以数据的冗余实现了变换的平移不变性,且具备较强的抗噪能力[17-18]

      假设h(n)为低通滤波器,因此信号a(n)的第j层非抽样连续近似可表示为:

      $$ {a_j}\left( n \right) = \sum\limits_{}^k {h\left( k \right){a_{j - 1}}\left( {n - {2^{j - 1}}k} \right)} ,j = 1,2 \cdots J $$ (1)

      通过高通滤波器g(n)可得到细节系数计算式为:

      $$ {d_j}\left( n \right) = \sum\limits_{}^k {g\left( n \right){a_{j - 1}}\left( {n - {2^{j - 1}}k} \right)} ,j = 1,2 \cdots J $$ (2)

      需要使用两个对偶滤波器h(n)和g(n)来精确地重构原始信号,且需满足正交镜像滤波器的条件。

      $$ \bar h\left( n \right) * h\left( n \right) + \bar g\left( n \right) * g\left( n \right) = \delta \left( n \right) $$ (3)

      式中,δ(n)是单位脉冲序列; *表示离散卷积。这种对偶滤波器的最简单选择就是h(n)= g(n)=δ(n),代入式(3)可得到:

      $$ g\left( n \right) = \delta \left( n \right) - h\left( n \right) $$ (4)

      由式(2)~(4)可得到细节系数计算公式:

      $$ {d_j}\left( n \right) = {a_{j - 1}}\left( n \right) - {a_j}\left( n \right),j = 1,2 \cdots J $$ (5)

      原始信号的重建可简单地通过式(6)得到:

      $$ a\left( n \right) = {a_J}\left( n \right) + \sum\limits_{j = 1}^J {{d_j}\left( n \right)} $$ (6)

      根据文献[19-22],在计算中选用了高阶Daubechies小波构造的低通滤波器, 由于没有对信号进行二次抽样,因此可以通过对低通滤波器进行增采样来代替。滤波器的增采样是通过在第j层每两个滤波器系数中间插入2j-1-1个零值实现的。

      多孔算法的编程实现步骤如下:

      1) 将j初始化为零并存储输入信号aj(n);

      2) j加1, 利用低通滤波器haj-1(n)进行离散卷积;

      3) 获取第j层的细节系数,dj(n)=aj-1(n)-aj(n);

      4) 如果j小于分解层数J,对低通滤波器h进行增采样,继续第2)步;

      5) 利用集合D={d0, d1dJ, aJ}表示数据的小波变换。

    • 假设S0是从一个均匀采样信号中获取的时间序列。利用Haar小波对S0进行分解后,从高频率区域获取的序列dm的能量可以认为是离散信号S0的频率指标。按照这一观点,任意信号S0≠0的频率指标可以定义为:

      $$ \omega \left( {{S_0}} \right) = \sum\limits_{m = a}^b {\frac{{{{\left\| {{d_m}} \right\|}^2}}}{{{{\left\| {{S_0}} \right\|}^2}}}} $$ (7)

      式中,ω(S0)反映了离散信号S0的频率特征。对于原始信号,由于存在大量高频信息,影响了对低频信息的反映,因此,对原始信号更注重其高频部分。一般通过3层分解(dm, m=1, 2…M)就能得到反映原始信号高频信息的最优结果,对于某些特定的应用需求可以考虑减小或者增加层数。

      假设(y1, y2ynyN)为一均匀采样信号,每个信号采样用下标n标记时间点。局部的频率指标可通过定义ωn=ω(bn)计算得到,其中bnyn前后长度为R的序列,这里R取8。

      $$ {b_n} = \left( {{y_{n - \left( {R/2 - 1} \right)}},{y_{n - \left( {R/2 - 2} \right)}} \cdots {y_{n + R/2}}} \right) $$ (8)

      式中,bn包括yn左边的3个采样和右边的4个采样以及yn; bn中未定义的采样用零代替,然后利用Haar小波变换计算得到bn的频率指标ωn

      对原始采样信号时间序列和对应的频率指标的时间序列进行异常模式搜索,下一步就是通过识别算法去发现异常模式。

    • 在时间序列ynωn中,所有长度为l的子序列可以通过滑动窗口的形式表示。异常模式的识别首先就是通过算法在时间序列ynωn中分别检测最短长度为l的相似向量,然后对搜索结果进行拼接,即可得到任意长度大于l的重复异常模式。

      首先建立矩阵 YW,分别代表时间序列ynωn中所有的长度为l的时间片段。

      $$ \mathit{\boldsymbol{Y}} = \left[ {\begin{array}{*{20}{c}} {{y_1}}& \cdots &{{y_l}}\\ \vdots&\ddots&\vdots \\ {{y_{N - l + 1}}}& \cdots &{{y_N}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{Y}}_1}}\\ \vdots \\ {{\mathit{\boldsymbol{Y}}_{N - l + 1}}} \end{array}} \right] $$ (9)
      $$ \mathit{\boldsymbol{W}} = \left[ {\begin{array}{*{20}{c}} {{w_1}}& \cdots &{{w_l}}\\ \vdots&\ddots&\vdots \\ {{w_{N - l + 1}}}& \cdots &{{w_N}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{W}}_1}}\\ \vdots \\ {{\mathit{\boldsymbol{W}}_{N - l + 1}}} \end{array}} \right] $$ (10)

      然后,分别对YW做相关性分析,得到相关性系数矩阵RYRW,以幅值搜索为例,频率指标的搜索方法相同, 有:

      $$ {\mathit{\boldsymbol{R}}_\mathit{\boldsymbol{Y}}} = \left[ {\begin{array}{*{20}{c}} {r\left( {{\mathit{\boldsymbol{Y}}_1},{\mathit{\boldsymbol{Y}}_1}} \right)}& \cdots &{r\left( {{\mathit{\boldsymbol{Y}}_1},{\mathit{\boldsymbol{Y}}_{N - l + 1}}} \right)}\\ \vdots &{}& \vdots \\ {r\left( {{\mathit{\boldsymbol{Y}}_i},{\mathit{\boldsymbol{Y}}_1}} \right)}& \cdots &{r\left( {{\mathit{\boldsymbol{Y}}_i},{\mathit{\boldsymbol{Y}}_{N - l + 1}}} \right)}\\ \vdots &{}& \vdots \\ {r\left( {{\mathit{\boldsymbol{Y}}_{N - l + 1}},{\mathit{\boldsymbol{Y}}_1}} \right)}& \cdots &{r\left( {{\mathit{\boldsymbol{Y}}_{N - l + 1}},{\mathit{\boldsymbol{Y}}_{N - l + 1}}} \right)} \end{array}} \right] $$ (11)

      由于只判断没有重叠部分的YiYj,仅当|i-j| >l时为有效的相似片段,且RY为对称矩阵,因此截取RY的1至N-2l+1行和l+1至N-l+1列为R′Y,这样可以使得R′Y为包含所有不重叠片段相关性信息的最小矩阵,有利于降低内存消耗和减少计算时间, 可表示为:

      $$ {{\mathit{\boldsymbol{R'}}}_\mathit{\boldsymbol{Y}}} = \left[ {\begin{array}{*{20}{c}} {r\left( {{\mathit{\boldsymbol{Y}}_1},{\mathit{\boldsymbol{Y}}_{l + 1}}} \right)}& \cdots &{r\left( {{\mathit{\boldsymbol{Y}}_1},{\mathit{\boldsymbol{Y}}_{N - l + 1}}} \right)}\\ \vdots &{}& \vdots \\ {r\left( {{\mathit{\boldsymbol{Y}}_{N - 2l + 1}},{\mathit{\boldsymbol{Y}}_{l + 1}}} \right)}& \cdots &{r\left( {{\mathit{\boldsymbol{Y}}_{N - 2l + 1}},{\mathit{\boldsymbol{Y}}_{N - l + 1}}} \right)} \end{array}} \right] $$ (12)

      r1r2分别为幅值和频率指标最小相关性,若r(Yi, Yj)>r1则认为观测值的两个片段幅值相似,同样,若r(Wi, Wj) >r2, 则认为两个片段频率相似。对于一组起点为ij的不同片段,若同时满足幅值相似和频率相似,则认为该组片段为相似片段。最后,对相似片段进行拼接,即可得到任意长度大于l的相似片段。算法流程图如图 1所示。

      图  1  算法流程图

      Figure 1.  Flow Chart of Algorithm

    • 多孔小波分解结果的不同分解层数对应的频率范围不同,且没有相位和采样率的变化,因此对不同频率范围的信号的搜索应在与之对应的分解层内进行。此外,不同的频率范围应对应着不同的最小长度l,在各层进行搜索时,应选择不同的l

      一个长度为T,采样间隔为dt的时间序列,其包含的频率范围为$\left[{\frac{1}{T}, \frac{1}{{2{\rm{d}}t}}} \right]$,对应的周期范围为[2dt, T],小波分解每次将高频部分的一半分离,如第1层结果对应的频率范围为$\left[{\frac{1}{2}, \left({\frac{1}{{2{\rm{d}}t}}-\frac{1}{T}} \right), \frac{1}{{2{\rm{d}}t}}} \right]$,第2层的结果对应的频率范围为$\left[{\frac{1}{{{\rm{8d}}t}}-\frac{3}{{4T}}, \frac{1}{{{\rm{4d}}t}}-\frac{1}{{2T}}} \right]$,对于高频部分,由于其周期短,且出现次数多,选择l为感兴趣的最小长度;对于低频部分,由于其周期长,出现次数少,最小长度l为其最小周期。根据经验选择l依据如下:对于周期小于20 min信号的搜索,定义l=20;对于周期大于20 min信号的搜索,定义l为该层对应的最小周期。

    • 为了验证本文算法的有效性,利用模拟数据进行了仿真实验。实验包括两个部分:第1部分模拟不同信噪比(signal noise ratio, SNR)信号的搜索结果,第2部分模拟在多孔小波分解结果中进行搜索。

    • 模拟一段长度为1 500 s的信号,在其中加入4个周期分别为16、32、64、128 s的正弦信号,另外还加入了周期分别为16、32、64 s的衰减正弦信号(见图 2),相同周期信号用相同颜色进行标注。然后加入一个随机生成的功率为1的高斯噪声。由于振幅为1、均值为0的正弦信号的平均功率为固定值0.707 1,通过调节其振幅实现对信噪比的调节,分别合成信噪比为1、2、4和8的合成信号进行实验,取r1=0.8,r2=0.9,得到结果如图 3表 1所示。

      图  2  模拟信号

      Figure 2.  Simulated Testing Data

      图  3  模拟信号搜寻结果

      Figure 3.  Search Results of Simulated Data

      表 1  不同模拟信号搜寻结果/%

      Table 1.  Search Results of Simulated Testing Data with Different SNR/%

      信噪比 周期/s 总比例
      16 32 64 128
      正弦 衰减 正弦 衰减 正弦 衰减 正弦
      1 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
      2 0.00 0.00 0.00 0.00 81.25 29.69 0.00 20.17
      4 0.00 32.81 100.00 60.94 96.09 50.78 43.55 58.66
      8 87.50 56.25 100.00 91.41 100.00 63.28 93.75 87.71

      图 3表 1可以看出,当信噪比为1时,模拟信号淹没在噪声中,在合成信号中搜索不到重复信号;随着信噪比增大,搜索结果逐渐增多,当信噪比为8时,模拟信号中87.71%被搜索出。

    • 选取图 3中信噪比为2的信号利用多孔小波分解方法对信号进行分解,每次分解为高频与低频两部分, 然后对低频部分进一步分解, 共分解5层, 结果如图 4所示, 并在不同分解结果中进行搜索,结果如图 5表 2所示。

      图  4  合成信号多孔小波分解结果

      Figure 4.  Wavelet Decomposition Results of Synthetic Signals

      图  5  合成信号多孔小波分解第2层和第3层搜索结果

      Figure 5.  The Second Layer and the Third Layer Search Results of Synthetic Signals

      表 2  模拟信号分层搜寻结果/%

      Table 2.  Search Results at Different Decomposition Result of Simulated Testing Data/%

      周期/s 总比例
      16 32 64 128
      正弦 衰减 正弦 衰减 正弦 衰减 正弦
      原始信号 0.00 0.00 0.00 0.00 81.25 29.69 0.00 20.17
      第1层 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
      第2层 60.94 90.63 100.00 100.00 99.22 72.66 58.59 77.63
      第3层 98.44 40.63 100.00 100.00 100.00 100.00 98.83 96.80

      图 5表 2搜索结果可以看出,原始信号由于噪声的影响,仅搜索到20.17%的模拟信号。多孔算法分解后的第1层主要为高频信号,噪声与信号频率十分接近,因此难以搜索到相似信号。而在第2层的搜索结果中,可以看到搜索的准确性明显提高,且对周期性信号识别率较高。

    • 为验证算法对实际观测数据的应用有效性,截取拉萨SG057以及武汉SG065和SG0073台超导重力仪2015-04-25尼泊尔地震前一天的秒采样数据,在扣除潮汐、气压、极移的影响后,采用3倍中误差法除去观测中的突跳,最后将3组数据依次拼接成一组数据(见图 6),并按照模拟数据的实验方法,选取r1=0.5,r2=0.9进行搜索, 结果如表 3

      图  6  拼接数据

      Figure 6.  Splicing Data

      表 3  拼接数据相似对

      Table 3.  Search Results of Splicing Data

      搜索结果 搜索结果
      1 262~297 3 096~3 131
      2 387~423 2 391~2 427
      3 749~780 2 879~2 910
      4 812~844 4 278~4 310
      5 950~982 3 385~3 417
      6 950~981 3 386~3 417
      7 951~983 3 385~3 417
      8 1 006~1 044 4 191~4 229
      9 1 137~1 168 3 001~3 032
      10 1 408~1 444 3 381~3 417
      11 1 408~1 443 3 382~3 417
      12 1 409~1 445 3 381~3 417
      13 1 409~1 442 3 384~3 417
      14 1 412~1 446 3 383~3 417
      15 1 464~1 498 2 904~2 938
      16 1 567~1 603 4 181~4 217
      17 1 851~1 885 3 291~3 325
      18 2 387~2 419 3 827~3 859

      表 3第2组结果和第18组结果可以得到,387~415段、2 391~2 419段和3 831~3 859段为相似段(见图 7),且分别属于3台仪器,由图 7(b)看出,3段数据的频率指标具有明显的相似性,而其幅值却很难找到该相似性,这也体现了综合频率指标和幅值作为判断依据的优点。从时间上看,2 391~2 419段和3 831~3 859段对应的是武汉两台仪器的相同时间,而比拉萨台晚了564 min,这可能与观测点的空间分布有关,拉萨超导重力仪距尼泊尔地震震中较近,而武汉两台超导重力仪距震中较远。

      图  7  拼接数据搜索结果

      Figure 7.  Part of Search Results of Splicing Data

      为了在不同分辨率下发现信号中的异常模式,利用多孔算法将拼接数据分解为5层,并针对每一层分解进行了异常模式搜索,结果见表 4。在第3层的搜索结果中,595~648段、1 847~1 900段和3 593~3 646段相似,同样属于3台仪器。虽然目前还不清楚造成这些异常重复模式的原因,但能够说明该方法可以作为一种有效的异常模式搜索手段,可从不同尺度识别连续重力观测数据中由仪器故障、环境变化和未知原因造成的异常。

      表 4  拼接数据分层搜索结果

      Table 4.  Search Results at Different Decomposition Result of Splicing Data

      搜索结果 搜索结果
      第1层 640~6741
      442~1 474
      2 608~2 6423
      528~3 560
      368~416 2 664~ 2712
      406~464 3 306~3 364
      第2层 526~575 3 426~3 475
      940~991 3 156~3 207
      1 131~1 178 2 990~3 037
      1 312~1 359 3 260~3 307
      26~75 3 040~3 089
      47~94 1 960~2 007
      232~282 3 293~3 343
      233~280 3 293~3 340
      514~568 3 085~3 139
      第3层 595~648 3 593~3 646
      672~725 2 799~2 852
      762~811 3 432~3 481
      1 023~1 081 2 974~3 032
      1 043~1 093 4 229~4 279
      1 847~1 913 3 593~3 659
      1 877~1 928 3 317~3 368
      17~83 1 963~2 029
      19~83 1 964~2 028
      19~97 2 899~2 977
      21~98 2 902~2 979
      65~131 1 957~2 023
      第4层 1 282~1 349 3 406~3 473
      1 400~1 466 3 136~3 202
      1 957~2 025 4 045~4 113
      1 958~2 034 4 045~4 121
      2 390~2 461 4 122~4 193
      2 616~2 679 4 116~4 179
      530~660 3 735~3 865
      1 340~1 504 3 033~3 197
      第5层 1 340~1 502 3 034~3 196
      1 340~1 485 3 035~3 180
      1 545~1 695 2 990~3 140
      1 545~1 695 2 991~3 141
    • 本文在小波多孔算法的基础上,提出了一种综合信号频率和幅值信息的连续重力观测数据多分辨率异常模式识别算法,并利用模拟数据和实际观测数据对算法开展了仿真与实际识别实验。仿真实验表明,通过不同小波分解结果的搜索,可以发现原始信号中较难识别的长周期信号;且利用频率指标和幅值结合的方法比仅通过幅值进行搜索有明显优势。在实际重力数据实验中,原始数据中搜索到3台仪器均存在的一段27 min的异常模式,其在频率指标上相关性在90%以上,幅值的相关性在50%以上。仿真实验和实际重力数据实验结果表明,该算法不仅能够准确地在带有噪声的模拟数据中识别人为添加的异常模式,同时也能够在实际观测数据中发现一些平时难以发现的微小异常。这种震前异常的普遍性仍有待进一步验证,实际观测数据搜索结果的更深层次物理解释需要进一步研究。

参考文献 (22)

目录

    /

    返回文章
    返回