三列不同相位原波形成的参量阵声场研究

杨德森，李中政，2，方尔正

(1.哈尔滨工程大学水声工程学院，黑龙江哈尔滨150001;2.海军92198部队，辽宁兴城125109)



三列不同相位原波形成的参量阵声场研究

杨德森1，李中政1，2，方尔正1

(1.哈尔滨工程大学水声工程学院，黑龙江哈尔滨150001;2.海军92198部队，辽宁兴城125109)

摘要:参量阵能生成低频、高指向性声束，在海底地层剖面测量、掩埋目标探测与识别、水下通信等水声工程领域有着广泛的应用。理论研究表明，当三列原波相互作用时，能形成具有两个差频分量的参量阵辐射系统，但尚未有文献对受原波相位影响的多频参量阵声场进行系统分析。基于此，本文综合计算精度和效率的需求，借鉴算子分离思想，利用时频域相结合的方法求解KZK方程，实现对三列不同相位原波相互作用形成的参量阵声场的描述。计算结果表明:参量阵差频信号的声压幅值随着原波相位的改变而改变，并且原波相位均为零时，两个差频分量的声压幅值最大。最后进行了试验研究，实验结果与仿真结果吻合较好。

关键词:参量阵;三波作用;原波相位;算子分离法;时频域相结合

1962年，Westervelt［1］提出了参量阵的最初模型:同向传播的两列高频原波利用媒质的非线性效应获得差频波的声发射装置。1965年，Berktay［2］引入调制“包络”的概念，给出了声学参量阵更精确和完整的理论解释，使得声源辐射信号不再局限于双频基波的情况。Lucilla Di Marcoberardin［3］，John A.Birken［4］对多列声波辐射形成系列差频分量的参量阵声场进行了理论及实验研究，为海底地质特性、管道或电缆的铺设提供技术指导。2002年，王润田等［5］用有多个差频档位的“堤防隐患监测声呐”对堤防的损毁程度进行探测和评估。显然，多频原波相互作用形成的参量阵在水声工程领域有非常广泛的应用前景。

在众多非线性声场计算的理论模型中，KZK方程［6－8］在计算参量阵近场声场时有着非常明显的优势。Masahiko Akiyama等［9－11］利用KZK方程对不同相位双频原波形成的参量阵声场特性进行计算，并进行了试验研究。Fenlon等［12］基于Burgers方程推导了多频有限振幅声波的Bessel－Fubini解，吕君等［13］根据Fenlon理论研究了二阶近似条件下的双频声源相互作用时的远场指向性。为研究有限振幅波在介质中传播的频域能量变化问题，杨德森等［14］基于频谱分解方法推导了单频及多频大振幅的传播规律。

到目前为止，并没有相关文献对三列声波相互作用形成的参量阵声场特性进行研究。本文借鉴算子分离思想［8］，综合考虑计算效率和精度，将KZK方程分解为线性(包含衍射、吸收效应)和非线性(包含非线性效应)两部分，线性部分采用二阶龙哥库塔有限差分法(DIRK2)和Crank－Nicolson有限差分法(CNFD)相结合［7，18］的方法;非线性部分则采用守恒型迎风格式［15－17］积分求解。最后将两部分进行整合，对不同相位三列原波相互作用形成的参量阵声场进行分析，旨在为参量阵的进一步工程应用提供相应的理论指导。

1　KZK方程的相关理论

综合考虑有限振幅声波传播过程中的吸收、衍射和非线性效应，建立适用于参量阵声场计算的KZK方程［6－7］:

式中:gn、hn是第n次谐波声压在空间位置(rj，zm)处的谐波系数。

1.1二阶龙格库塔法计算KZK方程的线性部分

利用算子分离思想对式(1)的吸收、衍射部分进行频域计算，得到如下表达式:

用IBFD法进行有限差分离散，化为矩阵有

式中:AJ×J=三对角矩阵，

将式(3)进行改写:

所以，利用IBFD法对KZK方程的线性部分进行离散得到的局部截断误差为O(dz)。为了提高计算精度，与CNFD法的局部截断误差相匹配，根据二阶龙格库塔思想［17］，将既定的轴向步长dz分为b1dz、b2dz(b1+b2=1)两部分，分别计算zm，zm+b1Δz、处的斜率s1、s2:

1.2守恒迎风格式计算KZK方程的非线性部分

采用频域方法对KZK方程进行数值计算，非线性部分形象、直观地体现了各阶谐波的耦合，但计算时间正比于N2(N为最大截断谐波阶次，衍射、吸收项的计算时间正比于N)，故计算量非常大。参量阵的形成过程属于弱非线性声学范畴，其非线性系数较小，可以将频域信号变换为时域波形，换取较大计算步长达到提高计算效率的目的。

1.1节的求解得到zm+1处声压幅值p(rj，zm+1)是考虑声波非线性传播衍射、吸收效应，接下来将其作为节点(rj，zm)处的初始条件对KZK方程的非线性部分(无粘滞Burgers方程)求解［15－17］:

方程(7)两边同时对时间进行积分，将声压转换到时域，有如下表达式:

对无粘滞Burgers方程进行守恒型迎风格式表示:

由式(10)、(11)就可以计算空间层zm→zm+1的时域波形了，然后将其通过FFT变换到频域，以其作为初始条件代入KZK方程的线性部分，求解轴向位置zm+1→zm+2的声场传播。依次类推，即可获知整个求解区域的声场空间分布特性。

2　仿真结果分析

按照上述KZK方程求解参量阵声场的方案，设定其计算区域的网格划分情况如下，轴向计算范围为zmax=3.5r0，径向计算范围是rmax=41a，考虑到计算的方便，在频域计算过程中，对谐波阶数进行截断，选定N=128。

图1　参量阵声场空间的网格划分示意图Fig.1 Mesh grids of the parametric array’s acoustic field

另外，为了尽量降低边界的反射，本文借鉴求解电磁场的PML吸收边界条件，将径向网格的外围边界进行处理，设定1倍半径长度的边界层作为PML边界［15］(如图1阴影部分所示)。

选定如下仿真参数:水的密度ρ=1 000 kg/m3、声速c=1 482 m/s，水的非线性系数为β=3.5，圆形活塞声源的表面峰值声压为P0=60 kPa，并且随着相位的改变而变化。半径为a=10 cm，水中声波衰减系数与频率的关系为:α/f2=2.17×10－13dB/m;声源辐射的声波:

2.1同相位原波形成的参量阵声场空间分布

本节设定三列原波相位满足θf1=θf2=θf3=0，对参量阵各个差频分量的声场特性进行研究。

由图2可知:

1)差频分量的声能量累积过程大致相同，都存在一个拐点。原波在近场相互作用时，非线性效应比吸收效应明显，差频波的声能量累积速率较快;到达拐点之后，非线性效应的声能量累积速率小于吸收效应引起的声衰减速率，声吸收效应在参量阵的轴向声场传播过程中占据主导地位，参量阵差频波声压幅值越来越小;

2)三列原波相互作用形成两个具有高指向性的参量阵差频波束，并且两个差频分量的波束宽度基本一致，但差频分量fd2=8 kHz的轴向声压幅值比fd1=8 kHz的声压幅值大。

图2　参量阵的声压幅值分布特性Fig.2 The amplitude distribution of the parametric array’s pressure

2.2差频分量fd1的声压幅值分布特性

对不同相位原波相互作用形成的参量阵差频分量fd1=4 kHz的声场进行分析，给出如下组图3，由图3可知:

1)当θf2=θf3=0时，随着原波f1=41 kHz相位θf1的增加(如图3(a)的上图所示)，同一轴向位置处差频分量fd1=4 kHz的参量阵声压幅值越来越低，当θf1=θf3=0时，随着原波f2=45 kHz相位θf2的增加(如图3(a)的下图所示)，同一轴向位置处差频分量fd1=4 kHz的参量阵声压幅值越来越低，当θf1=θf3=0时，随着原波f3=49 kHz相位θf3的增加(如图3(b)的上图所示)，同一轴向位置处差频分量fd1=4 kHz的参量阵声压幅值越来越低。这表明三列原波相互作用时，原波叠加形成的调制信号包络随着相位的改变而改变，进而影响到差频信号fd1的轴向声压幅值分布;

2)相较原波f1=41 kHz和f3=49 kHz相位改变引发的差频分量fd1=41 kHz轴向声压幅值变化，原波f2=45 kHz相位改变对参量阵差频分量fd1=41 kHz的影响更为明显，这是因为差频分量fd1的形成过程为fd1=f3－f2，fd1=f2－f1，显然，其非线性作用机制较为复杂。

图3　原波相位改变时，参量阵声压幅值变化规律Fig.3 The parametric array’s pressure amplitude characteristics due to the primary wave’s phase changed

为了对相同相位的原波相互作用形成的参量阵声场进行分析，给出仿真组图4。

图4　参量阵差频分量fd1的声压幅值特性Fig.4　The pressure amplitude characteristic of the difference frequency component fd1

由图4可知:

1)三列原波的相位满足θf1=θf2=θf3=0时，轴向声压幅值最大，并且原波f2=45 kHz的相位θf2改变时，同一轴向位置处参量阵差频分量fd1的声压幅值最低;

2)原波f1=41 kHz和f3=49 kHz具有相同的相位时，形成的参量阵差频分量fd1声压幅值相等。说明非载波原波的相位变化对参量阵声场空间分布具有相同的影响。

2.3差频分量fd2的声压幅值分布特性

本节对不同相位原波相互作用形成的参量阵差频分量fd2的声压幅值分布特性进行研究，给出如下仿真结果。

图5　原波相位改变形成的参量阵声压幅值变化特性Fig.5 The parametric array’s pressure amplitude characteristics due to the phase changed of the primary wave

相对fd1=41 kHz来讲，参量阵差频分量fd2=8 kHz的声压幅值分布受原波相位的影响较小，这是因为差频分量fd2的形成过程是fd2=f3－f1，其非线性作用机制更为简单。

另外，fd2的声压幅值有小范围的波动，并且同一轴向位置的声压幅值随着相位的增加而降低，这是由于三列原波叠加形成的调制信号峰值降低导致参量阵声源辐射系统的声源级降低。如果忽略声源级带来的影响，那么就可认为两列不同相位原波形成的参量阵差频波声压幅值分布受相位的影响较小，跟文献［9－12］的研究结果相吻合。

同样，为了对相同相位的不同原波相互作用形成的参量阵声场进行分析，给出仿真组图6。

图6　差频分量fd2的声压幅值特性Fig.6　The pressure amplitude characteristic of the difference frequency component fd2

图6表明，原波相位相同时，同一径向位置的差频分量幅值最大。随着原波相位的增加，参量阵差频分量fd2的声压幅值越来越低，并且原波f1，f3的相位相同时，声压幅值的分布相同，亦即非中心频率分量的原波相位改变具有相同的功效。

3　试验研究

为了对不同相位原波形成的参量阵声场进行研究，本文在哈尔滨工程大学水声工程学院的消声水池对轴向声压幅值进行测量的实验。

图7　参量阵轴向声压幅值仿真值与实验值的比较=60°Fig.7 Comparison of the experimental value of the parametric array with the simulation values when=60°

选定原波f1、f2、f3的相位分别为=90°、==0和=60°、==0°两组工况，对仿真结果和实验结果进行比对。

对图7进行分析可知，三列不同相位的原波相互作用可以生成具有两个高指向性、差频分量的参量阵辐射系统。另外，由于KZK方程在数值计算参量阵声场空间分布时，仅考虑传播媒介的粘滞性吸收，导致远场声场由于吸收效应引起的声衰减小于实验测量值。

图8给出参量阵差频分量fd1的声压幅值随轴向距离的变化图示，其中图8(a)表示f1=41 kHz分别为30°、90°，θf1=θf3=0时的声压幅值轴向分布图示，图8(b)表示f1=45 kHz，θf2分别为0°、60°，θf1=θf3=0时的声压幅值轴向分布图示。

图8　参量阵轴向声压幅值仿真值与实验值比较Fig.8　Comparison of the experimental value of the parametric array with the simulation values of axial sound pressure amplitude

图9　参量阵轴向声压幅值仿真值与实验值比较Fig.9 Comparison of the parametric array’s experimental value with the simulation values

图9表明在误差允许范围内，参量阵差频分量fd1随着轴向距离的增加逐渐降低，跟仿真结果吻合程度较高。

接下来选定原波频率为f3=49 kHz、其相位θf3=0°、30°、60°，θf1=θf2=0时的参量差频分量fd2的轴向声压幅值分布。

对实验结果分析可知，差频分量fd2随轴向距离的变化规律跟KZK方程仿真计算结果基本一致，并且在近场区域其声压幅值变化规律更为接近，亦即KZK方程描述参量阵的近场声场优势较为明显。

4　结论

本文借鉴KZK方程的时域算子分离思想，将时域有限差分与频域有限差分相结合，把求解参量阵非线性声场的KZK方程分为两部分:采用守恒迎风格式对KZK方程的非线性部分进行计算，其中由于得到的线性部分采用二阶龙格库塔法对KZK方程的衍射、吸收效应进行计算，最后进行整合，得到三列不同相位原波相互作用形成的参量阵声场分布特性，通过研究得出如下结论:

1)DIRK2是具有二阶计算精度的有限差分法，该算法可提高计算精度，并且能跟远场计算的Crank－Nicolson有限差分法计算精度相匹配;

2)相对频域有限差分计算KZK方程，采用守恒迎风格式变换到时域求解参量阵声场的非线性部分可有效提高计算效率，降低计算量;

3)当三列原波相互作用时，能获得两个差频分量的高指向性声束，并且三列原波的相位均为0时，其轴向声压幅值取最大值;

4)三列不同相位的原波相互作用时，频率较低的差频分量fd1由于其非线性作用机制较为复杂，其声压幅值随相位改变有明显的振荡。与此同时，频率较高的差频分量fd2的声压幅值随相位的改变仅有小范围的波动。

参考文献:

［1］WESTERVELT P J.Parametric acoustic array［J］.The Journal of the Acoustical Society of America，1963，35(4):535－537.

［2］BERKTAY H O.Possible exploitation of non－linear acoustics in underwater transmitting applications［J］.Journal of Sound and Vibration，1965，2(4):435－461.

［3］DI MARCOBERARDINO L，MARCHAL J，CERVENKA P.Nonlinear multi－frequency transmitter for sea－floor characterization［C］//Proceedings of the Acoustics Conference.Nantes，France，2012:2800－2805.

［4］BIRKEN J A.Empirical results from frequency－scanning nonlinear sonar in deep water［J］.The Journal of the A－coustical Society of America，1974，56(S1):S41.

［5］王润田.海底声学探测与底质识别技术的新进展［J］.声学技术，2002，21(1/2):96－98.WANG Runtian.Progress in detecting the geological formations and sediment properties by sound［J］.Technical A－coustics，2002，21(1/2):96－98.

［6］KAMAKURA T，HAMADA N，AOKI K，et al.Nonlinearly generated spectral components in the near－field of a directive sound source［J］.The Journal of the Acoustical Society of America，1989，85(6):2331－2337.

［7］NAZE TJØTTA J，TJØTTA S，VEFRING E H.Propagation and interaction of two collinear finite amplitude sound beams ［J］.The Journal of the Acoustical Society of America，1990，88(6):2859－2870.

［8］LEE Y S，HAMILTON M F.Time－domain modeling of pulsed finite－amplitude sound beams［J］.The Journal of the Acoustical Society of America，1995，97(2):906－917.

［9］KAMAKURA T，SAKAI S，NOMURA H，et al.Parametric audible sounds fields by phase－cancellation excitation of primary waves［J］.The Journal of the Acoustical Society of A－merica，2008，123(5):3694.

［10］KAMAKURA T，NOMURA H，SAKAI S.Spatial phaseinversion technique for parametric source with suppressed carrier［J］.The Journal of the Acoustical Society of America，2009，125(4):2717.

［11］KAMAKURA T，NOMURA H，AKIYAMA M，et al.Parametric sound fields formed by phase－inversion excitation of primary waves［J］.Acta Acustica United with Acustica，2011，97(2):209－218.

［12］FENLON F H.An extension of the Bessel－Fubini series for a mutihle－frequency CW acoustic source of finite amplitude ［J］.The Journal of the Acoustical Society of America，1972，51(1B):284－289.

［13］吕君，赵正予，周晨，等.有限振幅声波间的非线性相互作用对声源远场指向性的影响［J］.物理学报，2011，60(8):084301.LV Jun，ZHAO Zhengyu，ZHOU Chen，et al.Effect of finite－amplitude acoustic wave nonlinear interaction on farfield directivity of sound source［J］.Acta Physica Sinica，2011，60(8):084301.

［14］杨德森，董雷，时洁，等.大振幅波非线性传播的频率特性［J］.哈尔滨工程大学学报，2010，31(7):928－931.YANG Desen，DONG Lei，SHI Jie，et al.Frequency characteristics of nonlinear propagation of large amplitude waves［J］.Journal of Harbin Engineering University，2010，31(7):928－931.

［15］SONESON J E.A parametric study of error in the parabolic approximation of focused axisymmetric ultrasound beams ［J］.The Journal of the Acoustical Society of America，2012，131(6):EL481－EL486.

［16］SONESON J E.A user－friendly software package for HIFU simulation［C］//Proceedings of the 8th International Symposium on Therapeutic Ultrasound.Minneapolis，Minnesota:AIP，2009，1113:165－169.

［17］PINTON G F.Numerical methods for nonlinear wave propagation in ultrasound［D］.Durham:Department of Biomedical Engineering，Duke University，2007:37－46.

［18］LI Zhongzheng，FANG Erzheng，SHI Shengguo.Study on the calculation method of the KZK equation for parametric array［C］//Processings of the 12th International Conference on Signal Processing(ICSP).Hangzhou，China:IEEE，2014:139－143.

Study of the parametric array acoustic field yielded by three collimated primary waves with different phases

YANG Desen1，LI Zhongzheng1，2，FANG Erzheng1
(1.College of Underwater Acoustics Engineering，Harbin Engineering University，Harbin 150001，China;2.Unit 92198 of Naval，Xingcheng 125109，China)

Abstract:A parametric array can generate an acoustic beam of low frequency and high directivity.This has been widely used in acoustics engineering including seafloor profile measurement，buried underwater target detection and recognition，and underwater communications.Theoretical studies have shown that the interaction of three collimated primary waves can produce a parametric array radiation system with two difference frequencies.However，no systematic analysis has been done of a multi－frequency parametric array acoustic field affected by the primary wave phase.To ensure the precision and calculation efficiency of the computation，this study drew on ideas from operator separation，combining the time domain with the frequency domain to solve the KZK equation for the parametric array，thereby achieving a description of a parameter array acoustic field formed by the interaction of three collimated primary waves with different phases.The simulation results showed that the sound pressure amplitude of the parametric array changed with a change in the primary wave phases.The sound pressure amplitude of the parametric array difference frequency component was largest when the primary wave phase was zero.Finally，the experimental study and the simulation produced identical results.

Keywords:parametric array;three collimated primary waves;primary wave phase;operator splitting scheme;time－frequency domain

通信作者:杨德森，E－mail:dshyang@ hrbeu.edu.cn.

作者简介:杨德森(1957－)，男，中国工程院院士，教授，博士生导师.

基金项目:水声技术国防科技重点实验室基金资助项目(9140c2001 0613c20078).

收稿日期:2014－10－28.网络出版时间:2015－12－21.

中图分类号:O427.9

文献标志码:A

文章编号:1006－7043(2016)01－0007－07

doi:10.11990/jheu.201410082

网络出版地址:http://www.cnki.net/kcms/detail/23.1390.u.20151221.1555.022.html