cegao.m

机载双基雷达

%%========================================================
clear;clc;close all;
%%========================================================
%%参数——常量
C=3e8;                           %电磁波速度

%%参数——雷达规格
Fc=3e9;                          %载波频率
lambda=C/Fc;                     %波长 

%%参数——目标区域
X0=500;                          %方位向范围：[-X0,X0]
Y=3000;                         %测绘中心Y坐标
Y1=Y;
Y0=500;                          %成像宽度：2*Y0
Z0=500;                           %高度范围：[-Z0,Z0]
%%参数
theta=0;                         %平台波束斜视角
theta=theta*pi/180;              %rad
V=100;                           %平台速度：100 m/s
H=4000;                          %平台高度：4 km
Rb=sqrt(H^2+Y^2);                %区域中心斜距
Rc=Rb/cos(theta);                %目标区域距离向范围：[Rc-Y0,Rc+Y0]

%%参数——雷达天线
D=4;                             %天线方位向尺寸
Lsar=lambda*Rc/D;                %SAR合成孔径长度
Tsar=Lsar/V;                     %SAR合成孔径时间

%%参数——慢时间域
Ka=-2*(V*cos(theta))^2/lambda/Rc;            %doppler调频斜率
Ba=abs(Ka*Tsar);                %doppler调频频宽
PRF=Ba;                         %脉冲重复频率
PRT=1/PRF;                      %脉冲重复时间
ds=PRT;                         %慢时域采样间距
                                %为了进行FFT，对PRT、PRF、ds进行运算并更新
Nslow=ceil((2*X0+Lsar)/V/ds);   %慢时域采样点数
Nslow=2^nextpow2(Nslow);        %使采样点数为2的幂；N点

%Nslow=1024;

sn=linspace(-(X0+Lsar/2)/V,(X0+Lsar/2)/V,Nslow);%慢时域离散时间序列，1*N矩阵
PRT=(2*X0+Lsar)/V/Nslow;        %更新PRT、PRF、ds
PRF=1/PRT;
ds=PRT;
%%参数——快时域
Tr=5e-6;                        %Chirp脉冲持续时间：5 us
Br=50e6;                        %chirp信号调频带宽：50 MHz
Kr=Br/Tr;                       %chirp信号调频斜率
Fsr=2*Br;                       %快时域采样频率
dt=1/Fsr;                       %快时域采样间距
Rmin=sqrt((Y-Y0)^2+(H-Z0)^2);        %最小斜距
Rmax=sqrt((Y+Y0)^2+(H+Z0)^2+(Lsar/2)^2);     %最大斜距
Nfast=ceil(2*(Rmax-Rmin)/C/dt+Tr/dt);   %快时域采样点数
Nfast=2^nextpow2(Nfast);                %使采样点数为2的幂；M点
tm=linspace(2*Rmin/C,2*Rmax/C+Tr,Nfast);%快时域离散时间序列，1*M矩阵
dt=(2*Rmax/C+Tr-2*Rmin/C)/Nfast;        %更新Fsr、dt
Fsr=1/dt;

%%参数——分辨率
DY=C/2/Br;                          %距离向分辨率:DY=C/2/Br;
DX=D/2;                             %方位向分辨率:DX=D/2;

%%参数——点目标
Ntarget=1;                          %点目标个数

%格式[方位向, 距离向水平投影,高度，目标散射系数]
Z=-50;
Ptarget=[15*DY,Y+20*DX,Z,1];

%%========================================================
%%产生单、多点回拨回波信号，一个N*M矩阵
K=Ntarget;                       
N=Nslow;                         
M=Nfast;                        
T=Ptarget;                      
Srnm1=zeros(N,M);
for k=1:1:K                      
    sigma=T(k,4);                 %获取点目标散射系数
    Dslow=sn*V-T(k,1);            %慢时域中雷达与目标方位向距离，1*N矩阵
    R=sqrt(Dslow.^2+T(k,2)^2+(H-T(k,3))^2);%斜距，1*N矩阵
    tau=2*R/C;                    %时间延迟，1*N矩阵
    Dfast=ones(N,1)*tm-tau'*ones(1,M);                  %[t-2*R/C]项，N*M矩阵
    phase=pi*Kr*Dfast.^2-(4*pi/lambda)*(R'*ones(1,M));  %回波总相位，N*M矩阵
    Srnm1=Srnm1+sigma*exp(j*phase).*(0<Dfast&Dfast<Tr).*((abs(Dslow)<Lsar/2)'*ones(1,M));%回波
end
%%========================================================
%%========================================================
%%距离向压缩
tr=tm-2*Rc/C-Tr/2;                                 %快时域时间-最小延迟时间，1*M矩阵
Refr=exp(j*pi*Kr*tr.^2).*(0<tr&tr<Tr);             %距离压缩参考函数
Sr1=ifty(fty(Srnm1).*(ones(N,1)*conj(fty(Refr))));   %距离压缩后回波信号Sr
Gr1=abs(Sr1);                                        %Sr包络
%%========================================================
%%方位向压缩
ta=sn;                                             %慢时域时间，1*N矩阵
Refa=exp(j*pi*Ka*ta.^2).*(abs(ta)<Tsar/2);         %方位压缩参考函数
Sa1=iftx(ftx(Sr1).*(conj(ftx(Refa)).'*ones(1,M)));   %方位压缩后回波信号Sa（距离向压缩基础上）
Ga1=abs(Sa1);                                        %Sa包络
%%========================================================


%定位与显示
dr1=dt*C;
da1=ds*V;
[m,n]=find(Ga1==max(max(Ga1)));
r1=Rc+(n-668)*dr1/2;                                           
a1=(m-513)*da1;

disp('相对S1距离向位置：');disp(r1)
disp('相对S1方位向位置：');disp(a1)




%%参数——目标区域
Y=2000;                          %测绘中心Y坐标
Y2=Y;
%%参数
theta=0;                         %平台波束斜视角
theta=theta*pi/180;              %rad
Rb=sqrt(H^2+Y^2);                %区域中心斜距
Rc=Rb/cos(theta);                %目标区域距离向范围：[Rc-Y0,Rc+Y0]

%%参数——雷达天线
Lsar=lambda*Rc/D;                %SAR合成孔径长度
Tsar=Lsar/V;                     %SAR合成孔径时间

%%参数——慢时间域
Ka=-2*(V*cos(theta))^2/lambda/Rc;            %doppler调频斜率
Ba=abs(Ka*Tsar);                %doppler调频频宽
PRF=Ba;                         %脉冲重复频率
PRT=1/PRF;                      %脉冲重复时间
ds=PRT;                         %慢时域采样间距
                                %为了进行FFT，对PRT、PRF、ds进行运算并更新
Nslow=ceil((2*X0+Lsar)/V/ds);   %慢时域采样点数
Nslow=2^nextpow2(Nslow);        %使采样点数为2的幂；N点

%Nslow=1024;

sn=linspace(-(X0+Lsar/2)/V,(X0+Lsar/2)/V,Nslow);%慢时域离散时间序列，1*N矩阵
PRT=(2*X0+Lsar)/V/Nslow;        %更新PRT、PRF、ds
PRF=1/PRT;
ds=PRT;
%%参数——快时域
Fsr=2*Br;                       %快时域采样频率
dt=1/Fsr;                       %快时域采样间距
Rmin=sqrt((Y-Y0)^2+(H-Z0)^2);        %最小斜距
Rmax=sqrt((Y+Y0)^2+(H+Z0)^2+(Lsar/2)^2);     %最大斜距
Nfast=ceil(2*(Rmax-Rmin)/C/dt+Tr/dt);   %快时域采样点数
Nfast=2^nextpow2(Nfast);                %使采样点数为2的幂；M点
tm=linspace(2*Rmin/C,2*Rmax/C+Tr,Nfast);%快时域离散时间序列，1*M矩阵
dt=(2*Rmax/C+Tr-2*Rmin/C)/Nfast;        %更新Fsr、dt
Fsr=1/dt;

%格式[方位向, 距离向水平投影,高度，目标散射系数]
Ptarget=[15*DY,Y+20*DX,Z,1];
%%========================================================
%%产生单、多点回拨回波信号，一个N*M矩阵
K=Ntarget;                       
N=Nslow;                         
M=Nfast;                        
T=Ptarget;                      
Srnm2=zeros(N,M);
for k=1:1:K                      
    sigma=T(k,4);                 %获取点目标散射系数
    Dslow=sn*V-T(k,1);            %慢时域中雷达与目标方位向距离，1*N矩阵
    R=sqrt(Dslow.^2+T(k,2)^2+(H-T(k,3))^2);%斜距，1*N矩阵
    tau=2*R/C;                    %时间延迟，1*N矩阵
    Dfast=ones(N,1)*tm-tau'*ones(1,M);                  %[t-2*R/C]项，N*M矩阵
    phase=pi*Kr*Dfast.^2-(4*pi/lambda)*(R'*ones(1,M));  %回波总相位，N*M矩阵
    Srnm2=Srnm2+sigma*exp(j*phase).*(0<Dfast&Dfast<Tr).*((abs(Dslow)<Lsar/2)'*ones(1,M));%回波
end
%%========================================================
%%========================================================
%%距离向压缩
tr=tm-2*Rc/C-Tr/2;                                 %快时域时间-最小延迟时间，1*M矩阵
Refr=exp(j*pi*Kr*tr.^2).*(0<tr&tr<Tr);             %距离压缩参考函数
Sr2=ifty(fty(Srnm2).*(ones(N,1)*conj(fty(Refr))));   %距离压缩后回波信号Sr
Gr2=abs(Sr2);                                        %Sr包络
%%========================================================
%%方位向压缩
ta=sn;                                             %慢时域时间，1*N矩阵
Refa=exp(j*pi*Ka*ta.^2).*(abs(ta)<Tsar/2);         %方位压缩参考函数
Sa2=iftx(ftx(Sr2).*(conj(ftx(Refa)).'*ones(1,M)));   %方位压缩后回波信号Sa（距离向压缩基础上）
Ga2=abs(Sa2);                                        %Sa包络
%%========================================================


%定位与显示
dr2=dt*C;
da2=ds*V;
[m,n]=find(Ga2==max(max(Ga2)));
r2=Rc+(n-658)*dr2/2;                                
a2=(m-513)*da2;

disp('相对S2距离向位置：');disp(r2)
disp('相对S2方位向位置：');disp(a2)

Rt=6371300;             %地球平均半径
R1=H+Rt;             %雷达平台地心距
R2=H+Rt;
dY=Y1-Y2;
dR=r1-r2;
%costheta1=(r1^2+dY^2-r2^2)/(2*r1*dY);
%costheta2=dR/dY-dR^2/2/r1/dY+dY/2/r1;
thetavt=acos((r1^2+dY^2-r2^2)/(2*r1*dY));
thetav=pi/2-thetavt;

h=sqrt(r1^2+R1^2-2*r1*R1*cos(thetav))-Rt;

disp('高度位置：');disp(h)