Multiple Constraints Based Robust Matching of Poor-Texture Close-Range Images for Monitoring a Simulated Landslide

Landslides are one of the most destructive geo-hazards that can bring about great threats to both human lives and infrastructures. Landslide monitoring has been always a research hotspot. In particular, landslide simulation experimentation is an effective tool in landslide research to obtain critical parameters that help understand the mechanism and evaluate the triggering and controlling factors of slope failure. Compared with other traditional geotechnical monitoring approaches, the close-range photogrammetry technique shows potential in tracking and recording the 3D surface deformation and failure processes. In such cases, image matching usually plays a critical role in stereo image processing for the 3D geometric reconstruction. However, the complex imaging conditions such as rainfall, mass movement, illumination, and ponding will reduce the texture quality of the stereo images, bringing about difficulties in the image matching process and resulting in very sparse matches. To address this problem, this paper presents a multiple-constraints based robust image matching approach for poor-texture close-range images particularly useful in monitoring a simulated landslide. The Scale Invariant Feature Transform (SIFT) algorithm was first applied to the stereo images for generation of scale-invariate feature points, followed by a two-step matching process: feature-based image matching and area-based image matching. In the first feature-based matching step, the triangulation process was performed based on the SIFT matches filtered by the Fundamental Matrix (FM) and a robust checking procedure, to serve as the basic constraints for feature-based iterated matching of all the non-matched SIFT-derived feature points inside each triangle. In the following area-based image-matching step, the corresponding points of the non-matched features in each triangle of the master image were predicted in the homologous triangle of the searching image by using geometric constraints, followed by a refinement course with similarity constraint and robust checking. A series of temporal Single-Lens Reflex (SLR) and High-Speed Camera (HSC) stereo images captured during the simulated landslide experiment performed on the campus of Tongji University, Shanghai, were employed to illustrate the proposed method, and the dense and reliable image matching results were obtained. Finally, a series of temporal Digital Surface Models (DSM) in the landslide process were constructed using the close-range photogrammetry technique, followed by the discussion of the landslide volume changes and surface elevation changes during the simulation experiment.


Introduction
Landslides are one of the most destructive geo-hazards in mountain regions that can lead to great threats to human lives, as well as infrastructure damage in a very short time [1][2][3][4][5][6][7][8].This trend is likely to worsen in the future with the development of urbanization and economics, deforestation, and the increased regional rainfall in landslide-prone areas [9].Therefore, landslide monitoring and prediction has always been drawing great attention in both the engineering field and research community [10][11][12][13][14][15][16][17].
Generally speaking, landslides are usually triggered by some disastrous natural events (e.g., a great earthquake happening within a few seconds that brings tremendous damages, and a long period of heavy rainfall resulting from atrocious weather) or violent human activities (for example, excavating mountains using explosive when building a road or building) [18][19][20][21][22][23] with limited monitoring conditions, making it very difficult to get the real-time monitoring data for an in situ monitoring system during a typical landslide collapse process, thus reducing the feasibility in mechanism research of the landslide failure event [24].Sometimes an in situ monitoring system may record continuous data for very long time with no landslide happening, but the system could stop working under harsh weather conditions just before the landslide occurs.In this situation, the landslide simulation becomes a useful approach for obtaining critical parameters by allowing different types of monitoring sensors to work continuously during the landslide failure events in an controllable environment, being an important complement in landslide research [25][26][27][28][29].The landslide simulation platform also provides the possibility of employing photogrammetry, a non-contact and fast surface recording and reconstruction technology, in landslide research.The close-range photogrammetry in a landslide simulation can track and record the geometric deformation processes before and during the landslide failure event, providing an important input for understanding the landslide mechanism and evaluating the effects [30][31][32].
Image matching is a critical component in the stereo image processing to obtain the 3D information from image space, and this problem has long been studied in both photogrammetry and computer vision [33][34][35].A number of algorithms have been proposed on stereo image matching in different aspects [36][37][38][39][40][41][42][43].Lhuillier and Quan [44] proposed a quasi-dense matching algorithm between images based on the match propagation principle with discrete 2D gradient disparity limit and the uniqueness constraint.Furukawa and Ponce [45] implemented stereopsis as a match, expand, and filter procedure enforcing local photometric consistency and global visibility constraints.Zhu and Deng [46] used gradient orientation selective cross-correlation excluding the wrong points from correlation to image matching.Stentoumis et al. [47] provided an accurate dense matching using a local adaptive multi-cost approach, typical in hierarchical matching.Zhu et al. [48] and Song et al. [49] proposed a propagation-based stereo-matching algorithm.The former was under the dynamic triangle constraint; the latter constructed a line segment region for each pixel with local color and connectivity constraints.Stumpf et al. [50] used the algorithm implemented in Co-registration of Optically Sensed Images and Correlation (COSI-Corr) [51], which is based on phase correlation in the frequency domain, to achieve sub-pixel image correlation.These methods work very well on stereo images with good texture and can get dense matches, but only sparse conjugate points could be obtained for poor-texture images.
Recently, the Semi Global Matching (SGM) algorithm was proposed by Hirschmuller [52] and was first employed in HRSC Mars Express images, as well as in structured environments [53,54], for stereo matching before it was widely used in the computer vision society.As a dense image-matching algorithm, SGM has been applied to aerial full frame images and aerial pushbroom images in densely-populated regions for 3D reconstruction [55], high building roofs and facades for 3D models [56], as well as high-resolution satellite images for high mountain and glacier mapping [57].This algorithm has now been extended in many aspects including hardware implementation (CPU/GPU/FPGA) for real-time 3D mapping and navigation purposes [58,59].
Some image matching research has also been carried out for poor-texture stereo images.Wu et al. [60] presented an image matching method by integration of points and edges for dense image-matching on poor-texture images and good matches were obtained from space-borne, airborne, and terrestrial images with poor textures.Chen et al. [61] proposed a line-based matching method in low-texture area for high-resolution images with a narrow field-of-view camera or a short baseline.Bulatov et al. [62] used multi-view images and developed a dense matching method supported by triangular meshes, which is suitable for poor-texture images.
The above mentioned methods are effective for the stereo images with relatively poor textures, and edge information or multi-view images is essential for the reliable and dense matching in these methods.For the close-range stereo images obtained from landslide simulation observation (shown in Figure 1, two stereo image pairs captured by Single-Lens Reflex (SLR) and High-Speed Cameras (HSC), respectively), it is very difficult to find edge information on the landslide body.In these images, the landslide mass, itself, composed of soil, sand, and pebbles, shows similar surface characteristics, and other factors during the imaging process such as rainfall, mass movement, illumination, and ponding further reduce the texture quality in the regions to be matched, resulting in similar, homogeneous, low, or no textures.Figure 1a,b show similar or homogeneous textures in the landslide body collected by the SLR cameras before the landslide failure event.Figure 1c,d present the moment during the failure event when there are also many blank spots without textures induced by the reflection of water areas, which usually result in ambiguities during the matching process.Commercial software usually fails when doing dense image matching for this kind of images.For example, PMS (PhotoModeler Scanner, Vancouver, BC, Canada) [63] is a renowned software known for close-range image analysis and applications, but it requires images with good texture, and the matching result is not satisfactory for poor quality images.In this case, the exiting matching methods and software cannot generate good matches, and a new matching method is needed for effective image matching of the close-range landslide simulation images.
Remote Sens. 2016, 8, 396 3 of 24 airborne, and terrestrial images with poor textures.Chen et al. [61] proposed a line-based matching method in low-texture area for high-resolution images with a narrow field-of-view camera or a short baseline.Bulatov et al. [62] used multi-view images and developed a dense matching method supported by triangular meshes, which is suitable for poor-texture images.
The above mentioned methods are effective for the stereo images with relatively poor textures, and edge information or multi-view images is essential for the reliable and dense matching in these methods.For the close-range stereo images obtained from landslide simulation observation (shown in Figure 1, two stereo image pairs captured by Single-Lens Reflex (SLR) and High-Speed Cameras (HSC), respectively), it is very difficult to find edge information on the landslide body.In these images, the landslide mass, itself, composed of soil, sand, and pebbles, shows similar surface characteristics, and other factors during the imaging process such as rainfall, mass movement, illumination, and ponding further reduce the texture quality in the regions to be matched, resulting in similar, homogeneous, low, or no textures.Figure 1a,b show similar or homogeneous textures in the landslide body collected by the SLR cameras before the landslide failure event.Figure 1c,d present the moment during the failure event when there are also many blank spots without textures induced by the reflection of water areas, which usually result in ambiguities during the matching process.Commercial software usually fails when doing dense image matching for this kind of images.For example, PMS (PhotoModeler Scanner, Vancouver, BC, Canada) [63] is a renowned software known for close-range image analysis and applications, but it requires images with good texture, and the matching result is not satisfactory for poor quality images.In this case, the exiting matching methods and software cannot generate good matches, and a new matching method is needed for effective image matching of the close-range landslide simulation images.This paper presents a multiple constraints-based robust matching approach for poor-texture close-range images in simulated landslide monitoring, as well as imaged landslide surface changes caused by slope failure process.The proposed image matching approach includes two steps, feature-based image matching mainly constrained by triangulation, and area-based image matching confined by triangulation, regional Perspective Transformation (PT), and Epipolar Line (EL).The landslide simulation platform and experiment is first introduced in Section 2. The methodology is then presented in Section 3, including feature extraction and the two-step image matching process.Taking SLR stereo image pairs as an example, the overall workflow and the respective principles are depicted in this section to help readers understand the details.In Section 4, the SLR and HSC stereo image series are processed for image matching and the matching results are evaluated.In Section 5, the Digital Surface Model (DSM) series are first generated, then the experimental results obtained by the close-range photogrammetric systems are reported before and during the landslide slope failure process, including the landslide volume change analysis and the surface elevation changes.Section 6 gives the discussion.Finally, Section 7 describes the conclusions.This paper presents a multiple constraints-based robust matching approach for poor-texture close-range images in simulated landslide monitoring, as well as imaged landslide surface changes caused by slope failure process.The proposed image matching approach includes two steps, feature-based image matching mainly constrained by triangulation, and area-based image matching confined by triangulation, regional Perspective Transformation (PT), and Epipolar Line (EL).The landslide simulation platform and experiment is first introduced in Section 2. The methodology is then presented in Section 3, including feature extraction and the two-step image matching process.Taking SLR stereo image pairs as an example, the overall workflow and the respective principles are depicted in this section to help readers understand the details.In Section 4, the SLR and HSC stereo image series are processed for image matching and the matching results are evaluated.In Section 5, the Digital Surface Model (DSM) series are first generated, then the experimental results obtained by the close-range photogrammetric systems are reported before and during the landslide slope failure process, including the landslide volume change analysis and the surface elevation changes.Section 6 gives the discussion.Finally, Section 7 describes the conclusions.

Landslide Platform Set Up and Simulation Experiment
In this study, a scaled-down simulation platform was constructed on the campus of Tongji University in Shanghai, China, to reproduce a landslide-prone slope near the small town of Taziping in Sichuan Province, Western China, where the 2008 Wenchuan Earthquake left great threats of susceptible landslides due to the loose soil layers afterwards [27], shown as in Figure 2. The dimension of the landslide body was 6 m ˆ1.5 m ˆ3 m (length, width, and height), with three slope sections featuring inclinations of 30 ˝, 15 ˝, and 5 ˝.For real-time monitoring and early warning purposes, this platform was designed and implemented to include an artificial rainfall system, a sensor network, a subsystem for data collection and communication, a data server for storage, and a screen panel for visualization, for more details readers can refer to Qiao et al. [27] and Lu et al. [64].

Landslide Platform Set Up and Simulation Experiment
In this study, a scaled-down simulation platform was constructed on the campus of Tongji University in Shanghai, China, to reproduce a landslide-prone slope near the small town of Taziping in Sichuan Province, Western China, where the 2008 Wenchuan Earthquake left great threats of susceptible landslides due to the loose soil layers afterwards [27], shown as in Figure 2. The dimension of the landslide body was 6 m × 1.5 m × 3 m (length, width, and height), with three slope sections featuring inclinations of 30°, 15°, and 5°.For real-time monitoring and early warning purposes, this platform was designed and implemented to include an artificial rainfall system, a sensor network, a subsystem for data collection and communication, a data server for storage, and a screen panel for visualization, for more details readers can refer to Qiao et al. [27] and Lu et al. [64].The sensor network contains contact sensors that were installed in the landslide mass and used to record the environmental conditions to derive the geotechnical parameters of the landslide body, and the detailed research and discussion in this aspect can be found in [28,64].On the other hand, the non-contact/imaging sensors (cameras and video) capture the geometric changes in the slope surface during a simulated landslide deformation process.A stereo pair of NIKON D200 SLR cameras was deployed to collect the landslide surface changes during the entire simulation experiment at a relatively low frequency (a few seconds).To interpret the transient slope-failure process in detail, a HSC system, composed of a pair of synchronized DALSA Falcon 4M60 high-speed cameras that can capture images at high frequency (up to 62 Hz) with a synchronization accuracy of 0.1 ms, was employed.The cameras are designed mainly to capture the surface movement of the landslide body, and the locations of both the cameras and the waterway are restrained by the experiment site.In the landslide simulation experiment, the total amount of rainfall is huge, so the waterway is designed with holes (see Figure 2) to filter the water and guide the landslide mass to a designated area; thus, it has very little influence on the landslide process.
A set of well-distributed marked points (shown in Figure 1) that were fixed on the facility and measured by a total station were employed as GCPs to provide reference for recovery of the orientation parameters of the camera systems.
In this experiment, the SLR camera system was employed to record the entire landslide process from 12:25:00 to 14:50:00 at a frequency of 6 frames/min, and a total of 737 stereo image pairs were obtained.During the final failure event at 14:27:22, the low-frequency SLR camera system was unable to capture the detailed changes that occurred within only a few seconds, and the HSC system was started, and a total of 95 stereo image pairs were captured at a frequency of 20 frames/s.The experiment settings of the two stereo camera systems are listed in Table 1.The sensor network contains contact sensors that were installed in the landslide mass and used to record the environmental conditions to derive the geotechnical parameters of the landslide body, and the detailed research and discussion in this aspect can be found in [28,64].On the other hand, the non-contact/imaging sensors (cameras and video) capture the geometric changes in the slope surface during a simulated landslide deformation process.A stereo pair of NIKON D200 SLR cameras was deployed to collect the landslide surface changes during the entire simulation experiment at a relatively low frequency (a few seconds).To interpret the transient slope-failure process in detail, a HSC system, composed of a pair of synchronized DALSA Falcon 4M60 high-speed cameras that can capture images at high frequency (up to 62 Hz) with a synchronization accuracy of 0.1 ms, was employed.The cameras are designed mainly to capture the surface movement of the landslide body, and the locations of both the cameras and the waterway are restrained by the experiment site.In the landslide simulation experiment, the total amount of rainfall is huge, so the waterway is designed with holes (see Figure 2) to filter the water and guide the landslide mass to a designated area; thus, it has very little influence on the landslide process.
A set of well-distributed marked points (shown in Figure 1) that were fixed on the facility and measured by a total station were employed as GCPs to provide reference for recovery of the orientation parameters of the camera systems.
In this experiment, the SLR camera system was employed to record the entire landslide process from 12:25:00 to 14:50:00 at a frequency of 6 frames/min, and a total of 737 stereo image pairs were obtained.During the final failure event at 14:27:22, the low-frequency SLR camera system was unable to capture the detailed changes that occurred within only a few seconds, and the HSC system was started, and a total of 95 stereo image pairs were captured at a frequency of 20 frames/s.The experiment settings of the two stereo camera systems are listed in Table 1.Examples of the stereo images with poor-texture collected during the landslide simulation experiment are shown in Figure 1.

Methodology
Figure 3 illustrates the workflow of poor-texture close-range image processing proposed in this research.It is composed of two parts, the image geometry determination part before the simulation experiment (shown as green box in Figure 3) for 3D geo-referencing of the close-range photogrammetric systems, and the image processing and matching part (shown as blue box in Figure 3) for surface modeling during the landslide simulation experiment.Before the experiment, camera calibration is first performed and the GCPs are set up and measured for determination of Internal Orientation (IO) and External Orientation (EO) elements, which will be further employed in the process of image matching and 3D point calculation.After the landslide experiment, Scale Invariant Feature Transform (SIFT) algorithm is first applied to the poor-texture close-range stereo images for feature extraction, and then the images are processed for the two-step matching procedure, namely feature-based image-matching and area-based image matching; finally, the 3D points are obtained with orientation parameters (IOs and EOs) and a DSM is generated.Examples of the stereo images with poor-texture collected during the landslide simulation experiment are shown in Figure 1.

Methodology
Figure 3 illustrates the workflow of poor-texture close-range image processing proposed in this research.It is composed of two parts, the image geometry determination part before the simulation experiment (shown as green box in Figure 3) for 3D geo-referencing of the close-range photogrammetric systems, and the image processing and matching part (shown as blue box in Figure 3) for surface modeling during the landslide simulation experiment.Before the experiment, camera calibration is first performed and the GCPs are set up and measured for determination of Internal Orientation (IO) and External Orientation (EO) elements, which will be further employed in the process of image matching and 3D point calculation.After the landslide experiment, Scale Invariant Feature Transform (SIFT) algorithm is first applied to the poor-texture close-range stereo images for feature extraction, and then the images are processed for the two-step matching procedure, namely feature-based image-matching and area-based image matching; finally, the 3D points are obtained with orientation parameters (IOs and EOs) and a DSM is generated.At the very beginning, the camera calibration was carried out by PMS Software using the markers automatically identified from the calibration plates, and the initial IOs were then obtained using a self-calibrated bundle adjustment.The GCPs were pasted on the steel framework of the simulation platform, shown in Figure 1 as the red box.The GCP coordinates in object space were determined by a Sokkia total station with an accuracy of ±1 mm, and those in image space were At the very beginning, the camera calibration was carried out by PMS Software using the markers automatically identified from the calibration plates, and the initial IOs were then obtained using a self-calibrated bundle adjustment.The GCPs were pasted on the steel framework of the simulation platform, shown in Figure 1 as the red box.The GCP coordinates in object space were determined by a Sokkia total station with an accuracy of ˘1 mm, and those in image space were measured in PMS.The IOs and EOs were finally solved and refined simultaneously using PMS.In the following sub-sections, the multi-constraints based robust matching approach is described in detail based on a randomly selected SLR stereo image pair (Figure 1).

Feature Extraction
The feature points are needed for the image matching step and the feature extraction method is very important for the generation of dense point matching results.Here, we compared the performance of three commonly used feature extraction operators: SIFT, SURF, and STAR [37,[65][66][67][68].After trial and error, the maximum number of detected feature points of each method with the best settings is listed in Table 2, from which one can observe that the SIFT operator generates the most features.In fact, visual inspection also shows that the SIFT detector generated feature points everywhere, whereas other operators did not.So the SIFT operator was applied to the entire landslide body region and large amount of feature points (48,735 in left image and 21,548 in right image, respectively) were extracted from the stereo image pair, illustrated as Figure 4. Here, the difference of number of extracted features is mainly caused by the image quality.Compared with the left image, the right image is more blurred due to the imaging conditions.This large number of feature points is enough to provide the seeds for the image matching in the landslide simulation application.The next step is to match the feature points through the feature-to-area matching process.
detail based on a randomly selected SLR stereo image pair (Figure 1).

Feature Extraction
The feature points are needed for the image matching step and the feature extraction method is very important for the generation of dense point matching results.Here, we compared the performance of three commonly used feature extraction operators: SIFT, SURF, and STAR [37,[65][66][67][68].After trial and error, the maximum number of detected feature points of each method with the best settings is listed in Table 2, from which one can observe that the SIFT operator generates the most features.In fact, visual inspection also shows that the SIFT detector generated feature points everywhere, whereas other operators did not.So the SIFT operator was applied to the entire landslide body region and large amount of feature points (48,735 in left image and 21,548 in right image, respectively) were extracted from the stereo image pair, illustrated as Figure 4. Here, the difference of number of extracted features is mainly caused by the image quality.Compared with the left image, the right image is more blurred due to the imaging conditions.This large number of feature points is enough to provide the seeds for the image matching in the landslide simulation application.The next step is to match the feature points through the feature-to-area matching process.

Feature-Based Image-Matching
Figure 5 shows the flowchart of the feature-based matching process.In the first-level feature-based matching, the feature points derived from SIFT operator were matched by using the SIFT descriptor along the EL, followed by Fundemental Matrix (FM) filtering and robust checking process to distinguish the matched features and mismatches.Due to the poor and similar textures in the landslide body region, the SIFT descriptor is not capable of producing satisfying dense matches in a relatively large region [69,70], so more specific geometric constraints should be provided.A triangulation process was then performed based on the matched feature points to serve as the basic constraints, and for all the non-matched SIFT-derived feature points inside each triangle, the

Feature-Based Image-Matching
Figure 5 shows the flowchart of the feature-based matching process.In the first-level feature-based matching, the feature points derived from SIFT operator were matched by using the SIFT descriptor along the EL, followed by Fundemental Matrix (FM) filtering and robust checking process to distinguish the matched features and mismatches.Due to the poor and similar textures in the landslide body region, the SIFT descriptor is not capable of producing satisfying dense matches in a relatively large region [69,70], so more specific geometric constraints should be provided.A triangulation process was then performed based on the matched feature points to serve as the basic constraints, and for all the non-matched SIFT-derived feature points inside each triangle, the second-level feature-based matching result was obtained using the same matching and refinement strategy (SIFT operator, FM filtering, and robust checking) as that in the first-level process.The triangulation process was performed for all the refined matched feature points of the first-and second-levels feature-based matching, and a new level of feature-based matching will then be iterated.During this procedure, the number of refined matched feature points will be growing gradually to serve as the feature-based feature matching result, and this process will continue until the threshold of the triangle sizes is satisfied.
Remote Sens. 2016, 8, 396 7 of 24 second-level feature-based matching result was obtained using the same matching and refinement strategy (SIFT operator, FM filtering, and robust checking) as that in the first-level process.The triangulation process was performed for all the refined matched feature points of the first-and second-levels feature-based matching, and a new level of feature-based matching will then be iterated.During this procedure, the number of refined matched feature points will be growing gradually to serve as the feature-based feature matching result, and this process will continue until the threshold of the triangle sizes is satisfied.According to the EL principle, in stereo vision, for each point observed in one image, the corresponding point must be observed on the corresponding EL in the other image determined by the orientation parameters [71,72].For each of the SIFT feature points in left image, the SIFT descriptor was used for matching along the corresponding EL in the right image that was calculated by the IO and EO parameters.Here, the distance between the corresponding point and this EL was defined to be two pixels, considering the possible orientation errors induced by GCP measurement and camera calibration.After this process, the SIFT matches were obtained and the outliers were then filtered out by the FM that relates the corresponding points in stereo images.In computer vision, the FM F satisfies the condition that for any pair of corresponding points ↔ in the two images [72]: According to the EL principle, in stereo vision, for each point observed in one image, the corresponding point must be observed on the corresponding EL in the other image determined by the orientation parameters [71,72].For each of the SIFT feature points in left image, the SIFT descriptor was used for matching along the corresponding EL in the right image that was calculated by the IO and EO parameters.Here, the distance between the corresponding point and this EL was defined to be two pixels, considering the possible orientation errors induced by GCP measurement and camera calibration.After this process, the SIFT matches were obtained and the outliers were then filtered out by the FM that relates the corresponding points in stereo images.In computer vision, the FM F satisfies the condition that for any pair of corresponding points x Ø x 1 in the two images [72]: In this research, F can be estimated by the known IO and EO parameters of the stereo images [72].The result of SIFT matching and mismatches filtering is presented in Figure 6, where all the points in the green quadrangle (Figure 6a,b) were the initial matched features by the SIFT operator along EL.The pseudo corresponding points in Figure 6 were first filtered out by FM.In view of the orientation errors, the left part of Equation (1) x 1T Fx with normalized point coordinates for x and x' was set to be less than 0.01 during the filtering process.The SIFT matches that do not satisfy condition will be regarded as mismatches and will be filtered out, shown as the blue points in Figure 6 (LB and RB).In this research, F can be estimated by the known IO and EO parameters of the stereo images [72].The result of SIFT matching and mismatches filtering is presented in Figure 6, where all the points in the green quadrangle (Figure 6a,b) were the initial matched features by the SIFT operator along EL.The pseudo corresponding points in Figure 6 were first filtered out by FM.In view of the orientation errors, the left part of Equation (1) with normalized point coordinates for x and x' was set to be less than 0.01 during the filtering process.The SIFT matches that do not satisfy condition will be regarded as mismatches and will be filtered out, shown as the blue points in Figure 6 (LB and RB).The next step is to refine the matched feature points using more robust constraint, such as a selected matching cost function.A Normalized Correlation Coefficient (NCC) [71] was employed to examine the similarity of the matched feature points.Compared with other matching cost functions such as Sum of Absolute Difference s (SAD) and Census, NCC is statistically optimal for dealing with Gaussian noise [73] contained in the landslide simulation images.For each pair of matched feature points in the two images, NCC was calculated with a window size of 17-pixels by 17-pixels.The threshold of the NCC was set to 0.65, relatively low due to the fact that the candidates had been refined by FM.The yellow points in Figure 6 (LA and RA, LB and RB) are the pseudo corresponding points detected using NCC with the given threshold.It should be noted that, here, the purpose of FM filtering and robust checking is to remove all the pseudo corresponding points, leaving the real matches for first-level matching to serve as further constraints, so in this process some of the real matches may also be removed due to the rigid parameters.The remainder after this process are marked as red points in Figure 6.
The first-level refined matched feature points were then triangulated using the Delaunay criterion [74] for further matching.Figure 7 shows the Delaunay triangulation result of the first-level matched features.For each triangle satisfying the size requirement (for example, the yellow triangle in Figure 7), the second-level feature-based matching process was performed, during which all the none matched SIFT feature points were refined using the same matching-filtering-robust checking procedure as the first-level process.The triangle size requirement mainly relates to the lengths of its edges, for computation efficiency, we define two thresholds of the image coordinate differences in X and Y directions, respectively (here both five pixels), for every two of the three vertices.Only the triangle with all the coordinate differences larger than the corresponding thresholds will be employed for the next-level feature-based matching.The next step is to refine the matched feature points using more robust constraint, such as a selected matching cost function.A Normalized Correlation Coefficient (NCC) [71] was employed to examine the similarity of the matched feature points.Compared with other matching cost functions such as Sum of Absolute Difference s (SAD) and Census, NCC is statistically optimal for dealing with Gaussian noise [73] contained in the landslide simulation images.For each pair of matched feature points in the two images, NCC was calculated with a window size of 17-pixels by 17-pixels.The threshold of the NCC was set to 0.65, relatively low due to the fact that the candidates had been refined by FM.The yellow points in Figure 6 (LA and RA, LB and RB) are the pseudo corresponding points detected using NCC with the given threshold.It should be noted that, here, the purpose of FM filtering and robust checking is to remove all the pseudo corresponding points, leaving the real matches for first-level matching to serve as further constraints, so in this process some of the real matches may also be removed due to the rigid parameters.The remainder after this process are marked as red points in Figure 6.
The first-level refined matched feature points were then triangulated using the Delaunay criterion [74] for further matching.Figure 7 shows the Delaunay triangulation result of the first-level matched features.For each triangle satisfying the size requirement (for example, the yellow triangle in Figure 7), the second-level feature-based matching process was performed, during which all the none matched SIFT feature points were refined using the same matching-filtering-robust checking procedure as the first-level process.The triangle size requirement mainly relates to the lengths of its edges, for computation efficiency, we define two thresholds of the image coordinate differences in X and Y directions, respectively (here both five pixels), for every two of the three vertices.Only the triangle with all the coordinate differences larger than the corresponding thresholds will be employed for the next-level feature-based matching.
Remote Sens. 2016, 8, 396 9 of 24 The matching-filtering-robust checking process was iterated to obtain as many matched features as possible, at each iteration, all the refined matched feature points were used for triangulation for next level matching process.This process ended when no triangle meets the size requirement.In this way, the feature points in the landslide body can be sparsely matched, the result of which is presented in Figure 8.

Area-Based Image-Matching
Generally, a large collection of local feature points could be generated from an image by the SIFT algorithm, and these feature points are good candidates for feature matching to produce DSM.During the feature-based matching step, only part of the corresponding feature points from SIFT algorithm in the stereo pair were correctly matched.The area-based image matching step here is to try to match the remaining none-matched features by using a multiple-constraints assisted matching method.
Figure 9 shows the flowchart of area-based image matching process for both the left and right images (see the rhombuses).In this process, the sparsely-matched feature points derived from Section 3.2 were first triangulated using the Delaunay principle in the stereo image pair, respectively.Then, the area-based image-matching procedure was implemented for both the right image and left image, shown as the blue and red boxes in Figure 9.In each situation, the master image and searching image were first determined, and the corresponding points of the non-matched feature points in each triangle of the master image were then predicted in the homologous triangle of the searching image by using geometric constraint, followed by a refinement course with a similarity constraint and robust checking.Finally, the repeated matches were removed by using a redundancy checking procedure, and the area-based matching points were then obtained.The matching-filtering-robust checking process was iterated to obtain as many matched features as possible, at each iteration, all the refined matched feature points were used for triangulation for next level matching process.This process ended when no triangle meets the size requirement.In this way, the feature points in the landslide body can be sparsely matched, the result of which is presented in Figure 8.The matching-filtering-robust checking process was iterated to obtain as many matched features as possible, at each iteration, all the refined matched feature points were used for triangulation for next level matching process.This process ended when no triangle meets the size requirement.In this way, the feature points in the landslide body can be sparsely matched, the result of which is presented in Figure 8.

Area-Based Image-Matching
Generally, a large collection of local feature points could be generated from an image by the SIFT algorithm, and these feature points are good candidates for feature matching to produce DSM.During the feature-based matching step, only part of the corresponding feature points from SIFT algorithm in the stereo pair were correctly matched.The area-based image matching step here is to try to match the remaining none-matched features by using a multiple-constraints assisted matching method.
Figure 9 shows the flowchart of area-based image matching process for both the left and right images (see the rhombuses).In this process, the sparsely-matched feature points derived from Section 3.2 were first triangulated using the Delaunay principle in the stereo image pair, respectively.Then, the area-based image-matching procedure was implemented for both the right image and left image, shown as the blue and red boxes in Figure 9.In each situation, the master image and searching image were first determined, and the corresponding points of the non-matched feature points in each triangle of the master image were then predicted in the homologous triangle of the searching image by using geometric constraint, followed by a refinement course with a similarity constraint and robust checking.Finally, the repeated matches were removed by using a redundancy checking procedure, and the area-based matching points were then obtained.

Area-Based Image-Matching
Generally, a large collection of local feature points could be generated from an image by the SIFT algorithm, and these feature points are good candidates for feature matching to produce DSM.During the feature-based matching step, only part of the corresponding feature points from SIFT algorithm in the stereo pair were correctly matched.The area-based image matching step here is to try to match the remaining none-matched features by using a multiple-constraints assisted matching method.
Figure 9 shows the flowchart of area-based image matching process for both the left and right images (see the rhombuses).In this process, the sparsely-matched feature points derived from Section 3.2 were first triangulated using the Delaunay principle in the stereo image pair, respectively.Then, the area-based image-matching procedure was implemented for both the right image and left image, shown as the blue and red boxes in Figure 9.In each situation, the master image and searching image were first determined, and the corresponding points of the non-matched feature points in each triangle of the master image were then predicted in the homologous triangle of the searching image by using geometric constraint, followed by a refinement course with a similarity constraint and robust checking.Finally, the repeated matches were removed by using a redundancy checking procedure, and the area-based matching points were then obtained.This area-based image-matching algorithm is essentially an area-based correspondence approach assisted by multiple constraints for the matching of non-matched feature points in each triangle of the master image.Figure 10 shows an example of the area-based image matching process, in which the left image is the master image and the right image is the searching image, and here only the sparsely matched features are plotted out as red dots in the green boxes for visual effect, followed by a triangulation procedure that generates triangles as basic constraints [48].Figure 10A,B illustrate the details of the area-based matching process for the non-matched feature points within a triangle (with yellow edges).For a feature point Pm in the master triangle, the region containing the conjugate point in searching triangle is first predicted by using geometric constraints including EL and PT.Theoretically, the conjugate point Pc should lie on the EL L derived from the orientation parameters of the two images.Additionally, Pm and its corresponding point should also follow the PT principle [75] that can be realized by the existing sparsely-matched feature points.In this study, the PT matrix [72] between Pm and its corresponding point is determined by six corresponding This area-based image-matching algorithm is essentially an area-based correspondence approach assisted by multiple constraints for the matching of non-matched feature points in each triangle of the master image.Figure 10 shows an example of the area-based image matching process, in which the left image is the master image and the right image is the searching image, and here only the sparsely matched features are plotted out as red dots in the green boxes for visual effect, followed by a triangulation procedure that generates triangles as basic constraints [48].Figure 10A,B illustrate the details of the area-based matching process for the non-matched feature points within a triangle (with yellow edges).For a feature point Pm in the master triangle, the region containing the conjugate point in searching triangle is first predicted by using geometric constraints including EL and PT.Theoretically, the conjugate point Pc should lie on the EL L derived from the orientation parameters of the two images.Additionally, Pm and its corresponding point should also follow the PT principle [75] that can be realized by the existing sparsely-matched feature points.In this study, the PT matrix [72] between Pm and its corresponding point is determined by six corresponding feature points including the three vertices of the corresponding triangle and the vertices of the three triangles that connect with the corresponding triangle.The predicted conjugate point by PT is Ppt, the projection of which to the EL is Pp, shown in Figure 10C.The searching window for corresponding point of Pm is then determined as the rectangle in Figure 10C with Pp as its center, one side parallel to L, with a given side length (here 10 pixels).After applying the geometric constraint, the searching area for a matching point is restricted in a reasonable searching window, improving the matching reliability and efficiency, followed by the similarity constraint with NCC to determine the corresponding point.NCC is, thus, estimated between feature point Pm and each pixel within the searching window (rectangle) with a given window size (here also 17-pixels by 17-pixels), and the point corresponding to the maximum NCC value is then obtained.The corresponding point of Pm can be determined if the maximum NCC is larger than the given threshold (here is set to 0.78, larger than that in feature-based matching process according to the related reference [76] and manual interactive checking of the matching results), shown as Pc in Figure 10C.After that, a robust checking step was applied by bilateral matching that inversely determines the matching point of Pc in the previous master triangle (Figure 10A) using the same geometric and similarity constraints as described above, to check whether the matching result is Pm, which helped to remove the erroneous matches and enhance the robustness of the area-based matching result.
Remote Sens. 2016, 8, 396 11 of 24 one side parallel to L, with a given side length (here 10 pixels).After applying the geometric constraint, the searching area for a matching point is restricted in a reasonable searching window, improving the matching reliability and efficiency, followed by the similarity constraint with NCC to determine the corresponding point.NCC is, thus, estimated between feature point Pm and each pixel within the searching window (rectangle) with a given window size (here also 17-pixels by 17-pixels), and the point corresponding to the maximum NCC value is then obtained.The corresponding point of Pm can be determined if the maximum NCC is larger than the given threshold (here is set to 0.78, larger than that in feature-based matching process according to the related reference [76] and manual interactive checking of the matching results), shown as Pc in Figure 10C.After that, a robust checking step was applied by bilateral matching that inversely determines the matching point of Pc in the previous master triangle (Figure 10A) using the same geometric and similarity constraints as described above, to check whether the matching result is Pm, which helped to remove the erroneous matches and enhance the robustness of the area-based matching result.The area-based image-matching procedure was implemented in both the left image and right image by switching the master image and searching image which, on one hand, produced as many densely matched points as possible, and on the other, generated some repeated matches.A redundancy checking procedure was applied to remove the reduplicate matched points for those within one pixel difference.The final area-based matching result was then obtained, shown as in Figure 11.The area-based image-matching procedure was implemented in both the left image and right image by switching the master image and searching image which, on one hand, produced as many densely matched points as possible, and on the other, generated some repeated matches.A redundancy checking procedure was applied to remove the reduplicate matched points for those within one pixel difference.The final area-based matching result was then obtained, shown as in Figure 11.
The area-based image-matching procedure was implemented in both the left image and right image by switching the master image and searching image which, on one hand, produced as many densely matched points as possible, and on the other, generated some repeated matches.A redundancy checking procedure was applied to remove the reduplicate matched points for those within one pixel difference.The final area-based matching result was then obtained, shown as in Figure 11.

Results of Image Matching
In this paper, the landslide surface deformation description and analysis during the simulation process recorded by the stereo images using the proposed methodology was focused on.The landslide surface changes mainly appeared before and during the failure process, and there were no obvious changes occurred afterwards, so in this research we will only process and analyze the 568 SLR stereo image pairs acquired before the final failure event (12:25:00 to 14:27:20, hereafter the pre-failure stage) and 95 HSC stereo image pairs recorded during the failure process (14:27:22:000 to 14:27:26:700, hereafter the failure stage).

Image Processing for Stereo Image Series Recording Landslide Simulation Experiment
The low-quality landslide images were first preprocessed by using image enhancement techniques, such as a Wallis filter [77], in the landslide region to enhance and sharpen the texture patterns and increase the signal-to-noise ratio.The proposed image matching approach was then implemented on the stereo SLR and HSC image series to produce matched points.Figure 12 gives an example of the area-based image matching result for a pair of HSC stereo images.

Results of Image Matching
In this paper, the landslide surface deformation description and analysis during the simulation process recorded by the stereo images using the proposed methodology was focused on.The landslide surface changes mainly appeared before and during the failure process, and there were no obvious changes occurred afterwards, so in this research we will only process and analyze the 568 SLR stereo image pairs acquired before the final failure event (12:25:00 to 14:27:20, hereafter the pre-failure stage) and 95 HSC stereo image pairs recorded during the failure process (14:27:22:000 to 14:27:26:700, hereafter the failure stage).

Image Processing for Stereo Image Series Recording Landslide Simulation Experiment
The low-quality landslide images were first preprocessed by using image enhancement techniques, such as a Wallis filter [77], in the landslide region to enhance and sharpen the texture patterns and increase the signal-to-noise ratio.The proposed image matching approach was then implemented on the stereo SLR and HSC image series to produce matched points.Figure 12 gives an example of the area-based image matching result for a pair of HSC stereo images.The area-based image-matching method was employed to the simulated landslide stereo image series, and the resulting number of matched point pairs for both SLR and HSC images is shown in Figure 13, from which we can see that there are apparently more matched point pairs from SLR image series (around 32,000) than those from HSC images (around 15,000).The reason of this difference could be summarized as follows: (1) the texture in color images (SLR) are usually easier to be processed and recognized than that in grayscale images (HSC); (2) the larger focal length of SLR cameras (35.0 mm vs. 20.0mm, 75% larger) can help to capture more details in the landslide surface compared with HSC cameras, although the distance to the landslide platform is a bit longer for SLR cameras than that of HSC cameras; and (3) more importantly, the images taken by SLR cameras are mostly before the landslide failure event, with more distinct texture and less influence from debris flow, while the images captured by HSC cameras in landslide failure process are more inclined to be jeopardized by the debris flow and water reflection.This is also the reason that the number of The area-based image-matching method was employed to the simulated landslide stereo image series, and the resulting number of matched point pairs for both SLR and HSC images is shown in Figure 13, from which we can see that there are apparently more matched point pairs from SLR image series (around 32,000) than those from HSC images (around 15,000).The reason of this difference could be summarized as follows: (1) the texture in color images (SLR) are usually easier to be processed and recognized than that in grayscale images (HSC); (2) the larger focal length of SLR cameras (35.0 mm vs. 20.0mm, 75% larger) can help to capture more details in the landslide surface compared with HSC cameras, although the distance to the landslide platform is a bit longer for SLR cameras than that of HSC cameras; and (3) more importantly, the images taken by SLR cameras are mostly before the landslide failure event, with more distinct texture and less influence from debris flow, while the images captured by HSC cameras in landslide failure process are more inclined to be jeopardized by the debris flow and water reflection.This is also the reason that the number of matched points from SLR image series goes down towards the occurrence of the landslide failure.Nonetheless, the almost evenly-distributed matched point pairs generated by the proposed approach meet the requirement of landslide surface monitoring in the simulated experiment.

Evaluation of Image Matching Results
To assess the reliability of the proposed matching method on the stereo image series, the image matching results are evaluated both qualitatively (distribution and NCC of the matches, see Figure 14) and quantitatively (manual checking of the matching points).Figure 14 gives an example of the evaluation of the image matching result, it can be observed that most of the landslide body surface is covered by the densely-matched point pairs for the SLR and HSC image, respectively.The total number of matched point pairs for this SLR stereo images is 32,821, more than the number of feature points extracted in the right image due to the area-based matching step.It is noteworthy that there are some small regions or spots that are sparsely matched or non-matched in both SLR and HSC images, mainly caused by the shadows that produced by different illumination conditions, influence from cables of the underground sensors that are used to

Evaluation of Image Matching Results
To assess the reliability of the proposed matching method on the stereo image series, the image matching results are evaluated both qualitatively (distribution and NCC of the matches, see Figure 14) and quantitatively (manual checking of the matching points).Figure 14 gives an example of the evaluation of the image matching result, it can be observed that most of the landslide body surface is covered by the densely-matched point pairs for the SLR and HSC image, respectively.The total number of matched point pairs for this SLR stereo images is 32,821, more than the number of feature points extracted in the right image due to the area-based matching step.It is noteworthy that there are some small regions or spots that are sparsely matched or non-matched in both SLR and HSC images, mainly caused by the shadows that produced by different illumination conditions, influence from cables of the underground sensors that are used to collect other information related to landslide event [27,64], and occlusion of a borehole pressure gauge (the blue pipe in Figure 11).To further evaluate the image matching result, NCC is calculated for each pair of matched points with a 17-pixel by 17-pixel window size.The distribution and statistic of NCC are shown in Figure 14.We can see that the NCC ranges from 0.65 to 1.0, within the thresholds that were set in the feature-based and area-based matching processes.The NCC histograms of the SLR and HSC images present similar distribution, indicating the reliability of the area-based matching result.
To evaluate the image matching reliability quantitatively, 10% evenly-distributed matched point pairs were randomly selected as samples from all the matched points, and visual checking was applied to examine the correctness of each matched point pair with a threshold of two pixels for the corresponding points.For the exampled SLR image pair, 3282 (10% of 32,821) matched point pairs were checked and 51 mismatches were found, showing a correctness of 98.45%; for the HSC image pair, 1743 (10% of 17,426) matched point pairs were checked and 38 mismatches were countered, a correctness of 97.82%.This evaluation result shows the reliability of the image matching result in the experiment.

Evaluation of Image Matching Results
To assess the reliability of the proposed matching method on the stereo image series, the image matching results are evaluated both qualitatively (distribution and NCC of the matches, see Figure 14) and quantitatively (manual checking of the matching points).Figure 14 gives an example of the evaluation of the image matching result, it can be observed that most of the landslide body surface is covered by the densely-matched point pairs for the SLR and HSC image, respectively.The total number of matched point pairs for this SLR stereo images is 32,821, more than the number of feature points extracted in the right image due to the area-based matching step.It is noteworthy that there are some small regions or spots that are sparsely matched or non-matched in both SLR and HSC images, mainly caused by the shadows that produced by different illumination conditions, influence from cables of the underground sensors that are used to collect other information related to landslide event [27,64], and occlusion of a borehole pressure

Results of Landslide Surface and Volume Changes
This part will mainly analyze the changes of the landslide surface from two aspects, including landslide volume changes and surface elevation changes.

DSM Generation
On the basis of the matched point pairs and the corresponding image EOs and IOs, a forward intersection method [71] is used to calculate the space coordinates, followed by denoising [78] and natural neighbor interpolation [79] of the point cloud to generate landslide DSM series.Figure 15 gives examples of the generated DSMs from both SLR stereo images and HSC stereo images, respectively.Especially, a video clip based on the DSMs from HSC stereo images is generated to depict the simulated landslide surface changes during the approximate five-second slope failure process (see the Supplementary File).
Remote Sens. 2016, 8, 396 14 of 24 gauge (the blue pipe in Figure 11).To further evaluate the image matching result, NCC is calculated for each pair of matched points with a 17-pixel by 17-pixel window size.The distribution and statistic of NCC are shown in Figure 14.We can see that the NCC ranges from 0.65 to 1.0, within the thresholds that were set in the feature-based and area-based matching processes.The NCC histograms of the SLR and HSC images present similar distribution, indicating the reliability of the area-based matching result.
To evaluate the image matching reliability quantitatively, 10% evenly-distributed matched point pairs were randomly selected as samples from all the matched points, and visual checking was applied to examine the correctness of each matched point pair with a threshold of two pixels for the corresponding points.For the exampled SLR image pair, 3282 (10% of 32,821) matched point pairs were checked and 51 mismatches were found, showing a correctness of 98.45%; for the HSC image pair, 1743 (10% of 17,426) matched point pairs were checked and 38 mismatches were countered, a correctness of 97.82%.This evaluation result shows the reliability of the image matching result in the experiment.

Results of Landslide Surface and Volume Changes
This part will mainly analyze the changes of the landslide surface from two aspects, including landslide volume changes and surface elevation changes.

DSM Generation
On the basis of the matched point pairs and the corresponding image EOs and IOs, a forward intersection method [71] is used to calculate the space coordinates, followed by denoising [78] and natural neighbor interpolation [79] of the point cloud to generate landslide DSM series.Figure 15 gives examples of the generated DSMs from both SLR stereo images and HSC stereo images, respectively.Especially, a video clip based on the DSMs from HSC stereo images is generated to depict the simulated landslide surface changes during the approximate five-second slope failure process (see the Supplementary File).

Landslide Volume Changes
To analyze the volume changes at different regions caused by rainfall-induced landslide experiment, the landslide simulation platform was divided into three sections based on the slopes, sections 1-3.The landslide volume changes at each section was obtained by subtracting the initial DSMs at 12:25:00 and 14:27:22.000for SLR cameras and HSC cameras, respectively, shown in Figure 16.

Landslide Volume Changes
To analyze the volume changes at different regions caused by rainfall-induced landslide experiment, the landslide simulation platform was divided into three sections based on the slopes, Sections 1-3.The landslide volume changes at each section was obtained by subtracting the initial DSMs at 12:25:00 and 14:27:22.000for SLR cameras and HSC cameras, respectively, shown in Figure 16.The volume changes at each section before landslide failure are illustrated in Figure 16a.Due to power issues, there are two data loss periods (13:14:30-13:31:50 and13:52:40-14:02:50), yet this does not affect the analysis of the landslide change tendency.It can be observed that, in this stage, the landslide mass changed slightly: section 2 is reducing, section 3 increasing, and barely changes in section 1, indicating the major changes at the middle and toe parts of the landslide body during 12:25:00 to 14:27:20.Although no collapse occurs, this process has been accompanied by an increase of landslide sliding energy.With the persistent rainfall and landslide mass sliding, the energy reached its limit and triggered the collapse event which was recorded by the HSC system from 14:27:22.000 to 14:27:26.700,and distinct changes in all the three sections were presented in Figure 16b.In contrast to the pre-failure stage, in this stage the landslide mass in section 1 reduced rapidly, in section 3 increased significantly, and in section 2 endured a slowly increase manner.
The landslide surface volume at each section in the two stages was also estimated by setting the reference level at the bottom of the platform (2.5 m).Table 3 shows the landslide surface volume and volume difference recorded by SLR and HSC cameras.Surface volume changes at each section varied throughout different stages.Slight changes occurred in thepre-failure process (12:25:00-14:27:20) for all the three sections.Section 2, serving as the initialization zone that provided input to section 3, lost about 0.30 m 3 in volume.Significant volume changes substantially occurred in the slope failure process (14:27:22.000-14:27:26.700).Section 1 lost about 0.84 m 3 landslide mass, supplementing sections 2 and 3.In addition, the volume of section 3 increased 0.80 m 3 , indicating the major deposition zone.A noticeable feature of the surface volume is that there is an extraordinarily huge difference for the surface volumes of section 1 between 14:27:20 and 14:27:22.000,recorded by SLR and HSC cameras, respectively.This difference is mainly caused by the different Vertical Field Views (VFVs) of the two camera systems.As shown in Figure 17, the VFV of SLR cameras is smaller (18°) than that of HSC cameras (26°), resulting in a smaller imaging region and leading to significantly smaller surface volume in section 1 compared with the HSC cameras.The volume changes at each section before landslide failure are illustrated in Figure 16a.Due to power issues, there are two data loss periods (13:14:30-13:31:50 and13:52:40-14:02:50), yet this does not affect the analysis of the landslide change tendency.It can be observed that, in this stage, the landslide mass changed slightly: Section 2 is reducing, Section 3 increasing, and barely changes in Section 1, indicating the major changes at the middle and toe parts of the landslide body during 12:25:00 to 14:27:20.Although no collapse occurs, this process has been accompanied by an increase of landslide sliding energy.With the persistent rainfall and landslide mass sliding, the energy reached its limit and triggered the collapse event which was recorded by the HSC system from 14:27:22.000 to 14:27:26.700,and distinct changes in all the three sections were presented in Figure 16b.In contrast to the pre-failure stage, in this stage the landslide mass in Section 1 reduced rapidly, in Section 3 increased significantly, and in Section 2 endured a slowly increase manner.
The landslide surface volume at each section in the two stages was also estimated by setting the reference level at the bottom of the platform (2.5 m).Table 3 shows the landslide surface volume and volume difference recorded by SLR and HSC cameras.Surface volume changes at each section varied throughout different stages.Slight changes occurred in thepre-failure process (12:25:00-14:27:20) for all the three sections.Section 2, serving as the initialization zone that provided input to Section 3, lost about 0.30 m 3 in volume.Significant volume changes substantially occurred in the slope failure process (14:27:22.000-14:27:26.700).Section 1 lost about 0.84 m 3 landslide mass, supplementing Sections 2 and 3.In addition, the volume of Section 3 increased 0.80 m 3 , indicating the major deposition zone.A noticeable feature of the surface volume is that there is an extraordinarily huge difference for the surface volumes of section 1 between 14:27:20 and 14:27:22.000,recorded by SLR and HSC cameras, respectively.This difference is mainly caused by the different Vertical Field Views (VFVs) of the two camera systems.As shown in Figure 17, the VFV of SLR cameras is smaller (18 ˝) than that of HSC cameras (26 ˝), resulting in a smaller imaging region and leading to significantly smaller surface volume in section 1 compared with the HSC cameras.Although the change tendency of the landslide mass agrees with the actual process in the two stages, the calculated overall surface volumes of the landslide mass are not consistent during the experiment.Table 3 gives the calculated surface volume differences, and it can be seen that in the pre-failure stage the volume change was negative (−0.18 m 3 ), while in the failure stage it was positive (0.17 m 3 ).This phenomenon was also caused by the difference of the VFVs of the two camera systems, as illustrated in Figure 17.In the pre-failure stage, the view of the SLR cameras was obstructed by the waterway at the landslide toe, resulting in the blind spot in which the shaded region (Figure 17a) cannot be included to calculate the surface volume; therefore, the surface volume change was negative.In the failure stage, the reason for superfluous surface volume was that the shaded region in Figure 17b cannot be included into the surface volume calculation due to the VFV of the HSC system at the beginning (14:27:22.000)but, in the end (14:27:26:700), this shaded region was added to the surface volume calculation.Nonetheless, the small volume differences induced from VFVs are not significant to the surface volume comparison in the experiment.

Landslide Surface Elevation Changes
This section will describe the landslide surface elevation evolution during the sliding process by comparison of the DSMs at different times.Since there were no obvious changes in the pre-failure stage, here we only select two moments to demonstrate, and will pay more attention to the surface elevation changes in the failure stage.Figure 18 shows the elevation difference series of DSMs generated by ArcGIS 10.0 software at 12:58:10 and 13:14:00 compared to DSM at 12:25:00 in the pre-failure stage, and DSMs at 14:27:22.950(II), 14:27:23.950(III), 14:27:24.950(IV), 14:27:25.950(V) and 14:27:26.700(VI) compared to DSM at 14:27:22.000(I) in the failure stage.The results showed that in the pre-failure stage, the sliding surface was partly outcropping at the foot of the slope with slight landslide mass sliding.Major and rapid sliding processes and sudden surface collapse happened in the failure stage.The landslide mass in section 2 first slid into section 3, and then a large amount of the landslide mass slid down the slope; meanwhile, the sliding area continued to move upward to section 1.This process continued until the landside failure event stopped.Although the change tendency of the landslide mass agrees with the actual process in the two stages, the calculated overall surface volumes of the landslide mass are not consistent during the experiment.Table 3 gives the calculated surface volume differences, and it can be seen that in the pre-failure stage the volume change was negative (´0.18 m 3 ), while in the failure stage it was positive (0.17 m 3 ).This phenomenon was also caused by the difference of the VFVs of the two camera systems, as illustrated in Figure 17.In the pre-failure stage, the view of the SLR cameras was obstructed by the waterway at the landslide toe, resulting in the blind spot in which the shaded region (Figure 17a) cannot be included to calculate the surface volume; therefore, the surface volume change was negative.In the failure stage, the reason for superfluous surface volume was that the shaded region in Figure 17b cannot be included into the surface volume calculation due to the VFV of the HSC system at the beginning (14:27:22.000)but, in the end (14:27:26:700), this shaded region was added to the surface volume calculation.Nonetheless, the small volume differences induced from VFVs are not significant to the surface volume comparison in the experiment.

Landslide Surface Elevation Changes
This section will describe the landslide surface elevation evolution during the sliding process by comparison of the DSMs at different times.Since there were no obvious changes in the pre-failure stage, here we only select two moments to demonstrate, and will pay more attention to the surface elevation changes in the failure stage.Figure 18 shows the elevation difference series of DSMs generated by ArcGIS 10.0 software at 12:58:10 and 13:14:00 compared to DSM at 12:25:00 in the pre-failure stage, and DSMs at 14:27:22.950(II), 14:27:23.950(III), 14:27:24.950(IV), 14:27:25.950(V) and 14:27:26.700(VI) compared to DSM at 14:27:22.000(I) in the failure stage.The results showed that in the pre-failure stage, the sliding surface was partly outcropping at the foot of the slope with slight landslide mass sliding.Major and rapid sliding processes and sudden surface collapse happened in the failure stage.The landslide mass in Section 2 first slid into Section 3, and then a large amount of the landslide mass slid down the slope; meanwhile, the sliding area continued to move upward to Section 1.This process continued until the landside failure event stopped.To further investigate the surface elevation changes in the failure stage, the elevation profiles of the center line of the selected DSMs at the six key moments (same as Figure 18) are shown in Figure 19.All profiles show clearly that the boundary of sliding area and deposition area was 3.5 m far away from the front edge of platform.A loss of material up to 0.5 m in thickness occurred in the sliding area (upper part of section 1), and a large amount of material (up to 0.3 m in thickness) was accumulated downward the slope.In the end, the landslide toe crossed the barrier of the waterway and rushed out of the platform.To further investigate the surface elevation changes in the failure stage, the elevation profiles of the center line of the selected DSMs at the six key moments (same as Figure 18) are shown in Figure 19.All profiles show clearly that the boundary of sliding area and deposition area was 3.5 m far away from the front edge of platform.A loss of material up to 0.5 m in thickness occurred in the sliding area (upper part of Section 1), and a large amount of material (up to 0.3 m in thickness) was accumulated downward the slope.In the end, the landslide toe crossed the barrier of the waterway and rushed out of the platform.
the center line of the selected DSMs at the six key moments (same as Figure 18) are shown in Figure 19.All profiles show clearly that the boundary of sliding area and deposition area was 3.5 m far away from the front edge of platform.A loss of material up to 0.5 m in thickness occurred in the sliding area (upper part of section 1), and a large amount of material (up to 0.3 m in thickness) was accumulated downward the slope.In the end, the landslide toe crossed the barrier of the waterway and rushed out of the platform.

Comparison with Result of SGM Algorithm
In the landslide simulation experiment, the constantly changing landslide surface makes it very difficult to obtain transient ground truth of the landslide surface with high-precision using a laser scanner or other structured light sensor, such as the Microsoft Kinect [80].Commonly used ground laser scanners usually generate noise when scanning moving objects, especially for a non-rigid surface, such as that in the landslide failure process.In contrast, the Microsoft Kinect can record the sliding landslide surface in 3D with high frequency, but the poor positioning accuracy resulted from the low image resolution and the short ranging limit constrains its application for ground truth in the setup of the landslide simulation platform.Thus, we compare the image-matching result in this research with that from SGM, the commonly used algorithm in stereo computer vision that was employed to match glacier images with poor texture similar to the landslide surface [81].SGM performs pixel-wise matching based on a certain matching cost (such as mutual information or census) [73] and the approximation of a global smoothness constraint [52].The comparison is focused on the number and distribution of the matches, and the reliability of the matching result.
The SGM code employed in this research is from the LibTSgm library developed by the Institute for Photogrammetry at the University of Stuttgart [82], and can be downloaded freely from their website [83].Practical experience has shown that SGM is very robust in different applications and does not require parameter tuning [81], so here the SGM code is fully automated.The SGM algorithm was applied to the exampled SLR and HSC stereo image pairs and the image matching results were obtained.The numbers of matches from the SLR images and HSC images were 1,006,691 and 404,005, respectively, vastly outnumbering those from the proposed method due to the nature of per-pixel matching of SGM.The SGM matches are distributed evenly on the landslide surface of the SLR image, while there was a blank region on the bottom left side of the HSC image, as shown in Figure 20, mainly caused by the reflection of water areas.
Next, a similar matching reliability evaluation described in Section 4.2 was performed for the SGM results.Here we randomly selected 2000 evenly-distributed matched points from SLR and HSC stereo image pairs for visual inspection using the same threshold of two pixels.Results showed that the correctness of the matching result for SLR and HSC images was 75.70% (486 mismatches) and 70.50% (590 mismatches), respectively.Compared with the correctness rates evaluated in Section 4.2, we found that the proposed image matching algorithm significantly outweighed SGM method for both SLR (98.45% vs. 75.70%)and HSC (97.82% vs. 70.50%)stereo image pairs.
From the above analysis, we can see that although the SGM method generates matches much more as a result of per-pixel dense matching, our proposed method is more robust on the overall distribution of matched points and the reliability of the matching result for both SLR and HSC images in the landslide simulation experiment.
for Photogrammetry at the University of Stuttgart [82], and can be downloaded freely from their website [83].Practical experience has shown that SGM is very robust in different applications and does not require parameter tuning [81], so here the SGM code is fully automated.The SGM algorithm was applied to the exampled SLR and HSC stereo image pairs and the image matching results were obtained.The numbers of matches from the SLR images and HSC images were 1,006,691 and 404,005, respectively, vastly outnumbering those from the proposed method due to the nature of per-pixel matching of SGM.The SGM matches are distributed evenly on the landslide surface of the SLR image, while there was a blank region on the bottom left side of the HSC image, as shown in Figure 20, mainly caused by the reflection of water areas.Next, a similar matching reliability evaluation described in Section 4.2 was performed for the SGM results.Here we randomly selected 2000 evenly-distributed matched points from SLR and HSC stereo image pairs for visual inspection using the same threshold of two pixels.Results showed that the correctness of the matching result for SLR and HSC images was 75.70% (486 mismatches) and 70.50% (590 mismatches), respectively.Compared with the correctness rates evaluated in Section 4.2, we found that the proposed image matching algorithm significantly outweighed SGM method for both SLR (98.45% vs. 75.70%)and HSC (97.82% vs. 70.50%)stereo image pairs.
From the above analysis, we can see that although the SGM method generates matches much more as a result of per-pixel dense matching, our proposed method is more robust on the overall distribution of matched points and the reliability of the matching result for both SLR and HSC images in the landslide simulation experiment.

Accuracy Evaluation
The geometric accuracy of the detected feature points is very important for the reliable estimation of the landslide volume and surface elevation changes in the simulation experiment.Taking the HSC system as an example, the average ground position accuracy of the surface features is estimated in the following.
The coordinates of GCPs were measured by a Sokkia total station with an accuracy of ˘1 mm.The estimated accuracy of the EO parameters given by PhotoModeler software for the HSC system includes: (a) 0.3 mm in the X (horizontal) and Z (vertical) directions, and 0.4 mm in the Y (depth) direction for the camera center; and (b) 0.006 ˝, 0.021 ˝, and 0.018 ˝for the orientation angles about the X, Z, and Y coordinate axes, respectively.The calculated RMSEs of the ground coordinates were 0.9 mm in the X and Z directions, and 1.9 mm in the Y direction.Thus, the ground position accuracy was 2.3 mm.This estimated accuracy is considered as an internal accuracy that is caused by the EO parameter errors.
The detected surface features were generally non-structured image features.They were identified and matched at an estimated accuracy of 0.25 pixels.Given the camera baseline of 1.16 m for HSC system, and a depth of 5 m from the baseline to the middle of the slope, the computed ground position accuracy was 0.7 mm in the X and Z directions and 2.8 mm in the Y direction.Hence, the position error of surface features in the middle of the slope caused by the image feature measurement errors was 3.0 mm.Similarly, we can estimate that this error ranges from 1.1 mm for the closest surface features to 5.6 mm for the farthest features on the slope.Sicne the surface features changed constantly and were not accessible during the experiment, no ground truth was available to verify the accuracy of these moving surface features.Overall, the average ground position accuracy of the surface features was approximately 3.8 mm, as estimated using the error propagation law from the above-discussed position error components caused by errors of the EO parameters (2.3 mm), and image coordinate measurement (3.0 mm) under the assumption that these errors are independent.
In a similar manner, the average ground position accuracy of the surface features detected by the SLR camera system was approximately 3.2 mm.Considering the dimension of the landslide simulation platform (6 m ˆ1.5 m ˆ3 m), the estimated ground position accuracy of 3.2 mm and 3.8 mm is sufficient for reliable estimation of landslide volume and elevation changes.

Conclusions
This paper presents a novel robust image-matching approach for poor-texture, close-range images, composed of two steps: multiple-constraints assisted feature-based image matching and area-based image matching.The SLR and HSC stereo image series in the simulated landslide experiment were employed to illustrate the proposed method, and the reliable image matching results were obtained.Then the corresponding DSM series during the landslide process were constructed using the close-range photogrammetry technique, followed by the discussion of the landslide volume changes and surface elevation changes in the simulation experiment.
The research results support the following conclusions: (1) The proposed multiple-constraints based feature-to area-image matching methodology is capable of robustly matching the close-range, poor-texture images, obtaining almost evenly-and densely-distributed matches with sufficient matching accuracy.(2) The matching result of this method is relevant to the image quality that is usually affected by both the camera and capture settings, such as the image resolution and surface reflex of the object.For example, in the simulated landslide experiment, more matched points could be obtained from the color SLR images than from the grayscale HSC images due to the better imaging condition (e.g., higher resolution, less influence by water-pond regions).( 3) The proposed robust image-matching method can be successfully applied to the low-frequency SLR and high-frequency HSC stereo image series collected in the simulated landslide experiment for generation of sequential DSMs, which helps to reveal the landslide evolution process triggered by rainfall, especially based on the volume and surface elevation changes in the instantaneous failure event.
Despite the achievements in this research, there are currently several limitations that need further improvement in the future.For example, there are still very sparse, or even non-matched, regions affected by different factors, such as sensor cables and sliding fissures (see Figures 12 and 16).Due to the main objective in this research to observe the landslide surface changes, the accuracy of the generated DSMs is not discussed in detail in this paper.These limitations, which would be involved in imaged landslide surface deformation analysis, call for further comprehensive research in the future.

Figure 1 .
Figure 1.Examples of stereo images with poor texture in landslide simulation.Note: the red boxes represent the Ground Control Point (GCP) marks.(a) left image of the SLR stereo pair; (b) right image of SLR stereo pair; (c) left image of the HSC stereo pair; and (d) right image of the HSC stereo pair.

Figure 1 .
Figure 1.Examples of stereo images with poor texture in landslide simulation.Note: the red boxes represent the Ground Control Point (GCP) marks.(a) left image of the SLR stereo pair; (b) right image of SLR stereo pair; (c) left image of the HSC stereo pair; and (d) right image of the HSC stereo pair.

Figure 2 .
Figure 2. Landslide simulation platform used for collection of HSC and SLR stereo images.Note: (a) installation of the stereo camera systems on the landslide simulation platform; and (b) side view of the landslide simulation platform geometry.

Figure 2 .
Figure 2. Landslide simulation platform used for collection of HSC and SLR stereo images.Note: (a) installation of the stereo camera systems on the landslide simulation platform; and (b) side view of the landslide simulation platform geometry.

Figure 3 .
Figure 3. Flowchart of the proposed poor-texture close-range image processing approach.

Figure 3 .
Figure 3. Flowchart of the proposed poor-texture close-range image processing approach.

Figure 4 .
Figure 4. Feature points extracted by the SIFT operator in the SLR stereo image pair .Note: each point is a red dot, and the total numbers are 48,735 (a) and 21,548 (b), respectively.

Figure 4 .
Figure 4. Feature points extracted by the SIFT operator in the SLR stereo image pair .Note: each point is a red dot, and the total numbers are 48,735 (a) and 21,548 (b), respectively.

Figure 6 .
Figure 6.Filtering result of the first-level matched features.Note: blue features are filtered by FM, yellow ones by NCC, red ones are the real matched features.(a) left image; and (b) right image.

Figure 6 .
Figure 6.Filtering result of the first-level matched features.Note: blue features are filtered by FM, yellow ones by NCC, red ones are the real matched features.(a) left image; and (b) right image.

Figure 7 .Figure 8 .
Figure 7. Triangulation of the first-level matched features.Note: the red points are the refined matched features, and the blue lines are the triangulate edges.For visual effects, here we only plot the feature points within the yellow triangle.(a) left image; and (b) right image.

Figure 7 .
Figure 7. Triangulation of the first-level matched features.Note: the red points are the refined matched features, and the blue lines are the triangulate edges.For visual effects, here we only plot the feature points within the yellow triangle.(a) left image; and (b) right image.

Figure 7 .Figure 8 .
Figure 7. Triangulation of the first-level matched features.Note: the red points are the refined matched features, and the blue lines are the triangulate edges.For visual effects, here we only plot the feature points within the yellow triangle.(a) left image; and (b) right image.

Figure 8 .
Figure 8.The result of the feature-based matching.Note: the total number of matched features is 10,778.(a) left image; and (b) right image.

Figure 10 .
Figure 10.Triangulation generation and matching region prediction of area-based matching process.Note: (a) master image; and (b) searching image; A and B are two corresponding triangles, and C is the searching window in B for corresponding point of Pm.

Figure 10 .
Figure 10.Triangulation generation and matching region prediction of area-based matching process.Note: (a) master image; and (b) searching image; A and B are two corresponding triangles, and C is the searching window in B for corresponding point of Pm.

Figure 11 .
Figure 11.Area-based image matching result for SLR stereo image pair.The total number of matched points is 32,821.Note: (a) left image; and (b) right image.

Figure 11 .
Figure 11.Area-based image matching result for SLR stereo image pair.The total number of matched points is 32,821.Note: (a) left image; and (b) right image.

Figure 12 .
Figure 12.Example of area-based image matching for a pair of HSC stereo images.Note: (a) original HSC image pair; (b) SIFT features (23,250 points in the left image and 23,583 points in the right image); (c) feature-based image-matching result (4999 matched feature points); and (d) area-based image-matching result (17,426 matched feature points).

Figure 12 .
Figure 12.Example of area-based image matching for a pair of HSC stereo images.Note: (a) original HSC image pair; (b) SIFT features (23,250 points in the left image and 23,583 points in the right image); (c) feature-based image-matching result (4999 matched feature points); and (d) area-based image-matching result (17,426 matched feature points).

Figure 13 .
Figure 13.Number of matched point pairs for SLR stereo images (pre-failure stage) and HSC stereo images (failure stage).

Figure 14 .
Figure 14.Distribution and statistics of NCC of the image matching result for SLR and HSC stereo image pairs.Note: (a,b) are NCC distribution for SLR and HSC images, respectively, and (c,d) are the corresponding NCC statistics.

Figure 13 .
Figure 13.Number of matched point pairs for SLR stereo images (pre-failure stage) and HSC stereo images (failure stage).

Figure 13 .
Figure 13.Number of matched point pairs for SLR stereo images (pre-failure stage) and HSC stereo images (failure stage).

Figure 14 .
Figure 14.Distribution and statistics of NCC of the image matching result for SLR and HSC stereo image pairs.Note: (a,b) are NCC distribution for SLR and HSC images, respectively, and (c,d) are the corresponding NCC statistics.

Figure 14 .
Figure 14.Distribution and statistics of NCC of the image matching result for SLR and HSC stereo image pairs.Note: (a,b) are NCC distribution for SLR and HSC images, respectively, and (c,d) are the corresponding NCC statistics.

Figure 16 .
Figure 16.Landslide volume changes at each Section during the simulated experiment.Note: (a) pre-failure process recorded by SLR stereo images; (b) failure event captured by HSC stereo images.

Figure 16 .
Figure 16.Landslide volume changes at each Section during the simulated experiment.Note: (a) pre-failure process recorded by SLR stereo images; (b) failure event captured by HSC stereo images.

Figure 17 .
Figure 17.The geometry of VFV for SLR and HSC systems at the key moments (onset time, start, and end of failure process).(a) Zoomed geometry of VFV for SLR camerras in pre-failure stage; (b) Zoomed geometry of VFV for HSCS in failure stage.

Figure 17 .
Figure 17.The geometry of VFV for SLR and HSC systems at the key moments (onset time, start, and end of failure process).(a) Zoomed geometry of VFV for SLR camerras in pre-failure stage; (b) Zoomed geometry of VFV for HSCS in failure stage.

Figure 18 .
Figure 18.Landslide surface elevation changes at some critical moments in the pre-failure and failure stages.Note: the first row is the pre-failure stage, the second and third rows are the failure stage.The figures of surface elevation difference are generated via 3D Analyst Tools of ArcGIS 10.0 software using two DSMs.

Figure 18 .
Figure 18.Landslide surface elevation changes at some critical moments in the pre-failure and failure stages.Note: the first row is the pre-failure stage, the second and third rows are the failure stage.The figures of surface elevation difference are generated via 3D Analyst Tools of ArcGIS 10.0 software using two DSMs.

14: 27 :Figure 19 .
Figure 19.Landslide elevation profiles at six key moments in the failure stage.Figure 19.Landslide elevation profiles at six key moments in the failure stage.

Figure 19 .
Figure 19.Landslide elevation profiles at six key moments in the failure stage.Figure 19.Landslide elevation profiles at six key moments in the failure stage.

Figure 20 .
Figure 20.SGM result for the selected HSC stereo image pair (404,005 matches).Note: (a) left image, and (b) right image.

Figure 20 .
Figure 20.SGM result for the selected HSC stereo image pair (404,005 matches).Note: (a) left image, and (b) right image.

Table 1 .
Parameters of the stereo camera systems in the landslide simulation experiment.

Table 1 .
Parameters of the stereo camera systems in the landslide simulation experiment.

Table 2 .
Comparison of detected feature points using different feature detection methods.Note: the HSC stereo image pair here is the same as that in Section 4.1.

Table 2 .
Comparison of detected feature points using different feature detection methods.Note: the HSC stereo image pair here is the same as that in Section 4.1.

Table 3 .
Landslide surface volume and volume difference recorded by the SLR and HSC cameras.

Table 3 .
Landslide surface volume and volume difference recorded by the SLR and HSC cameras.