论文方法复现：A High Frequency Vibration Compensation Approach For Terahertz SAR Based on DSFMFT

---

1、概述

该文章是将振动误差的相位展开为傅里叶级数形式，并利用正弦调频傅里叶变换来估计相位傅里叶级数展开的系数，并根据所估计的系数进行相位补偿。
根据论文思路所设计的仿真算法流程为：
（1）对回波信号进行距离向压缩，RCMC和去斜操作
（2）提取回波数据中单个独立强散射点信号（由于所设计的仿真程序是单个点目标仿真，因此在仿真程序中将这一步简化了）
（3）提取单个独立强散射点相位并进行解缠，然后再利用最小二乘法去除一次相位项，所剩下的即为振动相位项
（4）进行正弦调频傅里叶变换，所得到的系数即为振动相位项按照傅里叶级数展开的系数
（5）根据估计得到的系数计算回波数据中所需补偿的振动相位（这一步在仿真程序中进行了简化）
（6）对补偿后的数据进行正常成像

2、仿真

关于太赫兹SAR振动分析的内容可参考太赫兹SAR振动分析_afgc223的博客-CSDN博客
对单个点目标进行仿真，设置振动斜距为

r_v=0.001*sin(2*pi*60*tm)+0.001*cos(2*pi*20*tm);

在进行完成步骤（2）后，步骤（3）中信号所提取的解缠后的相位为

利用去除一次相位项后，得到

而所设置的理论上的振动相位为

可见两者除了初始相位不同外，并没有太大的差别。

最后根据估计的系数进行补偿，所得到的补偿前（左）后（右）结果分别为

3、代码

生成回波数据的程序SAR_ECHO_Vibration.m ，以及sinc_interplotion.m代码参见https://blog.csdn.net/afgc223/article/details/133776470

主程序 SAR_RD.m

close all; clear; clc;
load('suv_point_sin.mat'); %加载回波数据                                                            %加载仿真数据
%% 参数设置
c = 2.9979e8;                                                               %光速  
fc=3e11;
Br=28.8e9;                                                                  %带宽
Vr=5;                                                                       %平台运动速度
PRT=pi*0.148/720/Vr;                                                        %方位向脉冲重复时间
theta=6.5/180*pi;                                                          
%% 计算相关参数
lambda = c/fc;                                                              %波长
[Nslow,Nfast] = size(suv);                                                  %慢时间，快时间个数
%deltaf=Br/(Nfast-1);                                                       %步进频率值
deltaR=c/2/Br;                                                              %步进距离
%R=0:deltaR:((Nfast-1)*deltaR);                                              %方位向距离
R=0:deltaR:((Nfast-1)*deltaR);
fa=1/PRT;                                                                   %方位向采样率
f =linspace(-fa/2,fa/2,Nslow).';                                            %多普勒频率分布
%x=((0:PRT:((Nslow-1)*PRT))*Vr).';                                           %方位向分布
x=linspace(-(Nslow-1)*PRT/2*Vr,(Nslow-1)*PRT/2*Vr,Nslow).';
%Fam=2*Vr/lambda*sin(theta/2);                                               %最大多普勒频率
%% rd算法处理
ECHO=THz_SAR_SFMFT(fc,Br,Vr,PRT,suv);
%% 绘图
Img=abs(ECHO)/max(max(abs(ECHO)));
figure,shading interp,colormap(jet);
imagesc(R,x,Img);axis image ,xlim([3.8,4.2]),ylim([-0.2,0.2]);
title('RD算法处理'),xlabel('距离(m)'),ylabel('方位(m)'),view(270,90);
set(gcf,'position',[0,0,1280,1280]);

实现振动补偿的函数 THz_SAR_SFMFT.m

function [ECHO] = THz_SAR_SFMFT(fc,Br,Vr,PRT,suv)
% 用于复现论文A High Frequency Vibration Compensation Approach For Terahertz SAR Based on 
% Sinusoidal Frequency Modulation Fourier Transform
%f0步进频率雷达中心频率，Br带宽，Vr平台速度，PRT方位向采样率，suv SAR回波信号
c = 2.9979e8;                                                               %光速  
%% 计算相关参数
lambda = c/fc;                                                              %波长
[Nslow,Nfast] = size(suv);                                                  %慢时间，快时间个数
%deltaf=Br/(Nfast-1);
deltaR=c/2/Br;                                                              %步进距离
R=0:deltaR:((Nfast-1)*deltaR);                                              %方位向距离
fa=1/PRT;                                                                   %方位向采样率
f =linspace(-fa/2,fa/2,Nslow).';                                            %多普勒频率分布
%x=(0:1:Nslow-1)*PRT*Vr;                                                   %计算方位向位置
%% 距离向逆傅里叶变换以及方位向傅里叶变换
ECHO=ifft(suv.').';
range=ones(Nslow,1)*exp(-1j*pi*(Nfast-1)*(0:Nfast-1)/Nfast);
ECHO=ECHO.*range;
suv=fftshift(fft(ECHO),1);
%% 距离徙动校正
N=6;                           %插值核长度
ECHO_RCMC=zeros([Nslow,Nfast]);
%抛物线近似模型
% RCM = (lambda*f/Vr).^2*R/8;    %距离徙动量
% RCMC=(ones(Nslow,1)*R+RCM)/deltaR+1;  
%精确的模型
Rrm=ones(Nslow,1)*R;
fam=f*ones(1,Nfast);
RCMC=(Rrm./sqrt(1-fam.^2*lambda^2/4/Vr^2))/deltaR+1;
h = waitbar(0,'插值中......');  %生成一个进度条
for i=1:Nslow 
    ECHO_RCMC(i,:)=sinc_interplotion(suv(i,:),RCMC(i,:),N);%  sinc插值
%     aai=(RCMC(i,:)<(Nfast-1))&(RCMC(i,:)>=0);
%     aa=aai.*RCMC(i,:)+(1-aai).*(Nfast-1);
%     ECHO_RCMC(i,:)=interp1(1:(Nfast),suv(i,:),aa,'linear');%线性插值
    waitbar(i/Nslow);
end
close(h);  %关闭进度条
suv=ifft(ifftshift(ECHO_RCMC,1));
%% 去斜处理与信号提取
tm=((0:1:Nslow-1)*PRT).';      %慢时间
ECHO_RCMC=suv.*exp(tm.^2*(1j*2*pi*fc/c*Vr^2./R));%.*exp(-1j*tm*(4*pi*fc/c*Vr^2*tm(end)/2./R));%去斜
[~,I_y]=find(abs(ECHO_RCMC)==max(max(abs(ECHO_RCMC))));
Sig=ECHO_RCMC(:,I_y);%提取信号
%% SFMFT
figure;plot(tm,phase(Sig));
P_Sig=phase(Sig);
%P_Sig = P_Sig - polyval(polyfit([1 Nslow],P_Sig([1 Nslow]),1),1:Nslow);
P_Sig = P_Sig - polyval(polyfit(tm,P_Sig,1),tm);
P_Sig=P_Sig-mean(P_Sig);
figure;plot(tm,P_Sig);
f_P_Sig=fftshift(fft(P_Sig));
figure;plot(f,abs(f_P_Sig));
Phi=exp(-1j*ifft(ifftshift(f_P_Sig)))*ones(1,Nfast);
suv=suv.*Phi;
ECHO_RCMC=fftshift(fft(suv),1);
%% 方位向匹配滤波
Ka = 2*Vr^2./(R*lambda); %多普勒调频率
for k=1:Nfast
    ECHO(:,k) = ifft(exp(-1j*pi*f.^2/Ka(k)).*ECHO_RCMC(:,k)); %方位向压缩后
end
end

3、总结

（1）虽然是对论文一个简单版的复现，并且没有讨论论文中所提及的相位解缠方面的内容，但基本上实现了论文的思路。
（2）该论文方法最关键的一点是需要找到回波数据中的独立强散射点，才能够使得提取的信号相位不至于存在太大的误差。而当场景中没有独立散射点时，那么该方法误差补偿的效果将不会特别理想。
（3）即使在存在独立散射点的情况下，利用单个散射点来估计整个场景的振动误差相位存在一定误差，因为单个散射点的信息只包含点合成孔径内的振动误差情况，而不能涵盖住整个场景的振动误差信息。

参考文献

S. Shi, C. Li, J. Hu, X. Zhang and G. Fang, “A High Frequency Vibration Compensation Approach For Terahertz SAR Based on Sinusoidal Frequency Modulation Fourier Transform,” in IEEE Sensors Journal, vol. 21, no. 9, pp. 10796-10803, 1 May1, 2021, doi: 10.1109/JSEN.2021.3056519.