2019, Vol. 39
2. 地理信息工程国家重点实验室，西安市雁塔路中段1号，710054;
3. 中国电子科技集团公司第20研究所，西安市白沙路1号，710054;
4. 西安测绘研究所，西安市雁塔路中段1号，710054
基于捷联惯导的海洋重力测量系统在作业时会受到惯性器件误差、仪器振动及海浪起伏等因素的干扰，一般通过数字低通滤波器消除噪声，并提取有用的重力异常信号[1-2]。数字低通滤波器的显著优势是拥有成熟的设计和实现方法，且能获得较为满意的滤波效果[3-4]，但在剔除高频噪声的同时会丢失部分有用的重力信号，降低结果的分辨率。
小波分析是在傅里叶变换基础上发展而来的，在信号处理、图像压缩、模式识别、数字分析等领域都有广泛的应用[5-6]，但在海洋重力测量方面研究相对较少。本文采用小波阈值滤波方法对捷联式海洋重力测量系统获得的原始数据进行滤波去噪实验，探讨小波阈值滤波在海洋重力测量数据处理中的适用性，并重点分析使用不同小波、阈值及施加方式的去噪效果。
1 小波分析小波变换是一种信号时频分析方法，具有多分辨分析的特点，在时频两域都具有表征信号局部特征的能力[7]。一般认为，信号经小波变换后产生的小波系数包含有信号的重要信息，其幅值较大、数目较少，而噪声对应的小波系数幅值小、数目多。通过在不同尺度上选取合适的阈值，并将小于该阈值的小波系数置零或弱化，保留大于阈值的小波系数，从而使信号中的噪声得到有效抑制，最后经小波逆变换就可以得到滤波后的重构信号[8-9]。
2 数值计算与分析沿指定测线航行得到重力异常实验值，所用仪器为捷联式海洋重力测量系统。图 1为原始未滤波测量数据及对应的频谱，图 2为数据处理系统输出的经截止频率为0.003 3 Hz和0.005 Hz(截止波长为300 s和200s)的Butterworth低通滤波器(简写为“BW”滤波器)滤波后的重力异常数据及对应的频谱。
|
图 1 原始重力异常测线数据及其频谱 Fig. 1 Raw gravity anomaly and spectrum |
|
图 2 Butterworth滤波后重力异常测线数据及其频谱 Fig. 2 Filtered gravity anomaly and spectrum by Butterworth |
从图 1和2可以看出，重力异常信号完全湮没在高频噪声中，经过Butterworth低通滤波后高频信号被剔除，处于低频带的真实重力异常信号突显出来。图 2(a)中矩形框部分为局部测线放大图，可以看出，截止频率为0.003 3 Hz的滤波结果要比截止频率为0.005 Hz的滤波结果平滑，去除的高频噪声更多。
考虑到重力异常信号为低频信号，且重力异常数据的采样频率为1 Hz，由小波多尺度(多分辨率)性质可知，小波分解至第8层时，对应的低频分量的频率范围为0~0.003 9 Hz，高频分量的频率范围为0.003 9~ 0.007 8 Hz；分解至第9层时，对应的低频分量的频率范围为0~0.002 Hz，而高频分量的频率范围为0.002~0.003 9 Hz，恰好有部分高频分量位于有用信号的频率范围内。
选取db、sym和coif类型小波分别采用4种阈值及3种阈值函数对图 1中原始重力异常信号进行降噪处理，分别对分解至第8层和第9层的高频部分进行阈值滤波处理，其他各层强制消噪(全部置零)，并将重构后的信号与截止频率为0.003 3 Hz和0.005 Hz的Butterworth低通滤波结果进行比较。
图 3为采用不同阈值及阈值函数的去噪结果与Butterworth低通滤波结果的对比(db7小波分解至第9层)，其中矩形框部分为局部测线放大后的效果。可以看出，db7小波去噪结果与Butterworth低通滤波结果整体吻合度较好，小波去噪结果较Butterworth低通滤波结果更加平滑，说明小波阈值滤波比Butterworth滤波更容易分解出噪声成分，可以更有效地消除重力畸变。另外，db7小波分解至第9层的阈值滤波结果与经截止频率为0.003 3 Hz的Butterworth低通滤波的结果更为吻合，这是由于第9层高频分量的频率范围为0~0.003 9 Hz，对高频系数进行阈值处理后，最终的滤波结果更接近截止频率为0.003 3 Hz的低通滤波结果。因此，在实际数据处理中要根据需要恢复的重力数据的分辨率要求来确定小波分解的层次，使恢复的重力结果更加精确。其他小波(包括db5、db6、db8、db9、db10，sym5、sym6、sym7、sym8、sym9、sym10、coif3、cpif4、coif5)分解至第9层的阈值滤波结果与db7滤波结果类似，限于篇幅不再给出结果对比图。值得注意的是，算例中采用史坦无偏风险阈值和启发式阈值公式得到的阈值相同。
|
1代表固定网值, 2代表史坦无偏风险阅值.3代表启发式园值, 4代表极大极小圆值, HT代表硬阀值处理.ST代表软阅值处理, IT代表改进阅值处理(调节因子a=0.4).如1-HT表示采用固定闽值且硬闽值处理 图 3 db7小波降噪结果 Fig. 3 db7 wavelet noise reduction result |
图 4为db7小波分解至第8层的阈值滤波结果与截止频率为0.005 Hz的Butterworth低通滤波结果的对比，其中矩形框部分为局部测线放大后的效果。不难看出，db7小波分解至第8层的去噪结果与Butterworth低通滤波结果吻合度较好，小波去噪结果较Butterworth低通滤波结果也更加平滑，说明前者更容易分解出噪声成分，可以更有效地消除重力畸变。对比局部测线放大图可以发现，硬阈值降噪曲线存在附加震荡，不具有软阈值和改进阈值降噪曲线的平滑性；软阈值降噪曲线较Butterworth低通滤波结果曲线在局部存在压缩现象，改进阈值降噪曲线在一定程度上缓解了这种信号压缩现象，其他小波(包括db5、db6、db8、db9、db10，sym5、sym6、sym7、sym8、sym9、sym10、coif3、cpif4、coif5)分解至第8层的阈值滤波结果与db7滤波结果类似。
|
图 4 db7小波降噪结果 Fig. 4 db7 wavelet noise reduction result |
小波分解至第8层时，高频分量的频率范围位于截止频率0.007 8 Hz(即滤波周期128 s)以内，足够提取海洋重力异常信号。因此，针对海洋重力异常信号滤波，小波分解至第8层或第9层即可。
为便于对比，统计各种小波在不同调节因子下降噪结果与Butterworth低通滤波结果(作为参考标准)差值的RMSE。图 5给出小波分解至第8层的阈值滤波结果与截止频率为0.005 Hz的Butterworth低通滤波结果差值的统计，可以看出，对于不同的小波，分别采用不同的阈值及阈值函数对第8层高频信号分量进行处理后，db小波中db7、db8、db9、db10的降噪效果较好，sym小波中sym7、sym8、sym9、sym10的降噪效果较好，coif小波中coif4、coif5的降噪效果较好，且这些小波的去噪结果与截止频率为0.005 Hz的Butterworth低通滤波结果差值的均方根误差均在0.25 mGal以内。对比图 5(d)~5(i)可以看出，除db5、db6、sym5、sym6、coif3小波外，其他小波在不同的阈值及阈值函数下对第8层高频分量进行阈值滤波处理后，最终的去噪效果相当，改进的阈值去噪效果并无明显提升。
|
(a)~(c)为固定阈值, (d)~(D为史坦无偏/启发式风险阈值，(g)~(i)为极大极小國值 图 5 小波在不同调节因子下的RMSE统计 Fig. 5 RMSE statistics of wavelets under different adjustment factors |
图 6为小波分解至第9层的阈值滤波结果与截止频率为0.003 3 Hz的Butterworth低通滤波结果差值的统计，可以看出，对于不同的小波，分别采用不同的阈值及阈值函数对第9层高频分量进行处理后，所有小波的去噪结果与截止频率为0.003 3 Hz的Butterworth低通滤波结果差值的均方根误差均在0.25 mGal以内。db小波中db6、db7、db8、db9、db10的降噪效果较好，sym小波中sym6、sym7、sym8、sym9、sym10的降噪效果较好，coif小波中coif3、coif4、coif5的降噪效果相当。另外，同一小波采用史坦无偏/启发式风险阈值和极大极小阈值的去噪效果要优于采用固定阈值的去噪效果，且采用史坦无偏/启发式风险阈值的去噪效果最优。
|
(a)~(c)为固定阈值, (d)~(f)为史坦无偏/启发式风险阈值，(g)~(i)为极大极小阈值 图 6 小波在不同调节因子下的RMSE统计 Fig. 6 RMSE statistics of wavelets under different adjustment factors |
除db5、sym5小波外，其他小波的硬阈值滤波结果与Butterworth低通滤波结果差值的RMSE均优于软阈值和改进阈值滤波结果，这是因为Butterworth低通滤波信号本身就具有很大的震荡性，而硬阈值滤波信号会产生附加震荡，软阈值和改进阈值估计得到的小波系数整体连续性较好，滤波信号不会产生附加震荡，因此硬阈值函数在均方根误差上优于软阈值和改进阈值。虽然硬阈值滤波结果的RMSE最小，但这并不意味着其滤波效果最好；相反，改进阈值的去噪效果要优于软阈值去噪效果，且采用固定阈值和极大极小阈值时，改进阈值的滤波效果更为明显。
小波分解至第8层时，最高层低频信号分量所对应的低通滤波器的截止频率为0.003 9 Hz，与截止频率为0.003 3 Hz的Butterworth低通滤波结果较为接近。为探究小波强制去噪的效果，直接对各层高频信号分量强制去噪，并重构信号，对比重构信号与截止频率为0.003 3 Hz的Butterworth低通滤波的结果，具体见图 7。可以看出，强制去噪结果与Butterworth低通滤波结果在波形上保持高度一致，但强制去噪结果在端部出现更多的无效数据，其中db小波强制去噪结果的边界效应最为严重，sym小波次之，coif小波最小，且同一类型小波随着小波数的增加，边界效应更加严重。
|
图 7 小波强制去噪结果 Fig. 7 Wavelet forced denoising results |
为对比不同小波强制去噪效果，先去除端部无效数据，再统计强制去噪结果与Butterworth低通滤波结果差值的RMSE，结果见图 8。可以看出，采用不同的小波进行强制去噪，其结果不尽相同，db小波中db6、db7、db8、db9、db10的去噪效果较好，sym小波中sym6、sym7、sym8、sym9、sym10的去噪效果较好，coif小波中coif3、coif4、coif5的去噪效果相当，这与阈值滤波得出的结论一致。
|
图 8 小波强制去噪结果RMSE统计 Fig. 8 Wavelet forced denoising results RMSE statistics |
对比阈值滤波和强制去噪的结果不难看出，强制去噪只能对固定频率段(小波分解后最高层低频分量对应的频率范围)的小波系数进行重构，而阈值滤波将小波分解后的低频分量和高频分量综合起来进行重构，可使频带更宽，且阈值滤波结果的边界效应要远小于强制去噪结果。因此，在实际处理过程中，阈值滤波比强制去噪更有优势，结果也更为精确。
3 结语本文结合捷联式海洋重力测量系统数据，选择db、sym和coif小波在不同阈值和阈值施加方式下的滤波结果进行对比分析。结果表明，同一阈值下，不同小波阈值滤波结果不尽相同，但差异不大；小波阈值滤波结果较Butterworth低通滤波结果更加平滑，说明小波阈值滤波比Butterworth滤波更容易分解出噪声成分，可更有效地消除重力畸变；阈值滤波中db7、db8、db9、db10、sym7、sym8、sym9、sym10、coif3、coif4和coif5小波去噪效果较好；采用固定阈值或极大极小阈值时，改进的阈值滤波结果明显优于软阈值滤波结果；硬阈值滤波结果与Butterworth低通滤波结果差值的均方根误差要优于软阈值和改进阈值的滤波结果；随着小波分解层次的增加，小波阈值滤波结果的端部效应也更加严重。
总之，只要选择合适的小波、小波阈值及阈值施加方式，用小波阈值滤波对海洋重力测量数据进行滤波是有效且可行的。本文是对事先采集好的重力异常测线信号进行处理，后续还需通过大量实验来验证小波阈值滤波在海洋重力测量实时数据处理中的适用性，并对在不同海况时滤波参数的选取作进一步研究。
致谢: 感谢西安测绘研究所提供的捷联式海洋重力实测数据。
| [1] |
孙中苗.航空重力测量理论、方法及应用研究[D].郑州: 信息工程大学, 2004 (Sun Zhongmiao. Theory, Methods and Applications of Airborne Gravimetry[D]. Zhengzhou: Information Engineering University, 2004)
(  0) |
| [2] |
韦建成, 肖云, 王利, 等. 海洋重力测量中低通滤波器的设计与应用[J]. 测绘科学与工程, 2017, 37(6): 11-16 (Wei Jiancheng, Xiao Yun, Wang Li, et al. Design and Application of Low-Pass Filter in Marine Gravimetry[J]. Surveying and Mapping Science and Engineering, 2017, 37(6): 11-16)
( 0) |
| [3] |
罗锋, 郭志宏, 骆遥, 等. 航空重力数据的等波纹FIR低通滤波试验[J]. 物探与化探, 2012, 36(5): 856-860 (Luo Feng, Guo Zhihong, Luo Yao, et al. Experimental Researches on FIR Lowpass Filter Based on Equiripple[J]. Geophysical and Geochemical Exploration, 2012, 36(5): 856-860)
( 0) |
| [4] |
李宜龙, 曾敏, 孙科, 等. 航空重力测量中FIR滤波与高斯滤波的比较[J]. 海洋测绘, 2014, 34(6): 40-42 (Li Yilong, Zeng Min, Sun Ke, et al. Comparison between FIR Filter and Gaussian Filter in Airborne Gravity Surveys[J]. Hydrographic Surveying and Charting, 2014, 34(6): 40-42 DOI:10.3969/j.issn.1671-3044.2014.06.010)
( 0) |
| [5] |
罗锋, 郭志宏, 王明, 等. 基于DB小波阈值去噪的航空重力数据试验[J]. 物探与化探, 2013, 37(4): 645-654 (Luo Feng, Guo Zhihong, Wang Ming, et al. Experimental Researches on the Threshold of Airborne Gravity Data Denoising Based on DB Wavelet Transform[J]. Geophysical and Geochemical Exploration, 2013, 37(4): 645-654)
( 0) |
| [6] |
郝燕玲, 姜鑫, 周广涛, 等. 基于小波分析的海洋重力测量滤波算法[J]. 测绘科学, 2014, 39(8): 116-119 (Hao Yanling, Jiang Xin, Zhou Guangtao, et al. Application of Wavelet Analysis in Marine Gravity Survey[J]. Science of Surveying and Mapping, 2014, 39(8): 116-119)
( 0) |
| [7] |
Mallat S G. A Theory for Multiresolution Signal Decomposition: The Wavelet Representation[J]. Communications on Pure and Applied Mathematics, 1989, 41(7): 674-693
( 0) |
| [8] |
刘涛, 曾祥利, 曾军. 实用小波分析入门[M]. 北京: 国防工业出版社, 2006 (Liu Tao, Zeng Xiangli, Zeng Jun. Introduction to Practical Wavelet Analysis[M]. Beijing: National Defense Industry Press, 2006)
( 0) |
| [9] |
董长虹, 高志, 余啸海. MATLAB小波分析工具箱原理与应用[M]. 北京: 国防工业出版社, 2004 (Dong Changhong, Gao Zhi, Yu Xiaohai. Principle and Application of MATLAB Wavelet Analysis Toolbox[M]. Beijing: National Defense Industry Press, 2004)
( 0) |
2. State Key Laboratory of Geographic Information Engineering, 1 Mid-Yanta Road, Xi'an 710054, China;
3. The 20th Research Institute of China Electronics Technology Group Co, 1 Baisha Road, Xi'an 710054, China;
4. Xi'an Research Institute of Surveying and Mapping, 1 Mid-Yanta Road, Xi'an 710054, China