水下声源方位估计

qianpeng897
0 阅读15分钟

水下声源方位估计

前言

水下目标定位问题是三维问题，可以将其转化为角度、距离和深度估计问题[1][2]，如下图所示为三维定位的原理图，在柱坐标系下，可以将声源坐标表示为(r,θ,d)(r,\theta,d)。均匀线阵和圆环阵进行波达方向估计（DOA估计）已被较为全面的研究[3]，波束形成的方法有常规波束形成（CBF）、最佳波束形成（MVDR）等[3]。距离和深度的估计可以使用匹配场方法（MFP）[4]，常见的方法有线性匹配场方法（CMFP）和自适应匹配场方法（AMFP），分别对应着DOA估计中的常规波束形成和最佳波束形成。

一、 波束形成原理

水下目标角度估计利用DOA估计，如下图所示，声源位于水平阵θ\theta角方向上，根据声源发射声波到达阵列的相位差，可以实现声源的方位估计。下面建立数学模型来描述这一过程。 浅海水平阵DOA估计原理图

1.1 阵列接收信号模型

首先阵列接收信号模型，通过该模型推导波束形成的物理原理。 在上图中选择海面中一点作为空间参考点，该点满足从该点向海底做垂线，垂足o′o'为阵列的几何中心，第m个阵元的位置可以表达为，

pm=[pxm,pym,pzm]T,m=1,...,M，(1)\boldsymbol{p}_m=\begin{bmatrix}p_{xm},p_{ym},p_{zm}\end{bmatrix}^\mathrm{T},m=1,...,M {，(1)}

将阵列的位置表示为矩阵形式，P=[p1,...,pm,...pM]H\mathcal{\boldsymbol{P}}=[\boldsymbol{p_1},...,\boldsymbol{p_m},...\boldsymbol{p_M}]^\mathrm{H}。由图1.1，可以定义信号传播方向单位向量为，

v(Ω)=−[sinϕcosθ,sinϕsinθ,cosϕ]T,m=1,...,M，(2)v(\Omega)=-[sin\phi cos\theta,sin\phi sin\theta,cos\phi]^\mathrm{T},m=1,...,M {，(2)}

可以得到阵元相对于参考点o′o'的时延为，

τm(Ω)=vT(Ω)pm/c,m=1,...,M，(3)\tau_m(\Omega)=\boldsymbol{v}^\mathrm{T}(\Omega)\boldsymbol{p}_m/\boldsymbol{c},\quad m=1,...,M {，(3)}

其中c为海水中参考声速。 假设参考点o′o'接收信号为s(t)，则阵元接收信号可以表达成，

sm(t)=s[t−τm(Ω)],m=1,...,M，(4)s_m(t)=s[t-\tau_m(\Omega)],m=1,...,M {，(4)}

转换为频域表达

Sm(ω)=S(ω)e−iωτm(Ω)，(5)S_m(\omega)=S(\omega)e^{-i\omega\tau_m(\Omega)} {，(5)}

式中，S(ω)S(\omega)为信号的频谱。定义波数向量

k(Ω)≜(ω/c)v(Ω)=kv(Ω)，(6)\boldsymbol{k}(\Omega)\triangleq(\omega/c)\boldsymbol{v}(\Omega)=k\boldsymbol{v}(\Omega) {，(6)}

式中，k=ω/c=2π/λk=\omega/c=2\pi/\lambda为波束，c为海洋参考声速，λ\lambda为信号波长。 于是，式(5)可以写成

Sm(ω)=S(ω)e−ikTpm，(7)S_m(\omega)=S(\omega)e^{-i\boldsymbol{k}^\mathrm{T}\boldsymbol{p}_m} {，(7)}

阵元接收信号可以表达成M×1列向量的形式，

xs(t)=[s1(t),...,sm(t),...,sM(t)]T，(8)\boldsymbol{x}_s(t)=[s_1(t),...,s_m(t),...,s_M(t)]^\mathrm{T} {，(8)}

根据式(7)和(8)。可以得到，

Xs(ω)=[S1(ω)...Sm(ω)...SM(ω)]=[exp⁡(−ikTp1)...exp⁡(−ikTpm)...exp⁡(−ikTpM)]S(ω)，(9)\boldsymbol{X}_s(\omega)=\begin{bmatrix}S_1(\omega)\\...\\S_m(\omega)\\...\\S_M(\omega)\end{bmatrix}=\begin{bmatrix}\exp(-i\boldsymbol{k}^\mathrm{T}\boldsymbol{p}_1)\\...\\\exp(-i\boldsymbol{k}^\mathrm{T}\boldsymbol{p}_m)\\...\\\exp(-i\boldsymbol{k}^\mathrm{T}\boldsymbol{p}_M)\end{bmatrix}S(\omega) {，(9)}

定义p(k,P)≜[p1(k,p1),...,pm(k,pm),...,pM(k,pM)]T≜[exp⁡(−ikTp1),...,exp⁡(−ikTpm),...exp⁡(−ikTpM)]T\boldsymbol{p}(\boldsymbol{k},\boldsymbol{P})\triangleq[p_1(\boldsymbol{k},\boldsymbol{p}_1),...,p_m(\boldsymbol{k},\boldsymbol{p}_m),...,p_M(\boldsymbol{k},\boldsymbol{p}_M)]^\mathrm{T}\triangleq[\exp(-i\boldsymbol{k}^{\mathrm{T}}\boldsymbol{p}_{1}),...,\exp(-i\boldsymbol{k}^{\mathrm{T}}\boldsymbol{p}_{m}),...\exp(-i\boldsymbol{k}^{\mathrm{T}}\boldsymbol{p}_{M})]^{\mathrm{T}}为方向导向向量。根据式(9)，并加入干扰和噪声，得到第m号阵元接收信号为

Xm(ω)=βS0(ω)e−iωτm(Ω0)+∑d=1DSd(ω)e−iωτm(Ωd)+Nm(ω)，(10)X_m(\omega)=\beta S_0(\omega)e^{-i\omega\tau_m(\Omega_0)}+\sum_{d=1}^DS_d(\omega)e^{-i\omega\tau_m(\Omega_d)}+N_m(\omega) {，(10)}

式中，β表示目标信号是否在接收信号中，β=1，信号存在；β=0，信号不存在。将式(10)表示为M×1维向量的形式，

X(ω)=βp(k0)S0(ω)+∑d=1Dp(kd)Sd(ω)+N(ω)，(11)\boldsymbol{X}(\omega)=\beta \boldsymbol{p}(\boldsymbol{k_0})S_0(\omega)+\sum_{d=1}^D\boldsymbol{p}(\boldsymbol{k_d})S_d(\omega)+\boldsymbol{N}(\omega) {，(11)}

式中，p(kd)(d=0,1,…,D)\boldsymbol{p}(\boldsymbol{k_d})(d=0,1,…,D)表示期望信号与干扰对应的方向响应向量。 假设信号、干扰、噪声互不相关，得到互谱矩阵为

Sx(ω)=E[X(ω)XH(ω)]=βSxs(ω)+Sxi(ω)+Sn(ω)=βp(k0)Ss0(ω)pH(k0)+∑d=1Dp(kd)Ssd(ω)pH(kd)+Sn(ω)ρn(ω)，(12)\boldsymbol{S}_\mathrm{x}(\omega)=E[\boldsymbol{X}(\omega)\boldsymbol{X}^\mathrm{H}(\omega)]=\beta\boldsymbol{S}_\mathrm{xs}(\omega)+\boldsymbol{S}_\mathrm{xi}(\omega)+\boldsymbol{S}_\mathrm{n}(\omega)\\=\beta\boldsymbol{p}(\boldsymbol{k}_0)S_{\mathrm{s}0}(\omega)\boldsymbol{p}^\mathrm{H}(\boldsymbol{k}_0)+\sum_{d=1}^D\boldsymbol{p}(\boldsymbol{k}_d)S_{\mathrm{s}d}(\omega)\boldsymbol{p}^\mathrm{H}(\boldsymbol{k}_d)+S_\mathrm{n}(\omega)\boldsymbol{\rho}_\mathrm{n}(\omega) {，(12)}

1.2 波束形成原理

波束形成方法的原理是，将阵列信号进行时延补偿，使信号同相相加，而噪声和干扰非同相相加，从而提高信噪比。 假设第m通道的传输函数为wm(ω)w_m(\omega)， w(ω)=[w1(ω),...wm(ω),...wM(ω)]T\boldsymbol{w}(\omega)=[w_1(\omega),...w_m(\omega),...w_M(\omega)]^\mathrm{T}（也称之为加权向量）频域波束形成的波束输出为，

Y(ω)=∑m=1Mwm(ω)Xm(ω)=wT(ω)X(ω)，(13)\boldsymbol{Y}(\omega)=\sum_{m=1}^Mw_m(\omega)X_m(\omega)=\boldsymbol{w}^\mathrm{T}(\omega)\boldsymbol{X}(\omega) {，(13)}

若需要得到输出功率，可将式(13)改为

σy2=E[Y(ω)YH(ω)]=wH(ω)Sx(ω)w(ω)，(14)\sigma_y^2=E[\boldsymbol{Y}(\omega)\boldsymbol{Y}^\mathrm{H}(\omega)]=\boldsymbol{w}^\mathrm{H}(\omega)\boldsymbol{S}_\mathrm{x}(\omega)\boldsymbol{w}(\omega) {，(14)}

常规波束形成器的加权向量为，

wc=p‾s/M，(15)\boldsymbol{w}_c=\overline{p}_s/M {，(15)}

式中，p‾s\overline{p}_s为信号的方向导向向量。 最佳波束形成器的加权向量为

wMVDR=Rx−1p‾sp‾sHRx−1p‾s，(16)\boldsymbol{w}_\mathrm{MVDR}=\frac{\boldsymbol{R}_\mathrm{x}^{-1}\boldsymbol{\overline{p}}_s}{\boldsymbol{\overline{p}}_s^\mathrm{H}\boldsymbol{R}_\mathrm{x}^{-1}\boldsymbol{\overline{p}}_s} {，(16)}

式中，分母 p‾sHRx−1p‾s\overline{\boldsymbol{p}}_{s}^{{\mathrm{H}}}R_{{\mathrm{x}}}^{-1}\overline{\boldsymbol{p}}_{s} 是一个标量系数。

1.3 波束响应和方位谱

波束响应是指某个方向的单位功率平面波在波束形成后的能量响应，用来考察波束形成器对不同方向来波信号的响应特性。省略频率宗量，波束响应公式为 B(Ω)=wHp(Ω),Ω∈Θ，(17)\boldsymbol{B}(\Omega)=\boldsymbol{w}^\mathrm{H}\boldsymbol{p}(\Omega),\quad\Omega\in\Theta {，(17)} 式中，Θ\Theta表示信号所有可能到达的方向集合（观察视区），对于三维观察空间，Θ\Theta为空间角；对于平面问题，此时ϕ=90°\phi=90°，Θ\Theta为xo'y平面角，θ∈[−180∘,180∘]\theta\in[-180^\circ,180^\circ]。 由式(17)可以看出，波束响应的加权向量不含自变量Ω\Omega，而方向响应向量中 含有自变量Ω，这一点要与方位谱公式区分开来。 方位谱是指在观察方向波束形成器的输出功率，其计算公式为 P(Ω)=wH(Ω)Rxw(Ω)，(18)\boldsymbol{P}(\Omega)=\boldsymbol{w}^{\mathrm{H}}(\Omega)\boldsymbol{R}_{\mathrm{x}}\boldsymbol{w}(\Omega) {，(18)} 由(18)可知，方位谱计算公式中加权向量含自变量ΩΩ，而RxR_x为阵元接收信号信号的协方差矩阵，不含自变量ΩΩ。 宽带波束形成方位谱的原理为窄带方位谱的加权求和 PB(Ω)=∑j=1NWjPN(fj,Ω)，(19)\boldsymbol{P}_B(\Omega)=\sum_{j=1}^NW_j\boldsymbol{P}_N(f_j,\Omega) {，(19)} 式中，PN(fj,Ω)P_N (f_j,Ω)为对应频率为f_j的窄带方位谱，PB(Ω)P_B (Ω)为宽带方位谱，WjW_j对应为加权系数，该式的求和方法为相干求和，若改为非相干求和，有 PB(Ω)=∑j=1NWj∣PN(fj,Ω)∣，(20)\boldsymbol{P}_B(\Omega)=\sum_{j=1}^NW_j|\boldsymbol{P}_N(f_j,\Omega)| {，(20)}

二、仿真

项目地址声源方位角估计

2.1. 水平均匀线阵波达方向估计基础仿真

考虑一个11元标准线列阵，坐标系采用笛卡尔坐标系，水平阵的位置从(0,-40m,42m)到(0,40m,42m)沿y轴方向均匀分布11个阵元，相邻阵元间隔d=8m。观察视区为xo'y平面，即ϕ=90°，Θ∈[-90°,90°]。对于水平均匀线阵，若无说明，阵列的布置与此阵列保持一致。 得到阵列坐标如下图所示。 假设期望信号方向为θ0=10°θ_0=10°，信号频率为50Hz，设置声速c=1500m/s，则平面波波长为λ=c/f=30mλ=c/f=30m，阵元间距d=8m<λ/2=15md=8m<\lambda/2=15m，满足空间采样定理，因此波束形成结果不会出现副极大。 得到波束图如下图所示。 如图所示，用对数表示的最高旁瓣值与期望方向主瓣值之差称为旁瓣级，用SLL表示。旁瓣级束宽（BWSLBW_{SL}），即波束主瓣功率下降到-3dB时的两方向间夹角。 得到该线阵常规波束旁瓣级SLL≈-13dB，旁瓣级束宽BWSL≈19°BW_{SL}≈19°。

2.1.1 常规波束图、常规波束扫描方位谱

考虑一个11元水平均匀线阵，假设一个功率为1的50Hz单频信号从10°方向入射到基阵，不考虑噪声，得到常规波束形成的常规波束图、常规波束扫描方位谱。 下面分别介绍各子图的意义。 (a)波束图的作用是指定目标方位，观测不同方向入射的信号对输出功率的影响。由图(a)可以看出在θ=10°时波束图出现峰值，说明对于该波束形成器来说10°方向的信号产生的响应最大。 (b)方位谱图被直接用来输出DOA估计结果。由图(b)可以看出，在θ=10°方位谱出现峰值，DOA估计的结果为10°，还可以看出CBF的旁瓣级较高，“能量泄露”较为严重。 (c )加权范数图，加权向量范数的大小反映波束形成器的灵敏度和稳健性，由图(c )可以看出，CBF方法在不同角度的加权向量范数均为-10dB，说明常规波束形成灵敏度低，稳健性好。 (d)阵增益图的作用是，研究波束形成的前后的信噪比增益与波束形成方向的关系。由图(d)可以看出，当波束对准信号方向时，即θ0=10°θ_0=10°，阵增益达到最大值10dB，当波束形成方向与目标方位由差别时，阵增益下降。

波达方向结果

主程序CBF_LA_NB.m代码如下（示例）：

% 清除命令窗口、工作区和关闭所有图形窗口
clc; clear all; close all;

% 仿真环境参数设置
c = 1500; % 声速
f_s = 50; % 采样频率
M = 11; % 阵元数
d_r = 8; % 阵元间隔

% 设置信号的俯仰角和入射角
phi = 90; % 信号俯仰角
theta = 10; % 信号入射角

% 设置信号幅度
sw = 1;

% 生成波束扫描角度范围
theta_B = [-90:0.1:90].';

% 创建阵列位置向量，P_r 的三个维度分别对应 x、y、z 坐标
P_r = [zeros(1,M); ((1:M) - (M +1)/2)*d_r; zeros(1,M)];

% 绘制阵列图像
figure;
% beaplot 函数用于绘制图形，这里传入的参数依次为 x 坐标、y 坐标、标记类型、标记大小、标题、x 轴标签、y 轴标签、是否显示图例、图例名称
beaplot(zeros(1,M), P_r(2,:),'o',1,'阵列流形','x(m)','y(m)',0,0);

% 计算波数
ksc = 2* pi* f_s/ c;

% 计算入射方向向量
pve = beamscp(ksc, P_r, theta, phi);

% 计算接收信号向量
X_s = pve* sw;

% 计算接收信号的自相关矩阵
Rx_s = X_s * X_s';

% 波束图

% 计算期望方向的权向量
w_c = pve/ M;

% 计算不同角度下的方向向量
pve_B = beamscp(ksc, P_r, theta_B, phi);

% 计算波束形成输出
B = w_c'* pve_B;

% 将波束形成输出转换为分贝形式
B_dB = 20* log10(abs(B));

% 绘制波束图
figure; 
beaplot(theta_B, B_dB, '-k', 1.5,'波束图','θ/(°)','20lgB(θ)/(dB)',0,0);
axis([-90 90 -70 0]);
grid off;

% 方位谱

% 计算方位谱的权向量
w_AziSpe = pve_B/ M;

% 初始化方位谱功率
P_AziSpe = zeros(length(theta_B),1);

% 采用 for 循环计算不同角度下的方位谱功率
for n =1 :length(theta_B)
    % 计算权向量与接收信号自相关矩阵的乘积，再与权向量转置相乘得到功率
    P_AziSpe(n,1) = w_AziSpe(:,n)'* Rx_s *w_AziSpe(:,n);
end

% 将方位谱功率转换为分贝形式
P_AziSpe_dB = 10*log10(abs(P_AziSpe));

% 绘制多个子图
figure;
set(gcf,'position',[100 100 800 600]);
subplot(221)
beaplot(theta_B, B_dB, '-k', 1.5,'(a)波束图','θ/(°)','20lgB(θ)/(dB)',0,0);
axis([-90 90 -70 5]);
grid on;

subplot(222)
beaplot(theta_B, P_AziSpe_dB, '-k', 1.5,'(b)纯信号时波束扫描方位谱'...
   ,'θ_0/(°)','10lgP(θ_0)/(dB)',0,0);
axis([-90 90 -70 5]);
grid on;

% % 设置信噪比并计算加噪声后的接收信号自相关矩阵
% SNR = 30;
% sigma2_s = 10^(SNR/10);
% Rx = sigma2_s * Rx_s + eye(M);

% % 初始化加噪声后的方位谱功率
% P_AziSpe2 = zeros(length(theta_B),1);
% for n =1 :length(theta_B)
%     P_AziSpe2(n,1) = w_AziSpe(:,n)'* Rx *w_AziSpe(:,n);
% end
% P_AziSpe_dB2 = 10*log10(abs(P_AziSpe2));

% 计算加权向量范数
w2 = sum(abs(w_AziSpe.^2),1);
w2_db = 10*log10(w2);

subplot(223)
beaplot(theta_B, w2_db, '-k', 1.5,'(c)加权向量范数'...
   ,'θ_0/(°)','10lgP(θ_0)/(dB)',0,0);
axis([-90 90 -15 30]);
grid on;

% 计算不同方向波束阵增益
Gain_w = zeros(length(theta_B),1);
for n =1 :length(theta_B)
    Gain_w(n,1) =abs(w_AziSpe(:,n)'* pve)^2/abs( w_AziSpe(:,n)'* w_AziSpe(:,n));
end
Gain_w_dB = 10*log10(abs(Gain_w));

subplot(224)
beaplot(theta_B, Gain_w_dB, '-k', 1.5,'(d)不同方向波束阵增益'...
   ,'θ_0/(°)','10lgGc(θ_0)/(dB)',0,0);
grid on;
axis([-90 90 -60 15]);

依赖方向响应向量计算函数beamscp.m代码如下（示例）：

function pve_B = beamscp(ksc, P_r, theta, phi)
% 函数说明：该函数输入为波数、阵元位置、入射方向角，输出方向响应向量

% 根据俯仰角 phi 和入射角 theta 计算单位方向向量 v_s
% sind 和 cosd 函数分别用于计算角度的正弦和余弦值（角度以度为单位）
% v_s 的三个分量分别对应 x、y、z 方向
v_s = - [sind(phi)* cosd(theta), sind(phi)* sind(theta), ones(length(theta),1)*cosd(phi)].';

% 根据波数 ksc 和单位方向向量 v_s 计算波矢 kve
kve = ksc * v_s;

% 计算方向响应向量 pve_B
% exp 函数用于计算指数，-1j * kve.' * P_r 表示对每个阵元位置与波矢的乘积取复指数
% 最后转置并输出
pve_B = exp(-1j * kve.' * P_r).';
end

依赖画图函数beaplot.m代码如下（示例）：

function h = beaplot(x,y,ltype,lwidth,Plottitle,xlab,ylab,reverse,xti)
% 函数说明：该函数用来修改 plot 函数的输出使其更美观。
% 注意事项：此函数未设计缺省值，泛化能力不强。如果需要修改，建议编写 default 语句以防止影响以前的程序运行。

% 设置图形背景颜色为白色
set(0,'defaultfigurecolor','w');

% 根据输出参数个数进行不同操作
if nargout == 0
    % 如果没有输出参数，直接绘制图形
    plot(x,y,ltype,'linewidth',lwidth);
elseif nargout == 1
    % 如果有一个输出参数，将绘制的图形句柄赋值给输出变量 h
    h = plot(x,y,ltype,'linewidth',lwidth);
end

% 设置图形标题
title(Plottitle,'fontsize',14,'fontname','黑体');

% 设置坐标轴字体为黑体
set(gca,'fontname','黑体');

% 设置 x 轴标签
xlabel(xlab);

% 设置 y 轴标签
ylabel(ylab);

% 设置 x 轴和 y 轴标签的字体大小和字体为黑体
set(get(gca,'xlabel'),'fontsize',12,'fontname','黑体');
set(get(gca,'ylabel'),'fontsize',12,'fontname','黑体');

% 根据参数 xti 的值设置 x 轴刻度
if xti==1
    set(gca,'xtick',x);
end

% 根据参数 reverse 的值进行坐标轴反转操作
if reverse
    axis ij;
end
end

2.1.2 MVDR波束图、MVDR波束扫描方位谱

考虑一个11元标准线列阵，假设基阵接收噪声是功率为0dB的高斯白噪声，即R_n=I。一个信噪比SNR=30dB的50Hz单频信号从θ_0=10°方向入射到基阵，得到最佳波束形成如下图所示。 (a)波束图，当θ_0=-10°，即期望方向为-10°，此时的波束图如图(a)黑色实线所示，可以看出波束响应在信号方向（θ=10°）有很明显的凹陷，这是由于当信号方向与期望方向不一致，波束形成时将信号当作干扰，抑制波束形成后功率大小。 (b)方位谱，可以看出MVDR方位谱呈现“尖峰状”，相对于CBF来说，MVDR有更好的分辨力。 (c)加权范数图，加权向量范数的大小反映波束形成器的稳健性，由图(c)可以看出，MVDR的加权向量范数在接近目标方位时增大，当θ=10°时突然减小至-10dB，此时MVDR退化为常规波束形成。MVDR在目标方位附近时加权向量范数较高，对应灵敏度高，失配对MVDR的影响较大，说明MVDR的稳健性较差[3]。 (d)阵增益图，由(d)图可以看出当θ_0=10°时，阵增益达到最大值10dB，在其他方向，阵增益非常小，说明MVDR相对CBF来说有更好的抗干扰能力。

MVDR波束图

% 清除命令窗口、工作区和关闭所有图形窗口
clc; clear all; close all;

% 设置声速、采样频率、阵元数和阵元间隔
c = 1500; % 声速
f_s = 50; % 采样频率
M = 11; % 阵元数
d_r = 8; % 阵元间隔

% 信号仿真

% 设置信号的俯仰角和入射角
phi = 90; % 信号俯仰角
theta = 10; % 信号入射角
% 设置信号幅度
sw = 1;
% 设置信噪比
SNR = 30;

% 生成波束扫描角度范围
theta_B = [-90:0.1:90].';


% 创建阵列位置向量，P_r 的三个维度分别对应 x、y、z 坐标
P_r = [zeros(1,M); ((1:M) - (M +1)/2)*d_r; zeros(1,M)];

% 绘制阵列图像
figure;
beaplot(zeros(1,M), P_r(2,:),'o',1,'阵列流形','x(m)','y(m)',0,0);

% 计算波数
ksc = 2* pi* f_s/ c;

% 计算入射方向向量
pve = beamscp(ksc, P_r, theta, phi);

% 计算接收信号向量
X_s = pve* sw;

% 计算接收信号的自相关矩阵
Rx_s = X_s * X_s';

% 计算噪声功率
sigma2_s = 10^(SNR/10);

% 计算加噪声后的接收信号自相关矩阵
Rx = sigma2_s * Rx_s + eye(M);

% 波束图

% 计算 MVDR 波束形成器的权向量（期望方向为 theta）
w_MVDR_1 = Rx\ pve/(pve'/ Rx* pve);

% 计算不同角度下的方向向量
pve_B = beamscp(ksc, P_r, theta_B, phi);

% 计算波束形成输出
B_1 = w_MVDR_1'* pve_B;

% 将波束形成输出转换为分贝形式
B_dB_1 = 20* log10(abs(B_1));

% 绘制波束图
figure;
set(gcf,'position',[100 100 800 600]);
subplot(221)
beaplot(theta_B, B_dB_1, '--r', 1.5,'(a)波束图','θ/(°)','20lgB(θ)/(dB)',0,0);
axis([-90 90 -70 0]);
hold on;

% 计算另一个入射角的方向向量和 MVDR 权向量
pve2 = beamscp(ksc, P_r, -10, phi);
w_MVDR_2 = Rx\ pve2/(pve2'/ Rx* pve2);
B_2 = w_MVDR_2'* pve_B;
B_dB_2 = 20* log10(abs(B_2));
beaplot(theta_B, B_dB_2, '-k', 1.5,'(a)波束图','θ/(°)','20lgB(θ)/(dB)',0,0);
legend('θ_0=10°','θ_0= - 10°');
grid on;

% 方位谱

% 初始化方位谱功率、加权向量范数和波束阵增益
P_AziSpe = zeros(length(theta_B),1);
w_AziSpe2 = zeros(length(theta_B),1);
Gain_w = zeros(length(theta_B),1);

% 采用 for 循环计算不同角度下的方位谱功率、加权向量范数和波束阵增益
for n =1 :length(theta_B)
    w_AziSpe = Rx\ pve_B(:,n)/(pve_B(:,n)'/ Rx* pve_B(:,n));
    w_AziSpe2(n,1) = w_AziSpe'*w_AziSpe;
    P_AziSpe(n,1) = w_AziSpe'* Rx *w_AziSpe;
    Gain_w(n,1) =abs(w_AziSpe'* pve)^2/abs( w_AziSpe'* w_AziSpe);
end

% 将方位谱功率转换为分贝形式
P_AziSpe_dB = 10*log10(abs(P_AziSpe));

subplot(222)
beaplot(theta_B, P_AziSpe_dB, '-k', 1.5,'(b)方位谱','θ_0/(°)','10lgP(θ_0)/(dB)',0,0);
axis([-90 90 -15 40]);
grid on;

% 将加权向量范数转换为分贝形式
w_AziSpe2_dB = 10*log10(abs(w_AziSpe2));

subplot(223)
beaplot(theta_B, w_AziSpe2_dB, '-k', 1.5,'(c)加权向量范数'...
  ,'θ_0/(°)','10lg||w(θ_0)||^2/(dB)',0,0);
axis([-90 90 -15 30]);
grid on;

% 将波束阵增益转换为分贝形式
Gain_w_dB = 10*log10(abs(Gain_w));

subplot(224)
beaplot(theta_B, Gain_w_dB, '-k', 1.5,'(d)不同方向波束阵增益'...
  ,'θ_0/(°)','10lgGc(θ_0)/(dB)',0,0);
grid on;
axis([-90 90 -150 15]);

总结

以上就是今天要讲的内容，本文介绍了波束形成算法的原理和仿真，讨论了CBF和MVDR的原理和仿真结果。

avatar
文章
阅读
粉丝