Mining area deformation monitoring and risk assessment method

Abstract

The invention provides a mining area deformation monitoring and risk assessment method, which comprises the following steps of: carrying out registration, de-leveling, filtering and phase unwrapping on acquired image data through InSAR (Interferometric Synthetic Aperture Radar), generating a coherence graph, a filtering interference graph and an unwrapping phase graph, carrying out phase decomposition by using a singular value algorithm, separating a deformation phase, an atmospheric delay phase and noise, and carrying out risk assessment. The method comprises the steps of expanding discrete deformation points into a continuous deformation field through Kriging interpolation, carrying out deformation grading by taking a deformation rate, an accumulated settlement amount and a deformation gradient as core indexes, identifying a deformation high-risk area, arranging measurement robot monitoring points in the high-risk area, and constructing a risk assessment model in combination with real-time deformation data of a measurement robot and geological parameters. And evaluating the risk grade of each part of the mining area through the dynamic risk index model, and triggering graded early warning. According to the method, the InSAR and measurement robot technology is combined to solve the data fusion problem of large-range deformation monitoring and local high-precision monitoring, and the dynamic risk assessment precision under the complex terrain of the mining area is improved.

Description

Technical Field

[0001] The present invention belongs to the technical fields of surveying and mapping and remote sensing, and particularly relates to a method for monitoring mining area deformation and risk assessment. Background Art

[0002] Surface deformation monitoring in mining areas is a core link to ensure safe production in mines and prevent geological disasters. Traditional monitoring means mainly rely on leveling and the Global Navigation Satellite System (GNSS), obtaining deformation data through discrete points. Its spatial resolution is low, the labor cost is high, and it is difficult to achieve large-scale continuous coverage. Slope radar can monitor large-scale deformation of mine slopes in real time, but the observation surface is limited. Interferometric Synthetic Aperture Radar (InSAR) technology can obtain millimeter-level deformation information of the entire mining area by analyzing the phase difference of radar satellite images, and has the ability of all-weather and large-scale monitoring. However, InSAR still faces problems such as atmospheric delay error, decorrelation, and time resolution limitation in mining area applications, and it is difficult to capture sudden deformations or local dynamic changes. To make up for the deficiencies of InSAR, total station robots, as a kind of high-precision automated dynamic monitoring equipment, are widely used in key areas. By networking total station instruments and prisms, it can achieve real-time monitoring of millimeter-level three-dimensional deformation and is sensitive to local deformation gradients. However, the coverage of total station robots is limited, and the deployment density is positively correlated with the cost, and it cannot independently support the monitoring requirements of the entire mining area.

[0003] After retrieval, it is found that existing mining area monitoring methods have problems such as imperfect monitoring systems and missing risk assessment models. The process integration of InSAR area monitoring and total station robot point monitoring has not been streamlined, and key factors such as deformation acceleration and geological parameters have not been fully incorporated into the mining area risk assessment system. Summary of the Invention

[0004] Aiming at the above deficiencies of the prior art, the present invention provides a method for monitoring mining area deformation and risk assessment, which effectively solves the data fusion problem of large-scale deformation monitoring and local high-precision monitoring by combining InSAR and total station robot technologies, and effectively improves the accuracy of dynamic risk assessment under complex terrain in mining areas.

[0005] The technical solution provided by the present invention: A method for monitoring mining area deformation and risk assessment, comprising the following steps:

[0006] (1) Obtain ascending and descending orbit SAR satellite images covering the mining area, introduce high-precision DEM data for terrain phase correction, convert the DEM into a standard format and crop it to the monitoring area range;

[0007] (2) Adopt the short baseline set strategy, select SAR images to form interferometric pairs, generate differential interferogram sets, perform registration, detrending, filtering, and phase unwrapping on each pair of master and slave images to generate coherence maps, filtered interferograms, and unwrapped phase maps;

[0008] (3) Use singular value decomposition to jointly solve each interferometric subset, obtain the time series deformation rate, separate the deformation phase, atmospheric delay phase, and noise, correct the atmospheric delay error in combination with the atmospheric model, expand the discrete deformation points into a continuous deformation field through Kriging interpolation, and calculate the cumulative deformation amount through the integral of the product of the deformation rate and time;

[0009] (4) Use the deformation rate, cumulative settlement amount, and deformation gradient as the core indicators for deformation grading, set the deformation rate threshold, cumulative settlement threshold, and deformation gradient threshold, trigger the grading response for the InSAR deformation value, and extract high-risk areas;

[0010] (5) In the high-risk areas, set up monitoring points for total station robots, and configure total station robots outside the high-risk areas to achieve automatic prism tracking with millimeter-level accuracy and real-time acquisition of three-dimensional coordinates;

[0011] (6) Normalize the deformation monitoring data of the total station robot and the geological parameters of the monitoring points, map the deformation rate, cumulative deformation amount, and geological vulnerability coefficient to the 0-1 interval through the range method, and establish a comprehensive risk assessment index formula:

[0012] R = α·V norm + β·D norm + γ·G norm

[0013] Among them, V norm , D norm , G norm are the normalized values of the deformation rate, cumulative deformation amount, and geological vulnerability coefficient respectively, and α, β, γ are the weight coefficients of the risk assessment index, which are determined by grey relational analysis or expert experience method.

[0014] According to the value of the comprehensive risk index R, confirm the risk level and generate a warning signal.

[0015] Furthermore, in step (1), the Sentinel-1A data downloaded from the European Space Agency platform is used, and its time span covers the mining cycle and the geological activity sensitive period. The single-scene coverage range reaches dozens of square kilometers, and the data interval period is 6-12 days. The precise orbit file downloaded from the European Space Agency is used to improve the registration accuracy.

[0016] Furthermore, in step (3), the deformation phase φ obs is separated from the interferometric phase φ def , the topographic phase φ topo, Atmospheric phase φ atm and noise phase φ noise

[0017] φ obs = φ def + φ topo + φ atn + φ noise

[0018] Parametrize the above formula to establish a function model between the interferometric phase, the deformation rate, and the topographic error:

[0019]

[0020] where v is the linear deformation rate, Δt k represents the time baseline of the k-th interferometric pair, λ is the radar wavelength, is the vertical spatial baseline of the k-th interferometric pair, Δh is the topographic error, R is the radar slant range, and θ is the radar incidence angle.

[0021] Furthermore, in step (3), the model equation is solved by the least squares method:

[0022] B·x = d

[0023] where B is an M×N dimensional design matrix, M is the number of interferometric pairs, N is the number of time nodes, x = [v1, v2,..., v N-1 , Δh] T is the vector of parameters to be solved, and d is the observed phase vector.

[0024] Introduce singular value decomposition to process the non-full rank matrix and optimize the result with the minimum norm solution

[0025] B = U·∑·V T

[0026] where U and V are unitary matrices, and ∑ is a diagonal matrix of singular values. Calculate the deformation rate vector through the above formula. Finally, the deformation time series is calculated by the following formula, and the cumulative deformation amount is calculated by the integral of the product of the deformation rate and time:

[0027]

[0028] In the formula, Δt is the time interval between adjacent images, and V i represents the rate within each monitoring interval.

[0029] Furthermore, in step (5), the three-dimensional coordinates of the monitoring points are obtained in real time by the triangulation method and the polar coordinate method. The coordinate difference formula is:

[0030] ΔX = S·cosα·sinβ, ΔY = S·cosα·cosβ, ΔZ = S·sinα

[0031] where S is the slant range, α is the vertical angle, and β is the horizontal angle.

[0032] The measuring robot collects deformation data every 10 - 60 minutes, improving the monitoring frequency to the hourly level. The cumulative deformation D represents the difference between the most recent collected value and the first collected data value. The deformation rate is calculated using time series integration:

[0033]

[0034] where ΔD i is the deformation amount of a single monitoring, λ is the time decay coefficient, n is the number of collections, and t i is the collection time.

[0035] Furthermore, the formula for mapping the deformation rate, cumulative deformation amount, and geological condition parameters to the 0 - 1 interval in step (6) is as follows:

[0036]

[0037] where X min and X max are the extreme values of the corresponding parameters, and the corresponding parameters include the deformation rate, cumulative deformation amount, and various geological parameters. The geological vulnerability coefficient includes the rock layer compressive strength, fracture zone density, groundwater activity intensity, and terrain slope. Its calculation formula is as follows:

[0038] G norm = aI norm岩性 + bI norm断裂 + cI norm水文 + dI norm坡度

[0039] where I norm is the normalized value of each geological parameter index, and a, b, c, and d are the weights of each geological parameter index, respectively, obtained by the analytic hierarchy process.

[0040] Furthermore, in step (6), according to the above - mentioned comprehensive risk index R value, four risk levels are divided for four - level early warning, and different countermeasures are taken for different risk levels.

[0041]

[0042] The present invention has the following beneficial effects:

[0043] (1) InSAR provides a wide-area periodic deformation trend, while the total station focuses on dynamic tracking of high-risk areas. By integrating the data of the two, the advantages of large-scale monitoring are retained, and the accuracy of local deformation details is improved, thus solving the limitations of a single technology.

[0044] (2) InSAR provides a long-term deformation trend (such as annual-scale cumulative settlement), while the total station captures short-period dynamic events (such as blasting vibration, rainfall-induced displacement). The two are combined to construct a "trend + event" dual-driven risk assessment model. After identifying the overall settlement trend of the mining area through InSAR, the total station monitors the real-time local crack expansion rate, and at the same time combines geological parameters (such as rock layer strength, slope, fracture zone density, and groundwater content) to construct a risk assessment model, triggering hierarchical warnings (such as yellow / red alerts) to achieve spatio-temporal dynamic quantification of disaster risks. Description of the Drawings

[0045] Figure 1 is the technical flow chart of the present invention;

[0046] Figure 2 is the mining area risk assessment plan diagram of the present invention. Detailed Embodiments

[0047] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts shall fall within the protection scope of the present invention.

[0048] As Figure 1 shown, a method for monitoring deformation and risk assessment in a mining area specifically includes the following steps:

[0049] Step 1: InSAR deformation monitoring and disaster census

[0050] Download Sentinel-1A data (SLC format, ≥20 scenes) using the European Space Agency platform. Attention should be paid to the time baseline coverage range and polarization mode (such as VV / VH). Download the precise orbit file (POD) from the European Space Agency to improve the registration accuracy. Download SRTM3V4 (90m resolution) or higher-precision SRTM1 V3 (30m resolution), and the DEM needs to be converted to the standard format and cropped to the monitoring area range.

[0051] The constraint space baseline is less than 100 meters and the time baseline is less than 50 days to maintain the interference coherence. For each pair of master and slave images, registration, detrending, filtering, and phase unwrapping are performed to generate a coherence map, a filtered interferogram, and an unwrapped phase map. Interferogram pairs are manually screened to eliminate pairs with low coherence and retain the results with good phase continuity.

[0052] The deformation phase (φ obs ), topographic phase (φ def ), atmospheric phase (φ topo ), and noise phase (φ atm ) are separated from the interferometric phase (φ noise ) through phase decomposition.

[0053] φ obs = φ def + φ topo + φ atm + φ noise (1)

[0054] The above formula is parameterized to establish a functional model between the interferometric phase and the deformation rate and topographic error as shown in formula (2), and the linear deformation rate is obtained.

[0055]

[0056] Among them, v is the linear deformation rate (mm / year), Δt k represents the time baseline of the kth interferogram pair, λ is the radar wavelength, is the vertical space baseline of the kth interferogram pair, Δh is the topographic error (m), R is the radar slant range, and θ is the radar incidence angle.

[0057] The model equation is solved by the least squares method:

[0058] B·x = d (3)

[0059] Among them, B is an M×N-dimensional design matrix (M is the number of interferogram pairs, N is the number of time nodes), x = [v1, v2,..., v N-1 , Δh] T is the parameter vector to be solved, and d is the observed phase vector.

[0060] The singular value decomposition (SVD) is introduced to process the non-full rank matrix, and the result is optimized with the minimum norm solution

[0061] B = U·∑·V T (4)

[0062] Among them, U and V are unitary matrices, and ∑ is a diagonal matrix of singular values. The deformation rate vector can be calculated through the above formula, and the final deformation time series is calculated by the following formula, and the cumulative deformation amount is calculated by the integral of the product of the deformation rate and time:

[0063]

[0064] Where Δt is the time interval between adjacent images, and V i represents the deformation rate within each monitoring interval and is applicable to the monitoring of slow deformation at the edge of the mining area.

[0065] By comparing with the on-site measured data of the total station, the monitoring error of SBAS needs to be less than 6 mm to meet the demand for millimeter-level deformation monitoring in the mining area. Set the deformation rate threshold (30 mm / month), cumulative settlement threshold (100 mm), and deformation gradient threshold (5 mm / m) to trigger a graded response to the InSAR deformation value, extract high-risk areas, and focus on other deformation areas.

[0066] Step 2: Total station real-time dynamic local deformation monitoring

[0067] Use the InSAR census results (deformation rate field, cumulative deformation, or deformation gradient) to delineate high-risk areas, and the total station is densely arranged according to the grid (spacing less than 50 m). The deformation gradient = change amount / distance, and the deformation gradient is calculated in the direction outward from the settlement center.

[0068] The three-dimensional coordinates of the monitoring points are obtained in real time by using the angle-side measurement method and polar coordinate method. The coordinate difference formula is:

[0069] ΔX = S·cosα·sinβ, ΔY = S·cosα·cosβ, ΔZ = S·sinα (6)

[0070] where S is the inclined distance, α is the vertical angle, and β is the horizontal angle.

[0071] The total station collects deformation data every 10 - 60 minutes, and the monitoring frequency is increased to the hourly level. The cumulative deformation D represents the difference between the most recent collected value and the first collected data value. The deformation rate is calculated using time series integration:

[0072]

[0073] where ΔD i is the deformation amount of a single monitoring, λ is the time decay coefficient (default is 0.05 / h), n is the number of collections, and t i is the collection time.

[0074] Step 3: Risk assessment model construction and risk level assessment

[0075] Normalize the deformation monitoring data of the total station and the geological parameters of the monitoring points, and map parameters such as deformation rate, cumulative deformation amount, and geological conditions to the 0 - 1 interval through the range method:

[0076]

[0077] Among them, X min and X max are the extreme values of the corresponding parameters.

[0078] On this basis, a comprehensive risk assessment index formula is established:

[0079] R = α·V norm + β·D norm + γ·G norm (9)

[0080] Among them, V norm , D norm , G norm are the normalized values of the deformation rate (mm / t), the cumulative deformation (mm), and the geological vulnerability coefficient (dimensionless), respectively, which are mapped to the interval of 0-1 through the range method. α, β, and γ are the weight coefficients of the risk assessment index, which are determined by grey relational analysis or expert experience method. For example, they are generally set as α = 0.4, β = 0.3, γ = 0.3, and can be dynamically adjusted later. The geological vulnerability coefficient includes parameters such as the rock compressive strength, fracture zone density, groundwater activity intensity, and terrain slope. The formula is

[0081] G norm = aI norm岩性 + bI norm断裂 + cI norm水文 + dI norm坡度 (10)

[0082] Among them, I norm is the normalized value of each geological parameter index. The weights of each geological vulnerability coefficient need to calculate the weights of each factor in combination with the analytic hierarchy process. Generally, they can be set as 0.3, 0.2, 0.2, 0.3, and the weights are optimized in combination with the actual conditions of the mining area (such as the rock burst tendency and goaf distribution of metal mines).

[0083] According to the above-mentioned comprehensive risk index γ value, four risk levels are divided, and four-level early warnings are carried out, and different countermeasures are taken for different risk levels.

[0084]

[0085]

[0086] A deformation monitoring and risk assessment system for a mining area according to the present invention includes the following parts: Part 1: InSAR monitoring module

[0087] (1) Data collection and preprocessing

[0088] Obtain ascending and descending orbit SAR satellite images covering the mining area (taking Sentinel-1A as an example). The time span should cover the mining cycle and the sensitive period of geological activities. The single-scene coverage area reaches dozens of square kilometers, and the data interval period is 6 - 12 days. Introduce high-precision DEM data (such as SRTM or generated by UAV aerial survey) for terrain phase correction to eliminate the influence of terrain undulation on the interference phase.

[0089] (2) Generation of interferometric pairs and time series analysis

[0090] Adopt the Small Baseline Subset (SBAS) strategy. Select SAR images with a spatial baseline less than 100 meters and a temporal baseline less than 50 days to form interferometric pairs, generate a differential interferogram atlas, and reduce the influence of spatio-temporal decorrelation. Through multi-master image registration and orbit refinement, eliminate the phase shift caused by orbit errors, and use singular value decomposition to jointly solve each interferometric subset to obtain the time series deformation rate (with an accuracy reaching the millimeter level).

[0091] (3) Deformation information extraction and accuracy verification

[0092] Use the singular value algorithm for phase decomposition to separate the deformation phase, atmospheric delay phase, and noise, and combine the atmospheric model to correct the atmospheric delay error. Expand the discrete deformation points into a continuous deformation field through Kriging interpolation, and compare with the measured data of the total station to verify the reliability of the SBAS monitoring results.

[0093] (4) Deformation grading and hidden danger identification

[0094] Use the deformation rate, cumulative settlement amount, and deformation gradient as the core indicators for deformation grading. Mark the areas with a monthly settlement amount greater than 30 mm, a cumulative settlement amount greater than 100 mm (usually the center of the settlement funnel), or a deformation gradient greater than 5 mm / m (such as the edge of the settlement funnel) as high-risk areas.

[0095] Content 2: Total station monitoring module

[0096] Total station target layout

[0097] In the high-risk areas identified by SBAS, layout total station monitoring points with a spacing less than 50 meters, focusing on covering fault zones, settlement centers, edges of goafs, and landslide boundaries. Configure total stations (such as Leica TM60) outside the high-risk areas to achieve automatic prism tracking with millimeter-level accuracy (±0.5 mm) and real-time three-dimensional coordinate acquisition, and increase the monitoring frequency to the hourly level. For fault zones, edges of goafs, and landslide boundaries, it is necessary to densely arrange monitoring points on both sides.

[0098] Content 3: Mining area risk assessment module

[0099] Construction of dynamic risk assessment model

[0100] Construct a risk assessment model by combining the real-time deformation data (deformation magnitude, deformation trend, etc.) of the total station with geological parameters (rock layer strength, slope, fracture zone density, groundwater content, etc.), and evaluate the risk levels of various parts of the mining area through a dynamic risk index model (deformation rate + cumulative settlement + geological vulnerability coefficient) to trigger hierarchical early warnings.

[0101] Although the embodiments of the present invention have been shown and described, those of ordinary skill in the art can understand that various changes, modifications, substitutions, and variations can be made to these embodiments without departing from the principles and spirit of the present invention. The scope of the present invention is defined by the appended claims and their equivalents.

Claims

1. A method for monitoring deformation and risk assessment in a mining area, characterized in that, It includes the following steps: (1) Obtain ascending and descending orbit SAR satellite images covering the mining area, introduce high-precision DEM data for terrain phase correction, convert the DEM to a standard format and crop it to the monitoring area range; (2) Adopt the short baseline set strategy, select SAR images to form interferometric pairs, generate differential interferogram sets, perform registration, detrending, filtering, and phase unwrapping on each pair of master and slave images to generate coherence maps, filtered interferograms, and unwrapped phase maps; (3) Use singular value decomposition to jointly solve each interferometric subset, obtain the time series deformation rate, separate the deformation phase, atmospheric delay phase, and noise, correct the atmospheric delay error in combination with the atmospheric model, expand the discrete deformation points into a continuous deformation field through Kriging interpolation, and calculate the cumulative deformation amount through the integral of the product of the deformation rate and time; (4) Take the deformation rate, cumulative settlement amount, and deformation gradient as the core indicators for deformation grading, set deformation rate thresholds, cumulative settlement thresholds, and deformation gradient thresholds, trigger a grading response to the InSAR deformation value, and extract high-risk areas; (5) In the high-risk area, set up monitoring points for total station robots, and configure total station robots outside the high-risk area to achieve automatic prism tracking with millimeter-level accuracy and real-time acquisition of three-dimensional coordinates; (6) Normalize the deformation monitoring data of the total station robot and the geological parameters of the monitoring points, map the deformation rate, cumulative deformation amount, and geological vulnerability coefficient to the 0-1 interval through the range method, and establish a comprehensive risk assessment index formula: R = α·V norm + β·D norm + γ·G norm Among them, V norm , D norm , G norm are the normalized values of the deformation rate, the cumulative deformation amount, and the geological vulnerability coefficient, respectively. α, β, and γ are the weight coefficients of the risk assessment index, which are determined by grey relational analysis or the expert experience method. According to the comprehensive risk assessment index R value, confirm the risk level and generate a warning signal.

2. The method for monitoring deformation and risk assessment in a mining area according to claim 1, wherein In step (1), Sentinel-1A data downloaded from the European Space Agency platform is used. Its time span covers the mining cycle and the sensitive period of geological activities. The single-scene coverage range reaches dozens of square kilometers, and the data interval period is 6-12 days. The precise orbit file downloaded from the European Space Agency is used to improve the registration accuracy.

3. A method for monitoring deformation and risk assessment in a mining area according to claim 1, characterized in that, In the step (3), the deformation phase φ obs is separated from the interference phase φ def , the topographic phase φ topo , the atmospheric phase φ atm and the noise phase φ noise φ obs = φ def + φ topo + φ atm + φ noise Parameterize the above formula, establish a function model between the interferometric phase, deformation rate, and terrain error, and obtain the linear deformation rate. where v is the linear deformation rate, and Δt k represents the time baseline of the k-th interferometric pair, λ is the radar wavelength, is the vertical spatial baseline of the k-th interferometric pair, Δh is the topographic error, R is the radar slant range, and θ is the radar incidence angle.

4. A method for monitoring deformation and risk assessment in a mining area according to claim 1, characterized in that, In step (3), the model equation is solved by the least squares method: B·x=d Among them, B is an M×N dimensional design matrix, M is the number of interference pairs, N is the number of time nodes, x = [v1, v2,..., v N-1 , Δh] T is the parameter vector to be solved, d is the observed phase vector, and singular value decomposition is introduced to process the non-full rank matrix to optimize the result with the minimum norm solution B = U·∑·V T where U and V are unitary matrices, ∑ is a diagonal matrix of singular values. The deformation rate vector is calculated through the above formula. The final deformation time series is calculated by the following formula, and the cumulative deformation amount is calculated through the integral of the product of the deformation rate and time: where Δt is the time interval between adjacent images, and V i represents the rate within each monitoring interval.

5. A method for monitoring deformation and risk assessment in a mining area according to claim 1, characterized in that, In step (5), the three-dimensional coordinates of the monitoring points are obtained in real time by the triangulation method and the polar coordinate method. The coordinate difference formula is: ΔX=S·cosα·sinβ, ΔY=S·cosα·cosβ, ΔZ=S·sinα where S is the slant range, α is the vertical angle, and β is the horizontal angle. The total station robot collects deformation data every 10-60 minutes, improves the monitoring frequency to the hourly level. The cumulative deformation D represents the difference between the most recent collected value and the first collected data value. The deformation rate is calculated using time series integration. Among them, ΔD i is the deformation amount measured once, λ is the time decay coefficient, n is the number of acquisitions, and t i is the acquisition time.

6. The deformation monitoring and risk assessment method for a mining area according to claim 1, characterized in that The formulas for mapping the deformation rate, cumulative deformation amount, and geological conditions to the 0-1 interval in step (6) are as follows: Among them, X min and X max are the extreme values of the corresponding parameters. The corresponding parameters include the deformation rate, the cumulative deformation amount, and various geological parameters. The geological vulnerability coefficient includes the rock compressive strength, the fracture zone density, the groundwater activity intensity, and the terrain slope. Its calculation formula is as follows: G norm = aI norm岩性 + bI norm断裂 + cI nomr水文 + dI norm坡度 Among them, I norm is the normalized value of each geological parameter index, and a, b, c, and d are the weights of each geological parameter index, respectively, which are calculated by the analytic hierarchy process.

7. A method for monitoring deformation and risk assessment in a mining area according to claim 1, characterized in that, In step (6), according to the comprehensive risk assessment index R value, four risk levels are divided for four-level early warning, and different countermeasures are taken for different risk levels.

Images