双站SAR模拟matlab代码

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Project: ThesisBiSimu_syn.m %% %% %% Author: K Yuan %% %% Date : May 16,2004 %% %% %% Functionality: 分析双基地SAR的几何特性，包括：照射区域，等距离线， %% %% 波束仿真，等多普勒线，卡西尼卵形线； %% %% 产生点目标数据，并使用RD算法进行成象 %% %% 可以模拟触发同步误差，相位同步误差，运动误差 %% %% 备注 : 用于飞行实验验证 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear all; % 清除所有的变量 clc; % 清屏 close all; % 关闭所有的图形窗口 DOPPLER = 1; % DOPPLER == 1则绘制多普勒历史，否则不绘。以下类似。 BEAM = 1; % BEAM == 1 则绘制收发波束对地照射区域 ELLIPSE = 1; % 等距离椭圆 EQUIDOPPLER = 1; % 等多普勒线绘制 CASSINI = 0; % 卡西尼卵形线，即等信噪比线 SIMULATION = 1; % 模拟产生原始数据 DISPLAY = 0; % 显示 TSEA = 1; % Time Synchronization Error Ajusment RD = 1; % RangeDoppler 算法 BP = 0; % BackProjection 算法 DEBUG = 0; % 调试 SEEDOPPLER = 0 % 观察多普勒频率 % % 选择参数 SCENARIO = 1; switch (SCENARIO) case 1; % 两平台平行飞行，以接收机为中心 % 暂时先考虑两机飞行速度一致的情况 % 假定接收平台飞过场景中心正上方的时刻为0 % 此为以接收机为中心的情况。 %Vt = [10.41,59.085,0].'; % 发射平台的飞行速度向量 Vt = [0,60,0].'; % 发射平台的飞行速度向量 T0 = 1e3*[-14.6970,0,6].'; % 发射平台 0 时刻的位置 Vr = [0,60,0].'; % 接收平台的飞行速度向量 R0 = 1e3*[0,0,3].'; % 接收平台 0 时刻的位置 CENTER = 1e3*[ 14.6969,0,0].'; % 场景中心 X= 1e3*[-30:50]; % 确定所绘的椭圆在x轴上的范围,该参数与收发平台的相对位置有较大的关系 [XColDem,XRowDem] = size(X);% 获取X的维数，后面计算的时候要用 RanCon = 60e3; % 椭球上的点到两焦点的距离和为 2RanCon. % 绘制等距离线用 RefRange = 2 * sqrt(sum((R0-CENTER).^2)); ; % 参考距离 Xc = CENTER(1,1); Yc = CENTER(2,1); X0 = 50; % 场景在X方向的测绘宽度的一半 Y0 = 100; % 场景在Y方向的测绘宽度的一半 CP = CENTER + [-X0,0,0].';% 场景中的最近距离点 ！！！注意，此最近距离点只在收发同侧的时候有效，如果收发平台位于测绘带两侧则还不能这样计算 %Rmin = 5000; % 收发平台的最短距离和 Rmin = sqrt(sum((T0 - CP).^2))+sqrt(sum((R0 - CP).^2)); lambda = 0.1; % 波长 case 2; % 非平行飞行，以接收机为中心 Vt = [20,sqrt(3600-400-25),5].'; % 发射平台的飞行速度向量 T0 = 1e3*[-14.6970,2,6].'; % 发射平台 0 时刻的位置 Vr = [0,60,0].'; % 接收平台的飞行速度向量 R0 = 1e3*[0,0,3].'; % 接收平台 0 时刻的位置 CENTER = 1e3*[ 14.6969,0,0].'; % 场景中心 X= 1e3*[-30:50]; % 确定所绘的椭圆在x轴上的范围,该参数与收发平台的相对位置有较大的关系 [XColDem,XRowDem] = size(X);% 获取X的维数，后面计算的时候要用 RanCon = 60e3; % 椭球上的点到两焦点的距离和为 2RanCon. % 绘制等距离线用 RefRange = 2 * sqrt(sum((R0-CENTER).^2)); ; % 参考距离 Xc = CENTER(1,1); Yc = CENTER(2,1); X0 = 500; % 场景在X方向的测绘宽度的一半 Y0 = 1000; % 场景在Y方向的测绘宽度的一半 CP = CENTER + [-X0,0,0].';% 场景中的最近距离点 ！！！注意，此最近距离点只在收发同侧的时候有效，如果收发平台位于测绘带两侧则还不能这样计算 %Rmin = 5000; % 收发平台的最短距离和 Rmin = sqrt(sum((T0 - CP).^2))+sqrt(sum((R0 - CP).^2)); lambda = 0.1; % 波长 otherwise error(' There is no such scenario!'); end; % %% 计算多普勒历史 delta_t = 0.001; % 引入delta_t主要是为了后面能够正确计算多普勒频率,此delta_t表示的是方位向上的时间（慢时间） t = [-15:delta_t:15]; % 方位向时间，取得比较长，主要是为了把趋势看得比较清楚些，实际的飞行时间是不会这么长的。 [ColuDem,TimeDem] = size(t); % 获取t 的维数，即RowDem % 点目标 pnum = 2; % 点目标数目 P = zeros(3,pnum); P(:,1) = CENTER ; P(:,2) = CENTER + 1e3*[2,6,2].'; P(:,3) = CENTER + 1e3*[0,-0.1,0].'; P(:,4) = CENTER + 1e3*[0.1,0,0].'; P(:,5) = CENTER + 1e3*[-0.1,0,0].'; % 产生运动轨迹 %MotionErrTra = zeros(3,TimeDem); % 发射平台运动误差 %MotionErrRec = zeros(3,TimeDem); % 接收平台运动误差 %% 注意:现在决定把运动误差加到后面Simulation部分来进行观察 T = T0*ones(1,TimeDem) + Vt*t ; % 发射机平台的位置，为空间中的一点，Rt0为其零时刻时的位置， % Vt为发射机平台的飞行方向。 R = R0*ones(1,TimeDem) + Vr*t ; % 接收机平台的位置，为空间中的一点，Rr0为其零时刻时的位置， % Vr为发射机平台的飞行方向。 Rt_abs = zeros(pnum,TimeDem); % 发射机相位中心到点目标P的距离向量 Rr_abs = zeros(pnum,TimeDem); % % 接收机相位中心到点目标P的距离向量 RanSum = zeros(pnum,TimeDem); % 点目标到收发平台的距离和 DopFre_t = zeros(pnum,TimeDem); % 由发射平台运动造成的多普勒频率 DopFre_r = zeros(pnum,TimeDem); % 由发射平台运动造成的多普勒频率 DopFre = zeros(pnum,TimeDem); % 由收发平台的共同运动造成的多普勒频率 DiffDopFre = zeros(pnum,TimeDem); % 对多普勒频率求导，即多普勒斜率，主要是想看一下多普勒频率的线性和非线性区 beta = zeros(pnum,TimeDem); %双基地角 theta_t = zeros(pnum,TimeDem); % 发射平台飞行方向与Rt间的夹角 theta_r = zeros(pnum,TimeDem); % 接收平台飞行方向与Rr间的夹角 for (count = 1:pnum) Rt = P(:,count)*ones(1,TimeDem) - T; % 发射机相位中心到点目标P的距离向量 Rr = P(:,count)*ones(1,TimeDem) - R; % 接收机相位中心到点目标P的距离向量 Rt_abs(count,:) = (sum(Rt.^2)).^0.5; % 点目标到发射平台的距离 Rr_abs(count,:) = (sum(Rr.^2)).^0.5; % 点目标到接收平台的距离 RanSum(count,:) = (sum(Rt.^2)).^0.5 + (sum(Rr.^2)).^0.5 ; % 点目标到收发平台的距离和 DopFre_t(count,:) = [NaN,-(diff(Rt_abs(count,:))/delta_t)/lambda]; % 由发射平台运动造成的多普勒频率 DopFre_r(count,:) = [NaN,-(diff(Rr_abs(count,:))/delta_t)/lambda]; % 由接收平台运动造成的多普勒频率 DopFre(count,:) = [NaN,-(diff(RanSum(count,:))/delta_t)/lambda]; % 由收发平台的共同运动造成的多普勒频率 % 经diff后DopFre_t的维数比t少了一维，所以要补一个NaN回来, % 否则下面在显示的时候会遇到问题 DiffDopFre(count,:) = [NaN,diff(DopFre(count,:))/delta_t]; % 对多普勒频率求导，即多普勒斜率，主要是想看一下多普勒频率的线性和非线性区 beta(count,:) = acos((sum(Rt.*Rr))./(Rt_abs(count,:).*Rr_abs(count,:))); %双基地角 beta(count,:) = beta(count,:)*180/pi; % 转换为角度 Vt_abs = (sum(Vt.^2)).^0.5; theta_t(count,:) = acos((sum((Vt*ones(1,TimeDem)).*Rt))./(Vt_abs.*Rt_abs(count,:))); % 发射平台飞行方向与Rt间的夹角 theta_t(count,:) = theta_t(count,:)*180/pi;% 转换为角度 Vr_abs = (sum(Vr.^2)).^0.5; theta_r(count,:) = acos((sum((Vr*ones(1,TimeDem)).*Rr))./(Vr_abs.*Rr_abs(count,:))); % 接收平台飞行方向与Rr间的夹角 theta_r(count,:) = theta_r(count,:)*180/pi;% 转换为角度 end; % 对应于for (count = 1:pnum) if (DOPPLER == 1) for(count = 1:pnum) figure(1); hold on; plot(t,Rt_abs(count,:),'red'); plot(t,Rr_abs(count,:),':green'); plot(t,RanSum(count,:),'-.blue'); hold off; xlabel('Azimuth Time,Second'); ylabel('Range,Meters'); title('Range Sum '); figure(2); hold on; plot(t,DopFre_t(count,:),'-red'); plot(t,DopFre_r(count,:),':green'); plot(t,DopFre(count,:),'-.blue'); plot(t,DiffDopFre(count,:),':red'); hold off; xlabel('Azimuth Time,Second'); ylabel('Doppler Frequency,Hz'); title('Doppler History'); figure(3); hold on; plot(t,beta(count,:),'red'); plot(t,theta_t(count,:),':green'); plot(t,theta_r(count,:),'-.blue'); hold off; xlabel('Azimuth Time,Second'); ylabel('Bistatic Angel,Degree'); title('Bistatic/Theta Ang