Monitoring and Analysis of Land Subsidence in Cangzhou Based on Small Baseline Subsets Interferometric Point Target Analysis Technology

Abstract

Cangzhou is located in the northeast part of the North China Plain; here, groundwater is the main water source for production and living. Due to the serious regional land subsidence caused by long-term overexploitation of groundwater, the monitoring of land subsidence in this area is significant. In this paper, we used the Small Baseline Subsets Interferometric Point Target Analysis (SBAS-IPTA) technique to process the Envisat-ASAR, Radarsat-2, and Sentinel-1A data and obtained the land subsidence of Cangzhou from 2004 to 2020. Additionally, we obtained winter wheat distribution information in Cangzhou using the Pixel Information Expert Engine (PIE-Engine) remote sensing cloud platform. On this basis, we analyzed the relationship between ground water level, winter wheat planting area, and the response of land subsidence according to the land use type and groundwater level monitoring data near the winter wheat growing area. The results show that during 2004–2020, the average annual subsidence rate of many places in Cangzhou was higher than 30 mm/year, and the maximum subsidence rate was 115 mm/year in 2012. From 2004 to 2020, the area of the subsidence funnel showed a trend of first increasing and then decreasing. In 2020, the subsidence funnel area reached 6.9 × 10³ km². The winter wheat planting area in the urban area showed a trend of first decreasing, then increasing and then decreasing, and it accounted for a large proportion in the funnel area. At the same time, we studied the relationship between the land subsidence rate and the water level at different burial depths and the response of winter wheat planting area. The results showed that the change of confined water level had a stronger response with the other two variables.

1. Introduction

Land subsidence is a slowly occurring environmental geological phenomenon caused by natural and human factors, which will seriously affect human production and lifestyles [1]. This kind of environmental phenomenon caused by multiple factors has the characteristics of a long formation time, wide influence, difficult prevention and treatment, and difficult repair [2]. It has become one of the geological engineering problems in major cities of the world because of the wide range of occurrence and the difficulty of in detecting it.

With the continuous increase in human water demand and consumption, the groundwater level in many countries and regions is constantly declining, which leads to the occurrence of different degrees of land subsidence disasters. For example, land subsidence caused by the intensive exploitation of groundwater for agricultural irrigation has been observed in various regions around the world. The San Joaquin Canyon [3] and Santa Clara Basin [4] in the United States, the Koyna River Basin [5] in India, the Rafsanjan Plain in Iran [6], the fore-alluvial plain at the eastern foot of the Taihang Mountains in China [7], and the southern margin of the Junggar Basin [8] are typical examples. In all of these areas, land subsidence caused by the extraction of groundwater for agricultural irrigation has been observed. Reeves et al. chose the irrigation circle of the central hub of the San Luis Valley as their research area [9]. They used European Remote Sensing Satellite (ERS) data to analyze the surface deformation field and inferred the changes in the water head of confined aquifers. They discovered that the deformation results exhibited the same seasonal periodicity as the monitored water head data. He et al. utilized Small Baseline Subset Interferometric Synthetic Aperture Radar (SBAS-InSAR) technology to monitor land subsidence in a loess irrigation area using Terra Synthetic Aperture Radar (TerraSAR) data [10]. They also conducted field investigations to determine the distribution of wells in the subsidence area; the main cause of the subsidence was the over-pumping of groundwater for agricultural irrigation. Wen et al. analyzed the InSAR measurement data of Heifangtai over a four-year period. They combined the InSAR data with the loess distribution and collapsibility characteristics in the area [11], and the results indicated that the collapsibility subsidence of the irrigation area of Heifangtai was gradually stable under long-term flood irrigation.

Cangzhou is one of the high-yield industrial and agricultural regions in the Beijing–Tianjin–Hebei Plain; it has relied on groundwater as its primary water source since the Beijing–Hangzhou Grand Canal (Cangzhou section) dried up in the 1960s [12]. Zhang et al. [13] proposed that the land subsidence in Cangzhou was mainly caused by the exploitation of deep groundwater, and the development of subsidence in the past ten years was basically consistent with the change in deep groundwater levels. Fang et al. [14] used the Fuzzy Analytical Hierarchy Process (FAHP) to evaluate the risk of land subsidence disaster in Cangzhou in 2016 based on the subsidence status, influencing factors, and disaster impact, and the results showed that Cangzhou was a high-risk area. Pan et al. [15] proposed that about 70% of the groundwater in the North China Plain is mainly used for irrigation to satisfy wheat cultivation and growth. Yan et al. [16] proposed that, in Cangzhou, there is a strong correlation coefficient between agricultural development and the depth of the groundwater. If agricultural development and utilization are reduced by 10%, the groundwater level will increase.

Therefore, we aimed to reveal the relationship between ground water level, land subsidence, and the response of winter wheat planting area in Cangzhou. In this study, the long-time series deformation characteristics of Cangzhou were obtained from Envisat-ASAR, Radarsat-2, and Sentinel-1A data using SBAS-IPTA technology. The content of this paper is as follows: Section 2 provides an overview of the geological background and research data of Cangzhou. Section 3 introduces the research methods used in this study. Section 4 includes the analysis of the temporal and spatial evolution characteristics of the land subsidence. Finally, in Section 5, we discuss the relationship between ground water level, land subsidence, and the response of winter wheat planting area in Cangzhou. The conclusion of this paper is finally presented in Section 6.

2. Study Area and Research Data

2.1. Study Area

Cangzhou is located in the central and eastern part of the North China Plain, adjacent to the Bohai Bay. Its longitude is 115°42′ E~117°50′ E, and its latitude is 37°29′ N~38°57′ N (Figure 1). It is an important channel connecting Beijing and Tianjin to the eastern coastal areas. Cangzhou is a prefecture-level administrative unit which comprises urban as well as agricultural areas. All the studied areas in this paper represent the agricultural areas of Cangzhou. The total area of Cangzhou is about 14.3 × 10³ km², and it is characterized by low and flat terrain, sloping from southwest to northeast, with a ground elevation of 2–15 m. Cangzhou is located in the mid-latitude region. The climate type is mainly warm temperate sub-humid continental monsoon, and it belongs to the semi-arid region. The temperature difference between the seasons is significant, and the four seasons are distinct. Sufficient light, along with spring drought, summer waterlogging, autumn coolness, and dry winter characteristics, make Cangzhou suitable for the development of agricultural production. In addition, the annual average temperature in Cangzhou ranges between 10 and 25 °C. The coldest month is January, with an average temperature of about 12.5 °C; the annual average precipitation is 581 mm, and the frost-free period lasts for 181 days [17,18]. The seasons vary in length, and the temperature varies greatly. The geomorphic genetic types of Cangzhou include alluvial, lacustrine, and marine deposits. According to the genetic types and geomorphic morphological characteristics, Cangzhou can be divided into three geomorphic types: an alluvial plain area, alluvial marine plain area, and marine plain area. There are many rivers in Cangzhou, but most of them are seasonal rivers and closed basins. There are 15 rivers in the territory, including Ziya New River, Zhangwei New River, South Canal, and Jiedijian River, with a total river length of 714 km. Due to its geographical location, the shallow groundwater in most areas of Cangzhou is saline water. As a result, the extraction of deep groundwater has become the main source for industrial and agricultural water [19]. Cangzhou is one of the major production areas for grain, cotton, and oil in Hebei Province. In 2022, the total area of land used for grain cultivation in Cangzhou reached 9.02 × 10³ km², which accounted for 63% of the city’s total area. Cangzhou is a traditional agricultural city. The main crops planted in Cangzhou follow a planting system of winter wheat and summer corn, with two crops a year [20]. Winter wheat is typically sown from late September to mid-October and emerges in late October. It enters the overwintering period from mid-December to late February of the following year. Winter wheat starts to return during the green period from late February to late March, and the jointing stage takes place in early to mid-April. The heading and flowering stage of winter wheat occurs in mid-May, while the milky maturity stage of winter wheat takes place in early and mid-June [21]. In this cycle, winter wheat utilizes more water during the sowing period (from late September to mid-October), returning green period (from late February to early March), jointing stage (from mid-April), heading stage (from late April to early May), and milky maturity stage (from mid-May) [22]. Moreover, the flat terrain, high soil fertility, abundant sunlight and heat resources, and consistent rainfall and temperatures during the same period create a regional environment that is highly conducive to the growth and development of winter wheat and other crops.

2.2. Research Data

2.2.1. Radar Remote Sensing Data

The radar remote sensing data include Envisat-ASAR data, Radarsat-2 data, and Sentinel-1A data, all of which are in C-band and VV polarization modes. The data coverage is shown in Figure 1. The Envisat-ASAR data include four kinds of orbit data from 10 December 2003, to 29 September 2010. The Radarsat-2 data mainly include 40 descending orbit datasets from 28 January 2012 to 21 October 2016. The Sentinel-1 data mainly include 68 ascending tracks from 14 January 2016, to 1 September 2020. However, the SAR data from 2011 are not included in this paper due to the absence of Radarsat-2 data in Hebei. The data list is shown in the

Table 1. Radar image information.

Satellite Identification	Data	Track	Frame	Beam Mode	Width Range (km)	Number of Images
Envisat-ASAR	10 December 2003– 29 September 2010	218	2835	Stripmap Mode	100 × 100	46
			2853			43
		447	2835			45
			2853			42
Radarsat-2	28 January 2012– 21 October 2016	21,528	4	Wide Mode	150 × 150	40
Sentinel-1A	14 January 2016– 1 September 2020	142	126	Interferometric Wideswath	250	68

.

2.2.2. Optical Remote Sensing

In order to determine the distribution of winter wheat planting in the study area from 2005 to 2020, this study utilized the Pixel Information Expert Engine (PIE-Engine) remote sensing cloud platform, in conjunction with phenological information of winter wheat in Cangzhou. Remote sensing images covering the entire winter wheat growing period in the study area (mid-October 2004 to mid-June 2020) were selected. The PIE-Engine can directly obtain the optical remote sensing datasets required for this study. These images included three remote sensing datasets: Landsat 5 TM, Landsat 7 TOA, and Sentinel-2 L2A [23,24]. The data are presented in

Table 2. Optical data source information.

Satellite Identification	Data	Spatial Resolution (m)	Return Period (Days)	Width Range (km)	Number of Band Classes
Landsat 5 TM	October 2004– June 2011	30	16	185	7
Landsat7 TOA	January 2012– June 2017	30	16	185	7
Sentinel-2 L2A	October 2017– June 2020	20	5	290	13

.

2.2.3. Other Data

The Digital Elevation Model (DEM) data were the Shuttle Radar Topography Mission (SRTM) data provided by the National Aeronautics and Space Administration (NASA) with a spatial resolution of 90 m. The data used to verify the inter-annual planting area of winter wheat monitored by the PIE-engine were sourced from the Cangzhou Statistical Yearbook from 2005 to 2022 (http://tj.cangzhou.gov.cn/, accessed on 17 October 2023) [25]. The data of the groundwater level used to overlay the subsidence rate and inter-annual planting area of winter wheat were extracted from the China Groundwater Yearbook from 2012 to 2020, published by the China Institute of Geological Environment Monitoring (https://geocloud.cgs.gov.cn/, accessed on 17 October 2023).

3. Research Method and Data Processing

3.1. Research Methods

3.1.1. Small Baseline Subsets Interferometric Point Target Analysis (SBAS-IPTA) Technology

IPTA technology interferes with Permanent Scatterers (PS) point sets that have stable spectral characteristics or high backscattering characteristics on Single Look Complex (SLC) data. The interference points, including the terrain phase, deformation phase, atmospheric phase, and phase noise are obtained. Compared to conventional PS-InSAR technology, this approach stores intermediate data in a vector format, which can significantly enhance the speed of data calculation [26,27]. However, its ability is often limited when processing wide-format data, such as Radarsat-2 data, for whole-format processing. However, SBAS-InSAR technology has low requirements for the monitoring range and is suitable for complex scenarios [12]. Therefore, this study utilized SBAS-InSAR technology and IPTA technology to accomplish the processing of SAR image data. The main technical process is shown in Figure 2.

(1)

Selection and registration of the reference image

The reference image was selected using the global correlation coefficient method, which is based on minimizing the time baseline, space baseline, and Doppler centroid frequency shift. The reference images of orbit data 218 and 447 from the Envisat-ASAR data in this paper are from 2 April 2008 and 14 March 2008, respectively. The reference images of the Radarsat-2 data and Sentinel-1A data are from 23 April 2014 and 25 June 2017, respectively. After determining the reference image, the other remaining images were registered to it using the same image coordinates.

(2)

Differential interferogram generation

We used MATLAB software (version: R2019)to select the time baseline and then GAMMA software (version: 20131203) to process the SAR data [28]. The N + 1 registered images were geocoded. The time baseline and spatial threshold were set. In this study, the time baseline and spatial baseline were set to 300 days and 300 m, respectively, for the three datasets. The final interference pairs for Envisat-ASAR were 148, 138, 143, and 132. The final interference pair for Radarsat-2 was 284, and the final interference pair for Sentinel-1A was 174. At the same time, the SRTM DEM with a resolution of 90 m was selected as the externally imported DEM data, and the terrain phase was removed using these data. The differential interference phase ${\phi }_{dint}$ in each interferogram consisted of several components, as shown in the following:

$${\phi }_{dint}={\phi }_{flat}+{\phi }_{topo}+{\phi }_{def}+{\phi }_{atm}+{\phi }_{noise}$$

(1)

where ${\phi }_{flat}$ represents the phase change caused by the change in curvature of the earth; ${\phi }_{topo}$ represents the phase change caused by topographic relief; ${\phi }_{def}$ is the deformation phase change caused by the radar target, which includes both linear deformation and nonlinear deformation phases; ${\phi }_{atm}$ is the phase change caused by the change in the electromagnetic wave path due to the non-uniform characteristics of atmospheric composition; and ${\phi }_{noise}$ is the phase change caused by residual noise.

(3)

Coherent point target extraction

The recognition and selection of highly stable coherent point targets are the basis and key of InSAR interferometry. The point target extraction method used in this study was the coherence coefficient threshold method. This method determines the appropriate coherence coefficient threshold by analyzing the image interference coherence map. It identifies the coherence point target as any highly coherent target point that exceeds the specified threshold [29]. The coherence coefficient threshold was set to 0.78 in this study.

(4)

Coherent point target deformation analysis

An iterative regression analysis was performed based on the selected reference points. This process involves successive iterations of linear deformation rate correction, elevation correction, a residual phase, and an unwinding interference phase. The linear deformation rate correction can be used to update the deformation model. Elevation correction results in an improved topographic phase model. At the same time, the residual phase is decomposed through spatiotemporal filtering, separating the atmospheric phase from the nonlinear deformation phase. The minimum cost flow phase unwrapping method is then used for spatial unwrapping. Finally, the linear deformation phase and nonlinear deformation phase are superimposed, and the deformation information on the time series is obtained through Singular Value Decomposition (SVD).

(5)

Merging the InSAR Monitoring Results

In this study, we only considered the difference in deformation in point targets that was caused by a different selection of reference points [30]. We first acquired the deformation for each SAR image using the IPTA technique. Secondly, we extracted the same interferometric point targets in overlapping regions. Thirdly, the deformation information fusion model based on the least square was constructed, and the deformation monitoring results of adjacent strip were fused to obtain the land subsidence information [28]. Meanwhile, the InSAR results derived by the Envisat-ASAR, Radarsat-2, and Sentinel-1A images were transformed from line of sight (LOS) to the vertical direction by (2) to ensure that the InSAR results have a common direction:

$$d_u=\frac{d_{LOS}}{\mathit{cos}\theta }$$

(2)

where $\theta$ is the central incidence angle, and $d_{LOS}$ is the deformation in the LOS direction. Finally, the long-term subsidence information of Cangzhou from 2004 to 2020 was derived.

3.1.2. Cangzhou Winter Wheat Planting Information Acquisition

In this study, three datasets were selected: Landsat 5 TM, Landsat 7 TOA, and Sentinel-2 L2A. These datasets cover the period from 2005 to 2020. Using the PIE-Engine remote sensing cloud platform, we processed remote sensing images of Cangzhou during the crucial period of winter wheat identification. The processing involved tasks such as image cropping, radiation calibration, cloud removal, Normalized Difference Vegetation Index (NDVI) calculation, image stitching, and image download. By utilizing NDVI remodel amplification and NDVI increase (decrease) slope thresholds, we established a universal remote sensing rapid mapping model for the winter wheat region. This model allowed us to map the inter-annual planting changes for winter wheat in Cangzhou from 2005 to 2020. The specific technical procedure is shown in Figure 3.

(1)

NDVI remodel amplification and NDVI increase (decrease) slope thresholds

The vegetation index can qualitatively and quantitatively evaluate vegetation biomass, growth activity, and vegetation coverage. It is an important indicator for characterizing vegetation cover through remote sensing. The common vegetation indices include the NDVI, Enhanced Vegetation Index (EVI), Ratio Vegetation Index (RVI), etc. In the research of index-based mapping, the NDVI reconstruction increment algorithm has been widely used to determine the index threshold [31].

This method determines a set of candidate thresholds by specifying the initial value, the end value, and the step size. Then, each candidate threshold is used to differentiate the target crop from the non-target crop, and the accuracy of the results is assessed. The threshold that corresponds to the highest accuracy is the optimal threshold. Based on the analysis of the phenological characteristics of winter wheat during its growth period, some scholars have established rules by observing the changes in NDVI thresholds of winter wheat during various key periods [32,33]. As a result, they have obtained extraction models for winter wheat with an accuracy of over 90%. This method fully utilizes the unique spectral features reflected in remote sensing images at various crucial stages during the winter wheat growth period. It enhances the distinction between winter wheat and other ground objects, thereby improving the accuracy of winter wheat extraction. The formula is shown in (3):

$$\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}=\frac{{\rho }_{Nir}-{\rho }_{Red}}{{\rho }_{Nir}+{\rho }_{Red}}$$

(3)

where ${\rho }_{Nir}$ and ${\rho }_{Red}$ represent the reflectance of the near-infrared band and the red band, respectively.

In this study, the NDVI index was calculated using the third and fourth bands of Landsat 5 TM, based on the band requirements of different datasets. Similarly, for Landsat 7 TOA, the third and fourth bands were selected to calculate the NDVI index. For Sentinel-2 L2A, Band 4 and Band 8 were selected for calculating the NDVI index. With pixels as the basic unit, the minimum NDVI value from 15 September to 15 November was determined for each pixel position. This resulted in an image representing the minimum NDVI value, which is referred to as ${\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_{\mathrm{m}\mathrm{i}\mathrm{n}}$. At each pixel position, the maximum NDVI value was filtered from 1 December to 31 March of the following year. This process results in an image displaying the maximum NDVI value, which is denoted as the maximum NDVI value (denoted as ${\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_{\mathrm{m}\mathrm{a}\mathrm{x}}$).

The period from 15 September to 15 November of the same year was the sowing and seedling period for winter wheat. During this period, the NDVI value of winter wheat was the lowest compared to other green vegetation crops. Therefore, the NDVI value during this period could effectively differentiate winter wheat from other crops. From 1 December to 31 March of the following year is the winter wheat overwintering and returning green stage. During this period, the NDVI value of winter wheat increases daily, allowing for effective differentiation between winter wheat and built-up lands or bare soil.

Since the NDVI vegetation index extraction alone may mistakenly include crops that were planted before winter wheat, the Normalized Burnt Ratio (NBR) vegetation index was used to filter them out. At the same time, the increase in NDVI (recorded as ${\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_{\mathrm{i}\mathrm{n}\mathrm{c}\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{s}\mathrm{e}}$) was used to minimize the spectral impact of neighboring non-winter wheat crops on the identification of winter wheat. The calculation formula for this method is shown in (4):

$${\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_{\mathrm{i}\mathrm{n}\mathrm{c}\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{s}\mathrm{e}}=\frac{{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_{\mathrm{m}\mathrm{a}\mathrm{x}}-{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}{\left|{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\right|}.$$

(4)

The results showed that ${\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_{\mathrm{i}\mathrm{n}\mathrm{c}\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{s}\mathrm{e}}$ and ${\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_{\mathrm{m}\mathrm{a}\mathrm{x}}$ were above 1.0 and 0.3, respectively, from sowing to 31 March of the following year. Therefore, this method was developed to extract winter wheat in Cangzhou.

(2)

Accuracy verification

The recognition results were evaluated, and the accuracy of the recognition was determined. In this study, the precision of the confusion matrix and the precision of the quantity were used to verify the accuracy of the extraction results of Cangzhou winter wheat. Among them, the accuracy of the confusion matrix was calculated by comparing the position of the verified sample pixel with the corresponding position in the recognized ground object image. The evaluation index for accuracy mainly included the Kappa coefficient and user accuracy.

The accuracy evaluation confusion matrix was constructed by combining the actual distribution vector data and the recognition results of winter wheat. The statistical data from the Cangzhou Municipal Bureau of Statistics from 2005 to 2020 were used as the true value to complete the quantitative accuracy evaluation [34,35,36].

4. Results and Analysis

4.1. Temporal and Spatial Evolution of Land Subsidence in Cangzhou

In this study, SBAS technology and IPTA technology were used to process Envisat-ASAR data from 2003 to 2010, Radarsat-2 data from 2012 to 2016, and Sentinel-1 data from 2016 to 2020. The monitoring results are shown in Figure 4.

The monitoring results showed that various degrees of uneven land subsidence occurred in Cangzhou between 2004 and 2020. In the past 16 years, Suning and Xian counties in the western part of Cangzhou and Dongguang County in the south of Cangzhou have been experiencing serious subsidence. The subsidence rate of Suning County increased from 54 mm/year in 2004 to 88 mm/year in 2020. The subsidence rate of Xian County increased from 60 mm/year in 2004 to 75 mm/year in 2020. Similarly, that of Dongguang County increased from 52 mm/year to 64 mm/year. On the other hand, Qing County’s subsidence showed a trend of decreasing year by year over the 16 years, decreasing from 48 mm/year in 2004 to 23 mm/year in 2020.

4.1.1. Center Transfer Law of Land Subsidence Funnel in Cangzhou

To provide a more detailed description of the spatial and temporal evolution patterns of land subsidence in Cangzhou, this study focused on extracting the center (the maximum point of land subsidence rate) of the land subsidence funnel (defined as the rate of land subsidence greater than 30 mm/year) in Cangzhou. The study also analyzed the transition characteristics of the center of the land subsidence funnel. The research results are shown in Figure 5 and

Table 3. The maximum annual subsidence rate and migration distance. (SFC means subsidence funnel center.)

Year	The Maximum Subsidence Rate of the SFC (mm/Year)	The Displacement of the SFC from the Previous Year (km)	Year	The Maximum Subsidence Rate of the SFC (mm/Year)	The Displacement of the SFC from the Previous Year (km)
2004	82	0	2013	108	852.55
2005	57	0.14	2014	105	0
2006	63	0.35	2015	105	0
2007	63	838.47	2016	92	1249.79
2008	66	0.01	2017	94	132.54
2009	53	0.17	2018	100	0
2010	55	25.58	2019	96	0.18
2012	115	783.58	2020	93	0.18

.

From 2004 to 2020, the center of the land subsidence funnel in Cangzhou migrated from Qing County to Suning County, moving spatially from the north of Cangzhou to the west. Additionally, the subsidence rate increased from 82 mm/year to 93 mm/year. In 2009, the subsidence rate at the center of the subsidence funnel was the lowest, which was 53 mm/year. In 2012, the subsidence rate at the center of the subsidence funnel was the highest, reaching 115 mm/year. From 2015 to 2016, the center of the subsidence funnel shifted from Cang County to Suning County in the west, covering a migration distance of up to 1250 km. The migration distances of 2012–2013 and 2006–2007 were also large, reaching 853 km and 838 km. In 2014 and 2015, the location of the subsidence funnel remained unchanged and was consistently located in Cang County, which is in the center of Cangzhou.

4.1.2. Variation Characteristics of Land Subsidence in the Funnel Area of Cangzhou

The temporal variation of the land subsidence funnel area in Cangzhou was further analyzed (

Table 4. The subsidence funnel area of Cangzhou from 2004 to 2020.

Year	The Area of Land Subsidence Funnel (×103 km2)	Year	The Area of Land Subsidence Funnel (×103 km2)
2004	2.8	2013	4.6
2005	3.1	2014	4.2
2006	3.0	2015	4.2
2007	3.9	2016	6.6
2008	3.3	2017	7.4
2009	2.4	2018	7.5
2010	3.0	2019	7.2
2012	7.6	2020	6.9

), and the spatial distribution variation characteristics of the land subsidence funnel in Cangzhou were analyzed (Figure 6). As can be seen from Table 4, the area of the land subsidence funnel in Cangzhou exhibited a trend of initially increasing and then decreasing from 2004 to 2020. From 2004 to 2010, the area of the land subsidence funnel showed a slow increase and then a decreasing trend, increasing from 2.8 × 10³ km² to 3.9 × 10³ km² and then decreasing to 3.0 × 10³ km². In 2012, the area of the land subsidence funnel reached the maximum (7.6 × 10³ km²). In 2020, the area of the land subsidence funnel was 6.9 × 10³ km², which indicated a decreasing trend but it was still higher than that in 2004.

5. Discussion

5.1. Relationship between Land Subsidence and Planting Distribution of Winter Wheat

Based on the PIE-Engine remote sensing cloud platform and three datasets, such as Landsat 5TM, this study utilized NDVI remodel amplification and NDVI increase (decrease) slope thresholds to extract the spatial distribution of winter wheat in Cangzhou. The study obtained the planting distribution of winter wheat in Cangzhou from 2005 to 2020. (The distribution infographic and comparison of accuracy are shown in the Appendix A.)

In order to calculate the winter wheat planting situation in the land subsidence funnel area of Cangzhou, the planting distribution information of winter wheat obtained from the cloud platform was spatially stacked with the subsidence funnel area. The proportion of the winter wheat area in the subsidence funnel area was then calculated, as shown in Figure 7 and Figure 8. As depicted in Figure 7, the distribution of winter wheat planting in the subsidence funnel area underwent a period of initial decline followed by expansion (2005–2010), and then entered a phase of gradual stability (2012–2020). From 2005 to 2010, the area affected by subsidence initially decreased and then gradually increased. In 2009, the subsidence funnel area was the smallest, only 2.4 × 10³ km², and the proportion of winter wheat planting in the subsidence funnel area was also the lowest at 74%. The area of the subsidence funnel gradually increased year by year, along with an increase in the proportion of winter wheat. In 2012, the maximum value was 7.8 × 10³ km², and the proportion of winter wheat planting in the subsidence funnel area was 97%. From 2012 to 2020, as mentioned above, with the adjustment of the agricultural planting structure in Cangzhou and the implementation of the policy to close self-provided wells, the proportion of winter wheat planting in the subsidence funnel area became increasingly stable year after year. From 2012 to 2020, the average proportion of winter wheat planting in the subsidence funnel area was 91%.

From the spatial distribution, it was found that the subsidence funnel distribution in Cangzhou was not uniform from 2005 to 2010. It was mainly concentrated in the western and northern regions of the city. The distribution of winter wheat in the subsidence funnel region, represented by Xian County in the west and Qing County in the north, showed a trend of initially decreasing and then increasing. In 2009, the proportion of winter wheat planting in the subsidence funnel area of Qing County, located in the north of Cangzhou, was significantly smaller than in previous years. From 2012 to 2020, the subsidence funnel area of Cangzhou was mainly concentrated in the western and southern regions, with the two parts showing a relatively concentrated and contiguous pattern. As can be seen from Figure 7, in 2012, winter wheat planting accounted for a significant proportion in Suning County in the west and Dongguang County in the south of Cangzhou. After 2012, the proportion of winter wheat planting in the subsidence funnel areas, represented by Xian County in the west and Dongguang County in the south, decreased slightly, but the distribution remained relatively stable.

Therefore, the winter wheat planting area in Cangzhou experienced a trend of initially decreasing, followed by an increase, and then slowly decreasing from 2005 to 2020. Spatially, the distribution of winter wheat planting in the south, north, and west of Cangzhou initially decreased, then increased, and finally stabilized. Therefore, based on this distribution difference, the trend of land subsidence in Cangzhou from 2005 to 2020 is similar to the distribution of winter wheat planting. It experienced a trend of initially decreasing, then expanding, then slowing down and becoming stable. Moreover, based on calculations and the analyses (Figure 8), it was observed that a greater amount of winter wheat was cultivated in the subsidence funnel area. This increased agricultural water usage in these areas and the surrounding regions may potentially contribute to land subsidence.

5.2. Response Characteristics of Different Land Use Types and Land Subsidence

In order to further analyze the spatio-temporal response process of land subsidence, winter wheat planting distribution, and groundwater levels in Cangzhou, we identified the soil utilization types of monitoring well locations, as shown in Figure 9. In light of the limited monitoring well data, the wells in the urban areas surrounding the Cangzhou sedimentation funnel were selected as a reference. As can be seen from

Table 5. Comparison of winter wheat planting areas obtained using PIE-Engine and the Cangzhou Statistical Yearbook from 2005 to 2020.

Well Number	Land Use Types	Subsidence Range (mm/Year)
Well 1	Agriculture and town residential land	9–112
Well 2	Agriculture and town residential land	1–96
Well 3	Agriculture, industrial, villages and town residential land	4–79

, the land use types of the three water level monitoring wells were similar. They are all located in the farmland area, with a large number of densely populated rural residential lands surrounding them. The maximum subsidence of the three wells was 112, 96, and 79 mm/year, respectively.

Well 1 (Figure 9a,b) is located near the subsidence funnel area of Qing County, Cangzhou. The area surrounding the well is primarily occupied by rubber manufacturing, non-ferrous metal processing, and new energy vehicle industries. Among these, the non-ferrous metal processing industry stands out as a high water consumption sector. In addition, the area is known for its thriving agricultural industry, with numerous farms cultivating crops such as corn, wheat, and soybeans. Consequently, the region’s agricultural water consumption is also substantial. According to historical documents, the drinking water for industry, agriculture, and people and animals in this area relied heavily on the extraction of deep groundwater resources [37]. This led to the formation of a large groundwater funnel, which played a significant role in the subsidence funnel area of Qing County, Cangzhou, by providing auxiliary support and accelerating the subsidence process.

Well 2 (Figure 9c,d) is located near the Suning subsidence funnel area of Cangzhou. This well is located in an area abundant in oil resources, with the surrounding region primarily focused on petroleum refining, new energy vehicle, and service textile industries. Among these, petroleum refining stands out as a high water consumption industry. In addition, in terms of agriculture, the region primarily cultivates crops such as wheat, corn, and cotton, which require significant amounts of water. Historical documents show that water consumption at this well location is primarily divided into industrial and domestic uses, which rely on deep groundwater. On the other hand, agricultural water consumption primarily relies on shallow groundwater [38]. Therefore, as a result of long-term over-exploitation of groundwater and high water demand from various sectors, a significant groundwater table depression funnel has formed in this region [39]. This has had a certain impact on the subsidence funnel area of Suning County, Cangzhou.

Well 3 (Figure 9e,f) is located in Nanpi County, Cangzhou. Nanpi County, as a traditional agricultural county, has been designated as a grain and cotton production base by the state. The main crops grown in Nanpi County are winter wheat, corn, cotton, and soybeans [40]. Among them, winter wheat and corn consume more water during their growth period. According to statistics, 73% of the well water in Nanpi County is used for agricultural production [41], making it a significant water source for agriculture. At the same time, the main industry in Nanpi County is the hardware and electrical industry, which has high water consumption and serves as the pillar industry [42]. Therefore, due to the coastal water shortage and salinization in Nanpi County, which is a representative area, the continuous exploitation of groundwater has led to a decrease in the depth of the groundwater. As a result, a groundwater funnel has formed, which has had a significant impact on the subsidence funnel area of Nanpi County.

5.3. Response Process of Subsidence Rate to Changes in Groundwater Level and Interannual Planting Area of Winter Wheat

In order to reveal the relationship between the interannual planting area of winter wheat, groundwater levels, and land subsidence rates, we selected data for superposition analysis from May 2012 to September 2020. The InSAR observations at each well point are the mean values of PS points in the buffer zone with a radius of 20 km from the well point.

As can be seen from Figure 10a, the temporal variation in the groundwater depth of confined well 1 from 2012 to 2020 can be divided into two stages: the groundwater level rising stage (May 2012 to December 2017) and the groundwater level stabilizing stage (January 2018 to September 2020). In the first stage, the groundwater level exhibited a fluctuating upward trend, with the groundwater depth decreasing from 74 m in 2012 to 50 m in 2017. At the same time, the subsidence rate of the well 1 showed a slowing trend, decreasing from 18 mm/year in 2012 to 4 mm/year in 2017. There was a strong correlation between the groundwater depth and the land subsidence and the correlation coefficient was 0.87. In the second stage, the groundwater depth exhibited a stable trend following a sharp increase. From January 2018 to September 2020, the groundwater depth increased to 4 m. In this stage, the change in the land subsidence rate was not significant and was approximately 4 mm/year. The correlation coefficient between the groundwater level change and land subsidence was 0.92. The trend of water level change in confined well 2 (Figure 10b) is consistent with that of confined well 1. From 2012 to 2020, the water level of confined well 2 can be divided into two main stages of development: the rising stage of the groundwater level (May 2012 to December 2016) and the stabilizing stage of the groundwater level (January 2018 to September 2020). In the first stage, the groundwater level exhibited a general upward trend, with the depth of the groundwater decreasing from 66 m in 2012 to 64 m in 2016. The subsidence rate of the well decreased from 6 mm/year in 2012 to 1 mm/year in 2016, while the correlation coefficient between them was 0.86. In the second stage, the fluctuation in the groundwater depth was small and was around 67 m. At the same time, the land subsidence rate also exhibited a consistent fluctuation pattern, with an average subsidence rate of approximately 2 mm/year. The correlation coefficient between them was 0.93. Through the analysis of the winter wheat planting situation, it was observed that the area of winter wheat declined from 2012 to 2020. Specifically, during the period from 2012 to 2017 (which coincided with the rising stage of groundwater levels), the area of winter wheat decreased by 23 × 10³ hm² (1 hm² = 0.01 km²). According to the analysis of the seasonal water use characteristics of winter wheat, the groundwater level of the well showed a downward trend, and the subsidence rate showed an upward trend during two periods: from the winter wheat greening period (from late February to late March) to the jointing period (mid-April), and from the heading period (from late April to early May) to the milky maturity stage (mid-May).

As can be seen from Figure 10c, the groundwater depth of phreatic well 1 exhibited a seasonal fluctuation trend from 2012 to 2020. Specifically, the groundwater depth at the well began to increase sharply in December 2017 and then stabilized. From May 2012 to December 2017, the depth of the groundwater fluctuated around 31 m, with variations of less than 1 m. In this stage, the land subsidence rate showed a decreasing trend, dropping from 16 mm/year in 2012 to 4 mm/year in 2017. Additionally, the change in groundwater depth exhibited a weak correlation with land subsidence, as indicated by a correlation coefficient of 0.46. From January 2018 to September 2020, the groundwater depth again showed stable fluctuation characteristics, and the variation amplitude was less than 4 m. The land subsidence rate showed a stable trend, and the subsidence rate was basically stable at 4 mm/year. The groundwater depth of submersible Well 2 (Figure 10d) is mainly divided into two development stages from 2012 to 2020, which are the initial rising and then decreasing stage (May 2012 to December 2016) and the stable stage (January 2018 to September 2020). In the first stage, the groundwater depth showed a decreasing trend, and decreased from 4 m in 2012 to 3 m in 2016. Simultaneously, the subsidence rate of the well exhibited a declining trend, decreasing from 6 mm/year in 2012 to 1 mm/year in 2016. The correlation coefficient between the two variables was 0.62. Overall, although the groundwater depth increased sharply in 2017, it generally exhibited a stable trend of fluctuation.

6. Conclusions

In this paper, we studied the spatial pattern and temporal evolution characteristics of land subsidence in Cangzhou by combining SBAS-IPTA results obtained from multi-source SAR datasets acquired from 2004 to 2020. The spatial–temporal distribution characteristics of winter wheat in Cangzhou from 2005 to 2020 were obtained using a PIE-Engine cloud remote sensing classification technique. On this basis, we analyzed the response characteristics of land subsidence evolution and winter wheat planting, and discussed the relationship between land subsidence rate and changes in the groundwater level at different levels as well as the response of winter wheat planting. The research conclusions are as follows:

• From 2004 to 2020, different degrees of land subsidence occurred in many places in Cangzhou, and the maximum annual subsidence rate increased from 60 mm/year to 82 mm/year. The center of the subsidence funnel moved from Qing County in the north of Cangzhou to Dongguang County in the south and then to Suning County in the west, during which, the maximum subsidence rate first decreased, then increased and then decreased, and the maximum migration distance reached 98 km.
• From 2004 to 2020, the area of land subsidence funnel in Cangzhou showed a trend of first increasing and then decreasing, among which, the area of the subsidence funnel reached the maximum in 2012, which was 7.6 × 10³ km². The planting distribution of winter wheat in the subsidence funnel area experienced a trend of first decreasing, expanding, then flattening and stabilizing year by year. In 2012, the proportion of winter wheat in the funnel area was the highest at 97%.
• The change in the groundwater level of the confined wells had a good correlation with land subsidence, and the correlation coefficient was always above 0.8. During the winter wheat irrigation period, the groundwater level of the confined water wells mostly showed a downward trend, while the subsidence rate showed an upward trend. The correlation between the changes in the groundwater level and land subsidence in phreatic wells was weak.