基于MATLAB的ERS SAR 原始数据提取和成像

沉迷推公式的猴子

已于 2024-07-18 16:03:17 修改

阅读量611

点赞数 11

文章标签： matlab 开发语言

于 2024-07-18 15:57:51 首次发布

本文链接：https://blog.csdn.net/qq_31482545/article/details/140525209

版权

概述
因为高分三号或者其他国产SAR卫星的原始数据获取的难度非常大。所以我们使用国外的ERS影像进行研究。
本文档阐述了如何从欧洲遥感（ERS）合成孔径雷达（SAR）系统参数和未聚焦的原始数据中提取信息，并使用范围迁移图像形成算法生成聚焦图像。我们将使用NASA的阿拉斯加卫星设施（ASF）发布的ERS数据集进行演示。
数据集下载

使用如下代码从给定的URL下载数据集：

outputFolder = pwd;
dataURL = ['https://ssd.mathworks.com/supportfiles/radar/data/' ...
'ERSData.zip'];
helperDownloadERSData(outputFolder,dataURL);

参数提取
SAR图像形成涉及从参数文件（E2_84699_STD_L0_F299.000.ldr）中提取参数。参数文件包含了形成图像所需的卫星和场景特定数据。像脉冲重复频率、波长、采样率、脉冲宽度、范围门延迟和传感器速度等参数都是从参数文件中提取的。其他参数如有效速度、最小距离、脉冲带宽和脉冲频率是通过从参数文件中提取的数据进行估计得到的。在这个例子中，我们使用辅助函数ERSParameterExtractor从输入文件 E2_84699_STD_L0_F299.000.ldr中提取系统参数。

c = physconst('LightSpeed');               
% Extract ERS system parameters
[fs,fc,prf,tau,bw,v,ro,fd] = ERSParameterExtractor('E2_84699_STD_L0_F299.000.ldr');

原始数据提取
原始数据从数据文件（E2_84699_STD_L0_F299.000.raw）中提取，该文件总共有28652行，每行11644字节。每行中，头部包含412字节，其余11232字节（5616复数）是每个回波的数据。对于ERS雷达，合成孔径长度为1296点。因此，辅助函数ERSDataExtractor读取2048行，因为这是2的幂。使用80%的波束宽度，辅助函数计算有效的方位点并取重叠的块。在这个例子中，我们使用辅助函数从输入文件E2_84699_STD_L0_F299.000.raw以及系统参数中提取未聚焦的原始数据。

% Extract raw data 
rawData = ERSDataExtractor('E2_84699_STD_L0_F299.000.raw',fs,fc,prf,tau,v,ro,fd).';

SAR图像形成
下一步是使用提取的系统参数和未聚焦的原始数据生成单视复杂度（SLC）图像。对于ERS-radar，雷达发射的是线性频率调制脉冲。phased.LinearFMWaveform系统对象使用提取的系统参数创建线性FM脉冲波形。斜向角使用接近于ERS系统的零的Doppler频率计算。范围迁移图像形成算法是一种频域算法，也被称为波数域处理方法。使用rangeMigrationLFM函数将原始数据聚焦到SLC图像。

% Create LFM waveform
waveform = phased.LinearFMWaveform('SampleRate',fs,'PRF',prf,'PulseWidth',tau,'SweepBandwidth',bw,'SweepInterval','Symmetric');
sqang = asind((c*fd)/(2*fc*v));        % Squint angle

% Range migration algorithm
slcimg = rangeMigrationLFM(rawData,waveform,fc,v,ro,'SquintAngle',sqang);

% Display image
figure(1)
imagesc(log(abs(slcimg)))
axis image
colormap('gray')
title('SLC Image')
ylabel('Range bin')
xlabel('Azimuth bin')

在这里插入图片描述
多视处理
通过进行多视处理来减少斑点的效果，从而在图像分辨率上做出权衡。在范围或方位方向，或者在两个方向上，后续的线路被平均以获得更好的图像并减少斑点。多视可以在时域或者频域进行，具体想了解的可以留言。

在这里插入图片描述

辅助函数

ERSParameterExtractor

function [fs,fc,prf,tau,bw,veff,ro,fdop] = ERSParameterExtractor(file)
% Open the parameter file to extract required parameters
fid = fopen(file,'r');

% Radar wavelength (satellite specific)
status = fseek(fid,720+500,'bof');
lambda = str2double(fread(fid,[1 16],'*char'));         % Wavelength (m)

% Pulse Repetition Frequency (satellite specific)
status = fseek(fid,720+934,'bof')|status;
prf = str2double(fread(fid,[1 16],'*char'));            % PRF (Hz)

% Range sampling rate (satellite specific)
status = fseek(fid,720+710,'bof')|status;
fs =str2double(fread(fid,[1 16],'*char'))*1e+06;        % Sampling Rate (Hz)

% Range Pulse length (satellite specific)
status = fseek(fid,720+742,'bof')|status;
tau = str2double(fread(fid,[1 16],'*char'))*1e-06;      % Pulse Width (sec)

% Range Gate Delay to first range cell
status = fseek(fid,720+1766,'bof')|status;
rangeGateDelay = str2double(fread(fid,[1 16],'*char'))*1e-03;   % Range Gate Delay (sec)

% Velocity X
status = fseek(fid,720+1886+452,'bof')|status;
xVelocity = str2double(fread(fid,[1 22],'*char'));    % xVelocity (m/sec)

% Velocity Y
status = fseek(fid,720+1886+474,'bof')|status;
yVelocity = str2double(fread(fid,[1 22],'*char'));    % yVelocity (m/sec)

% Velocity Z
status = fseek(fid,720+1886+496,'bof')|status;
zVelocity = str2double(fread(fid,[1 22],'*char'));    % zVelocity (m/sec)
fclose(fid);

% Checking for any file error
if(status==1)
    fs = NaN;
    fc = NaN;
    prf = NaN;
    tau = NaN;
    bw = NaN;
    veff = NaN;
    ro = NaN;
    fdop = NaN;
    return;
end

% Values specific to ERS satellites
slope = 4.19e+11;           % Slope of the transmitted chirp (Hz/s)
h = 790000;                 % Platform altitude above ground (m)
fdop = -1.349748e+02;       % Doppler frequency (Hz)

% Additional Parameters
Re = 6378144 ;              % Earth radius (m)

% Min distance
ro = time2range(rangeGateDelay);  % Min distance (m)  

% Effective velocity
v = sqrt(xVelocity^2+yVelocity^2+zVelocity^2);
veff = v*sqrt(Re/(Re+h));   % Effective velocity (m/sec)

% Chirp frequency
fc = wavelen2freq(lambda);  % Chirp frequency (Hz)     

% Chirp bandwidth
bw = slope*tau;             % Chirp bandwidth (Hz)
end
ERSDataExtractor

function rawData = ERSDataExtractor(datafile,fs,fc,prf,tau,v,ro,doppler)
c = physconst('LightSpeed');                    % Speed of light

% Values specific to data file
totlines = 28652;                                % Total number of lines
numLines = 2048;                                 % Number of lines
numBytes = 11644;                                % Number of bytes of data
numHdr = 412;                                    % Header size
nValid = (numBytes-numHdr)/2 - round(tau*fs);    % Valid range samples

% Antenna length specific to ERS
L = 10;

% Calculate valid azimuth points
range = ro + (0:nValid-1) * (c/(2*fs));             % Computes range perpendicular to azimuth direction 
rdc = range/sqrt(1-(c*doppler/(fc*(2*v))^2));       % Squinted range 
azBeamwidth = rdc * (c/(fc*L)) * 0.8;               % Use only 80%  
azTau = azBeamwidth / v;                            % Azimuth pulse length 
nPtsAz = ceil(azTau(end) * prf);                    % Use the far range value
validAzPts = numLines - nPtsAz ;                    % Valid azimuth points  

% Start extracting
fid = fopen(datafile,'r');
status = fseek(fid,numBytes,'bof');     % Skipping first line
numPatch = floor(totlines/validAzPts);  % Number of patches           


if(status==-1)
   rawData = NaN;
   return;
end
rawData=zeros(numPatch*validAzPts,nValid);
% Patch data extraction starts
for patchi = 1:numPatch      
    fseek(fid,11644,'cof');
    data = fread(fid,[numBytes,numLines],'uint8')'; % Reading raw data file
    
% Interpret as complex values and remove mean
data = complex(data(:,numHdr+1:2:end),data(:,numHdr+2:2:end));
data = data - mean(data(:));

rawData((1:validAzPts)+((patchi-1)*validAzPts),:) = data(1:validAzPts,1:nValid);
fseek(fid,numBytes + numBytes*validAzPts*patchi,'bof');
end
fclose(fid);
end
multilookProcessing

function image = multilookProcessing(slcimg,sx,sy)
[nx,ny] = size(slcimg);
nfx = floor(nx/sx);
nfy = floor(ny/sy);
image = (zeros(nfx,nfy));
for i=1:nfx
    for j = 1:nfy
        fimg=0;
        for ix = 1:sx
            for jy = 1:sy
                fimg = fimg+slcimg(((i-1)*sx)+ix,((j-1)*sy)+jy);
            end
        end
        image(i,j) = fimg/(sx*sy);
    end
end
end
helperDownloadERSData

function helperDownloadERSData(outputFolder,DataURL)
% Download the data set from the given URL to the output folder.

radarDataZipFile = fullfile(outputFolder,'ERSData.zip');

if ~exist(radarDataZipFile,'file')
    
    disp('Downloading ERS data (134 MiB)...');
    websave(radarDataZipFile,DataURL);
    unzip(radarDataZipFile,outputFolder);
end
end
主函数

c = physconst('LightSpeed');                % Speed of light

% Extract ERS system parameters
[fs,fc,prf,tau,bw,v,ro,fd] = ERSParameterExtractor('E2_84699_STD_L0_F299.000.ldr');

% Extract raw data 
rawData = ERSDataExtractor('E2_84699_STD_L0_F299.000.raw',fs,fc,prf,tau,v,ro,fd).';


% Create LFM waveform
waveform = phased.LinearFMWaveform('SampleRate',fs,'PRF',prf,'PulseWidth',tau,'SweepBandwidth',bw,'SweepInterval','Symmetric');
sqang = asind((c*fd)/(2*fc*v));        % Squint angle

% Range migration algorithm
slcimg = rangeMigrationLFM(rawData,waveform,fc,v,ro,'SquintAngle',sqang);

% Display image
figure(1)
imagesc(log(abs(slcimg)))
axis image
colormap('gray')
title('SLC Image')
ylabel('Range bin')
xlabel('Azimuth bin')


mlimg = multilookProcessing(abs(slcimg),4,20);

% Display Image
figure(2)
imagesc(log(abs(mlimg(1:end-500,:))))
axis image
colormap('gray')
title('Multi-look Image')
ylabel('Range bin')
xlabel('Azimuth bin')

以上代码适用于是2021a及以后的matlab版本。

注意：由于以上matlab代码是一次性把数据读到内存中，可能会出现内存不足的报错，所以你的电脑的内存要尽可能的大。