DOA观测模型

DOA即direction of arrival，包括目标到达观测站的方位角和仰角
下图所示为在三维空间直角坐标系下，DOA定位示意图

利用M个观测站进行定位
假设第i个观测站的坐标向量为：${\mathbf{s}}_i^o=[s_{ix}^o,s_{iy}^o,s_{iz}^o{]}^T$ ，目标的位置向量为：${\mathbf{u}}^o=[u_x^o,u_y^o,u_z^o{]}^T$ 。因此，方位角和仰角可以由

$${\theta }_i^o=arctan\left(\frac{u_y^o-s_{iy}^o}{u_x^o-s_{ix}^o}\right)$$
$${\varphi }_i^o=arctan\left(\frac{u_z^o-s_{iz}^o}{\sqrt{(u_x^o-s_{ix}^o{)}^2+(u_y^o-s_{iy}^o{)}^2}})\right)$$
其中${\theta }_i^o$ ，${\varphi }_i^o$分别表示目标到达第i个观测站的方位角和仰角的真实值。由于在实际中我们仅能够获得带有测量噪声的测量值，因此建立如下观测模型：
$$\mathbf{\theta }=[{\theta }_1^o,{\theta }_2^o,...,{\theta }_M^o{]}^T+[\delta {\theta }_1,\delta {\theta }_2,...,\delta {\theta }_M{]}^T={\mathbf{\theta }}^o+n_1$$
$$\mathbf{\varphi }=[{\varphi }_1^o,{\varphi }_2^o,...,{\varphi }_M^o{]}^T+[\delta {\varphi }_1,\delta {\varphi }_2,...,\delta {\varphi }_M{]}^T={\mathbf{\varphi }}^o+n_2$$
${\theta }^o$和${\varphi }^o$分别表示分别表示方位角和仰角的真实值向量。${\mathbf{n}}_1$,${\mathbf{n}}_2$为对应的噪声向量，其均服从零均值高斯分布，协方差矩阵分别为${\mathbf{Q}}_1$，${\mathbf{Q}}_2$。
因此可以得到DOA观测向量为：
$$\mathbf{z}={\mathbf{z}}^o+\mathbf{n}=[({\mathbf{\theta }}^o{)}^T,({\mathbf{\varphi }}^o{)}^T{]}^T+[{\mathbf{n}}_1^T,{\mathbf{n}}_2^T{]}^T$$
其中$n$的协方差矩阵为$\mathbf{Q}=blkdiag({\mathbf{Q}}_1,{\mathbf{Q}}_2)$。

CRLB

位置参数向量为${\mathbf{u}}^o=[u_x^o,u_y^o,u_z^o{]}^T$，测量向量为$\mathbf{z}$。因此$\mathbf{z}$关于${\mathbf{u}}^o$的对数似然函数为：
$$ln(p(\mathbf{z}|{\mathbf{u}}^o))=\kappa -(1/2)(\mathbf{z}-{\mathbf{z}}^o{)}^T{\mathbf{Q}}^{-1}(\mathbf{z}-{\mathbf{z}}^o)$$
其中$\kappa$表示常数项

因此，Fisher信息矩阵可以表示为：
$$F=E\left(\left(\frac{\partial ln(p(\mathbf{z}|{\mathbf{u}}^o))}{\partial ({\mathbf{u}}^o{)}^T}\right){\left(\frac{\partial ln(p(\mathbf{z}|{\mathbf{u}}^o))}{\partial ({\mathbf{u}}^o{)}^T}\right)}^T\right)=\left(\frac{\partial {\mathbf{z}}^o}{\partial ({\mathbf{u}}^o{)}^T}\right){\mathbf{Q}}^{-1}{\left(\frac{\partial {\mathbf{z}}^o}{\partial ({\mathbf{u}}^o{)}^T}\right)}^T$$
CRLB为Fisher矩阵的逆，因此CRLB表示为：
$$CRLB={\left(\left(\frac{\partial {\mathbf{z}}^o}{\partial ({\mathbf{u}}^o{)}^T}\right){\mathbf{Q}}^{-1}{\left(\frac{\partial {\mathbf{z}}^o}{\partial ({\mathbf{u}}^o{)}^T}\right)}^T\right)}^{-1}$$

仿真

clc;
close all;
clear all;
%%
% 测向站数量
M=6;
% 目标数量
N=1;
% 6个测向站的位置坐标,单位米
s1=[1200 1800 200];
s2=[-1500 -800 150];
s3=[1400 -600 -200];
s4=[-800 1200 120];
s5=[1300 -800 -250];
s6=[-1000 1600 -150];
% 目标位置,单位米
u1=[8000 6800 3000];
%%
for i=1:M
    %     组合字符串，得到接收站的变量名，以便循环
    str_si=['s',mat2str(i)];
    %    返回变量名对应的值
    si_value=eval(str_si);
    %    测向站到目标的距离
    distance(i)=sqrt(sum((u1-si_value).^2));
    %    测向站到目标的平均距离
    distance_average=sum(distance)/M;
    %     方位角真实值
    theta0(i)=atan((u1(1)-si_value(1))/(u1(2)-si_value(2)));
    %     仰角真实值
    beta0(i)=atan((u1(3)-si_value(3))/(sqrt((u1(1)-si_value(1))^2+(u1(2)-si_value(2))^2)));
end
%%
G0=zeros(2*M,3);
h0=zeros(2*M,1);
for i=1:1:M
    G0(2*(i-1)+1,:)=[cos(theta0(i)),-sin(theta0(i)),0];
    G0(2*(i-1)+2,:)=[sin(theta0(i))*sin(beta0(i)),cos(theta0(i))*sin(beta0(i)),-cos(beta0(i))];
    
    %     组合字符串，得到接收站的变量名，以便循环
    str_si=['s',mat2str(i)];
    %    返回变量名对应的值
    si_value=eval(str_si);
    h0(2*(i-1)+1,:)=si_value(1)*cos(theta0(i))-si_value(2)*sin(theta0(i));
    h0(2*(i-1)+2,:)=si_value(1)*sin(theta0(i))*sin(beta0(i))+si_value(2)*cos(theta0(i))*sin(beta0(i))-si_value(3)*cos(beta0(i));
end
result=G0*u1.'-h0;
%%
% 误差的标准差，弧度
deta_theta_more=(0.1:0.1:1).*pi/180;
big_loop_numer=length(deta_theta_more);

% 每一个测向标准差下的CRLB
CRLB=zeros(length(deta_theta_more));
%%
F_s_t0=zeros(2*M,3);
for i=1:M
    %     组合字符串，得到接收站的变量名，以便循环
    str_si=['s',mat2str(i)];
    %    返回变量名对应的值
    si_value=eval(str_si);
    
    F_s_t0(2*i-1,1)=cos(theta0(i))/sqrt(sum((u1-si_value).^2));
    F_s_t0(2*i-1,2)=-sin(theta0(i))/sqrt(sum((u1-si_value).^2));
    F_s_t0(2*i,1)=-(sin(theta0(i))*sin(beta0(i)))/sqrt(sum((u1-si_value).^2));
    F_s_t0(2*i,2)=-(cos(theta0(i))*sin(beta0(i)))/sqrt(sum((u1-si_value).^2));
    F_s_t0(2*i,3)=cos(beta0(i))/sqrt(sum((u1-si_value).^2));
end
%%
for big_loop=1:length(deta_theta_more)
    deta_theta=deta_theta_more(big_loop);
    cov_z=deta_theta^2*eye(2*M);
    % CRB界
    CRLB_sum(big_loop)=sqrt(trace((F_s_t0'*inv((cov_z))*F_s_t0)^-1));
end
figure
plot(deta_theta_more./(pi/180),CRLB_sum,'b.-');
xlabel('角度误差(°)')
ylabel('RMSE\m')

function [distance]=osjl(object,source)
% 输入均为列向量
distance=sqrt(sum((object-source).^2));
end

仿真结果如下所示：

欢迎私信交流无源定位、阵列信号处理