物理模型
固定深度
z
s
z_s
zs下,一个单频的点声源在水平方向上移动,当垂直阵与声源的水平距离为
r
i
r_i
ri时,记其接收到的一组声压列向量为
{
p
(
z
1
,
r
i
)
,
p
(
z
2
,
r
i
)
,
.
.
.
,
p
(
z
n
,
r
i
)
}
T
\{p(z_1,r_i),p(z_2,r_i),...,p(z_n,r_i)\}^T
{p(z1,ri),p(z2,ri),...,p(zn,ri)}T。改变水平距离,可以得到一个声压矩阵
P
{\rm P}
P:
简正波理论
根据简正波理论,理想波导的声场可以用如下式子表示:
p
(
r
,
z
)
=
∑
n
=
0
N
−
2
j
H
2
π
k
r
r
s
i
n
(
k
z
n
z
s
)
s
i
n
(
k
z
n
z
)
e
−
j
(
k
r
n
r
−
π
4
)
p(r,z)=\sum_{n=0}^{N}-\frac{2j}{H}\sqrt{\frac{2\pi}{k_rr}}sin(k_{zn}z_s)sin(k_{zn}z)e^{-j(k_{rn}r-\frac{\pi}{4})}
p(r,z)=n=0∑N−H2jkrr2πsin(kznzs)sin(kznz)e−j(krnr−4π)
根据声压的表达式形式,我们可以对声压矩阵进行分解:
P
=
e
j
π
4
Φ
Λ
R
{\rm P}=e^{j\frac{\pi}{4}}\Phi\Lambda{\rm R}
P=ej4πΦΛR
其中
Φ
\Phi
Φ对应各号简简正波的垂直模态函数,
Λ
\Lambda
Λ对应各号简正波的幅值,
R
{\rm R}
R对应
e
e
e指数项,下图为文献中的表示(文献中的简正波表达式与我们教材中的表达式存在系数上的差异,但不影响我们的仿真)
奇异值分解
我们现在已知的数据就是声压矩阵
P
{\rm P}
P,需要根据该矩阵提取得到深度模态函数,而深度模态函数的信息就在矩阵
Φ
\Phi
Φ,该矩阵的各个列向量是各号简正波深度模态函数在空间上的采样。
根据奇异值分解(SVD)的知识,一个矩阵可以进行如下方式的分解:
A
=
U
Λ
V
H
A=U\Lambda V^{H}
A=UΛVH
U
U
H
=
I
,
V
V
H
=
I
UU^{H}=I,VV^{H}=I
UUH=I,VVH=I
回到分解开来的声压矩阵:
1.
Φ
\Phi
Φ的各个列向量对应的是各号简正波的深度模态函数,故具有正交归一性,满足酉矩阵的定义;
2.
Λ
\Lambda
Λ是一个对角阵;
3.如果矩阵
R
{\rm R}
R也是个酉矩阵,那么就可以直接通过SVD分解得到
Φ
\Phi
Φ,
Λ
\Lambda
Λ和
R
{\rm R}
R。遗憾的是,
R
{\rm R}
R不是个酉矩阵,观察它的形式,发现如果把各列的
1
r
i
\sqrt{\frac{1}{r_i}}
ri1去掉,那么只剩下指数函数,各个行向量同样具有正交归一性;
所以引出接下来的操作,对声压矩阵
P
{\rm P}
P的各列进行能量归一化,即不同的列乘上一个不同的系数,使得各列的总能量(每一列的声压幅值的平方和)相同,满足:
∑
i
=
0
N
p
′
(
z
i
,
r
j
)
(
p
′
(
z
i
,
r
j
)
)
∗
=
1
\sum_{i=0}^{N}p^{'}(z_i,r_j)(p^{'}(z_i,r_j))^{*}=1
i=0∑Np′(zi,rj)(p′(zi,rj))∗=1
经过这个操作,那么
R
{\rm R}
R就变成了一个酉矩阵,此时我们就可以通过SVD的方法来分解得到
Φ
\Phi
Φ,
Λ
\Lambda
Λ和
R
{\rm R}
R了。matlab里的指令为svd()。至此,我们已经提取出了深度模态函数的信息,即矩阵
Φ
\Phi
Φ。
特征值分解
如果声源在水平距离上的取值过多的话(即列数过多),那么SVD将具有一个庞大的工作量,此时如果我们变换一下形式:
C
=
P
P
H
=
Φ
Λ
R
R
H
Λ
H
Φ
H
=
Φ
Λ
2
Φ
H
C={\rm P}{\rm P}^{H}=\Phi\Lambda{\rm R}{\rm R}^{H}\Lambda^{H}\Phi^{H}=\Phi\Lambda^{2}\Phi^{H}
C=PPH=ΦΛRRHΛHΦH=ΦΛ2ΦH
若
P
{\rm P}
P是一个
N
×
M
N\times M
N×M的矩阵
(
N
<
M
)
(N<M)
(N<M),则
C
C
C是一个
N
×
N
N\times N
N×N方阵,故可以进行特征值分解,特征值组成了矩阵
Λ
2
\Lambda^{2}
Λ2,特征向量组成了矩阵
Φ
\Phi
Φ。matlab里的指令为eig()。至此,同样能够提取出深度模态函数的信息,即矩阵
Φ
\Phi
Φ。
提取结果:
注意事项
1.水听器的个数
>
>
>简正波的号数,所以提取出来的特征值的个数
>
>
>简正波的号数,我们需要选取最大的那些特征值,这些才对应各个简正波;然后我们需要找出这些特征值对应的特征向量,这些对应的就是各号简正波的深度模态函数在垂直方向上的采样。
2.在作解析形式的深度模态函数与提取出来的深度模态函数的分析时:
- 如果遇到形状类似,但是存在一个正负号差异的情况,此时需要在该特征向量前面添上一个负号,因为在特征值分解时,特征向量之间仅需要满足正交归一性,所以方向可能会与解析形式的深度模态函数不同;
- 如果遇到幅值对不上的情况,注意解析形式的深度模态函数也是满足正交归一性的,其形式为 2 H s i n ( k z n z ) \sqrt{\frac{2}{H}}sin(k_{zn}z) H2sin(kznz),不要忘记了系数;
* 如何对应特征向量与简正波
首先需要明确的是,大的那些特征值对应的特征向量才对应不同号数的简正波,接近于0的那些特征值对应的特征向量无意义。
然后需要找到特征向量对应的那号简正波的模态函数,本来我是直接把特征向量和简正波都画出来,然后手动匹配,这样会导致程序的复用率低。通过和同学的讨论以及文献中的思路,设计了一种自动匹配的算法。
机理分析
- 观察分解得到的
Λ
\Lambda
Λ矩阵,是一个由特征值组成的对角阵,观察特征值,和简正波的幅值有着千丝万缕的联系;再观察能量归一化之前的矩阵
R
{\rm R}
R,前面乘上的
1
N
R
\frac{1}{\sqrt{N_R}}
NR1为了是保证能量归一化后的
R
{\rm R}
R各个行向量具有归一性,为了补偿,
Λ
\Lambda
Λ矩阵前面要乘上一个
N
R
\sqrt{N_R}
NR;
- Λ \Lambda Λ矩阵其他的系数来源于简正波的表达式,这里的 k i k_i ki是 i i i号简正波的水平波数;
- 如果仅仅想要消去矩阵 R {\rm R} R中的 1 r i \frac{1}{\sqrt{r_i}} ri1,仅需要对每列声压乘上 r i \sqrt{r_i} ri,但由于我们对声压矩阵的各列进行了能量归一化,所以实际上多乘了一个系数,需要把该系数乘到矩阵 Λ \Lambda Λ中去,系数 c c c: c = 1 ∑ i = 0 N z p ( z i , r j ) p ∗ ( z i , r j ) r i c=\frac{1}{\sum_{i=0}^{N_z} p(z_i,r_j)p^{*}(z_i,r_j)\sqrt{r_i}} c=∑i=0Nzp(zi,rj)p∗(zi,rj)ri1
- 提取出来的 Λ \Lambda Λ矩阵中的特征值= c 4 π N R H k r n s i n ( k z n z s ) c\sqrt{\frac{4\pi N_R}{Hk_{rn}}}sin(k_{zn}z_s) cHkrn4πNRsin(kznzs)
- 经过了一系列系数的补偿,特征值=简正波的幅值
算法思路
因为此时特征值与对应简正波幅值的位置还存在差异,所以需要进行位置匹配。
- t e m p temp temp=|特征值向量|-| i i i号简正波幅值|;
- 找到 t e m p temp temp向量中最小值出现的位置,该位置的特征值与第 i i i个位置的特征值交换位置,该位置的特征向量与第 i i i个位置的特征向量交换位置;
- 因为所有垂直模态在 0 + 0^{+} 0+的位置是正的,所以对位置匹配完成后的特征向量作判断,如果第一个元素<0,则该特征向量变号;
算法实现代码
%对提取出来的特征向量进行排序
%对特征向量的顺序按特征值与简正波幅值的差异进行排序
%g为有效的特征值组成的向量 t为各号简正波幅值组成的向量 V为各列特征向量组成的矩阵
for a=1:N+1
m=abs(abs(g)-abs(t(a)));
[value,index]=min(m);
if g(index)*t(a)<0%如果特征值和简正波幅值符号相反,那么对特征值取负号
g(index)=-g(index);
end
g([a,index])=g([index,a]);%交换特征值,使特征值与简正波号数对应
V(:,[a,index])=V(:,[index,a]);%交换特征向量
end
%这里的特征向量仅代表方向,在分解时可能会有负方向出现,故这里将其正过来
for a=1:N+1
if V(2,a)<0
V(:,a)=-V(:,a);%如果是负方向则翻转
end
end
Ref:
Extraction of acoustic normal mode depth functions using vertical line array data