A spatiotemporal scale landslide disaster inducing factor detection method and related device

Abstract

The application discloses a landslide disaster inducing factor detection method and related device of space-time scale, and relates to the technical field of disaster inducing factor detection. The method comprises the following steps: calculating a LOS deformation data set of a research area based on a SAR image set, DEM data, a multi-source remote sensing data set and a GACOS atmospheric data set of the research area; taking the LOS deformation data set as input, combining a three-dimensional solving model and an SPFM model to calculate a three-dimensional deformation rate data set of the research area; based on the three-dimensional deformation rate data set, detecting a first influencing factor from a first type of potential influencing factor in a time scale by using a correlation coefficient method; detecting a second influencing factor from a second type of potential influencing factor in a space scale by using a geographic detector; and taking the first influencing factor and the second influencing factor as disaster inducing factors. The application can improve the landslide disaster inducing factor detection precision.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of disaster inducing factor detection, in particular to a landslide disaster inducing factor detection method in space-time scale and related device. BACKGROUND

[0002] The Interferometric Synthetic Aperture Radar (InSAR) technology developed in recent years can implement continuous monitoring on a large-scale research area in all-weather and all-day high-precision, and has shown great advantages and values in many fields, but there is a problem of difficulty in disaster inducing factor (also known as disaster-causing factor) detection in the interpretation and analysis process of the monitoring results in the later stage, especially in the single LOS (Line-of-sight) deformation result. With the InSAR technology being applied in more and more scenarios, the geographical and topographical conditions of the monitored research area are becoming more and more complex, and the single LOS deformation result cannot meet the requirements of disaster inducing factor detection in some scenarios, and it is necessary to inverse the three-dimensional deformation of the ground to obtain more accurate disaster information, therefore, a quantitative analysis considering three-dimensional deformation is needed to complete the disaster inducing factor detection method.

[0003] Domestic and foreign researchers have done a lot of research on disaster inducing factor detection, and have carried out research on different landslide disaster inducing factor detection methods including statistical index (Statistical Index, SI), logistic regression (Logistic Regression, LR) model, field geological observation, random forest model, etc., but the results are not very ideal, and the operation complexity is high. Geographical detector is a new spatial statistics method for detecting spatial heterogeneity, which can show good quantitative effect on factor differentiation in space scale, and can accurately judge the linear and nonlinear interaction between each factor, but it is difficult for geographical detector to quantify the influence degree of potential disaster inducing factor in time scale, and it is necessary to supplement the disaster inducing factor detection method in time scale. SUMMARY

[0004] The purpose of the present application is to provide a landslide disaster inducing factor detection method in space-time scale and related device, which can solve the problems that the single LOS deformation result cannot meet the requirements of disaster inducing factor detection in some scenarios and it is difficult to complete the disaster inducing factor detection in time scale, and improve the landslide disaster inducing factor detection accuracy.

[0005] To achieve the above purpose, the present application provides the following solutions:

[0006] In a first aspect, the application provides a landslide disaster inducing factor detection method in a space-time scale, which comprises the following steps:

[0007] obtaining a SAR image set, DEM data, a multi-source remote sensing data set and a GACOS atmospheric data set of a research area; the SAR image set comprises SAR images at each first time in a preset time period, the multi-source remote sensing data set comprises multi-source remote sensing data at each second time in the preset time period, and the GACOS atmospheric data set comprises GACOS atmospheric data at each first time in the preset time period; the preset time period comprises a satellite ascending track process and a satellite descending track process;

[0008] calculating a LOS direction deformation data set of the research area based on the SAR image set, the DEM data, the multi-source remote sensing data set and the GACOS atmospheric data set; the LOS direction deformation data set comprises LOS direction deformation data at each first time in the preset time period;

[0009] calculating a three-dimensional deformation rate data set of the research area by taking the LOS direction deformation data set as input and combining a three-dimensional solving model and an SPFM model; the three-dimensional deformation rate data set comprises three-dimensional deformation rate data at each first time in the preset time period, the three-dimensional deformation rate data comprises a value of a three-dimensional deformation rate of each position point in the research area, and the three-dimensional deformation rate comprises a vertical deformation rate, a north-south direction deformation rate and an east-west direction deformation rate;

[0010] detecting a first influencing factor from a first type of potential influencing factor in a time scale and a second influencing factor from a second type of potential influencing factor in a space scale based on the three-dimensional deformation rate data set by using a correlation coefficient method and a geographic detector, and taking the first influencing factor and the second influencing factor as disaster inducing factors; the first type of potential influencing factor comprises precipitation and temperature, and the second type of potential influencing factor comprises at least two of DEM data, slope, slope direction, curvature, NDVI, land cover type, land surface temperature, water system and fault.

[0011] In a second aspect, the application provides a computer device, which comprises a memory, a processor and a computer program stored in the memory and executable on the processor, and the processor executes the computer program to implement the landslide disaster inducing factor detection method in a space-time scale.

[0012] In a third aspect, the application provides a computer readable storage medium, which stores a computer program, and the computer program is executed by a processor to implement the landslide disaster inducing factor detection method in a space-time scale.

[0013] In a fourth aspect, the present application provides a computer program product comprising a computer program which, when executed by a processor, implements the above-mentioned landslide disaster inducing factor detection method in the space-time scale.

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

[0015] The present application provides a landslide disaster inducing factor detection method in the space-time scale and related devices, based on the SAR image set, DEM data, multi-source remote sensing data set and GACOS atmospheric data set of the research area, the LOS deformation data set of the research area is calculated, the three-dimensional deformation rate data set of the research area is calculated by taking the LOS deformation data set as input, combining the three-dimensional solving model and the SPFM model, and the three-dimensional deformation rate data set is used to complete the landslide disaster inducing factor detection, by introducing the three-dimensional deformation rate, the problem that the single LOS deformation result cannot meet the disaster inducing factor detection requirement in some scenarios can be solved. And when the three-dimensional deformation rate data set is used to complete the landslide disaster inducing factor detection, based on the three-dimensional deformation rate data set, the first influencing factor is detected from the first type of potential influencing factor in the time scale by using the correlation coefficient method, the second influencing factor is detected from the second type of potential influencing factor in the space scale by using the geographic detector, and the first influencing factor and the second influencing factor are taken as the disaster inducing factor, by introducing the time scale, the landslide disaster inducing factor detection can be completed by fusing the space-time scale, and the problem that it is difficult to complete the disaster inducing factor detection in the time scale is solved. The present application can improve the landslide disaster inducing factor detection precision. BRIEF DESCRIPTION OF DRAWINGS

[0016] In order to more clearly illustrate the technical solutions in the embodiments of the present application or the prior art, the drawings needed in the embodiments will be briefly introduced below. Obviously, the drawings in the following description are only some embodiments of the present application, and other drawings can also be obtained by those skilled in the art without creative labor.

[0017] Figure 1 An application environment diagram of a landslide disaster inducing factor detection method in the space-time scale provided for Embodiment 1 of the present application.

[0018] Figure 2 A flowchart of a landslide disaster inducing factor detection method in the space-time scale provided for Embodiment 1 of the present application.

[0019] Figure 3 A detailed flowchart of a landslide disaster inducing factor detection method in the space-time scale provided for Embodiment 1 of the present application.

[0020] Figure 4SAR image coverage diagram of a research area in the south of Dingjie County, Shigatse City, Tibet Autonomous Region, provided for Embodiment 1 of the present application.

[0021] Figure 5 Three-dimensional deformation rate decomposition diagram provided for Embodiment 1 of the present application; wherein, Figure 5 (a) in the above is an east-west deformation rate, Figure 5 (b) in the above is a north-south deformation rate, Figure 5 (c) in the above is a vertical deformation rate.

[0022] Figure 6 Second type of potential influencing factor diagram provided for Embodiment 1 of the present application; wherein, Figure 6 (a) in the above is DEM (Digital Elevation Model, digital elevation model) data, Figure 6 (b) in the above is slope, Figure 6 (c) in the above is aspect, Figure 6 (d) in the above is curvature, Figure 6 (e) in the above is NDVI (Normalized Vegetation Index, normalized vegetation index), Figure 6 (f) in the above is land cover type, Figure 6 (g) in the above is land surface temperature (LST), Figure 6 (h) in the above is water system, Figure 6 (i) in the above is fault.

[0023] Figure 7 Disaster-inducing factor detection result diagram under a spatial scale, provided for Embodiment 1 of the present application.

[0024] Figure 8 Structural diagram of a computer device provided for Embodiment 2 of the present application. DETAILED DESCRIPTION

[0025] The technical solutions in the embodiments of the present application will be described clearly and completely below with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are only some of the embodiments of the present application, rather than all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those of ordinary skill in the art without creative work fall within the scope of protection of the present application.

[0026] Embodiment 1

[0027] The spatiotemporal scale landslide disaster-inducing factor detection method provided by the embodiments of the present application can be applied to, for example, Figure 1The application environment shown. Among them, the terminal communicates with the server through the network. The data storage system can store the data required by the server to process. The data storage system can be separately arranged, or integrated on the server, or placed on the cloud or other servers. The terminal can send a to-be-processed detection request to the server. After receiving the to-be-processed detection request, the server obtains a SAR image set of a research area, DEM data, a multi-source remote sensing data set and a GACOS atmospheric data set for the to-be-processed detection request; based on the SAR image set, the DEM data, the multi-source remote sensing data set and the GACOS atmospheric data set, the LOS deformation data set of the research area is calculated; taking the LOS deformation data set as input, combining the three-dimensional solving model and the SPFM model, the three-dimensional deformation rate data set of the research area is calculated; based on the three-dimensional deformation rate data set, the first influence factor is detected from the first type of potential influence factor in the time scale by using the correlation coefficient method, and the second influence factor is detected from the second type of potential influence factor in the spatial scale by using the geographic detector, and the first influence factor and the second influence factor are taken as the disaster inducing factor. The server can feed back the detection result of the disaster inducing factor for the detection request to the terminal.

[0028] In addition, in some embodiments, the spatiotemporal scale landslide disaster inducing factor detection method can also be implemented by the server or the terminal alone, such as being directly processed by the terminal for the to-be-processed detection request, or being acquired from the data storage system by the server and processed for the to-be-processed detection request.

[0029] Among them, the terminal can be, but is not limited to, various desktop computers, notebook computers, smart phones, tablet computers, Internet of Things devices and portable wearable devices, the Internet of Things devices can be smart speakers, smart televisions, smart air conditioners, smart vehicle-mounted devices, etc., and the portable wearable devices can be smart watches, smart bracelets, head-mounted devices, etc. The server can be implemented by an independent server or a server cluster composed of multiple servers, and can also be a cloud server.

[0030] In an exemplary embodiment, as shown in Figure 2 and Figure 3 A spatiotemporal scale landslide disaster inducing factor detection method is provided, which is executed by a computer device, specifically, can be executed by a terminal or a server alone, or can be executed by a terminal and a server together. In the embodiments of the present application, the method is applied to the server in Figure 1 for example, including the following steps.

[0031] Step S1, obtaining a SAR image set of a research area, DEM data, a multi-source remote sensing data set and a GACOS atmospheric data set; the SAR image set includes SAR images of each first time in a preset time period, the multi-source remote sensing data set includes multi-source remote sensing data of each second time in the preset time period, and the GACOS atmospheric data set includes GACOS atmospheric data of each first time in the preset time period; the preset time period includes a satellite ascending track process and a satellite descending track process.

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

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

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

[0035] By implementing steps S1 to S4 above, this embodiment, after determining the LOS-oriented deformation dataset of the study area, uses the LOS-oriented deformation dataset as input and combines the three-dimensional solution model and SPFM model to calculate the three-dimensional deformation rate dataset of the study area, thereby inverting the three-dimensional deformation of the land surface to obtain more accurate disaster information. Subsequently, based on the quantitative analysis of three-dimensional deformation, disaster-inducing factor detection is completed, solving the problem that the single LOS-oriented deformation result cannot meet the requirements of disaster-inducing factor detection in some scenarios. Furthermore, when completing disaster-inducing factor detection 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. The first and second influencing factors are used as disaster-inducing factors, thereby detecting disaster-inducing factors on a spatiotemporal scale and solving the problem of difficulty in quantifying the degree of influence of disaster-inducing factors on a time scale.

[0036] In this embodiment, SAR image sets, DEM data, multi-source remote sensing datasets, and GACOS atmospheric datasets for the study area are acquired. The SAR image set can be an L1-level time-series product acquired by a spaceborne synthetic aperture radar covering the study area, including SAR images at each first moment within a preset time period. The SAR images cover the study area, such as the SAR images acquired by Sentinel-1 covering the southern study area of ​​Dingjie County, Shigatse City, Tibet Autonomous Region. Figure 4 As shown, the multi-source remote sensing dataset can be a multispectral time-series dataset of the study area acquired using sensors such as Sentinel-2 MSI, MODIS, and Landsat MSS satellites. It includes multi-source remote sensing data for each second time point within a preset time period. The multi-source remote sensing data covers the study area. The first and second time points can be the same or different. Based on the acquisition time of each SAR image, GACOS atmospheric data covering the study area is downloaded. The GACOS atmospheric dataset includes GACOS atmospheric data for each first time point within the preset time period. The GACOS atmospheric data covers the study area, and the pixel value of each pixel in the GACOS atmospheric data is the atmospheric delay distance. The preset time period includes the satellite ascent and descent processes; therefore, the SAR image set includes multiple SAR images from the satellite ascent process and multiple SAR images from the satellite descent process.

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

[0038] The LOS deformation data set of the research region is calculated based on the SAR image set, DEM data, multi-source remote sensing data set and GACOS atmospheric data set, and specifically includes:

[0039] (1) The SAR images in the SAR image set are registered based on the DEM data to obtain a registered image data set of the research region. The registered image data set includes a registered SAR image at each first time in a preset time period, and the registered SAR image covers the research region.

[0040] In the actual application, the SAR image set obtained by Sentinel 1 is used, and the image registration process is assisted by the external DEM data to obtain the registered single-view complex image. The "three-step registration method" is mainly used to register the SAR image set in the TOPS mode, and specifically includes: in the first step, the geometric registration method is used for coarse registration based on the DEM data; in the second step, the intensity cross-correlation method is used for further registration; and in the third step, the enhanced spectral diversity method is used for fine registration. Through the three-step registration method, a registration accuracy of one thousandth can be achieved to obtain the registered image data set.

[0041] (2) The registered image data set and the DEM data are geocoded to obtain a fine lookup table for coordinate conversion between the radar coordinate system and the geographic coordinate system.

[0042] In the actual application, DEM data with a resolution basically consistent with the SAR image is introduced, and the DEM data and satellite orbit data are combined to generate a simulated SAR image at each first time in a preset time period. Then, the simulated SAR image and the real SAR image (i.e., the registered SAR image in the registered image data set) are registered to establish a registration relationship therebetween, realize the conversion of the radar image from slant range projection to orthographic projection, and obtain the fine lookup table. The "two-step registration method" is mainly used to obtain the fine lookup table, and specifically includes: in the first step, the DEM data and the satellite orbit data for SAR imaging are used to generate a simulated SAR image in combination with a SAR imaging model, and the simulated SAR image and the real SAR image are coarsely registered to estimate an 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 respectively the distance and azimuth pixel number in the radar coordinate system, g is the corresponding relationship between the radar coordinate system and the geographic coordinate system,

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

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

[0045] This embodiment is based on a multi-source remote sensing dataset and uses the GEE (Google Earth Engine) platform to retrieve remote sensing ecological indexes, including but not limited to NDVI, LST, NDSI (Normalized differential snow index), and other remote sensing ecological indexes of interest. The remote sensing ecological index dataset is obtained.

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

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

[0048] where GlacierArea is the glacier area; Ave represents the average value; S2 is a multi-source remote sensing dataset obtained by Sentinel-2 covering the study area in a preset time period (such as from January 1, 2021 to September 18, 2023); D1 is the first index; * indicates the threshold value calculated using the OSTU (maximum inter-class variance method), is the threshold value of the first index calculated using the OSTU; ∧ indicates and; D2 is the second index; ​is a threshold value of the second index calculated by using OSTU; NDSI is a normalized difference snow index; NDSI * is a threshold value of the normalized difference snow index calculated by using OSTU. When the first index is greater than or equal to the threshold value of the first index, the second index is greater than or equal to the threshold value of the second index, and the normalized difference snow index is greater than or equal to the threshold value of the normalized difference snow index, the pixel value of the pixel point is 1, which represents that the pixel point belongs to the glacier region.

[0049]

[0050] wherein Red is the reflectivity of the red band; SWIR1 is the reflectivity of the infrared band.

[0051]

[0052] wherein Blue is the reflectivity of the blue band.

[0053]

[0054] wherein Green is the reflectivity of the green band.

[0055] (4) Resample 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 including resampled index data of each second time in the preset time period, and the resampled atmospheric data set including resampled atmospheric data of each first time in the preset time period, the resampled index data and the resampled atmospheric data both covering the research region.

[0056] For the remote sensing ecological index data set and the GACOS atmospheric data set of the research region, resample them according to the geographic coded DEM data to obtain the range, corner point latitude and longitude, resolution and other parameters of the geographic coded DEM data, resample the remote sensing ecological index data and the GACOS atmospheric data to the above geographic parameter configuration to obtain multi-source remote sensing products consistent with the geographic parameter configuration of the geographic coded DEM data, specifically, first calculate the latitude and longitude of the center point of each grid in the DEM data according to the range, corner point latitude and longitude, resolution and other parameters of the geographic coded DEM data, the calculation formula being as follows:

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

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

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

[0060] Then, using the nearest neighbor, bilinear, B-spline, etc. resampling method, the remote sensing ecological index dataset and the GACOS atmospheric dataset are resampled according to the resampling formula to unify the remote sensing ecological index dataset and the GACOS atmospheric dataset to the geographic coordinate system and the parameter configuration of the DEM data after geographic coding. 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; Z(x, y) is the pixel value of the pixel point (x, y) of the original data.

[0063] (5) Based on the fine lookup table, the DEM data, the resampled index dataset and the resampled atmospheric dataset are respectively backward geocoded to map the DEM data, the resampled index dataset and the resampled atmospheric dataset in the geographic coordinate system to the radar coordinate system to obtain the converted DEM data, the converted index dataset and the converted atmospheric dataset in the radar coordinate system.

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

[0065]

[0066] wherein, is the pixel coordinate of a certain ground object in the radar coordinate system B; Interp is the convolution interpolation method: f = (Σ i w i ·f i ) / Σ i w i , w i is the weight of the i-th adjacent point, f iis the pixel value of the ith neighboring point, the neighboring point is a pixel point in the neighborhood of the pixel point where a certain ground object is located, the neighborhood can be determined according to requirements, such as an 8-neighborhood, and the weight can be determined by using a distance weighting method, that is, w i = 1 / d i , d i is the distance from the neighboring point to the pixel point where a certain ground object is located; L x (), L y () are the mapping relationships of the geographic coordinate system A and the radar coordinate system B in the rows and columns, respectively, and (xA, yA) is the pixel coordinate of the same ground object in the geographic coordinate system A.

[0067] (6) Combining the registered SAR images in the registered image data set to obtain a plurality of differential interference pairs, and performing differential interference on each differential interference pair to obtain a plurality of differential interference maps.

[0068] The embodiment combines differential interference pairs based on the small baseline set idea. The small baseline set idea is an idea of reducing the spatial decorrelation effect caused by a large baseline by selecting a group of interference maps with small time and spatial baselines. Specifically, in the combination mode of multiple master images, the differential interference pairs are combined under the restrictions that the time baseline is less than 90 days and the spatial baseline is less than 400 m, that is, each registered SAR image in the registered image data set is taken as a master image, and the registered SAR images in the registered image data set with a time baseline less than 90 days and a spatial baseline less than 400 m from the master image are selected as selected images, and the master image and each selected image are combined to form a differential interference pair.

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

[0070] (7) The undulating phase and flat ground phase are calculated by using the converted DEM data, and the undulating phase and flat ground phase are removed in each differential interference map to obtain a plurality of preliminary removed differential interference maps.

[0071] According to the converted DEM data in the radar coordinate system generated by geographic coding and satellite orbit parameters, the undulating phase and flat ground phase of the terrain are simulated, and the above two terrain-related undulating phases and flat ground phases are subtracted in each differential interference map to preliminarily remove the terrain-related phases, thereby obtaining a plurality of preliminary removed differential interference maps. The number of preliminary removed differential interference maps is the same as that of differential interference maps.

[0072] The embodiment also introduces Goldstein phase optimization. The generated differential interference map is iteratively phase-optimized by using a Goldstein filter to sequentially set the window sizes of 64*64, 32*32 and 16*16, so as to weaken the speckle noise in the differential interference map.

[0073] At this time, the relief phase and the flat phase are calculated by using the converted DEM data, and the relief phase and the flat phase are removed in each differential interferogram to obtain a plurality of preliminary removed differential interferograms, specifically including: calculating the relief phase and the flat phase by using the converted DEM data; for each differential interferogram, performing three times of phase optimization on the differential interferogram by using a Goldstein filter to obtain a denoised differential interferogram, wherein the windows used for the three times of phase optimization are 64*64, 32*32 and 16*16 respectively; removing the relief phase and the flat phase in each denoised differential interferogram to obtain a plurality of preliminary removed differential interferograms.

[0074] (8) Based on the converted exponential data set and the plurality of preliminary removed differential interferograms, a plurality of high-coherence points are determined, and based on all the high-coherence points, one high-coherence point is selected as a reference point.

[0075] Among them, the plurality of high-coherence points are determined specifically by using the amplitude deviation method and the coherence coefficient method.

[0076] Generally, the region of interest will be obtained based on the remote sensing ecological index of interest through the threshold value. In this embodiment, the glacier region extracted according to the NDSI index is further combined with the overlap, shadow mask and glacier region to generate a mask file. The mask file includes pixel points belonging to the overlap, shadow and glacier region. 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 all the preliminary removed differential interferograms are selected by using the amplitude deviation method and the coherence coefficient method. The high-coherence points can be the intersection of the high-coherence points selected by satisfying the preset amplitude deviation threshold value and the high-coherence points selected by satisfying the preset coherence coefficient threshold value, or the union of the high-coherence points selected by satisfying the preset amplitude deviation threshold value and the high-coherence points selected by satisfying the preset coherence coefficient threshold value. The high-coherence points selected by satisfying the preset amplitude deviation threshold value refer to the amplitude deviation of the pixel points being greater than the preset amplitude deviation threshold value, and the high-coherence points selected by satisfying the preset coherence coefficient threshold value refer to the coherence coefficient of the pixel points being greater than the preset coherence coefficient threshold value. The amplitude deviation and the coherence coefficient are used to obtain the high-coherence points, and the preset amplitude deviation threshold value and the preset coherence coefficient threshold value are 0.15 and 0.8 respectively.

[0077] The calculation formula of the amplitude deviation D A is:

[0078]

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

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

[0081]

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

[0083] (9) For each of the preliminary-removed differential interferograms, phase unwrapping is performed on each of the high-coherence points based on the reference point to obtain the unwrapped phase of each of the high-coherence points, and iterative regression analysis correction is performed on the unwrapped phase of each of the high-coherence points to remove the terrain phase in the preliminary-removed differential interferogram, thereby obtaining a re-removed differential interferogram.

[0084] wherein, phase unwrapping is performed on each of the high-coherence points based on the reference point, and the phase unwrapping specifically comprises: performing phase unwrapping on each of the high-coherence points based on the reference point by using the minimum cost flow method.

[0085] In this embodiment, the reference point is selected, and specifically, the reference point is selected from all the high-coherence points according to the center principle, the stability principle, the high-coherence principle and the non-island principle. For each of the high-coherence points, phase unwrapping is performed on the high-coherence point based on the reference point by using the minimum cost flow method in the spatial dimension, the unwrapped phase of the high-coherence point is obtained by solving, and iterative regression analysis correction is performed on the unwrapped phase of the high-coherence point until the height correction value of the iterative result is less than 0.1 m, thereby removing the terrain phase and obtaining a re-removed differential interferogram. The number of the re-removed differential interferograms is the same as that of the preliminary-removed differential interferograms.

[0086] In this embodiment, the iterative regression analysis correction is specifically performed by using an interference phase model, and the interference phase model is as follows:

[0087]

[0088] wherein, δφ i,j is the unwrapped phase of the high-coherence point; λ is the wavelength; v j is the deformation rate; Δt j =t2-t1, t2>t1, t2 is the time of the second scene-registered SAR image, t1 is the time of the first scene-registered SAR image, the first scene-registered SAR image and the second scene-registered SAR image forming a differential interference pair; δ i,j is the deformation residual; B j is the vertical baseline; Δz i is the DEM height error; and ρ​i is the slant range; θ i is the incidence angle; is the atmospheric delay phase; t 2,j is the time corresponding to the high-coherence point in the second scene-registered SAR image; t 1,j is the time corresponding to the high-coherence point in the first scene-registered SAR image; n i is the thermal phase noise.

[0089] (10) For each of the re-removed differential interferograms, based on the time of the differential interference pair corresponding to the re-removed differential interferogram, two pieces of 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 pieces of converted atmospheric data, and the atmospheric delay phase is removed from the re-removed differential interferogram to obtain a removed differential interferogram.

[0090] In this embodiment, the GACOS atmospheric data is differentially calculated according to the differential interference pair combination, and the relative value of each piece of GACOS atmospheric data is calculated based on the reference point. Specifically, the differential GACOS atmospheric data is subtracted by the atmospheric delay phase of the reference point, that is, The purpose of this is to unify the reference of the atmospheric delay phase of each piece of GACOS atmospheric data, and convert the GACOS atmospheric data into a two-way phase, that is, S M · 4π / λ, the obtained atmospheric delay phase is removed in the optimized unwrapped phase, so as to remove the atmospheric delay phase in the optimized unwrapped phase, and obtain a removed differential interferogram.

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

[0092]

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

[0094] (11) Based on all the removed differential interferograms, a LOS deformation data set of the research area is calculated.

[0095] In this embodiment, the LOS deformation data set is taken as input, and a three-dimensional deformation rate data set of the research region is calculated by combining a three-dimensional solving model and an SPFM model. The three-dimensional deformation rate data set includes three-dimensional deformation rate data of each first time 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 research region. The three-dimensional deformation rate includes a vertical deformation rate, a north-south deformation rate and an east-west deformation rate.

[0096] Since the SAR image set includes an ascending track data set and a descending track data set, the ascending track data set and the descending track data set are respectively processed as described above to obtain the LOS time sequence deformation of the ascending track process and the descending track process.

[0097] In this embodiment, the ascending and descending track LOS deformations are combined, the SPFM model is introduced, and the three-dimensional deformation rate results of the research region are solved by combining the three-dimensional solving model and the SPFM model to obtain the vertical, north-south and east-west deformation rates. The three-dimensional solving model is:

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

[0099] wherein, D Los is the LOS deformation; D v is the vertical deformation; θ is the radar incidence angle; D N is the north-south deformation; β is the radar heading angle; D E is the east-west deformation.

[0100] The SPFM model is:

[0101]

[0102] wherein, v U is the vertical deformation rate; is the first derivative in the north direction, H is the terrain height, X N is the north deformation component; v N is the north-south deformation rate; is the first derivative in the east direction, X E is the east deformation component; v E is the east-west deformation rate.

[0103] Taking the research region in the south 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 the embodiment, based on the three-dimensional deformation rate dataset, the first influencing factor is detected from the first type of potential influencing factors in the time scale by using the correlation coefficient method, the second influencing factor is detected from the second type of potential influencing factors in the space scale by using the geographic detector, and the first influencing factor and the second influencing factor are taken as the disaster inducing factors. The first type of potential influencing factors includes precipitation and temperature, and the second type of potential influencing factors includes at least two of DEM data, slope, slope direction, curvature, NDVI, land cover type, land surface temperature, water system and fault. Taking the southern research area of Dingjie County, Shigatse City, Tibet Autonomous Region as an example, the second type of potential influencing factors are shown in FIG. 1. Figure 6 .

[0105] In the embodiment, the geographic detector and the PCA analysis (Principal Component Analysis) are fused to detect the disaster inducing factors in the time and space scales, and specifically include: the geographic detector is used to detect the explanatory power of the second type of potential influencing factors on the deformation result, 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 the explanatory power; the PCA is used to extract the principal component of the vertical deformation rate of the research area, and the main component significant influence level of the first type of potential influencing factors such as temperature and precipitation on the landslide disaster in the research area is detected.

[0106] In the embodiment, based on the three-dimensional deformation rate dataset, the first influencing factor is detected from the first type of potential influencing factors in the time scale by using the correlation coefficient method, the second influencing factor is detected from the second type of potential influencing factors in the space scale by using the geographic detector, and specifically includes:

[0107] (1) Obtain the factor dataset of each first type of potential influencing factor and the value of each second type of potential influencing factor of the research area. The factor dataset includes the value of the first type of potential influencing factor at each first time in a preset time period. The value of the first type of potential influencing factor is the value of the first type of potential influencing factor at each position point in the research area. The value of the second type of potential influencing factor is the value of the second type of potential influencing factor at each position point in the research area.

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

[0109] The factor data set of each first type of potential influencing factor and the value of each second type of potential influencing factor of the research area are obtained, specifically including: extracting a glacier region of the research area by using a multi-source remote sensing data set, and removing the glacier region in the research area to obtain the factor data set of each first type of potential influencing factor of a non-glacier region of the research area and the value of each second type of potential influencing factor.

[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 including vertical deformation rate data of each first time in a preset time period, and the vertical deformation rate data including a value of a vertical deformation rate of each position point in the research area (after removing the glacier region).

[0111] (3) Taking the factor data set of each first type of potential influencing factor and the vertical deformation rate data set as input, a first influencing factor is detected from the first type of potential influencing factor in the time scale by using a correlation coefficient method.

[0112] (4) Taking the value of each second type of potential influencing factor and the value of the vertical deformation rate as input, a second influencing factor is detected from the second type of potential influencing factor in the spatial scale by using a geographic detector, and the value of the vertical deformation rate is vertical deformation rate data of a random first time in the vertical deformation rate data set or an average value of vertical deformation rate data of each first time in the vertical deformation rate data set.

[0113] At this time, taking the factor data set of each first type of potential influencing factor and the vertical deformation rate data set as input, a first influencing factor is detected from the first type of potential influencing factor in the time scale by using a 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 principal component in a plurality of principal components, the principal component data set including a value of the principal component of each first time in a preset time period, and the value of the principal component being a value of the principal component of each position point in the research area (after removing the glacier region).

[0115] In this embodiment, PCA is used to extract principal components of the vertical deformation rate data set in the research area, and a plurality of principal components with a cumulative probability of 98% of the principal components are obtained to participate in disaster inducing factor detection calculation in the time scale.

[0116] (2) For each first type of potential influencing factor, taking the factor data set corresponding to the first type of potential influencing factor and the principal component data set corresponding to each principal component as input, the correlation coefficient between the first type of potential influencing factor and each principal component is calculated.

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

[0118]

[0119] Wherein, r is the correlation coefficient; n is the number of samples, that is, the number of first time in the preset time period; X i is the ith sample X; is the mean of sample X; Y i is the ith sample Y; is the mean of sample Y; sample X and sample Y represent each principal component and the first type of potential influencing factor respectively.

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

[0121] Sort all correlation coefficients in descending order of correlation, and extract the first type of potential influencing factor located in the top N as the first influencing factor.

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

[0123] Table 1 Time scale disaster inducing factor detection principal component significance test table

[0124]

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

[0126] Using geographic detector to complete the detection process belongs to the existing mature technology, hereinafter only part of the steps are introduced: the second type of potential influencing factor is reclassified according to the parameter optimization principle, and the absolute value deviation of the vertical deformation rate of the spatial dimension is calculated, the second type of potential influencing factor is reclassified, and the quantile discretization method is used according to the parameter optimization 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, the absolute value deviation of the vertical deformation rate of the spatial dimension is calculated: x i is the vertical deformation rate of the ith pixel point, and it is assumed that the result obeys the standard normal distribution,

[0127] The explanatory power of the second type of potential influencing factor on the deformation result is detected by using the geographic detector, and a relative q value (q / p) is calculated, the q value representing the explanatory power of one variable on another variable, and the p value representing a confidence level of the explanatory power, when it is ensured that the confidence level represented by the p value is below 0.3, the corresponding q value and p value are involved in the calculation, and the calculation formula of the q value is:

[0128]

[0129] SST = Nσ 2 ;

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

[0131] The calculation result of the second influencing factor detected by using the geographic detector is shown in Table 1. Figure 7 .

[0132] The embodiment discloses a kind of multi-source remote sensing data fusion's time and space scale alpine region landslide disaster induced factor detection method, it is a kind of disaster induced factor detection method being carried out in time and space scale with considering three-dimensional deformation, improve the detection precision of disaster induced factor, this method includes: obtaining the L1 level time sequence product of satellite-borne synthetic aperture radar covering research area, and combining external DEM data auxiliary image registration, obtain the single view complex map after registration;DEM is encoded to SAR image, obtain the fine geographic-SAR lookup table of geographic coordinate system and radar coordinate system, ensure that geographic coding accuracy is better than 0.5Pixel;Using GEE platform, remote sensing ecological index is retrieved based on multi-source remote sensing data, including but not limited to NDVI, LST, NDSI etc., simultaneously obtain the GACOS atmospheric data of research area, according to the DEM data after geographic coding, remote sensing ecological index and GACOS atmospheric data are resampled, obtain the multi-source remote sensing product consistent with the geographic parameter configuration of DEM data after geographic coding;According to the fine geographic-SAR lookup table, the remote sensing ecological index of interest and GACOS atmospheric data are backwardly geographic coded, and the remote sensing ecological index of interest and GACOS atmospheric data under radar coordinate system are obtained;Based on the idea of small baseline set, difference interference pair is combined, and the single view complex map after registration is differentially interfered, to remove topographic-related phase, iteration is carried out Goldstein phase optimization;Combined with overlap, shadow mask and the interested remote sensing ecological index under suitable threshold, mask file is generated, and based on mask file, the high-coherence point of interest in difference interference map is selected by amplitude dispersion method and coherence coefficient method;Select reference point, and use least cost flow method to solve high-coherence point unwrapping phase, iteratively carry out regression analysis correction and remove terrain phase;Based on reference point, the relative value of each GACOS atmospheric data is calculated and converted into two-way phase, and the obtained atmospheric delay phase is removed in the optimized unwrapping phase;Joint lift-off track LOS deformation, introduce SPFM model to solve the three-dimensional deformation rate results of research area;According to the time and space dimension three-dimensional deformation rate results obtained, the disaster induced factor is detected from time and space scale by fusing spatial dimension geographic detector and time dimension PCA analysis.This multi-source remote sensing data fusion's time and space scale alpine region landslide disaster induced factor detection method provided in the embodiment solves the situation that disaster induced factor is difficult to detect due to LOS deformation, supplements the disaster induced factor detection method under time scale, and improves the level of landslide disaster monitoring and analysis.

[0133] The method provided by the embodiment provides a multi-source remote sensing data fusion time-space scale alpine landslide disaster inducing factor detection method. First, a potential disaster factor database is established by using GEE, and a region of interest is selected by fusing SAR and multispectral remote sensing images. Meanwhile, considering that single LOS deformation results cannot reflect the real ground deformation, a SPFM model is introduced to invert the ground three-dimensional deformation, obtain the time-space dimensional deformation data of the research region, and then fuse the geographic detector and the PCA analysis method to detect the disaster inducing factors in the time-space dimension, so as to obtain a comprehensive disaster inducing factor quantitative index system in the time-space dimension, effectively quantify the geological disaster inducing factor index, and supplement the deficiencies of related work.

[0134] The application also provides an application scenario of the time-space scale landslide disaster inducing factor detection method. Specifically, the time-space scale landslide disaster inducing factor detection method 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 for fusing the geographic detector and the principal component analysis technology to detect the disaster inducing factors in the time-space scale, and establish a landslide disaster influencing factor system of the research region. The interpretation link is used for interpreting the landslide based on the landslide disaster influencing factor system. The time-space scale landslide disaster inducing factor detection method belongs to the detection link.

[0135] Embodiment 2

[0136] In an exemplary embodiment, a computer device, which can be a server or a terminal, is provided, and an internal structure diagram of the computer device can be as shown in Figure 8 The computer device includes a processor, a memory, an input / output interface (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 operating system and the computer program in the non-volatile storage medium to run. 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 external devices. The communication interface of the computer device is used to communicate with external terminals through network connection. The computer program is executed by the processor to implement a time-space scale landslide disaster inducing factor detection method.

[0137] Those skilled in the art can understand that Figure 8It is to be noted that the structure shown in the figure is only a block diagram of part of the structure related to the scheme of the present application, and does not constitute a limitation on the computer device to which the scheme of the present application is applied. The specific computer device can 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, comprising a memory and a processor, the memory storing a computer program, and the processor implementing the method for detecting landslide disaster inducing factors in space-time scale in embodiment 1 when executing the computer program.

[0139] Embodiment 3

[0140] In an exemplary embodiment, a computer readable storage medium is provided, storing a computer program, and the computer program implementing the method for detecting landslide disaster inducing factors in space-time scale in embodiment 1 when executed by a processor.

[0141] Embodiment 4

[0142] In an exemplary embodiment, a computer program product is provided, comprising a computer program, and the computer program implementing the method for detecting landslide disaster inducing factors in space-time scale in embodiment 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 for analysis, stored data, displayed data, etc.) involved in the present application are all information and data authorized by the user or authorized by all parties, and the collection, use and processing of related data need to comply with relevant regulations.

[0144] The technical features of the above embodiments can be combined arbitrarily. To make the description concise, all possible combinations of the technical features in the above embodiments are not described, however, as long as the combination of the technical features does not exist contradictory, it should be considered as the scope of the present application.

[0145] The principles and implementation modes of the present application are described by using specific examples in the present application, and the above description of the embodiments is only used to help understand the method and its core idea of the present application; at the same time, for those skilled in the art, according to the idea of the present application, the specific implementation mode and application range will be changed. In conclusion, the content of the present application should not be understood as a limitation.

Claims

1. A method for detecting landslide disaster inducing factors in a spatiotemporal scale, characterized in that, The landslide disaster inducing factor detection method of the space-time scale comprises: SAR image sets, DEM data, multi-source remote sensing data sets and GACOS atmospheric data sets of a research region are acquired; the SAR image sets comprise SAR images of each first time in a preset time period, the multi-source remote sensing data sets comprise multi-source remote sensing data of each second time in the preset time period, and the GACOS atmospheric data sets comprise GACOS atmospheric data of each first time in the preset time period; the preset time period comprises a satellite ascending track process and a satellite descending track process; Based on the SAR image sets, the DEM data, the multi-source remote sensing data sets and the GACOS atmospheric data sets, LOS deformation data sets of the research region are calculated; the LOS deformation data sets comprise LOS deformation data of each first time in the preset time period; Based on the LOS deformation data sets, three-dimensional deformation rate data sets of the research region are calculated by combining a three-dimensional solving model and an SPFM model; the three-dimensional deformation rate data sets comprise three-dimensional deformation rate data of each first time in the preset time period, the three-dimensional deformation rate data comprise values of three-dimensional deformation rates of each position point in the research region, and the three-dimensional deformation rates comprise vertical deformation rates, north-south deformation rates and east-west deformation rates; Based on the three-dimensional deformation rate data sets, a first influencing factor is detected from first potential influencing factors in a time scale by using a correlation coefficient method, a second influencing factor is detected from second potential influencing factors in a space scale by using a geographic detector, and the first influencing factor and the second influencing factor are taken as disaster inducing factors; the first potential influencing factors comprise precipitation and temperature, and the second potential influencing factors comprise at least two of DEM data, slope, slope direction, curvature, NDVI, land cover type, land surface temperature, water system and fault. 2.The method according to claim 1, wherein, Based on the SAR image sets, the DEM data, the multi-source remote sensing data sets and the GACOS atmospheric data sets, the LOS deformation data sets of the research region are calculated, and specifically comprise: The SAR images in the SAR image sets are registered based on the DEM data to obtain registered image data sets of the research region; the registered image data sets comprise registered SAR images of each first time in the preset time period; The registered image data sets and the DEM data are geocoded to obtain a fine lookup table for coordinate conversion between a radar coordinate system and a geographic coordinate system; Remote sensing ecological index data sets of the research region are calculated by using the multi-source remote sensing data sets; the remote sensing ecological index data sets comprise remote sensing ecological index data of each second time in the preset time period; The remote sensing ecological index data sets and the GACOS atmospheric data sets are resampled based on the DEM data to obtain resampled index data sets and resampled atmospheric data sets; the resampled index data sets comprise resampled index data of each second time in the preset time period, and the resampled atmospheric data sets comprise resampled atmospheric data of each first time in the preset time period. The DEM data, the resampled exponential data set and the resampled atmospheric data set are respectively backward geocoded based on a fine lookup table to obtain converted DEM data, converted exponential data set and converted atmospheric data set in a radar coordinate system; The registered SAR image in the registered image data set is combined to obtain a plurality of differential interference pairs, and each differential interference pair is subjected to differential interference to obtain a plurality of differential interferograms; The undulating phase and the flat phase are calculated by using the converted DEM data, and the undulating phase and the flat phase are removed in each differential interferogram to obtain a plurality of preliminary removed differential interferograms; Based on the converted exponential data set and the plurality of preliminary removed differential interferograms, a plurality of high-coherence points are determined, and one high-coherence point is selected as a reference point based on all the high-coherence points; For each preliminary removed differential interferogram, the 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 subjected to iterative regression analysis correction to remove the terrain phase in the preliminary removed differential interferogram, thereby obtaining a second removed differential interferogram; For each second removed differential interferogram, two converted atmospheric data corresponding to the corresponding time of the differential interference pair of the second removed differential interferogram are selected from the converted atmospheric data set, and the atmospheric delay phase is calculated based on the two converted atmospheric data and removed in the second removed differential interferogram, thereby obtaining a removed differential interferogram; Based on all the removed differential interferograms, a LOS deformation data set of the research region is calculated. 3.The method according to claim 1, wherein, The three-dimensional solution model is: D Los = D v cos θ + D N sin θ sin β - D E sin θ cos β; where D Los is the LOS deformation; D v is the vertical deformation; θ is the radar incidence angle; D N is the north-south deformation; β is the radar heading angle; D E is the east-west deformation; The SPFM model is: where v U is the vertical deformation rate; is the first derivative in the north direction, H is the terrain height, X N is the north deformation component; v N is the north-south deformation rate; is the first derivative in the east direction, X E is the east deformation component; v E is the east-west deformation rate. 4.The method according to claim 1, wherein, Based on the three-dimensional deformation rate data set, a first influencing factor is detected from the first type of potential influencing factors in the time scale by using the correlation coefficient method, and a second influencing factor is detected from the second type of potential influencing factors in the spatial scale by using the geographic detector, which specifically includes: obtaining the factor data set of each first type of potential influencing factor and the value of each second type of potential influencing factor of the research region; the factor data set includes the value of the first type of potential influencing factor at each first time 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 the vertical deformation rate data at each first time in a preset time period; using the factor data set of each first type of potential influencing factor and the vertical deformation rate data set as input, a first influencing factor is detected from the first type of potential influencing factor in the time scale by using the correlation coefficient method; using the value of each second type of potential influencing factor and the value of the vertical deformation rate as input, a second influencing factor is detected from the second type of potential influencing factor in the spatial scale by using the geographic detector; the value of the vertical deformation rate is the vertical deformation rate data in the vertical deformation rate data set at a random first time or the average value of the vertical deformation rate data at each first time in the vertical deformation rate data set. 5.The method according to claim 4, wherein, The first influencing factor is detected from the first type of potential influencing factor in a time scale by using a correlation coefficient method, with the factor data set of each first type of potential influencing factor and the vertical deformation rate data set as inputs, specifically including: The principal component analysis method is used to analyze the vertical deformation rate data set, and a principal component data set corresponding to each of a plurality of principal components is obtained; the principal component data set includes the value of the principal component at each first time in a preset time period; For each first type of potential influencing factor, the correlation coefficient between the first type of potential influencing factor and each principal component is calculated, with the factor data set corresponding to the first type of potential influencing factor and the principal component data set corresponding to each principal component as inputs; The first influencing factor is detected from the first type of potential influencing factor based on all correlation coefficients. 6.The method according to claim 4, wherein, The factor data set of each first type of potential influencing factor and the value of each second type of potential influencing factor of the research area are obtained, specifically including: The glacier area of the research area is extracted by using a multi-source remote sensing data set, and the glacier area is removed in the research area to obtain the factor data set of each first type of potential influencing factor and the value of each second type of potential influencing factor of the non-glacier area of the research area; Wherein, the formula for extracting the glacier area is: GlacierArea = <Ave {Image E S2: [(D1 > D1 * ) A (D2 < D2 * ) A (NDSI > NDSI * )]} > 0.5>; where GlacierArea is the glacier area; S2 is a multi-source remote sensing dataset; D1 is a first index; is a threshold value for the first index; D2 is a second index; is a threshold value for the second index; NDSI is a normalized difference snow index; NDSI * is a threshold value for the normalized difference snow index; Wherein, Red is the reflectivity of the red band; SWIR1 is the reflectivity of the infrared band; Wherein, Blue is the reflectivity of the blue band. 7.The method according to claim 2, wherein, The relief phase and the flat phase are calculated by using the converted DEM data, and the relief phase and the flat phase are removed in each differential interferogram to obtain a plurality of preliminary removed differential interferograms, specifically including: The relief phase and the flat phase are calculated by using the converted DEM data; For each differential interferogram, the differential interferogram is phase-optimized three times by using a Goldstein filter to obtain a denoised differential interferogram; wherein the windows used for three times of phase optimization are 64*64, 32*32 and 16*16 respectively; The relief phase and the flat phase are removed in each denoised differential interferogram to obtain a plurality of preliminary removed differential interferograms; A plurality of high-coherence points are determined, specifically including: a plurality of high-coherence points are determined by using an amplitude dispersion method and a coherence coefficient method; Each high-coherence point is phase-unwrapped based on the reference point, specifically including: each high-coherence point is phase-unwrapped based on the reference point by using a minimum cost flow method; The formula for calculating the atmospheric delay phase is: wherein, is the atmospheric delay phase; is the difference data of the Mth converted atmospheric data and the Nth converted atmospheric data, i.e. the two converted atmospheric data; is the atmospheric delay phase of the difference data of the Mth converted atmospheric data and the Nth converted atmospheric data at the reference point (i,j); and λ is the wavelength of the SAR sensor.

8. A computer device comprising: Memory, processor, and computer program stored on the memory and executable on the processor, characterized in that the processor executes the computer program to implement the spatiotemporal scale landslide disaster inducing factor detection method of any one of claims 1-7.

9. A computer-readable storage medium having stored thereon a computer program, characterized in that, The computer program is executed by the processor to implement the spatiotemporal scale landslide disaster inducing factor detection method of any one of claims 1-7.

10. A computer program product comprising a computer program, characterized in that, The computer program is executed by the processor to implement the spatiotemporal scale landslide disaster inducing factor detection method of any one of claims 1-7.

Images