A landslide monitoring method based on SAR image and optical image fusion in complex mountainous areas

Abstract

The present invention discloses a landslide monitoring method for complex mountainous areas by fusing SAR images with optical images, including: using the R index method and GIS mountain shadow simulation method to quickly identify geometric distortion areas, analyzing the geometric relationship between the InSAR line of sight direction and the terrain, projecting the deformation results in the LOS direction to the slope direction, identifying the incoherent area by using SAR image intensity information, judging the InSAR identification results in combination with optical images, and fusing different remote sensing data. The method improves the resolution of landslide identification and the reliability of identification results.

Description

Technical Field

[0001] The present invention relates to the technical field of landslide identification, and more specifically to a landslide monitoring method by fusing SAR images with optical images in complex mountainous areas. Background Art

[0002] In recent years, geological disasters have occurred frequently in mountainous areas of my country. Among the many geological disasters, landslide disasters account for the highest proportion of geological disasters each year, accounting for nearly 49% of the total disasters in 2021, and landslides are one of the hot spots of social concern. According to the "Landslide Prevention Engineering Investigation Code" (GB / T 32864-2016) formulated by the China Geological Survey, landslides can be divided into five levels according to their volume: small, medium, large, extra-large and giant. With the development of multi-source data and multi-platform images, landslide identification in complex mountainous areas is mainly divided into two types of technologies: optical and InSAR. Optical mainly uses remote sensing interpretation methods, and InSAR technology mainly uses small baseline set technology (SBAS), permanent scatterer technology (PS), offset tracking technology (Offset-tracking) and interferogram stacking technology (Stacking). For landslide identification, using only SAR data will lead to problems such as geometric distortion, spatiotemporal incoherence, and long revisit periods; using only optical data will be affected by weather factors. Therefore, by combining optical and SAR data, landslides can be identified more completely. However, there is currently little research on comprehensive identification methods for complex mountain landslides using combined SAR data and optical data. Summary of the invention

[0003] The embodiment of the present invention provides a landslide monitoring method for complex mountainous areas by fusing SAR images with optical images, comprising:

[0004] Acquire satellite SAR image data in mountainous areas;

[0005] The terrain index R is calculated based on the digital elevation model DEM and SAR image data, and the active geometric distortion area of ​​the SAR image is obtained according to the R index change curve under different slopes;

[0006] Calculate the incident angle of the satellite line of sight LOS, based on the mountain shadow function in the GIS map and the relationship between the azimuth angle and the solar altitude angle, and calculate the passive geometric distortion area of ​​the SAR image;

[0007] Eliminate the active and passive geometric distortion areas in the SAR image, and use the Stacking technology to obtain the annual average deformation rate in the LOS direction of the eliminated geometric distortion area;

[0008] ArcGIS software was used to evaluate the accuracy of the annual average deformation rate in the LOS direction. The quality of the Stacking deformation monitoring results was determined based on the error in the annual average deformation rate in the LOS direction. The spatial analysis function of ArcGIS software was used to delineate the suspected landslide risk points in the Stacking deformation monitoring results.

[0009] Superimpose the suspected landslide potential points onto the optical image points to obtain the confirmed landslide potential points;

[0010] Project the monitoring results of InSAR line of sight LOS to the slope direction to obtain the deformation rate along the slope and monitor the movement characteristics of landslide potential points;

[0011] Through SBAS technology and optical image change analysis, the temporal change characteristics and optical image characteristics of the movement characteristics of the landslide hazard point are determined to determine whether the landslide hazard point should be continuously monitored in the future.

[0012] Preferably, the SAR image data is preprocessed, which includes:

[0013] Calculate the offset of the pixel to be registered in the SAR image relative to the corresponding pixel in the master image in the azimuth and range directions, and obtain the registered slave image through offset polynomial estimation;

[0014] Convert SAR coordinate system to geographic coordinate system.

[0015] Preferably, the method further comprises performing differential interferometry calculation on the pre-processed SAR image data, which comprises:

[0016] The temporal baseline is determined according to the coherence characteristics of the study area, and the spatial baseline is determined according to the degree of terrain undulation;

[0017] Calculate the phase integer ambiguity of each pixel in the interference pattern and perform phase unwrapping;

[0018] Reduce baseline error through baseline refinement;

[0019] To untangle.

[0020] Preferably, the terrain index R is calculated by the formula:

[0021]

[0022] In the formula: αα represents the satellite's incident angle, ββ represents the slope, θθ represents the satellite's heading angle, and γγ represents the slope direction. When the slope is greater than the satellite's incident angle, active overlap will occur. When the terrain index R is positive and less than sin(α)sin(α), perspective contraction will occur. When the terrain index R = sin(α)sin(α), the surface is flat.

[0023] Preferably, the LOS direction incident angle is calculated, and the calculation formula includes:

[0024]

[0025] In the formula, α i α i represents the i-th column incidence angle, H is the distance from the satellite to the center of the earth in m, h is the distance from the satellite to the geoid, S is the slant range resolution, and L is the near slant range.

[0026] Preferably, the deformation result in the LOS direction is projected onto the slope direction to obtain the deformation rate along the slope, including:

[0027] Calculate the LOS direction unit vector R;

[0028]

[0029] Re, Rn, Rz, and R represent the unit vectors in the east-west, north-south, vertical, and LOS directions, respectively;

[0030] The landslide is along the direction of the steepest slope unit vector U;

[0031] β mE β mE , β mN β mN , β mZ β mZ The unit vectors representing the steepest slope directions in the east-west, north-south, and vertical directions, respectively, with positive values ​​indicating upward and negative values ​​indicating downward;

[0032] According to the Euclid principle, the projection V of the landslide movement direction in the LOS direction is obtained. L V L ;

[0033]

[0034] V L -KU=0

[0035] V L =V E R E +V N R N +V Z R Z

[0036]

[0037] AX=V

[0038] X=(A) -1 V=(β mE V E+β mN V N +β mZ V Z ) -1 V

[0039] Define the angle between the LOS direction and the slope as A, and the apparent slope angle σσ to obtain the deformation rate along the slope:

[0040]

[0041] Where K represents the average velocity modulus along the slope, βmβm represents the deformation rate in the slope direction, and V L V L represents the deformation rate in the LOS direction, αα represents the slope direction, ββ represents the slope, and θθ represents the incident angle of the satellite. Indicates the angle between the satellite flight direction and the true north direction.

[0042] Preferably, after identifying the landslide hazard points in the SAR image, the SBAS technology is used to obtain the time series results.

[0043] Preferably, the deformation rate result of the sight line direction is obtained by stacking technology to prove the reliability of the monitoring result, which specifically includes:

[0044] Select the stable area in the InSAR monitoring results;

[0045] Use the spatial analysis function of ArcGIS software to determine the stable areas in the InSAR monitoring results;

[0046] Calculate the standard deviation of the stable area and eliminate the monitoring result values ​​that exceed 3 times the standard deviation;

[0047] According to the error statistics, the reliability of the monitoring results is calculated.

[0048] Preferably, the method further includes using an offset tracking technology to identify landslides for large-scale deformations, which includes:

[0049] Using the two registered intensity images, the total amount of displacement caused by surface deformation and satellite orbit offset is determined by the peak of the cross-correlation coefficient within the window.

[0050] The satellite orbit data is used to remove the offset caused by orbital deviation, thereby obtaining the two-dimensional deformation of the surface in range and azimuth.

[0051] The formula for the mutual correlation coefficient within the window is as follows:

[0052]

[0053] Where: ρ(x, y) is the normalized correlation coefficient, (x, y) is the reference pixel position, (xu, yu) is the search pixel position, r and s are the pixel values ​​of the reference pixel and the search pixel respectively, u r and u s are the average pixel values ​​of the reference and search windows, respectively.

[0054] The embodiment of the present invention provides a landslide monitoring method for complex mountainous areas by fusing SAR images with optical images. Compared with the prior art, the method has the following beneficial effects:

[0055] 1. Areas with large terrain undulations will cause obvious geometric distortion in SAR images. The identification of geometric distortion can help analyze the reasons why InSAR technology fails to identify geological disasters. It can also help determine whether the suspected hidden danger points identified can become real hidden danger points.

[0056] 2. There are certain differences between the deformation results of the InSAR line of sight direction and the actual movement direction of the geological disaster, and there are certain differences between the monitoring results of the line of sight direction and the results of the actual movement direction. By converting the line of sight direction into slope direction, it is helpful to determine the movement direction of the landslide.

[0057] 3. In order to determine whether a suspected hidden danger point can become a real hidden danger point, it is necessary to combine optical images to determine whether the hidden danger point poses a threat to residential areas and public facilities. Combining multiple optical images can also clearly show the deformation process of the landslide. BRIEF DESCRIPTION OF THE DRAWINGS

[0058] Figure 1 A technical scheme diagram of landslide comprehensive identification by fusion of optics and SAR for a landslide monitoring method by fusion of SAR images and optical images in a complex mountainous area provided by an embodiment of the present invention;

[0059] Figure 2 A landslide monitoring method for complex mountainous area by fusing SAR image and optical image provided by an embodiment of the present invention, and a sentinel data R index curve;

[0060] Figure 3 An active geometric distortion diagram of a landslide monitoring method for fusing SAR images and optical images in a complex mountainous area provided by an embodiment of the present invention;

[0061] Figure 4 A landslide monitoring method for complex mountainous area by fusing SAR images and optical images provided by an embodiment of the present invention is a graph showing changes in R index curves under different slopes;

[0062] Figure 5 A landslide monitoring method for complex mountainous area SAR image and optical image fusion provided by an embodiment of the present invention is a graph showing the relationship between the incident angle and the column;

[0063] Figure 6 A passive geometric distortion calculation diagram of a landslide monitoring method for fusing SAR images and optical images in a complex mountainous area provided by an embodiment of the present invention;

[0064] Figure 7 A data graph of the Sentinel orbit rise of a landslide monitoring method for complex mountainous areas by fusing SAR images and optical images provided by an embodiment of the present invention;

[0065] Figure 8 A slope aspect result diagram of a landslide monitoring method for complex mountainous areas by fusing SAR images with optical images provided by an embodiment of the present invention;

[0066] Fig. 9 An optical image change diagram of a landslide monitoring method for fusing SAR images and optical images in complex mountainous areas provided by an embodiment of the present invention. DETAILED DESCRIPTION

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

[0068] Example 1 See Figure 1 ;

[0069] This patent uses Sentinel-1 data as an example to introduce:

[0070] The acquired Sentinel-1 images are preprocessed, reasonable time-space thresholds are set to select image pairs, and differential interferometry calculations are performed. For the identification of geometric distortion, the geometric distortion of SAR images depends on the satellite LOS direction, ground slope and slope aspect, so the terrain index R method based on DEM and SAR system parameters can be used. The specific principles and steps are as follows:

[0071]

[0072] In the formula: αα represents the satellite's incident angle, ββ represents the slope, θθ represents the satellite's heading angle, and γγ represents the slope direction. When the slope is greater than the satellite's incident angle, active overlap will occur; when the terrain index R is positive and less than sin(α)sin(α), perspective contraction will occur; when the terrain index R = sin(α)sin(α), the surface is flat.

[0073] The above method only considers the influence of active geometric distortion, but does not calculate the influence of passive geometric distortion. Especially in mountainous areas, the visibility between slopes will be affected. The passive geometric distortion area can be calculated by the mountain shadow function in GIS. Assuming that the incident angle of the radar is a fixed value, the satellite is regarded as the sun at infinity, and the geometric distortion area is calculated based on the relationship between the azimuth angle and the solar altitude angle.

[0074]

[0075] In the formula, α i α i represents the i-th column incidence angle, H is the distance from the satellite to the center of the earth (m), h is the distance from the satellite to the geoid, S is the slant range resolution, and L is the near slant range.

[0076] The LOS turning slope direction deformation rate method assumes that Re, Rn, Rz, and R represent the unit vectors in the east-west, north-south, vertical, and LOS directions, respectively. That is:

[0077]

[0078] Assume that U represents the direction of the unit vector of the landslide along the steepest slope, β mE , β mN , β mZ Represent the unit vectors of the steepest slope in the east-west, north-south, and vertical directions respectively. Positive values ​​indicate upward and negative values ​​indicate downward. They are decomposed according to geometric relationships:

[0079]

[0080] Once the velocity vector in the coordinate system (E, N, Z) is defined, the landslide motion direction is projected in the LOS direction V L V L , according to Euclid's principle, it can be expressed as:

[0081]

[0082] V L -KU=0 (6)

[0083] V L =V E R E +V N R N +V Z R Z (7)

[0084]

[0085] AX=V (9)

[0086] X=(A)-1 V=(β mE V E +β mN V N +β mZ V Z ) -1 V (10)

[0087] Define the apparent slope angle σσ (the angle between the LOS direction and the slope) as A, and get the deformation rate along the slope, that is, Error: Reference source not found:

[0088]

[0089] Where K is the average velocity modulus along the slope; βmβm is the deformation rate in the slope direction (mm / a); V L V L represents the deformation rate in the LOS direction (mm / a); αα represents the slope direction; ββ represents the slope; θθ represents the incident angle of the satellite; Indicates the angle between the satellite's flight direction and true north.

[0090] From the analysis of formula 11, it can be found that when the value of σσ reaches 90°, the deformation rate of the slope will reach infinity, so the value of cosσcosσ needs to be limited. That is:

[0091]

[0092] Since the movement direction of the landslide on the slope is downward along the slope, V p V p At this time, it is a negative value; when V p V p When it is a positive value, it means that the movement direction of the landslide on the slope is upward.

[0093] Example 2, see Figure 2-9 :

[0094] This patent uses the Jiali Fault area as the study area to conduct experiments and analyze the technical effects of geometric distortion and LOS turning slope direction.

[0095] Figure 2 The R index variation curves at different slopes obtained by using the R index calculation method. Due to the large terrain undulations in the area, the maximum relative height difference in the study area exceeds 1 km, which causes geometric distortion in most areas on SAR images. Therefore, it is necessary to identify geometric distortion. Figure 4 This is the R index change curve under different slopes. As the slope increases, the change amplitude of the R index curve gradually increases.

[0096] The satellite parameters show that the distance from the satellite to the center of the earth is 7071741m, the distance from the satellite to the geoid is 6372337m, S is the slant range resolution of 2.33m, and L is the near slant range of 799135m. Using the above parameters, the incidence angle variation curve of the area is calculated. Figure 5 It can be found that for a complete SLC image, the incident angle range is (28.93°, 43.14°), and the incident angle difference between the near slant range and the far slant range exceeds 14°. In order to quickly calculate the passive geometric distortion area, the length and width of the SLC image are calculated to be 255×183km, and the interval is calculated in GIS at an interval of 1°, about 18km. Finally, the geometric distortion of the entire image can be estimated by combining the active and passive geometric distortion areas.

[0097] The Stacking technology is used to obtain the deformation rate results of the line of sight direction of the study area, in mm / a. Figure 7 A positive value indicates that the area is close to the satellite line of sight, while a negative value indicates that the area is far from the satellite line of sight. The terrain conditions in the study area are complex and there is a lack of field data. In order to evaluate the accuracy of the InSAR deformation monitoring results, the stable area ( Figure 7 The standard deviation of the calculated stable area is 2.47 mm / a. According to the statistical law of error, the reliability of the monitoring results is proved.

[0098] The deformation results of different track directions (ascending and descending) show opposite results in the deformation direction, which makes it impossible to effectively express the actual movement status of the potential hazard point, which is not conducive to the identification of the potential hazard point. Therefore, projecting the deformation results in the LOS direction to the slope direction can truly reflect the deformation and facilitate the identification of the potential hazard point ( Figure 8 ).

[0099] After the suspected landslide hazard points are delineated, remote sensing images of different time periods are selected to determine the deformation characteristics of the landslide. For example, the Longpu landslide is located near Wada Village, Basu County, Tibet, with geographic coordinates (97°7′10.96″E, 29°58′23.36″N). The landslide moves downward along the sliding surface of the landslide body. By comparing the optical images at different times, it can be found that from January 2016 to December 2018, the landslide body did not show obvious deformation in the optical images, but after December 2018, it was found that the landslide body had obvious deformation information, Figure (9).

[0100] Overall:

[0101] There will be obvious geometric distortion in areas with large terrain undulations. The identification of geometric distortion can help analyze the reasons why InSAR technology fails to identify geological hazards. It can also help determine whether the suspected hidden danger points identified can become real hidden danger points.

[0102] There is a certain difference between the deformation results of the InSAR line of sight direction and the actual movement direction of the geological disaster, and the monitoring results of the line of sight direction will be smaller than the results of the actual movement direction. By converting the line of sight direction into slope direction, it is helpful to determine the movement direction of the landslide.

[0103] In order to determine whether a suspected hidden danger point can become a real hidden danger point, it is necessary to combine optical images to determine whether the hidden danger point poses a threat to residential areas and public facilities. Combining multiple optical images can also clearly show the deformation process of the landslide.

[0104] Disclosed above are only a few specific embodiments of the present invention. Those skilled in the art may make various changes and modifications to the embodiments of the present invention without departing from the spirit and scope of the present invention. However, the embodiments of the present invention are not limited thereto. Any changes that can be conceived by those skilled in the art should fall within the scope of protection of the present invention.

Claims

1. A landslide monitoring method for complex mountainous areas by fusing SAR images with optical images, characterized in that: include: Acquire satellite SAR image data in mountainous areas; The terrain index R is calculated based on the digital elevation model DEM and SAR image data, and the active geometric distortion area of ​​the SAR image is obtained according to the R index change curve under different slopes; Calculate the incident angle of the satellite line of sight LOS, based on the mountain shadow function in the GIS map and the relationship between the azimuth angle and the solar altitude angle, and calculate the passive geometric distortion area of ​​the SAR image; Eliminate the active and passive geometric distortion areas in the SAR image, and use the Stacking technology to obtain the annual average deformation rate in the LOS direction of the eliminated geometric distortion area; ArcGIS software was used to evaluate the accuracy of the annual average deformation rate in the LOS direction. The quality of the Stacking deformation monitoring results was determined based on the error in the annual average deformation rate in the LOS direction. The spatial analysis function of ArcGIS software was used to delineate the suspected landslide risk points in the Stacking deformation monitoring results. Superimpose the suspected landslide potential points onto the optical image points to obtain the confirmed landslide potential points; Project the monitoring results of InSAR line of sight LOS to the slope direction to obtain the deformation rate along the slope and monitor the movement characteristics of landslide potential points; Through the SBAS technology and optical image change analysis, the temporal change characteristics and optical image characteristics of the movement characteristics of the landslide potential point are determined to determine whether the landslide potential point should be continuously monitored in the future; The monitoring result of the InSAR line of sight LOS is projected onto the slope direction to obtain the deformation rate along the slope, including: Calculate the LOS direction unit vector R; Re, Rn, Rz, and R represent the unit vectors in the east-west, north-south, vertical, and LOS directions, respectively; The landslide is along the direction of the steepest slope unit vector U; β mE , β mN , β mZ The unit vectors representing the steepest slope directions in the east-west, north-south, and vertical directions, respectively, with positive values ​​indicating upward and negative values ​​indicating downward; According to the Euclid principle, the projection V of the landslide movement direction in the LOS direction is obtained. L ; V L -KU=0 V L =V E R E +V N R N +V Z R Z AX=V X=(A) -1 V=(β mE V E +β mN V N +β mZ V Z ) -1 V Define the angle between the LOS direction and the slope as A, and the apparent slope angle σ, and obtain the deformation rate along the slope: Where K represents the average velocity modulus along the slope, βm represents the deformation rate in the slope direction, and V L represents the deformation rate in the LOS direction, α represents the slope direction, β represents the slope, θ represents the incident angle of the satellite, and θ represents the angle between the satellite flight direction and the true north direction.

2. The landslide monitoring method of complex mountainous area by fusing SAR image and optical image as claimed in claim 1, characterized in that: It also includes SAR image data preprocessing, which includes: Calculate the offset of the pixel to be registered in the SAR image relative to the corresponding pixel in the master image in the azimuth and range directions, and obtain the registered slave image through offset polynomial estimation; Convert SAR coordinate system to geographic coordinate system.

3. The landslide monitoring method of complex mountainous area by fusing SAR image and optical image as claimed in claim 2, characterized in that: It also includes differential interferometry calculation of pre-processed SAR image data, which includes: The temporal baseline is determined according to the coherence characteristics of the study area, and the spatial baseline is determined according to the degree of terrain undulation; Calculate the phase integer ambiguity of each pixel in the interference pattern and perform phase unwrapping; Reduce baseline error through baseline refinement; To untangle.

4. The landslide monitoring method of complex mountainous area by fusing SAR image and optical image as claimed in claim 1, characterized in that: The terrain index R is calculated by the following formula: In the formula: α represents the satellite's incident angle, β represents the slope, θ represents the satellite's heading angle, and γ represents the slope direction. When the slope is greater than the satellite's incident angle, active overlap will occur. When the terrain index R is positive and less than sin(α), perspective contraction will occur. When the terrain index R = sin(α, the surface is flat.

5. The landslide monitoring method of complex mountainous area by fusing SAR image and optical image as claimed in claim 1, characterized in that: The calculation formula for calculating the incident angle of the satellite line of sight direction LOS includes: In the formula, α i represents the i-th column incidence angle, H is the distance from the satellite to the center of the earth in m, h is the distance from the satellite to the geoid, S is the slant range resolution, and L is the near slant range.

6. The landslide monitoring method of complex mountainous area by fusing SAR image and optical image as claimed in claim 1, characterized in that: After the landslide potential points are determined, the SBAS technology is used to obtain the time series results.

7. The landslide monitoring method of complex mountainous area by fusing SAR image and optical image as claimed in claim 1, characterized in that: The stacking technology is used to obtain the annual average deformation rate in the LOS direction of the geometric distortion area, specifically including: Select the stable area in the InSAR monitoring results; Using the spatial analysis function of ArcGIS software, the stable areas in the InSAR monitoring results are analyzed. Calculate the standard deviation of the stable area and eliminate the monitoring result values ​​that exceed 3 times the standard deviation; According to the error statistics, the reliability of the monitoring results is calculated.

8. The landslide monitoring method of complex mountainous area by fusing SAR image and optical image as claimed in claim 1, characterized in that: It also includes the use of offset tracking technology for landslide identification for large-scale deformation, including: Using the two registered intensity images, the total amount of displacement caused by surface deformation and satellite orbit offset is determined by the peak of the cross-correlation coefficient within the window. The satellite orbit data is used to remove the offset caused by orbital deviation, thereby obtaining the two-dimensional deformation of the surface in range and azimuth. The formula for the mutual correlation coefficient within the window is as follows: Where: ρ(x,y) is the normalized correlation coefficient, (x,y) is the reference pixel position, (xu,yu) is the search pixel position, r and s are the pixel values ​​of the reference pixel and the search pixel respectively, u r and u s are the average pixel values ​​of the reference and search windows, respectively.

Images