SAR成像算法matab程序

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Bi-SAR %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear all; clc; close all; DOPPLER = 1; % DOPPLER == 1则绘制多普勒历史，否则不绘。以下类似。 BEAM = 1; % BEAM == 1 则绘制收发波束对地照射区域 ELLIPSE = 1; % 等距离椭圆 EQUIDOPPLER = 1; % 等多普勒线绘制 CASSINI = 1; % 卡西尼卵形线，即等信噪比线 SIMULATION = 1; % 模拟产生原始数据 DISPLAY = 1; % 显示 TSEA = 1; % Time Synchronization Error Ajusment RD = 1; % RangeDoppler 算法 DEBUG = 0; % 调试 SEEDOPPLER = 0; % 观察多普勒频率 %% 选择参数 SCENARIO = 1; switch (SCENARIO) case 1; % 两平台平行飞行，以接收机为中心，飞行速度一致 % 接收平台飞过场景中心正上方的时刻为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; % 椭球上的点到两焦点的距离和为 2*RanCon. % 绘制等距离线用 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 = 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 = sqrt(sum((T0 - CP).^2))+sqrt(sum((R0 - CP).^2)); lambda = 0.1; % 波长 otherwise error(' There is no such scenario!'); end; %% 计算多普勒历史 delta_t = 0.01; % 引入delta_t主要是为了后面能够正确计算多普勒频率,此delta_t表示的是方位向上的时间（慢时间） t = -15:delta_t:15; % 方位向时间，取得比较长，主要是为了把趋势看得比较清楚些，实际的飞行时间是不会这么长的。 [ColuDem,TimeDem] = size(t); % 获取t 的维数，即RowDem % 点目标 pnum = 1; % 点目标数目 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; 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; 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; 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 Angel Varying with Time '); end; % 对应for(count = 1:pnum) end; % 对应 if (DOPPLER == 1) % %% 绘制波束照射区域 if BEAM==1 % 计算接收波束照射区域 theta_r_ver = 11.68; % 接收波束在垂直向的角度,r代表reciever, ver代表 vertical theta_r_ver = theta_r_ver * pi/180;% 转化为弧度 %H_r = ;R_r = ; %phi_r = asin(H_r/R_r); % 接收波束的俯角,即从波束中心线向上转动直到波束水平时所转过的角度。区别于投射角 phi_r = 11.537; phi_r = phi_r*pi/180; % 转化为弧度 Height_r = R0(3,1); % 以接收平台在0时刻时的值为接收平台的高度，此时默认地面高度为0;如果考虑点目标的高度则还要修改程序； RBCR = Height_r*cot(phi_r+0.5*theta_r_ver); % Recieving Beam Closest Range, 注意：这是我造的词。Range 指地距。 % RBCR 代表了接收平台在地面的投影到接收波束较 近 边缘的距离。 RBFR = Height_r*cot(phi_r-0.5*theta_r_ver); % Recieving Beam Farest Range, 注意：这是我造的词。Range 指地距。 % RBFR 代表了接收平台在地面的投影到接收波束较