声传播理论模型综述文档格式.docx
《声传播理论模型综述文档格式.docx》由会员分享，可在线阅读，更多相关《声传播理论模型综述文档格式.docx（15页珍藏版）》请在冰点文库上搜索。
的声线形状是一致的，故
的上层可以用
的上层画出，
的下层可以用的
的下层画出；
在利用每条声线的周期性，可将其延拓到指定水平距离。
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'
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/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）;
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:
xmax=0;
j=1;
if（i<
=10）%向上发射
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:
-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:
xmax=xmax+x（i,2*number（i））;
j=j+1;
if（xmax>
=x0）
break;
while
（1）%向下发射
ifrem（j,2）==0%上面的周期延拓
-700+number（21-i）-1],[-700+number（21-i）-1:
holdon;
-700-number（i）+1],[-700-number（i）+1:
axis（[0,x0,-4000,1]）;
xlabel（'
distance/m'
）,ylabel（'
depth/m'
）,title（'
-pi/18-pi/18内的声线'
toc（T）
2、频散方程
clc;
clear,formatlong;
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:
%深度上分层
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）;
N%300hz,155hao
ifi==1
min=khn_min;
max=khn_max;
khn（i）=（min+max）/2;
max=khn（i-1）;
eq_r=2*i*pi;
while
（1）
zmin=interp1（cn（1:
n0）,zn（1:
n0）,w/khn（i））;
zmax=interp1（cn（n0:
length（zn））,w/khn（i））;
z_more=ceil（zmin）:
z_mid=z_more（2:
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
ifpluse<
max=khn（i）;
min=khn（i）;
khn（i）=（max+min）/2;
%
ifabs（max-min）<
=2*eps
sprintf（'
ξ（%d）=%f'
i,khn（i））
toc（T）
4、结论和问题
1、声线绘制，入射角大的声线翻转深度越深，跨度越大。
由于声速和深度在整体深度上并非严格线性关系，所以跨度随角度的变化率、翻转深度随角度的变化率没有找到规律的联系。
2、声线绘制，对称的分别向上下发射的角度，下翻转深度大于上翻转深度，说明了声源下层的声速变化率小于声源上层，这在声速-深度剖面图中得到验证；
这同时造成了上层的声线比下层的密，声波在上层传播时能量更集中。
3、声线绘制，声源深度离声速极小值深度较远的情况，如500m,1300m.会出现小掠射角度出现大跨度的情况。
事实上，只有严格在声速极小值点发生的各条声线其从声源到第一次翻转过程中声线凹凸性是一致的，其余的都会出现拐点（数学上的凹凸性分界点，不是反转点）。
只是当离声源较近时，这样的拐点不明显，当离声源较远时，如下图两种情况，会观察到明显的拐点，并造成上述现象。
4、频散方程，频散方程的实验结果验证了分析中得到的结论，同时发现
增大，同一号简正波的
也会变大，这是因为单纯
的改变不会影响上下翻转深度，但是
要保持不变，只有增大
5、频散方程，随着简正波号数的增大，
的变化率是在减小的。
这说明本征声线在较大掠射角时会变得更为密集。
6、频散方程，在程序刚开始调试过程中，只对方程两边相减的差设定了精度，发现了某些频率的特别号简正波陷入了死循环。
参考张爽同志的程序，发现其也存在此问题-在
时在第55号简正波处陷入死循环，这是一个很隐蔽的bug，想必学姐没发现。
因为学姐采用的是边矩形公式求积，而此程序采用更高精度的辛普生公式求积，此程序更容易出现死循环。
原因在于二分过程中当分到前后两个
的值之差小于eps,即matlab可以分辨的最小正数时，二分变在计算机上便无效了，但此时方程两边的差可能还没有达到精度要求，所以出现死循环。
数学表达为
，调试中找到了这样的例子，严格的数学证明不再给出。
所以在计算机上应用二分法求根时，只限定等式的精度是不合理的，而是应该设定二分的精度。
此程序中引入了这一点，死循环的情况消失。
升级会员 