基于Savitzky-Golay濾波的雷暴天氣大氣電場(chǎng)信號(hào)分析

李衛(wèi)平, 何靜, 曾祥平, 封隆永

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

0 引言

大氣電場(chǎng)是大氣電學(xué)的一個(gè)重要參數(shù)。觀測(cè)和研究表明,在雷暴等天氣時(shí)大氣電場(chǎng)場(chǎng)強(qiáng)受到擾動(dòng)而劇烈變化,呈現(xiàn)出不同于晴天大氣電場(chǎng)的顯著特征[1-2]。雷雨云內(nèi)電荷積累、電偶極子形態(tài)變化等均會(huì)導(dǎo)致地面大氣電場(chǎng)場(chǎng)強(qiáng)發(fā)生變化,雷電活動(dòng)時(shí)大量電荷在云內(nèi)、云—空、云—地間轉(zhuǎn)移,也會(huì)導(dǎo)致地面大氣電場(chǎng)發(fā)生突變甚至正負(fù)反轉(zhuǎn)。因此,可以通過(guò)地面大氣電場(chǎng)變化情況來(lái)反演雷雨云中電場(chǎng)變化,從而對(duì)目標(biāo)區(qū)域的雷電活動(dòng)進(jìn)行分析。雷電發(fā)生時(shí),由于云內(nèi)電荷積累變化或云閃、地閃,以及其他電磁干擾等造成地面大氣電場(chǎng)場(chǎng)強(qiáng)變化情況非常復(fù)雜,單純地從時(shí)域或頻域進(jìn)行分析不能完全掌握其全部情況,為了進(jìn)一步了解雷暴過(guò)程中大氣電場(chǎng)變化情況,本文基于重慶市南川區(qū)金佛山野外雷電觀測(cè)場(chǎng)2021年3—12月大氣電場(chǎng)數(shù)據(jù),首先利用Savitzky-Golay卷積平滑濾波(以下簡(jiǎn)稱S-G濾波)分離大氣電場(chǎng)信號(hào)高頻、低頻部分,而后使用快速傅里葉變換(FFT)等信號(hào)處理方法對(duì)濾波前后信號(hào)及殘差信號(hào)從時(shí)域、頻域兩方面進(jìn)行分析,并結(jié)合當(dāng)年ADTD閃電定位資料進(jìn)行對(duì)比研究。

1 資料和方法

1.1 大氣電場(chǎng)資料

大氣電場(chǎng)資料來(lái)源于重慶市南川區(qū)金佛山野外雷電觀測(cè)場(chǎng)的Pre-storm 2.0型場(chǎng)磨式大氣電場(chǎng)儀于2021年3—12月不間斷采集的數(shù)據(jù)。該電場(chǎng)儀工作基本原理:電機(jī)帶動(dòng)屏蔽金屬片旋轉(zhuǎn),使感應(yīng)金屬片交替暴露在電場(chǎng)中或被屏蔽,從而產(chǎn)生與外界電場(chǎng)強(qiáng)度成正比的感應(yīng)電荷,感應(yīng)片連接到放大處理電路及波形調(diào)整電路,輸出電壓信號(hào),經(jīng)過(guò)標(biāo)定,該電壓可以表征大氣電場(chǎng)的強(qiáng)度以及極性變化。其主要技術(shù)參數(shù)如表1所示。

表1 大氣電場(chǎng)儀主要技術(shù)參數(shù)

1.2 閃電定位數(shù)據(jù)

本文采用重慶市閃電定位系統(tǒng)2021年地閃定位資料。閃電定位系統(tǒng)由第二代ADTD閃電定位儀、中心數(shù)據(jù)處理站、用戶數(shù)據(jù)服務(wù)網(wǎng)絡(luò)及圖形顯示終端組成,其觀測(cè)網(wǎng)絡(luò)由重慶市內(nèi)的5個(gè)監(jiān)測(cè)站和鄰近省份的10個(gè)監(jiān)測(cè)站組成。系統(tǒng)實(shí)現(xiàn)了對(duì)地閃時(shí)間、位置(經(jīng)度、緯度)、雷電流峰值和極性不間斷自動(dòng)監(jiān)測(cè),探測(cè)效率85%以上,網(wǎng)內(nèi)探測(cè)定位精度小于300 m,閃電回?fù)舻奶幚頃r(shí)間在1 ms左右。

1.3 數(shù)字濾波方法

S-G濾波器又稱S-G卷積平滑器,它是一種特殊的低通濾波器,用來(lái)平滑噪聲數(shù)據(jù)。該方法最初由SAVITZKY A和GOLAY M于1964年提出,而后被廣泛地運(yùn)用于信號(hào)去噪,采用在時(shí)域內(nèi)基于多項(xiàng)式最小二乘法及窗口移動(dòng)實(shí)現(xiàn)最佳擬合的方法。與通常的濾波器要經(jīng)過(guò)時(shí)域-頻域-時(shí)域變換不同,S-G濾波直接處理時(shí)域數(shù)據(jù)進(jìn)行平滑,其平滑效果隨窗口寬度不同而不同。相對(duì)于均值平滑濾波,S-G濾波更能保留相對(duì)極大值、極小值和寬度等分布特征[3-4]。該濾波算法的另一優(yōu)點(diǎn)是其運(yùn)算量相對(duì)較小,對(duì)計(jì)算機(jī)的內(nèi)存及數(shù)據(jù)處理能力要求較低[5]。設(shè)濾波窗口的寬度n=2m+1,原始信號(hào)長(zhǎng)為N,窗口內(nèi)待平滑數(shù)據(jù)xi=(x-m,x-m+1,…,x0,x1,…,xm-1,xm),采用k-1次多項(xiàng)式對(duì)窗口內(nèi)的數(shù)據(jù)進(jìn)行多項(xiàng)式擬合:

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

這樣的n個(gè)方程構(gòu)成k元線性方程組(n>k),利用最小二乘法確定擬合參數(shù)a0,a1,a2,…,ak-1。即

其矩陣表達(dá)形式:

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

2 分析和檢驗(yàn)

2.1 窗口寬度及k值確定

利用S-G卷積平滑方法對(duì)大氣電場(chǎng)信號(hào)進(jìn)行濾波時(shí),首先要選取合適的濾波窗口寬度及k值。對(duì)一段長(zhǎng)度為N的大氣電場(chǎng)離散序列而言,濾波窗口寬度n及擬合階次k值大小對(duì)濾波效果影響顯著,n越大,平滑效果越好;k越大,擬合效果越好,保留的細(xì)節(jié)信息也越多。但過(guò)大的n會(huì)導(dǎo)致信號(hào)曲線過(guò)度平滑,對(duì)高頻部分去除較多,會(huì)導(dǎo)致部分信號(hào)的丟失;過(guò)大的k值,能保留信號(hào)的部分高頻特征但平滑效果又較差,同時(shí)也大幅增加了計(jì)算的復(fù)雜度。對(duì)窗口寬度與擬合階次的選取是對(duì)保持更多信號(hào)有效細(xì)節(jié)與平滑效果的綜合考慮,為了兼顧平滑效果及盡量保留信號(hào)特征,隨機(jī)選取一段長(zhǎng)度為N的雷暴天氣地面大氣電場(chǎng)信號(hào),選取不同的n、k值進(jìn)行S-G濾波,其效果如圖1所示?？梢钥闯?原始大氣電場(chǎng)信號(hào)由于含有較多的突變尖峰,信號(hào)重疊較為嚴(yán)重,很難看出電場(chǎng)變化的主要趨勢(shì),經(jīng)S-G濾波后的低頻信號(hào)則能夠較好地反映雷暴天氣時(shí)地面大氣電場(chǎng)變化的主要趨勢(shì),在n=51,k=3時(shí)平滑效果最好,但丟失了信號(hào)的一些細(xì)節(jié)特征。在n=33,k=3或n=51,k=5時(shí),濾波后的信號(hào)可以較好地反映地面大氣電場(chǎng)變化的主要趨勢(shì),同時(shí)也能較好地保留信號(hào)突變等細(xì)節(jié)特征,考慮到計(jì)算的復(fù)雜度,選取n=33,k=3。

圖1 不同窗口寬度、擬合階次時(shí)濾波效果對(duì)比

2.2 大氣電場(chǎng)信號(hào)頻譜特征分析

基于小波分析、傅里葉變換、EMD分解等信號(hào)處理方法對(duì)雷暴天氣大氣電場(chǎng)數(shù)據(jù)進(jìn)行頻譜分析方面的研究較多,主要集中在去噪效果研究、雷暴天氣與非雷暴天氣大氣電場(chǎng)信號(hào)頻譜對(duì)比、雷電預(yù)警等[6-9]。本文不再對(duì)比分析雷暴與非雷暴天氣時(shí)大氣電場(chǎng)信號(hào)頻譜區(qū)別,僅對(duì)雷暴天氣時(shí)大氣電場(chǎng)信號(hào)濾波前后及殘差信號(hào)作時(shí)域、頻域上的對(duì)比分析。對(duì)上述電場(chǎng)信號(hào)進(jìn)行FFT變換并繪制頻譜,如圖2所示。由圖2可以看出,雷暴天氣電場(chǎng)離散序列信號(hào)經(jīng)FFT變換后,其低頻部分及直流部分幅值較大,高頻部分幅值較小,且從0.05 Hz后幅值變化較為平穩(wěn),表明了電場(chǎng)能量主要存在于低頻部分。隨機(jī)選取20個(gè)雷暴天氣過(guò)程大氣電場(chǎng)信號(hào)進(jìn)行FFT變換,取其相應(yīng)頻率的均值后得到平均頻譜如圖3所示,亦可得出上述結(jié)論。此特征與陳紅兵等[10]的研究結(jié)論一致。設(shè)大氣電場(chǎng)原始信號(hào)序列為S0(N),則有下式:

圖2 濾波前后及殘差信號(hào)頻譜

圖3 雷暴天氣時(shí)大氣電場(chǎng)平均頻譜

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

定義Sh(N)為殘差,在濾波窗口n=33,擬合次數(shù)k=3時(shí),對(duì)濾波后的信號(hào)及殘差進(jìn)行FFT變換,可以看出大氣電場(chǎng)信號(hào)在經(jīng)過(guò)S-G濾波以后,信號(hào)的高、低頻部分實(shí)現(xiàn)了較好的分離,濾波后的低頻部分基本上在0.05 Hz以下,反映了電場(chǎng)信號(hào)的變化趨勢(shì)。高頻部分主要在0.04 Hz以上,反映了電場(chǎng)信號(hào)的快速變化特征。選取0.05 Hz作為高、低頻之間的分界線,則原始信號(hào)S0、濾波后信號(hào)Sl、殘差Sh的時(shí)域曲線如圖4所示。由圖4可見(jiàn),殘差Sh僅包含了微小的鋸齒狀波動(dòng)及尖峰脈沖。

圖4 濾波前后及殘差信號(hào)時(shí)域波形

通常閃電的持續(xù)時(shí)間為幾百毫秒,對(duì)大量地閃的觀測(cè)表明,閃電每一次放電過(guò)程都會(huì)引起大氣電場(chǎng)的突變,除先導(dǎo)和回?fù)敉?還有一系列更為細(xì)致的放電過(guò)程,也會(huì)引起電場(chǎng)變化[11]。由于大氣電場(chǎng)儀采樣頻率為1 Hz,根據(jù)奈奎斯特采樣定律,由樣值序列無(wú)失真恢復(fù)原信號(hào)的條件是采樣頻率fs≥2fh,其中fh為信號(hào)的最高頻率?？芍?該大氣電場(chǎng)儀離散數(shù)據(jù)序列能反映的電場(chǎng)信號(hào)的最高頻率fh為0.5 Hz,高于0.5 Hz的信號(hào)在經(jīng)FFT變換后的頻域不能得到體現(xiàn)。由此推論,經(jīng)濾波后得到的殘差信號(hào)Sh(N)不能簡(jiǎn)單地作為噪聲信號(hào)予以去除,相反,由于Sh(N)主要為信號(hào)的高頻部分,包含了較為豐富的閃電相關(guān)信息,因此,有必要對(duì)殘差信號(hào)作進(jìn)一步的分析。

為直觀展現(xiàn)閃電與殘差信號(hào)之間的關(guān)系,將2021年5月3日03—04時(shí),距大氣電場(chǎng)儀20 km范圍內(nèi)的ADTD閃電定位信息與大氣電場(chǎng)殘差信號(hào)進(jìn)行疊加,如圖5所示。其中次縱軸為地閃定位點(diǎn)到電場(chǎng)儀的距離,可以很明顯地看出,21個(gè)地閃點(diǎn)中有17個(gè)在時(shí)間上與殘差信號(hào)曲線的尖峰脈沖相重合,占比為81%,說(shuō)明距離較近的地閃是引起大氣電場(chǎng)突變的重要因素,在頻域中體現(xiàn)為信號(hào)的高頻分量。從這方面來(lái)考慮,濾波后的殘差信號(hào)包含了閃電信息,不能被作為噪聲去除。但也可看出,殘差信號(hào)尖峰脈沖的數(shù)量要多于地閃點(diǎn)數(shù)量,這可能由云閃、降水、云內(nèi)電荷轉(zhuǎn)移、電偶結(jié)構(gòu)變化、電磁干擾等多方面因素造成,下一步還需要作更深入的研究。

圖5 殘差與地閃點(diǎn)疊加對(duì)比

2.3 進(jìn)一步驗(yàn)證

雷暴天氣時(shí),大氣電場(chǎng)儀記錄的電場(chǎng)曲線波形不盡相同,為進(jìn)一步驗(yàn)證濾波后殘差與地閃之間的定性關(guān)系,選取閃電定位數(shù)據(jù)相對(duì)較多的3個(gè)時(shí)段的大氣電場(chǎng)信號(hào),即如圖6所示的3種不同形狀的信號(hào)波形進(jìn)行S-G濾波,得到相應(yīng)的殘差信號(hào)波形,并與同時(shí)段20 km范圍內(nèi)的ADTD定位數(shù)據(jù)進(jìn)行對(duì)比,查找時(shí)間上相重合的數(shù)量m并計(jì)算占定位數(shù)據(jù)總量的比例γ,得到表2?？梢钥闯?地閃定位數(shù)據(jù)與殘差信號(hào)脈沖時(shí)間重合比率在78%以上,進(jìn)一步驗(yàn)證了殘差信號(hào)包含閃電信息這一結(jié)論。

圖6 3種不同形狀的波形對(duì)比

表2 3種波形殘差脈沖與地閃點(diǎn)重合情況

3 總結(jié)

1) 使用S-G濾波器對(duì)雷暴天氣時(shí)大氣電場(chǎng)信號(hào)進(jìn)行濾波,選取窗口寬度n=33,擬合階次k=3時(shí),即可得到較為平滑的濾波曲線,同時(shí)也保留了必要的細(xì)節(jié)特征,相對(duì)于小波去噪等方法,其算法較為簡(jiǎn)便,如將多項(xiàng)式系數(shù)表預(yù)置在大氣電場(chǎng)儀探頭內(nèi),可實(shí)現(xiàn)探頭濾波分離高、低頻信號(hào),為后續(xù)的雷暴預(yù)警及雷電識(shí)別等工作提供依據(jù)。

2) 對(duì)濾波前后的大氣電場(chǎng)信號(hào)及殘差進(jìn)行FFT變換,則0.05 Hz為其高、低頻的分界線,直流及低頻部分的幅值較大,表明大氣電場(chǎng)信號(hào)能量主要存在于0～0.05 Hz的直流及低頻部分。大于0.05 Hz的高頻部分幅值分布則較為平均。

3) 地閃是引起大氣電場(chǎng)突變的重要因素,殘差信號(hào)不能被簡(jiǎn)單地作為噪聲予以去除,選取不同波形的大氣電場(chǎng)數(shù)據(jù)進(jìn)行濾波得到的殘差信號(hào)與ADTD定位數(shù)據(jù)作對(duì)比,驗(yàn)證了其中包含較為豐富的閃電信息,但引起殘差信號(hào)出現(xiàn)尖峰脈沖的因素很多,還需要作進(jìn)一步的觀測(cè)和研究。