留言板

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

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

小波多尺度分解和奇异谱分析在GNSS站坐标时间序列分析中的应用

戴海亮 孙付平 姜卫平 肖凯 朱新慧 刘婧

戴海亮, 孙付平, 姜卫平, 肖凯, 朱新慧, 刘婧. 小波多尺度分解和奇异谱分析在GNSS站坐标时间序列分析中的应用[J]. 武汉大学学报 ● 信息科学版, 2021, 46(3): 371-380. doi: 10.13203/j.whugis20190107
引用本文: 戴海亮, 孙付平, 姜卫平, 肖凯, 朱新慧, 刘婧. 小波多尺度分解和奇异谱分析在GNSS站坐标时间序列分析中的应用[J]. 武汉大学学报 ● 信息科学版, 2021, 46(3): 371-380. doi: 10.13203/j.whugis20190107
DAI Hailiang, SUN Fuping, JIANG Weiping, XIAO Kai, ZHU Xinhui, LIU Jing. Application of Wavelet Decomposition and Singular Spectrum Analysis to GNSS Station Coordinate Time Series[J]. Geomatics and Information Science of Wuhan University, 2021, 46(3): 371-380. doi: 10.13203/j.whugis20190107
Citation: DAI Hailiang, SUN Fuping, JIANG Weiping, XIAO Kai, ZHU Xinhui, LIU Jing. Application of Wavelet Decomposition and Singular Spectrum Analysis to GNSS Station Coordinate Time Series[J]. Geomatics and Information Science of Wuhan University, 2021, 46(3): 371-380. doi: 10.13203/j.whugis20190107

小波多尺度分解和奇异谱分析在GNSS站坐标时间序列分析中的应用

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

国家自然科学基金 41674042

国家自然科学基金 41804018

详细信息
    作者简介:

    戴海亮,硕士,主要从事导航时空基准研究。David_109@163.com

    通讯作者: 孙付平,博士,教授。sun.fp@163.com
  • 中图分类号: P228

Application of Wavelet Decomposition and Singular Spectrum Analysis to GNSS Station Coordinate Time Series

Funds: 

The National Natural Science Foundation of China 41674042

The National Natural Science Foundation of China 41804018

More Information
    Author Bio:

    DAI Hailiang, master, specializes in navigation space-time datum. E-mail: David_109@163.com

    Corresponding author: SUN Fuping, PhD, professor. E-mail: sun.fp@163.com
  • 摘要: 为了有效地提取GNSS(global navigation satellite system)站坐标时间序列中的有用信息,提高坐标时间序列的建模精度,提出一种小波多尺度分解与奇异谱分析相结合的非线性运动建模方法,并利用全球11个测站20年(1999―2018年)的GPS(global positioning system)垂向坐标时间序列对所提方法进行了验证。首先,通过小波分解将坐标时间序列分解到不同尺度上; 然后,对分解后的各层高频部分和低频部分进行奇异谱分析; 最后,通过叠加合成得到原始坐标时间序列的拟合值,并对所提方法的拟合效果进行评估。结果表明,与单纯的奇异谱分析方法相比,所提方法能够更加准确地从含噪声的有限尺度时间序列中提取趋势和周期等有用信息,降低了部分周期项如季节周期项、月周期项被当作噪声剔除的概率,并且建模精度有26%的提高。
  • 图  1  小波多分辨率分解示意图

    Figure  1.  Schematic Diagram of Wavelet Multi‐resolution Decomposition

    图  2  BJFS站垂向时间序列及其特征值贡献率

    Figure  2.  Vertical Time Series and Its Eigenvalue Contribution Ratio of BJFS Station

    图  3  前6阶特征值的重构成分

    Figure  3.  Reconstruction Components of the First 6 Order Eigenvalues

    图  4  前6阶重构成分的频率谱

    Figure  4.  Frequency Spectrums of the First 6 Order Reconstruction Components

    图  5  SSA的拟合结果及残差

    Figure  5.  Fitting Results and Residuals of Singular Spectrum Analysis

    图  6  残差序列的频率谱

    Figure  6.  Frequency Spectrum of the Residual Series

    图  7  小波分解后的各层低频部分和高频部分

    Figure  7.  Low-Frequency Part and High-Frequency Part of Each Layer After Wavelet Decomposition

    图  8  近似系数a3的特征值贡献率及拟合效果

    Figure  8.  Eigenvalue Contribution Ratio and Fitting Effect of Approximation Coefficient a3

    图  9  细节系数d3的特征值贡献率及拟合效果

    Figure  9.  Eigenvalue Contribution Ratio and Fitting Effect of Detail Coefficient d3

    图  10  细节系数d2的特征值贡献率及拟合效果

    Figure  10.  Eigenvalue Contribution Ratio and Fitting Effect of Detail Coefficient d2

    图  11  细节系数d1的特征值贡献率及拟合效果

    Figure  11.  Eigenvalue Contribution Ratio and Fitting Effect of Detail Coefficient d1

    图  12  小波分解和SSA的拟合结果及残差

    Figure  12.  Fitting Results and Residuals of Wavelet Decomposition and Singular Spectrum Analysis

    图  13  残差序列的频率谱

    Figure  13.  Frequency Spectrum of the Residual Series

    表  1  前14阶特征值贡献率统计

    Table  1.   Statistics of the First 14 Order Eigenvalue Contribution Ratio

    阶次 重构成分 特征值贡献率/% 阶次 重构成分 特征值贡献率/%
    1 RRC1 32.30 8 RRC8 1.20
    2 RRC2 30.45 9 RRC9 1.13
    3 RRC3 4.71 10 RRC10 1.06
    4 RRC4 4.68 11 RRC11 0.93
    5 RRC5 4.19 12 RRC12 0.83
    6 RRC6 1.55 13 RRC13 0.77
    7 RRC7 1.36 14 RRC14 0.77
    下载: 导出CSV

    表  2  常用小波基函数及其特性

    Table  2.   Common Wavelet Basis Functions and Their Characteristics

    小波基函数 支撑长度 消失矩阶数 对称性 特点
    Haar小波 1 1 对称 在时域上是不连续的,频率的局部性差
    dbN小波 2N N 近似对称 光滑性随着N的增加而增加,性能较好
    symN小波 2N N 近似对称 能够减少相位失真,更适合于图像处理领域
    Meyer小波 有限长度 对称 对信号进行分析会产生失真现象
    bior小波 重构:2Nr+1
    分解:2Nd+1
    Nr-1 不对称 具有线性相位性,在信号与图像的重构中应用比较普遍
    下载: 导出CSV

    表  3  不同分解层数下拟合精度对比

    Table  3.   Fitting Accuracy Comparison of Different Decomposition Layers

    拟合精度 分解层数
    2 3 4 5 6
    RMSE/mm 2.15 1.88 1.80 1.80 1.80
    MAE/mm 1.68 1.49 1.43 1.44 1.44
    下载: 导出CSV

    表  4  两种方法的垂直方向拟合精度对比

    Table  4.   Comparison of Fitting Accuracy Between Two Models in the Vertical Direction

    测站名称 纬度 SSA/mm 小波SSA/mm
    RMSE MAE RMSE MAE
    ADIS 9.2°N 2.59 1.98 1.90 1.44
    TUVA 8.3°S 2.94 2.17 2.23 1.71
    NAUR 0.3°S 2.82 2.09 2.09 1.61
    IISC 13.1°N 2.73 2.00 1.98 1.47
    WUHN 30.5°N 2.90 2.31 2.11 1.66
    STR1 35.2°S 2.50 1.81 1.83 1.38
    BJFS 39.6°N 2.48 1.93 1.88 1.49
    MAC1 54.3°S 2.12 1.65 1.65 1.27
    HOFN 64.2°N 2.13 1.67 1.66 1.30
    KELY 66.6°N 2.64 2.05 1.82 1.43
    MAW1 67.4°S 2.19 1.68 1.46 1.14
    下载: 导出CSV
  • [1] 马俊, 姜卫平, 邓连生, 等. GPS坐标时间序列噪声估计及相关性分析[J]. 武汉大学学报·信息科学版, 2018, 43 (10): 1 451-1 457 doi:  10.13203/j.whugis20160543

    Ma Jun, Jiang Weiping, Deng Liansheng, et al. Estimation Method and Correlation Analysis for Noise in GPS Coordinate Time Series[J]. Geomatics and Information Science of Wuhan University, 2018, 43 (10): 1 451-1 457 doi:  10.13203/j.whugis20160543
    [2] 马俊, 姜卫平, 周晓慧, 等. 联合小波和方差分量估计方法分析中国IGS测站时间序列变化特征[J]. 武汉大学学报·信息科学版, 2018, 43 (4): 629-636 doi:  10.13203/j.whugis20150731

    Ma Jun, Jiang Weiping, Zhou Xiaohui, et al. Coordinate Time Series of IGS Station in China Based on the Combination of Wavelet Spectral and Variance Component Estimation[J]. Geomatics and Information Science of Wuhan University, 2018, 43 (4): 629-636 doi:  10.13203/j.whugis20150731
    [3] Tregoning P, Watson C. Atmospheric Effects and Spurious Signals in GPS Analyses[J]. Journal of Geophysical Research Solid Earth, 2009, 114(B9), DOI: 10.1029/2009JB006344
    [4] 朱文耀, 符养, 李彦. GPS高程导出的全球高程振荡运动及季节变化[J]. 中国科学: 地球科学, 2003, 33(5): 470-481

    Zhu Wenyao, Fu Yang, Li Yan. Global Elevation Oscillation Movement and Seasonal Change Derived from GPS Elevation[J]. Science in China: Earth Science, 2003, 33(5): 470-481
    [5] Amiri-Simkooei A R, Tiberius C C J M, Teunissen P J G. Assessment of Noise in GPS Coordinate Time Series: Methodology and Results[J]. Journal of Geophysical Research Solid Earth, 2007, 112(B7): 1-19
    [6] 姜卫平, 马一方, 邓连生, 等. 毫米级地球参考框架的建立方法与展望[J]. 测绘地理信息, 2016, 41(4): 1-6

    Jiang Weiping, Ma Yifang, Deng Liansheng, et al. Establishment of mm-Level Terrestrial Reference Frame and Its Prospect[J]. Journal of Geomatics, 2016, 41(4): 1-6
    [7] Bogusz J, Klos A. On the Significance of Periodic Signals in Noise Analysis of GPS Station Coordinates Time Series[J]. GPS Solutions, 2016, 20(4): 655-664 doi:  10.1007/s10291-015-0478-9
    [8] 张鹏, 蒋志浩, 秘金钟, 等. 我国GPS跟踪站数据处理与时间序列特征分析[J]. 武汉大学学报·信息科学版, 2007, 32(3): 251-254 http://ch.whu.edu.cn/article/id/1848

    Zhang Peng, Jiang Zhihao, Bei Jinzhong, et al. Data Processing and Time Series Analysis for GPS Fiducial Stations in China[J]. Geomatics and Information Science of Wuhan University, 2007, 32(3): 251-254 http://ch.whu.edu.cn/article/id/1848
    [9] van Dam T M, Wahr J, Milly P C D, et al. Crustal Displacement due to Continent Water Loading[J]. Geophysical Research Letters, 2001, 28(4): 651‑654 doi:  10.1029/2000GL012120
    [10] 王敏, 沈正康, 董大南. 非构造形变对GPS连续站位置时间序列的影响和修正[J]. 地球物理学报, 2005, 48(5): 1 045-1 052

    Wang Min, Shen Zhengkang, Dong Danan. Effects of Non-tectonic Crustal Deformation on Continuous GPS Position Time Series and Correction to Them[J]. Chinese Journal of Geophysics, 2005, 48(5): 1 045-1 052
    [11] 张诗玉, 钟敏, 唐诗华, 等. 我国GPS基准站地壳垂直形变的大气负荷效应[J]. 武汉大学学报·信息科学版, 2006, 31(12): 1 090-1 093 http://ch.whu.edu.cn/article/id/2624

    Zhang Shiyu, Zhong Min, Tang Shihua, et al. Atmospheric Load Effect of Crustal Vertical Deformation at GPS Reference Stations in China[J]. Geomatics and Information Science of Wuhan University, 2006, 31(12): 1 090-1 093 http://ch.whu.edu.cn/article/id/2624
    [12] Jiang Weiping, Li Zhao, van Dam T, et al. Comparative Analysis of Different Environment Loading Methods and Their Impacts on the GPS Height Time Series[J]. Journal of Geodesy, 2013, 87(7): 687-703 doi:  10.1007/s00190-013-0642-3
    [13] Wyatt M, Kravtsov S, Tsonis A. Atlantic Multidecadal Oscillation and Northern Hemisphere's Climate Variability[J]. Climate Dynamics, 2012, 38(5): 929-949
    [14] Chen Q, Dam T V, Sneeuw N, et al. Singular Spectrum Analysis for Modeling Seasonal Signals from GPS Time Series[J]. Journal of Geodynamics, 2013, 72(12): 25-35
    [15] Zabalza J, Ren J, Zheng W, et al. Singular Spectrum Analysis for Effective Feature Extraction in Hyperspectral Imaging[J]. IEEE Geoscience and Remote Sensing Letters, 2014, 11(11): 1 886-1 890 doi:  10.1109/LGRS.2014.2312754
    [16] Kondrashov D, Berloff P. Stochastic Modeling of Decadal Variability in Ocean Gyres[J]. Geophysical Research Letters, 2015, 42(5): 1 543-1 553 doi:  10.1002/2014GL062871
    [17] 郭金运, 高文宗, 于红娟, 等. 基于奇异谱分析的静态相对重力观测重力固体潮提取[J]. 地球物理学报, 2018, 61(10) : 3 889-3 902

    Guo Jinyun, Gao Wenzong, Yu Hongjuan, et al. Gravity Tides Extracted from Relative Gravimetric Data with Singular Spectrum Analysis[J]. Chinese Journal of Geophysics, 2018, 61(10): 3 889-3 902
    [18] 文鸿雁. 基于小波理论的变形分析模型研究[D]. 武汉: 武汉大学, 2004

    Wen Hongyan. Research on Deformation Analysis Model Based on Wavelet Transform Theory[D]. Wuhan: Wuhan University, 2004
    [19] Ghaderpour E, Ince E S, Pagiatakis S. Least-Squares Cross-Wavelet Analysis and Its Applications in Geophysical Time Series[J]. Journal of Geodesy, 2018, 92(10): 1 223-1 236 doi:  10.1007/s00190-018-1156-9
    [20] 范朋飞. 高精度GPS站点坐标时间序列分析与应用[D]. 西安: 长安大学, 2013

    Fan Pengfei. Analysis and Application of High-Precision Coordinate Time Series of GPS Site[D].Xi'an: Chang'an University, 2013
    [21] 王解先, 连丽珍, 沈云中. 奇异谱分析在GPS站坐标监测序列分析中的应用[J]. 同济大学学报(自然科学版), 2013, 41(2): 282-288

    Wang Jiexian, Lian Lizhen, Shen Yunzhong. Application of Singular Spectral Analysis to GPS Station Coordinate Monitoring Series[J]. Journal of Tongji University(Natural Science), 2013, 41(2): 282-288
    [22] 周茂盛, 郭金运, 沈毅, 等. 基于多通道奇异谱分析的GNSS坐标时间序列共模误差的提取[J]. 地球物理学报, 2018, 61(11): 4 383-4 395

    Zhou Maosheng, Guo Jinyun, Shen Yi, et al. Extraction of Common Mode Errors of GNSS Coordinate Time Series Based on Multi-channel Singular Spectrum Analysis[J]. Chinese Journal of Geophysics, 2018, 61(11): 4 383- 4 395
  • [1] 吴浩, 卢楠, 邹进贵, 郭世泰.  GNSS变形监测时间序列的改进型3σ粗差探测方法 . 武汉大学学报 ● 信息科学版, 2019, 44(9): 1282-1288. doi: 10.13203/j.whugis20170338
    [2] 赵丹宁, 高蕊, 雷雨.  利用小波分解改进极移预报模型 . 武汉大学学报 ● 信息科学版, 2019, 44(12): 1797-1801. doi: 10.13203/j.whugis20180139
    [3] 刘立龙, 陈雨田, 黎峻宇, 田祥雨, 贺朝双.  活跃期区域电离层总电子短期预报及适用性分析 . 武汉大学学报 ● 信息科学版, 2019, 44(12): 1757-1764. doi: 10.13203/j.whugis20180145
    [4] 刘子维, 张晓彤, 张锐, 江颖, 韦进, 张坤.  连续重力观测异常模式的多分辨率识别算法 . 武汉大学学报 ● 信息科学版, 2018, 43(6): 840-846, 942. doi: 10.13203/j.whugis20160173
    [5] 王笑蕾, 张勤, 张双成.  基于EMD和WD联合算法的GPS水汽时间序列的周期性振荡分析 . 武汉大学学报 ● 信息科学版, 2018, 43(4): 620-628. doi: 10.13203/j.whugis20150596
    [6] 曹建农, 王平禄, 董昱威.  高分辨率遥感影像上居民地自动提取方法 . 武汉大学学报 ● 信息科学版, 2014, 39(7): 831-837.
    [7] 魏二虎, 李智强, 龚光裕, 张帅.  极移时间序列模型的拟合与预测 . 武汉大学学报 ● 信息科学版, 2013, 38(12): 1420-1424.
    [8] 裴媛媛, 廖明生, 王寒梅.  时间序列SAR影像监测堤坝形变研究 . 武汉大学学报 ● 信息科学版, 2013, 38(3): 266-269.
    [9] 何敏, 何秀凤.  利用时间序列干涉图叠加法监测江苏盐城地区地面沉降 . 武汉大学学报 ● 信息科学版, 2011, 36(12): 1461-1465.
    [10] 牟忠凯, 隋立芬, 张清华, 甘雨.  顾及线性化模型误差补偿的卡尔曼滤波算法 . 武汉大学学报 ● 信息科学版, 2011, 36(9): 1073-1076.
    [11] 赵红蕊, 唐中实, 李小文.  非线性不适定遥感反演中正则化参数的定量确定 . 武汉大学学报 ● 信息科学版, 2008, 33(6): 577-580.
    [12] 姚宜斌, 施闯.  IGS测站的非线性变化研究 . 武汉大学学报 ● 信息科学版, 2007, 32(5): 423-426.
    [13] 张鹏, 蒋志浩, 秘金钟, 党亚民.  我国GPS跟踪站数据处理与时间序列特征分析 . 武汉大学学报 ● 信息科学版, 2007, 32(3): 251-254.
    [14] 刘飞, 王新洲, 张鹏林, 余旭.  基于Lebesgue积分理论的时间序列信号中的背景重建 . 武汉大学学报 ● 信息科学版, 2006, 31(5): 462-465.
    [15] 李桂苓, 万剑华, 陶华学.  用差商代替导数的非线性最小二乘估计 . 武汉大学学报 ● 信息科学版, 2001, 26(2): 118-121.
    [16] 张晓东, 李德仁, 蔡东翔, 马洪超.  à trous小波分解在边缘检测中的应用 . 武汉大学学报 ● 信息科学版, 2001, 26(1): 29-33.
    [17] 程建权, 黄经南.  一种基于时序数据的动态聚类分析方法 . 武汉大学学报 ● 信息科学版, 1998, 23(3): 194-198.
    [18] 邵巨良, 李德仁.  小波变换及其在影像表示中的应用和其方向选择性算法的改进 . 武汉大学学报 ● 信息科学版, 1993, 18(2): 1-9.
    [19] 徐培亮.  应用时间序列方法作大坝变形预报 . 武汉大学学报 ● 信息科学版, 1988, 13(3): 23-31.
    [20] 徐培亮.  非线性函数的协方差传播公式 . 武汉大学学报 ● 信息科学版, 1986, 11(2): 92-99.
  • 加载中
图(13) / 表(4)
计量
  • 文章访问数:  85
  • HTML全文浏览量:  14
  • PDF下载量:  36
  • 被引次数: 0
出版历程
  • 收稿日期:  2019-05-29
  • 刊出日期:  2021-03-05

小波多尺度分解和奇异谱分析在GNSS站坐标时间序列分析中的应用

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

    国家自然科学基金 41674042

    国家自然科学基金 41804018

    作者简介:

    戴海亮,硕士,主要从事导航时空基准研究。David_109@163.com

    通讯作者: 孙付平,博士,教授。sun.fp@163.com
  • 中图分类号: P228

摘要: 为了有效地提取GNSS(global navigation satellite system)站坐标时间序列中的有用信息,提高坐标时间序列的建模精度,提出一种小波多尺度分解与奇异谱分析相结合的非线性运动建模方法,并利用全球11个测站20年(1999―2018年)的GPS(global positioning system)垂向坐标时间序列对所提方法进行了验证。首先,通过小波分解将坐标时间序列分解到不同尺度上; 然后,对分解后的各层高频部分和低频部分进行奇异谱分析; 最后,通过叠加合成得到原始坐标时间序列的拟合值,并对所提方法的拟合效果进行评估。结果表明,与单纯的奇异谱分析方法相比,所提方法能够更加准确地从含噪声的有限尺度时间序列中提取趋势和周期等有用信息,降低了部分周期项如季节周期项、月周期项被当作噪声剔除的概率,并且建模精度有26%的提高。

English Abstract

戴海亮, 孙付平, 姜卫平, 肖凯, 朱新慧, 刘婧. 小波多尺度分解和奇异谱分析在GNSS站坐标时间序列分析中的应用[J]. 武汉大学学报 ● 信息科学版, 2021, 46(3): 371-380. doi: 10.13203/j.whugis20190107
引用本文: 戴海亮, 孙付平, 姜卫平, 肖凯, 朱新慧, 刘婧. 小波多尺度分解和奇异谱分析在GNSS站坐标时间序列分析中的应用[J]. 武汉大学学报 ● 信息科学版, 2021, 46(3): 371-380. doi: 10.13203/j.whugis20190107
DAI Hailiang, SUN Fuping, JIANG Weiping, XIAO Kai, ZHU Xinhui, LIU Jing. Application of Wavelet Decomposition and Singular Spectrum Analysis to GNSS Station Coordinate Time Series[J]. Geomatics and Information Science of Wuhan University, 2021, 46(3): 371-380. doi: 10.13203/j.whugis20190107
Citation: DAI Hailiang, SUN Fuping, JIANG Weiping, XIAO Kai, ZHU Xinhui, LIU Jing. Application of Wavelet Decomposition and Singular Spectrum Analysis to GNSS Station Coordinate Time Series[J]. Geomatics and Information Science of Wuhan University, 2021, 46(3): 371-380. doi: 10.13203/j.whugis20190107
  • 随着GNSS(global navigation satellite system)观测数据的不断丰富和数据处理技术的不断发展,分析测站非线性运动规律可以获得各种地球物理现象和季节性影响,如果对这部分残差规律进行建模拟合从而消除其影响,那么可以进一步提高GPS(global positioning system)站坐标的精度[1]。尽管目前国际地球参考框架(international terrestrial reference frame,ITRF)中基准站的历元坐标和速度场已达到了毫米级,但几乎所有IGS(international GNSS service)站坐标时间序列(尤其是垂直方向)都呈现显著的非线性运动趋势,且振幅可达1~2 cm[2-5]。因此,地面点相对于参考框架的运动无法用线性速度进行完整描述,仅基于线性速率的模型来维持参考框架会有一定的局限性,而且其精度只能达到厘米级[6-7]

    要精确表达站点坐标,通常采用两种方式描述非线性运动:(1)不考虑引起测站坐标非线性变化的各种物理机制,只根据坐标时间序列本身的运动趋势建模;(2)从产生形变的物理机理入手,分析各因素影响,从而对测站坐标进行改正。国内外对IGS站坐标时间序列的分析已有不少相关研究,也取得较大进展。如文献[8]对中国IGS站坐标时间序列进行频谱分析,发现高程方向周年特性表现明显;文献[9-12]证实了大气负荷、水文负载、非潮汐海洋负载等负载变化是引起GPS测站垂向位移的主要因素。由于引起测站坐标时间序列非线性变化的各种机制比较复杂,目前,国际上还未建立包含多种机制影响的非线性变化理论改正模型,所以本文倾向于对站点本身的运动信息进行建模。

    奇异谱分析(singular spectrum analysis,SSA)是时间序列分析中一种常用的方法,在气候学[13-14]、测量学[15]和海洋学[16]等领域应用广泛,它的优越性在于无需任何先验信息和假设条件,能稳定识别和强化周期信号,把最可预报的分量聚集到若干个时间序列中,以此选择若干有意义的分量重建序列,降低噪声影响。因此,特别适合分析有周期震荡的时间序列数据[17]。根据贡献率较大的前几项来对时间序列进行重构,能有效地将半年及半年以上的周期项提取出来,而短周期项如季节周期项、月周期项容易被当作噪声忽略掉。因此,需要用一定方法削弱这种影响。

    为更好地从噪声信号中分离出时间序列的周期性、长期性变化信号等有用信息,提高建模精度,本文提出一种小波多尺度分解和SSA相结合的方法提取坐标时间序列中有用信息的本文思路,通过对全球11个测站20年(1999―2018年)的GPS垂向坐标时间序列进行实验分析,验证了利用本文方法能够有效地提取时间序列中的有用信息,提高GPS坐标时间序列的精度。

    • 多分辨率分析是小波分析的重要概念,又称为多尺度分析。它能将信号在不同分辨率上进行分解,得到低频部分和高频部分。其中,低频部分能够反映信号的概貌;而高频部分又叫细节部分,能够刻画信号的细节信息。小波多分辨率分析主要是通过空间分解,将信号分解在子空间列中,方便进行信号分析[18-20],其分解示意图如图 1所示。

      图  1  小波多分辨率分解示意图

      Figure 1.  Schematic Diagram of Wavelet Multi‐resolution Decomposition

      图 1S表示原始时间序列,A1、A2、A3表示小波分解后各层时间序列的低频部分,D1、D2、D3则表示小波分解后各层时间序列的高频部分。分解后的各部分与原始信号具有以下关系:S=A3+D3+D2+D1,以此类推,若要进行进一步的分解,可以把低频部分A3分解成低频部分A4和高频部分D4。因此,可以看出小波变换是几乎无损的信号分析方法。

    • 对给定的一维时间序列x=(x1x2x3xn),标准的SSA算法主要由如下4个步骤构成:

      1)构造时滞矩阵。选择适当的窗口长度$ L $,$ 1<L<N/2 $,一般情况下,$ L $的选取不宜超过整个数据长度N的1/3,如果根据事先经验大致确定数据的周期特征,则$ L $一般取周期的公倍数。根据向量x构建轨迹矩阵$ \mathit{\boldsymbol{X}} $为:

      $$ \mathit{\boldsymbol{X}}=({x}_{ij}{)}_{i, j=1}^{L, K}=\left[\begin{array}{cccc}{x}_{1}& {x}_{2}& \cdots & {x}_{K}\\ {x}_{2}& {x}_{3}& \cdots & {x}_{K+1}\\ ⋮& ⋮& \ddots & ⋮\\ {x}_{L}& {x}_{L+1}& \cdots & {x}_{N}\end{array}\right] $$ (1)

      式中,K=N$ - $L+1;$ \mathit{\boldsymbol{X}} $为$ M\times K $阶时滞矩阵,且$ \mathit{\boldsymbol{X}} $中副对角线上的元素相等,即对于$ \mathit{\boldsymbol{X}} $中的元素$ {x}_{i, j} $,有$ {x}_{i, j}={x}_{i-1, j+1} $,因此,$ \mathit{\boldsymbol{X}} $是一个Hankel矩阵。

      2)奇异值分解(singular value decomposition,SVD)。对时滞矩阵$ \mathit{\boldsymbol{X}} $进行奇异值分解:

      $$ \mathit{\boldsymbol{X}}=\mathit{\boldsymbol{U}}\cdot {\mathit{\pmb{\Lambda }}}^{\frac{1}{2}}\cdot {\mathit{\boldsymbol{V}}}^{\mathrm{T}} $$ (2)

      式中,$ {\mathit{\boldsymbol{U}}}_{1}, {\mathit{\boldsymbol{U}}}_{2}\cdots {\mathit{\boldsymbol{U}}}_{d} $分别为对应的特征向量;$ \mathit{\boldsymbol{\Lambda }}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}[{\lambda }_{1}, {\lambda }_{2}\cdots {\lambda }_{d}] $为矩阵$ \mathit{\boldsymbol{X}}{\mathit{\boldsymbol{X}}}^{\mathrm{T}} $和矩阵$ {\mathit{\boldsymbol{X}}}^{\mathrm{T}}\mathit{\boldsymbol{X}} $的前$ d $($ d=\mathrm{m}\mathrm{i}\mathrm{n}\left\{L, K\right\} $)个最大非零特征值构成的对角矩阵;$ {\mathit{\boldsymbol{\Lambda }}}^{\frac{1}{2}} $表示时滞矩阵$ \mathit{\boldsymbol{X}} $的奇异值。由此推知:

      $$ \left\{\begin{array}{l}{\mathit{\boldsymbol{X}}}_{i}=\sqrt{{\lambda }_{i}}{\mathit{\boldsymbol{U}}}_{i}{\mathit{\boldsymbol{V}}}_{i}^{\mathrm{T}}\\ {\mathit{\boldsymbol{V}}}_{i}={\mathit{\boldsymbol{X}}}^{\mathrm{T}}{\mathit{\boldsymbol{U}}}_{i}/\sqrt{{\lambda }_{i}}\end{array}\right. $$ (3)

      式中,$ \sqrt{{\lambda }_{i}} $为时滞矩阵$ \mathit{\boldsymbol{X}} $的奇异值,$ \sqrt{{\lambda }_{1}}\ge \sqrt{{\lambda }_{2}}\ge \cdots \ge \sqrt{{\lambda }_{d}}\ge 0 $为奇异谱。则时滞矩阵$ \mathit{\boldsymbol{X}} $的SVD可以表示为:

      $$ \mathit{\boldsymbol{X}}={\mathit{\boldsymbol{X}}}_{1}+{\mathit{\boldsymbol{X}}}_{2}+\cdots +{\mathit{\boldsymbol{X}}}_{d} $$ (4)

      3)分组。将矩阵$ {\mathit{\boldsymbol{X}}}_{i} $的下标$ \left\{\mathrm{1, 2}\cdots d\right\} $分成$ M $个互不相交的集合$ {I}_{1}, {I}_{2}\cdots {I}_{M} $,设$ I=\left[{i}_{1}, {i}_{2}\cdots {i}_{p}\right] $,那么与集合$ I $相关的矩阵$ {\mathit{\boldsymbol{X}}}_{I} $就可以表示为$ {\mathit{\boldsymbol{X}}}_{I}={\mathit{\boldsymbol{X}}}_{i1}+{\mathit{\boldsymbol{X}}}_{i2}+\cdots +{\mathit{\boldsymbol{X}}}_{ip} $,轨迹矩阵就可以表示成:

      $$ \mathit{\boldsymbol{X}}={\mathit{\boldsymbol{X}}}_{I1}+{\mathit{\boldsymbol{X}}}_{I2}+\cdots +{\mathit{\boldsymbol{X}}}_{IM} $$ (5)

      一般说来,在SVD中,每个矩阵$ {\mathit{\boldsymbol{X}}}_{I} $对轨迹矩阵$ \mathit{\boldsymbol{X}} $的贡献度与它的特征值有关,其关系可以用下面的式子表示:

      $$ {\eta }_{I}=\sum\limits _{i\in I}{\lambda }_{i}/\sum\limits _{i=1}^{L}{\lambda }_{i} $$ (6)

      4)对角平均化。对角平均化的目的就是把第3)步分解得到的矩阵$ {\mathit{\boldsymbol{X}}}_{IM} $重新转换成为长度为$ N $的新时间序列,称之为重建成分(reconstruction component,RC),所有RC之和等于原序列。假设$ \mathit{\boldsymbol{z}}=({z}_{1}, {z}_{2}\cdots {z}_{N}) $为$ \mathit{\boldsymbol{z}} $经过对角平均化得到的时间序列,则对角平均化的公式可表示为:

      $$ {z}_{n}=\left\{\begin{array}{l}\frac{1}{n}\sum\limits _{j=1}^{n}{z}_{j, n-j+1}, 1\le n\le L\\ \frac{1}{L}\sum\limits _{j=1}^{L}{z}_{j, n-j+1}, L<n<K\\ \frac{1}{N-n+1}\sum\limits _{j=n-K+1}^{L}{z}_{j, n-j+1}\begin{array}{c}\begin{array}{c}, \end{array}\end{array}K\le n\le N\end{array}\right. $$ (7)

      所有重建成分RC叠加后的和与原始序列$ x $相同,即:

      $$ \mathit{\boldsymbol{x}}={z}_{1}+{z}_{2}+\cdots +{z}_{N}=\sum\limits _{n=1}^{N}{z}_{n} $$ (8)

      截取前K个贡献大的成分近似表示原序列,则有[21]

      $$ \widehat{\mathit{\boldsymbol{x}}}={z}_{1}+{z}_{2}+\cdots +{z}_{K}=\sum\limits _{k=1}^{K}{z}_{k}, k=\mathrm{1, 2}\cdots N $$ (9)

      上述步骤中,有两个参数需自行选择,分别是窗口长度$ L $和重构阶次$ K $。

    • 首先,对原始坐标时间序列S进行小波分解和重构,得到低频部分和多层的高频部分。具体过程如下,对$ S=\left\{{S}_{1}, {S}_{2}\cdots {S}_{n}\right\} $进行小波多尺度分解和重构,得到下式:

      $$ S={A}_{B}+{D}_{1}+{D}_{2}+{D}_{3}+\cdots +{D}_{N} $$ (10)

      式中,$ {A}_{B} $为低频部分,$ {A}_{B}=\left\{{a}_{B1}, {a}_{B2}\cdots {a}_{Bn}\right\} $为S中低频信号的重构结果;$ {D}_{N} $为多层的高频部分,$ {D}_{1}=\left\{{d}_{11}, {d}_{12}\cdots {d}_{1n}\right\}, {D}_{2}=\left\{{d}_{21}, {d}_{22}\cdots {d}_{2n}\right\}\cdots {D}_{N}=\left\{{d}_{N1}, {d}_{N2}\cdots {d}_{Nn}\right\} $分别为原始坐标时间序列S中第1层到第N层的高频信号的重构结果。

      那么,$ {t}_{i} $时刻的原始坐标时间序列S可以表示为:

      $$ {S}_{i}={a}_{Bi}+{d}_{1i}+{d}_{2i}+\cdots +{d}_{Ni} $$ (11)

      式中,$ {a}_{Bi} $为低频信号在$ {t}_{i} $时刻的值;$ {d}_{1i}, {d}_{2i}\cdots {d}_{Ni} $分别为各层高频信号在$ {t}_{i} $时刻的值。

      利用SSA方法对$ {A}_{B}=\left\{{a}_{B1}, {a}_{B2}\cdots {a}_{Bn}\right\} $和$ {D}_{i}=\left\{{d}_{i1}, {d}_{i2}\cdots {d}_{in}\right\}, $i=1,2,3…N分别进行奇异谱分解和重建,得到低频信号在$ {t}_{i} $时刻的拟合值为$ {\widehat{a}}_{Bi} $,各层高频信号在$ {t}_{i} $时刻的拟合值分别为$ {\widehat{d}}_{1i}, {\widehat{d}}_{2i}\cdots {\widehat{d}}_{Ni} $。则原始坐标时间序列S在$ {t}_{i} $时刻的拟合序列为:

      $$ {\widehat{S}}_{i}={\widehat{a}}_{Bi}+{\widehat{d}}_{1i}+{\widehat{d}}_{2i}+\cdots +{\widehat{d}}_{Ni} $$ (12)
    • 实验采用的数据来自IGS网站提供的全球测站坐标时间序列文件,大部分测站时间序列包含1999―2018年近20 a的数据,采样间隔为一周。§1.2使用SSA算法对坐标时间序列数据进行处理时,有窗口长度L和重构阶次K两个参数需要自行选择。

      由于测站坐标时间序列中存在的已知周期有周年周期和半年周期,且本文所用数据的时间分辨率为一周,因此,在利用SSA方法提取坐标时间序列的信息时可以选择已知周期的最小公倍数52为窗口长度[22]

      参数K主要根据奇异值的贡献率来确定,若太小,则表示后面部分信号被当作噪声剔除掉;若太大,则会使得部分噪声被当作信号而提取出来。因此,确定窗口长度$ L $与重构阶次K是SSA的关键。

      鉴于篇幅,本文以BJFS站垂直方向时间序列为例,根据SSA方法对其进行分析,并将前14阶特征值的贡献率进行了统计,如图 2表 1所示。

      图  2  BJFS站垂向时间序列及其特征值贡献率

      Figure 2.  Vertical Time Series and Its Eigenvalue Contribution Ratio of BJFS Station

      表 1  前14阶特征值贡献率统计

      Table 1.  Statistics of the First 14 Order Eigenvalue Contribution Ratio

      阶次 重构成分 特征值贡献率/% 阶次 重构成分 特征值贡献率/%
      1 RRC1 32.30 8 RRC8 1.20
      2 RRC2 30.45 9 RRC9 1.13
      3 RRC3 4.71 10 RRC10 1.06
      4 RRC4 4.68 11 RRC11 0.93
      5 RRC5 4.19 12 RRC12 0.83
      6 RRC6 1.55 13 RRC13 0.77
      7 RRC7 1.36 14 RRC14 0.77

      图 2(b)表 1可以看出,有多对特征值近似相等,一般来说,近似相等的特征值表示周期与振幅相同的重构分量。其中RRC1RRC2的特征值贡献率较大,分别为32.30%和30.45%;RRC3RRC4RRC5贡献率次之,分别为4.71%、4.68%和4.19%;从第6阶开始,各重构序列的特征值贡献率较低,且阶次越高,贡献率越小。

      本文选择将特征值贡献率较大的前6阶特征值的重构成分进行快速傅里叶变换(fast Fourier transform, FFT)提取周期,结果如图 3图 4所示。

      图  3  前6阶特征值的重构成分

      Figure 3.  Reconstruction Components of the First 6 Order Eigenvalues

      图  4  前6阶重构成分的频率谱

      Figure 4.  Frequency Spectrums of the First 6 Order Reconstruction Components

      根据图 4可知,重构序列RRC1RRC2的周期相同,将其合并可以得到周期为1 a、振幅为4.76 mm的周期项;RRC5的周期分别为0.5 a(与RRC3的周期相同)和9 a(与RRC4的周期相同),将其合并得到周期为0.5 a、振幅为0.95 mm的周期项以及周期为9 a、振幅为0.88 mm的周期项;而RRC6对应的周期为0.3 a(即季节周期项)且振幅影响较小。上述结果与直接利用FFT提取的结果相吻合。

      综上可知,在利用SSA对BJFS站垂向时间序列进行分析时,原序列的趋势项和周期项主要分布在前5个重构序列中。因此,本文认为前5阶RC为主要信息成分,可以选择这5个特征值对时间序列进行重构,拟合结果如图 5所示。

      图  5  SSA的拟合结果及残差

      Figure 5.  Fitting Results and Residuals of Singular Spectrum Analysis

      图 5(a)可以看出,重构后的时间序列大致表现出了测站垂直方向时间序列的走势,但是拟合效果不太理想,图 5(b)中残差振幅在3 mm左右,而且还存在部分周期规律。

      在对残差序列进行频谱分析时发现,时间序列中还存在一些半年以下的周期项,如图 6所示。

      图  6  残差序列的频率谱

      Figure 6.  Frequency Spectrum of the Residual Series

      因此,基于SSA的非线性运动建模方法能有效地将半年及半年以上的周期项提取出来,而短周期项如季节周期项、月周期项则容易被忽略。

    • 在对原始坐标时间序列进行小波分解的过程中,还存在着小波基函数和分解层数如何确定的问题。对几种常用的小波基函数及其特性进行了统计,结果如表 2所示。

      表 2  常用小波基函数及其特性

      Table 2.  Common Wavelet Basis Functions and Their Characteristics

      小波基函数 支撑长度 消失矩阶数 对称性 特点
      Haar小波 1 1 对称 在时域上是不连续的,频率的局部性差
      dbN小波 2N N 近似对称 光滑性随着N的增加而增加,性能较好
      symN小波 2N N 近似对称 能够减少相位失真,更适合于图像处理领域
      Meyer小波 有限长度 对称 对信号进行分析会产生失真现象
      bior小波 重构:2Nr+1
      分解:2Nd+1
      Nr-1 不对称 具有线性相位性,在信号与图像的重构中应用比较普遍

      表 2可知,dbN小波具有良好的时间序列分析特性,在处理坐标时间序列数据方面优势明显。因此,本文一般选择dbN小波进行小波分解。

      理论上,随着分解层数的增加,测站坐标时间序列被划分得越来越细,各层信号的平稳性和平滑性也越来越好。由于分解过程中存在计算上的误差,以及分解层数越多,会带来计算时间的增加,相应的处理效率也会随之下降,因此,应选择合适的分解层数。

      以BJFS站垂直方向时间序列为例,在不同分解层数(2~6层)时,本文方法的均方根误差(root mean square error,RMSE)和平均绝对误差(mean absolute error,MAE)如表 3所示。根据表 3可知,当分解层数超过3层时,新本文方法的拟合精度随着分解层数的增加基本上没有太大的变化;当分解层数为3~4层时,本文方法的之间差异不大。考虑到测站坐标时间序列数据量很大,在保证处理效率的基础上,本文进行小波分解时,分解层数一般选择3层。

      表 3  不同分解层数下拟合精度对比

      Table 3.  Fitting Accuracy Comparison of Different Decomposition Layers

      拟合精度 分解层数
      2 3 4 5 6
      RMSE/mm 2.15 1.88 1.80 1.80 1.80
      MAE/mm 1.68 1.49 1.43 1.44 1.44

      根据表 2表 3确定的小波基函数以及分解层数,实验利用具有正交性和高度紧支撑性的db4小波将BJFS站垂向时间序列分解3层得到各层低频信息和高频信息,如图 7所示。

      图  7  小波分解后的各层低频部分和高频部分

      Figure 7.  Low-Frequency Part and High-Frequency Part of Each Layer After Wavelet Decomposition

      a3+d1+d2+d3之和即为原序列,a3为小波分解后的低频部分,称为近似系数,反映出了时间序列的概貌,表现出较强的规律性(主要包含了趋势项和周期项);d1、d2、d3分别为小波分解后的各层高频部分,称为细节系数,表现出的随机项成分较多。按照§2.1中的特征值选取方法,对小波分解后的低频部分a3和高频部分d1、d2、d3分别选择贡献率较大的前5阶特征值进行序列重构,结果如图 8~图 11所示。

      图  8  近似系数a3的特征值贡献率及拟合效果

      Figure 8.  Eigenvalue Contribution Ratio and Fitting Effect of Approximation Coefficient a3

      图  9  细节系数d3的特征值贡献率及拟合效果

      Figure 9.  Eigenvalue Contribution Ratio and Fitting Effect of Detail Coefficient d3

      图  10  细节系数d2的特征值贡献率及拟合效果

      Figure 10.  Eigenvalue Contribution Ratio and Fitting Effect of Detail Coefficient d2

      图  11  细节系数d1的特征值贡献率及拟合效果

      Figure 11.  Eigenvalue Contribution Ratio and Fitting Effect of Detail Coefficient d1

      图 8~图 11中可以看出,经过小波3层分解后的各层时间序列相对平稳了许多。从近似系数a3和细节系数d3中可以看出,重构序列和原始序列在大部分时间段基本上实现重合,而在细节系数d1和d2中表现出来的随机项成分较多,周期项较少。因此,各特征值贡献率差异不大,相应的拟合效果要小于其他层。将各层时间序列的拟合结果进行叠加得到最终的拟合结果,并将其与原始测站坐标时间序列进行比较,效果如图 12所示。

      图  12  小波分解和SSA的拟合结果及残差

      Figure 12.  Fitting Results and Residuals of Wavelet Decomposition and Singular Spectrum Analysis

      图 12(a)中可以看出,经过小波多尺度分解和SSA相结合的方法进行拟合后的时间序列不仅能够反映出测站坐标时间序列的走势,而且两者之间的差异明显减小。从剩余残差信息图 12(b)中可以看出,本文方法拟合后的残差振幅在2 mm左右,拟合精度有一定程度的提高。

      在对残差序列进行频谱分析时发现,一些短周期项如季节周期项、月周期项的影响明显减小,如图 13所示。

      图  13  残差序列的频率谱

      Figure 13.  Frequency Spectrum of the Residual Series

      与单纯的SSA方法相比,本文方法能更有效地对各层时间序列的周期项进行提取,减小部分短周期项如季节周期项、月周期项容易被忽略的影响。

    • 本文在全球范围内低、中、高纬度各选取了3~4个测站,利用§2.1和§2.2中的方法进行分析,通过计算剩余残差序列的MAE与RMSE,对两种方法在测站垂向时间序列的拟合效果进行了评估,结果如表 4所示。

      表 4  两种方法的垂直方向拟合精度对比

      Table 4.  Comparison of Fitting Accuracy Between Two Models in the Vertical Direction

      测站名称 纬度 SSA/mm 小波SSA/mm
      RMSE MAE RMSE MAE
      ADIS 9.2°N 2.59 1.98 1.90 1.44
      TUVA 8.3°S 2.94 2.17 2.23 1.71
      NAUR 0.3°S 2.82 2.09 2.09 1.61
      IISC 13.1°N 2.73 2.00 1.98 1.47
      WUHN 30.5°N 2.90 2.31 2.11 1.66
      STR1 35.2°S 2.50 1.81 1.83 1.38
      BJFS 39.6°N 2.48 1.93 1.88 1.49
      MAC1 54.3°S 2.12 1.65 1.65 1.27
      HOFN 64.2°N 2.13 1.67 1.66 1.30
      KELY 66.6°N 2.64 2.05 1.82 1.43
      MAW1 67.4°S 2.19 1.68 1.46 1.14

      剩余残差序列的RMSE和MAE分别为:

      $$ \left\{\begin{array}{l}{R}_{\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}=\sqrt{\frac{\sum\limits _{i=1}^{n}({Y}_{i}-{\widehat{Y}}_{i}{)}^{2}}{n}}\\ {M}_{\mathrm{M}\mathrm{A}\mathrm{E}}=\frac{1}{n}\sum\limits _{i=1}^{n}|{Y}_{i}-{\widehat{Y}}_{i}|\end{array}\right. $$ (13)

      式中,$ {\widehat{Y}}_{i} $是真实值$ {Y}_{i} $的拟合值;$ n $是样本容量或者数据个数。

      根据表 4图 5图 12可以看出,在测站垂直方向上,利用小波多尺度分解和SSA相结合的方法对坐标时间序列进行拟合时,其MAE与RMSE皆明显小于仅基于SSA方法的MAE与RMSE。

      经过计算,与SSA方法的拟合结果相比,本文方法的RMSE和MAE分别降低了约26.5%和25.5%。因此,可以得出与上述分析相同的结论,即在相同特征值的情况下,小波分解和SSA的拟合效果要优于仅基于SSA的拟合效果。

      由于经过小波分解后,坐标时间序列与原始序列相比在频率成分上相对单一、平滑,并且原始序列中的趋势项和周期项主要集中在低频部分中,而随机项则大部分存在于高频部分中。

      基于SSA将贡献率较大的前几项特征值进行时间序列重建的方法,会使一部分短周期项被当作噪声剔除掉。而本文方法可以将时间序列进行分层,在每一层上根据所选的特征值重建序列,在一定程度上减少了部分短周期项被当作噪声剔除的概率,从而进一步提高了建模精度。可以看出,采用小波分解和SSA相结合的方法进行非线性运动建模有很好的自适应性,并且建模精度较高。

    • 首先,本文提出一种小波多尺度分解和SSA相结合的非线性运动建模方法,通过小波分解将GPS坐标时间序列分解为各层高频和低频部分;然后,利用SSA对各层时间序列进行分析,截取贡献率较大的前几个特征值重构序列,并将各层时间序列的拟合结果进行叠加得到最终的拟合结果;最后,对本文方法的拟合结果进行了分析与评估,得出如下结论。

      1)与单纯的SSA方法相比,小波多尺度分解和SSA相结合的非线性运动建模方法能有效地融合两种方法的优点。本文方法将原始序列分解成频率成分相对单一、平滑的各层高频部分和低频部分,在每一层上根据所选的特征值重建序列,在一定程度上减少了部分周期项(如季节周期项、月周期项)被当作噪声剔除的概率。结果表明,本文方法的建模精度提高了约26%。

      2)小波分解虽然具有较大的优势,即可以将测站坐标时间序列分解到不同的尺度上,在不同尺度上分别进行SSA建模,但是要注意选择合适的分解层数。考虑到数据处理效率和精度要求,分解层数一般定为3层即可。

      3)小波分解和SSA相结合的方法对坐标时间序列建模具有很好的自适应性,由于测站所处地理位置不同,坐标时间序列所表现出来的规律也不尽相同,各测站前K阶特征值贡献率的大小就存在差异。因此,不宜选取相同的特征值个数对所有的测站进行处理。如何合理地选择特征值个数是今后需要进一步探讨的方向。

参考文献 (22)

目录

    /

    返回文章
    返回