遥感 SAR 遥感 SAR

【代码】ERS SAR 原始数据提取和成像

ytkz2024-04-092024-09-09

概述

因为高分三号或者其他国产SAR卫星的原始数据获取的难度非常大。所以我们使用国外的ERS影像进行研究。

本文档阐述了如何从欧洲遥感（ERS）合成孔径雷达（SAR）系统参数和未聚焦的原始数据中提取信息，并使用范围迁移图像形成算法生成聚焦图像。我们将使用NASA的阿拉斯加卫星设施（ASF）发布的ERS数据集进行演示。

以下代码均为matlab代码。

数据集下载

该示例使用的ERS数据集包含地球观测卫星委员会（CEOS）标准文件，包括：ldr文件（.ldr）、零级帧数据文件（.raw）、卷描述符文件（.vol）、处理信息文件（.pi）、元数据文件（.meta）和空文件（.nul）。

ldr文件包含了处理SAR数据（.raw）所需的相关元数据信息。

零级帧数据文件（.raw）包含了处理模拟SAR信号后收集的二进制SAR信号数据(主要数据)。

卷描述符是C包含了关于处理的简短总结。

处理信息文件提供了关于辅助信息。

元数据文件是一个ASCII文件，包含了使用ASF软件工具所需的大部分元数据。

空卷目录文件是CEOS帧分发的一部分，用于终止逻辑数据卷。

使用如下代码从给定的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');                % Speed of light

% 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).';

在matlab的工作区，可以查看rawData的格式为复数形式，数值类型为浮点型，如下图。

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')

多视处理

通过进行多视处理来减少斑点的效果，从而在图像分辨率上做出权衡。在范围或方位方向，或者在两个方向上，后续的线路被平均以获得更好的图像并减少斑点。多视可以在范围或方位或两个方向上进行。

辅助函数multilookProcessing在范围和方位方向上进行平均，分别查看4和20。

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')

总结

这个例子展示了如何提取系统参数如脉冲重复频率、波长、采样率、脉冲宽度、范围门延迟、传感器速度和原始数据。然后使用范围迁移图像形成算法聚焦原始数据。最后，应用多视处理来去除乘性噪声。

辅助函数

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')