基于小波变换的高边坡变形监测数据误差处理方法

2021-01-05雷孟飞周俊华汤金毅

导航定位学报 2020年6期

雷孟飞，周俊华，汤金毅

（湖南联智科技股份有限公司，长沙 410073）

0 引言

北斗卫星导航系统（BeiDou navigation satellite system, BDS）是我国自主研发的全天时、全天候提供高精度位置信息的导航定位系统。随着BDS的逐步完善，BDS 定位精度不断提高，BDS 高精度定位技术也逐渐被应用到各种工程领域施工期和维护期的变形监测中[1]。利用BDS 进行变形监测，可以克服人工测量所存在的受自然条件限制严重、不能实现自动化监测、监测成本高等弊端，因此在变形监测中的应用越来越广泛。

BDS 自动化监测系统采用静态相对定位的算法进行高精度解算，由于基准站和监测站的距离相对较近，解算结果水平方向的精度通常能够达到±3 mm，高程方向能达到±5 mm。监测结果中的大多误差已经通过差分算法消除，但是对于一些周期性的误差，比如多路径误差、周跳以及设备噪声造成的高频误差和异常跳动则无法剔除干净，特别对于一些高边坡监测项目，通常信号遮挡严重、卫星星座分布不均，周期性误差和高频误差明显，并不能得到有效剔除[2]。这些周期性误差和高频误差为边坡监测带来很大困难，技术人员很难从各项误差中提取有用信号来判断边坡是否发生形变以及是否需要形变预警，因此有必要对BDS变形监测的数据进行降噪处理，获得监测点的真实位移坐标时间序列。

小波变换是1 种新的变换分析方法，它克服了傅里叶变换窗口大小不能随频率变化等缺点，提供随频率改变的“时间-频率”窗口，可以探测信号中的瞬态成分及频率成分[3]。在BDS 高边坡变形监测过程中，可以利用小波变换将监测结果中的高频噪声信号分离出来并进行去除，提取有效的形变信息，提高监测精度和稳定性。国内外学者已经将小波变换较多地运用在全球定位系统（global positioning system, GPS）坐标时间序列分析中。文献[4]提出用小波变换方法提取 GPS 观测数据中的多路径信号，对观测结果进行修正；文献[5]利用小波分析对全球卫星导航系统（global navigation satellite system, GNSS）高频低频多路径误差进行滤波降噪，效果良好；文献[6]提出了1 种基于小波滤波与主成分分析相结合的GPS 去噪方法，构建GPS多路径误差的误差模型，并对GPS 多路径误差进行修正；文献[7]采用小波变换，提取了中国区域8 个国际GNSS 服务组织（International GNSS Service,IGS）站的坐标时间序列中的非线性运动信息，对季节性周期变化进行了分析研究。上述研究中所采用的GPS 坐标时序数据大多为IGS 站点数据，数据质量较高，没有考虑实际应用中卫星被遮挡、其他信号干扰的情况，因此小波变换在BDS 高边坡监测中的应用效果并不明确，并且上述研究中并没有针对性地对坐标时间序列中不同频率的误差项进行单独改正，监测结果容易丢失有用的细节部分。

针对上述问题，本文以广西某高速公路高边坡为实验测区，首先利用小波变换算法对监测结果进行分解，利用低频层恢复测试数据，剔除监测数据中的粗差数据，然后对监测数据进行2 次小波变换，并对变换后的各层系数分别进行改正，最后恢复监测数据，分析误差剔除效果。

1 高边坡变形监测数据误差处理方法

1.1 小波变换

小波变换主要分为连续小波变换和离散小波变换 2 类[8]。本文采用较常用的离散小波变换（discrete wavelet transform, DWT）进行处理，离散小波变换基本公式[9]为

式中：φ ( t)为小波母函数，对其进行平移和伸缩得到小波基函数集 φa,b(t )；a 为尺度因子；b 为平移因子；R 为实数集。通常取

式中：Z 为整数集； a0为尺度因子扩展步长；0b 为平移因子扩展步长；m 为尺度因子扩展系数；n 为平移因子扩展系数。

将式（2）代入式（1）可得[10]

离散小波变换公式为

式中： fDWT( a , b )为连续小波变换系数；*表示共轭[11]。

1.2 变形监测数据误差处理流程

基于小波变换的高边坡变形监测数据误差处理方法其基本思想是：首先利用小波变换对监测结果的粗差进行剔除；然后对处理后的时间序列进行2 次小波变换；最后将变换后的系数分为高频、中频、低频3 类，分别对每层细节系数采用不同的滤波算法进行滤波，利用各层滤波后的细节系数恢复监测时间序列。其方法流程如图1 所示，其中各层系数改正的处理流程如图2 所示。

图1 数据处理流程

1）数据预处理。选择监测点一段时间的监测数据进行处理，假设监测数据采样频率为s，监测数据时间序列为Y，数据预处理采用3 次样条插值算法对缺失数据进行插值；

2）小波变换。采用小波变换对监测点数据进行分解，分解后得到细节信号（监测结果的细节信息）以及趋势信号（包含监测点的位移趋势）；

3）分解系数分类。①分别计算第2 步的分层结果的频带范围，计算公式[12]为

图2 各层系数改正流程

式中：Y ' 为监测结果的采样频率，通用的单位采样频率为1 s/次，而本文中变形监测的采样频率较低，因此，本文中定义的单位采样频率为1 h/次；j 为分解后的第j 层；式（5）为细节系数的频带范围，式（6）为最后一层分解后的近似系数的频带范围。②根据不同频段的阈值，将上一步骤计算得到的各层的频带范围进行区分，其中频段的阈值为各分层系数高频、中频和低频的区分值，高频信号定义为频率大于24 h/次，中频信号约为24 h/次即周期约为24 h，低频信号频率小于24 h/次，其具体数值根据监测数据的长度和采样频率确定。

4）粗差识别修复。①利用第3 步得到的低频层和中频系数采用小波逆变换恢复监测数据，得到监测数据 Y '；②利用3 倍中误差原则识别粗差数据，具体为：|Yi－ Yi'| ＞ 3σ的数据为粗差数据，其中 Yi为第i 个原始监测数据， Yi'为第i 个低频系数恢复后的原始数据；③利用临近插值法对识别出的粗差数据进行修复，得到粗差修复后的监测数据 Y'' 。

5）高频系数滤波。将步第4 步2 次小波变换后得到的各层变换系数按照第3 步的方法进行分类，同样分为高频和低频2 部分，高频系数部分采用自适应卡尔曼滤波算法对高频系数进行滤波。

6）低频系数滤波。低频系数部分主要包含了监测点的趋势信息，因此低频系数部分采用滑动平均算法进行简单滤波，去除部分噪点，最大限度保留监测点的趋势信息，算法公式为

7）中频系数改正。中频系数改正采用中误差原则进行改正，首先计算中频系数的中误差Lσ ，然后按如下公式进行改正，即

式中：Mi为第i 个中频系数；Lσ 为中频系数中误差；为第i 个中频系数改正值。

8）监测数据恢复。利用处理后各层的系数恢复监测数据。

2 实验与结果分析

2.1 数据来源

本文采用的监测数据来自广西阳朔至鹿寨高速公路边坡2 号监测点（如图3 所示），该边坡高约60 m，五级边坡，植被茂密，边坡后方为陡坡向上。数据采集设备为本单位自己研发的BDS 高精度接收机，解算软件亦为本单位自己研发的高精度变形监测后处理软件，数据解算间隔为1 h，数据时间为2020-02-01—2020-05-19，另外在数据中加入3 个随机粗差来检验粗差剔除效果。

图3 监测点概况

2.2 处理结果

1）粗差修复。数据的小波变换函数选择“db08”，分解层数为6 层，a6 为第6 层分解后的近似系数，计算得到的各层频带范围如表1 所示。

表1 频带范围表单位：Hz

计算得到上文所述高频、中频和低频的阈值为0.041 7 Hz，高频层为d1、d2、d3，中频层为d4，低频层为d5、d6。利用低频层进行小波逆变换，恢复监测数据，并进行粗差识别和修复，修复效果如图4所示。

图4 粗差修复效果

从图4 可以看出，3 处粗差数据都能够得到较好修复，修复效果良好。

2）细节系数滤波。将上一步骤经过粗差修复后的数据进行2 次小波变换，然后对各项系数进行改正，改正效果如图5 所示。

图5 各层细节系数滤波效果

从图5 各层细节系数滤波效果图中可以看出：d1、d2、d3、d4、d5及d6各层的细节系数都有比较大的改正，其中d1、d2、d3及d4的改正效果最好，系数的频率和最大、最小值都有一定的减小；d4周期信号经过改正以后，振幅明显减小，周期性更加明显；d5和d6层改正前后细节信息差别不是很大，主要原因是为了保留监测点位移信息仅采用了滑动平均方法进行滤波，能够看到滤波前后d5和d62 层的整体趋势都没有改变，但是改正后的曲线更加平滑，细节抖动去除明显。

3）数据恢复。利用滤波后的各层系数进行小波逆变换，变换后的效果如图6 所示。

图6 误差削弱效果

从图6 的误差削弱效果中可以看出：处理后的观测数据时间序列的平稳性有很大提高，其中以天为周期的周期性位移信号有了较大改善，数据振幅更小；处理以后，噪声信息明显减少，可以清晰地看出点位的位移趋势，有利于更精确地分析监测站的位移趋势；在3 月30 日以后监测点发生20 mm 左右位移，经过小波变换处理以后监测点的位移曲线更加明显，位移量与处理前保持一致，位移开始时间并没有延迟；另外处理后的监测点数据3 月20—30 日之间的小幅度位移也得到了有效保留；统计了3 月20 日以前原始数据数据的中误差为3.58 mm，处理后数据的中误差降为0.92 mm；从上述分析可以看出本文小波变换方法处理后的监测数据时间序列的平稳性和数据质量都有明显提高。