Automatic Detection of Subsidence Troughs in SAR Interferograms Using Mathematical Morphology

Abstract

In this paper, an automatic algorithm for the detection of subsidence areas in SAR interferograms is proposed. It is based on the analysis of spatial distribution of the interferogram phase, and its coherence and entropy. The developed method was tested for differential interferograms generated on the basis of Sentinel-1 SAR images covering mining areas in South Poland. The obtained results were compared with those achieved using a method based on circular Gabor filters. Performed analysis revealed that the detection rate for the proposed method varied from 34% to 83%. It is an improved method based on Gabor filters that achieved a detection rate from 30% to 53%.

1. Introduction

Differential interferometric SAR (DInSAR ) analysis is a significant source of information in crisis management dealing with the problem of ground deformations [1,2,3]. DInSAR results can be used for monitoring landslides, volcanic activity [4], earthquakes [5], and mining-induced subsidence [6]. Nowadays, a very large amount of SAR data acquired by Sentinel-1 is freely available, making differential radar interferometry more accessible. However, to take real advantage of DInSAR in crisis management, it is necessary to ensure its automatic processing, which is a subject of many works [7,8,9]. A system that allows for the continuous monitoring of endangered areas and can alert about critical values of ground deformations is needed. Such a solution requires a technique for the automatic identification of subsidence patterns in SAR interferograms. There are many scientific papers on this problem that use different methods. In [10], the automatic detection of subsidence troughs in SAR interferograms based on circular Gabor filter was proposed. Its detection rate varied from 30% to 53% with a relatively low number of false alarms. In [11], circlet transform was compared with Hough transform [12,13]. There are also works that are based on a convolutional neural network [14,15]. Each of those methods has its limitations, which is why their efficiency in trough detection in noisy satellite images is low. Methods [10,11,12,13] use circular or elliptical troughs in the detection process. In many cases, the shape of the troughs on interferograms is deformed, which significantly reduces the effectiveness of these solutions. Moreover, the proposed method is much less computationally expensive, which is quite important due to the large size of the processed SAR images. However, most of the solutions [10,11,12,13,14,15] can be used to support the detection of subsidence in larger SAR images, but they are not suitable for fully automated operations.

Here, a new algorithm for the automatic detection of subsidence areas in radar interferograms is proposed. It identifies elliptical interferometric fringes that usually correspond to subsidence troughs. The developed method is intended for monitoring mining and postmining ground deformation. When elliptical fringes are caused by other deformation phenomena or atmospheric effects, they are also recognized. Therefore, the interpretation of the obtained results should be connected with a priori information about the studied region.

In the proposed algorithm, mathematical morphology methods are used for analyzing the phase, coherence, and entropy of interferograms, and for the identification of subsidence troughs. To speed up detection, the proposed algorithm consists of two stages. In the first stage, potential areas of subsidence troughs occurrence are identified. Afterwards, those areas are verified. The proposed method was tested on interferograms generated on the basis of Sentinel-1 SAR data.

2. Methodology

The proposed algorithm consists of two main steps. The first step eliminates from further analysis parts of the interferogram where, there are most probably no subsidence troughs. A map of the interferogram phase is used with maps of its coherence and entropy (calculated for data converted into unsigned int 8-bit types with Equation (1) in a moving window). For the phase map, lines with a sudden phase change (approximately from $-\pi$ to $+\pi$) are identified. Thos is motivated by the fact that such phase changes are characteristic for interferometric fringes that occur in the case of subsidence troughs represented by a full phase cycle. For this reason, subsidence troughs without full phase cycle may not be included in further analysis. This needs to be improved in the future.

(1)$Entropy=-\sum _{n=0}^{255}{p}_n·{log}_2{p}_n$

where $p_n$—sum of pixels with intensivity n divided by the number all pixels in window [16].

The identified lines that are very short are removed from the next steps of analysis. The continuity of the remained lines is improved by the closing operation with a linear structural element [17]. Areas inside those lines constitute a set of potential regions, where the centers of subsidence troughs may occur. An example of the obtained results is presented in the Figure 1 (obtained on the basis of the generated interferogram using Sentinel-1 data).

From the map presented in Figure 1, all regions with low coherence (noisy regions) and low entropy (homogeneous regions) were removed by considering two logical images presenting coherence larger than a defined threshold (${\gamma }_{th}=0.6$) and entropy above a defined threshold ($H_{th}=0.82$), respectively. To define the appropriate thresholds, tests were performed on several interferograms (different than those for which tests of the method were performed). Those interferograms, generated on the basis of Sentinel-1 images, differed from each other by time interval (6 and 12 days), and a number of subsidence troughs (from 1 to several). For the analyzed values of thresholds, the number of correctly detected troughs was divided by the number of total alarms. Thresholds that gave the highest value of this proportion were chosen as optimal. A procedure for the automatic selection of thresholds (${\gamma }_{th}$ and $H_{th}$ ) should be proposed in the future.

For a map with potential regions where the centers of subsidence troughs may occur, and logical maps of coherence and entropy, multiplication is applied. In Figure 2, the obtained results (areas surrounded by white lines) are presented on the phase image. Those areas represent more reliable regions (compared to those presented in Figure 1) where centers of the subsidence trough may occur (called perspective areas). The location and extent of perspective areas are an input for the second step of the detection procedure.

The second stage of the detection algorithm is much more time-consuming than the first one. For each pixel belonging to perspective areas, the center of the theoretical ellipse was set. The method of obtaining values the of theoretical ellipses radius is presented in the next section and Figure 3. After that, from this center, the rays in each direction (with 5^o step) were created. The length of the rays was set up while taking into account the resolution of the analyzed images and the information about the probable maximal spatial extent of the subsidence trough (in our case, for Sentinel-1 data, this was set to 75 pixels). The correlation between the phase values along those rays and a saw function was calculated, and is presented in Figure 4 and described in Equation (2).

(2)$\rho \left(r\right)=\pi -2\pi ·\frac{mod (r,\omega )}{\omega }$

where r—distance to the center of the theoretical ellipse, $\omega$—a period of saw function, $\mathrm{mod}$—modulo operation. The saw function represents the phase changes characteristic for the interferometric fringes. The shape and period of the saw function are based on a typical period in troughs in tested images.

After that, for each ray, distance $d_{max}$ (from the center of the theoretical ellipse) to the maximal value of the calculated correlation was determined. The line connecting the $d_{max}$ values for each ray determines the shape of the subsidence pattern (Figure 3—dotted line). The maximal and minimal values of dmax were used to define the minor and major axes, and orientation of the theoretical ellipse (Figure 3—solid line).

In the last step of the algorithm, for each pixel from perspective areas, the fitting error ($d_{max}\left(k\right)$) between the distance of the maximal correlation function ($d_{max}$) for each radius and theoretical ellipse ($d_{ellipse}$) was calculated according to Equation (3).

(3)$d_{error}=\sqrt{\sum _{k=1}^{Number of radius}{\left(d_{max}\left(k\right)-d_{ellipse}\left(k\right)\right)}^2}$

Values of $d_{error}$ were normalized and inverted to the range of <0,1>, where 1 means ${\forall }_k:d_{max}\left(k\right)=d_{ellipse}\left(k\right)$. Lastly, only those pixels for which the error function was greater than a threshold (in our case, 0.85) were identified as subsidence areas. The lower the value of the threshold was, the greater the number of false alarms. The greater value of the threshold can suggest that some subsidence areas could be undetected.

3. Experimental Analysis

The developed method was applied to a dataset covering a selected region of the Upper Silesian Coal Basin (USCB) in South Poland. USCB is located in southern Poland and in the Ostrava-Karvina region in the Czech Republic. Mining activity there has been conducted for over 200 years. The USCB region is Poland’s largest hard-coal basin and one of the largest hard-coal deposits in Europe. It is also a large metropolitan region with 37 towns and nearly 3 million residents. It is estimated that an area of around 600 km² already suffers from subsidence in USCB [18].

This is the same test site that was used in [10], where the detection method based on circular Gabor filters was presented and tested. For the selected region, the subsidence phenomenon is characteristic. The ground deformations are caused by extensive and long-lasting coal exploitation. The occurrence of subsidence troughs was confirmed many times using classical levelling or DInSAR analysis [19,20]. In the DInSAR method, a pair of SAR images that correspond to the same area and had been acquired at different times by the same SAR system were used. First, the analyzed SAR images need to be coregistrated, and the phase difference ($\Delta \varphi$) is then calculated for each pixel of the image. In this way, the flattened interferogram is obtained. To remove the component associated with the topography of the studied region, the Digital Elevation Model (DEM) is used. DEM is converted into the synthetic interferogram, which is then subtracted from the flattened interferogram. As a result, the differential interferogram is obtained. Its phase values can be directly linked to ground deformations. The subsidence phenomenon, especially related to mining activity, usually manifests itself on the differential interferograms as elliptical interferometric fringes. The density of the fringes is correlated with the size of the deformations.

In the presented work, as in [10], the efficiency of the presented subsidence trough detection algorithm was tested using three differential interferograms generated on the basis of Sentinel-1 SAR data. The interferogram concerns the following time spans: 10–22 October 2016 (interferogram I), 10 October–3 November 2016 (interferogram II) and 22 October–3 November 2016 (interferogram III). All interferograms were prepared using SNAP. The interferograms were filtered twice, using the Goldstein and Lee filtering methods. Goldstein filtering was applied before the removal of the topographic phase, whereas the Lee filtering method was applied to differential interferograms. The parameters of both methods are shown in Table 1. The multilooking procedure was not applied during processing step. More information about the used datasets and their processing parameters can be found in [10]. The main stages of the subsidence trough detection algorithm were performed using the Mathworks MATLAB 2019b computing package.

The proposed algorithm of the automatic detection of subsidence areas was applied to both interferograms computed using VV and VH channels of SAR images. Results are presented in Table 2 and Figure 5, Figure 6 and Figure 7. Table 2 $k1$ shows a number of undetected troughs that had a very irregular shape. $k2$ is a number of areas incorrectly classified as subsidence troughs (false alarms).

Table 2 shows that, for all analyzed cases, the best results were obtained for VV channels. In the case of VH channels, fewer subsidence troughs were detected, and more false alarms were noticed. Taking into account the results of both polarimetric channels (VH and VV), the number of false alarms could be significantly reduced (even more than sevenfold for interferogram I). However, in such cases, the number of correctly detected troughs can also decrease more than twice the time (interferogram II). The best results of subsidence area detection were obtained for interferogram I (Figure 5). In this case, for VV channels, 26 subsidence troughs (86.7%) were correctly identified, and 30 false alarms were created. Table 3 shows that, for interferogram I, the results for VH channels are only slightly worse than for VV channels. Therefore, it is beneficial to take into account both channels. VH and VV results are characterized by a high number of correctly detected subsidence troughs (25, which accounts for 83.3% of all troughs) and a low number of false alarms (only 4 cases).

The detection results for interferogram II (VV channel) are presented in Figure 6. The number of correctly identified subsidence troughs was high and totalled 25 (86.2% of all troughs). However, for interferogram II, the number of false alarms was very high, equal to 86 cases. Additionally, results for the VH channel were significantly worse, and false alarms decreased by using a logical product of VV- and VH-based results, compound detection accuracy.

For the third interferogram (VV channel), the proposed method allowed for correctly detecting 18 of 24 (75%) subsidence troughs (Figure 7). The number of false alarms was 25, which may be significantly reduced to 9 cases when the results of VH and VV are taken into account. In such a case, the number of undetected troughs increased from 6 to 7.

Table 2 also shows that, for all analyzed interferograms, the number of undetected troughs with irregular shapes ($k1$) and a number of false alarms with a shape similar to subsidence pattern ($k2$) were very low, and in both cases did not exceed 20%. In the case of a method based on a circular Gabor filter [10], the $k1$ and value of $k2$ reached 60%. For the algorithm proposed in this paper, the irregularity of the subsidence pattern is thus not a crucial problem in the detection procedure. Figure 8 shows the example of two undetected subsidence patterns. The subsidence pattern presented in Figure 8 (left image) was not detected because of its strong irregularity in different directions. This irregularity produced very high error function values. Figure 8 (right) shows the subsidence pattern, where an incomplete interferometric fringe pattern occurred (due to the relatively low ground subsidence magnitude—less than half of the radar wavelength). Figure 9 shows examples of false alarms from interferogram I. They appeared in regions where there were no similar patterns to those characteristic for subsidence troughs. Such regions should be removed from analysis in the first step. However, probably due to a not restrictive enough threshold on coherence and entropy, they passed to the second step of the analysis.

Table 3 shows the comparison between obtained results using the proposed algorithm and those from a method based on circular Gabor filters [10].

Table 3 shows that, for interferograms I and III, the number of correctly detected troughs by the proposed method was higher than that in the method based on circular Gabor filters. However, for interferogram III, the number of false alarms also increased. In the case of interferogram II, the proposed method gave worse results than those of the method based on circular Gabor filters. For this interferogram, the number of correctly detected troughs was similar for both methods, but the number of false alarms in the method proposed in this paper was significantly higher.

4. Conclusions

This work proposed a new algorithm for the automatic detection of subsidence areas in SAR interferograms using a method based on mathematical morphology. The proposed method was tested on three interferograms generated on the basis of Sentinel-1 SAR data. Detection results were compared with those obtained using a method based on circular Gabor filters. The detection rate for the proposed method varied from 34% to 83%. This is an improvement over the method based on circular Gabor filters that achieved detection rates from 30% to 53%. However, for both methods, different numbers of false alarms were reported. For the proposed method to be more operational, additional work should be devoted to the automatic selection of thresholds and dealing with subsidence patterns without a full phase cycle.

Author Contributions

Conceptualization, M.D., S.P.-S. and J.S.; methodology, M.D.; software, M.D., H.M. and R.M.; validation, S.P.-S., J.S., H.M. and R.M.; formal analysis, M.D., H.M. and R.M.; investigation, M.D., S.P.-S. and J.B.; resources, S.P.-S. and J.S.; data curation, M.D., A.F. and J.B.; writing—original draft preparation, S.P.-S. and M.D.; writing—review and editing, S.P.-S., M.D., A.F. and J.B.; visualization, M.D., H.M. and R.M.; supervision, A.F. and J.B.; project administration, A.F. and J.B.; funding acquisition, S.P.-S., M.D., A.F. and J.B. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by the AGH—University of Science and Technology, Faculty of Geology, Geophysics and Environmental Protection, as a part of statutory project. The work is a part of MineSAR Project.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data accessed upon request from any of the author.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Identified lines with potential (white) regions where centers of subsidence troughs may occur.

Figure 2. Location and extent of perspective areas presented in interferogram phase.

Figure 3. Theoretical ellipse and calculated shape of subsidence pattern.

Figure 4. Example of saw function with period [Forumla omitted. See PDF.] = 15 px.

Figure 5. Results of subsidence area detection for interferogram I (time period: 10–22 October 2016; used VV polarimetric channel).

Figure 6. Results of subsidence area detection for interferogram II (time period: 10 October–3 November 2016; used VV polarimetric channel).

Figure 7. Results of subsidence area detection for interferogram III (time period: 22 October–3 November 2016; used VV polarimetric channel).

Figure 8. Example of two undetected troughs from interferogram I.

Figure 9. Example of two false alarms from interferogram I.

Table 1

Parameters of filtering methods.

Goldstein Filtering	Lee Filtering
FFT size: 64 × 64	window size: 7 × 7
Window size: 3 × 3	target window size: 3 × 3
Coherence threshold: 0.2	sigma: 0.7
Adaptive filter exponent: 1.0

Table 2

Results of subsidence trough detection. $$k1$$—number of undetected troughs; $$k2$$—number of false alarms.

Interferogram I: 10–22 October 2016
Polarimetric	Detected	Undetected	False
Channel	Troughs	Troughs	Alarms
VV	26	4	30
VH	25	5	35
VH and VV	25	5 (k1 = 20%)	4 (k2 = 0%)
Interferogram II: 10 October 2016–3 November 2016
Polarimetric	Detected	Undetected	False
Channel	Troughs	Troughs	Alarms
VV	25	4	86
VH	10	19	102
VH and VV	10	19 (k1 = 10%)	39 (k2 = 7%)
Interferogram III: 22 October 2016–3 November 2016
Polarimetric	Detected	Undetected	False
Channel	Troughs	Troughs	Alarms
VV	18	6	25
VH	17	7	34
VH and VV	17	7 (k1 = 14%)	9 (k2 = 11%)

Table 3

Comparison of obtained results by the proposed algorithm and the method based on circular Gabor filters.

	Number of Troughs	Method Based on Circular Gabor Filters		Method Proposed in This Paper
		Detected Troughs	False Alarms	Detected Troughs	False Alarms
Int. I	30	16	21	25	4
Int. II	29	12	5	10	39
Int. III	24	7	0	17	9

References

1. Gupta, R.P. Remote Sensing Geology; 3rd ed. Springer: Berlin, Germany, 2017.

2. Sedighi, M.; Arabi, S.; Nankali, H.R.; Amighpey, M.; Tavakoli, F.; Soltanpour, A.; Motagh, M. Subsidence Detection Using InSAR and Geodetic Measurements in the North-West of Iran. Proceedings of the 2009 Fringe Workshop; Frascati, Italy, 30 November–4 December 2009.

3. Porzycka, S.; Strzelczyk, J. Preliminary results of ground deformations monitoring within mining area of “Prosper-Haniel” coal mine. Proceedings of the 12th International Multidisciplinary Scientific GeoConferences (SGEM); Balchik Municipality, Bulgaria, 17–22 June 2012; Volume 2, pp. 895-899.

4. Lu, Z.; Dzurisin, D. InSAR Imaging of Aleutian Volcanoes: Monitoring a Volcanic Arc from Space; Springer Science and Business Media: Chichester, UK, 2014.

5. Stramondo, S.; Chini, M.; Bignami, C.; Salvi, S.; Atzori, S. X-, C-, and L-band DInSAR investigation of the April 6, 2009, Abruzzi earthquake. IEEE Geosci. Remote Sens. Lett.; 2010; 8, pp. 49-53. [DOI: https://dx.doi.org/10.1109/LGRS.2010.2051015]

6. Nádudvari, A. Using radar interferometry and SBAS technique to detect surface subsidence relating to coal mining in Upper Silesia from 1993–2000 and 2003–2010. Environ. Socio-Econ. Stud.; 2016; 4, pp. 24-34. [DOI: https://dx.doi.org/10.1515/environ-2016-0003]

7. Deguchi, T. Automatic InSAR processing and introduction of its application studies. Proceedings of the 26th Asian Conference on Remote Sensing; Hanoi, Vietnam, 7–11 November 2005.

8. d’Oreye, N.; Celli, G. Automatic InSAR systematic processing and web based tool for efficient data mining: Application to volcano monitoring in Africa. Proceedings of the “Fringe 2009 Workshop”; Frascati, Italy, 30 November–4 December 2009.

9. Albino, F.; Biggs, J.; Yu, C.; Li, Z. Automated Methods for Detecting Volcanic Deformation Using Sentinel-1 InSAR Time Series Illustrated by the 2017–2018 Unrest at Agung, Indonesia. J. Geophys. Res. Solid Earth; 2020; 125, 017908. [DOI: https://dx.doi.org/10.1029/2019JB017908]

10. Porzycka-Strzelczyk, S.; Rotter, P.; Strzelczyk, J. Automatic detection of subsidence troughs in SAR interferograms based on circular Gabor filters. IEEE Geosci. Remote Sens. Lett.; 2018; 15, pp. 873-876. [DOI: https://dx.doi.org/10.1109/LGRS.2018.2815782]

11. Bała, J.; Dwornik, M.; Franczyk, A. Automatic subsidence troughs detection in SAR interferograms using circlet transform. Sensors; 2021; 21, 1706. [DOI: https://dx.doi.org/10.3390/s21051706] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33801252]

12. Klimczak, M.; Bała, J. Application of the Hough transform for subsidence troughs detection in SAR images. SGEM; 2017; 17, pp. 819-826.

13. Yip, R.K.; Tam, P.K.; Leung, D.N. Modification of Hough transform for circles and ellipses detection using a 2-dimensional array. Pattern Recognit.; 1992; 25, pp. 1007-1022. [DOI: https://dx.doi.org/10.1016/0031-3203(92)90064-P]

14. Rotter, P.; Muron, W. Automatic Detection of Subsidence Troughs in SAR Interferograms Based on Convolutional Neural Networks. IEEE Geosci. Remote Sens. Lett.; 2021; 18, pp. 82-86. [DOI: https://dx.doi.org/10.1109/LGRS.2020.2966079]

15. Pei, J.; Huang, Y.; Huo, W.; Zhang, Y.; Yang, J.; Yeo, T. SAR Automatic Target Recognition Based on Multiview Deep Learning Framework. IEEE Trans. Geosci. Remote Sens.; 2018; 56, pp. 2196-2210. [DOI: https://dx.doi.org/10.1109/TGRS.2017.2776357]

16. Gonzalez, R.C.; Woods, R.E.; Eddins, S.L. Digital Image Processing Using MATLAB; Prentice Hall: Hoboken, NJ, USA, 2003.

17. Dougherty, E.R. Mathematical Morphology in Image Processing; Marcel Dekker Inc.: New York, NY, USA, 1993.

18. Graniczny, M.; Kowalski, Z.; Tkowska, A.P.A.; Przyłucka, M. Observation of the Mining-Induced Surface Deformations Using C and L SAR Bands: The Upper Silesian Coal Basin (Poland) Case Study. Lecture Notes in Earth System Sciences; Springer: Berlin/Heidelberg, Gemany, 2013; pp. 249-255.

19. Perski, Z. Applicability of ERS-1 and ERS-2 InSAR for land subsidence monitoring in the Silesian Coal mining region, Poland. Int. Arch. Photogramm. Remote Sens.; 1998; 32, pp. 555-558.

20. Mirek, K.; Isakow, Z. Preliminary analysis of InSAR data from south-west part of Upper Silesian Coal Basin. Mineral Resour. Manag.; 2009; 25, pp. 239-246.