Model fitting by least squares

VESUVIUS PHASE UNWRAPPING

Figure 10 shows radarimages of Mt. Vesuviusin 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). 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 10. 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 10) but because these generally have spatially consistent phase, I would not describe the data as containing noise bursts.

vesuvio
Figure 10. 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 11. 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.

squeeze
Figure 11. Phase based on decimated data (left) and smoothed data (right).

From Figures 10 and 11 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 non-radar-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 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.

$$r e^{i \phi} = Z$$ (88)

$$\ln \vert r\vert + i \phi = \ln Z$$ (89)

$$\phi(x,y) = \Im \ln Z(x,y) ~+~ 2\pi N(x,y)$$ (90)

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:

$${\partial\phi \over \partial x} = \Im {1 \over Z}{\partial Z \over \partial x}$$ (91)

$${\partial\phi \over \partial x} = {\Im \bar Z {\partial Z \over \partial x} \over \bar Z Z }$$ (92)

For every point on the $y$-axis, equation (92) 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}$$ (93)

where $\nabla = ({\partial \over \partial x}, {\partial \over \partial y} )$. This is essentially the same application we solved flattening seismic data with the regression $\nabla \tau \approx {\bf d}$. Taking measurements to be phase differences between neighboring mesh points, it is more correct to interpret equation 93 as a difference equation than a differential equation. Since 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, they will 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:

$$\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]$$ (94)

ith this shorthand, the difference equation representation of the fitting goal (93) 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}$$ (95)

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

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

Model fitting by least squares