Drainage basin cascade reservoir dam group landslide monitoring method based on space-air-ground multi-source cooperation

Abstract

The invention discloses a drainage basin cascade reservoir dam group landslide monitoring method based on space-air-ground multi-source cooperation, and the method comprises the following steps: S1, obtaining drainage basin cascade reservoir dam group recognition data, and carrying out the preprocessing of the drainage basin cascade reservoir dam group recognition data; s2, calculating a surface deformation result, and obtaining a potential landslide mass of the drainage basin cascade reservoir dam group; s3, acquiring design data, and acquiring landslide mass classification; and S4, monitoring the landslide according to the landslide mass classification. According to the method, the drainage basin level, the water-power engineering level and the landslide mass level are analyzed by selecting data information based on space-air-ground multi-source data collaboration, different levels of landslide analysis monitoring methods are designed, and the landslide disaster monitoring and early warning capability of the drainage basin cascade reservoir dam group is effectively improved; the position, number, area and risk level of the landslide mass in the cascade reservoir dam group section of the drainage basin can be known more clearly from the level of the drainage basin, a detailed monitoring method, means and period are formulated through comprehensive analysis, research and judgment, and the landslide geological disaster monitoring and treatment capacity is improved.

Description

Technical Field

[0001] The present invention belongs to the technical field of landslide geological disaster identification and monitoring, and particularly relates to a landslide monitoring method for cascade reservoir dam groups in a basin based on multi-source collaboration of air, space and ground. Background Art

[0002] At present, China's geological disaster monitoring system has been basically formed, mainly adopting a monitoring method of integrating multi-source data of air, space and ground. Generally speaking, "air" refers to aerial remote sensing, "space" refers to space satellite remote sensing, and "ground" refers to the collection of ground-based remote sensing, field geological surveys, basic geological data and professional monitoring. Ground-based remote sensing usually uses devices such as lidar, panoramic cameras, and high-resolution satellite images to obtain data at specific locations on the earth's surface, and these data can be used for generating digital elevation models, three-dimensional terrain models, ground object classification, deformation monitoring, etc.

[0003] For example, the 2019 paper "Research on Dynamic Evaluation of Landslide Risk Based on Integrated Monitoring of Air, Space and Ground" carried out a research on the dynamic evaluation of landslide risk based on integrated monitoring of air, space and ground using 30m DEM data, Landsat 8 satellite data, geology, river network, basin, rainfall and other data. Taking Wanzhou District in the middle reaches of the Yangtze River as an example, it mainly analyzed and studied individual landslide monomers.

[0004] For example, the Chinese patent with the publication number CN118171918B discloses a method and system for landslide monitoring of air, space and ground based on multi-source data, which uses GNSS monitoring data, UAV aerial photography image data, ground real-time video data and environmental monitoring data to monitor landslides.

[0005] To sum up, the existing technologies mainly carry out research on the integrated monitoring technology of air, space and ground for local areas and monitoring means, lacking a monitoring method that systematically monitors different levels of cascade reservoirs, hydropower projects and individual landslide bodies from the perspective of cascade reservoir dam groups in a basin. Summary of the Invention

[0006] To solve the above technical problems, the present invention provides a landslide monitoring method for cascade reservoir dam groups in a basin based on multi-source collaboration of air, space and ground, effectively improving the monitoring and early warning ability of cascade reservoirs in the basin to landslide disasters.

[0007] The present invention is achieved through the following technical solutions.

[0008] A landslide monitoring method for cascade reservoir dam groups in a basin based on multi-source collaboration of air, space and ground provided by the present invention includes the following steps:

[0009] S1. Obtain the identification data of the cascade reservoir dam group in the basin, and preprocess the identification data of the cascade reservoir dam group. The identification data of the cascade reservoir dam group includes DEM data, InSAR data and high-resolution satellite data;

[0010] S2. Calculate the surface deformation result based on the identification data of the cascade reservoir-dam group in the basin, and obtain the potential landslide bodies of the cascade reservoir-dam group in the basin according to the surface deformation result;

[0011] S3. Obtain the design data, and obtain the landslide classification according to the potential landslide bodies of the cascade reservoir-dam group in the basin and the design data. The landslide classification includes landslide bodies with greater hazards, landslide bodies with threats, and landslide bodies with minor threats;

[0012] S4. Monitor the landslides according to the landslide classification.

[0013] Preferably, the InSAR data includes satellite position, velocity information, and SAR images;

[0014] The preprocessing of the identification data of the cascade reservoir-dam group in the basin includes the following steps:

[0015] Crop the DEM data;

[0016] Perform radiometric calibration, Doppler terrain correction, multi-look processing, multi-temporal registration, format conversion, and data cropping on the SAR images;

[0017] Perform radiometric calibration and geometric correction on the high-resolution satellite data.

[0018] Preferably, the S2. Calculate the surface deformation result based on the identification data of the cascade reservoir-dam group in the basin, and obtain the potential landslide bodies of the cascade reservoir-dam group in the basin according to the surface deformation result includes the following steps:

[0019] S21. Calculate the interference phase and coherence coefficient based on the InSAR data;

[0020] S22. Remove the flat-earth phase of the InSAR data, and remove the terrain phase of the InSAR data according to the DEM data;

[0021] S23. Perform differential interferometric processing on the InSAR data, and perform phase unwrapping on the InSAR data according to the coherence coefficient to obtain an unwrapping result map;

[0022] S24. Select ground control points according to the coherence coefficient and the unwrapping result map, and perform orbit refinement and re-flattening processing on the InSAR data according to the ground control points;

[0023] S25. Perform SBAS inversion according to the interference phase and the InSAR data to obtain an inversion result;

[0024] S26. Perform geocoding on the InSAR data according to the inversion result to obtain the surface deformation result. The surface deformation result includes the average surface deformation rate in a certain period of time in the area where the cascade reservoir-dam group is located;

[0025] S27. Set a deformation threshold, identify potential landslides based on high-resolution satellite data, the deformation threshold, and the surface deformation results, catalog the potential landslides, and obtain the potential landslide bodies of the cascade reservoir dams in the basin.

[0026] S28. Verify the surface deformation results.

[0027] Preferably, the step S24 of selecting ground control points according to the coherence coefficient and the unwrapping result map, and performing orbit refinement and re-flattening processing on the InSAR data according to the ground control points includes the following steps:

[0028] S241. Create ground control points according to the coherence coefficient.

[0029] S242. Select ground control points according to the unwrapping result map.

[0030] S243. Perform orbit refinement on the InSAR data according to the ground control points to obtain new orbit parameters.

[0031] S244. Re-flatten according to the new orbit parameters to obtain the re-flattening result, and update the InSAR data according to the re-flattening result.

[0032] Preferably, the step S26 of performing geocoding on the InSAR data according to the inversion result to obtain the surface deformation result includes the following steps:

[0033] S261. Convert the InSAR data from the spherical coordinate system to the northeast-up coordinate system according to the inversion result.

[0034] S262. Convert the InSAR data from the northeast-up coordinate system to the geocentric coordinate system.

[0035] Preferably, the step S4 of monitoring landslides according to the landslide body classification includes the following steps:

[0036] S41. Automatically monitor large-hazard landslide bodies with landslide GNSS to obtain monitoring data of large-hazard landslide bodies.

[0037] S42. Monitor potentially threatening landslide bodies by UAV aerial flight modeling to obtain monitoring data of potentially threatening landslide bodies.

[0038] S43. Obtain InSAR time series data, and obtain basin monitoring data according to the InSAR time series data.

[0039] Preferably, the step S41 of automatically monitoring large-hazard landslide bodies with landslide GNSS for large-hazard landslide bodies includes the following steps:

[0040] S411. Arrange monitoring stations according to large-hazard landslide bodies.

[0041] S412. Install the GNSS receiver at the monitoring site;

[0042] S413. Receive satellite signals through the GNSS receiver, collect GNSS data, and transmit the GNSS data to the data processing center;

[0043] S414. Process the GNSS data to obtain the monitoring data of the large - hazard landslide body.

[0044] Preferably, the step of S414. Processing the GNSS data to obtain the monitoring data of the large - hazard landslide body includes the following steps:

[0045] Convert the coordinate system of the GNSS data to obtain the coordinate conversion equation;

[0046] Calculate the conversion parameters according to the coordinate conversion equation to obtain the conversion equation;

[0047] Process the conversion equation through real - time kinematic differential to obtain the double - difference observations;

[0048] Calculate the positioning result according to the double - difference observations to obtain the monitoring data of the large - hazard landslide body.

[0049] Preferably, the step of S42. Using an unmanned aerial vehicle (UAV) for flight and modeling to monitor the threatening landslide body and obtaining the monitoring data of the threatening landslide body includes the following steps:

[0050] S421. Develop a UAV flight plan according to the scope of the threatening landslide body;

[0051] S422. Collect aerial photography data and image control point data according to the UAV flight plan;

[0052] S423. Conduct real - scene three - dimensional modeling based on the aerial photography data and image control point data to obtain a three - dimensional model;

[0053] S424. Obtain DOM and DSM from the three - dimensional model, and evaluate the accuracy of the three - dimensional model, DOM, and DSM;

[0054] S425. Extract the horizontal displacement according to the DOM;

[0055] S426. Extract the vertical displacement according to the horizontal displacement and DSM;

[0056] S427. Combine the horizontal displacement and vertical displacement to obtain the monitoring data of the threatening landslide body.

[0057] Preferably, the step of S425. Extracting the horizontal displacement according to the DOM includes the following steps:

[0058] Construct a Gaussian scale space based on the DOM to obtain the extreme points;

[0059] Locate the extreme points to obtain the key points;

[0060] Calculate the local gradient direction histogram based on key points;

[0061] Generate a feature description based on the key points and the local gradient direction histogram, and obtain the horizontal displacement;

[0062] The step S426 of extracting the vertical displacement according to the horizontal displacement and the DSM includes the following steps:

[0063] Obtain the elevation change according to the DSM;

[0064] Obtain the vertical displacement of the ground surface points according to the elevation change and the horizontal displacement;

[0065] Denoise the vertical displacement of the ground surface points to obtain the vertical displacement.

[0066] The beneficial effects of the present invention are as follows:

[0067] Select data based on the collaboration of multi-source data from space, air, and ground to analyze at the basin level, hydropower project level, and landslide level, and design monitoring methods for analyzing landslides at different levels, effectively improving the monitoring and early warning capabilities of cascade reservoir dams in the basin to respond to landslide disasters; it helps to more clearly understand the location, quantity, area, and risk level of landslides within the cascade reservoir dam section at the basin level, comprehensively analyze and judge to formulate detailed monitoring methods, means, and cycles, and improve the monitoring and governance capabilities of landslide geological disasters. BRIEF DESCRIPTION OF THE DRAWINGS

[0068] Figure 1 is the flowchart of the method provided by the embodiment of the present invention;

[0069] Figure 2 is the schematic diagram of DEM data provided by the embodiment of the present invention;

[0070] Figure 3 is the schematic diagram of InSAR data provided by the embodiment of the present invention;

[0071] Figure 4 is the schematic diagram of high-resolution satellite data provided by the embodiment of the present invention;

[0072] Figure 5 is the schematic diagram of the surface observation geometry by InSAR technology provided by the embodiment of the present invention;

[0073] Figure 6 is the schematic diagram of the surface deformation result provided by the embodiment of the present invention;

[0074] Figure 7 is the schematic diagram of the surface deformation result greater than the deformation threshold of 10 mm / a provided by the embodiment of the present invention;

[0075] Figure 8It is a schematic diagram of the geometric relationship between the LOS direction and the GNSS ENU direction provided by an embodiment of the present invention;

[0076] Figure 9 It is a schematic diagram of the surface deformation result of ascending orbit provided by an embodiment of the present invention;

[0077] Figure 10 It is a schematic diagram of extracting vertical displacement provided by an embodiment of the present invention. Specific embodiments

[0078] The technical solutions of the present invention will be further described below, but the scope of protection is not limited thereto.

[0079] As Figure 1 shown, a landslide monitoring method for cascade reservoir dams in a basin based on multi-source cooperation of space-air-ground includes the following steps:

[0080] S1. Obtain the identification data of the cascade reservoir dams in the basin, and preprocess the identification data of the cascade reservoir dams in the basin. The identification data of the cascade reservoir dams in the basin includes DEM data, InSAR data, and high-resolution satellite data;

[0081] The identification data of the cascade reservoir dams in the basin refers to the identification data of the area where the cascade reservoir dams are located, that is, the DEM data, InSAR data, and high-resolution satellite data of the area where the cascade reservoir dams are located. Among them, DEM is the digital elevation model, and InSAR is synthetic aperture radar interferometry.

[0082] The InSAR data includes satellite position, velocity information, and SAR images;

[0083] In this embodiment, landslide monitoring of the cascade reservoir dams in the experimental area is carried out. As Figures 2 to 4 shown, the identification data of the cascade reservoir dams in the basin obtained in this embodiment are: DEM data with a resolution of 30 meters of SRTM1 in the upper reaches of the Yellow River in the experimental area, that is, DEM data; C-band Sentine1-1A SAR images covering 79 + 102 scenes of ascending and descending orbits in the experimental area from 20210103 to 20240628, that is, InSAR data; high-resolution satellite data of Gaofen-1 and Gaofen-2 satellites, that is, high-resolution satellite data.

[0084] Among them, the basic principle of InSAR includes: When performing synthetic aperture radar interferometry, a synthetic aperture radar system carried by a satellite or an aircraft is used to obtain a pair of complex images of the same ground scene through simultaneous observation with two antennas or two nearly parallel observations. Due to the geometric relationship between the target and the positions of the two antennas, a phase difference is generated on the complex image, forming an interferogram. The interferogram contains precise information about the difference between the points along the slant range and the positions of the two antennas. Based on the phase difference information of the complex radar image, using the geometric relationship between the sensor altitude, radar wavelength, beam viewing direction, and antenna baseline distance, three-dimensional information of the ground target terrain is extracted through imaging processing, interferometric data processing, and geometric transformation, etc. The biggest difference between InSAR interferometry technology and other remote sensing means is that it is a phase-based measurement. Therefore, the data used is complex data, that is, single-look complex data.

[0085] As Figure 5 shown Figure 5 is a vertical cross-section of the satellite's observation of the ground, used to explain the principle of InSAR technology. Among them, S1 and S2 respectively represent the positions of the two antennas, the distance between them is represented by the baseline distance B, the angle between the baseline and the horizontal direction is α, and the baseline can be decomposed into a component B / / along the slant range direction and a component B ⊥ perpendicular to the slant range direction. H1 represents the altitude of satellite S1, H2 represents the viewing angle of satellite S2, R i represents the slant range from the satellite to a point P on the ground, where i = 1, 2, corresponding to satellite S1 and satellite S2 respectively, and the elevation of the point on the ground is represented by Z. Additionally, let R2 = R1 + ΔR.

[0086] Generally, assuming that the random phase contributions caused by the scattering characteristics are basically the same in the two acquired images, the phase of the interferogram is only related to the path difference of the radar signal propagation.

[0087]

[0088] Note: represents the phase, and λ represents the radar wavelength.

[0089] According to the cosine theorem, the following relational expressions can be calculated:

[0090]

[0091] Since the distance R i between the ground object and the SAR satellite is much larger than the spatial baseline B, and R1 is much larger than ΔR, the formula 1 and formula 1 can be combined to obtain the following formula:

[0092]

[0093] If θ1, α, λ, Given B and H1, the vertical height Z from the ground object to the reference ellipsoid can be derived through formula as follows:

[0094] Z = H1 - R1cosθ1

[0095]

[0096] During the actual solution process, the phase is converted to elevation as follows:

[0097]

[0098] The measurement accuracy of InSAR surface elevation is related to the vertical baseline B ⊥ , the distance r between the SAR sensor and the ground object, the radar wavelength λ, and the incident angle θ1.

[0099] Preprocessing the data for identifying the cascade reservoir dam group in the basin facilitates subsequent processing of the data for identifying the cascade reservoir dam group in the basin and realizes obtaining potential landslide bodies of the cascade reservoir dam group in the basin. The preprocessing specifically includes the following steps:

[0100] Clip the DEM data;

[0101] Perform radiometric calibration, Doppler terrain correction, multi-looking processing, multi-temporal registration, format conversion, and data clipping on the SAR image;

[0102] Perform radiometric calibration and geometric correction on the high-resolution satellite data.

[0103] S2. Calculate the surface deformation result based on the data for identifying the cascade reservoir dam group in the basin, and obtain the potential landslide bodies of the cascade reservoir dam group in the basin according to the surface deformation result;

[0104] Step S2 is to analyze the cascade reservoir dam group in the basin at the basin level. Among them, the potential landslide bodies of the cascade reservoir dam group refer to the potential landslide bodies in the area where the cascade reservoir dam group is located.

[0105] This step specifically includes the following steps:

[0106] S21. Calculate the interference phase and coherence coefficient based on the InSAR data;

[0107] The interference phase is the phase difference of corresponding pixels in two SAR images, which contains information on various factors such as surface displacement, earth curvature, terrain, atmospheric influence, and temporal variation of scatterers. The conjugate multiplication method usually refers to the method of multiplying a complex number by its conjugate complex number in complex number operations. For any complex number z = a + bi, where a and b are real numbers and i is the imaginary unit, its conjugate complex number is z - = a - bi. The rules of the conjugate multiplication method are as follows:

[0108] z·z- =(a + bi)(a - bi)=a 2 -bi 2 =a 2 +b 2 =|z| 2

[0109] The SAR coherence coefficient is an important parameter for evaluating the interference quality of two SAR images, which reflects the coherence between the two images. The coherence coefficient γ is a value between 0 and 1. The larger the γ value, the better the coherence between the two images and the smaller the noise.

[0110] Calculation formula of the coherence coefficient:

[0111] In practical applications, since the SAR sensor cannot obtain multiple observations at the same time, it is usually assumed that the equilibrium processes S1(t), S2(t) and S1(t)S2 * (t) are ergodic. The spatial average of multiple pixels in the estimation window is used to replace the time average, so as to obtain the estimator γ of the coherence magnitude γ ∧ :

[0112]

[0113] S1(L) and S2(L) represent the complex-valued backscattering coefficients of two SAR images at pixel point L.

[0114] S22, removing the flat-earth phase from the InSAR data, and removing the topographic phase from the InSAR data according to the DEM data;

[0115] The purpose of removing the flat-earth phase is to eliminate the phase changes caused by slant range variations and atmospheric disturbances, etc., and retain the phase changes related to terrain undulations.

[0116] The specific steps in this embodiment are as follows:

[0117] Using the relationship formula between phase and terrain elevation for elevation inversion: After obtaining the flat interferogram, combined with satellite orbit parameters, such as the height H of the satellite platform, the angle θ between the orbit and the earth's surface, the radar wavelength λ, and the baseline length B⊥, the digital elevation model, i.e., DEM, can be inverted. In this process, special attention needs to be paid to the problem of elevation ambiguity, and a prior DEM may be required for auxiliary correction.

[0118] Methods for eliminating the flat-earth effect: Calculating the flat-earth effect based on orbit parameters and the geodetic longitude and latitude of the center point of the imaging area; Calculating the flat-earth effect according to the image energy; Calculating the flat-earth effect by measuring the dominant interference fringe frequencies in the range and azimuth directions.

[0119] In this embodiment, the topographic phase is removed using DEM data. By subtracting the known DEM, a differential interferogram is obtained, thereby removing the influence of the topographic phase.

[0120] S23. Perform differential interferometric processing on InSAR data, perform phase unwrapping on the InSAR data according to the coherence coefficient, and obtain the unwrapping result map;

[0121] The unwrapping result map is used as a reference to select ground control points.

[0122] Differential interferometric processing is a technique for monitoring surface deformation, and time series analysis is a method for estimating the surface deformation rate using multiple interferograms. This method can identify and remove systematic errors, such as atmospheric delay and orbital errors.

[0123] In this embodiment, phase unwrapping is performed by minimum cost flow. Minimum cost flow is a commonly used unwrapping method, especially suitable for cases where large areas have low coherence or other limiting growth factors make unwrapping difficult. This method uses a square grid, considers all pixels on the image, and masks the pixels with coherence less than the threshold. The algorithm packages for the minimum cost flow problem include the Bellman-Ford algorithm, the SPFA algorithm, the Dijkstra algorithm, etc. These algorithms find the minimum cost path through different path relaxation techniques and perform flow augmentation until the minimum cost flow is found.

[0124] S24. Select ground control points according to the coherence coefficient and the unwrapping result map, and perform orbit refinement and re-flattening processing on the InSAR data according to the ground control points;

[0125] Orbit refinement means that in InSAR data processing, ground control points, abbreviated as GCPs, are used to redefine the baseline parameters to correct orbit errors. When the orbit parameters are not accurate enough, it will affect the conversion from the interferometric phase to the terrain height. The purpose of orbit refinement is to reduce the large fringes or orbit residuals on the interferogram caused by orbit errors. Through the artificially added GCP points, orbit recalculation can be performed to optimize the orbit parameters.

[0126] Re-flattening is a step performed after orbit refinement. The purpose is to re-modify the orbit parameters in the header file of the unwrapped image to calculate the phase offset, such as obtaining the absolute phase value. This step can reduce the orbit residuals on the interferogram and improve the accuracy of interferometric measurement.

[0127] The said S24. Select ground control points according to the coherence coefficient and the unwrapping result map, and perform orbit refinement and re-flattening processing on the InSAR data according to the ground control points includes the following steps:

[0128] S241. Create ground control points according to the coherence coefficient;

[0129] In the interferogram, that is, in the InSAR data, points with stable phase are selected as GCPs according to the coherence coefficient. These points usually have high coherence and no phase change, that is, flat points.

[0130] S242. Select ground control points according to the unwrapping result map;

[0131] In this embodiment, the filtered interferogram and the unwrapping result map are selected as references to assist in the selection of GCPs. Points with high coherence, no area with unremoved topographic phase, and no deformed area are selected as GCPs in the interferogram.

[0132] S243. Refine the orbit of the InSAR data according to the ground control points to obtain new orbit parameters;

[0133] S244. Re - flatten according to the new orbit parameters to obtain the re - flattened result, and update the InSAR data according to the re - flattened result.

[0134] S25. Perform SBAS inversion according to the interferometric phase and the InSAR data to obtain the inversion result;

[0135] The specific process of SBAS inversion is as follows:

[0136] Assume that within a period of time t0,..., t N N + 1 SAR images covering the same area are obtained. After registration of the master image and the slave image, and differential interferometry is performed by selecting appropriate temporal baselines and spatial baselines to obtain M differential interferograms. Assume that M is an odd number, and the relationship between M and N is shown in the following formula:

[0137]

[0138] Any t i The differential phase at time relative to the initial t0 time is unknown. The differential phase and the interferometric phase (k = 1,…,i…M) can be expressed as:

[0139]

[0140] Then, after removing the topographic phase, the interferometric phase of any pixel (x, r) in the k - th (k = 1,…,i…,M) differential interferogram is:

[0141]

[0142] λ represents the central wavelength of the radar signal emitted by the SAR satellite, d(t B ,x,r) and d(t A ,x,t) respectively represent t Aand t B The LOS deformation at time relative to the initial time t0. That is, the following relational expression can be deduced:

[0143]

[0144] d(t0, x, r) = 0

[0145] Assume that the main image sequence IE = [IE1…IE M and the auxiliary image sequence IS = [IS1…IS M are arranged in chronological order during data processing, and

[0146]

[0147] Then the differential interference phases of M differential interferograms can form the following observation equations:

[0148]

[0149] The above equation is a system of M equations with N unknowns, then:

[0150]

[0151] In the above equation, A represents an M*N-dimensional matrix. Let Then A can be expressed as:

[0152]

[0153] When all the generated interferograms are in a subset, then M≥N, the rank of matrix A is N. According to the above conditions, the estimated value of can be obtained

[0154]

[0155] In reality, the situation where all interference pairs are in the same subset is not common. Therefore, we need to face the situation where interference pairs belong to different subsets. When the interference pairs are divided into different subsets, the above equation shows rank deficiency, where A T A is a singular matrix. Assume that when M≥N, the interference pairs are divided into L subsets, then the rank deficiency number is N - L + 1, which will produce infinitely many solutions. Therefore, singular value decomposition is used to solve the minimum norm solution of matrix A.

[0156] A = USV T

[0157] In the above equation, U and V represent M×M and N×N orthogonal matrices respectively, and S represents an M×N matrix:

[0158]

[0159] In the above formula, r = N - L + 1, D = diag(σ1σ2...σ r ), σ i (i = 1, 2,..., r), representing the singular values of matrix A. Therefore, the parameter The solution of the least squares method with the minimum norm is expressed as:

[0160]

[0161] A + = VS + U +

[0162] That is

[0163]

[0164] In the above formula, u i and v i respectively represent the column vectors of U and V, A + , S + , U + respectively represent the generalized inverses of matrices A, S, and U.

[0165] The singular value decomposition has a minimum norm constraint on the obtained deformation. Therefore, this method must require the differential phase to approach 0, that is, it is required to be based on solving the minimum norm of the differential interferometric phase signal. However, it may introduce large discontinuities in the obtained cumulative deformation, resulting in physically meaningless results. Therefore, the average phase rate v between SAR images acquired at adjacent times is used to replace the unknown, expressed as:

[0166]

[0167] That is

[0168]

[0169] In the above formula, B is an M×N matrix, which can be expressed as:

[0170]

[0171] The SVD decomposition is used for matrix B to solve the cumulative deformation amount at the corresponding moment.

[0172] That is, B can be decomposed into the product of three matrices, the left singular vector matrix U, the singular matrix Σ, and the right singular vector matrix V T .

[0173] B = UΣV T

[0174] U is an m×m orthogonal matrix, and its column vectors are the left singular vectors of B.

[0175] Σ is an m×n diagonal matrix, and the non - negative real numbers on its diagonal are called singular values and are arranged in descending order, while other elements are zero.

[0176] V T is an n×n orthogonal matrix, and its column vectors are the right singular vectors of B.

[0177] S26. Geocode the InSAR data according to the inversion result to obtain the surface deformation result, and the surface deformation result includes the average surface deformation rate within a certain period of time in the area where the cascade reservoir dam group is located.

[0178] The S26. Geocoding the InSAR data according to the inversion result to obtain the surface deformation result includes the following steps:

[0179] Geocoding is the process of converting SAR data from the radar coordinate system (slant - range coordinate system) to the geographic coordinate system.

[0180] S261. Convert the InSAR data from the spherical coordinate system to the east - north - up coordinate system according to the inversion result.

[0181] The spherical coordinate system is the radar slant - range coordinate system.

[0182] The east - north - up coordinate system is also called the local - tangent plane coordinate system, abbreviated as the ENU coordinate system. The east - north - up coordinate system is a Cartesian coordinate system established with the observation station as the origin center and is a local coordinate system. The X - axis points to the local east direction, the Y - axis points to the local geographic north direction, and the Z - axis points vertically upward.

[0183] The spherical coordinate system (R, θ, ), the distance R, the distance from the target P to the origin, range [0, +∞]; the azimuth angle θ: the angle turned from the positive Z - axis to the projection of the line connecting the target P and the origin in the XOY plane in the counter - clockwise direction from the X - axis, range in radians [0, 2π], and the elevation angle the angle between the line connecting the target P and the origin and the positive Z - axis, range in radians [0, π].

[0184] The conversion between the coordinates (R, θ, ) of the target P in the radar - polar coordinate system and the coordinates (e, n, u) of the target P in the local - tangent plane coordinate system (ENU) with the antenna of the radar station as the origin is:

[0185]

[0186] S262. Convert the InSAR data from the east - north - up coordinate system to the geocentric coordinate system.

[0187] The coordinates of target P in the local geodetic coordinate system with the antenna of the radar station as the origin are (e, n, u). The longitude and latitude of the radar station are (λ, ), and its coordinates in the geocentric coordinate system are (X0, Y0, Z0). The coordinates of the target in the geocentric coordinate system are (X, Y, Z), that is:

[0188]

[0189] After geocoding, the surface deformation result is obtained. The surface deformation result obtained in this embodiment is as Figure 6 shown. Specifically, the surface deformation result includes: the average deformation rate of the surface in the experimental area from 2021 to 2024.

[0190] S27. Set the deformation threshold, identify potential landslides according to the high-resolution satellite data, deformation threshold and surface deformation result, catalog the potential landslides, and obtain the potential landslide bodies of the basin cascade reservoir dam group;

[0191] In this embodiment, the annual average deformation rate of 10 mm / a is used as the deformation threshold. According to the surface deformation result, combined with the terrain, high-resolution images and historical landslide data, the obtained image is as Figure 7 shown.

[0192] Catalog the potential landslide bodies of the basin cascade reservoir dam group according to the threshold, and mainly record parameters such as the number, name, slope, aspect, area and maximum deformation rate of the landslide body.

[0193] S28. Verify the surface deformation result.

[0194] After obtaining the time-series deformation of the experimental area through InSAR technology for surface deformation monitoring, in order to evaluate the reliability and accuracy of the deformation result, accuracy assessment is required. It can be analyzed by comparing different orbits and different types of SAR data, or by comparing InSAR results with ground monitoring data such as leveling data and GNSS data. However, since the ENU-direction deformation in GNSS measurement refers to the deformation amounts on the earth's surface in the eastward, northward and vertical directions, while the LOS-direction deformation obtained by InSAR measurement refers to the deformation amount along the radar line of sight. There is no simple direct conversion relationship between these two deformation amounts, and some calculations and conversions are required.

[0195] As Figure 8 shown, according to the geometric schematic diagram, the geometric relationship between the LOS-direction unit vector and the eastward, northward and vertical directions can be deduced, and its geometric relationship formula is as follows:

[0196]

[0197] Among them, since the surface deformation result of this embodiment is the descending orbit surface deformation result, it can also be verified with the ascending orbit surface deformation result as shown in Figure 9 shown.

[0198] S3. Obtain design data, and obtain the landslide classification based on the potential landslides of the cascade reservoir dam group in the basin and the design data. The landslide classification includes landslides with greater hazards, threatening landslides, and landslides with minor threats.

[0199] Step S3 is to analyze the cascade reservoir dam group in the basin at the hydropower project level.

[0200] The design data includes the management scope of the hydropower project, topography, geomorphology, and hub layout. Through the overlay analysis of the design data, the threat level of the landslide can be obtained, so as to classify the potential landslides of the cascade reservoir dam group in the basin and obtain the landslide classification. In this embodiment, for landslides with a slope greater than 25°, a vegetation coverage rate lower than 30%, and the distance between buildings and the landslide less than 200 m, which can cause damage to hydraulic structures and affect the operation of the hydropower project and directly endanger people's lives and property, they are defined as landslides with greater hazards; for landslides that can affect the external and on-site traffic of the hydropower project, they are defined as threatening landslides; and for those far from the management scope of the hydropower project, they are defined as landslides with minor threats.

[0201] S4. Monitor the landslides according to the landslide classification.

[0202] Step S4 is to analyze the cascade reservoir dam group in the basin at the landslide level and design corresponding monitoring methods to realize the method of analyzing landslide monitoring at different levels.

[0203] The step S4 of monitoring the landslides according to the landslide classification includes the following steps:

[0204] S41. Automatically monitor the landslides with greater hazards by using GNSS for the landslides with greater hazards and obtain the monitoring data of the landslides with greater hazards.

[0205] The step S41 of automatically monitoring the landslides with greater hazards by using GNSS for the landslides with greater hazards includes the following steps:

[0206] S411. Arrange monitoring stations according to the landslides with greater hazards.

[0207] According to factors such as the scope, shape, and geological conditions of the landslide, reasonably arrange monitoring stations on the landslide and the surrounding stable areas. Generally, key positions of the landslide, such as the top, shoulder, and toe of the slope, are mainly arranged, and the stations should be able to effectively reflect the deformation of the landslide. At the same time, set up a reference station, which is usually located at a stable location outside the influence range of the landslide and is used to provide accurate reference coordinates for the monitoring stations.

[0208] S412. Install the GNSS receiver at the monitoring site;

[0209] Install relevant equipment such as GNSS receivers at the monitoring site. Ensure that the equipment is firmly installed, the antenna can accurately receive satellite signals, and the equipment should be well protected to avoid damage from the natural environment and human factors. Connect the communication line so that the monitoring station can transmit data to the data processing center.

[0210] S413. Receive satellite signals through the GNSS receiver, collect GNSS data, and transmit the GNSS data to the data processing center;

[0211] The GNSS receiver automatically receives satellite signals and continuously and real - time collects data at the set sampling frequency. The collected data is mainly information such as the distance between the receiver and the satellite, which is used for subsequent coordinate calculation.

[0212] Use wired or wireless communication methods to transmit the collected data from the monitoring site to the data processing center or server. Ensure the integrity and accuracy of the data during the transmission process.

[0213] S414. Process the GNSS data to obtain monitoring data of large - hazard landslides.

[0214] At the data processing center, process the transmitted data, and calculate the three - dimensional coordinates of the monitoring site through professional software. Compare the coordinates at different times, analyze parameters such as the change amount and change rate of the coordinates, and use this to judge the displacement of the landslide. For example, if the elevation of a certain monitoring point continuously decreases and there is a continuous offset in the horizontal direction within a certain period of time, it may indicate that the landslide is active.

[0215] The step of S414. Process the GNSS data to obtain monitoring data of large - hazard landslides includes the following steps:

[0216] Convert the coordinate system of the GNSS data to obtain the coordinate conversion equation;

[0217] It is necessary to implement the conversion of the coordinate system of the GNSS collected data, and convert the WGS84 coordinates to the CGCS2000 plane coordinates. The conversion formula is as follows:

[0218]

[0219] Note: ΔX, ΔY, and ΔZ are the translation parameters in three directions, ωx, ωy, and ωz are the rotation parameters in three directions, and m represents the scaling ratio.

[0220] Calculate the conversion parameters according to the coordinate conversion equation to obtain the conversion equation;

[0221] As can be seen from the above formula, it is necessary to list the coordinate transformation equations using the coordinates of at least three or more common points in the survey area, and use the least squares principle to solve seven transformation parameters, and then obtain the transformation equation.

[0222] Principle of least squares: For a given set of data points (x i , y i )(i = 1, 2,..., n), assume that we want to use a function, such as a linear function y = ax + b, to fit these data points. Then for each data point (x i , y i ), the value predicted by the function is yy i . When the fitting function is y = ax + b, yy i = ax i + b. The true value y i of these points and the predicted value yy i have an error e i = y i - yy i . The goal of the least squares method is to find a set of parameters, such as a and b in the above linear function, to minimize the sum of the squares of all errors .

[0223] By processing the transformation equation through real-time kinematic differential, double-difference observations are obtained;

[0224] Real-time kinematic differential, abbreviated as RTK technology, is mainly based on the differential processing of carrier phase observations. Let the carrier phase observation value of satellite i observed by the reference station receiver be and the carrier phase observation value of satellite i observed by the user station receiver be Then the single-difference observation value can be expressed as:

[0225] Further, the double-difference observation value can be obtained. Let satellite j be another satellite, then the double-difference observation value between satellite i and satellite j at the reference station and the user station is:

[0226]

[0227] Calculate the positioning result according to the double-difference observation value to obtain the monitoring data of the landslide body with greater hazard.

[0228] In an ideal situation, ignoring errors such as measurement noise, the double-difference observation value is related to the baseline vector XX = (x, y, z) (i.e., the coordinate difference between the reference station and the user station), the satellite coordinates (x i , y i , z i ) and ((x j , y j , zj ) There are the following relationships with the carrier wavelength λ, etc.:

[0229]

[0230] In practical applications, the RTK system sends the observed data and coordinate information to the user station through the reference station. The user station uses this information combined with its own observed data to perform calculations according to the above mathematical relationships, thereby obtaining a high-precision positioning result.

[0231] By identifying and judging the monitoring data of the more dangerous landslide body, when the displacement exceeds a certain value, a geological disaster may occur to the landslide body, and the system automatically issues a warning message. Relevant personnel make decisions such as evacuating personnel and taking reinforcement measures based on the warning message and the data analysis results.

[0232] S42. The UAV flight and modeling monitoring targets the landslide body posing a threat, and obtains the monitoring data of the landslide body posing a threat.

[0233] The S42, where the UAV flight and modeling monitoring targets the landslide body posing a threat and obtains the monitoring data of the landslide body posing a threat, includes the following steps:

[0234] S421. Formulate a UAV flight plan according to the scope of the landslide body posing a threat;

[0235] For the landslide body posing a threat, according to the spatial layout, formulate a UAV flight plan, use a five-lens camera, set the flight route, and the cycle needs to be specified. Determine the scope of the acquisition area, resolution requirements, etc., and use professional mission planning software to design the flight route to ensure complete coverage of the area and sufficient overlap between adjacent flight routes and adjacent photos. Generally, the forward overlap is 60%-80%, and the side overlap is 30%-40%. Arrange the ground control points GCP in the survey area in advance.

[0236] S422. Collect aerial photography data and image control point data according to the UAV flight plan;

[0237] Fly under suitable weather conditions. Before taking off, check the status of the UAV and the camera to ensure the normal operation of the equipment. During the flight, take photos automatically according to the preset flight route and parameters, and pay attention to the stability of the flight attitude and altitude of the UAV. After the flight, check whether the collected data is complete. mainly check whether the number of photos meets the expectations and whether the photo quality is qualified. Re-collect the data that does not meet the requirements in a timely manner.

[0238] Collect within the survey area according to the already arranged image control points. The image control points are mainly stable and long-preserved marks in the survey area.

[0239] S423. Perform real - scene three - dimensional modeling based on aerial photography data and image control point data to obtain a three - dimensional model;

[0240] Real - scene three - dimensional modeling mainly includes processes such as feature extraction and matching, camera positioning and pose estimation, point cloud reconstruction, and model optimization and texture mapping.

[0241] Among them, the judgment method for feature matching is:

[0242] Let the image patches I1(x, y) and I2(x, y) around two feature points, and their normalized cross - correlation coefficient is:

[0243]

[0244] Among them, I1' and I2' are the means of the corresponding image patches respectively. By setting a threshold to judge whether they match, if r is greater than the threshold, it is considered a match.

[0245] The process of sparse point cloud reconstruction is:

[0246] Given the rotation matrices R1, R2 and translation vectors t1, t2 of two cameras, the internal parameter matrices K1, K2 of the two cameras, and the image points p1, p2 corresponding to the point P on the images of the two cameras, and the image point coordinates are normalized, then the coordinates X of the spatial point P can be calculated in the following way:

[0247] x1 = K1 -1 p1

[0248] x2 = K2 -1 p2

[0249] X = λx1

[0250]

[0251] The process of model optimization is:

[0252] The objective function is generally constructed based on the reprojection error. Let the spatial point X i , the camera pose ε i , which is uniformly represented including rotation and translation parameters, etc. The observed image point of the point X i under the camera ε i is x ij , and the reprojection error e ij can be expressed as: e ij = x ij - π(K, ε j , X i ), where π is the projection function, which projects the spatial point onto the image plane according to the camera internal parameters and pose to obtain the image point. The overall objective function is to minimize the sum of the squares of all reprojection errors min∑ i,j e ij2 , the camera pose and spatial point coordinates are optimized by solving a non - linear optimization algorithm.

[0253] S424. Obtain DOM and DSM from the three - dimensional model, and evaluate the accuracy of the three - dimensional model, DOM and DSM;

[0254] Among them, DOM represents the digital orthophoto map, and DSM represents the digital surface model.

[0255] In this embodiment, the process of obtaining DOM is as follows: First, aerial triangulation is carried out to obtain the exterior orientation elements of the image, and then the image is orthorectified according to the collinearity equation. During the rectification process, the influence of terrain undulation needs to be considered, and DEM data is used for height difference correction. For example, for each image point, the corresponding ground elevation is obtained according to the DEM data, the position difference between the image point and the ground projection point is calculated, and the position of the image point is adjusted to make the image an orthographic projection image. Finally, image mosaicking and cropping are carried out to obtain DOM.

[0256] Based on the collinearity equation in photogrammetry, assuming the image point coordinates are (x, y), the corresponding ground point coordinates are (X, Y, Z), the interior orientation elements of the camera (principal point coordinates (x0, y0), focal length f) and the exterior orientation elements (camera position (X s , Y s , Z s ), attitude angles ω, k), the collinearity equation is:

[0257]

[0258] In this embodiment, the process of obtaining DSM is as follows: Assuming the coordinates of a point in the three - dimensional model are (x, y, z), for numerous discrete points in the whole area, DSM can be expressed as a function Z = f(x, y) of the plane coordinates (x, y), where Z is the elevation value, that is, the height information, reflected by the ground surface and above - ground objects at the position of (x, y). When actually extracting from the real - scene three - dimensional model, it is often through analyzing and interpolating geometric elements such as each triangular patch in the model, if based on triangular mesh modeling or voxels, such as voxel - based modeling, etc., to determine the Z value corresponding to different (x, y) positions. For example, the weighted average interpolation method based on neighboring points is used to construct the relationship between Z and (x, y). Assuming there are points (x1, y1, z1), (x2, y2, z2) …… (x n , y n , z n ) these neighboring sampling points, the Z0 at a point to be solved (x0, y0), that is, the elevation of DSM at this point, can be expressed by an expression similar to the following. Taking simple distance - weighted average as an example:

[0259]

[0260] d i is the distance from the point to be determined (x0, y0) to the sampling point (x i , y i ), such as the Euclidean distance:

[0261]

[0262] In this embodiment, the process of evaluating the accuracy of the 3D model, DOM, and DSM is as follows:

[0263] The GCP dataset collected on-site is divided into two parts. One part is used for real-scene 3D reconstruction, and the other part is used as checkpoints to evaluate the accuracy of the 3D model, DOM, and DSM. The typical root mean square error is selected as the accuracy evaluation index. Identify each checkpoint from the produced DOM and extract the planar coordinates, and obtain the elevation of the checkpoint from the DSM based on the planar coordinates. The calculation formula is as follows:

[0264]

[0265] where n is the number of checkpoints, X DOM , Y DOM , Z DSM are the 3D coordinate data obtained from the DOM and DSM, X GNSS , Y GNSS , Z GNSS are the 3D coordinates of the checkpoints actually measured using RTK-CORS, and RMSE x , RMSE y , RMSE z are the root mean square errors in the X, Y, and Z directions.

[0266] S425. Extract the horizontal displacement according to the DOM;

[0267] It should be noted that here the DOM is a multi-temporal DOM, that is, multiple periods of DOM are required to extract the horizontal displacement according to the DOM.

[0268] The step S425 of extracting the horizontal displacement according to the DOM includes the following steps:

[0269] Construct a Gaussian scale space based on the DOM to obtain the extreme points;

[0270] In this embodiment, when constructing the Gaussian scale space, the Gaussian function is expressed as where (x, y) is the image pixel position and σ is the scale space factor. The scale space L(x, y, σ) of an image I(x, y) is obtained by convolving it with a Gaussian function: L(x, y, σ) = G(x, y, σ) × I(x, y). To detect the extreme points in the scale space, each pixel is compared with its 8 neighboring pixels in the current scale and the corresponding 18 pixels in the adjacent upper and lower scales, a total of 26 points. If the point is an extreme point among these points, it may be a feature point.

[0271] Locate the extreme points and obtain the key points;

[0272] For the detected extreme points, precise localization is required. The Taylor expansion is used to fit the scale space function D(x) (D(x) is the difference function of L(x, y, σ) in the scale space) where x = (x, y, σ) T . By taking the derivative and setting it to 0, the precise position tt of the key point is obtained:

[0273] At the same time, points with low contrast and edge response points are removed by calculating the Hessian matrix of this point. The determinant of the Hessian matrix Judges whether it is an edge response according to the ratio of its eigenvalues.

[0274] Calculate the local gradient direction histogram based on the key points;

[0275] For each key point, calculate its local gradient direction histogram to determine the main direction. In the neighborhood centered on the key point, the gradient magnitude m(x, y) and direction θ(x, y) of the pixel point (x, y) are calculated as follows:

[0276]

[0277] The neighborhood is divided into sub-intervals of 360° / 8 = 45°. The cumulative magnitude of the gradient direction in each sub-interval is counted, and the direction with the largest magnitude is the main direction of the key point.

[0278] Generate a feature description based on the key points and the local gradient direction histogram, and obtain the horizontal displacement.

[0279] Taking the key point as the center, a 16×16 neighborhood window is taken and divided into 4×4 sub-regions.

[0280] In each sub-region, the gradient histogram in 8 directions is calculated. In this way, each sub-region has 8 values, and a total of 4×4×8 = 128-dimensional feature description sub-vector is generated for the 4×4 sub-regions to describe the features of this key point.

[0281] After generating the feature description, it is necessary to remove the outliers in the feature description and then obtain the horizontal displacement.

[0282] Use a sliding window to obtain the horizontal displacement vectors in the neighborhood of the central pixel. Apply a Gaussian function to each displacement vector according to the modulus d of the vector difference between the displacement vector of the central pixel and the displacement vector of the central pixel i Assign the weight w i

[0283]

[0284] where N is the number of displacement vectors in the neighborhood of the central pixel, and σ g is the standard deviation of the Gaussian function. After weighting, calculate the displacement vector density ρ at the central pixel cv :

[0285]

[0286] S426. Extract the vertical displacement according to the horizontal displacement and DSM;

[0287] The step S426 of extracting the vertical displacement according to the horizontal displacement and DSM includes the following steps:

[0288] Obtain the elevation change according to the DSM;

[0289] Obtain the vertical displacement of the ground point according to the elevation change and the horizontal displacement;

[0290] As Figure 10 shown, due to the existence of horizontal movement, the elevation change and vertical displacement obtained from the multi-temporal DSM for the ground deformation points are different, and the greater the horizontal displacement and the original terrain slope, the greater the difference between the two. Before the ground movement, there is a ground point (x0, y0) in the image L1. After the ground movement, the horizontal displacement (Δx, Δy) at the ground point (x1, y1) is obtained by using the horizontal displacement detection result. It can be seen that the horizontal position of the point (x1, y1) after the ground movement in the image L2 after the ground movement is (x1 + Δx, y1 + Δy). The elevation H1(x1, y1) at the point (x1, y1) before the movement and the elevation H2(x1 + Δx, y1 + Δy) at the origin (x1, y1) after the movement (the point (x1 + Δx, y1 + Δy) in the image L2) are obtained through the DSM data before and after the ground movement. Finally, the vertical displacement ΔH at the ground point (x, y) can be obtained:

[0291] ΔH = H2(x1 + Δx, y1 + Δy) - H1(x1, y1)

[0292] Denoise the vertical displacement of the ground point to obtain the vertical displacement.

[0293] Adopt a denoising method based on the vertical displacement sequence. The elevation of ground point X before the surface movement is H1, and the elevation at X monitored in the i-th period is H i , where i ∈ [1, …, n], n is the total number of ground monitoring times, and the corresponding total number of monitoring periods is According to the elevation sequence {G1, G2, …, H n} at X, the vertical displacement Δ n-1 h i :

[0294] Δ n-1 h i = H n-i+j - H j , where j ∈ [1, 2, …, n - i]

[0295]

[0296] is the matrix composed of all period vertical displacements obtained from the elevation sequence at X. Theoretically, in the unaffected area, Δ n-i h j is 0, while in the subsidence area, Δ n-i h j , is less than 0. However, due to measurement errors, in the stable area, Δ n-i h j usually fluctuates within a certain range around 0, while in the subsidence area, Δ n-i h j may be greater than 0. Therefore, for Δ n-i h j , a tolerance value ε H is set, that is, when Δ n-i h j > ε H , it is an abnormal vertical displacement. If there is an abnormal displacement in one of all periods, it is determined that point X is an abnormal displacement point. ε H is determined according to the mean and standard deviation of the vertical displacements in the stable area of each period :

[0297]

[0298] S427. Obtain the monitoring data of the threatening landslide body by combining horizontal displacement and vertical displacement.

[0299] Similarly, by identifying and judging the monitoring data of the threatening landslide body, corresponding decisions are made, such as evacuating people and taking reinforcement measures.

[0300] S43. Obtain InSAR time series data and obtain the basin monitoring data according to the InSAR time series data.

[0301] InSAR time-series data refers to InSAR data in a recent time period. Through InSAR time-series data, recent surface deformation results, i.e., basin monitoring data, can be obtained, thereby enabling InSAR monitoring at the basin level. In this embodiment, the InSAR data is C-band Sentine1-1A SAR images of ascending and descending orbits covering 79 + 102 scenes in the experimental area from January 3, 2021 to June 28, 2024. Then the InSAR time-series data is SAR images covering the experimental area from July to August 2024, and basin monitoring data from July to August is obtained, thus realizing InSAR monitoring of the basin level from July to August in this embodiment.

[0302] Among them, the specific process of obtaining basin monitoring data through InSAR time-series data is the same as that in steps S21 to S26.

[0303] The potential landslide bodies of the cascade reservoir dam group in the basin identified and classified by the present invention can accurately determine the monitoring objects, effectively reduce the monitoring cost of large-scale areas, and implement the responsibilities of the monitoring objects. The monitoring data of landslide bodies with greater hazards and the monitoring data of landslide bodies posing threats obtained by monitoring are highly accurate. Combined with basin monitoring data, the landslide bodies of the cascade reservoir dam group in the basin can be monitored at different scales, further improving the monitoring accuracy and effectively enhancing the monitoring and early warning ability of the cascade reservoir dam group in the basin to respond to landslide disasters.

[0304] The present invention selects data materials based on the collaboration of multi-source data from space, air, and ground to analyze the basin level, hydropower project level, and landslide body level, designs monitoring methods for analyzing landslides at different levels, and effectively enhances the monitoring and early warning ability of the cascade reservoir dam group in the basin to respond to landslide disasters; it helps to more clearly know the location, quantity, area, and risk level of the landslide bodies in the cascade reservoir dam group section from the basin level, comprehensively analyze and judge to formulate detailed monitoring methods, means, and cycles, and improve the monitoring and treatment ability of landslide geological disasters.

Claims

1. A landslide monitoring method for cascade reservoirs and dams in a watershed based on multi-source coordination of air, ground and space, characterized in that: The following steps are involved: S1. Acquire the identification data of the cascade reservoir-dam group in the watershed, and pre-process the identification data of the cascade reservoir-dam group in the watershed, wherein the identification data of the cascade reservoir-dam group in the watershed includes InSAR data, DEM data and high-resolution satellite data; S2. Calculate the surface deformation results based on the identification data of the cascade reservoir-dam group in the basin, and obtain the potential landslide body of the cascade reservoir-dam group in the basin based on the surface deformation results; S3, obtaining design data, and obtaining landslide classification according to the potential landslides of the cascade reservoir dam group in the basin and the design data, wherein the landslide classification includes landslides with greater hazards, landslides with threats, and landslides with minor threats; S4. Monitor the landslide according to the classification of landslide bodies.

2. The method for monitoring landslides in a cascade reservoir-dam group in a watershed as claimed in claim 1, characterized in that: The InSAR data includes satellite position, velocity information and SAR images; The preprocessing of the identification data of the cascade reservoir and dam group in the watershed comprises the following steps: Clip the DEM data; Perform radiometric calibration, Doppler terrain correction, multi-view processing, multi-temporal registration, format conversion and data clipping on SAR images; Perform radiometric calibration and geometric correction on Gaofen satellite data.

3. The method for monitoring landslides in a cascade reservoir-dam group in a watershed as claimed in claim 1, characterized in that: The step S2, calculating the surface deformation results according to the identification data of the cascade reservoir dam group in the watershed, and obtaining the potential landslide body of the cascade reservoir dam group in the watershed according to the surface deformation results, comprises the following steps: S21, calculating the interference phase and coherence coefficient according to the InSAR data; S22, removing the flat ground phase of the InSAR data, and removing the terrain phase of the InSAR data according to the DEM data; S23, performing differential interferometry processing on the InSAR data, performing phase unwrapping on the InSAR data according to the coherence coefficient, and obtaining an unwrapping result graph; S24, selecting ground control points according to the coherence coefficient and the unwrapping result map, and performing track refinement and re-flattening processing on the InSAR data according to the ground control points; S25, performing SBAS inversion according to the interferometric phase and InSAR data to obtain an inversion result; S26, geocoding the InSAR data according to the inversion results to obtain surface deformation results, wherein the surface deformation results include an average surface deformation rate within a certain period of time in the area where the cascade reservoir-dam group is located; S27, setting deformation thresholds, identifying potential landslides based on high-resolution satellite data, deformation thresholds and surface deformation results, cataloguing potential landslides, and obtaining potential landslide bodies of cascade reservoir dam groups in the basin; S28. Verify the surface deformation results.

4. The method for monitoring landslides in a cascade reservoir-dam group in a watershed as claimed in claim 3, characterized in that: The step S24, selecting ground control points according to the coherence coefficient and the unwrapping result graph, and performing track refinement and re-flattening processing on the InSAR data according to the ground control points, comprises the following steps: S241, creating ground control points according to the coherence coefficient; S242, selecting a ground control point according to the unwrapping result map; S243, refining the orbit of the InSAR data according to the ground control points to obtain new orbit parameters; S244. Re-level the data according to the new orbit parameters, obtain the re-leveling result, and update the InSAR data according to the re-leveling result.

5. The method for monitoring landslides in a cascade reservoir-dam group in a watershed as claimed in claim 3, characterized in that: The S26, geocoding the InSAR data according to the inversion result to obtain the surface deformation result, comprises the following steps: S261, converting the InSAR data from the spherical coordinate system to the northeast celestial coordinate system according to the inversion results; S262. Convert the InSAR data from the northeast celestial coordinate system to the geocentric coordinate system.

6. The method for monitoring landslides in a cascade reservoir-dam group in a watershed according to claim 1, characterized in that: The S4, monitoring the landslide according to the landslide body classification, comprises the following steps: S41. Automatic monitoring of landslides with high hazard by GNSS, and acquisition of monitoring data of landslides with high hazard; S42, UAV flight modeling to monitor threatening landslides and obtain monitoring data of threatening landslides; S43. Acquire InSAR time series data, and acquire watershed monitoring data according to the InSAR time series data.

7. The method for monitoring landslides in a cascade reservoir-dam group in a watershed as claimed in claim 6, characterized in that: The S41, landslide GNSS automated monitoring of landslides with greater hazards, the landslides with greater hazards comprising the following steps: S411. Set up monitoring stations according to the landslide bodies with greater hazards; S412, install GNSS receivers at monitoring sites; S413, receiving satellite signals through a GNSS receiver, collecting GNSS data, and transmitting the GNSS data to a data processing center; S414. Process GNSS data to obtain monitoring data of landslides with greater hazards.

8. The method for monitoring landslides in a cascade reservoir-dam group in a watershed as claimed in claim 7, characterized in that: The S414, processing GNSS data to obtain monitoring data of landslides with greater hazards, comprises the following steps: Convert the GNSS data coordinate system and obtain the coordinate conversion equation; Calculate the conversion parameters according to the coordinate conversion equation to obtain the conversion equation; The double-difference observations are obtained by converting the equations through real-time dynamic difference processing; The positioning results are calculated based on the double-difference observations to obtain monitoring data of landslides with greater hazards.

9. The method for monitoring landslides in a cascade reservoir-dam group in a watershed as claimed in claim 6, characterized in that: The S42, UAV flight modeling monitoring of threatening landslide bodies, and obtaining monitoring data of threatening landslide bodies include the following steps: S421. Develop a UAV flight plan based on the scope of the threatening landslide body; S422, collecting aerial photography data and image control point data according to the UAV flight plan; S423, performing real scene three-dimensional modeling according to the aerial photography data and the image control point data to obtain a three-dimensional model; S424, obtaining DOM and DSM according to the three-dimensional model, and evaluating the accuracy of the three-dimensional model, DOM and DSM; S425, extracting horizontal displacement according to DOM; S426, extracting vertical displacement according to horizontal displacement and DSM; S427. Combine horizontal displacement and vertical displacement to obtain monitoring data of threatening landslides.

10. The method for monitoring landslides in a cascade reservoir-dam group in a watershed according to claim 9, characterized in that: The step S425 of extracting the horizontal displacement according to the DOM comprises the following steps: Construct Gaussian scale space according to DOM and obtain extreme points; Locate extreme points and obtain key points; Calculate the local gradient direction histogram based on the key points; Generate feature description based on key points and local gradient direction histogram to obtain horizontal displacement; The step S426 of extracting the vertical displacement according to the horizontal displacement and the DSM comprises the following steps: Obtain elevation changes based on DSM; Obtain vertical displacement of surface points based on elevation change and horizontal displacement; The vertical displacement of the surface points is denoised to obtain the vertical displacement.

Images