Disclosure of Invention
Aiming at the defects of the prior art, the underground mining subsidence monitoring method and system integrating time sequences and space-time deformation overcome the problem that interference phase aliasing is caused by overlarge subsidence gradient of underground mining and disengagement is caused by non-coherent subsidence centers, and other types of disentanglement problems (such as time correlation) are processed through an SBAS (system based architecture), a spatial interpolation method and a time model, so that the time sequence of the full subsidence surface of a mining area is accurately obtained.
In order to achieve the above object, the present invention provides the following solutions:
The underground mining subsidence monitoring method integrating the time sequence and the space-time deformation comprises the following steps:
Generating an InSAR interferogram pile based on the existing mining area SAR image;
Based on the InSAR interferogram pile, obtaining a mining area prediction settlement model and an unwrapped interferogram pile by an MA-PU auxiliary model phase unwrapping method;
Performing linear error removal and reference point removal processing on the unwrapped interferogram stack;
the pixel classification is carried out on the disentangled interferogram pile after the processing, and the disentangled interferogram pile is filled from the spatial relation and the time relation of pixel points based on pixel types;
Singular value decomposition is carried out on the filled disentangled interferogram pile, and a time deformation sequence of the sedimentation surface is obtained;
and based on the time deformation sequence, completing underground mining subsidence monitoring of fusion time sequence and space-time deformation.
Preferably, the method for obtaining the InSAR interferogram pile comprises the following steps:
performing image pairing according to the time base lines and the space base lines of all the mining area SAR images to obtain paired images;
oversampling is carried out on the paired images in the ranging direction, so as to obtain a sampled image;
and generating the InSAR interferogram stack by a differential interference technology based on the sampling images, wherein one pair of images generates one interferogram.
Preferably, the method for obtaining the mining area prediction settlement model by using the MA-PU auxiliary model phase unwrapping method comprises the following steps:
Based on a branch cutting method and a traditional MCF method, carrying out phase unwrapping on the interferograms in the InSAR interferogram stack to obtain an initial unwrapping phase;
Sub-sampling the initial phase unwrapping result based on a quadtree algorithm, and inverting the space model to obtain space model parameters;
Based on the space model parameters, obtaining a space model phase, namely a mining area prediction sedimentation model;
calculating a residual error phase between the spatial model phase and the initial unwrapping phase, and iterating;
And calculating the sum of the spatial model phases of all iterations and adding the sum with the unwrapped residual phase of the last iteration to obtain a final deformation, wherein the final deformation is extracted from the unwrapped interferogram.
Preferably, the pixel classification result comprises a full rank point, a far field point, a middle point and a near field point, and singular value decomposition is carried out on various pixels to obtain a time deformation sequence of a settlement surface;
And filling blank pixels of the intermediate point and the near field point in the unwrapped interferogram pile by using the space model and the time model, wherein the time model adopts a Logistic model.
Preferably, the method for filling the blank pixels of the intermediate points in the unwrapped interferogram stack comprises the following steps:
Inversion is carried out to obtain time model parameters of the time model based on a genetic algorithm and accumulated deformation of the exploitation date;
and calculating the pixel value of a blank pixel of the intermediate point in each unwrapping interferogram in the unwrapping interferogram stack based on the time model parameters.
The invention also provides an underground mining subsidence monitoring system integrating the time sequence and the space-time deformation, which is used for realizing the underground mining subsidence monitoring method and comprises the following steps:
The map pile acquisition module is used for generating an InSAR interferogram pile based on the existing mining area SAR image;
The unwrapping module is used for obtaining a mining area prediction settlement model and an unwrapping interferogram pile through an MA-PU auxiliary model phase unwrapping method based on the InSAR interferogram pile;
The unwrapping interferogram stack processing module is used for carrying out the processing of the unwrapping interferogram stack on the nonlinear error and the reference point;
The time deformation sequence acquisition module is used for classifying the pixels of the processed disentangled interferogram pile, filling the disentangled interferogram pile from the spatial relationship and the time relationship of pixel points based on pixel types, and carrying out singular value decomposition on the filled disentangled interferogram pile to acquire the time deformation sequence of the settlement surface;
and the monitoring module is used for completing underground mining settlement monitoring of fusion time sequence and space-time deformation based on the time deformation sequence.
Preferably, the map pile obtaining module includes:
the pairing unit is used for carrying out image pairing according to the time base lines and the space base lines of all the mining area SAR images to obtain paired images;
The sampling unit is used for oversampling the paired images in the ranging direction to obtain sampling images;
And the interferogram pile generation unit is used for generating the InSAR interferogram pile based on the sampling images through a differential interference technology, wherein one pair of images generate one interferogram.
Preferably, the unwrapping module includes:
The initial phase unwrapping unit is used for carrying out phase unwrapping on the interferograms in the InSAR interferogram stack based on a branch cutting method and a traditional MCF method to obtain an initial unwrapping phase;
The space model inversion unit is used for sub-sampling the initial phase unwrapping result based on a quadtree algorithm, and inverting the space model to obtain space model parameters;
The spatial model phase acquisition unit is used for acquiring a spatial model phase based on the spatial model parameters, wherein the spatial model phase is the mining area prediction sedimentation model;
The iteration unit is used for calculating the residual error phase between the space model phase and the initial unwrapping phase, and carrying out iteration;
And the model construction unit is used for calculating the sum of the spatial model phases of all iterations and adding the sum with the unwrapped residual phase of the last iteration to obtain a final deformation, wherein the final deformation is extracted from the unwrapped interferogram.
Compared with the prior art, the method has the beneficial effects that the time deformation sequence of the full settlement surface is obtained under the influence of factors such as large gradient deformation, time incoherence and the like based on the subsidence model of underground mining deformation, the spatial relationship and the time characteristics, so that more comprehensive deformation data is provided for mining area monitoring.
Detailed Description
The following description of the embodiments of the present invention will be made clearly and completely with reference to the accompanying drawings, in which it is apparent that the embodiments described are only some embodiments of the present invention, but not all embodiments. All other embodiments, which can be made by those skilled in the art based on the embodiments of the invention without making any inventive effort, are intended to be within the scope of the invention.
In order that the above-recited objects, features and advantages of the present invention will become more readily apparent, a more particular description of the invention will be rendered by reference to the appended drawings and appended detailed description.
Example 1
As shown in fig. 1, the underground mining subsidence monitoring method integrating time sequence and space-time deformation comprises the following steps:
S1, generating an InSAR interferogram pile based on an existing mining area SAR image;
A further embodiment is that the method for obtaining the InSAR interferogram stack comprises:
performing image pairing according to the time base lines and the space base lines of all the mining area SAR images to obtain paired images;
Oversampling is carried out on the paired images in the ranging direction, so as to obtain a sampled image;
and generating an InSAR interferogram stack based on the sampled images by a differential interference technology, wherein one pair of images generates an interferogram.
In this embodiment, as shown in part (a) of fig. 1, in order to obtain the time series information of the formation of the full-subsidence surface of a certain region, firstly, by publishing or purchasing the SAR images of the region in the corresponding time period, if the images are paired according to the time base lines and the space base lines of all the images, the paired images are oversampled in the ranging direction, and then, an interferogram pile is generated by the differential interference technology (one interferogram is generated by a pair of images).
S2, obtaining a mining area prediction subsidence model and an unwrapped interferogram pile by an MA-PU auxiliary model phase unwrapping method based on the InSAR interferogram pile, namely unwrapping each interferogram in the interferogram pile by adopting an MA-PU method, wherein the MA-PU method refers to the whole flow from the interferogram to the unwrapped image in the section (B) in fig. 1, and in the flow, the mining area prediction subsidence model (namely a space model subsidence part in the section C in fig. 1) can be obtained.
The further implementation mode is that the method for obtaining the mining area prediction settlement model through the MA-PU auxiliary model phase unwrapping method comprises the following steps:
(1) Based on a branch cutting method and a traditional MCF method (minimum cost flow unwrapping method), carrying out phase unwrapping on interferograms in an InSAR interferogram stack to obtain an initial unwrapping phase;
(2) Based on the quadtree algorithm, the initial phase unwrapping result is sub-sampled to reduce the number of far field measurement points. And inverting the space model to obtain space model parameters, inverting the space model parameters by using Genetic Algorithm (GA) based on the data obtained from the quadtree, and establishing the space model from the parameters and the formulas of the corresponding model
(3) Based on the space model parameters, obtaining a space model phase, wherein the space model phase is the mining area prediction sedimentation model;
(4) Calculating the residual phase between the spatial model phase and the initial unwrapping phase (subtracting the initial unwrapping phase from the model phase, i.e., the residual phase), and performing the iterations of (1) - (4), in this embodiment, by repeating (1) to (4) to iterate, taking the output of (4) as the input of (1), until no further improvement can be observed, i.e., no residual phase map with smaller phase value and smoother spatial distribution can be obtained.
And calculating the sum of the spatial model phases of all iterations and adding the sum with the unwrapped residual phase of the last iteration to obtain a final deformation, wherein the final deformation is extracted from the unwrapped interferogram. Specifically, the first iteration obtains a spatial model d1, the second iteration obtains d2, the final spatial model is d=d1+d2+, and the final deformation is d+the unwrapped residual phase of the last iteration.
Since the inversion of the deformation model is performed a plurality of times, the starting point of the parameter used for the inversion is not important, and a wide range of constraint boundaries can be used. Therefore, the influence of model imperfection and inaccurate geological parameters is greatly reduced.
The method was applied to Australia West CliffColliery, using 23 pieces of ALOS-1 data, and validated against GPS data. The results show that the method improves the accuracy and the correlation compared with the traditional method, and the measurement accuracy is improved from 30 mm and 0.97 correlation to 20 mm and 0.99 correlation respectively.
S3, carrying out linear error removal and reference point removal treatment on the unwrapped interferogram pile, removing linear trend through least square plane fitting in the embodiment, and removing reference points, namely selecting a region almost free of deformation in time, taking the region as a reference point, and subtracting the average phase of the reference point from the whole unwrapped interferogram. This is done for each map of the unwrapped map heap.
And S4, carrying out pixel classification on the processed disentangled interferogram pile, filling the disentangled interferogram pile from the spatial relationship and the time relationship of pixel points based on pixel types, and carrying out singular value decomposition on the filled disentangled interferogram pile to obtain a time deformation sequence of the settlement surface.
The invention aims to realize monitoring of the subsidence of the working surface of the whole mining area, and in practice, the generated interferogram has the problems that partial areas cannot be disentangled (i.e. deformation cannot be acquired) due to the large gradient subsidence and the destructive interference of the mining area, and the MA-PU method and pixel point classification are used for simulating the deformation value of the partial areas incapable of disentanglement according to the spatial relation (spatial model and spatial interpolation) and the time subsidence characteristic (time model) of the pixel points.
Specifically, time series sedimentation values are generated based on the mining area sedimentation point type. And classifying the space positions of the sedimentation points according to the full rank and the rank deficiency, and obtaining the time deformation sequence of the whole sedimentation surface.
The further implementation mode is that the pixel classification result comprises a full rank point, a far field point, an intermediate point and a near field point, singular value decomposition is carried out on various pixels to obtain a time deformation sequence of a settlement surface, and the specific flow design is shown in figure 1. Near field points refer to points located near the center of settlement, intermediate points refer to points located at the edges of settlement, and far field points refer to points away from the settlement region.
And filling blank pixels of intermediate points and near-field points in the unwrapped interferogram pile by using a space model and a time model, wherein the time model adopts a Logistic model.
In this embodiment, for each pixel, a design matrix is constructed from the successfully unwrapped interferograms for computing the deformation time series. Specifically, according to the interferogram successfully unwrapped, the time point is used for changing into a variable x, all equations ax=b are listed, the variable coefficient on the left side of the equations is the design matrix A, and the right side of the equations is the deformation b obtained by unwrapping each interference pair.
In the case of full rank pixels, the design matrix is full rank, and the deformation time sequence can be directly inverted. For far field points, the deformation time series is inverted using Singular Value Decomposition (SVD), similar to the SBAS algorithm. As described above, for ax=b, x=a -1 b is calculated, and since a is an m×n matrix, the non-invertible matrix needs to be pseudo-inverted by SVD (a -1), and for far field points, the distortion is small, and for the time zone distortion of the gap, the filling can be calculated by the adjacent time point distortion rate.
For intermediate and near field points, due to the nonlinear nature of mine subsidence, a spatial and temporal model is required to fill in blank pixels in the unwrapped interferograms before applying the time series inversion.
A further embodiment is a method for filling a blank pixel at an intermediate point in an unwrapped interferogram stack, comprising:
inverting to obtain time model parameters of the time model based on the genetic algorithm and accumulated deformation of the mining date;
based on the time model parameters, pixel values of blank pixels in each unwrapped interferogram at intermediate points within the unwrapped interferogram stack are calculated.
In this embodiment, for intermediate points, the deformation values missing in each individual unwrapped interferogram can be estimated from a time model, and a mining activity (e.g., a longwall mining face) can be obtained using the conventional form of a Logistic model:
Where D (t) is the displacement at t, D max represents the maximum displacement, k is the rate constant, c is the scale constant, and t 0 is the time at which the maximum sedimentation rate occurs. When a point is affected by multiple mining activities in time sequentially:
where D total (t) is the total displacement at time t and m represents the particular mth production event. The model allows the contribution of each longwall to the total displacement to be different, reflecting the uniqueness of the sedimentation at each site. Assuming InSAR paired data is available, the residual displacement of any observation point can be determined. The sum of squares of residuals of the X observation points is as follows:
Wherein, Representing InSAR derived satellite line of sight (LOS) deformation, t x is the acquisition time of the primary and secondary images. The main image is a reference image, and the sub image is a second image to be compared with the main image.
The inversion process of the present invention utilizes a Genetic Algorithm (GA) to obtain model parameters for a combined Logistic model. It should be noted that the Logistic model is not suitable for handling situations with lifting phenomena. To avoid incorrect inversion due to this problem, an accumulated deformation for each acquisition date is added to the inversion process, calculated by summing all consecutive pairs of data. And then, according to the inverted Logistic model parameters, calculating the pixel values of the vacant pixels of the intermediate points in the pile in each unwrapped interferogram, obtaining the deformation value of the intermediate points with full rank, and inverting the deformation time sequence. Specifically, according to the Logistic parameter and the formula thereof, deformation prediction values of all time points of the point can be obtained, if the deformation value of the point is missing (the deformation value of the point corresponding to a certain pair of unwrapping pictures is missing) in a time period of 368-414 days, the deformation value in the missing time period is obtained by subtracting the 368 th day deformation from the 414 th day deformation obtained by Logistic.
For near field points, since there is no data corresponding to the deformation peak period, relying on a time-based model alone may not accurately derive the correct model parameters. Therefore, spatial interpolation is needed to fill in these temporally empty data. In this study, the residual phase of each unwrapped interferogram is calculated first according to the corresponding spatial model in the spatial model (model of Okada or PIM, mogi, etc.) (the residual phase is obtained by subtracting the simulated deformation from the interferogram). And (3) applying an image filling technology to the unwrapped residual phase diagram, and then combining with the simulated deformation to generate a filled unwrapped interferogram. And updating the phase value of the near field point in the filling unwrapping interferogram to obtain a filling unwrapping interferogram stack. And then inverting the filled unwrapped image pile, and obtaining model parameters of the Logistic model from the filled unwrapped deformation data through a genetic algorithm. After obtaining the time model parameters, the deformation time series can be obtained by adopting the same method as the middle point (namely combining with successfully disentangled pixels).
And S5, completing underground mining subsidence monitoring of fusion time sequence and space-time deformation based on the time deformation sequence.
Example two
The invention also provides an underground mining subsidence monitoring system integrating the time sequence and the space-time deformation, which is used for realizing an underground mining subsidence monitoring method and comprises the following steps:
The map pile acquisition module is used for generating an InSAR interferogram pile based on the existing mining area SAR image;
The unwrapping module is used for obtaining a mining area prediction settlement model and an unwrapping interferogram pile through an MA-PU auxiliary model phase unwrapping method based on the InSAR interferogram pile;
The unwrapping interferogram stack processing module is used for carrying out the processing of the unwrapping interferogram stack on the nonlinear error and the reference point;
The time deformation sequence acquisition module is used for classifying the pixels of the processed disentangled interferogram pile, filling the disentangled interferogram pile from the spatial relationship and the time relationship of pixel points based on pixel types, and carrying out singular value decomposition on the filled disentangled interferogram pile to acquire the time deformation sequence of the settlement surface;
and the monitoring module is used for completing underground mining settlement monitoring of fusion time sequence and space-time deformation based on the time deformation sequence.
A further embodiment is that the pile acquisition module includes:
the pairing unit is used for carrying out image pairing according to the time base lines and the space base lines of all the mining area SAR images to obtain paired images;
The sampling unit is used for oversampling the paired images in the ranging direction to obtain sampling images;
And the interferogram pile generation unit is used for generating the InSAR interferogram pile based on the sampling images through a differential interference technology, wherein one pair of images generate one interferogram.
A further embodiment is that the unwrapping module comprises:
The initial phase unwrapping unit is used for carrying out phase unwrapping on the interferograms in the InSAR interferogram stack based on a branch cutting method and a traditional MCF method to obtain an initial unwrapping phase;
The space model inversion unit is used for sub-sampling the initial phase unwrapping result based on a quadtree algorithm, and inverting the space model to obtain space model parameters;
Wherein the quadtree algorithm divides the interferogram into four parts and calculates the Root Mean Square (RMS) in each quadrant, subdivides the quadrant further if RMS exceeds a given threshold, otherwise accepts and outputs the mean value of the quadrant. This process is performed in a loop until the size of the quadrant is equal to a given minimum root mean square threshold. The purpose is to downsample, reduce the follow-up space model inversion operation time.
The spatial model phase acquisition unit is used for acquiring a spatial model phase based on the spatial model parameters, wherein the spatial model phase is the mining area prediction sedimentation model;
The iteration unit is used for calculating the residual error phase between the space model phase and the initial unwrapping phase, and carrying out iteration;
The model construction unit is used for calculating the sum of the spatial model phases of all iterations and adding the sum with the unwrapped residual phase of the last iteration to obtain a final deformation quantity, specifically, the spatial model d1 is obtained in the first iteration, d2 is obtained in the second iteration, the final spatial model is d=d1+d2, and the final deformation quantity is d+the unwrapped residual phase of the last iteration. Wherein the final deformation is extracted from the unwrapped interferogram.
Example III
The area of investigation of this example is WEST CLIFF mining areas LW32, LW33, LW34 face
From fig. 2 ((a) the original interferogram after filtering, (b) the direct unwrapping of the image using MCF, (c) the sedimentation funnel image recovered by the method of the present invention, (d) the difference between (c) and (b)), it can be seen that the area pointed by the red arrow in sub-image (a) has a striped break due to noise, filtering and other factors. The maximum of the conventional direct use minimum cost flow method unwrapping is 25.8rad, while the maximum of the model iterative recovery method is 30.31rad. The difference between the two is concentrated in the broken range of the interference fringes, and the average difference is about 6.35rad, and the difference is obviously the broken range of the interference fringes. It is worth affirming that these factors can be effectively eliminated by the method, and correct results can be recovered. The method is improved in the accuracy of sedimentation monitoring, and has strong adaptability and robustness in practical application.
As shown in fig. 3, due to the influence of time decoherence and large gradient sedimentation in the mining area, the sedimentation surface information of a part of the area is partially lost in the mining process, and in other time periods, the area is not influenced by the large gradient sedimentation, and the sedimentation information of the time period is correctly acquired.
The experimental result shows that the time sequence deformation result obtained by the method is compared with GPS data, the average absolute error is 0.016m, the root mean square error is 0.02m, and the full settlement surface deformation time sequence diagram is shown in figure 4.
Under the condition that deformation data of a part of the subsidence area in the figure 3 are seriously lost, the method is applied to obtain the time sequence of the deformation of the full subsidence area of the mining area.
The above embodiments are merely illustrative of the preferred embodiments of the present invention, and the scope of the present invention is not limited thereto, but various modifications and improvements made by those skilled in the art to which the present invention pertains are made without departing from the spirit of the present invention, and all modifications and improvements fall within the scope of the present invention as defined in the appended claims.