三列不同相位原波形成的參量陣聲場研究

楊德森，李中政，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