[go: up one dir, main page]

Next Article in Journal
Towards a Digital Twin Prototype of Alpine Glaciers: Proposal for a Possible Theoretical Framework
Previous Article in Journal
Preliminary Results of the Three-Dimensional Plasma Drift Velocity at East Asian Low-Latitudes Observed by the Sanya Incoherent Scattering Radar (SYISR)
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Large-Scale Land Subsidence Monitoring and Prediction Based on SBAS-InSAR Technology with Time-Series Sentinel-1A Satellite Data

1
National Supercomputing Center in Zhengzhou, Zhengzhou University, Zhengzhou 450001, China
2
School of Geoscience and Technology, Zhengzhou University, Zhengzhou 450001, China
3
Henan Institute of Geological Survey, Zhengzhou 450001, China
4
National Engineering Laboratory Geological Remote Sensing Center for Remote Sensing Satellite Application, Zhengzhou 450001, China
*
Author to whom correspondence should be addressed.
Remote Sens. 2023, 15(11), 2843; https://doi.org/10.3390/rs15112843
Submission received: 18 April 2023 / Revised: 29 May 2023 / Accepted: 29 May 2023 / Published: 30 May 2023
(This article belongs to the Section Engineering Remote Sensing)
Figure 1
<p>The geographic location of the study area, the coverage of the Sentinel-1 image tiles, and the distribution of leveling benchmarks are superimposed on the ESRI satellite topographic map.</p> ">
Figure 2
<p>List of specific dates of SAR image data.</p> ">
Figure 3
<p>The method flowchart.</p> ">
Figure 4
<p>Partial SAR fusion experiment results for track 40 and track 142. (<b>a</b>) indicates the SAR experimental results of track 40 in the overlapping region, (<b>b</b>) indicates the SAR experimental results of track 142 in the overlapping region, and (<b>c</b>) represents the final fusion results after adding the offset.</p> ">
Figure 5
<p>Schematic diagram of the network unit structure of the LSTM model.</p> ">
Figure 6
<p>Architectural diagrams of the LSTM model (<b>a</b>), TCN model (<b>b</b>), and LSTM-TCN model (<b>c</b>).</p> ">
Figure 7
<p>The semiannual average spring and summer deformation rates of the land in the study area. The values are obtained by dividing the cumulative deformation values in spring and summer of 2020 and 2021 by 2. Positive values represent the rate of land uplift, and negative values represent the rate of land subsidence.</p> ">
Figure 8
<p>The semiannual average deformation rate of the land in the study area in autumn and winter. The values are obtained by dividing the cumulative deformation values in autumn and winter of 2020 and 2021 by 2. Positive values represent the rate of land uplift, and negative values represent the rate of land subsidence.</p> ">
Figure 9
<p>Correlation diagram between subsidence point data and leveling point measured data in the study area.</p> ">
Figure 10
<p>Average annual land deformation rates from November 2019 to February 2022. A positive number represents land uplift, and a negative number represents land subsidence.</p> ">
Figure 11
<p>Cumulative land deformation in the study area from November 2019 to February 2022. A positive number represents land uplift, and a negative number represents land subsidence.</p> ">
Figure 12
<p>Geographical distribution map of all experimental subsidence points.</p> ">
Figure 13
<p>Loss curves of 12 parameter combinations for three models. (<b>a</b>) Training Loss and (<b>b</b>) Validation Loss.</p> ">
Figure 14
<p>The correlation coefficient graph between the predictions of the three models and the results of SBAS-InSAR. (<b>a</b>–<b>c</b>) Correlation diagrams of the LSTM, TCN, and LSTM-TCN models, respectively.</p> ">
Figure 15
<p>Location map of the four validation points.</p> ">
Figure 16
<p>Comparison of the prediction results of points (<b>a</b>–<b>d</b>) with the monitoring results of SBAS-InSAR.</p> ">
Figure 17
<p>Experimental results of the three models with different numbers of samples. (<b>a</b>–<b>d</b>) represent the RMSE, <math display="inline"><semantics> <mrow> <msup> <mi>R</mi> <mn>2</mn> </msup> </mrow> </semantics></math>, MAE, and MAPE results of the three models, respectively.</p> ">
Figure 18
<p>Influence of different monthly sequence input matrices on experimental results. (<b>a</b>–<b>d</b>) represent the RMSE, <math display="inline"><semantics> <mrow> <msup> <mi>R</mi> <mn>2</mn> </msup> </mrow> </semantics></math>, MAE, and MAPE results of the three models, respectively.</p> ">
Review Reports Versions Notes

Abstract

:
Rapid urban development in China has aggravated land subsidence, which poses a potential threat to sustainable urban development. It is imperative to monitor and predict land subsidence over large areas. To address these issues, we chose Henan Province as the study area and applied small baseline subset-interferometric synthetic aperture radar (SBAS-InSAR) technology to obtain land deformation information for monitoring land subsidence from November 2019 to February 2022 with 364 multitrack Sentinel-1A satellite images. The current traditional time-series deep learning models suffer from the problems of (1) poor results in extracting a sequence of information that is too long and (2) the inability to extract the feature information between the influence factor and the land subsidence well. Therefore, a long short-term memory-temporal convolutional network (LSTM-TCN) deep learning model was proposed in order to predict land subsidence and explore the influence of environmental factors, such as the volumetric soil water layer and monthly precipitation, on land subsidence in this study. We used leveling data to verify the effectiveness of SBAS-InSAR in land subsidence monitoring. The results of SBAS-InSAR showed that the land subsidence in Henan Province was obvious and uneven in spatial distribution. The maximum subsidence velocity was −94.54 mm/a, and the uplift velocity was 41.23 mm/a during the monitoring period. Simultaneously, the land subsidence in the study area presented seasonal changes. The rate of land subsidence in spring and summer was greater than that in autumn and winter. The prediction accuracy of the LSTM-TCN model was significantly better than that of the individual LSTM and TCN models because it fully combined their advantages. In addition, the prediction accuracies, with the addition of environmental factors, were improved compared with those using only time-series subsidence information.

1. Introduction

Land subsidence is a geological issue in which land in a given region loses elevation under soil compression induced by natural or human activities. Humans continue to exploit coal, oil, gas, and groundwater for their own development needs [1,2,3,4,5,6]. Land subsidence has become a geological disaster that cannot be ignored in social and economic development [7,8]. Land subsidence creates large threats and losses for human productivity and life safety and severely hampers sustainable development [9,10]. Therefore, early monitoring and forecasting of land subsidence are of great significance to prevent and control geological disasters [11].
Traditional methods are mainly conventional geodetic observations based on discrete points, such as leveling, GNSS, and LiDAR technology, which are used to obtain highly accurate land deformation information to monitor land subsidence [12,13,14]. However, limited by small-scale monitoring areas, long measurement cycles, and high costs, these traditional approaches cannot meet the requirement of timely monitoring over large areas [15]. The monthly mass concentration change data from the GRACE gravity satellite cannot achieve high-precision monitoring of land subsidence due to its low resolution [16,17,18]. The time series InSAR technology can provide all-weather, high-precision high-resolution land deformation information and has been widely used for large-scale continuous land deformation monitoring [19,20,21,22,23]. Furthermore, InSAR technology has been widely used to detect land deformation caused by groundwater level fluctuations [24,25,26,27]. Differential InSAR (D-InSAR) is a technique for acquiring information about small land deformations by combining 2-view SAR imagery with external topographic data [28]. Time-series InSAR satellite remote sensing technology, such as persistent scatterer-InSAR (PS-InSAR) [29,30,31] and small baseline subset-InSAR (SBAS-InSAR) [32], has been used to estimate land deformation at high spatial resolutions with mm-scale accuracy. In particular, SBAS-InSAR can use all available interferometric pairs to achieve high accuracy even with shorter wavelength SAR data (such as C-band) over large areas, and consequently, it exhibits great potential applications to the field of large-scale and long time series land subsidence monitoring [33,34,35].
There are threemain approaches to land subsidence prediction: physical models [36], mathematical statistical models [37], and machine learning models [38,39,40,41,42]. The physical model predicts land subsidence based on a series of complex geological and hydrological characteristics of physical parameters through field investigations and experiments. However, these models cannot fit the land subsidence trend well because they require that assumptions be made about complex parameters in order to construct the corresponding models, which may fail at any time [43]. Mathematical statistical methods use numerical and empirical models, such as regression analysis and gray prediction, to predict land subsidence based on the analysis of internal relationships of a large amount of historical monitoring data [44]. However, it is difficult to generalize this method because it does not consider the internal structural relationship between underground rock and soil media. Machine learning methods can achieve high-accuracy prediction of land subsidence by mining the trend of historical land subsidence data and learning the internal nonlinear characteristics of land subsidence [45]. However, the current machine learning algorithms fail to adequately consider the information on the association characteristics between the influencing factors and land subsidence, and they do not show better prediction results for a particular study area. In fact, subsidence at any site is related to many factors, such as geological structure, soil volume water content, and rainfall distribution. If the correlation characteristics between these factors and land subsidence are not considered properly, the constructed prediction models may be invalid in local areas [46]. In particular, the long short-term memory (LSTM) is specifically designed for sequence prediction problems and has been extensively used in land subsidence prediction [47]. However, due to the serialized nature of LSTM data processing, only critical historical information can be obtained and retained by rolling in one direction, and only a limited length of data information can be obtained. A temporal convolutional network (TCN) is made up of parallelized convolutional networks that can process massive amounts of long sequence data more efficiently. Moreover, its dilated convolutional architecture in the network can flexibly control the size of the receptive field so that the network can process longer time series with certain parameters [48]. Compared with the LSTM model, the TCN model can extract the characteristic information between the influence factor and the land subsidence information and handle longer time-series information in order to realize the long time-series prediction of land subsidence. Therefore, a deep learning network architecture of LSTM-TCN was proposed in order to achieve a high-precision subsidence prediction, which combines the respective advantages of the LSTM and TCN models in this study.
In summary, this study considered Henan Province as a case study and aimed to investigate the application of SAR data in the monitoring and prediction of large-scale land subsidence. First, the Sentinel-1A image data of Henan Province from November 2019 to February 2022 were processed using SBAS-InSAR technology to obtain land deformation data. Then, we verified the accuracy of the SBAS-InSAR results by comparing them with the level data. Next, the spatiotemporal characteristics of land deformation were analyzed. We used the LSTM model to extract the time-series information of land subsidence and the TCN model to extract the feature information between the influencing factors and land subsidence. Finally, a deep learning approach based on the LSTM-TCN model was proposed in order to simulate the complex nonlinear relationship between time-series land subsidence data and related environmental factors to predict land subsidence.

2. Materials and Process

2.1. Study Area

Henan Province is located in the southern part of the North China Plain. Henan Province has a total area of 167,000 km2, with a geographical range of 31°23′N–36°22′N and 110°21′E–116°39′E (Figure 1). The central and eastern areas comprise the Huang-Huai-Hai alluvial plains, while the southwest is composed of the Nanyang Basin. Most of this study area is in a warm temperate zone, while the southern part is in a subtropical zone. The average annual precipitation ranges from 407 mm to 1295 mm (from 2017 to 2022) and is concentrated between June and August.

2.2. Sentinel-1A SAR Data

The Sentinel-1A satellite is an environmental monitoring satellite with a revisit period of 12 days. Its ground range resolution is 5 m, and its azimuth resolution is 20 m (single-look). This satellite has a fairly steep incident angle of approximately 34° with the normal orbit. The satellite’s angle of incidence ranges from 30.8° to 46.2°, providing considerable accuracy for SAR interferometric acquisition.
We collected 364 available Sentinel-1A images in the ascending track to cover the whole study area from November 2019 to February 2022 in order to derive the high-frequency land deformation results. The average time interval of these images was approximately one month. In addition, corresponding precision orbit files and 30 m resolution digital elevation model (DEM) data from the Space Shuttle Ordnance Survey Mission (SRTM) were used to correct the data in order to eliminate the effects of systematic errors and topographic residuals. Table 1 and Figure 2 show the details of the Sentinel-1A dataset used in this study.

2.3. Level Measurement Data

A total of 20 high-precision level measurement data points were used to evaluate the accuracy of the deformation results based on InSAR technology. The measurement period was from November 2019 to February 2022. The geographic locations of these monitoring points are indicated by red diamonds in Figure 1.

2.4. Environmental Influencing Factor Data

In this study, apart from the historical sequence of land deformation information, the volumetric soil water layer and monthly precipitation data as the environmental factor were used as the input variables for the model in order to predict the cumulative land subsidence. The volumetric soil water is related to the soil texture, soil depth, and underlying groundwater level, which was derived from ERA5 (fifth generation ECMWF reanalysis) monthly meteorological reanalysis data (ERA5 monthly average data from 1959 to present, available at https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-single-levels-monthly-means?tab=form, accessed on 20 March 2022). Its parameters represented the moisture content in soil layers 1, 2, 3, and 4 (surface is at 0 cm). Volumetric soil water data were obtained by using the monthly average values, and they were distributed over a 0.25° × 0.25° grid. Information on the four-layer volumetric soil aquifers used in this study is shown in Table 2. In this study, the volumetric soil water content data of ERA5 were preprocessed by downscaling [49] and then reprocessed by ordinary kriging linear interpolation in order to obtain data consistent with the InSAR spatial resolution. We obtained the corresponding precipitation data from the China Meteorological Administration. Furthermore, the monthly precipitation data were obtained by accumulating the daily precipitation data and then through global interpolation based on the ordinary kriging linear interpolation method for consistency with the spatial resolution of the InSAR monitoring data.

3. Methods

For land subsidence monitoring, we used SBAS-InSAR technology to obtain the Henan time-series cumulative subsidence data and annual average subsidence rate data. Specifically, we completed InSAR data processing through the InSAR scientific computing environment (ISCE) [50] and the open-source toolkit Miami InSAR Time-series Software in Python (MintPy) [51]. For land subsidence prediction, we used the volumetric soil water layer and precipitation data as environmental factor inputs combined with historical sequence land deformation information to predict land subsidence. Specifically, we first divided all of the time-series subsidence points and established the input matrix. Next, the three models of LSTM, TCN, and LSTM-TCN were constructed to predict the cumulative land subsidence in Henan Province. The overall flow diagram is shown in Figure 3.

3.1. SBAS-InSAR

The SBAS-InSAR method is a distributed scatterer (DS) method, and its main steps are described below. Assume that the N +1 SAR images of the study area were obtained at times t 0 , t 1 , , and t N . We selected the SAR image acquired in November 2019 as the reference image, and the remaining images were uniformly registered to the reference image. After setting the appropriate spatiotemporal baseline, M interferograms were generated. For the j -th interferogram generated by the combination of t a and t b ( t a < t b ), after removing the flat earth effect and terrain phase, the interference phase of any point in interferogram j can be represented as follows by Formula (1):
δ φ j = φ t b φ t a 4 π λ d t b d t a + Δ φ t o p , j + Δ φ a t m , j + Δ φ n o i s e , j
where λ represents the radar central wavelength; d t b and d t a represent the cumulative deformation corresponding to the line of sight (LOS) at t b and t a , respectively, with t 0 as the initial time; and Δ φ t o p , j , Δ φ a t m , j , and Δ φ n o i s e , j represent the residual terrain phase difference, atmospheric phase difference, and noise phase difference, respectively.
After removing the phase other than the deformation, the interferometric phase was simplified and represented by Equations (2) to (3):
Δ δ φ j = φ t b φ t a 4 π λ d t b d t a
d t b d t a = υ i t b t a
where υ i represents the average deformation rate in the direction of the radar LOS from t a to t b .
We used the differential interferogram for phase unrolling to generate the corresponding matrix, and then, we used the singular value decomposition (SVD) method to obtain the final land deformation information.
Figure 1 shows that the study area occupies four different SAR tracks: track 11, track 113, track 40, and track 142. We solved the problem of fusing SAR image data from different locations of the same orbit and different orbits by selecting common reference points in overlapping regions in this study. The black dots and green triangular dots in Figure 1 indicate the locations of the common reference points for the same and different orbits, respectively. The above common reference points are stable monitoring points with no displacement or deformation during the entire monitoring period of this study.
During the single-track SAR data processing, the black dots were used as reference points for the same track land deformation calculation to complete the data processing. Then, we compared the time-series displacement deformation at the level measurement points in the study area with the monitoring results and calculated the offset to calibrate the single-track InSAR monitoring results. For multitrack SAR data, the fusion of multitrack SAR data was essential because of the differences in monitoring results due to the different imaging angles of SAR images from different tracks, the capture time, and the position of reference points during data processing. In this study, we calibrated the orbit monitoring results by using the difference between the monitoring results of two different orbit overlapping areas as the offset in order to ensure the accuracy of the fusion results. Based on the single-track data processing, we used the green triangle point as the common reference point to obtain the second monitoring results again experimentally. Specifically, we took the calibrated track 40 deformation point as the reference. Then, we used the green triangular point in the overlapping area of track 40 and track 142 as the reference point to obtain the second monitoring results of track 40 and track 142. Next, we calibrated the first single-track monitoring results of track 142 based on the difference between the two experimental results for the track 40 and track 142 overlap area as the offset (Figure 4). Specifically, we took the SAR processing result of 142 for the first single track plus the offset as the result after fusion. Then, the results were stitched with the results of track 40 to complete the fusion processing of the SAR results of track 40 and track 142. Following the above method, the difference between the experimental results of the track 40 and track 113 overlapping area was calculated twice as the offset to calibrate the monitoring results of track 113 for the first single track. We selected the calibrated single-track 40 deformation point as the reference. Subsequently, we completed the data fusion of track 40, track 142, and track 113. Notably, when fusing the track 11 data, we took the results of the fused data of track 40, track 142, and track 113 as the reference. The difference between the second monitoring result of the overlapping area of track 11 and track 113 and the result of the overlapping area after the fusion of the track 113, track 142, and track 40 was calculated to complete the calibration of the first single-track monitoring result of track 11. Finally, we completed the data fusion of track 40, track 142, track 113, and track 11. In addition, we compared the final fusion results with the level measurement data for verification to ensure the accuracy of the InSAR monitoring results in the study area. InSAR monitoring data used in this study were obtained based on the World Geodetic System 1984 (WGS1984). All leveling data were obtained based on the China Geodetic Coordinate System 2000 (CGCS 2000). We used the specific parameter combination in the research area for the coordinate transformation and completed the data accuracy assessment in the same coordinate system.
In this study, we used ISCE to complete all SAR dataset registration, interference, and unwrapping in a Conda virtual environment, and 1014 interferograms were generated. Then, we used the open-source toolkit MintPy to perform a series of processing steps on the dataset in order to obtain the land deformation data. First, we registered and corrected the images with precise orbit files. Second, we removed the flat ground effect and terrain phase estimation with the DEM data of the SRTM satellite at a 30 m resolution. Next, the time baseline and spatial baseline were calculated, the coherence of pixels was calculated, and the points with high coherence were selected for expansion to generate an interferogram. Finally, we used a set of registered and expanded interferograms generated in the first step, and the land deformation data were generated. We used global atmospheric models for tropospheric delay correction in this study. Specifically, we processed the time-series data by introducing online ERA5 meteorological factor data. We also used the root mean square to evaluate various noises in the test process and obtained experimental results meeting the accuracy requirements. During the experiment, the method of automatically selecting multiple reference master images was used to calculate land deformation. Finally, the first image was used as a reference to obtain the final land deformation data.

3.2. LSTM-TCN Model

The LSTM model is a common prediction model for land subsidence, which has a strong ability to mine the nonlinear characteristics of geographical phenomena along with model time evolution and, thus, it has been extensively used in geosciences, including land subsidence, hydrology, and ecology [52]. The basic unit of the LSTM is composed of three gates (Figure 5). The LSTM models discard unimportant information and retain key information through gate control, which improves the studying ability of the models [53].
The main operations of this method are shown in Equations (4)–(8).
f t = σ W f h t 1 , x t + b f
i t = σ W i h t 1 , x t + b i
C t = f t C t 1 + i t tanh W c h t 1 , x t + b c
o t = σ W o h t 1 , x t + b o
h t = o t tanh C t
where f t represents the forget gate, σ and tan h represent the activation function, W f , W i , W c , and W o represent the input component of the weight matrix, h t 1 and C t 1 represent the hidden state value at the previous t 1 moment, x t represents the input value, i t represents the input unit and C t represents x t filtered values. b f , b i , b c , and b o represent the offset at time t , o t represents the output gate unit, and h t represents the hidden state value at the previous time t .
The operation between the dilated convolutional layer of the TCN model and the input sequence is given as follows:
F s = X d f s = i = 0 k 1 f i X s d * i
where the input sequence is X n , the filter is f : 0 ,   , k 1 , d is the size of the hole factor, k is the size of the filter, and s d i represents the previous neuron of the hole convolution layer. The meaning of the hole convolution is to control the interval size of each convolution layer neuron to increase or reduce the receptive field. When d is larger, the hole convolution layer has a larger receptive field, which can obtain more effective input sequence information.
From the above LSTM and TCN model architecture, the LSTM serialization calculation structure has certain limitations in the land subsidence prediction model. The TCN model can parallelize the input sequence to compensate for the lack of an LSTM model structure. Therefore, we propose an LSTM-TCN model that uses the extended convolution property of the TCN model to capture more critical information. Specifically, the LSTM-TCN model in this study combines the characteristics of the LSTM model to extract temporal information, the TCN convolutional layer to extract the characteristic information of influence factors and land subsidence, and the multidimensional residual connection layer to improve the accuracy of the experimental results. To better represent the input matrix, we divided the cumulative subsidence matrix and the corresponding environmental factor matrix. The cumulative subsidence of a single point m from November 2019 to February 2022 is shown in the following matrix L m .
L m = L m , 1 L m , 2 L m , 26 L m , 27
The time-series environmental factor data for a single point m from November 2019 to February 2022 are shown in the following matrix F m , where P m , 1 , , and P m , 27 , represent the monthly precipitation of point m from November 2019 to February 2022, V 1 m , 1 , , V 1 m , 27 , V 2 m , 1 , , V 2 m , 27 , V 3 m , 1 , , V 3 m , 27 , V 4 m , 1 , , and V 4 m , 27 represent the volumetric water content of the soil aquifer at point m in layers 1, 2, 3, and 4 from November 2019 to February 2022, respectively.
F m = P m , 1 V 1 m , 1 V 2 m , 1 V 3 m , 1 V 4 m , 1 P m , 2 V 1 m , 2 V 2 m , 2 V 3 m , 2 V 4 m , 2 P m , 26 V 1 m , 26 V 2 m , 26 V 3 m , 26 V 4 m , 26 P m , 27 V 1 m , 27 V 2 m , 27 V 3 m , 27 V 4 m , 27
To explore the influence of monthly time-series land subsidence on the experimental results, the above two matrices were combined to predict the cumulative land subsidence of time t + 1 according to the existing information and historical data obtained from time t . The input matrix X t , m of any subsidence point m is obtained as follows:
X t , m = L t , m P t , m V 1 t , m V 2 t , m V 3 t , m V 4 t , m L t 1 , m P t 1 , m V 1 t 1 , m V 2 t 1 , m V 3 t 1 , m V 4 t 1 , m L t d + 1 , m P t d + 1 , m V 1 t d + 1 , m V 2 t d + 1 , m V 3 t d + 1 , m V 4 t d + 1 , m
where the length d of the matrix represents the length used to control the historical time series, ranging from 0 to t .
The output matrix is Y t , m = L t + 1 , m .
We designed three models of LSTM, TCN, and LSTM-TCN for an accuracy comparison to verify the accuracy of the three models’ results. The architectures of the three models are shown in Figure 6.
The historical sequence of environmental factors and cumulative subsidence data of all subsidence points at time t was used as the two-dimensional input matrix, and the land-cumulative subsidence data at time t + 1 were used as output, as shown in Figure 6. The LSTM model used in the experiment is shown in Figure 6a, which was composed of three layers of LSTM hidden layers and a dense layer. The hidden layer of the TCN model is represented in Figure 6b, composed of four dilated convolution layers. Each layer’s sampling interval was increased exponentially by two, aiming to obtain a larger receptive field. A LSTM and TCN fusion model with a hidden layer consisting of a 1-layer LSTM model and 3-layer dilated convolutional layers is shown in Figure 6c, which combines the advantages of the above two models to predict cumulative land subsidence. The dilated convolution layers were connected by the residual connection module in Figure 6b,c, which obtained effective historical sequence information through feature extraction and normalized the information.

4. Results

4.1. Seasonal Characteristics of Land Subsidence and Its Causal Analysis

All subsidence points with annual average subsidence rates greater than −10 mm/a in the InSAR monitoring results were selected as rapid subsidence monitoring points in the study area for experiments to analyze the seasonal characteristics of land subsidence in Henan Province and to explore its causes in the context of environmental factors. Specifically, we calculated the average deformation rate of the land in the subsidence region in spring–summer and autumn–winter to explore the characteristics of seasonal land subsidence in the subsidence region. The sum of the cumulative spring–summer deformation variables of the subsidence area land in 2020 and 2021 was divided by two as the semiannual average spring–summer deformation rate of the subsidence area land, and the sum of the cumulative autumn–winter deformation variables was divided by two as the semiannual average autumn–winter deformation rate of the subsidence area land. If the deformation rate was positive, the land of the region was uplifting; otherwise, the land was subsiding. A total of 15,364 subsidence points were selected as characteristic points to explore the seasonal characteristics of land subsidence in Henan Province. Figure 7 and Figure 8 show the results of the seasonal characteristics of land subsidence. Figure 7 shows the semiannual average deformation rates in spring and summer in Henan Province. Anyang, Hebi, Shangqiu, Zhengzhou, Xuchang, Pingdingshan, Nanyang, Xinyang, Kaifeng, Luoyang and Sanmenxia all experienced subsidence in spring and summer in 2020 and 2021, as shown in Figure 7. Among them, the average subsidence rate of Luoyang exceeded 100 mm/semiannual. The average subsidence rates of Sanmenxia and Pingdingshan ranged from 80 mm/semiannual to 100 mm/semiannual. Figure 8 shows the semiannual average deformation rate of the land in Henan Province in autumn and winter. Anyang, Hebi, Shangqiu, Xuchang, Pingdingshan, Nanyang and Xinyang showed subsidence of the land in the autumn and winter of 2020 and 2021. The average subsidence rates in these regions were smaller than the subsidence rates in the corresponding regions in Figure 7. Some areas in Sanmenxia, Luoyang, Zhengzhou and Kaifeng experienced land uplift in the autumn and winter of 2020 and 2021. The average rate of land uplift in some areas of Luoyang exceeded 20 mm/semiannual. We also calculated the percentage of subsidence points in the study area where subsidence occurred in the autumn and winter of 2020 and 2021, which reached 93.34%. In general, the land in the subsidence area in the study area appeared to sink in spring and summer, the average rate of subsidence in the corresponding area was greater than that in autumn and winter, and the land ratio of some of the subsidence areas rose in autumn and winter compared with that in spring and summer, delaying the land subsidence phenomenon.
In this study, environmental factor data and subsidence point data of subsidence areas were combined to analyze the causes of seasonal characteristics of land subsidence in Henan Province, as shown in Table 3. In general, the semiannual average subsidence rate in the subsidence area increased by 1.19 mm/semiannual in spring–summer compared to autumn–winter. The semiannual average precipitation in the subsidence area reached 619.29 mm/semiannual in spring–summer and 229.61 mm/semiannual in autumn–winter, as shown in Table 3. The semiannual average precipitation in the subsidence area decreased significantly in autumn–winter compared to spring–summer. Combining the semiannual average subsidence rate data of spring–summer and autumn–winter in the subsidence region, the results showed that the influence of precipitation on land subsidence in the subsidence region has a certain hysteresis, which is manifested in the fact that abundant precipitation in spring and summer can replenish groundwater and delay land subsidence in autumn and winter. Therefore, in general, the subsidence phenomenon is more obvious in spring and summer in the subsidence area than in autumn and winter. In addition, we also found that the volumetric soil layer water content directly affects land subsidence. Specifically, the decrease in water content in the volumetric soil layer resulted in the subsidence of the land, and the increase in water content in the volumetric soil layer resulted in the delayed subsidence or uplift of the land. The semiannual average water content of the four volumetric soil layers in this study was significantly smaller in spring and summer than in autumn and winter, explaining the phenomenon that the semiannual average subsidence rate in the subsidence area was larger in spring and summer than in autumn and winter.

4.2. Land Subsidence Monitoring Results

To ensure the availability of InSAR data, we compared the annual mean deformation rate results between InSAR and leveling data and used the correlation coefficient graph to show the relationship between the two. According to Figure 9, the coefficient of determination ( R 2 ) of the two was 0.98, and the root mean square error (RMSE) value was 4.45 mm/a, which verified the availability of the monitoring data.
The cumulative land deformation data and annual average land deformation rate data are shown in Figure 10 and Figure 11 for the study area during the monitoring period, respectively. According to Figure 10, the land deformation of Anyang, Hebi, Zhengzhou, and Luoyang in the study area was relatively obvious. Among them, the land subsidence in Anyang was the most obvious, and the subsidence rate was more than 80 mm/a. In addition, the land uplift in Hebi, Zhengzhou, and Shangqiu was more than 20 mm/a. Figure 10 shows that there was an obvious subsidence center in Anyang, and the accumulated subsidence exceeded 150 mm. Adjacent to Anyang, in Hebi, the land of apparent uplift, the cumulative uplift reached 50 mm. Moreover, the accumulated subsidence in Zhengzhou, Kaifeng, Pingdingshan, and Xinxiang reached 150 mm, while the accumulated subsidence in Sanmenxia, Nanyang, Zhumadian, Xinyang, and Zhoukou reached 100 mm. According to the InSAR monitoring results, the accumulated land subsidence in Shangqiu city during the monitoring period was the largest, reaching −194.78 mm.
From the perspective of geographic location and topography, Anyang city has a mountainous area in the west, and the east is a plain. According to Figure 11, the western part of Anyang showed an upward trend, and the eastern part showed a downward trend, forming a subsidence funnel in Anyang County in the north. The Weihe River passes through Xunxian County in Hebi. There was a small subsidence funnel center in Xunxian County, which showed a clear upward trend around the center. The general topography of Zhengzhou is high in the southwest and low in the northeast, with a stepped descent. Figure 11 illustrates that the land of the northern part of Zhengzhou was uplifted, and the eastern part was decreased. Luoyang is in western Henan Province. The mountains and hills are connected, and the terrain is complex. There were two obvious subsidence centers in the eastern part of the city, and the cumulative subsidence exceeded 100 mm.

4.3. Land Subsidence Model Prediction Results

We tested and verified the input matrix of different monthly time series and found that the model accuracy was the highest when the length of input data was seventeen months. Therefore, the results of this section were based on the best performance over seventeen months. The influence of different input lengths on the accuracy of the prediction model is discussed in Section 5.3. A total of 11,189,314 deformation monitoring points were obtained, and the average density was 67/km2. For better monitoring and prediction of subsidence areas, 15,364 characteristic points with significant subsidence were selected for training and prediction of three models with an annual average subsidence rate exceeding −10 mm/a, and their geographical locations are shown in Figure 12. A total of 414,828 records for 15,364 subsidence points were used for model training, verification, and test evaluation. The input data were first normalized, and then, all the datasets were disrupted with random seeds. Notably, 70% of all datasets were used for training, 20% of all datasets were used to verify the trained model, and the remaining 10% of all datasets were used as a test set to evaluate the model.
We set different parameters and experimentally determined the parameter combinations for the optimal parameters. The loss curves of the experimental training and validation sets are shown in Figure 13, and the optimal parameter combination for the LSTM-TCN model is a learning rate (lr) of 0.0003, one LSTM layer, three dilated convolutional layers, 100 epochs, and 64 batch sizes. In addition, the optimal parameter combination for the LSTM model is a lr of 0.0003, three LSTM layers, 100 epochs, and 64 batch sizes. The optimal combination of parameters for the TCN model is a lr of 0.0003, 4 dilated convolutional layers, 100 epochs, and 64 batch sizes. We also observe that the training set loss for the four parameter combinations of the LSTM-TCN model is smaller than that of the other two models, and its convergence is faster than that of the LSTM and TCN models, as shown in Figure 13.
We constructed the experimental LSTM, TCN and LSTM-TCN models using the optimal combination of parameters of L4, T4, and LT4 in Table 4. We obtained the correlation results between the predicted data of the LSTM, TCN, and LSTM-TCN models and the SBAS-InSAR result data for a more visual observation of the experimental accuracy, as shown in Figure 14. We used the RMSE, mean absolute error (MAE), mean absolute percentage error (MAPE), prediction error less than 3 mm, 0~100 mm MAE (MAE values of InSAR monitoring results and model prediction results in the range of 0~100 mm), 100~200 mm MAE (MAE values of InSAR monitoring results and model prediction results in the range of 100~200 mm), and R 2 to evaluate the model’s accuracy. From the accuracy results, the RMSE of the LSTM-TCN model was reduced by 74.20% and 54.56% compared to the LSTM and TCN models, respectively. The results in Table 5 show that the prediction accuracy of the three models for subsidence results in the range of 0~100 mm is higher than that in the range of 100~200 mm. In addition, since the LSTM-TCN model can be processed in parallel, its training time is reduced by 40% and memory usage is reduced by 20% compared to the LSTM model.
Four significant subsidence points were selected to further verify the advantages of the proposed prediction models for the long time-series nonlinear variation in the subsidence centers. Figure 15 and Figure 16 show the geographical locations of the subsidence centers in Anyang, Hebi, Zhengzhou, and Luoyang, respectively, and the results of the fitting of the three models to the time-series subsidence variation in the subsidence points. According to Figure 15, all four validation points are located in the center of the subsidence funnel, which can reflect the future deformation trend of this subsidence center. Figure 16 shows that the LSTM-TCN model has better fitting results for the future deformation trend of the subsidence center compared with the LSTM model and the TCN model. Specifically, the LSTM-TCN model has a significant advantage over the LSTM and TCN models in predicting changes in D with lifting and sinking trends in the time-series deformation. The prediction accuracy R 2 of the LSTM-TCN model reached 0.78 in the fitting of the temporal variation at point D, which improved by 0.63 and 0.38 compared to the LSTM and TCN models, respectively. We also find that the LSTM model and TCN model have better prediction results for point A with a long-term subsidence trend, and their prediction accuracy RMSEs are 5.75 mm and 2.48 mm, respectively. However, the LSTM-TCN model achieves a prediction accuracy RMSE of 2.28 mm for the temporal variation in point A, which has better fitting results compared with the LSTM model and TCN model. For points B and C, where the rising and sinking phenomena occur frequently in the time-series deformation, the LSTM and TCN models have poor prediction results compared to point A. Their prediction accuracy R 2 is less than 0.6, while the prediction accuracy R 2 of the LSTM-TCN model is greater than 0.8. In general, the LSTM-TCN has better prediction results in predicting the time-series changes in subsidence points with long-term subsidence phenomena, indicating the prediction potential of the model in the future changes in deformation trends in subsidence centers.

5. Discussion

5.1. Influence of the Sample Division Ratio and Quantity on the Prediction Results

All sample points were divided and tested at 10%, 20%, 30%, 40%, 50%, 60%, 70%, and 80% to verify the sensitivity of the LSTM model, TCN model, and LSTM-TCN model to the number of subsidence points, and the results are shown in Figure 17. The results in Figure 17 show that the fluctuation range of the sample size change curve of the LSTM-TCN model is smaller than that of the other two models, indicating that the model is less sensitive to the number of subsidence points. When the sample size increases, the LSTM model prediction accuracy RMSE and MAE decrease, R 2 continues to increase, and the MAPE value decreases and then increases. The above results indicate that the LSTM model is more sensitive to the sample size and requires a sufficiently large sample size to improve the prediction accuracy. The TCN model has smaller fluctuations in the curve changes in RMSE, R 2 , MAE, and MAPE with increasing sample size prediction accuracy compared with the LSTM model, and the results of its MAPE curve changes are even close to those of the LSTM-TCN model. Collectively, the LSTM-TCN model can also obtain better prediction results when the sample size is small. When the sample size increases to the optimal threshold, the LSTM-TCN model can combine the advantages of the LSTM model and the TCN model to improve the accuracy of the prediction results. The prediction accuracy curve of the LSTM-TCN model remains smooth when the sample size reaches the optimal threshold, showing a low sensitivity to sample size.

5.2. Advantages of Adding Environmental Factors

The experimental results of the LSTM, TCN, and LSTM-TCN models without and with environmental factors are shown in Table 6. The results show that the accuracy of RMSE, R 2 , MAE, MAPE, and the prediction error is less than 3 mm, and the three models are improved after adding environmental factors. The R 2 of the LSTM model increased from 0.84 to 0.89, the R 2 of the TCN model increased from 0.94 to 0.96, and the R 2 of the LSTM-TCN model increased from 0.95 to 0.99, which greatly improved the accuracy. From the overall results, regardless of whether environmental factors were added, the accuracy ranking of the three models was LSTM-TCN > TCN > LSTM, which further demonstrated the effectiveness of the constructed LSTM-TCN model. After adding environmental factors, the LSTM-TCN model has more obvious accuracy improvement compared with the LSTM model and TCN model, which shows the superiority of the model for feature extraction between land subsidence and environmental factors.

5.3. Influence of the Different Lengths on the Prediction Results

In this study, we investigated the effect of different monthly time-series input matrices on the experimental results, which are shown in Figure 18. The results in Figure 18 show that the best experimental results were obtained with a two-dimensional matrix consisting of environmental factors and cumulative land subsidence for 17 consecutive months. In this study, two-dimensional input matrices, larger or smaller than 17 months for the 27-month time series of land-cumulative subsidence data, would reduce the accuracy results of the LSTM-TCN model. Choosing the appropriate two-dimensional matrix for the length of the monthly time series allows the model to fully extract the characteristic information between the land-cumulative subsidence data and the environmental factors and then to better fit the future deformation trend of the land subsidence, which provides us with a reference for timely monitoring and early warning of the occurrence of geological hazards.

5.4. Exploration of the LSTM-TCN-Based Model in Time-Series Land Subsidence Prediction

In the prediction of time-series land subsidence, researchers used traditional machine learning algorithms to model and predict land subsidence from a data-driven perspective. For example, the boosted regression trees and extreme gradient boosting algorithms model and predict land subsidence from a data-driven perspective, and this study achieved the modeling and prediction of land subsidence rates in agricultural areas [54]. There are also researchers who identified landslide-prone areas by mapping them with the help of machine learning models, which served to prevent and mitigate landslides. Specifically, this study used the logistic regression (LR) model and the random forest (RF) model to evaluate the sensitivity of landslides [55]. In addition, many researchers used LSTM models for modeling predictions of time-series land subsidence. The literature [56] compared the prediction results of the LSTM model with the recurrent neural network (RNN) model and a multi-layer perceptron model. The results showed that the LSTM model was more efficient in prediction than both models and can achieve effective short-term prediction of land subsidence. The literature [57] used a convolutional neural network to extract spatial feature information of land subsidence in mining areas and the LSTM model to extract temporal feature information of land subsidence in order to achieve the spatial and temporal prediction of land subsidence. The literature [58] divided the study area into subregions based on InSAR monitoring results and used LSTM models to capture the non-linear correlation characteristics of each subregion in order to achieve prediction of land subsidence. The literature [59] first used the SBAS-InSAR technique to extract land deformation information from the mining areas and then used a combination of LSTM models and grid search algorithms to predict the time-series land subsidence. However, the abovementioned traditional machine learning methods cannot extract and model the cumulative land subsidence information in a time series and cannot obtain comprehensive information on future land subsidence trends. Deep learning models, such as RNN and LSTM models, suffer from gradient explosion and gradient disappearance, which have limitations in modeling long time-series land subsidence prediction.
Compared with the above machine learning and deep learning models, the LSTM-TCN in this study relies on its scalable perceptual field to extract feature information on long sequences before and after, achieving the retention of longer sequence information and obtaining better prediction results. In addition, the LSTM-TCN model can also extract time-series feature information on influence factors and land subsidence, thus enabling the prediction of the time-series land subsidence by combining multiple influence factors. Notably, due to the parallelized sequence structure of the TCN model, the LSTM-TCN model can handle large amounts of data from a larger study area more efficiently. The LSTM-TCN model has promising applications in large-scale regional land subsidence in Henan Province.

5.5. Effect of Input Environmental Factors

Most researchers have used groundwater level data for modeling predictions of land subsidence, which was easy to achieve on a small scale. However, for large monitoring areas, groundwater level data for the entire area may be unavailable. In this study, in the absence of groundwater level data in part of the study area, soil volumetric water content data was used as an input factor in the prediction model in order to complete the prediction of the time-series land subsidence for the entire study area. We collected and collated data on groundwater levels in Zhengzhou City, as shown in Table 7. To assess the feasibility of the method in this study, we used groundwater levels and precipitation data, volumetric soil water content and precipitation data of Zhengzhou City as input factors for the LSTM-TCN model, respectively, and compared their prediction results, as shown in Table 8. Table 8 shows that the prediction accuracy of land subsidence was slightly better than that of volumetric soil water content and precipitation data when groundwater level and precipitation data were used as input factors for the LSTM-TCN model, indicating the importance of groundwater level data on land subsidence. However, in general, for the whole study area, the modelling experiments using volumetric soil water content and precipitation data for land subsidence prediction obtained better prediction results over a large study area, reflecting the feasibility of the experimental design of this study.

6. Conclusions

The Sentinel-1A multitrack image was processed by SBAS-InSAR technology to obtain the land deformation data of Henan Province in this study. Land subsidence was predicted by modeling land subsidence and environmental factor data. The conclusions are described as follows.
(1)
Henan Province has undergone significant deformation in the past two years. From November 2019 to February 2022, the maximum and minimum land deformation rates in Henan Province were −94.54 mm/a and 41.23 mm/a, respectively. The most severe land subsidence occurred in Shangqiu city, and the maximum cumulative subsidence was −194.78 mm.
(2)
Land subsidence in the study area has obvious seasonal characteristics. The semiannual average rate of land subsidence in the study area is greater in spring and summer than in autumn and winter. The effect of precipitation on land subsidence has a certain hysteresis, and abundant precipitation in spring and summer can replenish the lower groundwater and delay land subsidence in autumn and winter. The volumetric soil layer water content shows the same trend as land subsidence.
(3)
The LSTM-TCN model can combine the advantages of the LSTM and TCN models, which is better for exploring the nonlinear characteristics of land subsidence. Specifically, the LSTM-TCN model could effectively predict the cumulative land subsidence in Henan Province, and its prediction accuracy was better than that of the LSTM and TCN models in areas with obvious subsidence. When predicting the area with subsidence exceeding 100 mm, the MAE of the LSTM-TCN model was 1.97 mm, which was 68.35% and 63.98% higher than the 6.23 mm value of the LSTM model and the 5.48 mm value of the TCN model, respectively.
(4)
After introducing the environmental impact factors in the prediction model, the determination coefficient R 2 of the LSTM-TCN model increased from 0.95 to 0.99, indicating that environmental factors and geological characteristics are indispensable in exploring land subsidence.

Author Contributions

Conceptualization, H.G. and Y.Y.; methodology, J.W.; validation, R.Z., Q.C., J.L. and W.D.; formal analysis, B.Q.; investigation and resources and data curation, J.C.; writing—original draft preparation, H.G. and Y.Y.; writing—review and editing, D.Z.; visualization, H.B.; supervision and project administration, S.Z.; funding acquisition, H.G. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by Major Science and Technology Special Projects in Henan Province (221100210600); Major Science and Technology Special Projects in Henan Province (201400211000); Major Science and Technology Special Projects in Henan Province (201400210100); and Science and Technology Tackling Plan of Henan Province (222102320220).

Data Availability Statement

The data presented in this study are available upon request from the corresponding author.

Acknowledgments

The authors would like to thank the anonymous reviewers for their valuable feedback on the manuscript. We would also like to thank the ESA (European Space Agency) for providing the Sentinel-1A SAR dataset for this study, the ERA5 and China Meteorological Administration for providing environmental factor data, and the ISCE (2.5.3) and MintPy (1.3.3) development teams for providing the opensource code.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Carlson, G.; Shirzaei, M.; Werth, S.; Zhai, G.; Ojha, C. Seasonal and Long-Term Groundwater Unloading in the Central Valley Modifies Crustal Stress. J. Geophys. Res. Solid Earth 2020, 125, e2019JB018490. [Google Scholar] [CrossRef] [PubMed]
  2. Cianflone, G.; Vespasiano, G.; Tolomei, C.; De Rosa, R.; Dominici, R.; Apollaro, C.; Walraevens, K.; Polemio, M. Different Ground Subsidence Contributions Revealed by Integrated Discussion of Sentinel-1 Datasets, Well Discharge, Stratigraphical, and Geomorphological Data: The Case of the Gioia Tauro Coastal Plain (Southern Italy). Sustainability 2022, 14, 2926. [Google Scholar] [CrossRef]
  3. Cigna, F.; Tapete, D. Satellite InSAR survey of structurally-controlled land subsidence due to groundwater exploitation in the Aguascalientes Valley, Mexico. Remote Sens. Environ. 2021, 254, 112254. [Google Scholar] [CrossRef]
  4. Cigna, F.; Tapete, D. Urban growth and land subsidence: Multi-decadal investigation using human settlement data and satellite InSAR in Morelia, Mexico. Sci. Total Environ. 2022, 811, 152211. [Google Scholar] [CrossRef] [PubMed]
  5. Kamali Maskooni, E.; Naghibi, S.A.; Hashemi, H.; Berndtsson, R. Application of Advanced Machine Learning Algorithms to Assess Groundwater Potential Using Remote Sensing-Derived Data. Remote Sens. 2020, 12, 42. [Google Scholar] [CrossRef]
  6. Nguyen, M.; Lin, Y.N.; Tran, Q.C.; Ni, C.-F.; Chan, Y.-C.; Tseng, K.-H.; Chang, C.-P. Assessment of long-term ground subsidence and groundwater depletion in Hanoi, Vietnam. Eng. Geol. 2022, 299, 106555. [Google Scholar] [CrossRef]
  7. Guzy, A.; Malinowska, A. State of the Art and Recent Advancements in the Modelling of Land Subsidence Induced by Groundwater Withdrawal. Water 2020, 12, 2051. [Google Scholar] [CrossRef]
  8. van Natijne, A.L.; Bogaard, T.A.; van Leijen, F.J.; Hanssen, R.F.; Lindenbergh, R.C. World-wide InSAR sensitivity index for landslide deformation tracking. Int. J. Appl. Earth Obs. Geoinf. 2022, 111, 102829. [Google Scholar] [CrossRef]
  9. Cigna, F.; Tapete, D. Present-day land subsidence rates, surface faulting hazard and risk in Mexico City with 2014–2020 Sentinel-1 IW InSAR. Remote Sens. Environ. 2021, 253, 112161. [Google Scholar] [CrossRef]
  10. Peng, M.; Lu, Z.; Zhao, C.; Motagh, M.; Bai, L.; Conway, B.D.; Chen, H. Mapping land subsidence and aquifer system properties of the Willcox Basin, Arizona, from InSAR observations and independent component analysis. Remote Sens. Environ. 2022, 271, 112894. [Google Scholar] [CrossRef]
  11. Rezaei, A.; Mousavi, Z. Characterization of land deformation, hydraulic head, and aquifer properties of the Gorgan confined aquifer, Iran, from InSAR observations. J. Hydrol. 2019, 579, 124196. [Google Scholar] [CrossRef]
  12. Hu, L.; Navarro-Hernández, M.I.; Liu, X.; Tomás, R.; Tang, X.; Bru, G.; Ezquerro, P.; Zhang, Q. Analysis of regional large-gradient land subsidence in the Alto Guadalentín Basin (Spain) using open-access aerial LiDAR datasets. Remote Sens. Environ. 2022, 280, 113218. [Google Scholar] [CrossRef]
  13. Kumar Maurya, V.; Dwivedi, R.; Ranjan Martha, T. Site scale landslide deformation and strain analysis using MT-InSAR and GNSS approach–A case study. Adv. Space Res. 2022, 70, 3932–3947. [Google Scholar] [CrossRef]
  14. Shahbazi, S.; Mousavi, Z.; Rezaei, A. Constraints on the hydrogeological properties and land subsidence through GNSS and InSAR measurements and well data in Salmas plain, northwest of Urmia Lake, Iran. Hydrogeol. J. 2022, 30, 533–555. [Google Scholar] [CrossRef]
  15. Dong, J.; Lai, S.; Wang, N.; Wang, Y.; Zhang, L.; Liao, M. Multi-scale deformation monitoring with Sentinel-1 InSAR analyses along the Middle Route of the South-North Water Diversion Project in China. Int. J. Applied Earth Obs. Geoinf. 2021, 100, 102324. [Google Scholar] [CrossRef]
  16. Xie, X.; Xu, C.; Wen, Y.; Li, W. Monitoring Groundwater Storage Changes in the Loess Plateau Using GRACE Satellite Gravity Data, Hydrological Models and Coal Mining Data. Remote Sens. 2018, 10, 605. [Google Scholar] [CrossRef]
  17. Ojha, C.; Werth, S.; Shirzaei, M. Recovery of aquifer-systems in Southwest US following 2012–2015 drought: Evidence from InSAR, GRACE and groundwater level data. J. Hydrol. 2020, 587, 124943. [Google Scholar] [CrossRef]
  18. Vasco, D.W.; Kim, K.H.; Farr, T.G.; Reager, J.T.; Bekaert, D.; Sangha, S.S.; Rutqvist, J.; Beaudoing, H.K. Using Sentinel-1 and GRACE satellite data to monitor the hydrological variations within the Tulare Basin, California. Sci. Rep. 2022, 12, 3867. [Google Scholar] [CrossRef]
  19. Hou, J.; Xu, B.; Li, Z.; Zhu, Y.; Feng, G. Block PS-InSAR ground deformation estimation for large-scale areas based on network adjustment. J. Geod. 2021, 95, 111. [Google Scholar] [CrossRef]
  20. Lyu, M.; Ke, Y.; Guo, L.; Li, X.; Zhu, L.; Gong, H.; Constantinos, C. Change in regional land subsidence in Beijing after south-to-north water diversion project observed using satellite radar interferometry. GIScience Remote Sens. 2019, 57, 140–156. [Google Scholar] [CrossRef]
  21. Shi, X.; Yang, C.; Zhang, L.; Jiang, H.; Liao, M.; Zhang, L.; Liu, X. Mapping and characterizing displacements of active loess slopes along the upstream Yellow River with multi-temporal InSAR datasets. Sci. Total Environ. 2019, 674, 200–210. [Google Scholar] [CrossRef] [PubMed]
  22. Wu, S.; Yang, Z.; Ding, X.; Zhang, B.; Zhang, L.; Lu, Z. Two decades of settlement of Hong Kong International Airport measured with multi-temporal InSAR. Remote Sens. Environ. 2020, 248, 111976. [Google Scholar] [CrossRef]
  23. Zhang, Y.; Meng, X.M.; Dijkstra, T.A.; Jordan, C.J.; Chen, G.; Zeng, R.Q.; Novellino, A. Forecasting the magnitude of potential landslides based on InSAR techniques. Remote Sens. Environ. 2020, 241, 111738. [Google Scholar] [CrossRef]
  24. Goorabi, A.; Karimi, M.; Yamani, M.; Perissin, D. Land subsidence in Isfahan metropolitan and its relationship with geological and geomorphological settings revealed by Sentinel-1A InSAR observations. J. Arid. Environ. 2020, 181, 104238. [Google Scholar] [CrossRef]
  25. Li, J.; Wang, S.; Michel, C.; Russell, H.A.J. Surface deformation observed by InSAR shows connections with water storage change in Southern Ontario. J. Hydrol. Reg. Stud. 2020, 27, 100661. [Google Scholar] [CrossRef]
  26. Saowiang, K.; Giao, P.H. Numerical analysis of subsurface deformation induced by groundwater level changes in the Bangkok aquifer system. Acta Geotech. 2020, 16, 1265–1279. [Google Scholar] [CrossRef]
  27. Sharma, U.; Khan, A.; Dutta, V. Long-term sustainability of groundwater resources in the central Ganga Alluvial Plain, India: Study from Gomti River Basin. Environ. Dev. Sustain. 2021, 23, 16015–16037. [Google Scholar] [CrossRef]
  28. Wang, L.; Deng, K.; Zheng, M. Research on ground deformation monitoring method in mining areas using the probability integral model fusion D-InSAR, sub-band InSAR and offset-tracking. Int. J. Appl. Earth Obs. Geoinf. 2020, 85, 101981. [Google Scholar] [CrossRef]
  29. Fiorentini, N.; Maboudi, M.; Leandri, P.; Losa, M.; Gerke, M. Surface Motion Prediction and Mapping for Road Infrastructures Management by PS-InSAR Measurements and Machine Learning Algorithms. Remote Sens. 2020, 12, 3976. [Google Scholar] [CrossRef]
  30. Hussain, M.A.; Chen, Z.; Shoaib, M.; Shah, S.U.; Khan, J.; Ying, Z. Sentinel-1A for monitoring land subsidence of coastal city of Pakistan using Persistent Scatterers In-SAR technique. Sci. Rep. 2022, 12, 5294. [Google Scholar] [CrossRef]
  31. Zhang, Z.; Zeng, Q.; Jiao, J. Deformations monitoring in complicated-surface areas by adaptive distributed Scatterer InSAR combined with land cover: Taking the Jiaju landslide in Danba, China as an example. ISPRS J. Photogramm. Remote Sens. 2022, 186, 102–122. [Google Scholar] [CrossRef]
  32. Chang, F.; Dong, S.; Yin, H.; Wu, Z. Using the SBAS InSAR technique to monitor surface deformation in the Kuqa fold-thrust belt, Tarim Basin, NW China. J. Asian Earth Sci. 2022, 231, 105212. [Google Scholar] [CrossRef]
  33. Li, S.; Xu, W.; Li, Z. Review of the SBAS InSAR Time-series algorithms, applications, and challenges. Geod. Geodyn. 2022, 13, 114–126. [Google Scholar] [CrossRef]
  34. Li, X.; Yan, L.; Lu, L.; Huang, G.; Zhao, Z.; Lu, Z. Adjacent-Track InSAR Processing for Large-Scale Land Subsidence Monitoring in the Hebei Plain. Remote Sens. 2021, 13, 795. [Google Scholar] [CrossRef]
  35. Sabrian, P.G.; Saepuloh, A.; Kashiwaya, K.; Koike, K. Combined SBAS-InSAR and geostatistics to detect topographic change and fluid paths in geothermal areas. J. Volcanol. Geotherm. Res. 2021, 416, 107272. [Google Scholar] [CrossRef]
  36. Chen, M.; Tomás, R.; Li, Z.; Motagh, M.; Li, T.; Hu, L.; Gong, H.; Li, X.; Yu, J.; Gong, X. Imaging Land Subsidence Induced by Groundwater Extraction in Beijing (China) Using Satellite Radar Interferometry. Remote Sens. 2016, 8, 486. [Google Scholar] [CrossRef]
  37. Chitsazan, M.; Rahmani, G.; Ghafoury, H. Land subsidence susceptibility mapping using PWRSTFAL framework and analytic hierarchy process: Fuzzy method (case study: Damaneh-Daran Plain in the west of Isfahan Province, Iran). Environ. Monit. Assess. 2022, 194, 192. [Google Scholar] [CrossRef]
  38. Arabameri, A.; Saha, S.; Roy, J.; Tiefenbacher, J.P.; Cerda, A.; Biggs, T.; Pradhan, B.; Ngo, P.T.T.; Collins, A.L. A novel ensemble computational intelligence approach for the spatial prediction of land subsidence susceptibility. Sci. Total Environ. 2020, 726, 138595. [Google Scholar] [CrossRef]
  39. Chen, B.; Gong, H.; Chen, Y.; Lei, K.; Zhou, C.; Si, Y.; Li, X.; Pan, Y.; Gao, M. Investigating land subsidence and its causes along Beijing high-speed railway using multi-platform InSAR and a maximum entropy model. Int. J. Appl. Earth Obs. Geoinf. 2021, 96, 102284. [Google Scholar] [CrossRef]
  40. Hakim, W.; Achmad, A.; Lee, C.-W. Land Subsidence Susceptibility Mapping in Jakarta Using Functional and Meta-Ensemble Machine Learning Algorithm Based on Time-Series InSAR Data. Remote Sens. 2020, 12, 3627. [Google Scholar] [CrossRef]
  41. Li, H.; Zhu, L.; Dai, Z.; Gong, H.; Guo, T.; Guo, G.; Wang, J.; Teatini, P. Spatiotemporal modeling of land subsidence using a geographically weighted deep learning method based on PS-InSAR. Sci. Total Environ. 2021, 799, 149244. [Google Scholar] [CrossRef]
  42. Mohammady, M.; Pourghasemi, H.R.; Amiri, M.; Tiefenbacher, J.P. Spatial modeling of susceptibility to subsidence using machine learning techniques. Stoch. Environ. Res. Risk Assess. 2021, 35, 1689–1700. [Google Scholar] [CrossRef]
  43. Liu, Q.; Zhang, Y.; Deng, M.; Wu, H.; Kang, Y.; Wei, J. Time series prediction method of large-scale surface subsidence based on deep learning. Acta Geod. Cartogr. Sin. 2021, 50, 396–404. [Google Scholar] [CrossRef]
  44. Wang, Y.; Yang, G. Prediction of Composite Foundation Settlement Based on Multi-Variable Gray Model. Appl. Mech. Mater. 2014, 580–583, 669–673. [Google Scholar] [CrossRef]
  45. Zhou, D.; Zuo, X.; Zhao, Z. Constructing a Large-Scale Urban Land Subsidence Prediction Method Based on Neural Network Algorithm from the Perspective of Multiple Factors. Remote Sens. 2022, 14, 1803. [Google Scholar] [CrossRef]
  46. Hill, P.; Biggs, J.; Ponce-López, V.; Bull, D. Time-Series Prediction Approaches to Forecasting Deformation in Sentinel-1 InSAR Data. J. Geophys. Research. Solid Earth: JGR 2021, 126, e2020JB020176. [Google Scholar] [CrossRef]
  47. Ding, Q.; Shao, Z.; Huang, X.; Altan, O.; Zhuang, Q.; Hu, B. Monitoring, analyzing and predicting urban surface subsidence: A case study of Wuhan City, China. Int. J. Appl. Earth Obs. Geoinf. 2021, 102, 102422. [Google Scholar] [CrossRef]
  48. Bai, S.; Zico Kolter, J.; Koltun, V.J.a.e.-p. An Empirical Evaluation of Generic Convolutional and Recurrent Networks for Sequence Modeling. arXiv 2018, arXiv:1803.01271. Available online: https://ui.adsabs.harvard.edu/abs/2018arXiv180301271B (accessed on 20 March 2022).
  49. Jiang, Y.; Yang, K.; Shao, C.; Zhou, X.; Zhao, L.; Chen, Y.; Wu, H. A downscaling approach for constructing high-resolution precipitation dataset over the Tibetan Plateau from ERA5 reanalysis. Atmos. Res. 2021, 256, 105574. [Google Scholar] [CrossRef]
  50. Rosen, P.A.; Gurrola, E.M.; Agram, P.; Cohen, J.; Lavalle, M.; Riel, B.V.; Fattahi, H.; Aivazis, M.A.G.; Simons, M.; Buckley, S.M. The InSAR Scientific Computing Environment 3.0: A Flexible Framework for NISAR Operational and User-Led Science Processing. In Proceedings of the IGARSS 2018—2018 IEEE International Geoscience and Remote Sensing Symposium, Valencia, Spain, 22–27 July 2018; pp. 4897–4900. [Google Scholar] [CrossRef]
  51. Zhang, Y.; Heresh, F.; Falk, A. Small baseline InSAR time series analysis: Unwrapping error correction and noise reduction. Comput. Geosci. 2019, 133, 104331. [Google Scholar] [CrossRef]
  52. Li, H.; Zhu, L.; Gong, H.; Sun, H.; Yu, J. Land subsidence modelling using a long short-term memory algorithm based on time-series datasets. Proc. Int. Assoc. Hydrol. Sci. 2020, 382, 505–510. [Google Scholar] [CrossRef]
  53. Yu, Y.; Si, X.; Hu, C.; Zhang, J. A Review of Recurrent Neural Networks: LSTM Cells and Network Architectures. Neural Comput. 2019, 31, 1235–1270. [Google Scholar] [CrossRef] [PubMed]
  54. Naghibi, S.A.; Khodaei, B.; Hashemi, H. An integrated InSAR-machine learning approach for ground deformation rate modeling in arid areas. J. Hydrol. 2022, 608, 127627. [Google Scholar] [CrossRef]
  55. Zhang, H.; Song, Y.; Xu, S.; He, Y.; Li, Z.; Yu, X.; Liang, Y.; Wu, W.; Wang, Y. Combining a class-weighted algorithm and machine learning models in landslide susceptibility mapping: A case study of Wanzhou section of the Three Gorges Reservoir, China. Comput. Geosci. 2022, 158, 104966. [Google Scholar] [CrossRef]
  56. Chen, Y.; He, Y.; Zhang, L.; Chen, Y.; Pu, H.; Chen, B.; Gao, L. Prediction of InSAR deformation time-series using a long short-term memory neural network. Int. J. Remote Sens. 2021, 42, 6919–6942. [Google Scholar] [CrossRef]
  57. He, Y.; Yao, H.; Yang, W.; Yao, S.; Zhang, L.; Chen, Y.; Liu, T. Time-Series Analysis and Prediction of Surface Deformation in the Jinchuan Mining Area, Gansu Province, by Using InSAR and CNN–PhLSTM Network. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 2022, 15, 6732–6751. [Google Scholar] [CrossRef]
  58. Liu, Q.; Zhang, Y.; Wei, J.; Wu, H. and Deng, M. HLSTM: Heterogeneous Long Short-Term Memory Network for Large-Scale InSAR Ground Subsidence Prediction. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 2021, 14, 8679–8688. [Google Scholar] [CrossRef]
  59. Li, R.; Sun, J. Monitoring and prediction of tailings pond settlement based on integration of SBAS-InSAR and GS-LSTM. Met. Mine 2023, 1, 102–109. [Google Scholar] [CrossRef]
Figure 1. The geographic location of the study area, the coverage of the Sentinel-1 image tiles, and the distribution of leveling benchmarks are superimposed on the ESRI satellite topographic map.
Figure 1. The geographic location of the study area, the coverage of the Sentinel-1 image tiles, and the distribution of leveling benchmarks are superimposed on the ESRI satellite topographic map.
Remotesensing 15 02843 g001
Figure 2. List of specific dates of SAR image data.
Figure 2. List of specific dates of SAR image data.
Remotesensing 15 02843 g002
Figure 3. The method flowchart.
Figure 3. The method flowchart.
Remotesensing 15 02843 g003
Figure 4. Partial SAR fusion experiment results for track 40 and track 142. (a) indicates the SAR experimental results of track 40 in the overlapping region, (b) indicates the SAR experimental results of track 142 in the overlapping region, and (c) represents the final fusion results after adding the offset.
Figure 4. Partial SAR fusion experiment results for track 40 and track 142. (a) indicates the SAR experimental results of track 40 in the overlapping region, (b) indicates the SAR experimental results of track 142 in the overlapping region, and (c) represents the final fusion results after adding the offset.
Remotesensing 15 02843 g004
Figure 5. Schematic diagram of the network unit structure of the LSTM model.
Figure 5. Schematic diagram of the network unit structure of the LSTM model.
Remotesensing 15 02843 g005
Figure 6. Architectural diagrams of the LSTM model (a), TCN model (b), and LSTM-TCN model (c).
Figure 6. Architectural diagrams of the LSTM model (a), TCN model (b), and LSTM-TCN model (c).
Remotesensing 15 02843 g006
Figure 7. The semiannual average spring and summer deformation rates of the land in the study area. The values are obtained by dividing the cumulative deformation values in spring and summer of 2020 and 2021 by 2. Positive values represent the rate of land uplift, and negative values represent the rate of land subsidence.
Figure 7. The semiannual average spring and summer deformation rates of the land in the study area. The values are obtained by dividing the cumulative deformation values in spring and summer of 2020 and 2021 by 2. Positive values represent the rate of land uplift, and negative values represent the rate of land subsidence.
Remotesensing 15 02843 g007
Figure 8. The semiannual average deformation rate of the land in the study area in autumn and winter. The values are obtained by dividing the cumulative deformation values in autumn and winter of 2020 and 2021 by 2. Positive values represent the rate of land uplift, and negative values represent the rate of land subsidence.
Figure 8. The semiannual average deformation rate of the land in the study area in autumn and winter. The values are obtained by dividing the cumulative deformation values in autumn and winter of 2020 and 2021 by 2. Positive values represent the rate of land uplift, and negative values represent the rate of land subsidence.
Remotesensing 15 02843 g008
Figure 9. Correlation diagram between subsidence point data and leveling point measured data in the study area.
Figure 9. Correlation diagram between subsidence point data and leveling point measured data in the study area.
Remotesensing 15 02843 g009
Figure 10. Average annual land deformation rates from November 2019 to February 2022. A positive number represents land uplift, and a negative number represents land subsidence.
Figure 10. Average annual land deformation rates from November 2019 to February 2022. A positive number represents land uplift, and a negative number represents land subsidence.
Remotesensing 15 02843 g010
Figure 11. Cumulative land deformation in the study area from November 2019 to February 2022. A positive number represents land uplift, and a negative number represents land subsidence.
Figure 11. Cumulative land deformation in the study area from November 2019 to February 2022. A positive number represents land uplift, and a negative number represents land subsidence.
Remotesensing 15 02843 g011
Figure 12. Geographical distribution map of all experimental subsidence points.
Figure 12. Geographical distribution map of all experimental subsidence points.
Remotesensing 15 02843 g012
Figure 13. Loss curves of 12 parameter combinations for three models. (a) Training Loss and (b) Validation Loss.
Figure 13. Loss curves of 12 parameter combinations for three models. (a) Training Loss and (b) Validation Loss.
Remotesensing 15 02843 g013
Figure 14. The correlation coefficient graph between the predictions of the three models and the results of SBAS-InSAR. (ac) Correlation diagrams of the LSTM, TCN, and LSTM-TCN models, respectively.
Figure 14. The correlation coefficient graph between the predictions of the three models and the results of SBAS-InSAR. (ac) Correlation diagrams of the LSTM, TCN, and LSTM-TCN models, respectively.
Remotesensing 15 02843 g014
Figure 15. Location map of the four validation points.
Figure 15. Location map of the four validation points.
Remotesensing 15 02843 g015
Figure 16. Comparison of the prediction results of points (ad) with the monitoring results of SBAS-InSAR.
Figure 16. Comparison of the prediction results of points (ad) with the monitoring results of SBAS-InSAR.
Remotesensing 15 02843 g016
Figure 17. Experimental results of the three models with different numbers of samples. (ad) represent the RMSE, R 2 , MAE, and MAPE results of the three models, respectively.
Figure 17. Experimental results of the three models with different numbers of samples. (ad) represent the RMSE, R 2 , MAE, and MAPE results of the three models, respectively.
Remotesensing 15 02843 g017
Figure 18. Influence of different monthly sequence input matrices on experimental results. (ad) represent the RMSE, R 2 , MAE, and MAPE results of the three models, respectively.
Figure 18. Influence of different monthly sequence input matrices on experimental results. (ad) represent the RMSE, R 2 , MAE, and MAPE results of the three models, respectively.
Remotesensing 15 02843 g018
Table 1. Information about Sentinel-1A datasets.
Table 1. Information about Sentinel-1A datasets.
Orbital NumberQuantityPolarizations
Track 1184VV + VH
Track 113112VV + VH
Track 40112VV + VH
Track 14256VV + VH
Table 2. Volumetric soil aquifer information.
Table 2. Volumetric soil aquifer information.
NameDepthRange of ElevationSoil Type
Volume soil layer 10~7 cm23.2~2413.8 mSilty clay, soft soil
Volume soil layer 27~28 cm23.1~2414.5 mSilty clay, soft soil
Volume soil layer 328~100 cm22.9~2412.8 mSilty clay, soft soil
Volume soil layer 4100~289 cm22.3~2410.8 mSilty clay, soft soil
Table 3. The average cumulative subsidence and environmental factor information of all subsidence points in the study area in spring and summer, autumn and winter.
Table 3. The average cumulative subsidence and environmental factor information of all subsidence points in the study area in spring and summer, autumn and winter.
IndicatorsThe Semiannual Average in Spring and SummerThe Semiannual Average in Autumn and Winter
Accumulated subsidence−17.61 mm/semiannual−16.42 mm/semiannual
Precipitation619.29 mm/semiannual229.61 mm/semiannual
Volumetric soil layer water content (first layer)1.33 m3m−3/semiannual1.62 m3m−3/semiannual
Volumetric soil layer water content (second layer)1.29 m3m−3/semiannual1.64 m3m−3/semiannual
Volumetric soil layer water content (third layer)1.24 m3m−3/semiannual1.78 m3m−3/semiannual
Volumetric soil layer water content (fourth layer)1.20 m3m−3/semiannual1.65 m3m−3/semiannual
Table 4. Table of parameter combinations for the three different models.
Table 4. Table of parameter combinations for the three different models.
MethodslrBatch SizeEpochsLSTM LayerAtrous Convolution Layers
LSTML10.00023210020
L20.00026410030
L30.000212810030
L40.00036410030
TCNT10.00023210003
T20.00026410004
T30.000212810005
T40.00036410004
LSTM-TCNLT10.00023210013
LT20.00026410014
LT30.000212810015
LT40.00036410013
Table 5. The results of three models under other precision indexes. The 0~100 mm mean absolute error indicates MAE values of InSAR monitoring results and model prediction results in the range of 0~100 mm, and the 100~200 mm mean absolute error indicates MAE values of InSAR monitoring results and model prediction results in the range of 100~200 mm.
Table 5. The results of three models under other precision indexes. The 0~100 mm mean absolute error indicates MAE values of InSAR monitoring results and model prediction results in the range of 0~100 mm, and the 100~200 mm mean absolute error indicates MAE values of InSAR monitoring results and model prediction results in the range of 100~200 mm.
MethodsLSTMTCNLSTM-TCN
MAE (mm)6.163.511.58
MAPE (%)8.794.691.75
The prediction error is less than 3 mm (%)34.7155.4187.70
0~100 mm mean absolute error (mm)5.813.121.44
100~200 mm mean absolute error (mm)6.235.481.97
Table 6. Results of model experiments, depending on whether environmental factor variables are included.
Table 6. Results of model experiments, depending on whether environmental factor variables are included.
MethodsRMSE
(mm)
R 2 MAE
(mm)
MAPE
(%)
Prediction Error Is Less than 3 mm (%)
LSTM
(Univariate)
9.880.847.039.7229.35
LSTM8.400.896.168.7934.71
TCN
(Univariate)
6.270.944.546.4246.06
TCN4.770.963.514.6955.41
LSTM-TCN
(Univariate)
5.430.954.175.1950.52
LSTM-TCN2.170.991.581.7587.70
Table 7. Groundwater level data in Zhengzhou.
Table 7. Groundwater level data in Zhengzhou.
Soil TypeUnderground DepthSoil Properties
Powdery clay8~13.2 mYellowish-brown, high dry strength, high toughness.
Powdered clay5.4~8 mYellowish-brown, low dry strength, low toughness.
Powdered sand13.2~14.5 mBrownish-yellow, mainly chalky sand interspersed with a lot of chalky soil.
Fine Sand14.5~26.7 mBrownish-yellow, with a few mica fragments, and average grain gradation.
Table 8. Comparison of predicted land subsidence results for different input impact factors.
Table 8. Comparison of predicted land subsidence results for different input impact factors.
Study AreaInput FactorNumber of Subsidence PointsRMSE (mm)MAE (mm)
Zhengzhou CityVolumetric soil water content (1,2,3,4), precipitation.2011.921.21
Zhengzhou CityMonthly stabilization of groundwater level, precipitation.2011.481.01
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Guo, H.; Yuan, Y.; Wang, J.; Cui, J.; Zhang, D.; Zhang, R.; Cao, Q.; Li, J.; Dai, W.; Bao, H.; et al. Large-Scale Land Subsidence Monitoring and Prediction Based on SBAS-InSAR Technology with Time-Series Sentinel-1A Satellite Data. Remote Sens. 2023, 15, 2843. https://doi.org/10.3390/rs15112843

AMA Style

Guo H, Yuan Y, Wang J, Cui J, Zhang D, Zhang R, Cao Q, Li J, Dai W, Bao H, et al. Large-Scale Land Subsidence Monitoring and Prediction Based on SBAS-InSAR Technology with Time-Series Sentinel-1A Satellite Data. Remote Sensing. 2023; 15(11):2843. https://doi.org/10.3390/rs15112843

Chicago/Turabian Style

Guo, Hengliang, Yonghao Yuan, Jinyang Wang, Jian Cui, Dujuan Zhang, Rongrong Zhang, Qiaozhuoran Cao, Jin Li, Wenhao Dai, Haoming Bao, and et al. 2023. "Large-Scale Land Subsidence Monitoring and Prediction Based on SBAS-InSAR Technology with Time-Series Sentinel-1A Satellite Data" Remote Sensing 15, no. 11: 2843. https://doi.org/10.3390/rs15112843

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

Article Metrics

Back to TopTop