存在站址误差的DOA定位CRLB

站址误差下的DOA观测模型

在实际过程中，观测站可能处于运动状态，从而获得的观测站位置是带有误差的，将带有误差的观测站位置建模为：${\mathbf{s}}_{\mathbf{i}}={\mathbf{s}}_{\mathbf{i}}^{\mathbf{o}}+{\mathbf{n}}_{{\mathbf{s}}_{\mathbf{i}}}=[s_{i,x}^o,s_{i,y}^o,s_{i,z}^o{]}^T+[\delta s_{i,x},\delta s_{i,y},\delta s_{i,z}{]}^T$。其中${\mathbf{s}}_{\mathbf{i}}^{\mathbf{o}}$表示观测站真实位置向量，${\mathbf{n}}_{{\mathbf{s}}_{\mathbf{i}}}$表示随机误差。则将所有观测站向量组合得到：
$$\mathbf{s}={\mathbf{s}}^o+{\mathbf{n}}_{\mathbf{s}}=[({\mathbf{s}}_1^{\mathbf{o}}{)}^T,({\mathbf{s}}_2^{\mathbf{o}}{)}^T,...,({\mathbf{s}}_{\mathbf{M}}^{\mathbf{o}}{)}^T{]}^T+[{\mathbf{n}}_{{\mathbf{s}}_1}^{\mathbf{T}},{\mathbf{n}}_{{\mathbf{s}}_2}^{\mathbf{T}},...,{\mathbf{n}}_{{\mathbf{s}}_{\mathbf{M}}}^{\mathbf{T}}{]}^T$$
${\mathbf{n}}_{\mathbf{s}}$的协方差矩阵表示为${\mathbf{Q}}_{\mathbf{s}}$。

站址误差下的DOA定位的CRLB

存在站址误差时的未知向量为${\mathbf{\kappa }}^o=[({\mathbf{u}}^o{)}^T,({\mathbf{s}}^o{)}^T{]}^T$，测量向量为角度测量值，记为$\mathbf{z}$。则此时对数似然函数表示为：
$$ln(p(\mathbf{z}|{\mathbf{\kappa }}^o))=\kappa -(1/2)(\mathbf{z}-{\mathbf{z}}^o{)}^T{Q}^{-1}(\mathbf{z}-{\mathbf{z}}^o)-(1/2)(\mathbf{s}-{\mathbf{s}}^o{)}^T{Q}_{\mathbf{s}}^{-1}(\mathbf{s}-{\mathbf{s}}^o)$$
Fisher信息矩阵表示为：
$$\mathbf{F}=\begin{array}{c}\left[\begin{array}{cc}{\mathbf{F}}_1& {\mathbf{F}}_2\\ {\mathbf{F}}_2^T& {\mathbf{F}}_3\end{array}\right]\end{array}$$
其中：
$$\{\begin{array}{l}{\mathbf{F}}_1=\left(\left(\frac{\partial {\mathbf{z}}^{\mathbf{o}}}{\partial ({\mathbf{u}}^{\mathbf{o}}{)}^T}\right)Q^{-1}{\left(\frac{\partial {\mathbf{z}}^{\mathbf{o}}}{\partial ({\mathbf{u}}^{\mathbf{o}}{)}^T}\right)}^T\right)\\ {\mathbf{F}}_2=\left(\left(\frac{\partial {\mathbf{z}}^{\mathbf{o}}}{\partial ({\mathbf{u}}^{\mathbf{o}}{)}^T}\right)Q^{-1}{\left(\frac{\partial {\mathbf{z}}^{\mathbf{o}}}{\partial ({\mathbf{s}}^{\mathbf{o}}{)}^T}\right)}^T\right)\\ {\mathbf{F}}_3=\left(\left(\frac{\partial {\mathbf{z}}^{\mathbf{o}}}{\partial ({\mathbf{s}}^{\mathbf{o}}{)}^T}\right)Q^{-1}{\left(\frac{\partial {\mathbf{z}}^{\mathbf{o}}}{\partial ({\mathbf{s}}^{\mathbf{o}}{)}^T}\right)}^T\right)+{\mathbf{Q}}_{\mathbf{s}}^{-1}\end{array}$$

存在站址误差下的DOA定位CRLB表示为：
$$\mathbf{CRLB}=\begin{array}{c}\left[\begin{array}{cc}{\mathbf{X}}_1& {\mathbf{X}}_2\\ {\mathbf{X}}_2^T& {\mathbf{X}}_3\end{array}\right]\end{array}$$
其中：
$$\mathbf{CRLB}({\mathbf{u}}^o)={\mathbf{X}}_1={\left({\mathbf{F}}_1-{\mathbf{F}}_2{\mathbf{F}}_3^{-1}{\mathbf{F}}_2^T\right)}^{-1}={\mathbf{F}}_1^{-1}+{\mathbf{F}}_1^{-1}{\mathbf{F}}_2{\left({\mathbf{F}}_3-{\mathbf{F}}_2^T{\mathbf{F}}_1^{-1}{\mathbf{F}}_2\right)}^{-1}{\mathbf{F}}_2^T{\mathbf{F}}_1^{-1}$$

$$\mathbf{CRLB}({\mathbf{s}}^o)={\mathbf{X}}_3={\left({\mathbf{F}}_3-{\mathbf{F}}_2{\mathbf{F}}_1^{-1}{\mathbf{F}}_2^T\right)}^{-1}={\mathbf{F}}_3^{-1}+{\mathbf{F}}_3^{-1}{\mathbf{F}}_2{\left({\mathbf{F}}_1-{\mathbf{F}}_2^T{\mathbf{F}}_3^{-1}{\mathbf{F}}_2\right)}^{-1}{\mathbf{F}}_2^T{\mathbf{F}}_3^{-1}$$

matlab仿真

%CRLB随角度误差变化结果
clc;
close all;
clear all;
% 测向站数量
M=6;
% 目标数量
N=1;
%观测站位置
s1=[1500 -1900 1400].';
s2=[-1600 2000 1700].';
s3=[1800 -1500 -1600].';
s4=[-1200 1400 -2000].';
s5=[1700 1600 2000].';
s6=[-1800 -1400 -1800].';
for i=1:1:M
    str_si=['s',mat2str(i)];
    si_value=eval(str_si);
    s_true(3*(i-1)+1:3*(i-1)+3,1)=si_value;
end
% 目标位置
u1=[4000 4000 3000].';
for i=1:1:N
    str_ui=['u',mat2str(i)];
    ui_value=eval(str_ui);
    u_true(3*(i-1)+1:3*(i-1)+3,1)=ui_value;
end
z_error_more(1,:)=10.^[0:0.1:1.5];
z_error_more(2,:)=10.^[zeros(length(z_error_more),1)];
deta_more=[0.001,0.001,1,10];
% 站址误差功率不一样
power_s_error=[1 1 1 3 3 3 5 5 5 10 10 10 5 5 5 3 3 3].';
%% 存在站址误差
%% CRLB
for i=1:1:M
    eval(['s',mat2str(i),'(:,1)',' = ','s_true(3*(i-1)+1:3*(i-1)+3,1)',';']);
end
for i=1:1:N
    eval(['u',mat2str(i),'(:,1)',' = ','u_true(3*(i-1)+1:3*(i-1)+3,1)',';']);
end
%%     对目标位置的导数
% 方位角、仰角对目标位置的导数
diff_theta_ui=zeros(M,3,N);
diff_beta_ui=zeros(M,3,N);
%方位角、仰角对目标观测站的导数
diff_theta_si=zeros(M,3*M,N);
diff_beta_si=zeros(M,3*M,N);
%%
for j=1:1:N
    for i=1:1:M
        str_si=['s',mat2str(i)];
        si_value=eval(str_si);
        str_ui=['u',mat2str(j)];
        ui_value=eval(str_ui);
        
        ui_si=ui_value-si_value;
        ai2=sqrt(ui_si(1)^2+ui_si(2)^2);
        %     方位角对目标位置的导数
        diff_theta_ui_1=-ui_si(2)/ai2^2;
        diff_theta_ui_2=ui_si(1)/ai2^2;
        diff_theta_ui(i,:,j)=[diff_theta_ui_1,diff_theta_ui_2,0];
        %     仰角对目标位置的导数
        diff_beta_ui_1=-ui_si(1)*ui_si(3)/ai2/(osjl(ui_value,si_value))^2;
        diff_beta_ui_2=-ui_si(2)*ui_si(3)/ai2/(osjl(ui_value,si_value))^2;
        diff_beta_ui_3=ai2/(osjl(ui_value,si_value))^2;
        diff_beta_ui(i,:,j)=[diff_beta_ui_1,diff_beta_ui_2,diff_beta_ui_3];
        
        % 方位角、仰角对目标观测站位置的导数
        diff_theta_si(i,3*(i-1)+1:3*(i-1)+3,j)=-diff_theta_ui(i,:,j);
        diff_beta_si(i,3*(i-1)+1:3*(i-1)+3,j)=-diff_beta_ui(i,:,j);
    end
end

F_u1=[diff_theta_ui(:,:,1);diff_beta_ui(:,:,1)];
F_u=blkdiag(F_u1);

F_s1=[diff_theta_si(:,:,1);diff_beta_si(:,:,1)];
F_s=[F_s1.'].';

for big_big_loop=1:1:length(z_error_more)
    %方位角
    deta_theta=deta_more(1)*z_error_more(1,big_big_loop);
    R_theta=deta_theta^2*eye((M),(M));
    %仰角
    deta_beta=deta_more(2)*z_error_more(1,big_big_loop);
    R_beta=deta_beta^2*eye((M),(M));
    
    Rn_a1=blkdiag(R_theta,R_beta);
    Rn_a=blkdiag(Rn_a1);
    %站址误差
    deta_s=deta_more(4)*z_error_more(2,big_big_loop);
    power_s_error=[1 1 1 3 3 3 5 5 5 10 10 10 5 5 5 3 3 3].';
    Rs_a=deta_s^2*diag(power_s_error);
    %     观测站定位的先验均方根误差
    prior_cov_s(big_big_loop,1)=sqrt(trace(Rs_a));
    
    F1=F_u.'*inv(Rn_a)*F_u;
    F2=F_u.'*inv(Rn_a)*F_s;
    F3=F_s.'*inv(Rn_a)*F_s+inv(Rs_a);
    
    CRLB=inv([F1,F2;F2.',F3]);
    a=CRLB(1,1)+CRLB(2,2)+CRLB(3,3);
    
    crb_source1(big_big_loop,1)=sqrt(a);
    CRLB_sensor(big_big_loop,1)=sqrt(trace(CRLB)-(a));
end
CRLB_u_sensor_erros=[crb_source1];
%% 不存在站址误差
for big_big_loop=1:1:length(z_error_more)
    %方位角
    deta_theta=deta_more(1)*z_error_more(1,big_big_loop);
    R_theta=deta_theta^2*eye((M),(M));
    %仰角
    deta_beta=deta_more(2)*z_error_more(1,big_big_loop);
    R_beta=deta_beta^2*eye((M),(M));
    
    Rn_a1=blkdiag(R_theta,R_beta);
    Rn_a=blkdiag(Rn_a1);
    %站址误差
    deta_s=deta_more(4)*z_error_more(2,big_big_loop);
    power_s_error=[1 1 1 3 3 3 5 5 5 10 10 10 5 5 5 3 3 3].';
    Rs_a=deta_s^2*diag(power_s_error);
    %     观测站定位的先验均方根误差
    prior_cov_s(big_big_loop,1)=sqrt(trace(Rs_a));
    
    F1=F_u.'*inv(Rn_a)*F_u;
    
    CRLB=inv([F1]);
    a=CRLB(1,1)+CRLB(2,2)+CRLB(3,3);
    
    crb_source1(big_big_loop,1)=sqrt(a);
    CRLB_sensor(big_big_loop,1)=sqrt(trace(CRLB)-(a));
end
CRLB_u=[crb_source1];
%%
figure(1)
plot(20*log10(z_error_more(1,:)),CRLB_u_sensor_erros(:,1),'r.-')
hold on;
plot(20*log10(z_error_more(1,:)),CRLB_u(:,1),'b*-')
hold on;
legend('存在站址误差','不存在站址误差')
xlabel('角度误差')
ylabel('RMSE')

%%
%******CRLB随站址误差变化结果
%%
clc;
close all;
clear all;
% 测向站数量
M=6;
% 目标数量
N=1;
%观测站位置
s1=[1500 -1900 1400].';
s2=[-1600 2000 1700].';
s3=[1800 -1500 -1600].';
s4=[-1200 1400 -2000].';
s5=[1700 1600 2000].';
s6=[-1800 -1400 -1800].';
for i=1:1:M
    str_si=['s',mat2str(i)];
    si_value=eval(str_si);
    s_true(3*(i-1)+1:3*(i-1)+3,1)=si_value;
end
% 目标位置
u1=[4000 4000 3000].';
for i=1:1:N
    str_ui=['u',mat2str(i)];
    ui_value=eval(str_ui);
    u_true(3*(i-1)+1:3*(i-1)+3,1)=ui_value;
end
z_error_more(2,:)=10.^[0:0.1:1.5];
z_error_more(1,:)=10.^[zeros(length(z_error_more),1)];
deta_more=[0.001,0.001,1,10];
% 站址误差功率不一样
power_s_error=[1 1 1 3 3 3 5 5 5 10 10 10 5 5 5 3 3 3].';
%% 存在站址误差
%% CRLB
for i=1:1:M
    eval(['s',mat2str(i),'(:,1)',' = ','s_true(3*(i-1)+1:3*(i-1)+3,1)',';']);
end
for i=1:1:N
    eval(['u',mat2str(i),'(:,1)',' = ','u_true(3*(i-1)+1:3*(i-1)+3,1)',';']);
end
%%     对目标位置的导数
% 方位角、仰角对目标位置的导数
diff_theta_ui=zeros(M,3,N);
diff_beta_ui=zeros(M,3,N);
%方位角、仰角对目标观测站的导数
diff_theta_si=zeros(M,3*M,N);
diff_beta_si=zeros(M,3*M,N);
%%
for j=1:1:N
    for i=1:1:M
        str_si=['s',mat2str(i)];
        si_value=eval(str_si);
        str_ui=['u',mat2str(j)];
        ui_value=eval(str_ui);
        
        ui_si=ui_value-si_value;
        ai2=sqrt(ui_si(1)^2+ui_si(2)^2);
        %     方位角对目标位置的导数
        diff_theta_ui_1=-ui_si(2)/ai2^2;
        diff_theta_ui_2=ui_si(1)/ai2^2;
        diff_theta_ui(i,:,j)=[diff_theta_ui_1,diff_theta_ui_2,0];
        %     仰角对目标位置的导数
        diff_beta_ui_1=-ui_si(1)*ui_si(3)/ai2/(osjl(ui_value,si_value))^2;
        diff_beta_ui_2=-ui_si(2)*ui_si(3)/ai2/(osjl(ui_value,si_value))^2;
        diff_beta_ui_3=ai2/(osjl(ui_value,si_value))^2;
        diff_beta_ui(i,:,j)=[diff_beta_ui_1,diff_beta_ui_2,diff_beta_ui_3];
        
        % 方位角、仰角对目标观测站位置的导数
        diff_theta_si(i,3*(i-1)+1:3*(i-1)+3,j)=-diff_theta_ui(i,:,j);
        diff_beta_si(i,3*(i-1)+1:3*(i-1)+3,j)=-diff_beta_ui(i,:,j);
    end
end

F_u1=[diff_theta_ui(:,:,1);diff_beta_ui(:,:,1)];
F_u=blkdiag(F_u1);

F_s1=[diff_theta_si(:,:,1);diff_beta_si(:,:,1)];
F_s=[F_s1.'].';

for big_big_loop=1:1:length(z_error_more)
    %方位角
    deta_theta=deta_more(1)*z_error_more(1,big_big_loop);
    R_theta=deta_theta^2*eye((M),(M));
    %仰角
    deta_beta=deta_more(2)*z_error_more(1,big_big_loop);
    R_beta=deta_beta^2*eye((M),(M));
    
    Rn_a1=blkdiag(R_theta,R_beta);
    Rn_a=blkdiag(Rn_a1);
    %站址误差
    deta_s=deta_more(4)*z_error_more(2,big_big_loop);
    power_s_error=[1 1 1 3 3 3 5 5 5 10 10 10 5 5 5 3 3 3].';
    Rs_a=deta_s^2*diag(power_s_error);
    %     观测站定位的先验均方根误差
    prior_cov_s(big_big_loop,1)=sqrt(trace(Rs_a));
    
    F1=F_u.'*inv(Rn_a)*F_u;
    F2=F_u.'*inv(Rn_a)*F_s;
    F3=F_s.'*inv(Rn_a)*F_s+inv(Rs_a);
    
    CRLB=inv([F1,F2;F2.',F3]);
    a=CRLB(1,1)+CRLB(2,2)+CRLB(3,3);
    
    crb_source1(big_big_loop,1)=sqrt(a);
    CRLB_sensor(big_big_loop,1)=sqrt(trace(CRLB)-(a));
end
CRLB_u_sensor_erros=[crb_source1];
%% 不存在站址误差
for big_big_loop=1:1:length(z_error_more)
    %方位角
    deta_theta=deta_more(1)*z_error_more(1,big_big_loop);
    R_theta=deta_theta^2*eye((M),(M));
    %仰角
    deta_beta=deta_more(2)*z_error_more(1,big_big_loop);
    R_beta=deta_beta^2*eye((M),(M));
    
    Rn_a1=blkdiag(R_theta,R_beta);
    Rn_a=blkdiag(Rn_a1);
    %站址误差
    deta_s=deta_more(4)*z_error_more(2,big_big_loop);
    power_s_error=[1 1 1 3 3 3 5 5 5 10 10 10 5 5 5 3 3 3].';
    Rs_a=deta_s^2*diag(power_s_error);
    %     观测站定位的先验均方根误差
    prior_cov_s(big_big_loop,1)=sqrt(trace(Rs_a));
    
    F1=F_u.'*inv(Rn_a)*F_u;
    
    CRLB=inv([F1]);
    a=CRLB(1,1)+CRLB(2,2)+CRLB(3,3);
    
    crb_source1(big_big_loop,1)=sqrt(a);
    CRLB_sensor(big_big_loop,1)=sqrt(trace(CRLB)-(a));
end
CRLB_u=[crb_source1];
%%
figure(1)
plot(20*log10(z_error_more(2,:)),CRLB_u_sensor_erros(:,1),'r.-')
hold on;
plot(20*log10(z_error_more(2,:)),CRLB_u(:,1),'b*-')
hold on;
legend('存在站址误差','不存在站址误差')
xlabel('站址误差')
ylabel('RMSE')

参考：不存在站址误差的DOA定位CRLB
不存在站址误差的TDOA定位CRLB