矩形孔径参量阵相控非线性声场建模与实验研究

朱建军，李海森，魏玉阔，陈宝伟，周 天

（1．哈尔滨工程大学水声技术重点实验室，哈尔滨 150001；2．哈尔滨工程大学水声工程学院，哈尔滨 150001）

矩形孔径参量阵相控非线性声场建模与实验研究

朱建军1，2，李海森1，2，魏玉阔1，2，陈宝伟1，2，周 天1，2

（1．哈尔滨工程大学水声技术重点实验室，哈尔滨 150001；2．哈尔滨工程大学水声工程学院，哈尔滨 150001）

为了解决一阶抛物近似KZK方程无法直接对相控非线性声场进行建模的问题，依据原频声场准直特性构建等效参量阵原频自然指向性声场模型，建立其与原频相控声场的等效关系，将相控非线性声场的建模问题转换为等效参量阵声源条件的求解问题，同时提出算子分裂时域有限差分数值计算方法，实现了基于KZK抛物方程的矩形孔径参量阵相控非线性声场数值计算。以SES2000标准型参量声呐辐射相控非线性声场为研究对象，开展了计算机仿真和水池对比实验研究，研究结果验证了构建声场模型及其数值计算方法的有效性。

矩形孔径参量阵；相控非线性声场；声场建模

近年来，随着非线性声学理论的不断发展，声波的非线性效应被逐步应用于水声工程、生物医学和空气声学等领域［1－3］。当参量阵同时辐射两种频率的高强度有限振幅波时，声场中会产生声散射声—差频波、和频波与谐频波［4］。这一非线性效应有效解决了传统小孔径基阵无法相控发射低频、宽带、无旁瓣窄波束的问题［5］，使声学换能器突破了瑞利限的限制。基于这一技术特点，水声中参量声呐的差频探测技术在浅地层剖面探测［6］、浅海远程通信［7］、海底掩埋物探测［8］等领域获得了传统线性声呐无法比拟的优越性能。

参量声呐科学设计与合理应用均需以非线性声场（参量声场）研究为前提，涉及声场建模、计算以及相关实验验证研究。除Westervelt方程［9］、Burgers方程［10］等经典方程，KZK（Khokhlov Zabolotskaya Kuznetsov）方程［11］是目前为止最为精准的参量声场建模方程［12］。这是由于实际参量阵的孔径均具有一定的形状和大小，必须考虑声波衍射效应；所有介质均具有热传导性和黏滞性，介质对声波的吸收效应也不可忽略，上述三种建模方程中仅KZK方程全面考虑了声波的衍射效应、声吸收效应以及非线性效应。KZK方程在用于参量声场研究时既具有其自身特点，也存在一些尚未有效解决的问题：① 由于圆形孔径参量阵具有轴对称的声场分布特征，易于建模［13］，从研究之初就多针对这一声场结构开展研究，而矩形孔径参量阵辐射声场研究相对薄弱；② 以应用需求为导向，针对生物医学非线性谐波场的研究较为深入［14］，针对水声中广泛应用的差频声场研究相对较少；③ 受KZK方程一阶抛物近似仅能描述主轴附近±15°范围内参量声场的限制，无法采用KZK方程直接对相控参量声场进行描述［15］。

针对近年来水声领域对矩形孔径相控参量阵差频探测技术的迫切需求以及相控参量声场研究中存在的上述问题，从矩形孔径参量声场建模与数值计算、相控参量声场建模以及实验验证等几个方面开展重点研究，在对矩形孔径参量阵自然指向性参量声场建模的前提下，针对一级抛物近似KZK方程无法直接对相控参量声场进行建模的问题建立等效参量阵坐标系，通过基阵坐标变换将相控参量声场的建模问题转变为等效参量阵声源条件的求解问题，并在等效参量阵坐标系下实现相控参量声场的数值计算。

1 自然指向性参量声场建模与数值计算

1.1 非扩散参量声场KZK建模方程

由非线性声学基本方程推得近似至二阶微小项的非线性波动方程为［16］：

式中：p为声压，c0为小信号声速，ρ0为介质密度，ξ为吸收系数；β为非线性系数。由Muir平方根展开关系［17］，取一级抛物近似并做t′＝t－z／c0的参数变换，得到笛卡尔坐标系下的KZK方程为：

式中：δ为介质热粘吸收损耗因数，∇2⊥为笛卡尔坐标系下的拉普拉斯算子。公式中右侧三项因式依次表示声衍射、声吸收和非线性效应。为便于在笛卡尔坐标系中描述参量声场，定义矩形孔径参量阵坐标系见图1，其中基阵长、宽分别为2a和2b（a≥b），以辐射面中心点为坐标原点，基阵辐射面与xy平面重合，z为基阵法线方向，同时定义孔径纵横比R＝b／a≤1。

图1 矩形孔径参量阵坐标系Fig．1 Rectangularaperture parametric array coordinate

由于超出菲涅尔区后声波将以球面扩展方式传播，声场区域也随之扩大，从而声场计算量也会随之增加，声场建模需充分考虑这一问题，以获得优化的声场建模方案。Aanonsen等［18］提出参数去量纲化法实现扩散声场至非扩散声场的变换，达到了降低声场计算量的目的，KZK方程中参量的去量纲化形式为：σ＝z／R0，ux＝x／［a（1+σ）］，uy＝y／［a（1+σ］，P＝p（1+ σ）／p0，τ＝ω0t′－k0（x2+y2）／［2R0（1+σ）］，其中R0为瑞利距离，k0为原频波数，t′＝t－z／c0，p（x，y，z，t）为（x，y，z）场点处声压，p0为辐射面特征声压。将去量纲化参量代入式（2）并对等式两侧做时间积分，推得非扩散声场中的KZK抛物方程又称TBE（Transformed Beam E-quation）方程［19］：

式中：A＝α0R0为声吸收系数，α0为衰减系数，N＝R0／zs为非线性参数，zs为冲击距离。就基于式（3）在非扩散声场中对参量声场进行数值计算。

1.2 参量阵声源条件

参量阵辐射面上辐射声波的声压值为：

p（x，y，z＝0，t）＝p0·e（t）·s（t）·H（x，y）（4）式中：p0为初始声压幅值，e（t）为信号包络，s（t）为载频信号，H（x，y）为声源分布函数，其去量纲化表示形式为：

1.3 声场数值计算

由于KZK方程是高阶抛物方程，难以求出一般形式的封闭解析解，因此只能采用数值计算方法对声场进行计算。目前，KZK方程数值求解方法主要包括时域法和频域法［20］，频域法将声场表示为傅里叶级数展开各次谐波叠加的形式，但实际选取谐波阶次有限（最大谐波阶次称截断谐波次数），有限的截断谐波次数阻断了原频波向高次谐波以及高次谐波向低次谐波的能量传递，且在截断谐波频率处形成能量堆积，计算误差逐步增大。在计算短脉冲或宽带信号激励参量声场时，频域法必须保留大量频谱成分，而其计算量与截断谐波次数的平方成比例，从而计算量将很大，并且频域法难以对相控声场等具有非一致声源条件的声场进行建模。时域法不仅可以使用非一致性声源条件，而且基于波形的研究更加直观，因此采用时域法对参量声场进行数值计算。

定义τ时刻非扩散参量声场中场点（ux，uy，σ）处声压的离散表达形式为：

式中，Δux、Δuy分别为x，y向变化步长，Δτ时间变化步长，波束主轴方向传播距离σk＝∑Δσl，l∈（0，k－1），Δσl＝（1+σk）2Δσ0，σ0为初始步长。从而，数值积分约束区间为j∈［jmin，jmax］、m∈［mmin，mmax］和n∈［nmin，nmax］。同时定义声场的空间和时间零声压边界条件分别为：

由于矩形孔径参量阵自然指向性声场在x轴和y轴方向上均对称，因此计算1／4辐射声场区域即可对整个声场空间进行描述，此时引入附加边界条件：

取m∈［0，mmax］和n∈［0，nmax］。但即便如此，矩形孔径参量阵的声场计算仍是针对三维空间的数值计算，而圆形孔径参量阵的声场计算区域是基阵半径和声轴所在的二维平面，这正是矩形孔径参量阵声场研究的复杂性和难点所在。借鉴Lee［21］采用的算子分裂法求解思想，提出矩形孔径自然指向性参量声场时域有限差分数值计算方法，算法结构见图2。

图2 时域有限差分数值计算方法Fig．2 Time domain finite difference numerical calculation method

在声源附近（z≤zc，zc＝0．1R0），采用距离上一阶精度近似的隐式后向有限差分（Implicit Backward Finite Difference，IBFD）格式，以较小的距离变化步长消除声源条件中阶跃函数不连续引起的近场数值振荡，提高近场计算精度；在此区域外，采用绝对稳定的二阶精度近似Crank-Nicolson有限差分（Crank-Nicolson Finite Difference，CNFD）格式，在保证算法稳定性的同时可逐渐增大步长进一步减小声场计算量。

2 相控参量声场建模

2.1 等效参量阵坐标系

一级抛物近似KZK方程仅能有效描述声轴±15°范围内的参量声场，对具有超高指向性（通常＜5°）的参量阵来说，足以满足其自然指向性声场的建模要求，但在描述较大相控角度的声场时将失效［22］。

由参量阵基本理论可知，参量阵辐射面附近＜（L1＝4b2／λ0）的区域内原频波具有准直特性，该区域内垂直声轴的声场截面与基阵辐射面的大小和形状完全相同，而该区域至菲涅尔区外沿的声场中原频波做柱形扩展传播［23］，传播距离L2＝4（a2－b2）／λ0，矩形孔径参量阵原频声场分布见图3。

图3 矩形孔径参量阵原频声场分布Fig．3 Rectangular parametric array primary field distribution

对于相控声场，波束主极大偏离了基阵法线方向，此时，提出等效参量阵的概念，相控参量声场被认为由该等效参量阵产生，且此等效参量阵原频声场仍满足准直特性。建立等效参量阵二维坐标系xp－zp，使zp轴与等效参量阵法线重合（见图4），从而相控参量声场也可由该坐标系进行描述。

图4 等效参量阵坐标系Fig．4 Equivalent parametric array coordinate

由图4的几何结构推导出同一场点F在x－z和xp－zp两坐标系中坐标（x，z）和（xp，zp）间的变换关系为：

2.2 等效参量阵声源条件求解

由等效参量阵原频准直特性可知，xp－zp坐标系符合KZK方程声场建模的条件，只要根据x－z和xp－zp两坐标系间的坐标变换关系求解等效参量阵声源条件，即可在xp－zp坐标系中采用非扩散参量声场计算方法实现相控参量声场数值计算，将相控声场建模问题转换为等效参量阵声源条件的求解问题。由于声源附近非线性作用刚刚开始，声场中最主要的频率成分仍为原频，声散射声信号能量特别微弱，可忽略不计。此时声源附近原频声场可近似为线性声场，从而等效参量阵的声源条件可近似由线性声场计算方法求解得到。

假设参量阵辐射面周围为刚性障板，基于矩形增量区域划分模型计算原频声场，声场建模原理见图5，声场的离散形式表示为［24］：

其中，N为划分面元个数，Δs为面元面积，r为场点至面元的距离，un为第n个面元的表面复速度。对于一致激励声源，各面元表面复速度相同，而对于非一致激励的情况（如相控基阵），不同面元具有不同的un。

图5 原频声场计算模型Fig．5 Primary sound field calculation model

同时，推得相控角度为θp时等效参量阵辐射面上各场点在x－z坐标系中的坐标为：

式中：xL∈［－a，a］为横坐标，yw∈［－b，b］为纵坐标，yfield＝yw。因此，式（12）可用于等效参量阵声源条件求解。

相控角度为0°和15°时110 kHz原频波束上垂直声轴的声场剖面如图6所示，参量阵孔径为0．22 m× 0．22 m，dp＝6．7 cm（L1＝3．55 m）。0°相控声场剖面在x，y两维方向上均对称，而15°相控声场主轴在x轴向发生了相应偏移，但在y轴向仍对称。因此，相控参量声场数值计算时xp、yp两维方向上积分区间可定义为m∈［m′min，m′max］和n∈［0，nmax］，m′min、m′max对应等效参量阵xp轴向孔径，即相控参量声场可用yp轴方向上对称声场的1／2自然指向性声场来描述。

图6 原频准直区域内声场剖面Fig．6 Acoustic field profiles in thecollimate zone

3 水池实验验证分析

3.1 实验方案

水池实验在水声技术重点实验室的信道水池进行，水池大小为45 m×6 m×5 m，四面侧壁均铺设消声尖劈，吸声截止频率50～1 MHz，仅在19．2 kHz时存在一定的背景噪声干扰，但低于7 dB。基本实验系统由德国Innomar公司SES2000标准型参量声呐（原频中心频率100 kHz，差频4 kHz、5 kHz、6 kHz、8 kHz、10 kHz、12 kHz、15 kHz）、Reson TC4038标准水听器（带宽10～800 kHz，灵敏度－227dBre 1V／μPa）、TOPVIEW2000多通道采集器（最高采样频率25 MHz）和HY 1200声速剖面仪构成（用于水池中声速剖面测量，为声场建模提供准确的声速信息）。参量阵布放水深1．5 m，经刚性转接架装于可旋转的吊装杆上，水听器布放于小行车上，并处于同一深度。通过声轴上参量信号波形测试、原频和差频相控指向性声场测试以及轴向声场测试，与相同参数的计算机仿真结果进行比对，对本研究矩形孔径参量阵相控参量声场建模与数值计算方法的有效性进行验证，实验方案（见图7）。

3.2 实验结果及分析

图7 实验方案设计图Fig．7 Diagram of the experimental scheme

SES2000参量阵孔径为0．22 m×0．22 m，实测1．5 m处水中声速为1 474m／s。相控参量声场建模和数值计算在ThinkStation 4158－FB9工作站上进行，声场建模参量为：a＝0．11 m，b＝0．11 m，β＝3．6，原频声源级239 dB，Δx＝Δy＝0．03，ux＝uy＝10，IBFD计算时σ＝0．001（1+σ）2，CNFD计算时σ＝0．003（1+σ）2，fs＝3．075 MHz。水池实验测量数据处理时对标准水听器不同频率的灵敏度进行了一致性补偿处理。

12 m处仿真与实测参量信号波形及频谱如图8所示，差频信号频率为12 kHz。实测两原频信号的频谱幅度略微不同且信号波形出现了一定的振铃效应，分析原因主要是由于实际声呐系统响应难以做到完全一致且实测信号中混入了噪声干扰（电噪声和环境噪声等），而系统带宽有限是引起振铃效应的主要因素。但对比可以看出，总体上仿真与实测结果吻合得较好，初步验证了声场建模与数值计算方法的有效性。

图8 参量信号波形及频谱Fig．8Waveform and spectrum of parametric signal

同时，仿真了7 m和12 m处0°∶5°∶30°七个相控角度的差频指向性，并与实测结果进行对比（见图9）。由于SES2000最大相控角度为16°，因此仅测试了0°∶5°∶15°四个相控角度的指向性。图中差频主波束以外角度上噪声水平高于理论仿真结果，这是由于在水池环境噪声及实验系统电噪声等干扰因素影响下，随着测量角度增大信噪比逐渐降低造成的，在±10°范围外信噪比低于实验系统的检测阈，从而表现为主瓣外实测值（噪声）高于仿真指向性。

最后，对95 kHz原频信号以及12 kHz和15 kHz差频信号轴向声场的仿真与实测结果进行了比对（见图10），两者吻合较好，测量结果中差频曲线出现了微小的起伏，分析其原因是由于差频信号较弱，受环境噪声、电噪声以及小行车行走过程中引起的振动等瞬态干扰造成的。

图9 差频指向性图Fig．9 Directivity of differency frequency

图10 轴向声场Fig．10 Axialacoustic field

由对比实验可知，参量声场的仿真与实测结果整体上基本一致，充分验证了本文提出矩形孔径相控参量声场建模与数值计算方法的正确性和有效性。

4 结 论

系统地研究了矩形孔径参量阵相控参量声场的建模与数值计算方法，推导了矩形孔径参量声场的KZK建模方程，采用时域有限差分法实现了声场的数值计算，解决了基于KZK方程的相控参量声场的建模问题。参量阵原频声场的准直特性支持了参量阵相控参量声场等效于其等效参量阵自然指向性声场的结论，为利用一阶近似KZK抛物方程进行相控参量声场建模提供了理论依据。同时，声源条件求解问题是相控参量声场建模与数值计算的关键问题，相控参量阵非一致性声源条件的特点决定了时域有限差分法更适合于相控参量声场数值计算。计算机仿真结果与SES2000标准型参量声呐辐射相控参量声场吻合得很好，验证了声场建模和数值计算的正确性和有效性，研究成果可为相控参量声呐设计及差频海底探测研究提供理论依据和技术支持，对生物医学和空气声学等其它领域的相控参量声场研究也具有借鉴意义。

致谢：感谢水声技术重点实验室纪祥春老师、田宝晶老师以及李磊老师在开展水池实验过程中给予的大力支持和无私帮助。

［1］Olav R G，Kenneth G F，Johnny D，et al．Detecting atlantic herring by parametric sonar［J］．J．Acoust．Soc．Am，Express Letters，2010，127（4）：153－159．

［2］Egor V D，Bojan B G．Acoustic radiation force in tissue-like solids due to modulated sound field［J］．Journal of the Mechanics and Physics of Solids，2012，60：1791－1813．

［3］Akiyoshi I，Nobuo T，Hiroyuki I．Global active noise control using a delayed driving parametric array loudspeaker［J］．The Japan Society of Mechanical Engineers，2013，79（799）：604－617．

［4］Milan C，Michal B．Non-paraxialmodel for a parametric acoustic array［J］．J．Acoust．Soc．Am，2013，134（2）：933－938．

［5］Hideyuki N，Claes MH，Tomoo K．Numerical simulation of length-limited parametric sound beam［C］．International Congress on Ultrasonics，2012，1433：547－550．

［6］David G．B，Mark B M，William L K．Parametric acoustic array developmentat the USNavy's New London，Connecticut laboratory［C］．Acoustical Society of America，Physical A-coustics，2009，6：1－22．

［7］Karsten W，Thomas B，Tobias W．Parametric underwater communications［C］．11thEuropean Conference on Underwater Acoustics，2012，17：1－10．

［8］Akira A，TamakiU，Fumitaka M，etal．Sub-bottom synthetic aperture imaging sonar system using an AUV and an autonomous surface tracking vehicle for searching for buried shells of toxic chemicals［C］．Waterside Security conference，2010，1－3．

［9］Maxim S，Tony W H S．Simulation of nonlinear westervelt equation for the investigation of acoustic streaming and nonlinear propagation effects［J］．J．Acoust．Soc．Am，2013，134（5）：3931－3942．

［10］杨德森，兰朝凤，时胜国，等．水中声波非线性相互作用的声吸收研究［J］．振动与冲击，2012，31（8）：52－56．

YANG De-sen，LAN Chao-feng，SHI Sheng-guo，et al．Sound absorption of sound under interaction among underwater nonlinear acoustic variable parameters［J］．Journal of Vibration and Shock，2012，31（8）：52－56．

［11］LIUWei，YANG Jun，TIAN Jing．Time-domain modeling of finite-amplitude sound beams in three-dimensional Cartesian coordinate system［J］．Chinese Journal of Acoustics，2012，31（4）：408－422．

［12］杜宏伟．生物医学超声中若干非线性问题的研究［D］．合肥：中国科学技术大学，2007．

［13］Kenneth G F，David T IF，Philip R A．Calibration sphere for low-frequency parametric sonars［J］．J．Acoust．Soc．Am，2007，121（3）：1482－1490．

［14］Paul LMJvan N，MikhailG D，Martin D V，etal．Comparison of fundamental，second harmonic and superharmonic imaging：A simulation study［J］．J．Acoust．Soc．AM，2011，130（5）：3148－3157．

［15］Shi Chuang，Gan W S．Steerable parametric loudspeakerwith preprocessing methods［C］．Proceedings of Meetings on A-coustics，2013，19：1－6．

［16］李太宝．计算声学－声场的方程和计算方法［M］．北京：科学出版社，2005．

［17］Yang Jun，Sha Kan，Gan W S，etal．A fast field scheme for the parametric sound radiation from rectangular aperture source［J］．Chin．Phys．Lett，2004，21（1）：111－113．

［18］Aanonsen S I，Barkve T，Tjøtta J N，et al．Distortion and harmonic generation in the nearfield of a finite amplitude sound beam［J］．J．Acoust．Soc．Am，1984，75：749－768．

［19］Kaya O A，sahin A，Kaleci D．Pressure field of rectangular transducers at finite amplitude in three dimensions［J］．Ultrasound in Med．＆Biol，2006，32（2）：271－280．

［20］Egor V D，Bojan BG．On the KZK-type equation formodulated ultrasound fields［J］．Wave Motion，2013，50：763－775．

［21］Lee Y S．Numerical solution of the KZK equation for pulsed finite amplitude sound beams in thermoviscous fluids［D］．Texas：University of Texas，1993．

［22］MarcoMV．3D harmonic echocardiography［D］．Rotterdam：Erasmus University，2007．

［23］Lucilla D M，Jacques M，Pierre C．Nonlinearmulti-frequency generation for underwater application［J］．Applied Acoustic，2012，73：900－903．

［24］Kenneth BO，Leon A F．Sound field calculation for rectangular sources［J］．IEEE Transactions on Ultrasonics，1989，36（2）：242－248．

Modeling and experimental investigation on rectangular aperture parametric array phased nonlinear acoustic field

ZHU Jian-jun1，2，LIHai-sen1，2，WEIYu-kuo1，2，CHEN Bao-wei1，2，ZHOU Tian1，2
（1．Acoustic Science and Technology Laboratory，Harbin Engineering University，Harbin 150001，China；2．College of Underwater Acoustic Engineering，Harbin Engineering University，Harbin，150001，China）

In order to solve the problem of that phased nonlinear acoustic field could not bemodeled by one order parabolic approximation KZK equation，an equivalent parametric array primary frequency natural directivity acoustic field modelwas builtbased on the collimation characteristic of primary frequency acoustic field．The equivalent relation between this field and the primary frequency phased acoustic field was deduced，so the problem of phased nonlinear acoustic field modeling was transferred to the calculation of equivalent parametric array source condition，and an operator split time domain finite difference calculation method was derived．The rectangular aperture parametric array phased nonlinear field calculation based on KZK equation was realized．Taking the SES2000 parametric array acoustic field as object，the phased nonlinear acousticmodel and the corresponding field calculation method were examined by computer simulation and pool contrast experiment．

rectangular aperture parametric array；phased nonlinear acoustic field；acoustic field modeling

TB532；TB52+5

A

10．13465／j．cnki．jvs．2015．12．005

国家自然科学基金（41327004，41306182）；高等学校博士学科点专项科研基金（20112304130003）；中央高校基本科研业务费（HEUCF140501）

2014－03－25 修改稿收到日期：2014－05－14

朱建军 男，博士生，助教，1981年生

李海森 男，博士，教授，博士生导师，1962年生