Transcript D4P2a - SEOM

5th Advanced Training Course in Land
Remote Sensing
Synthetic Aperture Radar
Polarimetric Tomography
Practical session
Stefano Tebaldini
Dipartimento di Elettronica, Informazione e Bioingegneria – Politecnico di Milano
TomoSAR_Main.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%
% DEMONSTRATIVE TOMOGRAPHIC SAR PROCESSING FOR FOREST ANALYSIS
% AUTHOR: STEFANO TEBALDINI, POLITECNICO DI MILANO
% EMAIL: [email protected]
% TEL: +390223993614
%
% THE FOLLOWING SCRIPT AND ALL RELATED SCRIPTS/FUNCTIONS AND DATA ARE INTENDED AS
% MATERIAL FOR PRACTICAL SESSION D4P2a OF THE ESA'S 5TH ADVANCED TRAINING COURSE IN
% LAND REMOTE SENSING ESA, TO BE HELD IN VALENCIA, SPAIN, ON SEPTEMBER 08-12 2014.
%
% THIS SOFTWARE WAS DEVELOPED AND TESTED USING MATLAB R2011b
%
% ALL RELATED SCRIPTS AND FUNCTIONS WILL BE COMMENTED IN DETAILS DURING THE PRACTICAL
% SESSION
%
% SAR DATA USED IN THIS SCRIPT ARE PART OF THE SAR DATA-SET ACQUIRED BY DLR
% IN 2008 IN THE FRAME OF THE ESA CAMPAIGN BIOSAR 2008
% DATA FOCUSING, COREGISTRATION, PHASE FLATTENING, AND GENERATION OF KZ
% MAPS WERE CARRIED OUT BY DLR.
% DATA PHASE CALIBRATION WAS CARRIED OUT BY THE AUTHOR%
%
% TERRAIN ELEVATION AND FOREST HEIGHT DATA USED IN THIS SCRIPT ARE EXTRACTED FROM
% THE LIDAR DATA-SET ACQUIRED BY THE SWEDISH DEFENCE RESEARCH AGENCY (FOI)
% AND HILDUR AND SVEN WINQUIST'S FOUNDATION IN THE FRAME OF THE ESA
% CAMPAIGN BIOSAR 2008
% PROCESSING OF LIDAR DATA AND PROJECTION ONTO SAR GEOMETRY WAS CARRIED OUT
% BY THE AUTHOR
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5th Advanced Training Course in Land Remote Sensing
Data
Campaign
BioSAR 2008 - ESA
System
E-SAR - DLR
Site
Krycklan river catchment,
Northern Sweden
Scene
Boreal forest
Pine, Spruce, Birch, Mixed stand
Topography
Hilly
Tomographic
Tracks
6 + 6 – Fully Polarimetric
(South-West and North-East)
Carrier
Frequency
P-Band and L-Band
Slant range
resolution
1.5 m
Azimuth
resolution
1.6 m
Vertical resolution
(P-Band)
20 m (near range) to >80 m (far
range)
Vertical resolution
(L-Band)
6 m (near range) to 25 m (far
range)
5th Advanced Training Course in Land Remote Sensing
Forward model
Resolution is determined by pulse bandwidth along the slant range direction, and
by the lengths of the synthetic apertures in the azimuth and cross range directions
The SAR resolution cell is split into multiple layers, according to baseline aperture
Track N
c
r 
2B
Track n
v 
r
x 
2 Av
r
2 Ax
SAR Resolution Cell
Track 1
Tomographic Res. Cell
height
height
ground range
Δr
5th Advanced Training Course in Land Remote Sensing
For most systems:
Δv >> Δr, Δx
Δv
azimuth
B: pulse bandwidth
Av: baseline aperture
Ax: azimuth aperture
λ: carrier wavelength
ground range
z  v  sin  
Vertical wavenumber
Each focused SLC SAR image is obtained as the Fourier Transform of the scene
complex reflectivity along the cross-range coordinate
4


yn r, x    sr, x, v exp  j
bn v dv
r 

Change of variable from cross range to height
yn(r,x) : SLC pixel in the n-th image
s(r,x,v): average complex reflectivity of the scene
within the SAR 2D resolution cell at (r,x)
bn : normal baseline for the n-th image
λ : carrier wavelength
yn r , x    sr , x, z  exp jk z n   z dz
kz is usually referred to as vertical wavenumber
or phase to height conversion factor
k z n  
4 bn
r sin 
5th Advanced Training Course in Land Remote Sensing
height, z
z  v  sin 
Reference height
yn r , x    sr , x, z  exp jk z n   z   dz Note: z is always intended as height with respect to
a Digital Elevation Model (DEM)
 “zero” of the resulting Tomographic profiles
Track N
Track n
Track 1
height
θ
DEM
error
azimuth
ground range
Red = Reference terrain elevation
Orange = True terrain elevation
5th Advanced Training Course in Land Remote Sensing
DEM subtraction
The dependence on height is limited to the
phase terms kzz
yn r , x    sr , x, z  exp jk z n   z   dz
 Passing from one reference DEM to another  phase steering from Z1 to Z2
yn r , x; Z 2   yn r , x; Z1  exp  jk z n  Z1  Z 2 
Track N
Track n
Track 1
Tomographic profile from
Tomographic profile from
yn(r,x;Z1)
yn(r,x;Z2)
height
θ
azimuth
DEM
error
ground range
5th Advanced Training Course in Land Remote Sensing
Red = Z1 = Reference terrain elevation
Orange = Z2 = True terrain elevation
Phase calibration
Phase jitters in different passes result in signal defocusing
• Spaceborne: tropospheric and ionospheric phase screens
• Airborne: uncompensated platform motions on the order of a fraction of a
wavelength
Slant range
θ
Replicas
Height of
ambiguity
height
height
resolution cell
Replicas
Defocusing
azimuth
ground range
5th Advanced Training Course in Land Remote Sensing
azimuth
ground range
Phase calibration
Current navigational systems employed by airborne SARs do not provide, in general, subwavelength accuracy concerning the location of one flight line with respect to another
60
50
40
30
20
10
0
-10
TomoSAR at Remningstorp, Sweden – from BioSAR 2007
Original data
Phase calibrated data
height [m]
height [m]
Need for a data-driven Phase Calibration procedure
200
600
1000
1400
slant range [m]
1800
2200
60
50
40
30
20
10
0
-10
200
600
1000
1400
slant range [m]
1800
2200
TomoSAR at Kangerlussuaq, Greenland – from IceSAR 2012
200
100
0
-100
-200
Phase calibrated data
Height [m]
Height [m]
Original data
-4000 -3000 -2000 -1000
0
1000 2000 3000 4000 5000
Azimuth [m]
5th Advanced Training Course in Land Remote Sensing
200
100
0
-100
-200
-4000 -3000 -2000 -1000
0
1000 2000 3000 4000 5000
Azimuth [m]
Phase calibration
Current navigational systems employed by airborne SARs do not provide, in general, subwavelength accuracy concerning the location of one flight line with respect to another
60
50
40
30
20
10
0
-10
TomoSAR at Remningstorp, Sweden – from BioSAR 2007
Original data
Phase calibrated data
height [m]
height [m]
Need for a data-driven Phase Calibration procedure
200
600
1000
1400
slant range [m]
1800
2200
60
50
40
30
20
10
0
-10
200
600
1000
1400
slant range [m]
1800
2200
This data-set: phase calibration by PoliMi
TomoSAR at Kangerlussuaq, Greenland – from IceSAR 2012
200
100
0
-100
-200
Phase calibrated data
Height [m]
Height [m]
Original data
-4000 -3000 -2000 -1000
0
1000 2000 3000 4000 5000
Azimuth [m]
5th Advanced Training Course in Land Remote Sensing
200
100
0
-100
-200
-4000 -3000 -2000 -1000
0
1000 2000 3000 4000 5000
Azimuth [m]
TomoSAR_Main.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%
% LOAD DATA
if not(exist('I'))
load('BioSAR_2_L_Band_sample_data')
Master = 1
[Nr,Na,N] = size(I{1})
N_pol = length(I)
rem_dem_flag = 1
if rem_dem_flag % remove dem phases (optional)
for pol = 1:N_pol
for n = 1:N
dem_phase = kz(:,:,n).*(DEM - DEM_avg);
I{pol}(:,:,n) = I{pol}(:,:,n).*exp(1i*dem_phase);
end
end
end
Ch = {'HH','HV','VV'}
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%
Notes:
Data can be referenced to a flat DEM (DEM_avg) or to the Lidar DEM (DEM)
5th Advanced Training Course in Land Remote Sensing
TomoSAR_Main.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%
% Let's look at the data first....
for pol = 1:N_pol
figure, imagesc(az_ax,rg_ax,sum(abs(I{pol}),3)), colorbar
title(Ch{pol})
xlabel('azimuth [m]')
ylabel('range [m]')
end
figure, imagesc(az_ax,rg_ax,DEM), colorbar
title('DEM [m]')
xlabel('azimuth [m]')
ylabel('range [m]')
figure, imagesc(az_ax,rg_ax,FOR_H,[0 35]), colorbar
title('Forest height [m]')
xlabel('azimuth [m]')
ylabel('range [m]')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%
Notes:
DEM = Lidar DEM
FOR_H= Lidar forest height
5th Advanced Training Course in Land Remote Sensing
Results
4
HH
x 10
4400
2.5
4600
2
range [m]
4800
1.5
5000
1
5200
0.5
5400
-100
0
100
200 300 400
azimuth [m]
5th Advanced Training Course in Land Remote Sensing
500
600
700
Results
HV
7000
4400
6000
4600
5000
range [m]
4800
4000
5000
3000
5200
2000
5400
1000
-100
0
100
200 300 400
azimuth [m]
5th Advanced Training Course in Land Remote Sensing
500
600
700
Results
VV
4400
14000
4600
12000
10000
range [m]
4800
8000
5000
6000
5200
4000
5400
2000
-100
0
100
200 300 400
azimuth [m]
5th Advanced Training Course in Land Remote Sensing
500
600
700
Results
Forest height [m]
35
4400
30
4600
25
range [m]
4800
20
5000
15
5200
10
5400
5
-100
0
100
200 300 400
azimuth [m]
5th Advanced Training Course in Land Remote Sensing
500
600
700
0
Results
DEM [m]
4400
290
280
4600
270
260
range [m]
4800
250
240
5000
230
5200
220
210
5400
200
190
-100
0
100
200 300 400
azimuth [m]
5th Advanced Training Course in Land Remote Sensing
500
600
700
TomoSAR_Main.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%
% COHERENCE EVALUATION
% estimation window (in meters)
Wa_m = 30
Wr_m = 30
[COV_4D,a_sub,r_sub] = Generate_covariance_matrix(I{1},az_ax,rg_ax,Wa_m,Wr_m);
figure, InSAR_view(abs(COV_4D),[0 1]), colorbar
title('InSAR coherences')
figure, InSAR_view(angle(COV_4D),[-pi pi]), colorbar
title('InSAR phases')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%
Notes:
COV_4D is a 4D data structure representing the complex coherence as a function of each
interferometric pair, i.e.: ynm(r,x)
Generate_covariance_matrix.m = function to evaluate COV_4D from SLC
images
InSAR_view = function to view COV_4D as a big 2D matrix
5th Advanced Training Course in Land Remote Sensing
Generate_covariance_matrix.m
function [Cov,x_sub,y_sub] = Generate_covariance_matrix(F,x_ax,y_ax,Wx_m,Wy_m)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[Ny,Nx,N] = size(F);
% pixel sampling
dx = x_ax(2)-x_ax(1);
dy = y_ax(2)-y_ax(1);
% filter along x
Lx = round(Wx_m/2/dx);
filter_x = hamming(2*Lx+1);
% sub-sampling along x
x_sub = Lx+1:max(round(Lx/2),1):Nx-Lx;
% filter along y
Ly = round(Wy_m/2/dy);
filter_y = hamming(2*Ly+1);
% sub-sampling along y
y_sub = Ly+1:max(round(Ly/2),1):Ny-Ly;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5th Advanced Training Course in Land Remote Sensing
Generate_covariance_matrix.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Covariance matrix evaluation
Nx_sub = length(x_sub);
Ny_sub = length(y_sub);
Cov = ones(Ny_sub,Nx_sub,N,N);
for n = 1:N
In = F(:,:,n); % n-th image
% second-order moment
Cnn = filter_and_sub_sample(In.*conj(In),filter_x,filter_y,x_sub,y_sub);
for m = n:N
Im = F(:,:,m);
Cmm = filter_and_sub_sample(Im.*conj(Im),filter_x,filter_y,x_sub,y_sub);
Cnm = filter_and_sub_sample(Im.*conj(In),filter_x,filter_y,x_sub,y_sub);
% coherence
coe = Cnm./sqrt(Cnn.*Cmm);
Cov(:,:,n,m) = coe;
Cov(:,:,m,n) = conj(coe);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Cnm = filter_and_sub_sample(Cnm,filter_x,filter_y,x_sub,y_sub)
% filter and sub-sample
t = Cnm;
t = conv2(t,filter_x(:)','same');
t = t(:,x_sub);
t = conv2(t,filter_y(:),'same');
t = t(y_sub,:);
Cnm = t;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
th
5 Advanced Training Course in Land Remote Sensing
InSAR_view.m
function InSAR_view(DX,cax)
[Nx_out,Ny_out,N,a] = size(DX);
if a == N
flag_4D = 1;
else
flag_4D = 0;
end
DDX = zeros(N*Nx_out,N*Ny_out);
for n = 1:N
ind_n = [1:Nx_out] + Nx_out*(n-1);
for m = 1:N
ind_m = [1:Ny_out] + Ny_out*(m-1);
if flag_4D
DDX(ind_n,ind_m) = DX(:,:,n,m);
else
DDX(ind_n,ind_m) = DX(:,:,m) - DX(:,:,n);
end
end
end
if exist('cax')==1
if max(abs(cax-[-pi pi]))==0
disp('phase')
DDX = angle(exp(1i*DDX));
end
imagesc(DDX,cax)
else
imagesc(DDX)
end
axis off
5th Advanced Training Course in Land Remote Sensing
Results – InSAR coeherences – DEM subtracted
InSAR coherences
InSAR phases
1
3
0.9
0.8
2
0.7
1
0.6
0.5
0
0.4
-1
0.3
0.2
-2
0.1
0
5th Advanced Training Course in Land Remote Sensing
-3
Results – InSAR coeherences – DEM not subtracted
InSAR coherences
InSAR phases
1
3
0.9
0.8
2
0.7
1
0.6
0.5
0
0.4
-1
0.3
0.2
-2
0.1
0
Notes:
Noticeable topographic phases
Lower coherence magnitudes
5th Advanced Training Course in Land Remote Sensing
-3
TomoSAR_Main.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%
% TOMOGRAPHIC PROCESSING (3D focusing)
% vertical axis (in meters)
if rem_dem_flag % height w.r.t. DEM
dz = 0.5;
z_ax = [-20:dz:40];
else % % height w.r.t. average DEM
dz = 1;
z_ax = [-150:dz:150];
end
Nz = length(z_ax);
% half the number of azimuth looks to be processed
Lx = 10
% azimuth position to be processed (meters)
az_profile_m = 590;
az_profile_m = 678
az_profile_m = -92
% Focus in SAR geometry
TomoSAR_focusing
if rem_dem_flag == 0
% the following routines have been written assuming DEM phases are removed
return
end
% Geocode to ground geometry and compare to Lidar forest height
Geocode_TomoSAR
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%
5th Advanced Training Course in Land Remote Sensing
25
TomoSAR_focusing.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% pixel index
[t,a0] = min(abs(az_ax-az_profile_m));
az_ind = a0 + [-Lx:Lx];
% Focusing
for pol = 1:N_pol
Tomo_3D{pol} = zeros(Nz,Nr,length(az_ind));
for z = 1:Nz
t = I{pol}(:,az_ind,:).*exp(1i*kz(:,az_ind,:).*z_ax(z));
Tomo_3D{pol}(z,:,:) = mean(t,3);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Notes:
Just a discrete Fourier Transform
yn r , x    sr , x, z  exp jk z n   z   dz
5th Advanced Training Course in Land Remote Sensing
sˆr , x, z    yn r , x  exp jk z n   z 
n
Results – Tomographic profiles
Lidar DEM not subtracted (rem_dem_flag=0)
LIDAR DEM subtracted (rem_dem_flag=1)
 Reference height = DEM_avg = 200 m
 Reference height = LIDAR DEM
HH
HH
40
height [m]
height [m]
40
20
0
-20
20
0
-20
4400
4600
4800
5000
range [m]
5200
5400
4400
4600
4800
HV
5000
range [m]
5200
5400
5200
5400
5200
5400
HV
40
height [m]
height [m]
40
20
0
-20
20
0
-20
4400
4600
4800
5000
range [m]
5200
5400
4400
4600
4800
VV
5000
range [m]
VV
40
height [m]
height [m]
40
20
0
20
0
-20
-20
4400
4600
4800
5000
range [m]
5200
5400
5th Advanced Training Course in Land Remote Sensing
4400
4600
4800
5000
range [m]
Geocoding
Tomographic profiles have been generated in the coordinate system (r,z):
o
r = (Zero-Doppler) distance from the Master track
o
z = height w.r.t. the reference DEM
 A point at coordinates (Y,Z) in the ground range plane is found at
YMaster  Y 2  Z Master  Z 2
z  Z  z ref r 
r
Master Track
z
θ
Height, Z
r
Red = zref(r) = Reference
azimuth
ground range, Y
5th Advanced Training Course in Land Remote Sensing
terrain elevation
Geocode_TomoSAR.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% pixel index
[t,a0] = min(abs(az_ax-az_profile_m));
% Master position
Sy = interp1(S{Master}.x,S{Master}.y,az_profile_m);
Sz = interp1(S{Master}.x,S{Master}.z,az_profile_m);
% Terrain elevation
dem = DEM(:,a0)';
% Forest height
for_h = FOR_H(:,a0)';
% ground range as a function of slant range
y_of_r = sqrt(rg_ax.^2 - (Sz-dem).^2) + Sy;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% absolute ground range axis
dy = 1;
y_ax_abs = [min(y_of_r)-5:dy:max(y_of_r)+5];
% absolute height axis
z_ax_abs = [min(dem)-10:dz:max(dem)+30];
% ground range as a function of slant range
y_of_r = sqrt(rg_ax.^2 - (Sz-dem).^2) + Sy;
% resample lidar dem and lidar forest height from range to ground range
dem_gr = interp1(y_of_r,dem,y_ax_abs,'linear',nan);
for_h_gr = interp1(y_of_r,for_h,y_ax_abs,'linear',nan);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5th Advanced Training Course in Land Remote Sensing
Geocode_TomoSAR.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ground range, height coordinates
[Za,Ya] = ndgrid(z_ax_abs,y_ax_abs);
% slant range
R = sqrt( (Sy-Ya).^2 + (Sz-Za).^2 );
% reference dem
Z_ref = interp1(rg_ax,dem,R,'linear','extrap');
% height w.r.t. reference dem
Z = Za - Z_ref;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Geocode tomograms
for pol = 1:3
tomo_sar = Tomo_filt{pol};
tomo_sar = tomo_sar./max(tomo_sar(:));
% Geocoded tomogram
tomo_geo = interp2(rg_ax,z_ax,tomo_sar,R,Z);
% Geocoded tomogram - height w.r.t. Lidar
tomo_geo(isnan(tomo_geo)) = 0;
for y = 1:length(y_ax_abs)
tomo_geo_rel(:,y) = interp1(z_ax_abs,tomo_geo(:,y),z_ax + dem_gr(y));
end
%%%%%%%%%%%%%%%%% Draw pictures here%%%%%%%%%%%%%%%%%%%%%%5
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5th Advanced Training Course in Land Remote Sensing
Geocoding - results
VV SAR Geometry
40
Slant range – height w.r.t. reference
height [m]
SAR geometry
20
0
DEM
-20
4400
Note: Lidar forest height not matched
4600
4800
5000
range [m]
5200
5400
VV Ground Geometry
Ground geometry
Ground range – height
height [m]
280
260
240
220
200
180
2200
2400
2600
2200
2400
2600
2800
3000
3200
ground range [m]
VV Ground Geometry relative to DEM
3400
3600
3800
3400
3600
3800
40
DEM
Ground range – height w.r.t. reference
DEM
height [m]
Ground geometry w.r.t. reference
20
0
-20
Note: Lidar forest height well matched
5th Advanced Training Course in Land Remote Sensing
2800
3000
3200
ground range [m]
Red = Lidar terrain
Black = Lidar forest height
TomoSAR as a Spectral Estimation Problem
Performances are often limited by baseline sparseness and aperture
 SAR Tomography is commonly rephrased as a Spectral Estimation problem,
based on the analysis of the data covariance matrix among different tracks
General Procedure
Form the MB data
Evaluate the sample
Evaluate the cross-range distribution
vector [Nx1]
covariance matrix [NxN]
of the backscattered power through
by local multi-looking
some Spectral Estimator
 y1 r, x  
 y r, x 

y MB   2
  
 y r, x 
 N

ˆ  y MB y HMB
R
sˆr, x, v   Sˆ v 
2
L
Remark: it is customary to normalize R such that entries on the main diagonal are unitary
 R is the matrix of the interferometric coherences for
all baselines
5th Advanced Training Course in Land Remote Sensing
R
nm

E  yn ym* 
  
2
E yn E ym
2
  nm
TomoSAR as a Spectral Estimation Problem
• Beamforming:
inverse Fourier Transform; coarse spatial resolution; radiometrically consistent
ˆ az 
Sˆ z   aH z R
az   exp jkz 0z  exp jkz 1z   exp jkz N 1z 
T
Note: can also be obtained by averaging Tomograms
yn r , x    sr , x, z  exp jk z n   z   dz
Sˆ  z  
sˆr , x, z 

sˆr , x, z 
2
L
2
dz
L
5th Advanced Training Course in Land Remote Sensing
sˆr , x, z    yn r , x  exp jk z n   z 
n
TomoSAR as a Spectral Estimation Problem
• Capon Spectral Estimator:
spatial resolution is greatly enhanced, at the expense of radiometric accuracy;
Sˆ z  
1
ˆ 1az 
a H z R
Warning: many looks required to stabilize matrix inversion
Rationale: find a filter f = f(z) such that:
o S(z) = fHRf is minimized  optimal noise rejection
o fHa = 1 preservation of the signal component at height z
Result:
ˆ 1a z 
R
f H
ˆ 1a z 
a  z R
ˆ f
Sˆ  z   f H R
5th Advanced Training Course in Land Remote Sensing
1
ˆ 1a z 
a H  z R
TomoSAR_Main.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%
% TOMOGRAPHIC PROCESSING (cast as a spectral estimation problem)
% estimation window (in meters)
Wa_m = 10
Wr_m = 3
% coherence evaluation
for pol = 1:N_pol
[COV_4D_pol{pol},a_sub,r_sub] =
Generate_covariance_matrix(I{pol},az_ax,rg_ax,Wa_m,Wr_m);
end
for pol = 1:N_pol
figure, InSAR_view(abs(COV_4D_pol{pol}),[0 1]), colorbar
title(['Multi-baseline and multi-polarization InSAR coherences - ' Ch{pol}])
figure, InSAR_view(angle(COV_4D_pol{pol}),[-pi pi]), colorbar
title(['Multi-baseline and multi-polarization InSAR phases - ' Ch{pol}])
end
% Spectral estimation
TomoSAR_Spectral_Estimation
% Geocode to ground geometry and compare to Lidar forest height
Tomo_filt = Sp_est_i;
Geocode_TomoSAR
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%
5th Advanced Training Course in Land Remote Sensing
TomoSAR_Spectral_Estimation.m
% pixel index
az_ax_sub = az_ax(a_sub);
[t,a0] = min(abs(az_ax_sub-az_profile_m));
%
kz_sub = kz(r_sub,a_sub,:);
clear Sp_est
flag_Capon = 1
for r = 1:length(r_sub)
local_kz = squeeze(kz_sub(r,a0,:));
A = exp(-1i*z_ax(:)*local_kz');
%
for pol = 1:N_pol
cov = squeeze(COV_4D_pol{pol}(r,a0,:,:));
if flag_Capon % Capon
cov_i = inv(cov + 1e-2*eye(N));
Si = real(diag(A*cov_i*A'));
Sp_est{pol}(:,r) = 1./Si;
else % Fourier
M = A*cov*A';
Sp_est{pol}(:,r) = real(diag(A*cov*A'));
end
end
end
for pol = 1:N_pol
Sp_est_i{pol} = interp1(rg_ax(r_sub),Sp_est{pol}',rg_ax)';
end
%%%%%%%%%% draw pictures here %%%%%%%%%
5th Advanced Training Course in Land Remote Sensing
TomoSAR as a Spectral Estimation Problem
Beamforming:
Capon spectrum:
HH
40
20
20
height [m]
height [m]
HH
40
0
-20
4400
4600
4800
5000
range [m]
5200
0
-20
5400
4400
4600
4800
40
20
20
0
-20
4400
4600
4800
5000
range [m]
5200
-20
5400
4400
4600
4800
20
20
0
4600
4800
5000
range [m]
5000
range [m]
5200
5400
5200
5400
VV
40
height [m]
height [m]
VV
4400
5400
0
40
-20
5200
HV
40
height [m]
height [m]
HV
5000
range [m]
5200
5400
5th Advanced Training Course in Land Remote Sensing
0
-20
4400
4600
4800
5000
range [m]
TomoSAR as a Spectral Estimation Problem – Ground geometry
Beamforming:
Capon spectrum:
VV SAR Geometry
40
20
20
height [m]
height [m]
VV SAR Geometry
40
0
-20
0
-20
4400
4600
4800
5000
range [m]
5200
5400
4400
4600
4800
280
260
260
height [m]
height [m]
280
240
220
240
220
180
2200
2400
2600
2800 3000
3200
ground range [m]
3400
3600
3800
2200
2400
40
20
20
height [m]
40
0
2400
2600
2800 3000
3200
ground range [m]
3400
2800 3000
3200
ground range [m]
3400
3600
3800
3600
3800
0
-20
2200
2600
VV Ground Geometry relative to DEM
VV Ground Geometry relative to DEM
height [m]
5400
200
200
-20
5200
VV Ground Geometry
VV Ground Geometry
180
5000
range [m]
3600
3800
5th Advanced Training Course in Land Remote Sensing
2200
2400
2600
2800 3000
3200
ground range [m]
3400
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%% THANK YOU %%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FEEL FREE TO CONTACT ME AT:
% EMAIL: [email protected]
% TEL: +390223993614
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5th Advanced Training Course in Land Remote Sensing