Method for obtaining high-precision ground subsidence basin by fusing unmanned aerial vehicle DEM and InSAR phase

Abstract

The application discloses a method for obtaining high-precision surface subsidence basins by fusing unmanned aerial vehicle DEM and InSAR phases, and belongs to the field of surface subsidence extraction and geological disaster monitoring. The method uses a deformation monitoring method for fusing DSM (DEM) data and InSAR data obtained by an unmanned aerial vehicle, combines the large deformation of the unmanned aerial vehicle DSM (DEM) data and the high-precision characteristics of the InSAR data, and realizes the purpose of quickly, accurately and comprehensively obtaining deformation information of a mining area. The following steps are performed: S1, data acquisition; S2, data conversion; S3, phase fusion; S4, data solving; and finally, complete and accurate deformation of a subsidence basin in the mining area is obtained. The phase "whole number" in the image processing process is corrected by using the unmanned aerial vehicle to obtain the subsidence basin information in the mining area, so that the D-InSAR technology can be used to monitor the subsidence level of the mining area.

Description

TECHNICAL FIELD

[0001] The application relates to a method for obtaining a high-precision surface subsidence basin by fusing a UAV DEM and an InSAR phase, and belongs to the field of surface subsidence extraction and geological disaster monitoring. BACKGROUND

[0002] Coal mining not only meets the national energy demand, but also causes serious geological environment damage problems in the mining area, such as surface subsidence, ground fissures and other geological disasters. The rock stratum and surface subsidence caused by coal mining are the root cause of the mine geological environment damage. In particular, large-scale and high-intensity rapid mining is mainly used in the western region of China, and the surface movement damage caused by mining presents the characteristics of fast deformation speed, deep damage degree and wide range of influence. Under this form, how to quickly, accurately and comprehensively monitor the surface subsidence and environmental damage caused by high-intensity mining in the western region is the key to solving the problem.

[0003] According to the station arrangement form, the conventional surface movement deformation monitoring method can be divided into profile line observation stations and network observation stations. The former is widely used, which is to arrange "points" into "lines" on the working face strike / dip main section to form observation stations; the latter is to arrange some measurement points to form grid-shaped observation stations, and the latter is less used due to the limitation of terrain and object conditions (such as occupying more land). Both kinds of observation stations obtain limited "point" data. The conventional observation station mainly adopts conventional measurement methods including triangulation, traverse measurement (theodolite, total station), leveling (level), GPS and the like; the conventional observation station has high precision, is the most commonly used and effective means for obtaining the mine surface deformation at present, and plays an important role in the mining subsidence research and actual engineering. With the development of research, some shortcomings appear: 1) small scale, generally taking a single working face as the monitoring object, and the measurement method also determines that it is not convenient for large-scale (large-scale) operation; 2) only the surface movement deformation value is taken as the monitoring content, and the ecological environment damage information (such as vegetation spectrum, water body and the like) of the mining area cannot be obtained; 3) high cost, long cycle (poor timeliness), large workload, need to bury measurement points and difficult to be preserved for a long time and the like; 4) only "point" observation, the data and information amount is small, and the subsidence basin characteristics cannot be completely reflected. This also leads to that it has been unable to adapt to the surface damage monitoring task of high-intensity and large-scale mining in the western region under the new situation.

[0004] Unmanned Aerial Vehicle (UAV) photogrammetry and InSAR are currently widely used surveying and mapping new technologies. UAV photogrammetry is a platform based on unmanned aerial vehicles, integrating high-resolution (multi / hyperspectral) low-altitude optical sensors, GPS and IMU technologies, which can simultaneously obtain surface subsidence deformation and environmental information, adapt to multi-scale monitoring tasks, and have the advantages of flexibility, high efficiency, precision, low operation cost (outstanding timeliness and cost performance). InSAR technology obtains ground deformation through interference phase information, and its principle is to calculate the phase difference of SAR images of the same area through two or more passes to obtain DEM or separate ground deformation information, with an accuracy of centimeters to millimeters.

[0005] UAV and InSAR technology also have their own shortcomings in the application of mining deformation monitoring. The platform of UAV is relatively light, the photographic posture is unstable, and the image distortion is large, resulting in relatively low elevation accuracy. The actual measurement shows that the elevation accuracy is centimeters to decimeters. UAV has an advantage in monitoring large deformation in mining areas (such as the center of the subsidence basin), but it cannot accurately monitor the boundary area of the subsidence basin. The advantage of InSAR is that it can accurately monitor the deformation of the boundary area of the subsidence basin, but it is easily affected by the space-time incoherence, which makes it unable to reliably obtain the large deformation in the center of the subsidence basin. Therefore, UAV photogrammetry and InSAR are complementary, and through the fusion of the two, complete subsidence basin data can be obtained, providing a new technical approach to solve the problem of monitoring the geological environment damage caused by large-scale and high-intensity mining in mines.

[0006] Therefore, how to combine UAV and InSAR to quickly, accurately and comprehensively monitor the surface subsidence and environmental damage caused by high-intensity mining in the west is a technical problem that needs to be solved by those skilled in the art. SUMMARY

[0007] The present application is directed to the above problems, and proposes a method for obtaining high-precision surface subsidence basin by fusing UAV DEM and InSAR phase. The deformation monitoring method uses fused UAV obtained DSM (DEM) data and InSAR data, combines the characteristics of large deformation of UAV DSM (DEM) data and high precision of InSAR data, and realizes the purpose of quickly, accurately and comprehensively obtaining deformation information of mining areas.

[0008] The technical scheme of the present application is as follows:

[0009] S1, data acquisition;

[0010] N scenes of SAR images in the target time period are acquired by InSAR, and two DEM data in the target time period are acquired by UAV aerial photography in succession;

[0011] S2, data conversion;

[0012] The N scenes of SAR images obtained in step S1 are subjected to interference and unwrapping processing to obtain N-1 unwrapped phases, and through "stacking" and re-wrapping processing, a complex wrapped phase is obtained;

[0013] The two DEM data obtained in step S1 are subtracted to obtain the large deformation of the mining area subsidence; then the large deformation of the mining area subsidence is converted into absolute interference phase according to the radar wavelength, and is converted into wrapped phase to obtain the wrapped phase containing large deformation;

[0014] S3, phase fusion;

[0015] The wrapped phase containing large deformation obtained in step S2 is fused with the complex wrapped phase to correct the whole number in the SAR phase, and a fused phase containing large deformation is obtained;

[0016] S4, data solving;

[0017] Through solving processing, a complete and accurate deformation of the mining area subsidence basin is finally obtained.

[0018] Further, in the process of correcting the whole number in the SAR phase in step S3;

[0019] 1) Complex wrapped phase, after multi-period accumulation and de-trending of the difference phase, the phase wrapping operator is used for wrapping operation to obtain the wrapped phase;

[0020] 2) The wrapped phase of large deformation, first, the surface deformation W(x, y) obtained by the unmanned aerial vehicle, the satellite wavelength λ and the viewing angle θ are used to calculate the absolute phase Φ containing large deformation according to the relationship expression between the subsidence value of a certain point on the ground and the radar satellite wavelength, the incident angle and the phase of the point in the basic principle of D-InSAR, and then the wrapped operator is used for wrapping operation on the absolute phase to obtain the wrapped phase containing large deformation;

[0021] According to the basic principle of D-InSAR, the relationship between the subsidence value of a certain point on the ground and the radar satellite wavelength, the incident angle and the phase of the point can be expressed as:

[0022]

[0023] In the formula, △W is the subsidence value of the ground point; φ is the phase of the ground point; λ is the satellite wavelength; θ is the viewing angle.

[0024] The absolute phase of the ground point φ is between (-π, π), so the calculation formula of the phase integer number Zφ at the ground point (x, y) can be determined according to the subsidence obtained by the unmanned aerial vehicle:

[0025]

[0026] In the formula, Int[·] represents the integer operation, and the other symbols have the same meaning as formula (1). W(x, y) is the subsidence obtained by the unmanned aerial vehicle at the ground point (x, y);

[0027] 3) The absolute interference phase is the real phase value. After phase wrapping, the absolute interference phase mainly includes the wrapped phase with the principal value in (-π, π] and the integer number 2kπ.

[0028] Further, in the step S2 of processing the N SAR images, the SAR images are arranged in chronological order. Two-track differential interference processing is first performed on the adjacent two images, and then phase unwrapping is performed. The regions with a coherence less than 0.2 are masked to obtain N-1 unwrapped phase images. Thereafter, the phase values are interpolated in the masked regions of the unwrapped phase, and the N-1 unwrapped phase images are accumulated.

[0029] The stable regions in the accumulated unwrapped phase are selected as control points, a phase trend surface is fitted and removed, that is, phase tilt processing is performed to eliminate the slope effect caused by inaccurate baseline. The unwrapped phase is converted into complex wrapped phase again according to the phase cosine as the real part and the phase sine as the imaginary part.

[0030] Further, the specific process of the phase fusion in the step S3 is as follows: the wrapped phase containing large deformation obtained in the step S2 is deducted from the complex wrapped phase, and only the small deformation is retained to obtain the residual small deformation unwrapped phase. The absolute interference phase after the conversion of the DEM data in the step S2 is added back to the residual small deformation unwrapped phase.

[0031] Affected by the characteristics of large subsidence level settlement in the mining area, the phase "integer number" cannot be accurately determined in the unwrapping process, which causes a large deviation in the monitoring result and becomes the most important factor restricting the application of the D-InSAR technology in the mining area subsidence monitoring. The unmanned aerial vehicle can obtain large deformation in the mining area. The "integer number" of the phase in the image processing process is corrected by using the subsidence basin information obtained by the unmanned aerial vehicle in the mining area, so that the D-InSAR technology for monitoring the subsidence level in the mining area can be improved. BRIEF DESCRIPTION OF DRAWINGS

[0032] Figure 1 is the work flow chart of the case;

[0033] Figure 2 is the specific work flow chart of the case;

[0034] Figure 3a is DEM data acquired by UAV photogrammetry on June 9, 2018 in the embodiment of the present case;

[0035] Figure 3b is DEM data acquired by UAV photogrammetry on April 16, 2019 in the embodiment of the present case;

[0036] Figure 4 is the region multi-view image after the SAR image of the monitoring region is cropped in the embodiment of the present case;

[0037] Figure 5a is the differential interferogram after differential interferometric processing of SAR images on June 9, 2018 and July 27, 2018 in the embodiment of the present case;

[0038] Figure 5b is the differential interferogram after differential interferometric processing of SAR images on July 27, 2018 and August 20, 2018 in the embodiment of the present case;

[0039] Figure 5c is the differential interferogram after differential interferometric processing of SAR images on August 20, 2018 and November 24, 2018 in the embodiment of the present case;

[0040] Figure 5d is the differential interferogram after differential interferometric processing of SAR images on November 24, 2018 and January 11, 2019 in the embodiment of the present case;

[0041] Figure 5e is the differential interferogram after differential interferometric processing of SAR images on January 11, 2019 and February 4, 2019 in the embodiment of the present case;

[0042] Figure 5f is the differential interferogram after differential interferometric processing of SAR images on February 4, 2019 and February 28, 2019 in the embodiment of the present case;

[0043] Figure 5g is the differential interferogram after differential interferometric processing of SAR images on February 28, 2019 and March 24, 2019 in the embodiment of the present case;

[0044] Figure 5h is the differential interferogram after differential interferometric processing of SAR images on March 24, 2019 and April 17, 2019 in the embodiment of the present case;

[0045] Figure 6a is the unwrapped phase image after phase unwrapping of the differential interferogram on June 9, 2018 and July 27, 2018 in the embodiment of the present case;

[0046] Figure 6bis the unwrapped phase map after phase unwrapping of the difference interferograms in the embodiment of the present case on July 27, 2018 and August 20, 2018;

[0047] Figure 6c is the unwrapped phase map after phase unwrapping of the difference interferograms in the embodiment of the present case on August 20, 2018 and November 24, 2018;

[0048] Figure 6d is the unwrapped phase map after phase unwrapping of the difference interferograms in the embodiment of the present case on November 24, 2018 and January 11, 2019;

[0049] Figure 6e is the unwrapped phase map after phase unwrapping of the difference interferograms in the embodiment of the present case on January 11, 2019 and February 4, 2019;

[0050] Figure 6f is the unwrapped phase map after phase unwrapping of the difference interferograms in the embodiment of the present case on February 4, 2019 and February 28, 2019;

[0051] Figure 6g is the unwrapped phase map after phase unwrapping of the difference interferograms in the embodiment of the present case on February 28, 2019 and March 24, 2019;

[0052] Figure 6h is the unwrapped phase map after phase unwrapping of the difference interferograms in the embodiment of the present case on March 24, 2019 and April 17, 2019;

[0053] Figure 7 is the cumulative difference phase map after accumulation in the embodiment of the present case; Figures 6a-6b

[0054] Figure 8 is the surface deformation information map after subtraction of two period DEM data obtained by the unmanned aerial vehicle in the embodiment of the present case;

[0055] Figure 9 is the schematic diagram of the "small deformation" unwrapped phase in the embodiment of the present case;

[0056] Figure 10 is the schematic diagram after adding back the large deformation phase in the embodiment of the present case;

[0057] Figure 11 is the deformation information map after encoding in the embodiment of the present case;

[0058] Figure 12a is the three-dimensional rendering diagram of the deformation result in the embodiment of the present case;

[0059] Figure 12b is Figure 12a the enlarged view of the central region;​

[0060] Figure 13 is a schematic diagram of a three-dimensional laser scanning monitoring result in the embodiment. DETAILED DESCRIPTION

[0061] In order to clearly illustrate the technical features of the patent, the patent will be described in detail below with specific embodiments and in conjunction with the accompanying drawings.

[0062] In the embodiment, the subsidence basin of the mine area obtained by the UAV of Wangjiata Coal Mine is taken as a basis, the subsidence basin is resampled according to the ground resolution of the radar satellite image, and the subsidence values of the points are converted into the absolute phases of the satellite radar; then the “integer number of turns” of each point phase is determined, and the absolute phase of each point is added to the “integer number of turns” to obtain the true phase to achieve the purpose of unwrapping; finally, the subsidence basin is obtained by using the corrected phase. The surface deformation monitoring technical process of the mine area large deformation subsidence data obtained by subtracting the two period DEM of the UAV and the InSAR data is as follows. Figures 1-2

[0063] The analysis of the process is as follows:

[0064] Firstly, the acquired n SAR images are processed, and the phase unwrapping is performed to obtain n-1 differential phases,

[0065] Secondly, the UAV obtains the DEM of the mine area, and the large deformation of the mine area is obtained by subtracting the two period DEMs,

[0066] Thirdly, the large deformation obtained by the UAV is converted into the absolute interference phase according to the radar wavelength, and is converted into the wrapped phase.

[0067] Finally, the InSAR obtains the fusion of the differential wrapped phase and the wrapped phase containing the large deformation obtained by the UAV, corrects the “integer number of turns” in the SAR phase, obtains the fusion phase containing the large deformation, and finally obtains the complete and accurate deformation of the mine area subsidence basin through the solving and processing.

[0068] Specifically:

[0069] For the data obtained by the InSAR, N-1 differential interference images can be obtained; the phase unwrapping is performed on the N-1 differential interference images, the unwrapping results are accumulated, and the differential phase of the first period and the last period image is obtained; the differential phase is de-tilted and re-wrapped to obtain a complex wrapped phase;

[0070] For the data obtained by the UAV aerial photography, the first period and the last period DSM data are differentially processed, the results are filtered to make the obtained deformation information as smooth as possible; the obtained deformation information is converted into the absolute interference phase, and is converted into the wrapped phase to obtain the wrapped phase containing the large deformation; ​

[0071] Subtract the wrapped phase containing large deformation from the complex wrapped phase, and perform phase unwrapping to obtain the unwrapped phase;

[0072] Add the unwrapped phase back to the absolute interferometric phase to obtain the fusion phase containing large deformation; finally, the fusion phase containing large deformation is converted into deformation information and geocoded.

[0073] In the above, the DEM data obtained by two unmanned aerial photogrammetry on June 9, 2018 and April 16, 2019 is used for calculation, and the results are shown in Figure 3a 、 3b Since the DEM data obtained by the unmanned aerial vehicle covers a smaller range and cannot cover the area monitored by the SAR image, the 30m SRTM data is used as the DEM data for differential interference processing.

[0074] In addition, the calculation uses 9 RadarSat-2 SAR images from June 9, 2018 to April 17, 2019, as shown in Table 1; since the image range covers an area of 125km x 125km, in order to improve the data processing efficiency and reduce the processing difficulty, the monitoring area is cropped. The cropped multi-view image is shown in Figure 4 .

[0075] Since a RadarSat-2 SAR image is very large, but the area to be processed is relatively small, therefore, before processing the image, all the original images need to be cropped, and the cropped area is the range of the mining area to be studied, therefore, the same cropping processing is performed on the 9 images, since this step does not have other important information to display, it is only used to show the cropping effect, therefore Figure 4 is an example of the range of the area presented after cropping an image.

[0076] Table 1

[0077]

[0078] Differential interference processing of data obtained by InSAR.

[0079] Since the time interval between the first image and the last image is long, the ground object changes and the surface deformation are large, direct differential interference processing will result in serious decorrelation, therefore, the 9 SAR images are arranged in chronological order, two adjacent images are processed by two-track differential interference, a total of 8 interference pairs, and the DEM selects 30m SRTM data. The interference pair information is shown in Table 2 and Fig. 5.

[0080] Table 2

[0081] Number Primary image Secondary image Time baseline / day Vertical baseline / m 1 20180609 20180727 48 -27.00 2 20180727 20180820 24 -81.00 3 20180820 20181124 96 132.96 4 20181124 20190111 48 18.16 5 20190111 20190204 24 38.43 6 20190204 20190228 24 -250.08 7 20190228 20190324 24 -90.04 8 20190324 20190417 24 236.64

[0082] The phase unwrapping is performed on 8 interferometric pairs, and the areas with coherence less than 0.2 are masked, and the result is shown in Fig. 6.

[0083] The phase values in the masked areas of the unwrapped phase are interpolated, and the 8 unwrapped phase maps are accumulated to obtain the differential phase between 20180609-20190417. The obtained cumulative differential phase is shown in Figure 7

[0084] Figure 7 The settlement area is included in the differential phase. It can be concluded from the differential phase that the maximum subsidence value is about 20 cm, which obviously does not match the actual subsidence situation of the 2S201 working face of Wangjiata Coal Mine.

[0085] Therefore, in this case, the stable area in the cumulative unwrapped phase in the previous step is selected as the control point, the phase trend surface is fitted and removed, that is, the phase is de-skewed to eliminate the slope effect caused by inaccurate baseline. The unwrapped phase is converted to complex wrapped phase again according to the phase cosine as the real part and the phase sine as the imaginary part.

[0086] Specifically:

[0087] The DSM obtained on April 16, 2019 is subtracted from the DEM data obtained on June 9, 2018 to obtain "large deformation". Through threshold segmentation, the area with obviously large deformation value (non-mining influence) and the area with obviously large deformation gradient (DSM plane error leading to misalignment of mountain slopes, noise and other influences) are masked, and TIN is constructed in the masked area. The result is low-pass filtered and down-sampled to 5m to make the obtained deformation information as smooth as possible and consistent with the SAR image resolution. The obtained deformation information is shown in Figure 8

[0088] The obtained deformation information is "large deformation", which is converted to absolute interferometric phase according to the radar wavelength, and then converted to wrapped phase.

[0089] The complex wrapped phase includes "large deformation" which cannot be unwrapped due to large deformation gradient; the wrapped phase obtained in the previous step is also "large deformation". The wrapped phase obtained in the previous step is deducted from the complex wrapped phase to only retain "small deformation", so that the phase unwrapping can be performed. The unwrapped phase of the residual "small deformation" after deducting "large deformation" is shown in Figure 9

[0090] Since "large deformation" cannot be phase unwrapped, the absolute interferometric phase after conversion of "large deformation" is retained. The absolute interferometric phase is added back to the unwrapped phase with only residual small deformation, and the obtained result is shown in Figure 10

[0091] ​​​​The phase obtained in the last step is converted into deformation information and geocoded. The coded deformation information is as shown in FIG. 8: there are two lowest points in the center, which are the subsidence areas with a subsidence of-2.8 m. Figure 11

[0092] The three-dimensional rendering diagram of the obtained deformation result is as shown in FIG. 12.

[0093] As can be seen from the above figure, on April 15, 2019, the maximum subsidence value of the surface of Wangjiata Coal Mine was 2.8 m, which was basically consistent with the three-dimensional laser scanning monitoring result 2.725 m. The subsidence basin calculated after fusion is superimposed on the Wangjiata Coal Mine mining engineering plan, as shown in FIG. 13. Figure 13

[0094] It is shown from the embodiment that the method is applied in the Wangjiata Coal Mine in Inner Mongolia, and good results are obtained. The engineering case shows that the UAV / InSAR fusion technology has a significant advantage in the surface subsidence deformation monitoring of the western mining area, and the research results provide technical support for the surface damage monitoring of the western coal mining in China.

[0095] The present application has many specific implementation ways, and the above description is only the preferred embodiment of the present application. It should be pointed out that for ordinary skilled persons in the art, some improvements can be made without departing from the principles of the present application, and these improvements should also be regarded as the protection scope of the present application.​​

Claims

1. A method for obtaining high-precision surface subsidence basins by fusing UAV DEM and InSAR phase fusion, characterized in that, Follow these steps: S1, Data Acquisition; N SAR images within the target time period were acquired using InSAR, and two DEM data were acquired sequentially within the target time period using UAV aerial photography. S2, Data Conversion; The N SAR images obtained in step S1 are subjected to interferometric and unwrapping processing to obtain N-1 unwrapped phase images. Through superposition and rewrapping processing, complex wrapped phase images are obtained. Subtract the two DEM data obtained in step S1 to obtain the large deformation of the mining area subsidence; then convert the large deformation of the mining area subsidence into an absolute interference phase according to the radar wavelength, and then convert it into a tangled phase to obtain a tangled phase containing the large deformation. S3, Phase Fusion; The large-deformation winding phase obtained in step S2 is fused with the complex winding phase, and the integer number of cycles in the SAR phase is corrected to obtain the fused phase with large deformation. The specific process of phase fusion in step S3 is as follows: the entangled phase containing large deformation obtained in step S2 is subtracted from the complex entangled phase, and only small deformation is retained to obtain the unentangled phase of the residual small deformation; the absolute interference phase after DEM data conversion in step S2 is added back to the unentangled phase of the residual small deformation. S4, Data Solving; Through calculation and processing, the deformation of the mining area subsidence basin is finally obtained completely and accurately.

2. The method for obtaining high-precision surface subsidence basins by phase fusion of UAV DEM and InSAR according to claim 1, characterized in that, In step S2, when processing N SAR images, the SAR images are arranged in chronological order. First, two-track differential interferometry is performed on two adjacent images, and then phase unwrapping is performed. The region with coherence less than 0.2 is masked to obtain N-1 unwrapped phase images. After that, the phase value is interpolated in the masked region of the unwrapped phase, and the N-1 unwrapped phase images are accumulated. The stable region in the accumulated unwrapped phase is selected as the control point, the phase trend surface is fitted and removed, that is, the phase deskewing process is performed to eliminate the slope effect caused by the inaccurate baseline; the deskewing interference phase is converted back into a complex winding phase according to the phase cosine as the real part and the phase sine as the imaginary part.

Images