水田整地是水稻种植前的重要环节，精细的搅浆平地技术可以使稻田水层深浅均匀，节约农业灌溉用水，满足水稻种植“寸水不过田，寸水不漏泥”的要求；可以提高化肥、农药和除草剂的利用率，减少农业生产成本，促进水稻生产可持续发展；可以提高撒播、抛秧和插秧的作业质量，提高劳动生产率，是保证水稻高产、稳产的重要措施^[1-4]。

全球卫星导航系统(Global navigation satellite system, GNSS)控制的水田搅浆平地技术是目前先进的土地平整技术，相较于激光平地技术，GNSS技术精细平整，受风雪、雨雾等外界环境影响小，易于进行大面积水田作业，能够实现地形地势的测量^[5]。但是受多路径效应误差、拖拉机和搅浆机振动误差的影响，GNSS接收到的高程数据存在不可避免的偏差，导致水田地形的测量结果与真实情况不一致，影响基准面的准确度，进而降低水田搅浆平地机的作业质量。

为提高GNSS高程数据的精度与可信度，有必要对地形测量过程中采集到的数据进行降噪处理。夏友祥等^[6]分析了旱地农田平整系统GNSS定位数据误差源，分别利用卡尔曼和小波变换滤波对多路径效应和振动误差进行处理，提高了农田平整GNSS定位精度；梁友斌等^[7]采用小波硬阈值和卡尔曼融合滤波算法对水田旋耕平地机GNSS倾角信号进行去噪处理，减少了倾角信号中的噪声成分。虽然小波变换降噪和卡尔曼滤波在GNSS数据处理方面有一定的效果，但是小波基的选择不同会导致降噪效果不一样，小波基的选择原则需要进一步研究，且采用卡尔曼滤波算法需要建立精确的噪声模型才能取得较优的滤波效果，而对于水田平地机GNSS高程数据来说很难建立准确的噪声模型。胡炼等^[8]通过限幅加权递推平均滤波算法减少了GPS高程数据误差，滤波后的数据可以满足平地要求，但用于滤波的GPS高程定位数据的个数和权值的选择缺少理论支撑。由于经验模态分解(Empirical mode decomposition，EMD)算法在处理非线性、非平稳信号中具有优势，近年来其用于降噪的研究逐渐增多。戴吾蛟等^[9]根据高斯白噪声信号的本征模态函数(Intrinsic mode functions，IMFs)分量的能量密度与其平均周期的乘积为一常量这一特性来确定分解级数，并将其用于GPS多路径效应研究，有效削弱了随机误差，得到了更准确的多路径效应改正模型。罗飞雪等^[10]分别利用EMD低通滤波去噪法和小波分解对动态变形监测的GPS数据进行滤波处理，验证EMD算法在GPS动态变形监测数据去噪中的可行性与有效性。蒋卓芸^[11]利用EMD去噪法对高光谱数据去噪，效果良好。以上EMD降噪方法将信号分解为若干个频率由高到低的IMFs，其中，高频分量一般包含的是信号的噪声成分，将其剔除并将剩余IMF分量重构即完成对信号的滤波。然而中频分量中仍可能含有噪声成分，因此传统的EMD降噪方法在一定程度上会造成滤波不彻底，信号仍然含有少量噪声^[12]。

为解决水田平地机测量地形时GNSS高程数据存在噪声进而影响基准面建立准确性的问题，本文分析了基于GNSS控制的水田搅浆平地机测量高程数据时的误差来源，在传统EMD滤波的基础上提出了一种EMD与S-G联合降噪的方法，S-G滤波器同样适用于分析非平稳信号，且去除噪声的同时可以最大限度保持原始信号波形峰值的高度和宽度。联合降噪方法首先采用EMD将原始高程信号分解成若干IMFs，利用归一化自相关函数和相关系数将其细分为噪声IMFs、混合IMFs和有效IMFs，然后利用S-G算法对混合IMFs进行滤波，最后将S-G滤波后的IMFs与有效IMFs进行重构，得到最终降噪后的数据。通过静态试验与拉绳位移传感器测量数据进行对比，验证算法的可行性，并进行了田间试验进一步说明联合滤波算法的有效性，以期减小水田平地机在测量地势时存在的多路径效应和机器振动影响，从而提高水田平地作业质量。

1.   系统平台

基于GNSS的水田搅浆平地机的系统硬件组成如图1所示，系统主要由悬挂架、搅浆机、平地机构、液压系统、自动控制系统5部分组成。悬挂架和拖拉机三点悬挂装置挂接，平地机构设置于搅浆机后方，平地机构为平行四边形结构，液压系统固定于悬挂架上方，自动控制系统主要包括RTK GNSS系统、集成控制柜和触摸显示屏，其中，RTK GNSS系统包括基站、移动站接收机、主副卫星接收天线和差分信号天线。

图  1  系统硬件组成

1：拖拉机；2：触摸显示屏；3：差分信号天线；4：集成控制柜；5：液压集成阀块；6：悬挂架；7：搅浆机；8：平地机构；9：移动站卫星接收天线；10：基站

Figure  1.  System hardware components

1: Tractor; 2: Touch screen; 3: Differential signal antenna; 4: Integrated control cabinet; 5: Hydraulic integrated valve block; 6: Suspension frame; 7: Stirrer ; 8: Leveling mechanism; 9: Mobile station satellite receiving antenna; 10: Base station

下载: 全尺寸图片 幻灯片

使用120马力的拖拉机作为牵引动力，采用连云港双亚机械有限公司生产的搅浆机(作业幅宽为3 m)，试验中使用的导航系统是上海华测导航公司生产的P3-DT北斗高精度定位测向接收机，水平和垂直测量精度分别为(0.80±0.01)和(1.50±0.01) cm，接收数据频率为10 Hz。选用深圳市米兰特科技有限公司的MPS拉绳位移传感器，量程为1000 mm，输出信号为RS485数字信号，线性精度为 ± 0.1%，重复性精度为 ± 0.01%，在程序中设置GNSS和位移传感器同步接收数据。

2.   误差分析与算法描述

2.1   误差来源

影响全球卫星导航系统定位精度的原因有多种，如卫星星历误差、卫星钟差、电离层与对流层折射误差、多路径效应等。其中，前4种相关性较强的误差可以通过差分技术和误差纠正模型消除，但相关性较弱的多路径效应难以通过上述方法消除^[13-14]。因此，多路径效应是目前影响GNSS高定位精度的主要误差源^[15]。

基于GNSS自动控制的水田搅浆平地机，移动站卫星接收天线通过强磁铁吸附于平地机构两端，拖拉机和搅浆机启动后产生的无规则振动噪声会叠加到接收机采集的高程数据中，引起测量误差^[16]。为了提高系统的精度，需要对多路径效应引起的数据漂移和机械振动产生的随机噪声进行滤波处理。

2.2   算法及原理

2.2.1   EMD算法

经验模态分解(Empirical mode decomposition，EMD)是一种自适应的时频分析方法，该方法根据信号自身的特点将信号分解成若干个频率由高到低的本征模态函数(Intrinsic mode functions，IMFs)和1个残余分量，其中，高频IMFs分量一般为信号的尖锐部分，即噪声；低频IMFs分量一般为信号的有效成分。EMD算法适用于分析非线性、非平稳GNSS高程信号 ^[17-19]。

提取本征模态函数步骤如下^[20]：找出原始信号x(t)的所有极大值点和极小值点，利用3次样条函数拟合极大值点和极小值点分别得到上包络线e_max和下包络线e_min，计算均值包络m₁(t)、原信号x(t)与均值包络m₁(t)的差值h₁(t)。

判断h₁(t)是否满足以下2个条件：一是其极值点个数和过零点个数相差不大于1；二是其任意点处的上下包络均值为0。若不满足，则重复上述步骤，直至得到满足IMFs条件的模态分量c₁(t)。

用原信号x(t)减去c₁(t)，得到1个去除高频成分的新信号r₁(t)：

对r₁(t)重复上述过程，得到第2个IMF分量c₂(t)，反复进行，直到最后一阶分量(残余分量)至多存在1个极值点，经验模态分解结束。最终原始信号x(t)被分解为n个IMFs分量和1个残余分量：

2.2.2   S-G算法

S-G滤波法连续选取噪声信号$x[n]$中2M+1个数据(滤波窗口)，利用一个阶次为N的多项式进行拟合，多项式的系数根据最小二乘法拟合误差最小原则确定，将多项式曲线在滤波窗口中心点处的值作为降噪后的值，从左到右依次移动窗口，直至完成所有数据的滤波^[21-22]。构造拟合多项式：

拟合残差为：

令$\varepsilon $对${a_k}$的所有偏导数为0，即：

化简得：

引入矩阵$A = \left\{ {{a_{n,i}}} \right\}$，其中${a_{n,i}} = {n^i}, - M \leqslant n \leqslant M, 0 \leqslant i \leqslant N$，再引入辅助矩阵$B = {A^{\rm{T}}}A$，令：

可得：

式中，$H$第1行即为所求系数。

2.2.3   EMD与S-G联合滤波算法

EMD与S-G联合滤波算法如图2所示，主要分为3个步骤。

图  2  EMD与S-G联合滤波算法流程

Figure  2.  Flow chart of EMD and S-G combined filter algorithm

下载: 全尺寸图片 幻灯片

利用EMD算法将含有噪声的高程数据分解为若干IMFs和1个残余分量；通过求解IMFs的归一化自相关函数以及与原始信号的相关系数，将分解得到的IMFs分为噪声IMFs、混合IMFs和有效IMFs；将噪声IMFs舍弃，对混合IMFs进行S-G滤波，最后将S-G滤波后的“干净”信号和有效IMFs进行重构得到最终降噪高程信号。

本文利用随机噪声和一般信号的归一化自相关函数特性的不同来区分噪声信号和其他信号的IMFs分量界限^[23]。如图3所示，随机噪声信号在零点处归一化自相关函数值最大，在其他点迅速衰减至0附近，而对于一般信号，其最大值也出现在零点，但由于信号间存在相关性，在其他点并没有迅速衰减到很小的值。

图  3  不同信号归一化自相关函数

Figure  3.  Normalized autocorrelation function for different signals

下载: 全尺寸图片 幻灯片

本文利用相关系数(r) 区分有效IMFs与其他IMFs，通过计算每个IMF分量与原始信号的相关系数，把第1个r大于0.5及后面的低频IMFs作为有效IMFs^[24]。r表达式为：

式中，${{\rm{cov}}} ({\tilde x_i},{x_i})$表示$ {\tilde x_i} $与$ {x_i} $的协方差，$ {\sigma _{\tilde x}} $和$ {\sigma _x} $分别表示${\tilde x_i}$与$ {x_i} $的标准差。

3.   试验与结果分析

3.1   静态验证试验

3.1.1   试验方案

静态试验前，将GNSS基站放置于田间地头，移动站卫星接收天线对称磁吸于平地机构两侧，差分信号天线吸附于液压阀块顶部，位移传感器安装在调高油缸上，如图4所示。采用静态调高性能测试试验采集的数据进行验证，在调节过程中，GNSS和线性位移传感器同时记录数据，将线性位移传感器采集的数据经过转换后的高度作为真值，若经过联合滤波后的GNSS高程曲线和线性位移对应的高程吻合度高，则说明EMD与S-G联合滤波可以有效地剔除GNSS原始高程数据中的多路径效应和拖拉机振动的影响。

图  4  静态验证试验平台

Figure  4.  Static verification test platform

下载: 全尺寸图片 幻灯片

3.1.2   数据算法处理

根据“2.2”，首先将主副卫星接收天线的原始中间点高程数据采用EMD算法进行分解得到7个IMFs分量与1个趋势项，如图5所示；然后对EMD分解的7个IMFs 分量进行归一化自相关函数求解，结果如图6所示，从图6可以看出IMF1和IMF2为噪声分量，将其剔除；其次分别计算IMF1~IMF7分量与原始高程信号的相关系数，结果依次为0.1162、0.1221、0.2024、0.1876、0.2620、0.6159和0.6575，由计算结果可知，IMF6和IMF7为有效分量，对IMF3~IMF5进行S-G滤波，最后将滤波后的IMF3~IMF5分量与有效IMF分量及残余分量重构得到联合滤波后的高程信号。为更好地说明联合滤波的有效性，另外设置单独采用EMD滤波方案和单独采用S-G滤波作为对比。3种方案的降噪结果如图7所示。

图  5  原始高程数据进行EMD分解

Figure  5.  EMD decomposition of raw elevation data

下载: 全尺寸图片 幻灯片

图  6  IMFs的归一化自相关函数

Figure  6.  Normalized autocorrelation function of IMFs

下载: 全尺寸图片 幻灯片

图  7  3种方案降噪结果对比图

Figure  7.  Denoising results of three schemes

下载: 全尺寸图片 幻灯片

3.1.3   结果分析

由图7a、7b可以看出单独使用EMD和S-G滤波虽然可以抑制一部分多路径效应误差和振动误差，但数据仍有较多噪声，信号依然存在“毛刺”现象，误差较大。由图7c可以看出EMD与S-G联合滤波降噪较充分，数据最为平滑，大大提高了GNSS高程数据的精度。选用均方根误差(Root mean square error，RMSE)和信噪比(Signal to noise ratio，SNR)进行滤波效果评价，结果如表1所示。

表  1  不同算法的降噪效果统计

Table  1.  Statistics of noise reduction effect of different algorithms

算法 Algorithm	RMSE/mm	SNR/dB
滤波前 Before filtering	6.26	63.4562
EMD滤波 EMD filter	4.61	66.1109
S-G滤波 S-G filter	4.24	66.8486
联合滤波Combined filter	3.95	67.4653

下载: 导出CSV 

| 显示表格

RMSE是指滤波后信号与真实信号方差的平方根，用来衡量滤波后信号与真实信号之间的偏差，RMSE越小表示降噪效果越好，其表达式为：

式中，$ {x_i} $表示真实信号，$ {\tilde x_i} $表示滤波后信号。

SNR是指真实信号能量与噪声能量的比值，信噪比愈高表示降噪效果愈好，其表达式为：

式中，P_s、P_n分别表示真实信号功率和噪声信号功率。

由表1可知经EMD与S-G联合滤波后的均方根误差最小，信噪比最大，其次为S-G滤波方法，EMD相较于前2种滤波方法的降噪效果差。表1说明经联合降噪后的高程信号与拉绳位移传感器高程波形保持了较高的相似，联合滤波具有较优的效果。

3.2   田间试验结果

在江苏省灌云县五图河农场选取一块面积为20 m $ \times $ 50 m的水田，使用GNSS控制的水田搅浆平地机进行水田高程信息的采集。作业前，将设备按照“3.1.1”的方法安装，测量水田耕作层高程时，通过改变液压系统中电磁阀阀芯的位置将调高油缸设置为浮动模式，采用蛇形行走的方式进行离散点高程数据的采集，在地头拐弯处不采集数据。对高程数据滤波分幅宽进行，限于篇幅，只选择第1幅宽数据进行分析，结果如图8 所示。

图  8  第1幅宽原始高程数据和联合滤波处理效果图

Figure  8.  The raw elevation data and the effect of combined filtering of the first width

下载: 全尺寸图片 幻灯片

由图8 分析可知，经EMD与S-G滤波联合降噪后的高程数据得到较优平滑效果的同时保持了与原始数据几乎相同的变化趋势，说明噪声去除充分，且数据没有出现失真现象，信号中的有效信息得到保留。联合滤波后最大值与最小值均有所改善，最大值由17.107 m下降到17.100 m，最小值由17.006 m上升到17.011 m，高程差由10.1 cm下降到8.9 cm，数据波动范围减小了11.9%。通过对田间测量数据的分析进一步说明了 EMD和S-G滤波联合对水田平地机GNSS高程信号降噪的优越性。

4.   结论

本文基于水田平地机测量的GNSS实际高程数据，分析了GNSS高程数据误差来源，提出了一种EMD与S-G联合降噪的方法削弱多路径效应和振动对定位精度的影响，更加准确地反映水田信息。

静态验证试验联合滤波后的RMSE为3.95 mm，比滤波前降低了36.9%，SNR为67.4653 dB，比滤波前提高了6.3%。田间测量试验滤波后的数据波动范围减少了11.9%。试验结果表明，本文所提的EMD与S-G联合滤波算法有效削减了多路径效应误差和振动误差，改善了数据的平滑度，对于水田平地前基准面的建立及提高水田平地作业质量具有现实意义。