干涉合成孔径雷达（InSAR）工具箱.rar

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