Report

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 sr, 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 sr , 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 sr , 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 sr , 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 sr , 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 ˆ az Sˆ z aH z R az exp jkz 0z exp jkz 1z exp jkz N 1z T Note: can also be obtained by averaging Tomograms yn r , x sr , 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 ˆ 1az 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