留言板

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

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

基于改进累积和控制图的GNSS变形信息的识别与预警

吴昊 刘超 赵兴旺

吴昊, 刘超, 赵兴旺. 基于改进累积和控制图的GNSS变形信息的识别与预警[J]. 武汉大学学报 ● 信息科学版, 2020, 45(10): 1517-1525. doi: 10.13203/j.whugis20190003
引用本文: 吴昊, 刘超, 赵兴旺. 基于改进累积和控制图的GNSS变形信息的识别与预警[J]. 武汉大学学报 ● 信息科学版, 2020, 45(10): 1517-1525. doi: 10.13203/j.whugis20190003
WU Hao, LIU Chao, ZHAO Xingwang. Identification and Forewarning of GNSS Deformation Information Based on a Modified Cumulative Sum Control Chart[J]. Geomatics and Information Science of Wuhan University, 2020, 45(10): 1517-1525. doi: 10.13203/j.whugis20190003
Citation: WU Hao, LIU Chao, ZHAO Xingwang. Identification and Forewarning of GNSS Deformation Information Based on a Modified Cumulative Sum Control Chart[J]. Geomatics and Information Science of Wuhan University, 2020, 45(10): 1517-1525. doi: 10.13203/j.whugis20190003

基于改进累积和控制图的GNSS变形信息的识别与预警

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

国家自然科学基金 41474026

详细信息

Identification and Forewarning of GNSS Deformation Information Based on a Modified Cumulative Sum Control Chart

Funds: 

The National Natural Science Foundation of China 41474026

More Information
    Author Bio:

    WU Hao, postgraduate, specializes in early warning and forecast of deformation disaster. E-mail: haowu0505@gmail.com

    Corresponding author: LIU Chao, PhD, associate professor. E-mail: chaoliu0202@gmail.com
  • 摘要: 累积和(cumulative sum,CUSUM)控制图是一种较为成熟的突变点识别方法,可以较为有效地对全球导航卫星系统(global navigation satellite system,GNSS)坐标序列中的变形信息进行识别与预警。然而,对于长时间的连续监测,GNSS坐标序列中不可避免地存在粗差或离群值,其对CUSUM控制图的检测效果有一定的影响。鉴于此,提出了采用基于中位数构建的改进CUSUM控制图进行GNSS坐标序列中变形信息的识别与预警,并给出了该方法的基本原理和算法流程。结果表明:采用改进CUSUM控制图进行变形信息识别和预警的准确度优于传统方法,其灾害预警的误报率显著减小,说明了采用改进CUSUM控制图的必要性和优越性。
  • 图  1  偏移值Δ与最优决策量(h, k, ARL)的关系图

    Figure  1.  Relationship Between Offset Value(Δ) and Optimal Decision Quantity(h, k, ARL)

    图  2  实验采用的GNSS坐标序列

    Figure  2.  GNSS Coordinate Series Used in the Experiment

    图  3  含有不同偏移量的GNSS坐标序列

    Figure  3.  GNSS Coordinate Series with Different Offsets

    图  4  经典CUSUM控制图对不同标准差变形数据的识别结果

    Figure  4.  Recognition Results of CUSUM Control Chart for Different Standard Deviation Deformation Data

    图  5  改进CUSUM控制图对不同标准差变形数据的识别结果

    Figure  5.  Recognition Results of Modified CUSUM Control Chart for Different Standard Deviation Deformation Data

    图  6  含有突变型变形信息的GNSS坐标序列

    Figure  6.  GNSS Coordinate Series with an Abrupt Change

    图  7  经典CUSUM控制图对连续突变型变形数据的识别结果

    Figure  7.  Recognition Results of the Abrupt Changes in GNSS Coordinate Series by CUSUM

    图  8  改进CUSUM控制图对连续突变型变形数据的识别结果

    Figure  8.  Recognition Results of the Abrupt Changes in GNSS Coordinate Series by the Modified CUSUM

    表  1  GNSS坐标序列的正态性检验结果

    Table  1.   Normality Test Results of GNSS Coordinate Series

    参数 结果 参数 结果
    平均值 0.004 74 峰度 0.141 71
    标准误差 0.000 04 偏度 0.094 71
    中位数 0.004 61 W 0.997 21
    标准差 0.001 66 P 0.004 91
    下载: 导出CSV

    表  2  CUSUM及其改进方法对变形信息区间的识别结果/历元

    Table  2.   Recognition Results of CUSUM and Its Modified Method for Deformation Information/Epoch

    变形信息 识别结果
    经典CUSUM 改进CUSUM
    1倍标准差偏移 245~373 258~303
    2倍标准差偏移 244~443 252~345
    3倍标准差偏移 244~464 246~360
    区间真值 241~300
    下载: 导出CSV

    表  3  CUSUM及其改进方法对含有离散粗差的变形序列的识别结果/历元

    Table  3.   Recognition Results of CUSUM and Its Modified Method for Deformation Series With Discrete Gross Errors/Epoch

    变形信息 检验方向 识别结果
    经典CUSUM 改进CUSUM
    1倍标准差偏移 上偏移 16~20, 151~152, 245~373 258~303
    下偏移 45~47, 84, 87, 610~712, 1 243~1 245, 1 282, 1 285
    2倍标准差偏移 上偏移 244~446 252~345
    下偏移 45~47, 67~108, 361~370, 644~687, 1 243~1 245, 1 277~1 286
    3倍标准差偏移 上偏移 42~54, 244~491, 1 212~1 218 246~357
    下偏移 40, 41, 87, 639~650, 654~657, 660, 663, 664, 670~704, 1 238~1 240, 1 243~1 249, 1 253~1 256, 1 259, 1 262~1 263, 1 270~1 303
    下载: 导出CSV

    表  4  基于改进CUSUM控制图对含有连续粗差的变形序列的识别结果/历元

    Table  4.   Recognition Results of Deformation Series with Continuous Gross Error Based on the Modified CUSUM/Epoch

    样本容量 方案 改进CUSUM
    上偏移识别结果 下偏移识别结果
    3 1 252~342 -
    2 252~342 -
    3 252~345, 729~735 -
    4 252~342, 1 272~1 281 -
    5 1 250~385, 590 80, 85, 90, 390, 680~685, 1 240~1 245, 1 255
    2 250~385, 590 80, 85, 90, 390, 625, 680~685, 1 240~1 245, 1 255
    3 250~385, 590 80, 85, 90, 390, 625, 680~685, 1 240~1 245, 1 255
    4 250~385, 590, 1 274~1 295 80, 85, 90, 390, 625, 680~685, 1 240~1 245, 1 255
    7 1 252~378, 1 190 1 106, 1 225, 1 281~1 288
    2 252~378, 1 190 1 106, 1 225, 1 281~1 288
    3 252~378, 1 190 1 106, 1 225, 1 281~1 288
    4 252~378, 1 190 110, 612, 251, 281
    下载: 导出CSV

    表  5  CUSUM及其改进方法对连续突变型变形序列的识别结果/历元

    Table  5.   Recognition Results of the Abrupt Changes in GNSS Coordinate Series by CUSUM and Its Improved Method/Epoch

    变形信息 上偏移 下偏移
    经典CUSUM 1 209~1 500 34~149、387、633~635、638~748、986
    改进CUSUM 1 236~1 500
    区间真值 1 201~1 500
    下载: 导出CSV
  • [1] 陈永奇, James L.单历元GPS变形监测数据处理方法的研究[J].武汉测绘科技大学学报, 1998, 23(4): 324-328, 363 http://www.cnki.com.cn/Article/CJFDTotal-WHCH804.009.htm

    Chen Yongqi, James L. Development of the Methodo-logy for Single Epoch GPS Deformation Monitoring[J]. Journal of Wuhan University of Surveying and Mapping, 1998, 23(4): 324-328, 363 http://www.cnki.com.cn/Article/CJFDTotal-WHCH804.009.htm
    [2] Berardino P, Fornaro G, Lanari R, et al. A New Algorithm for Surface Deformation Monitoring Based on Small Baseline Differential SAR Interferograms[J]. IEEE Transactions on Geoscience and Remote Sensing, 2002, 40(11): 2 375-2 383 doi:  10.1109/TGRS.2002.803792
    [3] Monserrat O, Crosetto M. Deformation Measurement Using Terrestrial Laser Scanning Data and Least Squares 3D Surface Matching[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2008, 63(1): 142-154 doi:  10.1016/j.isprsjprs.2007.07.008
    [4] Im S B, Hurlebaus S, Kang Y J. Summary Review of GPS Technology for Structural Health Monitoring[J]. Journal of Structural Engineering, 2011, 139(10): 1 653-1 664 http://psr.kangnam.ac.kr/psr_bk_iss/vol15/V15_A_08.pdf
    [5] Leach M P, Hyzak M. GPS Structural Monitoring as Applied to Cable-Stayed Suspension Bridge[C]. Proceedings of the 20th FIG Congress, Melbourne, Australia, 1994
    [6] Lovse J W, Teskey W F, Lachapelle G, et al. Dynamic Deformation Monitoring of Tall Structure Using GPS Technology[J]. Journal of Surveying Engineering, 1995, 121(1): 35-40 doi:  10.1061/(ASCE)0733-9453(1995)121:1(35)
    [7] Guo J J. Research of Displacement and Frequency of Tall Building Under Wind Load Using GPS[C]. Proceedings of the 10th International Technical Meeting of Satellite Divistion of the Institute of Navigation, Missouri, USA, 1997
    [8] 姜卫平, 刘鸿飞, 刘万科, 等.西龙池上水库GPS变形监测系统研究及实现[J].武汉大学学报·信息科学版, 2012, 37(8): 949-952 http://ch.whu.edu.cn/article/id/283

    Jiang Weiping, Liu Hongfei, Liu Wanke, et al. CORS Development for Xilongchi Dam Deformation Monitoring[J]. Geomatics and Information Science of Wuhan University, 2012, 37(8): 949-952 http://ch.whu.edu.cn/article/id/283
    [9] 叶世榕, 赵乐文, 陈德忠, 等.基于北斗三频的实时变形监测数据处理[J].武汉大学学报·信息科学版, 2016, 41(6): 722-728 doi:  10.13203/j.whugis20140522

    Ye Shirong, Zhao Lewen, Chen Dezhong, et al. Real-Time Deformation Monitoring Data Processing Based on BDS Triple-Frequency Observations[J]. Geomatics and Information Science of Wuhan University, 2016, 41(6): 722-728 doi:  10.13203/j.whugis20140522
    [10] 胡东.利用BFAST算法分析GPS时间序列[D].成都: 西南交通大学, 2016

    Hu Dong. GPS Time Series Analysis Using BFAST Algorithm[D]. Chengdu: Southwest Jiaotong University, 2016
    [11] 伊廷华, 郭庆, 李宏男, 等. GPS异常监测数据的关联负选择分步识别算法[J].振动工程学报, 2015, 28(1): 1-8 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=zdgcxb201501001

    Yi Tinghua, Guo Qing, Li Hongnan, et al. Multi-step Identification Algorithm Based on Relational Negative Selection for GPS Abnormal Monitoring Data[J]. Journal of Vibration Engineering, 2015, 28(1): 1-8 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=zdgcxb201501001
    [12] 周纪芗, 茆诗松.质量管理统计方法[M].北京:中国统计出版社, 2008

    Zhou Jixiang, Mao Shisong. Statistical Methods for Quality Management[M]. Beijing: China Statistics Press, 2008
    [13] Miao Changqing, Wang Man, Tian Hongjin, et al. Damage Alarming of Long-span Suspension Bridge Based on GPS-RTK Monitoring[J]. Journal of Central South University, 2015, 22 (7): 2 800-2 808 doi:  10.1007/s11771-015-2811-4
    [14] Ogaja C, Wang J, Rizos C. Detection of Wind-Induced Response by Wavelet Transformed GPS Solutions[J]. Journal of Surveying Engineering, 2003, 129(3): 99-104 doi:  10.1061/(ASCE)0733-9453(2003)129:3(99)
    [15] Iz H B. Differencing Reveals Hidden Changes in Baseline Length Time-Series[J]. Journal of Geodesy, 2006, 80(5): 259-269 doi:  10.1007/s00190-006-0066-4
    [16] 郭庆. GPS结构监测数据失真分析与异常识别方法研究[D].大连: 大连理工大学, 2013

    Guo Qing. Research on Distortion Analysis and Abnormally Recognition of GPS Structural Monitoring Data[D]. Dalian: Dalian University of Technology, 2013
    [17] 伊廷华, 郭庆, 李宏男.基于控制图GPS异常监测数据检验方法研究[J].工程力学, 2013, 30(8): 133-141 http://www.cqvip.com/QK/95324X/201308/46715812.html

    Yi Tinghua, Guo Qing, Li Hongnan. The Research on Detection Methods of GPS Abnormal Monitoring Data Based on Control Chart[J]. Engineering Mechanics, 2013, 30(8): 133-141 http://www.cqvip.com/QK/95324X/201308/46715812.html
    [18] Quesenberry C P. On Properties of Binomial Q Charts for Variables[J]. Journal of Quality Technology, 1995, 27(3): 184-203 doi:  10.1080/00224065.1995.11979592
    [19] Mertikas S P, Damianidis K I. Monitoring the Quality of GPS Station Coordinates in Real Time[J]. GPS Solutions, 2007, 11(2): 119-128 doi:  10.1007/s10291-006-0044-6
    [20] Tran K P, Castagliola P, Celano G. Monitoring the Ratio of Population Means of a Bivariate Normal Distribution Using CUSUM Type Control Charts[J]. Statistical Papers, 2018, 59(1): 387-413 doi:  10.1007/s00362-016-0769-4
    [21] 胡雪龙, 周晓剑, 黄卫东, 等. CUSUM中位数控制图设计[J].数理统计与管理, 2018, 37(3): 469-477 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=sltjygl201803009

    Hu Xuelong, Zhou Xiaojian, Huang Weidong, et al. Design of the CUSUM Median Chart[J]. Journal of Applied Statistics and Management, 2018, 37(3): 469-477 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=sltjygl201803009
    [22] Shaphiro S S, Wilk M B. An Analysis of Variance Test for Normality[J]. Biometrika, 1965, 52(3): 591-611
    [23] Silverman B W. Using Kernel Density Estimates to Investigate Multimodality[J]. Journal of the Royal Statistical Society Series B (Methodological), 1981, 43(1): 97-99 doi:  10.1111/j.2517-6161.1981.tb01155.x
    [24] Parzen E. On Estimation of a Probability Density Function and Mode[J]. The Annals of Mathematical Statistics, 1962, 33(3): 1 065-1 076 doi:  10.1214/aoms/1177704472
    [25] Siegmund D. Error Probabilities and Average Sample Number of the Sequential Probability Ratio Test[J]. Journal of the Royal Statistical Society Series B (Methodological), 1975, 37(3): 394-401 doi:  10.1111/j.2517-6161.1975.tb01552.x
  • [1] 王晨辉, 郭伟, 孟庆佳, 刘炎炎, 毕逢东.  基于虚拟参考站的GNSS滑坡变形监测方法及性能分析 . 武汉大学学报 ● 信息科学版, 2022, 47(6): 990-996. doi: 10.13203/j.whugis20220102
    [2] 姜卫平, 梁娱涵, 余再康, 肖玉钢, 陈剡, 陈渠森.  卫星定位技术在水利工程变形监测中的应用进展与思考 . 武汉大学学报 ● 信息科学版, 2022, 47(10): 1625-1634. doi: 10.13203/j.whugis20220589
    [3] 党亚民, 王伟, 章传银, 白贵霞, 杨强.  综合GNSS和重力数据定量评价三峡地区地质环境稳定性 . 武汉大学学报 ● 信息科学版, 2020, 45(7): 1052-1057. doi: 10.13203/j.whugis20170417
    [4] 刘硕, 张磊, 李健.  GNSS系统间钟差辅助定位研究 . 武汉大学学报 ● 信息科学版, 2017, 42(7): 1015-1020. doi: 10.13203/j.whugis20150199
    [5] 张清华, 隋立芬, 贾小林, 田亮, 朱永兴.  利用高精度PPS测量进行GPS-GLONASS时差监测 . 武汉大学学报 ● 信息科学版, 2014, 39(11): 1347-1351.
    [6] 杨元喜, 曾安敏, 景一帆.  函数模型和随机模型双约束的GNSS数据融合及其性质 . 武汉大学学报 ● 信息科学版, 2014, 39(2): 127-131. doi: 10.13203/j.whugis20130378
    [7] 魏二虎, 万丽华, 金双根, 刘经南.  联合GNSS和SLR观测对地球自转参数的解算与分析 . 武汉大学学报 ● 信息科学版, 2014, 39(5): 581-585. doi: 10.13203/j.whugis20120213
    [8] 孙克文, 徐华敏, 刘伟, 丁志中.  卫星导航信号差分联合信道组合捕获策略 . 武汉大学学报 ● 信息科学版, 2013, 38(12): 1413-1419.
    [9] 章红平, 韩文慧, 黄玲, 耿长江.  地基GNSS全球电离层延迟建模 . 武汉大学学报 ● 信息科学版, 2012, 37(10): 1186-1189.
    [10] 李英冰, 闫景仙, 熊程波, 陈中新.  GNSS观测值的压缩方法研究 . 武汉大学学报 ● 信息科学版, 2012, 37(7): 818-822.
    [11] 闫伟, 欧吉坤, 袁运斌, DagobertoSalazar.  利用广播星历和区域参考网实现实时精密单点定位算法研究 . 武汉大学学报 ● 信息科学版, 2012, 37(10): 1190-1193.
    [12] 李振, 朱锋, 陈家君.  利用新息构造AKF估计和预报大坝变形 . 武汉大学学报 ● 信息科学版, 2012, 37(4): 454-457.
    [13] 张绍成, 郭际明, 孟祥广.  非差虚拟观测值的网络RTK实现 . 武汉大学学报 ● 信息科学版, 2011, 36(10): 1226-1230.
    [14] 蔡华, 赵齐乐, 孙汉荣, 胡志刚.  GNSS实时数据质量控制 . 武汉大学学报 ● 信息科学版, 2011, 36(7): 820-824.
    [15] 张飞舟, 杨东凯, 陈嘉, 杨伯钢.  GNSS软件接收机的实现与测试分析 . 武汉大学学报 ● 信息科学版, 2009, 34(5): 577-580.
    [16] 韩保民, 欧吉坤, 柴艳菊.  用拟准检定法探测和修复GPS数据中的粗差和周跳 . 武汉大学学报 ● 信息科学版, 2002, 27(3): 246-250.
    [17] 郑肇葆, 谭春健, 张润根.  按残差绝对值和最小原理的光束法区域网平差 . 武汉大学学报 ● 信息科学版, 1990, 15(2): 34-40.
    [18] 金新祥.  观测值残差和未知参数一次范数最小平差 . 武汉大学学报 ● 信息科学版, 1989, 14(3): 59-67.
    [19] 张方仁, 王勇.  主成份检验和粗差定位 . 武汉大学学报 ● 信息科学版, 1986, 11(4): 44-56.
    [20] 周东旭, 冯义楷, 张化疑, 付延光, 唐秋华.  联合卫星测高和GNSS观测的天津沿海近25年相对海平面变化分析 . 武汉大学学报 ● 信息科学版, 0, 0(0): 0-0. doi: 10.13203/j.whugis20210532
  • 加载中
图(8) / 表(5)
计量
  • 文章访问数:  641
  • HTML全文浏览量:  137
  • PDF下载量:  64
  • 被引次数: 0
出版历程
  • 收稿日期:  2019-12-27
  • 刊出日期:  2020-10-05

基于改进累积和控制图的GNSS变形信息的识别与预警

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

    国家自然科学基金 41474026

    作者简介:

    吴昊,硕士生,主要从事变形灾害的预警与预报。haowu0505@gmail.com

    通讯作者: 刘超,博士,副教授。chaoliu0202@gmail.com
  • 中图分类号: P228

摘要: 累积和(cumulative sum,CUSUM)控制图是一种较为成熟的突变点识别方法,可以较为有效地对全球导航卫星系统(global navigation satellite system,GNSS)坐标序列中的变形信息进行识别与预警。然而,对于长时间的连续监测,GNSS坐标序列中不可避免地存在粗差或离群值,其对CUSUM控制图的检测效果有一定的影响。鉴于此,提出了采用基于中位数构建的改进CUSUM控制图进行GNSS坐标序列中变形信息的识别与预警,并给出了该方法的基本原理和算法流程。结果表明:采用改进CUSUM控制图进行变形信息识别和预警的准确度优于传统方法,其灾害预警的误报率显著减小,说明了采用改进CUSUM控制图的必要性和优越性。

English Abstract

吴昊, 刘超, 赵兴旺. 基于改进累积和控制图的GNSS变形信息的识别与预警[J]. 武汉大学学报 ● 信息科学版, 2020, 45(10): 1517-1525. doi: 10.13203/j.whugis20190003
引用本文: 吴昊, 刘超, 赵兴旺. 基于改进累积和控制图的GNSS变形信息的识别与预警[J]. 武汉大学学报 ● 信息科学版, 2020, 45(10): 1517-1525. doi: 10.13203/j.whugis20190003
WU Hao, LIU Chao, ZHAO Xingwang. Identification and Forewarning of GNSS Deformation Information Based on a Modified Cumulative Sum Control Chart[J]. Geomatics and Information Science of Wuhan University, 2020, 45(10): 1517-1525. doi: 10.13203/j.whugis20190003
Citation: WU Hao, LIU Chao, ZHAO Xingwang. Identification and Forewarning of GNSS Deformation Information Based on a Modified Cumulative Sum Control Chart[J]. Geomatics and Information Science of Wuhan University, 2020, 45(10): 1517-1525. doi: 10.13203/j.whugis20190003
  • 高层建筑与大型桥梁变形、城市地表沉陷和高速公路边坡滑坡等风险,对人民生命财产安全构成了巨大的威胁。因此,有必要对变形体进行有效监测,获取变形信息,采用一定的数学方法进行数据处理与分析,进而进行有效的预警预报,减小灾害发生概率和影响范围[1-3]

    全球导航卫星系统(global navigation satellite system,GNSS)作为一种高精度定位技术,可以为变形体提供连续实时动态监测,已成为一种变形监测的重要技术手段,在相关变形监测领域得到了广泛应用[4-9]。采用GNSS进行动态变形监测,其主要目的之一是对变形灾害进行实时识别和预警。对于如何准确识别GNSS坐标序列中的变形信息,相关学者展开了广泛的研究。胡东[10]提出了利用添加季节和趋势分解与断点识别(breaks for additive season and trend,BFAST)算法分析GPS坐标序列,并利用蒙特卡洛方法对数据进行多尺度模拟,通过实验验证了此方法的有效性。伊廷华等[11]为了有效地判别GPS坐标序列的突变点,建立了GPS监测序列的关联模型,提出了基于关联负选择的坐标序列突变点的分步识别算法。在诸多研究突变点识别理论方法中,控制图理论以其较为成熟的理论和良好的检测效果,在坐标序列的突变点识别中得到了关注。常见的控制图理论包括休哈特控制图、累积和(cumulative sum,CUSUM)控制图和移动加权控制图等方法[12]。Miao等[13]应用均值控制图对桥梁的GNSS动态坐标异常变化进行统计模式识别,得到了很好的识别效果;Ogaja等[14]采用CUSUM控制图进行结构振动频率突变的检测,验证了该方法进行较小突变检测的可行性;Iz[15]采用CUSUM控制图对GPS和甚长基线干涉测量(very long baseline interferometry,VLBI)的差异进行突变检测,获取地震活动引起的基线的瞬时变化特征;郭庆[16]基于CUSUM控制图构建了GPS异常监测数据的预警模型,实现了不同标准差条件下突变点的有效探测;伊廷华等[17]对比分析了休哈特控制图和CUSUM控制图在标准差偏移方面的识别能力,发现两者各有优势,相互补充。

    虽然GNSS动态变形监测通常采用载波相对定位的方式可以消除其卫星钟差和接收机钟差,也可以将对流层误差和电离层误差等空间相关误差大大削弱,但多路径误差和衍射误差等非建模系统误差仍然影响监测结果,因此,GNSS坐标序列通常不服从正态分布,难以满足控制图相关研究规定的样本数据服从正态分布的要求。虽然相关学者提出了解决这一问题的途径,如Quesenberry等[18]提出的Q统计量的方法,对原始GNSS坐标序列进行转换,使转换后的坐标序列近似地服从正态分布,但考虑到GNSS长时间监测过程中可能出现的个别粗差或离群值,相关的正态化方法难以对其进行有效处理,其仍会影响采用相关控制图方法进行GNSS变形监测突变点的识别精度。鉴于此,本文以CUSUM控制图方法为例,提出采用基于中位数构建的改进CUSUM控制图进行GNSS坐标序列中变形灾害的识别,并与经典的CUSUM方法进行对比分析,得出了一些有益的结论。

    • 控制图是统计质量管理的一种重要手段和工具,其通常基于采集数据的样本均值和假设检验的原理所构造,用于监测生产过程是否处于控制状态。在变形监测领域,若被监测体处于安全状况或受到较小变形的影响,则GNSS监测数据具有相对稳定的均值和方差;若被监测体出现一定程度的变形,GNSS监测数据相应地出现一定程度的偏移,但其方差未发生改变。因此,可基于控制图理论建立变形灾害的识别模型。

      具体描述如下:对于GNSS监测的坐标序列Xt),t=1,2…nn为样本大小,认为其近似服从正态分布,即X(t)~N (μ0σ),其中,μ0为坐标序列的均值,σ为坐标序列的标准差。若监测体发生变形,其监测的坐标序列近似服从X(t)~N(μ0 + δσσ),其中,δσ为坐标序列的均值偏移量[19]

      μ1 = μ0 + δσ,进行假设检验。原假设H0为坐标序列未发生异常波动,其均值不变;备选假设H1为坐标序列在ss < n)时刻发生异常波动,即s时刻发生变形,认为在s时刻以后的监测数据的均值发生改变。

      在坐标序列服从正态分布的前提下,s时刻前后的Xt)的概率密度函数可以表示为:

      $$ \left\{ \begin{array}{l} {f_0}\left( x \right) = \frac{1}{{\sigma \sqrt {2{\rm{ \mathsf{ π} }}} }}{\rm{exp}}\left[ { - \frac{{{{\left( {x - {\mu _0}} \right)}^2}}}{{2{\sigma ^2}}}} \right], {H_0}成立\\ {f_1}\left( x \right) = \frac{1}{{\sigma \sqrt {2{\rm{ \mathsf{ π} }}} }}{\rm{exp}}\left[ { - \frac{{{{\left( {x - {\mu _1}} \right)}^2}}}{{2{\sigma ^2}}}} \right], {H_1}成立 \end{array} \right. $$ (1)

      构建对数似然比统计量:

      $$ \begin{array}{l} {\lambda _n} = \mathop \sum \limits_{t = 1}^n {\rm{ln}}\frac{{{f_1}\left( {{x_t}} \right)}}{{{f_0}\left( {{x_t}} \right)}} = \\ \;\;\;\;\;\;\;\mathop \sum \limits_{t = 1}^n \frac{1}{{2{\sigma ^2}}}\left[ {{{\left( {{x_t} - {\mu _0}} \right)}^2} - {{\left( {{x_t} - {\mu _1}} \right)}^2}} \right] \end{array} $$ (2)

      令均值偏移量μ1 - μ0 = δσ = Δk = Δ/2。则

      $$ {\lambda _n} = \mathop \sum \limits_{t = 1}^n \frac{\mathit{\Delta} }{{{\sigma ^2}}}\left( {{x_t} - {\mu _0} - k} \right) $$ (3)

      Δ > 0,则式(3)等价于$C_t^ + = {\rm{max}}\mathop \sum \limits_{j = 1}^t \left( {{x_j} - {\mu _0} - k} \right) $[17],其中,Ct+Δ > 0条件下的似然比统计量,其可以进一步表示为:

      $$ \begin{array}{*{20}{l}} {C_t^ + = C_{t - 1}^ + + \left( {{x_t} - {\mu _0} - k} \right) = }\\ {{\rm{max}}\left( {0, C_{t - 1}^ + + {x_t} - {\mu _0} - k} \right)} \end{array} $$ (4)

      式(4)即为CUSUM的上偏移统计量,同理,可得CUSUM的下偏移统计量为:

      $$C_t^ - = {\rm{min}}\left( {0, C_{t - 1}^ - + {x_t} - {\mu _0} + k} \right) $$ (5)

      为了削弱样本中可能出现的粗差或离群值对变形信息识别精度的影响,采用中位数代替原始样本进行似然比统计量的构建。中位数计算方法如下:

      $$ {{\tilde x}_i} = {\rm{median}}\left( {{x_{m \cdot \left( {i - 1} \right) + 1}}, {x_{m \cdot \left( {i - 1} \right) + 2}} \cdots {x_{m \cdot i}}} \right) $$ (6)

      式中,$ {{\tilde x}_i} $为原始样本的中位数;i为形成的中位数样本数目,i=1,2…INT(n/m),INT(·)为取整函数,m为设置的计算中位数的子样本容量;median(·)为中位数函数。

      采用式(6)得到的中位数样本构建改进的CUSUM(称为CUSUM中位数)的上、下偏移统计量:

      $$ \left\{ \begin{array}{l} {{\tilde C}_i}^ + = {\rm{max}}\left( {0, \tilde C_{i - 1}^ + + {{\tilde x}_i} - {{\tilde \mu }_0} - \tilde k} \right)\\ {{\tilde C}_i}^ - = {\rm{min}}\left( {0, \tilde C_{i - 1}^ - + {{\tilde x}_i} - {{\tilde \mu }_0} + \tilde k} \right) \end{array} \right. $$ (7)

      式中,$ {{{\tilde \mu }_0}}$为中位数样本的均值;$ \tilde k = \frac{{{{\tilde \mu }_1} - {{\tilde \mu }_0}}}{2} $,$ {{{\tilde \mu }_1}} $为发生变形后的中位数样本的均值。

      基于式(7)上、下统计量构建相应的变形灾害预警模型如下:

      $$ \left\{ \begin{array}{l} {{\tilde C}_i}{^ + } = {\rm{max}}\left( {0, {{\tilde C}}_{i - 1}^ + + {{\tilde x}_i} - {{\tilde \mu }_0} - \tilde k} \right)\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{{\tilde C}_i}{^ + } > h, 预警 \end{array} \right. $$ (8)
      $$ \left\{ \begin{array}{l} {{\tilde C}_i}^ - = {\rm{min}}\left( {0, \tilde C_{i - 1}^ - + {{\tilde x}_i} - {{\tilde \mu }_0} + \tilde k} \right)\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;{{\tilde C}_i}^ - < - h, 预警 \end{array} \right. $$ (9)

      式中,$ {{\tilde C}_0}^ + $和${{\tilde C}_0}^ - $分别为统计量的上、下偏移的初始值,通常取值为0;h为判断变形灾害发生的阈值,即当上统计量大于h或下统计量小于-h时,视为发生变形并进行预警。

      为了确定阈值h,引入平均运行链长(average run length,ARL),ARL为控制图发出失控状态之前的平均样本数。ARL的计算方法有很多,如常见的马尔可夫链法和积分方程法。本文拟采用Tran等提出的方法进行相关计算[20]

      当检测过程处于控制状态时,ARL被表示为ARL0;相反,当检测过程失控时,ARL被表示为ARL1。ARL0值越大,表示对异常数据判别中发生的误报率(第一类错误)越低;ARL1值越小,表示对异常数据判别中发生的漏报率(第二类错误)越低。故而,关于控制图最优参数的选取,可以认为是在ARL0固定不变的情况下,ARL1的值越小越好,即:

      $$ \left\{ \begin{array}{l} {\rm{AR}}{{\rm{L}}_1}\left( {h, k, \mathit{\Delta} } \right) = {\rm{min}}\\ \;\;\;\;\;\;\;\;\;{\rm{AR}}{{\rm{L}}_0} \ge \tau \end{array} \right. $$ (10)

      式中,τ为允许的最大误报率,依据传统的6σ原理,通常设τ = 370.4[21]。在此基础上,对于不同的偏移值Δ,通过模式搜索的方法获取控制限h的最优解。为了便于描述,在构建$ {{{\tilde x}_i}} $中位数统计量时将样本容量设为奇数,以样本容量为3、5和7为例,在Δ ∈ [0.2,2.0]时计算最优决策量(hk,ARL1),结果见图 1图 1中,在偏移取值不同的情况下,k在样本容量为3和5时取值近似,与偏移值Δ的关系可表示为k=Δ/2。样本容量为7时,在Δ≤1.6之前与样本容量是3和5时的取值近似,在Δ > 1.6之后逐渐稳定到0.8。控制限h与ARL的取值随着Δ的增加而逐渐下降,样本容量越大,则数值越低。

      图  1  偏移值Δ与最优决策量(h, k, ARL)的关系图

      Figure 1.  Relationship Between Offset Value(Δ) and Optimal Decision Quantity(h, k, ARL)

    • 由于CUSUM控制图理论是基于对服从正态分布的数据进行检验,因此,在进行变形数据的识别与预警之前需要对其进行正态性判别。判定方法大致可分为图形判别和构建检验统计量两类。本文选择的方法是构建统计量中的Shapiro-Wilks检验[22]。在本文中,通过R语言采用Shapiro-Wilks检验中p-value的验证方法。

    • 由于GNSS监测数据受到多路径误差及其他非建模误差的影响,其数据通常不服从正态性,因此,在运用控制图进行变形信息的检验和预警之前,需要对其进行转换。本文采用Quesenberry等[18]提出的方法,用来解决GNSS监测数据Xt)(t=1,2…n)不服从正态分布的问题[23-24]。设n个样本点的概率密度函数为:

      $$F\left( i \right) = \frac{1}{n}\mathop \sum \limits_{t = 1}^n {K_l}\left( {{x_i} - {x_t}} \right) = \frac{1}{{nl}}K\left( {\frac{{{x_i} - {x_t}}}{l}} \right) $$ (11)

      式中,Κ (·)为核函数;Κl (·)为缩放核函数,且Kl (x) = 1/(l·K (x)),其中l为一个平滑参数;i = 1,2…nt = 1,2…n。在此基础上,对原始GNSS监测数据进行Q统计量转换:

      $$ \tilde Q\left( i \right) = {\mathit{\Phi} ^{ - 1}}\left( {F\left( i \right)} \right) $$ (12)

      式中,Φ-1 (·)是标准正态分布的累积分布函数的逆函数。通过转换后的${\tilde Q} $统计量近似服从正态分布,已基本满足对GNSS变形信息的识别与预警算法的条件。

    • 综上所述,本文提出的基于CUSUM中位数控制图的GNSS变形信息识别与预警算法描述如下:

      1)获取监测体的GNSS变形监测数据,并计算其均值和标准差;

      2)对所采集的GNSS监测数据进行Shapiro-Wilks正态性检验,将不服从正态分布的数据使用Q统计量进行转换;

      3)将数据(或转换后的数据)按照式(6)计算中位数;

      4)统计步骤3)得到的中位数,根据式(8)、式(9)构建CUSUM中位数控制图的上、下偏移统计量;

      5)依据图 1,确定异常值识别模型的参数k和控制限h

      6)对比判断CUSUM中位数上、下偏移统计量是否超过控制限h,对于超过控制限的数据定义为异常数据并进行预警。

    • 数据来源为安徽理工大学某楼顶的GNSS连续跟踪站,采用其中相距约4 m的两个站点的观测数据进行相关算法的验证。GNSS接收机均采用天宝BD980板卡,天线型号为AT300,数据源为GPS/GLONASS/BDS三系统,采样频率设置为3 Hz。截取500 s的x方向坐标序列(1 500个观测历元),并减去其真实值(由长时间观测平均得到),如图 2所示。由于两个站点的距离较近,通过载波相位双差技术可以较为有效地削弱空间相关系统误差,可以认为其坐标序列中的系统误差主要为多路径误差等非建模系统误差;同时,站点架设在建筑物的高度仅约12 m,且相对稳固,通过长时间观测,未发现明显的变形信息。因此,可以将该段坐标序列作为不含有变形信息的GNSS相对定位的观测坐标序列,并考虑加入不同特征的变形信息,进行后续算法的测试。

      图  2  实验采用的GNSS坐标序列

      Figure 2.  GNSS Coordinate Series Used in the Experiment

      进行变形信息识别和预警,需先对原始GNSS坐标序列进行正态性检验,检验结果见表 1,其中,P值小于0.05,拒绝原假设,即认为原始GNSS坐标序列不服从正态分布。采用式(11)、式(12)对原始GNSS坐标序列进行数据转换,得到近似服从正态分布的坐标序列。在原始GNSS坐标序列中分别加入无粗差影响的不同标准差偏移、离散粗差影响的偏移、连续粗差影响的偏移和连续突变型变形等4种不同形式的变形信息,验证算法对不同条件下变形信息的识别和预警的性能。

      表 1  GNSS坐标序列的正态性检验结果

      Table 1.  Normality Test Results of GNSS Coordinate Series

      参数 结果 参数 结果
      平均值 0.004 74 峰度 0.141 71
      标准误差 0.000 04 偏度 0.094 71
      中位数 0.004 61 W 0.997 21
      标准差 0.001 66 P 0.004 91
    • 在原始GNSS坐标序列的241~300历元(80~100 s)处分别加入1~3倍的标准差的上偏移,形成3组突变型变形的测试数据,如图 3所示。

      图  3  含有不同偏移量的GNSS坐标序列

      Figure 3.  GNSS Coordinate Series with Different Offsets

      分别采用经典CUSUM控制图和基于中位数的改进CUSUM控制图两种方法对上述变形信息进行识别和预警。其中,经典CUSUM控制图的控制限的确定参考Siegmund[25]所提出的ARL的计算方法;基于中位数的改进CUSUM控制图的控制限的确定参见图 1,设样本容量为3,并取Δ=1。

      两种方法的变形信息识别结果如图 4图 5所示,可知:

      图  4  经典CUSUM控制图对不同标准差变形数据的识别结果

      Figure 4.  Recognition Results of CUSUM Control Chart for Different Standard Deviation Deformation Data

      图  5  改进CUSUM控制图对不同标准差变形数据的识别结果

      Figure 5.  Recognition Results of Modified CUSUM Control Chart for Different Standard Deviation Deformation Data

      1)无论采用哪种方法构建统计量,得到的统计量序列均有效地抑制了噪声的影响,变形信息得到了有效增强。以加入1倍标准差偏移为例,由于变形量较小,在原始坐标序列中,其变形信息被噪声淹没(见图 3),而在构建的统计量序列(见图 4(a)图 5(a))中,可以较为清晰地识别加入变形信息的区间。

      2)基于中位数的改进CUSUM控制图的识别精度高于经典CUSUM控制图。表 2给出了两种方法对不同变形信息区间的识别结果,其中,对于1~3倍标准差的变形,经典CUSUM控制图的误报历元数和漏报历元数分别为73、143、164和4、3、3,而对应基于中位数的改进CUSUM控制图的误报历元数和漏报历元数分别为3、45、60和17、11、5。新方法的漏报率有所增大,但其误报率得到了较大程度的减小,识别的变形区间明显收窄,与真实区间更为接近,表明在发生变形结束后,新方法构建的统计量具有更为迅速的回归速度。

      表 2  CUSUM及其改进方法对变形信息区间的识别结果/历元

      Table 2.  Recognition Results of CUSUM and Its Modified Method for Deformation Information/Epoch

      变形信息 识别结果
      经典CUSUM 改进CUSUM
      1倍标准差偏移 245~373 258~303
      2倍标准差偏移 244~443 252~345
      3倍标准差偏移 244~464 246~360
      区间真值 241~300

      3)虽然图 3中未加入任何下偏移变形信息,但经典CUSUM控制图仍检测出了少量下偏移变形信息。对于1~3倍标准差的变形,经典CUSUM控制图分别检测出21、39、45个下偏移变形历元,而基于中位数的改进CUSUM控制图未检测出任何下偏移变形信息,说明新方法可以有效地减小预报的误警率。

    • 在加入不同标准差偏移(见图 3)的基础上,随机选择不同的位置加入一定大小的粗差,用于进一步测试算法的性能。具体方案为:

      1)在1倍标准差偏移序列的16、151、610历元处分别加入10倍、5倍、-16倍标准差大小的粗差;

      2)在2倍标准差偏移序列的67、361、418历元处分别加入-5倍、-8倍、12倍标准差大小的粗差;

      3)在3倍标准差偏移序列的42、351、1 212历元处分别加入20倍、30倍、10倍标准差大小的粗差。

      经典CUSUM控制图和基于中位数的改进CUSUM控制图对含有离散粗差的变形信息的识别结果见表 3。由表 3可知:

      表 3  CUSUM及其改进方法对含有离散粗差的变形序列的识别结果/历元

      Table 3.  Recognition Results of CUSUM and Its Modified Method for Deformation Series With Discrete Gross Errors/Epoch

      变形信息 检验方向 识别结果
      经典CUSUM 改进CUSUM
      1倍标准差偏移 上偏移 16~20, 151~152, 245~373 258~303
      下偏移 45~47, 84, 87, 610~712, 1 243~1 245, 1 282, 1 285
      2倍标准差偏移 上偏移 244~446 252~345
      下偏移 45~47, 67~108, 361~370, 644~687, 1 243~1 245, 1 277~1 286
      3倍标准差偏移 上偏移 42~54, 244~491, 1 212~1 218 246~357
      下偏移 40, 41, 87, 639~650, 654~657, 660, 663, 664, 670~704, 1 238~1 240, 1 243~1 249, 1 253~1 256, 1 259, 1 262~1 263, 1 270~1 303

      1)加入粗差后,两种方法均可以识别出加入变形的区间(241~300历元),但识别的精度有所差异。主要原因为粗差的位置与变形区间临近,影响了变形信息的识别精度;而基于中位数的改进CUSUM控制图并未受到粗差的影响,对变形区间的识别结果与未加入粗差的结果基本一致。

      2)对于加入粗差的位置,经典CUSUM控制图均错误地识别为变形信息,而基于中位数的改进CUSUM控制图未识别出任何变形信息。如对于加入1倍标准差偏差的坐标序列,在其16、151和610历元分别加入了10倍、5倍、-16倍标准差大小的粗差,经典CUSUM控制图错误识别的区间分别为16~20、151~152和610~712历元,而基于中位数的改进CUSUM控制图未作出任何变形预警;其他两个变形序列的识别结果类同。

      3)统计两种方法进行变形预警的误报历元数和漏报历元数。对于1~3倍标准差的变形,经典CUSUM控制图的误报历元数和漏报历元数分别为193、258、319和4、3、3;而对应基于中位数的改进CUSUM控制图的误报历元数和漏报历元数分别为3、45、57和17、11、5。

    • 为了进一步分析在连续粗差影响下该方法的检验性能,在加入2倍标准差大小偏移(见图 3)的基础上,随机选择不同的位置加入一定大小的连续粗差,具体方案为:

      1)在50~51历元处分别加入3倍、5倍的标准差大小的粗差;

      2)在209~210历元处分别加入3倍、-5倍的标准差大小的粗差;

      3)在726~728历元处分别加入7倍、8倍、-9倍的标准差大小的粗差;

      4)在1 269~1 271历元处分别加入7倍、8倍、9倍的标准差大小的粗差。

      由于样本容量的选择直接影响中位数的生成,进而影响基于中位数的改进CUSUM控制图方法的识别精度,选择不同的样本容量(3、5和7)进行变形信息的识别,结果如表 4所示。由表 4可知:

      表 4  基于改进CUSUM控制图对含有连续粗差的变形序列的识别结果/历元

      Table 4.  Recognition Results of Deformation Series with Continuous Gross Error Based on the Modified CUSUM/Epoch

      样本容量 方案 改进CUSUM
      上偏移识别结果 下偏移识别结果
      3 1 252~342 -
      2 252~342 -
      3 252~345, 729~735 -
      4 252~342, 1 272~1 281 -
      5 1 250~385, 590 80, 85, 90, 390, 680~685, 1 240~1 245, 1 255
      2 250~385, 590 80, 85, 90, 390, 625, 680~685, 1 240~1 245, 1 255
      3 250~385, 590 80, 85, 90, 390, 625, 680~685, 1 240~1 245, 1 255
      4 250~385, 590, 1 274~1 295 80, 85, 90, 390, 625, 680~685, 1 240~1 245, 1 255
      7 1 252~378, 1 190 1 106, 1 225, 1 281~1 288
      2 252~378, 1 190 1 106, 1 225, 1 281~1 288
      3 252~378, 1 190 1 106, 1 225, 1 281~1 288
      4 252~378, 1 190 110, 612, 251, 281

      1)与表 3相比较,对于离散的粗差点,基于中位数的改进CUSUM控制图(样本容量为3)可以很好地避免粗差点的误识别,而对于连续的粗差(见表 4),该方法(样本容量为3)将连续粗差误识别为变形信息,如方案3和4的粗差被识别的区间分别为729~735历元和1 272~1 281历元。

      2)样本容量的增大虽然可以削弱连续粗差的误报率,但对应阈值不断减小,会造成其他位置误报率的增加。当样本容量设置为7时,粗差位置并未被识别为变形信息,但是1 190历元被误识别为上偏移变形,1 106、1 225等位置被误识别为下偏移变形。因此,对于样本容量的选择,还需要根据序列的具体特征综合考虑。

    • 为了进一步测试该算法对不同类型变形信息的检验性能,设置一组连续突变型的变形信息进行测试,具体方案为:在原始监测数据(见图 2)的基础上,选择在1 201~1 500历元处加入0.01 m的偏移量作为变形数据,如图 6所示。

      图  6  含有突变型变形信息的GNSS坐标序列

      Figure 6.  GNSS Coordinate Series with an Abrupt Change

      两种方法的结果分别如图 7图 8所示,其统计结果见表 5。由图 7图 8表 5可知,与经典CUSUM控制图方法相比,基于中位数的CUSUM控制图方法的误警率明显减小。如图 7所示,虽然在实验中未添加任何下偏移变形信息,但经典CUSUM控制图检验出的下偏移统计量超过控制限的数目高达230个,如34~149、387、638~748等历元均被误识别为下偏移变形信息;而基于中位数的改进CUSUM控制图方法未检测出任何的下偏移变形信息(见图 8),误警率得到了有效控制。另外,两种方法均有效地识别出了添加的偏移量区间(1 201~1 500历元),其中基于中位数的改进CUSUM控制图的漏报率略高于经典CUSUM控制图,但综合考虑误警率和漏警率,改进CUSUM控制图的识别精度高于经典CUSUM控制图。

      图  7  经典CUSUM控制图对连续突变型变形数据的识别结果

      Figure 7.  Recognition Results of the Abrupt Changes in GNSS Coordinate Series by CUSUM

      图  8  改进CUSUM控制图对连续突变型变形数据的识别结果

      Figure 8.  Recognition Results of the Abrupt Changes in GNSS Coordinate Series by the Modified CUSUM

      表 5  CUSUM及其改进方法对连续突变型变形序列的识别结果/历元

      Table 5.  Recognition Results of the Abrupt Changes in GNSS Coordinate Series by CUSUM and Its Improved Method/Epoch

      变形信息 上偏移 下偏移
      经典CUSUM 1 209~1 500 34~149、387、633~635、638~748、986
      改进CUSUM 1 236~1 500
      区间真值 1 201~1 500
    • 考虑到GNSS变形监测中可能出现的粗差会影响变形信息的识别与预警,提出采用基于中位数的CUSUM控制图进行变形信息的识别与预警。通过理论推导和数据处理与分析,得到如下结论:

      1)CUSUM控制图以及基于中位数的改进CUSUM控制图均可以有效地抑制噪声对变形信息的影响,增强监测序列的变形特征;

      2)对于未加入粗差的不同变形序列,基于中位数的改进CUSUM控制图对变形信息的识别精度高于经典CUSUM方法;

      3)对于加入离散粗差的变形序列,基于中位数的改进CUSUM控制图优势更为明显,其表现出了良好的稳健性,与采用该方法对未加入粗差的变形序列的识别结果基本相当;对于存在连续粗差的变形序列,采用基于中位数的改进CUSUM控制图进行变形信息的识别与预警,需要设置合适的样本容量,以达到最佳的检验效果。

参考文献 (25)

目录

    /

    返回文章
    返回