Karst Collapse Risk Zonation and Evaluation in Wuhan, China Based on Analytic Hierarchy Process, Logistic Regression, and InSAR Angular Distortion Approaches

: The current study presents a detailed assessment of risk zones related to karst collapse in Wuhan by analytical hierarchy process (AHP) and logistic regression (LR) models. The results showed that the LR model was more accurate with an area under the receiver operating characteristic (ROC) curve of 0.911 compared to 0.812 derived from the AHP model. Both models performed well in identifying high-risk zones with only a 3% discrepancy in area. However, for the medium-and low-risk classes, although the spatial distribution of risk zoning results were similar between two approaches, the spatial extent of the risk areas varied between ﬁnal models. The reliability of both methods were reduced signiﬁcantly by excluding the InSAR-based ground subsidence map from the analysis, with the karst collapse presence falling into the high-risk zone being reduced by approximately 14%, and karst collapse absence falling into the karst area being increased by approximately 6.5% on the training samples. To evaluate the practicality of using only results from ground subsidence maps for the risk zonation, the results of AHP and LR are compared with a weighted angular distortion (WAD) method for karst risk zoning in Wuhan. We ﬁnd that the areas with relatively large subsidence horizontal gradient values within the karst belts are generally spatially consistent with high-risk class areas identiﬁed by the AHP- and LR-based approaches. However, the WAD-based approach cannot be used alone as an ideal karst collapse risk assessment model as it does not include geological and natural factors into the risk zonation.


Introduction
Karst collapse is a significant geological hazard that occurs when subsurface physical and/or chemical dissolution of carbonate rocks form voids and cavities at depth that can propagate towards the surface and result in a collapse structure [1,2]. Environmental problems associated with karst collapse have a serious impact in many regions of the world, such as the Dead Sea area [3]; Florida and Texas, USA [4]; Hamburg, Germany [5]; Ebro Valley, Spain [6]; Elba Island, Italy [7]; Guangzhou, China [8]; Hamadan, Iran [9]. Similarly, in China, soluble rocks cover an area of 346.3 × 10 4 km 2 in the country, accounting for more than one-third of the total land area. The karst collapse incidents arising from this geological environment have affected more than 30 metropolises and 420 counties (districts) nationwide [10]. Among all cities, Wuhan, built on a well-developed karst ecosystem, is the most vulnerable mega-city (with people over 12.6 million) to karst collapse. From 1977 to 2018, 38 regions in Wuhan were affected by karst collapse, with the number of sinkholes reaching 102. These karst collapse incidents not only directly caused at least CNY 114.091 million in economic losses but also seriously jeopardized the safety of residential areas, transportation, engineering construction, and water supply [11]. Management of such unpredictable, but everlasting, hazards requires a detailed study of karst collapseinducing mechanisms and risk assessment of areas vulnerable to collapse.
In terms of the research on karst collapse-inducing mechanisms, previous studies have identified the main factors acting as triggers for karst collapse in Wuhan [10,12,13]. These include well-developed karst settings, the typical dualistic structure of the overlying soil, and the interactive runoff between the Yangtze River water, pore water, and karst water, as well as the extrinsic triggers such as engineering dewatering from large-scale municipal construction [14,15]. Moreover, localized ground subsidence was also observed at historical karst collapse sites in Wuhan, similar to that of Elba Island, Italy [7] and Alagoas, Brazil [16]. For the vulnerability assessment of karst collapse in Wuhan, all the existing studies adopted the analytical hierarchy process (AHP) method [17,18], of which the hierarchy modeling only included the essential controlling conditions. Neither the external predisposing factors nor the contemporary ground subsidence map of Wuhan was considered in the modeling [11,19,20].
In this paper, we extend previous susceptibility studies in Wuhan and present the outcomes of karst collapse risk zonation using two approaches commonly used in geohazard assessment, namely the qualitative AHP and the semi-quantitative logistic regression (LR) [21][22][23]. For the AHP and LR calculation, in addition to typical controlling factors included in the existing literature in Wuhan, we also consider the influence from urban land planning, the dynamic/static load factors from municipal traffic/high-rise buildings, and the ground deformation rates derived from the multitemporal InSAR processing of Sentinel-1A (April 2015 to June 2019). Moreover, we present results from Weighted Angular Distortion (WAD), derived from integration of InSAR-based ground subsidence rates, subsidence horizontal gradient (SHG), and municipal construction density for karst risk zoning in Wuhan. The WAD approach has been recently adopted as a simple, but practical and reliable, approach to evaluate surface faulting and subsidence risks in areas such as Beijing, China [24], Mexico City [25], and Alagoas, Brazil [16]. We compare the results from these different techniques in a common setting and evaluate the accuracy of AHP-and LR-based approaches, which provides important inferences on how to accurately assess karst collapse risk modeling in Wuhan.

Study Area
The study area covers seven central urban regions (I~VII in Figure 1a) of Wuhan city, with a total area of approximately 3260 km 2 . The plains, with a few low mountains and hills, form the basic geomorphological unit of the study area. Geologically, six karst belts (L1-L6 in Figure 1b) with individual lengths of about 35-63 km and widths of 0.5-15 km stretch throughout the city along a WNW-ESE direction. They cover a total area of 960.95 km 2 , accounting for 29.5% of the whole study area. Such karst geological conditions coupled with the intensive urban construction over the past 20 years, including nearly 12 subway lines and over 10,000 construction sites, resulted in 84 sinkhole incidents in the study area from 1977 to 2018 (Figure 1b). Moreover, 87% of these sinkholes occurred in the Baishazhou-Jiangdi zone, where the typical dualistic structure of the upper soft soils and lower sand is distributed (Figure 1b).

Materials and Data Layers Grading
For this study, we collected nine types of datasets in terms of geological environment intrinsic conditions and extrinsic trigger conditions as listed in Table 1. Geological and environmental intrinsic conditions include the following criterion layers: karst geology conditions (B1), overburden conditions (B2), hydrogeological conditions (B3), and karst surface subsidence conditions (B4). The extrinsic trigger condition includes the criterion layer of anthropological activities (B5). Below, we detail different components of the criterion layer.
Karst geology conditions (B1): these consist of stratigraphic lithology (C11) and the development degree of karst (C12). The soluble rock strata in the study area are mainly composed of dolomitic limestone of the Triassic Lower Middle System, flint tuberculous limestone of the Permian Lower System, thick-bedded spheroidal tuffs and dolomites of the Carboniferous Upper System, and clay-hosted tuffs of the Lower Carboniferous. They form two types of karst strata in Wuhan, i.e., covered and buried karst, which, together with non-carbonate areas, constituted the high, medium, and low susceptibility zones of karst collapse in Wuhan (Figure 2a). In addition, the development degree of karst determines the space for overlying soil storage and transport [1], which is an important predisposing factor for karst collapse. Therefore, we divided the development degree of karst in the study area into well-developed (κ > 10%), medium-developed (10% ≥ κ ≥ 3%), slight-developed (κ < 3%), and non-soluble areas (Figure 2b) according to the China karst collapse survey specification (1:50,000) [12].    Overburden conditions (B2): these include the overlying soil structure (C21) and thickness (C22). Most karst collapses occur in the dualistic structure overburden area of the upper soft soils and lower sand in the Baishazhou-Jiangdi zone (Figure 1b), with the soil thickness generally <15 m. This was closely followed by the multi-layered soft soil structure area (clay and sandy soil) with a thickness of 15-30 m, and the buried karst and single-layered soil structure area with a cover thickness of 30-40 m. A few karst collapses were recorded in areas with the soil thickness exceeding 40 m. Therefore, the overlying soil structure and thickness were classified into three and four categories, respectively, based on the recorded karst collapse incidents, as shown in Figure 2c,d.
Hydrogeological conditions (B3): Ref [14] revealed that there existed a frequent interactive runoff recharge of "Yangtze water-pore water-karst water" in the karst zones of the younger terrace on both sides of the Yangtze River, which promotes the formation and expansion of subsurface cavities. Therefore, we spatially analyzed the recorded karst collapses and ranked the proximity to the 4th class rivers (C31) as 0-1000 m, 1000-3000 m, 3000-5000 m, and >5000 m ( Figure 2e). In addition, the Quaternary pore water abundance (C32) ranked as >1000 m 3 /d, 100-1000 m 3 /d, and <100 m 3 /d according to the water yield data (i.e., cubic meters per day, m 3 /d) from 361 karst boreholes in the study area ( Figure 2f).
Karst surface subsidence conditions (B4): Ground subsidence can usually be regarded as an indicator of the localized sinkhole formation process, and this is the case in Wuhan. Spatial analyses of ground subsidence and sinkhole locations in karst zones showed that areas with dense historical sinkholes were often accompanied by severe localized ground subsidence [14,15]. To consider the influence of contemporary ground deformation in risk modeling, we exploited C-band Sentinel-1 SAR data, covering the period from April 2015 to June 2019, and analyzed them using the StaMPS-SBAS method [26] to derive the ground subsidence rates in Wuhan. Subsequently the InSAR-based ground subsidence rates (C42) in the study area [14,15] were divided into four categories using the natural break algorithm, i.e., −89.7-−5.8 mm/yr, −5.7-−1.3 mm/yr, −1.2-2.3 mm/yr, and 2.4-29 mm/yr, corresponding to the high-, medium-, low-, and non-prone risk areas, respectively ( Figure 2g).
Anthropological activities (B5): In addition to the four geological environment intrinsic conditions mentioned above, Ref [14] suggested that engineering dewatering from underconstruction subway lines and construction sites, and water depletion from industrial, commercial, and urban residential production and living, promotes the development of underground cavities and karst surface subsidence. Therefore, we conducted a distancebased multi-ring buffer analysis for all subway lines and construction sites with building heights over 250 m (C51) in the study area and classified them into four categories, i.e., <500 m, 500-1000 m, 1000-2000 m, and >2000 m ( Figure 2h). In addition, the urban land statutory planning provides a glimpse of the future engineering intensity in the karst area. Therefore, we sorted the urban planning map (C52) into four categories, i.e., manufacturing (M); transportation (T); municipal utilities (U); warehouse (W); land area, residential (R); commercial (C); land area, green space and agriculture (G); land area, and ecologically controlled land and water area (E) (Figure 2i).
The hierarchy, influence grading, and corresponding scale values of the above four geological environment intrinsic conditions (B1-B4) and anthropological activities (B5) with a total of nine sub-criteria are shown in Table 1.

. Decision Matrices and Consistency Test
We used the T.L. Satty 1-9 scaling method [27] for the pair-wise comparison of all factor layers under each criterion layer by relatively ranking their prevalence to karst sinkhole formation. The eigenvector associated with the principal eigenvalue of the pair-wise comparison matrix was used as the weight. As the comparison matrix is not necessarily a consistency matrix, we performed the consistency index (CI) and consistency random (RI) tests. Further, the consistency ratio (CR) could be defined as [27,28]: where λ max is the principal eigenvalue of the pairwise comparison matrix and n is the order of the matrix. The comparison matrix is entirely consistent if CI is equal to zero; a smaller value of CI corresponds to a better consistency. On the other hand, if CR is equal to 0.1, the comparison matrix has a good consistency, and the grading is reasonable. Otherwise, it does not meet the consistency principle and needs to be adjusted until the test is passed. The relative weights of the sub-criteria and data layers are provided in Table 2. Note where W C ij and W B i are the normalized weights of factor and criterion layer, respectively; W C ij −B i is the combined weights calculated by

Karst Collapse Susceptibility and Risk Assessment
We adopted the same method as above to assign the weight of karst collapse susceptibility (λ susc ) and the two triggering sub-criteria (λ C 51 , λ C 52 ). In the end, based on the comprehensive index evaluation model of hazards (Equation (2)), the above three types of layers (λ susc , λ C 51 , λ C 52 ) were rasterized and superimposed on each other for karst collapse risk zonation in Wuhan.

The LR-Based Approach to Risk Assessment
The logistic regression model is a common statistical analysis method for dichotomous problems, which explores the relationship between a binary dependent variable (usually 0 for the absence of geohazards and 1 for the presence of geohazards) and a set of independent variables [29,30]. The logistic regression function is given as: where p is the probability of the event (here sinkhole occurrence); b 0 is the intercept; b 1 , · · · b n are the coefficients of independent variables (X 1 , · · · X n ).
Taking the natural logarithm of both sides of the Equation (3) yields ln[p/(1 − p)], and using it as the dependent variable and all sub-criteria as the independent variables, the logistic regression equation Equation (4) is established: Since the magnitudes and scales of the sub-criteria vary in ranges, the logistic regression analysis requires unifying the sub-criteria under the same scale. To this end, based on the historical distribution of 84 sinkholes in the study area, we adopted the certain factors model [31] to calculate the relative weights (CF) of each grading under the factor layer.
where PP a is the conditional probability of a sinkhole occurring on the i−th level at a specific sub-criterion, denoted as the ratio of the number of sinkholes fallen in the i−th level to the area of this level; PP s is the ratio of the total number of karst sinkholes to the area of the study area.
The derived CF values CF ij {i = 1, 2, 3, 4, 5; j = 1, 2, . . . , 4} were used as the independent variables (X 1 , · · · X n ) for the LR model. The recorded sinkhole presence and the randomly generated sinkhole absence data were divided into training and test samples according to 80% and 20%, respectively. The CF values of the training sample points were input to Equation (4) to construct a binary LR model and derive the regression coefficients b 1 , · · · b n . The probability of occurrence of karst collapse in each unit was subsequently derived by substituting the b 1 , · · · b n back into Equation (3). At last, following the China karst collapse survey specification (1:50,000), the natural break algorithm in ArcGIS ® was adopted to perform the karst collapse risk zoning in Wuhan, and the rationality and accuracy of the constructed LR model were analyzed based on the reserved 20% sample.

The Weighted Angular Distortion Approach to Risk Assessment
After obtaining the mean ground subsidence rates using the StaMPS-SBAS technique, the risk zoning map for the cumulative period (April 2015 to June 2019) was produced based on the angular distortion β, which was calculated as the ratio of the subsidence horizontal gradient (SHG) between two adjacent pixels to the horizontal distance between them [32][33][34] using a slope algorithm available in ArcGIS ® . Compared to vertical ground deformation (V sub ), a large subsidence horizontal gradient (SHG) around deformation clusters contributes more significantly to risk zoning results, and the local municipal construction density (MC den ) represents, somehow, an indicator of vulnerability [33]. Therefore, we assigned three times more weight to the subsidence horizontal gradient than ground subsidence rates and introduced the municipal construction density calculated as the area ratio of subway lines, large construction sites, and the unit (90 m × 90 m) based on the results in Figure 2h. The proposed weighted angular distortion model (WAD) was defined as: Finally, we divided the risk map into four categories by using the natural break algorithm in ArcGIS ® .

Risk Zonation by AHP-Based Approach
According to the procedures described in Section 3.1, we first conducted the pairwise comparison for each sub-criterion under each criterion layer and obtained the judgment matrix ( Table 2). As seen from Table 2, the CI value of the four criterion layers is 0.0014, the CR ratio is 0.0016, and the consistency of the factor layer under each criterion layer is zero, which satisfies the consistency test. On this basis, normalization of the principal eigenvector associated with the judgement matrix for predisposing factors yielded the largest weight to overburden conditions (B2) (0.4486) followed by karst geology conditions (B1) (0.2347), karst surface subsidence conditions (B4) (0.2347), and hydrogeological conditions (B3) (0.082). The results are shown in the last two columns in Table 2.
Subsequently, based on the quantitative grading map of each sub-criterion in Figure 2a-i, the AHP-based karst collapse susceptibility evaluation model was achieved by using the weighted comprehensive index method: After obtaining the karst collapse susceptibility (Sus AHP ), we superimposed the influence of anthropological activities and used the same method in Section 3.1 to calculate the sub-criterion weights, as shown in Table 3.
Based on the natural break algorithm in ArcGIS ® , the areas with Risk AHP values between 2.4392 and 4.9089 were defined as high-risk, between 1.6566 and 2.4392 as mediumrisk, and between 0.8435 and 1.6566 as low-risk ( Figure 3).

Risk Zonation by LR-Based Approach
For the LR-based model, the created fishnet tool first generated a 90 m × 90 m grid in each sub-criterion layer. Sixty-eight sinkhole locations (80% of the total sinkholes) were randomly selected as training samples (we refer to them as sinkhole presence, and the opposite are sinkhole absence). In addition, in order to not lose generality, a certain percentage of sinkhole absence data were still needed to be involved in constructing the LR model. LR modeling using an unbalanced dataset (e.g., presence to absence ratio = 1:1.78) would result in a decreased root mean square error from 0.44 to 0.42 and an increase in model R 2 from 0.25 to 0.32 when compared to a balanced dataset [22,35,36]. Therefore, we generated sinkhole absence data by randomly locating 122 data points (1.78 × 68) within the karst zones according to the ratio in [22]. A total of 190 sinkhole presence and absence training data were used to construct the logistic regression model. Table 4 shows the CF values for each grading under each factor layer. We extracted the CF values of nine sub-criteria at 190 training samples and imported them into SPSS ® software [37] for binary logistic regression analysis. The correlation analyses of each sub-criterion using SPSS ® revealed that the correlations of overlying soil thickness with overlying soil structure and proximity to the fourth class rivers were −0.534 and −0.507, respectively. The correlations of proximity to subway lines and construction sites with the development degree of karst and urban planning map were −0.641 and 0.758, respectively. The correlations between other sub-criteria were less than 0.5. To avoid the accuracy loss of the LR model due to the involved redundant factors, overlying soil thickness and proximity to subway lines and construction sites were not included in the final LR model. Finally, we obtained the regression coefficients B for these seven factors and the intercepts, which were 5.640, 1.685, 2.343, −0.642, 2.606, 3.879, 0.104, and −1.961. The correlation matrix and LR parameter estimates for the determined evaluators are shown in Table 5 and Figure 4.  In the LR-based geohazard susceptibility assessment, B represents the weight of each sub-criterion while the sig reflects their significance in the model; the smaller the sig value, the more significant the sub-criterion is in the LR model. Moreover, the model is statistically significant only when the sig value is less than 0.05 [29]. In view of this, as seen from Figure 4, except for the proximity to the fourth class rivers and urban planning map, the sig values of other sub-criteria are less than 0.05. Eventually, this effective LR-based karst collapse risk assessment model ranked the seven predisposing factors to karst collapse risk (odds ratio, EXP(B)) as stratigraphic lithology (281.440), InSAR-based ground subsidence rates (48.383), Quaternary pore water abundance (13.539), overlying soil structure (10.414), development degree of karst (5.390), urban planning map (1.110), and proximity to the fourth class rivers (0.526) (Figure 4).
The regression coefficients of each sub-criterion were then substituted into Equation (4) to yield the logistic regression equation: where C11, . . . C52 are CF values of seven determined sub-criteria.
Referring to the China karst collapse survey specification (1:50,000) [12], we used the natural break algorithm in ArcGIS ® to categorize the areas with probability values between 0.2591 and 0.9999 as high-risk, between 0.0318 and 0.2591 as medium-risk, and between 0 and 0.0318 as low-risk. The other areas were less likely to experience karst collapse, and, thus, we obtained the karst collapse risk zoning map in Wuhan ( Figure 5).

Risk Zonation by InSAR-Based Angular Distortion Approach
A simple classification of the derived SHG from InSAR-based ground subsidence rates from April 2015 to June 2019 was conducted by accounting the geotechnical practice on practical limits for the allowable settlement of buildings (Figure 6a). For angular distortion level we took the lead from previous studies and considered four categories, low (β ≤ 1/3000), medium (1/3000 < β ≤ 1/1500), high (1/1500 < β ≤ 1/500), and very high (β > 1/500), indicating an increased likelihood of damage [32,34]. However, we adopted the risk level-I-II-III-IV (1/10 of the four categories) instead of the above four categories (very high, high, medium, low) because the ground subsidence rates in this study area during this period were generally low, and no grids with the derived SHG exceeded 0.002 rad, i.e., 1/500. As with the previous studies [14,15], Hankou and Wuchang (HK-WC) suffered the most severe ground subsidence with a maximum subsidence rate of −89 mm/year (Figure 2g). However, the geological conditions of non-karst strata and overlying multi-layered soft soils led to ground subsidence in Hankou and Wuchang, characterized mainly by consolidation mechanisms and with a large scale of uniform subsiding. Such an extent and distribution of ground subsidence resulted in moderate deformation gradients with a maximum SHG of 1/826 and 99.5% of the grids having SHG values less than 1/3000. In contrast, among the six karst belts, especially in the Tianxingzhou (L1) and Baishazhou karst belts (L3), localized ground subsidence led to relatively large deformation gradients at the edges of the deformed and non-deformed areas (a total of 339 grids with a cumulative area of 3.1 km 2 falling into the medium-high risk zone), which, in turn, caused severe tear damage to buildings (Figure 6a).
In addition, ground subsidence was observed along subway lines and around large construction sites (with building heights over 250 m) in karst zones. Therefore, we calculated the point and line densities MC den (Figure 6b) based on the results in Figure 2h and integrated them into Equation (6) to derive the calibrated risk values. In the end, the natural break algorithm was used to divide the calibrated risk values into four risk classes, corresponding to high-risk zone (−0.2271~−0.0386), medium-risk zone (−0.0386~−0.0099), low-risk zone (−0.0099~0.0014), and non-risk zone (0.0014~0.0120) (Figure 6c).

Comparison of Karst Collapse Risk Zonation Results
According to the GIS statistics of the three risk zonation results, approximately 17.7% (170.52 km 2 ) of the total karst area (960.95 km 2 ) in the risk map derived from the AHP-based approach fell into the high-risk class. The medium-and low-risk classes covered almost 51.6% (495.46 km 2 ) and 30.7% (294.97 km 2 ) of the total karst area. In contrast, the high-, medium-, and low-risk classes derived from the LR-based approach accounted for 14.4% (139.06 km 2 ), 36.0% (345.84 km 2 ), and 49.6% (476.05 km 2 ) of the total karst area, respectively.
The karst collapse high-risk areas identified from AHP-and LR-based approaches were located in (1) Jiefang Avenue, Houhu Avenue, Heping Avenue, and GongrenVillage in the L1 karst belt (zone 1 in Figures 3 and 5); (2) Moshui lake north road, Hanyang Avenue, and the area west of the Xunsi River in the L2 karst belt (zone 2 in Figures 3 and 5); (3) Hanyang Yingwu-Jiangang Road area, Wuchang Jiefang Bridge-Wutai Gate, Fenghuo Village, Maotan Port, and Baisha Road area in the L3 karst belt (zone 3 in Figures 3 and 5). The spatial distribution of high-risk areas identified by both approaches were highly correlated, with only a 3.3% area discrepancy. This suggested that the subjective judgment approach based on a priori information on karst collapse and the semi-quantitative approach based on the recorded sinkhole presence were comparable in terms of high-risk class zonation.
In terms of medium-and low-risk karst collapse classes, the zoning results of the AHP-and LR-based approaches were similar in spatial distribution; however, there were distinct discrepancies in the category areas. The AHP-based approach tended to allocate more areas to the medium-risk class (15.6% more) and less areas to the low-risk class (2.7% less). There are several reasons for this discrepancy in the results between AHP and LR. One could be the non-uniformity of rules in assigning weights (or scales) to each sub-criterion between the subjective judgment-based AHP model and the certain factorsbased semi-quantitative LR model. For instance, the AHP-based approach ranked the first four predominant predisposing factors to karst collapse as overlying soil structure (0.2991), proximity to subway lines and construction sites (0.2777), InSAR-based ground subsidence rates (0.2347), and development degree of karst (0.1565). However, for the LR-based approach, the stratigraphic lithology (281.440) had a dominant role followed by InSAR-based ground subsidence rates (48.383), Quaternary pore water abundance (13.539), and overlying soil structure (10.414). Alternatively, no rigid risk-grading criteria were available in karst collapse engineering practice. Instead, a natural break algorithm was used to subdivide the two types of derived risk values, using scalar values in AHP and using probability in LR, which results in different spatial extents of the risk areas in two methods.
In addition, the WAD-based approach, as an approach to risk zoning based directly on measured data, is essentially distinguished from the above two qualitative and semiquantitative models that rely on historical or empirical information. Although it was impossible to perform a strict grid-by-grid comparison, the areas with relatively large SHG values within the L1, L2, and L3 karst belts were generally spatially consistent with the high-risk class areas identified by the AHP-and LR-based approaches. For example, there were 217 grids with SHG values between 1/1500-1/500 (high susceptibility area) and a total area of 1.68 km 2 along the Houhu Avenue and Heping Avenue in the L1 karst belt, Moshui lake north road in the L2 karst belt, Hanyang Yingwu -Jiangang Road, and in the Fenghuo Village zone in the L3 karst belt (zones 1 , 2 , and 3 in Figures 3 and 5).
In order to evaluate the importance of including information from ground deformation into the risk models, we made a test by excluding InSAR-based ground subsidence rates (C41) from both the AHP-and LR-based approaches. As seen in Figure 7a,b the zoning results of the AHP and LR without the sub-criterion C41 vary significantly in the area of each risk zone compared to the model that considers ground deformation in Figure 5a,b. For instance, the area of high-, medium-, and low-risk classes varied by 6.9%, −4.7%, and −2.2% for the AHP-based approach, and −37.1%, 33.2%, and 3.9% for LR-based approach, respectively. Moreover, the new models reduced the karst collapse presence, falling into the high-risk zone by 15% (AHP) and 13%, (LR) respectively (Figure 7c), and increased the karst collapse absence, falling into the karst area by 5% (AHP) and 8% (LR) on the training samples (Figure 7d). This proves the necessity of including ground subsidence rates into karst collapse risk zoning. Note that the karst collapse risk situation reflected by the WAD-based approach was milder than the AHP-and LR-based results. This is because the WAD-based approach focuses on the tearing damage to buildings by inhomogeneous ground deformation. Although Wuhan showed different patterns of ground subsidence, the extents were generally small and there were no grids in the whole study area that exceeded the threshold of building tear resistance (with an SHG value of 0.002 radian) [32,34] to be included in very high-risk class, i.e., β > 1/500 ( Figure 6a). Therefore, the risk levels-I-II-III-IV in Figure 6a represent only the relative likelihood of collapse.

Zoning Results Test and Model Evaluation
We tested the reasonability of AHP-and LR-based approaches for karst collapse risk zoning and their accuracy based on the test samples (20% sinkhole presence and its 1.78 times sinkhole absence, i.e., 16 and 29, respectively). We conducted this test by (1) counting the proportion of sinkhole presence within each risk class to test model reasonability, and (2) using the receiver operating characteristic curve (ROC) to evaluate model accuracy [38][39][40]. To test the model reasonability, we calculated the ratio of sinkhole presences falling in each risk class to all sinkhole presences in the test samples (Perc i pres ), the percentage of the area of each risk class to the total area of karst zones (Perc i area ), and their ratio (Ratio i ). A reasonable karst collapse risk assessment model implies that, for the test samples, the percentage of sinkhole presences falling into the high-risk class should be the largest, and the percentage of the low-risk area over the whole study area should be the largest; further, Ratio non < Ratio low < Ratio medium < Ratio high .
From Table 6, the AHP-based approach predicted that 81.3% (13/16) of the sinkholes occurred in the high-risk zone, accounting for 17.7% of the total karst area. The LR-based approach performed better in mapping sinkhole risks, as it predicted 93.8% (15/16) of sink-holes within the high-risk zone, which accounted for 14.4% of the total karst area. Moreover, the zoning results derived from the LR-based approach satisfied the above three hypotheses, i.e., Perc high pres , Perc low area accounting for the maximum, and Ratio low < Ratio medium < Ratio high . However, the AHP-based approach didn't fulfill the hypothesis that Perc low area accounts for the maximum. Table 6. Verification results of AHP-and LR-based approaches for the karst collapse occurrence in Wuhan using the test samples. The ROC curve and the area under the ROC curve (AUC), plotted in a frame with the true positive rate (sensitivity) as the longitudinal axis and the false positive rate (1specificity) as the horizontal axis, and provides a more intuitive measure of the accuracy of the two models [38]. Therefore, we calculated the ROC curves for all the models using the test samples, and compared their AUC values to verify the accuracy of the AHP-and LR-based approaches. Note that, for the AHP model, to calculate its ROC, we normalized the derived risk values to approximate them as the probability of karst collapse (0 to 1). However, for the LR model, we directly calculated the karst collapse probability based on Equation (4), constructed from the training samples and the corresponding seven types of CF values at the 45 test points. Finally, the test samples' sinkhole presence and absence (1 and 0) were used to perform the ROC calculation for both models (Figure 8). From the results of the reasonability (Table 6) and accuracy tests (Figure 8), both modeling approaches performed well in karst collapse risk zoning. The AUC values derived from the AHP and LR models were 0.812 and 0.911, respectively, which both exceeded 0.8, indicating goodness of fit for karst collapse risk zoning for this study. However, the LR-based approach was relatively better than the AHP-based approach in terms of model accuracy, with the confusion matrix of the training data showing an overall accuracy of 98.4% for the model. This is due to the fact that (1) the Perc high pres and Perc low area of the LR-based approach were 12.5% and 18.9%, larger than those of the AHP-based approach, respectively, and (2) the AUC for the AHP-based approach was almost 10.9% smaller than those for the LR-based approach.

Analysis of the Applicability of Three Approaches on Karst Risk Zonation
Qualitative and semi-qualitative approaches to karst collapse risk zoning depend on correctly identifying predisposing factors and their relative contributions to sinkhole formation. The extensive work conducted by previous researchers on the recognition and contribution determination of predisposing factors to karst collapse and subsidence in Wuhan [13][14][15] provided valuable a priori information for selecting sub-criterion in the risk zoning modeling in this study, further ensuring the accuracy of karst collapse risk zoning by qualitative and semi-quantitative models (e.g., AHP and LR). However, for some karst collapse areas where a priori information is lacking (e.g., Guizhou, China), the biggest challenge of the AHP-or LR-based approach is the uncertainty in prioritizing predisposing factors to sinkhole formation. In particular, the AHP-based approach requires multiple evaluations through a trial-and-error process to generate logically consistent relative weights to predict existing sinkhole incidents better. It is worth noting that this trial and error approach is the introduction of bias towards sinkhole presence. Consequently, this approach may allocate more areas to the higher karst collapse risk classes (17.7% in this study), thus limiting its applicability in sinkhole risk mitigation strategies. For the LR-based approach, although the certain factors method could aid in the determination of the relative importance of sub-criteria on karst collapse occurrence, it remains a challenge to specify which factors are selected for modeling. The WAD-based approach proposed herein can reflect the vulnerability of surface structures the most objectively and quantitatively, particularly for localized collapse or deformation such as sinkholes. Therefore, it can aid the AHP-based approach in assigning scales and the grading of each sub-criterion in both models. However, the WAD-based approach cannot be used alone as an ideal karst collapse risk assessment model, as it does not consider geological and natural factors in the risk analysis, and is more suitable for cases where a single factor dominates ground subsidence. Moreover, if a localized subsidence pattern similar to that of a sinkhole exists in the non-karst area, it would lead to logical errors to stubbornly use this approach for karst collapse risk zoning.

Conclusions
The study presents a comparison between a subjective judgment-based AHP model, a certain factors-based semi-quantitative LR model, and the quantitative WAD-based approach to risk zonation related to karst collapse in Wuhan city. Our results suggest that the LR-based approach is, in general, superior to the AHP-based approach, as it (1) better satisfied the assumption that "the percentage of sinkholes falling into the high-risk class and the percentage of low-risk area are the largest", and (2) its area under the receiver operating characteristic (ROC) was 0.911 compared to the 0.812 obtained with the AHPbased approach. Nevertheless, both models performed well in identifying high-risk areas with a highly correlated spatial distribution and only a 3.3% discrepancy in area. For the medium-and low-risk classes, their spatial distribution of risk zoning results were similar; however, the category areas were varied. In particular, we found that inclusion of the InSAR-based ground subsidence velocities significantly improves the correct prediction of sinkhole presence and the ability to avoid incorrect prediction of sinkhole absence for both models. However, given the non-uniformity rules of the two models in assigning weights (or scales) to each sub-criterion, the LR-based approach ranked the stratigraphic lithology (281.440) as a dominant factor, followed by InSAR-based ground subsidence rates (48.383), Quaternary pore water abundance (13.539), and overlying soil structure (10.414), while the AHP-based approach identified overlying soil structure (0.2991) as the main factor for karst collapse, followed by proximity to subway lines and construction sites (0.2777), InSAR-based ground subsidence rates (0.2347), and the development degree of karst (0.1565).
In addition, the WAD-based method proposed in this study suggested that the regions with relatively large SHG values within the L1, L2, and L3 karst belts were generally spatially consistent with the high-risk class areas identified by the AHP-and LR-based approaches. Although it focused more on quantitatively deriving the tearing damage caused by local collapse or deformation (e.g., sinkholes), the WAD-based method could aid the AHP-based approach in assigning scales and the grading of each sub-criterion in both models, especially for some karst collapse areas where a priori information is lacking.