2-D phase unwrapping

Nizar Chemingui, Robert G. Clapp, and Jon Claerbout

nizar@sep.stanford.edu

ABSTRACT

We present a method for unwrapping the phase of a satellite data set that maps Mt.Vesuvius in Italy. The method finds a solution for a partial differential equation by redefining the problem as a system of regressions that we solve using linear inversion techniques with iterated reweighting. We show three different techniques for eliminating the noise bursts from the raw data and fitting it to a linear model. The first process is a two-step solution in which the data is first cleaned up and then fitted to a model by phase unwrapping. In the second approach, starting from a raw noisy data, we solve the inversion problem as a system of regressions using iterative weighting based on the residual. The third technique attempts to resolve a discontinuity in the data by solving a different system of regressions with weights based on the model. In the three processes, the inversion successfully unwrapped the phase of the radar image and produced an altitude map of the volcano region. We were able to accelerate the convergence of the solution by a smart preconditioning substitution. This transformation changes the problem-formulation variable using a leaky integration operator as a preconditioner. The substitution reduced the number of iterations in the least square-inversion by two orders of magnitude and provided a solution after very few iterations. The method, however, did not help resolve a presumed discontinuity in the model, and we may have to determine an optimal function for the iterated reweighting.

INTRODUCTION

We were presented with two radar signals, s₁(x,y,t) and s₂(x,y,t), from two satellite orbits 54 meters apart and 800 km above the surface. The signals were Fourier transformed and one was multiplied by the complex conjugate of the other, leading to the power spectrum $S_1(\omega) \overline {S_2}(\omega)$, whose square root gives the amplitude function, while the phase can be computed from the mathematical definition of a complex number. The straight mathematical definition for the phase function leads to a solution of a wrapped phase rather than a true phase. As a result, when looking at the phase map in Figure original, we do not see a hill-like structure, but rather something resembling a contour map. The contours of constant phase appear to be contours of constant altitude, which leads us to suppose that a study of radar theory would lead to a linear relation between phase and altitude. Actually, the altitude would be measured along the viewing angle, which happens to be 23 degrees from the vertical. The viewing angle is such that the mountain blocks some of the landscape behind it. Therefore, in that region, the contour lines join together, suggesting a jump in the altitude function over the hidden terrain.

The problem of uncovering the true phase from the complex data is known as phase unwrapping, for which we present a solution based on Claerbout's derivations. The solution reformulates the problem as a system of regressions that can be solved using linear least-square inversion techniques. We use iterated reweighting to estimate weighting functions based on the residuals or on the updated model of the previous iteration. Weighting with residual-based functions is a method of despiking that handles noise bursts in the residual. Similarly, model-based weighting functions handle noise bursts in data from a model that changes rapidly because it is high frequency or has jumps of discontinuity; such is the case with our altitude function.

For such a two-dimensional problem, convergence speed represents some problems. Our model contains 700 by 700 mesh points; so, theoretically, convergence is achieved after about half a million iterations. Even though in practice we only use a few thousand iterations, the convergence is still considerably slow. One way to accelerate the convergence is to make a change of variables and recast the problem in terms of an intermediate quantity. In a second stage we evaluate the solution for the original problem as a function of the intermediate solution.

In the next sections we will show how we estimate the weighting functions in the regressions and how we chose the preconditioning operator, but first we start by defining the phase unwrapping problem.

 

original
Figure 1 Radar image of Mt. Vesuvius. The amplitude is on the left, and the phase is on the right. (Umberto Spagnolini)

PHASE UNWRAPPING

A complex number is defined as:

$$\begin{eqnarray} Z & = & a+ib =re^{i\phi} \EQNLABEL{complex} \\ r & = & \sqrt{a^2+b^2}\\ \omega & = & \tan^{-1}{\left( b \over a \right)}\end{eqnarray}$$ (1) (2) (3)

The quantity $\omega$ in (3) describes an angular phase, but to solve for the true phase, $\phi$, we must consider two additional equations:

$$\begin{eqnarray} Z & = & re^{i(\omega+n\pi)} \quad \quad \quad \quad {\rm where}... ...\ \omega & = & \phi - \pi*{\rm int}{\left( \phi \over \pi \right)}\end{eqnarray}$$ (4) (5) (6)

Claerbout has proposed a method for retrieving the true phase $\phi$ instead of $\omega$. Starting from the definition of a complex number as $Z=re^{i\phi}$, we write:

$$\begin{eqnarray} \ln \vert r\vert + i \phi & = & \ln Z \\ \phi & = & \Im\ \ln Z ... ...over \partial x} \over \bar Z Z } \EQNLABEL{RHS} \EQNLABEL{x-comp}\end{eqnarray}$$ (7) (8) (9) (10)

where $\partial\phi \over \partial x$ is the x-component of the phase gradient which can be written for the 2-D case as $\nabla \phi= \bf d$.We can now formulate this as a 2-D inversion problem:

$${\bf L} \phi \approx \bf d \EQNLABEL{inversion}$$ (11)

where $\bf d$ is the right hand side of equation (10), hereinafter referred to as ``data phase gradient'', i.e, the x-component of $\bf d$ can be written as:

$${\bf d}_x = \Im\ { \bar Z {\partial Z \over \partial x} \over \bar Z Z } \EQNLABEL{inverse}$$ (12)

L is the transform (gradient operator i.e, $\nabla = ({\partial \over \partial x}, {\partial \over \partial y} )$ ) in which we propose to present the derivatives diagonally on a cartesian mesh:

$$\bold L \quad =\quad {1 \over 4}\ \begin{array} {\vert c\ve... ... & +1 \\ \hline \end{array}\end{array}\right) \EQNLABEL{OPER}$$ (13)

Notice that our transform operator is different from the original definition of (). Claerbout defines his modeling operator as ${\bf L}=(\bar Z Z)\nabla$ by moving the quantity $\bar Z Z$ to the left hand side of inversion in order to avoid division by zero. In our particular data set, however, the problems resulted from very high values in the amplitude term $\bar Z Z$; therefore, we chose to define $\bf L$ as in OPER and write the data phase gradient in inverse as:

$$d \approx \Im\ { \bar Z {\partial Z \over \partial x} \over \bar Z Z + \epsilon} \EQNLABEL{inverse2}$$ (14)

where $\epsilon$ is a small positive quantity used as a stabilizer to avoid zero division.

MODEL FITTING

Given very noisy raw data, our biggest challenge was handling the noise spikes and bursts while attempting the phase unwrapping. We used three different techniques to clean up the data and unwrap its phase:

1.

Clean up the noisy data phase gradient, then fit it to a model

2.

Fit the noisy data phase gradient to a linear model

3.

Fit the noisy data phase gradient to a model with discontinuities

De-bursting with Iterated reweighting

A method for creating cleaned up data $\bar {\bold d}$from the noisy data phase gradient $\bold d$ is solving the regression

$$\begin{eqnarray} 0 &\approx & \bold W (\bar {\bold d} - \bold d) \\ 0 &\approx & \sigma \nabla^2 \bar {\bold d} \end{eqnarray}$$ (15) (16)

where $\bold W$ is a diagonal matrix with weights sprinkled along the diagonal, and where $\nabla^2$ is a matrix with a roughener like the Laplacian or the gradient operator. Results showed a performance for $0 \approx \sigma \nabla \bar {\bold d}$ similar to that $0 \approx \sigma \nabla^2 \bar {\bold d}$.Throughout the examples that follow, we chose to use $\nabla^2$in the second regression.

Roughly, we want a weight w_i=|r_i|^-1/2 for medians and w_i=|r_i|^-1 for modes. To avoid zero division and to handle small values in the usual L₂ way we use weights such as those used in Claerbout :

$$\begin{eqnarray} \EQNLABEL{wmed} w_i^{\rm (median)} &=& {1 \over \sqrt{1 + \vert... ...de} w_i^{\rm (mode)} &=& {1 \over 1 + \vert r_i/\overline{r}\vert}\end{eqnarray}$$ (17) (18)

where $\overline{r}$ is an average residual magnitude chosen to be

$$\overline{r} = {\rm median}(\vert r_i\vert)$$ (19)

Figures data and deb-data show the data phase gradient (right hand side of inverse2) before and after smoothing with (15) and (16). The results show a definite increase in continuity. In addition, the number of outliers is greatly reduced by this process.

Once the data phase gradient has been filtered using regressions (15) and (16), the problem is reduced to performing the phase unwrapping by the solution of ${\bf L} \phi \approx \bar {\bold d}$. Figure method1 shows the results of the inversion. An altitude map of the mountain is uncovered and outlines the shape of the volcano.

Fitting noisy data to a linear model

In this approach, we fit the noisy data phase gradient $\bold d$to a linear model, say $\bold d \approx \bold \nabla \bold \phi$,using the following regressions:

$$\begin{eqnarray} 0 &\approx & \bold W (\bold \nabla \bold \phi -\bold d ) \\ 0 &\approx & \sigma \nabla^2 \bold \phi \EQNLABEL{mtd2}\end{eqnarray}$$ (20) (21)

We needed to experiment with the free parameter $\sigma$ in mtd2. A good result was obtained for $\sigma=3$.Notice that this approach is very similar to the first method of cleaning up the data phase gradient and then applying the phase unwrapping solution. The two operations are here combined in a single system of regressions. Figure method2 shows the results of the inversion. The method succeeded in uncovering the hill shape of the volcano.

 

data
Figure 2 Data phase gradient before smoothing with regressions (15) and (16). Each panel represents one component of the gradient along one diagonal direction on a cartesian mesh.

 

deb-data
Figure 3 Cleaned data phase gradient using regressions (15) and (16). Each panel represents one component of the gradient along a diagonal direction.

 

method1
Figure 4 Altitude map obtained by integrating the data phase gradient after disposing of the outliers (method-1).

 

method2
Figure 5 Altitude map obtained by method-2 using a value of $\sigma=2$ in regression (21).

Fitting noisy data to a discontinuous model

The previous approach handled noise bursts in the residual, which is not the same as noise bursts in the data. We then attempted to fit the data phase gradient to a model with discontinuity jumps using the following regression:

$$\begin{eqnarray} 0 &\approx & \bold W_d (\bold \nabla \bold \phi - \bold d ) \\ 0 &\approx & \bold W_h \nabla^2 \bold \phi\end{eqnarray}$$ (22) (23)

where $\bold W_d$ is a residual-based weighting function and $\bold W_h$is a model-based weighting function. In our model the discontinuity is present in the gradient of the the altitude map; therefore, we seek a weighting function based on the first derivative of the model rather than the model itself. Roughly we want a weighting function $\bold W_h$ that is inversely proportional to the gradient of the phase. We use a model-based weighting function that is similar to the one employed to handle the noise bursts in the residual, that is,

$$W_i^{h} = {1 \over \sqrt{1 + \vert h_i/\overline{h}\vert}} \;$$ (24)

where $h_i=\nabla \bold \phi_{i}$ and $\overline{h}$ is an average gradient magnitude chosen to be

$$\overline{h} = {\rm median}(\vert h_i\vert)\; .$$ (25)

Figure method3 shows the solution for the phase map. The image is not significantly different from those obtained by the previous approaches. Some sharp features are present along the slopes of the mountain and may have arisen from the important effect of $\sigma$.

 

method3
Figure 6 Altitude map obtained by method-3 using weights based on the model in regression (23).

PRECONDITIONING BY INTEGRATION

Convergence speed was something of a problem in this data set; about 2000 iterations were required to achieve a solution for the phase map. Resampling to every other mesh point, although reduced the computational time, did not actually decrease the number of iterations. Later we used a change of variables to accelerate the iterative inversion. The preconditioning substitution changes the problem formulation variable $\phi$in inverse to a new variable $\bf x$. We used a leaky integration operator () as a preconditioner so that the new variable $\bf x$ is like the derivative of the solution phase $\phi$.By the preconditioning transformation $\phi=\bf B \bf x$, we have recast the original inversion relation inverse into

$${\bf L} \bf Bx \approx \bf d \EQNLABEL{recast}$$ (26)

After solving this system for $\bf x$, we merely evaluate $\phi=\bf B \bf x$to get the phase map, solution of the original problem. Since B is essentially integration whereas our transform $\bf L$ is a differentiation operator, we were able to converge to a solution in about 30 iterations and therefore reducing the computational time by two orders of magnitude.

DISCUSSION

Each of the procedures discussed above has yielded a final phase function that describes an altitude map of a volcano. Methods (1) and (2) are quite similar with the exception of data cleaning and model fitting being done in two stages in (1) by solving the regressions:

$$\begin{eqnarray} 0 &\approx & \bold W (\bar {\bold d} - \bold d) \\ 0 &\approx & \sigma \nabla^2 \bar {\bold d} \end{eqnarray}$$ (27) (28)

and then integrating the data phase gradient by solving the system

$${\bf L} \phi \approx \bar {\bold d} \; .$$ (29)

Method (2) combines the two stages in a single step and solves the system

$$\begin{eqnarray} 0 &\approx & \bold W (\bold \nabla \bold \phi -\bold d ) \\ 0 &\approx & \sigma \nabla^2 \bold \phi \EQNLABEL{mtd2}\end{eqnarray}$$ (30) (31)

Most differences between the images arise from the choice of the free parameter $\sigma$in the second regression. The preconditioning transformation provided very rapid convergence to the solution so that differences in computational time between the three approaches were not significant.

To better illustrate the altitude map around the hidden terrain behind the mountain, we plotted the function $\rm {tan} (a * \bold \phi)$. The plot should look similar to that of the initial phase map. The frequency factor, a, controls the number of contours used in the display. Figure display-grad-bw shows a contour map of the results of method 3. As we notice, the contours are more continuous than in the original phase map. The technique, however, did not resolve the presumed discontinuity in the hidden section of the mountain. Although the contours come very close together around that region, they do not actually join, but rather remained parallel.

The inability to image the discontinuity illustrates that there is a significant room for improvement. In order to understand why the discontinuity was not imaged and to come up with a scheme that is capable of solving for the jump in the altitude map, we need to test our technique on a synthetic model that simulates the nature of the discontinuity in the Vesuvius data set. To obtain the correct altitude map we must migrate the data around the fault instead of straight integration across it. Future goals of this effort will target testing the method on canonical models and optimizing the weighting functions used in the iterated reweighting.

 

display-grad-bw
Figure 7 Plot of tan(a*phase). Results of method-3

CONCLUSIONS

We have presented a technique for unwrapping the phase of a radar image of Mt. Vesuvius in Italy. The method finds a solution for a partial differential equation by redefining the problem as a system of regressions that we solve using iterative linear inversion techniques. The biggest challenge was handling the noise spikes and bursts in the raw data set. Three different techniques were used for cleaning up the data and fitting it to a model. The techniques use iterated reweighting by functions that are based on the residual or the model. The results of the inversion were quite accurate in uncovering the hill shape of the volcano. The convergence of the solution was accelerated by a change of variable using a leaky integration operator as a preconditioner. This substitution reduced the number of iterations by two orders of magnitude and provided a very rapid solution. The method, however, did not resolve for the presence of a discontinuity in the altitude map. We suspect that the sharpness of the image will be improved with the use of optimal functions in the iterated reweighting.

ACKNOWLEDGEMENTS We thank Umberto Spagnolini of Politecnico di Milano for providing us with the data. He recently submitted a paper on the subject to IEEE Trans. Geoscience and Remote Sensing.

[MISC,paper]