基于Savitzky-Golay濾波的雷暴天氣大氣電場信號分析

李衛平, 何靜, 曾祥平, 封隆永

(重慶市防雷中心,重慶 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定位數據作對比,驗證了其中包含較為豐富的閃電信息,但引起殘差信號出現尖峰脈沖的因素很多,還需要作進一步的觀測和研究。