基于優(yōu)化MUSIC算法的云閃定位技術(shù)研究

盧 舟， 劉鐘中， 孫秀斌， 賈岸斌

（1.常德市氣象局，湖南 常德 415000；2.成都信息工程大學(xué) 電子工程學(xué)院，成都 610225；3.中國氣象局大氣探測重點(diǎn)實(shí)驗(yàn)室，成都 610225）

0 引言

極端天氣中的閃電活動(dòng)是一種災(zāi)害性天氣現(xiàn)象，以高電壓、大電流、強(qiáng)電磁輻射等為主要物理特征，常常造成嚴(yán)重的人員傷亡和經(jīng)濟(jì)損失［1］。隨著微電子、計(jì)算機(jī)、電力以及航空航天等技術(shù)的不斷發(fā)展，云閃引起的災(zāi)害越來越多，造成的危害也越來越大［2］。云閃過程中輻射的甚高頻段信號(hào)比地閃輻射的低頻段信號(hào)所包含的信息豐富很多，可利用云閃對地閃進(jìn)行預(yù)警。因此，研究精度高實(shí)時(shí)性好的云閃探測系統(tǒng)對云閃放電通道進(jìn)行定位就顯得非常緊要。

由于云閃多發(fā)生在云內(nèi)且具有瞬時(shí)性、劇烈性和頻繁性，要對其進(jìn)行準(zhǔn)確實(shí)時(shí)的監(jiān)測定位仍然有較大難度。近年來，隨著計(jì)算機(jī)技術(shù)、電子技術(shù)和相關(guān)理論的不斷發(fā)展，云閃定位隨之得到快速發(fā)展，許多理論方法在實(shí)際工程應(yīng)用中也逐步得到實(shí)現(xiàn)。根據(jù)張衛(wèi)杰等［3］所著文章中對數(shù)據(jù)采集系統(tǒng)的硬件設(shè)計(jì)可知，用于云閃定位的硬件系統(tǒng)可由云閃VHF信號(hào)接收天線、前端信號(hào)調(diào)理模塊、A/D采樣模塊和FPGA＋DSP數(shù)據(jù)處理模塊等組成。本文采用MUSIC測向算法作為DSP中的解算算法，對其在云閃定位上的應(yīng)用進(jìn)行了分析，提出了優(yōu)化MUSIC算法的方法，并用Matlab仿真了云閃空間譜圖。

1 閃電定位方法

閃電是發(fā)生于云間、云地間或者云體內(nèi)的瞬態(tài)強(qiáng)放電過程，大致可將其分為云閃和地閃兩類。閃電定位主要是通過測量閃電輻射的聲、光、電磁場信息，然后計(jì)算其相關(guān)參數(shù)來確定閃電的空間位置，進(jìn)而推測其活動(dòng)趨勢［4］。云閃定位和地閃定位原理相近，但前者是利用云閃電磁脈沖輻射源發(fā)出的甚高頻信號(hào)，后者則是利用閃電輻射的低頻和甚低頻電磁波。目前，閃電定位方法主要包括定向法、時(shí)差定位法以及時(shí)差測向聯(lián)合法，總體趨于向時(shí)差測向聯(lián)合法（TDOA）發(fā)展，該方法是通過多個(gè)子站既記錄閃電輻射源信號(hào)到達(dá)時(shí)間差又解算其來波信號(hào)方向信息來確定閃電空間位置。采用時(shí)差測向聯(lián)合法的多站云閃定位系統(tǒng)空間分布圖如圖1所示：

圖1 多站云閃定位系統(tǒng)空間分布圖

在多站云閃定位系統(tǒng)中，各探測子站的定位誤差主要來源于觀測時(shí)間誤差、布站位置誤差以及算法誤差。本文采用MUSIC算法作為時(shí)差測向聯(lián)合法中的測向算法來測算云閃輻射源信號(hào)的方位角和仰角，具有很高的分辨率和實(shí)時(shí)性，可以提高測向精度、減少定位誤差。

2 基于均勻圓陣的二維MUSIC算法

2.1 均勻圓陣模型

均勻圓陣（UCA）是一種具有全向探測能力的陣列，它能同時(shí)獲取多個(gè)云閃輻射源的方位角和俯仰角，從而對云閃進(jìn)行定位。它由N個(gè)相同的全向陣元組成，按等間隔分布在半徑為R的圓周上，圓周中心O為球坐標(biāo)系的原點(diǎn)和云閃入射信號(hào)的相位參考點(diǎn)，如圖 2 所示［5］［6］。 圖中，θi和 φi分別為云閃輻射源入射信號(hào)的仰角和方位角，i＝1，2，3…。若有P個(gè)遠(yuǎn)場窄帶信號(hào)源S（t）入射到此均勻圓陣中，則第 i個(gè)信源的入射角度為（θi，φi）。

圖2 探測系統(tǒng)中的均勻圓陣模型

方向矢量 a（θi，φi）是入射信號(hào)（θi，φi）對應(yīng)的陣列響應(yīng)，其表達(dá)式為［7］：

上式中：γm＝2πn/N（n＝0，1，…，N-1）；R 為圓陣半徑。

事實(shí)上，均勻圓陣所接收到的信號(hào)和均勻線陣所接收到的信號(hào)是相同的，但是二者對比起來均勻圓陣的陣因子 A＝［a （θ1，φ1） a （θ2，φ2.） … a（θP，φP）］中有入射信號(hào)的方位角和仰角兩個(gè)因子，比線陣的陣因子包含的方位信息多。采用MUSIC算法實(shí)現(xiàn)測向?qū)﹃嚵薪Y(jié)構(gòu)沒有嚴(yán)格要求，因此可以直接使用均勻圓陣對來波的方位角和仰角進(jìn)行測量。

2.2 二維MUSIC算法原理

MUSIC算法是通過對均勻圓陣天線所接收到的云閃輻射源信號(hào)的協(xié)方差矩陣進(jìn)行特征分解，然后由信號(hào)特征值、噪聲特征值分別得到信號(hào)子空間VS和噪聲子空間VN，再利用它們的正交特性對譜函數(shù)進(jìn)行譜峰搜索，最終完成對云閃信號(hào)的方位估計(jì)［8］。

假設(shè)P個(gè)窄帶信號(hào)入射到均勻圓陣天線中，陣元數(shù)為M，陣元間夾角為固定值，而且陣元數(shù)大于此時(shí)的信號(hào)數(shù)。若這些信號(hào)的波達(dá)方向分別為（θ1，φ1），（θ2，φ2），...，（θP，φP），陣因子 A＝［a（θ1，φ1），a（θ2，φ2），...，a（θP，φP）］M×P，n（t）為噪聲矢量，則陣列輸出信號(hào)可用下式表示［9］：

相應(yīng)的輸入陣列協(xié)方差矩陣為：

其中，RSS＝E［S（t）SH（t）］是信號(hào)相關(guān)矩陣，σn2是噪聲方差。

對2-3式中的RX進(jìn)行特征值分解得到M個(gè)特征值，對應(yīng) M 個(gè)特征向量 VK（K＝1，2，…，M）。 由較大的P個(gè)特征值對應(yīng)信號(hào)子空間，較小的M-P個(gè)特征值對應(yīng)噪聲子空間。把噪聲特征向量作為列構(gòu)造成噪聲子空間VN：

對應(yīng)于信號(hào)分量的導(dǎo)向矢量空間與VN正交，可得到空間譜搜索函數(shù)Pmusic（θ）如下式：

通過上式可知，分母中的 VN和 a（θ，φ）正交，所以Pmusic（θ，φ）的值在云閃信源處會(huì)有一個(gè)很尖的譜峰，此處即為云閃信源的位置。當(dāng)云閃輻射出多個(gè)信號(hào)時(shí)，可以通過搜索空間譜中的波峰來估計(jì)各信號(hào)的來波角度信息。二維MUSIC算法處理流程圖如圖3所示：

圖3 MUSIC算法簡要流程

3 算法優(yōu)化

在云閃定位系統(tǒng)中采用MUSIC算法來解算云閃輻射源的方位角和仰角等閃電參數(shù)具有很大的優(yōu)勢，該方法具有實(shí)時(shí)性好、測向精度高、能多信號(hào)同時(shí)測向等優(yōu)點(diǎn)。采用MUSIC算法實(shí)現(xiàn)云閃定位主要包括空間譜函數(shù)構(gòu)造、信源數(shù)估計(jì)以及譜峰搜索三個(gè)步驟，下面詳細(xì)介紹對這三個(gè)部分的優(yōu)化。

3.1 空間譜函數(shù)優(yōu)化

對均勻圓陣輸出信號(hào)的協(xié)方差矩陣進(jìn)行優(yōu)化處理，充分利用求得的特征值對應(yīng)的信號(hào)子空間和噪聲子空間包含的信息。

優(yōu)化空間譜函數(shù)構(gòu)造方法的實(shí)現(xiàn)步驟如下：

（1）采集云閃輸入信號(hào)樣本 X（i），i＝1，2，…N，計(jì)算輸入?yún)f(xié)方差矩陣Rx，如公式（2-3）所示，構(gòu)造新的協(xié)方差矩陣RX如式（3-1）：

上式中I為M×M的反向單位矩陣，R^*X為R^X的共軛矩陣。

（2）對3-1式中新構(gòu)造的協(xié)方差矩陣RX進(jìn)行特征分解重構(gòu)矩陣得：RX＝UΛH，可得對應(yīng)的信號(hào)子空間RXS1和噪聲子空間RXN1如下式：

（3）用低秩矩陣代替滿秩矩陣RX，再次構(gòu)造新的協(xié)方差矩陣RXX并對其進(jìn)行特征值分解，重復(fù)步驟（2），再一次構(gòu)造新的信號(hào)子空間和噪聲子空間［10］：

（4）對兩次得到的信號(hào)子空間和噪聲子空間進(jìn)行平均，如式（9）：

（5）根據(jù)US和UN構(gòu)造新的空間譜函數(shù)，如下式：

從上式可知，該空間譜函數(shù)充分利用了信號(hào)子空間和噪聲子空間所包含的信息，提高了MUSIC算法的穩(wěn)定性和在非理想情況下的分辨性能。根據(jù)上式進(jìn)行譜峰搜索，可以得到準(zhǔn)確性高、精度好的入射云閃輻射源信號(hào)方位角θ和仰角φ。

3.2 信號(hào)源數(shù)估計(jì)優(yōu)化

在MUSIC算法中估計(jì)信號(hào)源的來波方向之前需要先估計(jì)出信源數(shù)，否則會(huì)出現(xiàn)虛警或漏警。一般選用基于信息論準(zhǔn)則的信源數(shù)估計(jì)方法來完成該步驟，該方法在有高斯白噪聲的情況下能獲得較好的估計(jì)性能。但是在實(shí)際的云閃輻射源探測環(huán)境中，天線陣接收到的噪聲往往是空間色噪聲，這樣便不能準(zhǔn)確獲取陣列協(xié)方差矩陣的噪聲特征值，從而造成該方法的估計(jì)性能較差。為了解決這個(gè)問題，采用基于總體最小二乘擬合法來估計(jì)云閃信號(hào)源數(shù)。

基于總體最小二乘法（TLS）的信源數(shù)估計(jì)方法步驟如下：

（1） 算出均勻圓陣接收信號(hào)的輸出協(xié)方差矩陣：

（2）根據(jù)MUSIC算法原理可知，對協(xié)方差矩陣進(jìn)行特征值分解可以得到M個(gè)特征值，將其按從大到小的順序排列如下［11］：

（3）對式3-7中后面的k個(gè)特征值進(jìn)行擬合，擬合點(diǎn)數(shù)i（2≤i≤M）從所得特征值個(gè)數(shù)M到2逐漸減小，依次可以獲得對應(yīng)的擬合偏差 ε（1），ε（2），…，ε（M-1）。

（4）假設(shè)云閃入射信號(hào)源數(shù)為 n（n＝1，2，…，M），因?yàn)棣牛╧）為參與擬合數(shù)目的倒數(shù)擬合，所以待估計(jì)的信號(hào)源數(shù)目為。按式＝argmin{ε（k）＋p（k）}構(gòu)造目標(biāo)估計(jì)函數(shù)，式中罰函數(shù) p（k）＝1/2k（2M-k）logL，L為快拍數(shù)，能夠使該式取得最小值的k值就是所要估計(jì)的云閃信號(hào)源數(shù)。

3.3 譜峰搜索優(yōu)化

由于譜峰搜索在MUSIC算法中耗時(shí)較多，影響了云閃定位的實(shí)時(shí)性，所以在保證一定測向精度的情況下，提出了一種精度高實(shí)時(shí)性好的搜索方法。（1）采用一個(gè)較大的角度搜索步長（如 2°），在［0°，180°］的范圍里盡可能尋找出所有極大值出現(xiàn)的大概位置。假設(shè)粗略搜索到M個(gè)極大值，在沒有噪聲和干擾信號(hào)的情況下，M值等于信源數(shù)P。但在實(shí)際探測環(huán)境中存在無法避免的噪聲以及干擾信號(hào)，所以會(huì)有多于信源數(shù)P的極大值出現(xiàn)。此時(shí)需要對M和信號(hào)源估計(jì)數(shù)P的大小進(jìn)行比較，若M＝P則直接進(jìn)入細(xì)搜，否則，將所有峰值依次從大到小排序，選出其中P個(gè)最大的極大值，這些值就對應(yīng)著信號(hào)到達(dá)方向。 （2）在每個(gè)經(jīng)粗選的極大值（i，j）附近采用較小的角度搜索步長（如0.1°）進(jìn)行搜索，并分別與周圍八個(gè)點(diǎn)（i-1，j）、（i＋1，j）、（i，j-1）、（i，j＋1）、（i＋1，j＋1）、（i＋1，j-1）/（i-1，j-1）、（i-1，j＋1）比較；如果該極大值大于周圍八個(gè)點(diǎn)，則其為有效譜峰，（i，j）即為來波信號(hào)的方位角和仰角。（3）如果譜峰搜索結(jié)果還未達(dá)到精度要求，那么返回（2）粗略搜索出的峰值附近繼續(xù)減小角度搜索步長，重復(fù)進(jìn)行搜索，直到結(jié)果滿足云閃定位的精度要求［12］。

采用這種一次粗搜全方位多次細(xì)搜的優(yōu)化方法，可以避免在要求較高精度的爬山搜索算法中不必要的計(jì)算量，并且可以快速確定云閃信號(hào)來波方向的大概角度，然后再通過多次細(xì)搜找到滿足精度要求的方位角、仰角。優(yōu)化后的譜峰搜索方法可以極大的減少搜索運(yùn)算時(shí)間，同時(shí)提高了云閃定位結(jié)果的估計(jì)精度。搜索算法流程圖如圖4所示：

圖4 譜峰搜索流程圖

4 算法仿真及性能分析

假設(shè)均勻圓陣半徑滿足 r≤1/［4sin（π/M）］，陣元個(gè)數(shù) M＝6，快拍數(shù)：1024，信號(hào)源數(shù) P＝3；信號(hào)的方位角分別為（80°、100°、120°），俯仰角分別為（40°、50°、60°）， 分別在信噪比取 20dB、30dB 的情況下進(jìn)行二維MUSIC譜估計(jì)，如圖5（見彩頁）、圖6（見彩頁）所示：

假設(shè)均勻圓陣半徑滿足 r≤1/［4sin（π/M）］，陣元個(gè)數(shù) M＝8，快拍數(shù)：1024，信號(hào)源數(shù) P＝3；信號(hào)的方位角分別為（80°、100°、120°），俯仰角分別為（35°、40°、45°），信噪比均為 30dB，二維 MUSIC 算法優(yōu)化前和優(yōu)化后的二維MUSIC算法譜估計(jì)如圖7、圖8所示：

通過多次二維分級(jí)譜峰搜索獲得云閃輻射源的估計(jì)方位角、仰角平均值，如表一所示。

表1 兩種情況下二維估計(jì)仿真結(jié)果

通過對比圖六和圖七的空間譜圖可以看出，如果其他條件相同，增大信噪比可以使基于均勻圓陣的二維MUSIC算法DOA估計(jì)性能得到提高，對三個(gè)信號(hào)的分辨能力也得到增強(qiáng)；圖7（見彩頁）和圖8（見彩頁）的仿真條件考慮到云閃發(fā)生的連續(xù)性，針對信源之間的角度間隔可能會(huì)比較小進(jìn)行了設(shè)置；將優(yōu)化前和優(yōu)化后的二維MUSIC算法進(jìn)行仿真對比發(fā)現(xiàn)，優(yōu)化后的二維MUSIC算法譜峰更加陡峭，表明對來波方向估計(jì)的分辨力得到進(jìn)一步提升。從表一中的方位角和仰角估計(jì)數(shù)據(jù)可知，優(yōu)化后的二維MUSIC算法的DOA估計(jì)精度得到了很大的提高，并且對多個(gè)相隔較近的云閃輻射源信號(hào)能夠同時(shí)測向。

5 結(jié)論與展望

本文簡要介紹了閃電定位基本方法和均勻圓陣建模模型，重點(diǎn)概述了MUSIC算法中各個(gè)重要步驟的優(yōu)化方法；對算法中構(gòu)造空間譜函數(shù)的方法進(jìn)行了優(yōu)化，充分利用了信號(hào)和噪聲所包含的信息，提高了算法的分辨率；針對實(shí)際的空間噪聲環(huán)境，提出采用基于總體最小二乘擬合法來估計(jì)云閃信號(hào)源數(shù)；針對MUSIC算法中譜峰搜索耗時(shí)較多影響定位實(shí)時(shí)性的問題，提出了一種精度高實(shí)時(shí)性好的譜峰搜索方法。最后給出了不同信噪比情況下基于二維MUSIC算法的空間譜估計(jì)結(jié)果，對優(yōu)化前和優(yōu)化后的算法性能進(jìn)行了分析。通過對比分析，可以看到基于均勻圓陣二維優(yōu)化MUSIC算法作為云閃定位系統(tǒng)測向算法的優(yōu)越性，使用這種優(yōu)化的測向方法可以有效提高測向結(jié)果的實(shí)時(shí)性和精度。將這種方法應(yīng)用于估計(jì)云閃輻射源的來波方向，再結(jié)合多片并行高速數(shù)字信號(hào)處理器，可以實(shí)時(shí)解算多個(gè)云閃信號(hào)的方位角和仰角，準(zhǔn)確獲取云閃的空間位置，為云閃定位系統(tǒng)的實(shí)現(xiàn)提供了一種更優(yōu)的方法。

參考文獻(xiàn)：

［1］李云敏，孫秀斌，楊崧令，等.云閃定位算法及誤差分析［J］.電子測量與儀器學(xué)報(bào)，2012，26（10）：917-922.

［2］張義軍，周秀驥.雷電研究的回顧和進(jìn)展［J］.應(yīng)用氣象學(xué)報(bào)，2006，17（6）：829-834.

［3］張衛(wèi)杰，候孝民.高速大容量數(shù)據(jù)采集系統(tǒng)設(shè)計(jì)與實(shí)現(xiàn)［J］.電子測量與儀器學(xué)報(bào)，2005，19（4）：51-55.

［4］張義軍，孟青，馬明，等.閃電探測技術(shù)發(fā)展和資料應(yīng)用［J］.應(yīng)用氣象學(xué)報(bào)，2006，17（5）：611-620.

［5］ Mathews C P， ZoltowskiM D.Eigenstructure Techniques for 2-D Angle Estimation with Uniform Circular Array［J］.IEEE trans on Sp r，1994，42（9）：2395-2404.

［6］ 宋士瓊，陳旗，胡生亮.基于九元均勻圓陣的MUSIC測向方法研究［J］.計(jì)算機(jī)仿真，2008，25（3）：123-125.

［7］ 張小飛，汪飛，陳偉華.陣列信號(hào)處理的理論與應(yīng)用（第2 版）［M］.北京，國防工業(yè)出版社，2013.

［8］王永良，陳 輝，彭應(yīng)寧，等.空間譜估計(jì)理論與算法［M］.北京：清華大學(xué)出版社，2005：82-139.

［9］ SchmidtR.O.Multiple emitterlocation and signal parameterestimation ［J］.IEEE Trans.OnAP， 1986，34（3）：276-280.

［10］趙謙，董民，梁文娟.DOA估計(jì)算法的一種修正MUSIC算法的研究［J］.計(jì)算機(jī)工程與應(yīng)用，2012，48（10）：102-105.

［11］藍(lán)曉宇.提高空間譜估計(jì)分辨率的超分辨測向算法研究［D］.哈爾濱，哈爾濱工程大學(xué)，2012.

［12］郭躍，王宏遠(yuǎn)，陳思捷，等.陣列測向中的精確高速并行譜峰搜索算法 ［J］. 微電子學(xué)與計(jì)算機(jī)，2007，24（12）：50-54.

［13 潘軍，鐘一帆，蒙劍.雷電技術(shù)在現(xiàn)代建筑中的應(yīng)用［J］.氣象研究與應(yīng)用，2008，29（3）：49-54.

［14］鄒興奮，劉志輝，李可娟.建筑物防雷接地裝置結(jié)構(gòu)探析［J］.氣象研究與應(yīng)用，2010，31（S2）.

［15］黃劍釗.新型國家自動(dòng)站現(xiàn)場總線CAN的研究及維護(hù)方法［J］.氣象研究與應(yīng)用，2016，37（2）：87-89.

［16］張娟，林卓宏，陳巧淑，等.自動(dòng)氣象站觀測場防雷接地制式的技術(shù)分析［J］.氣象研究與應(yīng)用，2012，33（4）69-71.

［17］姚家釗，陳華宣.基層氣象臺(tái)站防雷問題淺析及解決方法［J］.氣象研究與應(yīng)用，2009，30（1）：78-81.

［18］黃仁立，羅曉軍，周開春.防城港天氣雷達(dá)樓雷擊防護(hù)等級(jí)及防護(hù)設(shè)計(jì)［J］.氣象研究與應(yīng)用，2014，35（4）：111-114.

［19］黎錦雷，韋菊，楊玉靜.新型自動(dòng)氣象站故障分析與排除［J］.氣象研究與應(yīng)用，2015，36（4）：102-104.

［20］傅俊霖，黃君健，何肖珍.防雷裝置接地電阻的認(rèn)識(shí)和探討［J］.氣象研究與應(yīng)用，2008，29（2）：68-69.

［21］韓建海，吳松.一次風(fēng)電場雷災(zāi)的特點(diǎn)及原因分析［J］.氣象研究與應(yīng)用，2016，37（1）：121-123.

［22］羅天龍.電涌保護(hù)器的簡易檢測方法［J］.氣象研究與應(yīng)用，2013，34（2）：92-94.

［23］楊仲江，盧燕，曹書華.用于防雷工程的電涌保護(hù)器的測試研究［J］.氣象研究與應(yīng)用，2007，28（4）：54-57.

［24］ 譚惠冰，杜建德，梁偉漢，等.電涌保護(hù)器（SPD）在低壓電氣系統(tǒng)中的設(shè)計(jì)［J］.氣象研究與應(yīng)用，2013，34（4）：92-97.

［25］林政，黎梓華，唐雷.淺談如何利用法拉第籠原理防護(hù)雷電電磁脈沖［J］.氣象研究與應(yīng)用，2009，30（1）：83-84.