DOA观测模型
DOA即direction of arrival，包括目标到达观测站的方位角和仰角
下图所示为在三维空间直角坐标系下，DOA定位示意图
利用M个观测站进行定位
假设第i个观测站的坐标向量为：
s
i
o
=
[
s
i
x
o
,
s
i
y
o
,
s
i
z
o
]
T
\ \boldsymbol s^o _i =[s^o_{ix},s^o_{iy},s^o_{iz}]^T
sio=[sixo,siyo,sizo]T ，目标的位置向量为：
u
o
=
[
u
x
o
,
u
y
o
,
u
z
o
]
T
\ \boldsymbol u^o =[u^o_{x},u^o_{y},u^o_{z}]^T
uo=[uxo,uyo,uzo]T 。因此，方位角和仰角可以由
θ
i
o
=
a
r
c
t
a
n
(
u
y
o
−
s
i
y
o
u
x
o
−
s
i
x
o
)
\theta^o_i =arctan \left( \frac{u^o_{y}-s^o_{iy}} {{u^o_{x}-s^o_{ix}}} \right)
θio=arctan(uxo−sixouyo−siyo)
ϕ
i
o
=
a
r
c
t
a
n
(
u
z
o
−
s
i
z
o
(
u
x
o
−
s
i
x
o
)
2
+
(
u
y
o
−
s
i
y
o
)
2
)
)
\phi ^o_i =arctan \left( \frac{u^o_{z}-s^o_{iz}} {\sqrt{(u^o_{x}-s^o_{ix})^2+(u^o_{y}-s^o_{iy})^2}} ) \right)
ϕio=arctan
(uxo−sixo)2+(uyo−siyo)2uzo−sizo)
其中
θ
i
o
\theta^o_i
θio ，
ϕ
i
o
\phi ^o_i
ϕio分别表示目标到达第i个观测站的方位角和仰角的真实值。由于在实际中我们仅能够获得带有测量噪声的测量值，因此建立如下观测模型：
θ
=
[
θ
1
o
,
θ
2
o
,
.
.
.
,
θ
M
o
]
T
+
[
δ
θ
1
,
δ
θ
2
,
.
.
.
,
δ
θ
M
]
T
=
θ
o
+
n
1
\boldsymbol \theta =[\theta^o_1,\theta^o_2,...,\theta^o_M]^T+[\delta \theta_1,\delta \theta_2,...,\delta \theta_M]^T=\boldsymbol \theta^o+n_1
θ=[θ1o,θ2o,...,θMo]T+[δθ1,δθ2,...,δθM]T=θo+n1
ϕ
=
[
ϕ
1
o
,
ϕ
2
o
,
.
.
.
,
ϕ
M
o
]
T
+
[
δ
ϕ
1
,
δ
ϕ
2
,
.
.
.
,
δ
ϕ
M
]
T
=
ϕ
o
+
n
2
\boldsymbol \phi =[\phi ^o_1,\phi ^o_2,...,\phi ^o_M]^T+[\delta \phi _1,\delta \phi _2,...,\delta \phi _M]^T=\boldsymbol \phi^o+n_2
ϕ=[ϕ1o,ϕ2o,...,ϕMo]T+[δϕ1,δϕ2,...,δϕM]T=ϕo+n2
θ
o
\theta^o
θo和
ϕ
o
\phi^o
ϕo分别表示分别表示方位角和仰角的真实值向量。
n
1
\boldsymbol n_1
n1,
n
2
\boldsymbol n_2
n2为对应的噪声向量，其均服从零均值高斯分布，协方差矩阵分别为
Q
1
\boldsymbol Q_1
Q1，
Q
2
\boldsymbol Q_2
Q2。
因此可以得到DOA观测向量为：
z
=
z
o
+
n
=
[
(
θ
o
)
T
,
(
ϕ
o
)
T
]
T
+
[
n
1
T
,
n
2
T
]
T
\boldsymbol z =\boldsymbol z^o+\boldsymbol n=[(\boldsymbol \theta^o)^T,(\boldsymbol \phi^o)^T]^T+[\boldsymbol n^T_1,\boldsymbol n^T_2]^T
z=zo+n=[(θo)T,(ϕo)T]T+[n1T,n2T]T
其中
n
n
n的协方差矩阵为
Q
=
b
l
k
d
i
a
g
(
Q
1
,
Q
2
)
\boldsymbol Q=blkdiag(\boldsymbol Q_1,\boldsymbol Q_2)
Q=blkdiag(Q1,Q2)。
CRLB
位置参数向量为
u
o
=
[
u
x
o
,
u
y
o
,
u
z
o
]
T
\ \boldsymbol u^o =[u^o_{x},u^o_{y},u^o_{z}]^T
uo=[uxo,uyo,uzo]T，测量向量为
z
\boldsymbol z
z。因此
z
\boldsymbol z
z关于
u
o
\ \boldsymbol u^o
uo的对数似然函数为：
l
n
(
p
(
z
∣
u
o
)
)
=
κ
−
(
1
/
2
)
(
z
−
z
o
)
T
Q
−
1
(
z
−
z
o
)
ln(p(\boldsymbol z|\boldsymbol u^o))=\kappa-(1/2)(\boldsymbol z-\boldsymbol z^o)^T\boldsymbol Q^{-1}(\boldsymbol z-\boldsymbol z^o)
ln(p(z∣uo))=κ−(1/2)(z−zo)TQ−1(z−zo)
其中
κ
\kappa
κ表示常数项
因此，Fisher信息矩阵可以表示为：
F
=
E
(
(
∂
l
n
(
p
(
z
∣
u
o
)
)
∂
(
u
o
)
T
)
(
∂
l
n
(
p
(
z
∣
u
o
)
)
∂
(
u
o
)
T
)
T
)
=
(
∂
z
o
∂
(
u
o
)
T
)
Q
−
1
(
∂
z
o
∂
(
u
o
)
T
)
T
F=E \left( \left( \frac {\partial ln(p(\boldsymbol z|\boldsymbol u^o))}{\partial (\boldsymbol u^o)^T} \right) \left( \frac {\partial ln(p(\boldsymbol z|\boldsymbol u^o))}{\partial (\boldsymbol u^o)^T} \right)^T \right)= \left( \frac {\partial \boldsymbol z^o}{\partial (\boldsymbol u^o)^T} \right) \boldsymbol Q^{-1} \left( \frac {\partial \boldsymbol z^o}{\partial (\boldsymbol u^o)^T} \right)^T
F=E((∂(uo)T∂ln(p(z∣uo)))(∂(uo)T∂ln(p(z∣uo)))T)=(∂(uo)T∂zo)Q−1(∂(uo)T∂zo)T
CRLB为Fisher矩阵的逆，因此CRLB表示为：
C
R
L
B
=
(
(
∂
z
o
∂
(
u
o
)
T
)
Q
−
1
(
∂
z
o
∂
(
u
o
)
T
)
T
)
−
1
CRLB=\left( \left( \frac {\partial \boldsymbol z^o}{\partial (\boldsymbol u^o)^T} \right) \boldsymbol Q^{-1} \left( \frac {\partial \boldsymbol z^o}{\partial (\boldsymbol u^o)^T} \right)^T \right)^{-1}
CRLB=((∂(uo)T∂zo)Q−1(∂(uo)T∂zo)T)−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
仿真结果如下所示：
欢迎私信交流无源定位、阵列信号处理
1306
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
