next up previous [pdf]

Next: Estimating the inverse gradient Up: Model fitting by least Previous: Gulf of Mexico Salt

VESUVIUS PHASE UNWRAPPING

Figure 11 shows radar images of Mt. Vesuvius[*]in Italy. These images are made from backscatter signals $s_1(t)$ and $s_2(t)$, recorded along two satellite orbits 800-km high and 54-m apart. The signals are very high frequency (the radar wavelength being 2.7 cm). The signals were Fourier transformed and one multiplied by the complex conjugate of the other, getting the product $Z=S_1(\omega) \bar S_2(\omega)$. The product's amplitude and phase are shown in Figure 11. Examining the data, you can notice that where the signals are strongest (darkest on the left), the phase (on the right) is the most spatially consistent. Pixel by pixel evaluation with the two frames in a movie program shows that there are a few somewhat large local amplitudes (clipped in Figure 11) but because these generally have spatially consistent phase, I would not describe the data as containing noise bursts.

vesuvio
vesuvio
Figure 11.
Radar image of Mt. Vesuvius. Left is the amplitude $\vert Z(x,y)\vert$. Nonreflecting ocean in upper-left corner. Right is the phase $\arctan( \Re Z(x,y), \Im Z(x,y)\;)$. (European Space Agency via Umberto Spagnolini)
[pdf] [png] [scons]

To reduce the time needed for analysis and printing, I reduced the data size two different ways, by decimation and local averaging, as shown in Figure 12. The decimation was to about 1 part in 9 on each axis, and the local averaging was done in $9\times 9$ windows giving the same spatial resolution in each case. The local averaging was done independently in the plane of the real part and the plane of the imaginary part. On the smoothed data the phase is less noisy.

squeeze
squeeze
Figure 12.
Phase based on decimated data (left) and smoothed data (right).
[pdf] [png] [scons]

From Figures 11 and 12, we see that contours of constant phase appear to be contours of constant altitude; this conclusion leads us to suppose that a study of radar theory would lead us to a relation like $Z(x,y)=e^{ih(x,y)}$, where $h(x,y)$ is altitude. We nonradar specialists often think of phase in $e^{i\phi} = e^{i\omega t_0(x,y)}$ as being caused by some time delay and being defined for some constant frequency $\omega$. Knowledge of this $\omega$ (as well as some angle parameters) would define the physical units of $h(x,y)$.

Because the flat land away from the mountain is all at the same phase (as is the altitude), the distance as revealed by the phase does not represent the distance from the ground to the satellite viewer. We are accustomed to measuring altitude along a vertical line to a datum; but here, the distance seems to be measured from the ground along a $23^\circ$ angle from the vertical to a datum at the satellite height.

Phase is a troublesome measurement, because we generally see it modulo $2\pi$. Marching up the mountain, we see the phase getting lighter and lighter until it suddenly jumps to black, which then continues to lighten as we continue up the mountain to the next jump. Let us undertake to compute the phase, including all its jumps of $2\pi$. Begin with a complex number $Z$ representing the complex-valued image at any location in the $(x , y )$-plane.

$\displaystyle r e^{i \phi}$ $\textstyle =$ $\displaystyle Z$ (84)
$\displaystyle \ln \vert r\vert + i \phi$ $\textstyle =$ $\displaystyle \ln Z$ (85)
$\displaystyle \phi(x,y)$ $\textstyle =$ $\displaystyle \Im \ln Z(x,y) ~+~ 2\pi N(x,y)$ (86)

Computers find the imaginary part of the logarithm with the arctan function of two arguments, atan2(y,x), which puts the phase in the range $-\pi < \phi \le \pi$, although any multiple of $2\pi$ could be added. We seem to escape the $2\pi N$ phase ambiguity by differentiating:
$\displaystyle {\partial\phi \over \partial x} \eq \Im {1 \over Z}{\partial Z \over \partial x} \eq
{\Im \bar Z {\partial Z \over \partial x} \over \bar Z Z }$     (87)

For every point on the $y$-axis, equation (87) is a differential equation on the $x$-axis. We could integrate them all to find $\phi(x,y)$. That sounds easy. On the other hand, the same equations are valid when $x$ and $y$ are interchanged, therefore we get twice as many equations as unknowns. Ideally either of these sets of equations is equivalent to the other; but for real data, we expect to be fitting this fitting goal:
\begin{displaymath}
\nabla \phi \quad \approx \quad {\Im \bar Z \nabla Z \over \bar Z Z}
\end{displaymath} (88)

where $\nabla = ({\partial \over \partial x}, {\partial \over \partial y} ) $. Mathematically, computing phase this way is like our previous seismic flattening with $\nabla \tau \approx {\bf d}$. Taking measurements to be phase differences between neighboring mesh points, it is more correct to interpret Equation (88) as a difference equation than a differential equation. Because we measure phase differences only over tiny distances (one pixel), we hope not to worry about phases greater than $2\pi$. But, if such jumps do occur, the jumps contribute to overall error.

Let us consider a typical location in the $(x , y )$ plane where the complex numbers $Z_{i,j}$ are given. Define a shorthand $a , b, c$, and $d$ as follows:


\begin{displaymath}
\left[
\begin{array}{ll}
a & b \\
c & d
\end{array} \r...
... & Z_{i,j+1} \\
Z_{i+1,j} & Z_{i+1,j+1}
\end{array} \right]
\end{displaymath} (89)

With this shorthand, the difference equation representation of the fitting goal (88) is:
\begin{displaymath}
\begin{array}{rcl}
\phi_{i+1,j} -\phi_{i,j} &\approx& \Del...
...\phi_{i,j+1} -\phi_{i,j} &\approx& \Delta\phi_{ab}
\end{array}\end{displaymath} (90)

Now, let us find the phase jumps between the various locations. Complex numbers $a$ and $b$ may be expressed in polar form, say $a=r_ae^{i\phi_a}$ and $b=r_be^{i\phi_b}$. The complex number $\bar a b = r_a r_b e^{i(\phi_b-\phi_a)}$ has the desired phase $\Delta \phi_{ab}$. To obtain it we take the imaginary part of the complex logarithm $\ln \vert r_a r_b\vert + i\Delta \phi_{ab}$:
\begin{displaymath}
\begin{array}{lllll}
\phi_b-\phi_a &=& \Delta \phi_{ab} &=...
...d-\phi_b &=& \Delta \phi_{bd} &=& \Im \ln \bar b d
\end{array}\end{displaymath} (91)

which gives the information needed to fill in the right side of (90), as done by subroutine igrad2init() from module unwrap [*].

The operator needed is igrad2, gradient with its adjoint, the divergence.



Subsections
next up previous [pdf]

Next: Estimating the inverse gradient Up: Model fitting by least Previous: Gulf of Mexico Salt

2014-12-01