An air-ground combined InSAR geological disaster monitoring method and device

Abstract

The present application belongs to the field of surveying and mapping engineering, and particularly relates to a kind of air-ground combined InSAR geological disaster monitoring method and device, method includes: respectively obtaining satellite digital orthophoto map of disaster area and first digital orthophoto map of disaster point, and satellite digital orthophoto map and first digital orthophoto map are matched with elevation model;Satellite digital orthophoto map and first digital orthophoto map are fused to obtain first fused digital orthophoto map;First digital elevation model and second digital elevation model are fused to obtain first fused digital elevation model;The mapping relationship between first fused digital orthophoto map and first fused digital elevation model is established to obtain air-ground combined three-dimensional terrain model;Ground-based SAR is fused with air-ground combined three-dimensional terrain model and is handled to obtain InSAR stereoscopic detection model.The present application solves the problems of low acquisition efficiency, limited monitoring range, insufficient data precision and poor real-time performance of traditional monitoring methods.

Description

TECHNICAL FIELD

[0001] The application belongs to the field of surveying and mapping engineering, and particularly relates to an air-ground combined InSAR geological disaster monitoring method and device. BACKGROUND

[0002] In recent years, various natural disasters occur frequently in China, showing the characteristics of wide disaster influence range, obvious disaster area characteristics, significant monsoon climate influence, and heavy geological disaster loss. In addition, engineering construction has changed the original geological environment, resulting in situations such as collapse, landslide, debris flow, ground subsidence, ground fissure and ground subsidence, which directly cause social and economic losses.

[0003] For geological disaster monitoring, domestic experts have proposed many monitoring methods, such as space-based monitoring technology. Space-based monitoring mainly relies on optical satellite images, synthetic aperture radar (SAR) images, satellite positioning, satellite-borne precipitation radar, etc. These methods are suitable for large-scale geological disaster screening and monitoring, have a long data acquisition period, are not high in accuracy, and are high in cost. Air-based monitoring mainly uses unmanned aerial vehicle low-altitude remote sensing technology to monitor various geological disaster macroscopic deformation signs, development characteristics and disaster rules, and is suitable for large-scale geological monitoring. The data period is higher than that of space-based monitoring, and the accuracy is affected by the attitude of the aircraft. In bad weather such as rain, snow and fog, unmanned aerial vehicles are difficult to obtain surface subsidence observation data, and there is a risk of crashing. Land-based monitoring mainly uses global navigation satellite system (GNSS) receivers and total station measurement methods to monitor various geological disaster changes. Although the data accuracy is high, it is a non-continuous monitoring method, and the observation data cannot directly represent the surface deformation value of the area far away from the observation pier. For landslides, debris flow and ground subsidence, GNSS and total station monitoring methods are not suitable, and the operation risk is high.

[0004] In view of the limitations of the prior art, the present application aims to provide an air-ground combined InSAR (Interferometric Synthetic Aperture Radar) geological disaster monitoring method, which can quickly establish a three-dimensional terrain model of the disaster area, realize large-scale, high-precision and real-time geological disaster monitoring, and improve the accuracy and response speed of disaster warning. SUMMARY

[0005] The present application provides an air-ground combined InSAR geological disaster monitoring method and device to solve the problems of low data acquisition efficiency, limited monitoring range, insufficient data accuracy, poor real-time performance and complex data processing of traditional methods.

[0006] To solve the above problems, the present application adopts the following scheme:

[0007] The application discloses an air-ground joint InSAR geological disaster monitoring method.

[0008] In a first aspect, the application provides an air-ground joint InSAR geological disaster monitoring method, comprising the following steps:

[0009] Satellite digital orthophoto images of a disaster area are acquired, and a first digital elevation model is matched with the satellite digital orthophoto images;

[0010] First digital orthophoto images of disaster points in the disaster area are acquired, and a second digital elevation model is matched with the first digital orthophoto images; the first digital orthophoto images are acquired by a UAV;

[0011] The satellite digital orthophoto images and the first digital orthophoto images are fused to obtain first fused digital orthophoto images; the first digital elevation model and the second digital elevation model are fused to obtain a first fused digital elevation model; a mapping relationship between the first fused digital orthophoto images and the first fused digital elevation model is established to obtain an air-ground joint three-dimensional terrain model;

[0012] Ground-based SAR is fused with the air-ground joint three-dimensional terrain model to obtain an InSAR stereoscopic detection model; and real-time monitoring of geological disasters in a disaster area is performed based on the InSAR stereoscopic detection model.

[0013] Preferably, the fusion of the ground-based SAR with the air-ground joint three-dimensional terrain model comprises the following steps:

[0014] InSAR orbit coordinates are acquired;

[0015] A radar deformation azimuth angle is determined;

[0016] Based on the InSAR orbit coordinates and the radar deformation azimuth angle, InSAR radar data is registered to the air-ground joint three-dimensional terrain model to obtain an initial InSAR stereoscopic detection model;

[0017] InSAR data observation points are established on the initial InSAR stereoscopic detection model to obtain an InSAR stereoscopic detection model.

[0018] Preferably, after the InSAR stereoscopic detection model is obtained, the InSAR observation data obtained from the InSAR data observation points are further optimized, and the optimization comprises:

[0019] In combination with the average correlation coefficient and amplitude deviation method, a correlation coefficient of adjacent images in the InSAR observation data is calculated, the correlation coefficient under a preset proportion is set as a threshold parameter, the InSAR data observation points are corrected, and the final InSAR data observation points are determined.

[0020] Preferably, the step of obtaining InSAR track coordinates comprises:

[0021] The coordinates of the InSAR device of the disaster point in the disaster area are obtained by aerial surveying and mapping by the unmanned aerial vehicle.

[0022] Based on the coordinates of the InSAR device, the coordinate values between the InSAR track end points are measured by using a single point coordinate measurement method.

[0023] Preferably, determining the radar deformation azimuth angle specifically comprises:

[0024] Obtaining InSAR original data from InSAR data observation points;

[0025] Obtaining cumulative displacement of each InSAR data observation point in the InSAR original data by differential interference method; obtaining deformation data of each InSAR data observation point;

[0026] Differential processing of deformation data obtained at different time points to obtain relative deformation data within a corresponding time period;

[0027] According to the deviation between the deformation data and the InSAR original data, the radar deformation azimuth angle is calculated.

[0028] Preferably, according to the deviation between the deformation data and the InSAR original data, the radar deformation azimuth angle is calculated specifically comprising:

[0029] According to the deformation data, the deformation coordinates of the InSAR track end points are obtained;

[0030] According to the InSAR original data, the original coordinates of the InSAR track end points are obtained;

[0031] The coordinates of the InSAR device are connected with the deformation coordinates of the InSAR track end points and the original coordinates of the InSAR track end points respectively, and the radar deformation azimuth angle is obtained.

[0032] Preferably, registering the InSAR observation data to the air-ground combined three-dimensional terrain model specifically comprises: registering according to the terrain point cloud data, track coordinates and InSAR observation data.

[0033] Preferably, it further comprises: after establishing the initial InSAR stereo detection model, verifying the spatial relationship to check whether the position of the coordinates of the InSAR radar observation data is consistent with the position on the air-ground combined three-dimensional terrain model, and if not, re-registering.

[0034] Preferably, it further comprises: by using a radar solving tool, solving the radar interference file in the format of DiffImage into a deformation map spot and registering it to the InSAR stereo detection model.

[0035] Preferably, the first digital orthographic image of the disaster point in the disaster area is obtained by using a UAV, the relative height of the UAV to the ground is controlled to be 30-40m, the image overlap rate of the heading and the lateral direction is 75%-85%, and the UAV flight area is greater than the disaster point by 50m.

[0036] Preferably, the image resolution of the first digital orthographic image is less than 5cm / pix, and pix represents a pixel.

[0037] Preferably, after obtaining the first digital orthographic image of the disaster point, the method further comprises calculating the obtained laser radar data into three-dimensional point cloud data.

[0038] Preferably, the construction method of the second digital elevation model comprises: obtaining laser radar data by using a UAV carrying a laser radar, calculating three-dimensional point cloud data according to the obtained laser radar data, and separating the point cloud of the disaster point by point cloud cropping.

[0039] Preferably, the separating of the point cloud of the disaster point by point cloud cropping specifically comprises:

[0040] By establishing a deep learning classification sample of point cloud data, the point cloud data is classified and displayed as ground point cloud data by using a Patchwork++ algorithm, the non-displayed part is non-ground point cloud data, the Patchwork++ algorithm is optimized by adjusting the deep learning iteration parameters and increasing the deep learning training sample, and accurate ground point cloud data is obtained.

[0041] Based on the same idea, an air-ground combined InSAR geological disaster monitoring device is also proposed, comprising at least one processor and a memory in communication connection with the at least one processor; the memory stores instructions executable by the at least one processor, and the instructions are executed by the at least one processor to enable the at least one processor to implement any one of the above-mentioned air-ground combined InSAR geological disaster monitoring methods.

[0042] Compared with the prior art, the beneficial effects of the present application are:

[0043] The application provides an air-ground combined InSAR geological disaster monitoring method, specifically comprising air-space three-dimensional terrain model fusion, ground InSAR and three-dimensional terrain fusion, PSInSAR (Persistent Scatterer Interferometric Synthetic Aperture Radar) time series data solving and multi-time series ground InSAR stereoscopic monitoring, ensuring the continuity of ground monitoring results, realizing real-time monitoring and early warning of geological disasters, improving the three-dimensional visual monitoring effect of geological disasters, and solving the problems of long data observation period, low acquisition efficiency, limited monitoring range, insufficient data precision, poor real-time performance and complex data processing of the traditional monitoring method. BRIEF DESCRIPTION OF DRAWINGS

[0044] Figure 1 It is a flow chart of an air-ground combined InSAR geological disaster monitoring method in embodiment 1 of the application.

[0045] Figure 2 It is an effect diagram of unmanned aerial vehicle, satellite, InSAR collected data fusion in embodiment 1 of the application.

[0046] Figure 3 It is an InSAR registration diagram in embodiment 1 of the application.

[0047] Figure 4 It is a ground surface overall subsidence diagram and subsidence amplitude change diagram in embodiment 2 of the application.

[0048] Figure 5 It is a diagram of change of ground surface overall subsidence with time in embodiment 2 of the application.

[0049] REFERENCE NUMERALS:

[0050] 1-satellite image, 2-unmanned aerial vehicle DOM image, 3-ground collapse area, 4-InSAR registration area, 5-InSAR position, 6-unmanned aerial vehicle DEM model. DETAILED DESCRIPTION

[0051] The application will be further described in detail in combination with test examples and specific embodiments. However, this should not be understood as limiting the scope of the above-mentioned subject matter of the application to the following examples, and any technology realized based on the content of the application belongs to the protection scope of the application.

[0052] In the description of specific embodiments of the present application, the orientation or position relationship indicated by the terms such as "upper", "lower", "left", "right", "center", "inner", "outer", "side" and the like appearing in the description of specific embodiments of the present application are expressed based on the orientation or position relationship shown in the drawings, or are the orientation or position relationship when the product / device / apparatus is usually used. These orientation or position relationship terms are only for the convenience of describing the present application or simplifying the description in specific embodiments, for the convenience of technicians to quickly understand the scheme, and therefore cannot be understood as indicating or implying that a specific device / component / element must have a specific orientation, or be constructed and operated in a specific position relationship, and therefore cannot be understood as limiting the present application.

[0053] In the description of specific embodiments of the present application, the technical terms "first", "second" and the like only distinguish one entity or operation from another entity or operation, and cannot be understood as indicating or implying relative importance or implicitly indicating the number, specific order or primary and secondary relationship of the indicated technical features. In the description of specific embodiments of the present application, the meaning of "multiple" is two and more than two, unless otherwise explicitly and specifically limited.

[0054] Reference herein to "an embodiment" means that a particular feature, structure, or characteristic described in connection with an embodiment can be included in at least one embodiment of the application. The appearance of the phrase in various places in the specification does not necessarily all refer to the same embodiment, nor is it necessarily independent or alternative embodiments to each other. It is explicitly and implicitly understood by those skilled in the art that the embodiments described herein can be combined with other embodiments.

[0055] Embodiment 1

[0056] An air-ground joint InSAR geological disaster monitoring method, a flowchart as shown in Figure 1 The method comprises the following steps:

[0057] Obtaining a satellite digital orthophoto map of a disaster area, and matching a first digital elevation model to the satellite digital orthophoto map;

[0058] Obtaining a first digital orthophoto map of a disaster point in the disaster area, and matching a second digital elevation model to the first digital orthophoto map; the first digital orthophoto map is obtained by a UAV;

[0059] Fusing the satellite digital orthophoto map and the first digital orthophoto map to obtain a first fused digital orthophoto map; fusing the first digital elevation model and the second digital elevation model to obtain a first fused digital elevation model; establishing a mapping relationship between the first fused digital orthophoto map and the first fused digital elevation model to obtain an air-ground joint three-dimensional terrain model;

[0060] The foundation SAR is fused with the three-dimensional terrain model of air-ground combination to obtain an InSAR stereo detection model; and real-time monitoring of geological disasters in the disaster area is performed based on the InSAR stereo detection model.

[0061] Taking an actual engineering project as an example, referring to Figures 2 to 4 The embodiment improves the geological disaster monitoring method, and proposes an air-ground combined InSAR geological disaster monitoring method, which specifically comprises the following steps:

[0062] I. Air-space three-dimensional terrain model fusion processing

[0063] (1) Satellite terrain rapid acquisition

[0064] In order to quickly establish a three-dimensional terrain model of the disaster area, satellite digital orthophoto map (DOM) of the disaster area is downloaded through BIG MAP software, and the pixel resolution is 2 m; DEM data is applied for from the surveying and mapping geographic information bureau, and the highest terrain precision is 12.5 m, and the data effect is as shown in Figure 2 The peripheral area in the figure is satellite image 1, the area close to the disaster area relative to the peripheral area is unmanned aerial vehicle DOM image 2, and the area close to the ground collapse area 3 relative to the unmanned aerial vehicle DOM image 2 is InSAR registration area 4, forming the characteristics that the image from the periphery to the ground collapse area 3 is clearer.

[0065] (2) Establish high-resolution DOM image

[0066] Here, the high-resolution digital orthophoto map (DOM) is obtained by reconstructing the unmanned aerial vehicle optical aerial photograph. Since the terrain of the collapse area is large, the unmanned aerial vehicle adopts full-automatic flight to collect images, which has the risk of crashing, therefore, the research group adopts a manual flight control mode. In order to improve the resolution of the DOM image, the relative height of the aircraft to the ground should be controlled at 30-40 m, and the image overlap rate in the heading and lateral directions is 75%-85%. At the same time, in order to quickly realize the monitoring of the ground collapse area 3, the unmanned aerial vehicle flight area is greater than the collapse area by about 50 m, the purpose is to reduce the unmanned aerial vehicle aerial photograph, quickly establish the image data of the rescue site, and provide services for the rescue command center. For the processing of the unmanned aerial vehicle DOM image 2, the DJI ZhiTu software is adopted, the multi-view unmanned aerial vehicle aerial photograph is imported, the coordinate reference is CGCS2000 (China Geodetic Coordinate System 2000), the local central meridian longitude is selected, the high-resolution digital orthophoto map (DOM) is established, and the image resolution is recommended to be less than 5 cm / pix (pix represents pixel), so as to facilitate the viewing of the ground cracking texture area.

[0067] (3) Establish high-precision DEM model

[0068] Since the high-precision digital elevation model (DEM) cannot be directly obtained, it needs to be obtained by processing point cloud data. The present embodiment adopts DJI unmanned aerial vehicle, which is equipped with lidar. In order to facilitate the viewing of the spatial positional relationship between the tunnel and the collapse body, the lidar point cloud data from the collapsed tunnel entrance to the ground surface in the collapse body area needs to be established. In order to improve the accuracy of the acquisition of the lidar point cloud data, the echo and true color mode are adopted, and a total of 5 flights are completed to collect the lidar point cloud data. After the data collection is completed, the original lidar point cloud data is calculated into three-dimensional point cloud by DJI ZhiTu, and the point cloud quantity established at the rescue site is 400497588.

[0069] In order to improve the efficiency of ground monitoring, the point cloud around the collapse body is segmented and retained by the point cloud cutting tool. By marking 30-50 feature points on the point cloud data, a point cloud deep learning classification sample is established by the Patchwork++ algorithm. The point cloud displayed by the algorithm classification is the ground point, and the non-displayed part is the non-ground point. PointNet++ is suitable for small-scale point cloud classification, directly processes unordered point sets, and extracts global and local features. In addition, DGCNN algorithm, Voxel-Based method and Transformer model can also be used. Dynamic Graph CNN (DGCNN, Dynamic Graph Convolutional Neural Network, a deep learning model based on graph convolutional neural network) captures the point cloud topology structure through dynamic graph convolution, which is suitable for complex terrain. Voxel-Based method (such as VoxelNet) is used to voxelize the point cloud and use 3D CNN for processing, which is suitable for high-density data. Transformer model (such as Point Transformer) can use self-attention mechanism to process long-range dependencies, which is suitable for large-scale scenarios. The above deep learning algorithms can output the segmented results in the classification task, such as distinguishing disaster areas (such as landslide boundary, subsidence area) and stable areas. Since single classification cannot accurately classify all ground points, the classification results can be optimized by adjusting the deep learning iteration parameters and increasing the number of deep learning training samples to improve the classification accuracy of the ground points. During iterative training, weighted cross-entropy loss (to solve the class imbalance, such as the scarcity of disaster points) is adopted, further segmented learning rate (initial 0.001, decay every 10 rounds) is adopted, and AdamW optimizer is also used to avoid overfitting. Since there are cliffs around the collapse body, the classification accuracy of the point cloud in the cliff area cannot be directly solved by optimizing the iteration parameters and training samples, and the ground points in the cliff area need to be classified by point cloud segmentation and resampling method, and then combined with the ground point cloud data.

[0070] Further, based on the curvature feature (such as PCA local plane fitting, Principal Component Analysis-based Local Plane Fitting, which is a principal component analysis algorithm) to segment the high curvature area, the cliff can be identified.

[0071] Delaunay triangulation (also called TIN triangulation) method, the ground point cloud data is converted into a triangular network, in order to reduce the number of triangular networks, ensure smooth computer operation, the edge length less than 0.1m, edge length greater than 30m above the triangular network is filtered out, at the same time, it can also improve the speed of establishing rescue three-dimensional scene. Then through the TIN triangulation DEM production method, set DEM interval parameters X=0.2m, Y=0.2m, to establish high-precision digital elevation model (DEM).

[0072] (4) Space, three-dimensional terrain model

[0073] First, satellite image 1 and unmanned aerial vehicle DOM; image fusion processing, using GIS software, through geometric correction to unify the coordinate reference of the two images, and then using IHS parameter fusion method to realize the fusion of satellite image 1 and unmanned aerial vehicle DOM, wherein, intensity (intensity) is set to 4-8, hue (hue) is set to 4-7, and saturation (saturation) is set to 8-15.

[0074] Second, to realize the fusion of DEM, the sparse linear expression method is used to realize the seamless fusion between satellite DEM and unmanned aerial vehicle DEM data, wherein the sparse line value is 25-40.

[0075] Finally, the mapping superposition relationship between the above fused DEM and DOM is established, that is, the three-dimensional terrain model of space-ground joint is obtained.

[0076] After fusing the digital elevation model (DEM) and the digital orthographic image (DOM), the two are combined to build a three-dimensional visual expression of the geographic scene through spatial mapping superposition relationship, and the specific relationship is as follows:

[0077] First step, determine the core relationship of mapping superposition

[0078] First, align the spatial coordinates, set the plane coordinates as (X, Y): DEM (terrain elevation data) and DOM (ground image data) are strictly aligned through the same plane coordinate system (such as UTM or WGS84), to ensure that each pixel (DOM) and the corresponding elevation point (DEM) are in the same position. Secondly, map the RGB information of DOM to the Z value (elevation value) of DEM to form a 2.5D surface model (i.e. each XY position corresponds to only one Z value). The data coupling mode is shown in Table 1:

[0079] Table 1 Data Coupling Methods

[0080]

[0081] Second Step: Implementation of Superimposition Relationship

[0082] The implementation of superimposition relationship specifically includes geometric correction and registration; texture mapping. The specific implementation of geometric correction and registration includes:

[0083] DOM processing: eliminate image distortion (such as lens distortion of UAV images, terrain displacement) through orthographic correction, or use resolution matching (such as 0.1m resolution of DOM needs to be aligned with 0.1m grid of DEM).

[0084] DEM processing: eliminate data voids (such as interpolate to complete missing areas of lidar point cloud data), or smoothing processing (such as Gaussian filtering to retain terrain trends and remove noise).

[0085] Texture mapping specifically includes: binding the RGB values of DOM as texture to the triangular mesh (TIN) or regular grid of DEM, and realizing seamless fitting through UV coordinate mapping.

[0086] Python code as follows:

[0087] # DEM + DOM Superimposition Based on GIS Library

[0088] import rasterio

[0089] from pyvista import dem_to_3d

[0090] dem = rasterio.open("dem.tif").read(1) # Read DEM

[0091] dom = rasterio.open("dom.tif").read([1,2,3]) # Read DOM (RGB three bands)

[0092] mesh = dem_to_3d(dem) # Convert DEM to 3D mesh

[0093] mesh.textures["dom"] = dom.transpose(1,2,0) # Bind DOM texture

[0094] Further, the resolution needs to be set, DOM resolution ≥ DEM resolution (to avoid texture blur), for example: if the DEM is 1 m grid, the corresponding DOM needs to be set to a resolution greater than or equal to 0.5 m resolution.

[0095] Further, the time should be consistent, and the time difference between the two acquisitions should be as small as possible (such as the post-landslide DOM needs to be synchronized with the post-disaster DEM). The coordinate system should be unified to ensure that the DEM and the DOM use the same projection (such as CGCS2000) and vertical datum (such as EGM96 elevation).

[0096] II. Ground-based SAR and three-dimensional terrain fusion processing

[0097] (1) Obtain InSAR orbit coordinates

[0098] Since ground-based InSAR does not have satellite positioning function, the observed data does not have coordinate values, therefore it needs to be converted to the coordinate system of unmanned aerial lidar, i.e. CGCS2000 coordinate system, for this purpose, the InSAR equipment needs to be installed before the aerial flight terrain data is collected, and then the flight survey operation is carried out. On the three-dimensional real scene model of the collapsed area, find the position of the InSAR equipment, and use the single point coordinate measurement method to measure the geodetic coordinate values (B, L, H) and the projection coordinate values (X, Y, H) of the left and right two end points of the InSAR orbit respectively, and save them as txt text.

[0099] (2) Determine the radar deformation azimuth angle

[0100] The InSAR original data can obtain the cumulative displacement of each target point at that time through differential interference technology, and the deformation data at T1 and T2 is processed by difference to obtain the relative deformation data in the time period. The radar deformation azimuth angle of this embodiment is 120°, and by adding the coordinate values of the two ends of the orbit on the three-dimensional terrain model, the radar deformation azimuth angle can be calculated.

[0101] Preferably, by adding the coordinate values of the two ends of the orbit on the three-dimensional terrain model, the radar geometric observation direction (i.e. azimuth direction) is essentially related to the geographic coordinate system. The specific steps of calculating the radar deformation azimuth angle are as follows:

[0102] First, prepare the radar orbit endpoint coordinates. Obtain the two endpoint coordinates of the radar satellite orbit on the ground projection, which usually comes from satellite ephemeris data or SAR image metadata. Use DEM data to construct the terrain surface to ensure that the coordinate system is consistent with the radar coordinate.

[0103] Second, calculate the radar azimuth vector.

[0104] Horizontal projection vector: ignore elevation (only use XY coordinates), calculate the horizontal projection vector of radar flight direction :

[0105] ;

[0106] where, and Two end point coordinates of radar satellite orbit projection on the ground.

[0107] Calculate the horizontal azimuth angle of radar flight direction ( ): Calculate the horizontal projection vector The angle between the north direction (clockwise), calculated using the arctangent function, formula:

[0108] ;

[0109] Convert the result to the angle system (°), range 0°~360°, example:

[0110] If the two end point coordinates of radar satellite orbit projection on the ground are M1(100,200), M2(150,180), then:

[0111] ;

[0112] ;

[0113] Third step, calculate the relationship between deformation azimuth angle and radar azimuth.

[0114] Deformation azimuth angle definition: deformation azimuth angle (such as 120°) refers to the angle between the ground deformation direction and the north direction (clockwise), which needs to be compared with the radar azimuth ( ) to analyze the deformation sensitivity.

[0115] Radar geometric constraints: radar is not sensitive to deformation along the azimuth direction ( ) (blind area), and is most sensitive to deformation perpendicular to the azimuth direction ( ).

[0116] Verification example:

[0117] If =111.8°, then:

[0118] Radar sensitive direction: 111.8°+90°=201.8° (i.e. southwest-northeast direction). The angle between the deformation azimuth angle 120° reported by the user and the radar sensitive direction (201.8°) is: |120°-201.8°|=81.8°, the angle is close to 90°, which indicates that the deformation direction is close to the radar sensitive direction, and the data is reliable.

[0119] The code for implementing radar deformation azimuth calculation with python is as follows:

[0120] import numpy as np

[0121] # Input track endpoint coordinates (x, y, z)

[0122] P1 = np.array([100, 200, 50])# Start point

[0123] P2 = np.array([150, 180, 55])# End point

[0124] # Calculate horizontal azimuth (degrees)

[0125] dx, dy = P2[0] - P1[0], P2[1] - P1[1]

[0126] azimuth_deg = np.degrees(np.arctan2(dx, dy)) % 360# Convert to 0~360°

[0127] print(f"Radar azimuth: {azimuth_deg:.1f}°")

[0128] (3) InSAR track coordinate registration

[0129] Use the radar coordinate registration tool to register the InSAR radar data to the three-dimensional terrain model to improve the three-dimensional stereo monitoring effect. In the registration tool, 5 files need to be prepared, of which 3 are terrain point cloud files, 1 is track coordinate file, and 1 is InSAR radar observation file;

[0130] In the InSAR registration interface, the first item is the point cloud file of the independent coordinate system after cutting, if there is no independent coordinate system, it is the point cloud file of the projection coordinate system; The second item is the radar track coordinate, which must be in the same coordinate system as the first item; The third item is the diffimage graph of the radar; The fourth item is the point cloud file of the CGCS2000 geodetic coordinate system after cutting; The fifth item is the point cloud file of the CGCS2000 projection coordinate system after cutting, if there is no independent coordinate system, it is the same as the first item. The InSAR registration diagram is shown in Figure 3 The red line area in the figure is the surface collapse area 3, the yellow area larger than the surface collapse area 3 is the InSAR registration area 4, and the red dot in the figure is the InSAR position 5. The area larger than the InSAR registration area 4 is the UAV DEM model 6.

[0131] (4) InSAR and three-dimensional terrain fusion

[0132] Radar position registration, after radar coordinate registration, select to verify the spatial relationship, check if the position of the radar coordinate is consistent with the position on the UAV tilt model, if not, it needs to be re-registered until the radar deformation map is fused into the three-dimensional terrain model. Secondly, fuse the radar deformation map, through the radar solving tool, calculate the radar interference file in the format of DiffImage into the deformation map spot, and register it to the DEM model and DOM image.

[0133] (5) Establish InSAR data observation

[0134] Using multi-source data fusion, creating monitoring work tasks on the monitoring and early warning platform, adding the ID of the radar to the equipment management, inputting the username and password on the login page of the slope radar monitoring and early warning system, clicking the login button to log in; Click the drop-down triangle on the right side of the application field to expand all groups under the application field, click the drop-down triangle in front of the group to expand all companies under the group, click the drop-down triangle in front of the company to expand all projects under the company, click the project name, and the webpage jumps to the three-dimensional scene to display the topography of the first project under the project. The three-dimensional earth displays the icons of all radars and sensors under the project.

[0135] III. PSInSAR (PS) time series data solving

[0136] In order to eliminate the influence of InSAR registration area 4 data on correlation coefficient, amplitude deviation, phase deviation and other factors, combined with the average correlation coefficient and amplitude deviation method, the correlation coefficient of adjacent images is calculated on the average correlation coefficient algorithm, at the same time, the correlation coefficient under the set proportion is taken as the threshold parameter, the final target point selection is carried out, and the real-time performance of ground-based SAR data processing is ensured. The main and auxiliary images are formed into an interference pair by transforming the adjacent two images, the number and quality of PS (Persistent Scatterer) point extraction are optimized, the adaptive threshold of average correlation coefficient is used, k ground-based SAR complex images are selected to form k interference pairs, and the correlation coefficient of the interference pair is calculated. Finally, the amplitude mean value and amplitude standard deviation of each pixel are calculated in turn by using the amplitude value of K ground-based InSAR images. In this example, a total of 17745 effective PS points are obtained, and the registration rate is 93.89%. In the rescue site, considering the time is relatively tense, when the PS effective point reaches more than 60%, it can be used for on-site observation.

[0137] PS, Persistent Scatterer technology can keep stable scattering characteristics of ground targets, such as natural stable rock surfaces, in a long time sequence of radar images. The stable backscattering and phase characteristics of rocks over time enable them to be used for high-precision ground deformation monitoring, achieving millimeter-level precision in ground deformation measurement.

[0138] Embodiment 2

[0139] The verification process of the air-ground combined InSAR geological disaster monitoring method is as follows:

[0140] (1) Upload observation period data

[0141] In this rescue example, the observation station of the ground-based SAR was established at 10:00 on May 5, and the observation data was started at 10:34. The first time series of dry radar images was obtained at 10:42. According to the requirements of the rescue command center, the observation station ended at 23:31 on May 7, and a total of 446 time series of radar dry images were obtained. It takes 8 minutes to observe one time series of radar data.

[0142] (2) Upload multi-time series InSAR data

[0143] Based on the geological disaster monitoring and early warning platform, by adding radar ID name, connecting computer and radar through wifi, realizing multi-time series radar data collection and uploading immediately, finding the long village geological monitoring project under the left directory tree, entering the disaster location, using mouse to click the radar icon, setting the interval time of radar data uploading to 8 minutes, consistent with the observation time of the radar.

[0144] (3) InSAR stereo monitoring and early warning

[0145] From the radar real-time data, 446 time series of radar data are superimposed simultaneously, and the overall subsidence of the entire collapse area is calculated by the radar multi-track interference phase difference method. In order to intuitively view the subsidence in the entire collapse area, the self-defined color spectrum is used, and the color is usually set to red when the subsidence value is larger. Based on the real-time monitoring platform of the slope radar, according to the first received data of the ground-based SAR, the monitoring and early warning area is selected by polygon drawing, and the subsidence of the ground collapse area 3 is monitored in real time. The early warning threshold is set to 5 cm / h, that is, when the subsidence value exceeds 5 cm per hour, the system will include the early warning, and the early warning report will be pushed to the rescue command center at the same time. Through the InSAR stereo detection model in the present application, real-time monitoring of geological disasters in the disaster area can be realized, and a plurality of subsidence values can be obtained, so as to perform early warning. The ground overall subsidence graph and the subsidence amplitude change graph obtained after analysis are shown in Figure 4 and the graph of the change of the overall ground subsidence with time is shown in Figure 5are shown.

[0146] Figure 4 In the figure, the abscissa represents the subsidence type, the ordinate represents the subsidence value (unit: cm), the black square line represents the maximum subsidence outside the surface collapse, the red square line represents the average subsidence of the surface, the blue square line represents the maximum subsidence inside the collapse, and the green square line represents the cumulative deformation value of the average subsidence.

[0147] Figure 5 In the figure, the abscissa represents the time, and the ordinate represents the subsidence value (unit: cm). Figure 5 In the figure, the cumulative deformation values of the maximum subsidence outside the surface collapse, the average subsidence of the surface, the maximum subsidence inside the collapse, and the average subsidence of the collapse pit are shown with the observation time.

[0148] Embodiment 3

[0149] The air-ground combined InSAR geological disaster monitoring device comprises at least one processor and a memory in communication connection with the at least one processor; the memory stores instructions executable by the at least one processor, and the instructions are executed by the at least one processor to implement the air-ground combined InSAR geological disaster monitoring method in Embodiment 1 or 2.

[0150] The above-described embodiments are only used to illustrate the technical solutions of the present application, rather than limit them; although the present application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that they can still modify the technical solutions recorded in the foregoing embodiments, or make equivalent replacements for part of the technical features; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the spirit and scope of the technical solutions of the embodiments of the present application, and should be included in the protection scope of the present application.

Claims

1. A combined air-ground InSAR method for monitoring geological hazards, characterized in that, Includes the following steps: Acquire satellite digital orthophotos of the disaster-stricken area and match the satellite digital orthophotos with a first digital elevation model; Acquire the first digital orthophoto of the disaster-stricken area and match it with a second digital elevation model; the first digital orthophoto is acquired by a drone. The satellite digital orthophoto and the first digital orthophoto are fused to obtain a first fused digital orthophoto; the first digital elevation model and the second digital elevation model are fused to obtain a first fused digital elevation model; Establish the mapping relationship between the first fused digital orthophoto and the first fused digital elevation model to obtain a three-dimensional terrain model combining air and ground features; The ground-based SAR and the air-ground combined three-dimensional terrain model are fused to obtain the InSAR stereo detection model. Real-time monitoring of geological disasters in disaster areas based on the InSAR three-dimensional detection model; The process of fusing the ground-based SAR with the air-ground combined three-dimensional terrain model specifically includes the following steps: Obtain InSAR orbit coordinates; Determine the radar deformation azimuth angle; Based on the InSAR orbital coordinates and the radar deformation azimuth angle, the InSAR radar data is registered to the air-ground joint three-dimensional terrain model to obtain the initial InSAR stereo detection model. InSAR data observation points are established on the initial InSAR stereo detection model; thus, the InSAR stereo detection model is obtained.

2. The air-ground combined InSAR geological disaster monitoring method as described in claim 1, characterized in that, After obtaining the InSAR stereo detection model, the InSAR observation data acquired from the InSAR data observation points are further optimized. The optimization process includes: By combining the average correlation coefficient and the amplitude deviation method, the correlation coefficient of adjacent images in the InSAR observation data is obtained. The correlation coefficient at a preset ratio is set as a threshold parameter to correct the InSAR data observation points and determine the final InSAR data observation points.

3. The air-ground combined InSAR geological disaster monitoring method as described in claim 1, characterized in that, The steps to obtain InSAR orbit coordinates include: The coordinates of InSAR devices at disaster-stricken locations in the disaster area were obtained through aerial surveying using unmanned aerial vehicles (UAVs). Based on the coordinates of the InSAR device, the coordinate values ​​between the endpoints of the InSAR orbit are measured using a single-point coordinate measurement method.

4. The air-ground combined InSAR geological hazard monitoring method as described in claim 1, characterized in that, Determining the radar deformable azimuth angle specifically includes: Obtain raw InSAR data from InSAR observation points; The cumulative displacement of each InSAR data observation point in the raw InSAR data is obtained by differential interferometry; thus, the deformation data of each InSAR data observation point is obtained. Deformation data acquired at different time points are differentially processed to obtain relative deformation data within the corresponding time period; the radar deformation azimuth angle is calculated based on the deviation between the deformation data and the original InSAR data.

5. The air-ground combined InSAR geological hazard monitoring method as described in claim 4, characterized in that, Calculating the radar deformation azimuth angle based on the deviation between the deformation data and the raw InSAR data specifically includes: The deformation coordinates of the InSAR orbit endpoints are obtained from the deformation data; The original coordinates of the InSAR orbit endpoints are obtained from the raw InSAR data; By connecting the coordinates of the InSAR device with the deformation coordinates of the InSAR orbit endpoint and the original coordinates of the InSAR orbit endpoint, the radar deformation azimuth angle can be obtained.

6. The air-ground combined InSAR geological disaster monitoring method as described in claim 1, characterized in that, Registering InSAR observation data to a combined air-ground 3D terrain model specifically involves: registering the data based on terrain point cloud data, orbital coordinates, and InSAR observation data.

7. The air-ground combined InSAR geological disaster monitoring method as described in claim 1, characterized in that, Also includes: After the initial InSAR stereo detection model is established, the spatial relationship is verified by checking whether the coordinates of the InSAR radar observation data are consistent with the positions on the air-ground joint three-dimensional terrain model. If they are inconsistent, they are re-registered.

8. The air-ground combined InSAR geological disaster monitoring method as described in claim 7, characterized in that, Also includes: Using radar resolution tools, radar interferometric files in DiffImage format are resolved into deformable patches and registered onto the InSAR stereo detection model.

9. The air-ground combined InSAR geological hazard monitoring method as described in claim 1, characterized in that, The first digital orthophoto of the disaster-stricken area was obtained by using a drone. The relative height of the drone to the ground was controlled at 30-40m, the overlap rate of the flight path and the side view was 75%-85%, and the drone's flight area was 50m larger than the disaster-stricken area.

10. The air-ground combined InSAR geological hazard monitoring method as described in claim 1, characterized in that, The image resolution of the first digital orthophoto is less than 5 cm / pix, where pix represents a pixel.

11. The air-ground combined InSAR geological hazard monitoring method as described in claim 1, characterized in that, After acquiring the first digital orthophoto of the disaster-stricken area, the process also includes converting the acquired lidar data into three-dimensional point cloud data.

12. The air-ground combined InSAR geological hazard monitoring method as described in claim 1, characterized in that, The method for constructing the second digital elevation model includes: acquiring lidar data by using a drone equipped with lidar, calculating a three-dimensional point cloud based on the acquired lidar data, and separating the point cloud of the disaster-stricken point by point cloud clipping.

13. The air-ground combined InSAR geological hazard monitoring method as described in claim 12, characterized in that, Point cloud trimming is used to separate the point clouds of disaster-stricken areas, specifically including: By establishing deep learning classification samples of point cloud data, the Patchwork++ algorithm is used to classify and display the point cloud data as ground point cloud data, while the undisplayed part is non-ground point cloud data. By adjusting the deep learning iteration parameters and increasing the deep learning training samples, the Patchwork++ algorithm is optimized to obtain accurate ground point cloud data.

14. An air-ground combined InSAR geological disaster monitoring device, characterized in that, It includes at least one processor and a memory communicatively connected to the at least one processor; the memory stores instructions executable by the at least one processor, which are executed by the at least one processor to enable the at least one processor to implement the air-ground combined InSAR geological hazard monitoring method according to any one of claims 1 to 13.

Images