remote sensing Spatially Varying Relationships between Land Subsidence and Urbanization: A Case Study in Wuhan, China

: Land subsidence has become an increasing global concern over the past few decades due to natural and anthropogenic factors. However, although several studies have examined factors affecting land subsidence in recent years, few have focused on the spatial heterogeneity of relationships between land subsidence and urbanization. In this paper, we adopted the small baseline subset-synthetic aperture radar interferometry (SBAS-InSAR) method using Sentinel-1 radar satellite images to map land subsidence from 2015 to 2018 and characterized its spatial pattern in Wuhan. The bivariate Moran’s I index was used to test and visualize the spatial correlations between land subsidence and urbanization. A geographically weighted regression (GWR) model was employed to explore the strengths and directions of impacts of urbanization on land subsidence. Our ﬁndings showed that land subsidence was obvious and unevenly distributed in the study area, the annual deformation rate varied from − 42.85 mm/year to +29.98 mm/year, and its average value was − 1.0 mm/year. A clear spatial pattern for land subsidence in Wuhan was mapped, and several apparent subsidence funnels were primarily located in central urban areas. All urbanization indicators were found to be signiﬁcantly spatially correlated with land subsidence at different scales. In addition, the GWR model results showed that all urbanization indicators were signiﬁcantly associated with land subsidence across the whole study area in Wuhan. The results of bivariate Moran’s I and GWR results conﬁrmed that the relationships between land subsidence and urbanization spatially varied in Wuhan at multiple spatial scales. Although scale dependence existed in both the bivariate Moran’s I and GWR models for land subsidence and urbanization indicators, a “best” spatial scale could not be conﬁrmed because the disturbance of factors varied over different sampling scales. The results can advance the understanding of the relationships between land subsidence and urbanization, and they will provide guidance for subsidence control and sustainable urban planning.


Introduction
Land subsidence (LS) is the gentle settlement of the ground surface due to the consolidation of compressible sediments or loss of regional earth materials as a result of water exploitation or extraction of oil and gas. Land subsidence is usually observed as a series of geological-environmental hazards, including severe destruction of buildings [1], roads [2], bridges, pipelines [3], railway tracks [4], and metro networks [5], and it increases the risk of urban flooding particularly in coastal regions experiencing sea-level rise [6]. In 1891, land subsidence was first observed in Mexico City while it has been detected and recorded in Shanghai (China) since 1921 [7,8]. Currently, land subsidence occurs mainly in regions with flat terrains where loose deposits accumulate in river deltas or coastal plains, especially densely populated areas, as well as urban or agricultural areas developed in temperate or arid climates that are characterized by long-term drought [9][10][11]. Although land subsidence affects 8% of the global terrestrial area, people at risk account for approximately 16% of the total world population comprising approximately 1.2 billion inhabitants [12]. Along with a growing population and the developments in urbanization and industrialization, the exposure of the population to land subsidence is expected to increase by 30% approximately to 1.6 billion inhabitants while predicted potential subsidence areas will increase by only 7% globally; 1596 of 7343 major world cities are predicted to be located in potential subsidence areas by 2040 [13]. Land subsidence is one of the considerable challenges facing the world, and poses a significant threat to the long-term sustainable development of the humankind. Therefore, it is crucial to effectively monitor and map land subsidence in real time [14].
Differential synthetic aperture radar interferometry (D-InSAR) spaceborne-based is a remotely sensed technology that enables investigation of widespread surface deformation across the earth. In contrast to conventional spirit levelling and GPS survey techniques, D-InSAR provides spatially dense displacement measurements that are updated periodically at relatively low cost. Advanced multitemporal InSAR (MT-InSAR) technology such as persistent scatterer interferometry (PS-InSAR) [15,16] and small baseline subset interferometry (SBAS-InSAR) [17,18] can effectively overcome spatial-temporal decorrelations and mitigate atmospheric delay effects, thus allowing the measurement of surface deformation with centimeter and even subcentimeter accuracy at very high spatial resolution [19]. To date, the MT-InSAR method has been widely applied to monitor deformation of the earth's surface [20], including volcanoes and seismic activity [21], landslides [22], glacial motion [23], mining-related subsidence [24], subsidence of urban or peri-urban area [25,26], and even large-scale land deformation nationwide [27].
Wuhan, as the largest city in central China, has recently been a relatively highly populous metropolis. In the last two decades, Wuhan has experienced rapid urbanization and industrialization combined with ongoing construction of high-rise buildings and massive underground space development. At present, a variety of efforts have been made to measure and map land subsidence derived from multisource space-borne remotely sensed data based on the MT-InSAR method in Wuhan, concentrating on wide-coverage urban surface settlement and building structure stability. It can be stated in detail as follows. Zhou monitored the spatial distribution of land subsidence from 2015 to 2016 using Sentinel-1 SAR images [28]. Han and Jiang [29,30] investigated the spatial pattern and temporal evolution characteristics of land subsidence, aiming to reveal the spatiotemporal variations in subsidence and the induction of subsidence. Bai and Zhang [31,32] qualitatively analyzed the relationship between land subsidence and influential factors from the aspects of natural conditions and human activities. In addition, Ding [33] predicted time-series surface subsidence based on long short-term memory (LSTM) model in selected key regions.
The aforementioned works have led to a substantial understanding of the spatial extent, magnitude, and temporal evolution of land deformation and qualitatively explored influential factors in Wuhan. However, local clustering patterns between land subsidence and urbanization have been ignored, and fewer studies have examined the quantitative impact of urbanization on land subsidence in Wuhan, especially for spatially varying relationships between land subsidence and urbanization at different scales. In fact, knowledge of spatially varying relationships is an important prerequisite to protect cities from damage due to surface subsidence. In addition, previous studies have neglected the interactions between land subsidence and its various impact factors, which exhibit scale dependence.
To address the gaps in the existing studies, we first adopted the SBAS-InSAR method using Sentinel-1 radar satellite images to map land subsidence and characterized its spatial pattern in Wuhan. Furthermore, the bivariate Moran's I index was used to test and visualize the spatial correlations between land subsidence and urbanization. In addition, a geographically weighted regression (GWR) model was employed to explore the strengths and directions of impacts of urbanization on land subsidence. Finally, we discussed scale effects on the spatially varying relationship between land subsidence and urbanization.

Study Area
Wuhan is located in the eastern region of the Jiang-Han Plain, with geographical coordinates between latitudes 29 • 58 N and 31 • 22 N and between longitudes 113 • 41 E and 115 • 05 E (see Figure 1). The Yangtze River and its largest branch, the Han River, have a confluence in central urban areas of Wuhan, which divides it into three main parts. The overall terrain of Wuhan is flat, with an average altitude of approximately 37 m above sea level. Soft clay layers with high compressibility and low strength are primarily spread along both banks of the Yangtze River in Wuhan. The belts of carbonate rock in Wuhan cover an area of more than 1100 km 2 , and they align in a trend with a nearly east-west orientation. Wuhan is in a northern humid subtropical monsoon climate zone, which is characterized by plentiful rainfall and abundant sunshine. The average annual temperature is approximately 16.6 • C and precipitation ranges from 1150 mm to 1450 mm. The rainfall is concentrated mainly in the rainy season from June to August every year, accounting for 41% of the total precipitation. visualize the spatial correlations between land subsidence and urbanization. In addition, a geographically weighted regression (GWR) model was employed to explore the strengths and directions of impacts of urbanization on land subsidence. Finally, we discussed scale effects on the spatially varying relationship between land subsidence and urbanization.

Study Area
Wuhan is located in the eastern region of the Jiang-Han Plain, with geographical coordinates between latitudes 29°58′ N and 31°22′ N and between longitudes 113°41′ E and 115°05′ E (see Figure 1). The Yangtze River and its largest branch, the Han River, have a confluence in central urban areas of Wuhan, which divides it into three main parts. The overall terrain of Wuhan is flat, with an average altitude of approximately 37 m above sea level. Soft clay layers with high compressibility and low strength are primarily spread along both banks of the Yangtze River in Wuhan. The belts of carbonate rock in Wuhan cover an area of more than 1100 km , and they align in a trend with a nearly east-west orientation. Wuhan is in a northern humid subtropical monsoon climate zone, which is characterized by plentiful rainfall and abundant sunshine. The average annual temperature is approximately 16.6 °C and precipitation ranges from 1150 mm to 1450 mm. The rainfall is concentrated mainly in the rainy season from June to August every year, accounting for 41% of the total precipitation. Wuhan, the capital city of Hubei Province, acts as an important industrial and economic center, cultural and educational base, and comprehensive transportation hub in central China. It covers a total area of 8494 km with permanent residents totaling approximately 11.08 million in 2018. Wuhan has witnessed unprecedented economic development and urban sprawl since the new millennium.The urban built-up area in 2000 was 221.01 km , while it reached 812.39 km at the end of 2018, and the average annual urban expansion rate was more than 7%. During the process of horizontal urban growth, a mass of high-rise buildings and underground space development promoted significant urban Wuhan, the capital city of Hubei Province, acts as an important industrial and economic center, cultural and educational base, and comprehensive transportation hub in central China. It covers a total area of 8494 km 2 with permanent residents totaling approximately 11.08 million in 2018. Wuhan has witnessed unprecedented economic development and urban sprawl since the new millennium.The urban built-up area in 2000 was 221.01 km 2 , while it reached 812.39 km 2 at the end of 2018, and the average annual urban expansion rate was more than 7%. During the process of horizontal urban growth, a mass of high-rise buildings and underground space development promoted significant urban expansion in the vertical direction, thus resulting in the continuous emergence of the subsidence phenomena in Wuhan.

Land Subsidence Extraction
Sentinel-1 mission SAR sensors operate in the C-band (a wavelength of approximately 5.6 cm) and two satellites consisting of Sentinel-1 A/B satellites observe earth's surface globally within 6 days and 12 days for a single satellite. The Sentinel-1 interferometric wide swath mode (IW) provides large swath width of~250 km images using the novel terrain observation by progressive scans (TOPSAR) imaging technique. The spatial resolutions of the IW mode are less than 20 m and 5 m in the azimuth and range directions (single look), respectively [34]. The Sentinel-1 satellite single look complex (SLC) data for interferometric applications can be accessed freely from European Space Agency (ESA) Copernicus Open Access Hub website (https://scihub.copernicus.eu (accessed on 8 December 2021).
In this work, we used 30 Sentinel-1A SAR images acquired in IW mode with VV (vertical-vertical) polarization between April 2015 and January 2019 to provide an assessment of land subsidence in Wuhan, China. All scenes were acquired along the descending orbits and average value of incidence angle is approximately 41.59 • . The SBAS algorithm [17] was employed to process multitemporal IW SLC level-1 data products to derive the land subsidence velocity. A combination of time-series images within thresholds of spatial and temporal baseline (smaller than 148 m and less than 365 days) was selected to generate a connected graph of differential interferograms, which allows maximization of geometric coherence [35]. After the removal of interferometric pairs with low coherence and poor unwrapping, a total of 99 differential interferograms were obtained. The average number of connections per scene is more than 5 to ensure sufficient interconnected redundancy. The flat-earth phase of interferograms can be determined and removed by the precise orbit determination (POD) data provided by the ESA. Shuttle Radar Topographic Mission (SRTM) DEM data with a resolution of 30 m obtained from the U.S. Geological Survey (https://lta.cr.usgs.gov (accessed on 8 December 2021) was used to simulate and eliminate topographic phases and geocode displacement results. We selected 30 stable pixels without displacement located in the study site as ground control points (GCPs) to perform orbital refinement and phase reflattening for interferometric pairs. A multilooking operation with a ratio of 1:4 in the azimuth and range directions was carried out to improve the phase performance of the differential interferograms, and the interferograms were processed with the adaptive Goldstein-Werner filter to further mitigate the effects of speckle noise. Then, phase unwrapping of each interferometric pair was implemented with minimum cost flow (MCF) network and Delaunay 3D method, setting the unwrapping coherence threshold ranging from 0.2 to 0.3 by trial and error. The singular value decomposition (SVD) method was employed to generate the minimum norm least square solution for the unwrapped phase for pixels exhibiting consistently high coherence levels from interferograms and to retrieve the deformation time series. In addition, the atmospheric phase signal and nonlinear displacement component were estimated and subtracted from the displacement time series through a low pass spatial filter combined with a high temporal pass filter. Finally, the line-of-sight (LOS) deformation was transformed into the vertical direction using the radar incident angle assuming that the displacement in the horizontal direction is negligible. Thus, positive values of the deformation rate indicate that the ground is moving upwards in the vertical direction (uplift), whereas negative values mean that the ground is moving downwards in the vertical direction (subsidence). The land subsidence velocity map was extracted using SARScape module version 5.2.1 in the ENVI software environment and is shown in Figure 2.

Urbanization Metric Quantification
Impervious surfaces, defined as artificial structures that prevent natural infiltration of water into the soil, are considered an indicator of urbanization [36,37]. The impervious surface area (ISA) data were extracted from the global artificial impervious area (GAIA) dataset [38], an annual product in raster format with a 30 × 30 m resolution between 1985 and 2018. The original ISA raster data in 2018 were resampled at a resolution of 500 × 500 m. Building and road network data in vector format were obtained from the Wuhan Nature Resource and Planning Bureau.
Night-time lights generated by anthropogenic activities correlate significantly with numerous urbanization and socioeconomic parameters at regional or global scales, which have been recorded by satellite sensors for a long time. Satellite-based artificial night-time light (NTL) observations provide a unique proxy measure for unveiling urbanization and regional development [39]. Night-time light satellite images are obtained from the extended time series (2000-2018) of global NPP-VIIRS-like night-time light data [40], which has a consistent temporal trend at both global and regional scales. NTL and ISA can be applied to represent the comprehensive degree of urbanization because they both belong to physical quantity.

Urbanization Metric Quantification
Impervious surfaces, defined as artificial structures that prevent natural infiltration of water into the soil, are considered an indicator of urbanization [36,37]. The impervious surface area (ISA) data were extracted from the global artificial impervious area (GAIA) dataset [38], an annual product in raster format with a 30 × 30 m resolution between 1985 and 2018. The original ISA raster data in 2018 were resampled at a resolution of 500 × 500 m. Building and road network data in vector format were obtained from the Wuhan Nature Resource and Planning Bureau.
Night-time lights generated by anthropogenic activities correlate significantly with numerous urbanization and socioeconomic parameters at regional or global scales, which have been recorded by satellite sensors for a long time. Satellite-based artificial nighttime light (NTL) observations provide a unique proxy measure for unveiling urbanization and regional development [39]. Night-time light satellite images are obtained from the extended time series (2000-2018) of global NPP-VIIRS-like night-time light data [40], which has a consistent temporal trend at both global and regional scales. NTL and ISA can be applied to represent the comprehensive degree of urbanization because they both belong to physical quantity.
Numerous high-rise buildings and roads densely concentrated within the limits of plane space transform the natural landscape, which is necessary content and one of the spatial manifestations of urbanization. Buildings and roads parallel the intensity of urbanization, and the changes in landscape characteristics can reflect the degree of human influences on the environment. The building load plays an important factor in land subsidence. The central area of the building group has larger subsidence and the subsidence superimposition effect is obvious [41]. The kernel density method is employed to quantify the component of urbanization related to buildings, called building kernel density (BKD), based on the weight of base areas by the number of floors. Similarly, the line density method is used to estimate the component urbanization related to roads named road line density (RLD), based on the weight of the type and grade of roads.

Exploratory Spatial Data Analysis
Tobler's first law of geography [42] pointed out that ubiquitous spatial dependence occurs widely in geographical phenomena. Spatial dependency is defined as an effect between the occurrence of a given geographical location and that of surrounding locations. The global and local Moran's I indicators were used to describe the global spatial dependence among variables. The global Moran's I statistic is defined as follows [43]: where n is equal to the total number of spatial units in the study area; X i and X j are the observed value of the variable for spatial units i and j (i = j), respectively; X denotes the mean of the variable; and W ij is the spatial weight between spatial units i and j defined by the inverse distance method, which is commonly used in row-standardized form. In general, the global Moran's I value ranges from −1 to 1. The Global Moran's I > 0 indicates that similar subsidence values are clustered together (positive spatial autocorrelation), whereas the global Moran's I < 0 indicates that dissimilar subsidence values are clustered together (negative spatial autocorrelation); and when the global Moran's I is zero, no spatial autocorrelation exists. The local indicator of spatial association (LISA; local Moran's I statistic) measures the degree of spatial autocorrelation in each sample unit. For each spatial unit i, the LISA is calculated as follows [44,45]: The LISA index can identify two spatial cluster types: a high-high cluster indicating a high value surrounded by higher value; and a low-low cluster indicating a low subsidence value surrounded by neighbors with lower values. Spatial outliers refer to those values that are significantly different from the values of their neighbors, including low-high (a low value surrounded by high value) and high-low (a high value surrounded by low values) outliers. The Monte Carlo simulation method (999 permutations) was used to test the statistical significance of Moran's I, and significance value for spatial autocorrelation was set at p < 0.05. The global Moran's I is regarded as the average LISA value of all spatial units. The land subsidence/uplift clusters were obtained using GeoDa software version 1.12 (GeoDa Press LLC, Chicago, IL, USA) and are shown in

Spatial Patterns of Urbanization
As shown in Figure 4, all urbanization indicators (ISA, NTL, BKD and RLD) gradually decrease from the city center to its outer periphery, and their spatial distributions are similar to each other globally. The ratio of the area for ISA across the study area is 44.40%, and the averages of NTL, BKD, and RLD are 16.31 n·W·cm · sr , 523.40, and 14.93, respectively. ISA and NTL are comprehensive indicators of urbanization, while BKD and RLD focus on a single aspect of urbanization. Therefore, the intragroup similarity between ISA and NTL, BKD and RDL is higher than the intergroup similarity in terms of spatial patterns. It is observed that water bodies, which are mainly composed of rivers, lakes, and reservoirs, play a crucial role in shaping the urban expansion of Wuhan according to the spatial distribution of urbanization indicators. In addition, BKD and RLD were both found to exhibit low values compared to the relatively values of high ISA and NTL in the southeastern part of the study area, which indicated local spatial disparities among urbanization indicators. The urban-rural gradient of NTL and BKD has a significant spatial variation as opposed to ISA and RLD. Although the spatial distribution of RLD is different from that of the original road, RLD also showed an obvious linear distribution characteristic. Local indicator of spatial association for land subsidence at block scales of 0.5 km × 0.5 km (a), 1 km × 1 km (b), 1.5 km × 1.5 km (c), and 2 km × 2 km (d) across the study area in Wuhan city.
We use the bivariate Moran's I index to measure the spatial association relationships between land subsidence and urbanization. The bivariate global Moran's I can detect whether a cluster or outlier exist in the study area, and bivariate local Moran's I is able to identify the exact location. The bivariate global Moran's I statistic is expressed as follows [46]: For each spatial unit p, the bivariate local Moran's I (bivariate LISA) is defined as follows: where n is the total number of spatial units; z p l is the standard value of land subsidence in spatial unit p; z q m is the standard value of the urbanization metric in spatial unit q; and W pq is the spatial weight between units p and q. The value of bivariate Moran's I ranges between −1 and 1, where a positive value suggests positive spatial correlation and a negative value indicates negative spatial correlation between two variables. When bivariate Moran's I is equal to 0, it signifies a random spatial pattern. In the case of statistical significance, the bivariate local Moran's I index divides the spatial relationship between land subsidence and urbanization in each sample unit into: "High-High (HH)", "Low-Low (LL)", "High-Low (HL)", and "Low-High (LH)". We used Monte Carlo randomization (9999 permutations) to assess the significance of the bivariate Moran's I. When the test is significant (p < 0.05), there is a clustered or dispersed pattern between two variables. The land subsidence/uplift clusters were obtained using GeoDa software version 1.12 (GeoDa Press LLC, Chicago, IL, USA; Anselin, Luc, 2006) and are shown in Figure 3.

Geographically Weighted Regression Model
GWR is a relatively simple but effective, technique that extends the traditional regression framework for exploring spatial nonstationarity. It allows different relationships to exist at different points in space, such that local rather than global parameters can be estimated. Brunsdon [47] described spatial heterogeneity as a condition where a global regression model cannot describe the relationship between the response variable and explanatory variables because of the variation in characteristics among the observation regions. Global regression models, such as ordinary least squares (OLS), assume constant relationships over space, ignoring the effects of spatial heterogeneity among the observations. The GWR model can capture spatial nonstationarity by allowing the variation in relationships across space, and the model can be defined as follows [48]: where y i represents the value of the dependent variable, β i0 is the constant term, x ik is the value of the independent variable k of unit I, β ik is the parameter estimate associated with x ik , and ε i is the random error. The local estimates for unit i using matrix representation are calculated as follows [47,48]:β where X is a (n × (k + 1)) independent variable matrix, the first column of which represents the intercept term and all of them are set as 1; Y denotes an n × 1 vector of dependent variables; W(i) is an n × n matrix with the element W ij indicating the spatial weight between units i and j, whose diagonal elements are spatial weights between two units, and the off-diagonal elements are set to zero. To obtain weights, we used the adaptive Gaussian function to define the spatial kernel: where b i represents the spatial distance between units i and j, i.e., the bandwidth of unit i, which determines whether the kernel function will be performed. The selection of bandwidth of spatial kernel function plays a critical role in GWR model performance [49], which is more important than the choice of spatial kernel function itself. In this work, we used the corrected Akaike information criterion (AICc) to determine the appropriate bandwidth of each kernel as it achieves a balance between goodness-of-fit and model complexity. Under the circumstance of GWR, the AICc is expressed as follows: AIC = 2n log e (σ 2 ) + n log e (2π) + n n + tr(S) n − 2 − tr(S) , where S denotes the hat matrix;σ 2 is defined as the variance in the error term; and tr(S) is the trace of the hat matrix.

Spatial Autocorrelations of Land Subsidence
A total of 2,013,600 pixels were ultimately identified as coherent targets (CTs) in the study area, the density of which was 1293 CTs per km 2 . In general, CT pixels were more densely distributed in the core of urban areas with plenty of buildings and roads than peri-urban areas. Figure 2 demonstrates the spatial distribution of the CT pixels and the corresponding vertical deformation velocity field map across the study area. During the whole observation period, the annual deformation rate of CT pixels varied from −42.85 mm/year to +29.98 mm/year, and its average value and standard deviation were −1.0 mm/year and 3.86 mm/year, respectively. As shown in Figure 2 The global Moran's I index was calculated to examine the spatial dependence of land subsidence across the whole study area (Figure 3). The value of global Moran's I at the four block scales ranged from 0.70 to 0.81 and passed the significance test at the 99% confidence level, which indicated a significantly positive global spatial autocorrelation of land subsidence across the entire study area. In general, the values of the global Moran's index for land subsidence have an increasing trend with increasing block scale, with the exception of 1.5 km × 1.5 km. The multiscale comparative analysis of the global Moran's I index suggests that the spatial autocorrelation of land subsidence is not an accidental phenomenon dependent on scale.
The LISA index was used to depict the local spatial correlation of land subsidence across the study area ( Figure 3). The cluster maps of LISA in Figure 3 exhibit the spatial aggregation state of land subsidence and distinct spatial patterns of the clusters at four block scales. The low-low clusters (serious land subsidence) were mainly concentrated in QK, JH, the south of Hanyang (HY), the southwest of HS, the north of Jiangxia (JX) and WC along the southern bank of the Yangtze River, while the high-high value agglomerations (land uplift) were mostly located north of HS and south of Huangpi (HP) along the Yangtze River. With an increase in block scale, the area of low-value and high-value agglomerations for land subsidence gradually shrank and some clusters disappeared. Nevertheless, the spatial distribution of low-value and high-value agglomerations for land subsidence became more concentrated. The LISA index can accurately delineate the funnel of land subsidence quantitatively. Thus, the multiscale LISA index can help guide corresponding actions at different administrative levels to address the consequences of land subsidence. The small number of spatial outliers (low-high or high-low) for land subsidence at four block scales were sporadically distributed within the study area. The high-low outlier refers to area where a unit with higher subsidence value is surrounded by other comparatively lower subsidence value units, it is likely to develop into an emerging subsiding area in the near future according to Tobler's first law of geography [42].

Spatial Patterns of Urbanization
As shown in Figure 4, all urbanization indicators (ISA, NTL, BKD and RLD) gradually decrease from the city center to its outer periphery, and their spatial distributions are similar to each other globally. The ratio of the area for ISA across the study area is 44.40%, and the averages of NTL, BKD, and RLD are 16.31 n·W·cm −2 ·sr −1 , 523.40, and 14.93, respectively. ISA and NTL are comprehensive indicators of urbanization, while BKD and RLD focus on a single aspect of urbanization. Therefore, the intragroup similarity between ISA and NTL, BKD and RDL is higher than the intergroup similarity in terms of spatial patterns. It is observed that water bodies, which are mainly composed of rivers, lakes, and reservoirs, play a crucial role in shaping the urban expansion of Wuhan according to the spatial distribution of urbanization indicators. In addition, BKD and RLD were both found to exhibit low values compared to the relatively values of high ISA and NTL in the southeastern part of the study area, which indicated local spatial disparities among urbanization indicators. The urban-rural gradient of NTL and BKD has a significant spatial variation as opposed to ISA and RLD. Although the spatial distribution of RLD is different from that of the original road, RLD also showed an obvious linear distribution characteristic.

Spatial Associations between Land Subsidence and Urban Development
The values of global bivariate Moran's I at four block scales are less than 0 (p < 0.01), indicating significantly negative spatial correlations between the four types of urbanization indicators and surface deformation ( Figure 5). That is, overall, urbanization is an important factor leading to land subsidence. However, the degree of negative spatial correlation varied with different urbanization factor types and scales. Among these types, the strongest negative correlation was found between land subsidence and BKD (Moran's I: ranging from −0.1911 to −0.1639), followed by that between land subsidence and NTL (Moran's I: ranging from −0.1793 to −0.1220) and that between land subsidence and RLD (Moran's I: ranging from −0.1224 to −0.0683). The weakest correlation was between land subsidence and ISA (Moran's I: ranging from −0.1121 to −0.0666). The negative spatial correlation between urbanization and land deformation manifested as a gradual increase in the block size.
The cluster maps of the bivariate local Moran's I in Figure 5 further illustrate the pattern of spatial heterogeneity in the relationships between land subsidence and urbanization at four block scales. We observed an obvious clustering pattern similarity between the spatial distributions of each of the urbanization indicators and land subsidence in the study area at four block scales. The low-high spatial outliers were mainly concentrated in the urban center of Wuhan, particularly in JA, JH, QK, and WC along the bank of the Yangtze River, and the southwest HS. A low-low cluster was significantly observed in the periphery of the study area, especially in the south HS, north JX. The high-low areas were

Spatial Associations between Land Subsidence and Urban Development
The values of global bivariate Moran's I at four block scales are less than 0 (p < 0.01), indicating significantly negative spatial correlations between the four types of urbanization indicators and surface deformation ( Figure 5). That is, overall, urbanization is an important factor leading to land subsidence. However, the degree of negative spatial correlation varied with different urbanization factor types and scales. Among these types, the strongest negative correlation was found between land subsidence and BKD (Moran's I: ranging from −0.1911 to −0.1639), followed by that between land subsidence and NTL (Moran's I: ranging from −0.1793 to −0.1220) and that between land subsidence and RLD (Moran's I: ranging from −0.1224 to −0.0683). The weakest correlation was between land subsidence and ISA (Moran's I: ranging from −0.1121 to −0.0666). The negative spatial correlation between urbanization and land deformation manifested as a gradual increase in the block size. clusters (except for low-high outliers) exhibited drastic declines in the number of clusters and the occupied areas and the low-high cluster appeared to have a concentrated distribution. The degree of global spatial correlation was the highest when low-high outliers became dominant agglomerates because clusters whose bivariate spatial autocorrelation type was opposing no longer canceled each other out. This suggests that the value of local bivariate Moran's I between land subsidence and urbanization indicators depends on the block scale to some degree.  The cluster maps of the bivariate local Moran's I in Figure 5 further illustrate the pattern of spatial heterogeneity in the relationships between land subsidence and urbanization at four block scales. We observed an obvious clustering pattern similarity between the spatial distributions of each of the urbanization indicators and land subsidence in the study area at four block scales. The low-high spatial outliers were mainly concentrated in the urban center of Wuhan, particularly in JA, JH, QK, and WC along the bank of the Yangtze River, and the southwest HS. A low-low cluster was significantly observed in the periphery of the study area, especially in the south HS, north JX. The high-low areas were mostly located across the region clustered in the northeast section of the study area. A high-high cluster was mainly distributed in the transition zone between East Lake and the Yangtze River and Hanyang (HY) district along the bank of the Han River, almost adjacent to the low-high regions. At block scale of 2 km × 2 km, the number of statistically significant clusters and their occupied areas both decreased compared to those at other scales in terms of the four types of urbanization indicators. In particular, the other three types of clusters (except for low-high outliers) exhibited drastic declines in the number of clusters and the occupied areas and the low-high cluster appeared to have a concentrated distribution. The degree of global spatial correlation was the highest when low-high outliers became dominant agglomerates because clusters whose bivariate spatial autocorrelation type was opposing no longer canceled each other out. This suggests that the value of local bivariate Moran's I between land subsidence and urbanization indicators depends on the block scale to some degree.

Impacts of Urbanization on Land Subsidence
To prevent the disturbance of potential multicollinearity on the parameter estimation of the model, each urbanization factor was independently analyzed with the land subsidence indicator in the GWR model due to high correlations existing among the urbanization indicators [50]. Thus, a total of 4 × 4 = 16 GWR models were generated and combined with four block scales. The standardized residual values of all GWR models range from −11.0298 to 6.7591 at most; more than 97% is in the range of −2.58-2.58. Therefore, the standardized residual values of all GWR models are randomly distributed at a 95% confidence level. To further examine whether the residuals from GWR models exhibit spatial randomness, a spatial autocorrelation analysis was performed on the residuals to obtain the global Moran's I statistics. As shown in Table 1, in general, low spatial autocorrelations of residuals from the GWR model are detected at the small block scale compared to the large scale, indicating that the variance in land subsidence over the study area is relatively random and exhibits spatial stationarity at the small block scale.  The local parameter estimates of the GWR model indicate the spatially varying relationships between the independent variable and response variable at different locations. The magnitude of the absolute value of the model regression coefficient denotes the degree of impact of an independent variable on land subsidence. In addition, the local adjusted determined coefficient (adjusted R 2 ) value from the GWR model is used to detect and assess the ability of the explanatory variable to explain the spatial variance in land subsidence, and a higher local adjusted R 2 value means better performance of the model. The ranges of local parameter estimate and adjusted R 2 between land subsidence and urbanization indicators obtained from GWR models are summarized in Table 2. Both positive and negative relationships are identified by local coefficients between urbanization indicators and land subsidence at different block scales, and the average of regression coefficients for NTL (except at the 1500 m scale) is less than 0. The local adjusted R 2 indicated that the urbanization indicators could explain more than 75% of the spatial variance in land subsidence on average at the four block scales (Table 2). In general, the explanatory power of urbanization indicators on land subsidence presents no significant difference from small to large block scales. However, the higher adjusted R 2 suggests that NTL has a stronger ability to explain the land subsidence changes than ISA at different block scales, which indicates that both ISA and NTL are relatively comprehensive indicators measuring the degree of urbanization, although NTL can reflect dynamic human activities and urban vitality better than ISA because artificial lights provide a direct signature of human activity [51]. As shown in Figure 6, the explanatory ability of GWR revealed by the local adjusted R 2 varies spatially. In general, GWR exhibits stronger explanatory power in WC, QS, central HS adjacent to East Lake, and northern JX. In contrast, the prediction ability of GWR appears to be lower in JA and around the periphery of the study area. The spatial patterns of coefficients of independent variables for urbanization indicators identified by GWR are clearly shown in Figure 7. In terms of the spatially varying regression coefficients, the directions (positive or negative) and strengths of the relationships between land deformation and urbanization indicators are not constant over the study area at different block scales. This result suggests that homogeneity and heterogeneity in the spatial relationships between land subsidence and urbanization indicators are sensitive to spatial scales. In addition, the spatial patterns of the regression coefficients of urbanization indicators tend to become more similar with increasing block scale. At the 2000 km block scale, negative relationships between urbanization indicators and land surface deformation are detected for JA, JH, QK, and WC along the Yangtze River and south HS, indicating that urban construction and anthropogenic activities resulted in the occurrence of land subsidence.
As shown in Figure 7, the negative values of the regression coefficients from the GWR model concentration region coincide with the distribution of soft soils or carbonate rocks, which have a high degree of urbanization. According to the definitions of urbanization indicators (ISA, NTL, BKD, and RLD), a region with a high degree of urbanization has large dynamic and static loads or is undergoing frequent construction and renewal activities. The urban construction of Wuhan city witnessed a stage of rapid development during the study period. In the construction process, groundwater extraction is required for operations around sites where buildings and subways are constructed. Previous studies have shown that the loss of groundwater results in consolidation of highly compressible soft soils and the dissolution of carbonate rocks, which thereby leads to land subsidence [30]. At the same time, the excavation of a subway tunnel inevitably disturbs the surrounding soil layers, followed by ground settlement. After the completion of project construction, i.e., in the process of building and subway operation, continuous dynamic and static loading act on the foundations of structures such as buildings, subways, and bridges. When the soil layer underneath a structure can no longer support the loading, settlement occurs within the structure and the surrounding area [32,41,52]. Remote Sens. 2022, 14, x FOR PEER REVIEW 14 of 19 Figure 6. Spatial patterns of local adjusted obtained from the GWR model for urbanization indicators at four block scales.
As shown in Figure 7, the negative values of the regression coefficients from the GWR model concentration region coincide with the distribution of soft soils or carbonate rocks, which have a high degree of urbanization. According to the definitions of urbanization indicators (ISA, NTL, BKD, and RLD), a region with a high degree of urbanization has large dynamic and static loads or is undergoing frequent construction and renewal activities. The urban construction of Wuhan city witnessed a stage of rapid development during the study period. In the construction process, groundwater extraction is required for operations around sites where buildings and subways are constructed. Previous studies have shown that the loss of groundwater results in consolidation of highly compressible soft soils and the dissolution of carbonate rocks, which thereby leads to land subsidence [30]. At the same time, the excavation of a subway tunnel inevitably disturbs the surrounding soil layers, followed by ground settlement. After the completion of project construction, i.e., in the process of building and subway operation, continuous dynamic and static loading act on the foundations of structures such as buildings, subways, and bridges. When the soil layer underneath a structure can no longer support the loading, settlement occurs within the structure and the surrounding area [32,41,52]. In contrast, some unexpected local relationships are also identified by the GWR model. For example, the southeastern and northeastern parts of the study area have concentrations of positive values. The density of deformation monitoring point pixels obtained in peri-urban areas is sparser than that in the center of the city due to dense vegetation and abundant waters. As a result, the land surface deformation monitoring data in the aggregated unit are easily influenced by random errors, and such unexpected results thus appear. Additionally, some omitted unknown variables in the GWR model may also contribute to the unexpected results.
1 Figure 7. Spatial distribution of regression coefficients between urbanization indicators and land subsidence.

Scale Effects of Relationships between Land Subsidence and Urbanization
The scale effect refers to smaller units being aggregated into larger units for spatial data, this likely provides inconsistent results [53]. Numerous studies have pointed out that the scale effect has been proven to be ubiquitous in the analysis of geographical phenomena and processes [54][55][56]. Neglecting the scale effect may lead to uncertainty in the results of spatial analyses and statistics and may even produce false conclusions. Thus, scale effects cannot be overlooked when analyzing the relationship between land subsidence and urbanization. In this work, the global and local bivariate Moran's I (Figure 5), the spatial patterns of the local adjusted R 2 (Figure 6), and the directions and strengths of the identified relationships between land subsidence and urbanization (Figure 7) all varied among the different block scales. In other words, the above results showed that the spatially varying relationship between land subsidence and urbanization was also scale dependent. First, the explanatory ability of the GWR models increased as the block scale increased, but approximately 20% of the spatial variance in land subsidence remained unexplained even at a block scale of 2 km. The reason for this may be that the larger block size was able to significantly reduce random disturbance for the estimation of land subsidence at the statistical unit due to the effect of spatial filtering. Nevertheless, the spatial autocorrelations of the residuals produced from the GWR models did not change accordingly with the variation in the block scale. Notably, significant positive spatial autocorrelations were found in the GWR model for land subsidence and ISA at the scales of 1.0 km (Moran's I = 0.1267, p < 0.01) and 2.0 km (Moran's I = 0.1150, p < 0.01), NTL (Moran's I = 0.1167, p < 0.01), and BKD (Moran's I = 0.1268, p < 0.01) at the 2.0 km scale. This finding indicated that although the results of the GWR model varied with changes in the block scale, certain bias occurred randomly in the calibration due to the interference of other factors that were not considered.
In fact, the interactions between land subsidence and urbanization indicators are very complicated, and thus it is difficult to design an appropriate method to sample spatial data. Therefore, to correctly analyze such interactions, a comparative analysis with multiscale data is very important and necessary. Although it is difficult to recommend the "best" spatial scales for the GWR model because the influences of factors vary over different sampling scales [57], the degrees of fitting for an observation measured by the adjusted R 2 and the spatial autocorrelation of the GWR model residuals should both be considered for the selection of a suitable block scale. In general, the suitable spatial scale in a GWR model should not only effectively prevent the interference of random factors but also reveal locally varying patterns of the identified relationships. Furthermore, the directions (positive or negative) and strengths of the spatial regression coefficients produced for the GWR models also varied across various scales and implied that spatial relationships between land subsidence and urbanization indicators gradually became global and that the spatial stationarity tended to be strong, which was consistent with the spatial resolution effects on the relationships between urban heat islands and their impact factors [58,59].
In this paper, an adaptive Gaussian kernel using the golden section search method was adopted to identify the optimal bandwidth size for a GWR model and mainly focused on the effects of block scale [60]. Future work is needed to consider a comparative analysis with multi-bandwidth models for GWR. In addition, we only analyzed spatially varying relationships between land subsidence and urbanization indicators at four block scales. As a result, it remains unknown whether there is a threshold within which scale effects are significant. It is therefore necessary to conduct a further comparative study of the GWR results obtained for a series of block scales [61].

Conclusions
In this study, the land deformation in Wuhan, the largest city in central China, obtained from Sentinel-1 SAR time-series datasets based on the SBAS-InSAR method was presented. We investigated the spatially varying relationships between land subsidence and urbanization in Wuhan by using bivariate Moran's I and GWR models. Our analysis concentrated on a set of empirical results that support the following conclusions. The derived deformation results showed that land subsidence was obvious and unevenly distributed in the study area, the annual deformation rate varies from −42.85 mm/year to +29.98 mm/year, and its average value was −1.0 mm/year. A clear spatial pattern for land subsidence in Wuhan was mapped, and several apparent subsidence funnels identified by the LISA index were primarily located in central urban areas. As indicated by bivariate global and local Moran's I, four types of urbanization indicators ISA, NTL, BKD, and RLD, were found to be significantly spatially correlated with land subsidence at different scales, which implied that the urbanization indicators could have an impact on the land subsidence of its surrounding neighbors. In addition, the GWR model results showed that all urbanization indicators were significantly associated with land subsidence across the whole study area in Wuhan, mainly controlled by the thickness of soft soil, but the relationships were not completely consistent among land subsidence and different urbanization indicators. The results of bivariate Moran's I and GWR confirmed that the relationships between land subsidence and urbanization varied spatially in Wuhan at multiple spatial scales. Moreover, scale dependence existed both in bivariate Moran's I and GWR models for land subsidence and all urbanization indicators, however, a "best" spatial scale could not be confirmed because the disturbances of factors vary over different sampling scales. We suggest that the results from our study advance the understanding of the spatially varying relationships between land subsidence and urbanization, and it is hoped that they will provide guidance for subsidence control and sustainable urban planning.