本发明涉及参量阵声场的计算，特别涉及基于准蒙特卡洛法的参量阵声场高精度快速计算方法。
背景技术：
1、当参量阵(pa,parametric array)在介质中辐射不同频率的两个超声波束时，由于二阶非线性，可以产生高度定向的差频声，其应用已相当广泛，例如空气中主动噪声控制、局部区域发声、水下通信和浅剖探测等，其中准确计算pa产生的差频声场非常重要。与传统声源相比，pa产生的差频声场更加复杂，在数学形式上涉及五重积分，计算成本相当巨大。而在数学和工程科学计算中，通常采用蒙特卡罗方法、准蒙特卡罗方法计算多重积分。但蒙特卡罗方法的误差是概率性的，得到解的精度不高；准蒙特卡罗方法是用确定性的超均匀分布序列代替蒙特卡罗方法中的随机数序列，可以得到确定性的误差估计和更高精度的解。
2、aanonsen指出，在不太靠近辐射表面的距离上，局部效应可以忽略不计，这使得二阶波动方程可以简化为westervelt方程(1aanonsen s i,barkve t,tjo/tta j n,etal.distortion and harmonic generation in the nearfield of a finite amplitudesound beam[j].the journal oftheacoustical society ofamerica,1984,75(3):749-768.)。由于多数情况下pa辐射的原波声压水平是有限大小的，可以使用准线性近似，这进一步简化了westervelt方程。使用pa表面上的二重瑞利积分(ri,rayleigh integral)来获得原波声场，然后构建无限大体积的虚源，源密度与原波压力的乘积成比例。最后使用整个空间上的三重体积积分获得差频声场。在弱非线性的情况下，除了在接近pa的位置之外，该解在整个空间中都是精确的，为了降低准线性解的计算成本，wen最先提出了高斯波束展开法(gbe,gaussian beam expansion)以简化瑞利积分和虚源体积积分的计算(2wen j j,breazeale ma.adiffraction beam field expressed as the superpositionofgaussian beams[j].the journal ofthe acoustical society ofamerica,1988,83(5):1752-1756.)，然而，gbe方法存在两个问题：一是当波束转向角超出近轴近似极限(约20°)时，轴上预测将不准确，二是在大转向角和低差频波频率下计算体积积分时，采用近轴近似会导致大的误差。后来，为了克服上述问题，zhuang等人提出了采用非近轴近似的高斯波束展开方法(ngbe,non-paraxial gaussian beam expansion)来计算可转向pa产生的声场(3zhuang t,zhong j,niu f,et al.a steerable non-paraxial gaussian beamexpansion for a steerable parametric array loudspeaker[j].the journaloftheacoustical society ofamerica,2023,153(1):124-136.)，然而，由于计算原理的限制，ngbe方法在近场或远离中心轴的地方仍可以观察到误差，另外，对于聚焦型换能器或波导环境等更复杂情况下，这些精度问题可能会更严重。zhang提出使用准蒙特卡罗方法(qmc,quasi-monte carlo)对超声换能器产生的声场进行建模(4zhang s,huangy,li x,etal.modeling of wave fields generated by ultrasonic transducers using a quasi-monte carlo method[j].the journal ofthe acoustical society ofamerica,2021,149(1):7-15.)，通过判断采样法和halton序列对声源进行伪随机采样，然后利用样本均值求解积分得到了声场分布，在一定程度上简化了声场计算，但其只考虑了单一频率下声源的辐射声场解，未考虑声场的非线性效应，仍属于线性声学范畴。
3、基于前人的理论，本发明提出了使用准蒙特卡洛方法计算pa非线性声场的方法。在准线性近似下，利用准蒙特卡洛方法对pa声源进行建模，求出原波声场解，然后考虑二阶非线效应，获得pa的差频声场解，有利于更快速地求解pa非线性声场的准线性解。其优点是在不采用额外近似，保证精度的前提下，提高了计算效率。
技术实现思路
1、针对现有技术中存在的技术问题，本发明提供了基于准蒙特卡洛法的参量阵声场高精度快速计算方法。所述方法通过对westervelt方程采用准线性近似和逐次逼近，参量阵差频声场的准线性解来源于无限大虚源的辐射，使用准蒙特卡洛方法求解原波声场，可以在不使用额外近似的条件下，准确快速地得出虚源密度的空间分布，最终获得参量阵声场分布。
2、本发明是通过以下技术方案实现的，本发明提出基于准蒙特卡洛法的参量阵声场高精度快速计算方法，所述方法包括以下步骤：
3、步骤1、对自由场均匀介质中的westervelt方程使用准线性近似和逐次逼近，从而分别给出稳态双频激励下原波声场和差频声场各自满足的helmholtz方程；
4、步骤2、对于平面活塞声源，原波声场满足的helmholtz方程的解利用green函数构造为瑞利表面积分ii；
5、步骤3、基于qmc方法，在活塞声源面s内，选取概率密度函数抽取服从的n个样本点并引入修正的green函数将瑞利表面积分ii近似成随机变量的数学期望的形式；
6、步骤4、对于振速均匀的平面活塞，取利用伪随机样本点来计算瑞利积分的近似值qi，从而得到原波声场解
7、步骤5、通过已求出的原波声场解得到虚源密度函数
8、步骤6、当活塞声源所处平面为刚性时，差频声场满足的helmholtz方程的解可利用green函数构造成关于虚源的体积积分形式，对该积分进行数值求解得到差频声场解
9、进一步地，在步骤1中，从自由场均匀介质中的westervelt方程出发
10、
11、其中p是声压，t是时间变量，c0是小信号下的绝热声速，等号右边第一项表征流体热粘度，δ是耗散系数，与频率ω处的声衰减系数α有关，关系式为等号右边第二项表征非线性，ρ0是环境流体密度，β是非线性参数；
12、半径为a的圆形pa产生频率为fi的原波，换能器表面的边界条件为
13、
14、其中j为虚数单位，vz表示垂直于换能器表面的振动速度，为声源表面某点到换能器中心的矢径，为处振速的幅度，ωi＝2πfi(i＝1,2)是原波的角频率；
15、基于准线性近似和逐次逼近方法，由式(1)可得原波声场满足以下齐次helmholtz方程
16、
17、其中ki＝ωi/c0+jαi和αi分别是频率为fi的原波的波数和声衰减系数；为场点位置，以pa中心为原点建立直角坐标系o-xyz，z轴垂直于换能器表面；
18、差频声场满足以下非齐次helmholtz方程
19、
20、其中kd＝ωd/c0+jαd和αd分别是频率为fd＝f1-f2的差频波的波数和声衰减系数，等号右边为虚源项且ωd＝ω1-ω2是差频角频率，为空间某虚源点位置。
21、进一步地，在步骤2中，对于平面活塞，式(3)的解利用green函数构造为瑞利表面积分ii
22、
23、其中是声源表面s上某点到空间场点的距离。
24、进一步地，在步骤3中，采用qmc方法，引入判断采样技术和halton序列，从目标声源区域获得有效数量的样本用于积分：首先在声源面s内选取概率密度函数且当时有引入修正的green函数
25、
26、从而声场可以被近似为随机变量的数学期望的形式
27、
28、进一步地，在步骤4中，选择足够数量的n个遵循概率密度函数的伪随机样本，然后将每个位置处的压力场近似表示成所有n个样本结果的平均值qi。
29、进一步地，在步骤5中，根据上一步解，可以得到空间处的虚源密度函数为
30、
31、其中是空间处的原波声压，*号代表复共轭。
32、进一步地，在步骤6中，在准线性近似下，差频声场可以被视为无限多个虚源在处产生的声音的叠加，因此式(4)的解可利用green函数构造成关于虚源的体积积分来求解
33、
34、其中是某虚源点到空间场点的距离，通过数值求解式(9)即可获得差频声场解。
35、本发明提出了使用准蒙特卡洛方法计算参量阵声场的方法。在不使用额外近似，只采用准线性近似的条件下，利用准蒙特卡洛方法对参量阵声源进行建模，从而快速求出原波声场解，有利于更高效地求解参量阵声场的准线性解。为今后更复杂的声源形式或者波导情况下参量阵声场的计算奠定基础，对水下参量阵的实际应用有积极的指导意义。
36、本发明的有益效果可以归纳为：
37、(1)利用准蒙特卡洛方法对参量阵声源进行建模，从而快速求出原波声场解，更高效地获得参量阵声场准线性解。
38、(2)可以计算出不同声源形式下的参量阵声场解，误差小，应用范围广；
39、(3)计算结果准确可靠，计算方法简单易行。
1.基于准蒙特卡洛法的参量阵声场高精度快速计算方法，其特征在于：所述方法包括以下步骤：
2.根据权利要求1所述的方法，其特征在于，在步骤1中，从自由场均匀介质中的westervelt方程出发
3.根据权利要求2所述的方法，其特征在于，在步骤2中，对于平面活塞，式(3)的解利用green函数构造为瑞利表面积分ii
4.根据权利要求3所述的方法，其特征在于，在步骤3中，采用qmc方法，引入判断采样技术和halton序列，从目标声源区域获得有效数量的样本用于积分：首先在声源面s内选取概率密度函数且当时有引入修正的green函数
5.根据权利要求4所述的方法，其特征在于，在步骤4中，选择足够数量的n个遵循概率密度函数的伪随机样本，然后将每个位置处的压力场近似表示成所有n个样本结果的平均值qi。
6.根据权利要求5所述的方法，其特征在于，在步骤5中，根据上一步解，可以得到空间处的虚源密度函数为
7.根据权利要求6所述的方法，其特征在于，在步骤6中，在准线性近似下，差频声场可以被视为无限多个虚源在处产生的声音的叠加，因此式(4)的解可利用green函数构造成关于虚源的体积积分来求解
