A method and system for monitoring deformation rate of large-scale landslides in mountainous areas using InSAR

Abstract

The present invention discloses a method and system for monitoring the deformation rate of a large-scale landslide in a mountainous area by InSAR, which relates to the field of synthetic aperture radar interferometry technology, and solves the problem of low accuracy of the existing method for calculating the deformation rate of landslide InSAR. The method comprises: using a mixed baseline strategy to generate a first interference phase map and a second interference phase map; using a time series phase optimization algorithm to optimize the first interference phase map to obtain a first target interference phase map; using a Gaussian filter algorithm to optimize the second interference phase map to obtain a second target interference phase map; combining the first target interference phase map and the second target interference phase map to obtain a mixed interference phase map; performing phase unwrapping processing and time series analysis on the mixed interference phase map to obtain the surface deformation rate of the target mountainous area. The method for monitoring the deformation rate of a large-scale landslide in a mountainous area provided by the present invention significantly improves the accuracy and coverage of InSAR deformation rate monitoring of landslides in mountainous areas.

Description

Technical Field

[0001] The invention relates to the field of synthetic aperture radar interferometry, and in particular to an InSAR deformation rate monitoring method and system for large-scale landslides in mountainous areas. Background Art

[0002] Landslide deformation monitoring can obtain the deformation of the landslide throughout the process. This information can be used to design landslide prevention and control projects to reduce the occurrence of landslide disasters. InSAR technologies such as PS-InSAR (Persistent Scatterer InSAR), SBAS-InSAR (Small Baseline Subset InSAR), and Distributed Scatterer InSAR (DS-InSAR) are commonly used for landslide monitoring. However, in complex mountainous areas, special terrain features will lead to interference factors such as terrain effects, the influence of unstable scatterers, and atmospheric delays in the area. These interference factors will cause phase distortion or dislocation of the collected InSAR images, thereby affecting the accuracy and reliability of the deformation signal. Therefore, the phase needs to be optimized during landslide deformation monitoring.

[0003] The DS-InSAR method is commonly used for phase optimization. However, the existing DS-InSAR method only obtains the interference phase image through a single baseline combination method, and then performs phase optimization based on the time series. This optimization method easily leads to data scarcity and difficulty in providing comprehensive deformation monitoring information, resulting in low accuracy of the landslide InSAR deformation rate calculation results, and unable to accurately carry out early warning and emergency response to geological disasters such as landslides.

[0004] Therefore, how to improve the application capability of InSAR in large-scale landslide deformation monitoring and overcome the problems of data scarcity and insufficient information has become one of the urgent problems to be solved in the current InSAR field. Summary of the invention

[0005] The purpose of the present invention is to provide a method and system for InSAR deformation rate monitoring of large-scale landslides in mountainous areas, which is used to solve the problems that the optimized phase data obtained by the existing landslide deformation monitoring methods are scarce and the information is insufficient, and the accuracy of the determined landslide InSAR deformation rate is low, so as to improve the accuracy and coverage of InSAR deformation rate monitoring of landslides in mountainous areas.

[0006] In order to achieve the above object, the present invention provides the following technical solutions:

[0007] In a first aspect, the present invention provides a method for monitoring deformation rate of large-scale landslides in mountainous areas using InSAR, comprising:

[0008] Obtain N SAR images of the target mountain area; N is an integer greater than 1;

[0009] The N-scene SAR images are interferometrically processed by a single main image combination method to obtain a first interferometric phase map; the N-scene SAR images are interferometrically processed by a short baseline combination method to obtain a second interferometric phase map;

[0010] The first interference phase image is optimized by using a sequential phase optimization algorithm to obtain a first target interference phase image;

[0011] Using a Gaussian filtering algorithm to optimize the second interference phase image to obtain a second target interference phase image;

[0012] Combining the first target interference phase image with the second target interference phase image to obtain a mixed interference phase image;

[0013] The mixed interference phase image is subjected to phase unwrapping processing and time series analysis to obtain the surface deformation rate of the target mountain area.

[0014] Optionally, the step of optimizing the second interference phase image by using a Gaussian filtering algorithm to obtain a second target interference phase image includes:

[0015] Using the formula:

[0016] ;

[0017] Performing Gaussian filtering on the real part of the second interference phase image to obtain a first Gaussian filtering result, and performing Gaussian filtering on the imaginary part of the second interference phase image to obtain a second Gaussian filtering result;

[0018] in, is the first Gaussian filtering result or the second Gaussian filtering result, , is the pixel value of the real part or the pixel value of the imaginary part of the second interference phase image input, For the current pixel The horizontal offset of For the current pixel The vertical offset of is the standard deviation of the Gaussian kernel;

[0019] The first Gaussian filtering result and the second Gaussian filtering result are compositely processed to obtain a second target interference phase map.

[0020] Optionally, the interferometric processing of the N SAR images by using a single main image combination method to obtain a first interferometric phase map includes:

[0021] The main image and the auxiliary image in the N SAR images are registered to obtain a plurality of first interference image pairs; the main image is the N / 2th image in the N SAR images, and the auxiliary image is an image other than the main image;

[0022] Using the formula:

[0023] ;

[0024] Calculating a first interference phase image of the first interference image pair;

[0025] in, is the first interference phase diagram, and are the two images of the first interference image pair, is the terrain phase compensation term of the first interferometric image pair expressed in complex form, is a mathematical constant, is the product of the amplitudes of the two images of the first interference image pair, For Take the complex conjugate.

[0026] Optionally, use the formula:

[0027] ;

[0028] Intercepting the N SAR images to obtain a target image;

[0029] in, is the minimum time baseline, is the maximum time baseline;

[0030] Combining images at adjacent time points in the target image in pairs to obtain a second interference image pair;

[0031] The second interference image pair is subjected to interference processing to obtain a second interference phase map.

[0032] Optionally, the step of optimizing the first interference phase image by using a timing phase optimization algorithm to obtain a first target interference phase image includes:

[0033] Selecting target homogeneous pixels in the first interferometric phase map according to the intensity data of the N SAR images;

[0034] Performing stack processing and matrix conversion processing on the target homogeneous pixels to obtain a covariance matrix and an initial coherence matrix;

[0035] Performing a positive definiteness test on the initial coherence matrix to obtain a target coherence matrix that meets the positive definiteness requirement;

[0036] Respectively calculating a first eigenvalue corresponding to the covariance matrix and a second eigenvalue corresponding to the target coherence matrix;

[0037] An eigenvector corresponding to the minimum eigenvalue of the first eigenvalue and the second eigenvalue is determined as a first target interference phase map.

[0038] Optionally, selecting target homogeneous pixels in the first interference phase map according to the intensity data of the N SAR images includes:

[0039] Using the formula:

[0040] ;

[0041] Determine the first confidence interval estimate;

[0042] in, is the time series average intensity of the pixel to be detected, is the time series average intensity of the reference pixel; intensity ratio The degree of freedom is F distribution; is the F distribution at the significance level The upper quantile of the lower is the F distribution at the significance level the lower quantile of the lower;

[0043] Confirming the pixels in the first window in the first confidence interval estimate as an initial homogeneous pixel set;

[0044] Calculating an average value of pixel intensities in the initial set of homogeneous pixels, and determining the average value as an estimated value of the reference pixel;

[0045] Using the formula:

[0046] ;

[0047] determining a second confidence interval estimate;

[0048] in, For the standard gamma distribution at a significant level The upper quantile of the lower For the standard gamma distribution at a significant level The lower quantile point, is the estimated value of the reference pixel;

[0049] The pixels in the second confidence interval estimate are determined as target homogeneous pixels.

[0050] Optionally, the stacking and matrix conversion processing is performed on the target homogeneous pixels to obtain a covariance matrix and an initial coherence matrix, including:

[0051] Using the formula:

[0052] ;

[0053] Loading the first interference phase image into a three-dimensional phase stack to obtain a three-dimensional phase image;

[0054] in, is the pixel position in the first interference phase image of scene N The phase of It is a three-dimensional phase stack;

[0055] extracting a homogeneous pixel stack within a third window at each pixel position in the three-dimensional phase image;

[0056] reshaping the homogeneous pixel stack into a homogeneous pixel matrix;

[0057] Using the formula:

[0058] ;

[0059] Converting the homogeneous pixel matrix into a covariance matrix;

[0060] in, is the covariance matrix, is a homogeneous pixel stack, is the number of target homogeneous pixel stacks in the third window;

[0061] The absolute values ​​in the covariance matrix are determined as an initial coherence matrix.

[0062] Optionally, performing positive definiteness detection on the initial coherence matrix to obtain a target coherence matrix that meets the positive definiteness includes:

[0063] By performing Cholesky decomposition on the initial coherence matrix, determining whether the initial coherence matrix meets the positive definiteness;

[0064] When the initial coherence matrix does not meet the positive definiteness, adding a target value on the diagonal of the initial coherence matrix to obtain a new coherence matrix;

[0065] It is determined whether the new coherence matrix meets the positive definiteness. If the new coherence matrix does not meet the positive definiteness, the target value is increased until the new coherence matrix meets the positive definiteness, thereby obtaining the target coherence matrix.

[0066] Optionally, performing phase unwrapping processing and time series analysis on the hybrid interference phase image to obtain the surface deformation rate of the target mountain area includes:

[0067] Perform atmospheric correction on the mixed interference phase diagram to obtain a corrected interference phase diagram;

[0068] Use the formula:

[0069] ;

[0070] Perform phase unwrapping on the corrected interference phase diagram to obtain the unwrapped phase;

[0071] where, is the unwrapped phase, is the corrected interference phase diagram, represents the number of phase jumps;

[0072] Use a linear regression model to fit the unwrapped phase to obtain a fitted interference phase function:

[0073] Use the least squares method to solve the fitted interference phase function to obtain the deformation rate and the initial phase.

[0074] Compared with the prior art, the present invention provides an InSAR deformation rate monitoring method for large-scale landslides in mountainous areas, comprising: acquiring N SAR images of a target mountainous area; interferometrically processing the N SAR images by a single main image combination method to obtain a first interference phase map; interferometrically processing the N SAR images by a short baseline combination method to obtain a second interference phase map; optimizing the first interference phase map by a time series phase optimization algorithm to obtain a first target interference phase map; optimizing the second interference phase map by a Gaussian filtering algorithm to obtain a second target interference phase map; combining the first target interference phase map and the second target interference phase map to obtain a mixed interference phase map; performing phase unwrapping processing and time series analysis on the mixed interference phase map to obtain the surface deformation rate of the target mountainous area, thereby realizing large-scale landslide deformation monitoring. By optimizing the time series phase of the first interferometric phase image obtained by the single main image combination method, the spatial coverage of deformation monitoring can be improved, and a higher point density and a wider coverage range are displayed in the low coherence area. At the same time, the short baseline combination method is introduced on the basis of the single main image combination method, and the number of interferometric phase images is increased. At the same time, the second interferometric phase image obtained by the short baseline combination method is Gaussian filtered, and a significantly higher number of measurement points and spatial density can be obtained, which fully captures the complex spatiotemporal deformation characteristics of mountain landslides, provides more accurate deformation information, and obtains a more reliable and accurate deformation rate. The present invention effectively solves the limitations of the traditional DS-InSAR method in the single baseline combination method and the time series phase optimization, improves the application ability of InSAR technology in large-scale landslide deformation monitoring, and overcomes the problems of data scarcity and insufficient information. It has important technical value and broad application scenarios, especially in surface deformation monitoring and geological disaster monitoring in complex mountainous environment areas, showing significant potential and advantages.

[0075] In a second aspect, the present invention further provides an InSAR deformation rate monitoring system for large-scale landslides in mountainous areas, comprising:

[0076] The SAR image acquisition module is used to acquire N SAR images of the target mountain area; N is an integer greater than 1;

[0077] A hybrid baseline processing module is used to perform interference processing on the N-scene SAR images by adopting a single main image combination method to obtain a first interference phase map; and to perform interference processing on the N-scene SAR images by adopting a short baseline combination method to obtain a second interference phase map;

[0078] A timing optimization module, used to optimize the first interference phase image by using a timing phase optimization algorithm to obtain a first target interference phase image;

[0079] A Gaussian filter algorithm optimization module, used to optimize the second interference phase image by using a Gaussian filter algorithm to obtain a second target interference phase image;

[0080] An interference phase combination module, used for combining the first target interference phase image with the second target interference phase image to obtain a mixed interference phase image;

[0081] The phase processing and surface deformation analysis module is used to perform phase unwrapping processing and timing analysis on the mixed interference phase image to obtain the surface deformation rate of the target mountain area.

[0082] Compared with the prior art, the beneficial effects of the InSAR deformation rate monitoring system for large-scale landslides in mountainous areas provided by the present invention are the same as the beneficial effects of the InSAR deformation rate monitoring method for large-scale landslides in mountainous areas described in the above technical solution, which will not be repeated here. BRIEF DESCRIPTION OF THE DRAWINGS

[0083] The drawings described herein are used to provide a further understanding of the present invention and constitute a part of the present invention. The exemplary embodiments of the present invention and their descriptions are used to explain the present invention and do not constitute an improper limitation of the present invention. In the drawings:

[0084] Figure 1 A flow chart of a method for monitoring deformation rate of large-scale landslides in mountainous areas provided by the present invention using InSAR;

[0085] Figure 2 The original interference phase map provided by the present invention is composed of a first interference phase map and a second interference phase map with image times of January 2, 2024 and February 19, 2024;

[0086] Figure 3 A schematic diagram of a first target interference phase diagram provided by the present invention;

[0087] Figure 4 A schematic diagram of a second target interference phase diagram provided by the present invention;

[0088] Figure 5 The present invention provides Figure 3 A magnified image of area A;

[0089] Figure 6 The present invention provides Figure 4 Enlarged view of area B;

[0090] Figure 7 A deformation rate map obtained by using the PS-InSAR monitoring method provided by the present invention;

[0091] Figure 8 A deformation rate map obtained by using the DS-InSAR monitoring method provided by the present invention;

[0092] Fig. 9 A deformation rate map obtained by using the SBAS-InSAR monitoring method provided by the present invention;

[0093] Fig.10 A deformation rate diagram obtained by using the method of the present invention provided by the present invention;

[0094] Fig.11 A structural schematic diagram of an InSAR deformation rate monitoring system for large-scale landslides in mountainous areas provided by the present invention. DETAILED DESCRIPTION

[0095] In order to clearly describe the technical solutions of the embodiments of the present invention, in the embodiments of the present invention, words such as "first" and "second" are used to distinguish the same items or similar items with basically the same functions and effects. For example, the first threshold and the second threshold are only used to distinguish different thresholds, and their order is not limited. Those skilled in the art can understand that words such as "first" and "second" do not limit the quantity and execution order, and words such as "first" and "second" do not necessarily limit them to be different.

[0096] It should be noted that, in the present invention, words such as "exemplary" or "for example" are used to indicate examples, illustrations or descriptions. Any embodiment or design described as "exemplary" or "for example" in the present invention should not be interpreted as being more preferred or more advantageous than other embodiments or designs. Specifically, the use of words such as "exemplary" or "for example" is intended to present related concepts in a specific way.

[0097] In the present invention, "at least one" means one or more, and "plurality" means two or more. "And / or" describes the association relationship of associated objects, indicating that three relationships may exist. For example, A and / or B can mean: A exists alone, A and B exist at the same time, and B exists alone, where A and B can be singular or plural. The character " / " generally indicates that the objects associated before and after are in an "or" relationship. "At least one of the following" or similar expressions refers to any combination of these items, including any combination of single or plural items. For example, at least one of a, b or c can mean: a, b, c, the combination of a and b, the combination of a and c, the combination of b and c, or the combination of a, b and c, where a, b, c can be single or multiple.

[0098] Before introducing the embodiments of the present invention, the following definitions are given for the relevant terms involved in the embodiments of the present invention:

[0099] Cholesky decomposition is a way to represent a symmetric positive definite matrix A as a lower triangular matrix L and its transpose LT. This method requires that all eigenvalues ​​of the matrix A must be greater than zero, thus ensuring that the diagonal elements of the decomposed lower triangular matrix L are also positive.

[0100] With the increasing number of SAR satellites, the interferometric measurement technology of synthetic aperture radar (InSAR) has been widely used in the fields of surface deformation monitoring and plays an important role. Especially in complex mountainous areas, especially in areas where landslides occur frequently, InSAR technology can provide large-scale, long-term and stable deformation monitoring data with its all-weather and all-time advantages, and has become an important monitoring method. However, in complex mountainous areas, due to the special terrain characteristics, InSAR technology faces many challenges, mainly including terrain effects, the influence of unstable scatterers and atmospheric delay, which can cause phase distortion or dislocation. For example, Dongchuan District, Kunming City, Yunnan Province, China, forms a "V" shape with Xiaojiang as the boundary. The highest peak is 4344.1 meters above sea level. Precipitation is mostly concentrated in the rainy period from May to October. Geological disasters occur frequently, and landslides and mudslides are particularly serious. Dongchuan has a complex geological structure, diverse strata and lithology, and dense and long-term active fault zones, resulting in the distribution of a large number of rock fracture zones and mylonites, making the surface deformation greatly affected by the complex geological structure and strong surface undulations. The area is known as a natural museum of debris flows, and the region is rich in mineral resources. Long-term mining has further aggravated the frequency of landslides. For this type of area with complex terrain, vegetation coverage and frequent geological disasters, traditional InSAR methods are difficult to fully capture the deformation information of the area, and cannot accurately monitor small-scale deformations in deeply cut alpine canyon areas. Therefore, it is urgent to improve InSAR technology and algorithms to adapt to such extremely complex terrain conditions.

[0101] The existing DS-InSAR method has improved the traditional InSAR technology. It adopts a single baseline combination strategy to obtain the interferogram and uses a time series to optimize the phase. This method easily leads to data scarcity, making it difficult to provide comprehensive deformation monitoring information, and the calculated deformation rate has low accuracy.

[0102] In order to solve the above problems, the present invention provides a method and system for monitoring the deformation rate of large-scale landslides in mountainous areas by InSAR, which will be described below with reference to the accompanying drawings.

[0103] See also Figure 1 The present invention provides a method for monitoring deformation rate of large-scale landslides in mountainous areas by InSAR, comprising the following steps:

[0104] Step S1: Acquire N SAR images of the target mountain area;

[0105] N is an integer greater than 1; N is an even number. N SAR images contain multiple radar images, which are acquired at different times.

[0106] Step S2: performing interference processing on the N-scene SAR images by adopting a single main image combination method to obtain a first interference phase map; performing interference processing on the N-scene SAR images by adopting a short baseline combination method to obtain a second interference phase map;

[0107] The single main image combination method is a strategy in which all auxiliary images interfere with the main image, wherein the N / 2th image is selected as the main image and the remaining images are used as auxiliary images.

[0108] The short baseline combination method is a strategy of selecting images at adjacent time points for interference based on the time baseline between images.

[0109] The interferometric processing of the N SAR images by using a single main image combination method to obtain a first interferometric phase map specifically includes:

[0110] The main image and the auxiliary image in N SAR images are registered to obtain multiple groups of first interference image pairs; the main image is the N / 2th image in the N SAR images, and the auxiliary image is the image other than the main image; each auxiliary image forms a first interference image pair with the main image.

[0111] Using formula (1):

[0112] (1)

[0113] Calculating a first interference phase image of a first interference image pair;

[0114] in, is the first interference phase diagram, and are the two images of the first interference image pair, is the terrain phase generated by the two images of the first interferometric image pair, is the terrain phase compensation term of the first interferometric image pair expressed in complex form. By multiplying it, the phase influence caused by the terrain can be removed, so that the final differential phase can focus more on the deformation information. is a mathematical constant, the base of natural logarithms, represents the imaginary part of a complex number, is the product of the amplitudes of the two images of the first interferometric image pair, which is used to normalize the phase difference result. For Take the complex conjugate.

[0115] The above-mentioned interferometric processing of the N SAR images by using the short baseline combination method to obtain the second interferometric phase map specifically includes:

[0116] Using formula (2):

[0117] (2)

[0118] Intercept N SAR images to obtain the target image; the target image has a time baseline of arrive Images within range.

[0119] in, is the minimum time baseline, is the maximum time baseline;

[0120] Combining images at adjacent time points in the target image in pairs to obtain a second interference image pair;

[0121] Using formula (3):

[0122] (3)

[0123] performing interference processing on the second interference image pair to obtain a second interference phase image;

[0124] in, is the second interference phase pattern, and are the two images of the second interference image pair, is the terrain phase generated by the two images of the second interferometric image pair, is the terrain phase compensation term of the second interferometric image pair expressed in complex form, is a mathematical constant, is the product of the amplitudes of the two images of the second interference image pair, For Take the complex conjugate.

[0125] Step S3: optimizing the first interference phase image by using a sequential phase optimization algorithm to obtain a first target interference phase image;

[0126] Step S3 can be implemented by the following steps:

[0127] Step S31: on the basis of the first interferometric phase map, selecting target homogeneous pixels in the first interferometric phase map according to the intensity data of N SAR images;

[0128] Assume that N SAR images obey the complex circular Gaussian distribution in the time dimension, and their intensity data obey the exponential distribution; according to the relationship between the exponential distribution and the gamma distribution, the probability density function of the intensity can be described by the gamma distribution. Within the first window, for the time series average intensity of the pixel to be detected and the time series average intensity of the reference pixel, the intensity ratio obeys the degree of freedom The F distribution of , so the first confidence interval estimate is as shown in formula (4):

[0129] (4)

[0130] Specifically, step S31 includes: using formula (4) to determine a first confidence interval estimate;

[0131] in, is the time series average intensity of the pixel to be detected, is the time series average intensity of the reference pixel; intensity ratio The degree of freedom is F distribution; is the F distribution at the significance level The upper quantile of the lower is the F distribution at the significance level the lower quantile of the lower;

[0132] Confirming the pixels in the first window in the first confidence interval estimation as the initial homogeneous pixel set; illustratively, the first window is a 7×7 window;

[0133] 5. Calculate the average pixel intensity of the pixels in the initial homogeneous pixel set, and determine the average value as the estimated value of the reference pixel; on this basis, further expand to the second window, for example: 15×15 window, and perform weighted gamma distribution estimation of the sample mean, so as to construct a more accurate confidence interval, as shown in formula (5):

[0134] (5)

[0135] Formula (5) is used to determine the second confidence interval estimate;

[0136] in, For the standard gamma distribution at a significant level The upper quantile of the lower For the standard gamma distribution at a significant level The lower quantile point, is the estimated value of the reference pixel;

[0137] The pixels within the second window in the second confidence interval estimate are determined as target homogeneous pixels.

[0138] Step S32: performing stack processing and matrix conversion processing on the target homogeneous pixels to obtain a covariance matrix and an initial coherence matrix;

[0139] Specifically, use formula (6):

[0140] (6)

[0141] Loading the first interference phase image into the three-dimensional phase stack to obtain a three-dimensional phase image;

[0142] in, is the pixel position in the first interference phase image of scene N The phase of It is a three-dimensional phase stack;

[0143] At each pixel position in the 3D phase image Select one A third window is formed, and a target homogeneous pixel stack within the third window is extracted;

[0144] Reshape the target homogeneous pixel stack into a target homogeneous pixel matrix , which is convenient for subsequent covariance matrix calculation; the target homogeneous pixel matrix It is in the form of the number of first interference phase images × the number of target homogeneous pixel stacks in the third window.

[0145] Using formula (7):

[0146] (7)

[0147] Convert homogeneous pixel matrix to covariance matrix;

[0148] in, is the covariance matrix, yes The conjugate transpose of is the number of target homogeneous pixel stacks in the third window, and the covariance matrix represents the phase consistency between target homogeneous pixels in different time series.

[0149] The absolute value is extracted from the covariance matrix, and the absolute value in the covariance matrix is ​​determined as the initial coherence matrix, as shown in formula (8):

[0150] (8)

[0151] in, is the initial coherence matrix; the initial coherence matrix is ​​used to evaluate the phase coherence between the first interference phase images, thereby providing a reference for subsequent phase optimization.

[0152] Step S33: Perform positive definiteness test on the initial coherence matrix to obtain a target coherence matrix that meets the positive definiteness requirement;

[0153] Specifically, by performing Cholesky decomposition on the initial coherence matrix, it is determined whether the initial coherence matrix meets the positive definiteness; if the decomposition is successful, the initial coherence matrix meets the positive definiteness, and if the decomposition fails, the initial coherence matrix does not meet the positive definiteness;

[0154] When the initial coherence matrix does not meet the positive definiteness, the target value is added to the diagonal of the initial coherence matrix to obtain a new coherence matrix; the new coherence matrix is ​​shown in formula (9):

[0155] (9)

[0156] in, is the new coherence matrix, is the identity matrix, is the target value, which is a very small positive number, such as 1, 2, 3, etc.

[0157] Determine whether the new coherence matrix meets the positive definiteness. If the new coherence matrix does not meet the positive definiteness, gradually increase the value of the target numerical value until the new coherence matrix meets the positive definiteness, and obtain the target coherence matrix to ensure the positive definiteness of the coherence matrix.

[0158] Step S34: Finally, the covariance matrix and the target coherence matrix are used to perform phase optimization on the first interference phase image to obtain a first target interference phase image.

[0159] Specifically, the eigenvector corresponding to the minimum eigenvalue is calculated using formula (10):

[0160] (10)

[0161] Solving formula (9) yields formula (11), and obtaining the first target interference phase diagram:

[0162] (11)

[0163] in, is the first target interference phase diagram, Indicates the selection of the eigenvector that minimizes the eigenvalue , represents the eigendecomposition, is the index of the main image.

[0164] Step S4: optimizing the second interference phase image by using a Gaussian filtering algorithm to obtain a second target interference phase image;

[0165] The principle of the Gaussian filtering algorithm is to select the size and standard deviation parameters of the Gaussian filtering kernel, smooth the phase noise, and finally obtain the interference phase map after Gaussian filtering optimization.

[0166] The core idea of ​​the Gaussian filter algorithm is to smooth the image by weighted averaging. Specifically, the new value of each pixel is the weighted average of all pixels in its neighborhood, where the weight is determined by a two-dimensional Gaussian function. The Gaussian filter convolution formula is shown in formula (12):

[0167] (12)

[0168] in, is a Gaussian filter, is the pixel value of the real part or the pixel value of the imaginary part of the second interference phase image input, is a two-dimensional Gaussian kernel function, representing the pixel position of the second interference phase The weight of the corresponding pixel, , For the current pixel The horizontal offset of For the current pixel The vertical offset of is the standard deviation of the Gaussian kernel, which is used to control the smoothness of the filter. It can be 0.9.

[0169] The second interference phase pattern is complex data, as shown in formula (13):

[0170] (13)

[0171] in, is the real part of the second interference phase pattern, is the imaginary part of the second interference phase pattern.

[0172] Therefore, the real part of the second interference phase image is subjected to Gaussian filtering using formula (12) to obtain the first Gaussian filtering result. The imaginary part of the second interference phase image is subjected to Gaussian filtering using formula (12) to obtain the second Gaussian filtering result. The first Gaussian filtering result and the second Gaussian filtering result are composited to obtain the second target interference phase image.

[0173] Step S5: combining the first target interference phase image and the second target interference phase image to obtain a mixed interference phase image;

[0174] It should be noted that the combination means that the first target interference phase image and the second target interference phase image are grouped into a set, and the set is called a mixed interference phase image.

[0175] Step S6: performing phase unwrapping processing and timing analysis on the hybrid interference phase image to obtain the surface deformation rate of the target mountain area.

[0176] Specifically, the atmospheric correction process of the mixed interferometric phase image is performed using the GACOS atmospheric correction product to obtain a corrected interferometric phase image;

[0177] Using formula (14):

[0178] (14)

[0179] Performing phase unwrapping processing on the corrected interference phase image to obtain an unwrapped phase;

[0180] in, is the unwrapped phase, To correct the interferometric phase pattern, Indicates the number of phase jumps, For time.

[0181] The unwrapped phase is fitted using a linear regression model to obtain a fitted interference phase function, which is shown in formula (15):

[0182] (15)

[0183] in, is the surface deformation rate, is the initial phase and t is the time.

[0184] The least square method is used to solve the fitted interference phase function to obtain the deformation rate and initial phase.

[0185] The goal of the least squares method is to minimize the optimization objective function, which is calculated as shown in formula (16):

[0186] (16)

[0187] right and Taking the derivative and setting it to zero, we get the optimal solution, as shown in formula (17):

[0188] (17)

[0189] in, is the average of the unwrapped phase, is the average value over time. Through the least squares method, the surface deformation rate of each pixel can be obtained, that is, , is the phase change, is the time variation.

[0190] pass Figure 1 It can be seen from the method that the present invention can improve the spatial coverage of deformation monitoring by optimizing the time series phase of the first interference phase map obtained by the single main image combination method, and show a higher point density and a wider coverage range in the low coherence area. At the same time, the short baseline combination method is introduced on the basis of the single main image combination method, and the number of interference phase maps is increased. At the same time, the second interference phase map obtained by the short baseline combination method is Gaussian filtered, which can obtain a significantly higher number of measurement points and spatial density, fully capture the complex spatiotemporal deformation characteristics of mountain landslides, provide more accurate deformation information, and obtain a more reliable and accurate deformation rate. The present invention effectively solves the limitations of the traditional DS-InSAR method in the single baseline combination method and the time series phase optimization, improves the application ability of InSAR technology in large-scale landslide deformation monitoring, and overcomes the problems of data scarcity and insufficient information. It has important technical value and broad application scenarios, especially in surface deformation monitoring and geological disaster monitoring in complex mountainous environment areas, showing significant potential and advantages. In addition, atmospheric correction of the mixed interference phase map can reduce the influence of atmospheric delay,

[0191] See also Figure 2-Figure 10 The present invention takes the monitoring of large-scale landslide deformation in Dongchuan District as an example to explain the effect of the present invention in detail. The specific implementation process includes the following steps:

[0192] First, the C-band Sentinel-1A SAR radar image data launched by the European Space Agency in 2014 was selected. The images were acquired in IW mode, and the time interval was 42 scenes of descending orbit data from May 7, 2023 to September 22, 2024. The 21st image was selected as the main image, and the remaining images were registered with it as auxiliary images. The single main image combination method and the short baseline combination method were used for interference processing of the registered SAR images. The original interference phase images of the image time January 2, 2024 and February 19, 2024 obtained by interference processing using the single main image combination method and the short baseline combination method are shown as follows: Figure 2 shown.

[0193] Based on the interference pattern generated by the single main image combination method, the intensity data of 42 SAR images are used to select homogeneous pixels, the covariance matrix of the selected homogeneous pixels is estimated and the coherence matrix is ​​calculated, the coherence matrix is ​​tested for positive definiteness, and finally the covariance matrix and the coherence matrix are used to optimize the interferometric phase in a sequential manner. The optimized interferometric phase patterns on January 2, 2024 and February 19, 2024 obtained by the single main image combination method using the sequential phase optimization method are shown in the figure below. Figure 3 shown.

[0194] Based on the interferogram generated by the short baseline combination method, the Gaussian filter algorithm is applied to the initial interferometric phase data in the interferogram, that is, the Gaussian filter kernel size is selected as 3 and the standard deviation parameter is selected as 0.9 to smooth the phase noise. Finally, the interferometric phase optimized by Gaussian filtering is obtained, as shown in Figure 4 As shown. The optimized single main baseline and short baseline interferometric phases are fused to construct a hybrid interferometric phase;

[0195] and Figure 2 In comparison, the optimized interference phase can effectively overcome the pseudo signal caused by spatial jump and make the phase more continuous. Both optimization methods can maintain spatial continuity, such as Figure 5 As shown in , the optimized interference phase of the covariance matrix has many discrete jump interference phase pseudo signals, whether in the area with low noise level or in the low coherence area, such as Figure 6 As shown in the figure, the optimized interference phase after Gaussian filtering is obviously more continuous, without obvious discrete jump interference phase pseudo signal, and the overall phase fringes are smoother and of better quality.

[0196] The optimized single main baseline and short baseline interferometric phases are fused to construct a hybrid interferometric phase.

[0197] The interferometric phases optimized by the single main baseline combination method and the short baseline combination method constitute a hybrid interferometric phase, and the atmospheric correction of the hybrid interferometric phase is performed using the GACOS atmospheric correction product to reduce the influence of atmospheric delay. The interferometric phase after atmospheric correction is phase unwrapped to obtain continuous phase information, and the surface deformation rate is estimated by the least squares method using the time series analysis method, thereby monitoring large-scale landslide deformation.

[0198] See also Figure 7-10 The application effect of the present invention in landslide monitoring in complex mountainous areas is evaluated by conducting a comparative analysis of the deformation rate results, the number of measurement points (MPs) and the spatial density of HGDS-InSAR, PS-InSAR, SBAS-InSAR and DS-InSAR. Among them, HGDS-InSAR is the method corresponding to the present invention, such as Figure 7 As shown in Figure 1, the extraction of temporal deformation information based on PS points of permanent scatterers with high coherence is suitable for high coherence areas such as cities. However, in complex mountainous areas, the density of PS points is low, the spatial coverage of deformation monitoring is limited, and it is difficult to fully characterize the activity characteristics of geological disasters such as landslides. Figure 8As shown in the figure, based on PS-InSAR, the temporal phase optimization is performed on the combination of multiple interferograms of a single main image to improve the spatial coverage of deformation monitoring. Although the DS-InSAR method shows higher point density and wider coverage in low coherence areas, it is limited by the baseline combination method of a single main image and the limited number of interferograms, which cannot fully capture the complex spatiotemporal deformation characteristics of mountain landslides. Fig. 9 As shown in the figure, increasing the number of interferometer pairs by short baseline combination significantly improves the coverage and monitoring accuracy of low coherence areas. However, since SBAS-InSAR does not perform phase optimization on the short baseline interferogram, the interferometric signal-to-noise ratio and phase accuracy in some areas are still insufficient, especially in areas with frequent landslide activities, where there is a certain degree of deformation omission. The limitation of this method is that, despite increasing the number of interferometer pairs, no further optimization can be made in terms of interferogram quality and temporal phase consistency. Fig.10 As shown, a short baseline combination is introduced based on the single main baseline combination method of DS-InSAR, and Gaussian filtering is performed on each pair of short baseline interferograms, and the coherence and phase accuracy are improved through time-series phase optimization. Compared with traditional methods, the method HGDS-InSAR provided by the present invention significantly improves the coverage and accuracy of deformation monitoring, especially in complex mountainous areas. HGDS-InSAR can obtain higher quality time-series deformation information in low-coherence areas (such as areas with dramatic terrain fluctuations). Fig.10 The deformation details of HGDS-InSAR in the area with concentrated landslide activity are shown, and its deformation pattern is more refined and the noise level is lower, which further verifies the reliability and accuracy advantages of this method in landslide monitoring in complex mountainous areas.

[0199] Referring to Table 1, it can be observed that, compared with the other three methods, the HGDS-InSAR method proposed in the present invention can obtain significantly higher number of measurement points and spatial density in complex mountainous environments, which fully proves that the HGDS-InSAR method can achieve excellent performance in the high-density monitoring point distribution in complex mountainous areas, and provide more accurate deformation information. Compared with the traditional SBAS-InSAR method, HGDS-InSAR not only increases the number of measurement points by about 86.03%, but also significantly improves the uniformity and density of spatial distribution, which provides more powerful data support for large-scale landslide deformation monitoring. In addition, the advantages of the method of the present invention in complex mountainous areas are not only reflected in the improvement of the density of monitoring points, but also can more accurately capture the subtle deformation of the mountain surface.

[0200] Table 1 The number of measurement points and their spatial density determined by different methods

[0201] method Number of measurement points (MPs) <![CDATA[Spatial density (number / km 2 )]]> PS 304103 36.17 DS 401647 47.77 SBAS 758860 90.26 HGDS 1411761 167.92

[0202] The embodiment of the present invention can divide the functional modules according to the above method example. For example, each functional module can be divided according to each function, or two or more functions can be integrated into one processing module. The above integrated module can be implemented in the form of hardware or in the form of software functional modules. It should be noted that the division of modules in the embodiment of the present invention is schematic and is only a logical function division. There may be other division methods in actual implementation.

[0203] In the case of dividing each functional module into corresponding functional modules, Fig.11 The figure shows a schematic diagram of the structure of a large-scale landslide InSAR deformation rate monitoring system in mountainous areas. Fig.11 As shown, the system includes:

[0204] The SAR image acquisition module 111 is used to acquire N SAR images of the target mountain area; N is an integer greater than 1;

[0205] The mixed baseline processing module 112 is used to perform interference processing on the N-scene SAR images by adopting a single main image combination method to obtain a first interference phase map; and perform interference processing on the N-scene SAR images by adopting a short baseline combination method to obtain a second interference phase map;

[0206] A timing optimization module 113, configured to optimize the first interference phase image by using a timing phase optimization algorithm to obtain a first target interference phase image;

[0207] A Gaussian filter algorithm optimization module 114 is used to optimize the second interference phase image by using a Gaussian filter algorithm to obtain a second target interference phase image;

[0208] An interference phase combination module 115, configured to combine the first target interference phase image with the second target interference phase image to obtain a mixed interference phase image;

[0209] The phase processing and surface deformation analysis module 116 is used to perform phase unwrapping processing and time series analysis on the mixed interference phase image to obtain the surface deformation rate of the target mountain area.

[0210] Optionally, the Gaussian filter algorithm optimization module 114 may be specifically used for:

[0211] Using the formula:

[0212] ;

[0213] Performing Gaussian filtering on the real part of the second interference phase image to obtain a first Gaussian filtering result, and performing Gaussian filtering on the imaginary part of the second interference phase image to obtain a second Gaussian filtering result;

[0214] in, is the first Gaussian filtering result or the second Gaussian filtering result, , is the pixel value of the real part or the pixel value of the imaginary part of the second interference phase image input, For the current pixel The horizontal offset of For the current pixel The vertical offset of is the standard deviation of the Gaussian kernel;

[0215] The first Gaussian filtering result and the second Gaussian filtering result are compositely processed to obtain a second target interference phase map.

[0216] Optionally, the mixed baseline processing module 112 may include:

[0217] A single main image combination processing unit is used to register the main image and the auxiliary image in the N SAR images to obtain multiple groups of first interference image pairs; the main image is the N / 2th image in the N SAR images, and the auxiliary image is an image other than the main image;

[0218] Using the formula:

[0219] ;

[0220] Calculating a first interference phase image of the first interference image pair;

[0221] in, is the first interference phase diagram, and are the two images of the first interference image pair, is the terrain phase compensation term of the first interferometric image pair expressed in complex form, is a mathematical constant, is the product of the amplitudes of the two images of the first interference image pair, For Take the complex conjugate.

[0222] The short baseline combination mode processing unit is used to adopt the formula:

[0223] ;

[0224] Intercepting the N SAR images to obtain a target image;

[0225] in, is the minimum time baseline, is the maximum time baseline;

[0226] Combining images at adjacent time points in the target image in pairs to obtain a second interference image pair;

[0227] The second interference image pair is subjected to interference processing to obtain a second interference phase map.

[0228] Optionally, the timing optimization module 113 may include:

[0229] a target homogeneous pixel determination unit, configured to select a target homogeneous pixel in the first interferometric phase map according to the intensity data of the N SAR images;

[0230] A covariance matrix and initial coherence matrix determination unit, used for performing stack processing and matrix conversion processing on the target homogeneous pixels to obtain a covariance matrix and an initial coherence matrix;

[0231] A positive definiteness detection unit, used for performing a positive definiteness detection on the initial coherence matrix to obtain a target coherence matrix that meets the positive definiteness;

[0232] an eigenvalue calculation unit, used to respectively calculate a first eigenvalue corresponding to the covariance matrix and a second eigenvalue corresponding to the target coherence matrix;

[0233] The first phase optimization unit is used to determine an eigenvector corresponding to the minimum eigenvalue of the first eigenvalue and the second eigenvalue as a first target interference phase map.

[0234] Optionally, the target homogeneous pixel determination unit may be specifically used for:

[0235] Using the formula:

[0236] ;

[0237] Determine the first confidence interval estimate;

[0238] in, is the time series average intensity of the pixel to be detected, is the time series average intensity of the reference pixel; intensity ratio The degree of freedom is F distribution; is the F distribution at the significance level The upper quantile of the lower is the F distribution at the significance level the lower quantile of the lower;

[0239] Confirming the pixels in the first window in the first confidence interval estimate as an initial homogeneous pixel set;

[0240] Calculating an average value of pixel intensities in the initial set of homogeneous pixels, and determining the average value as an estimated value of the reference pixel;

[0241] Using the formula:

[0242] ;

[0243] determining a second confidence interval estimate;

[0244] in, For the standard gamma distribution at a significant level The upper quantile of the lower For the standard gamma distribution at a significant level The lower quantile point, is the estimated value of the reference pixel;

[0245] The pixels in the second confidence interval estimate are determined as target homogeneous pixels.

[0246] Optionally, the covariance matrix and initial coherence matrix determination unit may be specifically used for:

[0247] Using the formula:

[0248]

[0249] Loading the first interference phase image into a three-dimensional phase stack to obtain a three-dimensional phase image;

[0250] in, is the pixel position in the first interference phase image of scene N The phase of It is a three-dimensional phase stack;

[0251] extracting a homogeneous pixel stack within a third window at each pixel position in the three-dimensional phase image;

[0252] reshaping the homogeneous pixel stack into a homogeneous pixel matrix;

[0253] Using the formula:

[0254] ;

[0255] Converting the homogeneous pixel matrix into a covariance matrix;

[0256] in, is the covariance matrix, is a homogeneous pixel stack, is the number of target homogeneous pixel stacks in the third window;

[0257] The absolute values ​​in the covariance matrix are determined as an initial coherence matrix.

[0258] Optionally, the positive detection unit may be specifically used for:

[0259] By performing Cholesky decomposition on the initial coherence matrix, determining whether the initial coherence matrix meets the positive definiteness;

[0260] When the initial coherence matrix does not meet the positive definiteness, adding a target value on the diagonal of the initial coherence matrix to obtain a new coherence matrix;

[0261] It is determined whether the new coherence matrix meets the positive definiteness. If the new coherence matrix does not meet the positive definiteness, the target value is increased until the new coherence matrix meets the positive definiteness, thereby obtaining the target coherence matrix.

[0262] Optionally, the phase processing and surface deformation analysis module 116 may include:

[0263] An atmospheric correction unit, used for performing atmospheric correction processing on the mixed interference phase image to obtain a corrected interference phase image;

[0264] Phase unwrapping unit, used to adopt the formula:

[0265] ;

[0266] Performing phase unwrapping processing on the corrected interference phase image to obtain an unwrapped phase;

[0267] in, is the unwrapped phase, To correct the interferometric phase pattern, Indicates the number of phase jumps;

[0268] The unwrapping phase fitting unit is used to fit the unwrapping phase using a linear regression model to obtain a fitting interference phase function:

[0269] The least square method solving unit is used to solve the fitting interference phase function by using the least square method to obtain the deformation rate and the initial phase.

[0270] The above mainly introduces the solution provided by the embodiment of the present invention from the perspective of the interaction between various modules. It can be understood that in order to realize the above functions, it includes hardware structures and / or software modules corresponding to the execution of various functions. It should be easy for those skilled in the art to realize that, in combination with the units and algorithm steps of each example described in the embodiments disclosed herein, the present invention can be implemented in the form of hardware or a combination of hardware and computer software. Whether a function is executed in the form of hardware or computer software driving hardware depends on the specific application and design constraints of the technical solution. Professional and technical personnel can use different methods to implement the described functions for each specific application, but such implementation should not be considered to exceed the scope of the present invention.

[0271] In the above embodiments, all or part of the embodiments may be implemented by software, hardware, firmware or any combination thereof. When implemented by software, all or part of the embodiments may be implemented in the form of a computer program product. The computer program product includes one or more computer programs or instructions. When the computer program or instruction is loaded and executed on a computer, the process or function described in the embodiment of the present invention is executed in whole or in part. The computer may be a general-purpose computer, a special-purpose computer, a computer network, a terminal, a user device or other programmable device. The computer program or instruction may be stored in a computer-readable storage medium or transmitted from one computer-readable storage medium to another computer-readable storage medium. For example, the computer program or instruction may be transmitted from one website, computer, server or data center to another website, computer, server or data center by wired or wireless means. The computer-readable storage medium may be any available medium that a computer can access or a data storage device such as a server or data center that integrates one or more available media. The available medium may be a magnetic medium, such as a floppy disk, a hard disk, or a tape; it may also be an optical medium, such as a digital video disc (DVD); it may also be a semiconductor medium, such as a solid state drive (SSD).

[0272] Although the present invention is described herein in conjunction with various embodiments, in the process of implementing the claimed invention, those skilled in the art may understand and implement other variations of the disclosed embodiments by viewing the drawings, the disclosure, and the appended claims. In the claims, the word "comprising" does not exclude other components or steps, and "one" or "an" does not exclude multiple situations. A single processor or other unit may implement several functions listed in a claim. Certain measures are recorded in different dependent claims, but this does not mean that these measures cannot be combined to produce good results.

[0273] Although the present invention has been described in conjunction with specific features and embodiments thereof, it is apparent that various modifications and combinations may be made thereto without departing from the spirit and scope of the present invention. Accordingly, this specification and the accompanying drawings are merely exemplary illustrations of the present invention as defined by the appended claims and are deemed to cover any and all modifications, variations, combinations or equivalents within the scope of the present invention. Obviously, those skilled in the art may make various modifications and variations to the present invention without departing from the spirit and scope of the present invention. Thus, the present invention is intended to include such modifications and variations if they fall within the scope of the claims of the present invention and their equivalents.

Claims

1. A method for monitoring deformation rate of large-scale landslides in mountainous areas by InSAR, characterized in that: include: Obtain N SAR images of the target mountain area; N is an integer greater than 1; The N-scene SAR images are interferometrically processed by a single main image combination method to obtain a first interferometric phase map; the N-scene SAR images are interferometrically processed by a short baseline combination method to obtain a second interferometric phase map; The first interference phase image is optimized by using a sequential phase optimization algorithm to obtain a first target interference phase image; Using a Gaussian filtering algorithm to optimize the second interference phase image to obtain a second target interference phase image; Combining the first target interference phase image with the second target interference phase image to obtain a mixed interference phase image; The mixed interference phase image is subjected to phase unwrapping processing and time series analysis to obtain the surface deformation rate of the target mountain area.

2. The InSAR deformation rate monitoring method for large-scale landslides in mountainous areas according to claim 1 is characterized in that: The adopting a Gaussian filtering algorithm to optimize the second interference phase image to obtain a second target interference phase image comprises: Using the formula: ； Performing Gaussian filtering on the real part of the second interference phase image to obtain a first Gaussian filtering result, and performing Gaussian filtering on the imaginary part of the second interference phase image to obtain a second Gaussian filtering result; in, is the first Gaussian filtering result or the second Gaussian filtering result, , is the pixel value of the real part or the pixel value of the imaginary part of the second interference phase image input, For the current pixel The horizontal offset of For the current pixel The vertical offset of is the standard deviation of the Gaussian kernel; The first Gaussian filtering result and the second Gaussian filtering result are compositely processed to obtain a second target interference phase map.

3. The InSAR deformation rate monitoring method for large-scale landslides in mountainous areas according to claim 1 is characterized in that: The interferometric processing of the N SAR images by using a single main image combination method to obtain a first interferometric phase map comprises: The main image and the auxiliary image in the N SAR images are registered to obtain a plurality of first interference image pairs; the main image is the N / 2th image in the N SAR images, and the auxiliary image is an image other than the main image; Using the formula: ； Calculating a first interference phase image of the first interference image pair; in, is the first interference phase diagram, and are the two images of the first interference image pair, is the terrain phase compensation term of the first interferometric image pair expressed in complex form, is a mathematical constant, is the product of the amplitudes of the two images of the first interference image pair, For Take the complex conjugate.

4. The InSAR deformation rate monitoring method for large-scale landslides in mountainous areas according to claim 1 is characterized in that: The interferometric processing of the N SAR images by using a short baseline combination method to obtain a second interferometric phase map comprises: Using the formula: ； Intercepting the N SAR images to obtain a target image; in, is the minimum time baseline, is the maximum time baseline; Combining images at adjacent time points in the target image in pairs to obtain a second interference image pair; The second interference image pair is subjected to interference processing to obtain a second interference phase map.

5. The InSAR deformation rate monitoring method for large-scale landslides in mountainous areas according to claim 1 is characterized in that: The adopting a sequential phase optimization algorithm to optimize the first interference phase image to obtain a first target interference phase image comprises: Selecting target homogeneous pixels in the first interferometric phase map according to the intensity data of the N SAR images; Performing stack processing and matrix conversion processing on the target homogeneous pixels to obtain a covariance matrix and an initial coherence matrix; Performing a positive definiteness test on the initial coherence matrix to obtain a target coherence matrix that meets the positive definiteness requirement; Respectively calculating a first eigenvalue corresponding to the covariance matrix and a second eigenvalue corresponding to the target coherence matrix; An eigenvector corresponding to the minimum eigenvalue of the first eigenvalue and the second eigenvalue is determined as a first target interference phase map.

6. The InSAR deformation rate monitoring method for large-scale landslides in mountainous areas according to claim 5 is characterized in that: The selecting target homogeneous pixels in the first interferometric phase map according to the intensity data of the N SAR images comprises: Using the formula: ； Determine the first confidence interval estimate; in, is the time series average intensity of the pixel to be detected, is the time series average intensity of the reference pixel; intensity ratio The degree of freedom is F distribution; is the F distribution at the significance level The upper quantile of the lower is the F distribution at the significance level the lower quantile of the lower; Confirming the pixels in the first window in the first confidence interval estimate as an initial homogeneous pixel set; Calculating an average value of pixel intensities in the initial set of homogeneous pixels, and determining the average value as an estimated value of the reference pixel; Using the formula: ； determining a second confidence interval estimate; in, For the standard gamma distribution at a significant level The upper quantile of the lower For the standard gamma distribution at a significant level The lower quantile point, is the estimated value of the reference pixel; The pixels within the second window in the second confidence interval estimate are determined as target homogeneous pixels.

7. The InSAR deformation rate monitoring method for large-scale landslides in mountainous areas according to claim 5 is characterized in that: The stacking and matrix conversion processing of the target homogeneous pixels to obtain the covariance matrix and the initial coherence matrix includes: Using the formula: ； Loading the first interference phase image into a three-dimensional phase stack to obtain a three-dimensional phase image; in, is the pixel position in the first interference phase image of scene N The phase of It is a three-dimensional phase stack; extracting a homogeneous pixel stack within a third window at each pixel position in the three-dimensional phase image; reshaping the homogeneous pixel stack into a homogeneous pixel matrix; Using the formula: ； Converting the homogeneous pixel matrix into a covariance matrix; in, is the covariance matrix, is a homogeneous pixel stack, is the number of target homogeneous pixel stacks in the third window; The absolute values ​​in the covariance matrix are determined as an initial coherence matrix.

8. The InSAR deformation rate monitoring method for large-scale landslides in mountainous areas according to claim 5 is characterized in that: The performing positive definiteness detection on the initial coherence matrix to obtain a target coherence matrix that meets the positive definiteness includes: By performing Cholesky decomposition on the initial coherence matrix, determining whether the initial coherence matrix meets the positive definiteness; When the initial coherence matrix does not meet the positive definiteness, adding a target value on the diagonal of the initial coherence matrix to obtain a new coherence matrix; It is determined whether the new coherence matrix meets the positive definiteness. If the new coherence matrix does not meet the positive definiteness, the target value is increased until the new coherence matrix meets the positive definiteness, thereby obtaining the target coherence matrix.

9. The InSAR deformation rate monitoring method for large-scale landslides in mountainous areas according to claim 1 is characterized by: The performing phase unwrapping processing and time series analysis on the hybrid interference phase image to obtain the surface deformation rate of the target mountain area includes: Performing atmospheric correction processing on the mixed interference phase image to obtain a corrected interference phase image; Using the formula: ； Performing phase unwrapping processing on the corrected interference phase image to obtain an unwrapped phase; in, is the unwrapped phase, To correct the interferometric phase pattern, Indicates the number of phase jumps; The unwrapped phase is fitted using a linear regression model to obtain a fitted interference phase function: The fitted interference phase function is solved by the least square method to obtain the deformation rate and the initial phase.

10. A large-scale landslide InSAR deformation rate monitoring system in mountainous areas, characterized in that: include: The SAR image acquisition module is used to acquire N SAR images of the target mountain area; N is an integer greater than 1; A hybrid baseline processing module is used to perform interference processing on the N-scene SAR images by adopting a single main image combination method to obtain a first interference phase map; and to perform interference processing on the N-scene SAR images by adopting a short baseline combination method to obtain a second interference phase map; A timing optimization module, used to optimize the first interference phase image by using a timing phase optimization algorithm to obtain a first target interference phase image; A Gaussian filter algorithm optimization module, used to optimize the second interference phase image by using a Gaussian filter algorithm to obtain a second target interference phase image; An interference phase combination module, used for combining the first target interference phase image with the second target interference phase image to obtain a mixed interference phase image; The phase processing and surface deformation analysis module is used to perform phase unwrapping processing and timing analysis on the mixed interference phase image to obtain the surface deformation rate of the target mountain area.

Images