CN113281742B - SAR landslide early warning method based on landslide deformation information and meteorological data - Google Patents

Abstract

The invention discloses an SAR landslide early warning method based on landslide deformation information and meteorological data, which comprises the following steps: s1, determining a rainfall threshold value of a landslide caused by effective rainfall; s2, utilizing the rainfall threshold value and the regional rainfall information to perform landslide early warning; s3, acquiring time sequence deformation data of the early warning area; s4, calculating a critical deformation threshold value of landslide occurrence; s5, checking the landslide early warning area obtained in the step S2 and the time sequence deformation data obtained in the step S3 by utilizing a critical deformation threshold value, and realizing the secondary early warning of the early warning area. Aiming at the complex characteristics of special sudden property, concealment, uncertainty and the like of landslide, the invention realizes the generalization of the relation among landslide, rainfall and deformation by researching the relativity between the time sequence rainfall and the time sequence deformation of the landslide, and can accurately determine whether the target is a landslide early warning area.

Description

A SAR Landslide Early Warning Method Based on Landslide Deformation Information and Meteorological Data

technical field

The invention relates to the technical field of synthetic aperture radar monitoring of geological disasters, in particular to a SAR landslide early warning method based on landslide deformation information and meteorological data.

Background technique

my country is a mountainous country, and the mountainous area accounts for 69% of the total land area. Common geological disasters in mountainous areas pose a great threat to road construction and operation. Landslide refers to the phenomenon that the rock and soil mass on the slope loses stability under the influence of natural or human factors, and slides down along the entire failure surface. As a common geological disaster, it is widely distributed and harmful. According to incomplete statistics, in addition to the casualties caused by landslide disasters, landslide disasters have caused thousands of kilometers of roads to suffer varying degrees of damage in the past 10 years, causing hundreds of millions of economic losses. With the development of our country's economy, the national highway construction will also continue to increase, and will continue to extend to the mountainous areas, which will mean that road disasters will also become increasingly serious. Through the experience and lessons of a series of landslide disaster events, it can be seen that early warning of landslide disasters and effective disaster analysis are the main ways to change "passive disaster avoidance and relief" to "active disaster prevention and treatment" and reduce the losses caused by disasters.

Traditional surface deformation monitoring methods mostly focus on surface monitoring technology and underground monitoring technology, such as geodetic method, GPS method, drilling point monitoring and geophysical detection technology. Traditional landslide monitoring methods are not continuous in time and space, and in the face of inaccessible and difficult-to-reach complex and dangerous mountainous areas, ground and underground monitoring technologies are faced with difficult communication, battery life, and difficult deployment and installation problems. In recent years, with the rapid development of communication technology and sensor technology, aerospace remote sensing monitoring technology has emerged as the times require. Optical remote sensing, spaceborne InSAR technology, aerial photogrammetry, and airborne LiDAR technology have transformed landslide disaster monitoring from the previous point-based monitoring method to the surface-based dynamic monitoring method. However, landslide disasters occur occasionally and are often accompanied by severe weather conditions, even at night. Optical remote sensing, aerial photogrammetry, and laser LiDAR technology are difficult to obtain image data with high resolution and high timeliness. The all-weather and all-weather radar can penetrate clouds and forests to obtain images, and then has the incomparable advantages of optical images.

Interferometric Synthetic Aperture Radar (InSAR) is a quantitative microwave remote sensing technology developed in the past half century. The three-dimensional information of the surface can be extracted by using the interferometric phase map and the attitude data of the platform equipped with radar sensors, and the change information of the land surface can be extracted by using the interferometric coherence analysis. From the perspective of surface deformation detection, the second difference technology can be used to remove the influence of terrain and other factors from the interferometric phase map, thereby extracting deformation information, which is the differential interferometry technology of synthetic aperture radar (D-InSAR). This makes InSAR technology have a broader application prospect in the field of landslide disaster warning and monitoring. At the same time, according to the statistics of sudden landslide disasters in China in recent years, continuous rainfall-induced landslides accounted for 65% of the total landslide occurrence, of which local rainfall-induced landslides accounted for 43% of the total occurrence, and continuous rainfall-induced landslides accounted for 66%. In other words, about two-thirds of sudden landslide disasters are closely related to atmospheric rainfall or meteorological factors. Therefore, combining InSAR technology with landslide inducing factors (rainfall) is an effective way to carry out landslide disaster early warning.

Contents of the invention

The purpose of the present invention is to provide a SAR landslide early warning method based on landslide deformation information and meteorological data, so as to overcome the defects in the prior art.

In order to achieve the above object, the technical scheme adopted in the present invention is as follows:

A SAR landslide early warning method based on landslide deformation information and meteorological data, comprising the following steps:

S1. Calculate the correlation coefficient between the number of historical landslide occurrences and the time-series rainfall before the landslide occurs, and determine the rainfall threshold that the effective rainfall causes the landslide;

S2. Correct the rainfall threshold obtained in step S1 according to the correlation coefficient between the time-series effective rainfall and the time-series deformation data;

S3, using the rainfall threshold and regional rainfall information to perform early warning of landslides;

S4. Obtain the time-series deformation data of the initial warning area;

S5. According to the relationship between the number of historical landslide occurrences and the time-series cumulative deformation data, determine the effective deformation amount that causes the landslide, and calculate the critical deformation threshold for the landslide accordingly;

S6. Use the critical deformation threshold to check the landslide initial warning area obtained in step S2 and the time-series deformation data obtained in step S3, so as to realize the secondary warning of the initial warning area.

Further, the calculation method of the effective rainfall in the step S1 is:

γ _n =γ ₀ +Tγ ₁ +T ² γ ₂ +T ³ γ ₃ +...+T ⁿ γ _n ;

In the formula, γ _n represents the effective rainfall; γ ₀ represents the rainfall of the day; γ ₁ , γ ₂ , γ ₃ , ..., γ _n represent the rainfall before the current day; n is the number of days passed;

Further, the step S3 specifically includes: judging the area where the rainfall threshold is greater than the set value as a landslide-prone area, and judging the area where the rainfall threshold is less than or equal to the set value as a stable area.

Further, the step S4 is to obtain the time-series deformation data of the landslide-prone area by performing time-series InSAR processing on the area.

Further, the time-series deformation data obtained by performing time-series InSAR processing on landslide-prone areas is specifically:

It is assumed that the jth interferometric phase map is expressed in the azimuth-distance pixel coordinate system (x, r) as:

In the formula, λ is the wavelength of the radar; d(t _B , x, r) and d(t _A , x, r) are the accumulated deformation of the line of sight defense corresponding to the reference time t ₀ at time t _B and t _A respectively, and d(t ₀ , x, r)≡0; d(t _i , x, r), (i=1, 2, 3...N) is used to represent the deformation time series, and the corresponding phase is There are: />

The M phase values corresponding to the deformation of each pixel point in the study area are represented by vectors

will be calculated from the differential interferogram values are represented as vectors:

Wherein, main image IE=[IE ₁ , IE ₂ , IE ₃ ,..., IE _M ], secondary image IS=[IS ₁ , IS ₂ , IS ₃ ,..., IS _M ];

On the basis of the above formula, the vector phase of the interferogram is expressed in the form of a matrix:

In the matrix, each row corresponds to each differential interferogram, and each column corresponds to SAR images at different times. The columns of the main image and secondary images in the matrix are ±1, and the remaining columns are 0.

like That is, the matrix G of order [M×N] can be expressed as

When a series of interference pairs are generated in the same small baseline subset, M≥N and N is the rank of G, the phase matrix can be obtained by the least square method:

For the phenomenon that it is impossible for the interference pairs of all images to be in the same baseline subset, the existing SARs are combined into a certain number of subsets according to the threshold, and then G ^T G is a matrix after rank reduction. If the matrix is not satisfied with the rank matrix, the singular value decomposition method is used to decompose the singular matrix into: G=USV ^T , where V ^T is the velocity of the average phase, U is the orthogonal matrix, and S is the diagonal matrix, and then the phase can be converted into the average phase velocity:

Then obtain the minimum norm solution of the velocity vector V, and integrate to obtain the estimated value of the phase, which is the accurate deformation phase vector δ.

Further, the step S5 is specifically: obtaining the deformation threshold by calculating the correlation coefficient between the time-series cumulative deformation amount when the landslide occurs and the number of occurrences of the landslide, and obtaining the corresponding deformation amount when the number of landslides jumps according to the relationship diagram between the two, using the deformation value corresponding to the first jump as the warning-level deformation threshold, and the deformation value corresponding to the second jump. The deformation threshold is the alarm-level deformation threshold.

Further, the deformation criterion model of the landslide area is trained according to the deformation data of the historical landslide area, and the accurate deformation threshold of the warning level and the deformation threshold of the warning level are obtained according to the model.

Further, the secondary warning of the preliminary warning area specifically includes determining the activity level of the preliminary warning area according to the relationship between the time series deformation data of the preliminary warning area and the deformation threshold of the warning level and the deformation threshold of the warning level.

Compared with the prior art, the advantage of the present invention is that: the present invention aims at complex characteristics such as suddenness, concealment and uncertainty unique to landslides, and realizes the induction of the relationship between landslides, rainfall and deformation by studying the correlation between landslide time-series rainfall and time-series deformation variables, and can accurately determine whether the target is a landslide early warning area.

Description of drawings

In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following will briefly introduce the accompanying drawings required in the description of the embodiments or prior art. Obviously, the accompanying drawings in the following description are only some embodiments of the present invention. For those of ordinary skill in the art, other accompanying drawings can also be obtained according to these drawings without creative work.

Fig. 1 is a flowchart of the SAR landslide early warning method based on landslide deformation information and meteorological data in the present invention.

Figure 2 is a scatter diagram of rainfall in the landslide area.

Fig.3 The relationship curve between effective rainfall and landslide occurrence times

Figure 4 is a preliminary early warning criterion diagram.

Fig. 5 is a flow chart of acquiring time-series deformation data of the target area.

Figure 6 is the deformation data map of the study area.

Figure 7 is a schematic diagram of two early warnings.

Detailed ways

The preferred embodiments of the present invention will be described in detail below in conjunction with the accompanying drawings, so that the advantages and features of the present invention can be more easily understood by those skilled in the art, so as to define the protection scope of the present invention more clearly.

Referring to shown in Fig. 1, it is a kind of typical SAR landslide early warning method based on landslide deformation information and meteorological data that the present embodiment discloses, comprises the following steps:

Step S1. Calculate the correlation coefficient between the number of historical landslide occurrences and the time-series rainfall before the landslide occurs, and determine the rainfall threshold for the effective rainfall to trigger the landslide. Specifically:

When the frequency of rainfall is positively correlated with the amount of rainfall, the possibility of landslides will increase. In this embodiment, the rainfall data of the rainfall observation station near the existing landslide area are collected, and the rainfall data of the day are multiplied by the effective rainfall coefficient of the day for a certain period of time before the landslide to obtain the effective rainfall.

γ _n =γ ₀ +Tγ ₁ +T ² γ ₂ +T ³ γ ₃ +...+T ⁿ γ _n (1);

In the formula, γ _n represents the effective rainfall; γ ₀ represents the rainfall of the day; γ ₁ , γ ₂ , γ ₃ , ..., γ _n represent the rainfall before the day; n is the number of days passed;

In this embodiment, 500 landslide points in Sichuan Province were obtained, and the time-series rainfall data of the rainfall monitoring station closest to each landslide point were collected, and the relationship between the number of occurrences of landslides and rainfall was statistically analyzed, and a scatter diagram of the relationship between the number of occurrences of landslides and different effective rainfall was formulated. The scatter diagram of rainfall in selected landslide areas is shown in Figure 2. According to the relationship curve between the effective rainfall and the number of landslides (as shown in Figure 3), it can be seen that when the effective rainfall reaches 60mm, the curve becomes steeper, the slope increases, and the number of landslides jumps. At this time, the correlation coefficient has also exceeded half. When the effective rainfall reaches 120mm and 200mm, similar phenomena also occur. Therefore, this paper sets the effective rainfall threshold D as 60mm, 120mm, and 200mm respectively, corresponding to the initial landslide warning level: attention level, warning level, and alarm level. Then obtain the preliminary warning criterion model of landslide hazards, as shown in Figure 4, which is divided into the lowest critical line and the middle critical line, and three areas A, B, and C. Below the lowest critical line of rainfall is the non-warning area, that is, A area; higher than the lowest critical line and lower than the middle critical line (ie, B area) is the candidate early warning area; higher than the middle critical line (ie, C area) is the key early warning area.

Step S2, calculate the correlation between the time-series rainfall and the time-series deformation according to the formula (2), correct the effective rainfall threshold obtained by S1, and obtain the contribution of the effective rainfall to the deformation of the landslide area according to the level of the correlation coefficient between the two, so as to verify the effectiveness of the effective rainfall threshold. In this embodiment, the correlation coefficient obtained according to formula (2) is consistent with the effective rainfall threshold value obtained by S1.

Step S3, using the rainfall threshold and regional rainfall information to perform preliminary warning of landslides, specifically: performing preliminary warning of the target area according to the effective rainfall threshold. Firstly, the time-series rainfall data of the target area is obtained, and the effective rainfall at each time point is calculated according to the above formula (1), and the preliminary early warning of the target area is carried out according to the formula (3), to determine whether the target is a landslide-prone area.

Step S4. Obtain the time-series deformation data of the initial warning area, specifically: determine the landslide-prone area through step S3, perform time-series InSAR processing on the SAR data of this area, and obtain the time-series deformation data of this area. This embodiment adopts SBAS-InSAR technology. First, set the threshold of space-time baseline to obtain interference image pairs; then, image registration is performed, and all images are registered to the super master image; then coherence calculation and interferogram filtering are performed to generate a coherence coefficient map and a filtered secondary differential interferogram; the unwrapping of the filtered secondary differential interferogram is achieved by minimizing the difference between the unwrapped phase gradient and the real phase gradient; Remove the residual orbital phase and horizon phase; on this basis, statistically analyze the phase of the high-coherence point obtained by the coherence coefficient threshold method, and model and solve each coherent point, and use the singular value decomposition method to obtain the linear deformation phase and the elevation error phase of the high-coherence area; finally separate the residual phase, perform secondary unwrapping on the residual phase, and use the spatial domain filter to obtain the atmospheric phase. So far, each component in the interference phase has been obtained. The linear deformation phase and the nonlinear deformation phase of .

The main processing flow of SBAS-InSAR technology is shown in Figure 5.

Step S5, according to the correlation coefficient between the time-series effective rainfall and time-series deformation data, determine the effective deformation amount that triggers the landslide, and calculate the critical deformation threshold for the occurrence of the landslide accordingly, specifically:

In this embodiment, the time-series SAR data of the existing landslide area is subjected to interference processing, and the time-series deformation data of various landslides are obtained using SBAS-InSAR technology. The specific operation process is as follows:

After the interference phase is removed, it can be assumed that the obtained interferogram does not contain residual terrain phase, atmospheric phase and noise phase. At this time, the jth interferogram in the azimuth-distance pixel coordinate system (x, r) can be expressed as:

In the formula, λ is the radar wavelength; d(t _B , x, r) and d(t _A , x, r) are respectively the accumulated deformation of the line of sight defense at time t _B and t _A corresponding to the reference time t ₀ , so d(t ₀ , x, r)≡0; d(t _i , x, r), (i=1, 2, 3...N (number of images)) is used to represent the deformation time series, and the corresponding phase is That is:

The M phase values corresponding to the deformation of each pixel point in the study area are represented by vectors:

will be calculated from the differential interferogram values are represented as vectors:

Wherein, main image IE=[IE ₁ , IE ₂ , IE ₂ , . . . , IE _M ], secondary image IS=[IS ₁ , IS ₂ , IS ₃ , . . . , IS _M ].

According to the above formula, the vector phase of the interferogram can be expressed in the form of matrix:

In the matrix, each row corresponds to each differential interferogram, and each column corresponds to a SAR image at a different time. In this matrix, the column of the main image and the secondary image is ±1, and the remaining columns are 0, such as That is, the matrix G of order [M×N] can be expressed as:

When a series of interference pairs are generated in the same small baseline subset, M≥N and N is the rank of G, the phase matrix can be obtained by the least square method:

For the phenomenon that it is impossible for the interference pairs of all images to be in the same baseline subset, the existing SARs are combined into a certain number of subsets according to the threshold, and then G ^T G is a reduced-rank matrix,

The effect of noise and coherence on the interfering pair can thus be reduced. Singular value decomposition (the core algorithm of SBAS-InSAR) is used for the phenomenon that the matrix is not a rank matrix, that is, the singular matrix is decomposed into:

G = USV ^T (12)

In the formula, V ^T is the velocity of the average phase, U is an orthogonal matrix, and S is a diagonal matrix. The phase can then be converted to an average phase velocity:

After obtaining the minimum norm solution of the velocity vector V, the integral is obtained to obtain the estimated value of the phase. A more accurate deformation phase vector δ can be obtained. Part of the time-series deformation rate diagram of this embodiment is shown in FIG. 6 , and the deformation criterion model of the landslide area is trained according to the deformation data of the landslide area.

The acquisition method of the time-series deformation rate criterion diagram is basically similar to that of the preliminary early warning criterion diagram. The deformation threshold is obtained by calculating the correlation coefficient between the time-series cumulative deformation and the number of landslide occurrences when the landslide occurs. According to the relationship diagram between the two, the corresponding deformation value is obtained when the number of landslides jumps. The deformation value corresponding to the first jump is the warning level deformation threshold, and the deformation value corresponding to the second jump is the alarm level deformation threshold. The secondary early warning criterion diagram of this embodiment is shown in FIG. 7 .

Step S6, using the critical deformation threshold to check the landslide initial warning area acquired in step S3 and the time-series deformation data acquired in step S4, so as to realize the secondary warning of the initial warning area. Specifically: comprehensively utilize the inversion deformation data before and after the landslide disaster occurs, and realize the accurate inversion of the time-series deformation variables and deformation rates in the landslide area. On this basis, step S5 establishes the critical deformation criterion diagram for the landslide occurrence. The critical deformation threshold of the landslide obtained in step S5 is used to check the landslide preliminary warning area and its time-series deformation data obtained in steps S3 and S4, and further determine the activity of the preliminary warning area according to the relationship between the time-series deformation data of the preliminary warning area and the deformation criterion map, and finally determine the type of the target area: attention level, warning level, and warning level.

Although the embodiment of the present invention has been described in conjunction with the accompanying drawings, the patent owner can make various deformations or modifications within the scope of the appended claims, as long as they do not exceed the protection scope described in the claims of the present invention, all should be within the protection scope of the present invention.

Claims (8)

1. A SAR landslide early warning method based on landslide deformation information and meteorological data, is characterized in that, comprises the following steps: S1, calculate the correlation coefficient between the number of times that the history landslide occurs and the effective rainfall before the landslide occurs, determine the rainfall threshold that the effective rainfall causes the landslide; S2. Correct the rainfall threshold obtained in step S1 according to the correlation coefficient between the time-series rainfall and the time-series deformation data; S3. Using the rainfall threshold value and regional rainfall information to carry out the initial warning of landslides, and obtain the initial warning area of landslides; S4. Obtain the time-series deformation data of the early warning area of the landslide; S5. According to the relationship between the number of landslide occurrences and the time-series cumulative deformation data, determine the effective deformation amount that causes the landslide, and calculate the critical deformation threshold of the landslide accordingly; S6. Check the landslide initial warning area obtained in step S3 and the time-series deformation data obtained in step S4 by using the critical deformation threshold, so as to realize the secondary warning of the landslide initial warning area. 2. the SAR landslide early warning method based on landslide deformation information and meteorological data according to claim 1, is characterized in that, the computing mode of the effective rainfall in described step S1 is: γ _n =γ ₀ +Tγ ₁ +T ² γ ₂ +T ³ γ ₃ +...+T ⁿ γ _n ; In the formula, γ _n represents the effective rainfall; γ ₀ represents the rainfall of the day; γ ₁ , γ ₂ , γ ₃ , ..., γ _n represent the rainfall before the current day; n is the number of days passed; 3. the SAR landslide early warning method based on landslide deformation information and meteorological data according to claim 1, it is characterized in that, described step S3 is specially: the judgment that effective rainfall is greater than rainfall threshold is the early warning area of landslide, the judgment that effective rainfall is less than or equal to effective rainfall threshold is stable region. 4. the SAR landslide early warning method based on landslide deformation information and meteorological data according to claim 3, is characterized in that, described step S4 is to obtain the time series deformation data of this area by carrying out time series InSAR processing to landslide early warning area. 5. the SAR landslide early warning method based on landslide deformation information and meteorological data according to claim 4, is characterized in that, described by carrying out time series InSAR processing to landslide early warning area to obtain the time series deformation data of this area specifically as: It is assumed that the jth interferometric phase map is expressed in the azimuth-distance pixel coordinate system (x, r) as: In the formula, λ is the wavelength of the radar; d(t _B , x, r) and d(t _A , x, r) are the accumulated deformation of the line of sight defense corresponding to the reference time t ₀ at time t _B and t _A respectively, and d(t ₀ , x, r)≡0; d(t _i , x, r), i=1, 2, 3...N is used to represent the deformation time series, and the corresponding phase is There are: /> The N phase value corresponding to the deformation of each pixel point in the study area is represented by a vector Represent the M values computed from the differential interferogram as a vector, where Wherein, main image IE=[IE ₁ , IE ₂ , IE ₃ ,..., IE _M ], secondary image IS=[IS ₁ , IS ₂ , IS ₃ ,..., IS _M ]; According to the above formula, the vector phase of the interferogram can be expressed in the form of matrix: In the matrix, each row corresponds to each differential interferogram, and each column corresponds to SAR images at different times. The columns of the main image and secondary images in the matrix are ±1, and the remaining columns are 0. like That is, the matrix G of order [M×N] can be expressed as When a series of interference pairs are generated in the same small baseline subset, M≥N and N is the rank of G, the phase matrix can be obtained by the least square method: For the phenomenon that it is impossible for the interference pairs of all images to be in the same baseline subset, the existing SAR is combined into a certain number of subsets according to the threshold, and then G ^T G is a matrix after rank reduction. If the matrix is not satisfied with the rank matrix, the singular value decomposition method is used to decompose the singular matrix into: G=USV ^T , where V ^T is the velocity of the average phase, U is the orthogonal matrix, and S is the diagonal matrix, and then the phase can be converted into the average phase velocity: Then obtain the minimum norm solution of the velocity vector V, and integrate to obtain the estimated value of the phase, which is the accurate deformation phase vector δ. 6. the SAR landslide early warning method based on landslide deformation information and meteorological data according to claim 1, it is characterized in that, described step S5 is specifically: obtain critical deformation threshold value by calculating the correlation coefficient between the time-series cumulative deformation amount and landslide occurrence number of times when landslide occurs, obtain the corresponding deformation amount when the number of times of landslide jumps according to the relationship diagram between the two, the deformation value corresponding to the first jump is the warning level deformation threshold value, and the deformation value corresponding to the second jump is the alarm level deformation threshold value. 7. the SAR landslide early warning method based on landslide deformation information and meteorological data according to claim 6, is characterized in that, obtains landslide area deformation criterion model according to the deformation data of historical landslide area and obtains early warning level deformation threshold and warning level deformation threshold according to this model. 8. the SAR landslide early warning method based on landslide deformation information and meteorological data according to claim 1, it is characterized in that, the secondary early warning of described landslide early warning area is specifically, according to the relationship between the time series deformation data of landslide early warning area and early warning level deformation threshold value, warning level deformation threshold value determines the active degree of landslide early warning area.