时间域航空电磁的天电噪声去除研究

贲放，黄威，路宁，韩飞，郑红闪，丁志强，李军峰

(1.自然资源部 地球物理电磁法探测技术重点实验室，河北 廊坊 065000； 2.中国地质科学院 地球物理地球化学勘查研究所，河北 廊坊 065000； 3.国家现代地质勘查工程技术研究中心，河北 廊坊 065000； 4.中水东北勘测设计研究有限责任公司，吉林 长春 130026)

0 引言

随着航空电磁正反演技术及仪器采集技术的进步，航空电磁法已广泛应用于诸多领域[1]，实现了数据密集多通道采集，采集的原始数据需要进行数据预处理为数据解释做准备。原始数据中包含了多种噪声，其中对数据影响较大的为天电噪声和运动噪声。降低天电噪声影响可以提高数据信噪比，并有利于增大勘探深度[2-4]和提高小异常的分辨能力[5]。

天电主要由大气中的带电电荷经过磁暴或者雷电作用发生电磁辐射而产生，并通过空气进行传播。磁暴引起的天电噪声幅值大且不可预测，而雷电引起的天电噪声有一定的规律。全球雷电发生的频率约在100次/秒，在电磁信号中表现为瞬时噪声(尖脉冲或振荡)，测量数据中表现为突跳现象。天电干扰的数量、频率、强度、持续时间及幅度随时间、位置和源距的不同均有所变化，对各频段数据的影响依赖于电离层情况。天电噪声的频带范围为5～100 kHz，在低频范围内幅值的变化和频率近似成反比，高频部分幅值迅速增长[6]，在500 Hz～2.5 kHz的范围内对数据影响较小，而在低于500 Hz和2.5 kHz～10 kHz的范围对数据影响较大[5]。赤道或热带地区的天电干扰比中纬度地区高200多倍，而中纬度地区的天电干扰夏天比冬天高100多倍；中午的天电干扰比早晚高10多倍[7]。在较近的勘探区，由于天电直接传播到测点导致数据受天电干扰较大，使得航空电磁测量能力下降甚至无法分辨；在较远的勘探区，天电通过电磁波在地面—电离层波导中多次反射到测点，相比较而言对数据影响小。经过反射的天电由多个闪电叠加而成[8]。天电噪声的水平极化导致其对数据垂直分量影响非常小，但由于作业过程中线圈会发生姿态变化，因此数据垂直磁场分量同样也受到天电噪声的影响[5]。在背景噪声较小地区，水平分量和垂直分量数据中分别有0.05%和 0.005%的噪声[9]，在天电噪声频繁情况下观测的数据，水平分量数据有2%含有天电噪声，垂直分量则有0.5%[4]含有天电噪声。

通过对天电特性的分析，研究去除天电噪声的方法。早在1984年，Macnae等利用裁剪法(pruning)或锥形叠加的方法对非平稳或相干噪声进行去除，达到提高信噪比的目的，其认为高频率采样的航空电磁数据量充足，且经过研究含有天电的数据量很小，裁剪掉含有天电影响的连续半周期数据不会对整体数据质量造成影响[9]。同年，Bednar和Watt提出利用α-trimmed均值/中值滤波的方法对噪声进行滤波，该方法可以对大规模数据进行快速滤波，对天电等突跳特点的噪声具有较好的滤除效果[10]。1989年，Sutarno和Vozoff[11]提出了回归M值—估计法，同时对最小平方法的最佳给出了证明。随后，Buselli和Cameron[4]利用该方法对天电噪声进行去除，同时对中值滤波、鲁棒统计法、均值滤波、M值—估计法进行了比较，经过研究当天电噪声符合高斯分布特性时，M值—估计法去除天电噪声较其他方法更好。1998年，Buselli等通过实测数据中天电和地面基站观测信号的相关性对天电噪声进行识别，阐述了局部噪声预测滤波和远噪声滤波，利用预测滤波对高频部分的天电噪声进行去除，去除过程中应用了神经网格。Cull[12]提出利用线性叠加对天电噪声进行去除，Lane等[13]通过对天电噪声的来源及特征进行分析后，选择使用非线性滤波对天电噪声进行去除。2000年，Leggatt等通过增大发射功率和利用更大计算能力的计算机来阻挡天电噪声对数据的干扰。2010年，Bouchedda等[5]利用小波提取和小波变换对天电噪声进行去除，小波提取技术是基于平稳小波变换，通过将时间域航空电磁数据分解为不同频率成分的细节和近似量[14]，在一阶小波细节系数中利用能量检测方法识别天电噪声，将数据中与天电噪声相关的小波细节系数进行去除，通过剩余小波细节系数利用反变换重构方法对信号进行重构，从而实现去除天电噪声的目的。虽然该方法可以有效地对天电噪声进行去除，但仍存在一定的不足，当天电噪声位于发射电流开、关期间内时，由于天电噪声中低频成分占主导降低了该方法对天电去除的有效性。为解决该问题，Bouchedda等提出了小波变换的方法，利用静态小波变换具有线性和位移不变性的特征，在小波域中进行小波细节系数叠加，其中利用了算术平均、中值平均及α-trimmed均值叠加技术，通过对合成数据和实测数据分析说明了该方法非常适合去除晚期道低频天电干扰。Reninger等[15]利用奇异值分解法对天电噪声进行有效去除。

对于天电噪声，去除的方法众多，各有优劣。裁剪法固然有效，但需要对天电噪声进行识别后方可剪裁，当天电噪声较多时，直接剪掉数据会影响叠加次数，从而对后期的数据预处理产生影响。直接对含有天电噪声数据进行中值滤波或均值滤波会使数据曲线中心产生偏移，例如当窗口滑动到含有天电成分的数据时，因数据窗内部含有尖峰信号，滤波后数值会远大于准确值，从而产生了因数据预处理而出现的假异常，且均值滤波是对信号进行平滑处理，对噪声去除效果不明显[16]。M值—估计法的叠加窗较短不适合天电噪声较多的情况，而实测数据中天电噪声含量是无法预测和估计的。利用神经网络对天电噪声进行去除时，需要利用不同时间、不同位置的大量数据对神经网络进行训练，严重降低了其天电噪声去除的速度[16]。小波变换和奇异值分解法去除天电噪声耗时较多。在实测数据中去除天电噪声，过程中不仅需要考虑去除效果，同时需要考虑处理速度等因素，而Bednar和Watt提出的α-trimmed均值/中值滤波方法不仅可以快速去除天电噪声，且去除效果较好。

综上，为保证原始数据量的不变性，参考α-trimmed均值/中值滤波方法进行改进，可称之为剪切均值滤波。原方法中描述不同的α(0～0.5)取值，对应了均值滤波或中值滤波，其中对[αN]进行判断和取整，α、N的取值对数据处理效果有较大影响，而剪切均值滤波直接给定滤波的窗宽(整数)、单边剪切长度(整数)，无需判断[αN]的情况进行取整，并通过对单窗口滤波参数(窗口长度，单边剪切窗宽)进行研究，提出对实测数据进行混合窗口滤波。

1 天电噪声去除原理

电磁工作的主要频率范围在5～25 kHz[8]，在实际数据采集过程中经常受到各种不同噪声的干扰，由于噪声的来源广泛，增加了实测数据预处理的难度，且不同于地面瞬变电磁，航空电磁[17]数据采集过程中多了天电噪声和运动噪声，是对数据影响较大的两种噪声。基于α-trimmed均值/中值滤波方法[10]改进的剪切均值滤波，取消了对[αN]的判断求整，直接给出滤波的窗宽、单边剪切窗宽:

其中：Y是经过剪切均值滤波后一个窗宽的数据；N是滤波的窗宽(奇数)；Lw是单边剪切窗宽，其中需要满足2Lw

1) 从原始数据中选取对应滤波窗宽(N)长度的数据；

2) 对窗口内数据进行升降排序；

3) 对排序数列的两端数据进行剪切，单边剪切长度为Lw；

4) 将剩余数据取均值，作为最后剪切均值滤波的结果；

5) 滑动窗口对所有待处理数据进行剪切均值滤波。

以上给出了剪切均值滤波的实现过程。当对测线数据进行滤波处理时，需要考虑边界数据的处理方法。本文考虑几种常用的处理方法，如周期扩边(图2a)、镜像扩边(图2b)，两种方法均可以较好地对边界数据进行处理，保持数据的边界变化趋势。

图1 剪切均值滤波过程示意Fig.1 The diagram of Trim-mean filter processing

图2 两种扩边方法Fig.2 Two extension methods

2 天电去除研究

剪切均值滤波涉及到两个参数：滤波的窗宽、单边剪切窗宽，两个参数的取值直接影响滤波效果。下面对两个参数的选取进行研究，文中实测数据是由基频25 Hz的时间域航空电磁系统采集，采样率为100 kHz，半正弦发射波形的占空比为1∶4。

2.1 滤波窗宽

当滤波窗口长度N取值过小时，会出现某一窗口内均为天电噪声数据，导致处理结果中天电噪声去除不彻底；当滤波窗口长度N取值过大时，不仅会增加计算时间，且处理结果中供电和断电区间内的数据过于平均，导致供电窗口数据提前和断电窗口数据延后的特征，甚至将断电处有用信号滤平。现截取部分含有天电噪声的数据，参数中单边剪切长度设置为8，改变滤波窗口大小来说明去噪效果。

图3中分别给出了不同滤波窗口长度(N=21、41、101)的去噪效果。当N值较小(N=21、41)时，去噪结果虽然保持待处理数据的变化趋势(图3中b、c处)，但天电噪声滤除不彻底(图3中a处蓝色实线)；当滤波窗口长度较大(N=101)时，可以较好地滤除天电噪声(图3中a处紫色实线)，但会使得供电处、断电处及峰值数据偏离原数据方向，出现数据过于平均而丢失有用信号(图3中c处紫色实线)或者供电数据提前(图3中b处紫色实线)、断电数据延后(图3中c处紫色实线)、峰值减小的情况。从本组处理结果可以看出小窗口滤波可以很好地保持供电处、峰值、断电处的数据变化情况，大窗口滤波可以较好地去除数据中的天电噪声。

图3 不同滤波窗口长度对噪声去除效果比较Fig.3 Comparison of noise removal effects of different filter window lengths

2.2 剪切长度

剪切长度的大小即为窗口两端去除数据量的大小，它决定了窗口内噪声数据去除的比例。剪切长度过小或过大，均可能导致窗口被平均化，或供电处、断电处滤波后数据偏移正确的方向。

图4中分别给出了不同单边剪切长度的去噪效果，其中滤波窗口N=101，单边剪切长度Lw=8、30、45。经过滤波后的供电处、峰值、断电处，Lw=8、30均出现了数据被严重平均化，导致了供电窗口提前、峰值变小、断电窗口延迟等问题(图4中b、c处深蓝和天蓝实线)，Lw=8时以上现象较严重。当Lw=45时，供电和断电处滤波效果相对较好，但峰值出现了过于平均化现象(图4中b处紫色实线)。图4中断电数据含有天电噪声，由于天电噪声位于断电数据中部且为大窗口剪切均值滤波，因此剪切长度为8和30时均可以较好地去除噪声，当剪切长度取值过大(Lw=45)时，滤波后数据波动明显(图4中a处紫色实线)。可以看出滤波窗口长度和剪切长度要相互参考适当取值。

2.3 混合窗口天电噪声滤波

天电噪声是不可预见的，可考虑利用阈值法先识别天电噪声再进行去除，具体可参考文献[8]、[9]和[16]。根据前面讨论的情况可知，剪切均值滤波法的小窗口滤波适合数据趋势变化明显部分，如供电处、峰值和断电处，大窗口滤波适合断电后数据的处理，因此本文提出采用混合窗口滤波来对天电噪声进行去除，即供电处、峰值和断电处数据采用小窗口剪切均值滤波，断电后衰减数据使用大窗口剪切均值滤波。

图4 不同单边剪切长度对噪声去除效果比较Fig.4 Comparison of noise removal effects of different trim lengths

图5给出了利用混合窗口对数据进行天电去除的效果，其中小窗口滤波参数为：N=11、Lw=3，大窗口滤波参数为N=101、Lw=30。从图中可以明显看出，小窗口滤波后数据保持了原数据供电处、峰值和断电处的变化趋势(图5中a处)；大窗口滤波可以较好的滤除天电噪声(图5中b处供电阶段、c处断电阶段)。采用混合窗口剪切均值滤波即可以保持数据变化趋势，又可以达到有效去除天电噪声的目的。

图5 混合窗口剪切均值滤波去噪效果Fig.5 Effects of hybrid window trim-mean filter

3 结论

本文利用剪切均值滤波对时间域航空电磁数据进行天电噪声去除。航空电磁测线数据中包含了快速变化的供电数据及断电后变化相对缓慢的数据，在研究过程中发现，如对测线数据进行同时间道剪切均值滤波，则需要先去除运动噪声，否则会出现数据整体趋势的均值化。在对测线数据进行剪切均值滤波时，滤波效果受滤波参数影响。经过研究，小窗口参数滤波会造成数据中天电噪声去除不彻底，大窗口参数滤波会造成数据过度平均化，损失原始信号严重。为此，本文利用混合窗口剪切均值滤波对天电噪声进行有效去除。混合窗口方法的使用不仅保留了供电处、峰值和断电处有用信号的变化趋势和稳定性，且有效去除了天电噪声，是一种快速、高效去除天电噪声的方法。天电噪声对水平分量影响大，因此希望通过本文研究，不仅可以有效地去除天电噪声提高数据信噪比和增大勘探深度，且可以实现数据水平和垂直分量的联合应用。

致谢：向中科院上海微系统所助理研究员裴易峰和吉林大学电磁“千人计划”研究团队成员提供的帮助表示感谢。