Method schwabisch

This method is described in [15]. It is a fast method that yields the radar coded heights. It is based on the idea to first compute the reference fase at a number of heights and then to compare the actual phase from the interferogram with these values to determine the height.

A problem is that the interferogram does not contain the reference phase anymore, so that has to be added to the estimated phi.

It uses a number of steps which are described below.

1. Compute reference phase at a NL locations (line,pixel) at NH (=3) heights (0, 2000 and 4000m).
   for h=0, 2000, 4000:
   • ellips.a=wgs84.a+h, ellips.b=wgs84.b+h
   • compute position of master satellite, corresponding point P on ellips and position of slave satellite, see annex D.
   • ${B_{\parallel}}= r_1 - r_2$
   • $\phi_R = -{4\pi\over\lambda} {B_{\parallel}}$
   • store the values and locations and heights.
   Note that the reference phase defined like this typically is something like 5000 (rad) even for h=0. Therefor the reference phase for h=0 is set to 0, (because in the unwrapped interferogram the reference phase is removed, if the phase of the unwrapped interferogram is 0, then this should yield a height of 0) and the reference phase for height h is set to refphaseh - refphase0. (This makes the computations as done later a little stupid, to estimate coefficients which are by definition equal to 0.) (An other possibility is not to do this here, but later when the functions are evaluated to add the reference phase to each pixel. I have tested this and the results are identical (+- .4m))
2. Compute for each location a polynomial (1d degree 1dD (= NH-1)) to describe the height as a function of reference phase at these points.
   For each location it holds ($\phi_i$ for height i):
   

   $$h = 0 = \alpha_0 + \alpha_1 \phi_0 + \alpha_2 \phi_0^2$$ (24)

   $$h = 2000 = \alpha_0 + \alpha_1 \phi_1 + \alpha_2 \phi_1^2$$ (25)

   $$h = 4000 = \alpha_0 + \alpha_1 \phi_2 + \alpha_2 \phi_2^2$$ (26)

   
   So it is easy to solve (exact) for $\alpha_i$ per location.
3. Compute (1dD+1=NH) polynomials to describe the coefficients of the previous step as a function of location. (now, for a random location, the coefficients of height as a function of the reference phase can be computed.) For example for $\alpha_0$ computed at NL locations (l,p) a 2d polynomial can be used:
   

   $$\alpha_{0_{lp}} = \beta_{00} + \beta_{10} l + \beta_{01}p + \beta{20} l^2 + ...$$

   $$= \sum_{i=0}^{d}\sum_{j=0}^{i} \beta_{i-j,j} l^{i-j}p^{j}$$ (27)

   
   A linear system can be easily set up and solved (least squares) by cholesky factorization. A rescaling needs to be applied to avoid instability. The system can be solved simultaneously for all alphas, because the normal matrix (and factorization) remains the same, but somehow our cholesky routine introduced an error (which is probably caused by using fortran in c) so we just solve 3 seperate times with the same factored normalmatrix.
4. Evaluate for all points at (line,pixel) with ok unwrapped phase the 2d polynomial to obtain the coefficients of the height(phi) function. Then evaluate the 1d polynomial to obtain the height.
   For all (l,p) with ok unwrapped phase:
   • Compute the alphas.
     For betas appropriate to alpha0:

     $$\alpha_0 = \sum_{i=0}^{d}\sum_{j=0}^{i} \beta_{i-j,j} l^{i-j}p^{j}$$ (28)

     
     

   • For betas appropriate to alpha1:

     $$\alpha_1 = \sum_{i=0}^{d}\sum_{j=0}^{i} \beta_{i-j,j} l^{i-j}p^{j}$$ (29)

     
     

   • repeat computing alphas until you have them all (1dD+1).
   • Compute the height.

     $$h = \sum_{i=0}^{1dD}\alpha_i\phi^i$$ (30)