Estimating the derivative of modulo-mapped phases
The determination of the unambiguous phase from noisy observations of complex
angularly modulated signals is an unsolved problem in general, especially if
phase and amplitude are mutually uncorrelated or even independent. This is
clearly the case for a complex SAR interferogram. In terms of signal theory a
SAR interferogram can be considered as a complex, simultaneously amplitude and
phase modulated 2D signal with non Gaussian error statistics. Usually the wanted
interferometric phase is obtained by a simple tan-1 operation
delivering phase values within the principal interval (e.g. -, , depending on
the kind of inverse function). These phases contain all the information needed
for the generation of digital terrain elevation maps of observed areas but they
do not contain that information in an unambiguous way as any absolute phase
offset (an integer multiple of 2) is lost. Furthermore they are subject to phase
noise coming from the superimposed amplitude noise in real and imaginary part of
the InSAR image. The process of resolving these phase ambiguities is usually
called phase unwrapping which in terms of signal theory is simply a two
dimensional phase demodulation problem. Nearly all classical approaches to phase
unwrapping, known from optical interferometry apply a sequence of
differentiating, taking the principal value of the discrete derivative and
integrating again along specified paths. The paper will show that such a
sequence of operations always yields more or less badly biased estimates of the
wanted derivative of the unambiguous phase. This especially applies to a
combination with Linear Least Squares techniques which are commonly used to
reduce stochastic phase errors. The paper will present the analysis, the
numerical, and the analytical evaluation of the phase bias depending on the
phase slope and on the signal to noise ratio. Finally some hints to circumvent
the problem will be given.
1. Theory and Concepts
1.1 1D Stochastic Analysis
Let the observed phase obtained from a tan-1-operation form real and imaginary part of the interferogram, be related to the unambiguous phase by the following mapping:
where (k) is the true unambiguous phase at time or point k, is the phase error and the bracket indicates the operation of taking the principal value of the argument phase term in a way that:
As a result of equations 1 and 2 the observed phase always lies within the base interval (-,. This effect, arising anytime when the interferometric phase is computed by a tan-1() operation from quadrature and inphase component of an interferogram, is well known as the wrap around effect.
Nearly all known phase unwrapping techniques with the exception of /LoKr1, LoKr2/ try to unwrap the mapped phases by a sequence of differentiating, taking the principal value of the discrete derivative and integrating again.
The operation of forming the discrete derivative yields:
In the sequence of operations in equation 3 we have only used the fact that adding any integer multiple of 2 does not change the result of a modulo-2 operation. are the phase errors at points k and k+1, respectively, mapped into the base interval (-,. The stochastic properties of these errors, namely distribution density and second order moments are known. Now again using the same identities in equation 3 as before we can further write:
where the last equality has exploited that along normal integration paths the true phase variation's modulus is always smaller than . Normal refers to those integration paths which do not cross a cutline. is the true discrete phase derivative.
Equation 4 clearly expresses the error which we commit when forming the discrete derivative from modulo-2 mapped noisy data. If there were no phase error present the result would be totally correct, but since phase errors always occur in normal interferograms, we commit a systematic error when 'differentiating' modulo-2 mapped data. The phase derivative will be - depending on the noise - more or less badly biased and so will be the unwrapped phase if the biased derivative is integrated again.
The further investigation of the bias error will be organized as follows. In a first step we will investigate the stochastic features of the phase error difference in the second term of equation 4. Then we will evaluate, how a modulo operation, which occurs twice, changes the known distribution density of random variable. Let us introduce the non-mapped phase error difference variable by:
Obviously the numerical values can vary between 2, as any of the terms in the difference can vary between . Assuming that the phase errors of two subsequent phase samples are independent random variables with identical and symmetric distributions the resulting distribution density of the phase difference is the correlation product of the individual phase distribution densities. Thus we have:
Later on we will need the distribution function rather than the density. This function is simply obtained by:
Conceptually equation 7 can be solved since the individual terms are known from /Mid1, LoKr1/. The phase error distribution density, given there, is:
where is the coherence at point k. The corresponding distribution function which we also need in equation 7, obtained from /Mid1/ is given by:
As indicated earlier the phase difference given in equation 5 and showing values between 2 is now mapped into the base interval between . The functional mapping is:
The distribution density of the -mapped phase error difference can be obtained by letting:
Using the correct functional mapping in each of the intervals and exploiting the fact the conditional probability density of a functionally mapped variable only consists of a Dirac impulse we have:
Now we return to the sum in equation 4 and introduce the non-mapped discrete difference by:
Further introducing the conditional density of conditioned on the fact that the true phase derivative takes on the value we may write:
In the last equality we have used the fact that the -mapped phase error difference is independent of the true phase derivative . Substituting equation 12 into 14 we obtain:
This is the conditional distribution density of conditioned on the fact that the true phase derivative takes on the value . This variable is -mapped again to yield (Equ. 4). Now utilizing the same arguments and reasoning as before, we get the final result for the conditional density of the 'mapped' phase derivative, conditioned on the fact that the true phase derivative takes on the value :
where the short hand expression has been utilized for convenience. This expression is the 2-cutout of the sum of three shifted replicas of the distribution density (cf. Equ. 15):
Figure 1 demonstrates the generation of for an arbitrary density :
To simplify the further derivation we introduce the bias error of the 'mapped' derivative:
If we now evaluate the conditional density of this error conditioned on , we can evaluate any stochastic measure of the bias conditioned on any value of the derivative, which we seek, e.g. the conditional mean of the bias. From probability theory we know that:
Substituting the identity of equ. 19 into equation 16 we obtain the wanted conditional density:
with given in equation 17. Figure 2 demonstrates the meaning of equation 20 graphically:
With the help of figure 2 the conditional expectation of the bias error is readily calculated:
1. As our first case we will consider the interval . For this case we can subdivide the integral into the following two parts:
2. The second case is given by: . Here the following sequence of operations is valid:
Since is an even density function symmetric around zero, the corresponding distribution will show the following symmetry:
From inspecting equations 22 and 23, respectively, we conclude that the conditional mean of the bias error is an odd function with respect to the nominal value :
1.2 Preliminary Observations
Summarizing equations 22-25 we note that:
a) Maximum Phase Slope
Inspecting equations 22 and 23 we note that the bias error clearly depends on the value of the true phase slope. Knowing that:
we conclude from equations 23 and 25 that:
Thus phase slopes of = will always - even for good coherence - be estimated with the maximum possible bias error of -sign()!
b) Ideal Case: Perfect Correlation
In this case we have:
and using equations 22 and 23 we note that:
c) Worst Case: No Correlation
In the zero correlation case we get a uniform distribution density and a ramp distribution:
and from equations 22 and 23 we get the final result:
so that in this case the estimated phase slope will always be zero, which is quite remarkable.
1.3 The General Case - Numerical Results
From equations 22, 23 we conclude that knowing F0() is completely sufficient for determining the bias error. If we furthermore restrict us to the case of phase slopes between , we can utilize equation 17 and write:
where in the last equality we have only used the usual symmetry properties (cf. equ. 24). For convenience we will furtheron restrict ourselves to the case of positive slopes so that we can substitute equation 32 into equation 23 and write:
The distribution density is periodic with respect to 4. This means that we can expand it in a Fourier series:
Substituting equation 34 into 33 we readily obtain:
where the Fourier coefficients are given by:
Since we do not want to solve the integral analytically we approximately calculate dm by FFT-techniques by:
is the sampled continuous density where:
The continuous density is, as indicated by equation 6, the continuous correlation:
The discrete equivalent employing the sampled versions of the individual densities is given by:
Realizing this discrete convolution as a cyclic convolution we carry over to FFT-techniques by writing:
The result of equation 41 can be easily obtained in the frequency domain by letting:
where the Fourier transforms are calulated by::
Then the rule for approximately evaluating the bias is:
Finally we obtain the solution for negative phase slopes by (equation 25):
Equations 44 and 45 provide the final result and form the basic framework for evaluating the bias error depending on the phase slope itself as well as on the form of the densities. These densities depend on the degree of coherence or on the SNR of the interferogram (the quality of the fringes). In the following we will give some quantitative results. We will assume identical distribution densities for two successive points. Figure 3 shows the wrapped phase error density for different coherence values.
Figure 4 shows the distribution density of the unwrapped phase slope error for different degrees of coherence. Figure 5 shows the outcoming bias over the true phase slope evaluated for different degrees of coherence.
2.0 Approaches to solve or circumvent the problem
The following chapter will give a short overview about how to avoid or circumvent the problem of biased phase derivative estimates.
Any filter operation to remove the stochastic influences from the phase derivative will produce an estimate which is 'nearer' to the conditional mean, which is not identical with the true phase slope. On the other hand the simple branch/cut methods which do not apply any filtering to the phase slope estimates provide noisy but unbiased phase estimates. If any filtering is to be applied it should be applied to complex data rather than to the phases or phase slopes.
Since the bias of the slope estimation clearly depends on the slope itself, one method to keep the bias small might be to keep the phase slope small. This can be achieved by successive flattening. The biased phase estimates obtained in the first run are utilized to 'demodulate' the complex interferogram. Then the residual phase slope is estimated again, this time with a smaller bias. The unambiguous phase is generated and used again to demodulate the interferogram. The whole loop is repeated until the bias has been reduced satisfactorily. Such a procedure (as reported in /For1 / clearly yields asymptotically unbiased phase estimates.
With the results given in equations 44 and 45 it should be possible to further utilize the phase slope estimator based on finite differences of wrapped phases and to eliminate the bias by simply subtracting it. This method would be advantageously applicable to homogeneous scenes with a sufficiently smooth coherence distribution, since then the densities would not have to be calculated and Fourier transformed for any individual pixel. In this case the additional computational burden would be moderate. The bias estimation of equations 44 and 45 would be the key to maintain Linear Least Squares phase unwrapping approaches.
Clearly the best solution to a problem is an approach which prevents the problem from arising. This can be achieved by applying unbiased phase slope estimators. All these estimators share the common property that they operate on complex data rather than on the phases. Either they use real and imaginary part of the interferogram and exploit the argument of a complex correlation kernel (known from Madsen's Correlation Doppler Estimator, as proposed by /Bam1/) or they operate in the power spectral density domain (such as the Local Spectral Mode Estimator, proposed by /KrLo2, LoKr2/). Another approach would be to use a nonlinear estimator in form of an Extended Kalman Filter which does not explicitly differentiate any mapped phases (as proposed in /LoKr1, KrLo1/). Recently a combination of local slope estimation and Kalman filtering techniques has been proposed in /KrLo2, LoKr2/. This combination seems to be the most powerful approach to phase unwrapping, yielding unbiased and nearly perfectly noisefree unwrapped phases down to coherence values of 0.3 without any prefiltering!
The paper has presented the analysis, the derivation and evaluation of the estimation bias if the phase slope is determined form modulo-2 mapped phases. A finite series representation for the phase slope bias has been given, where the coefficients can be easily determined from the FFT spectrum of the distribution density of the wrapped phase error. This kind of distribution density is known for a lot of special cases even if 'multi-look' prefiltering is applied. Thus the technique is widely applicable. The numerical results calculated with the approach are in perfect agreement with the expectation. Finally some hints to circumvent or solve the problem were presented.
The work reported in the paper is an integral part of our activities on phase
unwrapping in cooperation with DLR and has been funded by the DARA (project
„Rapid", grant number 50 EE 9431) which is greatly appreciated.
Keywords: ESA European Space Agency - Agence spatiale europeenne, observation de la terre, earth observation, satellite remote sensing, teledetection, geophysique, altimetrie, radar, chimique atmospherique, geophysics, altimetry, radar, atmospheric chemistry
|Copyright 2000 - European Space Agency. All rights reserved.|