駱瑤瑩， 卞宏志， 劉全楨， 劉保全， 傅正財， 張建勛， 劉亞坤

(1. 上海交通大學 電力傳輸與功率變換控制教育部重點實驗室， 上海 200030；2. 中國石化青島安全工程研究院 化學品安全控制國家重點實驗室， 山東 青島 266071；3. 國網(wǎng)福建省電力有限公司建設(shè)分公司， 福州 350001)

雷電作為一種長間隙強放電現(xiàn)象，產(chǎn)生時伴有強烈的電磁和聲波信號[1]，嚴重影響油罐儲運、電力設(shè)備、航空航天活動和森林安全等[2-4].雷電定位技術(shù)基于雷電產(chǎn)生時的電、磁、聲和光等信號特征，通過分析不同信號的傳播特性，反演計算得到雷電的位置信息和產(chǎn)生時間[5].雷電定位技術(shù)是目前用于研究雷電活動最主要的手段之一[6]，為關(guān)鍵系統(tǒng)的雷電防護和雷電災(zāi)害分析提供數(shù)據(jù)保障[7].

現(xiàn)有雷電定位系統(tǒng)多測量雷電放電時發(fā)出的電磁脈沖信號，其測量探頭可工作在甚高頻、甚低頻/低頻或極低頻等頻段，依據(jù)多個探頭測量得到的雷電電磁脈沖的波峰到達時間和各探頭間的時間差，計算確定雷電電磁輻射源的位置和產(chǎn)生時間.由該方法測得的雷電信息主要受探測子站探頭的時鐘精度、電磁環(huán)境和反演定位算法等影響.在雷電定位算法研究中，趙文光等[8]綜合利用到達時差定位法(Time Difference of Arrival, TDOA)和方向定位法，通過考慮平差計算剔除粗差的算法，提高了雷電定位的精度.胡志祥等[9]在雷電定位計算中引入粒子群優(yōu)化算法，得到了更精確的雷電發(fā)生位置.郭小紅等[10]提出優(yōu)化果蠅算法，實現(xiàn)了雷電定位結(jié)果的快速計算.李濤等[11]利用改進密度聚類算法對雷電位置進行聚類解算，提高了雷電定位的抗誤差干擾能力.梁麗等[12]綜合密度聚類算法和網(wǎng)格搜索法，降低了電磁探測數(shù)據(jù)誤差等對雷電定位的影響.在實際的TDOA雷電定位系統(tǒng)中，探測子站存在測量誤差，因此采用3個及3個以上探測子站的雷電定位系統(tǒng)會出現(xiàn)雷電位置信息難以求解或求解不唯一的問題.如何解決雷電定位計算中的求解問題，且減小探測站與雷電相對位置以及定位網(wǎng)絡(luò)幾何尺寸對計算結(jié)果的影響，具有實際應(yīng)用意義.

安裝于我國油罐等系統(tǒng)中的雷電定位系統(tǒng)具備同時采集聲波和電磁信息的能力.為充分利用聲波和電磁信號, 李皖等[13]利用人工識別雷擊電磁脈沖信號法得到了電磁信號與雷聲信號的時間差；湯偉等[14]將雷聲按點聲源處理，利用坐標搜索法進行雷聲方向定位.綜合已有研究發(fā)現(xiàn)，利用雷電放電時的聲波和電磁信號同時進行雷電定位可以提高定位精度，但在雷電電磁和聲波信號的融合交互方面還需進一步研究.因此，本文提出一種基于Levenberg-Marquardt算法(簡稱L-M算法)和聲波信息融合的TDOA雷電定位改進方法.根據(jù)已有探測站點的布站信息，將目標雷電監(jiān)測區(qū)域劃分為16個子區(qū)域；利用電磁信號的到達時差，在整個監(jiān)測區(qū)域內(nèi)計算雷電初定位結(jié)果，并將其作為依據(jù)，動態(tài)納入雷聲站點測得的聲波信息，避免聲波信號在長距離傳播時的大幅衰減和畸變.在計算求解雷電發(fā)生位置時，利用可融合聲波信息分析的L-M迭代算法，得到更精確的最優(yōu)解.所提方法在TDOA電磁信息中納入聲波信息，可以幫助解決雷電定位計算中的求解問題，提升雷電定位系統(tǒng)的定位精度和抗誤差擾動能力.

1 TDOA雷電定位方法

基于TDOA雷電定位系統(tǒng)，利用3個及3個以上探測子站測得的電磁信息，計算雷電電磁脈沖到達各探測子站的時間差，以確定電磁輻射源的位置.TDOA雷電定位系統(tǒng)一般由連接在同一個通信網(wǎng)絡(luò)上的若干探測站、中心站和用戶終端組成[15].其中，探測站主要包括電磁脈沖接收天線、信號處理單元、高精度時鐘單元、電源和通信等模塊.伴隨雷電產(chǎn)生的電磁信號到達接收天線，并由信號處理單元處理獲得同步的時間標簽后，通信口將電磁脈沖波峰到達探測站的精確時間實時發(fā)送到中心站.中心站的雷電位置分析儀根據(jù)電磁脈沖到達各探頭的時間差來確定雷電電磁輻射源的位置，并由中心站將定位結(jié)果輸出到各用戶工作站.在保證系統(tǒng)運行可靠性的基礎(chǔ)上，多個定位系統(tǒng)聯(lián)網(wǎng)可以實現(xiàn)數(shù)據(jù)共享，滿足雷電定位系統(tǒng)的擴展需求[16].

在雷電電磁信號測量中，雷電輻射源產(chǎn)生的電磁波到達探測站#1和探測站#2的時間差為Δt12，結(jié)合兩個探測站的地理位置信息，可以確定一條雙曲線，其計算式為

(1)

式中：(x1,y1)和(x2,y2)分別為TDOA探測子站T1和T2的位置坐標；c為電磁波傳播速度，且c=3×108m/s.根據(jù)雙曲線交匯原理，電磁波輻射源可以由交匯點位置坐標確定.在3個及3個以上探測站的情況中，不同雙曲線相互交匯并產(chǎn)生公共交點，據(jù)此可以確定雷電發(fā)生位置.

在TDOA雷電定位的實際計算中，雙曲線發(fā)散的幾何形狀特性導致整個雷電監(jiān)測區(qū)域的定位精度不均勻[17].采用3站TDOA定位可能會出現(xiàn)兩條雙曲線交于兩點的情況，導致雷電定位系統(tǒng)作出偽雷擊位置判斷[18].對于3站以上的TDOA定位系統(tǒng)，電磁波探測站在測量時存在誤差，因此利用探測站之間的時間差信息得到的多條雙曲線不易交于一點，從而影響基于到達時差信息的雷電定位方程組求解.對此，通常利用優(yōu)化算法來減小測量誤差的影響并解決方程組的求解問題[8-12].

2 基于L-M算法和聲波信息融合的TDOA雷電定位改進方法

2.1 L-M迭代算法

L-M算法是牛頓法的一種改進.在計算時，其能夠避免牛頓法由于Hessian矩陣奇異而導致算法無法繼續(xù)迭代的情況[19]，常用于非線性最小二乘問題的最優(yōu)化實現(xiàn)，能夠有效解決時差信息冗余時雷電位置最優(yōu)解的求解問題.

L-M算法用于雷電定位計算時的流程如圖1所示.其中，L0為雷電初定位結(jié)果，即L-M算法的迭代初值.在求解由多個探測子站的聲波信息得到的非線性方程組f(x)=0時，可轉(zhuǎn)換為最小值尋找問題：

圖1 利用L-M迭代算法進行雷電定位計算的流程圖Fig.1 Flow chart of lightning location calculation using L-M iterative algorithm

(2)

式中：x為需求解的雷電定位結(jié)果.高斯-牛頓法將迭代步長設(shè)定為

(3)

式中：I為單位矩陣；μ>0保證系數(shù)矩陣正定，從而保證迭代的下降方向，克服了因Hk不可逆而導致的迭代無法繼續(xù)問題.為保證算法在接近最優(yōu)解時能夠快速收斂[20]，當μ較大時，

當μ較小時，

hlm≈hgn

其中：τ為限定μ0在適當數(shù)值的系數(shù).μ的數(shù)值更新由系數(shù)?決定：

(4)

式中：F(xk)-F(xk+hlm)表示目標函數(shù)在步長hlm下的實際變化；L(0)-L(hlm)表示目標函數(shù)的預估變化值，其計算式為

(5)

(6)

‖xnew-x‖<ε2(‖x‖+ε2)

(7)

k>kmax

(8)

式中：ε1和ε2為迭代終止的條件參數(shù)；kmax為最大迭代步數(shù)，合適的kmax值可以防止迭代無限循環(huán)引起的算法崩潰.ε1、ε2和kmax主要根據(jù)L-M算法的具體應(yīng)用情況進行選取，對算法最終收斂結(jié)果影響不大.在本算法中，設(shè)ε1=ε2=10-12，kmax=200.在計算過程中，3個參數(shù)的取值滿足上述條件之一即終止算法，輸出雷電定位計算結(jié)果.當探測站測得的信號存在誤差或信息冗余時，會導致雷電定位計算無法得到非線性方程的唯一解，而L-M迭代算法可以為尋找雷電位置最優(yōu)解提供可靠的算法支持.

2.2 融合聲波信息的TDOA定位

TDOA雷電定位系統(tǒng)對測試精度要求較高，其中，探測站布局方式會引起時間誤差[21]，雷電電磁波峰值點因傳播路徑和距離的不同而發(fā)生漂移或畸變，導致時間測量誤差[22]，使得TDOA的定位誤差幅值達到幾百米或幾千米.此外，整個雷電監(jiān)測目標區(qū)域的TDOA雷電定位精度不均勻，存在局部定位誤差較大的區(qū)域塊.雷電發(fā)生時會發(fā)出強烈的電磁信號和聲波信號.聲波信號的傳播特性與信號頻率有關(guān)，在相同距離和環(huán)境背景下，信號頻率越低則聲波傳播性能越優(yōu)異[23].隨著傳播距離的增大，低于 100 Hz 的聲波信號幅度衰減速度較慢，可以作為有效的雷聲信息進行雷電定位研究.因此，可以考慮在進行雷電定位時納入雷聲聲波信息來減小局部區(qū)域塊的定位誤差，從而提升整個雷電監(jiān)測目標區(qū)域的雷電定位效果.

在TDOA雷電定位系統(tǒng)融合聲波信息時，利用波形延時壓縮處理技術(shù)可以有效發(fā)送聲波信號，減小聲波信號傳遞帶來的影響[14].為提高所納入雷聲聲波信息的準確性，綜合根據(jù)已有TDOA探測子站和雷聲探測站的位置信息，將整個雷電監(jiān)測區(qū)域進行子區(qū)域劃分處理，如圖2所示.其中，T3和T4為TDOA探測子站，N1～N4為4個雷聲探測站.共劃分為16個雷電監(jiān)測子區(qū)域.其中，子區(qū)域#1、#4、#13和#16為近時差站子區(qū)域，子區(qū)域#6、#7、#10 和#11為近雷聲站子區(qū)域，子區(qū)域#2、#3、#5、#8、#9、#12、#14和#15為時差-雷聲站點近似等距子區(qū)域.以TDOA探測子站與子站之間的連線(藍色線)為邊界，將線內(nèi)區(qū)域定義為探測網(wǎng)內(nèi)區(qū)域，線外區(qū)域即監(jiān)測邊緣區(qū)域定義為探測網(wǎng)外區(qū)域.在處理雷聲信號時，將極短時間內(nèi)連續(xù)的多點聲源近似處理為單個定位計算所需要的雷聲點聲源信號[13].在接收雷電放電時，聲波信號雷聲站中的聲波麥克風(M0～M4)陣列布置方式為十字陣列排布，如圖3所示.

圖2 TDOA探測子站和雷聲探測站布局Fig.2 Site location of TDOA lightning location sensors and acoustic sensors

圖3 十字麥克風陣列示意圖Fig.3 A diagram of crossed microphone array

當雷電發(fā)生在目標監(jiān)測區(qū)域內(nèi)時，4個TDOA探測子站會采集到一套相應(yīng)的雷電電磁信息，依據(jù)測得的雷電電磁時間信息，可列出由3個互相獨立方程組成的方程組：

(9)

式中：(x3,y3)和(x4,y4)分別為T3和T4的位置坐標.

利用式(9)計算得到雷電初定位結(jié)果(單一點)，記為L0=(x0,y0).根據(jù)初定位坐標判斷雷電發(fā)生位置所在的子區(qū)域.同時，利用子區(qū)域最小距離判定并選取合適的麥克風陣列，進行麥克風陣列所探雷聲信息的融合和求解.根據(jù)測得的雷聲聲波信息、麥克風陣列相對位置和聲波傳播特性，得到

(10)

式中：(xi,yi)和(xj,yj)分別為麥克風Mi和Mj的坐標；Δtij為雷聲到達時差；v為聲波傳播速度.

由此可知，對于包含多個麥克風的測量陣列系統(tǒng)，如由5個麥克風組成的測量陣列系統(tǒng)可以得到10個約束方程，即使在少數(shù)麥克風傳感器工作失效的特殊情況下，其仍能夠獲得足夠的雷聲信息.納入雷聲聲波信息后，方程組中的方程數(shù)目大于未知數(shù)數(shù)目，可以采用L-M算法進行迭代計算，求出雷電位置的最優(yōu)解.

基于L-M算法和聲波信息融合的TDOA雷電定位算法的整體流程如圖4所示.利用TDOA探測子站測得的電磁信號及其時間信息，判定雷電初定位結(jié)果.根據(jù)L0所屬子區(qū)域k′的判斷結(jié)果，調(diào)取距離子區(qū)域k′最近的雷聲站n′測得的聲波信息，并將其納入非線性方程組，利用L-M迭代算法計算更精確的雷電位置.

圖4 基于L-M算法和聲波信息融合的TDOA雷電定位算法流程Fig.4 Flow chart of TDOA lightning location approach considering acoustic information and L-M algorithm

2.3 算法驗證

在MATLAB中完成所提方法的算法編寫，并利用蒙特卡羅法[24]對雷電定位法進行精度評估.以邊長為40 km的正方形雷電監(jiān)測目標區(qū)域為例，據(jù)圖2中x軸和y軸建立坐標系4個時差法探測子站 T1～T4探頭位置的坐標分別為(15, 15)(-15, 15)(-15,-15)和(15,-15)；4個雷聲站N1～N4的中心麥克風位置的坐標分別為(8, 8)(-8, 8)(-8,-8)和(8,-8)；圖3中M1～M4與M0的距離均為20 m.

依據(jù)蒙特卡羅法，以2 km為空間步長將 40 km×40 km的雷電監(jiān)測區(qū)域進行網(wǎng)格劃分，將區(qū)域內(nèi)的441(21×21)個格點分別作為雷電發(fā)生位置進行雷電定位精度評估，具體步驟如下：

步驟1計算在某格點處產(chǎn)生的雷電電磁波和聲波信號到達每個探頭的準確時間.

步驟3分別利用TDOA定位的傳統(tǒng)法和所提改進法進行雷電位置計算，得到相應(yīng)的雷電定位結(jié)果.

步驟4評估雷電定位結(jié)果的誤差，即雷電發(fā)生的實際位置與定位計算結(jié)果之間的絕對差值.

步驟5重復步驟2至步驟4共100次，得到以該格點作為雷電發(fā)生位置的100次定位結(jié)果和定位誤差，并將100次定位誤差的均值作為該格點的最終定位誤差值.

圖5 雷電定位誤差分布Fig.5 Error distribution of lightning location

表1 不同子區(qū)域的雷電定位誤差均值

TDOA雷電定位傳統(tǒng)法和改進法在不同子區(qū)域內(nèi)的誤差均值分布如圖6所示.結(jié)合圖5可知，改進法利用雷電監(jiān)測區(qū)域測得的電磁時差信息開展雷電活動的初步定位，根據(jù)初步定位結(jié)果動態(tài)選取距離最近的雷聲站信息，同時融合電磁和聲波信息進行雷電活動精準定位，以改善整個監(jiān)測區(qū)域的雷電定位效果.傳統(tǒng)法在探測網(wǎng)內(nèi)的定位誤差為120～200 m，在探測網(wǎng)外的定位誤差接近300 m，在整個目標監(jiān)測區(qū)域內(nèi)的定位誤差均值為203.2 m，均方差為82.8；改進改法在探測網(wǎng)內(nèi)的定位誤差為90～110 m，在整個目標監(jiān)測區(qū)域內(nèi)的定位誤差均值為108.4 m，定位精度提高了46.6%，均方差為27.3.此外，利用傳統(tǒng)法定位時，位于監(jiān)測區(qū)域邊緣的雷電定位誤差明顯增大；而利用改進法進行定位時，位于監(jiān)測區(qū)域邊緣的雷電定位誤差得到改善，雷電定位精度比傳統(tǒng)法提高了51.2%.

圖6 不同子區(qū)域的雷電定位誤差均值分布Fig.6 Average lightning location error distribution in different sub-regions

3 結(jié)語

本文提出一種基于L-M算法和聲波信息融合的TDOA雷電定位改進方法，根據(jù)TDOA探測子站和雷聲探測站的位置信息，將目標監(jiān)測區(qū)域劃分為16個子區(qū)域.利用電磁時差信息完成雷電的初步定位計算，并根據(jù)初定位結(jié)果動態(tài)納入雷聲站探測得到的聲波信息.同時，利用L-M算法對融合的聲波信息進行計算，求解雷電位置最優(yōu)值.采用蒙特卡羅法分別對TDOA雷電定位傳統(tǒng)法和所提改進法進行精度評估，發(fā)現(xiàn)傳統(tǒng)TDOA雷電定位系統(tǒng)的定位誤差在整個監(jiān)測區(qū)域分布不均，在目標區(qū)域邊緣的雷電定位精度下降明顯，誤差均值為203.2 m，所算結(jié)果均方差為82.8；所提方法經(jīng)過聲波信息融合和L-M算法改進后，能夠有效降低區(qū)域誤差不均勻的程度并提高雷電監(jiān)測區(qū)域邊緣的定位精度，誤差均值為108.4 m，目標監(jiān)測區(qū)域整體定位精度提高了46.6%，所算結(jié)果均方差為27.3.所提方法能夠?qū)σ延欣纂姸ㄎ幌到y(tǒng)的改進和聲波信息的融合利用提供參考.