你的位置：首页 > 信息动态 > 新闻中心
信息动态
联系我们

中恒嘉业

地址：北京市海淀区苏州街大河庄苑底4号楼0105号

电话：15910357833

QQ：：823598168

邮件：823598168@qq.com

网站：http://www.zizl.cn

更多
更多

雷达信号处理---信源数估计算法总结（含MATLAB仿真）

2021/12/29 22:04:05

目录

一.引言

二.信源数估计算法

1.信息论方法

1）MDI准则

设置仿真环境：

​

2）AIC准则

设置仿真环境：

​

2.盖尔圆方法

设置仿真环境：

 三.总结

一.引言

        信源数估计在空间谱估计中，是一项十分关键的技术。对于大部分的超分辨算法来说，都需要已知空间中的信源数目。信源数目的过估计与欠估计都会影响算法的性能。以子空间类算法，比如MUSIC算法为例，信源数目估计错误，将导致噪声子空间的能量流入到信号子空间，导致最后角度估计结果出现较大误差，且出现虚警或漏警的可能。故，在本篇文章中，介绍几种常用的信源数估计方法。
 

二.信源数估计算法

1.信息论方法

信息论方法有一个统一的表达式

J(k)=L(k)+P(k)

其中，L(k)为对数似然函数，P(k)为罚函数。通过选择不同的对数似然函数与罚函数就可以得到不同的准则。譬如，有效检测准则（EDC准则）的表达式为

EDC(k)=L(M-k)ln\Lambda (k)+k(2M-k)C(L)

式中，k 为待估计的信源数，L为采样数，M为阵元个数，\Lambda (k)为似然函数，且

\Lambda (k)=\frac{\frac{1}{M-k}\sum_{i=k+1}^{M}\lambda_{i}}{(\prod_{i=k+1}^{M}\lambda_{i})^{\frac{1}{M-k}}}

式中\lambda_{i} 为特征值。

1）MDI准则

      当C(L) 为\frac{ln(L)}{2} 时，可以得到MDL准则。

MDL(k)=L(M-k)ln\Lambda (k)+\frac{1}{2}k(2M-k)lnL

算法思路为：

       Step1     计算阵元接收信号的协方差矩阵 R，并做特征值分解，将特征值按从大到小排序。

       Step2     将特征值以及阵元参数带入式中，若第d个特征值带入时，该式取最小值，则信源个数为k=d-1

补充：

       当 C(L)满足

 \left\{\begin{matrix} \lim_{k\rightarrow \infty }(C(L)/L)=0\\ \lim_{k\rightarrow \infty }(C(L)/lnlnL)=\infty \end{matrix}\right.

时，该准则具有估计一致性。可见，MDL准则具有估计一致性。

功能函数为

function [num]=MDL(R,N,M,D)
%MDL准则估计信源数目
%input R 阵元接收信号的协方差矩阵
%input N 采样数
%input M 阵元数
%input D 1：一维线阵 2：二维方阵
%output num 信源数目
if D==1
    [u,s,v] = svd(R);
    sd = diag(s);
    sd=sort(sd,'descend');%从大到小排列
    a = zeros(1,M);
    for m = 0:M-1
        negv = sd(m+1:M);
        Tsph = mean(negv)/((prod(negv))^(1/(M-m)));
        a(m+1) = N*(M-m)*log(Tsph) + m*(2*M-m)*log(N)/2;
    end
    [y,b] = min(a);
    num = b - 1;
else
    [u,s,v] = svd(R);
    sd = diag(s);
    sd=sort(sd,'descend');%从大到小排列
    a = zeros(1,prod(M)); %prod([M,N])=M*N
    for m = 0:prod(M)-1
        negv = sd(m+1:prod(M));
        Tsph = mean(negv)/((prod(negv))^(1/(prod(M)-m)));
        a(m+1) = N*(prod(M)-m)*log(Tsph) + m*(2*prod(M)-m)*log(N)/2;
    end
    [y,b] = min(a);
    num = b - 1;
end
end

设置仿真环境

两个相干信号入射到11*11的均匀面阵中，阵元间距半波长，其二维到达角度分别为[1.8° 4.6°]、[5.8° 4.6°]，其中，第一个角度为为方位角，第二个角度为俯仰角。方位角为5.8°的信号比方位角为1.8°的信号信噪比低3dB，快拍数32，信噪比从-5dB到20dB以1dB为间隔变化，每个信噪比点做300次蒙特卡洛实验。当估计信源数等于真是信源数时，判定为估计成功。

My=11;%y轴阵元数
Mz=11;%z轴阵元数
N=32;%快拍数
theta=[ 1.8 5.8];%方向角
phi=[ 4.6 5.6 ];%俯仰角
dis=[ 60 60 ];%信号源 
NUM=length(theta);%信号源数目
f0=3e9;%中频
c=3e8;
B=100e6;%带宽
lambda=c/f0;%波长
Te=20e-6;%脉冲宽度
PRI=120e-6;%脉冲周期
Fs=N/Te;%采样频率
N1=round(Fs*PRI);%一个脉冲周期内的采样点数
mu=B/Te;
d=lambda/2;%阵元间距
t1=0:1/Fs:PRI-1/Fs;
%%接收信号
Am=[1 10^(-0.15)];%幅值\功率
Alpha=[1 1 ];%相干系数
GROUP=NUM;
Ay0=zeros(My,NUM);
Az0=zeros(Mz,NUM);
X0=[];
s=zeros(NUM,N);
My0=6;%划分子阵
Mz0=6;%划分子阵
for i=1:NUM
    s(i,:)=Am(i)*Alpha(i)*exp(1i*pi*B/Te*((0:N-1)/Fs-(2*dis(i)/c)).^2);  %信号模型
end
for i=1:NUM
    Ay0(:,i)=exp(-1i*2*pi*(0:My-1)'*d/lambda*cos(phi(i)/180*pi)*sin(theta(i)*pi/180)); %y域空间导向矢量
end
for i=1:NUM
    Az0(:,i)=exp(-1i*2*pi*(0:Mz-1)'*d/lambda*sin(phi(i)/180*pi)); %z域空间导向矢量
end
for i=1:Mz
      X0=[X0;Ay0*diag(Az0(i,:))*s];
end
X1=X0;
snr=-5:1:20;
m=1;
success=zeros(1,length(snr));
tic;
for snr=-5:20
for time=1:300
X0=awgn(X1,snr,'measured');
%空间平滑
tic;
Rf=zeros(My0*Mz0,My0*Mz0);
x=zeros(My0*Mz0,N);
for i=1:My-My0+1
    for j=1:Mz-Mz0+1
        for k=1:Mz0
            x((k-1)*My0+1:(k-1)*My0+My0,:)=X0((k-1)*My+i+(j-1)*My:(k-1)*My+My0+i-1+(j-1)*My,:);
        end
        E=x*x'/N;
        Rf=Rf+E*E'/N;
    end
end
Rf=Rf/(My0-My0+1);
J=fliplr(eye(My0*Mz0,My0*Mz0));
Rb=J*conj(Rf)*J;
R=(Rf+Rb)/2;
num=MDL(R,N,[My0 Mz0],2);
if num==NUM
    success(m)=success(m)+1;
end
end
toc;
m=m+1;
end
success=success/300;

figure;
plot([-5:1:20],success,'-O');
title('不同信噪比下算法检测成功率');
xlabel('信噪比(dB)');
ylabel('检测成功');grid on;

2）AIC准则

当C(L) 为1时，就得到了AIC准则

AIC(k)=L(M-k)ln\Lambda (k)+\frac{1}{2}k(2M-k)

算法思路为：

       Step1     计算阵元接收信号的协方差矩阵R ，并做特征值分解，将特征值按从大到小排序。

       Step2     将特征值以及阵元参数带入式中，若第d个特征值带入时，该式取最小值，则信源个数为k=d-1

AIC准则并不是一致估计，即快拍数较大时，仍不能完全收敛于真实值.

 功能函数为

function [num]=AIC(R,N,M,D)
%AIC准则估计信源数目
%input R 阵元接收信号的协方差矩阵
%input N 采样数
%input M 阵元数
%input D 1：一维线阵 2：二维方阵
%output num 信源数目
if D==1
    [u,s,v] = svd(R);
    sd = diag(s);
    sd=sort(sd,'descend');
    a = zeros(1,M);
    for m = 0:M-1
        negv = sd(m+1:M);
        Tsph = mean(negv)/((prod(negv))^(1/(M-m)));
        a(m+1) = 2*N*(M-m)*log(Tsph) + m*(2*M-m)*2;
    end
    [y,b] = min(a);
    num = b - 1;
else
    [u,s,v] = svd(R);
    sd = diag(s);
    sd=sort(sd,'descend');
    a = zeros(1,prod(M));
    for m = 0:prod(M)-1
        negv = sd(m+1:prod(M));
        Tsph = mean(negv)/((prod(negv))^(1/(prod(M)-m)));
        a(m+1) = N*(prod(M)-m)*log(Tsph) + m*(2*prod(M)-m);
    end
    [y,b] = min(a);
    num = b - 1;
end
end

设置仿真环境

两个相干信号入射到11*11的均匀面阵中，阵元间距半波长，其二维到达角度分别为[1.8° 4.6°]、[5.8° 4.6°]，其中，第一个角度为为方位角，第二个角度为俯仰角。方位角为5.8°的信号比方位角为1.8°的信号信噪比低3dB，快拍数32，信噪比从-5dB到20dB以1dB为间隔变化，每个信噪比点做300次蒙特卡洛实验。当估计信源数等于真是信源数时，判定为估计成功。

My=11;
Mz=11;
N=32;%快拍数
theta=[ 1.8 5.8];%方向角
phi=[ 4.6 5.6 ];%俯仰角
dis=[ 60 60 ];%信号源 
NUM=length(theta);%信号源数目
f0=3e9;%中频
c=3e8;
B=100e6;%带宽
lambda=c/f0;%波长
Te=20e-6;%脉冲宽度
PRI=120e-6;%脉冲周期
Fs=N/Te;%采样频率
N1=round(Fs*PRI);%一个脉冲周期内的采样点数
mu=B/Te;
d=lambda/2;%阵元间距
t1=0:1/Fs:PRI-1/Fs;
%%接收信号
Am=[1 10^(-0.15)];%幅值\功率
Alpha=[1 1 ];%相干系数
GROUP=NUM;
Ay0=zeros(My,NUM);
Az0=zeros(Mz,NUM);
X0=[];
s=zeros(NUM,N);
My0=6;%划分子阵
Mz0=6;%划分子阵
for i=1:NUM
    s(i,:)=Am(i)*Alpha(i)*exp(1i*pi*B/Te*((0:N-1)/Fs-(2*dis(i)/c)).^2);  %信号模型
end
for i=1:NUM
    Ay0(:,i)=exp(-1i*2*pi*(0:My-1)'*d/lambda*cos(phi(i)/180*pi)*sin(theta(i)*pi/180)); %y域空间导向矢量
end
for i=1:NUM
    Az0(:,i)=exp(-1i*2*pi*(0:Mz-1)'*d/lambda*sin(phi(i)/180*pi)); %z域空间导向矢量
end
for i=1:Mz
      X0=[X0;Ay0*diag(Az0(i,:))*s];
end
X1=X0;
snr=-5:1:20;
m=1;
success=zeros(1,length(snr));
tic;
for snr=-5:20
for time=1:300
X0=awgn(X1,snr,'measured');
%空间平滑
tic;
Rf=zeros(My0*Mz0,My0*Mz0);
x=zeros(My0*Mz0,N);
for i=1:My-My0+1
    for j=1:Mz-Mz0+1
        for k=1:Mz0
            x((k-1)*My0+1:(k-1)*My0+My0,:)=X0((k-1)*My+i+(j-1)*My:(k-1)*My+My0+i-1+(j-1)*My,:);
        end
        E=x*x'/N;
        Rf=Rf+E*E'/N;
    end
end
Rf=Rf/(My0-My0+1);
J=fliplr(eye(My0*Mz0,My0*Mz0));
Rb=J*conj(Rf)*J;
R=(Rf+Rb)/2;
num=AIC(R,N,[My0 Mz0],2);
if num==NUM
    success(m)=success(m)+1;
end
end
toc;
m=m+1;
end
success=success/300;

figure;
plot([-5:1:20],success,'-O');
title('不同信噪比下算法检测成功率');
xlabel('信噪比(dB)');
ylabel('检测成功');grid on;

 观察结果不难发现，即使在高信噪比条件下，AIC准则依然无法收敛于真实值，这印证了AIC准则不是一致估计。

2.盖尔圆方法

Step1 求阵列接收信号的协方差矩阵R ,并对R 进行分块

R=\begin{bmatrix} R'& r\\ r^{H} & r_{MM} \end{bmatrix}

其中，R' 为R前M-1行和前M-1列构成的矩阵。对R' 做特征值分解，并将特征值从大到小排序，与其对应的特征向量V也按顺序重排。使用V构造酉矩阵

 T=\begin{bmatrix} V & 0\\ 0^{T} & 1 \end{bmatrix}

利用酉矩阵对协方差矩阵R做酉变换

R_{T}=T^{H}RT=\begin{bmatrix} \gamma _{1} &0 & \cdots & 0 & \rho _{1}\\ 0& \gamma _{2} & \cdots& 0 &\rho _{2} \\ \vdots & \vdots &\cdots & \vdots & \vdots \\ 0 & 0 & \cdots & \gamma _{M-1}&\rho _{M-1} \\ \rho_{1}^{*} & \rho_{2}^{*} & \cdots& \rho_{m-1}^{*} & r_{MM} \end{bmatrix}

Step2      计算下式

 GDE(k)=\rho _{k}-\frac{D(N)}{M-1}\sum_{i=1}^{M}\rho _{i}

其中D(N) 为调整因子，与快拍数负相关。当若第d个特征值带入时，该式出现负值，则信源个数为k=d-1

功能函数为

function [num]=Gerschgorin_disk_estimation (R,N)
%盖尔圆准则估计信源数目
%input R 协方差矩阵
%input N  快拍数
%output num 信源数
    [row_length col_length]=size(R);
    R_new=R(1:row_length-1,1:col_length-1);
    [V D]=eig(R_new);
    D =diag(D).';
    [D,I0] = sort(D);
    D=fliplr(D);
    V=fliplr(V(:,I0));
    U=[V zeros(row_length-1,1);zeros(1,col_length-1) 1];
    S=U'*R*U;
    k=1;
    while 1
        GDE=abs(S(k,col_length))-1/(2*N*(col_length-1))*sum(abs(S(1:row_length-1,col_length)));%调整因子取为1/N
        if GDE<0 ||k>col_length-2
            break;
        end
        k=k+1;
    end
    num=k-1;
end

设置仿真环境

两个相干信号入射到11*11的均匀面阵中，阵元间距半波长，其二维到达角度分别为[1.8° 4.6°]、[5.8° 4.6°]，其中，第一个角度为为方位角，第二个角度为俯仰角。方位角为5.8°的信号比方位角为1.8°的信号信噪比低3dB，快拍数32，信噪比从-5dB到20dB以1dB为间隔变化，每个信噪比点做300次蒙特卡洛实验。当估计信源数等于真是信源数时，判定为估计成功。

 

My=11;
Mz=11;
N=32;%快拍数
theta=[ 1.8 5.8];%方向角
phi=[ 4.6 5.6 ];%俯仰角
dis=[ 60 60 ];%信号源 
NUM=length(theta);%信号源数目
f0=3e9;%中频
c=3e8;
B=100e6;%带宽
lambda=c/f0;%波长
Te=20e-6;%脉冲宽度
PRI=120e-6;%脉冲周期
Fs=N/Te;%采样频率
N1=round(Fs*PRI);%一个脉冲周期内的采样点数
mu=B/Te;
d=lambda/2;%阵元间距
t1=0:1/Fs:PRI-1/Fs;
%%接收信号
Am=[1 10^(-0.15)];%幅值\功率
Alpha=[1 1 ];%相干系数
GROUP=NUM;
Ay0=zeros(My,NUM);
Az0=zeros(Mz,NUM);
X0=[];
s=zeros(NUM,N);
My0=6;%划分子阵
Mz0=6;%划分子阵
for i=1:NUM
    s(i,:)=Am(i)*Alpha(i)*exp(1i*pi*B/Te*((0:N-1)/Fs-(2*dis(i)/c)).^2);  %信号模型
end
for i=1:NUM
    Ay0(:,i)=exp(-1i*2*pi*(0:My-1)'*d/lambda*cos(phi(i)/180*pi)*sin(theta(i)*pi/180)); %y域空间导向矢量
end
for i=1:NUM
    Az0(:,i)=exp(-1i*2*pi*(0:Mz-1)'*d/lambda*sin(phi(i)/180*pi)); %z域空间导向矢量
end
for i=1:Mz
      X0=[X0;Ay0*diag(Az0(i,:))*s];
end
X1=X0;
snr=-5:1:20;
m=1;
success=zeros(1,length(snr));
tic;
for snr=-5:20
for time=1:300
X0=awgn(X1,snr,'measured');
%空间平滑
tic;
Rf=zeros(My0*Mz0,My0*Mz0);
x=zeros(My0*Mz0,N);
for i=1:My-My0+1
    for j=1:Mz-Mz0+1
        for k=1:Mz0
            x((k-1)*My0+1:(k-1)*My0+My0,:)=X0((k-1)*My+i+(j-1)*My:(k-1)*My+My0+i-1+(j-1)*My,:);
        end
        E=x*x'/N;
        Rf=Rf+E*E'/N;
    end
end
Rf=Rf/(My0-My0+1);
J=fliplr(eye(My0*Mz0,My0*Mz0));
Rb=J*conj(Rf)*J;
R=(Rf+Rb)/2;
num=Gerschgorin_disk_estimation (R,N);
if num==NUM
    success(m)=success(m)+1;
end
end
toc;
m=m+1;
end
success=success/300;

figure;
plot([-5:1:20],success,'-O');
title('不同信噪比下算法检测成功率');
xlabel('信噪比(dB)');
ylabel('检测成功');grid on;

 三.总结

        通过比较三种算法的性能可以知道。三种算法各有优势，MDL准则是一致性估计，在高信噪比条件下会有很好的估计效果。而AIC准则在低信噪比条件下，估计性能比MDL准则出色，在高信噪比条件下，由于自身的非一致性估计，无法收敛至真实值。盖尔圆准则相比于信息论准则，没有特征值分解的操作，因此就资源消耗来说，比两种信息论准则少。

        当然，以上仿真仅考虑高斯白噪声的条件，未考虑高斯色噪声以及阵元互耦等因素。在考虑上述因素的条件下，也有许多相应的改进算法被提出，这个在以后的篇章中会提到。敬请期待