Loess-based urban geological disaster assessment method and system

Abstract

The invention relates to the technical field of geological disaster assessment, and particularly discloses a loess-based urban geological disaster assessment method and system, and the method comprises the steps: collecting synthetic aperture radar images to construct a time sequence deformation data set; carrying out adaptive complete set empirical mode decomposition on the data set, and separating high-frequency and low-frequency components; performing Hilbert-Huang transform on the high-frequency component, identifying an irrigation time window and quantifying response intensity; eliminating irrigation period data in the low-frequency component, and interpolating and reconstructing a virtual tectonic settlement curve to obtain a tectonic settlement contribution value; taking the response intensity as an irrigation settlement contribution value, normalizing and weighting the response intensity and the irrigation settlement contribution value to obtain a comprehensive risk value, and outputting a zoning map; according to the method, the natural structure and the human irrigation are decoupled to the settlement contribution, and the evaluation accuracy is improved.

Description

Technical Field

[0001] This invention relates to the field of geological hazard assessment technology, specifically to a method and system for assessing geological hazards in loess cities. Background Technology

[0002] Loess cities refer to cities located in loess-covered areas. Cities like Lanzhou, Yan'an, Xi'an, and Xining in my country are typical examples. Loess possesses unique engineering geological characteristics such as well-developed vertical joints, high water sensitivity, and significant collapsibility, making loess urban areas highly susceptible to geological disasters. The main types of disasters include ground subsidence, landslides, collapses, and loess subsidence. Among these, the proportion of disasters caused by human engineering activities is increasing year by year. In particular, fluctuations in groundwater levels caused by human water use activities such as urban greening irrigation and seepage from pumping projects can easily induce deep loess subsidence, seriously threatening the safety of urban infrastructure and the planned use of land resources.

[0003] Existing geological hazard assessment methods, when processing InSAR time-series deformation data of loess cities, lack the ability to couple and analyze the dynamic indicator of human irrigation activity intensity. They rely solely on a single deformation rate threshold for hazard assessment. This results in an inability to effectively distinguish between natural creep caused by regional tectonic background and deep loess subsidence caused by groundwater level fluctuations due to concentrated summer irrigation when dealing with Yellow River terrace-type cities such as Lanzhou. This ambiguity in causal analysis significantly reduces the accuracy of hazard prediction models and easily leads to the misclassification of areas that can be managed by optimizing irrigation regulation as high-risk engineering management areas, causing unnecessary engineering relocation costs and resource allocation imbalances. Summary of the Invention

[0004] The purpose of this invention is to provide a method and system for assessing geological hazards in loess cities, in order to solve the problems mentioned above.

[0005] The objective of this invention can be achieved through the following technical solutions: A method for assessing geological hazards in loess-covered cities includes the following steps: S1: Collect surface deformation information of the target loess urban area in a continuous time series, construct the settlement-time curve of each monitoring point, and obtain the time series deformation dataset; S2: Perform adaptive full set empirical mode decomposition on the time-series deformation dataset, decompose it into several intrinsic mode function components representing different frequency characteristics, and extract the high-frequency components that characterize seasonal fluctuations and the low-frequency components that characterize long-term trends. S3: Perform Hilbert-Huang transform on the high-frequency components to calculate their instantaneous frequency and instantaneous amplitude. Identify the time window of irrigation activity based on the abrupt change point of the instantaneous frequency, and quantify the response intensity of irrigation events to loess subsidence based on the change in instantaneous amplitude within the corresponding time window. S4: The low-frequency components are reconstructed into deformation indices that characterize the regional tectonic background. Combined with the response intensity of irrigation events, tectonic-driven subsidence contribution values ​​and irrigation-driven subsidence contribution values ​​are generated respectively, so as to achieve mathematical decoupling of two different types of deformation. S5: Based on the decoupled tectonic-driven and irrigation-driven subsidence contribution values, the geological hazard risk of the target area is predicted, and the corrected risk zoning results are output to distinguish the different impacts of natural tectonic movements and human irrigation activities on land subsidence.

[0006] As a further aspect of the present invention: S1 specifically includes: Acquire a stack of synthetic aperture radar images covering the target area and extract distributed scatterers; The optimal phase is estimated by performing covariance matrix eigendecomposition on the interference phase of the distributed scatterer. Adaptive spatiotemporal filtering is applied to the optimal phase to obtain the temporal deformation phase of each monitoring point; The temporal deformation phase is converted into a settlement-time curve to form a temporal deformation dataset.

[0007] As a further aspect of the present invention: S2 specifically includes: Add auxiliary white noise to the time-series deformation dataset and perform extreme point envelope fitting to obtain multiple signals to be processed; For each signal to be processed, local mean filtering is performed, and fluctuation details at different scales are iteratively extracted to obtain a set of intrinsic mode function components; Calculate the zero-crossing rate of each intrinsic mode function component, and classify the components with a zero-crossing rate higher than the threshold as high-frequency components that characterize seasonal fluctuations, and classify the residual trend term as low-frequency components that characterize long-term trends. The high-frequency components and low-frequency components are reconstructed into high-frequency deformation signals and low-frequency deformation signals, respectively.

[0008] As a further aspect of the present invention: S3 specifically includes: Perform Hilbert transform on the high-frequency components to construct an analytic signal, and extract the instantaneous frequency and instantaneous amplitude at each sampling time from the analytic signal; Perform first-order difference operation on the instantaneous frequency sequence, detect the abrupt change when the difference value exceeds the set threshold, and define the continuous interval before and after the abrupt change as the time window for irrigation activities. Within the time window, the instantaneous amplitude is integrated, and the difference between the integral value and the average amplitude of the baseline period is calculated. The difference is used as the response intensity of the irrigation event to loess collapse. The time window and response intensity are associated and stored to form a record of the correspondence between irrigation events and loess subsidence.

[0009] As a further aspect of the present invention: the step of performing Hilbert transform on the high-frequency components to construct an analytic signal, and extracting the instantaneous frequency and instantaneous amplitude at each sampling moment from the analytic signal, specifically includes: The extreme point mirror extension is performed at both ends of the high-frequency component to suppress the endpoint effect caused by boundary processing. Perform Hilbert transform on the extended high-frequency components to obtain the real and imaginary parts of the analytic signal; Perform arctangent operation on the real and imaginary parts to obtain the instantaneous phase sequence, and perform differential calculation on the instantaneous phase sequence to obtain the instantaneous frequency at each sampling time; The sum of the square roots of the real and imaginary parts is taken to obtain the instantaneous amplitude at each sampling time, and the instantaneous frequency is stored in correspondence with the instantaneous amplitude.

[0010] As a further aspect of the present invention: S4 specifically includes: Based on the time window, the time period data corresponding to irrigation activities are removed from the low-frequency components to obtain the remaining low-frequency background segment; Spline interpolation was used to fit the remaining low-frequency background segment to reconstruct the virtual structural settlement curve without irrigation interference throughout the entire monitoring period. The cumulative amount of the virtual tectonic settlement curve at the end of the monitoring period is determined as the tectonic-driven settlement contribution value; The calculated response intensity is directly determined as the contribution value of irrigation-driven settlement, thus completing the mathematical decoupling of the two types of deformation.

[0011] As a further aspect of the present invention: the reconstruction process of the virtual settlement curve is as follows: The remaining low-frequency background segments are arranged in chronological order, and the time coordinates and sedimentation amounts of each remaining data point are extracted. Using the time coordinate as nodes, a piecewise polynomial is constructed between adjacent nodes using cubic spline functions, requiring the polynomial to have continuous first and second derivatives at the nodes; Based on the piecewise polynomial, point-by-point interpolation calculations are performed according to the complete time series of the original monitoring period to generate the estimated tectonic settlement value for each moment. By connecting the estimated tectonic settlement values ​​at all times in chronological order, a virtual tectonic settlement curve without irrigation interference is formed.

[0012] As a further aspect of the present invention: S5 specifically includes: The tectonic-driven settlement contribution value and the irrigation-driven settlement contribution value were normalized to obtain the dimensionless tectonic influence index and irrigation influence index. Multiply the construction impact index by the first preset weight to obtain the construction risk base, and multiply the irrigation impact index by the second preset weight to obtain the irrigation risk base; The risk baseline and irrigation risk baseline are summed at each monitoring point to obtain the comprehensive risk value for the corresponding point. Based on the comprehensive risk value from high to low, different levels of risk zones are delineated, and an overlay risk zoning map containing structural risk zones and irrigation risk zones is output.

[0013] A geological hazard assessment system based on loess cities includes: The time-series deformation data acquisition module collects surface deformation information of the target loess urban area over a continuous time series, constructs the settlement-time curve for each monitoring point, and obtains the time-series deformation dataset. The signal decomposition and component extraction module performs adaptive full set empirical mode decomposition on the time-series deformation dataset, decomposing it into several intrinsic mode function components representing different frequency characteristics, and extracting high-frequency components that characterize seasonal fluctuations and low-frequency components that characterize long-term trends. The irrigation response identification and quantification module performs Hilbert-Huang transform on high-frequency components to calculate their instantaneous frequency and instantaneous amplitude. It identifies the time window of irrigation activity based on the abrupt change point of instantaneous frequency and quantifies the response intensity of irrigation events to loess subsidence based on the change in instantaneous amplitude within the corresponding time window. The deformation causation decoupling module reconstructs low-frequency components into deformation indices characterizing the regional tectonic background, and combines them with the response intensity of irrigation events to generate tectonic-driven subsidence contribution values ​​and irrigation-driven subsidence contribution values ​​respectively, thereby achieving mathematical decoupling of two different causal deformations. The risk prediction and zoning output module predicts the geological hazard risk of the target area based on the decoupled tectonic-driven and irrigation-driven subsidence contribution values, and outputs the corrected risk zoning results to distinguish the different impacts of natural tectonic movements and human irrigation activities on land subsidence.

[0014] The beneficial effects of this invention are: (1) This invention combines adaptive complete set empirical mode decomposition with Hilbert-Huang transform to automatically extract high-frequency deformation signals caused by irrigation activities from time-series deformation data, and accurately identify the time window of irrigation events and quantify their response intensity. This method solves the technical problem of distinguishing between natural tectonic subsidence and human irrigation subsidence in traditional assessments, improves the accuracy of the analysis of the causes of urban land subsidence in loess areas, and provides a more reliable decision-making basis for the prediction of geological hazard risks.

[0015] (2) This invention achieves mathematical decoupling of the contribution of tectonic movement and irrigation activities to land subsidence by constructing a virtual tectonic settlement curve without irrigation interference, and performs weighted risk assessment and risk zoning based on this. This method can effectively avoid the waste of engineering relocation caused by misclassifying controllable irrigation areas as high-risk areas due to attribution errors, and at the same time provides targeted and differentiated guidance for urban land use planning, irrigation water scheduling and geological disaster prevention, with good social and economic benefits. Attached Figure Description

[0016] The invention will now be further described with reference to the accompanying drawings.

[0017] Figure 1 This is a flowchart of the method of the present invention; Figure 2 This is a system block diagram of the present invention. Detailed Implementation

[0018] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0019] Please see Figure 1 As shown, this invention is a method for assessing geological hazards in loess cities, comprising the following steps: S1: Collect surface deformation information of the target loess urban area in a continuous time series, construct the settlement-time curve of each monitoring point, and obtain the time series deformation dataset; S2: Perform adaptive full set empirical mode decomposition on the time-series deformation dataset, decompose it into several intrinsic mode function components representing different frequency characteristics, and extract the high-frequency components that characterize seasonal fluctuations and the low-frequency components that characterize long-term trends. S3: Perform Hilbert-Huang transform on the high-frequency components to calculate their instantaneous frequency and instantaneous amplitude. Identify the time window of irrigation activity based on the abrupt change point of the instantaneous frequency, and quantify the response intensity of irrigation events to loess subsidence based on the change in instantaneous amplitude within the corresponding time window. S4: The low-frequency components are reconstructed into deformation indices that characterize the regional tectonic background. Combined with the response intensity of irrigation events, tectonic-driven subsidence contribution values ​​and irrigation-driven subsidence contribution values ​​are generated respectively, so as to achieve mathematical decoupling of two different types of deformation. S5: Based on the decoupled tectonic-driven and irrigation-driven subsidence contribution values, the geological hazard risk of the target area is predicted, and the corrected risk zoning results are output to distinguish the different impacts of natural tectonic movements and human irrigation activities on land subsidence.

[0020] In S1, surface deformation information of the target loess urban area is collected over a continuous time series, and settlement-time curves are constructed for each monitoring point to obtain a time-series deformation dataset, which specifically includes: First, a stack of synthetic aperture radar (SAR) images covering the target loess urban area is acquired, and distributed scatterers (DSS) are extracted from this image stack. The specific implementation process is as follows: at least 20 SAR single-view complex images with the same orbit and incident angle covering the study area are selected. One of these images is used as the master image, and the remaining images are registered and resampled into the master image space. A homogeneous pixel identification algorithm based on a combination of amplitude deviation information and coherence coefficient is employed. A 7x7 rectangular estimation window is set for each pixel, and the amplitude statistical characteristics of the pixels within the window are calculated. Pixels with amplitude deviation values ​​less than 0.25 and maintaining high coherence over time (annual average coherence coefficient greater than 0.35) are selected as DSS candidate points, ultimately forming a DSS set.

[0021] Secondly, the covariance matrix eigenvalue decomposition is performed on the extracted interference phase of the distributed scatterers to estimate the optimal phase. Specifically, for each distributed scatterer, a 9x9 pixel block is selected centered on it. Using the complex values ​​of all pixels within the block on all interference pairs, a complex matrix with dimensions equal to the number of interference pairs multiplied by the number of pixels is constructed. This matrix is ​​then multiplied by its conjugate transpose to obtain the sample covariance matrix of the distributed scatterer. Eigenvalue decomposition is performed on this covariance matrix to extract the eigenvector corresponding to the largest eigenvalue. The phase component contained in this eigenvector is the optimal phase estimate for the distributed scatterer. This process effectively suppresses the interference of decoherent noise on the interference phase.

[0022] Next, adaptive spatiotemporal filtering is applied to the estimated optimal phase to obtain the temporal deformation phase of each monitoring point. In the spatial domain, an adaptive window based on distance weighting and coherence weighting is used for filtering. The window size is dynamically adjusted between 5 and 15 pixels according to the effective pixel density around the point to filter out atmospheric phase noise and high-frequency noise. In the temporal domain, a filtering method based on singular spectrum analysis is applied to the spatially filtered phase sequence. The phase sequence is embedded into a trajectory matrix, and singular value decomposition is performed on the matrix. The top principal components with a contribution rate of over 80% are retained for reconstruction, thereby separating the nonlinear deformation phase and obtaining the pure temporal deformation phase of each monitoring point.

[0023] Finally, the aforementioned temporal deformation phases are converted into settlement-time curves, forming a temporal deformation dataset. Specifically, the temporal deformation phases at each monitoring point are spatially unwrapped, and the true absolute phase value is recovered using a minimum cost flow algorithm. Based on the wavelength of the synthetic aperture radar sensor (e.g., 5.6 cm for a C-band sensor), the phase value is converted into deformation along the radar line-of-sight by multiplying the unwrapped phase value by the wavelength and then dividing by 4π. The deformation values ​​corresponding to the acquisition date of each image are arranged chronologically, with the date as the x-axis and the deformation value as the y-axis, constructing a complete settlement-time curve for each monitoring point. The settlement-time curves of all points together constitute a temporal deformation dataset covering the entire region.

[0024] In S2, an adaptive full ensemble empirical mode decomposition is performed on the time-series deformation dataset, decomposing it into several intrinsic mode function components representing different frequency characteristics. From these components, high-frequency components characterizing seasonal fluctuations and low-frequency components characterizing long-term trends are extracted, specifically including: First, auxiliary white noise is added to the time-series deformation dataset, and extreme point envelope fitting is performed to obtain multiple signals to be processed. Specifically, the overall average iteration is set to 100. During each iteration, a Gaussian white noise sequence with an amplitude of 0.2 times the standard deviation of the dataset is added to the original time-series deformation dataset, resulting in a noisy signal. For this noisy signal, all local maxima and local minima are identified. Cubic spline functions are used to interpolate and fit the maxima and minima respectively, forming upper and lower envelopes. The mean of the upper and lower envelopes is calculated to obtain the mean envelope. The noisy signal is subtracted from this mean envelope to obtain a preliminary component. The above process is repeated until the component satisfies two conditions of the intrinsic mode function: first, the absolute value of the difference between the number of extreme points and the number of zero-crossing points in the entire data sequence does not exceed 1; second, at any sampling time, the mean of the upper envelope fitted by the local maxima and the lower envelope fitted by the local minima is 0. At this point, the first intrinsic mode function (EMF) component is obtained. This component is subtracted from the original noisy signal to obtain the residual signal. The above filtering process is repeated on the residual signal until no new EMFs can be decomposed from it. The above complete decomposition process is then performed on noisy signals with 100 different sets of white noise added, resulting in a set of 100 EMFs.

[0025] Secondly, an overall average calculation is performed on the intrinsic mode functions of each order in each set of decomposition results to eliminate the influence of white noise. Specifically, the intrinsic mode functions of corresponding orders in all 100 sets of decomposition results are accumulated point by point, and then divided by 100 to obtain the final 1st, 2nd, up to nth order intrinsic mode function components. These components are arranged in descending order of frequency, that is, the 1st order component has the highest fluctuation frequency, and the fluctuation frequency of the component gradually decreases as the order increases over time.

[0026] Next, the zero-crossing rate of each intrinsic mode function (IMF) component is calculated, and the components are classified into high-frequency and low-frequency components based on the level of the zero-crossing rate. The zero-crossing rate is calculated as follows: for each IMF component, the total number of times the signal value changes from positive to negative or from negative to positive within a one-year sampling period is counted, and this number is divided by the total number of sampling points in a year to obtain the annual zero-crossing rate of the component. A frequency distinction threshold of 0.3 times per year is set, that is, if the annual zero-crossing rate of a component is higher than 0.3, the component is determined to be a high-frequency component representing seasonal fluctuations; if the annual zero-crossing rate is lower than or equal to 0.3, the component is determined to be a low-frequency component representing interannual variation. After the above screening, the components of order 1 to order m with an annual zero-crossing rate higher than 0.3 are classified into the high-frequency group, and the components of order m+1 to n, as well as the residual trend term that cannot be further decomposed, are classified into the low-frequency group.

[0027] Finally, all intrinsic mode function components of the selected high-frequency group are accumulated point by point to reconstruct a new time series, called the high-frequency deformation signal. Simultaneously, all intrinsic mode function components and residual trend terms of the low-frequency group are accumulated point by point to reconstruct another new time series, called the low-frequency deformation signal. This completes the separation of different frequency components from the original time-series deformation data, yielding two independent signals corresponding to seasonal fluctuations and long-term trends, respectively.

[0028] In S3, a Hilbert-Huang transform is applied to the high-frequency components to calculate their instantaneous frequency and amplitude. The time window for irrigation activities is identified based on the abrupt changes in instantaneous frequency, and the response intensity of irrigation events to loess subsidence is quantified based on the variation in instantaneous amplitude within the corresponding time window. Specifically, this includes: First, extreme point mirroring extension is performed on the high-frequency components to suppress end-point effects. Specifically, all local maxima and minima of the high-frequency components are identified. Three extreme points at the beginning (including maxima and minima) and three at the end of the signal are selected. Using the signal's two boundary lines as axes of symmetry, these extreme points are symmetrically mapped onto the outside of the signal, forming an extension sequence 10% the length of the original signal. The extended signal is then used as the new signal to be processed.

[0029] Then, a Hilbert transform is applied to the extended high-frequency components. The Hilbert transform is defined in the frequency domain as shifting the positive frequency portion of the signal by -90 degrees and the negative frequency portion by +90 degrees. In actual calculations, a Fast Fourier Transform is performed on the extended signal, its spectrum is multiplied by the frequency domain response function of the Hilbert transform (multiplied by the negative imaginary unit for positive frequencies and by the positive imaginary unit for negative frequencies), and then subjected to an inverse Fourier transform to obtain the Hilbert transform result. This result is the imaginary part of the analytic signal, with the original extended signal as the real part, thus constructing the analytic signal: ;in, To analyze the signal, For the extended high-frequency components, For its Hilbert transformation The unit is the imaginary unit. Instantaneous amplitude and instantaneous phase are extracted from the analytic signal. Instantaneous amplitude. We obtain it by taking the square root of the sum of the squares of the real and imaginary parts: Instantaneous phase It is obtained by arctangent operation of the real and imaginary parts: When calculating the instantaneous frequency, the instantaneous phase sequence is differentially divided. Let the sampling time interval be... (Unit: days), then the first Instantaneous frequency at each sampling time )for: ; in and These are the instantaneous phase values ​​at the current and previous moments, respectively. Since the difference cannot be calculated for the first sampling point, its instantaneous frequency is taken to be the same as that of the second point. The instantaneous frequency and instantaneous amplitude at each sampling moment are stored accordingly, forming an instantaneous frequency sequence and an instantaneous amplitude sequence.

[0030] Performing a first-order difference operation on the instantaneous frequency sequence yields the frequency change rate sequence: Set mutation detection threshold The threshold is set to twice the standard deviation of the frequency change rate sequence. Iterating through all time points, when... The absolute value is greater than When a candidate mutation point is reached, that moment is marked as a candidate mutation point. Since irrigation activities typically last for a period of time, the intervals covered by consecutive candidate mutation points are merged, that is, the continuous period between the appearance of the first mutation point and the disappearance of the last mutation point is defined as the time window for irrigation activities. If there are multiple separate time windows in a given year, the start and end times of each window are recorded separately.

[0031] For each identified time window, the integral value of the instantaneous amplitude within the window is calculated. The integral is approximated using the rectangular method, with the starting time of the window set as . The termination time is There are N sampling points in the window, and the sampling interval is... The integral value of the instantaneous amplitude within the window for: The baseline period is a continuous period preceding the specified time window that is free from irrigation disturbances. The length of the baseline period is the same as the window length, typically a 60-day period consisting of 15 days before and after the same calendar date in the previous year (excluding other irrigation periods). The average instantaneous amplitude within the baseline period is calculated. Multiply by the window length to obtain the baseline integral value: The difference between the window integral value and the baseline integral value: As the response intensity of this irrigation event to loess subsidence, this value reflects the cumulative effect of surface deformation caused by irrigation in the amplitude domain.

[0032] The start and end times of each identified irrigation activity time window, along with the corresponding response intensity value, are associated to form a record. For each monitoring point within the target area, a corresponding record of the relationship between the irrigation event and the loess subsidence response intensity is generated and stored in the database for subsequent decoupling analysis. If multiple time windows are detected at a point within a year, multiple records are recorded, with each record independently corresponding to one irrigation event.

[0033] In S4, low-frequency components are reconstructed into deformation indices characterizing the regional tectonic background. Combined with the response intensity of irrigation events, tectonic-driven and irrigation-driven subsidence contribution values ​​are generated separately, achieving mathematical decoupling between the two different types of deformation. Specifically, this includes: First, based on the irrigation activity time windows identified in step S3, the time period data corresponding to the irrigation activities are removed from the low-frequency components to obtain the residual low-frequency background segment. Specifically, the start and end times of each irrigation activity time window stored in step S3 are read. In the time series of the low-frequency components, all data points within the time periods covered by these windows are marked as invalid and removed. For time periods not marked as irrigation windows, their corresponding low-frequency component data points are retained. These retained discrete data points together constitute the residual low-frequency background segment. This background segment only contains long-period deformation information during non-irrigation periods, excluding the direct influence of irrigation activities on the low-frequency trend.

[0034] Secondly, cubic spline interpolation was performed on the remaining low-frequency background segment to reconstruct a virtual tectonic settlement curve without irrigation interference throughout the entire monitoring period. The specific construction process of cubic spline interpolation is as follows: First, all retained data points in the remaining low-frequency background segment are arranged in chronological order, and the corresponding time coordinates (in years) and settlement amounts (in millimeters) for each data point are extracted. These points are used as nodes for spline interpolation. Second, a cubic polynomial function is constructed between every two adjacent nodes, requiring that at each node, the function values ​​of the left and right polynomials are equal, the first derivative values ​​are equal, and the second derivative values ​​are equal. Simultaneously, the second derivative at both endpoints is specified as 0 as a natural boundary condition. Third, based on the above continuity and boundary conditions, a system of linear equations is established regarding the second derivative values ​​at each node. The coefficient matrix of the equation system is a tridiagonal matrix. The system of equations is solved using the pursuit method to obtain the second derivative values ​​at each node. The fourth step involves determining the coefficients of the cubic polynomial between every two adjacent nodes based on the node coordinates, function values ​​at the nodes, and the obtained second derivative values ​​at the nodes, thus obtaining the complete piecewise cubic spline function. The fifth step involves calculating the estimated tectonic settlement value for each time point by substituting the corresponding cubic polynomial for the complete time series within the original monitoring period (including missing moments within the irrigation window). Connecting the estimated tectonic settlement values ​​for all moments in chronological order creates a continuous, smooth, and irrigation-free virtual tectonic settlement curve. This curve fully covers the entire monitoring period and represents the tectonic settlement trend that would have occurred if irrigation had not taken place within the irrigation window.

[0035] Next, the cumulative amount of the virtual tectonic settlement curve at the end of the monitoring period is determined as the tectonic-driven settlement contribution value. Specifically, the settlement amount on the virtual tectonic settlement curve corresponding to the start time of the monitoring period is taken as the baseline value, and the settlement amount on the virtual tectonic settlement curve corresponding to the end time of the monitoring period is taken as the final value. The difference between the final value and the baseline value is calculated. This difference represents the total settlement amount caused solely by the regional tectonic background movement during the entire monitoring period, and is called the tectonic-driven settlement contribution value. The unit is the same as the original settlement amount, expressed in millimeters.

[0036] Finally, the response intensity calculated in step S3 is directly determined as the irrigation-driven settlement contribution value. Specifically, for each monitoring point, the response intensity value of that point in each irrigation event stored in step S3 is read. This value is the loess subsidence contribution caused by that irrigation event. If a monitoring point experiences multiple irrigation events within a year, all response intensity values ​​are summed to obtain the total settlement contribution of that point caused by irrigation activities in that year, which is called the irrigation-driven settlement contribution value, also expressed in millimeters. Thus, by constructing a virtual no-irrigation curve to isolate the tectonic contribution from the low-frequency components while retaining the irrigation contribution obtained from the high-frequency analysis, the mathematical decoupling of tectonic movement and irrigation subsidence—two different causes of deformation—is completed.

[0037] In S5, based on the decoupled tectonic-driven and irrigation-driven subsidence contribution values, the geological hazard risk of the target area is predicted, and the corrected risk zoning results are output to distinguish the different impacts of natural tectonic movements and human irrigation activities on land subsidence, specifically including: First, the decoupled tectonic-driven settlement contribution values ​​and irrigation-driven settlement contribution values ​​are normalized to obtain dimensionless tectonic influence indices and irrigation influence indices. The specific normalization method uses extreme value standardization: for all monitoring points within the target area, the minimum (denoted as the tectonic minimum) and maximum (denoted as the tectonic maximum) tectonic-driven settlement contribution values ​​are identified. For any given monitoring point, its tectonic-driven settlement contribution value is subtracted from the tectonic minimum, and then divided by the difference between the tectonic maximum and minimum; the result is the tectonic influence index for that point. Similarly, for all monitoring points, the minimum (denoted as the irrigation minimum) and maximum (denoted as the irrigation maximum) tectonic settlement contribution values ​​are identified. For any given monitoring point, its irrigation-driven settlement contribution value is subtracted from the irrigation minimum, and then divided by the difference between the irrigation maximum and minimum; the result is the irrigation influence index for that point. Both indices range from 0 to 1, where 0 represents the minimum contribution and 1 represents the maximum contribution.

[0038] Secondly, the structural impact index is multiplied by the first preset weight to obtain the structural risk base; the irrigation impact index is multiplied by the second preset weight to obtain the irrigation risk base. The first and second preset weights are determined as follows: Geological disaster records for the past 10 years in the city where the target area is located are collected, and the number of geological disasters caused by tectonic movement and the number of geological disasters caused by irrigation are counted. The proportions of each type of disaster to the total number of disasters are calculated, and these proportions are used as initial weights. The initial weights are normalized so that their sum is 1, resulting in the final first and second preset weights. For example, if the proportion of disasters caused by tectonic movement is 35% and the proportion of disasters caused by irrigation is 65%, then the first preset weight is 0.35 and the second preset weight is 0.65. The structural impact index of each location is multiplied by the first preset weight to obtain the structural risk base for that location; the irrigation impact index of each location is multiplied by the second preset weight to obtain the irrigation risk base for that location.

[0039] Next, the structural risk baseline and irrigation risk baseline are summed at each monitoring point to obtain the comprehensive risk value for that point. Specifically, for each monitoring point, the structural risk baseline and irrigation risk baseline are directly added together; the sum is the comprehensive risk value for that point. The comprehensive risk value ranges from 0 to 1; a higher value indicates a higher overall risk of geological disasters at that point.

[0040] Finally, risk zones of different levels are delineated according to the comprehensive risk value from high to low, and an overlay risk zoning map including structural risk zones and irrigation risk zones is output. The risk level classification adopts the natural discontinuity grading method: the comprehensive risk values ​​of all monitoring points are sorted, the difference between adjacent values ​​is calculated, and the location with the largest difference is selected as the grading discontinuity, dividing all points into four levels: low risk zone (comprehensive risk value less than 0.25), medium risk zone (comprehensive risk value between 0.25 and 0.50), high risk zone (comprehensive risk value between 0.50 and 0.75), and extremely high risk zone (comprehensive risk value greater than 0.75). Different colors are used to fill the areas of each risk level when outputting the risk zoning map. Meanwhile, the relative magnitudes of the structural risk base and the irrigation risk base are represented by superimposed dot symbols on the map. For example, the structural risk base is represented by a circle and the irrigation risk base is represented by a triangle. The size of the symbol is proportional to the value of the base, thus visually demonstrating the different impacts of natural tectonic movements and human irrigation activities on land subsidence on the same map.

[0041] Please see Figure 2 As shown, a geological hazard assessment system based on loess cities includes: The time-series deformation data acquisition module collects surface deformation information of the target loess urban area over a continuous time series, constructs the settlement-time curve for each monitoring point, and obtains the time-series deformation dataset. The signal decomposition and component extraction module performs adaptive full set empirical mode decomposition on the time-series deformation dataset, decomposing it into several intrinsic mode function components representing different frequency characteristics, and extracting high-frequency components that characterize seasonal fluctuations and low-frequency components that characterize long-term trends. The irrigation response identification and quantification module performs Hilbert-Huang transform on high-frequency components to calculate their instantaneous frequency and instantaneous amplitude. It identifies the time window of irrigation activity based on the abrupt change point of instantaneous frequency and quantifies the response intensity of irrigation events to loess subsidence based on the change in instantaneous amplitude within the corresponding time window. The deformation causation decoupling module reconstructs low-frequency components into deformation indices characterizing the regional tectonic background, and combines them with the response intensity of irrigation events to generate tectonic-driven subsidence contribution values ​​and irrigation-driven subsidence contribution values ​​respectively, thereby achieving mathematical decoupling of two different causal deformations. The risk prediction and zoning output module predicts the geological hazard risk of the target area based on the decoupled tectonic-driven and irrigation-driven subsidence contribution values, and outputs the corrected risk zoning results to distinguish the different impacts of natural tectonic movements and human irrigation activities on land subsidence.

[0042] The working principle of this invention is as follows: Synthetic aperture radar images of the target loess urban area are acquired. Through distributed scatterer extraction, covariance matrix eigenvalue decomposition, and spatiotemporal filtering, settlement-time curves for each monitoring point are constructed to obtain a time-series deformation dataset. Adaptive complete ensemble empirical mode decomposition is performed on the time-series deformation dataset to separate high-frequency components representing seasonal fluctuations and low-frequency components representing long-term trends. Hilbert-Huang transform is applied to the high-frequency components to extract instantaneous frequency and amplitude. Irrigation activity time windows are identified based on instantaneous frequency abrupt change points, and the difference between the integral value of the instantaneous amplitude within the window and the baseline integral value is used as the basis for calculation. The response intensity of irrigation events to loess subsidence is quantified. Based on the time window, the corresponding time period data in the low-frequency component is removed, and cubic spline interpolation is performed on the remaining low-frequency background segment to reconstruct the virtual tectonic settlement curve without irrigation interference. The cumulative amount at the end of the curve is used as the tectonic-driven settlement contribution value, and the response intensity is used as the irrigation-driven settlement contribution value, thus achieving mathematical decoupling of the two types of deformation. The two decoupled contribution values ​​are normalized and weighted to obtain the comprehensive risk value of each point. Based on the comprehensive risk value, risk zones of different levels are delineated, and a superimposed risk zoning map that simultaneously displays the contributions of tectonic risk and irrigation risk is output.

[0043] The foregoing has provided a detailed description of one embodiment of the present invention, but this description is merely a preferred embodiment and should not be construed as limiting the scope of the invention. All equivalent variations and modifications made within the scope of the claims of this invention should still fall within the patent coverage of this invention.

Claims

1. A method for assessing geological hazards in loess-covered cities, characterized in that, Includes the following steps: S1: Collect surface deformation information of the target loess urban area in a continuous time series, construct the settlement-time curve of each monitoring point, and obtain the time series deformation dataset; S2: Perform adaptive full set empirical mode decomposition on the time-series deformation dataset, decompose it into several intrinsic mode function components representing different frequency characteristics, and extract the high-frequency components that characterize seasonal fluctuations and the low-frequency components that characterize long-term trends. S3: Perform Hilbert-Huang transform on the high-frequency components to calculate their instantaneous frequency and instantaneous amplitude. Identify the time window of irrigation activity based on the abrupt change point of the instantaneous frequency, and quantify the response intensity of irrigation events to loess subsidence based on the change in instantaneous amplitude within the corresponding time window. S4: The low-frequency components are reconstructed into deformation indices that characterize the regional tectonic background. Combined with the response intensity of irrigation events, tectonic-driven subsidence contribution values ​​and irrigation-driven subsidence contribution values ​​are generated respectively, so as to achieve mathematical decoupling of two different types of deformation. S5: Based on the decoupled tectonic-driven and irrigation-driven subsidence contribution values, the geological hazard risk of the target area is predicted, and the corrected risk zoning results are output to distinguish the different impacts of natural tectonic movements and human irrigation activities on land subsidence.

2. The method for assessing geological hazards in loess-lined cities according to claim 1, characterized in that, S1 specifically includes: Acquire a stack of synthetic aperture radar images covering the target area and extract distributed scatterers; The optimal phase is estimated by performing covariance matrix eigendecomposition on the interference phase of the distributed scatterer. Adaptive spatiotemporal filtering is applied to the optimal phase to obtain the temporal deformation phase of each monitoring point; The temporal deformation phase is converted into a settlement-time curve to form a temporal deformation dataset.

3. The method for assessing geological hazards in loess-lined cities according to claim 1, characterized in that, S2 specifically includes: Add auxiliary white noise to the time-series deformation dataset and perform extreme point envelope fitting to obtain multiple signals to be processed; For each signal to be processed, local mean filtering is performed, and fluctuation details at different scales are iteratively extracted to obtain a set of intrinsic mode function components; Calculate the zero-crossing rate of each intrinsic mode function component, and classify the components with a zero-crossing rate higher than the threshold as high-frequency components that characterize seasonal fluctuations, and classify the residual trend term as low-frequency components that characterize long-term trends. The high-frequency components and low-frequency components are reconstructed into high-frequency deformation signals and low-frequency deformation signals, respectively.

4. The method for assessing geological hazards in loess-lined cities according to claim 1, characterized in that, S3 specifically includes: Perform Hilbert transform on the high-frequency components to construct an analytic signal, and extract the instantaneous frequency and instantaneous amplitude at each sampling time from the analytic signal; Perform first-order difference operation on the instantaneous frequency sequence, detect the abrupt change when the difference value exceeds the set threshold, and define the continuous interval before and after the abrupt change as the time window for irrigation activities. Within the time window, the instantaneous amplitude is integrated, and the difference between the integral value and the average amplitude of the baseline period is calculated. The difference is used as the response intensity of the irrigation event to loess collapse. The time window and response intensity are associated and stored to form a record of the correspondence between irrigation events and loess subsidence.

5. The method for assessing geological hazards in loess-lined cities according to claim 4, characterized in that, The process of performing a Hilbert transform on the high-frequency components to construct an analytic signal, and extracting the instantaneous frequency and instantaneous amplitude at each sampling moment from the analytic signal, specifically includes: The extreme point mirror extension is performed at both ends of the high-frequency component to suppress the endpoint effect caused by boundary processing. Perform Hilbert transform on the extended high-frequency components to obtain the real and imaginary parts of the analytic signal; Perform arctangent operation on the real and imaginary parts to obtain the instantaneous phase sequence, and perform differential calculation on the instantaneous phase sequence to obtain the instantaneous frequency at each sampling time; The sum of the square roots of the real and imaginary parts is taken to obtain the instantaneous amplitude at each sampling time, and the instantaneous frequency is stored in correspondence with the instantaneous amplitude.

6. The method for assessing geological hazards in loess-lined cities according to claim 1, characterized in that, S4 specifically includes: Based on the time window, the time period data corresponding to irrigation activities are removed from the low-frequency components to obtain the remaining low-frequency background segment; Spline interpolation was used to fit the remaining low-frequency background segment to reconstruct the virtual structural settlement curve without irrigation interference throughout the entire monitoring period. The cumulative amount of the virtual tectonic settlement curve at the end of the monitoring period is determined as the tectonic-driven settlement contribution value; The calculated response intensity is directly determined as the contribution value of irrigation-driven settlement, thus completing the mathematical decoupling of the two types of deformation.

7. The method for assessing geological hazards in loess-lined cities according to claim 6, characterized in that, The reconstruction process of the virtual settlement curve is as follows: The remaining low-frequency background segments are arranged in chronological order, and the time coordinates and sedimentation amounts of each remaining data point are extracted. Using the time coordinate as nodes, a piecewise polynomial is constructed between adjacent nodes using cubic spline functions, requiring the polynomial to have continuous first and second derivatives at the nodes; Based on the piecewise polynomial, point-by-point interpolation calculations are performed according to the complete time series of the original monitoring period to generate the estimated tectonic settlement value for each moment. By connecting the estimated tectonic settlement values ​​at all times in chronological order, a virtual tectonic settlement curve without irrigation interference is formed.

8. The method for assessing geological hazards in loess-lined cities according to claim 1, characterized in that, S5 specifically includes: The tectonic-driven settlement contribution value and the irrigation-driven settlement contribution value were normalized to obtain the dimensionless tectonic influence index and irrigation influence index. Multiply the construction impact index by the first preset weight to obtain the construction risk base, and multiply the irrigation impact index by the second preset weight to obtain the irrigation risk base; The risk baseline and irrigation risk baseline are summed at each monitoring point to obtain the comprehensive risk value for the corresponding point. Based on the comprehensive risk value from high to low, different levels of risk zones are delineated, and an overlay risk zoning map containing structural risk zones and irrigation risk zones is output.

9. A geological hazard assessment system based on loess urban areas, characterized in that, A method for implementing a geological hazard assessment method based on loess urban areas as described in any one of claims 1-8, comprising: The time-series deformation data acquisition module collects surface deformation information of the target loess urban area over a continuous time series, constructs the settlement-time curve for each monitoring point, and obtains the time-series deformation dataset. The signal decomposition and component extraction module performs adaptive full set empirical mode decomposition on the time-series deformation dataset, decomposing it into several intrinsic mode function components representing different frequency characteristics, and extracting high-frequency components that characterize seasonal fluctuations and low-frequency components that characterize long-term trends. The irrigation response identification and quantification module performs Hilbert-Huang transform on high-frequency components to calculate their instantaneous frequency and instantaneous amplitude. It identifies the time window of irrigation activity based on the abrupt change point of instantaneous frequency and quantifies the response intensity of irrigation events to loess subsidence based on the change in instantaneous amplitude within the corresponding time window. The deformation causation decoupling module reconstructs low-frequency components into deformation indices characterizing the regional tectonic background, and combines them with the response intensity of irrigation events to generate tectonic-driven subsidence contribution values ​​and irrigation-driven subsidence contribution values ​​respectively, thereby achieving mathematical decoupling of two different causal deformations. The risk prediction and zoning output module predicts the geological hazard risk of the target area based on the decoupled tectonic-driven and irrigation-driven subsidence contribution values, and outputs the corrected risk zoning results to distinguish the different impacts of natural tectonic movements and human irrigation activities on land subsidence.

Images