李卫平, 何静, 曾祥平, 封隆永

(重庆市防雷中心,重庆 401147)

0 引言

大气电场是大气电学的一个重要参数。观测和研究表明,在雷暴等天气时大气电场场强受到扰动而剧烈变化,呈现出不同于晴天大气电场的显著特征[1-2]。雷雨云内电荷积累、电偶极子形态变化等均会导致地面大气电场场强发生变化,雷电活动时大量电荷在云内、云—空、云—地间转移,也会导致地面大气电场发生突变甚至正负反转。因此,可以通过地面大气电场变化情况来反演雷雨云中电场变化,从而对目标区域的雷电活动进行分析。雷电发生时,由于云内电荷积累变化或云闪、地闪,以及其他电磁干扰等造成地面大气电场场强变化情况非常复杂,单纯地从时域或频域进行分析不能完全掌握其全部情况,为了进一步了解雷暴过程中大气电场变化情况,本文基于重庆市南川区金佛山野外雷电观测场2021年3—12月大气电场数据,首先利用Savitzky-Golay卷积平滑滤波(以下简称S-G滤波)分离大气电场信号高频、低频部分,而后使用快速傅里叶变换(FFT)等信号处理方法对滤波前后信号及残差信号从时域、频域两方面进行分析,并结合当年ADTD闪电定位资料进行对比研究。

1 资料和方法

1.1 大气电场资料

大气电场资料来源于重庆市南川区金佛山野外雷电观测场的Pre-storm 2.0型场磨式大气电场仪于2021年3—12月不间断采集的数据。该电场仪工作基本原理:电机带动屏蔽金属片旋转,使感应金属片交替暴露在电场中或被屏蔽,从而产生与外界电场强度成正比的感应电荷,感应片连接到放大处理电路及波形调整电路,输出电压信号,经过标定,该电压可以表征大气电场的强度以及极性变化。其主要技术参数如表1所示。

表1 大气电场仪主要技术参数

1.2 闪电定位数据

本文采用重庆市闪电定位系统2021年地闪定位资料。闪电定位系统由第二代ADTD闪电定位仪、中心数据处理站、用户数据服务网络及图形显示终端组成,其观测网络由重庆市内的5个监测站和邻近省份的10个监测站组成。系统实现了对地闪时间、位置(经度、纬度)、雷电流峰值和极性不间断自动监测,探测效率85%以上,网内探测定位精度小于300 m,闪电回击的处理时间在1 ms左右。

1.3 数字滤波方法

S-G滤波器又称S-G卷积平滑器,它是一种特殊的低通滤波器,用来平滑噪声数据。该方法最初由SAVITZKY A和GOLAY M于1964年提出,而后被广泛地运用于信号去噪,采用在时域内基于多项式最小二乘法及窗口移动实现最佳拟合的方法。与通常的滤波器要经过时域-频域-时域变换不同,S-G滤波直接处理时域数据进行平滑,其平滑效果随窗口宽度不同而不同。相对于均值平滑滤波,S-G滤波更能保留相对极大值、极小值和宽度等分布特征[3-4]。该滤波算法的另一优点是其运算量相对较小,对计算机的内存及数据处理能力要求较低[5]。设滤波窗口的宽度n=2m+1,原始信号长为N,窗口内待平滑数据xi=(x-m,x-m+1,…,x0,x1,…,xm-1,xm),采用k-1次多项式对窗口内的数据进行多项式拟合:

yi=a0+a1x+a2x2+…+ak-1xk-1

这样的n个方程构成k元线性方程组(n>k),利用最小二乘法确定拟合参数a0,a1,a2,…,ak-1。即

其矩阵表达形式:

Y(2m+1)×1=X(2m+1)×k·Ak×1+E(2m+1)×1

2 分析和检验

2.1 窗口宽度及k值确定

利用S-G卷积平滑方法对大气电场信号进行滤波时,首先要选取合适的滤波窗口宽度及k值。对一段长度为N的大气电场离散序列而言,滤波窗口宽度n及拟合阶次k值大小对滤波效果影响显著,n越大,平滑效果越好;k越大,拟合效果越好,保留的细节信息也越多。但过大的n会导致信号曲线过度平滑,对高频部分去除较多,会导致部分信号的丢失;过大的k值,能保留信号的部分高频特征但平滑效果又较差,同时也大幅增加了计算的复杂度。对窗口宽度与拟合阶次的选取是对保持更多信号有效细节与平滑效果的综合考虑,为了兼顾平滑效果及尽量保留信号特征,随机选取一段长度为N的雷暴天气地面大气电场信号,选取不同的n、k值进行S-G滤波,其效果如图1所示。可以看出,原始大气电场信号由于含有较多的突变尖峰,信号重叠较为严重,很难看出电场变化的主要趋势,经S-G滤波后的低频信号则能够较好地反映雷暴天气时地面大气电场变化的主要趋势,在n=51,k=3时平滑效果最好,但丢失了信号的一些细节特征。在n=33,k=3或n=51,k=5时,滤波后的信号可以较好地反映地面大气电场变化的主要趋势,同时也能较好地保留信号突变等细节特征,考虑到计算的复杂度,选取n=33,k=3。

图1 不同窗口宽度、拟合阶次时滤波效果对比

2.2 大气电场信号频谱特征分析

基于小波分析、傅里叶变换、EMD分解等信号处理方法对雷暴天气大气电场数据进行频谱分析方面的研究较多,主要集中在去噪效果研究、雷暴天气与非雷暴天气大气电场信号频谱对比、雷电预警等[6-9]。本文不再对比分析雷暴与非雷暴天气时大气电场信号频谱区别,仅对雷暴天气时大气电场信号滤波前后及残差信号作时域、频域上的对比分析。对上述电场信号进行FFT变换并绘制频谱,如图2所示。由图2可以看出,雷暴天气电场离散序列信号经FFT变换后,其低频部分及直流部分幅值较大,高频部分幅值较小,且从0.05 Hz后幅值变化较为平稳,表明了电场能量主要存在于低频部分。随机选取20个雷暴天气过程大气电场信号进行FFT变换,取其相应频率的均值后得到平均频谱如图3所示,亦可得出上述结论。此特征与陈红兵等[10]的研究结论一致。设大气电场原始信号序列为S0(N),则有下式:

图2 滤波前后及残差信号频谱

图3 雷暴天气时大气电场平均频谱

S0(N)=Sl(N)+Sh(N)

定义Sh(N)为残差,在滤波窗口n=33,拟合次数k=3时,对滤波后的信号及残差进行FFT变换,可以看出大气电场信号在经过S-G滤波以后,信号的高、低频部分实现了较好的分离,滤波后的低频部分基本上在0.05 Hz以下,反映了电场信号的变化趋势。高频部分主要在0.04 Hz以上,反映了电场信号的快速变化特征。选取0.05 Hz作为高、低频之间的分界线,则原始信号S0、滤波后信号Sl、残差Sh的时域曲线如图4所示。由图4可见,残差Sh仅包含了微小的锯齿状波动及尖峰脉冲。

图4 滤波前后及残差信号时域波形

通常闪电的持续时间为几百毫秒,对大量地闪的观测表明,闪电每一次放电过程都会引起大气电场的突变,除先导和回击外,还有一系列更为细致的放电过程,也会引起电场变化[11]。由于大气电场仪采样频率为1 Hz,根据奈奎斯特采样定律,由样值序列无失真恢复原信号的条件是采样频率fs≥2fh,其中fh为信号的最高频率。可知,该大气电场仪离散数据序列能反映的电场信号的最高频率fh为0.5 Hz,高于0.5 Hz的信号在经FFT变换后的频域不能得到体现。由此推论,经滤波后得到的残差信号Sh(N)不能简单地作为噪声信号予以去除,相反,由于Sh(N)主要为信号的高频部分,包含了较为丰富的闪电相关信息,因此,有必要对残差信号作进一步的分析。

为直观展现闪电与残差信号之间的关系,将2021年5月3日03—04时,距大气电场仪20 km范围内的ADTD闪电定位信息与大气电场残差信号进行叠加,如图5所示。其中次纵轴为地闪定位点到电场仪的距离,可以很明显地看出,21个地闪点中有17个在时间上与残差信号曲线的尖峰脉冲相重合,占比为81%,说明距离较近的地闪是引起大气电场突变的重要因素,在频域中体现为信号的高频分量。从这方面来考虑,滤波后的残差信号包含了闪电信息,不能被作为噪声去除。但也可看出,残差信号尖峰脉冲的数量要多于地闪点数量,这可能由云闪、降水、云内电荷转移、电偶结构变化、电磁干扰等多方面因素造成,下一步还需要作更深入的研究。

图5 残差与地闪点叠加对比

2.3 进一步验证

雷暴天气时,大气电场仪记录的电场曲线波形不尽相同,为进一步验证滤波后残差与地闪之间的定性关系,选取闪电定位数据相对较多的3个时段的大气电场信号,即如图6所示的3种不同形状的信号波形进行S-G滤波,得到相应的残差信号波形,并与同时段20 km范围内的ADTD定位数据进行对比,查找时间上相重合的数量m并计算占定位数据总量的比例γ,得到表2。可以看出,地闪定位数据与残差信号脉冲时间重合比率在78%以上,进一步验证了残差信号包含闪电信息这一结论。

图6 3种不同形状的波形对比

表2 3种波形残差脉冲与地闪点重合情况

3 总结

1) 使用S-G滤波器对雷暴天气时大气电场信号进行滤波,选取窗口宽度n=33,拟合阶次k=3时,即可得到较为平滑的滤波曲线,同时也保留了必要的细节特征,相对于小波去噪等方法,其算法较为简便,如将多项式系数表预置在大气电场仪探头内,可实现探头滤波分离高、低频信号,为后续的雷暴预警及雷电识别等工作提供依据。

2) 对滤波前后的大气电场信号及残差进行FFT变换,则0.05 Hz为其高、低频的分界线,直流及低频部分的幅值较大,表明大气电场信号能量主要存在于0～0.05 Hz的直流及低频部分。大于0.05 Hz的高频部分幅值分布则较为平均。

3) 地闪是引起大气电场突变的重要因素,残差信号不能被简单地作为噪声予以去除,选取不同波形的大气电场数据进行滤波得到的残差信号与ADTD定位数据作对比,验证了其中包含较为丰富的闪电信息,但引起残差信号出现尖峰脉冲的因素很多,还需要作进一步的观测和研究。