Geomorphological Change Detection Using an Integrated Method: A Case Study on the Taan River Watershed, Taiwan

The Chi-Chi earthquake of 1999 in Central Taiwan caused significant vertical uplift and ignited a dynamic fluvial erosion system in the Taan River. In order to understand the geomorphologic changes in the Taan River after the Chi-Chi earthquake, this study uses multi-temporal aerial photos to produce orthoimages and surface models to derive regional three dimensional movements from 2001 2009. The Particle Image Velocimetry (PIV) technique is applied to identify the horizontal movements from both orthoimages and digital surface models (DSMs). While PIV can provide estimates for the channel change direction and magnitude; vertical analysis such as river incision rate derivation could be performed using the height difference between multi-temporal surface models. The results from these two approaches are then integrated for river geomorphological interpretation. It is observed that the proposed scheme is capable of identifying three river evolution stages from the extracted movements in the Taan case. These stages are: years 2001 to 2003, 2004 to 2006, and 2007 to 2009, respectively. The total incision between 2001 and 2009 is about 10 m. This study demonstrates that both PIV and DSM subtraction are effective in river geomorphological change identification. The integration of these two approaches could provide more information when observing the evolution of river morphology.


InTRoDUCTIon
River landform is shaped by erosion and deposition processes, and influenced by external forces, such as climate change, tectonics, and human factors.Natural disasters such as earthquakes and typhoons frequently dramatically change the erosion and deposition mechanisms of a river.In 1999, the Chi-Chi earthquake (M w = 7.6) struck Central Taiwan, causing many surface ruptures.The Taan River experienced co-seismic displacement due to this earthquake.The displacement was induced by reverse fault across the main river and changed the local fluvial system equilibrium.The fault scarps across the rivers had a vertical slip (up to 10-m) (Lee et al. 2003;Huang et al. 2013) and initiated the morphological channel evolution.
Modern remote sensing techniques such as digital photogrammtery, LiDAR (Light Detection and Ranging), and InSAR (Interferometric Synthetic Aperture Radar), are effective tools for topographic survey.Likewise, multi-temporal topographic data facilitated quantified river morphology evolution analysis (Hervás et al. 2003;Demirkesen 2008;Miller et al. 2009;Tseng et al. 2009;Huang et al. 2013).Among the different analysis schemes, Particle Image Velocimetry (PIV) (Adrian 1991), a pixel analysis technique, has been applied to aerial images for geomorphologic studies (Lo et al. 2009;Tseng et al. 2009).The PIV concept assesses the spatial gray scale similarities in image sequences and identifies correspondence.The changes can then be derived from these correspondences among multi-temporal images.Meanwhile, Digital elevation model (DEM) subtraction, that is, computing the difference between two DEMs of the same site but at different times, is a common approach for assessing topographic changes (Dewitte and Demoulin 2005;Dewitte et al. 2008;Corsini et al. 2009;Prokop and Panholzer 2009;DeLong et al. 2012;Daehne and Corsini 2013).This study proposes an integrated method combining PIV and digital surface model (DSM) subtraction.The objectives of this study are: (1) to determine the 3D movements using an integrated method; and (2) to compare the PIV river morphology results from different data sources (i.e., aerial images, photogrammetrically derived DSM, and DSM collected with airborne LiDAR).

MeThoDoloGy
As shown in Fig. 1 the proposed scheme includes applying PIV to both orthoimages and DSMs independently, DSM subtractions and the integrated analysis.The details for this scheme are presented in the following sequence: (1) PIV analysis, (2) DSM subtraction, (3) verification with insitu data, (4) integrated analysis, and (5) river morphology evolution interpretation.

PIV Analysis
PIV was originally extensively applied in hydromechanics for flow field observation (Adrian 1991).The fundamental PIV principle is to determine the corresponding points between two different images with image matching algorithms such as cross-correlation (CC) to calculate the similarities to establish quantitative correspondence.PIV could be applied to either orthoimages or DSMs.The DSM heights are normalized to 0 -255 before CC calculation.The PIV results are horizontal displacement vectors.This normalization procedure is conducted due to the limitations in using existing software packages.A linear transformation (i.e., output = gain × input + offset) is applied to normalize the DSM height from the float to unsigned integer (i.e., 0 -255).The level of degradation resulting from this normalization procedure is not evaluated in this study, although it is expected to be minor in terms of the PIV result.
The CC, as shown in Eq. ( 1), estimates the degree of similarity between images with a moving window (Richards 2013).In order to obtain sub-pixel accuracy, Least Squares Gaussian Fitting is applied to interpolate the highest correlation peak.
Where, C(u, v) is the degree of correlation in displacement (u, v); f 1 and f 2 are the interrogation areas.The CC of f 1 (j, k) and f 2 (ju, kv) is evaluated.The resulting CC C(u, v) indicates the 'degree of similarity' between the f 1 and f 2 patches over the offset range in displacement (u, v).The window size is set as 64 × 64 pixels.

DSM Subtraction
The change in terrain could be observed directly from the height difference between two epochs, as shown in Eq. ( 2).This process is named DSM subtraction in this study.
Where, dDSM is the height difference; DSM(t 1 ) is the elevation at time t 1 ; and DSM(t 2 ) is the elevation at time t 2 .With a series of height differences resulting from DSM subtraction, the geological uplift and river incision balance could be analyzed.

Verification
Due to its importance, the Taan River is administrated by the Water Resources Agency of Taiwan.In order to monitor channel changes, a pair of monuments is established on both sides of the channel for repeated cross section surveying.There are three pairs of these survey monuments with four campaigns conducted between March 2000 and December 2008 available for this study.These cross section surveys are typically carried out with total-station instruments operated by a ground crew.In principle, all feature points on which the slope changes should be surveyed.The average sampling distance is about 4 m while most distances are about 10 m.This dataset is used as a reference for aerial survey verification (WRA 2000(WRA , 2008)).

Integrated Analysis
The proposed method retrieves 3D displacements using 2D displacement vectors from PIV and height changes from DSMs.The exact location of the changed area at different periods is identified with the horizontal movement from PIV.The height value at (x 0 , y 0 , time 0 ) is obtained from DSM at the initial start location (x 0 , y 0 ).The height value at (x 0 , y 0 , time 1 ) is obtained from DSM at (x 0 + dx, y 0 + dy), where (dx, dy) are the 2D displacement vectors from PIV.The 3D movement can then be calculated with Eq. ( 3) laid out below.The dz derived this way carries both land cover horizontal and vertical movement information.This is then compared with the vertical change from DSM subtraction.Where, DSM(x 0 , y 0 ) is the elevation at point (x 0 , y 0 ); DSM(x 1 , y 1 ) is the elevation of point (x 1 , y 1 ); and dz is the height difference with integrated analysis.

River Morphology evolution Interpretation
Both the PIV analysis and DSM subtraction provide information about the geomorphologic change in the site.Examining the nature of these two approaches, PIV mainly provides the horizontal displacement vectors; while DSM subtraction directly presents the vertical changes.For this case study, PIV would describe the sequence of river channel erosion changes on the sides and DSM subtraction would measure the rate of incision.These two information sources could be amalgamated.The height changes in the area with large displacement identified with PIV are analyzed in this study.

STUDy AReA
The Taan River is located in Central Taiwan, originating west of the Dabajian Mountain of the Hsuehshan Mountain Series.The drainage area is about 758 km 2 with a trunk length of about 95.76 km.This basin is featured by intensive precipitation caused by typhoons every year due to its geographical setting.The study area is the uplift reach river corridor located between 27.7 and 28.7 km upstream from the estuary, nearby Lan-Shih Bridge, as shown in Fig. 2. Due to the geological uplift across the river induced by the Chi-Chi earthquake of 1999, and the river erosion caused by typhoons, the phenomenal geomorphic changes have attracted great academic interest.
The study site area is about 8.7 km 2 .The landscape is a river channel without heavy vegetation and built-up neighborhoods.The stratum exposed is the Pliocene Cholan Formation, composed of sandstone, siltstone, mudstone, and shale in a monotonous alternating sequence.The young sediment rock is poorly cemented with low erosion resistance.The Tungshih anticline is the major geological structure in this study area.The surface rupture terminated after it passed the Taan River Valley along the NE-SW direction.Two parallel ruptures were along the fold scarps of the Tungshih anticline associated with the synclinal bends at the base of the fold limbs.Across the Taan River Valley these two surface ruptures created a pop-up structure with a longitudinal distance of 2 km (Chen et al. 2007).The three pairs of cross section monuments, profiles 44, 44-1, and 45 in Fig. 2, are established for monitoring this geological feature.Huang et al. (2013) proposed that there are four discriminative stages according to landform characteristics in the river morphology evolution of the uplift area.In stage 1 (years 1999 to 2001), sediment was deposited in the barrier lake in the upstream region and no clear sign of channel incision was observed in the downstream region.In stage 2 (years 2001 to 2004), intense incision occurred in the exposed anticline axis bedrock and several shallow and narrow incised channels competed to be the final main channel.In stage 3 (years 2004 to 2007), bedrock erosion was much more active than in the previous stages because of the concentrated discharge in the main channel and the effect of knickpoint migration.The maximum incision in the main channel was larger than 10 m.In stage 4 (years 2007 to 2010), the bedrock incision remained active towards the upstream scarp and eventually carved through the uplift along the main channel.

exPeRIMenTAl ReSUlTS
Aerial photos in nine blocks, one block a year from 2001 -2009, were collected in this study.The flight altitude and other parameters of these photos are listed in Table 1.In years 2001 to 2007, photos were taken with an analogue aerial frame camera.The film is scanned at 20 μm resolution with each pixel corresponding to 31 cm on ground.In years 2008 and 2009, the aerial images were taken using a medium format digital camera along with airborne LiDAR survey.The pixel sizes were 9 and 6.8 μm for 2008 and 2009, respectively.Photogrammetric processing, including aerial triangulation, DSM matching, and orthoimage generation, were conducted with the commercial software Socet GXP (BAE Systems 2010).Twelve ground control points (GCPs) were collected with static GPS in 2009 to support the aerial triangulation.These GCPs are all road-side features that can be identified from most images and spread over a much larger area than that of the study site.The road has existed for a long time, remaining static enough for certain markers to be used as control points.These control points can be used from years 2001 to 2009.That is, in order to have sufficient ground control, the aerial photo block was expanded.However, despite these markers, the 2003, 2004, and 2006 blocks all suffered from large uncertainties in aerial triangulation.The reason for the difference between (or within) the 2003 and 2004 blocks is likely the large deviation between the aerial photos that comprise the block, as they were taken in different seasons.The 2006 block also does not have sufficient coverage.Therefore, the DSM and orthoimage in years 2001, 2002, 2005, 2007, 2008, and 2009 were used in this study.The orthoimage spatial resolution and DSM spacing are 1 m.The image-derived DSM and orthoimage are produced from the same images.Therefore, these two data sets have the same coordinate system.For the LiDAR-derived DSM and its accompanied orthoimage, the orthoimage is generated from images taken during the Li-DAR mission with a medium camera.Therefore, these two data sets also use the same coordinates.The DSMs in 2008 and 2009 are derived from airborne LiDAR.The DSMs in 2001, 2002, 2005, and 2007 are derived from aerial photogrammetry.The RMS of aerial triangulation is better than 1 pixel.The DSMs derived from aerial photogrammetry were compared using control points and the average margin of error was less than 2 m in height and the DSM error derived from LiDAR were less than 1 m in both the horizontal and vertical directions.

PIV Analysis
Taking year 2001 as a base, the PIV results between 2001 and 2007 are shown in Fig. 4. The vector color indicates the magnitude of movement while the vector direction shows the direction of displacement.Comparing orthoimage and DSM, the displacement vector trends are similar.These vectors indicate the displacements, with most displacements located on the river channel.However, orthoimage (Fig. 4a) provides more detail than DSM (Fig. 4b).For example, more details of channel changes in zone B can be identified with orthoimage than with DSM.This result is likely from the richer information content of image spectral data for the land cover change.
Figures 5 and 6 show the quantitative results from the PIV result for years 2002 to 2009 with respect to 2001 and the increments of accumulated horizontal movements.The interval of magnitude is 5 pixels.Comparing orthoimage and DSM, these two have the same horizontal movement trend and the percentage of large movement increases from 2002 -2009.In Fig. 5, the maximum movement from PIV with orthoimage was 50 pixels in 2002 while the maximum movement in 2009 increased to 90 pixels.This indicates that the change from 2001 -2009 is much larger than from 2001 -2002.In fact, 79% of the horizontal movement is below 5 pixels in 2002 and 70% below 5 pixels in 2003.When comparing the results before and after 2004, the number of movements in the 10 -15 pixel range is 389 in year 2003 (Fig. 5a), but the number of movements of the same range has significantly increased to 643 in year 2004 (Fig. 5b).This indicates that year 2004 is a turning point.The other   In comparing the PIV results from orthoimage and DSM, both describe the same trend.However, the former shows more detailed features than the latter.Figure 7 shows the yearly PIV results with orthoimage, with the magnitude of PIV displacement vector color coded.From the hot spot, those colored in red and other displacement patterns, years 2002 and 2003 are similar, while the rest show other patterns.This observation can be explained by the channel stabilizing after 2004.The maximum magnitude of displacement for each year is: 47.4 (2002), 50.4 (2003), 61.3 (2004), 76.7 (2005), 41.6 (2006), 81.2 (2007), 173.9 (2008), 156.6 (2009).This provides evidence for identifying the start of the third stage, in 2008.
In general, the erosion base level of this river segment rose as the result of earthquake uplift.Zone A has a welldefined channel region well delineated by PIV.The erosion in the channel is clearly identified.Comparatively, zones B and C have a complex landscape and thus PIV is affected by mixed terrain without a well-defined eroded channel.
Profiles along the three ground surveyed are drawn from PIV with orthoimage.The x-axis is the point location while the y-axis is the magnitude (length of displacement) from PIV.All of the relative displacements are compared to the year 2001.Figure 8

DSM Subtraction
When compared to the year 2001, the changes in DSMs in years 2002, 2005, 2007, 2008, and 2009 are derived as shown in Fig. 11.In year 2002, Fig. 11a, the surrounding area outside the bank is generally elevated.However, the area within the bank is generally lowered.Observed in 2005, Fig. 11b, pronounced channel incision has developed.This incision trend continues in 2007, 2008, and 2009, with the major channel at almost in the same location.Elapsed with time, the width of the incised channel is reduced with an increase in depth.The maximum incision is 12 m in 2009 relative to 2001.The other observation is zone B, the popup area, where the width of the main channel is narrower than in the other two zones.
The above mentioned channel evolution, in terms of its location and width, and also its incision, is supported with the height profiles.From profile 44, Fig. 12a, a  The width/depth ratios were 10 in profile 44, 5 in profile 44-1, and 13 in profile 45 for maximum vertical river incision.Profile 44-1 has a relatively smaller ratio, which indicates a relatively intensive incision force.

Verification
Figure 13a shows the ground surveyed profile 44.The channel is located between 200 -450 m along the river.This segment incised 10 m from years 2000 to 2008.This agrees with those observed in DSM subtraction, but with less detail.Figure 13b shows profile 44-1.The deep channel is located at a distance of -450 to -350 m.The channel incised 20 m from years 2000 to 2008.Comparing this with profile 44, there are more incisions than side erosion.The width/depth ratio is also 5, the same as observed from DSM. Figure 13c shows profile 45.The channel is located at a distance of 50 -150 m, and the width increased from 2000 -2008.The incision is about 15 m.These observations also agree with the PIV and DSM subtraction analysis.As a concluding remark for the verification, both the flow path and incision trend from these three profiles and the result from these two analyses agree with the ground survey.

Integrated Analysis
As described in section 2.4 the integrated analysis derives three dimensional displacements from the two dimen-sional information obtained from PIV with the DSM subtraction integration.Based on the material used for PIV there are two approaches.Integrated method 1 uses the orthoimage and integrated method 2 uses DSM.The statistics for these two approaches are summarized in Tables 2 and 3 for method 1 and Tables 4 and 5 for method 2. In Tables 2 and 4, the statistics are derived from the entire area; while Tables 3 and 5 are limited to the channel area.As the entire area contains both channel and river bank areas the maximum and minimum differences are overestimated.
The average horizontal movements for Tables 2 and 3 range from 22.93 -27.08 and 22.95 -27.30pixels, respectively.The average height difference trend for integrated method 1 and DSM subtraction is the same.
In method 2, the average horizontal movement of the channel area is smaller than the average horizontal movement of the entire area.incises to the deeper side slope.The minimum difference is -14.35 m using the integrated method 2 and -16.15 m using DSM subtraction in year 2009.This corresponds with the maximum height difference by channel erosion.The integrated method 2 describes the magnitude of river incision from years 2001 to 2009 that responds to the influence of the horizontal movement to get the exact height difference.
To compare the integrated and DSM subtraction methods (Fig. 14), the height difference from DSM subtraction is slightly larger than the results from the integrated method.Therefore, subtracting the DSM directly will produce a larg-er height difference.

Discussion
River morphology evolution is essentially a dynamic system.The preceding context describes the horizontal movement and vertical change in river morphology in the uplift area.The integration method shows that the average horizontal movement is about 20 m and the channel incised intensely up to 15 m in the channel area.According to the horizontal and vertical change analysis we can separate the river

ConClUSIonS
The river morphology after earthquake uplift was analyzed with aerial photos and a digital elevation model in this study.Both the horizontal and vertical direction was integrated to describe the changing river and the exact 3D movement of the channel area.The maximum vertical slip of the fault scarps was about 10 m, disturbing the dynamic equilibrium of the fluvial system.Based on the multi-staged orthoimages and DSMs, the morphology evolution process in the uplift reach was divided into three distinct stages.The major action in the study area is river erosion after earthquake uplift, especially river incision.The surveyed crosssection profile verifies the result computed by algorithm and terrain data.The correctness of the horizontal movement and vertical changes were demonstrated.The river morphology shows that earthquake uplift causes river incision.The change detection integration method is helpful in interpreting the river morphology evolution after earthquake uplift.This method uses orthoimage and DSM to analyse river morphology evolution.The traditional approaches discuss the horizontal and vertical movements separately.The proposed integration method combined horizontal movement and height difference simultaneously to show real 3D displacement.The contribution of this study is to provide a scheme for detecting both horizontal and vertical movements for river evolution in an integrated manner.After the analysis and verification, the method proposed in this study has demonstrated that the detected movement agrees with the understanding based on river morphology.It is expected that higher frequency and longer period DSMs and orthoimages would be helpful for better interpretation.Regarding further studies, besides improvements in the quantitative accuracy validation with reference data, such as manually digitized from images or obtained from field surveys, the physical understanding of the matched features in the PIV analysis is also a fundamental issue for fully understanding this methodology.

Fig. 2 .
Fig. 2. Overview of the study area, including rupture line, anticline, and surveyed cross-section profile.The study area divides three zones based on different geological structures.(Color online only) 2000 and December 2007 by the Third River Management Office, and August 2008, December 2008 by Water Resources Planning Institute.These two offices are both managed by the Water Resources Agency of Taiwan.The study area is divided into three zones with the earthquake surface ruptures and Tungshih anticline (Fig. 2).The water flows from north to south.Zone A has surface ruptures and with one channel in year 2001.Zone B has no consistent channel and the terrain is very flat.Zone C has two branches in year 2001.After ten years, zone A still has one main channel, but the location has changed.Zone B has a regular flow path with a narrow channel and another deep channel.Zone C has one main channel.The channel width is narrow and meanders to the NW.Figures 3a and b show the orthoimage in years 2001 and 2009.Figures 3c and d show the DSM in years 2001 and 2009.From these figures the river channel changes significantly between 2001 and 2009.

Fig. 4 .Fig. 5 .
Fig. 4. Results of PIV between 2001 and 2007: (a) orthoimages, (b) DSM.The color bar is rainbow color and the maximum displacement is 50 pixels.(Color online only) shows profile 44.The river channel location is from 200 -500 m in the x-axis.The maximum displacement in year 2002 was 20 m at 500 m, 30 m in year 2004, and increased to 60 m at 400 m in year 2009.The channel area expanded in the x-axis from 200 -500 m after 2002.These displacements in time indicate significant changes in the river channel.Profile 44-1 is shown in Fig. 9.The distance of 150 m begins changing and has maximum movement of 20 m in year 2002.After year 2002, the channel area expanded via its x-axis to 150 -250 m and the maximum horizontal movement was about 30 m.The channel location stabilized after 2003.Profile 45 is shown in Fig. 10.There are two channel branches in this profile, those of 0 -200 and 300 -400 m.The maximum movement is 40 m at 100 m and 30 m at 300 -400 m, which stabilized after 2005.Due to the image limits an incomplete profile 45 is also covered.
deep channel developed between 2002 and 2005.That is, the channel incised 10 m at 300 m.But, no further incision occurred at this location after 2005.Likewise, there was a 10m vertical incision at 220 m in 2005, as shown in Fig. 12b, the profile 44-1.This incision is deepened to 20 m in year 2007.In profile 45, Fig. 12c, 10 m incision at 110 m is observed in 2005, furthered to 15 m in year 2009.
morphology process into three stages.At the beginning (stage 1), the basic erosion level rose and the erosive power of river flow increased in the uplift area.The channel did not have a stable flow path and the horizontal movement increased year by year.The landform was incised by water and changed rapidly.In year 2005 (stage 2), the horizontal movement was up to the maximum value and the channel eroded down to the weak bedrock.The channel moved to the W direction and the slope was smooth in zone A. The vertical and side erosion were in balance with the channel.The main channel has an unchanged flow path and incised vertical direction in zone B. The vertical incision was the most important factor that affected the surface.The slope of the main flow channel gradually increased each year.The maximum incision in the main channel was more than 10 m.The channel stabilized and the river incision acted on the landform in zone C. The channel slope in zone C was same as that of zone A but the channel width was larger than zone A.Comparing the behavior in the three zones the rate of increase in width and incision are alike in zones A and C; while the incision is more prominent in zone B.

Table 1 .
Parameters of aerial photo.

Table 2 .
Comparison of integrated method (PIV from images and DSMs) and DSM subtraction for entire area.

Table 3 .
Comparison of integrated method (PIV from images and DSMs) and DSM subtraction for channel area.

Table 4 .
Comparison of integrated method (PIV from DSMs and DSMs) and DSM subtraction for entire area.

Table 5 .
Comparison of integrated method (PIV from DSMs and DSMs) and DSM subtraction for channel area.