基于单水听器的匀速直线运动直升机三维参数估计方法

1.本发明涉及水下探空领域，具体涉及一种基于单水听器的匀速直线运动直升机三维参数估计方法。


背景技术：

2.反潜直升机对于水下潜艇具有很大威胁，为了提高水下潜艇对抗反潜机的能力，有必要尽快解决水下对空探测的技术问题。航空反潜多采用螺旋桨飞机或直升机，它们的辐射噪声是由离散谐波信号与宽带连续谱信号组成的，其中在15hz—200hz低频段内离散线谱的能量最为显著。线谱信号具有频率较低、功率集中、稳定性强的特点，可以传输较远距离，是水听器检测和识别空中目标的主要信息载体。当反潜机以一定高度匀速通过水面上方时，水听器接收的信号将会产生多普勒频率，频率变化曲线与目标固有频率、速度、距离等参数信息有关，这为目标运动参数的估计提供了可能。近年来，在水听器估计空中声源运动信息问题上，学者们已经做出了大量研究。buckingham(m.j.buckingham,e.m.giddens,f.simonet,t.r.hahn.propeller noise from a light aircraft for low-frequency measurements of the speed of sound in a marine sediment.journal of computational acoustics.2002；10(04):445—64.)实测直升机过顶飞行数据，分别获取了空气、海水和海底沉积层中的音频信号，实测数据表明过顶信号具有时变特性，验证了在水中可以有效接收到空气中声源产生的多普勒频率信号，为水下探空提供可能性。随后，学者们对水下探空的建模分析和理论研究主要在单个水听器或水听器阵列两个领域。在单个水听器方面，ferguson(b.g.ferguson.a ground-based narrow-band passive acoustic technique for estimating the altitude and speed of a propeller-driven aircraft.the journal of the acoustical society of america.1992；92(3):1403—1407.)使用射线理论分析了空气—水两种介质的声传播模型，并根据二维域的多普勒频率公式推导出空中目标过顶飞行参数的方法。随着水声探测设备的不断进步，具有“∞”型指向性的矢量水听器被应用于探测空中目标。与传统的标量水听器相比，矢量水听器不仅能获得目标信号标量场信息(如声压)，还可获得矢量场信息(如振速)。采用合理的数学模型对矢量信息进行处理，即可确定目标方位和其他飞行参数。在水听器阵列方面，ferguson和lo(k.w.lo,b.g.ferguson.broadband passive acoustic technique for target motion parameter estimation.ieee transactions on aerospace and electronic systems.2000jan；36(1):163—75.)采用宽孔径水听器阵列获取空中目标的方位信息，对空气中飞行的直升机进行参数估计和定位，实现了目标运动信息由二维域到三维域的扩展。彭朝晖等(刘凯悦,彭朝晖,张灵珊,王光旭.水下水平阵对空中运动声源的线谱探测.声学学报,2019；.44(4):567—575.)针对水下探测空中运动声源的干扰问题，提出一种稳健可靠的基于水听器阵列的线谱检测方法。在这些研究中，基于传统单水听器的水下探空算法仅能估计二维平面上的声源飞行参数问题，即假定直升机刚好从水听器顶部飞过，而没有考虑三维域上的偏航距离等信息。单矢量水听器可以估计方位信息，但存在指向性盲区，
实际应用中也更易受背景噪声干扰。而水听器阵列虽然能够估计出方位信息，但成本较高、布置不便，不适用于小体积的水下探测平台。


技术实现要素：

3.针对在水下估计反潜机的运动参数场景应用局限、布置成本较高的问题，本发明从两层等声速介质的声传播几何模型出发，利用直升机的螺旋桨谐波频率来计算水下声学的多普勒频率曲线。从三维空间几何关系上推导了直升机航行参数，包括固有频率、速度、高度、偏航距离，与多普勒之间的非线性函数关系。利用多普勒频率曲线不对称性，从拟合的多普勒曲线以及一阶和二阶导数上可计算上述参数。在不采用水听器阵列的前提下，该方法实现了对空中匀速直航目标的全部飞行参数解算。
4.一种基于单水听器的匀速直线运动直升机三维参数估计方法，该方法包括以下步骤：
5.步骤一：建立水听器接收匀速直线运动声源声波的几何模型；
6.定义静止的水听器节点记为点h，位于水下深度d处；恒定频率为f0的直升机辐射声源记为点s，它以恒定亚音速v、恒定高度h飞过点h斜上方，v《ca《cw，ca和cw分别是声波在空气和水介质中的传播速度；投射到静止节点h所在的水平面上的动声源记为s'，动点s'和节点h之间水平距离记为w(t)，即点声源s和节点h直线距离r(t)的水平投影；点声源s飞至距节点h的最近点处时，二者距离r
min
(t)的水平投影记为最短偏航距离w
min
；是水下视线角；空气介质和水介质的声波速度分别为ca、cw，从空气介质到水介质中的折射角分别为θi(t)、θ
t
(t)，折射率为
7.步骤二：根据声源与水听器相对运动的几何规律，推导出在t时刻多普勒频率fd(t)与声源频率f0、速度v、高度h、最短偏航距离w
min
之间的函数关系；
[0008][0009]
步骤三：利用水听器采集的多普勒频率信息，得到声源的频率和速度；利用多普勒频率fd(t)与声源频率f0、速度v的函数关系，通过下式计算点声源的固有频率f0和速度v
[0010][0011]
其中，fd(-∞)和fd( ∞)分别为近似的水听器与直升机水平距离大于二倍水听器
与直升机垂直距离的时刻t-2
与t
2
的多普勒频率；
[0012]
步骤四：利用步骤二中的多普勒频率fd(t)与高度h、最短偏航距离w
min
的函数关系，通过下式得到直升机辐射声源的高度和偏航距离：
[0013][0014]
其中，τ0是f
′d(t)极小值对应的时刻，τ-1
是f
″d(t)极小值对应的时刻。
[0015]
进一步地，所述步骤四中多普勒频率fd(t)的一阶导数为：
[0016][0017]
其中，
[0018]
这里，
[0019][0020][0021][0022][0023][0024][0025][0026][0027]
[0028][0029][0030][0031]
进一步地，多普勒频率fd(t)的二阶导数为
[0032][0033]
其中
[0034][0035][0036][0037][0038][0039]
与
[0040][0041]
其中，
[0042]
[0043][0044][0045][0046][0047][0048][0049][0050][0051][0052][0053][0054][0055][0056][0057]
[0058]
本发明针对小样本调频信号数据和未知概率分布噪声场景，设计了基于单水听器的匀速直线运动直升机三维参数估计技术，推导出一种准确性和稳定性水下探空技术。本发明与现有的水下探空方法相比具有如下有益效果：
[0059]
(1)本发明的方法能够适用于小样本调频信号数据和未知概率分布噪声场景；
[0060]
(2)针对在水下估计反潜机的运动参数问题，本发明的方法仅使用单水听器即能实现匀速圆弧运动声源的三维空间信息估计；
[0061]
(3)本发明的方法准确性高，稳定性强。
附图说明
[0062]
图1为空中沿直线运动的点声源与静止水听器节点的三维图。
[0063]
图2为声波传播路径的俯视图。
[0064]
图3为声波传播路径的侧视图。
[0065]
图4为多普勒频率曲线及其导数的理论结果。
[0066]
图5为多普勒频率曲线的导数。
具体实施方式
[0067]
下面根据附图和优选实施例详细描述本发明，本发明的目的和效果将变得更加明白，应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
[0068]
步骤一：建立水听器接收匀速直线运动声源声波的几何模型。
[0069]
图1是直升机声源沿定高直线ls飞过水听器节点的三维示意图。静止的水听器节点记为点h，位于水下深度d处。恒定频率为f0的直升机辐射声源记为点s，它以恒定亚音速v、恒定高度h飞过点h斜上方，v《ca《cw，ca和cw分别是声波在空气和水介质中的传播速度。投射到静止节点h所在的水平面上的动声源记为s'，该投影动点s'形成的直线轨迹记为ls'。动点s'和节点h之间水平距离记为w(t)，即点声源s和节点h直线距离r(t)的水平投影。声线在空气介质和水介质中传播距离分别为ra(t)、rw(t)。特别地，点声源s飞至距节点h的最近点(closest point of approach，cpa)处时，二者距离r
min
(t)的水平投影记为w
min
；是水下视线角，α(t)为偏向角，即直线轨迹ls与直线s'h的夹角，β(t)是空气中的视线角，与入射角θi(t)互为余角；从空气介质到水介质中的折射角分别为θi(t)、θ
t
(t)，折射率为
[0070]
点声源s的速度v可以在飞行水平面分解为两个相互垂直的分量，其中一个分量v1(t)与水听器处于同一垂面，其大小为
[0071]
v1(t)＝vcosα(t)
ꢀꢀꢀ
(1)
[0072]
进一步地，v1(t)可以分解为两个相互垂直的分量，其中分量v
11
(t)与入射声射线方向一致，其大小为
[0073]v11
(t)＝v1(t)cosβ(t)
ꢀꢀꢀ
(2)
[0074]
由于点声源s和水听器节点h相对速度随时间发生变化，因此水听器接收信号的瞬时频率也随时间发生变化。显然，水听器所接收到的多普勒频率曲线fd(t)由参数{f0,v,β
(t),α(t)}决定。对于固有频率为f0的动点声源s，水听器在t时刻所接收到的多普勒频率为
[0075][0076]
公式(3)给出多普勒频率fd(t)与声源运动参数之间的关系，下文利用点声源s和水听器节点h的几何关系，可以获得偏向角α(t)和空气中视线角β(t)。
[0077]
步骤二：根据声源与水听器相对运动的几何规律，推导出在不同时刻t多普勒频率fd与声源频率f0、速度v、高度h、最短偏航距离w
min
之间的函数关系；
[0078]
首先，根据图2所示的俯视图中动声源s与水听器节点h的位置关系，可以求解偏向角α(t)。设动点声源s在cpa的左侧或右侧时，速度v和时刻t为负或正，动点声源s和静止节点h之间的俯视位置关系如图3所示。根据勾股定理和三角函数，在时间t时刻，偏向角为
[0079][0080]
然后，根据图3所示的声线传播路径几何关系，可以求解空气中视线角β(t)。图3是动点声源s和静止节点h的侧视图，根据声源与水听器的几何关系，任意位置处的点声源有且仅有一条入射角为θi(t)、折射角为θ
t
(t)的声射线传播至水听器。
[0081]
在以∠sth为顶角δsth中，由余弦定理可知入射角θi(t)和折射角θ
t
(t)满足如下关系，
[0082][0083]
其中，声波在空气介质与水介质的传播距离分别为ra(t)和rw(t)，可写为
[0084][0085]
点声源s和水听器节点h的直线距离r(t)为
[0086][0087]
入射角θi(t)和折射角θ
t
(t)之间满足snell定理，且与空气视线角β(t)的互补关系和水下视线角的互补关系可表示为
[0088][0089]
将式(8)代入(5)后，整理可得水下视线角与点声源高度h、声源速度v、偏航距离w
min
、水听器深度d之间的关系
[0090][0091]
将式(4)、(8)和(9)代入式(3)中，可以消去水平角α(t)、空气视线角β(t)、折射角θi(t)和折射角θ
t
(t)，化简后可以得到方程组
[0092][0093]
对于方程组(10)，水下视线角作为可消参数，因此可以确定多普勒频率fd(t)与声源运动参数{f0,v,d,h,w
min
}之间的函数关系。
[0094]
步骤三：利用水听器采集的多普勒频率信息，得到声源的频率和速度。
[0095]
在建立直线飞行轨迹模型和分析三维多普勒频率曲线基础上，本发明提出一种可应用于三维空间内的单水听器估计飞行参数方法。若已知空气介质声速ca、水介质的声速cw、水听器深度d，本方法可依据式(10)使用单水听器所接收的多普勒频率曲线fd(t)估计直升机的固有频率f0、速度v、飞行高度h和最短偏航距离w
min
。
[0096]
声源s行至当前位置所处的时刻为
[0097][0098]
由表1可知，声源s在经过cpa点前后，cosα(t)分别为负值或正值，显然时刻t也是负或正。对于如式(10)所示的多普勒频率曲线，时刻的t正负性会导致空气中动声源在接近或远离水听器过程中所产生的多普勒频率曲线具有不对称性。
[0099]
表1各参数在穿过cpa点前后的变化规律
[0100][0101]
图4展示了在高度h＝151m处谐振频率为68hz的声源以速度v＝100m/s经过cpa点
前后，位于水下90m处水听器所接收多普勒频率曲线的理论值及其一阶、二阶导数。时刻t正负性造成了频率偏移量的不对称，声源s距声接收器h距离相同时，接近状态时的频率偏移量比远离状态时的频率偏移量大。
[0102]
根据表1可知，声源位置处于无穷时，此时式(10)所示的多普勒频率方程退化为
[0103][0104]
利用曲线因时刻t正负性带来的不对称以及声源处于无穷时多普勒频率方程，计算点声源的固有频率f0和速度v
[0105][0106]
设置固有频率f0对应的时刻t＝0。接着，联立多普勒频率fd(t)的一阶导数和二阶导数，并分别将极点值为零的时刻代入
[0107][0108]
其中，τ0是fd'(t)极小值对应的时刻，τ-1
是f
d”(t)极小值对应的时刻。于是，可以获得点声源的高度h和最短水平投影w
min
。
[0109]
式(14)中关于多普勒频率fd(t)的一阶导数和二阶导数的求解过程及结果具体如下：
[0110]
多普勒频率fd(t)的一阶导数为
[0111][0112]
其中，是水下视线角关于时间t的一阶导数。对式(9)所示的隐函数求导，整理后可得
[0113][0114]
这里，
[0115]
[0116][0117][0118][0119][0120][0121][0122][0123][0124][0125][0126][0127]
将式(a.2)代入式(a.1)后，可得多普勒频率fd(t)的一阶导数。
[0128]
接着，多普勒频率fd(t)的二阶导数为
[0129][0130]
其中
[0131][0132]
[0133][0134][0135][0136]
对式(a.2)所示的求导，整理后可得
[0137][0138]
这里，
[0139][0140][0141][0142][0143][0144][0145][0146][0147]
[0148][0149][0150][0151][0152][0153][0154][0155]
将式(a.2)、(a.4)代入式(a.3)后，可得多普勒频率fd(t)的二阶导数。
[0156]
美国声学学会提供了一段由水下90米处水听器记录的直升机通过水听器上方，总计120秒的音频数据，采样频率为44100hz。直升机在一定高度沿直线匀速航行，假定声波在空气和水中的速度为340m/s和1500m/s。根据单水听器采集的多普勒频率信号及其深度信息，使用本发明所提出的水下探空算法解算空中匀速直航直升机的全部飞行参数。
[0157]
根据本发明提出的基于单水听器的匀速直线运动估计直升机三维参数方法，可以实现最小偏航距离w
min
的估算。对于累积差分功率的时频估计结果，如图5所示的多普勒频率曲线，可近似认为f(-∞)＝73.17hz、f( ∞)＝64.37hz。于是，由式(13)可计算得
[0158][0159]
接着分别计算其一阶导数与二阶导数如图5所示，由图可知，一阶导数的极小值点在t0＝54.14sec处，二阶导数的极小值点在t-1
＝53.25sec处。图5中，一阶导数、二阶导数的横轴均为时间，纵轴分别代表多普勒频率随时间变换率以及多普勒频率随时间变换率的变换率，以此揭示多普勒频率随时间的变化关系。
[0160]
在传统二维方法中，可使用fd(t)的近似一阶导数式求解直升机航行高度，即
[0161]
[0162]
故直升机高度h为132米。而本发明提出的直升机三维参数估计方法通过联立fd(t)的一阶导数和二阶导数可得直升机航行高度与偏航距离，即
[0163][0164]
其中，τ0＝0秒，τ-1
＝t-1-t0＝-0.89秒。可得直升机高度h＝149.35m，最小偏航距离w
min
＝1.41m。
[0165]
使用窗内补零短时傅里叶变换的瞬时频率估计算法和无网格瞬时频率估计算法得出直升机螺旋桨的多普勒频移曲线，然后根据式(14)和(15)所示的二维参数估计算法和式(14)和(16)所示的三维参数估计算法可分别得到匀速运动直升机的飞行参数。按照美国声学学会公布的结果(www.tridipassio.com/asa2019)，直升机固有频率为68hz、速度约为100m/s、航行高度约为150m、偏航距离约为0m。比较传统二维参数估计方法和本发明的三维参数估计方法，后者更接近合理数值，因此本发明提出单水听器估计匀速直线运动直升机三维参数算法具备有效性。
[0166]
本领域普通技术人员可以理解，以上所述仅为发明的优选实例而已，并不用于限制发明，尽管参照前述实例对发明进行了详细的说明，对于本领域的技术人员来说，其依然可以对前述各实例记载的技术方案进行修改，或者对其中部分技术特征进行等同替换。凡在发明的精神和原则之内，所做的修改、等同替换等均应包含在发明的保护范围之内。