A method for automatically extracting a network of rills and calculating morphological indexes

Abstract

The application discloses a kind of tidal creek network automatic extraction and morphological index calculation method, it is related to data processing technical field, it includes the following steps: establishing scenario list and input data;The input data is dimensionless and is fused to generate fusion feature field;According to the statistical distribution of fusion feature field in effective analysis domain, determine the segmentation threshold, and according to the segmentation threshold, the fusion feature field is segmented, to obtain the initial tidal creek mask and optimization;Center line extraction is carried out to the optimized tidal creek mask, and tidal creek network data structure is generated;Local morphological index and network scale index are calculated.The application enables different tidal range and vegetation scene to be batch processed using the same spatial reference and calculation caliber, and forms a standardized result that can be directly used for wetland restoration effect evaluation.

Description

Technical Field

[0001] This invention belongs to the field of data processing technology, specifically relating to a method for automatic extraction and morphological index calculation of tidal channel networks. Background Technology

[0002] Tidal channels are the main water flow channels in coastal wetlands such as tidal flats, salt marshes, and mangroves. Their morphology and network organization directly affect tidal exchange, sediment transport, vegetation establishment, habitat connectivity, and wetland restoration efficiency. For coastal wetland restoration scenarios, the formation of tidal channel networks, the rationality of drainage paths, and how network morphology changes with tidal range, sediment, and vegetation conditions are important bases for evaluating restoration effects and optimizing engineering schemes. Numerical simulation is an important tool for the design, comparison, and evaluation of coastal wetland restoration schemes. To analyze the impact of factors such as tidal range, sediment conditions, and vegetation configuration on tidal channel development and wetland restoration effects, it is usually necessary to set up multiple sets of numerical simulation scenarios and compare and analyze the output results such as bed elevation, water depth, and flow velocity under different scenarios. As the number of simulation scenarios increases, it is necessary to use a unified spatial benchmark and calculation method to automate the simulation results of different scenarios to ensure the consistency and comparability of tidal channel extraction results and morphological indicators.

[0003] Existing methods for extracting tidal channel information are mostly geared towards remote sensing imagery, lidar, or digital elevation models, typically relying on spectral differences, topographic depression features, or manual thresholds to identify tidal channels. These methods are primarily used for mapping based on real surface observation data and are difficult to directly apply to the multiphysics output results generated by numerical simulations for coastal wetland restoration. They also fail to fully utilize the complex physical characteristics of tidal channels simultaneously exhibited in simulation data, such as relatively low topographic depressions, significant water depth, and strong flow velocities. Therefore, it is necessary to establish an automatic tidal channel extraction method oriented towards numerical simulation results, which can uniformly process and comprehensively identify physical fields such as bed elevation, water depth, and flow velocity.

[0004] On the other hand, numerical simulations for coastal wetland restoration typically require setting up and comparing multiple scenarios, while current post-processing methods often remain at the stage of single-scenario mapping or manual vectorization. Tidal channel masks, skeletons, nodes, links, and topographic rasters are stored in different data structures, and slight inconsistencies in scenario numbers, coordinates, resolutions, and effective analysis domains can lead to index mismatches. Even if the centerline of the tidal channel is obtained, it can usually only output single indices such as length or area, making it difficult to continuously calculate multi-scale indices such as width, depth, cross-sectional area, aspect ratio, drainage density, UPL (Unchanneled Path Length), OPL (Overmarsh Path Length), and geometric efficiency, and even more difficult to support automatic comparison of dozens or more restoration scenarios. Summary of the Invention

[0005] The purpose of this invention is to provide an automatic extraction and morphological index calculation method for tidal channel networks, which enables batch processing of different tidal ranges and vegetation scenarios using the same spatial benchmark and calculation caliber, and generates standardized results that can be directly used for wetland restoration effect evaluation.

[0006] To achieve the above objectives, the present invention adopts the following technical solution: A method for automatically extracting tidal channel networks and calculating morphological indicators includes the following steps: Establish a scenario list, read the input data for each scenario according to the scenario list, and align and crop all input data based on a unified spatial reference and effective analysis domain. The input data includes: bed elevation field, water depth field and / or flow velocity field. Within the effective analysis domain, the physical fields in the input data are dimensionless, and at least two dimensionless physical fields are fused to generate a fused feature field; the physical fields include: the bed elevation field, the water depth field, and / or the flow velocity field. A segmentation threshold is determined based on the statistical distribution of the fused feature field within the effective analysis domain, and the fused feature field is segmented according to the segmentation threshold to obtain an initial tidal channel mask; the initial tidal channel mask is then subjected to connectivity filtering and / or noise removal to obtain an optimized tidal channel mask. The optimized tidal channel mask is subjected to centerline extraction, and nodes and links are identified based on the extracted centerline to generate a tidal channel network data structure containing the attribute information of the nodes and links; the link attribute information includes link length, link width, and spatial coordinates. Based on the spatial correspondence between the tidal channel mask, the network data structure, and the bed surface elevation field, local morphological indices and network scale indices are calculated.

[0007] As a preferred embodiment of the method for automatic extraction and morphological index calculation of tidal channel networks, the step of generating the fused feature field includes: The bed elevation field is reverse normalized to obtain higher normalized characteristic values ​​in low-lying areas; the water depth field and / or the flow velocity field are forward normalized to obtain higher normalized characteristic values ​​in areas with large water depth and / or high flow velocity. The normalized bed elevation field, water depth field, and / or flow velocity field are summed pixel-by-pixel by weighted summation or by direct summation to obtain the fused feature field.

[0008] As a preferred method for automatic extraction and morphological index calculation of tidal channel networks, the segmentation threshold is dynamically determined based on the standard deviation of the fused feature field within the effective analysis domain: When the standard deviation is greater than the preset discriminant value, a lower quantile is used as the segmentation threshold; when the standard deviation is not greater than the preset discriminant value, a higher quantile is used as the segmentation threshold.

[0009] As a preferred embodiment of the automatic extraction and morphological index calculation method for tidal channel networks, the noise removal includes: deleting connected objects in the initial tidal channel mask whose number of pixels in the connected region is less than a preset pixel number threshold.

[0010] As a preferred embodiment of the automatic extraction and morphological index calculation method for tidal channel networks, the centerline extraction includes: performing skeletonization processing on the optimized tidal channel mask and pruning short branches on the skeletonization result; the nodes include: endpoints, branching points, and exit points; the links are continuous centerline segments between adjacent nodes.

[0011] As a preferred embodiment of the automatic extraction and morphological index calculation method for tidal channel networks, the local morphological index includes at least one of tidal channel width, tidal channel depth, cross-sectional area, and width-to-depth ratio. When calculating the local morphological index, the profile elevation value of the corresponding position is extracted from the bed elevation field based on the link centerline sampling point and the left and right bank endpoints in the network data structure.

[0012] As a preferred embodiment of the automatic extraction and morphological index calculation method for tidal channel networks, the network scale index includes at least one of drainage density, non-channelized path length, platform drainage characteristic path length, and geometric efficiency. Wherein, the drainage density is the ratio of the total length of all links to the area of ​​the effective analysis domain; The non-channelized path length is the distance from the non-tidal groove pixel to the nearest tidal groove pixel in the tidal groove mask; The characteristic path length of the platform drainage is obtained by fitting the exceedance probability distribution of the non-channelized path length; The geometric efficiency is the reciprocal of the product of the drainage density and the characteristic path length of the platform drainage.

[0013] As a preferred method for automatic extraction and morphological index calculation of tidal channel networks, the following steps are also included: The scenario identifier is associated with and stored in conjunction with the calculated local morphological index and network scale index, and a state tag for characterizing data integrity is generated.

[0014] As a preferred embodiment of the method for automatic extraction and morphological index calculation of tidal channel networks, the step of associating and storing the scenario identifier with the calculated local morphological index and network scale index, and generating a state marker for characterizing data integrity, further includes: Check whether the files in the scenario list exist, whether the coordinate reference system is consistent, whether the grid size is uniform, whether the units are correct, or whether the numbering mapping is complete; and generate a processing log based on the check results and the status markers.

[0015] As a preferred method for automatic extraction and morphological index calculation of tidal channel networks, the input data comes from the output results of the hydrodynamic-sediment-vegetation coupled numerical model at the target time, and the numerical model is the Delft3D model.

[0016] Beneficial effects: The automatic extraction and morphological index calculation method for tidal channel networks proposed in this invention utilizes multi-physics field fusion and the composite physical characteristics of tidal channels—namely, low-lying terrain, deep water, and strong flow velocity—for identification. Compared to a single elevation threshold method, this significantly reduces false identification of low-lying areas in tidal flats and improves the detection rate of narrow tributaries. Furthermore, the unified management of spatial reference, pixel resolution, and effective analysis domain through a scenario list ensures that all scenarios are processed under the same calculation caliber, resolving the issue of incomparable extraction results under different tidal ranges and vegetation conditions. Therefore, this invention enables batch processing of different tidal ranges and vegetation scenarios using the same spatial reference and calculation caliber, generating standardized results that can be directly used for wetland restoration effect evaluation. Attached Figure Description

[0017] Figure 1 This is a flowchart of a method for automatic extraction and morphological index calculation of tidal channel networks according to an embodiment of the present invention; Figure 2 Diagram showing the process of multi-physics fusion and tidal channel extraction; Figure 3 Diagram showing the construction process of the tidal channel centerline and network; Figure 4 A comparative diagram of tidal channel networks and drainage indicators in representative scenarios; Detailed Implementation 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.

[0018] Example: Reference Figure 1 This embodiment provides a method for automatic extraction and morphological index calculation of tidal channel networks, including the following steps: S1: Establish a scenario list, which records at least the scenario identifier, source file path and spatial parameters of each numerical simulation scenario; read the input data of each scenario according to the scenario list, and align and trim all input data based on a unified spatial reference and effective analysis domain, wherein the input data includes: bed elevation field, water depth field and / or flow velocity field. S2: Within the effective analysis domain, the physical fields in the input data are dimensionless, and at least two dimensionless physical fields are fused to generate a fused feature field; the physical fields include: the bed elevation field, the water depth field, and / or the flow velocity field; S3: Determine the segmentation threshold based on the statistical distribution of the fused feature field within the effective analysis domain, and segment the fused feature field according to the segmentation threshold to obtain an initial tidal channel mask; perform connectivity screening and / or noise removal on the initial tidal channel mask to obtain an optimized tidal channel mask; S4: Perform centerline extraction on the optimized tidal channel mask, and identify nodes and links based on the extracted centerline to generate a tidal channel network data structure containing the attribute information of the nodes and links; the link attribute information includes link length, link width, and spatial coordinates; S5: Based on the spatial correspondence between the tidal channel mask, the network data structure and the bed surface elevation field, calculate the local morphological index and the network scale index; S6: Associate and store the scenario identifier with the calculated local morphological index and network scale index, and generate a state tag to characterize data integrity.

[0019] Regarding the implementation environment and overall architecture of this embodiment, the present invention can be implemented based on a general-purpose computer or server system. This system includes at least: a central processing unit (CPU) or a graphics processing unit (GPU), system memory (for storing intermediate computation arrays and feature fields), non-volatile storage devices (for storing numerical simulation output files, scenario lists, and computation result reports), and input / output interfaces (for reading data files and outputting results). The software environment preferably employs a cross-platform combination of numerical computation and geographic information processing, such as MATLAB R2018b and later versions, Python 3.7 and later versions, and relies on open-source or commercial libraries such as NumPy (numerical computation library), SciPy (scientific computing library), GDAL (geographic raster processing library), and RivGraph (river network analysis library).

[0020] The following describes the invention in further detail using 32 standard scenarios output from an idealized Delft3D hydrodynamic-sediment-vegetation coupling model of near-natural restoration of coastal wetlands as examples. The model platform has a spatial range of 500m × 1200m, a grid resolution of 5m, and a corresponding effective analysis area of ​​600,000 m². The exported grid index ranges from 1 to 101 in the M direction and 62 to 301 in the N direction. The X and Y coordinates of the grid nodes exported from Delft3D are calculated using a 2×2 mean kernel to determine the corresponding grid center coordinates. The physical field values ​​retain the original exported values ​​and correspond one-to-one with the grid center coordinates. These specific parameters are for illustrative purposes only and do not constitute a limitation on the scope of application of the invention.

[0021] As described in step S1 above, regarding the scenario list data structure, a scenario list file (format can be CSV or Excel) is first created in a non-volatile storage device. Each row of this list corresponds to a numerical simulation scenario to be processed, and includes at least the following fields: standard_run_num (standard scenario number), source_path (source file root directory), tide_range (tidal range, in meters), veg_height (vegetation height, in meters), stem_diameter (stem diameter, in meters), max_density (maximum plant density, in plants / m²), and target_time_step (target simulation time step). The central processing unit reads this list and stores it in a table structure in system memory, while simultaneously creating a separate output subdirectory in the storage device based on standard_run_num.

[0022] Regarding physics field reading and coordinate registration: The processor locates the Delft3D output file for each scenario based on the path in the scenario list, typically in .dat or .nc format. Using the Delft3D QuickPlot interface or other equivalent NetCDF reading interfaces, the following physical raster fields at the target time are read: bed elevation field E (m); water depth field W (m); and flow velocity amplitude field V (m / s).

[0023] For negative water depth values ​​(representing exposed surfaces), the processor automatically corrects them to 0. After reading, the processor uses a 2×2 mean convolution kernel to map the physical quantity values ​​to the grid center coordinates based on the X and Y coordinate arrays of the original model grid, ensuring that all physical fields have the same spatial dimension and pixel position. If there are differences in the output grid resolution for different scenarios, the processor uniformly uses the nearest neighbor or bilinear interpolation method to resample to the standard resolution, which is 5 m in this embodiment.

[0024] Regarding the construction of the effective region mask: To avoid artificial boundaries of the model (e.g., raising the non-inlet boundary bed surface to 6 m to create a unique inlet) being misidentified as tidal channels, the processor executes boundary recognition logic: traversing the bed elevation field, pixels with elevation values ​​greater than 5.99 m are marked as invalid boundaries, and the feature values ​​of these locations are set to 0 or excluded in all subsequent calculations. Subsequently, pixels with water depths greater than 0.1616 m in the water depth field are used as valid water passage candidate areas, forming the effective region mask M_valid. This water depth threshold can be adjusted or calibrated in different embodiments according to the mesh resolution, model noise level, and representative scenario extraction effect. This threshold is only used to limit the possible range of tidal channels and is not directly used as the tidal channel boundary.

[0025] As described in step S2 above, multi-physics feature fusion is performed, firstly by field-by-field dimensionless processing: in system memory, the processor performs normalization operations on each physical grid field within the effective region M_valid. This invention employs three normalization strategies, as detailed below: (1) Formula for reverse normalization of bed surface elevation: in: The original bed elevation value (in meters) of any pixel is directly derived from the bed elevation field output by Delft3D. The maximum bed elevation value (unit: m) obtained statistically within the valid analysis domain M_valid of the current scenario; The minimum bed elevation value (unit: m) obtained statistically within the valid analysis domain M_valid of the current scenario; : The elevation feature value after reverse normalization (dimensionless, with a value range of [0,1]). The lower the elevation (the more depressed) the pixel, the closer this value is to 1, indicating that the pixel is more in line with the low-lying characteristics of a tidal ditch in terms of morphology.

[0026] Since tidal channels are relatively low-lying in shape, reverse normalization is used to make them appear as high values ​​in the fusion characteristics.

[0027] (2) Positive normalization formula for water depth: in: W: The original water depth value of any pixel (unit: m), directly from the water depth field output by Delft3D, with negative water depths corrected to 0; : The maximum water depth (in meters) obtained statistically within the effective analysis domain of the current scenario; The minimum water depth (in meters) obtained statistically within the effective analysis domain of the current scenario. : The water depth feature value after positive normalization (dimensionless, with a value range of [0,1]). The larger the water depth of a pixel, the closer this value is to 1, indicating that the pixel is more consistent with the characteristics of a tidal channel in terms of water depth.

[0028] (3) Positive normalization formula for flow velocity: in: V: The original velocity amplitude of any pixel (unit: m / s), which is synthesized from the u and v components of the depth-averaged velocity; The maximum flow velocity amplitude (unit: m / s) obtained statistically within the effective analysis domain of the current scenario. The minimum velocity amplitude (unit: m / s) obtained statistically within the effective analysis domain of the current scenario. : The positively normalized flow velocity characteristic value (dimensionless, with a value range of [0,1]). The higher the flow velocity of a pixel, the closer this value is to 1, indicating that the pixel is more consistent with the characteristics of a tidal channel in terms of hydrodynamic intensity.

[0029] To reduce the interference of outliers on the normalization results, other embodiments may employ quantile truncation normalization (e.g., truncating the 1% and 99th quantiles as upper and lower limits) or robust normalization (based on the median and interquartile range).

[0030] After feature field fusion and normalization, the processor merges the three types of feature fields pixel by pixel into a single fused feature field F(x, y). The general fusion formula is: in: The fusion feature value (dimensionless) of the pixel at coordinates (x,y). The higher the value, the more the pixel conforms to the characteristics of a tidal channel in terms of the comprehensive physical attributes of "low-lying terrain, greater water depth, and stronger flow velocity". The fusion weight of elevation features controls the proportion of importance of terrain depression features in the final judgment; The fusion weight of water depth features controls the proportion of importance of water depth features in the final judgment; The fusion weight of flow velocity features controls the proportion of importance of flow velocity features in the final judgment; The three weights satisfy the constraints: This ensures that the fused F value also falls within the [0,1] interval.

[0031] This embodiment employs a three-item equal-weight strategy, namely Since the weighted direct summation F = + + The weighted average differs from the weighted average by only a constant multiple and does not change the quantile-based sorting and segmentation results; therefore, the current embodiment uses direct summation to simplify the calculation. As an alternative implementation, the weights can be set based on verification samples or simulation phases of typical cross-sections: for example, the water depth weight can be increased during the high tide flooding phase (e.g., ...). =0.5), during the ebb tide return to the channel stage, the velocity weight can be increased (e.g., =0.5). If the input data lacks a velocity field or water depth field, the processor automatically resets the weights of the missing fields to 0 and renormalizes the remaining available features. In addition, vegetation density field, bed shear stress field, or suspended sediment concentration field can also be included in the fusion as additional features according to the same normalization logic.

[0032] As described in step S3, dynamic segmentation and optimization of the tidal trench mask. Adaptive threshold segmentation: The processor first calculates the standard deviation of the fused feature field F within the valid region M_valid. The calculation formula is as follows: in: : The fused feature value of the i-th pixel within the effective region; : The arithmetic mean of the fused feature values ​​of all pixels within the effective region; The total number of pixels within the effective area; The standard deviation of the fused feature field (dimensionless) reflects the feature contrast between the tidal channel and the tidal flat platform. The larger the value, the more distinct the distinction between the tidal channel and the background platform.

[0033] In this embodiment: when When the value is greater than 0.2, it indicates significant differences in the feature fields and a clear tidal groove structure. In this case, the processor uses the 8th percentile of the fused feature values ​​as the segmentation threshold T (i.e., ...). (in order to preserve as many fine tributaries as possible); when When the value is ≤0.2, it indicates that the difference in the feature field is small, and noise may dominate. In this case, the processor uses the 12th percentile as the segmentation threshold TT (i.e., ( ), to suppress false responses.

[0034] All pixels satisfying F(x,y)>T and located within M_valid constitute the initial tidal groove binary mask M_creek_initial (tidal groove pixels have a value of 1, and non-tidal groove pixels have a value of 0). The specific values ​​of the percentile thresholds (8% and 12%) and the discriminant value (0.2) in this embodiment are determined based on experiments using a 5 m resolution grid. As a more general alternative, the processor may also employ Otsu global thresholding, local adaptive thresholding, or a multi-threshold iterative scheme based on connectivity evaluation.

[0035] Post-mask processing optimization: To remove isolated false trenches caused by numerical noise or boundary effects, the processor performs connected component analysis on M_creek_initial (using the 8-neighborhood connectivity criterion). The number of pixels contained in each connected region is counted. , will satisfy < All connected objects are deleted (i.e., all cell values ​​in that region are set to 0). Among them: : The number of pixels contained in a connected region (unit: pixels); The preset area threshold is used to calculate the number of pixels (in units). In this embodiment... =10.

[0036] This pixel threshold should be converted to an actual area based on the grid resolution: for a 5 m resolution, 10 pixels correspond to 250 m². In practical applications, It should be based on the grid resolution Δx and the minimum true tidal channel area. According to the formula = / (Δ The conversion is determined. In this embodiment, considering that the model is an idealized regular mesh and the entrance is unique, the processor does not perform hole filling, morphological opening and closing operations, or entrance connectivity filtering to maintain the natural topology of the tidal channel. However, when dealing with real complex terrain or models with multiple entrances and exits, these post-processing operations can be selected as preferred enhancement methods.

[0037] The optimized trench mask M_creek_final is finally obtained and stored in system memory.

[0038] As described in step S4, centerline extraction and network construction Skeletonization and Short Branch Pruning: The processor performs centerline extraction on M_creek_final. In this embodiment, the binary mask is first iteratively skeletonized (skel) using the bwmorph function in the MATLAB environment until the skeleton width is a single pixel. Short branch pruning is then performed: the pruning intensity is set to 10 iterations, and skeleton connected branches with a length of less than 5 pixels are deleted after each iteration. After pruning, the final tidal channel mask is written to a non-volatile storage device in GeoTIFF format, with affine transformation parameters set according to the local Cartesian coordinates of the model's X and Y axes (raster origin set to (0,0), pixel size set to 5 m).

[0039] Network node and link generation: The processor calls the Python-based RivGraph network analysis library to read the GeoTIFF mask. RivGraph internally re-executes skeletonization and identifies three types of nodes based on the adjacency relationships of the skeleton cells: End points: Nodes with 1 adjacent cell; Bifurcations: Nodes with 3 or more adjacent pixels; Outlet: A node located at the mask boundary and with ≥1 adjacent pixel.

[0040] A continuous string of skeleton pixels between adjacent nodes is divided into a link. Each link records the following attributes in memory: link number. Starting node number End point node number Coordinate sequence sampled along the centerline {( Total link length (Accumulated by the Euclidean distance between adjacent sampling points, i.e.) (Where K is the total number of sampling points for the link). Simultaneously, RivGraph calculates the local width at each central sampling point by searching the left and right bank boundaries of the tidal channel along the normal direction of the link's centerline. The width value is bound to the corresponding sampling point. The network calculation results are output in GeoJSON format, containing complete geometric and attribute information for nodes and links.

[0041] It should be emphasized that this invention does not limit the use of the RivGraph library. In alternative embodiments, the same network construction function can be achieved using ArcGIS's ArcHydro toolset, GRASS GIS's r.thin and v.net modules, or an autonomous graph search algorithm based on morphological distance transformation.

[0042] As described in step S5, multi-scale index calculation The local cross-sectional morphology indicators are calculated as follows: The center sampling point and elevation binding processor reads the network GeoJSON file and extracts all center sampling points of each link, along with their corresponding left and right bank endpoint coordinates and local width. First, duplicate sampling records with completely identical center row and column coordinates are removed. Then, based on the X and Y coordinates of each sampling point, its row and column indices in the bed elevation raster are calculated (through the inverse operation of affine transformation). All covered bed elevation values ​​are sampled along the line connecting the left and right bank endpoints at pixel intervals, forming a set of discrete profile elevation sequence { , }, where m is the number of sampling points for this profile. Let be the bed surface elevation value at the j-th sampling point (unit: m).

[0043] If the elevation value of a certain endpoint in the cross-section parameter calculation is exactly equal to the artificial boundary elevation of the model (6m in this embodiment), then the endpoint is determined to be on the artificial boundary. The elevation of this endpoint is then replaced with the elevation value of the other valid endpoint to avoid boundary outliers affecting the profile statistics. For valid profiles, the following calculations are performed: (1) Local trench depth D (unit: m): The difference between the maximum and minimum elevations within the profile is taken as the calculation formula: in (.)and These are the maximum and minimum value operators, respectively.

[0044] (2) Cross-sectional area (Unit: m²): Based on the maximum elevation of the profile =max( Using the reference plane as a reference, the elevation difference at each sampling location is summed and then multiplied by the pixel size Δx (5 m in this embodiment). The calculation formula is as follows: Where Δx is the spatial resolution of the raster pixel (unit: m).

[0045] (3) Average width Average depth Average cross-sectional area : Take the arithmetic mean of the width, depth, and cross-sectional area of ​​all valid sampled profiles of the link.

[0046] (4) Width-to-depth ratio (Dimensionless): The calculation formula is: in This is the average width of the link (in meters). The average depth of the link (in meters).

[0047] If all elevation values ​​of a certain profile are invalid (such as -9999) or the sampling coordinates are out of bounds, the sample will be directly removed and not included in the statistics.

[0048] The aggregate statistics processor traverses all links, performing scene-level arithmetic mean calculations on the average width, depth, and cross-sectional area of ​​each link to obtain the overall average morphological index for that scene. Simultaneously, it counts the total number of pixels with a value of 1 in M_creek_final. Multiply by the single-pixel area Δx×Δy (25 m² in this example) to obtain the total area of ​​the tidal channel. (Unit: m²), the calculation formula is: The total length of the tidal channel is obtained by summing the lengths (len field) of all links. (Unit: m), the calculation formula is: Where K is the total number of links in this scenario. Let be the length of the k-th link.

[0049] The network and platform drainage indicators are calculated as follows: Drainage density calculation To effectively analyze area (600,000 m² in this example) is the denominator, representing the total length of the tidal channel. Calculate the displaced density for the molecule. (unit: ): in The total length of the tidal channel is given in meters. The total area of ​​the effective analysis domain (unit: m²).

[0050] The Unchannelized Path Length (UPL) calculation processor performs an Euclidean Distance Transform on M_creek_final. Specifically, for each non-ditch cell (value 0) in the raster, the Euclidean distance deucliddeuclid (unit: number of cells) to the nearest ditch cell (value 1) is calculated and multiplied by the cell scale Δx in the GeoTIFF affine transformation (5 m in this embodiment) to obtain the Unchannelized Path Length (UPL) in meters: UPL(x,y)= (x,y)×Δx Where UPL(x,y) is the non-channelized path length (in meters) of the non-tidal channel pixel at coordinates (x,y). (x, y) represents the Euclidean distance (in pixels) from the current pixel to the nearest tidal channel pixel, and Δx represents the pixel spatial resolution (in meters per pixel). The UPL value of the tidal channel pixel is set to 0. This calculation can be performed using GDAL's ComputeProximity function or SciPy's ndimage.distance_transform_edt, and both methods can be used interchangeably in this invention.

[0051] Platform drainage characteristic path length (OPL) calculation The processor extracts all pixel samples with positive UPL values ​​(UPL > 0), assuming a total sample size of n. These samples are then sorted in ascending order of UPL value. For the sample with sort index i (starting from 1), its empirical exceedance probability... The calculation is as follows: =1 in: i: Sample number sorted by UPL value from smallest to largest (i=1,2,…,n); n: UPL is the total number of positive (UPL>0) pixel samples; : The empirical exceedance probability of the i-th sample (dimensionless, with a value range of (0,1)), which means "the proportion of pixels whose UPL is greater than the current value".

[0052] Samples were selected within the exceedance probability window of 0.01 ≤ P ≤ 0.90, and a univariate linear regression was performed with ln(P) as the dependent variable and UPL as the independent variable. The model is as follows: ln(P) = a * UPL + b where: ln(P): The natural logarithm of the transcendental probability P (dimensionless). UPL: Non-channelized path length (in meters), used as the independent variable; a: Slope of the fitted line (unit: The physical meaning of α is the rate at which the exceedance probability decreases as the UPL increases; in a normal natural tidal channel system, a < 0. b: The intercept of the fitted line (dimensionless).

[0053] The slope a and intercept b are estimated using the least squares method, and the coefficient of determination is calculated. The platform drainage characteristic path length (OPL) is defined as the length of the fitted curve as it decays from the probability scale to the characteristic scale, and is calculated as follows: = Where a is the above fitting slope (unit: ... The unit of OPL is meters. The physical meaning of OPL is a characteristic scale of drainage path length. The smaller the value, the higher the drainage efficiency of the tidal channel network for the wetland platform.

[0054] When the number of valid positive UPL samples is less than 3, or the number of samples within the fitting window is less than 3, the processor determines that the OPL cannot be calculated in this scenario, records the result as NaN, and outputs a failure log. Preferably, the processor further checks whether the slope 'a' is less than 0 (to ensure physical meaning, i.e., the larger the UPL, the smaller the probability of exceeding the limit) and Whether the quality threshold (e.g., 0.8) is met; if not, it is marked as a low-quality result.

[0055] Finally, the processor calculates the geometric efficiency (GE) based on the drainage density. The geometric efficiency (GE) is a comprehensive index reflecting the water conveyance efficiency of the tidal channel network, calculated using the platform drainage characteristic path length (OPL) (GE is a dimensionless index). in Drainage density (unit: OPL is the characteristic path length of platform drainage (in meters). GE is a comprehensive measure of drainage efficiency; a larger value indicates higher drainage organization efficiency per unit drainage path length and unit tidal channel density. If Dd or OPL is 0 or an invalid value, GE is recorded as missing.

[0056] As described in step S6, results are summarized and quality control is performed. The main results table generator uses `standard_run_num` as the primary key and combines scenario parameters (tidal range, vegetation height, stem diameter, maximum density) and all calculated indicators (tidal channel area). Total length Average width Average depth Average cross-sectional area Aspect Ratio Drainage density The results are correlated with OPL and geometric efficiency (GE) and written to the master result table (in CSV or Excel format) in the storage device.

[0057] The data integrity verification and status marking processor performs the following five types of consistency checks for each scenario to ensure the reliability of batch processing results: File existence check: Check whether the input file pointed to by source_path in the scenario list actually exists in the storage device. If it is missing, record a "file missing" exception. Coordinate Reference System (CRS) Consistency Check: This check examines the affine transformation parameters and coordinate definitions in the GeoTIFF file output for each scenario, confirming that all scenarios use a unified local Cartesian coordinate system (X, Y) in meters. If an offset of the coordinate origin or inconsistent affine transformation parameters is detected, a warning message is recorded and subsequent processing of that scenario is stopped to prevent spatial misalignment from causing errors in index calculation.

[0058] Raster size consistency check: Read the GeoTransform parameter in the GeoTIFF file header of each scene, extract the cell width Δx and cell height Δy, compare them with the cell_size recorded in the scene list, and if the relative deviation exceeds 1%, trigger a warning and record the actual value; Unit consistency verification: All numerical indicators are checked for unit consistency, ensuring length is uniformly expressed as meters (m), area as square meters (m²), and density as density per meter (m²). If the original data unit is detected to be kilometers (km) or hectares (ha), the unit conversion is automatically performed (1 km = 1000 m, 1 ha = 10000 m²) and the conversion log is recorded; Number mapping integrity check: Check whether the standard_run_num in the scenario list corresponds one-to-one with the number in the output file name, check for records with skipped, duplicate or unmapped numbers, and record anomalies if mismatched numbers are found.

[0059] Based on the above verification results, the processor sets multiple existence flag fields for each scenario: has_area: Whether the trench mask file was successfully generated and its area. >0; has_links: Whether the network link data is not empty (total number of links K>0); has_opl: Whether the OPL fit is successful and of good quality (R²≥0.8 and a<0); has_results: This indicates that all three of the above items are valid.

[0060] Based on these markers, the overall data status is categorized into one of three types: "Complete" (all markers are True), "Partially Missing" (some markers are False), or "Missing" (critical files do not exist or all markers are False).

[0061] The exception log output processor captures all file read exceptions, coordinate out-of-bounds exceptions, and fitting calculation exceptions during processing, and records the exception type, occurrence time, associated scenario number, and detailed error information in the error.log text file. All verification results, unit conversion records, and warning messages are written to the qc_report.log file, which, together with status flags, constitute a complete quality control system.

[0062] Results Application and Scenario Comparison After completing all the above steps, the processor's output can be directly used to evaluate the restoration effect of coastal wetlands. For example, the 32 scenarios can be grouped according to tidal range (1 m, 2 m, 3 m, 5 m) and vegetation configuration, and plotted with... Figure 4 The tidal channel network morphology diagram and box plots for each indicator shown can intuitively identify the quantitative impact of increased tidal range or vegetation density on average channel width, width-to-depth ratio, drainage density, and geometric efficiency. Combined with the final vegetation cover, where P is the final plant density, vegetation cover is calculated based on P > 1 stem. The calculation of grid ratios can further differentiate between three ecological response types: "continuous vegetation colonization," "vegetation decline," and "no vegetation baseline," providing a quantitative basis for selecting tidal conditions and vegetation configuration schemes in restoration projects such as managed realignment.

[0063] The above embodiments are only used to illustrate the technical methods of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical methods of the present invention without departing from the spirit and scope of the technical methods of the present invention.

Claims

1. A method for automatic extraction and morphological index calculation of tidal channel networks, characterized in that, Includes the following steps: Establish a scenario list, read the input data for each scenario according to the scenario list, and align and crop all input data based on a unified spatial reference and effective analysis domain. The input data includes: bed elevation field, water depth field and / or flow velocity field. Within the effective analysis domain, the physical fields in the input data are dimensionless, and at least two dimensionless physical fields are fused to generate a fused feature field; the physical fields include: the bed elevation field, the water depth field, and / or the flow velocity field. A segmentation threshold is determined based on the statistical distribution of the fused feature field within the effective analysis domain, and the fused feature field is segmented according to the segmentation threshold to obtain an initial tidal channel mask; the initial tidal channel mask is then subjected to connectivity filtering and / or noise removal to obtain an optimized tidal channel mask. The optimized tidal channel mask is subjected to centerline extraction, and nodes and links are identified based on the extracted centerline to generate a tidal channel network data structure containing the attribute information of the nodes and links; the link attribute information includes link length, link width, and spatial coordinates. Based on the spatial correspondence between the tidal channel mask, the network data structure, and the bed surface elevation field, local morphological indices and network scale indices are calculated.

2. The method for automatic extraction and morphological index calculation of tidal channel networks according to claim 1, characterized in that, The step of generating the fused feature field includes: The bed elevation field is reverse normalized to obtain higher normalized characteristic values ​​in low-lying areas; the water depth field and / or the flow velocity field are forward normalized to obtain higher normalized characteristic values ​​in areas with large water depth and / or high flow velocity. The normalized bed elevation field, water depth field, and / or flow velocity field are summed pixel-by-pixel by weighted summation or by direct summation to obtain the fused feature field.

3. The method for automatic extraction and morphological index calculation of tidal channel networks according to claim 1, characterized in that, The segmentation threshold is dynamically determined based on the standard deviation of the fusion feature field within the effective analysis domain: When the standard deviation is greater than the preset discriminant value, a lower quantile is used as the segmentation threshold; when the standard deviation is not greater than the preset discriminant value, a higher quantile is used as the segmentation threshold.

4. The method for automatic extraction and morphological index calculation of tidal channel networks according to claim 1, characterized in that, The noise removal includes deleting connected objects in the initial tidal channel mask whose number of pixels in the connected region is less than a preset pixel number threshold.

5. The method for automatic extraction and morphological index calculation of tidal channel networks according to claim 1, characterized in that, The centerline extraction includes: skeletalizing the optimized tidal channel mask and pruning short branches on the skeletalized result; the nodes include: endpoints, branching points, and exit points; the link is a continuous centerline segment between adjacent nodes.

6. The method for automatic extraction and morphological index calculation of tidal channel networks according to claim 1, characterized in that, The local morphological indicators include at least one of the following: tidal channel width, tidal channel depth, cross-sectional area, and width-to-depth ratio; when calculating the local morphological indicators, the profile elevation values ​​at the corresponding locations are extracted from the bed elevation field based on the link centerline sampling points and the left and right bank endpoints in the network data structure.

7. The method for automatic extraction and morphological index calculation of tidal channel networks according to claim 1, characterized in that, The network metric includes at least one of drainage density, non-channelized path length, platform drainage characteristic path length, and geometric efficiency. Wherein, the drainage density is the ratio of the total length of all links to the area of ​​the effective analysis domain; The non-channelized path length is the distance from the non-tidal groove pixel to the nearest tidal groove pixel in the tidal groove mask; The characteristic path length of the platform drainage is obtained by fitting the exceedance probability distribution of the non-channelized path length; The geometric efficiency is the reciprocal of the product of the drainage density and the characteristic path length of the platform drainage.

8. The method for automatic extraction and morphological index calculation of tidal channel networks according to claim 1, characterized in that, It also includes the following steps: The scenario identifier is associated with and stored in conjunction with the calculated local morphological index and network scale index, and a state tag for characterizing data integrity is generated.

9. The method for automatic extraction and morphological index calculation of tidal channel networks according to claim 8, characterized in that, The step of associating and storing the scenario identifier with the calculated local morphological index and network scale index, and generating a state label for characterizing data integrity, further includes: Check whether the files in the scenario list exist, whether the coordinate reference system is consistent, whether the grid size is uniform, whether the units are correct, or whether the numbering mapping is complete; and generate a processing log based on the check results and the status markers.

10. The method for automatic extraction and morphological index calculation of tidal channel networks according to claim 1, characterized in that, The input data comes from the output of the hydrodynamic-sediment-vegetation coupled numerical model at the target time, and the numerical model is the Delft3D model.

Images