关键词：GNSS坐标时间序列 缺失数据 小波变换 信号提取

Dyadic wavelet transform and signal extraction of GNSS coordinate time series with missing data

JI Kunpu , SHEN Yunzhong

Corresponding author: SHEN Yunzhong,E-mail:yzshen@tongji.edu.cn.

Abstract: The GNSS position time series are often analyzed by using traditional dyadic wavelet transform, which requires that the time series must be complete. However, missing data inevitably occur in the GNSS position time series due to a variety of causes. In order to extract signals from the incomplete position time series, a modified dyadic wavelet transform algorithm is developed and the corresponding formulas are derived in this paper based on the principle that missing data can be reproduced by its wavelet coefficients. The equivalence between new algorithm and zero-padding algorithm is proved, which indicates that the zero-padding algorithm is essentially a least squares minimum norm solution. Finally, the real position time series of 27 based stations from Crustal Movement Observation Network of China (CMONOC) and simulated data are adopted to verify the validation of the new algorithm, the results show that the difference between the signals extracted by new algorithm and interpolation algorithm is small, with the differences of mean medium errors of 27 base stations are only 2.01%(North), 0.54%(East), 1.26%(Up) and the mean ratios of variance for difference of two signals to the variance for two signals are only 1.16%(North), 0.54%(East), 1.62%(Up).

Key words: GNSS coordinate time series missing data wavelet transform signals extraction

应用全球导航定位系统(GNSS)的连续跟踪数据解算的全球与区域基准站坐标已经积累了20多年的数据序列；利用这些坐标序列数据，可进行大量地壳运动变化的相关地学研究。由于受测站外部环境和观测技术本身的误差以及参考框架等因素影响，GNSS坐标序列中叠加了各类信号和噪声，表现出明显的季节性变化 [1] ，尤其在高程方向上更为明显 [2-4] 。季节性信号严重影响测站坐标序列的速度估计 [5] 和地球参考框架ITRF(international terrestrial reference frame)实现的精度 [6] ，因此，必须对季节性信号进行精确估计。传统的谐波模型 [7-9] 认为季节性信号是常振幅、常相位的，并用正余弦谐波表示。然而，已有研究表明，真实的季节性信号由于大气压、地表质量负载等因素呈显著的非线性变化 [10-11] ，其振幅和相位都是随时间变化的，传统的谐波模型不能很好地反映这种非线性变化。近年来，许多学者针对时变季节性信号的提取进行了相关研究。文献[12]提出了一种半参数模型提取季节性信号，但因其理论和计算太过复杂未能得到广泛应用；文献[13]采用Kalman滤波描述季节性信号的统计变化，但该模型仅适用于随机游走噪声，而实际的GNSS坐标序列主要包含闪烁噪声；文献[14]引入奇异谱分析来提取季节性信号，但奇异谱分析的滞后窗口选取具有主观性，不合适的滞后窗口甚至会得到错误的估计结果。

小波分析是20世纪80年代后期发展起来的一个新兴数学分支，因其具有优良的时频局部分析的能力，被广泛地应用到大地测量数据处理中 [15-17] 。作为一种时频域分析方法，小波分析能够将时间序列几乎无损地分解成低频部分和高频部分，再将分解后的信号映射到不同尺度上进行研究，以此得到信号在不同分辨率下的信息 [18] 。文献[19-20]分别基于小波多分辨率分析和小波包分析提取了变形监测数据序列的特征项; 文献[21-22]利用小波分析提取了GNSS坐标序列中季节性形变；文献[23]利用小波熵和小波变换系数模极大值组合法去除GNSS坐标序列中的白噪声和闪烁噪声。小波分析要求完整的时间序列数据 [24] ，然而GNSS站坐标序列必定存在缺值。传统的解决思路是通过事先插值补全缺失部分 [25-26] ，再对插值后的坐标序列进行处理，但插值仅适合于小范围的缺失，对于缺失较多的情形，插值不当可能会造成伪信息 [27] 。本文根据文献[28]对含缺值GNSS基准站网共模误差分析方法和文献[29]对含缺值的数据序列的奇异谱分析方法，提出了含缺值数据序列小波分析的非插值算法，并用中国地壳运动观测网络CMONOC(Crustal Movement Observation Network of China)实测数据和模拟数据验证本文算法的有效性。

1 含缺值数据的非插值二进小波变换方法1.1 二进小波变换方法

对于长度为 N 的离散时间序列信号 x=[ x 0 x 1 ··· x N －1 ] T ，采用二进小波进行变换计算时，其样本数必须是2的整次幂，即 N =2 J , J ∈ Z ，其中 Z 表示整数集， J 为二进小波分解的最大层数。第 j 层的尺度为 τ j =2 j , ( j =0, 1, 2, ···, J －1)，相应的小波系数有 n j = N /(2 τ j )个 [30] 。若母小波函数为 φ ( t )，则经偶数平移2 k 后的第 j 层二进小波可表示为 [30]

(1)

式中, t 为基小波函数的自变量。由于前 J －1层小波系数的总数为, 比样本数 N 少一个，因此必须加上最后一层，即第 J －1层的一个尺度系数才能重构原信号。若离散信号 x的第 j 层二进小波变换表示为 [30-31]

(2)

式中, wj 为第 j 层分解的小波系数，变换矩阵 Wj 可表示为

(3)

其中, 其余 k ≠0的系数 Wj , k 都可通过 Wj , 0 平移得到，即

(4)

式中，平移算子 T j 2k 是将 Wj , 0 的各元素平移2 j +1 k 位。若将小波系数按层数从低到高排列，并加上一个尺度系数 vJ －1 ，则二进小波变换系数可表示为

(5)

式中， VJ －1 为满足 VJ －1 ⊥ Wj 的单位行向量，根据正交条件可唯一确定，尺度系数 VJ －1 = VJ －1 x。变换矩阵 W是规范正交矩阵，满足 WT W= WWT = IN ，其中， IN 表示 N 阶单位阵。

实际数据处理时，需要分解的层次 K 要远远小于 J 。如图 1所示，每一层分解将空间分解成若干低频子空间和高频子空间，将信号向由这些高低频子空间构成的正交子空间投影展开，得到信号在不同分辨率下的信息 [32] 。若将0~ f 定义为空间 V 0 ，经过一层分解后，分成0~ f /2的低频子空间CA1( v 1 )和 f /2~ f 的高频子空间CD1( W 1 )，继续对低频子空间进行分解，最后共得到CAj, CDj, …, CD1( v j , W j , …, W 1 )，其中 j 为分解层数， W j 与 v j 为正交空间且各 W 子空间也相互正交, 因此经多层分解后，原空间分解为若干不同频率的正交子空间。文献[33]在多分辨率分析思想的基础上提出了塔形算法，可用于小波变换的快速计算。图 1为采样频率为 f 的信号多分辨分析，由于GNSS坐标序列的采样频率为1d，由图可知，重构的第 j 层信号的周期为2 j 。例如，半周年信号约为182 d，周年信号约为365 d，而182∈(2 7 , 2 8 )，365∈(2 8 , 2 9 )，因此重构的第7层和第8层分量即为半周年信号和周年信号 [22] ，故将分解层数 K 定为8。

图 1 小波多分辨分析 Fig. 1 Wavelet multi-resolution analysis

图选项

当时间序列数据存在缺值时，不能用式(5)计算 w，传统方法需要先插值后再进行变换。为了避免插值计算，需要对传统的二进小波变换方法进行改进。为推导方便，将式(5)中 W的第 j 列记为 cj -1 ，则有

(6)

式中， S 和 S分别表示[0, N -1]整数集内有效观测数据集和缺失数据集。根据式(5)，缺值数据可用小波系数表示为 x j = cj T w，代入式(6)右端第2项得

(7)

记，则式(7)可表示为

(8)

考虑到 W的正交性, 有

(9)

因此

(10)

即 A是一个对称幂等阵, 由幂等阵的性质可得

(11)

式中， N S 为数据长度，当含有缺值时 N S < N ，系数矩阵 A秩亏，其秩亏数为 N － N S ，因此式(8)存在无穷多组解。为求得唯一解，引入最小范数准则

(12)

则式(8)的最小范数解为

(13)

式中， A+ 为 A的彭罗斯-穆尔逆。显然， x 没有缺值时， A= In , L= Wx, 式(13)可退化为传统的二进小波变换式(5)。

根据式(5)或式(13)求得的小波系数 w可重构离散信号 x，即

(14)

记 dj +1 = Wj T wj , sJ = VJ －1 T vJ －1 ，其中 dj +1 表示第 j +1层小波的细节信号，反映了该层尺度上的信号变化， sJ 各元素均为 x的均值。因此，式(14)将离散信号 x表示成信号的均值与各尺度细节信号之和，即

(15)

若直接对 x中的缺值补零处理，由式(6)得

(16)

式中， wz 表示由补零法获得的小波系数。将式(16)减去式(13)得

(17)

考虑到 A为对称幂等阵且tr( A)= N S ，故 A有 N S 个数值为1的特征值和N－ N S 个数值为0的特征值。将 A进行奇异值分解得

(18)

式中, U∈ RN×N 为正交阵，则 A+ 为

(19)

即对称幂等阵的逆阵等于其本身。则式(17)可改写为

(20)

顾及式(9)得

(21)

即本文方法与补零处理方法等价，同时也证明了补零处理方法本质上是一种最小范数解法。

2 算例分析2.1 模拟试验

为验证算法的有效性，按式(22) [14] 模拟一段序列

(22)

式中， y 0 、 v 分别为初始位置和速度；Δ t i = t i －t 0 为相对历元， t i 、 t 0 分别为第 i 个观测历元和初始历元； s 和 ε 分别为时变季节项和噪声。采用式(23)模拟时变季节项

(23)

式中， a 、 b 和 d 、 e 分别为周年、半周年常振幅因子， c ( t i )= fe g sin( t i ) 为变振幅因子。采用白噪声+闪烁噪声的组合策略, 由Matlab函数mvnrnd按照不同噪声的协方差阵生成随机噪声。模拟所用的参数值见表 1。

表 1 模拟所用的参数值( σ w 和 σ f 分别为白噪声和闪烁噪声振幅)Tab. 1 Parameters adopted for simulated position time series( σ w and σ f denote amplitudes of white and flicker noise)

时间

线性趋势项

季节项

噪声

y 0 /mm

v /( mm/a )

a /mm

b /mm

d /mm

e /mm

f /mm

g /mm

σ w /mm

σ f /mm

2000—2005.5

4

0.5

1.1

1.0

0.7

0.5

0.1

0.9

3.16

1.02

表选项

图 2为模拟的坐标序列及其周期功率谱，其半年和年周期信号非常明显。采用coif5小波对扣除趋势项后的坐标序列进行8层分解，重构的各层分量如图 3所示。由图可见，越靠近顶层的分量频率越低，它更多的是反映原序列中的信号成分，反之越靠近底层的分量频率越高，代表着噪声。可通过相关系数法 [34] 确定信号与噪声的分界层，该方法通过计算原序列与各层分量的相关系数，从第1层起算，相关系数第1次出现局部极小值的层数即为信号与噪声的分界层 [34] 。相关系数的计算式如下

(24)

式中， R ( x, d i )表示原序列 x ( t )与第 i 层重构分量 d i ( t )的相关系数； N 为数据长度:和分别为 x 和 d i 的均值。各层分量与原坐标序列间的相关系数见表 2。由图 3和表 2可知， d 6 为信号与噪声的分界层，因此提取 d 7 + d 8 + a 8 层分量为无数据缺失时的真实信号 S ( t )。然后，对完整的坐标序列随机删除部分数据，其中2000—2005年，每年删除数据分别为1%至6%，然后采用线性插值算法和本文算法提取信号，记为 Ŝ I ( t )和 Ŝ N ( t )，利用式(25)计算提取的信号均方根误差

(25)

式中，RMSE下标 I 、 N 分别表示插值算法和本文算法； C 为有效观测数据集； M 表示有效观测历元数。显然求得的RMSE越小，提取的信号与真实信号就越接近。采用本文算法和插值算法提取的信号与真实信号对比如图 4所示，其中本文算法和插值算法的信号均方误差分别为0.080 7 mm和0.081 3 mm，两者并无明显差异。模拟500次试验，每次试验均按上述策略随机删除数据，分别采用本文算法、插值算法以及补零算法提取信号并按式(25)计算其均方根误差。图 5(a)为本文算法和插值算法提取信号的RMSE，其平均值分别为0.086 8 mm和0.081 8 mm，二者的精度相当。图 5(b)为本文算法与补零算法提取信号的RMSE之差值，显然二者的差值非常小，其平均值为-1.09×10 －8 mm, 进一步说明了本文算法与补零算法是等价的。

图 2 模拟的坐标序列及其周期功率谱 Fig. 2 Simulated position time series and its periodic power spectrum

图选项

图 3 模拟坐标序列的各层重构分量 Fig. 3 Reconstructed component of each layer of simulated position time series

图选项

表 2 原序列与各重构分量的相关系数Tab. 2 Correlation coefficients of original signal and reconstructed components

d

R

d

R

d 1

0.631 9

d 5

0.166 51

d 2

0.443 4

d 7

0.216 9

d 3

0.329 6

d 6

0.159 1

d 4

0.252 6

d 8

0.335 1

表选项

图 4 两种算法提取的信号与真实信号对 Fig. 4 Comparison of true signals with signals extracted by two algorithms

图选项

图 5 本文算法与插值算法提取信号的RMSE及其与补零算法提取信号RMSE之差 Fig. 5 RMSE of signals extracted by the proposed algorithm and interpolation algorithm and difference between signal extracted by new algorithm and zero-padding algorithms

图选项

利用中国地壳运动观测网络GNSS数据共享服务中心(http://www.cgps.ac.cn/)提供的中国大陆27个GNSS基准站1999—2019年坐标观测序列，所有基准站的坐标序列数据都是采用GAMIT/GLOBK10.4解算得到的，详细的处理方法与步骤可以参考ftp://ftp.cgps.ac.cn/doc/processing_manual.pdf。因仪器更换等原因，GNSS站坐标序列不可避免地存在缺失。考虑到绝大部分测站在1999—2002年段缺失较多 [28, 35] , 严重破坏信号的周期成分，因此本文选取2003—2019年的数据序列作为分析。图 6为27个基准站的位置分布及2003—2019年的数据缺值率，其中基准站TASH、YONG和SHAO的缺值率分别高达45.28%、35.73%和22.49%，基准站DXIN和LHAS的数据缺值率分别为11.71%和10.27%。除这5个基准站外，其余基准站的数据缺值率均低于10%。

图 6 27个基准站的位置分布及数据缺失率 Fig. 6 Geographic locations and percentage of data missing of 27 stations in CMONOC

图选项

以BJFS站观测数据为例，扣除观测数据的趋势项后，采用coif5小波通过式(13)计算小波系数，结合小波变换矩阵重构各层分量，其结果如图 7所示。将重构的各层分量相加即可完全重构原序列。图 8为本文算法重构的序列与对缺失数据进行线性插值后的序列的对比图，其中本文算法重构的原序列缺失数据全部为零，验证了上文所证明的等价性。表 3为原序列与各层重构分量的相关系数。由表可知，3个方向坐标序列中信号与噪声的分界层分别为 d 5 、 d 5 、 d 6 ，由此可提取原序列中的信号，记为 Ŝ N ( t )。对进行线性插值后的序列也进行同样的分解与重构并提取出信号，记为 Ŝ I ( t )。图 9为本文算法和插值算法提取的信号对比，显然两种算法提取的信号均能很好地描述原序列的非线性变化。分别从原序列中扣除 Ŝ N ( t )和 Ŝ I ( t )，求得残差 ê N ( t )和 ê I ( t )，并用式(26)计算其中误差

(26)

式中， S 为有效观测数据集， M 是有效观测数据的历元数。利用本文算法和插值算法对27个基准站3个方向坐标序列进行信号提取后，利用式(26)计算和，如图 10所示。除缺值率较大的基准站(如TASH和YONG)外，两种算法提取信号后的残差中误差无明显差异。表 4列出了27个基准站坐标序列的平均残差中误差,, 二者仅相差(由式(27)计算)2.01%(North)、0.54%(East)和1.26%(Up)

(27)

表 3 原序列与各重构分量的相关系数Tab. 3 Correlation coefficients between original time series and the reconstructed components of each layer

d i

R ( x, d i )

North

East

Up

d 1

0.301 2

0.242 2

0.358 8

d 2

0.230 9

0.178 0

0..290 8

d 3

0.185 7

0.151 6

0.255 3

d 4

0.153 1

0.117 4

0.208 1

d 5

0.149 2

0.111 6

0.192 9

d 6

0.207 6

0.127 2

0.191 8

d 7

0.281 6

0.223 5

0.242 2

d 8

0.373 5

0.272 3

0.587 9

表选项

图 7 BJFS站坐标序列的各层重构分量 Fig. 7 Reconstructed components of each layer of the position time series of BJFS station

图选项

图 8 本文算法恢复的值与线性插值对比 Fig. 8 Comparison of miss values recovered by the proposed algorithm with linear interpolation

图选项

图 9 两种算法提取的信号对比 Fig. 9 Comparison of signals extracted by two algorithms

图选项

图 10 27个基准站坐标序列残差的中误差 Fig. 10 Formal errors derived from the residual time series of 27 base stations

图选项

表 4 27个基准站坐标序列的平均残差中误差Tab. 4 Mean formal error of 27 base stations

项目

North

East

Up

1.247 1

1.411 9

4.382 6

1.272 7

1.419 5

4.438 3

Imp/(%)

2.01

0.54

1.26

表选项

记 E ( t )= Ŝ I ( t )－ Ŝ N ( t )为插值算法和本文算法提取信号的差异，则 Ŝ I ( t )、 Ŝ N ( t )和 E ( t )的方差为

(28)

E ( t )与 Ŝ I ( t )、 Ŝ N ( t )的方差比为

(29)

利用本文算法和插值算法提取27个基准站坐标序列的信号后，按式(29)计算 p I 和 p N , 如图 11所示。表 5为27个基准站坐标序列的平均方差比 pI 、 pN 。由图 11和表 5可知， pI 、 pN 非常接近，其差值，仅为1.16%(North)、0.54%(East)和1.62%(Up)。

图 11 27个基准站坐标序列的 p I 和 p N Fig. 11 p I and p N derived from signals of 27 base stations

图选项

表 5 27个基准站坐标序列的平均方差比Tab. 5 pI and pN of 27 base stations

项目

North

East

Up

pI /(%)

2.99

2.37

3.33

pN /(%)

4.15

2.91

4.95

Δ p/(%)

1.16

0.54

1.62

表选项

3 结语

GNSS基准站坐标序列中存在季节性信号，采用传统的二进小波变换提取信号时要求完整的坐标序列数据；当数据存在缺值时，再用传统二进小波变换公式计算，对缺失数据进行插值或补零处理。本文根据小波系数与坐标序列数据的重构关系，提出一种非插值的二进小波变换算法，并基于最小范数准则导出了计算公式。当数据没有缺值时，本文算法退化为传统算法；数据有缺失时，本文方法与传统补零算法等价，因此传统补零算法本质上是最小范数解法。最后通过模拟数据和实测数据验证了本文算法，且本文算法与插值算法提取的信号相差很小。

精

彩

回

顾