Abstract
Due to the lack of historical subsidence data on goaf sites alongside expressways, it is difficult to comprehensively evaluate the stability of goaf sites and accurately identify potential risk points. An investigation to study the stability of goaf sites in the Huaibei section of the Xu Huaifu Expressway is conducted through SBAS-InSAR technology and field investigation. Initially, the stability of the goaf site around the route is anticipated through field investigation. Afterwards, based on the extraction of Sentinel-1 satellite images, subsidence data of the goaf sites from 2017 to 2021 were extracted. Finally, the stability of risk points is analysed. Results demonstrate that although mining in the area stopped more than five years ago, the goaf site has continued to slowly sink. According to four years of goaf site subsidence data extracted by SBAS-InSAR technology, two subsidence risk areas, namely, Y1-05 and Y1-09, were discovered, and the subsidence areas displayed an obvious inverted “convex” shape. Of these two areas, Y1-05 is a potential risk point that was not found during the field investigation. The subsidence basin in the Y1-05 area has developed gradually, and the residual subsidence is expected to be close to 475 millimetres within the next 15 years. Although it can be considered a stable area, the subsidence value is quite large and may pose risks to route safety. Therefore, a flexible reinforced mesh can be adopted for the route. The subsidence basin of Y1-09 is the largest, but the expected residual subsidence is only 50 millimetres and will have insignificant influence on the route. Engineering practice shows that SBAS-InSAR technology can effectively analyse the surface deformation of goaf sites and has important practical significance in optimising the selection of special geological routes in expressways.
1. Introduction
Construction projects above goaf deposits may lead to secondary movement and the deformation of goaf, which in turn may lead to the subsidence, local cracking, and tilting of structures at ground level. Goaf has become a major geological disaster that restricts the environment of mining area, sustainable development of the mining industry, and social economy. With the continuous expansion of mining intensity and scope, the subsidence areas around mines have become larger, thus greatly limiting the economic development of mining regions. Mining subsidence is reflected in the contradiction between the supply and demand of coal resource exploitation and economic development. At present, along with the rapid economic development in Nanjing, Hefei, and other surrounding cities, the pressure of highway construction has continued to increase. This amplifies the contradiction between soil and mines and has even led to a series of social problems. Also, troubles are often encountered during the quest to sustainably develop mining towns and resource-rich regions. If there is an expressway through the goaf and then highway subsidence and uneven subgrade subsidence, the subgrade and line may experience deformation, thus resulting in disastrous consequences. However, due to a lack of subsidence history data in goaf areas, only empirical evaluation methods can currently be used to assess stability in goaf areas. Therefore, existing methods for the stability evaluation and residual subsidence prediction of goaf areas cannot meet the requirements of engineering practice. It is of great practical significance to evaluate goaf stability, which can be achieved by monitoring the subsidence of goaf.
InSAR technology combines synthetic aperture radar (SAR) and electromagnetic wave interference (EMI) technologies and is a method for the topographic survey and digital elevation inversion [1]. Since its creation, InSAR technology has been widely used for disaster warnings due to its ability to penetrate the atmosphere and obtain ground elevation and deformation information in all weather conditions and at any time of day. Besides, it possesses a deformation measurement accuracy up to the millimetre level [2, 3]. Ye et al. [4] applied InSAR technology to landslide monitoring in the Three Gorges reservoir area to obtain deformation rate field. The Sentinel-1 SAR data were used by Dai et al. [5] for postdisaster assessment of the Xinmo landslide in Maoxian County, Sichuan Province in 2017. This confirmed the applicability of space-borne InSAR technology under complex weather and topographic conditions in mountainous areas. Based on the improved SBAS-InSAR technology, Wang et al. [6] established a local incident angle model that analysed the invalid areas of InSAR results. This enabled them to determine the deformation rate of landslides and plot a distribution map of potential landslides. Besides, Li et al. [7] correlated landslide deformation with groundwater level and precipitation intensity, based on various InSAR data in Yizi village, Jin Shajiang. Xue and Lv [8] applied the unscented Kalman filter method and Sentinel-1 SAR data to predict landslide deformation in Maoxian County, Sichuan Province, China. Results indicated that the combined time series InSAR technology and unscented Kalman filter can be used to predict landslide deformation before large-scale landslides occur. InSAR technology also has great potential in surface deformation monitoring. Ma et al. [9] used Sentinel-1, COSMO-SkyMed, and TerraSAR-X data to monitor multi-scale subsidence in the Guangdong-Hong Kong-Macao Greater Bay Area. They discovered that sediment consolidation is the primary cause of subsidence, while groundwater extraction and artificial building loads are the trigger factors for deformation. Farolfi et al. [10] used the global navigation satellite system (GNSS) to correct the results of PSInSAR, which was applied to the monitoring of land subsidence in Ravenna and Ferrara, Italy. Also, Rateb and Kuo [11] utilised the Sentinel-1 SAR time series interferometry, hydrology, and weather experimental data to determine the relationship between land subsidence and groundwater storage decline near Baghdad, Iraq. Huang and He [12] obtained surface deformation field and velocity field data for the Hexi area in Nanjing by using SBAS technology, which confirms that this method is effective for extracting surface subsidence data. In the Hua Tugou oilfield in Northern Tibet, Li et al. [13] acquired a time series of surface deformation by using small baseline interference pairs and carried out accurate deformation monitoring on the surface. It can be seen that InSAR technology has already been widely used in surface subsidence monitoring, indicating that it also has strong practicability and generalization for surface subsidence inversion in coal mining areas.
In this study, we apply short baseline subset time series analysis technology (SBAS-InSAR) to monitor the history of goaf deformation in the Huaibei section of the Xu Huaifu Expressway. By minimising the influence of coherence removal and elevation error, we can obtain accurate data of goaf deformation and provide a basis for the stability assessment of goaf.
2. Project Overview and Field Investigation
2.1. Project Overview
The Xuzhou to Fuyang highway, which is about 215 km long, has helped to realise more convenient transportation between Anhui province and Jiangsu province and has enhanced the highway network in the region. The research area is located in the Huaibei section of Xu Huaifu Expressway. A diagram of the Huaibei section of the Xu Huaifu Expressway is shown in Figure 1. Huaibei City, Anhui Province, is one of the top ten coal bases in China, and coal mining has helped the urban economy to develop rapidly over the years. However, after long-term natural compaction of the underground goaf formed after coal mining, the separation layer and cracks in the rock mass are still under compaction, and the subsidence of goaf has continued to occur for a long time.

2.2. Investigation of Goaf Status
There are numerous mining areas along the Huaibei section of the Xu Huaifu Expressway. According to the field investigation, the Liuqiao 1 and Liudong mines, at a mileage of K48 + 000-K50 + 000, may influence the route stability. Goaf formations at Liuqiao 1 and Liudong mines affecting the route are labelled C1–C5 and the plane position diagram of expressway route and goaf is shown in Figure 2.

Three goaf sites in Liuqiao 1 may have an impact on the proposed route and are designated as C1, C2, and C3. These sites have been generated by mining in the No. 4 and No. 6 coal seams. Coal seam No. 4 is located in the Lower Permian Shihezi Formation. The mining depth is between −450 m and −670 m, and the thickness of the coal seam is between 0.9 m and 4.0 m. The goaf site is formed between 2008 and 2015. Coal seam No. 6 is situated in the Lower Permian Shanxi Formation. The mining depth is between −570 m and −770 m, and the thickness of the coal seam is between 0.9 m and 5.0 m. The goaf site was formed between 2011 and 2016.
Zone C1 is located to the west of the route. The surface of Zone C1 displays obvious mining subsidence, and serious cracks have formed on the surface. As Figure 3(a) illustrates, a subsidence basin has formed on the surface and a pond has appeared. The walls of the surrounding houses have cracked, and when the subsidence basin is formed, the surface tilted, the nearby village has been relocated, and the vegetation consists of mainly small plants.

(a)

(b)

(c)
Zone C2 and Zone C3 are formed by mining coal in the No. 6 coal seam. The formation time of goaf in Zone C2 is between 2009 and 2011, the mining level is from −150 m to −350 m, and the average thickness of the coal seam is 3.3 m. The formation time of goaf in Zone C3 is from 2009 to 2014, and the mining depth range is from −250 m to −400 m. According to the field investigation, there is an obvious surface collapse in Zone C2, and whole block subsidence has occurred in the C3 area near the provincial highway, as shown in Figures 3(b) and 3(c).
There are two goaf formations in the Liudong coal mine, which may have an impact on the proposed route, namely, C4 and C5. The C4 zone is created by mining the No. 10 coal seam, the mining time is from 2008 to 2014, the mining depth is −110 m to −370 m, and the thickness of the coal seam is 2.8 m–3.4 m. The mining times and range of Zone C5 are relatively complex. Between 2015 and 2016, it has a great impact on the formation of goaf. The mining depth is −110 m to −350 m, and the thickness of the coal seam is 2.8 m∼3.4 m. Field investigations revealed that obvious collapse has occurred on the surface in Zones C4 and C5, and a subsidence basin with ponds had formed, as shown in Figure 4.

(a)

(b)
2.3. Qualitative Analysis of Goaf Subsidence
The mining time of the No. 4 and No. 6 coal seams in Liuqiao 1 is 2008–2016. Therefore, the goaf mining time is between 4 and 12 years, meaning that a majority of the goaf formation time is more than five years. The coal mine belongs to a state-owned coal mining company, thus the coal seam recovery rate of the working face is higher and the collapse zone and fault zone are sufficiently developed. Compared with similar goaf sites and according to relevant experience, an empirical value of surface movement duration is shown in Figure 5(a), and we estimate that goaf subsidence has reached 80–90%. Because the active period of goaf subsidence has ended and it is now in the degeneration stage, the qualitative estimate of residual subsidence is smaller. However, due to the large buried depth of goaf and the long duration of surface deformation, residual deformation will continue for a long time in the future. Therefore, the deformation is occurring slowly and the possibility of large-scale subsidence is small. Similarly, in the No. 10 coal seams in the Liudong coal mine, which is mined between 2007 and 2016 but most intensively towards the year 2016, the mining completion time is shorter, but the coal seam depth is shallower. According to a field investigation and visit, it is presumed that 80–90% of goaf subsidence has already occurred and that the active period of goaf subsidence has ended. Based on data collection and combined with the subsidence coefficient of remaining mobile deformation, we constructed a contour map of current mining and residual subsidence, as displayed in Figure 5(b).

(a)

(b)
In summary, the active stage of surface deformation caused by goaf in the main coal seams has ended. The goaf in the shallow coal seams (mining depth <400 m) is fundamentally stable, while the goaf that appears in the deep coal seams is in a declining subsidence period with little residual subsidence. However, the residual deformation in Zone C2 will continue for a certain period. Because part of the Huaibei section of the Xu Huaifu Expressway is planned to pass over the mining face, we must investigate the subsidence of the goaf to determine its stability and prevent the occurrence of potential risks.
3. Time Sequence Analysis of Surface Subsidence
In this study, using the latest interferometric synthetic aperture radar time series analysis method (SBAS-InSAR) and based on recent Sentinel-1 satellite SAR images from the European Space Agency, we monitored goaf deformation in the main areas surrounding the Huaibei section of the Xu Huaifu Expressway. This provided enough data to support the subsequent analysis of deformation genesis. The core concept of SBAS-InSAR is to network all images according to their time and space baselines. If the interference pairs of small baselines can be connected into a strongly connected graph, the deformation time series can be estimated by using phase changes and the least-squares method. The temporal deformation series of the surface can be obtained while DInSAR processing can reduce the influence of decorrelation and errors of elevation and atmosphere.
The deformation extraction process of SBAS-InSAR technology is shown in Figure 6 and includes the following steps: (1) generate connection diagrams, register the input Sentinel-1 data single view complex images, conduct interference pair combination with the set time and space baselines as thresholds, and carry out differential interference processing for each group of image pairs; (2) use ENVI SARscape radar interference processing software to complete the generated interference figure, interferogram deflattening effect, adaptive filtering, and phase unwrapping operation steps. Each interferogram operation generates interference results and unwrapping images to realize a preliminary understanding of the deformation of the study area. At the same time, transform the differential interferogram from a SAR coordinate system to a geographic coordinate system to confirm the location of suspicious deformation points. (3) Access to Variables: remove phase changes due to noise, atmospheric error, ground effect, and other factors to preserve only the deformation phase. After analysing geometrical relationships, perform a conversion on the variables to obtain the average annual figure for the deformation rate of the entire area. Then, convert the annual average deformation rate grid file to vector points in the software to acquire the deformation rate and history of each vector point.

To accurately evaluate the stability of the study area, we extracted four consecutive years of deformation subsidence data from June 2017 to June 2021, which route station at K46 + 000 to K52 + 000. The SAR image data of goaf in the Huaibei section of the Xu Huaifu Expressway were preprocessed to obtain an interferogram of the study area, as shown in Figure 7. There are differences between satellite remote sensing images in summer and winter of a year. Through the coherence analysis of satellite images in the project area, the image coherence is continuous. Therefore, it can be seen that the settlement images of the study area are small in summer and winter.

4. Time Sequence Analysis of the Subsidence Risk Area
To facilitate analysis, the form variables of the obtained deformation map are accumulated to create a surface subsidence change map from June 2017 to June 2021, as shown in Figure 8. In the figure, the areas outlined in red designate zones with large subsidence deformation, while the blue-outlined areas show locations with smaller subsidence deformation. Figure 8 illustrates that the accumulative total goaf subsidence over the four years has generally been stable. Considering the extent of subsidence, subsidence range, and influence on the route, the subsidence risk of the Y1-05 and Y1-09 areas is relatively high. Significantly, a traditional theoretical analysis did not locate the potential risk point of the Y1-05 area, and the stability judgement of this area is crucial for the selection of the route. According to the document “Technical Rules for the Design and Construction of Goaf Highways”, the subsidence rate stability threshold at a goaf site is 60 mm/a. The maximum subsidence rate of Y1-05 is 54.7 mm/a–58.8 mm/a, while the subsidence rate at Y1-09 is 50.7 mm/a–54.8 mm/a. Based on the above data, our analysis will focus specifically on the Y1-05 and Y1-09 areas.

4.1. Time Sequence Analysis of the Deformation Area at Y1-05
The closest distance between the Y1-05 area and the east side of route station at K46 to K47 of the expressway is about 273 m. The area Y1-05 is close to the east side of Da Miaogou Village, and its subsidence centre is located at the midpoint between Zhangdian Village and Da Lianglou Village, with a shape variable over 200 mm, as shown in Figure 9.

The Y1-05 area is a potential risk point. To study the deformation trend of the entire area, eight monitoring points are selected to assess subsidence in the form of a “cross” distribution in the Y1-05 area. The distribution of the monitoring points and their deformation rates are presented in Figure 10. The subsidence at Points 1 and 2 are the most serious, and the subsidence at these points accelerated from August 2018. The subsidence at Point 1 slowed down in June 2020, exhibiting total subsidence of 191 mm from June 2017 to June 2021. The subsidence at Point 2 started to decelerate in December 2020, and there is no increase in subsidence at either point from the end of 2020. The subsidence at the outermost monitoring points (Points 5 and 8) is the smallest, and the other monitoring points present a similar trend, respectively experiencing the stage of acceleration—gentle—acceleration—gentle. The subsidence of all of the monitoring points tended to be stable in the final half-year since December 2020. It can be seen that the centre of the Y1-05 sedimentation basin is located around Point 1 and Point 2 and gradually radiates outward in a funnel shape.

(a)

(b)
The accumulated subsidence value at Point 1 is the greatest, but the subsidence in the past year is approximately 25 mm, and only about 10 mm in the most recent six months. In contrast, the subsidence at Point 2 reached 30 mm in the second half of 2020, which is at the threshold of goaf stability evaluation. Therefore, further analysis is still required to assess the overall stability of the area.
4.2. Time Sequence Analysis of the Deformation Area at Y1-09
The Y1-09 deformation area is located to the west of the route station at K49 + 200 to K50 + 100, and the nearest distance to the centre line is 205 m. There is discernible subsidence in the core of the area. According to the goaf distribution map, Y1-09 is situated in the C2 zone, which is consistent with the predictions of field investigations and theoretical analysis. Figure 11 shows the Y1-09 area subsidence distribution diagram.

By extracting the data from eight monitoring points in the Y1-09 deformation area for time-sequence analysis, we find that the subsidence at Points 6 and 8 is the smallest. Points 1, 5, and 7 demonstrate similar levels of subsidence, while the subsidence at Point 3 is the largest. When the monitoring points with similar subsidence are connected, it forms an obvious funnel shape with Point 3 at the core. In the past year, the subsidence rate at all monitoring points has become more stable. At Point 3, the subsidence is only 18 mm in the year from June 2020 to June 2021, and no obvious subsidence is observed at any other point from December 2020, indicating that this area is generally stable. Monitoring of subsidence in the Y1-09 area is shown in Figure 12.

(a)

(b)
5. Deformation Characteristics of the Subsidence Risk Area
To judge whether the study area may have an impact on the planned route, we conducted subsidence basin characteristic analysis and subsidence fitting of the study area according to the subsidence monitoring data, to predict potential subsidence.
5.1. Distribution Characteristics of the Subsidence Basin
The subsidence observation lines are arranged in the middle of the Y1-05 and Y1-09 areas, and the bearing of the observation lines is perpendicular to the direction of the route. The study area runs from left to right with the subsidence centre as the midpoint.
As Figure 13 illustrates, the subsidence curve shows symmetrical distribution with the subsidence core as the symmetry point and it presents an inverted “convex” shape overall. The subsidence curve can be divided into the accelerated subsidence zone, stable subsidence zone, and undisturbed zone from the centre to either side. As the distance from the subsidence centre grows, the subsidence decreases in a nonlinear manner. Figure 13(a) shows that from June 2017 to June 2018, a point 300 m away from the centre in the Y1-05 can be classified as the stable subsidence area. The subsidence boundary radius is 700 m–800 m, and the subsidence centre is about 1 km away from the centre line of the route. Therefore, the subsidence basin will have little influence on the road. Over time, the Y1-05 basin has shown a trend of decreasing subsidence after initial growth in every year. The squat of all points of the Y1-05 area observation lines reaches a maximum between June 2018 and June 2019, and the maximum subsidence point is offset to the left. The zones of accelerated subsidence and stable subsidence increased accordingly, and the range of accelerated subsidence zone approaches 800 m away from the centre of subsidence. Since June 2019, the subsidence of each point in the observation line has gradually reduced, resulting in minimal subsidence in the past year. Although the Y1-05 subsidence area is not exhibiting accelerated subsidence, it is still slowly sinking, and the subsidence basin is experiencing gradient-type development.

(a)

(b)
Figure 13(b) reveals that the Y1-09 subsidence basin is similar in shape to Y1-05. Subsidence has reduced year by year, and the rate of subsidence has displayed a decreasing trend; particularly, since June 2019, the subsidence curve has displayed a similar pattern over the last two years. And from June 2020 to June 2021, according to the subsidence curve, the stable subsidence area is significantly greater than from June 2019 to June 2020. This confirms that subsidence has entered a stage of stability, and because the expressway is about 500 m away from the subsidence centre, the subsidence basin will have a negligible influence on the route.
5.2. Residual Subsidence Risk Area Fitting Prediction
Residual subsidence is also an important index that evaluates the stability of goaf. Combined with an analysis of the subsidence time sequence, the residual subsidence in the Y1-05 and Y1-09 areas can be determined. Based on the centre point data of the Y1-05 and Y1-09 areas from June 2017 to June 2021, the subsidence amount can be fitted:
As Figure 14 illustrates, the similarity of the fitted curves is above 0.98 and the fitting results are very close to the measured values. Besides, the centre point subsidence of the two areas shows negative exponential attenuation. The Y1-05 area central subsidence amount is expected to be 720 mm, and the current subsidence is 245 mm. Over the next 15 years, the residual subsidence in the centre point of the Y1-05 area is predicted to be 475 mm, with annual average subsidence of 31.67 mm. This confirms that the field is stable, although the Y1-05 area subsidence basin has a trend of subsidence with no subsidence mutation, it still needs stronger monitoring. The subsidence of centre point in the Y1-09 area is 340 mm, and the total remaining subsidence is only 50 mm. From Figure 14, it is obvious that the subsidence rate of the Y1-09 area has entered a stable stage and the residual subsidence will have a negligible influence on the expressway.

5.3. Prediction of Subsidence Risk Area Site Expansion
In order to further evaluate the subsidence development range of the risk area sites, according to the shape of the subsidence basin, the model function is used to fit the shape of the subsidence basin in the Y1-05 area and Y1-09 area respectively. There are four unknown numbers, namely, , , , and , where determines the extreme value of function, determines the symmetry axis of the function, determines the opening amplitude of the function, and is the fine adjustment of the value of the function. The fitting results are shown in Figure 15 and Table 1.

(a)

(b)
From Table 1, we can see that the correlation coefficients of the fitting results are above 0.95, which indicates that the function model can describe the shape of the subsidence basin well. Therefore, this function model, combined with the predicted value of central point subsidence, can predict the development of future subsidence basins.
Taking the Y1-05 area as an example to illustrate the solution process, it is assumed that the shape function of the subsidence basin within the next 15 years is . As mentioned before, the basin has a significant range of accelerated subsidence, stable subsidence, and undisturbed areas. The subsidence basin in the Y1-05 area is sinking slowly, and the symmetry axis of the subsidence basin and the range of accelerated subsidence approaches a fixed value after June 2019. As shown in Figure 15(a), the subsidence slope is close to a constant in the stable subsidence area of the subsidence basin, and the slope of the stable subsidence area in different years is approximately parallel. Therefore, it can be assumed that the symmetry axis of the future subsidence basin and range of the accelerated subsidence area will remain unchanged, and the subsidence rate of the stable subsidence area will be close to a constant. Combined with the above fitting equation, we let the symmetry axis of the Y1-05 area subsidence basin is −30, and the following relationship is obtained:where x0 is the position of the symmetrical axis of the subsidence basin, is the subsidence rate of the stable subsidence area at the left of the subsidence basin, is the subsidence rate of the stable subsidence area at the right of the subsidence basin, and is the predicted value of the subsidence central point within the next 15 years.
Similarly, the shape function of the subsidence basin within the next 15 years in the Y1-09 area can be obtained. The shape functions of the subsidence basin in the Y1-05 area and Y1-09 area within the next 15 years, respectively, are as follows:
The prediction of development subsidence basin in Y1-05 area and Y1-09 area is shown in Figure 16. The stable subsidence area of the subsidence basin will be extended within the next 15 years. The most dangerous of these is the Y1-05 area, whose subsidence basin extends 100 metres from the edge of the line, despite the subsidence not affecting the stability of the route. However, in order to ensure the safety of the route and prevent uneven subsidence of the subsidence basin, which may lead to pavement tension shear failure in this section of the route, and it can be treated with flexible reinforced mesh.

6. Conclusion
(1)The SBAS-InSAR technology is utilised to rapidly obtain the ground surface subsidence of the goaf site in the Huaibei section of the Xu Huaifu Expressway from June 2017 to June 2021. Afterwards, the subsidence situation of the goaf site is accurately analysed, solving the problem of lacking historical data on surface subsidence in expressway route selection.(2)There is still residual subsidence on the surface, even after underground mining has ceased in the Huaibei section of the Xu Huaifu Expressway. Using SBAS-InSAR technology, we find that Y1-05 and Y1-09 are two areas with high subsidence risk. The Y1-05 area is a potential risk point, while the settlement of Y1-09 area is affected by goaf.(3)The prediction of the subsidence boundary radius of the Y1-05 area is between 300 m and 350 m, and the residual subsidence is expected to be 475 mm, highlighting the need to strengthen monitoring. The residual subsidence of the Y1-09 area is 50 mm and has entered a gentle stage.(4)Based on the subsidence value of the centre subsidence point, the development trend of subsidence basin within the next 15 years is predicted. The shape function of subsidence basin is highly consistent with the subsidence data of the subsidence basin, which provides strong support for the prediction of the subsidence development of the subsidence basin.Data Availability
The datasets used in the experiments and discussed in the paper are available from the corresponding author on reasonable request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research was funded by the National Natural Science Foundation of China Youth Project (51804006) and State Key Laboratory of Deep Coal Mining Response and Disaster Prevention and Control Independent Scientific Research (SKLMRDPC20ZZ04).