小波包消噪方法分析及改进

书书书第２８卷第４期２００８年８月大地测量与地球动力学ＪＯＵＲＮＡＬＯＦＧＥＯＤＥＳＹＡＮＤＧＥＯＤＹＮＡＭＩＣＳＶｏｌ．２８Ｎｏ．４　Ａｕｇ．，２００８　　文章编号：１６７１５９４２（２００８）０４０１０２０６小波包消噪方法分析及改进曲国庆１，２）　党亚民３）　章传银３）　苏晓庆２）１）山东科技大学，青岛　２６６５１０２）山东理工大学建筑工程学院，淄博　２５５０４９３）中国测绘科学研究院，北京　１０００３９摘　要　对小波包消噪中的两个关键点———阈值的选择和小波包系数的阈值量化进行了分析，并对阈值选择方法进行了改进。通过仿真试验，用信号消噪前后的信噪比和均方误差两个指标对各种方法进行对比分析。结果表明，改进后的阈值选择方法优于其他方法。关键词　变形信号；小波包；阈值；信号消噪；信噪比中图分类号：Ｐ２０７　　　　文献标示码：ＡＡＮＡＬＹＳＩＳＡＮＤＩＭＰＲＯＶＭＥＮＴＯＦＤＥＮＯＩＳＩＮＧＭＥＴＨＯＤＦＯＲＷＡＶＥＬＥＴＰＡＣＫＥＴＱｕＧｕｏｑｉｎｇ１，２），ＤａｎｇＹａｍｉｎ３），ＺｈａｎｇＣｈｕａｎｙｉｎ３）ａｎｄＳｕＸｉａｏｑｉｎｇ２）１）ＳｈａｎｇｄｏｎｇＵｎｉｖｅｒｓｉｔｙｏｆＳｃｉｅｎｃｅａｎｄＴｅｃｈｎｏｌｏｇｙ，Ｑｉｎｇｄａｏ　２６６５１０２）ＳｈａｎｄｏｎｇＵｎｉｖｅｒｓｉｔｙｏｆＴｅｃｈｎｏｌｏｇｙ，Ｚｉｂｏ　２５５０４９３）ＣｈｉｎｅｓｅＡｃａｄｅｍｙｏｆＳｕｒｖｅｙｉｎｇａｎｄＭａｐｐｉｎｇ，Ｂｅｉｊｉｎｇ　１０００３９Ａｂｓｔｒａｃｔ　Ｗａｖｅｌｅｔｐａｃｋｅｔａｎａｌｙｓｉｓｉｓａｐｒｅｃｉｓｅａｎａｌｙｓｉｓｍｅｔｈｏｄ．Ｔｈｅｒｅａｒｅｔｗｏｋｅｙｐｏｉｎｔｓ—ｃｈｏｏｓｉｎｇｔｈｒｅｓｈｏｌｄａｎｄｑｕａｎｔｉｔａｔｉｎｇｔｈｒｅｓｈｏｌｄｏｆｗａｖｅｌｅｔｐａｃｋｅｔｃｏｅｆｆｉｃｉｅｎｔｓ—ｉｎｗａｖｅｌｅｔｐａｃｋｅｔｄｅｎｏｉｓｉｎｇｍｅｔｈｏｄ．Ｉｎｔｈｉｓｐａｐｅｒ，ｔｈｅｔｗｏｐｏｉｎｔｓｉｎｄｉｆｆｅｒｅｎｔｍｅｔｈｏｄｓａｒｅａｎａｌｙｚｅｄａｎｄｃｏｍｐａｒｅｄ，ａｎｄｔｈｅｃｈｏｏｓｉｎｇｔｈｒｅｓｈｏｌｄｍｅｔｈｏｄｉｓｉｍｐｒｏｖｅｄ．Ｔｈｒｏｕｇｈｓｉｍｕｌａｔｉｖｅｔｅｓｔ，ｃｏｍｐａｒｉｎｇＳＮＲａｎｄＭＳＥｏｆｏｒｉｇｉｎａｌｓｉｇｎａｌａｎｄｔｈａｔｏｆｓｉｇｎａｌａｆｔｅｒｄｅｎｏｉｓｉｎｇｔｏｅｖａｌｕａｔｅｔｈｅｅｆｆｅｃｔｏｆｔｈｏｓｅｍｅｔｈｏｄｓ，ｉｔｉｓｐｒｏｖｅｄｔｈａｔｔｈｅｉｍｐｒｏｖｅｄｍｅｔｈｏｄａｔｔａｉｎｓｇｏｏｄｒｅｓｕｌｔｓａｎｄｉｓｓｕｐｅｒｉｏｒｔｏｏｔｈｅｒｍｅｔｈｏｄｓ．Ｋｅｙｗｏｒｄｓ：ｄｉｓｔｏｒｔｉｏｎｓｉｇｎａｌ；ｗａｖｅｌｅｔｐａｃｋｅｔ；ｔｈｒｅｓｈｏｌｄ；ｓｉｇｎａｌｄｅｎｏｉｓｉｎｇ；ＳＮＲ（ｓｉｇｎａｌｎｏｉｓｉｎｇｒａｔｉｏ）１　前言在现代变形测量技术中，由于全球定位系统技术的快速发展，提供了大范围、长时间、甚至不间断的动态观测，观测精度不断提高，观测环境愈来愈复杂，动态观测值中包含的系统误差与偶然误差相比数值较大，影响因素较多，函数关系复杂。在ＧＰＳ变形监测信号的多分辨分析中，小波方法只是将信号的低频部分进行分解，没有对高频部分进行剖分，信号的分析是一个不完整的树状结构。小波包分析则给出了一种更为精细的分析方法，不仅对信号的低频部分进行分解，同时也能够对高频部分进行细分，所以小波包分析是一个完整的树状结构，可得到更详细的信号特征。收稿日期：２００８０４０５基金项目：山东省自然科学基金（２００４ＸＺ３１）；中国科学院百人计划项目作者简介：曲国庆，男，１９６２年生，教授，山东科技大学博士研究生，从事测量数据处理教学与研究工作．Ｅ－ｍａｉｌ：[e-mail]　第４期曲国庆等：小波包消噪方法分析及改进许多学者利用小波分析法来检测ＧＰＳ相位观测值整周跳变［１］，利用小波分解不同层次间、信号之间等均具有一定的的相关性以及ＧＰＳ信号在传播过程中受到多种因素的影响，不同的因素呈现出一定的周期性等特点，提取信号特征［２－４］；利用小波分析方法分析处理ＧＰＳ信号多路径误差、相位周跳检测［５］、信号去噪［６，７］等。目前这些研究成果在一定程度上提高了数据处理的精度，取得了良好的效果。但在现代大地测量中，由于观测条件越来越复杂，要求的精度越来越高，信号消噪问题尚未得到很好的解决。小波包降噪方法有两个关键点：阈值的选择和小波包系数的阈值量化，从某种程度上说，它直接关系到信号消噪的质量。本文对小波包消噪中的关键点进行了分析比较，针对ＧＰＳ变形监测信号特点，对ｐｅｎａｌｔｙ阈值选择方法进行了改进，通过仿真试验进行了分析，并得到良好的效果。２　小波包分解与重构基本理论小波包降噪步骤为：信号的小波包分解、确定最佳小波包基、小波包分解系数的阈值量化、小波包信号重构。信号的ｆ（ｔ）正交小波分解的公式为［８］：Ｐｊ－１ｆ（ｔ）＝Ｐｊｆ（ｔ）＋Ｄｊｆ（ｔ）其中Ｐｊ－１ｆ（ｔ）＝∑ｋｘ（ｊ）ｋｊｋ（ｔ），Ｄｊｆ（ｔ）＝∑ｋｄ（ｊ）ｋφ（ｔ）正交小波的重建公式为：ｘ（ｊ）ｎ＝∑ｋｈ０（ｎ－２ｋ）ｘ（ｊ＋１）ｋ＋∑ｋｈ１（ｎ－２ｋ）ｄ（ｊ＋１）ｋ＝∑ｋｇ０（ｎ－２ｋ）ｘ（ｊ＋１）ｋ＋∑ｋｇ１（ｎ－２ｋ）ｄ（ｊ＋１）ｋ小波包的重构公式为：ｄｊ，ｎｌ＝∑ｋ［ｈ０（ｌ－２ｋ）ｄｊ＋１，２ｎｋ＋ｈ１（ｌ－２ｋ）ｄｊ＋１，２ｎ＋１ｋ］＝∑ｋｇ０（ｌ－２ｋ）ｄｊ＋１，２ｎｋ＋∑ｋｇ１（ｌ－２ｋ）ｄｊ＋１，２ｎｋ３　阈值选择准则分析比较及改进３．１　阈值选择准则的分析比较小波包降噪方法有两个关键点：阈值的选择和小波包系数的阈值量化。从某种程度上说，它直接关系到信号消噪的质量。如果阈值太小，去噪后的信号仍然有噪声；相反，阈值太大，重要的信号特征又被过滤掉，引起偏差。从直观上看，对于给定的小波系数，噪声越大，阈值也越大。１）通用阈值（ｓｑｔｗｏｌｏｇ准则）设含噪信号ｆ（ｔ）在尺度１～ｍ（ｍ＞Ｊ）上通过分解得到小波系数的个数总数和为ｎ，Ｊ为二进制尺度参数，附加噪声信号的标准差是σ，则通用阈值为：Ｔ１＝σ２ｌｎ（ｎ槡）。该方法的原理依据是Ｎ个具有独立同分布的标准高斯变量中的最大值小于Ｔ１的概率随着Ｎ的增大趋于１。若被测信号含有独立同分布的噪声时，经小波变换后，其噪声部分的小波系数也是独立同分布的。如果具有独立同分的噪声经小波分解后，它的系数序列长度Ｎ很大，则根据上述理论可知：该小波系数中最大值小于Ｔ的概率趋于１，即存在一个阈值Ｔ，使得该序列的所有小波系数都小于它。小波系数随着分解层次的加深，其长度也越来越短，根据Ｔ的计算公式，可知该阈值也越来越小。因此在假定噪声具有独立同分布特性的情况下，可通过设置简单的阈值来消除噪声。２）ｓｔｅｉｎ无偏风险阈值Ｔ２（ｒｉｇｒｓｕｒｅ准则）它是一种基于史坦的无偏似然估计（二次方程）原理的自适应阈值选择。对一个给定的阈值ｔ，得到它的似然估计，再将非似然ｔ最小化，得到所选的阈值，它是一种软件阈值估计器。具体的阈值选择规则为：设Ｗ为一向量，其元素为小波系数的平方并按由大到小的顺序排列，即Ｗ＝［ｗ１，ｗ２，…，ｗｎ］，且ｗ１≤ｗ２≤…ｗｎ，ｎ的含义同上。再设一风险向量Ｒ，其元素为：ｒｉ＝［ｎ－２ｉ－（ｎ－ｉ）ｗｉ＋∑ｉｋ＝１ｗｋ］／ｎ，ｉ－１，２，…，ｎ以Ｒ元素中的最小值ｒｂ作为风险值，由ｒｂ的下标变量ｂ求出对应的ｗｂ，则阈值Ｔ２为：Ｔ２＝σｗ槡ｂ３）启发式的ｓｔｅｉｎ无偏风险阈值Ｔ３（ｈｅｕｒｓｕｒｅ准则）启发式的ｓｔｅｉｎ无偏风险阈值是前两种阈值的综合，是一种良好的预测变量阈值选择方法。当满足某一条件时，选取阈值用固定阈值准则，否则，取无偏风险估计准则与通用阈值准则的较小者作为本准则的阈值。如果信噪比很小，估计有很大的噪声时，则采用这种固定的阈值。具体的阈值选取规则为：设Ｗ为ｎ个小波系数的平方和，令：η＝Ｗ－ｎｎ及μ＝（ｌｏｇ２ｎ）３２槡ｎ，则Ｔ３＝Ｔ１η＜μｍｉｎ（Ｔ１，Ｔ２）η＞{μ４）最大最小准则阈值Ｔ４（ｍｉｎｉｍａｘ准则）该准则采用的也是一种固定阈值，它产生一个最小均方误差的极值。在统计学上，极小极大原理常用来设计估计器，因为去噪信号可以看作与未知回归函数的估计式相似，则极小极大估计量可实现３０１大地测量与地球动力学２８卷在最坏条件下最大均方误差最小。具体的阈值选取规则为：Ｔ４＝σ（０．３９３６＋０．１８２９ｌｏｇ２ｎ）ｎ＞３２０ｎ{＜３２σ＝ｍｉｄｄｌｅ（Ｗ１，ｋ，０≤ｋ≤２ｊ－１－１）／０．６７４５式中，ｎ为小波系数的个数，σ为噪声信号的标准差；Ｗ１，ｋ表示尺度为１的小波系数。因此，式中σ的分子部分表示对分解出的第一级小波系数取绝对值后再取中值［１０］。当信号只有少量的高频系数位于噪声范围之内时，Ｍｉｎｉｍａｘ和ｒｉｇｒｓｕｒｅ阈值选择规则更加保守方便，仅将部分系数置为零，不容易丢失真实信号成分，因此在信号的高频信息有很少一部分在噪声范围内时，这两种阈值选择规则非常有用，可以将弱小的信号提取出来。另外两种阈值选取准则将系数全部置零，可以更有效地去除噪声。但也有可能将有用的高频信息当作噪声去除掉［１０］。５）Ｂｉｒｇｅｍａｓｓａｒｔ阈值，通过如下的规则求得：给定一个指定的分解层数ｊ，对ｊ＋１以及更高层，所有系数保留；对第ｉ层（１≤ｉ≤ｊ２Ｎ），保留绝对值最大的ｎｉ个系数，ｎｉ由下式确定：ｎｉ＝Ｍ（Ｊ＋２－ｉ）α。式中Ｍ和α为经验系数，缺省情况下Ｍ＝Ｌ（１），也就是第一层。一般情况下分解后系数的长度Ｍ满足Ｌ（１）≤Ｍ≤２Ｌ（１），α的取值因用途不同，在压缩情况下一般取α＝１．５，降噪情况下α＝３［１１］。６）采用Ｂｉｒｇｅｍａｓｓａｒｔ惩罚方法得到阈值，称ｐｅｎａｌｔｙ阈值，表示如下：ｃｒｉｔ＝－∑ｋ≤ｔｃ（ｋ２）＋２ｓｉｇｍａ２ｔ（Ａｌｐｈａ＋ｌｏｇ（ｎ／ｔ））式中，ｃ（ｋ）是小波包系数，它是按其绝对值递减的顺序存储的，ｎ是系数的个数，ｓｉｇｍａ＝Ｄｅｔ／０．６７４５，Ｄｅｔ为细节小波系数的绝对值中值。设ｔ是上式的极小值，那么小波包阈值即为ｃｔ，ｓｉｇｍａ值可以是零均值高斯白噪声的标准差，也可以根据第一层高频系数来估计噪声的标准差。Ａｌｐｈａ是调整参数，它必须是大于１的实数。其值越大，降噪信号的小波包表示越稀疏。小波包方法消噪采用的ｐｅｎａｌｔｙ阈值。为了避免小波系数计算中的边界效应，噪声水平标准差σ的估计是由细节小波系数的绝对值中值计算得到的，是一种鲁棒估计。当信号足够规则，尺度上的小波系数含有信号的细节都集中在少数小波系数上时，这种处理方式是比较适用的。３．２　阈值选取方法的改进在全球定位系统的精密定位、定轨中，观测量（卫星信号的相位）中除了有卫星到地面点的距离信息外，还包含电磁波信号在电离层和对流层的时延、多路径的影响等。由此得到的信号是不规则的，非平稳的，这样小波系数的绝对值中值在表征信号的特征时可能较差，如果选用这种方法来对噪声水平进行估计，就会出现一定的偏差。有必要估计一个适当的噪声水平，使阈值既能达到消噪的目的，又能保留信号的细节信息。对于ｐｅｎａｌｔｙ阈值法，在原有噪声估计因子中，加入小波系数绝对值平均值。小波系数绝对值平均值，在某种程度上反映了信号的特征，在描述信号的特点时，具有良好的参考价值，可以作为描述系数特性的一个因子。加入此因子后，可以避免中值在噪声估计时造成过于偏大或过于偏小的缺点。改正后的噪声估计计算如下：设ｘ中为第一层小波包系数的绝对值中值，ｘ平为第一层小波包系数的绝对值平均值，噪声估计表示为：ｓｉｇｍａ改＝（ｘ中－（ｘ中－ｘ平）／２）／０．６７４５阈值为：ｃｉｒｔ＝－∑ｋ≤ｔｃ（ｋ２）＋２ｓｉｇｍａ２改ｔ（Ａｌｐｈａ＋ｌｏｇ（ｎ／ｔ））式中：Ａｌｐｈａ一般取２或３。设ｔ是上式的极小值，那么小波包阈值即为ｃｔ，其它同ｐｅｎａｌｔｙ阈值。４　实例与分析本文采用的实验信号是ｂｕｍｐｓ含噪声信号，信噪比设为７，噪声为白噪声，数据点数为２１０（图１和图２）。４．１　最优小波包基的确定小波包基库是