基于電磁衛星的閃電哨聲波智能檢測算法的研究進展

袁靜, 王橋, 張學民, 楊德賀* , 王志國, 張樂, 申旭輝, 澤仁志瑪

1 防災科技學院， 河北 廊坊 062541 2 應急管理部國家自然災害防治研究院， 北京 100085 3 中國地震局地震預測研究所， 北京 100036 4 清華大學， 北京 100084 5 中國電信股份有限公司研究院， 北京 102209

0 引言

高空閃電發生時產生寬頻帶電磁波，電磁波中含有高頻成分和低頻成分.由于頻率高的群速度大，故高頻成分比低頻成分先到達，導致接收喇叭中出現下降的聲調，像哨聲一樣，故稱其為哨聲波.有的哨聲波，到達電離層后，沿著地磁力線傳播(如圖1所示)到達地球另一端共軛點，稱之為短哨聲波.還有一部分被電離層反射回到閃電發生地，稱之為長哨聲波(Chen et al., 1996).目前觀測閃電哨聲波主要有兩種方式：地面甚低頻(Very low frequency, VLF)觀測站和空中飛行器.最著名的地面VLF觀測站是電離層閃電全息成像陣列(張義軍和張陽, 2016; Johnson, 2000)，該陣列覆蓋了美國中西部閃電影響的絕大部分區域，通過觀測VLF信號，對閃電伴隨的電離層的擾動進行深入研究，其研究主要涉及到閃電引起的電離層變化特征、觀測方法、相互作用模式以及引起的中高層空間放電現象等方面(Johnson, 2000; Kucherov and Nikolaenko, 1979; Lay and Shao, 2011; Marshall et al., 2008a,b; Inan et al., 1991).因此，地面VLF觀測站能夠通過閃電哨聲波現象研究閃電引起的電離層擾動，其缺點是地面觀測站具有空間的局限性，且觀測信號容易受到地面干擾，比如地鐵、工地建設等人類活動的影響.衛星捕獲的閃電哨聲波避免了空間位置的限制且不容易受地面干擾影響，因此，基于電磁衛星的閃電哨聲波更適合開展電離層方面的研究和探索.

圖1 閃電哨聲波傳播路徑Fig.1 Propagation path of a lightning whistler

1993年，地震電離層擾動理論的出現推動了地震電磁衛星的發展(Parrot et al., 1993).法國2001年正式對外公布 DEMETER衛星計劃，2004年成功發射，2010年底停止運行.DEMETER衛星特定的軌道設計方案，使其成為全球第一顆專門用于地震和火山引起的電離層擾動監測的電磁衛星，是國際上地震電磁空間探測一個里程碑式的發展標志.隨后，在經歷了長期的科研論證后，中國于2013年正式啟動電磁衛星計劃，2018年2月2日，ZH-1衛星成功發射，2018年11月完成在軌測試，是國內首個地球物理場探測衛星，主要用于獲取全球空間電磁場、電磁波、電離層等離子體、高能粒子沉降等科學數據，是提取與地震前兆相關的空間環境信息的重要觀測手段.所獲取的探測資料有利于地震電離層及空間物理相關現象研究的開展，成為世界同行關注的熱點(申旭輝等, 2018; 張學民等, 2020).

ZH-1衛星的感應磁力儀(Wang et al., 2018)和電場探測儀(Huang et al., 2018a,b)均能捕獲到閃電哨聲波現象，其在VLF波段的電磁場數據的時頻圖中呈現的是色散狀，如圖2和圖3所示.產生該形態特征的物理機制是：高頻的電磁波傳播的速度快，先到達；低頻的電磁波的速度慢，后到達，反映在時頻圖上是頻率隨著時間逐漸減少，稱之為“色散譜”.當路徑長度較長或沿傳輸路徑的電子密度較密集時，色散譜變大.由此說明，色散形態與傳播路徑的長度和傳播路徑中的等離子體的電子密度分布有關.除此之外，閃電哨聲波在傳播的過程中攜帶了其他大量的地球空間物理的相關信息，因此，電磁衛星的閃電哨聲波是開展空間物理研究的重要媒介和手段.例如：Záhlava等(2018a)分析了DEMETER和RBSP對地球內部磁層哨聲模式波的縱向依賴性的測量結果，結果表明等離子體層內部的哨聲模式波的縱向依賴性很強.Bayupati等(2012)分析AKEBONO衛星觀測到的閃電哨聲波的色散形態，討論了閃電哨聲波沿軌道傳播的時間與電子密度分布的關系，表明了閃電哨聲波的色散趨勢是確定等離子體層中整體電子密度分布的有力方法.Oike等(2014)分析了AKEBONO衛星探測到的閃電哨聲波發生頻率與地面觀測到的閃電活動的空間分布和時間變化，表明了在電離層中閃電哨聲波的出現與閃電活動以及地球周圍的電子密度分布密切相關.Clilverd等(2002)使用全球閃電定位網(World-Wide Lightning Location Network, WWLLN)來確定閃電的來源位置，然后根據閃電哨聲波的傳播時間，遠程感測電子密度分布頻率，借助傳播理論導出電子密度沿傳播路徑的函數.Parrot等(2019)通過閃電哨聲波發現：全球溫度變化1°，閃電發生頻率變化5%～6%.Záhlava等(2018b)通過觀測閃電哨聲波色散形態的衰減程度發現晚上衰減慢，并得到晚間的空間電離層的濃度低的結論.此外，由于空中飛行的衛星的電子元件容易受高能電子的總劑量效應和單粒子效應影響而受損，借助哨聲波數據構造哨聲波參數統計模型，開展基于哨聲波模型的高能電子通量預測,對衛星設計和防護具有重要的意義(Horne et al., 2013; 羅旭東等, 2017).

圖2 張衡一號衛星閃電哨聲波形態圖例Fig.2 The time-frequency spectrogram of the lightning whistler observed by ZH-1 satellite

圖3 DEMETER衛星閃電形態圖例(Parrot et al., 2015)Fig.3 The time-frequency spectrogram of the lightning whistler observed by DEMETER (Parrot et al., 2015)

綜上，電磁衛星能夠捕獲到大量的閃電哨聲波，為借鑒閃電哨聲波的形狀特征和色散程度等開展地震電離層和空間物理研究提供了可靠的數據支撐.然而，目前通過人工識別閃電哨聲波現象的方式阻礙了相關研究的進展(Ferencz et al., 2001;Lichtenberger et al., 2008)，其主要原因是：

(1)電磁衛星高精度和全天候的工作模式導致電磁衛星存檔的數據體量很大，比如ZH-1衛星日產的數據量高達20 G.

(2)閃電頻發使得電磁衛星捕獲的閃電哨聲波現象的數量很大.閃電是在地球空間發生最為頻繁的一種自然現象，1925年，Brooks等(1925)使用雷暴日資料，估計全球在任一時刻存在1800個閃電.2003年Christian指出全球閃電每秒鐘45次.此外，據WWLLN(http:∥www.wwlln.com/)記載：僅2017年閃電發生了2.08億次.2019年7月28日，全球閃電分布如圖4所示，發生雷電約64萬個.

圖4 2019年7月28日全球閃電密度分布圖 (http:∥dudwlln.otago.ac.nz/)Fig.4 A global lightning density map on July 28, 2019 (http:∥dudwlln.otago.ac.nz/)

面臨數據層面的挑戰，近年來，隨著計算機視覺和深度學習等技術的發展，國內外學者借助人工智能技術開始了從電磁衛星數據中自動檢測閃電哨聲波的研究，閃電哨聲波傳播的物理機制決定了其在時頻圖中特征是具有明顯的色散形態，閃電哨聲波智能檢測的目標是從時頻圖中識別色散形態并進行定位.該問題的本質屬于計算機視覺中的目標檢測問題.目前發表的基于閃電哨聲波智能檢測算法涉及到空間物理、信號處理、圖像處理、模式識別和機器學習等多學科交叉，并取得了一定的效果.

本文將對閃電哨聲波智能檢測算法進行整理和對比分析，創新點主要體現在：首次綜述基于電磁衛星數據的閃電哨聲波智能檢測算法，對不同的算法進行分類整理和分析，為開展基于人工智能技術的海量電磁衛星數據處理提供參考.全文的內容組織安排如下：第一節介紹閃電哨聲波的時頻特征和類型；第二節介紹智能檢測算法的基本內容，分別從特征提取、分類(識別)和定位三個方面對不同的算法進行整理和對比分析；第三節介紹智能檢測算法的評價標準；第四節對部分算法進行復現；第五節總結當前閃電哨聲波智能檢測算法的現狀，展望閃電哨聲波智能檢測研究的發展趨勢和應用前景.

1 閃電哨聲波的時頻特征和類型

目前，人工識別電磁數據中的閃電哨聲波的過程通常包含三個環節：對電場數據/磁場數據進行預處理；采用重疊的快速傅里葉變換將預處理后的數據轉化成時頻圖；在時頻圖中，尋找色散形態區域.以此為依據，在DEMETER衛星中發現了大量的閃電哨聲波的降調形態(Helliwell, 1965; 胡云鵬等, 2020; Santolík et al., 2009)，如圖3所示.Helliwell(1965)根據色散程度將地面VLF觀測站觀測到的閃電哨聲波分為8種不同的類型，如表1所示.Santolík等(2009)指出地面VLF觀測站記錄的閃電形態與衛星觀測到的閃電形態極其相似且均具有頻散特性.不同類型的閃電哨聲波均與閃電源、傳播路徑的長短、電子密度等物理空間環境有關，是探測等離子層及輻射帶的重要媒介.值得注意的是目前的閃電哨聲波識別算法均適用于識別單條線的閃電哨聲波類型.

2 閃電哨聲波智能檢測算法

閃電哨聲波自動識別最早可追溯到2008年，Lichtenberger等(2008)指出通過哨聲波監測等離子體層電子密度變化遇到的瓶頸問題是需要大量人力從觀測數據中去提取和辨識哨聲波并分析哨聲波的軌跡.因此，作者在蒂豪尼15個地面觀測站的數據的基礎上，提出了一個自動檢測哨聲波的算法(Automatic whistler detector and analyzer system, AWD).該方法開辟了閃電哨聲波智能檢測算法的先例，并比較成功的應用到從地面VLF觀測站中自動識別閃電哨聲波，該方法對于識別分數跳(表1中的第8個類型)的閃電哨聲波效果較佳.

閃電哨聲波智能檢測算法的本質是屬于計算機視覺中的目標檢測，旨在從圖像中尋找感興趣的區域，提取該區域的特征，準確判斷該特征的類別，并進行定位.因此，完整的閃電哨聲波智能檢測算法的流程如圖5所示：(1)收集電磁場數據，通常是VLF波段的電場或磁場波形數據，如圖6所示；(2)對數據進行濾波等預處理；(3)采用快速傅里葉變換處理電磁場數據得到時頻圖，如圖7所示；(4)采用特征提取方法從時頻圖中提取能夠表示閃電哨聲波的特征；(5)對提取到的特征進行分類(識別)；(6)用邊界框對識別的目標進行定位.

圖5 閃電哨聲波智能檢測的流程Fig.5 The overall flow chart of the lightning whistlers detection

圖6 張衡一號感應式磁力儀測得的快速變化磁場數據Fig.6 The fast variation data of geomagnetic field obtained by Search Coil Magnetometer onboard ZH-1 satellite

圖7 張衡一號感應式磁力儀快速變化磁場數據的時頻圖Fig.7 The time-frequency spectrogram of the Search Coil Magnetometer data in Fig.6

在上述流程中，閃電哨聲波的特征提取、分類(識別)和定位是算法的核心環節.接下來分別從特征提取、分類(識別)和定位三個方面對現存的智能檢測識別算法進行梳理和對比.

2.1 特征提取

特征提取是計算機視覺智能識別過程中的關鍵環節.其主要作用是將智能模型不能識別的原始數據轉化為智能模型可以識別的數值特征的過程.比如說，圖像是由一系列像素點 (原始數據) 構成的，這些像素點本身無法被智能模型直接使用，需要將這些像素點轉化成角點、紋理等形式(數值特征)，供智能模型使用.特征的好壞對泛化性能有至關重要的影響(周志華, 2016).對閃電哨聲波進行特征提取旨在數值化閃電哨聲波的色散形態特征.本小節分別從人工提取特征和自動提取特征兩個方面進行整理和對比分析.

2.1.1 人工提取特征

所謂人工提取特征是根據人工分析并結合相關的處理手段提取特征的過程.由于閃電哨聲波具有明顯的色散形態，The Stanford VLF Group(2018)，Fiser等(2010)，Oike等(2014)，Ali Ahmad等(2019)和Dharma等(2014)采用圖像處理技術來提取閃電哨聲波的形態特征.The Stanford VLF Group(2018)針對DEMETER衛星的閃電哨聲波，提出的特征提取方法是：對時頻圖(如圖8所示)進行去噪處理和網格劃分，如圖9所示，然后計算每個網格中的平均幅度值，得到如圖10所示的特征表示.

圖8 時頻圖(The Stanford VLF Group, 2018)Fig.8 The time-frequency spectrogram (The Stanford VLF Group, 2018)

圖9 去噪和網格化(The Stanford VLF Group, 2018)Fig.9 Denoising and meshing (The Stanford VLF Group, 2018)

圖10 合并(The Stanford VLF Group, 2018)Fig.10 Binning (The Stanford VLF Group, 2018)

圖11 灰度化的時頻圖(Dharma et al., 2014)Fig.11 Grayscale image of the time-frequency spectrogram (Dharma et al., 2014)

圖12 二值化處理后的時頻圖(Dharma et al., 2014)Fig.12 Binary image of the time-frequency spectrogram (Dharma et al., 2014)

圖13 中值濾波和開運算處理后的時頻圖(Dharma et al., 2014)Fig.13 The results by median filter and open operation (Dharma et al., 2014)

圖14 連通域標記(Dharma et al., 2014)Fig.14 Labeling the connected domain (Dharma et al., 2014)

R=(x0,y0,w,h),

(2)

(3)

其中f(x,y)表示位置(x,y)的像素值，在本實例中其有兩種取值0和1.其目標區域是0，背景區域是1.(5)本方法檢測的是單條譜線的閃電哨聲波.當找到所有的連通區域后，發現不含閃電哨聲波的連通區域具有非常少的連通像素，因此，用下列篩選方法刪除非閃電哨聲波的連通區域.篩選方法是：統計每個連通區域的目標像素的數量，該值小于某閾值時，相應的連通區域和矩形框被刪除，剩下的矩形框標記的連通區域即為檢測到的閃電哨聲波區域.此時，若依然存在兩個重疊的矩形框，則意味著存在兩個連通區域，這兩個區域相互之間不連通且兩個連通區域的重心不同，表示是兩個不同的單條譜線的閃電哨聲波.

該方法適合檢測單條譜線的閃電哨聲波，其缺陷是容易將相近的兩個閃電哨聲波識別為一個閃電哨聲波.

上述閃電哨聲波智能檢測算法的特征提取依賴人工定義.這些方法與傳統的計算機視覺的目標檢測算法類似.在傳統的計算機視覺的目標檢測中，通常人工提取的特征包括尺度不變特征(Scale Invariant Feature Transform)(Lowe, 2004; Cai et al., 2015)，哈爾特征(Viola and Jones, 2001)等，針對提取的特征，再采用AdaBoost(Adaptive Boosting)(Zhang et al., 2016)和DPM(Deformable Parts Model)(Felzenszwalb et al., 2008)等方法對特征進行分類(識別)，其識別結果受特征的影響較大(Pei et al., 2016).因此，基于人工定義特征提取的閃電哨聲波智能檢測算法，最大的局限性是難以獲得魯棒性強的特征，對時頻噪聲十分敏感.

2.1.2 自動提取特征

自動提取特征旨在采用機器學習的方式由計算機自動從數據集中獲取特征.LeCun等(2015)在利用深度神經網絡自動提取圖像特征方面取得了一系列的突破.從此，基于深度卷積神經網絡的目標檢測越來越引起業界的重視(Liu et al., 2018).鑒于深度神經網絡具有強大的自動提取特征的能力，Konan等(2020)提出了兩種基于深度卷積神經網絡的閃電哨聲波的特征提取方法：基于滑動深度卷積神經網絡的閃電哨聲波算法(Sliding Deep Neural Convolutional Neural Network, SDNN)(Konan et al., 2020)和基于YOLOV3 (You Only Look Once and Version 3rd, YOLOV3)目標檢測網絡的閃電哨聲波算法(Konan et al., 2020)，并公開了其源碼(https:∥github.com/Kojey/MSc-whistler-waves-detector).

滑動深度卷積神經網絡的閃電哨聲波特征提取：主要的優勢是其卷積層具有強大的特征學習能力.由于閃電哨聲波具有簡單的形狀，提出了三個卷積層的卷積神經網絡對其進行特征提取.該網絡有兩個主要組成部分：卷積層和分類層，如圖15第1～3層所示為卷積層：第1層和第2層的卷積層具有卷積單元和池化單元：卷積單元負責提取圖像中的局部特征，第一層卷積單元有2個卷積核，第二層的卷積單元有4個卷積核；池化單元用來大幅降低參數量級(降維)防止過擬合；第3層的卷積層有8個卷積核，無池化層.全網絡結構共14個卷積核，主要作用是提取圖像的局部特征.

圖15 滑動的深度卷積神經網絡(Konan et al., 2020)Fig.15 Sliding deep convolutional neural networks (Konan et al., 2020)

圖16為卷積神經網絡模型學習到的卷積核及其相關的特征圖：圖16a是訓練得到卷積神經網絡模型后，顯示了該模型的第一層的兩個卷積核、第二層的四個卷積核以及第三層的8個卷積核；卷積核提取的特征都具有相同的功能(梯度遞減的線)，但在卷積核中的位置不同；此外，卷積核提取的特征具有明顯的特征相似性，其原因是哨聲波的形態特征除了色散狀沒有明顯的其他形狀特征；第三層的卷積核提取的特征的分辨率低，其主要原因是隨著網絡的深入，學習到的特征更加復雜且不容易理解.圖16b顯示的是應用卷積神經網絡模型提取閃電哨聲波特征，可以清楚地觀察到原始哨聲波的形狀.圖16c顯示的是應用卷積神經網絡模型提取非閃電哨聲波的特征，沒有觀察到明顯的形狀.

圖16 卷積層的輸出(Konan et al., 2020)(a) 訓練得到的卷積核; (b) 卷積核提取的閃電形態特征; (c) 卷積核提取的非閃電形態的特征.Fig.16 Output of convolutional layer (Konan et al., 2020)(a) Filters of the trained neural network; (b) Feature maps when detecting whistlers; (c) Feature Maps when detecting noise.

基于YOLOV3深度神經網絡的閃電哨聲波特征提取：借鑒YOLOV3深度神經網絡的特征學習能力，Konan等(2020)直接利用YOLOV3網絡訓練閃電哨聲波數據集，其過程包括：將閃電哨聲波波形數據預處理后轉化為時頻圖，對時頻圖進行圖像標注，在標注的數據上訓練神經網絡YOLOV3.YOLOV3的特征提取是目前最優的特征提取方式，由Redmon等(Redmon et al., 2016; Redmon and Farhadi, 2017, 2018)提出.網絡模型結構圖如圖17所示：此結構主要由75個卷積層構成，卷積層對于分析物體特征最為有效.由于沒有使用全連接層，該網絡可以適應任意大小的輸入圖像.此外，在YOLOV3中沒有出現池化層，取而代之的是將卷積層的步長設為2來達到下采樣的效果，同時將尺度不變特征傳送到下一層.其自動提取特征的優勢體現在YOLOV3中使用了殘差結構(對應著在上面的YOLOV3結構圖中的Residual Block).殘差結構是由一系列卷積層和一條短路徑組成，如圖18所示.隨著網絡越來越深，學習特征的難度也就越來越大，之前的卷積神經網絡是直接學習特征x，加入短路徑后變成在之前學習的特征的基礎上添加某些特征F(x)，以此獲得更為魯棒的特征H(x)=F(x)+x(Redmon and Farhadi, 2018).該策略使得機器學習的目標從學習完整的特征變成學習殘差特征，極大地降低了學習魯棒特征的難度.

圖17 YOLOV3的網絡結構(Redmon and Farhadi, 2018)Fig.17 The network structure of YOLOV3 (Redmon and Farhadi, 2018)

圖18 YOLOV3的殘差網絡(Redmon and Farhadi, 2018)Fig.18 The residual network in the network structure of YOLOV3 (Redmon and Farhadi, 2018)

基于SDNN網絡的閃電哨聲波檢測模型和基于YOLOV3網絡的閃電哨聲波檢測模型的調參過程將在第4.2.2小節進行闡述，模型的訓練過程將在第4.3小節進行介紹.

2.2 分類/識別

分類旨在通過設計算法為提取的特征分配類別標簽，通常也叫識別.閃電哨聲波分類算法的目標是對上節提取的特征進一步判斷其是閃電哨聲波或不是閃電哨聲波.本小節將現有的閃電哨聲波的分類算法分為兩類：基于模板匹配的分類算法和基于神經網絡的分類算法.分別從這兩個方面進行整理和對比分析.

2.2.1 基于模板匹配的分類

由于閃電哨聲波具有明顯的色散形態，The Stanford VLF Group(2018)，Fiser等(2010)，Ali Ahmad等(2019)和Oike等(2014)在提取的閃電哨聲波形態特征的基礎上，采用模板匹配的方式對閃電哨聲波的特征進行分類識別.首先：根據閃電哨聲波的色散形態制作模板，計算待識別區域與模板之間的相似度，當其達到某個閾值時，認為該區域是閃電哨聲波，否則視為沒有閃電哨聲波.

The Stanford VLF Group(2018)制作的閃電哨聲波的模板，是采用格子化的手段對圖像進行降采樣處理，得到如圖10相似的模板.用模板與待檢測的區域進行卷積計算得到相似度度量的值，以此閾值作為分類標準，高于該閾值，判斷待測區域為閃電哨聲波.

Fiser等(2010)提出閃電哨聲波的色散符合Eckersley公式(Gurnett et al., 1990)，t=Df-0.5，其中D是色散率，t是頻率f的波到達的時間.同時指出白天和晚上的色散狀態存在差異，則其制作兩個哨聲波模板.如圖19所示的是白天的哨聲波模板，時間寬度0.4 s，再按照同樣的方法得到晚上的哨聲波模板，模板記為b.然后，從電場波形數據中提取0.4 s的波形數據，做重疊短時傅里葉變換得到頻譜圖，用a表示，接著按照下列公式(4)計算互相關函數得到互相關值.最后采用閾值判斷互相關值最大的地方，即為檢測到的閃電哨聲波，如圖20所示.

(4)

圖19 閃電哨聲波模板(Fiser et al., 2010)Fig.19 Lightning whistler template (Fiser et al., 2010)

圖20 基于互相關熵的閃電哨聲波自動檢測(Fiser et al., 2010)Fig.20 Automatic detection of lightning whistler with cross-correlation entropy (Fiser et al., 2010)

當上述30個頻帶中的6個頻帶與其中一個直線匹配時，被判定為潛在的閃電哨聲波；當同時檢測到多個潛在哨聲波重疊時，從潛在哨聲中選出最匹配的一個，其結果如圖22所示.

圖21 時間變化與頻率之間的關系圖(Oike et al., 2014)Fig.21 Relationship between time and frequency (Oike et al., 2014)

圖22 檢測閃電哨聲波的結果(Oike et al., 2014)Fig.22 Detection results of lightning whistler (Oike et al., 2014)

Ali Ahmad等(2019)認為根據Eckersley公式制作的模板，能夠識別近似直線的閃電哨聲波，無法識別諸如鼻型、長型等類型的哨聲波，則提出了利用曲線檢測和模式規則結合的哨聲波識別方法.再用Bresenham線檢測算法檢測線(Bresenham, 1977)，利用如圖23所示的決策方案識別6種不同類型的閃電哨聲波.

圖23 決策規則(Ali Ahmad et al., 2019)Fig.23 Decision Rules (Ali Ahmad et al., 2019)

2.2.2 基于神經網絡的分類

基于SDNN的分類器(Konan et al., 2020)：如圖24所示，其分類環節是由最后兩層完成，包含128個神經元的全連接層和2個神經元的輸出層.為防止過度擬合，在分類環節增加了2個Dropout層，聯合前面的3個卷積層和2個池化層，整個網絡共有522112個可訓練參數.

圖24 滑動的三層深度卷積神經網絡 (Konan et al., 2020)Fig.24 Sliding three-layer deep convolutional neural network (Konan et al., 2020)

基于YOLOV3深度卷積神經網絡的分類方法：在傳統的深度卷積網絡中利用不同深度的特征直接分類，小的物體會在相對較大的特征圖中被識別出來，而大的物體會在相對較小的特征中被識別出來，從而達到對應不同尺度的物體的目的.然而在實際應用中，不同深度所對應的特征圖包含的信息不同，導致在不同級別的特征中進行目標檢測的精度不夠.為了克服該問題，YOLOV3采用了多尺度融合(Feature Pyramid Networks, FPN)技術，如圖25所示：采用上采樣(Upsample)和融合做法，融合了3個尺度(13×13、26×26和52×52)，在多個尺度的融合特征圖上分別獨立做檢測，使得對小目標的檢測效果得到明顯提升.此外，其輸出層的激活函數用邏輯回歸替換了傳統的softmax，輸出結果由類別標簽變為類別的概率.

圖25 FPN特征融合： (a) 13×13； (b) 26×26； (c) 52×52.Fig.25 FPN fusion

基于SDNN網絡的閃電哨聲波檢測模型和基于YOLOV3網絡的閃電哨聲波檢測模型的調參過程將在第4.2.2小節進行闡述，其訓練過程將在4.3小節進行介紹.

2.3 定位

閃電哨聲波智能檢測算法中的定位方法主要有兩種：滑動窗方法和智能判斷方法.The Stanford VLF Group(2018)，Fiser等(2010)，Ali Ahmad等(2019)和Oike等(2014)提出的定位方法與Konan等(2020)提出的滑動卷積神經網絡的閃電哨聲波檢測的定位方法均采用的是在時頻圖中滑動固定時間長度的窗，然后對窗內的區域進行特征提取和分類，該定位策略主要針對單條頻譜曲線的閃電哨聲波類型，且這些閃電哨聲波圖像區域的寬度相差不大，該方案無法檢測多個閃電源的閃電哨聲波和單個閃電源的多條頻譜曲線構成的閃電哨聲波的類型.通常，多個閃電源的閃電哨聲波類型包含多條頻譜曲線，若多條頻譜曲線均由每個閃電源的第一跳的哨聲波頻譜曲線組成，則每條頻譜曲線的強度都較大，反映在圖像上的特點是多條頻譜曲線的顏色/亮度差別不大；若多條頻譜曲線是由每個閃電源的多跳的哨聲波頻譜曲線組成，則頻譜曲線的頻譜強度存在差異和變化，這種差異會周期性出現，反映在圖像上的特點是顏色/亮度差異隨著時間的變化而周期性出現.而單個事件的多條頻譜曲線的閃電哨聲波類型也包括多條頻譜曲線，但頻譜曲線的頻譜強度隨著時間發生變化，時間滯后的頻譜曲線其頻譜強度逐漸變小，反映在圖像上的特點是其顏色/亮度逐漸變淡.針對頻譜曲線的強度隨時間變化的特點，可采用計算機視覺處理技術提取上述頻譜強度等信息，并進一步探索新方法，使其適用于檢測多閃電源的多頻譜曲線的閃電哨聲波和單閃電源的多頻譜曲線的閃電哨聲波.

基于YOLOV3深度神經網絡閃電哨聲波智能檢測算法 (Konan et al., 2020)采用的定位方法：其數學模型中融合了定位、特征提取和分類.位置的數學模型，如圖26所示.其中tx，ty，tw，th表示模型的預測輸出，cx和cy表示網格單元的坐標，pw和ph表示預測前邊界框的尺寸，bx，by，bw，bh就是預測得到的邊界框的中心坐標和尺寸，坐標的損失函數采用平方誤差.δ(tx)和δ(ty)表示在x軸和y軸上預測得到的邊界框的中心坐標相對于網絡單元坐標的距離.

圖26 YOLOV3 坐標模型 (Redmon and Farhadi, 2018)Fig.26 Coordinate representation in mathematical model of the YOLOV3 (Redmon and Farhadi, 2018)

3 評價指標

評估閃電哨聲波智能檢測算法通常需要用到基準數據，目前采用的基準數據主要是歐洲閃電探測網(http:∥www.euclid.org)和全球閃電定位網絡(http:∥www.wwlln.com/)的數據.在基準數據的基礎上，其評價指標主要有四種：錯分率、虛警率、漏檢率和F1值.四種指標從不同的角度評價了智能檢測模型的性能(楊德賀等, 2018).

設閃電哨聲波是正類樣本(Positive，P)，非閃電哨聲波是負類樣本(Negative，N)，接下來將闡述四種指標的定義.首先，符號定義如表2的混淆矩陣.

表2 矩陣Table 2 Matrix

NTP表示真正例的樣本的數量；NFN表示假負例的樣本的數量；NFP表示假正例的樣本的數量；NTN表示真負例的樣本的數量；四種評價指標的定義如下：

(1)錯分率=(NFP+NFN)/(NTN+NTP+NFP

+NFN) .

(5)

含義：在所有待檢測的數據中，誤識別的數據量所占的比例，通常該指標越小越好.

(2)虛警率=NFP/(NFP+NTP),

(6)

含義：在所有檢測到的閃電哨聲波的數據中，將非閃電哨聲波誤判為閃電哨聲波的數據所占的比例.

(3)漏檢率=NFN/(NTP+NFN),

(7)

含義：在所有已知的閃電哨聲波的樣本中，將閃電哨聲波誤判為非閃電哨聲波所占的比例.

(4)F1值(F1-Score)

=2/(1/Precision+1/Recall),

(8)

其中：Precision表示準確率，準確率=NTP/(NTP+NFP)：其含義是預測正確的圖片個數占總的正類預測個數的比例，該指標越大越好；Recall表示召回率，召回率=NTP/(NTP+NFN)：其含義是確定正類被預測為正類圖片占所有標注圖片的個數，該指標越大越好.通常情況下，準確率高，召回率就低，召回率高，準確率就低.指標F1-Score綜合考慮準確率和召回率的調和值.

目前評價閃電哨聲波智能檢測效果的定量指標如表3所示.通過表3發現，在閃電哨聲波智能檢測性能的評價指標方面，更注重的是錯分率和虛警率.同時，本文認為漏檢率亦重要.

表3 閃電哨聲波檢測算法的性能指標比較Table 3 Performance indexes for lightning whistler detection algorithm

虛警率：基于神經網絡學習到的閃電哨聲波的特征具有抽象和強魯棒的特點，強魯棒通常意味著學習到的特征具有旋轉不變性.比如，將圖27a分別旋轉了90°、180°和270°后則得到圖27(b—d)，再將圖27(a—d)分別輸入到YOLOV3網絡后獲取的特征是相同的.因此，若時頻圖中的非閃電哨聲波的空間現象呈現的形狀特征與閃電哨聲波呈現的形狀特征具有旋轉相似性，如圖28a所示：其矩形框內的圖像不是閃電哨聲波，但其色散狀特征旋轉180°后與閃電哨聲波色散狀極其相似，因此容易被誤識別為閃電哨聲波，致其虛警率升高，同時錯分率增加.

圖27 旋轉不變性(a) 原圖; (b) 逆時針旋轉90°; (c) 逆時針旋轉180°; (d) 逆時針旋轉270°.Fig.27 Rotation invariance(a) Original image; (b) 90° anticlockwise rotation; (c) 180° anticlockwise rotation; (d) 270° anticlockwise rotation.

圖28 實例(a) 虛警實例; (b) 漏檢實例.Fig.28 An example(a) False alarm; (b) Undetected.

漏檢率：若時頻圖中的閃電哨聲波比較密集，如圖28b所示最左邊的矩形框內有兩條閃電哨聲波，但由于密集，導致只識別出一條閃電哨聲波，致使漏檢率增加，同時錯分率增加.

除此之外，錯分率、虛警率和漏檢率具有以下意義：閃電哨聲波的檢測結果通常被用來做進一步的研究.閃電哨聲波智能檢測模型的錯分率綜合評價了采用該模型自動檢測閃電哨聲波的可靠性.錯分率低意味著模型檢測的效果比較可靠，即采用該模型檢測出電磁衛星數據中的閃電哨聲波后，可以直接使用檢測出的結果開展下一步研究.而若錯分率高，則意味著需要人工對識別結果進行復檢，以糾正錯分的數據，但錯分率難以指明復檢的方向.而虛警率和漏檢率恰恰可以指出復檢的方向.若虛警率高則意味著需要人工進一步從檢測到的閃電哨聲波的樣本中進行復檢以剔除將非閃電哨聲波識別為閃電哨聲波的數據.若漏檢率高意味需要人工進一步從原始數據中復檢被遺漏的閃電哨聲波.這些復檢工作對于進一步做閃電哨聲波的研究具有十分重要的作用.

4 現有的閃電哨聲波檢測算法在張衡衛星數據上的效果

閃電哨聲波智能檢測方法通過特征提取、分類和定位等過程實現了自動檢測閃電哨聲波的功能.由于特征提取是閃電哨聲波分類和定位的前提，本節將現有的閃電哨聲波智能檢測算法按照特征的提取方式分為兩類：基于人工定義特征的智能檢測算法和基于神經網絡自動提取特征的智能檢測算法.并從中選擇了三種具有代表性的算法在張衡一號衛星的數據上進行了復現.如圖29所示.在基于人工定義特征的智能檢測算法中，選擇基于形態學處理的智能檢測算法(Dharma et al., 2014)；在基于神經網絡自動提取特征的智能檢測算法中，分別選擇SDNN網絡(Konan et al., 2020)和YOLOV3網絡的智能檢測算法(Konan et al., 2020).接下來，將分別從數據介紹，代表性算法的參數設置，代表性算法的模型訓練和代表性算法的效果四個方面進行闡述.

圖29 三種代表性的算法Fig.29 Three typical algorithms

4.1 數據介紹

利用2018年8月份張衡衛星感應磁力儀VLF頻段Z分量的波形數據.以每2 s窗口且重疊90%的方式進行傅里葉變換生成時頻圖，用這些時頻圖創建數據集，僅關注時頻圖中的單個色散形態的閃電哨聲波.使用深度神經網絡常用的策略是構建較小的超參數訓練集，在其上盡量進行過擬合訓練尋找最佳的參數組合，當超參數選擇確定后再在模型訓練集上訓練模型，最后再在測試集上對該模型進行評價.因此，將張衡衛星的數據分為三類：超參數訓練集含有1000條閃電哨聲波數據(用來選擇模型的超參數)；模型訓練集含有2098條閃電哨聲波數據(超參數選定后用來訓練模型)；模型測試集含25834條閃電哨聲波(用來測試模型的效果)，均為單個頻譜線的閃電哨聲波.

4.2 代表性算法的參數設置

4.2.1 基于人工定義特征的智能檢測算法

本小節主要介紹基于形態學處理的閃電哨聲波檢測算法中每個環節中的參數設置情況.首先將原始時頻圖進行灰度化處理；再采用閾值分割進行二值化處理，該過程中的閾值分割算法采用自適應閾值分割算法(Otsu, 2007)；然后，進行中值濾波處理，該過程中采用的窗口大小是5×5；接著，進行開運算處理，該環節采用的結構算子是5×5的圓形結構算子；再提取連通域，用矩形框對其進行標記，用星號標記其重心位置；最后，統計連通域中的目標像素，當目標像素總量大于某閾值時，保留該連通域的矩形框，否則將刪除該連通區域的矩形框，該環節中的閾值是40.

4.2.2 基于神經網絡自動提取特征的智能檢測算法

目前，神經網絡具有自動地學習魯棒特征的能力且效果最優，但其檢測效果依賴于神經網絡模型的性能，該性能與訓練過程中的超參數的設置有關.其涉及的超參數主要有5種：學習率(learning rate或作lr)，批次大小(batch_size)，激活函數(activation)，輪次(epoch)和優化器(optimizer).在深度神經網絡中，上述5種超參數之間存在比較強的耦合性，使得尋找最優的超參數組合更困難.常用的方法是在超參數訓練集上進行訓練進而獲取最佳的超參數組合，當超參數選擇確定后再在模型訓練集上訓練模型，最后在測試集上對訓練好的模型性能進行評估.本小節將以基于SDNN網絡的閃電哨聲波檢測模型為例介紹超參數調整的過程.

學習率(lr)：指在優化算法中沿損失梯度的反方向更新網絡權重的幅度大小.學習率過大導致模型(在最優點附近)震蕩抑或不收斂；學習率過小則導致模型收斂速度偏慢，需要更長的時間訓練.學習率參數的調整需要結合優化器等參數.

輪次(epoch)：1輪次指的是完整的數據集通過了一次神經網絡的前向傳播和一次反向傳播.輪次數量過少容易導致神經網絡欠擬合，輪次過多容易導致神經網絡過擬合.因此，訓練輪次需要選擇合適的值.

批次大小(batch_size)：指的是每次送入神經網絡的樣本的數量.自2012年以來，其默認的策略是最大限度增大batch_size，直到“吃滿”內存為止.其原因是：batch_size太小導致訓練過程不容易收斂，而大的batch_size使得訓練過程中的振動幅度減小，且容易收斂.但是，2018年，Masters和Luschi (2018)通過開展大量的實驗指出過大的batch_size對測試集不利，即泛化能力較差.并指出最佳的實驗表現都是在batch_size處于2～32之間.因此，其選擇策略是：當有足夠算力時，選取batch_size為32或更小一些；算力不夠時，在效率和泛化性之間做折衷，盡量選擇更小的batch_size(Masters and Luschi, 2018).

優化器：深度神經網絡中的目標函數比較復雜且一般為自定義的損失函數，而神經網絡的優化目標是求得損失結果最小條件下的網絡參數狀態，這種優化問題一般不存在解析解，通常需要使用數值計算的理論和技術更新網絡參數，逼近最優解.將這種迭代更新求解網絡參數的工具稱為優化器.目前深度學習中有多種優化器可供選擇，其中最基本的優化器是隨機梯度下降法(Stochastic gradient descent, SGD)，SGD雖然能達到極小值，但是其他用的時間長，而且可能會被困在鞍點.類似這種的傳統優化器通常將學習率設置為常數或者根據訓練次數調節學習率，從而極大地忽略了學習率其他變化的可能性，進而嚴重影響了模型的性能.而自適應性學習率優化器通過采取一些策略更新學習率提高模型的訓練速度和模型性能，使得其廣泛應用于深度學習中.常用的自適應優化器包括Adadelta(Zeiler, 2012)，RMSprop(Root Mean Square Prop)(Graves et al., 2013)和Adam(Adaptive Moment Estimation)(Kingma and Ba, 2014)等.這些優化器均采用歷史梯度信息調整學習率.

激活函數：在神經網絡中用激活函數給神經網絡加入一些非線性因素，使得網絡可以更好地擬合復雜函數表達式.在深度神經網絡中，常用的激活函數是：Sigmoid，tanh，Relu和Softmax函數.激活函數的選擇通常按照下列規則：Softmax 常常是作為多分類神經網絡的輸出；Sigmoid函數常用于二分類，且適合淺層神經網絡且需要配備較小的初始化權重(防止輸入到激活函數的數值落入飽和區間)；tanh函數具有中心對稱性，適合于有對稱性的二分類；在深度學神經網絡中，ReLU函數被廣泛應用.其主要原因是：ReLU不僅僅使得梯度的計算變快，而且能夠緩解Sigmoid和tanh引起的梯度彌散和爆炸問題.

本文在復現SDNN網絡的閃電哨聲波智能檢測模型時，其調參過程如表4所示，詳細過程如下：

(1)學習率：設初始的參數組合(lr=0.1，optimizer=Adam，batch_size=64，activation=Relu，momentum=0.9，epoch=40).其loss曲線的變化如圖30所示，觀察前10個epoch的loss值在快速下降，從第10個輪次開始，loss值基本不再變化，因此，其最終的loss值是1.087.這說明初始的學習率太大，且從第10個輪次開始神經網絡幾乎停止“學習”.此時，繼續按照下列方式調整學習率：保持其他參數不變，將學習率減少到lr=0.01，得到如圖31所示的loss曲線，其中每個epoch的loss值均低于圖30的loss值，且從第10個epoch開始，loss值具有持續變化的趨勢，最終的loss值是1.051；繼續向下調整lr=0.001，得到如圖32所示的loss曲線：發現整個訓練過程中的loss值不僅僅下降速度快，而且隨著訓練輪次的增加，loss變化的趨勢也比較明顯，說明隨著訓練輪次的增加，神經網絡依然在“學習”，其最終的loss值是0.157.繼續調整lr=0.0001得到如圖33所示的loss曲線，此時發現其loss值下降速度快，神經網絡的”學習“趨勢依然存在，最終的loss值是0.384，將其與圖32進行比較發現：圖33所示的loss曲線整體明顯高于圖32所示的loss曲線，說明圖32所代表的神經網絡的性能優于圖33所代表的神經網絡的性能.此時，停止參數調整，將圖32對應的學習率設為SDNN網絡的初始學習率，即lr=0.001.

圖30 初始學習率0.1，輪次40Fig.30 Learning Rate=0.1 and epoch=40

圖31 初始學習率0.01，輪次40Fig.31 Learning Rate=0.01 and epoch=40

圖32 初始學習率0.001，輪次40Fig.32 Learning Rate=0.001 and epoch=40

(2)輪次：目前最優的超參數是圖32代表的參數組合(lr=0.001，optimizer=Adam，batch_size=64，activation=Relu，momentum=0.9，epoch=40)，其最終的loss值是0.157.繼續觀察圖32所示的loss曲線：發現在epoch=40的地方，loss依舊呈現下降趨勢，此時的參數調整策略如下：繼續增大epoch到70，得到圖34所示的loss曲線，其最終的loss值是0.0603，而且在epoch=70的附近，loss的變化趨勢已經逐漸平穩，說明此時的神經網絡已經停止了學習.因此，參數epoch的值無需再繼續增加，停止對該參數進行調整，設置SDNN網絡的epoch為70.

圖33 初始學習率0.0001，輪次40Fig.33 Learning Rate=0.0001 and epoch=40

圖34 初始學習率0.001，輪次70Fig.34 The learning rate=0.001 and epoch=70

(3)批次大小：在參數組合(lr=0.001，optimizer=Adam，batch_size=64，activation=Relu，momentum=0.9，epoch=70)的基礎上，為了得到最佳的batch_size，嘗試對其進行調整：如表4所示，保持其他參數不變，適當增大batch_size到128，發現其最終的loss值是0.0714；繼續增加batch_size到256，發現最終的loss值變為0.0742；此時意味著大的batch_size導致其loss變差，則嘗試batch_size=32，發現其最終的loss值變為0.0089；持續減少batch_size到16，發現其loss值較0.0089變大，為0.0187，此時，停止參數調整，選擇最低的loss值對應的批次大小為SDNN網絡的batch_size的值，即batch_size=32.

(4)優化器：在參數組合(lr=0.001，optimizer=Adam，batch_size=32，activation=Relu，momentum=0.9，epoch=70)的基礎上，保持其他參數不變，將優化器分別設置為：Adam，Adadelata和RMSprop，計算其最終的loss值，如表4所示，分別是0.0089，0.8011，0.9742.選擇最低的loss值所對應的優化器為SDNN網絡的優化器參數，即optimizer=Adam.

(5)激活函數：在參數組合(lr=0.001，optimizer=Adam，batch_size=32，activation=Relu，momentum=0.9，epoch=70)的基礎上，得到最終的loss值是0.0089.再將激活函數分別調整為tanh，Softmax，Sigmoid后進行模型訓練，得到最終的loss值，如表4所示，分別是1.703，1.847，1.912.選擇最低的loss值(0.0089)所對應的激活函數為SDNN網絡的激活函數，即activation=Relu.

表4 SDNN網絡的調參過程Table 4 Parameter adjustment of the SDNN model

按照類似的過程調試YOLOV3網絡的參數，如表5所示.得到SDNN網絡和YOLOV3網絡的最佳參數如表6所示.這說明按照上述的調參策略得到的超參數組合，能夠使得網絡參數在更新過程中獲得比較合適的方向和步長，達到學習的目的，并最終獲得較優的網絡參數.

表5 YOLOV3網絡的調參過程Table 5 Parameters adjustment of the YOLOV3 model

表6 SDNN網絡和YOLOV3網絡最終的參數設置Table 6 Final parameter settings of SDNN network and YOLOV3 network

4.3 代表性算法的模型訓練過程

基于形態學的閃電哨聲波智能檢測算法的模型，無訓練過程，無訓練時間成本.

基于SDNN網絡的閃電哨聲波智能檢測算法的模型主要有3層卷積神經網絡構成，其超參數的設置如表6所示.該模型的訓練環境：Windows10 64位系統，使用Nvidia GeForce RTX2080Ti顯卡，CUDA10.0和cuDNN調用GPU加速訓練，模型的訓練過程如圖35a所示：在前10個epoch內，其損失值急劇衰減，到第65個epoch后訓練過程基本處于穩定，訓練時間是0.16 h.

圖35 神經網絡的訓練過程(a) SDNN網絡; (b) YOLOV3網絡.Fig.35 Training the neural networks(a) SDNN; (b) YOLOV3.

基于YOLOV3網絡的閃電哨聲波智能檢測算法的模型主要有252層的神經網絡構成，其超參數的設置如表6所示，該模型的訓練環境與SDNN網絡的環境相同，其模型的訓練過程如圖35b所示：在前500epoch內，其損失值急劇衰減，直到第1800個epoch后損失值基本處于穩定，其訓練時間是4.5 h.

其中，神經網絡的訓練時間成本主要受三方面影響：硬件設施，模型復雜度和超參數中的優化器.

在硬件設施相同的情況下，模型越復雜、網絡層數越多，參數越多，其計算量就越大，同時其超參數“輪次”的值越高，因此，導致其訓練的時間越長.比如SDNN網絡由三層卷積網絡構成，其最佳的epoch為70，訓練時間是0.16 h；而YOLOV3神經網絡由252層網絡構成，其最佳的epoch為2000，訓練時間是4.5 h.

在網絡優化的過程中，訓練時間成本主要與優化器的選擇有關.最基本的優化器是SGD，但是SGD的優化速度較慢，容易駐留在鞍點.為了提高訓練速度和性能，大量的實驗表明通過設置一個根據迭代次數衰減的學習率具有明顯效果.其基本思想是：在訓練前期調整學習率的值變大，訓練后期調整其變小.常用的學習率衰減策略包括分段常數衰減，指數衰減，自然指數衰減，余弦衰減等.但對學習率超參數的初始設置具有較強的依賴，而自適應優化器主要利用歷史梯度信息自適應調整學習率，降低了對初始學習率參數的依賴性，同時加快了優化速度.常用的自適應優化器是Adadelta，RMSProp和Adam.在實際應用中，需結合具體問題有針對性地選擇合適的優化器.

4.4 代表性算法的復現效果

本小節主要從定量和定性的角度比較三種復現算法的效果：其中定量比較如表7所示，定性比較如圖36—38所示.基于形態學處理的閃電哨聲波智能檢測算法對噪聲敏感，雖然采用中值濾波、開運算等方式濾除了一部分干擾，如圖36(c,d,e)所示，但仍然難以克服噪聲對該算法產生的影響，如圖36f顯示的是連通域的標記結果，其中每個星號點是一個連通域的重心.對比圖36a和圖36f可知：圖36f中的兩個橢圓區域各有一條閃電哨聲波，即理論上每個橢圓區域內應該各有一個連通區域且每個區域有一個星號點所表示的連通域的重心.仔細觀察圖36f中的兩個橢圓區域發現：每個橢圓區域均存在多個星號點，說明每個橢圓區域存在多個連通域.說明連通特性被干擾點破壞，使得一個完整的哨聲波區域被標記成了若干個連通區域，導致連通域分析出現了較大的偏差，最終使得這兩條閃電哨聲波被漏檢.因此，該方法的漏檢率較高.而滑動卷積神經網絡(SDNN)通過特征學習獲得了較魯棒的哨聲波特征，并有效地克服了上述噪聲的干擾，但由于采用了固定大小的滑動窗，導致圖37最下邊的閃電哨聲波(橢圓標記區域)被漏檢了.YOLOV3網絡由于將自動特征提取、分類和定位統一整合在其損失函數所表示的數學模型中，其在定位方面比滑動卷積神經網絡更精準，如圖38所示，但模型的訓練時間長.總之，人工定義的特征不具有魯棒性容易受到噪聲等因素的干擾，自動化提取的特征較之更為魯棒，且自動定位效果亦更為精準.

表7 三種閃電哨聲波自動檢測算法在張衡一號衛星數據中的效果Table 7 Results of three different methods for detecting lightning whistler from the data captured by ZH-1 satellite

圖36 基于形態學處理的閃電哨聲波檢測算法(a) 原始時頻圖； (b) 灰度化； (c) 二值化； (d) 中值濾波； (e) 開運算； (f) 連通域分析/連通域標記； (g) 連通域篩選； (h) 最終結果.Fig.36 The Process of lightning whistler detection algorithm based on morphological processing(a) Original time-frequency diagram; (b) Gray image; (c) Binarization; (d) Median filtering; (e) Open operation; (f) Connected domain analysis/connected domain labeling; (g) Connected domain filtering; (h) Final result.

圖37 滑動卷積神經網絡的檢測的效果Fig.37 SDNN neural network

圖38 YOLOV3神經網絡的檢測的效果Fig.38 YOLOV3 neural network

5 總結與展望

近年來，電磁衛星的監測數據量急劇增長，閃電哨聲波智能檢測作為人工智能與地球物理交叉的研究方向，受到了廣泛關注.本文對基于電磁衛星數據的閃電哨聲波智能檢測進行了全面而層次化的綜述.通過綜述可以得到以下結論：

(1)電磁衛星：目前電磁衛星的閃電哨聲波的研究主要集中在DEMETER衛星、Arase 衛星和Akebono衛星，其目的是統計閃電哨聲波全球的分布特點，并通過衛星的閃電哨聲波反演電離層的高通電子量、等離子密度等信息.

(2)閃電哨聲波智能檢測方法在時頻圖的基礎上，借鑒計算機視覺和機器學習的相關技術，通過逐步實現特征提取、分類和定位等過程達到了自動檢測閃電哨聲波的目的.表8—10分別從特征提取、分類和定位三個方面對閃電哨聲波智能檢測算法進行了總結整理和對比.

表8總結和對比了閃電哨聲波智能檢測算法中的特征提取技術：①人工提取特征(The Stanford VLF Group, 2018; Fiser et al., 2010; Ali Ahmad et al., 2019; Dharma et al., 2014; Oike et al., 2014)：該方法是依賴人工觀察閃電哨聲波的時頻圖，結合閃電哨聲波的物理機理總結出閃電哨聲波的特征，比如色散形狀等，然后采用數字圖像處理技術設計預處理算法以增強該特征，最后再采用計算機視覺技術設計算法提取該特征.這種方式適合在數據量不大的情況下進行使用，且時間成本低.但閃電哨聲波的色散形狀亦是多種多樣，采用計算機視覺相關技術提取的人工定義的特征僅僅能夠代表一部分閃電哨聲波的特征，該特征不僅僅泛化性差且容易受到噪聲的干擾.說明人工提取的特征魯棒性差，在數據量大的情況下其檢測效果不佳.②自動提取特征(Konan et al., 2020)：首先在時頻圖上，由人工標注大量的單閃電源的單條頻譜線的閃電哨聲波，然后將這些標注圖片數據輸入到神經網絡去訓練得到一模型，該模型能夠自動提取閃電哨聲波特征.這種方式可以獲得較魯棒的閃電哨聲波特征，在數據量增大的情況下，其檢測到的單閃電源的單譜線的閃電哨聲波的效果最佳.

表8 閃電哨聲波智能檢測算法中的哨聲波特征提取方法的總結和對比Table 8 Summary of the feature extracting methods in the lightning whistler detection algorithms

表9總結和對比了閃電哨聲波智能檢測算法中的分類技術：①模板匹配(The Stanford VLF Group, 2018; Fiser et al., 2010; Ali Ahmad et al., 2019; Dharma et al., 2014; Oike et al., 2014)：該方法首先依賴人工提取閃電哨聲波的特征，對其進行歸一化處理后得到閃電哨聲波的特征模板，然后提取待識別的閃電哨聲波的特征，將該特征與特征模板進行匹配，匹配度高意味著該圖像是閃電哨聲波，否則認為該圖像是非閃電哨聲波圖像.該方法存在兩方面的局限性：第一個局限性來自于模板，閃電哨聲波雖然具有色散狀特征，但該特征有大量的不同的形態.因此，難以用一個標準化模板表示所有的閃電哨聲波的色散特征.第二個局限性來自于匹配，匹配策略可被視作一線性分類器，隨著閃電哨聲波數據量的增加以及色散特征的多樣化，線性可分性越來越難以滿足，導致分類結果變差.②神經網絡分類器(Konan et al., 2020)：神經網絡能夠很好地擬合閃電哨聲波圖像與識別結果之間的非線性關系，并克服線性分類器的缺點，總之，神經網絡分類器分類效果優于模板匹配分類器.

表9 閃電哨聲波智能檢測算法中的哨聲波分類(識別)方法的總結和對比Table 9 Summary of the classification methods in the lightning whistler detection algorithms

表10總結和對比了閃電哨聲波智能檢測算法中的定位技術：①滑動窗定位(The Stanford VLF Group, 2018; Fiser et al., 2010; Ali Ahmad et al., 2019; Dharma et al., 2014; Oike et al., 2014; Konan et al., 2020)：在時頻圖上，固定一個大小的時間窗，比如以步長2 s進行滑動，對每一個滑動區域進行特征提取和分類(識別)，該方法簡單易實現但最大的問題是閃電哨聲波持續的時間是不固定的，其時間尺度范圍從毫秒級到秒級，因此，滑動窗定位的方式只能定位出某個固定時間寬度的閃電哨聲波，定位精度低.②智能判斷定位(Konan et al., 2020)：采用YOLOV3神經網絡對閃電哨聲波圖像區域與其相鄰的上下文圖像區域進行建模，經過大量的訓練后得到的神經網絡模型能夠準確地定位到任意時間寬度的閃電哨聲波.

表10 閃電哨聲波智能檢測算法中的哨聲波定位方法的總結和對比Table 10 Summary of the location methods in the lightning whistler detection algorithms

總之，閃電哨聲波智能檢測的研究處于起步階段.目前的閃電哨聲波智能檢測算法僅在識別單條的閃電哨聲波方面取得一定的效果.通過綜述分析目前的電磁衛星的閃電哨聲波智能檢測算法，本文認為該領域未來的發展趨勢和應用前景將主要體現在以下四個方面.

(1)現有的智能檢測算法能夠檢測出大部分的單個閃電源的單條頻譜線的閃電哨聲波類型，如何檢測出表1中由單個閃電源的多條頻譜線組成的閃電哨聲波是探索空間物理現象機制亟需解決的問題.針對多條頻譜線的閃電哨聲波類型，其典型的特點是由多條色散狀的頻譜線組成，衰減速度隨著時間的變化逐漸變小，同時頻譜強度逐漸減少.可在單條頻譜線的哨聲波檢測的基礎上，通過計算機視覺處理技術提取頻譜曲線的斜率變化，用其表達該曲線的衰減信息，將提取得到斜率和頻譜強度隨時間變化的信息作為時間序列，進一步設計循環神經網絡以實現對單閃電源的多條頻譜線的閃電哨聲波進行智能檢測的功能.

(2)單條頻譜線的閃電哨聲波依然存在多種形態，在閃電哨聲波智能檢測結果的基礎上，能夠根據用戶提供的閃電哨聲波圖像自動地從海量的電磁衛星的數據池里找到與其形態相似的若干張哨聲波圖像并返回相關的時空信息，是未來研究的方向之一.其方案如圖39所示：首先通過對電磁衛星的數據進行預處理和傅里葉變換等方式建立電磁衛星數據池；然后，輸入待查詢的閃電哨聲波圖像，提取待查詢的閃電哨聲波圖像的形狀特征、位置信息、頻譜強度特征以及斜率特征等創建待檢索向量；然后借用哈希技術對待檢索向量進行哈希編碼；接下來，在電磁衛星數據池中采用閃電哨聲波智能檢測算法檢測識別單條頻譜線的閃電哨聲波，再對每個檢測到的閃電哨聲波采用計算機視覺等技術提取其形狀特征、位置信息、頻譜強度特征以及斜率特征等信息，利用這些信息創建向量表示該閃電哨聲波，同時對該向量進行哈希編碼.計算待檢索數據的哈希編碼與數據池中每個數據向量的哈希編碼之間的漢明距離，選擇前K(K是非負整數)個最小的漢明距離所對應的閃電哨聲波圖像即為檢索到的閃電哨聲波圖像，且檢索到的圖像共有K張.

圖39 基于哈希的哨聲波圖像智能檢索Fig.39 The image retrieval technology of lightning whistler based on Hashing

(3)目前的閃電哨聲波檢測算法沒有公開的龐大的數據庫，構造閃電哨聲波電磁衛星的數據庫亦是目前亟需解決的問題，可通過收集ZH-1衛星VLF波段感應磁力儀數據和電場探測儀數據，將其轉化成時頻圖后，通過人工處理方式建立統一的電磁衛星的閃電哨聲波數據庫，以期大力推動閃電哨聲波智能檢測算法的發展.

(4)目前的閃電哨聲波智能檢測算法能夠檢測出部分閃電哨聲波，在此基礎上，根據其形態特征開展哨聲波的軌跡跟蹤是進行電子密度估計的關鍵，值得研究和探索.可通過跟蹤其軌跡，利用軌跡的幅值和色散度反演電子密度.

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