A shoreline health assessment method based on spatiotemporal coupling diagnosis

Abstract

This invention relates to the field of riverbank health assessment and discloses a shoreline health assessment method based on spatiotemporal coupling diagnosis. The method includes: acquiring raw data on riverbank status to generate a multi-source dataset of riverbank status; extracting temporal indicators of shoreline characteristics from the multi-source dataset to construct a spatiotemporal characteristic data cube; extracting features from the spatiotemporal characteristic data cube for temporal causal relationship diagnosis and spatial heterogeneity analysis of vegetation shoreline stabilization effects, predicting temporal change trends, delineating risk warning zones, and achieving dynamic risk warning; and generating assessment conclusions based on the dynamic risk warning, covering scoring, diagnosis, warning, and countermeasures. According to the above technical solution, efficient querying and coupled analysis can be supported to identify shoreline sections that are currently in good condition but pose future risks, enabling early intervention and overcoming the limitations of traditional post-evaluation methods. Ultimately, this addresses the industry pain points of low efficiency, narrow coverage, shallow analysis, and weak support in the assessment of natural riverbank conditions.

Description

Technical Field

[0001] This invention relates to the field of river shoreline health assessment, and more specifically, to a shoreline health assessment method based on spatiotemporal coupling diagnosis. Background Technology

[0002] Establishing national river and lake health records and conducting periodic river and lake health assessments are among the core tasks of water conservancy work in the new era. As a key component of river and lake ecosystems, the natural condition of shorelines directly determines the basic health status of rivers and lakes, making their assessment a core element of the river and lake health assessment system. The core objective of river and lake health assessment is not simply to output grading conclusions, but to accurately identify the root causes of problems in shoreline ecosystems, predict the overall changing trends of shorelines and rivers and lakes, provide a scientific basis for shoreline protection and restoration, comprehensive river and lake management, and routine supervision, and promote the continuous and positive development of river and lake ecosystems.

[0003] The current "Technical Guidelines for River and Lake Health Assessment (Trial)" has clearly defined the evaluation indicators and calculation methods for the natural conditions of riverbanks, laying a standardized foundation for river and lake health monitoring. However, current practical assessment work still has significant shortcomings: the calculation of various indicators mainly relies on manual field surveys or static analysis of single remote sensing data, which is difficult to adapt to the modern needs of dynamic river and lake monitoring. On the one hand, manual surveys suffer from high costs, significant terrain limitations, and limited coverage, making it impossible to achieve full coverage of long-distance and complex river sections; on the other hand, traditional assessment cycles last for 1-2 months, resulting in serious timeliness deficiencies and failing to meet the requirements of routine dynamic monitoring and timely completion of the assessment of rivers and lakes in the census list.

[0004] Furthermore, while evaluation models based on single remote sensing data have improved efficiency to some extent, they have fundamental limitations: First, they treat shoreline indicators such as slope and vegetation cover in isolation, ignoring the spatiotemporal correlations and synergistic effects among various factors, resulting in evaluation results that lack systematicity and scientific rigor. Second, relying on single-point-of-time data for analysis fails to reflect the dynamic evolution of shoreline ecology and makes it difficult to identify hidden problems such as "gradual degradation" of the shoreline. Third, the evaluation results only reflect grading conclusions, lacking diagnosis of the "causes" of shoreline degradation and early warning of deterioration risks, failing to accurately align with actual governance needs, and making it difficult to effectively support the river and lake chief system assessment and forward-looking, precise management decisions.

[0005] Therefore, it is imperative to construct an automated, intelligent, and in-depth shoreline natural condition assessment system. This system will effectively address the pain points of traditional shoreline assessments, such as low efficiency, narrow coverage, superficial analysis dimensions, and weak decision support. It will achieve the goals of multi-source data fusion, spatiotemporal dynamic analysis, precise risk early warning, and targeted governance, shifting river and lake shoreline assessment from passive result output to proactive risk early warning. This will provide solid scientific support for the systematic protection and restoration of river and lake ecological shorelines, comprehensively strengthening the basic defense line for river and lake health. Summary of the Invention

[0006] To achieve the above objectives, this application provides a shoreline health assessment method based on spatiotemporal coupling diagnosis, comprising the following steps: The raw data of river shoreline status is acquired, and the raw data of river shoreline status is preprocessed and calibrated to generate a multi-source dataset of river shoreline status. The raw data of river shoreline status is a three-dimensional monitoring system consisting of satellite remote sensing data, UAV LiDAR data and ground verification data. The dataset is a comparable multi-source dataset generated after preprocessing and calibration. Based on a multi-source dataset of river shoreline conditions, time-series indicators of shoreline characteristics were extracted. These indicators were then combined with annual precipitation anomaly Pt, flood season flow level Ft, and human activity variable Ht to construct a spatiotemporal characteristic data cube. The time-series shoreline characteristics indicators included: average slope score. Vegetation coverage score Elevation difference score between the top and bottom of the slope Matrix characteristic spectral score Slope toe elevation difference score and the baseline year vegetation cover score Generate time-series features, represented as: , , , , , , , and ; Features are extracted from spatiotemporal symmetry data cubes to diagnose temporal causal relationships and analyze the spatial heterogeneity of vegetation bank stabilization effects, predict temporal change trends, delineate risk warning zones, and achieve dynamic risk warning. Evaluation conclusions are generated based on dynamic risk warnings, and these conclusions cover scoring, diagnosis, warnings, and countermeasures.

[0007] The preprocessing includes: performing radiometric and geometric corrections on satellite remote sensing data, and converting UAV LiDAR data into a digital elevation model after point cloud denoising. The calibration process includes: achieving three-level spatiotemporal alignment of satellite remote sensing data, UAV LiDAR data, and ground verification data through affine transformation; unifying the coordinate systems of satellite remote sensing data, UAV LiDAR point clouds, and ground verification data; aligning the high-frequency data of ground verification data with the low-frequency data of satellites and UAVs to the same time node; identifying outliers, removing only equipment error anomalies, and retaining risky extreme values.

[0008] The construction of the spatiotemporal vital sign data cube includes the following steps: The target shoreline is divided into evaluation units, and each evaluation unit is numbered i (i=1,2,...,n). Standardized scores of shoreline vital signs time series indicators were extracted one by one based on the time series t. Extract time-series indicators of shoreline vital signs and watershed environmental variables; The temporal characteristics of shoreline vital signs and watershed environmental variables are generated, and the comprehensive value of shoreline stability is calculated. ; Constructing the matrix of independent variables ; Construct a spatiotemporal vital sign data cube.

[0009] Among them, the extracted shoreline vital signs time series indicators include: Average slope score of the bank slope The average slope of the bank was calculated based on the digital elevation model and assigned a value according to the guideline classification standard. Vegetation Coverage Score The assignment methods include: calculating vegetation cover using a pixel-based binary model, and calculating the coefficient of variation based on three years of cover data. According to the coefficient of variation The value of the vegetation coverage score Assignment; Elevation difference score between the top and bottom of the slope The elevation difference H between the top and bottom of the bank slope is extracted and assigned based on the digital elevation model; Matrix characteristic spectral score The assignment method includes: extracting matrix feature spectra from hyperspectral data of satellite remote sensing data, classifying and identifying matrix types through support vector machines, wherein the matrix types include bedrock, soil and rock, clay and non-clay, and assigning values ​​based on matrix types; Slope toe elevation difference score The assignment method includes: calculating the slope toe elevation difference ΔZ between the two digital elevation models, and scoring according to the range of the slope toe elevation difference ΔZ; Baseline year vegetation cover score The vegetation coverage value is assigned based on the calculated baseline year.

[0010] The construction of the spatiotemporal vital sign data cube includes the following steps: Construct a time-series vital signs data matrix for each shore segment i. , represented as: ; All sections of the shore Together with the corresponding basic information, time series information, and quality control information, they constitute a spatiotemporal characteristic data cube of the target river; The associated data of the spatiotemporal characteristic data cube is stored in the PostGIS spatiotemporal database, and the DEM raster and satellite image data are stored in MinIO cloud storage and associated with the PostGIS spatiotemporal database through file paths.

[0011] Furthermore, the diagnosis of temporal causality includes the following steps: Five univariate time series were generated based on five hysteresis stability indices for each shore segment i. , where j=1...5; Using the ADF test model Implement the ADF test; Determine the lag order p; By combining the F-test and Bonferroni correction, the significance level was ensured to reach α=0.005; Generating a shore-level causal network includes: defining the nodes of the shore-level causal network as hysteresis stability indices, wherein the hysteresis stability indices are... , , , and Draw directed edges between nodes based on significant Granger causes; The nodes and directed edges are shown in the causal network graph of shoreline stability.

[0012] Furthermore, the spatial heterogeneity analysis of vegetation bank stabilization effect was achieved through a geographically and temporally weighted regression model, including the following steps: Define the comprehensive value of shoreline stability Define the independent variable matrix as the dependent variable. All variables are independent variables; calculate the latitude and longitude of the center point of each shore segment i. As spatial coordinates, the standardized time variable t is used as a time coordinate; Preprocess the set of independent variables; Construct a spatiotemporal weight matrix to quantify the composite proximity of geographical distance and temporal distance, and determine the sample weight allocation for local regression; Based on cross-validation, the bandwidth parameters of the geographic and time-weighted regression model are optimized. Local regression coefficients are iteratively calculated and t-tests are performed to estimate the strength of spatiotemporal nonstationary relationships and test the statistical significance of local estimates; this includes: for each shore segment i, based on the spatiotemporal weight matrix... We construct a weighted least squares regression model, which is expressed as: ,in, For the dependent variable vector, Design a matrix for the independent variable. This is a vector of local regression coefficients. Let be the residual vector; where, Let be the comprehensive value of the dynamic stability of the shoreline for the i-th shore segment at time t. For the local intercept coefficient of shore segment i, for Local coefficients Local coefficients, ... for Local coefficients The coefficient iterative solution is expressed as: ,in, Let be the estimated local regression coefficient value for shore segment i. The spatiotemporal diagonal weight matrix corresponding to shore segment i; calculate the t-statistic of the k-th independent variable at shore segment i. ,according to Determine significance; Based on a geographically and temporally weighted regression model, a heat map of the spatial distribution of vegetation bank stabilization effect is generated to predict the temporal trend of change.

[0013] Furthermore, generating a spatial distribution heatmap of the vegetation bank stabilization effect includes the following steps: For all sections of the shore Standardization processing The target riverbanks are classified into high-effect zones, medium-effect zones, and low-effect zones. Using the spatial vector map of the target river shoreline as the base map, the classification results of each shore segment are colored, auxiliary information is overlaid, and key areas are marked to realize the creation of a spatial distribution heat map of the vegetation shoreline stabilization effect.

[0014] Furthermore, in implementing the cross-validation, a test set is constructed sequentially using the dependent and independent variables of each shoreline segment, while the dependent and independent variables of the remaining shoreline segments serve as the training set. A geographically and temporally weighted regression model is trained based on the training set to predict the results of the test set. Value, calculate the prediction error CV; The error formula for cross-validation is defined as follows: ,in, The total cross-validation error is given by bandwidth h. Let i be the measured comprehensive value of shoreline stability for shoreline segment i. For the bandwidth h of the shore segment i Predicted value, where n is the total number of shore segments; Define a candidate range for bandwidth h, calculate the CV value for each different h, and select the h with the smallest CV value as the optimal bandwidth.

[0015] The implementation of dynamic risk early warning includes the following steps: Identifying dominant factors based on causal network graphs of shoreline stability; The changing trends of the dominant factors were analyzed using the Mann-Kendall trend test and Sen's Slope rate of change estimation. Determine the triggering conditions for multi-level dynamic risk warnings and perform spatial mapping to achieve dynamic risk warnings; including: determining the triggering conditions for multi-level dynamic risk warnings and performing spatial mapping based on the warning level of the shoreline; Output an early warning report, which includes: a dynamic risk warning list for the shoreline, a spatial distribution map of the risk warnings, and specific intervention measures.

[0016] According to this invention, the causal relationships and spatial heterogeneity among indicators reflecting the natural condition of riverbanks can be demonstrated. Based on characteristic indicators, a spatiotemporal characteristic data cube combining time, space, and indicators is constructed to support efficient querying and coupled analysis. A precise risk early warning mechanism is established to identify dominant factors, conduct trend testing and quantification thresholds of dominant factors, identify riverbank sections that are currently in good condition but pose future risks, and achieve early intervention, breaking through the limitations of traditional ex-post evaluation. Ultimately, this addresses the industry pain points of low efficiency, narrow coverage, shallow analysis, and weak support in the evaluation of the natural condition of riverbanks. Attached Figure Description

[0017] Figure 1 This is a schematic diagram of the steps of a shoreline health assessment method based on spatiotemporal coupling diagnosis provided in an embodiment of the present invention. Detailed Implementation

[0018] This invention is based on multi-source time-series remote sensing data from air, space, and ground. First, it constructs a dual-matrix spatiotemporal characteristic data cube to achieve data standardization and unified management. Then, through the coupled analysis of Granger causality test (time dimension) and GTWR model (spatial dimension), it finds the temporal driving relationship and spatial difference pattern of shoreline degradation. Finally, it identifies the dominant factors, delineates dynamic risk warning zones, and outputs a comprehensive evaluation conclusion that can be implemented, achieving a breakthrough from static post-event evaluation to dynamic pre-event warning.

[0019] The specific implementation of the present invention will now be described in detail with reference to the accompanying drawings.

[0020] Figure 1 A schematic diagram illustrating the steps involved in evaluating the natural condition of riverbanks is provided, including the following steps: Step S100: Obtain raw data on river shoreline status, preprocess and calibrate the raw data on river shoreline status, and generate a multi-source dataset of river shoreline status.

[0021] The raw data on river shoreline status is comprised of a three-dimensional monitoring system integrating satellite remote sensing data, UAV LiDAR data, and ground verification data. After preprocessing to eliminate errors and calibration to achieve spatiotemporal alignment, a comparable, high-precision multi-source dataset is generated. Specifically, different data collection requirements exist when acquiring raw river shoreline status data: 1) Satellite remote sensing data: Priority should be given to multispectral images from Gaofen-6 and Gaofen-1, and hyperspectral images from Gaofen-5, supplemented by 10m resolution images from Sentinel-2; establish a data acquisition cycle (e.g., 5 days, 10 days, determined according to the actual situation of the study area) to ensure data continuity; the data acquisition period should cover the evaluation baseline year and the previous 3 consecutive years, with 2-3 data acquisitions each year during the vegetation growing season (May-September), cloud cover ≤15%, and no shoreline obstruction; hyperspectral images should cover the 400-900nm matrix-sensitive band.

[0022] 2) UAV LiDAR data: The data acquisition equipment is a UAV equipped with a Livox Avia LiDAR sensor and a ranging accuracy of approximately ±2cm. Data is collected by flying once before the flood season each year (e.g., in June) and once after the flood season (October). The flight altitude is 120-150m, and the point cloud density is ≥150 points / m². 2 The heading overlap rate is ≥80%, and the lateral overlap rate is ≥60%. The data collection route includes the main route parallel to the shoreline and the verification route perpendicular to the shoreline, ensuring full coverage of the top, toe, and surface of the shore slope.

[0023] 3) Ground verification data: Monitoring points are adaptively deployed according to watershed morphology and shoreline complexity. The monitoring range of the monitoring points covers the upstream-middle-downstream of the riverbank, steep slopes-gentle slopes, and natural shorelines-artificially intervened shorelines. Along the distribution route, multiple (e.g., 3-5) verification points are set up per 1km shoreline section to determine the distribution density. The total number of monitoring points is determined based on the watershed area, such as ≥1000km. 2 The total number of monitoring points in the entire watershed area is ≥30 or ≥50; vegetation quadrat (1m×1m) data is collected once a month; slope morphology data is measured once a quarter; and matrix sample data (including soil, rock type, clay content, and erosion resistance coefficient) is collected once every six months.

[0024] Furthermore, the raw data on river shoreline status is preprocessed to eliminate error sources such as atmospheric interference and data noise; and calibration processing such as spatiotemporal registration strategy is adopted to achieve accurate alignment of satellite imagery, UAV and ground data, so that the horizontal error is ≤1m and the vertical error is ≤0.3m, ensuring strict comparability of data from different sources and at different times on the same shoreline.

[0025] Preprocessing includes: 1) Satellite remote sensing data preprocessing, including radiometric and geometric correction: The hyperspectral data were denoised using Savitzky-Golay filtering with a denoising window size of 5-11 points and a polynomial order of 2. Atmospheric correction was then performed using the radiative transfer equation method, with the correction formula as follows: , in To correct the surface radiance, Atmospheric transmittance, For surface reflectance, This refers to the solar irradiance outside the atmosphere. For sky optical path radiation; For satellite remote sensing data with ≥10 ground control points (such as bridges and road intersections) within the target watershed, the quadratic polynomial transformation method is used to correct image distortion, so that the error after correction is ≤0.5 pixels.

[0026] 2) Preprocessing of UAV LiDAR data, including point cloud denoising and conversion into a digital elevation model (DEM). When denoising point clouds, statistical filtering is used to remove outliers, as shown below: , among which, among which Let be the three-dimensional distance from the i-th LIDAR point to the neighborhood center. , , ) represents the coordinates of the i-th LiDAR point. ) is the mean of the coordinates of all points in the neighborhood, when ( When the distance to neighboring points is equal to the standard deviation, the point is identified as an outlier and removed.

[0027] When converting to a digital elevation model (DEM), Kriging interpolation is used to transform the denoised point cloud data into a DEM with a resolution ≤0.5m, represented as: ,in, For the elevation of the interpolation point, Given the elevation of a point, is the interpolation weight, and n is the number of known points participating in the interpolation.

[0028] Calibration processes include spatiotemporal registration, coordinate system unification, time alignment, and outlier handling. 1) Time registration refers to the three-level spatiotemporal alignment of satellite remote sensing data, UAV LiDAR data, and ground verification data through affine transformation; spatiotemporal alignment is achieved through the affine transformation formula. The implementation is as follows: (x, y) represents the original data coordinates, and (x', y') represents the registered coordinates. , , , , , The transformation parameters are obtained by solving for ≥10 ground control points within the target watershed.

[0029] 2) When unifying the coordinate system, satellite remote sensing data, UAV LiDAR point clouds and ground verification data are converted to the CGCS2000 National Geodetic Coordinate System / Gauss-Kruger Projection 3° zone; the images are in GeoTIFF format, the point clouds are in LAS 1.4 format, and the monitoring and analysis results are stored in the PostGIS spatiotemporal database.

[0030] 3) Time alignment is achieved using linear interpolation: Data from different acquisition frequencies are time-matched using linear interpolation, aligning high-frequency data from ground verification data with low-frequency data from satellites and UAVs to the same time node, thus meeting the time series analysis requirements of the target watershed.

[0031] 4) Outlier handling refers to: using median absolute deviation (MAD) to identify outliers, removing only outliers caused by equipment errors, and retaining extreme values ​​of risks such as shoreline erosion and collapse to improve the robustness of the model.

[0032] The raw data on river shoreline status, after the above processing, constitute a multi-source dataset of river shoreline status.

[0033] Step S110: Extract time-series indicators of shoreline characteristics based on multi-source datasets of river shoreline status, and construct a spatiotemporal characteristic data cube by combining annual precipitation anomaly Pt, flood season flow level Ft, and human activity variable Ht. Extracting time-series indicators of shoreline vital signs includes the following steps: 1) Divide the target shoreline into evaluation units: Referring to the current guidelines (i.e., the Technical Guidelines for River and Lake Health Assessment (Trial)), the requirements are as follows: straight shoreline section 1km / segment, curved / steep slope shoreline section 500m / segment, artificial shoreline section 200-500m / segment. Each evaluation unit corresponds to a shoreline section numbered i (i=1,2,...,n). 2) Define the time series: Extract standardized scores (0-100 points) of indicators one by one based on the time series t. The time span of the time series t must include at least 3 years and include data before and after the flood season. 3) Extract time-series indicators of shoreline characteristics and watershed environmental variables; The watershed environmental variables include annual precipitation anomaly Pt and flood season flow level Ft extracted from satellite remote sensing data; the human activity variable Ht includes the intensity of shoreline human intervention based on hyperspectral identification, specifically composed of data such as the proportion of dike construction and sand mining activities; This step extracts two independent vegetation indicators, used for dynamic stability calculations and guideline compliance scoring, respectively. The dynamic time-series indicators include the average slope score. Vegetation coverage score Elevation difference score between the top and bottom of the slope Matrix characteristic spectral score Slope toe elevation difference score Static benchmark indicators include the vegetation cover score for the benchmark year. .

[0034] The method for extracting time-series indicators of shoreline vital signs is as follows: ① Average slope score of the bank slope The average slope of the bank was calculated based on the digital elevation model and assigned a value according to the guideline grading standard. The slope thresholds were: 100 points for 15°, 75 points for 30°, 25 points for 45°, and 0 points for 60°. For slopes above 15°, the score was calculated using linear interpolation based on the corresponding threshold range. ②Vegetation coverage score First, the vegetation cover is calculated using a pixel-based binary model, represented as follows: ,in For vegetation coverage, For the current pixel normalized vegetation index, The NDVI value of bare soil. For dense vegetation NDVI value, and =0.7; then, the coefficient of variation was calculated based on the 3-year coverage data. Vegetation coverage score based on CV value The scoring is conducted, with specific CV thresholds as follows: 5% = 100 points, 25% = 75 points, 50% = 25 points, and 75% = 0 points; for CVs above 5%, scores are calculated using linear interpolation based on the corresponding threshold range. ③Elevation difference score between the top and bottom of the slope The elevation difference H between the top and bottom of the slope is extracted based on the digital elevation model and assigned according to the guide grading standard. The threshold values ​​for the elevation difference H are: 100 points for 1m, 75 points for 2m, 25 points for 3m, and 0 points for 5m. When the elevation difference H is greater than 1m, the score is calculated using linear interpolation. ④ Matrix characteristic spectral score Matrix characteristic spectra were extracted from hyperspectral data based on satellite remote sensing data. Support Vector Machine (SVM) was used to classify and identify matrix types, and the accuracy of classification and identification was calibrated using ground validation data to ensure a Kappa coefficient ≥ 0.85. Matrix types included bedrock, soil-rock, clay, and non-clay. The stability of each matrix type was as follows: bedrock > soil-rock > clay > non-clay. Matrix characteristic spectral scores were assigned to matrix types according to the guidelines. Scoring: Bedrock is 100 points, soil and rock is 75 points, clay is 25 points, and non-clay is 0 points.

[0035] ⑤ Slope toe elevation difference score First, calculate the slope toe elevation difference ΔZ between the two digital elevation models. Refer to the guidelines and score the slope toe elevation difference ΔZ. The threshold values ​​for elevation difference ΔZ are: 0.1m = 100 points, 0.3m = 75 points, 0.5m = 25 points, and >0.7m = 0 points. When the elevation difference ΔZ is greater than 0.1m, linear interpolation is used to calculate the score. ⑥ Baseline year vegetation cover score : Adopting the vegetation coverage score A consistent pixel-based binary model is used to calculate the base year vegetation cover. The values ​​are assigned according to the guideline grading standards. The base year vegetation cover thresholds are: 75% = 100 points, 50% = 75 points, 25% = 50 points, 5% = 25 points, 5% = 5 points, and 0% = 0 points. When the coverage is between the grading thresholds, the score is calculated using linear interpolation.

[0036] 4) The temporal characteristics of the generated shoreline feature time series indicators are expressed as follows: , , , , and Where i is the shore segment number and t is the time sequence number; The comprehensive value of shoreline stability is calculated and expressed as: =( + + + + ) / 5; The time-series characteristics of the annual precipitation anomaly P, flood season flow level F, and human activity variable H are expressed as follows: , and .

[0037] 5) Construct an independent variable matrix based on time series characteristics. ; Independent variable matrix The independent variables include core independent variables, lagged stability indicators, watershed environmental variables, and human activity variables; the core independent variables include Hysteresis stability indicators include , , , , ; Furthermore, the dataset of independent variables is standardized to eliminate the influence of dimensions, preparing for subsequent multicollinearity diagnosis.

[0038] Independent variable matrix Represented as: 6) Construct a spatiotemporal vital sign data cube; First, construct the temporal characteristic data matrix for each shore segment i. , represented as: ; Then, all the shore sections Together with the corresponding basic information (shore section number, geographical location, length, shoreline type), time series information, and quality control information (data accuracy, completeness), it constitutes the spatiotemporal characteristic data cube of the target river; A hybrid architecture combining a spatiotemporal database and cloud storage is adopted. The associated data of the spatiotemporal characteristic data cube is stored in the PostGIS spatiotemporal database, and a composite index is established using the shoreline ID and timestamp to achieve millisecond-level queries. Large-capacity data such as DEM raster and satellite imagery are stored in MinIO cloud storage and associated with the PostGIS spatiotemporal database through file paths.

[0039] At this point, the spatiotemporal symmetry data cube is completed and can be used for subsequent diagnosis of temporal causal relationships between indicators and spatial heterogeneity analysis of vegetation shoreline stabilization effects, to determine the natural condition of the river shoreline.

[0040] Step S120: Extract features from the spatiotemporal symptom data cube, perform temporal causal relationship diagnosis, analyze the spatial heterogeneity of vegetation bank stabilization effect to predict temporal change trends, delineate risk warning zones, and realize dynamic risk warning; 1. Diagnosis of temporal causality The temporal causality diagnosis is based on the temporal characteristics of five hysteresis stability indices for each shore segment i, generating five sets of univariate time series. Where j=1...5, corresponding to , , , , ; t=1...n, n≥6, i.e., 3 years × 2 seasonal nodes). The diagnosis of temporal causality includes the following steps: 1) Using an ADF test model that includes a constant term and a trend term to test... Implement the ADF test; whereby the ADF test model is expressed as: , in, For first-order difference operators, For univariate time series, For constant terms, For trend term coefficients, For the lagged term coefficient, Let p be the iteration ordinal number of the lag term, and p be the lag order. The coefficients of the difference lag term, This is the random error term; Set the criteria for determining the null hypothesis H0 and the alternative hypothesis H1. The criteria include: calculating the ADF test statistic under significance levels of 1%, 5%, and 10%, and comparing the ADF test statistic under each condition with the preset critical value. If the ADF test statistic is less than the critical value and the probability value P is less than 0.05, then H0 is rejected and the sequence is determined to be stationary.

[0041] If a time series fails the ADF stationarity test, a first-order differencing transformation is performed on it, which is represented as: The ADF test was re-executed on the series after the first-order difference transformation. After all series were stationary, the Johansen cointegration test was used to verify the long-term equilibrium relationship between variables and avoid spurious regression.

[0042] 2) Determine the lag order p based on the AIC / BIC criterion; First, construct the VAR model: using 5 stationary time series for each shore segment as variables, construct the VAR(p) model, which is expressed as: , in, This constitutes a lagged stability index, A k is a 5×5 lag coefficient matrix, where p is the lag order; Secondly, based on the length of the time series data (n≥6), the range of the candidate lag order p is set. For example, to avoid p being too large and causing insufficient degrees of freedom, which would affect the test power, p is set to 1-3.

[0043] Then, the Akaike Information Criterion (AIC) is defined as follows: Where k = 5 × 5 × p + 5, and L is the model likelihood function value; Define the Bayesian Information Criterion (BIC) as follows: , where n is the length of the time series data.

[0044] Calculate the AIC and BIC values ​​for p=1, 2, and 3 respectively, and select the p value corresponding to the minimum AIC and BIC values ​​as the lag order p; if the AIC and BIC values ​​conflict, prioritize the lag order p with the minimum BIC.

[0045] 3) Combine the F-test and Bonferroni correction to ensure that the significance level reaches α=0.005; When performing the F-test, the null hypothesis H0 is set: X is not a Granger cause of Y, and the alternative hypothesis H1 is set: X is a Granger cause of Y, where X and Y are any two indicators from the lagged stability indices. The constrained model (with H0 constraint) and the unconstrained model are constructed to calculate the F-statistic, expressed as: ,in, To constrain the sum of squared residuals in the model, is the sum of squared residuals of the unconstrained model, p is the lag order, and n is the time series length; Set the significance level α = 0.05, and obtain the critical value from the F-distribution table. ; If the calculated F-statistic > If the probability value P < 0.05, then reject H0 and determine that X is a Granger cause of Y; Bonferroni was used for multiple test correction, including pairwise tests of hysteresis stability indices, with a certain number of tests. ; Corrected significance level ; A significant causal relationship between X and Y is confirmed only when the probability value P of the F-test is less than 0.005.

[0046] 4) Generate shoreline-level causal networks: Define the nodes of the shore-segment-level causal network as hysteresis stability indices ( , , , and When defining nodes, label each node with an indicator name and a shoreline ID; Define the edges of the shore-level causal network: If index X is a significant Granger cause of Y (P < 0.005), then draw a directed edge from X to Y; the thickness of the edge is graded according to the probability value P (thick line when P < 0.001, thin line when 0.001 ≤ P < 0.005). Furthermore, the riverbank sections are summarized and classified, namely: causal relationships are summarized according to riverbank section type (natural riverbank section / artificial riverbank section) and watershed division (upstream / midstream / downstream), and the proportion of riverbank sections with each type of causal relationship is statistically analyzed.

[0047] A causal network diagram of shoreline stability is generated according to the type of shoreline segment. The node (indicator) and directed edge (causal relationship) are displayed in the causal network diagram. A statistical table of the proportion of causal relationships is generated to clarify the prevalence of various driving relationships and provide a basis for subsequent identification of dominant factors.

[0048] 2. Spatial Heterogeneity Analysis of Vegetation Bank Stabilization Effect Spatial heterogeneity analysis of vegetation bank stabilization effect was performed using a geographically and temporally weighted regression model (GTWR), with all independent variables directly taken from the independent variable matrix. Specifically, it includes the following steps: 1) Define the comprehensive value of shoreline stability Define the independent variable matrix as the dependent variable. All variables are independent variables, namely: the baseline year vegetation cover score (PCr), the lagged stability index, watershed environmental variables, and human activity variables are independent variables; the latitude and longitude of the center point of each shoreline i are obtained based on the geometric center of the evaluated shoreline. As spatial coordinates, the standardized time variable t is used as a time coordinate; 2) Preprocess the set of independent variables, including: calculating the variance inflation factor (VIF) of all independent variables and performing VIF verification. If VIF > 10, use stepwise regression to remove variables with high VIF, or use PCA to transform highly correlated variables into uncorrelated principal components so that the set of independent variables is free from serious multicollinearity. 3) Construct a spatiotemporal weight matrix to quantify the composite proximity of geographical distance and temporal distance, and determine the sample weight allocation for local regression; In practice, the spatiotemporal distance is first calculated using a weighted summation method. , represented as: , in, The spatial weights are adaptive (automatically optimized in the range of 0.6-0.8 based on the spatial heterogeneity of the watershed). Let i be the spatial distance between shore segment i and shore segment j. The maximum spatial distance across the entire shoreline. This represents the time difference between shoreline segment i and shoreline segment j. The maximum time difference; The spatial distance between shore segment i and shore segment j is calculated using the spherical distance formula, and is expressed as follows: , where R is the Earth's radius (6371km). , Longitude of shoreline segment i; The method for calculating the time difference between shore segment i and shore segment j is as follows: .

[0049] Secondly, a Gaussian kernel function is used to construct the spatiotemporal weight matrix. Spatiotemporal weight matrix Each element Let the weight of shore segment j with respect to i be expressed as: Where h is the bandwidth parameter, and .

[0050] 4) Based on cross-validation, optimize the bandwidth parameters of GTWR to balance its goodness of fit and complexity, avoiding overfitting or underfitting. Specifically, construct a test set for the dependent and independent variables of each shore segment sequentially, and use the dependent and independent variables of the remaining shore segments as the training set. Train the GTWR model based on the training set and predict the results for the test set. The value is used to calculate the prediction error CV, in order to achieve cross-validation based on data from all shorelines; Define the error of cross-validation The formula is, ,in, The total cross-validation error is given by bandwidth h. Let i be the measured comprehensive value of shoreline stability for shoreline segment i. For the bandwidth h of the shore segment i The predicted value is n, which is the total number of shoreline segments. Based on this, the candidate range of bandwidth h is set (e.g., standardized to 0.1-5.0 according to the total shoreline length). The CV value corresponding to different h is calculated one by one. The h corresponding to the smallest CV value is selected as the optimal bandwidth. For example, the optimal bandwidth of a certain river is h=1.2, and the CV is 128.3, which is the minimum value.

[0051] 5) Perform iterative calculation and t-test of local regression coefficients to estimate the strength of spatiotemporal nonstationary relationships and test the statistical significance of local estimates; specifically including: For each shore segment i, based on the spatiotemporal weight matrix We construct a weighted least squares regression model, which is expressed as: , in, For the dependent variable vector, Let be the comprehensive value of the dynamic stability of the shoreline for the i-th shore segment at time t; The independent variable design matrix for shore segment i is directly constructed from the independent variable matrix Vi.

[0052] This is a vector of local regression coefficients; For the local intercept coefficient of shore segment i, for Local coefficients for Local coefficients, ... for Local coefficients for Local coefficients for Local coefficients The residual vector; The coefficient iterative solution is expressed as: ; in, Let be the estimated local regression coefficient value for shore segment i. For the spatiotemporal diagonal weight matrix corresponding to shore segment i, during the iteration process, the coefficient difference between two adjacent iterations is ≤10. -6 When the iteration ends, proceed. Calculate the t-statistic of the k-th independent variable at shore segment i. , represented as: ,in, This represents the estimated local regression coefficient of the k-th independent variable. Let be the standard error of the local regression coefficient of the k-th independent variable; according to Determine the local effect of the independent variable on the dependent variable and assess its significance, including: setting α = 0.05, if... >t 0.025 If the coefficient is (n-7) and the probability value P < 0.05, then the coefficient is considered significant.

[0053] At this point, GTWR can output a spatial distribution heat map of the vegetation bank stabilization effect, predicting the temporal trend of change.

[0054] 6) When outputting the spatial distribution heat map of vegetation bank stabilization effect, the following steps are included: First, for all sections of the shore Standardization was performed, and the target riverbanks were classified based on statistical analysis of the entire sample: high-efficiency areas. The first 20%, the medium-effect zone The middle 60% is the low-efficiency zone. The last 20%); Then, a spatial heat map is drawn: using the spatial vector map of the target river shoreline as the base map, each 1km shoreline segment is divided into... The classification results are colored (e.g., high-effect areas in red, medium-effect areas in yellow, and low-effect areas in blue); auxiliary information such as shoreline type and watershed zoning is overlaid to mark key areas. Based on this, supporting analyses can be conducted to explain the causes of spatial differences in vegetation shoreline stabilization effects (e.g., topographic slope, matrix type, human activity disturbance, etc.) and to propose targeted suggestions (e.g., prioritizing vegetation ecological restoration in high-effect areas and combining matrix reinforcement projects in low-effect areas).

[0055] 3. By comprehensively judging the dominant factors, risk warning zones are delineated to achieve dynamic risk warning; By using a heat map of the spatial distribution of vegetation bank stabilization effects, the time and area of ​​risk requiring early warning are determined. Based on this, dynamic risk warnings are issued according to a comprehensive assessment of dominant factors, specifically including the following steps: 1) Identifying dominant factors based on the causal network diagram of shoreline stability: In the Granger causality test, indicators that establish a significant causal relationship (probability P < 0.05) with more than two shoreline characteristic time-series indicators, and whose causal relationship accounts for ≥ 50% across all shoreline segments, are designated as core dominant factors; in the GTWR model, for The absolute value of the regression coefficient is the largest ( Indicators with a value ≥0.6 and significant (P<0.05) are set as key dominant factors; for each shoreline segment, the intersection of the core dominant factor and the key dominant factor is taken. If the intersection is empty, the key dominant factor is taken. If there are multiple indicators in the intersection, they are sorted by the absolute value of the coefficient, and the first two elements are taken as dominant factors.

[0056] 2) The changing trends of dominant factors were analyzed using the Mann-Kendall trend test and Sen's slope rate of change estimation: In the Mann-Kendall trend test analysis, the time series data Z of the dominant factor are used. t (t=1 to n, n≥6) are used as inputs. The null hypothesis H0 is that the sequence has no trend, and the alternative hypothesis H1 is that the sequence has a significant upward or downward trend. Calculate the rank of the sequence S, denoted as: ,in The sign function (positive is 1, negative is -1, zero is 0). The dominant factor is time series data; Calculate variance , Standardized statistics To determine the trend, set α=0.05. If Z>1.96 (and P<0.05), it is considered a significant increase; if Z<-1.96 (P<0.05), it is considered a significant decrease; otherwise, there is no significant trend.

[0057] When estimating the rate of change of Sen's Slope, the slope of all pairwise data pairs in the time series data of the dominant factor is calculated. , represented as: , where i and j are the position indices of the time series data points.

[0058] Furthermore, Sen's Slope estimate , which is the median of all slopes, is robust and unaffected by outliers; where Q is in points per year, and if Q=-6.2, it means that the core dominant factor has decreased by an average of 6.2 points per year.

[0059] 3) Determine the triggering conditions for multi-level dynamic risk early warning and perform spatial mapping to realize dynamic risk early warning; First, determine the triggering conditions for multi-level dynamic risk warnings (as shown in Table 1), and then perform spatial plotting based on the warning levels of the shoreline sections. Table 1 Multi-level dynamic risk warning information

[0060] When performing spatial mapping, all shore segments are traversed according to the triggering conditions in Table 1 to filter out the warning levels of each shore segment and clarify the core information and risk sources of each shore segment. The core information includes shore segment number, warning level, core dominant factor, downward trend (Z value, Sen's Slope), current score and historical quantile. The risk source is the reason for the risk caused by the decline of this indicator based on the causal network.

[0061] Simultaneously, an early warning report is output, which may include a table file containing detailed information on all warning shorelines, an image file containing an explanation of the warning level and a priority ranking of intervention recommendations, and a document containing specific intervention measures for warning shorelines with different dominant factors.

[0062] Step S130: Generate evaluation conclusions based on dynamic risk warnings, the evaluation conclusions covering scoring, diagnosis, warnings and countermeasures.

[0063] When generating evaluation conclusions, the scoring formula is first determined, including: Overall score for shoreline stability: =( + + + + ) / 5; Shoreline Natural Condition (BH) Final Score: BH = ×0.4 + ×0.6.

[0064] The calculated results are then compared with the manual scoring results to ensure that the consistency error is ≤3 points.

[0065] Furthermore, dynamic evaluation conclusions are generated, including: 1) Core findings include: total score of shoreline natural condition, warning level, main problems, and number of high-risk shoreline sections; 2) Evaluation scope and data description, including shoreline length, number of evaluated shorelines, data sources, and preprocessing accuracy; 3) Overall compliance score, for example: BH / / Total score and grade, shoreline section scoring table, BH score spatial distribution map; 4) Dynamic diagnostic analysis results, including attribution analysis of each bank segment (e.g., scores, dominant factors, causal relationships), spatial distribution map of dominant factors, and analysis of vegetation bank stabilization effect; 5) Risk warning list, including details of the warning shoreline (ID, level, warning indicators, downward trend, recommended intervention time) and risk spatial distribution map; 6) Remediation recommendations, including remediation measures based on dominant factors (matrix-dominant / vegetation-dominant / slope-dominant) and remediation priority ranking; 7) Other content, such as data quality reports, model parameter settings, and test result details.

[0066] Based on guidelines, this invention defines multiple characteristic indicators. Through Granger causality tests reflecting the temporal dimension and GTWR coupling analysis reflecting the spatial dimension, it constructs causal relationships and spatial heterogeneity among indicators reflecting the natural condition of riverbanks. Building upon these characteristic indicators, it constructs a spatiotemporal characteristic data cube combining time, space, and indicators to support efficient querying and coupled analysis. A precise risk warning mechanism is established, identifying dominant factors, conducting trend tests and quantifying thresholds for these factors, and identifying riverbank sections that are currently in good condition but pose future risks, enabling early intervention and overcoming the limitations of traditional ex-post evaluation. Ultimately, it addresses the industry pain points of low efficiency, narrow coverage, superficial analysis, and weak support in the evaluation of the natural condition of riverbanks.

[0067] The above-disclosed embodiments are merely a few specific examples of the present invention. However, the present invention is not limited thereto, and any variations that can be conceived by those skilled in the art should fall within the protection scope of the present invention.

Claims

1. A method for evaluating the health of a shoreline based on spatiotemporal coupling diagnosis, characterized in that, Includes the following steps: The raw data of river shoreline status is acquired, and the raw data of river shoreline status is preprocessed and calibrated to generate a multi-source dataset of river shoreline status. The raw data of river shoreline status is a comparative multi-source dataset generated by a three-dimensional monitoring system consisting of satellite remote sensing data, UAV LiDAR data and ground verification data, which is preprocessed and calibrated. Based on the aforementioned multi-source dataset of river shoreline status, time-series indicators of shoreline characteristics are extracted. These are then combined with annual precipitation anomaly Pt, flood season flow level Ft, and human activity variable Ht to construct a spatiotemporal characteristic data cube. The time-series shoreline characteristics indicators include lag stability indicators and a baseline year vegetation cover score. The hysteresis stability index includes: the average slope score of the bank slope. Vegetation coverage score Elevation difference score between the top and bottom of the slope Matrix characteristic spectral score Slope toe elevation difference score Based on the lag stability index and the baseline year vegetation cover score. The generation time series characteristics of annual precipitation anomaly Pt, flood season flow level Ft, and human activity variable Ht are respectively expressed as: , , , , , , , and Where i is the shore segment number and t is the time sequence number; Features are extracted from the spatiotemporal symmetry data cube to perform temporal causal relationship diagnosis and spatial heterogeneity analysis of vegetation shoreline stabilization effect, generate shoreline stability causal network diagram, predict temporal change trends, delineate risk warning zones, and realize dynamic risk warning. Evaluation conclusions are generated based on dynamic risk warnings, and these conclusions cover scoring, diagnosis, warnings, and countermeasures.

2. The shoreline health assessment method based on spatiotemporal coupling diagnosis according to claim 1, characterized in that, The preprocessing includes: performing radiometric and geometric correction on satellite remote sensing data, and converting UAV LiDAR data into a digital elevation model after point cloud denoising. The calibration process includes: achieving three-level spatiotemporal alignment of satellite remote sensing data, UAV LiDAR data, and ground verification data through affine transformation; unifying the coordinate systems of satellite remote sensing data, UAV LiDAR point clouds, and ground verification data; aligning the high-frequency data of ground verification data with the low-frequency data of satellites and UAVs to the same time node; identifying outliers, removing only equipment error anomalies, and retaining risky extreme values.

3. The shoreline health assessment method based on spatiotemporal coupling diagnosis according to claim 1, characterized in that, The construction of the spatiotemporal vital sign data cube includes the following steps: The target shoreline is divided into evaluation units, each of which is numbered i, i=1,2,...,n; Standardized scores of shoreline vital signs time series indicators were extracted one by one based on the time series t. Extract time-series indicators of shoreline vital signs and watershed environmental variables; The temporal characteristics of shoreline vital signs and watershed environmental variables are generated, and the comprehensive value of shoreline stability is calculated. ; Constructing the matrix of independent variables ; Construct a spatiotemporal vital sign data cube.

4. The shoreline health assessment method based on spatiotemporal coupling diagnosis according to claim 3, characterized in that, The extracted shoreline vital signs time-series indicators include: The average slope score of the bank slope The average slope of the bank was calculated based on the digital elevation model and assigned a value according to the guideline classification standard. The vegetation coverage score The assignment methods include: calculating vegetation cover using a pixel-based binary model, and calculating the coefficient of variation based on three years of cover data. According to the coefficient of variation The value of the vegetation coverage score Assignment; Elevation difference score between the top and bottom of the slope The elevation difference H between the top and bottom of the slope is extracted and assigned based on the digital elevation model; Matrix characteristic spectral score The assignment method includes: extracting matrix feature spectra from hyperspectral data of satellite remote sensing data, classifying and identifying matrix types through support vector machines, wherein the matrix types include bedrock, soil and rock, clay and non-clay, and assigning values ​​based on matrix types; Slope toe elevation difference score The assignment method includes: calculating the slope toe elevation difference ΔZ between the two digital elevation models, and scoring according to the range of the slope toe elevation difference ΔZ; Baseline year vegetation cover score The vegetation coverage value is assigned based on the calculated baseline year.

5. The shoreline health assessment method based on spatiotemporal coupling diagnosis according to claim 3, characterized in that, The construction of the spatiotemporal vital sign data cube includes the following steps: Construct a time-series vital signs data matrix for each shore segment i. , is represented as: ； All sections of the shore Together with the corresponding basic information, time series information, and quality control information, they constitute a spatiotemporal characteristic data cube of the target river; The associated data of the spatiotemporal characteristic data cube is stored in the PostGIS spatiotemporal database, and the DEM raster and satellite image data are stored in MinIO cloud storage and associated with the PostGIS spatiotemporal database through file paths.

6. The shoreline health assessment method based on spatiotemporal coupling diagnosis according to claim 3, characterized in that, The temporal causality diagnosis includes the following steps: A univariate time series is generated based on the hysteresis stability index of each shore segment i. Where j is the index number of the hysteresis stability index and j=1...5; Using the ADF test model Implement the ADF test; Determine the lag order p; By combining the F-test and Bonferroni correction, the significance level was ensured to reach α=0.005; Generate a shore-level causal network, including: defining the nodes of the shore-level causal network as hysteresis stability indices; and drawing directed edges between nodes based on significant Granger causes. The nodes and directed edges are shown in the causal network graph of shoreline stability.

7. The shoreline health assessment method based on spatiotemporal coupling diagnosis according to claim 6, characterized in that, The spatial heterogeneity analysis of the vegetation bank stabilization effect was achieved through a geographically and temporally weighted regression model, including the following steps: Define the comprehensive value of shoreline stability Define the independent variable matrix as the dependent variable. All variables are independent variables; calculate the latitude and longitude of the center point of each shore segment i. As spatial coordinates, the standardized time variable t is used as a time coordinate; Preprocess the set of independent variables; Construct a spatiotemporal weight matrix to quantify the composite proximity of geographical distance and temporal distance, and determine the sample weight allocation for local regression; Based on cross-validation, the bandwidth parameters of the geographic and time-weighted regression model are optimized. Local regression coefficients are iteratively calculated and t-tests are performed to estimate the strength of spatiotemporal nonstationary relationships and test the statistical significance of local estimates; this includes: for each shore segment i, based on the spatiotemporal weight matrix... We construct a weighted least squares regression model, which is expressed as: ,in, For the dependent variable vector, Design a matrix for the independent variable. This is a vector of local regression coefficients. Let be the residual vector; where, Let be the comprehensive value of the dynamic stability of the shoreline for the i-th shore segment at time t. For the local intercept coefficient of shore segment i, for Local coefficients for Local coefficients, ... for Local coefficients; the coefficient iterative solution is expressed as: ,in, Let be the estimated local regression coefficient value for shore segment i. The spatiotemporal diagonal weight matrix corresponding to shore segment i; calculate the t-statistic of the k-th independent variable at shore segment i. ,according to Determine the significance; Based on a geographically and temporally weighted regression model, a heat map of the spatial distribution of vegetation bank stabilization effect is generated to predict the temporal trend of change.

8. The shoreline health assessment method based on spatiotemporal coupling diagnosis according to claim 7, characterized in that, The output of the spatial distribution heat map of vegetation bank stabilization effect includes the following steps: For all sections of the shore Standardization processes were implemented to classify the target riverbanks into high-efficiency, medium-efficiency, and low-efficiency zones. Using the spatial vector map of the target river shoreline as the base map, the classification results of each shore segment are colored, auxiliary information is overlaid, and key areas are marked to realize the creation of a spatial distribution heat map of the vegetation shoreline stabilization effect.

9. The shoreline health assessment method based on spatiotemporal coupling diagnosis according to claim 7, characterized in that, To implement the cross-validation, a test set is constructed sequentially using the dependent and independent variables of each shoreline segment. The dependent and independent variables of the remaining shoreline segments are used as the training set. A geographically and temporally weighted regression model is trained based on the training set to predict the results of the test set. Value, calculate the prediction error CV; The error formula for cross-validation is defined as follows: ,in, Let h be the total cross-validation error with bandwidth h. Let i be the measured comprehensive value of shoreline stability for shoreline segment i. For the bandwidth h of the shore segment i Predicted value, where n is the total number of shore segments; Define a candidate range for bandwidth h, calculate the CV value for each different h, and select the h with the smallest CV value as the optimal bandwidth.

10. The shoreline health assessment method based on spatiotemporal coupling diagnosis according to claim 1, characterized in that, The implementation of dynamic risk early warning includes the following steps: Identifying dominant factors based on causal network graphs of shoreline stability; The changing trends of the dominant factors were analyzed using the Mann-Kendall trend test and Sen's Slope rate of change estimation. Determine the triggering conditions for multi-level dynamic risk warnings and perform spatial mapping to achieve dynamic risk warnings; including: determining the triggering conditions for multi-level dynamic risk warnings and performing spatial mapping based on the warning level of the shoreline; Output an early warning report, which includes: a dynamic risk warning list for the shoreline, a spatial distribution map of the risk warnings, and specific intervention measures.

Images