Geological catastrophe monitoring method for diversion canal of non-pressure diversion type hydropower station

Abstract

The invention discloses a geological catastrophe monitoring method for a diversion canal of a non-pressure diversion type hydropower station, and relates to the technical field of geological catastrophe monitoring. Determining a region of interest in the range of the target diversion canal according to the change of the slope ratio of the terrain along the way of the target diversion canal; expanding the direction of the window in the region of interest along the trend of the diversion canal, determining homogeneous pixels and performing phase optimization, and obtaining candidate distributed target pixels from the homogeneous pixels according to a phase consistency discrimination method; fusing the permanent scatterer candidate points in the time sequence synthetic aperture radar image with the candidate distributed target pixels, performing phase compensation according to the self-statistical point locations, and calculating the deformation value of each position in the range of the target diversion canal to determine a suspected geological catastrophe region; and determining the geological catastrophe area of the target diversion canal according to the high-resolution optical and thermal infrared images of the suspected geological catastrophe area and the deformation value. The method improves the recognition accuracy of the geological catastrophe along the diversion canal of the non-pressure diversion type hydropower station.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of geological disaster monitoring, and particularly relates to a geological disaster monitoring method for a water diversion channel of a non-pressure diversion type hydropower station. BACKGROUND

[0002] At present, the water diversion channel of the non-pressure diversion type hydropower station is an important part of the hydropower station. Different from the main project of the power station, the water diversion channel has the characteristics of long distribution distance, complex geological conditions along the line and long-term operation. Since the hydropower station needs water level potential energy for power generation, it is generally located in an area with large topographic changes. The corresponding area has complex and variable terrain, large geomorphologic span, complex geological background conditions and relatively intense tectonic activity, which may cause the water diversion channel to be potentially threatened by various geological disasters.

[0003] Traditional geological disaster monitoring and identification of the water diversion channel mainly relies on artificial field investigation and fixed-point monitoring of instruments and equipment. Not only a large amount of manpower and material resources are consumed, but also the method is limited by complex field terrain and weather conditions. Some areas are difficult to reach, and it is difficult to comprehensively and accurately collect data along the water diversion channel, resulting in low data collection efficiency and limited coverage. It is even more difficult to achieve high-frequency dynamic monitoring. In addition, the traditional method of relying solely on on-site manual identification is excessively dependent on the identification experience of relevant personnel, and is prone to misjudgment and omission. Moreover, the initial characteristics of the geological disasters along the water diversion channel are not significant, and it is also difficult for on-site personnel to interpret, so that the geological disaster area along the line cannot be found in time, and the best treatment stage of the geological disaster is missed. High-resolution optical remote sensing is also an important method for monitoring geological disasters along the water diversion channel. With the characteristics of high resolution and rapid large-scale coverage, it is mainly used for identification and monitoring evaluation of significant geological disaster-related features along the water diversion channel, and then for analysis and interpretation of geological disaster information along the water diversion channel to serve the safe and stable operation of the water diversion channel. However, it is difficult to image in cloudy and rainy days and at night, and it is difficult to ensure the timeliness and regularity of monitoring along the water diversion channel. Moreover, the water diversion channel is often located in an area with complex and variable terrain, and the atmospheric convection activity is relatively intense. The cloud and fog coverage is more serious than in plain areas. In addition, optical remote sensing cannot detect small deformations on the ground, and its detection capability for the initial characteristics of geological disasters is poor. It is generally used for loss evaluation after the occurrence of geological disasters, and it is difficult to meet the needs of detection and monitoring of the initial signals of geological disasters along the water diversion channel. The change of vegetation coverage along the water diversion channel caused by seasonal changes also seriously affects the reliability of optical remote sensing in identifying geological disaster features.

[0004] In summary, the traditional artificial field verification and high-resolution optical remote sensing survey methods face various challenges in monitoring geological disasters along the water diversion channel, and it is difficult to meet the needs of accurate identification of geological disasters along the water diversion channel. SUMMARY

[0005] Therefore, it is necessary to provide a geological disaster monitoring method for a diversion channel of a non-pressure diversion type hydropower station to improve the identification accuracy of geological disasters along the diversion channel of the non-pressure diversion type hydropower station.

[0006] The application adopts the following technical scheme: The application provides a geological disaster monitoring method for a diversion channel of a non-pressure diversion type hydropower station, comprising: Obtaining time-series synthetic aperture radar images in a range where the target diversion channel is located, and registering all the synthetic aperture radar images to the same space-time reference; According to the change of the terrain slope ratio along the target diversion channel, a region of interest is determined; in the region of interest, the window direction is extended along the direction of the diversion channel, the homogenous pixels are determined from the registered synthetic aperture radar images, and the phase is optimized; according to the phase consistency discrimination method, the candidate distributed target pixels are determined from the phase-optimized homogenous pixels; Obtaining the candidate points of the permanent scatterers in the time-series synthetic aperture radar images, and fusing the candidate points of the permanent scatterers with the candidate distributed target pixels to obtain target monitoring points; According to the distribution of the target monitoring points, a self-statistical point is selected, the phase stability of each target monitoring point is calculated according to the self-statistical point, and based on the phase stability, the atmospheric phase of the target monitoring points is compensated by using the weighted least square method; According to the phase of the target monitoring points after compensation, a suspected geological disaster area in the range where the target diversion channel is located is determined; According to the high-resolution optical and thermal infrared images of the suspected geological disaster area, the temperature data of the suspected geological disaster area are determined, and according to the temperature data and the deformation value of the suspected geological disaster area, the disaster anomaly coefficient of the suspected geological disaster area is determined; According to the disaster anomaly coefficient of the suspected geological disaster area, the geological disaster area of the target diversion channel is determined.

[0007] Optionally, according to the change of the terrain slope ratio along the target diversion channel, the region of interest is determined, comprising: According to the terrain slope ratio along the target diversion channel, the width along the line of the target diversion channel is determined; According to the width along the line of the target diversion channel, the region of interest is determined.

[0008] Optionally, the calculation formula of the width along the line of the target diversion channel is: ; wherein, the width along the line is W, the terrain slope ratio is S, and the basic width is B.

[0009] Optionally, the phase optimization of the homogeneous pixels comprises: constructing a sample complex coherence matrix according to all the homogeneous pixels; optimizing the time-series phase sequence of each homogeneous pixel by solving the maximum likelihood estimation phase based on the eigenvalue decomposition method according to the sample complex coherence matrix.

[0010] Optionally, the sample complex coherence matrix is constructed in the following manner: ; wherein, is the sample complex coherence matrix, is is the amplitude normalized result vector of the complex observation vector of the amplitude synthetic aperture radar image at a distributed pixel, is is the conjugate transpose of is a set of homogeneous pixels of the distributed target, is the corresponding number of homogeneous pixels.

[0011] Optionally, the eigenvalue decomposition method is as follows: ; wherein, is a scale factor, is the time-series phase sequence optimized by the EMI method, and the solution of the eigenvalue decomposition method is is the eigenvector corresponding to the minimum eigenvalue after the eigenvalue decomposition, is the time-series phase sequence to be optimized.

[0012] Optionally, the suspected geological disaster area in the range of the target diversion channel is determined according to the compensated phase of the target monitoring point, comprising: calculating the deformation value of each target monitoring point according to the compensated phase of the target monitoring point; calculating the local Moran coefficient at each target monitoring point according to the deformation value at each target monitoring point; calculating the deformation threshold based on the target monitoring point with the local Moran coefficient greater than 0, and determining the target monitoring point with the deformation value greater than the deformation threshold as the suspected geological disaster area.

[0013] Optionally, the calculation formula of the disaster anomaly coefficient of the suspected geological disaster area is as follows: ; wherein, denotes the deformation value at the i th pixel in the deformation statistical window of the suspected geological disaster area; denotes the area corresponding to the i th pixel in the deformation statistical window. ​​This indicates the number of pixels within the deformation statistics window; The temperature statistics window indicating the suspected geological disaster area is the first one. The area corresponding to each pixel; Indicates the number of times the temperature statistics window is filled. Temperature of each pixel; This represents the average temperature of the pixels within the temperature statistics window. This indicates the number of pixels within the temperature statistics window.

[0014] Optionally, based on the catastrophic anomaly coefficient of the suspected geological hazard area, the geological hazard area of ​​the target water diversion canal is determined, including: Areas suspected of geological disasters with a disaster anomaly coefficient greater than the safety coefficient threshold are identified as geological disaster areas for the target water diversion canal.

[0015] Optionally, the method further includes: After identifying the geological disaster area of ​​the target water diversion channel, a geological disaster distribution map is drawn for on-site disaster verification.

[0016] This invention provides a geological disaster monitoring device for the diversion canal of an unpressurized water diversion hydropower station, comprising: The acquisition module is used to acquire time-series synthetic aperture radar images within the range of the target water diversion canal and register all synthetic aperture radar images to the same spatiotemporal reference. The cropping module is used to determine the region of interest based on the change in the slope ratio of the terrain along the target irrigation canal; within the region of interest, it expands the direction of its window along the irrigation canal, identifies homogeneous pixels from the registered synthetic aperture radar image and performs phase optimization, and determines candidate distributed target pixels from the phase-optimized homogeneous pixels according to the phase consistency discrimination method. The compensation module is used to acquire candidate permanent scatterer points in temporal synthetic aperture radar images, fuse the candidate permanent scatterer points with candidate distributed target pixels to obtain target monitoring points; select self-statistical points according to the distribution of target monitoring points, calculate the phase stability of each target monitoring point according to the self-statistical points, and use the weighted least squares method to compensate for the atmospheric phase of the target monitoring points based on the phase stability. The calculation module is used to calculate the deformation value of each location within the target water diversion canal based on the compensated phase of the target monitoring point, and to determine the suspected geological disaster area within the target water diversion canal based on the deformation value of each location. The first determining module is used to determine the temperature data of the suspected geological disaster area based on the high-resolution optical and thermal infrared images of the suspected geological disaster area, and to determine the disaster anomaly coefficient of the suspected geological disaster area based on the temperature data and deformation value of the suspected geological disaster area. The second determining module is configured to determine a geological disaster area of the target diversion channel according to a disaster anomaly coefficient of the suspected geological disaster area.

[0017] The application provides a computer readable storage medium, which stores a computer program, and the computer program is executed by a processor to realize the geological disaster monitoring method of the non-pressure diversion type hydropower station diversion channel.

[0018] The application provides a computer device, which comprises a memory, a processor and a computer program stored in the memory and capable of running on the processor, and the processor realizes the geological disaster monitoring method of the non-pressure diversion type hydropower station diversion channel when executing the program.

[0019] The application adopts the above at least one technical scheme to achieve the following beneficial effects: In the application, the method of shearing data with a fixed width along the length direction of the traditional diversion channel is improved according to the linear characteristics of the distribution of the diversion channel and the difference in topography and geomorphology along the line, the data filtering region of interest along the line of the diversion channel is selected based on the topographic slope ratio, the omission of remote geological disasters along the line of the diversion channel caused by the fixed width is avoided, meanwhile, the complex atmospheric correction method for the linear diversion channel is specifically proposed according to the atmospheric delay problem faced by the small deformation extraction technology of the time series InSAR along the line of the diversion channel, the phase stability discrimination method at the corresponding position is constructed through the sampling points randomly distributed along the line of the diversion channel, the phase stability is obtained, the basic assumption that the phase of the geological disaster is high frequency and the atmospheric phase is low frequency is used, the above phase stability is used, the atmospheric phase is compensated by using the weighted least square estimation, and the accuracy of the small deformation monitoring is improved; on the basis of accurately obtaining the surface subsidence along the line of the diversion channel, the geological disaster index along the line of the diversion channel is constructed by combining the high-resolution unmanned aerial vehicle obtained thermal infrared surface temperature data and the sensitivity of the seepage induced by the deformation of the diversion channel, the distribution of the geological disaster area can be more accurately indicated, and the identification accuracy of the geological disaster along the line of the non-pressure diversion type hydropower station diversion channel is improved. BRIEF DESCRIPTION OF DRAWINGS

[0020] The drawings described herein are used to provide further understanding of the application, and form a part of the application.

[0021] Figure 1 A geological disaster monitoring method flow chart of a non-pressure diversion type hydropower station diversion channel is provided. Figure 2 A topographic complexity adaptive width along the line of the diversion channel is provided. Figure 3A schematic diagram of spatial scale difference of atmospheric delay phase and geological disaster related phase is provided for the present application. Figure 4 A schematic diagram of calculation and application of geological disaster index along the diversion channel is provided for the present application. Figure 5 Another flowchart of the geological disaster monitoring method of the diversion channel of the non-pressure diversion type hydropower station is provided for the present application. DETAILED DESCRIPTION

[0022] In order to make the purpose, technical scheme and advantages of the present application clearer, the technical scheme of the present application will be described clearly and completely below in combination with specific embodiments of the present application and corresponding drawings. Obviously, the described embodiments are only some of the embodiments of the present application, not all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor fall within the scope of protection of the present application.

[0023] How to realize early identification and rapid monitoring of the early characteristics of geological disasters along the diversion channel of the hydropower station by fusing multi-source effective remote sensing data and various technical methods is a problem to be solved by the present application. Therefore, the present application provides a geological disaster monitoring method for the diversion channel of the non-pressure diversion type hydropower station.

[0024] In view of the problems faced by traditional artificial field monitoring and the shortcomings of high-resolution optics, the present application focuses on accurate detection of early characteristics such as deformation and water seepage related to geological disasters along the diversion channel. The present application makes full use of the characteristics of all-weather and all-day of space-borne synthetic aperture radar (SAR) image and the characteristics of unmanned aerial vehicle thermal infrared remote sensing not dependent on solar radiation. Through the time-series InSAR deformation extraction technology optimized for the spatial distribution form of the diversion channel, on the basis of obtaining the micro-deformation of the ground surface in the corresponding area, combined with the ground temperature anomaly information extracted by the unmanned aerial vehicle thermal infrared, the geological disaster index along the diversion channel is constructed, and the corresponding mapping is completed, which is used to indicate the distribution and type of geological disasters along the diversion channel of the hydropower station, so as to realize timely detection and identification of potential geological disasters along the diversion channel.

[0025] Unlike conventional time-series InSAR technology, the present application is aimed at the linear characteristics presented by the distribution of the diversion channel and the differences in the topography along the line, improves the method of fixed width shearing data along the length of the diversion channel of the traditional, and proposes a data selection method along the line of the diversion channel based on the complexity of the terrain, which not only ensures the data processing efficiency, but also avoids the omission of remote geological disasters along the line of the diversion channel caused by the fixed width; at the same time, in view of the atmospheric delay problem faced by the small deformation extraction technology of time-series InSAR along the diversion channel, a complex atmospheric correction method for linear diversion channels is proposed, a phase stability discrimination method is constructed at the corresponding position through the sampling points randomly distributed along the diversion channel, the phase stability is obtained, and by means of the basic assumption of high frequency of geological disaster related phase and low frequency of atmospheric phase, the above phase stability is used to compensate the corresponding atmospheric phase by using weighted least squares estimation, and the accuracy of small deformation monitoring is improved; on the basis of accurately obtaining the ground subsidence along the diversion channel, combined with the high-resolution unmanned aerial vehicle obtained thermal infrared surface temperature data, by means of its sensitivity to seepage induced by diversion channel deformation, a geological disaster index along the diversion channel is constructed, which can more accurately indicate the distribution of geological disaster area, and the acquisition of corresponding spaceborne SAR and unmanned aerial vehicle thermal infrared remote sensing image data has good repeatability and regularity, which can support the long-term monitoring of geological disasters along the diversion channel, so as to complete the mapping of the distribution of geological disasters along the diversion channel, and then serve the on-site verification of the diversion channel geological disasters, and lay a good data foundation for the rapid and accurate acquisition of diversion channel geological disasters.

[0026] The technical solutions provided by the embodiments of the present application will be described in detail below with reference to the drawings.

[0027] Figure 1 It is a flowchart of a geological disaster monitoring method for a non-pressure diversion type hydropower station diversion channel in the present application, which specifically includes the following steps: S101, acquiring time-series synthetic aperture radar images in the range where the target diversion channel is located, and registering all synthetic aperture radar images to the same space-time reference.

[0028] The time-series synthetic aperture radar images are acquired based on the position information of the on-orbit SAR satellite and the research area; the time-series synthetic aperture radar images include L or C band time-series synthetic aperture radar images, and the time-series spaceborne synthetic aperture radar images include M amplitude synthetic aperture radar images.

[0029] According to the time baseline, the spatial baseline and the spectrum center of the synthetic aperture radar image, based on the basic principle of minimum mean of the spatial baseline and the time baseline and minimum difference of the spectrum center, a super master image is selected from all time series of the synthetic aperture radar image. Generally, the synthetic aperture radar image obtained at the middle time period of the observation period is taken as the basis, and then the synthetic aperture radar image obtained at the imaging time which makes the spatial baseline and the spectrum center difference minimum is selected as the master image. Meanwhile, the influence of extreme weather such as heavy rain and snow should be considered, and the synthetic aperture radar image under the influence should not be taken as the master image. Finally, the images corresponding to the other time series of the spaceborne synthetic aperture radar image are registered to the master image.

[0030] S102, according to the change of the terrain slope ratio along the target diversion channel, a region of interest is determined; in the region of interest, the direction of the window is extended along the diversion channel, the homogeneous pixels are determined from the registered synthetic aperture radar image, and the phase optimization is performed, and the candidate distributed target pixels are determined from the phase-optimized homogeneous pixels according to the phase consistency discrimination method.

[0031] According to the change of the terrain slope ratio along the target diversion channel, the region of interest is determined, including: according to the terrain slope ratio along the target diversion channel, the width along the line of the target diversion channel is determined; and according to the width along the line of the target diversion channel, the region of interest is determined.

[0032] Specifically, according to the position distribution of the target diversion channel, the width along the line is adaptively set according to the terrain complexity (see formula (1)), and with the change of the terrain slope ratio along the diversion channel, the width is narrow in the flat area and wide in the steep terrain area, so that the surrounding terrain and its possible influence can be better compatible, and the adaptive ability of the terrain is possessed, and the limitation of the insufficient adaptability of the traditional fixed width to the surrounding complex terrain is broken; the calculation formula of the width along the line of the target diversion channel is: (1); wherein, is the width along the line, with the unit of km; is the terrain slope ratio, is the basic width, which is generally taken as 1, that is, the terrain change within the initial 1km width is obtained.

[0033] In the actual scene, the diversion channel has a certain width range along its extension direction. Taking the target diversion channel itself as the core, the region covered by the calculated width along the line is defined as the region of interest which needs to be focused on and analyzed.

[0034] For example, assuming that the width along the line of the target diversion channel is 5 meters, then on the image, taking the center line of the diversion channel as the reference, extending a certain distance (such as 2.5 meters) to both sides, the belt-shaped region formed is the region of interest. For example, Figure 2As shown, Figure 2 Adaptive width diagram of diversion channel along the terrain complexity.

[0035] Among the selection of homogeneous pixels, the traditional method is a square window, but the present application considers the linear characteristics of the diversion channel, and takes its main direction as a constraint. In the determined region of interest, the window is expanded in the direction of the target diversion channel, for example, a rectangle or a strip-shaped distribution shape is expanded along the target diversion channel, and then HTCI (Homogeneous Test Confidence Interval) method is used to select homogeneous pixels: (2); Among them, is a probability operator, and are the time series amplitude mean of the reference sample and the sample to be tested, respectively, and are the quantile points of the gamma distribution at the confidence level , is the number of samples.

[0036] In one embodiment, the phase of the homogeneous pixel is optimized, including: constructing a sample complex coherence matrix according to all homogeneous pixels; according to the sample complex coherence matrix, the maximum likelihood estimation phase is solved based on the eigenvalue decomposition method, and the time series phase sequence of each homogeneous pixel is optimized.

[0037] Phase optimization and distributed target selection. Based on the homogeneous pixel screening result, the sample covariance matrix and the complex coherence matrix are calculated, and the statistical characteristics of the distributed target can be described by the homogeneous pixel to construct the sample complex coherence matrix (CCM). The construction method of the sample complex coherence matrix is as follows:

[0038] (3); Among them, is the sample complex coherence matrix, is the complex observation vector of the amplitude synthetic aperture radar image at a certain distributed pixel after amplitude normalization, is the conjugate transpose of is the homogeneous pixel set of the distributed target, is the corresponding number of homogeneous pixels.

[0039] The maximum likelihood estimation phase is solved by using the eigenvalue decomposition method, and the eigenvalue decomposition method is: (4); wherein, is a scale factor, is a time series phase sequence optimized by the EMI method, is a time series phase sequence to be optimized. The solution of formula (4) is the eigenvector corresponding to the minimum eigenvalue after eigenvalue decomposition.

[0040] For each homogeneous pixel, according to the sample complex coherence matrix and the original time series phase sequence of the homogeneous pixel , solve formula (4) to obtain the optimized time series phase sequence of the homogeneous pixel.

[0041] According to the optimized time series phase sequence of all homogeneous pixels, the candidate homogeneous pixels are screened by the phase consistency discrimination method, that is, on this basis, the quality of phase optimization is described by using time coherence, and then high-quality homogeneous pixel points that meet the requirements, that is, candidate distributed target pixels, are screened. The corresponding method is to define the phase optimization time coherence , the higher the value, the better the consistency of the phase before and after optimization, and the higher the optimization quality. It can be expressed by the following formula: (5); wherein, and are the differential interferometric phases before and after optimization, respectively, , and are the phase values of the first and aperture synthetic aperture radar image, , and are the phase values of the first and aperture synthetic aperture radar image, , by setting a specific threshold for time coherence, use formula (5) to calculate the time coherence of each homogeneous pixel, and by comparing with the time coherence threshold, the homogeneous pixels higher than the threshold are considered as reliable homogeneous pixels, that is, candidate distributed target pixels.

[0042] S103, obtain the permanent scatter candidate points in the time series synthetic aperture radar image, and fuse the permanent scatter candidate points into candidate distributed target pixels, to obtain target monitoring points; according to the distribution of the target monitoring points, select a self-statistical point, calculate the phase stability of each target monitoring point according to the self-statistical point, and based on the phase stability, use the weighted least squares method to compensate the atmospheric phase of the target monitoring points.

[0043] The permanent scatterer (PS) candidate points are selected from the time series synthetic aperture radar image by using a gray scale discrete index, the permanent scatterer candidate points and candidate distributed scatterer (DS) pixels are fused through a set operation to obtain multiple observation points, and the stable observation points are further optimized and screened based on phase consistency to obtain target monitoring points.

[0044] The phase of any one interference figure can be linearly represented by the phases of other two related interference figures, as shown in the following formula: (6) Wherein, represents the differential interference phase generated by the SAR image with the serial number and after phase optimization; and respectively represent the differential interference phases of the main image and the image , the image after phase optimization; represents a phase winding operator.

[0045] For any observation point, according to formula (6), if the two sides of the formula are equal, or the difference between the two sides of the formula is within a certain range, it indicates that the observation point is a stable observation point, and then becomes a target monitoring point.

[0046] According to the point self-statistic phase stability, the atmospheric phase estimation and compensation are carried out. According to the general geometric characteristics of the water diversion channel, the channel presents the characteristics of being narrow and long, and the width-length ratio is extremely low. The atmosphere changes greatly along the water channel direction (10-50 kilometers), and changes little in the direction perpendicular to the water channel (1-5 kilometers). The direction perpendicular to the water channel is low frequency, and the direction along the water channel is relatively high frequency (the deformation characteristics of the channel along the line are very high frequency). Considering the local and discrete characteristics of geological disasters, the statistical characteristics of the data itself and the polynomial fitting method are used to estimate and remove the atmospheric phase. Specifically, first, according to the distribution of the target monitoring point positions obtained based on formulas (4) and (5) in the region of interest, a part of the positions, i.e., the self-statistic positions, are selected from the above effective target monitoring point positions by means of random sampling, and the atmospheric phase is estimated and removed at the corresponding self-statistic positions by using the polynomial fitting method. 17 The phase stability W at the corresponding point is statistically analyzed using a window of 17 or larger, as shown in Formula (7). Then, the low-frequency atmospheric phase change is fitted with statistical points scattered in the region and a polynomial. In order to improve the accuracy of the estimation, the stability is weighted, as shown in Formula (8), to avoid the influence of phase noise on atmospheric noise estimation, and at the same time, to reduce the interference to high-frequency signals related to geological disasters.

[0047] Phase stability based on its own statistical characteristics mainly depends on the neighborhood phase variance. (7); in, Indicates the self-statistical point i Phase stability, Based on self-statistical points i Sampling points in the neighborhood within the center window, This refers to the number of sampling points in the window, which is typically 17×17. The average phase within the corresponding window. As a stable term, it is generally set as to , For the window number j The phase of each sampling point.

[0048] The phase stability of each self-statistical point can be calculated based on formula (7).

[0049] The polynomial coefficients weighted by the phase stability of the self-statistical points are estimated using the least squares method as follows: (8); in, and These are the location information and phase information of the self-statistical points, respectively. The coefficient matrix of the polynomial. The corresponding weight matrix has its element values ​​calculated using formula (7). The coefficients are based on the least squares optimization estimation. In addition to the location information of the target monitoring points, the atmospheric phase at the corresponding target monitoring point can be calculated. Subtracting the atmospheric phase estimated using the least squares method from the optimized phase can reduce the adverse effects of the atmospheric phase and improve the reliability of subsequent deformation monitoring. For any target monitoring point, the location information of the target monitoring point is multiplied by a coefficient. The atmospheric phase of the target monitoring point is obtained. The corresponding atmospheric phase is subtracted from the phase of the target monitoring point optimized in S102 to obtain the phase of the target monitoring point after compensation.

[0050] like Figure 3 As shown, Figure 3Schematic diagram of spatial scale difference between atmospheric delay phase and related phase of geological disasters.

[0051] S104, determine a suspected geological disaster area in a range where the target diversion channel is located according to the compensated phase of the target monitoring point.

[0052] In one embodiment, determining a suspected geological disaster area in a range where the target diversion channel is located according to the compensated phase of the target monitoring point includes: calculating a deformation value at each target monitoring point according to the compensated phase of the target monitoring point; calculating a local Moran's coefficient of each target monitoring point according to the deformation value at each target monitoring point; calculating a deformation threshold based on the target monitoring point whose local Moran's coefficient is greater than 0, and determining the target monitoring point whose deformation value is greater than the deformation threshold as the suspected geological disaster area.

[0053] According to the compensated phase of the target monitoring point, the deformation rate of each target monitoring point is calculated based on the time series synthetic aperture radar image by means of the discrete point unwrapping method, and then the cumulative deformation variable, i.e. the deformation value, of each target monitoring point is determined according to the time span of the time series synthetic aperture radar image.

[0054] The deformation rate and the time series cumulative deformation distribution of the region of interest are calculated. The deformation rate is, for example, 10 mm / y, i.e. the annual average rate is 10 mm; the time span of the time series SAR image is 2 years, and the time series cumulative deformation is, for example, 20 mm, which means that the cumulative deformation in 2 years is 20 mm.

[0055] The calculation formula of the local Moran's coefficient is: (9); wherein, is the local Moran's coefficient at position , denotes the deformation value at position , is the mean value of all position deformation values, is the spatial weight matrix element, indicating the spatial relationship between position and position , adjacent is 1 and non-adjacent is 0, is the number of target monitoring points. The corresponding deformation threshold and the local Moran's coefficient for calculating the deformation distribution are set to distinguish the heterogeneity of the deformation area, and then the suspected geological disaster area is screened out according to the deformation information.

[0056] According to the deformation value at each target monitoring point, the local Moran's coefficient of each target monitoring point is calculated based on formula (9).

[0057] According to the deformation value at each target monitoring point, the local Moran's coefficient of each target monitoring point is calculated based on formula (9). ​

[0058] The deformation values of the target monitoring points with the local Morlet coefficient greater than 0 are obtained, and the average value of the deformation values is determined as a deformation threshold value.

[0059] In S105, the temperature data of the suspected geological disaster area is determined according to the high-resolution optical and thermal infrared images of the suspected geological disaster area, and the disaster anomaly coefficient of the suspected geological disaster area is determined according to the temperature data and the deformation value of the suspected geological disaster area.

[0060] In the suspected geological disaster area, the high-resolution optical and thermal infrared images of the suspected geological disaster area are obtained by using the optical and thermal infrared cameras carried by the unmanned aerial vehicle, and the temperature data of the suspected geological disaster area is obtained; the geological disaster anomaly coefficient is calculated based on the temperature and deformation, as shown in formula (10); according to the disaster anomaly coefficient, the geological disaster situation of the corresponding area is indicated, and the greater the value is, the higher the risk is. The calculation formula of the disaster anomaly coefficient of the suspected geological disaster area is:

[0061] (10); Wherein, represents the deformation value of the i-th pixel in the deformation statistical window of the suspected geological disaster area; represents the area corresponding to the i-th pixel in the deformation statistical window; represents the number of pixels in the deformation statistical window; represents the area corresponding to the i-th pixel in the temperature statistical window of the suspected geological disaster area; represents the temperature of the i-th pixel in the temperature statistical window; represents the average temperature of the pixels in the temperature statistical window; represents the number of pixels in the temperature statistical window. It should be noted that the calculation of the disaster anomaly coefficient is not for a single point, but for the statistical information of all points in a certain window, so the statistical window is the window size when calculating the disaster anomaly coefficient, which can be 11x11, 15x15, etc., and can be set as needed. As shown in FIG. 1,

[0062] It should be noted that the calculation of the disaster anomaly coefficient is not for a single point, but for the statistical information of all points in a certain window, so the statistical window is the window size when calculating the disaster anomaly coefficient, which can be 11x11, 15x15, etc., and can be set as needed.

[0063] As shown in FIG. 1, Figure 4 Figure 4 is a schematic diagram for calculating and applying the geological disaster index along the water diversion channel.

[0064] In S106, the geological disaster area of the target water diversion channel is determined according to the disaster anomaly coefficient of the suspected geological disaster area.

[0065] ​​​In one embodiment, a suspected geological disaster area with a catastrophe anomaly coefficient greater than a safety coefficient threshold is determined as a geological disaster area of the target diversion channel.

[0066] And, after determining the geological disaster area of the target diversion channel, a geological disaster distribution map is drawn for on-site disaster verification.

[0067] In one embodiment, as shown in Figure 5 the present application also provides a flowchart of a geological disaster monitoring method for a diversion channel of a non-pressure diversion type hydropower station, which includes the following steps: S501, acquiring time-series SAR data according to the area where the diversion channel is located and the monitoring period.

[0068] S502, based on the acquired time-series SAR data, selecting a super master image using the principle of minimum space-time baseline and minimum spectral center difference, and realizing SAR data registration.

[0069] S503, self-adaptive cropping of the region of interest according to the geometric distribution of the diversion channel and the terrain complexity along the DEM.

[0070] S504, generation of an optimal space-time baseline map.

[0071] Among them, considering the local climate characteristics, avoiding the influence of snow season images, generating a space-time baseline map, ensuring the maximum availability and coherence of SAR images, avoiding the dependence on the connection of a single image pair, and guaranteeing and improving the quality of subsequent inversion results.

[0072] S505, based on the registered SAR images, using HTCI and EMI methods to carry out homogeneous pixel extraction and phase optimization, and then determining DS candidate points based on the phase consistency principle.

[0073] S506, based on the intensity dispersion index, selecting PS candidate points, and fusing with DS candidate points to form a monitoring target point candidate dataset, and then based on the phase consistency principle, carrying out phase stability estimation, and selecting points with high stability as the final target monitoring points.

[0074] S507, according to the distribution of the monitoring point positions in the cropped region, selecting part of the point positions by random sampling, and calculating the statistical stability of the phase change within a certain window range, and then based on it, carrying out weighted least squares parameter optimization estimation, obtaining a polynomial fitting the atmospheric related phase, and subtracting from the corresponding phase, thereby realizing atmospheric estimation and compensation.

[0075] S508, based on the atmospheric corrected differential interferometric phase and monitoring point position information, constructing a delaunnay triangular network, phase unwrapping and time series analysis, and finally obtaining the deformation rate and time series cumulative deformation of the corresponding region.

[0076] S509, on the basis of setting a deformation threshold, and in combination with a local model index based on deformation, jointly determining a deformation suspected disaster area.

[0077] S510, carrying out unmanned aerial vehicle thermal infrared data acquisition in the suspected disaster area.

[0078] S511, using deformation information and thermal infrared data, calculating a geological disaster anomaly coefficient of the corresponding area according to deformation rate and time-series cumulative deformation.

[0079] S512, based on the obtained geological disaster anomaly coefficient, carrying out geological disaster distribution mapping along the diversion canal.

[0080] S513, according to the geological disaster distribution map and the corresponding disaster coefficient size, carrying out artificial on-site verification.

[0081] The present application fully integrates the advantages of satellite-borne SAR and unmanned aerial vehicle thermal infrared data in the monitoring of geological disasters along the diversion canal, obtains high-precision deformation information of the ground surface along the diversion canal through the time-series InSAR technology based on satellite-borne SAR, overcomes the omission problem caused by traditional fixed-width screening and cutting and the problem of accurate atmospheric estimation, and on this basis, in combination with high-resolution thermal infrared remote sensing data of the unmanned aerial vehicle, constructs a geological disaster index along the diversion canal, realizes quantitative analysis and mapping of the geological disasters of the diversion canal, and provides a new method for identification and monitoring of geological disasters along the diversion canal.

[0082] The present application has the following advantages: 1. The present application overcomes the shortcomings of traditional fixed-width cutting along the diversion canal, establishes a region cutting method that is adaptive to terrain changes, and can effectively avoid and reduce the omission of remote geological disasters; 2. The present application carries out atmospheric phase estimation and compensation based on the phase stability weighting of the diversion canal itself, can specifically weaken the influence of the atmosphere on the accurate extraction of deformation along the linearly distributed diversion canal, and realizes high-precision extraction of geological disaster-related ground deformation; 3. By means of the seepage phenomenon that may be accompanied by geological disasters of the diversion canal, on the basis of deformation extraction, the present application uses the ground temperature information obtained from high-resolution thermal infrared images of the unmanned aerial vehicle, combines ground deformation and temperature information, constructs an index for indicating geological disasters along the diversion canal, makes up for the shortcomings of single deformation information, can more intuitively and quantitatively describe the possibility of geological disasters along the diversion canal, and on this basis, carries out geological disaster mapping of the diversion canal, which can more intuitively and accurately serve the detection and identification of geological disasters along the diversion canal.

[0083] In the application, the method for monitoring geological disasters of the water diversion channel of the non-pressure diversion type hydropower station can be executed without considering the order of the steps shown in the method, and the order of the steps can be determined according to actual needs, and the application does not limit the order of the steps. Figure 1 The order of the execution of the steps is shown in the method, and the order of the execution of the steps can be determined according to actual needs, and the application does not limit the order of the steps.

[0084] The method for monitoring geological disasters of the water diversion channel of the non-pressure diversion type hydropower station provided by one or more embodiments of the application is based on the same idea, and the application further provides a corresponding device for monitoring geological disasters of the water diversion channel of the non-pressure diversion type hydropower station. The acquisition module is configured to acquire time-series synthetic aperture radar images in a range where the target water diversion channel is located, and register all the synthetic aperture radar images to the same space-time reference. The cropping module is configured to determine a region of interest according to changes in the terrain slope ratio along the target water diversion channel, determine homogeneous pixels from the registered synthetic aperture radar images in the direction of extending the window along the strike of the water diversion channel, and perform phase optimization, determine candidate distributed target pixels from the phase-optimized homogeneous pixels according to the phase consistency discrimination method. The compensation module is configured to acquire candidate points of permanent scatterers in the time-series synthetic aperture radar images, fuse the candidate points of permanent scatterers with the candidate distributed target pixels to obtain target monitoring points, select self-statistical points according to the distribution of the target monitoring points, calculate the phase stability of each target monitoring point according to the self-statistical points, and compensate the atmospheric phase of the target monitoring points by using the weighted least squares method based on the phase stability. The calculation module is configured to calculate the deformation value of each position in the range where the target water diversion channel is located according to the phase of the target monitoring points after compensation, and determine the suspected geological disaster area in the range where the target water diversion channel is located according to the deformation value of each position. The first determination module is configured to determine the temperature data of the suspected geological disaster area according to the high-resolution optical and thermal infrared images of the suspected geological disaster area, and determine the disaster anomaly coefficient of the suspected geological disaster area according to the temperature data and the deformation value of the suspected geological disaster area. The second determination module is configured to determine the geological disaster area of the target water diversion channel according to the disaster anomaly coefficient of the suspected geological disaster area.

[0085] The specific limitation of the geological disaster monitoring device for the diversion channel of the non-pressure diversion type hydropower station can refer to the limitation of the geological disaster monitoring method for the diversion channel of the non-pressure diversion type hydropower station, which will not be repeated here. Each module in the above-mentioned geological disaster monitoring device for the diversion channel of the non-pressure diversion type hydropower station can be realized by software, hardware and their combination. The above-mentioned modules can be embedded in or independent of the processor in the computer device in hardware form, or can be stored in the memory in the computer device in software form, so as to call and execute the operation corresponding to each module by the processor.

[0086] The application further provides a computer readable storage medium, which stores a computer program, and the computer program can be used to execute the above-mentioned Figure 1 The application provides a geological disaster monitoring method for a diversion channel of a non-pressure diversion type hydropower station.

[0087] The application further provides a computer device, which comprises a processor, an internal bus, a network interface, a memory and a non-volatile memory at the hardware level, and can further comprise other hardware required by business. The processor reads the corresponding computer program from the non-volatile memory into the memory and then runs, so as to realize the above-mentioned Figure 1 The application provides a geological disaster monitoring method for a diversion channel of a non-pressure diversion type hydropower station.

[0088] Those skilled in the art can understand that all or part of the above-mentioned embodiment methods can be completed by a computer program to instruct related hardware. The computer program can be stored in a non-volatile computer readable storage medium, and when the computer program is executed, the computer program can include the flow of the above-mentioned embodiment method. In each embodiment of the application, any reference to the memory, storage, database or other medium can include at least one of the non-volatile and volatile memory. The non-volatile memory can include a read-only memory (ROM), a magnetic tape, a floppy disk, a flash memory or an optical memory. The volatile memory can include a random access memory (RAM) or an external cache memory. As an illustration but not limitation, the RAM can be in various forms, such as static random access memory (SRAM) or dynamic random access memory (DRAM).

[0089] The technical features of the above embodiments can be combined in any manner. For the sake of brevity, not all possible combinations of the technical features in the above embodiments are described, however, as long as the combinations of the technical features do not contradict each other, they should be considered to be within the scope of the present application.

Claims

1. A method for monitoring geological disasters in a water diversion channel of a non-pressure diversion type hydropower station, characterized in that, include: Acquire temporal synthetic aperture radar images of the area where the target water diversion canal is located, and register all synthetic aperture radar images to the same spatiotemporal reference. Based on the changes in the terrain slope ratio along the target irrigation canal, the region of interest is determined; within the region of interest, the window is expanded along the direction of the irrigation canal, and homogeneous pixels are identified from the registered synthetic aperture radar image and phase optimization is performed; candidate distributed target pixels are identified from the phase-optimized homogeneous pixels based on the phase consistency discrimination method. Acquire candidate points of permanent scatterers in temporal synthetic aperture radar images, and fuse the candidate points of permanent scatterers with candidate distributed target pixels to obtain the target monitoring points; Based on the distribution of target monitoring points, self-statistical points are selected, and the phase stability of each target monitoring point is calculated based on the self-statistical points. Then, based on the phase stability, the atmospheric phase of the target monitoring points is compensated using the weighted least squares method. Based on the compensated phase of the target monitoring points, the suspected geological disaster area within the range of the target water diversion canal is determined; Based on high-resolution optical and thermal infrared images of suspected geological disaster areas, temperature data of suspected geological disaster areas are determined, and based on temperature data and deformation values ​​of suspected geological disaster areas, disaster anomaly coefficients of suspected geological disaster areas are determined. Based on the catastrophic anomaly coefficient of the suspected geological disaster area, the geological disaster area of ​​the target water diversion canal is determined.

2. The method of claim 1, wherein, Based on the changes in the topographic slope ratio along the target irrigation canal, the region of interest is determined, including: The width of the target water diversion canal along the route is determined based on the slope ratio of the terrain along the route. The region of interest is determined based on the width of the target water diversion channel.

3. The method of claim 2, wherein, The formula for calculating the width of the target irrigation canal is as follows: ； wherein, is the line width, is the terrain slope ratio, is the base width.

4. The method of claim 1, wherein, Phase optimization for homogeneous pixels includes: Construct the sample complex coherence matrix based on all homogeneous pixels; Based on the sample complex coherence matrix, the maximum likelihood estimation phase is solved using the eigenvalue decomposition method, and the temporal phase sequence of each homogeneous pixel is optimized.

5. The method of claim 4, wherein, The method for constructing the sample complex coherence matrix is ​​as follows: ； wherein, is a sample complex coherence matrix, is is a vector of amplitude normalized complex observations of a synthetic aperture radar image at a distributed pixel, is is a conjugate transpose of is a set of homogeneous pixels of the distributed target, is a number of corresponding homogeneous pixels.

6. The method of claim 5, wherein, The eigenvalue decomposition method is as follows: ； wherein, is a scale factor, is a sequence of timing phases optimized by the EMI method, the solution of the eigen decomposition method is is a sequence of timing phases optimized by the EMI method, the solution of the eigen decomposition method is is a sequence of timing phases optimized by the EMI method, the solution of the eigen decomposition method is 7. The method of claim 1, wherein, Based on the compensated phase of the target monitoring points, the suspected geological disaster area within the scope of the target water diversion canal is determined, including: Based on the compensated phase of the target monitoring point, calculate the deformation value of each target monitoring point; Calculate the local Moran coefficient at each target monitoring point based on the deformation value at each target monitoring point; Deformation thresholds are calculated for target monitoring points with local Moran coefficients greater than 0, and target monitoring points with deformation values ​​greater than this deformation threshold are identified as suspected geological disaster areas.

8. The method of claim 1, wherein, The formula for calculating the catastrophic anomaly coefficient in a suspected geological disaster area is as follows: ； wherein, represents the deformation value at the i-th pixel in the deformation statistical window of the suspected geologic disaster area; represents the area corresponding to the i-th pixel in the deformation statistical window; represents the number of pixels in the deformation statistical window; represents the area corresponding to the i-th pixel in the temperature statistical window of the suspected geologic disaster area; represents the temperature of the i-th pixel in the temperature statistical window; represents the average temperature of the pixels in the temperature statistical window; represents the number of pixels in the temperature statistical window.​​​​ 9. The method of claim 1, wherein, Based on the catastrophic anomaly coefficient of the suspected geological disaster area, the geological disaster area of ​​the target water diversion canal is determined, including: Areas suspected of geological disasters with a disaster anomaly coefficient greater than the safety coefficient threshold are identified as geological disaster areas for the target water diversion canal.

10. The method of claim 1, wherein, The method further includes: After identifying the geological disaster area of ​​the target water diversion channel, a geological disaster distribution map is drawn for on-site disaster verification.

Images