基于正交传播算子的闪电宽带甚高频辐射源定位方法研究*

李书磊 邱实 石立华 李云 段艳涛

(陆军工程大学,电磁环境效应与光电工程国家级重点实验室,南京 210007)

1 引 言

闪电甚高频(VHF)辐射源定位技术可以实现闪电放电通道时空演变过程的高分辨率成像分析,对闪电放电机理的认识及其电磁环境效应防护具有重要的意义[1-4].

自Shao等[5,6]提出宽带干涉仪后,基于短基线(基线长度通常几米到几十米)宽带VHF定位的各类系统在闪电物理研究中发挥了重要作用.Ushio等[7]采用分段触发的方式对一次空中触发闪电的下行负先导产生的VHF辐射进行观测和特征分析,但没有探测到上行正先导产生的辐射.董万胜等[8]利用三阵元正交基线宽带干涉仪阵列基于干涉法对人工触发闪电进行了观测研究,首次探测到人工触发闪电上行正先导的辐射.针对阵列天线间的信号时延估计,邱实等[9]将时域互相关技术应用于宽带VHF辐射源定位.Sun等[10]提出基于小波变换的广义互相关时延估计,利用基小波在不同尺度下的频谱对互功率谱密度进行加权,降低了噪声的影响.Stock等[11]设计了连续采集宽带干涉仪,并利用广义互相关算法提高了弱辐射源的探测能力.

闪电起始的双向先导理论认为强电场环境下的闪电通道将分别以正、负先导的方式同时双向发展[12,13].然而,正击穿的辐射功率比负击穿至少低一个量级,传统VHF辐射源定位方法往往导致负击穿掩盖掉同时发展的正击穿.正先导的发展过程观测和特征分析是澄清闪电发展过程的一个亟待解决的问题[14].此外,反冲先导是否同样存在未被VHF探测到的正先导部分仍尚待观测[15].这对定位算法的弱辐射源定位和多源定位能力提出了较高要求.针对弱辐射源定位,通过多阵元对辐射信号进行聚焦,具有较高的空间分辨能力、较高的信号增益及较强的抗干扰能力.Wang等[16,17]基于七阵元“L”型天线阵,首次采用时间反转技术对人工触发闪电放电通道进行定位,并提出了频域时间反转(FDTR)成像方法,证明了该方法在弱辐射源探测中的有效性.该方法实质上是谱估计理论在空域的扩展形式,但阵列的角度分辨力受到空域“傅里叶限”的限制,在一定程度上降低了定位精度,导致位于一个波束宽度内的多个空间目标不可分辨.另外,Wang等[17]对该方法的双源定位能力进行了实验验证,但并未在实际闪电定位中应用.Marcos等[18,19]提出传播算子算法进行波达方向(direct of arrival,DOA)估计,通过线性运算构造传播算子,利用子空间的正交性获得空间谱峰,提高了算法的空间分辨力,同时避免矩阵特征分解所带来的计算量.该方法因其计算复杂度低、易于工程实现等优点广泛应用于雷达目标定位、移动通信领域.

针对时间反转波束宽度较宽、角度分辨力低的缺点,本文将正交传播算子(OPM)方法应用于闪电宽带VHF辐射源定位,针对宽带VHF信号,采用非相干子空间处理方法(ISM)将带宽内的有效频点进行平均,以减小噪声干扰.并将定位结果与FDTR算法进行了对比,表明该算法定位弱辐射源的有效性和同窗双源事件定位的优势.

2 宽带OPM闪电辐射源定位方法的原理

正交传播算子方法是一种基于空间谱估计理论的经典DOA估计方法.利用宽带VHF阵列接收数据线性构造特征子空间,突破了瑞利限的限制,可实现高分辨率的DOA估计[20,21].该方法最初针对窄带DOA估计提出,而闪电辐射的VHF信号是一种典型的宽带信号,频率的不同导致了信号子空间的差异,因此需要通过ISM方法将有效频点加权,然后进行闪电辐射源DOA估计.

假设K个闪电辐射信号s1(t),s2(t),…,sk(t)(k=1,2,…,K)分别来自〈θ1,φ1〉,〈θ2,φ2〉,…,〈θk,φk〉方向,其中θ和φ分别表示辐射源的方位角和仰角,阵元噪声假设为不相关的高斯噪声.则对于M(要求M＞K)元VHF阵列,输出矢量可表示为

其中,fj表示对VHF信号进行J点离散傅里叶变换后的第j个频率点(j=0,1,…,J—1).X(fj)=[X1(fj),X2(fj),···,XM(fj)]T表示频率点fj处的阵列接收数据,A(fj,〈θ,φ〉)为VHF天线阵列流型,其形式为

假设闪电辐射源信号相互独立,则天线阵列流型矩阵A是列满秩的,则矩阵A有K行线性无关,其他M—K行可由这K行线性表示.因此可将天线阵列流形矩阵分解为

由于矩阵A1为非奇异矩阵,因此矩阵A2可以用矩阵A1的线性表达,即

其中,矩阵P∈CK×(M-K)称为传播算子.令则有

其中,IM-K和0分别为M-K阶单位矩阵和(M-K)×K阶零矩阵.上式表明矩阵Q的列向量和VHF天线阵列方向矢量正交,则矩阵Q的列向量所张成的子空间即为噪声子空间.将矩阵Q正交化即为正交传播算子Q0:

定义窄带OPM定位算法的空间谱如下:

最后,针对闪电宽带VHF信号,根据ISM算法的原理,将频带内有效频点进行平均,从而得到宽带VHF信号的OPM算法空间谱:

在实际应用中,由于VHF天线阵列的流型矩阵是未知的,因此一般从阵列快拍数据协方差矩阵来估计Q0,基本步骤为:

1)利用VHF天线阵列接收数据计算协方差矩阵R,表达式为

2)对矩阵R进行分块,

3)由于VHF天线接收数据中噪声的存在,传播算子的估计由最小化问题得到此问题的解为并进一步得到

4)根据(8)式构造宽带VHF空间谱,搜索空间谱获得闪电辐射源DOA估计结果; 实际上,通过上述步骤构造的“倒谱”是利用子空间正交性使得原来在辐射源位置零值的谱变成趋于无限大,而非辐射源点处的值为相对有限值,这就将辐射源在区域中与非辐射源位置上的谱幅度区别开来; 而实际由于噪声、阵列互耦、辐射源相干等因素,使得空间谱中出现非源点的伪峰; 为了实现辐射源同窗双源定位,需要首先判断是否存在双源; 实测数据分析表明,搜索次峰值与主峰值的比大于0.5时,可以认为空间中同时存在两个辐射源,并分别进行定位;

5)质量控制方法,首先,采用密度聚类分析[22]的方法将离散的伪点去除.设置一组邻域参数(簇的最小点个数和临近半径,此处分别设置为20和0.03)判断辐射源点是否紧密相连,从而去除离散伪点; 然后,使用能量比(ER)作为辅助判据,进一步滤除聚成一簇的伪点,其表达式如(9)式所示,即信号方向的能量与所有可能方向的总能量的比值,

3 数值仿真分析

为了分析算法性能,通过数值仿真研究了不同信噪比时算法的定位误差、以及空间谱半峰值宽度、角度分辨力等表征双源定位性能的参数,并与FDTR算法进行了对比.

参考实际宽带VHF信号波形特征,仿真时采用正弦高斯调制信号作为辐射源信号,其数学表达式如下:

式中,A0表示信号幅度,实际处理时进行归一化处理;f0表示信号的中心频率,此处设置为100 MHz;τ1和τ2用于调节信号的带宽,参考VHF天线实际采集信号带宽,此处均设置为20 ns.

天线阵列的设置型式,参考实验室野外观测试验的设置方式,采用7天线“L”型平面阵列,阵元间距设置为9 m,详细情况参见4.1节.仿真时辐射源的方位角和仰角事先指定,然后根据辐射源与各天线之间的空间坐标关系获得各天线接收信号的延迟,最后利用定位算法计算辐射源的位置.

3.1 单一辐射源的定位误差分析

设置单一辐射源S的位置为θ=120.25°,φ=30.50°.为了定量评估定位结果的准确性,定义定位误差为DOA估计结果与辐射源实际位置的球面距离,以角度表示.此处规定定位误差小于1°时,即视为算法正确识别目标.针对不同信噪比(SNR)分别进行10000次独立重复试验,平均定位误差与信噪比的关系如图1所示.

图1 定位误差与信噪比的关系Fig.1.Relationship between locating error and SNR.

由图1可得,仅对单一辐射源而言,宽带OPM定位算法和宽带FDTR算法的定位性能基本一致,二者均能在信噪比大于—12 dB时对辐射源进行有效定位,且其定位误差均在0.3°以下,在一定程度上反映了算法的鲁棒性.当信噪比小于—12 dB时,两种算法仅能以一定的概率正确识别目标,识别概率如表1所列.随信噪比下降,识别概率相应降低,同一信噪比下两种算法的识别概率相近.

表1 低信噪比条件下两种算法的识别概率的比较Table 1. Comparison of recognition probabilities of two algorithms under low SNR.

3.2 双辐射源的分辨能力分析

闪电发展过程中多分支通道同时发展以及闪电起始等双向先导事件,增加了双源同窗事件发生的概率.

为了定量说明宽带OPM辐射源定位方法对于双辐射源的定位性能,分别定义空间谱半峰值宽度(HPW)和角度分辨力.其中,空间谱HPW为空间谱峰值半功率点的宽度,主要反映空间谱的尖锐程度,如图2(a)所示.图2(b)和图2(c)分别为方位角和仰角所对应的HPW与SNR的关系.

由图2可得,在信噪比大于—12 dB时,宽带OPM定位方法的半峰值宽度比FDTR方法小,表明其谱峰更加尖锐,旁瓣相对平缓,能量较为集中.并且随信噪比的增大,其HPW呈减小趋势,而FDTR方法的HPW几乎不随SNR而变化.当信噪比小于—12 dB时,两种方法的HPW都随着信噪比的减小而增大,表明算法的分辨力下降,此时两种方法的HPW几乎相同.

角度分辨力则是指定位算法分辨两个相近辐射源的能力.在进行空间谱搜索时,若相近的两个辐射源可分辨,那么其对应的主峰与次峰要大于最小可分辨变角度.设置两个等幅不相干辐射源,SNR设为0dB.分别利用两种算法对不同入射方向的辐射源进行DOA估计,其最小可分辨角度如表2所列.其中,最小可分辨变角度通过球面距离表示.

图2 方位角和仰角的HPW与SNR的关系(a)HPW示意;(b)方位角;(c)仰角Fig.2.Relationship between HPW and SNR of azimuth and elevation:(a)Sketch map of HPW;(b)azimuth;(c)elevation.

表2 不同入射角度时两种方法的最小可分辨角度Table 2. Minimum distinguishable angles of two methods at different incident angles.

由表2可得,宽带OPM方法的最小可分辨角度小于宽带FDTR方法,表明基于子空间正交构造空间谱的定位方法具有更好的角度分辨能力,从而具有更好的双源定位能力和定位精度,表明OPM方法在双源定位方面的优势.

4 人工触发闪电试验观测结果分析

4.1 观测系统

自 2016年 起,SLOT(Jiangsu Lightning Observation Team)团队持续在江苏苏北地区开展人工触发闪电试验,该地区属江苏省内雷电多发区,且视野开阔,具备闪电的长期观测条件.试验场包括两个闪电观测站点和试验站点,其中观测站点A处布置了一套多天线辐射源连续观测系统(multiple-antennas radiation continuous obser vation system,MARCOS),采用 7个天线呈“L”型平面阵列,阵元间距设置为9 m,阵列设置示意图如图3所示.同时利用8通道LeCroy高速数字存储示波器采集阵元的VHF信号,采样率为500 MHz,存储深度为250 Mpts/ch,可连续记录500 ms闪电数据.该系统通过增加阵元数量、增大动态接收范围、采用连续记录模式,提高了对弱辐射源的探测能力.此外,该站点还布置有同步记录电磁场的快、慢天线、磁场天线等,以及进行光学观测的高速摄像等.

图3 MARCOS系统天线布置示意图Fig.3.Schematic diagram of MARCOS system antennas layout.

4.2 观测结果分析

2017年8月1日23时左右,雷暴云自东南向西北运动,于23时05分53秒采用经典触发方式成功触发闪电(记为Trig230553).触发时地面大气电场为—13.5 kV/m(大气电学符号规定),据此判断此次触发为负极性触发.触发闪电的快电场波形记录如图4所示,0时刻代表记录信号的触发时刻.此次触发闪电共有5次直窜先导-回击过程,整个放电过程持续时间约600 ms.

图4 Trig230553的快电场波形Fig.4.Fast electric field waveform of Trig230553.

利用宽带OPM闪电VHF辐射源定位算法,对此次触发闪电的放电通道进行了定位,其中滑动窗的窗长设置为512个采样点(1.024 µs),重叠率为75%.整个记录时段的闪电成像结果如图5(a)所示,并将其与同站设置的高速摄像拍摄的结果(图5(b))进行了对比.两者在闪电通道形状上具有较高的一致性,但是高速摄像受限于视野范围,没有获得该次闪电其他的发展过程,而基于VHF系统的成像可以对高速摄像视野以外的通道有较好的描述.此次触发闪电的上行正先导起始于A点,且放电通道分支较少,当上行正先导发展到仰角为45°左右时出现明显分支通道,之后仅在正先导发展后期可以观察到两条细小分支通道.

在上行正先导发展末期,在正先导主通道和分支通道末端分别引发了数次反冲先导过程,该阶段持续时间大约为100 ms.大约210 ms时,引发了三次沿正先导主通道反向传输的直窜先导,其持续时间分别为3,1,0.7 ms,到达地面后均引发回击,其中第二次直窜先导过程如图6所示.这三次直窜先导过程均起始于相同的放电区域B处,沿之前上行正先导开辟的通道向地面方向发展,定位通道连续、清晰,且没有其他分支通道.

图5 Trig230553放电通道的定位结果(a)VHF阵列定位结果;(b)高速摄像结果Fig.5.Locating results of discharge channels of Trig230553:(a)Mapping results of VHF arrays;(b)HSV results.

图6 第二次直窜先导过程Fig.6.The second dart-leader process.

4.3 宽带OPM闪电定位方法和FDTR方法的对比分析

为了比较宽带OPM闪电定位方法和FDTR方法对弱辐射源的定位性能,选取该触发闪电上行正先导阶段的数据,分别利用两种方法进行通道成像,结果如图7所示.

由图7可得,两种方法均可以清晰地分辨出闪电发展通道的结构,且具有较高的一致性,对通道发展过程的描述基本类似.需要注意的是,由于两种方法空间谱能量集中的程度不同,为了便于比较,所以在滤除伪点时只采用CR阈值法,即两种方法设定相同阈值(此处设置为0.3).针对上行正先导阶段,对定位点数进行统计,OPM方法和FDTR方法的定位辐射源数分别为27015和24446.宽带OPM算法成功解算的定位点数量更多,这主要是由于宽带OPM方法的角度分辨力高于FDTR方法,对于多通道同时发展的情形,其成像结果的细节表现更优.

图7 Trig230553上行正先导成像结果比较(a)宽带OPM方法定位结果;(b)宽带FDTR方法定位结果Fig.7.Comparison of mapping results for upward positive leader of Trig230553:(a)Mapping results of OPM method;(b)mapping results of FDTR method.

4.4 双源定位性能分析

闪电发展过程中,多分支通道同时发展时使得VHF天线接收信号发生交错,因此滑动窗口时可能存在双源同窗的情况.传统干涉法或互相关方法通过设置高重叠率滑动窗可在一定程度上提高时间分辨率,但对于时间上相互交错的多辐射源信号无法处理.宽带OPM定位方法在空域对信号进行处理,其空间谱表示辐射源在空间各个方向上的能量分布.在搜索空间谱最优解时,按照次峰和主峰能量比大于0.5时判定双源同时存在,否则为单一辐射源.Trig230553发展到—5 ms时,上行正先导出现明显分支(图7中以P1和P2标识),采用OPM方法对5 ms时间内的两分支进行了双源自动识别定位,结果如图8所示.

图8 宽带OPM方法双源定位结果Fig.8.Mapping results of broad band OPM method for two sources of Trig230553.

图8反映了该时间内仰角随时间的变化,P1和P2分别对应于图7中的上行正先导的分支.红色菱形表示每个窗口空间谱主峰所对应的定位结果,黄色正方形表示次峰所对应的结果,两者为同一时间窗对应的辐射源,如图8中局部放大图所示.在5 ms时间段内,宽带OPM方法可定位到有效双源的窗口数为37.由于估计偏保守,实际同时存在双源的窗口个数可能更多.自VHF定位技术应用闪电探测以来,均只能实现同一时间窗口单个辐射源的定位,图8给出的定位结果首次实现了对同一时间窗的双源同时定位.双辐射源的定位结果证实了闪电各个分支上的一些击穿事件存在同时发展的现象,这对于揭示闪电放电发生发展的物理本质具有重要的参考价值和意义,尤其对于闪电起始等双向先导发展模式的研究以及模型的建立提供了新的工具和思路.

利用双源定位算法,对此次触发闪电的上行正先导阶段进行了双源定位,定位到有效点数为32943.但通道形状基本与单源定位结果基本类似,这主要是由于此次触发闪电双源定位的通道单源基本也可以交替发展的方式进行定位.另外,此处双源定位只能给出幅度比小于2的两个辐射源的位置,需要增加天线阵列数目等方法进一步提高物理孔径,增强弱辐射源的定位能力.

由于OPM方法的角度分辨力优于FDTR方法,具有较小的最小分辨角度,旁瓣抑制能力更强,信源方向的空间谱更为突出,从而更适于空间相距较近的分支通道同时发展时的定位.

5 结 论

本文将基于空间谱估计理论的OPM方法应用于闪电宽带VHF辐射源定位,实现了闪电放电通道时空演变过程的成像.该方法首先通过VHF阵列接收数据协方差矩阵的线性分解形成正交化传播算子,并利用子空间正交构造空间谱,然后进行空间谱搜索确定辐射源的位置.针对宽带VHF信号,采用ISM方法将带宽内的有效频点进行平均,减小噪声干扰.

为了验证算法的有效性,对OPM定位方法和FDTR方法的定位误差、空间谱半峰值宽度、角度分辨力进行了数值仿真.结果表明,针对单一辐射源,宽带OPM定位算法的定位误差与FDTR算法相近,但OPM算法的角度分辨力优于FDTR算法,对于分辨双辐射源性能更优.利用宽带OPM定位算法分析了一次人工触发闪电放电通道的时空发展过程,结果表明该算法可以较高的时空分辨率清晰地描绘出闪电通道的基本结构及放电通道的发展过程,并且OPM算法对双辐射源的定位能力有较大优势.下一步将利用该算法双源定位的优势,进一步开展对闪电起始先导发展过程的定位研究,揭示正、负击穿同时存在时的发展规律.