雷达导引头制导仿真

雷达导引头搜索目标的仿真过程
雷达制导分为搜索和跟踪。
首先是参数设置，雷达工作参数；导弹位置坐标、姿态角、速度；目标距离、方位、俯仰；相控阵天线阵列
第一步，将目标转换到导弹坐标系下，得出弹目连线对应的角度

function [xx,yy,zz] = zuobiao(theta,phi,gamma,target)
%角度转弧度
theta = deg2rad(theta);
phi = deg2rad(phi);
gamma = deg2rad(gamma);
%坐标转换
%地面坐标至弹体坐标
L = [cos(phi) * cos(theta) , sin(phi) , -cos(phi)*sin(theta)
   -sin(phi)*cos(theta)*cos(gamma)+sin(theta)*sin(gamma),  cos(phi)*cos(gamma),  sin(phi)*sin(theta)*cos(gamma)+cos(theta)*sin(gamma) 
   sin(phi)*cos(theta)*sin(gamma)+sin(theta)*cos(gamma),  -cos(phi)*sin(gamma),  -sin(phi)*sin(theta)*sin(gamma)+cos(theta)*cos(gamma)];
Pos= L * target.';
xx = Pos(1);
yy = Pos(2);
zz = Pos(3);
end

第二步，计算相控阵天线图的和波束、俯仰差波束、防卫差波束

lam = 0.03;
M_ant = 21;
N_ant = 21;
dy = 0.015;
dz = 0.015;
slope = 0;  %天线阵面与弹体夹角

drz = 2*pi*dz/lam;
dry = 2*pi*dy/lam;
M = 300;
N =300;
theta0 = 40;
phi0 = 30;
theta = linspace(-180,180,M);
phi = linspace(-90,90,N);
F_the = zeros(M,N);
F_phi = zeros(M,N);
for iz = 1:M_ant
    for iy = 1:N_ant
        F_the = F_the + exp(1j*(iz-1)*drz*(cosd(phi).'*sind(theta)-cosd(90-theta0)));
        F_phi = F_phi + exp(1j*(iy-1)*dry*(sind(phi).' * ones(1,M)-cosd(90-phi0)));
    end
end
F = F_the .* F_phi;
[phi2,theta2]=meshgrid(phi,theta);
figure;
subplot(1,2,1);mesh(theta2,phi2,abs(F));shading interp;title('球坐标系');
[x,y,z]=sph2cart(deg2rad(theta2),deg2rad(phi2),abs(F));  %球面到笛卡尔坐标系转换
subplot(1,2,2);surf(x,y,z);axis equal;shading interp;title('直角坐标系');
%xlim([0,100])
figure;
subplot(1,2,1);polarplot(theta/180*pi,abs(F(M/2,:)));title('方位40度');
subplot(1,2,2);polarplot(phi/180*pi,abs(F(:,N/2)));title('俯仰30度');

第三步，循环扫描，建立回波模型
第四步，脉冲压缩、选取距离波们、MTI、MTD、CFAR

%% 回波生成
%波束扫描循环
for j = 1:5
    for i = 1:5
        %         (j-1)*5+i
        theta1 = theta_start +(i-1) * theta_i;
        phi1  = phi_start + (j-1) * phi_i;
        M = 180;
        N = 180;
        % 和波束
        im =ceil( (theta1+180)/360*M);
        in = ceil((phi1+90)/180*N);
        F1 =abs(F_he(im,in))/F0;
        % 方位差波束
        F2 = abs(F_cha(im,in))/F_cha0;
        % 俯仰差波束
        F3 = abs(F_cha_fu(im,in))/F_cha_fu0;

        T = 1/fs:1/fs:N*PRT;
        MM = length(T);
        s = zeros(1,MM);
        %和波束、方位差波束、俯仰差波束
if  F1 < 0.5
        disp(num2str((j-1)*5+i))
        disp('距离：0m')
        disp('速度：0m/s')
        disp(' ')
else


        Gr = F1;
        M = floor(PRT*fs);
        t = (1:M)/fs;
        %17个脉冲
        for n = 1:N
            tmp = zeros(1,M);
            ind_tianxian = (i-1)*5+j;
            Rt = R0-(V_dao+V_tar)*(ind_tianxian-1)*T_n;
            Pr = Pt*Gt*Gr*lam^2*sigma/(4*pi)^3/R0^4/Loss1/Loss2;

            %散射点累加
            for k = 1:numpoint
                tau = 2*Rt/c;
                tmp = tmp + Pr * 1e12 * rectpuls((t-tau),Tr) *exp(-1j*2*pi*f0*tau) .* exp(1j*pi*Kr*(t-tau).^2);
            end
            s(1,(n-1)*MM/N+1:(n-1)*MM/N+M) = tmp ;
        end
        noise = sqrt(K*B) *1e-3* wgn(1,length(T),K*Tn*B,'complex');
        s = Gain*(s+noise);
        s_echo = s;

        for ik = 1:length(T)
            uu = real(s(1,ik));
            if abs(uu) >ui
                uu = (uu/abs(uu))*k*log(uu/ui);
            else
                uu = k*uu;
            end
            s_echo(1,ik) = uu;
        end


        %% 距离波门建立，起始81 结束800
        s_ad = s_echo;
        n_gate =8000-81+1;
        N = 17;
        L = ceil(MM/N);
        for ik = 1:N
            s_rg((ik-1)*n_gate+1:ik*n_gate) = s_ad((ik-1)*L+81:(ik-1)*L+8000);
        end


        %脉冲压缩模块
        tr = (1:PRT*fs)'/fs;
        s_ref = rectpuls(tr,Tr) .* exp(1j*pi*Kr*tr.^2);
        ht = conj(fliplr(s_ref));
        s_rd = conv(ht,s_rg);
        s_rd = s_rd(M/2:end-M/2);


        %%     MTD
        vd = lam/PRT/4;%最大不模糊速度
        s_mti = s_rd(1,n_gate+1:end)- s_rd(1,1:end-n_gate);
        s_mti = s_rd(1,n_gate+1:end);

        N = 16;
        win = taylorwin(N);
        for m = 1:N
            ww(m,:) = exp(-1i*2*pi*m*(0:N-1)/N) .* win';
        end
        hd = fftshift(fft(ww,n_gate,2),2);


        s_mat = reshape(s_mti,n_gate,N);
        [~,ir] = find(abs(s_mat(1,:)) == max(abs(s_mat(1,:))));
        r = ir * c/2/B;
        f = (-N/2:N/2-1)/N/PRT;
        s_mat2 = fftshift(fft(s_mat,N,2),2);


        H = hd.' .* s_mat;
        s_mtd = reshape(H,1,N*n_gate);

        %% 对MTD和通道信号处理，取模值最大点的位置

        s_new = [zeros(1,(pro_len+win_len)),s_mtd,zeros(1,(pro_len+win_len))];
        s_cfar = s_mtd;
        m_cfarRate = 0.8;
        cfarRate = 10*sqrt(-2*log10(m_cfarRate));
        cfarRate = 1e-3          ;
        for jk = (pro_len+win_len)+1:length(s_mtd)+(pro_len+win_len)
            x1 =  sum(abs(s_new(jk-(pro_len+win_len):jk-pro_len-1)))/win_len;
            x2 =  sum(abs(s_new(jk+pro_len+1:jk+(pro_len+win_len))))/win_len;
            if (x1+x2)/2 > cfarRate
                s_cfar(jk-(pro_len+win_len)) = s_mtd(jk-(pro_len+win_len));
            else
                s_cfar(jk-(pro_len+win_len)) = 0;
            end
        end


        H_cfar = reshape(s_cfar,n_gate,N);



        %搜索模块
        [ind, ind_v] = find(H_cfar == max(max(H_cfar)) );


%         A = abs(H_cfar(ind,ind_v));
%         B = abs(H_cfar(ind,ind_v-1));
%         C = abs(H_cfar(ind,ind_v+1));
%         if B == C
%             ind_v2 = ind_v;
%         else
%             if B > C
%                 ind_v2 = ind_v - 0.5 - 0.5 *(B-A)/(B+A);
%             else
%                 ind_v2 = ind_v + 0.5 + 0.5 *(C-A)/(C+A);
%             end
%         end

        fd =  ind_v/N/PRT;
        v = lam*fd/2;

        disp(num2str((j-1)*5+i))
        disp(['距离：',num2str(r),'m'])
        disp(['速度：',num2str(v),'m/s'])
        disp(' ')
    end


        %{
        % 绘图
         figure;
         subplot(4,3,1);plot(real(s_echo));title('接收机回波');xlabel('t/s');ylabel('幅值');
         subplot(4,3,2);plot(real(s_rg));title('距离波门选取');xlabel('t/s');ylabel('幅值');
         subplot(4,3,3);plot(abs(s_rd));title('脉压');xlabel('t/s');ylabel('幅值');
         subplot(4,3,4);plot(abs(s_mtd));title('MTD后');xlabel('t/s');ylabel('幅值');
         subplot(4,3,5);plot((1:n_gate)/n_gate,abs(hd.'));title('MTD滤波器组频率响应');xlabel('归一化频率/Hz');ylabel('幅值')
         subplot(4,3,6);plot(abs(s_cfar));title('cfar后');

         
         subplot(4,1,3);imagesc(f,1:n_gate,abs(s_mat2));title('距离-多普勒图');ylabel('距离维');xlabel('多普勒频率');
         subplot(4,2,7);surf(abs(H));shading interp;title('MTD结果');
         subplot(4,2,8);mesh(abs(H_cfar));shading interp;title('CFAR结果');
        %}
    end

end

参考文献：
郭明明，雷达导引头系统建模仿真与性能评估[D]西安电子科技大学