Land Subsidence Monitoring and Dynamic Prediction of Reclaimed Islands with Multi-Temporal InSAR Techniques in Xiamen and Zhangzhou Cities, China

: Artiﬁcial islands and land reclamation are one of the most important ways to expand urban space in coastal cities. Long-term consolidation of reclaimed material and compaction of marine sediments can cause ground subsidence, which may threaten the buildings and infrastructure on the reclaimed lands. Therefore, it is crucial to monitor the land subsidence and predict the future deformation trend to mitigate the damage and take measures for the land reclamation and any infrastructure. In this paper, a total of 125 SAR images acquired by the C-band Sentinel-1A satellite between June 2017 and September 2021 are collected. The small baseline subsets (SBAS) SAR interferometry (InSAR) method is ﬁrst conducted to detect the land deformation in Xiamen and Zhangzhou cities of Fujian Province, China, and the distributed scatterers (DS)-InSAR method is used to recover the complete deformation history of some typical areas including Xiamen Airport in Dadeng Island and Shuangyu Island. Then, the sequential estimation and the geotechnical model are jointly applied to demonstrate the current and future evolution of land subsidence of the constructed roads on Shuangyu Island. The results show that the maximum cumulative deformation reaches 425 mm of Xiamen Xiang’an Airport and 626 mm of Shuangyu Island, and the maximum deformation is predicted to be as large as 1.1 m by 2026 of Shuangyu Island. This research will provide important guidelines for the design and construction of Xiamen Xiang’an Airport and Shuangyu Island to prevent and control land subsidence.


Introduction
With the expansion of coastal reclamation, economic growth is inevitably accompanied by negative environmental and ecological issues such as biodiversity loss, water pollution, and wetland degradation [1,2]. At present, coastal reclamation projects have been carried out in many countries, such as in the USA [3], China [4,5], and Singapore [6]. Furthermore, many airports are built on the reclaimed or partially reclaimed land for public foundations, such as Changi Airport in Singapore [6] and Hong Kong International Airport in China [7]. As reclaimed land is usually on some unconsolidated marine sediment, the compaction of the soil can cause significant land subsidence over the time [8]. After land reclamation is complete, the serious land subsidence will lead to irreversible damage to buildings, roads, airport runways, and underground pipelines, which causes significant economic losses and high maintenance costs [5].
Land reclamation-induced deformation can be measured by traditional high precision technologies such as GNSS and leveling. However, these techniques are time consuming, costly, and operated on sparse benchmarks [9]. Over the last 30 years, repeat-pass spaceborne synthetic aperture radar interferometry (InSAR) has become one of the most powerful geodetic technology for high-resolution and large-coverage deformation measurements [10]. To overcome the limitations of conventional differential InSAR technologies, including DEM errors, satellite orbital errors, and atmospheric artefacts, Multi-Temporal InSAR (MT-InSAR) technologies have been developed to obtain the high precision deformation time series including Persistent Scatterers InSAR (PSInSAR) [11], Small Baseline Subset (SBAS) InSAR [12,13], Temporally Coherent Point InSAR (TCPInSAR) [14], Interferometric Point Target Analysis (IPTA) [15], Intermittent SBAS (ISBAS) [16], and Distributed Scatterers-InSAR(DS-InSAR) [17]. These advanced InSAR technologies have facilitated ground subsidence monitoring with land reclamation. For instance, Liu et al. used Sentinel-1A data to investigate the ground subsidence caused by reclamation at Xiamen Xiang'an Airport [18]. Zhao et al. used Envisat ASAR, COSMO-SkyMed, and Sentinel-1 A SAR data to investigate the surface deformation associated with reclamation at Lingang New City, Shanghai from 2007 to 2017 [4]. Wu et al. used ERS-1/2, ENVISAT ASAR, COSMO-SkyMed (CSK), and Sentinel-1A to link the ground deformation time series caused by two decades of reclamation at Hong Kong Airport [7]. Aslan et al. investigated the spatial extent and ground deformation rate related to land reclamation in the megacity of Istanbul, Turkey by using the ERS-1/2, Envisat, and Sentinel-1 datasets [19]. Shi et al. combined InSAR observations and Terzaghi consolidation theory to study the reclamation-related subsidence of Jinzhou Bay International Airport in Dalian, China [20].
For reclamation projects, temporally changed man-made structures, vegetation, and ground seepage result in the incoherence of SAR images. DS-InSAR techniques can attenuate the effects of incoherence and provide more monitoring targets. For the DS-InSAR algorithms, the SqueeSAR method was first proposed by Ferretti et al. in 2010, including homogeneous image element identification and phase optimization, and reconstructed the deformation of DS image elements [17]. In 2015, Cao et al. revealed the mathematical principles of the phase optimization process and induced two optimization methods, equal-weighted and coherent-weighted [21]. In the same year, Jiang et al. converted the hypothesis testing problem of the homogeneous image element identification process into a confidence interval finding problem, and proposed the fast homogeneous image element selection algorithm (FaSHPS), which improved the computational efficiency compared with the KS and other testing methods [22,23]. In 2019, Zhao et al. used the covariance matrix decomposition method to reduce the noise of interferograms and increase the density of coherent targets [24].
With the arrival of SAR big data era, the near real-time processing of SAR data can facilitate the development of disaster mitigation. At present, many algorithms for dynamic processing of SAR data have been developed, for example the sequential least squares algorithm [25][26][27][28], which helps us improve the efficiency of SAR data processing and reduce the consumption of storage space. In this paper, a geotechnical model is used to predict land subsidence caused by reclamation. It was first proposed by Yang et al. [29] and then applied by Zhao et al. for the prediction of future subsidence trends on reclaimed land using InSAR deformation time series [30]. In addition, the model was also used by Wu et al. to link multi-platform InSAR long-term reclaimed land deformation time series [7].
Xiamen Dadeng Island is undergoing a large-scale land reclamation project to build one of the most important airports in China, Xiamen Xiang'an Airport. At present, the first and second stages of reclamation work have been completed. The airport is planned for an area of 43 square kilometers on Dadeng Island and about 65% of the land is reclaimed [18]. The construction of Xiamen Xiang'an Airport experienced significant and uneven subsidence on the reclaimed land due to soil consolidation of the marine sediments. Zhangzhou Shuangyu Island is located in the Zhangzhou China Merchants Economic and Technological Development Zone and is the largest man-made island in China. The planned area of Shuangyu Island is approximately 221.67 hectares, 82% of which was completed by reclamation (182.30 hectares). Due to the soil consolidation of the marine sediments, significant land subsidence occurred after the completion of the reclamation.
It is therefore important to carry out the continuous ground deformation monitoring for the Shuangyu Island and the Dadeng Island Xiang'an Airport to determine the potential subsidence during their construction and future operation.
In this paper, the ground subsidence caused by land reclamation is characterized, monitored, and predicted to better understand the spatial and temporal evolution of land reclamation-induced deformation. Firstly, the SBAS-InSAR method is used to investigate the large area of subsidence in the coastal areas of Xiamen and Zhangzhou cities, Fujian Province. Then, the DS-InSAR method is considered to obtain a more complete deformation field in the typical reclaimed land areas (i.e., Shuangyu Island and Dadeng Island). Finally, the sequential estimation and geotechnical model are combined to predict the future deformation trends of the roads on Shuangyu Island.

Study Area
The study area shown in Figure 1 is located within the cities of Xiamen and Zhangzhou in southern of China, where Dadeng Island, with an area of about 13 km 2 , is located in the southeast of Xiang'an District, Xiamen city, about 25 km west of Xiamen city center and about 15 km south of Jinmen County, Taiwan, China [31]. Xiamen Xiang'an Airport, shown in Figure 2a,c, will be built in the southeast of Dadeng Island, with a total planned area of approximately 31 km 2 , where the reclamation works account for approximately 84% of the total airport area (26 km 2 ). Currently, the first and second stages of the project have been completed. The area is about 3 km 2 in the first stage project, about 7.58 km 2 in the second stage, and about 14.06 km 2 is designed in the third stage [18,32]. Shuangyu Island is located in the Zhangzhou China Merchants Economic and Technological Development Zone, southwest of Xiamen city center, in the shape of a circle surrounded by two fishes, with a radius of 840 m and a planned area of 221.67 hectares. The reclaimed land area is 182.30 hectares, and the shoreline is 11.7 km. The Shuangyu Island project is a key construction project in Fujian Province, which will be developed into an ecological island with a combination of marine tourism and holidays, special cultural and artistic leisure activities, and residences in the future.    The reclamation of Dadeng Island started in March 2014 and the second stage of land reclamation was completed in February 2017, as Google Earth's historical imagery can show us. The airport is expected to be completed by 2045 [32]. It can also be seen that the Shuangyu Island project was officially started in 2010 by comparing the optical historical images shown in Figure 2b,d. Reclamation of the entire project was completed by the end of 2017. It is clear to see from Figure 2d that the road has been constructed. The full development of Shuangyu Island will be completed by 2024 [35]. Xiamen Xiang'an Airport belongs to a coastal landform unit, where shallow silt (including silt-sand and sand-mixed silt) in the reclamation area accounts for 76% with a thickness of less than 5 m, 14% with a silt thickness of 5 to 8 m, and 10% with a silt less than 8 m. According to a survey of the silt thickness distribution, there is a deep N-S trench in the middle of the reclamation area with the maximum silt thickness of 30 m and a width of 400 m. During the reclamation process, the first stage of the reclamation area was filled with sand blowing to form the land area at an elevation of 5 m. The second stage of the reclamation area was filled with dredged mud blowing to form the land area at an elevation of 4 m. Shuangyu Island is located in a shallow inshore area, with the natural mud surface elevation ranging from −5 to −2 m in the western section and −7 to −5 m in the eastern section. The soft and weak soil layers of Shuangyu Island are mainly silt and mixed sand, and the thickness of it at the east and south sides of Shuangyu Island is less than 5 m, while the one at the central and southeast sides is from 5 to 20 m, with the maximum one being 23.5 m [33]. Shuangyu Island is filled with a combination of conventional landform materials, i.e., marine sand and dredged soil, where the northern section was directly filled with marine sand, and the southern section is blown with dredged soil [34].
The reclamation of Dadeng Island started in March 2014 and the second stage of land reclamation was completed in February 2017, as Google Earth's historical imagery can show us. The airport is expected to be completed by 2045 [32]. It can also be seen that the Shuangyu Island project was officially started in 2010 by comparing the optical historical images shown in Figure 2b,d. Reclamation of the entire project was completed by the end of 2017. It is clear to see from Figure 2d that the road has been constructed. The full development of Shuangyu Island will be completed by 2024 [35].

Data
We selected 125 ascending Sentinel-1A TOPS mode Interferometric Wide Swath (IWS) images acquired from 13 June 2017 to 21 October 2021 with VV polarization, C-band, and an incidence angle of 34 degrees. During this period, the first and second stages of the Xiamen Xiang'an Airport project on Dadeng Island and the Shuangyu Island reclamation project were largely completed. The topographic phase was also removed using 30 m SRTM DEM, where most areas are covered by water. For the processing of data, we used a multi-look ratio of 4:1 in the range and azimuth direction to obtain a spatial resolution of approximately 16 m and selected the interference pairs shown in Figure 3 to calculate the deformation rates and time series. For the dynamic prediction of ground subsidence, we employed the sequential estimation and the geotechnical model, where we divided the interferometric pairs into two groups. The first group is the archived interferometric pairs, which are used for the initialization of the deformation parameters (indicated by blue lines). The second group is the interferometric pairs generated from the newly acquired SAR images to link older archived SAR images (indicated by the red lines). As studied by Wang et al. [26], the sequential estimation can get the identical deformation time series results to the conventional SBAS method.
We selected 125 ascending Sentinel-1A TOPS mode Interferometric Wide Swath (IWS) images acquired from 13 June 2017 to 21 October 2021 with VV polarization, Cband, and an incidence angle of 34 degrees. During this period, the first and second stages of the Xiamen Xiang'an Airport project on Dadeng Island and the Shuangyu Island reclamation project were largely completed. The topographic phase was also removed using 30 m SRTM DEM, where most areas are covered by water. For the processing of data, we used a multi-look ratio of 4:1 in the range and azimuth direction to obtain a spatial resolution of approximately 16 m and selected the interference pairs shown in Figure 3 to calculate the deformation rates and time series. For the dynamic prediction of ground subsidence, we employed the sequential estimation and the geotechnical model, where we divided the interferometric pairs into two groups. The first group is the archived interferometric pairs, which are used for the initialization of the deformation parameters (indicated by blue lines). The second group is the interferometric pairs generated from the newly acquired SAR images to link older archived SAR images (indicated by the red lines). As studied by Wang et al. [26], the sequential estimation can get the identical deformation time series results to the conventional SBAS method.  Table 1 shows the SAR data coverage and techniques involved in different regions. In terms of data processing, the SBAS-InSAR method is used rather than the DS-InSAR method to obtain large scale deformation rates because of the slow processing speed and high computer storage and computational performance requirements of the DS-InSAR method. To overcome the effects of incoherence caused by vegetation and surface changes in Shuangyu and Dadeng Islands, the DS-InSAR method is used to obtain the land subsidence rate and recovered the complete deformation field of the nonlinear subsidence due to the reclamation in two typical artificial islands. In order to update the prediction model in near real time, a sequential least squares estimation technique was added to the Shuangyu Island data processing to quickly update the deformation time series dynamically.   Table 1 shows the SAR data coverage and techniques involved in different regions. In terms of data processing, the SBAS-InSAR method is used rather than the DS-InSAR method to obtain large scale deformation rates because of the slow processing speed and high computer storage and computational performance requirements of the DS-InSAR method. To overcome the effects of incoherence caused by vegetation and surface changes in Shuangyu and Dadeng Islands, the DS-InSAR method is used to obtain the land subsidence rate and recovered the complete deformation field of the nonlinear subsidence due to the reclamation in two typical artificial islands. In order to update the prediction model in near real time, a sequential least squares estimation technique was added to the Shuangyu Island data processing to quickly update the deformation time series dynamically.

Method
To investigate the spatial and temporal evolution of land deformation during the land reclamation in the coastal areas of Xiamen and Zhangzhou City, we use SBAS-InSAR and DS-InSAR techniques and remote sensing historical images to investigate the deformation procedure.
Firstly, pre-processing of the acquired sentinel-1A SAR data is conducted, which includes radiometric correction, selection of the required bursts, fine alignment to 1/1000th Remote Sens. 2022, 14, 2930 6 of 19 of a pixel, azimuthal spectral de-skewing, and bursts stitching. As for the SBAS-InSAR method, assuming N + 1 SAR images selected for the same area, interferometric pairs are generated based on a spatiotemporal baseline thresholds, and the number of interferograms M should satisfy the following conditions.
Then, interferogram filtering and unwrapping are performed. The unwrapped phase of any pixel (x, r) of the m-th differential interferogram is expressed as where t A and t B represent the acquisition date of the SAR image, respectively, corresponding to the m-th differential interferogram; δψ (de f ,x,r) is the deformation phase component in the line of sight (LOS) direction between t A and t B ; δψ (ε,x,r) is the topographic phase error; δψ (a,x,r) is the atmospheric phase error; δψ (n,x,r) is the noise phase. The atmospheric errors related to elevation are removed by polynomial fitting. Finally, the deformation time series on the high-coherence pixels are obtained by the least squares or singular value decomposition (SVD) criteria. As the land deformation in the coastal region is mainly due to compaction in the vertical direction, we convert the LOS deformation to the vertical direction with the following equation.
where d vertical is the vertical deformation; d los is the deformation in the LOS direction; θ is the incidence angle.
For the deformation monitoring of the two typical study areas (i.e., Dadeng Island and Shuangyu Island), we used the DS-InSAR method proposed by Zhao et al. for robust estimation [24], where the covariance matrix is estimated based on homogeneous points, followed by the eigenvalue decomposition and phase optimization to the covariance matrix to achieve filter the noise of the interferogram.
As for the dynamic land subsidence monitoring of Shuangyu Island, the sequential estimation proposed by Wang et al. [26][27][28]36] and the correlation model between ground subsidence and time adopted by Zhao et al. [30] are considered. As the sequential estimation does not process all SAR data repeatedly when new SAR data added, the deformation parameters can be dynamically updated without loss of accuracy, improving the computational efficiency of InSAR processing. As a result, the deformation time series can be updated rapidly along with each new acquisition of SAR image. At the same time, the best-fit geotechnical model is dynamically updated.
According to previous investigations [37,38], there are two main stages of ground subsidence associated with land reclamation, i.e., primary consolidation and long-term secondary compression. Most of the reclamation-induced deformation will last for decades.
Here, a geotechnical model is used to predict the future evolution of ground subsidence on Shuangyu Island, which shows the ground subsidence as a function of time and can be expressed by the following equation [29]: where t represents the number of days relative to the start time of 20 May 2017; L(t) represents the cumulative ground subsidence at time t; S m denotes asymptotic (total) subsidence assuming infinite time; and k and λ are the model parameters that affect the curvature of the model.

Large Area Identification of Land Subsidence
First, the deformation areas in Xiamen and Zhangzhou cities were detected. The obtained deformation rate maps (Figure 4) show that a total of 23 deformation areas were detected, 17 of which are on the coastal regions. Note that most areas were stable with a vertical deformation rate between −10 mm/year and 10 mm/year, and most areas suffered uneven ground subsidence caused by land reclamation, especially in Dadeng and Shuangyu Islands. The maximum deformation rate of the subsidence center reached 100 mm/a in Dadeng and Shuangyu Islands.
where t represents the number of days relative to the start time of 20 May 2017; L(t) represents the cumulative ground subsidence at time t; S denotes asymptotic (total) subsidence assuming infinite time; and k and λ are the model parameters that affect the curvature of the model.

Large Area Identification of Land Subsidence
First, the deformation areas in Xiamen and Zhangzhou cities were detected. The obtained deformation rate maps (Figure 4) show that a total of 23 deformation areas were detected, 17 of which are on the coastal regions. Note that most areas were stable with a vertical deformation rate between −10 mm/year and 10 mm/year, and most areas suffered uneven ground subsidence caused by land reclamation, especially in Dadeng and Shuangyu Islands. The maximum deformation rate of the subsidence center reached 100 mm/a in Dadeng and Shuangyu Islands.

Dadeng Island
We calculated the deformation rate and cumulative deformation time series in Dadeng Island as shown in Figures 5a and 6. It can be seen that there are four main subsidence zones, of which A1 is located within the future planned airstrip of Xiamen Xiang'an Airport, A2 and A3 are located in the west and northwest of Dadeng Island, and A4 is located in the north. To investigate the spatial deformation of the subsidence zones, four profiles, namely AA', BB', CC', and DD', were extracted, and shown in Figure 5b-d. The maximum vertical subsidence in the four subsidence centers exceeded 90 mm/a. As the new airport has been completed in the first and second stages, the deformation in most areas was less than 10 mm/a, indicating a stable stage. Figure 6 shows the cumulative deformation time series maps, where it is noticeable to see the increase of the magnitude

Dadeng Island
We calculated the deformation rate and cumulative deformation time series in Dadeng Island as shown in Figures 5a and 6. It can be seen that there are four main subsidence zones, of which A1 is located within the future planned airstrip of Xiamen Xiang'an Airport, A2 and A3 are located in the west and northwest of Dadeng Island, and A4 is located in the north. To investigate the spatial deformation of the subsidence zones, four profiles, namely AA', BB', CC', and DD', were extracted, and shown in Figure 5b-d. The maximum vertical subsidence in the four subsidence centers exceeded 90 mm/a. As the new airport has been completed in the first and second stages, the deformation in most areas was less than 10 mm/a, indicating a stable stage. Figure 6 shows the cumulative deformation time series maps, where it is noticeable to see the increase of the magnitude and extent of land subsidence, and the maximum cumulative subsidence was over 425 mm. and extent of land subsidence, and the maximum cumulative subsidence was over 425 mm.

Shuangyu Island
We also calculated the deformation rate of Shuangyu Island and its surrounding area, as shown in Figure 7a. There are four subsidence centers, where two centers are on Shuangyu Island and two are in the coastal area adjacent to Shuangyu Island. Verified by the historical remote sensing images, four subsidence centers were caused by land reclamation. The profiles crossing four subsidence centers (EE', FF', GG', HH') were also extracted to reveal the spatial characteristics of the land subsidence. Two areas on Shuangyu Island experienced severe ground subsidence, with the maximum vertical deformation rate of 120 mm/a. The subsidence center with profile GG' suffered the maximum vertical deformation rate of 130 mm/a.

Shuangyu Island
We also calculated the deformation rate of Shuangyu Island and its surrounding area, as shown in Figure 7a. There are four subsidence centers, where two centers are on Shuangyu Island and two are in the coastal area adjacent to Shuangyu Island. Verified by the historical remote sensing images, four subsidence centers were caused by land reclamation. The profiles crossing four subsidence centers (EE', FF', GG', HH') were also extracted to reveal the spatial characteristics of the land subsidence. Two areas on Shuangyu Island experienced severe ground subsidence, with the maximum vertical deformation rate of 120 mm/a. The subsidence center with profile GG' suffered the maximum vertical deformation rate of 130 mm/a.

Deformation Prediction of Land Subsidence at Shuangyu Island
The magnitude and velocity of land subsidence after land reclamation depends to a large extent on the type and thickness of the reclamation material, the thickness of the

Deformation Prediction of Land Subsidence at Shuangyu Island
The magnitude and velocity of land subsidence after land reclamation depends to a large extent on the type and thickness of the reclamation material, the thickness of the underlying soil materials, and the duration of time after the reclamation is completed [39]. According to previous studies [40,41], most predictions related to land subsidence are based on the consolidation theory, for example, Terzaghi theory of consolidation [37,42].
Shuangyu Island has already experienced severe land subsidence after the completion of reclamation. Figure 8a shows the cumulative vertical deformation map, and the maximum cumulative vertical deformation reached 626 mm. After the completion of the reclamation, land subsidence will take place immediately and continue for decades [37,38]. Therefore, it is necessary to predict the future development of land subsidence on the artificial island. To this end, we use the sequential estimation to dynamically update the deformation time series of the Shuangyu Island, where we mainly focus on the future deformation of the newly constructed roads. We predict the deformation of the road by using the geotechnical model [29,30], and verify it by the InSAR sequential deformation estimation [28]. underlying soil materials, and the duration of time after the reclamation is completed [39]. According to previous studies [40,41], most predictions related to land subsidence are based on the consolidation theory, for example, Terzaghi theory of consolidation [37,42]. Shuangyu Island has already experienced severe land subsidence after the completion of reclamation. Figure 8a shows the cumulative vertical deformation map, and the maximum cumulative vertical deformation reached 626 mm. After the completion of the reclamation, land subsidence will take place immediately and continue for decades [37,38]. Therefore, it is necessary to predict the future development of land subsidence on the artificial island. To this end, we use the sequential estimation to dynamically update the deformation time series of the Shuangyu Island, where we mainly focus on the future deformation of the newly constructed roads. We predict the deformation of the road by using the geotechnical model [29,30], and verify it by the InSAR sequential deformation estimation [28]. We use the best-fit curve model to predict the deformation at a given time, where three parameters , k, and λ of Equation (4) are estimated for each pixel. Then, the deformation time series is updated with the sequential estimation at the given time. Accordingly, the best-fit curve model is updated iteratively. As shown in Figure 9, we divided the data into two groups, i.e., archived SAR data are from 13 June 2017 to 11 January 2021 and newly acquired SAR data are from 13 January 2021 to 20 September 2021. The prediction model is iteratively updated as more observations are added. The 20 new observations matched the archived prediction model accurately. As shown in Figure 10, the standard deviation of the difference between the archived model predictions and the true values calculated from the sequential estimates at four points is within 3.7 mm, which greatly verifies the effectiveness of the geotechnical model. Therefore, the updated prediction model can be applied to predict the future land subsidence trend. We use the best-fit curve model to predict the deformation at a given time, where three parameters S m , k, and λ of Equation (4) are estimated for each pixel. Then, the deformation time series is updated with the sequential estimation at the given time. Accordingly, the best-fit curve model is updated iteratively. As shown in Figure 9, we divided the data into two groups, i.e., archived SAR data are from 13 June 2017 to 11 January 2021 and newly acquired SAR data are from 13 January 2021 to 20 September 2021. The prediction model is iteratively updated as more observations are added. The 20 new observations matched the archived prediction model accurately. As shown in Figure 10, the standard deviation of the difference between the archived model predictions and the true values calculated from the sequential estimates at four points is within 3.7 mm, which greatly verifies the effectiveness of the geotechnical model. Therefore, the updated prediction model can be applied to predict the future land subsidence trend. Figure 9. The geotechnical model predicted vertical deformation time series at four points S1 to S4 (a-d), whose positions are marked in Figure 8b. Figure 10. The difference between the prediction by archived deformation time sereis and the estimation by the sequential method at four points S1 to S4 (a-d).
Therefore, we make the preliminary prediction of the cumulative vertical deformation of the roads from 2021 to 2026 by using an updated prediction model calculated from the first SAR acquisition date on 13 June 2017. In Figure 11, the vertical deformation time series maps from 13 June 2017 to 3 August 2021 were obtained by using sequential estimation method, while the ones from 30 October 2021 to 18 March 2026 shown in the black box are predicted based on the updated prediction model. The quality of the nonlinear curve fitting was verified by the standard deviation of the difference between the InSAR time series and the best-fitting model. The average standard deviation was 4.5 mm, indicating the high accuracy of the prediction model. The results show that the maximum cumulative vertical deformation of the road on Shuangyu Island will exceed 1.1 m by 2026. Figure 9. The geotechnical model predicted vertical deformation time series at four points S1 to S4 (a-d), whose positions are marked in Figure 8b.
Remote Sens. 2022, 14, 2930 12 of 18 Figure 9. The geotechnical model predicted vertical deformation time series at four points S1 to S4 (a-d), whose positions are marked in Figure 8b. Figure 10. The difference between the prediction by archived deformation time sereis and the estimation by the sequential method at four points S1 to S4 (a-d).
Therefore, we make the preliminary prediction of the cumulative vertical deformation of the roads from 2021 to 2026 by using an updated prediction model calculated from the first SAR acquisition date on 13 June 2017. In Figure 11, the vertical deformation time series maps from 13 June 2017 to 3 August 2021 were obtained by using sequential estimation method, while the ones from 30 October 2021 to 18 March 2026 shown in the black box are predicted based on the updated prediction model. The quality of the nonlinear curve fitting was verified by the standard deviation of the difference between the InSAR time series and the best-fitting model. The average standard deviation was 4.5 mm, indicating the high accuracy of the prediction model. The results show that the maximum cumulative vertical deformation of the road on Shuangyu Island will exceed 1.1 m by 2026. Figure 10. The difference between the prediction by archived deformation time sereis and the estimation by the sequential method at four points S1 to S4 (a-d).
Therefore, we make the preliminary prediction of the cumulative vertical deformation of the roads from 2021 to 2026 by using an updated prediction model calculated from the first SAR acquisition date on 13 June 2017. In Figure 11, the vertical deformation time series maps from 13 June 2017 to 3 August 2021 were obtained by using sequential estimation method, while the ones from 30 October 2021 to 18 March 2026 shown in the black box are predicted based on the updated prediction model. The quality of the non-linear curve fitting was verified by the standard deviation of the difference between the InSAR time series and the best-fitting model. The average standard deviation was 4.5 mm, indicating the high accuracy of the prediction model. The results show that the maximum cumulative vertical deformation of the road on Shuangyu Island will exceed 1.1 m by 2026.

Spatial Evolution of Artificial Islands
Ten-year historical remote sensing images from 2011 to 2021 of Xiamen Dadeng Island and Zhangzhou Shuangyu Island were acquired. The reclamation process of Shuangyu Island is shown in Figure 12a-d, where it can be divided into two main stages. The first stage is shown with the white line and was substantially completed on 1 January 2012. The second stage is shown with the yellow line and was largely completed on 12 January 2014. By now, the roads on the island have been constructed. Meanwhile, the reclamation process of the Xiamen Xiang'an Airport on Dadeng Island is shown in Figure  12e

Spatial Evolution of Artificial Islands
Ten-year historical remote sensing images from 2011 to 2021 of Xiamen Dadeng Island and Zhangzhou Shuangyu Island were acquired. The reclamation process of Shuangyu Island is shown in Figure 12a-d, where it can be divided into two main stages. The first stage is shown with the white line and was substantially completed on 1 January 2012. The second stage is shown with the yellow line and was largely completed on 12 January 2014. By now, the roads on the island have been constructed. Meanwhile, the reclamation process of the Xiamen Xiang'an Airport on Dadeng Island is shown in Figure 12e

Secondary Landfill and Land Subsidence
As investigated, two land subsidence centers in the western and northern p Dadeng Island, Xiamen, were first filled in 2017. Figure 13 shows four optical remot ing images acquired on 18 June 2019 (a), 30 April 2020 (b), 26 October 2020 (e), and 2019 (f) and two deformation rate maps (c) and (d) from May 2017 to Novembe where the maximum vertical deformation rate reached 120 mm/a. The two reclaim eas show significant seepage via (a) and (b), indicating that they are still in an u state. We can also see from Figure 13e,f that these two areas were filled secondl 2019 to 2020. Figure 14 shows the deformation time series in the vertical direction at four points, P1 to P4, after the secondary landfilling, where the locations of four poi shown in Figure 13a,f, i.e., P1 and P2 are in area A and P3 and P4 are in area B. The mation of four points began to accelerate after April 2020, which can be explained following three possible reasons: the consolidation of the new filled material, the c

Secondary Landfill and Land Subsidence
As investigated, two land subsidence centers in the western and northern parts of Dadeng Island, Xiamen, were first filled in 2017. Figure 13 Figure 14 shows the deformation time series in the vertical direction at four typical points, P1 to P4, after the secondary landfilling, where the locations of four points are shown in Figure 13a,f, i.e., P1 and P2 are in area A and P3 and P4 are in area B. The deformation of four points began to accelerate after April 2020, which can be explained as the following three possible reasons: the consolidation of the new filled material, the consolidation of the previous filled material, and the compaction of the unconsolidated marine sediment after secondary land reclamation. Figure 15 shows the deformation time series in the vertical direction at two typical points P5 and P6 after the secondary landfilling of Xiamen Xiang'an Airport, which are shown in the last frame of Figure 6. Green lines in Figure 15a,b indicate the fitting lines calculated by Equation (3). It can be seen that the same acceleration phenomenon also occurred at Xiamen Xiang'an Airport.

Analysis of the Causes of Uneven Subsidence on Reclaimed Land
According to the seabed geological conditions of Shuangyu Island, the thickness of it at the east and south sides of Shuangyu Island is less than 5 m, while the one at the central and southeast sides is from 5 to 20 m, with the maximum one as 23.5 m [33]. Figure 16 is a correlation of the fill elevation and cumulative deformation. It can be seen that there are two uplifted hills on the southeast side of Shuangyu Island, with elevations of 18.5 m and 21.5 m. Two parts of the northwest side are uplifted to hills with elevations of 12.5 and 8.2 m [34].
The cumulative subsidence map shows that only two deformation centers are formed on the southeast side, and not on the northwest side. It can be speculated that the thickness of the soft seabed soil should be the main factor that causes the surface subsidence, and the reclamation elevation is the secondary factor. A similar situation exists in the reclamation area of Xiamen Xiang'an Airport on Dadeng Island, where the subsidence center is a deep ditch area with a thickness of 30 m of soft soil.    Figure 15 shows the deformation time series in the vertical direction at two typical points P5 and P6 after the secondary landfilling of Xiamen Xiang'an Airport, which are shown in the last frame of Figure 6. Green lines in Figure 15a,b indicate the fitting lines calculated by Equation (3). It can be seen that the same acceleration phenomenon also occurred at Xiamen Xiang'an Airport.

Analysis of the Causes of Uneven Subsidence on Reclaimed Land
According to the seabed geological conditions of Shuangyu Island, the thickness of it at the east and south sides of Shuangyu Island is less than 5 m, while the one at the central and southeast sides is from 5 to 20 m, with the maximum one as 23.5 m [33]. Figure 16 is a correlation of the fill elevation and cumulative deformation. It can be seen that there are two uplifted hills on the southeast side of Shuangyu Island, with elevations of 18.5 m and 21.5 m. Two parts of the northwest side are uplifted to hills with elevations of 12.5 and 8.2 m [34]. The cumulative subsidence map shows that only two deformation centers are formed on the southeast side, and not on the northwest side. It can be speculated that the thickness of the soft seabed soil should be the main factor that causes the surface subsidence, and the reclamation elevation is the secondary factor. A similar situation exists in the reclamation area of Xiamen Xiang'an Airport on Dadeng Island, where the subsidence center is a deep ditch area with a thickness of 30 m of soft soil. Xiang'an Airport, whose positions are marked in Figure 6.
Remote Sens. 2022, 14, 2930 16 of 18 Figure 15 shows the deformation time series in the vertical direction at two typical points P5 and P6 after the secondary landfilling of Xiamen Xiang'an Airport, which are shown in the last frame of Figure 6. Green lines in Figure 15a,b indicate the fitting lines calculated by Equation (3). It can be seen that the same acceleration phenomenon also occurred at Xiamen Xiang'an Airport.

Analysis of the Causes of Uneven Subsidence on Reclaimed Land
According to the seabed geological conditions of Shuangyu Island, the thickness of it at the east and south sides of Shuangyu Island is less than 5 m, while the one at the central and southeast sides is from 5 to 20 m, with the maximum one as 23.5 m [33]. Figure 16 is a correlation of the fill elevation and cumulative deformation. It can be seen that there are two uplifted hills on the southeast side of Shuangyu Island, with elevations of 18.5 m and 21.5 m. Two parts of the northwest side are uplifted to hills with elevations of 12.5 and 8.2 m [34]. The cumulative subsidence map shows that only two deformation centers are formed on the southeast side, and not on the northwest side. It can be speculated that the thickness of the soft seabed soil should be the main factor that causes the surface subsidence, and the reclamation elevation is the secondary factor. A similar situation exists in the reclamation area of Xiamen Xiang'an Airport on Dadeng Island, where the subsidence center is a deep ditch area with a thickness of 30 m of soft soil.   [34]. Background is the cumulative subsidence of Shuangyu Island.

Conclusions
In this paper, the SBAS-InSAR and DS-InSAR methods are used to detect the land subsidence areas along the coast of Zhangzhou City and Xiamen City, China, and to monitor the spatiotemporal land subsidence evolution of reclaimed lands with Sentinel-1 SAR images from 13 June 2017 to 20 September 2021. A total of 23 deformation areas were detected, of which 17 land subsidence areas are on the coastal regions including two typical deformation areas, i.e., Dadeng Island and Shuangyu Island.
The maximum vertical deformation rate of Dadeng Island is 120 mm/a, and it reaches 130 mm/a for Shuangyu Island, which is mainly correlated with the land reclamation. Specifically, the secondary reclamation accelerates the land subsidence significantly.
The sequential estimation method is applied to update the land subsidence dynamically and the geotechnical model is used iteratively to predict the deformation in the near future. Results show that the maximum cumulative deformation will reach 1.1 m by 2026 alone. Accordingly, special measures should be taken for the safety of infrastructure construction.
These research results will have important guiding significance for the design of Xiamen Xiang'an Airport and Shuangyu Island, the next stage of construction, and the prevention and control of disasters.