siminterf.m

荷兰Delft大学开发的insar（干涉合成孔径雷达）图像处理部分源代码

function [interf, DEM, noise, coherence, watermask] = siminterf(Bperp,fracdim,water,maxheight,numlines,numpixels,doplot,donoise,dorefpha)
% SIMINTERF -- SIMulate phase of INTERFerogram (not amplitude).
%
%   Function to simulate an (unwrapped) interferogram by 'radarcoding'
%   a terrain model (DEM). 'Radarcoding' means converting the terrain
%   coordinates to the radar system [azimuth lines, range pixels]. 
%   The terrain model either is a geometric body (i), a fractal DEM (ii),
%   or an input matrix (iii). The range pixel spacing of the DEM is 80 meters.
%
%   The algorithm used radarcodes the DEM per line by computing the slant
%   range r to master and slave satellite for each point of the DEM, yielding
%   a function phase(r). This function is interpolated to a regular grid
%   (radar range pixel coordinate system). Zero Doppler data is assumed,
%   thus the azimuth coordinates are equal in both systems. The phase of 
%   the 'flat Earth' is computed using a far field approximation, and
%   subtracted by default. To compute r, the position of master and slave
%   satellite are fixed by using a certain Bperp and baseline orientation.
%   The range, incidence angle to the first terrain point, and the wavelength
%   are also fixed. The orbits are assumed not to converge (easy to simulate.)
% 
%   Input parameters:
%     INT = SIMINTERF by itself simulates a DEM of a cone and the
%     corresponding interferogram. It is verbose and makes plots.
%
%     SIMINTERF(Bperp) uses the specified perpendicular baseline. A larger
%     baseline corresponds with a smaller height ambiguity (more fringes).
%     A crude approximation of the height ambiguity is ha = 10000/Bperp.
%
%     SIMINTERF(Bperp,D) where D is the fractal dimension of the simulated
%     DEM. Smaller D implies smoother surface. Earth topography can be simulated
%     by a fractal dimension of approximately 2.3.
%       D == DEM     Use this input matrix as DEM; 
%                     Input arguments: water, height, lines, pixels
%                     are not assigned if this option is used.
%             0      Simulate a cone;
%            -1      Simulate a ramp;
%            -2      Simulate pyramid;
%             other  Use this fractal dimension to simulate a DEM.
%
%     SIMINTERF(Bperp,D,WATER) uses the specified percentage to 
%     create water areas of height 0 (approximately). The phase of 
%     these areas is uniform noise, the coherence is gaussian below 0.2.
%
%     SIMINTERF(Bperp,D,water,HEIGHT) simulates a DEM with the
%     specified maximum HEIGHT. (The minimum height always equals 0.)
%
%     SIMINTERF(Bperp,D,water,height,LINES) creates an interferogram
%     with the specified number of azimuth lines.
%
%     SIMINTERF(Bperp,D,water,height,lines,PIXS) creates an
%     interferogram with the specified number of range pixels.
%
%     SIMINTERF(Bperp,D,water,height,lines,pixs,DOPLOT) where
%     DOPLOT=0 prevents the generation of plots. Figures 1:6 are used.
%
%     SIMINTERF(Bperp,D,water,height,lines,pixs,doplot,DONOISE)
%     where DONOISE = 0 disables noise generation.
%
%     SIMINTERF(Bperp,D,W,H,lines,pixs,doplot,donoise,DOREFPHA)
%     where DOREFPHA is 0 if reference phase does not have to be
%     subtracted.
%
%   Defaults:
%     Bperp    = 150   (meters)
%     D        = 0     (i.e. a cone)
%     WATER    = 20    (percentage)
%     HEIGHT   = 700   (meters)
%     LINES    = 512   (number of azimuth lines)
%     PIXELS   = 512   (number of range bins in range)
%     DOPLOT   = 1     (do plot results)
%     DONOISE  = 1     (do add noise based on terrain slopes)
%     DOREFPHA = 1     (do subtract reference phase of 'flat Earth')
%
%   Additional output:
%     [INT,DEM] = SIMINTERF(...) optionally outputs the simulated DEM (in 
%     terrain coordinates).
%     [INT,DEM,NOISE] = SIMINTERF(...) optionally outputs the noise matrix
%     that is added to the INTerferogram based on the geometrical decorellation.
%     [INT,DEM,NOISE,COHERENCE] = SIMINTERF(...) optionally outputs the
%     coherence map computed from the terrain slopes (geometric decorrelation).
%     [INT,DEM,NOISE,COHERENCE,WATERMASK] = SIMINTERF(...) optionally 
%     outputs the waterarea index vector as returned by FIND.
%
%   EXAMPLES:
%   To simulate an interferogram with a baseline of 100 meter,
%   rough terrain, and 30% water area:
%     Bperp = 100; D=2.4; water=30;
%     siminterf(Bperp,D,water);
%
%   To fastly simulate some interferograms to test this script:
%     Bperp=100;  D=2.1;    water=0; H=500;
%     nlines=100; npix=200; doplot=1; donoise=1; doflat=1;
%     siminterf(Bperp,D,water,H,nlines,npix,doplot,donoise,doflat);
%
%   To radarcode an input DEM:
%     Bperp=50; X=256.*peaks(256);
%     siminterf(Bperp,X);
%
%   To view the unwrapped interferogram, and obtain the coherence map:
%     Bperp=200;  D=2.3;    water=30; H=100;
%     nlines=256; npix=256; doplot=0; donoise=1; doflat=1;
%     [INT,DEM,NOISE,COH] = ...
%       siminterf(Bperp,D,water,H,nlines,npix,doplot,donoise,doflat);
%     figure(1);
%       subplot(221), imagesc(DEM); axis 'image'; title 'DEM'; colorbar;
%       subplot(222), imagesc(INT); axis 'image'; title 'INT'; colorbar;
%       subplot(223), imagesc(wrap(INT)); axis 'image'; title 'W(INT)'; colorbar;
%       subplot(224), imagesc(COH); axis 'image'; title 'COH'; colorbar;
%
%   To make the interferogram complex, use e.g.,  
%     cint   = complex(cos(interf),sin(interf));
%   or:
%     ampl   = ones(size(interferogram));
%     cint   = ampl.*complex(cos(interf),sin(interf));
%
%   See also FRACTOPO, ANAFRACDEMO, ANGLE, COMPLEX, WRAP, CONE, PYRAMID, HEIGHTAMB
%
%   The DEOS FRACTAL and INSAR toolboxes are used
%   (www.geo.tudelft.nl/doris/)
%

% Created by Bert Kampes 05-Oct-2000
% Tested by Erik Steenbergen
%
%   TODO:
%     -more geometrical figure if fracdim=-1,-2, etc., 
%     -demo for people new to insar
%


%%% better switch more, but how, get(0)?
more off

%%% Set defaults for general parameters.
%%% Handle input; could be done smarter...
% (Bperp,dem,water,maxheight,numlines,numpixels,doplot,donoise,dorefpha)
%   1     2   3      4         5         6         7     8       9
if (nargin<9) dorefpha  = 1;    end;
if (nargin<8) donoise   = 1;    end;
if (nargin<7) doplot    = 1;    end;
if (nargin<6) numpixels = 512;  end;
if (nargin<5) numlines  = 512;  end;
if (nargin<4) maxheight = 700; end;
if (nargin<3) water     = 20;   end;
if (nargin<2) fracdim   = 0;    end;%  cone
if (nargin<1) Bperp     = 150;  end;

%%% Set output parameters.
% [interferogram, DEM, noise, coherence, watermask]
%     1            2    3       4          5
if (nargout>=3) donoise = 1;  end;
%if (nargout<1) siminterf(); end;



%%% Check if input fracdim (D) was a matrix (use it as DEM).
if (prod(size(fracdim))>1) 
  DEM                  = 999;			% 999 identifier use input matrix
  [DEM,fracdim]        = deal(fracdim,DEM);	% swap
  [numlines,numpixels] = size(DEM);
  maxheight            = max(DEM(:));
  if (nargin<3)  water = 0;   end;
end;


%%% Fixed variables (do not change w/o reason).
alpha           = deg2rad(10.);%	[rad] baseline orientation
lambda          = 0.05666;%		[m]   wavelength
theta           = deg2rad(19.);%	[rad] looking angle to first pixel
R1              = 830000.;%		[m]   range to first point
pi4divlam       = (-4.*pi)./lambda;


%%% Provide some info.
disp(['Lines (azimuth):        ', num2str(numlines)]);
disp(['Pixels (range bins):    ', num2str(numpixels)]);
disp(['Perpendicular baseline: ', num2str(Bperp), ' m']);
disp(['Height ambiguity:       ', num2str(round(heightamb(Bperp))), ' m']);
disp(['Fractal dimension DEM:  ', num2str(fracdim)]);	
disp(['Water area:             ', num2str(water), '%']);	
disp(['Minimum height DEM:     ', num2str(0), ' m']);
disp(['Maximum height DEM:     ', num2str(maxheight), ' m']);
disp(['Make plots:             ', num2str(doplot)]);	
disp(['Compute noise:          ', num2str(donoise)]);	
disp(' ');


%%% Simulation fractal DEM or cone.
switch fracdim
  case  999
    disp('Method DEM:             Using input matrix.');
  case  0
    disp('Method DEM:             Creating cone');
    DEM   = cone(numlines,numpixels);
  case -1
    disp('Method DEM:             Creating ramp');
    DEM   = ones(numlines,1) * lying(linspace(0,maxheight,numpixels));
  case -2
    disp('Method DEM:             Creating pyramid');
    DEM   = pyramid(numlines,numpixels);
  otherwise
    disp('Method DEM:             Creating fractal');
    %%% Create DEM of dimensions (numlines)x(numpixels)
    %%% See script fractopo for improvements & water areas.
    %%% and 3D visualization.
    beta = 7 - 2*fracdim;	% valid for 2D area [pr.corr.RH]
    DEM = fracsurf(numlines,beta,'n',1000);
end


%%% Add water and rescale to [0:maxheight].
if ( fracdim ~= 999 )		% input matrix
  disp(['Rescaling DEM:          [0:',num2str(maxheight),']']);
  DEM = DEM(1:numlines,1:numpixels);
  minDEM         = min(DEM(:));
  maxDEM         = max(DEM(:));
  waterlevel     = minDEM + 0.01*water.*(maxDEM-minDEM);
  DEM            = (DEM-waterlevel) .* (maxheight./(maxDEM-waterlevel));
  DEM(DEM<0)     = 0.;
end


%%% The actual radarcoding by interp1.
% use nearest, certainly for DEM fractal (?)
% due to problems with layover areas with cubic method.
%// BK 21-Nov-2000
method = 'nearest';
%method = 'linear';
%method = 'cubic';
disp(['Radarcoding DEM:        ', method ,' interpolation method']);
numpixelsdem = size(DEM,2);%			in (oversampled) DEM
dx           = 80;%    				[m] DEM resolution
%dx           = scenewidth/(numpixelsdem-1);%	[m] DEM resolution
scenewidth   = dx.*(numpixelsdem-1);%		[m]   ground range
x0           = sin(theta) .* R1;%		x coord. first DEM point
sat1_x       = 0.;%				x coord. of master satellite
sat1_y       = cos(theta) .* R1 + DEM(1,1);%	y coord. (height)
Margin       = maxheight;% [m] be on save side
maxrange     = sqrt((x0+(numpixelsdem-1)*dx).^2+sat1_y.^2)-Margin;
R1extra      = R1+Margin;
totalrange   = maxrange-R1extra;
rangebinsize = totalrange/numpixels;
rangegrid    = R1extra:rangebinsize:maxrange-rangebinsize;

%%% Give some more info.
disp(['Wavelength:             ', num2str(round(lambda*1000)./10), ' cm']);
disp(['Resolution of DEM:      ', num2str(round(dx)), ' m']);
disp(['Width of scene:         ', num2str(round(scenewidth/1e3)), ' km']);
disp(['Resolution of interf:   ', num2str(round(rangebinsize)), ' m']);
disp(['Satellite height:       ', num2str(round(sat1_y./1e3)), ' km']);
disp(['Looking angle:          ', num2str(rad2deg(theta)), ' deg']);