Flat-earth phase removal algorithm improved with frequency information of interferogram a. Bin Ai a, Kai Liu b∗, Xia Li a, Dong hai, Lia School of Geography science and Planning, Sun Yat-sen University, 135 West Xingang Road, Guangzhou, P.R. China, 510275, b. Guangzhou Institute of Geography, Guangzhou, P.R. China, 510070 ABSTRACT As is well known, Interferometric synthetic aperture radar (InSAR) has been widely used in remote sensing field, which can reflect actual topographic trend or possible surface deformation. The precision of interferometric phase is critical to the final measurement. Due to the orbit attitude influence, such phase difference between the scattering elements on the same height level, which is named as flat-earth phase, usually causes the complex interferogram dense and difficult to be used in further procedures. Before phase unwrapping, interferogram must be flattened to derive accurate topographic or deformation information. Traditional methods pose problem to retrieve accurate flat-earth phase, which finally lead to inaccurate elevation or deformation information. A new algorithm of flat-earth phase removal is proposed in this paper based on SAR satellite system geometry and spectrum information of actual interferogram. The basic procedure of the method is firstly introduced and then the test results are following listed. From the comparison between the new algorithm and conventional ones, some advantages can be easily shown: The whole calculation can be easily understood and applied; accurate flat-earth phase can be retrieved and removed not only airborne SAR image but also satellite SAR image, which will improves the quality of complex interferogram to be well unwrapped; Keywords: Flat-earth effect; Complex Interferogram; interferometric synthetic aperture radar (InSAR); Frequency shift 1. Introduction Now InSAR has been an important hotspot and widely used in radar imaging area for its high ability to retrieve elevation and surface deformation information, which can reflect such information according to the phase difference between two or more images. InSAR has been proved to be applicable to many fields such as surveying topography, estimating ocean’s currents, detecting and locating moving targets, and so on. Basic procedures of InSAR processing can be concluded as follows [1]: selecting an image couple, co-registering the images, generating interferogram, flattening interferogram, filtering noise, unwrapping phase, reconstructing and correcting digital elevation model (DEM). But during the data processing, for the imaging geometry, the phase difference will exist between the two scattering elements locating on the same height, and phases in the interferogram often parallel to the azimuth or slant aspect, which is defined as flat-earth phase effect. It is very important to remove flat-earth phase, which usually increases phase density and causes the complexity of unwrapping [2, 3, 4]. And accuracy of flat-earth removal will influence the further process; especially make the phase to DEM transformation more difficult for the nonlinear function and Ground Control Points (GCP) requirement. Main goal of flattening is to wipe out the flat earth phase which has nothing with topographic elevation or surface deformation [4]. In the past years, many researchers have done much more on the flattening algorithm; and many effective algorithms are also proved to well solve the flat earth removal problem. Conventional methods can be mainly summarized into categories such as: one is based on the satellite orbit information, which calculates the flat-earth phase of each pixel according to the geometry between the two satellites during the image survey [5]. But it is most based on assumption baseline not varying with time, the calculation errors often cause the inaccuracy in flat-earth phase removal, and also it is difficult to realize the evaluation for the complex calculation; The other one uses the coarse external DEM to remove flat-earth phase, which is also difficult to carry out for the lack of external DEM; Another common algorithm is to use the frequency information of original interforgram [2]. When flat-earth phase is wrapped in the interferogram, ∗ Corresponding author, Email: liukai@gdas.ac.cn. Geoinformatics 2008 and Joint Conference on GIS and Built Environment: Classification of Remote Sensing Images, edited by Lin Liu, Xia Li, Kai Liu, Xinchang Zhang, Proc. of SPIE Vol. 7147, 71471A © 2008 SPIE · CCC code: 0277-786X/08/$18 · doi: 10.1117/12.813247 Proc. of SPIE Vol. 7147 71471A-1 the spectrum peak will not appear on the zero-frequency. However, common frequency-shift method only transforms the spectrum peak to zero frequency point constrainedly without considering actual characters of interferogram. Generally speaking, traditional methods cannot flatten the interferogram accurately enough, though the density of fringes can be remarkably decreased, which makes the subsequent phase filtering and unwrapping convenient, the inaccurate flat-earth phase is then added back into the interferogram before DEM reconstruction [1], [4]. It obviously wastes time and resource. An improved flattening algorithm is proposed in the view of spectrum information and furthered tested with the sample data. This paper is organized as follows. In Section I, the general analysis of flat-earth effect in InSAR is discussed. Then the problems of common frequency shift flattening methods are analyzed in Section II-A, some improvements on common methods are further provided in Section III-B. In Section IV, the improved algorithm is tested by InSAR images. One is Shuttle Radar Topography Mission (SRTM) data imaged from a mountain area in Ötztal / Austria. The other is JERS-1 data imaged from a volcano from volcano island Unzen in Japan. Compared to common methods, the flattening result shows that sidelong fringes are basically eliminated and interferometric fringes are clear and consistent with topography; finally, some conclusions are drawn in section 2. Principle of Flat-earth Phase SAR images are acquired at slant range direction; even different points without any topography variation, there will still be regular fringes appeared in interferogram for the range difference [4, 9]. In the case of repeat-pass interferometry, the orbits usually show yaw angle difference during imaging on different time, which will result some regular phase stripes on azimuth aspect. Such data as JERS-1 and Radarsat, unstable orbits leads to flat-earth also appearing in the azimuth direction. It is well known the original interferogram is composed of several different parts including topographic phase, the possible deformation phase, the atmosphere, the track error and so on [6]. When there is no deformation, original phase difference can be expressed as: φ = φtopo + φ flat + φnoise (1) If surface deformation takes place during the imaging times span, original interferometric phase will be characterized as[7, : 8] φ = φdeformation + φtopo + φ flat + φatm + φnoise (2) The purpose of flattening is to remove φ flat in the interferogram. To explain the formation of flat-earth phase, repeat-pass imaging geometry and spatial relationship is analyzed in this paper as shown in fig.1. Fig.1 Geometry of Flat-earth Phase Proc. of SPIE Vol. 7147 71471A-2 Where A1, A2 represents antenna’s positions of two times imaging respectively, the distance between the two antennas is defined as baseline B, α represents the angle between baseline and horizontal surface, P1 and P2 are the points on the same height level with different range distance to the imaging satellite, r11 and r12 are the distances between the two satellites and an arbitrary target on a terrain respectively, θ is the look angle, variation of incidence angle between P1 and P2 is equal to ∆θ , and h0 is the height of the platform. Different representations for the baseline are figure derived: 1) parallel B∥ and perpendicular B⊥; 2) horizontal Bh and vertical Bv; Consequently the phase of P1 and P2 can be given as follows: Where, λ φ=− 4π λ B sin(θ − α ) φ' = − 4π λ B sin(θ + ∆θ − α ) (3) is the wavelength of radar. The phase difference between P1 and P2 can further be given as ∆φR = − 4π B cos(θ − α )∆R 4πB⊥ ∆R =− λ r11 tan θ λr11 tan θ (4) From formula (3), it can be understood that phase difference between the two points on the same height level is related to range distance difference, and when the orbits are not parallel to each other or big slope angle appearing in the imaging area, regular flat-earth phase parallel to the range aspect will also expose to the interferogram. Flat-earth effect cannot be ignored during InSAR processing.. As is discussed above, though without any altitude changes between the two points P1 and P2, regular stripes will still expose in the interferometric phase due to the different look angles during the repeat-pass imaging, And also due to the orbit error, the flat-earth effect will appear both in the range and azimuth direction [4]. The dense stripes caused by flat-earth effect can overshadow the phase variety of topographical changes, and moreover bring much burden for subsequent phase filtering and unwrapping. 3. Flatten Algorithm Based on improved Spectrum Information of Interferogram In the view of image information, frequency spectrum of interferogram is characterized by some typical shape, which generated by flat-earth phase is generally sharp and narrow. Then the whole frequency spectrum peak in original interferogram will not centralize to zero-frequency. This information can provide ideas for flat-earth removal, and many researches proposed different algorithms. The flattened interferogram φtopo will be basically consistent with the shape of topography which does not vary greatly. So the spectrum information of a flattened cint will always produce a sharper pulse at the zero-frequency position [2]. Frequency shift method is commonly used for its convenience without geometry parameters. 3.1 Problems of Conventional Frequency Shift methods To remove flat-earth phase, common frequency-shift method usually locates the maximum fringe frequency in the range direction where the flat-earth phase dominates, and removes it with corresponding linear phase compensation in the time domain (by complex multiplication) or the frequency domain (by circular shift). Although relatively simple, conventional frequency-shift method often fails to remove flat-earth phase accurately for the reasons following: only integer samples of frequency deviation can be easily detected and compensated to the zero frequency. But when the deviation shows decimal numbers, the residual flat earth phase is neglected and will still remain in the interferogram; this method usually make some assumption that fringes variation show the same trend in the whole interferogram, which is often not corresponding to actual situation, especially in cross-track interferometry, spanning of flat-earth phase at far range is wider than that in near range. 3.2 An improvement on Conventional Frequency Shift Algorithm In order to overcome the limitations above, an improved algorithm is proposed in this paper. Three key modules are discussed and implemented. (1)Detection of non-integer frequency deviation Actually, frequency deviation caused by flat-earth phase can be divided into inter part and decimal part as formula (5): Proc. of SPIE Vol. 7147 71471A-3 φ flat = ( N d + N i )2π (5) Where, N d is the integer deviation part, N i represents decimal part. Since the spectrum of original complex interferogram is distributed at discrete samples by FFT technique, only the integer frequency deviation N d can be detected by conventional spectrum analysis. To solve this problem, the cubic spline function is here introduced to inspect the spectrum with more details, which can be realized with Matlab toolbox. For a spectrum array with length of N, it is interpolated as Cubic spline data interpolation. So the frequency spectrum is further discreted, finer spectrum information can be measured and processed. (2)Block division To adjust our method referring to the different phase density variation along the slant range, original interferogram is divided into several blocks in both range and azimuth direction. Each block can overlap to a conjoint one. This technique would be effective when the parallel base line changes with time. (3)Evaluation Index Next, the definition of entropy is also introduced here to determine the exact frequency deviation. As we know, in most cases, if the flat-earth phase can be removed from the interferogram accurately, the remaining interferogram contains a minimum number of strips. Entropy can represent the strip density of the interferogram, by definition as formula (6): H = ∑ - Pi log Pi (6) i Where, Pi denotes the probability of different phase value. For every test frequency deviation, the value of entropy is calculated. When the entropy achieves a minimum value, the flat-earth phase is absolutely removed. According to the analysis above, the flow chart of the flattening algorithm can be realized on the following aspects: 1) Firstly, original interferogram is divided into several blocks in azimuth aspect, the spectral information info1×N of a complex interferogram cintM×N (M-point range pixels and N-point azimuth lines) in every block is correspondingly defined as: Info1× N = ∑ Block _ azi i =1 FFT (cint1× N ) (7) Where, the FFT operator (along each range line) means 1-D fast Fourier transformation and the Σ operator (cross each azimuth line) avoids the influence by noises. Detect the integer frequency deviation of (7) 2) Interpolate (7) and search a more detailed value (non-integer) 3) Calculate the entropy. If it is a maximum value, then continue; else back to 2) 4) Shift the Interpolated interferogram, and resample it like the original one. The spectral information will always produce a sharp pulse at the zero-frequency when the flat-earth phase is absolutely removed from the interferogram. The whole procedure is listed in fig 2 Proc. of SPIE Vol. 7147 71471A-4 Initial interferogram and geometry parameters Calculating frequency spectrum info of cint Resampling frequency spectrum e.g. 0.01 sample Blocking in azimuth Blocking in range Obtaining the deviation between the peak and zero frequency Obtaining the deviation between the peak and zero frequency Moving gradually Moving gradually No No 0 freq position? 逐步 移动 Entropy min? 0 freq position? Entropy min? Yes Yes Calculating flat earth phase in azimuth 是 Calculating flat earth phase in azimuth Removal of flat earth phase Fig.2 Flow Chart of the Improved Algorithm 4. Application In this section, two pairs of sample data are used to validate the algorithm. The results are furthermore compared with those derived from conventional methods. One test data pair is SRTM data, the baseline of which is very short, only tens of meters. The other one is InSAR pair of Mt. Unzen, Japan, the SAR data is ©METI/NASDA (1992-1993), which retains ownership of the JERS-1 data. The baseline is relatively longer, up to hundreds of meters. 4.1 Results of SRTM Test Data The Amplitude of SLC is shown as fig.3 (a), which can well display the topographic trend of the imaging area. After coregistering the slave to the master, original interferogram with size of 2048*2048 is formed with the complex matrix conjugation as shown in fig 3(b). There exists significant flat Earth phase. The dense phase stripes are regular parallel to some aspect. According to frequency spectrum curve analysis, the peak is obviously not located on zero frequency. Actual topographic trend is obviously shaded for the flat-earth phase, which will burden the following process such as unwrapping and phase to height transformation. Proc. of SPIE Vol. 7147 71471A-5 Raw Interferogram of SRTM Data 0 3 200 2 400 600 1 Azimuth 800 1000 0 1200 -1 1400 1600 -2 1800 -3 2000 0 Fig.3 Amplitude of SLC 200 400 600 800 1000 1200 1400 1600 1800 2000 Slant Range Fig.3 (b) Raw Interferogram (1) Result of flat-earth removal with conventional frequency moving method Flat-earth phase will be estimated and removed by the circular shift principle within common frequency shift method. The results are shown as fig.4 (a) and fig.4 (b). The method takes the whole interferogram as equal pattern from the near range to the far range, which makes the assumption of fixed parallel baseline. The flat-earth phase stripes are evenly distributed unlike the pattern of raw interferogram. Therefore, the method succeeds in retrieving accurate flat-earth phase in the middle of imaging area but fails at the far range and near range area. And only the integer numbers of difference between the peak and the zero frequency point can be moved. Flat-Earth Phase Derived with Conventional Frequency Shift Method Interferogram Flattened with Conventional Frequency Shift Method 3 3 200 50 2 400 2 100 150 600 1 1 200 1000 0 1200 -1 1400 Azimuth Azimuth 800 250 0 300 -1 350 400 1600 -2 1800 -2 450 2000 -3 200 400 600 800 1000 1200 1400 1600 1800 2000 500 -3 50 100 Slant Range 150 200 250 300 350 400 450 500 Slant Range Fig.4 (a) Flat-earth phase Derived with Common Method Fig.4 (b) Interferogram Flattened with Common Method (2)Result of flat-earth removal with the method proposed in the paper Based on the analysis of frequency spectrum information, the interferogram is then flattened with the technique proposed in this paper. Firstly, different sizes of blocks along the azimuth direction are divided; frequency spectrum information is then derived with formula (7), which is further interpolated with the precision of 0.1 samples; with the evaluating criterion as listed in formula (6), appropriate deviation between the peak and zero frequency will be well assessed after iterative calculation. Results are listed as fig.5. Therein, the whole flat-earth phase is shown as fig.5 (a), and the corresponding flattened interferogram is shown as fig.5 (b). From the figures shown below, the flat-earth phase, whose frequency spectrum shows quite sharp shape, is unequal in the raw interferogrram; the stripes are dense in the near range, while sparse in the far range, which is similar with the Proc. of SPIE Vol. 7147 71471A-6 pattern of raw interferogram. The interferogram flattened then matches the SLC image very well, which can be further used to retrieve the actual topographic trend, though the value of which is reversed, it won’t affect the further processing; and the peak of the frequency spectrum is well moved to the zero frequency point. It shows that the method proposed can be well used to flatten interferogram with unequal stripes. Compared the results derived from the two algorithms with the raw interferogram and SLC image, big difference can be easily found. Nearly even interval between the flat-earth phase stripes is exposed to us as shown in fig.4 (a), which is obviously inconsistent with the pattern of raw interferogram. Much more residual flat-earth phase is still remained especially in the far range as shown in fig.4 (b). Inaccurate flat earth will be added back to the phase unwrapped and ground reference points will be required for phase to height transformation. While, the results obtained by the method improved conform to the SLC image well, the uneven interval between the flat earth phases demonstrates the variation of parallel baseline during imaging survey, and the relation between the flattened interferogram and topography phase can be distinctly built, which will make the following procedures easier. Flat-Earth Phase Flattened and Multilooked Interferogram 3 3 200 50 2 600 1 Azimuth 800 1000 2 100 Azimuth 400 150 1 200 250 0 1200 0 300 -1 1400 -1 350 400 1600 -2 -2 1800 450 2000 200 400 600 500 -3 800 1000 1200 1400 1600 1800 2000 Slant Range Fig.5 (a) Flat-Earth Phase in the Interferogram -3 50 100 150 200 250 300 350 400 450 500 Slant Range Fig.5 (b) Interferogram Flattened and Multilooked 4.2 Results of JERS-1 Test Data The method is also tested with the other satellite JERS-1 InSAR image pair. Amplitude image is displayed as fig.6 (a), and Raw interferogram is also shown as fig.6 (b). From the figure, flat-earth phases increase the density of the stripes in the interferogram, little useful information can be derived. Unlike the SRTM data, flat-earth phase appears both in the slant range direction and azimuth direction, which will make the flattening more complex. Furthermore, it shows that the peak of the frequency spectrum, showing unsymmetrical shape, locates close to coordinate 700 in the slant range direction clearly deviating from the zero frequency point. Raw Interferogram 3 2 500 1 Azimuth 1000 0 1500 -1 2000 -2 2500 -3 200 400 600 800 1000 1200 1400 Slant Range Fig.6 (a) Amplitude of SLC Fig.6 (b) Raw interferogram of JERS-1 data Proc. of SPIE Vol. 7147 71471A-7 (1) Result of flat-earth removal with conventional frequency moving method Firstly, to compare the method proposed in this paper with the conventional methods, results derived from common frequency moving method are listed in fig.7. As is shown, there are still much more residual flat-earth phase, especially under the condition of flat-earth phase appearing both in the slant range and azimuth direction, which reflects the common method cannot remove that from the interferogram accurately. Interferogram Flattened with Common Technique 3 100 2 200 Azimuth 300 1 400 0 500 600 -1 700 -2 800 -3 900 200 400 600 800 1000 1200 1400 Slant Range Fig.7 Interferogram Flattened with Common Method (Near range is on the left) (2) Result of flat-earth removal with the method proposed in the paper Then, conformed to the procedures of SRTM test data, flat-earth phases in the slant range and azimuth direction are respectively wiped off. Firstly, several blocks are divided in the azimuth based on the frequency spectrum information of raw interferogram. According to the procedures listed in the fig.2, flat-earth phases in the slant range direction can be well derived. The same process to the azimuth direction is then carried out, only dividing the interferogram without flatearth phases in the slant range.Accurate results are also obtained as shown in fig.7, which can further prove that it is more applicable into the flat-earth phase removal field. Flat-Earth Phase Obtained with Improved Method Interferogram Flattened with the improved algorithm 3 3 100 2 Azimuth 1 1000 2 200 300 1 Azimuth 500 400 0 1500 0 500 -1 600 -1 2000 700 -2 -2 800 2500 -3 200 400 600 800 1000 1200 1400 -3 900 200 Slant Range Fig8 (a) Flat-earth Obtained with the Improved Technique 400 600 800 1000 1200 1400 Slant Range Fig.8 (b) Corresponding Flattened Interferogram Proc. of SPIE Vol. 7147 71471A-8 Obviously, the pattern of flat-earth phase is also consistent with the raw interferogram, distributing regularly and densely as shown in fig.8 (a). The density of stripes decreases distinctly with the improved flattening algorithm, which basically conforms to the actual satiation. From the results shown above, difference can be easily seen between the two methods. There are still much more flatearth phase in the interferogram after calculation with the conventional technique, especially in the far range for short baseline data such as SRTM image or in the azimuth aspect for relative long baseline such as JERS-1data. However, the phases in the interferogram can basically demonstrate the trend of actual topography with the method improved in this paper, flat-earth phases are thoroughly removed from original interferogram between the master and slave image, which can well move the peak of frequency spectrum to the zero frequency point consistent with conventional pattern. 5. Conclusion and Discussion The principle of flat-earth phase is discussed, and the generally flattening techniques are also summarized in this paper. Two categories of flattening algorithms have been introduced in section V. The first type determines the flat-earth phase in terms of geometrical parameters, otherwise fit the expression by a polynomial from known DEM. However in many cases such reliable parameters or reference DEM are usually absent. The second type computes the flat-earth phase by measuring the dominant fringe frequency of the interferogram and introducing into the corresponding linear phase compensation. This algorithm may be more reasonable if it considered the following two: one is that only an integer of dominant fringes can be measured, and the other is that the linear phase compensation has comparatively big errors corresponding to the nonlinearity in formula (2). Obviously the first type is only dependent on the external information without the assistance from the interferogram itself, while the second type is in reverse. Aiming to the limitation of common frequency shift technique, an improved algorithm is then proposed based on the frequency spectrum information, which is further validated with two pairs of sample data. It is shown that the method proposed in this paper can wipe off flat-earth phase accurately, not only for short baseline data but for relative long baseline InSAR pair, more effective than the conventional ones, which will make phase unwrapping more easily. Compared with conventional frequency shift methods, the method can be concluded as: (1) Only information of original interferogram is required, which can avoid some errors for lack of parameters or inaccurate external referring information; (2) The blocks can avoid the same trend of flat-earth phase in the raw interferogram, which is more consistent with imaging geometry, especially in cross-track interferometry. (3) The interpolation and iterative process can obtain more detail information, and can better cope with the situation moving decimal sampling points to zero frequency point. ACKNOWLEDGEMENT This work was support by Guangdong Academy of sciences outstanding young (2007.12). REFERENCES 1 J. Moreira, M. Schwäbisch, G. Fornaro, R. et al, (1995). “X-SAR interferometry: First results,” IEEE Trans. Geosci. Remote Sensing, 33 (4): 950-956. 2 Yang Lei, Zhao Yongjun1, Wang Zhigang. (2004). Frequency Shifted - based Approach for InSAR Flat Earth Effect Removal. JOURNAL OF ELECTRONIC MEASUREMENT AND INSTRUMENT, 18 (4): 15-20. (in Chinese) 3 ZHANG Yanjie. (2005). A Comparison of the Different Models Used For Interferograms Flattening. International Geoscience and Remote Sensing Symposium (IGARSS) 5: 5194-5196. 4 ZENG qiming, LI xiaofan, GAO liang, (2006) An Improvement to Flattening in Interferometric SAR Processing. Proc. of SPIE. 6200: 62000D-1—62000D-6. 5 K. Ren, V. Prinet, and X. Shi, (2003). “Spaceborne InSAR flat earth removal processing based on geolocation method,” Geomatics and Information Science of Wuhan University, 28 (3): 326-329. (in Chinese) Proc. of SPIE Vol. 7147 71471A-9 6 YU Jingtao, CHEN Ying, (2005). “A Discussion on Some Key Techniques of INSAR Data Processing”, Remote Sensing Information, 2005.1: 24~27. 7 F. Gatelli, A.M. Guarnieri, F. Parizzi, et al, (1994). “The wavenumber shift in SAR interferometry,” IEEE Trans. Geosci. Remote Sensing, 32 (4):. 855-865. 8 R. Siegmund, M. Bao S. Lehner, et al, (2004). “First demonstration of surface currents imaged by hybrid along- and cross-track interferometric SAR,” IEEE Trans. Geosci. Remote Sensing, 42 (5): 511-519. 9 D. Geudtner, M. Schwabisch (1996). An algorithm for precise reconstruction of InSAR imaging geometry: Application to “flat earth” phase removal, phase-to-height conversion, and geocoding of InSAR-derived DEMs[J]. EUSAR’96, Konigswinter, Germany: 249-252. Proc. of SPIE Vol. 7147 71471A-10