声传播理论模型综述.docx

1、声传播理论模型综述声道中的声线绘制与频散方程数值求解作者：王翰卓1、模型内容给定大洋某处深度上离散采集的声速值，某深度处存在声速的极小值点，因为存在声道。认为在大洋一定范围内声速在水平方向上没有变化。解决以下问题-1、画出连续的声速剖面曲线2、当在适当深度上放置声源，使得当声波掠射角在之间时，所有的声线都会上下发生反转而不发生海底或海面的反射。应用射线声线的模型，画出其在给定距离上的声线。3、将射线声线与波动声学相联系，在只存在声线反转的情况下，给定声波频率，应用频散方程，求解特征波数。2、算法分析1、给定深度-声速曲线上的离散的点，可应用分段样条插值的方法得到深度上更多点的声速信息，同时保证

2、了曲线的平滑。2、算法如下所示 （1）输入传播距离，设置声源深度。 （2）对深度进行离散，对每个声层上的声速进行线性插值。 （3）计算声源处的声速值，用找出每条声线翻转处的声速，从声源处分开上下两层，分别用线性插值找出不同声线的上下翻转深度。 （4）对每条声线进行计算。先求出声线每层分层处的掠射角，在计算其所走的水平距离；每个角度值记录其从出发到其到达第一次翻转位置处的每层的水平位移。 （5）根据同一深度上某一声线的略射角相同，即声线的轴对称性，得到声线从第一次反转点到回到声源高度位置上的水平位移信息，最终得到一个跨度上的水平位移信息。 （6）进行画图。与的声线形状是一致的，故的上层可以用的上

3、层画出，的下层可以用的的下层画出；在利用每条声线的周期性，可将其延拓到指定水平距离。3、问题分析 波矢的方向即射线的某点的切线，水平方向的波矢具有不变性，一定时，由于，由声源极小值处发射的声线在向上下传播的过程中由于变小，相应与的夹角变小，直到即达到射线的反转点。为了使得所有的声线都反转，而且海底声速大于海面声速，故；同时因为，所以。越小，反转点离声源越远。当使得成立时，其为本征波数，显然，一定时，越小，越大，且积分上下限变长，简正波号数增大。增大时，简正波的最大号数也会增加。编程的算法如下：（1）给定，计算分层后的声速以及波数。（2）计算声速最小值，给出的范围。（3）用计算最大的简正波号数。

4、（4）用二分法在中找的根，积分的计算选用复化辛普生公式，最终确定每一号水平波数。3、源代码及运行结果1、声线绘制clc,clear all,format longT=tic;depth1=0.0,150.0,305.0,533.0,610.0,680.0;depth2=762.0,1372.0,1829.0,3048.0,4000.0;depth=depth1 depth2;c1=1507.2,1498.1,1491.7,1480.7,1478.9,1478.0;c2=1478.6,1483.2,1488.6,1507.5,1523.0;c=c1 c2;z0=700;%声源位置c0=inter

5、p1(depth,c,z0);x0=input(input the distance of propagation); z1=0:700;z2=700:4000;z=z1 z2;%深度上进行分层figure(1)plot(spline(depth,c,z),z,r);%样条差值后画出声速-深度曲线set(gca,yDir,reverse);set(gca,XAxisLocation,top);%y轴反转，x轴取有效部分title(sound speed section);xlabel(speed(m/s);ylabel(depth(m);grid on; c1n=interp1(depth1,

6、z0,c1,c0,z1);%线性差值，求封层后的声速c2n=interp1(z0,depth2,c0,c2,z2);theta=-pi/18:pi/180:-pi/180,pi/180:pi/180:pi/18;%入射角度cr=c0./cos(theta);reverse_up=interp1(c1n,z1,cr(1:10);%计算上下翻转深度reverse_down=interp1(c2n,z2,cr(11:20);x=zeros(20,2500);%存储每条声线的离散水平分量number=ceil(700-reverse_up,reverse_down-700)+1;%每条声线记录的基本x

7、点的个数,同时也是声线的高度x(:,1)=0;for i=1:20 if i=10 beta=acos(c1n(701:-1:701-number(i)+2)*cos(theta(i)/c0);%记录掠射角 beta(length(beta)+1)=0; else beta=acos(c2n(1:number(i)-1)*cos(theta(i)/c0); beta(length(beta)+1)=0; end for j=2:number(i) x(i,j)=x(i,j-1)+1/tan(beta(j-1)+beta(j)/2); end x(i,j+1:2*number(i)=2*x(i,

8、j)-fliplr(x(i,1:j);%x(i,number(i)+x(i,2*number(i)-j);%对称延拓两倍endfigure(2);%用对称性和周期性画图 for i=1:20 xmax=0; j=1; if(i=x0) break; end end else while(1)%向下发射 if rem(j,2)=0%上面的周期延拓 plot(xmax+x(21-i,1:2*number(21-i),-700:-700+number(21-i)-1,-700+number(21-i)-1:-1:-700); hold on; xmax=xmax+x(21-i,2*number(21

9、-i); else plot(xmax+x(i,1:2*number(i),-700:-1:-700-number(i)+1,-700-number(i)+1:-700); hold on; xmax=xmax+x(i,2*number(i); end j=j+1; if(xmax=x0) break; end end end axis(0,x0,-4000,1);grid on; hold on; end xlabel(distance/m),ylabel(depth/m),title(-pi/18-pi/18内的声线);toc(T) 2、频散方程clc;clear,format long;

10、T=tic;f=input(input the frequency);w=2*pi*f;%设定角频率%depth=0.0 150.0 305.0 533.0 610.0 680.0 762.0 1372.0 1829.0 3048.0 4000.0;%已知深度值c=1507.2 1498.1 1491.7 1480.7 1478.9 1478.0 1478.6 1483.2 1488.6 1507.5 1523.0;%已知声速值zn=0:4000;%深度上分层cn=interp1(depth,c,zn);%声速上线性插值c0=min(c);%找到最小声速n0=find(cn=c0);z0=zn

11、(n0);%最小声速对应的深度kz=w./cn;%波数分层khn_min=w/cn(1);%水平波数的最大最小值khn_max=w/c0; %找对应频率的最大简正波号数khn=khn_min;%最小水平波数时，积分等式左边值越大zmin=zn(1);%上翻转高度zmax=interp1(cn(n0:length(cn),zn(n0:length(zn),w/khn);%下翻转高度z_more=zmin:0.5:floor(zmax);%复化辛普生公式求积分值z_mid=z_more(2:2:length(z_more)-1);cnmid=interp1(zn,cn,z_mid);summ=4*

12、sum(sqrt(w./cnmid).2-khn2);summ=summ+2*sum(sqrt(kz(2:floor(zmax).2-khn2)+sqrt(kz(1)2-khn2)+sqrt(kz(floor(zmax)+1)2-khn2);eq_l=summ/6+trapz(sqrt(kz(floor(zmax)+1)2-khn2),0);%等式左边N=floor(eq_l/pi+1/2); %二分法找水平波数khn=zeros(1,N);for i=1:N%300hz,155hao if i=1 min=khn_min;max=khn_max; khn(i)=(min+max)/2; el

13、se min=khn_min;max=khn(i-1); khn(i)=(min+max)/2; end eq_r=2*i*pi; while (1) zmin=interp1(cn(1:n0),zn(1:n0),w/khn(i);%上翻转高度 zmax=interp1(cn(n0:length(cn),zn(n0:length(zn),w/khn(i);%下翻转高度 z_more=ceil(zmin):0.5:floor(zmax);%复化辛普生公式求积分值 z_mid=z_more(2:2:length(z_more)-1); cnmid=interp1(zn,cn,z_mid); sum

14、m=4*sum(sqrt(w./cnmid).2-khn(i)2); summ=summ+2*sum(sqrt(kz(ceil(zmin)+2:floor(zmax).2-khn(i)2)+sqrt(kz(floor(zmax)+1)2-khn(i)2)+sqrt(kz(ceil(zmin)+1)2-khn(i)2); eq_l=2*(summ/6+trapz(sqrt(kz(floor(zmax)+1)2-khn(i)2),0)+trapz(sqrt(kz(ceil(zmin)+1)2-khn(i)2),0)+pi;%等式左边 pluse=eq_l-eq_r; if abs(pluse)=1

15、0-8 break; else if pluse0 max=khn(i); else min=khn(i); end end khn(i)=(max+min)/2;% if abs(max-min)=2*eps break; end end sprintf(%d)=%f,i,khn(i)endtoc(T)4、结论和问题1、声线绘制，入射角大的声线翻转深度越深，跨度越大。由于声速和深度在整体深度上并非严格线性关系，所以跨度随角度的变化率、翻转深度随角度的变化率没有找到规律的联系。2、声线绘制，对称的分别向上下发射的角度，下翻转深度大于上翻转深度，说明了声源下层的声速变化率小于声源上层，这在声速-

16、深度剖面图中得到验证；这同时造成了上层的声线比下层的密，声波在上层传播时能量更集中。3、声线绘制，声源深度离声速极小值深度较远的情况，如500m,1300m. 会出现小掠射角度出现大跨度的情况。事实上，只有严格在声速极小值点发生的各条声线其从声源到第一次翻转过程中声线凹凸性是一致的，其余的都会出现拐点（数学上的凹凸性分界点，不是反转点）。只是当离声源较近时，这样的拐点不明显，当离声源较远时，如下图两种情况，会观察到明显的拐点，并造成上述现象。 4、频散方程，频散方程的实验结果验证了分析中得到的结论，同时发现增大，同一号简正波的也会变大，这是因为单纯的改变不会影响上下翻转深度，但是要保持不变，只

17、有增大。5、频散方程，随着简正波号数的增大，的变化率是在减小的。这说明本征声线在较大掠射角时会变得更为密集。6、频散方程，在程序刚开始调试过程中，只对方程两边相减的差设定了精度，发现了某些频率的特别号简正波陷入了死循环。参考张爽同志的程序，发现其也存在此问题-在时在第55号简正波处陷入死循环，这是一个很隐蔽的bug，想必学姐没发现。因为学姐采用的是边矩形公式求积，而此程序采用更高精度的辛普生公式求积，此程序更容易出现死循环。原因在于二分过程中当分到前后两个的值之差小于eps,即matlab可以分辨的最小正数时，二分变在计算机上便无效了，但此时方程两边的差可能还没有达到精度要求，所以出现死循环。数学表达为，调试中找到了这样的例子，严格的数学证明不再给出。所以在计算机上应用二分法求根时，只限定等式的精度是不合理的，而是应该设定二分的精度。此程序中引入了这一点，死循环的情况消失。