CN103091676A - Mining area surface subsidence synthetic aperture radar interferometry monitoring and calculating method - Google Patents

Abstract

A mining area surface subsidence synthetic aperture radar interferometry monitoring and calculating method belongs to mining area surface subsidence monitoring and calculating method. The method comprises the following steps: having format conversion, calibration, pre-filtering and interference to interferometric synthetic aperture radar (InSAR) data to obtain an InSAR interference phase, eliminating flat ground effect, an terrain phase and orbit errors of the interference phase by means of precise orbit data and an external dynamic effect model (DEM), obtaining phase values which only contain deformation information of a ground surface after filtering to residual phases, on the condition of large deformation gradient, through phase unwrapping, obtaining deflection of the edge of a mining subsidence basin, fusing the deflection with a few ground measured data, reversely deducing probability integral method parameters of the subsidence basin through a genetic algorithm, calculating sinking values of an arbitrary point of the whole surface subsidence basin through required probability integral method parameters and geological mining data, and therefore mining subsidence deformation field is produced. The mining area surface subsidence synthetic aperture radar interferometry monitoring and calculating method has the advantages of being high in monitoring accuracy, large in monitoring range, convenient to operate, low in cost and large in technical content.

Description

Monitoring and resolving method for mining area surface mining subsidence synthetic aperture radar interferometry

Technical Field

The invention relates to a monitoring and resolving method for mining area surface mining subsidence, in particular to a monitoring and resolving method for mining area surface mining subsidence synthetic aperture radar interferometry.

Technical Field

The synthetic aperture radar interferometry is used as an all-weather, large-visible-range, stable-track and high-precision geodetic surveying technology which cannot be compared with optical remote sensing, the application field is continuously expanded, and the synthetic aperture radar interferometry is always a popular research subject at home and abroad. With the continuous updating of the space-based radar system and the great abundance of SAR data, the technology has obtained unprecedented development and application. The synthetic aperture radar interferometry technique is expressed in english as: InSAR, Interferometric Synthetic Aperture Radar.

The "two-track difference" method in the dinsar (differential Interferometric Synthetic Aperture radar) is a new space-to-ground observation technology which takes two SAR images in the same area as basic processing data, obtains an interference image by solving the phase difference of the two SAR images, removes a terrain phase by an external dem (digital Elevation model) and obtains surface deformation information. The DInSAR has great application potential in the fields of volcano monitoring, ground settlement, seismic deformation field acquisition, landslide monitoring and the like.

In the case of a repetitive orbit, if the earth surface deformation moves within the time interval between the two image acquisitions, the interference fringes formed by the two images mainly contain the following phase information (fig. 1):

in the formula,

is the terrain phase;

is a satellite sight squarePhase of deformation to the earth's surface;is a phase resulting from atmospheric delay or the like;

is the phase due to the reference plane;

is the phase caused by noise.

Flat ground phase calculation formula:

${\mathrm{\φ}}_{\mathrm{fat}}=-\frac{4\mathrm{\π}}{\mathrm{\λ}}{B}_1=-\frac{4\mathrm{\π}}{\mathrm{\λ}}B\sin ({\mathrm{\θ}}_0-a)$

where λ is the radar wavelength, B is the spatial baseline, and θ₀The radar incidence angle is defined, and alpha is an included angle between a space baseline and the horizontal direction;

topographic phase calculation formula:

in the formula,

is a vertical baseline;

the view angle difference of the radar sight line caused by the elevation of the ground point;

atmospheric retardation phase: the residual phase and an external DEM are subjected to linear regression resolving to obtain the phase-locked loop;

noise phase: the cancellation is performed by Goldstein frequency domain filtering.

Deformation phase: after removing the phase component, the final deformation phase is obtained, and the formula is as follows:

maximum deformation phase gradient D available by InSAR technology_LThe method comprises the following steps:

$D_L=\frac{\mathrm{\λ}}{2\mathrm{\μ}}$

in the formula, mu is the side length of the pixel, and lambda is the wavelength.

For C-band SAR images of ERS and Envisat satellites, the wavelength is 56 mm, the resolution is 20 m, and D_LIs 0.0028; l-band SAR image of ALOS satellite with wavelength of 230 mm and resolution of 10 m and D_L0.023; the X wave band SAR image of the Terras SAR-X satellite has the wavelength of 32 mm and the resolution of 2 m, D_LIs 0.016.

The InSAR technology has self limitation, is easily influenced by temporal and spatial incoherence, and has quite harsh application conditions. Particularly, the ground surface vegetation cover of most mining areas in China is large, the relief of the ground surface is large, the subsidence caused by resource exploitation has very large deformation gradient, the correct solution is difficult to adopt the existing InSAR technology, and the application of the InSAR technology in the ground surface subsidence monitoring of the mining areas is greatly limited. And mineral resources in China are rich, mining modes are relatively extensive, illegal mining and illegal digging sometimes occur, so that disasters such as ground surface settlement of a mining area are serious, and the disasters cannot be completely monitored one by adopting a conventional geodetic measurement means.

Disclosure of Invention

The invention aims to provide a monitoring and resolving method for mining area surface mining subsidence synthetic aperture radar interferometry, which solves the problem that the real settlement of the surface cannot be accurately obtained due to phase loss coherence when the mining area surface subsidence is monitored by the synthetic aperture radar interferometry technology InSAR, and also solves the problem that more surface control points are needed when the probability integration method parameter inversion is carried out in the mining subsidence prediction theory.

The aim of the invention is achieved by the following technical solution:

a monitoring and resolving method for mining area surface mining subsidence synthetic aperture radar interferometry comprises the following steps:

1) generating an interference pattern: the method comprises the steps of selecting an InSAR image, converting a format, registering the image, pre-filtering and generating an interference pattern so as to obtain an interference phase pattern of a monitored area; the registration process is as follows: performing rough rail registration and fine pixel-level and sub-pixel-level registration, calculating the offset of the secondary image relative to the primary image, performing polynomial fitting, and completing resampling of the secondary image according to a polynomial fitting coefficient; the interference process is as follows: conjugate multiplying the corresponding pixels of the main image and the resampled auxiliary image to obtain an interference fringe pattern;

2) resolving a deformation phase: eliminating a flat ground effect, an orbit error and a terrain phase in an interference phase by means of precise orbit data and an external DEM (digital elevation model), and filtering a residual phase to obtain a phase value only containing earth surface deformation information; the DEM is an external Digital Elevation Model, and English is expressed as a Digital Elevation Model;

the interference image contains the following phase information:

in the formula,

is the terrain phase;

is the earth surface deformation phase in the satellite sight direction;

is a phase resulting from atmospheric delay or the like;

is the phase due to the reference plane;is the phase caused by noise;

the land phase is calculated using the following equation:

${\mathrm{\φ}}_{\mathrm{fat}}=-\frac{4\mathrm{\π}}{\mathrm{\λ}}{B}_1=-\frac{4\mathrm{\π}}{\mathrm{\λ}}B\sin ({\mathrm{\θ}}_0-a)$

wherein λ is radar wavelength, B is space baseline, and θ₀The radar incidence angle is defined, and alpha is an included angle between a space baseline and the horizontal direction;

the terrain phase is calculated using the following equation:

wherein,

is a vertical baseline;

the view angle difference of the radar sight line caused by the elevation of the ground point;

atmospheric retardation phase: the residual phase and an external DEM are subjected to linear regression resolving to obtain the phase-locked loop;

noise phase: elimination is carried out through Goldstein frequency domain filtering;

deformation phase: after removing the phase component, the final deformation phase is obtained, and the formula is as follows:

please explain two unknowns in the formula.

Wherein, λ is radar wavelength, and Δ r is deformation in radar sight direction;

3) phase unwrapping and geocoding: we unwrap the phase of deformation of the winding using a least cost-stream method and convert the unwrapped phase to a surface subsidence. That is, the final surface sedimentation amount formula is:

in the formula,

is the deformation phase after unwrapping;

with the assistance of an external DEM, projecting the earth surface settlement under a WGS-84 coordinate system;

4) selecting reliable sinking points on the image: the step is to unwind the winding deformation phase by adopting a minimum cost flow method, convert the unwinding phase into the ground surface subsidence, and finally select the subsidence point position with more reliable subsidence basin edge according to the coherence of the point position.

5) Calculating a probability integration method parameter: fusing a small amount of measured data of the earth surface mobile observation station with image reliable point location sinking data to form an earth surface control reference point set with predicted parameters, and continuously and circularly solving through steps of crossing, variation and the like based on a genetic optimization algorithm to calculate final probability integration method parameters;

the most important step of the probability integration method parameter inversion based on the genetic optimization algorithm is to determine an encoding rule and a fitness function, and the method specifically comprises the following steps:

A. and (3) encoding: calculating the parameters of the probability integration method by adopting a genetic algorithm, wherein each parameter needs to be subjected to coding of a chromosome structure, and the coding rule can be selected from binary coding and real number coding; because the probability integration method parameter inversion is a complex nonlinear optimization problem, the adoption of binary coding can affect the calculation precision and the calculation efficiency of the evolutionary algorithm, and the real number coding is suitable for genetic algorithms with large range and high precision, so the real number coding is selected; arranging the probability integration method parameters into a chromosome string of a genetic algorithm according to the sequence, wherein each gene in the string represents one parameter of the probability integration method;

B. fitness function: combining the phase unwrapping characteristic of the InSAR technology, selecting the expected subsidence value of the earth surface control reference point set to be consistent with the actual measured subsidence values corresponding to the point positions, determining a fitness function, evaluating the fitness of each individual of the group in the calculation, sequencing the fitness from large to small, and eliminating the individual with small fitness;

C. selecting the individual with the maximum fitness value, and directly transmitting the individual to the next generation; carrying out genetic operation on the current generation population by using operation operators such as crossover, mutation and the like to generate a next generation population;

D. repeating the steps B-C to continuously optimize the calculation result of the probability integration method parameters until the termination condition is met; the obtained parameter genetic code is decoded to obtain a group of parameter sequences of a probability integration method;

6) resolving a ground surface subsidence basin: and (3) jointly resolving the subsidence value of any point of the whole surface subsidence basin by utilizing the probability integration method parameters and geological mining data, thereby generating a mining subsidence deformation field in the mining area.

Six parameters of the surface deformation predicted by the probability integration method are a sinking coefficient q, a main influence radius r, a main influence angle tangent tg beta and an inflection point offset S₀Horizontal migration coefficient b, mining influence propagation angle theta₀Here, the subsidence value of any point of the subsidence basin is solved, and according to the basic principle of the probability integration method model, the subsidence value of any point (x, y) on the earth surface caused by underground mining can be expressed as follows:

$W(x,y)=\mathrm{mq}\cos \mathrm{\α}\∫\∫W_{\mathrm{oci}}(x,y)\mathrm{dxdy}$

$W_{\mathrm{coi}}(x,y)=\frac{\mathrm{<figure-callout\; id="2"\; label="l\; r"\; filenames="HDA00002759719200011.png"\; state="\{\{state\}\}">l</figure-callout>}}{r^{}}\mathrm{epx}\left(\frac{-\mathrm{\&pi;}{(x-x_i)}^2}{r^2}\right)-\exp \left(\frac{-\mathrm{\&pi;}{(y-y_i+H_i\mathrm{ctg\&theta;})}^2}{r^2}\right)$

in the formula: q is a sinking coefficient; r is the major radius of influence, r = H₀/tgβ;H₀Average mining depth; theta is a mining influence propagation angle; tg β is the primary influence angle β tangent; (xi, yi) is a plane coordinate of the central point of the mining unit i; (x, y) is the coordinates of an arbitrary point on the earth's surface.

The method has the advantages that due to the adoption of the scheme, the method sufficiently fuses relevant theories of InSAR and mining subsidence, overcomes the problem that the traditional InSAR method cannot accurately acquire the surface subsidence under the large deformation gradient of the mining area, and solves the problem that more surface control points are needed during the parameter inversion of the probability integration method in the mining subsidence prediction theory to a certain extent.

The advantages are that: the method has the advantages of clear structure of the whole process, high monitoring precision, large range, simple realization of operation process, low cost and large technical content, overcomes the problem that the traditional InSAR method cannot accurately acquire the surface subsidence under the large deformation gradient of the mining area, solves the problem that more surface control points are needed during the probability integration method parameter inversion in the mining subsidence prediction theory to a certain extent, and provides a solution for the application and popularization of the InSAR technology in the mining area.

Drawings

FIG. 1 is a schematic diagram of interferometric phase composition analysis according to an embodiment of the present invention.

Fig. 2 is a flow chart of an InSAR monitoring and resolving method for mining area surface mining subsidence according to an embodiment of the invention.

FIG. 3 is a face mining simulated subsidence basin map.

FIG. 4 is a phase diagram of simulated subsidence unwrapping for face mining using conventional methods.

FIG. 5 is a plot of the results of the resolution of the face mining simulation subsidence of the present invention.

FIG. 6 is a diagram of the results of a conventional method of resolving actual work surface subsidence.

FIG. 7 is a graph of the results of the method of the present invention in resolving actual work surface subsidence.

FIG. 8 is a comparison graph of the calculation result of the present invention and the measured data of the ground control point.

Detailed Description

The invention will be described in further detail below with reference to the figures and specific implementation processes:

example 1: a monitoring and resolving method for mining area surface mining subsidence synthetic aperture radar interferometry comprises the following steps:

1) generating an interference pattern: the method comprises the steps of selecting an InSAR image, converting a format, registering the image, pre-filtering and generating an interference pattern so as to obtain an interference phase pattern of a monitored area; the registration process is as follows: performing rough rail registration and fine pixel-level and sub-pixel-level registration, calculating the offset of the secondary image relative to the primary image, performing polynomial fitting, and completing resampling of the secondary image according to a polynomial fitting coefficient; the interference process is as follows: conjugate multiplying the corresponding pixels of the main image and the resampled auxiliary image to obtain an interference fringe pattern;

2) resolving a deformation phase: eliminating a flat ground effect, an orbit error and a terrain phase in an interference phase by means of precise orbit data and an external DEM (digital elevation model), and filtering a residual phase to obtain a phase value only containing earth surface deformation information; the DEM is as follows: an external Digital Elevation Model, denoted Digital Elevation Model in English;

the interference image contains the following phase information:

in the formula,is the terrain phase;

is a satellite sight squarePhase of deformation to the earth's surface;

is a phase resulting from atmospheric delay or the like;

is the phase due to the reference plane;

is the phase caused by noise;

the land phase is calculated using the following equation:

${\mathrm{\&phi;}}_{\mathrm{fat}}=-\frac{4\mathrm{\&pi;}}{\mathrm{\&lambda;}}{B}_1=-\frac{4\mathrm{\&pi;}}{\mathrm{\&lambda;}}B\sin ({\mathrm{\&theta;}}_0-a)$

wherein λ is radar wavelength, B is space baseline, and θ₀The radar incidence angle is defined, and alpha is an included angle between a space baseline and the horizontal direction;

the terrain phase is calculated using the following equation:

wherein,

is a vertical baseline;

the view angle difference of the radar sight line caused by the elevation of the ground point;

atmospheric retardation phase: the residual phase and an external DEM are subjected to linear regression resolving to obtain the phase-locked loop;

noise phase: elimination is carried out through Goldstein frequency domain filtering; the Goldstein frequency domain filtering is a frequency domain filtering method, an interferogram is converted into a frequency domain from a space domain by utilizing Fourier transform, power spectrum smoothing is carried out, and denoising in a signal region is realized by smoothing the frequency domain;

deformation phase: after removing the phase component, the final deformation phase is obtained, and the formula is as follows:

please explain two unknowns in the formula.

Wherein, λ is radar wavelength, and Δ r is deformation in radar sight direction;

3) phase unwrapping and geocoding: we unwrap the phase of deformation of the winding using a least cost-stream method and convert the unwrapped phase to a surface subsidence. That is, the final surface sedimentation amount formula is:

in the formula,

is the deformation phase after unwrapping;

with the assistance of an external DEM, projecting the earth surface settlement under a WGS-84 coordinate system;

4) selecting reliable sinking points on the image: the step is to unwind the winding deformation phase by adopting a minimum cost flow method, convert the unwinding phase into the ground surface subsidence, and finally select the subsidence point position with more reliable subsidence basin edge according to the coherence of the point position.

5) Calculating a probability integration method parameter: fusing a small amount of measured data of the earth surface mobile observation station with image reliable point location sinking data to form an earth surface control reference point set with predicted parameters, and continuously and circularly solving through steps of crossing, variation and the like based on a genetic optimization algorithm to calculate final probability integration method parameters;

the most important step of solving the probability basic method parameters by adopting the genetic optimization algorithm is to determine a coding rule and a fitness function, and the method specifically comprises the following steps:

A. and (3) encoding: calculating the parameters of the probability integration method by adopting a genetic algorithm, wherein each parameter needs to be subjected to coding of a chromosome structure, and the coding rule can be selected from binary coding and real number coding; because the probability integration method parameter inversion is a complex nonlinear optimization problem, the adoption of binary coding can affect the calculation precision and the calculation efficiency of the evolutionary algorithm, and the real number coding is suitable for genetic algorithms with large range and high precision, so the real number coding is selected; arranging the probability integration method parameters into a chromosome string of a genetic algorithm according to the sequence, wherein each gene in the string represents one parameter of the probability integration method;

B. fitness function: combining the phase unwrapping characteristic of the InSAR technology, selecting the expected subsidence value of the earth surface control reference point set to be consistent with the actual measured subsidence values corresponding to the point positions, determining a fitness function, evaluating the fitness of each individual of the group in the calculation, sequencing the fitness from large to small, and eliminating the individual with small fitness;

C. selecting the individual with the maximum fitness value, and directly transmitting the individual to the next generation; carrying out genetic operation on the current generation population by using operation operators such as crossover, mutation and the like to generate a next generation population;

D. repeating the steps B-C to continuously optimize the calculation result of the probability integration method parameters until the termination condition is met; the obtained parameter genetic code is decoded to obtain a group of parameter sequences of a probability integration method;

6) resolving a ground surface subsidence basin: calculating the subsidence value of any point of the whole surface subsidence basin by utilizing the probability integration method parameters and geological mining data in a combined manner, thereby generating a mining subsidence deformation field of the mining area;

a method for calculating the prediction of the earth surface subsidence basin by adopting a probability integration method, wherein six parameters of the earth surface deformation predicted by the probability integration method are a subsidence coefficient q, a main influence radius r, a main influence angle tangent tg beta and a inflection point offset S₀Horizontal migration coefficient b, mining influence propagation angle theta₀Here, the subsidence value of any point of the subsidence basin is solved, and according to the basic principle of the probability integration method model, the subsidence value of any point (x, y) on the earth surface caused by underground mining can be expressed as follows:

$W(x,y)=\mathrm{mq}\cos \mathrm{\&alpha;}\&Integral;\&Integral;W_{\mathrm{oci}}(x,y)\mathrm{dxdy}$

$W_{\mathrm{coi}}(x,y)=\frac{\mathrm{<figure-callout\; id="2"\; label="l\; r"\; filenames="HDA00002759719200011.png"\; state="\{\{state\}\}">l</figure-callout>}}{r^{}}\mathrm{epx}\left(\frac{-\mathrm{\&pi;}{(x-x_i)}^2}{r^2}\right)-\exp \left(\frac{-\mathrm{\&pi;}{(y-y_i+H_i\mathrm{ctg\&theta;})}^2}{r^2}\right)$

in the formula: q is a sinking coefficient; r is the major radius of influence, r = H₀/tgβ;H₀Average mining depth; theta is a mining influence propagation angle; tg β is the primary influence angle β tangent; (xi, yi) is a plane coordinate of the central point of the mining unit i; (x, y) is the coordinates of an arbitrary point on the earth's surface.

In fig. 2, a flow of a method for monitoring and resolving an InSAR for mining subsidence on the surface of a mining area includes the following steps:

(1) selection and format conversion of SAR data. The preconditions for phase unwrapping using the probabilistic integration model are: each phase value is true and reliable, so that the probability integration method can be used for predicting the mining area more ideally. This prediction is a static prediction, i.e. the final subsidence basin formed by the surface after face extraction, which takes a certain amount of time, typically extending over 1 year. The real SAR image is influenced by a time baseline and a space baseline, so that interference is difficult to form; secondly, the phase information of the pixel is difficult to reach an ideal state due to the influence of thermal noise, atmospheric delay and the like; in addition, the working face is generally under the farmland, and the ground surface subsidence deformation gradient caused by the mining of the working face is large, so that the pixel coherence can not be effectively ensured. Therefore, there is a great limitation in obtaining images.

In reality, the amount of subsidence of the ground surface of the mining area is very large (for example, the surface of some working surfaces sinks by 3m in 11 days), and the problems of image cost and the like are also considered, and SAR data with high resolution, short satellite revisit period and long wavelength is selected to monitor large deformation of the mining area during data selection. After the SAR image is selected, the original data format of the SAR needs to be converted into a single-view complex image (SLC).

(2) And (5) registering the SAR images. The currently commonly used registration method is a three-step registration method: coarse registration, wherein the registration precision is about 30 pixels; secondly, pixel level registration; and thirdly, subpixel level registration. The method comprises fitting a biquadratic model, resampling and the like, generally speaking, the sub-pixel level registration accuracy reaches 1/8 pixels, when the registration accuracy is superior to 1/8 pixels, the caused decorrelation is small (about 4%), and the accuracy requirement of SAR interference processing is met.

Optionally, (3) pre-filtering the single-view complex image. Since the SAR image has spectral shifts in both the distance direction and the azimuth direction, phase noise is introduced into the interferogram, and therefore, it is necessary to perform pre-filtering in the distance direction and the azimuth direction before generating the interferogram. The azimuth filtering is a filtering process performed on the master-slave image in the azimuth direction so as to retain the same doppler spectrum. Distance pre-filtering refers to the process of removing the local spectral shift between the master and slave images from the local interferogram and then filtering out the intra-spectral noise using a band-pass filter.

(4) An interferogram is generated. And carrying out conjugate multiplication on corresponding pixels of the master image and the slave image so as to obtain an interference pattern. The result of conjugate multiplication is a complex number whose modulus values are called the interference intensity map and whose phase values are called the fringe map or interferogram. The phase values are here wound, the absolute values of which are not greater than pi.

(5) And resolving the deformation phase. The method comprises the steps of eliminating the flat ground effect, the orbit error and the terrain phase in the interference phase by means of precise orbit data and an external DEM, and obtaining a phase value only containing earth surface deformation information after filtering the residual phase.

The land phase is calculated using the following equation:

${\mathrm{\&phi;}}_{\mathrm{fat}}=-\frac{4\mathrm{\&pi;}}{\mathrm{\&lambda;}}{B}_1=-\frac{4\mathrm{\&pi;}}{\mathrm{\&lambda;}}B\sin ({\mathrm{\&theta;}}_0-a)$

wherein lambda is the radar wavelength, B is a space baseline, theta 0 is the radar incident angle, and alpha is the included angle between the space baseline and the horizontal direction;

the terrain phase is calculated using the following equation:

is a vertical baseline;

the view angle difference of the radar sight line caused by the elevation of the ground point;

atmospheric retardation phase: the residual phase and an external DEM are subjected to linear regression resolving to obtain the phase-locked loop;

noise phase: the cancellation is performed by Goldstein frequency domain filtering.

Deformation phase: after removing the phase component, the final deformation phase is obtained, and the formula is as follows:

(6) phase unwrapping and geocoding: we unwrap the phase of deformation of the winding using a least cost-stream method and convert the unwrapped phase to a surface subsidence. That is, the final surface sedimentation amount formula is:

in the formula,

the phase of deformation after unwrapping.

With the aid of an external DEM, the surface subsidence is again projected under the WGS-84 coordinate system.

(7) Selecting reliable sinking points on the image: in practical application, because of the influence of a plurality of incoherent factors of the SAR image, the large ground deformation gradient, the large vegetation coverage and other factors existing in a mining area, the ground subsidence calculated by the InSAR technology is not completely correct, and sometimes the subsidence in a basin is even completely wrong. The more reliable sinking point locations at the edges of the sinking basins need to be selected according to the coherence of the point locations, wherein the larger the coherence of the point locations is, the higher the reliability of the sinking information of the sinking point locations is, and the more accurate the resolved sinking basins are.

(8) Calculating a probability integration method parameter: a small amount of measured data of the earth surface mobile observation station and image reliable point location sinking data are fused to form an earth surface control reference point set with predicted parameters, and final probability integration method parameters are calculated through continuous loop iterative solution by utilizing a genetic optimization algorithm and through steps of crossing, variation and the like. The most important step of the parameter calculation algorithm is to determine a coding rule and a fitness function, and the method specifically comprises the following steps:

A. and (5) encoding. The method adopts a genetic algorithm to calculate the parameters of the probability integration method, and each parameter needs to be subjected to coding of a chromosome structure, and the coding rule can select binary coding and real number coding. Because the probability integration method parameter inversion is a complex nonlinear optimization problem, the adoption of binary coding can affect the calculation precision and the calculation efficiency of the evolutionary algorithm, and the real number coding is suitable for genetic algorithms with large range and high precision, so the real number coding is selected. The probability integration method parameters are arranged in order as a chromosome string of a genetic algorithm, and each gene in the string represents one parameter of the probability integration method.

B. A fitness function. And (3) selecting the expected sinking values of the earth surface control reference point set to be consistent with the actual measurement sinking values corresponding to the point positions by combining the phase unwrapping characteristic of the InSAR technology to determine a fitness function, evaluating the fitness of each individual of the group in the calculation, sequencing the fitness from large to small, and eliminating the individuals with small fitness.

C. And selecting the individual with the maximum fitness value, and directly transmitting the individual to the next generation. And carrying out genetic operation on the current generation population by using operation operators such as crossover, mutation and the like to generate a next generation population.

D. And repeating the steps B-C to continuously optimize the calculation result of the probability integration method parameters until the termination condition is met. And decoding the obtained parameter genetic code to obtain a group of parameter sequences of a probability integration method.

[0081](9) Resolving a ground surface subsidence basin: calculating the subsidence value of any point of the whole earth surface subsidence basin by utilizing the probability integration method parameters and geological mining data in a combined manner, thereby generating oresAnd (5) mining the subsidence deformation field. Six parameters of the surface deformation predicted by the probability integration method are a sinking coefficient q, a main influence radius r, a main influence angle tangent tg beta and an inflection point offset S₀The horizontal movement coefficient b and the mining influence propagation angle theta 0 are mainly used for calculating the subsidence value of any point of a subsidence basin, and according to the basic principle of a probability integration method model, the subsidence value of any point (x, y) on the earth surface caused by underground mining can be expressed as follows:

$W(x,y)=\mathrm{mq}\cos \mathrm{\&alpha;}\&Integral;\&Integral;W_{\mathrm{oci}}(x,y)\mathrm{dxdy}$

$W_{\mathrm{coi}}(x,y)=\frac{\mathrm{<figure-callout\; id="2"\; label="l\; r"\; filenames="HDA00002759719200011.png"\; state="\{\{state\}\}">l</figure-callout>}}{r^{}}\mathrm{epx}\left(\frac{-\mathrm{\&pi;}{(x-x_i)}^2}{r^2}\right)-\exp \left(\frac{-\mathrm{\&pi;}{(y-y_i+H_i\mathrm{ctg\&theta;})}^2}{r^2}\right)$

in the formula: q is a sinking coefficient; r is the major radius of influence, r = H₀/tgβ；H₀Average mining depth; theta is a mining influence propagation angle; tg β is the primary influence angle β tangent; (xi, yi) is a plane coordinate of the central point of the mining unit i; (x, y) is the coordinates of an arbitrary point on the earth's surface.

The maximum sinking value under the full mining condition is

W₀＝qmcosα

In the formula, q is a sinking coefficient; m is the mining thickness of the coal bed; alpha is the coal seam dip angle (in degrees).

In fact, the surface subsidence prediction area is an insufficient mining area in most cases, and the surface subsidence prediction area can be regarded as superposition of two semi-infinite mining according to the probability integration theory without loss of generality. In semi-infinite mining in inclined coal seam mining, mining of a unit [ coordinate (s, t) ] causes the subsidence value of any point [ coordinate (x, y) ] of the earth surface to be:

$W(x,y)={\mathrm{<figure-callout\; id="0"\; label="W"\; filenames="HDA00002759719200021.png"\; state="\{\{state\}\}">W</figure-callout>}}_0{\&Integral;}_0^{+\&infin;}\frac{1}{{\mathrm{<figure-callout\; id="2"\; label="r"\; filenames="HDA00002759719200011.png"\; state="\{\{state\}\}">r</figure-callout>}}^2}\exp (-\mathrm{\&pi;}\frac{{(x-s)}^2}{r^2})\mathrm{ds}-{\&Integral;}_0^{+\&infin;}\frac{1}{{\mathrm{<figure-callout\; id="2"\; label="r"\; filenames="HDA00002759719200011.png"\; state="\{\{state\}\}">r</figure-callout>}}^2}\exp (-\mathrm{\&pi;}\frac{{(y-r+l)}^2}{r^2})\mathrm{dt}$

wherein r is the major influence radius; h is the depth of cut, 1 ═ H-cot theta

The mining of the coal seam is equivalent to the geometric superposition of two semi-infinite mining during the limited mining, and the calculation formula of the surface subsidence value is as follows:

in the formula, D₁、D₂The working face has long tendency and long trend respectively; s_{On the upper part}、S_{Lower part}、S_{Left side of}、S_{Right side}The inflection point offset distances of the upper, lower, left and right boundaries are respectively.

Compared with the existing InSAR surface subsidence monitoring method, the embodiment of the invention can solve the problem that the traditional InSAR method cannot accurately acquire the surface subsidence under the large deformation gradient of the mining area, and simultaneously solves the problem that more surface control points are needed during the probability integration method parameter inversion in the mining subsidence prediction theory to a certain extent.

In order to better illustrate the effectiveness and superiority of the technical solution of the embodiment of the present invention, the simulation of the mining subsidence data and the actual mining subsidence data are now performed. As shown in fig. 3, (a) is the ground surface subsidence caused by the horizontal face mining with the subsidence coefficient of 0.8 simulated by the present invention, the maximum subsidence is about 1350mm, the simulated subsidence is converted into the winding phase according to the SAR signal parameter of the C-band, the phase unwrapping by the conventional InSAR method as shown in (b) is difficult to obtain the real ground surface deformation phase, and the method of the present invention as shown in (C) can correctly calculate the ground surface subsidence. As shown in fig. 4, it is a comparison graph of the results of resolving the surface subsidence of the actual working surface of a certain mine, and it is difficult to obtain the true phase of the surface deformation and even cause errors by the phase unwrapping using the conventional InSAR method as shown in (a), while the method of the present invention as shown in (b) can correctly resolve the surface subsidence. As shown in FIG. 5, the measured data of the surface monitoring point of the working surface is well matched with the calculation result of the invention, and the effectiveness and the superiority of the invention are also illustrated.

Claims (4)

1. A monitoring and resolving method for mining subsidence synthetic aperture radar interferometry of surface mining in a mining area is characterized in that,

the method comprises the following steps:

1) generating an interference pattern: the method comprises the steps of selecting an InSAR image, converting a format, registering the image, pre-filtering and generating an interference pattern so as to obtain an interference phase pattern of a monitored area; the registration process is as follows: performing rough rail registration and fine pixel-level and sub-pixel-level registration, calculating the offset of the secondary image relative to the primary image, performing polynomial fitting, and completing resampling of the secondary image according to a polynomial fitting coefficient; the interference process is as follows: conjugate multiplying the corresponding pixels of the main image and the resampled auxiliary image to obtain an interference fringe pattern;

2) resolving a deformation phase: eliminating a flat ground effect, an orbit error and a terrain phase in an interference phase by means of precise orbit data and an external DEM (digital elevation model), and filtering a residual phase to obtain a phase value only containing earth surface deformation information;

the interference image contains the following phase information:

in the formula,

is the terrain phase;is the earth surface deformation phase in the satellite sight direction;is a phase resulting from atmospheric delay or the like;

is the phase due to the reference plane;

is the phase caused by noise;

flat ground phase calculation formula:

${\mathrm{\&phi;}}_{\mathrm{fat}}=-\frac{4\mathrm{\&pi;}}{\mathrm{\&lambda;}}{B}_1=-\frac{4\mathrm{\&pi;}}{\mathrm{\&lambda;}}B\sin ({\mathrm{\&theta;}}_0-a)$

wherein λ is radar wavelength, B is space baseline, and θ₀The radar incidence angle is defined, and alpha is an included angle between a space baseline and the horizontal direction;

topographic phase calculation formula:

wherein,

is a vertical baseline;the view angle difference of the radar sight line caused by the elevation of the ground point;

atmospheric retardation phase: the residual phase and an external DEM are subjected to linear regression resolving to obtain the phase-locked loop;

noise phase: elimination is carried out through Goldstein frequency domain filtering;

deformation phase: after removing the phase component, the final deformation phase is obtained

The formula is as follows:

wherein, λ is radar wavelength, and Δ r is deformation in radar sight direction;

3) phase unwrapping and geocoding: the deformation phase of the winding is unwound by adopting a minimum cost flow method, and the unwinding phase is converted into the ground subsidence; that is, the final surface sedimentation amount formula is:

in the formula,is the deformation phase after unwrapping;

with the assistance of an external DEM, projecting the earth surface settlement under a WGS-84 coordinate system;

4) selecting reliable sinking points on the image: the step is to unwind the winding deformation phase by adopting a minimum cost flow method, convert the unwinding phase into the ground surface subsidence, and finally select the subsidence point position with more reliable subsidence basin edge according to the coherence of the point position.

5) Calculating a probability integration method parameter: fusing a small amount of measured data of the earth surface mobile observation station with the image reliable point phase to form an earth surface control reference point set with predicted parameters, and continuously performing loop iteration solution to calculate final probability integration method parameters by using a genetic optimization algorithm and through the steps of crossing, variation and the like;

the most important step of utilizing the probability integration method parameter inversion based on the genetic optimization algorithm is to determine the encoding rule and the fitness function, and the method specifically comprises the following steps:

[0044] a, encoding: calculating the parameters of the probability integration method by adopting a genetic algorithm, wherein each parameter needs to be subjected to coding of a chromosome structure, and the coding rule can be selected from binary coding and real number coding; because the probability integration method parameter inversion is a complex nonlinear optimization problem, the adoption of binary coding can affect the calculation precision and the calculation efficiency of the evolutionary algorithm, and the real number coding is suitable for genetic algorithms with large range and high precision, so the real number coding is selected; arranging the probability integration method parameters into a chromosome string of a genetic algorithm according to the sequence, wherein each gene in the string represents one parameter of the probability integration method;

[0045] fitness function: combining the phase unwrapping characteristic of the InSAR technology, selecting the expected subsidence value of the earth surface control reference point set to be consistent with the actual measured subsidence values corresponding to the point positions, determining a fitness function, evaluating the fitness of each individual of the group in the calculation, sequencing the fitness from large to small, and eliminating the individual with small fitness;

[0046] c, selecting the individual with the maximum fitness value, and directly transmitting the individual to the next generation; carrying out genetic operation on the current generation population by using operation operators such as crossover, mutation and the like to generate a next generation population;

[0047] d, repeating the steps B-C to continuously optimize the calculation result of the probability integration method parameters until the termination condition is met; the obtained parameter genetic code is decoded to obtain a group of parameter sequences of a probability integration method;

6) resolving a ground surface subsidence basin: and (3) jointly resolving the subsidence value of any point of the whole surface subsidence basin by utilizing the probability integration method parameters and geological mining data, thereby generating a mining subsidence deformation field in the mining area.

2. The method for monitoring and calculating interferometry of synthetic aperture radar for mining area surface mining subsidence according to claim 1, wherein the method comprises the following steps: when the sinking point is selected on the image in the step 4), the deformation quantity obtained by the conventional InSAR cannot be completely used as the final mining area ground surface sinking quantity, wherein a large amount of phase unwrapping errors exist, and the point position of the edge of a partial sinking basin is selected as a reliable sinking point on the image according to coherence.

3. The method for monitoring and calculating interferometry of synthetic aperture radar for mining area surface mining subsidence according to claim 1, wherein the method comprises the following steps: in the step 5), when solving the probability integration method parameters, combining a small amount of measured data of the earth surface observation station with the point subsidence calculated by the InSAR technology, and continuously and circularly iterating and solving to calculate final probability integration method parameters by using a genetic optimization algorithm through steps of crossing, mutation and the like; providing a range of predicted parameters of a probability integration method, generating a group of initial chromosome population, and predicting the subsidence corresponding to the actual measurement point position of the earth surface; using the surface subsidence and a small number of surface actual measurement points resolved by InSAR as a reference, and comparing the actual measurement values with the predicted values to establish a fitness function; calculating the fitness of the fitter according to the survival principle of the fitter, namely according to a fitness function, and selecting a chromosome which is more suitable for the environment from the fitness function to carry out crossing and mutation so as to generate a new generation chromosome group which is more suitable for the environment, namely a probability integration method parameter; and continuously iterating and evolving, and finally converging and acquiring a chromosome which is most suitable for the environment, wherein the population is the optimal probability integration method parameter.

4. The method for monitoring and calculating interferometry of synthetic aperture radar for mining area surface mining subsidence according to claim 1, wherein the method comprises the following steps: when the surface subsidence basin is solved in the step 6), the subsidence value of any point (x, y) on the surface caused by underground mining is solved according to a probability integration method model by adopting a mining subsidence prediction theory which is most widely applied and perfect:

$W(x,y)=\mathrm{mq}\cos \mathrm{\&alpha;}\&Integral;\&Integral;W_{\mathrm{oci}}(x,y)\mathrm{dxdy}$

$W_{\mathrm{coi}}(x,y)=\frac{l}{r^2}\mathrm{epx}\left(\frac{-\mathrm{\&pi;}{(x-x_i)}^2}{r^2}\right)-\exp \left(\frac{-\mathrm{\&pi;}{(y-y_i+H_i\mathrm{ctg\&theta;})}^2}{r^2}\right)$

in the formula: q is a sinking coefficient; r is the major radius of influence, r = H₀/tgβ;H₀Average mining depth; theta is a mining influence propagation angle; tg β is the primary influence angle β tangent; (xi, yi) is a plane coordinate of the central point of the mining unit i; (x, y) is the coordinates of an arbitrary point on the earth's surface.

Images