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

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


背景技术：

2.潜艇利用自身水声接收设备来探测空中旋转翼发动机的辐射噪声，以实现反潜直升机的识别任务。urick
1.用射线法分析了空气中声源正下方附近的水下声场，并指出空气中声波以四种途径到达水听器：直接折射波、海底海面反射波、渐逝波以及海水散射波。这说明了在水中可以有效接收到空气中声源产生的多普勒频移信号，这满足了水下探空的前提条件。目前现有理论模型仅能估计直升机沿直线航行的飞行参数。例如在二维域方面，即空中目标过顶直线飞行模型，已有大量文献分析了空气—水两种介质的声传播模型，并推导出该场景声源运动参数的方法。在三维域方面，一种方法是使用单矢量水听器和恰当数学模型来估计目标方位和其他飞行参数；另一种方法是采用宽孔径水听器阵列获取空中目标的方位信息，对空气中飞行的直升机进行参数估计和定位。然而，相比直升机直线航行，沿圆弧航行时直升机与声源之间距离函数更为复杂，导致以上算法在估计圆弧航行的飞行参数时会失效。


技术实现要素：

3.针对空中匀速圆弧运动声源所激发的水声信号特征，本发明采用单水听器解决了空中快速飞行目标状态三维参数估计问题。本发明首先以直升机离散线谱为声源特征，在空气—水介质中建立声源线谱特征在匀速圆弧运动下三维多普勒传播模型。然后根据多普勒频移曲线、声源运动模型以及声线传播几何关系，选取三个时间观测点计算目标多普勒频移值，推导了单水听器估计空中匀速圆弧运动声源的三维参数估计算法。最后，利用提出的瞬时频率估计算法获取精确地多普勒频移曲线以及匀速圆弧场景下的水下探空算法，即可获取声源飞行参数。
4.本发明的目的通过如下的技术方案来实现：
5.一种基于单水听器的匀速圆弧运动直升机三维参数估计方法，包括以下过程：
6.步骤一：建立水听器接收匀速圆周运动声源声波的几何模型；
7.直升机声源在水听器节点上空以
⊙
o为圆心做匀速圆周运动，静止的水听器节点记为点 h，位于水下深度d处；恒定频率为f0的直升机辐射声源记为点s，它以恒定亚音速v、恒定高度飞过点h上方附近，其中，v＜ca＜cw，ca和cw分别是声波在空气和水介质中的传播速度；投射到静止节点h所在的水平面上的动声源记为s
′
，该投影动点s
′
形成的直线轨迹记为
⊙
o'；动点s'和节点h之间距离记为ρ(t)，即点声源s和节点h直线距离r(t)的水平投影，ρ(t)随动点s'的位置改变而变化；是水下视线角；从空气介质到水介质中的折射h角分
别为θi(t)、θ
t
(t)，折射率为点声源s的速度v在飞行水平面上分解为两个相互垂直的分量，一个分量v1(t)指向水听器方向，其大小为
8.v1(t)＝vcosα(t)
9.其中，α(t)为偏向角，是轨迹在当前位置处切线与直线s'h的夹角；
10.v1(t)在垂面上分解为两个相互垂直的分量，一个分量v
11
(t)与入射声射线方向一致，其大小为
11.v
11
(t)＝v1(t)cosβ(t)
12.其中，空气中视线角β(t)与入射角θi(t)互为余角。
13.步骤二：根据声源与水听器相对运动的几何规律，推导出在不同时刻t多普勒频移fd与声源频率f0、速度v、高度h、偏航距离w
min
之间的函数关系；
[0014][0015]
步骤三，应用声源和水听器的几何关系、多普勒频移曲线，估计出声源的运动参数。
[0016]
进一步地，所述步骤三中，
[0017]
声源s在时刻ta、tb和tc依次经过不共线的三点a、b、c，且相邻两点时间差为δt
ab
＝t
b-ta、δt
bc
＝t
c-tb；声源s在点a或b发出的声波传播至水听器h具有一定时间差δt
ab
或δt
bc
，水听器h接收到声波信号时，声源s已运动至点b或c处，产生多普勒频移；点m为弧线的中点，点{a,b,c,m}的坐标分别为
[0018][0019][0020][0021][0022]
其中，
[0023][0024]
以圆心点o为顶点、圆弧上的任二点为底边所构成的三角形应用余弦定理，有
[0025][0026]
对于∠poq，有如下关系式
[0027][0028]
以节点投影h'为顶点、圆弧上的任二点为底边所构成的三角形应用余弦定理，有
[0029][0030]
对于∠ph'q，有如下关系式
[0031][0032]
对于ρ
p
，有如下关系式
[0033][0034]
以节点投影o为顶点、h'p为底边所构成的三角形应用余弦定理，有
[0035][0036]
声源位于p∈{a,b,c,m}时，以折射点t
p
为顶点、声源s与水听器h连线为底边所构成的三角形，应用余弦定理，有
[0037][0038]
此外，声源位于p∈{a,b,c,m}时，声线由声源s传至水听器节点h的传播时间为
[0039][0040]
根据弧长公式，声源s一定时间内经过的弧长与其弧度有如下关系
[0041]
v(τ
b-τm)＝rσ
[0042]
v(τ
a-τm)＝rω
[0043]
声源位于p∈{a,b,c,m}时，其瞬时多普勒频率为
[0044][0045]
联立以上共58个方程，将x
p
、y
p
、ρ
p
、β
p
、τ
p
、r、ω、ω'、ω”、σ、σ'、υ、υ'因子视为可消参数，求解出直升机的固有频率f0、速度v、飞行高度h。
[0046]
本发明的有益效果如下：
[0047]
(1)本发明的方法能够适用于小样本调频信号数据和未知概率分布噪声场景。
[0048]
(2)针对在水下估计反潜机的运动参数问题，本发明的方法仅使用单水听器即能实现匀速圆弧运动声源的三维空间信息估计。
[0049]
(3)本发明的方法准确性高，稳定性强。
附图说明
[0050]
图1为空中沿曲线运动的点声源与静止水听器节点的三维图。
[0051]
图2为空中动点声源与静止水听器节点的俯视图。
[0052]
图3为声波传播路径的垂面图。
[0053]
图4为声波传播路径的几何关系图，其中，图(a)为点b位于弧线内，图(b)为点b位于弧线外。
具体实施方式
[0054]
下面根据附图和优选实施例详细描述本发明，本发明的目的和效果将变得更加明白，应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
[0055]
步骤一：建立水听器接收匀速圆周运动声源声波的几何模型
[0056]
在空气—水介质环境中，如图1所示，直升机声源在水听器节点上空以
⊙
o为圆心、r 为半径做匀速圆弧运动。静止的水听器节点记为点h，位于水下深度d处。恒定频率为f0的直升机辐射声源记为点s，它以恒定亚音速v、恒定高度h飞过点h上方附近，其中，v＜ca＜cw，ca和cw分别是声波在空气和水介质中的传播速度；
[0057]
投射到静止节点h所在的水平面上的动声源记为s'，该投影动点s'形成的轨迹记为
⊙
o'。动点s'和节点h之间距离记为ρ(t)，即点声源s和节点h直线距离r(t)的水平投影，ρ(t)随动点s'的位置改变而变化。入射角为θi(t)的声源发出的声线在点t处发生折射进入水下，到达水听器节点h。
[0058]
点声源s的速度v可以在飞行水平面上分解为两个相互垂直的分量，一个分量v1(t)指向水听器方向，其大小为
[0059]
v1(t)＝vcosα(t)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
[0060]
其中，α(t)为偏向角，是轨迹在当前位置处切线与直线s'h的夹角。进一步地，v1(t)可以在垂面上分解为两个相互垂直的分量，一个分量v
11
(t)与入射声射线方向一致，其大小为
[0061]v11
(t)＝v1(t)cosβ(t)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2)
[0062]
其中，空气中视线角β(t)与入射角θi(t)互为余角。
[0063]
由于点声源s和水听器节点h相对速度随时间发生变化，因此水听器接收信号的瞬时频率也随时间发生变化。显然，水听器所接收到的多普勒频移曲线fd(t)由参数{f0,v,β(t),α(t)} 决定。对于固有频率为f0的动点声源s，水听器在t时刻所接收到的多普勒频率为
[0064][0065]
公式(3)给出多普勒频移fd(t)与声源运动参数之间的关系。由于难以直接求解偏向角α(t)与空气视线角β(t)随时间t的变化关系，因此需要将该角度函数转化为声源s与水听器节点h之间位置关系，以方便求解多普勒频移。
[0066]
步骤2：根据声源与水听器相对运动的几何规律，推导出多普勒频移与声源频率、速度、高度之间的函数关系。
[0067]
首先，根据图2所示的俯视图中动声源s与水听器节点h的位置关系，可以求解偏向角α(t)。在时间t时刻，声源s与节点h的坐标分别为(xs(t),ys(t))、(xh,yh)，因此偏向角可由声源s与水听器节点h的坐标位置求解
[0068][0069]
然后，根据图3所示的声线传播路径几何关系，可以求解空气中视线角β(t)。动点声源s和静止节点h之间的侧视图，展示了入射角为θi(t)和折射角为θ
t
(t)的声线传播路径关系。
[0070]
在以∠sth为顶角δsth中，由余弦定理可知入射角θi(t)和折射角θ
t
(t)满足如下关系，
[0071][0072]
其中，声波在空气介质与水介质的传播距离分别为ra(t)和rw(t)，可写为
[0073][0074]
点声源s和水听器节点h的直线距离r(t)为
[0075]
r2(t)＝(d h)2 v2t2 ρ2(t)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(7)
[0076]
入射角θi(t)和折射角θ
t
(t)之间满足的snell定理，且与空气视线角β(t)的互补关系和水下视线角的互补关系可表示为
[0077][0078]
将式(8)代入(5)后，整理可得水下视线角与点声源高度h、声源速度v、偏航距离 w
min
、水听器深度d之间的关系
[0079][0080]
最后，将式(4)、(8)和(9)代入式(3)中，可以消去偏向角α(t)、空气视线角β(t)、折射角θi(t)和折射角θ
t
(t)，化简后可以得到方程组
[0081][0082]
对于方程组(10)，水下视线角作为可消参数，因此可以确定多普勒频移fd(t)与声源运动参数{f0,v,d,h,ρ(t)}之间的函数关系。
[0083]
由于声源s在水平方向沿曲线航迹运动，静止水听器h观测声波，水听器h观测声源s 的偏向角α(t)与空气视线角β(t)会随声源s与水听器h相对位置改变而变化。在图2所示的坐标系中，设定声源s与水听器h之间距离变小的阶段，速度v为正；设定声源s与水听器h 之间距离变大的阶段，速度v为负。
[0084]
由式(4)可知，声源s在运动过程中，偏向角α(t)变化范围为[0,π]rad。其中，声源s朝向水听器h运动时，偏向角α(t)∈[0,π/2]rad，远离水听器h运动时，偏向角α(t)∈[π/2,
π]rad。特别地，若声源s相距水听器h无限远，偏向角为或速度v的瞬时方向垂直于声源s与水听器h连线时，偏向角α(t)为π/2rad。由式(9)可知，空气视线角β(t)变化范围为[0,arctanh/ρ
min
(t)]rad。其中，若声源s相距水听器h无限远，空气视线角为声源s经至水听器h最近点时，空气视线角β(t)为arctanh/ρ
min
(t)rad弧度。声源s朝向水听器h运动过程中，观测到的多普勒频率大于声源s的声波频率，即fd(t)＞f0(t)，声源s远离水听器h运动过程中，观测到的多普勒频率小于声源s的声波频率，即fd(t)＜f0(t)，速度v的瞬时方向垂直于声源s与水听器h连线时，观测到的多普勒频率等于声源s的声波频率，即fd(t)＝f0(t)。
[0085]
步骤三，应用声源和水听器的几何关系、多普勒频移曲线，估计出声源的运动参数。
[0086]
在建立圆周飞行轨迹模型和分析三维多普勒频移曲线基础上，若已知空气介质声速ca、水介质的声速cw、水听器深度d，本节使用单水听器所接收的水声数据来估计多普勒频移曲线fd(t)，并根据式(10)、声源与水听器几何关系可以估计直升机的固有频率f0、速度v、飞行高度h。
[0087]
图4展示了声波传播几何关系俯视图，声源s在时刻ta、tm和tc依次经过不共线的三点a、 m、c，其中点m为弧线的中点，即由于声源s在点a发出的声波传播至水听器h具有一定时间差δt
ab
，水听器h接收到声波信号时，声源s已运动至点b处。若时间差δt
ac
大于声波传至水听器时间，则点b位于弧上，如图4(a)所示；若时间差δt
ac
小于声波传至水听器时间，则点b位于弧外，如图4(b)所示。
[0088]
在图4所示的坐标系中，声源经过点{a,m,c}是等间隔的，令点{a,m,c}的坐标分别为
[0089][0090][0091][0092]
这里，
[0093][0094]
为求解声源飞行参数信息，本节需确定辅助点b的信息，其为声源由点a处航行δt
ab
时间后到达的位置，然后根据点{a,b,c,m}的几何关系来求解参数。
[0095]
首先，令点b坐标为
[0096][0097]
这里，
[0098][0099]
其中，sgn是符号函数。
[0100]
声源位于点a时，声线由声源s传至水听器节点h的传播时间为
[0101][0102]
根据弧长公式，声源s一定时间内经过的弧长与其弧度有如下关系
[0103]
v|t
b-tm|＝rσ
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(16)
[0104]
其中，声源到达点b的时刻为
[0105]
tb＝ta τaꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(17)
[0106]
然后，根据点{a,b,c,m}的位置关系，建立几何关系。以圆心点o为顶点、选取圆弧上的不同的二点连线(即)为底边所构成的三角形应用余弦定理，有
[0107][0108]
对于∠poq，有如下关系式
[0109][0110]
以节点投影h'为顶点、圆弧上的任二点为底边所构成的三角形应用余弦定理，有
[0111][0112]
对于∠ph'q，有如下关系式
[0113][0114]
对于ρ
p
，有如下关系式
[0115][0116]
以节点投影o为顶点、h'p为底边所构成的三角形应用余弦定理，有
[0117][0118]
声源位于p∈{a,b,c,m}时，以折射点t
p
为顶点、声源s与水听器h连线为底边所构成的三角形，其几何关系如图3所示，应用余弦定理可以得到
[0119][0120]
根据弧长公式，声源s一定时间内经过的弧长与其弧度有如下关系
[0121]
v|t
a-tm|＝rω
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(25)
[0122]
声源位于p∈{a,b,c,m}时，其瞬时多普勒频率为
[0123][0124]
联立式(11)—(26)共56个方程，并将(6个)、(6个)、(6个)、(4个)、x
p
(4 个)、y
p
(4个)、ρ
p
(4个)、β
p
(4个)、(4个)、y0、τa、r、ω、ω'、ω”、σ、σ'、υ、υ'等52个因子视为可消参数，可以求解出直升机固有频率f0、速度v、飞行高度h。
[0125]
选定直升机运行轨迹中等间隔三点，根据式(11)、(12)设定三点坐标，并利用式(13)—(17) 设定辅助点b的位置。
[0126]
然后，由式(18)所示的余弦定理可得，
[0127][0128]
由式(20)所示的余弦定理可得，
[0129][0130]
其中，
[0131][0132]
由式(23)所示的余弦定理可得，
[0133][0134]
由式(24)所示的余弦定理可得，
[0135]
[0136][0137]
在点{a,b,c,m}处，多普勒频率与固有频率、速度、水下视线角有如下关系
[0138][0139]
matlab中fsolve函数采用最小二乘法，可以实现非线性方程组的求解。在图4中，选取3个等间隔时刻及其多普勒频移信息，即(ta,f(ta))、(tm,f(tm))和(tc,f(tc))。联立以上方程组，可求解直升机的飞行参数。
[0140]
本领域普通技术人员可以理解，以上所述仅为发明的优选实例而已，并不用于限制发明，尽管参照前述实例对发明进行了详细的说明，对于本领域的技术人员来说，其依然可以对前述各实例记载的技术方案进行修改，或者对其中部分技术特征进行等同替换。凡在发明的精神和原则之内，所做的修改、等同替换等均应包含在发明的保护范围之内。