基于Mallat小波的厦门CORS基准站时间序列分析.pdf

第 11 卷 第 3 期 导航定位学报 Vol.11，No.3 2023 年 6 月 Journal of Navigation and Positioning Jun.，2023 引文格式：杨海鹏,赵利江,省天琛.基于 Mallat 小波的厦门 CORS 基准站时间序列分析J.导航定位学报,2023,11(3):179-184.（YANG Haipeng,ZHAO Lijiang,SHENG Tianchen.Time series analysis of Xiamen CORS reference station based on Mallat waveletJ.Journal of Navigation and Positioning,2023,11(3):179-184.）DOI:10.16547/ki.10-1096.20230324.基于 Mallat 小波的厦门 CORS 基准站时间序列分析 杨海鹏，赵利江，省天琛（青海省基础测绘院，西宁 810000）摘要：为了进一步分析连续运行参考站系统（CORS）的稳定性，提出一种CORS 基准站时间序列分析方法：以福建省厦门市为例，采用马拉特（Mallat）小波分析方法，对厦门区域基准站（XMCORS）的稳定性进行分析；利用小波方差图和小波功率谱图对时间序列中的周期项进行分析，确定周期和分解尺度，然后进行多尺度分解，设置阈值去噪并重构全球卫星导航系统（GNSS）观测数据的时间序列。结果表明，GNSS 观测数据的时间序列具备多周期项特性，软阈值法去噪效果比硬阈值去噪效果更好。关键词：小波分析；全球卫星导航系统；时间序列；周期性；多尺度分解；去噪 中图分类号：P228 文献标志码：A 文章编号：2095-4999(2023)03-0179-06 Time series analysis of Xiamen CORS reference station based on Mallat wavelet YANG Haipeng,ZHAO Lijiang,SHENG Tianchen （Qinghai Institute of Basic Surveying and Mapping,Xining 810000,China）Abstract：In order to further analyze the stability of continuously operating reference stations(CORS),the paper proposed a time series analysis method of CORS:taking Xiaman,Fujian as the example,Mallat wavelet analysis method was used to analyze the stability of XMCORS;the wavelet square difference diagram and wavelet power spectrum were used to analyze the period term in the time series to determine the period and decomposition scale,and the multi-scale decomposition was carried out,then the threshold was set for denoising,and the time series of global navigation satellite system(GNSS)observation were reconstructed finally.Result showed that GNSS time series would have the characteristics of multi-period term,and the denoising effect of soft threshold method could be better than that of hard threshold method.Keywords:wavelet analysis;global navigation satellite system(GNSS);time series;periodicity;multiscale decomposition;denoising 0 引言 全球卫星导航定位系统（global navigation satellite system，GNSS）地基增强系统被广泛应用在测绘、无人机、自动驾驶等行业。省级、市级的连 续 运 行 参 考 站 系 统（continuously operating reference stations，CORS）在城市精细管理、智能城市中提供时空基准。CORS 基准站的稳定分析对地区的地壳运动、坐标框架维护至关重要。为了更好地分析福建省连续运行卫星定位服务系统（厦门区域）基准站（XMCORS）稳定性，提升实时和事后定位服务质量，通过建立 XMCORS 基准站的时间序列和速度场，保障厦门市现代大地坐标框架的时效性，分析厦门 CORS 基准站的坐标时间序列，获得长期、季节性、随机性等多种地球物理量并分析其机理1。信号分析采用时-频域分析方法，通过傅里叶变化探测平稳信号的频率特征。然而，地球科学中 收稿日期：2022-03-02 第一作者简介：杨海鹏（1979），男，河北永清人，本科学历，工程师，研究方向为卫星导航精密定位技术和应用。180 导航定位学报 2023 年 6 月 的大多数信号（如 GNSS 时间序列）都不是平稳信号，而是具有趋势项、周期性等特点。文献2提出一种时-频多分辨功能的小波分析（wavelet analysis）法，能够更好地研究地球科学时间序列问题；该方法能清晰地揭示出隐藏在时间序列中的多种变化周期，以及反映不同时间尺度的变化趋势。文献3提出了小波多尺度分析的方法，提供小波函数的构造，通过小波变换对信号的不同频率成分进行分解，为信号滤波、信噪分离和特征提取提供了有效途径。基准站稳定性分析包括数据质量分析、获取基准站坐标时间序列、分析坐标时间序列的噪声模型和获取基准站水平速度场和垂向速度场，综合评价基准站稳定性情况。本文在研究小波分析理论的基础上，采用小波分析应用于厦门 CORS基准站分析，重点对时间序列周期、噪声等特性进行研究。1 小波分析方法 小波分析的核心思想与傅里叶变换相同，都是利用一组函数去表示或者无限逼近信号，这一组函数称为小波函数系。类似于向量基的概念，小波函数系是通过小波基函数平移和伸缩构成的一族函数。设小波基函数为?（?表示平方可积的实数空间），则小波函数系为4 (),a btbta baaa=其中R1，0（1）在 GNSS 时间序列分析中?为尺度因子，反映小波的周期长度；b 为平移因子，反映时间上的平移，在本文研究中选择 morlet 小波函数。实数空间中任意给定的函数?，其连续小波变换是?与小波基函数?,?的内积，即()()(),()()da btbW a bf ttf ttaa=1R（2）式中：?,?为小波变换系数；()tba是()tba的复共轭函数；t 是时间（自变量）。实际观测信号都是离散非连续的，将积分改为求和即可。1.1 小波方差和小波功率谱 小波方差4定义为小波变换系数的平方值在平移因子 b 值域上积分，即 ()()Var|,|d=2aW a bb（3）小波方差Var?是尺度因子?的函数，由式（3）可知，小波方差能反映信号能量随尺度?的分布；小波方差随?变化过程图形，称为小波方差图。根据函数特性，小波方差图反映出信号在不同种尺度下的能量波动，即反映信号中的主周期性。由于 morlet 小波函数是个复数形式，GNSS 时间序列?经过小波变换之后得到的小波系数?,?也是复数形式：有实部?,?和虚部?,?；依据实部、虚部计算小波系数的振幅|?,?|和相位()(),tan(),W a bW a b1IR，小波功率谱4定义为振幅的平方，有 ()()()|,|,W a bW a bW a b=+222RI（4）1.2 小波多尺度分解和重构 GNSS 测站空间坐标可以分为 3 个一维分量形式（如 E、N、U），因此本文采用一维小波多尺度分析。设有?空间的子空间序列：?012VVV（5）尺度函数?和小波函数?经过伸缩、平移得到（?为实数空间）：?2?;?,?（6）?2?;?,?（7）式（6）、式（7）构成（?、?），是?的正交基函数。?和?满足关系式（二尺度方程）：()()()nxl nxn=22Z（8）()()()nxh nxn=22Z（9）式中：?是低通滤波器；?是高通滤波器，且?1?1?。由此，得到 GNSS 信号的一维多尺度分解（比如高程方向）为 ()()()()JJJjjnkkkknkjkf xcxncxdx=+01ZZZ （10）式（10）最右边的第一项()()JJJkkkfxcx=Z是()f x在尺度 J下的逼近，J值越大逼近效果越差，其系数为()f x的离散细节；右边第二项()=jgx()jjkkkdxZ是?的频率在2?到2?之间的连续细节部分。第 3 期 杨海鹏，等.基于 Mallat 小波的厦门 CORS 基准站时间序列分析 181 ()()jjkknjjkknccl nkdch nk=1122ZZ（11）式中：?是尺度系数；?是小波系数。式（11）即为马拉特（Mallat）小波多尺度分解和重构算法。若已知双尺度方程中的滤波器系数，可由?快速计算出小波多尺度的逼近和细节5。1.3 多尺度去噪 小波多尺度分析是将信号进行分解、处理后再重构，能够将信号中某些特征分层显示6。对于GNSS时间序列资料的分析，通常考虑长周期运动、同震变形以及震后变形。GNSS单站、单分量坐标序列可用数学模型拟合，即()()()()()()sincossin cos=+12242giiiiniijigjijy tabtctdtetftg H tTv（12）式中：?和?表示变量时间和地震发生时刻，单位为年；?、?表示初始位置和形变速率；?、?、?、?表示周年和半周年系数；?表示同震位移；?为阶跃函数；?为残差值。GNSS 时间序列中存在噪声，去噪操作如下：1）小波分解。使用 morlet 小波函数对 GNSS时间序列进行多层小波分解。2）阈值设定。对不同层小波细节系数用合适的阈值处理。3）小波重构。重建阈值调整后的小波系数，重构新的GNSS时间序列。式（12）操作的核心是步骤 2）阈值设定；本文研究硬阈值（式（13）和软阈值（式（14）2 种去噪效果，有：,if(),otherwiseiiTixxTx=00（13）()max(,)=10 1TiTxxx（14）式中 T 为阈值，根据软件 GAMIT/GLOBK 解算的GNSS 数据 3 倍中误差设定 T7，即 T=3。最后，采用信噪比8-10（signal-noise ratio，SNR）指标来反映小波去噪的效果，即 sNRzlog()=1010pSp（15）式 中：s()/=2npfnn为 原 始 信 号 功 率；()z()/=?2npf nf nn为噪声功率，其中()?f n为去噪后的信号，()f n为原始信号。信噪比越大，去噪效果越好。2 实验与结果分析 采用高精度的 GNSS 数据处理软件 GAMIT/GLOBK 17.4 对 XMCORS 的 6 个基准站（如图 1所示）2013 年 1 月 1 日2016 月 12 月 31 日期间的数据进行处理。通过引入中国周边 11 个国际GNSS 服务组织（International GNSS Service，IGS）核心站，将 IGS 站点坐标和基准框架引入到基线解算的过程。GMAIT 基线解算的重复性较高。图 1 厦门 XMCORS 基准站分布 表 1 中?表示基线长度，基线重复性北向（N）为 1.55 mm+1.4610-9?，东向（E）为 1.66 mm+1.4610-9?，垂向（U）为 5.61 mm+3.3410-9?。GLOBK 平差结果显示：各基准站 N 方向的归一化均方根（normalized root mean square，NRMS）值小于 1.66 mm，平均 NRMS 值为 0.84 mm；E 方向的 NRMS 值都小于 1.44 mm，平均 NRMS 值为0.75 mm；U 方向的 NRMS 值都小于 0.96 mm，平均NRMS值 为0.60 mm；同 时，根 据GAMIT/GLOBK 软件解算结果，确定=0.73 mm。基于上述研究，采用小波分析 GNSS 时间序列流程如图 2 所示7,11。其中|?|?表示小波功率谱，DWT 表示离散小波分解（discrete wavelet transform，DWT）。由于地壳运动导致 GNSS 时间序列存在长期趋势项，因此在小波分析前需要进行“去趋势”操作得到零均值的时间序列（如图 3 182 导航定位学报 2023 年 6 月 表 1 GAMIT/GLOBK 解算结果 mm 精度指标 北向 东向 垂向 基线重复性 1.55+1.4610-9D 1.66+1.4610-9D 5.61+3.3410-9D 平均 NRMS 值 0.84 0.75 0.60 图 2 GNSS 时间序列小波分析流程 图 3 XMJY 去趋势项时间序列 所示），能够看出基准站的时间序列存在明显周期波动，为后续谱分析和小波分析提供了数据支撑。依据上述理论对厦门 CORS 基准站时间序列数据进行小波分析。图 4 展示了其中的 XMJY 站小波方差和小波功率谱。小波方差和功率谱能够反映数据中的周期性质，图中在 12 a 间出现显著波峰，这说明 GNSS 时间序列中存在显著的周年项。式（12）的数学模型中也体现出，GNSS 时间序列就是线性项和周期项的叠加信号。GNSS时间序列中同时存在多个周期因素，依据图 4 显示，存在 4 个明显波动，并依据文献8研究结果，尝试将 GNSS 时间序列进行 4 层小波分解（如图5 所示），并根据小波去噪原理对每一层的高频部分进行去噪，之后再重构时间序列。表 2 统计了利用 CATS（create and analyze time series）软件对 6 个厦门 CORS 基准站坐标三分量所有时间序列进行小波去噪处理前后的噪声振幅估值和信噪比值。图 6 和表 2 显示软阈值和硬阈值去噪效果，结果显示软阈值法去噪效果更好，将 GNSS 时间序列中大部分有色噪声过滤掉了。以往研究发现，图 4 XMJY 站 GPS 时间序列小波方差和功率谱 第 3 期 杨海鹏，等.基于 Mallat 小波的厦门 CORS 基准站时间序列分析 183 图 5 GNSS 时间序列的原始序列和 4 层分解的高频、低频序列 图 6 GNSS 时间序列与硬阈值法、软阈值法去噪对比 表 2 去闪烁噪声前后 GNSS 时间序列统计信息基准站 基准站 状态 垂向白 噪声振 幅/mm 垂向闪烁噪声振幅/mm SNRXMJY 去噪前原始数据 3.9 17.1 硬阈值去噪后 3.0 6.0 7.8软阈值去噪后 2.5 1.5 6.0（续）基准站状态 垂向白 噪声振 幅/mm 垂向闪烁噪声振幅/mm SNRXMZC去噪前原始数据 4.3 10.5 硬阈值去噪后 3.1 3.1 9.8软阈值去噪后 2.9 2.5 10.9XMJM去噪前原始数据 3.8 16.0 硬阈值去噪后 2.8 5.9 8.3软阈值去噪后 2.7 2.0 10.2XMDF去噪前原始数据 4.6 21.0 硬阈值去噪后 3.1 5.4 9.8软阈值去噪后 2.2 1.1 10.5XMDM去噪前原始数据 3.1 18.6 硬阈值去噪后 2.9 7.4 7.3软阈值去噪后 2.4 2.2 8.6XMDD去噪前原始数据 3.8 13.9 硬阈值去噪后 3.5 6.3 7.7软阈值去噪后 2.0 1.8 11.3在评估 GNSS 基准站速度场时，有色噪声会干扰结果精度，根据本文研究结果发现软阈值去噪可 184 导航定位学报 2023 年 6 月 以解决该问题。另一方面，时间序列中的噪声与GNSS 建站方式有很大关系，对于研究建站区域的地质以及浅地表的地球物理因素，需要保留真实的有色噪声，则建议采用硬阈值去噪法。3 结束语 小波分析方法可以对时域信号的局部特征进行放大，能将信号分解到不同的尺度上进行多尺 度分析。本文介绍了小波分析的基本理论，并将小波技术运用在厦门 CORS 基准站时间序列分析中，对基准站坐标时间序列的周期性质和多尺度分解、去噪及重构进行研究。结果表明，软阈值比硬阈值去噪效果更好。2 种去噪方法的应用场景不同：1）估计 GNSS 基准站速度场时，采用软阈值去噪法；2）研究 GNSS 建站区域的地质以及浅地表的地球物理现象，则建议采用硬阈值去噪法。参考文献 1 李婧.“陆态网络”基准站坐标时间序列变化特性分析D.信息工程大学,2013.2 MORLET J,ARENS G,FOURGEAU I,et al.Wave propagation and sampling theoryJ.Geoplysics,1982,47:203-236.3 MALLAT S.A theory of multiresolution signal decomposition:the wavelet representationJ.IEEE Transactions on Pattern Analysis and Machine Intelligence,1989,11(7):674-693.4 TORRENCE C,COMPO G P.A practical guide to wavelet analysisJ.Bulletin of the American Meteorological Society,1998,79(1):61-78.5 汪海洪.小波多尺度分析在地球重力场中的应用研究D.武汉大学,2005.6 黄声享,刘经南.GPS 变形监测系统中消除噪声的一种有效方法J.测绘学报,2002(2):104-107.7 范朋飞.高精度 GPS 站点坐标时间序列分析与应用D.长安大学,2013.8 陈强,黄声享,王韦.小波去噪效果评价的另一指标J.测绘信息与工程,2008(5):13-14.9 马俊,曹成度,姜卫平,等.利用小波包系数信息熵去除 GNSS 站坐标时间序列有色噪声J.武汉大学学报(信息科学版),2021,46(9):1309-1317.10 周亚,王立峰,张思慧,等.IGS 连续运行参考站高程时间序列功率谱分析J.太赫兹科学与电子信息学报,2014,12(1):103-107.11 田亮,孙付平.基于 GPS 测站坐标残差序列的小波工具应用与分析J.测绘工程,2013,22(1):44-46,50.