Landslide deformation analysis method based on Offset-Tracking technology

Abstract

The present invention provides a landslide deformation analysis method based on Offset-Tracking technology, which forms multiple groups of image pairs by acquiring time series images covering the same landslide monitoring area; performing matching search on each group of image pairs, determining the integer pixel point pair with the largest mutual correlation coefficient, the sub-pixel point pair with the largest mutual correlation coefficient around the integer pixel point pair, and the target sub-pixel point pair around the sub-pixel point pair; then obtaining the deformation observation matrix of each group of image pairs; establishing the observation equations of the deformation observation value matrix, deformation time matrix and deformation rate matrix of the image pair, using the least squares criterion or the weighted singular value decomposition method to solve according to different situations, and obtaining the rates of two images in each group of image pairs; finally generating a two-dimensional deformation time series image. The present invention can effectively form reliable two-dimensional deformation time series images on different software according to the landslide deformation offset, providing researchers with an effective analysis basis.

Description

Technical Field

[0001] The invention belongs to the technical field of landslide monitoring, and in particular relates to a landslide deformation analysis method based on Offset-Tracking technology. Background Art

[0002] Landslide is a very serious natural disaster that can cause great property losses to the country and the people. Due to its serious harm to society, landslide monitoring has become an important means and method in landslide research.

[0003] Among the existing landslide monitoring technologies, GPS monitoring is favored by researchers because of its advantages such as all-weather, high-precision and fully automatic. However, the measurement data of GPS monitoring is complex to solve, the monitoring signal is easily blocked, and it is expensive and has high measurement costs, making it unsuitable for large-scale monitoring. Another monitoring method is drone photogrammetry, which is easily affected by climatic conditions, has poor continuity in all-weather operation, and cannot obtain long-term surface deformation processes. There is also an offset tracking technology based on synthetic aperture radar (SAR) that can measure surface displacement with high precision, but this method requires the interference process to be obtained through interference phase, which is prone to produce decoherence results, resulting in reduced measurement accuracy and failure of phase unwrapping.

[0004] With the development of landslide monitoring technology, researchers have begun to introduce optical imaging to identify landslides. In optical image recognition, most of the visual interpretation is combined with the landslide deformation area obtained by other monitoring methods mentioned above for comparative analysis to determine the accurate landslide deformation area. However, in actual practice, landslides are a process with precursors and continuous development, and this method cannot provide a reliable analysis process. Summary of the invention

[0005] In order to solve the above problems existing in the prior art, the present invention provides a landslide deformation analysis method based on Offset-Tracking technology. The technical problem to be solved by the present invention is achieved through the following technical solutions:

[0006] The present invention provides a landslide deformation analysis method based on Offset-Tracking technology, which includes:

[0007] Obtain time series images covering the same landslide monitoring area;

[0008] According to the similarity of the solar altitude angle and the azimuth angle, the time series images are combined into image pairs to obtain multiple groups of image pairs;

[0009] Perform a matching search on each image pair to determine the integer pixel pair with the largest mutual correlation coefficient, the sub-pixel pair with the largest mutual correlation coefficient around the integer pixel pair, and the target sub-pixel pair around the sub-pixel pair;

[0010] The integer displacement between two pixels in an integer pixel pair, the fractional displacement between two sub-pixels in a sub-pixel pair, and the fractional offset between two target sub-pixels in a target sub-pixel pair are accumulated to obtain the deformation measurement matrix of each image pair;

[0011] Establish observation equations of deformation observation value matrix, deformation time matrix and deformation rate matrix of image pairs;

[0012] According to the elements of the deformation time matrix, the observation equation is solved using the least squares criterion or the weighted singular value decomposition method to obtain the rates of the two images in each image pair;

[0013] Based on the deformation rates of the two images in each image pair, a two-dimensional deformation time series image is generated.

[0014] Optionally, the time series images are multiple images sorted in time series; the number of image pairs and the number of images sorted in time series satisfy:

[0015]

[0016] Where N represents the number of images sorted in time series, and M represents the number of image pairs.

[0017] Optionally, performing a matching search on each group of image pairs to determine the integer pixel pair with the largest mutual correlation coefficient, the sub-pixel pair with the largest mutual correlation coefficient around the integer pixel pair, and the target sub-pixel pair around the sub-pixel pair includes:

[0018] One image in each image pair is used as the reference image and the other image is used as the search image. An integer pixel search is performed on each image pair to determine the integer pixel pair at the location where the maximum value of the mutual correlation coefficient is located.

[0019] Perform interpolation calculation on the sub-pixel points around the whole pixel point pair to determine the sub-pixel point pair at the position where the maximum value of the mutual correlation coefficient is located;

[0020] Gaussian regression analysis is performed on the sub-pixel points around the sub-pixel point pair to determine the target sub-pixel point pair.

[0021] Optionally, one image in each image pair is used as a reference image and the other image is used as a search image. An integer pixel search is performed on each image pair, and the integer pixel pair at the location where the maximum value of the mutual correlation coefficient is determined includes:

[0022] One image in each image pair is used as the reference image and the other as the search image;

[0023] Use a first sliding window and a second sliding window of the same size to slide in the target area of ​​the reference image and the search area of ​​the search image respectively, and calculate the first mutual correlation coefficient in turn to obtain a first mutual correlation coefficient matrix;

[0024] Determine the integer pixel point pair at the position where the maximum value of the first mutual correlation coefficient is located;

[0025] Interpolation calculation is performed on sub-pixel points around the whole pixel point pair to determine the sub-pixel point pair at the position where the maximum value of the mutual correlation coefficient is located, including:

[0026] Meshing is performed at the location where the maximum value of the first mutual correlation coefficient is located;

[0027] Using two third sliding windows and a fourth sliding window of the same size to slide in the grid, and using a bilinear interpolation algorithm to sequentially calculate the second mutual correlation coefficients of the third sliding window and the fourth sliding window at the gridded position, to obtain a second mutual correlation coefficient matrix;

[0028] Gaussian regression analysis is performed on the sub-pixel points around the sub-pixel point pair to determine that the target sub-pixel point pair includes:

[0029] Gaussian regression analysis is performed on the sub-pixel point at the position of the second cross-correlation coefficient peak and the 9 surrounding sub-pixels to obtain the optimal target sub-pixel point of the Gaussian regression model.

[0030] Optionally, according to the elements of the deformation time matrix, the observation equation is solved using the least squares criterion or the weighted singular value decomposition method to obtain the rates of the two images in each image pair, including:

[0031] When the deformation time matrix is ​​full rank, the observation equation is solved using the least squares criterion to obtain the rates of the two images in each image pair;

[0032] When the deformation time matrix is ​​rank-deficient, the observation equation is solved using the weighted singular value decomposition (SVD) method to obtain the rates of the two images in each image pair.

[0033] Optionally, when the deformation time matrix is ​​rank-deficient, the observation equation is solved using the weighted singular value decomposition (SVD) method to obtain the rates of the two images in each image pair including:

[0034] When the deformation time matrix is ​​rank-deficient, the observation equation is transformed using a weighted matrix to obtain the transformed observation equation;

[0035] Perform singular value decomposition on the weight matrix and deformation time matrix in the transformed observation equation;

[0036] Based on the singular value decomposition results of the weighted matrix and the deformation matrix, the transformed observation equation is solved to obtain the rates of the two images in each image pair.

[0037] Optionally, the deformation observation matrix includes deformation observation values ​​of all image pairs, the deformation time matrix includes the time interval between two images in each image pair, and the deformation rate matrix includes: the rates of two images in each image pair;

[0038] The observation equation is expressed as:

[0039] Bv=d

[0040] Wherein, B is an M×N coefficient matrix, and the elements in the matrix B consist of the time intervals between adjacent images in the image pair; m i represents the main image in the image pair, and the slave image is s i , the elements in matrix B can be expressed as:

[0041]

[0042] d is an M×1 deformation measurement matrix, the matrix elements are composed of the deformation observation values ​​of all cross-correlated image pairs, d = [d1d2 … d M ] T ; v is an N×1 matrix whose elements consist of the rates between adjacent images in an image pair, v = [v1y2… v n ] T .

[0043] Optionally, the observation equation is solved using the least squares criterion, and the rates of the two images in each image pair are expressed as:

[0044] v=(B T PB) -1 B T Pd

[0045] Where P is an M×M weighted matrix, P=diag(ω1,ω2,…,ω M ),ω i The standard deviation σ of each image pair i The reciprocal of is an indicator to measure the weight of each image pair;

[0046] Optionally, the transformed observation equation is expressed as:

[0047]

[0048] The results of singular value decomposition of the weight matrix and coefficient matrix are expressed as:

[0049]

[0050] Among them, svd is the singular value decomposition function; U is an M×M orthogonal matrix, which is the ZZT eigenvector; V is an N×M matrix, which is the Z T The eigenvector of Z; S is an M×M diagonal matrix;

[0051] The transformed observation equation is solved using the singular value decomposition result, and the rate of the two images in each image pair is expressed as:

[0052]

[0053] Among them, S + is the pseudo-inverse of S, given by the singular value σ i Composition, S + =diag(1 / σ1,1 / σ2,...,1 / σ N-L+1 , 0, 0, …, 0), the number of singular values ​​that are 0 is (L-1).

[0054] Optionally, after obtaining the deformation observation matrix of each set of image pairs, the landslide deformation analysis method further includes:

[0055] Decompose the deformation observation matrix in the east-west direction and the north-south direction to obtain the deformation observation values ​​in the east-west direction and the north-south direction;

[0056] Calculate the signal-to-noise ratio of the image pairs used to obtain the deformation observations in the east-west and north-south directions;

[0057] Use signal-to-noise ratio to verify the deformation measurement matrix;

[0058] Filtering the deformation observation values ​​in the deformation observation matrix to obtain smooth deformation observation values;

[0059] The smoothed deformation observations are decomposed into east-west and north-south directions to obtain the deformation observations in the east-west and north-south directions.

[0060] The present invention provides a landslide deformation analysis method based on Offset-Tracking technology, which forms multiple groups of image pairs by acquiring time series images covering the same landslide monitoring area; performing matching search on each group of image pairs, determining the integer pixel point pair with the largest mutual correlation coefficient, the sub-pixel point pair with the largest mutual correlation coefficient around the integer pixel point pair, and the target sub-pixel point pair around the sub-pixel point pair; then obtaining the deformation observation matrix of each group of image pairs; establishing the observation equations of the deformation observation value matrix, deformation time matrix and deformation rate matrix of the image pair, using the least squares criterion or the weighted singular value decomposition method to solve according to different situations, and obtaining the rates of two images in each group of image pairs; and finally generating a two-dimensional deformation time series image. The present invention can effectively form reliable two-dimensional deformation time series images on different software according to the landslide deformation offset, providing researchers with an effective analysis basis.

[0061] The present invention will be further described in detail below with reference to the accompanying drawings and embodiments. BRIEF DESCRIPTION OF THE DRAWINGS

[0062] Figure 1 It is a flow chart of a landslide deformation analysis method based on Offset-Tracking technology provided by an embodiment of the present invention;

[0063] Figure 2 is a schematic diagram of a process for generating a two-dimensional temporal deformation image provided by an embodiment of the present invention;

[0064] Figure 3 is a schematic diagram of the coverage of the study area provided by an embodiment of the present invention;

[0065] Figure 4 is a cumulative time series surface deformation result map provided by an embodiment of the present invention;

[0066] Figure 5 It is a two-dimensional deformation result diagram of the 10·10 landslide provided by an embodiment of the present invention;

[0067] Figure 6 This is a result diagram of the east-west deformation of the 11.3 landslide provided by an embodiment of the present invention;

[0068] Figure 7 It is a two-dimensional deformation result diagram after the Baige landslide provided by an embodiment of the present invention;

[0069] Figure 8 is a rear profile analysis diagram of the Baige landslide provided by an embodiment of the present invention;

[0070] Fig. 9 is a schematic diagram of the time series offset of the COSI-Corr Baige landslide provided by an embodiment of the present invention;

[0071] Fig.10It is a schematic diagram of the time series offset of the CARST Baige landslide provided by an embodiment of the present invention. DETAILED DESCRIPTION

[0072] The present invention is further described in detail below with reference to specific embodiments, but the embodiments of the present invention are not limited thereto.

[0073] Existing technologies mainly analyze and study the evolution of landslides through field monitoring and remote sensing monitoring. Since the field monitoring process requires a lot of manpower and material resources and has a high risk index, it is generally used for auxiliary monitoring. Remote sensing monitoring is divided into optical remote sensing monitoring and radar remote sensing monitoring according to the source of the image. Radar remote sensing monitoring obtains the deformation process by studying the phase and amplitude information in the image; optical remote sensing monitoring was originally based on visual interpretation. With the increasing maturity of satellite launch technology and aerial measurement technology, the accuracy of optical image acquisition has been continuously improved. The acquisition of multi-phase images through image cross-correlation algorithms to monitor the movement of landslides has been widely cited.

[0074] like Figure 1 As shown, the present invention provides a landslide deformation analysis method based on Offset-Tracking technology, which includes:

[0075] S1, obtain time series images covering the same landslide monitoring area;

[0076] Among them, the time series images are multiple images sorted in time series;

[0077] S2, according to the similarity of the solar altitude angle and the azimuth angle, the time series images are combined into image pairs to obtain multiple groups of image pairs;

[0078] It is worth noting that optical image time series deformation monitoring is evolved from short baseline differential interferometry (SBAS-InSAR). N+1 optical images covering the same area are arranged in order of time acquisition (t_0,…,t_N), and image pairs are formed according to the similarity of solar altitude angle and solar azimuth angle to obtain M cross-correlated image pairs. The relationship between the number of image pairs M and the number of images (N+1) satisfies:

[0079]

[0080] Where N represents the number of images sorted in time series, and M represents the number of image pairs.

[0081] S3, performing a matching search on each group of image pairs to determine the integer pixel pair with the largest mutual correlation coefficient, the sub-pixel pair with the largest mutual correlation coefficient around the integer pixel pair, and the target sub-pixel pair around the sub-pixel pair;

[0082] As an optional implementation manner of the present invention, S3 includes:

[0083] S31, using one image in each image pair as a reference image and the other image as a search image, performing integer pixel search on each image pair to determine the integer pixel pair at the location where the maximum value of the mutual correlation coefficient is located;

[0084] S32, performing interpolation calculation on sub-pixel points around the integer pixel point pair to determine the sub-pixel point pair at the position where the maximum value of the mutual correlation coefficient is located;

[0085] S33, performing Gaussian regression analysis on sub-pixel points around the sub-pixel point pair to determine the target sub-pixel point pair.

[0086] As an optional implementation manner of the present invention, S31 includes:

[0087] S311, using one image in each image pair as a reference image and the other image as a search image;

[0088] S312, using two first sliding windows and second sliding windows of the same size to slide in the target area of ​​the reference image and the search area of ​​the search image respectively, and calculating first mutual correlation coefficients in sequence to obtain a first mutual correlation coefficient matrix;

[0089] S313, determining the integer pixel point at which the maximum value of the first mutual correlation coefficient is located;

[0090] As an optional implementation of the present invention, step S32 includes:

[0091] S321, gridding is performed at the position where the maximum value of the first mutual correlation coefficient is located;

[0092] S322, using two third sliding windows and a fourth sliding window of the same size to slide in the grid, and using a bilinear interpolation algorithm to sequentially calculate second cross-correlation coefficients of the third sliding window and the fourth sliding window at the gridded position, to obtain a second cross-correlation coefficient matrix;

[0093] S323, performing Gaussian regression analysis on sub-pixel points around the sub-pixel point pair to determine that the target sub-pixel point pair includes:

[0094] S324, performing Gaussian regression analysis on the sub-pixel point at the position where the second cross-correlation coefficient peak is located and the 9 surrounding sub-pixel points to obtain the optimal target sub-pixel point of the Gaussian regression model.

[0095] S4, accumulating the integer displacement between two pixels in the integer pixel pair, the decimal displacement between two sub-pixels in the sub-pixel pair, and the decimal offset between two target sub-pixels in the target sub-pixel pair, to obtain a deformation measurement matrix for each image pair;

[0096] Among them, the deformation observation matrix includes the deformation observation values ​​of all image pairs, the deformation time matrix includes the time interval between two images in each image pair, and the deformation rate matrix includes: the rate of two images in each image pair;

[0097] It is worth noting that the present invention uses the Offset-Tracking technology to solve the deformation observation value, that is, the offset. The basic principle of the Offset-Tracking technology is to calculate the offset between pixels by using a correlation matching algorithm for optical images. The core of the technology is the correlation matching algorithm. The present invention mainly introduces the implementation of the frequency domain-based cross-correlation algorithm. The frequency domain cross-correlation algorithm is based on the principle of Fourier transform, which converts the offset value in the spatial domain into a phase difference value in the frequency domain, thereby obtaining the displacement difference between the two images. The calculation formula of the frequency domain cross-correlation algorithm is:

[0098]

[0099]

[0100] Wherein F1(u, v) and F2(u, v) are Fourier transform functions of the two original images f1(x, y) and f2(x, y), * represents the complex conjugate, and (x0, y0) represents the relative displacement value of the two original images.

[0101] The image offset value (x0, y0) can be solved by formula (2, 3) in two ways. One is to solve it in the frequency domain by fitting the phase plane, and the other is to obtain the peak value of the Dirac function δ(x-x0, y-y0) by inverse Fourier transform, and then use the surface fitting algorithm to solve the offset value (x0, y0), that is, the deformation observation value.

[0102] In practice, the steps of obtaining the two-dimensional offset of the image using the Offset-Tracking technology are mainly divided into the following steps:

[0103] (1) Obtain two images before and after the landslide, and use them as the reference image and search image, respectively. Use two sliding windows s1 and s2 of the same size to slide in the target area r1 of the reference image f1 and the search area r2 of the search image f2, respectively, and calculate the mutual correlation coefficient in turn to obtain the mutual correlation coefficient matrix. Calculate the integer displacement d1 between the two images from the position where the maximum value of the mutual correlation coefficient is located.

[0104] (2) Gridding is performed using surrounding pixels at the location of the maximum value of the cross-correlation coefficient. Two other sliding windows s3 and s4 of the same size are selected, and the cross-correlation coefficients of s3 and s4 at the gridding location are calculated in turn according to the bilinear interpolation algorithm to obtain another cross-correlation coefficient matrix. After Gaussian filtering, the decimal offset d2 is calculated based on the location of the maximum value of the cross-correlation coefficient.

[0105] (3) Gaussian regression analysis is performed using the peak value of the mutual correlation coefficient obtained in step 2) and its surrounding points 9 to obtain a decimal offset d3 with an accuracy less than d2.

[0106] (4) The total displacement d is obtained from d1, d2, and d3.

[0107] It is worth noting that the total displacement d is the offset between pixels. For the image pair, an offset matrix, namely the deformation observation matrix, is finally formed, and each element of the matrix is ​​the value of the offset.

[0108] refer to Figure 2 After obtaining the deformation observation matrix of each group of image pairs, the landslide deformation analysis method of the present invention further includes:

[0109] Step a: Decompose the deformation observation matrix in the east-west direction and the north-south direction to obtain the deformation observation values ​​in the east-west direction and the north-south direction;

[0110] Step b: Calculate the signal-to-noise ratio of the image pairs used to obtain the deformation observation values ​​in the east-west direction and the north-south direction;

[0111] Step c: Use the signal-to-noise ratio to verify the deformation measurement matrix;

[0112] Step d: filtering the deformation observation values ​​in the deformation observation matrix to obtain smooth deformation observation values;

[0113] Step e: Decompose the smoothed deformation observations into east-west and north-south directions to obtain deformation observations in the east-west and north-south directions.

[0114] S5, establishing observation equations of deformation observation value matrix, deformation time matrix and deformation rate matrix of the image pair;

[0115] The deformation rate between adjacent image pairs is regarded as an unknown number, and the deformation time and deformation observation value of all image pairs are regarded as known numbers, and the corresponding observation equation is constructed.

[0116] The observation equation is expressed as:

[0117] Bv=d (4)

[0118] Wherein, B is an M×N coefficient matrix, and the elements in the matrix B consist of the time intervals between adjacent images in the image pair; m i represents the main image in the image pair, and the slave image is s i , the elements in matrix B can be expressed as:

[0119]

[0120] d is an M×1 deformation measurement matrix, the matrix elements are composed of the deformation observation values ​​of all cross-correlated image pairs, d = [d1d2 … d M ] T ; v is an N×1 matrix whose elements consist of the rates between adjacent images in an image pair, v = [v1v2…v n ] T .

[0121] S6, according to the elements of the deformation time matrix, the observation equation is solved using the least squares criterion or the weighted singular value decomposition method to obtain the rates of the two images in each image pair;

[0122] The present invention solves the unknown number v, where v is an N×1 matrix, and the matrix elements are composed of the rates between adjacent images in each image pair, that is, v=[v1v2 … v n ] T .

[0123] S7, generating a two-dimensional deformation time-series image based on the deformation rates of the two images in each image pair.

[0124] As an optional implementation manner of the present invention, S6 includes:

[0125] S61, when the deformation time matrix is ​​full rank, the observation equation is solved using the least squares criterion to obtain the rates of the two images in each image pair;

[0126] The solution process is divided into two cases. One is when B is full rank, and the least squares criterion is used to solve the observation equation. The rate of the two images in each image pair is expressed as:

[0127] v=(B T PB) -1 B T Pd (6)

[0128] Where P is an M×M weighted matrix, P=diag(ω1,ω2,…,ω M ),ω i The standard deviation σ of each image pair i The reciprocal of is an indicator to measure the weight of each image pair;

[0129]

[0130] S62, when the deformation time matrix is ​​rank-deficient, the observation equation is solved using a weighted singular value decomposition (SVD) method to obtain the rates of the two images in each image pair.

[0131] Another case is when the matrix B is rank deficient and the unknowns are solved by weighted singular value decomposition (SVD) method.

[0132] This step S62 can be implemented by the following steps:

[0133] S621, when the deformation time matrix is ​​rank-deficient, transform the observation equation using a weighted matrix to obtain a transformed observation equation;

[0134] First, the observation equation needs to be transformed. The transformed observation equation is expressed as:

[0135]

[0136] S622, performing singular value decomposition on the weighted matrix and the deformation time matrix in the transformed observation equation;

[0137] Among them, the results of singular value decomposition of the weight matrix and coefficient matrix are expressed as:

[0138]

[0139] Among them, svd is the singular value decomposition function; U is an M×M orthogonal matrix, which is ZZ T Eigenvector; V is an N×M matrix, which is Z T The eigenvector of Z; S is an M×M diagonal matrix;

[0140] S623, based on the singular value decomposition results of the weighted matrix and the deformation matrix, the transformed observation equation is solved to obtain the rates of the two images in each image pair.

[0141] Finally, the unknown v is solved, that is, the transformed observation equation is solved using the singular value decomposition result, and the rate of the two images in each image pair is expressed as:

[0142]

[0143] Among them, S + is the pseudo-inverse of S, given by the singular value σ i Composition, S + =diag(1 / σ1,1 / σ2,…,1 / σ N-L+1 , 0, 0, …, 0), the number of singular values ​​that are 0 is (L-1).

[0144] The present invention provides a landslide deformation analysis method based on Offset-Tracking technology, which forms multiple groups of image pairs by acquiring time series images covering the same landslide monitoring area; performing matching search on each group of image pairs, determining the integer pixel point pair with the largest mutual correlation coefficient, the sub-pixel point pair with the largest mutual correlation coefficient around the integer pixel point pair, and the target sub-pixel point pair around the sub-pixel point pair; then obtaining the deformation observation matrix of each group of image pairs; establishing the observation equations of the deformation observation value matrix, deformation time matrix and deformation rate matrix of the image pair, using the least squares criterion or the weighted singular value decomposition method to solve according to different situations, and obtaining the rates of two images in each group of image pairs; and finally generating a two-dimensional deformation time series image. The present invention can effectively form reliable two-dimensional deformation time series images on different software according to the landslide deformation offset, providing researchers with an effective analysis basis.

[0145] Offset-Tracking technology, also known as offset tracking technology, is a new technology for monitoring surface deformation based on image cross-correlation algorithm. It is widely used in monitoring landslide deformation, co-seismic deformation, glacier flow rate, dune movement and other surface deformation. The present invention applies CARST (Cryosphere And Remote Sensing Toolkit) software to the landslide deformation analysis method of the present invention to monitor landslide surface deformation. The software COSI-Corr (Co-registration of Optically Sensed Images and Correlation) uses offset tracking technology to obtain surface deformation information, providing more comprehensive technical support for landslide surface deformation monitoring and making deformation results more reliable. Both softwares enable optical remote sensing images to be used in Windows and Linux systems to monitor landslide surface deformation through offset tracking technology, solving the limitation of obtaining landslide deformation results in a single software. The present invention uses optical images for offset tracking (Offset-Tracking), obtains the landslide time series deformation process and sliding rate, and analyzes the inducing factors and creep process.

[0146] like Figure 3 As shown in the figure, the Baige landslide is located on the right bank of the Jinsha River in Baige Village, Jiangda County, Qamdo City, Tibet Autonomous Region, with coordinates of 31°04′55.8″N, 98°42′14.69″. It is located in the transition zone between the Qinghai-Tibet Plateau and the Sichuan Basin - the Jinsha River tectonic zone. Due to the large terrain difference, landform erosion and weathering, and abnormal climate and environmental changes, the area experienced two consecutive landslides on October 11, 2018 and November 3, 2018, which were referred to as the "10·10" landslide and the "11·3" landslide respectively.

[0147] The present invention uses 20 Sentinel-2 images acquired on Google Earth Engine before and after the Baige landslide disaster, covering a time span from November 13, 2015 to November 11, 2021. The data details are shown in Table 1. In the present invention, two images acquired on November 13, 2015 and February 5, 2018 before the "10·10" landslide, two images acquired on October 28, 2018 and November 2, 2018 before the "11·3" landslide, and two images acquired on November 12, 2018 and November 11, 2021 after the Baige landslide are selected as the original data for offset tracking and comparison by the COSI-Corr and CARST methods. 16 images covering November 2015 to June 2018 are selected for offset tracking and cumulative displacement extraction of surface deformation time series.

[0148] Table 1 Sentinel-2 imaging parameters

[0149]

[0150] In the present invention, optical images with a month number not greater than 3 and a low cloud content in the white grid area are selected for time series analysis. Image pairs are selected based on the solar altitude angle and solar azimuth angle of the screened images, and each scene image with a similar solar altitude angle and solar azimuth angle is combined into an image pair to obtain 65 groups of image pairs. The offset of each group of image pairs is tracked, and image cross-correlation processing is performed according to the parameters set in Table 2 to obtain the two-dimensional deformation information of each group of image pairs. The east-west displacement and north-south displacement obtained for each pair of image pairs are used as known deformation values, and weights are assigned to image pairs with different time intervals to obtain the cumulative time series surface deformation results. Figure 4 shown.

[0151] In the present invention, the selected 6 images are subjected to image cross-correlation processing by setting step size and window parameters, as shown in Table 2, and the surface deformation results in the east-west and north-south directions are obtained and compared and analyzed.

[0152] Table 2 Detailed information on software parameters

[0153]

[0154] Both methods use the offset-tracking technology to obtain the two-dimensional deformation of the "10·10" Baige landslide from November 2015 to February 2018. Figure 5It can be seen that the two results have similar spatial distribution ranges. In the past three years, the maximum two-dimensional deformation value obtained by COSI-Corr is 31.53m, and the maximum two-dimensional deformation value obtained by CARST is 25.75m. From the known results, it can be seen that the "10·10" Baige landslide slides in an east-west direction as a whole. Therefore, the profile line AA' is selected, and four profile points P1, P2, P3, and P4 are selected along the profile line, such as Figure 5 shown.

[0155] The trend of the two-dimensional cumulative deformation values ​​of the east-west profile line AA' in COSI-Corr and CARST is consistent, and the average two-dimensional deformation along the profile line obtained by the two methods is 27.066m and 21.828m, respectively. This shows that the deformation value obtained by CARST is generally smaller than the deformation result of COSI-Corr, but both methods can reflect the overall deformation range of the "10·10" landslide.

[0156] According to the existing results, the deformation area of ​​the "11·3" Baige landslide event mainly occurred in the west and upper part of the "10·10" landslide event. Therefore, two profile lines BB' and CC' are selected, such as Figure 6 The east-west cumulative displacement values ​​obtained by COSI-Corr and CARST range from -14.95m to 17.63m and -13.75m to 18m, respectively, where positive means moving eastward and negative means moving westward, indicating that the movement rate of the landslide was close to 5 meters per day in the past 5 days. Figure 6 It can be seen that the COSI-Corr results are closer to the actual situation than the CARST results. The deformation area in COSI-Corr is divided into two directions: moving from west to northeast and from south to northeast, while the deformation area in CARST tends to be more overall deformation, indicating that the COSI-Corr results are more reliable.

[0157] The changes in the two-dimensional deformation values ​​of COSI-Corr and CARST along the profiles BB' and CC' show that the central areas of Region 3 and Region 1 are the points with higher deformation values, the transition zone from Region 3 to Region 1 is the lowest point of deformation value, and the transition from Region 2 to Region 1 also shows a trend of first decreasing and then increasing. The average deformation values ​​of the two-dimensional deformation results obtained by CARST along the profiles BB' and CC' are 3.68m and 7.56m respectively, indicating that although the deformation results obtained by CARST fail to show the local deformation results, there is still a certain fluctuation trend on the profile line, and the fluctuation range is large.

[0158] In order to understand whether the Baige landslide continued to deform after it occurred, images from January 2019 to November 2021 were selected for deformation monitoring. The initial search window and the final search window selected by COSI-Corr for cross-correlation matching are both 32*32, with a step size of 1*1. The maximum displacement of the obtained two-dimensional deformation value is about 20m, as shown in Figure 2. Figure 7 As shown in (left), there are three obvious deformation areas; the initial search window selected by CARST during the cross-correlation matching process is 64*64, the final search window is 8*8, and the step size is 1*1. The results are as follows Figure 7 As shown in (right), the maximum displacement value in the two-dimensional deformation result is 20.153m, which is similar to the maximum displacement value obtained by COSI-Corr. The deformation area appears in the figure in the form of sheets, and the spatial distribution range is larger and wider than the specific deformation position in COSI-Corr.

[0159] To further understand the deformation area after the Baige landslide, two profiles DD' and EE' were drawn along the three deformation areas obtained by COSI-Corr. The distribution positions of the two profiles in the results obtained by CARST are shown in the figure below. Figure 7 As shown in the figure, the two-dimensional deformation values ​​along the section lines DD' and EE' change in the two methods, as shown in the figure. Figure 8 As shown, it shows that the deformation value along the profile line DD' in area 1 does not produce obvious deformation by COSI-Corr, while CARST produces larger deformation values; in area 2, the deformation values ​​obtained by the two methods are consistent, and the change trends are also the same, indicating that area 2 is an obvious deformation area; the deformation value along the profile line EE' in area 3 has larger deformation values ​​in both methods, with similar deformation trends and no obvious fluctuations, indicating that area 3 is also a deformation area.

[0160] Taking the "10·10" Baige landslide as an example, the displacement of the Baige landslide and its evolution over time are studied. By processing and analyzing 16 Sentine-2 images from 2015-11 to 2018-06 (Table 1), the singular value decomposition algorithm (SVD) is used to calculate the time series surface deformation of the study area. In the cross-correlation results of the "10·10" Baige landslide, the north-south horizontal displacement value has less influence on the two-dimensional deformation results than the east-west displacement value. Therefore, the total deformation time series offset results of the Baige landslide are shown in this paper. Fig. 9 It shows that when using COSI-Corr, the two-dimensional temporal deformation displacement field of the study area during this period; Fig.10 It shows the results of the 2D temporal deformation migration of the study area when using CARST;

[0161] Since the offsets of the image pairs obtained by COSI-Corr and CARST are different when tracking the offsets, and there will be errors to varying degrees, the image pairs selected by the present invention are not exactly the same in the two methods. The offset results of the image pairs with large errors in the two methods are eliminated. Finally, COSI-Corr uses 65 image pairs for singular value decomposition, and CARST uses 59 image pairs for singular value decomposition. Fig. 9 It can be seen that the maximum cumulative displacement deformation value calculated by COSI-Corr is 43.5m. Fig.10 It can be seen that the maximum cumulative displacement value calculated by CARST is 35m, and the two-dimensional deformation result obtained by CARST is smaller than the two-dimensional deformation result of COSI-Corr.

[0162] The above contents are further detailed descriptions of the present invention in combination with specific preferred embodiments, and it cannot be determined that the specific implementation of the present invention is limited to these descriptions. For ordinary technicians in the technical field to which the present invention belongs, several simple deductions or substitutions can be made without departing from the concept of the present invention, which should be regarded as falling within the protection scope of the present invention.

Claims

1. A landslide deformation analysis method based on Offset-Tracking technology, characterized in that: include: Obtain time series images covering the same landslide monitoring area; According to the similarity of the solar altitude angle and the azimuth angle, the time series images are combined into image pairs to obtain multiple groups of image pairs; Perform a matching search on each image pair to determine the integer pixel pair with the largest mutual correlation coefficient, the sub-pixel pair with the largest mutual correlation coefficient around the integer pixel pair, and the target sub-pixel pair around the sub-pixel pair; The integer displacement between two pixels in an integer pixel pair, the fractional displacement between two sub-pixels in a sub-pixel pair, and the fractional offset between two target sub-pixels in a target sub-pixel pair are accumulated to obtain the deformation measurement matrix of each image pair; Establish observation equations of deformation observation value matrix, deformation time matrix and deformation rate matrix of image pairs; According to the elements of the deformation time matrix, the observation equation is solved by using the least squares criterion or the weighted singular value decomposition method to obtain the rates of the two images in each image pair; Based on the deformation rates of the two images in each group of image pairs, a two-dimensional deformation time-series image is generated.

2. The landslide deformation analysis method according to claim 1, characterized in that: The time series images are multiple images sorted in time series; the number of image pairs and the number of images sorted in time series satisfy: Where N represents the number of images sorted in time series, and M represents the number of image pairs.

3. The landslide deformation analysis method according to claim 1, characterized in that: The matching search is performed on each group of image pairs to determine the integer pixel pair with the largest mutual correlation coefficient, the sub-pixel pair with the largest mutual correlation coefficient around the integer pixel pair, and the target sub-pixel pair around the sub-pixel pair, including: One image in each image pair is used as the reference image and the other image is used as the search image. An integer pixel search is performed on each image pair to determine the integer pixel pair at the location where the maximum value of the mutual correlation coefficient is located. Performing interpolation calculation on sub-pixel points around the integer pixel point pair to determine the sub-pixel point pair at the position where the maximum value of the mutual correlation coefficient is located; Gaussian regression analysis is performed on sub-pixel points around the sub-pixel point pair to determine the target sub-pixel point pair.

4. The landslide deformation analysis method according to claim 3, characterized in that: The method of using one image in each image pair as a reference image and the other image as a search image, performing integer pixel search on each image pair, and determining the integer pixel pair at the location where the maximum value of the mutual correlation coefficient is located comprises: One image in each image pair is used as the reference image and the other as the search image; Use a first sliding window and a second sliding window of the same size to slide in the target area of ​​the reference image and the search area of ​​the search image respectively, and calculate the first mutual correlation coefficient in turn to obtain a first mutual correlation coefficient matrix; Determine the integer pixel point pair at the position where the maximum value of the first mutual correlation coefficient is located; The interpolation calculation is performed on the sub-pixel points around the integer pixel point pair to determine the sub-pixel point pair at the position where the maximum value of the mutual correlation coefficient is located, comprising: Meshing is performed at the location where the maximum value of the first mutual correlation coefficient is located; Using two third sliding windows and a fourth sliding window of the same size to slide in the grid, and using a bilinear interpolation algorithm to sequentially calculate the second mutual correlation coefficients of the third sliding window and the fourth sliding window at the gridded position, to obtain a second mutual correlation coefficient matrix; Performing Gaussian regression analysis on sub-pixel points around the sub-pixel point pair to determine the target sub-pixel point pair includes: Gaussian regression analysis is performed on the sub-pixel point at the position of the second cross-correlation coefficient peak and the 9 surrounding sub-pixels to obtain the optimal target sub-pixel point of the Gaussian regression model.

5. The landslide deformation analysis method according to claim 1, characterized in that: According to the elements of the deformation time matrix, the observation equation is solved using the least squares criterion or the weighted singular value decomposition method to obtain the rates of the two images in each image pair, including: When the deformation time matrix is ​​of full rank, the observation equation is solved using the least squares criterion to obtain the rates of the two images in each image pair; When the deformation time matrix is ​​rank-deficient, the observation equation is solved using a weighted singular value decomposition (SVD) method to obtain the rates of the two images in each image pair.

6. The landslide deformation analysis method according to claim 4, characterized in that: When the deformation time matrix is ​​rank-deficient, the observation equation is solved using the weighted singular value decomposition (SVD) method to obtain the rates of the two images in each image pair including: When the deformation time matrix is ​​rank-deficient, the observation equation is transformed using a weighted matrix to obtain a transformed observation equation; Perform singular value decomposition on the weight matrix and deformation time matrix in the transformed observation equation; Based on the singular value decomposition results of the weighted matrix and the deformation matrix, the transformed observation equation is solved to obtain the rates of the two images in each image pair.

7. The landslide deformation analysis method according to claim 6, characterized in that: The deformation observation matrix includes deformation observation values ​​of all image pairs, the deformation time matrix includes the time interval between two images in each image pair, and the deformation rate matrix includes: the rate of two images in each image pair; The observation equation is expressed as: Bv=d Wherein, B is an M×N coefficient matrix, and the elements in the matrix B consist of the time intervals between adjacent images in the image pair; m i represents the main image in the image pair, and the slave image is s i , the elements in matrix B can be expressed as: d is an M×1 deformation measurement matrix, the matrix elements are composed of the deformation observation values ​​of all cross-correlated image pairs, d = [d1 d2… d M ] T ; v is an N×1 matrix whose elements consist of the rates between adjacent images in an image pair, v = [v1 v2 … v n ] T .

8. The landslide deformation analysis method according to claim 7, characterized in that: The observation equation is solved using the least squares criterion, and the rates of the two images in each image pair are expressed as: v=(B T PB) -1 B T Pd Where P is an M×M weighted matrix, P=diag(ω1,ω2,…,ω M ),ω i The standard deviation σ of each image pair i The reciprocal of is an indicator to measure the weight of each image pair; 9. The landslide deformation analysis method according to claim 7, characterized in that: The transformed observation equation is expressed as: The results of singular value decomposition of the weight matrix and coefficient matrix are expressed as: Among them, svd is the singular value decomposition function; U is an M×M orthogonal matrix, which is ZZ T Eigenvector; V is an N×M matrix, which is Z T The eigenvector of Z; S is an M×M diagonal matrix; The transformed observation equation is solved using the singular value decomposition result, and the rate of the two images in each image pair is expressed as: Among them, S + is the pseudo-inverse of S, given by the singular value σ i Composition, S + =diag(1 / σ1,1 / σ2,…,1 / σ N-L+1 , 0, 0, …, 0), the number of singular values ​​that are 0 is (L-1).

10. The landslide deformation analysis method according to claim 1, characterized in that: After obtaining the deformation observation matrix of each group of image pairs, the landslide deformation analysis method further includes: Decomposing the deformation observation matrix in the east-west direction and the north-south direction to obtain deformation observation values ​​in the east-west direction and the north-south direction; Calculate the signal-to-noise ratio of the image pairs used to obtain the deformation observations in the east-west and north-south directions; Using the signal-to-noise ratio to verify the deformation measurement matrix; Filtering the deformation observation values ​​in the deformation observation matrix to obtain smooth deformation observation values; The smoothed deformation observations are decomposed into east-west and north-south directions to obtain the deformation observations in the east-west and north-south directions.

Images