Open AccessArticle

Terrain Point Cloud Assisted GB-InSAR Slope and Pavement Deformation Differentiate Method in an Open-Pit Mine

Xiangtian Zheng

Xiufeng He

^1,*

Xiaolin Yang

²,

Haitao Ma

²,

Zhengxing Yu

²,

Guiwen Ren

²,

Jiang Li

²,

Hao Zhang

² and

Jinsong Zhang

School of Earth Science and Engineering, Hohai University, Nanjing 211100, China

China Academy of Safety Science and Technology, Beijing 100012, China

Beijing Institute of Surveying and Mapping, Beijing 100038, China

Author to whom correspondence should be addressed.

Sensors 2020, 20(8), 2337; https://doi.org/10.3390/s20082337

Submission received: 6 March 2020 / Revised: 9 April 2020 / Accepted: 14 April 2020 / Published: 20 April 2020

(This article belongs to the Section Remote Sensors)

Download

Browse Figures

Figure 1
(a) Optical image of the target open pit slope. (b) It is difficult to distinguish the pavement surface from the slope surface using the mapping result generated by nearest-neighbor interpolation. The orange color areas consist of several pixels across multiple slope steps. A ground-based synthetic aperture radar interferometry (GB-InSAR) cumulative deformation map was produced in a month, and was mapped on the terrain model. The orange-colored area shows great deformation but is difficult to interpret by geotechnical engineers. "> Figure 2
Logical experiment scheme of the data fusion between GB-InSAR interferograms and terrain point cloud data using nearest-neighbor interpolation and geometric-scattering model separately. "> Figure 3
The synthetic aperture model was observed at a fixed-station, overhead perspective. (a) The key process of GB-InSAR echo data processing. (b) S-SAR radar system which contains transmitting and receiving antennas, radar sensor and a linear rail. (c) Ideal state: "One-step, one-stop" SFCW signal transmitting and receiving mode of straight-line repeated-track controlled by a pulse-width modulation (PWM). "> Figure 4
(a) The GB-InSAR image grid. The spatial resolution of the grid is lower than terrestrial laser scanning (TLS) points cloud resolution. Each grid may correspond to multiple terrain surface model points. (b) GB-InSAR image pixel responds to the same slope area after 2D image co-registration. "> Figure 5
(a) The geometric mapping relationship between normal center rectangular coordinate system GB-InSAR and 3D terrain data. (b) Terrain point cloud azimuth angle top view. "> Figure 6
Terrain point cloud introduces greater detail so that radar wave incident angles define geometric and scattering models. The model helps in solving the one-to-many problem during the data fusion process. For the road surface and slope surface of the slope, they have completely different relative incident angles. "> Figure 7
Assignment weights according to point cloud incident angle geometric-scattering model. "> Figure 8
(a) Surface color of the road-slope model indicating <math display="inline"><semantics> <mrow> <msub> <mi>θ</mi> <mrow> <mi>s</mi> <mi>l</mi> <mi>o</mi> <mi>p</mi> <mi>e</mi> </mrow> </msub> </mrow> </semantics></math> variance extracted by Equation (5). (b) Geocoded model with the actual collected terrain point cloud to a local coordinate system via linear transformation. (c) Normal vectors of the pavement and slope of the rectangle. (d) <math display="inline"><semantics> <mrow> <mi>γ</mi> <mo stretchy="false">(</mo> <msub> <mi>θ</mi> <mi>i</mi> </msub> <mo stretchy="false">)</mo> </mrow> </semantics></math> the variance of each point. The deep blue color area is the pavement. In this case, they have a bad incident angle. "> Figure 9
(a) Simulated sub-areas with small deformation 2D deformation map. (b) Geometrical mapping 3D result of deformation in (a) using a nearest-neighbor interpolation method. "> Figure 10
(a) 2D deformation diagram shows a slight caving-in phenomenon. (b) The nearest-neighbor method creates a 3D result which has no obvious discrimination between slope and road surfaces. (c) Simulated deformation diagram of Region D after accumulation. (d) Accumulated 3D matching produces obvious banding. A similar phenomenon occurs in adjacent regions to D. (e) Nearest-neighbor method result. (f) The result generated by the geometric and scattering weight model, where deformation results distinguish road surface from the slope surface. The pavement deformation disappears because the incident angle is not good. If the pavement surface is rugged as shown in <a href="#sensors-20-02337-f006" class="html-fig">Figure 6</a>, the pavement deformation would be reserved. "> Figure 11
(a) Radar local incident angle and (b) normalized results, local zoom in <a href="#sensors-20-02337-f011" class="html-fig">Figure 11</a>a. The arrow represents the direction of the incident vector (Equation (11)). (c) Weight varies when the scattering model changes. (d) The arrow represents the direction of local most likely deformed by gravitation or explosion fluctuation. "> Figure 12
The cumulative deformation map from GMT+8 13:03 10 October 2018 to 17:23 was used in the mapping experiment. Three sub-areas were also selected. (a) Nearest-neighbor interpolation mapping result. (b) Paper method mapping result. "> Figure 13
The deformation of the transition zone is removed to a certain extent by the model weighting method of the paper. These edge boundary anomalies, however, are often marginal areas with poor data quality, but are affected by strong coherence scattering points. (a) Nearest-neighbor interpolation mapping result of sub-area A. (b) Paper method mapping result of sub-area A’. "> Figure 14
Multiple terrain surface deformed points may locate in one radar pixel spatial resolution element which crosses pavement and slope. The local layover distortion is hard to identify without geometric and scattering information. (a) Nearest-neighbor interpolation mapping result of sub-area B. (b) Paper method mapping result of sub-area B’. "> Figure 15
A deformation region with large local undulations. (a) Nearest-neighbor interpolation mapping result of sub-area C. (b) Paper method mapping result of sub-area C’. "> Figure 16
The deformation in one radar pixel generated by two methods. (a) Nearest-neighbor interpolation deformation result. The whole pixel may equal one vital scattering target. (b) Deformation result using the geometric and scattering model. The terrain model point cloud may equal several distributed targets. ">

Versions Notes

Abstract

Ground-based synthetic aperture radar interferometry (GB-InSAR) is a valuable tool for deformation monitoring. The 2D interferograms obtained by GB-InSAR can be integrated with a 3D terrain model to visually and accurately locate deformed areas. The process has been preliminarily realized by geometric mapping assisted by terrestrial laser scanning (TLS). However, due to the line-of-sight (LOS) deformation monitoring, shadow and layover often occur in topographically rugged areas, which makes it difficult to distinguish the deformed points on the slope between the ones on the pavement. The extant resampling and interpolation method, which is designed for solving the scale difference between the point cloud and radar pixels, does not consider the local scattering characteristics difference of slope. The scattering difference information of road surface and slope surface in the terrain model is deeply weakened. We propose a differentiated method with integrated GB-InSAR and terrain surface point cloud. Local geometric and scattering characteristics of the slope were extracted, which account for pavement and slope differentiating. The geometric model is based on a GB-InSAR system with linear repeated-pass and the topographic point cloud relative observation geometry. The scattering model is based on k-nearest neighbor (KNN) points in small patches varies as radar micro-wave incident angle changes. Simulation and a field experiment were conducted in an open-pit mine. The results show that the proposed method effectively distinguishes pavement and slope surface deformation and the abnormal area boundary is partially relieved.

Keywords:

ground-based synthetic aperture radar interferometry (GB-InSAR); deformation monitoring; geometric and scattering model; slope deformation monitoring

1. Introduction

Deformation monitoring can reveal unstable slope surfaces to prevent serious safety hazards. It is deployed to determine changes in the shape and displacements of a deformable body in temporal and spatial domains [1]. Advancements in signal processing, integrated circuit, microwave, and antenna technologies have made synthetic aperture radar (SAR) a notable innovation in the radar technology field. Interferometry is a technique for extracting deformation information by comparing the observed phases before and after the deformation of an object. SAR interferometry technology can be used for all-weather situations and the high-precision monitoring of ground surface deformations across large areas [2]. Ground-based synthetic aperture radar interferometry (GB-InSAR) has been tested in various fields in recent years including the monitoring of landslides, open-pit mines, ground subsidence, architectural structures, bridges, sinkholes, and glacier displacements [3,4,5,6,7,8]. Despite these achievements, there are yet obstacles hindering the practical application of GB-InSAR deformation monitoring. Traditional SAR can obtain 2D images in steep and complex terrain, but 3D objects may be severe distorted on the 2D image which obscures several targets [9,10,11]. SAR sensors can resolve different monitored targets based on their slant ranges from the targets to the radar center and the angles away from the centerline of the radar antenna middle beam. Unlike the orthographic projection of a traditional topographic map and the central projection of photogrammetry, SAR-based image projection generally requires manual interpretation to identify a deformed target and its position in the planar coordinate system of a radar image. The interpretation of radar images requires understanding and judgment of radar observation geometry and imaging principles, matching with 3D scenes, and accurately locating deformation areas and variations. The analyst must have sufficient knowledge to ensure appropriate such identifications. The issue of accurately interpreting 3D measurement results is key regarding successful GB-InSAR implementation.

Previous research has explored GB-InSAR mapping technology in the GB-InSAR geometric mapping context. Tapet et al., for example, integrated terrestrial laser scanning (TLS) and GB-InSAR to monitor heritage building deformations and accurately located deformed areas [5]. Yue et al. fused GB-radar and survey robot data to determine the spatial relationship between the robot monitored 3D deformations and the GB-radar coordinate system [12]. The method was used for 1D deformation and was not expanded in high dimensions. GB-InSAR and InSAR have unique observation geometry and measurement modes, so applying the InSAR processing method to GB-InSAR produces a substantial error. The integration of high-precision external auxiliary measurements may be an effective solution to this problem. Zou, of Wuhan University, fused ground-based SAR and 3D laser scanning data for the 3D visualization of GB-InSAR in a hydro-dam deformation monitoring experiment using point cloud, coordinate transformation and projection, and 2D four-parameter transformation techniques [13]. Yang et al. proposed a geometric 3D matching method based on both range and azimuth constraints for frontal and lateral strip SAR and terrain data, [14]. Wang et al. used 3D laser scanning data to assist the 3D coordinate transformation of GB-InSAR image pixels; they also applied a scaling factor method based on the parameter estimation of a similar transformation. They put forward the idea of mismatching correction based on control points [9]. In general, inconsistencies among systems, variations in the measurement environment, and scene image changes remain problematic. In an experiment on ground-based SAR deformation monitoring, Yang et al. used a positioning method to determine the rail endpoint coordinate and artificial targets [4]. Lombardi et al. [15] integrated GB-InSAR and 3D laser scanning for landslide monitoring. Zheng et al. investigated data fusion and deformation comparison with GB-InSAR and 3D laser scanning in a controlled field test and later they integrated GB-InSAR, TLS and unmanned aerial vehicle photogrammetry (UAV) in a rock slide residual mass emergency monitoring case [1,16].

The geometric mapping of the GB-InSAR image and 3D terrain data is an important step in the deformation monitoring process. During a monitoring case in an open-pit mine, we encountered a problem: mapping results cannot differentiate slope steps differences. The problem confused geo-engineers in analyzing the deduction of spatial and temporal deformation patterns. Based on the above analysis, we can conclude that extant geometric mapping methods are deficient in the mapping relation of multiple spatially distributed target points of the radar pixel with the point cloud, as well as the 3D features (e.g., structure, texture, and occlusion) in reflecting targets. The points cloud may help inversion of the scattering difference of GB-InSAR, which would help differentiate the pavement and slope. Conventional mapping methods based on nearest-neighbor interpolation cannot be effectively applied to SAR pixels. Based on the basic cognition that SAR image is a kind of target scattering characteristic difference representation map. We develop a method that is comprised of a geometric mapping and scattering variance weighting model. The model is introduced in setting the interpolation weight of points in one radar pixel. The method may help in identifying shadows and layover in topographically rugged areas.

The paper is organized as follows. Section 2 discusses the problem we encountered and solving scheme by an improved linear GB-InSAR and TLS data fusion method. Section 3 further discusses the GB-InSAR system and the rough surface slope topographic data model. The GB-InSAR image and spatial point cloud 3D mapping method, which is based on geometric and scattering models, is also delineated in Section 3. Section 4 presents a data fusion experiment in an open-pit mine in China, which was conducted to validate the model. Section 5 discusses the influence factors and possible solutions. Section 6 provides a summary and recommendations for future research.

2. Linear GB-InSAR and TLS Data Fusion

2.1. Problem and Solving Scheme

A time-series of 2D images was fused with GB-InSAR and the 3D terrain surface point cloud data with TLS as shown in Figure 1. This fusion is a step-wise process. The key processes are as follow:

Calculate relative ranges and azimuth angle according to GB-InSAR monitoring geometry.
Nearest-neighbor interpolation for one radar pixel-to-multiple terrain surface points mapping.
3D visualization.

Figure 1b shows the mapped result, where the gray-scale color marks the terrain model generated by the TLS point clouds, which consists of several 1 m² patches. The colored layer is the GB-InSAR deformation 3D projection “cloud diagram”. However, the mapped result shows no differences in the pavement and slope. Geotechnical experts cannot intuitively evaluate the stability of a given slope according to the "cloud diagram" of deformation changes, and the results are shown as a band-like pattern. These bands do not conform to any constitutive model of slope deformation and would not help geotechnical destroy deduction. Additionally, the open pit examined here was excavated in strict accordance with the engineering design. There is a resting angle in the monitored slope; the safety factor is higher when the slope and road surface are at an angle. The result cannot verify the safety angle, so the mapping process should be improved.

The study flow of the paper is as shown in Figure 2. It depends on the assumption that the existing nearest-neighbor interpolation does not consider the incident angle and surface information. We further research to provide a solution to this problem and add the two factors in our geometric and scattering model.

2.2. GB-InSAR and Point Cloud

2.2.1. GB-InSAR Images

The proposed method was validated in an experiment at an open-pit mine in China using a Ku-band GB-InSAR system developed by the China Academy of Safety Science and Technology (CASSAT) and named as slope SAR (S-SAR). The system is shown in Figure 3.

The radar sensor gathered data based on the principle of independent propagation and superposition of electromagnetic waves. The data were coherent after processing. A synthetic-aperture radar is generally an antenna array of N numbers of antennas. The echo signals were received by each antenna and then coherently superimposed. The common GB-InSAR may use an stepped frequency continuous wave (SFCW) signal or linear frequency modulation continuous wave (LFMCW) signal. The SFCW system is easier than the LFMCW to implement, hardware-wise. The linear rail driver adopts the pulse-width modulation (PWM) signal to trigger and control the sliding platform of the linear track "step-by-step and stop". In a single scan from start to end, the original data as-obtained were processed and saved in a complex array. Each element

C_{1}, C_{2}, \dots, C_{n}

in the array corresponds to a specific location on the rail, allowing the data of a certain frequency band to be obtained. The term “repeat-track”, here, refers to passing the same area repeatedly. In this approach, the system acquired SAR data sets by passing the same area twice, covering it with a slightly different viewing geometry. The imaging geometry of the repeat-track SAR interferometry is shown in Figure 3c. On the first pass, the radar wave was transmitted from the transmitting antenna every 1 s. After interaction with the terrain, the backscattered return was also recorded by receiving antenna every 1 s. The signal was then processed to a complex SAR image, with 2 s on the second pass; the signal was processed to another complex SAR image. The system made repeated observations of specific scenes at different time points to obtain radar images in a time series. The phase changes of the image sequence can be analyzed to obtain high-precision deformation information of the observed scene. The elevation in the monitored scene fluctuates as 2D images were gathered in 3D space. The GB-InSAR used an overhead perspective to monitor unstable slopes.

Figure 4 shows the characteristics of the smallest element in a radar image, which is the key parameter in the matching process. The GB-InSAR image grid is the projected image of the observed scene on the slant range plane. After the image co-registration of the GB-InSAR and single look complex image (SLC) in the precise orbit mode, images with the same (M, N) indexed pixel, coherence, and deformation characteristics continuously correspond to the same side of the target slope in the ideal state. In monitoring unstable slopes, the local area observed by GB-InSAR usually has its own 3D coordinate system as well as the topographic data obtained by TLS, total stations, GPS, or interferometric SAR. The TLS is fast and accurate; its monitoring scope is similar to that of the GB-InSAR. (Rigel’s vz4000 model has a mapping range of 4000 m, while the VZ6000 can reach 6000 m in a better environment [17]). In most cases, the digital surface model (DSM) can be rapidly constructed via TLS. The digital elevation model (DEM) was generated after elevation correction. Automatic deformation monitoring with GB-InSAR works via 2D grid imaging methods, which are stable, rapid, and well-tested. This section focuses on the preprocessing of the topographic point cloud as 2D polar coordinate information. This stepwise process included point cloud preprocessing followed by the acquisition of the distance and perspective of spatial points relative to the GB-InSAR.

2.2.2. Point Cloud

As shown in Figure 5a, the distributed ideal targets were selected as the reference points

A_{1} (x_{1}, y_{1}, z_{1})

and

A_{2} (x_{2}, y_{2}, z_{2})

on the 3D terrain slope surface located on both sides of the central axis of the synthetic aperture.

{\begin{cases} r_{1} = \sqrt{{(x_{1} - x_{s})}^{2} + {(y_{1} - y_{s})}^{2} + {(z_{1} - z_{s})}^{2}} \\ = {‖ S - A_{1} ‖}_{2} \\ r_{2} = \sqrt{{(x_{2} - x_{s})}^{2} + {(y_{2} - y_{s})}^{2} + {(z_{2} - z_{s})}^{2}} \\ = {‖ S - A_{2} ‖}_{2} \end{cases}

(1)

where

{‖ \cdot ‖}_{2}

represents the Euclidean 2-norm of the space vector. As shown in Figure 5b, the local azimuth angle in the 2D plane coordinate system can be calculated as follows:

∠ R_{0} S A = \arccos (\frac{| A_{k} A |}{| A_{k} S |}),

(2)

where

R_{0}

represents the center of the monitoring area,

∠ R_{0} S A

is shown in Figure 5b which denotes for relative view angle of target

A

. The angle needs a further conversion according to the specified radar image coordinates when geometric mapping.

3. Methods

3.1. Geometric Mapping Between Image Space and Terrain Space

Let

I

denote for the deformation map 2D matrix. The size of

I

(M, N)

and is coded in polar coordinates

(r, θ)

or pseudo-polar coordinates

(r, \sin (θ))

. Let

I_{p c l}

represent the 2D matrix formed by the relative range and relative azimuth angle derived by Equations (1) and (2). There is no existing deformation value corresponding to the internal coordinates of

I_{p c l}

I

. Then, the mapping relation we have to calculate is

f : I \to I_{p c l}

. In our case, the spatial resolution of the point cloud is higher than that of radar, the appropriate method of resampling and interpolation is needed to assign deformation value to the point cloud. Mapping relations are as follows:

{\begin{array}{l} \begin{array}{l} f 1 : (r, \sin (θ)) \to (m, n) \\ f 2 : (r, θ) \to (X_{1}, Y_{1}, Z_{1}) \\ f 3 : (X_{1}, Y_{1}, Z_{1}) \to (X_{t}, Y_{t}, Z_{t}) \end{array} \\ f = f 1 \circ f 2 \circ f 3 \end{array}

(3)

where

f 1

represents the mapping of pseudo-polar coordinates (e.g., plane coordinates) to pixel coordinates,

(m, n)

represents the 2D image pixel index.

f 2

is the mapping of polar coordinates (e.g., plane coordinates) to the local coordinate system of the point cloud, and

(X_{1}, Y_{1}, Z_{1})

denotes local coordinates.

f 3

is the mapping of the local coordinate system of the point cloud to the geodetic coordinate system,

(X_{t}, Y_{t}, Z_{t})

stands for the geodetic coordinates and symbol ‘

\circ

’ represents a cascade.

3.2. Geometric and Scattering Weight Model

As shown in Figure 6, slope surface nodes and their nearest neighbor discrete points in the DEM consist of small surface element patches. Traditional matching methods mostly work based on the geometric model of the radar and terrain. Here, a geometric and scattering model was applied for the patches to solve the one-to-many problem.

The scattering model is commonly used in SAR image simulations. The problem to be solved, as discussed in this section, is similar to that of SAR image simulation: The terrain point cloud data imported externally have higher dimensions and resolution than the radar data. Deformation results with high dimensionality and resolution can be retrieved based on the existing GB-InSAR monitoring results. The most commonly used scattering model is the small surface element model [18], each small surface element is tangent to the actual plane. Each surface must be small enough to effectively approximate the surface fluctuation. The spatial geometry parameters of a small surface element are determined by its central position vector and normal vector. The open-pit mine in this case study was subject to daily blasting operations which also shook the command center (building in Figure 6). Slope stability in the lower part of the command center is the primary target of mine safety applications.

An algorithm flow was established (Figure 7) to explicitly show how to determine the weights assigned to the point cloud by the radar pixel based on the geometric-scattering model. K values need to be set according to the difference in resolution between the point cloud and the radar image. The normal vector can be computed using the method introduced by Hoppe et al. based on the directed distance function (signed short function) of the scattered point cloud the normal vector [19]. Assuming that the terrain point cloud plane is continuous after sampling, the most similar K adjacent points can be determined for each point p by the k-nearest neighbor classification method, which has been widely used, and as such there is no need to describe it in this paper. The local incident angle can be obtained quickly via Equation (4). The focal point coordinates

S (x_{1}, y_{1}, z_{1})

and the phase center coordinates

A (x_{2}, y_{2}, z_{2})

of the small surface element allow the local incident angle

θ_{i}

for each small surface element to be computed as:

θ_{i} = \arccos \frac{- \vec{n_{i}} \cdot \vec{S A_{i}}}{| {\vec{n}}_{i} | | \vec{S A_{i}} |},

(4)

where

{\vec{n}}_{i}

represents the normal vector of the ith point in the terrain cloud.

The steepness in the k-neighborhood of each point on the surface can be calculated with a differential operation. A smaller value indicates flatter terrain in the neighborhood and vice versa. The slope of the k-neighborhood at each point can be calculated as follows:

θ_{s l o p e} = \arctan \sqrt{f_{x}^{2} + f_{y}^{2}},

(5)

where

θ_{s l o p e}

is the radian value of the slope angle in the k neighborhood,

f_{x}

is the change rate of elevation in the x direction, and

f_{y}

is the change rate of elevation in the y direction. The direction vector of the small face element is calculated based on the slope of the k-neighborhood as:

θ_{a s p e c t} = \arctan \sqrt{\frac{f_{y}}{f_{x}}},

(6)

where

θ_{a s p e c t}

is the orientation of the k-neighborhood. The orientation represents the maximum change direction of the elevation change at any point in the terrain point cloud. The standard deviation of surface heights was used to evaluate the degree to which the surface fluctuates and to judge whether the k-nearest neighbor surface element scattering model can be used. Suppose a surface is in the x–y plane and the height of some point

(x, y)

above this plane is

z (x, y)

. Assuming that the center of the surface is at the origin, then the average height of the surface can be calculated by Equations (7) and (8).

\bar{z} = \frac{1}{L_{x} L_{y}} \int_{- L_{x} / 2}^{L_{x} / 2} \int_{- L_{y} / 2}^{L_{y} / 2} z (x, y) d x d y,

(7)

The corresponding secondary moment is:

\bar{z^{2}} = \frac{1}{L_{x} L_{y}} \int_{- L_{x} / 2}^{L_{x} / 2} \int_{- L_{y} / 2}^{L_{y} / 2} z^{2} (x, y) d x d y,

(8)

The standard deviation of the surface heights

σ

σ = {(\bar{z^{2}} - {\bar{z}}^{2})}^{\frac{1}{2}},

(9)

The Fraunhofer criterion of the antenna’s far-field condition is also useful here. If the maximum phase difference between the center and edge of the antenna and the receiving point is greater than or equal to

π / 8

, it can serve as the criterion of the scattering model.

σ \geq \frac{λ}{32 \cos θ_{i}},

(10)

where

λ

is the wavelength. The direction vector of orientation can replace the normal vector of k-nearest neighbors in Equation (4) to calculate the point’s k-neighborhood relative to the radar incidence angle:

θ_{i} = \arccos \frac{- n_{a s p e c t} \cdot \vec{S A}}{| n_{a s p e c t} | | \vec{S A} |},

(11)

where

\vec{S A}

is explicated in Figure 5a. The equation of the upper hemisphere using Lambert’s law is the weight assigned to the cosine graph according to the distribution of Equation (12):

γ (θ_{i}) = γ (0) \cos θ_{i},

(12)

where

γ (0)

is the initial value of the weight can be determined by other methods with distinct geometric features. The scattering observed here is caused by the projected area. The spheres are independent of each other and separated far enough so that they do not overlap in any way. Therefore, the re-radiation directivity diagram of the isotropic spheres remains unchanged.

The next section shows how to use the proposed process with simulated data, and then verifies the proposed method with actual data.

4. Results

A slope-pavement model was established (Figure 8a). The step model consists of slope and pavement components. The deformable regions were simulated across the road and slope on a deformation image based on the simulated monitored slope. The nearest-neighbor interpolation method and the proposed method were used separately for mapping to observe whether the projected deformation was consistent with the shape of the step. The model has been geocoded with the target area based on block stretching and linear transformation (Figure 8b).

Normal $\vec{n_{i}}$ : Normal vector computing is integrated into many point cloud processing software and C++ library, for example: Cloud compare, MeshLab, and point cloud library (PCL). We programmed with the PCL library and referenced the Cloud compare built-in minimum spanning tree method to extract normal vectors of the sampled step model. $\vec{n_{i}}$ is as shown in Figure 8c. The pavement and the slope can be differentiated by $\vec{n_{i}}$ .
Orientation: The orientation of each point $n_{a s p e c t}$ in the model is easy to extract. As the step model coded, it has a certain horizontal coordinate system plane. The angle between $\vec{n_{i}}$ and horizontal plane normal can be treated as an orientation vector $n_{a s p e c t}$ , and the $θ_{a s p e c t}$ of each point in Equation (6) is identified.
Line-of-sight (LOS): LOS is the vector that connects each point to radar Station S.
Relative incident angle $θ_{i}$ : $θ_{i}$ can be determined by Equation (11).
Weight $γ (θ_{i})$ : $γ (θ_{i})$ is the normalization of $\cos (θ_{i})$ . The result is shown in Figure 8d.

One rather important thing to note, the simulated stepped model is more like a “smooth surface” consist of many particles. However, the target area in a natural scene is more like a “rough surface”. Figure 8d shows the weight difference of the stepped model. The pavement weight was much lower than the slope surface through SAR beam incidents in an overlooking perspective. The weight of pavement only means the incident angle is not good, which may lead to a disappearance of deformation in the geometric mapping process. The topographically rugged surface of the pavement will not be affected by the deformation disappearance.

Sub-area deformation positioning: The mapping relation table of the 3D point cloud and radar image was established via geometric mapping. Several subregions were identified in the diagram that are composed of groups of pixels that appear as blocks on the surface. As shown in Figure 9a, these blocks have a global displacement along the LOS direction; this simulates the displacement detected by the radar rather than changes in the 3D model. The simulation results can then be mapped to a 3D space as shown in Figure 9b. This case features sub-area overall deformation. The GB-InSAR image as-collected from the open-pit mine was used to match the real terrain point cloud and simulated terrain point cloud using the method described in Section 2. Several sub-areas that cover the pavement and slope surface were then created as shown in Figure 9a–c were the index of the sub-area. However, in the real world, the deformation of the road surface and slope surface is always different, so simulation case 1 only shows the spatial positioning of each sub-area.

Continuous deformation generated through simulation in Region C: This was a simulation of a slope surface collapse, where the flow of water in the slope was caused by internal empty surface sag. The deformation field was simulated via a quadratic function model, as shown in Figure 10a. At this point, the lower part of the slope moved down relative to the radar resulting in a negative displacement relative to the radar. This was then projected into 3D space using the nearest-neighbor method (Figure 10b). If the area produces persistent, identical depressions and is detected by GB-InSAR, then the subsequent deformation image is the same as the image shown in Figure 10c. Then the deformation values of each pixel can be accumulated to form the image. At this point, there is a gradual and non-linear deterioration of the slope. If a pixel is across the slope and road surface, the assignment of the point cloud within the pixel set should be distinguishable. Scattering model analysis can solve this problem. If the radar wave incident angle of the road is poor but that of the slope is acceptable, then the nearest-neighbor method results of the road surface and slope are the same. The proposed method makes the road result more distinguishable by reducing its weight and strengthening slope results. Based on the above assumptions, the normal vector of each point was calculated followed by the slope angle of the small surface element adjacent to each point.

Target area real application: Firstly, the pavement and slope scattering variance can be extracted by the S-SAR microwave incident angle and is colored as shown in Figure 11a. The bad incident areas have no color. The example figure of normal and calculated

θ_{a s p e c t}

is shown in Figure 11b,d. The

θ_{a s p e c t}

shows that the sub-area would most probably move towards a given direction when affected by gravitation. However, when they are affected by other kinds of outside forces, the model should couple with the force analysis. Since the target area is a permanent slope which is only fluctuated by daily mining explosion, it can be treated as deforms along

n_{a s p e c t}

vector. The weight

γ (θ_{i})

variance is shown in Figure 11c. The areas with bad incident angles are shown as deep blue and brighter yellow areas are areas with good incident angles. As the road surface height went lower, the incident situation got better. The dark color of each road edge is the soil retaining wall which is to protect trucks from falling. The wall area weights also vary by incident angle changes. The nearest-neighbor mapping result and paper method result is as shown in Figure 12.

Three sub-areas A, A’, B, B’, C. C’ were selected, corresponding to the simulated step model sub-areas. Area A is as shown in Figure 13a, the bright yellow boundary shows a great deformation, with the highest reach near 56 mm, which is quite large. The situation is quite similar to the controlled field test from in our previous work [1], which are abnormal boundaries caused by the persistent scattering points that affect the neighborhood area. The area itself is not in a good coherence situation and the phase error is significant. Then, the deformation extracted in these boundary areas is also abnormal. In the A’ area, as is shown in Figure 13b, the situation is relieved, especially for the pavement areas: For the red point, from 55.5943 mm to 2.4283 mm; for the black point from 47.4978 mm to 18.6112 mm. For the slope area, the real deformation still exists, a geologist can introduce their analysis method for further destruction deduction.

GB-InSAR monitoring is flexible, and monitoring geometry should be well designed to balance between layover distortion and shadow. The local incident angle and K-nearest neighbor patches aspect angle in the geometric and scattering model may help the layover distortion analysis (the case is shown in Figure 14). The layover distortion would mostly happen when the SAR beam incident angle is near perpendicular. The area B and B’ show the same radar spatial resolution element. The B area shows no difference when across the road and slope. The sub-area deforms about 20 mm during the experiment. However, for the paper method, the B’ showed a difference as the incident angle changes. The road and slope were differentiated, and the inside small geo-elements deform according to their local shape. The result would help small area (less < 200 m²) slide analysis, which often happens in open-pit mines.

As shown in Figure 15, the areas C and C’ show an area with an undulated surface. The case is quite similar to the road and slope crossing area. These areas are commonly seen in an open pit where roads go two directions, or the roads are a regional rotary; these are often unstable when trucks pass. The abnormal boundary was relieved by the proposed method in the paper. The nearest-neighbor method result is confused by abnormal deformation and the real deformed area is concealed.

5. Discussion

This section concentrates on the following issues:

How the weight influences the deformation mapping result.
How can we get closer to the real deformation.

In a strict sense, the weight calculated based on Equation (12) in the paper may not conform to the scattering differences of local areas in any natural scene. However, it is widely used in SAR image simulation. It comes from the hypothesis that the surface is a collection of “small spherical” scatterers. Some scattering model is simulated based on

\cos θ

function and a lot of researchers adopt the

\cos^{2} θ

function to simulate the scattering difference. In addition, if possible, we need to design our own compromise function to accommodate the local geological condition. The factor will affect deformation mapping. For example, the result in Figure 13 shows sudden changes from 40–50 mm to 20–30 mm. The situation may confuse us if the mapping result is correct or near the real deformation. Figure 16a,b is the one radar pixel interpolated deformation surf plot of A and A’ (Figure 13). The figure would help to understand the influences when using a geometric-scattering model.

According to Figure 16, the following understanding can be concluded.

The GB-InSAR image pixel contains several targets. Due to the limited spatial resolution, these targets can be equivalent to one vital “scattering target”. The terrain point cloud model spatial resolution is higher than radar. The local area with a good incident angle in the radar pixel plays a major role in forming "sub-targets", and the shape variables are distributed on these distributed sub-targets. In this study, the deformation on the road surface is not eliminated, and their deformation is concentrated on the distributed sub-targets with good incident angles on the road surface.
If the precise deformation in the radar pixel is known through the ground control points or global navigation satellite system (GNSS) and other high time resolution measuring instruments, a more accurate model of mapping can be established. The optimizing objective function the least square solution of the temporal sequence deformation of each sub-target and the deformation function of the radar pixel.

From the perspective of geotechnical application, the difference between slope and road surface is highlighted, and the instability can be analyzed from a macroscopic perspective. Moreover, the model proposed in this paper can be integrated with the existing classical methods, rather than abandoning the previous methods. Fusion analysis is not the focus of this paper, and it will be further studied in the future.

6. Conclusions

In this study, geometric and scattering models were established to solve the “one-to-many” problem and pavement, slope differentiation. A scattering and geometric weighting model was proposed during the interpolation process of mapping. A simulation case and a target area field test were conducted to test the method and compared it with the nearest-neighbor interpolation method. The following problems were resolved:

The pavement and slope surface deformation were differentiated.
The parameters can be adjusted to avoid band-like phenomena in the experiment.
The abnormal deformed boundaries were relieved to a certain extent.

The geometric and scattering models are not without defects. In the future, more adaptable models should be introduced according to regional conditions of the slope under analysis so that the deformation measured by GB-InSAR can be best integrated with other geographical models. This will allow geologists to analyze the safety of unstable slopes.

Author Contributions

X.Z. carried out the simulation, experiment, and data analyses. He also wrote the main portion of the manuscript. X.H. reviewed and contributed to the text. X.Y. was responsible for the methodology and imaging process of GB-InSAR. H.M. proposed the correction conceptualization. Z.Y. validated the whole integration scheme. G.R. managed the data curation and supervised the GB-InSAR monitoring. J.L. processed the point cloud and performed TLS modeling. H.Z. wrote the SAR image preprocessing algorithm. J.Z. proposed advice and took part in the writing of the text. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by [Natural Science Foundation of China] grant number [41830110]. And the APC was funded by [the Fundamental Research Funds for the Central Universities] grant number [B200202011]. And [Basic Scientific Research Project of China Academy of Safety Science and Technology] grant number 2020JBKY07.

Conflicts of Interest

The authors declare no conflict of interest.

References

Zheng, X.; Yang, X.; Ma, H.; Ren, G.; Zhang, K.; Yang, F.; Li, C. Integrated Ground-Based SAR Interferometry, Terrestrial Laser Scanner, and Corner Reflector Deformation Experiments. Sensors 2018, 18, 4401. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Pieraccini, M.; Miccinesi, L. Ground-Based Radar Interferometry: A Bibliographic Review. Remote Sens. 2019, 11, 1029. [Google Scholar] [CrossRef] [Green Version]
Carla, T.; Tofani, V.; Lombardi, L.; Raspini, F.; Bianchini, S.; Bertolo, D.; Thuegaz, P.; Casagli, N. Combination of GNSS, satellite InSAR, and GBInSAR remote sensing monitoring to improve the understanding of a large landslide in high alpine environment. Geomorphology 2019, 335, 62–75. [Google Scholar] [CrossRef]
Xiaolin, Y.; Yanping, W.; Yaolong, Q.; Weixian, T.; Wen, H. Experiment study on deformation monitoring using ground-based SAR. Synthetic Aperture Radar. In Proceedings of the Asia-Pacific Conference on Synthetic Aperture Radar (APSAR), Tsukuba, Japan, 23–27 September 2013. [Google Scholar]
Tapete, D.; Casagli, N.; Luzi, G.; Fanti, R.; Gigli, G.; Leva, D. Integrating radar and laser-based remote sensing techniques for monitoring structural deformation of archaeological monuments. J. Archaeol. Sci. 2013, 40, 176–189. [Google Scholar] [CrossRef] [Green Version]
Pieraccini, M.; Miccinesi, L. An Interferometric MIMO Radar for Bridge Monitoring. IEEE Geosci. Remote Sens. Lett. 2019, 16, 1383–1387. [Google Scholar] [CrossRef]
Intrieri, E.; Gigli, G.; Nocentini, M.; Lombardi, L.; Mugnai, F.; Fidolini, F.; Casagli, N. Sinkhole monitoring and early warning: An experimental and successful GB-InSAR application. Geomorphology 2015, 241, 304–314. [Google Scholar] [CrossRef] [Green Version]
Luzi, G.; Pieraccini, M.; Mecatti, D.; Noferini, L.; Macaluso, G.; Tamburini, A.; Atzeni, C. Monitoring of an Alpine Glacier by Means of Ground-Based SAR Interferometry. IEEE Geosci. Remote Sens. Lett. 2007, 4, 495–499. [Google Scholar] [CrossRef]
Wang, P.; Xing, C. Research on Coordinate Transformation Method of GB-SAR Image Supported by 3D Laser Scanning Technology. Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci. 2018, XLII-3, 1757–1763. [Google Scholar] [CrossRef] [Green Version]
Hu, C.; Deng, Y.; Wang, R.; Tian, W.; Zeng, T. Two-Dimensional Deformation Measurement Based on Multiple Aperture Interferometry in GB-SAR. IEEE Geosci. Remote Sens. Lett. 2016, 14, 208–212. [Google Scholar] [CrossRef]
Sansosti, E.; Berardino, P.; Manunta, M.; Serafino, F.; Fornaro, G. Geometrical SAR image registration. IEEE Trans. Geosci. Remote Sens. 2006, 44, 2861–2870. [Google Scholar] [CrossRef]
Yue, J.; Yue, S.; Wang, X.; Guo, L. Research on multi-source data integration and the extraction of three-dimensional displacement field based on GBSAR. In Proceedings of the International Workshop on Thin Films for Electronics, Electro-Optics, Energy and Sensors, Suzhou, China, 4–6 July 2015. [Google Scholar]
Zou, J.; Tian, J.; Chen, Y.; Mao, Q.; Li, Q. Research on Data Fusion Method of Ground-Based SAR and 3D Laser Scanning. J. Geomat. 2015, 40, 26–30. [Google Scholar]
Yang, J.; Qi, Y.L.; Tan, W.X.; Yanping, W.; Wen, H. Three-dimensional matching algorithm for geometric mapping between GB-SAR image and terrain data. J. Univ. Chin. Acad. Sci. 2015, 32, 422–427. (In Chinese) [Google Scholar]
Lombardi, L.; Nocentini, M.; Frodella, W.; Nolesini, T.; Bardi, F.; Intrieri, E.; Carlà, T.; Solari, L.; Dotta, G.; Ferrigno, F.; et al. The Calatabiano landslide (southern Italy): Preliminary GB-InSAR monitoring data and remote 3D mapping. Landslides 2017, 14, 685–696. [Google Scholar] [CrossRef] [Green Version]
Zheng, X.; Yang, X.; Ma, H.; Ren, G.; Yu, Z.; Yang, F.; Zhang, H.; Gao, W. Integrative Landslide Emergency Monitoring Scheme Based on GB-INSAR Interferometry, Terrestrial Laser Scanning and UAV Photography. J. Phys. Conf. Ser. 2019, 1213. [Google Scholar] [CrossRef]
Available online: http://www.riegl.com/nc/products/terrestrial-scanning/ (accessed on 28 March 2020).
Ulaby, F.; Long, D.; Blackwell, W.; Elachi, C.; Fung, A.K.; Ruf, C.; Sarabandi, K.; van Zyl, J.; Zebker, H. Microwave Radar and Radiometric Remote Sensing; University of Michigan Press: Ann Arbor, MI, USA, 2014. [Google Scholar]
Hoppe, H.; Derose, T.; Duchamp, T.; McDonald, J.; Stuetzle, W. Surface reconstruction from unorganized points. In Proceedings of the 19th Annual Conference on Computer Graphics and Interactive Techniques, Chicago, IL, USA, 26–31 July 1992; Volume 26, pp. 71–78. [Google Scholar]

Figure 1. (a) Optical image of the target open pit slope. (b) It is difficult to distinguish the pavement surface from the slope surface using the mapping result generated by nearest-neighbor interpolation. The orange color areas consist of several pixels across multiple slope steps. A ground-based synthetic aperture radar interferometry (GB-InSAR) cumulative deformation map was produced in a month, and was mapped on the terrain model. The orange-colored area shows great deformation but is difficult to interpret by geotechnical engineers.

Figure 2. Logical experiment scheme of the data fusion between GB-InSAR interferograms and terrain point cloud data using nearest-neighbor interpolation and geometric-scattering model separately.

Figure 3. The synthetic aperture model was observed at a fixed-station, overhead perspective. (a) The key process of GB-InSAR echo data processing. (b) S-SAR radar system which contains transmitting and receiving antennas, radar sensor and a linear rail. (c) Ideal state: "One-step, one-stop" SFCW signal transmitting and receiving mode of straight-line repeated-track controlled by a pulse-width modulation (PWM).

Figure 4. (a) The GB-InSAR image grid. The spatial resolution of the grid is lower than terrestrial laser scanning (TLS) points cloud resolution. Each grid may correspond to multiple terrain surface model points. (b) GB-InSAR image pixel responds to the same slope area after 2D image co-registration.

Figure 5. (a) The geometric mapping relationship between normal center rectangular coordinate system GB-InSAR and 3D terrain data. (b) Terrain point cloud azimuth angle top view.

Figure 6. Terrain point cloud introduces greater detail so that radar wave incident angles define geometric and scattering models. The model helps in solving the one-to-many problem during the data fusion process. For the road surface and slope surface of the slope, they have completely different relative incident angles.

Figure 7. Assignment weights according to point cloud incident angle geometric-scattering model.

Figure 8. (a) Surface color of the road-slope model indicating

θ_{s l o p e}

variance extracted by Equation (5). (b) Geocoded model with the actual collected terrain point cloud to a local coordinate system via linear transformation. (c) Normal vectors of the pavement and slope of the rectangle. (d)

γ (θ_{i})

the variance of each point. The deep blue color area is the pavement. In this case, they have a bad incident angle.

Figure 8. (a) Surface color of the road-slope model indicating

θ_{s l o p e}

γ (θ_{i})

the variance of each point. The deep blue color area is the pavement. In this case, they have a bad incident angle.

Figure 9. (a) Simulated sub-areas with small deformation 2D deformation map. (b) Geometrical mapping 3D result of deformation in (a) using a nearest-neighbor interpolation method.

Figure 10. (a) 2D deformation diagram shows a slight caving-in phenomenon. (b) The nearest-neighbor method creates a 3D result which has no obvious discrimination between slope and road surfaces. (c) Simulated deformation diagram of Region D after accumulation. (d) Accumulated 3D matching produces obvious banding. A similar phenomenon occurs in adjacent regions to D. (e) Nearest-neighbor method result. (f) The result generated by the geometric and scattering weight model, where deformation results distinguish road surface from the slope surface. The pavement deformation disappears because the incident angle is not good. If the pavement surface is rugged as shown in Figure 6, the pavement deformation would be reserved.

Figure 11. (a) Radar local incident angle and (b) normalized results, local zoom in Figure 11a. The arrow represents the direction of the incident vector (Equation (11)). (c) Weight varies when the scattering model changes. (d) The arrow represents the direction of local most likely deformed by gravitation or explosion fluctuation.

Figure 12. The cumulative deformation map from GMT+8 13:03 10 October 2018 to 17:23 was used in the mapping experiment. Three sub-areas were also selected. (a) Nearest-neighbor interpolation mapping result. (b) Paper method mapping result.

Figure 13. The deformation of the transition zone is removed to a certain extent by the model weighting method of the paper. These edge boundary anomalies, however, are often marginal areas with poor data quality, but are affected by strong coherence scattering points. (a) Nearest-neighbor interpolation mapping result of sub-area A. (b) Paper method mapping result of sub-area A’.

Figure 14. Multiple terrain surface deformed points may locate in one radar pixel spatial resolution element which crosses pavement and slope. The local layover distortion is hard to identify without geometric and scattering information. (a) Nearest-neighbor interpolation mapping result of sub-area B. (b) Paper method mapping result of sub-area B’.

Figure 15. A deformation region with large local undulations. (a) Nearest-neighbor interpolation mapping result of sub-area C. (b) Paper method mapping result of sub-area C’.

Figure 16. The deformation in one radar pixel generated by two methods. (a) Nearest-neighbor interpolation deformation result. The whole pixel may equal one vital scattering target. (b) Deformation result using the geometric and scattering model. The terrain model point cloud may equal several distributed targets.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Zheng, X.; He, X.; Yang, X.; Ma, H.; Yu, Z.; Ren, G.; Li, J.; Zhang, H.; Zhang, J. Terrain Point Cloud Assisted GB-InSAR Slope and Pavement Deformation Differentiate Method in an Open-Pit Mine. Sensors 2020, 20, 2337. https://doi.org/10.3390/s20082337

AMA Style

Zheng X, He X, Yang X, Ma H, Yu Z, Ren G, Li J, Zhang H, Zhang J. Terrain Point Cloud Assisted GB-InSAR Slope and Pavement Deformation Differentiate Method in an Open-Pit Mine. Sensors. 2020; 20(8):2337. https://doi.org/10.3390/s20082337

Chicago/Turabian Style

Zheng, Xiangtian, Xiufeng He, Xiaolin Yang, Haitao Ma, Zhengxing Yu, Guiwen Ren, Jiang Li, Hao Zhang, and Jinsong Zhang. 2020. "Terrain Point Cloud Assisted GB-InSAR Slope and Pavement Deformation Differentiate Method in an Open-Pit Mine" Sensors 20, no. 8: 2337. https://doi.org/10.3390/s20082337

APA Style

Zheng, X., He, X., Yang, X., Ma, H., Yu, Z., Ren, G., Li, J., Zhang, H., & Zhang, J. (2020). Terrain Point Cloud Assisted GB-InSAR Slope and Pavement Deformation Differentiate Method in an Open-Pit Mine. Sensors, 20(8), 2337. https://doi.org/10.3390/s20082337

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

Terrain Point Cloud Assisted GB-InSAR Slope and Pavement Deformation Differentiate Method in an Open-Pit Mine

Abstract

1. Introduction

2. Linear GB-InSAR and TLS Data Fusion

2.1. Problem and Solving Scheme

2.2. GB-InSAR and Point Cloud

2.2.1. GB-InSAR Images

2.2.2. Point Cloud

3. Methods

3.1. Geometric Mapping Between Image Space and Terrain Space

3.2. Geometric and Scattering Weight Model

4. Results

5. Discussion

6. Conclusions

Author Contributions

Funding

Conflicts of Interest

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI