Denoising Method for Deformation Monitoring Data Based on ICEEMD-ICA and MDP Principle
-
摘要: 针对经验模态分解(empirical mode decomposition,EMD)方法存在信噪分离不准确的缺陷,以及独立分量分析(independent component analysis,ICA)存在不确定性的问题,提出了一种改进完备集成经验模态分解(improved complete ensemble empirical mode decomposition, ICEEMD)、ICA与最小失真准则(minimal distortion principle,MDP)相结合进行变形数据去噪的方法。首先,使用ICEEMD方法对变形监测数据进行有效分解,并以此构建虚拟噪声信号;其次,对虚拟噪声进行二次ICEEMD分解,提取更接近真实噪声的二次虚拟噪声信号,再以二次虚拟噪声和原变形数据组成输入观测通道,使用ICA进行处理;然后,通过计算ICA处理后的独立分量与输入信号的相关系数,解决独立分量的排序不确定性与相位不确定性问题;最后,使用MDP准则有效解决了独立分量的幅值不确定性。对加噪仿真数据和实际桥梁GNSS变形监测数据进行详细分析,结果表明,所提方法可取得良好的去噪效果,有效提升去噪的性能指标,充分验证了所提方法在变形监测数据去噪中具备的可行性和有效性。
-
关键词:
- 改进完备集成经验模态分解 /
- 独立分量分析 /
- 二次虚拟噪声 /
- 最小失真准则 /
- 变形监测数据去噪
Abstract:Objectives Considering the inaccurate separation of signal and noise of empirical mode decomposition (EMD) method and the uncertainty of independent component analysis (ICA), a new method for denoising deformation data with improved complete ensemble empirical mode decomposition (ICEEMD), independent component analysis (ICA) and minimal distortion principle (MDP) is proposed.Methods Firstly, ICEEMD method is used to decompose the deformation monitoring data effectively, and the virtual noise signal is constructed. Secondly, ICEEMD decomposition of virtual noise is carried out to extract twice virtual noise signal which is closer to real noise. The input observation channel is composed of twice virtual noise and original deformation data and processed by ICA. Then, by calculating the correlation coefficient between the independent components and the input signal after ICA processing, the sorting uncertainty and phase uncertainty of independent components can be solved. Finally, the MDP criterion is used to effectively solve the amplitude uncertainty of independent components.Results Through the detailed analysis of noisy simulation data and actual bridge GNSS deformation monitoring data, the results show that the proposed method has achieved good denoising effect and can effectively improve the performance of denoising.Conclusions It also fully verified the feasibility and effectiveness of the proposed method indenoising of deformation monitoring data. -
随着北斗、GPS、GLONASS、伽利略等卫星导航系统星座的逐步完善,目前可用卫星数已达78颗,GNSS(global navigation satellite system)除了应用于传统精密测绘领域,还在手机定位服务、移动测量技术、机械控制、产量监测、智能交通系统和无人机等新兴领域得到了越来越广泛的应用[1]。与高精度高可靠性定位需求的传统精密测绘相比,这些新兴应用领域对定位精度要求相对较低,但对接收机的成本、体积及在复杂环境下的可用性等因素较为敏感。传统的精密测绘一般采用大地测量型接收机,定位精度高,但价格较为昂贵、体积相对较大。用户由于使用成本的限制更多地会采用单频多模接收机,为满足大众化的应用需求,低成本的GNSS单频精密定位方法具有重要的研究价值。而单频周跳探测与修复方法又是GNSS单频精密定位亟待解决的关键问题。
与传统多站相对定位数据处理模式相比,PPP(precise point positioning)、PPP-AR(PPP ambiguity resolution)、URTK(undifferenced network RTK)等单站处理模式虽然在理论基础上与多站模式等价,但由于其在数据处理时的灵活性,可以满足大规模网、大区域的分层并行解算,因此越来越被人们所接受[2-8]。与多站处理模式相比,单站处理模式下无法利用周边测站的观测数据消除观测值中电离层、对流层、卫星钟差等误差影响。此外,单频数据预处理无法采用MW(melbourne-wubbena)、GF(geometry free)、IF(ionosphere free)组合观测值辅助周跳探测,因此单站单频数据预处理较多站模式以及单站双频/多频模式要更为困难。
目前,针对单频数据预处理主要有高次差法、多项式拟合法、Kalman滤波法、小波变换法、多普勒观测值辅助探测法等。多项式拟合法探测周跳算法的基础是假设观测值随时间的变化可以用一个高阶多项式来表示,这一假设很容易被接收机自身的运动所打破,因此该方法通常不适用于动态定位中周跳的探测[9]。3阶多项式模型的Kalman滤波探测周跳,仅适用于高采样率的观测数据[10]。基于双差的单频单站数据预处理方法,可以有效探测出大部分的周跳观测值,但如果参考卫星的观测值发生周跳,组成双差观测方程时,会把参考卫星的周跳值引入到每一个双差观测值中,这就为确定发生周跳的卫星号带来一定的困难, 甚至导致误判[11]。综上所述,针对单站单频数据预处理,现有方法还存在一定的局限性,为此本文提出一种历元间单差观测值的单站单频周跳探测与修复方法。
1 周跳探测原理
1.1 探测周跳发生历元
假定测站i在历元t的位置为基准站p,在历元t+1的位置为流动站q, 组成历元间单差观测方程式(1)。相邻历元间电离层延迟和对流层延迟的变化量有限,卫星钟变化微小,历元间单差观测值可基本消除电离层、对流层及卫星钟差的影响,故式(1) 不考虑这些误差。残余的电离层、对流层及卫星钟差被吸收至Δε中。
$$ \lambda \Delta \varphi = \rho {s_q}-\rho {s_p}-c\Delta t-\lambda \Delta N + \Delta \varepsilon $$ (1) 式中,Δφ为历元间单差观测值;ρsp为历元t的卫地距;ρsq为历元t+1的卫地距;Δt为接收机钟变化量;ΔN为周跳值;Δε为观测噪声;c为光速;λ为载波波长。
即使相邻历元的观测值都是同一台接收机观测值,历元间单差观测值也消除不了接收机钟差影响,需作为参数估计。假定观测值未发生周跳,则前后历元的整周模糊度相等可通过历元间单差消除,因此共存在3个坐标参数和1个接收机钟差参数,如需确定发生周跳的观测值或进行周跳修复,需保证至少存在4颗卫星的观测数据未发生周跳。采用相对定位模式可解出测站i在历元t+1处的观测值验后单位权中误差。观测值未发生周跳,如图 1所示,验后单位权中误差一般在毫米级。如果发生周跳即ΔN不为零,图 2(b)和图 3(b)中验后单位权中误差会显著增大,因此可根据验后单位权中误差探测发生周跳的历元。
1.2 确定周跳发生卫星
在测量数据服从正态分布的情形下,最小二乘估计具有最优统计性质,即它是最优线性无偏估计[12]。但是由于最小二乘具有良好均衡误差的特性,观测值中有粗差出现时,估值必将迁就粗差,因此个别的粗差会对估值产生较大影响,导致估计严重失实,因此本文采用抗差最小二乘法确定发生周跳的异常卫星,等价权函数的选取如式(2)[13]:
$$ {\bar P_i} = \left\{ \begin{array}{l} {P_i}{, _i}|{{\hat V}_i}| \le {k_0}\\ \frac{{{k_0}{P_i}}}{{|{{\hat V}_i}|}}{\left( {\frac{{{k_1}-|{{\hat V}_i}|}}{{{k_1}-{k_0}}}} \right)^2}, {k_0} < |{{\hat V}_i}| \le {k_1}\\ 0, |{{\hat V}_i}| > {k_1} \end{array} \right. $$ (2) 式中,$ {{\hat V}_i} $为标准化残差; Pi采用高度角定权;本文k0取1.0, k1取2.5。
采用抗差最小二乘法解算法方程,逐步进行迭代计算,直至前后两次解的坐标差值符合限差要求为止。此时可获得每个观测值的残差vi,当标准化残差$|{{\hat V}_i}| = \frac{{|{v_i}|}}{{\sqrt {{Q_{vv}}} }} $大于某一阈值(本文取0.05 m)时,可确定该观测值有周跳发生。Qvv只与图形结构及观测值的权阵有关,权阵采用先验权即高度角定权。
1.3 周跳修复
当确定某颗卫星观测值发生周跳时,在式(1) 中增加一个周跳值待估参数ΔN。循环发生周跳的异常卫星,每次循环确定一颗卫星的周跳值,其余发生周跳的异常卫星观测值权降为0,使其观测值对参数估计不起作用。直至修复所有周跳值,结束循环。
2 有效性测试分析
采用2010年5月21日山西省连续运行参考站系统(SXCORS)一组1 s采样间隔的实测静态和动态数据验证本文方法的有效性。实测数据为双频观测数据,为模拟单频实时周跳探测与修复的过程,仅采用L1、CA观测值及卫星广播星历。
2.1 实测静态数据测试
为验证本文方法对实测静态观测数据的有效性,选取SXCORS中某一基准站的观测数据进行试验,GPS接收机和天线类型分别为LEICA GRX1200GGPRO和LEIAR25 LEIT,能够采集C1、L1、P2和L2共4种观测值。本文测试了前一小时共3 600历元连续无周跳的观测数据,平均可视卫星数为9.80。按照本文方法逐历元计算出历元间单差观测值的单位权中误差(见图 1)。各历元的单位权中误差均为毫米级,验证了所采用的实验数据并不存在周跳。
文献[11]中验证了基于双差的单频单站数据预处理方法相对多项式拟合的优越性,是单频数据预处理中较为有效的一种方法,以下进一步验证本文方法相对该方法的有效性。在历元100处对参考卫星G15添加1周周跳,图 2显示了两种方案的处理结果,其中A方案为基于双差的单频单站数据预处理方法,B方案为本文方法。图 2(a)显示了观测值未发生周跳时的抗差估计验后残差,两种方案得到的验后残差基本一致。参考卫星发生周跳时,图 2(b)中两种方案的验后单位权都被显著的放大,A方案和B方案都可以探测出周跳发生的历元。然而如图 2(c)所示,A方案中参考卫星的周跳值被星间单差观测值吸收,每颗卫星的验后残差被不正常地放大从而导致周跳探测失效。而本文方法由于无需设置参考卫星,观测值间不存在相关性,如图 2(d)所示,只有卫星G15的验后残差被显著的放大。因此,任意卫星发生周跳均不会影响本文周跳探测方法的有效性。
本文方法能否正确定位周跳的位置取决于发生周跳的个数。在复杂观测环境下,发生周跳的异常卫星数、周跳值大小存在多样性,周跳探测与修复的有效性可能会有所降低。为系统测试分析本文周跳探测与修复方法在不同观测条件下的抗差性,每个历元分别添加n个不同的周跳,周跳值大小分别为1,2,…,n,参与统计的历元在观测时段内都至少存在4颗卫星的观测数据未发生周跳,历元数总计3 600,统计结果见表 1。
表 1 静态数据周跳探测与修复成功率统计Table 1. Statistics on the Rate of the Successful Cycle Slip Detection and Repair of Static Data(发生周跳卫星数/可视卫星数)/% 周跳探测成功率/% 探测到的周跳修复成功率/% 0~10 98.88 100 10~15 99.25 100 15~20 99.11 100 20~25 96.83 100 25~30 95.86 100 30~35 70.50 100 35~40 58.11 100 40~45 34.31 100 45~50 0.94 100 从表 1可见,当发生周跳的异常卫星数小于可视卫星数的30%时,周跳探测成功率在95%以上;当发生周跳的异常卫星数大于可视卫星数的30%时,随着发生周跳异常卫星数的增加,周跳探测的成功率逐渐下降。由表 1还可以看出,如果周跳探测成功,则周跳修复成功率高达100%。
2.2 实测动态数据测试
为了进一步验证本文周跳探测与修复方法的有效性和可靠性,利用SXCORS网中的一组实测动态数据进行验证分析,GPS接收机和天线类型分别为Trimble 5800和TRM5800,能够采集C1、L1、L2和P2共4种类型的观测类型。选取时段10:50:40~10:55:21共280个历元连续无周跳的观测数据,平均可视卫星数为9,利用本文方法逐历元计算出历元间单差观测值的单位权中误差(图 3(a)),实测动态数据所得的中误差相比实测静态数据略高,但基本在同一量级,验证了所采用的实验数据并不存在周跳。
对G15卫星在历元100处添加1周周跳,图 3(b)显示了观测值添加周跳后解算得到的验后单位权中误差。从图 3(b)中可以看出,在历元100处的验后单位权中误差显著增大,可探测出该历元有周跳发生,采用与静态测试相同的处理流程便可实现周跳探测与修复。为测试本文方法对动态数据的抗差性,选取10:50:41~11:09:01时段中无周跳发生的历元为统计时段,平均可视卫星数为9。每个历元分别添加n个不同的周跳,周跳值大小分别为1,2,…,n,参与统计的历元在观测时段内都至少存在4颗卫星的观测数据未发生周跳,历元数总计1 035,统计结果见表 2。
表 2 动态数据周跳探测与修复成功率统计Table 2. Statistics on the Rate of the Successful Cycle Slip Detection and Repair of Dynamic Data(发生周跳卫星数/可视卫星数)/% 周跳探测成功率/% 探测到的周跳修复成功率/% 0~15 98.74 100 15~20 98.84 100 20~25 99.61 100 25~30 99.61 100 30~35 42.03 100 35~40 34.49 100 40~45 17.29 100 45~50 0.48 100 由表 2可以看出,即使采用实测动态数据,仍可得到与静态数据类似的测试分析结果,这进一步验证了本文方法在进行单站单频周跳探测与修复时的有效性。
3 结语
基于实测静态和动态数据的测试结果表明,按照本文方法,根据历元间差分观测值的验后单位权中误差,能够100%探测出发生周跳的历元。当至少4颗卫星未发生周跳时,如异常卫星数小于可视卫星数的30%,则在95%以上的情况下可以有效确定异常卫星。当异常卫星过多时,本文方法确定异常卫星的成功率会有所下降。当无法确定异常卫星时,建议将所有可视卫星重新初始化。整体而言,只要能够探测出发生周跳的异常卫星,本文方法对其观测值进行修复的成功率高达100%。本文提出的单站单频周跳探测与修复方法,对于GNSS单频精密定位技术发展,以及将相对廉价、轻便的GNSS单频接收机设备应用于城市定位服务、精准农业、智能交通系统等新兴领域具有一定的推动作用。
-
表 1 4种方法的去噪指标
Table 1 Denoising Indexes of Four Methods
去噪方法 SNR/dB RMSE ICEEMD 10.727 3 0.209 8 EMD-ICA 9.514 7 0.241 3 ICEEMD-ICA 10.152 1 0.224 2 本文方法 12.765 7 0.165 9 表 2 实例数据中4种方法的去噪指标
Table 2 Denoising Indexes of Four Methods in Instance Data
去噪方法 SNR/dB RMSE/ m ICEEMD 24.964 2 7.768 8 EMD-ICA 23.263 4 9.449 1 ICEEMD-ICA 25.518 0 7.288 9 本文方法 29.129 8 4.809 2 -
[1] 章浙涛, 朱建军, 匡翠林, 等. 一种小波包混合滤波方法及其应用[J]. 武汉大学学报·信息科学版, 2014, 39(4): 471-475 doi: 10.13203/j.whugis20120141 Zhang Zhetao, Zhu Jianjun, Kuang Cuilin, et al. A Hybrid Filter Method Based on Wavelet Packet and Its Application[J]. Geomatics and Information Science of Wuhan University, 2014, 39(4): 471- 475 doi: 10.13203/j.whugis20120141
[2] 栾元重, 栾亨宣, 李伟, 等. 桥梁变形监测数据小波去噪与Kalman滤波研究[J]. 大地测量与地球动力学, 2015, 35(6): 1 041-1 045 https://www.cnki.com.cn/Article/CJFDTOTAL-DKXB201506031.htm Luan Yuanzhong, Luan Hengxuan, Li Wei, et al. Research on Wavelet Denoising and Kalman Filter in Bridge Deformation Data[J]. Journal of Geodesy and Geodynamics, 2015, 35(6): 1 041-1 045 https://www.cnki.com.cn/Article/CJFDTOTAL-DKXB201506031.htm
[3] Yi T H, Li H N, Gu M. Wavelet Based Multi-step Filtering Method for Bridge Health Monitoring Using GPS and Accelerometer[J]. Smart Structures and Systems, 2013, 11(4): 331-348 doi: 10.12989/sss.2013.11.4.331
[4] Huang N E, Shen Z, Long S R, et al. The Empirical Mode Decomposition and the Hilbert Spectrum for Nonlinear and Non-stationary Time Series Analysis[J]. Proceeding of the Royal Society: A Mathematical Physical & Engineering Sciences, 1998, 454 (1): 903-995 http://www.ccpo.odu.edu/~klinck/Reprints/PDF/huangPRSLA1998.pdf
[5] 罗亦泳, 黄城, 张静影. 基于变分模态分解的变形监测数据去噪方法[J]. 武汉大学学报·信息科学版, 2020, 45(5): 784-790 doi: 10.13203/j.whugis20180437 Luo Yiyong, Huang Cheng, Zhang Jingying. Denoising Method of Deformation Monitoring Data Based on Variational Mode Decomposition[J]. Geomatics and Information Science of Wuhan University, 2020, 45(5): 784-790 doi: 10.13203/j.whugis20180437
[6] 范千. 单历元GPS变形信息特征提取的EMD方法[J]. 山东科技大学学报(自然科学版), 2011, 30(5): 78-82 https://www.cnki.com.cn/Article/CJFDTOTAL-SDKY201105016.htm Fan Qian. Empirical Mode Decomposition Method for Single Epoch GPS Deformation Information Characteristic Extracting[J]. Journal of Shandong University of Science and Technology, 2011, 30 (5): 78-82 https://www.cnki.com.cn/Article/CJFDTOTAL-SDKY201105016.htm
[7] 罗飞雪, 戴吾蛟. 小波分解与EMD在变形监测应用中的比较[J]. 大地测量与地球动力学, 2010, 30(3): 137-141 https://www.cnki.com.cn/Article/CJFDTOTAL-DKXB201003034.htm Luo Feixue, Dai Wujiao. Comparison of EMD with Wavelet Decomposition for Dynamic Deformation Monitoring[J]. Journal of Geodesy and Geodynamics, 2010, 30(3): 137-141 https://www.cnki.com.cn/Article/CJFDTOTAL-DKXB201003034.htm
[8] Wu Z H, Huang N E. Ensemble Empirical Mode Decomposition: A Noise Assisted Data Analysis Method[J]. Advances in Adaptive Data Analysis, 2009, 1 (1): 1-41 doi: 10.1142/S1793536909000047
[9] Torres M E, Colominas M A, Schlotthauer G, et al. A Complete Ensemble Empirical Mode Decom position with Adaptive Noise[C]// International Conference on Acoustics, Speech and Signal Processing(ICASSP), Prague, Czech Republic, 2011
[10] Colominas M A, Schlotthauer G, Torres M E. Improved Complete Ensemble EMD: A Suitable Tool for Biomedical Signal Processing[J]. Biomedical Signal Processing and Control, 2014, 14(1): 19-29 http://download.xuebalib.com/n2fU9SncaWt.pdf
[11] 吴金斌, 周世健. 经验模式分解联合独立分量分析降噪研究[J]. 测绘科学, 2016, 41(7): 197-201 https://www.cnki.com.cn/Article/CJFDTOTAL-CHKD201607035.htm Wu Jinbin, Zhou Shijian. The De-noising Study Based on Empirical Mode Decomposition and Independent Component Analysis[J]. Science of Surveying and Mapping, 2016, 41(7): 197-201 https://www.cnki.com.cn/Article/CJFDTOTAL-CHKD201607035.htm
[12] 于帅, 刘超, 李盟盟, 等. EMD-WaveletICA耦合模型及其在GPS坐标序列降噪中的应用[J]. 测绘科学技术学报, 2016, 33(2): 139-144 https://www.cnki.com.cn/Article/CJFDTOTAL-JFJC201602006.htm Yu Shuai, Liu Chao, Li Mengmeng, et al. EMDWavelet-ICA Coupled Model and Its Application in GPS Coordinate Series De-noising[J]. Journal of Geomatics Science and Technology, 2016, 33(2): 139-144 https://www.cnki.com.cn/Article/CJFDTOTAL-JFJC201602006.htm
[13] 赫彬, 张雅婷, 白艳萍. 基于ICA-CEEMD小波阈值的传感器信号去噪[J]. 振动与冲击, 2017, 36(4): 226-232 https://www.cnki.com.cn/Article/CJFDTOTAL-ZDCJ201704036.htm He Bin, Zhang Yating, Bai Yanping. A Method for Sensor Signal De-noising Based on ICA-CEEMD Wavelet Threshold[J]. Journal of Vibration and Shock, 2017, 36(4): 226-232 https://www.cnki.com.cn/Article/CJFDTOTAL-ZDCJ201704036.htm
[14] 王胜, 李磊. 加权优选方法在信号分解中的应用与软件研制[J]. 测绘地理信息, 2020, 45(5): 54-58 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXG202005013.htm Wang Sheng, Li Lei. Application and Software Development of Weighted Optimization Method Used on Signal Decomposition[J]. Journal of Geomatics, 2020, 45(5): 54-58 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXG202005013.htm
[15] 徐信, 李志华, 杨越. 基于EMD-ICA的激电数据降噪处理方法[J]. 计算机应用研究, 2017, 34(6): 1 737-1 741 https://www.cnki.com.cn/Article/CJFDTOTAL-JSYJ201706031.htm Xu Xin, Li Zhihua, Yang Yue. De-noising of IP Data Based on EMD-ICA[J]. Application Research of Computers, 2017, 34(6): 1 737-1 741 https://www.cnki.com.cn/Article/CJFDTOTAL-JSYJ201706031.htm
[16] 谢将剑, 张军国, 董方舟, 等. 基于虚拟观测的FastICA的坠砣位移的风荷载分量去除[J]. 铁道学报, 2014, 36(7): 44-50 https://www.cnki.com.cn/Article/CJFDTOTAL-TDXB201407008.htm Xie Jiangjian, Zhang Junguo, Dong Fangzhou, et al. FastICA Remove of Windload Component from Balance Weight Displacement Based on Virtual Observation[J]. Jouranl of the China Railway Society, 2014, 36(7): 44-50 https://www.cnki.com.cn/Article/CJFDTOTAL-TDXB201407008.htm
[17] Hyvarinen A. Fast and Robust Fixed-Point Algorithms for Independent Component Analysis[J]. IEEE Transactions on Neural Networks, 1999, 10 (3): 626-634 http://www.cs.helsinki.fi/u/ahyvarin/papers/TNN99_reprint.pdf
[18] 刘东瀛, 邓艾东, 刘振元, 等. 基于EMD与相关系数原理的故障声发射信号降噪研究[J]. 振动与冲击, 2017, 36(19): 71-77 https://www.cnki.com.cn/Article/CJFDTOTAL-ZDCJ201719011.htm Liu Dongying, Deng Aidong, Liu Zhenyuan, et al. De-noising Method for Fault Acoustic Emission Signals Based on the EMD and Correlation Coefficient[J]. Journal of Vibration and Shock, 2017, 36 (19): 71-77 https://www.cnki.com.cn/Article/CJFDTOTAL-ZDCJ201719011.htm
[19] Matsuoka K. Minimal Distortion Principle for Blind Source Separation[C]// The 41st SICE Annual Conference, New York, USA, 2002
[20] Meng X L, Nguyen D T, Xie Y L, et al. Design and Implementation of a New System for Large Bridge Monitoring-GeoSHM[J]. Sensors, 2018, 18(1): 775-797 http://www.mdpi.com/1424-8220/18/3/775/pdf