Geostatistical Data Fusion of Multiple Type Observations to Improve Land Subsidence Monitoring Resolution in the Choushui River Fluvial Plain , Taiwan

Precise leveling, global positioning system (GPS), and interferometric synthetic aperture radar (InSAR) have had marked influences on the geodesy field. Many studies have applied these techniques to land subsidence monitoring and predictions. Land subsidence in the Choushui River Fluvial Plain (CRFP) in western central Taiwan is severe, primarily because of excessive groundwater pumping. Exploiting various observational techniques this study employed a geostatistical cokriging algorithm to integrate multiple types of observations for mapping land subsidence in the CRFP between 1993 and 2008 to effectively reduce the regional effects and interpolation biases in InSAR observations. Precise leveling data and persistent scatterer InSAR results were first assessed through variogram analysis, the results of which revealed similar directional variograms and no nugget effect for the experimental variograms. All observed sill values decreased with increasing temporal intervals and the ranges of the same type of observations tended to be stable for time intervals of different lengths. The accuracy of a cokriged land subsidence map was verified using continuous GPS data at the same cell locations and was significantly improved compared to when the precise leveling data were used. The cokriging interpolations of land subsidence in the CRFP indicated that the severe subsidence areas had moved from the coastal area to the central CRFP. An obvious critical subsidence migration point for the southern CRFP occurred in 1998. Overall subsidence rate trends increased in the early periods (1993 1998 and 1998 2003) and decreased in the most recent period (2005 2008).


InTRODuCTIOn
Managing groundwater resources is imperative for areas with limited surface water supplies or rapid urban development.Long-term and unrestricted groundwater withdrawal can lead to serious land subsidence.Severe land subsidence has been recognized in many cities such as Shanghai, Las Vegas and Bangkok (Hoffmann et al. 2001;Phien-wej et al. 2006;Liu and Helm 2008).The consequences of land subsidence can be significant because non-uniform compaction and sinking lowlands might damage public constructions and increase the risk of flooding (Abidin et al. 2001).Numerous geodetic techniques with different scale resolutions are usually employed to monitor and quantify the behavior, degree and deformation speed of land subsidence in a specified area.Over the years, precise leveling, global positioning system (GPS) and interferometric synthetic aperture radar (InSAR) have proven to be powerful techniques in the field of geodesy.Many studies have applied these techniques for land subsidence monitoring and prediction (Reilinger et al. 1977;Dixon 199l;Yu and Chen 1994;Galloway et al. 1998;Massonnet and Feigl 1998;Ferretti et al. 2001;Sato et al. 2003;Chang et al. 2004;Hooper et al. 2004Hooper et al. , 2007;;Galloway and Burbey 2011;Yen et al. 2011;Wu et al. 2013).
The Choushui River Fluvial Plain (CRFP) in central Taiwan has been affected by severe land subsidence since the 1970s.The extent of accumulated land subsidence in the CRFP varies from a few centimeters to several meters, depending on the locations of aquifer materials and groundwater withdrawal patterns.The Water Resources Agency (WRA) and Central Geological Survey (CGS) have conducted long-term hydrologic and geologic studies for aquifer characterization and land subsidence monitoring (CGS 1999;Liu et al. 2001).The collection of geodetic data, including precise leveling, GPS, and recently installed compaction monitoring well (CMW), can provide detailed information on the temporal and spatial variations in surface deformation in the CRFP (Liu et al. 2004;Hwang et al. 2008).However, using precise leveling to monitor a large area regularly is time-consuming and expensive.Installing and maintaining even a limited number of continuous GPS and CMW stations to resolve the spatial resolution of local land subsidence is costly.Differential interferometric synthetic aperture radar (DInSAR) is an alternative for mapping the spatial distribution of land subsidence at a relatively high resolution (Galloway et al. 1998;Hoffmann et al. 2001Hoffmann et al. , 2003;;Calderhead et al. 2011).This approach has been applied in quantifying land subsidence in many areas of Taiwan (e.g., Chang et al. 2004;Hung et al. 2010).
Previous investigations have recognized that the errors in DInSAR observations can be temporal and spatial decorrelations such as atmospheric delay, baseline errors and digital elevation model (DEM) errors.These regional effects might limit the extraction of meaningful information for identifying surface deformations (Hooper et al. 2004).Ferretti et al. (2000Ferretti et al. ( , 2001) ) proposed the identification of a Permanent Scatterer (PS) for any SAR image to improve the accuracy of DInSAR in temporal and spatial decorrelations.The PS is considered stable and continuous in temporal series and can return strong signals to dominate a pixel for improving spatial resolution.Hooper et al. (2004Hooper et al. ( , 2007) ) and Hooper (2008) successfully obtained deformation signals from agricultural and vegetational SAR data, improving the accuracy of estimated displacements in areas with nonlinear deformation behavior.Tung and Hu (2012) modified the PS approach and employed a Delaunay triangulation algorithm developed by Liu et al. (2009) to extract the PS pixels in non-urban areas in the southern CRFP.The average pixel density was increased by nearly three PS pixels km -2 in their study.
Recently developed computational technologies can integrate multiple types of observations to improve land subsidence monitoring resolution.This is particularly important for areas where the observations from precise leveling, GPS and InSAR images are sparse and limited in spatial and temporal scale.Hung et al. (2011) integrated precise leveling and PS observations and identified PS pixels and extracted more than 100 pixels km -2 in the CRFP from 2006 -2008.From the differences between persistent scat-terer InSAR (PSI) and precise leveling data at the leveling benchmarks, a kriging algorithm revealed the distribution of the differences between the PSI and precise leveling data in the CRFP.The distribution of the differences was then added back into the original PSI observations to improve the land subsidence monitoring accuracy.The accuracy of such data integration approaches requires sufficient leveling benchmarks and PSI because the kriging weightings depend heavily on the leveling benchmark number and distribution.In the CRFP, however, the leveling benchmarks varied considerably before and after 2000 (see Fig. 1).Before 2000 the subsidence areas were located mainly along the shoreline and the leveling benchmarks were installed near the subsidence areas.Numerous leveling benchmarks were installed in the CRFP after 2000 and the inconsistency of the leveling benchmarks between the two periods (before and after 2000) might lead to logical and technical difficulties when implementing the approach proposed in Hung et al. (2011).
The present study employed a geostatistical cokriging interpolation algorithm to integrate multiple types of land subsidence observations in the CRFP between 1993 and 2008.The precise leveling data and PSI results were first assessed through a variogram analysis to obtain representative geostatistical structures.A nonlinear least squares model was developed to obtain the best fit of the geostatistical structures for the precise leveling data and PSI results.From such geostatistical structures cokriging interpolation can be applied to map land subsidence in different periods in the CRFP.The cokriging interpolations were verified with continuous GPS data at the same locations.The CRFP surface deformations during different time periods were then systematically quantified with detailed comparisons.

HyDROGeOLOGICAL DeSCRIPTIOn AnD
SOuRCeS OF ObSeRvATIOnS

Hydrogeological Description
Taiwan is located at the boundary between the Philippine Sea Plate (PSP) and Eurasian Plate (EUP), where the currently ongoing collision started at approximately 5 -7 Ma.The PSP moves northwest toward the EUP at an average rate of 82 mm yr -1 (Yu et al. 1997).The tectonic compression divides Taiwan into four physiographic regions, which from west to east are the Coastal Plain, the Western Foothills, the Central Range, and the Coastal Range (Suppe 1984;Teng 1990).In central Taiwan the CRFP is part of the Coastal Plain and the CRFP area is approximately 2500 km 2 .The Changhua and Chiuchiunkeng Faults are defined as two segments of the geological boundary between the Coastal Plain and the Western Foothills.These two faults run along Pakua Tableland and Taolao Hill and are treated as the eastern CRFP boundary (see Fig. 1).The other CRFP boundaries are the Taiwan Strait to the west, the Wu River to the north and the Putzu River to the south.
The Choushui River is the longest river in Taiwan.Its flow areas cover the Central Range, the Western Foothills and the Coastal Plain.The CRFP sedimentary thickness varies from 750 -3000 m and the sediment grain sizes decrease from the east to the west (Lin et al. 1992).The sedimentary environment is frequently altered by factors such as rising sea levels, river flooding and river channel migration.The materials deposited on the riverbed include clay, fine sand, coarse sand, and gravel (CGS 1999).
The WRA developed a groundwater monitoring network in 1990.The CRFP has over 200 groundwater monitoring wells.The depths of these wells vary from 200 -300 m depending on the well location.The CRFP aquifer system can be classified into three to four aquifer layers and several aquitard layers based on well logs.The average annual rainfall in the CRFP is 1200 -2000 mm.However, over 80% of the annual rainfall is concentrated from May to September.Such rainfall patterns provide a large amount of groundwater that can be pumped out during dry seasons (from October to April of the next year).Most pumping wells pump groundwater from Aquifer 2 and 3 (CGS 1999;Liu et al. 2001Liu et al. , 2004;;Chen et al. 2010).The CRFP has long been characterized by high-density agriculture, aquaculture and industrial parks.Surface water storage capacity is limited.The groundwater resources constitute the main water re-source in the CRFP.Excessive groundwater pumping has caused land subsidence since the 1970s.

Precise Leveling
Since 1975 land subsidence surveys have recorded precise leveling data for the CRFP.The WRA collects such data based on the first-or second-order standard of leveling procedures.The leveling procedure specifications require that any loop misclosure should be less than 3 mm K , where K is the distance between two neighboring benchmarks in kilometers.Precise leveling is best accomplished during the wet season or rainy season because of minimal groundwater withdrawal and soil compaction in those periods (Hung et al. 2010;Galloway and Burbey 2011).
The main advantage of precise leveling is its high accuracy (the reliability value is 2 mm for two leveling benchmarks within 1 km), but monitoring a wide area regularly is extremely time-consuming.Figure 2 shows the precise leveling benchmark count between 1990 and 2010.Before 1998 the benchmark numbers were less than 100.After 2000 the benchmark numbers increased dramatically.Benchmarks cover a large area in the CRFP and resolve the areas where land subsidence is severe.According to WRA records the CRFP precise leveling surveys run from 1992 -2008.We divided these 17 years into three periods: Period 1 (1992 -1998) contains 34 benchmarks, Period 2 (1998 -2003) contains 99 benchmarks, andPeriod 3 (2005 -2008) contains 298 benchmarks.We further defined three overlapping intervals for Period 3: Interval I (2005 -2006), Interval II (2005-2007), and Interval III (2005-2008).Between 1992 and 2008 the mean NB061 velocity (precise leveling benchmark) was -1.27 mm yr -1 .This leveling benchmark is the most stable one of all precise leveling benchmarks during these years.Therefore, the location of NB061 was considered the reference point for this study.

Persistent Scatterers InSAR (PSI)
The CRFP area covers two SAR image modes (Track 232,Frame 3123;and Track 232,Frame 3141) (Hooper et al. 2007;Hooper 2008).The topographic effect was removed using the DEM from NASA's SRTM mission and the precise orbit parameters from the Delft Institute for Earth-Oriented Space Research were used to correct the orbit errors.  2 shows the spatiotemporal distribution of all interferograms that were used in this study.

MeTHODOLOGy
The spatial resolution of PSI observations was 500 times, which was higher than that of CRFP leveling observations (Hung et al. 2011).To simplify the analysis of the leveling and PSI data we generated 9867 active cells for the CRFP (see Fig. 1).The uniform cell size was 500 × 500 m.Different cell sizes and boundaries might lead to slight differences in the results.The geostatistical parameters required for the cokriging analysis were estimated through a variogram analysis.The FORTRAN routines gamv and cokb3d in GSLIB (Deutsch and Journel 1998) were employed to develop experimental variograms and perform cokriging interpolation, respectively.The estimated land subsidence was then compared with the available continuous GPS data from several GPS stations in the CRFP.In the following sections we briefly describe the algorithms used for the data analysis.

experimental variogram
The experimental variogram analysis is the key procedure in obtaining geostatistical structures for a site-specific condition (Deutsch and Journel 1998).A variogram describes the spatial dependency of referenced observations in a one-or multidimensional space.An experimental semivariogram, ( ) h cl , is defined as half of the average squared difference between two univariate values, approximately separated by vector h, and can be expressed as: where N(h) is the number of pairs, and z(u i ) and z(u i + h) are the respective values of precise leveling or PSI at locations u i and u i + h of pair i. Vector h can be specified with particular directions and lag distance (Deutsch and Journel 1998).

The cross semivariogram ( )
is defined as half of the h-increment average product relative to two variables (bivariate) and can be computed as follows: where z L (u i ) and z L (u i +h) are the respective precise leveling values at locations u i and u i + h of pair i, and z P (u i ) and  (cross semivariogram) versus the incremental distance h.The experimental variograms must be described with continuous functions for most applications such as Gaussian, exponential, or spherical models.Based on the leveling and PSI data analysis in the CRFP the Gaussian model variogram was used for the continuous function in these analyses.The Gaussian model has the following mathematical form (Deutsch and Journel 1998): where C is the sill and a is the range.Equation ( 3) can be applied for multidimensional problems.The C and a values were obtained according to the best fit between the experimental variogram and Eq. ( 3).However, the parameter adjustments are highly dependent on both the specific site conditions and user experience.
A nonlinear least squares algorithm was employed in this study to develop a fitting program for estimating parameters C and a.For the observation data spatial anisotropy analysis Fig. 3 illustrates how vector h was applied to search the observation pairs with the defined parameters in the gamv routine in GSLIB.We used 1 km for the lag distance and 0.5 km for the lag tolerance, and considered six directions when analyzing the variograms.The azimuth angles were North, N30°E, N60°E, N90°E, N60°W, and N30°W.The tolerance angles for different directions were fixed to 15° and the bandwidth for each angle was fixed to 5 km.

Cokriging Interpolation
The cokriging interpolation is a linear combination of primary (i.e., the precise leveling) and secondary (i.e., the PSI) data.It minimizes the variance in the estimation error by exploiting the cross-correlation between two variables  (Isaaks and Srivastava 1990).Through various observational techniques when a subsidence area involves multiple types of observations (such as precise leveling and PSI) for the same period.The observations are usually considered cross-correlated and can be integrated to improve the spatial resolution and land subsidence monitoring accuracy.
The precise leveling observations in this study were defined as the primary variables and are denoted as ( ) L u 1 a , whereas the PSI observations were defined as the secondary variables, denoted as ( ) u P 2 a .The cokriging estimations of Z COK (u) for a desired variable u can be expressed as: where 1 m a is the weighting applied to the 1 a th precise leveling observation, and u 1 a is the precise leveling benchmark at the location 1 a ; and 2 m a l is the weighting applied to PSI observation 2 a , and u 2 a l is the PSI pixel at location 2 a .To reduce the influences on the weightings from using different observation types we employed the standardized ordinary cokriging, which involves creating new secondary variables with the same means as the primary variables.Subsequently, all of the weightings are constrained and the weightings summation is equal to 1. From this modification Eq. ( 4) yields with the condition where m L and m P are the stationary means of precise leveling and PSI observations, respectively.

ReSuLTS AnD DISCuSSIOn
This section first presents the variogram analysis of different types of observations and the cokriging interpolations for Period 3 (sections 4.1 and 4.2).The GPS data are employed to verify the accuracy of the integrated observations (section 4.3).Subsequently, the cokriging interpolations for Periods 1 and 2 are introduced using the same procedure (section 4.4), and the cokriged land subsidence results are compared with those obtained from other studies (section 4.5).

Different Types of Observations and variogram Analyses
Fifteen Envisat interferograms from Track 232, Frame 3123 and 17 Envisat interferograms from Track 232, Frame 3141 were used to estimate the radar line-of-sight (LOS) displacement for Period 3. In the same period 298 precise leveling data were recorded for CREP vertical deformation.The radar LOS was 23° from the vertical deformation and the outcome was an accumulation of three-dimensional movement.According to Ching et al. (2011) the greatest horizontal GPS velocity in the CRFP between 1995 and 2005 was close to 5 mm yr -1 .This was calculated by comparing against the reference point.Assuming that all horizontal vectors and the azimuth (the flight direction of satellite) are perpendicular, the 5 mm yr -1 of horizontal velocity contributed only 1.95 mm yr -1 in LOS velocity.Because of the negligible effect of the horizontal velocity in the LOS velocity, the PSI results were relatively sensitive to vertical motion.The first row in Fig. 4 shows the vertical displacement rates of the precise leveling data and PSI results in the CRFP from Intervals I to III.The cold colors represent the land subsidence relative to the reference point and the warm colors show the surface uplift.Both the precise leveling data and PSI results clearly indicated two local subsidence areas by the north and south sides of the Choushui River.The highest subsidence rates derived from the precise leveling data from Intervals I to III were 84, 78, and 76 mm yr -1 , whereas those derived from the PSI results for the corresponding periods were 94, 73, and 64 mm yr -1 .The shapes of the subsidence areas, locations of the subsidence centers and patterns of the subsidence rates between these two types of observations were similar in the CRFP area.However, a plot of all observations on the same scale would not provide a clear quantification of their similarities and differences.To quantify the correlations between these observations, semivariograms, and cross semiveriograms were employed to analyze the spatial correlations.
The other three rows in Fig. 4 show the leveling variogram, PSI variogram, and cross variograms in six directions.The red dotted line and triangles represent the north direction, the green dotted line and squares represent N30°E, the blue dotted line and stars represent N60°E, the yellow dotted line and diamonds represent N90°E, the purple dotted line and cycles indicate N60°W, and the gray dotted line and hexagon represent N30°W.The black solid line is a summary of six directional semivariograms fitted with a Gaussian model and the associated parameters including the sills and ranges in the different intervals.Table 2 lists the detailed parameters for the spatial analyses.
In Interval I the fitted sill and range of the precise leveling data were 741 and 20.395 km.The optimally fitted sill and range of the PSI results were 1074 and 14.154 km; and the fitted sill and range of the PSI and leveling cross variogram were 393 and 12.703 km, respectively.In Interval II the fitted sill and range of the precise leveling data were 654 and 21.837 km.The fitted sill and range of the PSI results were 636 and 14.003 km; and the fitted sill and range of the PSI and leveling cross variogram were 314 and 11.993 km, respectively.In Interval III the fitted sill and range of the precise leveling data were 526 and 20.172 km.The fitted sill and range of the PSI results were 493 and 14.039 km; and the fitted sill and range of the PSI and leveling cross variogram were 278 and 11.964 km, respectively.No nugget effect was introduced into the experimental variograms in this study.
The sill values for the different types of observations generally decreased from Intervals I to III because of the retarded subsidence rate.The seimivariogram plots showed that the differences among the sills in different directions were nonsignificant.The PSI variograms appeared to be relatively smooth compared with those of the precise leveling data because the PSI results were averaged beforehand in every active cell.The precise leveling data sill, however, was smaller than the PSI results for the different time intervals, indicating that the data variations among the leveling observations were relatively small.For the same observation types the ranges in the experimental variograms were nearly constant among the different time intervals.The range of variograms reflected the spatial correlation and subsidence area patterns in the CRFP.In the context of the correlation between the precise leveling data and the PSI results, the cross variograms displayed similar sill temporal patterns and similar spatial patterns in the ranges.The sill patterns reflected the variations in subsidence slowing gradually from Intervals I to III.However, the range patterns approximately indicated the radii of the subsidence areas, which were similar throughout the three time intervals.

Cokriging Interpolations for the Three Intervals in Period 3
We divided Period 3 (2005 -2008) into three intervals: Interval I (2005Interval I ( -2006)), Interval II (2005-2007), and Interval III (2005-2008), to produce contrasting analyses from the different time-averaged lengths.We employed the fitted parameters in Table 2 for cokriging interpolations to obtain the land subsidence in the spatiotemporal variations.Note that we used a low nugget effect value of 1 to maintain the numerical stability in the Gaussian variogram model (Deutsch and Journel 1998).
The cokriging interpolations in the three intervals are shown in the upper row of Fig. 5.The negative rate values (shown in cold colors) represent the areas with land subsidence.The largest land subsidence rates in Intervals I, II, and III were -99.98, -93.17, and 80.97 mm yr -1 , respectively.In Interval I the land subsidence rates exceeding -50 mm yr -1 covered 1760 cells with an area of 440 km 2 in the CRFP.The same subsidence rates for Interval II covered 1360 cells with an area of 340 km 2 and those for Interval III covered only 1266 cells with an area of 316.5 km 2 .The average subsidence rate or area decreased over time.The bottom row of Fig. 5 summarizes the statistics regarding the differences between the cokriging interpolation and PSI results for these three intervals.The distribution of the differences was close to a normal distribution with the respective skewness values of -0.302, -0.238, and -0.277.The root mean square (RMS) values in these three intervals were 28.88, 21.89, and 23.71 mm yr -1 .

Cokriged Results Comparisons with Continuous GPS Observations
The five continuous GPS stations in the CRFP were acquired from WRA.Two stations, Hsi-Kang (HK) and Hsin-Hsing (HH), recorded the vertical displacement from 2005 -2008 and the other stations, Tu-Ku (TK), Lin-Nei (LN), and Ko-Tso (KT), recorded the vertical displacement from 2006 -2008.Period 3 in this study was artificially divided into three intervals.Therefore, the mean velocity of the GPS data from 2006 -2007 was used for comparison with the cokriging interpolation for Interval II, and the GPS data from 2006 -2008 were used for a comparison with the cokriging interpolation for Interval III.Table 3 lists detailed information for the GPS data.Figure 6 shows the direct comparisons of the GPS data with the cokriking interpolations, precise leveling data and PSI results.The correlation coefficient for the cokriging interpolation and GPS data at the same cell is 0.9603, and the RMS is 10.56 mm yr -1 (marked with a red circle in Fig. 6).The correlation coefficient and RMS for the precise leveling data and GPS data are 0.9495 and 4.95 mm yr -1 , respectively (marked with a green circle in Fig. 6).For the PSI results and GPS data, the correlation coefficient and RMS are 0.3256 and 31.49mm yr -1 , respectively (marked with a blue circle in Fig. 6). Figure 6 clearly shows that the precise leveling data gives the highest correlation with the GPS data (close to   -72.463 (2005 -2006) a 45° dotted line).The PSI results are less accurate because large variations were obtained for the specified GPS stations.Integrating the precise leveling data and PSI results improved land subsidence monitoring for the specified points.However, the estimation method consistently underestimated the land subsidence.
Notable advantages to integrating the precise leveling data and PSI results based on cokriging interpolation are described as follows.The precise leveling technique has been recognized as an accurate approach for monitoring land subsidence.However, the precise leveling approach demands extensive resources over long periods, which limits the availability of precise leveling data for large areas.The PSI technique provides high-resolution spatial observations through an identified PS for a large area.This study employed cokriging interpolation to exploit the precise leveling and PSI techniques to obtain high-resolution data with greater accuracy for monitoring land subsidence.

Cokriging Interpolations for Periods 1 and 2
We verified the cokriging interpolation results with GPS data for Period 3. In this section we follow same analytical process and estimate the cokriging interpolations for Periods 1 and 2. Period 1 covers 34 precise leveling data, 26 ERS interferograms from Track 232, Frame 3123 and 26 ERS interferogram from Track 232, Frame 3141. Figure 7a shows the precise leveling data and PSI results.Near the shoreline some precise leveling data is shaded black to indicate that the subsidence rate was higher than 100 mm yr -1 .These subsidence centers, however, do not appear in the PSI results.The variogram analysis results are displayed in  7b.Subsidence rates greater than -100 mm yr -1 are marked in black; they are distributed over the map coastal area.One hundred ninety cells with subsidence rates in excess of -100 mm yr -1 cover a total area of 47.5 km 2 .The total CRFP area with subsidence rates between -50 and -100 mm yr -1 is 451.25 km 2 (1805 cells).Figure 7d shows the statistics for the differences between the cokriging interpolation and PSI results.The skewness is -1.757 and RMS is 30.89mm yr -1 .
For Period 2, 28 ERS interferograms from Track 232, Frame 3123, 23 ERS interferogram form Track 232, Frame 3141, and 99 precise leveling data were used for the cokriging interpolation.Figure 8a shows the PSI results and precise leveling data.Some precise leveling data, shaded in black, show subsidence rates in excess of 100 mm yr -1 near the northern side of the Choushui River mouth.Although the subsidence centers indicated by the precise leveling data do not appear clearly in the PSI results, the shapes and patterns of the subsidence areas in these two types of observations are similar.The variogram analyses are shown in Fig. 8c and the fitted parameters are listed in Table 4.The sill and range values of the precise leveling variogram are 884 and 16.151 km, those of the PSI variogram are 168 and 30.892 km, and those of the PSI and leveling cross variogram are 236 and 25.534 km, respectively.Figure 8b shows the cokriging interpolation results for Period 2. Two hundred nineteen cells with subsidence rates in excess of -100 mm yr -1 are shaded in black, covering a total area of 54.75 km 2 .In the northern part of the CRFP the subsidence center is located in the coastal area.However, in the southern part, the subsidence center is located in the central CRFP.The total area with subsidence rates between -50 and -100 mm yr -1 covers 561.5 km 2 (2246 cells).Figure 8d presents the statistics regarding the differences between the cokriging interpolations and PSI results.For Period 2 the skewness is -1.414 and RMS is 30.30mm yr -1 .These results are similar to those for Period 1.
The results for these two periods clearly indicate that severe land subsidence moved from the coastal area to the central CRFP.Additionally, regardless of whether the subsidence rate exceeded -100 mm yr -1 or was between -50 and -100 mm yr -1 , the total subsidence area increased from Periods 1 to 2. The negative skewness values (-1.757 and -1.414) for these two periods indicate that the cokriging interpolations generally estimated the land subsidence rate to be greater than the estimates derived from the PSI results.According to the cokriging interpolation procedures the precise leveling data were used for the primary variables and the interpolation integrates the leveling and PSI observations.We

Results Comparison with Observations from Previous Studies
The cokriging interpolation presents two subsidence centers located in the northern and southern CRFP.The subsidence migrations for both centers were from the coastal area to the central CRFP.However, their migration times appear to differ.The primary subsidence migration time span for the northern CRFP was between Periods 2 and 3.For the southern CRFP the subsidence migration occurred between Periods 1 and 2. Liu et al. (2004) reported on four precise leveling benchmarks (SC1, SC2, SC3, and SC4) installed along the coastal area that recorded vertical deformations from 1976 to 2001.The SC1 leveling benchmark is the northernmost of these four leveling benchmarks.The primary subsidence at SC1 occurred between 1976 and 1986.The SC2 leveling benchmark is located proximal to the subsidence center in the northern CRFP in Period 1.The land subsidence at SC2 became notable after 1991 and continued through 2001.During that 10-yr period the subsidence rate reached 150 mm yr -1 .The SC3 and SC4 benchmarks are located in the southern CRFP.Their subsidence rates decreased gradually after 1996.The historical data of these four leveling surveys confirm the cokriged subsidence behavior and verify the subsidence migration patterns in Periods 1 and 2.
Tung and Hu (2012) used their PSInSAR algorithm to process ERS images taken between 1996 to 1999.The results showed two subsidence centers in the southern coastal CFRP area and the southern central CRFP.The estimated subsidence rate reached 80 mm yr -1 .Although their study period (1996 -1999) overlaps with Periods 1 and 2 in the present study, the PSInSAR observations captured the subsidence characteristics for these two periods.The year 1998, between Periods 1 and 2, was a turning point during which the subsidence areas migrated to the southern CRFP.
For Period 3 our cokriging interpolation yielded results matching the data fusion of Hung et al. (2011), who combined leveling and PSI data for the 2006 -2008 period.In contrast to the study in Hung et al. (2011), we quantified the variations in geostatistical structure for different periods from 1993 -2008.The spatial correlations among the different types of observations were implemented in the cokriging interpolations for land subsidence mapping.Each cokriging interpolation was interpolated once with the spatial correlation between the precise leveling data and PSI results.Our procedure can reduce both the errors from regional effects and the biases of interpolated calculations from SAR images.Yen et al. (2008) adopted the DInSAR technique to process ERS images from between 1996 and 2000 and concluded that the Meishan Fault has a 10 mm yr -1 uplift rate.They contested the accuracy with which earthquake potential could be evaluated and the ability to identify surface deformations in the CRFP area.In Period 1 the southeast part of the CRFP, where the Meishan Fault and Chiuchiunkeng Fault are located, presented an uplifted deformation.Because there were no precise leveling data, the uplift rate was dominated by the PSI results and was calculated to reach 40 mm yr -1 .These findings should be incorporated with the land subsidence findings to clarify the deformation and earthquake probabilities in the CRFP area.

COnCLuSIOnS
Using a geostatistical method this study showed that cokriging interpolation is a systematic and effective procedure for exploiting multiple types of observations.The results can capture the spatiotemporal variations of different types of observations and map land subsidence behaviors in the CRFP.Based on our analysis results we argue the following: (1) Precise leveling and PSI observations were subjected to geostatistical analysis in this study.The procedure reduced the errors from regional effects and interpolation calculations.The results clearly delineated the historical surface deformations in the CRFP between 1993 and 2008.
(2) The variogram results for different types of observations for the three intervals for Period 3 showed that the six directional variograms were similar and no nugget effects were observed.The sill values of all observations decreased with the increase in temporal averaging (i.e., from Intervals I to III).For observations of the same type, the ranges of the variograms were nearly constant for all three intervals because all of them had similar subsidence areas.The optimally fitted parameters for geostatistical structures in the CRFP were the basis for the cokriging interpolations.
(3) In Period 3 the cokriging interpolation showed that the subsidence rates and areas decreased from Intervals I to III.The correlation coefficient for the cokriging interpolation and GPS data at the same cells was 0.9603 and the RMS was 10.56 mm yr -1 .The accuracy of the cokriging interpolation was notably improved compared with that of the PSI results and approximated the precise leveling data.(4) The cokriging interpolations for the three periods indicated that the severe subsidence areas migrated from the coastal area to the central CRFP and the turning point of subsidence migration for the southern CRFP occurred in 1998.There is an increasing subsidence trend from Periods 1 to 2 and a decreasing trend from Periods 2 to 3.

Fig. 1 .
Fig. 1.Study area.The black circle encloses the leveling benchmarks established before 2000, and the white circle encloses those set up after 2000.The red triangle is a continuous GPS station.The empty square shows the cell size adopted in this study.The red lines indicate the active fault locations.(Color online only)

Fig. 2 .
Fig. 2. SAR image spatiotemporal distribution and a histogram of the leveling benchmark numbers.(Color online only)

Fig. 3 .
Fig. 3. Searching condition for vector h in a variogram (adapted from Deutsch and Journel 1998).

Fig. 4 .
Fig. 4. Vertical deformation and variogram analyses for the precise leveling data and PSI results obtained between 2005 and 2008.The first row shows the vertical displacement rate of the precise leveling data and PSI results in the CRFP.The cold colors represent the land subsidence relative to the reference point and the warm colors show the surface uplift.The other three rows display the leveling variogram, PSI variogram, and cross variogram in six directions, North, N30°E, N60°E, N90°E, N60°W, and N30°W.(Color online only)

Fig
Fig. 7c and the fitted parameters are listed in Table 4.The sill and range values of the precise leveling variogram are 1152 and 7.775 km, those of the PSI variogram are 145 and 18 km and those of the PSI and leveling cross variogram are 204 and 18.347 km, respectively.The cokriging interpolation results for Period 1 are shown in Fig.7b.Subsidence rates greater than -100 mm yr -1 are marked in black; they are distributed over the map coastal area.One hundred ninety cells with subsidence rates in excess of -100 mm yr -1 cover a total area of 47.5 km 2 .The total CRFP area with subsidence rates between -50 and -100 mm yr -1 is 451.25 km 2 (1805 cells).Figure7dshows the statistics for the differences between the cokriging interpolation and PSI results.The skewness is -1.757 and RMS is 30.89mm yr -1 .For Period 2, 28 ERS interferograms from Track 232, Frame 3123, 23 ERS interferogram form Track 232, Frame 3141, and 99 precise leveling data were used for the cokriging interpolation.Figure8ashows the PSI results and precise leveling data.Some precise leveling data, shaded in black, show subsidence rates in excess of 100 mm yr -1 near the northern side of the Choushui River mouth.Although the subsidence centers indicated by the precise leveling data do not appear clearly in the PSI results, the shapes and patterns of the subsidence areas in these two types of observations are similar.The variogram analyses are shown in Fig.8cand the fitted parameters are listed in Table4.The sill and range values of the precise leveling variogram are 884 and 16.151 km, those of the PSI variogram are 168 and 30.892 km, and those of the PSI and leveling cross variogram are 236 and 25.534 km, respectively.Figure8bshows the cokriging interpolation results for Period 2. Two hundred nineteen cells with subsidence rates in excess of

Fig. 6 .
Fig. 6.Comparison of the GPS data and other observations.Abbreviations of the GPS stations are as follows: HK (Hsi-Kang), HH (Hsin-Hsing), TK (Tu-Ku), LN (Lin-Nei), and KT (Ko-Tso).The numbers 1, 2, and 3 represent Intervals I, II, and III.The red circles are comparisons of the cokriging interpolations and GPS values.The green circles are comparisons of the precise leveling data and GPS values.The blue circles are comparisons of PSI results and GPS values.The 45° dotted line indicates that the subsidence rates of two observations are identical.(Color online only) Fig. 7. (a) Vertical displacement and variogram analyses for the precise leveling and PSI data obtained between 1993 and 1998.(b) Cokriging interpolation for Period 1. (c) The variogram analysis includes the leveling variogram, PSI variogram, and cross variogram in six directions.(d) The error statistics for the differences between the cokriging interpolation and PSI results.(Color online only) Fig. 8. (a) Vertical displacement and variogram analyses for the precise leveling and PSI data obtained between 1998 and 2003.(b) Cokriging interpolation for Period 2. (c) The variogram analysis includes the leveling variogram, PSI variogram, and cross variogram in six directions.(d) The error statistics for the differences between the cokriging interpolation and PSI results.(Color online only) Based on the time frames used in the precise leveling, Period 1 includes 26 pairs of ERS interferograms (master image on 27 March 1997) from Track 232, Frame 3123 and 26 pairs of ERS interferograms (master image on 27 March 1997) from Track 232, Frame 3141.Period 2 includes 28 pairs of ERS interferograms (master image on 3 September 1998) from Track 232, Frame 3123 and 23 pairs of ERS interferograms (master image on 3 September 1998) from Track 232, Frame 3141.Period 3 includes 15 pairs of Envisat interferograms (master image on 15 November 2007) from Track 232, Frame 3123 and 17 pairs of Envisat interferograms (master image on 15 March 2007) from Track 232, Frame 3141.Table 1 lists the detailed information for all SAR images.Figure

Table 1 .
SAR image information for this study.
Note: a: Length of perpendicular baseline component.b: Time difference relative to the master image.A negative sign means the time before the master image; a positive sign means the time after the master image.

Table 2 .
Theoretical variogram parameters fitted using a Gaussian variogram model in the three time intervals.

Table 4 .
Theoretical variogram parameters fitted using Gaussian variogram model in Periods 1 and 2.