Intelligent data processing method and system for geological disaster early warning

Abstract

The invention discloses an intelligent data processing method and system for geological disaster early warning, and relates to geological disaster early warning, and the method comprises the steps: collecting multi-source heterogeneous data; extracting topographic features of the geological disaster area according to the optical satellite image data and the laser radar data; performing differential interference processing on the synthetic aperture radar data to obtain differential interference phase data of the geological disaster area; performing phase unwrapping by using an improved minimum cost flow algorithm to obtain a continuous absolute phase value # imgabs0 #; performing nonlinear estimation on the estimation result of the surface deformation quantity to obtain state parameter estimation of the geological disaster area; and performing early warning detection by using a support vector machine (SVM). In order to solve the problem of low geological disaster prediction precision caused by difficult phase unwrapping in the prior art, the unwrapping problem is converted into a space-time network flow graph, space-time continuity constraints are introduced, the incidence relation between data in different periods is considered, a multi-scale strategy is adopted for solving and the like, and the prediction precision is improved.

Description

Technical Field

[0001] The present application relates to the field of geological disaster early warning, and in particular to an intelligent data processing method and system for geological disaster early warning. Background Art

[0002] With the impact of global climate change and human engineering activities, the frequency and intensity of geological disasters are increasing, posing a serious threat to the safety of people's lives and property and social and economic development. Timely and accurate monitoring and early warning of geological disasters have become an important issue in the field of disaster prevention and mitigation.

[0003] However, the formation mechanism of geological disasters is complex, involving multiple factors such as geological structure, topography, hydrology and meteorology. Traditional single monitoring methods are difficult to meet the increasingly severe disaster prevention needs. Synthetic aperture radar interferometry (InSAR) technology has shown great potential in surface deformation monitoring and geological disaster early warning due to its large range, high precision, and all-day and all-weather characteristics. Differential interferometric synthetic aperture radar (DInSAR) technology can extract small displacement change information on the surface by performing differential interferometry processing on radar data acquired at different times, and has become one of the important means of geological disaster monitoring and early warning.

[0004] However, in practical applications, DInSAR technology still faces many challenges. First, when the time interval between two radar data acquisitions is long, factors such as surface cover changes and vegetation growth will cause the temporal decorrelation of radar signals, affecting the quality of the differential interferometer phase, thereby reducing the accuracy of displacement estimation. Secondly, the differential interferometer phase is a fuzzy phase with a period of 2π, and the absolute phase needs to be restored through a phase unwrapping algorithm. However, the phase unwrapping algorithm has certain limitations and is prone to unwrapping errors, especially in areas with large displacement gradients such as faults and landslides, which affects the reliability of displacement estimation. In addition, DInSAR technology obtains the displacement component in the direction of the radar line of sight. To convert it into a real three-dimensional surface displacement vector, it is also necessary to consider the influence of factors such as radar imaging geometry and terrain undulation. Summary of the invention

[0005] In view of the low accuracy of geological disaster prediction caused by the difficulty of phase unwrapping in the prior art, the present application provides an intelligent data processing method and system for geological disaster early warning, which converts the unwrapping problem into a space-time network flow graph, introduces space-time continuity constraints, considers the correlation between data in different periods, and adopts a multi-scale strategy for solving, thereby improving the prediction accuracy.

[0006] The purpose of this application is achieved through the following technical solutions.

[0007] One aspect of the present application provides an intelligent data processing method for geological disaster early warning, including: S1, collecting multi-source heterogeneous data, wherein the multi-source heterogeneous data includes synthetic aperture radar data, optical satellite image data and laser radar data; S2, extracting terrain features of the geological disaster area based on the optical satellite image data and the laser radar data, wherein the terrain features include slope, slope aspect and elevation; S3, performing differential interference processing on the synthetic aperture radar data to obtain differential interference phase data of the geological disaster area; S4, performing phase unwrapping based on the differential interference phase data using an improved minimum cost flow algorithm to obtain a continuous absolute phase value φ; the improved minimum cost flow algorithm introduces spatiotemporal continuity constraints and considers the correlation between differential interference phase data of different periods; converting the absolute phase value φ into a deformation variable in the line of sight direction, and estimating the surface deformation variable of the geological disaster area based on the geometric parameters of the synthetic aperture radar data; S5, performing nonlinear estimation on the surface deformation variable estimation result to obtain a state parameter estimation of the geological disaster area; S6, using the state parameter estimation as input and performing early warning detection using a support vector machine SVM.

[0008] Further, S4, based on the differential interferometry phase data, the phase unwrapping is performed using the improved minimum cost flow algorithm to obtain continuous absolute phase values. , including: S41, converting the differential interference phase data into a space-time network flow graph, each node in the space-time network flow graph represents a pixel point at a space-time position, and the phase difference between adjacent pixels in space and the phase change between adjacent pixels in time are represented as edge costs; S42, setting source nodes and sink nodes in the space-time network flow graph, and setting a flow equal to the total number of pixels of the differential interference phase data; S43, performing multi-scale decomposition on the differential interference phase data to obtain phase graphs at different spatial scales; S44, at the coarsest scale, obtaining the minimum cost path in the space-time network flow graph through the minimum cost maximum flow algorithm as the initial untangling result; S45, at each scale, using the untangling result of the previous scale as the initial value, and adjusting the edge cost according to the terrain characteristics and space-time continuity constraints, and obtaining the untangling result of the current scale through the minimum cost maximum flow algorithm; S46, fusing the untangling results of each scale to obtain the final absolute phase value .

[0009] Among them, geological disasters usually occur in areas with complex and changeable terrain, such as mountainous areas and hills. Complex terrain can cause a large number of phase jumps and nonlinear changes in differential interferometry phase data. The minimum cost flow algorithm is based on the assumption of linearization, and it is difficult to accurately model the phase changes caused by complex terrain, which may lead to unwrapping errors. The occurrence and development of geological disasters have spatiotemporal correlations, and there is a certain correlation between differential interferometry phase data in different periods. The traditional minimum cost flow algorithm only considers data from a single period and does not make full use of time series information, which may lead to poor temporal consistency of the unwrapping results.

[0010] In this application, on the one hand, when constructing a network flow graph, in addition to considering the phase difference between adjacent pixels as the cost of the edge, terrain information is also introduced for adaptive adjustment. For areas with complex terrain, the weight of the phase difference is increased to emphasize the continuity of the phase; for areas with flat terrain, the weight of the phase difference is reduced to allow more phase jumps. Through adaptive edge cost calculation, phase untangling under different terrain conditions can be better adapted.

[0011] On the other hand, the spatiotemporal correlation of geological hazards is used to introduce spatiotemporal continuity constraints in the network flow graph. The differential interferometric phase data of different periods are combined to construct a spatiotemporal network flow graph, in which each node represents a pixel point at a spatiotemporal position. When calculating the cost of the edge, not only the phase difference between spatially adjacent pixels but also the phase change between temporally adjacent pixels is considered. By introducing spatiotemporal continuity constraints, the consistency of the disentanglement results in the time dimension can be improved.

[0012] Further, S42, setting source nodes and sink nodes in the space-time network flow graph, and setting a flow equal to the total number of pixels of the differential interference phase data, including: in the space-time network flow graph, setting the first pixel point of the first time step as the source node, and setting the last pixel point of the last time step as the sink node; counting the total number of pixels N of the differential interference phase data, using N as the flow of the source node, and using -N as the flow of the sink node; setting the flow of other nodes except the source node and the sink node to zero.

[0013] Further, S45, at each scale, the disentanglement result of the previous scale is used as the initial value, and the cost of the edge is adjusted according to the terrain characteristics and the space-time continuity constraint, and the disentanglement result of the current scale is obtained by the minimum cost maximum flow algorithm, including: interpolating the disentanglement result of the previous scale to the current scale as the initial value of the disentanglement of the current scale; according to the terrain characteristics, calculating the local slope value at each node, dividing the area where the node is located into a first area and a second area according to a preset threshold, the first area represents a complex terrain area, and the second area represents a flat terrain area; for the nodes in the first area, increasing The proportion of the phase difference on the edge between the node and the adjacent node in the edge cost is reduced for the nodes in the second area; the phase change between each node and the adjacent node in the time series in the spatiotemporal network flow graph after the terrain feature constraint is calculated; the phase change is multiplied by the time continuity weight to obtain the time continuity cost, and the time continuity weight is set according to the time baseline length; the time continuity cost is superimposed on the cost of the corresponding edge to obtain the adjusted edge cost; the minimum cost maximum flow algorithm is used to obtain the minimum cost path in the spatiotemporal network flow graph after the edge cost is adjusted as the untangling result of the current scale.

[0014] Furthermore, the time continuity weight is set according to the following formula: , where T represents the length of the time baseline, α controls the overall magnitude of the time continuity weight, and β controls the rate at which the time continuity weight increases as the time baseline length increases. In this application, the time continuity weight It is logarithmically related to the time baseline length T. When the time baseline length is short, the time continuity weight is small, and the phase change between adjacent time points has little effect on the edge cost; when the time baseline length is long, the time continuity weight is large, and the phase change between adjacent time points has a greater impact on the edge cost, so as to ensure the continuity and consistency of the unwrapping results on the long time series. By introducing the preset parameters α and β, the changing trend of the time continuity weight can be flexibly adjusted to adapt to different time baseline lengths and application scenarios. Among them, α controls the overall amplitude of the time continuity weight, and β controls the rate at which the time continuity weight increases with the increase of the time baseline length.

[0015] Furthermore, the absolute phase value The deformation variable in the line of sight is converted into the deformation variable in the line of sight direction, and the surface deformation variable of the geological disaster area is estimated according to the geometric parameters of the synthetic aperture radar data, including: according to the wavelength λ of the synthetic aperture radar data, the absolute phase value φ is converted into the deformation variable d_los in the line of sight direction, and the conversion formula is: ; Calculate the incident angle θ and azimuth α of the synthetic aperture radar signal; According to the incident angle θ and azimuth α, the deformation in the line of sight direction Convert to three-dimensional surface shape ; Combined with the topographic characteristics of the geological disaster area, the least squares equations are constructed to solve the three-dimensional surface deformation variables ; Using the Kriging interpolation algorithm, the three-dimensional surface deformation variables are spatially interpolated to obtain the three-dimensional surface deformation field in the geological disaster area.

[0016] Furthermore, the incident angle θ and azimuth angle α of the synthetic aperture radar signal are calculated, including: extracting orbital parameters of the satellite platform from the synthetic aperture radar data, the orbital parameters including satellite position, velocity and attitude; calculating the geometric relationship between the satellite platform and the geological disaster area based on the orbital parameters and the geographic coordinates of the geological disaster area; and calculating the incident angle θ and azimuth angle α of the synthetic aperture radar signal in the geological disaster area based on the geometric relationship.

[0017] Furthermore, the deformation in the sight direction is calculated by the following formula: Convert to three-dimensional surface shape : ,in, It represents the deformation in the line of sight direction and the surface displacement component in the direction of radar signal propagation, with the unit of meter (m). Represents the east-west (x-axis) surface deformation variable, indicating the displacement component of the surface in the east-west direction, in meters (m). Positive values ​​indicate movement to the east, and negative values ​​indicate movement to the west. Represents the surface deformation in the north-south direction (y-axis), indicating the displacement component of the surface in the north-south direction, in meters (m). A positive value indicates movement to the north, and a negative value indicates movement to the south. Represents the vertical (z-axis) surface deformation, which represents the displacement component of the surface in the vertical direction, in meters (m). A positive value indicates upward movement, and a negative value indicates downward movement. α represents the azimuth of the radar signal, which represents the angle between the projection direction of the radar signal on the horizontal plane and the true north direction, in radians (rad). The azimuth angle ranges from [0, 2π]. θ represents the incident angle of the radar signal, which represents the angle between the incident direction of the radar signal and the direction of the surface normal, in radians (rad). The incident angle ranges from (0, π / 2). β represents the angle between the surface deformation direction and the radar signal incident plane, which represents the angle between the projection of the surface deformation direction on the radar signal incident plane and the line of sight, in radians (rad). The β ranges from [0, π].

[0018] This application introduces the angle β between the surface deformation direction and the radar signal incident plane to take into account the situation that the deformation direction of the geological disaster area may not be completely consistent with the radar line of sight. When β=0, it means that the deformation direction is completely consistent with the line of sight; when β=π / 2, it means that the deformation direction is perpendicular to the line of sight. and ) is combined into the horizontal contribution, that is, , and then multiply by , and obtain its projection component in the line of sight direction. This can more accurately describe the distribution characteristics of surface deformation in the horizontal direction.

[0019] Further, S5, nonlinear estimation is performed on the surface deformation variable estimation result to obtain the state parameter estimation of the geological disaster area, including: S51, taking the surface deformation variable as the observation quantity, and taking the state parameter of the geological disaster area as the state quantity, wherein the state parameter includes the surface deformation rate and the surface deformation acceleration; S52, establishing the nonlinear state equation and observation equation between the observation quantity and the state quantity; S53, initializing the estimated value and covariance matrix of the state parameter; S54, predicting the state parameter at the current moment according to the estimated value and state equation of the state parameter at the previous moment; S55, updating the estimated value of the state parameter at the current moment according to the observation quantity and observation equation at the current moment by using the extended Kalman filter algorithm; S56, repeating steps S53 to S55 until the estimated value of the state parameter at all moments is updated to obtain the state parameter estimation of the geological disaster area.

[0020] The traditional Kalman filter assumes that the system state and the observation equation satisfy a linear relationship, but the actual geological disaster process often has nonlinear characteristics. For example, the movement speed of the landslide body is nonlinearly related to the sliding distance. The traditional Kalman filter is difficult to accurately describe this nonlinear dynamic. Secondly, the Kalman filter assumes that the system noise and the observation noise obey the Gaussian distribution, but the actual geological disaster data may contain non-Gaussian noise, such as outliers, outliers, etc. The traditional Kalman filter is sensitive to these noises, resulting in a decrease in estimation accuracy.

[0021] By establishing nonlinear state equations and observation equations between surface deformation variables and state parameters, this application can more accurately characterize the dynamic evolution of geological disasters, such as the nonlinear relationship between the movement speed of landslides and the sliding distance. Unlike the traditional Kalman filter that assumes a linear model, the nonlinear model of this application can better approximate the actual geological disaster process and improve the accuracy of state estimation.

[0022] This application is aimed at nonlinear state equations and observation equations, and traditional Kalman filtering cannot be directly applied. This application introduces an extended Kalman filter algorithm, which converts nonlinear problems into linear problems for processing by locally linearizing nonlinear models. The extended Kalman filter predicts and updates the state parameters at each moment, realizes the recursive calculation of state estimation, retains the excellent properties of the Kalman filter, and can process nonlinear systems, improving the estimation accuracy and computational efficiency.

[0023] Another aspect of the present application also provides an intelligent data processing system for geological disaster warning, which is used to execute an intelligent data processing method for geological disaster warning of the present application.

[0024] Compared with the prior art, the advantages of this application are:

[0025] When traditional differential interferometric synthetic aperture radar (DInSAR) technology calculates the changes in surface displacement between different periods, when the time interval between two radar data acquisitions is long, factors such as surface cover changes and vegetation growth will cause temporal decorrelation of radar signals, affecting the quality of differential interferometric phase, thereby reducing the accuracy of displacement estimation. This application introduces spatiotemporal continuity constraints, considers the correlation between data from different periods during phase unwrapping, and reduces the impact of temporal decorrelation. At the same time, a multi-source data fusion strategy is used to extract terrain features in combination with optical images and lidar data to provide auxiliary information for InSAR data processing, further improving the quality of differential interferometric phase and the accuracy of displacement estimation.

[0026] The differential interference phase is a fuzzy phase with a period of 2π, and the absolute phase needs to be restored through a phase unwrapping algorithm. However, the phase unwrapping algorithm has certain limitations and is prone to unwrapping errors, especially in areas with large displacement gradients such as faults and landslides, which affects the reliability of displacement estimation. This application uses an improved minimum cost flow algorithm to perform phase unwrapping, transforms the problem into a spatiotemporal network flow graph, introduces terrain feature constraints and spatiotemporal continuity constraints, and adopts a multi-scale unwrapping strategy to effectively overcome the limitations of traditional phase unwrapping algorithms. In areas with large displacement gradients such as faults and landslides, the cost of the edge is adjusted through terrain feature constraints to reduce unwrapping errors; the spatiotemporal continuity constraints are used to consider the phase change relationship between adjacent pixels to improve the unwrapping accuracy; through a multi-scale unwrapping strategy, optimization solutions are performed at different spatial scales to reduce the risk of unwrapping error propagation, and ultimately obtain a more reliable absolute phase value.

[0027] The traditional Kalman filter assumes that both the system state and the observation equation satisfy a linear relationship, but the actual geological disaster process often has nonlinear characteristics, such as the movement speed of the landslide body and the nonlinear relationship between the sliding distance, and the traditional Kalman filter is difficult to accurately describe this nonlinear dynamic. This application introduces nonlinear state equations and observation equations, and uses the extended Kalman filter algorithm to dynamically estimate the state parameters of the geological disaster area, which can effectively characterize the nonlinear dynamic characteristics of the geological disaster process. The extended Kalman filter predicts and updates the state parameters at each moment by locally linearizing the nonlinear model, and realizes the recursive calculation of state estimation. Compared with the traditional Kalman filter, the extended Kalman filter can handle nonlinear systems and improve the accuracy of state estimation. At the same time, by setting reasonable state parameters, such as surface deformation rate, acceleration and deformation direction, the dynamic evolution process of geological disasters can be fully characterized, providing a more accurate and reliable basis for disaster warning. BRIEF DESCRIPTION OF THE DRAWINGS

[0028] The present application will be further described in the form of exemplary embodiments, which will be described in detail by the accompanying drawings. These embodiments are not restrictive, and in these embodiments, the same number represents the same structure, wherein:

[0029] Figure 1 is an exemplary flow chart of an intelligent data processing method for geological disaster early warning according to some embodiments of the present application;

[0030] Figure 2 is an exemplary flow chart of calculating the absolute phase value according to some embodiments of the present application;

[0031] Figure 3 This is an exemplary flow chart for calculating the state parameter estimation of a geological disaster area according to some embodiments of the present application. DETAILED DESCRIPTION

[0032] The method and system provided in the embodiments of the present application are described in detail below with reference to the accompanying drawings.

[0033] like Figure 1 As shown, multi-source heterogeneous data are collected, and the multi-source heterogeneous data include synthetic aperture radar data, optical satellite image data and laser radar data; the terrain characteristics of the geological disaster area are extracted according to the optical satellite image data and the laser radar data, and the terrain characteristics include slope, slope direction and elevation; the synthetic aperture radar data are differentially interfered to obtain differential interference phase data of the geological disaster area; according to the differential interference phase data, the improved minimum cost flow algorithm is used to perform phase unwrapping to obtain continuous absolute phase values The improved minimum cost flow algorithm introduces space-time continuity constraints and considers the correlation between differential interference phase data of different periods; the absolute phase value It is converted into a deformation variable in the line of sight direction, and the surface deformation variable of the geological disaster area is estimated based on the geometric parameters of the synthetic aperture radar data; the surface deformation variable estimation results are nonlinearly estimated to obtain the state parameter estimation of the geological disaster area; the state parameter estimation is used as input and support vector machine (SVM) is used for early warning detection.

[0034] S1, collect multi-source heterogeneous data, select appropriate SAR satellite platforms, such as TerraSAR-X, Sentinel-1, etc., and obtain multi-temporal SAR data according to the geographical location and time span of the study area. The collected SAR data should cover the study area and meet certain temporal resolution and spatial resolution requirements, such as a time interval of no more than 1 month and a spatial resolution better than 10 meters. Obtain the orbital parameters and metadata of the SAR data for subsequent data processing and analysis.

[0035] Select appropriate optical satellite platforms, such as Landsat and Sentinel-2, to obtain multi-temporal optical image data according to the geographical location and time span of the study area. The collected optical image data should cover the study area and meet certain temporal resolution and spatial resolution requirements, such as a time interval of no more than 1 month and a spatial resolution better than 10 meters. Obtain metadata for the optical image data, including imaging time, sensor type, solar altitude angle, etc., for subsequent data processing and analysis.

[0036] According to the geographical location and scope of the study area, select a suitable LiDAR platform, such as airborne LiDAR or ground-based LiDAR system. Develop a LiDAR data acquisition plan, including parameters such as flight altitude, speed, and point cloud density, to ensure that the acquired point cloud data has sufficient density and accuracy, such as a point cloud density of no less than 10 points per square meter. Perform LiDAR data acquisition and record auxiliary information such as GPS time and inertial navigation data for subsequent point cloud data processing.

[0037] S2, based on the optical satellite image data and lidar data, extract the terrain features of the geological disaster area, perform radiometric calibration on the optical satellite image data, and convert the DN value into a radiometric brightness value or a reflectivity value. Perform atmospheric correction to remove the effects of atmospheric scattering and absorption on the image, and obtain a surface reflectivity image. Perform geometric correction to correct the geometric distortion of the image so that the image corresponds to the true position on the ground. Methods such as RPC orthorectification and GCP geometric correction can be used. According to the geographic coordinate range of the geological disaster area, such as the coordinates of the four corner points of a rectangular area, the corrected optical image is cropped to extract the image data covering the target area. The cropped optical image data is saved in a commonly used data format, such as GeoTIFF format, and the metadata information of the image is recorded.

[0038] Obtain high-precision DEM data of the target area, such as SRTM DEM or DEM obtained by airborne LiDAR. Use DEM data and auxiliary data of the image, such as RPC coefficients or GCP, to orthorectify the optical image of the target area. During the orthorectification process, consider the resolution and accuracy of the DEM and select an appropriate resampling method, such as bilinear interpolation or cubic convolution interpolation. Generate an orthorectified image and perform a quality check to ensure that the geometric position of the image is consistent with the DEM.

[0039] De-noising the original LiDAR point cloud data, removing abnormal height value points and low point density areas, and improving the quality of point cloud data. Perform filtering, such as using statistical filtering or morphological filtering methods, to remove discrete noise points in the point cloud. Extract the LiDAR point cloud data covering the target area based on the geographic coordinate range of the geological disaster area.

[0040] The orthorectified optical image and the preprocessed LiDAR point cloud data are registered to ensure that they are spatially aligned. A high-resolution DSM is generated using LiDAR point cloud data. Irregular point cloud data can be interpolated into regular grid data using methods such as TIN triangulation and Kruger. The generated DSM is evaluated for quality, such as checking the DSM's elevation accuracy, spatial resolution and other indicators to ensure that it meets the requirements for terrain feature extraction.

[0041] Calculate the slope of DSM data to get the slope value of each pixel. You can use the maximum descending slope method, polynomial fitting and other methods. Calculate the aspect of DSM data to get the aspect value of each pixel. You can use the inverse tangent function method, polynomial fitting and other methods. Interpolate and smooth the DSM data to reduce the impact of data noise and discrete points and obtain smooth elevation data. Store the calculated slope, aspect and elevation data in a certain data format to form a terrain feature data set for subsequent analysis and application.

[0042] S3, differential interferometry processing is performed on synthetic aperture radar data to obtain differential interferometry phase data of geological disaster areas. First, precise orbit data of synthetic aperture radar data, such as ephemeris data or GPS positioning data, is obtained. Using orbit data, the position and velocity of the satellite platform at the time of imaging are calculated to establish the orbit model of the satellite. According to the orbit model, orbit correction is performed on radar data to eliminate orbit errors caused by the movement of the satellite platform, such as perturbation and clock error. Through orbit correction, the spatial consistency of radar data acquired at different times is ensured, laying the foundation for subsequent interferometry processing.

[0043] Obtain high-precision DEM data of the study area, such as SRTM DEM or DEM obtained by airborne LiDAR. Use DEM data to perform terrain correction on synthetic aperture radar data to remove the impact of terrain undulations on radar signals. The main steps of terrain correction include: calculating the incident angle and azimuth of the radar signal; calculating the slant range and ground distance of each pixel based on the incident angle and azimuth; using DEM data, calculating the terrain phase of each pixel; removing the terrain phase from the radar signal to obtain radar data without terrain effects.

[0044] Select radar data at one time as the master image, and data at other times as slave images. Use image registration algorithms, such as cross-correlation algorithms or feature point matching algorithms, to register slave images to the geometric coordinate system of the master image. During the registration process, consider the characteristics of radar data, such as oblique imaging, azimuth compression, etc., and select appropriate registration parameters and strategies. Through data registration, ensure that radar data acquired at different times are accurately aligned in space, and provide accurate pixel correspondence for subsequent interference processing.

[0045] According to the time range and frequency of geological disasters, such as the occurrence time and evolution period of landslides, select a suitable time baseline. The selection of the time baseline should take into account the detection capability and time resolution of the deformation signal, and generally select a time interval close to or slightly smaller than the deformation period. According to the orbital parameters and imaging mode of the radar data, such as orbital altitude, incident angle, etc., select a suitable spatial baseline. The selection of the spatial baseline should take into account coherence and phase sensitivity, and generally select a spatial baseline combination that meets a certain coherence threshold and deformation detection accuracy.

[0046] The master and slave images of each data pair are conjugate multiplied to obtain a complex interferogram, which contains phase information caused by factors such as surface deformation and topography. The complex interferogram is processed using multi-view processing techniques, such as slant range averaging or multi-view filtering, to improve phase quality and reduce speckle noise and phase discontinuity. Combined with orbital parameters and DEM data, the differential interferometric phase caused by surface deformation is separated from the complex interferogram, and interference factors such as terrain residual phase and atmospheric effects are removed.

[0047] Adaptive filtering algorithms, such as Goldstein filtering or non-local mean filtering, are used to filter the differential interferometry phase. Adaptive filtering can adaptively adjust the filtering parameters according to the local statistical characteristics of the phase, such as variance, gradient, etc., and effectively suppress phase noise while retaining deformation information. Through phase filtering, the signal-to-noise ratio and spatial continuity of the differential interferometry phase are improved, providing high-quality input data for subsequent phase unwrapping.

[0048] Phase unwrapping algorithms such as the least squares method and the branch-cut method are used to convert the differential interferometric phase from the wrapped phase modulo 2π to a continuous phase value. During the phase unwrapping process, factors such as phase gradient, residual, and coherence are considered to select appropriate unwrapping paths and strategies to avoid unwrapping errors as much as possible. Through phase unwrapping, continuous phase values ​​are obtained to represent the deformation of the surface along the radar line of sight, providing basic data for subsequent deformation inversion and analysis.

[0049] like Figure 2 As shown in S4, based on the differential interference phase data, the improved minimum cost flow algorithm is used to perform phase unwrapping to obtain continuous absolute phase values , including: S41, converting the differential interferometer phase data into a spatiotemporal network flow graph: representing the differential interferometer phase data as a three-dimensional matrix, the dimension of the matrix is ​​(number of time steps, number of rows, number of columns). Each pixel in the space-time network flow graph G is represented as a node, and the coordinates of the node are (t, m, n), corresponding to the pixel in the matrix For example, the node Represents the pixel at the t-th time step, m-th row, and n-th column.

[0050] Calculate the phase difference between adjacent pixels in space, that is, calculate the phase difference between each pixel and its four adjacent pixels above, below, left and right in the same time step. , calculate it and the node Phase difference between: Phase difference between adjacent pixels above: ; Phase difference between adjacent pixels below: ; Phase difference between adjacent pixels on the left: ; Phase difference between adjacent pixels on the right: .

[0051] Calculate the phase change between adjacent pixels in time, that is, calculate the phase change between the two previous and next time steps for each pixel at the same position. , calculate it and the node , Phase change between: Phase change from the previous time step: ; Phase change in the next time step: .

[0052] In the spatiotemporal network flow graph G, the phase difference between adjacent pixels in space and the phase change between adjacent pixels in time are taken as the cost of the edge: With Node The edge cost between is: ;node With Node The edge cost between is: ;node With Node The edge cost between is: ;node With Node The edge cost between is: ;node With Node The edge cost between is: ;node With Node The edge cost between them is: Through the above steps, the differential interference phase data is converted into a spatiotemporal network flow graph G, where nodes represent pixels, edges represent phase differences or phase changes between adjacent pixels, and the cost of edges is the absolute value of the phase difference or phase change.

[0053] S42, set the source node, sink node and flow rate. In the spatiotemporal network flow graph G, set the first pixel point of the first time step (i.e., the node with coordinates (0, 0, 0)) as the source node s, that is, . Set the last pixel point of the last time step (i.e., the node with coordinates (T-1, M-1, N-1)) as the sink node t, that is, . Count the total number of pixels N of the differential interferometric phase data, that is, N = T × M × N. For example, if the dimension of the differential interferometric phase data is (10, 100, 100), the total number of pixels N = 10 × 100 × 100 = 100,000. Let N be the flow of the source node s, indicating that the flow out of the source node s is N. In the spatiotemporal network flow graph G, add a directed edge from the source node s to the super source node ss, the edge capacity is N, and the edge cost is 0, that is, f (s, ss) = N, c (s, ss) = 0.

[0054] Let -N be the flow of sink node t, indicating that the flow into sink node t is N. In the spatiotemporal network flow graph G, add a directed edge from super sink node tt to sink node t, with the capacity of the edge being N and the cost of the edge being 0, i.e., f(tt, t) = N, c(tt, t) = 0. For other nodes except source node s and sink node t , setting their flow to 0 means that the net flow of these nodes is 0, that is, the inflow and outflow are equal. In the spatiotemporal network flow graph G, for each node (In addition to s and t), add a slave node The directed edge to the super sink node tt has an edge capacity of ∞ and an edge cost of 0, that is, , . In the spatiotemporal network flow graph G, source nodes, sink nodes and flow are set to transform the phase unwrapping problem of differential interferometry phase into a minimum cost flow problem. Among them, the source node s represents the starting point of phase unwrapping, the sink node t represents the end point of phase unwrapping, the flow of the source node and the sink node is the total number of pixels N, and the net flow of other nodes is 0 to ensure the flow balance constraint.

[0055] S43, perform multi-scale decomposition on the differential interferometric phase data, and use wavelet transform to perform multi-scale decomposition on the differential interferometric phase data. Commonly used wavelet basis functions include Haar wavelet, Daubechies wavelet, etc. Taking Haar wavelet as an example, perform two-dimensional wavelet decomposition on the matrix Φ: perform one-dimensional wavelet decomposition on each row of the matrix Φ to obtain low-frequency subband L and high-frequency subband H. Perform one-dimensional wavelet decomposition on each column of L and H to obtain four subbands: LL, LH, HL, and HH. LL represents the low-frequency subband, which retains the low-frequency information of the original data; LH, HL, and HH represent high-frequency subbands, which retain the high-frequency information in the horizontal, vertical, and diagonal directions, respectively. Downsample the low-frequency subband LL to obtain the phase map of the next scale. The downsampling operation usually uses row and column interval extraction to reduce the number of rows and columns of the LL subband to 1 / 2 of the original. For example, if the dimension of LL is (M / 2, N / 2), the downsampled phase map The dimension of is (M / 4, N / 4). Repeat steps 2 and 3 to get the downsampled phase map Perform wavelet decomposition and downsampling to obtain a lower resolution phase map Φ_2. And so on, until the preset coarsest scale L is reached. Phase diagram at the level scale The dimension is .

[0056] The coarsest scale L is usually determined based on the spatial resolution of the differential interferometric phase data and the size of the target area, and is generally an integer between 3 and 5. Multi-scale decomposition obtains a series of phase maps with different spatial resolutions. ,in, Represents the original differential interferometry phase data, with the highest resolution; Represents the phase map at the coarsest scale, with the lowest resolution. Phase map The original phase data is retained in The low-frequency information at the level scale reflects the overall change trend of the phase.

[0057] Phase images at different scales capture the multi-scale features of phase information step by step, from coarse to fine. Wavelet transform is used to perform multi-scale decomposition of differential interferometric phase data to obtain a series of phase images with different spatial resolutions, providing input data for subsequent multi-scale phase unwrapping. Multi-scale decomposition can effectively extract the multi-scale features of phase information, reduce the impact of phase noise, and improve the accuracy and robustness of phase unwrapping.

[0058] S44, obtain the initial disentanglement result at the coarsest scale, and construct the spatiotemporal network flow graph at the coarsest scale L obtained by multi-scale decomposition .picture The nodes in the graph represent the phase diagram at the coarsest scale. The edge represents the phase difference or phase change between adjacent pixels. In step S42, the flow rates of the source node s, the sink node t, and the pixel node are set.

[0059] Define the objective function of the minimum cost flow problem, that is, to find the minimum cost flow from the source node s to the sink node t: Maximum flow constraint: The flow outflow from the source node s is equal to the flow inflow to the sink node t, which is equal to the total number of pixels N. Flow balance constraint: For nodes other than the source node s and the sink node t, the inflow and outflow are equal. Capacity constraint: The flow on each edge does not exceed the capacity of the edge. Minimum cost: Under the premise of satisfying the above constraints, the total edge cost is minimized. Use the minimum cost maximum flow algorithm to solve the above minimum cost flow problem and obtain the minimum cost flow f_L from the source node s to the sink node t. Commonly used minimum cost maximum flow algorithms include the continuous shortest path algorithm, the primal dual algorithm, etc.

[0060] According to the minimum cost flow , extract the minimum cost path . Minimum cost path is a path from source node s to sink node t, where the total cost of the edges passing through it is the minimum. The nodes on the grid are connected in chronological order to obtain the initial disentanglement result at the coarsest scale L. .path Each node v on the graph corresponds to the phase diagram at the coarsest scale The edge between adjacent nodes represents the phase difference or phase change. The phase difference between adjacent pixels can be determined based on the direction and weight of the edge.

[0061] Starting from the source node s, along the path The phase differences are accumulated in sequence to obtain the absolute phase value of each pixel. The initial unwrapping result is obtained at the coarsest scale L. , as the initial value for subsequent fine-scale unwrapping. Since the phase image at the coarsest scale has the lowest resolution, smaller phase noise and ambiguity, the initial unwrapping result is usually highly reliable and can provide a good starting point for subsequent iterative optimization.

[0062] S45, iteratively optimize the disentanglement results at each scale, starting from the coarsest scale L, and iteratively optimize the disentanglement results step by step towards the fine scale. Assume that the current processing is the lth scale, and the previous scale has been obtained The untangling result . The unwrapping result of the previous level Interpolate to the current scale l to get the initial unwrapping result .

[0063] Use bilinear interpolation or spline interpolation to The size is enlarged to the phase diagram at the current scale l During the interpolation process, The phase value of the smooth transition is used to reduce the phase difference between adjacent pixels. According to the terrain characteristics at the current scale l, the spatiotemporal network flow graph is The side costs in the calculation are adjusted.

[0064] The local slope value at each pixel can be calculated using digital elevation model (DEM) data or other terrain data. According to the preset slope threshold, the area where the pixel is located is divided into a complex terrain area (the first area) and a flat terrain area (the second area).

[0065] For pixels in complex terrain areas, the proportion of the phase difference between the pixels and adjacent pixels in the edge cost is increased, so that phase unwrapping is more dependent on local phase gradient information. The edge cost can be adjusted by introducing a terrain complexity factor α: ,in, represents the cost of the edge between pixels m and n, Represents the phase difference between two pixels. Represents the original edge cost. The value range of α is [0, 1]. The larger the α is, the higher the proportion of phase difference in the edge cost is.

[0066] For pixels in areas with gentle terrain, the proportion of phase difference in edge cost is reduced to make phase unwrapping smoother and reduce the impact of local phase noise. The edge cost can be adjusted by introducing the terrain gentleness factor β: , where the value range of β is [0, 1]. The smaller β is, the lower the proportion of phase difference in edge cost is.

[0067] In space-time network flow graph In the above example, the phase change between each node and its adjacent nodes in the time series is calculated. , calculate its time difference with the previous node and the next moment node The phase change between: ,in, Representation Node The phase value at and Respectively represent the phase value of the node at the same position at the previous moment and the next moment.

[0068] The phase change Multiply by the time continuity weight , get the time continuity cost . Time continuity weight According to the time baseline length T, the weight calculation formula is: α controls the weight amplitude, and β controls the rate at which the weight increases over the time baseline. The values ​​of α and β can be adjusted according to actual conditions. Usually, α takes a value between 0.11 and β takes a value between 0.010.1. Time continuity cost It reflects the continuity of the phase in the time series. The larger it is, the more drastic the phase change over time and the worse the time continuity.

[0069] Time Continuity Cost Superimposed on the space-time network flow graph The cost of the corresponding side The adjusted side cost is For the node Its temporal neighbors and The cost of updating the edge between is: , adjusted edge cost The spatial phase difference and temporal phase change are comprehensively considered, making the phase unwrapping result more continuous in time.

[0070] After the adjustment, the spatiotemporal network flow graph In the example, the minimum cost maximum flow algorithm is used to obtain the minimum cost path. ,Will The nodes on the left and right are connected in sequence to obtain the disentanglement result at the current scale l. Repeat until the finest scale 0 is reached, and the final unwrapping result is obtained. The unwrapping results are iteratively optimized at each scale, gradually improving the accuracy and reliability of phase unwrapping. By introducing terrain features and time continuity constraints, phase noise and ambiguity can be effectively suppressed, and more accurate and continuous unwrapping results can be obtained.

[0071] S46, fusion of disentanglement results at different scales. In the process of multi-scale disentanglement, a series of disentanglement results at different scales are obtained. ,in, Represents the unwrapping result at the lth scale. The unwrapping results at different scales are fused using the weighted average method to obtain the final absolute phase value φ. For each pixel point (m, n), calculate its absolute phase value at different scales .

[0072] Calculate the weight of the disentanglement result at each scale , the calculation of weights can take into account the following factors: the higher the scale level, the lower the resolution of the unwrapping result, the smaller the phase noise and ambiguity, and the larger the weight should be. The weight can be set proportional to the scale level, for example The higher the terrain complexity, the lower the reliability of the disentanglement result, and the smaller the weight should be. The terrain complexity factor at each scale can be used to calculate the To adjust the weights, for example The better the time continuity, the more reliable the disentanglement result, and the greater the weight should be. To adjust the weights, for example .

[0073] The disentanglement results at different scales The corresponding weight Multiply and sum to get the final absolute phase value at the pixel point (m, n) : , for all pixels, repeat to get the complete absolute phase value matrix . For the absolute phase value matrix Post-processing, including phase filtering, phase conversion and other operations, is performed to obtain the final unwrapping result.

[0074] Phase filtering can use methods such as Gaussian filtering and median filtering to further remove phase noise and outliers. The results of multi-scale unwrapping are fused to obtain the final absolute phase value. Compared with the unwrapping results at a single scale, the results of multi-scale fusion can better balance local details and global trends, and improve the accuracy and reliability of the unwrapping results.

[0075] like Figure 3As shown, S5, nonlinear estimation is performed on the surface deformation variable estimation result to obtain the state parameter estimation of the geological disaster area, including: S51, taking the surface deformation variable as an observation quantity and taking the state parameter of the geological disaster area as a state quantity: the observation quantity Represents the surface deformation at time k, including the displacement of the landslide block , surface tilt angle and surface crack width Etc., can be obtained through differential interferometry SAR, GPS monitoring and other methods. Represents the state parameters of the geological disaster area at time k, including the displacement of the landslide block ,speed Etc., reflecting the movement state of the landslide body.

[0076] S52, establish the nonlinear state equation and observation equation between the observed quantity and the state quantity; state equation: , f(*) is the nonlinear state transfer function, which represents the state quantity The state quantity at the previous moment , control amount and process noise The non-linear relationship between them. Represents external factors that affect landslide stability, such as rainfall and groundwater level wait. represents the model error and external disturbance in the state equation, which is usually assumed to be Gaussian white noise.

[0077] The specific form of f(*) can be constructed based on the landslide stability analysis model and the hydrological model, for example: ,in, Represent the nonlinear state transfer functions of landslide displacement, velocity and acceleration respectively.

[0078] Observation equation: , h(*) is a nonlinear observation function, which represents the observed quantity With state quantity and observation noise The non-linear relationship between them. represents the measurement error in the observation equation, which is usually assumed to be Gaussian white noise. The specific form of h(*) can be used to relate the state quantity to the observation quantity, for example: ,in, They represent the nonlinear observation functions of landslide displacement, surface inclination angle, and surface crack width respectively.

[0079] S53, initialize the estimated value and covariance matrix of the state parameter; set the estimated value X(0|0) of the state parameter at the initial moment (k=0), which can usually be set based on prior knowledge or experience, such as initializing the displacement, velocity, and acceleration to 0. Set the state parameter covariance matrix P(0|0) at the initial moment, which represents the uncertainty of the estimated value of the state parameter, and can usually be set to a diagonal matrix, with the diagonal elements being the prior variance of the state parameter.

[0080] S54, predict the state parameters at the current moment based on the estimated values ​​of the state parameters at the previous moment and the state equation; according to the state equation, the estimated values ​​of the state parameters at the previous moment (k-1) are and control volume Substitute the state transfer function f(*) to get the predicted value of the state parameter at the current time (k) : , calculate the covariance matrix of the predicted values ​​of the state parameters , which represents the uncertainty of the predicted value: ,in, is the state transfer function f(*) in The Jacobian matrix at , is the covariance matrix of the process noise.

[0081] S55, based on the observed quantity and observation equation at the current moment, using the extended Kalman filter algorithm to update the estimated value of the state parameter at the current moment; based on the observation equation, the predicted value of the state parameter at the current moment Substitute the observation function h(*) to get the predicted value of the observation at the current moment : , calculate the residual between the predicted value of the observation and the actual observation and the covariance matrix of the residuals : , ,in, is the observation function h(*) in The Jacobian matrix at , is the covariance matrix of the observation noise. Calculate the Kalman gain matrix : , update the state parameter estimates and the covariance matrix : , , where I is the identity matrix.

[0082] S56, repeat steps S53 to S55 until the estimated values ​​of the state parameters at all times are updated to obtain the estimated state parameters of the geological disaster area. Set the time step k from 1 to N, where N is the total number of observation moments. For each moment k, repeat steps S53 to S55 to update the estimated state parameters. and the covariance matrix .

[0083] The extended Kalman filter algorithm is used to perform nonlinear estimation on the surface shape variable estimation results to obtain the state parameter estimation of the geological disaster area. This method takes into account the nonlinear relationship between the observed quantity and the state quantity, establishes a dynamic model through the state equation and the observation equation, and uses the extended Kalman filter algorithm to recursively update the state parameter estimation value, which can effectively integrate multi-source observation information and improve the accuracy and reliability of state parameter estimation. At the same time, the uncertainty of the state parameter estimation value can be quantified through the propagation and update of the covariance matrix.

[0084] S6, collect historical geological disaster data, including the estimated state parameter values, environmental factors, disaster level and other information when the geological disaster occurred. The estimated state parameter values ​​of the geological disaster area are used as the input features of the SVM model, for example: the displacement of the landslide block ,speed Estimated values ​​of state parameters. External factors affecting landslide stability, such as rainfall and groundwater level Etc. According to the severity of geological disasters, geological disaster warnings are divided into multiple levels as the output categories of the SVM model, for example: Blue warning: low landslide risk, need to continue monitoring. Yellow warning: high landslide risk, need to strengthen monitoring and preventive measures. Orange warning: high landslide risk, need to activate emergency plans, make evacuation preparations. Red warning: extremely high landslide risk, need to evacuate the danger zone immediately.

[0085] Select the kernel function type and kernel function parameters of the SVM model, evaluate the performance of different kernel function types and parameter combinations through cross-validation and other methods, and select the optimal kernel function setting. For different warning levels, design corresponding warning information content, for example: Blue warning: indicates that the landslide risk is low, and it is recommended to continue monitoring and pay attention to weather changes and human activities. Yellow warning: indicates that the landslide risk is high, and it is recommended to strengthen monitoring, formulate emergency plans, and prepare necessary protective measures. Orange warning: indicates that the landslide risk is very high, it is recommended to activate the emergency plan, make evacuation preparations, and pay close attention to the changes in the landslide body. Red warning: indicates that the landslide risk is extremely high, it is recommended to evacuate the dangerous area immediately, ensure the safety of personnel, and coordinate with relevant departments to carry out rescue work. The SVM model is used to estimate the state parameters of the geological disaster area for early warning detection, determine the warning level according to the detection results, and generate corresponding warning information. This method makes full use of historical geological disaster data, establishes a nonlinear mapping relationship between state parameters and disaster risks through machine learning, can automatically identify the early signs of landslides, and issue early warning information in a timely manner, providing a scientific basis for the prevention and response of geological disasters.

Claims

1. An intelligent data processing method for geological disaster early warning, characterized in that: include: S1, collecting multi-source heterogeneous data, where the multi-source heterogeneous data includes synthetic aperture radar data, optical satellite image data, and lidar data; S2, extracting the terrain features of the geological disaster area based on the optical satellite image data and the laser radar data, wherein the terrain features include slope, aspect and elevation; S3, performing differential interferometry processing on the synthetic aperture radar data to obtain differential interferometry phase data of the geological disaster area; S4, based on the differential interferometry phase data, the improved minimum cost flow algorithm is used to perform phase unwrapping to obtain continuous absolute phase values ; The improved minimum cost flow algorithm introduces space-time continuity constraints and considers the correlation between differential interferometry phase data of different periods; The absolute phase value Convert it into deformation in the line of sight direction, and estimate the surface deformation of the geological disaster area based on the geometric parameters of synthetic aperture radar data; S5, performing nonlinear estimation on the surface shape variable estimation results to obtain the state parameter estimation of the geological disaster area; S6, taking the state parameter estimation as input, uses the support vector machine SVM for early warning detection.

2. The intelligent data processing method for geological disaster early warning according to claim 1 is characterized in that: S4, based on the differential interferometry phase data, the improved minimum cost flow algorithm is used to perform phase unwrapping to obtain continuous absolute phase values ,include: S41, converting the differential interference phase data into a spatiotemporal network flow graph, wherein each node in the spatiotemporal network flow graph represents a pixel point at a spatiotemporal position; the phase difference between adjacent pixel points in space and the phase change between adjacent pixel points in time are represented as edge costs; S42, setting a source node and a sink node in the spatiotemporal network flow graph, and setting a flow rate equal to the total number of pixels of the differential interferometric phase data; S43, performing multi-scale decomposition on the differential interferometry phase data to obtain phase maps at different spatial scales; S44, at the coarsest scale, obtain the minimum cost path in the space-time network flow graph through the minimum cost maximum flow algorithm as the initial disentanglement result; S45, at each scale, using the disentanglement result of the previous scale as the initial value, and adjusting the cost of the edge according to the terrain characteristics and the space-time continuity constraint, and obtaining the disentanglement result of the current scale through the minimum cost maximum flow algorithm; S46, fuse the unwrapping results of each scale to obtain the final absolute phase value .

3. The intelligent data processing method for geological disaster early warning according to claim 2 is characterized in that: S42, setting a source node and a sink node in the spatiotemporal network flow graph, and setting a flow rate equal to the total number of pixels of the differential interferometric phase data, including: In the spatiotemporal network flow graph, the first pixel point of the first time step is set as the source node, and the last pixel point of the last time step is set as the sink node; Count the total number of pixels N of the differential interferometric phase data, take N as the flow of the source node, and take -N as the flow of the sink node; Set the flow rates of all nodes except source and sink nodes to zero.

4. The intelligent data processing method for geological disaster early warning according to claim 2 is characterized in that: S45, at each scale, the disentanglement result of the previous scale is used as the initial value, and the cost of the edge is adjusted according to the terrain characteristics and the space-time continuity constraint, and the disentanglement result of the current scale is obtained through the minimum cost maximum flow algorithm, including: Interpolate the unwrapping result of the previous scale to the current scale as the initial value of the unwrapping of the current scale; According to the terrain features, the local slope value at each node is calculated, and the area where the node is located is divided into a first area and a second area according to a preset threshold, wherein the first area represents an area with complex terrain, and the second area represents an area with flat terrain; For nodes in the first area, the proportion of the phase difference on the edge between the node and the adjacent node in the edge cost is increased, and for nodes in the second area, the proportion is reduced; Calculate the phase change between each node and its adjacent nodes in the time series in the spatiotemporal network flow graph constrained by terrain features; Multiplying the phase change by a time continuity weight to obtain a time continuity fee, wherein the time continuity weight is set according to the length of the time baseline; Add the time continuity fee to the fee of the corresponding edge to get the adjusted edge fee; The minimum cost maximum flow algorithm is used to obtain the minimum cost path in the space-time network flow graph after adjusting the edge costs as the disentanglement result of the current scale.

5. The intelligent data processing method for geological disaster early warning according to claim 4 is characterized in that: Set the time continuity weight according to the following formula : ； Among them, T represents the length of the time baseline, α controls the overall amplitude of the time continuity weight, and β controls the rate at which the time continuity weight increases with the increase of the time baseline length; T represents the length of the time baseline.

6. The intelligent data processing method for geological disaster early warning according to claim 1 is characterized in that: The absolute phase value The deformation variables in the line of sight are converted into deformation variables, and the surface deformation variables of the geological disaster area are estimated based on the geometric parameters of the synthetic aperture radar data, including: According to the wavelength λ of the synthetic aperture radar data, the absolute phase value Converted to the deformation in the direction of sight , the conversion formula is: ; Calculate the incident angle θ and azimuth angle α of the synthetic aperture radar signal; According to the incident angle θ and the azimuth angle α, the deformation in the line of sight is Convert to three-dimensional surface shape ; Combined with the topographic characteristics of the geological disaster area, the least squares equations are constructed to solve the three-dimensional surface deformation variables. ; The Kriging interpolation algorithm is used to perform spatial interpolation on the three-dimensional surface deformation variables to obtain the three-dimensional surface deformation field in the geological disaster area.

7. The intelligent data processing method for geological disaster early warning according to claim 6 is characterized in that: Calculate the incident angle θ and azimuth angle α of the synthetic aperture radar signal, including: Extracting orbital parameters of a satellite platform from synthetic aperture radar data, wherein the orbital parameters include satellite position, velocity and attitude; Calculate the geometric relationship between the satellite platform and the geological disaster area based on the orbital parameters and the geographical coordinates of the geological disaster area; According to the geometric relationship, the incident angle θ and azimuth angle α of the synthetic aperture radar signal in the geological disaster area are calculated.

8. The intelligent data processing method for geological disaster early warning according to claim 6 is characterized in that: The deformation in the sight direction is calculated by the following formula Convert to three-dimensional surface shape : ； in, Represents the deformation in the sight direction; Represents the x-axis surface shape variable; Represents the surface shape variable on the y-axis; Represents the surface deformation variable along the z axis; Indicates the azimuth of the radar signal; represents the incident angle of the radar signal; It represents the angle between the ground deformation direction and the incident plane of the radar signal.

9. The intelligent data processing method for geological disaster early warning according to any one of claims 1 to 8, characterized in that: S5, performing nonlinear estimation on the surface shape variable estimation results to obtain the state parameter estimation of the geological disaster area, including: S51, taking the surface deformation variable as an observation quantity and the state parameter of the geological disaster area as a state quantity, wherein the state parameter includes the surface deformation rate and the surface deformation acceleration; S52, establishing a nonlinear state equation and an observation equation between the observed quantity and the state quantity; S53, initializing the estimated values ​​and covariance matrix of the state parameters; S54, predicting the state parameters at the current moment according to the estimated values ​​of the state parameters at the previous moment and the state equation; S55, based on the observed quantity and observation equation at the current moment, using the extended Kalman filter algorithm to update the estimated value of the state parameter at the current moment; S56, repeating steps S53 to S55 until the estimated values ​​of the state parameters at all times are updated to obtain the estimated state parameters of the geological disaster area.

10. An intelligent data processing system for geological disaster early warning, characterized in that: include: At least one processing unit; used to execute instructions to implement the intelligent data processing method for geological disaster warning as described in any one of claims 1 to 9.

Images