提示：文章写完后，目录可以自动生成，如何生成可参考右边的帮助文档
论文方法复现：Compensation for High-Frequency Vibration of Platform in SAR Imaging Based on Adaptive Chirplet Decomposition
1、Chirplet 分解
该论文是以Chirplet 分解为基础实现太赫兹SAR振动补偿的参数估计，Chirplet分解是将一个信号分解为多个线性调频信号叠加的形式，即
其中
为高斯包络的线性调频信号。
根据论文的思路设计如下Chirplet 分解算法：
(1)设置初始固定窗函数长度，寻找信号Chirplet变换的幅值最大的时刻和调频率。
(2)根据选定的调频率和信号时刻计算使得Chirplet变换幅值最大的窗函数长度。
(3)选择窗函数形状使得信号与具有窗函数包络的chirp信号内积最大
(4)从原始信号中减去分解得到的chirp信号，返回（1）继续迭代，直到达到所设定的迭代次数。
下面进行信号Chirplet 分解的仿真，所设置的仿真信号为
s_t=exp(1j*pi*4*sin(B*t)).*(2*exp(1j*2*pi*100*t)+3*exp(-1j*2*pi*300*t));
设置初始固定窗函数长度为原始信号长度的1/10，迭代次数为14，仿真信号的时频关系图为
随着迭代进行，原始信号不断被分解为chirp信号，分解的过程如下
代码如下：
主程序 main.m
clc;clear;close all
%% 参数设置
B=100;
T=0.1;
K=-B/T;
fs=2200;
t=(0:1/fs:T);
%% 信号采样
%s_t=exp(1j*pi*K*t.^2);
%s_t=exp(1j*pi*K*t.^2).*(2*exp(1j*2*pi*100*t)+3*exp(-1j*2*pi*300*t));
%s_t=exp(1j*pi*4*sin(B*t)); /length(Index);/length(Index);/length(Index);/length(Index);/length(Index);/length(Index);/length(Index);/length(Index);/length(Index);/length(Index);/length(Index);/length(Index);/length(Index);/length(Index);
s_t=exp(1j*pi*4*sin(B*t)).*(2*exp(1j*2*pi*100*t)+3*exp(-1j*2*pi*300*t));
%% 自适应Chirp分解
Adaptive_Chirplet_Decomposition(s_t,t,fs,round(length(s_t)/10),14)
用于实现chirplet分解的函数 Adaptive_Chirplet_Decomposition.m
function [A,f] = Adaptive_Chirplet_Decomposition(Sig,t,SampFreq,d0,N)
%ADAPTIVE_CHIRPLET_DECOMPOSITION 用于实现论文Compensation for High-Frequency Vibration of Platform in SAR
%Imaging Based on Adaptive Chirplet Decomposition 信号的Chirplet分解
% Sig 被分解的原始信号
% t 时间轴
% SampFreq 信号采样率
% d0 初始设置的高斯窗函数宽度
% N 算法迭代次数，迭代完成后，算法被分解为N个chirp信号
% f 估计出的振动频率
% A 估计出的振动幅度
t_c=zeros(1,N);%储存中心时间数据
alpha_c=zeros(N,1);%储存调频率
for j=1:N
N_S=length(Sig);%信号长度
% 首先估计最佳调频率alpha_optimal和时间频率轴位置I_t,I_f
N_alpha=1000;%设置调频率扫描个数
alpha=linspace(-SampFreq*2/(t(end)-t(1)),SampFreq*2/(t(end)-t(1)),N_alpha);%设置调频率扫描范围
maximum=0;
index=0;
h = waitbar(0,'调频率扫描中......'); %生成一个进度条
for i=1:N_alpha
[Spec,Freq,~] = Chirplet_Transform(Sig,length(Sig),d0,SampFreq,alpha(i));
b=max(max(abs(Spec)));
if maximum<b
Max_Spec=Spec;
maximum=b;
index=i;
end
waitbar(i/N_alpha);
end
close(h);
figure;colormap(jet);imagesc(t,Freq,abs(Max_Spec));ylim([-1000,1000]);
[I_f,I_t]=find(abs(Max_Spec)==maximum);
alpha_optimal=alpha(index);
t_c(j)=t(I_t);
alpha_c(j)=alpha_optimal;
% 估计信号最佳持续时间d_optimal
d=5:N_S;%设置持续时间扫描范围
maximum=0;
index=0;
for i=1:length(d)
[Spec,~,~] = Chirplet_Transform_d(Sig,length(Sig),d(i),SampFreq,alpha_optimal,I_t,6);
b=abs(Spec(I_f,I_t));
if maximum<b
maximum=b;
index=i;
end
end
d_optimal=d(index);
%估计信号系数
win_c=linspace(0,8,100);
minimum=norm(Sig);
index=0;
for i=1:length(win_c)
[~,Chirp,iIndex1] = Chirplet_Transform_d(Sig,length(Sig),d_optimal,SampFreq,alpha_optimal,I_t,win_c(i));
chirp_=Chirp.*exp(-1j*2*pi*t(iIndex1)*Freq(I_f));
C=chirp_*Sig(iIndex1).'/(chirp_*chirp_');
b=norm(Sig(iIndex1)-C*conj(chirp_));
if minimum>b
index=i;
minimum=b;
chirp_optimal=chirp_;
C_optimal=C;
end
end
win_c(index)
% [~,Chirp,iIndex1] = Chirplet_Transform_d(Sig,length(Sig),d_optimal,SampFreq,alpha_optimal,I_t,6);
% chirp_=Chirp.*exp(-1j*2*pi*t(iIndex1)*Freq(I_f));
% chirp_optimal=chirp_;
% C_optimal=chirp_*Sig(iIndex1).'/(chirp_*chirp_');
Sig(iIndex1)=Sig(iIndex1)-C_optimal*conj(chirp_optimal);
end
[t_cc,I]=sort(t_c);
f=(N-1)/4/sum(diff(t_cc));
alpha_cc=abs(alpha_c(I));
t_cc=t_cc-t_cc(1);
%A=-diff(t_cc)*abs(alpha_cc(2:N))/(N-1)/2/pi/f*2*pi;
A=-2*t_cc*abs(alpha_cc)/(N)/(N-1)/2/pi/f*2*pi;
end
Chirplet变换 Chirplet_Transform.m 代码参见 https://blog.csdn.net/afgc223/article/details/133710620?spm=1001.2014.3001.5501
Chirplet_Transform_d.m
function [Spec,Chirp,iIndex1] = Chirplet_Transform_d(Sig,fLevel,WinLen,SampFreq,alpha,I_t,win_c)
% --------------------INPUT------------------%
% Sig 一维原始信号 为横向量
% fLevel 频谱长度
% WinLen 高斯窗函数长度
% SampFreq 信号采样频率
% alpha 信号调频率（chirp rate）
% I_t 只对时间轴位置为I_t进行计算
% --------------------OUTPUT------------------%
% Spec Chirplet变换结果
% Chirp 估计的Chirp信号
% t 水平时间轴
%--------------------------------------------------------%
if (nargin < 1)
error('At least one parameter required!');
end
if (nargin < 4) % 默认采样率 1000Hz
SampFreq = 1000;
end
if (nargin < 3) % 默认窗函数长度64
WinLen = 64;
end
if (nargin < 2) % 默认频谱长度512
fLevel = 512;
end
%% data preparation
SigLen = length(Sig); % 信号长度
t = (0:1:SigLen-1)*(1/SampFreq); % 计算时间轴
%% Frequency axis and its half-axis points
fLevel = ceil(fLevel/2) * 2+1; % 将频谱长度转化为最近的奇数
Lf = (fLevel - 1)/2; % 频谱长度的一半
%% Generate Gauss window functions
WinLen = ceil(WinLen/2) * 2+1; % 将窗函数长度转化为最接近的奇数
WinFun = exp(-win_c* linspace(-1,1,WinLen).^2); % 高斯窗函数, 固定时间宽度 [-1, 1], 数据点长度为WinLen
Lw = (WinLen - 1)/2; % 窗函数长度的一半
%% CT spectrum result array (initialized to 0)
Spec = zeros(fLevel,SigLen); %用于存储Chirplet变换结果
%% Sliding window signal,data segmentation preparation
for iLoop = I_t
% Determine the upper and lower limits of the left and right signal index subscripts (note that to prevent the edge width from exceeding the time domain, the retrieval error!)
iLeft = min([iLoop-1, Lw, Lf]);
iRight = min([SigLen-iLoop, Lw, Lf]); %左右两边的滑动窗口均取这些数的最小值，以防止数组越界
iIndex = -iLeft:iRight;%滑动窗口索引值
iIndex1 = iIndex + iLoop; % 原始信号的索引
iIndex2 = iIndex + Lw + 1; % 窗函数的索引
Index = iIndex + Lf +1; % 频率轴索引
R_Operator = exp(-1i * 2*pi*alpha * t(iIndex1).^2 / 2); % 频率旋转算子
S_Operator = exp(1i * 2*pi*alpha * t(iLoop) * t(iIndex1)); % 频率移位算子
Sig_Section = Sig(iIndex1).*R_Operator.*S_Operator; % 信号滑动窗片段乘以旋转和移位算子，正常的Chirplet变换还应乘以相位exp(-1i * 2*pi*alpha * t(iLoop).^2 / 2),
%这里给省略了，因为对于每个给定的iLoop,该相位项为常数项。
Spec(Index, iLoop) = Sig_Section.* conj(WinFun(iIndex2))/norm(WinFun);%/length(Index) % 乘以窗函数
end
%% STFT
Spec = fftshift(fft(Spec),1); % 最后进行STFT
%Freq = linspace(-0.5*SampFreq,0.5*SampFreq,fLevel); % 计算频率轴
Chirp=R_Operator.*S_Operator.*WinFun(iIndex2);
end
2、将Chirplet 分解用于THz SAR单频谐波振动补偿的点目标仿真
关于太赫兹SAR振动分析的内容可参考https://blog.csdn.net/afgc223/article/details/133776470
所设置的单频谐波振动模型为
r_v=0.003*sin(2*pi*20*tm);%单次谐波振动斜距模型
仿真的过程描述为：
（1）根据振动模型生成回波数据，并进行距离压缩和RCMC
（2）在距离和方位向时域进行去斜操作
（3）找到回波数据中散射强度最大的信号进行自适应Chirplet分解
（4）根据分解所得到的信号计算振动谐波的幅度和频率，所依据的是论文中的如下公式
（5）利用估计得到的振动模型进行相位补偿
（6）最后获得振动补偿后的结果
经过去斜后，所得到散射强度最大的信号的时频分析图为
对该信号Chirplet分解的结果为
估计得到振动谐波的振幅与频率分别为A=0.0034m和f=20.4833Hz，这与仿真设置的参数很接近，这里分别给出不进行补偿，根据估计参数进行补偿以及根据真实仿真参数进行补偿的结果，分别如下
根据仿真参数进行补偿的结果
可以发现，即使估计的参数和仿真参数如此接近，进行补偿后的点目标的聚焦效果仍不是特别理想。
代码如下：
生成回波数据的程序SAR_ECHO_Vibration.m ，以及sinc_interplotion.m 代码参见 https://blog.csdn.net/afgc223/article/details/133776470?spm=1001.2014.3001.5501
主程序 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_VC_ACD(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_VC_ACD.m
function [ECHO] = THz_SAR_VC_ACD(fc,Br,Vr,PRT,suv)
% 用于复现论文Compensation for High-Frequency Vibration of Platform in SAR
% Imaging Based on Adaptive Chirplet Decomposition
%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));%去斜
[~,I_y]=find(abs(ECHO_RCMC)==max(max(abs(ECHO_RCMC))));
Sig=ECHO_RCMC(:,I_y).';%提取信号
%% 信号分解
%angle=diff(unwrap(phase(Sig)))/(tm(2)-tm(1))/2/pi/20;
[A,ff]=Adaptive_Chirplet_Decomposition(Sig,tm.',fa,round(length(Sig)/10),7);
% C=sin(pi/3-acos(1/4));
% A=4*pi/lambda*0.003*C;
% ff=20;
% A=0;
% A=-11.4201;
% ff=20.4833;
Phi=exp(1j*A*sin(2*pi*ff*tm))*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
总结
（1）进行Chirp分解所需的计算量较大，特别是对于长度较长的信号
（2）该论文所估计的单次振动谐波受Chirp分解结果的影响而存在一定的误差
（3）要想得到很好的图像聚焦精度对估计参数的精确度要求较高，利用改论文的方法进行太赫兹SAR单频振动补偿很难达到这一点
参考
[1] Y.Wang, Z. Wang, B. Zhao and L. Xu, “Compensation for High-Frequency Vibration of Platform in SAR Imaging Based on Adaptive Chirplet Decomposition,” in IEEE Geoscience and Remote Sensing Letters, vol. 13, no. 6, pp. 792-795, June 2016, doi: 10.1109/LGRS.2016.2544945.
[2] J. C. O’Neill and P. Flandrin, “Chirp hunting,” Proceedings of the IEEE-SP International Symposium on Time-Frequency and Time-Scale Analysis (Cat. No.98TH8380), Pittsburgh, PA, USA, 1998, pp. 425-428, doi: 10.1109/TFSA.1998.721452.
595
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
