Spatial-temporal scale landslide disaster induction factor detection method and related device

Abstract

The invention discloses a spatio-temporal scale landslide disaster induction factor detection method and a related device, and relates to the technical field of disaster induction factor detection, and the method comprises the steps: calculating an LOS direction deformation data set of a research region based on an SAR image set, DEM data, a multi-source remote sensing data set and a GACOS atmosphere data set of the research region, taking the LOS direction deformation data set as the input, and carrying out the detection of a landslide disaster induction factor. Calculating a three-dimensional deformation rate data set of the research area in combination with a three-dimensional calculation model and an SPFM model, and detecting a first impact factor from the first type of potential impact factors under a time scale by using a correlation coefficient method based on the three-dimensional deformation rate data set; and detecting a second impact factor from the second type of potential impact factors under the spatial scale by using a geographic detector, and taking the first impact factor and the second impact factor as disaster induction factors. The detection precision of the landslide disaster induction factors can be improved.

Description

Technical Field

[0001] The present application relates to the technical field of disaster inducing factor detection, and in particular to a method and related device for detecting landslide disaster inducing factors on a spatiotemporal scale. Background Art

[0002] The synthetic aperture radar (InSAR) technology developed in recent years can continuously monitor a large-scale research area all day, all weather and with high precision, and has shown great advantages and value in many fields. However, in the later interpretation and analysis of the monitoring results, there is a problem of difficulty in detecting disaster-inducing factors (also called disaster-causing factors), especially under the single LOS (Line-of-sight, radar line of sight) deformation results. With the increasing application scenarios of InSAR technology, the geographical and topographic conditions of the monitored research areas are becoming more and more complex. The single LOS deformation results cannot meet the requirements of disaster-inducing factor detection in some scenarios. It is necessary to invert the three-dimensional deformation of the surface to obtain more accurate disaster information. Therefore, a method of disaster-inducing factor detection that takes into account the quantitative analysis of three-dimensional deformation is needed.

[0003] Domestic and foreign researchers have conducted a lot of research on the detection of disaster inducing factors, and have carried out research on different methods of detecting landslide disaster inducing factors, including statistical index (SI), logistic regression (LR) model, on-site geological observation, random forest model, etc., but the results are not particularly ideal and the operation complexity is relatively high. Geographic detector is a new spatial statistical method for detecting spatial heterogeneity. It can show good quantitative effect on the heterogeneity of factors under spatial scale and can accurately judge the linear and nonlinear interactions between each factor. However, it is difficult for geographic detector to quantify the influence of potential disaster inducing factors under time scale, and it is necessary to supplement the detection method of disaster inducing factors under time scale. Summary of the invention

[0004] The purpose of this application is to provide a method and related device for detecting landslide disaster inducing factors on a spatiotemporal scale, which can solve the problem that a single LOS deformation result cannot meet the requirements for disaster inducing factor detection in some scenarios and is difficult to complete the detection of disaster inducing factors on a time scale, thereby improving the accuracy of landslide disaster inducing factor detection.

[0005] To achieve the above objectives, this application provides the following solutions:

[0006] In a first aspect, the present application provides a method for detecting landslide hazard inducing factors at a spatiotemporal scale, the method comprising:

[0007] Acquire a SAR image set, DEM data, a multi-source remote sensing data set and a GACOS atmospheric data set of the study area; the SAR image set includes a SAR image at each first moment in a preset time period, the multi-source remote sensing data set includes multi-source remote sensing data at each second moment in the preset time period, and the GACOS atmospheric data set includes GACOS atmospheric data at each first moment in the preset time period; the preset time period includes a satellite orbit raising process and a satellite orbit lowering process;

[0008] Based on the SAR image set, DEM data, multi-source remote sensing data set and GACOS atmospheric data set, the LOS deformation data set of the study area is calculated; the LOS deformation data set includes the LOS deformation data of each first moment in a preset time period;

[0009] Taking the LOS deformation data set as input, the three-dimensional deformation rate data set of the study area is calculated in combination with the three-dimensional solution model and the SPFM model; the three-dimensional deformation rate data set includes the three-dimensional deformation rate data at each first moment in the preset time period, the three-dimensional deformation rate data includes the value of the three-dimensional deformation rate of each position point in the study area, and the three-dimensional deformation rate includes the vertical deformation rate, the north-south deformation rate and the east-west deformation rate;

[0010] Based on the three-dimensional deformation rate dataset, the correlation coefficient method is used to detect the first influencing factor from the first type of potential influencing factors on a time scale, and the geographic detector is used to detect the second influencing factor from the second type of potential influencing factors on a spatial scale, and the first influencing factor and the second influencing factor are used as disaster inducing factors; the first type of potential influencing factors include precipitation and temperature, and the second type of potential influencing factors include at least two of DEM data, slope, aspect, curvature, NDVI, land cover type, surface temperature, water system and fault.

[0011] In a second aspect, the present application provides a computer device, comprising: a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement the above-mentioned method for detecting landslide hazard inducing factors at a spatiotemporal scale.

[0012] In a third aspect, the present application provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the above-mentioned method for detecting landslide hazard inducing factors at a spatiotemporal scale.

[0013] In a fourth aspect, the present application provides a computer program product, including a computer program, which, when executed by a processor, implements the above-mentioned method for detecting landslide hazard inducing factors at a spatiotemporal scale.

[0014] According to the specific embodiments provided in this application, this application has the following technical effects:

[0015] The present application provides a method and related device for detecting landslide disaster inducing factors at a spatiotemporal scale. Based on the SAR image set, DEM data, multi-source remote sensing data set and GACOS atmospheric data set of the study area, the LOS deformation data set of the study area is calculated. The LOS deformation data set is used as input, and the three-dimensional deformation rate data set of the study area is calculated in combination with the three-dimensional solution model and the SPFM model. The three-dimensional deformation rate data set is subsequently used to complete the detection of landslide disaster inducing factors. By introducing the three-dimensional deformation rate, the problem that a single LOS deformation result cannot meet the requirements of disaster inducing factor detection in some scenarios can be solved. When using the three-dimensional deformation rate data set to complete the detection of landslide disaster inducing factors, based on the three-dimensional deformation rate data set, the first influencing factor is detected from the first type of potential influencing factors at the time scale using the correlation coefficient method, and the second influencing factor is detected from the second type of potential influencing factors at the spatial scale using the geographic detector, and the first influencing factor and the second influencing factor are used as disaster inducing factors. By introducing the time scale, the spatiotemporal scale can be integrated to complete the detection of landslide disaster inducing factors, solving the problem that it is difficult to complete the detection of disaster inducing factors at the time scale. The present application can improve the accuracy of detecting landslide disaster inducing factors. BRIEF DESCRIPTION OF THE DRAWINGS

[0016] In order to more clearly illustrate the embodiments of the present application or the technical solutions in the prior art, the drawings required for use in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present application. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0017] Figure 1 This is an application environment diagram of a method for detecting landslide hazard inducing factors at a spatiotemporal scale provided in Example 1 of the present application.

[0018] Figure 2 A schematic flow chart of a method for detecting landslide hazard inducing factors at a spatiotemporal scale provided in Example 1 of the present application.

[0019] Figure 3 A detailed flow chart of a method for detecting landslide hazard inducing factors at a spatiotemporal scale provided in Example 1 of the present application.

[0020] Figure 4This is a schematic diagram of SAR image coverage of the southern study area of ​​Dingjie County, Shigatse City, Tibet Autonomous Region provided in Example 1 of the present application.

[0021] Figure 5 A schematic diagram of three-dimensional deformation rate decomposition provided in Example 1 of the present application; wherein, Figure 5 (a) is the east-west deformation rate, Figure 5 (b) is the north-south deformation rate, Figure 5 (c) in the figure is the vertical deformation rate.

[0022] Figure 6 A schematic diagram of the second type of potential influencing factors provided in Example 1 of the present application; wherein, Figure 6 (a) in the figure is DEM (Digital Elevation Model) data. Figure 6 (b) is the slope, Figure 6 (c) in the equation is the slope direction, Figure 6 (d) is the curvature, Figure 6 (e) in the figure is NDVI (Normalized Vegetation Index). Figure 6 (f) in the equation is the land cover type. Figure 6 (g) in the equation is the Land Surface Temperature (LST). Figure 6 (h) in the equation is the water system. Figure 6 (i) in the figure is a fault.

[0023] Figure 7 A schematic diagram of the disaster inducing factor detection results at the spatial scale provided in Example 1 of the present application.

[0024] Figure 8 A schematic diagram of the structure of a computer device provided in Example 2 of the present application. DETAILED DESCRIPTION

[0025] The following will be combined with the drawings in the embodiments of the present application to clearly and completely describe the technical solutions in the embodiments of the present application. Obviously, the described embodiments are only part of the embodiments of the present application, not all of the embodiments. Based on the embodiments in the present application, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of this application.

[0026] Example 1

[0027] The method for detecting landslide disaster inducing factors at a spatiotemporal scale provided in the present application can be applied to Figure 1In the application environment shown. Among them, the terminal communicates with the server through the network. The data storage system can store the data that the server needs to process. The data storage system can be set up separately, integrated on the server, or placed on the cloud or other servers. The terminal can send the detection request to be processed to the server. After the server receives the detection request to be processed, for the detection request to be processed, the server obtains the SAR image set, DEM data, multi-source remote sensing data set and GACOS atmospheric data set of the study area; based on the SAR image set, DEM data, multi-source remote sensing data set and GACOS atmospheric data set, the LOS deformation data set of the study area is calculated; using the LOS deformation data set as input, the three-dimensional deformation rate data set of the study area is calculated in combination with the three-dimensional solution model and the SPFM model; based on the three-dimensional deformation rate data set, the correlation coefficient method is used to detect the first influencing factor from the first type of potential influencing factors at the time scale, and the geographic detector is used to detect the second influencing factor from the second type of potential influencing factors at the spatial scale, and the first influencing factor and the second influencing factor are used as disaster inducing factors. The server can feedback the detection results of the disaster inducing factors for the detection request to the terminal.

[0028] In addition, in some embodiments, the method for detecting landslide hazard inducing factors at a spatiotemporal scale can also be implemented independently by a server or a terminal. For example, the terminal can directly process the detection request to be processed, or the server can obtain the detection request to be processed from the data storage system and process the detection request to be processed.

[0029] The terminals may be, but are not limited to, various desktop computers, laptops, smart phones, tablet computers, IoT devices and portable wearable devices. IoT devices may be smart speakers, smart TVs, smart air conditioners, smart car-mounted devices, etc. Portable wearable devices may be smart watches, smart bracelets, head-mounted devices, etc. The server may be implemented as an independent server or a server cluster consisting of multiple servers, or a cloud server.

[0030] In an exemplary embodiment, Figure 2 and Figure 3 As shown, a method for detecting landslide hazard inducing factors at a spatiotemporal scale is provided. The method is executed by a computer device, and specifically can be executed by a computer device such as a terminal or a server alone, or can be executed by a terminal and a server together. In the embodiment of the present application, the method is applied to Figure 1 The following steps are used to illustrate the server in the example.

[0031] Step S1, obtaining a SAR image set, DEM data, a multi-source remote sensing data set and a GACOS atmospheric data set of a study area; the SAR image set includes a SAR image at each first moment in a preset time period, the multi-source remote sensing data set includes multi-source remote sensing data at each second moment in the preset time period, and the GACOS atmospheric data set includes GACOS atmospheric data at each first moment in the preset time period; the preset time period includes a satellite orbit raising process and a satellite orbit lowering process.

[0032] Step S2, based on the SAR image set, DEM data, multi-source remote sensing data set and GACOS atmospheric data set, calculate the LOS deformation data set of the study area; the LOS deformation data set includes the LOS deformation data of each first moment in a preset time period.

[0033] Step S3, taking the LOS deformation data set as input, combining the three-dimensional solution model and the SPFM model to calculate the three-dimensional deformation rate data set of the study area; the three-dimensional deformation rate data set includes the three-dimensional deformation rate data at each first moment in a preset time period, the three-dimensional deformation rate data includes the value of the three-dimensional deformation rate of each position point in the study area, and the three-dimensional deformation rate includes the vertical deformation rate, the north-south deformation rate and the east-west deformation rate.

[0034] Step S4, based on the three-dimensional deformation rate data set, uses the correlation coefficient method to detect the first influencing factor from the first category of potential influencing factors on a time scale, and uses the geographic detector to detect the second influencing factor from the second category of potential influencing factors on a spatial scale, and uses the first influencing factor and the second influencing factor as disaster inducing factors; the first category of potential influencing factors includes precipitation and temperature, and the second category of potential influencing factors includes at least two of DEM data, slope, aspect, curvature, NDVI, land cover type, surface temperature, water system and fault.

[0035] The above steps S1 to S4 are implemented. After determining the LOS deformation dataset of the study area, this embodiment uses the LOS deformation dataset as input, combines the three-dimensional solution model and the SPFM model to calculate the three-dimensional deformation rate dataset of the study area, thereby inverting the three-dimensional deformation of the surface to obtain more accurate disaster information. Subsequently, the disaster inducing factor detection is completed based on the quantitative analysis of the three-dimensional deformation, so as to solve the problem that the single LOS deformation result cannot meet the requirements of disaster inducing factor detection in some scenarios. When completing the disaster inducing factor detection based on the three-dimensional deformation rate dataset, based on the three-dimensional deformation rate dataset, the correlation coefficient method is used to detect the first influencing factor from the first type of potential influencing factors at the time scale, and the geographic detector is used to detect the second influencing factor from the second type of potential influencing factors at the spatial scale, and the first influencing factor and the second influencing factor are used as disaster inducing factors, so as to detect the disaster inducing factors at the time and space scale, and solve the problem that it is difficult to quantify the degree of influence of the disaster inducing factors at the time scale.

[0036] In this embodiment, a SAR image set, DEM data, a multi-source remote sensing data set, and a GACOS atmospheric data set are obtained for the study area. The SAR image set may be an L1-level time series product obtained by a spaceborne synthetic aperture radar covering the study area, including a SAR image at each first moment in a preset time period. The SAR image covers the study area, and the SAR image obtained by Sentinel-1 covering the southern study area of ​​Dingjie County, Shigatse City, Tibet Autonomous Region is as follows: Figure 4 As shown, the multi-source remote sensing data set can be a multi-spectral time series data set of the study area obtained by using sensors such as Sentinel-2 MSI, MODIS, and Landsat satellite MSS, including multi-source remote sensing data at each second moment in a preset time period, and the multi-source remote sensing data covers the study area. The first moment and the second moment can be the same or different. According to the acquisition time of each SAR image, the GACOS atmospheric data covering the study area is downloaded. The GACOS atmospheric data set includes the GACOS atmospheric data at each first moment in the preset time period, and the GACOS atmospheric data covers the study area. The pixel value of the pixel point of the GACOS atmospheric data is the atmospheric delay distance. The preset time period includes the satellite orbit raising process and the satellite orbit lowering process, and the SAR image set includes multiple SAR images of the satellite orbit raising process and multiple SAR images of the satellite orbit lowering process.

[0037] In this embodiment, based on the SAR image set, DEM data, multi-source remote sensing data set and GACOS atmospheric data set, the LOS deformation data set of the study area is calculated. The LOS deformation data set includes the LOS deformation data at each first moment in a preset time period. The LOS deformation data covers the study area, including the value of the LOS deformation at each location point in the study area.

[0038] Among them, based on SAR image sets, DEM data, multi-source remote sensing data sets and GACOS atmospheric data sets, the LOS deformation data sets of the study area are calculated, including:

[0039] (1) Based on the DEM data, the SAR images in the SAR image set are registered to obtain a registered image dataset of the study area. The registered image dataset includes the registered SAR images at each first moment in a preset time period, and the registered SAR images cover the study area.

[0040] This embodiment combines external DEM data to assist the image registration process of the SAR images in the SAR image set, and obtains the registered single-view complex map (i.e., the registered SAR image). In practical applications, the SAR image set acquired by Sentinel-1 is used, and the external DEM data is combined to assist the image registration process, and the registered single-view complex map is obtained. The "three-step registration method" is mainly used to register the SAR image set in the TOPS mode, which specifically includes: the first step is based on DEM data, and the geometric registration method is used for rough registration; the second step is to use the intensity cross-correlation method (which calculates the correlation between intensities) for further registration; the third step is to use the enhanced spectral diversity method for fine registration. The three-step registration method can achieve a registration accuracy of one thousandth, and obtain a registered image data set.

[0041] (2) Geocoding the registered image dataset and DEM data to obtain a detailed lookup table for coordinate conversion between the radar coordinate system and the geographic coordinate system.

[0042] In this embodiment, the registered image dataset and the DEM data are geocoded to obtain a refined geographic-SAR lookup table (i.e., a refined lookup table) for coordinate conversion between the geographic coordinate system and the radar coordinate system. In practical applications, DEM data having a resolution substantially consistent with that of the SAR image is introduced, and the simulated SAR image at each first moment in a preset time period is generated by combining the DEM data and the satellite (satellite used to obtain the SAR image) orbit data. The simulated SAR image and the real SAR image (i.e., the registered SAR image in the registered image dataset) are then registered to establish a registration relationship between the two, and the conversion of the radar image from the slant range projection to the orthographic projection is realized to obtain a refined lookup table. The refined lookup table is mainly obtained by using a "two-step registration method", which specifically includes: the first step is to generate a simulated SAR image using the DEM data and the satellite orbit data of the SAR imaging, and to combine the SAR imaging model, and to roughly register the simulated SAR image and the real SAR image, and to estimate the initial offset to obtain an initial lookup table. The initial lookup table includes the corresponding relationship between the SAR imaging coordinate system (i.e., the radar coordinate system) and the WGS84 coordinate system (i.e., the geographic coordinate system). r, a are the pixel numbers in the range and azimuth directions in the radar coordinate system, respectively; g is the correspondence between the radar coordinate system and the geographic coordinate system. They are the latitude, longitude and elevation in the geographic coordinate system respectively; the second step is to use the intensity cross-correlation method to cross-correlate the simulated SAR image and the real SAR image, and then estimate the offset of the simulated SAR image and the real SAR image according to the matching results and the initial lookup table, that is, to select the control points with the same name in the simulated SAR image and the real SAR image respectively, and to fit the control points with the same name using polynomials, to calculate the offset polynomial and to obtain the refined lookup table.

[0043] When refining the offset estimation of the simulated SAR image and the real SAR image according to the intensity cross-correlation method, the registration window is set to 128*128 and the offset polynomial order is set to 4.

[0044] (3) A remote sensing ecological index dataset of the study area is calculated using a multi-source remote sensing dataset. The remote sensing ecological index dataset includes remote sensing ecological index data at every second moment in a preset time period. The remote sensing ecological index data covers the study area, including the value of the remote sensing ecological index of each location point in the study area.

[0045] This embodiment is based on a multi-source remote sensing data set and uses the GEE (Google Earth Engine) platform to invert remote sensing ecological indices, including but not limited to NDVI, LST, NDSI (Normalized differential snow index), etc. Of course, other remote sensing ecological indices of interest can also be selected to obtain a remote sensing ecological index data set.

[0046] This embodiment also uses the NDSI index to extract the glacier area, where the extraction formula is:

[0047] GlacierArea= <Ave{Image∈S2:[(D 1 ≥D 1 * )∧(D 2 ≤D 2 * )∧(NDSI≥NDSI * )]}≥0.5>;

[0048] Where GlacierArea is the glacier area; Ave represents the calculated average value; S2 is the multi-source remote sensing dataset acquired by Sentinel-2 covering the study area during a preset time period (e.g., from January 1, 2021 to September 18, 2023); D 1 is the first index; * indicates the threshold value calculated by OSTU (maximum between-class variance method), is the threshold of the first index calculated by OSTU; ∧ indicates that; D 2 is the second index; is the threshold of the second index calculated by OSTU; NDSI is the normalized difference snow index; NDSI * is the threshold of the normalized difference snow index calculated by OSTU. When the first index is greater than or equal to the threshold of the first index, the second index is greater than or equal to the threshold of the second index, and the normalized difference snow index is greater than or equal to the threshold of the normalized difference snow index, the pixel value of the pixel is 1, indicating that the pixel belongs to the glacier area.

[0049]

[0050] Among them, Red is the reflectivity of the red band; SWIR 1 is the reflectivity in the infrared band.

[0051]

[0052] Among them, Blue is the reflectivity of the blue band.

[0053]

[0054] Among them, Green is the reflectivity of the green band.

[0055] (4) Based on the DEM data, the remote sensing ecological index dataset and the GACOS atmospheric dataset are resampled to obtain a resampled index dataset and a resampled atmospheric dataset. The resampled index dataset includes the resampled index data at each second moment in the preset time period, and the resampled atmospheric dataset includes the resampled atmospheric data at each first moment in the preset time period. Both the resampled index data and the resampled atmospheric data cover the study area.

[0056] For the remote sensing ecological index dataset and GACOS atmospheric dataset of the study area, resample them according to the geocoded DEM data to obtain the range, corner latitude and longitude, resolution and other parameters of the geocoded DEM data. Resample the remote sensing ecological index data and GACOS atmospheric data to the above geographic parameter configuration to obtain multi-source remote sensing products consistent with the geographic parameter configuration of the geocoded DEM data. Specifically, first calculate the longitude and latitude of the center point of each grid in the DEM data according to the range, corner latitude and longitude, resolution and other parameters of the geocoded DEM data. The calculation formula is as follows:

[0057] Lat i =Lat min +(i+0.5)×ΔLat,i=0,1,2,...,Nlat -1;

[0058] Lon j =Lon max -(j+0.5)×ΔLon,j=0,1,2,...,N lon -1;

[0059] Among them, Lat i is the latitude of the center point of the grid (i, j); Lat min is the minimum latitude of the grid; ΔLat is the pixel latitude interval; N lat is the number of grid rows; Lon j is the longitude of the center point of the grid (i, j); Lon max is the maximum longitude of the grid; ΔLon is the pixel longitude interval; N lon is the number of grid columns.

[0060] Then, the nearest neighbor, bilinear, B-spline and other resampling methods are used to resample the remote sensing ecological index dataset and the GACOS atmospheric dataset according to the resampling formula, and the remote sensing ecological index dataset and the GACOS atmospheric dataset are unified to the geographical coordinate system and parameter configuration of the geocoded DEM data. The resampling formula is:

[0061] Z′(x′,y′)=R(Z(x,y));

[0062] Wherein, Z′(x′, y′) is the pixel value of the pixel point (x′, y′) of the expected resampled data; R is the sampling method; and Z(x, y) is the pixel value of the pixel point (x, y) of the original data.

[0063] (5) Based on the refined lookup table, the DEM data, the resampled index data set, and the resampled atmosphere data set are backward geocoded, and the DEM data, the resampled index data set, and the resampled atmosphere data set in the geographic coordinate system are mapped to the radar coordinate system to obtain the transformed DEM data, the transformed index data set, and the transformed atmosphere data set in the radar coordinate system.

[0064] When mapping the data in the geographic coordinate system to the radar coordinate system according to the fine lookup table, the mapping formula used is:

[0065]

[0066] in, is the pixel coordinate of a certain object in the radar coordinate system B; Interp is the convolution interpolation method: f = (∑ i w i ·f i ) / ∑ i wi , w i is the weight of the i-th neighboring point, f i is the pixel value of the ith neighboring point. The neighboring point is the pixel point in the neighborhood of the pixel point where a certain feature is located. The neighborhood can be determined according to needs, such as 8 neighborhoods. The weight can be determined by the distance weighted method, that is, w i =1 / d i , d i L is the distance from the neighboring point to the pixel where a certain feature is located; x (),L y () are the mapping relationships between geographic coordinate system A and radar coordinate system B in rows and columns, and (xA, yA) are the pixel coordinates of the same object in geographic coordinate system A.

[0067] (6) The registered SAR images in the registered image data set are combined to obtain a plurality of differential interferometry pairs, and differential interferometry is performed on each differential interferometry pair to obtain a plurality of differential interferograms.

[0068] This embodiment combines differential interferometry pairs based on the idea of ​​a small baseline set. The idea of ​​a small baseline set is to reduce the spatial decorrelation effect caused by a large baseline by selecting a group of interferograms with a small baseline in time and space. Specifically, the differential interferometry pairs are combined in a combination manner of multiple main images under the restriction that the time baseline is less than 90 days and the space baseline is less than 400 meters. That is, each registered SAR image in the registered image data set is used as the main image, and the registered SAR images in the registered image data set with a time baseline less than 90 days and a space baseline less than 400 meters with the main image are selected as selected images, and the main image and each selected image are respectively combined into a differential interferometry pair.

[0069] Then, for each differential interferometer pair, differential interferometry is performed on the two registered SAR images included in the differential interferometer pair to obtain a differential interferogram.

[0070] (7) The undulating phase and the flat ground phase are calculated using the converted DEM data, and the undulating phase and the flat ground phase are removed from each differential interferogram to obtain a plurality of differential interferograms after preliminary removal.

[0071] According to the converted DEM data and satellite orbit parameters in the radar coordinate system generated by geocoding, the terrain undulation phase and flat ground phase are simulated, and the above two terrain-related undulation phases and flat ground phases are subtracted from each differential interferogram to preliminarily remove the terrain-related phase, and obtain multiple differential interferograms after preliminary removal. The number of the preliminary removed differential interferograms is the same as that of the differential interferograms.

[0072] This embodiment also introduces Goldstein phase optimization, uses Goldstein filter to sequentially set 64*64, 32*32, and 16*16 window sizes, iteratively performs phase optimization on the generated differential interferogram, and weakens the coherent speckle noise in the differential interferogram.

[0073] At this time, the undulation phase and the flat ground phase are calculated using the converted DEM data, and the undulation phase and the flat ground phase are removed in each differential interferogram to obtain multiple preliminary differential interferograms after removal, specifically including: using the converted DEM data to calculate the undulation phase and the flat ground phase; for each differential interferogram, using the Goldstein filter to perform three phase optimizations on the differential interferogram to obtain a denoised differential interferogram, wherein the windows used for the three phase optimizations are 64*64, 32*32 and 16*16 respectively; removing the undulation phase and the flat ground phase in each denoised differential interferogram to obtain multiple preliminary differential interferograms after removal.

[0074] (8) Based on the converted index data set and multiple differential interference patterns after preliminary removal, multiple high coherence points are determined, and based on all the high coherence points, a high coherence point is selected as a reference point.

[0075] The step of determining a plurality of high coherence points specifically includes: determining a plurality of high coherence points using an amplitude deviation method and a coherence coefficient method.

[0076] Generally, the area of ​​interest is obtained through a threshold based on the remote sensing ecological index of interest. In this embodiment, the glacier area extracted by the NDSI index is further combined with the overlap, shadow mask and glacier area to generate a mask file. The mask file includes pixel points belonging to the overlap, shadow and glacier areas. Based on the mask file, that is, in the remaining pixel points after removing the pixel points in the mask file, the high coherence points of interest in the differential interference map after all preliminary removals are selected through the amplitude deviation method and the coherence coefficient method. The high coherence points can be high coherence points that meet the preset amplitude deviation threshold at the same time. The intersection of the point and the high coherence point selected by satisfying the preset coherence coefficient threshold may also be the union of the high coherence point selected by satisfying the preset amplitude deviation threshold and the high coherence point selected by satisfying the preset coherence coefficient threshold. The high coherence point selected by satisfying the preset amplitude deviation threshold means that the amplitude deviation of the pixel point is greater than the preset amplitude deviation threshold, and the high coherence point selected by satisfying the preset coherence coefficient threshold means that the coherence coefficient of the pixel point is greater than the preset coherence coefficient threshold, so as to obtain the high coherence point by adopting the amplitude deviation and the coherence coefficient, and the preset amplitude deviation threshold and the preset coherence coefficient threshold are 0.15 and 0.8 respectively.

[0077] Amplitude deviation D A =The calculation formula is:

[0078]

[0079] Among them, σ A and m A They represent the standard deviation and mean of the amplitude of the pixel in the time series respectively.

[0080] The calculation formula of the coherence coefficient is:

[0081]

[0082] Wherein, γ is the coherence coefficient; m is the number of image rows; n is the number of image columns; M and S are the signals of high coherence points in the first scene registered SAR image and the second scene registered SAR image, respectively. The first scene registered SAR image and the second scene registered SAR image form a differential interferometer pair; (i, j) represents the row and column number of the high coherence point; * represents conjugate multiplication.

[0083] (9) For each differential interferogram after preliminary removal, phase unwrapping is performed on each high coherence point based on the reference point to obtain the unwrapped phase of each high coherence point, and the unwrapped phase of each high coherence point is corrected by iterative regression analysis to remove the terrain phase from the differential interferogram after preliminary removal to obtain a differential interferogram after secondary removal.

[0084] Wherein, phase unwrapping is performed on each high coherence point based on the reference point, specifically including: based on the reference point, phase unwrapping is performed on each high coherence point using the minimum cost flow method.

[0085] In this embodiment, a reference point is selected, and a reference point is selected from all high coherence points according to the central principle, the stability principle, the high coherence principle, and the non-island area principle; for each high coherence point, based on the reference point, the minimum cost flow method of the space dimension is used to perform phase unwrapping on the high coherence point, and the unwrapped phase of the high coherence point is obtained, and the unwrapped phase of the high coherence point is corrected by iterative regression analysis until the elevation correction value of the iterative result is Less than 0.1m, in order to remove the terrain phase, to obtain the differential interferogram after the second removal, the number of the differential interferograms after the second removal is the same as the number of the differential interferograms after the initial removal.

[0086] This embodiment specifically adopts an interference phase model for iterative regression analysis and correction. The interference phase model is:

[0087]

[0088] Among them, δφ i,j is the unwrapping phase of the high coherence point; λ is the wavelength; v j is the deformation rate; Δt j =t 2 -t 1 ,t 2>t 1 , t 2 is the time of the second scene registered SAR image, t 1 is the time of the first scene registered SAR image, the first scene registered SAR image and the second scene registered SAR image form a differential interferometer pair; δ i,j is the deformation residual; B j is the vertical baseline; Δz i is the DEM height error; ρ i is the slope distance; θ i is the angle of incidence; is the atmospheric delay phase; t 2,j is the time corresponding to the high coherence point in the SAR image after the second scene registration; t 1,j is the time corresponding to the high coherence point in the SAR image after the first scene registration; n i is the thermal phase noise.

[0089] (10) For each second-removed differential interferogram, based on the time of the differential interference pair corresponding to the second-removed differential interferogram, two converted atmospheric data at the corresponding time are selected from the converted atmospheric data set, and the atmospheric delay phase is calculated based on the two converted atmospheric data, and the atmospheric delay phase is removed from the second-removed differential interferogram to obtain the removed differential interferogram.

[0090] This embodiment performs differential calculation on the GACOS atmospheric data according to the differential interference pair combination, and calculates the relative value of each GACOS atmospheric data based on the reference point. Specifically, the atmospheric delay phase of the reference point is subtracted from the differential GACOS atmospheric data, that is, The purpose of this is to unify the atmospheric delay phase benchmark of each GACOS atmospheric data and convert the GACOS atmospheric data into a two-way phase, namely S M ·4π / λ, the obtained atmospheric delay phase is removed from the optimized unwrapped phase, so as to remove the atmospheric delay phase from the optimized unwrapped phase and obtain a differential interferogram after removal.

[0091] The calculation formula of the atmospheric delay phase is:

[0092]

[0093] in, is the atmospheric delay phase; is the Mth converted atmospheric data S M The difference data with the Nth converted atmospheric data, the Mth converted atmospheric data and the Nth converted atmospheric data are two converted atmospheric data; is the atmospheric delay phase of the difference data between the Mth converted atmospheric data and the Nth converted atmospheric data at the reference point (i, j); λ is the wavelength of the SAR sensor (i.e., the sensor used to acquire SAR images).

[0094] (11) Based on all the removed differential interference patterns, the LOS deformation dataset of the study area is calculated.

[0095] In this embodiment, the LOS deformation data set is used as input, and the three-dimensional deformation rate data set of the study area is calculated in combination with the three-dimensional solution model and the SPFM model. The three-dimensional deformation rate data set includes the three-dimensional deformation rate data at each first moment in a preset time period. The three-dimensional deformation rate data includes the value of the three-dimensional deformation rate of each position point in the study area. The three-dimensional deformation rate includes the vertical deformation rate, the north-south deformation rate, and the east-west deformation rate.

[0096] Since the SAR image set includes an orbit raising dataset and an orbit descending dataset, the above processing is performed on the orbit raising dataset and the orbit descending dataset respectively to obtain the temporal deformation in the LOS direction of the orbit raising process and the orbit descending process.

[0097] This embodiment combines the LOS deformation of the lifting rail, introduces the SPFM model, and combines the three-dimensional solution model and the SPFM model to solve the three-dimensional deformation rate results of the study area to obtain the deformation rates in the vertical, north-south, and east-west directions. The three-dimensional solution model is:

[0098] D Los =D v cosθ+D N sinθsinβ-D E sinθcosβ;

[0099] Among them, D Los is the deformation in the LOS direction; D v is the vertical deformation; θ is the radar incident angle; D N is the north-south deformation; β is the radar heading angle; D E It is an east-west deformation.

[0100] The SPFM model is:

[0101]

[0102] Among them, v U is the vertical deformation rate; is the first derivative to the north, H is the terrain height, X N is the northward deformation component; ν N is the north-south deformation rate; is the first derivative in the east direction, X E is the eastward deformation component; vE is the east-west deformation rate.

[0103] Taking the southern study area of ​​Dingjie County, Shigatse City, Tibet Autonomous Region as an example, the generated three-dimensional deformation rate decomposition results are shown in Figure 5 .

[0104] In this embodiment, based on the three-dimensional deformation rate data set, the correlation coefficient method is used to detect the first influencing factor from the first type of potential influencing factors at the time scale, and the geographic detector is used to detect the second influencing factor from the second type of potential influencing factors at the spatial scale, and the first influencing factor and the second influencing factor are used as disaster inducing factors. Among them, the first type of potential influencing factors include precipitation and temperature, and the second type of potential influencing factors include at least two of DEM data, slope, aspect, curvature, NDVI, land cover type, surface temperature, water system and fault. Taking the southern study area of ​​Dingjie County, Shigatse City, Tibet Autonomous Region as an example, the second type of potential influencing factors can be found in Figure 6 .

[0105] This embodiment integrates geographic detectors and PCA analysis (Principal Component Analysis) to detect disaster inducing factors from a temporal and spatial scale, specifically including: using geographic detectors to detect the explanatory power of the second type of potential influencing factors on deformation results, and calculating the relative q value (q / p), where the q value represents the explanatory power of one variable on another variable, and the p value represents the confidence level of this explanatory power; using PCA to extract the principal components of the vertical deformation rate in the study area, and detecting the significant influence level of the first type of potential influencing factors such as temperature and precipitation in the study area on the principal components of landslide disasters in the study area.

[0106] Among them, based on the three-dimensional deformation rate data set, the correlation coefficient method is used to detect the first influencing factor from the first type of potential influencing factors at the time scale, and the geographic detector is used to detect the second influencing factor from the second type of potential influencing factors at the spatial scale, specifically including:

[0107] (1) Obtain a factor data set of each first-category potential influencing factor and a value of each second-category potential influencing factor in the study area. The factor data set includes the value of the first-category potential influencing factor at each first moment in a preset time period. The value of the first-category potential influencing factor is the value of the first-category potential influencing factor at each location point in the study area, and the value of the second-category potential influencing factor is the value of the second-category potential influencing factor at each location point in the study area.

[0108] The factor data set can be obtained by interpolation. The precipitation data inverted by the GEE platform is selected according to the acquisition time of the SAR image, and the temperature data is interpolated using the B-spline method. The temperature data corresponding to the acquisition time of the SAR image is selected from the interpolated full-time temperature data to form a new precipitation data set and temperature data set.

[0109] Among them, obtaining a factor dataset of each first-category potential influencing factor and a value of each second-category potential influencing factor in the study area specifically includes: using a multi-source remote sensing dataset to extract the glacier area in the study area, and removing the glacier area in the study area, obtaining a factor dataset of each first-category potential influencing factor and a value of each second-category potential influencing factor in the non-glacier area of ​​the study area.

[0110] (2) Based on the three-dimensional deformation rate data set, a vertical deformation rate data set is determined. The vertical deformation rate data set includes the vertical deformation rate data at each first moment in a preset time period. The vertical deformation rate data includes the value of the vertical deformation rate at each location point in the study area (after removing the glacier area).

[0111] (3) Taking the factor dataset and vertical deformation rate dataset of each first-category potential influencing factor as input, the correlation coefficient method is used to detect the first influencing factor from the first-category potential influencing factors on a time scale.

[0112] (4) Taking the value of each second-category potential influencing factor and the value of the vertical deformation rate as input, the geographic detector is used to detect the second influencing factor from the second-category potential influencing factors at a spatial scale. The value of the vertical deformation rate is a vertical deformation rate data at a random first moment in the vertical deformation rate data set or the average value of each vertical deformation rate data at the first moment in the vertical deformation rate data set.

[0113] At this time, the factor data set and vertical deformation rate data set of each first-class potential influencing factor are used as input, and the first influencing factor is detected from the first-class potential influencing factors on a time scale using the correlation coefficient method, specifically including:

[0114] (1) The principal component analysis method is used to analyze the vertical deformation rate data set to obtain a principal component data set corresponding to each of the multiple principal components. The principal component data set includes the value of the principal component at each first moment in a preset time period. The value of the principal component is the value of the principal component at each location point in the study area (after removing the glacier area).

[0115] This embodiment uses PCA to extract the principal components of the vertical deformation rate data set in the study area, obtains multiple principal components with a cumulative probability of 98%, and participates in the detection and calculation of disaster inducing factors on a time scale.

[0116] (2) For each first-category potential influencing factor, the factor data set corresponding to the first-category potential influencing factor and the principal component data set corresponding to each principal component are used as input to calculate the correlation coefficient between the first-category potential influencing factor and each principal component.

[0117] This embodiment uses the Pearson correlation coefficient to detect and calculate the quantitative impact evaluation index of precipitation, temperature and other time series data on different principal components. The calculation formula is:

[0118]

[0119] Among them, r is the correlation coefficient; n is the number of samples, that is, the number of the first moment in the preset time period; X i is the i-th sample X; is the mean of sample X; Y i is the i-th sample Y; is the mean of sample Y; sample X and sample Y represent the principal components and the first type of potential influencing factors respectively.

[0120] (3) Detect the first influencing factor from the first category of potential influencing factors based on all correlation coefficients.

[0121] All correlation coefficients are sorted in descending order of correlation, and the first type of potential influencing factors ranked in the top N are extracted as the first influencing factors.

[0122] The calculation results of the first type of potential impact factors are shown in Table 1.

[0123] Table 1. Principal component significance test table of time scale disaster inducing factor detection

[0124]

[0125] In Table 1, PC1, PC2, PC3, PC4, and PC5 represent the five principal components determined, respectively; * represents a relatively high correlation; ** represents a significant correlation.

[0126] The use of geographic detectors to complete the detection process is an existing mature technology. Only some of the steps are introduced here: the second type of potential influencing factors are reclassified according to the parameter optimal principle, and the absolute value deviation of the vertical deformation rate in the spatial dimension is calculated. Specifically, the second type of potential influencing factors are reclassified and the quantile discretization method is used according to the parameter optimal principle: P k =n·k / m,P k is the position of the kth quantile, n is the total number of data, k is the quantile number, k = 1, 2, ..., m, m is the total number of quantiles, and the absolute value deviation of the vertical deformation rate in the spatial dimension is calculated: x i is the vertical deformation rate of the ith pixel, assuming that the result follows a standard normal distribution,

[0127] The geographic detector is used to detect the explanatory power of the second type of potential influencing factors on the deformation results, and the relative q value (q / p) is calculated. The q value represents the explanatory power of one variable on another variable, and the p value represents the confidence level of this explanatory power. When the confidence level represented by the p value is guaranteed to be below 0.3, the corresponding q value and p value are involved in the calculation. The calculation formula of the q value is:

[0128]

[0129] SST=Nσ 2 ;

[0130] Wherein, h=1,...,L is the stratification, classification or partition of factor X (i.e., one or more second-category potential influencing factors); N h and N are the number of units in layer h and the whole area, respectively; and σ 2 are the variances of the Y values ​​of layer h and the entire area (i.e., the absolute value deviation of the vertical deformation rate), respectively; SSW and SST are the sum of the variances within the layer and the total variance of the entire area, respectively.

[0131] The calculation results of the second impact factor detected by the geographic detector can be found in Figure 7 .

[0132] The present embodiment discloses a method for detecting landslide disaster inducing factors in high-cold areas at a spatiotemporal scale by fusing multi-source remote sensing data, which relates to the field of microwave remote sensing measurement technology. The method is a method for detecting disaster inducing factors that takes into account three-dimensional deformation and is performed at a spatiotemporal scale, and improves the detection accuracy of disaster inducing factors. The method comprises: obtaining L1-level time series products of a satellite-borne synthetic aperture radar covering a study area, and combining external DEM data to assist in image registration to obtain a registered single-view complex map; geocoding the DEM to SAR image to obtain a refined geo-SAR lookup table of the geographic coordinate system and the radar coordinate system to ensure that the geocoding accuracy is better than 0.5 Pixel; utilizing the GEE platform to invert remote sensing ecological indices based on multi-source remote sensing data, including but not limited to NDVI, LST, NDSI, etc., and simultaneously obtaining GACOS atmospheric data of the study area, resampling the remote sensing ecological index and the GACOS atmospheric data according to the geocoded DEM data, and obtaining a multi-source remote sensing product with a consistent geographic parameter configuration with the geocoded DEM data; and performing geocoding on the remote sensing ecological index of interest according to the refined geo-SAR lookup table. Backward geocoding is performed with GACOS atmospheric data to obtain the remote sensing ecological index of interest and GACOS atmospheric data in the radar coordinate system; differential interference pairs are combined based on the idea of ​​small baseline sets, and differential interference is performed on the registered single-view complex image to preliminarily remove the terrain-related phase, and iterative Goldstein phase optimization is performed; mask files are generated by combining overlapping masks, shadow masks and remote sensing ecological indexes of interest under appropriate thresholds, and based on the mask files, high coherence points of interest in the differential interferogram are selected by the amplitude deviation method and the coherence coefficient method; reference points are selected, and the unwrapped phase of high coherence points is solved by the minimum cost flow method, and regression analysis is iteratively performed to correct and remove the terrain phase; the relative value of each GACOS atmospheric data is calculated based on the reference point and converted into a two-way phase, and the obtained atmospheric delay phase is removed from the optimized unwrapped phase; the SPFM model is introduced to solve the three-dimensional deformation rate results of the study area in combination with the LOS deformation of the ascending and descending orbits; according to the obtained three-dimensional deformation rate results in the spatiotemporal dimension, the geographical detector in the spatial dimension and the PCA analysis in the temporal dimension are integrated to detect disaster-inducing factors from the spatiotemporal scale. The method for detecting landslide disaster causes in high-cold areas at a spatiotemporal scale by fusing multi-source remote sensing data provided in this embodiment solves the problem that disaster inducing factors are difficult to detect due to LOS deformation, supplements the method for detecting disaster inducing factors at a time scale, and improves the level of landslide disaster monitoring and analysis.

[0133] The method for detecting landslide disaster causes in high-cold areas at a spatiotemporal scale by fusing multi-source remote sensing data provided in this embodiment first uses GEE to establish a potential disaster factor database, and fuses SAR and multispectral remote sensing images to select areas of interest; at the same time, taking into account the situation that a single LOS deformation result cannot reflect the actual surface deformation, the SPFM model is introduced to invert the three-dimensional surface deformation, and the spatiotemporal deformation data of the study area are obtained. Then, the geographic detector and PCA analysis method are integrated to detect disaster inducing factors in the spatiotemporal dimension, and a comprehensive disaster factor quantitative index system in the spatiotemporal dimension is obtained, which effectively quantifies the geological disaster disaster factor indicators and supplements the deficiencies of related work.

[0134] The present application also provides an application scenario, which applies the above-mentioned method for detecting landslide disaster inducing factors at a spatiotemporal scale. Specifically, the method for detecting landslide disaster inducing factors at a spatiotemporal scale provided in this embodiment can be applied in a landslide disaster interpretation scenario. The landslide disaster interpretation scenario includes a detection link and an interpretation link. The detection link is used to integrate geographic detectors and principal component analysis technology to detect disaster inducing factors at a spatiotemporal scale, and establish a landslide disaster influencing factor system in the study area. The interpretation link is used to interpret landslides based on the landslide disaster influencing factor system. The method for detecting landslide disaster inducing factors at a spatiotemporal scale provided in this embodiment belongs to the detection link.

[0135] Example 2

[0136] In an exemplary embodiment, a computer device is provided. The computer device may be a server or a terminal. The internal structure diagram thereof may be as follows: Figure 8 As shown. The computer device includes a processor, a memory, an input / output interface (Input / Output, referred to as I / O) and a communication interface. The processor, the memory and the input / output interface are connected through a system bus, and the communication interface is connected to the system bus through the input / output interface. The processor of the computer device is used to provide computing and control capabilities. The memory of the computer device includes a non-volatile storage medium and an internal memory. The non-volatile storage medium stores an operating system, a computer program and a database. The internal memory provides an environment for the operation of the operating system and the computer program in the non-volatile storage medium. The database of the computer device is used to store data. The input / output interface of the computer device is used to exchange information between the processor and an external device. The communication interface of the computer device is used to communicate with an external terminal through a network connection. When the computer program is executed by the processor, a method for detecting landslide hazard inducing factors at a spatiotemporal scale is implemented.

[0137] Those skilled in the art will understand that Figure 8The structure shown in the figure is only a block diagram of a part of the structure related to the solution of the present application, and does not constitute a limitation on the computer device to which the solution of the present application is applied. The specific computer device may include more or fewer components than those shown in the figure, or combine certain components, or have a different arrangement of components.

[0138] In an exemplary embodiment, a computer device is provided, including a memory and a processor, wherein a computer program is stored in the memory, and when the processor executes the computer program, the method for detecting landslide hazard inducing factors at a spatiotemporal scale in Example 1 is implemented.

[0139] Example 3

[0140] In an exemplary embodiment, a computer-readable storage medium is provided, which stores a computer program. When the computer program is executed by a processor, the method for detecting landslide hazard inducing factors at a spatiotemporal scale in Example 1 is implemented.

[0141] Example 4

[0142] In an exemplary embodiment, a computer program product is provided, including a computer program, which implements the method for detecting landslide hazard inducing factors at a spatiotemporal scale in Example 1 when executed by a processor.

[0143] It should be noted that the user information (including but not limited to user device information, user personal information, etc.) and data (including but not limited to data used for analysis, stored data, displayed data, etc.) involved in this application are all information and data authorized by the user or fully authorized by all parties, and the collection, use and processing of relevant data must comply with relevant regulations.

[0144] The technical features of the above embodiments may be combined arbitrarily. To make the description concise, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0145] This article uses specific examples to illustrate the principles and implementation methods of this application. The description of the above embodiments is only used to help understand the method and core ideas of this application. At the same time, for those skilled in the art, according to the ideas of this application, there will be changes in the specific implementation methods and application scope. In summary, the content of this specification should not be understood as limiting this application.

Claims

1. A method for detecting landslide hazard inducing factors at a spatiotemporal scale, characterized in that: The method for detecting landslide disaster inducing factors at a temporal and spatial scale includes: Acquire a SAR image set, DEM data, a multi-source remote sensing data set and a GACOS atmospheric data set of the study area; the SAR image set includes a SAR image at each first moment in a preset time period, the multi-source remote sensing data set includes multi-source remote sensing data at each second moment in the preset time period, and the GACOS atmospheric data set includes GACOS atmospheric data at each first moment in the preset time period; the preset time period includes a satellite orbit raising process and a satellite orbit lowering process; Based on the SAR image set, DEM data, multi-source remote sensing data set and GACOS atmospheric data set, the LOS deformation data set of the study area is calculated; the LOS deformation data set includes the LOS deformation data of each first moment in a preset time period; Taking the LOS deformation data set as input, the three-dimensional deformation rate data set of the study area is calculated in combination with the three-dimensional solution model and the SPFM model; the three-dimensional deformation rate data set includes the three-dimensional deformation rate data at each first moment in the preset time period, the three-dimensional deformation rate data includes the value of the three-dimensional deformation rate of each position point in the study area, and the three-dimensional deformation rate includes the vertical deformation rate, the north-south deformation rate and the east-west deformation rate; Based on the three-dimensional deformation rate dataset, the correlation coefficient method is used to detect the first influencing factor from the first type of potential influencing factors on a time scale, and the geographic detector is used to detect the second influencing factor from the second type of potential influencing factors on a spatial scale, and the first influencing factor and the second influencing factor are used as disaster inducing factors; the first type of potential influencing factors include precipitation and temperature, and the second type of potential influencing factors include at least two of DEM data, slope, aspect, curvature, NDVI, land cover type, surface temperature, water system and fault.

2. The method for detecting landslide hazard inducing factors at a spatiotemporal scale according to claim 1, characterized in that: Based on SAR image sets, DEM data, multi-source remote sensing data sets and GACOS atmospheric data sets, the LOS deformation data sets of the study area are calculated, including: Based on the DEM data, the SAR images in the SAR image set are registered to obtain a registered image dataset of the study area; the registered image dataset includes the registered SAR images at each first moment in a preset time period; Geocoding the registered image dataset and DEM data to obtain a detailed lookup table for coordinate conversion between radar coordinate system and geographic coordinate system; A remote sensing ecological index data set of the study area is calculated using a multi-source remote sensing data set; the remote sensing ecological index data set includes remote sensing ecological index data at each second moment in a preset time period; Resampling the remote sensing ecological index data set and the GACOS atmospheric data set based on the DEM data to obtain a resampled index data set and a resampled atmospheric data set; the resampled index data set includes the resampled index data at each second moment in a preset time period, and the resampled atmospheric data set includes the resampled atmospheric data at each first moment in the preset time period; Based on the refined lookup table, the DEM data, the resampled index data set and the resampled atmosphere data set are back-geocoded to obtain the transformed DEM data, the transformed index data set and the transformed atmosphere data set in the radar coordinate system; The registered SAR images in the registered image data set are combined to obtain a plurality of differential interferometry pairs, and differential interferometry is performed on each differential interferometry pair to obtain a plurality of differential interferograms; The undulating phase and the flat ground phase are calculated by using the converted DEM data, and the undulating phase and the flat ground phase are removed from each differential interferogram to obtain a plurality of differential interferograms after preliminary removal; Based on the converted index data set and the plurality of differential interference patterns after preliminary removal, a plurality of high coherence points are determined, and based on all the high coherence points, a high coherence point is selected as a reference point; For each preliminary removed differential interferogram, phase unwrapping is performed on each high coherence point based on a reference point to obtain an unwrapped phase of each high coherence point, and an iterative regression analysis correction is performed on the unwrapped phase of each high coherence point to remove the terrain phase from the preliminary removed differential interferogram to obtain a second removed differential interferogram; For each second-removed differential interferogram, based on the time of the differential interference pair corresponding to the second-removed differential interferogram, two converted atmospheric data at the corresponding time are selected from the converted atmospheric data set, and the atmospheric delay phase is calculated based on the two converted atmospheric data, and the atmospheric delay phase is removed from the second-removed differential interferogram to obtain the removed differential interferogram; Based on all the removed differential interference patterns, the LOS deformation dataset of the study area was calculated.

3. The method for detecting landslide hazard inducing factors at a spatiotemporal scale according to claim 1, characterized in that: The three-dimensional solution model is: D Los =D v cosθ+D N sinθsinβ-D E sinθcosβ; Among them, D Los is the deformation in the LOS direction; D v is the vertical deformation; θ is the radar incident angle; D N is the north-south deformation; β is the radar heading angle; D E It is an east-west deformation; The SPFM model is: Among them, v U is the vertical deformation rate; is the first derivative to the north, H is the terrain height, X N is the northward deformation component; v N is the north-south deformation rate; is the first derivative in the east direction, X E is the eastward deformation component; v E is the east-west deformation rate.

4. The method for detecting landslide hazard inducing factors at a spatiotemporal scale according to claim 1, characterized in that: Based on the three-dimensional deformation rate dataset, the correlation coefficient method is used to detect the first influencing factor from the first type of potential influencing factors at the time scale, and the geographic detector is used to detect the second influencing factor from the second type of potential influencing factors at the spatial scale, including: Obtaining a factor data set of each first-category potential influencing factor and a value of each second-category potential influencing factor in the study area; the factor data set includes the value of the first-category potential influencing factor at each first moment in a preset time period; Based on the three-dimensional deformation rate data set, a vertical deformation rate data set is determined; the vertical deformation rate data set includes vertical deformation rate data at each first moment in a preset time period; Taking the factor data set and vertical deformation rate data set of each first-category potential influencing factor as input, the first influencing factor is detected from the first-category potential influencing factors in time scale using the correlation coefficient method; Taking the value of each second-category potential influencing factor and the value of the vertical deformation rate as input, the geographic detector is used to detect the second influencing factor from the second-category potential influencing factors at a spatial scale; the value of the vertical deformation rate is a random vertical deformation rate data at the first moment in the vertical deformation rate data set or the average value of each vertical deformation rate data at the first moment in the vertical deformation rate data set.

5. The method for detecting landslide hazard inducing factors at a spatiotemporal scale according to claim 4 is characterized in that: Taking the factor dataset and vertical deformation rate dataset of each first-class potential influencing factor as input, the correlation coefficient method is used to detect the first influencing factor from the first-class potential influencing factors on a time scale, specifically including: The vertical deformation rate data set is analyzed by using a principal component analysis method to obtain a principal component data set corresponding to each principal component of the plurality of principal components; the principal component data set includes the value of the principal component at each first moment in a preset time period; For each first-class potential influencing factor, the factor data set corresponding to the first-class potential influencing factor and the principal component data set corresponding to each principal component are used as input to calculate the correlation coefficient between the first-class potential influencing factor and each principal component; The first impact factor is detected from the first type of potential impact factors based on all correlation coefficients.

6. The method for detecting landslide hazard inducing factors at a spatiotemporal scale according to claim 4, characterized in that: Obtain the factor data set of each first-category potential influencing factor and the value of each second-category potential influencing factor in the study area, including: The glacier area in the study area is extracted using multi-source remote sensing datasets, and the glacier area is removed from the study area to obtain the factor dataset of each first-category potential influencing factor and the value of each second-category potential influencing factor in the non-glacier area of ​​the study area; The formula for extracting the glacier area is: GlacierArea= <Ave{Image∈S2:[(D1≥D1 * )∧(D2≤D2 * )∧(NDSI≥NDSI * )]}≥0.5>; Among them, GlacierArea is the glacier area; S2 is the multi-source remote sensing dataset; D1 is the first index; is the threshold of the first index; D2 is the second index; is the threshold of the second index; NDSI is the normalized difference snow index; NDSI * is the threshold of the normalized difference snow cover index; Among them, Red is the reflectivity of the red band; SWIR1 is the reflectivity of the infrared band; Among them, Blue is the reflectivity of the blue band.

7. The method for detecting landslide hazard inducing factors at a spatiotemporal scale according to claim 2, characterized in that: The undulating phase and the flat ground phase are calculated using the converted DEM data, and the undulating phase and the flat ground phase are removed from each differential interferogram to obtain multiple differential interferograms after preliminary removal, including: The ups and downs phase and the flat ground phase are calculated using the converted DEM data; For each differential interferogram, the Goldstein filter is used to perform three phase optimizations on the differential interferogram to obtain the denoised differential interferogram; the windows used for the three phase optimizations are 64*64, 32*32 and 16*16 respectively; The undulating phase and the flat phase are removed from each denoised differential interferogram to obtain a plurality of preliminary removed differential interferograms; Determining a plurality of high coherence points, specifically comprising: determining a plurality of high coherence points using an amplitude deviation method and a coherence coefficient method; Based on the reference point, phase unwrapping is performed on each high coherence point respectively, specifically including: based on the reference point, using the minimum cost flow method to phase unwrapping is performed on each high coherence point respectively; The calculation formula of the atmospheric delay phase is: in, is the atmospheric delay phase; is the difference data between the Mth converted atmospheric data and the Nth converted atmospheric data, and the Mth converted atmospheric data and the Nth converted atmospheric data are two converted atmospheric data; is the atmospheric delay phase of the difference data between the Mth converted atmospheric data and the Nth converted atmospheric data at the reference point (i, j); λ is the wavelength of the SAR sensor.

8. A computer device comprising: A memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement the method for detecting landslide hazard inducing factors at a spatiotemporal scale according to any one of claims 1 to 7.

9. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the method for detecting landslide hazard inducing factors at a spatiotemporal scale according to any one of claims 1 to 7 is implemented.

10. A computer program product, comprising a computer program, characterized in that When the computer program is executed by a processor, the method for detecting landslide hazard inducing factors at a spatiotemporal scale according to any one of claims 1 to 7 is implemented.

Images