Disaster hazard InSAR rapid identification and accurate positioning method and system

Abstract

The application relates to the technical field of geological disaster monitoring, in particular to a disaster hidden danger InSAR rapid identification and accurate positioning method and system, the method comprising the following steps: obtaining original polarized interference data of a target region through a dual-band SAR receiving antenna and performing enhancement processing to obtain enhanced polarized interference data; performing terrain self-adaptive radiation correction and image registration based on terrain constraints on the enhanced polarized interference data to obtain pretreated interference image data; obtaining initial deformation field data through phase unwrapping; determining spatial strain gradient and time sequence evolution acceleration based on the initial deformation field data, and constructing a three-dimensional deformation tensor feature space according to the spatial strain gradient and the time sequence evolution acceleration; determining a disaster probability prediction result based on the three-dimensional deformation tensor feature space; and constructing a visual disaster probability prediction product view and displaying the same according to the disaster probability prediction result. The application has the effect of improving the identification timeliness and positioning accuracy of geological disaster hidden dangers.

Description

Technical Field

[0001] This application relates to the field of geological disaster monitoring technology, and in particular to an InSAR method and system for rapid identification and precise location of potential disaster hazards. Background Technology

[0002] Currently, Synthetic Aperture Radar Interferometry (InSAR) is a remote sensing technique that utilizes the phase information of synthetic aperture radar to monitor surface deformation. This technology can provide high-precision displacement measurement data over a wide area and in geological disaster early warning, ground subsidence analysis, and deformation monitoring of major engineering projects, all around the clock.

[0003] Existing disaster monitoring solutions mostly employ single-band radar data combined with conventional interferometry procedures, identifying potential hazards by calculating the numerical values ​​of single deformation displacements at surface points. This approach lacks efficient parallel computing support in the preprocessing stage and often focuses solely on displacement magnitude in risk assessment, neglecting the spatially uneven gradient of deformation and its accelerated temporal evolution. This results in the timeliness and accuracy of identifying small-scale disaster hazards in complex terrain failing to meet practical early warning requirements. Therefore, a comprehensive solution capable of rapid InSAR data processing, intelligent hazard identification, and precise location is urgently needed. Summary of the Invention

[0004] To improve the timeliness and accuracy of geological hazard identification and location, this application provides a method and system for rapid InSAR identification and precise location of geological hazards.

[0005] The above-mentioned objective of this application is achieved through the following technical solution:

[0006] A method for rapid InSAR identification and precise location of potential disaster hazards, the method comprising:

[0007] The raw polarization interferometric data of the target area is obtained by using a dual-band SAR receiving antenna;

[0008] The original polarimetric interferometric data is enhanced to obtain enhanced polarimetric interferometric data;

[0009] The enhanced polarimetric interferometric data is subjected to terrain-adaptive radiometric correction, and the enhanced polarimetric interferometric data after terrain-constrained quadratic surface fitting algorithm is used for image registration to obtain preprocessed interferometric image data.

[0010] The preprocessed interferometric image data is subjected to phase unwrapping to obtain initial deformation field data;

[0011] Based on the initial deformation field data, the spatial strain gradient and temporal evolution acceleration are determined, and a three-dimensional deformation tensor feature space is constructed according to the spatial strain gradient and the temporal evolution acceleration.

[0012] Based on the aforementioned three-dimensional deformation tensor feature space, the disaster probability prediction results are determined;

[0013] Based on the disaster probability prediction results, a visual disaster probability prediction product view is constructed and presented through the display interface of the disaster probability prediction platform.

[0014] By adopting the above technical solutions, and through dual-band SAR acquisition and enhancement processing, it is possible to achieve multi-source information complementarity and signal gain enhancement in disaster areas at different frequency bands. This solves the problem of insufficient perception capability of a single band in complex environments and provides high-quality raw data for subsequent calculations. Through terrain-adaptive radiometric correction and terrain-constrained quadratic surface fitting image registration, it is possible to eliminate the uneven radiation distribution caused by complex terrain and achieve precise alignment of main and auxiliary images. This effectively corrects image geometric distortion and significantly improves the coherence of multi-phase images. By constructing a three-dimensional deformation tensor feature space containing spatial strain gradient and temporal evolution acceleration, it is possible to transform displacement values ​​into multi-dimensional feature vectors reflecting the dynamic evolution of the Earth's surface. This allows for the accurate capture of critical signals for the accelerated instability and transformation of hidden dangers. Through visualization, it is possible to achieve automated probability prediction and three-dimensional intuitive positioning of disaster hazards. This overcomes the shortcomings of existing technologies, such as long data processing cycles, low accuracy in identifying small-scale hazards, and high false alarm rates, and enables rapid identification and accurate positioning of disaster hazards.

[0015] In a preferred embodiment, this application can be further configured such that: the enhancement processing of the original polarimetric interferometric data to obtain enhanced polarimetric interferometric data specifically includes:

[0016] Calculate the signal arrival time difference between each array element in the dual-band SAR receiving antenna, and perform phase weighting compensation on the original polarization interferometric data based on the signal arrival time difference;

[0017] The original polarimetric interferometric data after phase weighting compensation is denoised to obtain high signal-to-noise ratio polarimetric interferometric data.

[0018] The high signal-to-noise ratio polarized interferometric data is compressed in the range direction using the ω-k algorithm to obtain the range-compressed data.

[0019] The range-compressed data is subjected to azimuth self-focusing processing based on the phase gradient estimation method to obtain the enhanced polarization interferometry data.

[0020] By adopting the above technical solutions, and by performing phase weighted compensation and using the ω-k algorithm and phase gradient estimation method for range compression and azimuth autofocus, it is possible to achieve in-depth optimization of the original polarimetric interferometric data in the time domain, frequency domain, and spatial domain, thereby significantly improving the signal-to-noise ratio and geometric focusing accuracy of the original image.

[0021] In a preferred embodiment, this application can be further configured as follows: The enhanced polarimetric interferometric data undergoes terrain-adaptive radiometric correction, and a terrain-constrained quadratic surface fitting algorithm is used to perform image registration on the terrain-adaptive radiometric corrected enhanced polarimetric interferometric data to obtain preprocessed interferometric image data. Specifically, this includes:

[0022] Obtain the external DEM model data of the target area, and calculate the initial geometric offset of the auxiliary image relative to the main image in the enhanced polarimetric interferometry data based on the external DEM model data;

[0023] Using the DEM model of the target area, the local slope and azimuth of each sampling point in the DEM model are calculated;

[0024] Based on the relationship between the local slope, azimuth angle and satellite incident angle, the radiation normalization coefficient is calculated. The enhanced polarization interferometric data is then processed using the radiation normalization coefficient to eliminate the difference in terrain brightness, thereby obtaining terrain adaptive radiometric correction data.

[0025] The search window is determined using the initial geometric offset in the cross-correlation peak region of the main image and the auxiliary image, and a quadratic surface equation about the cross-correlation coefficient is constructed by least squares fitting.

[0026] Based on the quadratic surface equation, an offset vector is determined, and the offset vector is used for image resampling and alignment to align the auxiliary image with the main image, thereby obtaining the preprocessed interferometric image data.

[0027] By adopting the above technical solution, using external DEM model data for terrain adaptive radiometric correction and registration based on cross-correlation coefficient quadratic surface fitting, it is possible to offset the false brightness differences caused by the relationship between the hillside orientation and the satellite incident angle and lock the offset vector with extremely high precision. This ensures that the main and auxiliary images achieve complete physical overlap in terms of ground feature details and significantly reduces the interference phase noise caused by registration deviation.

[0028] In a preferred embodiment, this application can be further configured such that: performing phase unwrapping processing on the preprocessed interferometric image data to obtain initial deformation field data specifically includes:

[0029] Image pairs whose temporal and vertical baselines are both within the corresponding preset thresholds are selected, and an interferogram stack consisting of at least a preset number of the preprocessed interferometric image data is constructed.

[0030] The phase map in the interferogram stack is divided into multiple sub-blocks, and the sub-blocks are processed in parallel based on the statistical cost flow unwrapping algorithm to obtain continuous phase stack data.

[0031] Based on the continuous phase stack data, a deformation observation equation is constructed. The structure matrix of the deformation observation equation is solved by generalized inverse solution using singular value decomposition to obtain the initial value of the deformation rate of each monitoring point.

[0032] The atmospheric delay phase in the initial value of the deformation rate is identified using a preset turbulent atmospheric model, and the atmospheric delay phase is removed from the initial value of the deformation rate to obtain the initial deformation field data.

[0033] By adopting the above technical solutions, accelerating the unwrapping of statistical cost flow and using the singular value decomposition method to solve the deformation equation in a generalized inverse, it is possible to solve the matrix rank deficiency problem caused by data discontinuity while ensuring the timeliness of massive data processing. By introducing a turbulent atmospheric model to remove the water vapor delay phase, it is possible to ensure that the calculated initial deformation field data excludes atmospheric interference, thereby obtaining the true surface line-of-sight displacement values.

[0034] In a preferred embodiment, this application can be further configured as follows: determining the spatial strain gradient and temporal evolution acceleration based on the initial deformation field data, and constructing a three-dimensional deformation tensor feature space according to the spatial strain gradient and the temporal evolution acceleration, specifically includes:

[0035] Extract the line-of-sight cumulative deformation and deformation time series of each monitoring point from the initial deformation field data;

[0036] Based on the cumulative deformation along the line of sight, the first spatial partial derivative between the target monitoring point and the neighboring monitoring points is calculated to obtain the spatial strain gradient.

[0037] Perform a second-order difference operation on the time series to calculate the time series evolution acceleration;

[0038] The three-dimensional deformation tensor feature space is obtained by vector mapping the cumulative deformation along the line of sight, the spatial strain gradient, and the temporal evolution acceleration.

[0039] By adopting the above technical solution, by calculating the first-order spatial partial derivative of the surface displacement to obtain the spatial strain gradient and performing second-order difference operations to obtain the temporal evolution acceleration, a deep correlation between the spatial non-uniform settlement characteristics and the temporal nonlinear evolution characteristics can be established, thereby forming a digital tensor feature that comprehensively depicts the dynamic evolution of disaster hazards, greatly enhancing the system's ability to identify the evolution stage of hazards.

[0040] In a preferred embodiment, this application can be further configured such that: determining the disaster probability prediction result based on the three-dimensional deformation tensor feature space specifically includes:

[0041] The three-dimensional deformation tensor feature space is input into a pre-trained ResNet-50 convolutional neural network, and the spatial morphological features of the disaster-prone area are obtained through residual convolutional layers.

[0042] Based on the spatial morphological features, an attention mechanism is used to assign weights to the temporal evolution sequence in the three-dimensional deformation tensor to obtain a temporal fusion feature vector.

[0043] The time-series fusion feature vector is input into a random forest classifier for nonlinear classification processing to output the disaster probability prediction result for each monitoring point.

[0044] By adopting the above technical solution, extracting spatial morphological features through the ResNet-50 convolutional neural network and introducing an attention mechanism to assign weights to key temporal features, the model can automatically focus on the abnormal moments in the deformation sequence that contribute the most to disaster instability, thereby effectively identifying and filtering out false deformation signals caused by seasonal environmental factors and improving the confidence of disaster probability prediction results.

[0045] In a preferred embodiment, this application can be further configured as follows: The step of constructing a visualized disaster probability prediction product view based on the disaster probability prediction results, and presenting the disaster probability prediction product view through the display interface of the disaster probability prediction platform, specifically includes:

[0046] The WebGL engine is used to load the three-dimensional point cloud data of the target area, and a heat map layer is rendered on the three-dimensional point cloud data according to the disaster probability prediction results.

[0047] Based on the geographical range corresponding to the heat map, GIS analysis tools are used to calculate the slope parameters of disaster hazard points and the range of threatened watersheds.

[0048] Based on the three-dimensional geographic coordinates of the disaster hazard points, the slope parameters, and the threatened watershed range, an early warning report is generated;

[0049] Based on preset layer display rules, the early warning report is processed by layer synthesis to generate the disaster probability prediction product view.

[0050] By adopting the above technical solutions, and by using the WebGL rendering engine to perform volume rendering and fusion, and by using GIS analysis tools to calculate slope parameters and the scope of threatened watersheds, an intuitive mapping relationship between risk prediction probability and digital geographic scene can be established, and structured early warning reports can be automatically generated. This provides disaster prevention decision-making departments with a fast, accurate and highly interactive disaster location product view.

[0051] The second objective of this invention is achieved through the following technical solution:

[0052] A method for rapid InSAR identification and precise location of potential disaster hazards, the method comprising:

[0053] The data acquisition module is used to acquire the raw polarimetric interferometric data of the target area through a dual-band SAR receiving antenna;

[0054] The enhancement processing module is used to enhance the original polarization interference data to obtain enhanced polarization interference data;

[0055] The data preprocessing module is used to perform terrain-adaptive radiometric correction on the enhanced polarimetric interferometric data, and to perform image registration on the enhanced polarimetric interferometric data after terrain-adaptive radiometric correction using a terrain-constrained quadratic surface fitting algorithm to obtain preprocessed interferometric image data.

[0056] The deformation field extraction module is used to perform phase unwrapping processing on the preprocessed interferometric image data to obtain initial deformation field data.

[0057] The feature space construction module is used to determine the spatial strain gradient and temporal evolution acceleration based on the initial deformation field data, and to construct a three-dimensional deformation tensor feature space according to the spatial strain gradient and the temporal evolution acceleration.

[0058] The probability prediction module is used to determine the disaster probability prediction result based on the three-dimensional deformation tensor feature space.

[0059] The visualization module is used to construct a visualized disaster probability prediction product view based on the disaster probability prediction results, and to present the disaster probability prediction product view through the display interface of the disaster probability prediction platform.

[0060] By adopting the above technical solutions, and through dual-band SAR acquisition and enhancement processing, it is possible to achieve multi-source information complementarity and signal gain enhancement in disaster areas at different frequency bands. This solves the problem of insufficient perception capability of a single band in complex environments and provides high-quality raw data for subsequent calculations. Through terrain-adaptive radiometric correction and terrain-constrained quadratic surface fitting image registration, it is possible to eliminate uneven radiation distribution caused by complex terrain and achieve precise alignment of primary and auxiliary images. This effectively corrects image geometric distortion and significantly improves the coherence of multi-phase images. By constructing a three-dimensional deformation tensor feature space containing spatial strain gradient and temporal evolution acceleration, displacement values ​​can be transformed into multi-dimensional feature vectors reflecting the dynamic evolution of the Earth's surface. This accurately captures the critical signals of accelerated instability and transformation of hidden dangers. Through visualization, it is possible to achieve automated probabilistic prediction and intuitive three-dimensional positioning of disaster hazards. This overcomes the shortcomings of existing technologies, such as long data processing cycles, low accuracy in identifying small-scale hazards, and high false alarm rates, enabling rapid identification and accurate positioning of disaster hazards.

[0061] In summary, this application includes at least one of the following beneficial technical effects:

[0062] 1. By acquiring and enhancing dual-band SAR, it is possible to achieve multi-source information complementarity and signal gain enhancement in disaster areas at different frequency bands, thereby solving the problem of insufficient perception capability of a single band in complex environments and providing high-quality raw data for subsequent calculations. Through terrain-adaptive radiometric correction and terrain-constrained quadratic surface fitting image registration, it is possible to eliminate the uneven radiation distribution caused by complex terrain and achieve precise alignment of main and auxiliary images, thereby effectively correcting image geometric distortion and significantly improving the coherence of multi-phase images. By constructing a three-dimensional deformation tensor feature space containing spatial strain gradient and temporal evolution acceleration, it is possible to transform displacement values ​​into multi-dimensional feature vectors reflecting the dynamic evolution of the Earth's surface, thereby accurately capturing the critical signal of accelerated instability transformation of hidden dangers. Through visualization, it is possible to achieve automated probability prediction and three-dimensional intuitive positioning of disaster hazards, thereby overcoming the shortcomings of existing technologies such as long data processing cycles, low accuracy of small-scale hazard identification, and high false alarm rate, and achieving rapid identification and accurate positioning of disaster hazards.

[0063] 2. By using the WebGL rendering engine to perform volume rendering and fusion, and by using GIS analysis tools to calculate slope parameters and the scope of threatened watersheds, an intuitive mapping relationship between risk prediction probability and digital geographic scene can be established, and structured early warning reports can be automatically generated. This provides disaster prevention decision-making departments with a fast, accurate and highly interactive disaster location product view. Attached Figure Description

[0064] Figure 1This is a flowchart illustrating the implementation of an InSAR method for rapid identification and precise location of disaster hazards in one embodiment of this application.

[0065] Figure 2 This is a principle block diagram of an InSAR rapid identification and precise positioning system for disaster hazards according to one embodiment of this application. Detailed Implementation

[0066] The following embodiments will help those skilled in the art to further understand the function of this application, but do not limit this application in any way. It should be noted that those skilled in the art can make several modifications and improvements without departing from the concept of this application. These all fall within the protection scope of this application.

[0067] In the following description, specific details such as particular system architectures and techniques are set forth for illustrative purposes and not for limitation, in order to provide a thorough understanding of the embodiments of this application. However, those skilled in the art will understand that this application may also be implemented in other embodiments without these specific details. In other instances, detailed descriptions of well-known systems, apparatuses, circuits, and methods have been omitted so as not to obscure the description of this application with unnecessary detail.

[0068] It should be understood that, when used in this application specification and the appended claims, the term "comprising" indicates the presence of the described features, integrals, steps, operations, elements and / or components, but does not exclude the presence or addition of one or more other features, integrals, steps, operations, elements, components and / or a collection thereof.

[0069] The present application will be further described in detail below with reference to the accompanying drawings.

[0070] In one embodiment, such as Figure 1 As shown, this application discloses a method for rapid InSAR identification and accurate location of disaster hazards, which specifically includes the following steps:

[0071] S10. Obtain the original polarization interferometric data of the target area through the dual-band SAR receiving antenna.

[0072] Specifically, a dual-band SAR receiving antenna refers to an induction device capable of synchronously responding to electromagnetic echoes in different frequency bands. For example, it can simultaneously utilize the X-band, which has high spatial resolution, and the L-band, which has strong vegetation penetration capability, for coordinated observation. Raw polarization interferometric data refers to the original electrical signal sequence containing complex amplitude and initial phase information formed after the radar transmits and receives specific polarization signals to target areas such as mountains, cities, or mining areas. Where A is the complex amplitude and Φ is the initial phase (in rad). By acquiring this multi-source, multi-band data set, we can provide an initial data foundation with rich ground object scattering characteristics and phase depth information for subsequent identification of small-scale disaster hazards, thereby ensuring that the system can cover the real surface deformation state under different materials and coverage at the perception level.

[0073] S20. Enhance the original polarization interferometric data to obtain enhanced polarization interferometric data.

[0074] Specifically, enhancement processing refers to a series of primary processing actions performed on the raw signals acquired in the preceding steps to improve data quality and geometric fidelity. By mathematically compensating and correcting phase deviations in the raw polarimetric interferometric data caused by propagation path interference, satellite attitude fluctuations, or environmental thermal noise, the chaotic raw echo stream can be transformed into a set of digital images with clear focus, accurate geometric topology, and significantly improved signal-to-noise ratio. This enhanced polarimetric interferometric data, as the sole input source for subsequent precise deformation calculations, directly determines the clarity of the final disaster identification and the sensitivity to capturing minute displacement features.

[0075] S30. Perform terrain-adaptive radiometric correction on the enhanced polarimetric interferometric data, and use a terrain-constrained quadratic surface fitting algorithm to perform image registration on the enhanced polarimetric interferometric data after terrain-adaptive radiometric correction to obtain preprocessed interferometric image data.

[0076] Specifically, terrain-adaptive radiometric correction refers to dynamically adjusting the echo energy distribution of each sampling point in the image based on the topographic relief of the target area to eliminate radiometric distortion caused by mountain shadows or overlapping phenomena. The image alignment criterion uses a terrain-constrained mathematical fitting model to find the geometric correlation between the main image and the auxiliary image at the sub-pixel level. Through this precise alignment operation, enhanced polarimetric interferometric data acquired at different times can achieve complete overlap in geospatial coordinates, thereby generating preprocessed interferometric image data with extremely high coherence and eliminating false features of terrain brightness and darkness, providing a robust pixel alignment benchmark for subsequent extraction of millimeter-level deformation phase.

[0077] S40. Perform phase unwrapping on the preprocessed interferometric image data to obtain the initial deformation field data.

[0078] Specifically, phase unwrapping refers to the process of converting the wrapped phase values ​​extracted from the preceding aligned image within a single period into continuous phase values ​​that reflect the actual physical displacement through mathematical integration or path search algorithms. By transforming the abstract phase difference information into a deformation scalar with clear physical meaning, it is possible to construct initial deformation field data covering the entire target area. This data intuitively represents the cumulative displacement value of each monitoring point on the ground relative to the satellite line of sight during the observation period, and is the core numerical basis for subsequent disaster evolution trend analysis.

[0079] S50. Based on the initial deformation field data, determine the spatial strain gradient and temporal evolution acceleration, and construct a three-dimensional deformation tensor feature space according to the spatial strain gradient and temporal evolution acceleration.

[0080] Specifically, spatial strain gradient refers to the physical index used to quantitatively describe whether there is local uneven settlement or shear failure on the ground surface by calculating the difference in deformation displacement between adjacent spatial units in the initial deformation field. Temporal evolution acceleration is a dynamic index extracted by analyzing the rate of change of deformation over time, reflecting whether the ground surface movement is in an accelerated instability state. By encapsulating these parameters reflecting spatial abrupt change characteristics and temporal acceleration characteristics into the same high-dimensional mathematical matrix, a three-dimensional deformation tensor feature space can be constructed to comprehensively and multidimensionally characterize the evolution of disaster risks. For example, when a landslide risk point generates a shear gradient in space and exhibits accelerated subsidence in time, this tensor space can more comprehensively characterize its dangerous situation.

[0081] S60. Based on the three-dimensional deformation tensor feature space, determine the disaster probability prediction results.

[0082] Specifically, determining the disaster probability prediction result refers to using intelligent recognition algorithms to perform pattern matching and risk assessment on the previously constructed feature tensors. By comprehensively weighing the abnormal fracture features in the spatial gradient and the abrupt landslide signal in the temporal acceleration, the risk level of geological disasters at each monitoring point in the target area can be quantitatively scored, generating a confidence score or probability distribution map representing the severity of the hidden danger at that location. This prediction method based on multidimensional tensor space can effectively filter out false deformations caused by environmental interference, thereby outputting accurate identification conclusions for small-scale landslides, collapses and other hidden dangers.

[0083] S70. Based on the disaster probability prediction results, construct a visual disaster probability prediction product view and present the disaster probability prediction product view through the display interface of the disaster probability prediction platform.

[0084] Specifically, constructing a visualized disaster probability prediction product view refers to the multi-level fusion rendering of complex mathematical prediction probabilities with the three-dimensional geographic environment of the target area. By dynamically overlaying the disaster probability prediction results onto the real terrain model in the form of color heat maps or risk labels, and displaying them in an interactive interface, relevant personnel can intuitively and quickly locate the specific latitude and longitude, slope risk, and threatened range of high-risk hazard points, thereby achieving intuitive inspection and accurate location of disaster hazards.

[0085] In one embodiment, step S20, which involves enhancing the original polarization interferometric data to obtain enhanced polarization interferometric data, specifically includes:

[0086] S21. Calculate the signal arrival time difference between each array element in the dual-band SAR receiving antenna, and perform phase weighting compensation on the original polarization interferometric data based on the signal arrival time difference.

[0087] Specifically, the signal arrival time difference between each array element refers to the time delay determined by the antenna element spacing d (in meters) and the target incident angle θ (in degrees). Where c is the speed of light, and its value is 3 × 10⁻⁶. 8 The unit is m / s; the raw polarization interferometric data S received by each array element i (t) multiplied by the compensated phase weight , where f is the carrier frequency (in Hz), which enables the echo signal from the target direction to be aligned in phase, thereby achieving constructive interference to enhance the target signal strength and suppress clutter from non-target directions.

[0088] S22. The original polarization interferometric data after phase weighting compensation is subjected to noise reduction processing to obtain high signal-to-noise ratio polarization interferometric data.

[0089] Specifically, noise reduction refers to smoothing the phase-weighted complex signal using a mean filter operator. By calculating the coherence coefficient γ (ranging from 0 to 1) of pixels within a local window, such as a 5×5 window, and removing noise components with coherence values ​​below a preset threshold, the phase jump caused by random thermal noise in the echo can be significantly reduced. Thus, high signal-to-noise ratio polarized interferometric data with stronger phase continuity can be obtained while preserving the edge details of ground object interference fringes.

[0090] S23. Range-compressed data is obtained by using the ω-k algorithm to perform range compression on high signal-to-noise ratio polarized interferometric data.

[0091] Specifically, range compression using the ω-k algorithm refers to performing coordinate mapping and matched filtering on radar echoes in the wavenumber domain, through the mapping relationship. Resampling is performed on high signal-to-noise ratio polarimetric interferometric data, where k x and k y These are the azimuth wavenumber and range wavenumber (in rad / m), respectively, and ω is the angular frequency (in rad / s). This processing can correct the geometric coupling between range and azimuth caused by large squint observations, and refocus the pulse energy to the corresponding range sampling unit, thereby producing range compressed data.

[0092] S24. Based on the phase gradient estimation method, the range-compressed data is subjected to azimuth self-focusing processing to obtain enhanced polarization interferometric data.

[0093] Specifically, azimuth self-focusing based on the phase gradient estimation method refers to extracting the phase gradient $ΔΦ (in rad) of the signal in the azimuth direction, and calculating the residual phase error caused by motion instability using the least squares criterion. And perform phase correction operation on the range-compressed data. This eliminates the defocusing caused by flight platform attitude disturbances, ultimately obtaining enhanced polarization interferometric data with focusing accuracy reaching the theoretical limit.

[0094] In one embodiment, step S30 involves performing terrain-adaptive radiometric correction on the enhanced polarimetric interferometric data, and then using a terrain-constrained quadratic surface fitting algorithm to perform image registration on the terrain-adaptive radiometric corrected enhanced polarimetric interferometric data to obtain preprocessed interferometric image data. Specifically, this includes:

[0095] S31. Obtain the external DEM model data of the target area, and calculate the initial geometric offset of the auxiliary image relative to the main image in the enhanced polarimetric interferometry data based on the external DEM model data.

[0096] Specifically, external DEM model data refers to a gridded digital elevation matrix H with units of 30m. By mapping the geographic coordinates of the main and auxiliary images to the WGS84 coordinate system and solving the instantaneous satellite orbit vector, the initial geometric offset vector (Δx0, Δy0) between the two images caused by terrain undulation and baseline length B (in meters) can be calculated. This vector is measured in pixels and provides a reliable search reference point for subsequent refined sub-pixel registration.

[0097] S32. Using the DEM model of the target area, calculate the local slope and azimuth of each sampling point in the DEM model.

[0098] Specifically, calculating the local slope α (in °) and azimuth β (in °) refers to performing a first-order spatial difference operation on the elevation points in the DEM model, i.e., based on the formula... Determining the tilt of ground sampling points, these topographic geometric parameters reflect the real-time attitude of the ground surface relative to the direction of radar wave illumination, and are the core basis for subsequent radiation accuracy compensation.

[0099] S33. Based on the relationship between local slope, azimuth angle and satellite incident angle, calculate the radiation normalization coefficient. Use the radiation normalization coefficient to process the enhanced polarization interferometric data to eliminate the difference in terrain brightness and obtain terrain adaptive radiation correction data.

[0100] Specifically, the radiation normalization coefficient K is determined by the local incident angle θ. loc The dimensionless correction factor, determined by the cosine value (in degrees), is calculated by performing an operation. Weighted adjustment of image brightness values ​​can counteract overbrightness caused by hillsides facing the radar beam or underbrightness caused by hillsides facing away from the beam, so that the corrected terrain adaptive radiometric correction data can more realistically reflect the electromagnetic scattering characteristics of the ground objects themselves.

[0101] S34. Determine the search window in the cross-correlation peak region of the main image and auxiliary image using the initial geometric offset, and construct a quadratic surface equation about the cross-correlation coefficient through least squares fitting.

[0102] Specifically, constructing the quadratic surface equation involves selecting a pixel window within the neighborhood of the cross-correlation peaks after initial alignment, collecting the cross-correlation coefficient ρ (ranging from 0 to 1) of each sampling point as observation samples, and fitting a smooth mathematical surface model using the least squares method. This model can characterize the geometric similarity distribution of the main and auxiliary images at the sub-pixel level as a continuous function.

[0103] S35. Based on the quadratic surface equation, determine the offset vector and use the offset vector to perform image resampling and alignment so that the auxiliary image is aligned with the main image, thus obtaining preprocessed interferometric image data.

[0104] Specifically, determining the offset vector involves performing differentiation on the quadratic surface equation and solving the system of equations. This allows for the precise sub-pixel coordinates (x) of the point where the cross-correlation coefficient is maximized. peak ,y peak The offset vector determined in this way can achieve an accuracy of 0.05 pixels. Then, the Singer function is used to perform pixel resampling interpolation on the auxiliary image to ensure that each ground point in the two images is perfectly aligned in coordinates, thereby obtaining preprocessed interferometric image data that meets the requirements of subsequent interferometric calculations.

[0105] In one embodiment, step S40, which involves performing phase unwrapping on the preprocessed interferometric image data to obtain initial deformation field data, specifically includes:

[0106] S41. Select image pairs whose time baseline and vertical baseline are both within the corresponding preset thresholds, and construct an interferogram stack consisting of at least a preset number of preprocessed interferometric image data.

[0107] Specifically, the screening baseline refers to the decision made by the execution logic that the time baseline is within a preset threshold, i.e., B. temp ≤90 days and all vertical baselines are within the preset threshold, i.e., B perp For images with a depth of ≤300m, by selecting digital image sequences that meet the above conditions, it is possible to effectively suppress decoherent noise caused by seasonal growth of surface vegetation or excessive satellite orbit span, thereby constructing a digital interferogram stack containing at least 30 images with stable phase quality.

[0108] S42. Divide the phase map in the interferogram stack into multiple sub-blocks, and process the sub-blocks in parallel based on the statistical cost flow unwrapping algorithm to obtain continuous phase stack data.

[0109] Specifically, dividing the phase map into multiple sub-blocks allows for the simultaneous processing of different local regions of a wide-swath image using parallel computing resources. Within each sub-block, a statistical cost flow algorithm is applied to find the optimal path for phase integration; that is, a network flow model is constructed and the cost function is minimized. , where w ij The weight is ΔΦ, which is the phase difference between adjacent pixels. The wrapped phase is restored to an absolute phase value (in rad) that is not limited by the 2π period, thereby obtaining continuous phase stack data that reflects the continuous displacement trend of the earth's surface.

[0110] S43. Based on continuous phase stack data, construct deformation observation equations. Use singular value decomposition to solve the generalized inverse of the structure matrix of the deformation observation equations to obtain the initial values ​​of deformation rate at each monitoring point.

[0111] Specifically, inverse solution using singular value decomposition refers to solving the deformation observation equations. Decompose the time interval matrix A into UΣV T Through calculation The initial value of the deformation rate V (in millimeters per year) for each monitoring point is solved, so that even if some observation data is missing, the initial value of the deformation rate of each monitoring point throughout the entire observation period can be robustly calculated.

[0112] S44. Use a preset turbulent atmospheric model to identify the atmospheric delay phase in the initial value of deformation rate, and remove the atmospheric delay phase from the initial value of deformation rate to obtain the initial deformation field data.

[0113] Specifically, the atmospheric random phase component Φ with specific spatial autocorrelation characteristics in the initial values ​​is extracted using a turbulent atmospheric model.atm (Unit: rad), by performing subtraction. By removing the delay error caused by tropospheric water vapor fluctuations, initial deformation field data that can truly reflect millimeter-level displacement changes on the Earth's surface can be obtained.

[0114] In one embodiment, step S50, namely, determining the spatial strain gradient and temporal evolution acceleration based on the initial deformation field data, and constructing a three-dimensional deformation tensor feature space based on the spatial strain gradient and temporal evolution acceleration, specifically includes:

[0115] S51. Extract the cumulative line-of-sight deformation and deformation time series of each monitoring point in the initial deformation field data.

[0116] Specifically, extracting cumulative deformation and sequence refers to reading the line-of-sight displacement value L(t) (in millimeters, mm) of each sampling point at different time points t, forming an array sequence that characterizes the displacement evolution trajectory of that point from the monitoring start date to the current date. This digital sequence records the displacement evolution process of the slope over several months or even several years.

[0117] S52. Based on the cumulative deformation along the line of sight, calculate the first-order spatial partial derivative between the target monitoring point and the neighboring monitoring points to obtain the spatial strain gradient.

[0118] Specifically, calculating the first-order spatial partial derivative refers to performing a difference operation. By calculating the difference in displacement between the target point and its neighboring points within a unit distance (in meters), the spatial strain gradient (in mm / m) reflecting the characteristics of abrupt changes on the ground surface is obtained. This index can accurately identify the discontinuous displacement areas caused by landslide boundaries.

[0119] S53. Perform a second-order difference operation on the time series to calculate the time series evolution acceleration.

[0120] Specifically, performing second-order difference operations refers to using the formula Calculate the second-order rate of change of displacement over time. For example, if the displacement of a certain slope is 2mm, 5mm, and 15mm over three consecutive months, the temporal evolution acceleration (in mm / d²) extracted by this calculation can reflect whether the surface movement is in a dangerous critical period of transition from constant creep to exponential acceleration.

[0121] S54. By mapping the line of sight to the cumulative deformation, spatial strain gradient, and temporal evolution acceleration, a three-dimensional deformation tensor feature space is obtained.

[0122] Specifically, vector mapping refers to mapping the cumulative deformation L and spatial strain gradient of each monitoring point. And the temporal evolution acceleration 'a' is combined into an eigenvector. This data is then mapped to a predefined three-dimensional mathematical vector space. By constructing this carrier that integrates multi-dimensional physical features, a structurally unified feature data source is provided for subsequent recognition algorithms.

[0123] In one embodiment, step S60, which determines the disaster probability prediction result based on the three-dimensional deformation tensor feature space, specifically includes:

[0124] S61. Input the three-dimensional deformation tensor feature space into the pre-trained ResNet-50 convolutional neural network, and obtain the spatial morphological features of the disaster-prone area through the residual convolutional layer.

[0125] Specifically, extracting spatial morphological features refers to using convolutional kernels to perform multi-layer dimensionality reduction and feature mapping on the input three-dimensional deformation tensor, and learning the implicit spatial topological relationships in the deformation data through the residual network y=F(x)+x, thereby automatically identifying spatial geometric morphological features that are highly correlated with landslide hazards.

[0126] S62. Based on spatial morphological features, the temporal evolution sequence in the three-dimensional deformation tensor is weighted using an attention mechanism to obtain a temporal fusion feature vector.

[0127] Specifically, using the attention mechanism to assign weights is the basis for calculating scores. Different importance weight factors are assigned to the deformation characteristics at different observation times, so that the model can prioritize the key time nodes of acceleration mutation and generate a time-series fusion feature vector that reflects the abnormal evolution trend.

[0128] S63. Input the time-series fusion feature vector into a random forest classifier for nonlinear classification processing to output the disaster probability prediction result for each monitoring point.

[0129] Specifically, nonlinear classification using a random forest classifier refers to constructing an ensemble learning model composed of multiple decision trees, and calculating the confidence score (ranging from 0% to 100%) of each monitoring point for potential disaster risks based on the voting results of each feature component under different decision paths, as the result of disaster probability prediction.

[0130] In one embodiment, step S70, namely, constructing a visualized disaster probability prediction product view based on the disaster probability prediction results and presenting the disaster probability prediction product view through the display interface of the disaster probability prediction platform, specifically includes:

[0131] S71. Load the 3D point cloud data of the target area using the WebGL engine, and render a heat map on the 3D point cloud data based on the disaster probability prediction results.

[0132] Specifically, rendering a heatmap refers to mapping probability scores to specific color mapping values ​​(such as red for high risk) based on a digital shader program, and merging color attributes into a 3D point cloud model using a volume rendering algorithm, thereby achieving intuitive visual identification of high-risk areas in a digital twin scene.

[0133] S72. Based on the geographical range corresponding to the heat map, use GIS analysis tools to calculate the slope parameters of disaster hazard points and the range of threatened watersheds.

[0134] Specifically, calculating slope parameters and threatened range refers to conducting spatial analysis on the area covered by the heat map in the digital geographic information system, extracting the maximum tangent plane angle of the DEM within this range and performing digital hydrological runoff simulation calculations to automatically define the affected range of the digital watershed that may be affected after the hazard becomes unstable.

[0135] S73. Generate early warning reports based on the three-dimensional geographic coordinates, slope parameters, and threatened watershed range of disaster hazard points.

[0136] Specifically, generating an early warning report refers to automatically extracting key attributes such as predicted probability distribution, latitude and longitude coordinates, average slope, and digital disaster area, and packaging them into a digital report file containing statistical charts and early warning conclusions according to a preset standard format. For example, generating XML or JSON format documents that conform to the CAP protocol standard, aiming to provide emergency management departments with directly usable digital decision-making basis.

[0137] S74. Based on preset layer display rules, perform layer synthesis processing on the early warning report to generate a disaster probability prediction product view.

[0138] Specifically, layer compositing refers to superimposing the 3D terrain model layer, thermal cloud layer, and digital early warning briefing layer in a multi-dimensional space according to the preset digital layer overlay logic. This ensures that information at different levels can be clearly and distinctly presented in the platform's display interface, thereby generating a highly integrated digital disaster probability prediction product view.

[0139] It should be understood that the sequence number of each step in the above embodiments does not imply the order of execution. The execution order of each process should be determined by its function and internal logic, and should not constitute any limitation on the implementation process of the embodiments of this application.

[0140] In one embodiment, an InSAR rapid identification and precise location system for disaster hazards is provided, which corresponds one-to-one with the InSAR rapid identification and precise location method for disaster hazards described in the above embodiment. For example... Figure 2As shown, this InSAR rapid identification and precise positioning system includes a data acquisition module, an enhancement processing module, a data preprocessing module, a deformation field extraction module, a feature space construction module, a probability prediction module, and a visualization module. Detailed descriptions of each functional module are as follows:

[0141] The data acquisition module is used to acquire the raw polarimetric interferometric data of the target area through a dual-band SAR receiving antenna;

[0142] The enhancement processing module is used to enhance the original polarization interference data to obtain enhanced polarization interference data;

[0143] The data preprocessing module is used to perform terrain-adaptive radiometric correction on the enhanced polarimetric interferometric data, and to perform image registration on the enhanced polarimetric interferometric data after terrain-adaptive radiometric correction using a terrain-constrained quadratic surface fitting algorithm to obtain preprocessed interferometric image data.

[0144] The deformation field extraction module is used to perform phase unwrapping processing on the preprocessed interferometric image data to obtain initial deformation field data.

[0145] The feature space construction module is used to determine the spatial strain gradient and temporal evolution acceleration based on the initial deformation field data, and to construct a three-dimensional deformation tensor feature space according to the spatial strain gradient and the temporal evolution acceleration.

[0146] The probability prediction module is used to determine the disaster probability prediction result based on the three-dimensional deformation tensor feature space.

[0147] The visualization module is used to construct a visualized disaster probability prediction product view based on the disaster probability prediction results, and to present the disaster probability prediction product view through the display interface of the disaster probability prediction platform.

[0148] Specific limitations regarding the InSAR rapid identification and precise location system for disaster hazards can be found in the limitations of the InSAR rapid identification and precise location method for disaster hazards described above, and will not be repeated here. Each module in the aforementioned InSAR rapid identification and precise location system for disaster hazards can be implemented entirely or partially through software, hardware, or a combination thereof. These modules can be embedded in or independent of the processor in a computer device, or stored in the memory of a computer device as software, so that the processor can call and execute the corresponding operations of each module.

[0149] Those skilled in the art will clearly understand that, for the sake of convenience and brevity, the above-described division of functional units and modules is used as an example. In practical applications, the above functions can be assigned to different functional units and modules as needed, that is, the internal structure of the device can be divided into different functional units or modules to complete all or part of the functions described above.

[0150] The above-described embodiments are only used to illustrate the technical solutions of this application, and are not intended to limit them. Although this application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of this application, and should all be included within the protection scope of this application.

Claims

1. A method for rapid InSAR identification and precise location of potential disaster hazards, characterized in that, The InSAR rapid identification and accurate positioning method includes: The raw polarization interferometric data of the target area is obtained by using a dual-band SAR receiving antenna; The original polarimetric interferometric data is enhanced to obtain enhanced polarimetric interferometric data; The enhanced polarimetric interferometric data is subjected to terrain-adaptive radiometric correction, and the enhanced polarimetric interferometric data after terrain-constrained quadratic surface fitting algorithm is used for image registration to obtain preprocessed interferometric image data. The preprocessed interferometric image data is subjected to phase unwrapping to obtain initial deformation field data; Based on the initial deformation field data, the spatial strain gradient and temporal evolution acceleration are determined, and a three-dimensional deformation tensor feature space is constructed according to the spatial strain gradient and the temporal evolution acceleration. Based on the aforementioned three-dimensional deformation tensor feature space, the disaster probability prediction results are determined; Based on the disaster probability prediction results, a visual disaster probability prediction product view is constructed and presented through the display interface of the disaster probability prediction platform.

2. The InSAR rapid identification and accurate positioning method according to claim 1, characterized in that, The enhancement processing of the original polarization interferometric data to obtain enhanced polarization interferometric data specifically includes: Calculate the signal arrival time difference between each array element in the dual-band SAR receiving antenna, and perform phase weighting compensation on the original polarization interferometric data based on the signal arrival time difference; The original polarimetric interferometric data after phase weighting compensation is denoised to obtain high signal-to-noise ratio polarimetric interferometric data. The high signal-to-noise ratio polarized interferometric data is compressed in the range direction using the ω-k algorithm to obtain the range-compressed data. The range-compressed data is subjected to azimuth self-focusing processing based on the phase gradient estimation method to obtain the enhanced polarization interferometry data.

3. The InSAR rapid identification and accurate positioning method according to claim 2, characterized in that, The process of performing terrain-adaptive radiometric correction on the enhanced polarimetric interferometric data and then using a terrain-constrained quadratic surface fitting algorithm to perform image registration on the terrain-adaptive radiometric correction-corrected enhanced polarimetric interferometric data to obtain preprocessed interferometric image data specifically includes: Obtain the external DEM model data of the target area, and calculate the initial geometric offset of the auxiliary image relative to the main image in the enhanced polarimetric interferometry data based on the external DEM model data; Using the DEM model of the target area, the local slope and azimuth of each sampling point in the DEM model are calculated; Based on the relationship between the local slope, azimuth angle and satellite incident angle, the radiation normalization coefficient is calculated. The enhanced polarization interferometric data is then processed using the radiation normalization coefficient to eliminate the difference in terrain brightness, thereby obtaining terrain adaptive radiometric correction data. The search window is determined using the initial geometric offset in the cross-correlation peak region of the main image and the auxiliary image, and a quadratic surface equation about the cross-correlation coefficient is constructed by least squares fitting. Based on the quadratic surface equation, an offset vector is determined, and the offset vector is used for image resampling and alignment to align the auxiliary image with the main image, thereby obtaining the preprocessed interferometric image data.

4. The InSAR rapid identification and accurate positioning method according to claim 1, characterized in that, The step of performing phase unwrapping processing on the preprocessed interferometric image data to obtain initial deformation field data specifically includes: Image pairs whose temporal and vertical baselines are both within the corresponding preset thresholds are selected, and an interferogram stack consisting of at least a preset number of the preprocessed interferometric image data is constructed. The phase map in the interferogram stack is divided into multiple sub-blocks, and the sub-blocks are processed in parallel based on the statistical cost flow unwrapping algorithm to obtain continuous phase stack data. Based on the continuous phase stack data, a deformation observation equation is constructed. The structure matrix of the deformation observation equation is solved by generalized inverse solution using singular value decomposition to obtain the initial value of the deformation rate of each monitoring point. The atmospheric delay phase in the initial value of the deformation rate is identified using a preset turbulent atmospheric model, and the atmospheric delay phase is removed from the initial value of the deformation rate to obtain the initial deformation field data.

5. The InSAR rapid identification and accurate positioning method according to claim 1, characterized in that, The process of determining the spatial strain gradient and temporal evolution acceleration based on the initial deformation field data, and constructing a three-dimensional deformation tensor feature space based on the spatial strain gradient and the temporal evolution acceleration, specifically includes: Extract the line-of-sight cumulative deformation and deformation time series of each monitoring point from the initial deformation field data; Based on the cumulative deformation along the line of sight, the first spatial partial derivative between the target monitoring point and the neighboring monitoring points is calculated to obtain the spatial strain gradient. Perform a second-order difference operation on the deformation time series to calculate the time series evolution acceleration; The three-dimensional deformation tensor feature space is obtained by vector mapping the cumulative deformation along the line of sight, the spatial strain gradient, and the temporal evolution acceleration.

6. The InSAR rapid identification and accurate positioning method according to claim 1, characterized in that, The determination of disaster probability prediction results based on the three-dimensional deformation tensor feature space specifically includes: The three-dimensional deformation tensor feature space is input into a pre-trained ResNet-50 convolutional neural network, and the spatial morphological features of the disaster-prone area are obtained through residual convolutional layers. Based on the spatial morphological features, an attention mechanism is used to assign weights to the temporal evolution sequence in the three-dimensional deformation tensor to obtain a temporal fusion feature vector. The time-series fusion feature vector is input into a random forest classifier for nonlinear classification processing to output the disaster probability prediction result for each monitoring point.

7. The InSAR rapid identification and accurate positioning method according to claim 1, characterized in that, The step of constructing a visualized disaster probability prediction product view based on the disaster probability prediction results and presenting the disaster probability prediction product view through the display interface of the disaster probability prediction platform specifically includes: The WebGL engine is used to load the three-dimensional point cloud data of the target area, and a heat map layer is rendered on the three-dimensional point cloud data according to the disaster probability prediction results. Based on the geographical range corresponding to the heat map, GIS analysis tools are used to calculate the slope parameters of disaster hazard points and the range of threatened watersheds. Based on the three-dimensional geographic coordinates of the disaster hazard points, the slope parameters, and the threatened watershed range, an early warning report is generated; Based on preset layer display rules, the early warning report is processed by layer synthesis to generate the disaster probability prediction product view.

8. A rapid InSAR system for identifying and accurately locating potential disaster hazards, characterized in that, The InSAR rapid identification and precise location method for disaster hazards as described in claims 1-7, wherein the InSAR rapid identification and precise location system comprises: The data acquisition module is used to acquire the raw polarimetric interferometric data of the target area through a dual-band SAR receiving antenna; The enhancement processing module is used to enhance the original polarization interference data to obtain enhanced polarization interference data; The data preprocessing module is used to perform terrain-adaptive radiometric correction on the enhanced polarimetric interferometric data, and to perform image registration on the enhanced polarimetric interferometric data after terrain-adaptive radiometric correction using a terrain-constrained quadratic surface fitting algorithm to obtain preprocessed interferometric image data. The deformation field extraction module is used to perform phase unwrapping processing on the preprocessed interferometric image data to obtain initial deformation field data. The feature space construction module is used to determine the spatial strain gradient and temporal evolution acceleration based on the initial deformation field data, and to construct a three-dimensional deformation tensor feature space according to the spatial strain gradient and the temporal evolution acceleration. The probability prediction module is used to determine the disaster probability prediction result based on the three-dimensional deformation tensor feature space. The visualization module is used to construct a visualized disaster probability prediction product view based on the disaster probability prediction results, and to present the disaster probability prediction product view through the display interface of the disaster probability prediction platform.

Images