An Improved Adaptive Template Size Pixel-Tracking Method for Monitoring Large-Gradient Mining Subsidence

The monitoring of large-gradient deformation caused by coal mining is of great significance to the prevention and management of disasters in mining areas. The interferometric synthetic aperture radar (InSAR) method captures the small-gradient ground deformation on the edge of the subsidence basin accurately but is unreliable for capturing large-gradient deformation.The intensitybased pixel-tracking method (e.g., the normalized cross-correlation (NCC) method) can overcome the limitations of InSAR’s maximum detectable displacement gradient and incoherence. However, the pixel-tracking method is sensitive to template size. It is difficult to estimate ground subsidence accurately by the conventional pixel-tracking method with fixed template size. In this paper, the signal-to-noise ratio (SNR) is redefined and an improved locally adaptive template sizemethod is proposed by identifying optimal template adaptively based onmaximization of the redefined SNR.The constraint radius is used to constrain the search area in this improved method. The frequency of misrepresentation is reduced by finding the peak of the correlation coefficient surface within the search area. Both simulation data and real ground subsidence data are used to test this algorithm.The results show that this method can improve monitoring accuracy compared with the traditional pixel-tracking method for fixed template size.


Introduction
China is not only the world's largest coal consumer but also the largest coal producer.Accounting for 42% of global production, coal accounts for more than 74% of China's domestic energy consumption.Over 96% of the coal output in China comes from underground coal mining.Underground coal extraction inevitably causes ground subsidence, which severely disturbs the ground surface [1][2][3][4].Subsidence in coalfields has already become one of the most prominent social and environmental problems in China.
Mining-induced subsidence has been studied using numerous ground-based methods, including levelling, total station surveys, and GPS field surveys, which are primarily based on point-by-point measurements.These methods are time-consuming, labor-intensive, costly, and limited in spatial coverage and precision.With the rapid progress in Earth observation techniques in recent years, different interferometric synthetic aperture radar (D-InSAR) and time-series InSAR methods such as PS-InSAR [5][6][7][8], SBAS-InSAR [9][10][11][12][13][14][15][16], TCP-InSAR [17][18][19][20], and ITPA-InSAR [21,22] have been proven to be effective methods with an accuracy of millimeters to centimeters for mapping ground deformation induced by various natural and anthropogenic causes.However, those methods are based on phase unwrapping and have significant limitations based on deformation gradient.For large deformation gradients, the phase unwrapping will fail [23][24][25].The ground deformation caused by coal mining has the characteristics of a large deformation gradient.The InSAR technique based on phase unwrapping can only obtain the deformation information of edge of the subsidence basin.
Although NCC has been documented as being reliable for measuring mass surface displacements [37,44,45], NCC is sensitive enough to monitor subsidence related to mining.A small template size results in misrepresentation due to the amount of noise [46].A larger template size results in greater computation time and deformation distortion [47].The optimum template size strikes a balance between noise reduction and distortion reduction.In order to obtain accurate monitoring results, template sizes need to be locally adaptive as the deformation varies spatially.Based on the image signal-to-noise ratio (SNR), a locally adaptive template size method is proposed for matching repeat images of Earth surface mass movements by Debella-Gilo and Kääb [46].Debella-Gilo and Kääb's method has a drawback for monitoring large-gradient deformation by SAR image pair: a lot of gaps appear in displacement vectors due to the presence of a large amount of noise in the SAR image.
In this paper, an improved locally adaptive template sizes pixel-tracking method is proposed for monitoring largegradient subsidence caused by mining.In order to reduce the appearance of misrepresentation, a constraint radius determined by the maximum of subsidence caused by mining is used to limit the search area.The peak of the crosscorrelation coefficient surface is found in the search area.The SNR is redefined and the optimal template is determined based on the maximization of redefined SNR.Simulated data and real ground subsidence data are used to test this improved algorithm.

Calculation of Normalized Cross-Correlation (NCC)
Coefficients Surface.The NCC algorithm is a measurement of similarity widely used in image matching to measure the similarity of templates between a master image and slave image.A higher cross-correlation coefficient means a higher similarity [44].The deformation can be determined based on the position of the maximum value of the correlation coefficient [45].As shown in Figure 1, the master image over the mining area is taken in time 1 and the slave image over the same area is taken in 2.The slave image is registered to master image.The ground has been deformed sometime between 1 and 2; the pixel (, ) moved to ( − ,  − V).In order to acquire the displacement in  direction, the twodimensional cross-correlation coefficient matrix should be obtained by cross-correlation coefficient calculation between the reference temporary templates, which takes every pixel in search template as central pixel.The reference template size is defined as 2 + 1 pixels, where  represents the number of pixels on either side of the central pixel (, ) in both dimensions.The NCC coefficient  between the reference and temporary templates in the slave image is calculated by the following equation: ; here,  is the mean intensity in the reference template and  is the mean intensity in the temporary template.The sizes of the temporary and reference templates are equal.The value of  ranges between 0 (when the matching entities are not exactly the same) and 1 (when the matching entities are exactly the same).

Definition of SNR Ratio and Constraint
Radius.The cross-correlation coefficient is mainly dependent on template size [46,47].For small template size, it is difficult to distinguish the signal correlation from noise; as the template size increases, more and more pixels are included, the noise is suppressed, and the cross-correlation coefficient began to increase.When noise is optimally suppressed, the template is optimal.With further increasing template size, the maximum NCC coefficient starts to decrease because additional pixels mean more dissimilarity as a result of increased geometric distortion.Therefore, the optimal template exists within a certain interval.A template that is too large or too small is not suitable for pixel tracking.In order to obtain the optimal template size, the SNR is introduced as a judgment index.
Greater SNR value means that noise is suppressed.SNR is calculated as follows: where  max is the maximum cross-correlation coefficient in the areas centering on the center pixel of the corresponding template and  pixels for the constraint radius; and  is the mean value of the cross-correlation coefficients in 5 × 5 pixel areas centered on  max .SAR intensity images are well known to be noisy due to radar speckle noise.The use of large template does not optimally remove the influence of this noise.This large amount of noise is likely to cause significant misrepresentation in the SAR imagery.Because SAR images are registered and the maximum ground deformation caused by mining can be estimated using geological data, misrepresentation may be reduced by constraining the radius to limit the search area for the maximum correlation coefficient.A subsidence of more than twenty meters is extremely rare in coalmines.In order to reduce the frequency of misrepresentation, for RADARSAT-2 SAR images with a pixel size of 2.66 m, constraining the radius to five pixels is optimal.For different SAR images with different pixel sizes, the constraining radius should change according to the size of SAR image and the computing equation (3) as follows: where  is the constraining radius;  is the maximum value of vertical ground subsidence; it can be calculated by the probability integral method according to geological data; generally speaking, 20 m is enough for mining subsidence;  is the incidence angle, and  range is the size of the SAR image in the range direction; int( ) is an integer operation.

Implementation of Adaptive Template Sizes.
A greater SNR means a more accurate location of the peak of correlation coefficient can be acquired.The SNR is related to the template size, but the relationship between SNR and template size is difficult to model mathematically.Perhaps the most effective way to acquire the optimal template size is to calculate the SNR of different template sizes and find the maximum value of SNR.In this section, a method based on locally adaptive template sizes is described.For every pixel in the monitoring area, assume that the template size interval is [, ], where  is the smallest value of the template size and  is the largest.In general, for ground deformation caused by coal mining, the template interval from 21 to 121 is sufficient.The algorithm is implemented as follows.
Step 1.The SNR is computed using the following equation: where  is the initial template size step (usually 8),  = (0, 1, 2, . . ., ), + represents the template sizes,  max (+) is the maximum value of the correlation surface in areas centering on the center of the corresponding template with a size of  +  with a radius of five pixels, and  (+) is the mean value of the 5 × 5 pixel template centered on the maximum correlation coefficient in the correlation surface.Then, the maximum value of the initial SNR is calculated as follows: where SNR1 is the maximum value of the initial SNR corresponding to the template size temp1 and max{⋅ ⋅ ⋅ } represents the maximum value of the components in the braces.
Step 2. The template size step is decreased to /2; only two SNRs with template sizes of (temp1 − /2) and (temp1 + /2) are computed using (4), and the maximum value of the correlation surface is calculated by the following equation: where SNR2 corresponds to the template size temp2.
Step 3. The template size step is decreased to /4; the two SNRs with template sizes of (temp2 − /4) and (temp2 + /4) are computed by (4), and the maximum value of the correlation surface is calculated by the following equation: where SNR3 corresponds to the template size temp3, which is the optimal template size.Once these three steps have been completed, a method converting oversampling of the crosscorrelation surface into a higher resolution surface using twodimensional interpolation algorithms and fitting of a secondorder polynomial to the correlation surface around the peak is used to match the exact positions of two corresponding pixels [46].Assume that the pixel level (integer) position of the peak of cross-correlation surface has been calculated as ( 0 ,  0 ); this pixel has two neighbors in the range direction: ( 0 − 1,  0 ) and ( 0 + 1,  0 ).In order to find the exact position of the peak, ( 8) is used to calculate the noninteger location of the peak, which will be added to the integer pixel location: where Δ is the noninteger location of the peak, ( 0 − 1,  0 ), ( 0 ,  0 ), and ( 0 + 1,  0 ) are the cross-correlation coefficients corresponding to the pixel positions of ( 0 −1,  0 ), ( 0 ,  0 ), and ( 0 + 1,  0 ), respectively.After determining the exact position of the peak of cross-correlation coefficient surface, the displacement in the range direction can be expressed as follows: where  range is the displacement in the range direction,  exact is the exact position of peak of cross-correlation coefficient surface,  central is the position of central pixel in reference template, and  range is the size of SAR image in range direction.

Study Area and Data
In order to test the performance of this improved algorithm, both simulated image pair named simulated data with known displacements by probability integral method for each pixel and real ground deformation caused by mining in the Daliuta mining area named real ground subsidence data are used in this paper.The simulated SAR image was made by resampling the master RADARSAT-2 SAR image with pixel spacing of 2.66 m in the range direction and of 2.88 m in the azimuthal direction.The simulated ground subsidence is predicted by probability integral method (shown in Figure 2).According to the incident angle and pixel size of RADARSAT-2, the vertical subsidence is converted into subsidence in the range direction and transformed into offsets measured in pixels.Based on the offsets, the bicubic convolution interpolation function is used to resample the master RADARSAT-2 image.The resampled image (slave image) contains the complete information of ground subsidence predicted by probability integral method.The simulated RADARSAT-2 data (slave image) is only a resampling of master RADARSAT-2 image, so there is no time-related noise in simulated data.
The real ground subsidence data composed of two RADARSAT-2 SAR images generated on 27 November 2012 and 27 March 2013 was collected from an area located in Daliuta mining area, Shenmu County, Yulin City, in the north of Shaanxi Province, which is one of the largest areas of coal exploitation in China.It belongs to a desert-foothills region with a dry climate, long winters with little snow or rain, and short summers.The two working faces in the study area, 52304 and 22307, are shown in Figure 3(a).Working face 22307 was completely exhausted by the pillar-and-room method in 2005, and residual ground subsidence was small.Therefore, this working face made little contribution to the ground subsidence captured by RADARSAT-2 images between 27 November 2012 and 27 March 2013.A fully mechanized caving technique was used on working face 52304, and the rate and spatial extent of subsidence in this face were large.The mining conditions of working face 52304 were as follows: altitude of the ground surface was between 1154.8 and 1269.9 m, mining depth was 235 m, and average thickness of the coal seam was 6.94 m.There are 45 GPS monitoring points in the strike direction and 27 in the dip direction (shown in Figure 3), measured between 27 November 2012 and 27 March 2013.GPS points were used to test this improved locally adaptive template sizes pixeltracking method.

4.1.
Results of Simulated Data.Two registered simulated RADARSAT-2 SAR images were analyzed using the improved locally adaptive template size pixel-tracking method.The results of this analysis are shown in Figure 4.The simulated RADARSAT-2 SAR images are processed by the fixed template size pixel-tracking method with template sizes of 31, 61, 121, 181, and 241.In order to evaluate the performance of this method, two hatched lines in the strike and dip directions were extracted (Figure 5); the root mean squared error (RMSE), the mean absolute value of difference (MAVD), the maximum absolute value of difference, and the minimum absolute value of difference between the simulated subsidence and monitoring results for different fixed template sizes are computed and shown in Table 1.
Summaries within the results of the fixed template size pixel-tracking method for simulated data are as follows.(i) When the template size is small (i.e., 31), monitoring results show a serrated profile, especially where the subsidence is large.(ii) With an increased template size, the profile becomes   smoother and closer to the simulated subsidence profile (i.e., 61), but when the template is too large (i.e., 121 or 181), deformation compression is present.The larger the template is, the more severe deformation compression is.(iii) Compared to the fixed template size pixel-tracking method, the improved locally adaptive template size pixel-tracking method can better monitor ground subsidence.

Results of Real Ground Subsidence.
The improved locally adaptive template size pixel-tracking method and three fixed template sizes of 61, 91, and 121, were used to examine the subsidence caused by mining (Figure 6).The RMSE and the MAVD between the monitoring results and the GPS results were computed in both the strike and dip directions (Table 2).In order to standardize the comparison, 45 GPS results in strike direction and 27 GPS results in dip direction were converted to the range direction according to the RADARSAT-2 parameters.As shown in Table 2, the locally adaptive template size pixel-tracking method was more accurate than fixed template size pixel-tracking method.Overall, the locally adaptive template size pixel-tracking method was more accurate in the dip direction as opposed to the strike direction.From point 11 to point 18 in the strike direction, the MAVD is large, and it is hard to monitor the subsidence exactly using the pixeltracking method (Figure 6(a)).This is because the extensive cracks in the ground surface (Figure 7) lead to apparent changes in the ground surface which will cause dramatic changes of intensity in the slave SAR image.The dramatic changes of intensity cause significant reduction of the crosscorrelation coefficient in ground crack affected areas, leading to a significant reduction in monitoring accuracy and even to mismatch.

Effect of Template Size on the Pixel-Tracking Method.
Template size is a crucial factor in accurately analyzing subsidence with a given area.Small template sizes commonly misrepresent the total amount of subsidence, while large  template sizes cause compression of deformation.The larger the template size is, the more pronounced the compression of deformation is.In order to analyze the influence of template size on pixel tracking, take the pixel in the center of subsidence basin as an example; the correlation coefficient surfaces with the template sizes of 11, 31, 61, 121, 181, and 241 are calculated, respectively, and shown in Figure 8.
When the template size is small (e.g., 11), multiple peaks appear in the correlation coefficient surface (Figure 8(a)), which will produce error within the pixel tracking.This is due to the fact that small template sizes cannot effectively suppress the noise within the SAR image.As the template size increases to 31, a clear single peak in the correlation coefficient surface appears (Figure 8(b)), which is mainly due to the noise being effectively suppressed.When the template size increases to 61, the single peak is still clear but a small uplift is present around the peak.As the template size continues to increase to 121, the shape begins to broaden and phenomenon of passivation appears (Figure 8(d)).With the further increase of template size, more and more pixels are included, which makes the correlation coefficient around the peak increase, so that the peak shape shows the passivation phenomenon.As the increase of the template size continues (e.g., 181 and 241), the passivation phenomenon is more obvious (Figures 8(e) and 8(f)), and correlation coefficients around the central pixel of the reference template are closer to the peak, so that the location of the quadratic polynomial fitting peak is closer to the center of the reference template, resulting in deformation compression.
The optimal template size produces the surface with an obvious correlation coefficient peak with the noise effectively suppressed (Figure 8(b)).In this paper, SNR is used as an important indicator to find the optimal template size.According to the definition of SNR (in Section 2.2), the larger the SNR value is, the more obvious the correlation coefficient peak is, which is conducive to the precise monitoring of deformation.The relationships between SNR and template size in the simulation data and real data are shown in Figure 9.The relationship between SNR and template size is difficult to express in mathematical terms, but it exhibits some patterns.In the beginning, the SNR increases with the increase of template size, when the SNR reaches the maximum value; the SNR begins to decrease as template size continues to increase.After the SNR decreases to a certain value, the SNR value remained stable.The optimal template size exists in a range, according to the processors' experiences; the template interval from 21 to 121 is sufficient for mining subsidence.The template size larger than 121 will cause compression of deformation, while template size smaller than 21 will cause low level of SNR due to the fact that the noise cannot be effectively suppressed in SAR image pair.

Accuracy Analysis of Simulated and Real Data.
The improved locally adaptive template size pixel-tracking method is used to monitor the simulated data and the real ground subsidence from the Daliuta mining area.The monitoring accuracy of the simulated data is higher than that of the real data (Table 3).This is mainly because the simulated data is only a resampling of the master image according to the predicted results of deformation by the probability integral method, in addition to the speckle noise contained within the SAR image.In the real data, the time interval between the capture of the two images is 120 days, which will result in a large time-related noise.Also the attitudes of the satellite cannot be exactly the same when both images were captured, introducing some geometric-related noise into the     calculation.This noise will significantly affect the accuracy of the pixel-tracking method.Rapid excavation of the coalmines resulted in serious ground surface ruptures (Figure 7) in the strike direction.In areas where ground surface ruptures are severe, the intensity of the SAR image changes drastically, and even if the improved locally adaptive template size pixel-tracking method is used, surface deformation cannot be obtained correctly.Therefore, the monitoring accuracy in the dip direction is significantly better than that in the strike direction in 52304 working face of Daliuta mining area (Table 3).
According to the statistical analysis of the optimal template size of the simulated data and the real deformation data (Figure 10), the optimal template size for the real data is significantly larger than the optimal template size for the simulated data.This is mainly due to the different levels of noise in simulated data and the real deformation data.In the real data, it is difficult to obtain the appropriate template size  due to the presence of a large amount of noise.Increasing the template size is the most effective way to suppress this noise, but larger template size will result in significant deformation compression (e.g., the template size of 121 both in simulated data and in real data shown in Figures 5 and 6).The strategy based on maximization of redefined SNR can better balance noise and deformation compression, resulting in more accurate acquisition of deformation monitoring results (e.g., locally adaptive algorithm shown in Figures 5 and  6).
This improved locally adaptive template size pixeltracking method increases the accuracy of the monitoring results and avoids the selection of a fixed template size according to processor's experience; however, this algorithm contains some defects.The main drawback is that it is computationally expensive to find the optimum template size, being more than thirty times slower than using a fixed template size.For example, to monitor the subsidence in 140× 140 pixel area on a computer with 6 GB RAM and a 2.4 GHz processor speed, the fixed template size of 61 took about 6 minutes, while the improved locally adaptive template size algorithm took about 3 hours in MATLAB.The process can take up to several hours depending on computer performance and the number of templates to be processed.In the future, we should attempt to optimize this algorithm and improve the calculation efficiency.

Conclusions
The classic pixel-tracking method based on normalized crosscorrelation can overcome the limitations of the maximum detectable displacement gradient that InSAR is subject to, but the accuracy of fixed template size method is severely affected by template size, especially in areas with a small area of subsidence and large subsidence gradient.In this paper, an improved locally adaptive template size pixel-tracking method based on maximization of redefined SNR is proposed.By constraining the radius, the search region where the correlation coefficient peak appears is constrained, which greatly reduces the frequency of misrepresentation.According to experiments carried out in the Daliuta mining area and using simulated data, the accuracy of simulated data in the strike direction is 0.063 m and is 0.047 m in the dip direction; the accuracy of the real ground subsidence data in the strike direction is 0.381 m and is 0.151 m in the dip, which are significantly higher than the fixed template size pixel-tracking method.In general, a higher spatial resolution means higher monitoring accuracy for SAR images by pixel tracking.The improved locally adaptive template size pixeltracking method should perform better when monitoring large-gradient mining subsidence using higher spatial resolution SAR images.

Figure 1 :
Figure 1: Schematic diagram of pixel tracking.The slave image is registered to the master image; 2 + 1 represents the size of template. direction indicates the range direction (the direction in which the displacement occurs), and  direction represents the azimuth direction.

Figure 2 :
Figure 2: Simulated subsidence by probability integral method (maximum value is 2.6 m).

Figure 3 :
Figure 3: Location of the study areas.In (a), the red frame represents the areas covered by the RADARSAT-2 images, green points are GPS monitoring points, and the white arrow represents the direction of mining.(b) is the interferogram generated by two RADARSAT-2 SAR images on 27 November 2012 and 27 March 2013.Because of the limitations of deformation gradient and incoherence, the D-InSAR technique based on phase unwrapping can only obtain the deformation information in edge of the subsidence basin.

Figure 4 :
Figure 4: Results of improved locally adaptive template size pixeltracking method for monitoring simulated ground subsidence.Two red hatched lines in the strike and dip directions are used to evaluate the performance of improved method.

Figure 5 :
Figure 5: Profiles of results from pixel tracking with different template sizes: (a) profiles of different template sizes results and real subsidence in the strike direction; (b) profiles of different template sizes results and real subsidence in the dip direction.

Figure 6 :
Figure 6: Comparison of different template sizes and locally adaptive template sizes with GPS data: (a) in the strike direction (points from 11 to 18 are shown in the black dashed box; the RMSE and MAVD are larger due to ground crack); (b) in the dip direction.

Figure 7 :
Figure 7: Results of improved locally adaptive template size pixel-tracking method in SAR geometry.Black points represent GPS monitoring locations and there are 45 GPS monitoring locations in strike direction and 27 GPS monitoring locations in dip direction.The distance between two adjacent GPS points is about 20 meters.The black frame denotes an area with numerous ground cracks; even improved adaptive template size pixel-tracking method in this area cannot get accurate monitoring results.

Figure 9 :
Figure 9: Relationship between template sizes and SNR on the central pixel of subsidence basin in simulated data and real data.

Figure 10 :
Figure 10: Histograms of optimal template sizes in monitoring real and simulated subsidence: (a) real subsidence; (b) simulated subsidence.

Table 1 :
Accuracy of different template sizes in the pixel-tracking method for monitoring simulated ground subsidence in the strike and dip directions.

Table 2 :
Accuracy of different template sizes in the strike and dip directions.

Table 3 :
Accuracy of simulated data and real data in the strike and dip directions.