CN112269176B - Early recognition and monitoring method for mine surface subsidence - Google Patents

Early recognition and monitoring method for mine surface subsidence Download PDF

Info

Publication number
CN112269176B
CN112269176B CN202011094157.7A CN202011094157A CN112269176B CN 112269176 B CN112269176 B CN 112269176B CN 202011094157 A CN202011094157 A CN 202011094157A CN 112269176 B CN112269176 B CN 112269176B
Authority
CN
China
Prior art keywords
monitoring
amplitude
point
phase
surface subsidence
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202011094157.7A
Other languages
Chinese (zh)
Other versions
CN112269176A (en
Inventor
王慧妮
李康伦
涂美义
王卓羽
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Wuhan Institute of Technology
Original Assignee
Wuhan Institute of Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Wuhan Institute of Technology filed Critical Wuhan Institute of Technology
Priority to CN202011094157.7A priority Critical patent/CN112269176B/en
Publication of CN112269176A publication Critical patent/CN112269176A/en
Application granted granted Critical
Publication of CN112269176B publication Critical patent/CN112269176B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
    • G01S13/9021—SAR image post-processing techniques
    • G01S13/9023—SAR image post-processing techniques combined with interferometric techniques
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
    • G01S13/9094—Theoretical aspects
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
    • G01S7/41—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00 using analysis of echo signal for target characterisation; Target signature; Target cross-section
    • G01S7/418—Theoretical aspects

Abstract

The invention discloses an early recognition and monitoring method for mine ground surface subsidence, which comprises the steps of firstly selecting PS by using an amplitude deviation method, then carrying out DS point recognition by combining classification information and a mathematical statistical method, connecting selected measuring points by using a Delaunay triangulation network and adding a distance limiting principle by combining the PS and the DS, removing an atmospheric phase by using one of an atmospheric model method, high-pass filtering of a time domain and low-pass filtering of a space domain to reduce an atmospheric phase APS effect, DEM residual errors and orbit errors, and finally estimating DEM errors and deformation rates; the deformation estimation technology based on the robustness of the PS and the DS can perform high-density early identification of the surface subsidence in a large area of the mine, can quickly and accurately acquire the surface subsidence information in a large range of the mine, improves the identification rate, the accuracy and the monitoring density of a subsidence area, and has important significance for the early identification and monitoring of the surface subsidence area of the mine.

Description

Early recognition and monitoring method for mine surface subsidence
Technical Field
The invention relates to the technical field of surface deformation identification and monitoring, in particular to an early identification and monitoring method for mine surface subsidence.
Background
Along with the rapid development of economy in China, the demand for exploiting underground energy is more and more large, the exploitation of underground coal usually causes continuous ground settlement, which can cause great threat to infrastructure and safety of coal mine areas, monitoring and investigation on the settlement of the coal mine areas are urgently needed, knowledge is provided for a disaster early warning system, and the traditional leveling and GPS data can generate reliable ground settlement measurement, however, the field investigation is time-consuming and labor-consuming, and a deformation map with high spatial sampling density cannot be provided;
synthetic Aperture Radar (SAR) interferometry (InSAR) is a powerful technique for obtaining high-precision elevation angle and surface deformation along the sight direction by using phase information of different SAR images, a DInSAR technique has been used for monitoring surface movement of a coal mine area, however, the conventional DInSAR is generally influenced by decorrelation of time, geometry and atmosphere, and due to the reasons of low density of PSs or coherent points, uneven distribution and the like, the conventional method has certain limitation when drawing a mining subsidence map, and most coal mine areas are located in suburban areas, artificial objects are too sparse to obtain enough PSI measurement points, so that the invention provides an early identification and monitoring method for mine surface subsidence to solve the problems in the prior art.
Disclosure of Invention
In view of the above problems, the present invention aims to provide an early identification and monitoring method for mine surface subsidence, which utilizes the technology of performing deformation estimation by using the robustness of PS and DS, can perform high-density early identification of surface subsidence in a mine in a large area, can quickly and accurately acquire surface subsidence information in a large area of the mine, and improves the identification rate, accuracy and monitoring density of a subsidence area, and has important significance for early identification and monitoring of the mine surface subsidence area.
In order to realize the purpose of the invention, the invention is realized by the following technical scheme: a method for early recognition and monitoring of mine surface subsidence comprises the following steps:
selecting and registering main images, selecting N SAR images according to a minimum sum criterion of three baselines, a comprehensive correlation model method or a minimum time baseline criterion, selecting the main images from the N SAR images, carrying out combined registration on the selected main images, and selecting amplitude/coherent images as a selected image of a PS point and a DS point identification image;
selecting PS points, measuring the phase noise level by utilizing time sequence amplitude information, particularly in a high signal-to-noise ratio pixel, utilizing pixel amplitude stability characteristics to replace the phase noise level to obtain a target point, and if Gaussian noise with standard deviation of sigma n exists in a real part and an imaginary part of the SAR complex image respectively, subjecting the amplitude A to screening distribution, namely:
Figure BDA0002723120100000021
in the formula (1), g is the energy reflected by the ground object target, and the signal-to-noise ratio is expressed as:
Figure BDA0002723120100000022
when in use
Figure BDA0002723120100000023
The following relationships exist:
Figure BDA0002723120100000024
Figure BDA0002723120100000025
in the formula sigmavIs the standard deviation of the phase, σAIs the standard deviation of the amplitude of the time sequence, σnRIs the standard deviation of the real part of the complex number, σnIIs the imaginary standard deviation of the complex number, DAIs an amplitude dispersion index.
And taking the standard deviation of the time sequence amplitude and the standard deviation of the phase as a stability analysis tool of the PS point to perform stability analysis on the detected PS point.
Step three, potential DS selection, namely calibrating the SAR image, averaging the calibrated and calibrated SAR image to obtain an average amplitude map, identifying a potential DS area by using the average amplitude map, determining a statistical homogeneity point for each pixel point by using a KS (K-link quality) test method, eliminating a related target and a continuous scattering target through classification information to obtain a rough DS candidate set, and then performing a refinement test according to the statistical characteristics of each candidate object in the candidate set;
step four, optimizing DS selection, and carrying out AD statistics A on the candidate set object subjected to the thinning test by utilizing AD visual inspection and the pair of pixels a and b2Is composed of
Figure BDA0002723120100000031
Wherein N represents the number of SAR images, Fa (x) and Fb (x) are empirical cumulative distribution functions of amplitudes x of pixels a and b, Fa, b (x) are empirical joint cumulative distribution functions of combined pixel values amplitude x, if AD statistic is less than threshold value, it is assumed that two pixels belong to the same uniform area to obtain adjacent points, and when the candidate object and the adjacent points are distributed in the same area, the candidate object is selected to obtain an optimized DS selection point set.
Step five, calculating deformation rate and DEM error, selecting selected PS points and DS points, connecting the selected pixel points to remove atmospheric phase to reduce atmospheric phase APS effect, using Delaunay triangulation network and adding distance limiting principle to facilitate connecting the selected measuring points, estimating triangulation network deformation parameters and elevation correction values by an indirect adjustment observation method, and performing deformation estimation and DEM error estimation by formula (5)
Figure BDA0002723120100000041
T in formula (5)k、
Figure BDA0002723120100000042
R and theta respectively represent a time base line, a normal base line, a tilting range and an incidence angle, lambda is a wavelength, delta v and delta epsilon are deformation rate change and elevation error change of adjacent pixels, and a final deformation rate and DEM error are calculated by an equation (6)
Figure BDA0002723120100000043
In equation (6), N is the number of interferograms, and k-phase is the phase difference between two adjacent points in the kth interferogram.
And sixthly, identifying a landslide area in a large range by setting the deformation to be larger than 3 cm/year, marking out a key area, and monitoring by combining rail lifting and rail lowering, so that the monitoring certainty reaches full-area coverage.
The further improvement lies in that: in the second step, a point selection method of amplitude deviation index threshold and intensity threshold in series is used for selecting the PS point, wherein the formula (3) represents the point selection method of the amplitude deviation index threshold.
The further improvement lies in that: and the KS test in the third step is based on a fitting goodness test method of general distribution, an empirical distribution function is adopted to determine whether two pixels have consistency, and the potential DS is selected by a method of combining classification information and a mathematical statistics method in the third step.
The further improvement lies in that: and in the third step, the statistical characteristics of the candidate objects in the candidate set are subjected to the refining test by adopting an adaptive sample selection strategy.
The further improvement lies in that: the neighboring points selected in step four constitute an estimation window, the size and shape of which may vary, which estimation window can be used both for estimation of complex coherence values and for adaptive filtering.
The further improvement lies in that: and in the fifth step, the atmospheric phase is removed by using one of an atmospheric model method, high-pass filtering in a time domain and low-pass filtering in a space domain to reduce atmospheric phase APS, DEM residual and orbit errors.
The further improvement lies in that: in the fifth step, formula (5) represents the phase difference between two pixels x and y adjacent to each edge of the network in the kth interference pattern.
The invention has the beneficial effects that: the technology for deformation estimation by using the robustness of the PS and the DS can perform high-density early identification of the ground surface subsidence in a large-area mine, can quickly and accurately acquire ground surface subsidence information in a large range of the mine, improves the identification rate, the accuracy and the monitoring density of a subsidence area, has important significance for the early identification and monitoring of the ground surface subsidence area of the mine, can detect candidate objects subjected to a refinement test by using an AD test, is more sensitive to the change of a distribution tail, can give an optimal detection rate, is convenient for connecting selected measuring points by using a Delaunay triangulation network and adding a distance limiting principle, and can effectively increase the effective monitoring range of a monitoring area by combined monitoring of ascending and descending tracks.
Drawings
FIG. 1 is a flow chart of the present invention.
Detailed Description
In order to further understand the present invention, the following detailed description will be made with reference to the following examples, which are only used for explaining the present invention and are not to be construed as limiting the scope of the present invention.
Example 1
According to fig. 1, the embodiment provides a method for early identification and monitoring of mine surface subsidence, which includes the following steps:
selecting and registering main images, selecting N SAR images according to a minimum sum criterion of three baselines, a comprehensive correlation model method or a minimum time baseline criterion, selecting the main images from the N SAR images, carrying out combined registration on the selected main images, and selecting amplitude/coherent images as a selected image of a PS point and a DS point identification image;
selecting PS points, measuring the phase noise level by utilizing time sequence amplitude information, particularly in high signal-to-noise ratio pixels, and acquiring target points by utilizing pixel amplitude stability characteristics to replace the phase noise level, wherein for SAR complex images, the real part and the imaginary part respectively have standard deviations of sigmanThe amplitude a follows a filter distribution, i.e.:
Figure BDA0002723120100000061
in the formula (1), g is the energy reflected by the ground object target, and the signal-to-noise ratio is expressed as:
Figure BDA0002723120100000062
when in use
Figure BDA0002723120100000063
The following relationships exist:
Figure BDA0002723120100000064
Figure BDA0002723120100000065
in the formula sigmavIs the standard deviation of the phase, σAIs the standard deviation of the amplitude of the time sequence, σnRIs the standard deviation of the real part of the complex number, σnIIs the imaginary standard deviation of the complex number, DAIs an amplitude dispersion index.
And taking the standard deviation of the time sequence amplitude and the standard deviation of the phase as a stability analysis tool of the PS point to perform stability analysis on the detected PS point.
Step three, potential DS selection, namely calibrating the SAR image, averaging the calibrated and calibrated SAR image to obtain an average amplitude map, identifying a potential DS area by using the average amplitude map, determining a statistical homogeneity point for each pixel point by using a KS (K-link quality) test method, eliminating a related target and a continuous scattering target through classification information to obtain a rough DS candidate set, and then performing a refinement test according to the statistical characteristics of each candidate object in the candidate set;
step four, optimizing DS selection, and carrying out AD statistics A on the candidate set object subjected to the thinning test by utilizing AD visual inspection and the pair of pixels a and b2Is composed of
Figure BDA0002723120100000071
Wherein N represents the number of SAR images, Fa (x) and Fb (x) are empirical cumulative distribution functions of amplitudes x of pixels a and b, Fa, b (x) are empirical joint cumulative distribution functions of combined pixel values amplitude x, if AD statistic is less than threshold value, it is assumed that two pixels belong to the same uniform area to obtain adjacent points, and when the candidate object and the adjacent points are distributed in the same area, the candidate object is selected to obtain an optimized DS selection point set.
Step five, calculating deformation rate and DEM error, selecting selected PS points and DS points, connecting the selected pixel points to remove atmospheric phase to reduce atmospheric phase APS effect, using Delaunay triangulation network and adding distance limiting principle to facilitate connecting the selected measuring points, estimating triangulation network deformation parameters and elevation correction values by an indirect adjustment observation method, and performing deformation estimation and DEM error estimation by formula (5)
Figure BDA0002723120100000072
T in formula (5)k、
Figure BDA0002723120100000073
R and theta respectively represent a time base line, a normal base line, a tilting range and an incidence angle, lambda is a wavelength, delta v and delta epsilon are deformation rate change and elevation error change of adjacent pixels, and a final deformation rate and DEM error are calculated by an equation (6)
Figure BDA0002723120100000074
In equation (6), N is the number of interferograms, and k-phase is the phase difference between two adjacent points in the kth interferogram.
And sixthly, identifying a landslide area in a large range by setting the deformation to be larger than 3 cm/year, marking out a key area, and monitoring by combining rail lifting and rail lowering, so that the monitoring certainty reaches full-area coverage.
In the second step, a point selection method of amplitude deviation index threshold and intensity threshold in series is used for selecting the PS point, wherein the formula (3) represents the point selection method of the amplitude deviation index threshold.
And the KS test in the third step is based on a fitting goodness test method of general distribution, an empirical distribution function is adopted to determine whether two pixels have consistency, and the potential DS is selected by a method of combining classification information and a mathematical statistics method in the third step.
And in the third step, the statistical characteristics of the candidate objects in the candidate set are subjected to the refining test by adopting an adaptive sample selection strategy.
The neighboring points selected in step four constitute an estimation window, the size and shape of which may vary, which estimation window can be used both for estimation of complex coherence values and for adaptive filtering.
And in the fifth step, an atmospheric model method is used for removing the atmospheric phase so as to reduce the atmospheric phase APS, the DEM residual error and the orbit error.
In the fifth step, formula (5) represents the phase difference between two pixels x and y adjacent to each edge of the network in the kth interference pattern.
Example 2
According to fig. 1, the embodiment provides a method for early identification and monitoring of mine surface subsidence, which includes the following steps:
selecting and registering main images, selecting N SAR images according to a minimum sum criterion of three baselines, a comprehensive correlation model method or a minimum time baseline criterion, selecting the main images from the N SAR images, carrying out combined registration on the selected main images, and selecting amplitude/coherent images as a selected image of a PS point and a DS point identification image;
selecting PS points, measuring the phase noise level by utilizing time sequence amplitude information, particularly in high signal-to-noise ratio pixels, and acquiring target points by utilizing pixel amplitude stability characteristics to replace the phase noise level, wherein for SAR complex images, the real part and the imaginary part respectively have standard deviations of sigmanThe amplitude a follows a filter distribution, i.e.:
Figure BDA0002723120100000091
in the formula (1), g is the energy reflected by the ground object target, and the signal-to-noise ratio is expressed as:
Figure BDA0002723120100000092
when in use
Figure BDA0002723120100000093
The following relationships exist:
Figure BDA0002723120100000094
Figure BDA0002723120100000095
in the formula sigmavIs the standard deviation of the phase, σAIs the standard deviation of the amplitude of the time sequence, σnRIs the standard deviation of the real part of the complex number, σnIIs the imaginary standard deviation of the complex number, DAIs an amplitude dispersion index.
And taking the standard deviation of the time sequence amplitude and the standard deviation of the phase as a stability analysis tool of the PS point to perform stability analysis on the detected PS point.
Step three, potential DS selection, namely calibrating the SAR image, averaging the calibrated and calibrated SAR image to obtain an average amplitude map, identifying a potential DS area by using the average amplitude map, determining a statistical homogeneity point for each pixel point by using a KS (K-link quality) test method, eliminating a related target and a continuous scattering target through classification information to obtain a rough DS candidate set, and then performing a refinement test according to the statistical characteristics of each candidate object in the candidate set;
step four, optimizing DS selection, and carrying out AD statistics A on the candidate set object subjected to the thinning test by utilizing AD visual inspection and the pair of pixels a and b2Is composed of
Figure BDA0002723120100000101
Wherein N represents the number of SAR images, Fa (x) and Fb (x) are empirical cumulative distribution functions of amplitudes x of pixels a and b, Fa, b (x) are empirical joint cumulative distribution functions of combined pixel values amplitude x, if AD statistic is less than threshold value, it is assumed that two pixels belong to the same uniform area to obtain adjacent points, and when the candidate object and the adjacent points are distributed in the same area, the candidate object is selected to obtain an optimized DS selection point set.
Step five, calculating deformation rate and DEM error, selecting selected PS points and DS points, connecting the selected pixel points to remove atmospheric phase to reduce atmospheric phase APS effect, using Delaunay triangulation network and adding distance limiting principle to facilitate connecting the selected measuring points, estimating triangulation network deformation parameters and elevation correction values by an indirect adjustment observation method, and performing deformation estimation and DEM error estimation by formula (5)
Figure BDA0002723120100000102
T in formula (5)k、
Figure BDA0002723120100000103
R and theta respectively represent a time base line, a normal base line, a tilting range and an incidence angle, lambda is a wavelength, delta v and delta epsilon are deformation rate change and elevation error change of adjacent pixels, and a final deformation rate and DEM error are calculated by an equation (6)
Figure BDA0002723120100000104
In equation (6), N is the number of interferograms, and k-phase is the phase difference between two adjacent points in the kth interferogram.
And sixthly, identifying a landslide area in a large range by setting the deformation to be larger than 3 cm/year, marking out a key area, and monitoring by combining rail lifting and rail lowering, so that the monitoring certainty reaches full-area coverage.
In the second step, a point selection method of amplitude deviation index threshold and intensity threshold in series is used for selecting the PS point, wherein the formula (3) represents the point selection method of the amplitude deviation index threshold.
And the KS test in the third step is based on a fitting goodness test method of general distribution, an empirical distribution function is adopted to determine whether two pixels have consistency, and the potential DS is selected by a method of combining classification information and a mathematical statistics method in the third step.
And in the third step, the statistical characteristics of the candidate objects in the candidate set are subjected to the refining test by adopting an adaptive sample selection strategy.
The neighboring points selected in step four constitute an estimation window, the size and shape of which may vary, which estimation window can be used both for estimation of complex coherence values and for adaptive filtering.
And fifthly, removing the atmospheric phase by using a time domain high-pass filtering method to reduce the atmospheric phase APS, the DEM residual error and the orbit error.
In the fifth step, formula (5) represents the phase difference between two pixels x and y adjacent to each edge of the network in the kth interference pattern.
Example 3
According to fig. 1, the embodiment provides a method for early identification and monitoring of mine surface subsidence, which includes the following steps:
selecting and registering main images, selecting N SAR images according to a minimum sum criterion of three baselines, a comprehensive correlation model method or a minimum time baseline criterion, selecting the main images from the N SAR images, carrying out combined registration on the selected main images, and selecting amplitude/coherent images as a selected image of a PS point and a DS point identification image;
selecting PS points, measuring the phase noise level by utilizing time sequence amplitude information, particularly in high signal-to-noise ratio pixels, and acquiring target points by utilizing pixel amplitude stability characteristics to replace the phase noise level, wherein for SAR complex images, the real part and the imaginary part respectively have standard deviations of sigmanThe amplitude a follows a filter distribution, i.e.:
Figure BDA0002723120100000121
in the formula (1), g is the energy reflected by the ground object target, and the signal-to-noise ratio is expressed as:
Figure BDA0002723120100000122
when in use
Figure BDA0002723120100000123
The following relationships exist:
Figure BDA0002723120100000125
Figure BDA0002723120100000124
in the formula sigmavIs the standard deviation of the phase, σAIs the standard deviation of the amplitude of the time sequence, σnRIs the standard deviation of the real part of the complex number, σnIIs the imaginary standard deviation of the complex number, DAIs an amplitude dispersion index.
And taking the standard deviation of the time sequence amplitude and the standard deviation of the phase as a stability analysis tool of the PS point to perform stability analysis on the detected PS point.
Step three, potential DS selection, namely calibrating the SAR image, averaging the calibrated and calibrated SAR image to obtain an average amplitude map, identifying a potential DS area by using the average amplitude map, determining a statistical homogeneity point for each pixel point by using a KS (K-link quality) test method, eliminating a related target and a continuous scattering target through classification information to obtain a rough DS candidate set, and then performing a refinement test according to the statistical characteristics of each candidate object in the candidate set;
step four, optimizing DS selection, and carrying out AD statistics A on the candidate set object subjected to the thinning test by utilizing AD visual inspection and the pair of pixels a and b2Is composed of
Figure BDA0002723120100000131
Wherein N represents the number of SAR images, Fa (x) and Fb (x) are empirical cumulative distribution functions of amplitudes x of pixels a and b, Fa, b (x) are empirical joint cumulative distribution functions of combined pixel values amplitude x, if AD statistic is less than threshold value, it is assumed that two pixels belong to the same uniform area to obtain adjacent points, and when the candidate object and the adjacent points are distributed in the same area, the candidate object is selected to obtain an optimized DS selection point set.
Step five, calculating deformation rate and DEM error, selecting selected PS points and DS points, connecting the selected pixel points to remove atmospheric phase to reduce atmospheric phase APS effect, using Delaunay triangulation network and adding distance limiting principle to facilitate connecting the selected measuring points, estimating triangulation network deformation parameters and elevation correction values by an indirect adjustment observation method, and performing deformation estimation and DEM error estimation by formula (5)
Figure BDA0002723120100000132
T in formula (5)k、
Figure BDA0002723120100000133
R and theta respectively represent a time base line, a normal base line, a tilting range and an incidence angle, lambda is a wavelength, delta v and delta epsilon are deformation rate change and elevation error change of adjacent pixels, and a final deformation rate and DEM error are calculated by an equation (6)
Figure BDA0002723120100000134
In equation (6), N is the number of interferograms, and k-phase is the phase difference between two adjacent points in the kth interferogram.
And sixthly, identifying a landslide area in a large range by setting the deformation to be larger than 3 cm/year, marking out a key area, and monitoring by combining rail lifting and rail lowering, so that the monitoring certainty reaches full-area coverage.
In the second step, a point selection method of amplitude deviation index threshold and intensity threshold in series is used for selecting the PS point, wherein the formula (3) represents the point selection method of the amplitude deviation index threshold.
And the KS test in the third step is based on a fitting goodness test method of general distribution, an empirical distribution function is adopted to determine whether two pixels have consistency, and the potential DS is selected by a method of combining classification information and a mathematical statistics method in the third step.
And in the third step, the statistical characteristics of the candidate objects in the candidate set are subjected to the refining test by adopting an adaptive sample selection strategy.
The neighboring points selected in step four constitute an estimation window, the size and shape of which may vary, which estimation window can be used both for estimation of complex coherence values and for adaptive filtering.
And fifthly, removing the atmospheric phase by using a low-pass filtering method of a spatial domain to reduce the atmospheric phase APS, the DEM residual error and the orbit error.
In the fifth step, formula (5) represents the phase difference between two pixels x and y adjacent to each edge of the network in the kth interference pattern.
The technology for deformation estimation by using the robustness of the PS and the DS, which is provided by the mine surface subsidence early-stage identification and monitoring method, can carry out high-density surface subsidence early-stage identification in a large-area mine, can quickly and accurately acquire surface subsidence information in a large range of the mine, improves the identification rate, the accuracy and the monitoring density of a subsidence area, has important significance for the early identification and monitoring of the mine surface subsidence area, can detect candidate objects after a thinning test by using an AD test, is more sensitive to the change of a distribution tail, can give the optimal detection rate, is convenient for connecting selected measuring points by using a Delaunay triangulation network and adding a distance limiting principle, and can effectively increase the effective monitoring range of the monitoring area by combined monitoring of ascending and descending tracks.
The foregoing illustrates and describes the principles, general features, and advantages of the present invention. It will be understood by those skilled in the art that the present invention is not limited to the embodiments described above, which are described in the specification and illustrated only to illustrate the principle of the present invention, but that various changes and modifications may be made therein without departing from the spirit and scope of the present invention, which fall within the scope of the invention as claimed. The scope of the invention is defined by the appended claims and equivalents thereof.

Claims (7)

1. The early-stage identification and monitoring method for the mine surface subsidence is characterized by comprising the following steps of:
selecting and registering main images, selecting N SAR images according to a minimum sum criterion of three baselines, a comprehensive correlation model method or a minimum time baseline criterion, selecting the main images from the N SAR images, carrying out combined registration on the selected main images, and selecting amplitude/coherent images as a selected image of a PS point and a DS point identification image;
selecting PS points, measuring the phase noise level by utilizing time sequence amplitude information, and in the high signal-to-noise ratio pixel, utilizing the pixel amplitude stability characteristic to replace the phase noise level to obtain a target point, wherein for the SAR complex image, the real part and the imaginary part respectively have standard deviation sigmanThe amplitude a follows a filter distribution, i.e.:
Figure FDA0003164670540000011
in the formula (1), g is the energy reflected by the ground object target, and the signal-to-noise ratio is expressed as:
Figure FDA0003164670540000012
when in use
Figure FDA0003164670540000013
The following relationships exist:
Figure FDA0003164670540000014
Figure FDA0003164670540000015
in the formula sigmavIs the standard deviation of the phase, σAIs the standard deviation of the amplitude of the time sequence, σnRIs the standard deviation of the real part of the complex number, σnIIs the imaginary standard deviation of the complex number, DAIs an amplitude dispersion index;
taking the standard deviation of the time sequence amplitude and the standard deviation of the phase as a stability analysis tool of the PS point to perform stability analysis on the detected PS point;
step three, potential DS selection, namely calibrating the SAR image, averaging the calibrated and calibrated SAR image to obtain an average amplitude map, identifying a potential DS area by using the average amplitude map, determining a statistical homogeneity point for each pixel point by using a KS (K-link quality) test method, eliminating a related target and a continuous scattering target through classification information to obtain a rough DS candidate set, and then performing a refinement test according to the statistical characteristics of each candidate object in the candidate set;
step four, optimizing DS selection, and carrying out AD statistics A on the candidate set object subjected to the thinning test by utilizing AD visual inspection and the pair of pixels a and b2Is composed of
Figure FDA0003164670540000021
Wherein N represents the number of SAR images, Fa (x) and Fb (x) are empirical cumulative distribution functions of amplitudes x of pixels a and b, Fa, b (x) are empirical joint cumulative distribution functions of combined pixel values and amplitudes x, if AD statistic is less than threshold, the two pixels are assumed to belong to the same uniform area to obtain adjacent points, and if the candidate object and the adjacent points are distributed in the same area, the candidate object is selected to obtain an optimized DS selection point set;
step five, calculating deformation rate and DEM error, selecting selected PS points and DS points, connecting the selected pixel points to remove atmospheric phase to reduce atmospheric phase APS effect, using Delaunay triangulation network and adding distance limiting principle to facilitate connecting selected measuring points, estimating triangulation network deformation parameters and elevation correction values by an indirect adjustment observation method, and performing DEM error estimation by formula (5)
Figure FDA0003164670540000022
T in formula (5)k、
Figure FDA0003164670540000023
R and theta respectively represent a time base line, a normal base line, a tilting range and an incidence angle, lambda is a wavelength, delta v and delta epsilon are deformation rate change and elevation error change of adjacent pixels, and a final deformation rate (6) is calculated
Figure FDA0003164670540000024
N in the formula (6) is the number of interferograms, and the k phase is the phase difference of two adjacent points in the kth-channel interferogram;
and sixthly, identifying a landslide area in a large range by setting the deformation to be larger than 3 cm/year, marking out a key area, and monitoring by combining rail lifting and rail lowering, so that the monitoring certainty reaches full-area coverage.
2. The method for early identification and monitoring of mine surface subsidence according to claim 1, wherein: in the second step, a point selection method of amplitude deviation index threshold and intensity threshold in series is used for selecting the PS point, wherein the formula (3) represents the point selection method of the amplitude deviation index threshold.
3. The method for early identification and monitoring of mine surface subsidence according to claim 1, wherein: and the KS test in the third step is based on a fitting goodness test method of general distribution, an empirical distribution function is adopted to determine whether two pixels have consistency, and the potential DS is selected by a method of combining classification information and a mathematical statistics method in the third step.
4. The method for early identification and monitoring of mine surface subsidence according to claim 1, wherein: and in the third step, the statistical characteristics of the candidate objects in the candidate set are subjected to the refining test by adopting an adaptive sample selection strategy.
5. The method for early identification and monitoring of mine surface subsidence according to claim 1, wherein: the neighboring points selected in step four constitute an estimation window, the size and shape of which may vary, which estimation window can be used both for estimation of complex coherence values and for adaptive filtering.
6. The method for early identification and monitoring of mine surface subsidence according to claim 1, wherein: and in the fifth step, the atmospheric phase is removed by using one of an atmospheric model method, high-pass filtering in a time domain and low-pass filtering in a space domain to reduce atmospheric phase APS, DEM residual and orbit errors.
7. The method for early identification and monitoring of mine surface subsidence according to claim 1, wherein: in the fifth step, formula (5) represents the phase difference between two pixels x and y adjacent to each edge of the network in the kth interference pattern.
CN202011094157.7A 2020-10-14 2020-10-14 Early recognition and monitoring method for mine surface subsidence Active CN112269176B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011094157.7A CN112269176B (en) 2020-10-14 2020-10-14 Early recognition and monitoring method for mine surface subsidence

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011094157.7A CN112269176B (en) 2020-10-14 2020-10-14 Early recognition and monitoring method for mine surface subsidence

Publications (2)

Publication Number Publication Date
CN112269176A CN112269176A (en) 2021-01-26
CN112269176B true CN112269176B (en) 2021-09-14

Family

ID=74338184

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011094157.7A Active CN112269176B (en) 2020-10-14 2020-10-14 Early recognition and monitoring method for mine surface subsidence

Country Status (1)

Country Link
CN (1) CN112269176B (en)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113253271A (en) * 2021-07-12 2021-08-13 中国测绘科学研究院 Deformation rate estimation method and system
CN116148853B (en) * 2023-02-20 2024-02-20 中国矿业大学 Mining area subsidence interference phase filtering method and device

Citations (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104111457A (en) * 2014-07-23 2014-10-22 中国国土资源航空物探遥感中心 Mutual-inspection and temporal fusion method for surface subsidence monitoring result of PSInSAR for lifting track
CN105487072A (en) * 2015-12-29 2016-04-13 武汉工程大学 Method and system of joint location based on T2/R time difference and Doppler shift
CN106023157A (en) * 2016-05-10 2016-10-12 电子科技大学 Mountain area surface micro deformation information extraction method based on SAR images
WO2018027332A1 (en) * 2016-08-08 2018-02-15 Comercial E Industrial Gesecology Limitada Method and system for the analysis and generation of early or predictive alerts concerning the stability of slopes in open-pit mines
WO2018229485A1 (en) * 2017-06-15 2018-12-20 The University Of Nottingham Land deformation measurement
CN109118520A (en) * 2018-07-27 2019-01-01 中国电力科学研究院有限公司 A kind of minery power grid shaft tower displacement monitoring method and system
CN110412574A (en) * 2019-09-05 2019-11-05 河海大学 A kind of distributed object InSAR timing sequence process method and apparatus of temporal and spatial coherence enhancing
CN110456346A (en) * 2019-06-28 2019-11-15 深圳市水务规划设计院股份有限公司 A kind of electric power pylon inclination monitoring method based on InSAR technology
CN110568440A (en) * 2019-09-10 2019-12-13 四川省地质工程勘察院集团有限公司 method for monitoring deformation of complex mountain area based on DS-InSAR technology
CN110763187A (en) * 2019-09-30 2020-02-07 中国科学院测量与地球物理研究所 Stable ground settlement monitoring method based on radar distributed target
CN111142119A (en) * 2020-01-10 2020-05-12 中国地质大学（北京） Mine geological disaster dynamic identification and monitoring method based on multi-source remote sensing data

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8718940B2 (en) * 2010-11-30 2014-05-06 Halliburton Energy Services, Inc. Evaluating surface data

Patent Citations (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104111457A (en) * 2014-07-23 2014-10-22 中国国土资源航空物探遥感中心 Mutual-inspection and temporal fusion method for surface subsidence monitoring result of PSInSAR for lifting track
CN105487072A (en) * 2015-12-29 2016-04-13 武汉工程大学 Method and system of joint location based on T2/R time difference and Doppler shift
CN106023157A (en) * 2016-05-10 2016-10-12 电子科技大学 Mountain area surface micro deformation information extraction method based on SAR images
WO2018027332A1 (en) * 2016-08-08 2018-02-15 Comercial E Industrial Gesecology Limitada Method and system for the analysis and generation of early or predictive alerts concerning the stability of slopes in open-pit mines
WO2018229485A1 (en) * 2017-06-15 2018-12-20 The University Of Nottingham Land deformation measurement
CN109118520A (en) * 2018-07-27 2019-01-01 中国电力科学研究院有限公司 A kind of minery power grid shaft tower displacement monitoring method and system
CN110456346A (en) * 2019-06-28 2019-11-15 深圳市水务规划设计院股份有限公司 A kind of electric power pylon inclination monitoring method based on InSAR technology
CN110412574A (en) * 2019-09-05 2019-11-05 河海大学 A kind of distributed object InSAR timing sequence process method and apparatus of temporal and spatial coherence enhancing
CN110568440A (en) * 2019-09-10 2019-12-13 四川省地质工程勘察院集团有限公司 method for monitoring deformation of complex mountain area based on DS-InSAR technology
CN110763187A (en) * 2019-09-30 2020-02-07 中国科学院测量与地球物理研究所 Stable ground settlement monitoring method based on radar distributed target
CN111142119A (en) * 2020-01-10 2020-05-12 中国地质大学（北京） Mine geological disaster dynamic identification and monitoring method based on multi-source remote sensing data

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
"A new algorithm for processing interferometric Data-stacks: SqueeSAR";Alessandro Ferretti等;《IEEE transactions on geosciences and remote sensing》;20110930;第49卷(第9期);第3460-3470页 *

Also Published As

Publication number Publication date
CN112269176A (en) 2021-01-26

Similar Documents

Publication Publication Date Title
Pacheco et al. Retrieval of nearshore bathymetry from Landsat 8 images: A tool for coastal monitoring in shallow waters
CN112269176B (en) Early recognition and monitoring method for mine surface subsidence
Wang et al. Retrieval of phase history parameters from distributed scatterers in urban areas using very high resolution SAR data
CN110174044B (en) Bridge longitudinal displacement deformation monitoring method based on PSI technology
CN112198511A (en) Integrated geological disaster census method based on starry sky and ground
JP2003500658A (en) Procedures for radar measurements of urban area and landslide motion.
CN113960595A (en) Surface deformation monitoring method and system
CN110473187B (en) Object-oriented line scanning three-dimensional pavement crack extraction method
CN112284332B (en) High-rise building settlement monitoring result three-dimensional positioning method based on high-resolution INSAR
CN110763187A (en) Stable ground settlement monitoring method based on radar distributed target
CN116047519A (en) Point selection method based on synthetic aperture radar interferometry technology
CN113281749B (en) Timing sequence InSAR high coherence point selection method considering homogeneity
Zhao et al. A temporal phase coherence estimation algorithm and its application on DInSAR pixel selection
Cenci et al. Describing the quality assessment workflow designed for DEM products distributed via the Copernicus Programme. Case study: The absolute vertical accuracy of the Copernicus DEM dataset in Spain
CN113687353A (en) DS target phase optimization method based on homogeneous pixel time sequence phase matrix decomposition
Tao et al. Ground deformation retrieval using quasi coherent targets DInSAR, with application to suburban area of Tianjin, China
CN115236658B (en) Road surface crack three-dimensional form monitoring method based on active radar remote sensing cooperation
Korzeniowska et al. Mapping gullies using terrain surface roughness
CN113625241A (en) Differential settlement monitoring and early warning method
Narayan et al. Persistent scatter identification and look-angle error estimation using similar time-series interferometric pixels
Yan et al. Determining suitable spaceborne SAR observations and ground control points for surface deformation study in rugged terrain with InSAR technique
CN113758445A (en) Detection method for tire tread depth and embedded foreign matter fragments
Friedl et al. Using SAR Interferograms and Coherence Images for Object-Based Delineation of Unstable Slopes
CN112859073B (en) Road damage assessment method based on PSInSAR technology
He et al. An improved method for phase triangulation algorithm based on the coherence matrix eigen-decomposition in time-series SAR interferometry

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant