袁靜,王子杰,澤仁志瑪,王志國,豐繼林,申旭輝,吳鵬,王橋,楊德賀,王統(tǒng)領(lǐng),周樂

1 防災科技學院，河北三河 065201 2 應急管理部國家自然災害防治研究院，北京 100085 3 清華大學，北京 100084 4 湖州師范學院體育學院，浙江湖州 313000

0 引言

閃電是高空頻發(fā)的自然災害，全球范圍平均每秒發(fā)生約44次，全年累計約14億次(Christian et al.,2003).閃電產(chǎn)生寬頻帶的電磁脈沖，能夠傳播到電離層并激發(fā)起電磁哨聲波.哨聲波是由于電磁波在傳播過程中高低頻成分之間存在相速度差,通常情況下高頻相速度快提前到達衛(wèi)星高度，低頻相速度慢后續(xù)到達，導致其在衛(wèi)星記錄的電磁場時頻圖中呈現(xiàn)頻率隨時間下降的“L”形態(tài)的色散狀(Barkhausen,1930;Storey,1953;Helliwell,1965).我國首顆電磁監(jiān)測試驗衛(wèi)星張衡一號衛(wèi)星記錄的典型閃電哨聲波如圖1所示.當閃電哨聲波傳播的路徑較長、電子密度較高或磁場強度較強時，色散變大(Carpenter and Anderson,1992).由于閃電哨聲波的形態(tài)攜帶了大量的空間環(huán)境信息，被廣泛用于空間環(huán)境監(jiān)測中，是揭開圈層耦合機理的重要研究手段(Chen et al.,2017;Carpenter and Anderson,1992;Singh et al.,2018;Oike et al.,2014;Bayupati et al.,2012;Clilverd et al.,2002;Kishore et al.,2016;Záhlava et al.,2018；Parrot et al.,2019；Horne et al.,2013;羅旭東等,2017).

圖1 張衡一號衛(wèi)星記錄的VLF頻段磁場X分量上的閃電哨聲波(Bx分量)Fig.1 The lightning whistler observed at VLF magnetic field′s X component by ZH-1 satellite (Bx component)

當前，基于閃電哨聲波的空間物理研究主要針對衛(wèi)星記錄的單個閃電哨聲波事件進行深入分析并反演其空間物理環(huán)境的相關(guān)參數(shù)，然而通常閃電哨聲波事件淹沒在衛(wèi)星觀測的海量電磁場數(shù)據(jù)中，完全依賴人工識別效率低下，難度極高，導致了閃電哨聲波全球時空分布規(guī)律和相關(guān)參數(shù)的研究甚少.

2008年，國內(nèi)外學者開始借助人工智能技術(shù)克服難點，初步發(fā)展了閃電哨聲波圖像自動識別算法.目前閃電哨聲波識別算法的流程是首先對觀測波形數(shù)據(jù)進行帶通濾波，然后利用快速傅里葉變換將波形數(shù)據(jù)轉(zhuǎn)化成時頻圖，最后借助機器學習或計算機視覺技術(shù)等自動識別時頻圖中的L色散形態(tài).比如Lichtenberger等(2008)在反演電子密度時提出依賴人工處理大量地面甚低頻(Very Low Frequency，VLF)觀測數(shù)據(jù)將導致相關(guān)研究面臨巨大的技術(shù)瓶頸，并提出了基于滑動模板匹配技術(shù)的閃電哨聲波自動檢測方法，其模板制作符合Bernard提出的閃電哨聲波的形態(tài)特征(Bernard,1973).該算法已經(jīng)大規(guī)模應用于Marion和SANAE的VLF地面觀測站數(shù)據(jù)處理，其存在的缺陷是需要事先從時頻圖中移除閃電脈沖、地面電力系統(tǒng)諧波輻射和人工發(fā)射源等引起的干擾現(xiàn)象,具有較高虛警率和漏檢率(Lichtenberger et al.,2008).Zhou等(2020)針對武漢VLF地面觀測數(shù)據(jù)中隱藏的閃電吱聲(tweek)現(xiàn)象，通過設置能量譜閾值和時間寬度閾值的方式提出了簡單快捷地自動識別算法.但該方法不適用于從張衡衛(wèi)星數(shù)據(jù)中識別閃電哨聲波，因為張衡一號衛(wèi)星電磁場觀測載荷具有高靈敏度的特點，導致時頻圖的背景噪聲的能量譜強度與閃電哨聲波軌跡的能量譜強度相差不大，難以通過設置能量譜閾值的方式對閃電哨聲波進行粗定位.斯坦福大學的VLF電波研究小組Stanford VLF Group(2009)最早對電磁衛(wèi)星的閃電哨聲波現(xiàn)象開展自動識別:截取固定時間窗的時頻圖，結(jié)合去噪處理、網(wǎng)格劃分、計算平均幅度值等計算機視覺技術(shù)完成特征提取功能，最后采用模板匹配的分類策略實現(xiàn)閃電哨聲波粗定位.Dharma等(2014)認為該特征的性能受制于網(wǎng)格劃分的數(shù)量，利用閃電哨聲波區(qū)域的顏色變化較小且具有明顯連通的特點，提出了基于連通域分析的閃電哨聲波粗定位方法.其缺陷是特征魯棒性不高，算法效果受背景噪聲影響很大.Oike等(2014)和Fiser等(2010)以Eckersley公式(Eckersley,1935)為基礎，借助觀測數(shù)據(jù)分別制作了白天和夜間的閃電哨聲波模板，再采用互相關(guān)熵的模板匹配策略完成閃電哨聲波識別和粗定位.Ali Ahmad等(2019)認為上述閃電哨聲波模板不符合實際情況，借助邊緣提取等計算機視覺技術(shù)提出了表達多種色散形態(tài)的特征提取方法，最后利用基于決策樹規(guī)則的分類算法完成識別.袁靜等(2021a)認為按照上述方法提取的閃電哨聲波特征魯棒性差且分類器的泛化性能低，則根據(jù)閃電哨聲波的色散形態(tài)設計了能夠增強閃電哨聲波L色散特征的卷積核并采用支持向量機做分類器獲得了較高的識別性能.

鑒于深度神經(jīng)網(wǎng)絡在提取圖像特征和擬合非線性函數(shù)等方面獲得的巨大突破(LeCun et al.,2015;Liu et al.,2018)，Konan等(2020)提出了兩種基于深度神經(jīng)網(wǎng)絡的閃電哨聲波粗定位算法：基于滑動深度卷積神經(jīng)網(wǎng)絡(Sliding Deep Neural Convolutional Neural Network,SDNN)算法和基于YOLOv3 (You Only Look Once version 3rd)神經(jīng)網(wǎng)絡算法.SDNN神經(jīng)網(wǎng)絡主要包含兩部分：3個卷積層和2個分類層.該算法的實施過程：截取某固定時間寬度的時頻圖，利用卷積層提取圖像特征，再利用分類層進行識別，從而實現(xiàn)閃電哨聲波粗定位(Konan et al.,2020).該算法具有高魯棒性特征的提取功能，泛化能力強的分類功能，但漏檢率高，原因是基于固定時間寬度的定位策略容易漏檢其他時間寬度的閃電哨聲波(Konan et al.,2020；袁靜等，2021b).YOLOv3神經(jīng)網(wǎng)絡的閃電哨聲波檢測算法包含兩個主要組成部分：YOLOv3主體網(wǎng)絡和非極大值抑制(Non-Maximum Suppression，NMS)算法.YOLOv3主體網(wǎng)絡主要由75個卷積層構(gòu)成，無全連接層，適應任意大小的輸入圖像；無池化層，尺度不變特征能傳送到下一層；采用殘差結(jié)構(gòu)極大地降低了學習魯棒特征的難度；利用上下文圖像信息構(gòu)造目標定位的數(shù)學模型；通過上述過程主體網(wǎng)絡將輸出多個粗定位的預測框，最后再采用NMS算法對主體網(wǎng)絡輸出的預測框進行過濾和優(yōu)化,從而實現(xiàn)閃電哨聲波區(qū)域粗定位；以上優(yōu)勢使得YOLOv3深度神經(jīng)網(wǎng)絡較其他粗定位算法具有更高的精度、更快的速度和更高的效率(Konan et al.,2020；袁靜等,2021b)，但需要配置高性能的GPU設備，且消耗的內(nèi)存資源高達233M.

總之，目前主流的閃電哨聲波識別算法需要將原始波形數(shù)據(jù)轉(zhuǎn)化為時頻圖，對算力和存儲設備的要求較為苛刻，適合離線數(shù)據(jù)，因而無法直接應用于星載.張衡一號感應磁力儀(Search Coil Magnetometer,SCM)載荷探測到的閃電哨聲波形數(shù)據(jù)以語音方式播放出來后能夠非常清晰地聽到類似口哨的聲音，意味著其頻率正好在人耳可聽到的范圍之內(nèi)(Wicks et al.,2016).梅爾頻率倒譜系數(shù)(MFCCs)恰恰模擬了人耳的聽覺特性(Davis and Mermelstein,1980)能有效提取閃電哨聲波的聽覺特征,LSTM神經(jīng)網(wǎng)絡引入了時間維信息適合處理聽覺特征中的聲音事件(Hochreiter and Schmidhuber,1997),這使得采用語音識別技術(shù)識別閃電哨聲波成為可能.因此，本文以張衡一號衛(wèi)星SCM觀測的磁場波形數(shù)據(jù)為研究對象，首次提出了基于智能語音技術(shù)的閃電哨聲波自動識別算法.

1 數(shù)據(jù)處理流程

2018年2月，我國首顆電磁衛(wèi)星張衡一號(ZH-1)發(fā)射成功，具備了天基觀測閃電哨聲波的能力.ZH-1衛(wèi)星覆蓋南北緯65°，在中國大陸及周邊1000 km區(qū)域及全球兩個地震帶(太平洋地震帶和歐亞地震帶)進行詳查模式的觀測，其他區(qū)域為巡查模式(Zhang et al.2018).ZH-1衛(wèi)星在軌飛行高度約507 km，其位置接近電離層頂部和等離子層邊界，這個區(qū)域有豐富的ELF/VLF頻段波動事件(Yan et al.2018;Wang et al.,2018)，如閃電哨聲波、準周期輻射等(Zhima et al.,2020).ZH-1軌道傾角97.4°，屬于太陽同步軌道，降交點地方時為下午2∶00；軌道回歸周期為5天，即每5天星下點軌跡相同；在一個回歸周期內(nèi)能夠?qū)崿F(xiàn)全球約500 km 空間分辨率的觀測.衛(wèi)星繞地球飛行一圈約94 min，大部分載荷在±65°緯度范圍開機工作，觀測數(shù)據(jù)按升軌(夜晚)和降軌(白天)分別存儲，每半軌(升/降軌)觀測約34 min；在同一天內(nèi)相鄰的升軌(或降軌)空間分辨率約2000 km.所搭載的感應磁力儀載荷(SCM)，通過法拉第電磁感應定律獲得電離層的感應磁場數(shù)據(jù)，能夠捕獲全球閃電哨聲波信號，其在巡查模式下僅獲得功率譜數(shù)據(jù)(Zhou et al.,2018;Huang et al.,2019).到目前為止，ZH-1已經(jīng)在軌觀測3年多，采集了大量的全球電磁場的波形和功率譜數(shù)據(jù)，其中SCM的三分量X/Y/Z包含3個頻段 ULF/ELF/VLF，頻點范圍ULF:10 Hz～200 Hz，ELF:200 Hz～2.2 kHz，VLF:12.5 Hz～20 kHz，原始波形數(shù)據(jù)的采樣率為51.2 kHz，功率譜數(shù)據(jù)的頻點間隔ULF:0.25 Hz，ELF:2.5 Hz，VLF：12.5 Hz，詳查模式VLF波形數(shù)據(jù)80 ms包含4096個點，每天產(chǎn)生大約10 G的數(shù)據(jù)量(Wang et al.,2018;Cao et al.,2018;Shen et al.,2018a,b).面對海量數(shù)據(jù)的挑戰(zhàn)，亟需實現(xiàn)基于原始波形的閃電哨聲波自動識別.

基于原始波形的閃電哨聲波識別方案主要由三部分組成，如圖2所示：數(shù)據(jù)整理、數(shù)據(jù)預處理和基于智能語音的閃電哨聲波自動識別算法.

圖2 自動識別ZH-1衛(wèi)星記錄閃電哨聲波的數(shù)據(jù)處理流程圖Fig.2 The diagram of automatic recognition of the lightning whistler recorded by ZH-1 satellite

(1)數(shù)據(jù)整理

本文的數(shù)據(jù)收集主要來自ZH-1衛(wèi)星2018年8月SCM載荷VLF波段的詳查數(shù)據(jù).首先，以0.16 s的時間滑動窗從原始波形數(shù)據(jù)中截取數(shù)據(jù)，該數(shù)據(jù)含有8192個點，將其轉(zhuǎn)化為音頻片段；再對截取數(shù)據(jù)進行重復短時傅里葉變換得到其時頻圖；然后根據(jù)時頻圖是否存在L色散形態(tài)特征進行人工標注；最終獲得10200段音頻數(shù)據(jù)集(閃電哨聲波數(shù)據(jù)5100段，非閃電哨聲波數(shù)據(jù)5100段).請注意本文中的時頻圖僅僅是為了查看是否存在閃電哨聲波，并不參與識別算法的計算.

(2)數(shù)據(jù)預處理

為了有效避免由于噪聲和信號的不穩(wěn)定導致的干擾，增強閃電哨聲波的波形特征，首先對原波形數(shù)據(jù)進行去趨勢處理，如式(1)所示：

(1)

其中，s(n)為原始信號，S(k)為去趨勢后的信號.結(jié)果如圖3所示.圖3a是含有閃電哨聲波的原波形，對其進行趨勢處理后的結(jié)果如圖3b所示，圖3c是不含閃電哨聲波的波形數(shù)據(jù)，對其進行去趨勢處理的結(jié)果如圖3d所示.

圖3 VLF磁場的原始波形和去趨勢處理(a)含閃電哨聲波的原始波形；(b)對圖(a)去趨勢處理；(c)不含閃電哨聲波的原始波形；(d)對圖(c)進行去趨勢處理.Fig.3 The raw and detrended waveform of the VLF magnetic field data(a)The raw wave containing the lightning whistler;(b)The result by processing the (a)with detrended method;(c)The raw wave not containing the lightning whistler;(d)The result by processing the (c)with detrended method.

(3)閃電哨聲波MFCCs音頻特征提取

由于人耳能明顯地聽到閃電哨聲波的嘶嘶聲，依據(jù)人耳的聽覺機理所設計的MFCCs能夠提取閃電哨聲波的聲音特征，將在第二節(jié)詳細介紹其提取過程.

(4)LSTM神經(jīng)網(wǎng)絡分類器

該環(huán)節(jié)主要包含訓練神經(jīng)網(wǎng)絡和應用神經(jīng)網(wǎng)絡兩個過程.其中，訓練神經(jīng)網(wǎng)絡指的是在訓練樣本集上提取MFCCs特征，利用該特征訓練LSTM神經(jīng)網(wǎng)絡；應用神經(jīng)網(wǎng)絡指的是在測試集上提取MFCCs，將其輸入訓練好的LSTM網(wǎng)絡，得到最終的識別結(jié)果,將在第三節(jié)詳細介紹其實現(xiàn)過程.

2 閃電哨聲波MFCCs音頻特征提取算法

MFCCs特征提取過程見圖4所示，主要包括預加權(quán)重，分幀加窗，快速傅里葉變換，Mel濾波組，對數(shù)運算，離散余弦變換(DCT)和動態(tài)差分.

圖4 MFCCs特征參數(shù)提取流程示意圖Fig.4 Flow chart for extracting MFCCs parameters

2.1 預加重、分幀、加窗和快速傅里葉變換

預加重處理：目的是對語音的高頻部分進行加重，增加高頻部分的分辨率.

s(n)=sn-μsn-1,(2)

式中sn是原始信號，s(n)為處理后的信號，參數(shù)μ的值介于0.9～1.0之間，由于SCM采樣率是51.2 kHz，則本文選取μ=0.97.

分幀處理：先將N個采樣點集合成一個觀測單位，稱為幀，本項目中N為512.為了避免相鄰兩幀的變化過大，設置兩相鄰幀之間存在重疊區(qū)域，此重疊區(qū)域包含了M個取樣點，通常M的值約為N的1/2或1/3.對應的時間長度是:

512/51200×1000=10 ms,(3)

加窗處理：窗長40 ms，窗移8 ms，對信號加窗以避免短時語音段邊緣的影響(Jibbs效應).加窗的定義如下：

sω(n)=s(n)×ω(n),(4)

式中ω(n)為窗函數(shù)，sω(n)為加窗后的信號，本文選用Hamming窗來進行加窗處理，ω定義如式(5)所示：

(5)

不同的α值會產(chǎn)生不同的漢明窗，默認選取0.46.

快速傅里葉變換：由于信號在時域上的變換通常很難看出信號的特性，所以通常將它轉(zhuǎn)換為頻域上的能量分布來觀察，不同的能量分布代表不同語音的特性.對分幀，加窗后的各幀信號進行快速傅里葉變換得到各幀的頻譜，并對語音信號的頻譜取模平方得到語音信號的功率譜，如式6所示：

(6)

其中，sω(n)為加窗后的信號，X(k)為快速傅里葉變換后得到的信號，N表示傅里葉變換的點數(shù).

2.2 Mel濾波器組

將能量譜通過一組Mel尺度的三角形濾波器組，采用一個有M個濾波器的三角濾波器，中心頻率為f(m).M為濾波器個數(shù)，本文M選默認值26，Hm(k)表示能量譜權(quán)重.其中Hm(k)：

Hm(k)=

(7)

其中f(m)滿足：

2Mel(f(m))=Mel(f(m-1))+Mel(f(m+1)),(8)

Mel標度頻率域提取出來的倒譜參數(shù)與頻率成非線性對應關(guān)系，見圖5所示，用式(9)近似表示為

圖5 Mel頻率與線性頻率的關(guān)系圖Fig.5 Relationship between the Mel frequency and the linear frequency

(9)

其中，f為頻率.

2.3 對數(shù)運算、離散余弦變換DCT

對數(shù)運算：將FFT得到的頻譜系數(shù)X(k)用順序三角濾波器進行濾波處理得到一組能量系數(shù)m1,m2,m3….濾波器組中每三個濾波器的跨度在Mel刻度上是相等的.所有的濾波器總體上覆蓋的范圍從0 Hz到采樣頻率的二分之一，計算能量系數(shù)s(m)的公式如下：

(10)

其中，X(k)為快速傅里葉變換后得到的信號，Hm(k)表示能量譜權(quán)重，M為濾波器的個數(shù).

計算濾波器組輸出能量系數(shù)的對數(shù)能量，其公式為

s′(m)=lns(m),0≤m≤M

(11)

其中，s(m)為能量系數(shù)，s′(m)為對數(shù)能量系數(shù).

離散余弦變換：目的是去除各維信號之間的相關(guān)性，將信號映射到低維空間.將上述的對數(shù)能量帶入離散余弦變換，求出L階的Mel-scale Cepstrum參數(shù)，如式(12)所示：

(12)

其中，c(n)為倒譜系數(shù)，L為階數(shù)在MFCC中通常選取8～13，本文L選取13.

2.4 動態(tài)差分

標準的倒譜參數(shù)MFCC只反映了語音參數(shù)的靜態(tài)特性，語音的動態(tài)特性可以用這些靜態(tài)特征的差分譜來描述，把語音的動、靜態(tài)特征結(jié)合起來能有效提高系統(tǒng)的識別性能.其差分譜的計算公式如式(13)所示：

(13)

其中,dt表示第t個一階差分，Ct表示第t個倒譜系數(shù)，L表示倒譜系數(shù)的階數(shù)，K表示一階導數(shù)的時間差，可取1或2.

最后將c(n)、dt(K=1)和dt(K=2)拼接得到一個16×39的二維張量，其每行表示一幀的能量值.幀能量的組成是由39維MFCC參數(shù)(13維MFCC倒譜系數(shù)+13維一階差分參數(shù)+13維二階差分參數(shù)).按照上述方法分別提取圖3中各子圖的MFCCs特征，將其繪制成如圖6所示的幀能量圖.橫坐標表示MFCC倒譜的個數(shù)，縱坐標代表時間.圖6a是含有閃電哨聲波的原始波形的MFCCs特征，對圖6a的原始波形進行去趨勢處理后的MFCCs特征如圖6b所示，圖6c是非含有閃電哨聲波的原始波形的MFCCs特征，對圖6c的原始波形進行去趨勢處理后的MFCCs特征如圖6d所示.觀察發(fā)現(xiàn)：閃電哨聲波和非閃電哨聲波在MFCCs特征上具有較強的鑒別性，比如MFCCs特征圖的第三列存在明顯差異.

圖6 MFCC-Time幀能量圖(a)含有閃電哨聲波；(b)對圖(a)去趨勢；(c)不含有閃電哨聲波；(d)對圖(c)去趨勢.Fig.6 The energy map of MFCC-Time Frame(a)Containing the lightning whistler;(b)The result by processing the (a)with detrended method;(c)Not containing the lightning whistler;(d)The result by processing the (c)with detrended method.

3 LSTM神經(jīng)網(wǎng)絡分類器算法

磁場原始數(shù)據(jù)是時序信號，當出現(xiàn)閃電哨聲波時，通過語音播放器能夠聽到明顯“哨聲”，提取其MFCCs特征發(fā)現(xiàn)其具有明顯的序列關(guān)聯(lián)性.LSTM 恰恰是一類專門用于處理序列關(guān)聯(lián)數(shù)據(jù)的神經(jīng)網(wǎng)絡.LSTM網(wǎng)絡中存在的細胞單元使該網(wǎng)絡具備“記住”長時間歷史信息的能力，從而可以通過學習捕捉到當前時刻的信息與歷史信息之間的依賴關(guān)系，以此提升對序列數(shù)據(jù)建模的能力，最終獲得較好的分類效果.鑒于以上優(yōu)勢，本文采用LSTM網(wǎng)絡對閃電哨聲波的MFCCs特征進行分類建模，其基本結(jié)構(gòu)見圖7所示，主要由4部分組成：細胞狀態(tài)、遺忘門f、輸入門i和輸出門o三種門控單元(Hochreiter and Schmidhuber,1997).

圖7 LSTM單元結(jié)構(gòu)示意圖Fig.7 The unit structure of LSTM

遺忘門：用來決定當前需要丟棄之前的哪些磁場信息.LSTM網(wǎng)絡通過學習決定讓網(wǎng)絡記住哪些磁場數(shù)據(jù)內(nèi)容，遺忘門定義如下：

ft=δ(Wt·[ht-1,xt]+bf),(14)

其中，xt是網(wǎng)絡在t時刻的輸入，其數(shù)值是MFCCs序列中的第t個數(shù)值；δ表示取Sigmoid函數(shù)值；Wt表示t時刻的權(quán)重矩陣；ht-1表示將1到t-1時刻的MFCCs序列數(shù)據(jù)(x1,x2,…,xt-1)輸入到 LSTM 神經(jīng)網(wǎng)絡后的輸出的抽象特征信息；bf表示t時刻的偏置量.ft的元素取值范圍是0到1，表示t時刻的遺忘程度，0表示全忘，1表示全記住.遺忘門具有過濾掉與閃電哨聲波無關(guān)的特征信息的功能.

輸入門：處理MFCCs序列數(shù)據(jù)中當前位置的數(shù)據(jù).圖7中間部分的左側(cè)是sigmoid激活函數(shù)，這個函數(shù)是用來決定哪些輸入的特征會被記住，其定義如式(15)所示.

it=δ(Wi·[ht-1,xt]+bi).

(15)

細胞狀態(tài)：細胞狀態(tài)是隱狀態(tài)，類似一個存儲信息的容器，保存序列數(shù)據(jù)中的關(guān)聯(lián)性信息，具有歷史記憶功能，其定義如式(16)所示.

Ct=ft×Ct-1+it×Rt,(16)

其中，Rt是輸入門右側(cè)的tanh部分，該部分的主要作用是利用上一次的輸出和本次輸入計算新的候選的細胞狀態(tài)信息，其定義是：Rt=tanh(Wc·[ht-1,xt]+bc)；ft×Ct-1意味著通過前一個時刻的細胞狀態(tài)Ct-1與遺忘門ft結(jié)合，忘記舊細胞狀態(tài)中一些信息；it×Rt意味著通過輸入門it和候選細胞狀態(tài)Rt結(jié)合，記住新候選細胞狀態(tài)中一些信息.

輸出門：控制當前隱狀態(tài)的輸出信息，其基本計算過程為

ot=δ(Wo·[ht-1,xt]+bi),(17)

其中，ot是輸出門，通過將前一個時刻的輸出信息ht-1和當前的數(shù)據(jù)xt輸入到Sigmoid函數(shù)得到的.

最后輸出當前時刻的重要信息ht，其定義如下：

ht=ot×tanh(Ct),(18)

該信息是通過聯(lián)合輸出門和當前的細胞狀態(tài)Ct進行計算而得到.最后將ht輸入到sigmoid函數(shù)實現(xiàn)自動識別.ht含有序列之間的關(guān)聯(lián)性信息，這使得序列數(shù)據(jù)的識別效果得到了大幅提升.

4 實驗和分析

4.1 實驗流程及LSTM神經(jīng)網(wǎng)絡模型參數(shù)設置

實驗流程包括數(shù)據(jù)整理、MFCCs特征提取、LSTM模型訓練和評估指標值，并進行1000次重復實驗.具體詳細步驟如下：

(1)數(shù)據(jù)集：含有5100個閃電哨聲波波形樣本集WD和5100個非閃電哨聲波的樣本集NWD.

(2)訓練集：分別從樣本集WD和NWD中各隨機選取50%作為訓練樣本，構(gòu)建訓練集.

(3)測試集：將樣本集WD和NWD剩下的樣本組建成測試集.

(4)特征提取：用4種不同的特征提取方法提取訓練集的音頻特征，4種特征分別是：原始波形數(shù)據(jù)特征，用Original表示；對原始波形數(shù)據(jù)進行去趨勢處理后的特征，用Original_Detrend表示；對原始波形數(shù)據(jù)采用MFCCs處理后的特征，用Original_MFCC表示；先對原始波形數(shù)據(jù)進行去趨勢處理，再采用MFCCs處理后的特征，用Original_Detrend_MFCC表示.

(5)訓練過程：用(4)中提及的4種特征分別訓練LSTM分類模型，得到4種不同的LSTM分類器.

(6)測試過程：在測試集上先采用步驟(4)中4種不同的特征提取方式提取特征，再將4種特征分別輸入到4種不同的LSTM分類模型中進行識別，輸出識別結(jié)果.

(7)評估：對識別效果采用4種指標進行評估：精確度(Precision)、召回率(Recall)、F1值和ROC面積(AUC-ROC)(袁靜等，2021b).

采用不同的輸入特征訓練LSTM神經(jīng)網(wǎng)絡所需的超參數(shù)也不盡相同，對每一種特征分類器，均采用十折交叉的方式獲得LSTM神經(jīng)網(wǎng)絡模型的超參數(shù)，見表1所示.其中優(yōu)化器選擇的是自適應梯度優(yōu)化器(Adaptive Gradient Algorithm,Adagrad)(Duchi et al.,2011)；損失函數(shù)選擇的是適用于二分類的二元交叉熵損失函數(shù)(Binary_CrossEntropy)(Goodfellow et al.,2017).

表1 基于4種不同特征的LSTM神經(jīng)網(wǎng)絡的參數(shù)Table 1 The parameters setting of four LSTM neural networks

為了全面地評估本文算法的有效性，開展1000次實驗：每一次實驗的訓練集和測試集均使用4種不同的特征提取方法提取圖像特征，并以此訓練得到4個不同的分類器，最后采用精度、召回率、F1和AUC-ROC評估每一個分類器的性能，4種指標的詳細定義請參考袁靜等(2021a)的文獻資料.由于每次的訓練集和測試集不同，單次的4個評估指標難以充分評價本文提出的閃電哨聲波識別算法的效果，因此，開展實驗1000次，并在4個評估指標的基礎上制定了如下的評價策略：

(1)訓練集和測試集的表現(xiàn)；

(2)部分識別結(jié)果展示；

(3)總體識別精度評價：對1000次實驗的評價指標進行求平均值的評價策略；

(4)穩(wěn)定性評價和差異性評價：對1000次實驗的評價指標采用盒形圖評估分類的穩(wěn)定性.為了評價不同的特征分類器之間是否具有明顯差異，采用 T 檢驗進行差異性評價.閾值p=0.05，即小于0.05為具有明顯差異，若大于0.05表明不具有明顯差異.

4.2 訓練集和測試集的表現(xiàn)

利用表1的超參數(shù)開展1000次實驗，得到每次實驗的測試集和訓練集在精度指標上的表現(xiàn)，對其進行求平均得到如表2所示的結(jié)果，說明本文提出的算法無論在訓練集和測試集的表現(xiàn)均為最佳.

表2 基于4種不同特征的LSTM神經(jīng)網(wǎng)絡在訓練集和測試集上的表現(xiàn)Table 2 The performance of four LSTM neural network models on training set and test set

4.3 部分識別結(jié)果展示

部分識別結(jié)果的波形圖和相應的時頻圖繪制如圖8和圖9所示，其中波形圖是識別結(jié)果，此處的時頻圖僅僅是為了可視化波形中是否存在閃電哨聲波.圖8是閃電哨聲波的識別結(jié)果，圖8a是正確識別出的閃電哨聲波，圖8b是未識別出的閃電哨聲波，識別不成功的原因是：閃電哨聲波的能量較弱且背景干擾強，導致閃電哨聲波的趨勢特征不明顯，采用去趨勢處里后哨聲波趨勢被干擾淹沒，造成了不存在閃電哨聲波的假象，如圖8c所示.圖9是非閃電哨聲波的部分識別結(jié)果，(a)是正確識別出的非閃電哨聲波，(b)是誤識別的非閃電哨聲波，誤識別的原因是：原始波形數(shù)據(jù)中存在強烈的干擾信號且出現(xiàn)近似閃電哨聲波的趨勢特征，如圖9b的時頻圖的黑框處，對相應的波形進行去趨勢處理后的結(jié)果如圖9c所示.

圖8 閃電哨聲波(a)正確識別(右子圖是原始波形，左子圖是其時頻圖)；(b)未識別(右邊是原始波形，左邊是其時頻圖)；(c)對圖(b)進行去趨勢.Fig.8 Results of recognition of lightning whistler(a)Accurate (The right panel is the wave and the left panel is the time-frequency plot of the right panel);(b)Wrong (The right panel is the wave and the left panel is the time-frequency plot of the right panel);(c)The result by processing the (b)with detrended method.

圖9 非閃電哨聲波(a)正確識別(右子圖是波形數(shù)據(jù)，左子圖是其時頻圖)；(b)錯誤識別(右子圖是波形數(shù)據(jù)，左子圖是其時頻圖)；(c)去趨勢后的(b).Fig.9 Results of recognition of non lightning whistler(a)Accurate (The right panel is the wave and the left panel is the time-frequency plot of the right panel);(b)Wrong (The right panel is the wave and the left panel is the time-frequency plot of the right panel);(c)The result by processing the (b)with detrended method.

4.4 總體識別精度評價

1000次實驗后分別獲得1000個精確度(Precision)、召回率(Recall)、F1值(F1socre)、ROC面積(AUC-ROC)值、時間消耗(Cost time)和內(nèi)存消耗(Cost memory)，分別對其進行均值計算以評估基于智能語音技術(shù)的哨聲波識別效果，如表3所示.

表3中Original+LSTM表示用原始波形訓練LSTM分類器，Original_Detrend+LSTM表示對原始波形進行去趨勢處理后再訓練LSTM分類器，Original_MFCC+LSTM表示對原始波形提取MFCCs特征，Original_Detrend_MFCC+LSTM表示對原始波形進行去趨勢處理后再提取其MFCCs特征，最后用該特征訓練LSTM分類器.通過觀察表3發(fā)現(xiàn)：直接用原始波形數(shù)據(jù)訓練LSTM分類器的識別算法(Original+LSTM),具有最少的時間消耗和內(nèi)存消耗，分別是2.08 s和82.790 MB，但該算法在分類精度、召回率、F1值和AUC-ROC四個指標上的表現(xiàn)最差.本文提出識別算法(Original_Detrend_MFCC+LSTM)在4個指標上的表現(xiàn)效果最佳，分別達到0.967，0.842，0.900和0.907，且由于采用了MFCCs特征，將每段0.16 s的音頻數(shù)據(jù)量從8192減少到684(16×39)，使得其時間消耗和內(nèi)存消耗與Original+LSTM相近，達到2.24 s和82.89 MB.而Original_Detrend+LSTM算法為了得到較好的分類結(jié)果使用了雙層LSTM網(wǎng)絡，導致該算法損失了更多的時間和和內(nèi)存.值得注意的是目前最佳的基于時頻圖的閃電哨聲波識別算法采用的是YOLOV3深度卷積神經(jīng)網(wǎng)絡(袁靜等，2021b)，其在CPU上消耗的時間成本是6.71 s，消耗的內(nèi)存資源是233 MB.

表3 1000次實驗后平均效果Table 3 Statistical results of 1000 experiments

總之，在基于原始波形的閃電哨聲波識別中，聯(lián)合MFCCs音頻特征提取和LSTM神經(jīng)網(wǎng)絡技術(shù)的閃電哨聲波識別算法的分類效果最優(yōu).與基于時頻圖的識別算法比較，其消耗的時間成本和內(nèi)存資源最小.

4.5 穩(wěn)定性和差異性評價

本小節(jié)將針對不同的LSTM分類器的分類效果進行穩(wěn)定性和差異性評價.

(1)穩(wěn)定性評價：對每種指標的1000個數(shù)據(jù)，繪制其箱型圖，如圖10所示.閃電哨聲波的識別精度(Precision)的1000組數(shù)據(jù)分布圖如圖10的Precision圖所示：其橫軸是不同的特征分類器，縱軸是精度.可發(fā)現(xiàn)采用Original_Detrend_MFCC特征分類器的Precision箱體的高度低于Original、Original_Detrend和Original_MFCC特征分類器，說明該特征分類器在Precision指標上的表現(xiàn)更穩(wěn)定；該箱體的位置高于Original、Original_Detrend和Original_MFCC特征分類器，說明該特征分類器在Precision指標上的表現(xiàn)更優(yōu).用上述方法觀察圖10的Recall、F1score和AUC-ROC的箱體，均能得出上述同樣結(jié)論.總之，本文提出的分類器在4個評價指標上均最優(yōu)且最穩(wěn)定性.

圖10 評估LSTM分類器性能的箱型圖Fig.10 Box plot for evaluating the performance of four LSTM network models

(2)差異性評價：為檢驗不同分類器的性能是否存在明顯的差異，采用兩兩獨立樣本的T檢驗方法對顯著性差異進行定量評價，其顯著性P值越高，表明差異性越小，通常采用的閾值是0.05，其含義是若差異性小于0.05則認為存在明顯差異；若大于0.05則認為兩組實驗不存在明顯差異.結(jié)果見表4所示.

觀察表4的Precision 的T值檢驗表中的第一行第二列的值是0，說明采用Original特征分類器與Original_Detrend特征分類器在精度指標上存在明顯差異，繼續(xù)觀察發(fā)現(xiàn)，Original_Detrend_MFCCs特征分類器的精度與其他兩個特征分類器的精度也存在明顯的差異；但相同特征分類器在精度指標上的T檢驗值是1，說明相同特征分類器在精度指標上未表現(xiàn)出顯著性差異，比如第三行第三列的值是1.用上述方法觀察表4的Recall的T值檢驗表、F1score的T值檢驗表和AUC-ROC的T值檢驗表，可以得到同等結(jié)論.總之，本文提出的識別算法與其他識別算法在精度、召回率、F1值和AUC-ROC 4個評價指標均表現(xiàn)出了明顯差異，說明本文提出的算法顯著地改善了識別效果.

表4 T檢驗Table 4 The T test

5 討論

上述實驗表明本文提出的閃電哨聲波自動識別算法具有一定的效果.算法方案中的原始波形特征提取和LSTM神經(jīng)網(wǎng)絡對閃電哨聲波自動識別結(jié)果有非常重要的影響，本章將對其產(chǎn)生的影響進行較深入的討論和分析.

5.1 提取波形特征

為了分析4種特征的時間軌跡，本小節(jié)隨機選用10個閃電哨聲波樣本和 10 個非閃電哨聲波樣本的音頻數(shù)據(jù)，繪制4種波形特征的時間軌跡如圖11所示.

圖11a是原始波形數(shù)據(jù)的時間序列，圖11b是對原始波形數(shù)據(jù)進行去趨勢處理后的時間序列，圖11c是對原始波形數(shù)據(jù)進行MFCCs特征提取后的時間序列，對圖11b進行MFCCs特征提取，得到16×39維的特征圖，參考圖6的效果，將該特征圖的第三列特征按照時間順序展開，得到如圖11d所示的結(jié)果；其中W表示閃電哨聲波波形樣本，NonW表示非閃電哨聲波波形樣本.觀察圖11a發(fā)現(xiàn)，含有閃電哨聲波原始波形的時間軌跡和非閃電哨聲波的軌跡雜糅在一起，增加了分類的難度；對原始波形數(shù)據(jù)進行去趨勢處理后，含有閃電哨聲波的樣本軌跡的類內(nèi)差變小，如圖11b所示，但分類難度依然較大；對原始波形數(shù)據(jù)提取MFCCs特征后，含有閃電哨聲波的樣本軌跡具有可分性，但有部分雜糅在一起且內(nèi)類差較大，如圖11c所示；對原始波形進行去趨勢處理并提取MFCCs特征，繪制其時間軌跡如圖11c所示，發(fā)現(xiàn)：含有閃電哨聲波的時間軌跡和非閃電哨聲波的軌跡具有明顯的可分性.

圖11 波形數(shù)據(jù)中的語音特征的時間軌跡(a)原數(shù)據(jù);(b)原數(shù)據(jù)去趨勢；(c)原數(shù)據(jù)+MFCCs特征；(d)原數(shù)據(jù)+去趨勢+MFCCs特征；橙色和藍色曲線分別代表閃電哨聲波和非閃電哨聲波.Fig.11 The visualization of the audio feature,the orange and blue curves represent the lightning induced whistler waves and non-whistler waves(a)The audio feature from the raw waveform (Marked as Original);(b)The audio feature from the detrended raw waveform (Marked as Original_detrend);(c)The audio feature extracted from MFCCs of the raw wave (Marked as Original_MFCCs)；(d)The audio feature extracted from MFCCs of the detrended wave (Marked as Original_Detrend_MFCCs).

5.2 LSTM神經(jīng)網(wǎng)絡的抽象映射

LSTM神經(jīng)網(wǎng)絡的輸出門的最后時刻的隱藏信息特征ht含有時間序列的抽象特征，該抽象特征包括該時間序列的歷史信息和趨勢信息等，對最終的分類結(jié)果產(chǎn)生關(guān)鍵影響.本小節(jié)在測試集中隨機選擇200個閃電哨聲波波形樣本(WD)和200個非閃電哨聲波波形樣本(NWD)，將這些樣本通過4種不同特征的LSTM分類器的抽象特征繪制成圖12，其中W表示閃電哨聲波波形樣本，NonW表示非閃電哨聲波波形樣本.

根據(jù)這些樣本數(shù)據(jù)，計算不同的LSTM網(wǎng)絡的抽象特征的類內(nèi)差異度和類間差異度，得到見表5所示的結(jié)果.

通過觀察圖12和表5可以發(fā)現(xiàn)，基于MFCCs的LSTM分類器的抽象特征，其閃電哨聲波類內(nèi)差是0.06609，非閃電哨聲波的類內(nèi)差是0.00024，這說明同類樣本的聚集性強；同時，閃電哨聲波和非閃電哨聲波的類間差是0.26357，說明不同類間的差異性強.總之，本文提出的算法具有類內(nèi)差小、類間差大的特點，意味著該算法更容易實現(xiàn)準確分類.

圖12 不同LSTM網(wǎng)絡中的抽象特征(a)原數(shù)據(jù)+LSTM網(wǎng)絡;(b)原數(shù)據(jù)去趨勢+LSTM網(wǎng)絡；(c)原數(shù)據(jù)+MFCCs+LSTM網(wǎng)絡；(d)原數(shù)據(jù)+去趨勢+MFCCs特征+LSTM網(wǎng)絡.Fig.12 The abstract features from four LSTM networks(a)Original+LSTM;(b)Original+Detrend+LSTM；(b)Original+MFCCs+LSTM；(d)Original+Detrend+MFCCs+LSTM.

表5 4種LSTM網(wǎng)絡中的抽象特征的差異Table 5 Comparison of the abstract features among four LSTM networks

5.3 網(wǎng)絡結(jié)構(gòu)對LSTM分類結(jié)果的影響

從實驗結(jié)果發(fā)現(xiàn)不同的網(wǎng)絡結(jié)構(gòu)對LSTM神經(jīng)網(wǎng)絡的分類效果產(chǎn)生不同的影響，接下來將從LSTM網(wǎng)絡的性能和抽象特征可分性兩個方面進行討論.

(1)性能評估

本小節(jié)對不同網(wǎng)絡結(jié)構(gòu)的LSTM網(wǎng)絡進行十折交叉驗證并計算交叉驗證得分和訓練所需的時間消耗，結(jié)果見表6所示.

表6 不同超參數(shù)的LSTM網(wǎng)絡Table 6 LSTM network with different parameters

通過觀察表6發(fā)現(xiàn)，使用加入Dropout的單層LSTM網(wǎng)絡(LSTM網(wǎng)絡C)比另外兩個LSTM網(wǎng)絡在十折交叉驗證中得分高，達到0.947，其分類器的平均時間消耗也是最少，為43.673 s.由于神經(jīng)網(wǎng)絡會發(fā)生過擬合現(xiàn)象，加入Dropout能起到一種Vote的作用和減少神經(jīng)元之間的共適應性，能提升了網(wǎng)絡的精度和泛化能力.

(2)特征可分角度分析

為了定性的比較不同超參數(shù)下的LSTM網(wǎng)絡的抽象特征是否具有可分性，本小節(jié)將60個閃電哨聲波和60個非閃電哨聲波分別輸入到不同超參數(shù)的LSTM網(wǎng)絡中提取隱藏信息特征ht，并將其隨著時間變化的軌跡繪制如圖13所示.其中W表示閃電哨聲波波形樣本數(shù)據(jù)，NonW表示非閃電哨聲波波形樣本數(shù)據(jù).

圖13a的時間軌跡是來自加入失活層的雙層LSTM網(wǎng)絡，圖13b的時間軌跡來自加入失活層的單層LSTM網(wǎng)絡，圖13c的時間軌跡來自去掉失活層的單層LSTM網(wǎng)絡.觀察圖13發(fā)現(xiàn)，LSTM網(wǎng)絡-A的W和NonW的特征存在重疊，總體區(qū)分度不是很高；LSTM網(wǎng)絡-B的W和NonW的特征區(qū)分度相對較高，仍然存在不同類別的特征交錯的情況；相比之下，LSTM網(wǎng)絡-C的W和NonW的特征主要分布在兩個不同區(qū)域，特征交錯較少.該現(xiàn)象說明LSTM網(wǎng)絡-C能夠提高閃電哨聲波和非閃電哨聲波的區(qū)分度，具有較強的閃電哨聲波識別能力.

圖13 不同LSTM網(wǎng)絡隱藏特征的時間序列(a)LSTM網(wǎng)絡A；(b)LSTM網(wǎng)絡B；(c)LSTM網(wǎng)絡C.Fig.13 The plots of the time series of hidden features in four LSTM networks(a)LSTM A；(b)LSTM B；(c)LSTM C.

6 結(jié)論

在ZH-1衛(wèi)星運行過程中，其搭載的SCM每天產(chǎn)生大約10GB的數(shù)據(jù)量，其中絕大部分是不存在閃電哨聲波的數(shù)據(jù).如何實現(xiàn)星載實時識別閃電哨聲波、傳回更有意義的數(shù)據(jù)、減少存儲壓力變得尤為重要.鑒于閃電哨聲波的頻率范圍在人耳聽覺范圍之內(nèi)，基于語音識別技術(shù)的閃電哨聲波自動識別已經(jīng)成為可能.

本文在ZH-1號衛(wèi)星的SCM數(shù)據(jù)上開展星載閃電哨聲波自動識別算法的探索和研究.根據(jù)閃電哨聲波能被人耳聽到的特性，采用了MFCCs特征提取方式增強閃電哨聲波的聽覺特征，并采用淺層長短期記憶(LSTM)回歸神經(jīng)網(wǎng)絡對特征進行分類，其分類結(jié)果在精度，F(xiàn)1socre以及AUC值指標上均高于90%，同時將此方法與基于YOLOv3神經(jīng)網(wǎng)絡的閃電哨聲波檢測算法進行對比發(fā)現(xiàn)：準確率相當，卻能夠節(jié)省150.11 MB的存儲空間以及4430 ms的時間消耗，極大地增加了星載實時識別閃電哨聲波的可能性.MFCCs是模擬人耳聽覺特性所設計的，但由于Hz-Mel頻率非線性的對應關(guān)系，使得在低頻區(qū)域使用的濾波器數(shù)量較多，分布密集，而中高區(qū)域使用的濾波器較少，分布稀疏.使得MFCCs隨著頻率的提高其計算的精度就隨之下降.而閃電會產(chǎn)生強烈的寬帶無線電波，尤其是在300 Hz至20 kHz的甚低頻(VLF)頻帶中，因此后續(xù)將通過增加中高頻的濾波器的方式進一步改善MFCCs特征，從而提高其高頻部分的計算精度.

致謝本工作使用了中國國家航天局和中國地震局支持的張衡一號衛(wèi)星的觀測數(shù)據(jù)(http:∥leos.ac.cn),特別感謝來自應急管理部國家自然災害防治研究院的張衡一號衛(wèi)星團隊的所有成員為本文研究數(shù)據(jù)提供的技術(shù)服務支持.