三元次声阵定位火箭发射的方法

殷 昊，庞新良，张震川

（国民核生化灾害防护国家重点实验室，北京 102205）

0 引 言

次声源定位是一种典型的被动式探测定位技术，根据利用的信号信息和参数的不同，可大致将定位技术分为四类：利用信号的到达时间（Time of Arrival,TOA）定位、利用信号的到达时间差（Time Difference of Arrival,TDOA）定位、利用信号的到达方位角（Azimuth of Arrival,AOA）定位、利用接收信号的强度（Received Signal Strength,RSS）定位。其中，基于到达时间差的定位方法定位精度较高，实用性较强，是目前应用与研究最广泛的一种声源定位方法。

常用的是基于互相关的时延估计方法，其中最常用的方法主要是广义互相关法（Generalized Cross Correlation,GCC）。该方法被Knapp等[1]于1976年提出，其基本原理是在求取信号互相关函数之前对功率谱进行加权滤波，突出信号并抑制噪声干扰部分，从而突出相关函数在时延处的峰值，取其峰值进行时延估计。广义互相关方法采用了众多加权函数[2-6]，如由 Roth提出的Roth加权函数（The Roth Processor,ROTH）、平滑相关变换（The Smoothed Coherence Transform,SCOT）、相位变换（The Phase Transform,PHAT）、最大似然（Maximum Likelihood,ML）估计/HT（Hannan-Thomson）以及一些改善函数等，起到锐化峰值、抑制噪声的作用。近几年来，众多学者致力于 GCC时延估计方法的改进，小波变换[7]、时频分析、希尔伯特变换[8]、神经网络[9]、经验模态分解（Empirical Mode Decomposition,EMD）[10]等方法被引入到时延估计方法中，能有效地进行非平稳环境下的时延估计。为了更进一步精确时延量，杨亦春等[11]针对震前异常次声波一次GCC得到的时延量远远大于信号持续时间的情况，对信号采用二步法时延估计：第一步先根据能量分布粗略估计时延；第二步再用 GCC法计算精确时延，最终的时延量为两次时延之和。考虑到一次互相关法受噪声影响较大，文献[12]提出了二次相关法，可在更低的信噪比环境下较准确地估计时延。

因为火箭发射的时间与发射点的位置明确，对于检测次声时延估计、定位算法等具有重要意义，具有重要的研究价值[13]。火箭发射是一种运动次声源，其产生的次声波信号主要由火箭飞行时激发空气产生，因此近距离采集到的次声信号持续时间可长达 3～4 min，甚至更久。文献[14]考虑到次声监测网络孔径较大且信号相关性较小，获取的时延估计值准确度不高，因此采用最小方差法对阵列周边的区域进行扫描，进而得出声源位置的估计值，与实际位置相差数十千米。唐伟等[13]对朝鲜4.13“光明星3号”卫星发射次声信号进行分析，分别采用单台、两个次声台站两种方法对火箭爆炸区域进行估算。

由于次声波的远距离传播特性，大多数学者考虑的是次声源的远程定位。文献[15]中考虑在小于100 km的范围内，次声定位应用的缺乏，开发了一种 srcLoc方法，该方法是一种近场、严格基于TDOA的声学定位方法，其优势主要在于没有限制性的大气假设，试验结果表明该方法在96%的站点上比传统的定位方法的定位精度提高了2倍以上。而火箭次声信号的持续时间较长，且声源并非静止的，难以识别判定不同发射阶段所产生的次声信号。因此，本文提出一种基于短时能量突变的时延估计方法对近场火箭发射次声信号进行定位。

1 火箭发射信号时延估计

基于 TDOA的声源定位方法是利用传感器阵列各个阵元接收声信号因其传输过程中传播路径不同而引起的时间差，然后结合阵列的几何关系和时延量估计出声源的位置[16]，其包含时延估计和定位解算两部分，该定位算法的关键在于获取精确的时间差，利用时间差可以进一步确定声源的速度、方位和距离等参数。

1.1 信号特征分析

对某次三元阵采集到的火箭次声信号先进行去除大气本底压力处理后，再采用0.01～20 Hz的带通滤波器对信号进行去噪，得到的信号波形如图1所示，其中图1（a）、1（b）、1（c）分别表示1、2、3号传感器采集到的次声信号进行滤波后的波形图，横坐标表示的是次声信号到达传感器的时刻。从图1可以看到，随着时间的增加，次声信号波形逐渐增强，有几个起伏，然后慢慢消失，波形持续时间长、强度高，这符合典型的液体燃料发动机产生的次声信号波形特征[17]，与实际情况相符。

次声信号的频谱图如图2所示，此次火箭发射次声信号的频率范围是 4～12 Hz。次声源到达三个传感器的传播路径存在差异，在传播过程中的衰减程度不一致，且火箭发射次声源一直处于运动状态，多普勒效应会引起次声信号的幅值和相位的偏移，每个传感器接收到的次声信号的频带范围存在差异。

图2 相应的次声信号的频谱图Fig.2 The corresponding spectrum diagrams of the infrasound signals

由于火箭次声信号持续时间长达3～4 min，且声源位置相对于传感器阵列位置一直在发生变化，继续采用峰值检测或者广义互相关法求时间差的方法误差较大，所以本文提出一种基于信号短时能量突变的时延估计算法，首先获得接收信号的各个传感器之间的时间差，再对其进行位置解算。

1.2 火箭发射次声信号短时能量突变

当火箭发射产生的次声信号到达各个传感器时，传感器输出的次声信号能量会发生突变。一般地，无次声事件发生时，次声传感器接收到的信号能量幅度较低且波动起伏不大。因此，当传感器接收到火箭发射产生的次声信号时，传感器输出的信号能量幅值会发生突变且幅值较大，即火箭次声信号到达传感器的时间点就是信号能量发生突变的位置。实际上，传感器接收到的次声信号不仅包含火箭发射产生的次声信号，还包含其他噪声信号，如风噪声等，有用信号与噪声信号的叠加导致不能准确地确定火箭发射次声信号能量突变的位置，即不能准确确定次声信号到达传感器的时间点。次声探测站的传感器接收到的信号在一定程度上是经过玫瑰状降噪管过滤了一部分噪声后的信号，因此虽然还有其他无用信号的存在，但仍然可以认为火箭发射次声信号到达各个传感器的时间大约在信号能量发生突变的位置。信号能量发生突变位置的附近区域相对更容易找到，因此首先可以先找到包含信号到达传感器时间点的该区域。火箭发射产生的次声信号持续时间较长，通过对次声信号进行加窗计算，计算窗内的信号能量（即短时能量），先找出次声信号能量发生突变的区域[18]。

设第i个时间窗内次声信号的短时能量计算公式[18]为

式中，s（·）表示火箭发射产生的次声信号；w（·）为窗函数，此处窗函数选为矩形窗。本文采集到的火箭发射次声信号数据较长，因此选取窗长为1 000。定义矩形窗函数的窗长为N，其函数表达式为

根据式（2）计算得到次声信号的短时能量序列{Ei,i=0,1,2,…}后，可以大致看出信号的能量变化情况。假设次声信号首先在第I个时间窗内发生能量突变，计算该时间窗内的次声信号的瞬时能量：

图3为三个传感器接收到的火箭发射产生的次声信号根据上述过程计算的短时能量的结果。

图3 相应的次声信号的短时能量曲线Fig.3 The corresponding short-time energy curves of the infrasound signals

图4为对火箭发射次声信号的短时能量值进行相邻能量点之间计算差值后的能量突变幅度的曲线。可以明显地看出，次声信号的能量变化存在多个极值，即火箭发射过程中，有多个阶段产生了次声信号。

图4 短时能量突变幅度曲线Fig.4 The amplitude curves of short-time energy mutation

图5为传感器接收到的信号首次发生能量突变的时间窗内的次声信号能量变化幅值情况的曲线。即火箭发射次声信号到达传感器的时间点包含在该时间窗内。

图5 短时能量突变时间窗内能量变化曲线Fig.5 Energy change curves within the time window of short-time energy mutation

1.3 火箭次声信号时间差估计

设两两传感器之间的时间差为tij（i,j=1,2,3），若tij＜0，说明次声波信号先到达第i号传感器，后到达第j号传感器。理想情况下，闭环三元阵的时间差闭环和等于0，即传感器之间的延时量满足

而在实际情况中，由于噪声和大气传播路径的复杂性，闭环系统延时量之和只能趋于0[19-20]，即：

根据图5可知，信号能量突变幅度在整个信号持续时间范围内出现多个极大值，根据各路信号能量突变幅度设置峰值检测的条件获取时间差。本文设置1号传感器最小峰值高度为1.5，2号传感器最小峰值高度为1.5，3号传感器最小峰值高度为0.4，检测峰值的时间间隔为 25个采样点。可得到多个满足条件的位置，采用第一个位置为信号到达时间点，进而求解出本文中的次声事件根据上述方法得到的时间差为

计算t12+t23+t31=0，满足条件。

2 近场次声源的三站定位方法

对采集到的次声信号进行频谱分析，火箭次声波频率集中在4～12 Hz之间，相应波长λ约为27～81 m，由于声源S与阵列之间的距离d满足：

此时声源S处于近场范围，需要按照球面波来计算。

2.1 基于时间差的三站定位原理

基于时间差的次声定位技术主要是利用到达任意两个传感器之间的时间差，并结合传感器阵列之间的几何关系[21]，进而确定声源位置。

三元阵中传感器的位置关系如图6所示。根据三元阵的经纬度坐标可计算出三个传感器之间的距离：

图6 三元阵位置示意图Fig.6 The position diagram of tripartite array

建立如图7所示的次声三元阵，以其中一个阵元为原点，正北方向为Y轴正向、正东方向为X轴正向建立平面直角坐标系。本文中，A、B、C三点表示3号、1号、2号传感器所在位置。

图7 三元次声阵定位原理图Fig.7 Locating principle diagram of the ternate infrasound array

假设某一火箭发射事件发生后，三个传感器接收到次声信号的时间为：t1、t2和t3，若A点最先接收到次声信号，以B点为圆心，声速vt13为半径画圆，其中v为声速，从A点向该圆画两条切线，连接B和两个切点做法线，这两个法线方向即为可能的来波方向；同理，以C点为圆心，vt23为半径画圆，从A点向该圆做两条切线，连接C和两个切点做法线，这两个法线方向即为可能的另外两个来波方向；两两来波方向相交于四个点S1、S2、S3和S4，这4个点即为可能的发射点。

根据式（8）的阵元间距，再结合正弦定理即可得到三元阵的夹角以及相对于正北方向的夹角：∠ABC、∠ACB、∠BAC、LAB、LBC、LAC、α和β。假如次声波最先到达3号传感器即A点，A点和B点的时间差获得来波方向为BS1和BS4，A点和C点间的时间差获得来波方向为CS1和CS4。AB的来波方向和AC的来波方向相交得到4个点：S1、S2、S3和S4，这四个点即为可能的发射点（v为声速）：

由此可得S4的方位角为：φ=-1 80°+∠S4AB+α。此处φ表示目标点与正北方向的夹角，设γ为目标点与正东方向的夹角，得到最终的方位角为：γ=90°-φ。

同理，可得到其余可能的发射点的方位角和距离。从实际工程应用方面考虑，可排除在阵内的发射点S1。再将得到的其余三个可能发射点到三个阵元的距离长短与次声到达时间进行逻辑判断，若某发射点到1号传感器的距离大于到2号传感器的距离，而t1＜t2，从逻辑上排除该发射点。从排除后剩下的点中计算发射点声波传输到B点和C点的实测波程差与计算得到的波程差中选择误差最小的，即为实际发射点。

2.2 定位算法实现框图

采用编程语言进行定位算法的实现，框图如图8所示。

图8 三元次声阵的定位算法流程图Fig.8 Flowchart of locating algorithm of the ternate infrasound array

2.3 定位结果

采用第1节得出的时延量代入上述的定位方法进行试验，根据最小波程差得到声速 v=325.5 m·s-1，进行定位计算距离A点为54.412 km，与正东方向的夹角为118.984 3°。与实际位置对比，定向误差在2°以内，距离误差约为3.17%。

3 结 论

从火箭发射次声信号波形看出，火箭发射次声信号持续时间较其他次声事件波形较长，且由于火箭发射是一种运动次声源，采用峰值检测或广义互相关等方法进行时延估计误差较大。因此本文采用了一种基于短时能量突变进行时延估计，并采用近场的三站定位方法进行定位解算，求解结果与实际情况方位角相差在2°以内，距离误差约为3.17%，定位精度较高。