单星单天线仅测信号到达时间定位模型方法

参考《空间电子侦察定位原理》第6章

算法原理

代码实现

%% start
clc;
close all;
clear all;

%% 参数设置
% 通用参数
c = 3e8 ;  % 信号传播速度 单位：m/s

% 单颗卫星参数
v = 500e3 ; % 卫星速度 单位：m/s

V = [v,0,0] ; % 速度矢量  单位：m/s 
H = 7e3 ; % 卫星运行高度 单位：m

Sa_loc = [0,0,H] ; % 卫星的初始位置 依次对应坐标的X Y Z ,Z为高度


% 射频源参数
Tar_loc = [37e3,37e3,0];
fc   = 500e9 ; % 射频源载频
Tr   = 500e-3; % 脉冲重复周期 单位：s
PRF  = 1/Tr ;  % 脉冲重复周期 单位：Hz

%% ==============================多次测量=================================
% 设置测量时间
t_m = [0,1,10,100,1000,10000,1000e3] ; % 设置测量时间点
n_m = length(t_m) ;
A = zeros(n_m,3) ; % 卫星加速度 单位:m/s2
Sa = zeros(n_m,3);
% 生成测量位置坐标矩阵
for i = 1:1:n_m
    Sa(i,:)=[Sa_loc(1)+v*t_m(i),0,H]; % 卫星沿着X轴方向飞行
end

% 计算距离矢量矩阵
B = repmat(Tar_loc, n_m, 1);

R = Sa - B ; % 距离矢量矩阵

% 对每行进行平方和开根号操作求出距离标量
r = sqrt(sum(R.^2, 2)) ;

% 生成单位方向矢量
U = R ./ repmat(r, 1, size(R, 2));

% 相对速度计算
V = repmat(V, n_m, 1);
V_1 = abs(U .* V) ; % 相对速度矩阵
R_1 = sum(V_1, 2) ; % 相对速度标量列

% 相对加速度计算
R_2 = zeros(n_m,1);
for i = 1:1:n_m
    R_2(i,:) = (V(i,:)*V(i,:)' + R(i,:)*A(i,:)')/r(i)-(R(i,:)*V(i,:)')/(r(i)^3) ;
end

% 生成各个脉冲的TOA
K = 5 ; % 连续脉冲个数
TK = zeros(n_m,K);
a  = zeros(n_m,1);
b  = zeros(n_m,1);
for i = 1:1:n_m  
    for k = 1:1:K
        TK(i,k) = k*Tr/(1-R_1(i)/c) + ((R_2(i)/c)*(k*Tr)^2)/(2*(1-R_1(i)/c)^3);
    end
    a(i) = 1/(1-R_1(i)/c);
    b(i) = (R_2(i)/c)/(2*(1-R_1(i)/c)^3);
end


% 解方程
KT1 = 1:1:K;
KT = [KT1.^2;KT1];
a_c  = zeros(n_m,1);  % 方程求解出的a参数列
b_c  = zeros(n_m,1);  % 方程求解出的b参数列
R_2_C = zeros(n_m,1); % 方程求解出的加速度参数列
for j = 1:1:n_m
    theta = pinv(KT*KT')*(KT)*(TK(j,:)');
    % 方程求解出的a b 加速度参数
    b_c(j) = theta(1)/(Tr^2) ;
    a_c(j) = theta(2)/Tr ;
    R_2_C(j) = 2*b_c(j)*c/(a_c(j)^3); 
end
p=0;
for x = 30e3: 1000 : 50e3
    p = p + 1 ;
    q = 0 ;
    for y = 30e3: 1000 : 50e3

        q = q + 1 ;

        Tar = [x,y,0]; % 假设目标就是在地面 z = 0
        % 计算距离矢量矩阵
        B = repmat(Tar, n_m, 1);

        R_c = Sa - B ; % 距离矢量矩阵

        % 对每行进行平方和开根号操作求出距离标量
        r_c = sqrt(sum(R_c.^2, 2)) ;


        % 相对加速度计算
        R_2_C1 = zeros(n_m,1);
        for i = 1:1:n_m
            R_2_C1(i,:) = (V(i,:)*V(i,:)' + R_c(i,:)*A(i,:)')/r_c(i)-(R_c(i,:)*V(i,:)')/(r_c(i)^3) ;
        end
        % 频率变化率计算

       
        % 计算损失函数
        Loss(p,q)= abs(sum((R_2_C - R_2_C1).^2));
    end
end


figure;mesh(Loss);
figure;mesh(max(max(Loss))-Loss );
[a,b] = find(Loss == min(min(Loss)))

仿真结果分析

• 没有使用经纬度转换，假设的是在一个三维坐标系中，目标在坐标系的地面上（Z=0）
• 损失误差图不会出现很明显的峰值但是确实是有损失最小的点
•