进入新世纪以来，国家的综合国力不断增强，国民生产总值和国家经济水平快速提升，各种大坝、桥梁、隧道、高铁等利国利民工程相继竣工使用或开工建设，这些大型土木水利建设工程由于空间跨越较大，结构较为复杂，在施工或使用阶段受到降雨、大风等外界因素的影响时会产生微小的结构形变，当这种微小的形变经过日积月累超过一定限制时，就会影响建筑物的正常使用，并带来安全隐患。因此，对建筑物进行持续、可靠和高精度的变形监测从而提前预防事故发生具有重要的现实意义^[1-3]。

当前常用的变形监测方法为利用安装在建筑物所处环境的专业测量仪器设备（如电子水准仪、GPS、倾斜仪等）对时空域中各种与建筑物变形有关的环境因素（如位移、温度、湿度等）进行采集，然后采用一定的方法对采集到的数据进行处理和分析，最后给出当前阶段建筑物安全性的评估和未来变形发展趋势的预报。受环境因素、观测手段及测量仪器自身的影响，采集到的初始变形监测数据中往往会存在较多的测量误差，即噪声。由于建筑物的变形是一个缓变过程，得到的观测数值较小，噪声的存在会影响对变形信息的提取，严重时会造成对建筑物安全性的错误评估。因此在对变形监测数据分析处理时，需要提前将数据中包含的噪声分量滤除。

针对变形监测数据的噪声抑制问题，国内外学者进行了深入的研究，文献[4-6]提出采用中值滤波、维纳滤波和卡尔曼滤波等基于时域、空间域和频域滤波方法对线性平稳时间序列中的噪声进行抑制，取得了较好的噪声抑制效果，但是该类方法在面对非线性非平稳的变形监测数据时不能获得理想的噪声抑制性能；小波变换^[7]、经验模态分解（Empirical Mode Decomposition，EMD）^[8]等方法通过对观测序列进行多尺度分解，将噪声分量和信号分量分解到不同的小波层或本征模函数，最后通过信号重构实现噪声抑制，通常可以获得较为理想的噪声抑制结果，被广泛应用于对非线性非平稳变形监测数据的分析处理中^[9]；随着人工智能和机器学习技术研究的不断深入，利用神经网络方法对变形监测数据进行处理和分析成为了当前的研究热点^[10-13]，文献[14]利用RBF神经网络的最佳逼近、全局最优和快速迭代训练的优势对混杂白噪声的变形监测数据进行处理和预测，并采用实际数据对所提方法进行了验证。上述研究都是针对变形监测数据中存在的随机误差（高斯白噪声）进行分析，而实际工程应用过程中，除了随机误差外，变形监测数据中也会存在例如多径误差、测量粗差等色噪声，色噪声的存在同样会对变形监测数据分析和处理产生严重影响，因此需要提前滤除。

本文在前述研究的基础上，针对变形监测数据中随机误差（白噪声）和多径误差（色噪声）同时存在的问题，将CLEAN算法^{[15, 16]}引入GPS变形监测数据处理领域，并根据变形监测数据中白噪声、色噪声和有用信号的分布特性差异，提出一种基于广义匹配滤波（generalized matched filtering，GMF）^[17]和CLEAN的噪声抑制算法，在实现噪声抑制的同时较好的保留了变形监测数据中的有用信息，从而提高变形监测精度，通过仿真实验和实际GPS变形监测数据的噪声抑制试验对所提方法的性能进行验证，并与传统基于小波和EMD的去噪方法进行比较。

1 算法原理

CLEAN算法最初是用来对综合孔径成像雷达的旁瓣信号进行抑制从而提升图像质量，作为一种经典的信号去耦方法，一经提出便得到了广泛关注，特别是在从复杂信号中提取特定频率分量信号领域，充分发挥了其精度高，运算量小等优势。利用CLEAN算法对信号s(t)中进行处理的步骤为^{[15, 17]}：

1）将s(t)作为算法的初始输入信号；

2）对s(t)进行快速傅里叶变换（fast Fourier transformation，FFT），将其由时域变换至频域，得到s(t)的频谱s(f)=FFT(s(t))；

3）寻找并记录s(f)中最大值点对应的幅度A_d，频率f_d和相位φ_d；

4）利用步骤（3）的结果构建信号s(t)中能量最大的谐波信号x_d(t)

$ x_{d}(t)=A_{d} \exp \left(j 2 \pi f_{d} t+\varphi_{d}\right) $

(1)

5）利用s(t)减去x_d(t)，得到去除最大谐波的剩余信号s_r(t)；

6）判断s_r(t)是否满足终止条件，不满足则令s(t)=s_r(t)，重复进行步骤2）~步骤5），直到满足迭代终止条件。

根据上述对CLEAN算法迭代步骤的描述可知，一次迭代可以提取信号中的最大频率分量，通过多次迭代可以依次将信号中的谐波分量进行提取，从而实现对信号中随机噪声的抑制。

正如前文分析，GPS变形监测数据中除了随机噪声外，多径效应、测量粗差等因素会造成监测数据中存在色噪声，因此需要对CLEAN算法进行优化，以适应白噪声和色噪声同时存在的GPS监测数据处理问题。

考虑到短时间内，GPS接收机所处环境通常不会发生剧烈变换，因此多径误差、测量粗差等色噪声在一定时间范围内是较为稳定且相关的，因此可以采用广义匹配滤波^[17]的方式对色噪声进行滤波，首先将其转变为白噪声，然后利用CLEAN算法进行滤除。

存在色噪声的GPS变形监测序列u(t)可以用以下形式表示：

$ u(t)=s(t)+g(t) $

(2)

式中，s(t)表示包含变形信息的有用信号；g(t)为多径噪声、测量粗差等色噪声。观测序列的自相关矩阵可以表示为：

$ \boldsymbol{B}_{u}=\frac{1}{N} \sum\limits_{n=1}^{N} \boldsymbol{u}_{n} \boldsymbol{u}_{n}^{\mathrm{T}}=\frac{1}{N} \sum\limits_{n=1}^{N} \boldsymbol{s}_{n} \boldsymbol{s}_{n}^{\mathrm{T}}+\frac{1}{N} \sum\limits_{n=1}^{N} \boldsymbol{g}_{n} \boldsymbol{g}_{n}^{\mathrm{T}}=\boldsymbol{B}_{s}+\boldsymbol{B}_{{g}} $

(3)

式（3）的推导过程中认为s(t)与g(t)相互独立，其中，B_u为监测信号的协方差矩阵；u_n为协方差矩阵对应的特征向量；B_s和B_g分别为s(t)和g(t)的协方差矩阵；s_n和g_n分别为对应的特征向量。利用cholesky分解对B_g进行分析可得：

$ \boldsymbol{B}_{g}=\boldsymbol{V}^{\mathrm{T}} \boldsymbol{V} $

(4)

其中，V为cholesky分解得到的下三角矩阵，将式（4）带入式（3）可得：

$ \boldsymbol{B}_{u}=\boldsymbol{B}_{s}+\boldsymbol{B}_{g} \Rightarrow \boldsymbol{V}^{-\mathrm{T}} \boldsymbol{B}_{u} \boldsymbol{V}^{-1}=\boldsymbol{V}^{-\mathrm{T}} \boldsymbol{B}_{s} \boldsymbol{V}^{-1}+\boldsymbol{I} $

(5)

式中，I为单位矩阵。

经过上述处理可以看出GPS监测数据中的色噪声被白化，其自相关矩阵被转变为单位阵。此时包含变形信息的有用信号s(t)的自相关矩阵可以表示为：

$ \overline{\boldsymbol{B}_{s}}=\sum\limits_{n=1}^{K}\left(\boldsymbol{V}^{-\mathrm{T}} \boldsymbol{s}_{n}\right)\left(\boldsymbol{V}^{-\mathrm{T}} \boldsymbol{s}_{n}\right)^{\mathrm{T}} $

(6)

通过上述分析可知，对于同时含有色噪声和白噪声的GPS变形监测数据进行噪声抑制，可以采取联合使用GMF和CLEAN的方法，首先利用GMF对观测序列进行处理，将色噪声转化为白噪声，然后利用CLEAN算法对其进行分析，提取其中的有用信息从而实现噪声抑制。

2 仿真试验

为了能够准确评估所提方法对白噪声和色噪声的抑制性能，需要构建一组真值已知的数据作为标准数据进行分析。在仿真实验中，利用MATLAB仿真产生一组振幅、频率已知的谐波信号，然后对其叠加高斯白噪声，在此基础上再利用MATLAB自带的imnoise()函数添加乘性色噪声构建测试信号，利用所提方法对其进行噪声抑制。由于变形监测数据具备典型的非平稳、非线性特征，而小波方法和EMD方法是两种经典的非平稳非线性信号噪声抑制方法，因此对于同样的含噪数据分别采用小波方法和EMD方法对其噪声抑制，将结果与所提方法进行对比，能够更好的评估所提方法的噪声抑制性能。

仿真实验中涉及的参数设置为：

1）仿真序列长度为800；

2）谐波信号：

$ s(t)=A_{1} \cos \left(2 {\rm{ \mathsf{ π}}} f_{1} t\right)+A_{2} \cos \left(2 {\rm{ \mathsf{ π}}} f_{2} t\right)+A_{3} \cos \left(2 {\rm{ \mathsf{ π}}} f_{3} t\right) $

其中f₁=1 000 Hz, f₂=900 Hz, f₃=1 200 Hz为谐波信号频率，A₁=3, A₂=1, A₃=2为谐波信号幅度；

3）利用MATLAB自带awgn()函数对谐波信号叠加高斯白噪声n(t)，得到信噪比（signal to noise ratio，SNR）为5dB的含噪信号x(t)=s(t)+n(t)，本文中信噪比定义如下：

$ S N R=10 \times \log _{10} \frac{\bar{P}_{s}}{P_{n}}=10 \times \log _{10} \frac{\frac{1}{D} \sum\limits_{t=1}^{D} P_{s(t)}}{P_{n}} $

(7)

式中，P_s为测试信号s(t)的平均功率；P_s(t)=|s(t)|²为测试信号s(t)中第t个样本的功率；P_n表示高斯白噪声的功率谱密度；D表示仿真序列长度；

4）利用MATLAB自带的imnoise()函数对含噪信号x(t)叠加乘性色噪声w(t)，得到信噪比为5 dB，信杂比（signal to clutter ratio，SCR）为5 dB的测试信号y(t)=x(t)+w(t)=s(t)+n(t)+w(t)，本文中信杂比定义如下：

$ S C R=10 \times \log _{10} \frac{\bar{P}_{s}}{P_{w}}=10 \times \log _{10} \frac{\frac{1}{D} \sum\limits_{t=1}^{D} P_{s(t)}}{\frac{1}{D} \sum\limits_{t=1}^{D} P_{w(t)}} $

(8)

式中，P_w(t)=|w(t)|²为乘性噪声w(t)中第t个样本的功率。

图 1给出了仿真得到的谐波信号s(t)，信噪比为5 dB的含噪信号x(t)和信噪比为5 dB，信杂比为5 dB的测试信号y(t)。从图 1可以看出，当谐波信号叠加高斯白噪声后，原本平滑信号被叠加了很多“毛刺”，但是从信号中仍能看到信号的周期性，再叠加色噪声后，已经不能看出信号的周期性，即信号中的有用信息被噪声污染，无法从中提取有效的信息。图 2给出了利用所提方法，小波方法和EMD方法进行噪声抑制得到的结果，每蝠图中第一行子图为噪声抑制后信号，第二行子图为噪声抑制后信号和谐波信号的残差。从图 2结果可以看出，所提方法在有效抑制高斯白噪声和乘性色噪声的同时几乎不对原始谐波信号产生任何影响，保留了信号中的细节信息，原始信号与重构信号的残差在整个采样范围内很小，噪声抑制性能远远优于小波方法和EMD方法。究其原因在于小波方法需要预设小波基函数，而基函数是脱离信号的，并且分解层数和阈值的设置也会影响噪声抑制性能；而EMD方法提取的IMF虽然是基于信号的，但是EMD方法由于模式混叠影响，通常会收敛到局部最优解。

为了定量的描述三种方法在不同信噪比和信杂比条件下的噪声抑制性能，采用式（9）所示的相关系数和式（10）所示均方误差作为评价指标，相关系数从波形的相似性维度评估噪声抑制性能，相关系数越接近于1表明滤噪声抑制后的信号与原始序列波形越相似，噪声抑制性能越好；均方误差从能量的维度评估噪声抑制性能，均方误差越小说明噪声抑制后信号与原始信号的能量越接近，噪声抑制性能越好。

表 1给出了信噪比为0 dB~20 dB，信杂比为0 dB~20 dB条件下利用三种方法进行噪声抑制得到去噪信号与原始信号的相关系数值，从中可以直观看出，在信噪比和信杂比相同的条件下，所提方法相对于其他两种方法可以获得最优的噪声抑制性能。假设能够达到实际应用所需的相关系数要优于0.85，所提方法所需实测数据的最低信噪比为5 dB，最低信杂比为0 dB即可，而小波方法所需最低信噪比为5 dB，最低信杂比为5 dB，EMD方法所需最低信噪比为10 dB，最低信杂比为5 dB，即实际应用过程中，所提方法对实测数据信噪比，信杂比的要求最低，最容易满足。

表 2给出了信噪比为0~20 dB，信杂比为0~20 dB条件下利用三种方法进行噪声抑制得到去噪信号与原始信号的重构误差，从中可以直观的看出，在信噪比和信杂比相同的条件下，所提方法相对于其他两种方法可以获得最优的噪声抑制性能。假设能够达到实际应用所需的相关系数要小于0.2，所提方法所需实测数据的最低信噪比为5 dB，最低信杂比为0 dB即可，而小波方法所需最低信噪比为5 dB，最低信杂比为5 dB，EMD方法所需最低信噪比为10 dB，最低信杂比为5 dB，即实际应用过程中，所提方法对实测数据信噪比，信杂比的要求最低，最容易满足。同时从表 1~表 2可以看出，重构误差维度的分析结果与相关系数维度一致。

$ x {cor}=\frac{\frac{1}{D} \sum\limits_{t=1}^{D} \bar{s}(t) \cdot s(t)}{\sqrt{\frac{1}{D} \sum\limits_{t=1}^{D} \bar{s}^{2}(t)} \cdot \sqrt{\frac{1}{D} \sum\limits_{t=1}^{D} s^{2}(t)}} $

(9)

$ r m s e=\sqrt{\frac{1}{D} \sum\limits_{t=1}^{D}|\bar{s}(t)-s(t)|^{2}} $

(10)

3 工程实例

实际工程应用中，对微弱变形监测数据中的测量随机误差和多径噪声进行抑制是后续数据分析与处理的基础，具有重要意义。本文采用的实测数据为某高层建筑的GPS动态监测数据，选择X方向上800期累计变形监测数据进行分析。图 3给出了实测数据的时域波形和对应的频谱，并在图中对测量仪器、方法等因素引起的白噪声和多径因素引起的色噪声进行了标注。

图 4给出了利用所提方法、小波方法和EMD方法对图 3所示实测数据进行噪声抑制的结果。每蝠插图中，第一行子图为噪声抑制后的时域信号，第二行子图为噪声抑制后的频域信号。可以看出，利用所提方法进行噪声抑制得到的时域信号更加平滑，频域信号更加“干净”，表明所提方法可以有效的实现对实测数据中随机噪声和多径噪声的抑制，并且最大程度的保留了信号中细节部分和趋势，试验结果要明显优于另外两种对比方法。同时对比图 4（b）和图 4（c）的时域信号可以看出，小波方法相对于EMD方法可以获得更好的随机噪声抑制性能，得到的时域信号“毛刺”要少一些，但是从频域信号可以看出，EMD方法的多径噪声抑制性能要更好一些，得到的频谱中多径噪声被较好抑制。

4 结束语

利用GPS对建筑物的形变进行监测时，监测数据中不可避免的会混杂着由测量仪器等产生的随机噪声和由多径效应、测量粗差等产生的色噪声，这些噪声的存在会对监测数据的测量精度产生影响，因此在对数据进行处理和分析前需要提前对噪声进行抑制。针对该问题，本文提出一种基于GMF-CLEAN的变形监测数据噪声抑制方法，通过仿真实验从相关系数和重构均方误差两个维度对GMF-CLEAN方法的噪声抑制性能进行验证，并将试验结果与传统小波方法和EMD方法进行对比，结果表明在不同信噪比的信杂比条件下，所提方法均能获得最好的噪声抑制性能，并且要获取满足实际应用要求的噪声抑制能力，所提方法对实测数据的信噪比和信噪比要求最低。最后结合某工程实际案例对所提方法的噪声抑制性能进行了验证，结果表明所提方法可以获得较好的噪声抑制性能，能够满足实际工程应用需求。