基于Savitzky-Golay滤波的雷暴天气大气电场信号分析.pdf

Microcomputer Applications Vol.39,No.10,2023文章编号：10 0 7-7 57 X（2 0 2 3)10-0 0 14-0 4基于Savitzky-Golay滤波的雷暴天气大气电场信号分析基金项目微型电脑应用2 0 2 3年第39 卷第10 期李卫平，何静，曾祥平，封隆永（重庆市防雷中心，重庆40 1147）摘要：基于重庆市南川区金佛山野外雷电观测场2 0 2 1年大气电场数据及重庆市ADTD闪电定位资料，利用Savitzky-Golay卷积平滑滤波、快速傅里叶变换等信号处理方法对大气电场信号从时域、频域两方面进行分析。分析结果表明：选取窗口宽度n=33，拟合阶次k=3时，滤波效果较为理想；大气电场信号能量主要存在于0 0.0 5Hz的低频部分；地闪是引起大气电场突变的重要因素，残差信号不能简单地作为噪声予以去除，选取雷暴天气时不同波形的大气电场数据进行滤波得到的残差信号与ADTD定位数据作对比，验证了其中包含较为丰富的闪电信息。关键词：大气电场；S-G滤波；频谱分析；残差中图分类号：TP23文献标志码：AAnalysis of Atmospheric Electric Field Signal in Thunderstorm WeatherBased on Savitzky-Golay FilteringLI Weiping,HE Jing,ZENG Xiangping,FENG Longyong(Chongqing Lightning Protection Center,Chongqing 401147,China)Abstract:Based on the 2021 atmospheric electric field data of Jinfoshan lightning observation field in Nanchuan District and theADTD data of Chongqing,we analyze the atmospheric electric field signal in time domain and frequency domain by signal pro-cessing methods such as Savitzky-Golay filtering and fast Fourierian transform.The results show that when the window widthn=33 and the fitting order k=3,The filtering effect is relative ideal.The signal energy of atmospheric electric field mainly ex-ists in the low frequency part of 00.05 Hz;CG lightning is an important factor causing sudden change of atmospheric electricfield,and the residual signal can not be simply removed as noise.Compared the residual signal obtained by filtering the atmos-pheric electric field data of different waveforms in thunderstorm weather and the ADTD data,it is proved that the signals con-tain rich information of CG lightning.Key words:atmospheric electric field;S-G filtering;spectrum analysis;residuals卷积平滑滤波（以下简称S-G滤波）分离大气电场信号高频、0引言低频部分，而后使用快速傅里叶变换（FFT)等信号处理方法大气电场是大气电学的一个重要参数。观测和研究表对滤波前后信号及残差信号从时域、频域两方面进行分析，明，在雷暴等天气时大气电场场强受到扰动而剧烈变化，呈并结合当年ADTD闪电定位资料进行对比研究。现出不同于晴天大气电场的显著特征-2。雷雨云内电荷积1资料和方法累、电偶极子形态变化等均会导致地面大气电场场强发生变化，雷电活动时大量电荷在云内、云一空、云一地间转移，也会导致地面大气电场发生突变甚至正负反转。因此，可以通过地面大气电场变化情况来反演雷雨云中电场变化，从而对目标区域的雷电活动进行分析。雷电发生时，由于云内电荷积累变化或云闪、地闪，以及其他电磁干扰等造成地面大气电场场强变化情况非常复杂，单纯地从时域或频域进行分析不能完全掌握其全部情况，为了进一步了解雷暴过程中大气电场变化情况，本文基于重庆市南川区金佛山野外雷电观测场2 0 2 1年312 月大气电场数据，首先利用Savitzky-Golay作者简介：李卫平（198 1一），男，硕士，高级工程师，研究方向为雷电科学与防护技术；曾祥平（197 8 一），男，本科，高级工程师，研究方向为雷电科学与防护技术；封隆永（198 7 一），女，本科，工程师，研究方向为雷电科学与防护技术。通信作者：何静（198 6 一），女，本科，工程师，研究方向为雷电科学与防护技术。141.1大气电场资料大气电场资料来源于重庆市南川区金佛山野外雷电观测场的Pre-storm2.0型场磨式大气电场仪于2 0 2 1年3一12月不间断采集的数据。该电场仪工作基本原理：电机带动屏蔽金属片旋转，使感应金属片交替暴露在电场中或被屏蔽，从而产生与外界电场强度成正比的感应电荷，感应片连接到放大处理电路及波形调整电路，输出电压信号，经过标定，该电压可以表征大气电场的强度以及极性变化。其主要技术参数如表1所示。Microcomputer Applications Vol.39,No.10,20233W采样速率5ms测量范围300kV/m灵敏度5V/m分辨率1 V/m有效探测范围1020km测量误差5%转速误差1%1.2闪电定位数据本文采用重庆市闪电定位系统2 0 2 1年地闪定位资料。闪电定位系统由第二代ADTD闪电定位仪、中心数据处理站、用户数据服务网络及图形显示终端组成，其观测网络由重庆市内的5个监测站和邻近省份的10 个监测站组成。系统实现了对地闪时间、位置（经度、纬度）、雷电流峰值和极性不间断自动监测，探测效率8 5%以上，网内探测定位精度小：(0）1(c）o1：m-11(m-1)1ymL1(m）1其矩阵表达形式：Y(2m+1)x1=X(2m+1)X AkxI+E(2m+1)XIA的最小二乘解A为A=(XT,X)-1.XT.Y经滤波后，输出值或预测值Y为Y=X.A=X.(XT.X)-I.XT.Y=B.Y2分析和检验2.1窗口宽度及k值确定利用S-G卷积平滑方法对大气电场信号进行滤波时，首先要选取合适的滤波窗口宽度及k值。对一段长度为N的大气电场离散序列而言，滤波窗口宽度n及拟合阶次k值大小对滤波效果影响显著，n越大，平滑效果越好；k越大，拟合效果越好，保留的细节信息也越多。但过大的n会导致信号曲线过度平滑，对高频部分去除较多，会导致部分信号的丢失；过大的k值，能保留信号的部分高频特征但平滑效果又较差，同时也大幅增加了计算的复杂度。对窗口宽度与拟合阶次的选取是对保持更多信号有效细节与平滑效果的综合考虑，为了兼顾平滑效果及尽量保留信号特征，随机选取一段长度为N的雷暴天气地面大气电场信号，选取不同的n、k值进行S-G滤波,其效果如图1所示。可以看出，原始大气电场信号由于含有较多的突变尖峰，信号重叠较为严重，很基金项目表1大气电场仪主要技术参数于30 0 m，闪电回击的处理时间在1ms左右。名称主要技术参数电压范围24 V(DC)功耗(-m)1mJm+1：微型电脑应用2 0 2 3年第39 卷第10 期1.3数字滤波方法S-G滤波器又称 S-G卷积平滑器，它是一种特殊的低通滤波器,用来平滑噪声数据。该方法最初由 SAVITZKY A和GOLAYM于19 6 4年提出，而后被广泛地运用于信号去噪，采用在时域内基于多项式最小二乘法及窗口移动实现最佳拟合的方法。与通常的滤波器要经过时域一频域一时域变换不同，S-G滤波直接处理时域数据进行平滑，其平滑效果随窗口宽度不同而不同。相对于均值平滑滤波，S-G滤波更能保留相对极大值、极小值和宽度等分布特征3-41。该滤波算法的另一优点是其运算量相对较小，对计算机的内存及数据处理能力要求较低5。设滤波窗口的宽度n=2m十1，原始信号长为N，窗口内待平滑数据a,=（-m，-m+1，ao，1，,c m-1，a m）,采用k一1次多项式对窗口内的数据进行多项式拟合：yi=ao+ai+a2a+ak-1a-1这样的n个方程构成k元线性方程组（nk），利用最小二乘法确定拟合参数aoa1,a2,ak-1。即(-m）(-m）k-11(-m+1)1：(m）(cm)k-1难看出电场变化的主要趋势，经S-G滤波后的低频信号则能够较好地反映雷暴天气时地面大气电场变化的主要趋势，在n=51,k=3时平滑效果最好，但丢失了信号的一些细节特征。在n=33,k=3或n=51,k=5时，滤波后的信号可以较好地反映地面大气电场变化的主要趋势，同时也能较好地保留信号突变等细节特征，考虑到计算的复杂度，选取n=33，k=3。2.2大气电场信号频谱特征分析基于小波分析、傅里叶变换、EMD分解等信号处理方法对雷暴天气大气电场数据进行频谱分析方面的研究较多，主要集中在去噪效果研究、雷暴天气与非雷暴天气大气电场信号频谱对比、雷电预警等6-9。本文不再对比分析雷暴与非雷暴天气时大气电场信号频谱区别，仅对雷暴天气时大气电场信号滤波前后及残差信号作时域、频域上的对比分析。对上述电场信号进行FFT变换并绘制频谱，如图2 所示。由图2 可以看出，雷暴天气电场离散序列信号经FFT变换后，其低频部分及直流部分幅值较大，高频部分幅值较小，且从0.05Hz后幅值变化较为平稳，表明了电场能量主要存在于低频部分。随机选取2 0 个雷暴天气过程大气电场信号进行FFT变换，取其相应频率的均值后得到平均频谱如图3所示，亦可得出上述结论。此特征与陈红兵等10 的研究结论一致。设大气电场原始信号序列为S。（N），则有下式：S。(N)=S(N)+Sh(N)15.ao(a-m+1)(a-m+1)k-1：(co)-1：(acm-1)-1厂emale-m+1：ak-2+0：ak-2emr-1Lak-1emMicrocomputer Applications Vol.39,No.10,2023基金项目微型电脑应用2 0 2 3年第39 卷第10 期20100-10078F-10020100-10078F-1001.0A/0.50A/1.00.501.00.50图2滤波前后及残差信号频谱0.8r0.70.60.5A/0.40.30.20.10图3雷暴天气时大气电场平均频谱定义S.（N)为残差，在滤波窗口n=33，拟合次数k=3时，对滤波后的信号及残差进行FFT变换，可以看出大气电场信号在经过S-G滤波以后，信号的高、低频部分实现了较好的分离，滤波后的低频部分基本上在0.0 5Hz以下，反映了电场信号的变化趋势。高频部分主要在0.0 4Hz以上，反映了电场信号的快速变化特征。选取0.0 5Hz作为高、低频之间的分界线，则原始信号S。、滤波后信号S、残差S的时域曲线如图4所示。由图4可见，残差S仅包含了微小的锯齿状波动及尖峰脉冲。200100-10500100015002000250030005001000 15002000250030005001000150020002500300050010001500200025003000图1不同窗口宽度、拟合阶次时滤波效果对比0.040.08f/Hz0.040.08f/Hz0.040.08f/Hz0.040.08f/Hz078F0-1002010-10078F-1000.120.160.120.160.120.160.120.18500500100015002000 2500 300050010001500 2000 250030005001000 15002000 25003000200.20-200500100015002000 2500300020-200.200500100015002000 25003000500-50-1000.20050010001500 200025003000图4滤波前后及残差信号时域波形通常闪电的持续时间为几百毫秒，对大量地闪的观测表明，闪电每一次放电过程都会引起大气电场的突变，除先导和回击外，还有一系列更为细致的放电过程，也会引起电场变化1。由于大气电场仪采样频率为1Hz，根据奈奎斯特采样定律，由样值序列无失真恢复原信号的条件是采样频率f2fh,其中f为信号的最高频率。可知，该大气电场仪离散数据序列能反映的电场信号的最高频率f为0.5Hz，高于0.5Hz的信号在经FFT变换后的频域不能得到体现。由此推论，经滤波后得到的残差信号S,（N)不能简单地作为0.20噪声信号予以去除，相反，由于S,（N)主要为信号的高频部分，包含了较为丰富的闪电相关信息，因此，有必要对残差信号作进一步的分析。为直观展现闪电与残差信号之间的关系，将2 0 2 1年5月3日0 3一0 4时，距大气电场仪2 0 km范围内的ADTD闪电定位信息与大气电场残差信号进行叠加，如图5所示。其中次纵轴为地闪定位点到电场仪的距离，可以很明显地看出，2 1个地闪点中有17 个在时间上与残差信号曲线的尖峰脉冲相重合，占比为8 1%，说明距离较近的地闪是引起大气电场突变的重要因素，在频域中体现为信号的高频分量。从这方面来考虑，滤波后的残差信号包含了闪电信息，不能被1610001500200025003000Microcomputer Applications Vol.39,No.10,2023作为噪声去除。但也可看出，残差信号尖峰脉冲的数量要多于地闪点数量，这可能由云闪、降水、云内电荷转移、电偶结构变化、电磁干扰等多方面因素造成，下一步还需要作更深人的研究。40/-20-30402.3进一步验证雷暴天气时，大气电场仪记录的电场曲线波形不尽相同，为进一步验证滤波后残差与地闪之间的定性关系，选取闪电定位数据相对较多的3个时段的大气电场信号，即如图6所示的3种不同形状的信号波形进行S-G滤波，得到相应的残差信号波形，并与同时段2 0 km范围内的ADTD定位数据进行对比，查找时间上相重合的数量m并计算占定位数据总量的比例，得到表2。可以看出，地闪定位数据与残差信号脉冲时间重合比率在7 8%以上，进一步验证了残差信号包含闪电信息这一结论。20A/-20020/-20020A/0-200图6 3种不同形状的波形对比表2 3种波形残差脉冲与地闪点重合情况波形时长/s220 km内地闪总量重合数量占比/%波形13000波形23000波形330003总结1）使用S-G滤波器对雷暴天气时大气电场信号进行滤基金项目波，选取窗口宽度n=33，拟合阶次k=3时，即可得到较为平滑的滤波曲线，同时也保留了必要的细节特征，相对于小波去噪等方法，其算法较为简便，如将多项式系数表预置在大气电场仪探头内，可实现探头滤波分离高、低频信号，为后续的雷暴预警及雷电识别等工作提供依据。一残差值725地闪点205时间/s图5残差与地闪点叠加对比500100015002200025003000时间/s50010001500200025003000时间/s50010001500200025003000时间/s777160微型电脑应用2 0 2 3年第39 卷第10 期2）对滤波前后的大气电场信号及残差进行FFT变换，则0.0 5Hz为其高、低频的分界线，直流及低频部分的幅值较大，表明大气电场信号能量主要存在于0 0.0 5Hz的直流及低频部分。大于0.0 5Hz的高频部分幅值分布则较为平均。3）地闪是引起大气电场突变的重要因素，残差信号不能被简单地作为噪声予以去除，选取不同波形的大气电场数据进行滤波得到的残差信号与ADTD定位数据作对比，验证了其中包含较为丰富的闪电信息，但引起残差信号出现尖峰脉冲的因素很多，还需要作进一步的观测和研究。1秦微，张其林，姜苏，等.基于大气电场资料的雷电临近预警研究J.南京信息工程大学学报（自然科学版),2 0 16,8(3):2 47-2 51.2 徐栋璞，王振会，曾庆锋，等.近地面大气电场数据EMD方法分析JI.气象科技，2 0 13,41(1)：17 0-17 6.3吴水龙，郭阳，曲广浩，等.基于Savitzky-Golay滤波器的多变量监控报警方法研究J.计算机与应用化学，2019,36(5):478-482.4胡顺石，黄春晓，杨斌，等。自适应加权Savitzky-Golay滤波重构MODIS植被指数时间序列JI.测绘科学，2020,45(4):105-116.5蔡天净，唐瀚.Savitzky-Golay平滑滤波器的最小二乘拟合原理综述J.数字通信，2 0 11，38（1)：6 3-6 8.6 李颖，王振会，肖稳安，等.大气电场的FFT频谱分析及雷暴预报研究J.气象科学，2 0 13，33（1）：6 6-7 0.7刘鹭，姚慧茹.雷电过程中大气电场信号处理方法探究J.沙漠与绿洲气象，2 0 2 1，15（5）：12 3-131.8 行鸿彦，张强，季鑫源.基于地面电场资料的雷暴临近预报研究J.大气科学学报，2 0 17，40(1)：111-117.9余海，张廷龙，刘文杰，等.基于小波分析对地面电场信号的去噪及闪电信息提取J.海南大学学报（自然科学版）2 0 2 0,38（1)：6 7-7 4.6685.715678.875490.00参考文献10陈红兵，徐文，王振会，等.基于小波分析的地面大气电场观测数据处理技术研究J.大气科学学报，2012,35(6)762-767.11陈渭民.雷电学原理M.2版.北京：气象出版社，2006.（收稿日期：2 0 2 2-0 3-18）17.