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

: 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 inﬂuence 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 inﬂuence 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 signiﬁcantly 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.


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 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.

Study Area
Henan Province is located in the southern part of the North China Plain. Henan Province has a total area of 167,000 km 2 , 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.
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.

Study Area
Henan Province is located in the southern part of the North China Plain. Henan Province has a total area of 167,000 km 2 , 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.

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 (singlelook). 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.

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 Shu le 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.

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

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.

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-singlelevels-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.

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.

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): 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 ∆ϕ top,j , ∆ϕ atm,j , and ∆ϕ noise,j represent the residual terrain phase difference, atmospheric phase difference, and noise phase difference, respectively.

SBAS-InSAR
The SBAS-InSAR method is a distributed sca erer (DS) method, and its main steps are described below. Assume that the +1 SAR images of the study area were obtained at times , , ⋯, and . 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 se ing the appropriate spatiotemporal baseline, interferograms were generated. For the -th interferogram generated by the combination of and ( < ), after removing the flat earth effect and terrain phase, the interference phase of any point in interferogram can be represented as follows by Formula (1): where represents the radar central wavelength; and represent the cumulative After removing the phase other than the deformation, the interferometric phase was simplified and represented by Equations (2) to (3): 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 Remote Sens. 2023, 15, 2843 8 of 28 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.
Remote Sens. 2023, 15, x FOR PEER REVIEW 8 of 29 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

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 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].
where represents the forget gate, and tanh represent the activation function, , , , and represent the input component of the weight matrix, ℎ and represent the hidden state value at the previous − 1 moment, represents the input value, represents the input unit and represents filtered values. , , , and represent the offset at time , represents the output gate unit, and ℎ represents the hidden state value at the previous time .
The operation between the dilated convolutional layer of the TCN model and the input sequence is given as follows: The main operations of this method are shown in Equations (4)- (8).
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: where the input sequence is X ∈ R n , the filter is f : {0, . . . , k − 1} → R , 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 .
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: 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 was used as the two-dimensional input matrix, and the land-cumulative subsidence data at time + 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 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.

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. Figures 7 and 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.

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 Figures 10 and 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.  The cumulative land deformation data and annual average land deformation rate data are shown in Figures 10 and 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.  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.

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/km 2 . 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. Figures 15 and 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 timeseries 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.

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.

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, 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, , 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.

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.

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.

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 recur-rent 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. In the prediction of time-series land subsidence, researchers used traditional mach learning algorithms to model and predict land subsidence from a data-driven perspecti For example, the boosted regression trees and extreme gradient boosting algorith 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.

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.

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.

Data Availability Statement:
The data presented in this study are available upon request from the corresponding author.