VESUVIUS PHASE UNWRAPPING

Figure 2 shows radar images of Mt. Vesuvius in Italy. These images are made from backscatter signals s₁(t) and s₂(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). They 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 2. 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 2) but because these generally have spatially consistent phase I would not describe the data as containing noise bursts.

 

vesuvio90
Figure 2 Radar image of Mt. Vesuvius. Left is the amplitude. Non-reflecting ocean in upper left corner. Right is the phase. (Umberto Spagnolini)

To reduce the time needed for analysis and printing, I reduced the data size two different ways, by decimation and by local averaging, as shown in Figure 3. 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. Putting them back together again showed that the phase angle of the averaged data behaves much more consistently. This adds evidence that the data is not troubled by noise bursts.

 

squeeze90
Figure 3 Phase based on decimated data (left) and smoothed data (right).

From Figures 2 and 3 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=e^ih where h is altitude (in units unknown to us nonspecialists). 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 of its jumps of $2\pi$.Begin with a complex number Z representing the complex-valued image at any location in the (x,y)-plane.

$$\begin{eqnarray} r e^{i \phi} &=& Z \\ \ln \vert r\vert + i \phi &=& \ln Z \\ \phi &=& \Im\ \ln Z \end{eqnarray}$$ (68) (69) (70)

A computer will find the imaginary part of the logarithm with the arctan function of two arguments, atan2(y,x), which will put 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:

$$\begin{eqnarray} {\partial\phi \over \partial x}&=& \Im\ {1 \over Z}{\partial Z ... ...x}&=& {\Im\ \bar Z {\partial Z \over \partial x} \over \bar Z Z }\end{eqnarray}$$ (71) (72)

For every point on the y-axis, equation (72) is a differential equation on the x-axis, and 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, so we get twice as many equations as unknowns. For ideal data, either of these sets of equations should be equivalent to the other, but for real data we expect to be fitting the fitting goal  

$$\nabla \phi \quad \approx \quad {\Im\ \bar Z \nabla Z \over \bar Z Z}$$ (73)

where $\nabla = ({\partial \over \partial x}, {\partial \over \partial y} )$.

We will be handling the differential equation as a difference equation using an exact representation on the data mesh. By working with the phase difference of neighboring data values, we do not have to worry about phases greater than $2\pi$(except where phase jumps that much between mesh points). Thus we solve (73) with finite differences instead of differentials. Module igrad2 is a linear operator for the difference representation of the gradient operator.  

#$=head1 NAME
#$
#$igrad2 - 2D gradient operator 
#$  
#$=head1 SYNOPSIS
#$
#$C<ierr=igrad2_init(n1,n2)>
#$      
#$C<ierr=igrad2_lop(adj,add,xx,yy)>
#$              
#$=head1 PARAMETERS
#$              
#$=over 4
#$
#$=item n1,n2  - integer   
#$
#$      dimensions of dataset
#$     
#$=item adj,add,xx,yy -
#$
#$      standard operator interface
#$  
#$=back
#$
#$=head1 DESCRIPTION
#$
#$ Calculate 2-D gradient
#$
#$=head1 SEE ALSO
#$
#$L<igrad1>
#$
#$=head1 LIBRARY
#$  
#$B<geef90>
#$ 
#$=cut
#$

module igrad2 {                       # 2-D gradient with adjoint,  r= grad( p)
integer :: n1, n2
#%_init   (n1, n2)
#%_lop ( p(n1, n2),  r(n1,n2,2))
integer i,j
do i= 1, n1-1 {
do j= 1, n2-1 {
	if( adj) {
		p(i+1,j  )  +=  r(i,j,1) 
                p(i  ,j  )  -=  r(i,j,1) 
                p(i  ,j+1)  +=  r(i,j,2) 
                p(i  ,j  )  -=  r(i,j,2)
		}
	else {  r(i,j,1)  +=  ( p(i+1,j) - p(i,j)) 
                r(i,j,2)  +=  ( p(i,j+1) - p(i,j)) 
		}
	}}
}

To do the least-squares fitting (73) we pass the igrad2 module to the Krylov subspace solver. (Other people might prepare a matrix and give it to Matlab.)

The difference equation representation of the fitting goal (73) is:  

$$\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}$$ (74)

where we still need to define the right-hand side. Define the parameters a, b, c, and d as follows:

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

Arbitrary 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{array} {lllll} \phi_b-\phi_a &=& \Delta \phi_{ab} &=... ...d-\phi_b &=& \Delta \phi_{bd} &=& \Im \ln \bar b d \end{array}$$ (76)

which gives the information needed to fill in the right-hand side of (74), as done by subroutine igrad2init() from module unwrap .