ERS SAR 生データ抽出とイメージ形成

R2022a 以降
ライブ スクリプトを開く

この例では、欧州リモート センシング (ERS) 合成開口レーダー (SAR) システムのパラメーターと集束していない生データを抽出し、レンジ移行イメージ形成アルゴリズムを使用して、集束したイメージを生データから生成する手順を示します。

ERS データセットは、欧州宇宙機関の ERS 衛星 (ERS-1 および ERS-2) から取得した SAR データを NASA が提供したものです。Earth observation Heritage Data Program (LTDP+) の一環として、ERS ミッションは、科学者に対し、地球のダイナミクスをより深く理解するのに役立つ、歴史的に正確で容易にアクセスできる情報を提供します。ERS データセットには、集束していない生データと共に、シーンの集束イメージの生成に使用できるシステム パラメーター ファイルが含まれています。

このワークフローを説明するために、NASA のアラスカ衛星施設 (ASF) [1] によって公開されている ERS データセットを使用します。データセットはこちらからダウンロードできます。ここでの目標は、集束したイメージを生データから生成するモデルを開発することです。

データセットのダウンロード

この例では、地球観測衛星委員会 (CEOS) 標準ファイル、すなわちリーダー ファイル (.ldr)、レベル ゼロ フレーム データ ファイル (.raw)、ボリューム記述子ファイル (.vol)、処理情報ファイル (.pi)、メタデータ ファイル (.meta)、およびヌル ファイル (.nul) を含む ERS データセットを使用します。CEOS リーダー ファイルには、付属の SAR データ (.raw) を処理するために必要な関連メタデータ情報が含まれています。レベル ゼロ フレーム データ ファイル (.raw) には、アナログ SAR 信号を処理した後に収集されたバイナリ SAR 信号データが含まれています。ボリューム記述子は CEOS フレーム配布の一部であり、処理に関する簡単な概要が含まれています。処理情報ファイルには、データ セット自体とそのシステム上の場所を記述した CEOS 変換に関する情報が含まれています。メタデータ ファイルは、ASF ソフトウェア ツールを使用するために必要なメタデータのほとんどを含む ASCII ファイルです。ヌル ボリューム ディレクトリ ファイルは、CEOS フレーム配布の一部です。これは論理データ ボリュームを終了させるために使用されます。

helperDownloadERSData 補助関数を使用して、指定された URL からデータセットをダウンロードします。データセットのサイズは 134 MB です。ERSData.zip には、リーダー ファイル、生ファイル、ライセンス ファイルの 3 つのファイルが含まれています。

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

インターネット接続の速度によっては、ダウンロード プロセスに時間がかかることがあります。このコードは、ダウンロード プロセスが完了するまで、MATLAB® の実行を一時停止します。または、Web ブラウザーを使用してデータ セットをローカル ディスクにダウンロードし、ファイルを抽出することもできます。その場合は、コード内の変数 outputFolder を、ダウンロードしたファイルの場所に変更します。

パラメーター抽出

SAR イメージ形成には、生データ ファイルに付属するパラメーター ファイル (E2_84699_STD_L0_F299.000.ldr) からパラメーターを抽出することが含まれます。パラメーター ファイルには、イメージを形成するために必要な衛星およびシーン固有のデータが含まれています。パルス繰り返し周波数、波長、サンプル レート、チャープ信号のパルス幅、レンジ ゲート遅延、センサーの速度など、いくつかのパラメーターは、パラメーター ファイルの説明に記載されている特定のアドレスに従って、パラメーター ファイルから抽出されます。有効速度、最小距離、チャープ帯域幅、チャープ周波数などの他のパラメーターは、パラメーター ファイルから抽出されたデータを使用して推定されます。

この例では、E2_84699_STD_L0_F299.000.ldr を入力ファイルとして受け取る補助関数 ERSParameterExtractor を使用して、システム パラメーターを抽出します。

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

生データ抽出

生データは、合計 28652 行、1 行あたり 11644 バイトのデータ ファイル (E2_84699_STD_L0_F299.000.raw) から抽出されます。各行において、ヘッダーは 412 バイトで構成され、残りの 11232 バイト (5616 個の複素数) が各エコーのデータとなります。ERS レーダーの場合、合成開口は 1296 点の長さです。したがって、補助関数 ERSDataExtractor は 2 のべき乗である 2048 行を読み取ります。補助関数は、ビーム幅を 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 レーダーの場合、レーダーから線形周波数変調されたチャープ信号が発信されます。phased.LinearFMWaveform System object は、抽出されたシステム パラメーターを使用して線形 FM パルス波形を作成します。スクイント角は、ドップラー周波数を用いて計算され、ERS システムではほぼゼロとなります。レンジ移行イメージ形成アルゴリズムは、周波数領域アルゴリズムであり、波数領域処理法としても知られています。rangeMigrationLFM 関数を使用して、生データを SLC イメージに集束させます。

SAR SLC イメージは、乗法性ノイズとしてモデル化できるスペックルによって特徴付けられ、結果として SAR イメージにごま塩のような外観が生じます。

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

Figure contains an axes object. The axes object with title SLC Image, xlabel Azimuth bin, ylabel Range bin contains an object of type image.

マルチルック処理

スペックルの影響は、イメージ解像度とのトレードオフを行うマルチルック処理を実行することによって軽減されます。レンジ方向、方位角方向、またはその両方において、後続のラインを平均化することで、スペックルを低減したより鮮明なイメージが得られます。マルチルックは、レンジ方向、方位角方向、またはその両方で実行できます。

補助関数 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')

Figure contains an axes object. The axes object with title Multi-look Image, xlabel Azimuth bin, ylabel Range bin contains an object of type image.

まとめ

この例では、パルス繰り返し周波数、波長、サンプル レート、チャープ信号のパルス幅、レンジ ゲート遅延、センサーの速度、生データなどのシステム パラメーターを抽出する方法を示します。次に、レンジ移行イメージ形成アルゴリズムを使用して、生データを集束させます。最後に、乗法性ノイズを除去するためにマルチルック処理を適用します。

補助関数

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

参考文献

[1] Dataset: ERS-2, ESA 2011. Retrieved from ASF DAAC 10 November 2021.