留言板

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

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

改进的启发式分割算法在GNSS坐标时间序列阶跃探测中的应用

姚宜斌 冉启顺 张豹

姚宜斌, 冉启顺, 张豹. 改进的启发式分割算法在GNSS坐标时间序列阶跃探测中的应用[J]. 武汉大学学报 ( 信息科学版), 2019, 44(5): 648-654. doi: 10.13203/j.whugis20170322
引用本文: 姚宜斌, 冉启顺, 张豹. 改进的启发式分割算法在GNSS坐标时间序列阶跃探测中的应用[J]. 武汉大学学报 ( 信息科学版), 2019, 44(5): 648-654. doi: 10.13203/j.whugis20170322
YAO Yibin, RAN Qishun, ZHANG Bao. Application of Improved Heuristic Segmentation Algorithm to Step Detection of GNSS Coordinate Time Series[J]. Geomatics and Information Science of Wuhan University, 2019, 44(5): 648-654. doi: 10.13203/j.whugis20170322
Citation: YAO Yibin, RAN Qishun, ZHANG Bao. Application of Improved Heuristic Segmentation Algorithm to Step Detection of GNSS Coordinate Time Series[J]. Geomatics and Information Science of Wuhan University, 2019, 44(5): 648-654. doi: 10.13203/j.whugis20170322

改进的启发式分割算法在GNSS坐标时间序列阶跃探测中的应用

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

国家自然科学基金 41574028

详细信息
    作者简介:

    姚宜斌, 博士, 教授, 主要从事GNSS数据处理研究。ybyao@whu.edu.cn

    通讯作者: 冉启顺, 博士生。RQS0905@126.com
  • 中图分类号: P228

Application of Improved Heuristic Segmentation Algorithm to Step Detection of GNSS Coordinate Time Series

Funds: 

The National Natural Science Foundation of China 41574028

More Information
    Author Bio:

    YAO Yibin, PhD, professor, specializes in GNSS data processing and algorithm research. E⁃mail:ybyao@whu.edu.cn

    Corresponding author: RAN Qishun, PhD candidate. E⁃mail:RQS0905@126.com
  • 摘要: 在启发式分割算法的基础上,引入z检验和标准正态均一性检验(standard normal homogeneity test,SNHT),提出了一种新的阶跃探测法,并将其应用于20个陆态网(Crustal Movement Observation Network ofChina,CMONOC)和10个国际全球卫星导航系统服务(International Global Navigation Satellite System Ser-vice,IGS)参考站近5 a的坐标时间序列。结果显示,对于IGS站东(east,E)、北(north,N)、垂直(up,U)3个方向分量已知阶跃的平均准确探测率分别为68.5%、71.3%、64.8%,且利用阶跃点前后25 d初始拟合残差的平均值之差得到修复后的坐标时间序列保持了良好的连续性。基于剔除粗差后的坐标时间序列,利用坐标时间序列分析软件(coordinate time series analysis software,CATS)估计出的速度显示:对于由已知阶跃与改进算法探测出的阶跃历元的组合方式,所有测站E、N、U方向的估计速度与CMONOC官网发布的速度的平均偏差分别为2.86 mm/a、1.14 mm/a、2.31 mm/a,明显小于由已知阶跃和肉眼判断明显阶跃历元的组合方式获得的平均速度偏差。上述结果表明,改进的启发式分割算法可应用于坐标时间序列的阶跃探测。
  • 图  1  局域研究测站分布图

    Figure  1.  Distribution of Local Research Stations

    图  2  测站阶跃探测情况

    Figure  2.  Step Detection of Stations

    图  3  坐标序列阶跃修复结果

    Figure  3.  Correction Results of Offset in Coordinate Series

    表  1  误差改正模型/策略

    Table  1.   Error Correction Model/Strategy

    误差类型 改正模型/策略
    固体潮 IERS Convention 2010
    海洋负荷潮 otl_FES2004.grid
    极移 IERS Convention 2010
    相对论效应 文献[11]
    卫星和接收机天线相位中心误差 IGS08.atx
    天线相位缠绕 文献[12]
    地球自转 文献[13]
    对流层延迟 随机游走估计(20 cm+2 cm/$\sqrt {\rm{h}} $)
    对流层映射函数 NMF[14]
    电离层延迟 消电离层组合,消除一阶项延迟,忽略高阶项
    距离观测值误差 相位0.003 m,伪距0.6 m
    接收机钟差 作为参数进行估计
    下载: 导出CSV

    表  2  IGS站各分量RMSE/mm

    Table  2.   Each Component's RMSE of IGS Stations/mm

    测站 E方向 N方向 U方向
    AIRA 31.20 6.67 23.08
    BJFS 11.72 5.10 10.08
    CHAN 9.66 5.35 12.12
    GUAO 12.78 3.75 6.62
    LHAZ 15.42 4.30 9.00
    SHAO 8.85 6.18 8.56
    SUWN 14.30 5.54 10.60
    ULAB 12.30 5.86 8.91
    URUM 10.67 6.25 8.37
    WUHN 27.14 6.68 42.37
    均值 15.40 5.57 13.97
    下载: 导出CSV

    表  3  IGS站各分量已知阶跃探测情况/%

    Table  3.   IGS Stations's Detection Result of Known Step for Each Component/%

    探测指标 E方向 N方向 U方向
    准确率 68.5 71.3 64.8
    漏判率 31.5 28.7 35.2
    下载: 导出CSV

    表  4  测站估计速度与发布速度的平均偏差/(mm·a-1)

    Table  4.   Average Deviation Between the Estimated and Release Velocity of All Stations/(mm·a-1)

    阶跃组成 E方向 N方向 U方向
    与ITRF08相比(仅IGS站) 与CMONOC相比(所有测站) 与ITRF08相比(仅IGS站) 与CMONOC相比(所有测站) 与ITRF08相比(仅IGS站) 与CMONOC相比(所有测站)
    已知阶跃+探测结果 2.53 2.86 1.02 1.14 2.62 2.31
    已知阶跃+肉眼判断 3.11 3.33 1.37 1.10 3.21 3.97
    注:探测结果所包含的阶跃历元是快速人工干预(增加漏判、剔除误判)后剩余的历元
    下载: 导出CSV
  • [1] 黄立人. GPS基准站坐标分量时间序列的噪声特性分析[J].大地测量与地球动力学, 2006, 26(2):31-33, 38 http://d.old.wanfangdata.com.cn/Periodical/dkxbydz200602006

    Huang Liren. Noise Properties in Time Series of Coordinate Component at GPS Fiducial Stations[J]. Journal of Geodesy and Geodynamics, 2006, 26(2):31-33, 38 http://d.old.wanfangdata.com.cn/Periodical/dkxbydz200602006
    [2] Gazeaux J, Williams S, King M, et al. Detecting Offsets in GPS Time Series: First Results from the Detection of Offsets in GPS Experiment[J]. Journal of Geophysical Research, 2013, 118(5):2397-2407 http://cn.bing.com/academic/profile?id=80e541ad414ce1c40005c340777f5585&encoded=0&v=paper_preview&mkt=zh-cn
    [3] Teunissen P J G. Testing Theory: An Introduction, Series on Mathematical Geodesy and Postioning[M]. Delft: Delft University Press, 2009
    [4] Bernaola-Galván P, Ivanov P C, Nunes Amaral L A, et al. Scale Invariance in the Nonstationarity of Human Heart Rate[J]. Physical Review Letters, 2001, 87(16): 1-4 doi:  10.1103-PhysRevLett.87.168105/
    [5] 符淙斌, 王强.气候突变的定义和检测方法[J].大气科学, 1992, 16(4):482-493 doi:  10.3878/j.issn.1006-9895.1992.04.11

    Fu Chongbin, Wang Qiang.Definition and Detection the Sudden Change of Climate[J].Chinese Journal of Atmospheric Sciences, 1992, 16(4):482-493 doi:  10.3878/j.issn.1006-9895.1992.04.11
    [6] 封国林, 龚志强, 董文杰, 等.基于启发式分割算法的气候突变检测研究[J].物理学报, 2005, 54(11):5494-5499 http://d.old.wanfangdata.com.cn/Periodical/wlxb200511087

    Feng Guolin, Gong Zhiqiang, Dong Wenjie, et al. Abrupt Climate Change Detection Based on Heuristic Segmentation Algorithm[J]. Acta Physica Sinica, 2005, 54(11):5494-5499 http://d.old.wanfangdata.com.cn/Periodical/wlxb200511087
    [7] 汪丽娜, 陈晓宏, 李粤安, 等.水文时间序列突变点分析的启发式分割方法[J].人民长江, 2009, 40(9): 15-17, 106 doi:  10.3969/j.issn.1001-4179.2009.09.006

    Wang Lina, Chen Xiaohong, Li Yuean, et al. Heuristic Segmentation Method for Change-Point Analysis of Hydrological Time Series[J]. Yangtze River, 2009, 40(9): 15-17, 106 doi:  10.3969/j.issn.1001-4179.2009.09.006
    [8] Ducré-Robitaille J F, Vincent L A, Boulet G. Comparison of Techniques for Detection of Discontinuities in Temperature Series[J]. International Journal of Climatology, 2003, 23(9): 1087-1101 doi:  10.1002/joc.v23:9
    [9] Alexandersson H. A Homogeneity Test Applied to Precipitation Data[J]. Journal of Climatology, 1986, 6(6):661-675 doi:  10.1002/joc.v6:6
    [10] Williams S D P. CATS: GPS Coordinate Time Series Analysis Software[J]. GPS Solutions, 2008, 12(2): 147-153 doi:  10.1007/s10291-007-0086-4
    [11] Kouba J, Pierre Héroux. Precise Point Positioning Using IGS Orbit and Clock Products[J]. GPS Solutions, 2001, 5(2):12-28 doi:  10.1007/PL00012883
    [12] Wu J T, Wu S C, Hajj G A, et al. Effects of Antenna Orientation on GPS Carrier Phase[J]. Manuscripta Geodaetica, 1993, 18(2):91-98 http://cn.bing.com/academic/profile?id=01cde011ff31e1c8def187fbb6771e9a&encoded=0&v=paper_preview&mkt=zh-cn
    [13] 易重海.实时精密单点定位理论与应用研究[D].长沙: 中南大学, 2011 http://cdmd.cnki.com.cn/Article/CDMD-10533-1011177780.htm

    Yi Chonghai. Research on Theory and Application of Real Time Precise Point Positioning[D]. Changsha: Central South University, 2011 http://cdmd.cnki.com.cn/Article/CDMD-10533-1011177780.htm
    [14] Niell A E. Global Mapping Functions for the Atmosphere Delay at Radio Wavelengths[J]. Journal of Geophysical Research Solid Earth, 1996, 101(B2): 3227-3246 doi:  10.1029/95JB03048
    [15] Khaliq M N, Ouarda T B M J. On the Critical Values of the Standard Normal Homogeneity Test (SNHT)[J]. International Journal of Climatology, 2007, 27(5):681-687 doi:  10.1002/(ISSN)1097-0088
    [16] 伍吉仓, 邓康伟, 陈永奇.地心坐标系与站心坐标系中的速度转换及误差传播[J].大地测量与地球动力学, 2005, 25(3):13-18 http://d.old.wanfangdata.com.cn/Periodical/dkxbydz200503003

    Wu Jicang, Deng Kangwei, Chen Yongqi. Velocity Transformation and Error Propagation Between Geocentric Coordinate System and Site-centric Coordinate System[J].Journal of Geodesy and Geodynamics, 2005, 25(3):13-18 http://d.old.wanfangdata.com.cn/Periodical/dkxbydz200503003
  • [1] 马俊, 曹成度, 姜卫平, 周吕.  利用小波包系数信息熵去除GNSS站坐标时间序列有色噪声 . 武汉大学学报 ( 信息科学版), 2021, 46(9): 1309-1317. doi: 10.13203/j.whugis20190353
    [2] 吴昊, 刘超, 赵兴旺.  基于改进累积和控制图的GNSS变形信息的识别与预警 . 武汉大学学报 ( 信息科学版), 2020, 45(10): 1517-1525. doi: 10.13203/j.whugis20190003
    [3] 赵齐乐, 许小龙, 马宏阳, 刘经南.  GNSS实时精密轨道快速计算方法及服务 . 武汉大学学报 ( 信息科学版), 2018, 43(12): 2157-2166. doi: 10.13203/j.whugis20180374
    [4] 刘帅, 贾小林, 孙大伟.  GNSS星载原子钟性能评估 . 武汉大学学报 ( 信息科学版), 2017, 42(2): 277-284. doi: 10.13203/j.whugis20150344
    [5] 陈俊平, 周建华, 严宇, 陈倩, 王彬.  GNSS数据处理时空参数的相关性 . 武汉大学学报 ( 信息科学版), 2017, 42(11): 1649-1657. doi: 10.13203/j.whugis20170278
    [6] 陈冠旭, 刘焱雄, 柳响林, 戴礼文, 范士杰.  船载GNSS探测海洋水汽信息的影响因子分析 . 武汉大学学报 ( 信息科学版), 2017, 42(2): 270-276. doi: 10.13203/j.whugis20150028
    [7] 卢立果, 刘万科, 李江卫.  一种有效的LLL规约算法 . 武汉大学学报 ( 信息科学版), 2016, 41(8): 1118-1124. doi: 10.13203/j.whugis20140484
    [8] 吴寒, 姚宜斌, 陈鹏, 肖勇.  GNSS电离层层析成像算法研究 . 武汉大学学报 ( 信息科学版), 2013, 38(12): 1405-1408.
    [9] 查峰, 许江宁, 李京书, 何泓洋.  一类捷联惯导系统模糊内阻尼算法的改进 . 武汉大学学报 ( 信息科学版), 2013, 38(6): 705-709.
    [10] 章红平, 韩文慧, 黄玲, 耿长江.  地基GNSS全球电离层延迟建模 . 武汉大学学报 ( 信息科学版), 2012, 37(10): 1186-1189.
    [11] 李英冰, 闫景仙, 熊程波, 陈中新.  GNSS观测值的压缩方法研究 . 武汉大学学报 ( 信息科学版), 2012, 37(7): 818-822.
    [12] 闫伟, 欧吉坤, 袁运斌, DagobertoSalazar.  利用广播星历和区域参考网实现实时精密单点定位算法研究 . 武汉大学学报 ( 信息科学版), 2012, 37(10): 1190-1193.
    [13] 冯威, 黄丁发, 严丽, 李萌.  GNSS双频整周关系约束模糊度算法研究 . 武汉大学学报 ( 信息科学版), 2012, 37(8): 945-948.
    [14] 蔡华, 赵齐乐, 孙汉荣, 胡志刚.  GNSS实时数据质量控制 . 武汉大学学报 ( 信息科学版), 2011, 36(7): 820-824.
    [15] 刘志平, 何秀凤, 郭广礼, 查剑锋.  GNSS模糊度降相关算法及其评价指标研究 . 武汉大学学报 ( 信息科学版), 2011, 36(3): 257-261.
    [16] 刘经南, 赵莹, 张小红.  GNSS无线电掩星电离层反演技术现状与展望 . 武汉大学学报 ( 信息科学版), 2010, 35(6): 631-635.
    [17] 孙保琪, 欧吉坤, 盛传贞, 刘吉华.  一种适于Compass周跳探测的三频数据优化组合 . 武汉大学学报 ( 信息科学版), 2010, 35(10): 1157-1160.
    [18] 张飞舟, 杨东凯, 陈嘉, 杨伯钢.  GNSS软件接收机的实现与测试分析 . 武汉大学学报 ( 信息科学版), 2009, 34(5): 577-580.
    [19] 李征航, 陈锴, 刘万科, 王文丽.  GNSS电离层延迟模型的数学统一与方法扩展 . 武汉大学学报 ( 信息科学版), 2007, 32(8): 699-703.
    [20] 陈祥, 杨志强, 田镇, 杨兵, 梁沛.  GA-VMD与多尺度排列熵结合的GNSS坐标时序降噪方法 . 武汉大学学报 ( 信息科学版), 0, 0(0): 0-0. doi: 10.13203/j.whugis20210215
  • 加载中
图(3) / 表(4)
计量
  • 文章访问数:  996
  • HTML全文浏览量:  67
  • PDF下载量:  223
  • 被引次数: 0
出版历程
  • 收稿日期:  2018-02-28
  • 刊出日期:  2019-05-05

改进的启发式分割算法在GNSS坐标时间序列阶跃探测中的应用

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

    国家自然科学基金 41574028

    作者简介:

    姚宜斌, 博士, 教授, 主要从事GNSS数据处理研究。ybyao@whu.edu.cn

    通讯作者: 冉启顺, 博士生。RQS0905@126.com
  • 中图分类号: P228

摘要: 在启发式分割算法的基础上,引入z检验和标准正态均一性检验(standard normal homogeneity test,SNHT),提出了一种新的阶跃探测法,并将其应用于20个陆态网(Crustal Movement Observation Network ofChina,CMONOC)和10个国际全球卫星导航系统服务(International Global Navigation Satellite System Ser-vice,IGS)参考站近5 a的坐标时间序列。结果显示,对于IGS站东(east,E)、北(north,N)、垂直(up,U)3个方向分量已知阶跃的平均准确探测率分别为68.5%、71.3%、64.8%,且利用阶跃点前后25 d初始拟合残差的平均值之差得到修复后的坐标时间序列保持了良好的连续性。基于剔除粗差后的坐标时间序列,利用坐标时间序列分析软件(coordinate time series analysis software,CATS)估计出的速度显示:对于由已知阶跃与改进算法探测出的阶跃历元的组合方式,所有测站E、N、U方向的估计速度与CMONOC官网发布的速度的平均偏差分别为2.86 mm/a、1.14 mm/a、2.31 mm/a,明显小于由已知阶跃和肉眼判断明显阶跃历元的组合方式获得的平均速度偏差。上述结果表明,改进的启发式分割算法可应用于坐标时间序列的阶跃探测。

English Abstract

姚宜斌, 冉启顺, 张豹. 改进的启发式分割算法在GNSS坐标时间序列阶跃探测中的应用[J]. 武汉大学学报 ( 信息科学版), 2019, 44(5): 648-654. doi: 10.13203/j.whugis20170322
引用本文: 姚宜斌, 冉启顺, 张豹. 改进的启发式分割算法在GNSS坐标时间序列阶跃探测中的应用[J]. 武汉大学学报 ( 信息科学版), 2019, 44(5): 648-654. doi: 10.13203/j.whugis20170322
YAO Yibin, RAN Qishun, ZHANG Bao. Application of Improved Heuristic Segmentation Algorithm to Step Detection of GNSS Coordinate Time Series[J]. Geomatics and Information Science of Wuhan University, 2019, 44(5): 648-654. doi: 10.13203/j.whugis20170322
Citation: YAO Yibin, RAN Qishun, ZHANG Bao. Application of Improved Heuristic Segmentation Algorithm to Step Detection of GNSS Coordinate Time Series[J]. Geomatics and Information Science of Wuhan University, 2019, 44(5): 648-654. doi: 10.13203/j.whugis20170322
  • 由于地震、接收机天线更换等已知因素以及一些未知原因,全球卫星导航系统(Global Navigation Satellite System, GNSS)连续运行观测站的点位会发生突变[1],这种突变以阶跃的形式反映在测站的坐标时间序列中。如果阶跃不能得到有效的探测和改正,则GNSS坐标时间序列势必会对速度场估计、坐标参考框架维持和地球动力学研究等带来不利影响。

    文献[2]通过给坐标时间序列中人为添加随机大小的阶跃,分别将目前主流的自动、半自动探测法与直接人工判断法的结果进行对比,结果发现:直接肉眼判断的探测结果往往更好,大多数探测法对于阶跃值较大的历元探测效果良好,但对于小阶跃探测能力较差,并且容易出现误判。如文献[3]利用Neyman-Pearson检验法,对坐标时间序列历元间一次差分构成的新序列进行粗差探测,认为粗差点对应原始序列中的阶跃,该方法不能有效探测小尺度的阶跃;澳大利亚地球科学局提出的GA法,对每个历元将坐标序列分为左右两个子序列,利用线性函数进行拟合,通过对比拟合残差判断阶跃,该方法不能有效判断序列起始和结尾部分的数据中是否存在阶跃[3]。同时,肉眼判别具有较强的主观性和局限性,在噪声水平较高时,无法识别较小的阶跃。

    针对这些问题,本文利用精密单点定位(precise point positioning, PPP)技术解算出的GNSS连续运行观测站多年测站坐标,在启发式分割算法[4-7]的基础上,引入z检验以及气象学领域常用的标准正态均一性检验(standard normal homogeneity test,SNHT)[8-9]探测阶跃。利用坐标时间序列分析软件(coordinate time series analysis software, CATS)[10]估计出的测站速度分别与陆态网(Crustal Movement Observation Network of China, CMONOC)和国际GNSS服务(International GNSS Service, IGS)发布的速度进行对比,根据对应速度间的偏差进一步验证本文方法的可靠性。

    • 本文选取了华北平原上均匀分布的20个CMONOC连续运行观测站和中国及周边10个IGS参考站作为研究对象,时间跨度为2011-01-01—2015-12-31。研究区域内的各测站概略分布如图 1所示。

      图  1  局域研究测站分布图

      Figure 1.  Distribution of Local Research Stations

      为了确保数据解算的准确性和稳定性,舍弃单天观测时段总长小于3 h的观测文件。利用笔者团队研发的精密单点定位程序PASSION对每个测站坐标进行解算,获取单天解。程序解算过程中,利用IGS提供的事后精密星历以及钟差文件来改正卫星的轨道误差和钟误差,测站坐标、接收机钟差和对流层延迟模型化后的残余误差均通过卡尔曼滤波进行估计。表 1给出了本文定位程序中使用的部分误差改正模型和处理策略。

      表 1  误差改正模型/策略

      Table 1.  Error Correction Model/Strategy

      误差类型 改正模型/策略
      固体潮 IERS Convention 2010
      海洋负荷潮 otl_FES2004.grid
      极移 IERS Convention 2010
      相对论效应 文献[11]
      卫星和接收机天线相位中心误差 IGS08.atx
      天线相位缠绕 文献[12]
      地球自转 文献[13]
      对流层延迟 随机游走估计(20 cm+2 cm/$\sqrt {\rm{h}} $)
      对流层映射函数 NMF[14]
      电离层延迟 消电离层组合,消除一阶项延迟,忽略高阶项
      距离观测值误差 相位0.003 m,伪距0.6 m
      接收机钟差 作为参数进行估计

      PPP解算坐标的框架与精密星历的框架一致,本文数据处理采用IGS05/IGS08/IGb08星历,解算结果是对应ITRF05/ITRF08框架下的坐标。由于国际地球自转服务(International Earth Rotation Service, IERS)最新发布的国际地球参考框架ITRF14处理数据的时间段包含了本文研究的时段,且其发布的IGS站速度由整网平差估计得到,同时顾及了年周期、半年周期项以及震后形变,相对于ITRF08框架下仅考虑线性趋势估计而得到的速度更加准确、可靠。所以本文采取的策略是首先利用IGS官网发布的转换参数,将所有坐标统一转换到ITRF08框架下,同时将IGS站在ITRF14框架下的速度转换到ITRF08框架下。然后分别选取各站对应的地心坐标序列中2011年第1个有PPP解算坐标的历元的三维坐标作为站心,将对应的地心坐标转换成站心坐标序列,结果标记为S0。最后对任一测站,分别对其原始的各站心坐标分量应用拉依达准则,将满足式(1)的历元点剔除:

      $$\left| {{X_j} - \bar X} \right| > a\sqrt {\frac{1}{{n - 1}}\mathop { \sum\limits_{j = 1} }\limits^n {{({X_j} - \bar X)}^2}} $$ (1)

      式中,Xj表示某一测站在特定分量中第j个历元的坐标值;X表示该测站在指定分量上所有历元坐标的均值;n表示该分量中测站5 a的有效观测历元总数;a表示中误差的倍数(通常选取3),根据剔除粗差后的坐标时间序列图判断粗差剔除效果,并适当调整其大小。将剔除粗差后的各站站心坐标序列标记为S1

    • 首先在http://sopac-ftp.ucsd.eduasures/网站下载本文选取的10个IGS站在ITRF08框架下的地心坐标序列,并将之转换成该框架下的站心坐标序列(转换方法与§1.1中对应测站相同);然后与本文解算得到的10个IGS站的S1序列对比,将测站各分量上两个时间序列残差的均方根误差(root mean square error, RMSE)作为评定PPP解算质量的依据。

      $${\rm{RMSE}} = \sqrt {\frac{1}{k}\mathop {\mathop \sum\limits_{i = 1} }\limits^k {{\left( {{X_{{\rm{ppp}}, i}} - {X_{{\rm{sopac}}, i}}} \right)}^2}} $$ (2)

      式中,k表示测站在指定分量下的坐标总数;Xppp, i表示PPP解算的第i个历元的站心坐标值;Xsopac, i表示斯克里普斯轨道和永久阵列中心(Scripps Orbit and Permanent Array Center, SOPAC)发布的相同历元的站心坐标值。统计结果如表 2所示。其中,E、N、U分别表示东、北、垂直分量。

      表 2  IGS站各分量RMSE/mm

      Table 2.  Each Component's RMSE of IGS Stations/mm

      测站 E方向 N方向 U方向
      AIRA 31.20 6.67 23.08
      BJFS 11.72 5.10 10.08
      CHAN 9.66 5.35 12.12
      GUAO 12.78 3.75 6.62
      LHAZ 15.42 4.30 9.00
      SHAO 8.85 6.18 8.56
      SUWN 14.30 5.54 10.60
      ULAB 12.30 5.86 8.91
      URUM 10.67 6.25 8.37
      WUHN 27.14 6.68 42.37
      均值 15.40 5.57 13.97

      表 2可以看出,对于所选的10个IGS站,除AIRA、WUHN测站的E分量和U分量外,本文利用PPP解算出的其他测站各坐标分量与SOPAC发布的坐标差值的RMSE均在15.40 mm内,所有IGS站的各分量平均RMSE未超过15.40 mm,其中各站N分量的RMSE最小。总体来说,PPP解算的坐标与SOPAC发布的坐标差异较小,可用于后续分析。

    • 在本节的数据处理中,基于测站各分量序列S1,将去除趋势项、年周期、半年周期的残差序列标准化,并将标准化后的序列作为改进算法的输入项,最终探测出测站各分量坐标序列发生阶跃的历元。

    • 对于测站各分量坐标序列S1采用式(3)进行拟合(不考虑已知阶跃,便于验证改进算法的探测能力),得到拟合残差序列,并将其标准化,标记为R

      $$\begin{array}{*{20}{c}} {y\left( {{t_i}} \right) = a + b{t_i} + c{\rm{sin}}\left( {2{\rm{ \mathsf{ π} }}{t_i}} \right) + d{\rm{cos}}\left( {2{\rm{ \mathsf{ π} }}{t_i}} \right) + }\\ {e{\rm{sin}}\left( {4{\rm{ \mathsf{ π} }}{t_i}} \right) + f{\rm{cos}}\left( {4{\rm{ \mathsf{ π} }}{t_i}} \right) + {v_i}} \end{array}$$ (3)

      式中,a表示参考时刻的位置;b表示线性速率;cd为年周期项系数;ef表示半年周期项系数;ti表示特定分量第i个历元的时间值,为了便于方程的拟合求解,本文使用小数点年表示,即由年积日与该年总天数的比值加上对应年份数值而得;y(ti)表示ti时间测站对应坐标分量下,经预处理得到的站心坐标;vi表示残差。

    • z检验法常用来检验两个数据量较大且服从标准正态分布的随机样本均值的显著性差异。将长度为n的序列R分成两个子序列X1X2,用XiSini分别表示子序列Xi的均值、方差以及序列长度(i=1, 2),构建统计量${z_0} = \frac{{{{\bar X}_1} - {{\bar X}_2}}}{{{{\left( {{s_1}/{n_1} + {s_2}/{n_2}} \right)}^{1/2}}}}$。本文选择显著性水平α=0.05,当|z0| > zα/2时,认为两个子序列的样本均值有显著性差异,即认为对应的序列S1在点n1处可能存在阶跃。

    • 将长度为n的序列R分成两个子序列X1X2,构建统计量T0=n1X12+(n-n1)X22,其中n1表示X1的长度,X1X2分别表示子序列X1X2的均值。当T0大于某个临界值时,即认为对应的序列S1n1处可能存在阶跃。本文选择显著性水平α=0.05,根据序列长度n在文献[15]中选择合适的临界值。

    • 1)对于由n个点构成的序列R,计算第i(1 < i < n)个点左边部分和右边部分子序列的均值ul(i)、ur(i)以及方差sl(i)、sr(i),构建z检验统计量${z_0} = \frac{{{u_l}\left( i \right) - {u_r}\left( i \right)}}{{{{[{s_l}\left( i \right)/{n_l} + {s_r}\left( i \right)/{n_r}]}^{1/2}}}}$。

      2)若|z0| > zα/2,表明序列S1在点i处可能存在位置突变,利用R中点i的左边部分和右边部分子序列执行SNHT检验;反之,i=i+1,重新执行步骤1)。

      3)当点i在SNHT检验中也被判断为阶跃点时,认为对应测站的坐标分量序列S1在点i处存在位置突变,并将R在点i处分割为左右两个子序列;否则,i=i+1,重新执行步骤1)。

      将分割而得的子序列替换成R重复执行上述操作,直至所有的子序列都不可分割或不存在阶跃为止。通常,当子序列的长度小于或等于l1时(l1为最小的子序列长度,常选25 d),不再对其进行分割。通过改变l1z检验以及SNHT检验的显著性水平值,可对序列不同尺度的阶跃进行探测。

    • 本文仅以WUHN、ULAB测站N分量的S1为例,在对二者阶跃探测时,z检验与SNHT检验的显著性水平均选择0.05。自动探测的结果如图 2所示,其中绿色竖线表示正确探测出的阶跃点,橘黄色竖线表示误判的阶跃点,灰色竖线表示漏判的阶跃点,蓝色竖线表示已知阶跃的历元,横轴表示小数点年。

      图  2  测站阶跃探测情况

      Figure 2.  Step Detection of Stations

      图 2可以看出,两站含有较多且难以直接判断的阶跃,采用改进后的方法,最终WUHN测站共探测出10个阶跃点,即绿色和橘黄色竖线标记的历元,从左至右分别是2011.225、2011.408、2011.997、2013.008、2013.488、2013.712、2014.049、2014.537、2014.997、2015.068,其中正确判断和误判的历元各5个,蓝色竖线表示2个已知阶跃点2011.222、2013.879;与WUHN测站类似,ULAB测站总共探测出10个阶跃点,从左至右分别是2011.184、2011.605、2011.874、2012.514、2012.852、2012.959、2013.104、2013.918、2013.173、2014.595,其中正确判断和误判的历元各5个,灰色竖线表示在历元2012.997的漏判点,4个已知阶跃点分别为2011.189、2011.718、2012.855、2013.019。

      对于测站WUHN,改进算法探测出的第1个阶跃点2011.225与已知阶跃点2011.222仅相隔1 d,可认为探测出了1个已知阶跃点。对于测站ULAB,改进算法探测出的第1个阶跃点2011.184与已知阶跃点2011.189仅相隔2 d;第2个阶跃点2011.605与已知阶跃点2011.718虽然相差的时间较长,但相隔的这段时间缺失观测数据,从时间序列图来看属于正确探测的历元;第5个阶跃点2012.852与已知阶跃点2012.855仅相差1 d。所以改进算法对于ULAB测站的N分量序列S1,探测出了3个已知阶跃点。

      统计所有IGS站E、N、U这3个分量对于已知阶跃的平均准确探测率与漏判率,如表 3所示,其中准确率表示探测出的已知阶跃个数与该站已知阶跃总个数的比值,漏判率与准确率的和为1。

      表 3  IGS站各分量已知阶跃探测情况/%

      Table 3.  IGS Stations's Detection Result of Known Step for Each Component/%

      探测指标 E方向 N方向 U方向
      准确率 68.5 71.3 64.8
      漏判率 31.5 28.7 35.2

      上述结果表明,改进算法对于序列的已知阶跃探测情况表现良好。利用所有正确判断和漏判阶跃点前后n0天(通常选取25 d)初始拟合残差的平均值之差修复S1的阶跃,初步修复后的序列保持了良好的连续性,不存在明显的突变,表明本文提出的阶跃探测方法是可行的。其中,WUHN和ULAB两测站N分量坐标序列S1初步修复前后的对比结果如图 3所示。

      图  3  坐标序列阶跃修复结果

      Figure 3.  Correction Results of Offset in Coordinate Series

    • 为了进一步验证改进算法的可靠性,首先基于所有测站各分量的序列S1,假设其噪声组合为白噪声与闪烁噪声[1],将已知阶跃历元分别与改进算法探测出的阶跃历元、肉眼判断明显的阶跃历元组合作为CATS软件的输入项,基于趋势项、年周期、半年周期的函数模型估计其测站速度及精度;然后将IERS发布的IGS站在ITRF14框架下的地心坐标速度(已经转换到ITRF08框架)转换成站心速度[16];同时以与本文处理数据的时间段最接近为原则,从CMONOC官网收集采用GIPSY或GAMIT软件处理得到的对应测站站心速度,以便对比分析。

      统计两种阶跃组合方案估计出的各站速度分别与ITRF14发布的IGS站速度(已转换到ITRF08框架下)、CMONOC发布的测站速度的平均偏差,结果如表 4所示。

      表 4  测站估计速度与发布速度的平均偏差/(mm·a-1)

      Table 4.  Average Deviation Between the Estimated and Release Velocity of All Stations/(mm·a-1)

      阶跃组成 E方向 N方向 U方向
      与ITRF08相比(仅IGS站) 与CMONOC相比(所有测站) 与ITRF08相比(仅IGS站) 与CMONOC相比(所有测站) 与ITRF08相比(仅IGS站) 与CMONOC相比(所有测站)
      已知阶跃+探测结果 2.53 2.86 1.02 1.14 2.62 2.31
      已知阶跃+肉眼判断 3.11 3.33 1.37 1.10 3.21 3.97
      注:探测结果所包含的阶跃历元是快速人工干预(增加漏判、剔除误判)后剩余的历元

      表 4可以看出,对于已知阶跃与改进算法探测出的阶跃历元组合方式,CATS的估计结果显示,所有测站E、N、U方向的估计速度与CMONOC官网发布速度的平均偏差分别为2.86 mm/a、1.14 mm/a、2.31 mm/a;IGS站的估计速度与ITRF14发布的测站速度(已转换到ITRF08下)的平均偏差更小。对于已知阶跃与肉眼判断的阶跃历元的组合方式,估计结果显示,所有测站E、U方向的估计速度与CMONOC官网发布速度的平均偏差分别为3.33 mm/a、3.97 mm/a,明显高于第1种方案估计得到的对应的平均速度偏差,但N方向的平均速度偏差与第1种方案估计得到的平均速度偏差相当,主要是因为原始站心坐标序列中N方向的外符合精度最好,阶跃极少;IGS站E、N、U方向的估计速度与IERS发布的ITRF14框架下的站心速度(已转换到ITRF08下)平均偏差分别为3.11 mm/a、1.37 mm/a、3.21 mm/a,也明显高于第1种方案估计得到的平均速度偏差。

    • 本文基于启发式分割算法的原理,引入z检验和SNHT检验,提出了一种新的阶跃探测方法。该方法在双重统计检验的基础上,进行快速人工干预,能够较为准确地找到序列中存在的阶跃;同时通过改变检验的显著性水平,可对不同尺度的阶跃进行探测,比直接肉眼判断更加可靠,最大限度地降低了主观因素的影响。但需要注意的是,在观测质量较差、噪声水平较大的时段,可能会出现较多误判,不可避免地增加了人工干预的工作量。

      测站坐标时间序列图直观地反映了探测结果,初步修复阶跃后的坐标序列保持了良好的连续性。统计各站通过CATS估计出的速度显示,相对于仅考虑已知阶跃和肉眼判断明显阶跃历元的组合方式,利用改进算法探测出的历元与已知阶跃历元的组合方式估计得到的速度与ITRF14、CMONOC发布的测站速度的平均偏差更小,进一步表明了本文改进的阶跃探测法是一种可行的方法,能够有效地探测坐标序列中存在的阶跃,可应用于GNSS连续运行观测站的速度估计等研究中。

参考文献 (16)

目录

    /

    返回文章
    返回