基于散射声场分离算法在浅海信道中测量方法

1.本发明涉及声全息领域中浅海信道中目标识别的测量方法，特别涉及基于散射声场分离算法在浅海信道中测量方法。


背景技术：

2.声全息是将全息照相理论引用到声学领域而形成的声成像技术，根据成像距离的不同分为：近场声全息、远场声全息、常规声全息。与传统的全息技术相比，近场声全息由于测量面位置的选取，使其既包含了远场传播波成分，又包含了近场倏逝波成分，从而突破瑞利距离的限制,可以精确的识别噪声源和可视化空间声场。由于这样的优点和特点，使得近场声全息在声源的分析方面具有广阔的应用前景，成为噪声源识别和定位的重要方法之一。近场声全息技术的核心是其全息变换算法。其基本原理是利用基阵测量得到的目标辐射声场的声压或振速等声学量，利用不同的重构算法，进行声场的可视化重建或预测。
3.常规的近场声全息对环境要求比较严格，要求声源都在测量面的同一侧，另一侧为自由场环境。而在实际工程的测量环境下，这一点是很难保障的。此时进行重建的结果会因为存在其他目标的反射、散射等干扰，重建精度难以保障，甚至会出现完全影响判断。因此，对于非自由场的nah技术研究日益成为关键问题。
4.对于目标声散射问题，最原始的是积分方程法，其基本原理是helmholtz积分公式，理论上可计算任意复杂形状目标的散射声场，但积分方程本身的求解比较困难，但是其有明显的缺点，解的不稳定性和巨大的计算量。rayleigh提出分离变量，结合边界条件求解简正级数，该方法可对无限长圆柱、球等给出严格解析解，仅适用于规则形状目标，此外，该方法的解为无穷级数形式，导致解的收敛性比较差。已经提出的严格理论解都只适用于一些简单目标。于是，对于复杂目标的散射声场计算，出现了多种数值解法和近似解法。比如：t矩阵法、物理声学计算方法、亮点模型和板块元理论等方法。有限元法是一种十分常见的数值计算方法，通过对区域离散化近似，线性叠加来完成计算，但在单独处理散射声场问题时，尤其是复杂、大尺寸目标和高频计算时，计算速度慢，对计算机硬件的要求高。相较于有限元方法，边界元法计算精度高，计算量减小，但面临大阶次矩阵的求逆计算，因而对具有复杂形状的目标建模十分困难。


技术实现要素：

5.本发明的目的在于提供基于散射声场分离算法在浅海信道中测量方法，通过声场分离算法，有效的分离得到入射声和散射声，重构散射声场。为解决目标散射声场特性分析和测量问题提供新的研究途径，以解决上述背景技术中提出的问题。
6.为实现上述目的，本发明提供如下技术方案：
7.基于散射声场分离算法在浅海信道中测量方法，包括以下步骤：
8.步骤一：建立基于弹性结构浅海信道散射声场模型，采用双层柱面阵列进行声压数据采样，提取浅海信道中两同轴柱形全息面上的复声压数据；
9.步骤二：提取两柱形全息面数据，计算全息面声压角谱；
10.步骤三：采用柱面声场分离技术对两柱形测量全息面声压数据进行分解，分离得到散射声场；
11.步骤四：重构散射声场，计算声场重构误差。对重建面的声压角谱进行波数域加窗。对加窗后的重建声压角谱进行傅里叶逆变换获得重建面声压。
12.进一步地，步骤一结合待分析的模型及有限空间三维尺寸，建立弹性结构浅海信道声场模型，并进行网格划分，水域网格按照一个波长内不少于六个点的规则，采用自由四面体网格进行网格划分。
13.进一步地，步骤二对于理想流体介质定义为连续并且在运动的过程中不存在能耗问题的介质，三维环境下，声波在理想流体介质中的传播规律方程，表示为三个基本方程：
[0014][0015][0016][0017]
式(1)、(2)中，ρ0在声学中表示传播介质密度；式(3)中ρ
′
是密度逾量，表示介质中有声场时密度和无声场时的密度之差，是一个与时间和空间相关的物理变量。式(3)中c0表示介质内的声音传播速度；式(1)、(2)中v和p分别表示声场中的质点振速和声压。
[0018]
进一步地，步骤三声源激发的散射场和边界反射激发的散射声场视为目标的散射声场，声源具有声源1和声源2，两声源分布在两个全息面两侧。
[0019]
进一步地，步骤四对连续的声压信号进行加窗处理，降低在柱面近场声全息的重构过程中由于对连续的声压信号进行有限截断而产生的频谱泄漏误差。
[0020]
进一步地，步骤四的重建的计算过程中对波数域进行加窗，滤除高波数误差，其中，采用的滤波器为低通圆对称指数滤波器：
[0021][0022]
其中，α为窗函数的陡峭系数；k
c
为截止波束。
[0023]
与现有技术相比，本发明的有益效果是：本发明提出的基于散射声场分离算法在浅海信道中测量方法，通过有限元软件对声场进行仿真，获取更为精准的两柱形全息面复声压数据，通过声场分离算法，有效的分离得到入射声和散射声，重构散射声场。为解决目标散射声场特性分析和测量问题提供新的研究途径。
附图说明
[0024]
图1为本发明的方法流程图；
[0025]
图2为本发明步骤一中所建立的点声源
‑
弹性球浅海信道散射声场模型；
[0026]
图3为本发明步骤二中所声源和两测量面之间的关系；
[0027]
图4(a)为步骤二中声压数据提取主视图；
[0028]
图4(b)为步骤二中声压数据提取x
‑
z面视图；
[0029]
图5为本发明中选取的声压值对比轴位置；
[0030]
图6(a)为全息面距离目标z
h
＝0.17m时，目标在全息面1声压幅值对比；
[0031]
图6(b)为全息面距离目标z
h
＝0.2m时,目标在全息面2上声压幅值对比；
[0032]
图6(c)为全息面距离目标z
h
＝0.17m时,声源在全息面1声压幅值对比；
[0033]
图6(d)为全息面距离目标z
h
＝0.2m时，声源在全息面2上声压幅值对比；
[0034]
图7(a)为全息面距离目标z
h
＝0.17m时,目标在全息面1声压幅值对比；
[0035]
图7(b)为全息面距离目标z
h
＝0.2m时,目标在全息面2上声压幅值对比；
[0036]
图7(c)为全息面距离目标z
h
＝0.17m时,声源在全息面1声压幅值对比；
[0037]
图7(d)为全息面距离目标z
h
＝0.2m时,声源在全息面2上声压幅值对比。
具体实施方式
[0038]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0039]
基于散射声场分离算法在浅海信道中测量方法，如图1所示，包括以下步骤：
[0040]
步骤一：建立基于点声源
‑
弹性球浅海信道散射声场模型，提取浅海信道中全息面上的复声压数据。
[0041]
结合待分析的模型及有限空间三维尺寸，建立点声源
‑
弹性球浅海信道声场模型，建立网格模型时，可采用有限元仿真软件comsol建立模型，如图2所示，并进行网格划分，水域网格按照一个波长内不少于六个点的规则，采用自由四面体网格进行网格划分，图2中a表示海面无限软边界，b表示海底无限硬边界，c表示边界反射，d表示入射声，e表示散射声，f表示边界反射的散射声，g表示声源，m表示测量矩阵。
[0042]
步骤二：提取两柱形全息面数据，计算全息面声压角谱。
[0043]
对于理想流体介质一般可以定义为连续并且在运动的过程中不存在能耗问题的介质。三维环境下，声波在理想流体介质中的传播规律方程，可以表示为三个最近本方程，它们的表达式如下：
[0044][0045][0046][0047]
式(1)、(2)中，ρ0在声学中表示传播介质密度；式(3)中ρ
′
是密度逾量，表示介质中有声场时密度和无声场时的密度之差，是一个与时间和空间相关的物理变量。式(3)中c0表示介质内的声音传播速度；式(1)、(2)中v和p分别表示声场中的质点振速和声压。
[0048]
联立上述(1)、(2)、(3)三个方程可以得到均匀、静止理想流体中小振幅波波动方程：
[0049][0050]
上式(4)中的
▽2为laplace算子，基于直角坐标系中的函数关系表达式如下：
[0051][0052]
若令x＝r cosθ，y＝r sinθ，则可以将helmholtz方程从直角坐标系转换到柱面坐标系对应于柱面坐标系的laplace算子的表达式为：
[0053][0054]
利用分离变数法对柱面坐标系下的helmholtz方程求解，首先将偏微分方程改写为几个常微分方程乘积的形式，假设该方程解的形式可以写成：
[0055]
p(r,θ,z)＝r(r)θ(θ)z(z)
ꢀꢀꢀ
(7)
[0056]
可以得到柱面坐标系下的helmholtz方程式：
[0057][0058]
式(8)中，第一个括号内的项只与变量r和θ有关，第二个括号内的项只与变量z有关，并且，自变量r只出现在r中，自变量θ只出现在θ中，自变量z只出现在z中，也就是说，自变量r、θ、z是相互独立的三个变量。若等式想要成立，则必须满足两个括号内的项均为常数。于是，把这两个常数分别记作k
z
和k
r
设
[0059][0060][0061]
其中，常数k
r
满足
[0062][0063]
同时，可以改写为
[0064][0065]
若上式(12)成立，也需满足方程左右两边的项均为常数的条件，因此，选择常数n2，令其满足如下方程
[0066][0067]
将上述表达式(13)代入方程，并在方程两边同时除以r2/r，可以得到
[0068][0069]
该式为贝塞尔方程。
[0070]
由数学物理方程知识可知，贝塞尔方程的解包含第一类贝塞尔函数j
n
(k
r
r)和第二
类贝塞尔函数y
n
(k
r
r)，其中，第二类贝塞尔函数也被称作纽曼函数，可以用符号n
n
进行表示。故解可以由第一类贝塞尔方程和第二类贝塞尔方程的线性组合
[0071]
r(r)＝c1j
n
(k
r
r)+c2y
n
(k
r
r)
ꢀꢀꢀ
(15)
[0072]
进行表示，该式(15)表示贝塞尔方程的驻波解，其中，c1、c2为任意常数。
[0073]
贝塞尔方程的行波通解表达式需要由第一类汉克尔函数和第二类汉克尔函数进行线性表示
[0074][0075]
其中，第一类汉克尔函数对应于向外发散的波，第二类汉克尔函数对应于向内收敛的波。
[0076]
在中，选取k
z
时并无任何条件限制，但是，在上述讨论中，默认即有
[0077]
而当k
z
＞k时，k
r
应写作是一个虚数，此时，应改写为修正的贝塞尔方程。
[0078][0079]
该方程的解可以由第一类修正贝塞尔函数i
n
(k
′
r
r)和第二类修正贝塞尔函数k
n
(k
′
r
r)线性表示，其中应该注意的是，第二类修正的贝塞尔函数k
n
(k
′
r
r)对应于第一类汉克尔函数的修正结果，第一类修正的贝塞尔函数对应于第二类汉克尔函数的修正结果。
[0080]
总结以上分析推导，可以得到柱面坐标系下的helmholtz方程的通解表达式为
[0081][0082]
式(18)中，d1(k
z
z)和d2(k
z
z)可以是任意常数。
[0083]
将p
n
(r,k
z
)定义为声压量在位置r处关于周向变量θ和轴向变量z的空间fourier变换，即为柱面波数谱：
[0084][0085]
步骤三：采用柱面声场分离技术对测量全息面声压数据进行分解，分离得到散射声场。
[0086]
双柱面声场分离，顾名思义，应用该技术时需要两个全息测量柱面，全息面与声源的空间分布情况如下图3所示。
[0087]
图3中，s
h1
和s
h2
分别是半径为r
h1
和r
h2
的全息柱面，在分析空间内存在两个声源，分别为a和b，其中a表示声源2，b表示声源1，两声源分布在两个全息面两侧，其中，声源1分布在全息柱面s
h1
内部，声源2位于全息柱面s
h2
之外，根据声场叠加原理，可知全息面上总声场应为两声源辐射声场之和，可以记作
[0088]
p(r
h1
,θ,z)＝p
i
(r
h1
,θ,z)+p
ii
(r
h1
,θ,z)
ꢀꢀꢀ
(20)
[0089]
p(r
h2
,θ,z)＝p
i
(r
h2
,θ,z)+p
ii
(r
h2
,θ,z)
ꢀꢀꢀ
(21)
[0090]
式(20)和(21)中，p
i
(r
h1
,θ,z)和p
ii
(r
h1
,θ,z)分别代表了声源1和声源2在全息柱面s
h1
上产生的声压值，同理，p
i
(r
h2
,θ,z)和p
ii
(r
h2
,θ,z)分别代表了声源1和声源2在全息柱面s
h2
上产生的声压值。
[0091]
对空间域声压进行空间fourier变换可以得到波数域声压，即声波的柱面波谱，做空间傅里叶变换可得：
[0092][0093][0094]
根据近场声全息理论，不同的柱面之间的声压柱面波谱具有一定的转换关系，因此，可以得到如下的表达式，该表达式描述了同一声源在全息面s
h1
和全息面s
h2
上产生的声压柱面波谱的关系：
[0095][0096][0097]
联立上述(24)和(25)方程式，便可解得声源1在全息面s
h1
和s
h2
产生的声压波谱以及声源2在全息面s
h1
和s
h2
产生的声压波谱，具体公式如下：
[0098][0099][0100][0101][0102]
步骤四：重构散射声场，计算声场重构误差。对重建面的声压角谱进行波数域加窗。对加窗后的重建声压角谱进行傅里叶逆变换获得重建面声压。
[0103]
因为全息测量面的实际操作过程中，全息孔径上往往只能选取全息面上一定面积的离散点处的声压信号，也就是相当于是对连续的声压信号进行了一个有限截断，所以在近场声全息的重构过程中就会产生频谱泄漏误差，这种情况通常会被称之为“有限孔径效应”或者“窗效应”。“有限孔径效应”是全息面上的声压角谱p(k
x
,k
y
,z
h
)产生计算误差的主要原因。
[0104]
离散fourier变换存在波数域采样的过程，在这个波数域采样的过程中会产生“卷绕误差”，通常又被称之为孔径重复效应。以上误差的卷入，都会导致重建结果的失真。
[0105]
重建过程对于空间频率的误差非常敏感，因此必须采用适当的措施来抑制高空间频率成分误差的干扰，一种行之有效的方法就是在计算过程中对波数域进行加窗，滤除高波数误差。其中，最常被采用的滤波器要属低通圆对称指数滤波器：
[0106][0107]
其中，α为窗函数的陡峭系数；k
c
为截止波束。
[0108]
对本发明的有益效果如下方式得以验证：
[0109]
在comsol数值计算和matlab仿真条件下，对该方法进行仿真实验：
[0110]
仿真参数如下：
[0111]
信道塑料球壳验证仿真：
[0112]
目标球壳采取pvc材料，密度1380kg/m3，泊松比0.31，杨氏模量3.5
×
109pa
[0113]
球半径r＝0.15m，厚度0.01m球壳中心位于原点处。
[0114]
点声源位置：(0.5，0，0)。
[0115]
信道环境设置：12.8m*12.8m*5m。
[0116]
流体介质参数：水，密度1000kg/m3，声速1500m/s。
[0117]
模型边界条件：模型上边界为水面和空气交界面，为绝对软边界，反射系数为
‑
1，下边界为绝对硬硬边界，反射系数为1。
[0118]
计算频率：f＝200hz
‑
2khz。
[0119]
全息面测量参数：矩柱形全息面，全息面轴向高度l＝4m,周向角度间隔10
°
，轴向采样间隔dz＝0.04m,全息面距离目标球中心距离0.17m、0.20m、0.25m、0.30m、0.35m、0.40m、0.45m。
[0120]
重建面距离目标球中心距离：0.17m、0.20m。
[0121]
声压数据提取如图4，图4(a)为声压数据提取主视图，图4(b)为声压数据提取x
‑
z面视图，图5为选取的声压值对比轴位置，轴向声压仿真结果如图6所示，其中图6(a)为全息面距离目标z
h
＝0.17m时，目标在全息面1声压幅值对比。图6(b)为全息面距离目标z
h
＝0.2m时,目标在全息面2上声压幅值对比。图6(c)为全息面距离目标z
h
＝0.17m时,声源在全息面1声压幅值对比。图6(d)为全息面距离目标z
h
＝0.2m时，声源在全息面2上声压幅值对比。
[0122]
其中图7为三维效果对比图，图7(a)为全息面距离目标z
h
＝0.17m时,目标在全息面1声压幅值对比。图7(b)为全息面距离目标z
h
＝0.2m时,目标在全息面2上声压幅值对比。图7(c)为全息面距离目标z
h
＝0.17m时,声源在全息面1声压幅值对比。图7(d)为全息面距离目标z
h
＝0.2m时,声源在全息面2上声压幅值对比。
[0123]
以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明披露的技术范围内，根据本发明的技术方案及其发明构思加以等同替换或改变，都应涵盖在本发明的保护范围之内。