闪电是通常发生于雷暴过程中的一种瞬态大电流放电现象，研究闪电通道发展结构对于认识闪电的物理过程和雷暴起电、放电特征有重要意义^[1-4]。研究表明，闪电过程会产生两种类型的声音辐射:一是由闪电放电热通道的急速膨胀引起的，属于可闻雷声；二是闪电使云中电场快速减小时由存储在雷暴云静电场中的能量转换产生的声音，属于次声^{[1, 5]}。本文讨论的是前者，可闻雷声的探测可通过麦克风阵列实现，并通过采集设备进行同步采集，再利用合理的算法进行定位，得到的闪电过程中不同位置的声源即可重建闪电通道的三维空间分布。由于雷声的传播距离一般不会超过30 km，所以雷声定位可用于对近距离的闪电活动进行探测，且声音传播速度较慢，所以利用雷声对雷电距离进行估计误差不会太大，因此雷声定位技术在局地区域的雷电探测方面可作为其他雷电定位技术^[6-7](如基于电磁波到达时间差的技术) 的补充手段。

20世纪70年代，Few^[8]首次提出了基于雷声观测重建闪电通道的方法，实验中架设了4个麦克风组成的Y型结构阵列，数据采样率为1 KS/s，基线长度为30 m。Few利用雷声信号到达阵列的时间差和雷声信号与电磁信号到达的时间差，大体刻画出了闪电通道结构。MacGorman^[9]在架设麦克风阵列时，仪器的数据采样率为500 S/s，4个麦克风组成的边长为50 m的正方形阵列，并考虑了雷声在大气中传播时受到温度、风速和风切变等因素的影响，进一步改善了定位结果。Akiyama等^[10]考虑了地面噪声对信号的影响，将其放置在距离地面60 cm高度。采用由4个麦克风组成的阵列，分别被架设在边长为3.5 m某一顶角为60°的菱形的4个顶点上。由于采用的采集设备的采样率达到了100 KS/s，所以Akiyama架设的麦克风阵列的基线长度有了明显的缩短。Arechiga等^[11]在触发闪电附近，放置频率范围为3.3~500 Hz的麦克风阵列，采样率为1 KS/s，并与能够准确定位闪电通道的LMA (Lightning Mapping Array, 闪电定位阵列) 系统比较，确定了该方法的精度：在近距离 ( < 5 km) 处误差为100 m，在远距离 (>6 km) 处误差为650 m。Johnson等^[12]使用4个宽频麦克风组成的阵列，结合缓慢的搜索和距离测距用来区别三维雷声区域 (大约12 km)，采样率为1 KS/s。章涵等^[13]采用4个麦克风组成阵列，3个在xy平面，1个在z轴方向，基线长为1 m，采样率为100 KS/s，记录时间为35 s，大大提高了阵列的便携性。Qiu等^[14]首次采用了二维VHF干涉仪和单站三维雷声定位相结合的系统对闪电通道的研究。采用不在同一直线上的3个麦克风，基线长12.5 m，采样率为100 KS/s，持续时间为30 s。VHF定位系统着重定位先导分支，而雷声定位系统则着重定位主通道，他结合两种方法的优缺点弥补了各自的不足，较完整地重构了闪电通道。

Few^[8]、MacGorman等^[9]、章涵等^[13]和Qiu等^[14]都采用互相关函数法对雷声信号进行处理，以得到雷声到达不同麦克风的时间差。Few等^[15]讨论了雷声重建方法的精度，提出了互相关函数法在雷声定位中的局限性：互相关函数法不能区分出同一时刻到达麦克风阵列的不同闪电分支。Akiyama等^[10]提出了互功率谱法能够区分不同通道产生的不同声波频率，弥补了互相关函数法的不足。本文在文献[10]研究的基础上考虑了大气介质对雷声传播的影响，提出了改进的互功率谱相位法，增强了雷声信号的抗噪能力，提高了定位的精度。

从2006年开始，中国气象科学研究院每年夏季都联合广东省气象部门开展“广东野外雷电综合观测试验”，包括人工触发闪电和自然闪电综合观测试验^[16-17]。本文使用的定位系统架设在广东省气象局楼顶，对击中高建筑物的自然闪电进行观测。

麦克风阵列的结构和性能指标在文献[13]中已详细描述，本文不再赘述。几个关键的指标如下：麦克风响应的频率为15~20000 Hz，基线长度为1 m，采样率为100 KS/s，记录长度为35 s，预触发时间为5 s。高速摄像系统架设在阵列附近，能够同时记录闪电通道发展过程，本文将用拍摄的闪电照片与麦克风定位的结果比较，对雷声定位结果进行检验。

本系统中，麦克风阵列的基线长度为1 m，对于几百米或更远的闪电，通道产生的声音信号传播到麦克风阵列时可以被近似看成是平面波。通过章涵等^[13]的分析，对于一个由多个麦克风组成的阵列，雷声信号到达每个麦克风的时间不同，利用到达时间差可以计算得到声源对于麦克风阵列的方位角和仰角。闪电产生的电磁辐射会使声波观测系统产生感应信号，利用声波和电磁波到达麦克风的时间差可以计算得到声源的距离，结合得到的方向信息就可得到声源的三维位置。

图 1是多分叉闪电通道导致的雷声重叠现象示意图，其中通道至麦克风阵列的3条虚线长度相同。由于任意时间点阵列的实测信号是所有同时传播到阵列的源信号的总和，而互相关函数法对信号振幅具有选择性，只能研究时间窗内幅值最大的声源点^[13]，因此用互相关函数法计算不同麦克风之间的时延很难区分出如图 1所示的闪电通道结构。

图 1

图 1. 多分叉闪电通道导致的雷声重叠现象示意图 Fig 1. Schematic of the thunder overlapping caused by multi-branch lightning channel

考虑到重建如图 1所示的闪电通道结构，Akiyama等^[10]提出了用频率区分声源点的方法。假设不同通道产生的声波频率不同，每个频率点的时延τ_ij即可求得。时延τ_ij可由互功率谱S_ij(ω) 的虚部和实部求得，表达式如下：

(1)

式 (1) 中，ω是雷声的角频率。因此雷声信号在时间窗内可获得不同位置的声源点。

时延τ_ij可由式 (1) 计算，每个频率点对应一个值，因此可分辨出不同的闪电通道。但互功率谱相位法也有一定局限性，两个麦克风之间的相位差有2nπ的不确定性。为了避免这个不确定性，麦克风之间的距离必须满足L < 2π|v_p|/ωcoaα，其中n和v_p分别是整数和相速度。

基于时延估计的声源定位方法是近几年发展起来的一种传感器阵列声源定位方法，在运算量上远远小于可控波束法和谱估计法，且适当改进后，即使在一定的噪声和混响环境下也有较好的定位精度，并适合于实时定位的应用^[18-20]。基于时延估计方法的声源定位流程分为时延估计和声源定位两部分^[21]，其中时延估计方法的精度和鲁棒性是关系到声源定位精确与否的关键因素。本文先用时延方法计算出各麦克风之间的到达时间时延，再通过时延计算出声源的方向信息。

时延估计方法的种类很多，常用的有互相关函数法、最小均方自适应滤波、互功率谱相位法等^[22-23]。本文首先介绍互相关函数法和互功率谱相位法，考虑了在闪电发生过程中会产生高频噪声、雷声到达地面反射时产生的混响以及雷声信号上叠加的各种噪声 (如采集电路本身的噪声) 等，提出了改进的互功率谱相位时延估计方法，该方法对噪声和混响均具有很强的鲁棒性。

假设两路接收信号的信号模型^[24]：

(2)

(3)

其中，S(t) 为雷声信号，n₁(t)，n₂(t) 为互不相关的高斯白噪声，S(t) 和n₁(t)，n₂(t) 也互不相关。τ₁，τ₂分别是雷声源点传播到麦克风1和麦克风2的时间，τ₁₂=τ₁-τ₂即为麦克风1、麦克风2间的时延。

x₁(t) 和x₂(t) 互相关函数可表示为

(4)

由于S(t), n₁(t) 和n₂(t) 三者互不相关，上式可简化为

(5)

由自相关函数的性质可知，当τ-(τ₁-τ₂)=0时R₁₂最大，此时R₁₂(τ) 对应的最大值即为麦克风1、麦克风2的时延。

利用互相关函数法进行时延估计的主要特点是方法简单。但该方法有两个明显的不足之处，一是其假定了信号和噪声及噪声之间均不相关，这在有些情况下不一定能得到满足；二是相关函数是在一种严格数学意义上的统计平均或在平稳遍历条件下替代统计平均的无穷时间平均。而在实际应用中，严格数学意义上的统计平均或无穷平均是不可能做到的，而只能用有限的时间平均替代无穷平均或统计平均，这种短时处理导致噪声对相关函数的影响不容忽略，从而使R₁₂(τ) 的峰值不明显，降低了时延估计的精度。

由于互相关函数与互功率谱满足傅里叶关系^[25]，

(6)

其中，G₁₂(ω) 是两路信号的互功率谱，Ψ₁₂(ω) 是加权函数。图 2是互功率谱相位法时延估计原理框图，针对不同环境，可选择不同加权函数，使R₁₂(τ) 具有比较尖锐的峰值。R₁₂(τ) 峰值处为两个麦克风间的时延，但低信噪比环境使这种分析很不稳定，因此选择合适的加权函数是本文的难点。

图 2

图 2. 互功率谱相位法时延估计原理框图 Fig 2. Block diagram of the time delay estimation using cross-power spectrum phase method

雷声在大气中传播时，由于气溶胶及大气层之间的折射和反射，雨和水汽对雷声的反射和吸收等^[26-27]，会产生类似声波在水中传播的多途效应，严重影响了雷声信号处理的性能。赵真等^[28]研究表明：在多途条件下，幅度平方加权能够削弱次峰的干扰，使R₁₂(τ) 的峰值更加尖锐，进而能大大减少由于信噪比降低及信号起伏引起的估计误差，提高定位的精度。

为了检验这两种方法，在仿真实验中计算机模拟两段多频率声源信号S₁(t)，S₂(t)，采样频率为2 kHz，采样长度为1024，并加上加性高斯白噪声。如图 3所示，S₂(t) 是将S₁(t) 时延时间τ后的信号。

图 3

图 3. 模拟的两段多频率声音信号 (a) 及其中截取的一段信号 (b) Fig 3. Two simulated sound signals with multiple frequency (a) and a section of the signal (b)

仿真方案：分别取3种常用的加权函数, 具体见表 1。其中，G₁₂(ω) 是接收信号S₁(t) 和S₂(t) 互功率谱。

表 1

表 1 常用的加权函数 Table 1 Common weighting functions 编号 名称 加权函数 1 互相关函数法 (CCF) 1 2 互功率谱相位法 (CSP) Ψ12(ω)=1 3 幅度平方加权互功率谱相位法 (ASW-CSP) Ψ12(ω)=|G12(ω)|2

表 1 常用的加权函数 Table 1 Common weighting functions

仿真结果：信噪比为-10 dB (噪声信号很大) 时3种方法的比较。

由图 4看出，在信噪比很低 (-10 dB) 时，互相关函数法的互相关函数的峰值很不明显，几乎淹没在干扰峰之中，互功率谱相位法的互功率谱函数虽然峰值明显，但也出现了一些较强的干扰，而幅度平方加权互功率谱相位法的互功率谱函数峰值明显，对干扰信号有较好的抑制能力。

图 4

图 4. 信噪比为-10 dB时3种方法对干扰的抑制能力 Fig 4. Three methods for interference suppression ability at signal to noise ratio of-10 dB

选取2012年5月5日08:01:00(北京时，下同) 广州发生的一次自然闪电 (编号F1203)。图 5a是麦克风阵列的实测信号，其中包括电磁信号、雷声信号和噪声干扰。图 5b是实测信号经过低通滤波器，滤掉叠加在实测信号中的噪声信号，以减少对雷声定位的干扰。

图 5

图 5. 麦克风探测到的F1203雷声信号 (a) 和经过滤波之后的信号 (b) Fig 5. The original thunder signal of lightning flash F1203 obtained by a microphone (a) and the filtered signal (b)

用互相关函数法 (CCF) 和改进的互功率谱相位法 (ASW-CSP) 对雷声声源定位，为了检验雷声定位方法的结果，与通道的照片信息进行了比较。

为便于对比，图 6c的视野范围与图 6a相同，将高速摄像照片和麦克风阵列定位点的方位角统一为以观测点正北方为0°，顺时针旋转0°~360°。在图 6c中，定位点与照片描绘的通道有大约1°~2°的偏差，可能是阵列架设时造成的。由高速摄像资料可知，此次闪电有两个接地点，左边接地点的方位角约为112°，右边接地点的方位角约为132°。由图 6c可知，利用ASW-CSP在这两个方位角附近均定位到了一些声源点，说明基于ASW-CSP的声源定位方法是可靠的。CCF虽然能够定位出右边闪电通道的大致形状，但定位点较离散，ASW-CSP的定位点集中，能够更好地描绘出右边闪电通道的细节 (如拐角处)，可能是因为CCF对噪声的抑制能力较弱；CCF在左边接地点附近和左边闪电通道均没有定位点，ASW-CSP定位点与左边闪电通道形状大致吻合，可以看到，因为CCF对信号振幅的选择性而不能分辨出此处的闪电通道分支，而ASW-CSP弥补了这种不足^{[8, 13]}。相比之下，ASW-CSP对闪电通道的定位效果明显优于CCF。

图 6

图 6. 闪电F1203的定位结果 (a) 高速摄像拍摄到的图像 (镜头焦距14 mm)，(b) 高速摄像拍摄到的图像 (镜头焦距24 mm)，(c) 两种方法定位的声源点与二维闪电通道的对比 Fig 6. The location results of lightning flash F1203(a) the high speed images of lightning F1203 (lens focal length is 14 mm), (b) the high speed images of lightning F1203(lens focal length is 24 mm), (c) the comparison between thunder source locations using two methods and two-dimensional channels

先通过麦克风阵列采集4组雷声数据，再通过改进的互功率谱相位法得到声源的方位信息，通过声电差估算出声源点到阵列的距离，根据距离和方位信息即可得到声源点的三维信息。

声音是以物质波动的形式来传播的，其大小与物体的密度和弹性有关，进而与温度和压强有关。在干燥空气中，声速的经验公式^[29]如下：

(7)

式 (7) 中，T为温度，单位：℃。潮湿空气的声速略有增加，但是增加幅度不到0.5%，大多数场合可以忽略不计。

雷声伴随闪电产生，一般认为可以听到的雷声被认为是加热的闪电通道迅速扩张而引起的。设雷声速度为v₁，闪电回击产生的电磁脉冲速度为v₂，计算声源点距离

(8)

由于v₂»v₁，可得到距离S=v₁Δt，其中Δt是电磁信号和雷声到达阵列的时间差。

图 7所示的闪电 (编号F0904) 发生于2009年8月24日19:08：04，由高速摄像照片确定的接地点为距离观测点约1420 m、高度为105 m的建筑物顶部。由图 8可以看出，该次自然闪电包含一次回击，发生在5 s (触发时刻) 处幅值很大的脉冲是由闪电回击产生的电磁信号，与伴随闪电产生的雷声信号到达麦克风的时间相差4.18 s，可估算出闪电过程距离探测点大约1410 m，与由高速摄像照片确定的位置的距离1420 m很接近，说明根据回击产生的电磁信号和雷声信号到达的时间差来估算闪电的距离具有较好的精度。

图 7

图 7. 闪电F0904的高速摄像照片 Fig 7. A high speed image of lightning flash F0904

图 8

图 8. 麦克风1实测得到的F0904雷声信号 Fig 8. The thunder record of Microphone 1 for F0904

本文用CCF和ASW-CSP对该过程也进行了通道三维重建，结果如图 9所示。由图 9a可以看出在闪电接地主通道处CCF和ASW-CSP均吻合较好，在分叉1和分叉2处ASW-CSP定位结果能够与高速摄像二维照片在闪电通道的分支部分有较好的对应，而CCF在其附近没有定位点，这是因为CCF算法基于振幅相关性^[10]不能区分同时到达阵列的不同声源点，且由于CCF抗噪能力较弱，定位点离散不能较好地描绘出通道形状。结果表明：在低信噪比环境和多分支闪电的情况下，ASW-CSP比CCF能够表现出更好的定位效果。

图 9

图 9. 闪电F0904的定位结果 (a) 两种方法定位的声源点与二维闪电通道的对比, (b) 雷声声源点的三维定位结果 Fig 9. The location results of lightning flash F0904 (a) the comparison between thunder source locations using two methods and two-dimensional lightning channels, (b) three-dimensional locations of entire thunder sources

本文针对原有雷声定位方法在多分支和抗噪能力方面的不足，提出了一种改进的互功率谱相位法，改善了定位效果，得到以下初步结果：

1) 互相关函数法和改进的互功率谱相位法分别得到雷声信号的到达时间差，并与高速摄像记录的二维照片对比，验证了两种方法的可靠性。

2) 比较了两种定位方法在闪电通道多分支的情况下，改进的互功率谱相位法能够区分出同时到达阵列的不同声源点。

3) 利用闪电回击产生的电磁信号与雷声信号到达阵列的时间差得出声源点的距离，从而反演出雷声声源点的三维信息，并描绘出闪电通道结构。

通过架设单站麦克风阵列对闪电通道进行三维定位，降低了对架设环境的要求和布站的成本。虽然定位结果精度不是很高，但对于近距离、多分叉、多接地点的闪电观测来说，单站麦克风阵列提供了一种简单可行的三维观测方案，对于小范围的闪电监测、预警以及雷电物理研究来讲有一定应用前景。目前本文采用的定位方法在实时性和定位能力上都有较好的效果，在下一步工作中，将与硬件结合研发出一套能够实时处理的、基于单站麦克风阵列的雷声声源三维定位系统，使其在闪电三维结构的相关研究中发挥其特有的作用。