Source and pattern identification of ground deformation based on the framework of Blind 1 Source Separation--Nonnegative Matrix Factorization: a case study in a long-term GPS 2 monitoring mine

29 Ground deformation caused by mining affects the safety of underground mining and buildings. 30 High-resolution characterization of ground deformation respond to latent sources is a first step 31 toward improved hazard forecasting. We rely on long-term GPS displacement data and an improved 32 NMF algorithm of BSS to identify the sources driving ground deformation. The NMF identifies 33 three sources (S1, S2, and S3) with distinct spatial and temporal surface deformation patterns and 34 quantitatively reveals the contribution of each source to ground deformation. S1 captures horizontal 35 Y-displacement related to horizontal tectonic stress σ 1 . S2 dominates vertical Z-displacement with 36 a widespread isotropous deformation driven by self-weight body force. S3 controls horizontal X- 37 displacement related to horizontal tectonic stress σ 2 . Besides, we find that independent 38 contributions of these three sources can be resolved from the GPS data. The results show that the 39 sharp change of source contribution and even transformation of dominate source are highly related 40 to the severe deformation belt. This phenomenon is time independent, and helps us to select sites 41 and find the potential risk area at its early stage.


Introduction 45
Ground deformation is an common phenomenon in terms of civil engineering, groundwater 46 pumping, earthquake, landslide and underground mining operations (Cao et  there are multiple sources (factors) that can cause rock movement, such as topography, type of rock, 50 structure of rock mass, tectonic stress, mining technology, and co-seismic vibration effects, etc. The 51 combination of these sources and the dominant source vary significantly depending on the local 52 geological conditions. Without greater confidence in how these proposed factors affect ground 53 deformation near surface, we cannot reliably use the results of geodetic analyses to infer shallow 54 deformation, nor can we formulate models to predict deformation in future events. 55 In recent years, these sources have been widely studied by various methods to have a deep 56 understanding of the mechanism of ground deformation. There are four broadly categories. 57 First, numerical model based on finite element or discrete element method. It is a forward 58 method that takes some sources (e.g., stress condition, joint distribution, and mining operations) and fixing other parameters at the same time. The second category is classical statistical methods (e.g., 65 correlation analysis, factor analysis, regression analysis, etc.). These methods try to find the law of 66 ground deformation from observed data (Hui et al. 2018). They also can help to find out parts of the 67 sources leading to ground deformation. However, the hidden sources of massive data with high 68 dimensions cannot be easily interpreted by using these methods (Sahu et al. 2017; Song et al. 2011). 69 The third category is analytical solution based on empirical methods, profile function, and influence 70 function(Díaz-Fernández et al. 2010). These methods take the geometrical and mechanical 71 parameters as input to predict the ground deformation. The type and value of parameters are highly 72 dependent on the prior knowledge of researchers. However, it is unreasonable to predict before we 73 know the mechanism of each potential source causing ground deformation. The forth category is 74 machine learning consists of a series of supervised learning algorithms (Chen et  identified by prior knowledge and then to fit the observed data. Therefore, these methods are not 77 self-constraint by data and much dependent on the knowledge of researchers. 78 However, none of these studies try to objectively and quantitatively find out that 1) what kind 79 of and how many sources leading to the ground deformation? 2) what is the independent pattern of 80 each source? 3) how much is the proportion of each source contributing to the deformation at a 81 specific site? and 4) how the dominant source transfers with time among these sources? 82 Mechanism and pattern vary significantly among different sources when causing ground 83 deformation. As a result, the characteristics of ground deformation induced by different sources can 84 be also different in terms of horizontal and vertical displacement. In plain words, even though we 85 did not know who they (sources) are, how many of them there and how they did (causing ground 86 deformation), but what they have done was recorded in the data. Based on this logic, it is possible 87 to find a backward way to separate the sources from the long-term monitored displacement data. If 88 sources can be successfully decomposed from observed data, hidden pattern of each source can be 89 revealed independently with knowledge of petrology, tectonic geology, fracture mechanics, and rock 90 mechanics, etc. Most importantly, it is the transformation of dominant sources with time in a 91 concerned area that is essential to the prediction of ground deformation and safety management. 92 Indeed, sources identification causing temporal and spatial development of ground 93 deformation is still challenging. In recent years, source identification has received extensive 94 attention in areas such as signal processing, biomedical engineering, data analysis and data mining. 95 Nonetheless, the theory and applications are still being developed. Source identification can be 96 complex because 1) some of source signals may mutually interfere, 2)there may be attenuation 97 (linearly or nonlinearly) of signals when transmitting through the medium or discontinuity, 3) the 98 ratios of signal to noise vary in real word, 4) some of sources may be partly collinear or have similar 99

patterns. 100
The observed displacement data can be seen as the mixture signal of different displacement 101 patterns with some unknown mixing processes. From this point of view, we can introduce Blind 102 Source Separation (BSS) in the field of signal processing to find the hidden sources and their 103 patterns. The aim of BSS is to process observations acquired by sensor arrays in such a way that the 104 original unknown source signals are extracted by various algorithms without knowing or with 105 limited information about the characteristics of the transmission channels through which the sources 106 propagate to the sensors. Therefore, BSS can be used to reveal hidden sources and patterns in a 107 mixed data set. 108 The basic framework of BSS can be described as a linear mixture model by equation 1. 109 Where ∈ × is observed data (mixtures). ∈ × is an unknown mixing matrix 111 representing the linear combination of the sources and each row of ∈ × is a source. 112 represents additive noise or interference introduced during mixing and transmission. 113 Typically, the only known information in BSS is the mixtures . We need to determine the 114 number of sources and the solution of mixing matrix and , i.e., mathematically, we need to 115 solve equation 2 116 = min || − || 2 (2) 117 It is noteworthy that such a problem is ill-posed and has an infinite number of solutions. Thus, 118 extra imposed constraints (e.g., statistical independence, non-negativity, sparsity, maximum mutual 119 information, etc.) and regularized optimization procedures are often required to minimize the 120 objective function and to find an optimal and robust solution. Different constraints are usually 121 imposed on and based on particular cases and may lead to different results. Many methods 122 are included in BSS (e.g., independent component analysis (ICA), sparse component analysis (SCA), 123 non-negative matrix/tensor factorization (NMF/NTF), etc.). They all can successfully identify the 124 hidden sources and patterns from observed data set. However, it is the physical meaning and 125 interpretation of particular case that curial to determine which method to be used. and temporal surface deformation components (sources). However, the performance of ICA is not 131 good in some areas (e.g., in bio-signal and genomic signal) because the latent sources are partially 132 or entirely correlated. In other words, it is hard and often impossible to verify the statistical 133 assumptions on the sources (Tangirala et al. 2007). Under this situation, much effort has been 134 devoted to find an alternative approach with weaker assumptions. SCA started to emerge at the end 135 of the 1990s and became prominent during the 2000s (Chang et al. 2006; ).Unlike 136 ICA, SCA is suited to correlated sources and tries to explain data as a mixture of sparse 137 components (Karoui et al. 2012). In practice, both non-negativity and sparsity are often desirable or 138 necessary when the hidden sources have physical meaning. However, results of ICA and SCA often 139 lead to the subtraction in order to reconstruction of observed data. As a result, it is difficult to 140 interpret the physical meaning of sources. 141 In contrast to ICA, NMF enforces a non-negative constraint exhibiting some degree of natural 142 sparsity on the mixing matrix and source matrix (Lee and Seung 1999). In other words, NMF 143 does not allow subtraction and combinations, and, therefore, it is often used to quantitatively 144 describe the parts that comprise the entire entity. Furthermore, matrix factorization methods that 145 exploit non-negativity and sparsity constraints usually lead to estimation of the hidden sources with 146 specific structures and physical interpretations, in contrast to other BSS methods. For example, it is 147 the sources that contribute to the development (rather than reduce) of ground deformation that gains 148 our most concerns. Sparsity constraint significantly limits the number of sources to represent the 149 observed ground deformation data. Non-negativity constraint gives the physical interpretation of 150 each source on the deformation at monitoring sites. 151 In this study, the highlight is that without knowing other subsurface conditions and other prior 152  Here, we give the mathematical form of this particular problem. 176 Where × is observed data by m GPS sites and n monitoring campaigns. × is an 178 unknown mixing matrix representing the linear combination of p sources. × is a source 179 contribution at each discrete time moment.
× represents the measurement errors or noises. The 180 number of hidden sources is definitely less than the number of GPS monitoring sites (i.e., ≤ ). 181

182
In order to estimate matrix and , a cost function is chosen to quantify the similarity 183 measure (also referred to as divergence of dissimilarity, distance depending on the probability 184 distribution of estimated sources and the structure of data) between reconstructed data and observed 185 data. || • || 2 represents the Frobenius norm which is the squared Euclidean distance. and 186 are nonnegative regularization parameters (typically, In order to impose sparsity, we choose to use ℓ 1 -norm to minimize the cost function 192 To minimize the cost function , we choose the multiplicative update algorithm (equation 6, The iterative algorithm is stopped when there is little or no improvement of cost function 203 between successive iterations, = 10 −10 . 204

206
Obviously, the final solution is almost a local minima because the global minima is seldom 207 achievable. In order to achieve a reasonable and optimal solution, another three issues (the first issue 208 is to develop a problem dependent cost function) in NMF should be stressed: 1) alleviation of the 209 rotation indeterminacy, 2) determination of the number of sources, 3) initialization of and . 210

Alleviation of the rotation indeterminacy: 211
The eq. 3 can also be presented as It is easy to find many 212 rotational matrixes such that × =̂ and − × =̂ are also the solutions of eq. 3. 213 Therefore, the solution of NMF is non-unique. To overcome this, we should keep the input data , 214 decomposed matrixes and sufficiently sparse (Rickard and Cichocki 2008). One common 215 way is to normalize the input data to make and zero-grounded. Besides, it is a necessary 216 step to find a meaningful approximation of the hidden pattern for ground deformation neglecting 217 the scaling ambiguity. 218

Number of sources: 219
Remember that the number of sources (rank of , p) and mixing process is unknown in NMF. In this work, a heuristic method proposed by Minka (2001) was used to calculate the 230 smoothness index defined as: 231 is the eigenvalue of the covariance matrix for observed data in a 233 descending order ( 1 ≥ 2 ≥ ⋯ ≥ > 0). The sample variance is defined as: 234 Finally, the number of sources p can be selected by the following criterion: 236 = arg min =1,2,…, −2 ( ) (11) 237

Initialization: 238
The intrinsic of multiplicative update algorithms is nonconvex, even though the cost function 239 is strictly convex with respect to one set of variables (i.e., either or ). Therefore, it is impossible 2) Run the proposed algorithm in this paper for each set of initial matrices. 250

251
All the possible solutions are evaluated on the basis of criteria in terms of relative error (e), 252 orthogonality (o) and sparsity (u) of identified sources (Hoyer 2004). Each of these criteria reveals 253 the capability of the demixing method from a specific aspect. The final optimal solution is a tradeoff 254 among these three criteria, and this can be defined as a Euclidean distance in a 3D rectangular 255 coordinate system. 256 = min(|| || 2 ) , = 1,2, … , * ( + 1) (12) 257 Where is a vector composed of the three parameters( , , ). , and are defined as 258 the following equation 13 to 15 with a monotonic property. The value of , and ranges 259 [0, +∞), [0, +∞), and [1, +∞), respectively. A lower value of each parameter represents a better 260 performance in terms of a specific aspect. 261

Data processing 299
GPS monitoring network contains a reference net and a deformation monitoring net on the 300 ground surface. Reference net containing 7 benchmarks is set far away from the deformed area. 301 During the long-term period, some of points are destroyed. Therefore, a total of 78 monitoring points 302 with successive time series out of 101 monitoring points are selected in this study. These points 303 distribute along the exploratory lines with a spacing of 50 m between two points on the same line. 304 Besides, an additional monitoring line (H line) is laid out on the upper surface of the ore body, which 305 is perpendicular to the exploratory lines (Fig. 2b). Original data is received by Z-12-type GPS 306 receiver and antenna (Ashtech, USA) (Fig. 2c). The nominal accuracies of the horizontal and 307 vertical displacements are 3 mm+0.5 ppm and 5 mm + 1 ppm, respectively. Square error of 308 processed displacement data is about ±1.96 mm. For a single monitoring point, displacement data 309 comprises of the increment of X, Y, Z (local coordinate system) and each of them has 30 records. 310 Besides, displacement data is offset without changing the structure of data to ensure there is no 311 negative value in the matrix. 312

Results 313
The GPS monitoring data is a mixture of a series of unknown sources with dynamic mixing 314 ratios. Our goal is to separate the spatial and temporal patterns of dominate long-term and common-315 mode. The NMF does this by minimizing cost function with non-negative and sparsity constraints. 316 To seek a trade-off between interpretability and statistical fidelity, the formulated evaluation criteria 317 helps to select the optimal solution from all possible solutions that has the highest resolution of each 318 source pattern. 319 After 100 tries of NMF and a series of optimization criteria by using the parameters of Table  320 1, we obtain the 100 solutions with different accuracies (Fig. 3). The optimal solution is determined 321 by error (e), orthogonality (o) and sparsity (u) according to eq. 12 with the number of iteration steps 322 of 1245. The value of e, o, and u is 0.02, 11.14, and 1.06, respectively. Although the NMF cannot 323 give the accurate solution, the 100 tries show a validation in statistics. The reconstructed data with 324 less noises (Fig. 4c) fits well with original data (Fig. 4a). 325 The final solution consists of three sources (Fig. 4b), namely source 1 (S1), source 2 (S2), and 326 source 3 (S3), which optimally isolate the major spatial and temporal trends underlying the ground 327 deformation time series. For comparison, the amplitude is normalized by the relative maximum. It 328 is worth noting that these three sources are not orthogonal because NMF does not assume the 329 orthogonality of sources. Therefore, the sum of source contribution is not equal to one according to 330 the equation 16. The source contribution to original data (Table 2) shows that S1 has the highest 331 contribution of 0.53, and S3 has the secondary highest contribution of 0.51, and S2 has the lowest 332 contribution of 0.28. 333 S1 is a source that has a high impact on Y displacement and a low impact on X and Z 334 displacement. In contrast, S3 shows a high contribution to X displacement and a low contribution 335 to Y and Z displacement. S2 primarily contributes to the Z displacement and has a low contribution 336 to X and Y displacement (Fig. 4b). 337

Spatial pattern of source 338
The mixing matrix A represents contribution of the three sources at different monitoring point. 339 Therefore, each column of the mixing matrix A can be regarded as the spatial pattern of one source. 340 As the element of A is fixed, the spatial pattern of each source is stable among different monitoring 341 campaigns. This means that these three spatial patterns are time-independent and not to be affected 342 by mining intensity and mining methods. The value of displacement in different monitoring 343 campaign is controlled by the magnitude of source transient. Fig. 5 shows the contour of each 344 column of matrix A that represent the spatial pattern of S1 (Fig. 5a), S2 (Fig. 5b), and S3 (Fig. 5c), 345 respectively. 346 S1 captures the feature of Y displacement. In general, the value of Y displacement decreases 347 from northwest to the southeast. There are two displacement centers that is nearly divided by line

Physical meaning of identified source 362
In study area, high tectonic stress has been empirically considered as an important factor to the 363 ground deformation (Ma et al. 2015). However, the mechanism of ground deformation induced by 364 tectonic stress has never been revealed and quantitatively supported. Herein, according to the 365 identified patterns by NMF, we can make a forward step to deeply understand the process about 366 how tectonic stress affects ground deformation. 367 In natural conditions, the maximum principal stress σ 1 is perpendicular to the strike of ore 368 body (Ma et al. 2015). The long axis of strain ellipse is NW-SE and the relationship of stress is σ 1 369 ＞σ 2 (Fig. 6b). In post-mining condition, we find that the Z displacement at the center of ground 370 deformation in S2 pattern can be simplified as two sources controlled by S1 and S3 because the 371 contribution of S2 is nearly zero. Therefore, the ellipse shape at the center region (Fig. 5b) can be 372 thought as strain ellipse and the long axis is NE-SW. The rotation of ellipse indicates that mining 373 operations has caused the stress rotation in a larger area. At the same time, the relationship of stress 374 changes to σ 2 ＞σ 1 . This can be supported by deformation traces. In detail, by overlying S1 and S3 375 patterns, we find a set of "x" shape strong deformation belts on the surface (Fig. 6a). Through the 376 deformation traces, it is evident that the deformation belt by S3 (left wing) shows a shear property, 377 and the slide direction can be obtained to reveal σ 2 ＞σ 1 . The two quasi-circular closed region at 378 each end of shear belt may be caused by stress concentration (Fig. 5c) (Guo et al. 2020). Deformation 379 belt by S1 (right wing) shows an extension property because the displacement contour line in cross 380 section of "x" is straight and parallel to each other. In fact, as the mining tunnel is perpendicular to 381 σ 1 , broader free faces cause a larger stress release of σ 1 than σ 2 . This results in a series of extension 382 fractures perpendicular to σ 1 on surface (field investigation can be seen in Ma et al. (2015)) and 383 constraint of Y displacement by increasing the normal stress on joint planes. As a result, σ 1 causes 384 a large scale of X displacement at the beginning, and, then, σ 2 becomes the dominate stress in the 385 disturbed stress field which dominates Y displacement. This results in that S1 and S3 explain 386 different volume information of original data (Table 2) with the explanation ratios of ξ 1 /ξ 3 =1.05. It 387 is worth noting that the sum of explanation ξ 1 , ξ 2 , and ξ 3 is larger than 1 because these three sources 388 are not orthogonal causing the existence of partial correlation. Therefore, we conclude that S1 is 389 driven horizontal principal stress σ 2 , and S3 is driven horizontal principal stress σ 1 . 390 S2 describes a homogeneously and broadly distributed pattern that represents the self-weight 391 body force caused by gravity. In the scale of study area, density of rock mass is the only factor affect 392 the value and homogeneity of self-weight body force. We find that the mechanism of self-weight 393 body force on ground deformation is significantly different from empirical knowledge. Specifically, 394 on one hand, surface deformation induced by self-weight body force explains the least information 395 of original data compared to S1 and S3. On the other hand, self-weight body force has a high 396 contribution to Z displacement far away from center while low contribution near the subsidence 397 center (Fig. 5b). Compression state of rock mass is found at the center of ground deformation in 398 many studies Xia et al. 2016). In this study, high confining pressure caused by σ 1 399 and σ 2 produces enough force of friction on joint planes with high dip angle that alleviating the Z 400 displacement. At the same time, Z displacement is dominated by σ 1 and σ 2 because of the 401 existence of the joint planes with lower dip angle. This can be supported by ratios of source 402 contribution on single point such as 14~7 and H~11 (Fig. 9c). Besides, we propose that the high 403 contribution of S2 away from the center may be caused by the extension-subsidence after the stress 404 release of σ 1 and σ 2 . 405

Temporal pattern of source 406
Temporal evolution of displacement is critical for disaster early warning. However, previous 407 studies focus on the displacement evolution responding to the combined actions of ambient 408 sources (Li et al. 2020b). In contrast, in this study, temporal evolution of each displacement 409 component is related to separate source. Fig. 7 shows the temporal pattern of three sources. The 410 temporal pattern does not show periodicity or seasonality related to precipitation indicating that the 411 precipitation induced consolidation deformation can be ignored in this arid region. Another 412 important conclusion from the analyses is that the time lag is not obvious by comparing the curve 413 of three sources (Fig. 7). This indicates these three sources are independent in physical mechanism. 414 The deformation rate of three sources slowly increases with time. In detail, S1 induces the 415 largest increment of deformation rate of 0.0094, and S3 induces a medium value of 0.0086, and S2 416 causes a lowest increment of deformation rate of 0.0056 (Fig. 7). The different increment of 417 deformation rate may be driven by the different physical process. The stress concentration is more 418 remarkable with mining depth increase, and, thus, the displacement controlled by σ 1 and σ 2 419 increases slowly. In contrast, the self-weight body force has no obvious change because it is 420 controlled by gravity and the volume and density of overlying rock mass. Therefore, the increment 421 of deformation rate caused by S2 would be a stable value ideally. However, the discontinuities, 422 especially the joint and other structure planes, contribute to a component of Z displacement when 423 suffering horizontal stress. As the dynamic of stress redistribution, the normal stress may decrease 424 on structure planes at the tensile zone mainly distributed at the edge of deformation area (Fig. 5b)  425 can favor the Z displacement driven by self-weight body force. 426 By comparing the source contribution to X, Y, and Z displacement, we can divide the temporal 427 pattern into three stages (Fig. 8). These three stages reflect the evaluation of ground deformation. A 428 "vibration-acceleration" phenomenon is revealed. During the vibration period, the stress 429 redistribution and adjustment of discontinuities occur in rock mass. The acceleration is a relatively 430 faster process, which is driven by the stress concentration and the degradation of rock mass. Another 431 phenomenon is that each source has a more distinct contribution to its controlling displacement 432 component while keeps a continuous decrease contribution to another two components. A 433 reasonable explanation is that the overlying rock mass is becoming loose, and the confining pressure 434 is decreasing on the plane of discontinuities. Therefore, S3 contributes more X displacement and 435 less Y, and Z displacement along the plane (S1 and S2 have the same mechanism with S3). This  these studies make a forward to reveal what drives the deformation. As we all know, ground 447 deformation is always drove by multi sources in reality. Prediction may be more valuable and 448 reliable on the basis of physical mechanism of each source. In this study, we assume that the GPS 449 data is a mixture of a series of unknown sources with dynamic mixing ratios. As a result, we separate 450 two horizontal tectonic stresses and self-weight body force controlling the characteristics and 451 tendency of ground deformation. 452 On the basis of separated sources, we, firstly, propose a new conception "controlling pattern" 453 that is defined as a combination of sources ordered by their contribution at a certain time. It is 454 dynamic with time and any source can become the dominate source. It is known that different value 455 of confining pressure can cause different shape and location of rupture in the triaxial compression 456 test. Like that, different controlling patterns will lead to different features of ground deformation. In 457 addition, the location of micro deformation by first controlling pattern will play an important role 458 in subsequent deformation process like "buckets effect". Fig. 9 shows the contribution ratio of S1, 459 S2, and S3 at each point. The relationship of three ratios is stable generally except for several points 460 located on strong deformation belts. For X displacement, the dominate source is S3, and the 461 secondary and the third is S1 and S2, respectively. In the strong deformation belts, the source 462 contribution of S3 decreases sharply while the contribution of S1 increases at points, such as 8~3, 463 10~9, 14-7, 18~6, H-11, and Q~1 (Fig. 9a, 9b). The trend of Y displacement is similar to X 464 displacement. Transformation of dominate source is found in Z displacement (Fig. 9c). The 465 dominate sources S2 transfers to S1 or S3 in strong deformation belts. In addition, these changes 466 are the same trend in each monitoring campaigns, which means we can identify it at an early phase 467 of ground deformation helping to prevent risk like the severe damage of ventilation shaft on line 14 468 factor in a regional scale because we even not understand the mechanic of a single factor, much less 483 the coupling mechanics of all factors. However, it is more meaningful for safety management to 484 know the comprehensive effect of these influence factors. As a result, NMF successfully identify 485 three driving sources and the corresponding spatial and temporal patterns that is a reflection of the 486 comprehensive effect of these influence factors. It is worth noting that driving sources identification 487 is a core step to analyze the mechanism and evolution of ground deformation because they are 488 carriers of these influence factors. 489 From the analysis in section 4, we find that discontinuities such as joints and faults, may 490 dominant effects on the deformation behaviors of the rock mass in study area. Although NMF cannot 491 give the precise effect of discontinuities, a conclusion can be proposed that discontinuities 492 significantly contribute to rock movement in study area. In fact, field investigations find that there 493 are four groups of dominate joints (Fig. 10a, c) with an average of dip direction is 7.2° (percent: joints with E and S dip direction can significantly affect the ground deformation. The consistency 496 between direction of principal stress (Fig. 6a) and E and S dip direction favors the shear movement 497 of rock, therefore, causing the shear belt and extension belt on the ground (Fig. 6a). In addition, S1 498 and S3 have a high contribution of Z displacement (Fig. 4) while the pattern of S2 shows that S2 499 has a small contribution to X and Y displacement (Fig. 4), which suggests that the discontinuities 500 with high dip angle are highly developed and, therefore, the rock has a small slide component in X 501 and Y direction. Field investigations also reveal the highly developed of joints with high dip angle. 502 As shown in Fig. 10c, the average dip of four groups of joints is 59.6°, 54.2°, 64.2°, 63°, respectively. 503 The temporal pattern of three sources are also affected by discontinuity. The characteristic of 504 each stage in Fig. 8 is different. The amplitude of source continuously increases from the first stage 505 to the third stage. For example, S3 has a normalized amplitude of 0.63 in the first stage and increases 506 to 0.72 in the second stage, and 0.82 in the third stage (Fig. 11a). S1 and S2 have the similar changes 507 in amplitude (Fig. 11b, c). In addition, the repeated mining operations can disturb tectonic stress 508 causing the concentration or relaxation of stress but the self-weight body force remains stable. This 509 leads to the oscillation amplitude of S2 is smaller than that of S1 and S3. These temporal patterns 510 reveal that the overlaying rock mass is becoming increasingly sensitive to stress perturbations, 511 which is consistent with the process closing to the volumetric expansion point in uniaxial 512 compression test (Fig.12). During this process, the discontinuity becomes unstable and deforms 513 easily under small stress perturbation, especially shear stress perturbation. However, NMF cannot 514 predict when the ground deformation will evolve to failure (collapse of overlaying strata). 515

Conclusions 516
Our analyses demonstrate the applicability of a new inverse method for analysis of ground 517 deformation based on NMF algorithm of BSS. The solution is optimized by newly formulated 518 evaluation criteria to reduce the interference among sources and improve the resolution of each 519 pattern. The unknown sources are identified from a set of GPS time series data without any 520 information about the sources, underground mining conditions, and the physical processes 521 impacting the displacement propagation through rock mass. 522 The NMF based method identifies three unique sources causing the observed data. These are 523 listed in a descending order in terms of source contribution to original data: S1 (driven by horizontal 524 principal stress σ 2 ), S3(driven by horizontal principal stress σ 1 ), and S2 (driven by self-weight 525 body force caused by gravity). They appear to be proportionally manifested at the monitoring point. 526 Relative contribution of each source remains the same order as S1>S3>S2 except for the monitoring 527 points on strong deformation belts where the source contribution changes distinctly and even 528 dominate source transformation occurs. Besides, each source has stable spatial pattern which is 529 time-independent. This allows us to identify the potential risk area at an earlier stage. Finally, the 530 NMF based method can be applied to any real problem where temporal system behavior is 531 monitored at multiple locations. 532 Although NMF is useful to long-term and cyclical driving source, NMF cannot identify some 533 instant driving sources like blasting vibrating load and seismic load. The data sampling frequency 534 (interval: 0.5 year) affects the result of NMF and high frequency sampling may overcome this 535 problem and show more details of sources. 536

Acknowledgements 537
We gratefully acknowledge our financial support from the National Science

Declaration of Competing Interest 554
None. 555