非线性独立分量分析（论文推荐嵇昆浦）

非线性独立分量分析（论文推荐嵇昆浦）收稿日期：2019-05-08；修回日期：2020-01-20嵇昆浦 沈云中同济大学测绘与地理信息学院 上海 200092【测绘学报的个人主页】长按复制此条消息，长按复制打开抖音查看TA的更多作品##7NsBSynuc88##[抖音口令]【引文格式】嵇昆浦 沈云中. 含缺值GNSS基准站坐标序列的非插值小波分析与信号提取. 测绘学报，2020，49(5)：537-546. DOI: 10.11947/j.AGCS.2020.20190163含缺值GNSS基准站坐标序列的非插值小波分析与信号提取

《测绘学报》

构建与学术的桥梁 拉近与权威的距离

《测绘学报》抖音自开通以来，聚焦于测绘地理信息学术前沿进展，受到了广大专家学者的大力支持，播放量数万，粉丝1.7万。

复制链接，关注我们哦！

【测绘学报的个人主页】长按复制此条消息，长按复制打开抖音查看TA的更多作品##7NsBSynuc88##[抖音口令]

【引文格式】嵇昆浦 沈云中. 含缺值GNSS基准站坐标序列的非插值小波分析与信号提取. 测绘学报，2020，49(5)：537-546. DOI: 10.11947/j.AGCS.2020.20190163

含缺值GNSS基准站坐标序列的非插值小波分析与信号提取

嵇昆浦 沈云中

同济大学测绘与地理信息学院 上海 200092

收稿日期：2019-05-08；修回日期：2020-01-20

基金项目：国家自然科学基金（41731069）

第一作者简介：嵇昆浦(1995-) 男 硕士生 研究方向为GNSS数据处理理论与方法。E-mail:1575540259@qq.com

通信作者：沈云中，E-mail:yzshen@tongji.edu.cn

摘要：受多种因素影响，GNSS基准站坐标序列通常都含有缺值，传统小波分析需要对缺值数据进行内插或补零处理。本文基于小波系数与时间序列观测数据的重构关系，提出了一种非插值的二进小波变换的最小范数解法，导出了相应的计算式，并严格证明了传统的补零处理算法与本文的最小范数解法等价。最后利用中国地壳运动观测网络一期27个基准站实测数据以及模拟数据进行了验证分析。结果表明，本文的非插值算法与插值算法提取的信号差异较小，27个基准站坐标序列的平均残差中误差仅相差2.01%（North），0.54%（East）和1.26%（Up），两种算法提取的信号之差与信号平均方差比仅相差1.16%（North），0.54%（East）和1.62%（Up）。

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

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

JI Kunpu SHEN Yunzhong

College of Surveying and Geoinformatics Tongji University Shanghai 200092 China

Foundation support: The National Natural Science Foundation of China (No.41731069)

First author: JI Kunpu (1995-) male postgraduate majorsinthe theories and methods of GNSS data proG cessing.E-mail:1575540259@qq.com.

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₀x₁···xN－1]^T，采用二进小波进行变换计算时，其样本数必须是2的整次幂，即N=2J J∈Z，其中Z表示整数集，J为二进小波分解的最大层数。第j层的尺度为τ_j=2j (j=0 1 2 ··· J－1)，相应的小波系数有n_j=N/(2τ_j)个^[30]。若母小波函数为φ(t)，则经偶数平移2k后的第j层二进小波可表示为^[30]

(1)

式中 t为基小波函数的自变量。由于前J－1层小波系数的总数为

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

(2)

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

(3)

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

(4)

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

(5)

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

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

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

图选项

1.2 含缺值数据的改进二进小波变换

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

(6)

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

(7)

记

(8)

考虑到W的正交性 有

(9)

因此

(10)

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

(11)

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

(12)

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

(13)

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

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

(14)

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

(15)

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

(16)

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

(17)

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

(18)

式中 U∈R^N×N为正交阵，

则A 为

(19)

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

(20)

顾及式(9)得

(21)

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

2 算例分析2.1 模拟试验

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

(22)

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

(23)

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

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

时间	线性趋势项		季节项						噪声
	y0/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₆为信号与噪声的分界层，因此提取d₇ d₈ a₈层分量为无数据缺失时的真实信号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^－8mm 进一步说明了本文算法与补零算法是等价的。

图 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
d1	0.631 9	d5	0.166 51
d2	0.443 4	d7	0.216 9
d3	0.329 6	d6	0.159 1
d4	0.252 6	d8	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

图选项

2.2 实测试验

利用中国地壳运动观测网络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₅、d₅、d₆，由此可提取原序列中的信号，记为Ŝ_N(t)。对进行线性插值后的序列也进行同样的分解与重构并提取出信号，记为Ŝ_I(t)。图 9为本文算法和插值算法提取的信号对比，显然两种算法提取的信号均能很好地描述原序列的非线性变化。分别从原序列中扣除Ŝ_N(t)和Ŝ_I(t)，求得残差ê_N(t)和ê_I(t)，并用式(26)计算其中误差

(26)

式中，S为有效观测数据集，M是有效观测数据的历元数。利用本文算法和插值算法对27个基准站3个方向坐标序列进行信号提取后，利用式(26)计算

(27)

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

di	R(x di)
	North	East	Up
d1	0.301 2	0.242 2	0.358 8
d2	0.230 9	0.178 0	0..290 8
d3	0.185 7	0.151 6	0.255 3
d4	0.153 1	0.117 4	0.208 1
d5	0.149 2	0.111 6	0.192 9
d6	0.207 6	0.127 2	0.191 8
d7	0.281 6	0.223 5	0.242 2
d8	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个基准站坐标序列的平均方差比p_I、p_N。由图 11和表 5可知，p_I、p_N非常接近，其差值，仅为1.16%(North)、0.54%(East)和1.62%(Up)。

图 11 27个基准站坐标序列的p_I和p_NFig. 11p_Iandp_Nderived from signals of 27 base stations

图选项

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

人物 | 测绘学报主编杨元喜院士：珠峰高程测量最高精度从何而来

招聘启事｜武汉大学测绘学院诚聘海外优秀人才（2020年）

通知 | 自然资源部地理国情监测重点实验室2020年开放基金申请指南

长图｜用生命测量珠峰的中国故事

测绘图书 | 《新型地理计算模式及其在双评价中的应用》实现版权输出

孙文舟 殷晓冬 暴景阳 I《测绘学报（英文版）》（JGGS）精选论文

权威 | 专业 | 学术 | 前沿微信、抖音小视频投稿邮箱 | song_qi_fan@163.com

微信公众号中搜索「测绘学报」，关注我们，长按上图二维码，关注学术前沿动态。

欢迎加入《测绘学报》作者QQ群： 751717395

进群请备注：姓名 单位 稿件编号