A method for obtaining surface subsidence sequence by fusing spatiotemporal simulation phase and DInSAR phase

Abstract

The present invention discloses a method for obtaining a surface subsidence sequence by fusing a spatiotemporal simulation phase and a DInSAR phase, wherein the method comprises: obtaining a probability integral method model parameter and a time function model parameter of a mining monitoring area, and constructing a mining mining subsidence spatiotemporal model according to the probability integral method model parameter and the time function model parameter; obtaining a DInSAR phase of the mining monitoring area, and obtaining a simulation phase of the mining monitoring area according to the mining mining subsidence spatiotemporal model; obtaining a fused phase according to the DInSAR phase and the simulation phase, and unwrapping the fused phase to obtain a surface subsidence sequence of the mining monitoring area. Thus, the surface subsidence sequence of the mining monitoring area is obtained according to the fused phase, the accuracy of the phase unwrapping is improved, so as to improve the resolution of the surface subsidence sequence, thereby reducing the mining monitoring cost while improving the mining monitoring efficiency and the reliability of the monitoring results.

Description

Technical Field

[0001] The invention relates to the technical field of mining area surface subsidence monitoring, and in particular to a method for fusing a spatiotemporal simulation phase and a DInSAR phase to obtain a surface subsidence sequence. Background Art

[0002] At present, in the process of coal resource collection, it is usually necessary to obtain accurate surface subsidence information in mining areas to provide a basis for subsequent mining damage prediction, control and comprehensive management of subsidence areas. Among them, relevant technologies usually use leveling technology or DInSAR (differential interferometry of synthetic aperture radar) to obtain surface subsidence information in mining areas.

[0003] However, the problem with the related technology is that, on the one hand, in the process of using leveling technology to obtain surface subsidence information in mining areas, it is affected by the monitoring range of the mining area, the sparse distribution of leveling points, the long observation period, etc., which makes it difficult for the leveling technology to obtain large-area surface deformation information. On the other hand, in the process of using DInSAR technology to obtain surface subsidence information in mining areas, it is affected by the complex surface of the coal mining subsidence area, the dense surface vegetation, the fast surface subsidence speed, and the relatively small mining subsidence range, which leads to the low accuracy of DInSAR technology in mining surface subsidence monitoring. Therefore, there is an urgent need for a mining area surface subsidence information monitoring method that can obtain large-area surface deformation information with high accuracy. Summary of the invention

[0004] The present invention aims to solve at least one of the technical problems in the related art to a certain extent. To this end, the first purpose of the present invention is to propose a method for obtaining a surface subsidence sequence by fusing a spatiotemporal simulation phase and a DInSAR phase, which can obtain a surface subsidence sequence in a mining monitoring area according to the fused phase, improve the accuracy of phase unwrapping, and improve the resolution of the surface subsidence sequence, thereby reducing the mining monitoring cost while improving the mining monitoring efficiency and the reliability of the monitoring results.

[0005] In order to achieve the above-mentioned purpose, the method for fusing spatiotemporal simulation phase and DInSAR phase to obtain surface subsidence sequence proposed in an embodiment of the present invention includes: obtaining probability integration method model parameters and time function model parameters of a mining monitoring area, and constructing a spatiotemporal model of mining subsidence according to the probability integration method model parameters and the time function model parameters; obtaining the DInSAR phase of the mining monitoring area, and obtaining the simulation phase of the mining monitoring area according to the spatiotemporal model of mining subsidence; obtaining a fused phase according to the DInSAR phase and the simulation phase, and unwrapping the fused phase to obtain a surface subsidence sequence of the mining monitoring area.

[0006] According to the method for obtaining a surface subsidence sequence by fusing a spatiotemporal simulation phase and a DInSAR phase proposed in an embodiment of the present invention, the probability integral method model parameters and the time function model parameters of the mining monitoring area are obtained, and a spatiotemporal model of mining subsidence is constructed based on the probability integral method model parameters and the time function model parameters, and then the DInSAR phase of the mining monitoring area is obtained, and the simulation phase of the mining monitoring area is obtained based on the spatiotemporal model of mining subsidence, and the fused phase is obtained based on the DInSAR phase and the simulation phase, and the fused phase is unwrapped to obtain the surface subsidence sequence of the mining monitoring area. Thus, the surface subsidence sequence of the mining monitoring area is obtained based on the fused phase, the accuracy of the phase unwrapping is improved, so as to improve the resolution of the surface subsidence sequence, thereby reducing the mining monitoring cost while improving the mining monitoring efficiency and the reliability of the monitoring results.

[0007] In addition, the method for obtaining a surface subsidence sequence by fusing the spatiotemporal simulation phase and the DInSAR phase according to the above embodiment of the present invention may also have the following additional technical features:

[0008] According to one embodiment of the present invention, the method of obtaining the probability integration method model parameters of the mining monitoring area includes: obtaining level monitoring time series data and constructing an inversion objective function of the probability integration method model parameters; obtaining multiple groups of the probability integration method model parameters using multiple preset algorithms based on the level monitoring time series data and the inversion objective function; using a weighted combination prediction model to weightedly fuse the multiple groups of the probability integration method model parameters, and using the weighted fused combined prediction parameters as the probability integration method model parameters of the mining monitoring area.

[0009] According to one embodiment of the present invention, the method of obtaining the time function model parameters of the mining monitoring area includes: determining the monitoring point timing start time of the level monitoring time series data, and adjusting the level monitoring time series data according to the monitoring point timing start time; according to the monitoring point timing start time and the adjusted level monitoring time series data, using the iterative function extreme value LM algorithm to obtain the time function model parameters of the mining monitoring area.

[0010] According to one embodiment of the present invention, obtaining the DInSAR phase of the mining monitoring area includes: obtaining historical radar monitoring data of the mining monitoring area; generating an interference image pair corresponding to the historical radar monitoring data according to the DInSAR imaging time; and obtaining the DInSAR phase of the mining monitoring area according to the interference image pair.

[0011] According to one embodiment of the present invention, obtaining the DInSAR phase of the mining monitoring area according to the interference image pair includes: performing image cropping, precise orbit data loading, flat ground phase removal, terrain phase removal, interference, phase filtering, and geocoding on the interference image pair in sequence to obtain the DInSAR phase of the mining monitoring area.

[0012] According to one embodiment of the present invention, the simulated phase of the mining monitoring area is obtained according to the mining subsidence spatiotemporal model, including: calculating the surface vertical settlement corresponding to the DInSAR imaging time according to the mining subsidence spatiotemporal model; converting the surface vertical settlement into a slant range change in the radar line of sight; performing phase conversion on the slant range change to obtain the simulated phase of the mining monitoring area.

[0013] According to one embodiment of the present invention, obtaining the fused phase based on the DInSAR phase and the simulated phase includes: spatially resampling the DInSAR phase and the simulated phase; calculating the phase difference between the spatially resampled DInSAR phase and the simulated phase; filtering the phase difference using a Kalman filter algorithm to obtain a filtered phase; and summing the simulated phase and the filtered phase to obtain the fused phase.

[0014] According to an embodiment of the present invention, after acquiring the fused phase, the method further includes: correcting the fused phase to an interval of [-π,π].

[0015] According to an embodiment of the present invention, the unwrapping of the fused phase includes: unwrapping the fused phase using a fast two-dimensional phase unwrapping algorithm with non-continuous path reliability sorting.

[0016] According to an embodiment of the present invention, after unwrapping the fused phase, the method further comprises: performing Loess smoothing on the unwrapping result to obtain a smoothed surface subsidence sequence of the mining monitoring area of ​​the mining area.

[0017] Additional aspects and advantages of the present invention will be given in part in the following description and in part will be obvious from the following description, or will be learned through practice of the present invention. BRIEF DESCRIPTION OF THE DRAWINGS

[0018] Figure 1 It is a flow chart of a method for obtaining a surface subsidence sequence by fusing a spatiotemporal simulation phase and a DInSAR phase according to an embodiment of the present invention;

[0019] Figure 2 is a schematic diagram of a process for obtaining the probability integral method model parameters of a mining monitoring area according to an embodiment of the present invention;

[0020] Figure 3 It is a schematic diagram of a process for obtaining the probability integral method model parameters of a mining monitoring area according to a specific embodiment of the present invention;

[0021] Figure 4 is a schematic diagram of a process for obtaining time function model parameters of a mining monitoring area according to an embodiment of the present invention;

[0022] Figure 5 It is a schematic diagram of a process for obtaining time function model parameters of a mining monitoring area according to a specific embodiment of the present invention;

[0023] Figure 6 is a schematic diagram of a process for obtaining a DInSAR phase of a mining monitoring area according to an embodiment of the present invention;

[0024] Figure 7 It is a schematic diagram of a process of obtaining a simulated phase of a mining monitoring area according to a mining subsidence spatiotemporal model according to an embodiment of the present invention;

[0025] Figure 8 It is a schematic diagram of a process of obtaining a simulation phase of a mining monitoring area according to a mining subsidence spatiotemporal model according to a specific embodiment of the present invention;

[0026] Fig. 9 is a schematic diagram of a process of obtaining a fused phase according to a DInSAR phase and a simulated phase according to an embodiment of the present invention;

[0027] Fig.10 It is a schematic diagram of a process of obtaining a fused phase based on a DInSAR phase and a simulated phase according to a specific embodiment of the present invention. DETAILED DESCRIPTION

[0028] Embodiments of the present invention are described in detail below, examples of which are shown in the accompanying drawings, wherein the same or similar reference numerals throughout represent the same or similar elements or elements having the same or similar functions. The embodiments described below with reference to the accompanying drawings are exemplary and are intended to be used to explain the present invention, and should not be construed as limiting the present invention.

[0029] The following describes a method for acquiring a surface subsidence sequence by fusing a spatiotemporal simulation phase and a DInSAR phase according to an embodiment of the present invention with reference to the accompanying drawings.

[0030] Figure 1 It is a flow chart of a method for obtaining a surface subsidence sequence by fusing a spatiotemporal simulation phase and a DInSAR phase according to an embodiment of the present invention.

[0031] like Figure 1As shown, the method for fusing the spatiotemporal simulation phase and the DInSAR phase to obtain the surface subsidence sequence includes the following steps:

[0032] S101, obtaining probability integration method model parameters and time function model parameters of the mining monitoring area of ​​the mining area, and constructing a mining subsidence spatiotemporal model of the mining area according to the probability integration method model parameters and the time function model parameters.

[0033] Specifically, in some embodiments of the present invention, the probability integral method model parameter par of the mining monitoring area can be obtained by PIM , construct a probability integral method model (PIM model) of the corresponding mining monitoring area, and, by obtaining the time function model parameters of the mining monitoring area, construct a time function model (for example, Gompertz model) of the corresponding mining monitoring area, and then, construct a mining subsidence spatiotemporal model (for example, Gompertz-PIM model) of the corresponding mining monitoring area based on the probability integral method model and the time function model.

[0034] S102, obtaining the DInSAR phase of the mining monitoring area of ​​the mining area, and obtaining the simulation phase of the mining monitoring area of ​​the mining area according to the mining subsidence spatiotemporal model of the mining area.

[0035] Specifically, in some embodiments of the present invention, the DInSAR (Differential Interferometry Synthetic Aperture Radar) phase position of the mining monitoring area can be obtained based on radar monitoring data. In addition, the simulated phase of the mining monitoring area can be obtained based on the level monitoring time series data and the mining subsidence spatiotemporal model.

[0036] S103, obtaining a fused phase according to the DInSAR phase and the simulated phase, and unwrapping the fused phase to obtain a surface subsidence sequence in the mining monitoring area.

[0037] It should be understood that in the above-mentioned embodiment of the present invention, first, the probability integral method model parameters and the time function model parameters of the mining monitoring area are obtained, and the mining subsidence spatiotemporal model is constructed according to the probability integral method model parameters and the time function model parameters, and then the DInSAR phase of the mining monitoring area is obtained, and the simulation phase of the mining monitoring area is obtained according to the mining subsidence spatiotemporal model, and then the fusion phase is obtained according to the DInSAR phase and the simulation phase, and the fusion phase is unwrapped to obtain the surface subsidence sequence of the mining monitoring area. Thus, the surface subsidence sequence of the mining monitoring area is obtained according to the fusion phase, and the accuracy of the phase unwrapping is improved to improve the resolution of the surface subsidence sequence, thereby reducing the mining monitoring cost while improving the mining monitoring efficiency and the reliability of the monitoring results.

[0038] Furthermore, if Figure 2 As shown, the probability integration method model parameters of the mining monitoring area are obtained, including:

[0039] S201, obtaining level monitoring time series data and constructing an inversion objective function of probability integration method model parameters.

[0040] Specifically, in some embodiments of the present invention, level monitoring time series data can be obtained, and an inversion objective function of the probability integration method model parameters can be constructed, so that the probability integration method model parameters can be inverted using the inversion objective function and the level monitoring time series data.

[0041] S202, according to the level monitoring time series data and the inversion objective function, using a variety of preset algorithms, obtain multiple groups of probability integration method model parameters.

[0042] Specifically, in some embodiments of the present invention, multiple sets of probability integration method model parameters can be obtained based on the level monitoring time series data and the inversion objective function using a variety of preset algorithms (for example, genetic algorithm GA, particle swarm algorithm PSO, differential evolution algorithm DE, whale optimization algorithm WOA).

[0043] S203, using a weighted combination prediction model to weightedly fuse multiple groups of probability integration method model parameters, and using the weighted fused combination prediction parameters as the probability integration method model parameters of the mining monitoring area.

[0044] It should be noted that, due to the difference in parameter correlation and preset algorithms in the inversion process of the probability integration method model parameters, there are problems of different inversion results, that is, multiple sets of probability integration method model parameters are obtained. Therefore, in some embodiments of the present invention, a weighted combination prediction model (OWCM model) can be used to weightedly fuse multiple sets of probability integration method model parameters, and the weighted fused combined prediction parameters are used as the probability integration method model parameters of the mining monitoring area, thereby improving the accuracy and reliability of constructing the probability integration method model of the corresponding mining monitoring area.

[0045] The following is combined with Figure 3 In accordance with the specific embodiment of the present invention, a specific process of obtaining the probability integral method model parameters of the mining monitoring area is described accordingly.

[0046] Specifically, obtaining the probability integration method model parameters of the mining monitoring area includes the following steps:

[0047] S21, based on the PIM principle, calculate the settlement of any point (x, y) on the mining surface.

[0048] Specifically, the settlement of any point (x, y) on the mining surface can be obtained by the following calculation formula:

[0049]

[0050]

[0051] Among them, W0 is the maximum sinking value, W 0 (x) is the sinking value of the projection point on the main section, W 0 (y) The subsidence value of the projection point on the main dip section, u is the mining unit of the probability integral, m is the mining thickness, H0 is the average mining depth, H1 is the downhill mining depth, H2 is the uphill mining depth, q is the subsidence coefficient, tanβ is the tangent value of the main influencing angle, s4 is the offset distance of the inflection point of the left boundary, s3 is the offset distance of the inflection point of the right boundary, s2 is the offset distance of the uphill inflection point, s1 is the offset distance of the downhill inflection point, θ is the mining influence propagation angle, D1 is the length of the strike working face, and D3 is the length of the dip working face.

[0052] It should be noted that going uphill refers to the tunnel that is excavated above the mining level and along the inclination of the coal seam (that is, the tunnel that is inclined upward). Going downhill refers to the tunnel that is excavated below the mining level and along the inclination of the coal seam (that is, the tunnel that is inclined downward). The strike refers to the intersection of the rock layer and the horizontal plane. The directions indicated by the two ends of the strike line are the strike of the rock layer. There are countless parallel strike lines, but there are only two strikes, which differ by 180°. The dip refers to the direction of the projection of the dip line on the horizontal plane. The straight line on the layer that is perpendicular to the strike line and drawn downward along the inclination is called the dip line. The inflection point offset distance refers to the distance between the actual mining boundary and the calculated boundary.

[0053] Based on the above formula, it can be known that the probability integral method model parameters of the mining monitoring area that needs to be inverted may include: q, tanβ, θ, s4, s3, s2, s1.

[0054] Therefore, in a specific embodiment of the present invention, the probability integral method parameter inversion objective function f can be constructed vv , where f vv The expression is as follows:

[0055]

[0056] Where n is the number of monitoring points of the level monitoring time series data involved in the parameter inversion, W is the PIM calculation value, It is the time series data value of level monitoring.

[0057] S22, using a variety of preset algorithms (for example, genetic algorithm, particle swarm algorithm, differential evolution algorithm, whale optimization algorithm), looping iteratively to make fvv Reach the minimum.

[0058] Specifically, the actual values ​​of the probability integral method model parameters based on the mining monitoring area Due to the different parameter correlation and preset algorithm selection, taking the preset algorithms as the above genetic algorithm GA, particle swarm algorithm PSO, differential evolution algorithm DE, and whale optimization algorithm WOA as examples, the probability integration method model parameters of different mining monitoring areas can be obtained: (par1, par2, par3, par4),

[0059] Therefore, in a specific embodiment of the present invention, a weighted combination prediction model is also used to perform weighted fusion on multiple groups of probability integration method model parameters.

[0060] Specifically, the weights in the weighted combination prediction model are w = (w1, w2, w3, w4) T , the predicted value of the model parameter par of the probability integral method of the mining monitoring area using the jth preset algorithm j =[q j ,b j ,tanβ j ,θ j ,(s4,s3,s2,s1) j ], where j=1,2,3,4.

[0061] Furthermore, the fitting deviation can be calculated based on the actual value and predicted value of the probability integral method model parameter in the mining monitoring area of ​​the mining area, where the calculation formula is as follows:

[0062]

[0063] Among them, e t is the fitting deviation, is the actual value of the probability integration method model parameter in the mining monitoring area, par j It is the predicted value of the model parameters of the probability integration method in the mining monitoring area using the jth preset algorithm.

[0064] In addition, the sum of squares of the prediction errors of the weighted combination prediction model is:

[0065]

[0066]

[0067] The objective function expression is:

[0068]

[0069] R=(1,1,1…,1) T

[0070] The Lagrange multiplier algorithm is used to solve the objective function, and the Lagrange multiplier λ is introduced. Q can be expressed as:

[0071] Q=w T Ew+2λ(R T w-1)

[0072] Depend on have to:

[0073] R T w=1

[0074] R T λE -1 R=1

[0075]

[0076] Therefore, the optimal weight vector is:

[0077]

[0078] The minimum value of the objective function is:

[0079]

[0080] Construct a weighted combination prediction model:

[0081]

[0082] S23, using a weighted combination prediction model to perform weighted fusion on multiple groups of probability integration method model parameters, to obtain the probability integration method model parameters of the mining monitoring area after weighted fusion, and based on the probability integration method model parameters of the mining monitoring area, to construct a probability integration method model corresponding to the mining monitoring area.

[0083] It should be understood that in the above embodiments of the present invention, a variety of preset algorithms (for example, genetic algorithm GA, particle swarm algorithm PSO, differential evolution algorithm DE, whale optimization algorithm WOA) can be used to invert multiple groups of probability integration method model parameters of mining monitoring areas, and the weighted combination prediction model is used to weightedly fuse the multiple groups of probability integration method model parameters to obtain the weighted fusion probability integration method model parameters par of the mining monitoring area. PIM , so that the probability integral method model parameter par PIM Construct a probability integration method model corresponding to the mining monitoring area.

[0084] Furthermore, if Figure 4 As shown, the time function model parameters of the mining monitoring area are obtained, including:

[0085] S301, determining a monitoring point timing start time of the level monitoring time series data, and adjusting the level monitoring time series data according to the monitoring point timing start time.

[0086] Specifically, in some embodiments of the present invention, the time when the working face is mined to the point directly below the corresponding monitoring point can be determined as the monitoring point timing start time of the level monitoring time series data, and the level monitoring time series data can be adjusted according to the monitoring point timing start time, thereby ensuring the consistency of the level monitoring time series data of each monitoring point, which is conducive to improving the accuracy and reliability of constructing the time function model.

[0087] S302, according to the start time of the monitoring point sequence and the adjusted level monitoring sequence data, using the iterative function extreme value LM algorithm to obtain the time function model parameters of the mining monitoring area of ​​the mining area.

[0088] Specifically, in some embodiments of the present invention, an iterative function extreme value LM algorithm can be used to perform iterative calculations in combination with the monitoring point timing start time and the adjusted level monitoring timing data to obtain the time function model parameters of the mining monitoring area.

[0089] The following is combined with Figure 5 The specific process of obtaining the time function model parameters of the mining monitoring area of ​​the mining area is described accordingly with the specific embodiment of the present invention.

[0090] Specifically, obtaining the time function model parameters of the mining monitoring area of ​​the mining area includes the following steps:

[0091] S31, determining the time function of any monitoring point.

[0092] Among them, the Gompertz time function expression of any monitoring point is:

[0093]

[0094] Among them, W max is the maximum settlement, k is the growth coefficient, and b is the inflection point. The inflection point is specifically used to describe the dynamic deformation of the mining surface, which is affected by the starting settlement time and the settlement speed respectively.

[0095] Furthermore, in a specific embodiment of the present invention, the time when the working face is mined to the point directly below the corresponding monitoring point is determined as the monitoring point timing start time of the level monitoring timing data, and the level monitoring timing data is adjusted according to the monitoring point timing start time.

[0096] S32, taking the start time of the monitoring point sequence of the level monitoring time series data as the start time of the Gompertz time function G(t), and using the iterative function extreme value LM algorithm combined with the adjusted level monitoring time series data, the Gompertz time function is cyclically operated to obtain the time function model parameters (k, b) of the mining monitoring area, thereby constructing the time function model of the corresponding mining monitoring area according to the time function model parameters (k, b).

[0097] It should be understood that in the above-mentioned embodiments of the present invention, the level monitoring time series data can be used in combination with the iterative function extreme value (LM, Levenberg-Marquardt) algorithm to invert the time function model parameters (k, b) of the mining monitoring area of ​​the mining area, so that the time function model parameters (k, b) can be used to construct the time function model of the corresponding mining monitoring area.

[0098] Furthermore, in the above-mentioned probability integral method model parameter par PIM After constructing the probability integration method model W(x, y) of the corresponding mining monitoring area and the time function model G(t) of the corresponding mining monitoring area using the time function model parameters (k, b) as mentioned above, the specific process of constructing the mining subsidence spatiotemporal model of the corresponding mining monitoring area according to the probability integration method model and the time function model is explained.

[0099] Specifically, the average daily working face advancement distance [l1,l2,…,l n ], and divide the mining face into n blocks in units of days, where the size of n is determined by the mining speed of the working face. At this time, the settlement amount of each mining monitoring area can be obtained according to the probability integral method model [W1, W2, …, W n ], and the dynamic settlement of each mining monitoring area and the surface dynamic settlement can be obtained according to the time function model: W i (t) = G(t)W i , where W i (t) is the dynamic settlement of each mining monitoring area; Among them, W(t) is the dynamic surface settlement in the mining monitoring area of ​​each mining area.

[0100] Furthermore, if Figure 6 As shown, the DInSAR phase of the mining monitoring area is obtained, including:

[0101] S401, obtaining historical radar monitoring data of the mining monitoring area.

[0102] Specifically, in some embodiments of the present invention, the historical radar monitoring data of the mining monitoring area may be radar monitoring data obtained by monitoring the mining monitoring area within a preset time period using a C-band synthetic aperture radar satellite.

[0103] S402, generating an interferometric image pair corresponding to the historical radar monitoring data according to the DInSAR imaging time.

[0104] Specifically, in some embodiments of the present invention, the historical radar monitoring data may be differentially processed according to the DInSAR imaging time, thereby generating an interferometric image pair corresponding to the historical radar monitoring data.

[0105] S403, obtaining the DInSAR phase of the mining monitoring area of ​​the mining area according to the interference image pair.

[0106] Specifically, in some embodiments of the present invention, corresponding image processing may be performed on the interference image pair to obtain the DInSAR phase of the mining monitoring area.

[0107] The specific process of obtaining the DInSAR phase of the mining monitoring area is described below in conjunction with the specific embodiments of the present invention.

[0108] Specifically, InSAR can obtain phase information through the phase difference between two echo signals of the sensor and the parameters of the satellite orbit, where the phase information can be expressed as:

[0109]

[0110] in, is the InSAR phase information, is the flat ground phase, is the terrain phase, is the surface deformation phase, is the atmospheric phase caused by the ionospheric and tropospheric delays, and is the noise phase.

[0111] Furthermore, the round-trip slant distances correspond to the slant distance phases ψ1 and ψ2, which can be expressed as:

[0112]

[0113]

[0114] Among them, n1 and n2 are the phase integer cycles of the radar along the slant range R1 and R2, respectively. and are the non-integer phases of the radar's round-trip propagation along the slant ranges R1 and R2, and λ is the radar wavelength.

[0115] Furthermore, in a specific embodiment of the present invention, after removing the atmospheric phase Flat Earth Phase and noise phase After that, the surface deformation phase The expression is:

[0116]

[0117] Among them, R1-R2 is the phase change affected by terrain undulation, and Δr is the slant range change of DInSAR monitoring surface deformation along the radar line of sight. At this time, the terrain phase can be simulated by an external digital elevation model (DEM). and will If the term is removed, the surface deformation phase The expression can be further updated as:

[0118]

[0119] Furthermore, in a specific embodiment of the present invention, the historical radar monitoring data of the mining monitoring area can be used in combination with the above-mentioned surface deformation phase The interferometric image pairs are obtained by expression, and the interferometric image pairs are processed accordingly using SNAP software to obtain the DInSAR phase coordinates of the mining monitoring area.

[0120] It should be understood that in the above-mentioned embodiment of the present invention, the historical radar monitoring data of the mining monitoring area can be obtained from the historical radar monitoring data, and the interference image pair corresponding to the historical radar monitoring data is generated according to the DInSAR imaging time, so that the DInSAR phase of the mining monitoring area can be obtained according to the interference image pair.

[0121] Furthermore, the DInSAR phase of the mining monitoring area is obtained according to the interference image pair, including: image cropping, precise orbit data loading, flat ground phase removal, terrain phase removal, interference, phase filtering, and geocoding of the interference image pair in sequence to obtain the DInSAR phase of the mining monitoring area.

[0122] It should be understood that in the above-mentioned embodiments of the present invention, the SNAP software can be used to perform image cropping, precise orbit data loading, flat ground phase removal, terrain phase removal, interference, phase filtering, geocoding and other data processing on the interference image pairs in sequence, so as to obtain the DInSAR phase of the mining monitoring area.

[0123] Furthermore, if Figure 7 As shown, the simulated phase of the mining monitoring area is obtained according to the mining subsidence spatiotemporal model, including:

[0124] S501, calculating the surface vertical settlement amount corresponding to the DInSAR imaging time according to the mining subsidence spatiotemporal model of the mining area.

[0125] Specifically, in some embodiments of the present invention, the constructed mining subsidence spatiotemporal model can be used to calculate the vertical settlement of the ground surface corresponding to the DInSAR imaging time.

[0126] S502, converting the vertical settlement of the ground surface into a slant range change in the radar line of sight.

[0127] Specifically, in some embodiments of the present invention, the vertical settlement of the ground surface can be converted into the slant distance change in the radar line of sight, so as to obtain the simulated phase of the mining monitoring area according to the slant distance change in the radar line of sight.

[0128] S503, performing phase conversion on the slant distance variation to obtain a simulated phase of the mining monitoring area.

[0129] Specifically, in some embodiments of the present invention, the slant range variation of the radar line of sight can be phase-converted to obtain the simulated phase of the mining monitoring area.

[0130] The following is combined with Figure 8 In accordance with the specific embodiment of the present invention, the specific process of obtaining the simulated phase of the mining monitoring area is described accordingly.

[0131] S51, the surface vertical settlement corresponding to DInSAR is calculated by constructing the Gompertz-PIM mining subsidence spatiotemporal model.

[0132] S52, converting the vertical settlement of the ground surface into a change in the slant range from the radar line of sight to the LOS.

[0133] The calculation formula for the change in slant range from the radar line of sight to the LOS is as follows:

[0134] Δr(t)=W(t)cosθ;

[0135] Where θ is the radar incident angle, Δr(t) is the change in slant range from the radar line of sight to the LOS, and W(t) is the vertical settlement of the ground surface.

[0136] S53, converting the slant range variation of the radar line of sight into a simulated phase value of the mining monitoring area.

[0137] Among them, the calculation formula of the simulated phase value of the mining monitoring area is as follows:

[0138]

[0139] in, It is the simulated phase value of the mining monitoring area.

[0140] Based on the above formula, the simulated phase of the radar line of sight (i.e., the simulated phase of the mining monitoring area) can be generated from the dynamic surface settlement W(t). ), the calculation formula is as follows:

[0141]

[0142] It should be understood that in the above-mentioned embodiment of the present invention, the constructed mining subsidence spatiotemporal model is used to calculate the surface vertical settlement W(t) corresponding to the DInSAR imaging time, so that according to the phase generation principle, the simulated phase of the mining monitoring area is calculated by W(t).

[0143] Furthermore, if Fig. 9 As shown, the fusion phase is obtained according to the DInSAR phase and the simulated phase, including:

[0144] S601, spatially resample the DInSAR phase and the simulated phase.

[0145] Specifically, in some embodiments of the present invention, the DInSAR phase and the simulated phase may be spatially resampled to make the spatial resolutions of the DInSAR phase and the simulated phase consistent, so that the DInSAR phase and the simulated phase may be subsequently phase-fused.

[0146] S602, calculating the phase difference between the spatially resampled DInSAR phase and the simulated phase.

[0147] Specifically, in some embodiments of the present invention, since the spatial resolutions of the spatially resampled DInSAR phase and the simulated phase are consistent, the difference between the spatially resampled DInSAR phase and the simulated phase can be used as the phase difference between the DInSAR phase and the simulated phase.

[0148] S603, using a Kalman filter algorithm to filter the phase difference to obtain a filtered phase.

[0149] Specifically, in some embodiments of the present invention, a Kalman filter algorithm may be used to filter the phase difference between the spatially resampled DInSAR phase and the simulated phase, thereby obtaining a filtered phase.

[0150] S604, summing the simulated phase and the filtered phase to obtain a fused phase.

[0151] Optionally, in some other embodiments of the present invention, other fusion algorithms (for example, Helmert variance component estimation, neural network, etc.) may be used to fuse the DInSAR phase and the simulated phase to obtain a fused phase.

[0152] Specifically, in some embodiments of the present invention, after obtaining the filtered phase, the simulated phase and the filtered phase may be added to obtain the fused phase.

[0153] The following is combined with Fig.10 The specific process of obtaining the fusion phase is described accordingly with the specific embodiment of the present invention.

[0154] S61, resample the DInSAR phase space to the same spatial resolution as the simulated phase.

[0155] It can be understood that, in some embodiments of the present invention, the DInSAR phase space is resampled to be consistent with the spatial resolution of the simulated phase, thereby preparing for subsequent phase fusion.

[0156] S62, calculating the phase difference between the simulated phase and the DInSAR phase.

[0157] The calculation formula of the phase difference between the simulated phase and the DInSAR phase is as follows:

[0158]

[0159] in, To simulate the phase, is the DInSAR phase, is the phase difference.

[0160] S63, using a Kalman filter algorithm to filter the phase difference to obtain a filtered phase.

[0161] Specifically, in a specific embodiment of the present invention, the Kalman filter algorithm is used to The 1 / 3 standard deviation of the observed noise variance R for all pixels of the DInSAR phase With analog phase Perform filtering and perform the following calculations:

[0162] (1) Calculate the state forecast value at time t+1:

[0163] (2) Calculate the covariance matrix of the forecast value: P(t+1|t)=φP(t|t)φ T +ΓQΓ T ;

[0164] (3) Calculate the gain matrix at time t+1: K(t+1)=P(t+1|t)H T [HP(t+1|t)H T +R] -1 ;

[0165] (4) Calculate the state estimate at time t+1:

[0166] (5) Calculate the state estimate to obtain the covariance matrix: P(t+1|t+1)=[I n -K(t+1)H]P(t+1|t);

[0167] Thus, through the above loop operation, the filtering phase is obtained

[0168] S64, adding the simulated phase and the filtered phase to obtain a fused phase.

[0169] The calculation formula of the fusion phase is as follows:

[0170]

[0171] in, is the filtering phase, The fusion phase.

[0172] It should be understood that in the above embodiments of the present invention, the Kalman filter algorithm can be used to fuse the simulated phase Phase with DInSAR Get fusion phase Therefore, the surface subsidence sequence in the mining monitoring area is obtained according to the fused phase, and the accuracy of phase unwrapping is improved to improve the resolution of the surface subsidence sequence, thereby reducing the mining monitoring cost while improving the mining monitoring efficiency and the reliability of monitoring results.

[0173] Furthermore, after acquiring the fused phase, the method further includes: correcting the fused phase to the interval [-π,π].

[0174] It should be understood that in the above-mentioned embodiment of the present invention, the fusion phase Corrected to the interval [-π,π], where the calculation formula is as follows:

[0175]

[0176] Furthermore, the fused phase is unwrapped, including: unwrapping the fused phase using a fast two-dimensional phase unwrapping algorithm based on non-continuous path reliability sorting.

[0177] Optionally, in some other embodiments of the present invention, the fused phase may be unwrapped using an LP minimum norm method, a minimum cost flow method, or the like.

[0178] It should be understood that in the above-mentioned embodiment of the present invention, a fast two-dimensional phase unwrapping algorithm of non-continuous path reliability sorting can be used to unwrap the fused phase. The unwrapping is performed to obtain the surface subsidence sequence in the mining monitoring area.

[0179] Furthermore, after unwrapping the fused phase, the method further includes: performing Loess smoothing on the unwrapping result to obtain a smoothed surface subsidence sequence in the mining monitoring area of ​​the mining area.

[0180] Optionally, in some other embodiments of the present invention, the unwrapping result may be smoothed using a Bisquare algorithm, an Adjacent-Averaging algorithm, or the like.

[0181] It should be understood that in the above-mentioned embodiments of the present invention, the unwrapping results can be subjected to Loess smoothing processing to obtain a smoothed surface subsidence sequence in the mining monitoring area of ​​the mining area, so as to improve the accuracy of the fused phase unwrapping results, thereby improving the mining monitoring efficiency and the reliability of the monitoring results, while reducing the mining monitoring cost.

[0182] In summary, according to the method for obtaining a surface subsidence sequence by fusing a spatiotemporal simulation phase and a DInSAR phase proposed in an embodiment of the present invention, the probability integral method model parameters and the time function model parameters of the mining monitoring area are obtained, and a spatiotemporal model of mining subsidence is constructed based on the probability integral method model parameters and the time function model parameters, and then the DInSAR phase of the mining monitoring area is obtained, and the simulation phase of the mining monitoring area is obtained based on the spatiotemporal model of mining subsidence, and the fused phase is obtained based on the DInSAR phase and the simulated phase, and the fused phase is unwrapped to obtain the surface subsidence sequence of the mining monitoring area. Thus, the surface subsidence sequence of the mining monitoring area is obtained based on the fused phase, and the accuracy of the phase unwrapping is improved to improve the resolution of the surface subsidence sequence, thereby reducing the mining monitoring cost while improving the mining monitoring efficiency and the reliability of the monitoring results.

[0183] It should be noted that the logic and / or steps represented in the flowchart or otherwise described herein, for example, can be considered as a sequenced list of executable instructions for implementing logical functions, and can be specifically implemented in any computer-readable medium for use by an instruction execution system, device or apparatus (such as a computer-based system, a system including a processor, or other system that can fetch instructions from an instruction execution system, device or apparatus and execute instructions), or in combination with these instruction execution systems, devices or apparatuses. For the purposes of this specification, "computer-readable medium" can be any device that can contain, store, communicate, propagate or transmit a program for use by an instruction execution system, device or apparatus, or in combination with these instruction execution systems, devices or apparatuses. More specific examples of computer-readable media (a non-exhaustive list) include the following: an electrical connection portion with one or more wirings (electronic device), a portable computer disk box (magnetic device), a random access memory (RAM), a read-only memory (ROM), an erasable and programmable read-only memory (EPROM or flash memory), an optical fiber device, and a portable compact disk read-only memory (CDROM). In addition, the computer-readable medium may even be paper or other suitable medium on which the program is printed, since the program may be obtained electronically, for example, by optically scanning the paper or other medium and then editing, interpreting or processing in other suitable ways if necessary, and then stored in a computer memory.

[0184] It should be understood that the various parts of the present invention can be implemented by hardware, software, firmware or a combination thereof. In the above-mentioned embodiments, a plurality of steps or methods can be implemented by software or firmware stored in a memory and executed by a suitable instruction execution system. For example, if implemented by hardware, as in another embodiment, it can be implemented by any one of the following technologies known in the art or their combination: a discrete logic circuit having a logic gate circuit for implementing a logic function for a data signal, a dedicated integrated circuit having a suitable combination of logic gate circuits, a programmable gate array (PGA), a field programmable gate array (FPGA), etc.

[0185] In the description of this specification, the description with reference to the terms "one embodiment", "some embodiments", "examples", "specific examples", or "some examples" means that the specific features, structures, materials or characteristics described in conjunction with the embodiment or example are included in at least one embodiment or example of the present invention. In this specification, the schematic representation of the above terms does not necessarily refer to the same embodiment or example. Moreover, the specific features, structures, materials or characteristics described may be combined in any one or more embodiments or examples in a suitable manner.

[0186] In the description of the present invention, it is to be understood that the terms “center”, “longitudinal”, “lateral”, “length”, “width”, “thickness”, “up”, “down”, “front”, “back”, “left”, “right”, “vertical”, “horizontal”, “top”, “bottom”, “inside”, “outside”, “clockwise”, “counterclockwise”, “axial”, “radial”, “circumferential”, etc., indicating orientations or positional relationships based on the orientations or positional relationships shown in the accompanying drawings, are only for the convenience of describing the present invention and simplifying the description, and do not indicate or imply that the referred device or element must have a specific orientation, be constructed and operated in a specific orientation, and therefore should not be understood as limiting the present invention.

[0187] In addition, the terms "first" and "second" are used for descriptive purposes only and should not be understood as indicating or implying relative importance or implicitly indicating the number of the indicated technical features. Therefore, the features defined as "first" and "second" may explicitly or implicitly include at least one of the features. In the description of the present invention, the meaning of "plurality" is at least two, such as two, three, etc., unless otherwise clearly and specifically defined.

[0188] In the present invention, unless otherwise clearly specified and limited, the terms "installed", "connected", "connected", "fixed" and the like should be understood in a broad sense, for example, it can be a fixed connection, a detachable connection, or an integral connection; it can be a mechanical connection or an electrical connection; it can be a direct connection or an indirect connection through an intermediate medium, it can be the internal connection of two elements or the interaction relationship between two elements, unless otherwise clearly defined. For ordinary technicians in this field, the specific meanings of the above terms in the present invention can be understood according to specific circumstances.

[0189] In the present invention, unless otherwise clearly specified and limited, a first feature being "above" or "below" a second feature may mean that the first and second features are in direct contact, or the first and second features are in indirect contact through an intermediate medium. Moreover, a first feature being "above", "above" or "above" a second feature may mean that the first feature is directly above or obliquely above the second feature, or simply means that the first feature is higher in level than the second feature. A first feature being "below", "below" or "below" a second feature may mean that the first feature is directly below or obliquely below the second feature, or simply means that the first feature is lower in level than the second feature.

[0190] Although the embodiments of the present invention have been shown and described above, it is to be understood that the above embodiments are exemplary and are not to be construed as limitations of the present invention. A person skilled in the art may change, modify, replace and vary the above embodiments within the scope of the present invention.

Claims

1. A method for obtaining a surface subsidence sequence by fusing a spatiotemporal simulation phase and a DInSAR phase, characterized in that: The method comprises: Obtaining probability integration method model parameters and time function model parameters of the mining monitoring area, and constructing a mining subsidence spatiotemporal model according to the probability integration method model parameters and the time function model parameters; Acquire the DInSAR phase of the mining monitoring area of ​​the mining area, and acquire the simulated phase of the mining monitoring area of ​​the mining area according to the mining subsidence spatiotemporal model of the mining area; Acquire a fused phase according to the DInSAR phase and the simulated phase, and unwrap the fused phase to obtain a surface subsidence sequence in the mining monitoring area of ​​the mining area; The acquiring a fused phase according to the DInSAR phase and the simulated phase comprises: spatially resampling the DInSAR phase and the simulated phase; Calculate the phase difference between the spatially resampled DInSAR phase and the simulated phase; Using a Kalman filter algorithm to filter the phase difference to obtain a filtered phase; The simulation phase and the filtering phase are summed to obtain the fusion phase.

2. The method for obtaining a surface subsidence sequence by fusing a spatiotemporal simulation phase and a DInSAR phase according to claim 1, characterized in that: The method of obtaining the probability integration method model parameters of the mining monitoring area includes: Acquire the level monitoring time series data and construct the inversion objective function of the probability integration method model parameters; According to the level monitoring time series data and the inversion objective function, using a plurality of preset algorithms, a plurality of groups of the probability integration method model parameters are obtained; A weighted combination prediction model is used to weightedly fuse multiple groups of probability integration method model parameters, and the weighted fused combination prediction parameters are used as the probability integration method model parameters of the mining monitoring area of ​​the mining area.

3. The method for obtaining a surface subsidence sequence by fusing a spatiotemporal simulation phase and a DInSAR phase according to claim 2, characterized in that: The method of obtaining the time function model parameters of the mining monitoring area includes: Determine a monitoring point timing start time of the level monitoring time series data, and adjust the level monitoring time series data according to the monitoring point timing start time; According to the start time of the monitoring point sequence and the adjusted level monitoring sequence data, the time function model parameters of the mining monitoring area of ​​the mining area are obtained by using the iterative function extreme value LM algorithm.

4. The method for obtaining a surface subsidence sequence by fusing a spatiotemporal simulation phase and a DInSAR phase according to claim 1, characterized in that: The step of acquiring the DInSAR phase of the mining monitoring area comprises: Obtaining historical radar monitoring data of the mining monitoring area of ​​the mining area; Generate an interferometric image pair corresponding to the historical radar monitoring data according to the DInSAR imaging time; The DInSAR phase of the mining monitoring area of ​​the mining area is obtained according to the interference image pair.

5. The method for obtaining a surface subsidence sequence by fusing a spatiotemporal simulation phase and a DInSAR phase according to claim 4, characterized in that: The step of obtaining the DInSAR phase of the mining monitoring area according to the interference image pair includes: The interference image pair is subjected to image cropping, precise orbit data loading, flat ground phase removal, terrain phase removal, interference, phase filtering, and geocoding in sequence to obtain the DInSAR phase of the mining monitoring area of ​​the mining area.

6. The method for obtaining a surface subsidence sequence by fusing a spatiotemporal simulation phase and a DInSAR phase according to claim 4, characterized in that: Acquiring a simulated phase of the mining monitoring area of ​​the mining area according to the mining subsidence spatiotemporal model of the mining area includes: Calculate the vertical settlement of the ground surface corresponding to the DInSAR imaging time according to the mining subsidence spatiotemporal model of the mining area; Converting the vertical settlement of the ground surface into a change in slant range in the radar line of sight; The slant distance variation is phase-converted to obtain a simulated phase of the mining monitoring area of ​​the mining area.

7. The method for obtaining a surface subsidence sequence by fusing a spatiotemporal simulation phase and a DInSAR phase according to claim 1, characterized in that: After acquiring the fusion phase, the method further includes: The fused phase is corrected to the interval [-π,π].

8. The method for obtaining a surface subsidence sequence by fusing a spatiotemporal simulation phase and a DInSAR phase according to claim 1, characterized in that: The unwrapping of the fused phase includes: The fused phase is unwrapped using a fast two-dimensional phase unwrapping algorithm based on non-continuous path reliability sorting.

9. The method for obtaining a surface subsidence sequence by fusing a spatiotemporal simulation phase and a DInSAR phase according to claim 8, characterized in that: After unwrapping the fused phase, the method further includes: The unwrapping results are subjected to Loess smoothing to obtain a smoothed surface subsidence sequence in the mining monitoring area of ​​the mining area.

Images