Detecting the Amount of Eroded and Deposited Sand Using DInSAR

DInSAR techniques can be used to monitor changes in desert height due to sandstorms. Hunshandake Sandy Land was used as the test area and was a main source of sand in the vicinity of Beijing. In order to study the source of sandstorm and its impact, a pair of EnviSat ASAR images for January 20 and October 26, 2004 were processed. After the image configuration, flat earth effect correction, data filtering, phase unwrapping and geo-coding, a deformation model over Hunshandake was built. The deformation in the test area ranges from -160 to +120 cm, and the height increased in most areas and decreased in a few locations, which is generally consistent with the sandstorm appearing at Beijing during 2004.


INTRODUCTION
Located in the middle of Xilinguole Plateau, Hunshandake Sandy Land, is among one of the 10 largest desert areas in China, and is one of the largest windy and sandy areas with semi-arid weather in the mid-temperate zone.The flux of sand and dust over Hunshandake was calculated using observed data from 2004 (Zhang et al. 2007).The fluxes of sand and dust under the non-dust weather condition, flying dust weather condition, and sandstorm weather condition were ±5, ±30, and -200 to +300 μg m -2 s -1 , respectively.Even in non-dust weather conditions, a certain amount of sand and dust are transported.The critical wind speed which can transport sand is 6 m s -1 (Zhang et al. 2007).In Hunshandake Sandy Land, the length of time with winds stronger than 6 m s -1 is about 1500 ~ 3381 h per year (Liu and Wang 2007), which implies that Hunshandake can act as a main source of blowing sand.In the past, Hunshandake was a sandy grassland with many lakes and rivers.But recently, grasslands, rivers, and lakes have shrunk, and the desert area has expanded due to the arid conditions and over-grazing.Especially in the past few decades, human activity has damaged the vegetation cover, and the desertification is becoming more serious.Up to now, Hunshandake Sandy Land has already been one of the main sources of sand transported to Beijing, China's capital city (Yue et al. 2008).With the rapid economic development, the deterioration of the environment has occurred and desertification has become more obvious, which damages efforts towards the sustainable development of China's economy and society (Wang et al. 2002;Li 2003).
Recently, both the government and the public in China have paid more attention to Hunsandake Sandy Land due to the increasing frequency of sandstorms in the northern China.The wind is frequent and strong, with a mean wind speed of about 3.5 to 5.3 m s -1 , and winds greater than 6 m s -1 occurring for 1500 to 3381 h annually (Liu and Wang 2007).Sand, which is distributed on the surface and soils can be eroded by wind.In response, the height of the desert surface changes much due to the frequent replacement of sand.According to the monitoring data for active dunes of Liu and Wang (2007), the height of active dune could decrease by about 70 cm between 1 October 2003 and 10 May 2004.
Complex image pairs contain information about amplitude and phase from synthetic aperture radar (SAR) observations.Only the amplitude information is used for ordinary SAR data processing, but the phase information is fully considered when using interferometric SAR (InSAR) techniques.Due to the phase difference between the two complex SAR images caused by the relative geometric relationship between the target and the two antennae, we can derive the interference fringes by the interference data processing, which are used to calculate the distance difference between the ground point and the two antennae along the oblique direction.The three-dimensional change of each point in the image can be measured precisely using the relationship among the sensor height, radar wavelength, beam direction, and antenna baseline (Massonnet and Feigl 1998).Differential InSAR (DInSAR) can process SAR image pairs over the same area with the mode of repeated orbits, to obtain the radar phase difference resulting from the change in distance between the ground surface and the satellite, and then compute the distance change of each point.Gabriel and Goldstein (1988) discussed the theory of centimeter level land surface deformation detection based on DInSAR technique, and calculated the land surface deformation of the watering area in Inpeirier Valley in Southeast California (USA) using SeaSat L-band SAR images (Gabriel and Goldstein 1988).Massonnet et al. (1993) calculated the deformation field of Lander Earthquake occurring in 1992 using ERS 2/1-SAR data processing.The comparison showed that the results of DInSAR were consistent with measurements from other techniques and the elastic deformation model (Massonnet et al. 1993).Their study initialized the land surface deformation detection with DInSAR technique all over the world.In 2002, Wang et al. for the first time in China, extracted the deformation field of Zhangbei-Shangyi Earthquake (6.2 M) occurring in 1998 using ERS 2/1 radar data by DInSAR technique (Wang et al. 2001).With the theoretic observation precision of centimeter and even to millimeter level, InSAR could be widely used in the earthquake deformation field, volcano deformation, glacier movement and city subsidence monitoring (Hanssen 2001;Wang 2007) In order to detect the source of the sandstorm and its influence, EnviSat ASAR satellite radar data are processed using two-pass differential interferometric measurement techniques to monitor desert surface deformation in Hunshandake Sandy Land in the paper.

THE PRINCIPLE OF INSAR DEFORMATION MONITORING
As shown in Fig. 1, { is the phase difference of the reflection signal from the same land surface which is received by two antennae.{ is caused from the distance difference δ between two antennae (Hanssen 2001).In two-pass inter-ferometric measurements, the distance difference of signal transmission is 2δ, and so According to the geometric relation shown in Fig. 1, we can find where, Z(y) is the topographical height, h is the satellite altitude, λ is the wavelength of radar wave, { is the phase difference, B is the baseline length, θ is the angle of view, and α is the angle of inclination.According to Eq. ( 2), we can calculate the topographical height as long as we can calculate the related parameters.Meanwhile, the azimuth pixel position is correlated with the imaging epoch, and the distance pixel position is correlated with the distance from the antenna to the land surface (Rantakokko and Rosenholm 1999).

DINSAR MONITORING OF HEIGHT CHANGE IN DESERT AREA
The Hunshandake Sandy Land area is located at longitude 112 to 118°E and latitude 42 to 44°N as shown in Fig. 2. Most of the Hunshandake Sandy Land is covered by sand with only a small area vegetated.The goal of this study is to detect the sand dune movement and monitor changes of the sand surface.
EnviSat (ENVIronmental SATellite) is equipped with ASAR, which can acquire high-resolution SAR images of the globe under any weather conditions with an active phase  Table 1.Basic parameters of the interferometric image pair.
array at incidence angles of 15 -45° (http://envisat.esa.int/).In this study, considering their temporal-spatial baseline and interferometric quality, two ASAR images of January 20 and October 26, 2004 were selected for preprocessing and analysis of their interferometric capability, shown in Table 1.Then we can obtain the deformation between January 20 and October 26, 2004.The center of the research area was located at 114.5°E and 42.7°N.
The two-pass algorithm was employed for raw data interferometric processing using the Doris software (Delft object-oriented radar intererometric software).Doris, as one free software with open source codes, is designed to obtain three-dimensional information and deformation of land surfaces (Kampes and Usai 1999).SRTM (Shuttle Radar Topography Mission) DEM issued by NASA (National Aeronautics and Space Administration, USA) was used to delete topography effects.SRTM DEM, covering about 80 percent of the surface of the Earth from 60°N to 54°S, is built with WGS-84 reference ellipsoid as the horizontal datum and EGM96 geoid as the vertical datum, according to  (Farr et al. 2007) The deformation error in the sight direction due to the DEM model error is not larger than 5 mm (Massonnet and Feigl 1998).
In our data processing, the orbital ephemeris is from the precise orbit of DORIS (Doppler Orbitography and Radiopositioning Integrated by Satellite) positioning system released by ESA (European Space Agency).The filtering was implemented with the Goldstein method.The Snaphu program (Chen and Zebker 2000) developed by Stanford University was employed in the phase unwrapping.Snaphu computes the initial position of each point according to the orbit parameters based on a statistical method, and then, considering the interferometric capability between points, calculates the best solution from the point with the highest interferometry to the surrounding points according to maximum a posteriori.In this way, the phase unwrapping of the full interferometric map can be obtained.A flowchart of the data processing is shown in Fig. 3.
First a coarse co-registration was performed using the precise orbit from DEOS (Delft Institute for Earth-Oriented Space Research, Netherlands) to find the initial displacement (38, -1) of the image center.The precise co-registration is a key problem to obtain the interferogram.Theoretically, 1/10 pixel accuracy is essential.Once the control points are determined, the relative rectification (pixel re-sampling) is also important.The coherence coefficient was used due to the initial displacement of 38 for azimuth direction and -1 for distance direction.Then a sub-pixel fine co-registration was performed using the coherence coefficient to obtain the geometric displacement fitting parameters of the image by polynomial fitting, shown in Table 2.After re-sampling the slave image, interference processing, and flat earth effect removal processing, we obtain the interferogram including the topography phase and displacement phase shown in Fig. 4, which is from 10 : 2 (azimuth : range) multi-look processing.
We computed the topographic phase correction shown in Fig. 5 according to SRTM DEM, and performed differential interference processing to arrive at the differential interferogram shown in Fig. 6.On the whole, the fringes with high coherence in the interferogram are clear, continuous, and of high resolution, indicating that the interferometirc capability of these two ASAR images was acceptable under the conditions of dry weather, low vegetative cover, and good temporal baseline.The ranging error of InSAR is mainly from the atmospheric delay, DEM, orbit error, etc., among which atmospheric wet delay errors are the largest.In the study area, the water vapor content was low.According to Wright et al. (Wright et al. 2004), errors from changes in water vapor content are smaller than 1 mm a -1 .In addition, the precision of DORIS initial orbit is better than 0.2 m, and the radial precision of the post orbit determination can approach 5cm (Scharroo and Visser 1998).Deformation errors in the sight line direction from DEM with a horizontal precision of 20 m and a vertical precision of 15 m are smaller than 5 mm (Massonnet and Feigl 1998).The rang-  ing error from DEM and orbit error will be lower because a perpendicular baseline was selected.After the phase unwrapping and geo-coding, we obtain the deformation field of the desert surface in the sandstorm source, as shown in Fig. 7, which is a three-dimensional deformation map of the red rectangle shown in Fig. 6.
The three-dimensional desert surface deformation field is shown in Fig. 7, in which X is in latitude direction, and Y is in longitude direction, and the colors show different values of deformation (unit in centimeter).From Fig. 7, we can find that the desert surface of the research area was obviously deformed during the 280 days.Area A markedly decreased in height, and the ascending surface occurred mainly in area B, where the largest rise and subsidence were greater than 1 m.In order to analyze the deformation in detail, we constructed the vertical profiles in latitude and longitude directions of the area showing the largest increase in height in area A and the largest subsidence in area B, shown in Fig. 8.In the profiles, deformation values and trends in the two directions were obvious.For further visual analysis, we computed deformation isolines for area A and B as shown in Fig. 9, in which the main isolines are -160, -145, -130, -115 and 15, 30, 45, 60, 75, etc. (cm).
In order to compare the monitoring results with results from the DInSAR technique, we made two GPS RTK surveys of the desert land surface in the study area within two years.Positioning precision with the relative GPS positioning technique can be up to the mm level (Guo et al. 2005).In each survey, we set two GPS reference stations and made continuous 24-hour observations to obtain precise coordinates.At the same time, observations were made on the moving stations.Considering the topography of the observation area, the sample time interval was 1 second.According to the results of data post-processing using Bernese software (http://www.bernese.unibe.ch/), the precision in height was as good as 5 cm.We built the desert surface height model using curve fitting.Comparing the height difference from the model with the results from the DInSAR, the biggest difference was 0.6 m, and the average difference was 0.13 m.
The coherence of complex image pairs was not good enough to some extent due to long time interval between the image pairs.Besides temporal and spatial factors, the interference quality is related with the scattering characteristics of the image area.In the study area, the interferometric capability should be excellent because of very low vegetative cover, but atmospheric effects, temperature difference, and other factors will affect the coherence.More details on ways  to improve the image coherence over desert areas will require further research.

ANALYSIS AND DISCUSSION
Sandstorms were frequent over the northern China in the spring of 2004.In Beijing, there were 8 sandstorms, and they became stronger and stronger, especially the sandstorm occurring on March 26 -29, bringing Beijing and surrounding 11 provinces the worst sandstorm in the last several years.During that sandstorm, the city of Beijing was covered by sand clouds with winds of level 5 -7 and wind gusts to level 8.The visibility in part of the city was only several meters.The air was seriously polluted by the sandstorm.The Hunshandake Desert was the main source of the sandstorm, which was related to a Mongolian cyclone (Zhang 2004).In 2004, sandstorms in Inner Mongolia came early, frequently and strongly due to the basic conditions for sandstorms: abundant sand sources, strong wind, and unstable thermal structure of the atmosphere.This makes Hunshandake and Keerqin Sandy Land the important route of sand and dust transportation (Zhi 2004).
As shown in Fig. 6, from January 20 to October 26, 2004, the height of the test area changed greatly, ranging from -160 to +120 cm.The height in part of the test area increased (as shown in area B and C of Fig. 6), but in most of the area, the height decreased.On one hand, some sand, which is in a small amount, from Mongolia and west Inner Mongolia was brought to Hunshandake by the wind.On the other hand, as a sand and dust source, Hunshandake offers a large amount of sand for dusty weather and sandstorms, which results in the decrease of the desert land surface in most of the study area.The largest height change was -160 cm because of the very serious sandstorm in 2004.
By using longer wavelength microwaves, which can penetrate dry desert and soil with low permittivity, InSAR technology was used to obtain more detail about the land surface and underground targets than other earth observation tools, including visible light and near-infrared, could provide.In desert areas, accurate co-registration was easy to realize due to the sparse vegetation.But it was difficult to obtain stable scattering data due to the movement of dunes, and the interferogram, including the plane displacement and the height change, will show some distortions.

CONCLUSIONS
In order to detect the sandstorm source and monitor changes in the height of the land surface, InSAR techniques were employed for the first time to monitor the changes in Hunshandake Sandy Land.Two EnviSat ASAR images for January 26 and October 20, 2004 were selected for DInSAR processing because of the strong and frequent sandstorms over northern China in the spring of 2004.A three-dimensional map and isoline map of land surface deformation of the study area were obtained.According to the results, the height change of the study area ranged from -160 to +120 cm.In most of the study area, the height decreased, indicating the sandstorm source.In some parts of the study area, the desert land surface height increased due to sand brought by the wind from Mongolia and other places, and the movement of the sand dune.Our results indicate that DInSAR technique is useful for monitoring height changes in desert areas, and has the potential to become an important tool in studies on sandstorm erosion and deposition.
This research was supported in part by the National Natural Science Foundation of China under grant 40974016 and 40774009, the National Hi-tech R&D Program of China under grant 2009AA121402, the Special Project Fund of Taishan Scholars of Shandong Province, China under grant TSXZ0502, and the Research & Innovation Team Support Program of SDUST, China.
SAR data from February 11 to 22, 2000 by a US space shuttle (http://srtm.csi.cgiar.org/).The spatial resolution of this model is 3m × 3m, and the horizontal and vertical precisions are 20 and 16 m respectively

Table 2 .
Coregistration fitting parameters of the interference image pair.