声传播理论模型综述.docx
《声传播理论模型综述.docx》由会员分享，可在线阅读，更多相关《声传播理论模型综述.docx（15页珍藏版）》请在冰豆网上搜索。
声传播理论模型综述
声道中的声线绘制与频散方程数值求解
作者：
王翰卓
1、模型内容
给定大洋某处深度上离散采集的声速值，某深度处存在声速的极小值点，因为存在声道。
认为在大洋一定范围内声速在水平方向上没有变化。
解决以下问题-
1、画出连续的声速剖面曲线
2、当在适当深度上放置声源，使得当声波掠射角在
之间时，所有的声线都会上下发生反转而不发生海底或海面的反射。
应用射线声线的模型，画出其在给定距离上的声线。
3、将射线声线与波动声学相联系，在只存在声线反转的情况下，给定声波频率，应用频散方程，求解特征波数。
2、算法分析
1、给定深度-声速曲线上的离散的点，可应用分段样条插值的方法得到深度上更多点的声速信息，同时保证了曲线的平滑。
2、算法如下所示
（1）输入传播距离，设置声源深度。
（2）对深度进行离散，对每个声层上的声速进行线性插值。
（3）计算声源处的声速值
，用
找出每条声线翻转处的声速，从声源处分开上下两层，分别用线性插值找出不同声线的上下翻转深度。
（4）对每条声线进行计算。
先求出声线每层分层处的掠射角
，在计算其所走的水平距离
；每个角度值记录其从出发到其到达第一次翻转位置处的每层的水平位移。
（5）根据同一深度上某一声线的略射角相同，即声线的轴对称性，得到声线从第一次反转点到回到声源高度位置上的水平位移信息，最终得到一个跨度上的水平位移信息。
（6）进行画图。
与
的声线形状是一致的，故
的上层可以用
的上层画出，
的下层可以用的
的下层画出；在利用每条声线的周期性，可将其延拓到指定水平距离。
3、问题分析
波矢的方向即射线的某点的切线，水平方向的波矢具有不变性，
一定时，由于
，由声源极小值处发射的声线在向上下传播的过程中由于
变小，相应
与
的夹角变小，直到
即达到射线的反转点。
为了使得所有的声线都反转，而且海底声速大于海面声速，故
；同时因为
，所以
。
越小，反转点离声源越远。
当
使得
成立时，其为本征波数，显然，
一定时，
越小，
越大，且积分上下限变长，简正波号数增大。
增大时，简正波的最大号数也会增加。
编程的算法如下：
（1）给定
，计算分层后的声速以及波数
。
（2）计算声速最小值
，给出
的范围
。
（3）用
计算最大的简正波号数
。
（4）用二分法在
中找
的根
，积分的计算选用复化辛普生公式
，最终确定每一号水平波数。
3、源代码及运行结果
1、声线绘制
clc,clearall,formatlong
T=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=[depth1depth2];
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=[c1c2];
z0=700;%声源位置
c0=interp1（depth,c,z0）;
x0=input（'inputthedistanceofpropagation'）;
z1=0:
700;z2=700:
4000;z=[z1z2];%深度上进行分层
figure
（1）
plot（spline（depth,c,z）,z,'r'）;%样条差值后画出声速-深度曲线
set（gca,'yDir','reverse'）;set（gca,'XAxisLocation','top'）;%y轴反转，x轴取有效部分
title（'soundspeedsection'）;
xlabel（'speed（m/s）'）;ylabel（'depth（m）'）;
gridon;
c1n=interp1（[depth1,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点的个数,同时也是声线的高度
x（:
1）=0;
fori=1:
20
ifi<=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
forj=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,j）-fliplr（x（i,1:
j））;%x（i,number（i））+x（i,2*number（i）-j）;%对称延拓两倍
end
figure
（2）;%用对称性和周期性画图
fori=1:
20
xmax=0;
j=1;
if（i<=10）%向上发射
while
（1）
ifrem（j,2）==0%下面周期的延拓
plot（xmax+x（21-i,1:
2*number（21-i））,[[-700:
-1:
-700-number（21-i）+1],[-700-number（21-i）+1:
-700]]）;holdon;
xmax=xmax+x（21-i,2*number（21-i））;
else%上面周期的延拓
plot（xmax+x（i,1:
2*number（i））,[[-700:
-700+number（i）-1],[-700+number（i）-1:
-1:
-700]]）;holdon;
xmax=xmax+x（i,2*number（i））;
end
j=j+1;
if（xmax>=x0）
break;
end
end
else
while
（1）%向下发射
ifrem（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]]）;holdon;
xmax=xmax+x（21-i,2*number（21-i））;
else
plot（xmax+x（i,1:
2*number（i））,[[-700:
-1:
-700-number（i）+1],[-700-number（i）+1:
-700]]）;holdon;
xmax=xmax+x（i,2*number（i））;
end
j=j+1;
if（xmax>=x0）
break;
end
end
end
axis（[0,x0,-4000,1]）;gridon;
holdon;
end
xlabel（'distance/m'）,ylabel（'depth/m'）,title（'-pi/18-pi/18内的声线'）;
toc（T）
2、频散方程
clc;clear,formatlong;
T=tic;
f=input（'inputthefrequency'）;
w=2*pi*f;%设定角频率%
depth=[0.0150.0305.0533.0610.0680.0762.01372.01829.03048.04000.0];%已知深度值
c=[1507.21498.11491.71480.71478.91478.01478.61483.21488.61507.51523.0];%已知声速值
zn=0:
4000;%深度上分层
cn=interp1（depth,c,zn）;%声速上线性插值
c0=min（c）;%找到最小声速
n0=find（cn==c0）;
z0=zn（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*sum（sqrt（（w./cnmid）.^2-khn^2））;
summ=summ+2*sum（sqrt（kz（2:
floor（zmax））.^2-khn^2））+sqrt（kz
（1）^2-khn^2）+sqrt（kz（floor（zmax）+1）^2-khn^2）;
eq_l=summ/6+trapz（[sqrt（kz（floor（zmax）+1）^2-khn^2）,0]）;%等式左边
N=floor（eq_l/pi+1/2）;
%二分法找水平波数
khn=zeros（1,N）;
fori=1:
N%300hz,155hao
ifi==1
min=khn_min;max=khn_max;
khn（i）=（min+max）/2;
else
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）;
summ=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;
ifabs（pluse）<=10^-8
break;
else
ifpluse<0
max=khn（i）;
else
min=khn（i）;
end
end
khn（i）=（max+min）/2;%
ifabs（max-min）<=2*eps
break;
end
end
sprintf（'ξ（%d）=%f',i,khn（i））
end
toc（T）
4、结论和问题
1、声线绘制，入射角大的声线翻转深度越深，跨度越大。
由于声速和深度在整体深度上并非严格线性关系，所以跨度随角度的变化率、翻转深度随角度的变化率没有找到规律的联系。
2、声线绘制，对称的分别向上下发射的角度，下翻转深度大于上翻转深度，说明了声源下层的声速变化率小于声源上层，这在声速-深度剖面图中得到验证；这同时造成了上层的声线比下层的密，声波在上层传播时能量更集中。
3、声线绘制，声源深度离声速极小值深度较远的情况，如500m,1300m.会出现小掠射角度出现大跨度的情况。
事实上，只有严格在声速极小值点发生的各条声线其从声源到第一次翻转过程中声线凹凸性是一致的，其余的都会出现拐点（数学上的凹凸性分界点，不是反转点）。
只是当离声源较近时，这样的拐点不明显，当离声源较远时，如下图两种情况，会观察到明显的拐点，并造成上述现象。
4、频散方程，频散方程的实验结果验证了分析中得到的结论，同时发现
增大，同一号简正波的
也会变大，这是因为单纯
的改变不会影响上下翻转深度，但是
要保持不变，只有增大
。
5、频散方程，随着简正波号数的增大，
的变化率是在减小的。
这说明本征声线在较大掠射角时会变得更为密集。
6、频散方程，在程序刚开始调试过程中，只对方程两边相减的差设定了精度，发现了某些频率的特别号简正波陷入了死循环。
参考张爽同志的程序，发现其也存在此问题-在
时在第55号简正波处陷入死循环，这是一个很隐蔽的bug，想必学姐没发现。
因为学姐采用的是边矩形公式求积，而此程序采用更高精度的辛普生公式求积，此程序更容易出现死循环。
原因在于二分过程中当分到前后两个
的值之差小于eps,即matlab可以分辨的最小正数时，二分变在计算机上便无效了，但此时方程两边的差可能还没有达到精度要求，所以出现死循环。
数学表达为
，调试中找到了这样的例子，严格的数学证明不再给出。
所以在计算机上应用二分法求根时，只限定等式的精度是不合理的，而是应该设定二分的精度。
此程序中引入了这一点，死循环的情况消失。
升级会员 