【简正波作业】——深度模态函数提取

物理模型

　固定深度$z_s$下，一个单频的点声源在水平方向上移动，当垂直阵与声源的水平距离为$r_i$时，记其接收到的一组声压列向量为 { 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。改变水平距离，可以得到一个声压矩阵$\mathrm{P}$：
　

简正波理论

  根据简正波理论，理想波导的声场可以用如下式子表示：
$$p(r,z)=\sum _{n=0}^N-\frac{2j}{H}\sqrt{\frac{2\pi }{k_{r}r}}sin(k_{zn}{z}_s)sin(k_{zn}z)e^{-j(k_{rn}r-\frac{\pi }{4})}$$
　根据声压的表达式形式，我们可以对声压矩阵进行分解：
　$$\mathrm{P}=e^{j\frac{\pi }{4}}\Phi \Lambda \mathrm{R}$$
　其中$\Phi$对应各号简简正波的垂直模态函数，$\Lambda$对应各号简正波的幅值，$\mathrm{R}$对应$e$指数项，下图为文献中的表示(文献中的简正波表达式与我们教材中的表达式存在系数上的差异，但不影响我们的仿真)
　

奇异值分解

  我们现在已知的数据就是声压矩阵$\mathrm{P}$，需要根据该矩阵提取得到深度模态函数，而深度模态函数的信息就在矩阵$\Phi$，该矩阵的各个列向量是各号简正波深度模态函数在空间上的采样。
　根据奇异值分解(SVD)的知识，一个矩阵可以进行如下方式的分解：
　$$A=U\Lambda V^H$$$$U{U}^H=I,V{V}^H=I$$
　回到分解开来的声压矩阵：
　1.$\Phi$的各个列向量对应的是各号简正波的深度模态函数，故具有正交归一性，满足酉矩阵的定义；
　2.$\Lambda$是一个对角阵；
　3.如果矩阵$\mathrm{R}$也是个酉矩阵，那么就可以直接通过SVD分解得到$\Phi$，$\Lambda$和$\mathrm{R}$。遗憾的是，$\mathrm{R}$不是个酉矩阵，观察它的形式，发现如果把各列的$\sqrt{\frac{1}{r_i}}$去掉，那么只剩下指数函数，各个行向量同样具有正交归一性；
　所以引出接下来的操作，对声压矩阵$\mathrm{P}$的各列进行能量归一化，即不同的列乘上一个不同的系数，使得各列的总能量(每一列的声压幅值的平方和)相同，满足：
$$\sum _{i=0}^N{p}^{{}^{\prime }}(z_i,r_j)(p^{{}^{\prime }}(z_i,r_j){)}^{\ast }=1$$
　经过这个操作，那么$\mathrm{R}$就变成了一个酉矩阵，此时我们就可以通过SVD的方法来分解得到$\Phi$，$\Lambda$和$\mathrm{R}$了。matlab里的指令为svd()。至此，我们已经提取出了深度模态函数的信息，即矩阵$\Phi$。

特征值分解

  如果声源在水平距离上的取值过多的话(即列数过多)，那么SVD将具有一个庞大的工作量，此时如果我们变换一下形式：
$$C=\mathrm{P}{\mathrm{P}}^H=\Phi \Lambda \mathrm{R}{\mathrm{R}}^H{\Lambda }^H{\Phi }^H=\Phi {\Lambda }^2{\Phi }^H$$
　若$\mathrm{P}$是一个$N\times M$的矩阵$(N<M)$，则$C$是一个$N\times N$方阵，故可以进行特征值分解，特征值组成了矩阵${\Lambda }^2$，特征向量组成了矩阵$\Phi$。matlab里的指令为eig()。至此，同样能够提取出深度模态函数的信息，即矩阵$\Phi$。

提取结果：

注意事项

1.水听器的个数$>$简正波的号数，所以提取出来的特征值的个数$>$简正波的号数，我们需要选取最大的那些特征值，这些才对应各个简正波；然后我们需要找出这些特征值对应的特征向量，这些对应的就是各号简正波的深度模态函数在垂直方向上的采样。
2.在作解析形式的深度模态函数与提取出来的深度模态函数的分析时：

• 如果遇到形状类似，但是存在一个正负号差异的情况，此时需要在该特征向量前面添上一个负号，因为在特征值分解时，特征向量之间仅需要满足正交归一性，所以方向可能会与解析形式的深度模态函数不同；
• 如果遇到幅值对不上的情况，注意解析形式的深度模态函数也是满足正交归一性的，其形式为$\sqrt{\frac{2}{H}}sin(k_{zn}z)$，不要忘记了系数；

* 如何对应特征向量与简正波

  首先需要明确的是，大的那些特征值对应的特征向量才对应不同号数的简正波，接近于0的那些特征值对应的特征向量无意义。
　然后需要找到特征向量对应的那号简正波的模态函数，本来我是直接把特征向量和简正波都画出来，然后手动匹配，这样会导致程序的复用率低。通过和同学的讨论以及文献中的思路，设计了一种自动匹配的算法。

机理分析

1. 观察分解得到的$\Lambda$矩阵，是一个由特征值组成的对角阵，观察特征值，和简正波的幅值有着千丝万缕的联系；再观察能量归一化之前的矩阵$\mathrm{R}$，前面乘上的$\frac{1}{\sqrt{N_R}}$为了是保证能量归一化后的$\mathrm{R}$各个行向量具有归一性，为了补偿，$\Lambda$矩阵前面要乘上一个$\sqrt{N_R}$；
   
   
2. $\Lambda$矩阵其他的系数来源于简正波的表达式，这里的$k_i$是$i$号简正波的水平波数；
3. 如果仅仅想要消去矩阵$\mathrm{R}$中的$\frac{1}{\sqrt{r_i}}$，仅需要对每列声压乘上$\sqrt{r_i}$，但由于我们对声压矩阵的各列进行了能量归一化，所以实际上多乘了一个系数，需要把该系数乘到矩阵$\Lambda$中去，系数$c$：$$c=\frac{1}{{\sum }_{i=0}^{N_z}p(z_i,r_j)p^{\ast }(z_i,r_j)\sqrt{r_i}}$$
4. 提取出来的 $\Lambda$矩阵中的特征值=$c\sqrt{\frac{4\pi N_R}{H{k}_{rn}}}sin(k_{zn}{z}_s)$
5. 经过了一系列系数的补偿，特征值=简正波的幅值

算法思路

因为此时特征值与对应简正波幅值的位置还存在差异，所以需要进行位置匹配。

1. $temp$=|特征值向量|-|$i$号简正波幅值|；
2. 找到$temp$向量中最小值出现的位置，该位置的特征值与第$i$个位置的特征值交换位置，该位置的特征向量与第$i$个位置的特征向量交换位置；
3. 因为所有垂直模态在$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