一种多元阵次声时延估计方法与流程

1.本发明涉及一种多元阵次声时延估计方法，属于声学定位技术领域，用于如火箭发射、 泥石流、火山爆发、地震等持续时间较长的次声事件定位或定向。


背景技术：

2.在现在的次声定位方法中，基于广义互相关时延估计、单阵定向多阵定位的次声定位技 术得到了广泛研究，针对阵元间时间差的时延量获取提出了大量的改进方法。
3.火箭发射事件的时间与地点明确，是核爆次声探测中一类具有重要参考研究价值的次声 源，对于次声时延估计、定位算法等具有重要参考意义。火箭发射是一种运动次声源，其产 生的次声波信号主要由火箭飞行时激发空气产生，因此采集到的火箭发射次声信号持续时间 可长达3～4分钟，甚至更久。由于次声波的远距离传播特性，大多数学者考虑的是次声源的 远程定位，且火箭次声信号的持续时间较长，声源并非静止的，难以识别判定不同发射阶段 所产生的次声信号，采用常用的广义互相关时延估计定位技术时间差的计算精度较差，依据 广义互相关法进行时延估计并定位结算时，甚至出现南辕北辙的事情，难以实现对次声事件 的定位。


技术实现要素：

4.本发明的目的是解决次声事件定位中阵元间时延估计不准导致定位不准的问题，提供一 种多元阵次声时延估计方法，通过该方法可以较为准确地确定火箭发射点次声源的来波方向， 并采用三元次声阵实现对火箭发射点的定位。
5.本发明解决上述问题采用的技术方案：多元阵次声阵时延估计方法步骤如下：
6.步骤一、通过对次声传感器采集到的信号先去除大气本底压力，进行[0.1 12]hz带通或 者12hz低通滤波处理后得到的次声信号波形，对次声信号波形进行加矩形窗处理，矩形窗 长度n＝1000；
[0007]
步骤二、计算第i个时间窗内的次声信号短时能量：
[0008][0009]
s(
·
)表示次声信号，w(
·
)为矩形窗函数；
[0010]
步骤三、通过步骤二计算得到第i个时间窗内的短时能量值，将短时能量值的后一个值 减去前一个值的相邻瞬时能量点之间的差值绘制成能量突变幅度曲线；设置最小峰值检测高 度为噪声峰峰值的10倍和最小峰值检测距离为采样周期25倍记录第一个峰值；记录第一个 峰值的采样时间为传感器感应时间，将阵内传感器的能量突变幅度曲线的第一个峰值时间两 两相减获得次声信号到达阵元间的两两时延量。
[0011]
步骤四、以次声监测阵中一个阵元所在位置为原点，正北方向为y轴正向、正东方向为 x轴正向建立平面直角坐标系；
[0012]
步骤五、根据两点的经纬度坐标a(longa，lata)、b(longb，latb)按照
[0013][0014]
计算出相邻阵元间距，再运用三角正弦定理即得到三元阵δabc的夹角以及相对于正北方 向的夹角；当次声波最先到达a点，以a点和b点的时间差获得两个来波方向，a点和c 点间的时间差获得两个来波方向；a点、b点的来波方向和a点、c点的来波方向相交得到 4个点，这四个点即为发射点；
[0015]
步骤六、根据实际工程应用和次声信号到达传感器的先后顺序，排除在阵内的发射点， 再按照实测波程差与步骤三得到的时延量和次声速度计算得到的波程差中误差最小的进行选 择，选择误差最小的点为实际发射点；
[0016]
步骤七、当次声波先到达b点，以b点和a点的时间差得到两个来波方向，以b点和 c点的时间差得到两个来波方向；当次声波先到达c点，以c点和a点的时间差得到两个来 波方向，以c点和b点的时间差得到两个来波方向，分别得到四个发射点，根据步骤六得到 当b点或c点最先接收到次声信号时的定位算法；
[0017]
步骤八、将短时能量突变得到的时间差作为时延量输入到步骤四至步骤七中，得到次声 事件的来波方向和次声传感器与发射点之间的距离，实现对次声事件的定位。
[0018]
对次声信号进行滤波处理，滤波器的设计为[0.01 20]hz的带通滤波器或20hz的低通滤 波器，带通滤波器的低端变化范围为0.001～10hz、高端为2～100hz，带通滤波器的通带截 止频率高至100hz。
[0019]
对去噪后的次声波波形进行加矩形窗处理，n选择范围为10～10000。
[0020]
最小峰值检测高度为噪声峰峰值的3～100倍，最小峰值检测距离为采样周期10～100 倍。
[0021]
延估计方法的处理对象是火箭发射、中远区核爆探测、泥石流的次声事件。
[0022]
本发明的有益效果：对持续时间较长或由一系列连续事件组成的次声事件进行定位时， 采用目前常用的互相关时延估计算法，通常能得到南辕北辙的定向或定位效果，该算法采用 能量突变的方法进行时延估计，可由采用常用的基本相关性算法计算时延量误差从25～120％ 提高至6％以内，从而实现次声定位距离误差约3％左右、定向误差小于2
°
，大大提高定向和 定位精度。
附图说明
[0023]
图1 卫星发射事件滤波后的三元阵次声信号
[0024]
图中：a为1号站滤波波形，b为2号站滤波波形，c为3号站滤波波形；
[0025]
纵坐标为声压，单位为pa，横坐标为时间。
[0026]
图2 三路传感器接收到的火箭次声信号计算短时能量的结果
[0027]
图中：细实线为1号站短时能量曲线，虚线为2号站短时能量曲线，粗实线为3号站短 时能量曲线；
[0028]
纵坐标为短时能量，单位为j，横坐标为时间，单位为s。
[0029]
图3 火箭次声信号短时能量突变幅度的曲线
[0030]
图中：细实线为1号站短时能量突变曲线，虚线为2号站短时能量突变曲线，粗实线为 3号站短时能量突变曲线；
[0031]
纵坐标为短时能量突变，单位为j，横坐标为时间，单位为s。
[0032]
图4 短时能量突变时间窗内能量变化曲线
[0033]
图中：细实线为1号站短时能量突变曲线，虚线为2号站短时能量突变曲线，粗实线为 3号站短时能量突变曲线；
[0034]
纵坐标为短时能量突变，单位为j，横坐标为时间，单位为s。
[0035]
图5 三元次声阵坐标系示意图
[0036]
图6 时延估计方法进行定位算法流程图
具体实施方式
[0037]
下面结合附图和实施例对本发明作进一步描述。
[0038]
实例1
[0039]
本发明的可用于次声源定位的多元阵时延估计方法，所述的次声源定位是通过短时能量 突变进行时延估计，再利用时延量结合多元阵的几何关系实现对次声源的定位。
[0040]
以某次卫星发射事件次声信号为例，所述的短时能量突变进行时延估计的具体步骤如下：
[0041]
步骤1)通过对传感器采集到的信号进行去噪、滤波处理后得到的三元次声信号波形如 图1所示，对火箭发射次声信号进行加矩形窗处理，定义矩形窗的窗长为n，其函数表达式 为
[0042][0043]
本实例选择的矩形窗窗长n＝1000。
[0044]
步骤2)计算第i个时间窗内的次声信号短时能量：
[0045][0046]
其中，s(
·
)为火箭发射次声信号；w(
·
)为矩形窗函数。
[0047]
三元次声传感器接收到的火箭次声信号计算短时能量的结果如图2所示。
[0048]
步骤3)计算次声信号能量突变情况，观察能量变化幅度与快慢情况。根据步骤2计算 得到的时间窗的瞬时能量，计算时间窗内的相邻瞬时能量点之间的差值，得到时间窗内的能 量突变幅度曲线。图3为火箭次声信号短时能量突变幅度的曲线，可以明显看到出现多个局 部能量突变峰值。图4为发生短时能量突变的时间窗内次声信号能量突变幅度的曲线，火箭 次声信号到达传感器时间点即包含在该时间窗内。
[0049]
根据各路传感器采集到的信号得到的能量情况，设置最小峰值检测高度约为噪声峰峰值 的10倍和最小峰值检测距离约为采样周期25倍记录第一个峰值；记录第一个峰值的采样时 间为传感器感应时间，将阵内多个传感器的能量突变幅度曲线的第一个峰值时间两两相减获 得次声信号到达阵元间的两两时延量。
[0050]
设两两传感器之间的时间差为t
ij
(i，j＝1，2，3)，当t
ij
＜0，说明次声波信号先到达第i号传 感器，后到达第j号传感器。理想情况下，闭环三元阵的时间差闭环和等于0，即传感器之间 的延时量满足
[0051]
t
12
t
23
t
31
＝0
[0052]
根据图4可知，信号能量突变幅度在整个信号持续时间范围内出现多个极大值，根
据各 路信号能量突变幅度设置峰值检测的条件获取时间差。本例设置1号传感器最小峰值高度为1.5，2号传感器最小峰值高度为1.5，3号传感器最小峰值高度为0.4，检测峰值距离为25。 可得到多个满足条件的位置，采用第一个位置为信号到达时间点，进而求解出本文中的次声 事件根据上述方法得到的时间差为：
[0053]
t
12
＝-7.05s，t
13
＝7.34s，t
23
＝14.39s
[0054]
计算t
12
t
23
t
31
＝0，满足条件。
[0055]
采用短时能量突变算法计算时延量误差如下表所示，误差值在6％以内。
[0056][0057]
而采用常用的基本互相关算法时，计算得到的时延量如下表所示，误差值在25％-120％ 之间。
[0058][0059]
所述的利用时延量结合三元阵的几何关系确定声源位置的具体步骤如下：
[0060]
步骤4)以其中一个阵元为原点，正北方向为y轴正向、正东方向为x轴正向建立平面 直角坐标系。假设某一火箭发射事件发生后，1号、2号、3号传感器接收到次声信号的时间 分别为：t1、t2和t3，声速为v，1、2号传感器的时间差为：t
21
＝t2-t1，1号、3号传感器之 间的时间差为：t
31
＝t3-t1，2号、3号传感器之间的时间差为：t
32
＝t3-t2，并设1号、2号和3 号传感器所在位置分别为a点、b点和c点。若a点最先接收到次声信号，以b点为圆心， (v*t
21
)为半径画圆，从a点向该圆画两条切线，连接b点和两个切点做法线，这两个法线方 向即为可能的来波方向；同理，以c点为圆心，(v*t
31
)为半径画圆，从a点向该圆做两条切 线，连接c点和两个切点做法线，这两个法线方向即为可能的另外两个来波方向；两两来波 方向相交于四个点s1、s2、s3和s4，这4个点即为可能的火箭发射点，如图4所示。
[0061]
步骤(5)根据所建立的三角阵的经纬度坐标计算出阵元间距，再运用正弦定理即可得到 三元阵的夹角以及相对于正北方向的夹角：∠abc、∠acb、∠bac、lab、lbc、lac、 α和β。假如次声波最先到达a点，a点和b点的时间差获得来波方向为bs2和bs4，a点 和c点间的时间差获得来波方向为cs1和cs4。ab的来波方向和ac的来波方向相交得到4 个点：s1、s2、s3和s4，这四个点即为可能的发射点：
[0062]
∠s4cb＝∠acb arccos(t
31
*v/l
ac
)
[0063]
∠s4ba＝arccos(t
21
*v)
[0064]
过点b作bd
⊥
cs4交s4c于点d，
[0065]
∠s4bd＝∠abc-(90
°‑
∠dcb) ∠s4ba＝∠abc ∠s4ba ∠s4cb-90
°
[0066][0067][0068][0069]
根据以上公式可得s4的正北方向的夹角为：设γ为目标点与正东方 向的夹角，得到最终的方位角为：
[0070]
∠s3cb＝∠s4cb＝∠4cb arccos(t
31
*v/l
ac
)
[0071]
∠s3ba＝∠s4ba＝arccos(t
21
*v/l
ab
)
[0072]
∠s3bd＝∠abc-∠s3ba-∠dbc＝∠abc-∠s3ba-90
°
∠s3cb
[0073][0074][0075][0076]
根据以上公式可得到s3的正北方向的夹角为：设γ为目标点与正东 方向的夹角，得到最终的方位角为：
[0077]
∠s2ba＝∠s4ba＝arccos(t
21
*v/l
ab
)
[0078]
∠s2ca＝arccos(t
31
*v/l
ac
)
[0079]
∠s2cb＝∠acb-∠s2ca
[0080]
延长s2b，过点c作ce
⊥
s2b于点e，
[0081]
∠cbe＝180
°‑
∠abc-∠s2ba
[0082]
l
ce
＝l
bc
*sin(∠cbe)
[0083]
∠bce＝90
°‑
∠cbe
[0084]
l
be
＝l
bc
*sin(∠bce)
[0085]
∠s2ce＝∠s2cb ∠bce
[0086][0087][0088][0089][0090]
根据以上公式可得到s2的正北方向的夹角为：设γ为目标点与正东 方向的夹角，得到最终的方位角为：
[0091]
∠s1ba＝∠s3ba
[0092]
∠s1bc＝∠abc-∠s1ba
[0093]
过点c作cf
⊥
s1b于点f，
[0094]
∠bcf＝90
°‑
∠s1bc
[0095]
l
cf
＝l
bc
*cos(∠bcf)
[0096]
∠s1cf＝∠s2cb-∠bcf
[0097]
∠s1ca＝∠s2ca
[0098][0099]
[0100][0101]
根据以上公式可得到s1的正北方向的夹角为：设γ为目标点与正东 方向的夹角，得到最终的方位角为：
[0102]
步骤(6)从实际工程应用和次声信号到达传感器的先后顺序两方面考虑，可进一步排除 在阵内的发射点以及距离和次声到达时间不合逻辑的发射点s1。将s2、s3和s4按照实测波 程差与计算得到的波程差中误差最小的进行选择，即误差最小的点为计算得到的发射点。
[0103]
步骤(7)当次声波先到达b点，以b点和a点的时间差得到两个来波方向，以b点和 c点的时间差得到两个来波方向；当次声波先到达c点，以c点和a点的时间差得到两个来 波方向，以c点和b点的时间差得到两个来波方向，可分别得到四个发射点，根据步骤六得 到当b点或c点最先接收到次声信号时计算得到的发射点。
[0104]
步骤(8)将步骤3采用短时能量突变算法得到的阵元间时延量，输入到步骤5的三元阵 定位算法中，算出s2、s3和s4的点位置，通过步骤6得到波程差最小的点，即为火箭发射点， 从发射点到次声传感器的布设点的方向为来波方向，距离为发射点到次声传感器布设点的距 离，最终实现对火箭发射事件的定位。
[0105]
以上步骤流程如图6所示，采用以上步骤进行定位计算，某次卫星发射事件的次声传感 器理卫星发射塔距离为54.412km，与实际来波方向与正东方向的夹角为118.9843
°
。与实际位 置对比，采用短时能量突变算法计算的来波方向为120.1655
°
，距离为52.7399km，定向误差 在2
°
以内，距离误差为3.17％。如下表所示。
[0106]